text:
using CSV
using DelimitedFiles
using DataFrames
using Dates
using Plots
using TimerOutputs
using BenchmarkTools
using Revise
using PDFmerger

include("read_in_cases.jl")
include("generative_model.jl")
include("processing_sim.jl")
include("plot_forecasts.jl")

function simulate_all_states(file_date, states_to_run, nsims)
    """
    This runs the branching process for each of the states in states_to_run and then
    merges and saves the files.
    """

    # set seed for consistent plots (NOTE: this is not useful when multithreading
    # enabled as we use separate seeds but the simulation pool should handle that)
    rng = Random.Xoshiro(2022)

    simulation_start_dates = Dict{String, String}(
        "NSW" => "2021-06-23",
        "QLD" => "2021-11-01",
        "SA" => "2021-11-01",
        "TAS" => "2021-11-01",
        "WA" => "2021-11-01",
        "ACT" => "2021-08-01",
        "NT" => "2021-11-01",
        "VIC" => "2021-08-01",
    )

    # date we want to apply increase in cases due to Omicron
    omicron_dominant_date = Dates.Date("2021-12-10")

    # get the latest onset date
    latest_start_date = Dates.Date(maximum(v for v in values(simulation_start_dates)))

    pop_sizes = Dict{String, Int}(
        "NSW" => 8189266,
        "QLD" => 5221170,
        "SA" => 1773243,
        "TAS" => 541479,
        "VIC" => 6649159,
        "WA" => 2681633,
        "ACT" => 432266,
        "NT" => 246338,
    )

    initial_conditions = Dict{String, NamedTuple{(:S, :A, :I), Tuple{Int64, Int64, Int64}}}(
        "NSW" => (S = 5, A = 8, I = 0),
        "QLD" => (S = 0, A = 0, I = 0),
        "SA" => (S = 0, A = 0, I = 0),
        "TAS" => (S = 0, A = 0, I = 0),
        "VIC" => (S = 20, A = 20, I = 0),
        "WA" => (S = 0, A = 0, I = 0),
        "ACT" => (S = 0, A = 0, I = 0),
        "NT" => (S = 0, A = 0, I = 0),
    )

    (local_case_dict, import_case_dict) = read_in_cases(file_date, rng)
    dates = local_case_dict["date"]
    last_date_in_data = dates[end]
    forecast_end_date = last_date_in_data + Dates.Day(35)

    # create vector for dates
    onset_dates = latest_start_date:Dates.Day(1):forecast_end_date

    for state in states_to_run

        forecast_start_date = Dates.Date(simulation_start_dates[state])
        local_cases = local_case_dict[state]
        import_cases = import_case_dict[state]

        # named tuple for initial condi
        D0 = initial_conditions[state]
        N = pop_sizes[state]

        # get the observed cases
        cases_pre_forecast = sum(local_case_dict[state][dates .< forecast_start_date])
        local_cases = local_case_dict[state][dates .>= forecast_start_date]
        import_cases = import_case_dict[state]

        (D, U, TP_local) = simulate_branching_process(
            D0,
            N,
            nsims,
            local_cases,
            import_cases,
            cases_pre_forecast,
            forecast_start_date,
            file_date,
            omicron_dominant_date,
            state,
        )

        save_simulations(D,TP_local,state,file_date,onset_dates,rng)

    end

    # merge all the simulation and TP files for states into single CSV
    merge_simulation_files(file_date)
    merge_TP_files(file_date)
    plot_all_forecast_intervals(file_date, states_to_run, local_case_dict)

    return nothing

end

function plot_all_forecast_intervals(file_date, states, local_case_dict)
    """
    Simple wrapper function to plot the forecasts with various zoom and confidence levels.
    """

    plot_all_forecasts(file_date, states_to_run, local_case_dict, confidence_level="both")
    plot_all_forecasts(file_date, states_to_run, local_case_dict, zoom=true, confidence_level="both")
    plot_all_forecasts(file_date, states_to_run, local_case_dict, confidence_level="50")
    plot_all_forecasts(file_date, states_to_run, local_case_dict, confidence_level="95")

    file_name_tmp = "UoA_forecast_"
    dir_name = "figs/case_forecasts/"*file_date*"/"

    pdf_filenames = [
        dir_name*file_name_tmp*file_date*"_zoomed_both_intervals.pdf",
        dir_name*file_name_tmp*file_date*"_both_intervals.pdf",
    ]
    # merge the pdfs and delete the files
    merge_pdfs(
        pdf_filenames,
        dir_name*file_name_tmp*file_date*".pdf",
        cleanup=true,
    )

    return nothing
end

function simulate_single_state(file_date,state,nsims)
    """
    This runs the branching process for each of the states in states_to_run and then
    merges and saves the files.
    """

    simulation_start_dates = Dict{String, String}(
        "NSW" => "2021-06-23",
        "QLD" => "2021-11-01",
        "SA" => "2021-11-01",
        "TAS" => "2021-11-01",
        "WA" => "2021-11-01",
        "ACT" => "2021-08-01",
        "NT" => "2021-11-01",
        "VIC" => "2021-08-01",
    )

    # get the latest onset date
    latest_onset_date = Dates.Date(maximum(v for v in values(simulation_start_dates)))
    onset_dates = latest_onset_date:Dates.Day(1):latest_onset_date+Dates.Day(size(D,1)-1)

    pop_sizes = Dict{String, Int}(
        "NSW" => 8189266,
        "QLD" => 5221170,
        "SA" => 1773243,
        "TAS" => 541479,
        "VIC" => 6649159,
        "WA" => 2681633,
        "ACT" => 432266,
        "NT" => 246338,
    )

    initial_conditions = Dict{String, NamedTuple{(:S, :A, :I), Tuple{Int64, Int64, Int64}}}(
        "NSW" => (S = 5, A = 5, I = 0),
        "QLD" => (S = 0, A = 0, I = 0),
        "SA" => (S = 0, A = 0, I = 0),
        "TAS" => (S = 0, A = 0, I = 0),
        "VIC" => (S = 20, A = 20, I = 0),
        "WA" => (S = 0, A = 0, I = 0),
        "ACT" => (S = 0, A = 0, I = 0),
        "NT" => (S = 0, A = 0, I = 0),
    )

    (local_case_dict, import_case_dict) = read_in_cases(file_date)
    dates = local_case_dict["date"]

    forecast_start_date = Dates.Date(simulation_start_dates[state])
    local_cases = local_case_dict[state]
    import_cases = import_case_dict[state]

    # named tuple for initial condi
    D0 = initial_conditions[state]
    N = pop_sizes[state]

    # get the observed cases
    cases_pre_forecast = sum(local_case_dict[state][.!plot_ind])
    local_cases = local_case_dict[state][plot_ind]
    import_cases = import_case_dict[state]

    (D, U, TP_local) = simulate_branching_process(
        D0,
        N,
        nsims,
        local_cases,
        import_cases,
        cases_pre_forecast,
        forecast_start_date,
        file_date,
        state,
    )

    save_simulations(D,state,file_date,onset_dates)

    # merge all the simulation and TP files for states into single CSV
    merge_simulation_files(file_date)
    merge_TP_files(file_date)

    # now we can plot everything
    plot_all_forecasts(file_date,state,local_case_dict,confidence_level="both")
    plot_all_forecasts(file_date,state,local_case_dict,confidence_level="50")
    plot_all_forecasts(file_date,state,local_case_dict,confidence_level="95")

    return nothing

end

function visualise_data()

    # choose file date
    file_date = "2022-01-04"

    # states to simulate
    states_to_run = [
        "NSW",
        "QLD",
        "SA",
        "TAS",
        "VIC",
        "WA",
        "ACT",
        "NT",
    ]

    simulation_start_dates = Dict{String, String}(
        "NSW" => "2021-06-23",
        "QLD" => "2021-11-01",
        "SA" => "2021-11-01",
        "TAS" => "2021-11-01",
        "WA" => "2021-11-01",
        "ACT" => "2021-08-01",
        "NT" => "2021-11-01",
        "VIC" => "2021-08-01",
    )

    # get the latest onset date
    latest_onset_date = Dates.Date(maximum(v for v in values(simulation_start_dates)))
    onset_dates = latest_onset_date:Dates.Day(1):latest_onset_date+Dates.Day(size(D,1)-1)

    pop_sizes = Dict{String, Int}(
        "NSW" => 8189266,
        "QLD" => 5221170,
        "SA" => 1773243,
        "TAS" => 541479,
        "VIC" => 6649159,
        "WA" => 2681633,
        "ACT" => 432266,
        "NT" => 246338,
    )

    initial_conditions = Dict{String, NamedTuple{(:S, :A, :I), Tuple{Int64, Int64, Int64}}}(
        "NSW" => (S = 5, A = 5, I = 0),
        "QLD" => (S = 0, A = 0, I = 0),
        "SA" => (S = 0, A = 0, I = 0),
        "TAS" => (S = 0, A = 0, I = 0),
        "VIC" => (S = 20, A = 20, I = 0),
        "WA" => (S = 0, A = 0, I = 0),
        "ACT" => (S = 0, A = 0, I = 0),
        "NT" => (S = 0, A = 0, I = 0),
    )

    (local_case_dict, import_case_dict) = read_in_cases(file_date)
    dates = local_case_dict["date"]

    # indices to actually plot over (simulation horizon)
    plot_ind = dates .>= forecast_start_date

    plot_stuff = false
    if plot_stuff
        # plot the actual cases
        fig = plot(legend=false)
        bar!(
            fig,
            dates[plot_ind],
            local_cases[plot_ind]+import_cases[plot_ind],
            lc="gray",
            fc="gray",
            fillalpha=0.4,
            linealpha=0.4
        )
        bar!(
            fig,
            dates[plot_ind],
            import_cases[plot_ind],
            lc="gray",
            fc="gray",
            fillalpha=0.4,
            linealpha=0.4
        )
    end

end


text:
module ModelTest

