text:
using BVHFiles
using Test

@testset "graph" begin
    println("Starting test of 'graph'")

    g = BVHGraph(3)
    @test add_edge!(g, 1, 2) == true
    @inferred add_edge!(g, 2, 3, offset = [1.0, 3.0, 5.0])
    @inferred add_vertex!(g, name = "JOINT L3", sequence = :ZXY)

    @inferred has_vertex(g, 2)
    @test has_vertex(g, 2) == true
    @test has_vertex(g, 5) == false

    @inferred has_edge(g, 2, 3)
    @test has_edge(g, 2, 3) == true
    @test has_edge(g, 3, 2) == false

    @inferred props(g, 3)

    @inferred AbstractString name(g)
    @inferred AbstractString name(g, 3)
    @inferred AbstractString name!(g, "Test.bvh")
    @inferred AbstractString name!(g, 1, "JOINT hip")
    @inferred AbstractString name!(g, 2, "JOINT L5")

    @inferred sequence(g)
    @inferred sequence(g, 3)
    sequence!(g, :ZXY)
    sequence!(g, 1, :ZXY)

    @inferred frames(g)
    @test frames(g) == 1

    @inferred frametime(g)
    @test frametime(g) == 0.01

    @inferred Matrix{Float64} rotations(g, 1)
    @inferred Matrix{Float64} rotations!(g, 1, rand(Float64, 10, 3))
    @inferred Matrix{Float64} rotations(g, 1, 5)

    @inferred Matrix{Float64} positions(g)
    @inferred Matrix{Float64} positions(g, 1)
    @inferred Matrix{Float64} positions!(g, rand(Float64, 10, 3))
    @inferred Matrix{Float64} positions!(g, 1, rand(Float64, 10, 3))
    @inferred Matrix{Float64} positions(g, 1, 5)

    @inferred Vector{Float64} offset(g)
    @inferred Vector{Float64} offset(g, 1, 2)
    @inferred Vector{Float64} offset!(g, rand(Float64, 3))
    @inferred Vector{Float64} offset!(g, 1, 2, rand(Float64, 3))

    @inferred find(g, "JOINT L5")
    @test find(g, "JOINT hip") == 1
    @test find(g, "JOINT L5") == 2

    @inferred eltype(g)
    @inferred edgetype(g)
    @test is_directed(g) == true
    @inferred print(g)

    @inferred nv(g)
    @test nv(g) == 4
    @inferred ne(g)
    @test ne(g) == 2

    @inferred vertices(g)
    @inferred edges(g)

    @inferred inneighbors(g, 3)
    @test inneighbors(g, 3) == [2]
    @inferred outneighbors(g, 1)
    @test outneighbors(g, 1) == [2]

    @inferred rem_edge!(g, 2, 3)
    @inferred rem_vertex!(g, 3)
    @inferred fadj(g)
    @inferred badj(g)

    @inferred BVHGraph(3)
end

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

abstract type BATDataVector{T} <: DenseVector{T} end

function Base.merge!(X::BATDataVector, Xs::BATDataVector...)
    for Y in Xs
        append!(X, Y)
    end
    X
end

Base.merge(X::BATDataVector, Xs::BATDataVector...) = merge!(deepcopy(X), Xs...)


text:
module day2

readInput() = readlines(joinpath(@__DIR__, "input.txt"))

function checkerP1(s::String)
    s_lb, s_ub, s_char, pwd = split(s, ['-', ':', ' '], keepempty=false)
    lb, ub, char = parse(Int, s_lb), parse(Int, s_ub), s_char[1]
    return lb <= count(==(char), pwd) <= ub
end

part1(input) = count(checkerP1, input)

function checkerP2(s::String)
    s_lb, s_ub, s_char, pwd = split(s, ['-', ':', ' '], keepempty=false)
    lb, ub, char = parse(Int, s_lb), parse(Int, s_ub), s_char[1]
    return (pwd[lb] == char) ⊻ (pwd[ub] == char)
end

part2(input) = count(checkerP2, input)

end


text:
function update_constant!(s::MOI.LessThan{T}, val::T) where {T}
    return MOI.LessThan{T}(s.upper - val)
end

function update_constant!(s::MOI.GreaterThan{T}, val::T) where {T}
    return MOI.GreaterThan{T}(s.lower - val)
end

function update_constant!(s::MOI.EqualTo{T}, val::T) where {T}
    return MOI.EqualTo{T}(s.value - val)
end

function update_constant!(s::MOI.Interval{T}, val::T) where {T}
    return MOI.Interval{T}(s.lower - val, s.upper - val)
end

# Affine
function update_parameter_in_affine_constraints!(model::Optimizer)
    for (F, S) in keys(model.affine_constraint_cache.dict)
        affine_constraint_cache_inner = model.affine_constraint_cache[F, S]
        affine_added_cache_inner = model.affine_added_cache[F, S]
        if !isempty(affine_constraint_cache_inner)
            update_parameter_in_affine_constraints!(
                model.optimizer,
                model.parameters,
                model.updated_parameters,
                affine_constraint_cache_inner,
                affine_added_cache_inner,
            )
        end
    end
    return model
end

function update_parameter_in_affine_constraints!(
    optimizer::OT,
    parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
    updated_parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
    affine_constraint_cache_inner::MOI.Utilities.DoubleDicts.DoubleDictInner{
        F,
        S,
        V1,
    },
    affine_added_cache_inner::MOI.Utilities.DoubleDicts.DoubleDictInner{F,S,V2},
) where {OT,T,F,S,V1,V2}
    for (ci, param_array) in affine_constraint_cache_inner
        update_parameter_in_affine_constraints!(
            optimizer,
            affine_added_cache_inner[ci],
            param_array,
            parameters,
            updated_parameters,
        )
    end
    return optimizer
end

function update_parameter_in_affine_constraints!(
    optimizer::OT,
    ci::CI,
    param_array::Vector{MOI.ScalarAffineTerm{T}},
    parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
    updated_parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
) where {OT,T,CI}
    param_constant = zero(T)
    for term in param_array
        if !isnan(updated_parameters[p_idx(term.variable)])
            param_constant +=
                term.coefficient * (
                    updated_parameters[p_idx(term.variable)] -
                    parameters[p_idx(term.variable)]
                )
        end
    end
    if param_constant != zero(Float64)
        set = MOI.get(optimizer, MOI.ConstraintSet(), ci)
        set = update_constant!(set, param_constant)
        MOI.set(optimizer, MOI.ConstraintSet(), ci, set)
    end
    return ci
end

function update_parameter_in_affine_objective!(model::Optimizer)
    if !isempty(model.affine_objective_cache)
        objective_constant = zero(Float64)
        for term in model.affine_objective_cache
            if !isnan(model.updated_parameters[p_idx(term.variable)])
                param_old = model.parameters[p_idx(term.variable)]
                param_new = model.updated_parameters[p_idx(term.variable)]
                aux = param_new - param_old
                objective_constant += term.coefficient * aux
            end
        end
        if !iszero(objective_constant)
            F = MOI.get(model.optimizer, MOI.ObjectiveFunctionType())
            f = MOI.get(model.optimizer, MOI.ObjectiveFunction{F}())
            MOI.modify(
                model.optimizer,
                MOI.ObjectiveFunction{F}(),
                MOI.ScalarConstantChange(f.constant + objective_constant),
            )
        end
    end
    return model
end

function update_parameter_in_quadratic_constraints_pc!(model::Optimizer)
    for (ci, fparam) in model.quadratic_constraint_cache_pc
        param_constant = zero(Float64)
        for term in fparam
            if !isnan(model.updated_parameters[p_idx(term.variable)])
                param_constant +=
                    term.coefficient * (
                        model.updated_parameters[p_idx(term.variable)] -
                        model.parameters[p_idx(term.variable)]
                    )
            end
        end
        if param_constant != zero(Float64)
            set = MOI.get(
                model.optimizer,
                MOI.ConstraintSet(),
                model.quadratic_added_cache[ci],
            )
            set = update_constant!(set, param_constant)
            MOI.set(
                model.optimizer,
                MOI.ConstraintSet(),
                model.quadratic_added_cache[ci],
                set,
            )
        end
    end
end

function update_parameter_in_quadratic_objective_pc!(model::Optimizer)
    if !isempty(model.quadratic_objective_cache_pc)
        objective_constant = zero(Float64)
        for term in model.quadratic_objective_cache_pc
            if !isnan(model.updated_parameters[p_idx(term.variable)])
                param_old = model.parameters[p_idx(term.variable)]
                param_new = model.updated_parameters[p_idx(term.variable)]
                aux = param_new - param_old
                objective_constant += term.coefficient * aux
            end
        end
        if !iszero(objective_constant)
            F = MOI.get(model.optimizer, MOI.ObjectiveFunctionType())
            f = MOI.get(model.optimizer, MOI.ObjectiveFunction{F}())
            MOI.modify(
                model.optimizer,
                MOI.ObjectiveFunction{F}(),
                MOI.ScalarConstantChange(f.constant + objective_constant),
            )
        end
    end
    return model
end

function update_parameter_in_quadratic_constraints_pp!(model::Optimizer)
    for (ci, fparam) in model.quadratic_constraint_cache_pp
        param_constant = zero(Float64)
        for term in fparam
            if !isnan(model.updated_parameters[p_idx(term.variable_1)]) &&
               !isnan(model.updated_parameters[p_idx(term.variable_2)])
                param_constant +=
                    term.coefficient * (
                        (
                            model.updated_parameters[p_idx(term.variable_1)] *
                            model.updated_parameters[p_idx(term.variable_2)]
                        ) - (
                            model.parameters[p_idx(term.variable_1)] *
                            model.parameters[p_idx(term.variable_2)]
                        )
                    )
            elseif !isnan(model.updated_parameters[p_idx(term.variable_1)])
                param_constant +=
                    term.coefficient *
                    model.parameters[p_idx(term.variable_2)] *
                    (
                        model.updated_parameters[p_idx(term.variable_1)] -
                        model.parameters[p_idx(term.variable_1)]
                    )
            elseif !isnan(model.updated_parameters[p_idx(term.variable_2)])
                param_constant +=
                    term.coefficient *
                    model.parameters[p_idx(term.variable_1)] *
                    (
                        model.updated_parameters[p_idx(term.variable_2)] -
                        model.parameters[p_idx(term.variable_2)]
                    )
            end
        end
        if param_constant != zero(Float64)
            set = MOI.get(
                model.optimizer,
                MOI.ConstraintSet(),
                model.quadratic_added_cache[ci],
            )
            set = update_constant!(set, param_constant)
            MOI.set(
                model.optimizer,
                MOI.ConstraintSet(),
                model.quadratic_added_cache[ci],
                set,
            )
        end
    end
    return model
end

