text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds - pODEs Discrete
# A package for discretize and relax methods for bounding pODEs.
# See https://github.com/PSORLab/DynamicBoundspODEsDiscrete.jl
#############################################################################
# src/DiscretizeRelax/utilities/relax.jl
# Defines the relax! routine for the DiscretizeRelax integrator.
#############################################################################

function DBB.relax!(d::DiscretizeRelax{M,T,S,F,K,X,NY}) where {M <: AbstractStateContractor, T <: Number, S <: Real, F, K, X, NY}

    fill!(d.time, 0.0)

    storage_len = length(d.storage)
    for i=1:(1000-storage_len)
        push!(d.storage, zeros(T, d.nx))
        push!(d.storage_apriori, zeros(T, d.nx))
        push!(d.time, 0.0)
    end

    # reset relax!
    for i = 1:length(d.storage)
        fill!(d.storage[i], zero(T))
    end
    for i = 1:length(d.storage_apriori)
        fill!(d.storage_apriori[i], zero(T))
    end
    empty!(d.relax_t_dict_indx)
    empty!(d.relax_t_dict_flt)
    d.step_result.time = 0.0

    set_P!(d)         # Functor set P and P - p values for calculations
    compute_X0!(d)     # Compute initial condition values

    # Get initial time and integration direction
    t = d.tspan[1]
    tmax = d.tspan[2]
    sign_tstep = copysign(1, tmax - t)
    d.time[1] = t
    d.contractor_result.times[1] = t

    # Computes maximum step size to take (either hit support or max time)
    support_indx = 1
    next_support = Inf
    tsupports = d.tsupports
    if !isempty(tsupports)
        if (tsupports[1] == 0.0)
            support_indx += 1
        end
        next_support = tsupports[support_indx]
    end

    # initialize QR type storage
    set_Δ!(d.contractor_result.Δ, d.storage)
    fill!(d.contractor_result.A_Q[1], 0.0)
    fill!(d.contractor_result.A_inv[1], 0.0)
    for i = 1:d.nx
        d.contractor_result.A_Q[1][i,i] = 1.0
        d.contractor_result.A_inv[1][i,i] = 1.0
    end

    d.exist_result.status_flag = RELAXATION_NOT_CALLED

    # Begin integration loop
    hlast = 0.0
    d.step_count = 0
    is_adaptive = d.exist_result.hj <= 0.0
    d.exist_result.hj = !is_adaptive ? d.exist_result.hj : 0.01*(tmax - d.contractor_result.times[1])
    d.exist_result.predicted_hj = d.exist_result.hj
    stored_value_count = length(d.relax_t_dict_indx)
    for step_number = 2:(d.step_limit+2)
        if (sign_tstep*d.step_result.time <= sign_tstep*tmax) &&
           (d.exist_result.predicted_hj != 0.0)

            delT = abs(sign_tstep*d.step_result.time - sign_tstep*tmax)

            # max step size is min of predicted, when next support point occurs,
            # or the last time step in the span
            tv = d.step_result.time
            d.exist_result.hj = min(d.exist_result.hj, next_support - tv, tmax - tv)
            hj_limit = min(next_support - tv, tmax - tv)
            d.exist_result.hj_max = tmax - tv
            is_support_pnt = (hj_limit == next_support - tv)

            d.contractor_result.steps[1] = d.exist_result.hj
            d.contractor_result.step_count = step_number

            # perform step size calculation and update bound information
            single_step!(d.exist_result, d.contractor_result, d.step_params,
                         d.step_result, d.method_f!, step_number-1, hj_limit,
                         delT, d.constant_state_bounds, d.polyhedral_constraint)

            # unpack storage
            if step_number - 1 > length(d.time)
                push!(d.storage, copy(d.contractor_result.X_computed))
                push!(d.storage_apriori, copy(d.exist_result.Xj_apriori))
                push!(d.time, d.contractor_result.times[1])
            else
                copyto!(d.storage[step_number], d.contractor_result.X_computed)
                copyto!(d.storage_apriori[step_number], d.exist_result.Xj_apriori)
                d.time[step_number] = d.step_result.time
                d.contractor_result.times[1] = d.step_result.time
            end
            if is_support_pnt
                stored_value_count += 1
                d.relax_t_dict_indx[stored_value_count] = step_number
                d.relax_t_dict_flt[next_support] = step_number
                support_indx += 1
                if support_indx <= length(tsupports)
                    next_support = tsupports[support_indx]
                else
                    next_support = Inf
                end
            end

            # throw error if limit exceeded
            if d.exist_result.status_flag !== RELAXATION_NOT_CALLED
                d.error_code = d.exist_result.status_flag
                break
            end
            d.step_count += 1
        else
            break
        end
    end

    if (d.step_count > d.step_limit) && (sign_tstep*d.time[d.step_count + 1] < sign_tstep*tmax)
        d.error_code = LIMIT_EXCEEDED
    end

    # cut out any unnecessary array elements
    resize!(d.storage, d.step_count + 1)
    resize!(d.storage_apriori, d.step_count + 1)
    resize!(d.time, d.step_count + 1)

    if d.error_code === RELAXATION_NOT_CALLED
        d.error_code = COMPLETED
    end

    return nothing
end


text:
module ScikitSpatial

using LinearAlgebra
using StaticArrays

export

    # Composite types (structs)
    Line,
    Plane,
    Circle,
    Sphere,

    # Measurement
    cosine_similarity,
    angle_between,
    distance,

    # Comparison
    is_zero,
    are_parallel,
    are_perpendicular,
    are_coplanar,
    on_surface,

    # Projection
    project,

    # Intersection
    intersect,

    # Transformation
    Vector,
    unit,
    to_point,
    centroid,
    mean_center

include("_util.jl")
include("types/base.jl")
include("computations/base.jl")

end

ext: jl


text:
# module Replanning
#
# using Parameters
# using LightGraphs
# using GraphUtils
# using CRCBS
# using ..TaskGraphs

using Printf: @sprintf

export
    ProjectRequest

"""
    ProjectRequest

Encodes a "request" that arrives in the factory command center.
    `t_request` encodes the time at which the request reaches the command center
    `t_arrival` is the time at which the project materials will be available

TODO more flexibility with t_arrival by allowing different materials to arrive
at different times
"""
@with_kw struct ProjectRequest
    schedule::OperatingSchedule = OperatingSchedule()
    t_request::Int              = -1
    t_arrival::Int              = -1
end
function project_request(args...)
    req = ProjectRequest(args...)
    set_t0!(req,req.t_request)
    process_schedule!(req.schedule)
    return req
end

export
    RepeatedAbstractPC_TAPF,
    RepeatedPC_TAPF,
    RepeatedC_PC_TAPF

abstract type RepeatedAbstractPC_TAPF end
CRCBS.get_env(prob::RepeatedAbstractPC_TAPF) = prob.env

"""
    RepeatedPC_TAPF{S<:SearchEnv} <: RepeatedAbstractPC_TAPF

Encodes a Repeated PC_TAPF problem.
Usage:
    `prob = RepeatedPC_TAPF(env::SearchEnv,requests::Vector{ProjectRequest})`
"""
struct RepeatedPC_TAPF{S<:SearchEnv} <: RepeatedAbstractPC_TAPF
    env::S
    requests::Vector{ProjectRequest}
end
RepeatedPC_TAPF(env::SearchEnv) = RepeatedPC_TAPF(env,Vector{ProjectRequest}())
struct RepeatedC_PC_TAPF{S<:SearchEnv} <: RepeatedAbstractPC_TAPF
    env::S
    requests::Vector{ProjectRequest}
end
RepeatedC_PC_TAPF(env::SearchEnv) = RepeatedC_PC_TAPF(env,Vector{ProjectRequest}())

for T in [:RepeatedPC_TAPF,:RepeatedC_PC_TAPF]
    @eval $T(prob::$T,env::SearchEnv) = $T(env,prob.requests)
end

construct_routing_problem(prob::RepeatedPC_TAPF,env) = PC_TAPF(env)
construct_routing_problem(prob::RepeatedC_PC_TAPF,env) = C_PC_TAPF(env)

"""
    trim_route_plan(search_env, route_plan, T)

Construct a trimmed route_plan that stops at a certain time step
"""
function trim_route_plan(search_env, route_plan, T)
    trim_solution!(build_env(search_env),route_plan,T)
end

export
    get_active_and_fixed_vtxs,
    get_active_vtxs,
    get_fixed_vtxs,
    split_active_vtxs!,
    fix_precutoff_nodes!,
    break_assignments!,
    prune_schedule,
    prune_project_schedule,
    splice_schedules!

"""
    get_active_and_fixed_vtxs(sched::OperatingSchedule,t)

"active" vertices "straddle" the query time t
"fixed" vertices finish before the query time t
"""
function get_active_and_fixed_vtxs(sched::OperatingSchedule,t)
    active_vtxs = Set{Int}()
    fixed_vtxs = Set{Int}()
    for v in vertices(get_graph(sched))
        if get_tF(sched,v) + minimum(get_local_slack(sched,v)) <= t
            push!(fixed_vtxs, v)
    elseif get_t0(sched,v) <= t < get_tF(sched,v) + minimum(get_local_slack(sched,v))
            push!(active_vtxs, v)
        end
    end
    @assert all(map(v->get_tF(sched,v) + minimum(get_local_slack(sched,v)), collect(active_vtxs)) .>= t)
    active_vtxs, fixed_vtxs
end
get_active_vtxs(sched,t) = get_active_and_fixed_vtxs(sched,t)[1]
get_fixed_vtxs(sched,t) = get_active_and_fixed_vtxs(sched,t)[2]

function split_action_node!(sched::OperatingSchedule,problem_spec::ProblemSpec,node::ScheduleNode,x,t)
    @assert isa(node.id,ActionID)
    pred1,pred2 = split_node(node.node,x)
    node1 = replace_in_schedule!(sched,pred1,node.id)
    set_t0!(node1,get_t0(node))
    set_tF!(node1,t)
    node2 = add_node!(sched,make_node(sched,problem_spec,pred1))
    set_t0!(node2,t)
    # set_tF!(sched,id2,get_tF(node))
    for v in outneighbors(sched,get_vtx(sched,node1.id))
        rem_edge!(sched,node1,v)
        add_edge!(sched,node2,v)
    end
    add_edge!(sched,node1,node2)
    node1,node2
end

"""
    split_active_vtxs!(sched::OperatingSchedule,
        problem_spec::ProblemSpec,t;

Split all GO nodes that "straddle" the cutoff time.
"""
function split_active_vtxs!(sched::OperatingSchedule,problem_spec::ProblemSpec,t;
    robot_positions::Dict{RobotID,ROBOT_AT}=Dict{RobotID,ROBOT_AT}()
    )
    G = get_graph(sched)
    # identify nodes "cut" by timestep t
    active_vtxs, fixed_vtxs = get_active_and_fixed_vtxs(sched,t)
    # split active nodes
    for v in active_vtxs
        node = get_node(sched,v)
        if isa(node,BOT_GO) # split TODO why aren't we splitting CARRY, COLLECT, and DEPOSIT as well?
            x = robot_positions[node.r].x
            node1,node2 = split_action_node!(sched,problem_spec,node,x,t)
            set_path_spec!(node1,PathSpec(get_path_spec(node1),fixed=true,plan_path=false))
            set_path_spec!(node2,PathSpec(get_path_spec(node1),tight=true))
        end
    end
    set_leaf_operation_vtxs!(sched)
    process_schedule!(sched)
    sched
end

"""
    fix_precutoff_nodes!(sched::OperatingSchedule,
        problem_spec::ProblemSpec,t)

Identify all nodes that end before the cutoff time, and change their path spec
    so that the route planner will not actually plan a path for them.
"""
function fix_precutoff_nodes!(sched::OperatingSchedule,problem_spec::ProblemSpec,t)
    # active_vtxs = Set{Int}()
    active_vtxs, fixed_vtxs = get_active_and_fixed_vtxs(sched,t)
    # set all fixed_vtxs to plan_path=false
    for v in fixed_vtxs
        set_path_spec!(sched,v,PathSpec(get_path_spec(sched,v), plan_path=false, fixed=true))
    end
    # verify that all vertices following active_vtxs have a start time > 0
    @assert all(map(v->get_t0(sched,v), collect(fixed_vtxs)) .<= t)
    @assert all(map(v->get_tF(sched,v) + minimum(get_local_slack(sched,v)), collect(active_vtxs)) .>= t)
    sched
end
function fix_precutoff_nodes!(env::SearchEnv,t=minimum(map(length, get_paths(get_route_plan(env)))))
    fix_precutoff_nodes!(get_schedule(env),get_problem_spec(env),t)
    return env
end

function break_assignments!(sched::OperatingSchedule,problem_spec,v)
    G = get_graph(sched)
    node_id = get_vtx_id(sched,v)
    node = get_node_from_id(sched,node_id)
    if isa(node, AbstractRobotAction)
        new_node = replace_robot_id(node,RobotID(-1)) # TODO Why is this line here? I don't think the robot id should be replaced in this node--just it's successors
        if isa(node,BOT_GO)
            for v2 in outneighbors(G,v)
                if isa(get_node_from_vtx(sched,v2),Union{BOT_COLLECT,TEAM_COLLECT})
                    rem_edge!(G,v,v2)
                    new_node = replace_destination(new_node,LocationID(-1))
                end
            end
        end
        replace_in_schedule!(sched,problem_spec,new_node,node_id)
    elseif isa(node,TEAM_ACTION)
        for i in 1:length(node.instructions)
            n = node.instructions[i]
            node.instructions[i] = replace_robot_id(node,RobotID(-1))
        end
    end
    return sched
end

"""
    break_assignments!(sched::OperatingSchedule,problem_spec::ProblemSpec)

Break all assignments that are eligible for replanning
"""
function break_assignments!(sched::OperatingSchedule,problem_spec::ProblemSpec)
    G = get_graph(sched)
    for v in vertices(G)
        path_spec = get_path_spec(sched,v)
        if path_spec.fixed == true
            continue
        end
        break_assignments!(sched,problem_spec,v)
    end
    propagate_valid_ids!(sched,problem_spec)
    sched
end

"""
    prune_schedule(sched::OperatingSchedule,
        problem_spec::ProblemSpec,t)

    remove nodes that don't need to be kept around any longer
"""
function prune_schedule(sched::OperatingSchedule,
        problem_spec::ProblemSpec,
        t,
        )
    G = get_graph(sched)

    # identify nodes "cut" by timestep
    active_vtxs, fixed_vtxs = get_active_and_fixed_vtxs(sched,t)
    # construct set of all nodes to prune out.
    remove_set = Set{Int}()
    for v in active_vtxs
        if isa(get_node_from_vtx(sched,v),Operation)
            push!(remove_set, v)
        end
        for e in edges(bfs_tree(G,v;dir=:in))
            if isa(get_node_from_vtx(sched, e.dst),Operation)
                push!(remove_set, e.dst)
            end
        end
    end
    for v in collect(remove_set)
        for e in edges(bfs_tree(G,v;dir=:in))
            if !(isa(get_node_from_vtx(sched, e.dst),ROBOT_AT))
                push!(remove_set, e.dst)
            end
        end
    end
    # Construct new graph
    new_sched = OperatingSchedule()
    keep_vtxs = setdiff(Set{Int}(collect(vertices(G))), remove_set)
    # add all non-deleted nodes to new project schedule
    for v in keep_vtxs
        add_node!(new_sched,get_node(sched,v))
    end
    # add all edges between nodes that still exist
    for e in edges(get_graph(sched))
        add_edge!(new_sched, get_vtx_id(sched, e.src), get_vtx_id(sched, e.dst))
    end
    # Initialize new cache
    # draw new ROBOT_AT -> GO edges where necessary
    for v in vertices(new_sched)
        node = get_node(new_sched, v)
        pred = node.node
        if isa(pred,BOT_GO) && indegree(new_sched,v) == 0
            robot_id = get_robot_id(pred)
            robot_node = replace_in_schedule!(new_sched, BOT_AT(robot_id, pred.x1), robot_id)
            set_t0!(robot_node,get_t0(node))
            add_edge!(new_sched, robot_node, node)
        end
        if isa(pred,Operation) && indegree(new_sched,v) < num_required_predecessors(pred)
            # add deleted objects back in
            input_ids = Set(map(v2->get_object_id(get_node_from_vtx(new_sched,v2)), inneighbors(new_sched,v)))
            for (object_id,o) in preconditions(pred)
                # if !(get_object_id(o) in input_ids)
                if !(object_id in input_ids)
                    o_node = add_node!(new_sched,make_node(sched,problem_spec,o))
                    add_edge!(new_sched, o_node.id, node)
                    set_t0!(o_node,get_t0(node))
                    set_tF!(o_node,get_t0(node))
                end
            end
        end
    end
    set_leaf_operation_vtxs!(new_sched)
    process_schedule!(new_sched)
    # init planning cache with the existing solution
    @assert sanity_check(new_sched)

    new_sched
end
prune_schedule(env::SearchEnv,args...) = prune_schedule(get_schedule(env),get_problem_spec(env),args...)
prune_schedule(solver,env::SearchEnv,args...) = prune_schedule(env,args...)

"""
    `prune_project_schedule`

Remove all vertices that have already been completed. The idea is to identify
all `Operation`s that are completed before `t`, remove all nodes upstream of
them (except for ROBOT_AT nodes), and create new edges between the ROBOT_AT
nodes and their first GO assignments.
"""
function prune_project_schedule(sched::OperatingSchedule,problem_spec::ProblemSpec,t;
        robot_positions::Dict{RobotID,ROBOT_AT}=Dict{RobotID,ROBOT_AT}()
    )
    new_sched = prune_schedule(sched,problem_spec,t)
    # split active nodes
    new_sched = split_active_vtxs!(new_sched,problem_spec,t;robot_positions=robot_positions)
    # freeze nodes that terminate before cutoff time
    fix_precutoff_nodes!(new_sched,problem_spec,t)
    # Remove all "assignments" from schedule
    break_assignments!(new_sched,problem_spec)
    new_sched
end

"""
    splice_schedules!(sched::P,next_sched::P) where {P<:OperatingSchedule}

Merge next_sched into sched
"""
function splice_schedules!(sched::P,next_sched::P,enforce_unique=true) where {P<:OperatingSchedule}
    for v in vertices(next_sched)
        node = get_node(next_sched,v)
        if !has_vertex(sched,node)
            add_node!(sched, node)
        elseif enforce_unique
            throw(ErrorException(string("Vertex ",v," = ",
                string(node.node),
                " already in sched."
                )))
        end
    end
    for e in edges(get_graph(next_sched))
        node1 = get_node(next_sched, e.src)
        node2 = get_node(next_sched, e.dst)
        add_edge!(sched, node1, node2)
    end
    set_leaf_operation_vtxs!(sched)
    sched
end
splice_schedules!(sched,next_sched) = sched

export
    ReplannerModel,
    ReplannerConfig,
    DeferUntilCompletion,
    ReassignFreeRobots,
    MergeAndBalance,
    ConstrainedMergeAndBalance,
    Oracle,
    NullReplanner,
    get_commit_time,
    replan!

function get_timeout_buffer end
function get_route_planning_buffer end
function get_commit_time end
function break_assignments! end
function set_time_limits! end
function split_active_vtxs! end
function fix_precutoff_nodes! end
function replan! end

"""
    ReplannerModel

Abstract type. Concrete subtypes currently include `DeferUntilCompletion`,
`ReassignFreeRobots`, `MergeAndBalance`, `Oracle`, `NullReplanner`
"""
abstract type ReplannerModel end

"""
    ReplannerConfig

Stores parameters for configuring an instance of `ReplannerModel`.
Fields:
* real_time::Bool - if `true`, adjust solver runtime limits to meet real time
    operation constraints.
* timeout_buffer::Float64 - defines the minimum amount of time that must be for
    replanning the route
* route_planning_buffer::Float64 - defines the minimum amount of time that must
    be allocated for replanning the route
* commit_threshold::Float64 - defines where the ``freeze time'', or the time
    step from which the solver may change the existing plan.
"""
@with_kw mutable struct ReplannerConfig
    real_time::Bool                 = true
    timeout_buffer::Float64         = 1
    route_planning_buffer::Float64  = 2
    commit_threshold::Int           = 10
    max_time_limit::Int             = 100
end
get_replanner_config(config::ReplannerConfig) = config
get_replanner_config(model) = model.config

export
    get_real_time_flag,
    get_timeout_buffer,
    get_route_planning_buffer,
    get_commit_threshold,
    get_max_time_limit,
    set_real_time_flag!,
    set_timeout_buffer!,
    set_route_planning_buffer!,
    set_commit_threshold!,
    set_max_time_limit!

get_real_time_flag(model)           = get_replanner_config(model).real_time
get_timeout_buffer(model)           = get_replanner_config(model).timeout_buffer
get_route_planning_buffer(model)    = get_replanner_config(model).route_planning_buffer
get_commit_threshold(model)         = get_replanner_config(model).commit_threshold
get_max_time_limit(model)           = get_replanner_config(model).max_time_limit
function set_real_time_flag!(model,val) get_replanner_config(model).real_time = val end
function set_timeout_buffer!(model,val) get_replanner_config(model).timeout_buffer = val end
function set_route_planning_buffer!(model,val) get_replanner_config(model).route_planning_buffer = val end
function set_commit_threshold!(model,val) get_replanner_config(model).commit_threshold = val end
function set_max_time_limit!(model,val) get_replanner_config(model).max_time_limit = val end

"""
    DeferUntilCompletion <: ReplannerModel

Allow work to begin on the new project only after all other work is completed.
"""
@with_kw struct DeferUntilCompletion <: ReplannerModel
    config::ReplannerConfig = ReplannerConfig()
end
"""
    ReassignFreeRobots   <: ReplannerModel

Allow robots to begin working on the new project as soon as they have finished
their current assignments.
"""
@with_kw struct ReassignFreeRobots   <: ReplannerModel
    config::ReplannerConfig = ReplannerConfig()
end
"""
    MergeAndBalance      <: ReplannerModel

Allow replanning from scratch for all assignments and routes except for those
that will take place
"""
@with_kw struct MergeAndBalance      <: ReplannerModel
    config::ReplannerConfig = ReplannerConfig()
end
"""
    ConstrainedMergeAndBalance <: ReplannerModel

Identical to `MergeAndBalance`, except that the replanning commit time may be
pushed farther into the future to limit the total problem size at each replanning
stage.
"""
@with_kw struct ConstrainedMergeAndBalance      <: ReplannerModel
    config::ReplannerConfig = ReplannerConfig()
    max_problem_size::Int   = 400 # max permissible number of variables in sub problem
end
@with_kw struct Oracle               <: ReplannerModel
    config::ReplannerConfig = ReplannerConfig(
        timeout_buffer         = -110,
        route_planning_buffer   = 10,
        commit_threshold        = 0,
    )
end
@with_kw struct NullReplanner        <: ReplannerModel
    config::ReplannerConfig = ReplannerConfig()
end

export
    ReplanningProfilerCache,
    FullReplanner

"""
    ReplanningProfilerCache

Stores information during replanning
* `schedule` - maintains (non-pruned) overall project schedule
* `stage_ids` - maps id of each node in schedule to stage index
* `stage_results` - vector of result dicts
* `features` - vector of `::FeatureExtractors` to extract after each stage
* `final_features` - vector of `::FeatureExtractors` to extract after final stage
"""
@with_kw struct ReplanningProfilerCache
    schedules::Vector{OperatingSchedule}    = Vector{OperatingSchedule}()
    project_ids::Dict{AbstractID,Int}       = Dict{AbstractID,Int}()
    stage_results::Vector{Dict{String,Any}} = Vector{Dict{String,Any}}()
    final_results::Dict{String,Any}         = Dict{String,Any}()
    features::Vector{FeatureExtractor}      = Vector{FeatureExtractor}()
    final_features::Vector{FeatureExtractor}= Vector{FeatureExtractor}()
end
function reset_cache!(cache::ReplanningProfilerCache)
    empty!(cache.schedules)
    empty!(cache.project_ids)
    empty!(cache.stage_results)
    empty!(cache.final_results)
    cache
end

"""
    FullReplanner{R,S}

Wrapper for a `ReplannerModel` and a `PC_TAPF` solver, as well as a
`ReplanningProfilerCache`.
"""
@with_kw struct FullReplanner{R,S}
    solver::S                       = NBSSolver()
    replanner::R                    = MergeAndBalance()
    cache::ReplanningProfilerCache  = ReplanningProfilerCache()
end
reset_cache!(planner::FullReplanner) = reset_cache!(planner.cache)
get_replanner_config(planner::FullReplanner) = get_replanner_config(planner.replanner)
for op in [
        :get_commit_time,
        :break_assignments!,
        :set_time_limits!,
        :split_active_vtxs!,
        :fix_precutoff_nodes!,
        ]
    @eval $op(planner::FullReplanner,args...) = $op(planner.replanner,args...)
end
CRCBS.get_logger(planner::FullReplanner) = CRCBS.get_logger(planner.solver)

export ReplannerWithBackup

"""
    ReplannerWithBackup{A,B}

Consists of a primary and a backup solver. The backup solver is supposed to be
fast--it computes a "fallback" plan in case the primary solver fails to find a
good enough solution before time out.
"""
struct ReplannerWithBackup{A,B}
    primary_planner::A
    backup_planner::B
end
function set_real_time_flag!(model::ReplannerWithBackup,val)
    set_real_time_flag!(model.backup_planner,val)
    set_real_time_flag!(model.primary_planner,val)
end
function reset_cache!(planner::ReplannerWithBackup)
    reset_cache!(planner.primary_planner)
    reset_cache!(planner.backup_planner)
end

function get_commit_time(replan_model, search_env, request, commit_threshold=get_commit_threshold(replan_model))
    request.t_request + commit_threshold
end
get_commit_time(replan_model::Oracle, search_env, request, args...) = request.t_request
get_commit_time(replan_model::DeferUntilCompletion, search_env, request, commit_threshold) = max(request.t_request + commit_threshold, maximum(get_tF(get_schedule(search_env))))
get_commit_time(replan_model::NullReplanner,args...) = get_commit_time(DeferUntilCompletion(),args...)
function get_commit_time(replan_model::ReassignFreeRobots, search_env, request, commit_threshold)
    free_time = makespan(get_schedule(search_env))
    for v in vertices(get_schedule(search_env))
        # node = get_node_from_vtx(get_schedule(search_env),v)
        node = get_node(get_schedule(search_env),v)
        if isa(node.node,GO)
            if get_id(get_destination_location_id(node.node)) == -1
                free_time = min(free_time, get_t0(node))
            end
        end
    end
    max(request.t_request + commit_threshold,free_time)
end
function get_commit_time(replan_model::ConstrainedMergeAndBalance, search_env, request, commit_threshold)
    nv_max = replan_model.max_problem_size - nv(request.schedule)
    t_commit = request.t_request + commit_threshold
    if 0 < nv_max < nv(get_schedule(search_env))
        t_commit = max(t_commit, sort(get_tF(get_schedule(search_env)); rev=true)[nv_max])
    end
    return t_commit
end

break_assignments!(replan_model::ReplannerModel,args...) = break_assignments!(args...)
break_assignments!(replan_model::ReassignFreeRobots,args...) = nothing
break_assignments!(replan_model::DeferUntilCompletion,args...) = nothing
break_assignments!(replan_model::NullReplanner,args...) = nothing

function set_time_limits!(replan_model,solver,t_request,t_commit)
    real_time = get_real_time_flag(replan_model)
    if real_time
        set_time_limits!(real_time,replan_model,solver,t_request,t_commit)
    end
end

function set_time_limits!(flag::Bool,replan_model,solver,t_request,t_commit)
    # set_runtime_limit!(solver, (t_commit - t_request) - get_timeout_buffer(replan_model))
    # # set_runtime_limit!(assignment_solver(solver), solver.time_limit - get_route_planning_buffer(replan_model))
    # set_runtime_limit!(solver.assignment_model, runtime_limit(solver) - get_route_planning_buffer(replan_model))
    # @assert runtime_limit(solver) > 0.0
    # solver
    set_runtime_limit!(solver, (t_commit - t_request) - get_timeout_buffer(replan_model))
    set_runtime_limit!(solver, min(runtime_limit(solver),get_max_time_limit(replan_model)))
    set_runtime_limit!(assignment_solver(solver), runtime_limit(solver) - get_route_planning_buffer(replan_model))
    @assert runtime_limit(solver) > 0.0
    set_deadline!(solver, time() + (t_commit - t_request))
    solver
end
# function set_time_limits!(flag::Bool,replan_model::DeferUntilCompletion,solver,t_request,t_commit)
#     set_runtime_limit!(solver, (t_commit - t_request) - get_timeout_buffer(replan_model))
#     set_runtime_limit!(solver, min(runtime_limit(solver),get_max_time_limit(replan_model)))
#     set_runtime_limit!(assignment_solver(solver), runtime_limit(solver) - get_route_planning_buffer(replan_model))
#     @assert runtime_limit(solver) > 0.0
#     solver
# end
function set_time_limits!(flag::Bool,replan_model::NullReplanner,solver,t_request,t_commit)
    set_runtime_limit!(solver,-1.0)
    return solver
end

split_active_vtxs!(replan_model::ReplannerModel,new_sched,problem_spec,t_commit;kwargs...) = split_active_vtxs!(new_sched,problem_spec,t_commit;kwargs...)

fix_precutoff_nodes!(replan_model,new_sched,problem_spec,t_commit) = fix_precutoff_nodes!(new_sched,problem_spec,t_commit)

"""
    replan!(solver, replan_model, search_env, env_graph, problem_spec, route_plan,
        next_sched, t_request, t_arrival; commit_threshold=5,kwargs...)

Combine an existing solution with a new project request `next_sched` that
is requested at time `t_request`, and whose raw materials become available at
`t_arrival`.
The exact behavior of this function depends on `replan_model <: ReplannerModel`
"""
# function replan!(solver, replan_model, search_env, env_graph, problem_spec, route_plan, next_sched, t_request, t_arrival; commit_threshold=5,kwargs...)
function replan!(solver, replan_model, search_env, request;
        commit_threshold=get_commit_threshold(replan_model),
        kwargs...
        )
    sched               = get_schedule(search_env)
    cache               = get_cache(search_env)
    route_plan          = get_route_plan(search_env)
    problem_spec        = get_problem_spec(search_env)
    env_graph           = get_graph(search_env)
    t_request           = request.t_request
    next_sched          = request.schedule

    @assert sanity_check(sched," in replan!()")
    # Freeze route_plan and schedule at t_commit
    t_commit = get_commit_time(replan_model, search_env, request, commit_threshold)
    t_final = minimum(map(length, get_paths(get_route_plan(search_env))))
    t_split = min(t_commit,t_final)
    robot_positions=get_env_snapshot(search_env,t_split)
    reset_solver!(solver)
    set_time_limits!(replan_model,solver,t_request,t_commit)
    # Update operating schedule
    new_sched = prune_schedule(replan_model,search_env,t_split)
    @assert sanity_check(new_sched," after prune_schedule()")
    # split active nodes
    new_sched = split_active_vtxs!(replan_model,new_sched,problem_spec,t_split;robot_positions=robot_positions)
    @assert sanity_check(new_sched," after split_active_vtxs!()")
    # freeze nodes that terminate before cutoff time
    new_sched = fix_precutoff_nodes!(replan_model,new_sched,problem_spec,t_split)
    # Remove all "assignments" from schedule
    break_assignments!(replan_model,new_sched,problem_spec)
    @assert sanity_check(new_sched," after break_assignments!()")
    # splice projects together!
    set_t0!(next_sched,t_commit)
    splice_schedules!(new_sched,next_sched)
    process_schedule!(new_sched)
    @assert sanity_check(new_sched," after splice_schedules!()")
    base_search_env = construct_search_env(
        solver,
        new_sched,
        search_env,
        initialize_planning_cache(new_sched)
        )
    base_route_plan = initialize_route_plan(search_env,get_cost_model(base_search_env))
    trimmed_route_plan = trim_route_plan(base_search_env, base_route_plan, t_commit)
    SearchEnv(base_search_env, route_plan=trimmed_route_plan)
end
# replan!(solver, replan_model::NullReplanner, search_env, args...;kwargs...) = search_env
function replan!(solver::FullReplanner,args...;kwargs...)
    replan!(solver.solver,solver.replanner,args...;kwargs...)
end

# end


text:
module MPCC2DPlot

using Relaxation
using PyPlot

"""
Plot les contours de la relaxation puis les itérées pour une contrainte de complémentarité
z : tableau de dimension 2
"""
function RelaxationPlot(z::Any,r::Float64,s::Float64,t::Float64,a::Float64,b::Float64)

 #Step 1 : plot des contraintes de relaxations
 x = linspace(-r,b,1000); y = Relaxation.psi(x,r,s,t)
 PyPlot.plot(x, y, color="red", linewidth=2.0, linestyle="-")
 PyPlot.plot(y, x, color="blue", linewidth=2.0, linestyle="-")

 #Step 2 : plot des contraintes de positivité
 x = linspace(-r,b,1000); y = -r*ones(length(x))
 PyPlot.plot(x, y, color="green", linewidth=2.0, linestyle="-")
 PyPlot.plot(y, x, color="green", linewidth=2.0, linestyle="-")

 #Step 3: tracer les itérés
 PyPlot.plot(z[:,1], z[:,2],marker="o",markerfacecolor="black", linestyle="")
 return true
end

#end of module
end


text:
struct GosperIterator
    n::Int64
    k::Int64
end

"""
    gosper(n, k)

Returns iterator over numbers whose binary representation correspond to the
indicator vector of a subset of {1,...,n} of size k.

Example:
   collect(gosper(3,2)) == [3, 5, 6]
"""
function gosper(n::Int, k::Int)
    @assert 0 < k <= n <= 62
    return GosperIterator(n, k)
end

const GosperState = Tuple{Int64, Int64}

eltype(it::GosperIterator) = Int64
length(it::GosperIterator) = binomial(it.n, it.k)

function start(it::GosperIterator)
    (2^(it.k) - 1, 2^(it.n) -1)
end

function iterate(it::GosperIterator, state::GosperState = start(it))

    if state[1] > state[2]
        return nothing
    end

        # Gosper's hack: Finds the next smallest number with exactly k bits
    # set to 1 in its binary representation.
    o = state[1]
    u = state[1] & -state[1]
    v = u + state[1]
    y = v + (div(xor(v, state[1]), u) >> 2)
    return o, (y, state[2])
end

struct BinaryIndexIterator
    d::Int64#  dimension
    indices::Array{Int64}# indices of 1's
    init_value::Int64# value that the iterator add the subnumbers to it
end

"""
For some a in {0,1}^d, that have 1's in the places "indices", iterate over all substes
of indices in "indices" and return the number of it's binary representation
"""
function BinaryItr(d::Int64, indices::Array{Int64}, init_value::Int64 = 0)
    @assert 0 < d <= 62

    BinaryIndexIterator(d, indices, init_value)
end

"""
Now the indices of the 1's are the number itself with those indices
"""
function BinaryItr(d::Int64, indices::Int64, init_value::Int64 = 0)
    @assert 0 < d <= 62

    _, arr_idx = norm_1_with_indices(indices)

    BinaryIndexIterator(d, Array{Int64}(arr_idx), init_value)
end

const BinaryIndexState = Int64

eltype(it::BinaryIndexIterator) = Int64
length(it::BinaryIndexIterator) = 2^(length(it.indices))

function start(it::BinaryIndexIterator)
   0
end

function iterate(it::BinaryIndexIterator, state::BinaryIndexState = start(it))

    if state == length(it)
        return nothing
    end

    sub_binary_number = it.init_value

    temp_state = state

    for i=1:length(it.indices)
        if temp_state == 0
            break
        end
        if temp_state%2 == 1
            sub_binary_number += 1<<(it.indices[i]-1)
        end
        temp_state >>= 1
    end

    state += 1

    return sub_binary_number, state

end

"""
Calculate the norm1 of binary number (the number of 1's in the number).
return: -norm1(alpha)
        -array with the indices of the 1's
"""
function norm_1_with_indices(alpha::Int64)

    @assert alpha >= 0

    count = 0
    index = 1
    idx = []
    while alpha > 0
        if alpha % 2 == 1
            count += 1
            push!(idx, index)
        end
        alpha >>= 1
        index += 1
    end
    count, Array{Int64}(idx)
end

"""
Calculate the norm1 of binary number (the number of 1's in the number).
return: -norm1(alpha)
        -array with the indices of the 1's
"""
function norm_1(alpha::Int64)

    @assert alpha >= 0

    count = 0
    while alpha > 0
        if alpha % 2 == 1
            count += 1
        end
        alpha >>= 1
    end
    count
end


text:
# This script is licensed under the MIT License (MIT).

jupyter_cmds = ["jupyter"]
python_cmds = ["python3", "python"]

function detect_application(
    application_name::AbstractString,
    envvar_name::AbstractString,
    julia_pkg_name::AbstractString,
    python_pkg_name::AbstractString,
    command_candidates::AbstractVector{<:AbstractString}
)
    if get(ENV, envvar_name, nothing) == ""
        @info "\$$envvar_name is set to \"$(ENV[envvar_name])\", $julia_pkg_name will create a Julia-internal Conda installation for $python_pkg_name."
    else
        if haskey(ENV, envvar_name)
            cmd_name = ENV[envvar_name]
        else
            cmd_idxs = findall(!isequal(nothing), Sys.which.(command_candidates))
            if !isempty(cmd_idxs)
                cmd_name = command_candidates[first(cmd_idxs)]
            else
                @error "No $(join(command_candidates, " or ")) command available, consider installing Anaconda and adding the Anaconda \"bin\" directory to your \$PATH"
                exit(1)
            end
        end

        try
            cmd_version = chomp(read(`$cmd_name --version`, String))
            # ENV[envvar_name] = cmd_name
            @info "Using command \"$cmd_name\" for $application_name, $application_name version $cmd_version"
        catch err
            @error "Trying to use $application_name command \"$cmd_name\", but \"$cmd_name --version\" fails."
            exit(1)
        end

        if application_name == "Python"
            try
                read(`$cmd_name -c "import matplotlib.pyplot"`, String)
                @info "Python can import matplotlib.pyplot"
            catch err
                @error "Python cannot import matplotlib.pyplot, consider installing Anaconda and adding the Anaconda \"bin\" directory to your \$PATH"
                exit(1)
            end
        end

        if !haskey(ENV, envvar_name)
            println("Please set environment variable $envvar_name=\"$cmd_name\"")
        end
    end
end

detect_application("Jupyter", "JUPYTER", "PyCall.jl/PyPlot.jl", "Jupyter", jupyter_cmds)
detect_application("Python", "PYTHON", "PyCall.jl/PyPlot.jl", "matplotlib/PyPlot", python_cmds)


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ a70cdae2-14c4-11eb-0f58-3bf672d63847
begin
	using SpecialFunctions
	using SpecialPolynomials
	using HypergeometricFunctions
	using PyPlot
	using Polynomials
	using PyCall
end

# ╔═╡ 814bda9c-14c4-11eb-0ba8-a71a2d9c234b
md"### Playing with the discrete polynomials definitions"

# ╔═╡ bb8453e6-18b7-11eb-22dd-519ac2f2d0f6
mp = pyimport("mpmath")

# ╔═╡ e4232312-14c6-11eb-2578-d775876aa9f4
begin
	rc("text", usetex = true)
	ion()
end

# ╔═╡ 65817e78-14c6-11eb-2016-9d7cf09a2801
begin
	j = 16
	N = 2*j
	const p = 2
end

# ╔═╡ c6c1b848-18b7-11eb-03ab-ef48179d3b5d
function kravp_mp(n,q)
	convert(Float64,mp.hyp2f1(-n,-q,-N,p,zeroprec = 250))
end

# ╔═╡ 8f1e8a68-18b9-11eb-2888-7d8b6effde6a
function nbinom(n,k)
    1/((n+1)*beta(k+1, n-k+1))
end

# ╔═╡ faa722e2-18b7-11eb-1899-e5bdcf897db3
function kravf_mp(n,q)
	((-1)^n/(2^j))*sqrt(nbinom(2*j,n)*nbinom(2*j,j+q))*kravp_mp(n,q+j)
end

# ╔═╡ 11d8fc38-18b8-11eb-19fe-c739310c212c
begin
	K1v_mp = zeros(Float64,(N+1,N+1))
	for l in -j:j
		for m in 0:N
			K1v_mp[l+j+1,m+1] = kravf_mp(m,l)
		end
	end
end

# ╔═╡ d6f280d4-18b8-11eb-3240-1fb563728533
imshow(K1v_mp,cmap="gray")

# ╔═╡ 75b2d10e-18bb-11eb-0f1e-0fc1d7d51b07
begin
	Kn = zeros(Float64,(N+1,N+1))
	for l in -j:j
		for m in 0:N
			Kn[m+1,l+j+1] = kravf_mp(m,l)
		end
	end
end

# ╔═╡ a4fcc568-18bc-11eb-0f1d-a172c97d48d2
begin
	imshow(Kn,cmap="gray")
	xlabel("q")
	xticks([0,j,N],[-j,0,j])
	ylabel("n")
	yticks([0,j,N],[N,j,0])
end

# ╔═╡ f1ed4446-18bd-11eb-27b3-6da15e29c6d5
plot(rotl90(Kn)[N+1,:])

# ╔═╡ 65e59e8e-18be-11eb-3e41-cb3799c1cc46
plot(Kn[N+1,:])

# ╔═╡ Cell order:
# ╟─814bda9c-14c4-11eb-0ba8-a71a2d9c234b
# ╠═a70cdae2-14c4-11eb-0f58-3bf672d63847
# ╟─bb8453e6-18b7-11eb-22dd-519ac2f2d0f6
# ╠═e4232312-14c6-11eb-2578-d775876aa9f4
# ╠═65817e78-14c6-11eb-2016-9d7cf09a2801
# ╠═c6c1b848-18b7-11eb-03ab-ef48179d3b5d
# ╠═8f1e8a68-18b9-11eb-2888-7d8b6effde6a
# ╠═faa722e2-18b7-11eb-1899-e5bdcf897db3
# ╠═11d8fc38-18b8-11eb-19fe-c739310c212c
# ╠═d6f280d4-18b8-11eb-3240-1fb563728533
# ╠═75b2d10e-18bb-11eb-0f1e-0fc1d7d51b07
# ╠═a4fcc568-18bc-11eb-0f1d-a172c97d48d2
# ╠═f1ed4446-18bd-11eb-27b3-6da15e29c6d5
# ╠═65e59e8e-18be-11eb-3e41-cb3799c1cc46


text:
## roctypes.jl  Types used for ROC computations
## (c) 2013--2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

## TNT: a target-nontarget tuple
"""
`TNT(tar, non`) is a type that holds target and non-target scores in a test for a two-class
classifier.
"""
type TNT{T<:Real}
    """target scores (scofres for which, actually, the one hypothesis is true)"""
    tar::Vector{T}              # target scores
    """non-target scores (score for which, actually, the oter hypothesis is true)"""
    non::Vector{T}              # non-target scores
end

## Roc: essential data to store Receiver Operating Characteristics
"""
`Roc` is a type that stores the essential performance information that can be extracted from a
set of supervised trials, i.e., target and non-target scores from a two-class classifier.  Apart from
the (minimalized) arrays for probabilities of false-alarm and miss---the coordinates of the ROC
curve---, they are the threshold for these, a boolean whether-or-not this point lies on the
convex hull, and the associated optimal log-likelihood-ratio associated to the line segment.
"""
type Roc{T<:Real}
    pfa::Vector{Float64}        # probability of false alarm
    pmiss::Vector{Float64}      # probability of miss
    ch::BitArray                # point lies on convex hull
    θ::Vector{T}                # threshold
    llr::Vector{Float64}        # optimal log likelihood ratio
end

## A traditional decision cost function
"""
`DCF` is a type that holds one or more cost functions with which the performance of a two-class
classifier can be assessed.
"""
type DCF{PTT,CFT,CMT}
    ptar::PTT               # target prior
    cfa::CFT                # cost of a false alarm
    cmiss::CMT              # cost of a miss
    function DCF(ptar, cfa, cmiss)
        lengths = map(length, (ptar, cfa, cmiss))
        lmax = maximum(lengths)
        for l in lengths
            l in [1,lmax] || error("Inconsistent vector lengths")
        end
        new(ptar, cfa, cmiss)
    end
end
if VERSION < v"0.4.0-dev"
    typealias ArrayOrReal{T<:Real} Union(Array{T}, Real)
else
    typealias ArrayOrReal{T<:Real} Union{Array{T}, Real}
end
DCF(ptar::ArrayOrReal, cfa::ArrayOrReal, cmiss::ArrayOrReal) = DCF{typeof(ptar),typeof(cfa),typeof(cmiss)}(ptar, cfa, cmiss)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(sdata::AbstractData; variables=nothing)
  # retrieve underlying domain
  sdomain = domain(sdata)

  # valid variables
  validvars = sort([var for (var, V) in GeoStatsBase.variables(sdata)])

  # plot all variables by default
  variables == nothing && (variables = validvars)
  @assert variables ⊆ validvars "invalid variable name"

  # shared plot specs
  layout --> length(variables)

  for (i, var) in enumerate(variables)
    # retrieve valid values
    vals = map(1:npoints(sdata)) do ind
      if isvalid(sdata, ind, var)
        v = sdata[ind,var]
        v isa Number ? v : get(v)
      else
        NaN
      end
    end
    @series begin
      subplot := i
      title --> string(var)
      legend --> false
      sdomain, vals
    end
  end
end


text:
#!/bin/env julia

const filename = "julia-vim-L2U-table"
const sourcefile = joinpath(JULIA_HOME, "..", "..", "doc", "manual", "unicode-input-table.rst")

isreadable(sourcefile) || error("unable to find source file $sourcefile")

function fix_compose_chars(uni)
  u = uni[1]
  isc = ('\u0300' ≤ u ≤ '\u036F') ||
        ('\u1DC0' ≤ u ≤ '\u1DFF') ||
        ('\u20D0' ≤ u ≤ '\u20FF') ||
        ('\uFE20' ≤ u ≤ '\uFE2F')
  return isc ? "\u25CC" * uni : uni
end

open("$filename.txt","w") do f
    print(f, """
      $filename.txt  LaTeX-to-Unicode reference table

      ===================================================================
      LATEX-TO-UNICODE REFERENCE TABLE    *L2U-ref* *julia-vim-L2U-reference*

        Note: This file is autogenerated from the script '$(basename(Base.source_path()))'
        The symbols are based on the documentation of Julia version $VERSION
        See |julia-vim| for the LaTeX-to-Unicode manual.

      """)

    col_headers = ["Code point(s)", "Character(s)", "Tab completion sequence(s)", "Unicode name(s)"]

    code = String[]
    unicode = String[]
    latex = String[]
    desc = String[]

    open(sourcefile) do sf
        for l in eachline(sf)
            m = match(r"
                ^(U\+[0-9A-F]{5}(\ \+\ U\+[0-9A-F]{5})*)
                \s+
                (\X)
                \s+
                (\\\\[^[:space:],]+(,\ \\\\[^[:space:],]+)*)
                \s+
                (.*)$
                "x, l)
            m == nothing && continue
            c, _, u, l, _, d = m.captures
            push!(code, c)
            push!(unicode, fix_compose_chars(u))
            push!(latex, replace(l, "\\\\", "\\"))
            push!(desc, d)
        end
    end

    cw = max(length(col_headers[1]), maximum(map(length, code)))
    uw = max(length(col_headers[2]), maximum(map(length, unicode)))
    lw = max(length(col_headers[3]), maximum(map(length, latex)))
    dw = max(length(col_headers[4]), maximum(map(length, desc)))

    print_padded(c, u, l, d) = println(f, rpad(c, cw), " ", rpad(u, uw), " ", rpad(l, lw), " ", d)

    print_padded(col_headers[1:3]..., col_headers[4] * "~")
    print_padded("-"^cw, "-"^uw, "-"^lw, "-"^dw)

    for (c, u, l, d) in zip(code, unicode, latex, desc)
        print_padded(c, u, l, d)
    end
    print_padded("-"^cw, "-"^uw, "-"^lw, "-"^dw)

    println(f)
    println(f, "  vim", # separated to avoid an error from vim which otherwise tries to parse this line
               ":tw=$(cw+uw+lw+dw+3):et:ft=help:norl:")
end

max_stars_repo_name: mzdravkov/dotfiles

max_issues_repo_name: mzdravkov/dotfiles

max_forks_repo_name: mzdravkov/dotfiles


text:
"""
Copyright (c) 2014: Brian J Smith and other contributors:

https://github.com/brian-j-smith/Mamba.jl/contributors

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""

########################### Demo #########################
# x, d = 5, Gamma(1,1)
# transformedx = link(d, x)
# invlink(d, transformedx)
# invlink(d, -1)

#################### Fallback Methods ####################

link(d::Distribution, x, transform::Bool=true) = x
invlink(d::Distribution, x, transform::Bool=true) = x

function logpdf(d::UnivariateDistribution, x, transform::Bool)
  all(insupport(d, x)) ? logpdf(d, x) : -Inf
end

function logpdf(d::MultivariateDistribution, x, transform::Bool)
  all(insupport(d, x)) ? logpdf(d, x) : -Inf
end

function logpdf(d::MatrixDistribution, x, transform::Bool)
  all(insupport(d, x)) ? logpdf(d, x) : -Inf
end

#################### Discrete Support Grids ####################

typealias GridUnivariateDistribution
          Union{Bernoulli, Binomial, Categorical, DiscreteUniform,
                Hypergeometric, NoncentralHypergeometric}

grid(d::GridUnivariateDistribution) = Float64[minimum(d):maximum(d)]
grid(d::Distribution) =
  error("discrete grid not available for ", typeof(d), " distributions")

#################### PDMatDistribution ####################

typealias PDMatDistribution Union{InverseWishart, Wishart}

function link(D::PDMatDistribution, X::Matrix, transform::Bool=true)
  n = dim(D)
  value = similar(X, int(n * (n + 1) / 2))
  k = 1
  if transform
    U = chol(X)
    for i in 1:n
      value[k] = log(U[i,i])
      k += 1
    end
    for i in 1:n, j in (i+1):n
      value[k] = U[i,j]
      k += 1
    end
  else
    for i in 1:n, j in i:n
      value[k] = X[i,j]
      k += 1
    end
  end
  value
end

function invlink(D::PDMatDistribution, x::Vector, transform::Bool=true)
  n = dim(D)
  value = zeros(VariateType, n, n)
  k = 1
  if transform
    for i in 1:n
      value[i,i] = exp(x[k])
      k += 1
    end
    for i in 1:n, j in (i+1):n
      value[i,j] = x[k]
      k += 1
    end
    return At_mul_B(value, value)
  else
    for i in 1:n, j in i:n
      value[i,j] = value[j,i] = x[k]
      k += 1
    end
    return value
  end
end

function logpdf{T<:Real}(D::PDMatDistribution, X::Matrix{T}, transform::Bool)
  value = logpdf(D, X)
  if transform && isfinite(value)
    U = chol(X)
    n = dim(D)
    for i in 1:n
      u = U[i,i]
      value += log(2 * u) + (n - i + 1) * log(u)
    end
  end
  value
end

#################### TransformDistribution ####################

typealias TransformDistribution{T<:ContinuousUnivariateDistribution}
  Union{T, Truncated{T}}

function link(d::TransformDistribution, x, transform::Bool=true)
  if transform
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
      return logit((x - a) / (b - a))
    elseif lowerbounded
      return log(x - a)
    elseif upperbounded
      return log(b - x)
    else
      return x
    end
  else
     return x
  end
end

function invlink(d::TransformDistribution, x, transform::Bool=true)
  if transform
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
      return (b - a) * invlogit(x) + a
    elseif lowerbounded
      return exp(x) + a
    elseif upperbounded
      return b - exp(x)
    else
      return x
    end
  else
    return x
  end
end

function logpdf(d::TransformDistribution, x::Real, transform::Bool)
  insupport(d, x) || return -Inf
  value = logpdf(d, x)
  if transform
    a, b = minimum(d), maximum(d)
    lowerbounded, upperbounded = isfinite(a), isfinite(b)
    if lowerbounded && upperbounded
      y = (x - a) / (b - x)
      value += log((b - a) * y / (y + 1.0)^2)
    elseif lowerbounded
      value += log(x - a)
    elseif upperbounded
      value += log(b - x)
    end
  end
  value
end

# Special case for NormalInverseGamma
function link(d :: NormalInverseGamma, x :: Tuple{Real,Real})
  (mean,var) = x
  return (mean, log(var))
end

function invlink(d :: NormalInverseGamma, x :: Tuple{Real,Real})
  (mean,logvar) = x
  return (mean, exp(logvar))
end

function logpdf(d :: NormalInverseGamma, x :: Tuple{Real,Real}, t :: Bool)
  (mean,var) = x
  var < 0 ? -Inf : logpdf(d,mean,var)
end


text:
AgeDistribution = zeros(NBYears, 52)#[zeros(52) for i in 1:NBYears] # 52 = 18 -> 70

# Compute Age Distribution for Initial Manpower
TempActualVar = 1
for imp in 1:length(InitManpower)
    for i in 1:length(InitMPPartCP[imp])
        for j in (1 + InitManpower[imp].Seniority):GuyCareerPaths[InitMPPartCP[imp][i]].Length
            AgeDistribution[j - InitManpower[imp].Seniority, GuyCareerPaths[InitMPPartCP[imp][i]].RecruitmentAge + j + InitManpower[imp].ActualSeniority - InitManpower[imp].Seniority - 18] += sol.sol[TempActualVar]
        end
        TempActualVar += 1
    end
end

# New recruits age distribution
for y in 1:NBYears
    for i in 1:length(GuyCareerPaths)
        for j in 1:GuyCareerPaths[i].Length # j => year for const
            if (y + j - 1) <= NBYears
                AgeDistribution[y + j - 1, GuyCareerPaths[i].RecruitmentAge + j - 18] +=  (1 - AttritionAtYear(GuyCareerPaths[i], j))* sol.sol[ InitMPDivisionNb + (y - 1)*length(GuyCareerPaths) + i]
            end
        end
    end
end

using Plots
plotly()

gui( surface([18:70],  [SimulationYear:(SimulationYear+NBYears)], AgeDistribution, xlabel = "Age (y)", ylabel = "Year", zlabel = "Amount", size=(1200,900), title = "Age distribution", color = :lightrainbow))

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Transformer2W <: ACBranch
        name::String
        available::Bool
        activepower_flow::Float64
        reactivepower_flow::Float64
        arc::Arc
        r::Float64
        x::Float64
        primaryshunt::Float64
        rate::Union{Nothing, Float64}
        services::Vector{Service}
        ext::Dict{String, Any}
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

The 2-W transformer model uses an equivalent circuit assuming the impedance is on the High Voltage Side of the transformer. The model allocates the iron losses and magnetizing susceptance to the primary side.

# Arguments
- `name::String`
- `available::Bool`
- `activepower_flow::Float64`
- `reactivepower_flow::Float64`
- `arc::Arc`
- `r::Float64`: System per-unit value, validation range: (-2, 4), action if invalid: error
- `x::Float64`: System per-unit value, validation range: (-2, 4), action if invalid: error
- `primaryshunt::Float64`: System per-unit value, validation range: (0, 2), action if invalid: error
- `rate::Union{Nothing, Float64}`, validation range: (0, nothing), action if invalid: error
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct Transformer2W <: ACBranch
    name::String
    available::Bool
    activepower_flow::Float64
    reactivepower_flow::Float64
    arc::Arc
    "System per-unit value"
    r::Float64
    "System per-unit value"
    x::Float64
    "System per-unit value"
    primaryshunt::Float64
    rate::Union{Nothing, Float64}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function Transformer2W(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    Transformer2W(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, services, ext, forecasts, InfrastructureSystemsInternal(), )
end

function Transformer2W(; name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, services=Device[], ext=Dict{String, Any}(), forecasts=InfrastructureSystems.Forecasts(), )
    Transformer2W(name, available, activepower_flow, reactivepower_flow, arc, r, x, primaryshunt, rate, services, ext, forecasts, )
end

# Constructor for demo purposes; non-functional.
function Transformer2W(::Nothing)
    Transformer2W(;
        name="init",
        available=false,
        activepower_flow=0.0,
        reactivepower_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        r=0.0,
        x=0.0,
        primaryshunt=0.0,
        rate=nothing,
        services=Device[],
        ext=Dict{String, Any}(),
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::Transformer2W) = value.name
"""Get Transformer2W available."""
get_available(value::Transformer2W) = value.available
"""Get Transformer2W activepower_flow."""
get_activepower_flow(value::Transformer2W) = value.activepower_flow
"""Get Transformer2W reactivepower_flow."""
get_reactivepower_flow(value::Transformer2W) = value.reactivepower_flow
"""Get Transformer2W arc."""
get_arc(value::Transformer2W) = value.arc
"""Get Transformer2W r."""
get_r(value::Transformer2W) = value.r
"""Get Transformer2W x."""
get_x(value::Transformer2W) = value.x
"""Get Transformer2W primaryshunt."""
get_primaryshunt(value::Transformer2W) = value.primaryshunt
"""Get Transformer2W rate."""
get_rate(value::Transformer2W) = value.rate
"""Get Transformer2W services."""
get_services(value::Transformer2W) = value.services
"""Get Transformer2W ext."""
get_ext(value::Transformer2W) = value.ext

InfrastructureSystems.get_forecasts(value::Transformer2W) = value.forecasts
"""Get Transformer2W internal."""
get_internal(value::Transformer2W) = value.internal


text:
###
### Mer specific specializations of src/biosequence/predicates.jl
###

Base.cmp(seq1::T, seq2::T) where {T<:AbstractMer} = cmp(encoded_data(seq1), encoded_data(seq2))
Base.:(==)(x::T, y::T) where {T<:AbstractMer} = encoded_data(x) == encoded_data(y)
Base.isless(x::T, y::T) where {T<:AbstractMer} = isless(encoded_data(x), encoded_data(y))

function Base.hash(x::Mer{<:NucleicAcidAlphabet{2},K}, h::UInt) where {K}
    return Base.hash_uint64(encoded_data(x) ⊻ K ⊻ h)
end


text:
using DataFrames: DataFrames, DataFrame
using Dates: Dates, Date
using TimeZones
using TypedTables: TypedTables, Table

@testset "store" begin
    @testset "dataframe" begin
        @system SStoreDataFrame(DataFrameStore, Controller) begin
            a(s) => s[:a] ~ track::int
            b(s) => s[:b] ~ track
        end
        a = [1, 2, 3]
        b = [4.0, 5.0, 6.0]
        df = DataFrame(; a, b)
        n = DataFrames.nrow(df)
        r = simulate(SStoreDataFrame, config=:0 => (:df => df, :ik => :a), stop=n-1)
        @test r.a == a
        @test r.b == b
    end

    @testset "day" begin
        @system SStoreDay(DayStore, Controller) begin
            a(s) => s.a ~ track
        end
        r = mktemp() do f, io
            write(io, "day (d),a\n0,0\n1,10\n2,20")
            close(io)
            c = (:Clock => :step => 1u"d", :0 => :filename => f)
            simulate(SStoreDay, config=c, stop=2)
        end
        @test r.a[1] == 0
        @test r.a[2] == 10
        @test r.a[end] == 20
    end

    @testset "date" begin
        @system SStoreDate(DateStore, Controller) begin
            a(s) => s.a ~ track
        end
        r = mktemp() do f, io
            write(io, "date (:Date),a\n2020-12-09,0\n2020-12-10,10\n2020-12-11,20")
            close(io)
            config = (
                :Clock => :step => 1u"d",
                :Calendar => :init => ZonedDateTime(2020, 12, 9, tz"UTC"),
                :0 => :filename => f,
            )
            simulate(SStoreDate; config, stop=2u"d")
        end
        @test r.a[1] == 0
        @test r.a[2] == 10
        @test r.a[end] == 20
    end

    @testset "time" begin
        @system SStoreTime(TimeStore, Controller) begin
            a(s) => s.a ~ track
        end
        r = mktemp() do f, io
            write(io, "date (:Date),time (:Time),a\n2020-11-15,01:00,0\n2020-11-15,02:00,10\n2020-11-15,03:00,20")
            close(io)
            config = (
                :Clock => :step => 1u"hr",
                :Calendar => :init => ZonedDateTime(2020, 11, 15, 1, tz"America/Los_Angeles"),
                :0 => (:filename => f, :tz => tz"America/Los_Angeles"),
            )
            simulate(SStoreTime; config, stop=2u"hr")
        end
        @test r.a[1] == 0
        @test r.a[2] == 10
        @test r.a[end] == 20
    end

    @testset "table" begin
        @system SStoreTable(TableStore, Controller) begin
            a(s) => s[:a] ~ track::int
            b(s) => s[:b] ~ track
        end
        a = [1, 2, 3]
        b = [4.0, 5.0, 6.0]
        tb = Table(; a, b)
        n = length(tb)
        r = simulate(SStoreTable, config=:0 => :tb => tb, stop=n-1)
        @test r.a == a
        @test r.b == b
    end
end


text:
# Constructors:
###############
"""
    GBVector{T}(n = libgb.GxB_INDEX_MAX)
"""
function GBVector{T}(n = libgb.GxB_INDEX_MAX) where {T}
    return GBVector{T}(libgb.GrB_Vector_new(toGBType(T),n))
end

GBVector{T}(dims::Dims{1}) where {T} = GBVector{T}(dims...)
GBVector{T}(nrows::Base.OneTo) where {T} =
    GBVector{T}(nrows.stop)
GBVector{T}(nrows::Tuple{Base.OneTo,}) where {T} = GBVector{T}(first(nrows))
"""
    GBVector(I::Vector, X::Vector{T})

Create a GBVector from a vector of indices `I` and a vector of values `X`.
"""
function GBVector(I::AbstractVector{U}, X::AbstractVector{T}; dup = BinaryOps.PLUS, nrows = maximum(I)) where {U<:Integer, T}
    x = GBVector{T}(nrows)
    build(x, I, X, dup = dup)
    return x
end

#iso valued constructors.
"""
    GBVector(I, x; nrows = maximum(I))

Create an `n` length GBVector `v` such that `M[I[k]] = x`.
The resulting vector is "iso-valued" such that it only stores `x` once rather than once for
each index.
"""
function GBVector(I::AbstractVector{U}, x::T;
    nrows = maximum(I)) where {U<:Integer, T}
    A = GBVector{T}(nrows)
    build(A, I, x)
    return A
end

"""
    GBVector(n, x)

Create an `n` length dense GBVector `v` such that M[I[k]] = x.
The resulting vector is "iso-valued" such that it only stores `x` once rather than once for
each index.
"""
function GBVector(n::Integer, x::T) where {T}
    v = GBVector{T}(n)
    v[:] = x
    return v
end

# Some Base and basic SparseArrays/LinearAlgebra functions:
###########################################################
Base.unsafe_convert(::Type{libgb.GrB_Vector}, v::GBVector) = v.p

function Base.copy(v::GBVector{T}) where {T}
    return GBVector{T}(libgb.GrB_Vector_dup(v))
end

clear!(v::GBVector) = libgb.GrB_Vector_clear(v)

function Base.size(v::GBVector)
    return (Int64(libgb.GrB_Vector_size(v)),)
end

SparseArrays.nnz(v::GBVector) = Int64(libgb.GrB_Vector_nvals(v))
Base.eltype(::Type{GBVector{T}}) where{T} = T

function Base.similar(
    ::GBVector{T}, ::Type{TNew},
    dims::Dims{1}
) where {T, TNew}
    return GBVector{TNew}(dims...)
end

function Base.deleteat!(v::GBVector, i)
    libgb.GrB_Vector_removeElement(v, libgb.GrB_Index(i))
    return v
end

function Base.resize!(v::GBVector, n)
    libgb.GrB_Vector_resize(v, n)
    return v
end

function LinearAlgebra.diag(A::GBMatrix{T}, k::Integer = 0; desc = C_NULL) where {T}
    return GBVector{T}(libgb.GxB_Vector_diag(A, k, desc))
end

#We need these until I can get a SparseArrays.nonzeros implementation
function Base.show(io::IO, ::MIME"text/plain", v::GBVector)
    gxbprint(io, v)
end

function Base.show(io::IO, v::GBVector)
    gxbprint(io, v)
end

function Base.show(io::IOContext, v::GBVector)
    gxbprint(io, v)
end

# Type dependent functions build, setindex, getindex, and findnz:
for T ∈ valid_vec
    if T ∈ gxb_vec
        prefix = :GxB
    else
        prefix = :GrB
    end

    # Build functions
    func = Symbol(prefix, :_Vector_build_, suffix(T))
    @eval begin
        function build(v::GBVector{$T}, I::Vector, X::Vector{$T}; dup = BinaryOps.PLUS)
            nnz(v) == 0 || throw(libgb.OutputNotEmptyError("Cannot build vector with existing elements"))
            length(X) == length(I) || DimensionMismatch("I and X must have the same length")
            libgb.$func(v, Vector{libgb.GrB_Index}(I) .- 1, X, length(X), dup[$T])
        end
    end
    # Setindex functions
    func = Symbol(prefix, :_Vector_setElement_, suffix(T))
    @eval begin
        function Base.setindex!(v::GBVector{$T}, x, i::Integer)
            x = convert($T, x)
            return libgb.$func(v, x, libgb.GrB_Index(i) - 1)
        end
    end
    # Getindex functions
    func = Symbol(prefix, :_Vector_extractElement_, suffix(T))
    @eval begin
        function Base.getindex(v::GBVector{$T}, i::Integer)
            x = Ref{$T}()
            result = libgb.$func(x, v, libgb.GrB_Index(i) - 1)
            if result == libgb.GrB_SUCCESS
                return x[]
            elseif result == libgb.GrB_NO_VALUE
                return nothing
            else
                throw(ErrorException("Invalid extractElement return value."))
            end
        end
    end
    # findnz functions
    func = Symbol(prefix, :_Vector_extractTuples_, suffix(T))
    @eval begin
        function SparseArrays.findnz(v::GBVector{$T})
            nvals = Ref{libgb.GrB_Index}(nnz(v))
            I = Vector{libgb.GrB_Index}(undef, nvals[])
            X = Vector{$T}(undef, nvals[])
            libgb.$func(I, X, nvals, v)
            nvals[] == length(I) == length(X) || throw(DimensionMismatch("length(I) != length(X)"))
            return I .+ 1, X
        end
        function SparseArrays.nonzeros(v::GBVector{$T})
            nvals = Ref{libgb.GrB_Index}(nnz(v))
            X = Vector{$T}(undef, nvals[])
            libgb.$func(C_NULL, X, nvals, v)
            nvals[] == length(X) || throw(DimensionMismatch(""))
            return X
        end
        function SparseArrays.nonzeroinds(v::GBVector{$T})
            nvals = Ref{libgb.GrB_Index}(nnz(v))
            I = Vector{libgb.GrB_Index}(undef, nvals[])
            libgb.$func(I, C_NULL, nvals, v)
            nvals[] == length(I) || throw(DimensionMismatch(""))
            return I .+ 1
        end
    end
end

function build(v::GBVector{T}, I::Vector, x::T) where {T}
    nnz(v) == 0 || throw(libgb.OutputNotEmptyError("Cannot build vector with existing elements"))
    x = GBScalar(x)
    return libgb.GxB_Vector_build_Scalar(
            v,
            Vector{libgb.GrB_Index}(I),
            x,
            length(I)
        )
end

# Indexing functions:
#####################
function _outlength(u, I)
    if I == ALL
        wlen = size(u)
    else
        wlen = length(I)
    end
    return wlen
end

"""
    extract!(w::GBVector, u::GBVector, I; kwargs...)::GBVector

Extract a subvector from `u` into the output vector `w`. Equivalent to the matrix definition.
"""
function extract!(
    w::GBVector, u::GBVector, I;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    I, ni = idx(I)
    libgb.GrB_Vector_extract(w, mask, getaccum(accum, eltype(w)), u, I, ni, desc)
    return w
end

function extract!(
    w::GBVector, u::GBVector, ::Colon;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    return extract!(w, u, ALL; mask, accum, desc)
end

"""
    extract(u::GBVector, I; kwargs...)::GBVector

Extract a subvector from `u` and return it. Equivalent to the matrix definition.
"""
function extract(
    u::GBVector, I;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    wlen = _outlength(u, I)
    w = similar(u, wlen)
    return extract!(w, u, I; mask, accum, desc)
end

function extract(u::GBVector, ::Colon; mask = C_NULL, accum = nothing, desc=DEFAULTDESC)
    extract(u, ALL; mask, accum, desc)
end

function Base.getindex(
    u::GBVector, I;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    return extract(u, I; mask, accum, desc)
end

function Base.getindex(u::GBVector, ::Colon; mask = C_NULL, accum = nothing, desc = DEFAULTDESC)
    return extract(u, :)
end

function Base.getindex(
    u::GBVector, i::Union{Vector, UnitRange, StepRange};
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    return extract(u, i; mask, accum, desc)
end
"""
    subassign(w::GBVector, u::GBVector, I; kwargs...)::GBVector

Assign a subvector of `w` to `u`. Return `u`. Equivalent to the matrix definition.
"""
function subassign!(
    w::GBVector, u, I;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    I, ni = idx(I)
    u isa Vector && (u = GBVector(u))
    if u isa GBVector
        libgb.GxB_Vector_subassign(w, mask, getaccum(accum, eltype(w)), u, I, ni, desc)
    else
        libgb.scalarvecsubassign[eltype(u)](w, mask, getaccum(accum, eltype(w)), u, I, ni, desc)
    end
    return nothing
end

"""
    assign(w::GBVector, u::GBVector, I; kwargs...)::GBVector

Assign a subvector of `w` to `u`. Return `u`. Equivalent to the matrix definition.
"""
function assign!(
    w::GBVector, u, I;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    I, ni = idx(I)
    u isa Vector && (u = GBVector(u))
    if u isa GBVector
        libgb.GrB_Vector_assign(w, mask, getaccum(accum, eltype(w)), u, I, ni, desc)
    else
        libgb.scalarvecassign[eltype(u)](w, mask, getaccum(accum, eltype(w)), u, I, ni, desc)
    end
    return nothing
end

function Base.setindex!(
    u::GBVector, x, ::Colon;
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    subassign!(u, x, ALL; mask, accum, desc)
    return nothing
end
function Base.setindex!(
    u::GBVector, x, I::Union{Vector, UnitRange, StepRange};
    mask = C_NULL, accum = nothing, desc = DEFAULTDESC
)
    subassign!(u, x, I; mask, accum, desc)
    return nothing
end


text:
using Test
using POMDPModels

let
    rng = MersenneTwister(41)
    p = SimpleLightDark()
    dp = DSimpleLightDark()

    @test discount(p) == 0.95
    @test discount(dp) == 0.95

    s0 = rand(rng, initialstate(p))
    s0 = rand(rng, initialstate(dp))

    s0 = 0
    s0, _, r = @gen(:sp, :o, :r)(p, s0, +1, rng)
    @test s0 == 1
    @test r == -1.0
    s1, _, r = @gen(:sp, :o, :r)(p, s0, 0, rng)
    @test isterminal(p, s1) == true
    @test r == -100.0
    s2 = 10
    obs = rand(rng, observation(p, s2))
    @test abs(obs-10.0) <= 0.1
    obs = rand(rng, observation(dp, s2))
    @test pdf(observation(dp, s2), 71) == pdf(observation(dp, s2), 72)
end


text:
###########################################################################
# Gauss-Kronrod integration-weight computation for arbitrary floating-point
# types and precision, implemented based on the description in:
#
#    Dirk P. Laurie, "Calculation of Gauss-Kronrod quadrature rules,"
#    Mathematics of Computation, vol. 66, no. 219, pp. 1133-1145 (1997).
#
# for the Kronrod rule, and for the Gauss rule from the description in
#
#    Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997).
#
# Arbitrary-precision eigenvalue (eignewt & eigpoly) and eigenvector
# (eigvec1) routines are written by SGJ, independent of the above sources.
#
# Since we only implement Gauss-Kronrod rules for the unit weight function,
# the diagonals of the Jacobi matrices are zero and certain things simplify
# compared to the general case of an arbitrary weight function.

# Given a symmetric tridiagonal matrix H with H[i,i] = 0 and
# H[i-1,i] = H[i,i-1] = b[i-1], compute p(z) = det(z I - H) and its
# derivative p'(z), returning (p,p').
function eigpoly(b,z,m=length(b)+1)
    d1 = z
    d1deriv = d2 = one(z)
    d2deriv = zero(z)
    for i = 2:m
        b2 = b[i-1]^2
        d = z * d1 - b2 * d2
        dderiv = d1 + z * d1deriv - b2 * d2deriv
        d2 = d1
        d1 = d
        d2deriv = d1deriv
        d1deriv = dderiv
    end
    return (d1, d1deriv)
end

# compute the n smallest eigenvalues of the symmetric tridiagonal matrix H
# (defined from b as in eigpoly) using a Newton iteration
# on det(H - lambda I).  Unlike eig, handles BigFloat.
function eignewt(b,m,n)
    # get initial guess from eig on Float64 matrix
    H = SymTridiagonal(zeros(m), Float64[ b[i] for i in 1:m-1 ])
    lambda0 = sort(eigvals(H))

    lambda = Array{eltype(b)}(undef, n)
    for i = 1:n
        lambda[i] = lambda0[i]
        for k = 1:1000
            (p,pderiv) = eigpoly(b,lambda[i],m)
            lambda[i] = (lamold = lambda[i]) - p / pderiv
            if abs(lambda[i] - lamold) < 10 * eps(lambda[i]) * abs(lambda[i])
                break
            end
        end
        # do one final Newton iteration for luck and profit:
        (p,pderiv) = eigpoly(b,lambda[i],m)
        lambda[i] = lambda[i] - p / pderiv
    end
    return lambda
end

# given an eigenvalue z and the matrix H(b) from above, return
# the corresponding eigenvector, normalized to 1.
function eigvec1(b,z::Number,m=length(b)+1)
    # "cheat" and use the fact that our eigenvector v must have a
    # nonzero first entries (since it is a quadrature weight), so we
    # can set v[1] = 1 to solve for the rest of the components:.
    v = Array{eltype(b)}(undef, m)
    v[1] = 1
    if m > 1
        s = v[1]
        v[2] = z * v[1] / b[1]
        s += v[2]^2
        for i = 3:m
            v[i] = - (b[i-2]*v[i-2] - z*v[i-1]) / b[i-1]
            s += v[i]^2
        end
        rmul!(v, 1 / sqrt(s))
    end
    return v
end

"""
    gauss([T,] N, a=-1, b=1)

Return a pair `(x, w)` of `N` quadrature points `x[i]` and weights `w[i]` to
integrate functions on the interval `(a, b)`,  i.e. `sum(w .* f.(x))`
approximates the integral.  Uses the method described in Trefethen &
Bau, Numerical Linear Algebra, to find the `N`-point Gaussian quadrature
in O(`N`²) operations.

`T` is an optional parameter specifying the floating-point type, defaulting
to `Float64`. Arbitrary precision (`BigFloat`) is also supported.
"""
function gauss(::Type{T}, N::Integer) where T<:AbstractFloat
    if N < 1
        throw(ArgumentError("Gauss rules require positive order"))
    end
    o = one(T)
    b = T[ n / sqrt(4n^2 - o) for n = 1:N-1 ]
    x = eignewt(b,N,N)
    w = T[ 2*eigvec1(b,x[i])[1]^2 for i = 1:N ]
    return (x, w)
end

gauss(N::Integer) = gauss(Float64, N) # integration on the standard interval (-1,1)

# re-scaled to an arbitrary interval:
gauss(N::Integer, a::Real, b::Real) = gauss(typeof(float(b-a)), N, a, b)
function gauss(::Type{T}, N::Integer, a::Real, b::Real) where T<:AbstractFloat
    x, w = gauss(T, N)
    s = T(b-a)/2
    x .= a .+ (x .+ 1) .* s
    w .*= abs(s)
    return (x, w)
end

"""
    kronrod([T,] n)

Compute `2n+1` Kronrod points `x` and weights `w` based on the description in
Laurie (1997), appendix A, simplified for `a=0`, for integrating on `[-1,1]`.
Since the rule is symmetric, this only returns the `n+1` points with `x <= 0`.
The function Also computes the embedded `n`-point Gauss quadrature weights `gw`
(again for `x <= 0`), corresponding to the points `x[2:2:end]`. Returns `(x,w,wg)`
in O(`n`²) operations.

`T` is an optional parameter specifying the floating-point type, defaulting
to `Float64`. Arbitrary precision (`BigFloat`) is also supported.

Given these points and weights, the estimated integral `I` and error `E` can
be computed for an integrand `f(x)` as follows:

    x, w, wg = kronrod(n)
    fx⁰ = f(x[end])                # f(0)
    x⁻ = x[1:end-1]                # the x < 0 Kronrod points
    fx = f.(x⁻) .+ f.((-).(x⁻))    # f(x < 0) + f(x > 0)
    I = sum(fx .* w[1:end-1]) + fx⁰ * w[end]
    if isodd(n)
        E = abs(sum(fx[2:2:end] .* wg[1:end-1]) + fx⁰*wg[end] - I)
    else
        E = abs(sum(fx[2:2:end] .* wg[1:end])- I)
    end
"""
function kronrod(::Type{T}, n::Integer) where T<:AbstractFloat
    if n < 1
        throw(ArgumentError("Kronrod rules require positive order"))
    end
    o = one(T)
    b = zeros(T, 2n+1)
    b[1] = 2*o
    for j = 1:div(3n+1,2)
        b[j+1] = j^2 / (4j^2 - o)
    end
    s = zeros(T, div(n,2) + 2)
    t = zeros(T, div(n,2) + 2)
    t[2] = b[n+2]
    for m = 0:n-2
        u = zero(T)
        for k = div(m+1,2):-1:0
            l = m - k + 1
            k1 = k + n + 2
            u += b[k1]*s[k+1] - b[l]*s[k+2]
            s[k+2] = u
        end
        s,t = t,s
    end
    for j = div(n,2):-1:0
        s[j+2] = s[j+1]
    end
    for m = n-1:2n-3
        u = zero(T)
        for k = m+1-n:div(m-1,2)
            l = m - k + 1
            j = n - l
            k1 = k + n + 2
            u -= b[k1]*s[j+2] - b[l]*s[j+3]
            s[j+2] = u
        end
        k = div(m+1,2)
        if 2k != m
            j = n - (m - k + 2)
            b[k+n+2] = s[j+2] / s[j+3]
        end
        s,t = t,s
    end
    for j = 1:2n
        b[j] = sqrt(b[j+1])
    end

    # get negative quadrature points x
    x = eignewt(b,2n+1,n+1) # x <= 0

    # get quadrature weights
    w = T[ 2*eigvec1(b,x[i],2n+1)[1]^2 for i in 1:n+1 ]

    # Get embedded Gauss rule from even-indexed points, using
    # the method described in Trefethen and Bau.
    for j = 1:n-1
        b[j] = j / sqrt(4j^2 - o)
    end
    gw = T[ 2*eigvec1(b,x[i],n)[1]^2 for i = 2:2:n+1 ]

    return (x, w, gw)
end

kronrod(N::Integer) = kronrod(Float64, N)

###########################################################################
# Type-stable cache of quadrature rule results, so that we don't
# repeat the kronrod calculation unnecessarily.

# precomputed n=7 rule in double precision (computed in 100-bit arithmetic),
# since this is the common case.
const xd7 = [-9.9145537112081263920685469752598e-01,
             -9.4910791234275852452618968404809e-01,
             -8.6486442335976907278971278864098e-01,
             -7.415311855993944398638647732811e-01,
             -5.8608723546769113029414483825842e-01,
             -4.0584515137739716690660641207707e-01,
             -2.0778495500789846760068940377309e-01,
             0.0]
const wd7 = [2.2935322010529224963732008059913e-02,
             6.3092092629978553290700663189093e-02,
             1.0479001032225018383987632254189e-01,
             1.4065325971552591874518959051021e-01,
             1.6900472663926790282658342659795e-01,
             1.9035057806478540991325640242055e-01,
             2.0443294007529889241416199923466e-01,
             2.0948214108472782801299917489173e-01]
const gwd7 = [1.2948496616886969327061143267787e-01,
              2.797053914892766679014677714229e-01,
              3.8183005050511894495036977548818e-01,
              4.1795918367346938775510204081658e-01]

# cache of T -> n -> (x,w,gw) Kronrod rules, to avoid recomputing them
# unnecessarily for repeated integration.   We initialize it with the
# default n=7 rule for double-precision calculations.  We use a cache
# of caches to allow us to evaluate the cache in a type-stable way with
# a generated function below.
const rulecache = Dict{Type,Dict}(
    Float64 => Dict{Int,NTuple{3,Vector{Float64}}}(7 => (xd7,wd7,gwd7)),
    Float32 => Dict{Int,NTuple{3,Vector{Float32}}}(7 => (xd7,wd7,gwd7)))

# for BigFloat rules, we need a separate cache keyed by (n,precision)
const bigrulecache = Dict{Tuple{Int,Int}, NTuple{3,Vector{BigFloat}}}()

cachedrule(::Type{T}, n::Integer) where {T<:Number} = cachedrule(typeof(float(real(one(T)))), n::Integer)

function cachedrule(::Type{BigFloat}, n::Integer)
    key = (n, precision(BigFloat))
    haskey(bigrulecache, key) ? bigrulecache[key] : (bigrulecache[key] = kronrod(BigFloat, n))
end

# use a generated function to make this type-stable
@generated function cachedrule(::Type{T}, n::Integer) where {T<:AbstractFloat}
    cache = haskey(rulecache, T) ? rulecache[T] : (rulecache[T] = Dict{Int,NTuple{3,Vector{T}}}())
    :(haskey($cache, n) ? $cache[n] : ($cache[n] = kronrod($T, n)))
end


text:
using ImagineInterface, ImagineFormat

ais = parse_ai("t.ai")

# Grab one particular signal
piezo = getname(ais, "axial piezo monitor")
# piezo is just a reference, the data are loaded on-demand. This lets you
# work with long recordings.
# Extract values in physical units, which here represent the position
data = get_samples(piezo)
# We can also get them in their original voltage units...
datav = get_samples(piezo; sampmap=:volts)
# ...or even in the raw Int16 format
dataraw = get_samples(piezo; sampmap=:raw)

# Let's display this signal
using ImaginePlots   # you have to install this and its dependencies manually
using Plots
plot(piezo)

stimuli = getname(ais, "stimuli")
# This channel is just noise, but if it had had a sequence of TTL pulses
# this would give us a list of scan #s at which the pulses start
stimstarts = find_pulse_starts(stimuli; sampmap=:volts)

max_forks_repo_name: HolyLab/ImagineInterface


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

abstract type AbstractMCMCTunerConfig end
export AbstractMCMCTunerConfig

abstract type AbstractMCMCTuner end
export AbstractMCMCTuner

function mcmc_init end
export mcmc_init

function mcmc_tune_burnin! end
export mcmc_tune_burnin!

function isviable end

@with_kw struct MCMCInitStrategy
    ninit_tries_per_chain::ClosedInterval{Int64} = 8..128
    max_nsamples_pretune::Int64 = 25
    max_nsteps_pretune::Int64 = 250
    max_time_pretune::Float64 = Inf
end

export MCMCInitStrategy

MCMCInitStrategy(tuner_config::AbstractMCMCTunerConfig) =
    MCMCInitStrategy()

gen_tuners(ids::Range{<:Integer},
    chainspec::MCMCSpec,
    exec_context::ExecContext,
    tuner_config::AbstractMCMCTunerConfig,
    ) = [tuner_config(chainspec(id, exec_context), init_proposal = true) for id in ids]

function mcmc_init(
    chainspec::MCMCSpec,
    nchains::Int,
    exec_context::ExecContext = ExecContext(),
    tuner_config::AbstractMCMCTunerConfig = AbstractMCMCTunerConfig(chainspec.algorithm),
    convergence_test::MCMCConvergenceTest = BGConvergence(),
    init_strategy::MCMCInitStrategy = MCMCInitStrategy(tuner_config);
    ll::LogLevel = LOG_INFO
)
    @log_msg ll "Trying to generate $nchains viable MCMC chain(s)."

    min_nviable = minimum(init_strategy.ninit_tries_per_chain) * nchains
    max_ncandidates = maximum(init_strategy.ninit_tries_per_chain) * nchains

    ncandidates = zero(Int64)

    tuners = similar([tuner_config(chainspec(zero(Int64), exec_context), init_proposal = false)], zero(Int))
    cycle = one(Int)
    while length(tuners) < min_nviable && ncandidates < max_ncandidates
        n = min(min_nviable, max_ncandidates - ncandidates)
        @log_msg ll+1 "Generating $n $(cycle > 1 ? "additional " : "")MCMC chain(s)."
        new_tuners = gen_tuners(ncandidates + (one(Int64):n), chainspec, exec_context, tuner_config)
        ncandidates += n

        @log_msg ll+1 "Testing $(length(new_tuners)) MCMC chain(s)."

        chains = map(x -> x.chain, new_tuners)

        run_tuning_iterations!(
            (), new_tuners, exec_context;
            max_nsamples = max(5, div(init_strategy.max_nsamples_pretune, 5)),
            max_nsteps =  max(50, div(init_strategy.max_nsteps_pretune, 5)),
            max_time = init_strategy.max_time_pretune / 5,
            ll = ll+2
        )

        filter!(isviable, new_tuners)
        @log_msg ll+1 "Found $(length(new_tuners)) viable MCMC chain(s)."

        if !isempty(new_tuners)
            run_tuning_iterations!(
                (), new_tuners, exec_context;
                max_nsamples = init_strategy.max_nsamples_pretune,
                max_nsteps = init_strategy.max_nsteps_pretune,
                max_time = init_strategy.max_time_pretune,
                ll = ll+2
            )

            nsamples_thresh = floor(Int, 0.8 * median([nsamples(t.chain.state) for t in new_tuners]))

            filter!(t -> nsamples(t.chain.state) >= nsamples_thresh, new_tuners)
            @log_msg ll+1 "Found $(length(new_tuners)) MCMC chain(s) with at least $(nsamples_thresh) samples."

            append!(tuners, new_tuners)
        end

        cycle += 1
    end

    length(tuners) < min_nviable && error("Failed to generate $min_nviable viable MCMC chains")

    m = nchains
    tidxs = linearindices(tuners)
    n = length(tidxs)

    mode_1 = tuners[1].stats.mode
    modes = Array{eltype(mode_1)}(length(mode_1), n)
    for i in tidxs
        modes[:,i] = tuners[i].stats.mode
    end

    final_tuners = similar(tuners, zero(Int))

    if 2 <= m < size(modes, 2)
        clusters = kmeans(modes, m)
        clusters.converged || error("k-means clustering of MCMC chains did not converge")

        mincosts = fill(Inf, m)
        tuner_sel = fill(zero(Int), m)

        for i in tidxs
            j = clusters.assignments[i]
            if clusters.costs[i] < mincosts[j]
                mincosts[j] = clusters.costs[i]
                tuner_sel[j] = i
            end
        end

        @assert all(j -> j in tidxs, tuner_sel)

        for i in sort(tuner_sel)
            push!(final_tuners, tuners[i])
        end
    else
        @assert length(tuners) == nchains
        resize!(final_tuners, nchains)
        copy!(final_tuners, tuners)
    end

    @log_msg ll "Selected $(length(final_tuners)) MCMC chain(s)."

    final_tuners
end

@with_kw struct MCMCBurninStrategy
    max_nsamples_per_cycle::Int64 = 1000
    max_nsteps_per_cycle::Int64 = 10000
    max_time_per_cycle::Float64 = Inf
    max_ncycles::Int = 30
end

export MCMCBurninStrategy

MCMCBurninStrategy(tuner_config::AbstractMCMCTunerConfig) =
    MCMCBurninStrategy()

MCMCBurninStrategy(tuner::AbstractMCMCTuner) =
    MCMCBurninStrategy()

function mcmc_tune_burnin!(
    callbacks,
    tuners::AbstractVector{<:AbstractMCMCTuner},
    convergence_test::MCMCConvergenceTest,
    burnin_strategy::MCMCBurninStrategy,
    exec_context::ExecContext;
    strict_mode::Bool = false,
    ll::LogLevel = LOG_INFO
)
    @log_msg ll "Begin tuning of $(length(tuners)) MCMC chain(s)."

    chains = map(x -> x.chain, tuners)
    nchains = length(chains)

    user_callbacks = mcmc_callback_vector(callbacks, chains)

    cycles = zero(Int)
    successful = false
    while !successful && cycles < burnin_strategy.max_ncycles
        cycles += 1
        run_tuning_cycle!(
            user_callbacks, tuners, exec_context;
            max_nsamples = burnin_strategy.max_nsamples_per_cycle,
            max_nsteps = burnin_strategy.max_nsteps_per_cycle,
            max_time = burnin_strategy.max_time_per_cycle,
            ll = ll+2
        )

        stats = [x.stats for x in tuners] # ToDo: Find more generic abstraction
        ct_result = check_convergence!(convergence_test, chains, stats, ll = ll+2)

        ntuned = count(c -> c.tuned, chains)
        nconverged = count(c -> c.converged, chains)
        successful = (ntuned == nconverged == nchains)

        for i in eachindex(user_callbacks, tuners)
            user_callbacks[i](1, tuners[i])
        end

        @log_msg ll+1 "MCMC Tuning cycle $cycles finished, $nchains chains, $ntuned tuned, $nconverged converged."
    end

    if successful
        @log_msg ll "MCMC tuning of $nchains chains successful after $cycles cycle(s)."
    else
        msg = "MCMC tuning of $nchains chains aborted after $cycles cycle(s)."
        if strict_mode
            error(msg)
        else
            @log_msg LOG_WARNING msg
        end
    end

    successful
end

struct NoOpTunerConfig <: BAT.AbstractMCMCTunerConfig end
export NoOpTunerConfig

(config::NoOpTunerConfig)(chain::MCMCIterator; kwargs...) =
    NoOpTuner(chain)

struct NoOpTuner{C<:MCMCIterator} <: AbstractMCMCTuner
    chain::C
end

export NoOpTuner

MCMCInitStrategy(tuner_config::NoOpTunerConfig) =
    MCMCInitStrategy(1..1, zero(Int), zero(Int), Inf)

MCMCBurninStrategy(tuner_config::NoOpTunerConfig) =
    MCMCBurninStrategy(zero(Int), zero(Int), Inf, zero(Int))

isviable(tuner::NoOpTuner) = true

function mcmc_init(
    chainspec::MCMCSpec,
    nchains::Int,
    exec_context::ExecContext,
    tuner_config::NoOpTunerConfig,
    convergence_test::MCMCConvergenceTest,
    init_strategy::MCMCInitStrategy;
    ll::LogLevel = LOG_INFO
)
    @log_msg ll "NoOpTuner generating $nchains MCMC chain(s)."

    [tuner_config(chainspec(id, exec_context), init_proposal = true) for id in one(Int64):Int64(nchains)]
end

function mcmc_tune_burnin!(
    callbacks,
    tuners::AbstractVector{<:NoOpTuner},
    convergence_test::MCMCConvergenceTest,
    burnin_strategy::MCMCBurninStrategy,
    exec_context::ExecContext;
    ll::LogLevel = LOG_INFO,
    kwargs...
)
    @log_msg ll "Tune/Burn-In with NoOpTuner doing nothing."
end


text:
varname_cplx2real(varname::String) = (varname*"_Re",varname*"_Im")

function Base.show(io::IO, x::Variable)
    print(io, x.name)
end

function Base.show(io::IO, d::Degree)
    print(io, "(", d.explvar, ",", d.conjvar, ")")
end

function Base.print(io::IO, exp::Exponent)
  if exp == Exponent()
    print(io, "1")
  else
    expo = exp.expo
    sortedCollec = sort(collect(expo), by=(x)->x[1].name)
    i=length(sortedCollec)
    for (var, deg) in sortedCollec
      if var.kind <: Complex && deg.conjvar>0
        print(io, "conj(", var, ")")
        if deg.conjvar > 1
          print(io, "^", deg.conjvar)
        end
        if deg.explvar > 0
          print(io, " * ")
        end
      end
      if deg.explvar == 1
        print(io, var)
      elseif deg.explvar > 1
        print(io, var, "^", deg.explvar)
      end
      if i > 1
        print(io, " * ")
      end
      i -= 1
    end
  end
end

function Base.print(io::IO, P::Polynomial)
  poly = P.poly
  i = length(poly)
  sorted_keys = sort(collect(keys(poly)))
  for expo in sorted_keys
    λ = poly[expo]
    if λ != 0
      print(io, "(", λ, ")")
    end
    if expo.degree != Degree(0,0)
      print(io, "*")
      print(io, expo)
      if i > 1
        print(io, " + ")
      end
      i -= 1
    end
  end
end

function Base.print(io::IO, pt::Point)
  for (var, val) in sort(collect(pt), by=x->x[1].name)
    println(io, var, " ", val)
  end
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
export ckpt_to_jld2, load_pretrained, list_pretrains, @pretrained_str

using JLD2, Requires, Pkg.TOML, DataDeps
using GoogleDrive: google_download
using Flux: loadparams!

const configs = open(TOML.parse, joinpath(@__DIR__, "pretrains.toml"))

"""     load_model(weights, config)
Loads the weights into a Flux model
"""
function load_model(weights, config)
    lstm_units = config["lstm_units"]
    input_dims = config["input_dims"]

    layers = []
    for i in 1:config["num_layers"]
        layer = BasicLSTM(
                     (i == 1) ? input_dims : lstm_units,
                     lstm_units
                )
        push!(layers, layer)
    end
    push!(layers, Dense(lstm_units, input_dims))

    model = Chain(layers...)

    weight_names = keys(weights)
    rnn_weights = filter(name -> occursin("rnn", name), weight_names)
    dense_weights = filter(name -> occursin("fully_connected", name), weight_names)

    for i in 1:length(rnn_weights)
        cell_weights = filter(name -> occursin("cell_$(i-1)", name), rnn_weights)
        for j in cell_weights
            if occursin("kernel", j)
                kernel = weights[j]
                if i == 1
                    Wi = @view kernel[:, 1:input_dims]
                    Wh = @view kernel[:, input_dims+1:end]
                else
                    Wi = @view kernel[:, 1:lstm_units]
                    Wh = @view kernel[:, lstm_units+1:end]
                end
                loadparams!(model[i], [Wi, Wh])
            elseif occursin("bias", j)
                loadparams!(model[i].cell.b, [weights[j]])
            end
        end
    end

    for j in dense_weights
        if occursin("weights", j)
            loadparams!(model[end].W, [weights[j]])
        elseif occursin("biases", j)
            loadparams!(model[end].b, [weights[j]])
        end
    end

    model
end

"""     ckpt2_to_jld2(ckptpath::String, ckptname::String, savepath::String)
Loads a pre-trained model from a tensorflow checkpoint and saves to jld2
"""
function ckpt_to_jld2(ckptpath::String; ckptname::String="perfrnn.ckpt", savepath::String="./")
    files = readdir(ckptpath)
    ckptname ∉ files && error("The checkpoint file $ckptname is not found")
    ckptname*".meta" ∉ files && error("The checkpoint meta file is not found")
    weights = readckpt(joinpath(ckptpath, ckptname))
    jld2name = normpath(joinpath(savepath, ckptname[1:end-5]*".jld2"))
    @info "Saving the model weights to $jld2name"
    JLD2.@save jld2name weights
end

# From Transformers.jl
"""     readckpt(path)
Load weights from a tensorflow checkpoint file into a Dict.
"""
readckpt(path) = error("readckpt require TensorFlow.jl installed. run `Pkg.add(\"TensorFlow\"); using TensorFlow`")

@init @require TensorFlow="1d978283-2c37-5f34-9a8e-e9c0ece82495" begin
    import .TensorFlow
    function readckpt(path)
        weights = Dict{String, Array}()
        TensorFlow.init()
        ckpt = TensorFlow.pywrap_tensorflow.x.NewCheckpointReader(path)
        shapes = ckpt.get_variable_to_shape_map()

        for (name, shape) ∈ shapes
            weight = ckpt.get_tensor(name)
            if length(shape) == 2
                weight = collect(weight')
            end
            weights[name] = weight
        end

        weights
    end
end

"""     load_pretrain(path)
Loads a pre-trained performance rnn model.
"""
function load_pretrained(model_name::String)
    if model_name ∉ keys(configs)
        error("""Invalid model.
               Please try list_pretrains() to check the available pre-trained models""")
    end

    model_path = @datadep_str("$model_name/$model_name.jld2")
    if !endswith(model_path, ".jld2")
        error("""Invalid file. A jld2 file is required to load the model.
                 If this is a tensorflow checkpoint file, run ckpt_to_jld2 to convert""")
    end

    JLD2.@load model_path weights
    model = load_model(weights, configs[model_name])
    num_velocitybins = configs[model_name]["velocity_bins"]
    perfencoder = PerformanceOneHotEncoding(num_velocitybins=num_velocitybins)
    perfrnn = PerfRNN(model, perfencoder, num_velocitybins)
end

"""     list_pretrains()
List all the available pre-trained models.
"""
function list_pretrains()
    println.(keys(configs))
    nothing
end

function register_configs(configs)
    for (model_name, config) in pairs(configs)
        model_desc = description(config["description"], config["host"], config["link"])
        checksum = config["checksum"]
        url = config["url"]
        dep = DataDep(model_name, model_desc, url, checksum;
                      fetch_method=google_download)
        DataDeps.register(dep)
    end
end

# From Transformers.jl
macro pretrained_str(name)
    :(load_pretrained($(esc(name))))
end

function description(description::String, host::String, link::String, cite=nothing)
  """
  $description
  Released by $(host) at $(link).
  $(isnothing(cite) ? "" : "\nCiting:\n$cite")"""
end

max_stars_repo_name: VasanthManiVasi/PerformanceRNN

max_issues_repo_name: VasanthManiVasi/PerformanceRNN

max_forks_repo_name: VasanthManiVasi/PerformanceRNN


text:
"""
Return a list of parameter names based on user supplied values.
"""
function parameter_names end

parameter_names(p, pnames::Nothing) = ["p$i" for i in eachindex(p)]

function parameter_names(p, pnames)
    if length(p) == length(pnames)
        return [string(pname) for pname in pnames]
    else
        throw(ArgumentError("Number of parameters does not match number of parameter names)"))
    end
end