text:
module SolarLEDs

export main, spaceship_dj, sun_wind_switch

using LibSerialPort, Blink, Interact, COBS, ProgressMeter, Flatten, StaticArrays

#remove this
#=function get_port_list(;nports_guess::Integer=64)
    ports = sp_list_ports()
    port_list = String[]
    for port in unsafe_wrap(Array, ports, nports_guess, own=false)
        port == C_NULL && return port_list
        push!(port_list, sp_get_port_name(port))
    end
    sp_free_port_list(ports)
    return port_list
end=#

const N_UPLOAD_ATTEMPTS = 10
const BUTTON_LABELS = vcat(string.(0:9), "#", "*", "→", "↑", "←", "↓", "OK")
const N_BUTTONS = length(BUTTON_LABELS)
const BAUD = 115200

function tobytes(x::Integer)
    low = UInt8(x & 0xFF)
    high = UInt8((x >> 8) & 0xFF)
    return low, high
end
tobytes(x::Number) = tobytes(Int(x))
struct Sun
    int::UInt8
    low::UInt8
    high::UInt8
    siz::UInt8
end
function Sun(card, pos, rad, int, n_leds_per_strip)
    pos -= 1 + rad
    if card == :NS
        pos += n_leds_per_strip
    end
    low, high = tobytes(pos)
    return Sun(int, low, high, 2rad + 1)
end
function guiaxes(setup)
    card = radiobuttons(setup.cardinals)
    pos = slider(setup.elevations, value = first(setup.elevations))
    rad = slider(setup.radii, value = first(setup.radii))
    int = slider(setup.intensities, value = first(setup.intensities))
    on(pos) do p
        a = p - 1
        if rad[] > a
            rad[] = a
        end
        a = setup.n_leds_per_strip - p[]
        if rad[] > a
            rad[] = a
        end
    end
    on(rad) do r
        a = r + 1
        if pos[]  < a
            pos[] = a
        end
        a = setup.n_leds_per_strip - r
        if pos[] > a
            pos[] = a
        end
    end
    output = map(Sun, card, pos, rad, int, setup.n_leds_per_strip)
    wdg = Widget{:sun}(["card" => card, "pos" => pos, "rad" => rad, "int" => int], output = output)
    @layout! wdg vbox( pad(1em, hbox("Cardinal axes", :card)), pad(1em, hbox("Position", :pos)), pad(1em, hbox("Radius", :rad)), pad(1em, hbox("Intensity", :int)))
end

function uploaded(sp, msg)
    encode(sp, msg)
    @showprogress 0 "Uploading..." for i in 1:100
        if bytesavailable(sp) > 0
            return true
        else
            sleep(0.1)
        end
    end
    return false
end
function attemptupload(sp, msg)
    for attempt = 1:N_UPLOAD_ATTEMPTS
        if uploaded(sp, msg)
            msg2 = decode(sp)
            if msg2 == msg
                return nothing
            end
        end
        println("attempt #$attempt to upload data failed...")
    end
    error("failed to upload data after $N_UPLOAD_ATTEMPTS attempts")
end

struct Setup
    nsuns::Int
    n_leds_per_strip::Int
    cardinals::Vector{Symbol}
    elevations::Vector{Int}
    radii::Vector{Int}
    intensities::Vector{Int}
end

function main(; n_leds_per_strip::Int = 150, cardinals = [:EW, :NS], elevations = 1:150, radii = 0:25, intensities = 0:255)
    if 0 ∉ intensities
        push!(intensities, 0)
    end
    if !issorted(elevations)
        sort!(elevations)
    end
    if !issorted(radii)
        sort!(radii)
    end
    if !issorted(intensities)
        sort!(intensities)
    end
    nsuns = 4
    @assert 0 < nsuns "number of suns must be larger than zero"
    @assert 0 < n_leds_per_strip "number of LEDs per strip must be larger than zero"
    main(Setup(nsuns, n_leds_per_strip, cardinals, elevations, radii, intensities))
end

sun_wind_switch() = main(n_leds_per_strip = 150, elevations = [5, 17, 36, 48, 59, 63, 65, 66, 68, 69, 71, 73, 74, 76, 77, 79, 83, 94, 106, 125, 137], radii = [0], intensities = [255])
spaceship_dj() = main(n_leds_per_strip = 73, elevations = [12, 29, 46, 64], radii = [0,1,2], intensities = 0:255, cardinals = [:NA])

function returnport()
    for port in get_port_list()
        sp = LibSerialPort.sp_get_port_by_name(port)
        if occursin(r"arduino", LibSerialPort.sp_get_port_usb_manufacturer(sp))
            LibSerialPort.sp_free_port(sp)
            return port
        end
        LibSerialPort.sp_free_port(sp)
    end
    return nothing
end

function main(setup)

    port = returnport()
    isnothing(port) && error("no ports were detected...")
    serialport = open(port, BAUD)

    bottoms = Dict(l => guisuns(serialport, setup) for l in BUTTON_LABELS)
    top = tabs(BUTTON_LABELS)
    bottom = map(top) do l
        bottoms[l]
    end
    #=download = button("Download")
    on(download) do _
    l = top[]
    i = UInt8(findfirst(isequal(l), BUTTON_LABELS) - 1)
    uploaded([0x02, i])
    msg = decode(serialport)
    setgui(bottoms[l], msg)
    # bottom[] = bottom[]
    end=#
    upload = button("Upload")
    on(upload) do _
        l = top[]
        i = UInt8(findfirst(isequal(l), BUTTON_LABELS) - 1)
        msg = vcat(0x01, i, bottom[][])
        attemptupload(serialport, msg)
    end
    uploadall = button("Upload all")
    on(uploadall) do _
        for (i,l) in enumerate(BUTTON_LABELS)
            b = bottoms[l]
            msg = vcat(0x01, i - 1, b[])
            attemptupload(serialport, msg)
        end
    end
    reset = button("Reset")
    on(reset) do _
        flush(serialport)
        uploaded(serialport, [UInt8(3)])
        @assert Bool(decode(serialport)[])
    end
    w = Window()
    body!(w, dom"div"(hbox(pad(1em, uploadall), pad(1em, upload), pad(1em, reset)), top, bottom))

    return nothing
end

function guisuns(sp, setup)
    suns = [guiaxes(setup) for i in 1:setup.nsuns]
    output = map(suns...) do ss...
        SVector(flatten(ss))
    end
    on(output) do pl
        encode(sp, [0x00; pl])
    end
    d = Dict("sun$i" => s for (i,s) in enumerate(suns))
    wgt = Widget{:suns}(d, output = output)
    @layout! wgt hbox(:sun1, :sun2, :sun3, :sun4)
end

### TODO
# implement downloading
# maybe prepare hosting this to a server for mobile interface

end # module


text:
#=
bachelorbench:
- Julia version:
- Author: Dan
- Date: 2021-10-31
=#

push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL

pdb_filepath = "1AKE.pdb"

comp = System(pdb_filepath)

function counterComp()
    res = collectResidues(comp)
    filter!(x->x.name_ == "ALA", res)
end

counterComp()
elapsed = @elapsed counterComp()

println(elapsed)

ext: jl


text:
"""
    mergeSubtrees(original::Union{CGPGenotype, STGPGenotype}, subtree::Array{Node},
                  pointIndex::Integer, removedLen::Integer)

Merges two subtrees at the given merging point.
"""
function mergeSubtrees(original::Union{CGPGenotype, STGPGenotype}, subtree::Array{Node},
                       pointIndex::Integer, removedLen::Integer)

    originalTree = original._representation
    originalLen = length(originalTree)
    subtreeLen = length(subtree)
    newLen = originalLen - removedLen + subtreeLen
    childTree = Array{Node}(undef, newLen)

    for i=1:newLen
        if i < pointIndex
            childTree[i] = originalTree[i]
        elseif pointIndex <= i < (subtreeLen+pointIndex)
            childTree[i] = subtree[i-pointIndex+1]
        else
            childTree[i] = originalTree[i+(removedLen-subtreeLen)]
        end
    end

    return childTree
end

"""
    chooseAnotherNode(node::Node, terminalSet::Array{TerminalNode},
                      functionSet::Array{FunctionNode}, rng::Random.AbstractRNG,
                      allowedType::Union{DataType, Union} = Any)

Chooses, if possible, another node of the same arity and type of the given one.
Returns `nothing` if such node does not exist.
"""
function chooseAnotherNode(node::Node, terminalSet::Array{TerminalNode},
                           functionSet::Array{FunctionNode}, rng::Random.AbstractRNG,
                           allowedType::Union{DataType, Union} = Any)
    local newSet
    if allowedType == Any
        if typeof(node) <: TerminalNode
           newSet = filter(x->x!=node, terminalSet)
        elseif typeof(node) <: FunctionNode
           arity = getArity(node)
           newSet = filter(x->getArity(x)==arity && x!=node, functionSet)
        end
    else
        if typeof(node) <: TerminalNode
           newSet = filter(x->getType(x)<:allowedType && x!=node, terminalSet)
        elseif typeof(node) <: FunctionNode
           arity = getArity(node)
           newSet = filter(x->getType(x)<:allowedType && getArity(x)==arity && x!=node, functionSet)
        end
    end

    return !isempty(newSet) ? rand(rng, newSet) : nothing
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using Random; Random.seed!(1234);

include("src/data.jl")
using .DataUtil
include("src/stepwise/SelectiveInference.jl")
using .SelectiveInference
include("src/ci.jl")
using .ConfidenceInterval
using LinearAlgebra:I
using Distributions
using IntervalSets
using Parameters
using ProgressMeter
using Plots, StatsPlots, Statistics, DataFrames

function selective_CI_rand(test_prob::SelectiveHypothesisTest; α=0.05, isOC=false)
    @unpack alg, alg_result, test_stats, 𝐗, 𝐲_obs, 𝚺 = test_prob
    M_obs = collect(alg_result.model.M)
    isempty(M_obs) && return nothing
    test_stat = test_stats[rand(1:length(M_obs))]
    z_obs, σ²_z = test_stat.z_obs, test_stat.σ²_z

    # compute Z = {z | selection_event}
    z_min = -abs(z_obs) - 10 * √(σ²_z); z_max = abs(z_obs) + 10 * √(σ²_z)
    Z = isOC ? SelectiveInference.region(𝐗, 𝐲_obs, 𝚺, test_stat, alg)[2] ∩ [z_min..z_max] : SelectiveInference.truncated_interval(alg_result.model, 𝐗, 𝚺, test_stat, alg; z_min=z_min, z_max=z_max)

    return confidence_interval(z_obs, √(σ²_z), Z; α=α)
end

# CI experiment
function CI_experiment(param::DataParam, alg::StepwiseFeatureSelection; α=0.05, isOC=false, iter_size=5)
    CI_lengths = fill(NaN, iter_size)
    @showprogress for i = 1:iter_size
        try
            𝐗, 𝐲_obs = make_dataset(param)
            alg_result = stepwise(𝐗, 𝐲_obs, I, alg)
            test_prob = SelectiveHypothesisTest(alg, alg_result, 𝐗, 𝐲_obs, I)
            CI = selective_CI_rand(test_prob; α=α, isOC=isOC)
            isnothing(CI) && continue
            CI_lengths[i] = CI.right - CI.left
        catch e
            println(e)
            continue
        end
    end
    return CI_lengths
end

using JLD
function main(alg::StepwiseFeatureSelection)
    # parameters
    N = 100
    β₀ = 0.0
    𝛃 = randn(10)
    α = 0.05
    param = DataParam(N=N, β₀=β₀, 𝛃=𝛃)

    iter_size = 100
    CI_lengths = CI_experiment(param, alg; α=α, isOC=false, iter_size=iter_size)
    CI_lengths_OC = CI_experiment(param, alg; α=α, isOC=true, iter_size=iter_size)

    # plot result
    p = plot()
    violin!(p, fill(1, length(CI_lengths)), CI_lengths; c=:red2, label=nothing)
    violin!(p, fill(2, length(CI_lengths_OC)), CI_lengths_OC; c=:blue2, label=nothing)
    plot!(p; xticks=(1:2, ["homotopy", "over-conditioning"]))
    plot!(p; yscale=:log10)
    plot!(p; xlabel="method", ylabel="length of confidence interval")
    plot!(p; guidefontsize=16, legendfontsize=12, tickfontsize=10)
    savefig(p, "CI_length_forward_backward.pdf")
end

main(StepwiseFeatureSelection())


text:
using Pkg
import REPL
using REPL.TerminalMenus

using DataStructures
using HTTP
using JSON

# TOPICS.txt parser
function parsetopics(topicstxt)
    # category -> topic
    # e.g. Basic Concepts -> abstract_classes
    topics = OrderedDict{String, Vector{String}}()

    category = ""
    for (i, l) in enumerate(topicstxt)
        if isempty(category)
            category = l
            topics[category] = String[]
            continue
        elseif isempty(l)
            # new line means that the next line will be a new heading
            category = ""
            continue
        elseif all(c -> c == '-', l)
            # underline, ignore it
            continue
        end

        push!(topics[category], l)
    end

    topics
end

# [Yes/No] selector
# The following function has been derived from BinaryBuilder.jl's `yn_prompt`:
# https://github.com/JuliaPackaging/BinaryBuilder.jl/blob/60124d0532ae689227d6c707ab899841e6306c10/src/wizard/utils.jl#L118-L141
# BinaryBuilder.jl license info: Copyright (c) 2017: Elliot Saba [MIT License]
# Full license text: https://github.com/JuliaPackaging/BinaryBuilder.jl/blob/60124d0532ae689227d6c707ab899841e6306c10/LICENSE.md
function yn_select(question, default = :y)
    @assert default in (:y, :n)
    ynstr = default == :y ? "[Y/n]" : "[y/N]"
    while true
        print(question, " ", ynstr, ": ")
        answer = lowercase(strip(readline()))
        if isempty(answer)
            return default == :y
        elseif answer == "y" || answer == "yes"
            return true
        elseif answer == "n" || answer == "no"
            return false
        else
            println("Unrecognized answer. Answer `y` or `n`.")
        end
    end
end

# slug
function choose_slug()
    print("exercise slug: ")
    readline()
end

# core
iscore(slug) = yn_select("Is $slug a core exercise?", :n)

# unlocked_by
function choose_unlocked_by(slug, core, core_exercises)
    if !core
        menu = RadioMenu(core_exercises)
        choice = request("choose the core exercise that $slug is unlocked by:", menu)
        core_exercises[choice]
    else
        nothing
    end
end

# difficulty
function choose_difficulty(slug)
    while true
        print("choose the difficulty of $slug [1-10]: ")
        d = parse(Int, readline())
        if d ∈ 1:10
            return d
        else
            @error "difficulty must be a number from 1-10"
        end
    end
end

function shared_topics()
    url = "https://raw.githubusercontent.com/exercism/problem-specifications/master/TOPICS.txt"
    r = HTTP.request("GET", url)
    split(String(r.body), '\n')
end

# topics
function choose_topics(slug)
    println("Please select all of the following topics that are covered by this exercise.")
    alltopics = merge(parsetopics(shared_topics()), parsetopics(readlines(open("TOPICS.txt"))))
    topics = String[]
    for category in keys(alltopics)
        # MultiSelectMenu requires at least two options
        # if a category only has one entry, add an empty entry
        if length(alltopics[category]) < 2
            push!(alltopics[category], "")
        end

        menu = MultiSelectMenu(alltopics[category])
        choices = request("$category:", menu)
        for i in choices
            push!(topics, alltopics[category][i])
        end
        println()
    end
    topics
end

const configlet = joinpath("bin", "configlet")

# check if configlet has been installed and the script is running from the root dir of the repo
try
    read(`$configlet version`)
catch e
    isa(e, Base.IOError) && error("could not run configlet\n\nMake sure to run the script from the root directory of the track repo\nand that the configlet has been fetched using `bin/fetch-configlet`.\n\n")
end

# load config.json
const cfg = JSON.parsefile("config.json")

# load exercises from config
const core_exercises = [ex["slug"] for ex in cfg["exercises"] if ex["core"]]

# show prompts
const uuid = strip(String(read(`$configlet uuid`)))
const slug = choose_slug()
const core = iscore(slug)
const unlocked_by = choose_unlocked_by(slug, core, core_exercises)
const difficulty = choose_difficulty(slug)
const topics = choose_topics(slug)

const cfg_entry = OrderedDict(
    "uuid" => uuid,
    "slug" => slug,
    "core" => core,
    "unlocked_by" => unlocked_by,
    "difficulty" => difficulty,
    "topics" => topics,
)

println(json(cfg_entry))


text:
struct Device
    handle::nvmlDevice_t

    function Device(index::Integer)
        handle_ref = Ref{nvmlDevice_t}()
        nvmlDeviceGetHandleByIndex_v2(index, handle_ref)
        return new(handle_ref[])
    end

    function Device(uuid::Base.UUID)
        handle_ref = Ref{nvmlDevice_t}()
        nvmlDeviceGetHandleByUUID("GPU-$uuid", handle_ref)
        return new(handle_ref[])
    end

    function Device(serial::String)
        handle_ref = Ref{nvmlDevice_t}()
        nvmlDeviceGetHandleBySerial(serial, handle_ref)
        return new(handle_ref[])
    end
end

Base.unsafe_convert(::Type{nvmlDevice_t}, dev::Device) = dev.handle

# iteration

struct DeviceIterator end

devices() = DeviceIterator()

Base.eltype(::DeviceIterator) = Device

function Base.iterate(iter::DeviceIterator, i=1)
    i >= length(iter) + 1 ? nothing : (Device(i-1), i+1)
end

function Base.length(::DeviceIterator)
    count_ref = Ref{Cuint}()
    nvmlDeviceGetCount_v2(count_ref)
    return count_ref[]
end

Base.IteratorSize(::DeviceIterator) = Base.HasLength()

# queries

function name(dev::Device)
    buf = Vector{Cchar}(undef, NVML_DEVICE_NAME_BUFFER_SIZE)
    nvmlDeviceGetName(dev, pointer(buf), length(buf))
    return unsafe_string(pointer(buf))
end

function brand(dev::Device)
    ref = Ref{nvmlBrandType_t}()
    nvmlDeviceGetBrand(dev, ref)
    return ref[]
end

function uuid(dev::Device)
    buf = Vector{Cchar}(undef, NVML_DEVICE_UUID_BUFFER_SIZE)
    nvmlDeviceGetUUID(dev, pointer(buf), length(buf))
    uuid_str = unsafe_string(pointer(buf))
    @assert startswith(uuid_str, "GPU-")
    return Base.UUID(uuid_str[5:end])
end

function serial(dev::Device)
    buf = Vector{Cchar}(undef, NVML_DEVICE_SERIAL_BUFFER_SIZE)
    nvmlDeviceGetSerial(dev, pointer(buf), length(buf))
    return unsafe_string(pointer(buf))
end

# watt
function power_usage(dev::Device)
    ref = Ref{Cuint}()
    nvmlDeviceGetPowerUsage(dev, ref)
    return ref[] / 1000
end

# joules
function energy_consumption(dev::Device)
    ref = Ref{Culonglong}()
    nvmlDeviceGetTotalEnergyConsumption(dev, ref)
    return ref[] / 1000
end

# bytes
function memory_info(dev::Device)
    ref = Ref{nvmlMemory_t}()
    nvmlDeviceGetMemoryInfo(dev, ref)
    return (total=Int(ref[].total), free=Int(ref[].free), used=Int(ref[].used))
end

function compute_mode(dev::Device)
    ref = Ref{nvmlComputeMode_t}()
    nvmlDeviceGetComputeMode(dev, ref)
    return ref[]
end

function compute_capability(dev::Device)
    major = Ref{Cint}()
    minor = Ref{Cint}()
    nvmlDeviceGetCudaComputeCapability(dev, major, minor)
    return VersionNumber(major[], minor[])
end

function compute_processes(dev::Device)
    count_ref = Ref{Cuint}(0)
    res = unsafe_nvmlDeviceGetComputeRunningProcesses(dev, count_ref, C_NULL)
    if res == NVML_SUCCESS
        return nothing
    elseif res !== NVML_ERROR_INSUFFICIENT_SIZE
        throw_api_error(res)
    end

    infos = Vector{nvmlProcessInfo_t}(undef, count_ref[])
    nvmlDeviceGetComputeRunningProcesses(dev, count_ref, infos)

    return Dict(Int(info.pid) => (used_gpu_memory=Int(info.usedGpuMemory),)
                for info in infos)
end


text:
using JLD

include("instruction.jl")
include("map.jl")

function build_dict(instructions)
    d = Dict{AbstractString, Int}()

    for ins in instructions
        for w in ins.text
            for t in split(w, "-")
                get!(d, t, 1+length(d))
            end
        end
    end

    return d
end

function build_char_dict(instructions)
    d = Dict{Char, Int}()
    get!(d, ' ', 1+length(d))

    for ins in instructions
        for w in ins.text
            for t in split(w, "-")
                for c in t
                    get!(d, c, 1+length(d))
                end
            end
        end
    end

    return d
end

#converts tokens to onehots
function ins_arr(d, ins)
    arr = Any[]
    vocablength = length(d)

    for w in ins
        for t in split(w, "-")
            indx = haskey(d, t) ? d[t] : vocablength + 1
            onehot = zeros(Float32, 1, vocablength + 1)
            onehot[1, indx] = 1.0
            push!(arr, onehot)
        end
    end

    return arr
end

#converts tokens to embeddings
function ins_arr_embed(embeds, d, ins)
    arr = Any[]
    vocablength = length(d)

    for w in ins
        for t in split(w, "-")
            vec = haskey(d, t) ? transpose(embeds[t]) : transpose(embeds["unk"])
            push!(arr, vec)
        end
    end

    return arr
end

#converts chars to onehots
function ins_char_arr(d, ins)
    arr = Any[]
    vocablength = length(d)

    for i=1:length(ins)
        t = i==length(ins) ? ins[i] : "$(ins[i]) "
        for c in t
            indx = haskey(d, c) ? d[c] : vocablength + 1
            onehot = zeros(Float32, 1, vocablength + 1)
            onehot[1, indx] = 1.0
            push!(arr, onehot)
        end
    end
    return arr
end

#=
builds the agent's view
agent centric
up: front
right hand side: right of the agent
down: back
left hand side: left of the agent
(20, 20) is the agent curent location and it is a node
neighbors of a node are edges
=#
function state_agent_centric(map, loc; vdims = [39, 39])
    #lfeatvec = length(Items) + length(Floors) + length(Walls) + 3
    lfeatvec = length(Items) + length(Floors) + length(Walls) + length(MapColors) + 3
    view = zeros(Float32, vdims[1], vdims[2], lfeatvec, 1)
    mid = [round(Int, vdims[1]/2), round(Int, vdims[2]/2)]

    if loc[3] == 0
        ux = 0; uy = -1;
        rx = 1; ry = 0;
        dx = 0; dy = 1;
        lx = -1; ly = 0;
    elseif loc[3] == 90
        ux = 1; uy = 0;
        rx = 0; ry = 1;
        dx = -1; dy = 0;
        lx = 0; ly = -1;
    elseif loc[3] == 180
        ux = 0; uy = 1;
        rx = -1; ry = 0;
        dx = 0; dy = -1;
        lx = 1; ly = 0;
    else
        ux = -1; uy = 0;
        rx = 0; ry = -1;
        dx = 1; dy = 0;
        lx = 0; ly = 1;
    end

    current = loc[1:2]

    i, j = mid

    view[i, j, map.nodes[(loc[1], loc[2])]] = 1.0
    view[i,j, lfeatvec-2, 1] = 1.0

    i = i - 1
    #up
    while i > 0
        next = (current[1] + ux, current[2] + uy)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[i, j, length(Items) + floor, 1] = 1.0
            view[i, j, length(Items) + length(Floors) + wall, 1] = 1.0

            for c in ColorMapping[floor]
                view[i, j, length(Items) + length(Floors) + length(Walls) + MapColors[c], 1] = 1.0
            end

            view[i,j, lfeatvec-1, 1] = 1.0
            i = i - 1
            view[i, j, map.nodes[(next[1], next[2])], 1] = 1.0
            view[i,j, lfeatvec-2, 1] = 1.0
            current = next
            i = i - 1
        else
            #view[i,j, lfeatvec-1, 1] = 1.0
            view[i, j, lfeatvec, 1] = 1.0
            break
        end
    end

    current = loc[1:2]
    i, j = mid
    j = j + 1
    #right
    while j <= vdims[2]
        next = (current[1] + rx, current[2] + ry)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[i, j, length(Items) + floor, 1] = 1.0
            view[i, j, length(Items) + length(Floors) + wall, 1] = 1.0
            view[i,j, lfeatvec-1, 1] = 1.0

            for c in ColorMapping[floor]
                view[i, j, length(Items) + length(Floors) + length(Walls) + MapColors[c], 1] = 1.0
            end

            j = j + 1
            view[i, j, map.nodes[(next[1], next[2])], 1] = 1.0
            view[i,j, lfeatvec-2, 1] = 1.0
            current = next
            j = j + 1
        else
            #view[i,j, lfeatvec-1, 1] = 1.0
            view[i, j, lfeatvec, 1] = 1.0
            break
        end
    end

    current = loc[1:2]
    i, j = mid
    i = i + 1
    #down
    while i <= vdims[1]
        next = (current[1] + dx, current[2] + dy)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[i, j, length(Items) + floor] = 1.0
            view[i, j, length(Items) + length(Floors) + wall] = 1.0
            view[i,j, lfeatvec-1, 1] = 1.0

            for c in ColorMapping[floor]
                view[i, j, length(Items) + length(Floors) + length(Walls) + MapColors[c], 1] = 1.0
            end

            i = i + 1
            view[i, j, map.nodes[(next[1], next[2])]] = 1.0
            view[i,j, lfeatvec-2, 1] = 1.0
            current = next
            i = i + 1
        else
            #view[i,j, lfeatvec-1, 1] = 1.0
            view[i, j, lfeatvec] = 1.0
            break
        end
    end

    current = loc[1:2]
    i, j = mid
    j = j - 1
    #left
    while j > 0
        next = (current[1] + lx, current[2] + ly)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[i, j, length(Items) + floor] = 1.0
            view[i, j, length(Items) + length(Floors) + wall] = 1.0
            view[i,j, lfeatvec-1, 1] = 1.0

            for c in ColorMapping[floor]
                view[i, j, length(Items) + length(Floors) + length(Walls) + MapColors[c], 1] = 1.0
            end

            j = j - 1
            view[i, j, map.nodes[(next[1], next[2])]] = 1.0
            view[i,j, lfeatvec-2, 1] = 1.0
            current = next
            j = j - 1
        else
            #view[i,j, lfeatvec-1, 1] = 1.0
            view[i, j, lfeatvec] = 1.0
            break
        end
    end

    nv = zeros(Float32, 5, 20, size(view, 3), 1)
    nv[1, :, :, 1] = view[1:20, 20, :, :]
    nv[2, :, :, 1] = view[20, 20:end, :, :]
    nv[3, :, :, 1] = view[20:end, 20, :, :]
    nv[4, :, :, 1] = view[20, 1:20, :, :]
    nv[5, :, :, 1] = view[1:20, 20, :, :]
    #return view

    return nv
end

function state_agent_centric_multihot(map, loc)
    lfeatvec = length(Items) + length(Floors) + length(Walls) + 1
    view = zeros(Float32, 1, length(Items) + 4 * lfeatvec)

    if loc[3] == 0
        ux = 0; uy = -1;
        rx = 1; ry = 0;
        dx = 0; dy = 1;
        lx = -1; ly = 0;
    elseif loc[3] == 90
        ux = 1; uy = 0;
        rx = 0; ry = 1;
        dx = -1; dy = 0;
        lx = 0; ly = -1;
    elseif loc[3] == 180
        ux = 0; uy = 1;
        rx = -1; ry = 0;
        dx = 0; dy = -1;
        lx = 1; ly = 0;
    else
        ux = -1; uy = 0;
        rx = 0; ry = -1;
        dx = 1; dy = 0;
        lx = 0; ly = 1;
    end

    current = loc[1:2]

    view[1, map.nodes[(loc[1], loc[2])]] = 1.0
    walkable = false

    #up
    while true
        next = (current[1] + ux, current[2] + uy)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[1, length(Items) + length(Items) + floor] = 1.0
            view[1, length(Items) + length(Items) + length(Floors) + wall] = 1.0

            view[1, length(Items) + map.nodes[(next[1], next[2])]] = 1.0
            current = next
            walkable = true
        else
            if !walkable
                view[1, length(Items) + lfeatvec] = 1.0
            end
            break
        end
    end

    current = loc[1:2]
    walkable = false

    #right
    while true
        next = (current[1] + rx, current[2] + ry)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[1, length(Items) + lfeatvec + length(Items) + floor] = 1.0
            view[1, length(Items) + lfeatvec + length(Items) + length(Floors) + wall] = 1.0

            view[1, length(Items) + lfeatvec + map.nodes[(next[1], next[2])]] = 1.0
            current = next
            walkable = true
        else
            if !walkable
                view[1, length(Items) + 2*lfeatvec] = 1.0
            end
            break
        end
    end

    current = loc[1:2]
    walkable = false

    #down
    while true
        next = (current[1] + dx, current[2] + dy)
        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[1, length(Items) + 2*lfeatvec + length(Items) + floor] = 1.0
            view[1, length(Items) + 2*lfeatvec + length(Items) + length(Floors) + wall] = 1.0

            view[1, length(Items) + 2*lfeatvec + map.nodes[(next[1], next[2])]] = 1.0
            current = next
            walkable = true
        else
            if !walkable
                view[1, length(Items) + 3*lfeatvec] = 1.0
            end
            break
        end
    end

    current = loc[1:2]
    walkable = false

    #left
    while true
        next = (current[1] + lx, current[2] + ly)

        if haskey(map.edges[(current[1], current[2])], (next[1], next[2]))#check the wall existence
            wall, floor = map.edges[(current[1], current[2])][(next[1], next[2])]
            view[1, length(Items) + 3*lfeatvec + length(Items) + floor] = 1.0
            view[1, length(Items) + 3*lfeatvec + length(Items) + length(Floors) + wall] = 1.0

            view[1, length(Items) + 3*lfeatvec + map.nodes[(next[1], next[2])]] = 1.0
            current = next
            walkable = true
        else
            if !walkable
                view[1, length(Items) + 4*lfeatvec] = 1.0
            end
            break
        end
    end

    return view
end

function action(curr, next)
    a = 0
    if curr[1] != next[1] || curr[2] != next[2]#move
        a = 1
    elseif !(next[3] == 270 && curr[3] == 0) && (next[3] > curr[3] || (next[3] == 0 && curr[3] == 270))#right
        a = 2
    elseif !(next[3] == 0 && curr[3] == 270) && (next[3] < curr[3] || (next[3] == 270 && curr[3] == 0))#left
        a = 3
    else
        a = 4
    end
    return a
end

function getactions(path)
    actions = Any[]
    for i=1:length(path)
        curr = path[i]
        next = i == length(path) ? curr : path[i+1]
        push!(actions, action(curr, next))
    end
    return actions
end

function build_instance(instance::Instruction, map::Map, vocab; vdims=[39, 39], emb=nothing, encoding="grid")
    words = emb == nothing ? ins_arr(vocab, instance.text) : ins_arr_embed(emb, vocab, instance.text)

    states = Any[]
    Y = zeros(Float32, length(instance.path), 4)

    for i=1:length(instance.path)
        curr = instance.path[i]
        next = i == length(instance.path) ? curr : instance.path[i+1]
        Y[i, action(curr, next)] = 1.0
        push!(states, encoding == "grid" ? state_agent_centric(map, curr) : state_agent_centric_multihot(map, curr))
    end

    return (words, states, Y)
end

function minibatch(data; bs=100)
    sort!(data; by=t->length(t[1]))
    batches = []

    vocab = size(data[1][1][1], 2)
    vdims = size(data[1][2][1])

    for i=1:bs:length(data)
        l = i + bs - 1
        l = l > length(data) ? length(data) : l
        words = Any[]
        views = Any[]
        ys = Any[]
        maskouts = Any[]

        maxe = maximum(map(ind->length(data[ind][1]), i:l))
        maxd = maximum(map(ind->length(data[ind][2]), i:l))

        for enc=1:maxe
            word = zeros(Float32, (l-i+1), vocab)

            for j=i:l
                t = j-i+1

                if length(data[j][1]) >= enc
                    word[t, :] = data[j][1][enc]
                end
            end
            push!(words, word)
        end

        multihot = length(vdims) == 2

        for dec=1:maxd
            view = multihot ? zeros(Float32, (l-i+1), vdims[2]) : zeros(Float32, vdims[1], vdims[2], vdims[3], (l-i+1))
            y = zeros(Float32, (l-i+1), 4)
            maskout = ones(Float32, (l-i+1), 1)

            for j=i:l
                t = j-i+1
                if length(data[j][2]) >= dec
                    if !multihot
                        view[:, :, :, t] = data[j][2][dec]
                    else
                        view[t, :] = data[j][2][dec]
                    end

                    y[t, :] = data[j][3][dec, :]
                else
                    maskout[t, 1] = 0.0
                end
            end

            push!(views, view)
            push!(ys, y)
            push!(maskouts, maskout)
        end

        push!(batches, (words, views, ys, maskouts))
    end

    return batches
end

function build_data(trn_ins, tst_ins, outfile; charenc=false, encoding="grid")
    fname = "data/maps/map-grid.json"
    grid = getmap(fname)

    fname = "data/maps/map-jelly.json"
    jelly = getmap(fname)

    fname = "data/maps/map-l.json"
    l = getmap(fname)

    maps = Dict("Grid" => grid, "Jelly" => jelly, "L" => l)

    voc_ins = copy(trn_ins)
    append!(voc_ins, tst_ins)

    println("Building the vocab...")
    vocab = !charenc ? build_dict(voc_ins) : build_char_dict(voc_ins)

    println("Converting data...")
    trn_data = map(x -> build_instance(x, maps[x.map], vocab; encoding=encoding), trn_ins)

    println("Saving...")

    save(outfile, "vocab", vocab, "maps", maps, "data", trn_data)
    println("Done!")
end

max_stars_repo_name: ozanarkancan/navi-app

max_issues_repo_name: ozanarkancan/navi-app

max_forks_repo_name: ozanarkancan/navi-app


text:
using GrowthMaps, GeoData, ModelParameters, Dates, Test, StaticArrays
using Unitful: °C, K, hr, d, mol, cal

dimz = Y(10:10), X(100:10:110)

# Set up series data
stressdata = GeoArray.([[1. 2.], [1. 2.],
                        [2. 3.], [3. 4.], [2.5 3.5],
                        [4. 5.], [5. 6.],
                        [6. 7.], [6. 7.], [6. 7.]], Ref(dimz); name=:stress)
# TODO set up tempdata
tempdata = GeoArray.([[270. 280.], [270. 280.],
                      [270. 280.], [270. 280.], [270. 280.],
                      [270. 280.], [270. 280.],
                      [270. 280.], [270. 280.], [270. 280.]], Ref(dimz); name=:tempdata)

# Build a GeoSeries
stacks = [GeoStack(NamedTuple{(:stress, :temp)}((stressdata[i], tempdata[i]))) for i in 1:length(stressdata)]
timedim = (Ti([DateTime(2016, 1, 3, 9),
               DateTime(2016, 1, 6, 15),
               DateTime(2016, 2, 3, 10),
               DateTime(2016, 2, 3, 14),
               DateTime(2016, 2, 18, 10),
               DateTime(2016, 3, 3, 3),
               DateTime(2016, 3, 3, 8),
               DateTime(2016, 4, 3, 14),
               DateTime(2016, 4, 4, 10),
               DateTime(2016, 4, 16, 14)
              ]; mode=Sampled(; span=Regular(Hour(3)))),)
series = GeoSeries(stacks, timedim)
@test series[At(DateTime(2016, 1, 3, 9))][:stress] == [1. 2.]

@testset "mapgrowth float" begin
    # Set up models
    lowerthreshold = 5
    lowermortalityrate = -1
    lower = Layer(:stress, LowerStress(lowerthreshold, lowermortalityrate))

    upperthreshold = 5K
    uppermortalityrate = -1/K
    upper = Layer(:stress, K, UpperStress(upperthreshold, uppermortalityrate))

    # Lower
    output = mapgrowth(lower;
        series=series,
        tspan=DateTime(2016, 1, 3):Month(1):DateTime(2016, 4, 3)
    );

    # Test were are not touching the original arrays
    @test series[At(DateTime(2016, 1, 3, 9))][:stress] == [1. 2.]

    @test output[Ti(1)] == [-4.0 -3.0]
    @test output[Ti(2)] == [-2.5 -1.5]
    @test output[Ti(At(DateTime(2016, 3, 3)))] == [-0.5 0.0]
    @test output[Ti(At(DateTime(2016, 4, 3)))] == [0.0 0.0]

    @test typeof(dims(output)) <: Tuple{Y,X,Ti}
    @test length(val(dims(output, Ti))) == 4

    # Upper
    output = mapgrowth(upper;
        series=series,
        tspan=DateTime(2016, 1, 3):Month(1):DateTime(2016, 4, 3)
    );

    @test output[Ti(1)] == [0. 0.]
    @test output[Ti(2)] == [0. 0.]
    @test output[Ti(At(DateTime(2016, 3, 3)))] == [0.0 -0.5 ]
    @test output[Ti(At(DateTime(2016, 4, 3)))] == [-1.0 -2.0]

    # Lower and Uppera, in a Model with an extra period
    output = mapgrowth(Model((lower, upper));
        series=series,
        tspan=DateTime(2016, 1, 3):Month(1):DateTime(2016, 5, 3)
    );

    # Test were are still not touching the original arrays
    @test series[At(DateTime(2016, 1, 3, 9))][:stress] == [1. 2.]

    @test output[Ti(1)] == [-4.0 -3.0]
    @test output[Ti(2)] == [-2.5 -1.5]
    @test output[Ti(At(DateTime(2016, 3, 3)))] == [-0.5 -0.5]
    @test output[Ti(At(DateTime(2016, 4, 3)))] == [-1.0 -2.0]

    @test_logs (:warn,"No files found for the 1 month period starting 2016-05-03T00:00:00")
    mapgrowth(
        Model((lower, upper));
        series=series,
        tspan=DateTime(2016, 1, 3):Month(1):DateTime(2016, 5, 3)
    );
end

struct MultiRate <: RateModel end

@inline GrowthMaps.rate(l::MultiRate, x) = SA[x/K, 2x/K, 3x/K]

@testset "mapgrowth returning static array" begin
    model = Layer(:stress, K, MultiRate())

    output = mapgrowth(model;
        series=series,
        tspan=DateTime(2016, 1, 3):Month(1):DateTime(2016, 4, 3),
        initval=SA[0.0, 0.0, 0.0]
    );
    @test output[1, 1, 1] == SA[1.0, 2.0, 3.0]
    @test output[1, 2, 1] == SA[2.0, 4.0, 6.0]