using Test
using Elsa
using Elsa: TightbindingModel, OnsiteTerm, HoppingTerm, padtotype

@testset "onsite" begin
    r = SVector(0.0, 0.0)
    os = (1, 1.0, @SMatrix[1 0; 0 1], I, r -> 1, r -> I, r -> SMatrix{3,3}(I))
    ss = (missing, :A, (:A,), (:A, :B))
    cs = (1, 1.0, 1.0f0)
    ts = (Float32, Float64, SMatrix{3,3}, SMatrix{3,2,Float64}, SMatrix{1,4,Float32})
    for o in os, s in ss, c in cs
        ons = c * onsite(o, sublats = s)
        @test ons isa
            TightbindingModel{1,Tuple{OnsiteTerm{typeof(o),S,typeof(c)}}} where S
        term = first(ons.terms)
        for t in ts
            @test padtotype(term(r, r), t) isa t
        end
    end
end

@testset "hopping" begin
    r = SVector(0.0, 0.0)
    hs = (1, 1.0, @SMatrix[1 0; 0 1], I, (r, dr) -> 1, (r, dr) -> I, (r, dr) -> SMatrix{3,3}(I))
    ss = (missing, :A, (:A,), (:A, :B), ((:A,:B), :C), ((:A,:B), (:C,)), ((:A,:B), (:C, :D)))
    cs = (1, 1.0, 1.0f0)
    ts = (Float32, Float64, SMatrix{3,3}, SMatrix{3,2,Float64}, SMatrix{1,4,Float32})
    dns = (missing, (1,), ((1,2), (3,4)))
    rs = (1, 2.0, Inf)
    for h in hs, s in ss, c in cs, d in dns, r in rs
        hop = c * hopping(h, sublats = s, dn = d, range = r)
        @test hop isa
            TightbindingModel{1,Tuple{HoppingTerm{typeof(h),S,D,typeof(float(r)),typeof(c)}}} where {S,D}
        term = first(hop.terms)
        for t in ts
            @test padtotype(term(r, r), t) isa t
        end
    end
end

@testset "term algebra" begin
    r = SVector(0, 0)
    model = onsite(1) + hopping(2I)
    @test (t -> t(r, r)).(model.terms) == (1, 2I)
    model = onsite(1) - hopping(2I)
    @test (t -> t(r, r)).(model.terms) == (1, -2I)
    model = -onsite(@SMatrix[1 0; 1 1]) - 2hopping(2I)
    @test (t -> t(r, r)).(model.terms) == (-@SMatrix[1 0; 1 1], -4I)
    @test model(r, r) == @SMatrix[-5 0; -1 -5]
end

end # module

ext: jl


text:
module NCP_MC

include("../Misc.jl")
include("BackwardKernel.jl")
using .Misc, .BackwardKernel

using Distributions, SpecialFunctions
import Base.Iterators: product

export Particle, smc_weight, filtering_proposal, resample

mutable struct Particle{ℜ <: Real}
    S_R::Array{ℜ}
    S_I::Array{Array{ℜ}}
    function Particle(X::Array, R::Int; a::ℜ=0.0) where {ℜ <: Real}
        return new{ℜ}(fill(a/R,R), [fill(a/(R*Iₙ),R,Iₙ) for Iₙ ∈ size(X)])
    end
    function Particle(S_R::Array{ℜ}, S_I::Array{Array{ℜ}}) where {ℜ <: Real}
        return new{ℜ}(S_R, S_I)
    end
end

function resample(W::Array,Π::Array{Particle{ℜ}},u=rand()) where {ℜ <: Real} # systematic resampling
    P, N = length(W), length(Π[1].S_I)
    j = 0
    cum_Wⱼ = cum_Wᵢ = -u
    for i ∈ 1:P
        rᵢ = ceil(cum_Wᵢ + P*W[i]) - ceil(cum_Wᵢ) # number of replicas for ith particle
        for _ ∈ 2.0:rᵢ
            j+=1
            while ceil(cum_Wⱼ+ P*W[j]) - ceil(cum_Wⱼ) > 0 # find next j to be replaced
                cum_Wⱼ += P*W[j]
                j+=1
            end
            cum_Wⱼ += P*W[j]

            # replace j by i
            Π[j].S_R .= Π[i].S_R
            for n ∈ 1:N
                Π[j].S_I[n] .= Π[i].S_I[n]
            end
        end
        cum_Wᵢ += P*W[i]
    end
end

function smc_weight(X::Array{ℜ,N}, R::Ƶ, P::Ƶ=1; a=1.0, b=a/nansum(X), resampling=true) where {ℜ<:Real, Ƶ<:Integer,N}
    T = Ƶ(round(sum(X)))
    ESS = P

    Π = [Particle(X, R; a=a) for p=1:P]
    log_Z = a*log(b) - (a+T)*log(b + 1) - sum(lgamma, X .+ 1)

    log_w, W = fill(log_Z,P), fill(1.0/P,P)
    log_q, q = zeros(R), zeros(R)

    for i ∈ EventQueue(X)
        for p ∈ 1:P
            log_q .= (1-N) .* log.(Π[p].S_R)
            for  n ∈ 1:N
                log_q .+= log.(Π[p].S_I[n][:,i[n]])
            end

            log_ν = logsumexp(log_q)
            log_q .-= log_ν
            q .= exp.(log_q)
            r = rand(Categorical(q))

            Π[p].S_R[r] += 1.0
            for n ∈ 1:N
                Π[p].S_I[n][r,i[n]] += 1.0
            end

            log_w[p] += log_ν
        end

        log_Z = logmeanexp(log_w)
        W .= exp.(log_w .- logsumexp(log_w))
        ESS = 1.0/sum(W .* W)

        if resampling && ESS < P/2
            resample(W,Π)
            log_w .= log_Z
        end
    end
    p = rand(Categorical(W))
    return log_Z, ESS, Π[p]
end

function filtering_proposal(X::Array{ℜ,N}, R::Ƶ, P::Ƶ=1; a=1.0, b=a/nansum(X)) where {ℜ<:Real, Ƶ<:Integer,N}
    I = size(X)
    α_I = (a/R) ./ I
    α_R = a/R

    S, S_R, S_I = zeros(ℜ,R,I...), zeros(ℜ,R), [zeros(ℜ,R,Iₙ) for Iₙ ∈ size(X)]

    S₊ = sum(X)
    log_Z = a*log(b) - (a+S₊)*log(b + 1) + lgamma(a + S₊) - lgamma(a)
    log_w, log_q, q = fill(log_Z,P), zeros(R), zeros(R)

    for p ∈ 1:P
        S .= 0
        S_R .= 0
        for n ∈ 1:N
            S_I[n] .= 0
        end

        for (t,i) in enumerate(EventQueue(X))
            log_q .=  sum(n -> log.(S_I[n][:,i[n]] .+ α_I[n]),1:N) - (N-1)*log.(S_R .+ α_R) .- log(t+a-1)
            q .= exp.(log_q .- logsumexp(log_q))
            r = rand(Categorical(q))

            S_R[r] += 1
            for n ∈ 1:N
                S_I[n][r,i[n]] += 1
            end
            S[r,i...] += 1

            log_w[p] += log_q[r]
        end
        log_w[p] -= sum(lgamma,S .+ 1)
    end
    return log_w
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
################################################################################
#
#                       arf.jl: wrapper for arf
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015 Tommy Hofmann
#
################################################################################
#
#  TODO:
#  - Sensible way to deal with rounding modes (at least documentation)
#
################################################################################

export arf, ArfField

export pos_inf, neg_inf, nan, isposinf, isneginf, isnan, isnormal, max, min

################################################################################
#
#  Types and memory management
#
################################################################################

const ArfFieldID = ObjectIdDict()

mutable struct ArfField <: Field
  prec::Int
  rndmode::Cint

  function ArfField(p::Int = 256, r::Cint = Cint(4))
    try
      return ArfFieldID[p,r]::ArfField
    catch
      z = new(p,r)
      ArfFieldID[p,r] = z
      return z
    end
  end
end

mutable struct arf
  exp::Int # fmpz
  size::UInt64 # mp_size_t
  d1::Int64 # mantissa_struct
  d2::Int64
  parent::ArfField

  function arf()
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(i::Int)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_si, :libarb), Void, (Ptr{arf}, Int), &z, i)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(i::UInt)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_ui, :libarb), Void, (Ptr{arf}, UInt), &z, i)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(i::fmpz)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_fmpz, :libarb), Void, (Ptr{arf}, Ptr{fmpz}), &z, &i)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(x::BigFloat)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_mpfr, :libarb), Ptr{arf}, (Ptr{arf}, Ptr{BigFloat}), &z, &x)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(x::Cdouble)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_d, :libarb), Void, (Ptr{arf}, Cdouble), &z, x)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(x::mag)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_mag, :libarb), Void, (Ptr{arf}, Ptr{mag}), &z, &x)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(x::arf, p::Int, r::Cint)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_round, :libarb), Void,
                (Ptr{arf}, Ptr{arf}, Int, Cint), &z, &x, p, r)
    finalizer(z, _arf_clear_fn)
    return z
  end

  function arf(x::Int, p::Int, r::Cint)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_round_si, :libarb), Void,
                  (Ptr{arf}, Int, Int, Cint), &z, x, p, r)
    finalizer(z, _arf_clear_fn)
    return z
  end

 function arf(x::UInt, p::Int, r::Cint)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_round_ui, :libarb), Void,
                  (Ptr{arf}, UInt, Int, Cint), &z, x, p, r)
    finalizer(z, _arf_clear_fn)
    return z
  end

 function arf(x::fmpz, p::Int, r::Cint)
    z = new()
    ccall((:arf_init, :libarb), Void, (Ptr{arf}, ), &z)
    ccall((:arf_set_round_fmpz, :libarb), Void,
                  (Ptr{arf}, Ptr{fmpz}, Int, Cint), &z, &x, p, r)
    finalizer(z, _arf_clear_fn)
    return z
  end