function update_parameter_in_quadratic_objective_pp!(model::Optimizer)
    if !isempty(model.quadratic_objective_cache_pp)
        objective_constant = zero(Float64)
        for term in model.quadratic_objective_cache_pp
            if !isnan(model.updated_parameters[p_idx(term.variable_1)]) &&
               !isnan(model.updated_parameters[p_idx(term.variable_2)])
                objective_constant +=
                    term.coefficient * (
                        (
                            model.updated_parameters[p_idx(term.variable_1)] *
                            model.updated_parameters[p_idx(term.variable_2)]
                        ) - (
                            model.parameters[p_idx(term.variable_1)] *
                            model.parameters[p_idx(term.variable_2)]
                        )
                    )
            elseif !isnan(model.updated_parameters[p_idx(term.variable_1)])
                objective_constant +=
                    term.coefficient *
                    model.parameters[p_idx(term.variable_2)] *
                    (
                        model.updated_parameters[p_idx(term.variable_1)] -
                        model.parameters[p_idx(term.variable_1)]
                    )
            elseif !isnan(model.updated_parameters[p_idx(term.variable_2)])
                objective_constant +=
                    term.coefficient *
                    model.parameters[p_idx(term.variable_1)] *
                    (
                        model.updated_parameters[p_idx(term.variable_2)] -
                        model.parameters[p_idx(term.variable_2)]
                    )
            end
        end
        if objective_constant != zero(Float64)
            F = MOI.get(model.optimizer, MOI.ObjectiveFunctionType())
            f = MOI.get(model.optimizer, MOI.ObjectiveFunction{F}())
            MOI.modify(
                model.optimizer,
                MOI.ObjectiveFunction{F}(),
                MOI.ScalarConstantChange(f.constant + objective_constant),
            )
        end
    end
    return model
end

function update_parameter_in_quadratic_constraints_pv!(model::Optimizer)
    for (ci, quad_aff_vars) in model.quadratic_constraint_cache_pv
        # We need this dictionary because we could have terms like
        # p_1 * v_1 + p_2 * v_1 and we should add p_1 and p_2 on the update
        new_coeff_per_variable = Dict{MOI.VariableIndex,Float64}()
        for term in quad_aff_vars
            # Here we use the convention that the parameter always comes as the first variables
            # in the caches
            if !isnan(model.updated_parameters[p_idx(term.variable_1)])
                coef = term.coefficient
                param_new = model.updated_parameters[p_idx(term.variable_1)]
                if haskey(new_coeff_per_variable, term.variable_2)
                    new_coeff_per_variable[term.variable_2] += param_new * coef
                else
                    new_coeff_per_variable[term.variable_2] = param_new * coef
                end
            end
        end
        if haskey(
            model.quadratic_constraint_variables_associated_to_parameters_cache,
            ci,
        )
            for aff_term in
                model.quadratic_constraint_variables_associated_to_parameters_cache[ci]
                coef = aff_term.coefficient
                if haskey(new_coeff_per_variable, aff_term.variable)
                    new_coeff_per_variable[aff_term.variable] += coef
                else
                    new_coeff_per_variable[aff_term.variable] = coef
                end
            end
        end
        old_ci = model.quadratic_added_cache[ci]
        changes = Vector{MOI.ScalarCoefficientChange}(
            undef,
            length(new_coeff_per_variable),
        )
        i = 1
        for (vi, value) in new_coeff_per_variable
            changes[i] = MOI.ScalarCoefficientChange(vi, value)
            i += 1
        end
        # Make multiple changes at once.
        MOI.modify(
            model.optimizer,
            fill(old_ci, length(new_coeff_per_variable)),
            changes,
        )
    end
    return model
end

function update_parameter_in_quadratic_objective_pv!(model::Optimizer)
    if !isempty(model.quadratic_objective_cache_pv)
        # We need this dictionary because we could have terms like
        # p_1 * v_1 + p_2 * v_1 and we should add p_1 and p_2 on the update
        new_coeff_per_variable = Dict{MOI.VariableIndex,Float64}()
        for term in model.quadratic_objective_cache_pv
            # Here we use the convention that the parameter always comes as the first variables
            # in the caches
            if !isnan(model.updated_parameters[p_idx(term.variable_1)])
                coef = term.coefficient
                param_new = model.updated_parameters[p_idx(term.variable_1)]
                if haskey(new_coeff_per_variable, term.variable_2)
                    new_coeff_per_variable[term.variable_2] += param_new * coef
                else
                    new_coeff_per_variable[term.variable_2] = param_new * coef
                end
            end
        end

        for aff_term in
            model.quadratic_objective_variables_associated_to_parameters_cache
            coef = aff_term.coefficient
            if haskey(new_coeff_per_variable, aff_term.variable)
                new_coeff_per_variable[aff_term.variable] += coef
            else
                new_coeff_per_variable[aff_term.variable] = coef
            end
        end

        F_pv = MOI.get(model.optimizer, MOI.ObjectiveFunctionType())
        changes = Vector{MOI.ScalarCoefficientChange}(
            undef,
            length(new_coeff_per_variable),
        )
        i = 1
        for (vi, value) in new_coeff_per_variable
            changes[i] = MOI.ScalarCoefficientChange(vi, value)
            i += 1
        end
        MOI.modify(model.optimizer, MOI.ObjectiveFunction{F_pv}(), changes)
    end
    return model
end

# Vector Affine
function update_parameter_in_vector_affine_constraints!(model::Optimizer)
    for (F, S) in keys(model.vector_constraint_cache.dict)
        vector_constraint_cache_inner = model.vector_constraint_cache[F, S]
        if !isempty(vector_constraint_cache_inner)
            update_parameter_in_vector_affine_constraints!(
                model.optimizer,
                model.parameters,
                model.updated_parameters,
                vector_constraint_cache_inner,
            )
        end
    end
    return model
end

function update_parameter_in_vector_affine_constraints!(
    optimizer::OT,
    parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
    updated_parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
    vector_constraint_cache_inner::MOI.Utilities.DoubleDicts.DoubleDictInner{
        F,
        S,
        V,
    },
) where {OT,T,F,S,V}
    for (ci, param_array) in vector_constraint_cache_inner
        update_parameter_in_vector_affine_constraints!(
            optimizer,
            ci,
            param_array,
            parameters,
            updated_parameters,
        )
    end

    return optimizer
end

function update_parameter_in_vector_affine_constraints!(
    optimizer::OT,
    ci::CI,
    param_array::Vector{MOI.VectorAffineTerm{T}},
    parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
    updated_parameters::MOI.Utilities.CleverDicts.CleverDict{
        ParameterIndex,
        T,
        typeof(MOI.Utilities.CleverDicts.key_to_index),
        typeof(MOI.Utilities.CleverDicts.index_to_key),
    },
) where {OT,T,CI}
    cf = MOI.get(optimizer, MOI.ConstraintFunction(), ci)

    n_dims = length(cf.constants)
    param_constants = zeros(T, n_dims)

    for term in param_array
        vi = term.scalar_term.variable

        if !isnan(updated_parameters[p_idx(vi)])
            param_constants[term.output_index] =
                term.scalar_term.coefficient *
                (updated_parameters[p_idx(vi)] - parameters[p_idx(vi)])
        end
    end

    if param_constants != zeros(T, n_dims)
        MOI.modify(
            optimizer,
            ci,
            MOI.VectorConstantChange(cf.constants + param_constants),
        )
    end

    return ci
end

function update_parameters!(model::Optimizer)
    update_parameter_in_affine_constraints!(model)
    update_parameter_in_affine_objective!(model)
    update_parameter_in_quadratic_constraints_pc!(model)
    update_parameter_in_quadratic_objective_pc!(model)
    update_parameter_in_quadratic_constraints_pv!(model)
    update_parameter_in_quadratic_objective_pv!(model)
    update_parameter_in_quadratic_constraints_pp!(model)
    update_parameter_in_quadratic_objective_pp!(model)
    update_parameter_in_vector_affine_constraints!(model)

    # Update parameters and put NaN to indicate that the parameter has been updated
    for (parameter_index, val) in model.updated_parameters
        if !isnan(val)
            model.parameters[parameter_index] = val
            model.updated_parameters[parameter_index] = NaN
        end
    end

    return model
end


text:
using LinearAlgebra
using TensorOperations
import Strided: UnsafeStridedView

@testset "Arrays" begin

    idx = [ (1,1,1,1), (2,2,2,2), (3,3,3,3)]
    val = rand(3)
    A = FermiSparse(idx, val)
    @test ndims(A) == 4
    @test length(A) == 3
    @test eltype(A) == Float64

end

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Probabilistic <: Forecast
        label::String
        percentiles::Vector{Float64}
        data::TimeSeries.TimeArray
    end

A Probabilistic forecast for a particular data field in a Component.

# Arguments
- `label::String`: label of component parameter forecasted
- `percentiles::Vector{Float64}`: Percentiles for the probabilistic forecast
- `data::TimeSeries.TimeArray`: timestamp - scalingfactor
"""
mutable struct Probabilistic <: Forecast
    "label of component parameter forecasted"
    label::String
    "Percentiles for the probabilistic forecast"
    percentiles::Vector{Float64}
    "timestamp - scalingfactor"
    data::TimeSeries.TimeArray
end

function Probabilistic(; label, percentiles, data, )
    Probabilistic(label, percentiles, data, )
end

"""Get Probabilistic label."""
get_label(value::Probabilistic) = value.label
"""Get Probabilistic percentiles."""
get_percentiles(value::Probabilistic) = value.percentiles
"""Get Probabilistic data."""
get_data(value::Probabilistic) = value.data


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Actors, Printf
import Actors: spawn

const maxitems = 10

mutable struct Store
    capacity::Int
    items::Array{Any,1}
    prod::Array{Link,1}
    cons::Array{Link,1}
    count::Int
    Store(capacity::Int) = new(capacity, Any[], Link[], Link[], 0)
end

struct Prod
    name::String
    store::Link
end

struct Cons
    name::String
    store::Link
end

@msg Put Full Done Ok Take Empty Notify

available(s::Store) = 0 < length(s.items) < s.capacity
isfull(s::Store) = length(s.items) ≥ s.capacity
Base.isempty(s::Store) = isempty(s.items)

function (s::Store)(::Put, prod, item)
    if isfull(s)
        send(prod, Full(), item)
        push!(s.prod, prod)
    elseif s.count < maxitems
        push!(s.items, item)
        s.count += 1
        s.count == maxitems ?
            send(prod, Done()) :
            send(prod, Ok(), item)
        !isempty(s.cons) && send(popfirst!(s.cons), Notify())
    else
        send(prod, Done())
    end
end
function (s::Store)(::Take, cons)
    if isempty(s)
        send(cons, Empty())
        push!(s.cons, cons)
    else
        send(cons, popfirst!(s.items))
        !isempty(s.prod) && send(popfirst!(s.prod), Notify())
    end
end

function prod_start(p::Prod, start)
    become(producing, p)
    send(self(), start+1)
    send(prn, "producer $(p.name) started")
end
function producing(p::Prod, item)
    sleep(rand())
    send(p.store, Put(), self(), item)
end
function producing(p::Prod, ::Ok, item)
    send(prn, "producer $(p.name) delivered item $item")
    send(self(), item+1)
end
function producing(p::Prod, ::Full, item)
    send(prn, "producer $(p.name) stalled with item $item")
    become(stalled, p, item)
end
function producing(p::Prod, ::Done)
    send(prn, "producer $(p.name) done")
    stop()
end
function stalled(p::Prod, item, ::Notify)
    send(p.store, Put(), self(), item)
    become(producing, p)
end

function cons_start(c::Cons)
    become(buying, c)
    send(c.store, Take(), self())
    send(prn, "consumer $(c.name) started")
end
function buying(c::Cons, item)
    become(consuming)
    send(self(), c)
    send(prn, "consumer $(c.name) got item $item")
end
function buying(c::Cons, ::Empty)
    become(waiting, c)
    send(prn, "consumer $(c.name) found store empty")
end
function consuming(c)
    sleep(rand())
    become(buying, c)
    send(c.store, Take(), self())
end
function waiting(c::Cons, ::Notify)
    become(buying, c)
    send(c.store, Take(), self())
end

prn = spawn(s->print(@sprintf("%s\n", s)))
st = spawn(Store(5))
pr1 = spawn(prod_start, Prod("A", st), 100)
pr2 = spawn(prod_start, Prod("B", st), 200)
pr3 = spawn(prod_start, Prod("C", st), 300)
cs1 = spawn(cons_start, Cons("U", st))
cs2 = spawn(cons_start, Cons("V", st))

foreach(x->send(x), (pr1,pr2,pr3,cs1,cs2))


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# Depth-first search
export DFS
"""
  `DFS(fst::WFST, initial; filter=arc->true)`

Creates an iterator that performs a depth-first search (DFS) of the `fst` starting from the state `initial`.
The function filter can be specify in order to perform the search ignoring arcs with user-defined properties.

For each iterator the tuple `(p,s,n,d,e)` is returned, where:

* `p` is the parent parent state (p==0 if there is no parent)
* `s` is the current state
* `n` next state: `n==0` if `s` has no arcs or if completely explored. These will give you standard DFS iterations.
* `d` is a `Bool`: if `true` means a new state is discovered, can be used to pick up DFS its (see example below)
* `e` is a `Bool`, if `true` means all arcs have been inspected
* `a` inspected arc (empty if `d == false`)

The following cases/colors are possible:
* `d == false && n ==0` -> black state, state has been checked completely
* `d == false && n !=0` -> `n` is already gray
* `d == true -> n` state is colored gray, i.e. is visited for the first time

# Example

```julia
julia> fst = WFST(Dict(1=>1));

julia> add_arc!(fst, 1, 2, 1, 1, 1.0);

julia> add_arc!(fst, 1, 5, 1, 1, 1.0);

julia> add_arc!(fst, 1, 3, 1, 1, 1.0);

julia> add_arc!(fst, 2, 6, 1, 1, 1.0);

julia> add_arc!(fst, 2, 4, 1, 1, 1.0);

julia> add_arc!(fst, 5, 4, 1, 1, 1.0);

julia> add_arc!(fst, 3, 4, 1, 1, 1.0);

julia> add_arc!(fst, 3, 7, 1, 1, 1.0);

julia> initial!(fst,1)
WFST #states: 7, #arcs: 8, #isym: 1, #osym: 1
|1/0.0f0|
1:1/1.0f0 → (2)
1:1/1.0f0 → (5)
1:1/1.0f0 → (3)
(2)
1:1/1.0f0 → (6)
1:1/1.0f0 → (4)
(3)
1:1/1.0f0 → (4)
1:1/1.0f0 → (7)
(4)
(5)
1:1/1.0f0 → (4)
(6)
(7)

julia> dfs = FiniteStateTransducers.DFS(fst,1); # DFS starting from state 1

julia> for (p,s,n,d,e,a) in dfs
         println("\$p,\$s,\$n")
         if d == true
           println("visiting first time \$n (Gray)")
         else
           if e
             println("completed state \$s (Black)")
           else
             println("node \$n already visited")
           end
         end
       end
0,1,2
visiting first time 2 (Gray)
1,2,6
visiting first time 6 (Gray)
2,6,0
completed state 6 (Black)
1,2,4
visiting first time 4 (Gray)
2,4,0
completed state 4 (Black)
1,2,0
completed state 2 (Black)
0,1,5
visiting first time 5 (Gray)
1,5,4
node 4 already visited
1,5,0
completed state 5 (Black)
0,1,3
visiting first time 3 (Gray)
1,3,4
node 4 already visited
1,3,7
visiting first time 7 (Gray)
3,7,0
completed state 7 (Black)
1,3,0
completed state 3 (Black)
0,1,0
completed state 1 (Black)

```
"""
struct DFS{W, A, D, F <: WFST{W,A,D}}
  fst::F
  is_visited::Vector{Bool}
  is_completed::Vector{Bool} #TODO possibly this can be rm
  initial::D
  filter::Function
  empty_arc::Arc{W,D}
  function DFS(fst::F, initial::D; filter=arc->true) where {W,A,D,F <: WFST{W,A,D}}
    is_visited = zeros(Bool, length(fst))
    is_completed = zeros(Bool, length(fst))
    empty_arc = Arc{W,D}(D(0),D(0),W(0),D(0))
    new{W,A,D,F}(fst, is_visited, is_completed, initial, filter, empty_arc)
  end
end

function init(d::DFS{W,A,D,F}) where {W,A,D,F}
  s = d.initial
  q = Stack{Tuple{D,D}}()
  d.is_visited[s] = true
  push!(q,(s,1)) # tuple having (parent state,arc iterator state)
end

function Base.iterate(d::DFS{W,A,D}, q = init(d)) where {W,A,D}
  if isempty(q)
    return nothing
  else
    s, arc_it_state = pop!(q)
    arcs = d.fst[s]
    state_arc = iterate(arcs,arc_it_state)
    while state_arc != nothing && (!d.filter(state_arc[1]))
      state_arc = iterate(arcs,state_arc[2])
    end

    if !isempty(q)
      p, = first(q)
    else
      p = D(0)
    end
    if state_arc == nothing
      d.is_completed[s] = true
      # s is completed (black)
      return (p,s,D(0),false,true,d.empty_arc), q
    else
      arc, a = state_arc
      push!(q,(s,a))
      n = get_nextstate(arc)
      if !(d.is_visited[n])
        d.is_visited[n] = true
        push!(q,(n,1))
        # n is newly visit and represent DFS indices
        # labeled as gray
        return (p,s,n,true ,false,arc), q
      else
        return (p,s,n,false,false,arc), q
      end
    end
  end
end

# simple implementation (just for testing pourposes)
function recursive_dfs(fst, i = get_initial(fst;single=true))
  is_visited = zeros(Bool, length(fst))
  recursive_dfs!(is_visited, fst, i)
  return is_visited
end

function recursive_dfs!(is_visited, fst, i)
  # is_visited must be zeros(Bool,length(fst))
  is_visited[i] = true
  for arc in fst[i]
    w = get_nextstate(arc)
    if !is_visited[w]
      recursive_dfs!(is_visited, fst, w)
    end
  end
end


text:
using HDF5

"""
    get_group(h5::Union{HDF5.File,HDF5.Group}, name) -> HDF5.Group

Returns an `HDF5.Group` object for `name` in `h5`.  If `name` does not exist in
`h5`, a group named `name` will be created.  If `name` exists in `h5`, but it is
not as group, an error will be thrown.
"""
function get_group(h5::Union{HDF5.File,HDF5.Group}, name)::HDF5.Group
    haskey(h5, name) ? h5[name] : create_group(h5, name)
end

"""
    write_dataset(h5::Union{HDF5.File,HDF5.Group}, name, data; opts=())

Write `data` to dataset `name` in `h5`.  If `name` does not already exist it
will be created with any kwargs that are given in `opts`
"""
function write_dataset(h5::Union{HDF5.File,HDF5.Group}, name::AbstractString, data;
                     opts=()
                    )
    created = false
    if !haskey(h5, name)
        create_dataset(h5, name, data; opts...)
        created = true
    end
    write(h5[name], data)

    created
end

"""
    isstoreworthy(x)::Bool

Determines whether `x` is worthy of being strored in a BeamformerRecipe file.
"""
function isstoreworthy(x)
    @warn "unsupprted type $(typeof(x)) is not storeworthy" _module=nothing _file=nothing
    false
end
isstoreworthy(x::Nothing) = false
isstoreworthy(x::AbstractString) = length(x) > 0
isstoreworthy(x::AbstractFloat) = !isnan(x)
# The only integer scalar fields in BeamformerRecipe files are positive integers
isstoreworthy(x::Integer) = x > 0
isstoreworthy(x::AbstractArray) = sizeof(x) > 0

"""
    groupname(info::Info) -> String

Name of HDF5 group that holds datasets for `typeof(info)` structs.
"""
groupname(info::Info) = "/$(lowercase(string(typeof(info))))"

"""
    to_hdf5(group::HDF5.Group, info::Info; allopts...)

Store the data in `info` in the HDF5 datasets within `group`.  When splatted,
`allopts` should be pairs mapping symbols of `info` field names to NamedTuples
of HDF5 dataset creation options.
"""
function to_hdf5(group::HDF5.Group, info::Info; allopts...)
    if HDF5.name(group) != groupname(info)
        @warn "adding $typeof(info) to group $groupname(info)"
    end
    for name in propertynames(info)
        data = getproperty(info, name)
        isstoreworthy(data) || continue
        write_dataset(group, string(name), data; opts=get(allopts, name, ()))
    end
end

function to_hdf5(h5::HDF5.File, info::Info)
    group = get_group(h5, groupname(info))
    to_hdf5(group, info)
end

function to_hdf5(::HDF5.File, ::Nothing)
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sin(x)^5, x) == :(-1//3*c^-1*cos(x)^3+(1/2)*b*c^-2*cos(x)^2+-1*c^-3*(b^2+-1*c*(a+2c))*cos(x)+c^-4*(b^2+-4*a*c)^-1//2*(b^4+2*c^2*(a+c)^2+-2*c*b^2*(c+2a))*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*cos(x)))+(1/2)*b*c^-4*(b^2+-2*c*(a+c))*log(a+b*cos(x)+c*cos(x)^2))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sin(x)^3, x) == :(c^-1*cos(x)+-1//2*b*c^-2*log(a+b*cos(x)+c*cos(x)^2)+-1*c^-2*(b^2+-4*a*c)^-1//2*(b^2+-2*c*(a+c))*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*cos(x))))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sin(x), x) == :(2*(b^2+-4*a*c)^-1//2*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*cos(x))))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*csc(x), x) == :((2a+2b+2c)^-1*log(1+-1*cos(x))+-1*(-2b+2a+2c)^-1*log(1+cos(x))+b*(a+b+c)^-1*(-2b+2a+2c)^-1*log(a+b*cos(x)+c*cos(x)^2)+-1*(b^2+-4*a*c)^-1//2*(a+b+c)^-1*(a+c+-1b)^-1*(b^2+-2*c^2+-2*a*c)*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*cos(x))))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*csc(x)^3, x) == :(-1//4*(a+c+-1b)^-2*(a+-2b+3c)*log(1+cos(x))+1//4*(a+b+c)^-2*(a+2b+3c)*log(1+-1*cos(x))+(b^2+-4*a*c)^-1//2*(a^2+c^2+-1*b^2+2*a*c)^-2*(b^4+2*c^2*(a+c)^2+-2*c*b^2*(c+2a))*arctanh((b^2+-4*a*c)^-1//2*(b+2*c*cos(x)))+(a+b+c)^-1*(-2b+2a+2c)^-1*csc(x)^2*(b+-1*(a+c)*cos(x))+-1//2*b*(a^2+c^2+-1*b^2+2*a*c)^-2*(b^2+-2*c*(a+c))*log(a+b*cos(x)+c*cos(x)^2))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sin(x)^2, x) == :(-1*x*c^-1+c^-1*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2b+2*(b^2+-4*a*c)^-1//2*(b^2+-2*c*(a+c)))*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+c^-1*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2b+-2*(b^2+-4*a*c)^-1//2*(b^2+-2*c*(a+c)))*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*csc(x)^2, x) == :((1+cos(x))^-1*(-2b+2a+2c)^-1*sin(x)+-1*(1+-1*cos(x))^-1*(2a+2b+2c)^-1*sin(x)+-2*b*c*(a+b+c)^-1*(a+c+-1b)^-1*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(1+-1*b^-1*(b^2+-4*a*c)^-1//2*(b^2+-2*c*(a+c)))*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+-2*b*c*(a+b+c)^-1*(a+c+-1b)^-1*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(1+b^-1*(b^2+-4*a*c)^-1//2*(b^2+-2*c*(a+c)))*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((-2+cos(x)^2+cos(x))^-1*sin(x), x) == :(-1//3*log(1+-1*cos(x))+1//3*log(2+cos(x)))
@test integrate((4+cos(x)^2+-5*cos(x))^-1*sin(x), x) == :(-1//3*log(4+-1*cos(x))+1//3*log(1+-1*cos(x)))
@test integrate((3+cos(x)^2+-2*cos(x))^-1*sin(x), x) == :((1/2)*2^(1/2)*arctan((1/2)*2^(1/2)*(1+-1*cos(x))))
@test integrate((13+cos(x)^2+-4*cos(x))^-2*sin(x), x) == :(-1//54*arctan(-2//3+1//3*cos(x))+(234+-72*cos(x)+18*cos(x)^2)^-1*(2+-1*cos(x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*cos(x)^4, x) == :((1/2)*x*c^-1+x*c^-3*(b^2+-1*a*c)+(1/2)*c^-1*cos(x)*sin(x)+-1*b*c^-2*sin(x)+-1*c^-3*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2*b^3+2*(b^2+-4*a*c)^-1//2*(b^4+2*a^2*c^2+-4*a*c*b^2)+-4*a*b*c)*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+-1*c^-3*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2*b^3+-2*(b^2+-4*a*c)^-1//2*(b^4+2*a^2*c^2+-4*a*c*b^2)+-4*a*b*c)*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*cos(x)^3, x) == :(c^-1*sin(x)+-1*b*x*c^-2+c^-2*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2*b^2+-2*a*c+2*b^3*(b^2+-4*a*c)^-1//2+-6*a*b*c*(b^2+-4*a*c)^-1//2)*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+c^-2*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2*b^2+-2*a*c+-2*b^3*(b^2+-4*a*c)^-1//2+6*a*b*c*(b^2+-4*a*c)^-1//2)*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*cos(x)^2, x) == :(x*c^-1+-1*c^-1*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2b+2*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c))*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+-1*c^-1*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2b+-2*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c))*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*cos(x), x) == :((b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2+2*b*(b^2+-4*a*c)^-1//2)*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(2+-2*b*(b^2+-4*a*c)^-1//2)*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1, x) == :(-4*c*(b^2+-4*a*c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+4*c*(b^2+-4*a*c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sec(x), x) == :(a^-1*arctanh(sin(x))+-2*c*a^-1*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(1+-1*b*(b^2+-4*a*c)^-1//2)*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+-2*c*a^-1*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(1+b*(b^2+-4*a*c)^-1//2)*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sec(x)^2, x) == :(a^-1*tan(x)+-1*b*a^-2*arctanh(sin(x))+2*b*c*a^-2*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(1+-1*b^-1*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c))*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+2*b*c*a^-2*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(1+b^-1*(b^2+-4*a*c)^-1//2*(b^2+-2*a*c))*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))
@test integrate((a+b*cos(x)+c*cos(x)^2)^-1*sec(x)^3, x) == :((1/2)*a^-1*arctanh(sin(x))+a^-3*(b^2+-1*a*c)*arctanh(sin(x))+(1/2)*a^-1*sec(x)*tan(x)+-1*b*a^-2*tan(x)+-2*c*a^-3*(b^2+-4*a*c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*(b^3+(b^2+-4*a*c)^(1/2)*(b^2+-1*a*c)+-3*a*b*c)*arctan((b+-1*(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+-1*(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x))+2*c*a^-3*(b^2+-4*a*c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+-2c)^-1//2*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*(b^3+-1*(b^2+-4*a*c)^(1/2)*(b^2+-1*a*c)+-3*a*b*c)*arctan((b+(b^2+-4*a*c)^(1/2)+-2c)^(1/2)*(b+(b^2+-4*a*c)^(1/2)+2c)^-1//2*tan((1/2)*x)))

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.9 trig^m (a+b cos^n+c cos^(2 n))^p.jl


text:
using GeneticsMakie
using Test
using CairoMakie
using CSV
using DataFrames
using SnpArrays
using Statistics

@testset "Plotting genes/isoforms" begin
    gencode = CSV.read("data/gencode.gtf", DataFrame)
    GeneticsMakie.parsegtf!(gencode)
    @test ncol(gencode) == 14

    gene = "KMT2E"
    chr, start, stop = GeneticsMakie.findgene(gene, gencode)
    @test ("7", 104581390, 104755466) == (chr, start, stop)
    range1 = start - 1e6
    range2 = stop + 1e6

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, chr, range1, range2, gencode)
    GeneticsMakie.labelgenome(f[1, 1, Bottom()], chr, range1, range2)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, "6", range1, range2, gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, chr, start, gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, gene, gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, chr, range1, range2, (["KMT2E", "SETD1A"], ["red", "blue"]), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, chr, start, (["KMT2E", "SETD1A"], ["red", "blue"]), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, gene, (["KMT2E", "SETD1A"], ["red", "blue"]), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, chr, range1, range2, ("KMT2E", "red"), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, "6", range1, range2, ("KMT2E", "red"), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, chr, start, ("KMT2E", "red"), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgenes!(ax, gene, ("KMT2E", "red"), gencode)
    save("gene.png", f)
    @test isfile("gene.png")
    rm("gene.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotisoforms!(ax, gene, gencode; orderby = ["ENST00000667857", "ENST00000482560", "Random"])
    save("isoform.png", f)
    @test isfile("isoform.png")
    rm("isoform.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotisoforms!(ax, gene, gencode; text = :l)
    save("isoform.png", f)
    @test isfile("isoform.png")
    rm("isoform.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotisoforms!(ax, gene, gencode; text = :b)
    save("isoform.png", f)
    @test isfile("isoform.png")
    rm("isoform.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotisoforms!(ax, gene, gencode; text = :r)
    save("isoform.png", f)
    @test isfile("isoform.png")
    rm("isoform.png")
end

@testset "Munging summmary stats" begin
    gwas = CSV.read("data/sumstats.csv", DataFrame)
    GeneticsMakie.mungesumstats!(gwas)
    @test ncol(gwas) == 7

    gwas = CSV.read("data/sumstats.csv", DataFrame)
    select!(gwas, Not("BETA"))
    GeneticsMakie.mungesumstats!(gwas)
    @test ncol(gwas) == 6

    gwas = CSV.read("data/sumstats.csv", DataFrame)
    select!(gwas, "#CHROM", :POS, :PVAL)
    GeneticsMakie.mungesumstats!(gwas)
    @test ncol(gwas) == 4

    df = GeneticsMakie.findgwasloci(gwas)
    @test nrow(df) == 1
    df = GeneticsMakie.findgwasloci([gwas, gwas])
    @test nrow(df) == 1
end

@testset "Plotting GWAS" begin
    gwas = DataFrame(
        CHR = rand(string.(collect(1:22)), 1000),
        BP = rand(1:1000, 1000),
        P = rand(1000)
    )

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotgwas!(ax, gwas)
    save("manhattan.png", f)
    @test isfile("manhattan.png")
    rm("manhattan.png")
end

@testset "Plotting LD" begin
    LD = rand(10, 10)
    LD = LD + transpose(LD)
    for i in 1:size(LD, 1)
        LD[i, i] = 1
    end
    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotld!(ax, LD)
    save("ld.png", f)
    @test isfile("ld.png")
    rm("ld.png")
end

@testset "Plotting LocusZoom" begin
    kgp = SnpData("data/kgp")
    gwas = CSV.read("data/locus.csv", DataFrame)
    gwas.CHR = string.(gwas.CHR)

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotlocus!(ax, "7", 103581390, 105755466, gwas)
    save("locuszoom.png", f)
    @test isfile("locuszoom.png")
    rm("locuszoom.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotlocus!(ax, "7", 105755466, 103581390, gwas)
    save("locuszoom.png", f)
    @test isfile("locuszoom.png")
    rm("locuszoom.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotlocus!(ax, "7", 103581390, 105755466, gwas; ld = kgp)
    save("locuszoom.png", f)
    @test isfile("locuszoom.png")
    rm("locuszoom.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotlocus!(ax, "7", 103581390, 105755466, gwas; ld = (kgp, "rs11764361"))
    save("locuszoom.png", f)
    @test isfile("locuszoom.png")
    rm("locuszoom.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotlocus!(ax, "7", 103581390, 105755466, gwas; ld = (kgp, "rs111931861"))
    save("locuszoom.png", f)
    @test isfile("locuszoom.png")
    rm("locuszoom.png")

    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotlocus!(ax, "7", 103581390, 105755466, gwas; ld = (kgp, "rs11764361"), showsnplabel = false)
    save("locuszoom.png", f)
    @test isfile("locuszoom.png")
    rm("locuszoom.png")
end

@testset "Plotting QQ plot" begin
    P = rand(1000)
    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotqq!(ax, P)
    save("qq.png", f)
    @test isfile("qq.png")
    rm("qq.png")

    df = DataFrame(P = P)
    f = Figure()
    ax = Axis(f[1, 1])
    GeneticsMakie.plotqq!(ax, P)
    save("qq.png", f)
    @test isfile("qq.png")
    rm("qq.png")
end

@testset "Plotting correlation" begin
    f = Figure()
    ax = Axis(f[1, 1])
    n = 10
    GeneticsMakie.plotrg!(ax, (rand(n, n) .* 2 .- 1), string.(1:n), circle = true)
    colsize!(f.layout, 1, Aspect(1, 1))
    rowsize!(f.layout, 1, 18 * n)
    save("cor.png", f)
    @test isfile("cor.png")
    rm("cor.png")
end

@testset "Plotting TWAS" begin
end


text:
"""
    optimize!(CGP, VarDervs, OptParams, elocs, forces_w, γ[, iterative])

Optimize the network parameters in-place with Stochastic Reconfiguration by
MATRIX-FREE inversion of the S-matrix.

# Arguments
- `CGP::CGPARAMS`: A conjugate gradient parameter struct for when Conjugate
Gradient is performed instead of explicit inversion of the S-matrix.
- `VarDervs::VARDERIVS`: contains sampled variational derivatives and
regularisation parameter.
- `OptParams::OPTPARAMSMFREE`: contains optimization variables: forces,
wave function parameter update vector ΔW (NOT the S-matrix).
- `elocs`: Each column stores the sampled local energies from a worker.
- `forces_w`: Stores the result of adjoint(O) * E (both standardized) of each
worker in a column, the summed forces are stored in `OptParams.f`.
- `γ::Float64`: learning parameter
- `NetSettings`: struct containing all parameters of the RBM network

"""
function optimize!(CGP::CGPARAMS, VarDervs::VARDERIVS, OptParams::OPTPARAMSMFREE,
    elocs, forces_w, γ, NetSettings)
    @unpack o_tot = VarDervs
    @unpack f = OptParams
    o_tot .-=  mean(mean(o_tot, dims=2), dims=3) # Mean over all vardervs on all workers.
    elocs .-= mean(elocs)  # Mean over all eloc samples on all workers.
    mul!(f, forces_w, VarDervs, elocs, NetSettings)

    cg_iters = my_cg!(CGP, VarDervs, OptParams, γ, NetSettings, tol=1e-3)
    return cg_iters
end

"""
    optimize!(CGP, VarDervs, OptParams, elocs, forces_w, γ[, iterative])

Optimize the network parameters in-place with Stochastic Reconfiguration using
the EXPLICITLY BUILT S-matrix.

# Arguments

- `CGP::CGPARAMS`: A conjugate gradient parameter struct for when Conjugate
Gradient is performed instead of explicit inversion of the S-matrix.
- `VarDervs::VARDERIVS`: contains sampled variational derivatives and
regularisation parameter.
- `OptParams::OPTPARAMS`: contains optimization variables: S, S^(-1), forces,
wave function parameter update vector ΔW.
- `elocs`: Each column stores the sampled local energies from a worker.
- `forces_w`: Stores the result of adjoint(O) * E (both stanardized) of each
worker in a column, the summed forces are stored in `OptParams.f`.
- `γ::Float64`: learning parameter
- `NetSettings`: network (optimization/sampling) settings

"""
function optimize!(CGP::CGPARAMS, VarDervs::VARDERIVS, OptParams::OPTPARAMS,
    elocs, forces_w, γ, NetSettings)
    @unpack iterative_inverse, mc_trials, dim_rbm, regulator = NetSettings
    @unpack f, s, dw_tot = OptParams
    @unpack o_tot = VarDervs
    o_tot .-= mean(mean(o_tot, dims=2), dims=3)
    elocs .-= mean(elocs)

    mul!(f, forces_w, VarDervs, elocs, NetSettings)

    mul!(s, reshape(o_tot, dim_rbm, mc_trials * nworkers()),
        adjoint(reshape(o_tot, dim_rbm, mc_trials * nworkers())) )
    lmul!(1 / (mc_trials * nworkers()), s)

    if iterative_inverse
        my_cg!(CGP,
            VarDervs.reg * Matrix{Complex{Float64}}(I, dim_rbm, dim_rbm) + s,
            OptParams,
            γ,
            NetSettings,
            tol=1e-3)
    else
        inv_s = zeros(Complex{Float64}, 0, 0)
        try
            inv_s = inv(Hermitian(VarDervs.reg *
                Matrix{Complex{Float64}}(I, dim_rbm, dim_rbm) + s))
        catch SingularException
            println("increased λ")
            VarDervs.reg = 2 * VarDervs.reg
            inv_s = inv(Hermitian(VarDervs.reg *
                Matrix{Complex{Float64}}(I, dim_rbm, dim_rbm) + s))
            VarDervs.reg = regulator  # Revert to original in the case λ was raised.
        end
        dw_tot .= -γ * inv_s * f
    end
    return ""
end

"""
    revert_netparams!(NetParams, OptParams, NetSettings)

Revert the network parameters using the previous shift ΔW saved in
OptParams.dw_tot. Run this function after NaNs occur.
"""
@views function revert_netparams!(NetParams, OptParams, NetSettings)
    @unpack a, b, w = NetParams
    @unpack dw_tot = OptParams
    @unpack n, m = NetSettings
    a .-= dw_tot[1 : n]
    b .-= dw_tot[n + 1: n + m]
    w .-= reshape(dw_tot[n + m + 1: end], n, m)
end

"""
    update_netparams!(NetParams, OptParams, NetSettings)

Update the network parameters using the proposed network shift ΔW saved in
OptParams.dw_tot. Run this function after optimization.
"""
@views function update_netparams!(NetParams, OptParams, NetSettings)
    @unpack a, b, w = NetParams
    @unpack dw_tot = OptParams
    @unpack n, m = NetSettings
    a .+= dw_tot[1 : n]
    b .+= dw_tot[n + 1: n + m]
    w .+= reshape(dw_tot[n + m + 1: end], n, m)
end

"""
    init_optim(OptParamStruct)

Takes a `OPTPARAMS` struct as an argument. Consequently initializes the optimization
parameters (WITH S-matrix).
"""
function init_optim(OptParamStruct, NetSettings)
    @unpack dim_rbm = NetSettings
    s = zeros(Complex{Float64}, dim_rbm, dim_rbm)
    f = zeros(Complex{Float64}, dim_rbm)
    dw_tot = zeros(Complex{Float64}, dim_rbm)
    optimparameters_inst = OptParamStruct(s, f, dw_tot, inv_s)
    return optimparameters_inst
end

"""
    init_optim_mfree(OptParamStruct)

Takes an `OPTPARAMSMFREE` struct as argument. Consequently initializes the
optimization parameters (WITHOUT S-matrix).
"""
function init_optim_mfree(OptParamStruct, NetSettings)
    @unpack dim_rbm = NetSettings
    f = zeros(Complex{Float64}, dim_rbm)
    dw_tot = zeros(Complex{Float64}, dim_rbm)
    optparam_inst = OptParamStruct(f, dw_tot)
end

"""
    distributed_o_eloc(forces_w, VarDervs, elocs)

Writes the part of adjoint(O)⋅E that is sampled by worker number `id`
to a SharedArray `forces_w`. See also `mul!`.
"""
function distributed_o_eloc(forces_w, VarDervs, elocs, id)
    mul!(view(forces_w, :, id),
        view(VarDervs.o_tot, :, :, id),
        view(elocs, :, id))
end

"""
    mul!(forces, forces_w, VarDervs, elocs)

Overloaded function for the calculation of
F = (mc_trials * nworkers)^(-1) * adjoint(O) ⋅ Eloc,  sums over all workers.
"""
function LinearAlgebra.mul!(forces, forces_w, VarDervs::VARDERIVS,
    elocs::SharedArray{Complex{Float64}, 2}, NetSettings)
    @unpack mc_trials = NetSettings
    if nworkers() > 30
        @sync for (id, worker) in enumerate(workers())
            @async remotecall_wait(distributed_o_eloc,
                worker,
                forces_w,
                VarDervs,
                elocs,
                id)
            end
    else
        for id=1:nworkers()
            mul!(view(forces_w, :, id),
                view(VarDervs.o_tot, :, :, id),
                view(elocs, :, id))
            end
    end
    forces .= dropdims(sum(forces_w, dims=2), dims=2)
    # overwrites v with adjoint(O) * (O * b)
    lmul!(1 / (mc_trials * nworkers()), forces)
    # overwrites v with scaled v
end


text:
# Test extensions to basic Base methods
@testset "Base Extensions" begin
    # initialize models and references
    m = InfiniteModel()
    m2 = InfiniteModel()
    ivref = InfiniteVariableRef(m, 1)
    pvref = PointVariableRef(m, 2)
    hvref = HoldVariableRef(m, 3)
    pref = ParameterRef(m, 1)
    # variable compare
    @testset "(==)" begin
        @test ivref == ivref
        @test pvref == pvref
        @test hvref == hvref
        @test ivref == InfiniteVariableRef(m, 1)
        @test pvref == PointVariableRef(m, 2)
        @test hvref == HoldVariableRef(m, 3)
        @test !(ivref == InfiniteVariableRef(m, 2))
        @test !(ivref == InfiniteVariableRef(m2, 1))
        @test !(ivref != InfiniteVariableRef(m, 1))
        @test !(pref == ivref)
    end
    # copy(v)
    @testset "copy(v)" begin
        @test copy(ivref) == ivref
        @test copy(pvref) == pvref
        @test copy(hvref) == hvref
    end
    # copy(v, m)
    @testset "copy(v, m)" begin
        @test copy(ivref, m2) == InfiniteVariableRef(m2, 1)
        @test copy(pvref, m2) == PointVariableRef(m2, 2)
        @test copy(hvref, m2) == HoldVariableRef(m2, 3)
    end
    # broadcastable
    @testset "broadcastable" begin
        @test isa(Base.broadcastable(ivref), Base.RefValue{InfiniteVariableRef})
        @test isa(Base.broadcastable(pvref), Base.RefValue{PointVariableRef})
        @test isa(Base.broadcastable(hvref), Base.RefValue{HoldVariableRef})
    end
end

# Test core JuMP methods
@testset "Core JuMP Extensions" begin
    # initialize models and references
    m = InfiniteModel()
    m2 = InfiniteModel()
    ivref = InfiniteVariableRef(m, 1)
    pvref = PointVariableRef(m, 2)
    hvref = HoldVariableRef(m, 3)
    pref = ParameterRef(m, 1)
    # isequal_canonical
    @testset "JuMP.isequal_canonical" begin
        @test isequal_canonical(ivref, ivref)
        @test isequal_canonical(pvref, pvref)
        @test isequal_canonical(hvref, hvref)
        @test !isequal_canonical(ivref, InfiniteVariableRef(m2, 1))
        @test !isequal_canonical(ivref, InfiniteVariableRef(m, 2))
    end
    # variable_type(m)
    @testset "JuMP.variable_type(m)" begin
        @test variable_type(m) == GeneralVariableRef
    end
    # variable_type(m, t)
    @testset "JuMP.variable_type(m, t)" begin
        @test variable_type(m, Infinite) == InfiniteVariableRef
        @test variable_type(m, Point) == PointVariableRef
        @test variable_type(m, Hold) == HoldVariableRef
        @test variable_type(m, Parameter) == ParameterRef
        @test_throws ErrorException variable_type(m, :bad)
    end
end

# Test precursor functions needed for add_parameter
@testset "Basic Reference Queries" begin
    # initialize model and infinite variable
    m = InfiniteModel()
    ivref = InfiniteVariableRef(m, 1)
    info = VariableInfo(false, 0, false, 0, false, 0, false, 0, false, false)
    param = InfOptParameter(IntervalSet(0, 1), Number[], false)
    pref = add_parameter(m, param, "test")
    m.vars[1] = InfiniteVariable(info, (pref, ))
    # JuMP.index
    @testset "JuMP.index" begin
        @test JuMP.index(ivref) == 1
    end
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(ivref) == m
    end
    # JuMP.is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, ivref)
        @test !is_valid(InfiniteModel(), ivref)
        @test !is_valid(m, InfiniteVariableRef(m, 5))
    end
end


text:
using Documenter
using InMemoryDatasets

# DocMeta.setdocmeta!(InMemoryDatasets, :DocTestSetup, :(using InMemoryDatasets); recursive=true)

# Build documentation.
# ====================

makedocs(
    # options
    # modules = [InMemoryDatasets],
    doctest = false,
    clean = false,
    sitename = "InMemoryDatasets",
    # format = Documenter.HTML(
    #     canonical = "https://sl-solution.github.io/InMemoryDataset.jl/stable/",
    #     edit_link = "main"
    # ),
    pages = Any[
        "Introduction" => "index.md",
        "First Steps" => "man/basics.md",
        "Tutorial" => "man/tutorial.md",
        "User Guide" => Any[
            "Missing Values" => "man/missing.md",
            "Formats" => "man/formats.md",
            "Call functions on each observation" => "man/map.md",
            "Row-wise operations" => "man/byrow.md",
            "Transform columns" => "man/modify.md",
            "Filter observations" => "man/filter.md",
            "Sort" => "man/sorting.md",
            "Group observations" => "man/grouping.md",
            "Aggregation" => "man/aggregation.md",
            "Transpose data" => "man/transpose.md",
            "Joins" => "man/joins.md"
        ],
        "Gallery" => "man/gallery.md",
        "Performance tips" => "man/performance.md"
        # "API" => Any[
        #     "Functions" => "lib/functions.md"
        # ]
    ],
    strict = true
)

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    # options
    repo = "github.com/sl-solution/InMemoryDatasets.jl",
    target = "build",
    deps = nothing,
    make = nothing,
    devbranch = "master"
)


text:
##########################################################################################
# Project directory
##########################################################################################
export projectdir, datadir, srcdir, plotsdir, scriptdir, papersdir
export projectname
export findproject, quickactivate

"""
    projectdir()
Return the directory of the currently active project. Ends with `"/"`.

```julia
projectdir(folder::String) = joinpath(projectdir(), folder)*"/"
```
Return the directory of the `folder` in the active project.
"""
projectdir() = dirname(Base.active_project())*"/"
projectdir(folder::String) = joinpath(projectdir(), folder)*"/"

datadir() = projectdir()*"data/"
srcdir() = projectdir()*"src/"
plotsdir() = projectdir()*"plots/"
scriptdir() = projectdir()*"scripts/"
papersdir() = projectdir()*"papers/"
testdir() = projectdir()*"test/"

"""
    projectname()
Return the name of the currently active project.
"""
projectname() = Pkg.REPLMode.promptf()[2:end-7]

"""
    findproject(path = pwd()) -> project_path
Recursively search `path` and its parents for a valid Julia project file
(anything in `Base.project_names`).
If it is found return its path, otherwise issue a warning and return
`nothing`.

The function stops searching if it hits either the home directory or
the root directory.
"""
function findproject(dir::AbstractString = pwd())
    # look for project file in current dir and parents
    home = homedir()
    while true
        for proj in Base.project_names
            file = joinpath(dir, proj)
            Base.isfile_casesensitive(file) && return dir
        end
        # bail at home directory
        dir == home && break
        # bail at root directory
        old, dir = dir, dirname(dir)
        dir == old && break
    end
    @warn "Could not find find a project file by recursively checking "*
    "given `path` and its parents. Returning `nothing` instead."
    return nothing
end

"""
    quickactivate(path [, name::String])
Activate the project found by [`findproject`](@ref) of the `path`, which
recursively searches the `path` and its parents for a valid Julia
project file.

Optionally check if `name` is the same as the activated project's name.
If it is not, throw an error.

This function is _first_ activating the project and _then_ checking if
it matches the `name`.

!!! warning
    Note that to access `quickactivate` you need to be `using DrWatson`.
    For this to be possible `DrWatson` must be already added in the
    existing global environment. The version of `DrWatson` loaded therefore
    will be the one of the global environment, and not of the activated project.
    To avoid unexpected behavior take care so that these two versions coincide.

    **In addition please be very careful to not write:**
    ```julia
    using DrWatson, Package1, Package2
    quickactivate(@__DIR__)
    # do stuff
    ```
    **but instead load packages after activating the project:**
    ```julia
    using DrWatson
    quickactivate(@__DIR__)
    using Package1, Package2
    # do stuff
    ```
    This ensures that the packages you use will all have the versions dictated
    by your activated project (besides `DrWatson`, since this is impossible
    to do using `quickactivate`).
"""
function quickactivate(path, name = nothing)
    projectpath = findproject(path)
    projectpath === nothing && return nothing
    Pkg.activate(projectpath)
    if !(name === nothing) && projectname() != name
        error(
        "The activated project did not match asserted name. Current project "*
        "name is $(projectname()) while the asserted name is $name."
        )
    end
    return nothing
end

##########################################################################################
# Project setup
##########################################################################################
export initialize_project

const DEFAULT_PATHS = [
"_research", "src/", "scripts/",
"plots/", "notebooks/",
"papers/",
"test",
"data/sims/",
"data/exp_raw/",
"data/exp_pro/",
]
const DEFAULT_README = """
This is an awesome new scientific project that uses `DrWatson`!\n
"""

"""
    initialize_project(path [, name]; kwargs...)
Initialize a scientific project expected by `DrWatson` inside the given `path`.
If its `name` is not given, it is assumed to be the folder's name.

The new project remains activated for you to immidiately add packages.

## Keywords
* `readme = true` : adds a README.md file.
* `authors = nothing` : if a string or container of strings, adds the authors in the
  Project.toml file.
* `force = false` : If the `path` is _not_ empty then throw an error. If however `force`
  is `true` then recursively delete everything in the path and create the project.
* `git = true` : Make the project a Git repository.
"""
function initialize_project(path, name = basename(path);
    force = false, readme = true, authors = nothing,
    git = true)

    mkpath(path)
    rd = readdir(path)
    if length(rd) != 0
        if force
            for d in rd
                rm(joinpath(path, d), recursive = true, force = true)
            end
        else
            error("Project path is not empty!")
        end
    end

    if git; repo = LibGit2.init(path); end
    git && LibGit2.commit(repo, "Initial commit")
    Pkg.activate(path)
    try
        Pkg.add("DrWatson")
    catch
        @warn "Could not add DrWatson to project. Adding Pkg instead..."
        Pkg.add("Pkg")
    end
    # Default folders
    for p in DEFAULT_PATHS
        mkpath(joinpath(path, p))
    end

    git && LibGit2.add!(repo, "Project.toml")
    git && LibGit2.add!(repo, "Manifest.toml")
    git && LibGit2.add!(repo, DEFAULT_PATHS...)
    git && LibGit2.commit(repo, "Folder setup by DrWatson")

    # Default files
    cp(joinpath(@__DIR__, "defaults", "gitignore.txt"), joinpath(path, ".gitignore"))
    cp(joinpath(@__DIR__, "defaults", "intro.jl"), joinpath(path, "scripts/intro.jl"))
    cp(joinpath(@__DIR__, "defaults", "runtests.jl"), joinpath(path, "test/runtests.jl"))
    files = vcat(".gitignore", "/scripts/intro.jl","/test/runtests.jl")
    if readme
        write(joinpath(path, "README.md"), DEFAULT_README)
        push!(files, "README.md")
    end
    pro = read(joinpath(path, "Project.toml"), String)
    w = "name = \"$name\"\n"
    if !(authors === nothing)
            w *= "authors = "*sprint(show, vecstring(authors))*"\n"
    end
    write(joinpath(path, "Project.toml"), w, pro)
    push!(files, "Project.toml")

    git && LibGit2.add!(repo, files...)
    git && LibGit2.commit(repo, "File setup by DrWatson")
    return path
end

vecstring(a::String) = [a]
vecstring(a::Vector{String}) = a
vecstring(c) = [string(a) for a in c]

function greet()
    println(
    """
    DrWatson is currently in beta.
    Help us make it better by opening
    issues on GitHub or submitting feature requests!

    Currently active project is: $(projectname())

    Have fun with your new project!
    """
    )
end


text:
module InteractBulma

import InteractBase, Sass

struct BulmaTheme<:InteractBase.WidgetTheme
    path::String
end

const examplefolder = joinpath(@__DIR__, "..", "examples")
const main_css = joinpath(@__DIR__, "..", "assets", "main.min.css")
const main_confined_css = joinpath(@__DIR__, "..", "assets", "main_confined.min.css")
const font_awesome = InteractBase.font_awesome

function InteractBase.libraries(b::BulmaTheme)
    bulmalib = joinpath(b.path, InteractBase.isijulia() ? "main_confined.min.css" : "main.min.css")
    vcat(font_awesome, InteractBase.style_css, bulmalib)
end

const main_scss = joinpath(@__DIR__, "..", "assets", "main.scss")
const main_confined_scss = joinpath(@__DIR__, "..", "assets", "main_confined.scss")
const _overrides = joinpath(@__DIR__, "..", "assets", "_overrides")
const _variables = joinpath(@__DIR__, "..", "assets", "_variables")

function copy_or_empty(::Nothing, dest)
    dest = dest*".scss"
    open(io -> nothing, dest, "w")
    return dest
end
function copy_or_empty(src::AbstractString, dest)
    extension = split(src, '.')[end]
    dest = "$dest.$extension"
    cp(realpath(abspath(src)), dest, force = true)
    return dest
end

function compile_theme(output = mktempdir(); overrides = nothing, variables = nothing)
    asset_dir = joinpath(@__DIR__, "..", "assets")
    for file in readdir(asset_dir)
        startswith(file, "_") && rm(joinpath(asset_dir, file))
    end
    copy_or_empty(overrides, _overrides)
    copy_or_empty(variables, _variables)
    Sass.compile_file(main_scss, joinpath(output, "main.min.css"), output_style = Sass.compressed)
    Sass.compile_file(main_confined_scss, joinpath(output, "main_confined.min.css"), output_style = Sass.compressed)
    return BulmaTheme(output)
end

end # module


text:
let sys = DenseSystem(test_params_flat...),
    max_iter = 100,
    decomp = IterativeDecomposition(sys, max_iter),
    sys_decomp = DiagonalSystem(decomp)

    @test decomp.degens == [1]

    # Remove overall phase.
    vs_decomp = decomp.vs[1,1] > 0 ? decomp.vs : -decomp.vs
    vs_expected = reshape([0.5833168784, 0.5773295720, 0.5713422656], 3, 1)
    @test isapprox(vs_decomp, vs_expected)

    expected_energy = reshape([-2.880208314], 1, 1)
    expected_freq = reshape([1.0, 2.0], test_M, 1)
    expected_lin = reshape([0.0149574448, 0.0179572297], test_M, 1, 1)
    expected_quad = reshape([0.0179572297 0.0209570146; 0.0209570146 0.0239567995], test_M, test_M, 1, 1)
    sys_expected = DiagonalSystem(expected_energy, expected_freq, expected_lin, expected_quad)
    @test isapprox(sys_decomp, sys_expected)
end

let sys = DenseSystem(test_params_sym...),
    max_iter = 100,
    decomp = IterativeDecomposition(sys, max_iter)

    @test decomp.degens == [1, 2]
end

let sys = DenseSystem(test_params_sym3...),
    max_iter = 100

    @test_throws DegeneracyException IterativeDecomposition(sys, max_iter)
end


text:
include("QuadForm/Basic.jl")
include("QuadForm/Lattices.jl")
include("QuadForm/Genus.jl")


text:
using OrthogonalPolynomialsQuasi, BandedMatrices, LazyArrays, Test

@testset "Ultraspherical" begin
    @testset "Transforms" begin
        U = Ultraspherical(1)
        x = axes(U,1)
        Un = U[:,Base.OneTo(5)]
        @test factorize(Un) isa ContinuumArrays.TransformFactorization
        @test (Un \ x) ≈ [0,0.5,0,0,0]
        @test (U * (U \ exp.(x)))[0.1] ≈ exp(0.1)
    end

    @testset "Operators" begin
        @testset "Lowering" begin
            λ = 1
            wC1 = WeightedUltraspherical(λ)
            wC2 = WeightedUltraspherical(λ+1)
            L = wC1 \ wC2
            @test L isa BandedMatrix
            @test bandwidths(L) == (2,0)
            @test L[1,1] ≈ 2λ*(2λ+1)/(4λ*(λ+1))
            @test L[3,1] ≈ -2/(4λ*(λ+1))
            u = [randn(10); zeros(∞)]
            @test wC1[0.1,:]'*(L*u) ≈ wC2[0.1,:]'*u
        end

        @testset "Weighted Derivative" begin
            T = Chebyshev()
            wC1 = WeightedUltraspherical(1)
            wC2 = WeightedUltraspherical(2)
            x = axes(wC2,1)
            D = Derivative(x)

            @test wC1 \ (D*wC2) isa BandedMatrix
            @test (wC1 \ (D*wC2))[1:5,1:5] == BandedMatrix(-1 => [-1.5,-4,-7.5,-12])

            u = wC2 * [randn(5); zeros(∞)]
            @test (D*u)[0.1] ≈ ((D*T) * (T\u))[0.1]
        end

        @testset "Interrelationships" begin
            @testset "Chebyshev–Ultrashperical" begin
                T = Chebyshev()
                U = ChebyshevU()
                C = Ultraspherical(2)
                D = Derivative(axes(T,1))

                @test C\C === pinv(C)*C === Eye(∞)
                D₀ = U\(D*T)
                D₁ = C\(D*U)
                @test D₁ isa BandedMatrix
                @test (D₁*D₀)[1:10,1:10] == diagm(2 => 4:2:18)
                @test D₁*D₀ isa MulMatrix
                @test bandwidths(D₁*D₀) == (-2,2)

                S₁ = (C\U)[1:10,1:10]
                @test S₁ isa BandedMatrix{Float64}
                @test S₁ == diagm(0 => 1 ./ (1:10), 2=> -(1 ./ (3:10)))
            end
            @testset "Legendre" begin
                @test Ultraspherical(0.5) \ (UltrasphericalWeight(0.0) .* Ultraspherical(0.5)) == Eye(∞)
                @test Legendre() \ (UltrasphericalWeight(0.0) .* Ultraspherical(0.5)) == Eye(∞)
            end
        end
    end

    @testset "test on functions" begin
        T = Chebyshev()
        U = Ultraspherical(1)
        D = Derivative(axes(T,1))
        f = T*Vcat(randn(10), Zeros(∞))
        @test (U*(U\f)).args[1] isa Ultraspherical
        @test (U*(U\f))[0.1] ≈ f[0.1]
        @test (D*f)[0.1] ≈ ForwardDiff.derivative(x -> (ChebyshevT{eltype(x)}()*f.args[2])[x],0.1)
    end

    @testset "Evaluation" begin
        C = Ultraspherical(2)
        @test @inferred(C[0.1,Base.OneTo(0)]) == Float64[]
        @test @inferred(C[0.1,Base.OneTo(1)]) == [1.0]
        @test @inferred(C[0.1,Base.OneTo(2)]) == [1.0,0.4]
        @test @inferred(C[0.1,Base.OneTo(3)]) == [1.0,0.4,-1.88]
    end

end


text:
macro f_str(str)
    inbetweens = String[]
    arguments  = String[]
    formatstrs = String[]
    formatstr = ""
    nextidx = 1
    while formatstr !== nothing
        nextidx, inbetween, formatstr = JuFormatting.iterate_braces_pair(str, nextidx)

        push!(inbetweens, inbetween)
        if formatstr !== nothing
            temp = split(formatstr, ":")
            if temp[1] == ""
                error("f_str: empty expression not allowed")
            end
            push!(arguments, temp[1])
            if length(temp) > 1
                push!(formatstrs, temp[2])
            else
                push!(formatstrs, "")
            end
        end
    end

    parsed_args = esc.(Meta.parse.(arguments))

    :(fstr_combine($inbetweens, $formatstrs, $(parsed_args...)))
end

function fstr_combine(inbetweens, fmtstrs, arguments...)
    result = ""
    for i in 1:length(fmtstrs)
        result *= inbetweens[i] * fmt(fmtstrs[i], arguments[i])
    end
    result * inbetweens[end]
end


text:
using Pkg
Pkg.activate(@__DIR__)
using Pluto
Pluto.run()

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
#!/usr/bin/env julia

# https://github.com/JuliaEditorSupport/julia-emacs/blob/master/make-julia-latexsubs.jl
@assert VERSION >= v"1"

import REPL

"""
Create latex symbols formatted for elisp as either abbrev or hash table.
- ds      : elisp output data structure - either "abbrev" or "hash"
- varname : name of hashtable if "hash" output
- dest    : out stream
"""
function generate_unicode(ds::AbstractString, varname::AbstractString, dest::IO)
    if dest != stdout
        println(dest, ";; -*- coding: utf-8; mode: emacs-lisp; -*-")
    end
    if ds == "abbrev"
        println(dest, "(define-abbrev-table 'unicode-latex-abbrev-table\n  '(")
    elseif ds == "hash"
        println(dest, "(defvar $(varname) (make-hash-table :test 'equal))\n")
    end

    for (k, v) in sort!(collect(REPL.REPLCompletions.latex_symbols), by=last)
        ks = escape_string(k)
        vs = escape_string(v)

        if ds == ""
            println(dest, "\"$ks\" => \"$vs\"")
        elseif ds == "hash" && occursin(r"^\\U[0-9A-Fa-f]+$", vs)
            cp = vs[3:end]
            println(dest,
                    "(let ((c (decode-char 'ucs #x$cp)))\n",
                    "  (if c (puthash \"$ks\" (char-to-string c) $(varname))))")
        else
            if ds == "hash"
                println(dest, "(puthash \"$ks\" \"$vs\" $(varname))")
            elseif ds == "abbrev"
                println(dest, "    (\"$ks\" \"$vs\" nil :system t)")
            else
                println(dest, "\"$ks\" => \"$vs\"")
            end
        end
    end

    if ds == "abbrev"
        println(dest,
                "    )\n",
                "  \"Unicode latex abbrevs (julia generated).\"\n",
                "  :regexp \"\\\\(\\\\\\\\[\\\\A-Za-z0-9^]+\\\\)\")")
    end
end

# CLI

if length(ARGS) < 1
    generate_unicode("", "", stdout)
else
    dest = length(ARGS) > 2 ? open(ARGS[3], "w") : stdout
    if ARGS[1] == "abbrev"
        generate_unicode("abbrev", "", dest)
    elseif ARGS[1] == "hash"
        generate_unicode(ARGS[1], ARGS[2], dest)
    else
        println(dest, "Unknown type")
    end
end

max_stars_repo_name: nverno/nvp

max_issues_repo_name: nverno/nvp

max_forks_repo_name: nverno/nvp


text:
facts("Query parsing.") do

    context("Functions.") do

        @fact @query(f) --> Query([f, :f], 0)
        @fact @query(g) --> Query([g, :g], 0)

        @fact @query(f, 1) --> Query([f, :f], 1)
        @fact @query(g, 1) --> Query([g, :g], 1)

        @fact @query(A.f) --> Query([A.f, :f], A, 0)
        @fact @query(A.g) --> Query([A.g, :g], A, 0)

        @fact @query(A.f, 1) --> Query([A.f, :f], A, 1)
        @fact @query(A.g, 1) --> Query([A.g, :g], A, 1)

        @fact @query(A.B.f) --> Query([A.B.f, :f], A.B, 0)
        @fact @query(A.B.g) --> Query([A.B.g, :g], A.B, 0)

        @fact @query(A.B.f, 1) --> Query([A.B.f, :f], A.B, 1)
        @fact @query(A.B.g, 1) --> Query([A.B.g, :g], A.B, 1)

    end

    context("Methods.") do

        @fact @query(f()) --> Query([@which(f())], 0)
        @fact @query(g()) --> Query([@which(g())], 0)

        @fact @query(f(1)) --> Query([@which(f(1))], 0)
        @fact @query(g(1)) --> Query([@which(g(1))], 0)

        @fact @query(f(), 1) --> Query([@which(f())], 1)
        @fact @query(g(), 1) --> Query([@which(g())], 1)

        @fact @query(f(1), 1) --> Query([@which(f(1))], 1)
        @fact @query(g(1), 1) --> Query([@which(g(1))], 1)

        @fact @query(A.f()) --> Query([@which(A.f())], A, 0)
        @fact @query(A.g()) --> Query([@which(A.g())], A, 0)

        @fact @query(A.f(1)) --> Query([@which(A.f(1))], A, 0)
        @fact @query(A.g(1)) --> Query([@which(A.g(1))], A, 0)

        @fact @query(A.f(), 1) --> Query([@which(A.f())], A, 1)
        @fact @query(A.g(), 1) --> Query([@which(A.g())], A, 1)

        @fact @query(A.f(1), 1) --> Query([@which(A.f(1))], A, 1)
        @fact @query(A.g(1), 1) --> Query([@which(A.g(1))], A, 1)

        @fact @query(A.B.f()) --> Query([@which(A.B.f())], A.B, 0)
        @fact @query(A.B.g()) --> Query([@which(A.B.g())], A.B, 0)

        @fact @query(A.B.f(1)) --> Query([@which(A.B.f(1))], A.B, 0)
        @fact @query(A.B.g(1)) --> Query([@which(A.B.g(1))], A.B, 0)

        @fact @query(A.B.f(), 1) --> Query([@which(A.B.f())], A.B, 1)
        @fact @query(A.B.g(), 1) --> Query([@which(A.B.g())], A.B, 1)

        @fact @query(A.B.f(1), 1) --> Query([@which(A.B.f(1))], A.B, 1)
        @fact @query(A.B.g(1), 1) --> Query([@which(A.B.g(1))], A.B, 1)

    end

    context("Macros.") do

        @fact @query(@m)      --> Query([getfield(LexiconTests, symbol("@m"))], 0)
        @fact @query(@m(), 1) --> Query([getfield(LexiconTests, symbol("@m"))], 1)

        @fact @query(A.@m)      --> Query([getfield(A, symbol("@m"))], A, 0)
        @fact @query(A.@m(), 1) --> Query([getfield(A, symbol("@m"))], A, 1)

        @fact @query(A.B.@m)      --> Query([getfield(A.B, symbol("@m"))], A.B, 0)
        @fact @query(A.B.@m(), 1) --> Query([getfield(A.B, symbol("@m"))], A.B, 1)

    end

    context("Types.") do

        @fact @query(T)     --> Query([T, :T], 0)
        @fact @query(A.T)   --> Query([A.T, :T], A, 0)
        @fact @query(A.B.T) --> Query([A.B.T, :T], A.B, 0)

        @fact @query(T, 1)     --> Query([T, :T], 1)
        @fact @query(A.T, 1)   --> Query([A.T, :T], A, 1)
        @fact @query(A.B.T, 1) --> Query([A.B.T, :T], A.B, 1)

        @fact @query(S)     --> Query([S, :S], 0)
        @fact @query(A.S)   --> Query([A.S, :S], A, 0)
        @fact @query(A.B.S) --> Query([A.B.S, :S], A.B, 0)

        @fact @query(S, 1)     --> Query([S, :S], 1)
        @fact @query(A.S, 1)   --> Query([A.S, :S], A, 1)
        @fact @query(A.B.S, 1) --> Query([A.B.S, :S], A.B, 1)

    end

    context("Globals.") do

        @fact @query(K)     --> Query([K, :K], 0)
        @fact @query(A.K)   --> Query([A.K, :K], A, 0)
        @fact @query(A.B.K) --> Query([A.B.K, :K], A.B, 0)

        @fact @query(K, 1)     --> Query([K, :K], 1)
        @fact @query(A.K, 1)   --> Query([A.K, :K], A, 1)
        @fact @query(A.B.K, 1) --> Query([A.B.K, :K], A.B, 1)

    end
end


text:
abstract type AbstractObject end

struct Agent <: AbstractObject end
const AGENT = Agent()

struct Wall <: AbstractObject end
const WALL = Wall()

struct Goal <: AbstractObject end
const GOAL = Goal()

struct Door <: AbstractObject end
const DOOR = Door()

struct Key <: AbstractObject end
const KEY = Key()

struct Gem <: AbstractObject end
const GEM = Gem()

struct Obstacle <: AbstractObject end
const OBSTACLE = Obstacle()

struct Box <: AbstractObject end
const BOX = Box()

struct Target <: AbstractObject end
const TARGET = Target()

struct Target1 <: AbstractObject end
const TARGET1 = Target1()

struct Target2 <: AbstractObject end
const TARGET2 = Target2()

struct Body <: AbstractObject end
const BODY = Body()

struct Food <: AbstractObject end
const FOOD = Food()

struct Basket <: AbstractObject end
const BASKET = Basket()

struct Ball <: AbstractObject end
const BALL = Ball()


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_energy_market_unit_priority" begin

    function create_instance_with_expensive_wind()
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 150.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_2", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                0., 100.,
                                                0., 15.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_2", PSCOPF.Networks.PILOTABLE,
                                                0., 100.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", TS[1], "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_2", TS[1], "S1", 25.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", 55.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(TS[1] => PSCOPF.FREE),
                                "wind_1_2" => SortedDict(TS[1] => PSCOPF.FREE),
                                "prod_1_1" => SortedDict(TS[1] => PSCOPF.FREE),
                                "prod_1_2" => SortedDict(TS[1] => PSCOPF.FREE), )
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)

        return TS, ech, firmness, context
    end

    #=
    By default, Energy markets picks Limitable units first. Then, the cheapest units.

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=150   |
      S1: 20            |
                        |
    (limitable) wind_1_2|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |
      S1: 25            |
                        |
    (pilotable) prod_1_1|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=15    |
                        |
    (pilotable) prod_1_2|
     Pmin=0, Pmax=100   |
     Csta=0, Cprop=10   |
                        |

    Demand : 55
    Available limitable power : 45

    The units prod_1_1 and prod_1_2 are cheaper than wind_1_2.

    The market selects the limitable units first.
    So, it uses wind_1_2 before prod_1_1 and prod_1_2 even if their proportional cost is lower.

    => wind_1_1 = 20 and wind_1_2 = 25
    remaining demand : 10MW
    => prod_1_1 = 10 (because prod_1_1 is cheaper than prod_1_2)
    =#
    @testset "energy_market_picks_renewables_first" begin
        TS, ech, firmness, context = create_instance_with_expensive_wind()

        market = PSCOPF.EnergyMarket()

        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test 25. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_2", TS[1], "S1")
        @test 10. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09
        # pmin = 0 for prod_1_1 & prod_1_2
        @test ismissing( PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1") )
        @test ismissing( PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1") )
    end

    #=
    In S1:
      load = 55, wind=45 => missing 10
      prod_1_1 has a pmin=20 => we produce 20 costing 200 (20MW*10)
                                + we cap 10MW of wind
      prod_1_2 has a pmin=5 => we produce 10 costing 150 (10MW*15)
    => use prod_1_2
    In S2:
      load = 60, wind=45 => missing 15
      prod_1_1 has a pmin=20 => we produce 20 costing 200 (20MW*10)
                                + we cap 5MW of wind
      prod_1_2 has a pmin=5 => we produce 15 costing 225 (15MW*15)
    => use prod_1_1
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load_2
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=1     | S2: 60
      S1: 45            |
      S2: 45            |
                        |
    (pilotable) prod_1_1|
    Pmin=20, Pmax=100   |
    Csta=0, Cprop=10    |
                        |
    (pilotable) prod_1_2|
     Pmin=5, Pmax=100   |
     Csta=0, Cprop=15   |
                        |
    =#
    @testset "energy_market_may_cap_due_to_pilotables_pmin" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables : have a Pmin but no start cost
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_2", PSCOPF.Networks.PILOTABLE,
                                                5., 100.,
                                                0., 15.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 45.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S2", 45.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S2", 60.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            ),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                    )
        # initial generators state
        generators_init_state = SortedDict(
            "prod_1_1" => PSCOPF.ON,
            "prod_1_2" => PSCOPF.ON
        )

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        market = PSCOPF.EnergyMarket()
        result = PSCOPF.run(market, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_market_schedule!(context, ech, result, firmness, market)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        #S1
        @test 45. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S1")
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S1")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S1") < 1e-09 #cheapest but pmin=20
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test 10. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S1")
        @test value(result.limitable_model.p_global_capping[TS[1], "S1"]) < 1e-09

        #S2
        @test 40. ≈ PSCOPF.get_prod_value(context.market_schedule, "wind_1_1", TS[1], "S2")
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_1", TS[1], "S2")
        @test 20. ≈ PSCOPF.get_prod_value(context.market_schedule, "prod_1_1", TS[1], "S2")
        @test PSCOPF.OFF == PSCOPF.get_commitment_value(context.market_schedule, "prod_1_2", TS[1], "S2")
        @test PSCOPF.get_prod_value(context.market_schedule, "prod_1_2", TS[1], "S2") < 1e-09
        @test 5. ≈ value(result.limitable_model.p_global_capping[TS[1], "S2"])

        @test value(result.objective_model.start_cost) < 1e-09
        @test value(result.objective_model.prop_cost) ≈ (
              (0. * 10. + 10. * 15) #S1
            + (20. * 10. + 0. * 15) #S2
        )
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
@with_kw mutable struct OneIntent
    playerintent::PLAYERINTENT
    keybutton::AbstractPlotting.Keyboard.Button
    scale::Float64 = 1.0
    intended::Bool = false
end

struct NoInputHandler <: AbstractInputHandler end

struct InputHandler <: AbstractInputHandler
    currentintent::Vector{OneIntent} # does this need to be a parameter for speed?
    playerintentindex::Dict{PLAYERINTENT, Int64}
    callback::Observables.Observable{Nothing}
end

numtype = Float64
movementscale = numtype(0.2)
zoomscale = numtype(1.0)
rotscale = numtype(1.0)

defaultkeyboardmapping() = [
                            OneIntent(playerintent=moveforward, keybutton=Makie.Keyboard.w, scale=movementscale),
                            OneIntent(playerintent=movebackward, keybutton=Makie.Keyboard.s, scale=movementscale),
                            OneIntent(playerintent=moveleft, keybutton=Makie.Keyboard.a, scale=movementscale),
                            OneIntent(playerintent=moveright, keybutton=Makie.Keyboard.d, scale=movementscale),
                            OneIntent(playerintent=zoomout, keybutton=Makie.Keyboard.j, scale=zoomscale),
                            OneIntent(playerintent=zoomin, keybutton=Makie.Keyboard.k, scale=zoomscale),
                            OneIntent(playerintent=rotateleft, keybutton=Makie.Keyboard.q, scale=rotscale),
                            OneIntent(playerintent=rotateright, keybutton=Makie.Keyboard.e, scale=rotscale),
                            OneIntent(playerintent=plantlighthouse, keybutton=Makie.Keyboard.l),
                            OneIntent(playerintent=stoprenderloop, keybutton=Makie.Keyboard.p),
                            ]

function generate_input_handler(events::AbstractPlotting.Events, currentintent::Vector{OneIntent})

    # input processor
    callback = Makie.lift(events.keyboardbuttons) do but
        #@show but
        for i in eachindex(currentintent)
            @unpack keybutton = currentintent[i]
            currentintent[i].intended = keybutton in but
        end
        #@show currentintent
        nothing
    end

    playerintentindex = Dict([oneindex.playerintent => i for (i, oneindex) in enumerate(currentintent)])

    @show typeof(callback), callback
    handler = InputHandler(currentintent, playerintentindex, callback)

    handler
end

true

#=

    Makie.lift(image.events.keyboardbuttons) do but
        #@show but

        wdown = Makie.Keyboard.w in but
        adown = Makie.Keyboard.a in but
        sdown = Makie.Keyboard.s in but
        ddown = Makie.Keyboard.d in but

        jdown = Makie.Keyboard.j in but
        kdown = Makie.Keyboard.k in but

        udown = Makie.Keyboard.u in but
        idown = Makie.Keyboard.i in but

        endscene = Makie.Keyboard.left_control in but

        nothing
    end
        expscale = numtype(2.0 ^ -scalefactor)
        yrangeextent = numtype(expscale)
        xrangeextent = numtype(expscale * aspectratio)
        xstart = numtype(-xrangeextent / 2)
        ystart = numtype(-yrangeextent / 2)

        cosrot = numtype(cos(rotation))
        sinrot = numtype(sin(rotation))

        if wdown
            centerxlocal += -sinrot * movementscale * expscale * deltatime
            centerylocal += cosrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end
        if adown
            centerxlocal -= cosrot * movementscale * expscale * deltatime
            centerylocal -= sinrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end
        if sdown
            centerxlocal -= -sinrot * movementscale * expscale * deltatime
            centerylocal -= cosrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end
        if ddown
            centerxlocal += cosrot * movementscale * expscale * deltatime
            centerylocal += sinrot * movementscale * expscale * deltatime
            #@show centerxlocal
            #@show centerylocal
        end

        # note: scale factor for Float32 is approximately capped at 32
        if jdown
            scalefactor -= zoomscale * deltatime
            expscale = numtype(2.0 ^ -scalefactor)
            yrangeextent = numtype(expscale)
            xrangeextent = numtype(expscale * aspectratio)
            xstart = numtype(-xrangeextent / 2)
            ystart = numtype(-yrangeextent / 2)
            #@show scalefactor
        end
        if kdown
            scalefactor += zoomscale * deltatime
            expscale = numtype(2.0 ^ -scalefactor)
            yrangeextent = numtype(expscale)
            xrangeextent = numtype(expscale * aspectratio)
            xstart = numtype(-xrangeextent / 2)
            ystart = numtype(-yrangeextent / 2)
            #@show scalefactor
        end

        if udown
            rotation -= rotscale * deltatime
            cosrot = numtype(cos(rotation))
            sinrot = numtype(sin(rotation))
            #@show rotation
        end
        if idown
            rotation += rotscale * deltatime
            cosrot = numtype(cos(rotation))
            sinrot = numtype(sin(rotation))
            #@show rotation
        end

        if curtime - mobiuswalktime > 3
            lambdaangletarget, adisktarget = randomdiskmobius()
            mobiuswalktime = curtime
        end

=#

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta