text:
# avoiding-global-variables.jl

# Let's start with a little demonstration of the differences in speed
# between two basic things. We're going to be using @time to show this
# off.

# BUT FIRST.

# We need to talk about the fuckery that is... Global scope.

########################################################################

println("------------------------")

# Julia's global scoping is freaking weird. For example, this code works
# as-is.

x = 1

for i in 1:3
    println(x, "\t", i)
end

println("----------------")

# # *This* code, meanwhile, will wig out on you if you uncomment it. It'll
# # give you an "UndefVarError: x not defined", at least at the time of
# # writing.
#
# x = 2
#
# for i in 1:3
#     x = x + 2
#     println(x, "\t", i)
# end
#
# # Fucking. Infuriating.

# What's going on here? It turns out that when you open up a loop like a
# for loop, Julia's scope rules get weirdly exclusionary about
# everything. If you want to use that x, you need to first pull it in
# with an explicit "global x" line. Like so:

x = 3

for i in 1:3
    global x
    x = x + 3
    println(x, "\t", i)
end

# And that works just fine. It's one of the most annoying gotchas on the
# planet.

# Except you NEVER EVER EVER want to do that, because global variables
# 💩 the 🛌 in Julia. Let's watch!

println("------------------------")

########################################################################

println("-.-.-.-.-.-.-.-.-.-.-.-.")

# Here is a loop that does something trivially simple: It adds 1 to p, 1,000,000 times.

begin
    p = 1
    for i in 1:1000000
        global p           # REMEMBER THIS PART OR SUFFER <==
        p += 1
    end
end

# And here is that exact same code, wrapped inside a function call.

function f()
    begin
        p = 1
        for i in 1:1000000  # NO GLOBAL NEEDED ANYMORE, GOOD RIDDANCE.
            p += 1
        end
    end
end

f()

# You would expect these two to have roughly the same performance,
# right? Oh, you poor, naive fool. If only you knew.

println("-.-.-.-.-.-.-.")

# Let's bust out the @time macro, which has been nicely baked into the
# Base package, and see how the two code blocks *actually* turn out.

@time begin
    p = 1
    for i in 1:1000000
        global p           # REMEMBER THIS PART OR SUFFER <==
        p += 1
    end
end

function f()
    begin
        p = 1
        for i in 1:1000000  # NO GLOBAL NEEDED ANYMORE, GOOD RIDDANCE.
            p += 1
        end
    end
end

@time f()

println("-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.")

# On my machine, the second one is almost an *order of magnitude* faster
# than the first.

# Not convinced? Let's ramp things up to the point where we human beings
# can actually feel the difference. This code runs 10,000,000 times.

@time begin
    p = 1
    for i in 1:10000000
        global p           # REMEMBER THIS PART OR SUFFER <==
        p += 1
    end
end

function f()
    begin
        p = 1
        for i in 1:10000000  # NO GLOBAL NEEDED ANYMORE, GOOD RIDDANCE.
            p += 1
        end
    end
end

@time f()

println("-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.")
# And this code? 100,000,000 times.

@time begin
    p = 1
    for i in 1:100000000
        global p           # REMEMBER THIS PART OR SUFFER <==
        p += 1
    end
end

function f()
    begin
        p = 1
        for i in 1:100000000  # NO GLOBAL NEEDED ANYMORE, GOOD RIDDANCE.
            p += 1
        end
    end
end

@time f()

# Take a look at the @time performance for all those little
# global assignments. They keep going up. On my machine,
# they're going up from 0.05 seconds, to 0.5 seconds, to 5
# whole seconds. The allocations go up from 1 M allocations,
# to 10 M allocations, to 100 M allocations.

# Meanwhile, the function takes almost the same amount of time *and* allocations every single time.

# Global variables are WACK. Avoid them if you at all can.

########################################################################

println("====================================")
println("TIME TO IMPORT PLOTS and DIFF EQ ::")
# Alright, big woop. What does that have to do with what we're doing?

# I'm glad you asked. Let's look at some example code.

@time using Plots
@time using DifferentialEquations

########################################################################

println("❤❤❤❤❤❤❤❤❤❤❤❤❤❤❤❤❤")
# Here's a Lorenz system. As it is, it's taking me about 1 second to run
# on my system, with @time.

function lorenz(du,u,p,t)
     du[1] = 10.0*(u[2]-u[1])
     du[2] = u[1]*(28.0-u[3]) - u[2]
     du[3] = u[1]*u[2] - (8/3)*u[3]
     return du
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,1000.0)
prob = ODEProblem(lorenz,u0,tspan)
@time sol = solve(prob)

plot(sol)
png("./techniques-to-increase-speed/plots/lorenz-1")

println("====================")

# Here's a different way to code it without any assignment.

function lorenz(u,p,t)
     return [10.0*(u[2]-u[1]),
             u[1]*(28.0-u[3]) - u[2],
             u[1]*u[2] - (8/3)*u[3]]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,1000.0)
prob = ODEProblem(lorenz,u0,tspan)
@time sol = solve(prob)

plot(sol)
png("./techniques-to-increase-speed/plots/lorenz-2")

# Here, we have NO ASSIGNMENTS WHATSOEVER. And the code is almost 20
# times faster as a result. If you compare the two PNGs side by side,
# you can see that they are the same thing.

# Let's try this same technique with some of the code we wrote just
# earlier today, and see if we can't wring some performance benefits out
# of the fuckers.

########################################################################

# Here's our original KuramotoModel. Very stateful.

function KuramotoModel(N, K, omega)
    @assert(N==size(omega)[1], "number of oscillators != N")

    function Kuramoto(dtheta, theta, p, t)
        for i in 1:N
            theta[i] = mod2pi(theta[i])

            stabilizers = map(theta_j -> sin(theta_j - theta[i]), theta)
            stabilizing_term = (K / N) * sum(stabilizers)

            dtheta[i] = omega[i] + stabilizing_term

            theta[i] = mod2pi(theta[i])
        end
    end

    return Kuramoto
end
# This model has an enormous amount of overhead. For 2 measly
# oscillators, we ended up with 728 MB and over 14 million allocations.
# We can do better. We *must* do better.

# But how? First, let's generate a picture just so we know we aren't
# messing anything up. Fundamentally, this is a refactoring project.

println(":::::::::::::::::::::::::::::::::")

tspan = (0.0, 60.0)

K2 = KuramotoModel(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_1")

# On a virgin run, plotting all that data takes almost 20 seconds, and
# solving takes almost 8. Let's see if we can do better by paring down
# the state stuff.

########################################################################

# First off, let's try to pare out the stabilizers stuff in there, as well as one of the mod2pi assignments.

function KuramotoModel2(N, K, omega)
    @assert(N==size(omega)[1], "number of oscillators != N")

    function Kuramoto(dtheta, theta, p, t)
        for i in 1:N
            dtheta[i] = omega[i] + (K / N) * sum(map(theta_j -> sin(theta_j - theta[i]), theta))
            theta[i] = mod2pi(theta[i])
        end
    end

    return Kuramoto
end

println(":::::::::::::::::::::::::::::::::")

tspan = (0.0, 60.0)

K2 = KuramotoModel2(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_2")

# Any luck? Em... Yea, it seems like we got a slight performance
# increase. From 15 to 9 seconds. Not earth-shattering, but reasonable.

########################################################################

# Okay, let's rerun the whole script virgin and see what happens.

# Oh, wait. Actually, the second one seems to run with slightly *worse*
# time. And a whole lot more garbage collection going on. Well... Let's
# see if we can't keep going, aha.

function KuramotoModel3(N, K, omega)
    @assert(N==size(omega)[1], "number of oscillators != N")
    # Let's try initializing dtheta within the function, instead of
    # having it as an explicit argument.

    function Kuramoto(theta, p, t)
        dtheta = Array{Float64}(undef, N)
        for i in 1:N
            dtheta[i] = omega[i] + (K / N) * sum(map(theta_j -> sin(theta_j - theta[i]), theta))
            theta[i] = mod2pi(theta[i])
        end
        return dtheta
    end

    return Kuramoto
end

println(":::::::::::::::::::::::::::::::::")

tspan = (0.0, 60.0)

K2 = KuramotoModel3(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_3")

#######################################################

# How does our virgin run look now?

# Ah, christ. It looks like we're actually getting even WORSE now as
# time goes on. What gives? Well... At least our plots all look the
# same.

#######################################################

# Let's try moving the dtheta out and see if that makes any difference.

function KuramotoModel4(N, K, omega)
    @assert(N==size(omega)[1], "number of oscillators != N")
    # Let's try initializing dtheta within the function, instead of
    # having it as an explicit argument.

    dtheta = Array{Float64}(undef, N)

    function Kuramoto(theta, p, t)
        for i in 1:N
            dtheta[i] = omega[i] + (K / N) * sum(map(theta_j -> sin(theta_j - theta[i]), theta))
            theta[i] = mod2pi(theta[i])
        end
        return dtheta
    end

    return Kuramoto
end

println(":::::::::::::::::::::::::::::::::")

tspan = (0.0, 60.0)

K2 = KuramotoModel4(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_4")

#######################################################

# I think I might have misunderstood what he meant by "non-allocating".
# Let me try a different approach.

function KuramotoModel5(N, K, omega)
    @assert(N==size(omega)[1], "number of oscillators != N")

    function Kuramoto!(dtheta, theta, p, t)
        for i in 1:N
            dtheta[i] = omega[i] + (K / N) * sum(map(theta_j -> sin(theta_j - theta[i]), theta))
            theta[i] = mod2pi(theta[i])
        end
        return dtheta
    end

    return Kuramoto!
end

println(":::::::::::::::::::::::::::::::::")

tspan = (0.0, 60.0)

K2 = KuramotoModel5(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_5")

#######################################################

# There is a bit of a speed-up! 🙂 Plotting takes
# longer, however. Let's try another virgin run.

# Yes. According to this, the change we made just now puts us
# comfortable back at around the same place we were with KuramotoModel2.
# KuramotoModel the original is still the fastest though, lol.

function KuramotoModel6(N, K, omega)
    @assert(N==size(omega)[1], "number of oscillators != N")

    function Kuramoto(dtheta, theta, p, t)
        for i in 1:N
            theta[i] = mod2pi(theta[i])

            stabilizers = map(theta_j -> sin(theta_j - theta[i]), theta)
            stabilizing_term = (K / N) * sum(stabilizers)

            dtheta[i] = omega[i] + stabilizing_term

            theta[i] = mod2pi(theta[i])
        end
    end

    return Kuramoto
end
# This model has an enormous amount of overhead. For 2 measly
# oscillators, we ended up with 728 MB and over 14 million allocations.
# We can do better. We *must* do better.

# But how? First, let's generate a picture just so we know we aren't
# messing anything up. Fundamentally, this is a refactoring project.

println(":::::::::::::::::::::::::::::::::")

tspan = (0.0, 60.0)

K2 = KuramotoModel6(2, 1.0, [1.0, 1.0])

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_6")
########################################################################

# Let's actually try hard-coding it, and see what happens when we do.

function Kuramoto(du, u, p, t)
    u[1] = mod2pi(u[1])
    u[2] = mod2pi(u[2])
    du[1] = 1.0 + (1.0/2.0) * sin(u[2] - u[1])
    du[2] = 1.0 + (1.0/2.0) * sin(u[1] - u[2])
end

tspan = (0.0, 60.0)

K2 = Kuramoto

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_7")

# Yea, I hate to admit it but this in-place form is the way to fuckin'
# go. OOF.

########################################################################

# Let's move on to the next tip: STACK allocation.

using StaticArrays

function Kuramoto_Static(u, p, t)
    dx = 1.0 + (1.0/2.0) * sin(u[2] - u[1])
    dy = 1.0 + (1.0/2.0) * sin(u[1] - u[2])
    @SVector [dx, dy]
end

tspan = (0.0, 60.0)

K2 = Kuramoto_Static

theta0 = [0.0, pi]
prob = ODEProblem(K2, theta0, tspan)

@time sol = solve(prob, Euler(), dt=0.00001) # Yes, Euler() is important.
@time plot(sol)

png("./techniques-to-increase-speed/plots/naive-2-osc-solution_8")


text:
module ForsakenCitySatellite

using ..Ahorn, Maple

function satelliteFinalizer(entity::Maple.ForsakenCitySatellite)
    x, y = Ahorn.position(entity)

    entity.data["nodes"] = [
        (x + 64, y),
        (x + 48, y),
    ]
end

const placements = Ahorn.PlacementDict(
    "Forsaken City Satellite" => Ahorn.EntityPlacement(
        Maple.ForsakenCitySatellite,
        "point",
        Dict{String, Any}(),
        satelliteFinalizer
    )
)

Ahorn.nodeLimits(entity::Maple.ForsakenCitySatellite) = 2, 2

birdSprite = "scenery/flutterbird/idle00.png"
gemSprite = "collectables/heartGem/0/00.png"

dishSprite = "objects/citysatellite/dish.png"
lightSprite = "objects/citysatellite/light.png"
computerSprite = "objects/citysatellite/computer.png"
computerScreenSprite = "objects/citysatellite/computerscreen.png"

# Offsets from game
computerOffsetX = 32
computerOffsetY = 24

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ForsakenCitySatellite)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ((x, y), (x, y)))
    birdX, birdY = nodes[1]
    gemX, gemY = nodes[2]

    Ahorn.drawSprite(ctx, gemSprite, gemX, gemY)
    Ahorn.drawSprite(ctx, birdSprite, birdX, birdY, jx=0.5, jy=1.0)

    for node in nodes
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, x, y, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

function Ahorn.selection(entity::Maple.ForsakenCitySatellite)
    x, y = Ahorn.position(entity)

    nodes = get(entity.data, "nodes", ((x, y), (x, y)))
    birdX, birdY = nodes[1]
    gemX, gemY = nodes[2]

    return [
        Ahorn.coverRectangles([
            Ahorn.getSpriteRectangle(dishSprite, x, y, jx=0.5, jy=1.0),
            Ahorn.getSpriteRectangle(lightSprite, x, y, jx=0.5, jy=1.0),
            Ahorn.getSpriteRectangle(computerSprite, x + computerOffsetX, y + computerOffsetY),
            Ahorn.getSpriteRectangle(computerScreenSprite, x + computerOffsetX, y + computerOffsetY)
        ]),
        Ahorn.getSpriteRectangle(birdSprite, birdX, birdY, jx=0.5, jy=1.0),
        Ahorn.getSpriteRectangle(gemSprite, gemX, gemY)
    ]
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ForsakenCitySatellite, room::Maple.Room)
    x, y = Ahorn.position(entity)

    Ahorn.drawSprite(ctx, dishSprite, x, y, jx=0.5, jy=1.0)
    Ahorn.drawSprite(ctx, lightSprite, x, y, jx=0.5, jy=1.0)
    Ahorn.drawSprite(ctx, computerSprite, x + computerOffsetX, y + computerOffsetY)
    Ahorn.drawSprite(ctx, computerScreenSprite, x + computerOffsetX, y + computerOffsetY)
end

end

ext: jl


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct ScenariosMetadata <: ForecastMetadata
        name::String
        resolution::Dates.Period
        initial_timestamp::Dates.DateTime
        interval::Dates.Period
        scenario_count::Int64
        count::Int
        time_series_uuid::UUIDs.UUID
        horizon::Int
        scaling_factor_multiplier::Union{Nothing, Function}
        internal::InfrastructureSystemsInternal
    end

A Discrete Scenario Based time series for a particular data field in a Component.

# Arguments
- `name::String`: user-defined name
- `resolution::Dates.Period`
- `initial_timestamp::Dates.DateTime`: time series availability time
- `interval::Dates.Period`: time step between forecast windows
- `scenario_count::Int64`: Number of scenarios
- `count::Int`: number of forecast windows
- `time_series_uuid::UUIDs.UUID`: reference to time series data
- `horizon::Int`: length of this time series
- `scaling_factor_multiplier::Union{Nothing, Function}`: Applicable when the time series data are scaling factors. Called on the associated component to convert the values.
- `internal::InfrastructureSystemsInternal`
"""
mutable struct ScenariosMetadata <: ForecastMetadata
    "user-defined name"
    name::String
    resolution::Dates.Period
    "time series availability time"
    initial_timestamp::Dates.DateTime
    "time step between forecast windows"
    interval::Dates.Period
    "Number of scenarios"
    scenario_count::Int64
    "number of forecast windows"
    count::Int
    "reference to time series data"
    time_series_uuid::UUIDs.UUID
    "length of this time series"
    horizon::Int
    "Applicable when the time series data are scaling factors. Called on the associated component to convert the values."
    scaling_factor_multiplier::Union{Nothing, Function}
    internal::InfrastructureSystemsInternal
end

function ScenariosMetadata(name, resolution, initial_timestamp, interval, scenario_count, count, time_series_uuid, horizon, scaling_factor_multiplier=nothing, )
    ScenariosMetadata(name, resolution, initial_timestamp, interval, scenario_count, count, time_series_uuid, horizon, scaling_factor_multiplier, InfrastructureSystemsInternal(), )
end

function ScenariosMetadata(; name, resolution, initial_timestamp, interval, scenario_count, count, time_series_uuid, horizon, scaling_factor_multiplier=nothing, internal=InfrastructureSystemsInternal(), )
    ScenariosMetadata(name, resolution, initial_timestamp, interval, scenario_count, count, time_series_uuid, horizon, scaling_factor_multiplier, internal, )
end

"""Get [`ScenariosMetadata`](@ref) `name`."""
get_name(value::ScenariosMetadata) = value.name
"""Get [`ScenariosMetadata`](@ref) `resolution`."""
get_resolution(value::ScenariosMetadata) = value.resolution
"""Get [`ScenariosMetadata`](@ref) `initial_timestamp`."""
get_initial_timestamp(value::ScenariosMetadata) = value.initial_timestamp
"""Get [`ScenariosMetadata`](@ref) `interval`."""
get_interval(value::ScenariosMetadata) = value.interval
"""Get [`ScenariosMetadata`](@ref) `scenario_count`."""
get_scenario_count(value::ScenariosMetadata) = value.scenario_count
"""Get [`ScenariosMetadata`](@ref) `count`."""
get_count(value::ScenariosMetadata) = value.count
"""Get [`ScenariosMetadata`](@ref) `time_series_uuid`."""
get_time_series_uuid(value::ScenariosMetadata) = value.time_series_uuid
"""Get [`ScenariosMetadata`](@ref) `horizon`."""
get_horizon(value::ScenariosMetadata) = value.horizon
"""Get [`ScenariosMetadata`](@ref) `scaling_factor_multiplier`."""
get_scaling_factor_multiplier(value::ScenariosMetadata) = value.scaling_factor_multiplier
"""Get [`ScenariosMetadata`](@ref) `internal`."""
get_internal(value::ScenariosMetadata) = value.internal

"""Set [`ScenariosMetadata`](@ref) `name`."""
set_name!(value::ScenariosMetadata, val) = value.name = val
"""Set [`ScenariosMetadata`](@ref) `resolution`."""
set_resolution!(value::ScenariosMetadata, val) = value.resolution = val
"""Set [`ScenariosMetadata`](@ref) `initial_timestamp`."""
set_initial_timestamp!(value::ScenariosMetadata, val) = value.initial_timestamp = val
"""Set [`ScenariosMetadata`](@ref) `interval`."""
set_interval!(value::ScenariosMetadata, val) = value.interval = val
"""Set [`ScenariosMetadata`](@ref) `scenario_count`."""
set_scenario_count!(value::ScenariosMetadata, val) = value.scenario_count = val
"""Set [`ScenariosMetadata`](@ref) `count`."""
set_count!(value::ScenariosMetadata, val) = value.count = val
"""Set [`ScenariosMetadata`](@ref) `time_series_uuid`."""
set_time_series_uuid!(value::ScenariosMetadata, val) = value.time_series_uuid = val
"""Set [`ScenariosMetadata`](@ref) `horizon`."""
set_horizon!(value::ScenariosMetadata, val) = value.horizon = val
"""Set [`ScenariosMetadata`](@ref) `scaling_factor_multiplier`."""
set_scaling_factor_multiplier!(value::ScenariosMetadata, val) = value.scaling_factor_multiplier = val
"""Set [`ScenariosMetadata`](@ref) `internal`."""
set_internal!(value::ScenariosMetadata, val) = value.internal = val


text:
using Arpack, SparseArrays, LinearAlgebra

#include("file.jl")

# defining vectors
n = 10
a = Matrix(I, n, n)

b = [1 0 0; 0 1 0; 0 0 1]
c = zeros(10,1)
d = ones(10,1)

e = sparse(a)
eigs(e)

# defining functions
function f(x)
    return x + 1
end

# g(1) will fail, but g(1.0) will work
function g(x::Float64)
    return x+1
end

# h(1) will work, but h("1") will not
function h(x::T) where {T <: Number}
    return x+1
end

function main(x)
    f(x)
    #g(x)
    h(x)

    return nothing
end

main(1)

@time main(1)

using BenchmarkTools

@benchmark main(1)

using Plots

a = rand(10,10)

heatmap(a)

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
module AdventOfCode

using HTTP, Dates, JSON
using Mocking
export setup_files

_base_url(year, day) = "https://adventofcode.com/$year/day/$day"

function _get_cookies()
    if "AOC_SESSION" ∉ keys(ENV)
        error("Session cookie in ENV[\"AOC_SESSION\"] needed to download data.")
    end
    return Dict("session" => ENV["AOC_SESSION"])
end

function _download_data(year, day)
    result = @mock HTTP.get(_base_url(year, day) * "/input", cookies = _get_cookies())
    if result.status == 200
        return result.body
    end
    error("Unable to download data")
end

function _template(year, day; include_year = true)
    rel_path = joinpath("data", "day_$day.txt")
    include_year && (rel_path = joinpath(string(year), rel_path))

    """
    # $(_base_url(year, day))
    using AdventOfCode

    input = readlines("$rel_path")

    function part_1(input)
        nothing
    end
    @info part_1(input)

    function part_2(input)
        nothing
    end
    @info part_2(input)
    """
end

function _setup_data_file(year, day; include_year = true)
    rel_path = joinpath("data", "day_$day.txt")
    include_year && (rel_path = joinpath(string(year), rel_path))

    data_path = joinpath(pwd(), rel_path)
    if isfile(data_path)
        @warn "$data_path already exists. AdventOfCode.jl will not redownload it"
        return nothing
    end
    time_req = HTTP.get("http://worldclockapi.com/api/json/est/now")
    current_datetime = JSON.parse(String(time_req.body))["currentDateTime"]
    current_date = Date(current_datetime[1:10])
    if _is_unlocked(year, day)
        data = _download_data(year, day)
        mkpath(splitdir(data_path)[1])
        open(data_path, "w+") do io
            write(io, data)
        end
    end
end

function _is_unlocked(year, day)
    time_req = HTTP.get("http://worldclockapi.com/api/json/est/now")
    current_datetime = JSON.parse(String(time_req.body))["currentDateTime"]
    current_date = Date(current_datetime[1:10])
    is_unlocked = current_date >= Date(year, 12, day)
    if !is_unlocked
        @warn "Advent of Code for year $year and day $day hasn't unlocked yet."
    end
    is_unlocked
end

"""
    setup_files(year, day; force = false, include_year = true)
    setup_files(; force = false, include_year = true)

Downloads the input file for the specified year and date and stores that file in
`{year}/data/day_{day}.txt` or `data/day_{day}.txt` if `include_year` is set to `false`.
Also sets up a template file in `{year}/src/day_{day}.jl` for your script
(`src/day_{day}.jl` if `include_year` is set to `false`).
`force=true` will recreate the `src` file even if it already exists.
The `data` file will not be re-downloaded even with `force=true`
since it's a static file and to reduce load on AdventOfCode's servers.

If `year` and `day` are not provided, the setup defaults to today's date.
"""
function setup_files(year, day; force = false, include_year = true)
    is_unlocked = _is_unlocked(year, day)
    rel_path = joinpath("src", "day_$day.jl")
    include_year && (rel_path = joinpath(string(year), rel_path))
    code_path = joinpath(pwd(), rel_path)
    is_unlocked &&  _setup_data_file(year, day, include_year = include_year)
    if !force && isfile(code_path)
        @warn "$code_path already exists. To overwrite, re-run with `force=true`"
    else
        mkpath(splitdir(code_path)[1])
        open(code_path, "w+") do io
            write(io, _template(year, day, include_year=include_year))
        end
    end
    return code_path
end

setup_files(; force = false, include_year = true) = setup_files(year(today()), day(today()), force = force, include_year = include_year)

end


text:
using Test
using Main.CGP
@assert MyInt == UInt16   # Must be set to MyInt16 in CGP.jl

# Does not test comp_experiment(), rand_comp_circuit(), rand_gate_circuit()
@testset "compositional" begin  # @testset makes random_chromosome(p) always return same chromosome

p = Parameters( 2, 1, 3, 3 )
funcs = default_funcs(p.numinputs)
ctx = construct_context(p.numinputs)

cg3 = CGP.Cg( NAND, [1,2], 1 )
cg4 = CGP.Cg( XOR, [2,3], 1 )
cg5 = CGP.Cg( AND, [4,3], 1 )
cg6 = CGP.Cg( NOR, [3,5], 1 )
cg7 = CGP.Cg( NOR, [5,4], 1 )
cc7 = CGP.CC( CGP.CC[cg3, cg4, cg5, cg6, cg7], [1,2], 1)
ch7 = circuit((1,2), ((3,NAND,1,2), (4,XOR,2,3), (5,AND,4,3), (6,NOR,3,5), (7,NOR,5,4)))
recurse_step = [0]
@test execute!(cg3,ctx,recurse_step,active_only=true) == execute!(cg3,ctx,recurse_step,active_only=false) == 0x0007
recurse_step = [0]
@test execute!(cg4,[0x000c, 0x000a, 0x0007],recurse_step,active_only=true) == execute!(cg4,[0x000c, 0x000a, 0x0007],recurse_step,active_only=false) == 0x000d
recurse_step = [0]
@test execute!(cg5,[0x000c, 0x000a, 0x0007, 0x000d],recurse_step,active_only=true) == execute!(cg5,[0x000c, 0x000a, 0x0007, 0x000d], recurse_step,active_only=false) == 0x0005
recurse_step = [0]
@test execute!(cg6,[0x000c, 0x000a, 0x0007, 0x000d, 0x0005],recurse_step,active_only=true) == execute!(cg6,[0x000c, 0x000a, 0x0007, 0x000d, 0x0005],recurse_step,active_only=false) == 0x0008
recurse_step = [0]
@test execute!(cg7,[0x000c, 0x000a, 0x0007, 0x000d, 0x0005, 0x0008],recurse_step,active_only=true) == execute!(cg7,[0x000c, 0x000a, 0x0007, 0x000d, 0x0005, 0x0008],recurse_step,active_only=false) == 0x0002
recurse_step = [0]
@test output_values(ch7)[1] == execute!(cc7,ctx,recurse_step) == execute!(cg7,[0x000c, 0x000a, 0x0007, 0x000d, 0x0005, 0x0008],recurse_step,active_only=true)

cc1 = CGP.CC( CGP.CC[cg3, cg4], [1,2], 1 )
recurse_step = [0]
@test execute!(cc1,ctx,recurse_step,active_only=true) == execute!(cc1,ctx,recurse_step,active_only=false) == 0x000d
cc4 = CGP.CC( CGP.CC[cg4, cc1], [2,1], 1 )
recurse_step = [0]
@test execute!(cc4,ctx,recurse_step,active_only=true) == execute!(cc4,ctx,recurse_step,active_only=false) == 0x000b
cc6 = CGP.CC( CGP.CC[cc1,cc4], [1,2], 1 )
recurse_step = [0]
@test execute!(cc6,ctx,recurse_step,active_only=true) == execute!(cc6,ctx,recurse_step,active_only=false) == 0x000b

numinputs = 3
p = Parameters(numinputs,1,4,3)
ctx = construct_context(numinputs)
ccg = CGP.Cg(AND,[2,3],1)
recurse_step = [0]
@test execute!(ccg,ctx,recurse_step) == 0x0088
recurse_step = [0]
cc2 = CGP.CC(CGP.CC[CGP.Cg(AND,[2,3],1),CGP.Cg(XOR,[1,4],1)],[1,2,3],1)
@test execute!(cc2,ctx,recurse_step) == 0x0078

ch = random_chromosome(p)
print_circuit(ch)
println("output_values(ch)[1]: ",@sprintf("0x00%x",output_values(ch)[1]))
println("chromosome_to_circuit(ch): ",chromosome_to_circuit(ch))
recurse_step = [0]
@test execute!(chromosome_to_circuit(ch),ctx,recurse_step) == output_values(ch)[1]

cg4 = CGP.Cg( NAND, [1,2], 1 ) #4
cg5 = CGP.Cg( XOR, [4,3], 1 )  #5
cg6 = CGP.Cg( AND, [4,3], 1 )  #6
cg7 = CGP.Cg( NOR, [3,5], 1 )  #7
cg8 = CGP.Cg( OR, [5,6] ,1 )  #8
cc1 = CGP.CC( CGP.CC[cg4, cg5 ], [2,1,3], 1 )
cc2 = CGP.CC( CGP.CC[cg5, cg4, cg7], [2,3,1], 1 )
cc3 = CGP.CC( CGP.CC[cg4, cg5, cg6, cg7, cg8], [3,1,2], 1 )
cc4 = CGP.CC( CGP.CC[cc1], [1,3,2], 1 )
cc5 = CGP.CC( CGP.CC[cc4, cc2], [1,3,2], 1 )
recurse_step = [0]
@test execute!(cg4,ctx,recurse_step) == 0x003f
recurse_step = [0]
@test execute!(cc1,ctx,recurse_step) == 0x0095 == execute_chromosome(circuit_to_chromosome(cc1),ctx)[1]
recurse_step = [0]
@test execute!(cc2,ctx,recurse_step) == 0x0008 == execute_chromosome(circuit_to_chromosome(cc2),ctx)[1]
recurse_step = [0]
@test execute!(cc3,ctx,recurse_step) == 0x00df == execute_chromosome(circuit_to_chromosome(cc3),ctx)[1]
recurse_step = [0]
@test execute!(cc4,ctx,recurse_step) == 0x0093 == execute!(cc1,ctx[[1,3,2]],recurse_step)
recurse_step = [0]
@test execute!(cc5,ctx,recurse_step) == 0x0008

recurse_step = [0]
@test execute_chromosome(ch,ctx)[1] == execute!(chromosome_to_circuit(ch),ctx,recurse_step)
recurse_step = [0]
@test execute_chromosome(circuit_to_chromosome(cc1),ctx)[1] == execute!(cc1,ctx,recurse_step)
cc9 = random_gate_circuit( 3, 4, funcs)
ch9 = circuit_to_chromosome(cc9); print_circuit(ch9)
recurse_step = [0]
@test execute_chromosome(ch9,ctx)[1] == execute!(chromosome_to_circuit(ch9),ctx,recurse_step)
recurse_step = [0]
@test execute_chromosome(circuit_to_chromosome(cc9),ctx,permute_context=true)[1] == execute!(cc9,ctx,recurse_step)
recurse_step = [0]
@test execute_chromosome(circuit_to_chromosome(cc9),ctx)[1] == execute!(cc9,ctx,recurse_step)

Random.seed!(1);
needs_list = [([0x000a],2)]
numinputs=2
numcompositions=2
numcircuits=4
funcs = [AND,OR,XOR,NAND,NOR,IN1,IN2]
(needs,ac) =comp_experiment( numinputs, numcompositions, needs_list, numcircuits, funcs );
@test needs[[0x000a]].goal[1] == 0x000a  # Might not work if random number generation changes

end  # @testset


text:
## This code was suggested by Doug Bates, August 20, 2018
function compute_maxlam(xty::AbstractVector{<:Real},
                         penalty_factor::AbstractVector{<:Real})
    # compute largest value for lambda.
    # we need to divide by penalty_factor (making sure
    # not to divide by zero) to account for
    # differential penalization
    lambda_max = maximum(x / (iszero(pf) ? one(pf) : pf) for (x, pf) in zip(xty, penalty_factor))
end

function make_steps_A(xtx::AbstractMatrix{T}) where T <: BLAS.BlasFloat
    pp, p  = size(xtx)

    # compute largest eigenvalue
    d, v   = eigs(xtx, nev = 1, which = :LR)
    d      = reinterpret(Float64, d)[1]

    # make a vector of stepsizes based on d
    steps  = repeat([d], outer = p)

    # Calculate A from oem algorithm
    A      = Diagonal(steps) - xtx
    steps, A, d
end

function compute_xtx_y(x::AbstractMatrix{<:Real},
                       y::AbstractVector{<:Real},
                       intercept::Bool = true,
                       standardize::Bool = true)

    n, p = size(x)

    intscale = 1.0

    if standardize
        colsquare = sqrt.(sum(x .^ 2, dims = 1) / (n - 1))
        colsums   = sum(x, dims = 1) / n
    else
        colsquare = ones(p)
        colsums   = ones(p)
    end

    colsquare_inv = 1 ./ colsquare

    if intercept & standardize
        xtx = zeros(p + 1, p + 1)
        xty = zeros(p + 1)
        intscale = 1 / n

        colsums .*= colsquare_inv

        xtx[2:(p+1),2:(p+1)] = Diagonal(colsquare_inv[1,:]) * Symmetric(x'x) * Diagonal(colsquare_inv[1,:]) / n
        xtx[1, 2:(p+1)] = colsums * intscale
        xtx[2:(p+1), 1] = xtx[1, 2:(p+1)]'
        xtx[1, 1] = 1.0 * intscale
        xty[2:(p+1)] = ((x' * y) .* colsquare_inv') / n'
        xty[1] = mean(y) * intscale

    elseif intercept
        xtx = zeros(p + 1, p + 1)
        xty = zeros(p + 1)

        xtx[2:(p+1), 2:(p+1)] = Symmetric(x'x) / n
        xtx[1, 2:(p+1)] = colsums
        xtx[2:(p+1), 1] = xtx[1, 2:(p+1)]'
        xtx[1, 1] = 1.0

        xty[2:(p+1)] = x' * y / n'
        xty[1] = mean(y)
    elseif standardize
        xtx = Diagonal(colsquare_inv[1,:]) * Symmetric(x'x) * Diagonal(colsquare_inv[1,:]) / n
        xty = (((x' * y) .* colsquare_inv') / n')[:,1]
    else
        xtx = Symmetric(x'x) / n
        xty = x' * y / n'
    end
    Symmetric(xtx), xty, colsquare_inv, intscale
end

function create_steps(x::StridedMatrix{T},
                      group_col::StridedVector{Int}) where T <: BLAS.BlasFloat
    grps   = unique(group_col)
    n_grps = length(grps)
    n, p   = size(x)

    steps  = ones(p)
    p_cur = 0

    xtx_list = Any[]
    grp_idx_list = Any[]
    grp_bool_list = Any[]
    steps_list = Any[]
    eigens = zeros(n_grps)
    for j in 1:n_grps
        cur_grp = group_col .== grps[j]
        grp_len = sum(cur_grp)
        xc      = x[:,cur_grp]
        xtx_cur = Symmetric(xc'xc) / n

        d, v   = eigs(xtx_cur, nev = 1)

        push!(steps_list, repeat(d, outer = grp_len))
        push!(xtx_list, xtx_cur)
        push!(grp_bool_list, cur_grp)
        push!(grp_idx_list, findall(cur_grp))
        eigens[j] = d[1]
    end
    eigens, steps_list, xtx_list, grp_idx_list, grp_bool_list, grps, n_grps
end


text:
module Fridge

include("recipeWebscraper.jl")

using .recipeWebscraper, JLD2

export checkIngredients, greedyFindCombo, findBestRecipe, randomCombo, Neighbour, SAFindCombo, scrapeRecipe, loadRecipeDBCSV

#==================================================
          CHECK INGREDIENTS FUNCTIONS
==================================================#

function createIngredientDatabase(recipeDict)
    # create an unique vector of all ingredients used in the recipe database
    ingredientList = []
    for ingredients in values(recipeDict)
        append!(ingredientList,ingredients)
    end
    ingredientList = unique(ingredientList)

    return ingredientList
end

"""
    checkIngredients(fridgeList, ingredientList)

This function checks if the foods in your fridge are also found in the ingredient overview
of the recipe database. If they are not found in the database, regex is used to find possible alternatives.
For instance cheese may be replaced by swiss cheese.

## Input:
- fridgeList: A list containing the different foods in your fridge as a string.
- ingredientList: A list containing all the different ingredients that are used in the recipe database.
- test: An optional Boolean argument. When true the function simulates user input.

## Output:
- fridgeList: The (adapted) given fridgeList

"""
function checkIngredients(fridgeList,ingredientList; test=false)

    print("Checking if the food in your fridge is found in our database.\n\n")

    # check if the ingredients in the fridge are found in the database
    for food in fridgeList
        if food in ingredientList
            print("Found $food in the ingredient database.\n")
        else
            print("Did not find $food in the ingredient database.\n\n")
            if any(occursin.(food,ingredientList))
                print("Possible alternatives in database:\n")
                alternatives = ingredientList[occursin.(food,ingredientList)]
                for (indexNum, alternative) in enumerate(alternatives)
                    print("[$indexNum] $alternative\n")
                end
                print("If you want to take an alternative type its number, else type no.\n")
                if !test
                    # read wanted alternative number
                    answer = readline()
                else
                    # simulate user input
                    answer = string(rand(1:length(alternatives)))
                    print("The computer typed $answer.")
                end
                if answer != "no"
                    correctInput = false
                    while !correctInput
                        try
                            answer = parse(Int64,answer)
                            correctInput = true
                        catch
                            print("Please only type the number of the alternative, eg. 1\n")
                            answer = readline()
                        end
                    end
                    fridgeList[fridgeList .== food] .= alternatives[answer]
                    print("Replaced $food with $(alternatives[answer])\n")
                end
            else
                print("Did not find any alternative in the ingredient database.\n")
            end
        end
    end
    print("Done checking ingredients.\n\n")

    return fridgeList
end

#==================================================
            OBJECTIVE FUNCTIONS
==================================================#

# objective function used to score a single recipe
fridgeObjective(x::Array{Int64}) = sum(x[1:end-1]) + 6*sum(x[1:end-1] .== 0) + 2*x[end]

# objective function used to score a combination of recipes
fridgeObjective(x) = compatible(x) ? sum(sum(x) .== 0)*6 + sum(x)[end]*2 : Inf

#==================================================
                GREEDY ALGORITHM
==================================================#

"""
    greedyFindCombo(fridgeList, recipeDict, numRecipes)

This function uses greedy search to find a good combination of recipes that match your fridge content.
It ranks all recipes based on the following formula

``score = (ingredients from fridge used) + 6*(ingredients in fridge remaining) + 2*(extra ingredients needed)``

## Input:
- fridgeList: A list containing the different foods in your fridge as a string.
- recipeDict: A dictionary in which the keys are the recipe names and the responding values are a list of the needed ingredients.
- numRecipes: The max amount of recipes that a combo should contain.

## Output:
- bestCombo: A dictionary containing the best found combination of recipes.

"""
function greedyFindCombo(fridgeList, recipeDict, numRecipes)

    bestCombo = Dict()
    ingredientsArray = []
    namesArray = []

    for (name,ingredients) in recipeDict
        push!(namesArray,name)
        push!(ingredientsArray, recipeToNumVector(fridgeList, ingredients))
    end

    for i = 1:numRecipes

        bestOrder = sortperm(ingredientsArray, by=fridgeObjective)
        ingredientsArray = ingredientsArray[bestOrder]
        namesArray = namesArray[bestOrder]

        # break if all recipes that are left, don't use anything from the fridge
        if all(ingredientsArray[1][1:end-1] .== 0)
            break
        end

        tempRecipeName = namesArray[1]
        bestCombo[tempRecipeName] = ingredientsArray[1]

        # break if all ingredients are used
        if all(isone.(sum(values(bestCombo))))
            break
        end

        namesArray = [name for (name, ingredientList) in zip(namesArray, ingredientsArray) if sum(ingredientList[1:end-1] .& bestCombo[tempRecipeName][1:end-1] ) == 0]
        ingredientsArray = [ingredientList for ingredientList in ingredientsArray if sum(ingredientList[1:end-1] .& bestCombo[tempRecipeName][1:end-1] ) == 0]

        # break if there are no recipes left
        if isempty(namesArray)
            break
        end

    end

    return bestCombo
end

#==================================================
                NEIGBOURHOODS
==================================================#

"""
    randomCombo(fridgeList, recipeDict, numRecipes)

This function gives a random combination of recipes from the provided recipe dictionary.

## Input:
- fridgeList: A list containing the different foods in your fridge as a string.
- recipeDict: A dictionary in which the keys are the recipe names and the responding values are a list of the needed ingredients.
- numRecipes: The max amount of recipes that a combo should contain.

## Output:
- randCombo: A dictionary containing a random combination of recipes.
"""
function randomCombo(fridgeList, recipeDict, numRecipes)

    randCombo = Dict()
    ingredientsArray = []
    namesArray = []

    for (name,ingredients) in recipeDict
        push!(namesArray,name)
        push!(ingredientsArray, recipeToNumVector(fridgeList, ingredients))
    end

    for i = 1:numRecipes

        randIndex = rand(1:length(namesArray))
        tempRecipeName = namesArray[randIndex]
        randCombo[tempRecipeName] = ingredientsArray[randIndex]

        # break if all ingredients are used
        if all(isone.(sum(values(randCombo))))
            break
        end

        namesArray = [name for (name, ingredientList) in zip(namesArray, ingredientsArray) if sum(ingredientList[1:end-1] .& randCombo[tempRecipeName][1:end-1] ) == 0]
        ingredientsArray = [ingredientList for ingredientList in ingredientsArray if sum(ingredientList[1:end-1] .& randCombo[tempRecipeName][1:end-1] ) == 0]

        # break if there are no recipes left
        if isempty(namesArray)
            break
        end

    end

    return randCombo
end

"""
    Neighbour(curSolution, fridgeList, recipeDict, numRecipes, tabuList, randRecipe)

This is a function that looks for a neighbour of the current solution. This function is used in the simulated annealing algorithm.
It also uses a tabulist to stimulate the use of new solutions.

## Input:
- curSolution: The current best combination of recipes, given as a dictionary in which the keys are the recipe names and the responding values are a list of the needed ingredients.
- fridgeList: A list containing the different foods in your fridge as a string.
- recipeDict: A dictionary in which the keys are the recipe names and the responding values are a list of the needed ingredients.
- numRecipes: The max amount of recipes that a combo should contain.
- tabuList: A list of recipes that should not be used in the found neighbour.
- randRecipe: A Boolean `true` or `false` value. When `true`, random recipes are used for the neighbour.

## Output:
- neighbour: A dictionary containing a combination of recipes.

"""
function Neighbour(curSolution, fridgeList, recipeDict, numRecipes, tabuList, randRecipe)

    toRemove = rand(curSolution)[1]
    # adapt the fridgeList so that only ingredients from the removed ingredient are available
    tempFridgeList = copy(fridgeList)
    for recipe in keys(curSolution)
        if recipe != toRemove
            tempFridgeList = [i for i in fridgeList if !in(i,recipeDict[recipe])]
            numRecipes -= 1
        end
    end
    # adapt the recipeDict and use greedy search to find a new solution
    tempRecipeDict = copy(recipeDict)
    for recipe in keys(curSolution)
        delete!(tempRecipeDict,recipe)
    end

    for recipe in tabuList
        try delete!(tempRecipeDict,recipe)
        catch e
        end
    end

    if randRecipe
        neighbour = randomCombo(tempFridgeList, tempRecipeDict, numRecipes)
    else
        neighbour = greedyFindCombo(tempFridgeList, tempRecipeDict, numRecipes)
    end

    # correct recipe vectors
    for recipe in keys(neighbour)
        neighbour[recipe] = recipeToNumVector(fridgeList,recipeDict[recipe])
    end

    # here combine the two dictionaries
    tempCurSolution = copy(curSolution)
    delete!(tempCurSolution,toRemove)
    neighbour = merge(neighbour,tempCurSolution)

    return neighbour
end

#==================================================
        SIMULATED ANNEALING ALGORITHM
==================================================#

"""
    SAFindCombo(curSolution,  fridgeList, recipeDict, numRecipes, randRecipe; kT=100, r=0.75, Tmax=4, Tmin=1, tabuLength=3)

This function uses simulated annealing to find a better combination of recipes that match your fridge content.
It starts with the current solution and tries to improve this.

## Input:
- curSolution: The current best combination of recipes, given as a dictionary in which the keys are the recipe names and the responding values are a list of the needed ingredients.
- fridgeList: A list containing the different foods in your fridge as a string.
- recipeDict: A dictionary in which the keys are the recipe names and the responding values are a list of the needed ingredients.
- numRecipes: The max amount of recipes that a combo should contain.
- tabuList: A list of recipes that should not be used in the found neighbour.
- randRecipe: A Boolean `true` or `false` value. When `true`, random recipes are used for the neighbour.

## Optional Inputs:
- kT: repetitions per temperature
- r: cooling rate
- Tmax: maximal temperature to start
- Tmin: minimal temperature to end
- tabuLength: Number of cycli that recipe needs to be blocked

## Output:
- solution: A dictionary containing the best found combination of recipes.

"""
function SAFindCombo(curSolution,  fridgeList, recipeDict, numRecipes, randRecipe;
    kT=100, # repetitions per temperature
    r=0.75, # cooling rate
    Tmax=4, # maximal temperature to start
    Tmin=1, # minimal temperature to end
    tabuLength=3) # number of cycli that recipe needs to be blocked

    @assert 0 < Tmin < Tmax "Temperatures should be positive"
	@assert 0 < r < 1 "cooling rate is between 0 and 1"

	solution = curSolution
    if isempty(solution)
        print("There are no recipes in the database matching your fridge.")
        T = Tmin-1
    else
        obj = fridgeObjective([i for i in values(solution)])
	    T = Tmax
    end

    tabuList = String[i for i in keys(curSolution)]

	while T > Tmin
        if isempty(solution)
            print("There are no recipes in the database matching your fridge.")
            break
        end
        print("T = $T \n")
		# repeat kT times
		for i in 1:kT
			sn = Neighbour(solution, fridgeList, recipeDict, numRecipes, tabuList, randRecipe)  # random neighbor
			obj_sn = fridgeObjective([i for i in values(sn)])
			# if the neighbor improves the solution, keep it
			# otherwise accept with a probability determined by the
			# Metropolis heuristic
			if obj_sn < obj || rand() < exp(-(obj_sn-obj)/T)
				solution = sn
				obj = obj_sn
			end

            for recipe in keys(sn)
                if !in(recipe, tabuList)
                    if length(tabuList) < tabuLength
                        pushfirst!(tabuList,recipe)
                    else
                        pop!(tabuList)
                        pushfirst!(tabuList,recipe)
                    end
                end
            end
		end

		# decay temperature
		T *= r
	end

    return solution
end

#==================================================
                OVERVIEW FUNCTION
==================================================#

"""
    findBestRecipe(fridgeList, csvPath; numRecipes=3, randRecipe=false)

This function combines all other functions. This function checks if your ingredients are in the database,
if not it offers possible alternatives. Next it uses simulated annealing to find a better recipe combination.

## Input:
- fridgeList: A list containing the different foods in your fridge as a string.
- dataPath: A relative or absolute path to a .csv or .jld2 file containing the recipe database.

## Optional Inputs:
- numRecipes: The max amount of recipes that a combo should contain.
- randRecipe: A Boolean `true` or `false` value. When `true`, random recipes are used to find the neighbour in simulated annealing.
- testmode: An optional Boolean argument. When true the function simulates user input.

## Output:
- SASolution: A dictionary containing the best found combination of recipes.

"""
function findBestRecipe(fridgeList, dataPath; numRecipes=3, randRecipe=false, testmode=false)

    # load the recipe dictionary from the db file
    recipeDict = dataPath[end-3:end] == ".csv" ? loadRecipeDBCSV(dataPath) : load(dataPath)

    # create a list of all ingredients in your database
    ingredientList = createIngredientDatabase(recipeDict)

    # check for every food in your fridge if it's in the database. If not check if their are alternatives.
    fridgeList = checkIngredients(fridgeList, ingredientList, test=testmode)

    # find the best greedy recipe
    greedySolution = greedyFindCombo(fridgeList, recipeDict, numRecipes)
    print("greedySolution = $greedySolution\n")

    # find the best recipe with SA
    SASolution = SAFindCombo(greedySolution,  fridgeList, recipeDict, numRecipes, randRecipe)

    # print the solution
    print("The best recipes to make are:\n\n")
    for recipeName in keys(SASolution)
        print("$(recipeName) : $(recipeDict[recipeName])\n")
    end

    return SASolution
end

#==================================================
            SUPPORTING FUNCTIONS
==================================================#

"""
    recipeToNumVector(fridgeList,ingredientList)

This function changes an ingredients list of a recipe to a numeric vector. This vector has a length equal
to the number of foods in your fridge plus one. If a certain food from the fridge is used in the recipe, that position in the vector is a 1.
If not it is a 0. The last index of the vector contains the amount of extra ingredients needed.

## Input:
- fridgeList: A list containing the different foods in your fridge as a string.
- ingredientList: A list containing all the ingredients used in a specific recipe.

## Output:
- numVector: This vector has a length equal to the number of foods in your fridge plus one. If a certain food from the fridge is used in the recipe,
that position in the vector is a 1. If not it is a 0. The last index of the vector contains the amount of extra ingredients needed.

"""
function recipeToNumVector(fridgeList,ingredientList)

    numVector = zeros(Int64,length(fridgeList)+1)
    for i in 1:length(fridgeList)
        numVector[i] = fridgeList[i] in ingredientList ? 1 : 0
    end
    numVector[end] = length(ingredientList) - sum(numVector)
    return numVector
end

compatible(x) = !any(sum(x)[1:end-1] .>= 2) # checks if two recipes use a same ingredient

end # module


text:
######################## Simulation Definations #######################
function create_uc_template()
    service = Dict(
        :ReserveUp =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveUp}, PSI.RangeReserve),
        :ReserveDown =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveDown}, PSI.RangeReserve),
    )
    devices = Dict(
        :Generators =>
            PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalStandardUnitCommitment),
        :Ren => PSI.DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch),
        :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
        :HydroROR =>
            PSI.DeviceModel(PSY.HydroEnergyReservoir, PSI.HydroDispatchRunOfRiver),
        :RenFx => PSI.DeviceModel(PSY.RenewableFix, PSI.FixedOutput),
        :Batt => PSI.DeviceModel(PSY.GenericBattery, PSI.BookKeepingwReservation),
    )

    template_uc = PSI.template_unit_commitment(
        network = PSI.PM.NFAPowerModel,
        devices = devices,
        services = service,
    )
    return template_uc
end

function create_ed_template()
    service = Dict(
        :ReserveUp =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveUp}, PSI.RangeReserve),
        :ReserveDown =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveDown}, PSI.RangeReserve),
    )
    devices = Dict(
        :Generators => PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch), # ramp limited creats infeasibilities
        :Ren => PSI.DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch),
        :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
        :HydroROR =>
            PSI.DeviceModel(PSY.HydroEnergyReservoir, PSI.HydroDispatchRunOfRiver),
        :RenFx => PSI.DeviceModel(PSY.RenewableFix, PSI.FixedOutput),
        :Batt => PSI.DeviceModel(PSY.GenericBattery, PSI.BookKeeping),
    )
    template_ed = PSI.template_economic_dispatch(
        network = PSI.PM.NFAPowerModel,
        devices = devices,
        services = service,
    )

    return template_ed
end

function create_continuous_template()
    service = Dict(
        :ReserveUp =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveUp}, PSI.RangeReserve),
        :ReserveDown =>
            PSI.ServiceModel(PSY.VariableReserve{PSY.ReserveDown}, PSI.RangeReserve),
    )
    devices = Dict(
        :Generators => PSI.DeviceModel(PSY.ThermalStandard, PSI.ThermalRampLimited),
        :Ren => PSI.DeviceModel(PSY.RenewableDispatch, PSI.RenewableFullDispatch),
        :Loads => PSI.DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad),
        :HydroROR =>
            PSI.DeviceModel(PSY.HydroEnergyReservoir, PSI.HydroDispatchRunOfRiver),
        :RenFx => PSI.DeviceModel(PSY.RenewableFix, PSI.FixedOutput),
        :Batt => PSI.DeviceModel(PSY.GenericBattery, PSI.BookKeeping),
    )

    template_uc = PSI.template_economic_dispatch(
        network = PSI.PM.NFAPowerModel,
        devices = devices,
        services = service,
    )
    return template_uc
end

function create_stage(
    template::PSI.OperationsProblemTemplate,
    sys::PSY.System,
    solver::JuMP.MOI.OptimizerWithAttributes;
    kwargs...,
)
    constraint_duals = get(kwargs, :constraint_duals, nothing)
    if isnothing(constraint_duals)
        stage = PSI.Stage(
            PSI.GenericOpProblem,
            template,
            sys,
            solver;
            balance_slack_variables = true,
            services_slack_variables = true,
        )
        return stage
    else
        stage = PSI.Stage(
            PSI.GenericOpProblem,
            template,
            sys,
            solver;
            balance_slack_variables = true,
            services_slack_variables = true,
            export_pwl_vars = true,
            constraint_duals = constraint_duals,
        )
        return stage
    end
end

function create_sequence(stage_1::AbstractString, stage_2::AbstractString)
    sequence = PSI.SimulationSequence(
        step_resolution = Dates.Hour(24),
        order = Dict(1 => stage_1, 2 => stage_2),
        feedforward_chronologies = Dict(
            (stage_1 => stage_2) => PSI.Synchronize(periods = 24),
        ),
        horizons = Dict(stage_1 => 24, stage_2 => 1),
        intervals = Dict(
            stage_1 => (Dates.Hour(24), PSI.Consecutive()),
            stage_2 => (Dates.Hour(1), PSI.Consecutive()),
        ),
        feedforward = Dict(
            (stage_2, :devices, :Generators) => PSI.SemiContinuousFF(
                binary_source_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        cache = Dict(
            (stage_1,) => PSI.TimeStatusChange(PSY.ThermalStandard, PSI.ON),
            (stage_1,) => PSI.StoredEnergy(PSY.GenericBattery, PSI.ENERGY),
        ),
        ini_cond_chronology = PSI.InterStageChronology(),
    )
    return sequence
end

function change_load_profile!(rawsys::PSY.PowerSystemTableData, multiplier::Float64)
    rt_dir = joinpath(
        dirname(rawsys.directory),
        "timeseries_data_files/Load/REAL_TIME_regional_Load.csv",
    )
    da_dir = joinpath(
        dirname(rawsys.directory),
        "timeseries_data_files/Load/DAY_AHEAD_regional_Load.csv",
    )
    rt_load = CSV.read(rt_dir)
    da_load = CSV.read(da_dir)
    rt_load[!, [Symbol(1), Symbol(2), Symbol(3)]] =
        rt_load[!, [Symbol(1), Symbol(2), Symbol(3)]] .* multiplier
    da_load[!, [Symbol(1), Symbol(2), Symbol(3)]] =
        da_load[!, [Symbol(1), Symbol(2), Symbol(3)]] .* multiplier
    CSV.write(rt_dir, rt_load)
    CSV.write(da_dir, da_load)
    return
end

function build_training_list(sys::PSY.System)
    category_dict = Dict()
    capacity_dict = Dict()
    reg_capacity_dict = Dict()
    technology = [
        PSY.PrimeMovers.ST,
        PSY.PrimeMovers.CT,
        PSY.PrimeMovers.CC,
        PSY.PrimeMovers.WT,
        PSY.PrimeMovers.HY,
        PSY.PrimeMovers.PVe,
        PSY.PrimeMovers.BA,
    ]
    devices = collect(PSY.get_components(PSY.Generator, sys))
    storage_devices = collect(PSY.get_components(PSY.Storage, sys))
    devices = vcat(devices, storage_devices)
    buses = collect(PSY.get_components(PSY.Bus, sys))
    for tech in technology
        _devices = filter(x -> x.primemover == tech, devices)
        if isempty(_devices)
            category_dict[Symbol(tech)] = []
            capacity_dict[Symbol(tech)] = 0.0
        else
            category_dict[Symbol(tech)] = PSY.get_name.(_devices)
            capacity_dict[Symbol(tech)] = sum(get_capacity.(_devices))
        end
        for bus in buses
            if !isempty(_devices)
                reg_devices = filter(x -> PSY.get_name(x.bus) == PSY.get_name(bus), _devices)
                if isempty(reg_devices)
                    reg_capacity_dict[(Symbol(tech), PSY.get_name(bus))] = 0.
                else
                    reg_capacity_dict[(Symbol(tech), PSY.get_name(bus))] = sum(get_capacity.(reg_devices))
                end
            else
                reg_capacity_dict[(Symbol(tech), PSY.get_name(bus))] = 0.
            end
        end
    end
    return category_dict, capacity_dict, reg_capacity_dict
end

get_capacity(v::PSY.ThermalGen) = v.activepowerlimits.max
get_capacity(v::PSY.HydroGen) = v.activepowerlimits.max
get_capacity(v::PSY.RenewableGen) = v.rating
get_capacity(v::PSY.Storage) = v.capacity.max

function create_simulation(
    file_path::AbstractString,
    days::Int64,
    pruned_unit::Dict{Type{<:PSY.Component}, Array{AbstractString}},
    device_list::Dict,
    time_series_directory::AbstractString,
    rts_dir = nothing,
    load_multiplier = 1.0,
    ;
    kwargs...,
)
    rts_dir = get_rts_data(; rts_dir = rts_dir)
    rawsys = read_rts_data(rts_dir)
    simulation_name = get(kwargs, :simulation_name, "surrogate_epec")
    optimizer = get(kwargs, :optimizer, Cbc_optimizer)
    sys_UC = PSY.System(
        rawsys;
        forecast_resolution = Dates.Hour(1),
        time_series_directory = time_series_directory,
    )

    add_copies!(sys_UC, device_list)

    prune_system_devices!(sys_UC, pruned_unit)
    template_uc = create_uc_template()
    template_ed = create_ed_template()

    uc_stage = create_stage(template_uc, sys_UC, optimizer)
    ed_stage = create_stage(
        template_ed,
        sys_UC,
        optimizer;
        constraint_duals = [
            Symbol("nodal_balance_active__PowerSystems.Bus"),
            Symbol("requirement__PowerSystems.VariableReserve_PowerSystems.ReserveDown"),
            Symbol("requirement__PowerSystems.VariableReserve_PowerSystems.ReserveUp"),
        ],
    )

    stages_definition = Dict("UC" => uc_stage, "ED" => ed_stage)

    sequence = create_sequence("UC", "ED")

    sim = PSI.Simulation(
        name = simulation_name,
        steps = days,
        stages = stages_definition,
        stages_sequence = sequence,
        simulation_folder = file_path,
        initial_time = Dates.DateTime("2020-01-01T00:00:00"),
    )
    PSI.build!(sim)

    sim_results = PSI.execute!(sim;)

    # res_uc = load_simulation_results(sim_results, "UC");
    res_ed = PSI.load_simulation_results(sim_results, "ED")
    PSI.write_to_CSV(res_ed)
    category_dict, capacity_dict, reg_capacity_dict = build_training_list(sys_UC)
    bus_mapping = Dict(
        PSY.get_name(d) => PSY.get_name(PSY.get_bus(d))
        for d in PSY.get_components(PSY.StaticInjection, sys_UC) if (typeof(d) <: PSY.Generator) || (typeof(d) <: PSY.Storage)
    )
    capacity_mapping = Dict(
        PSY.get_name(d) => get_capacity(d)
        for d in PSY.get_components(PSY.StaticInjection, sys_UC) if (typeof(d) <: PSY.Generator) || (typeof(d) <: PSY.Storage)
    )
#     reg_capacity_dict = Dict(
#         PSY.get_name(d) => get_capacity(d) for d in PSY.get_components(PSY.Generator, sys_UC)
#     )
    res_mapping = Dict(
        PSY.get_name(d) => PSY.get_name.(PSY.get_services(d))
        for d in PSY.get_components(PSY.StaticInjection, sys_UC) if (typeof(d) <: PSY.Generator) || (typeof(d) <: PSY.Storage)
    )
    results_dir = sim_results.results_folder
    training_data, reg_training_data, cap_price, cap_accepted_bid = calculate_revenue(
        res_ed,
        category_dict,
        capacity_dict,
        capacity_mapping,
        bus_mapping,
        res_mapping,
        reg_capacity_dict,
        optimizer,
    )
    CSV.write(results_dir * "/training_data.csv", training_data)
    CSV.write(results_dir * "/reg_training_data.csv", reg_training_data)
    CSV.write(results_dir * "/Capacity_market_price.csv", Dict("price" => cap_price))
    CSV.write(results_dir * "/Capacity_market_accepted_bids.csv", cap_accepted_bid)
    return
end

function create_continuous_simulation(
    file_path::AbstractString,
    days::Int64,
    pruned_unit::Dict{Type{<:PSY.Component}, Array{AbstractString}},
    device_list::Dict,
    time_series_directory::AbstractString,
    rts_dir = nothing,
    load_multiplier = 1.0;
    kwargs...,
)

    rts_dir = get_rts_data(; rts_dir = rts_dir)
    rawsys = read_rts_data(rts_dir)
    simulation_name = get(kwargs, :simulation_name, "surrogate_epec")
    optimizer = get(kwargs, :optimizer, Cbc_optimizer)
    sys = PSY.System(
        rawsys;
        forecast_resolution = Dates.Hour(1),
        time_series_directory = time_series_directory,
    )
    add_copies!(sys, device_list)
    prune_system_devices!(sys, pruned_unit)
    template = create_continuous_template()
    stage = create_stage(
        template,
        sys,
        optimizer;
        constraint_duals = [
            Symbol("nodal_balance_active__PowerSystems.Bus"),
            Symbol("requirement__PowerSystems.VariableReserve_PowerSystems.ReserveDown"),
            Symbol("requirement__PowerSystems.VariableReserve_PowerSystems.ReserveUp"),
        ],
    )
    stages_definition = Dict("ED" => stage)
    sequence = PSI.SimulationSequence(
        step_resolution = Dates.Hour(24),
        order = Dict(1 => "ED"),
        horizons = Dict("ED" => 24),
        intervals = Dict("ED" => (Dates.Hour(24), PSI.Consecutive())),
        ini_cond_chronology = PSI.IntraStageChronology(),
    )
    sim = PSI.Simulation(
        name = simulation_name,
        steps = days,
        stages = stages_definition,
        stages_sequence = sequence,
        simulation_folder = file_path,
        initial_time = Dates.DateTime("2020-01-01T00:00:00"),
    )
    PSI.build!(sim)
    sim_results = PSI.execute!(sim;)
    res_ed = PSI.load_simulation_results(sim_results, "ED")
    PSI.write_to_CSV(res_ed)
    category_dict, capacity_dict, reg_capacity_dict = build_training_list(sys)
    bus_mapping = Dict(
        PSY.get_name(d) => PSY.get_name(PSY.get_bus(d))
        for d in PSY.get_components(PSY.Generator, sys)
    )
    capacity_mapping = Dict(
        PSY.get_name(d) => get_capacity(d) for d in PSY.get_components(PSY.Generator, sys)
    )
#     reg_capacity_dict = Dict(
#         PSY.get_name(d) => get_capacity(d) for d in PSY.get_components(PSY.Generator, sys)
#     )
    res_mapping = Dict(
        PSY.get_name(d) => PSY.get_name.(PSY.get_services(d))
        for d in PSY.get_components(PSY.Generator, sys)
    )
    results_dir = sim_results.results_folder
    training_data, reg_training_data, cap_price, cap_accepted_bid = calculate_revenue(
        res_ed,
        category_dict,
        capacity_dict,
        capacity_mapping,
        bus_mapping,
        res_mapping,
        reg_capacity_dict,
        optimizer,
    )
    CSV.write(results_dir * "/training_data.csv", training_data)
    CSV.write(results_dir * "/reg_training_data.csv", reg_training_data)
    CSV.write(results_dir * "/Capacity_market_price.csv", Dict("price" => cap_price))
    CSV.write(results_dir * "/Capacity_market_accepted_bids.csv", cap_accepted_bid)
    return
end


text:
struct Node
    left::Int
    right::Int
end

struct TreePartition{N,T} <: AbstractBoxPartition{Box{N,T}}
    domain::Box{N,T}
    nodes::Vector{Node}
    regular_partitions::Vector{BoxPartition{N,T}}
end

Base.copy(partition::TreePartition) = TreePartition(partition.domain, copy(partition.nodes), copy(partition.regular_partitions))

depth(partition::TreePartition) = length(partition.regular_partitions) - 1
keytype(::Type{<:TreePartition}) = Tuple{Int,Int}

function Base.show(io::IO, partition::TreePartition)
    print(io, "TreePartition of depth $(depth(partition))")
end

function keys_all(partition::TreePartition)
    if depth(partition) != 0
        error("not implemented")
    end

    return [(0, 1)]
end

TreePartition(domain::Box) = TreePartition(domain, [Node(0, 0)], [BoxPartition(domain)])

dimension(partition::TreePartition{N,T}) where {N,T} = N

function key_to_box(partition::TreePartition, key::Tuple{Int,Int})
    return key_to_box(partition.regular_partitions[key[1] + 1], key[2])
end

#=
function unsafe_point_to_ints(partition::BoxPartition, point) # should maybe be put into partition_regular.jl
    x = (point .- partition.left) .* partition.scale    # counts how many boxes x is away from left (componentwise)
    return map(xi -> Base.unsafe_trunc(Int, xi), x)
end
=#

function tree_search(tree::TreePartition{N,T}, point) where {N,T}
    regular_partitions = tree.regular_partitions

    # start at root
    node_idx = 1
    tree_depth = depth(tree)
    current_depth = 0
    ints = zeros(SVector{N,Int})

    while current_depth+1 <= tree_depth
        ints_next = unsafe_point_to_ints(regular_partitions[current_depth+2], point)

        point_is_left = iseven(ints_next[(current_depth % N) + 1])

        node = tree.nodes[node_idx]
        node_idx_next = ifelse(point_is_left, node.left, node.right)

        if node_idx_next == 0
            break
        end

        node_idx = node_idx_next

        ints = ints_next
        current_depth += 1
    end

    key = (current_depth, sum(ints .* regular_partitions[current_depth+1].dimsprod) + 1)

    return key, node_idx
end

function point_to_key(partition::TreePartition, point)
    if point_to_key(partition.regular_partitions[1], point) === nothing
        return nothing
    end

    return tree_search(partition, point)[1]
end

function subdivide!(tree::TreePartition{N, T}, key::Tuple{Int,Int}) where {N, T}
    search_key, node_idx = tree_search(tree, key_to_box(tree, key).center)

    @assert key == search_key

    node = tree.nodes[node_idx]

    @assert node.left == 0 && node.right == 0

    new_node = Node(length(tree.nodes) + 1, length(tree.nodes) + 2)
    push!(tree.nodes, Node(0, 0))
    push!(tree.nodes, Node(0, 0))

    tree.nodes[node_idx] = new_node

    if key[1] == depth(tree)
        push!(tree.regular_partitions, BoxPartition(tree.domain, depth=depth(tree)+1 #=ntuple(i -> depth(tree)+1, N)=#))
    end

    return tree
end


text:
normalize(arr) = arr ./ sum(arr)

function generateAdjacencyMatrix(R_t::Array{Int,2}, C::Function)
	assert(size(R_t,2) == 2)
	A = zeros(Int,max(R_t[:,1]),max(R_t[:,2]))
	for tag = 1:size(R_t,1)
		u = R_t[tag,1]
		d = R_t[tag,2]
		A[A[:,d] .>= 1,d] += 1
		A[u,d] = 1
	end
	return sparse(C(A))
end

function runSpear(A::SparseMatrixCSC{Float64,Int64}, k = 10)
	M = size(A,1)    # number of unique users
	N = size(A,2)    # number of unique documents
	E = ones(M)'     # experience scors of users
	Q = ones(N)'     # qality scors of documents

	for i = 1:k # calculation of expertise and quality stores
		E = normalize(Q * A')
		Q = normalize(E * A)
	end

	return E, Q

end


text:
module GeoCommon

  type SimpleLocation
    name::String
    latitude::Float64
    longitude::Float64
  end

end

max_stars_repo_name: astrieanna/geojl

max_issues_repo_name: astrieanna/geojl

max_forks_repo_name: astrieanna/geojl


text:
# Autogenerated wrapper script for cuQuantum_jll for aarch64-linux-gnu
export libcustatevec, libcutensornet

JLLWrappers.@generate_wrapper_header("cuQuantum")
JLLWrappers.@declare_library_product(libcustatevec, "libcustatevec.so.0")
JLLWrappers.@declare_library_product(libcutensornet, "libcutensornet.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcustatevec,
        "lib/libcustatevec.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libcutensornet,
        "lib/libcutensornet.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Node

TODO: make this immutable

Record path in a maze
"""
mutable struct Node{P}
point::P #Tuple{Int64,Int64}  # MazeLocation
parent::Union{Node{P},Nothing}
cost::Float64
heuristic::Float64

# Node() = ( # incomplete initialization
#     me = new();
#     #me.cost = 0.0;
#     #me.heuristic = 0.0;
#     me
# )
function Node{P}() where {P}
    me = new()
    return me
end
# Node(point) = (
#     me = new();
#     me.point = point;
#     me.parent = nothing;
#     me
# )
function Node{P}(point::P) where {P}
      me = new();
      me.point = point;
      me.parent = nothing;
      me
end
# Node(point,parent) = (
#     me = new();
#     me.point = point;
#     me.parent = parent;
#     me
# )
function Node{P}(point::P, parent::Node{P}) where {P}
      me = new();
      me.point = point;
      me.parent = parent;
      me
end
# Node(point,parent,cost,heuristic) = (
#     me = new();
#     me.point = point;
#     me.parent = parent;
#     me.cost = cost;
#     me.heuristic = heuristic;
#     me
# )
end
#Base.:(<)(x::Node, y::Node) = x.id < y.id
Base.:(==)(x::Node,y::Node) = x.point == y.point

"""
DijkstraNode

Node for use in Dijkstra algorithm
"""
struct DijkstraNode
    index::Int64
    distance::Float64
    # function DijkstraNode(index, distance)
    #     new(index, distance)
    # end
end
# make sure not mis-use these 2 functions
Base.:(<)(x::DijkstraNode, y::DijkstraNode) = error("Unsupported operation")# x.distance < y.distance
Base.:(==)(x::DijkstraNode, y::DijkstraNode) = error("Unsupported operation") #x.distance == y.distance #error("Unsupported operation") #x.index == y.index
# to make Dict and PriorityQueue behavior the way we wanted
Base.isequal(x::DijkstraNode, y::DijkstraNode) = isequal(x.index, y.index) #x.index == y.index #error("Unsupported operation")
Base.hash(x::DijkstraNode, h::UInt64=UInt64(13)) = hash(x.index,h)

"""
    node_to_path(node::Node{P})::Vector{P} where {P}

Collect `point` from `Node` and its parents to show the path
from the `node` all the way back until parent is `nothing`.
"""
function node_to_path(node::Node{P})::Vector{P} where {P}
    path = Vector{P}()
    push!(path, node.point)
    while !isnothing(node.parent)
        node = node.parent
        #push!(path, node.point)
        insert!(path, 1, node.point)
    end

    return path
end


text:
"""
	pvaluelocal(d::ContinuousUnivariateDistribution, x::Number ; tail::Symbol=:both)::Float64
	pvaluelocal{T<:Number}(d::Vector{T}, x::Number ; tail::Symbol=:both, as::Bool=false)::Float64

Obtain the p-value associated with the inputs.
d::ContinuousUnivariateDistribution -> d is the distribution under the null and x is the test statistic.
d::Vector{T} -> d is a bootstrapped vector of test statistics and x is the value of the test statistic under the null.
"""
function pvaluelocal(d::ContinuousUnivariateDistribution, x::Number ; tail::Symbol=:both)::Float64
	if tail == :both
		lpv = cdf(d, x)
		rpv = ccdf(d, x)
		lpv < rpv && return min(2*lpv, 1.0)
		return min(2*rpv, 1.0)
	elseif tail == :left
	    return cdf(d, x)
	elseif tail == :right
	    return ccdf(d, x)
	else
		error("Keyword argument tail=$(tail) is invalid")
	end
end
function pvaluelocal(xVec::Vector{T}, x::Number ; tail::Symbol=:both, as::Bool=false)::Float64 where {T<:Number}
	!as && (xVec = sort(xVec))
	i = searchsortedlast(xVec, x)
	if tail == :both
		lpv = Float64(i / length(xVec))
		rpv = 1.0 - lpv
		lpv < rpv && return min(2*lpv, 1.0)
		return min(2*rpv, 1.0)
	elseif tail == :left
		return Float64(i / length(xVec))
	elseif tail == :right
		return Float64(1.0 - (i / length(xVec)))
	else
		error("Keyword argument tail=$(tail) is invalid")
	end
end


text:
function test_progress_messenger(messenger; LES=false, model_kwargs...)
    model = create_model(; model_kwargs...)
    test_progress_messenger(model, messenger)
end

function test_progress_messenger(model, messenger; LES=false)
    simulation = Simulation(model; Δt=1e-2,
                            stop_iteration=10,
                            iteration_interval=1,
                            progress=messenger,
                            )

    run!(simulation)

    return nothing
end


text:
function _current_exceptions(codestate, parms)
    if isempty(parms) || parms[1] == current_task()
        return codestate.interpstate.exceptions
    else
        return Compat.current_exceptions(parms...)
    end
end
function _catch_backtrace(codestate, parms)
    if isempty(codestate.interpstate.exceptions)
        return []
    else
        return last(codestate.interpstate.exceptions).backtrace
    end
end
function _rethrow(codestate, parms)
    if isempty(parms)
        if isempty(codestate.interpstate.exceptions)
            # parms = (ErrorException("rethrow() not allowed outside a catch block"),)
            parms = Any[ErrorException("rethrow() not allowed outside a catch block")]
        else
            # parms = (pop!(codestate.interpstate.exceptions).exception,)
            parms = Any[pop!(codestate.interpstate.exceptions).exception]
        end
    elseif isempty(codestate.interpstate.exceptions)
        # parms = (ErrorException("rethrow(exc) not allowed outside a catch block"),)
        parms = Any[ErrorException("rethrow(exc) not allowed outside a catch block")]
    end
    return Base.throw(parms...)
end
# https://github.com/JuliaLang/julia/issues/43556
function _iolock_begin(codestate, parms)
    # @show :_iolock_begin codestate.interpstate.iolock codestate.src
    # if !codestate.interpstate.iolock
    #     Base.iolock_begin(parms...)
    #     codestate.interpstate.iolock = true
    # end
    return nothing
end
function _iolock_end(codestate, parms)
    # @show :_iolock_end codestate.interpstate.iolock codestate.src
    # if codestate.interpstate.iolock
    #     codestate.interpstate.iolock = false
    #     Base.iolock_end(parms...)
    # end
    return nothing
end
function _sigatomic_begin(codestate, parms)
    # @show :_sigatomic_begin codestate.interpstate.sigatomic codestate.src
    # if !codestate.interpstate.sigatomic
        Base.sigatomic_begin(parms...)
        codestate.interpstate.sigatomic = true
    # end
    return nothing
end
function _sigatomic_end(codestate, parms)
    # @show :_sigatomic_end codestate.interpstate.sigatomic codestate.src
    if codestate.interpstate.sigatomic
        codestate.interpstate.sigatomic = false
        Base.sigatomic_end(parms...)
    end
    return nothing
end
function _llvmcall(codestate, parms)
    # @show :_llvmcall parms
    funname = gensym("llvmcall")
    argnames = ((Symbol(:(_), i) for i = 1:length(parms) - 3)...,)
    ir = parms[1]
    rt = parms[2]
    at = parms[3]
    return eval_ast(codestate.interpstate, moduleof(codestate.mod_or_meth),
        quote
            function $funname($(argnames...))
                Base.llvmcall($ir, $rt, $at, $(argnames...))
            end
            $funname($(parms[4:end]...))
        end)
end
function _eval(codestate, parms)
    wc1 = Base.get_world_counter()
    @show :_eval Int(wc1)
    if length(parms) == 1
        ans = interpret_ast(moduleof(codestate.mod_or_meth), parms[1],
            codestate.interpstate.options)
    else
        ans = interpret_ast(parms[1], parms[2],
            codestate.interpstate.options)
    end
    wc2 = Base.get_world_counter()
    @show :_eval Int(wc2)
    return ans
end

struct InterpOptions
    debug::Bool
    stats::Bool
    modules::Vector{Module}
    limiteds::Vector{Module}
    budget::UInt
end

function options(debug, stats, modules, limiteds, budget)
    InterpOptions(debug, stats, modules, limiteds, budget)
end

const Callsite = Union{Ptr{Nothing}, Base.InterpreterIP}

const ModuleOrMethod = Union{Module, Method}

moduleof(meth::Method) = meth.module
moduleof(mod::Module) = mod

const WorldType = Tuple{UInt, Type}

typeof_lower(type::Type) = Type{type}
typeof_lower(val) = typeof(val)

function world_type(callable, parms)
    wc = Base.get_world_counter()
    t = Base.to_tuple_type(typeof_lower.(parms))
    tt::Type = Base.signature_type(callable, t)
    (wc, tt)
end

function _which(wt)
    #= try
        return which(wt[2])
    catch exception
        @show :_which exception
        if !(exception isa ErrorException)
            rethrow()
        else
            return nothing
        end
    end =#
    @static if VERSION >= v"1.8.0-DEV"
        # Adapted Base._which
        min_valid = Base.RefValue{UInt}(typemin(UInt))
        max_valid = Base.RefValue{UInt}(typemax(UInt))
        match = ccall(:jl_gf_invoke_lookup_worlds, Any,
            (Any, UInt, Ptr{Csize_t}, Ptr{Csize_t}),
            wt[2], wt[1], min_valid, max_valid)
        return match
        #= if match isa Core.MethodMatch
            return match.method
        else
            return nothing
        end =#
    else
        return ccall(:jl_gf_invoke_lookup, Any, (Any, UInt), wt[2], wt[1])
    end
end

mutable struct Statistics
    calls::UInt
    elapsed::UInt
end

mutable struct InterpState
    options::InterpOptions
    wc::UInt
    intercepts::Dict{Base.Callable, Base.Callable}
    meths::Dict{WorldType, Tuple{Method, Core.SimpleVector, Core.CodeInfo}}
    stats::Dict{ModuleOrMethod, Statistics}
    exceptions::Vector{NamedTuple{(:exception, :backtrace), Tuple{Any, Vector{<:Callsite}}}}
    iolock::Bool
    sigatomic::Bool
end

function interp_state(mod, options)
    wc = Base.get_world_counter()
    interpstate = InterpState(options, wc, Dict(), Dict(), Dict(), [], false, false)
    @static if !isdefined(Base, :current_exceptions)
        interpstate.intercepts[Base.catch_stack] = _current_exceptions
    else
        interpstate.intercepts[Base.current_exceptions] = _current_exceptions
    end
    interpstate.intercepts[Base.catch_backtrace] = _catch_backtrace
    interpstate.intercepts[Base.rethrow] = _rethrow
    interpstate.intercepts[Base.iolock_begin] = _iolock_begin
    interpstate.intercepts[Base.iolock_end] = _iolock_end
    interpstate.intercepts[Base.sigatomic_begin] = _sigatomic_begin
    interpstate.intercepts[Base.sigatomic_end] = _sigatomic_end
    interpstate.intercepts[Base.llvmcall] = _llvmcall
    interpstate.intercepts[Core.Intrinsics.llvmcall] = _llvmcall
    # interpstate.intercepts[Core.eval] = _eval
    # interpstate.intercepts[Base.eval] = _eval
    # interpstate.intercepts[mod.eval] = _eval
    # union!(interpstate.passthroughs, intrinsics())
    # union!(interpstate.passthroughs, builtins())
    #= push!(interpstate.passthroughs, Core.eval)
    push!(interpstate.passthroughs, Base.eval)
    push!(interpstate.passthroughs, getfield(mod, :eval))
    push!(interpstate.passthroughs, Base.lock)
    push!(interpstate.passthroughs, Base.unlock) =#
    interpstate
end

function update_interp_statistics(interpstate::InterpState, mod_or_meth, time)
    if !haskey(interpstate.stats, mod_or_meth)
        interpstate.stats[mod_or_meth] = Statistics(UInt(1), time)
    else
        interpstate.stats[mod_or_meth].calls += UInt(1)
        interpstate.stats[mod_or_meth].elapsed += time
    end
end

function check_interp_statistics(interpstate::InterpState, mod_or_meth)
    interpstate.options.budget == 0 ||
    !haskey(interpstate.stats, mod_or_meth) ||
    interpstate.stats[mod_or_meth].elapsed < interpstate.options.budget
end

mutable struct CodeState
    interpstate::InterpState
    wt::WorldType
    mod_or_meth::ModuleOrMethod
    sparam_vals::Core.SimpleVector
    src::Core.CodeInfo
    pc::Int
    cc::Int
    ssavalues::Vector{Any}
    childstates::Vector{Vector{Tuple{Base.Callable, Union{Nothing, CodeState}}}}
    childstats::Vector{Statistics}
    slots::Vector{Any}
    handlers::Vector{Int}
end

const ChildState = Tuple{Base.Callable, Union{Nothing, CodeState}}

function code_state_from_thunk(interpstate, mod, src)
    # @show :code_state_from_thunk
    wc = Base.get_world_counter()
    CodeState(interpstate, (wc, Nothing), mod, Core.svec(), src, 1, 1,
        Vector{Any}(undef, length(src.code)),
        Vector{Vector{ChildState}}(undef, length(src.code)),
        Vector{Statistics}(undef, length(src.code)),
        Vector{Any}(undef, length(src.slotnames)),
        Int[])
end

generate_lower(codestate::CodeState, expr::Expr) = Meta.lower(moduleof(codestate.mod_or_meth), expr)
generate_lower(codestate::CodeState, val) = val

function code_state_from_call(codestate::CodeState, wt)
    @nospecialize wt
    # @show :code_state_from_call
    if !haskey(codestate.interpstate.meths, wt)
        @static if VERSION >= v"1.8.0-DEV"
            match = _which(wt)
            if match === nothing
                @show :code_state_from_call wt match
                return nothing
            end
            meth = match.method
            # mi = Core.Compiler.specialize_method(match)
            # sparam_vals = mi.sparam_vals
            (ti, sparam_vals) = Core.Compiler.normalize_typevars(meth, wt[2], match.sparams)
        else
            meth = _which(wt)
            if meth === nothing
                @show :code_state_from_call wt meth
                return nothing
            end
            (ti, sparam_vals) = ccall(:jl_type_intersection_with_env, Any, (Any, Any), wt[2], meth.sig)
        end
        if !isdefined(meth, :generator)
            src = Base.uncompressed_ast(meth)
            # src = Base.uncompressed_ir(meth)
        else
            local expr_or_src
            try
                generator = meth.generator
                if codestate.mod_or_meth isa Module
                    expr_or_src = Base.invokelatest(generator, sparam_vals..., generator.argnames...)
                else
                    # expr_or_src = generator(sparam_vals..., generator.argnames...)
                    expr_or_src = Base.invokelatest(generator, sparam_vals..., generator.argnames...)
                end
            catch exception
                @show :code_state_from_call exception
                return nothing
            end
            if expr_or_src isa Core.CodeInfo
                src = expr_or_src
            elseif expr_or_src isa Expr
                expr = expr_or_src
                src = generate_lower(codestate, expr)
            else
                @show expr_or_src
                @assert false
            end
        end
        codestate.interpstate.meths[wt] = (meth, sparam_vals, src)
    else
        meth, sparam_vals, src = codestate.interpstate.meths[wt]
    end
    CodeState(codestate.interpstate, wt, meth, sparam_vals, src, 1, 1,
        Vector{Any}(undef, length(src.code)),
        Vector{Vector{ChildState}}(undef, length(src.code)),
        Vector{Statistics}(undef, length(src.code)),
        Vector{Any}(undef, length(src.slotnames)),
        Int[])
end

function update_code_state(codestate::CodeState, callable, parms)
    codestate.pc = 1
    codestate.ssavalues = Vector{Any}(undef, length(codestate.src.code))
    # codestate.childstats = Vector{Statistics}(undef, length(codestate.src.code))
    codestate.slots = Vector{Any}(undef, length(codestate.src.slotnames))
    empty!(codestate.handlers)
    meth = codestate.mod_or_meth::Method
    codestate.slots[1] = callable
    if !meth.isva
        @views codestate.slots[2:meth.nargs] .= parms[1:meth.nargs - 1]
    else
        @views codestate.slots[2:meth.nargs - 1] .= parms[1:meth.nargs - 2]
        # codestate.slots[meth.nargs] = parms[meth.nargs - 1:end]
        codestate.slots[meth.nargs] = (parms[meth.nargs - 1:end]...,)
    end
    codestate
end

function update_code_statistics(codestate::CodeState, time)
    if !isassigned(codestate.childstats, codestate.pc)
        codestate.childstats[codestate.pc] = Statistics(UInt(1), time)
    else
        codestate.childstats[codestate.pc].calls += UInt(1)
        codestate.childstats[codestate.pc].elapsed += time
    end
end

function check_code_statistics(codestate::CodeState)
    codestate.interpstate.options.budget == 0 ||
    !isassigned(codestate.childstats, codestate.pc) ||
    codestate.childstats[codestate.pc].elapsed < codestate.interpstate.options.budget
end


text:
function exch(::Scenario, from::C, to::C) :: Float64 where {C<:Currencies.Currency}
    # same currency
    return 1.0
end

function exch(p::Scenario, from::C1, to::C2) :: Float64 where {C1<:Currencies.Currency, C2<:Currencies.Currency}

    if !haskey(p, from) # ex.: BRL -> USD
        @assert haskey(p, to)
        @assert Currencies.cashcurrency(p[to]) == from
        return 1.0 / Currencies.cashvalue(p[to])

    elseif Currencies.cashcurrency(p[from]) == to # ex.: USD -> BRL
        return Currencies.cashvalue(p[from])

    else
        @assert haskey(p, to)
        @assert Currencies.cashcurrency(p[from]) == Currencies.cashcurrency(p[to]) # same functional currency
        return Currencies.cashvalue(p[from]) / Currencies.cashvalue(p[to])
    end
end

@inline function exch(p::Scenario, from_cash::Currencies.Cash{C1, T}, to_currency::C2) :: T where {C1<:Currencies.Currency, C2<:Currencies.Currency, T<:Real}
    return exch(p, C1(), C2()) * Currencies.cashvalue(from_cash)
end

function exchcash(p::Scenario, from::C1, to::C2) :: Currencies.Cash{C2, Float64} where {C1<:Currencies.Currency, C2<:Currencies.Currency}
    return exch(p, from, to) * to
end


text:
import LinearAlgebra: norm_sqr, ⋅

import ..Constraints: FixedDistanceConstraint
import ..Types: move!
import ..Utils: acceleration

import ..Verlet

@inline function _move_half_step!(
    r::Vector3s, v::Vector3s,
    dt::Float64, tolerance::Float64, m::Function, dist::Function,
    max_iter::Int, last_r::Vector3s, constraints::Vector{FixedDistanceConstraint}
)
    N = length(r)
    K = length(constraints)
    moved = ones(Bool, N)
    moving = zeros(Bool, N)
    iter = 0
    correcting = true
    while correcting && iter < max_iter
        iter += 1
        correcting = false
        @inbounds for ki in 1:K
            constraint = constraints[ki]
            i, j, d = constraint.i, constraint.j, constraint.distance
            !(moved[i] || moved[j]) && continue
            d² = d ^ 2
            ξ = 2 * tolerance * d²
            rij = dist(last_r[i], last_r[j])
            s = dist(r[i], r[j])
            rmi, rmj = 1 / m(i), 1 / m(j)
            c = norm_sqr(s) - d²
            if abs(c) ≥ ξ
                correcting = true
                moving[i] = moving[j] = true
                g = c / 2dt / (s ⋅ rij) / (rmi + rmj)
                dvi = -rij * rmi * g
                dvj = rij * rmj * g
                v[i] += dvi
                v[j] += dvj
                r[i] += dvi * dt
                r[j] += dvj * dt
            end
        end
        moved .= moving
        moving .= false
    end
end

@inline function _move_full_step!(
    r::Vector3s, v::Vector3s,
    tolerance::Float64, m::Function, dist::Function,
    max_iter::Int, constraints::Vector{FixedDistanceConstraint}
)
    N = length(r)
    K = length(constraints)
    moved = ones(Bool, N)
    moving = zeros(Bool, N)
    iter = 0
    correcting = true
    while correcting && iter < max_iter
        iter += 1
        correcting = false
        @inbounds for ki in 1:K
            constraint = constraints[ki]
            i, j, d = constraint.i, constraint.j, constraint.distance
            !(moved[i] || moved[j]) && continue
            d² = d ^ 2
            ξ = tolerance * d²
            rij = dist(r[i], r[j])
            vij = v[i] - v[j]
            rmi, rmj = 1 / m(i), 1 / m(j)
            c = rij ⋅ vij
            if abs(c) ≥ ξ
                correcting = true
                moving[i] = moving[j] = true
                k = c / d² / (rmi + rmj)
                v[i] += -rij * rmi * k
                v[j] += rij * rmj * k
            end
        end
        moved .= moving
        moving .= false
    end
end

function move!(integrator::RattleIntegrator)
    r = integrator.positions
    v = integrator.velocities
    dt = integrator.timestep
    m = integrator.mass_function
    f = integrator.force_function
    dist = integrator.distance_function
    fv = integrator.forces
    max_iter = integrator.max_iterations
    constraints = integrator.constraints
    tol = integrator.tolerance

    last_r = copy(r)

    a = acceleration(fv, m)
    Verlet._move_half_step!(r, v, a, dt)
    _move_half_step!(r, v, dt, tol, m, dist, max_iter, last_r, constraints)

    fv .= f(r)
    a = acceleration(fv, m)
    Verlet._move_full_step!(r, v, a, dt)
    _move_full_step!(r, v, tol, m, dist, max_iter, constraints)

    integrator, dt
end


text:
export AgentsIO

"""
    AgentsIO

Submodule containing functionality for serialization and deserialization of model data
to and from files.
"""
module AgentsIO
using Agents

include("csv_integration.jl")
include("jld2_integration.jl")
end