text:
# Usually, running tests work with either
# 1) inside VSCode with ⌘⇧T (command+shift+T) (but plots will show up externally)
# 2) in shell=>Julia with `]test` (which work the same as above)
# 3) by launching the VSCode REPL and then `]test` (but now no plots are displayed)
#
# To overcome the plotting issue and have the plots displayed in the plots panel:
# - I added all the dependencies for the Package in the Package.toml of the test package
# - Replaced `using ChainPlot` with `include("../src/ChainPlot.jl")`
# - Open up `runtests.jl` in VSCode
# - "Julia: Change to This Directory"
# - "Julia: Activate This Environment" (dismiss popup warning if it appears)
# - "Julia: Execute File in REPL"
# voilá!

# using ChainPlot
using Flux
using Plots

include("../src/ChainPlot.jl")

gr()
theme(:default)

themes = [
    :default
    :dark
    :ggplot2 # somehow throwing error
    :juno
    :lime
    :orange
    :sand
    :solarized
    :solarized_light
    :wong
    :wong2
    :gruvbox_dark
    :gruvbox_light
]

using Colors
using Cairo
using Compose
using LightGraphs
using MetaGraphs
using GraphPlot

nnr = Chain(Dense(2,5,σ), RNN(5,4,relu), LSTM(4,4), GRU(4,4), Dense(4,3))
mg_nnr = ChainPlot.chaingraph(nnr)
@show get_prop(mg_nnr, 3, :layer_type)
@show get_prop(mg_nnr, 9, :index_in_layer)
@show get_prop(mg_nnr, 12, :neuron_color)
@show collect(edges(mg_nnr))
locs_x = [get_prop(mg_nnr, v, :loc_x) for v in vertices(mg_nnr)]
locs_y = [get_prop(mg_nnr, v, :loc_y) for v in vertices(mg_nnr)]
nodefillc = [parse(Colorant, get_prop(mg_nnr, v, :neuron_color)) for v in vertices(mg_nnr)]
draw(PNG("img/mg_nnr.png", 600, 400), gplot(mg_nnr, locs_x, locs_y, nodefillc=nodefillc))

dl = Dense(2,3)
display(plot(dl, title="$dl", titlefontsize=12))
savefig("img/dl.png")

rl = RNN(3,5)
display(plot(rl, title="$rl", titlefontsize=12))
savefig("img/rl.png")

llstm = LSTM(4,7)
display(plot(llstm, title="$llstm", titlefontsize=12))
savefig("img/llstm.png")

lgru = GRU(5,7)
display(plot(lgru, title="$lgru", titlefontsize=12))
savefig("img/lgru.png")

nnd = Chain(Dense(2,5), Dense(5,7,σ), Dense(7,2,relu),Dense(2,3))
display(plot(nnd, title="$nnd", titlefontsize=10, xaxis=nothing))
savefig("img/nnd.png")

nnr = Chain(Dense(2,5,σ), RNN(5,4,relu), LSTM(4,4), GRU(4,4), Dense(4,3))
display(plot(nnr, title="$nnr", titlefontsize=7))
savefig("img/nnr.png")

dx(x) = x[2:end]-x[1:end-1]
x³(x) = x.^3
nna = Chain(Dense(2,5,σ), dx, RNN(4,6,relu), x³, LSTM(6,4), GRU(4,4), Dense(4,3))
display(plot(nna, title="$nna", titlefontsize=7))
savefig("img/nna.png")

nnx = Chain(x³, dx, LSTM(5,10), Dense(10,5))
input_data = rand(6)
display(plot(nnx, input_data, title="$nnx", titlefontsize=9))
savefig("img/nnx.png")

nnrlwide = Chain(Dense(5,8), RNN(8,20), LSTM(20,10), Dense(10,7))
display(plot(nnrlwide, title="$nnrlwide", titlefontsize=9))
savefig("img/nnrlwide.png")

reshape6x1x1(a) = reshape(a, 6,  1, 1)
slice(a) = a[:,1,1]
nnrs = Chain(x³, Dense(3,6), reshape6x1x1, Conv((2,), 1=>1), slice, Dense(5,4))
display(plot(nnrs, Float32.(rand(3)), title="$nnrs", titlefontsize=9))
savefig("img/nnrs.png")

reshape3x3x1x1(a) = reshape(a, 3, 3, 1, 1)
nnrs2d = Chain(x³, Dense(4,9), reshape3x3x1x1, Conv((2,2), 1=>1), slice)
display(plot(nnrs2d, Float32.(rand(4)), title="$nnrs2d", titlefontsize=9))
savefig("img/nnrs2d.png")

#= nncg = Chain(Conv((3,3), 1=>8, leakyrelu, pad = 1),GroupNorm(8,4))
display(plot(nncg, Float32.(rand(6,6,1,1)), title="$nncg", titlefontsize=10))
savefig("img/nncg.png") =#

hdf5()
plot(nnr, title="$nnr with HDF5", titlefontsize=7)
Plots.hdf5plot_write("img/nnrhdf5.hdf5")
gr()
plthdf5_read = Plots.hdf5plot_read("img/nnrhdf5.hdf5")
display(plthdf5_read)

gr()
for t in themes
    theme(t)
    try
        display(plot(nnr, title="With theme $t", titlefontsize=10))
        savefig("img/nnr_$t.png")
    catch err
        println("Error in chain plot with theme $t: $err")
    end
end

theme(:default)

nothing


text:
using ApproxFun, DualNumbers, LinearAlgebra, Plots
    import ApproxFun: DualFun, jacobian

x=Fun(-π..π)
N = (T,u) -> [u(-π)-1;u(π)-1;u'(π)-1;u'' - u - sin(T*x/π)]

T0=Tk=0.
    u0=uk=one(x)

    (T,u)=(T0,u0)
    for k=1:4
        global T,u
        J1=jacobian.(N(DualFun(Fun(T)),u))
        J2=jacobian.(N(T,DualFun(u)))
        J=Operator[J1 J2]
        Tk,uk=J\N(T,u)
        @show Tk
        T,u=T-Number(Tk),u-uk
    end

plot(u)
norm(N(T,u)[end])  #9e-17

ext: jl

max_stars_repo_path: ODEs/Multivariate newton iteration.jl

max_issues_repo_path: ODEs/Multivariate newton iteration.jl

max_forks_repo_path: ODEs/Multivariate newton iteration.jl


text:
export request_guild_members,
    update_voice_status,
    update_status,
    update_presence

const conn_properties = Dict(
    "\$os"      => string(Sys.KERNEL),
    "\$browser" => "Discord.jl",
    "\$device"  => "Discord.jl",
)

const EMPTY = ErrorException("Empty")

const OPCODES = Dict(
    0 =>  :DISPATCH,
    1 =>  :HEARTBEAT,
    2 =>  :IDENTIFY,
    3 =>  :STATUS_UPDATE,
    4 =>  :VOICE_STATUS_UPDATE,
    6 =>  :RESUME,
    7 =>  :RECONNECT,
    8 =>  :REQUEST_GUILD_MEMBERS,
    9 =>  :INVALID_SESSION,
    10 => :HELLO,
    11 => :HEARTBEAT_ACK,
)

# Connection.

"""
    open(c::Client; delay::Period=Second(7))

Connect a [`Client`](@ref) to the Discord gateway.

The `delay` keyword is the time between shards connecting. It can be increased from its
default if you are frequently experiencing invalid sessions upon connection.
"""
function Base.open(c::Client; resume::Bool=false, delay::Period=Second(7))
    isopen(c) && throw(ArgumentError("Client is already connected"))
    c.ready = false

    # Clients can only identify once per 5 seconds.
    resume || sleep(c.shard * delay)

    logmsg(c, DEBUG, "Requesting gateway URL")
    resp = HTTP.get("$DISCORD_API/v$(c.version)/gateway")
    data = JSON.parse(String(resp.body))
    url = "$(data["url"])?v=$(c.version)&encoding=json"
    logmsg(c, DEBUG, "Connecting to gateway"; url=url)
    v = isdefined(c, :conn) ? c.conn.v + 1 : 1
    c.conn = Conn(opentrick(HTTP.WebSockets.open, url), v)

    logmsg(c, DEBUG, "receiving HELLO"; conn=c.conn.v)
    data, e = readjson(c.conn.io)
    e === nothing || throw(e)
    op = get(OPCODES, data["op"], data["op"])
    op === :HELLO || error("Expected opcode HELLO, received $op")
    hello(c, data)

    data = if resume
        Dict("op" => 6, "d" => Dict(
            "token" => c.token,
            "session_id" => c.state.session_id,
            "seq" => c.hb_seq,
        ))
    else
        d = Dict("op" => 2, "s" => c.hb_seq, "d" => Dict(
            "token" => c.token,
            "properties" => conn_properties,
        ))
        if !isempty(c.state.login_presence)
            d["d"]["presence"] = c.state.login_presence
        end
        if c.shards > 1
            d["d"]["shard"] = [c.shard, c.shards]
        end
        d
    end

    op = resume ? :RESUME : :IDENTIFY
    logmsg(c, DEBUG, "Writing $op"; conn=c.conn.v)
    writejson(c.conn.io, data) === nothing || error("Writing $op failed")

    logmsg(c, DEBUG, "Starting background maintenance tasks"; conn=c.conn.v)
    @async heartbeat_loop(c)
    @async read_loop(c)

    c.ready = true
end

"""
    isopen(c::Client) -> Bool

Determine whether the [`Client`](@ref) is connected to the gateway.
"""
Base.isopen(c::Client) = c.ready && isdefined(c, :conn) && isopen(c.conn.io)

"""
    close(c::Client)

Disconnect the [`Client`](@ref) from the Discord gateway.
"""
function Base.close(c::Client; statuscode::Int=1000)
    c.ready = false
    isdefined(c, :conn) && close(c.conn.io; statuscode=statuscode)
end

"""
    wait(c::Client)

Wait for an open [`Client`](@ref) to close.
"""
function Base.wait(c::Client)
    while isopen(c)
        wait(c.conn.io.cond)
        # This is an arbitrary amount of time to wait,
        # but we want to wait long enough to potentially reconnect.
        sleep(30)
    end
end

# Gateway commands.

"""
    request_guild_members(
        c::Client,
        guilds::Union{Integer, Vector{<:Integer};
        query::AbstractString="",
        limit::Int=0,
    ) -> Bool

Request offline guild members of one or more [`Guild`](@ref)s. [`GuildMembersChunk`](@ref)
events are sent by the gateway in response.
More details [here](https://discordapp.com/developers/docs/topics/gateway#request-guild-members).
"""
function request_guild_members(
    c::Client,
    guild::Integer;
    query::AbstractString="",
    limit::Int=0,
)
    return request_guild_members(c, [guild]; query=query, limit=limit)
end

function request_guild_members(
    c::Client,
    guilds::Vector{<:Integer};
    query::AbstractString="",
    limit::Int=0,
)
    throw_if_closed(c)

    return writejson(c.conn.io, Dict("op" => 8, "d" => Dict(
        "guild_id" => guilds,
        "query" => query,
        "limit" => limit,
    ))) === nothing
end

"""
    update_voice_state(
        c::Client,
        guild::Integer,
        channel::Union{Integer, Nothing},
        mute::Bool,
        deaf::Bool,
    ) -> Bool

Join, move, or disconnect from a voice channel. A [`VoiceStateUpdate`](@ref) event is sent
by the gateway in response.
More details [here](https://discordapp.com/developers/docs/topics/gateway#update-voice-state).
"""
function update_voice_state(
    c::Client,
    guild::Integer,
    channel::Union{Integer, Nothing},
    mute::Bool,
    deaf::Bool,
)
    throw_if_closed(c)

    return writejson(c.conn.io, Dict("op" => 4, "d" => Dict(
        "guild_id" => guild,
        "channel_id" => channel,
        "self_mute" => mute,
        "self_deaf" => deaf,
    ))) === nothing
end

"""
    update_status(
        c::Client,
        since::Union{Int, Nothing},
        activity::Union{Activity, Nothing},
        status::Union{PresenceStatus, AbstractString},
        afk::Bool,
    ) -> Bool

Indicate a presence or status update. A [`PresenceUpdate`](@ref) event is sent by the
gateway in response.
More details [here](https://discordapp.com/developers/docs/topics/gateway#update-status).
"""
function update_status(
    c::Client,
    since::Union{Int, Nothing},
    game::Union{Dict, NamedTuple, Activity, Nothing},
    status::Union{PresenceStatus, AbstractString},
    afk::Bool,
)
    throw_if_closed(c)

    return writejson(c.conn.io, Dict("op" => 3, "d" => Dict(
        "since" => since,
        "game" => game,
        "status" => status,
        "afk" => afk,
    ))) === nothing
end

# Client maintenance.

function heartbeat_loop(c::Client)
    v = c.conn.v
    try
        sleep(rand(1:round(Int, c.hb_interval / 1000)))
        heartbeat(c) || logmsg(c, ERROR, "Writing HEARTBEAT failed"; conn=v)

        while c.conn.v == v && isopen(c)
            sleep(c.hb_interval / 1000)
            if c.last_hb > c.last_ack && isopen(c) && c.conn.v == v
                logmsg(c, DEBUG, "Encountered zombie connection"; conn=v)
                reconnect(c)
            elseif !heartbeat(c) && c.conn.v == v && isopen(c)
                logmsg(c, ERROR, "Writing HEARTBEAT failed"; conn=v)
            end
        end
        logmsg(c, DEBUG, "Heartbeat loop exited"; conn=v)
    catch e
        logmsg(c, ERROR, "Heartbeat loop exited unexpectedly:\n$(catchmsg(e))"; conn=v)
    end
end

function read_loop(c::Client)
    v = c.conn.v
    try
        while c.conn.v == v && isopen(c)
            data, e = readjson(c.conn.io)
            if e !== nothing
                if e == EMPTY
                    continue
                elseif c.conn.v == v
                    handle_read_error(c, e)
                else
                    logmsg(c, DEBUG, "Read failed, but the connection is outdated"; conn=v, e=e)
                end
            elseif haskey(HANDLERS, data["op"])
                HANDLERS[data["op"]](c, data)
            else
                logmsg(c, WARN, "Unkown opcode"; op=data["op"])
            end
        end
        logmsg(c, DEBUG, "Read loop exited"; conn=v)
    catch e
        logmsg(c, ERROR, "Read loop exited unexpectedly:\n$(catchmsg(e))"; conn=v)
    end
end

# Event handlers.

function dispatch(c::Client, data::Dict)
    c.hb_seq = data["s"]

    T = get(EVENT_TYPES, data["t"], UnknownEvent)

    handlers = allhandlers(c, T)
    # If there are no handlers to call, don't bother parsing the event.
    isempty(handlers) && return

    evt = if T === UnknownEvent
        UnknownEvent(data)
    else
        val, e = tryparse(c, T, data["d"])
        if e === nothing
            val
        else
            T = UnknownEvent
            UnknownEvent(data)
        end
    end

    for (tag, handler) in handlers
        @async try
            handler.f(c, evt)
        catch e
            logmsg(c, ERROR, catchmsg(e); event=T, handler=tag)
            push!(c.state.errors, evt)
        finally
            # TODO: There are race conditions here.
            if handler.expiry isa Int && handler.expiry != -1
                handler.expiry -= 1
            end
        end
    end
end

function heartbeat(c::Client, ::Dict=Dict())
    e = writejson(c.conn.io, Dict("op" => 1, "d" => c.hb_seq))
    e === nothing && (c.last_hb = now())
    return e === nothing
end

function reconnect(c::Client, ::Dict=Dict(); resume::Bool=true)
    logmsg(c, INFO, "Reconnecting"; resume=resume)
    close(c; statuscode=resume ? 4000 : 1000)
    open(c; resume=resume)
end

function invalid_session(c::Client, data::Dict)
    logmsg(c, WARN, "Received INVALID_SESSION"; resumable=data["d"])
    sleep(rand(1:5))
    reconnect(c; resume=data["d"])
end

hello(c::Client, data::Dict) = c.hb_interval = data["d"]["heartbeat_interval"]

heartbeat_ack(c::Client, ::Dict) = c.last_ack = now()

# Gateway opcodes => handler function.
const HANDLERS = Dict(
    0   => dispatch,
    1   => heartbeat,
    7   => reconnect,
    9   => invalid_session,
    10  => hello,
    11  => heartbeat_ack,
)

# Error handling.

const CLOSE_CODES = Dict(
    1000 => :NORMAL,
    4000 => :UNKNOWN_ERROR,
    4001 => :UNKNOWN_OPCODE,
    4002 => :DECODE_ERROR,
    4003 => :NOT_AUTHENTICATED,
    4004 => :AUTHENTICATION_FAILED,
    4005 => :ALREADY_AUTHENTICATED,
    4007 => :INVALID_SEQ,
    4008 => :RATE_LIMITED,
    4009 => :SESSION_TIMEOUT,
    4010 => :INVALID_SHARD,
    4011 => :SHARDING_REQUIRED,
)

function handle_read_error(c::Client, e::Exception)
    logmsg(c, DEBUG, "Handling a $(typeof(e))"; e=e, conn=c.conn.v)
    c.ready || return
    if e isa HTTP.WebSockets.WebSocketError
        handle_close(c, e.status)
    else
        logmsg(c, ERROR, sprint(showerror, e))
        isopen(c) || reconnect(c)
    end
end

function handle_close(c::Client, status::Integer)
    err = get(CLOSE_CODES, status, :UNKNOWN_ERROR)
    if err === :NORMAL
        close(c)
    elseif err === :UNKNOWN_ERROR
        reconnect(c)
    elseif err === :UNKNOWN_OPCODE  # Probably a library bug.
        reconnect(c)
    elseif err === :DECODE_ERROR  # Probably a library bug.
        reconnect(c)
    elseif err === :NOT_AUTHENTICATED  # Probably a library bug.
        reconnect(c)
    elseif err === :AUTHENTICATION_FAILED
        logmsg(c, ERROR, "Authentication failed")
        close(c)
    elseif err === :ALREADY_AUTHENTICATED  # Probably a library bug.
        reconnect(c)
    elseif err === :INVALID_SEQ  # Probably a library bug.
        reconnect(c)
    elseif err === :RATE_LIMITED  # Probably a library bug.
        reconnect(c)
    elseif err === :SESSION_TIMEOUT
        reconnect(c)
    elseif err === :INVALID_SHARD
        logmsg(c, ERROR, "Invalid shard")
        close(c)
    elseif err === :SHARDING_REQUIRED
        logmsg(c, ERROR, "Sharding required")
        close(c)
    end
end

# Helpers.

function readjson(io)
    return try
        data = readavailable(io)
        if isempty(data)
            nothing, EMPTY
        else
            JSON.parse(String(data)), nothing
        end
    catch e
        nothing, e
    end
end

function writejson(io, body)
    return try
        write(io, json(body))
        nothing
    catch e
        e
    end
end

function throw_if_closed(c::Client)
    isopen(c) || throw(ArgumentError("Client is not connected"))
end


text:
#-------------------------------------------------------------------------------
# File: repair.jl
# Description: This files contains all function for the "repair" phase.
# Date: November 03, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

"""
    find_first_violation(solution::Solution, instance::Instance)

Return the index of the first car violating the paint_batch_limit constraint in the sequence,
return -1 if such a car does not exist
"""
function find_first_violation(solution::Solution, instance::Instance)
    first_violation = -1
    position = instance.nb_late_prec_day+1
    while position <= solution.length
        if solution.colors[position].width > instance.nb_paint_limitation
            first_violation = position + instance.nb_paint_limitation
            break
        end
        position += solution.colors[position].width
    end
    return first_violation
end

"""
    first_strategy_repair!(solution::Solution, instance::Instance)

First strategy
"""
function first_strategy_repair!(solution::Solution, instance::Instance)
    # First strategy
    b0 = instance.nb_late_prec_day+1

    position = b0
    counter = 1
    current_color = instance.color_code[solution.sequence[position]]
    position += 1
    first_violation = 0
    while position <= instance.nb_cars
        if instance.color_code[solution.sequence[position]] == current_color
            counter += 1
        else
            counter = 1
            current_color = instance.color_code[solution.sequence[position]]
        end
        if counter > instance.nb_paint_limitation
            if position >= b0
                if first_violation == 0
                    first_violation = position
                end
                car_RC_value = instance.RC_keys[solution.sequence[position]]
                car_pos = 1
                len = length(instance.same_RC[car_RC_value])
                while (car_pos <= len) && (instance.color_code[instance.same_RC[car_RC_value][car_pos]] == current_color)
                    car_pos += 1
                end
                if car_pos <= len
                    move_exchange!(solution, position, car_pos, instance)
                    counter = first_violation - position
                    if first_violation >= position
                        position = first_violation - 1
                    else
                        position -= 1
                    end
                    first_violation = 0
                elseif first_violation >= position
                    position -= 2
                    if counter == 2*instance.nb_paint_limitation
                        # This strategy can't repair
                        position = first_violation
                    end
                end
            else
                if first_violation != 0
                    position = first_violation
                else
                    # pcc are not possible in nb_late_prec_day so we keep the counter
                end
            end
        end
        position += 1
    end

    position = b0
    counter = 1
    current_color = instance.color_code[solution.sequence[position]]
    position += 1
    first_violation = 0
    while position <= instance.nb_cars
        if instance.color_code[solution.sequence[position]] == current_color
            counter += 1
        else
            counter = 1
            current_color = instance.color_code[solution.sequence[position]]
        end
        if counter > instance.nb_paint_limitation
            if position >= b0
                if first_violation == 0
                    first_violation = position
                end
                car_HPRC_value = instance.HPRC_keys[solution.sequence[position]]
                car_pos = 1
                len = length(instance.same_HPRC[car_HPRC_value])
                while (car_pos <= len) && (instance.color_code[instance.same_HPRC[car_HPRC_value][car_pos]] == current_color)
                    car_pos += 1
                end
                if car_pos <= len
                    move_exchange!(solution, position, car_pos, instance)
                    counter = first_violation - position
                    if first_violation >= position
                        position = first_violation - 1
                    else
                        position -= 1
                    end
                    first_violation = 0
                elseif first_violation >= position
                    position -= 2
                    if counter == 2*instance.nb_paint_limitation
                        # This strategy can't repair
                        position = first_violation
                    end
                end
            else
                if first_violation != 0
                    position = first_violation
                else
                    # pcc are not possible in nb_late_prec_day so we keep the counter
                end
            end
        end
        position += 1
    end
end

"""
    second_repair!(solution::Solution, instance::Instance)

Second strategy
"""
function second_strategy_repair!(solution::Solution, instance::Instance)
    # Second strategy
    first_violation = find_first_violation(solution, instance)
    b0 = instance.nb_late_prec_day+1
    while first_violation != -1
        # Compute the best_insertion index
        solution_value = weighted_sum(solution, instance)
        cost_insertion = zeros(instance.nb_cars)
        for i in b0:instance.nb_cars
            if i <= instance.nb_late_prec_day || i == first_violation # avoiding insterting in sequence from the precedent day and at the same position
                cost_insertion[i] = Inf
            else
                #checking for validity of insertion :  if it does not add an other violation of Paint batch constraint
                batch_size = 1
                batch_color = instance.color_code[solution.sequence[first_violation]]
                if instance.color_code[solution.sequence[i]] == batch_color # if inserting create a batch of length greater than one
                    batch_size += 1
                    j = 1
                    same_batch_sup = true
                    same_batch_inf = true
                    while same_batch_inf || same_batch_sup
                        if i-j > 0 && instance.color_code[solution.sequence[i-j]] == batch_color && same_batch_inf
                            batch_size += 1
                        else
                            same_batch_inf = false
                        end

                        if i+j <= instance.nb_cars && instance.color_code[solution.sequence[i+j]] == batch_color && same_batch_sup
                            batch_size += 1
                        else
                            same_batch_sup = false
                        end
                        j += 1
                    end
                end
                if batch_size > instance.nb_paint_limitation
                    cost_insertion[i] = Inf
                else
                    solution_copy = deepcopy(solution)
                    move_insertion!(solution, first_violation, i, instance)
                    cost_insertion[i] = weighted_sum(solution_copy, instance) - solution_value
                end
            end
        end
        best_insertion = argmin(cost_insertion)[1]
        move_insertion!(solution, first_violation, best_insertion, instance)

        first_violation = find_first_violation(solution, instance)
    end
end

"""
    repair!(solution::Solution, instance::Instance)

Apply 2 repair strategies on `solution`.
"""
function repair!(solution::Solution, instance::Instance)
    if solution.colors === nothing
        initialize_batches!(solution, instance)
    end

    first_strategy_repair!(solution, instance)
    second_strategy_repair!(solution, instance)

    @assert find_first_violation(solution, instance) == -1
    return solution
end


text:
# This is Julia code for discretizing AR(1) process by using Tauchen method (1986)
# written on 2017/10/17
# z' = rho*z+e e~N(0,(sig_e)^2)
#
# INPUTS
# N:     Number of grids
# m:     Max number of std devs from mean
# rho:   Coeffecient of AR(1) process z'=mew+rho*z+e, e~N(0,sigmasq)
# sig_e: Std of innovation, e, in AR(1) process
#
# OUTPUT
# vZ:    grids
# mPI:   transition prob. matrix

## Preliminary

###############################################################
## 0. Define CDF function (from https://www.johndcook.com/blog/cpp_phi/)
###############################################################

function phi(x::Real)
    # constants
    const a1 =  0.254829592;
    const a2 = -0.284496736;
    const a3 =  1.421413741;
    const a4 = -1.453152027;
    const a5 =  1.061405429;
    const p  =  0.3275911;

    # Save the sign of x
    sign = 1;
    if (x < 0)
        sign = -1;
    end
    x = abs(x)/sqrt(2.0);

    # A&S formula 7.1.26
    t = 1.0/(1.0 + p*x);
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*exp(-x*x);

    return 0.5*(1.0 + sign*y);
end

##############################################################
##############################################################
##                          MAIN PART
##############################################################
##############################################################

###########################################################
## 1. Define Parameters
###########################################################

# Ex. Standard quarterly calibration for productivity process
# of US economy

const N = 5;            # Number of grids, N
const m = 3;            # Max number of std devs from mean
const rho = 0.95;       # AR(1) coefficient
const sig_e = 0.007;    # std of innovation

###########################################################
## 2. Generating grids
###########################################################

sig_z = sqrt((sig_e^2)/(1-rho^2));  # std of AR(1) r.v.
z_max = m*sig_z;                    # Max grid
z_min = -z_max;                     # Min grid
d = (z_max-z_min) / (N-1);          # incriment of grids
vZ = linspace(z_min, z_max, N);     # Grid vector 1*N

###########################################################
## 3. Transition matrix
###########################################################

mPI = zeros(N,N); # Transition matrix N*N

for i in 1:N
    for j in 1:N
        if j == 1
            mPI[i,j] = phi((vZ[1]+d/2-rho*vZ[i])/sig_e);       # If z' is min grid
        elseif j == N-1
            mPI[i,j] = 1 - phi((vZ[N]-d/2-rho*vZ[i])/sig_e);   # z' is the max grid
        else
            mPI[i,j] = phi((vZ[j]+d/2-rho*vZ[i])/sig_e) - phi((vZ[j]-d/2-rho*vZ[i])/sig_e);
        end
    end
end

###########################################################
## 4. Display result with 4 decimals
###########################################################
print("Grid points are \n")
for i in 1:N
    @printf "%6.4f" vZ[i]
    print(" ")
end
print("\n")
print("Transition prob. matrix is\n")
for i in 1:N
    for j in 1:N
        @printf "%6.4f" mPI[i,j]
        print(" ")
    end
    print("\n")
end


text:
using DecayFit

irf = simulate_irf(0.047,256,44.7,100000,0.25,1.5)

dk = simulate_dk(0.097,256,44.7,irf,zeros(Int,size(irf)),
[0.1, 2.3, 0, 0, 0, 0],1,0)


text:
using Sockets
using Base.Threads

@enum Status reception gambling
@enum Game black_jack slot_machine roulette

mutex_dict = Dict()

mutable struct Player
    balance :: Int64
    status :: Status
    current_game :: Game
    socket :: TCPSocket
    dimension :: Int
    msg :: String
    token :: String
    diarrhea :: Bool
end

#julia + async + docker makes some problems...
#therefore we try use a mutex
function open_file_try(path, mode, max_trys=5)
    global mutex_dict
    if !haskey(mutex_dict, path)
        mutex_dict[path] = Mutex()
    end
    lock(mutex_dict[path])
    f = open(path, mode)
    unlock(mutex_dict[path])
    return f
end


text:
"""
Struct for Asian Floating Strike Option

		asOption=AsianFloatingStrikeOption(T::num1,isCall::Bool=true) where {num1 <: Number}

Where:\n
		T	=	Time to maturity of the Option.
		isCall  = true for CALL, false for PUT.
"""
struct AsianFloatingStrikeOption{num <: Number} <: AsianPayoff{num}
    T::num
    isCall::Bool
    function AsianFloatingStrikeOption(T::num, isCall::Bool = true) where {num <: Number}
        if T <= 0.0
            error("Time to Maturity must be positive")
        else
            return new{num}(T, isCall)
        end
    end
end

export AsianFloatingStrikeOption;

function payout(S::abstractArray, opt_::AsianFloatingStrikeOption) where {abstractArray <: AbstractArray{num_}} where {num_ <: Number}
    iscall = opt_.isCall ? 1 : -1
    zero_typed = zero(eltype(S))
    return max(iscall * (S[end] - mean(S)), zero_typed)
end


text:
## Save the small wikipedia graph
using SparseArrays
A = [0 1 1 1 1 0 1 0 0 0
     1 0 0 0 0 0 0 0 0 0
     0 0 0 0 1 1 1 0 0 0
     1 1 0 0 1 1 0 1 0 0
     1 1 1 1 0 0 0 0 1 0
     0 0 1 0 0 0 1 0 1 1
     0 0 0 0 0 1 0 0 0 1
     0 0 0 0 0 0 0 0 1 0
     0 0 0 0 0 0 0 1 0 0
     0 0 0 0 0 0 0 0 0 0]
A = sparse(A)
pages = ["PageRank", "Google", "Adjacency matrix",
"Markov chain", "Eigenvector", "Directed graph",
"Graph", "Linear system", "Vector space", "Multiset"]
ei,ej = findnz(A)
xy = [1.5 6.5;0 5.5; 0 2.25; -1.5 6.5; 0 4; -1.5 0.85; 1.5 0.85; -1 4.5; -1 3; 0 0.25]
##
using DelimitedFiles
writedlm("wiki-simple.edges", [ei ej])
writedlm("wiki-simple.xy", xy)
writedlm("wiki-simple.nodes", pages)
##
using GraphRecipes
#graphplot(ei,ej, x=xy[:,1], y=xy[:,2])
B = A.*A'
U = A - B # we need A' so that the ei, ej plots

graphplot(findnz(triu(B,1))[1:2]..., x = xy[:,1], y = xy[:,2],
  markercolor=nothing, markerstrokecolor=:white, curves=false,
  markersize=0,
  linecolor=1, linealpha=0.8, linewidth=0.7, arrow=:both, shorten=0.85,
  axis_buffer=0.5, background=nothing)

graphplot!(findnz(U)[1:2]..., x = xy[:,1], y = xy[:,2],
  markercolor=nothing, markerstrokecolor=:white, curves=false, curvature_scalar = -0.4,
  markersize=0,
   linecolor=2, linealpha=0.8, linewidth=0.7, arrow = :tail, shorten=0.85,
  axis_buffer=0.5, background=nothing)

xlims!(xlims().*1.2)
for i=1:10
  annotate!(xy[i,1],xy[i,2],text(pages[i],9,:white))
end
plot!(size=(300,300))
savefig("wiki-simple.pdf")

##
using GraphRecipes
#graphplot(ei,ej, x=xy[:,1], y=xy[:,2])
B = A.*A'
U = A - B # we need A' so that the ei, ej plots

graphplot(findnz(triu(B,1))[1:2]..., x = xy[:,1], y = xy[:,2],
  markercolor=nothing, markerstrokecolor=:white, curves=false,
  markersize=0, markerstrokewidth=0,
  linecolor=1, linealpha=0.8, linewidth=0.7, arrow=:both, shorten=0.85,
  axis_buffer=0.5, background=nothing)

graphplot!(findnz(U)[1:2]..., x = xy[:,1], y = xy[:,2],
  markercolor=nothing, markerstrokecolor=:white, curves=false, curvature_scalar = -0.4,
  markersize=0, markerstrokewidth=0,
   linecolor=2, linealpha=0.8, linewidth=0.7, arrow = :tail, shorten=0.85,
  axis_buffer=0.5, background=nothing)

xlims!(xlims().*1.2)
for i=1:10
  annotate!(xy[i,1],xy[i,2],text("$i",9,:white))
end
plot!(size=(300,300))
savefig("wiki-simple-numbers.pdf")

##
#= Tikz graph
\Vertex[L={PageRank},x=1.5,y=6.5]{1}
\Vertex[L={Google},x=0,y=5.5]{2}
\Vertex[L={Adjacency matrix},x=0,y=2.25]{3}
\Vertex[L={Markov chain},x=-1.5,y=6.5]{4}
\Vertex[L={Eigenvector},x=0,y=4]{5}
\Vertex[L={Directed graph},x=-1.5,y=0.85]{6}
\Vertex[L={Graph},x=1.5,y=0.85]{7}
\Vertex[L={Linear system},x=-1,y=4.5]{8}
\Vertex[L={Vector space},x=-1,y=3]{9}
\Vertex[L={Multiset},x=0,y=0.25]{10}

\Edge(4)(2)
\Edge(5)(2)
\Edge(4)(6)
\Edge(1)(7)
\Edge(3)(7)
\Edge(4)(8)
\Edge(5)(9)
\Edge(6)(9)
\Edge(6)(10)
\Edge(7)(10)
\Edge(2)(1)
\Edge(4)(1)
\Edge(5)(1)
\Edge(1)(2)
\Edge(5)(3)
\Edge(6)(3)
\Edge(1)(4)
\Edge(5)(4)
\Edge(1)(5)
\Edge(3)(5)
\Edge(4)(5)
\Edge(3)(6)
\Edge(7)(6)
\Edge(6)(7)
\Edge(9)(8)
\Edge(8)(9)
\Edge(1)(3)
=#


text:
@testset "life expectancy" begin

    @testset "parametric" begin

        # ALMC 2.5
        m = MortalityTables.Gompertz(a=0.0003,b=log(1.07))

        @test life_expectancy(m,0) ≈ 71.938 atol=0.001
        @test life_expectancy(m,50) ≈ 26.691 atol=0.001
    end

    @testset "vector" begin
        # if no distribution of deaths given, assume curtate otherwise complete
        t = MortalityTables.table("1980 CSO Basic Table – Male, ANB")

        # calculate sum of tpx in Excel
        @test life_expectancy(t.ultimate,55) ≈ 22.16469212 atol=1e-6
        @test life_expectancy(t.ultimate,100) ≈ 0.0 atol=1e-6

        @test_throws ArgumentError life_expectancy(t,100)

        # relation of curtate to complete, ALMC 2.6.1
        @test life_expectancy(t.ultimate,55,MortalityTables.Uniform()) ≈ 22.16469212 + 0.5 atol=1e-3
    end

end


text:
import HTTP.URIs: escapeuri

function post_on_slack_channel(msg, token, channel)
    try
        HTTP.get("https://slack.com/api/chat.postMessage?token=$token&channel=$channel&as_user=true&text=$(escapeuri(msg))")
    catch ex
        @info("Error while posting to slack channel: $(sprint(Base.showerror, ex, catch_backtrace()))")
    end
end


text:
function main(args)
    return args
end
# Load inline arguments for worker unit
num_workers, this = main(ARGS)
num_workers = parse(Int64, num_workers)
this = parse(Int64, this)

println("Worker " * string(this) * " of " * string(num_workers))
#
# Load packages and functions
include("../src/init.jl")
include("../src/fct.jl")
include("../src/worker_fct.jl")
data_type = "cal" # phy
dataset = "v07.01"

if data_type == "phy"
    plots_base_path = "../../waveforms/data/plots/raw_" * dataset * "/"
    base_path_raw   = "../../waveforms/data/raw_" * dataset * "/"
    base_path       = "../../waveforms/data/" * dataset * "/"
elseif data_type == "cal"
    plots_base_path = "../../waveforms/calib/plots/raw_" * dataset * "/"
    base_path_raw   = "../../waveforms/calib/raw_" * dataset * "/"
    base_path       = "../../waveforms/calib/" * dataset * "/"
end
number_of_pulses = 1e4
bl_range = 1:1:200

files = glob(base_path_raw * "*/*.h5");
files = get_share_for_worker(files, num_workers, this)

data_cal = Dict()
for ch in 0:1:36
    data_cal[string(ch)] = Dict()
    data_cal[string(ch)]["data"] = Table(energy = [], run = [], channel = [], AoEvetoed = [], datasetID = [], AoEclassifier = [], A = [], E = [], waveform = [])
    data_cal[string(ch)]["counter"] = 1
end

pro = Progress(length(files), dt=0.5,
                barglyphs=BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
                barlen=10);
for file in files
    if stat(file).size > 1000
        data = HDF5.h5open(file, "r") do h5f
            LegendHDF5IO.readdata(h5f, "data")
        end;
        for ch in unique(data.channel)
            ch_str = lpad(ch, 2, "0");
            temp = data |> @filter(_.channel == ch) |> Table
#             cal_factor = cal[string(ch)]["value"]
            waveforms = []
            A = []
            E = []
            for wf in temp.waveform
                pulse   = -1 .*float.(wf.value)
                pulse .-= sum(pulse[bl_range])/length(bl_range)
#                 pulse .*= cal_factor
                push!(A, maximum(multi_mwa(diff(pulse),5,3)))
                push!(E, maximum(multi_mwa(diff(pulse),201,6)))
                push!(waveforms, RDWaveform(wf.time, pulse))
            end

            append!(data_cal[string(ch)]["data"], Table(energy   = temp.energy,
                                                        run      = temp.run,
                                                        channel  = temp.channel,
                                                        AoEvetoed= temp.AEvetoed,
                                                        datasetID= temp.datasetID,
                                                        AoEclassifier = temp.AEclassifier,
                                                        A        = A,
                                                        E        = E,
                                                        waveform = waveforms))
            if size(data_cal[string(ch)]["data"],1) >= number_of_pulses
                filename = base_path *  ch_str * "-" * channel_to_name[ch] * "/"
                !isdir(filename) ? mkpath(filename) : "path exists"
                filename *= lpad(data_cal[string(ch)]["counter"], 4, "0") * "-" * string(this) * ".h5"
                data_cal[string(ch)]["counter"] += 1
                HDF5.h5open(filename, "w") do h5f
                    LegendHDF5IO.writedata( h5f, "data", Table( energy       = float.(data_cal[string(ch)]["data"].energy),
                                                                run          = Int.(data_cal[string(ch)]["data"].run),
                                                                channel      = Int.(data_cal[string(ch)]["data"].channel),
                                                                AoEvetoed    = Int.(data_cal[string(ch)]["data"].AoEvetoed),
                                                                datasetID    = Int.(data_cal[string(ch)]["data"].datasetID),
                                                                AoEclassifier= float.(data_cal[string(ch)]["data"].AoEclassifier),
                                                                A            = float.(data_cal[string(ch)]["data"].A),
                                                                E            = float.(data_cal[string(ch)]["data"].E),
                                                                waveform     = StructArray{RDWaveform}(Array{RDWaveform,1}(data_cal[string(ch)]["data"].waveform))))
                end
                data_cal[string(ch)]["data"] = Table(energy = [], run = [], channel = [], AoEvetoed = [], datasetID = [], AoEclassifier = [], A = [], E = [], waveform = [])
            end
        end
    end
    next!(pro)
end
for ch in 0:1:36
    ch_str = lpad(ch, 2, "0");
    if size(data_cal[string(ch)]["data"],1) > 0
        filename = base_path *  ch_str * "-" * channel_to_name[ch] * "/"
        !isdir(filename) ? mkpath(filename) : "path exists"
        filename *= lpad(data_cal[string(ch)]["counter"], 4, "0") * "-" * string(this) * ".h5"
        data_cal[string(ch)]["counter"] += 1
        HDF5.h5open(filename, "w") do h5f
            LegendHDF5IO.writedata( h5f, "data", Table( energy       = float.(data_cal[string(ch)]["data"].energy),
                                                        run          = Int.(data_cal[string(ch)]["data"].run),
                                                        channel      = Int.(data_cal[string(ch)]["data"].channel),
                                                        AoEvetoed     = Int.(data_cal[string(ch)]["data"].AoEvetoed),
                                                        datasetID    = float.(data_cal[string(ch)]["data"].datasetID),
                                                        AoEclassifier= float.(data_cal[string(ch)]["data"].AoEclassifier),
                                                        A            = float.(data_cal[string(ch)]["data"].A),
                                                        E            = float.(data_cal[string(ch)]["data"].E),
                                                        waveform     = StructArray{RDWaveform}(Array{RDWaveform,1}(data_cal[string(ch)]["data"].waveform))))
        end
    end
end

max_stars_repo_path: calib/02 - Data Calibration A and E determination.jl


text:
include("define_constants.jl")
# Include linear algebara:
if VERSION >= v"0.7"
  using LinearAlgebra
end

# Include definition of Transit structure type:
include("transit_structure.jl")
# Include code which computes linear limb-darkening term:
include("s2.jl")
# Include code which computes M_n and N_n:
include("Mn_compute.jl")
include("Nn_compute.jl")

"""
    compute_uniform(t)

Computes the coefficient for the uniform disk case, `s_0`, given
a `TransitStruct` instance `t`.
"""
function compute_uniform!(t::Transit_Struct{T}) where {T <: Real}
  r=t.r; b=t.b; r2=r*r; b2=b*b
  # Uniform disk case:
  # Compute sn[1] and its derivatives:
  if b <= 1-r  # k^2 > 1
    t.sn[1] = pi*(1-r2)
    if t.grad
      t.dsndr[1] = -2*pi*r
      t.dsndb[1] = 0.
    end
    t.kap0 = convert(T,pi); t.kck = zero(T)
  else
    # Twice area of kite-shaped region connecting centers of circles & intersection points:
  #  t.kite_area2 = sqrt(sqarea_triangle(one(T),b,r))
    t.kite_area2 = sqrt(t.sqarea)
    # Angle of section for occultor:
    t.kap0  = atan(t.kite_area2,(r-1)*(r+1)+b2)
    # Angle of section for source:
    t.pimkap1 = atan(t.kite_area2,(r-1)*(r+1)-b2)
    # Flux of visible uniform disk:
    t.sn[1] = t.pimkap1 - r2*t.kap0 + .5*t.kite_area2
    if t.grad
      t.dsndr[1]= -2*r*t.kap0
      t.dsndb[1]= t.kite_area2/b
    end
    t.kck = t.kite_area2*t.fourbrinv
  end
  return
end

"""
    sqarea_triangle(a,b,c)

Function which computes sixteen times the square of the area
of a triangle with sides a, b and c using Kahan method.
How to compute (sixteen times the) area squared of triangle with
high precision (Goldberg 1991).
"""
function sqarea_triangle(a::T,b::T,c::T) where {T <: Real}

  # First, do a quick sort of three numbers:
  if c > b
    tmp = c; c = b; b = tmp
  end
  if c > a
    tmp = a; a = c; c = tmp
  end
  if b > a
    tmp = b; b = a; a = tmp
  end
  area = (a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c))
  return area

end

"""
    transit_poly_g(t)

Given a `TransitStruct` instance `t`, computes a limb-darkened transit light
curve without gradient.

"""
function transit_poly_g(t::Transit_Struct{T}) where {T <: Real}
  # Number of limb-darkening components to include (beyond 0 and 1):
  # We are parameterizing these with the function:
  # g_n = g_n [(n+2) z^n - n z^{n-2}] for n >= 2
  # while g_{0} = g_0 z^0 (uniform source) and g_{1} = g_1 z^1 (linear limb-darkening)
  # which gives a Green's integral of:
  # P(G_n) = \int_{\pi-\phi}^{2\pi+\phi} (1-r^2-b^2-2br s_\varphi)^{n/2} (r+b s_\varphi) d\varphi
  # for which we have a solution in terms of M_n.
  # The variable "t" is a structure which contains transit parameters
  # and intermediate quantities computed from these:
  r=t.r; b=t.b; n = t.n; r2=r*r; b2 = b*b
  # Set up a vector for storing results of P(G_n)-Q(G_n); note that
  # this is a different vector than the Starry case:

  # Check for different cases:
  if b >= 1+r
    # unobscured - return one:
    return one(T)
  end
  if r >= 1+b
    # full obscuration - return zero:
    return zero(T)
  end
  if b == 0.0
    # Annular eclipse - integrate around the full boundary of both bodies:
    onemr2 = 1-r2
    flux = zero(T); t.sqrt1mr2 = sqrt(onemr2)
    flux = (t.g_n[1]*onemr2+t.twothird*t.g_n[2]*t.sqrt1mr2^3)
    fac= 2r2*onemr2
    @inbounds for i=2:t.n
      flux += -t.g_n[i+1]*fac
      fac *= t.sqrt1mr2
    end
    return flux*pi*t.den
  else
  # Next, compute k^2 = m:
    t.onembmr2=(r+1-b)*(1-r+b); t.fourbr = 4b*r; t.fourbrinv = inv(t.fourbr)
    t.sqbr = sqrt(b*r); t.sqbrinv = inv(t.sqbr)
    t.onembmr2inv=inv(t.onembmr2); t.sqonembmr2 = sqrt(t.onembmr2)
    t.onembpr2 = (1-r-b)*(1+b+r)
    t.sqarea = sqarea_triangle(one(T),r,b)
  #  t.k2 = t.onembmr2*t.fourbrinv
    t.k2 = t.onembpr2*t.fourbrinv+1
    t.onemr2mb2 = (1.0-r)*(1.0+r)-b2
    if t.k2 > 0
      t.k = sqrt(t.k2)
    else
      println("negative k2: ",t.k2," r: ",r," b: ",b)
      t.k2 = 0.0; t.k = 0.0
    end
    if t.k2 > 1
      if t.k2 > 2.0
        t.kc2 = 1.0-inv(t.k2)
        t.kc = sqrt(t.kc2)
      else
        t.kc2 = t.onembpr2/((1-b+r)*(1-r+b))
        t.kc = sqrt(t.kc2)
      end
    else
      if t.k2 > 0.5
        t.kc2 = (r-1+b)*(b+r+1)*t.fourbrinv
        t.kc = sqrt(t.kc2)
      else
        t.kc2 = 1.0-t.k2
        t.kc = sqrt(t.kc2)
      end
    end
  end

  # Compute uniform case, sn[1]:
  compute_uniform!(t)
  # If uniform model, then return:
  if t.n == 0
    return t.g_n[1]*t.sn[1]*t.den
  end

  # Compute linear case, sn[2]:
  #s2!(t)
  t.sn[2],t.Eofk,t.Em1mKdm = s2_ell(r,b)
  #sn2,Eofk,Em1mKdm = s2_ell(r,b)
  # Now, compare results:
  #if abs(t.sn[2]-sn2) > 1e-8*abs(sn2)
  #  println("r: ",r," b: ",b," sn[2]: ",t.sn[2]," sn2: ",sn2)
  #end
  # Special case of linear limb-darkening:
  if t.n == 1
    flux = t.g_n[1]*t.sn[1]+t.g_n[2]*t.sn[2]
    flux *= t.den  # for g_2 and above, the flux is zero.
    return flux
  end

  if t.n >= 2
  # Special case of quadratic limb-darkening:
  # Transformed expressions from Mandel & Agol:
    eta2 = r2*(r2+2*b2)
    if t.k2 >= 1
      four_pi_eta = 2pi*(eta2-1.0)
    else
      four_pi_eta = 2*(-t.pimkap1+eta2*t.kap0-0.25*t.kite_area2*(1.0+5r2+b2))
    end
    t.sn[3] = 2*t.sn[1]+four_pi_eta
    if t.n == 2
      flux = t.g_n[1]*t.sn[1]+t.g_n[2]*t.sn[2]+t.g_n[3]*t.sn[3]
      flux *= t.den  # for g_2 and above, the flux is zero.
      return flux
    end
  end

  # Compute the M_n functions:
  if t.k2 > 0
  #  if (t.k2 < 0.5 || t.k2 > 20.0) && t.n > 3
    if (t.k2 < 0.5) && t.n > 3
  # This computes M_n for the largest four m, and then works down to smaller values:
      Mn_lower!(t)
    else
  # This computes Mn for m=0 to 3, and then works upward to larger m:
      Mn_raise!(t)
    end
  end

  # Add up first three terms in flux numerator:
  flux = t.g_n[1]*t.sn[1]+t.g_n[2]*t.sn[2]+t.g_n[3]*t.sn[3]
  # Next, loop over the Green's function components:
  @inbounds for n=3:t.n
  #  pofgn_M = (1+(r-b)*(r+b))*t.Mn[n+1]-t.Mn[n+3]
  #  pofgn_M = t.onemr2mb2*t.Mn[n+1]-t.Mn[n+3]
  #  pofgn_M = 2*r^2*t.Mn[n+1]-n/(n+2)*((1-r^2-b^2)*t.Mn[n+1]+(1-(b-r)^2)*((b+r)^2-1)*t.Mn[n-1])
  #  pofgn_M = 2*r^2*t.Mn[n+1]-n/(n+2)*((1-r^2-b^2)*t.Mn[n+1]+sqarea_triangle(one(T),r,b)*t.Mn[n-1])
  #  pofgn_M = 2*r^2*t.Mn[n+1]-n/(n+2)*((1-r^2-b^2)*t.Mn[n+1]+t.sqarea*t.Mn[n-1])
  #  pofgn_M = 2*r2*t.Mn[n+1]-n/(n+2)*(t.onemr2mb2*t.Mn[n+1]+t.sqarea*t.Mn[n-1])
    pofgn_M = 2*r2*t.Mn[n+1]-n*t.ninv[n+2]*(t.onemr2mb2*t.Mn[n+1]+t.sqarea*t.Mn[n-1])
  # Q(G_n) is zero in this case since on limb of star z^n = 0 at the stellar
  # boundary for n > 0.
  # Compute sn[n]:
    t.sn[n+1] = -pofgn_M
      flux += t.g_n[n+1]*t.sn[n+1]
  end
  # That's it!
  # flux = sum(t.g_n.*t.sn)/(pi*(t.g_n[1]+t.twothird*t.g_n[2]))  # for g_2 and above, the flux is zero.
  # Divide by denominator:
  flux *= t.den  # for g_2 and above, the flux is zero.
  return flux
end

"""
    transit_poly!(r,b,u_n,dfdrbu)

Given a radius-ratio, impact parameter, vector of limb-darkening coefficients
of size N, and pre-allocated vector for derivatives of size N+2, returns
the flux (normalized to one for unocculted star) for a limb-darkened transit
light curve and its gradient with respect to the input parameters is returned.

# Arguments
- `r::Real`: The radius of the occultor in units of the radius of the occulted body.
- `b::Real`: The (initial) impact parameter of the occultation.
- `u_n::Array{Real,1}`: The array of limb darkening coefficients.  Size N.
- `dfdrbu::Array{Real,1}`: Gradient of the flux with respect to all input parameters. Size is N+2
"""
function transit_poly!(r::T,b::T,u_n::Array{T,1},dfdrbu::Array{T,1}) where {T <: Real}
  # Initialize `TransitStruct`:
  t = transit_init(r,b,u_n,true)
  # Pass g_n (without last two dummy values):
  flux = transit_poly_g!(t)
  # Now, transform derivaties from g to u:
  fill!(dfdrbu,zero(T))
  # u_n derivatives:
  dfdrbu[1] = t.dfdrb[1]
  dfdrbu[2] = t.dfdrb[2]
  @inbounds for i=1:t.n, j=0:t.n
    dfdrbu[i+2] += t.dfdg[j+1]*t.dgdu[j+1,i]
  end
  return flux
end

"""
    transit_poly(r,b,u_n)

Given a radius-ratio, impact parameter, vector of limb-darkening coefficients
of size N, returns the flux (normalized to one for unocculted star) for a
limb-darkened transit.

# Arguments
- `r::Real`: The radius of the occultor in units of the radius of the occulted body.
- `b::Real`: The (initial) impact parameter of the occultation.
- `u_n::Array{Real,1}`: The array of limb darkening coefficients.  Size N.
"""
function transit_poly(r::T,b::T,u_n::Array{T,1}) where {T <: Real}
  t = transit_init(r,b,u_n,false)
  # Pass g_n (without last two dummy values):
  return transit_poly_g(t)
end

"""
    transit_poly!(t)

Given a `TransitStruct` instance `t`, with Green's coefficients already initialized,
computes a limb-darkened transit light curve with optional gradient with
respect to r, b and u coefficients.

"""
function transit_poly!(t::Transit_Struct{T}) where {T <: Real}
  # This function assumes that g_n has already been computed from u_n
  # (this can be used to save compute time when limb-darkening is fixed for
  # a range of radii/impact parameters).
  # Pass transit structure, and compute flux:
  if t.grad
    flux = transit_poly_g!(t)
    # Now, transform derivaties from g to u:
    fill!(t.dfdu,zero(T))
    BLAS.gemv!('T',1.0,t.dgdu,t.dfdg,0.0,t.dfdu)
  #  @inbounds for i=1:t.n, j=0:t.n
  #    t.dfdu[i] += t.dfdg[j+1]*t.dgdu[j+1,i]
  #  end
    return flux
  else
    return transit_poly_g(t)
  end
  return
end

"""
    transit_poly_g!(t)

Given a `TransitStruct` instance `t`, with Green's coefficients already initialized,
computes a limb-darkened transit light curve including gradient with
respect to r, b and g coefficients.

"""
function transit_poly_g!(t::Transit_Struct{T}) where {T <: Real}
  r = t.r; b=t.b; n = t.n; r2 =r*r; b2=b*b; bcut = 1e-3
  @assert((length(t.g_n)) == length(t.dfdg))
  @assert(r > 0)
  # Number of limb-darkening components to include (beyond 0 and 1):
  # We are parameterizing these with the function:
  # g_n = g_n [(n+2) z^n - n z^{n-2}] for n >= 2
  # while g_{0} = g_0 z^0 (uniform source) and g_{1} = g_1 z^1 (linear limb-darkening)
  # which gives a Green's integral of:
  # P(G_n) = \int_{\pi-\phi}^{2\pi+\phi} (1-r^2-b^2-2br s_\varphi)^{n/2} (r+b s_\varphi) d\varphi
  # for which we have a solution in terms of M_n.
  # Compute the derivative of the flux with respect to the different coefficients.

  # Set up a vector for storing results of P(G_n)-Q(G_n); note that
  # this is a different vector than the Starry case:
  # Check for different cases:
  if b >= 1+r || r ==  0.0
    # unobscured - return one, and zero derivatives:
    fill!(t.dfdrb,zero(T))
    fill!(t.dfdg,zero(T))
    return one(T)
  end
  if r >= 1+b
    # full obscuration - return zero, and zero derivatives:
    fill!(t.dfdrb,zero(T))
    fill!(t.dfdg,zero(T))
    return zero(T)
  end
  if b == 0.0
    # Annular eclipse - integrate around the full boundary of both bodies:
    flux = zero(T); onemr2 = 1-r2; t.sqrt1mr2 = sqrt(onemr2)
    fill!(t.dfdg,zero(T))
    flux = (t.g_n[1]*onemr2+t.twothird*t.g_n[2]*t.sqrt1mr2^3)*pi*t.den
    fac  = 2r2*onemr2*pi*t.den
    facd = -2r*pi*t.den
    t.dfdrb[1] = t.g_n[1]*facd + t.g_n[2]*facd*t.sqrt1mr2
    @inbounds for i=2:t.n
      flux -= t.g_n[i+1]*fac
      t.dfdrb[1] += t.g_n[i+1]*facd*(2*onemr2-i*r2)
      t.dfdg[i+1] -= fac
      fac *= t.sqrt1mr2
      facd *= t.sqrt1mr2
    end
    #  dfdrb[2]=0 since the derivative with respect to b is zero.
    t.dfdg[1] = (onemr2-flux)*pi*t.den
    t.dfdg[2] = t.twothird*(t.sqrt1mr2^3-flux)*pi*t.den
    return flux
  else
  # Next, compute k^2 = m:
    t.onembmr2=(r-b+1)*(1-r+b); t.fourbr = 4b*r; t.fourbrinv = inv(t.fourbr)
    t.sqbr = sqrt(b*r); t.sqbrinv = inv(t.sqbr)
    t.onembmr2inv=inv(t.onembmr2); t.sqonembmr2 = sqrt(t.onembmr2)
    t.onembpr2 = (1-r-b)*(1+b+r)
    t.sqarea = sqarea_triangle(one(T),r,b)
    t.k2 = t.onembmr2*t.fourbrinv;
    t.onemr2mb2 = (1.0-r)*(1.0+r)-b2
    if t.k2 > 0
      t.k = sqrt(t.k2)
    else
      println("negative k2: ",t.k2," r: ",r," b: ",b)
      t.k2 = 0.0; t.k = 0.0
    end
    if t.k2 >= 1
      if t.k2 > 2.0
        t.kc2 = 1.0-inv(t.k2)
        t.kc = sqrt(abs(t.kc2))
      else
        t.kc2 = t.onembpr2/((1+r-b)*(1-r+b))
        t.kc = sqrt(abs(t.kc2))
      end
    else
      if t.k2 > 0.5
        t.kc2 = (r-1+b)*(b+r+1)*t.fourbrinv
        t.kc = sqrt(abs(t.kc2))
      else
        t.kc2 = (r-1+b)*(b+r+1)*t.fourbrinv
        t.kc = sqrt(abs(t.kc2))
      end
    end
  end

  # Compute uniform case:
  compute_uniform!(t)

  # Compute sn[2] and its derivatives:
  #if t.n >= 1
  #  s2!(t)
  #  # debug:
  #  s2_grad = zeros(T,2)
  #  sn2,Eofk,Em1mKdm = s2!(r,b,s2_grad)
  #  if (abs(t.sn[2]-sn2) > 1e-8) || (abs(t.s2_grad[1]-s2_grad[1]) > 1e-8) || (abs(t.s2_grad[2]-s2_grad[2]) > 1e-8)
  #    println("r: ",r,"b: ",b," s2 new: ",t.sn[2]," s2 old: ",sn2)
  #    println("ds2/dr new: ",t.s2_grad[1]," ds2/dr old: ",s2_grad[1])
  #    println("ds2/db new: ",t.s2_grad[2]," ds2/db old: ",s2_grad[2])
  #    read(STDIN,Char)
  #  end
  #end
  if t.n >= 1
    t.sn[2],t.Eofk,t.Em1mKdm = s2!(r,b,t.s2_grad)
    t.dsndr[2] = t.s2_grad[1]
    t.dsndb[2] = t.s2_grad[2]
  end

  # Special case of quadratic limb-darkening:
  if t.n >= 2
  #  if t.n == 2
  # Transformed expressions from Mandel & Agol for n=2:
    r2pb2 = r2+b2
    eta2 = r2*(r2pb2+b2)
    deta2dr =  2*r*r2pb2
    deta2db = 2*b*r2
    if t.k2 > 1
      four_pi_eta = 2pi*(eta2-1.0)
      detadr = 4pi*deta2dr
      detadb = 4pi*deta2db
    else
      four_pi_eta = 2*(-t.pimkap1+eta2*t.kap0-0.25*t.kite_area2*(1.0+5r2+b2))
      detadr = 8r*(r2pb2*t.kap0-t.kite_area2)
      detadb = 2/b*(4*b2*r2*t.kap0-(1+r2pb2)*t.kite_area2)
    end
    t.sn[3] = 2*t.sn[1]+four_pi_eta
    t.dsndr[3] = 2*t.dsndr[1]+detadr
    t.dsndb[3] = 2*t.dsndb[1]+detadb
  #  else
    if t.n > 2
      # Compute the M_n functions:
      if t.k2 > 0
  #      if (t.k2 < 0.5 || t.k2 > 20.0) && t.n > 3
        if (t.k2 < 0.5) && t.n > 3
      # This computes Mn for largest four m, and then works down to smaller values:
          Mn_lower!(t)
          if b < bcut
            Nn_lower!(t)
          end
        else
      # This computes Mn and then works upward to larger m:
          Mn_raise!(t)
          if b < bcut
            Nn_raise!(t)
          end
        end
      end
      # Next, loop over the Green's function components:
      binv = inv(b)
  #    @inbounds for n=2:t.n
      @inbounds for n=3:t.n
  #      pofgn_M = (1+(r-b)*(r+b))*t.Mn[n+1]-t.Mn[n+3]
  #      pofgn_M = t.onemr2mb2*t.Mn[n+1]-t.Mn[n+3]
  #      pofgn_M = 2*r^2*t.Mn[n+1]-n/(n+2)*((1-r2-b2)*t.Mn[n+1]+t.kite_area2^2*t.Mn[n-1])
  #      pofgn_M = 2*r^2*t.Mn[n+1]-n/(n+2)*(t.onemr2mb2*t.Mn[n+1]+sqarea_triangle(one(T),r,b)*t.Mn[n-1])
  #      pofgn_M = 2*r2*t.Mn[n+1]-n/(n+2)*(t.onemr2mb2*t.Mn[n+1]+t.sqarea*t.Mn[n-1])
        pofgn_M = 2*r2*t.Mn[n+1]-n*t.ninv[n+2]*(t.onemr2mb2*t.Mn[n+1]+t.sqarea*t.Mn[n-1])
      # Q(G_n) is zero in this case since on limb of star z^n = 0 at the stellar
      # boundary for n > 0.
      # Compute sn[n]:
        t.sn[n+1] = -pofgn_M
        dpdr_M = 2*r*((n+2)*t.Mn[n+1]-n*t.Mn[n-1])
  #      dpdr_M = 2*r*((2-(n+2)*(b-r)^2)*t.Mn[n-1]-4*(n+2)*b*r*t.Nn[n-1])
  #      dpdr_M = 2*r*((2+n-n/t.onembmr2)*t.Mn[n+1]-n/t.k2*t.Nn[n+1])
        t.dsndr[n+1] = -dpdr_M
        if b < bcut
        # For small b, use function which doesn't involve division by b:
          dpdb_M = n*(t.Mn[n-1]*(2*r^3+b^3-b-3*r2*b)+b*t.Mn[n+1]-4*r^3*t.Nn[n-1])
        else
          dpdb_M = n*binv*((t.Mn[n+1]-t.Mn[n-1])*(r2+b2)+(b2-r2)^2*t.Mn[n-1])
        end
        t.dsndb[n+1] = -dpdb_M
      end
    end
  end
  # That's it!
  # Compute derivatives with respect to the coefficients:
  flux = zero(T)
  t.dfdrb[1]=zero(T)  # Derivative with respect to r
  t.dfdrb[2]=zero(T)  # Derivative with respect to b
  @inbounds for i=0:t.n
    # derivatives with respect to the coefficients:
    t.dfdg[i+1]= t.sn[i+1]*t.den
    # total flux:
    flux += t.g_n[i+1]*t.dfdg[i+1]
    # derivatives with respect to r and b:
    t.dfdrb[1] += t.g_n[i+1]*t.dsndr[i+1]*t.den
    t.dfdrb[2] += t.g_n[i+1]*t.dsndb[i+1]*t.den
  end
  # Include derivatives with respect to first two g_n parameters:
  t.dfdg[1] -= flux*t.den*pi
  t.dfdg[2] -= flux*t.den*pi*t.twothird
  return flux
end

max_stars_repo_name: dfm/limbdark

max_issues_repo_name: dfm/limbdark

max_forks_repo_name: dfm/limbdark


text:
using Documenter, PhyloNetworks

makedocs()

# Versions of `mkdocs` and `mkdocs-material` specified manually to avoid conflicts.
# See here: https://discourse.julialang.org/t/mkdocs-material-in-documenter/13764/3
# To be kept in mind: those versions might evolve in the future.
deploydocs(
    deps   = Deps.pip("pygments", "mkdocs==0.17.5", "mkdocs-material==2.9.4", "python-markdown-math"),
    repo = "github.com/crsl4/PhyloNetworks.jl.git",
    julia  = "0.6",
    osname = "linux"
)


text:
using Flux, Statistics, BSON, Dates, Random, Images, FileIO, JLD2
using Flux: onehotbatch, onecold, crossentropy, throttle, @epochs
using Base.Iterators: partition, take
using CuArrays
using CUDAnative: tanh, log, exp

@load "./train.jld2" train_imgs train_labels
@load "./test.jld2" test_imgs test_labels

batch_size = 128
test_set = ((gpu(cat((test_imgs[batch]/255f0)..., dims=4)), gpu(onehotbatch(test_labels[batch],1:3755))) for batch in partition(1:length(test_imgs), batch_size))

function conv_unit(chanel, nb_filters, mp=false)
  conv_bn = Chain(Conv((3, 3), chanel => nb_filters, leakyrelu, pad=1, stride=1),BatchNorm(nb_filters))
  mp ? Chain(conv_bn..., x -> maxpool(x, (3, 3), stride=2, pad=1)) : conv_bn
end

model = Chain(
  # input size 64 x 64 x 1 x batch_szie (WHCN order)
  Conv((3, 3), 1 => 64, leakyrelu, pad=1, stride=2),
  BatchNorm(64),
  conv_unit(64, 128)...,
  conv_unit(128, 128, true)...,
  conv_unit(128, 256)...,
  conv_unit(256, 256, true)...,
  conv_unit(256, 384)...,
  conv_unit(384, 384)...,
  conv_unit(384, 384, true)...,
  conv_unit(384, 512)...,
  conv_unit(512, 512)...,
  conv_unit(512, 512, true)...,
  x -> reshape(x, :, size(x, 4)),
  Dense(2*2*512, 1024, leakyrelu),
  BatchNorm(1024),
  Dense(1024, 256, leakyrelu),
  BatchNorm(256),
  Dense(256, 3755),
  softmax) |> gpu

loss(x, y) = crossentropy(model(x), y)
predict = mapleaves(Flux.Tracker.data, model)
accuracy(x, y) = mean(cpu(onecold(predict(x))) .== onecold(cpu(y)))
opt = ADAM()
# uncomment to apply inverse time decay
# Flux.Optimise.Optimiser(InvDecay(), opt)

# callback function to show accuracy and save model every 600s
evalcb = throttle(600) do
  Flux.testmode!(predict)
  acc = mean(accuracy(x...) for x in test_set)
  @show(acc)
  model_cpu = cpu(model)
  rightnow = (Dates.format(now(), "yyyy-mm-ddTHH-MM-SS"))
  BSON.@save "model-$(rightnow)-A$(round(UInt16,10000acc)).bson" model_cpu acc
end

@epochs 2 begin
  train_set = ((gpu(cat((train_imgs[batch]/255f0)..., dims=4)), gpu(onehotbatch(train_labels[batch],1:3755))) for batch in partition(randperm(length(train_imgs)), batch_size))
  Flux.train!(loss, params(model), train_set, opt, cb=evalcb)
end

max_stars_repo_name: afternone/HCCR

max_issues_repo_name: afternone/HCCR

max_forks_repo_name: afternone/HCCR


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

using Markdown
using InteractiveUtils

# ╔═╡ c52e11fe-83fb-11eb-29b7-77aeaf7ebc83
using Images,TestImages

# ╔═╡ 8228eedc-8416-11eb-0daa-81aacd64af55
begin
using Base: OneTo, tail
using OffsetArrays

struct PaddedView{T,N,I,A} <: AbstractArray{T,N}
    fillvalue::T
    data::A
    indices::I
    function PaddedView{T,N,I,A}(fillvalue,data,indices::NTuple{N,AbstractUnitRange}) where {T,N,I,A}
        new{T,N,I,A}(fillvalue, data, indices)
    end
end

function PaddedView(fillvalue::FT, data::AbstractArray{T}, args...) where {FT, T}
    PaddedView{filltype(FT, T)}(fillvalue, data, args...)
end

function PaddedView{FT}(fillvalue,data::AbstractArray{T,N},
			sz::Tuple{Integer,Vararg{Integer}},dims=1) where {FT,T,N}
	if(dims==0)
	  val=copy(sz[1])
	 new_tuple=(size(data)[2],val)
	elseif(dims==1)
	  val=copy(sz[2])
	 new_tuple=(val,size(data)[1])
	end
    inds = map(OneTo, new_tuple)
	dims=1
    PaddedView{FT,N,typeof(inds),typeof(data)}(convert(FT, fillvalue), data, inds)
end

filltype(::Type, ::Type{T}) where T = T
filltype(::Type{FT}, ::Type{T}) where {FT<:Union{Nothing, Missing}, T} = Union{FT, T}
filltype(::Type{FT}, ::Type{T}) where {FT, T<:Union{Nothing, Missing}} = Union{FT, T}
# ambiguity patch
filltype(::Type{FT}, ::Type{T}) where {FT<:Union{Nothing, Missing}, T<:Union{Nothing, Missing}} = Union{FT, T}

Base.axes(A::PaddedView) = A.indices
@inline Base.axes(A::PaddedView, d::Integer) = d <= ndims(A) ? A.indices[d] : default_axes(A.indices)
default_axes(::NTuple{N,I}) where {N,I<:AbstractUnitRange} = convert(I, OneTo(1))
default_axes(::Tuple{}) = OneTo(1)
default_axes(::Any) = OneTo(1)

Base.size(A::PaddedView) = map(length, axes(A))

Base.parent(A::PaddedView) = A.data

Base.@propagate_inbounds function Base.getindex(A::PaddedView{T,N}, i::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, i...)
    if Base.checkbounds(Bool, A.data, i...)
        return convert(T, A.data[i...])
    end
    return A.fillvalue
end

function Base.showarg(io::IO, A::PaddedView, toplevel)
    print(io, "PaddedView(", A.fillvalue, ", ")
    Base.showarg(io, parent(A), false)
    print(io, ", (", join(A.indices, ", "))
    print(io, ndims(A) == 1 ? ",))" : "))")
    toplevel && print(io, " with eltype ", eltype(A))
end

end

# ╔═╡ e02da704-8416-11eb-2529-99ef5a991fbf
begin
img_source = testimage("woman_darkhair")
pyramid = gaussian_pyramid(img_source, 5, 2, 1)
end

# ╔═╡ 6e54a776-8417-11eb-0635-f7686eaabd70
PaddedView(-1, pyramid[2], (512,512),0) # 0 if we want update in rows direction, and 1 if we want to update in columns direction

# ╔═╡ f7565dc2-8416-11eb-3a3e-0130d446dd84

# ╔═╡ Cell order:
# ╠═c52e11fe-83fb-11eb-29b7-77aeaf7ebc83
# ╠═8228eedc-8416-11eb-0daa-81aacd64af55
# ╠═6e54a776-8417-11eb-0635-f7686eaabd70
# ╠═e02da704-8416-11eb-2529-99ef5a991fbf
# ╠═f7565dc2-8416-11eb-3a3e-0130d446dd84

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
using Plots

Rotation(R, θ) = [cos(θ) -sin(θ);sin(θ) cos(θ)] * R

function DragonFractal(;order, θ, points)
    θ₁ = θ[1]; θ₂ = θ[2]
    for num in 1:order
        i = 1
        while i < length(points)
            point = (points[i+1] - points[i]) * sin(abs(θ₁))
            point = Rotation(point, θ₁)  + points[i]
            insert!(points,i+1,point)
            i += 2
            point = (points[i+1] - points[i]) * sin(abs(θ₂))
            point = Rotation(point, θ₂)  + points[i]
            insert!(points,i+1,point)
            i += 2
        end
    end
    return points
end

for i in [1,5,10,15]
    Parameters = Dict(:order => i,
                    :θ => [pi/4, -pi/4],
                    :points => [[0.0; 0.0],[1.0; 0.0],[1.0; 1.0]])

    points = DragonFractal(;Parameters...)
    half = Int(floor(length(points)/2))

    plot(hcat(points...)[1,1:half], hcat(points...)[2,1:half], legend = false, color =:blue, border=:none)
    plot!(hcat(points...)[1,half:end], hcat(points...)[2,half:end], legend = false, color =:red, border=:none)
    savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet1\\Figs\\Q2\\Heighway-dragon-O$i.png")

    Parameters = Dict(:order => i,
                    :θ => [-pi/4, -pi/4],
                    :points => [[0.0; 0.0],[1.0; 0.0],[1.0; 1.0]])
    points = DragonFractal(;Parameters...)
    plot(hcat(points...)[1,:], hcat(points...)[2,:], legend = false, border=:none, color =:black, linewidth=0.4)
    savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet1\\Figs\\Q2\\Lévy-C-curve-O$i.png")
end


text:
const LOG_FORMAT = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sss")

function log_stream(stream::Requests.ResponseStream, log_path::String; verbose::Bool = true)
    while true
        tick = next(stream)
        open(log_path,"a") do f
            write(f, Dates.format(tick.recd, LOG_FORMAT), " ")
            write(f, Dates.format(tick.sent, LOG_FORMAT), " ")
            if verbose
                write(f, string(tick.bid), " ")
                write(f, string(tick.ask), " ")
            end
            write(f, tick.pair, "\n")
        end
    end
end

function log_tick(tick::Tick, log_path::String)
    open(log_path,"a") do f
        write(f, JSON.json(tick), "\n")
    end
end

function log_json(stream::Requests.ResponseStream, log_path::String)
    while true
        log_tick(next(stream), log_path)
    end
end

function log_db(stream::Requests.ResponseStream, db::SQLite.DB, table::String)
    while true
        SQLite.query(db, "INSERT INTO $table VALUES $(string(next(stream))[11:end])");
    end
end

max_stars_repo_name: cadwallader/Oanda

max_issues_repo_name: cadwallader/Oanda

max_forks_repo_name: cadwallader/Oanda


text:
module TestEnsembles

# using Revise
using Test
using Random
using MLJ
using MLJBase
import MLJModels
using CategoricalArrays
using Distributions

@load KNNRegressor

## HELPER FUNCTIONS

@test MLJ._reducer([1, 2], [3, ]) == [1, 2, 3]
@test MLJ._reducer(([1, 2], [:x, :y]), ([3, 4, 5], [:z, :w, :a])) ==
    ([1, 2, 3, 4, 5], [:x, :y, :z, :w, :a])

pair_vcat(p, q) = (vcat(p[1], q[1]), vcat(p[2], q[2]))

## WRAPPED ENSEMBLES OF FITRESULTS

# target is :deterministic :multiclass false:
atom = MLJModels.DeterministicConstantClassifier()
L = ['a', 'b', 'j']
L2 = categorical(L)
ensemble = [L2[1], L2[3], L2[3], L2[2]]
n=length(ensemble)
atomic_weights = fill(1/n, n) # ignored by predict below
wens = MLJ.WrappedEnsemble(atom, ensemble)
X = MLJ.table(rand(3,5))
@test predict(wens, atomic_weights, X) == categorical(vcat(['j','j','j'],L))[1:3]

# target is :deterministic :continuous false:
atom = MLJModels.DeterministicConstantRegressor()
ensemble = Float64[4, 7, 4, 4]
atomic_weights = [0.1, 0.5, 0.2, 0.2]
wens = MLJ.WrappedEnsemble(atom, ensemble)
@test predict(wens, atomic_weights, X) ≈ [5.5, 5.5, 5.5]

# target is :probabilistic :multiclass false:
atom = ConstantClassifier()
L = categorical(['a', 'b', 'j'])
d1 = UnivariateFinite(L, [0.1, 0.2, 0.7])
d2 = UnivariateFinite(L, [0.2, 0.3, 0.5])
ensemble = [d2,  d1, d2, d2]
atomic_weights = [0.1, 0.5, 0.2, 0.2]
wens = MLJ.WrappedEnsemble(atom, ensemble)
X = MLJ.table(rand(2,5))
d = predict(wens, atomic_weights, X)[1]
@test pdf(d, 'a') ≈ 0.15
@test pdf(d, 'b') ≈ 0.25
@test pdf(d, 'j') ≈ 0.6

# target is :probabilistic :continuous false:
atom = ConstantRegressor()
d1 = Distributions.Normal(1, 2)
d2 = Distributions.Normal(3, 4)
ensemble = [d2,  d1, d2, d2]
atomic_weights = [0.1, 0.5, 0.2, 0.2]
wens = MLJ.WrappedEnsemble(atom, ensemble)
X = MLJ.table(rand(2,5))
d = predict(wens, atomic_weights, X)[1]

## ENSEMBLE MODEL

# target is :deterministic :multiclass false:
atom=MLJModels.DeterministicConstantClassifier()
X = MLJ.table(ones(5,3))
y = categorical(collect("asdfa"))
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
atomic_weights = rand(10)
atomic_weights = atomic_weights/sum(atomic_weights)
ensemble_model.atomic_weights = atomic_weights
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
p = predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)
@test MLJBase.target_scitype(ensemble_model) == MLJBase.target_scitype(atom)

# target is :deterministic :continuous false:
atom = MLJModels.DeterministicConstantRegressor()
X = MLJ.table(ones(5,3))
y = Float64[1.0, 2.0, 1.0, 1.0, 1.0]
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
@test reduce(* , [x ≈ 1.0 || x ≈ 1.25 for x in fitresult.ensemble])
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
ensemble_model.bagging_fraction = 1.0
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
@test unique(fitresult.ensemble) ≈ [1.2]
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
atomic_weights = rand(10)
atomic_weights = atomic_weights/sum(atomic_weights)
ensemble_model.atomic_weights = atomic_weights
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)

# target is :deterministic :continuous false:
atom = MLJModels.DeterministicConstantRegressor()
Random.seed!(1234)
X = MLJ.table(randn(10,3))
y = randn(10)
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom, rng=1)
ensemble_model.out_of_bag_measure = [MLJ.rms,MLJ.rmsp]
ensemble_model.n = 2
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
# TODO: the following test fails in distributed version (because of multiple rng's ?)
@test abs(report.oob_measurements[1] - 1.0834) < 0.001
ensemble_model = MLJ.DeterministicEnsembleModel(atom=atom,rng=Random.MersenneTwister(1))
ensemble_model.out_of_bag_measure = MLJ.rms
ensemble_model.n = 2
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)

# target is :probabilistic :multiclass false:
atom = ConstantClassifier()
X = MLJ.table(ones(5,3))
y = categorical(collect("asdfa"))
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.ProbabilisticEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
fitresult.ensemble
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
ensemble_model.bagging_fraction = 1.0
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
fitresult.ensemble
d = predict(ensemble_model, fitresult, MLJ.selectrows(X, test))[1]
@test pdf(d, 'a') ≈ 2/5
@test pdf(d, 's') ≈ 1/5
@test pdf(d, 'd') ≈ 1/5
@test pdf(d, 'f') ≈ 1/5
@test mode(d) == 'a'
atomic_weights = rand(10)
atomic_weights = atomic_weights/sum(atomic_weights)
ensemble_model.atomic_weights = atomic_weights
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)
# test sample weights
w = [1,100,1,1,1]
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y, w)
p2 = predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
@test mode(p2[1] ) == 's'

# target is :probabilistic :continuous false:
atom = ConstantRegressor()
X = MLJ.table(ones(5,3))
y = Float64[1.0, 2.0, 2.0, 1.0, 1.0]
train, test = partition(1:length(y), 0.8);
ensemble_model = MLJ.ProbabilisticEnsembleModel(atom=atom)
ensemble_model.n = 10
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
d1 = fit(Distributions.Normal, [1,1,2,2])
d2 = fit(Distributions.Normal, [1,1,1,2])
# @test reduce(* , [d.μ ≈ d1.μ || d.μ ≈ d2.μ for d in fitresult.ensemble])
# @test reduce(* , [d.σ ≈ d1.σ || d.σ ≈ d2.σ for d in fitresult.ensemble])
d=predict(ensemble_model, fitresult, MLJ.selectrows(X, test))[1]
for dc in d.components
    @test pdf(dc, 1.52) ≈ pdf(d1, 1.52) || pdf(dc, 1.52) ≈ pdf(d2, 1.52)
end
ensemble_model.bagging_fraction = 1.0
fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
d = predict(ensemble_model, fitresult, MLJ.selectrows(X, test))[1]
d3 = fit(Distributions.Normal, y)
@test pdf(d, 1.52) ≈ pdf(d3, 1.52)
atomic_weights = rand(10)
atomic_weights = atomic_weights/sum(atomic_weights)
ensemble_model.atomic_weights = atomic_weights
predict(ensemble_model, fitresult, MLJ.selectrows(X, test))
MLJBase.info_dict(ensemble_model)
# @test MLJBase.output_is(ensemble_model) == MLJBase.output_is(atom)

# test generic constructor:
@test EnsembleModel(atom=ConstantRegressor()) isa Probabilistic
@test EnsembleModel(atom=MLJModels.DeterministicConstantRegressor()) isa Deterministic

@testset "further test of sample weights" begin
    N = 20
    X = (x = rand(3N), );
    y = categorical(rand("abbbc", 3N));
    atom = @load KNNClassifier
    ensemble_model = MLJ.ProbabilisticEnsembleModel(atom=atom,
                                                    bagging_fraction=1,
                                                    n = 5)
    fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
    @test predict_mode(ensemble_model, fitresult, (x = [0, ],))[1] == 'b'
    w = map(y) do η
        η == 'a' ? 100 : 1
    end
    fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y, w)
    @test predict_mode(ensemble_model, fitresult, (x = [0, ],))[1] == 'a'

    ensemble_model.rng = 1234 # always start with same rng
    ensemble_model.bagging_fraction=0.6
    ensemble_model.out_of_bag_measure = [BrierScore(), cross_entropy]
    fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
    m1 = report.oob_measurements[1]
    fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y)
    m2 = report.oob_measurements[1]
    @test m1 == m2
    # supplying sample weights should change the oob meausurements for
    # measures that support weights:
    fitresult, cache, report = MLJ.fit(ensemble_model, 1, X, y, w)
    m3 = report.oob_measurements[1]
    @test !(m1 ≈ m3)
end

    ## MACHINE TEST (INCLUDES TEST OF UPDATE)

N =100
X = (x1=rand(N), x2=rand(N), x3=rand(N))
y = 2X.x1  - X.x2 + 0.05*rand(N)

atom = KNNRegressor(K=7)
ensemble_model = EnsembleModel(atom=atom)
ensemble = machine(ensemble_model, X, y)
train, test = partition(eachindex(y), 0.7)
fit!(ensemble, rows=train)
@test length(ensemble.fitresult.ensemble) == ensemble_model.n
ensemble_model.n = 15
@test_logs((:info, r"Training"),
          fit!(ensemble))
@test length(ensemble.fitresult.ensemble) == 15
ensemble_model.n = 20
@test_logs((:info, r"Updating"),
           (:info, r"Building"),
           fit!(ensemble))
@test length(ensemble.fitresult.ensemble) == 20
ensemble_model.n = 5
@test_logs((:info, r"Updating"),
           (:info, r"Truncating"),
           fit!(ensemble))
@test length(ensemble.fitresult.ensemble) == 5

@test !isnan(predict(ensemble, MLJ.selectrows(X, test))[1])

end
true


text:
# A linked list node used to implement a stack trace within Salsa. This trace is used for
# debug tracing (to understand incrementality) and is included in Salsa wrapped Exceptions
# for improved error information.
#
# We use a linked list for the stack since it is immutable and so will be thread-safe, and
# because the linked list can be cheaply branched (becoming a "tree") for derived functions
# called concurrently with the same parent.
struct SalsaStackFrame
    dp::DependencyKey
    next::Union{Nothing,SalsaStackFrame}
end

function Base.collect(stack::SalsaStackFrame)
    out = DependencyKey[stack.dp]
    while stack.next !== nothing
        stack = stack.next
        push!(out, stack.dp)
    end
    return out
end

Base.length(stack::SalsaStackFrame) = stack.next === nothing ? 1 : length(stack.next) + 1

# For cycle detection in Debug-mode:
# NOTE: The quadratic performance hazard here; this is to be disabled in production.
function stack_has_key(stack::Union{SalsaStackFrame,Nothing}, key)
    # NOTE: Important that this is implemented as a while-loop, not a recursive function
    # because it may be called already deep in the stack, and this could effectively double
    # the current stack height. (It seems that the tail-recursion optimization is missed
    # with juila -O0, so it's better to do this manually.)
    while stack !== nothing
        if stack.dp == key
            return true
        end
        stack = stack.next
    end
    return false
end


text:
module NeuralPDELogging

using NeuralPDE
using TensorBoardLogger

"""This function overrides the empty function in NeuralPDE in order to use TensorBoardLogger in that package
This is light type piracy but it should be alright since this is a subpackage of NeuralPDE"""
function NeuralPDE.logvector(logger::TBLogger, vector::AbstractVector{R}, name::AbstractString, step::Integer) where R <: Real
    for j in 1:length(vector)
        log_value(logger, "$(name)/$(j)", vector[j], step=step)
    end
    nothing
end

"""This function overrides the empty function in NeuralPDE in order to use TensorBoardLogger in that package.
This is light type piracy but it should be alright since this is a subpackage of NeuralPDE"""
function NeuralPDE.logscalar(logger::TBLogger, scalar::R, name::AbstractString, step::Integer) where R <: Real
    log_value(logger, "$(name)", scalar, step=step)
    nothing
end

end


text:
struct SDFFile{P,B,C,PC}
    name::String
    header::Header
    blocks::B
    param::Ref{P}
    cache::Ref{C}
    particle_cache::Ref{PC}
end

include("read.jl")
include("utils.jl")
include("angular_momentum.jl")

Base.keys(sdf::SDFFile) = keys(sdf.blocks)
Base.getindex(sdf::SDFFile, idx::Vararg{AbstractString, N}) where N = sdf[Symbol.(idx)...]
Base.haskey(sim::SDFFile, key) = haskey(sim.param[], key)
Base.haskey(sim::SDFFile, block, key) = haskey(sim.param[][block], key)

# show
function Base.show(io::IO, ::MIME"text/plain", sdf::SDFFile)
    t = si_round(get_time(sdf))
    n = length(keys(sdf))
    print(io, "SDFFile with $n entries at t = " * t)
end


text:
function labyrinth(r)
    sides = []

    # функция, которая выводит робота в левый верхний угол
    function go_to_angle()
        while (!isborder(r, Nord) || !isborder(r, West))
            while !isborder(r, Nord)
                move!(r, Nord)
                push!(sides, Nord)
            end
            while !isborder(r, West)
                move!(r, West)
                push!(sides, West)
            end
        end
    end

    function put_marker_to_angle(side)
        while !isborder(r, side)
            move!(r, side)
        end
        putmarker!(r)
    end

    function solve()
        go_to_angle()
        for i in [Sud, Ost, Nord, West]
            put_marker_to_angle(i)
        end

        while length(sides) > 0
            move!(r, inverse_side(pop!(sides)))
        end
    end

    solve()
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
abstract type AbstractRegulationFormulation <: AbstractDeviceFormulation end
struct ReserveLimitedRegulation <: AbstractRegulationFormulation end
struct DeviceLimitedRegulation <: AbstractRegulationFormulation end

"""
This function add the upwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{DeltaActivePowerUpVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.RegulationDevice}
    return AddVariableSpec(;
        variable_name = make_variable_name(DeltaActivePowerUpVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

"""
This function add the downwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{DeltaActivePowerDownVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.RegulationDevice}
    return AddVariableSpec(;
        variable_name = make_variable_name(DeltaActivePowerDownVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

"""
This function add the upwards scheduled regulation variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{AdditionalDeltaActivePowerUpVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.RegulationDevice}
    return AddVariableSpec(;
        variable_name = make_variable_name(AdditionalDeltaActivePowerUpVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

"""
This function add the variables for power generation output to the model
"""
function AddVariableSpec(
    ::Type{AdditionalDeltaActivePowerDownVariable},
    ::Type{U},
    psi_container::PSIContainer,
) where {U <: PSY.RegulationDevice}
    return AddVariableSpec(;
        variable_name = make_variable_name(AdditionalDeltaActivePowerDownVariable, U),
        binary = false,
        lb_value_func = x -> 0.0,
    )
end

function add_constraints!(
    psi_container::PSIContainer,
    ::Type{RangeConstraint},
    ::Type{DeltaActivePowerUpVariable},
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, DeviceLimitedRegulation},
    ::Type{AreaBalancePowerModel},
    ::Nothing,
) where {T <: PSY.StaticInjection}
    parameters = model_has_parameters(psi_container)
    var_name_up = make_variable_name(DeltaActivePowerUpVariable, T)
    var_up = get_variable(psi_container, var_name_up)

    names = [PSY.get_name(g) for g in devices]
    time_steps = model_time_steps(psi_container)

    up = Symbol("regulation_limits_up_$(T)")
    container_up = add_cons_container!(psi_container, up, names, time_steps)

    constraint_infos = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, "get_max_active_power")
        constraint_info = DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_max_active_power(x),
            ts_vector,
            x -> PSY.get_active_power_limits(x),
        )
        constraint_infos[ix] = constraint_info
    end

    if parameters
        base_points_param =
            get_parameter_container(psi_container, make_variable_name(ACTIVE_POWER, T))
        multiplier = get_multiplier_array(base_points_param)
        base_points = get_parameter_array(base_points_param)
    end

    for d in constraint_infos
        name = get_component_name(d)
        limits = get_limits(d)
        for t in time_steps
            rating = parameters ? multiplier[name, t] : d.multiplier
            base_point = parameters ? base_points[name, t] : get_timeseries(d)[t]
            container_up[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                var_up[name, t] <= limits.max - base_point * rating
            )
        end
    end
    return
end

function add_constraints!(
    psi_container::PSIContainer,
    ::Type{RangeConstraint},
    ::Type{DeltaActivePowerDownVariable},
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, DeviceLimitedRegulation},
    ::Type{AreaBalancePowerModel},
    ::Nothing,
) where {T <: PSY.StaticInjection}
    parameters = model_has_parameters(psi_container)
    var_name_dn = make_variable_name(DeltaActivePowerDownVariable, T)
    var_dn = get_variable(psi_container, var_name_dn)

    names = [PSY.get_name(g) for g in devices]
    time_steps = model_time_steps(psi_container)

    dn = Symbol("regulation_limits_dn_$(T)")
    container_dn = add_cons_container!(psi_container, dn, names, time_steps)

    constraint_infos = Vector{DeviceTimeSeriesConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        ts_vector = get_time_series(psi_container, d, "get_max_active_power")
        constraint_info = DeviceTimeSeriesConstraintInfo(
            d,
            x -> PSY.get_max_active_power(x),
            ts_vector,
            x -> PSY.get_active_power_limits(x),
        )
        constraint_infos[ix] = constraint_info
    end

    if parameters
        base_points_param =
            get_parameter_container(psi_container, make_variable_name(ACTIVE_POWER, T))
        multiplier = get_multiplier_array(base_points_param)
        base_points = get_parameter_array(base_points_param)
    end

    for d in constraint_infos
        name = get_component_name(d)
        limits = get_limits(d)
        for t in time_steps
            rating = parameters ? multiplier[name, t] : d.multiplier
            base_point = parameters ? base_points[name, t] : get_timeseries(d)[t]
            container_dn[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                var_dn[name, t] <= base_point * rating - limits.min
            )
        end
    end
    return
end

function add_constraints!(
    psi_container::PSIContainer,
    ::Type{RangeConstraint},
    ::Type{DeltaActivePowerUpVariable},
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, ReserveLimitedRegulation},
    ::Type{AreaBalancePowerModel},
    ::Nothing,
) where {T <: PSY.StaticInjection}
    var_name_up = make_variable_name(DeltaActivePowerUpVariable, T)
    var_up = get_variable(psi_container, var_name_up)

    names = [PSY.get_name(g) for g in devices]
    time_steps = model_time_steps(psi_container)

    up = Symbol("regulation_limits_up_$(T)")
    container_up = add_cons_container!(psi_container, up, names, time_steps)

    for d in devices
        name = PSY.get_name(d)
        limit_up = PSY.get_reserve_limit_up(d)
        for t in time_steps
            container_up[name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, var_up[name, t] <= limit_up)
        end
    end
    return
end

function add_constraints!(
    psi_container::PSIContainer,
    ::Type{RangeConstraint},
    ::Type{DeltaActivePowerDownVariable},
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, ReserveLimitedRegulation},
    ::Type{AreaBalancePowerModel},
    ::Nothing,
) where {T <: PSY.StaticInjection}
    var_name_dn = make_variable_name(DeltaActivePowerDownVariable, T)
    var_dn = get_variable(psi_container, var_name_dn)

    names = [PSY.get_name(g) for g in devices]
    time_steps = model_time_steps(psi_container)

    dn = Symbol("regulation_limits_dn_$(T)")
    container_dn = add_cons_container!(psi_container, dn, names, time_steps)

    for d in devices
        name = PSY.get_name(d)
        limit_up = PSY.get_reserve_limit_dn(d)
        for t in time_steps
            container_dn[name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, var_dn[name, t] <= limit_up)
        end
    end
    return
end

function ramp_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, DeviceLimitedRegulation},
    ::Type{AreaBalancePowerModel},
    ::Nothing,
) where {T <: PSY.ThermalStandard}
    R_up = get_variable(psi_container, DeltaActivePowerUpVariable, T)
    R_dn = get_variable(psi_container, DeltaActivePowerDownVariable, T)

    resolution = Dates.value(Dates.Second(model_resolution(psi_container)))
    names = [PSY.get_name(g) for g in devices]
    time_steps = model_time_steps(psi_container)

    container_up = add_cons_container!(psi_container, :ramp_limits_up, names, time_steps)
    container_dn = add_cons_container!(psi_container, :ramp_limits_dn, names, time_steps)

    for d in devices
        ramp_limits = PSY.get_ramp_limits(d)
        isnothing(ramp_limits) && continue
        scaling_factor = resolution * SECONDS_IN_MINUTE
        name = PSY.get_name(d)
        for t in time_steps
            container_up[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                R_up[name, t] <= ramp_limits.up * scaling_factor
            )
            container_dn[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                R_dn[name, t] <= ramp_limits.down * scaling_factor
            )
        end
    end
    return
end

function participation_assignment!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, <:AbstractRegulationFormulation},
    ::Type{AreaBalancePowerModel},
    ::Nothing,
) where {T <: PSY.StaticInjection}
    time_steps = model_time_steps(psi_container)
    R_up = get_variable(psi_container, DeltaActivePowerUpVariable, T)
    R_dn = get_variable(psi_container, DeltaActivePowerDownVariable, T)
    R_up_emergency = get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, T)
    R_dn_emergency = get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, T)
    area_reserve_up = get_variable(psi_container, DeltaActivePowerUpVariable, PSY.Area)
    area_reserve_dn = get_variable(psi_container, DeltaActivePowerDownVariable, PSY.Area)

    component_names = [PSY.get_name(d) for d in devices]
    participation_assignment_up = JuMPConstraintArray(undef, component_names, time_steps)
    participation_assignment_dn = JuMPConstraintArray(undef, component_names, time_steps)
    assign_constraint!(
        psi_container,
        "participation_assignment_up",
        participation_assignment_up,
    )
    assign_constraint!(
        psi_container,
        "participation_assignment_dn",
        participation_assignment_dn,
    )

    expr_up = get_expression(psi_container, :emergency_up)
    expr_dn = get_expression(psi_container, :emergency_dn)
    for d in devices
        name = PSY.get_name(d)
        services = PSY.get_services(d)
        if length(services) > 1
            device_agc = (a for a in PSY.get_services(d) if isa(a, PSY.AGC))
            area_name = PSY.get_name.(PSY.get_area.(device_agc))[1]
        else
            device_agc = first(services)
            area_name = PSY.get_name(PSY.get_area(device_agc))
        end
        p_factor = PSY.get_participation_factor(d)
        for t in time_steps
            participation_assignment_up[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                R_up[name, t] ==
                (p_factor.up * area_reserve_up[area_name, t]) + R_up_emergency[name, t]
            )
            participation_assignment_dn[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                R_dn[name, t] ==
                (p_factor.dn * area_reserve_dn[area_name, t]) + R_dn_emergency[name, t]
            )
            JuMP.add_to_expression!(expr_up[area_name, t], -1 * R_up_emergency[name, t])
            JuMP.add_to_expression!(expr_dn[area_name, t], -1 * R_dn_emergency[name, t])
        end

    end

    return
end

function regulation_cost!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{PSY.RegulationDevice{T}},
    ::DeviceModel{PSY.RegulationDevice{T}, <:AbstractRegulationFormulation},
) where {T <: PSY.StaticInjection}
    time_steps = model_time_steps(psi_container)
    R_up = get_variable(psi_container, DeltaActivePowerUpVariable, T)
    R_dn = get_variable(psi_container, DeltaActivePowerDownVariable, T)
    R_up_emergency = get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, T)
    R_dn_emergency = get_variable(psi_container, AdditionalDeltaActivePowerUpVariable, T)

    for d in devices
        cost = PSY.get_cost(d)
        p_factor = PSY.get_participation_factor(d)
        up_cost =
            isapprox(p_factor.up, 0.0; atol = 1e-2) ? SERVICES_SLACK_COST : 1 / p_factor.up
        dn_cost =
            isapprox(p_factor.dn, 0.0; atol = 1e-2) ? SERVICES_SLACK_COST : 1 / p_factor.dn
        for t in time_steps
            JuMP.add_to_expression!(
                psi_container.cost_function,
                R_up_emergency[PSY.get_name(d), t],
                up_cost,
            )
            JuMP.add_to_expression!(
                psi_container.cost_function,
                R_dn_emergency[PSY.get_name(d), t],
                dn_cost,
            )
        end
    end
    return
end

function NodalExpressionSpec(
    ::Type{<:PSY.RegulationDevice{T}},
    ::Type{AreaBalancePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.StaticInjection}
    return NodalExpressionSpec(
        "get_max_active_power",
        make_variable_name(ActivePowerVariable, T),
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        1.0,
        JuMP.VariableRef,
    )
end


text:
quoted(x) = typof(x) === Quote || typof(x) === Quotenode
unquoted(x) = typof(x) === UnaryOpCall && typof(x.args[1]) === OPERATOR && kindof(x.args[1]) == CSTParser.Tokens.EX_OR

function get_ids(x, q = false, ids = [])
    if quoted(x)
        q = true
    end
    if q && unquoted(x)
        q = false
    end
    if isidentifier(x)
        !q && push!(ids, x)
    elseif x.args !== nothing
        for i in 1:length(x.args)
            get_ids(x.args[i], q, ids)
        end
    end
    ids
end

function collect_bindings_refs(x::EXPR, bindings = [], refs = [])
    if bindingof(x) !== nothing
        push!(bindings, x)
    end
    if StaticLint.hasref(x)
        push!(refs, x)
    end
    if x.args !== nothing
        for a in x.args
            collect_bindings_refs(a, bindings, refs)
        end
    end
    return bindings, refs
end

function remove_ref(x::EXPR)
    if hasref(x) && refof(x) isa Binding
        for ia in enumerate(refof(x).refs)
            if ia[2] == x
                deleteat!(refof(x).refs, ia[1])
                setref!(x, nothing)
                return
            end
        end
        error()
    end
end

function clear_binding(x::EXPR)
    if bindingof(x) isa Binding
        for r in bindingof(x).refs
            if r isa EXPR
                setref!(r, nothing)
            elseif r isa Binding
                clear_binding(r)
            end
        end
        x.binding = nothing
    end
end
function clear_scope(x::EXPR)
    if x.scope isa Scope
        setparent!(x.scope, nothing)
        empty!(x.scope.names)
        if typof(x) === FileH && x.scope.modules isa Dict && haskey(x.scope.modules, "Base") && haskey(x.scope.modules, "Core")
            m1 = x.scope.modules["Base"]
            m2 = x.scope.modules["Core"]
            empty!(x.scope.modules)
            x.scope.modules["Base"] = m1
            x.scope.modules["Core"] = m2
        else
            x.scope.modules = nothing
        end
    end
end

function clear_ref(x::EXPR)
    if refof(x) isa Binding
        for i in 1:length(refof(x).refs)
            if refof(x).refs[i] == x
                deleteat!(refof(x).refs, i)
                break
            end
        end
        setref!(x, nothing)
    elseif refof(x) !== nothing
        setref!(x, nothing)
    end
end
function clear_meta(x::EXPR)
    clear_binding(x)
    clear_ref(x)
    clear_scope(x)
    if x.args !== nothing
        for a in x.args
            clear_meta(a)
        end
    end
end

function get_root_method(b, server)
    return b
end

function get_root_method(b::Binding, server, b1 = nothing, bs = Binding[])
    if b.overwrites === nothing || b == b.overwrites || !(b.overwrites isa Binding)
        return b
    elseif b.overwrites.t == getsymbolserver(server)["Core"].vals["Function"]
        return get_root_method(b.overwrites, server, b, bs)
    elseif b.overwrites.t == getsymbolserver(server)["Core"].vals["DataType"]
        return b.overwrites
    else
        return b
    end
end

function retrieve_delayed_scope(x)
    if (CSTParser.defines_function(x) || CSTParser.defines_macro(x)) && scopeof(x) !== nothing
        if parentof(scopeof(x)) !== nothing
            return parentof(scopeof(x))
        else
            return scopeof(x)
        end
    elseif typof(x) === Export
        return retrieve_scope(x)
    elseif parentof(x) !== nothing
        return retrieve_delayed_scope(parentof(x))
    end
    return nothing
end

function retrieve_scope(x)
    if scopeof(x) != nothing
        return scopeof(x)
    elseif parentof(x) isa EXPR
        return retrieve_scope(parentof(x))
    end
    return
end