end


text:
"""
    MIOTree_defaults(kwargs...)

Contains default MIOTree parameters, and modifies them with kwargs.
"""
function MIOTree_defaults(;kwargs...)
    d = Dict(:max_depth => 5,
        :cp => 1e-6,
        :hypertol => 0.005, # hyperplane separation tolerance
        :minbucket => 0.01,
        :regression_sparsity => :all,
        :hyperplane_sparsity => :all,
        :regression => false)
    set_params!(d; kwargs...)
    return d
end

"""
    $(TYPEDSIGNATURES)

Returns nodes of a dense binary tree starting from root.
"""
function generate_binary_tree(root::BinaryNode, max_depth::Int)
    levs = Dict{Int64, Array{BinaryNode}}(0 => [root])
    nodes = [root]
    for i = 1:max_depth
        idxs = levs[i-1][end].idx .+ collect(1:2^i)
        levs[i] = BinaryNode[]
        for j = 1:length(levs[i-1])
            lc = BinaryNode(idxs[2*j-1])
            leftchild(levs[i-1][j], lc)
            push!(levs[i], lc)
            push!(nodes, lc)
            rc = BinaryNode(idxs[2*j])
            rightchild(levs[i-1][j], rc)
            push!(levs[i], rc)
            push!(nodes, rc)
        end
    end
    return nodes
end

mutable struct MIOTree
    model::JuMP.Model
    root::BinaryNode
    params::Dict
    classes::Union{Nothing, Array{Any}}
    solver

    function MIOTree(solver; kwargs...)
        root = BinaryNode(1)
        mt = new(JuMP.Model(solver),
                 root,
                 MIOTree_defaults(),
                 nothing,
                 solver)
        for (key, val) in kwargs
            if key in keys(mt.params)
                set_param!(mt, key, val)
            else
                throw(ErrorException("Bad kwarg with key $(key) and value $(val) in MIOTree constructor."))
            end
        end
        return mt
    end
end

get_param(mt::MIOTree, sym::Symbol) = get_param(mt.params, sym)

set_param!(mt::MIOTree, sym::Symbol, val::Any) = set_param!(mt.params, sym, val)

set_params!(mt::MIOTree; kwargs...) = set_params!(mt.params; kwargs...)

"""
    $(TYPEDSIGNATURES)

Clones an MIOTree via deepcopy.
"""
function clone(mt::MIOTree)
    return deepcopy(mt)
end

""" Returns all BinaryNodes of MIOTree. """
allnodes(mt::MIOTree) = [mt.root, alloffspring(mt.root)...]

""" Returns all leaf BinaryNodes of MIOTree. """
allleaves(mt::MIOTree) = [nd for nd in allnodes(mt) if is_leaf(nd)]

"""
    $(TYPEDSIGNATURES)

Returns the leaf nodes in which the data X fall.
"""
function apply(mt::MIOTree, X::Matrix)
    vals = BinaryNode[] # TODO: initialize empty array instead, based on types of indices in MIOTree.
    for i = 1:size(X, 1)
        row = X[i,:]
        nd = mt.root
        while !is_leaf(nd)
            lhs = sum(nd.a .* row)
            if lhs ≤ nd.b
                nd = nd.left
            else
                nd = nd.right
            end
        end
        push!(vals, nd)
    end
    return vals
end

"""
    $(TYPEDSIGNATURES)

Makes predictions using a tree, based on data X.
"""
function predict(mt::MIOTree, X::Matrix)
    vals = [] # TODO: initialize empty array instead, based on types of labels in MIOTree.
    if get_param(mt, :regression)
        for i = 1:size(X, 1)
            row = X[i,:]
            nd = mt.root
            while !is_leaf(nd)
                lhs = sum(nd.a .* row)
                if lhs ≤ nd.b
                    nd = nd.left
                else
                    nd = nd.right
                end
            end
            push!(vals, nd.label[1] + sum(nd.label[2] .* row))
        end
    else
        for i = 1:size(X, 1)
            row = X[i,:]
            nd = mt.root
            while !is_leaf(nd)
                lhs = sum(nd.a .* row)
                if lhs ≤ nd.b
                    nd = nd.left
                else
                    nd = nd.right
                end
            end
            push!(vals, nd.label)
        end
    end
    return vals
end

apply(mt::MIOTree, X::DataFrame) = apply(mt, Matrix(X))
# TODO: improve this by making sure that the DataFrame labels are in the right order.

"""
    $(TYPEDSIGNATURES)

Returns the prediction accuracy of MIOTree.
For classification: misclassification accuracy.
For regression: R^2.
"""

function score(mt::MIOTree, X, Y)
    preds = predict(mt, X)
    if get_param(mt, :regression)
        return 1 - sum((preds .- Y).^2) / sum((preds .- sum(Y)/length(Y)*ones(length(Y))).^2)
    else
        return sum(preds .== Y)/length(Y)
    end
end

"""
    $(TYPEDSIGNATURES)

Returns the number of nonzero hyperplane coefficients of the MIOTree.
"""
function complexity(mt::MIOTree)
    scor = 0
    for node in allnodes(mt)
        if !is_leaf(node) && !isnothing(node.a)
            nonzeros = length(node.a) - sum(isapprox.(node.a, zeros(length(node.a))))
            scor += nonzeros
        end
    end
    return scor
end

"""
    $(TYPEDSIGNATURES)

Cleans all variables and constraints from a MIOTree.
"""
function clean_model!(mt::MIOTree)
    mt.model = JuMP.Model(mt.solver)
    return
end

"""
    $(TYPEDSIGNATURES)

Populates the nodes of the MIOTree using optimal solution of the MIO problem.
"""
function populate_nodes!(mt::MIOTree)
    termination_status(mt.model) == MOI.OPTIMAL ||
        throw(ErrorException("MIOTree must be trained before it can be pruned."))
    m = mt.model
    for nd in allnodes(mt)
        if !is_leaf(nd)
            aval = getvalue.(m[:a][nd.idx, :])
            if sum(isapprox.(aval, zeros(length(aval)); atol = 1e-10)) != length(aval)
                nd.a = aval
                nd.b = getvalue.(m[:b][nd.idx])
            end
        end
    end
    if get_param(mt, :regression)
        for lf in allleaves(mt)
            regr_coeffs = getvalue.(m[:beta][lf.idx, :])
            regr_const = getvalue(m[:beta0][lf.idx])
            set_classification_label!(lf, (regr_const, regr_coeffs))
        end
    else
        for lf in allleaves(mt) # Then populate the class values...
            class_values = [isapprox(getvalue.(m[:ckt][i, lf.idx]), 1; atol=1e-4) for i = 1:length(mt.classes)]
            if sum(class_values) == 1
                set_classification_label!(lf, mt.classes[findall(class_values)[1]])
            elseif sum(class_values) > 1
                throw(ErrorException("Multiple classes assigned to node $(lf.idx)."))
            end
        end
    end
    return
end

"""
    prune!(mt::MIOTree)

Prunes MIOTree depending on hyperplane coefficients and class values.
See ```populate_nodes''' for
how to populate nodes using optimal solution data.
"""
function prune!(mt::MIOTree)
    queue = allleaves(mt)
    while !isempty(queue) # A bottom-up approach for pruning.
        node = popfirst!(queue)
        par = node.parent
        if is_leaf(node) && isnothing(node.label) && !isnothing(node.parent)
            offspring = alloffspring(par) # TODO: speed this up
            labels = [of.label for of in offspring if !isnothing(of.label)]
            if length(labels) == 1
                delete_children!(par)
                set_classification_label!(par, labels[1])
            elseif length(labels) == 0
                delete_children!(par)
                push!(queue, par)
            end
        end
    end
    leaves = allleaves(mt)
    if any(isnothing(leaf.label) for leaf in leaves)
        @warn "A subset of leaves is unlabeled. Must label for correct predictions. "
    end
    return
end

""" Prunes the MIOTree from the root down, so that its offspring are garbage-collected. """
chop_down!(mt::MIOTree) = delete_children!(mt.root) # TODO: check garbage collection.

function generate_binary_tree(mt::MIOTree)
    isempty(alloffspring(mt.root)) || throw(ErrorException("The MIOTree must be ungrown/chopped down to generate a dense binary tree. "))
    generate_binary_tree(mt.root, get_param(mt, :max_depth))
    return
end


text:
# generic groupby, combine, apply
# grpCombApply(df, grp::String, comb::String, f::Function)
#		Args:
#			`df:: dataframe`
#			`grp:: String`: Column name to group by
# 			`comb:: String`: Column name to combine
#			`f:: Function`: Function to apply to the combined column
function grpCombApply(df, grp::String, comb::String, f::Function)
    dfg = groupby(df, grp);
    combine(dfg, comb => f);
end


text:
using Dates
using DataFrames
using Test

using Stonks: SchemaValidationError, UpdatableSymbol, to_dataframe, last_workday
using Stonks.Stores: FileStore, load, save, update
using Stonks.Stores:
  find_update_candidates,
  list_partition_nesting,
  list_nested_files,
  validate_schema,
  symbols_from_df
using Stonks.Models: AssetPrice, AssetInfo, ExchangeRate

include("test_utils.jl")

@testset "Store updates" begin
  dest = joinpath(@__DIR__, "data/test_stonks")
  symbols = ["AAPL", "MSFT"]

  ds = FileStore{AssetPrice}(;
    path=dest, ids=[:symbol], partitions=[:symbol], time_column="date"
  )

  @testset "Generate Symbol types from datastore ids" begin
    ref_date = last_workday() - Day(7)
    prices = to_dataframe(fake_price_data(7, ref_date, symbols))
    upd_smb = symbols_from_df(prices, ds.ids; time_column=ds.time_column)
    @test upd_smb == map(x -> (x, ref_date + Day(1), last_workday()), symbols)
  end

  @testset "Update candidates, updates available" begin
    df = to_dataframe(fake_price_data(7, last_workday() - Dates.Day(7), symbols))
    candidates = find_update_candidates(df, symbols; ids=ds.ids, time_column=ds.time_column)
    @test length(candidates) == length(symbols)
    @test isa(maximum(map(x -> x.from, candidates)), Date)
  end

  @testset "Update candidates, updates available + new symbols" begin
    df = to_dataframe(fake_price_data(7, last_workday() - Dates.Day(7), symbols))
    symbols_new = vcat(symbols, ["IBM", "GOOG"])
    candidates_new = find_update_candidates(
      df, symbols_new; ids=ds.ids, time_column=ds.time_column
    )
    @test length(candidates_new) == length(symbols_new)
    @test @chain filter(x -> ismissing(x.from), candidates_new) map(x -> x.ticker, _) ==
      ["IBM", "GOOG"]
  end

  @testset "Update candidates, no updates available" begin
    df = to_dataframe(fake_price_data(7, last_workday(), symbols))
    candidates = find_update_candidates(df, symbols; ids=ds.ids, time_column=ds.time_column)
    @test isempty(candidates)
  end

  @testset "Update candidates, updates from new tickers only" begin
    df = to_dataframe(fake_price_data(7, today(), symbols))
    symbols_new = vcat(symbols, ["IBM", "GOOG"])
    candidates = find_update_candidates(
      df, symbols_new; ids=ds.ids, time_column=ds.time_column
    )
    @test length(candidates) == 2
    @test map(x -> x.ticker, candidates) == ["IBM", "GOOG"]
  end

  @testset "Update candidates exchange_rate, updates available" begin
    symbols = ["EUR/USD", "USD/CAD"]
    df = to_dataframe(fake_exchange_data(7, last_workday() - Dates.Day(7), symbols))
    candidates = find_update_candidates(
      df, symbols; ids=["base", "target"], time_column=ds.time_column
    )
    @test length(candidates) == length(symbols)
    @test isa(maximum(map(x -> x.from, candidates)), Date)
    @test all(map(c -> !isnothing(c.fx_pair), candidates))
  end

  @testset "Update candidates exchange_rate, updates available + new_symbols" begin
    symbols = ["EUR/USD", "USD/CAD"]
    df = to_dataframe(fake_exchange_data(7, last_workday() - Dates.Day(7), symbols))
    symbols_new = vcat(symbols, ["USD/JPY"])
    candidates = find_update_candidates(
      df, symbols_new; ids=["base", "target"], time_column="date"
    )
    @test length(candidates) == length(symbols_new)
    @test @chain filter(x -> ismissing(x.from), candidates) map(x -> x.ticker, _) ==
      ["USDJPY=X"]
    @test all(map(c -> !isnothing(c.fx_pair), candidates))
  end

  @testset "Update candidates exchange_rate, updates from new symbols only" begin
    symbols = ["EUR/USD", "USD/CAD"]
    df = to_dataframe(fake_exchange_data(7, last_workday() - Dates.Day(7), symbols))
    symbols_new = ["USD/JPY", "EUR/CZK"]
    candidates = find_update_candidates(
      df, symbols_new; ids=["base", "target"], time_column="date"
    )
    @test length(candidates) == length(symbols_new)
    @test @chain filter(x -> ismissing(x.from), candidates) map(x -> x.ticker, _) ==
      ["USDJPY=X", "EURCZK=X"]
  end

  @testset "Update candidates exchange_rate, no updates available" begin
    symbols = ["EUR/USD", "USD/CAD"]
    df = to_dataframe(fake_exchange_data(7, today(), symbols))
    candidates = find_update_candidates(
      df, symbols; ids=["base", "target"], time_column="date"
    )
    @test isempty(candidates)
  end
end


text:
using Ipopt
using MathProgBase
using Random

struct Quadratic <: MathProgBase.AbstractNLPEvaluator
    V::Matrix{Float64}
    m::Vector{Float64}
    A::Matrix{Float64}
end

function MathProgBase.initialize( ::Quadratic, requested_features::Vector{Symbol} )
    unimplemented = setdiff( requested_features, [:Grad,:Jac] )
    if !isempty( unimplemented )
        error( "The following features aren't implemented: " * join( string.(unimplemented), ", " ) )
    end
end

MathProgBase.features_available( ::Quadratic ) = [:Grad]

MathProgBase.eval_f( q::Quadratic, x ) = (x' * 0.5 * q.V + q.m') * x

function MathProgBase.eval_g( q::Quadratic, g, x )
    g[:] = q.A * x
end

function MathProgBase.eval_grad_f( q::Quadratic, grad_f, x )
    grad_f[:] = q.V * x + q.m
end

MathProgBase.jac_structure( q::Quadratic ) =
    ((m,n) -> ([(i % m) + 1 for i in 1:m*n], vcat([fill(i,m) for i in 1:n]...)))(size(q.A)...)

function MathProgBase.eval_jac_g( q::Quadratic, J, x )
    J[:] = q.A
end

solver = IpoptSolver()

model = MathProgBase.NonlinearModel(solver)

Random.seed!(1)
S = rand( 3, 3 )
V = S'*S
m = rand(3)
A = rand(1,3)
b = rand(1)
q = Quadratic(V,m,A)
x0 = rand(3)
MathProgBase.loadproblem!(model, 3, 1, fill(-Inf,3), fill(Inf,3), fill(-Inf,1), b, :Min, q)
MathProgBase.setwarmstart!( model, rand(3) )
MathProgBase.optimize!(model)
@assert( MathProgBase.status(model) == :Optimal )
x = MathProgBase.getsolution(model)

@assert( maximum(abs.(x + q.V\q.m)) < 1e-8 )

if hasmethod( MathProgBase.freemodel!, Tuple{typeof(model)} )
    MathProgBase.freemodel!( model )
end

b = [-5.0]
model = MathProgBase.NonlinearModel(solver)
MathProgBase.loadproblem!(model, 3, 1, fill(-Inf,3), fill(Inf,3), fill(-Inf,1), b, :Min, q)
x0 = rand(3)
MathProgBase.setwarmstart!( model, x0 )
MathProgBase.optimize!(model)
@assert( MathProgBase.status(model) == :Optimal )
x = MathProgBase.getsolution(model)

@assert( q.A * x <= b .+ 1e-6 )

lambda = -((q.A * inv(q.V) * q.m + b)/(q.A * inv(q.V) * q.A'))[1]
@assert( maximum(abs.(x - inv(V)*(- q.m' - lambda * q.A)')) < 1e-6 )

if hasmethod( MathProgBase.freemodel!, Tuple{typeof(model)} )
    MathProgBsae.freemodel!( model )
end

# the following optimization takes many iterations
# I'm not sure why
#A = rand(2,3)
#b = rand(2)
#q = Quadratic( V, m, A )
#model = MathProgBase.NonlinearModel(solver)
#MathProgBase.loadproblem!(model, 3, 2, fill(-Inf,3), fill(Inf,3), fill(-Inf,2), b, :Min, q)
#x0 = rand(3)
#MathProgBase.setwarmstart!( model, x0
#MathProgBase.optimize!(model)

#if hasmethod( MathProgBase.freemodel!, Tuple{typeof(model)} )
#    MathProgBsae.freemodel!( model )
#end

# looseing the constraints does work
Ap = A[1:1,:]
bp = b[1:1]
q = Quadratic( V, m, Ap )
model = MathProgBase.NonlinearModel(solver)
MathProgBase.loadproblem!(model, 3, 1, fill(-Inf,3), fill(Inf,3), fill(-Inf,1), bp, :Min, q)
MathProgBase.setwarmstart!( model, rand(3) )
MathProgBase.optimize!(model)
x = MathProgBase.getsolution(model)
@assert( all(Ap * x .<= bp .+ 1e-6) )


text:
using Test
using DFTK
using DFTK: interpolate_density, transfer_blochwave, compute_transfer_matrix
using LinearAlgebra

include("testcases.jl")

@testset "Interpolation of density" begin
    lattice = Array{Float64}(I, 3, 3)
    N = 20
    f(n) = cos(2π*(n-1)/N)
    grid_in = (N, N, N)
    grid_out = (N+20, N-5, N+2)
    ρin = [f(i)*f(j)+f(k) for i = 1:N, j = 1:N, k = 1:N]

    # Test interpolation
    ρout = interpolate_density(ρin, grid_in, grid_out, lattice, lattice)
    ρin2 = interpolate_density(ρout, grid_out, grid_in, lattice, lattice)
    @test maximum(abs.(ρin2-ρin)) < .01

    # Test supercell
    ρout = interpolate_density(ρin, grid_in, grid_in, lattice, 2lattice)
    ρout1 = interpolate_density(ρout, grid_in, grid_in, 2lattice, 4lattice)
    ρout2 = interpolate_density(ρin, grid_in, grid_in, lattice, 4lattice)
    @test maximum(abs.(ρout2 - ρout1)) < sqrt(eps())

    ρout = interpolate_density(ρin, grid_in, grid_in, lattice, 2lattice)
    ρout1 = interpolate_density(ρout, grid_in, grid_out, 2lattice, 4lattice)
    ρout2 = interpolate_density(ρin, grid_in, grid_out, lattice, 4lattice)
    @test maximum(abs.(ρout2 - ρout1)) < .01
end

@testset "Transfer of blochwave" begin
    tol = 1e-7
    Ecut = 5

    model  = model_LDA(silicon.lattice, silicon.atoms, silicon.positions)
    kgrid  = [2, 2, 2]
    kshift = [1, 1, 1] / 2
    basis  = PlaneWaveBasis(model; Ecut, kgrid, kshift)

    ψ = self_consistent_field(basis; tol=tol, callback=identity).ψ

    ## Testing transfers from basis to a bigger_basis and backwards

    # Transfer to bigger basis then same basis (both interpolations are
    # tested then)
    bigger_basis = PlaneWaveBasis(model; Ecut=(Ecut + 5), kgrid, kshift)
    ψ_b = transfer_blochwave(ψ, basis, bigger_basis)
    ψ_bb = transfer_blochwave(ψ_b, bigger_basis, basis)
    @test norm(ψ - ψ_bb) < eps(eltype(basis))

    T    = compute_transfer_matrix(basis, bigger_basis)  # transfer
    Tᵇ   = compute_transfer_matrix(bigger_basis, basis)  # back-transfer
    Tψ   = [Tk * ψk   for (Tk, ψk)   in zip(T, ψ)]
    TᵇTψ = [Tᵇk * Tψk for (Tᵇk, Tψk) in zip(Tᵇ, Tψ)]
    @test norm(ψ - TᵇTψ) < eps(eltype(basis))

    # TᵇT should be the identity and TTᵇ should be a projection
    TᵇT = [Tᵇk * Tk  for (Tk, Tᵇk) in zip(T, Tᵇ)]
    @test all(M -> maximum(abs, M-I) < eps(eltype(basis)), TᵇT)

    TTᵇ = [Tk  * Tᵇk for (Tk, Tᵇk) in zip(T, Tᵇ)]
    @test all(M -> M ≈ M*M, TTᵇ)

    # Transfer between same basis (not very useful, but is worth testing)
    bigger_basis = PlaneWaveBasis(model; Ecut, kgrid, kshift)
    ψ_b = transfer_blochwave(ψ, basis, bigger_basis)
    ψ_bb = transfer_blochwave(ψ_b, bigger_basis, basis)
    @test norm(ψ-ψ_bb) < eps(eltype(basis))

end


text:
# -------------------- FIXED EFFECTS --------------------
abstract ModelEffect
abstract FixedEffect <: ModelEffect

function getFixedFormula(idf::DataFrame, df_name::AbstractString="df_data")
    v_intercept= idf[idf[:x] .== "(Intercept)" , :estimate][1]
    v_group1 = idf[idf[:x] .== "group" , :estimate][1]
    v_out = string(v_intercept)
    for fixedeff in setdiff(idf[:x], ["(Intercept)","group"])
        v_coef=  idf[idf[:x] .== fixedeff, :estimate][1]
        v_out *= " + ("*string(v_coef)*"*"*df_name*"[:"*fixedeff*"])"
    end
    #if length(logvar) > 0
    #    v_out*="+log( "*df_name*"[:"*logvar*"] + 1) "
    #end
    v_out2=v_out*"+"*string(v_group1)
    return v_out, v_out2
end

# JStack.getFixedFormula(JStack.readFixedFile(dfx,"occ"))

function genFixedTotals(feff::FixedEffect, df_data::DataFrame)
    model=feff.v_model
    src=feff.src
    #logvar=feff.logvar
    intercept= src[src[:x] .== "(Intercept)" , :estimate][1]
    group1 = src[src[:x] .== "group" , :estimate][1]
    coefs=Float64[]
    cols=Symbol[]
    for fixedeff in setdiff(src[:x], ["(Intercept)","group"])
        v_coef =  src[src[:x] .== fixedeff, :estimate][1]
        push!(coefs,v_coef)
        push!(cols,Symbol(fixedeff))
    end
    #if length(logvar) > 0
    #    push!(cols,Symbol(logvar))
    #end
    println(coefs)
    println(cols)

    function calc_feff(row::DataFrameRow, ftype::Int64)
        tot=intercept
        ccnt=length(row)
        #if length(logvar) > 0
        #    ccnt=ccnt-1
        #end
        for x in 1:ccnt
            tot=tot+(row[x]*coefs[x])
        end
        #if length(logvar) > 0
        #    tot=tot+log(row[Symbol(logvar)]+1)
       # end
        if ftype == 1
            tot=tot+group1
        end
        return tot
    end
    df_data[Symbol("pre_"*model*"_score0")] = map(x->calc_feff(x,0), eachrow(df_data[cols]) )
    df_data[Symbol("pre_"*model*"_score1")] = map(x->calc_feff(x,1), eachrow(df_data[cols]) )

    df_data[Symbol("pre_"*model*"_score0")]=map(Float64,df_data[Symbol("pre_"*model*"_score0")])
    df_data[Symbol("pre_"*model*"_score1")]=map(Float64,df_data[Symbol("pre_"*model*"_score1")])
    #df_data[Symbol("pre_occ_score0")]=eval(parse(mocc.feff.fmula0))
    #df_data[Symbol("pre_occ_score1")]=eval(parse(mocc.feff.fmula1))
end
#genFixedTotals(mocc.feff,df_data)
#df_data[[:pre_occ_score1,:xpre_occ_score1]]

type FOcc <: FixedEffect
    sdf::DataFrame
    src::DataFrame
    fmula0::AbstractString
    fmula1::AbstractString
    v_model::AbstractString
    #logvar::AbstractString
    function FOcc(df_data::DataFrame, idf::DataFrame) this=new(); this.v_model="occ"; #this.logvar="trps_pre_p1";
        this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model);
        this.src = idf
        this.fmula0, this.fmula1 = getFixedFormula(this.src)
        #println(this.fmula0,"\n\n",this.fmula1)
        #df_data[Symbol("pre_"*this.v_model*"_score0")]=eval(parse(this.fmula0))
        #df_data[Symbol("pre_"*this.v_model*"_score1")]=eval(parse(this.fmula1))
        this.sdf[:B] = idf[idf[:x] .== "group", :estimate][1]
        this.sdf[:SE] = idf[idf[:x] .== "group", :std_error][1]
        this.sdf[:P] = idf[idf[:x] .== "group", :pr_z_][1]
        genFixedTotals(this,df_data)
        return this
    end
end

type FDolOcc   <: FixedEffect
    sdf::DataFrame
    src::DataFrame
    fmula0::AbstractString
    fmula1::AbstractString
    v_model::AbstractString
    #logvar::AbstractString
    function FDolOcc(df_data::DataFrame, idf::DataFrame) this=new(); this.v_model="dolocc"; #this.logvar="dol_per_trip_pre_p1";
        this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model);
        this.src = idf
        this.fmula0, this.fmula1 = getFixedFormula(this.src)
        #println(this.fmula0,"\n\n",this.fmula1)
        #df_data[Symbol("pre_"*this.v_model*"_score0")]=eval(parse(this.fmula0))
        #df_data[Symbol("pre_"*this.v_model*"_score1")]=eval(parse(this.fmula1))
        this.sdf[:B] = idf[idf[:x] .== "group", :estimate][1]
        this.sdf[:SE] = idf[idf[:x] .== "group", :std_error][1]
        this.sdf[:P] = idf[idf[:x] .== "group", :pr_z_][1]
        genFixedTotals(this,df_data)
        return this
    end
end

type FPen   <: FixedEffect
    sdf::DataFrame
    src::DataFrame
    fmula0::AbstractString
    fmula1::AbstractString
    v_model::AbstractString
    #logvar::AbstractString
    function FPen(df_data::DataFrame, idf::DataFrame) this=new(); this.v_model="pen"; #this.logvar="buyer_pre_p1";
        this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model);
        this.src = idf
        this.fmula0, this.fmula1 = getFixedFormula(this.src)
        #println(this.fmula0,"\n\n",this.fmula1)
        #df_data[Symbol("pre_"*this.v_model*"_score0")]=eval(parse(this.fmula0))
        #df_data[Symbol("pre_"*this.v_model*"_score1")]=eval(parse(this.fmula1))
        this.sdf[:B] = idf[idf[:x] .== "group", :estimate][1]
        this.sdf[:SE] = idf[idf[:x] .== "group", :std_error][1]
        this.sdf[:P] = idf[idf[:x] .== "group", :pr_z_][1]
        genFixedTotals(this,df_data)
        return this
    end
end

type FDolHH
    sdf::DataFrame
    v_model::AbstractString
    function FDolHH() this=new(); this.v_model="dolhh"; this.sdf=pushSDFrow!(deepcopy(SDF),"Total Campaign",this.v_model); return this  end
end

function Base.show(io::IO, fdolhh::FDolHH)
    show(io,fdolhh.sdf)
end
function Base.show(io::IO, fixedeffect::FixedEffect)
    show(io,fixedeffect.sdf)
end

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
# export test_adjacency, test_laplacian, test_accessors, test_arithmetic, test_other
using ArnoldiMethod

@testset "Graph matrices" begin
    function converttest(T::Type, var)
        @test typeof(T(var)) == T
    end

    function constructors(mat)
        adjmat = CombinatorialAdjacency(mat)
        stochmat = StochasticAdjacency(adjmat)
        adjhat = NormalizedAdjacency(adjmat)
        avgmat = AveragingAdjacency(adjmat)
        return adjmat, stochmat, adjhat, avgmat
    end

    function test_adjacency(mat)
        adjmat, stochmat, adjhat, avgmat = constructors(mat)
        @test adjmat.D == vec(sum(mat, dims=1))
        @test adjmat.A == mat
        @test isa(sparse(mat), SparseMatrixCSC)
        @test isa(sparse(stochmat), SparseMatrixCSC)
        @test isa(sparse(adjhat), SparseMatrixCSC)
        @test isa(sparse(avgmat), SparseMatrixCSC)
        @test isa(convert(CombinatorialAdjacency, adjmat), CombinatorialAdjacency)
        @test isa(convert(CombinatorialAdjacency, avgmat), CombinatorialAdjacency)
        @test prescalefactor(adjhat) == postscalefactor(adjhat)
        @test postscalefactor(stochmat) == prescalefactor(avgmat)
        @test prescalefactor(adjhat) == postscalefactor(adjhat)
        @test prescalefactor(avgmat) == Noop()
    end

    function test_laplacian(mat)
        lapl = CombinatorialLaplacian(CombinatorialAdjacency(mat))
        @test typeof(lapl) <: Laplacian
        #constructors that work.
        @test adjacency(lapl).A == mat
        adj = adjacency(lapl)
        @test typeof(StochasticAdjacency(adj)) <: StochasticAdjacency
        @test typeof(NormalizedAdjacency(adj)) <: NormalizedAdjacency
        @test typeof(AveragingAdjacency(adj)) <: AveragingAdjacency

        @test typeof(adjacency(lapl)) <: CombinatorialAdjacency
        # converttest(SparseMatrix{Float64}, lapl)

        adjmat, stochmat, adjhat, avgmat = constructors(mat)
        @test typeof(adjacency(lapl))  <: CombinatorialAdjacency
        stochlapl = StochasticLaplacian(StochasticAdjacency(adjmat))
        @test typeof(adjacency(stochlapl)) <: StochasticAdjacency
        averaginglapl = AveragingLaplacian(AveragingAdjacency(adjmat))
        @test typeof(adjacency(averaginglapl)) <: AveragingAdjacency

        normalizedlapl = NormalizedLaplacian(NormalizedAdjacency(adjmat))
        @test typeof(adjacency(normalizedlapl)) <: NormalizedAdjacency
        @test !(typeof(adjacency(normalizedlapl)) <: CombinatorialAdjacency)

        #constructors that fail.
        @test_throws MethodError CombinatorialAdjacency(lapl)
        @test_throws MethodError StochasticLaplacian(lapl)
        @test_throws MethodError NormalizedLaplacian(lapl)
        @test_throws MethodError AveragingLaplacian(lapl)
        @test_throws MethodError convert(CombinatorialAdjacency, lapl)

        L = sparse(lapl)

        @test sum(abs, (sum(L, dims=1))) == 0
    end

    function test_accessors(mat, n)
        adjmat, stochmat, adjhat, avgmat = constructors(mat)
        dv = degrees(adjmat)
        @test degrees(StochasticLaplacian(stochmat)) == dv
        @test degrees(NormalizedLaplacian(adjhat)) == dv
        @test degrees(AveragingLaplacian(avgmat)) == dv

        for m in (adjmat, stochmat, adjhat, avgmat)
            @test degrees(m) == dv
            @test eltype(m) == eltype(m.A)
            @test size(m) == (n, n)
            #@fact length(m) --> length(adjmat.A)
        end
    end

    function test_arithmetic(mat, n)
        adjmat, stochmat, adjhat, avgmat = constructors(mat)
        lapl = CombinatorialLaplacian(adjmat)
        onevec = ones(Float64, n)
        v = adjmat * ones(Float64, n)
        @test sum(abs, (adjmat * onevec)) > 0.0
        @test sum(abs, ((stochmat * onevec) / sum(onevec))) ≈ 1.0
        @test sum(abs, (lapl * onevec)) == 0
        g(a) = sum(abs, (sum(sparse(a), dims=1)))

        @test g(lapl) == 0
        @test g(NormalizedLaplacian(adjhat)) > 1e-13
        @test g(StochasticLaplacian(stochmat)) > 1e-13

        @test eigs(adjmat, which=LR())[1][1] > 1.0
        @test eigs(stochmat, which=LR())[1][1] ≈ 1.0
        @test eigs(avgmat, which=LR())[1][1] ≈ 1.0
        @test eigs(lapl, which=LR())[1][1] > 2.0
        @test eigs(sparse(lapl), which=SR())[1][2] > 0.0
        @test eigs(sparse(lapl), which=SR())[1][1] < 1e-7

        lhat = NormalizedLaplacian(adjhat)
        @test eigs(lhat, which=LR())[1][1] < 2.0 + 1e-9
    end

    function test_other(mat, n)
        adjmat = CombinatorialAdjacency(mat)
        lapl = CombinatorialLaplacian(CombinatorialAdjacency(mat))
        @test size(lapl, 1) == n
        @test size(lapl, 2) == n
        @test size(lapl, 3) == 1

        @test_throws MethodError symmetrize(StochasticAdjacency(adjmat))
        @test_throws MethodError symmetrize(AveragingAdjacency(adjmat))
        @test !issymmetric(AveragingAdjacency(adjmat))
        @test !issymmetric(StochasticAdjacency(adjmat))
        @test_throws MethodError symmetrize(NormalizedAdjacency(adjmat)).A # --> adjmat.A

        begin
            @test CombinatorialAdjacency(mat) == CombinatorialAdjacency(mat)
            S = StochasticAdjacency(CombinatorialAdjacency(mat))
            @test S.A == S.A
            @test sparse(S) != S.A
            @test adjacency(S) == S.A
            @test NormalizedAdjacency(adjmat) != adjmat
            @test StochasticLaplacian(S) != adjmat
            @test_throws MethodError StochasticLaplacian(adjmat) # --> not(adjmat)
            @test !issymmetric(S)
        end
    end

    function test_symmetry(mat, n)
        adjmat = CombinatorialAdjacency(mat)
        lapl = CombinatorialLaplacian(CombinatorialAdjacency(mat))
        @test size(lapl, 1) == n
        @test size(lapl, 2) == n
        @test size(lapl, 3) == 1

        @test_throws MethodError symmetrize(StochasticAdjacency(adjmat))
        @test_throws MethodError symmetrize(AveragingAdjacency(adjmat))
        @test_throws MethodError symmetrize(NormalizedAdjacency(adjmat)).A # --> adjmat.A
        @test symmetrize(adjmat).A == adjmat.A
        # these tests are basically the code
        @test symmetrize(adjmat, :triu).A == triu(adjmat.A) + triu(adjmat.A)'
        @test symmetrize(adjmat, :tril).A == tril(adjmat.A) + tril(adjmat.A)'
        @test symmetrize(adjmat, :sum).A == adjmat.A + adjmat.A

        @test issymmetric(StochasticAdjacency(adjmat)) == false
        @test issymmetric(AveragingAdjacency(adjmat)) == false

        normalized_adjmat = NormalizedAdjacency(adjmat)
        @test issymmetric(normalized_adjmat) == issymmetric(normalized_adjmat.A)

        @test_throws ArgumentError symmetrize(adjmat, :fake)
    end

    function test_punchedmatrix(mat, n)
        adjmat = CombinatorialAdjacency(mat)
        ahatp  = PunchedAdjacency(adjmat)
        y = ahatp * perron(ahatp)
        @test dot(y, ahatp.perron) ≈ 0.0 atol = 1.0e-8
        @test sum(abs, y) ≈ 0.0 atol = 1.0e-8
        eval, evecs = eigs(ahatp, which=LM())
        @test eval[1] - (1 + 1.0e-8)  <= 0
        @test dot(perron(ahatp), evecs[:, 1]) ≈ 0.0 atol = 1e-8
        ahat = ahatp.A
        @test isa(ahat, NormalizedAdjacency)

        z = ahatp * perron(ahat)
        @test norm(z) ≈ 0.0 atol = 1e-8
    end

    n = 10

    mat = Float64.(sprand(Bool, n, n, 0.3))

    test_adjacency(mat)
    test_laplacian(mat)
    test_accessors(mat, n)

    mat = symmetrize(Float64.(sprand(Bool, n, n, 0.3)))
    test_arithmetic(mat, n)
    test_other(mat, n)
    test_symmetry(mat, n)
    test_punchedmatrix(mat, n)

    """Computes the stationary distribution of a random walk"""
    function stationarydistribution(R::StochasticAdjacency; kwargs...)
        er = eigs(R, nev=1, which=LR(); kwargs...)
        l1 = er[1][1]
        abs(l1 - 1) < 1e-8 || error("failed to compute stationary distribution") # TODO 0.7: should we change the error type to InexactError?
        p = real(er[2][:, 1])
        if p[1] < 0
            for i in 1:length(p)
                p[i] = -p[i]
            end
        end
        return p
    end

    function stationarydistribution(A::CombinatorialAdjacency; kwargs...)
        R = StochasticAdjacency(A)
        stationarydistribution(R; kwargs...)
    end

    # Random walk demo
    n = 100
    p = 16 / n
    M = sprand(n, n, p)
    M.nzval[:] .= 1.0
    A = CombinatorialAdjacency(M)
    sd = stationarydistribution(A; mindim=10)
    @test all(sd .>= 0)

    @testset "Noop" begin
        @testset "Noop broadcasted with * is identity" begin
            testobjects = [1, "string", 'c', cycle_graph(3), [1,2,3], sin, Noop()]
            for object in testobjects
                @test Noop() .* object === object
            end
        end

        @testset "Noop * is identity" begin
            testobjects = [1, "string", 'c', cycle_graph(3), [1,2,3], sin, Noop()]
            for object in testobjects
                @test Noop() * object === object
            end
        end
    end

end


text:
"""
    Layers

 Defines the `Layer` type for a substrate layer.
"""
module Layers

using ..PSSFSSLen  # For length units and ustrip only
using StaticArrays: SVector

export Layer, Gblock, TEorTM, TE, TM

@enum TEorTM TE=1 TM=2

"""
    Layer <: Any

An instance of the `Layer` type represents a single dielectric layer of the physical structure.
It contains the electrical properties of the dielectric layer. For layers not included
in a Gblock, an instance of Layer also specifies the periodicity (via the
reciprocal lattice vectors) and stores the mode constants for the Floquet modes present
in the layer.

    Layer(width::0mm, ϵᵣ=1.0, tanδ=0.0, μᵣ=1.0, mtanδ=0.0)
    Layer(width::0mm, epsr=1.0, tandel=0.0, mur=1.0, mtandel=0.0)

Create a `Layer` instance with the specified electrical properties.  All arguments
are optional keyword arguments with default values as shown above. They can be
supplied in any order.  Typically the first and last `Layer` in a composite FSS or PSSFSS
structure are generated by a call to `Layer()` to represent the semi-infinite vacuum
regions surrounding the structure.  Using the default `width` of `0mm` means that the phase
reference planes are located at the surfaces just adjacent to these semi-infinite end `Layer`s.

# Arguments
- `width`: The layer width (i.e. thickness) expressed as a
        [`Unitful`](https://github.com/timholy/Unitful.jl)
        length quantity. For convenience the following unit
        suffixes are exported by this module: `m`, `cm`, `mil`, `inch`,
        so one can specify, e.g., `width=20mil`.  Note that `width` can
        be negative for the first and/or final layer of the composite structure,
        which has the effect of shifting the phase reference plane towards
        the interior of the composite structure.  This is sometimes needed
        when interfacing with other programs, such as TEP file generation for
        Ticra's `GRASP` program.
- `ϵᵣ` or `epsr`: Relative permittivity of the dielectric.
- `tanδ` or `tandel`: Loss tangent (electrical) of the dielectric.
- `μᵣ` or `mur`: Relative permeability of the dielectric.
- `mtanδ` or `mtandel`: Loss tangent (magnetic) of the dielectric.

"""
mutable struct Layer
    name::String
    ϵᵣ::ComplexF64  # Relative permittivity
    μᵣ::ComplexF64  # Relative permeability
    user_width::Unitful.Length # Layer thickness in Unitful units
    width::Float64  # (thickness in meters without attached units)
    #  Mode indices, sorted in order of increasing real part, or decreasing
    #  imaginary part, of γ, the mode attenuation constant.  M is the x
    #  index, N is the y index, and P = TE or TM:
    P::Vector{TEorTM}
    M::Vector{Int}
    N::Vector{Int}
    β₁::SVector{2,Float64}  # Reciprocal lattice vector (1/m)
    β₂::SVector{2,Float64}  # Reciprocal lattice vector (1/m)
    # Transverse portion of modal propagation vectors (radians/meter):
    β::Vector{SVector{2,Float64}}
    γ::Vector{ComplexF64} # Modal attenuation constants (nepers/meter)
    Y::Vector{ComplexF64} # Modal admittances (Siemens)
    c::Vector{ComplexF64} # Modal normalization constants (volts/meter)
    tvec::Vector{SVector{2,Float64}} # Modal unit electric field vectors (unitless)

    # Interior constructor:
    function Layer(;name="Layer", width::Unitful.Length=0u"mm", ϵᵣ::Real=1.0, epsr::Real=1.0,
                   tanδ::Real=0.0, tandel::Real=0.0, μᵣ::Real=1.0, mur::Real=1.0,
                   mtanδ::Real=0.0, mtandel::Real=0.0)
        epsr ≠ 1.0 && (ϵᵣ=epsr)
        mur ≠ 1.0 && (μᵣ=mur)
        tandel ≠ 0.0 && (tanδ=tandel)
        mtandel ≠ 0.0 && (mtanδ=mtandel)
        cϵᵣ = ϵᵣ * complex(1.0, -tanδ)
        cμᵣ = μᵣ * complex(1.0, -mtanδ)
        new(name, cϵᵣ, cμᵣ, width, float(ustrip(u"m", width)),
            TEorTM[], Int[], Int[], SVector{2}([0.,0.]), SVector{2}([0.,0.]),
            SVector{2,Float64}[], ComplexF64[], ComplexF64[], ComplexF64[],
            SVector{2,Float64}[])
    end # function
end # struct

Base.:(==)(l1::Layer, l2::Layer) =
             all(f -> getfield(l1, f) == getfield(l2, f), 1:nfields(l1))

Base.show(io::IO, ::MIME"text/plain", l::Layer) =
    print(io, l.name, ": width=", l.user_width, ", ϵᵣ=", real(l.ϵᵣ), ", tanδ=", -imag(l.ϵᵣ)/real(l.ϵᵣ),
            ", μᵣ=", real(l.μᵣ), ", mtanδ=", -imag(l.μᵣ)/real(l.μᵣ), ", ", length(l.P), " modes" )

end # module


text:
"""
Various utility functions for using across models, controllers and views.
"""
module Helpers

using Genie, Genie.Router, URIParser, Genie.Loggers, HTTP, Genie.Flax, Genie.Sessions

export request, response, flash, wsclient, flash_has_message

"""
    request() :: HTTP.Request

Returns the `Request` object associated with the current HTTP request.
"""
function request() :: HTTP.Request
  request(@params)
end

"""
    request(params::Dict{Symbol,Any}) :: HTTP.Request

Returns the `Request` object associated with the current HTTP request.
"""
function request(params::Dict{Symbol,Any}) :: HTTP.Request
  if haskey(params, Genie.PARAMS_REQUEST_KEY)
    return params[Genie.PARAMS_REQUEST_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_REQUEST_KEY) key"
    log(msg, :err)
    error(msg)
  end
end

"""
    response() :: HTTP.Response

Returns the `Response` object associated with the current HTTP request.
"""
function response() :: HTTP.Response
  response(@params)
end

"""
    response(params::Dict{Symbol,Any}) :: HTTP.Response

Returns the `Response` object associated with the current HTTP request.
"""
function response(params::Dict{Symbol,Any}) :: HTTP.Response
  if haskey(params, Genie.PARAMS_RESPONSE_KEY)
    return params[Genie.PARAMS_RESPONSE_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_RESPONSE_KEY) key"
    log(msg, :err)
    error(msg)
  end
end

"""
    flash()

Returns the `flash` dict object associated with the current HTTP request.
"""
function flash()
  flash(@params)
end

"""
    flash(params::Dict{Symbol,Any})

Returns the `flash` dict object associated with the current HTTP request.
"""
function flash(params::Dict{Symbol,Any})
  if haskey(params, Genie.PARAMS_FLASH_KEY)
    return params[Genie.PARAMS_FLASH_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_FLASH_KEY) key"
    log(msg, :err)
    error(msg)
  end
end

"""
    flash(value::Any) :: Nothing

Stores `value` on the `flash`.
"""
function flash(value::Any) :: Nothing
  flash(value, @params)
end

"""
    flash(value::Any, params::Dict{Symbol,Any}) :: Nothing

Stores `value` on the `flash`.
"""
function flash(value::Any, params) :: Nothing
  Sessions.set!(Sessions.session(params), Genie.PARAMS_FLASH_KEY, value)
  params[Genie.PARAMS_FLASH_KEY] = value

  nothing
end

"""
    flash_has_message() :: Bool

Checks if there's any value on the flash storage
"""
function flash_has_message() :: Bool
  ! isempty(flash())
end

"""
    wsclient(params::Dict{Symbol,Any}) :: HTTP.WebSockets.WebSocket

Returns the `WebSocket` object associated with the current WS request.
"""
function wsclient(params::Dict{Symbol,Any}) :: HTTP.WebSockets.WebSocket
  if haskey(params, Genie.PARAMS_WS_CLIENT)
    return params[Genie.PARAMS_WS_CLIENT]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_WS_CLIENT) key"
    log(msg, :err)
    error(msg)
  end
end

end


text:
function StepFunction(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return StepFunction((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function value(obj::StepFunction, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end


text: incx = 1


text:
using PyPlot

d = 0:0.01:1

function kp(d, A...)
    s = d * (1-d)
    r = 0.0
    for (i, a) in enumerate(A)
        r += a * (1-d)^(i-1)
    end
    return s * r
end

fig, axs = subplots()
plot(d, kp.(d, 0.0, 1.0), "-")
plot(d, kp.(d, 1.0, 0.0), "--")
plot(d, kp.(d, 0.1, 1.0), "-.")
plot(d, kp.(d, 1.0, -1.0), ":")
gcf()


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule rdma_core_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("rdma_core")
JLLWrappers.@generate_main_file("rdma_core", UUID("69dc3629-5c98-505f-8bcd-225213cebe70"))
end  # module rdma_core_jll


text:
## Bessel
@registerN Bessel AbstractCCOP1 α

"""
    Bessel{α}

Implements the [Bessel](https://dlmf.nist.gov/18.34) polynomials, introduced by [Krall and Frink](https://www.ams.org/journals/tran/1949-065-01/S0002-9947-1949-0028473-1/S0002-9947-1949-0028473-1.pdf) (with `b=2`). The  case `a=2` corresponds to the
[Bessel](https://en.wikipedia.org/wiki/Bessel_polynomials) polynomials of Wikipedia. The Bessel  polynomials are not orthogonal over  a domain of the real  line, rather over an arbitray curve in the complex plane enclosing the  origin.  The weight  function is `ρ(x)=(2πi)^(-1)∑Γ(α)/Γ(α+n-1)(-β/x)^n`,   where `β=2`.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> 𝐐 = Rational{Int}
Rational{Int64}

julia> x = variable(Polynomial{𝐐})
Polynomials.Polynomial(x)

julia> [basis(Bessel{3//2, 𝐐}, i)(x) for i in 0:5]
6-element Vector{Polynomial{Rational{Int64}, :x}}:
 Polynomials.Polynomial(1//1)
 Polynomials.Polynomial(1//1 + 3//4*x)
 Polynomials.Polynomial(1//1 + 5//2*x + 35//16*x^2)
 Polynomials.Polynomial(1//1 + 21//4*x + 189//16*x^2 + 693//64*x^3)
 Polynomials.Polynomial(1//1 + 9//1*x + 297//8*x^2 + 1287//16*x^3 + 19305//256*x^4)
 Polynomials.Polynomial(1//1 + 55//4*x + 715//8*x^2 + 10725//32*x^3 + 182325//256*x^4 + 692835//1024*x^5)
```

"""
Bessel
export Bessel
basis_symbol(::Type{<:Bessel{α}}) where {α} = "Cᵅ"
Polynomials.domain(::Type{<:Bessel}) = Polynomials.Interval(0, Inf)

abcde(::Type{<:Bessel{α}}) where {α} = NamedTuple{(:a, :b, :c, :d, :e)}((1, 0, 0, α, 2))

function k1k0(::Type{P}, k::Int) where {α,P<:Bessel{α}}
    k < 0 && return zero(eltype(P)) / 1
    iszero(k) && return (one(eltype(P)) * α) / 2

    val = one(eltype(P))
    val *= (2k + α) * (2k + α - 1)
    val /= (k + α - 1) * 2
    val
end

norm2(::Type{<:Bessel{α}}, n) where {α} =
    -1^(n + α - 1) * Γ(1 + n) * 2^(α - 1) / (Γ(n + α - 2) * (2n + α - 1))

# From https://www.ams.org/journals/tran/1949-065-01/S0002-9947-1949-0028473-1/S0002-9947-1949-0028473-1.pdf we use
# kn = (n + α - 1)_n / 2^n not (n+α+1)_n/2^n
# Koepf suggests kn = (n + α + 1)_n/2^n
function leading_term(P::Type{<:Bessel{α}}, n::Int) where {α}
    one(eltype(P)) / 2^n * Pochhammer(n + α - 1, n)
end

weight_function(::Type{<:Bessel{α}}) where {α} =
    z -> (2π * im)^(-1) * z^(α - 2) * exp(-2 / z)
generating_function(::Type{<:Bessel}) = (t, x) -> error("XXX")
function classical_hypergeometric(::Type{<:Bessel{α}}, n, x) where {α}
    as = (-n, n + α - 1)
    bs = ()
    pFq(as, bs, -x / 2) #/ kn
end

## Overrides XXX fails wih 1 and 2
function B̃n(P::Type{<:Bessel{α}}, n::Int) where {α}
    val = one(eltype(P))
    (iszero(n) && α == 2) && return val
    val *= 2 * (α - 2)
    val /= (2 * n + α) * (2 * n + α - 2)
    val
end

function C̃n(P::Type{<:Bessel{α}}, n::Int) where {α}
    val = one(eltype(P))
    n == 1 && return -(val * 4) / (α^2 * (α + 1))

    val *= -4 * n * (n + α - 2)
    val /= (2 * n + α - 3) * (2 * n + α - 2)^2 * (2 * n + α - 1)
    val
end

#Bn(::Type{<:Bessel{1}}, n::Int, ::Type{S}) where {S} = error("α=1 is not correct")
#B̃n(P::Type{<:Bessel{2}}, ::Val{0}) where {α} =  one(eltype(P))  # 0  otherwise
#iszero(N) #?  one(eltype(P)) : zero(eltype(P)))
#C̃n(P::Type{<:Bessel{α}}, ::Val{1}) where {α} =  -(one(eltype(P))*4)/(α^2*(α + 1))

b̂̃n(::Type{<:Bessel{2}}, n::Int) = (one(eltype(P)) * 2) / (n * (2n + 2))
b̂̃n(::Type{<:Bessel{2}}, ::Val{0}) = one(eltype(P)) * Inf


text:
#### Worker Structs ####

struct SubproblemRecord
    problem::AbstractSubproblem
    branch_vars::Vector{VariableID}
    leaf_vars::Vector{VariableID}
    subproblem_callbacks::Vector{SubproblemCallback}
end

mutable struct WorkerRecord
    id::Int
    warm_start::Bool
    subproblems::Dict{ScenarioID,SubproblemRecord}
    lb_subproblems::Dict{ScenarioID,SubproblemRecord}
end

function WorkerRecord(id::Int)
    return WorkerRecord(id,
                        false,
                        Dict{ScenarioID,SubproblemRecord}(),
                        Dict{ScenarioID,SubproblemRecord}(),
                        )
end

#### Worker Functions ####

function worker_loop(id::Int,
                     input::RemoteChannel,
                     output::RemoteChannel,
                     )::Int

    running = true
    my_record = WorkerRecord(id)

    try

        while running
            msg = take!(input)
            running = process_message(msg, my_record, output)
        end

    catch e

        while isready(input)
            take!(input)
        end

        rethrow()

    finally
        ; # intentional no-op
    end

    return 0
end

function process_message(msg::Message,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool
    error("Worker unable to process message type: $(typeof(msg))")
    return false
end

function process_message(msg::Abort,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool
    return false
end

function process_message(msg::DumpState,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool
    put!(output, WorkerState(record.id, record))
    return true
end

function process_message(msg::Initialize,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool

    record.warm_start = msg.warm_start

    # Master process needs the variable map messages to proceed. Initial solves are
    # not needed until much later in the process. So get all variable maps back first
    # then start working on the initial solves and the objective augmentation

    _build_subproblems(output,
                       record,
                       msg.scenarios,
                       msg.scenario_tree,
                       msg.create_subproblem,
                       msg.create_subproblem_args,
                       msg.create_subproblem_kwargs,
                       msg.subproblem_callbacks
                       )

    _initial_solve(output, record)

    _penalty_parameter_preprocess(output, record, msg.r)

    return true
end

function process_message(msg::InitializeLowerBound,
                         record::WorkerRecord,
                         output::RemoteChannel,
                         )::Bool

    for scen in msg.scenarios
        lb_sub = msg.create_subproblem(scen,
                                       msg.create_subproblem_args...;
                                       msg.create_subproblem_kwargs...
                                       )
        var_map = report_variable_info(lb_sub, msg.scenario_tree)
        (branch_ids, leaf_ids) = _split_variables(msg.scenario_tree,
                                                  collect(keys(var_map))
                                                  )
        add_lagrange_terms(lb_sub, branch_ids)
        record.lb_subproblems[scen] = SubproblemRecord(lb_sub,
                                                       branch_ids,
                                                       leaf_ids,
                                                       SubproblemCallback[]
                                                       )
    end

    return true

end

function process_message(msg::PenaltyInfo,
                         record::WorkerRecord,
                         output::RemoteChannel,
                         )::Bool

    if !haskey(record.subproblems, msg.scen)
        error("Worker $(record.id) received penalty parameter values for scenario $(msg.scen). This worker was not assigned this scenario.")
    end

    sub = record.subproblems[msg.scen]

    add_ph_objective_terms(sub.problem,
                           sub.branch_vars,
                           msg.penalty)

    return true
end

function process_message(msg::Ping,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool
    put!(output, Ping())
    return true
end

function process_message(msg::Solve,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool

    if !haskey(record.subproblems, msg.scen)
        error("Worker $(record.id) received solve command for scenario $(msg.scen). This worker was not assigned this scenario.")
    end

    sub = record.subproblems[msg.scen]

    update_ph_terms(sub.problem, msg.w_vals, msg.xhat_vals)

    _execute_subproblem_callbacks(sub, msg.niter, msg.scen)

    start = time()
    sts = solve_subproblem(sub.problem)
    stop = time()

    if record.warm_start == true
        warm_start(sub.problem)
    end

    var_vals = report_values(sub.problem, sub.branch_vars)

    put!(output, ReportBranch(msg.scen,
                              sts,
                              objective_value(sub.problem),
                              stop - start,
                              var_vals)
         )

    return true

end

function process_message(msg::SolveLowerBound,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool

    if !haskey(record.lb_subproblems, msg.scen)
        error("Worker $(record.id) received lower-bound solve command for scenario $(msg.scen). This worker was not assigned this scenario.")
    end

    lb_sub = record.lb_subproblems[msg.scen]
    update_lagrange_terms(lb_sub.problem, msg.w_vals)

    start = time()
    sts = solve_subproblem(lb_sub.problem)
    stop = time()

    if record.warm_start == true
        warm_start(sub.problem)
    end

    put!(output, ReportLowerBound(msg.scen,
                                  sts,
                                  objective_value(lb_sub.problem),
                                  stop - start)
         )

    return true

end

function process_message(msg::ShutDown,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool

    for (s, sub) in pairs(record.subproblems)
        leaf_vals = report_values(sub.problem, sub.leaf_vars)
        put!(output, ReportLeaf(s, leaf_vals))
    end

    return false

end

function process_message(msg::SubproblemAction,
                         record::WorkerRecord,
                         output::RemoteChannel
                         )::Bool

    sub = record.subproblems[msg.scen].problem
    msg.action(sub, msg.args...; msg.kwargs...)

    return true

end

#### Internal Functions ####

function _build_subproblems(output::RemoteChannel,
                            record::WorkerRecord,
                            scenarios::Set{ScenarioID},
                            scen_tree::ScenarioTree,
                            create_subproblem::Function,
                            create_subproblem_args::Tuple,
                            create_subproblem_kwargs::NamedTuple,
                            subproblem_callbacks::Vector{SubproblemCallback}
                            )::Nothing
    for scen in scenarios

        # Create subproblem
        sub = create_subproblem(scen,
                                create_subproblem_args...;
                                create_subproblem_kwargs...)

        # Assign variable ids (scenario, stage, index) to variable names and send
        # that map to the master process -- master process expects names to match
        # for variables that have anticipativity constraints
        var_map = report_variable_info(sub, scen_tree)
        put!(output, VariableMap(scen, var_map))

        # Break variables into branch and leaf
        (branch_ids, leaf_ids) = _split_variables(scen_tree,
                                                  collect(keys(var_map)))

        # Save subproblem and relevant data
        record.subproblems[scen] = SubproblemRecord(sub,
                                                    branch_ids,
                                                    leaf_ids,
                                                    subproblem_callbacks
                                                    )
    end

    return

end

function _execute_subproblem_callbacks(sub::SubproblemRecord,
                                       niter::Int,
                                       scen::ScenarioID
                                       )::Nothing

    for spcb in sub.subproblem_callbacks
        spcb.h(spcb.ext, sub.problem, niter, scen)
    end

    return

end

function _initial_solve(output::RemoteChannel,
                        record::WorkerRecord
                        )::Nothing

    for (scen, sub) in record.subproblems

        start = time()
        sts = solve_subproblem(sub.problem)
        stop = time()

        if record.warm_start == true
            warm_start(sub.problem)
        end

        var_vals = report_values(sub.problem, sub.branch_vars)
        put!(output, ReportBranch(scen,
                                  sts,
                                  objective_value(sub.problem),
                                  stop - start,
                                  var_vals)
             )

    end

    return

end

function _penalty_parameter_preprocess(output::RemoteChannel,
                                       record::WorkerRecord,
                                       r::Type{P},
                                       ) where P <: AbstractPenaltyParameter

    if is_subproblem_dependent(r)
        for (scen, sub) in record.subproblems
            pen_map = report_penalty_info(sub.problem, sub.branch_vars, r)
            put!(output, PenaltyInfo(scen, pen_map))
        end
    end

    return

end

function _split_variables(scen_tree::ScenarioTree,
                          vars::Vector{VariableID},
                          )::NTuple{2,Vector{VariableID}}

    branch_vars = Vector{VariableID}()
    leaf_vars = Vector{VariableID}()

    for vid in vars
        if is_leaf(scen_tree, vid)
            push!(leaf_vars, vid)
        else
            push!(branch_vars, vid)
        end
    end

    return (branch_vars, leaf_vars)

end


text:
import Base: read

using GZip
using LightXML

"`Format` is used for write special `parse` (and `read`) methods on it."
abstract Format

"""Returns the `filename`.
Throws an `ErrorException` if the file doesn't exist, or a warning if the file is empty."""
function check_file(filename)
  if !isfile(filename)
    throw(ErrorException(string(filename, " doesn't exist!")))
  elseif filesize(filename) == 0
    warn(string(filename, " is empty!"))
  end
  filename
end

"Returns `true` if the file exists and isn't empty."
isnotemptyfile(filename) = isfile(filename) && filesize(filename) > 0

function _read(completename, filename, format, args...; kargs...) # for using with download, since filename doesn't have file extension
  check_file(filename)
  if endswith(completename, ".xml.gz") || endswith(completename, ".xml")
    document = parse_file(filename)
    try
      parse(document, format, args...; kargs...)
    finally
      free(document)
    end
  else
    io = endswith(completename, ".gz") ? gzopen(filename, "r") : open(filename, "r")
    try
      parse(io, format, args...; kargs...)
    finally
      close(io)
    end
  end
end

"""
```read(pathname, Format [, Type [, … ] ] ) -> Type```

This function opens a file in the `pathname` and calls `parse(io, ...)`for the given `Format` and `Type` on it.
If the  `pathname` is an HTTP or FTP URL, the file is downloaded with `download` in a temporal file.
Gzipped files should end on `.gz`.
"""
function read{T<:Format}(completename::AbstractString, format::Type{T}, args...; kargs...)
  if startswith(completename, "http://") || startswith(completename, "https://") || startswith(completename, "ftp://")
    filename = download(completename)
    try
      _read(completename, filename, format, args...; kargs...)
    finally
      rm(filename)
    end
  else
    _read(completename, completename, format, args...; kargs...) # completename and filename are the same
  end
end


text:
mutable struct Triangle{T} <: Intersectable{T}
    o::SVector{3,T}
    e₁::SVector{3,T}
    e₂::SVector{3,T}
    na::SVector{3,T}
    nb::SVector{3,T}
    nc::SVector{3,T}
end

Triangle(ABC::NTuple{3,SVector{3,T}}, normals::NTuple{3,SVector{3,T}}) where T=
    Triangle(ABC[1], ABC[2]-ABC[1], ABC[3]-ABC[1], normals...)

Triangle(ABC::NTuple{3,SVector{3,T}},
         normal::SVector{3,T}=cross(ABC[2]-ABC[1], ABC[3]-ABC[1])) where T =
             Triangle(ABC, (normal, normal, normal))

function extents(tri::Triangle)
    a,b,c=vertices(tri)
    [extrema([a[i],b[i],c[i]]) for i=1:3]
end

function add_tris!(fun::Function, acc::A, faces::Vector{NTuple{3}}) where {A<:Accelerator}
    for face in faces
        add!(acc, fun(face...))
    end
end

function add_tris!(acc::A, vertices::Vector{SVector{3,T}},
                   faces::Vector{NTuple{3}}) where {A<:Accelerator, T}
    add_tris!(acc, faces) do (i,j,k)
        Triangle((vertices[i], vertices[j], vertices[k]))
    end
end

function add_tris!(acc::A, vertices::Vector{SVector{3,T}}, normals::Vector{SVector{3,T}},
                   faces::Vector{NTuple{3}}) where {A<:Accelerator, T}
    add_tris!(acc, faces) do (i,j,k)
        Triangle((vertices[i], vertices[j], vertices[k]),
                 (normals[i], normals[j], normals[k]))
    end
end

#=
Tomas Möller & Ben Trumbore (1997). Fast, minimum
storage ray-triangle intersection. Journal of Graphics Tools, 2(1),
21–28. http://dx.doi.org/10.1080/10867651.1997.10487468
=#

function Base.intersect(ray::Ray{T}, tri::Triangle{T}) where T
    p = cross(ray.dir, tri.e₂)
    det = dot(tri.e₁, p)

    (det < eltype(tri.o)(1e-6)) && return false

    tv = ray.pos - tri.o

    u = dot(tv, p)

    (u < 0 || u > det) && return false

    v = dot(ray.dir, cross(tv, tri.e₁))
    (v < 0 || u+v > det) && return false

    true
end

function Base.intersect!(i::Intersection{T}, tri::Triangle{T}) where T
    p = cross(i.ray.dir, tri.e₂)
    det = dot(tri.e₁, p)

    (det < eltype(tri.o)(1e-6)) && return

    tv = i.ray.pos - tri.o

    u = dot(tv, p)

    (u < 0 || u > det) && return

    q = cross(tv, tri.e₁)
    v = dot(i.ray.dir, q)
    (v < 0 || u+v > det) && return

    inv_det = inv(det)
    t = dot(tri.e₂, q) * inv_det
    (t > i.t) && return

    u *= inv_det
    v *= inv_det

    i.o = tri
    i.t = t
    i.u = u
    i.v = v
end

normal(tri::Triangle{T}, p::SVector{3,T}, i::Intersection{T}) where T =
    normalize((1 - i.u - i.v)*tri.na + i.u * tri.nb + i.v * tri.nc)

vertices(tri::Triangle{T}) where T= (tri.o, tri.o + tri.e₁, tri.o + tri.e₂)

function aabb(t::Triangle{T}) where T
    a,b,c = vertices(t)
    AABB(compmin(compmin(a,b),c), compmax(compmax(a,b),c))
end

export Triangle, extents, add_tris!, normal, aabb