text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# print.jl
# All "pretty printers" for JuMP types.
# - Delegates to appropriate handler methods for REPL or IJulia.
# - These handler methods then pass the correct symbols to use into a
#   generic string builder. The IJulia handlers will also wrap in MathJax
#   start/close tags.
# - To find printing code for a type in this file, search for `## TypeName`
# - Code here does not need to be fast, in fact simplicity trumps speed
#   within reason as this code is thorny enough as it is.
# - Corresponding tests are in test/print.jl, although test/operator.jl
#   is also testing the constraint/expression code extensively as well.
# - Base.print and Base.string both delegate to Base.show, if they are not
#   separately defined.
#############################################################################

using Printf

# Used for dispatching
"""
    PrintMode

An abstract type used for dispatching printing.
"""
abstract type PrintMode end

"""
    REPLMode

A type used for dispatching printing. Produces plain-text representations.
"""
abstract type REPLMode <: PrintMode end

"""
    IJuliaMode

A type used for dispatching printing. Produces LaTeX representations.
"""
abstract type IJuliaMode <: PrintMode end

Base.show(io::IO, model::AbstractModel) = _print_summary(io, model)

struct _LatexModel{T<:AbstractModel}
    model::T
end

"""
    latex_formulation(model::AbstractModel)

Wrap `model` in a type so that it can be pretty-printed as `text/latex` in a
notebook like IJulia, or in Documenter.

To render the model, end the cell with `latex_formulation(model)`, or call
`display(latex_formulation(model))` in to force the display of the model from
inside a function.
"""
latex_formulation(model::AbstractModel) = _LatexModel(model)

Base.show(io::IO, model::_LatexModel) = _print_latex(io, model.model)
Base.show(io::IO, ::MIME"text/latex", model::_LatexModel) = show(io, model)

function Base.print(model::AbstractModel)
    for d in Base.Multimedia.displays
        if Base.Multimedia.displayable(d, "text/latex") &&
           startswith("$(typeof(d))", "IJulia.")
            return display(d, "text/latex", latex_formulation(model))
        end
    end
    return _print_model(stdout, model)
end
Base.print(io::IO, model::AbstractModel) = _print_model(io, model)

# Whether something is zero or not for the purposes of printing it
# oneunit is useful e.g. if coef is a Unitful quantity.
_is_zero_for_printing(coef) = abs(coef) < 1e-10 * oneunit(coef)
# Whether something is one or not for the purposes of printing it.
_is_one_for_printing(coef) = _is_zero_for_printing(abs(coef) - oneunit(coef))
_sign_string(coef) = coef < zero(coef) ? " - " : " + "

# Helper function that rounds carefully for the purposes of printing
# e.g.   5.3  =>  5.3
#        1.0  =>  1
function _string_round(f::Float64)
    iszero(f) && return "0" # strip sign off zero
    str = string(f)
    return length(str) >= 2 && str[end-1:end] == ".0" ? str[1:end-2] : str
end
_string_round(f) = string(f)

# REPL-specific symbols
# Anything here: https://en.wikipedia.org/wiki/Windows-1252
# should probably work fine on Windows
function _math_symbol(::Type{REPLMode}, name::Symbol)
    if name == :leq
        return Sys.iswindows() ? "<=" : "≤"
    elseif name == :geq
        return Sys.iswindows() ? ">=" : "≥"
    elseif name == :eq
        return Sys.iswindows() ? "==" : "="
    elseif name == :times
        return "*"
    elseif name == :sq
        return "²"
    elseif name == :ind_open
        return "["
    elseif name == :ind_close
        return "]"
    elseif name == :for_all
        return Sys.iswindows() ? "for all" : "∀"
    elseif name == :in
        return Sys.iswindows() ? "in" : "∈"
    elseif name == :open_set
        return "{"
    elseif name == :dots
        return Sys.iswindows() ? ".." : "…"
    elseif name == :close_set
        return "}"
    elseif name == :union
        return Sys.iswindows() ? "or" : "∪"
    elseif name == :infty
        return Sys.iswindows() ? "Inf" : "∞"
    elseif name == :open_rng
        return "["
    elseif name == :close_rng
        return "]"
    elseif name == :integer
        return "integer"
    elseif name == :succeq0
        return " is semidefinite"
    elseif name == :Vert
        return Sys.iswindows() ? "||" : "‖"
    elseif name == :sub2
        return Sys.iswindows() ? "_2" : "₂"
    else
        error("Internal error: Unrecognized symbol $name.")
    end
end

# IJulia-specific symbols.
function _math_symbol(::Type{IJuliaMode}, name::Symbol)
    if name == :leq
        return "\\leq"
    elseif name == :geq
        return "\\geq"
    elseif name == :eq
        return "="
    elseif name == :times
        return "\\times "
    elseif name == :sq
        return "^2"
    elseif name == :ind_open
        return "_{"
    elseif name == :ind_close
        return "}"
    elseif name == :for_all
        return "\\quad\\forall"
    elseif name == :in
        return "\\in"
    elseif name == :open_set
        return "\\{"
    elseif name == :dots
        return "\\dots"
    elseif name == :close_set
        return "\\}"
    elseif name == :union
        return "\\cup"
    elseif name == :infty
        return "\\infty"
    elseif name == :open_rng
        return "\\["
    elseif name == :close_rng
        return "\\]"
    elseif name == :integer
        return "\\in \\mathbb{Z}"
    elseif name == :succeq0
        return "\\succeq 0"
    elseif name == :Vert
        return "\\Vert"
    elseif name == :sub2
        return "_2"
    else
        error("Internal error: Unrecognized symbol $name.")
    end
end

_wrap_in_math_mode(str) = "\$\$ $str \$\$"
_wrap_in_inline_math_mode(str) = "\$ $str \$"

_plural(n) = (isone(n) ? "" : "s")

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

"""
    _print_summary(io::IO, model::AbstractModel)

Print a plain-text summary of `model` to `io`.

An `AbstractModel` subtype should implement `show_objective_function_summary`,
`show_constraints_summary` and `show_backend_summary` for this method to work.
"""
function _print_summary(io::IO, model::AbstractModel)
    println(io, "A JuMP Model")
    sense = objective_sense(model)
    if sense == MOI.MAX_SENSE
        print(io, "Maximization")
    elseif sense == MOI.MIN_SENSE
        print(io, "Minimization")
    else
        print(io, "Feasibility")
    end
    println(io, " problem with:")
    # TODO: Use MOI.Name for the name of a JuMP model.
    println(
        io,
        "Variable",
        _plural(num_variables(model)),
        ": ",
        num_variables(model),
    )
    if sense != MOI.FEASIBILITY_SENSE
        show_objective_function_summary(io, model)
    end
    show_constraints_summary(io, model)
    show_backend_summary(io, model)
    names_in_scope = sort(collect(keys(object_dictionary(model))))
    if !isempty(names_in_scope)
        println(io)
        print(
            io,
            "Names registered in the model: ",
            join(string.(names_in_scope), ", "),
        )
    end
end

"""
    show_backend_summary(io::IO, model::Model)

Print a summary of the optimizer backing `model`.

`AbstractModel`s should implement this method.
"""
function show_backend_summary(io::IO, model::Model)
    model_mode = mode(model)
    println(io, "Model mode: ", model_mode)
    if model_mode == MANUAL || model_mode == AUTOMATIC
        println(io, "CachingOptimizer state: ", MOIU.state(backend(model)))
    end
    # The last print shouldn't have a new line
    print(io, "Solver name: ", solver_name(model))
    return
end

"""
    _print_model(io::IO, model::AbstractModel)

Print a plain-text formulation of `model` to `io`.

An `AbstractModel` subtype should implement `objective_function_string`,
`constraints_string` and `_nl_subexpression_string` for this method to work.
"""
function _print_model(io::IO, model::AbstractModel)
    sense = objective_sense(model)
    if sense == MOI.MAX_SENSE
        print(io, "Max ")
        println(io, objective_function_string(REPLMode, model))
    elseif sense == MOI.MIN_SENSE
        print(io, "Min ")
        println(io, objective_function_string(REPLMode, model))
    else
        println(io, "Feasibility")
    end
    println(io, "Subject to")
    for constraint in constraints_string(REPLMode, model)
        println(io, " ", replace(constraint, '\n' => "\n "))
    end
    # TODO: Generalize this when similar functionality is needed for
    # `AbstractModel`.
    nl_subexpressions = _nl_subexpression_string(REPLMode, model)
    if !isempty(nl_subexpressions)
        println(io, "With NL expressions")
    end
    for expr in nl_subexpressions
        println(io, " ", expr)
    end
    return
end

"""
    _print_latex(io::IO, model::AbstractModel)

Print a LaTeX formulation of `model` to `io`.

An `AbstractModel` subtype should implement `objective_function_string`,
`constraints_string` and `_nl_subexpression_string` for this method to work.
"""
function _print_latex(io::IO, model::AbstractModel)
    println(io, "\$\$ \\begin{aligned}")
    sense = objective_sense(model)
    if sense == MOI.MAX_SENSE
        print(io, "\\max\\quad & ")
        println(io, objective_function_string(IJuliaMode, model), "\\\\")
    elseif sense == MOI.MIN_SENSE
        print(io, "\\min\\quad & ")
        println(io, objective_function_string(IJuliaMode, model), "\\\\")
    else
        println(io, "\\text{feasibility}\\\\")
    end
    constraints = constraints_string(IJuliaMode, model)
    if !isempty(constraints)
        print(io, "\\text{Subject to} \\quad")
    end
    for constraint in constraints
        println(io, " & ", constraint, "\\\\")
    end
    # TODO: Generalize this when similar functionality is needed for
    # `AbstractModel`.
    nl_subexpressions = _nl_subexpression_string(IJuliaMode, model)
    if !isempty(nl_subexpressions)
        print(io, "\\text{With NL expressions} \\quad")
    end
    for expr in nl_subexpressions
        println(io, " & ", expr, "\\\\")
    end
    return print(io, "\\end{aligned} \$\$")
end

# TODO(odow): deprecate this?
function model_string(print_mode, model::AbstractModel)
    if print_mode == IJuliaMode
        return sprint(_print_latex, model)
    else
        return sprint(_print_model, model)
    end
end

_nl_subexpression_string(print_mode, ::AbstractModel) = String[]

function _nl_subexpression_string(print_mode, model::Model)
    strings = String[]
    if model.nlp_data !== nothing
        num_subexpressions = length(model.nlp_data.nlexpr)::Int
        for k in 1:num_subexpressions
            ex = model.nlp_data.nlexpr[k]
            expr_string = _tape_to_expr(
                model,
                1, # start index in the expression
                ex.nd,
                adjmat(ex.nd),
                ex.const_values,
                [], # parameter_values (not used)
                [], # subexpressions (not needed because !splat_subexpressions)
                model.nlp_data.user_operators,
                false, # generic_variable_names
                false, # splat_subexpressions
                print_mode,
            )
            if print_mode == IJuliaMode
                expr_name = "subexpression_{$k}"
            else
                expr_name = "subexpression[$k]"
            end
            push!(strings, "$expr_name: $expr_string")
        end
    end
    return strings
end

"""
    show_objective_function_summary(io::IO, model::AbstractModel)

Write to `io` a summary of the objective function type.
"""
function show_objective_function_summary(io::IO, model::Model)
    nlobj = _nlp_objective_function(model)
    print(io, "Objective function type: ")
    if nlobj === nothing
        println(io, objective_function_type(model))
    else
        println(io, "Nonlinear")
    end
end

"""
    objective_function_string(print_mode, model::AbstractModel)::String

Return a `String` describing the objective function of the model.
"""
function objective_function_string(print_mode, model::Model)
    nlobj = _nlp_objective_function(model)
    if nlobj === nothing
        return function_string(print_mode, objective_function(model))
    else
        return nl_expr_string(model, print_mode, nlobj)
    end
end

struct _SolutionSummary
    verbose::Bool
    solver::String
    # Status
    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    dual_status::MOI.ResultStatusCode
    raw_status::String
    result_count::Int
    has_values::Bool
    has_duals::Bool
    # Candidate solution
    objective_value::Float64
    objective_bound::Union{Missing,Float64}
    dual_objective_value::Union{Missing,Float64}
    primal_solution::Union{Missing,Dict{String,Float64}}
    dual_solution::Union{Missing,Dict{String,Float64}}
    # Work counters
    solve_time::Float64
    barrier_iterations::Union{Missing,Int}
    simplex_iterations::Union{Missing,Int}
    node_count::Union{Missing,Int}
end

"""
    solution_summary(model::Model; verbose::Bool = false)

Return a struct that can be used print a summary of the solution.

If `verbose=true`, write out the primal solution for every variable and the
dual solution for every constraint, excluding those with empty names.

## Examples

When called at the REPL, the summary is automatically printed:
```julia
julia> solution_summary(model)
[...]
```

Use `print` to force the printing of the summary from inside a function:
```julia
function foo(model)
    print(solution_summary(model))
    return
end
```
"""
function solution_summary(model::Model; verbose::Bool = false)
    return _SolutionSummary(
        verbose,
        solver_name(model),
        termination_status(model),
        primal_status(model),
        dual_status(model),
        raw_status(model),
        result_count(model),
        has_values(model),
        has_duals(model),
        objective_value(model),
        _try_get(objective_bound, model),
        _try_get(dual_objective_value, model),
        verbose ? _get_solution_dict(model) : missing,
        verbose ? _get_constraint_dict(model) : missing,
        solve_time(model),
        _try_get(simplex_iterations, model),
        _try_get(barrier_iterations, model),
        _try_get(node_count, model),
    )
end

"""
    Base.show([io::IO], summary::SolutionSummary; verbose::Bool = false)

Write a summary of the solution results to `io` (or to `stdout` if `io` is not
given).
"""
function Base.show(io::IO, summary::_SolutionSummary)
    println(io, "* Solver : ", summary.solver)
    println(io)
    _show_status_summary(io, summary)
    _show_candidate_solution_summary(io, summary)
    _show_work_counters_summary(io, summary)
    return
end

function _show_status_summary(io::IO, summary::_SolutionSummary)
    println(io, "* Status")
    println(io, "  Termination status : ", summary.termination_status)
    println(io, "  Primal status      : ", summary.primal_status)
    println(io, "  Dual status        : ", summary.dual_status)
    if summary.verbose
        println(io, "  Result count       : ", summary.result_count)
        println(io, "  Has duals          : ", summary.has_duals)
    end
    println(io, "  Message from the solver:")
    println(io, "  \"", summary.raw_status, "\"")
    println(io)
    return
end

function _show_candidate_solution_summary(io::IO, summary::_SolutionSummary)
    println(io, "* Candidate solution")
    println(io, "  Objective value      : ", summary.objective_value)
    _print_if_not_missing(
        io,
        "  Objective bound      : ",
        summary.objective_bound,
    )
    _print_if_not_missing(
        io,
        "  Dual objective value : ",
        summary.dual_objective_value,
    )
    if summary.verbose && summary.has_values
        println(io, "  Primal solution :")
        for variable_name in sort(collect(keys(summary.primal_solution)))
            println(
                io,
                "    ",
                variable_name,
                " : ",
                summary.primal_solution[variable_name],
            )
        end
    end
    if summary.verbose && summary.has_duals
        println(io, "  Dual solution :")
        for constraint_name in sort(collect(keys(summary.dual_solution)))
            println(
                io,
                "    ",
                constraint_name,
                " : ",
                summary.dual_solution[constraint_name],
            )
        end
    end
    println(io)
    return
end

function _show_work_counters_summary(io::IO, summary::_SolutionSummary)
    println(io, "* Work counters")
    println(io, "  Solve time (sec)   : ", @sprintf("%.5f", summary.solve_time))
    _print_if_not_missing(
        io,
        "  Simplex iterations : ",
        summary.simplex_iterations,
    )
    _print_if_not_missing(
        io,
        "  Barrier iterations : ",
        summary.barrier_iterations,
    )
    _print_if_not_missing(io, "  Node count         : ", summary.node_count)
    return
end

function _get_solution_dict(model)
    dict = Dict{String,Float64}()
    if has_values(model)
        for x in all_variables(model)
            variable_name = name(x)
            if !isempty(variable_name)
                dict[variable_name] = value(x)
            end
        end
    end
    return dict
end

function _get_constraint_dict(model)
    dict = Dict{String,Float64}()
    if has_duals(model)
        for (F, S) in list_of_constraint_types(model)
            for constraint in all_constraints(model, F, S)
                constraint_name = name(constraint)
                if !isempty(constraint_name)
                    dict[constraint_name] = dual(constraint)
                end
            end
        end
    end
    return dict
end

function _try_get(f, model)
    try
        return f(model)
    catch
        return missing
    end
end

_print_if_not_missing(io, header, ::Missing) = nothing
_print_if_not_missing(io, header, value) = println(io, header, value)

#------------------------------------------------------------------------
## VariableRef
#------------------------------------------------------------------------

function function_string(::Type{REPLMode}, v::AbstractVariableRef)
    var_name = name(v)
    if !isempty(var_name)
        return var_name
    else
        return "noname"
    end
end
function function_string(::Type{IJuliaMode}, v::AbstractVariableRef)
    var_name = name(v)
    if isempty(var_name)
        return "noname"
    end
    # We need to escape latex math characters that appear in the name.
    # However, it's probably impractical to catch everything, so let's just
    # escape the common ones:
    # Escape underscores to prevent them being treated as subscript markers.
    var_name = replace(var_name, "_" => "\\_")
    # Escape carets to prevent them being treated as superscript markers.
    var_name = replace(var_name, "^" => "\\^")
    # Convert any x[args] to x_{args} so that indices on x print as subscripts.
    m = match(r"^(.*)\[(.+)\]$", var_name)
    if m !== nothing
        var_name = m[1] * "_{" * m[2] * "}"
    end

    return var_name
end

#------------------------------------------------------------------------
## GenericAffExpr
#------------------------------------------------------------------------

function function_string(mode, a::GenericAffExpr, show_constant = true)
    # If the expression is empty, return the constant (or 0)
    if length(linear_terms(a)) == 0
        return show_constant ? _string_round(a.constant) : "0"
    end

    term_str = Array{String}(undef, 2 * length(linear_terms(a)))
    elm = 1

    for (coef, var) in linear_terms(a)
        pre = _is_one_for_printing(coef) ? "" : _string_round(abs(coef)) * " "

        term_str[2*elm-1] = _sign_string(coef)
        term_str[2*elm] = string(pre, function_string(mode, var))
        elm += 1
    end

    if elm == 1
        # Will happen with cancellation of all terms
        # We should just return the constant, if its desired
        return show_constant ? _string_round(a.constant) : "0"
    else
        # Correction for very first term - don't want a " + "/" - "
        term_str[1] = (term_str[1] == " - ") ? "-" : ""
        ret = join(term_str[1:2*(elm-1)])
        if !_is_zero_for_printing(a.constant) && show_constant
            ret = string(
                ret,
                _sign_string(a.constant),
                _string_round(abs(a.constant)),
            )
        end
        return ret
    end
end

#------------------------------------------------------------------------
## GenericQuadExpr
#------------------------------------------------------------------------

function function_string(mode, q::GenericQuadExpr)
    length(quad_terms(q)) == 0 && return function_string(mode, q.aff)

    # Odd terms are +/i, even terms are the variables/coeffs
    term_str = Array{String}(undef, 2 * length(quad_terms(q)))
    elm = 1
    if length(term_str) > 0
        for (coef, var1, var2) in quad_terms(q)
            pre =
                _is_one_for_printing(coef) ? "" : _string_round(abs(coef)) * " "

            x = function_string(mode, var1)
            y = function_string(mode, var2)

            term_str[2*elm-1] = _sign_string(coef)
            term_str[2*elm] = "$pre$x"
            if x == y
                term_str[2*elm] *= _math_symbol(mode, :sq)
            else
                term_str[2*elm] *= string(_math_symbol(mode, :times), y)
            end
            if elm == 1
                # Correction for first term as there is no space
                # between - and variable coefficient/name
                term_str[1] = coef < zero(coef) ? "-" : ""
            end
            elm += 1
        end
    end
    ret = join(term_str[1:2*(elm-1)])

    aff_str = function_string(mode, q.aff)
    if aff_str == "0"
        return ret
    else
        if aff_str[1] == '-'
            return string(ret, " - ", aff_str[2:end])
        else
            return string(ret, " + ", aff_str)
        end
    end
end

#------------------------------------------------------------------------
## Constraints
#------------------------------------------------------------------------

"""
    show_constraints_summary(io::IO, model::AbstractModel)

Write to `io` a summary of the number of constraints.
"""
function show_constraints_summary(io::IO, model::Model)
    for (F, S) in list_of_constraint_types(model)
        n_constraints = num_constraints(model, F, S)
        println(
            io,
            "`$F`-in-`$S`: $n_constraints constraint",
            _plural(n_constraints),
        )
    end
    if !iszero(num_nl_constraints(model))
        println(
            io,
            "Nonlinear: ",
            num_nl_constraints(model),
            " constraint",
            _plural(num_nl_constraints(model)),
        )
    end
end

"""
    constraints_string(print_mode, model::AbstractModel)::Vector{String}

Return a list of `String`s describing each constraint of the model.
"""
function constraints_string(print_mode, model::Model)
    strings = String[]
    for (F, S) in list_of_constraint_types(model)
        for cref in all_constraints(model, F, S)
            push!(
                strings,
                constraint_string(print_mode, cref, in_math_mode = true),
            )
        end
    end
    if model.nlp_data !== nothing
        for nl_constraint in model.nlp_data.nlconstr
            push!(
                strings,
                nl_constraint_string(model, print_mode, nl_constraint),
            )
        end
    end
    return strings
end

## Notes for extensions
# For a `ConstraintRef{ModelType, IndexType}` where `ModelType` is not
# `JuMP.Model` or `IndexType` is not `MathOptInterface.ConstraintIndex`, the
# methods `JuMP.name` and `JuMP.constraint_object` should be implemented for
# printing to work. If the `AbstractConstraint` returned by `constraint_object`
# is not `JuMP.ScalarConstraint` nor `JuMP.VectorConstraint`, then either
# `JuMP.jump_function` or `JuMP.function_string` and either `JuMP.moi_set` or
# `JuMP.in_set_string` should be implemented.
function Base.show(io::IO, ref::ConstraintRef)
    return print(io, constraint_string(REPLMode, ref))
end
function Base.show(io::IO, ::MIME"text/latex", ref::ConstraintRef)
    return print(io, constraint_string(IJuliaMode, ref))
end

"""
    function_string(print_mode::Type{<:JuMP.PrintMode},
                    func::Union{JuMP.AbstractJuMPScalar,
                                Vector{<:JuMP.AbstractJuMPScalar}})

Return a `String` representing the function `func` using print mode
`print_mode`.
"""
function function_string end

function Base.show(io::IO, f::AbstractJuMPScalar)
    return print(io, function_string(REPLMode, f))
end
function Base.show(io::IO, ::MIME"text/latex", f::AbstractJuMPScalar)
    return print(io, _wrap_in_math_mode(function_string(IJuliaMode, f)))
end

function function_string(print_mode, vector::Vector{<:AbstractJuMPScalar})
    return "[" * join(function_string.(print_mode, vector), ", ") * "]"
end

function function_string(
    ::Type{REPLMode},
    A::AbstractMatrix{<:AbstractJuMPScalar},
)
    str = sprint(show, MIME"text/plain"(), A)
    lines = split(str, '\n')
    # We drop the first line with the signature "m×n Array{...}:"
    lines = lines[2:end]
    # We replace the first space by an opening `[`
    lines[1] = '[' * lines[1][2:end]
    for i in 1:length(lines)
        lines[i] = lines[i] * (i == length(lines) ? ']' : ';')
    end
    return join(lines, '\n')
end

function function_string(
    print_mode::Type{IJuliaMode},
    A::AbstractMatrix{<:AbstractJuMPScalar},
)
    str = "\\begin{bmatrix}\n"
    for i in 1:size(A, 1)
        line = ""
        for j in 1:size(A, 2)
            if j != 1
                line *= " & "
            end
            if A isa Symmetric && i > j
                line *= "\\cdot"
            else
                line *= function_string(print_mode, A[i, j])
            end
        end
        str *= line * "\\\\\n"
    end
    return str * "\\end{bmatrix}"
end

"""
    function_string(print_mode::{<:JuMP.PrintMode},
                    constraint::JuMP.AbstractConstraint)

Return a `String` representing the function of the constraint `constraint`
using print mode `print_mode`.
"""
function function_string(print_mode, constraint::AbstractConstraint)
    f = reshape_vector(jump_function(constraint), shape(constraint))
    return function_string(print_mode, f)
end

"""
    in_set_string(print_mode::Type{<:PrintMode}, set)

Return a `String` representing the membership to the set `set` using print mode
`print_mode`.
"""
function in_set_string end

function in_set_string(print_mode, set::MOI.LessThan)
    return string(_math_symbol(print_mode, :leq), " ", set.upper)
end

function in_set_string(print_mode, set::MOI.GreaterThan)
    return string(_math_symbol(print_mode, :geq), " ", set.lower)
end

function in_set_string(print_mode, set::MOI.EqualTo)
    return string(_math_symbol(print_mode, :eq), " ", set.value)
end

function in_set_string(print_mode, set::MOI.Interval)
    return string(
        _math_symbol(print_mode, :in),
        " ",
        _math_symbol(print_mode, :open_rng),
        set.lower,
        ", ",
        set.upper,
        _math_symbol(print_mode, :close_rng),
    )
end

in_set_string(print_mode, ::MOI.ZeroOne) = "binary"
in_set_string(print_mode, ::MOI.Integer) = "integer"

in_set_string(::Type{IJuliaMode}, ::MOI.ZeroOne) = "\\in \\{0, 1\\}"
in_set_string(::Type{IJuliaMode}, ::MOI.Integer) = "\\in \\mathbb{Z}"

function in_set_string(print_mode, set::Union{PSDCone,MOI.AbstractSet})
    # Use an `if` here instead of multiple dispatch to avoid ambiguity errors.
    if print_mode == REPLMode
        return _math_symbol(print_mode, :in) * " $(set)"
    else
        set_str = replace(replace(string(set), "{" => "\\{"), "}" => "\\}")
        return "\\in \\text{$(set_str)}"
    end
end

"""
    in_set_string(print_mode::Type{<:PrintMode}, constraint::AbstractConstraint)

Return a `String` representing the membership to the set of the constraint
`constraint` using print mode `print_mode`.
"""
function in_set_string(print_mode, constraint::AbstractConstraint)
    set = reshape_set(moi_set(constraint), shape(constraint))
    return in_set_string(print_mode, set)
end

function constraint_string(print_mode, constraint_object::AbstractConstraint)
    func_str = function_string(print_mode, constraint_object)
    in_set_str = in_set_string(print_mode, constraint_object)
    if print_mode == REPLMode
        lines = split(func_str, '\n')
        lines[1+div(length(lines), 2)] *= " " * in_set_str
        return join(lines, '\n')
    else
        return func_str * " " * in_set_str
    end
end

function constraint_string(
    print_mode,
    constraint_name::String,
    constraint_object::AbstractConstraint;
    in_math_mode::Bool = false,
)
    prefix = isempty(constraint_name) ? "" : constraint_name * " : "
    constraint_str = constraint_string(print_mode, constraint_object)
    if print_mode == IJuliaMode
        if in_math_mode
            return constraint_str
        elseif isempty(prefix)
            return _wrap_in_math_mode(constraint_str)
        else
            return prefix * _wrap_in_inline_math_mode(constraint_str)
        end
    else
        return prefix * constraint_str
    end
end

"""
    constraint_string(
        print_mode,
        ref::ConstraintRef;
        in_math_mode::Bool = false)

Return a string representation of the constraint `ref`, given the `print_mode`.

`print_mode` should be `IJuliaMode` or `REPLMode`.
"""
function constraint_string(print_mode, ref::ConstraintRef; in_math_mode = false)
    return constraint_string(
        print_mode,
        name(ref),
        constraint_object(ref),
        in_math_mode = in_math_mode,
    )
end

#------------------------------------------------------------------------
## _NonlinearExprData
#------------------------------------------------------------------------

"""
    nl_expr_string(model::Model, print_mode, c::_NonlinearExprData)

Return a string representation of the nonlinear expression `c` belonging to
`model`, given the `print_mode`. `print_mode` should be `IJuliaMode` or
`REPLMode`.
"""
function nl_expr_string(model::Model, print_mode, c::_NonlinearExprData)
    ex = _tape_to_expr(
        model,
        1,
        c.nd,
        adjmat(c.nd),
        c.const_values,
        [],
        [],
        model.nlp_data.user_operators,
        false,
        false,
        print_mode,
    )
    if print_mode == IJuliaMode
        ex = _latexify_exponentials(ex)
    end
    return string(ex)
end

# Change x ^ -2.0 to x ^ {-2.0}
# x ^ (x ^ 2.0) to x ^ {x ^ {2.0}}
# and so on
_latexify_exponentials(ex) = ex
function _latexify_exponentials(ex::Expr)
    for i in 1:length(ex.args)
        ex.args[i] = _latexify_exponentials(ex.args[i])
    end
    if length(ex.args) == 3 && ex.args[1] == :^
        ex.args[3] = Expr(:braces, ex.args[3])
    end
    return ex
end
#------------------------------------------------------------------------
## _NonlinearConstraint
#------------------------------------------------------------------------
const NonlinearConstraintRef = ConstraintRef{Model,NonlinearConstraintIndex}

function Base.show(io::IO, c::NonlinearConstraintRef)
    return print(
        io,
        nl_constraint_string(
            c.model,
            REPLMode,
            c.model.nlp_data.nlconstr[c.index.value],
        ),
    )
end

function Base.show(io::IO, ::MIME"text/latex", c::NonlinearConstraintRef)
    constraint = c.model.nlp_data.nlconstr[c.index.value]
    return print(
        io,
        _wrap_in_math_mode(
            nl_constraint_string(c.model, IJuliaMode, constraint),
        ),
    )
end

# TODO: Printing is inconsistent between regular constraints and nonlinear
# constraints because nonlinear constraints don't have names.
"""
    nl_constraint_string(model::Model, mode, c::_NonlinearConstraint)

Return a string representation of the nonlinear constraint `c` belonging to
`model`, given the `print_mode`. `print_mode` should be `IJuliaMode` or
`REPLMode`.
"""
function nl_constraint_string(model::Model, mode, c::_NonlinearConstraint)
    s = _sense(c)
    nl = nl_expr_string(model, mode, c.terms)
    if s == :range
        out_str =
            "$(_string_round(c.lb)) " *
            _math_symbol(mode, :leq) *
            " $nl " *
            _math_symbol(mode, :leq) *
            " " *
            _string_round(c.ub)
    else
        if s == :<=
            rel = _math_symbol(mode, :leq)
        elseif s == :>=
            rel = _math_symbol(mode, :geq)
        else
            rel = _math_symbol(mode, :eq)
        end
        out_str = string(nl, " ", rel, " ", _string_round(_rhs(c)))
    end
    return out_str
end

#------------------------------------------------------------------------
## Opaque nonlinear objects
#------------------------------------------------------------------------
# TODO: This could be pretty printed.

function function_string(::Type{<:PrintMode}, p::NonlinearExpression)
    return "Reference to nonlinear expression #$(p.index)"
end

function function_string(::Type{<:PrintMode}, p::NonlinearParameter)
    relevant_parameters = filter(
        i -> i[2] isa NonlinearParameter && i[2].index == p.index,
        p.m.obj_dict,
    )
    if length(relevant_parameters) == 1
        par_name = first(relevant_parameters)[1]
        return "Reference to nonlinear parameter $(par_name)"
    else
        return "Reference to nonlinear parameter #$(p.index)"
    end
end

function Base.show(io::IO, ex::Union{NonlinearExpression,NonlinearParameter})
    return Base.show(io, function_string(REPLMode, ex))
end

function Base.show(
    io::IO,
    ::MIME"text/latex",
    ex::Union{NonlinearExpression,NonlinearParameter},
)
    return print(io, function_string(IJuliaMode, ex))
end

# TODO: Print the status of the NLPEvaluator, features available, etc.
function Base.show(io::IO, evaluator::NLPEvaluator)
    return Base.show(io, "A JuMP.NLPEvaluator")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## 08.jl
# 暗号文

function cipher(sentence)
  s = []
  for c in sentence
    islower(c) && push!(s, Char(219 - Int(c)) )
    !islower(c) && push!(s, c)
  end
  return join(s,"")
end

sentence = "Now I need a drink, alcoholic of course,
 after the heavy lectures involving quantum mechanics."

enc = cipher(sentence)
println(enc)
dec = cipher(enc)
println(dec)
println(dec == enc)

ext: jl


text:
module DiffEqMonteCarlo

using DiffEqBase, RecursiveArrayTools, StaticArrays, Distributed, Statistics
import DiffEqBase: MonteCarloSummary

include("solve.jl")
include("analysis.jl")

export MonteThreads, MonteDistributed, MonteSplitThreads, MonteSerial

export get_timestep, get_timepoint, apply_timestep, apply_timepoint,
       componentwise_vectors_timestep, componentwise_vectors_timepoint

export componentwise_mean, componentwise_meanvar

export timestep_mean, timestep_median, timestep_quantile, timestep_meanvar,
       timestep_meancov, timestep_meancor, timestep_weighted_meancov

export timeseries_steps_mean, timeseries_steps_median, timeseries_steps_quantile,
       timeseries_steps_meanvar, timeseries_steps_meancov,
       timeseries_steps_meancor, timeseries_steps_weighted_meancov

export timepoint_mean, timepoint_median, timepoint_quantile,
       timepoint_meanvar, timepoint_meancov,
       timepoint_meancor, timepoint_weighted_meancov

export timeseries_point_mean, timeseries_point_median, timeseries_point_quantile,
       timeseries_point_meanvar, timeseries_point_meancov,
       timeseries_point_meancor, timeseries_point_weighted_meancov

end # module


text:
"""
This module provides a set of helper macros and helper functions

The exported macros are:

$(EXPORTS)

"""
module Helpers
using DocStringExtensions

import Base: +, -, *, /, transpose, adjoint
import LinearAlgebra: cross, tr, dot

export @abstractmethod
export @notimplemented
export @notimplementedif
export @unreachable
export @check
export tfill
export get_val_parameter
export first_and_tail
export GridapType
#export operate

include("Macros.jl")

include("HelperFunctions.jl")

include("GridapTypes.jl")

end # module Helpers


text:
function myfunc(x)
    x1, x2 = x[1], x[2]
    return x1^4 - 2*x2*x1^2 + x2^2 + x1^2 - 2*x1 + 5
end

lb = [-3.0, -1.0]
ub = [2.0, 6.0];

particles = Particles(100, lb, ub)
pso(particles, myfunc, minstep=1e-3)
# ou posso fazer a função criar um e depor retornar
particles2 = pso(myfunc, lb, ub, minstep=1e-3);

mínimo em [1,1] com valor 4


text:
@testset "Events/Processing Instructions" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    DC = E.DataContent
    ME = E.MarkupError
    PI = E.ProcessingInstruction

    @testset "Events/Processing Instructions (Positive)" begin
        @test evaluate("<?x?>")   == [ PI("x", "", L.Location("a buffer", -1)) ]
        @test evaluate("<?xxx?>") == [ PI("xxx", "", L.Location("a buffer", -1)) ]

        @test evaluate("<?x?><?y?>") == [ PI("x", "", L.Location("a buffer", -1)),
                                          PI("y", "", L.Location("a buffer", -1)) ]

        @test evaluate("<?é?>") == [ PI("é", "", L.Location("a buffer", -1)) ]

        # The specification allows for trailing white space after the PI target, but it's considered pure syntax by the
        # parser, so make sure it doesn't appear in the generated event. Try various types and numbers of white space.
        #
        @test evaluate("<?target ?>") == [ PI("target", "", L.Location("a buffer", -1)) ]
        @test evaluate("<?target   ?>") == [ PI("target", "", L.Location("a buffer", -1)) ]
        @test evaluate("<?target \n ?>") == [ PI("target", "", L.Location("a buffer", -1)) ]

        # The specification excludes all case variants of "xml" for the PI target, but the parser doesn't care. So
        # ensure they get captured at this level.
        #
        events = evaluate("<?XML?>")
        @test events == [ ME("ERROR: A PI target cannot be any case variant of 'XML'.", L.Location("a buffer", -1)),
                          PI("XML", "", L.Location("a buffer", -1)) ]

        events = evaluate("<?xml?>")
        @test events == [ ME("ERROR: A PI target cannot be any case variant of 'XML'.", L.Location("a buffer", -1)),
                          PI("xml", "", L.Location("a buffer", -1)) ]

        events = evaluate("<?xmL?>")
        @test events == [ ME("ERROR: A PI target cannot be any case variant of 'XML'.", L.Location("a buffer", -1)),
                          PI("xmL", "", L.Location("a buffer", -1)) ]

        @test evaluate("<?target value?>") == [ PI("target", "value", L.Location("a buffer", -1)) ]
        @test (evaluate("<?target value value value?>") == [ PI("target", "value value value", L.Location("a buffer", -1)) ])
        @test (evaluate("<?target value\nvalue\nvalue?>") == [ PI("target", "value\nvalue\nvalue", L.Location("a buffer", -1)) ])

        @test evaluate("<?target ééé?>") == [ PI("target", "ééé", L.Location("a buffer", -1)) ]

        # The PI value is anything following the first white space, including any white space trailing at the end of the
        # PI. Make sure that any trailing white space is captured.
        #
        @test evaluate("<?target value ?>") == [ PI("target", "value ", L.Location("a buffer", -1)) ]

        # A '?' can appear in a PI value.
        #
        @test evaluate("<?Are you there? Hello, World!?>") == [ PI("Are", "you there? Hello, World!", L.Location("a buffer", -1)) ]
    end

    @testset "Events/Processing Instructions (Negative ... no target)" begin
        # Check that a missing PI target is caught.
        #
        events = evaluate("<??>")
        @test (events == [ ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)),
                           DC("?", false, L.Location("a buffer", -1)),
                           DC(">", false, L.Location("a buffer", -1)) ])

        # Check that EOI is caught.
        #
        events = evaluate("<?")
        @test (events == [ ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)) ])

        # Check that a random token is caught. But careful ... the parser keeps going, so we only check the FIRST
        # event. (Otherwise we're testing the results of some other part of the parser.)
        #
        events = evaluate("<?<")
        @test (events == [ ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])
    end

    @testset "Events/Processing Instructions (Negative ... no terminator)" begin
        # Check that EOI is caught.
        #
        events = evaluate("<?target")
        @test (events == [ ME("ERROR: Expecting '?>' to end a processing instruction.", L.Location("a buffer", -1)) ])

        events = evaluate("<?target value")
        @test (events == [ ME("ERROR: Expecting '?>' to end a processing instruction.", L.Location("a buffer", -1)) ])

        events = evaluate("<?target value ")
        @test (events == [ ME("ERROR: Expecting '?>' to end a processing instruction.", L.Location("a buffer", -1)) ])

        # Check that a random token is caught. Note the random token goes into the PI value.
        #
        events = evaluate("<?target<")
        @test (events == [ ME("ERROR: Expecting '?>' to end a processing instruction.", L.Location("a buffer", -1)) ])
    end
end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
# piecewise homogenous Poisson process

# n intervals - n + 1 bounds
# (b_1, b_2]
# (b_2, b_3]
# ..
# (b_n, b_{n+1}]

function compute_total(bounds, rates)
    num_intervals = length(rates)
    if length(bounds) != num_intervals + 1
        error("Number of bounds does not match number of rates")
    end
    total = 0.
    bounds_ascending = true
    for i=1:num_intervals
        lower = bounds[i]
        upper = bounds[i+1]
        rate = rates[i]
        len = upper - lower
        if len <= 0
            bounds_ascending = false
        end
        total += len * rate
    end
    (total, bounds_ascending)
end

struct PiecewiseHomogenousPoissonProcess <: Distribution{Vector{Float64}} end
const piecewise_poisson_process = PiecewiseHomogenousPoissonProcess()

function Gen.logpdf(::PiecewiseHomogenousPoissonProcess, x::Vector{Float64}, bounds::Vector{Float64}, rates::Vector{Float64})
    cur = 1
    upper = bounds[cur+1]
    lpdf = 0.
    for xi in sort(x)
        if xi < bounds[1] || xi > bounds[end]
            error("x ($xi) lies outside of interval")
        end
        while xi > upper
            cur += 1
            upper = bounds[cur+1]
        end
        lpdf += log(rates[cur])
    end
    (total, bounds_ascending) = compute_total(bounds, rates)
    if bounds_ascending
        lpdf - total
    else
        -Inf
    end
end

function Gen.random(::PiecewiseHomogenousPoissonProcess, bounds::Vector{Float64}, rates::Vector{Float64})
    x = Vector{Float64}()
    num_intervals = length(rates)
    for i=1:num_intervals
        lower = bounds[i]
        upper = bounds[i+1]
        rate = (upper - lower) * rates[i]
        n = random(poisson, rate)
        for j=1:n
            push!(x, random(uniform_continuous, lower, upper))
        end
    end
    x
end


text:
using DelaySSAToolkit
using Test, DiffEqJump, StaticArrays
#σ_off: Gon -> Goff
#σ_on: Goff -> Gon
#ρ_on: Gon -> Gon + N, triggers N => τ 0
#d: P -> 0
# 1. Gon, 2. Goff, 3. P
params = [0.1, 0.1, 10., 0.1, 10., 50.]
σ_off, σ_on, ρ_on, d, τ, tf = params
rates = [σ_off, σ_on, ρ_on,  d]
react_stoich = [[1=>1],[2=>1],[1=>1],[3=>1]]
net_stoich = [[1=>-1,2=>1],[1=>1,2=>-1],[3=>1],[3=>-1]]
mass_action_jump = MassActionJump(rates, react_stoich, net_stoich; scale_rates = false)
jumpset = JumpSet((),(),nothing,mass_action_jump)
delay_trigger = Dict(3=>[1=>τ])
delay_complete = Dict(1=>[3=>-1])
delay_affect! = function (integrator, rng)
    i = rand(rng, 1:length(integrator.de_chan[1]))
    deleteat!(integrator.de_chan[1],i)
end
delay_interrupt = Dict(4=>delay_affect!)

delayjumpset = DelayJumpSet(delay_trigger,delay_complete,delay_interrupt)
u0 =@SVector [0,1,0]
de_chan0 = [[]]
tspan = (0.,tf)
dprob = DiscreteProblem(u0, tspan)

djprob = DelayJumpProblem(dprob, DelayDirect(), jumpset, delayjumpset, de_chan0, save_positions = (true,true),save_delay_channel = true)
sol = @benchmark solve(djprob, SSAStepper())

algo_list =[DelayDirect(),DelayRejection(),DelayMNRM(),DelayDirectCR()]
meanlist=[]
medianlist=[]
minlist=[]
maxlist=[]
stdlist=[]
algo_list[1]

djprob = DelayJumpProblem(dprob, DelayMNRM(), jumpset, delayjumpset, de_chan0, save_positions = (true,true),save_delay_channel = true)
    a=@benchmark solve(djprob, SSAStepper())
    push!(meanlist,copy(mean(a).time/1e9))
    push!(medianlist,copy(median(a).time/1e9))
    push!(minlist,copy(minimum(a).time/1e9))
    push!(maxlist,copy(maximum(a).time/1e9))
    push!(stdlist,copy(std(a).time/1e9))
    print(djprob.aggregator," OK","\n")

1:length(integrator.de_chan[1])
1:1
length([])
@time testalgo(algo_list)

meanlist
medianlist
minlist
maxlist
stdlist
df=DataFrame(mean=meanlist,median=medianlist,min=minlist,max=maxlist,std=stdlist)
CSV.write("C:/Users/86158/Desktop/algotest/save_delay_channel_save_positions=T.csv",df)


text:
using Distributed
using TimerOutputs
using SparseArrays

to = TimerOutput()

if nprocs() == 1
    addprocs(3)
end
# Number of logical CPU cores available in the system.
# println("CPU cores: ", Sys.CPU_CORES)
# Number of available processes.
# Sys.cpu_summary()
println("nprocs: ", nprocs())
println("nworkers: ", nworkers())
println("positive numbers are the faces, negative simulate the end of block")

@everywhere fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)

ch = RemoteChannel(()->Channel{Int}(32));
@everywhere function process_block(channel, i)
    print(".")
    for j=1:rand(1:4)
        # simulate extracting usefull faces from the brick
        sleep(rand())
        # fib(rand(10:40))
        put!(channel, i*10 + j)
    end
    put!(channel, -i)
end

number_of_blocks = 10

@distributed for i = 1:number_of_blocks
    process_block(ch, i)
end

Fchar = spzeros(Int64, 1000)
# Read all the values
n = 0
while n < number_of_blocks
    global n
    val = take!(ch)
    print("$val ")
    if val < 0
        n += 1
    else
        Fchar[val] = 1
    end
end


text:
abstract type AbstractData{N,T<:Real} end

mutable struct Data{N,T<:Real} <: AbstractData{N,T}
	pos::SVector{N,T}
	idx::Int64
	hsml::T
	mass::T
end

Data{N,T}() where {N,T} = Data{N,T}(zero(SVector{N,T}),0,0,0)

function assign_node_data!(n::D, old::D, new::D) where {D<:AbstractData}
	masssum = old.mass + new.mass
	inv_masssum = 1.0 / masssum
	n.pos  = (old.mass .* old.pos .+ new.mass .* new.pos) .* inv_masssum
	n.mass = masssum #this line has to be put after insertpart!, otherwise we're inserting back the wrong node.mass!!!
	n.hsml = old.hsml > new.hsml ? old.hsml : new.hsml
	assign_additional_node_data!(n, old, new)
end

function assign_additional_node_data!(n::Data, old::Data, new::Data)
end


text:
module V
const NAMES = [
    "pERK"
    "pERK_DUSP"
    "pRSK"
    "cFos"
    "Cd"
    "uRb"
    "tE2f"
    "E2f"
    "tCe"
    "CeP21"
]

#name2idx
for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module


text:
module Example

end # module


max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Base.Threads
using Base.Threads.Atomic
import Base: push!, shift!, isempty, getindex

using UnsafeAtomics

export ThreadQueue, LevelSynchronous, bfs

"""
A Queue with the head and tail guarded by Atomics.
"""
immutable ThreadQueue{T}
    data::Vector{T}
    head::Atomic{Int}
    tail::Atomic{Int}
end

function ThreadQueue(T::Type, maxlength::Int)
    q = ThreadQueue(Vector{T}(maxlength), Atomic{Int}(1), Atomic{Int}(1))
    return q
end

function push!{T}(q::ThreadQueue{T}, val::T)
    # TODO: check that head > tail
    offset = atomic_add!(q.tail, 1)
    q.data[offset] = val
    return offset
end

function shift!{T}(q::ThreadQueue{T})
    # TODO: check that head < tail
    offset = atomic_add!(q.head, 1)
    return q.data[offset]
end

function isempty(q::ThreadQueue)
    return ( q.head[] == q.tail[] ) && q.head != 1
    # return q.head == length(q.data)
end

function getindex{T}(q::ThreadQueue{T}, iter)
    return q.data[iter]
end

abstract type BFSAlgorithm end
type LevelSynchronous <: BFSAlgorithm end

function bfskernel(
        alg::LevelSynchronous, s::Stinger, next::ThreadQueue, parents::Array{Int64},
        level::Array{Int64}
    )
    @threads for src in level
        foralledges(s, src) do edge, src, etype
            direction, neighbor = edgeparse(edge)
            if (direction != 1)
                parent = UnsafeAtomics.unsafe_atomic_cas!(parents, neighbor+1, -2, src)
                if parent==-2
                    push!(next, neighbor) #Push onto queue
                end
            end
        end
    end
end

"""
    bfs(alg::LevelSynchronous, s::Stinger, source::Int64, nv::Int64)

Uses the `LevelSynchronous` algorithm to perform a parallel BFS with the
specified source node.
"""
function bfs(alg::LevelSynchronous, s::Stinger, source::Int64, nv::Int64)
    next = ThreadQueue(Int, nv)
    parents = fill(-2, nv)
    bfs(alg, s, next, source, parents)
end

function bfs(
        alg::LevelSynchronous, s::Stinger, next::ThreadQueue, source::Int64,
        parents::Array{Int64}
    )
    parents[source+1]=-1 #Set source to -1
    push!(next, source)
    while !isempty(next)
        level = next[next.head[]:next.tail[]-1]
        next.head[] = next.tail[] #reset the queue
        bfskernel(alg, s, next, parents, level)
    end
    return parents
end


text:
macro fastcall(expr)
    if !(expr isa Expr) || expr.head != :call || length(expr.args) < 1 || !(expr.args[1] isa Symbol)
        error("@fastcall macro expects an expression like f(x, y, ...; a=A, b=B, ...). Got $expr")
    end

    fname = expr.args[1]

    simplecall = true
    args = expr.args[2:end]
    for arg ∈ args
        if iskw(arg)
            simplecall = false
            break
        end
    end

    if simplecall
        return esc(expr)
    else
        newargs = make_ordered_kwlist(args)
        return Expr(:call, esc(fname), Expr(:tuple, newargs...))
    end
end

function make_ordered_kwlist(args)
    std_args = Vector{Any}()
    kw_args = Vector{Any}()
    kwlist = Vector{Symbol}()
    for arg ∈ args
        if iskw(arg)
            push!(kw_args, :(@kw($(esc(arg.args[1])) = $(esc(arg.args[2])))))
            push!(kwlist, arg.args[1])
        else
            push!(std_args, esc(arg))
        end
    end

    perm = sortperm(kwlist)
    kw_args = kw_args[perm]

    return (Expr(:tuple, KWSentinal(), kw_args...), std_args...)
end


text:
using Genie.Configuration, Logging
server_sett = TOML.parsefile("config/env/settings.toml")
port = server_sett["port"]
host = server_sett["host"]
print(string("https://", host, ":", port))
const config = Settings(
  server_port                     = parse(Int64, port),
  server_host                     = host,
  log_level                       = Logging.Info,
  log_to_file                     = false,
  server_handle_static_files      = true,
  path_build                      = "build",
  format_julia_builds             = true,
  format_html_output              = true
)

ENV["JULIA_REVISE"] = "auto"


text:
import Pkg

@inline function _with_temp_dir(f::F) where F <: Base.Callable
    original_directory = pwd()
    my_tmp_dir = mktempdir(; cleanup=false)
    cd(my_tmp_dir)
    result = f(my_tmp_dir)
    cd(original_directory)
    rm(my_tmp_dir; force = true, recursive = true)
    return result
end

# @inline function _with_temp_depot(f::F) where F <: Base.Callable
#     return _with_temp_dir() do my_temp_depot
#         original_depot_path = deepcopy(Base.DEPOT_PATH)
#         empty!(Base.DEPOT_PATH)
#         push!(Base.DEPOT_PATH, my_temp_depot)
#         result = f(my_temp_depot)
#         empty!(Base.DEPOT_PATH)
#         for x in original_depot_path
#             push!(Base.DEPOT_PATH, x)
#         end
#         return result
#     end
# end

@inline function _with_temp_project(f::F) where F <: Base.Callable
    return _with_temp_dir() do my_temp_project
        original_active_project = Base.active_project()
        Pkg.activate(my_temp_project; shared = false)
        result = f(my_temp_project)
        Pkg.activate(original_active_project; shared = false)
        return result
    end
end

@inline function _activate_possibly_shared_project(project::AbstractString;
                                                   shared_project::Bool = false)
    if length(project) > 0 && project[1] == '@'
        Pkg.activate(project[2:end]; shared = true)
    else
        Pkg.activate(project; shared = shared_project)
    end
    Pkg.instantiate()
    return nothing
end

@inline function _with_different_project(f::F,
                                         project::AbstractString;
                                         shared_project::Bool = false) where F <: Base.Callable
    original_active_project = Base.active_project()
    _activate_possibly_shared_project(project; shared_project = shared_project)
    result = f(Base.active_project())
    Pkg.activate(original_active_project; shared = false)
    return result
end


text:
function ghactions(; kwargs...)
    username = ENV["COMPAT_USERNAME"]
    spec = RepoSpec(
        reponame = ENV["GITHUB_REPOSITORY"],
        username = username,
        token = ENV["COMPAT_TOKEN"],
        masterbranch = get(ENV, "COMPAT_MASTER_BRANCH", MASTER_BRANCH),
        compatbranch = get(ENV, "COMPAT_COMPAT_BRANCH", COMPAT_BRANCH),
        gitconfig = Dict{String,String}("user.name" => username),
    )
    update_tomls!(spec; kwargs...)
    return nothing
end

function update_tomls!(rspec::RepoSpec; update_from_test::Bool = true, overwrite::Bool = true, kwargs...)
    masterbranch = rspec.masterbranch
    compatbranch = rspec.compatbranch

    gitcmd = create_git_cmd(rspec.gitconfig)
    @debug "Git command: $gitcmd"

    with_tempdir() do _
        auth = authenticate(rspec.token)
        ghrepo = repo(rspec.reponame, auth = auth)

        url = gitauthurl(rspec.username, rspec.token, rspec.reponame)
        run(`$gitcmd clone $url repo`)
        cd("repo")

        run(`$gitcmd checkout $masterbranch`)
        run(`$gitcmd pull`)
        if overwrite
            run(`$gitcmd checkout -B $compatbranch`)
        else
            run(`$gitcmd checkout -b $compatbranch`)
        end

        if update_from_test
            rs, msg = update_tomls_from_test!(pwd(); kwargs...)
        else
            rs, msg = update_tomls!(pwd(), kwargs...)
        end

        for r in rs
            project_updated(r) && safe_git_add(gitcmd, r.project_file)
            manifest_updated(r) && safe_git_add(gitcmd, manifestfile_path(dirname(r.project_file)))
        end

        if success(`$gitcmd diff-index --quiet --cached HEAD --`)
            @info "No changes"
        else
            title = "New compat entries"
            body = string(msg)

            run(`$gitcmd commit -m $title`)
            run(`$gitcmd push --force -u origin $compatbranch`)

            existing_pr = find_existing_pr(ghrepo, auth, compatbranch, masterbranch)
            if existing_pr === nothing
                params = Dict(
                    "title" => title,
                    "body" => body,
                    "labels" => PR_LABELS,
                    "head" => compatbranch,
                    "base" => masterbranch,
                )
                pr = create_pull_request(ghrepo; params = params, auth = auth)
                # NOTE: Adding labels when PR is created doesn't appear to work so
                # we edit it after creation.
                edit_issue(ghrepo, pr, params = Dict("labels" => PR_LABELS), auth = auth)
            else
                params = Dict("title" => title, "body" => body, "labels" => PR_LABELS)
                edit_issue(ghrepo, existing_pr, params = params, auth = auth)
            end
        end

        return rs
    end
end

function safe_git_add(gitcmd, path)
    if success(`$gitcmd check-ignore $path`)
        # path is in .gitignore
        @warn "$path is ignored by git. Skipping."
    else
        run(`$gitcmd add $path`)
    end
end

function find_existing_pr(repo::Repo, auth::Authorization, head::String, base::String)
    params = Dict(
        "state" => "open",
        "head" => "$(repo.owner.login):$head",
        "base" => base,
        "per_page" => 50,
        "page" => 1,
    )
    prs = Vector{PullRequest}()

    page_prs, page_data = pull_requests(repo, params = params, auth = auth)
    while true
        append!(prs, page_prs)
        if haskey(page_data, "next")
            page_prs, page_data = pull_requests(repo, auth = auth, start_page = page_data["next"])
        else
            break
        end
    end

    if isempty(prs)
        @debug "No PR found"
        return nothing
    elseif length(prs) == 1
        pr = first(prs)
        if occursin(COMPAT_UUID, pr.body)
            @info "Found existing PR"
            return pr
        else
            error("PR already exists but it wasn't created by Compat. Aborting.")
        end
    else
        error("More than one compat PR found. This shouldn't happen.")
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using DocOpt

#CompilerTools.OptFramework.set_debug_level(3)
#CompilerTools.LambdaHandling.set_debug_level(3)
#ParallelAccelerator.set_debug_level(3)
#ParallelAccelerator.DomainIR.set_debug_level(3)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(3)

@acc function compute(src, dst, N)
  runStencil(dst, src, N, :oob_skip) do b, a
    b[0,0,0] = (a[0,0,-1] + a[0,0,1] + a[0,-1,0] + a[0,1,0] + a[-1,0,0] + a[1,0,0]) / 6
    return a, b
  end
  return isodd(N) ? dst : src
end

function initialize(x, y, z)
  src = Array{Float32}(x, y, z)
  rand!(src)
  dst = copy(src)
  return src, dst
end

function checksum(src)
  local z, y, x
  (x, y, z) = size(src)
  sum(src) / (x * y * z)
end

function main()
    doc = """laplace-3d.jl

Laplace 6-point 3D stencil.

Usage:
  laplace-3d.jl -h | --help
  laplace-3d.jl [--size=<size>] [--iterations=<iterations>]

Options:
  -h --help                  Show this screen.
  --size=<size>              Specify a 3d array size (<size> x <size> x <size>) [default: 300].
  --iterations=<iterations>  Specify a number of iterations [default: 100].
"""
    arguments = docopt(doc)

    size = parse(Int, arguments["--size"])
    iterations = parse(Int, arguments["--iterations"])

    srand(0)

    function laplace_3d(iterations)
        local src, dst
        println("Run laplace-3d with size ", size, "x", size, "x", size, " for ", iterations, " iterations.")
        (src, dst) = initialize(size, size, size)
        tic()
        compute(src, dst, 0)
        println("SELFPRIMED ", toq())

        tic()
        dst = compute(src, dst, iterations)
        println("SELFTIMED ", toq())
        println("checksum: ", checksum(dst))
    end

    laplace_3d(iterations)

end

main()


text:
"""
$(TYPEDEF)
Optimizer for finding a sparse implicit relationship via alternating the left hand side of the problem and
solving the explicit problem, as introduced [here](https://royalsocietypublishing.org/doi/10.1098/rspa.2020.0279).

```math
\\argmin_{x} \\|x\\|_0 ~s.t.~Ax= 0
```

# Fields
$(FIELDS)

# Example
```julia
ImplicitOptimizer(STLSQ())
ImplicitOptimizer(0.1f0, ADMM)
```
"""
mutable struct ImplicitOptimizer{T} <: AbstractSubspaceOptimizer{T}
    """Explicit Optimizer"""
    o::AbstractOptimizer{T}

    function ImplicitOptimizer(threshold = 1e-1, opt = STLSQ)
        return new{typeof(threshold)}(opt(threshold))
    end

    function ImplicitOptimizer(opt::AbstractOptimizer{T}) where T
        return new{T}(opt)
    end

end

Base.summary(opt::ImplicitOptimizer) = "Implicit Optimizer using "*summary(opt.o)

get_threshold(opt::ImplicitOptimizer) = get_threshold(opt.o)

function (opt::ImplicitOptimizer{T})(X, A, Y, λ::V = first(opt.o.λ);
    maxiter::Int64 = maximum(size(A)), abstol::V = eps(eltype(T)),
    rtol::V = zero(eltype(T)) ,progress = nothing,
    f::Function = F(opt),
    g::Function = G(opt),
    scale_coefficients::Bool = false)  where {T, V}

    exopt = opt.o

    n,m = size(A)
    ny, my = size(Y)
    nx, mx = size(X)
    nq, mq = 0,0

    # Closure for the pareto function
    fg(x, A, y) = (g∘f)(x, A, y)
    fg(x, A) = (g∘f)(x,A)

    xzero = zero(eltype(X))
    xone = one(eltype(X))

    obj = xzero
    sparsity = xzero
    conv_measure = xzero

    iters = 0
    converged = false

    # Build a quadratic matrix
    x_tmp = zeros(eltype(X), m, m) # The temporary solution
    x_opt = zeros(eltype(X), m, m) # All solutions
    inds = [false for _ in 1:m]

    _progress = isa(progress, Progress)

    # Set progress -1 and
    initial_prog = _progress ? progress.counter : 0

    # Iterate over all columns of A ( which represent the lhs )
    @views for j in 1:m
        inds .= true
        inds[j] = false
        x_tmp[j, j] = -one(eltype(X)) # We set this to -1 for the lhs
        # Solve explicit problem
        x_tmp[inds, j:j] .= init(exopt, A[:, inds], A[:, j:j])

        # Use optimizer
        λs = sparse_regression!(x_tmp[inds, j:j], A[:, inds], A[:, j:j],exopt, maxiter = maxiter, abstol = abstol)

        # Normalize
        if scale_coefficients
            x_tmp[inds,j:j] ./= maximum(abs.(x_tmp[inds, j:j]))
        end

        if _progress
            sparsity, obj = f(x_tmp[inds, :], A[:, inds], A[:, j:j], λs[1])

            ProgressMeter.next!(
            progress;
            showvalues = [
                (:Threshold, λ), (:Objective, obj), (:Sparsity, sparsity)
                ]
                )
        end
    end

    # Reduce the solution size to linear independent columns
    @views x_tmp = linear_independent_columns(x_tmp, rtol)

    # Indicate if already used
    _included = zeros(Bool, my, size(x_tmp, 2))
    for i in 1:my, j in 1:size(x_tmp, 2)
        # Check, if already included
        any(_included[:, j]) && continue
        # Selector
        inds .= true; inds[j] = false

        if fg(x_tmp[inds, j], A[:, inds], A[:, j]) < fg(X[inds, i], A[:, inds], A[:, j])
            X[:, i] .= x_tmp[:, i]
            _included[i,j] = true
        end

    end

    if rank(X'X) < my
        @warn "$opt @ $λ has found illconditioned equations. Vary the threshold or relative tolerance."
    end

    return
end


text:
module Solution
  chromosomeCheck(sperm) =
    'Y' ∈ sperm ? "Congratulations! You're going to have a son." : "Congratulations! You're going to have a daughter."
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using Test
using ThreeBodyDecay

@testset "Sum over polarization" begin
    mΛb = 5.61960
    mπ = 0.14
    # intermediate resonances
    mΣb = 5.81065; ΓΣb = 0.005
    mΣb_x = 5.83032; ΓΣb_x = 0.0094

    mΛb2S = 6.3; # just a peak of the plot

    tbs = ThreeBodySystem(mπ,mΛb,mπ; m0=mΛb2S,
        two_js = ThreeBodySpins(0, 1, 0; two_h0=1))  # 1/2+ 0- 0- 1/2+
    dpp = randomPoint(tbs)
    σs = dpp.σs

    # lineshape
    dc_Σb1 = DecayChainLS(1, σ->BW(σ,mΣb,ΓΣb); tbs=tbs, two_s = 1, parity='+', Ps=['-','+','-','+'])
    dc_Σb3 = DecayChainLS(3, σ->BW(σ,mΣb,ΓΣb); tbs=tbs, two_s = 1, parity='+', Ps=['-','+','-','+'])

    dc_Σb_x1 = DecayChainLS(1, σ->BW(σ,mΣb_x,ΓΣb_x); tbs=tbs, two_s = 3, parity='+', Ps=['-','+','-','+'])
    dc_Σb_x3 = DecayChainLS(3, σ->BW(σ,mΣb_x,ΓΣb_x); tbs=tbs, two_s = 3, parity='+', Ps=['-','+','-','+'])

    full_a(σs,two_λs) = sum(amplitude(σs,two_λs, ch) for ch in [dc_Σb1, dc_Σb_x1, dc_Σb3, dc_Σb_x3])
    full_a(dpp) = full_a(dpp.σs, dpp.two_λs)
    #
    total_I(σs) = sum(abs2(full_a(σs,two_λs)) for two_λs in itr(tbs.two_js))

    test_I = summed_over_polarization((σs,two_λs)->abs2(full_a(σs,two_λs)),tbs.two_js)
    @test test_I(dpp.σs) == total_I(dpp.σs)
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
using Base.Iterators
using Test
using Arbitrary

using MPSIMD

# Define type
const IntN128 = IntN{Float64, Int32, 4}
const maxval = big(2)^120       # don't go near overflow yet
const maxval2 = typeof(maxval)(sqrt(maxval))

# Generate arbitrary values
isinbounds(x) = abs(x) <= maxval
xs = collect(take(Base.Iterators.filter(isinbounds, arbitrary(BigInt)), 100))
ys = collect(take(Base.Iterators.filter(isinbounds, arbitrary(BigInt)), 100))
zs = collect(take(Base.Iterators.filter(isinbounds, arbitrary(BigInt)), 100))

@testset "Internal representation" begin
    z = IntN128()
    @test Tuple(z.digits) == (0, 0, 0, 0)
    i0 = IntN128(Int32(0))
    @test Tuple(i0.digits) == (0, 0, 0, 0)
    i1 = IntN128(Int32(1))
    @test Tuple(i1.digits) == (1, 0, 0, 0)
    @test Tuple(normalize(i1).digits) == (1, 0, 0, 0)
    im1 = IntN128(Int32(-1))
    @test Tuple(im1.digits) == (-1, 0, 0, 0)
    @test Tuple(normalize(im1).digits) == (-1, 0, 0, 0)
end

@testset "Conversions" begin
    @test BigInt(zero(IntN128)) == 0
    @test BigInt(one(IntN128)) == 1
    for x in xs
        @test BigInt(IntN128(x)) == x
    end
end

@testset "Comparisons" begin
    for x in xs
        @test signbit(IntN128(x)) === signbit(x)
        @test BigInt(abs(IntN128(x))) == abs(x)
        @test cmp(IntN128(x), IntN128(x)) === 0
        @test (IntN128(x) == IntN128(x)) === true
        @test (IntN128(x) != IntN128(x)) === false
        @test (IntN128(x) < IntN128(x)) === false
        @test (IntN128(x) <= IntN128(x)) === true
        @test (IntN128(x) > IntN128(x)) === false
        @test (IntN128(x) >= IntN128(x)) === true
        @test isequal(IntN128(x), IntN128(x)) === true
        @test isless(IntN128(x), IntN128(x)) === false
    end
    for (x,y) in zip(xs, ys)
        @test cmp(IntN128(x), IntN128(y)) === cmp(x, y)
        @test (IntN128(x) == IntN128(y)) === (x == y)
        @test (IntN128(x) != IntN128(y)) === (x != y)
        @test (IntN128(x) < IntN128(y)) === (x < y)
        @test (IntN128(x) <= IntN128(y)) === (x <= y)
        @test (IntN128(x) > IntN128(y)) === (x > y)
        @test (IntN128(x) >= IntN128(y)) === (x >= y)
        @test isequal(IntN128(x), IntN128(y)) == isequal(x, y)
        @test isless(IntN128(x), IntN128(y)) == isless(x, y)
    end
end

@testset "Addition and subtraction" begin
    for x in xs
        @test BigInt(+IntN128(x)) == +x
        @test BigInt(-IntN128(x)) == -x
    end
    for (x,y) in zip(xs, ys)
        @test BigInt(IntN128(x) + IntN128(y)) == x + y
        @test BigInt(IntN128(x) - IntN128(y)) == x - y
    end
    s = IntN128(0)
    for x in xs
        s += IntN128(x)
    end
    @test BigInt(s) == sum(xs)
end

@testset "Scaling" begin
    for (x,y) in zip(xs, ys)
        y1 = y % Int32
        @test BigInt(IntN128(x) * y1) == x * y1
    end
end

@testset "Multiplication" begin
    for (x,y) in zip(xs, ys)
        if abs(x * y) <= maxval
            @test BigInt(IntN128(x) * IntN128(y)) == x * y
        end
        x1 = x % maxval2
        y1 = y % maxval2
        @test BigInt(IntN128(x1) * IntN128(y1)) == x1 * y1
        if x != 0
            x2 = x
            y2 = y ÷ x
            @test BigInt(IntN128(x2) * IntN128(y2)) == x2 * y2
        end
        if x1 != 0
            x3 = x1
            y3 = y ÷ x1
            @test BigInt(IntN128(x3) * IntN128(y3)) == x3 * y3
        end
    end
end

# +z
# -z
#
# i0 = zero(IntN128)
# i1 = IntN128(1)
#
# @test i0 === normalize(i0)
# @test i1 === normalize(i1)
#
# @show normalize(i0)
# @show normalize(i1)
# @show normalize(-i1)
#
# @test i1 - i1 === i0
# @test i1 + -i1 === i0
#
# using InteractiveUtils
# # @code_native MPSIMD.shift_up(i1.digits)
# # @code_native normalize1(i1)
# @code_native i0 + i0

using BenchmarkTools

@benchmark x * y setup=((x,y)=rand(Int16,2)*big(2)^(128÷2-32))
@benchmark x * y setup=((x,y)=rand(Int16,2)*big(2)^(256÷2-32))
@benchmark x * y setup=((x,y)=rand(Int16,2)*big(2)^(512÷2-32))
@benchmark x * y setup=((x,y)=rand(Int16,2)*big(2)^(1024÷2-32))

using BitIntegers

@benchmark Int128(x) * Int128(y) setup=((x,y)=rand(Int,2))
@benchmark Int256(x) * Int256(y) setup=((x,y)=rand(Int,2))
@benchmark Int512(x) * Int512(y) setup=((x,y)=rand(Int,2))
@benchmark Int1024(x) * Int1024(y) setup=((x,y)=rand(Int,2))

const IntN256 = IntN{Float64, Int32, 8}
const IntN512 = IntN{Float64, Int32, 16}
const IntN1024 = IntN{Float64, Int32, 32}

@benchmark IntN128(x) * IntN128(y) setup=((x,y)=rand(Int,2))
# @benchmark IntN256(x) * IntN256(y) setup=((x,y)=rand(Int,2))
# @benchmark IntN512(x) * IntN512(y) setup=((x,y)=rand(Int,2))
# @benchmark IntN1024(x) * IntN1024(y) setup=((x,y)=rand(Int,2))


text:
### FlatS2 types
# (spin-2 fields are just special FieldTuples)
const FlatQUMap{P,T,M}     = FieldTuple{QUMap,     NamedTuple{(:Q,:U),NTuple{2,FlatMap{P,T,M}}}, T}
const FlatQUFourier{P,T,M} = FieldTuple{QUFourier, NamedTuple{(:Q,:U),NTuple{2,FlatFourier{P,T,M}}}, Complex{T}}
const FlatEBMap{P,T,M}     = FieldTuple{EBMap,     NamedTuple{(:E,:B),NTuple{2,FlatMap{P,T,M}}}, T}
const FlatEBFourier{P,T,M} = FieldTuple{EBFourier, NamedTuple{(:E,:B),NTuple{2,FlatFourier{P,T,M}}}, Complex{T}}
# some handy Unions
const FlatS2{P,T,M}=Union{FlatEBMap{P,T,M},FlatEBFourier{P,T,M},FlatQUMap{P,T,M},FlatQUFourier{P,T,M}}
const FlatQU{P,T,M}=Union{FlatQUMap{P,T,M},FlatQUFourier{P,T,M}}
const FlatEB{P,T,M}=Union{FlatEBMap{P,T,M},FlatEBFourier{P,T,M}}
const FlatS2Map{P,T,M}=Union{FlatQUMap{P,T,M},FlatEBMap{P,T,M}}
const FlatS2Fourier{P,T,M}=Union{FlatQUFourier{P,T,M},FlatEBFourier{P,T,M}}

spin(::Type{<:FlatS2}) = S2

### convenience constructors
for (F,F0,(X,Y),T) in [
        (:FlatQUMap,     :FlatMap,     (:Qx,:Ux), :T),
        (:FlatQUFourier, :FlatFourier, (:Ql,:Ul), :(Complex{T})),
        (:FlatEBMap,     :FlatMap,     (:Ex,:Bx), :T),
        (:FlatEBFourier, :FlatFourier, (:El,:Bl), :(Complex{T}))
    ]
    doc = """
        # main constructor:
        $F(
            $X::AbstractArray, $Y::AbstractArray; $(F0==:FlatFourier ? "\n        Nside, # required, size of the map in pixels" : "")
            θpix,  # optional, resolution in arcmin (default: 1)
            ∂mode, # optional, fourier∂ or map∂ (default: fourier∂)
        )

        # more low-level:
        $F{P}($X::AbstractArray, $Y::AbstractArray) # specify pixelization P explicilty
        $F{P,T}($X::AbstractArray, $Y::AbstractArray) # additionally, convert elements to type $T
        $F{P,T,M<:AbstractArray{$T}}($X::M, $Y::M) # specify everything explicilty

    Construct a `$F` object. The top form of the constructor is most convenient
    for interactive work, while the others may be more useful for low-level code.
    """
    @eval begin
        @doc $doc $F
        $F($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array; kwargs...) =
            $F(($F0($X; kwargs...), $F0($Y; kwargs...)))
        $F{P}($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array) where {P} =
            $F{P}(($F0{P}($X), $F0{P}($Y)))
        $F{P,T}($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array) where {P,T} =
            $F{P,T}($F0{P,T}($X), $F0{P,T}($Y))
        $F{P,T,M}($X::AbstractRank2or3Array, $Y::AbstractRank2or3Array) where {P,T,M} =
            $F{P,T,M}($F0{P,T,M}($X), $F0{P,T,M}($Y))
    end
end

### properties
function propertynames(f::FlatS2)
    (:fs, propertynames(f.fs)...,
     (Symbol(string(k,(f isa FlatMap ? "x" : "l"))) for (k,f) in pairs(f.fs))...)
end
getproperty(f::FlatQUMap,     ::Val{:Qx}) = getfield(f,:fs).Q.Ix
getproperty(f::FlatQUMap,     ::Val{:Ux}) = getfield(f,:fs).U.Ix
getproperty(f::FlatQUFourier, ::Val{:Ql}) = getfield(f,:fs).Q.Il
getproperty(f::FlatQUFourier, ::Val{:Ul}) = getfield(f,:fs).U.Il
getproperty(f::FlatEBMap,     ::Val{:Ex}) = getfield(f,:fs).E.Ix
getproperty(f::FlatEBMap,     ::Val{:Bx}) = getfield(f,:fs).B.Ix
getproperty(f::FlatEBFourier, ::Val{:El}) = getfield(f,:fs).E.Il
getproperty(f::FlatEBFourier, ::Val{:Bl}) = getfield(f,:fs).B.Il
getproperty(f::FlatS2,        ::Val{:P})  = f
function getindex(f::FlatS2, k::Symbol; full_plane=false)
    maybe_unfold = (full_plane && k in [:El,:Bl,:Ql,:Ul]) ? x->unfold(x,fieldinfo(f).Ny) : identity
    B = @match k begin
        (:P)         => identity
        (:E  || :B)  => @match f begin
            _ :: FlatQUMap     => EBMap
            _ :: FlatQUFourier => EBFourier
            _                  => identity
        end
        (:Q  || :U)  => @match f begin
            _ :: FlatEBMap     => QUMap
            _ :: FlatEBFourier => QUFourier
            _                  => identity
        end
        (:Ex || :Bx) => EBMap
        (:El || :Bl) => EBFourier
        (:Qx || :Ux) => QUMap
        (:Ql || :Ul) => QUFourier
        _ => throw(ArgumentError("Invalid FlatS2 index: $k"))
    end
    maybe_unfold(getproperty(B(f),k))
end

function getindex(D::DiagOp{<:FlatEBFourier}, k::Symbol)
    @unpack El, Bl = diag(D)
    @unpack sin2ϕ, cos2ϕ, P = fieldinfo(diag(D))
    f = @match k begin
        (:QQ)        => FlatFourier{P}(@. Bl*sin2ϕ^2 + El*cos2ϕ^2)
        (:QU || :UQ) => FlatFourier{P}(@. (El - Bl) * sin2ϕ * cos2ϕ)
        (:UU)        => FlatFourier{P}(@. Bl*cos2ϕ^2 + El*sin2ϕ^2)
        _            => getproperty(D.diag, k)
    end
    Diagonal(f)
end

### basis conversion

QUFourier(f::FlatQUMap) = FlatQUFourier(map(Fourier,f.fs))
QUFourier(f::FlatEBMap) = f |> EBFourier |> QUFourier
QUFourier(f::FlatEBFourier{P,T}) where {P,T} = begin
    @unpack sin2ϕ, cos2ϕ = fieldinfo(f)
    Ql = @. - f.El * cos2ϕ + f.Bl * sin2ϕ
    Ul = @. - f.El * sin2ϕ - f.Bl * cos2ϕ
    FlatQUFourier(Q=FlatFourier{P}(Ql), U=FlatFourier{P}(Ul))
end

QUMap(f::FlatQUFourier)  = FlatQUMap(map(Map,f.fs))
QUMap(f::FlatEBMap)      = f |> EBFourier |> QUFourier |> QUMap
QUMap(f::FlatEBFourier)  = f |> QUFourier |> QUMap

EBFourier(f::FlatEBMap) = FlatEBFourier(map(Fourier,f.fs))
EBFourier(f::FlatQUMap) = f |> QUFourier |> EBFourier
EBFourier(f::FlatQUFourier{P,T}) where {P,T} = begin
    @unpack sin2ϕ, cos2ϕ = fieldinfo(f)
    El = @. - f.Ql * cos2ϕ - f.Ul * sin2ϕ
    Bl = @.   f.Ql * sin2ϕ - f.Ul * cos2ϕ
    FlatEBFourier(E=FlatFourier{P}(El), B=FlatFourier{P}(Bl))
end

EBMap(f::FlatEBFourier) = FlatEBMap(map(Map,f.fs))
EBMap(f::FlatQUMap)     = f |> QUFourier |> EBFourier |> EBMap
EBMap(f::FlatQUFourier) = f |> EBFourier |> EBMap

QUFourier(f′::FlatQUFourier, f::FlatQUMap) = (map(Fourier,f′.fs,f.fs); f′)
QUMap(f′::FlatQUMap, f::FlatQUFourier) = (map(Map,f′.fs,f.fs); f′)

Map(f::FlatQUFourier) = QUMap(f)
Map(f::FlatEBFourier) = EBMap(f)
Map(f::FlatS2Map) = f
Fourier(f::FlatQUMap) = QUFourier(f)
Fourier(f::FlatEBMap) = EBFourier(f)
Fourier(f::FlatS2Fourier) = f

### simulation and power spectra
function white_noise(rng::AbstractRNG, ::Type{F2}) where {F2<:FlatS2}
    F = (((::Type{<:FieldTuple{B,NamedTuple{Names,NTuple{2,F}}}}) where {B,Names,F}) -> F)(F2)
    FlatEBMap(E=white_noise(rng,F), B=white_noise(rng,F))
end
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S2}, CℓEE::InterpolatedCℓs, CℓBB::InterpolatedCℓs; kwargs...) where {P,T,M}
    Diagonal(FlatEBFourier(E=diag(Cℓ_to_Cov(P,T,S0,CℓEE;kwargs...)), B=diag(Cℓ_to_Cov(P,T,S0,CℓBB;kwargs...))))
end
function Cℓ_to_Cov(::Type{P}, ::Type{T}, ::Type{S2}, (CℓEE, ℓedges, θname)::Tuple, CℓBB::InterpolatedCℓs; units=fieldinfo(P).Ωpix) where {P,T}
    C₀ = Cℓ_to_Cov(P, T, S2, CℓEE, CℓBB, units=units)
    Cbins = Diagonal.(FlatEBFourier.(MidPasses(ℓedges) .* [diag(C₀).E], [zero(diag(C₀).B)]))
    BinRescaledOp(C₀,Cbins,θname)
end

function get_Cℓ(f1::FlatS2, f2::FlatS2=f1; which=(:EE,:BB), kwargs...)
    Cℓ = (;[Symbol(x1*x2) => get_Cℓ(getindex(f1,Symbol(x1)),getindex(f2,Symbol(x2)); kwargs...) for (x1,x2) in split.(string.(ensure1d(which)),"")]...)
    which isa Symbol ? Cℓ[1] : Cℓ
end

function ud_grade(f::FlatS2{P}, args...; kwargs...) where {P}
    f′ = FieldTuple(map(f->ud_grade(f, args...; kwargs...), f.fs))
    B′ = (f′[1] isa FlatMap) ? (f isa FlatQU ? QUMap : EBMap) : (f isa FlatQU ? QUFourier : EBFourier)
    FieldTuple{B′}(f′)
end


text:
using Revise
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters
using Distributed
using SharedArrays
using JLD2
using RCall
using DataFrames
using Dates: today

using ShaleDrillingLikelihood: _ψ2, royalty, produce, num_i, value_function, EVobj, cost, drill
using Base.Iterators: flatten, product, OneTo

using Formatting: generate_formatter
num_nonzeros(x) = length(nonzeros(x))
intstr = generate_formatter("%'d")

# ------------------- number of simulations ----------------------

pargs = parse_commandline_counterfactuals()
print_parsed_args_counterfactuals(pargs)

DO_PAR         = !pargs["noPar"]
JLD2FILE       = pargs["jld2"]
DATE_STOP      = pargs["dateStop"]
TECH_YEAR_ZERO = pargs["techYearZero"]
RFILEDIR       = pargs["rFileDir"]

# --------------- create data ---------------

if "SLURM_JOBID" in keys(ENV)
    DATADIR = "/home/magerton/haynesville/intermediate_data"
    SLURM_JOBID = ENV["SLURM_JOBID"]
else
    DATADIR = "E:/projects/haynesville/intermediate_data"
    SLURM_JOBID = "999999"
end

println_time_flush("Loading results from $JLD2FILE")
file = jldopen(JLD2FILE, "r")
    DATAPATH = file["DATAPATH"]
    M        = file["M"]
    ddmnovf  = file["ddm_novf"]
    theta    = file["theta1"]
close(file)
println_time_flush("Done!")

REWARD = reward(ddmnovf)
ddm = DynamicDrillModel(ddmnovf, REWARD)

# --------------- create data ---------------

# load in data from disk
rdatapath = joinpath(DATADIR, DATAPATH)
println_time_flush("loading $rdatapath")

data_royalty    = DataRoyalty(REWARD, rdatapath)
data_produce    = DataProduce(REWARD, rdatapath)
data_drill_prim = DataDrillPrimitive(REWARD, rdatapath)
data_drill      = DataDrill(data_drill_prim, ddm)

evdims = size(EVobj(value_function(ddm)))
println_time_flush("EV dimension is $evdims, implying $(intstr(prod(evdims))) states")

# full dataset
dataset_full = DataSetofSets(data_drill, data_royalty, data_produce, CoefLinks(data_drill))

println_time_flush("Data created")

let thts = split_thetas(dataset_full, theta)
    println("theta_drill = $(round.(thts[1]; digits=3))")
    println("theta_roy   = $(round.(thts[2]; digits=3))")
    println("theta_pdxn  = $(round.(thts[3]; digits=3))")
end

# ------------------- parallel startup -----------------------

if DO_PAR
    pids = start_up_workers(ENV)
    @everywhere using ShaleDrillingLikelihood
    println_time_flush("Library loaded on workers")
end

# ------------------- simulations -----------------------

N = num_i(dataset_full)
posteriors = SharedPosterior(M, N; pids = pids)
sim = SimulationDraws(M, dataset_full)

@eval @everywhere set_g_BaseDataSetofSets($dataset_full)
@eval @everywhere set_g_SimulationDrawsMatrix($sim)
@eval @everywhere set_g_SharedPosterior($posteriors)
@eval @everywhere ShaleDrillingLikelihood.update!(get_g_BaseDataSetofSets(), $theta)
println_time_flush("starting simulations")

# map(i -> simloglik_posterior!(i, theta, sim, posteriors, dataset_full), OneTo(N))
pmap(i -> simloglik_posterior!(i, theta), CachingPool(pids), OneTo(N))
rmprocs(pids)

unitidx = [Int32(n) for (m,n) in product(OneTo(M), OneTo(N))]
psi2    = _ψ2(sim)
post_d  = sdata(drill(posteriors))
post_r  = sdata(royalty(posteriors))
post_p  = sdata(produce(posteriors))

posterior_df = DataFrame(
    uniti      = vec(unitidx),
    psi2       = vec(psi2),
    post_drill = vec(post_d),
    post_roy   = vec(post_r),
    post_pdxn  = vec(post_p)
)

# ------------------- save -----------------------

dpath_no_rdata = replace(DATAPATH, ".RData" => "")
filenm = "posterior-$(SLURM_JOBID)-$(dpath_no_rdata)-$(today()).RData"
FILEPATH = joinpath(RFILEDIR, filenm)

println("saving posterior to $FILEPATH")
flush(stdout)

@rput posterior_df FILEPATH

R"""
library(data.table)
setDT(posterior_df)
save(posterior_df, file=FILEPATH)
rm(posterior_df, FILEPATH)
gc()
"""

println("all done! :)")
flush(stdout)


text:
"""
`trivial_allocation!`: Trivially allocate each particle to have one offspring.

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
"""
function trivial_allocation!(E::TE, B::TB) where {TE<:Ensemble, TB<:Bins}
   @. E.o = 1;
   @. B.target = B.n;
   E, B
end

"""
`minimal_bin_allocation!`: Allocates a single particle to be spawned within each
nonempty bin

### Arguments
* `B` - bin data structure
"""
function minimal_bin_allocation!(B::TB) where {TB<:Bins}

   # increment the number of offspring by one for each nonempty bin
   @. B.target += (B.n > 0)
   B
end

"""
`targeted_bin_allocation!`: Targeted allocation of particles amongst bins using
a specified function, `G:(p, E, B, t) → [0,∞)` for bin `p`. Falls back to
uniform allocation amongst the non-empty bins in the event that this `G` fails
to normalize.

### Arguments
* `B` - bin data structure
* `E` - particle ensemble
* `G` - target function
* `t` - t-th seletion step
* `n_allocate` - number of particles to allocate
  ### Optional Arguments
* `allocation_resampler=systematic` - resampling scheme amongst bins
"""
function targeted_bin_allocation!(B::TB, E::TE, G::F, t::Int, n_allocate::Int; allocation_resampler = systematic) where {TE<:Ensemble,TB<:Bins,F<:Function}

   n_bins = length(B)

   # identify nonempty bins
   non_empty_bins = findall(n -> n > 0, B.n);
   Ñ = zeros(n_bins);
   ρ = zeros(n_bins);

   for p in non_empty_bins
      Ñ[p] = G(p, E, B, t);
   end

   # compute probabilities when this can be normalized
   ρ .= Ñ / sum(Ñ);
   B.target .+= allocation_resampler(n_allocate, ρ);

   B
end

"""
`optimal_bin_allocation!`: Optimally particles according to the bins,
using a value function to approximate mutation variance.

### Arguments
* `B` - bin data structure
* `E` - particle ensemble
* `v²` - v² variance function estimator
* `t` - t-th seletion step
* `n_allocate` - number of particles to allocate
### Optional Arguments
* `allocation_resampler=systematic` - resampling scheme amongst bins
"""
function optimal_bin_allocation!(B::TB, E::TE, v²::F, t::Int, n_allocate::Int; allocation_resampler = systematic) where {TE<:Ensemble,TB<:Bins,F<:Function}

   function G(p, E, B, t)
      particle_ids = findall(isequal(p), E.b)
      bin_target = sqrt(B.ν[p] * (E.ω[particle_ids] ⋅ v².(E.ξ[particle_ids], t)))
      return bin_target
   end

   targeted_bin_allocation!(B, E, G, t, n_allocate, allocation_resampler = allocation_resampler)

   B
end

"""
`uniform_bin_allocation!`: Uniformly allocate particles amongst bins.

### Arguments
* `B` - bin data structure
* `E` - particle ensemble
* `n_allocate` - number of particles to allocate
### Optional Arguments
* `allocation_resampler=systematic` - resampling scheme amongst bins
"""
function uniform_bin_allocation!(B::TB, E::TE, n_allocate::Int; allocation_resampler = systematic) where {TE<:Ensemble,TB<:Bins}
   G = (p, E, B, t) -> 1.0
   targeted_bin_allocation!(B, E, G, 0, n_allocate, allocation_resampler = allocation_resampler)
   B
end

"""
`within_bin_allocation!`: Once the number of offspring within each bin are set,
allocate them amongst the particles within the bin.  This assumes that the bin
allocations of the bins have completed.

### Arguments
* `E` - particle ensemble
* `B` - bin data structure
### Optional Arguments
* `within_bin_resampler=multinomial` - resampling scheme within bins
"""
function within_bin_allocation!(E::TE, B::TB; within_bin_resampler = multinomial) where {TE<:Ensemble,TB<:Bins}
   n_particles = length(E)
   non_empty_bins = findall(n -> n > 0, B.n)

   # compute number of offspring of each particle bin by bin
   for p in non_empty_bins
      # get particle indices for bin p
      particle_ids = findall(isequal(p), E.b)
      @inbounds E.o[particle_ids] .+= within_bin_resampler(B.target[p], E.ω[particle_ids] / B.ν[p])
   end
   E
end


text:
"""
A `struct` for individuals that keeps individual-specific variables.
"""
mutable struct Ind{B<:AbstractFloat, C<:AbstractArray, D<:AbstractArray, E<:AbstractArray} <: AbstractAgent
  id::Int  # the individual ID
  pos::Tuple{Int, Int}  # the individuals position
  species::Int  # the species ID the individual belongs to
  biotic_phenotype::Vector{B}
  abiotic_phenotype::Vector{B}
  epistasisMat::C  # epistasis matrix
  pleiotropyMat::D  # pleiotropy matrix
  q::E  # expression array
  age::Int
  sex::Bool
  interaction_history::MArray{S, Int} where S # records the last interaction with all species
  energy::B  # determines whether need to feed (energy=0) or not.
  W::B  # survival probability
  isalive::Bool
end

struct Params{F<:AbstractFloat, I<:Int, N<:AbstractString}
  ngenes::Vector{I}
  nphenotypes::Vector{I}
  growthrates::Vector{F}
  selectionCoeffs::Vector{F}
  ploidy::Vector{I}
  optvals::Vector{Vector{Vector{Matrix{F}}}}
  optinds::Vector{Vector{I}}
  mutProbs::Vector{Vector{DiscreteNonParametric{Bool, Float64, Vector{Bool}, Vector{Float64}}}}
  mutMagnitudes::Vector{Vector{UnivariateDistribution{S} where S<:ValueSupport}}
  N::Vector{Vector{I}}
  E::Vector{Normal{F}}
  generations::I
  nspecies::I
  new_N::Vector{Vector{I}}
  migration_traits::Vector{I}
  vision_radius::Vector{I}
  check_fraction::Vector{F}
  migration_thresholds::Vector{F}
  step::MVector{1, Int64}
  nodes::Matrix{Tuple{I, I}}
  biotic_phenotypes::Vector{Vector{I}}
  abiotic_phenotypes::Vector{Vector{I}}
  max_ages::Vector{I}
  names::Dict{I, N}
  food_sources::Matrix{F}
  interactions::Matrix{F}
  resources::Matrix{I}
  resources_org::Matrix{I}
  recombination::Vector{Poisson{F}}
  initial_energy::Vector{F}
  bottlenecks::Vector{BitMatrix}
  repro_start::Vector{Int}
  repro_end::Vector{Int}
  seed::Int
end

const variance = 1.0

"""
    model_initiation(dd)

Innitializes the model.
"""
function model_initiation(dd)

  if !isnothing(dd[:model]["seed"])
    Random.seed!(dd[:model]["seed"])
  else
    rnd = Int(rand(1:1e6))
    dd[:model]["seed"] = rnd
    Random.seed!(rnd)
  end

  properties, species_arrays = create_properties(dd)
  epistasisMat, pleiotropyMat, expressionArrays = species_arrays

  space = dd[:model]["space"]

  if isnothing(space)
    fspace = GridSpace((1, 1))
  elseif typeof(space) <: Tuple
    fspace = GridSpace(space, periodic=dd[:model]["periodic"], metric=dd[:model]["metric"])
  end

  indtype = EvoDynamics.Ind{typeof(0.1), eltype(epistasisMat), eltype(pleiotropyMat), eltype(expressionArrays)}

  model = ABM(indtype, fspace, properties=properties, scheduler=Schedulers.randomly)

  # create and add agents
  for sp in 1:model.nspecies
    for (pos, n) in enumerate(model.N[sp])
      for ind in 1:n
        abiotic_ph = get_abiotic_phenotype(sp, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], model)
        biotic_ph = get_biotic_phenotype(sp, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], model)
        W = abiotic_fitness(abiotic_ph, sp, pos, model)
        W = adjust_fitness(W, sp, model)
        sex = false
        if model.ploidy[sp] == 2
          sex = rand((true, false))
        end
        interaction_history = MVector{model.nspecies, Int}(fill(0, model.nspecies))
        initial_energy = model.initial_energy[sp]
        add_agent!(model.nodes[pos], model, sp, biotic_ph, abiotic_ph, epistasisMat[sp], pleiotropyMat[sp], expressionArrays[sp], 0, sex, interaction_history, initial_energy, W, true)
      end
    end
  end

  return model
end

nnodes(x::GridSpace) = prod(size(x))
nnodes(x::ABM) = nnodes(x.space)

function create_properties(dd)
  nspecies = length(dd[:species])

  Ed = [Normal(0.0, dd[:species][i]["environmental noise"]) for i in 1:nspecies]
  Mdists = [[DiscreteNonParametric([true, false], [i, 1-i]) for i in dd[:species][arr]["mutation probabilities"]] for arr in 1:nspecies]  # μ (probability of change)
  Ddists = [[Normal(0, dd[:species][ar]["mutation magnitudes"][1]), DiscreteNonParametric([true, false], [dd[:species][ar]["mutation magnitudes"][2], 1-dd[:species][ar]["mutation magnitudes"][2]]), Normal(0, dd[:species][ar]["mutation magnitudes"][3])] for ar in 1:nspecies]  # amount of change in case of mutation

  # make single-element arrays 2D so that linAlg functions will work
  newA = Array{Array{Float64}}(undef, nspecies)
  newQ = Array{Array{Float64}}(undef, nspecies)
  for i in 1:nspecies
    if length(dd[:species][i]["epistasis matrix"]) == 1
      newA[i] = reshape(dd[:species][i]["epistasis matrix"], 1, 1)
      newQ[i] = reshape(dd[:species][i]["expression array"], 1, 1)
    else
      newA[i] = dd[:species][i]["epistasis matrix"]
      newQ[i] = dd[:species][i]["expression array"]
    end
  end

  epistasisMatS = [MArray{Tuple{size(newA[i])...}}(newA[i]) for i in eachindex(newA)]
  pleiotropyMatS = [MArray{Tuple{size(dd[:species][i]["pleiotropy matrix"])...}}(dd[:species][i]["pleiotropy matrix"]) for i in 1:nspecies]
  expressionArraysS = [MArray{Tuple{size(newQ[i])...}}(newQ[i]) for i in eachindex(newQ)]

  ngenes = [dd[:species][i]["number of genes"] for i in 1:nspecies]
  nphenotypes = [dd[:species][i]["number of phenotypes"] for i in 1:nspecies]
  growthrates = [dd[:species][i]["growth rate"] for i in 1:nspecies]
  selectionCoeffs = [dd[:species][i]["selection coefficient"] for i in 1:nspecies]
  ploidy = [dd[:species][i]["ploidy"] for i in 1:nspecies]
  optvals = [dd[:species][i]["optimal phenotype values"] for i in 1:nspecies]
  optinds = [dd[:species][i]["optimal phenotypes"] for i in 1:nspecies]
  Ns = [dd[:species][i]["N"] for i in 1:nspecies]
  migration_traits = [dd[:species][i]["migration phenotype"] for i in 1:nspecies]
  vision_radius = [dd[:species][i]["vision radius"] for i in 1:nspecies]
  check_fraction = [dd[:species][i]["check fraction"] for i in 1:nspecies]
  migration_thresholds = [dd[:species][i]["migration threshold"] for i in 1:nspecies]
  generations = dd[:model]["generations"]
  step = MVector{1, Int}(undef)
  step[1] = 0
  nnodes = Matrix{Tuple{Int, Int}}(undef, dd[:model]["space"]...)
  for col in 1:size(nnodes, 2)
    for row in 1:size(nnodes, 1)
      nnodes[row, col] = (row, col)
    end
  end
  biotic_phenotyps = [dd[:species][i]["biotic phenotypes"] for i in 1:nspecies]
  abiotic_phenotyps = [dd[:species][i]["abiotic phenotypes"] for i in 1:nspecies]
  max_ages = [dd[:species][i]["age"] for i in 1:nspecies]
  names = Dict(i => dd[:species][i]["name"] for i in 1:nspecies)
  recombination = [Poisson(dd[:species][i]["recombination"]) for i in 1:nspecies]
  initial_energy = [AbstractFloat(dd[:species][i]["initial energy"]) for i in 1:nspecies]
  bottlenecks = [dd[:species][i]["bottleneck times"] for i in 1:nspecies]
  repro_start = [dd[:species][i]["reproduction start age"] for i in 1:nspecies]
  repro_end = [dd[:species][i]["reproduction end age"] for i in 1:nspecies]

  properties = Params(ngenes, nphenotypes, growthrates, selectionCoeffs, ploidy, optvals, optinds, Mdists, Ddists, Ns, Ed, generations, nspecies, Ns, migration_traits, vision_radius, check_fraction, migration_thresholds, step, nnodes, biotic_phenotyps, abiotic_phenotyps, max_ages, names, dd[:model]["food sources"], dd[:model]["interactions"], dd[:model]["resources"], deepcopy(dd[:model]["resources"]), recombination, initial_energy, bottlenecks, repro_start, repro_end, dd[:model]["seed"])

  return properties, (epistasisMatS, pleiotropyMatS, expressionArraysS)
end

function return_opt_phenotype(species::Int, generation::Int, site::Int, model::ABM)
  nabiotic = length(model.abiotic_phenotypes[species])
  output = Array{typeof(0.1)}(undef, nabiotic)
  for trait in 1:nabiotic
    output[trait] = model.optvals[species][model.optinds[species][generation+1]][trait][site]
  end
  return output
end

function return_opt_phenotype(species::Int, generation::Int, site::Tuple{Int, Int}, model::ABM)
  nabiotic = length(model.abiotic_phenotypes[species])
  output = Array{typeof(0.1)}(undef, nabiotic)
  for trait in 1:nabiotic
    output[trait] = model.optvals[species][model.optinds[species][generation+1]][trait][site[1],site[2]]
  end
  return output
end

function model_step!(model::ABM)
  model.step[1] += 1
  model.resources .= model.resources_org
end

function agent_step!(agent::Ind, model::ABM)
  # update age
  agent.age += 1
  # abiotic survive
  abiotic_survive!(agent, model)  # the agent first survives then feeds, reproduces, and interacts.
  if !agent.isalive
    return
  end
  # use food
  burn_energy!(agent)
  # consume basic energy if agent can
  consume_food!(agent, model)
  # interact with other species
  interact!(agent, model)
  # Kill the agent if it doesn't have energy
  if agent.isalive && agent.energy < 0
    remove_agent!(agent, model)
    return
  end
  # survive
  survive!(agent, model)
  # migrate
  migrate!(agent, model)
  # reproduction for the haploid
  reproduce!(agent, model)
  # bottleneck
  if agent.isalive && model.bottlenecks[agent.species][coord2vertex(agent.pos, model), model.step[1]]
    if EvoDynamics.bottleneck(agent, model)
      remove_agent!(agent, model)
    end
  end
end

function bottleneck(agent, model)
  # @info "Package bottleneck function"
  return false
end

function remove_agent!(agent, model)
  agent.isalive = false
  kill_agent!(agent, model)
end

"""
    survive!(agent::Ind, model::ABM)

Kills the agent if it has no energy, or is too old.
"""
function survive!(agent::Ind, model::ABM)
  if !agent.isalive
    return
  elseif agent.energy < 0
    remove_agent!(agent, model)
  elseif agent.age ≥ model.max_ages[agent.species]
    remove_agent!(agent, model)
  # elseif rand() > agent.W
  #   remove_agent!(agent, model)
  end
end

"""
Kills the agent by chance given its fitness `W`
"""
function abiotic_survive!(agent::Ind, model::ABM)
  if agent.isalive && rand() > agent.W
    remove_agent!(agent, model)
  end
end

function adjust_fitness!(agent::Ind, model::ABM)
  W = agent.W < 0 ? 0.0 : agent.W
  newW = 1.0 - ( (1.0 - W) * model.selectionCoeffs[agent.species])
  agent.W = newW
end

function adjust_fitness(W, species, model::ABM)
  W2 = W < 0 ? 0.0 : W
  newW = 1.0 - ( (1.0 - W2) * model.selectionCoeffs[species])
  return newW
end

function burn_energy!(agent)
  agent.energy -= 1
end

function consume_food!(agent::Ind, model::ABM)
  environmental_consumption = model.food_sources[agent.species, agent.species]
  if environmental_consumption > 0 && model.resources[agent.pos...] >= environmental_consumption
    model.resources[agent.pos...] -= environmental_consumption
    agent.energy += environmental_consumption
  end
end

"Mutate an agent."
function mutate!(agent::Ind, model::ABM)
  mutated = false
  # mutate gene expression
  if rand(model.mutProbs[agent.species][1])
    agent.q .+= rand(model.mutMagnitudes[agent.species][1], model.ngenes[agent.species] * model.ploidy[agent.species])
    mutated = true
  end
  # mutate pleiotropy matrix
  if rand(model.mutProbs[agent.species][2])
    randnumbers = rand(model.mutMagnitudes[agent.species][2], size(agent.pleiotropyMat))
    agent.pleiotropyMat[randnumbers] .= .!agent.pleiotropyMat[randnumbers]
    mutated = true
  end
  # mutate epistasis matrix
  if rand(model.mutProbs[agent.species][3])
    agent.epistasisMat .+= rand(model.mutMagnitudes[agent.species][3], size(agent.epistasisMat))
    mutated = true
  end
  # update biotic and abiotic phenotypes and W
  if mutated
    update_fitness!(agent, model)
  end
end

coord2vertex(c, model) = c[1] + (size(model.space)[1] * (c[2]-1))

function update_fitness!(agent::Ind, model::ABM)
  abiotic_phenotype = get_abiotic_phenotype(agent, model)
  biotic_phenotype = get_biotic_phenotype(agent, model)
  W = abiotic_fitness(abiotic_phenotype, agent.species, agent.pos, model)
  W = adjust_fitness(W, agent.species, model)
  agent.biotic_phenotype .= biotic_phenotype
  agent.abiotic_phenotype .= abiotic_phenotype
  agent.W = W
end


text:
@testset "Persistent Image" begin

    intrs = diagram(0.0=>Inf, 0.0=>1.0, 1.0=>1.0, 1.0=>2.0)

    pimg = PersistentImage(intrs...)
    @test length(pimg.intervals) == length(intrs)
    @test pimg(0,0) ≈ 0.0627 atol=1e-3
    @test length(vec(pimg, 5, 6)) == 30

end