text:
module PkgButlerEngine

import Mustache
import Pkg

function configure_pkg(path::AbstractString; channel=:auto, template=:auto)
    channel in (:auto, :stable, :dev) || error("Invalid value for channel.")
    template in (:auto, :default, :bach) || error("Invalid value for template.")

    path_for_butler_workflows_folder = joinpath(path, ".github", "workflows")

    path_for_main_butler_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-butler-workflow.yml")
    path_for_main_butler_dev_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-butler-dev-workflow.yml")

    mkpath(path_for_butler_workflows_folder)

    if channel == :auto
        channel = isfile(path_for_main_butler_dev_workflow) ? :dev : :stable
    end

    if channel == :stable
        if isfile(path_for_main_butler_dev_workflow)
            rm(path_for_main_butler_dev_workflow, force=true)
        end

        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-butler-workflow.yml"), path_for_main_butler_workflow, force=true)
    elseif channel == :dev
        if isfile(path_for_main_butler_workflow)
            rm(path_for_main_butler_workflow, force=true)
        end

        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-butler-dev-workflow.yml"), path_for_main_butler_dev_workflow, force=true)
    end

    path_for_config_file = joinpath(path, ".jlpkgbutler.toml")

    if template !== :auto
        if isfile(path_for_config_file)
            config_content = Pkg.TOML.parsefile(path_for_config_file)

            if haskey(config_content, "template")
                if template !== :auto
                    config_content["template"] = string(template)

                    open(path_for_config_file, "w") do f
                        Pkg.TOML.print(f, config_content)
                    end
                end
            elseif template !== :default
                config_content["template"] = string(template)

                open(path_for_config_file, "w") do f
                    Pkg.TOML.print(f, config_content)
                end
            end
        elseif template !== :default
            open(path_for_config_file, "w") do f
                Pkg.TOML.print(f, Dict{String,Any}("template" => string(template)))
            end

        end
    end
end

function cp_with_mustache(src, dest, vals)
    content = read(src, String)

    open(dest, "w") do file
        Mustache.render(file, content, vals, tags=("\$[[", "]]"))
    end
end

function ensure_project_has_julia_compat(path)
    proj_file = isfile(joinpath(path, "JuliaProject.toml")) ? joinpath(path, "JuliaProject.toml") : joinpath(path, "Project.toml")

    pkg_toml_content = Pkg.TOML.parsefile(proj_file)

    if !haskey(pkg_toml_content, "compat")
        pkg_toml_content["compat"] = Dict{String,String}()
    end

    if !haskey(pkg_toml_content["compat"], "julia")
        pkg_toml_content["compat"]["julia"] = "1"

        open(proj_file, "w") do f
            Pkg.TOML.print(f, pkg_toml_content)
        end
    end
end

function ensure_project_uses_new_enough_documenter(path)
    doc_proj_file = isfile(joinpath(path, "docs", "JuliaProject.toml")) ? joinpath(path, "docs", "JuliaProject.toml") : joinpath(path, "docs", "Project.toml")

    if isfile(doc_proj_file)
        pkg_toml_content = Pkg.TOML.parsefile(doc_proj_file)

        if haskey(pkg_toml_content, "compat") && haskey(pkg_toml_content["compat"], "Documenter")
            documenter_compat_bound = pkg_toml_content["compat"]["Documenter"]
            version_bound = Pkg.Types.semver_spec(documenter_compat_bound)

            # This is the list of versions that don't work
            invalid_versions = Pkg.Types.semver_spec("0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.20,0.21,0.22,0.23")

            if !isempty(intersect(invalid_versions, version_bound))
                # TODO This is a bit crude. Ideally we would try to see whether a version >= 0.24 is listed in the compat section, and if so,
                # only remove the offending versions.
                pkg_toml_content["compat"]["Documenter"] = "~0.24"

                open(doc_proj_file, "w") do f
                    Pkg.TOML.print(f, pkg_toml_content)
                end

            end
        end

    end
end

function construct_version_matrix(path)
    proj_file = isfile(joinpath(path, "JuliaProject.toml")) ? joinpath(path, "JuliaProject.toml") : joinpath(path, "Project.toml")

    pkg_toml_content = Pkg.TOML.parsefile(proj_file)

    julia_compat_bound = pkg_toml_content["compat"]["julia"]

    version_spec = Pkg.Types.semver_spec(julia_compat_bound)

    versions = [v"1.0.5" => "'1.0'", v"1.1.1" => "'1.1'", v"1.2.0" => "'1.2'", v"1.3.1" => "'1.3'", v"1.4.2" => "'1.4'", v"1.5.3" => "'1.5'", v"1.6.5" => "'1.6'", v"1.7.2" => "'1.7'"]

    compat_versions = filter(i -> i[1] in version_spec, versions)

    return join(map(i -> i[2], compat_versions), ", ")
end

function construct_matrix_exclude_list(path)
    path_for_config_file = joinpath(path, ".jlpkgbutler.toml")

    if isfile(path_for_config_file)
        config_content = Pkg.TOML.parsefile(path_for_config_file)

        if haskey(config_content, "strategy-matrix-exclude")
            option_value = config_content["strategy-matrix-exclude"]

            line_ending = Sys.iswindows() ? "\r\n" : "\n"

            exclude_configs = split(option_value, ";", keepempty=false)
            exclude_configs = strip.(exclude_configs)

            ret = ""

            for ec in exclude_configs

                lines = split(ec, ",", keepempty=false)
                lines = strip.(lines)

                ret *= line_ending * " "^10 * "- " * lines[1] * ( length(lines) > 1 ? line_ending * join(string.(" "^12, lines[2:end]), line_ending) : "" )
            end

            return ret
        end
    end

    return ""
end

function configure_tagbot!(path, view_vals)
    path_for_config_file = joinpath(path, ".jlpkgbutler.toml")

    if isfile(path_for_config_file)
        config_content = Pkg.TOML.parsefile(path_for_config_file)

        if haskey(config_content, "custom-registry")
            option_value = config_content["custom-registry"]

            view_vals["include_custom_registry"] = "true"
            view_vals["JL_CUSTOM_REGISTRY"] = option_value
        end
    end
end

function add_compathelper(path)
    path_for_butler_workflows_folder = joinpath(path, ".github", "workflows")
    path_for_compathelper_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-compathelper-workflow.yml")

    cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-compathelper-workflow.yml"), path_for_compathelper_workflow, force=true)
end

function detect_template(path)
    path_for_config_file = joinpath(path, ".jlpkgbutler.toml")

    if isfile(path_for_config_file)
        config_content = Pkg.TOML.parsefile(path_for_config_file)

        if haskey(config_content, "template")
            return lowercase(config_content["template"])
        end
    end

    return "default"
end

function update_pkg_bach(path)
    if isfile(joinpath(path, ".travis.yml"))
        rm(joinpath(path, ".travis.yml"), force=true)
    end

    if isfile(joinpath(path, "appveyor.yml"))
        rm(joinpath(path, "appveyor.yml"), force=true)
    end

    if isfile(joinpath(path, ".appveyor.yml"))
        rm(joinpath(path, ".appveyor.yml"), force=true)
    end
end

function update_pkg(path::AbstractString)
    configure_pkg(path)

    template = detect_template(path)

    path_for_butler_workflows_folder = joinpath(path, ".github", "workflows")

    path_for_ci_master_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-ci-master-workflow.yml")
    path_for_ci_pr_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-ci-pr-workflow.yml")
    path_for_docdeploy_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-docdeploy-workflow.yml")
    path_for_codeformat_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-codeformat-pr-workflow.yml")
    path_for_tagbot_workflow = joinpath(path_for_butler_workflows_folder, "jlpkgbutler-tagbot-workflow.yml")

    path_for_docs_make_file = joinpath(path, "docs", "make.jl")

    mkpath(path_for_butler_workflows_folder)

    ensure_project_has_julia_compat(path)

    view_vals = Dict{String,Any}()
    view_vals["JL_VERSION_MATRIX"] = construct_version_matrix(path)
    # if template == "bach"
    #     view_vals["include_codeformat_lint"] = "true"
    # end
    view_vals["ADDITIONAL_MATRIX_EXCLUDES"] = construct_matrix_exclude_list(path)
    configure_tagbot!(path, view_vals)

    cp_with_mustache(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-ci-master-workflow.yml"), path_for_ci_master_workflow, view_vals)
    cp_with_mustache(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-ci-pr-workflow.yml"), path_for_ci_pr_workflow, view_vals)
    cp_with_mustache(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-tagbot-workflow.yml"), path_for_tagbot_workflow, view_vals)

    if isfile(path_for_docs_make_file)
        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-docdeploy-workflow.yml"), path_for_docdeploy_workflow, force=true)
    else isfile(path_for_docdeploy_workflow)
        rm(path_for_docdeploy_workflow, force=true)
    end

    ensure_project_uses_new_enough_documenter(path)

    add_compathelper(path)

    isfile(path_for_codeformat_workflow) && rm(path_for_codeformat_workflow, force=true)

    if template == "bach"
        update_pkg_bach(path)

        cp(joinpath(@__DIR__, "..", "templates", "jlpkgbutler-codeformat-pr-workflow.yml"), path_for_codeformat_workflow, force=true)
    end
end

end # module


text:
using Modia
using Base.Test

@testset "RunSimulations" begin

Modia.ModiaLogging.resetTestStatus()

include("models/TestVariableTypes.jl")
include("models/TestEquations.jl")

include("models/TestElectrical.jl")
include("models/TestFilter.jl")
include("models/TestArrayOfComponents.jl")
include("models/TestSpatialDiscretization.jl")

#include("models/TestCoupledInertias.jl")
#include("models/TestPendulum.jl")
#include("models/TestStateSelection.jl")
#include("models/TestRetranslationOmega.jl")

Modia.ModiaLogging.printTestStatus()

end # testset
nothing

ext: jl


text:
import Base.@kwdef

abstract type AbstractSpline end

@kwdef struct Spline
    points::Vector{Point}
    duration::Real = 10000.
    sharpness::Real = 0.85
    centers = []
    controls = [[]]
    stepLength::Real = 60
    length::Real = length(points)
    delay::Real = 0
    steps = []

    function Spline(; points=points::Vector{Point}, duration=duration::Real=10000., sharpness=sharpness::Real=0.85,
        centers=centers=[],
        controls=controls=[[]],
        stepLength=stepLength::Real=60,
        length=length::Real = length(points),
        delay=delay::Real= 0,
        steps=steps= [])

        for point in points
            point.coordinates[3] = length(point) > 2 ? point.coordinates[3] : 0
        end

        for i in eachindex(points)
            p1 = points[i]
            p2 = points[i + 1]

            push!(centers,
                Point([(p1.coordinates[1] + p2.coordinates[1]) / 2,
                (p1.coordinates[2] + p2.coordinates[2]) / 2,
                (p1.coordinates[3] + p2.coordinates[3]) / 2]))
        end
        push!(controls, [points[1], points[1]])

        for j in eachindex(centers)
            c1 = centers[j]
            c2 = centers[j + 1]
            dx = points[j + 1].coordinates[1] - (c1.coordinates[1] + c2.coordinates[1]) / 2
            dy = points[j + 1].coordinates[2] - (c1.coordinates[2] + c2.coordinates[2]) / 2
            dz = points[j + 1].coordinates[3] - (c1.coordinates[3] + c2.coordinates[3]) / 2

            push!(controls, [
                Point([(1. - sharpness) * points[j + 1].coordinates[1] + sharpness * c1.coordinates[1] + dx,
                (1. - sharpness) * points[j + 1].coordinates[2] + sharpness * c1.coordinates[2] + dy,
                (1. - sharpness) * points[j + 1].coordinates[3] + sharpness * c1.coordinates[3] + dz]),
                Point([(1. - sharpness) * points[j + 1].coordinates[1] + sharpness * c2.coordinates[1] + dx,
                (1. - sharpness) * points[j + 1].coordinates[2] + sharpness * c2.coordinates[2] + dy,
                (1. - sharpness) * points[j + 1].coordinates[3] + sharpness * c2.coordinates[3] + dz])
                ])
        end

        push!(controls, [points[length - 1], points[length - 1]])
        steps = cacheSteps(stepLength)

        return new(points, duration, sharpness, centers, controls, stepLength, length, delay, steps)
    end

end # struct

function cacheSteps(minDist::Real)
    steps = Real[]
    lastStep = pos(0)
    push!(steps, 0)

    for i in 1:10:duration
        step = pos(i)
        dist = sqrt(
            (step.coordinates[1] - lastStep.coordinates[1])^2 +
            (step.coordinates[2] - lastStep.coordinates[2])^2 +
            (step.coordinates[3] - lastStep.coordinates[3])^2)

        if dist > minDist
            push!(steps, i)
            lastStep = step
        end
    end

    return steps

end

function pos(time::Real)
    t = time - delay
    t < 0 && (t = 0)
    t > duration && (t = duration - 1)

    t2 = t / duration
    t2 >= 1 && return points[length - 1]

    n = floor((length(points) - 1) * t2)
    t1 = (length - 1) * t2 - n

    return bezier(t1, points[n], controls[n][1], controls[n + 1][0], points[n + 1])
end

function bezier(; t::Real, p1::Point, c1::Point, c2::Point, p2::Point)
    b = [t^3, (3 * t^2 * (1-t)), (3 * t * (1 - t) * (1 - t)), ((1 - t) * (1 - t) * (1 - t))]

    return Point([
        p2.coordinates[1] * b[1] + c2.coordinates[1] * b[2] + c1.coordinates[1] * b[3] + p1.coordinates[1] * b[4],
        p2.coordinates[2] * b[1] + c2.coordinates[2] * b[2] + c1.coordinates[2] * b[3] + p1.coordinates[2] * b[4],
        p2.coordinates[3] * b[1] + c2.coordinates[3] * b[2] + c1.coordinates[3] * b[3] + p1.coordinates[3] * b[4]])
end


text:
function get_bb(obj::FieldDSSTZonalContext)
    return jcall(obj, "getBB", RealFieldElement, ())
end

function get_bo_a(obj::FieldDSSTZonalContext)
    return jcall(obj, "getBoA", RealFieldElement, ())
end

function get_bo_a_bpo(obj::FieldDSSTZonalContext)
    return jcall(obj, "getBoABpo", RealFieldElement, ())
end

function get_cxo2n2a2(obj::FieldDSSTZonalContext)
    return jcall(obj, "getCXO2N2A2", RealFieldElement, ())
end

function get_hk(obj::FieldDSSTZonalContext)
    return jcall(obj, "getHK", RealFieldElement, ())
end

function get_k2mh2(obj::FieldDSSTZonalContext)
    return jcall(obj, "getK2MH2", RealFieldElement, ())
end

function get_k2mh2o2(obj::FieldDSSTZonalContext)
    return jcall(obj, "getK2MH2O2", RealFieldElement, ())
end

function get_m2ao_a(obj::FieldDSSTZonalContext)
    return jcall(obj, "getM2aoA", RealFieldElement, ())
end

function get_m_co2ab(obj::FieldDSSTZonalContext)
    return jcall(obj, "getMCo2AB", RealFieldElement, ())
end

function get_mean_motion(obj::FieldDSSTZonalContext)
    return jcall(obj, "getMeanMotion", RealFieldElement, ())
end

function get_muoa(obj::FieldDSSTZonalContext)
    return jcall(obj, "getMuoa", RealFieldElement, ())
end

function get_oo_ab(obj::FieldDSSTZonalContext)
    return jcall(obj, "getOoAB", RealFieldElement, ())
end

function get_oon2a2(obj::FieldDSSTZonalContext)
    return jcall(obj, "getOON2A2", RealFieldElement, ())
end

function get_roa(obj::FieldDSSTZonalContext)
    return jcall(obj, "getRoa", RealFieldElement, ())
end

function get_x(obj::FieldDSSTZonalContext)
    return jcall(obj, "getX", RealFieldElement, ())
end

function get_x2on2a2xp1(obj::FieldDSSTZonalContext)
    return jcall(obj, "getX2ON2A2XP1", RealFieldElement, ())
end

function get_x3on2a(obj::FieldDSSTZonalContext)
    return jcall(obj, "getX3ON2A", RealFieldElement, ())
end

function get_xon2a2(obj::FieldDSSTZonalContext)
    return jcall(obj, "getXON2A2", RealFieldElement, ())
end

function get_xx(obj::FieldDSSTZonalContext)
    return jcall(obj, "getXX", RealFieldElement, ())
end

function get_xxx(obj::FieldDSSTZonalContext)
    return jcall(obj, "getXXX", RealFieldElement, ())
end


text:
#!/usr/bin/env julia

# publication11

using Base.Test
using MathProgBase
using Logging
using Drs

A = Float64[1 1; 1 2; 3 1]
b = Float64[3, 5, 6]
c = Float64[-1, -2]

s = linprog(c, A, '<', b, -Inf, Inf, DrsMathProgSolver())

@test s.status == :Optimal
@test s.objval == -5
@test s.sol == [0, 5/2]

max_stars_repo_name: abogushov/lplib

max_issues_repo_name: abogushov/lplib

max_forks_repo_name: abogushov/lplib


text:
using MultiAgentReinforcementLearning
using Test

@testset "MultiAgentReinforcementLearning.jl" begin
    # Write your own tests here.
end


text:
"""
    $(SIGNATURES)

Creates an axis array with given range and units, as well as spin polarization if
`polarized` is `true`. The spin dimension is the last. The arguments pairs indicate the name
and ranges (or categories) of the axes.
"""
function zeros_axisarray(T::Type, polarized::Bool, axes::Vararg{Tuple{Symbol, Any}};
                         kwargs...)
    if polarized
        zeros_axisarray(T, axes..., kwargs..., (:spin, [:↑, :↓ ]))
    else
        zeros_axisarray(T, axes...; kwargs...)
    end
end
function zeros_axisarray(T::Type, axes::Vararg{Tuple{Symbol, Any}}; kwargs...)
    @argcheck all(typeof(axis[1]) <: Symbol for axis in axes)
    underlying = zeros(T, ((length(axis[2]) for axis in chain(axes, kwargs))...))
    AxisArray(underlying, (Axis{n}(r) for (n, r) in chain(axes, kwargs))...)
end

"""
    $(SIGNATURES)

Creates an array with named axis and with correct units for an electronic density. The first
argument can be either an actual quantity with the correct dimensions, or a unitless subtype
of number. In the latter case, the units default to Hartree atomic units. In the former
case, one should take care to give a fully qualified quantity, including actual units (at
least until triangular dispatch makes its way into Julia).
"""
function density_array{Q <: Number}(T::Type{Q}, spin::Bool, args...; kwargs...)
    @argcheck dimension(Q) == NoDims
    density_array(HartreeUnits.ρ{Q}, spin, args...; kwargs...)
end
function density_array{Q <: DFTUnits.Ρ}(T::Type{Q}, spin::Bool, args...; kwargs...)
    zeros_axisarray(T, spin, args...; kwargs...)
end

@traitfn function similar_potential_from_rho(U::Type, ρ::::(!HasSpinDim))
    similar(ρ, U)
end
@traitfn function similar_potential_from_rho(U::Type, ρ::::HasSpinDim)
    similar(view(ρ, Axis{:spin}(:↑)), U)
end

"""
    $(SIGNATURES)

Creates an axis arrays for the potential, guessing the exact type from the type of ρ.
"""
function similar_potential{Q <: DFTUnits.Ρ}(ρ::AxisArray{Q})
    @argcheck :radius ∈ axisnames(ρ)
    const i = findfirst(axisnames(ρ), :radius)
    const axis = axisvalues(ρ)[i]
    T = typeof(one(eltype(ρ)) / one(eltype(axis)))
    U = HartreeUnits.ϵ{T}
    fill!(similar_potential_from_rho(U, ρ), zero(U))
end


text:
"""
    mutable struct SourceTreeData{N,T,V}

Struct used to store data associated with a node of the `SourceTree`. The
parametric type `N` represents the ambient dimension, `T` represents the
primitive type used to represent the points (e.g. `Float32` or `Float64` for
single/double precision), and `V` is the type of the density that will be
interpolated.

# Fields:
- `msh`: uniform cartesian mesh of the interpolation space,
    indexable by a `CartesianIndex`
- `conedatadict`: dictionary mapping a cone tag `I::CartesianIndex` to an
  `Array` storing its interpolation values/coefficient.
- `farlist`: nodes on the target tree that can be evaluated through
    interpolation
- `nearlist`: nodes on the target tree that must be evaluated by a direct sum.
"""
mutable struct SourceTreeData{N,T}
    msh::UniformCartesianMesh{N,T}
    conedatadict::Dict{CartesianIndex{N},UnitRange{Int64}}
    farlist::Vector{TargetTree{N,T}}
    nearlist::Vector{TargetTree{N,T}}
end

function SourceTreeData{N,T}() where {N,T}
    domain       = HyperRectangle{N,T}(ntuple(i->0,N),ntuple(i->0,N))
    msh          = UniformCartesianMesh(domain, ntuple(i->1,N))
    conedict     = Dict{CartesianIndex{N},UnitRange{Int64}}()
    farlist      = Vector{TargetTree{N,T}}()
    nearlist     = Vector{TargetTree{N,T}}()
    return SourceTreeData{N,T}(msh, conedict, farlist, nearlist)
end

"""
    const SourceTree{N,T,V}

Type alias for a `ClusterTree` of points (represented as `SVector{N,T}`) storing
data of type [`SourceTreeData{N,T,V}`](@ref). See the documentation of
[`ClusterTree`](@ref) for more details.
"""
const SourceTree{N,T} = ClusterTree{Vector{SVector{N,T}},HyperRectangle{N,T},SourceTreeData{N,T}}

# getters for SourceTree
msh(t::SourceTree)                            = t.data.msh
conedatadict(t::SourceTree)                   = t.data.conedatadict
active_cone_idxs(t::SourceTree)               = keys(conedatadict(t))
farlist(t::SourceTree)                        = t.data.farlist
nearlist(t::SourceTree)                       = t.data.nearlist
center(t::SourceTree)            = t |> container |> center
cone_domains(t::SourceTree)                   = msh(t) |> ElementIterator # iterator of HyperRectangles
cone_domain(t::SourceTree, I::CartesianIndex) = cone_domains(t)[I]
conedata(t::SourceTree,I::CartesianIndex)    = conedatadict(t)[I]

# setters
_addtofarlist!(s::SourceTree, t::TargetTree)  = push!(farlist(s), t)
_addtonearlist!(s::SourceTree, t::TargetTree) = push!(nearlist(s), t)

"""
    _addnewcone!(node::SourceTree,I::CartesianIndex,p)

If `node` does not contain cone `I`, add it to the `conedict`.
"""
function _addnewcone!(node::SourceTree, idxcone, ::Val{P}) where {P}
    cd = conedatadict(node)
    if !haskey(cd, idxcone)
        push!(cd, idxcone => -1:-1)
    end
    return node
end

"""
    _init_msh!(node,domain,ds)

Create a uniform cartesian mesh of the interpolation domain of `node` with
element spacing `ds`.
"""
function _init_msh!(node::SourceTree, domain, ds)
    node.data.msh = UniformCartesianMesh(domain; step = ds)
    return node
end

"""
    interp2cart(s,xc,h)
    interp2cart(s,source::SourceTree)

Map interpolation coordinate `s` into the corresponding cartesian coordinates.

Interpolation coordinates are similar to `polar` (in 2d) or `spherical` (in 3d)
coordinates centered at the `xc`, except that the
radial coordinate `r` is replaced by `s = h/r`.

If passed a `SourceTree` as second argument, call `interp2cart(s,xc,h)` with
`xc` the center of the `SourceTree` and `h` its radius.
"""
function interp2cart(si, source::SourceTree)
    bbox = container(source)
    xc   = center(bbox)
    h    = radius(bbox)
    interp2cart(si, xc, h)
end

function interp2cart(si, xc::SVector{N}, h::Number) where {N}
    if N == 2
        s, θ = si
        @assert -π ≤ θ < π
        r = h / s
        sin_θ,cos_θ = sincos(θ)
        x,y = r*cos(θ), r*sin(θ)
        return SVector(x, y) + xc
    elseif N == 3
        s, θ, ϕ     = si
        r           = h / s
        sin_ϕ,cos_ϕ = sincos(ϕ)
        sin_θ,cos_θ = sincos(θ)
        x           = r * cos_ϕ * sin_θ
        y           = r * sin_ϕ * sin_θ
        z           = r * cos_θ
        return SVector(x, y, z) + xc
    else
        notimplemented()
    end
end

"""
    cart2interp(x,c,h)
    cart2interp(x,source::SourceTree)

Map cartesian coordinates `x` into interpolation coordinates `s` centered at `c`
and scaled by `h`. Interpolation coordinates are similar to hyperspherical
coordinates (e.g. polar in two-dimensions and spherical in three-dimensions),
except that the radial coordinates is replaced by `s=h/r`.

If passed a `SourceTree`, use its `center` and radius as `c` and `h` respectively.
"""
@inline function cart2interp(x, source::SourceTree)
    bbox = container(source)
    c    = center(bbox)
    h    = radius(bbox)
    cart2interp(x, c, h)
end

@inline function cart2interp(x, c::SVector{N}, h) where {N}
    xp = x - c
    if N == 2
        x, y = xp
        r = sqrt(x^2 + y^2)
        θ = atan(y, x)
        s = h / r
        return SVector(s, θ)
    elseif N == 3
        x, y, z = xp
        ϕ = atan(y, x)
        a = x^2 + y^2
        θ = atan(sqrt(a),z)
        r = sqrt(a + z^2)
        s = h / r
        # s = h * invsqrt(a + z^2)
        return SVector(s, θ, ϕ)
    else
        notimplemented()
    end
end

"""
    initialize_source_tree(;points,datatype,splitter,Xdatatype)

Create the tree-structure for clustering the source `points` using the splitting
strategy of `splitter`, returning a `SourceTree` built through the empty
constructor (i.e. `data = SourceTreeData()`).

# Arguments
- `points`: vector of source points.
- `splitter`: splitting strategy.
- `Xdatatype`: type container of target points.
- `Vdatatype`: type of value stored at interpolation nodes
"""
function initialize_source_tree(; Ypoints::Vector{SVector{N,T}}, splitter) where {N,T}
    source_tree_datatype = SourceTreeData{N,T}
    source_tree = ClusterTree{source_tree_datatype}(Ypoints, splitter)
    return source_tree
end

"""
    compute_interpolation_domain(source::SourceTree)

Compute the minimal domain (as a `HyperRectangle`) in interpolation space for which
`source` must provide a covering through its cone interpolants.
"""
function compute_interpolation_domain(source::SourceTree{N,Td},::Val{P}) where {N,Td,P}
    lb = svector(i->typemax(Td),N) # lower bound
    ub = svector(i->typemin(Td),N) # upper bound
    # nodes on far_list
    for far_target in farlist(source)
        for x in elements(far_target) # target points
            s  = cart2interp(x,source)
            lb = min.(lb,s)
            ub = max.(ub,s)
        end
    end
    refnodes = cheb2nodes(P)
    if !isroot(source)
        source_parent = parent(source)
        for I in active_cone_idxs(source_parent)
            rec  = cone_domain(source_parent,I)
            lc,hc = low_corner(rec), high_corner(rec)
            c0 = (lc+hc)/2
            c1 = (hc-lc)/2
            cheb  = map(x-> c0 + c1 .* x,refnodes)
            for si in cheb
                # interpolation node in physical space
                xi = interp2cart(si,source_parent)
                s  = cart2interp(xi,source)
                lb = min.(lb,s)
                ub = max.(ub,s)
            end
        end
    end
    # create the interpolation domain
    domain = HyperRectangle(lb,ub)
    return domain
end

"""
    cone_index(x::SVector,tree::SourceTree)

Given a point `x` in cartesian coordinates, return the index `I` of the cone in
`source` to whcih `x` belongs, as well as its parametric coordinate `s`.

## See also: [`cart2interp`](@ref)
"""
function cone_index(x::SVector, tree::SourceTree)
    m = msh(tree)  # UniformCartesianMesh
    s = cart2interp(x, tree)
    I = element_index(s,m)
    return I, s
end

"""
    element_index(s::SVector,m::UniformCartesianMesh)

Given a point `s ∈ m`, return the index `I` of the element in `m` containing
`s`.
"""
function element_index(s::SVector{N},m::UniformCartesianMesh{N}) where {N}
    els = ElementIterator(m)
    sz = size(els)
    Δs = step(m)
    lc = low_corner(m)
    uc = high_corner(m)
    # assert that lc <= s <= uc?
    # @assert all(lc .<= s .<= uc)
    I = ntuple(N) do n
        q = (s[n] - lc[n]) / Δs[n]
        i = ceil(Int, q)
        clamp(i, 1, sz[n])
    end
    return CartesianIndex(I)
end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
"""
SARSA-λ
-------

See learner/sarsaLlearner.jl for the core implementation.

get_all_puzzles(), get_all_episodes_in_puzzle(), get_data_by_puzzle_episode()
are in rl_data_toolsforalgos.jl
"""
function getLL_SarsaLambda(data::Array{Int,2}, alpha::Float64, gamma::Float64,
                            lambda::Float64, temp::Float64;
                            reward_magnitude=1., action_cost = 0., init_Q=0.,
                            initQSA=nothing, use_replacing_trace=true)
    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_SarsaLambda supports only one subject")
    end

    LL_unbiased = 0.

    action_selection_probs = zeros(Float64, nr_samples)

    # signals:
    RPE_record = zeros(Float64, nr_samples)
    Qsa_record = zeros(Float64, nr_samples) # trial-by-trial action value
    Vs_record = zeros(Float64, nr_samples)  # trial-by-trial max(Q(s,:))
    finalQSA_record = Dict{Int,Any}() # end of Episode, Q(sa) estimate
    RPE_QV_record = zeros(Float64, nr_samples) # the RPE defined as  Q(s_t, a_t)  -  gamma*V(s_t+1), where V(s_t+1) = max(Q(s_t+1, :))

    # Record various values (not used eventually in paper)
    state_on_counter_V1 = 1; state_on_counter_V2 = 1; state_on_counter_V3 = 1
    total_nr_of_state_on = nr_samples + data_stats.nr_episodes  # actions + nr_goal states (no action at goal)
    StateValue_V1_record = zeros(Float64, total_nr_of_state_on)  # value = max(Q(s,:)), 0 at goal state
    StateValue_V2_record = zeros(Float64, total_nr_of_state_on)  # r + gamma*V_next
    StateValue_V3_record = zeros(Float64, total_nr_of_state_on)  # V AFTER update. Start state: V before update

    PolicyParametersImbalance_record = zeros(Float64, nr_samples)
    PolicyImbalance_record = zeros(Float64, nr_samples)
    PolicyParametersImbalanceAllstates_record = Array{Array{Float64,1}, 1}(undef, nr_samples)
    [PolicyParametersImbalanceAllstates_record[i] = zeros(Float64, nr_states) for i in 1:nr_samples]
    PolicyParametersDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)
    PolicyDiffChosenVsUnchosen_record = zeros(Float64, nr_samples)

    current_puzzle_id = -1
    # Puzzle = independent block of the experiment.
    # The code was made for general purposes, in this experiment we
    # have only ONE puzzle.
    is_start_state = true
    isendofexperiment = false
    modulatingfactor = 0.
    learner = []
    for i = 1:nr_samples
        if current_puzzle_id != data[i,Int(c_puzzle)]
            current_puzzle_id = data[i,Int(c_puzzle)]
            # if new puzzle -> reset learning
            if initQSA == nothing
                learner = SarsaL(nr_states = nr_states,
                                    nr_actions = nr_actions,
                                    α = alpha, γ = gamma,
                                    λ = lambda, initvalue = init_Q,
                                    use_replacing_trace = use_replacing_trace
                                    )
            else
                learner = SarsaL(nr_states = nr_states,
                                    nr_actions = nr_actions,
                                    α = alpha,
                                    γ = gamma,
                                    λ = lambda,
                                    QSA = copy(initQSA[current_puzzle_id]),
                                    use_replacing_trace = use_replacing_trace
                                    )
            end
        end
        state = data[i, Int(c_state)]
        action = data[i, Int(c_action)]
        reward::Float64= data[i, Int(c_reward)]
        reward *= reward_magnitude
        next_state = data[i, Int(c_nextState)]

        Qsa_record[i] = copy(learner.QSA[state, action])
        Vs_record[i] = copy(maximum(learner.QSA[state, :]))

        StateValue_V1_record[state_on_counter_V1] = maximum(learner.QSA[state, :])
        state_on_counter_V1 += 1
        StateValue_V2_record[state_on_counter_V2] = learner.γ*maximum(learner.QSA[state, :])
        state_on_counter_V2 += 1

        # and if goal state:
        if reward>0
            StateValue_V1_record[state_on_counter_V1]  = 0
            state_on_counter_V1 += 1
            StateValue_V2_record[state_on_counter_V2] = reward #+ action_cost
            state_on_counter_V2 += 1
        end

        if is_start_state
            StateValue_V3_record[state_on_counter_V3] = maximum(learner.QSA[state, :])
            state_on_counter_V3 += 1
        end

        p_actions = getSoftmaxProbs(learner.QSA[state,:], temp, [0., 0.])
        action_selection_probs[i] = copy(p_actions[action])
        LL_unbiased += log(p_actions[action])

        PolicyParametersImbalance_record[i] = abs(learner.QSA[state, 1] - learner.QSA[state, 2])
        PolicyImbalance_record[i] = abs(p_actions[1] - p_actions[2])
        PolicyParametersImbalanceAllstates_record[i] = abs.(learner.QSA[:, 1] - learner.QSA[:, 2])
        PolicyParametersDiffChosenVsUnchosen_record[i] = learner.QSA[state, action] - learner.QSA[state, mod(action,2)+1]
        PolicyDiffChosenVsUnchosen_record[i] = p_actions[action] - p_actions[mod(action,2)+1]

        # Sort out if end of episode or end of experiment
        next_action = 0
        if reward==0 #non goal state
            rpe_qv = learner.γ * maximum(learner.QSA[next_state, :]) - learner.QSA[state, action]
            if nr_samples > i
                next_action = data[i+1,Int(c_action)]
            else
                # println("i:$i \nS:$state | a:$action | R:$reward | S':$next_state")
                isendofexperiment = true
            end
        end

        update!(learner, state, action, next_state, next_action, reward,
                isendofexperiment, modulatingfactor)

        if reward>0; rpe_qv = copy(learner.rpe[1]); end

        RPE_record[i] = copy(learner.rpe[1])
        RPE_QV_record[i] = copy(rpe_qv)

        finalQSA_record[current_puzzle_id] = copy(learner.QSA)

        StateValue_V3_record[state_on_counter_V3] = maximum(learner.QSA[state, :])
        state_on_counter_V3 += 1
        is_start_state = reward>0 # next state will be a start state
    end

    maxLL = LL_unbiased
    trial_by_trial_p_sa = (action_selection_probs)

    trace_type = use_replacing_trace ? "ReplTrc" : "EligTrc"
    algo = "Sarsa-Lambda ($trace_type)"

    param_v = [alpha, gamma, lambda, temp,reward_magnitude, action_cost, init_Q]
    param_str = "alpha=$(round(alpha,digits=4)), gamma=$(round(gamma,digits=4)),
                lambda=$(round(lambda,digits=4)), temp=$(round(temp,digits=4)),
                reward_magnitude=$(round(reward_magnitude,digits=4)),
                action_cost=$(round(action_cost,digits=4)),
                initQSA=$(repr(initQSA)), use_replacing_trace=$use_replacing_trace"

    results_dict = Dict{String,Any}("trialByTrial_Qsa" => Qsa_record,
                                    "final_QSA" => finalQSA_record,
                                    "trialByTrial_Vs" => Vs_record,
                                    "RPE_QV" => RPE_QV_record,
                                    "StateValue_V1_record" => StateValue_V1_record,
                                    "StateValue_V2_record" => StateValue_V2_record,
                                    "StateValue_V3_record" => StateValue_V3_record,
                                    "PolicyParametersImbalance" => PolicyParametersImbalance_record,
                                    "PolicyImbalance" => PolicyImbalance_record,
                                    "PolicyParametersImbalanceAllstates" => PolicyParametersImbalanceAllstates_record,
                                    "PolicyParametersDiffChosenVsUnchosen" => PolicyParametersDiffChosenVsUnchosen_record,
                                    "PolicyDiffChosenVsUnchosen" => PolicyDiffChosenVsUnchosen_record)

    signals = Step_by_step_signals(data_stats.all_subject_ids[1],
    RPE_record, #RPE
    [], #SPE
    trial_by_trial_p_sa,
    maxLL,
    data_stats,
    algo,
    param_v,
    param_str,
    99, #biased_action_idx, dummy value
    results_dict)
    return maxLL, signals
end


text:
function empty_xaamdi_table()
    (E                = empty_col(uenergy ),
     TotalAttenuation = empty_col(umassatt),
     EnergyLoss       = empty_col(umassatt))
end

const XAAMDITable = typeof(empty_xaamdi_table())
struct XAAMDIData <: DataSource
    element_tables::Vector{XAAMDITable}
    compound_tables::Dict{Symbol, XAAMDITable}
end

emptytable(::Type{XAAMDIData}) = empty_xaamdi_table()

const XAAMDI = load(XAAMDIData, Zs=1:92, compounds=[:H2O], dir=datapath("XAAMDI"))

default_data_source(mat, E, ::EnergyLoss) = XAAMDI

function get_table(s::XAAMDIData, key::Symbol)::XAAMDITable
    res = get(s.compound_tables, key, nothing)
    if res === nothing
        get_element_table(s, key)
    else
        res
    end
end


text:
# Main file for solving the J-SAA and J+-SAA problems

using Distributions, StatsFuns, StatsBase, JuMP, Gurobi, GLMNet

const gurobi_env = Gurobi.Env()

# read these as inputs from the corresponding "instance file"
modRepNum = parse(ARGS[1])
degNum = parse(ARGS[2])
dataRepNum = parse(ARGS[3])
sampSizeNum = parse(ARGS[4])

#*========= MODEL INFORMATION ==========
const caseNum = 2
const paramsFile = "params_case" * string(caseNum) * "_jsaa.jl"
#*======================================

#*========= INCLUDE FILES ====================
include("maxParameters.jl")
include("randomSeeds.jl")

include(paramsFile)
# the next two data files are defined in paramsFile
include(modelFile)
include(dataFile)

include("regressionMethods.jl")
include("genJSAAScenarios.jl")
include("solveModel.jl")
include("estimateSolnCost.jl")
include("evaluateSoln.jl")
#*======================================

# set seed for reproducibility
srand(startingSeed)

# directory name for storing results
const baseDirName = "case" * string(caseNum) * "_jsaa/" * "mod_" * string(modRepNum) * "/" * "deg_" * string(degNum) * "/"
const subDirName = baseDirName * "rep_" * string(dataRepNum) * "/"
const subDirName2 = subDirName * regressionMethod * "/" * "samp_" * string(sampSizeNum) * "/"
mkpath(subDirName2)

#*========= PRINT OPTIONS ====================
const storeResults = true

const infoFile = "ddsp.txt"
const modelDataFile = "model_data.txt"
const numSampFile = "num_samples_" * regressionMethod * ".txt"

const jsaaObjFile = "jsaa_obj.txt"
const jsaaDDObjFile = "jsaa_ddobj.txt"
const jpsaaObjFile = "jpsaa_obj.txt"
const jpsaaDDObjFile = "jpsaa_ddobj.txt"
const jsaaTimeFile = "jsaa_time.txt"

const covRealFile = "covariate_obs.txt"
#*======================================

#*========= GENERATE MODEL PARAMETERS ====================

srand(randomSeeds_models[modRepNum])

covariate_mean, covariate_covMat, coeff_true = generateBaseDemandModel(numCovariates)

#*========= STORE RESULTS ====================
if(storeResults)

	# write details to text file, including some key details about the test instance
	details_file = baseDirName * infoFile
	open(details_file, "w") do f
		write(f,"case number: $caseNum \n")
		write(f,"numResources: $numResources \n")
		write(f,"numCustomers: $numCustomers \n")
		write(f,"model replicate number: $modRepNum \n")
		write(f,"degree: $(degree[degNum]) \n")
		write(f,"sample sizes: $numDataSamples \n")
		write(f,"demand_errors_scaling: $demand_errors_scaling \n")
		write(f,"regressionMethod: $regressionMethod \n\n")

		write(f,"randomSeeds_MC: $randomSeeds_MC \n")
		write(f,"randomSeeds_models: $randomSeeds_models \n")
		write(f,"randomSeeds_data: $randomSeeds_data \n")
	end

	mod_data_file = baseDirName * modelDataFile
	open(mod_data_file, "w") do f
		write(f,"covariate_mean = $covariate_mean \n")
		write(f,"covariate_covMat = $covariate_covMat \n")
		write(f,"coeff_true = $coeff_true \n")
	end

	samp_size_file = baseDirName * numSampFile
	open(samp_size_file, "w") do f
		for i = 1:length(numDataSamples)
			write(f,"$(numDataSamples[i]) \n")
		end
	end

end
#*============================================

#*========= GENERATE DATA REPLICATE ====================

srand(randomSeeds_data[dataRepNum])

covariate_obs = generateCovariateReal(numCovariates,covariate_mean,covariate_covMat)

demand_data, covariate_data = generateDemandData(numCovariates,numDataSamples[sampSizeNum],degree[degNum],covariate_mean,covariate_covMat,coeff_true)

#*========= STORE RESULTS ====================
if(storeResults)

	cov_real_file = subDirName * covRealFile
	open(cov_real_file, "w") do f
		write(f,"covariate_obs = $covariate_obs")
	end

end
#*============================================

#*========= SOLVE THE J-SAA and J+-SAA MODELS ====================

tic()

demand_scen_jsaa, demand_scen_jpsaa = generateJSAAScenarios(demand_data,covariate_data,covariate_obs,regressionMethod)

z_soln_jsaa, objDDJSAA = solveSAAModel(demand_scen_jsaa)

jsaaObjEstimates = estimateSolnQuality(z_soln_jsaa,covariate_obs,degree[degNum],numMCScenarios,numMCReplicates,coeff_true)

z_soln_jpsaa, objDDJpSAA = solveSAAModel(demand_scen_jpsaa)

jpsaaObjEstimates = estimateSolnQuality(z_soln_jpsaa,covariate_obs,degree[degNum],numMCScenarios,numMCReplicates,coeff_true)

jsaaTime = toq()

#*========= STORE RESULTS ====================
if(storeResults)

	# write data to text file
	obj_est_file1 = subDirName2 * jsaaObjFile
	open(obj_est_file1, "w") do f
		for i = 1:length(jsaaObjEstimates)
			write(f,"$(jsaaObjEstimates[i]) \n")
		end
	end

	obj_est_file2 = subDirName2 * jpsaaObjFile
	open(obj_est_file2, "w") do f
		for i = 1:length(jpsaaObjEstimates)
			write(f,"$(jpsaaObjEstimates[i]) \n")
		end
	end

	ddobj_est_file1 = subDirName2 * jsaaDDObjFile
	open(ddobj_est_file1, "w") do f
		write(f,"$objDDJSAA")
	end

	ddobj_est_file2 = subDirName2 * jpsaaDDObjFile
	open(ddobj_est_file2, "w") do f
		write(f,"$objDDJpSAA")
	end

	time_obj_file = subDirName2 * jsaaTimeFile
	open(time_obj_file, "w") do f
		write(f,"$jsaaTime")
	end

end
#*============================================

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
module TestRadialODE

export compare_to_qnm

const tolerance = 1e-3 ## tolerance we compare to

include("../src/CustomTypes.jl")
include("../src/RadialODE.jl")
include("ReadQNM.jl")
using .CustomTypes
using .ReadQNM
import .RadialODE as RO
import Test: @test

##============================================================

"""
Compare against values computed by Leo Stein's qnm code,

J.Open Source Softw. 4 (2019) 42, 1683

This is written in python so we need to bind to the library
"""
function compare_to_qnm(
      nr::myI,
      s::myI,
      n::myI,
      l::myI,
      m::myI,
      avals::Vector{myF}
   )

   println("Comparison to qnm library: s=$s, n=$n, l=$l, m=$m")

   for a in avals
      om, la = qnm(n,s,m,l,a)

      ls_c, vs_c = RO.eig_vals_vecs_c(nr, s, m, a, tomyC(om))

      #println("testing: a=$a, ω=$om\nΛ=$la, Cheb Λ=$(ls_c[1])")

      @test abs(ls_c - la)/max(1,abs(la))<tolerance
   end
end

end


text:
module SBML

using EzXML, MathML, ModelingToolkit, Symbolics

include("parse.jl")
include("sbml2reactionsystem.jl")

export sbml_to_sysinfo, sbml_lists

export SbmlModel
export build_map, build_par_map, build_comp_map, build_spec_map,
       promotelocalparameters, getkineticlaw, getmath, build_reactions, to_table, uniquekeys, get_name, make_nice

end


text:
struct AggressivenessPOWNodeBelief{P<:MLPOMDP}
    model::P
    b::AggressivenessBelief
    r::Float64
end

rand(rng::AbstractRNG, b::AggressivenessPOWNodeBelief) = (rand(rng, b.b), b.r)

struct AggressivenessPOWFilter
    params::WeightUpdateParams
end

POMCPOW.belief_type(::Type{AggressivenessPOWFilter}, ::Type{P}) where P<:MLPOMDP = AggressivenessPOWNodeBelief{P}

function POMCPOW.init_node_sr_belief(f::AggressivenessPOWFilter, p::MLPOMDP, s, a, sp, o, r)
    particles = [Vector{Float64}() for c in o.cars]
    cweights = [Vector{Float64}() for c in o.cars]
    gen = CorrelatedIDMMOBIL(p.dmodel.behaviors)
    maybe_push_one!(particles, cweights, f.params, p.dmodel.phys_param, gen, sp, o)
    ab = AggressivenessBelief(gen, o, particles, cweights)
    return AggressivenessPOWNodeBelief(p, ab, r)
end

function POMCPOW.push_weighted!(b::AggressivenessPOWNodeBelief, f::AggressivenessPOWFilter, s, sp, r)
    maybe_push_one!(b.b.particles, b.b.cweights, f.params, b.model.dmodel.phys_param, b.b.gen, sp, b.b.physical)
end

actions(p::MLPOMDP, b::AggressivenessPOWNodeBelief) = actions(p, b.b.physical)

function actions(p::MLPOMDP, h::POWTreeObsNode)
    if isroot(h)
        return actions(p, belief(h))
    else
        return actions(p, sr_belief(h))
    end
end

struct BehaviorPOWNodeBelief{P<:MLPOMDP,G}
    model::P
    b::BehaviorParticleBelief{G}
    r::Float64
end

rand(rng::AbstractRNG, b::BehaviorPOWNodeBelief) = (rand(rng, b.b), b.r)

struct BehaviorPOWFilter
    params::WeightUpdateParams
end

POMCPOW.belief_type(::Type{BehaviorPOWFilter}, ::Type{P}) where P<:MLPOMDP = BehaviorPOWNodeBelief{P, gen_type(P)}

function POMCPOW.init_node_sr_belief(f::BehaviorPOWFilter, p::MLPOMDP, s, a, sp, o, r)
    particles = [Vector{IDMMOBILBehavior}() for c in o.cars]
    cweights = [Vector{Float64}() for c in o.cars]
    gen = p.dmodel.behaviors
    maybe_push_one!(particles, cweights, f.params, p.dmodel.phys_param, gen, sp, o)
    ab = BehaviorParticleBelief(gen, o, particles, cweights)
    return BehaviorPOWNodeBelief(p, ab, r)
end

function POMCPOW.push_weighted!(b::BehaviorPOWNodeBelief, f::BehaviorPOWFilter, s, sp, r)
    maybe_push_one!(b.b.particles, b.b.cweights, f.params, b.model.dmodel.phys_param, b.b.gen, sp, b.b.physical)
end

actions(p::MLPOMDP, b::BehaviorPOWNodeBelief) = actions(p, b.b.physical)


text:
"""
    worker(balancer<:AbstractLoadBalancer)

Start a worker that receives work messages and performs computation asynchronously.
"""
function worker(balancer::T) where {T<:AbstractLoadBalancer}

    local_chl = Channel{Message}(64)

    @sync begin
        # MESSAGE HANDLER SUBTASK
        @async _msg_handler(balancer, local_chl)
        # SUBTASK 1
        @async do_work(balancer, local_chl)
    end
    put!(balancer.stat_chl, Message(:done, myid()))
end

"""
    do_work(balancer<:AbstractLoadBalancer, local_chl::Channel{Message})

Pull work from local_chl to compute, and send status updates to msg_chl.

This worker receives two message kinds (<| `local_chl`):
- `:work`
- `:end` - This worker will exit upon receiving this message

This worker also produces two message kinds (|> `balancer.msg_chl`)
- `:_idle` - Signal that this worker is idle (`local_chl` empty)
- `:_nonidle` - Signal that this worker is nonidle
"""
function do_work(balancer::T,
                 local_chl::Channel{Message}) where {T<:AbstractLoadBalancer}

    idle::Bool = true
    my_msg_chl = get_msg_chl(balancer)

    info("worker ", string(myid()), " do_work")

    while true
        # wait for an item to be *appended* (going from i -> i+1, not necessary from 0 -> 1)
        # REMARK: Could this block indefinitely if :end arrives before this call to wait() and,
        #         since no messages are expected to arrive after :end, this waits for another
        #         value to be appended. Should we assert that the channel is empty at this point?
        # REMARK: One way of solving this could be to have the msg_handler signal on the local
        #         channels put_cond... perhaps
        @assert length(local_chl.data) == 0
        wait(local_chl)

        last_signal = time()

        # while there are messages in the local channel
        while isready(local_chl)

            local msg = take!(local_chl)

            # if this message is end, exit
            if msg.kind == :end
                put!(balancer.res_chl, Message(:end, myid()))
                return
            elseif msg.kind == :work
                # if this worker has been idle, signal it is no longer idle
                if idle
                    info("worker ", string(myid()), " working")
                    idle = false
                    last_signal = time()
                    put!(my_msg_chl, Message(:_nonidle, myid()))
                # if it has been 1 second since the last nonidle message

                elseif time() - last_signal > 1
                    put!(my_msg_chl, Message(:_nonidle, myid()))
                    last_signal = time()
                end

                # do work
                local work = msg._data2

                if work.units > 1
                    # split the work
                    work_1, work_2 = split_work(work)

                    # place the extra work back in the local channel
                    put!(local_chl, Message(:work, myid(), work_2))
                else
                    work_1 = work
                end

                # simulate calculation
                sleep(work_1.unitcost)

                # take 1 down, pass it around
                work_1.units -= 1
                if work_1.units > 0
                    # place back in local channel
                    put!(local_chl, Message(:work, myid(), work_1))
                end

                put!(balancer.res_chl, Message(:work, myid(), work_1.unitcost))
            else
                err("_worker received unrecognized message! ", string(msg))
                error("do_work error")
            end
        end

        # local_chl is now empty
        if !idle
            idle = true
            put!(my_msg_chl, Message(:_idle, myid()))
        end
    end
end # end _worker function

"""
    send_work(balancer<:AbstractLoadBalancer, local_chl::Channel{Message}, msg::Message)

Attempt to send a piece of work from this worker to another.
"""
function send_work(balancer::T,
                   local_chl::Channel{Message},
                   msg::Message) where {T<:AbstractLoadBalancer}
    # attempt to move some local work to the remote worker
    other_wid = msg.data
    other_msg_chl = get_msg_chl(balancer, other_wid)
    if isready(local_chl)
        if fetch(local_chl).kind == :end
            return
        end
        work = take!(local_chl)
        put!(other_msg_chl, work)
    else
        put!(other_msg_chl, Message(:nowork, myid()))
    end
end

max_stars_repo_name: dillondaudert/jlance

max_issues_repo_name: dillondaudert/jlance

max_forks_repo_name: dillondaudert/jlance


text:
#=
The "Resynthesizer" algorithm for texture transfer in images.

Copyright 2021 Lloyd Konneker

Original algorithm by Paul Harrison.
=#

#=
This project is an exercise in translating and refactoring algorithms.
Code derived from https://github.com/bootchk/resynthesizer.git, written in C.
Specifically from the  /lib directory, the library implementing the algorithm.

One main lesson is that the original C code is very hard to understand.
(Even though I wrote the C code, but I didn't devise the original algorithm)
-The separate declarations of variables is confusing.
-The large functions are hard to understand.
-Better to name functions properly than to use comments.
-Hard to understand data flow when pointers are used for OUT variables
(although Julia doesn't ease that when you pass mutable objects.)
-Better to pass structs of related values

The derivation is mostly a translation from language C to Julia.
But also refactoring:
- renamed things to better denote what they are, using Smalltalk like naming
- separated functions into smaller files

Retains the original programming strategy of a reentrant C library,
where there are no globals and everything is passed in function arguments.
Except that parameters are global constants.

Also retains the architecture of the original:
a basic, parameterized algorithm, wrapped in use cases of the algorithm.
For example, the basic algorithm takes two images,
but a use case may take one image and derive two images
to pass to the basic algorithm.

Uses many functions built-in to Julia e.g. findall and colordiff.
The original implemented many of those functions from scratch.

Also uses Julia's bounds checking to good advantage,
instead of implementing clipping algorithms from scratch.

The method of development: test based.
After each small migration of code, retest.
>cd tests
>julia julia --project=. resynthesizer.jl
When progress is made and the produced image changes from prior tests,
update reference image interactively:
>julia
julia> include("resynthesizer.jl")
Expect a dialog to let you update the reference test, i.e. the gold standard result.
=#

#=
Another lesson is in "generic" algorithms.
While refactoring, I attempt to generalize to multidimension.
The existing algorithm is specialized to 2D array of color.
Julia is all about writing generic algorithms.
While refactoring, I attempt to maintain genericity.
So that the algorithm could work with  multidimension array AKA tensor,
of differing element type than "color".
IOW the algorithm might now work for the applications:
synthesize missing text, synthesize video, synthesize music.

The algorithm is generic on:
- ValueType: the type of the points being compared in the innermost loops of the algorithm
- DimensionCount: the dimensions of the tensors (e.g. Images are 2D)
So these type variables flow through the code.
When you open a specific tensor file (e.g. an image),
that specifies these type variables, and concretizes many abstract types.
Then Julia compiles the algorithm for those concrete types.

TODO Substitute "tensor" or "MDA multidimensional array" for "image"
Don't use the word "image" except as an example.

TODO In many places, iterate over the dimensions of the tensor,
instead of using just two dimensions.

TODO test harness for text, video, music

My experience has been:
it is relatively easy to specify as generic,
but hard to get performance of the original implementation.
I struggled with too many allocations and with abstract versus concrete types.
I started without type annotations, hoping they would be inferred,
then added many type annotations.
The takeaway might be: it is better to supply type annotations
(even if they could be inferred)
because it is not easy to understand what the inferred types are
(without a lot of testing and harnessing the code
and using code_warntype)
=#

include("image.jl")
include("passes/passes.jl")
include("result/synthResult.jl")
include("result/searchResult.jl")

#include("scatterPatch.jl")
include("parameters.jl")

# global
const parameters = ResynthesizerParameters()

const metric = DE_94()

#=

Note that the original implementation used the word "source" for the target image.
Also, the original used the word "target" ambiguously to mean the masked region of the target.
=#

#=
The applications accepts many abstract types,
but internally uses more concrete types.
=#
function resynthesize(
    targetImage::MaskedImage{ValueType, DimensionCount},  # AbstractArray{T,N},
    corpusImage::MaskedImage{ValueType, DimensionCount}    # AbstractArray{Bool,N},
    ) where {ValueType, DimensionCount}

    #=
    The mask on the target defines the synth region to be possibly mutated.
    The mask on the corpus defines the selected region to sample from.
    =#

    @debug "Image length, size" length(targetImage.image) size(targetImage.image)

    #=
    Seed the RNG so that the random number sequence is reproducible i.e. deterministic.

    Not strictly necessary, but it aids in debugging so that tests are more reproducible.
    =#
    Random.seed!(123)

    # uninitialized mutable struct
    synthPatch = ScatterPatch(targetImage)

    # initial fundamental, overall result
    synthResult = initialSynthResult(targetImage)

    # Initial, singleton searchResult, allocated once and reused by passes
    searchResult = SearchResult(targetImage.image)

    # offsets to span the targetImage
    sortedOffsets = SortedOffsets(targetImage.image)

    # TODO use toggled_asserts
    # This is temporary, while chasing performance
    @assert isconcretetype(typeof(corpusImage))
    @assert isconcretetype(typeof(synthPatch))
    @assert isconcretetype(typeof(synthResult))
    @assert isconcretetype(typeof(sortedOffsets))

    #=
    Nullify the synth region.
    For visual debugging when images, when animated, or at the end too much black is bad.
    The algorithm does not require this.
    =#
    nullifyMaskedRegion(targetImage)

    makePassesUntilGoodEnough(targetImage, corpusImage, synthPatch, synthResult, searchResult, sortedOffsets)
    # assert targetImage is mutated and synthResult is mutated

    #=
    The result is the possibly mutated original (sans mask).
    Possibly: depends on parameters.withReplacement.
    Assert that only the masked region of the original is mutated.
    =#

    # TODO return the synthResult
    return targetImage.image
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using Gadfly

samplesize = 100
coord = Coord.cartesian(xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0)
xsc = Scale.x_continuous(minvalue=0.0,maxvalue=0.0)
ysc = Scale.y_continuous(minvalue=0.0,maxvalue=1.0)
colsc = Scale.color_continuous(minvalue=-0,maxvalue=1)
Data = DataFrame(X=map(sin,map(acos, session_similarity(1:samplesize))) , Y=session_similarity(1:samplesize))
Data[:slope] = Data[:Y] ./ Data[:X] # rise over run

layerAB = layer(Data, x=repeat([0], inner=samplesize), y=repeat([0], inner=samplesize), xend=:X, yend=:Y, color=:Y, Geom.segment)

p = plot(layerAB, xsc, ysc, colsc, coord)
draw(PNG("cossim.png", 6inch, 6inch), p)


text:
module HTMLWriterTests

using Test

import JSON
import Documenter.Writers.HTMLWriter: generate_version_file, expand_versions

function verify_version_file(versionfile, entries)
    @test isfile(versionfile)
    content = read(versionfile, String)
    idx = 1
    for entry in entries
        i = findnext(entry, content, idx)
        @test i !== nothing
        idx = last(i)
    end
end

@testset "HTMLWriter" begin
    mktempdir() do tmpdir
        versionfile = joinpath(tmpdir, "versions.js")
        versions = ["stable", "dev",
                    "2.1.1", "v2.1.0", "v2.0.1", "v2.0.0",
                    "1.1.1", "v1.1.0", "v1.0.1", "v1.0.0",
                    "0.1.1", "v0.1.0"] # note no `v` on first ones
        cd(tmpdir) do
            for version in versions
                mkdir(version)
            end
        end

        # expanding versions
        versions = ["stable" => "v^", "v#.#", "dev" => "dev"] # default to makedocs
        entries, symlinks = expand_versions(tmpdir, versions)
        @test entries == ["stable", "v2.1", "v2.0", "v1.1", "v1.0", "v0.1", "dev"]
        @test symlinks == ["stable"=>"2.1.1", "v2.1"=>"2.1.1", "v2.0"=>"v2.0.1",
                           "v1.1"=>"1.1.1", "v1.0"=>"v1.0.1", "v0.1"=>"0.1.1",
                           "v2"=>"2.1.1", "v1"=>"1.1.1", "v2.1.1"=>"2.1.1",
                           "v1.1.1"=>"1.1.1", "v0.1.1"=>"0.1.1"]
        generate_version_file(versionfile, entries)
        verify_version_file(versionfile, entries)

        versions = ["v#"]
        entries, symlinks = expand_versions(tmpdir, versions)
        @test entries == ["v2.1", "v1.1"]
        @test symlinks == ["v2.1"=>"2.1.1", "v1.1"=>"1.1.1", "v2"=>"2.1.1", "v1"=>"1.1.1",
                           "v2.0"=>"v2.0.1", "v1.0"=>"v1.0.1", "v0.1"=>"0.1.1",
                           "v2.1.1"=>"2.1.1", "v1.1.1"=>"1.1.1", "v0.1.1"=>"0.1.1"]
        generate_version_file(versionfile, entries)
        verify_version_file(versionfile, entries)

        versions = ["v#.#.#"]
        entries, symlinks = expand_versions(tmpdir, versions)
        @test entries == ["v2.1.1", "v2.1.0", "v2.0.1", "v2.0.0", "v1.1.1", "v1.1.0",
                          "v1.0.1", "v1.0.0", "v0.1.1", "v0.1.0"]
        @test symlinks == ["v2.1.1"=>"2.1.1", "v1.1.1"=>"1.1.1", "v0.1.1"=>"0.1.1",
                           "v2"=>"2.1.1", "v1"=>"1.1.1", "v2.1"=>"2.1.1",
                           "v2.0"=>"v2.0.1", "v1.1"=>"1.1.1", "v1.0"=>"v1.0.1", "v0.1"=>"0.1.1"]
        generate_version_file(versionfile, entries)
        verify_version_file(versionfile, entries)

        versions = ["v^", "devel" => "dev", "foobar", "foo" => "bar"]
        entries, symlinks = expand_versions(tmpdir, versions)
        @test entries == ["v2.1", "devel"]
        @test ("v2.1" => "2.1.1") in symlinks
        @test ("devel" => "dev") in symlinks
        generate_version_file(versionfile, entries)
        verify_version_file(versionfile, entries)

        versions = ["stable" => "v^", "dev" => "stable"]
        @test_throws ArgumentError expand_versions(tmpdir, versions)
    end
end

end


text:
# Autogenerated wrapper script for MKL_jll for i686-linux-gnu
export libmkl_core, libmkl_rt

using IntelOpenMP_jll
JLLWrappers.@generate_wrapper_header("MKL")
JLLWrappers.@declare_library_product(libmkl_core, "libmkl_core.so.1")
JLLWrappers.@declare_library_product(libmkl_rt, "libmkl_rt.so.1")
function __init__()
    JLLWrappers.@generate_init_header(IntelOpenMP_jll)
    JLLWrappers.@init_library_product(
        libmkl_core,
        "lib/libmkl_core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmkl_rt,
        "lib/libmkl_rt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
using DelimiteFiles
=#

# NOTE: These simulation setting are identical to that of test-sims-6-5 EXCEPT
# that W_true is different. The W_true here has low abundance (< 10% each) for
# 2 pairs (i.e. 4) features. Thus, they should be harder to recover, and
# uncertainty in R_i should be greater.

# Simulation Name
simname = basename("$(@__DIR__)")
println(simname); flush(stdout)

# NOTE: write to scratchdir
function outdir_suffix(maxtemp, ntemps, degree)
  return "maxtemp$(maxtemp)-ntempts$(ntemps)-degree$(degree)-N500" # TODO:N2000
end

# NOTE: modify
settings = let
  Z = Matrix{Bool}(readdlm(joinpath(@__DIR__, "Z.txt"), Int, comments=true))
  W = readdlm(joinpath(@__DIR__, "W.txt"), comments=true)
  [Dict(:simname => simname,
        :repfam_dist_scale => 1.0,
        :N => [500, 500], # TODO: [2000, 2000],
        :Z => Z,
        :W => W,
        :thin_samps => 2,
        :nburn => 10000,
        :nsamps => 2000,
        :Lmcmc => Dict(0 => 2, 1 => 2),
        :Kmcmc => 10,
        :maxtemp => maxtemp,
        :ntemps => ntemps,
        :ncores => ntemps,
        :degree => degree,
        :dataseed => 1,
        :mcmcseed => 1,
        :outdir_suffix => outdir_suffix(maxtemp, ntemps, degree))
   for maxtemp in [4]
   for ntemps in [20]
   for degree in [1]]
end


text:
using Statistics,Serialization,StatsPlots, DataFrames

# paras = ["no_pol","dist","dist_ho"]

filename_prefix = "../figures/"
path_to_data = "../data/mutation/"

badpol = true # do you have data for badpol and switches?
show_capacity = false # for infected plot

# set intial values as in coivd par ini
datat = 7 # interval of data collection
T = 744 #582  # time periods

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 102
# set periods for covid plots
inf_start_day = 1
inf_end_day = 600 # cut of after period 400
inf_start_week = 1
inf_end_week = 86 # 58*7 = 406 and 42*7 = 294

datapoint = fld(T,datat)
include("../data/mutation/beta_l10/covid_par_ini.jl") # necessary for some basic parameter values (nhh)

# parameter values
paras = ["beta_l0","beta_l5","beta_l30","beta_l50", "beta_l100"]

# scatterplot
gdploss_means = []
totcas_means = []
pubacc_means = []

for (i,p) in enumerate(paras)
    #filename_prefix = string("adap_start//", p , "//")
    prefix = string(path_to_data, "/", p , "/")
    # store data scatterplot
    worker_results = deserialize(open("$(prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means,gdploss)
    push!(totcas_means,totcas)
    push!(pubacc_means,pubacc)
    if p == "beta_l5"
        global point_A_x = gdploss
        global point_A_y = totcas
    end
    if p == "beta_l50"
            global point_D_x = gdploss
            global point_D_y = totcas
        end
end
plfull = plot(gdploss_means, totcas_means, linewidth = 2, linecolor = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means, totcas_means, linewidth = 2, seriestype = :scatter, color = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
plarrows = quiver([gdploss_means[ll-1]],[totcas_means[ll-1]],quiver=([gdploss_means[ll]-gdploss_means[ll-1]],[totcas_means[ll]-totcas_means[ll-1]]), linecolor = [:black], linewidth = 2, xlabel = "average GDP loss [%]", ylabel = "mortality [%]")
gdploss_means = reverse(gdploss_means)
totcas_means = reverse(totcas_means)
ll = length(paras)
quiver!(plarrows, gdploss_means[2:ll-1],totcas_means[2:ll-1],quiver=(gdploss_means[3:ll]-gdploss_means[2:ll-1],totcas_means[3:ll]-totcas_means[2:ll-1]), linecolor = [:black], linewidth = 2)
scatter!(gdploss_means, totcas_means, markersize = 3, markercolor = [:black], legend=false)

# parameter values
#paras = ["p10_050","p30_050","p50_050"]
#
# scatterplot
#gdploss_means = []
#totcas_means = []
#pubacc_means = []
#
#for (i,p) in enumerate(paras)
#    #filename_prefix = string("adap_start//", p , "//")
#    filename_prefix = string( p , "//")
#    # store data scatterplot
#    worker_results = deserialize(open("$(filename_prefix)batchdata.dat"))
#    gdploss = mean((results -> results[:togdploss]).(worker_results))
#    totcas = mean(100*(results -> results[:totcas]).(worker_results))
#    pubacch = (results -> results[:pubacctraj]).(worker_results)
#    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
#    global nn = size(gdpmean)[1]
#    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
#    push!(gdploss_means,gdploss)
#    push!(totcas_means,totcas)
#    push!(pubacc_means,pubacc)
#end
#plot!(plfull, gdploss_means, totcas_means, linewidth = 2, linestyle = [:dot], linecolor = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
#plot!(plfull, gdploss_means, totcas_means, linewidth = 2, seriestype = [:scatter], color = [:black], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
#
#ll = length(paras)
#quiver!(plarrows, gdploss_means[1:ll-1],totcas_means[1:ll-1],quiver=(gdploss_means[2:ll]-gdploss_means[1:ll-1],totcas_means[2:ll]-totcas_means[1:ll-1]), linecolor = [:black], linestyle = [:dot], linewidth = 2)
#scatter!(gdploss_means, totcas_means, markersize = 3, markercolor = [:black], legend=false)

paras = ["beta_l5","beta_l5_alpha_o025","beta_l5_alpha_o050","beta_l5_alpha_o075","beta_l5_alpha_o100"]

# scatterplot
gdploss_means1 = []
totcas_means1 = []
pubacc_means1 = []

for (i,p) in enumerate(paras)
    #filename_prefix = string("adap_start//", p , "//")
    prefix = string(path_to_data, "/", p , "/")
    # store data scatterplot
    worker_results = deserialize(open("$(prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means1,gdploss)
    push!(totcas_means1,totcas)
    push!(pubacc_means1,pubacc)
end
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linecolor = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, seriestype = [:scatter], color = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, gdploss_means1[1:ll-1],totcas_means1[1:ll-1],quiver=(gdploss_means1[2:ll]-gdploss_means1[1:ll-1],totcas_means1[2:ll]-totcas_means1[1:ll-1]), linecolor = [:blue], linewidth = 2)
scatter!(gdploss_means1, totcas_means1, markersize = 3, markercolor = [:blue], legend=false)

paras = ["beta_l5_alpha_l050","beta_l5_alpha_l050_alpha_o025","beta_l5_alpha_l050_alpha_o050"]

# scatterplot
gdploss_means1 = []
totcas_means1 = []
pubacc_means1 = []

for (i,p) in enumerate(paras)
    #filename_prefix = string("adap_start//", p , "//")
    prefix = string(path_to_data, "/", p , "/")
    # store data scatterplot
    worker_results = deserialize(open("$(prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means1,gdploss)
    push!(totcas_means1,totcas)
    push!(pubacc_means1,pubacc)
    if p == "beta_l5_alpha_l050_alpha_o050"
        global point_E_x = gdploss
        global point_E_y = totcas
    end
end
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linecolor = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, seriestype = [:scatter], color = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, gdploss_means1[1:ll-1],totcas_means1[1:ll-1],quiver=(gdploss_means1[2:ll]-gdploss_means1[1:ll-1],totcas_means1[2:ll]-totcas_means1[1:ll-1]), linecolor = [:blue], linewidth = 2)
scatter!(gdploss_means1, totcas_means1, markersize = 3, markercolor = [:blue], legend=false)

paras = ["beta_l50","beta_l50_alpha_o025","beta_l50_alpha_o050","beta_l50_alpha_o075","beta_l50_alpha_o100"]

# scatterplot
gdploss_means1 = []
totcas_means1 = []
pubacc_means1 = []

for (i,p) in enumerate(paras)
    #filename_prefix = string("adap_start//", p , "//")
    prefix = string(path_to_data, "/", p , "/")
    # store data scatterplot
    worker_results = deserialize(open("$(prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means1,gdploss)
    push!(totcas_means1,totcas)
    push!(pubacc_means1,pubacc)
    if p == "beta_l50_alpha_o050"
        global point_C_x = gdploss
        global point_C_y = totcas
    end
end
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linestyle= [:dot], linecolor = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linestyle= [:dot], seriestype = [:scatter], color = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, gdploss_means1[1:ll-1],totcas_means1[1:ll-1],quiver=(gdploss_means1[2:ll]-gdploss_means1[1:ll-1],totcas_means1[2:ll]-totcas_means1[1:ll-1]), linecolor = [:blue], linestyle= [:dot], linewidth = 2)
scatter!(gdploss_means1, totcas_means1, markersize = 3, markercolor = [:blue], linestyle= [:dot], legend=false)

# paras = ["5","5_025","5_050","5_075"]
#
# # scatterplot
# gdploss_means1 = []
# totcas_means1 = []
# pubacc_means1 = []
#
# for (i,p) in enumerate(paras)
#     #filename_prefix = string("adap_start//", p , "//")
#     filename_prefix = string( p , "//")
#     # store data scatterplot
#     worker_results = deserialize(open("$(filename_prefix)batchdata.dat"))
#     gdploss = mean((results -> results[:togdploss]).(worker_results))
#     totcas = mean(100*(results -> results[:totcas]).(worker_results))
#     pubacch = (results -> results[:pubacctraj]).(worker_results)
#     gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
#     global nn = size(gdpmean)[1]
#     pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
#     push!(gdploss_means1,gdploss)
#     push!(totcas_means1,totcas)
#     push!(pubacc_means1,pubacc)
#     if p == "50_050"
#         global point_C_x = gdploss
#         global point_C_y = totcas
#     end
# end
# plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linestyle= [:dot], linecolor = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
# plot!(plfull, gdploss_means1, totcas_means1, linewidth = 2, linestyle= [:dot], seriestype = [:scatter], color = [:blue], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
#
# ll = length(paras)
# quiver!(plarrows, gdploss_means1[1:ll-1],totcas_means1[1:ll-1],quiver=(gdploss_means1[2:ll]-gdploss_means1[1:ll-1],totcas_means1[2:ll]-totcas_means1[1:ll-1]), linecolor = [:blue], linestyle= [:dot], linewidth = 2)
# scatter!(gdploss_means1, totcas_means1, markersize = 3, markercolor = [:blue], linestyle= [:dot], legend=false)

paras = ["beta_l5_alpha_l025","beta_l5_alpha_l050","beta_l5_alpha_l075","beta_l5","beta_l5_alpha_l125"]

# scatterplot
gdploss_means2 = []
totcas_means2 = []
pubacc_means2 = []

for (i,p) in enumerate(paras)
    #filename_prefix = string("adap_start//", p , "//")
    prefix = string(path_to_data, "/", p , "/")
    # store data scatterplot
    worker_results = deserialize(open("$(prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means2,gdploss)
    push!(totcas_means2,totcas)
    push!(pubacc_means2,pubacc)
    if p == "beta_l5_alpha_l050"
        global point_B_x = gdploss
        global point_B_y = totcas
    end
end
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, linecolor = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, seriestype= [:scatter], color = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, [gdploss_means2[ll-1]],[totcas_means2[ll-1]],quiver=([gdploss_means2[ll]-gdploss_means2[ll-1]],[totcas_means2[ll]-totcas_means2[ll-1]]), linecolor = [:red], linewidth = 2)
gdploss_means2 = reverse(gdploss_means2)
totcas_means2 = reverse(totcas_means2)
quiver!(plarrows, gdploss_means2[2:ll-1],totcas_means2[2:ll-1],quiver=(gdploss_means2[3:ll]-gdploss_means2[2:ll-1],totcas_means2[3:ll]-totcas_means2[2:ll-1]), linecolor = [:red], linewidth = 2)
scatter!(gdploss_means2, totcas_means2, markersize = 3, markercolor = [:red], legend=false)

paras = ["beta_l50_alpha_l025","beta_l50_alpha_l050","beta_l50_alpha_l075","beta_l50","beta_l50_alpha_l125"]

# scatterplot
gdploss_means2 = []
totcas_means2 = []
pubacc_means2 = []

for (i,p) in enumerate(paras)
    #filename_prefix = string("adap_start//", p , "//")
    prefix = string(path_to_data, "/", p , "/")
    # store data scatterplot
    worker_results = deserialize(open("$(prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    pubacch = (results -> results[:pubacctraj]).(worker_results)
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    pubacc = mean([pubacch[i][datapoint+1]*datat/(365*gdpmean[2]) for  i = 1:size(pubacch)[1]])
    push!(gdploss_means2,gdploss)
    push!(totcas_means2,totcas)
    push!(pubacc_means2,pubacc)
end
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, linestyle= [:dot], linecolor = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)
plot!(plfull, gdploss_means2, totcas_means2, linewidth = 2, linestyle= [:dot], seriestype= [:scatter], color = [:red], legend=false) #, linewidth = [2], linecolor = [:black], xlabel = "Average GDP Loss", ylabel = "Total Casualties", group = paras, legend =:top)

ll = length(paras)
quiver!(plarrows, [gdploss_means2[ll-1]],[totcas_means2[ll-1]],quiver=([gdploss_means2[ll]-gdploss_means2[ll-1]],[totcas_means2[ll]-totcas_means2[ll-1]]), linecolor = [:red], linestyle= [:dot], linewidth = 2)
gdploss_means2 = reverse(gdploss_means2)
totcas_means2 = reverse(totcas_means2)
quiver!(plarrows, gdploss_means2[2:ll-1],totcas_means2[2:ll-1],quiver=(gdploss_means2[3:ll]-gdploss_means2[2:ll-1],totcas_means2[3:ll]-totcas_means2[2:ll-1]), linecolor = [:red], linestyle= [:dot], linewidth = 2)
scatter!(gdploss_means2, totcas_means2, markersize = 3, markercolor = [:red], linestyle= [:dot], legend=false)

### add specific points A B C
scatter!([point_A_x], [point_A_y], markersize = 4, markercolor = [:green])
scatter!([point_A_x+ 0.15], [point_A_y], markersize = 0, markercolor = [:white], series_annotations = text.(["B"], :bottom))

scatter!([point_B_x], [point_B_y], markersize = 4, markercolor = [:green])
scatter!([point_B_x], [point_B_y-0.02], markersize = 0, markercolor = [:white], series_annotations = text.(["C"], :bottom))

#scatter!([point_C_x], [point_C_y], markersize = 4, markercolor = [:green])
#scatter!([point_C_x-0.15], [point_C_y], markersize = 0, markercolor = [:white], series_annotations = text.(["D"], :right))

scatter!([point_D_x], [point_D_y], markersize = 4, markercolor = [:green])
scatter!([point_D_x+0.15], [point_D_y+0.001], markersize = 0, markercolor = [:white], series_annotations = text.(["A"], :bottom))

scatter!([point_E_x], [point_E_y], markersize = 4, markercolor = [:green])
scatter!([point_E_x+.15], [point_E_y], markersize = 0, markercolor = [:white], series_annotations = text.(["D"], :left))

ylims!(plarrows,0.0,0.25)
xlims!(plarrows,0,6)

savefig(plarrows,"$(filename_prefix)variation_policy_mutation.pdf")


text:
using Individual.SchemaBase
using Catlab.ACSetInterface
using Test

@testset "construct model with integer initial states" begin

    initial_states = [1,2,3,1,2,3]
    state_labels = ["S", "I", "R"]

    SIR = IBM{String}()
    initialize_states(SIR, initial_states, state_labels)

    @test nparts(SIR, :Person) == length(initial_states)
    @test nparts(SIR, :State) == length(state_labels)
    @test subpart(SIR, :statelabel) == state_labels
    @test subpart(SIR, :state) == initial_states
    @test subpart(SIR, :state_update) == zeros(Int64, length(initial_states))

    @test npeople(SIR) == length(initial_states)
    @test npeople(SIR, "S") == count(x -> x == 1, initial_states)
    @test npeople(SIR, "I") == count(x -> x == 2, initial_states)
    @test npeople(SIR, "R") == count(x -> x == 3, initial_states)
    @test nstate(SIR) == length(state_labels)
    @test statelabel(SIR) == state_labels
    @test get_index_state(SIR) == 1:length(initial_states)
    @test get_index_state(SIR, "S") == findall(x -> x == 1, initial_states)
    @test get_index_state(SIR, "I") == findall(x -> x == 2, initial_states)
    @test get_index_state(SIR, "R") == findall(x -> x == 3, initial_states)

end

@testset "construct model with string initial states" begin

    initial_states = ["S", "I", "R", "S", "I", "R"]
    state_labels = ["S", "I", "R"]

    SIR = IBM{String}()
    initialize_states(SIR, initial_states, state_labels)

    @test nparts(SIR, :Person) == length(initial_states)
    @test nparts(SIR, :State) == length(state_labels)
    @test subpart(SIR, :statelabel) == state_labels
    @test subpart(SIR, :state) == indexin(initial_states, state_labels)
    @test subpart(SIR, :state_update) == zeros(Int64, length(initial_states))

    @test npeople(SIR) == length(initial_states)
    @test npeople(SIR, "S") == count(x -> x == "S", initial_states)
    @test npeople(SIR, "I") == count(x -> x == "I", initial_states)
    @test npeople(SIR, "R") == count(x -> x == "R", initial_states)
    @test nstate(SIR) == length(state_labels)
    @test statelabel(SIR) == state_labels
    @test get_index_state(SIR) == 1:length(initial_states)
    @test get_index_state(SIR, "S") == findall(x -> x == "S", initial_states)
    @test get_index_state(SIR, "I") == findall(x -> x == "I", initial_states)
    @test get_index_state(SIR, "R") == findall(x -> x == "R", initial_states)

end

@testset "construct model errors properly" begin

    initial_states = [1,2,3,1,2,3,4]
    state_labels = ["S", "I", "R"]

    SIR = IBM{String}()
    @test_throws ArgumentError initialize_states(SIR, initial_states, state_labels)

    initial_states = ["S", "I", "R", "S", "I", "R", "X"]
    @test_throws ArgumentError initialize_states(SIR, initial_states, state_labels)

end

@testset "queuing and applying state updates works" begin

    initial_states = ["S", "I", "R", "S", "I", "R"]
    state_labels = ["S", "I", "R"]

    SIR = IBM{String}()
    initialize_states(SIR, initial_states, state_labels)

    @test_throws ArgumentError queue_state_update(SIR, 1, "X")
    queue_state_update(SIR, 1, "R")
    @test findfirst(x -> x == 3, subpart(SIR, :state_update)) == 1

    queue_state_update(SIR, 6, "S")
    @test findfirst(x -> x == 1, subpart(SIR, :state_update)) == 6

    apply_state_updates(SIR)

    new_state = initial_states
    new_state[1] = "R"
    new_state[6] = "S"
    @test subpart(SIR, :state) == indexin(new_state, state_labels)
    @test subpart(SIR, :state_update) == zeros(Int64, length(initial_states))

    queue_state_update(SIR, 1:length(initial_states), "I")
    reset_states(SIR, initial_states)
    @test subpart(SIR, :state) == indexin(initial_states, state_labels)
    @test subpart(SIR, :state_update) == zeros(Int64, length(initial_states))

end


text:
function debug_variable(io::IO, variable_to_index, BT, args...)
    for C in added_constrained_variable_types(BT, args...)
        i = get(variable_to_index, C, nothing)
        if i != nothing
            MOIU.print_with_acronym(io, "  [$i] constrained variables in `$(C[1])` are not supported\n")
        end
    end
end
function _print_unsupported_constraint(io::IO, constraint_to_index, C)
    i = get(constraint_to_index, C, nothing)
    if i != nothing
        MOIU.print_with_acronym(io, "  ($i) `$(C[1])`-in-`$(C[2])` constraints are not supported\n")
    end
end
function debug_constraint(io::IO, constraint_to_index, BT, args...)
    for C in added_constraint_types(BT, args...)
        _print_unsupported_constraint(io, constraint_to_index, C)
    end
end
function debug_objective(io::IO, objective_to_index, BT, args...)
    F = set_objective_function_type(BT, args...)
    i = get(objective_to_index, (F,), nothing)
    if i != nothing
        MOIU.print_with_acronym(io, "  |$i| objective function of type `$F` is not supported\n")
    end
end
function debug(b::LazyBridgeOptimizer, variable_to_index::OrderedDict,
               constraint_to_index::OrderedDict, objective_to_index::OrderedDict;
               io::IO = Base.stdout)
    for (v, i) in variable_to_index
        S = v[1]
        MOIU.print_with_acronym(io, "[$i] constraint variables in `$S` are not supported because")
        no_bridge = true
        for BT in b.variable_bridge_types
            if Variable.supports_constrained_variable(BT, S)
                if no_bridge
                    println(io, ":")
                    no_bridge = false
                end
                MOIU.print_with_acronym(io, "  Cannot use `$(Variable.concrete_bridge_type(BT, S))` because:\n")
                debug_variable(io, variable_to_index, BT, S)
                debug_constraint(io, constraint_to_index, BT, S)
            end
        end
        if no_bridge
            println(io, " no added bridge supports bridging it.")
        end
        F = MOIU.variable_function_type(S)
        println(io, "  Cannot add free variables and then constrain them because:")
        _print_unsupported_constraint(io, constraint_to_index, (F, S))
    end
    for (c, i) in constraint_to_index
        F, S = c
        MOIU.print_with_acronym(io, "($i) `$F`-in-`$S` constraints are not supported because")
        no_bridge = true
        for BT in b.constraint_bridge_types
            if MOI.supports_constraint(BT, F, S)
                if no_bridge
                    println(io, ":")
                    no_bridge = false
                end
                MOIU.print_with_acronym(io, "  Cannot use `$(Constraint.concrete_bridge_type(BT, F, S))` because:\n")
                debug_variable(io, variable_to_index, BT, F, S)
                debug_constraint(io, constraint_to_index, BT, F, S)
            end
        end
        if no_bridge
            println(io, " no added bridge supports bridging it.")
        end
    end
    for (o, i) in objective_to_index
        F = o[1]
        MOIU.print_with_acronym(io, "|$i| objective function of type `$F` is not supported because")
        no_bridge = true
        for BT in b.objective_bridge_types
            if Objective.supports_objective_function(BT, F)
                if no_bridge
                    println(io, ":")
                    no_bridge = false
                end
                MOIU.print_with_acronym(io, "  Cannot use `$(Objective.concrete_bridge_type(BT, F))` because:\n")
                debug_variable(io, variable_to_index, BT, F)
                debug_constraint(io, constraint_to_index, BT, F)
                debug_objective(io, objective_to_index, BT, F)
            end
        end
        if no_bridge
            println(io, " no added bridge supports bridging it.")
        end
    end
end
function debug(b::LazyBridgeOptimizer, variables::OrderedSet, constraints::OrderedSet, objectives::OrderedSet; kws...)
    variable_to_index = OrderedDict{Tuple{DataType}, Int}()
    for (i, key) in enumerate(variables)
        variable_to_index[key] = i
    end
    constraint_to_index = OrderedDict{Tuple{DataType, DataType}, Int}()
    for (i, key) in enumerate(constraints)
        constraint_to_index[key] = i
    end
    objective_to_index = OrderedDict{Tuple{DataType}, Int}()
    for (i, key) in enumerate(objectives)
        objective_to_index[key] = i
    end
    debug(b, variable_to_index, constraint_to_index, objective_to_index; kws...)
end
function debug(b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
               S::Type{<:MOI.AbstractSet}; io::IO = Base.stdout)
    if MOI.supports_constraint(b, F, S)
        if F == MOIU.variable_function_type(S)
            # This may be thanks to a variable bridge so `F`-in-`S` constraints
            # are maybe not supported but constraint variables in `S` are
            # definitely supported.
            MOIU.print_with_acronym(io, "Constrained variables in `$S` are supported.\n")
        else
            MOIU.print_with_acronym(io, "`$F`-in-`$S` constraints are supported.\n")
        end
    else
        message = " are not supported and cannot be bridged into supported" *
            " constrained variables and constraints. See details below:"
        if F == MOIU.variable_function_type(S)
            MOIU.print_with_acronym(io, "Constrained variables in `$S`")
            println(io, message)
            required_tuples = required(b, (S,))
            fill_required!(required_tuples..., b, F, S)
            debug(b, required_tuples...; io = io)
        else
            MOIU.print_with_acronym(io, "`$F`-in-`$S` constraints")
            println(io, message)
            debug(b, required(b, (F, S))...; io = io)
        end
    end
end

"""
    debug_supports_constraint(
        b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
        S::Type{<:MOI.AbstractSet}; io::IO = Base.stdout)

Prints to `io` explanations for the value of [`MOI.supports_constraint`](@ref)
with the same arguments.
"""
function debug_supports_constraint(
    b::LazyBridgeOptimizer, F::Type{<:MOI.AbstractFunction},
    S::Type{<:MOI.AbstractSet}; kws...)
    debug(b, F, S; kws...)
end
function debug(b::LazyBridgeOptimizer,
               F::Type{<:MOI.AbstractScalarFunction}; io::IO = Base.stdout)
    MOIU.print_with_acronym(io, "Objective function of type `$F` is")
    if supports_bridging_objective_function(b, F)
        MOIU.print_with_acronym(io, " supported.\n")
    else
        MOIU.print_with_acronym(io, " not supported and cannot be bridged" *
            " into a supported objective function by adding only supported" *
            " constrained variables and constraints. See details below:\n")
        debug(b, required(b, (F,))...; io = io)
    end
end

"""
    debug_supports(b::LazyBridgeOptimizer, ::MOI.ObjectiveFunction{F}; io::IO = Base.stdout) where F

Prints to `io` explanations for the value of [`MOI.supports`](@ref) with the
same arguments.
"""
function debug_supports(
    b::LazyBridgeOptimizer, ::MOI.ObjectiveFunction{F}; kws...) where F
    debug(b, F; kws...)
end


text:
using Reachability.ReachSets: tune_δ

#===== Projectile model =====
We test the line plot!(x->x, x->-24*x+375, 0., 50., line=2., color="red", linestyle=:solid, legend=:none)
=#
A = [0. 0.5 0. 0. ; 0. 0. 0. 0. ; 0. 0. 0. 0.7 ; 0. 0. 0. 0.]
X0 = Singleton([0.,5.,100.,0])
U = ConstantInput(Singleton([0.,0.,0.,-9.81]))
S = IVP(CLCDS(A, Matrix(1.0I, 4, 4), nothing, U), X0)
time_horizon = 20.0
prec = 1e-4
initial_δ = 0.5
property = is_contained_in(LinearConstraint([24., 0., 1, 0], 375.))
algorithm(N, δ) =
    solve(S,
          Options(:mode => "check", :plot_vars => [1, 3], :T => time_horizon,
                  :property => property),
          op=BFFPSV18(:δ => δ, :vars => [1, 3], :partition=>[1:2, 3:4])
         ).satisfied

tune_δ(algorithm, time_horizon, prec, initial_δ)


text:
@testset "Mismatches" begin
    function test_mismatches(a, b)
        count = 0
        for (x, y) in zip(a, b)
            count += x != y
        end
        @test mismatches(a, b) === mismatches(b, a) === count
    end

    for len in 1:64, _ in 1:10
        a = random_dna_kmer(len)
        b = random_dna_kmer(len)
        test_mismatches(DNAKmer(a), DNAKmer(b))
        test_mismatches(Kmer{DNAAlphabet{4}}(a), Kmer{DNAAlphabet{4}}(b))

        a = random_rna_kmer(len)
        b = random_rna_kmer(len)
        test_mismatches(RNAKmer(a), RNAKmer(b))
        test_mismatches(Kmer{RNAAlphabet{4}}(a), Kmer{RNAAlphabet{4}}(b))

        a = AAKmer(random_aa(len))
        b = AAKmer(random_aa(len))
        test_mismatches(a, b)
    end
end

@testset "Matches" begin
    function test_matches(a, b)
        count = 0
        for (x, y) in zip(a, b)
            count += x == y
        end
        @test matches(a, b) === matches(b, a) === count
    end

    for len in 1:64, _ in 1:10
        a = random_dna_kmer(len)
        b = random_dna_kmer(len)
        test_matches(DNAKmer(a), DNAKmer(b))
        test_matches(Kmer{DNAAlphabet{4}}(a), Kmer{DNAAlphabet{4}}(b))

        a = random_rna_kmer(len)
        b = random_rna_kmer(len)
        test_matches(RNAKmer(a), RNAKmer(b))
        test_matches(Kmer{RNAAlphabet{4}}(a), Kmer{RNAAlphabet{4}}(b))

        a = AAKmer(random_aa(len))
        b = AAKmer(random_aa(len))
        test_matches(a, b)
    end
end


text:
using Base.Test

using Rifraf

@testset "RifrafSequence" begin
    @testset "empty" begin
        bandwidth = 10
        scores = Scores(-1., -2., -3., -4., -5.)
        rseq = RifrafSequence(DNASeq(), LogProb[], bandwidth, scores)
        @test length(rseq) == 0
        @test length(RifrafSequence()) == 0
    end

    @testset "scores" begin
        seq = DNASeq("ACGT")
        error_log_p = LogProb[-1., -2., -3., -4.]
        bandwidth = 10
        scores = Scores(-1., -2., -3., -4., -5.)
        rseq = RifrafSequence(seq, error_log_p, bandwidth, scores)

        @test rseq.match_scores == log10.(1.0 - exp10.(error_log_p))
        @test rseq.mismatch_scores == error_log_p + scores.mismatch
        @test rseq.ins_scores == error_log_p + scores.insertion
        @test rseq.del_scores == LogProb[-1., -1., -2., -3., -4.] + scores.deletion
        @test rseq.codon_ins_scores == LogProb[-1., -2.] + scores.codon_insertion
        @test rseq.codon_del_scores == LogProb[-1., -1., -2., -3., -4.] + scores.codon_deletion
    end

    @testset "phred" begin
        seq = DNASeq("ACGT")
        phreds = Int8[3, 50, 10, 70]
        bandwidth = 10
        scores = Scores(-1., -2., -3., -4., -5.)

        # error log p
        rseq = RifrafSequence(seq, phreds, bandwidth, scores)
        # TODO: test
    end

    @testset "update scores" begin
        seq = DNASeq("ACGT")
        error_log_p = LogProb[-1., -2., -3., -4.]
        bandwidth = 10
        scores = Scores(-1., -2., -3., -4., -5.)
        rseq = RifrafSequence(seq, error_log_p, bandwidth, scores)

        new_scores = Scores(-1., -1., -1., -1., -1.)
        new_rseq = RifrafSequence(rseq, new_scores)
        @test new_rseq.ins_scores == new_rseq.mismatch_scores
    end
end