end

function _arf_clear_fn(x::arf)
  ccall((:arf_clear, :libarb), Void, (Ptr{arf}, ), &x)
end

parent(x::arf) = x.parent

################################################################################
#
#  Parent object overloads
#
################################################################################

function (r::ArfField)()
  z = arf()
  z.parent = r
  return z
end

function (r::ArfField)(x::Int)
  z = arf(x, r.prec, r.rndmode)
  z.parent = r
  return z
end

function (r::ArfField)(x::UInt)
  z = arf(x, r.prec, r.rndmode)
  z.parent = r
  return z
end

function (r::ArfField)(x::fmpz)
  z = arf(x, r.prec, r.rndmode)
  z.parent = r
  return z
end

function (r::ArfField)(x::Float64)
  z = arf(x)
  z = arf(z, r.prec, r.rndmode)
  z.parent = r
  return z
end

function (r::ArfField)(x::BigFloat)
  z = arf(x)
  z = arf(z, r.prec, r.rndmode)
  z.parent = r
  return z
end

function (r::ArfField)(x::mag)
  z = arf(x)
  z = arf(z, r.prec, r.rndmode)
  z.parent = r
  return z
end

################################################################################
#
#  Special values
#
################################################################################

for (s,f) in (("zero", "arf_zero"), ("one", "arf_one"),
              ("pos_inf", "arf_pos_inf"), ("neg_inf", "arf_neg_inf"),
              ("nan", "arf_nan"))
  @eval begin
    function($(Symbol(s)))(r::ArfField)
      z = r()
      ccall(($f, :libarb), Void, (Ptr{arf}, ), &z)
      return z
    end
  end
end

################################################################################
#
#  String I/O
#
################################################################################

# this function is crap
function show(io::IO, a::arf)
  #cstr = ccall((:arf_get_str, :libarb), Ptr{UInt8}, (Ptr{arf}, ), a.data)
  #print(io, bytestring(cstr))
  #ccall((:flint_free, :libflint), Void, (Ptr{UInt8},), cstr)
  return show(io, BigFloat(a))
end

function show(io::IO, a::ArfField)
  print(io, "Field of arf's with precision $(a.prec)")
  print(io, " and rounding mode $(a.rndmode)")
end

################################################################################
#
#  Predicates
#
################################################################################

for (s,f) in (("iszero", "arf_iszero"), ("isone", "arf_is_one"),
              ("isposinf", "arf_is_pos_inf"),
              ("isneginf", "arf_is_neg_inf"),
              ("isinf", "arf_is_inf"),
              ("isnan", "arf_is_nan"),
              ("isnormal", "arf_is_normal"),
              ("isfinite", "arf_isfinite"),
              ("isspecial", "arf_is_special"))
  @eval begin
    function($(Symbol(s)))(x::arf)
      return Bool(ccall(($f, :libarb), Cint, (Ptr{arf},), x.data))
    end
  end
end

################################################################################
#
#  Comparison
#
################################################################################

function ==(x::arf, y::arf)
  r = ccall((:arf_equal, :libarb), Cint, (Ptr{arf}, Ptr{arf}), &x, &y)
  return Bool(r)
end

function compare(x::arf, y::arf)
  r = ccall((:arf_cmp, :libarb), Cint, (Ptr{arf}, Ptr{arf}), &x, &y)
  return r
end

function <(x::arf, y::arf)
  r = compare(x,y)
  r < 0 ? true : false
end

function >(x::arf, y::arf)
  r = compare(x,y)
  r > 0 ? true : false
end

function max(x::arf, y::arf)
  z = parent(x)()
  ccall((:arf_max, :libarb), Void, (Ptr{arf}, Ptr{arf}, Ptr{arf}), &z, &x, &y)
  return z
end

function min(x::arf, y::arf)
  z = parent(x)()
  ccall((:arf_min, :libarb), Void, (Ptr{arf}, Ptr{arf}, Ptr{arf}), &z, &x, &y)
  return z
end

################################################################################
#
#  Unary operations
#
################################################################################

function -(x::arf)
  z = parent(x)()
  ccall((:arf_neg, :libarb), Void, (Ptr{arf}, Ptr{arf}), &z, &x)
  return z
end

################################################################################
#
#  Binary operations
#
################################################################################

function +(x::arf, y::arf)
  check_parent(x,y)
  z = parent(x)()
  ccall(("arf_add", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

function +(x::arf, y::UInt)
  z = parent(x)()
  ccall(("arf_add_ui", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, UInt, Int, Cint),
              &z, &x, y, parent(x).prec, parent(x).rndmode)
  return z
end

+(x::UInt, y::arf) = +(y, x)

function +(x::arf, y::Int)
  z = parent(x)()
  ccall(("arf_add_si", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Int, Int, Cint),
              &z, &x, y, parent(x).prec, parent(x).rndmode)
  return z
end

+(x::Int, y::arf) = +(y, x)

function +(x::arf, y::fmpz)
  z = parent(x)()
  ccall(("arf_add_fmpz", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{fmpz}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

+(x::fmpz, y::arf) = +(y, x)

function *(x::arf, y::arf)
  z = parent(x)()
  ccall(("_arf_mul", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

function *(x::arf, y::UInt)
  z = parent(x)()
  ccall(("arf_mul_ui", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, UInt, Int, Cint),
              &z, &x, y, parent(x).prec, parent(x).rndmode)
  return z
end

*(x::UInt, y::arf) = *(y, x)

function *(x::arf, y::Int)
  z = parent(x)()
  ccall(("arf_mul_si", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Int, Int, Cint),
              &z, &x, y, parent(x).prec, parent(x).rndmode)
  return z
end

*(x::Int, y::arf) = *(y, x)

function *(x::arf, y::fmpz)
  z = parent(x)()
  ccall(("arf_mul_fmpz", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{fmpz}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

*(x::fmpz, y::arf) = *(y, x)

function -(x::arf, y::arf)
  check_parent(x,y)
  z = parent(x)()
  ccall(("arf_sub", :libarb), Void,
                (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
                &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

function -(x::arf, y::UInt)
  z = parent(x)()
  ccall(("arf_sub_ui", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, UInt, Int, Cint),
              &z, &x, y, parent(x).prec, parent(x).rndmode)
  return z
end

-(x::UInt, y::arf) = -(y - x)

function -(x::arf, y::Int)
  z = parent(x)()
  ccall(("arf_sub_si", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Int, Int, Cint),
              &z, &x, y, parent(x).prec, parent(x).rndmode)
  return z
end

-(x::Int, y::arf) = -(y - x)

function -(x::arf, y::fmpz)
  z = parent(x)()
  ccall(("arf_sub_fmpz", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{fmpz}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

-(x::fmpz, y::arf) = -(y - x)

function /(x::arf, y::arf)
  check_parent(x,y)
  z = parent(x)()
  ccall((:arf_div, :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
  return z
end

################################################################################
#
#  Unsafe arithmetic
#
################################################################################

function add!(x::arf, y::arf)
  ccall(("arf_add", :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
end

function div!(z::arf, x::arf, y::arf)
  ccall((:arf_div, :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
end

function sub!(x::arf, y::arf)
  ccall((:arf_sub, :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
end

function mul!(x::arf, y::arf)
  ccall((:_arf_mul, :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, &y, parent(x).prec, parent(x).rndmode)
end

################################################################################
#
#  Some functions
#
################################################################################

function sign(x::arf)
  r = ccall((:arf_sgn, :libarb), Cint, (Ptr{arf}, ), &x)
  return r
end

function abs(x::arf)
  z = parent(x)()
  ccall((:arf_abs, :libarb), Void, (Ptr{arf}, Ptr{arf}), &z, &x)
  return z
end

function Base.sqrt(x::arf)
  z = parent(x)()
  ccall((:arf_sqrt, :libarb), Void,
              (Ptr{arf}, Ptr{arf}, Int, Cint),
              &z, &x, parent(x).prec, parent(x).rndmode)
  return z
end

################################################################################
#
#  Conversion
#
################################################################################

function Cdouble(a::arf, rnd::Cint = 4)
  z = ccall((:arf_get_d, :libarb), Cdouble, (Ptr{arf}, Cint), a.data, rnd)
  return z
end

function BigFloat(x::arf)
  old_prec = get_bigfloat_precision()
  set_bigfloat_precision(parent(x).prec)
  z = BigFloat(0)
  r = ccall((:arf_get_mpfr, :libarb), Cint,
                (Ptr{BigFloat}, Ptr{arf}, Cint), &z, &x, Cint(0))
  set_bigfloat_precision(old_prec)
  return z
end


text:
# transition band is cosine. This allows satisfies the power complementary condition.
# cos( (yp-y)/(yp-ys) * π/2 ) for a falling edge. yp is pass band, ys is stop band.
function cosinetransitionlowpassfunction(y::T, yp::T, ys::T)::T where T
    abs_y = abs(y)

    if abs_y <= yp
        return one(T)
    elseif yp < abs_y <= ys
        #return cos( (yp-y)/(ys-yp)*π/2 )
        return cos( (yp-abs_y)/(ys-yp)*π/2 )
    end

    return zero(T)
end

function cosinetransitionhighpassfunction(x::T, xp::T, xs::T)::T where T
    abs_x = abs(x)

    if abs_x <= xs
        return zero(T)
    elseif xs < abs_x <= xp
        #return cos( (xp-x)/(xp-xs)*π/2 )
        return cos( (xp-abs_x)/(xp-xs)*π/2 )
    end

    return one(T)
end

# rs is rising edge's stopband.
# rp is rising edge's passband.
# fp is falling edge's passband.
# fs is falling edge's stopband.
function cosinetransitionbandpassfunction(y::T, rs, rp, fp, fs)::T where T
    abs_y = abs(y)

    # rising edge.
    if abs_y <= rs
        return zero(T)
    elseif rs < abs_y <= rp
        #return cos( (rp-y)/(rp-rs)*π/2 )
        return cos( (rp-abs_y)/(rp-rs)*π/2 )
    end

    # falling edge.
    if rp <= abs_y <= fp
        return one(T)
    elseif fp < abs_y <= fs
        #return cos( (fp-y)/(fs-fp)*π/2 )
        return cos( (fp-abs_y)/(fs-fp)*π/2 )
    end

    return zero(T)
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
"""
    AbstractChain{C, R}

A chain is a formal linear combination of `C`-type elements in `R`-module.

Given a `C`-type element set ``E``, we can construct a free `R`-module ``M`` that
has ``E ⊆ M`` as a basis ``M``, the module ``M`` is the module of the formal linear combinations
of elements of ``E``, or free module over ``E``, ``R^{(E)}``, s.t. given a finite subset
``{X_1, ..., X_n}`` of ``E``, a *formal linear combination* of ``X_1, ..., X_n`` is

``a_1 X_1 + ··· + a_n X_n``

where the ``a_i`` ∈ `R`.
"""
abstract type AbstractChain{C,R} end

# Chain interface
"""
    dim(ch::AbstractChain)

Returns the dimension of the chain `ch`.
"""
dim(ch::AbstractChain) = throw(MethodError(dim, (typeof(ch),)))

"""
    length(ch::AbstractChain)

Returns the length of the chain `ch`.
"""
length(ch::AbstractChain) = throw(MethodError(length, (typeof(ch),)))

"""
    keys(ch::AbstractChain)

Returns elements of the chain `ch`.
"""
keys(ch::AbstractChain) = throw(MethodError(keys, (typeof(ch),)))

"""
    values(ch::AbstractChain)

Returns coefficients of the chain `ch`.
"""
values(ch::AbstractChain) = throw(MethodError(values, (typeof(ch),)))

"""
    copy(ch::AbstractChain)

Returns the copy of the chain `ch`.
"""
copy(ch::AbstractChain) = throw(MethodError(copy, (typeof(ch),)))

"""
    setindex!(ch::AbstractChain{C,R}, v::R, e::C)

Sets an element `e` coefficient `v` in the chain `ch`.
"""
setindex!(ch::AbstractChain{C,R}, v::R, e::C) where {C,R} = throw(MethodError(setindex!, (typeof(ch), R, C)))

"""
    getindex(ch::AbstractChain{C,R}, e::C)

Returns a coefficient of the element `e` of the chain `ch`.
"""
getindex(ch::AbstractChain{C,R}, e::C) where {C,R} = throw(MethodError(getindex, (typeof(ch), C)))

"""
    push!(ch:: AbstractChain{C,R}, e::Pair{C,R})

Add element with coefficient `e` to the chain `ch`. If the element is already in the chain then the coefficient
value of `e` will be added to the chain element coefficient.

Use `setindex!` to reset an element coefficient in the chain.
"""
push!(ch:: AbstractChain{C,R}, e::Pair{C,R}) where {C,R} = throw(MethodError(push!, (typeof(ch),Tuple{C,R})))

"""
    iterate(ch::AbstractChain)

Return a element/coefficient iterator of the chain `ch`
"""
iterate(ch::AbstractChain, state...) = throw(MethodError(simplify, (typeof(ch),Any)))

"""
    simplify(ch::AbstractChain)

Returns a simplified (without 0 elements) copy of the chain `ch`.
"""
simplify(ch::AbstractChain) = throw(MethodError(simplify, (typeof(ch),)))

# Auxiliary methods

function show(io::IO, ch::AbstractChain)
    print(io, "$(dim(ch)): ")
    if iszero(ch)
        print(io, "0")
    else
        for (i,(id,v)) in enumerate(ch.cells)
            i != 1 && print(io, " + ")
            print(io, "$v[0x$(lpad(string(id,base=16),16,'0'))]")
        end
    end
end

iszero(ch::AbstractChain) = length(ch) == 0
keytype(ch::AbstractChain{C,R}) where {C,R} = C
valtype(ch::AbstractChain{C,R}) where {C,R} = R
push!(ch::AbstractChain{C,R}, e::C, c::R) where {C,R} = push!(ch, e=>c)
push!(ch::AbstractChain{C,R}, e::Tuple{C,R}) where {C,R} = push!(ch, e[1]=>e[2])

"""
    map(f, ch::AbstractChain)

Apply function `f` to every coefficient of the chain `ch`.
"""
map(f, ch::AbstractChain) = map!(f, copy(ch), ch)

function map!(f, dest::AbstractChain, src::AbstractChain)
    for (k,v) in src
        dest[k] = f(v)
    end
    return dest
end

"""
    in(e, ch::AbstractChain)

Check if the element `e` in the chain `ch`.
"""
in(e::C, ch::AbstractChain{C,R}) where {C,R} = in(e, keys(ch))

"""
    append!(a::AbstractChain, b::AbstractChain)

Perform mutable addition of the chain `b` to the chain `a`.
"""
function append!(a::AbstractChain, b::AbstractChain)
    for (k,v) in b
        if k ∈ a
            a[k] += v
        else
            push!(a, k=>v)
        end
    end
    return a
end

"""
    +(ch::AbstractChain{C,R}, e::Tuple{C,R})

Perform non-mutable addition of the element `e` with a coefficient to the chain `ch`.
"""
(+)(ch::AbstractChain{C,R}, e::Tuple{C,R}) where {C,R} = push!(copy(ch), e)

"""
    +(a::AbstractChain, b::AbstractChain)

Perform non-mutable addition of the chain `b` to the chain `a`.
"""
(+)(a::AbstractChain, b::AbstractChain) = append!(copy(a), b)

"""
    +(a::AbstractChain, b::AbstractChain)

Perform non-mutable addition of the chain `b` to the chain `a`.
"""
function (-)(a::AbstractChain, b::AbstractChain)
    c = copy(a)
    for (k,v) in b
        push!(c, k => -v)
    end
    return c
end

"""
    +(ch::AbstractChain{C,R}, r::R)

Perform non-mutable addition of the coefficient `r` to every element of the chain `ch`.
"""
(+)(ch::AbstractChain{C,R}, r::R) where {C,R} = map(e->e+r, ch)
(+)(r::R, ch::AbstractChain{C,R}) where {C,R} = ch+r

"""
    *(ch::AbstractChain{C,R}, r::R)

Perform non-mutable addition of the coefficient `r` to every element of the chain `ch`.
"""
(*)(ch::AbstractChain{C,R}, r::R) where {C,R} = map(e->e*r, ch)
(*)(r::R, ch::AbstractChain{C,R}) where {C,R} = ch*r


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Permuted congruential generator

"""
    mutable struct PCG <: AbstractRNG

Random number generator that implement the Permuted Congruential Generator, a simple fast space-efficient
statistically good algorithms for random number generation. See [O'Neill (2014)](https://www.pcg-random.org/paper.html).

# Fields

- `state::UInt64`: the state of the generator.
- `inc::UInt64`: sequence identifier.
"""
mutable struct PCG <: AbstractRNG
    state::UInt64
    inc::UInt64

    function PCG(state::UInt64 = UInt64(42), inc::UInt64 = UInt64(54))
        self = new(UInt64(0), (inc << UInt64(1)) | UInt64(1))
        rand(self)
        self.state += state
        rand(self)
        self
    end
end

@doc """
    PCG(state::UInt64 = UInt64(42), inc::UInt64 = UInt64(54))

Constructor for a [`PCG`](@ref) instance.

If no parameter is specified, the generated instance will have a state of 42 and a sequence identifier of 54.
""" PCG(::UInt64, ::UInt64)

function Base.rand(r::PCG)
    rand(r, Sampler(r, UInt32, Val(1)))
end

function Base.rand(r::PCG, ::Random.SamplerType{UInt32})
    oldstate = r.state
    r.state = UInt64(oldstate * UInt64(6364136223846793005) + r.inc)
    xorshifted = UInt32(((oldstate >> UInt64(18)) ⊻ oldstate) >> UInt64(27) & typemax(UInt32))
    rot = oldstate >> UInt64(59)
    UInt32((xorshifted >> rot) | (xorshifted << ((-rot) & UInt32(31))))
end

function Base.rand(r::PCG, ::Random.SamplerTrivial{Random.CloseOpen01{Float64}, Float64})
    convert(Float64, rand(r)/typemax(UInt32))
end

function Base.rand(r::PCG, ::Random.SamplerType{T}) where {T <: Integer}
    convert(T, rand(r) % typemax(T))
end

function Base.rand(r::PCG, ::Type{T}) where {T <: AbstractFloat}
    convert(T, rand(r, Sampler(r, Float64, Val(1))))
end


text:
using LIBLINEAR
using CTMCEnsemble
using StatsBase

function loaddata(file, nfeature)
    data = zeros(nfeature, 0)
    label = Int[]
    for l in eachline(file)
        items = split(l)
        push!(label, parse(Int, items[1]))
        f = zeros(nfeature)
        for i in items[2:end]
            j, v = split(i, ":")
            f[parse(Int, j)] = parse(Float64, v)
        end
        data = hcat(data, f)
    end

    (data, label)
end

function vote(p...)
    m = countmap([p...])
    collect(keys(m))[indmax(values(m))]
end

vote_all(preds) = vote.((getindex.(preds,1))...)

onehot(i, n, eps=0) = setindex!(ones(n) * eps / (n - 1), 1 - eps, i)

function pred_all(preds, method)
    results = Int[]
    n = size(preds[1][1], 1)
    for i = 1:n
        p = []
        for (pr, ran) in preds
            push!(p, (onehot(findfirst(ran, pr[i]), length(ran), 1e-5), ran))
        end
        push!(results, indmax(method(p)))
    end
    results
end

function cv(features, labels; nfold=10, C=1, bias=1)
    reorder = shuffle(1:length(labels))
    d = features[:, reorder]
    l = labels[reorder]
    nclass = maximum(labels)

    results = []
    voteacc = 0
    ctmcacc = 0
    meanacc = 0
    prodacc = 0
    for f = 1:nfold
        trainmask  = trues(length(l))
        trainmask[f:nfold:end] = false
        traindata  = d[:, trainmask]
        trainlabel = l[trainmask]
        testdata   = d[:, !trainmask]
        testlabel  = l[!trainmask]

        preds = []
        for i=1:nclass, j=i+1:nclass
            localcls = [i, j]
            trainmask = in.(trainlabel, [localcls])
            # testmask = in.(testlabel, [localcls])
            m = linear_train(trainlabel[trainmask], traindata[:, trainmask], C=C, bias=bias)
            p, _ = linear_predict(m, testdata)
            push!(preds, (p, localcls))
        end

        voteacc += sum(vote_all(preds) .== testlabel)
        ctmcacc += sum(pred_all(preds, ctmc) .== testlabel)
        meanacc += sum(pred_all(preds, average) .== testlabel)
        prodacc += sum(pred_all(preds, product) .== testlabel)
    end

    println("> C=$C, bias=$bias")
    println("vote: ", voteacc / length(l))
    println("ctmc: ", ctmcacc / length(l))
    println("mean: ", meanacc / length(l))
    println("prod: ", prodacc / length(l))
end


text:
# Leonardo Lozano, Daniel Duque, Andrés L. Medaglia (2016) An Exact Algorithm for the Elementary Shortest Path Problem with Resource Constraints. Transportation Science 50(1):348-357. https://doi.org/10.1287/trsc.2014.0582

function should_rollback(p::Label, pg::ESPPRC_Instance)
    # Section 4.3 Rollback Pruning
    if length(p.path) < 3
        return false
    end

    v_i = p.path[end-2]
    v_k = p.path[end-1]
    v_j = p.path[end]

    # p: ...... v_i -> v_k -> v_j
    # pp: ..... v_i -> v_j

    cost_p = pg.cost[v_i, v_k] + pg.cost[v_k, v_j]
    cost_pp = pg.cost[v_i, v_j]

    path_p = p.path
    path_pp = p.path[1:end-2]
    push!(path_pp, v_j)

    time_p = calculate_path_time(path_p, pg)
    time_pp = calculate_path_time(path_pp, pg)

    if cost_pp <= cost_p && time_pp <= time_p
        # dominated, should rollback
        return true
    else
        # non-dominated, no need to rollback
        return false
    end
end

function bounding_time_index(current_time::Float64, btimes::Vector{Float64})
    # time_values (sorted from greatest to smallest)
    Δ = btimes[1] - btimes[2]
    time_ub = btimes[1]
    k = Int(ceil((time_ub - current_time) / Δ)) + 1
    k = max(k, 1)
    # if k <= length(btimes)
    #     @assert current_time >= btimes[k]
    # end
    return k
end

function isbounded(p::Label, primal_bounds::Vector{Label}, lower_bounds::Matrix{Float64}, btimes::Vector{Float64}, pg::ESPPRC_Instance, bounding::Bool; root=pg.origin)
    if isempty(primal_bounds) || isempty(primal_bounds)
        @warn("isempty primal_bounds or primal_bounds")
        return false
    end

    v_i = p.path[end]
    bound = bounding ? min_lower_bound = minimum(lower_bounds[root, :]) : primal_bounds[root].cost

    # must use from the bound matrix B the lower closest value to τ available.
    # time_values (sorted from greatest to smallest)
    # τ <= p.time
    k = bounding_time_index(p.time, btimes)
    if k > length(btimes)
        # In this case, p.time < btimes[end]; no lower_bounds info available.
        return false
    elseif lower_bounds[v_i, k] < Inf
        # @assert btimes[k] <= p.time
        # The condition below should be strict inequality.
        # If it is set >=, then some problems cannot be solved optimally.
        bounded = p.cost + lower_bounds[v_i, k] > bound
        return bounded
    else
        return false
    end
end

function bounding_scheme!(btimes::Vector{Float64}, neg_cost_routes::Vector{Label}, pg::ESPPRC_Instance)
    n_nodes = length(pg.service_time)
    set_N = 1:n_nodes

    primal_bounds = initialize_label.(set_N, n_nodes; cost=Inf)
    # best_labels = Array{Label, 1}(undef, n_nodes)
    # for i in eachindex(best_labels)
    #     best_labels[i] = initialize_label(i, n_nodes; cost=Inf)
    # end

    # lower_bounds = Matrix{Label}(undef, n_nodes, length(btimes))
    # for i in 1:n_nodes, k in 1:length(btimes)
    #     lower_bounds[i, k] = initialize_label(i, n_nodes; cost=Inf)
    # end
    lower_bounds = fill(Inf, n_nodes, length(btimes))

    bounding_iteration = 0
    for k in 1:length(btimes)
        cur_time = btimes[k]

        for v_i in 1:n_nodes
            if v_i != pg.destination
                p = initialize_label(v_i, n_nodes)
                if cur_time <= pg.late_time[v_i]
                    p.time = max(cur_time, pg.early_time[v_i])
                    pulse_procedure!(p, primal_bounds, lower_bounds, btimes, neg_cost_routes, pg;
                                        init_time=cur_time, bounding=true, root=v_i)
                end
                bounding_iteration += 1
            end
        end
    end

    # @show btimes
    # for i in 1:n_nodes, k in 1:length(btimes)
    #     @show k, i, btimes[k], lower_bounds[i, k].cost
    # end

    # @show bounding_iteration

    return primal_bounds, lower_bounds
end # bounding_scheme

function pulse_procedure!(p::Label, primal_bounds::Vector{Label}, lower_bounds::Matrix{Float64}, btimes::Vector{Float64}, neg_cost_routes::Vector{Label}, pg::ESPPRC_Instance; init_time=0.0, bounding=false, root=pg.origin)
    if length(neg_cost_routes) >= pg.max_neg_routes
        return
    end

    # current node
    v_i = p.path[end]
    if p.time > pg.late_time[v_i]
        return
    end
    p.time = max(p.time, pg.early_time[v_i])

    if v_i == pg.destination   # Arrived at the destination. Update the best bounds
        if p.cost < primal_bounds[root].cost
            primal_bounds[root] = Label(p.time, p.load, copy(p.flag), p.cost, copy(p.path))
            # deepcopy(p)
        end
        if bounding
            k = bounding_time_index(init_time, btimes)
            if k <= length(btimes)
                if p.cost < lower_bounds[root, k]
                    for kk in k:length(btimes)
                        # lower_bounds[root, kk] = deepcopy(p)
                        lower_bounds[root, kk] = p.cost
                    end
                end
            end
        end
        if root == pg.origin && p.cost < 0
            push!(neg_cost_routes, p)
        end
        return
    end

    # If it didn't arrive at the destination yet
    if isbounded(p, primal_bounds, lower_bounds, btimes, pg, bounding, root=root)
        return
    end
    if should_rollback(p, pg)
        return
    end

    # Create a new pulse
    for v_j in pg.forward_star[v_i]
        if p.flag[v_j] == 0
            isreachable, new_time, new_load = forward_reach(p, v_i, v_j, pg)
            if isreachable
                global counter += 1

                pp = Label(new_time, new_load, copy(p.flag), p.cost + pg.cost[v_i, v_j], copy(p.path))
                push!(pp.path, v_j)
                pp.flag[v_j] = 1

                # pp = deepcopy(p)
                # push!(pp.path, v_j)
                # pp.cost += pg.cost[v_i, v_j]
                # pp.load = new_load
                # pp.time = new_time
                # pp.flag[v_j] = 1

                pulse_procedure!(pp, primal_bounds, lower_bounds, btimes, neg_cost_routes, pg;
                            init_time=init_time, bounding=bounding, root=root)
            end
        end
    end

    return
end

function solveESPPRCpulse(org_pg::ESPPRC_Instance; step=-1, max_neg_routes=MAX_INT::Int)
    pg = deepcopy(org_pg)
    graph_reduction!(pg)
    pg.max_neg_routes = max_neg_routes

    n_nodes = length(pg.late_time)

    # Discrete Time Steps for Bounding
    time_ub = calculate_max_T(pg)
    time_lb = 0.1 * time_ub
    if step == -1
        step = Int(floor(time_ub*0.9 / 10)) + 1
    end
    btimes = collect(time_ub-step : -step : time_lb)

    # Container for negative reduced cost solutions
    neg_cost_routes = Vector{Label}(undef,0)

    # Bounding Scheme
    primal_bounds, lower_bounds = bounding_scheme!(btimes, neg_cost_routes, pg)

    p = initialize_label(pg.origin, n_nodes)
    pulse_procedure!(p, primal_bounds, lower_bounds, btimes, neg_cost_routes, pg)

    # @show counter
    if max_neg_routes < MAX_INT && !isempty(neg_cost_routes)
        return find_best_label!(neg_cost_routes), neg_cost_routes
    else
        return primal_bounds[pg.origin], neg_cost_routes
    end
end


text:
#@require UnicodePlots begin
"""
Description

      The uplot() method plots a time series in REPL, using the UnicodePlots package.

Usage

      using MarketData
      uplot(cl)
      uplot(ohlc, 4) # the fourth column ("Close") will be plotted

Method Signature(s)

      uplot(t::TimeSeries.TimeArray, n::Int=1)

Details

      UnicodePlots provides simple plots inside of REPL. The uplot() method is a wrapper for the UnicodePlots.lineplot()
      method that takes a TimeSeries.TimeArray as the first argument. The second argument defines the column to be
      plotted and defaults to n=1.

References

      https://github.com/Evizero/UnicodePlots.jl

See Also

      uhist() which plots a histogram of a time series
      plot() which plots a time series in a separate window
"""

function uplot(t::TimeSeries.TimeArray, n::Int=1)
    s = t.timestamp
    v = t.values[:,n]
    c = t.colnames[n]
    isa(t.meta, AbstractString) ? m = t.meta :  m = ""
    lineplot(s, v, title=string(m, " ", c))
end
#end # of requires macro block


text:
```
The gas laws to calculate volume by
temperature and/or pressure changes respectively :
Charles, Boyles, Gay Lussacs and Avogadros laws
```

"""
Calculate the volume by using charles law, stating that gas volume
increases or decreases directly relative to change in temprature.

curentVolume = Curent c.c of gas
curenttemperature = The Curent temperature in either C or F
newtemperature = temperature the gas is being heated, or cooled to
temperatureUnit = (c)entigrade, (f)ahrenheit, or (a)bsolute
"""
function charles_law(curentVolume, curenttemperature, newtemperature, temperatureUnit)
    newVolume = 0;
    curenttemperature = absolute_temperature(curenttemperature, temperatureUnit)
    newtemperature = absolute_temperature(newtemperature, temperatureUnit)
    newVolume = curentVolume * (newtemperature / curenttemperature)
    return newVolume
end

"""
When temperature remains contstant but pressure changes we can calculate new volume by using
Boyles law when solving for pressure.
Which simlpy states : New volume = old volume * old pressure devided by new pressure

curentPressure = The curent pressure in mm
curentVolume = The curent volume of the gas
desiredPressure = The pressure we want to be at
"""
function boyles_law_solve_pressure(curentPressure, curentVolume, desiredPressure)
    newVolume = 0;
    newVolume = (curentVolume * curentPressure) / desiredPressure
    return newVolume
end

"""
When temperature remains contstant but volume changes we can calculate new pressure by using
Boyles law when solving for volume.
Which simlpy states : New pressure = old pressure * curent volume devided by new volume

curentPressure = The curent pressure in mm
curentVolume = The curent volume of the gas
desiredVolume = The volume we want to have
"""
function boyles_law_solve_volume(curentPressure, curentVolume, desiredVolume)
    newPressure = 0
    newPressure = (curentPressure * curentVolume) / desiredVolume
    return newPressure
end

"""
The pressure of a gas of fixed mass and fixed volume is directly
proprtional to the gas's absolute temperature.

If a gas's temperature increases, so does its pressure if the mass and volume are
held constant.

pressure = Curent pressure in atm
temperature = Curent temperature of container
newTemp = temperature at which we want to know the pressure
temperatureUnit = (c)entigrade, (f)ahrenheit, or (a)bsolute
"""
function gay_lussacs_law(pressure, temperature, newTemp, temperatureUnit)
    newTemp = absolute_temperature(newTemp, temperatureUnit)
    temperature = absolute_temperature(temperature, temperatureUnit)
    newPressure = (pressure * newTemp) / temperature
    return newPressure
end

"""
Avrogados law under the assumption of ideal gases teaches us that
any gas, while having the same temperature, pressure and volume will have the
same number of molecules. or : V1|N1 = V2|N2
"""
function avogadros_law()
    return null
end

"""
One of the use cases of grahams law
is that we can calculate the ratio of diffusion rates of two gasses.

the rate of effusion or diffusion of a gass in inversely proprtional to the square
root of the molar mass of the gas.

using g/mol give the molar masses of both gasses.
"""
function grahams_law_velocity_ratio(molar_mass_a,molar_mass_b)
        ratio = sqrt(molar_mass_a/molar_mass_b)
    return ratio
end


text:
module Core2CoreLatency

using ..ThreadPinning: pinthread, pinthreads, getcpuids, @tspawnat

const State = Int
const Preparing = 0
const Ready = 1
const Ping = 2
const Pong = 3
const Finish = 4

Base.@kwdef struct Sync
    state::Threads.Atomic{State} = Threads.Atomic{State}(Preparing)
end

state(S::Sync) = S.state[]
function wait_until(S::Sync, expected_state::State)
    while state(S) != expected_state
    end
    return nothing
end
function set(S::Sync, state::State)
    S.state[] = state
    return nothing
end
function wait_as_long_as(S::Sync, wait_state::State)
    loaded_state = state(S)
    while loaded_state == wait_state
        loaded_state = state(S)
    end
    return loaded_state
end

function run_bench(cpu1::Integer, cpu2::Integer; nsamples::Integer=100, mode::Symbol=:min)
    cpu1 == cpu2 && return zero(Float64)
    Threads.nthreads() >= 2 || @error("Need at least two Julia threads.")

    S = Sync()
    pinthread(cpu1)

    second_thread = @tspawnat 2 begin
        pinthread(cpu2)
        set(S, Ready)

        state = wait_as_long_as(S, Ready)
        while state != Finish
            if state == Ping
                set(S, Pong)
                state = wait_as_long_as(S, Pong)
            end
        end
    end

    wait_until(S, Ready)

    Δts = zeros(typeof(time_ns()), nsamples)
    @inbounds for i in 1:nsamples
        Δts[i] = begin
            t = time_ns()
            set(S, Ping)
            wait_until(S, Pong)
            time_ns() - t
        end
    end

    if mode == :avg
        Δt = sum(Float64, Δts) / nsamples
    elseif mode == :min || mode == :minimum
        Δt = Float64(minimum(Δts))
    else
        throw(ArgumentError("Unkown mode $mode."))
    end

    set(S, Finish)
    fetch(second_thread)
    return Δt
end

end # module


text:
get_inner_vars_count(::PSY.DynamicGenerator) = GEN_INNER_VARS_SIZE
get_inner_vars_count(::PSY.DynamicInverter) = INV_INNER_VARS_SIZE
get_inner_vars_count(::PSY.PeriodicVariableSource) = 0
get_inner_vars_count(::PSY.SingleCageInductionMachine) = 0
get_inner_vars_count(::PSY.SimplifiedSingleCageInductionMachine) = 0

index(::Type{<:PSY.TurbineGov}) = 1
index(::Type{<:PSY.PSS}) = 2
index(::Type{<:PSY.Machine}) = 3
index(::Type{<:PSY.Shaft}) = 4
index(::Type{<:PSY.AVR}) = 5

index(::Type{<:PSY.DCSource}) = 1
index(::Type{<:PSY.FrequencyEstimator}) = 2
index(::Type{<:PSY.OuterControl}) = 3
index(::Type{<:PSY.InnerControl}) = 4
index(::Type{<:PSY.Converter}) = 5
index(::Type{<:PSY.Filter}) = 6

"""
Wraps DynamicInjection devices from PowerSystems to handle changes in controls and connection
status, and allocate the required indexes of the state space.
"""
struct DynamicWrapper{T <: PSY.DynamicInjection}
    device::T
    system_base_power::Float64
    system_base_frequency::Float64
    static_type::Type{<:PSY.StaticInjection}
    bus_category::Type{<:BusCategory}
    connection_status::Base.RefValue{Float64}
    V_ref::Base.RefValue{Float64}
    ω_ref::Base.RefValue{Float64}
    P_ref::Base.RefValue{Float64}
    Q_ref::Base.RefValue{Float64}
    inner_vars_index::Vector{Int}
    ix_range::Vector{Int}
    ode_range::Vector{Int}
    bus_ix::Int
    global_index::Base.ImmutableDict{Symbol, Int}
    component_state_mapping::Base.ImmutableDict{Int, Vector{Int}}
    input_port_mapping::Base.ImmutableDict{Int, Vector{Int}}
end

function state_port_mappings(
    dynamic_device::T,
    device_states,
) where {T <: Union{PSY.DynamicGenerator, PSY.DynamicInverter}}
    component_state_mapping = Dict{Int, Vector{Int}}()
    input_port_mapping = Dict{Int, Vector{Int}}()
    for c in PSY.get_dynamic_components(dynamic_device)
        ix = index(typeof(c))
        component_state_mapping[ix] = _index_local_states(c, device_states)
        input_port_mapping[ix] = _index_port_mapping!(c, device_states)
    end
    return (component_state_mapping, input_port_mapping)
end

function state_port_mappings(dynamic_device::PSY.DynamicInjection, device_states)
    component_state_mapping = Dict{Int, Vector{Int}}()
    input_port_mapping = Dict{Int, Vector{Int}}()
    return (component_state_mapping, input_port_mapping)
end

function DynamicWrapper(
    device::T,
    bus_ix::Int,
    ix_range,
    ode_range,
    inner_var_range,
    sys_base_power,
    sys_base_freq,
) where {T <: PSY.Device}
    dynamic_device = PSY.get_dynamic_injector(device)
    @assert dynamic_device !== nothing
    device_states = PSY.get_states(dynamic_device)

    component_state_mapping, input_port_mapping =
        state_port_mappings(dynamic_device, device_states)

    return DynamicWrapper{typeof(dynamic_device)}(
        dynamic_device,
        sys_base_power,
        sys_base_freq,
        T,
        BUS_MAP[PSY.get_bustype(PSY.get_bus(device))],
        Base.Ref(1.0),
        Base.Ref(PSY.get_V_ref(dynamic_device)),
        Base.Ref(PSY.get_ω_ref(dynamic_device)),
        Base.Ref(PSY.get_P_ref(dynamic_device)),
        Base.Ref(PSY.get_reactive_power(device)),
        inner_var_range,
        ix_range,
        ode_range,
        bus_ix,
        Base.ImmutableDict(
            sort!(device_states .=> ix_range, by = x -> x.second, rev = true)...,
        ),
        isempty(component_state_mapping) ? Base.ImmutableDict{Int, Vector{Int}}() :
        Base.ImmutableDict(component_state_mapping...),
        isempty(input_port_mapping) ? Base.ImmutableDict{Int, Vector{Int}}() :
        Base.ImmutableDict(input_port_mapping...),
    )
end

function DynamicWrapper(
    device::PSY.Source,
    bus_ix::Int,
    ix_range,
    ode_range,
    inner_var_range,
    sys_base_power,
    sys_base_freq,
)
    dynamic_device = PSY.get_dynamic_injector(device)
    @assert dynamic_device !== nothing
    device_states = PSY.get_states(dynamic_device)
    IS.@assert_op PSY.get_X_th(dynamic_device) == PSY.get_X_th(device)
    IS.@assert_op PSY.get_R_th(dynamic_device) == PSY.get_R_th(device)

    return DynamicWrapper{typeof(dynamic_device)}(
        dynamic_device,
        sys_base_power,
        sys_base_freq,
        PSY.Source,
        BUS_MAP[PSY.get_bustype(PSY.get_bus(device))],
        Base.Ref(1.0),
        Base.Ref(0.0),
        Base.Ref(0.0),
        Base.Ref(0.0),
        Base.Ref(0.0),
        collect(inner_var_range),
        collect(ix_range),
        collect(ode_range),
        bus_ix,
        Base.ImmutableDict(Dict(device_states .=> ix_range)...),
        Base.ImmutableDict{Int, Vector{Int}}(),
        Base.ImmutableDict{Int, Vector{Int}}(),
    )
end

function _index_local_states(component::PSY.DynamicComponent, device_states::Vector{Symbol})
    component_state_index = Vector{Int}(undef, PSY.get_n_states(component))
    component_states = PSY.get_states(component)
    for (ix, s) in enumerate(component_states)
        component_state_index[ix] = findfirst(x -> x == s, device_states)
    end
    return component_state_index
end

function _index_port_mapping!(
    component::PSY.DynamicComponent,
    device_states::Vector{Symbol},
)
    ports = Ports(component)
    index_component_inputs = Vector{Int}()
    for i in ports[:states]
        tmp = [(ix, var) for (ix, var) in enumerate(device_states) if var == i]
        isempty(tmp) && continue
        push!(index_component_inputs, tmp[1][1])
    end
    return index_component_inputs
end

get_device(wrapper::DynamicWrapper) = wrapper.device
get_device_type(::DynamicWrapper{T}) where {T <: PSY.DynamicInjection} = T
get_bus_category(wrapper::DynamicWrapper) = wrapper.bus_category
get_inner_vars_index(wrapper::DynamicWrapper) = wrapper.inner_vars_index
get_ix_range(wrapper::DynamicWrapper) = wrapper.ix_range
get_ode_ouput_range(wrapper::DynamicWrapper) = wrapper.ode_range
get_bus_ix(wrapper::DynamicWrapper) = wrapper.bus_ix
get_global_index(wrapper::DynamicWrapper) = wrapper.global_index
get_component_state_mapping(wrapper::DynamicWrapper) = wrapper.component_state_mapping
get_input_port_mapping(wrapper::DynamicWrapper) = wrapper.input_port_mapping

get_system_base_power(wrapper::DynamicWrapper) = wrapper.system_base_power
get_system_base_frequency(wrapper::DynamicWrapper) = wrapper.system_base_frequency

get_P_ref(wrapper::DynamicWrapper) = wrapper.P_ref[]
get_Q_ref(wrapper::DynamicWrapper) = wrapper.Q_ref[]
get_V_ref(wrapper::DynamicWrapper) = wrapper.V_ref[]
get_ω_ref(wrapper::DynamicWrapper) = wrapper.ω_ref[]

set_P_ref(wrapper::DynamicWrapper, val::Float64) = wrapper.P_ref[] = val
set_Q_ref(wrapper::DynamicWrapper, val::Float64) = wrapper.Q_ref[] = val
set_V_ref(wrapper::DynamicWrapper, val::Float64) = wrapper.V_ref[] = val
set_ω_ref(wrapper::DynamicWrapper, val::Float64) = wrapper.ω_ref[] = val

# PSY overloads for the wrapper
PSY.get_name(wrapper::DynamicWrapper) = PSY.get_name(wrapper.device)
PSY.get_ext(wrapper::DynamicWrapper) = PSY.get_ext(wrapper.device)
PSY.get_states(wrapper::DynamicWrapper) = PSY.get_states(wrapper.device)
PSY.get_n_states(wrapper::DynamicWrapper) = PSY.get_n_states(wrapper.device)
PSY.get_base_power(wrapper::DynamicWrapper) = PSY.get_base_power(wrapper.device)

PSY.get_machine(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator} =
    wrapper.device.machine
PSY.get_shaft(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator} =
    wrapper.device.shaft
PSY.get_avr(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator} =
    wrapper.device.avr
PSY.get_prime_mover(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator} =
    wrapper.device.prime_mover
PSY.get_pss(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicGenerator} =
    wrapper.device.pss

PSY.get_converter(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter} =
    wrapper.device.converter
PSY.get_outer_control(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter} =
    wrapper.device.outer_control
PSY.get_inner_control(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter} =
    wrapper.device.inner_control
PSY.get_dc_source(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter} =
    wrapper.device.dc_source
PSY.get_freq_estimator(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter} =
    wrapper.device.freq_estimator
PSY.get_filter(wrapper::DynamicWrapper{T}) where {T <: PSY.DynamicInverter} =
    wrapper.device.filter

# PSY overloads of specific Dynamic Injectors

PSY.get_V_ref(wrapper::PSY.SingleCageInductionMachine) = 1.0
PSY.get_V_ref(wrapper::PSY.SimplifiedSingleCageInductionMachine) = 1.0
PSY.get_P_ref(wrapper::PSY.SingleCageInductionMachine) = PSY.get_τ_ref(wrapper)
PSY.get_P_ref(wrapper::PSY.SimplifiedSingleCageInductionMachine) = PSY.get_τ_ref(wrapper)
PSY.get_ω_ref(wrapper::PSY.SingleCageInductionMachine) = 1.0
PSY.get_ω_ref(wrapper::PSY.SimplifiedSingleCageInductionMachine) = 1.0

function get_local_state_ix(
    wrapper::DynamicWrapper,
    ::Type{T},
) where {T <: PSY.DynamicComponent}
    return wrapper.component_state_mapping[index(T)]
end

function get_input_port_ix(
    wrapper::DynamicWrapper,
    ::Type{T},
) where {T <: PSY.DynamicComponent}
    return wrapper.input_port_mapping[index(T)]
end

function get_local_state_ix(wrapper::DynamicWrapper, ::T) where {T <: PSY.DynamicComponent}
    return get_local_state_ix(wrapper, T)
end

function get_input_port_ix(wrapper::DynamicWrapper, ::T) where {T <: PSY.DynamicComponent}
    return get_input_port_ix(wrapper, T)
end

struct StaticWrapper{T <: PSY.StaticInjection, V}
    device::T
    connection_status::Base.RefValue{Float64}
    V_ref::Base.RefValue{Float64}
    θ_ref::Base.RefValue{Float64}
    P_ref::Base.RefValue{Float64}
    Q_ref::Base.RefValue{Float64}
    bus_ix::Int
end

function DynamicWrapper(device::T, bus_ix::Int) where {T <: PSY.Device}
    bus = PSY.get_bus(device)
    StaticWrapper{T, BUS_MAP[PSY.get_bustype(bus)]}(
        device,
        Base.Ref(1.0),
        Base.Ref(PSY.get_magnitude(bus)),
        Base.Ref(PSY.get_angle(bus)),
        Base.Ref(PSY.get_active_power(device)),
        Base.Ref(PSY.get_reactive_power(device)),
        bus_ix,
    )
end

function StaticWrapper(device::T, bus_ix::Int) where {T <: PSY.Source}
    bus = PSY.get_bus(device)
    return StaticWrapper{T, BUS_MAP[PSY.get_bustype(bus)]}(
        device,
        Base.Ref(1.0),
        Base.Ref(PSY.get_internal_voltage(device)),
        Base.Ref(PSY.get_internal_angle(device)),
        Base.Ref(PSY.get_active_power(device)),
        Base.Ref(PSY.get_reactive_power(device)),
        bus_ix,
    )
end

# get_bustype is already exported in PSY. So this is named this way to avoid name collisions
get_bus_category(::StaticWrapper{<:PSY.StaticInjection, U}) where {U} = U

# TODO: something smart to forward fields
get_device(wrapper::StaticWrapper) = wrapper.device
get_bus_ix(wrapper::StaticWrapper) = wrapper.bus_ix

get_P_ref(wrapper::StaticWrapper) = wrapper.P_ref[]
get_Q_ref(wrapper::StaticWrapper) = wrapper.Q_ref[]
get_V_ref(wrapper::StaticWrapper) = wrapper.V_ref[]
get_θ_ref(wrapper::StaticWrapper) = wrapper.θ_ref[]

set_P_ref(wrapper::StaticWrapper, val::Float64) = wrapper.P_ref[] = val
set_Q_ref(wrapper::StaticWrapper, val::Float64) = wrapper.Q_ref[] = val
set_V_ref(wrapper::StaticWrapper, val::Float64) = wrapper.V_ref[] = val
set_θ_ref(wrapper::StaticWrapper, val::Float64) = wrapper.θ_ref[] = val

PSY.get_bus(wrapper::StaticWrapper) = PSY.get_bus(wrapper.device)
PSY.get_active_power(wrapper::StaticWrapper) = PSY.get_active_power(wrapper.device)
PSY.get_reactive_power(wrapper::StaticWrapper) = PSY.get_reactive_power(wrapper.device)
PSY.get_name(wrapper::StaticWrapper) = PSY.get_name(wrapper.device)
PSY.get_ext(wrapper::StaticWrapper) = PSY.get_ext(wrapper.device)

mutable struct ExpLoadParams
    P_exp::Float64
    P_coeff::Float64
    Q_exp::Float64
    Q_coeff::Float64
end

mutable struct StaticLoadWrapper
    bus::PSY.Bus
    V_ref::Float64
    θ_ref::Float64
    P_power::Float64
    P_current::Float64
    P_impedance::Float64
    Q_power::Float64
    Q_current::Float64
    Q_impedance::Float64
    exp_params::Vector{ExpLoadParams}
    exp_names::Dict{String, Int}
    bus_ix::Int
end

function StaticLoadWrapper(bus::PSY.Bus, loads::Vector{PSY.ElectricLoad}, bus_ix::Int)
    P_power = 0.0
    P_current = 0.0
    P_impedance = 0.0
    Q_power = 0.0
    Q_current = 0.0
    Q_impedance = 0.0

    # Add ZIP Loads
    for ld in loads
        if isa(ld, PSY.PowerLoad)
            if PSY.get_available(ld) &&
               PSY.get_model(ld) == PSY.LoadModels.ConstantImpedance
                P_impedance += PSY.get_active_power(ld)
                Q_impedance += PSY.get_reactive_power(ld)
            elseif PSY.get_available(ld) &&
                   PSY.get_model(ld) == PSY.LoadModels.ConstantCurrent
                P_current += PSY.get_active_power(ld)
                Q_current += PSY.get_reactive_power(ld)
            elseif PSY.get_available(ld) &&
                   PSY.get_model(ld) == PSY.LoadModels.ConstantPower
                P_power += PSY.get_active_power(ld)
                Q_power += PSY.get_reactive_power(ld)
            else
                error(
                    "Not supported load model in $(PSY.get_number(bus)) named $(PSY.get_name(ld))",
                )
            end
        end
    end

    # Add Exponential Loads
    exp_loads = filter(x -> isa(x, PSY.ExponentialLoad) && PSY.get_available(x), loads)
    exp_params = Vector{ExpLoadParams}(undef, length(exp_loads))
    dict_names = Dict{String, Int}()
    if !isempty(exp_loads)
        for (ix, ld) in enumerate(exp_loads)
            dict_names[PSY.get_name(ld)] = ix
            exp_params[ix] = ExpLoadParams(
                PSY.get_active_power(ld),
                PSY.get_active_power_coefficient(ld),
                PSY.get_reactive_power(ld),
                PSY.get_reactive_power_coefficient(ld),
            )
        end
    end

    return StaticLoadWrapper(
        bus,
        PSY.get_magnitude(bus),
        PSY.get_angle(bus),
        P_power,
        P_current,
        P_impedance,
        Q_power,
        Q_current,
        Q_impedance,
        exp_params,
        dict_names,
        bus_ix,
    )
end

PSY.get_bus(wrapper::StaticLoadWrapper) = wrapper.bus
PSY.get_name(wrapper::StaticLoadWrapper) = PSY.get_name(wrapper.bus)

get_V_ref(wrapper::StaticLoadWrapper) = wrapper.V_ref
get_θ_ref(wrapper::StaticLoadWrapper) = wrapper.θ_ref
get_P_power(wrapper::StaticLoadWrapper) = wrapper.P_power
get_P_current(wrapper::StaticLoadWrapper) = wrapper.P_current
get_P_impedance(wrapper::StaticLoadWrapper) = wrapper.P_impedance
get_Q_power(wrapper::StaticLoadWrapper) = wrapper.Q_power
get_Q_current(wrapper::StaticLoadWrapper) = wrapper.Q_current
get_Q_impedance(wrapper::StaticLoadWrapper) = wrapper.Q_impedance
get_exp_params(wrapper::StaticLoadWrapper) = wrapper.exp_params
get_exp_names(wrapper::StaticLoadWrapper) = wrapper.exp_names
get_bus_ix(wrapper::StaticLoadWrapper) = wrapper.bus_ix

set_V_ref!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.V_ref = val
set_θ_ref!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.θ_ref = val
set_P_power!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.P_power = val
set_P_current!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.P_current = val
set_P_impedance!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.P_impedance = val
set_Q_power!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.Q_power = val
set_Q_current!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.Q_current = val
set_Q_impedance!(wrapper::StaticLoadWrapper, val::Float64) = wrapper.Q_impedance = val

function set_connection_status(wrapper::Union{StaticWrapper, DynamicWrapper}, val::Int)
    if val == 0
        @debug "Generator $(PSY.get_name(wrapper)) status set to off"
    elseif val == 1
        @debug "Generator $(PSY.get_name(wrapper)) status set to on"
    else
        error("Invalid status $val. It can only take values 1 or 0")
    end
    wrapper.connection_status[] = Float64(val)
end

get_connection_status(wrapper::Union{StaticWrapper, DynamicWrapper}) =
    wrapper.connection_status[]


text:
function sprintfrange{T<:Integer}(a::Array{T,1})
    len = length(a)
    0 < len || return ""
    dropme = falses(len)
    dropme[2:end-1] = Bool[a[i-1]==a[i]-1 && a[i+1]==a[i]+1 for i in 2:(len-1)]
    s = [string(i) for i in a]
    s[dropme] = "X"
    s = join(s, ",")
    replace(s, r",[,X]+,", "-")
end

testa = [ 0,  1,  2,  4,  6,  7,  8, 11, 12, 14,
         15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
         25, 27, 28, 29, 30, 31, 32, 33, 35, 36,
         37, 38, 39]

println("Testing range-style formatting.")
println("   ", testa, "\n       =>\n   ", sprintfrange(testa))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module BenchStructureFields
using BenchmarkTools
using Statistics
using Gadfly
using MLStyle
using DataFrames
using MacroTools: @capture
using ..ArbitrarySampler
using ..Utils

lineno = @spec LineNumberNode(::Int isa (x -> mod(x, 100) + 1))
field = @spec :($(::Symbol)::$(::Symbol)) || ::Symbol

specs = [
    Symbol("N(fields)=2") => @spec(quote
        $$lineno
        struct $(::Symbol)
            $([$lineno{0, 2}..., $field{2}...]...)
        end
    end),
    Symbol("N(fields)=5") => @spec(quote
        $$lineno
        struct $(::Symbol)
            $([$lineno{0,2}..., $field{5}...]...)
        end
    end),
    Symbol("N(fields)=10") => @spec(quote
        $$lineno
        struct $(::Symbol)
            $([$lineno{0,2}..., $field{10}...]...)
        end
    end),
    Symbol("N(fields)=50") => @spec(quote
        $$lineno
        struct $(::Symbol)
            $([$lineno{0,2}..., $field{50}...]...)
        end
    end),
]

struct StructField end
@active StructField(x) begin
    @match x begin
        :($a::$b) => Some((a, b))
        a::Symbol => Some((a, Any))
        _ => nothing
    end
end

using MLStyle.AbstractPatterns: effect
struct PushTo end
function MLStyle.pattern_uncall(::Type{PushTo}, self, _, _, args)
    length(args) === 1 || error("PushTo accepts 1 arg.")
    container = args[1]
    effect() do target, _, _
        :(push!($container, $target))
    end
end

implementations = [
    :MLStyle => function (ex)
        fields = Tuple[]
        @match ex begin
            quote
                $(::LineNumberNode)
                struct $typename
                    $(Many[::LineNumberNode||StructField(PushTo(fields))]...)
                end
            end => (typename, fields)
        end
    end,
    :MacroTools => function (ex)
        @capture(ex, struct T_
            fields__::Any
        end)
        (T, fields)
    end,
]

records = NamedTuple{(:time_mean, :implementation, :case)}[]
for (spec_id, spec) in specs
    # group_key = string(spec_id)
    # suite[group_key] = BenchmarkGroup()
    for (impl_id, impl_fn) in implementations
        bench′ =
            @benchmark $impl_fn(sample) setup = (sample = $generate($spec)) samples = 2000
        time′ = mean(bench′.times)
        @info :bench (spec_id, impl_id, time′)
        push!(records, (time_mean = time′, implementation = impl_id, case = spec_id))
    end
end

df = DataFrame(records)

@info df

theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size = 5px,
)

report_meantime, df_time = report(
    df, Guide.title("Extracting Expr(:struct, ...)");
    benchfield = :time_mean
)

open("stats/bench-structfields.txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/bench-structfields.svg", 14inch, 6inch), report_meantime)

end


text:
# Exports

export Shape, PostOpObj, create, combine, unpack, repack, repack!

abstract type Shape end

mutable struct PostOpObj
    objs::Vector{Any}
    operations::Array{Function}
    function PostOpObj(obj, func)
        if ~isa(obj, Vector)
            obj = [obj]
        end
        if ~isa(func, Vector)
                func = [func]
        end
        new(obj, func)
    end
end

function Base.copy(x::T) where T<:Union{Shape, PostOpObj}
    T(copy.([getfield(x, fn) for fn in fieldnames(T)])...)
end

"""
    create(shape::T; resolution=nothing, rand_=0.0, type::Int64=1) where T <: Shape

Abstact function for creating **Shape** objects.

## Returns
    - X : Initial reference position
    - V : Initial velocity
    - Y : Initial position
    - volume : Volume per particle point
    - type: Type of particle point

"""
function create(shape::T; resolution=nothing, rand_=0.0, type::Int64=1) where T <: Shape
    error("Not implemented for type **$(typeof(shape))** yet.")
end

function mycat(a, b)
    if length(size(a))==1
        return vcat(a, b)
    else
        return hcat(a, b)
    end
end

function create(pobj::PostOpObj, args...; kwargs...)
    temp = nothing
    out = nothing
    for obj in pobj.objs
        temp = create(obj, args...; kwargs...)
        if isa(out, Nothing)
            out = temp
        else
            out = mycat.(out, temp)
        end
    end
    for func in pobj.operations
        out = func(out)
    end
    return out
end

function combine(obj1::T1, obj2::T2) where {T1<:Union{Shape, PostOpObj}, T2<:Union{Shape, PostOpObj}}
    PostOpObj([obj1, obj2], [])
end

function unpack(d::Dict)
    return d[:x], d[:v], d[:y], d[:volume], d[:type]
end

function repack(args...; keys_ = (:x, :v, :y, :volume, :type))
    d = Dict()
    for i in 1:5
        d[keys_[i]] = args[i]
    end
    d
end

function repack!(d::Dict, keys_, vals)
    if length(keys_)==length(vals)
        for i in 1:length(vals)
            d[keys_[i]] = vals[i]
        end
    else
        error("Lengths are not same.")
    end
    d
end

include("./cuboid.jl")
include("./sphere.jl")
include("./disk.jl")
include("./cylinder.jl")
include("./cone.jl")

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh