text:
module Expl

using DBInterface
using SQLite
using Dates
using TimeZones
using Unicode
using StringEncodings
using ArgParse
using HTTP

using ...Klio
using ..Mattermost

const MAX_UTF16_LENGTH_ITEM = 50
const MAX_UTF16_LENGTH_EXPL = 500
const MAX_EXPL_COUNT = 20
const MAX_FIND_EXPL_COUNT = 20
const MAX_FIND_ITEM_COUNT = 100
const MAX_TOPEXPL = 100

db_initialized = false

# get the SQLite.DB, creating/updating it if required
function init_db()
    db = SQLite.DB(Klio.settings.expl_sqlite_file)

    # perform idempotent (!) database initialization once per execution
    global db_initialized
    if !db_initialized
        # id must be AUTOINCREMENT because monotonicity is required for some queries
        # nick is NULL for some old entries
        # datetime (unix timestamp) is NULL for some old entries
        # item_norm is a normalized variant of item
        SQLite.execute(db, """
            CREATE TABLE IF NOT EXISTS t_expl (
                id INTEGER PRIMARY KEY AUTOINCREMENT,
                nick TEXT,
                item TEXT NOT NULL,
                item_norm TEXT NOT NULL,
                expl TEXT NOT NULL,
                datetime INTEGER,
                enabled INTEGER NOT NULL
            )
            """)

        SQLite.createindex!(db, "t_expl", "idx_expl_item_norm", "item_norm", unique = false, ifnotexists = true)

        db_initialized = true
    end

    # this should happen automatically, according to the SQLite.jl docs
    SQLite.register(db, SQLite.regexp, nargs = 2, name = "regexp")

    return db
end

# common query to determine the different indexes
# returns all entries (including disabled ones) unordered
const QUERY_BY_ITEM_NORM = """
	SELECT rowid rowid, id, nick, item, item_norm, expl, datetime, enabled,
		   CASE WHEN enabled <> 0 THEN ROW_NUMBER() OVER (PARTITION BY enabled <> 0 ORDER BY id) END normal_index,
		   ROW_NUMBER() OVER (ORDER BY id) permanent_index,
           CASE WHEN enabled <> 0 THEN ROW_NUMBER() OVER (PARTITION BY enabled <> 0 ORDER BY id DESC) END tail_index
	FROM t_expl WHERE item_norm = :item_norm
"""

const QUERY_BY_ITEM_NORM_PARAM = :item_norm

# normalize a string (expl item) for easy searchability
item_normalize(item) = Unicode.normalize(item,
    compat = true,
    casefold = true,
    stripignore = true,
    stripcc = true,
    stable = true)

# number of 16-bit words in the UTF-16 encoding of the given string
# string(s) is required because StringEncodings doesn't support SubString
utf16_length(s) = length(encode(string(s), enc"UTF-16BE")) >> 1

abstract type ExplIndex end

struct NormalExplIndex <: ExplIndex
    index::Int64
end

struct TailExplIndex <: ExplIndex
    index::Int64
end

struct PermanentExplIndex <: ExplIndex
    index::Int64
end

# string representation, used by join and string interpolation
Base.print(io::IO, x::NormalExplIndex) = print(io, string(x.index))
Base.print(io::IO, x::PermanentExplIndex) = print(io, 'p', string(x.index))
Base.print(io::IO, x::TailExplIndex) = print(io, '-', string(x.index))

abstract type ExplIndexSelector end

struct SingleExplIndexSelector <: ExplIndexSelector
    index::ExplIndex
end

struct RangeExplIndexSelector <: ExplIndexSelector
    start::ExplIndex
    stop::ExplIndex
end

struct AllExplIndexSelector <: ExplIndexSelector
end

function Base.tryparse(t::Type{NormalExplIndex}, s::AbstractString)
    index = tryparse(fieldtype(t, :index), s)
    if isnothing(index) || index <= 0
        return nothing
    end
    return NormalExplIndex(index)
end

function Base.tryparse(t::Type{TailExplIndex}, s::AbstractString)
    if !startswith(s, '-')
        return nothing
    end
    index = tryparse(fieldtype(t, :index), s[2:end])
    if isnothing(index) || index <= 0
        return nothing
    end
    return TailExplIndex(index)
end

function Base.tryparse(t::Type{PermanentExplIndex}, s::AbstractString)
    if !startswith(s, 'p')
        return nothing
    end
    index = tryparse(fieldtype(t, :index), s[2:end])
    if isnothing(index) || index <= 0
        return nothing
    end
    return PermanentExplIndex(index)
end

function Base.tryparse(::Type{ExplIndex}, s::AbstractString)
    for t in [NormalExplIndex, TailExplIndex, PermanentExplIndex]
        i = tryparse(t, s)
        if !isnothing(i)
            return i
        end
    end
    return nothing
end

function Base.tryparse(::Type{ExplIndexSelector}, s::AbstractString)
    p = split(s, ':', limit = 3)
    if length(p) > 2
        return nothing
    end

    start = tryparse(ExplIndex, p[1])
    if isnothing(start)
        return nothing
    end

    if length(p) == 1
        return SingleExplIndexSelector(start)
    end

    stop = tryparse(ExplIndex, p[2])
    if isnothing(stop)
        return nothing
    end

    return RangeExplIndexSelector(start, stop)
end

struct ExplEntry
    rowid::Int64
    item::AbstractString
    indexes::Vector{ExplIndex}
    text::AbstractString
end

# string representation, used by join and string interpolation
Base.print(io::IO, e::ExplEntry) = print(io, "$(e.item)[" * join(e.indexes, '/') * "]: $(e.text)")

struct ItemEntry
    item::AbstractString
    count::Int64
end

# string representation, used by join and string interpolation
Base.print(io::IO, e::ItemEntry) = print(io, "$(e.item)($(e.count))")

# unique ExplIndex subtypes used by an ExplIndexSelector
indextypes(s::SingleExplIndexSelector) = [typeof(s.index)]
indextypes(s::RangeExplIndexSelector) = unique([typeof(s.start), typeof(s.stop)])
indextypes(s::AllExplIndexSelector) = [NormalExplIndex]

sqlify(i::NormalExplIndex, op = "=") = "$(i.index) $op normal_index"
sqlify(i::PermanentExplIndex, op = "=") = "$(i.index) $op permanent_index"
sqlify(i::TailExplIndex, op = "=") = "-$(i.index) $op -tail_index"

sqlify(s::SingleExplIndexSelector) = sqlify(s.index)
sqlify(s::RangeExplIndexSelector) = '(' * sqlify(s.start, "<=") * ") AND (" * sqlify(s.stop, ">=") * ')'
sqlify(s::AllExplIndexSelector) = "1 = 1"

sqlify(ss::Vector{<:ExplIndexSelector}) = '(' * join(sqlify.(ss), ") OR (") * ')'

function convert_expl_row(nt, index_types)
    text = replace(nt.:expl, r"[[:space:]]" => " ")

    metadata = []
    if !ismissing(nt.:nick)
        push!(metadata, nt.:nick)
    end
    if !ismissing(nt.:datetime)
        datetime = Dates.format(ZonedDateTime(Dates.epochms2datetime(nt.:datetime), Klio.settings.expl_time_zone, from_utc = true), Klio.settings.expl_datetime_format)
        push!(metadata, datetime)
    end
    if !isempty(metadata)
        text = "$text (" * join(metadata, ", ") * ')'
    end

    indexes = []
    if NormalExplIndex in index_types
        push!(indexes, NormalExplIndex(nt.:normal_index))
    end
    if PermanentExplIndex in index_types
        push!(indexes, PermanentExplIndex(nt.:permanent_index))
    end
    if TailExplIndex in index_types
        push!(indexes, TailExplIndex(nt.:tail_index))
    end

    return ExplEntry(nt.:rowid, nt.:item, indexes, text)
end

function convert_item_row(nt)
    return ItemEntry(nt.:item, nt.:cnt)
end

function add(req)
    parts = split(rstrip(req.text), limit = 3)
    if length(parts) !== 3
        return OutgoingWebhookResponse("Syntax: $(parts[1]) <Begriff> <Erklärung>")
    end
    _, item, expl = parts

    if utf16_length(item) > MAX_UTF16_LENGTH_ITEM
        return OutgoingWebhookResponse("Tut mir leid, der Begriff ist leider zu lang.")
    end
    if utf16_length(expl) > MAX_UTF16_LENGTH_EXPL
        return OutgoingWebhookResponse("Tut mir leid, die Erklärung ist leider zu lang.")
    end

    item_norm = item_normalize(item)

    db = init_db()

    DBInterface.execute(db, "INSERT INTO t_expl(nick, item, item_norm, expl, datetime, enabled) VALUES (:nick, :item, :item_norm, :expl, :datetime, :enabled)",
        (; Dict(
            :nick => req.user_name,
            :item => item,
            :item_norm => item_norm,
            :expl => expl,
            :datetime => Dates.datetime2epochms(Dates.now(Dates.UTC)),
            :enabled => 1
        )...))

    entries = []
    for nt in DBInterface.execute(db, "SELECT * FROM ($QUERY_BY_ITEM_NORM) WHERE enabled <> 0 AND rowid = last_insert_rowid()",
                        (; Dict(QUERY_BY_ITEM_NORM_PARAM => item_norm)...))
        push!(entries, convert_expl_row(nt, [NormalExplIndex, PermanentExplIndex]))
    end

    count = length(entries)
    if count == 1
        text = "Ich habe den folgenden neuen Eintrag hinzugefügt:\n```\n$(entries[1])\n```"
    elseif count == 0
        text = "Ich konnte den eben hinzugefügten Eintrag leider nicht mehr finden."
    else
        error("more than one entry added by request \"$(req.text)\"")
    end

    return OutgoingWebhookResponse(text)
end

function expl(req)
    parts = split(req.text)
    selectors = tryparse.(ExplIndexSelector, parts[3:end])
    if length(parts) < 2 || any(selectors .== nothing)
        return OutgoingWebhookResponse("Syntax: $(parts[1]) <Begriff> { <Index> | <VonIndex>:<BisIndex> }")
    end

    # default selector (all)
    if isempty(selectors)
        selectors = [AllExplIndexSelector()]
    end

    # convert selectors to SQL
    selectors_sql = sqlify(selectors)

    # determine index types used in selectors
    index_types = union(map(indextypes, selectors)...)

    item = parts[2]
    item_norm = item_normalize(item)

    db = init_db()

    entries = []
    for nt in DBInterface.execute(db, "SELECT * FROM ($QUERY_BY_ITEM_NORM) WHERE enabled <> 0 AND ($selectors_sql) ORDER BY id",
                            (; Dict(QUERY_BY_ITEM_NORM_PARAM => item_norm)...))
        push!(entries, convert_expl_row(nt, index_types))
    end

    count = length(entries)
    if count == 0
        text = "Ich habe leider keinen Eintrag gefunden."
    else
        if count == 1
            text = "Ich habe den folgenden Eintrag gefunden:"
        elseif count <= MAX_EXPL_COUNT
            text = "Ich habe die folgenden $count Einträge gefunden:"
        else
            text = "Ich habe $count Einträge gefunden (https://expl.klio.s2000.at/$(HTTP.URIs.escapeuri(item))), das sind die letzten $MAX_EXPL_COUNT:"
            entries = entries[end-MAX_EXPL_COUNT+1:end]
        end

        text = "$text\n```\n" * join(entries, '\n') * "\n```"
    end

    return OutgoingWebhookResponse(text)
end

function www_expl(req)
    # default selector (all)
    selectors = [AllExplIndexSelector()]

    # convert selectors to SQL
    selectors_sql = sqlify(selectors)

    # determine index types used in selectors
    index_types = union(map(indextypes, selectors)...)

    item = replace(req.target, "/expl/" => "")
    item = HTTP.URIs.unescapeuri(item)
    item_norm = item_normalize(item)

    db = init_db()

    entries = []
    for nt in DBInterface.execute(db, "SELECT * FROM ($QUERY_BY_ITEM_NORM) WHERE enabled <> 0 AND ($selectors_sql) ORDER BY id",
                            (; Dict(QUERY_BY_ITEM_NORM_PARAM => item_norm)...))
        push!(entries, convert_expl_row(nt, index_types))
    end

    count = length(entries)
    if count == 0
        text = "Ich habe leider keinen Eintrag gefunden."
    else
        if count == 1
            text = "Ich habe den folgenden Eintrag gefunden:"
        else
            text = "Ich habe die folgenden $count Einträge gefunden:"
        end

        text = "$text\n" * join(entries, '\n')
    end

    text = HTTP.Strings.escapehtml(text)
    text = "<!doctype html><meta charset=\"utf-8\">
$text
"

    return HTTP.Response(text)
end

function del(req)
    parts = split(req.text)
    if length(parts) != 3
        return OutgoingWebhookResponse("Syntax: $(parts[1]) <Begriff> <Index>")
    end

    index = tryparse(ExplIndex, parts[3])
    if isnothing(index)
        return OutgoingWebhookResponse("Syntax: $(parts[1]) <Begriff> <Index>")
    end

    # convert index to SQL
    index_sql = sqlify(index)

    # determine index type used
    index_types = [typeof(index)]

    item = parts[2]
    item_norm = item_normalize(item)

    db = init_db()

    entries = []
    for nt in DBInterface.execute(db, "SELECT * FROM ($QUERY_BY_ITEM_NORM) WHERE enabled <> 0 AND ($index_sql)",
                        (; Dict(QUERY_BY_ITEM_NORM_PARAM => item_norm)...))
        push!(entries, convert_expl_row(nt, index_types))
    end

    changes = 0
    if length(entries) == 1
        DBInterface.execute(db, "UPDATE t_expl SET enabled = 0 WHERE enabled <> 0 AND rowid = :rowid",
                        (; Dict(:rowid => entries[1].rowid)...))

        for nt in DBInterface.execute(db, "SELECT changes() changes")
            changes = nt.:changes
        end
    end

    if changes == 0
        text = "Ich habe leider keinen Eintrag zum Löschen gefunden."
    elseif changes == 1
        text = "Ich habe den folgenden Eintrag gelöscht:\n```\n$(entries[1])\n```"
    else
        error("more than one entry matched request \"$(req.text)\"")
    end

    return OutgoingWebhookResponse(text)
end

# cannot use QUERY_BY_ITEM_NORM, because the indexes need to be partitioned by item_norm
const QUERY_BY_EXPL_REGEXP = """
    SELECT t.rowid rowid, t.id, t.nick, t.item, t.expl, t.datetime,
           (SELECT COUNT(1) FROM t_expl WHERE item_norm = t.item_norm AND enabled <> 0 AND id <= t.id) normal_index,
           (SELECT COUNT(1) FROM t_expl WHERE item_norm = t.item_norm AND id <= t.id) permanent_index,
           (SELECT COUNT(1) FROM t_expl WHERE item_norm = t.item_norm AND enabled <> 0 AND id >= t.id) tail_index
    FROM t_expl t WHERE t.expl REGEXP :regexp AND t.enabled <> 0 ORDER BY t.id
"""

# use the last item to represent a matching item_norm
const QUERY_BY_ITEM_REGEXP = """
    SELECT (SELECT item FROM t_expl WHERE item_norm = t.item_norm ORDER BY id DESC LIMIT 1) item,
           COUNT(1) cnt
    FROM t_expl t WHERE t.item_norm REGEXP :regexp AND t.enabled <> 0 GROUP BY t.item_norm ORDER BY MAX(t.id)
"""

const QUERY_REGEXP_PARAM = :regexp

const FIND_ARG_PARSE_SETTINGS = ArgParseSettings(exc_handler = (s, err) -> rethrow(err))
@add_arg_table! FIND_ARG_PARSE_SETTINGS begin
    "--keys", "-k"
        nargs = 0
        help = "find items whose unicode-normalized item key matches the unicode-normalized regex"
    "regex"
        help = "the regular expression to use"
        required = true
end

function find(req)
    parts = split(req.text)

    local args
    try
        args = parse_args(parts[2:end], FIND_ARG_PARSE_SETTINGS)
    catch err
        if isa(err, ArgParseError)
            return OutgoingWebhookResponse("Syntax: $(parts[1]) [ -k | --keys ] <Suchbegriff>")
        end
        rethrow(err)
    end

    if args["keys"]
        query = QUERY_BY_ITEM_REGEXP
        query_params = Dict(QUERY_REGEXP_PARAM => item_normalize(args["regex"]))
        conv_func = convert_item_row
        max_count = MAX_FIND_ITEM_COUNT
        sep = ", "
    else
        query = QUERY_BY_EXPL_REGEXP
        query_params = Dict(QUERY_REGEXP_PARAM => args["regex"])
        conv_func = nt -> convert_expl_row(nt, [NormalExplIndex, PermanentExplIndex])
        max_count = MAX_FIND_EXPL_COUNT
        sep = "\n"
    end

    db = init_db()

    entries = []
    for nt in DBInterface.execute(db, query, (; query_params...))
        push!(entries, conv_func(nt))
    end

    count = length(entries)
    if count == 0
        text = "Ich habe leider keinen Eintrag gefunden."
    else
        if count == 1
            text = "Ich habe den folgenden Eintrag gefunden:"
        elseif count <= max_count
            text = "Ich habe die folgenden $count Einträge gefunden:"
        elseif args["keys"]
            text = "Ich habe $count Einträge gefunden, das sind die letzten $max_count:"
            entries = entries[end-max_count+1:end]
        else
            text = "Ich habe $count Einträge gefunden (https://find.klio.s2000.at/$(HTTP.URIs.escapeuri(args["regex"]))), das sind die letzten $max_count:"
            entries = entries[end-max_count+1:end]
        end

        text = "$text\n```\n" * join(entries, sep) * "\n```"
    end

    return OutgoingWebhookResponse(text)
end

function www_find(req)
    item = replace(req.target, "/find/" => "")
    item = HTTP.URIs.unescapeuri(item)

    query = QUERY_BY_EXPL_REGEXP
    query_params = Dict(QUERY_REGEXP_PARAM => item)
    conv_func = nt -> convert_expl_row(nt, [NormalExplIndex, PermanentExplIndex])
    sep = "\n"

    db = init_db()

    entries = []
    for nt in DBInterface.execute(db, query, (; query_params...))
        push!(entries, conv_func(nt))
    end

    count = length(entries)
    if count == 0
        text = "Ich habe leider keinen Eintrag gefunden."
    else
        if count == 1
            text = "Ich habe den folgenden Eintrag gefunden:"
        else
            text = "Ich habe die folgenden $count Einträge gefunden:"
        end

        text = "$text\n" * join(entries, sep)
    end

    text = HTTP.Strings.escapehtml(text)
    text = "<!doctype html><meta charset=\"utf-8\">
$text
"

    return HTTP.Response(text)
end

const QUERY_TOPEXPL = """
    SELECT (SELECT item FROM t_expl WHERE item_norm = t.item_norm ORDER BY id DESC LIMIT 1) item,
           COUNT(1) cnt
    FROM t_expl t WHERE t.enabled <> 0 GROUP BY t.item_norm ORDER BY 2 DESC, MAX(t.id) LIMIT :limit
"""

const QUERY_TOPEXPL_LIMIT_PARAM = :limit

function topexpl(req)
    parts = split(req.text)
    if length(parts) > 1
        return OutgoingWebhookResponse("Syntax: $(parts[1])")
    end

    db = init_db()

    entries = []
    for nt in DBInterface.execute(db, QUERY_TOPEXPL,
                            (; Dict(QUERY_TOPEXPL_LIMIT_PARAM => MAX_TOPEXPL)...))
        push!(entries, convert_item_row(nt))
    end

    count = length(entries)
    if count == 0
        text = "Ich habe leider keinen Eintrag gefunden."
    else
        if count == 1
            text = "Ich habe den folgenden wichtigsten Eintrag gefunden:"
        else
            text = "Ich habe die folgenden $count wichtigsten Einträge gefunden:"
        end

        text = "$text\n```\n" * join(entries, ", ") * "\n```"
    end

    return OutgoingWebhookResponse(text)
end

end


text:
# Day 14: Extended Polymerization
function addinsert!(d::Dict, key, val)
    if haskey(d, key)
        d[key] += val
    else
        d[key] = val
    end
end

function load(::Day{14}, path)
    lines = readlines(path)
    polymer = first.(split(lines[1], ""))

    pairs = Dict{Vector{Char},Int}()
    for i in 1:(length(polymer) - 1)
        pair = polymer[i:(i + 1)]
        addinsert!(pairs, pair, 1)
    end

    rules = Dict(
        map(lines[3:end]) do l
            ins, out = split(l, " -> ")
            return first.(split(ins, "")), first(out)
        end,
    )

    return pairs, rules, polymer[end]
end

# Part 1 & 2 are the same – at least after rewriting part 1
function pair_insert(pairs, rules)
    new_pairs = Dict{Vector{Char},Int}()
    for (pair, count) in pairs
        elem = rules[pair]
        addinsert!(new_pairs, [pair[1], elem], count)
        addinsert!(new_pairs, [elem, pair[2]], count)
    end
    return new_pairs
end

function score_polymer(pairs, last_elem)
    element_count = Dict{Char,Int}()
    for (pair, count) in pairs
        addinsert!(element_count, pair[1], count)
    end
    addinsert!(element_count, last_elem, 1)
    min, max = extrema(values(element_count))
    return max - min
end

function polymerize(pairs, rules, steps)
    for _ in 1:steps
        pairs = pair_insert(pairs, rules)
    end
    return pairs
end

function solve1(::Day{14}, data)
    pairs, rules, last_element = data
    return score_polymer(polymerize(pairs, rules, 10), last_element)
end
function solve2(::Day{14}, data)
    pairs, rules, last_element = data
    return score_polymer(polymerize(pairs, rules, 40), last_element)
end

testresult1(::Day{14}) = 1588
testresult2(::Day{14}) = 2188189693529

ext: jl


text:
export
    OrbitSystemNode,
    set!

using Printf

# import ..Nodes.draw

mutable struct OrbitSystemNode <: Ranger.AbstractNode
    base::Nodes.NodeData
    transform::Nodes.TransformProperties{Float64}

    # Collection of nodes.
    children::Array{Ranger.AbstractNode,1}

    color::Rendering.Palette

    polygon::Geometry.Polygon

    angular_motion::Animation.AngularMotion{Float64}
    anchor_motion::Animation.AngularMotion{Float64}
    tri_motion::Animation.AngularMotion{Float64}

    anchor::Custom.AnchorNode
    triangle::Custom.OutlinedTriangle

    detection::Nodes.Detection
    inside::Bool

    aabb::Geometry.AABB{Float64}

    function OrbitSystemNode(world::Ranger.World, name::String, parent::Ranger.AbstractNode)
        o = new()

        o.base = Nodes.NodeData(name, parent, world)
        o.transform = Nodes.TransformProperties{Float64}()
        o.children = Array{Ranger.AbstractNode,1}[]
        o.color = Rendering.White()
        o.polygon = Geometry.Polygon{Float64}()
        o.angular_motion = Animation.AngularMotion{Float64}()
        o.anchor_motion = Animation.AngularMotion{Float64}()
        o.tri_motion = Animation.AngularMotion{Float64}()

        o.aabb = Geometry.AABB{Float64}()
        o.detection = Nodes.Detection(Rendering.Lime(), Rendering.Red())
        o.inside = false

        o
    end
end

function build(node::OrbitSystemNode, world::Ranger.World)
    Geometry.add_vertex!(node.polygon, 0.0, 0.0)
    Geometry.add_vertex!(node.polygon, 0.0, 0.0)
    Geometry.add_vertex!(node.polygon, 0.0, 0.0)
    Geometry.add_vertex!(node.polygon, 0.0, 0.0)

    Geometry.build!(node.polygon)

    # amgle is measured in angular-velocity or "degrees/second"
    node.angular_motion.rate = -45.0    # degrees/second

    # The anchor will rotate and we want that rotation to propagate to the
    # children.
    # Another way to say this is:
    # The anchor has a position, rotation and scale. We want to "filter-out" the
    # from any children as the children will determine their own scale.
    node.anchor = Custom.AnchorNode(world, "AnchorNode", node)
    Nodes.set_scale!(node.anchor, 20.0)
    node.anchor.color = RangerGame.lightblue
    node.anchor_motion.rate = 45.0

    # The child filter of the anchor node above.
    anchor_filter = Filters.TransformFilter(world, "AnchorTransformFilter", node.anchor)
    # Change the default rotation flag from "true" to "false".
    anchor_filter.exclude_rotation = false
    push!(node.anchor.children, anchor_filter);

    # Now we add children of the filter.
    node.triangle = Custom.OutlinedTriangle(world, "YellowTriangle", anchor_filter)
    Nodes.set_scale!(node.triangle, 50.0)
    Nodes.set_position!(node.triangle, 200.0, 0.0)
    # Nodes.set_scale!(node.triangle, 2.0)
    # Nodes.set_position!(node.triangle, 3.0, 0.0)
    node.triangle.color = RangerGame.yellow
    push!(anchor_filter.children, node.triangle);
    node.tri_motion.rate = -90.0    # degrees/second
end

# --------------------------------------------------------
# Timing
# --------------------------------------------------------
function Nodes.update(node::OrbitSystemNode, dt::Float64)
    # println("OrbitSystemNode::update : ", layer)
    Animation.update!(node.angular_motion, dt);
    Animation.update!(node.anchor_motion, dt);
    Animation.update!(node.tri_motion, dt);

    Nodes.update!(node.detection, node)

    Nodes.update(node.triangle, dt)
end

function Nodes.interpolate(node::OrbitSystemNode, interpolation::Float64)
    value = Animation.interpolate!(node.angular_motion, interpolation)
    Nodes.set_rotation_in_degrees!(node, value)

    value = Animation.interpolate!(node.anchor_motion, interpolation)
    Nodes.set_rotation_in_degrees!(node.anchor, value)

    value = Animation.interpolate!(node.tri_motion, interpolation)
    Nodes.set_rotation_in_degrees!(node.triangle, value)
end

# --------------------------------------------------------
# Life cycle events
# --------------------------------------------------------
function Nodes.enter_node(node::OrbitSystemNode, man::Nodes.NodeManager)
    println("enter ", node);
    node.inside = Nodes.check!(node.detection, node, node.base.world)
    # Register node as a timing target in order to receive updates
    Nodes.register_target(man, node)
    Nodes.register_event_target(man, node);
end

function Nodes.exit_node(node::OrbitSystemNode, man::Nodes.NodeManager)
    println("exit ", node);
    Nodes.unregister_target(man, node);
    Nodes.unregister_event_target(man, node);
end

function Nodes.draw(node::OrbitSystemNode, context::Rendering.RenderContext)
    if Nodes.is_dirty(node)
        Rendering.transform!(context, node.polygon)
        Nodes.set_dirty!(node, false)
    end

    Rendering.set_draw_color(context, node.color)
    Rendering.render_outlined_polygon(context, node.polygon, Rendering.CLOSED);

    # Nodes.draw(node.detection, context)
    aabb_color = Nodes.highlight_color(node.detection, node.inside)
    Rendering.set_draw_color(context, aabb_color)

    Geometry.set!(node.aabb, Nodes.get_bucket(node))
    Rendering.render_aabb_rectangle(context, node.aabb)
end

function set!(node::OrbitSystemNode, minx::Float64, miny::Float64, maxx::Float64, maxy::Float64)
    Geometry.set!(node.polygon.mesh.vertices[1], minx, miny)
    Geometry.set!(node.polygon.mesh.vertices[2], minx, maxy)
    Geometry.set!(node.polygon.mesh.vertices[3], maxx, maxy)
    Geometry.set!(node.polygon.mesh.vertices[4], maxx, miny)

    Nodes.set_dirty!(node, true)
end

# --------------------------------------------------------
# Events
# --------------------------------------------------------
function Nodes.io_event(node::OrbitSystemNode, event::Events.MouseMotionEvent)
    # println("io_event ", event, ", node: ", node)
    Nodes.set_device_point!(node.detection, Float64(event.x), Float64(event.y))
    node.inside = Nodes.check!(node.detection, node, node.base.world)

    Nodes.io_event(node.triangle, event)
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function Nodes.get_children(node::OrbitSystemNode)
    node.children
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# SimLynx/src/statistics.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

# Thanks to SimLynx.jl for some of this code.

abstract type Stats{T<:Real} end

"""
    AccumulatedStats{T<:Real} <: Stats{T}

    AccumulatedStats()
    AccumulatedStats(value::T)

Collection of statistics to describe the DAG.

# Fields
- `min::T` - the minimum for an array of values
- `max::T` - the maximum for an array of values
- `n::Float64` - the weighted number of an array of values
- `sum::Float64` - the weighted sum for an array of values
- `sum_squares::Float64` - the weighted sum of squares for an array of values
"""
mutable struct Statistics{T<:Real} <: Stats{T}
    min::T
    max::T
    n::Float64
    sum::Float64
    sum_squares::Float64
    Statistics{T}() where {T<:Real} =
        new(zero(T), zero(T), 0.0, 0.0, 0.0)
        Statistics{T}(value::T) where {T<:Real} =
        new(value, value, 0.0, 0.0, 0.0)
end

# Stats Operations
mean(stats::Stats) = stats.sum / stats.n
mean_square(stats::Stats) = stats.sum_squares / stats.n
variance(stats::Stats) = mean_square(stats) - mean(stats)^2
stddev(stats::Stats) = sqrt(variance(stats))

"Overloaded statistics properties for mean, variance, and stddev."
function Base.getproperty(stats::Stats, name::Symbol)
    if name === :mean
        return mean(stats)
    elseif name === :variance
        return variance(stats)
    elseif name === :stddev
        return stddev(stats)
   else # fallback to getfield
       return getfield(stats, name)
   end
end

"""
Print the statistics with optional format, title, and output stream
"""
function print_stats(stats::Stats;
                     format="%f",
                     title="Statistics",
                     output=stdout)
    println(output, title)
    println(output, "     min = $(stats.min)")
    println(output, "     max = $(stats.max)")
    println(output, "       n = $(stats.n)")
    println(output, "    mean = $(stats.mean)")
    println(output, "variance = $(stats.variance)")
    println(output, "  stddev = $(stats.stddev)")
    flush(output)
end


text:
module FailedInterface

import DBAPI
using FactCheck
using DataStructures

type BadInterface <: DBAPI.DatabaseInterface end
type BadConnection <: DBAPI.DatabaseConnection{BadInterface} end
type BadCursor <: DBAPI.DatabaseCursor{BadInterface} end
type BadFixedLengthCursor <: DBAPI.FixedLengthDatabaseCursor{BadInterface} end
type BadQuery <: DBAPI.DatabaseQuery end

facts("Failed Interface") do
    # data structures for testing
    empty_data_structures = (
        Array{Any}[Array{Any}(0)],
        Array{Any}[Array{Any}(0, 0)],
        Dict{Any, Any}[Dict{Any, Any}()],
        PriorityQueue[PriorityQueue()],
        Dict{Any,Array{Any}}(1=>Array{Any}(0)),
        Dict{Any,Array{Any}}(1=>Array{Any}(0, 0)),
        Array{Nullable{Any}}[Array{Nullable{Any}}(0)],
        Array{Nullable{Any}}[Array{Nullable{Any}}(0, 0)],
        Dict{Any, Nullable{Any}}[Dict{Any, Nullable{Any}}()],
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(0)),
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(0, 0)),
    )

    empty_2d_data_structures = (
        Array{Any}(0,0),
        Dict{Any, Any}(),
        PriorityQueue(),
        Array{Nullable{Any}}(0,0),
        Dict{Any, Nullable{Any}}(),
    )

    filled_pq = PriorityQueue()
    filled_pq[1] = 1

    filled_data_structures = (
        Array{Any}[Array{Any}(1)],
        Array{Any}[Array{Any}(1, 1)],
        Dict{Any, Any}[Dict{Any, Any}(1=>5)],
        PriorityQueue[filled_pq],
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(1)),
        Array{Nullable{Any}}[Array{Nullable{Any}}(1)],
        Array{Nullable{Any}}[Array{Nullable{Any}}(1, 1)],
        Dict{Any, Nullable{Any}}[Dict{Any, Nullable{Any}}(1=>5)],
        Dict{Any,Array{Nullable{Any}}}(1=>Array{Nullable{Any}}(1)),
    )

    filled_2d_pq = PriorityQueue()
    filled_2d_pq[1, 1] = 1

    filled_2d_data_structures = (
        Array{Any}(1),
        Array{Any}(1, 1),
        Dict{Any, Any}((1, 1)=>5),
        filled_2d_pq,
        Array{Nullable{Any}}(1),
        Array{Nullable{Any}}(1, 1),
        Dict{Any, Nullable{Any}}((1, 1)=>Nullable{Any}(5)),
    )

    dummy_io = IOBuffer()

    @fact_throws DBAPI.NotImplementedError DBAPI.connect(BadInterface)
    @fact_throws DBAPI.NotImplementedError DBAPI.connect(BadInterface, "foobar")
    @fact_throws DBAPI.NotImplementedError DBAPI.connect(BadInterface, "foobar"; port=2345)
    @fact_throws DBAPI.NotImplementedError DBAPI.connect(identity, BadInterface)

    connection = BadConnection()

    @fact_throws DBAPI.NotImplementedError Base.show(connection)
    @fact_throws DBAPI.NotImplementedError DBAPI.commit(connection)
    @fact_throws DBAPI.NotImplementedError DBAPI.isopen(connection)
    @fact_throws DBAPI.NotImplementedError DBAPI.close(connection)
    @fact_throws DBAPI.NotImplementedError DBAPI.isopen(connection)
    @fact_throws DBAPI.NotImplementedError DBAPI.commit(connection)
    @fact_throws DBAPI.NotSupportedError DBAPI.rollback(connection)
    @fact_throws DBAPI.NotImplementedError DBAPI.cursor(connection)

    # no length
    cursor = BadCursor()
    @fact isa(cursor, DBAPI.DatabaseCursor) --> true
    @fact_throws DBAPI.NotSupportedError length(cursor)

    for i in [12, :twelve]
        for j in [12, :twelve]
            @fact_throws DBAPI.NotSupportedError cursor[i, j]
        end
    end

    @fact_throws DBAPI.NotSupportedError cursor["far", Set([1,2,3])]

    for ds in empty_2d_data_structures
        @fact_throws DBAPI.NotSupportedError DBAPI.fetchintoarray!(ds, cursor)
    end

    for ds in empty_data_structures
        @fact_throws DBAPI.NotSupportedError DBAPI.fetchintoarray!(ds, cursor)
    end

    for ds in empty_data_structures
        @fact_throws DBAPI.NotSupportedError DBAPI.fetchintoarray!(ds, cursor)
    end

    for ds in filled_2d_data_structures
        @fact_throws DBAPI.NotSupportedError DBAPI.fetchintoarray!(ds, cursor)
    end

    for ds in filled_data_structures
        @fact_throws DBAPI.NotSupportedError DBAPI.fetchintorows!(ds, cursor)
    end

    for ds in filled_data_structures
        @fact_throws DBAPI.NotSupportedError DBAPI.fetchintocolumns!(ds, cursor)
    end

    # fixed length
    cursor = BadFixedLengthCursor()
    @fact isa(cursor, DBAPI.DatabaseCursor) --> true
    @fact isa(cursor, DBAPI.FixedLengthDatabaseCursor) --> true

    @fact_throws DBAPI.NotImplementedError DBAPI.connection(cursor)
    @fact_throws DBAPI.NotImplementedError Base.show(cursor)

    @fact_throws DBAPI.NotImplementedError DBAPI.rows(cursor)
    @fact_throws DBAPI.NotSupportedError DBAPI.columns(cursor)

    @fact_throws DBAPI.NotSupportedError DBAPI.execute!(cursor, "foobar", (1, "d"))

    for i in [12, :twelve]
        for j in [12, :twelve]
            @fact_throws DBAPI.NotImplementedError cursor[i, j]
        end
    end

    @fact_throws DBAPI.NotImplementedError cursor["far", Set([1,2,3])]

    for ds in empty_2d_data_structures
        @fact DBAPI.fetchintoarray!(ds, cursor) --> exactly((ds, 0))
    end

    for ds in empty_data_structures
        # in this case there is a row, it's just empty
        @fact DBAPI.fetchintorows!(ds, cursor) --> exactly((ds, 1))
    end

    for ds in empty_data_structures
        @fact DBAPI.fetchintocolumns!(ds, cursor) --> exactly((ds, 0))
    end

    for ds in filled_2d_data_structures
        @fact_throws DBAPI.NotImplementedError DBAPI.fetchintoarray!(ds, cursor)
    end

    for ds in filled_data_structures
        @fact_throws DBAPI.NotImplementedError DBAPI.fetchintorows!(ds, cursor)
    end

    for ds in filled_data_structures
        @fact_throws DBAPI.NotImplementedError DBAPI.fetchintocolumns!(ds, cursor)
    end

    # testing that these methods exist and run
    # @fact_throws does not do that, unfortunately
    Base.showerror(dummy_io, DBAPI.NotImplementedError{BadInterface}())
    Base.showerror(dummy_io, DBAPI.NotSupportedError{BadInterface}())
    Base.showerror(dummy_io, DBAPI.DatabaseQueryError(BadInterface, BadQuery()))

    for empty_ds in empty_data_structures
        @fact_throws DBAPI.NotImplementedError begin
            for ds in DBAPI.DatabaseFetcher(:rows, empty_ds, cursor)
                @fact true --> false  # should never be reached
            end
        end
    end

    for empty_ds in empty_data_structures
        @fact_throws DBAPI.NotImplementedError begin
            for ds in DBAPI.DatabaseFetcher(:columns, empty_ds, cursor)
                @fact true --> false  # should never be reached
            end
        end
    end

    for empty_ds in empty_data_structures
        @fact_throws DBAPI.NotImplementedError begin
            for ds in DBAPI.DatabaseFetcher(:array, empty_ds, cursor)
                @fact true --> false  # should never be reached
            end
        end
    end
end

end


text:
using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

# 1. read in and reshape
cd("C:/Users/jakek/Documents/econ_phd/structural_modeling_course/forked_class/ProblemSets/PS5-ddc/")
include("create_grids.jl")
df = CSV.read("busdataBeta0.csv")
df = @transform(df, bus_id = 1:size(df,1))
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

##### 2. static model.
beta_glm = glm(@formula(Y ~ Odometer+Branded ), df_long, Binomial(), LogitLink())
println(beta_glm)

###### 3.
beta = 0.9
data = CSV.read("busdata.csv")
Y = @select(data, :Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20)
Y = convert(Matrix, Y)
O = @select(df, :Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
O = convert(Matrix, O)
X = data[:, 43:62]
Z = Vector(df[:,:RouteUsage])
X = convert(Matrix, X)
Zst = convert(Vector, data[:,:Zst])
B = convert(Vector, data[:,:Branded])
# create xst

## b
zval,zbin,xval,xbin,xtran = create_grids()

### c.
@views @inbounds function zurcher(theta, Y, X, O,Zst,B,beta,zval,zbin,xval,xbin,xtran)
	T = promote_type(eltype(X),eltype(theta), eltype(Y), eltype(O), eltype(B), eltype(Zst),eltype(beta))
	Fut = repeat(zeros(T,size(xtran,1),2 ),1,1,21)
	for t=reverse(1:20)
		for b=0:1
			for z=1:zbin
				for x=1:xbin
					row = x+(z-1)*xbin
					v1 = theta[1]+theta[2]*xval[x] + theta[3]*b+(xtran[row,:]'*Fut[(z-1)*xbin+1:z*xbin,b+1,t+1])
					v0 = xtran[1+(z-1)*xbin,:]'*Fut[(z-1)*xbin+1:z*xbin,b+1,t+1]
					Fut[row,b+1,t] = beta*log(exp(v1)+exp(v0))
				end
			end
		end
	end
	loglik = 0
	for i=1:size(Y,1)
		for t=1:20
			row0 = 1+(Zst[i]-1)*xbin
			row1 = X[i,t]+(Zst[i]-1)*xbin
			vflow = theta[1]+theta[2]*O[i,t]+ theta[3]*B[i]
			vdiff = vflow + (xtran[row1,:].-xtran[row0,:])'*Fut[row0:row0+xbin-1,B[i]+1,t+1]
			loglik = loglik + (Y[i,t]==1)*vdiff - log(1+exp(vdiff))
		end
	end
	return -loglik
end

starting = vec(coef(beta_glm))
td = TwiceDifferentiable(b -> zurcher(b, Y, X, O, Zst,B,beta,zval,zbin,xval,xbin,xtran), θ_true; autodiff = :forward)
beta_zurcher = optimize(td,θ_true, Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true))
H  = Optim.hessian!(td, beta_zurcher.minimizer)
se= sqrt.(diag(inv(H)))
res = vcat(beta_zurcher.minimizer', se' )
CSV.write("q5.csv", convert(DataFrame,res'))


text:
module JuliaReachDevDocs

# this project is used for documentation purposes only

end # module

max_stars_repo_name: JuliaReach/JuliaReachDevDocs

max_issues_repo_name: JuliaReach/JuliaReachDevDocs

max_forks_repo_name: JuliaReach/JuliaReachDevDocs


text:
include_dependency("core.jl")

#=
# NOTE: This function needs to be adjusted.
"""
Find a longest path on DD.

# Input
- `dd::DecisionDiagram`: The DD.
- `func::Array{Function}`: An array of weight functions (each component represents the separable function for that variable)

# Output
- The longest path (composed of arc labels)
- The length of the longest path
"""
function longest_path(dd::DecisionDiagram, func::Array{Function})

    n = nvars(dd)
    @assert(n == length(func))       #makes sure the size of the fractional point matches the dimension of variables represented by dd

    node_num = nnodes(dd)
    node_obj_vals = Array{Float64}(  node_num)   #stores the maximum objective value computed from the source to each node
    node_obj_paths = Array{Int}(  node_num)   #at each node, stores the index of the previous node of longest path from source to that node
    node_obj_labels = Array{Int}(  node_num)  #stores the label value of the last arc of the longest path from source to each node

    #setting the root objective value equal to zero
    node_obj_vals[root(dd)] = 0

    #computing the longest path from source to a node, by searching through all nodes and their incoming arcs, layer by layer
    for i=2:n+1, node in dd.layers[i]
        temp_v = -Inf
        temp_p = 0
        temp_l = 0
        for parent in inneighbors(dd, node), label in get_arc_labels(dd, parent, node)
            val = node_obj_vals[parent] + func[i-1](label)    #computes the head node value based on the tail and arc label
            if val > temp_v
                temp_v = val
                temp_p = parent
                temp_l = label
            end
        end
        node_obj_vals[node] = temp_v
        node_obj_paths[node] = temp_p
        node_obj_labels[node] = temp_l
    end

    terminal_node = terminal(dd)
    lpv = node_obj_vals[terminal_node]   #the length of the longest path

    lp = Array{Int}(  0)
    p = terminal_node
    #computing the arc labels on the longest path through a backtracking
    for i=n:-1:1
        # push!(lp, get_node_obj_label(dd, p))
        push!(lp, node_obj_labels[p])
        # p = get_node_obj_path(dd, p)
        p = node_obj_paths[p]
    end

    return reverse(lp), lpv
end
=#

"""
Second method for function Longest_Path, where weight functions are linear for all variables.

# Input
- `dd::DecisionDiagram`: The DD.
- `cf::Array{T} where T<:Float64`: An array of coefficients of the variables in the linear objective function.
"""

#=
function longest_path(dd::DecisionDiagram, cf::Array{T} where T<:Float64)
    func = Array{Function}(  length(cf))
    for (i, c) in enumerate(cf)                 #creates an array of separable functions for each variable
        func[i] = (x::Int64) -> c*x
    end
    return longest_path(dd, func)
end
=#

#This version is for DD structure where we store arc layers, and use them directly to compute the longest path
function longest_path(dd::DecisionDiagram, cf::Array{Float64})::Tuple{Array{Float64, 1}, Float64}

    n::Int = nvars(dd)
    @assert(n == length(cf))       #makes sure the size of the fractional point matches the dimension of variables represented by dd

    aux_val::Array{Array{Float64,1}, 1} = [Array{Float64,1}(0) for i=1:n+1]

    #setting the root objective value equal to zero
    push!(aux_val[1], 0)

    #computing the longest path from source to a node, by searching through all nodes and their incoming arcs, layer by layer
    for i::Int=1:n
        aux_val[i+1] = fill(-Inf, get_node_layer_size(dd, i+1))
        for j::Int=1:get_arc_layer_size(dd, i)
            val::Float64 = aux_val[i][get_arc_tail(dd, i, j)] + cf[i]*get_arc_label(dd, i, j)    #computes the head node value based on the tail and arc label
            if val > aux_val[i+1][get_arc_head(dd, i, j)]
                aux_val[i+1][get_arc_head(dd, i, j)] = val
            end
        end
    end

    # Computing the max obj for the terminal node
    lpv::Float64 = aux_val[n+1][1]   #the length of the longest path

    lp::Array{Float64, 1} = Array{Float64, 1}(0)
    p::Int = 1                   #terminal node id at the last layer
    # Computing the arc labels on the longest path through a backtracking
    for i::Int=n+1:-1:2
        for tail in values(inneighbors(dd, i, p)), in_arc::Int in tail
        #for in_arc::Int in inneighbors(dd, i, p)                       # this is for inneighbor that is array of int
            val::Float64 = aux_val[i-1][get_arc_tail(dd, i-1, in_arc)] + cf[i-1]*get_arc_label(dd, i-1, in_arc)    #computes the head node value based on the tail and arc label
            if val == aux_val[i][p]
                push!(lp, get_arc_label(dd, i-1, in_arc))
                p = get_arc_tail(dd, i-1, in_arc)
                break
            end
        end
    end

    return reverse(lp), lpv
end

max_stars_repo_name: ctjandra/ddopt-ecp

max_issues_repo_name: ctjandra/ddopt-ecp

max_forks_repo_name: ctjandra/ddopt-ecp


text:
# Command: Get Page Source
"""
    source(session::Session)::String

The Get Page Source command returns a string serialization of the DOM of the current browsing context active document.
"""
function source(session::Session)::String
    @unpack addr, id = session
    response =
        HTTP.get("$addr/session/$id/source", [("Content-Type" => "application/json; charset=utf-8")])
    @assert response.status == 200
    JSON3.read(response.body).value
end


text:
abstract type Event end

struct StreamStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    encoding::AbstractString
end

struct StreamEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
end

struct DocumentStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    explicit::Bool
    version::Union{Tuple, Nothing}
    tags::Union{Dict{AbstractString, AbstractString}, Nothing}

    function DocumentStartEvent(start_mark::Mark,end_mark::Mark,
                                explicit::Bool, version=nothing,
                                tags=nothing)
        new(start_mark, end_mark, explicit, version, tags)
    end
end

struct DocumentEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
    explicit::Bool
end

struct AliasEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::Union{AbstractString, Nothing}
end

struct ScalarEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::Union{AbstractString, Nothing}
    tag::Union{AbstractString, Nothing}
    implicit::Tuple
    value::AbstractString
    style::Union{Char, Nothing}
end

struct SequenceStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::Union{AbstractString, Nothing}
    tag::Union{AbstractString, Nothing}
    implicit::Bool
    flow_style::Bool
end

struct SequenceEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
end

struct MappingStartEvent <: Event
    start_mark::Mark
    end_mark::Mark
    anchor::Union{AbstractString, Nothing}
    tag::Union{AbstractString, Nothing}
    implicit::Bool
    flow_style::Bool
end

struct MappingEndEvent <: Event
    start_mark::Mark
    end_mark::Mark
end


text:
"""
    getmask(ls::Container{<:Container})

get the mask for batched data.
"""
function getmask(ls::Container{<:Container})
    lens = map(length, ls)
    m = zeros(Float32, maximum(lens), length(lens))

    for (i, l) ∈ enumerate(ls)
        selectdim(selectdim(m, 2, i), 1, 1:length(l)) .= 1
    end
    reshape(m, (1, size(m)...))
end

"""
    getmask(m1::A, m2::A) where A <: Abstract3DTensor

get the mask for the covariance matrix.
"""
getmask(m1::A, m2::A) where A <: Abstract3DTensor = permutedims(m1, [2,1,3]) .* m2

"""
    Embed(size::Int, Vocab::Vocabulary)

The Embedding Layer, `size` is the hidden size. `Vocab` is vocabulary set.
"""
struct Embed{F ,W <: AbstractArray{F}, T}
    Vocab::Vocabulary{T}
    embedding::W
end

@treelike Embed

Embed(size::Int, vocab::Vocabulary) = Embed(vocab, param(randn(Float32, size, length(vocab))))

(e::Embed)(x::AbstractArray{Int}) = e.embedding * onehotarray(length(e.Vocab), x)
(e::Embed{F,W,T})(x::Container{<:Union{T, Container{T}}}) where {F,W,T} = e.embedding * onehot(e, x)
(e::Embed{F})(x, scale) where {F} = e(x) .* convert(F, scale)

Flux.onehot(e::Embed{F, W, T}, x::Container{<:Union{T, Container{T}}}) where {F,W,T} = onehot(e.Vocab, x)
Flux.onehot(v::Vocabulary{T}, x::Container{<:Union{T, Container{T}}}) where T = onehot(v, v(x))
Flux.onehot(e::Embed, x::AbstractArray{Int}) = onehot(e.Vocab, x)
Flux.onehot(v::Vocabulary, x::AbstractArray{Int}) = onehotarray(length(v), x)

Flux.onecold(e::Embed, p) = Flux.onecold(e.Vocab, p)
Flux.onecold(v::Vocabulary{T}, p) where T = decode(v, _onecold(p))

function _onecold(p)
    p = collect(p)
    y = Array{Int}(undef, Base.tail(size(p)))
    for i = 1:length(y)
        ind = Tuple(CartesianIndices(y)[i])
        y[i] = argmax(@view(p[:, ind...]))
    end
    y
end

function Base.show(io::IO, e::Embed)
    print(io, "Embed($(size(e.embedding, 1)), ")
    show(io,  e.Vocab)
    print(io, ")")
end


text:
#TODO: Newick parser

"""
    to_df(root::Node)::DataFrame

This function returns a matrix representation of the tree structure. The matirx
is returned as a DataFrame so that the names of the columns are the names of the
tips in the tree. The entry `df[i,j]` is the length of the edge connecting node
`i` with node `j`.

Returns a Dataframe.

* `root` : root of tree used to create matrix represenation.
"""
function to_df(root::Node)::DataFrame

    post_order_iteration = post_order(root)

    name_list = [i.name for i in post_order_iteration]

    temp_ar = zeros(Float64, (length(post_order_iteration), length(post_order_iteration)))

    for i in post_order_iteration
        if i.nchild != 0
            ind = indexin(i.name, name_list)
            for j in i.child
                ind2 = indexin(j.name, name_list)
                temp_ar[ind[1], ind2[1]] = j.inc_length
            end # end for
        end # end if
    end # end for
    df = convert(DataFrame, temp_ar)
    names!(df, [Symbol(i) for i in name_list])
    return df
end # end function to_df

"""
    from_df(df::DataFrame)::Node

This function takes a DataFrame and turns it into a tree. It assumes a rooted
binary tree is stored in the matrix. No checks are performed.

Returns the root node of the tree.

* `df` : Dataframe used to create a tree.
"""
function from_df(df::DataFrame)::Node

    node_list::Vector{Node} = [Node(String(i) , [0.0], Node[], 0, true, 0.0, "0") for i in names(df)]

    for (col_index, col) in enumerate(eachcol(df))
        for (row_index, entry) in enumerate(col)
            if entry != 0
                # there is a branch connecting these nodes
                node_list[col_index].inc_length = entry
                add_child!(node_list[row_index], node_list[col_index])
            end # end if
        end # end for
    end # end for
    i::Int = 0
    # find the root
    for (ind, n) in enumerate(node_list)
        if n.root == true
            i = ind
            break
        end # end if
    end # for
    node::Node = node_list[i]

    # do some bookeeping here and set the binary representation of the nodes
    set_binary!(node)
    return node
end # function from_df

"""
    newick(node::Node)::String
Creates a newick representation of the tree.

Returns a properly formatted newick String.

* `node` : root node of tree used to create the newick string.
"""
function newick(root::T)::String  where T<:GeneralNode
    # get the newickstring
    newickstring = newick(root, "")

    # Some polishing.
    newickstring = chop(newickstring)
    newickstring = string(newickstring, ";")
    return newickstring
end

"""
    newick(root::T, newickstring::AbstractString) where T<:GeneralNode

Do the newick recursion. It is meant as the internal iterator function.
"""
function newick(root::T, newickstring::AbstractString) where T<:GeneralNode
    if root.nchild != 0
        # internal node
        newickstring = string(newickstring, "(")
        for child in root.children
            newickstring = string(newick(child,newickstring))
        end # for
        newickstring = chop(newickstring)
        return string(newickstring,")", root.name, ":", root.inc_length,",")

    else
        # leave
        return string(newickstring, root.name, ":", root.inc_length, ",")
    end # if
end

#################### Covariance wrapper ####################
function to_covariance(tree::TreeStochastic)::Array{T,2} where T <: Real
    blv::Vector{T} = get_branchlength_vector(tree)
    to_covariance(tree.value, blv)
end # end to_covariance

function to_covariance(tree::N) where {N<:GeneralNode, T <: Real}
    blv = get_branchlength_vector(tree)
    to_covariance(tree, blv)
end # end to_covariance

function to_covariance(tree::TreeStochastic, blv::Vector{T})::Array{T,2} where T<: Real
    to_covariance(tree.value, blv)
end # end to_covariance

"""
    to_covariance_ultra(tree::Node)::Array{T,2} where T<: Real

Get the covariance matrix of the ultrametric version of `tree` with height 1.

Returns an Array of Real numbers.

* `tree` : root of tree used to perform calculation.

"""
function to_covariance_ultra(tree::N)::Array{T,2} where {T<: Real, N <:GeneralNode}
    # scale the branchlength between 0 and 1
    blv = get_branchlength_vector(tree)
    blv ./= tree_height(tree)

    # copy the tree, to maintain the original one
    root = deepcopy(tree)
    set_branchlength_vector!(root, blv)
    # make the tree ultrametric
    force_ultrametric!(root)

    # calculate variance-covariance matrix
    to_covariance(root, blv)
end # end function to_covariance_ultra

#################### To Matrix convertes ####################
"""
    to_distance_matrix(tree::T)::Array{Float64,2} where T <:GeneralNode

Calculate the distance matrix over the set of leaves.

Returns an Array of Floats.

* `tree` : root node of tree used to perform caclulcation.
"""
function to_distance_matrix(tree::T)::Array{Float64,2} where T <:GeneralNode
    leaves::Vector{T} = get_leaves(tree)
    ll = length(leaves)
    distance_mat = zeros(Float64, ll, ll)
    for i in 1:ll
        for j in 1:ll
            if i>j
                d = node_distance(tree, leaves[i], leaves[j])
                distance_mat[i,j] = d
                distance_mat[j,i] = d
            end # if
        end # for
    end #for
    distance_mat
end # function to_distance_matrix

"""
    to_covariance(tree::N, blv::Array{T})::Array{T,2} where {N<:GeneralNode,T<: Real}

Calcualte the variance-covariance matrix from `tree`. An entry (i,j) of the matrix
is defined as the length of the path connecting the latest common ancestor
of i and j with the root of the tree.

Returns an Array of Real numbers.

* `tree` : Node in tree of interest.

* `blv` : branchlength vector of tree.

"""
function to_covariance(tree::N, blv::Vector{T})::Array{T,2} where {N<:GeneralNode,T<: Real}
    leaves::Vector{N} = get_leaves(tree)
    ll = length(leaves)
    covmat = zeros(T, ll, ll)
    #@inbounds for ((ind,itm),(jnd,jtm)) in Iterators.product(enumerate(leaves), enumerate(leaves))
    @inbounds for ind = 1:ll, jnd = 1:ind
        itm = leaves[ind]
        if ind == jnd
            covmat[ind,jnd] = covmat[ind,jnd] + reduce(+, @view blv[get_path(tree, itm)])
        else

            lca = find_lca(tree, itm, leaves[jnd])
            if !lca.root
                tmp = reduce(+, @view blv[get_path(tree, lca)])
                covmat[ind,jnd] = covmat[ind,jnd] + tmp
                covmat[jnd,ind] = covmat[jnd,ind] + tmp
            end # if
        end # if

    end # for
    covmat
end# function to_covariance

function to_covariance_func(tree::N)::Array{Function,2} where {N<: GeneralNode}
    leaves = get_leaves(tree)
    ll = length(leaves)
    covmat = Array{Function, 2}(undef, ll, ll)
    @inbounds for ind = 1:ll, jnd = 1:ind
        itm = leaves[ind]
        if ind == jnd
            sympath = get_path(tree, itm)
            covmat[ind,jnd] = y -> sum(y[sympath])
        else

            lca = find_lca(tree, itm, leaves[jnd])
            if !lca.root
                sympath = get_path(tree, lca)
                covmat[ind,jnd] = y -> sum(y[sympath])
                covmat[jnd,ind] = y -> covmat[ind,jnd](y)
            else
                covmat[ind,jnd] = y -> 0.0
                covmat[jnd,ind] = y -> 0.0
            end # if

        end # if

    end # for
    covmat
end# function to_covariance

function to_covariance_ultra(tree::Node) where T<: Real
    mv = tree_height(tree)
    blv = get_branchlength_vector(tree)
    blv ./= mv
    root = deepcopy(tree)
    set_branchlength_vector!(root, blv)
    force_ultrametric(root)

    leaves = get_leaves(root)
    ll = length(leaves)
    covmat = zeros(ll, ll)
    for i in 1:ll
        node1 = leaves[i]
        for j in 1:ll
            if i == j

                path = path_length(root, node1)
                covmat[i,j] = path
            elseif i>j
                lca = find_lca(root, node1, leaves[j])
                if lca.root != true
                    path = path_length(root, lca)
                    covmat[i,j] = path
                    covmat[j,i] = path

                end
            end
        end

    end
    covmat
end


text:
module PyTinkerforge

using Sockets
using PyCall
using Conda
using DocStringExtensions

function __init__()
  try
    tinkerforge = pyimport("tinkerforge")
  catch e
    Conda.pip_interop(true)
    Conda.pip("install", "Tinkerforge")
  end
end

include("ip_connection.jl")
include("devices/device_display_names.jl")
include("devices/device_factory.jl")

end


text:
push!(LOAD_PATH,"../src/") # necessary for unregistered packages

using Documenter
using StartUpDG

makedocs(
    sitename = "StartUpDG.jl",
    repo = "https://github.com/jlchan/StartUpDG.jl",
    modules=[StartUpDG],
    pages = [
        "Home" => "index.md",
        "Background and conventions" => "conventions.md",
        "`RefElemData`" => "RefElemData.md",
        "`MeshData`" => "MeshData.md",
        "Example: computing DG derivatives" => "ex_dg_deriv.md",
        "Timestepping" => "tstep_usage.md",
        "Reference" => "index_refs.md",
        "Authors" => "authors.md",
    ]
)

deploydocs(
    repo = "github.com/jlchan/StartUpDG.jl",
)


text:
module CoulombIntegrals

using LinearAlgebra

using LazyArrays
import LazyArrays: materialize
import ContinuumArrays: axes
import ContinuumArrays.QuasiArrays: AbstractQuasiArray, AbstractQuasiMatrix, MulQuasiArray, QuasiAdjoint
using IntervalSets

using CompactBases
import CompactBases: locs
using BlockBandedMatrices

using SpecialFunctions

const RadialOrbital{T,B<:AbstractQuasiMatrix} = Applied{<:Any, typeof(*), <:Tuple{B,<:AbstractVector{T}}}

include("analytical.jl")
include("lazy-coulomb.jl")
include("poisson.jl")
include("coulomb_repulsion_potential.jl")

end # module


text:
#ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля без внутренних перегородок и маркеров.
#РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки.

#r = Robot(animate = false, "1.sit")

function putmarker(r,side)
    while isborder(r,side) == false
        move!(r,side)
        putmarker!(r)
    end
end
function moverobot(r,side)
    while ismarker(r) == true
        move!(r,side)
    end
end
function reverse(side)
    for i=0:3
        if side == HorizonSide(i)
            return HorizonSide(mod(i+2, 4))
        end
    end
end
function krest(r)
    for side in (HorizonSide(i) for i=0:3)
        putmarker(r,side) #функция передвигает и ставит марки
        moverobot(r,reverse(side)) #функция перемещает в первоначальное положение, revers делает сторону противоположной
    end
    putmarker!(r)
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
function Jsingle!(J::Matrix{Float64},x::Vector{Float64},beta::Vector{Float64},rho::Float64)
    # J = zeros(4,4);
    J[1,1] = 1.;
    J[2,1] = rho.*x[1]./(x[2].^2);
    J[3,1] = 1./x[2];

    J[2,2] = (0.5.*beta[1].*x[2].^-0.5 .- rho.*x[1].^2.*x[2].^-3);
    J[3,2] = -x[1].*x[2].^-2 .+ sqrt.(0.5.*beta[1]./rho).*x[2].^-0.75;

    J[1,3] = -1.;
    J[2,3] = -rho.*x[3]./(x[4].^2);
    J[4,3] = 1./x[4];

    J[2,4] = (-0.5.*beta[2].*x[4].^-0.5 .+ rho.*x[3].^2.*x[4].^-3);
    J[4,4] = -x[3].*x[4].^-2 .- sqrt.(0.5.*beta[2]./rho).*x[4].^-0.75;

    # return J
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# This file contains the type and functions related to the Ricker model, also sets
# up the enviorment

# load packages
using Distributions
using DataFrames
using StatsBase
using Optim
#using Lasso
using GLM
using JLD
using HDF5

# load functions
#include("/home/samwiq/project 1 accalerated delated accaptance mcmc/gpmodel/gp_model.jl")
#include("/home/samwiq/project 1 accalerated delated accaptance mcmc/adaptive updating algorithms/adaptiveupdate.jl")
#include("/home/samwiq/project 1 accalerated delated accaptance mcmc/select case/selectcase.jl")

include(pwd()[1:end-13]*"/gpmodel/gp_model.jl")
include(pwd()[1:end-13]*"/adaptive updating algorithms/adaptiveupdate.jl")
include(pwd()[1:end-13]*"/select case/selectcase.jl")

# types:

"Type for prior distribution"
type PriorDistribution
  dist::String
  prior_parameters::Array{Float64}
end

"Type for the data"
type Data
  y::Array{Float64}
end

"Pararameters for the model"
type ModelParameters
  theta_true::Array{Float64} # [log(r) log(phi) log(sigma)]
  theta_known::Float64 # NaN
  theta_0::Array{Float64} # [log(r_0) log(phi_0) log(sigma_0)]
end

"Type for algorithm parameters for PMCMC"
type AlgorithmParameters
  R::Int64
  N::Int64
  burn_in::Int64
  alg::String
  pf_alg::String
  print_interval::Int64
end

"Type for algorithm parameters for gpPMCMC"
type AlgorithmParametersgpPMCMC
  R::Int64
  N::Int64
  burn_in::Int64
  length_training_data::Int64
  alg::String
  pf_alg::String
  print_interval::Int64
  est_method::String
  lasso::Bool
  pred_method::String
  independet_sampling::Bool
  noisy_est::Bool
  compare_GP_and_PF::Bool
  beta_MH::Float64
end

"Type for the problem (including algorithm parameters) for the PMCMC algorithm"
type Problem
  data::Data
  alg_param::AlgorithmParameters
  model_param::ModelParameters
  adaptive_update::AdaptationAlgorithm
  prior_dist::PriorDistribution
end

# this should be removed...
"Type for the problem (including algorithm parameters) for the gpPCMC algorithm"
type gpProblem
  data::Data
  alg_param::AlgorithmParametersgpPMCMC
  model_param::ModelParameters
  adaptive_update::AdaptationAlgorithm
  prior_dist::PriorDistribution
end

"Type for the results"
type Result
  Theta_est::Array{Float64}
  loglik_est::Array{Float64}
  accept_vec::Array{Float64}
  prior_vec::Array{Float64}
end

"Type for the results for the gpPCMC algorithm"
type gpResult
  Theta_est::Array{Float64}
  loglik_est::Array{Float64}
  accept_vec::Array{Float64}
  prior_vec::Array{Float64}
  compare_GP_PF::Array{Int64}
  data_gp_pf::Array{Float64}
  nbr_early_rejections::Int64
end

"type for the algorithm parameters for the ABC-MCMC algorithm"
type AlgorithmParametersABCMCMC
  R::Int64
  burn_in::Int64
  alg::String
  print_interval::Int64
  nbr_summary_stats::Int64
  eps::Vector
  start::String
end

"Type for the problem (including algorithm parameters) for the PMCMC algorithm"
type ProblemABCMCMC
  data::Data
  alg_param::AlgorithmParametersABCMCMC
  model_param::ModelParameters
  adaptive_update::AdaptationAlgorithm
  prior_dist::PriorDistribution
end

"Type for the results for the ABC-MCMC algorithm"
type ResultABCMCMC
  Theta_est::Array{Float64}
  accept_vec::Array{Float64}
  prior_vec::Array{Float64}
end

"Type for the parameters of the DA algorithm"
type DAParameters
  R::Int64 # nbr of DA iterations
  N::Int64 # length of training data
end

# load functions
include("mcmc.jl")

# set-up functions:

doc"""
    set_up_problem(use_sim_data::Bool=true;nbr_of_unknown_parameters::Int64=2,
  prior_dist::String = "Uniform", adaptiveon::Bool = true, ploton::Bool = false,
  nbr_of_cores::Int64 = 8, pf_algorithm::String = "bootstrap", T::Int64 = 50,x0::Float64 = 7.,
  print_interval::Int64 = 500, alg::String="PMCMC")

set_up_problem sets the parameters for the problem defined by the inputs.

# Output
* `problem`: An instance of the type Problem cointaining all information for the problem
"""
function set_up_problem(use_sim_data::Bool=true;nbr_of_unknown_parameters::Int64=2,
  prior_dist::String = "Uniform", adaptiveon::Bool = true, ploton::Bool = false,
  nbr_of_cores::Int64 = 8, pf_algorithm::String = "bootstrap", T::Int64 = 50,x0::Float64 = 7.,
  print_interval::Int64 = 500, alg::String="PMCMC")

  # set algorithm parameters
  theta_true = log([44.7; 10; 0.3])
  theta_0 =  log([3;3;10])
  Theta_parameters = [0 10; 0 4;-10 1]
  theta_known = NaN

  # create instance of AlgorithmParameters (set parameters to default values)
  alg_param = AlgorithmParameters(10000,500,5000,alg,pf_algorithm,print_interval)

  # create instance of ModelParameters, all theta paramters are on log-scale
  model_param = ModelParameters(theta_true,theta_known,theta_0)

  # set data
  if use_sim_data
    y = generate_data(T, exp(theta_true[1]),exp(theta_true[2]),exp(theta_true[3]),x0,ploton)
  end

  # create instance of Data
  data = Data(y)

  # create instance of PriorDistribution
  prior = PriorDistribution(prior_dist, Theta_parameters)

  # create instance of AdaptiveUpdate
  adaptive_update = AMUpdate(eye(length(theta_0)), 2.4/sqrt(length(theta_0)), 1., 0.7, 50)

  # return the an instance of Problem
  return Problem(data, alg_param, model_param, adaptive_update, prior)

end

doc"""
    set_up_gp_problem(use_sim_data::Bool=true;nbr_of_unknown_parameters::Int64=2,
  prior_dist::String = "Uniform", adaptiveon::Bool = true, ploton::Bool = false,
  nbr_of_cores::Int64 = 8, pf_algorithm::String = "bootstrap", est_method::String = "ml",lasso::Bool=true,
  independet_sampling::Bool = false, compare_GP_and_PF::Bool = true,noisy_est::Bool = true,
  pred_method::String = "sample",T::Int64 = 50,  x0::Float64 = 7., print_interval::Int64 = 500, nbr_predictions::Int64 = 25,
  selection_method::String="max_loglik", alg::String="PMCMC")

set_up_problem sets the parameters for the problem defined by the inputs.

# Output
* `problem`: An instance of the type gpProblem cointaining all information for the problem
"""
function set_up_gp_problem(use_sim_data::Bool=true;nbr_of_unknown_parameters::Int64=2,
  prior_dist::String = "Uniform", adaptiveon::Bool = true, ploton::Bool = false,
  nbr_of_cores::Int64 = 8, pf_algorithm::String = "bootstrap", est_method::String = "ml",lasso::Bool=true,
  independet_sampling::Bool = false, compare_GP_and_PF::Bool = true,noisy_est::Bool = true,
  pred_method::String = "sample",T::Int64 = 50,  x0::Float64 = 7., print_interval::Int64 = 500, alg::String="PMCMC", beta_MH::Float64=0.1)

  # set algorithm parameters
  theta_true = log([44.7; 10; 0.3])
  theta_0 =  log([3;3;10])
  Theta_parameters = [0 10; 0 4;-10 1]
  theta_known = NaN
  # create instance of AlgorithmParametersgpPMCMC (set parameters to default values)
  alg_param = AlgorithmParametersgpPMCMC(10000,500,2000,2000,alg,pf_algorithm,print_interval,est_method,lasso,pred_method,independet_sampling,noisy_est,compare_GP_and_PF,beta_MH)

  # create instance of ModelParameters, all theta paramters are on log-scale
  model_param = ModelParameters(theta_true,theta_known,theta_0)

  # set data
  if use_sim_data
    # load data from "data.csvexport_generated_data
    y = generate_data(T, exp(theta_true[1]),exp(theta_true[2]),exp(theta_true[3]),x0,ploton)
  end

 # create instance of Data
 data = Data(y)

  # create instance of PriorDistribution
  prior = PriorDistribution(prior_dist, Theta_parameters)

  # create instance of AdaptiveUpdate
  adaptive_update = AMUpdate(eye(length(theta_0)), 2.4/sqrt(length(theta_0)), 1., 0.7, 50)

  # return the an instance of gpProblem
  return gpProblem(data, alg_param, model_param, adaptive_update, prior)

end

#=
doc"""
    set_up_problem(use_sim_data::Bool=true;nbr_of_unknown_parameters::Int64=2,
  prior_dist::String = "Uniform", adaptiveon::Bool = true, ploton::Bool = false,
  nbr_of_cores::Int64 = 8, pf_algorithm::String = "bootstrap", T::Int64 = 50,x0::Float64 = 7.,
  print_interval::Int64 = 500, alg::String="PMCMC")

set_up_problem sets the parameters for the problem defined by the inputs.

# Output
* `problem`: An instance of the type Problem cointaining all information for the problem
"""
function set_up_abcmcmc_problem(use_sim_data::Bool=true;nbr_of_unknown_parameters::Int64=2,
  prior_dist::String = "Uniform", adaptiveon::Bool = true, ploton::Bool = false,
  nbr_of_cores::Int64 = 8, T::Int64 = 50,x0::Float64 = 7.,
  print_interval::Int64 = 500, alg::String="abcmcmc", nbr_summary_stats::Int64=5,
  eps::Float64=0.1)

  # set algorithm parameters
  theta_true = log([44.7; 10; 0.3])
  theta_0 =  log([3;3;10])
  Theta_parameters = [0 10; 0 4;-10 1]
  theta_known = NaN

  # create instance of AlgorithmParametersABCMCMC (set parameters to default values)
  alg_param = AlgorithmParametersABCMCMC(10000,2000,alg,print_interval,nbr_summary_stats,eps)

  # create instance of ModelParameters, all theta paramters are on log-scale
  model_param = ModelParameters(theta_true,theta_known,theta_0)

  # set data
  if use_sim_data
    y = generate_data(T, exp(theta_true[1]),exp(theta_true[2]),exp(theta_true[3]),x0,ploton)
  end

  # create instance of Data
  data = Data(y)

  # create instance of PriorDistribution
  prior = PriorDistribution(prior_dist, Theta_parameters)

  # create instance of AdaptiveUpdate
  adaptive_update = AMUpdate(eye(length(theta_0)), 2.4/sqrt(length(theta_0)), 1., 0.7, 50)

  # return the an instance of Problem
  return ProblemABCMCMC(data, alg_param, model_param, adaptive_update, prior)

end
=#

doc"""
    generate_data(T::Int65=50,r::Float64,phi::Float64,sigma::Float64,x0::Float64,ploton::Bool)

Generates data from the Ricker model.
"""
function generate_data(T::Int64,r::Float64,phi::Float64,sigma::Float64,x0::Float64=7.0,ploton::Bool=false)

  e = rand(Normal(0,sigma),T)
  y = zeros(T)
  x = zeros(T)

  # first iteration

  x[1] = r*x0*exp(-x0+e[1])

  y[1] = rand(Poisson(phi*x[1]))

  @simd for t = 2:T
    @inbounds x[t] = r*x[t-1]*exp(-x[t-1]+e[t])
    y[t] = rand(Poisson(phi*x[t]))
  end

  return y
end

max_stars_repo_path: Ricker model/rickermodel.jl

max_stars_repo_name: SamuelWiqvist/adamcmcpaper

max_issues_repo_path: Ricker model/rickermodel.jl

max_issues_repo_name: SamuelWiqvist/adamcmcpaper

max_forks_repo_path: Ricker model/rickermodel.jl

max_forks_repo_name: SamuelWiqvist/adamcmcpaper


text:
@testset "build" begin
    mktempdir() do dir
        files = map(1:2) do i
            without_extension = "file$i"
            file = "$(without_extension).jl"
            content = pluto_notebook_content("x = begin sleep(3); x = 3000 + $i; end")
            path = joinpath(dir, file)
            write(path, content)
            return file
        end
        parallel_build!(dir)

        html_file = joinpath(dir, "file1.html")
        @test contains(read(html_file, String), "3001")

        html_file = joinpath(dir, "file2.html")
        @test contains(read(html_file, String), "3002")
    end
end

@testset "invalid_notebook" begin
    try
        mktempdir() do dir
            content = pluto_notebook_content("@assert false")
            file = "file.jl"
            path = joinpath(dir, file)
            write(path, content)
            parallel_build!(dir, [file])
        end
        error("Test should have failed")
    catch AssertionError
        @test true # Success.
    end
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
using DataFrames
using CSV
using Turing
using StatsBase
using StatsPlots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %% 5.1, 5.2
d = CSV.read(datadir("exp_raw/WaffleDivorce.csv"), DataFrame)
d.D = zscore(d.Divorce)
d.M = zscore(d.Marriage)
d.A = zscore(d.MedianAgeMarriage)

std(d.MedianAgeMarriage)

# %% 5.3
# age of marriage
@model function divorce_A(A, D)
    a ~ Normal(0, 0.2)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA)
    D ~ MvNormal(μ, σ)
end

m5_1 = divorce_A(d.A, d.D)
prior = sample(m5_1, Prior(), 50) |> DataFrame

# %% 5.4
x = -2:0.1:2
plot()
for r in eachrow(prior)
    p = lin(r.a, x, r.bA)
    plot!(x, p, color = :black, alpha = 0.4)
end
plot!(legend = false)
plot!(xlabel = "Observed divorce", ylabel = "Predicted divorce")

# %% 5.5
q5_1 = quap(m5_1)
post = DataFrame(rand(q5_1.distr, 1000)', q5_1.params)

A_seq = range(-3, 3.2, length = 30)
mu = lin(post.a', A_seq, post.bA') |> meanlowerupper

scatter(d.A, d.D, alpha = 0.4, legend = false)
plot!(A_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))

# %% 5.6
# marriage rate
@model function divorce_M(M, D)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, M, bM)
    D ~ MvNormal(μ, σ)
end

q5_2 = quap(divorce_M(d.M, d.D))
post = DataFrame(rand(q5_2.distr, 1000)', q5_2.params)

M_seq = range(-3, 3.2, length = 30)
mu = lin(post.a', M_seq, post.bM') |> meanlowerupper

scatter(d.M, d.D, alpha = 0.4, legend = false)
plot!(M_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))

# %% 5.7 - 5.9
missing  # TODO

# %% 5.10
# both age at and rate of marriage
@model function divorce_A_M(A, M, D)
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA, M, bM)
    D ~ MvNormal(μ, σ)
end

m5_3 = divorce_A_M(d.A, d.M, d.D)
q5_3 = quap(m5_3)
DataFrame(rand(q5_3.distr, 1000)', q5_3.params) |> precis

# %% 5.11
# Well, this ain't pretty. But before I have something reasonable, I at least wanted to
# something, even if it's messily thrown together. TODO

x = [
    q5_1.coef.bA,       # bA
    NaN,
    q5_3.coef.bA,
    NaN,                # Spaceholder
    NaN,                # bM
    q5_2.coef.bM,
    q5_3.coef.bM,
]
xerr = sqrt.([
    q5_1.vcov[2, 2],    # bA
    NaN,
    q5_3.vcov[3, 3],
    NaN,                # Spaceholder
    NaN,                # bM
    q5_2.vcov[2, 2],
    q5_3.vcov[2, 2],
])
ylab = [
    "bA m5.1",
    "bA m5.2",
    "bA m5.3",
    "",
    "bM m5.1",
    "bM m5.2",
    "bM m5.3",
]

scatter(x, ylab, xerr = xerr, legend = false)
vline!([0])

# %% 5.12
N = 50
age = randn(N)
mar = rand.(Normal.(-age))
div = rand.(Normal.(age))

# %% 5.13, 5.14
@model function divorce_AM(A, M)
    a ~ Normal(0, 0.2)
    bAM ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bAM)
    M ~ MvNormal(μ, σ)
end

sort!(d, :A)  # for nice looking graphs

q5_4 = quap(divorce_AM(d.A, d.M))
post = DataFrame(rand(q5_4.distr, 1000)', q5_4.params)
mu = lin(post.a', d.A, post.bAM') |> meanlowerupper
resid = d.M .- mu.mean

scatter(d.A, d.M, legend = false)
plot!(d.A, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))
plot!(xlabel = "Age at marriage (std)", ylabel = "Marriage rate (std)")

# %% 5.15, 5.16
post = DataFrame(rand(q5_3.distr, 1000)', q5_3.params)
μ = lin(post.a', d.A, post.bA', d.M, post.bM') |> meanlowerupper

scatter(d.D, μ.mean, yerror = (μ.mean .- μ.lower, μ.upper .- μ.mean), legend = false)
plot!(identity, range(extrema(d.D)..., length = 10))
plot!(xlabel = "Observed divorce", ylabel = "Predicted divorce")

# %% 5.17
xy = [d.D μ.mean]
for loc in ("ID", "UT", "ME", "RI")
    coord = xy[d.Loc .== loc, :]
    annotate!(coord..., loc, :bottom)
end
plot!()

# %% 5.18
N = 100
x_real = randn(N)
x_spur = rand.(Normal.(x_real))
y = rand.(Normal.(x_real))
d = DataFrame((; y, x_real, x_spur))  # the semicolon turns it into a NamedTuple which then
                                      # gives the DataFrame the names of the columns

@df d corrplot([:y :x_real :x_spur])  # either
corrplot(Matrix(d), label = names(d)) # or

# %% 5.19
d = DataFrame(CSV.File(datadir("exp_raw/WaffleDivorce.csv")))
d.D = zscore(d.Divorce)
d.M = zscore(d.Marriage)
z_age, unz_age = zscore_transform(d.MedianAgeMarriage)
d.A = z_age(d.MedianAgeMarriage)

@model function divorce_19(A, M, D)
    # A -> M
    aM ~ Normal(0, 0.2)
    bAM ~ Normal(0, 0.5)
    σ_M ~ Exponential(1)
    μ_M = lin(aM, A, bAM)
    M ~ MvNormal(μ_M, σ_M)
    # A -> D <- M
    a ~ Normal(0, 0.2)
    bM ~ Normal(0, 0.5)
    bA ~ Normal(0, 0.5)
    σ ~ Exponential(1)
    μ = lin(a, A, bA, M, bM)
    D ~ MvNormal(μ, σ)
end

q5_3_A = quap(divorce_19(d.A, d.M, d.D))

# %% 5.20 - 5.22
A_seq = range(-2, 2, length = 30)
post = DataFrame(rand(q5_3_A.distr, 1000)', q5_3_A.params)

μM = lin(post.aM', A_seq, post.bAM')
M = rand.(Normal.(μM, post.σ_M'))
μD = lin(post.a', A_seq, post.bA', M, post.bM')
D = rand.(Normal.(μD, post.σ')) |> meanlowerupper

plot(A_seq, D.mean, ribbon = (D.mean .- D.lower, D.upper .- D.mean))
plot!(legend = false, xlabel = "manipulated A", ylabel = "counterfactual D")

# %% 5.23
# I use zscore_transform from the tools.jl file so I don't have to hardcode the mean and
# stddev of d.MedianAgeMarriage.
M = rand.(Normal.(lin(post.aM', (20.0, 30.0) |> z_age, post.bAM'), post.σ_M'))
μD = lin(post.a', (20.0, 30.0) |> z_age, post.bA', M, post.bM') |> meanlowerupper
μD.mean[2] - μD.mean[1]
# Eh, yes, 4.5 stddev are kinda big but so is a change in the median age of marriage from
# 20 to 30 years, clocking in at 8 stddevs. I mean, think about it, a median age of
# marriage of 20 years?!

# %% 5.24
M_seq = range(-2, 2, length = 30)
A_seq = zeros(length(M_seq))
M = rand.(Normal.(lin(post.aM', A_seq, post.bAM'), post.σ_M'))
μD = lin(post.a', A_seq, post.bA', M_seq, post.bM')
D = rand.(Normal.(μD, post.σ')) |> meanlowerupper

plot(M_seq, D.mean, ribbon = (D.mean .- D.lower, D.upper .- D.mean))
plot!(legend = false, xlabel = "manipulated M", ylabel = "counterfactual D")

# %% 5.25 - 5.27
# This is a little superfluous since we aren't hiding any details but here we go.
A_seq = range(-2, 2, length = 30)
post = DataFrame(rand(q5_3_A.distr, 1000)', q5_3_A.params)
M_sim = rand.(Normal.(post.aM' .+ A_seq .* post.bAM', post.σ_M'))
D_sim = rand.(Normal.(post.a' .+ A_seq .* post.bA', post.σ'))

max_stars_repo_path: scripts/5_1 Spurious association.jl

max_issues_repo_path: scripts/5_1 Spurious association.jl

max_forks_repo_path: scripts/5_1 Spurious association.jl


text:
function place_init_grid_order(client, config)
    val = parse(Float64, config["total_value"]) * parse(Float64, config["init_state"]["init_fraction"])/100
    price = client.get_ticker(symbol = config["symbol"])["lastPrice"]
    current = parse(Float64, price)
    qty = floor(val / current, digits=config["precision"])

    if config["testing"] == "1"
        order = try client.create_test_order(
            symbol=config["symbol"],
            side=client.SIDE_BUY,
            type=client.ORDER_TYPE_LIMIT,
            timeInForce=client.TIME_IN_FORCE_GTC,
            quantity=qty,
            price=price
        );
        catch
            nothing
        end
        if order !== nothing
            order = Dict(
                "price" => price,
                "orderId" => 42,
                "executedQty" => "0.00",
                "origQty" => string(qty),
                "side" => client.SIDE_BUY,
                "type" => "limit",
                "symbol" => config["symbol"],
                "status" => "NEW",
                "cummulativeQuoteQty" => "0.00000000",
            )
        else
            return nothing
        end
    else
        order = try client.order_limit_buy(
            symbol=config["symbol"],
            quantity=qty,
            price=price
        );
        catch
            nothing
        end
        if order === nothing
            return nothing
        end
    end
    order["time"] = string(now())
    return order
end
export place_init_grid_order

function check_init_grid_order(client, config)
    if config["testing"] == "1"
        ticker_file = joinpath(Pinguu.base_dir, "wp_data/ticker/" * config["symbol"] * "_kline_1m_candle.json")
        ticker_data = JSON.parsefile(ticker_file);
        ticker = Pinguu.price2string(ticker_data[end][end], config["sell_precision"])
        tmp = deepcopy(config["init_state"]["order"])
        @info("Current: " * ticker)
        @info("Buy: " * tmp["price"])
        if parse(Float64, tmp["price"]) >= parse(Float64, ticker)
            tmp["status"] = client.ORDER_STATUS_FILLED
            tmp["executedQty"] = tmp["origQty"]
            tmp["cummulativeQuoteQty"] = round(parse(Float64, tmp["price"]) * parse(Float64, tmp["origQty"]), digits=config["sell_precision"])
        end
    else
        tmp = try client.get_order(symbol=config["symbol"], orderId=config["init_state"]["order"]["orderId"])
        catch
            return nothing
        end
    end
    return tmp
end
export check_init_grid_order


text:
########################## NIGP (Noisy Input GP) ###############################
# IDEA: training procedure for input variances (chain derivatives through corrective term)
# Laplace approximation to treating measurement points as multivariate Gaussians?

# conditional with two covariance inputs inspired by approximation put forth in
# "Gaussian Process Training with Input Noise" by McHutchon and Rasmussen 2011
# generalizes the work by allowing non-constant input noise variances
# ASSUMPTION: x and y are vectors of data, not data vectors
function conditional(G::Gaussian, Gx::Gaussian{<:AbstractVector, <:Diagonal},
                                  Gy::Gaussian{<:AbstractVector, <:Diagonal})
    x, y = mean(Gx), mean(Gy)
    Σx::Diagonal, Σy::Diagonal = cov(Gx), cov(Gy)
    C = conditional(G, x, y, Σy) # canonical conditioning
    ∇(x::Real) = ForwardDiff.derivative(mean(C), x) # differentiate posterior mean
    ∇(x::AbstractVector) = ForwardDiff.gradient(mean(C), x)
    Σ = Diagonal(zeros(length(x)))
    for i in eachindex(x)
        ∇x = ∇(x[i])
        Σ[i, i] = max(dot(∇x, Σx[i, i], ∇x), 0)
    end
    Σ += Σy
    conditional(G, x, y, Σ)
end

# WARNING: Σx is the d by d noise covariance for all inputs x,
# similar to σy² being the noise variance for all outputs y
function conditional(G::Gaussian,
                     x::AbstractVector, Σx::Union{Real, AbstractMatrix},
                     y::AbstractVector, σy²::Real)
    n = length(x)
    Σx = Diagonal(fill(Σx, n))
    Gx = Gaussian(x, Σx)
    Σy = σy² * I(n)
    Gy = Gaussian(y, Σy)
    conditional(G, Gx, Gy)
end


text:
# All functions related to dispersal:
# - dispersal
# - establishment
# - viability checking

"""
    disperse!(world, static)

Dispersal of individuals within the world.
"""
function disperse!(world::Array{Patch,1}, static::Bool = true)
    # TODO: additional border conditions, refactor
    for patch in world
        idx = 1
        while idx <= size(patch.seedbank,1)
            dispmean = patch.seedbank[idx].traits["dispmean"]
            dispshape = patch.seedbank[idx].traits["dispshape"]
            xdir = rand([-1,1]) * rand(Logistic(dispmean,dispshape))/sqrt(2) # scaling so that geometric mean...
            ydir = rand([-1,1]) * rand(Logistic(dispmean,dispshape))/sqrt(2) # ...follows original distribution
            xdest = patch.location[1] + Int(round(xdir))
            ydest = patch.location[2] + Int(round(ydir))
            !patch.isisland ? target = checkborderconditions(world, xdest, ydest) : target = (xdest, ydest)
            possdest = findall(x -> in(x.location, [target]), world)
            static && filter!(x -> world[x].isisland, possdest) # disperse only to islands
            if !static || patch.isisland
                indleft = splice!(patch.seedbank,idx) # only remove individuals from islands!
                idx -= 1
            end
            if length(possdest) > 0 # if no viable target patch, individual dies
                if static && !patch.isisland
                    indleft = patch.seedbank[idx]
                end
                destination = rand(possdest) # currently there is only one possible destination
                push!(world[destination].community, indleft)
            end
            idx += 1
        end
    end
end

"""
    establish!(patch, nniches)

Establishment of individuals in patch `p`: Sets the adaptation parameters (~fitness)
according to an individual's adaptation to the niches of the surrounding environment.

A maximum of two niches (temperature and "precipitation") is currently supported.
"""
function establish!(patch::Patch, nniches::Int=1)
    temp = patch.temp
    idx = 1
    while idx <= size(patch.community,1)
        if patch.community[idx].marked
            opt = patch.community[idx].traits["tempopt"]
            tol = patch.community[idx].traits["temptol"]
            fitness = gausscurve(opt, tol, temp, 0.0)
            fitness > 1 && (fitness = 1) # should be obsolete
            fitness < 0 && (fitness = 0) # should be obsolete
            patch.community[idx].tempadaptation = fitness
            if nniches >= 2
                opt = patch.community[idx].traits["precopt"]
                tol = patch.community[idx].traits["prectol"]
                fitness = gausscurve(opt, tol, patch.prec, 0.0)
                fitness > 1 && (fitness = 1) # should be obsolete
                fitness < 0 && (fitness = 0) # should be obsolete
                patch.community[idx].precadaptation = fitness
            end
            patch.community[idx].marked = false
        end
        idx += 1
    end
end

"""
    establish!(world, nniches, static)

Carry out establishment for each patch in the world.
"""
function establish!(world::Array{Patch,1}, nniches::Int=1, static::Bool = true)
    for patch in world
        (patch.isisland || !static) && establish!(patch, nniches) # pmap(!,patch) ???
    end
end

"""
    checkviability!(community, settings)

Check whether all individuals in the passed community conform to a basic set of
constraints (i.e. all traits are present and certain properties are >= 0).
Individuals that fail the test are removed from the community.
"""
function checkviability!(community::Array{Individual, 1}, settings::Dict{String, Any})
    idx=1
    while idx <= size(community,1)
        reason = ""
        dead = false
        community[idx].size <= 0 && (dead = true) && (reason *= "size ")
        any(collect(values(community[idx].traits)) .< 0) && (dead = true) && (reason *= "traitvalues ")
        community[idx].traits["repsize"] <= community[idx].traits["seedsize"] && (dead = true) && (reason *= "seed/rep ")
        community[idx].tempadaptation < 0 && (dead = true) && (reason *= "fitness ")
        community[idx].precadaptation < 0 && (dead = true) && (reason *= "fitness ")
        community[idx].traits["selfing"] > 1 && (dead = true) && (reason *= "selfing ")
        community[idx].traits["seqsimilarity"] > 1 && (dead = true) && (reason *= "seqsimilarity ")
        !traitsexist(community[idx].traits, settings) && (dead = true) && (reason *= "missingtrait ")
        if dead
            simlog("Individual not viable: $reason. Being killed.", 'w')
            splice!(community,idx)
            continue
        end
        idx += 1
    end
end

"""
    checkviability(world, settings)

Check the viability of all individuals.
"""
function checkviability!(world::Array{Patch,1}, settings::Dict{String, Any})
    for patch in world
        #XXX pmap(checkviability!,patch) ???
        checkviability!(patch.community, settings)
    end
end

"""
    traitsexist(traits, settings)

Check a trait dict to make sure it contains the full set of traitnames required
by the model (as defined in the settings).
"""
function traitsexist(traits::Dict{String, Float64}, settings::Dict{String, Any})
    missingtraits = setdiff(settings["traitnames"], keys(traits))
    if length(missingtraits) > 0
        simlog("Missing trait $missingtraits. Individual might be killed.", 'w')
        return false
    end
    true
end

"""
    traitsexist(individual, settings)

Make sure an individual organism has the full set of traits required by the model
(as defined in the settings).
"""
function traitsexist(ind::Individual, settings::Dict{String, Any})
    traitnames = settings["traitnames"]
    for trait in traitnames
        if !haskey(ind.traits, trait)
            simlog("Individual is missing trait $trait. Might be killed.", 'e')
            return false
        end
    end
    true
end

"""
    gausscurve(b, c, x, a=1.0)

Calculate the value of the Gauss function ("bell curve") at point x; with
a being the maximum height of the curve, b the position of the curve center and
c the standard deviation ("width").
"""
function gausscurve(b, c, x, a = 1.0)
    if c != 0 && a != 1.0
        a = 1 / (c * sqrt(2 * pi))
        y = a * exp(-(x-b)^2/(2*c^2))
    elseif c != 0
        y = a * exp(-(x-b)^2/(2*c^2))
    else
        y = 0.0
    end
end

"""
    findisland(w)

within world `w`, find out in which direction from the continent the island(s) lie(s).
"""
function findisland(world::Array{Patch,1})
    xmin = minimum(map(x->x.location[1],world))
    xmax = maximum(map(x->x.location[1],world))
    ymin = minimum(map(x->x.location[2],world))
    ymax = maximum(map(x->x.location[2],world))
    westernborder = filter(x->x.location[1]==xmin,world)
    northernborder = filter(x->x.location[2]==ymax,world)
    easternborder = filter(x->x.location[1]==xmax,world)
    southernborder = filter(x->x.location[2]==ymin,world)
    if all(map(x->x.isisland,westernborder))
        return "west"
    elseif all(map(x->x.isisland,northernborder))
        return "north"
    elseif all(map(x->x.isisland,easternborder))
        return "east"
    elseif all(map(x->x.isisland,southernborder))
        return "south"
    else
        return "none"
    end
end

"""
    checkborderconditions!(w, x, y)

check if coordinates `x` and `y` lie within world `w` and correct if not,
considering defined border conditions.
"""
function checkborderconditions(world::Array{Patch,1}, xdest::Int, ydest::Int)
    xmin = minimum(map(x->x.location[1],world))
    xmax = maximum(map(x->x.location[1],world))
    ymin = minimum(map(x->x.location[2],world))
    ymax = maximum(map(x->x.location[2],world))
    xrange = xmax - xmin + 1 # we're counting cells!
    yrange = ymax - ymin + 1 # we're counting cells!
    xshift = xdest - xmin + 1 # 1-based count of cells
    yshift = ydest - ymin + 1 # 1-based count of cells
    xshift > 0 ? outofx = abs(xshift) : outofx = abs(xshift) + 1
    while outofx > xrange
        outofx -= xrange
    end
    outofx -= 1
    yshift > 0 ? outofy = abs(yshift) : outofy = abs(yshift) + 1
    while outofy > yrange
        outofy -= yrange
    end
    outofy -= 1
    islanddirection = findisland(world::Array{Patch,1})
    if islanddirection == "west"
        xdest > xmax && (xdest = xmax - outofx) # east: reflective
        ydest < ymin && (ydest = ymax - outofy) # south: periodic
        ydest > ymax && (ydest = ymin + outofy) # north: periodic
    elseif islanddirection == "north"
        ydest < ymin && (ydest = ymin + outofy) # south: reflective
        xdest < xmin && (xdest = xmax - outofx) # west: periodic
        xdest > xmax && (xdest = xmin + outofx) # east: periodic
    elseif islanddirection == "east"
        xdest < xmin && (xdest = xmin + outofx) # west: reflective
        ydest < ymin && (ydest = ymax - outofy) # south: periodic
        ydest > ymax && (ydest = ymin + outofy) # north: periodic
    elseif islanddirection == "south"
        ydest > ymax && (ydest = ymax + outofy) # north: reflective
        xdest < xmin && (xdest = xmax - outofx) # west: periodic
        xdest > xmax && (xdest = xmin + outofx) # east: periodic
    else
        ydest > ymax && (ydest = ymin + outofy) # north: periodic
        xdest > xmax && (xdest = xmin + outofx) # east: periodic
        ydest < ymin && (ydest = ymax - outofy) # south: periodic
        xdest < xmin && (xdest = xmax - outofx) # west: periodic
    end
    xdest, ydest
end

max_stars_repo_name: lleiding/gemm

max_issues_repo_name: lleiding/gemm

max_forks_repo_name: lleiding/gemm


text:
# Creates a tex file with 35 columns of state_total
# for use in plotCorner_ttv.py:

using JLD2
using DelimitedFiles

@load "../../../data/T1_hmc_total_02212020.jld2" state_total

writedlm("../../../data/state_total.txt",state_total[1:35,:]')


text:
#sinekernel.jl

n=200
x=(0:.025:1)*3*pi
include("hermitekernel.jl")
K=hermitekernel(n-1,x/sqrt(2*n))*pi/sqrt(2*n)

## Plot
using Winston
module Plot3D
    using Winston
    # Warning: the 3D interface is very new at the time
    # of writing this and is likely to not be available
    # and/or have changed
    if Winston.output_surface == :gtk
        include(Pkg.dir("Winston","src","canvas3d_gtk.jl"))
    else
        include(Pkg.dir("Winston","src","canvas3d.jl"))
    end
end
xx = Float64[xx for xx = x, yy = x]
yy = Float64[yy for xx = x, yy = x]
Plot3D.plot3d(Plot3D.surf(xx,K,yy))

Plot3D.plot3d(Plot3D.surf(
    (u,v)->u,
    (u,v)->sin(u-v)./(u-v),
    (u,v)->v,
    x, x))


text:
#!/usr/bin/julia

function base10ToAny(in_dig,base)
    if base >= 21
        println("Warning: There're not special characters for digits in base 21 or more")
    end
    digits_b10=[]
    q=in_dig
    while true
        if q % base <= 9
            push!(digits_b10,q % base)
        elseif q % base == 10
            push!(digits_b10,'A')
        elseif q % base == 11
            push!(digits_b10,'B')
        elseif q % base == 12
            push!(digits_b10,'C')
        elseif q % base == 13
            push!(digits_b10,'D')
        elseif q % base == 14
            push!(digits_b10,'E')
        elseif q % base == 15
            push!(digits_b10,'F')
        elseif q % base == 16
            push!(digits_b10,'G')
        elseif q % base == 17
            push!(digits_b10,'H')
        elseif q % base == 18
            push!(digits_b10,'I')
        elseif q % base == 19
            push!(digits_b10,'J')
        else
            r=(q % base)
            push!(digits_b10,"($r)")
        end
        q=div(q,base)
        if q == 0
            break
        end
    end
    l=length(digits_b10)
    to_ret=[]
    for idx in 1:l
        push!(to_ret,pop!(digits_b10))
    end
    if base <= 10
        return parse(Int,join(to_ret))
    else
        return join(to_ret)
    end
end

function anyBaseTo10(iterDig,base)
    if base >= 21
        println("Error: Only works for a base 20 or less")
        return
    end
    if !isa(iterDig,AbstractString)
        println("Error: Input n should be a string")
        return
    end
    summation=0
    l=sizeof(iterDig)
    for idx in 1:l
        digit = iterDig[ l - (idx-1)]
        if digit == 'A'
            digit = 10
        elseif digit == 'B'
            digit = 11
        elseif digit == 'C'
            digit = 12
        elseif digit == 'D'
            digit = 13
        elseif digit == 'E'
            digit = 14
        elseif digit == 'F'
            digit = 15
        elseif digit == 'G'
            digit = 16
        elseif digit == 'H'
            digit = 17
        elseif digit == 'I'
            digit = 18
        elseif digit == 'J'
            digit = 19
        else
            digit = parse(Int64, digit)
        end
        summation = digit*(base^(idx-1)) + summation
    end
    return summation
end

function binAddition(in_bin_1,in_bin_2)
    in_bin_1=base10ToAny(in_bin_1,2)
    in_bin_2=base10ToAny(in_bin_2,2)
    ib1=string(in_bin_1)
    ib2=string(in_bin_2)
    l1=sizeof(ib1)
    l2=sizeof(ib2)
    if l1 > l2
        AddZeros = l1 - l2
        ib2=("0" ^ AddZeros) * ib2
        larger=l1
    elseif l2 > l1
        AddZeros = l2 - l1
        ib1=("0" ^ AddZeros) * ib1
        larger=l2
    else
        larger=l1
    end
    c=0
    all_s=[]
    for idx in 1:larger
        dig1= parse(Int64, ib1[ larger - (idx-1) ])
        dig2= parse(Int64, ib2[ larger - (idx-1) ])
        d=div( dig1 + dig2 + c, 2)
        s=dig1 + dig2 + c - 2d
        c=d
        push!(all_s,s)
    end
    push!(all_s,c)
    justZero = true
    catResult=[]
    for idx in 1:length(all_s)
        if length(all_s) == 2 && join(all_s) == "00"
            return 0
        end
        if justZero && all_s[ length(all_s) - (idx -1 ) ] == 0
            continue
        end
        if all_s[ length(all_s) - (idx -1 ) ]  == 1
            justZero = false
        end
        push!(catResult,all_s[ length(all_s) - (idx -1 ) ])
    end
    return anyBaseTo10(join(catResult),2)
end

function binMultiplication(a,b)
    a=base10ToAny(a,2)
    b=base10ToAny(b,2)
    a=string(a)
    b=string(b)
    la=sizeof(a)
    lb=sizeof(b)
    if la > lb
        AddZeros = la - lb
        b=("0" ^ AddZeros) * b
        larger=la
    elseif lb > la
        AddZeros = lb - la
        a=("0" ^ AddZeros) * a
        larger=lb
    else
        larger=la
    end
    c_partials=[]
    for idx in 1:larger
        if b[larger - (idx -1)] == '1'
            c = a * ("0" ^ (idx - 1))
            push!(c_partials,c)
        else
            c = "0"
            push!(c_partials,c)
        end
    end
    p = 0
    for partial in c_partials
        partial = parse(Int64, partial)
        p = binAddition(p,anyBaseTo10(string(partial),2))
    end
    return p
end

# Division Algorithm
# a = d(q) + r
function intDivision(a::Int,d::Int)
    q = 0
    r = abs(a)
    while r >= d
        r = r - d
        q = q + 1
    end
    if a < 0 && r > 0
        r = d - r
        q = -(q + 1)
    end
    return [q,r]
end

# Modular Exponentiation
# x = b ^ n mod m
function modExp(b::Int, n::Int, m::Int)
    n=base10ToAny(n,2) #convert to binary
    x=1
    power=mod(b,m)
    ns=string(n)
    for idx in 1:sizeof(ns)
        if ns[sizeof(ns) - (idx-1)] == '1'
            x = mod(x * power,m)
        end
        power = mod((power * power),m)
    end
    return x
end

install.packages("plotly")
# The Euclidean Algorithm
function gcd(a::Int,b::Int)
    x=a
    y=b
    while y != 0
        r = mod(x,y)
        x = y
        y = r
    end
    return x
end

function sieveEratosthenes(n::Int) # accepts one integer argument
    isprime = ones(Bool, n) # n-element vector of true-s
    isprime[1] = false
    # 1 is not a prime
    for i in 2:round(Int, sqrt(n)) # loop integers from 2 to sqrt(n)
        if isprime[i]
            # conditional evaluation
            for j in (i*i):i:n # sequence with step i
            isprime[j] = false
            end
        end
    end
    return filter(x -> isprime[x], 1:n) # filter using anonymous function
end

max_stars_repo_name: elbecerrasoto/tutti-frutti

max_issues_repo_name: elbecerrasoto/tutti-frutti

max_forks_repo_name: elbecerrasoto/tutti-frutti


text:
@polly function kernel_symm(alpha, beta, C, A, B)
    m,n = size(C)
    for i = 1:m, j = 1:n
        temp2 = zero(eltype(C))
        for k = 1:(i-1)
            C[k,j] += alpha * B[i,j] * A[i,k]
            temp2 += B[k,j] * A[i,k]
        end
        C[i,j] = beta * C[i,j] + alpha * B[i,j] * A[i,i] + alpha * temp2
    end
end

let
    m = 1000
    n = 1200

    C = zeros(Float32,m,n)
    A = zeros(Float32,m,m)
    B = zeros(Float32,m,n)

    alpha = Float32(1.5)
    beta = Float32(1.2)
    for i = 1:m, j = 1:n
        C[i,j] = ((i+j) % 100) / m
        B[i,j] = ((n+i-j) % 100) / m
    end
    for i = 1:m
        for j = 1:i
            A[i,j] = ((i+j) % 100) / m
        end
        for j = i+1:m
            A[i,j] = -999 # regions of arrays that should not be used
        end
    end

    SUITE["symm"] = @benchmarkable kernel_symm($alpha, $beta, C, A, B) setup = (C = copy($C); A = copy($A); B = copy($B))
end


text:
using Fermi.HartreeFock

import Base: show

export RMP2

"""
    Fermi.HartreeFock.RHFAlgorithm

Abstract type for RMP2 implementations.
"""
abstract type RMP2Algorithm end

"""
    Fermi.MollerPlesset.get_rmp2_alg()

Returns a singleton type corresponding to a RMP2 implementation based on the options.
"""
function get_rmp2_alg()
    implemented = [RMP2a()]
    N = Options.get("mp2_alg")
    try
        return implemented[N]
    catch BoundsError
        throw(FermiException("implementation number $N not available for RMP2."))
    end
end

"""
    RMP2{T} <: AbstractMPWavefunction

    TODO
"""
struct RMP2{T} <: AbstractMPWavefunction
    correlation::T
    energy::AbstractFloat
end

# For each implementation a singleton type must be create
struct RMP2a <: RMP2Algorithm end
include("RMP2a.jl")

function RMP2(x...)
    if !any(i-> i isa RMP2Algorithm, x)
        RMP2(x..., get_rmp2_alg())
    else
        # Print the type of arguments given for a better feedback
        args = "("
        for a in x[1:end-1]
            args *= "$(typeof(a)), "
        end
        args = args[1:end-2]*")"
        throw(FermiException("invalid arguments for RMP2 method: $args"))
    end
end

## MISCELLANEOUS
# Pretty printing
function string_repr(X::RMP2)
    out = ""
    out = out*" ⇒ Fermi Restricted MP2 Wave function\n"
    out = out*" ⋅ Correlation Energy:     $(X.correlation)\n"
    out = out*" ⋅ Total Energy:           $(X.energy)"
    return out
end

function show(io::IO, ::MIME"text/plain", X::RMP2)
    print(io, string_repr(X))
end