text:
using UCIData
using Test

using DataFrames

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

expected_sizes = Dict(
    "classification" => Dict(
        "acute-inflammations-1" => (120, 8, 5, 0),
        "acute-inflammations-2" => (120, 8, 5, 0),
        "arrhythmia" => (452, 281, 0, 408),
        "balance-scale" => (625, 6, 0, 0),
        "balloons-a" => (20, 6, 4, 0),
        "balloons-b" => (20, 6, 4, 0),
        "balloons-c" => (20, 6, 4, 0),
        "balloons-d" => (16, 6, 4, 0),
        "banknote-authentication" => (1372, 6, 0, 0),
        "blood-transfusion-service-center" => (748, 6, 0, 0),
        "breast-cancer" => (286, 11, 8, 9),
        "breast-cancer-wisconsin-diagnostic" => (569, 32, 0, 0),
        "breast-cancer-wisconsin-prognostic" => (198, 34, 0, 4),
        "breast-cancer-wisconsin-original" => (699, 11, 0, 16),
        "car-evaluation" => (1728, 8, 6, 0),
        "chess-king-rook-vs-king-pawn" => (3196, 38, 36, 0),
        "chess-king-rook-vs-king" => (28056, 8, 6, 0),
        "climate-model-simulation-crashes" => (540, 20, 0, 0),
        "cnae-9" => (1080, 858, 0, 0),
        "congressional-voting-records" => (435, 18, 16, 392),
        "connectionist-bench-sonar" => (208, 62, 0, 0),
        "connectionist-bench" => (990, 12, 0, 0),
        "contraceptive-method-choice" => (1473, 11, 1, 0),
        "credit-approval" => (690, 17, 9, 67),
        "cylinder-bands" => (540, 41, 19, 999),
        "dermatology" => (366, 36, 0, 8),
        "echocardiogram" => (132, 9, 0, 103),
        "ecoli" => (336, 9, 0, 0),
        "fertility" => (100, 11, 2, 0),
        "flags" => (194, 30, 6, 0),
        "glass-identification" => (214, 11, 0, 0),
        "haberman-survival" => (306, 5, 0, 0),
        "hayes-roth" => (132, 6, 0, 0),
        "heart-disease-cleveland" => (303, 15, 4, 6),
        "heart-disease-hungarian" => (294, 15, 4, 782),
        "heart-disease-switzerland" => (123, 15, 4, 273),
        "heart-disease-va" => (200, 15, 4, 698),
        "hepatitis" => (155, 21, 0, 167),
        "hill-valley-noise" => (606, 102, 0, 0),
        "hill-valley" => (606, 102, 0, 0),
        "horse-colic" => (300, 24, 9, 1604),
        "image-segmentation" => (210, 21, 0, 0),
        "indian-liver-patient" => (583, 12, 1, 4),
        "ionosphere" => (351, 36, 0, 0),
        "iris" => (150, 6, 0, 0),
        "lenses" => (24, 6, 4, 0),
        "letter-recognition" => (20000, 18, 0, 0),
        "libras-movement" => (360, 92, 0, 0),
        "lung-cancer" => (32, 58, 56, 5),
        "magic-gamma-telescope" => (19020, 12, 0, 0),
        "mammographic-mass" => (961, 7, 2, 162),
        "monks-problems-1" => (124, 8, 6, 0),
        "monks-problems-2" => (169, 8, 6, 0),
        "monks-problems-3" => (122, 8, 6, 0),
        "mushroom" => (8124, 24, 22, 2480),
        "nursery" => (12960, 10, 8, 0),
        "optical-recognition-handwritten-digits" => (3823, 66, 0, 0),
        "ozone-level-detection-eight" => (2534, 74, 0, 14937),
        "ozone-level-detection-one" => (2536, 74, 0, 14938),
        "parkinsons" => (195, 23, 0, 0),
        "pen-based-recognition-handwritten-digits" => (7494, 18, 0, 0),
        "planning-relax" => (182, 14, 0, 0),
        "poker-hand" => (25010, 12, 0, 0),
        "post-operative-patient" => (90, 10, 7, 3),
        "qsar-biodegradation" => (1055, 43, 0, 0),
        "seeds" => (210, 9, 0, 0),
        "seismic-bumps" => (2584, 20, 4, 0),
        "shuttle-landing-control" => (15, 8, 6, 26),
        "skin-segmentation" => (245057, 5, 0, 0),
        "soybean-large" => (307, 37, 35, 712),
        "soybean-small" => (47, 37, 35, 0),
        "spambase" => (4601, 59, 0, 0),
        "spect-heart" => (80, 24, 22, 0),
        "spectf-heart" => (80, 46, 0, 0),
        "statlog-project-german-credit" => (1000, 22, 13, 0),
        "statlog-project-landsat-satellite" => (4435, 38, 0, 0),
        "teaching-assistant-evaluation" => (151, 7, 4, 0),
        "thoracic-surgery" => (470, 18, 13, 0),
        "thyroid-disease-allbp" => (2800, 31, 22, 4556),
        "thyroid-disease-allhyper" => (2800, 31, 22, 4556),
        "thyroid-disease-allhypo" => (2800, 31, 22, 4556),
        "thyroid-disease-allrep" => (2800, 31, 22, 4556),
        "thyroid-disease-ann-thyroid" => (3772, 23, 0, 0),
        "thyroid-disease-dis" => (2800, 31, 22, 4556),
        "thyroid-disease-new-thyroid" => (215, 7, 0, 0),
        "thyroid-disease-sick-euthyroid" => (3163, 27, 18, 5329),
        "thyroid-disease-sick" => (2800, 31, 22, 4556),
        "thyroid-disease-thyroid-0387" => (9172, 31, 22, 14629),
        "tic-tac-toe-endgame" => (958, 11, 9, 0),
        "trains" => (10, 34, 22, 51),
        "wall-following-robot-navigation-2" => (5456, 4, 0, 0),
        "wall-following-robot-navigation-24" => (5456, 26, 0, 0),
        "wall-following-robot-navigation-4" => (5456, 6, 0, 0),
        "wine" => (178, 15, 0, 0),
        "yeast" => (1484, 10, 0, 0),
        "zoo" => (101, 18, 15, 0),
    ),
    "regression" => Dict(
        "abalone" => (4177, 10, 1, 0),
        "ailerons" => (13750, 42, 0, 0),
        "airfoil-self-noise" => (1503, 7, 0, 0),
        "airline-costs" => (31, 11, 0, 0),
        "auto-mpg" => (398, 9, 1, 6),
        "automobile" => (205, 27, 8, 59),
        "beer-aroma" => (23, 9, 0, 0),
        "blog-feedback" => (52397, 282, 200, 0),
        "communities-and-crime-2" => (2215, 126, 0, 41406),
        "communities-and-crime" => (1994, 124, 0, 36851),
        "computer-hardware" => (209, 10, 1, 0),
        "concrete-slump-test-compressive" => (103, 9, 0, 0),
        "concrete-slump-test-flow" => (103, 9, 0, 0),
        "concrete-slump-test-slump" => (103, 9, 0, 0),
        "construction-maintenance" => (33, 6, 0, 0),
        "cpu-act" => (8192, 23, 0, 0),
        "cpu-small" => (8192, 14, 0, 0),
        "elevators" => (16599, 20, 0, 0),
        "forest-fires" => (517, 14, 2, 0),
        "geographic-origin-music" => (1059, 70, 0, 0),
        "home-mortgage" => (18, 8, 0, 0),
        "housing" => (506, 15, 0, 0),
        "hybrid-price" => (153, 5, 0, 0),
        "immigrant-salaries" => (35, 5, 0, 0),
        "insurance-company-benchmark" => (5822, 87, 0, 0),
        "japan-emmigration" => (45, 7, 0, 0),
        "kin8nm" => (8192, 10, 0, 0),
        "lpga-2008" => (157, 8, 0, 0),
        "lpga-2009" => (146, 13, 0, 0),
        "online-news-popularity" => (39644, 61, 0, 0),
        "parkinsons-telemonitoring-motor" => (5875, 18, 0, 0),
        "parkinsons-telemonitoring-total" => (5875, 18, 0, 0),
        "pyrim" => (74, 29, 0, 0),
        "texas-jan-temp" => (16, 5, 0, 0),
        "triazines" => (186, 62, 0, 0),
        "tv-sales" => (31, 10, 0, 0),
        "vote-for-clinton" => (2704, 11, 0, 0),
        "wiki4he" => (913, 54, 8, 1995),
        "wine-quality-red" => (1599, 13, 0, 0),
        "wine-quality-white" => (4898, 13, 0, 0),
        "yacht-hydrodynamics" => (308, 8, 0, 0),
        "year-prediction-msd" => (515345, 92, 0, 0),
    ),
)

@testset "$datasettype" for datasettype in UCIData.list_dataset_types()
  @testset "$dataset" for dataset in UCIData.list_datasets(datasettype)
    df = UCIData.dataset(dataset)
    n_categoric = sum(map(x -> string(x)[1] == 'C', names(df)))
    n_missing = sum(count(ismissing, c) for c in eachcol(df))
    @test (size(df, 1), size(df, 2), n_categoric, n_missing) ==
          expected_sizes[datasettype][dataset]
  end
end

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
module HarmonicPowerModels

    # using pkgs
    using ProgressMeter

    # import pkgs
    import JuMP
    import PowerModels
    import SignalDecomposition
    import InfrastructureModels
    import Interpolations
    import Memento

    # import types
    import PowerModels: AbstractPowerModel, AbstractIVRModel
    import PowerModels: ids, ref, var, con, sol, nw_ids, nws
    import InfrastructureModels: replicate, sol_component_value_edge

    # pkg constants
    const _PMs = PowerModels
    const _HPM = HarmonicPowerModels
    const _SDC = SignalDecomposition
    const _IMs = InfrastructureModels
    const _INT = Interpolations

    function __init__()
        global _LOGGER = Memento.getlogger(PowerModels)
    end

    # const
    const freq = 50.0
    const nw_id_default = 1

    # funct
    sorted_nw_ids(pm) = sort(collect(nw_ids(pm)))

    # paths
    const BASE_DIR = dirname(@__DIR__)

    # include
    include("core/constraint_template.jl")
    include("core/data.jl")
    include("core/variable.jl")

    include("form/iv.jl")

    include("prob/hopf_iv.jl")
    include("prob/hpf_iv.jl")

    # export
    export BASE_DIR

    export replicate

    export run_hopf_iv, run_hpf_iv

end

ext: jl


text:
const UnitTriangle{T} = UnitSimplex{2,T,:closed}

ClassicalOrthogonalPolynomials.checkpoints(d::UnitTriangle{T}) where T = [SVector{2,T}(0.1,0.2), SVector{2,T}(0.2,0.3)]

# expansion in OPs orthogonal to
# x^a*y^b*(1-x-y)^c
# defined as
# P_{n-k}^{2k+b+c+1,a}(2x-1)*(1-x)^k*P_k^{c,b}(2y/(1-x)-1)

struct JacobiTriangle{T,V} <: BivariateOrthogonalPolynomial{T}
    a::V
    b::V
    c::V
end

JacobiTriangle(a::T, b::T, c::T) where T = JacobiTriangle{float(T),T}(a, b, c)
JacobiTriangle() = JacobiTriangle(0,0,0)
==(K1::JacobiTriangle, K2::JacobiTriangle) = K1.a == K2.a && K1.b == K2.b && K1.c == K2.c

axes(P::JacobiTriangle{T}) where T = (Inclusion(UnitTriangle{T}()),blockedrange(oneto(∞)))

copy(A::JacobiTriangle) = A

Base.summary(io::IO, P::JacobiTriangle) = print(io, "JacobiTriangle($(P.a), $(P.b), $(P.c))")

"""
    TriangleWeight(a, b, c)

is a quasi-vector representing `x^a * y^b * (1-x-y)^c`.
"""
struct TriangleWeight{T,V} <: Weight{T}
    a::V
    b::V
    c::V
end

TriangleWeight(a::T, b::T, c::T) where T = TriangleWeight{float(T),T}(a, b, c)

const WeightedTriangle{T} = WeightedBasis{T,<:TriangleWeight,<:JacobiTriangle}

WeightedTriangle(a, b, c) = TriangleWeight(a,b,c) .* JacobiTriangle(a,b,c)

axes(P::TriangleWeight{T}) where T = (Inclusion(UnitTriangle{T}()),)

Base.summary(io::IO, P::TriangleWeight) = print(io, "x^$(P.a)*y^$(P.b)*(1-x-y)^$(P.c) on the unit triangle")

@simplify function *(Dx::PartialDerivative{1}, P::JacobiTriangle)
    a,b,c = P.a,P.b,P.c
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockBroadcastArray(hcat,
        ((k .+ (b-1)) .* (n .+ k .+ (b+c-1)) ./ (2k .+ (b+c-1))),
        ((n .+ k .+ (a+b+c)) .* (k .+ (b+c)) ./ (2k .+ (b+c-1)))
        )
    JacobiTriangle(a+1,b,c+1) * _BandedBlockBandedMatrix(dat', axes(k,1), (-1,1), (0,1))
end

@simplify function *(Dy::PartialDerivative{2}, P::JacobiTriangle)
    a,b,c = P.a,P.b,P.c
    k = mortar(Base.OneTo.(oneto(∞)))
    T = promote_type(eltype(Dy), eltype(P)) # avoid bug in convert
    JacobiTriangle(a,b+1,c+1) * _BandedBlockBandedMatrix((k .+ convert(T, b+c))', axes(k,1), (-1,1), (-1,1))
end

# @simplify function *(Δ::Laplacian, P)
#     _lap_mul(P, eltype(axes(P,1)))
# end

@simplify function *(Ac::QuasiAdjoint{<:Any,<:JacobiTriangle}, B::JacobiTriangle)
    A = parent(Ac)
    @assert A == B == JacobiTriangle(0,0,0)
    a,b,c = A.a,A.b,A.c
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    _BandedBlockBandedMatrix((@. exp(loggamma(n+k+b+c)+loggamma(n-k+a+1)+loggamma(k+b)+loggamma(k+c)-loggamma(n+k+a+b+c)-loggamma(k+b+c)-loggamma(n-k+1)-loggamma(k))/((2n+a+b+c)*(2k+b+c-1)))',
                                axes(k,1), (0,0), (0,0))
end

function Wy(a,b,c)
    k = mortar(Base.OneTo.(oneto(∞)))
    _BandedBlockBandedMatrix((-k)', axes(k,1), (1,-1), (1,-1))
end
function Wx(a,b,c)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockVcat(
        ((k .+ (c-1)) .* ( k .- n .- 1 ) ./ (2k .+ (b+c-1)))',
        (k .* (k .- n .- a) ./ (2k .+ (b+c-1)))'
        )
    _BandedBlockBandedMatrix(dat, axes(k,1), (1,-1), (1,0))
end

 function Rx(a,b,c)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = BlockHcat(
        BroadcastVector((n,k,bc1,abc) -> (n + k +  bc1) / (2n + abc), n, k, b+c-1, a+b+c),
        BroadcastVector((n,k,abc) -> (n + k +  abc) / (2n + abc), n, k, a+b+c)
        )
    _BandedBlockBandedMatrix(dat', axes(k,1), (0,1), (0,0))
end
function Ry(a,b,c)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = PseudoBlockArray(Vcat(
        ((k .+ (c-1) ) .* (n .+ k .+ (b+c-1)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ((k .- n .- a ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ((k .+ (c-1) ) .* (k .- n .- 1) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ((n .+ k .+ (a+b+c) ) .* (k .+ (b+c)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))'
        ), (blockedrange(Fill(2,2)), axes(n,1)))
    _BandedBlockBandedMatrix(dat, axes(k,1), (0,1), (0,1))
end

function Lx(a,b,c)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = PseudoBlockArray(Vcat(
        ((n .- k .+ a) ./ (2n .+ (a+b+c)))',
        ((n .- k .+ 1) ./ (2n .+ (a+b+c)))'
        ), (blockedrange(Fill(1,2)), axes(n,1)))
    _BandedBlockBandedMatrix(dat, axes(k,1), (1,0), (0,0))
end
function Ly(a,b,c)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = PseudoBlockArray(Vcat(
        ((k .+ (b-1) ) .* (n .+ k .+ (b+c-1)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        (k .* (k .- n .- a ) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ((k .+ (b-1) ) .* (k .- n .- 1) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ( k .* (n .+ k .+ (a+b+c) ) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))'
        ), (blockedrange(Fill(2,2)), axes(n,1)))
    _BandedBlockBandedMatrix(dat, axes(k,1), (1,0), (1,0))
end
function Lz(a,b,c)
    n = mortar(Fill.(oneto(∞),oneto(∞)))
    k = mortar(Base.OneTo.(oneto(∞)))
    dat = PseudoBlockArray(Vcat(
        ((k .+ (c-1) ) .* (n .+ k .+ (b+c-1)) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        (k .* (n .- k .+ a ) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ((k .+ (c-1) ) .* (k .- n .- 1) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))',
        ( (-k) .* (n .+ k .+ (a+b+c) ) ./ ((2n .+ (a+b+c)) .* (2k .+ (b+c-1) )))'
        ), (blockedrange(Fill(2,2)), axes(n,1)))
    _BandedBlockBandedMatrix(dat, axes(k,1), (1,0), (1,0))
end

function \(w_A::WeightedTriangle, w_B::WeightedTriangle)
    wA,A = w_A.args
    wB,B = w_B.args

    @assert wA.a == A.a && wA.b == A.b && wA.c == A.c
    @assert wB.a == B.a && wB.b == B.b && wB.c == B.c

    if A.a + 1 == B.a && A.b == B.b && A.c == B.c
        Lx(B.a, B.b, B.c)
    elseif A.a == B.a && A.b + 1 == B.b && A.c == B.c
        Ly(B.a, B.b, B.c)
    elseif A.a == B.a && A.b == B.b && A.c + 1 == B.c
        Lz(B.a, B.b, B.c)
    else
        error("not implemented for $A and $wB")
    end
end

\(w_A::JacobiTriangle, w_B::WeightedTriangle) =
    (TriangleWeight(0,0,0) .* w_A) \ w_B

function \(A::JacobiTriangle, B::JacobiTriangle)
    if A.a == B.a && A.b == B.b && A.c == B.c
        Eye((axes(B,2),))
    elseif A.a == B.a + 1 && A.b == B.b && A.c == B.c
        Rx(B.a, B.b, B.c)
    elseif A.a == B.a && A.b == B.b + 1 && A.c == B.c
        Ry(B.a, B.b, B.c)
    # elseif A.a == B.a && A.b == B.b && A.c == B.c + 1
    #     Rz(B.a, B.b, B.c)
    else
        error("not implemented for $A and $B")
    end
end

## Jacobi matrix

jacobimatrix(::Val{1}, P::JacobiTriangle) = Lx(P.a+1, P.b, P.c) * Rx(P.a, P.b, P.c)
jacobimatrix(::Val{2}, P::JacobiTriangle) = Ly(P.a, P.b+1, P.c) * Ry(P.a, P.b, P.c)

## Evaluation

"""
Represents B⁺, the pseudo-inverse defined in forward recurrence.
"""
struct TriangleRecurrenceA{T} <: AbstractMatrix{T}
    bˣ::Vector{T}
    d::NTuple{3,T}  # last row entries
end

"""
Represents -B⁺*[Aˣ; Aʸ].
"""
struct TriangleRecurrenceB{T} <: AbstractMatrix{T}
    aˣ::Vector{T}
    bˣ::Vector{T}
    d::NTuple{2,T}
end

"""
Represents B⁺*[Cˣ; Cʸ].
"""
struct TriangleRecurrenceC{T} <: AbstractMatrix{T}
    bˣ::Vector{T}
    cˣ::Vector{T}
    d::T  # last row entries
end

function TriangleRecurrenceA(n, X, Y)
    bˣ = view(X,Block(n+1,n))[band(0)]
    Bʸ = view(Y,Block(n+1,n))'
    d₂ = inv(Bʸ[n,n+1])
    d₁ = -Bʸ[n,n] * d₂/bˣ[n]
    # max avoids bounds error, result unused
    d₀ = -Bʸ[n,max(1,n-1)] * d₂/bˣ[max(1,n-1)]
    TriangleRecurrenceA(bˣ, (d₀, d₁, d₂))
end

function TriangleRecurrenceB(n, X, Y)
    aˣ = view(X,Block(n,n))[band(0)]
    bˣ = view(X,Block(n+1,n))[band(0)]
    Aʸ = view(Y,Block(n,n))'
    Bʸ = view(Y,Block(n+1,n))'
    d₂ = inv(Bʸ[n,n+1])
    d₁ = Bʸ[n,n] * d₂/bˣ[n]*aˣ[n]-Aʸ[n,n]*d₂
    # max avoids bounds error, result unused
    d₀ = Bʸ[n,max(1,n-1)] * d₂ * aˣ[max(1,n-1)]/bˣ[max(1,n-1)] - Aʸ[n,max(1,n-1)]*d₂
    TriangleRecurrenceB(aˣ, bˣ, (d₀, d₁))
end

function TriangleRecurrenceC(n, X, Y)
    cˣ = view(X,Block(n-1,n))[band(0)]
    bˣ = view(X,Block(n+1,n))[band(0)]
    Cʸ = view(Y,Block(n-1,n))'
    Bʸ = view(Y,Block(n+1,n))'
    d₂ = inv(Bʸ[n,n+1])
    d = -Bʸ[n,n-1]/bˣ[n-1]*cˣ[n-1]*d₂ + Cʸ[n,n-1]*d₂
    TriangleRecurrenceC(bˣ, cˣ, d)
end

function size(A::TriangleRecurrenceA)
    n = length(A.bˣ)
    (n+1,2n)
end

function size(BA::TriangleRecurrenceB)
    n = length(BA.bˣ)
    (n+1, n)
end

function size(BA::TriangleRecurrenceC)
    n = length(BA.bˣ)
    (n+1, n-1)
end

function getindex(A::TriangleRecurrenceA{T}, k::Int, j::Int) where T
    n,m = size(A)
    if k < n && j == k
        return inv(A.bˣ[k])
    elseif k == n
        d₀, d₁, d₂ = A.d
        if j == m
            return d₂
        elseif j == n-1
            return d₁
        elseif j == n-2
            return d₀
        end
    end
    return zero(T)
end

function getindex(A::TriangleRecurrenceB{T}, k::Int, j::Int) where T
    n,m = size(A)
    if k < n && j == k
        return -A.aˣ[k] / A.bˣ[k]
    elseif k == n
        d₀, d₁ = A.d
        if j == m
            return d₁
        elseif j == m-1
            return d₀
        end
    end
    return zero(T)
end

function getindex(A::TriangleRecurrenceC{T}, k::Int, j::Int) where T
    n,m = size(A)
    if k ≤ m && j == k
        return A.cˣ[k] / A.bˣ[k]
    elseif k == n && j == m
        return A.d
    end
    return zero(T)
end

function xy_muladd!((x,y), A::TriangleRecurrenceA, v::AbstractVector, β, w::AbstractVector)
    n = size(A,1)
    d₀, d₁, d₂ = A.d
    w_1 = view(w, 1:n-1)
    w_1 .= A.bˣ .\ x .* v .+ β .* w_1
    if n > 2
        w[n] = d₀*x*v[n-2] + (d₁*x+d₂*y)*v[n-1] + β*w[n]
    else
        w[n] = (d₁*x+d₂*y)*v[n-1] + β*w[n]
    end
    w
end

function xy_muladd!(xy, Ac::Adjoint{<:Any,<:TriangleRecurrenceA}, v::AbstractVector, β, w::AbstractVector)
    x,y = xy
    A = parent(Ac)
    n = size(A,1)
    d₀, d₁, d₂ = A.d
    w .= A.bˣ .\ x .* view(v,1:n-1) .+ β .* w
    if n > 2
        w[n-2] += d₀*x*v[n]
        w[n-1] += (d₁*x+d₂*y)*v[n]
    else
        w[n-1] += (d₁*x+d₂*y)*v[n]
    end
    w
end

function ArrayLayouts.muladd!(α::T, BA::TriangleRecurrenceB{T}, c::AbstractVector{T}, β::T, dest::AbstractVector{T}) where T
    @boundscheck (size(BA,2) == length(c) && size(BA,1) == length(dest)) || throw(DimensionMismatch())
    aˣ,bˣ = BA.aˣ,BA.bˣ
    Ba_1, Ba_2 = BA.d

    n = length(bˣ)
    w = view(dest,1:n)
    w .= (-α) .* (bˣ .\ aˣ .* c) .+ β .* w
    dest[n+1] = α*(n > 1 ? Ba_1 * c[end-1] + Ba_2 * c[end] : Ba_2 * c[end]) + β*dest[n+1]
    dest
end

function ArrayLayouts.muladd!(α::T, BAc::Adjoint{T,<:TriangleRecurrenceB{T}}, c::AbstractVector{T}, β::T, dest::AbstractVector{T}) where T
    @boundscheck (size(BAc,2) == length(c) && size(BAc,1) == length(dest)) || throw(DimensionMismatch())
    BA = parent(BAc)
    aˣ,bˣ = BA.aˣ,BA.bˣ
    Ba_1, Ba_2 = BA.d

    n = length(bˣ)
    dest .= (-α) .* (bˣ .\ aˣ .* view(c,1:n)) .+ β .* dest
    if n > 1
        dest[end-1] += α*Ba_1*c[n+1]
    end
    dest[end] += α*Ba_2*c[n+1]
    dest
end

function ArrayLayouts.muladd!(α::T, BA::TriangleRecurrenceC{T}, u::AbstractVector{T}, β::T, dest::AbstractVector{T}) where T
    @boundscheck (size(BA,2) == length(u) && size(BA,1) == length(dest)) || throw(DimensionMismatch())
    bˣ,cˣ = BA.bˣ,BA.cˣ
    d = BA.d

    n = length(bˣ)
    w = view(dest,1:n-1)
    w .= α .* (view(bˣ,1:n-1) .\ cˣ .* u) .+ β .* w
    dest[n] *= β
    dest[n+1] = α*d * u[end] + β*dest[n+1]
    dest
end

function ArrayLayouts.muladd!(α::T, BAc::Adjoint{T,TriangleRecurrenceC{T}}, u::AbstractVector{T}, β::T, w::AbstractVector{T}) where T
    @boundscheck (size(BAc,2) == length(u) && size(BAc,1) == length(w)) || throw(DimensionMismatch())
    BA = parent(BAc)
    bˣ,cˣ = BA.bˣ,BA.cˣ
    d = BA.d

    n = length(bˣ)
    w .= α .* (view(bˣ,1:n-1) .\ cˣ .* view(u,1:n-1)) .+ β .* w
    w[n-1] += α*d*u[n+1]

    w
end

# following also resizes
function LinearAlgebra.lmul!(BAc::Adjoint{T,TriangleRecurrenceC{T}}, u::AbstractVector{T}) where T
    @boundscheck size(BAc,2) == length(u) || throw(DimensionMismatch())
    BA = parent(BAc)
    bˣ,cˣ = BA.bˣ,BA.cˣ
    d = BA.d

    n = length(bˣ)
    v = view(u,1:n-1)
    v .= view(bˣ,1:n-1) .\ cˣ .* v
    u[n-1] += d*u[n+1]

    resize!(u, n-1)
end

function LinearAlgebra.mul!(dest::AbstractVector{T}, BA::Union{TriangleRecurrenceB,TriangleRecurrenceC,Adjoint{<:Any,<:TriangleRecurrenceB},Adjoint{<:Any,<:TriangleRecurrenceC}}, c::AbstractVector) where T
    ArrayLayouts.zero!(dest) # avoid NaN
    muladd!(one(T), BA, c, zero(T), dest)
end

##
# TODO: almost generic, but left specialised for now because we need
# to make block indexing of infinite Jacobi operators fast.
# This has to do with recognizing a block-row vector with banded blocks
# times a block-col vec with banded blocks is banded.
# The functionality will need to be added to LazyBandedMatrices.jl
function tri_forwardrecurrence(N::Int, X, Y, x, y)
    T = promote_type(eltype(X),eltype(Y))
    ret = PseudoBlockVector{T}(undef, 1:N)
    N < 1 && return ret
    ret[1] = 1
    N < 2 && return ret
    ret[2] = x/X[1,1]-1
    ret[3] = -Y[2,1]/(Y[3,1]*X[1,1]) * (x-X[1,1]) + (y-Y[1,1])/Y[3,1]

    X_N = X[Block.(1:N), Block.(1:N-1)]
    Y_N = Y[Block.(1:N), Block.(1:N-1)]

    for n = 2:N-1
        # P[n+1,xy] == (A[n]*[x*Eye(n); y*Eye(n)] + B[n])*P[n,xy] - C[n]*P[n-1,xy]
        A = TriangleRecurrenceA(n, X_N, Y_N)
        B = TriangleRecurrenceB(n, X_N, Y_N)
        C = TriangleRecurrenceC(n, X_N, Y_N)

        p_0 = view(ret,Block(n-1))
        p_1 = view(ret,Block(n))
        p_2 = view(ret,Block(n+1))

        mul!(p_2, B, p_1)
        xy_muladd!((x,y), A, p_1, one(T), p_2)
        muladd!(-one(T), C, p_0, one(T), p_2)
    end
    ret
end

getindex(P::JacobiTriangle, xy::SVector{2}, JR::BlockOneTo) =
    tri_forwardrecurrence(Int(JR[end]), jacobimatrix(Val(1), P), jacobimatrix(Val(2), P), xy...)

###
# Clenshaw
###

function getindex(f::Expansion{T,<:JacobiTriangle}, xy::SVector{2}) where T
    P,c∞ = arguments(f)
    c = paddeddata(c∞)
    N = blocksize(c,1)

    N == 1 && return c[1]
    X = jacobimatrix(Val(1), P)
    Y = jacobimatrix(Val(2), P)
    X_N = X[Block.(1:N), Block.(1:N-1)]
    Y_N = Y[Block.(1:N), Block.(1:N-1)]
    γ1 = Array{T}(undef, N-1)
    γ2 = Array{T}(undef, N)
    A = TriangleRecurrenceA(N-1, X_N, Y_N)
    B = TriangleRecurrenceB(N-1, X_N, Y_N)
    copyto!(γ2, view(c,Block(N)))
    copyto!(γ1, view(c,Block(N-1)))
    muladd!(one(T), B', γ2, one(T), γ1)
    xy_muladd!(xy, A', γ2, one(T), γ1)

    for n = N-2:-1:1
        # P[n+1,xy] == (A[n]*[x*Eye(n); y*Eye(n)] + B[n])*P[n,xy] - C[n]*P[n-1,xy]
        A,B,C = TriangleRecurrenceA(n, X_N, Y_N),TriangleRecurrenceB(n, X_N, Y_N),TriangleRecurrenceC(n+1, X_N, Y_N)
        # some magic! C can be done in-place, otherwise
        # we would need a secondary buffer.
        # this also resizes γ2
        lmul!(C', γ2)
        # Need to swap sign since it should have been -C
        γ2 .= (-).(γ2) .+ view(c,Block(n))
        γ2,γ1 = γ1,γ2
        muladd!(one(T), B', γ2, one(T), γ1)
        xy_muladd!(xy, A', γ2, one(T), γ1)
    end

    γ1[1]
end

getindex(f::Expansion{T,<:JacobiTriangle}, xys::AbstractArray{<:SVector{2}}) where T =
    [f[xy] for xy in xys]

# getindex(f::Expansion{T,<:JacobiTriangle}, x::AbstractVector{<:Number}) where T =
#     copyto!(Vector{T}(undef, length(x)), view(f, x))

# function copyto!(dest::AbstractVector{T}, v::SubArray{<:Any,1,<:Expansion{<:Any,<:JacobiTriangle}, <:Tuple{AbstractVector{<:Number}}}) where T
#     f = parent(v)
#     (x,) = parentindices(v)
#     P,c = arguments(f)
#     clenshaw!(paddeddata(c), recurrencecoefficients(P)..., x, Fill(_p0(P), length(x)), dest)
# end

lgamma(x) = logabsgamma(x)[1]

# FastTransforms uses normalized Jacobi polynomials so this corrects the normalization
function _ft_trinorm(n,k,a,b,c)
    if a == 0 && b == c == -0.5
        k == 0 ? sqrt((2n+1)) /sqrt(2π) : k*sqrt(2n+1)*exp(lgamma(k)-lgamma(k+0.5))
    else
        sqrt((2n+b+c+a+2)*(2k+c+b+1))*exp(k*log(2)+(lgamma(n+k+b+c+a+2)+lgamma(n-k+1)-log(2)*(2k+a+b+c+2)-
                lgamma(n+k+b+c+2)-lgamma(n-k+a+1) + lgamma(k+c+b+1)+lgamma(k+1)-
                    log(2)*(c+b+1)-lgamma(k+c+1)-lgamma(k+b+1))/2)
    end
end

function tridenormalize!(F̌,a,b,c)
    for n = 0:size(F̌,1)-1, k = 0:n
        F̌[n-k+1,k+1] *= _ft_trinorm(n,k,a,b,c)
    end
    F̌
end

function _trigrid(N::Integer)
    M = N
    x = [sinpi((2N-2n-1)/(4N))^2 for n in 0:N-1]
    w = [sinpi((2M-2m-1)/(4M))^2 for m in 0:M-1]
    [SVector(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
end
trigrid(N::Block{1}) = _trigrid(Int(N))

function grid(Pn::SubQuasiArray{T,2,<:JacobiTriangle,<:Tuple{Inclusion,BlockSlice{<:BlockRange{1}}}}) where T
    kr,jr = parentindices(Pn)
    trigrid(maximum(jr.block))
end

struct TriTransform{T}
    tri2cheb::FastTransforms.FTPlan{T,2,FastTransforms.TRIANGLE}
    grid2cheb::FastTransforms.FTPlan{T,2,FastTransforms.TRIANGLEANALYSIS}
    a::T
    b::T
    c::T
end

TriTransform{T}(F::AbstractMatrix{T}, a, b, c) where T =
    TriTransform{T}(plan_tri2cheb(F, a, b, c), plan_tri_analysis(F), a, b, c)

TriTransform{T}(N::Block{1}, a, b, c) where T = TriTransform{T}(Matrix{T}(undef, Int(N), Int(N)), a, b, c)

*(T::TriTransform, F::AbstractMatrix) = DiagTrav(tridenormalize!(T.tri2cheb\(T.grid2cheb*F),T.a,T.b,T.c))

function factorize(V::SubQuasiArray{T,2,<:JacobiTriangle,<:Tuple{Inclusion,BlockSlice{BlockOneTo}}}) where T
    P = parent(V)
    _,jr = parentindices(V)
    N = maximum(jr.block)
    TransformFactorization(grid(V), TriTransform{T}(N, P.a, P.b, P.c))
end


text:
"""
Concatenate strings.
"""
function cat_str end

"""
    cat_str(a, b)

Concatenate strings.
"""
function cat_str(a, b)
 string(a, b)
end

"""
    cat_str(a, b, c)

Concatenate strings.
"""
function cat_str(a, b, c)
  string(a, b, c)
end


text:
"""
  `SingerAlg.PrimesGAP.jl` provides a few functions based on GAP
  that are faster than the corresponding functions from the Julia package
  `Primes.jl`.
  Since one of the aims of GAP's `SingerAlg` package is
  to compare implementations in GAP with Julia implementations,
  we use GAP's factorization code in both situations.
"""
module PrimesGAP

import GAP

isprime(n::T) where {T<:Integer} = GAP.Globals.IsPrimeInt(GAP.julia_to_gap(n))

nextprime(n::T) where {T<:Integer} = GAP.Globals.NextPrimeInt(GAP.julia_to_gap(n-1))

factors(n::T) where {T<:Integer} = Vector{Int}(GAP.Globals.Factors(GAP.julia_to_gap(n)))
#T We would better provide a function `factor` as in `Primes.jl`,
#T to whose result one can apply `sort(collect(Set, .))`.

end

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
# module Bukdu

function plug(func::Function, arg::Any)
    func(arg)
end

module Plug

import ..plug

include("plug/router.jl")
include("plug/static.jl")
include("plug/logger.jl")

end # module Bukdu.Plug


text:
insts = OrderedDict()
insts["minimal"] = [
    # rootdet
    ((3, MatNegRtdet()),),
    ((3, MatNegRtdetEFExp()),),
    ((3, MatNegRtdetEFPow()),),
    # tr neglog
    ((2, MatLogdetCone()),),
    ((2, MatNegLog()),),
    ((2, MatNegLogDirect()),),
    # tr negentropy
    ((3, MatNegEntropy()),),
    ((3, MatNegEntropyEigOrd()),),
    # tr negsqrt
    ((3, MatNegSqrt()),),
    ((3, MatNegSqrtEigOrd()),),
    # tr negpower01
    ((3, MatNegPower01(0.7)),),
    ((3, MatNegPower01EigOrd(0.7)),),
    # tr power12
    ((3, MatPower12(1.3)),),
    ((3, MatPower12EigOrd(1.3)),),
    ]
insts["fast"] = [
    ((50, MatNegRtdet()),),
    ((12, MatNegRtdetEFExp()),),
    ((8, MatNegRtdetEFPow()),),
    ((40, MatNegSqrt()),),
    ((5, MatNegSqrtEigOrd()),),
    ((6, MatNegSqrtConjEigOrd()),),
    ((25, MatNegSqrtConjDirect()),),
    ((100, MatLogdetCone()),),
    ((75, MatNegLog()),),
    ((7, MatNegLogEigOrd()),),
    ((20, MatNegLogDirect()),),
    ((60, MatNegEntropy()),),
    ((8, MatNegEntropyEigOrd()),),
    ((50, MatPower12(1.5)),),
    ((5, MatPower12EigOrd(1.5)),),
    ((50, MatNegPower01(0.4)),),
    ((40, MatNegPower01Conj(0.4)),),
    ((5, MatNegPower01EigOrd(0.4)),),
    ]
insts["various"] = [
    ((400, MatNegRtdet()),),
    ((70, MatNegRtdetEFExp()),),
    ((40, MatNegRtdetEFPow()), nothing, (default_tol_relax = 100,)),
    ((200, MatNegSqrt()),),
    ((10, MatNegSqrtConjEigOrd()),),
    ((125, MatNegSqrtConjDirect()),),
    ((300, MatLogdetCone()),),
    ((150, MatNegLog()),),
    ((12, MatNegLogEigOrd()),),
    ((75, MatNegLogDirect()),),
    ((150, MatNegEntropy()),),
    ((10, MatNegEntropyEigOrd()),),
    ((80, MatNegPower01(0.4)),),
    ((40, MatNegPower01Conj(0.6)),),
    ((15, MatPower12(1.5)),),
    ((9, MatPower12EigOrd(1.5)), nothing, (default_tol_relax = 100,)),
    ]
return (ExperimentDesignJuMP, insts)


text:
## manipulate to go here
using Mustache

## make_control
##
## maps a tuple to a control
## The tuple has the form:
## (:nm, values..., {:label=>"label", :initial=>"initial", :control=>"slider"})
## where the ending `Dict` is optional
##
## By default, the initial symbol, `:nm`, represents an unbound variable in
## the manipulated expression. During the evaluation of the
## expression, its value comes from the control.
##
## The `values...` is used to specify the range of acceptable values or type of value. The

## (:nm, Bool) -> checkbox
make_control(parent::Container, x::Bool) = checkbox(parent, x)

## (:nm, Range) -> slider
function make_control(parent::Container, x::Union(Range, Range1, Ranges))
    obj = slider(parent, x)
    obj[:size] = [100, 20]
    obj
end

## (:nm, Range, Range) -> 2dslider
make_control(parent::Container, x::Union(Range, Range1,Ranges), y::Union(Range, Range1,Ranges)) = slider2d(parent, x, y)

## (:nm, Vector) -> length <= 3 -> radio else combo. (Can be buttongroup, slider, listbox, ...
function make_control(parent::Container, x::Vector)
    if length(x) <= 3
        radiogroup(parent, x, orientation=:vertical)
    else
        cb = combobox(parent, x)
        cb[:value] = cb.model.items[1] # initialize
        cb
    end
end

## (:nm, String) -> edit
make_control(parent::Container, x::String) = lineedit(parent, x)

## (:nm, Real) -> edit with coercion
make_control(parent::Container, x::Real) = lineedit(parent, x, coerce=parsefloat)
## (:nm, Int) -> edit with coercion
make_control(parent::Container, x::Int) = lineedit(parent, x, coerce=parseint)

## evaluation context
## context to store dynamic values
module ManipulateContext
end

function dict_to_module(d::Dict, context) ## stuff values into Manipulate Context
    for (k,v) in d
        eval(context, :($(symbol(k)) = $v))
    end
end

## The main manipulate object and some methods
type ManipulateObject
    expr
    controls
    window
    output_area
    toolkit
end

manipulate_signal(x) = "valueChanged"
manipulate_signal(x::LineEdit) = "editingFinished"

setsize(object::ManipulateObject, sz::Vector) = setsize(object.toolkit, object.window, sz)
function getValue(object::ManipulateObject)
    d = Dict()
    for (k, v) in object.controls
        if !isa(v, Button)
            d[k] = getValue(v)
        end
    end
    d
end

## set GUI value for symbol to value
function setValue(object::ManipulateObject, key::Symbol, value)
    ## get id from controls
    ctrl = filter(u -> isa(u, MControlType) && u.nm == key, object.controls)
    if length(ctrl) > 0
        setValue(object.toolkit, ctrl[1], value)
    end
end
## Give [] interface to values
import Base: getindex, setindex!
function getindex(object::ManipulateObject, x::Symbol)
    ctrl = filter(u -> isa(u, MControlType) && u.nm == x, object.controls)
    length(ctrl) > 0 ? getValue(object.toolkit, ctrl[1]) : nothing
end
setindex!(object::ManipulateObject, value,  x::Symbol) = setValue(object, x, value)

ClearDisplay(self::ManipulateObject) = ClearDisplay(self.toolkit, self)

## Manipulate an expression with controls specified by args...
##
## `expr`: an expression to manipulate. Last value computed is
## displayed.  Variables in the expression can refer to symbols
## specified through the args.... These are replaced with the values
## generated by the GUI. Expressions are easily created with `quote`
## blocks, though there are some subtle issues, such as loading of
## modules.
##
## Expressions which return `Winston` objects are plotted, others are
## displayed as text objects after coercion by the `string` method.
##
## `args...`: A collection of control specifications. These can be
## simple specification, such as (:n, 1:10) which would map to a
## slider for the symbol n. There are other familiar control types
## without a shortcut, such as `MInput` for text input.
##
## Keyword arguments:
##
## `toolkit`: MIME specification of toolkit,
## eg. `MIME("application/x-tcltk")`
##
## `title`: title for window
##
## `control_placement`: a symbol, `:left` or `top` to indicate where
## output window should go.
##
## `modules`: A vector of symbols listing needed modules to load. One
## can't call `using ModuleName` within the expression to manipulate.
##
## Output:
##
## Creates a window with controls to manipulate the expression. The
## return value is of type `ManipulateObject`.
function manipulate(expr, args...;
                    toolkit::MIME=default_toolkit, ##MIME("application/x-tcltk"),
                    title::String="Manipulate",
                    control_placement::Symbol=:left, ## or :left
                    modules::Vector=[],              # eg. [:Winston, :SymPy]
                    width::Int=600,
                    height::Int=400,
                    kwargs...)

    controls = Any[]
    self = ManipulateObject(expr, controls, nothing, nothing, toolkit)

    ## should provide per instance context
    context = ManipulateContext

    for i in modules
        f = "manipulate-$(string(i)).jl"
        if isfile(f)
            require(f)
        end
        eval(context, Expr(:using, i))
    end

    eval(context, :($(:self) = $self))

    function update_gui(value)
        ## need to push these into context
        d = getValue(self)
        dict_to_module(d, context)
        p = eval(context, expr)
        if !isa(p, Nothing)
            Display(toolkit, self, p, context=context)
        end
    end

    set_these = Dict()

    controls = Dict()

    ## layout
    self.window = w = window(toolkit=toolkit, title=title, size=[width, height], visible=false)
    istk() && Tk.pack_stop_propagate(w[:widget])

    f = hbox(w)
    push!(w, f)

    lb = vbox(f)
    lb[:sizepolicy] = (:fixed, :fixed)
    lb[:alignment] = (:left, :top)

    self.output_area = rb = vbox(f)
    rb[:sizepolicy] = (:expand, :expand)
    isqt() && push!(rb, pyplotgraphic(rb))

    push!(f, lb)
    push!(f, rb)

    fl = formlayout(lb)
    fl[:sizepolicy] = (:expand, :fixed)
    for ctrl in args
        if isa(ctrl, Tuple)
            nm = ctrl[1]
            isa(nm, Symbol) || error("Initial value of a control is a symbol")
            ## do we have a dict? does it have a control?
            d = ctrl[end]
            if isa(ctrl[end], Dict) && haskey(d, :control)
                widget = d[:control](fl, ctrl[2:end-1]...)
            elseif isa(ctrl[end], Dict) && !haskey(d, :control)
                widget = make_control(fl, ctrl[2:end-1]...)
            else
                widget = make_control(fl, ctrl[2:end]...)
            end

            ## initial
            if isa(d, Dict) && haskey(d, :initial)
                set_these[nm] = d[:initial]
            end
            ## label {:label=>"string"} {:label=>nothing} or skip to get string(:nm)
            label = (isa(d, Dict) && haskey(d, :label)) ? d[:label] : string(ctrl[1])

            connect(widget, manipulate_signal(widget), update_gui)
            push!(fl, widget, label)
            controls[nm] = widget
        elseif isa(ctrl, String)
            if ctrl == "----"
                push!(fl, separator(fl), nothing)
            else
                push!(fl, label(fl, ctrl), nothing)
            end
        else

            error("What is this doing here?")
        end
    end

    push!(lb, fl)

    self.controls = controls
    [setValue(self, k, v) for (k, v) in set_these] # initialize any values
    raise(w)
    update_gui(nothing)

    self.window[:visible] = true
    self
end


text:
module Style
using Printf
using Markdown
using Cairo
using ComplexVisual
@ComplexVisual.import_huge

import Main.DocGenerator: DocSource, DocCreationEnvironment, DocContext,
        Document, substitute_marker_in_markdown, create_doc_icon, append_md

"""
# [![./Style_docion.png]({image_from_canvas: get_doc_icon()}) Style](./Style.md)

Styles (i.e. subtypes of `CV_ContextStyle`) are used to govern the appearance
of painting and drawing actions. They are used to "bundle" typical
(Cairo-context) changes in order to make them reusable and easily combinable
(with `→:Tuple{T, S} where {T<:CV_ContextStyle, S<:CV_ContextStyle}`).

## Quick links

* `cv_color` (`cv_black`  `cv_white`)
* `cv_linewidth`
* `cv_antialias` (`cv_antialias_best`  `cv_antialias_none`)
* `cv_opmode` (`cv_op_source`  `cv_op_over`)
* `cv_fillstyle` (`cv_fill_winding`  `cv_fill_even_odd`)
* `cv_fontface`  `cv_fontsize`
* `CV_CombiContextStyle`   `→:Tuple{T, S} where {T<:CV_ContextStyle, S<:CV_ContextStyle}`
* `cv_create_context:Tuple{Function, CV_Canvas, Union{Nothing, CV_ContextStyle}}`

## How styles work

Before the painting and/or drawing operation(s) the function `cv_prepare` for
subtypes of `CV_ContextStyle` is called:

```
cv_prepare(context::C, style::S)
    context     C <: CV_Context          (often C <: CV_CanvasContext)
    style       S <: CV_ContextStyle     (often S <: CV_CanvasContextStyle)
```

Several styles can be combined
(with `→:Tuple{T, S} where {T<:CV_ContextStyle, S<:CV_ContextStyle}`) to
a single style:

```
new_style = cv_black → cv_linewidth(3) → cv_antialias_best
```

Styles can be attached to painters with `↦:Tuple{CV_ContextStyle, CV_Painter}`:

```
cv_color(0.7, 0.4, 0.4) ↦ CV_FillPainter()
```

Use parenthesis to combine more styles and attach them to a painter
(parens are needed because arrows are right-associative):

```
( cv_black → cv_linewidth(2) )  ↦ CV_LinePainter(...)
```

## `doc: CV_ContextStyle`

## `doc: CV_CanvasContextStyle`

"""
style_intro() = nothing

"""
## `doc: cv_create_context:Tuple{Function, CV_Canvas, Union{Nothing, CV_ContextStyle}}`
"""
help_create_context() = nothing

"""
## `doc: CV_CombiContextStyle`

## `doc: →:Tuple{T, S} where {T<:CV_ContextStyle, S<:CV_ContextStyle}`
"""
help_combi_style() = nothing

"""
## `doc: cv_color`

## `doc: cv_black`

## `doc: cv_white`
"""
help_color() = nothing

"""
## `doc: cv_linewidth`
"""
help_linewidth() = nothing

"""
## `doc: cv_antialias`

## `doc: cv_antialias_best`

## `doc: cv_antialias_none`
"""
help_antialias() = nothing

"""
## `doc: cv_operatormode`

## `doc: cv_opmode`

## `doc: cv_op_source`

## `doc: cv_op_over`
"""
help_operator() = nothing

"""
## `doc: cv_fillstyle`

## `doc: cv_fill_winding`

## `doc: cv_fill_even_odd`
"""
help_fillstyle() = nothing

"""
## `doc: cv_fontface`

## `doc: cv_fontsize`
"""
help_font() = nothing

"""
## `doc: CV_MathCoorStyle`
"""
help_coor_style() = nothing

function get_doc_icon()
    layout = CV_2DLayout()

    math_can = CV_Math2DCanvas(0.0 + 1.0im, 1.0 + 0.0im, 70)

    bg_painter = cv_white ↦ CV_FillPainter()
    line_painters = Vector{CV_StyledPainter}()
    for x in -0.5:0.2:0.9
        line_segs = [ [x+0.0im, x+0.5+1.0im ] ]
        style = cv_color(0.3, 0.3, min(0.2+0.5+x, 1.0)) → cv_linewidth(round( 5*(x+1.0) ))
        push!(line_painters, style ↦ CV_LinePainter(line_segs))
    end

    cv_create_context(math_can) do canvas_context
        cv_paint(canvas_context, bg_painter)
        for lp in line_painters
            cv_paint(canvas_context, lp)
        end
    end

    can_l = cv_add_canvas!(layout, math_can, cv_anchor(math_can, :center), (0, 0))

    stext = cv_text("Style", cv_black → cv_fontface("serif") → cv_fontsize(25))
    ptext = cv_add_canvas!(layout, stext, cv_anchor(stext, :center), (0, 0))

    can_layout = cv_canvas_for_layout(layout)
    cv_create_context(can_layout) do con_layout
        can_l(con_layout)
        ptext(con_layout)
    end

    icon = create_doc_icon(can_layout)
    return icon
end

function create_document(doc_env::DocCreationEnvironment)
    doc_source = DocSource("Style", @__MODULE__)
    context = DocContext(doc_env, doc_source)

    md = Markdown.MD()
    for part in (style_intro, help_create_context, help_combi_style,
            help_color, help_linewidth, help_antialias, help_operator,
            help_fillstyle, help_font, help_coor_style)
        part_md = Base.Docs.doc(part)
        substitute_marker_in_markdown(context, part_md)

        append_md(md, part_md)
    end

    doc = Document(doc_source, md)
    return doc
end

end
# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using NgSpice
using Test

rootdir = joinpath(@__DIR__, "..")
netpath = joinpath(rootdir, "inputs", "mosfet.cir")

@testset "interface" begin
    p = convert(Ptr{Nothing}, 0)
    @test ngSpice_Init(p, p, p, p, p, p, p) == 0
    @test ngSpice_Command("source $netpath") == 0
    @test ngSpice_CurPlot() != C_NULL
    @test ngSpice_Command("run") == 0
    @test typeof(ngSpice_AllPlots()) == Ptr{Ptr{UInt8}}
    @test ngSpice_running() == 0
end

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
@testset "Tables" begin
    X = (
        x = rand(5),
        y = rand(Int, 5),
        z = categorical(collect("asdfa")),
        w = rand(5)
    )
    s = schema(X)
    @test info(X) == schema(X)
    @test s.scitypes == (Continuous, Count, Multiclass{4}, Continuous)
    @test s.types == (Float64, Int64, CategoricalValue{Char,UInt32}, Float64)
    @test s.nrows == 5

    @test_throws ArgumentError schema([:x, :y])

    t = scitype(X)
    @test t <: Table(Continuous, Finite, Count)
    @test t <: Table(Infinite, Multiclass)
    @test !(t <: Table(Continuous, Union{Missing, Count}))

    @test MLJScientificTypes._nrows(X) == 5
    @test MLJScientificTypes._nrows(()) == 0
    @test MLJScientificTypes._nrows((i for i in 1:7)) == 7
end

@testset "csvfile" begin
    X = (x = rand(4), )
    CSV.write("test.csv", X)
    file = CSV.File("test.csv")
    @test scitype(file) == scitype(X)
    rm("test.csv")
end


text:
using LightGraphs
using MetaGraphs
using Test

using LabelledGraphs

@testset "LabelledEdges compare equall iff their source and destination are equal" begin
    @test LabelledEdge((1, 1), (2, 2)) == LabelledEdge((1, 1), (2, 2))
    @test LabelledEdge((2, 2), (1, 1)) ≠ LabelledEdge((1, 1), (2, 2))
    @test LabelledEdge((1, 1), (2, 2)) ≠ LabelledEdge((1, 1), (3, 3))
end

for graph_type ∈ [SimpleGraph, SimpleDiGraph]
@testset "Initializing LabelledGraph with another graph of type $graph_type" begin

    @testset "fails if number of labels is different than source graph's number of vertices" begin
        g = graph_type(5)
        labels = [1, 4, 5, 3]
        @test_throws ArgumentError LabelledGraph(labels, g)
    end

    @testset "fails if labels are not unique" begin
        g = path_digraph(4)
        labels = [1, 5, 10, 5]
        @test_throws ArgumentError LabelledGraph(labels, g)
    end

    @testset "gives a graph isomorphic to the source graph" begin
        g = graph_type(4)
        for (i, j) ∈ [(1, 2), (3, 2), (2, 4)] add_edge!(g, i, j) end
        labels = [20, 4, 5, 6]

        lg = LabelledGraph(labels, g)

        @testset "number of vertices in both graphs is the same" begin
            @test nv(lg) == nv(g)
        end

        @testset "sequence of vertices of LabelledGraph is equal to the labels used" begin
            @test collect(vertices(lg)) == labels
        end

        @testset "presence and absence of vertices is correctly reported" begin
            for v ∈ (20, 4, 5, 6)
                @test has_vertex(lg, v)
            end

            @test !has_vertex(lg, 1)
        end

        @testset "number of edges in both graph sis the same" begin
            @test ne(lg) == ne(g) == length(edges(lg))
        end

        # For non-directed LightGraphs, edges are reported lex-ordered, regardless of
        # the ordering of the original edge (i.e. edge (2, 1) is always reported as (1, 2)).
        # Since this ordering propagates to our LabelledGraphs, we need to define set of
        # expected edges accordingly.
        expected_edges = is_directed(g) ?
            [(i, j) for (i, j) ∈ [(20, 4), (5, 4), (4, 6)]] :
            [(i, j) for (i, j) ∈ [(20, 4), (4, 5), (4, 6)]]
        @testset "set of edges of LabelledGraph comprises source graph's edges with vertices mapped to labels" begin
            @test Set([(src(e), dst(e)) for e ∈ edges(lg)]) == Set(expected_edges)
        end

        @testset "presence and absence of edges is correctly reported" begin
            for (u, v) ∈ expected_edges
                @test has_edge(lg, LabelledEdge(u, v))
                @test has_edge(lg, u, v)
            end

            @test !has_edge(lg, 1, 2)
            @test !has_edge(lg, LabelledEdge(1, 2))
        end

        @testset "LabelledGraph is directed iff source graph is also directed" begin
            @test is_directed(lg) == is_directed(g)
        end
    end
end
end

for graph_type ∈ [SimpleGraph, SimpleDiGraph]
@testset "Initializing LabelledGraph{$graph_type} with only labels" begin
    lg = LabelledGraph{graph_type}([1, 7, 3, 5])

    @testset "gives graph with vertices corresponding to the labels" begin
        @test vertices(lg) == [1, 7, 3, 5]
    end

    @testset "and an empty set of edges" begin
        @test ne(lg) == length(edges(lg)) == 0
    end
end
end

@testset "Adding edge (s, d) to undirected LabelledGraph" begin
    lg = LabelledGraph{SimpleGraph}([4, 5, 0])

    add_edge!(lg, 4, 5)
    add_edge!(lg, LabelledEdge(4, 0))

    @testset "makes both (s, d) and (d, s) present in the graph begin" begin
        @test has_edge(lg, 4, 5) && has_edge(lg, 5, 4)
        @test has_edge(lg, 4, 0) && has_edge(lg, 0, 4)
    end

    @testset "does not add any other edge" begin
        @test ne(lg) == length(edges(lg)) == 2
    end

    @testset "makes s and d both ingoing and outgoing neighbors" begin
        @test 4 ∈ neighbors(lg, 5) && 4 ∈ outneighbors(lg, 5) && 4 ∈ inneighbors(lg, 5)
        @test 5 ∈ neighbors(lg, 4) && 5 ∈ outneighbors(lg, 4) && 5 ∈ inneighbors(lg, 4)

        @test 4 ∈ neighbors(lg, 0) && 4 ∈ outneighbors(lg, 0) && 4 ∈ inneighbors(lg, 0)
        @test 0 ∈ neighbors(lg, 4) && 0 ∈ outneighbors(lg, 4) && 0 ∈ inneighbors(lg, 4)
    end
end

@testset "Adding (s, d) edge to directed LabelledGraph" begin
    lg = LabelledGraph{SimpleDiGraph}(["a", "b", "c", "d"])

    add_edge!(lg, "a", "c")
    add_edge!(lg, "d", "a")

    @testset "does not add (d, s) edge" begin
        @test has_edge(lg, "a", "c") && !has_edge(lg, "c", "a")
        @test has_edge(lg, "d", "a") && !has_edge(lg, "a", "d")
    end

    @testset "does not add any other edge" begin
        @test ne(lg) == length(edges(lg)) == 2
    end

    @testset "makes s and d neighbors (as reported by all_neighbors)" begin
        @test "a" ∈ all_neighbors(lg, "c") && "c" ∈ all_neighbors(lg, "a")
        @test "d" ∈ all_neighbors(lg, "a") && "a" ∈ all_neighbors(lg, "d")
    end

    @testset "makes s an ingoing neighbor of d and d outgoing neighbor of s" begin
        @test "a" ∈ inneighbors(lg, "c") && "c" ∈ outneighbors(lg, "a")
        @test "d" ∈ inneighbors(lg, "a") && "a" ∈ outneighbors(lg, "d")
    end
end

for source_graph ∈ (path_digraph(5), path_graph(5))
@testset "Adding vertex to LabelledGraph{$source_graph}" begin

    lg = LabelledGraph(["a", "b", "c", "d", "e"], source_graph)

    @testset "is not possible if the label is duplicated" begin
        @test_throws ArgumentError add_vertex!(lg, "a")
    end

    add_vertex!(lg, "f")

    @testset "makes new vertex present in vertices list" begin
        @test "f" ∈ vertices(lg)
    end

    @testset "increases number of vertices by one" begin
        @test nv(lg) == length(vertices(lg)) == 6
    end

    @testset "makes it possible to connect new vertex to previously existing one" begin
        add_edge!(lg, "f", "b")
        @test has_edge(lg, "f", "b")
    end
end
end

for source_graph ∈ (path_digraph(5), path_graph(5))
@testset "Adding multiple vertices to LabelledGraph{$(typeof(source_graph))}" begin

    lg = LabelledGraph(["a", "b", "c", "d", "e"], source_graph)

    @testset "is not possible if any labels are duplicated, in which case no verts are added" begin
        @test_throws ArgumentError add_vertices!(lg, ["f", "g", "h", "a"])
        @test nv(lg) == 5
    end

    new_labels = ["u", "v", "w"]
    add_vertices!(lg, new_labels)

    @testset "makes new vertices present in vertices list" begin
        @test all(label ∈ vertices(lg) for label ∈ new_labels)
    end

    @testset "increases number of vertices by the number of added vertices" begin
        @test nv(lg) == length(vertices(lg)) == 8
    end

    @testset "makes it possible to connect any two of new vertices" begin
        add_edge!(lg, "u", "w")
        add_edge!(lg, "w", "v")
        @test has_edge(lg, "u", "w")
        @test has_edge(lg, "w", "v")
    end

    @testset "makes it possible to connet new vertex and previously existing one" begin
        add_edge!(lg, "w", "a")
        add_edge!(lg, "b", "u")
        @test has_edge(lg, "w", "a")
        @test has_edge(lg, "b", "u")
    end
end
end

@testset "LabelledGraph can be constructed using default-type aliases" begin
    lg = LabelledGraph(["a", "b", "c"])
    @test !is_directed(lg)
    @test nv(lg) == 3

    lg = LabelledDiGraph(["a", "b", "c"])
    @test is_directed(lg)
    @test nv(lg) == 3
end

for graph_type ∈ (MetaGraph, MetaDiGraph)
@testset "LabelledGraph backed up by $graph_type can store metainformation" begin
    lg = LabelledGraph{graph_type}([2, 5, 10])
    add_edge!(lg, 2, 5)
    set_prop!(lg, 2, :x, 10.0)
    set_prop!(lg, 2, 5, :y, "test")
    set_prop!(lg, LabelledEdge(2, 5), :z, [1, 2,3])
    set_prop!(lg, :name, "the ising model")

    @test get_prop(lg, 2, :x) == 10.0
    @test get_prop(lg, LabelledEdge(2, 5), :z) == [1, 2, 3]
    @test get_prop(lg, 2, 5, :y) == "test"
    @test get_prop(lg, :name) == "the ising model"
end
end

for graph_type ∈ (MetaGraph, MetaDiGraph)
@testset "LabelledGraph backed up by $graph_type can store multiple metainformation at once" begin
    lg = LabelledGraph{graph_type}([2, 5, 10])

    vertex_props = Dict(:x => 10.0, :y => -1.0)
    sd_props = Dict(:a => "test", :b => "Fortran")
    edge_props = Dict(:u => [1, 2, 3], :v => 0)
    graph_props = Dict(:name => "the", :title => "ising model")

    add_edge!(lg, 2, 5)
    set_props!(lg, 2, vertex_props)
    set_props!(lg, 2, 5, sd_props)
    set_props!(lg, LabelledEdge(2, 5), edge_props)
    set_props!(lg, graph_props)

    @test props(lg, 2) == vertex_props
    @test props(lg, LabelledEdge(2, 5)) == merge(sd_props, edge_props)
    @test props(lg, 2, 5) == merge(sd_props, edge_props)
    @test props(lg) == graph_props
end
end

for graph_type ∈ (MetaGraph, MetaDiGraph)
@testset "LabelledGraph backed up by $graph_type should support querying for property existence" begin
    lg = LabelledGraph{graph_type}([2, 5, 10])
    add_edge!(lg, 2, 5)
    set_prop!(lg, 2, :x, 10.0)
    set_prop!(lg, 2, 5, :y, "test")
    set_prop!(lg, LabelledEdge(2, 5), :z, [1, 2,3])
    set_prop!(lg, :name, "the ising model")

    @test has_prop(lg, 2, :x)
    @test !has_prop(lg, 2, :z)
    @test has_prop(lg, LabelledEdge(2, 5), :z)
    @test has_prop(lg, 2, 5, :y)
    @test !has_prop(lg, 2, 5, :zenek)
    @test has_prop(lg, :name)
    @test !has_prop(lg, :title)
end
end

for graph_type ∈ (SimpleGraph, SimpleDiGraph, MetaGraph, MetaDiGraph)
@testset "Vertices labels are preserved when constructing induced subgraph of LabelledGraph{$graph_type}" begin
    lg = LabelledGraph{graph_type}([10, 20, 40, 50])
    add_edge!(lg, 40, 10)
    add_edge!(lg, 50, 10)
    add_edge!(lg, 20, 40)

    sub_lg, vmap = induced_subgraph(lg, [10, 20, 40])

    @test nv(sub_lg) == 3
    @test vertices(sub_lg) == [10, 20, 40]

    @test ne(sub_lg) == 2
    @test has_edge(sub_lg, 20, 40)
    @test has_edge(sub_lg, 40, 10)
    @test vmap == [10, 20, 40]
end
end

for graph_type ∈ (MetaGraph, MetaDiGraph)
@testset "Metainformation are correctly preserved when constructing induced subgraph of LabelledGraph{$graph_type}" begin
    lg = LabelledGraph{graph_type}([5, 10, 15])
    add_edge!(lg, 5, 10)

    set_prop!(lg, :name, "The Ising model")
    set_prop!(lg, 5, :x, 20)
    set_prop!(lg, 5, 10, :y, 30)

    sub_ig, vmap = induced_subgraph(lg, [10, 5])
    @test get_prop(sub_ig, :name) == "The Ising model"
    @test get_prop(sub_ig, 5, :x) == 20
    @test get_prop(sub_ig, 5, 10, :y) == 30
end
end


text:
@testset "PSO" begin

    RANDOM_SEED = 909231
    rng = MersenneTwister(RANDOM_SEED)

    f_sphere(x) = sum(x.^2)

    # * Single-run tests for TestFunctions, Easom
    @test begin
        assert_results = []
        ground_truth = @SVector [π, π]
        for k = 1:10
            val = PSO(
                Easom(),
                Population(35, 2, -100.0, 100.0, rng),
                5_000,
                rng
            )
            result = isapprox(val.x, ground_truth)
            push!(assert_results, result)
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Single-run tests for user-defined functions
    @test begin
        assert_results = []
        ground_truth = zeros(30)
        for k = 1:10
            val = PSO(
                f_sphere,
                Population(35, 30, -10.0, 10.0, rng),
                15_000,
                rng
            )
            push!(assert_results, isapprox(val.x, ground_truth, atol=1e-8))
        end
        # if at least 80% of the time converges, the test passes
        if count(assert_results) >= 8
            true
        else
            println(count(assert_results))
            false
        end
    end

    # * Check that kwargs are handled correctly
    @test begin
        ground_truth = zeros(30)
        val = PSO(
            f_sphere,
            Population(35, 30, -10.0, 10.0, rng),
            15_000,
            rng;
            c1=1.0,
            c2=0.5
        )

        isapprox(val.x, ground_truth, atol=1e-8)
    end
end


text:
module TestApp

Base.eval(Main, :(const UserApp = TestApp))

include("genie.jl")

Base.eval(Main, :(const Genie = TestApp.Genie))
Base.eval(Main, :(using Genie))

using Genie, Genie.Router, Genie.Renderer, Genie.AppServer

end

max_stars_repo_name: tetsugps/julia

max_issues_repo_name: tetsugps/julia

max_forks_repo_name: tetsugps/julia


text:
"""
    NamedDimsStyle{S}
This is a `BroadcastStyle` for NamedDimsArray's
It preserves the dimension names.
`S` should be the `BroadcastStyle` of the wrapped type.
"""
struct NamedDimsStyle{S<:BroadcastStyle} <: AbstractArrayStyle{Any} end
NamedDimsStyle(::S) where {S} = NamedDimsStyle{S}()
NamedDimsStyle(::S, ::Val{N}) where {S,N} = NamedDimsStyle(S(Val(N)))
NamedDimsStyle(::Val{N}) where {N} = NamedDimsStyle{DefaultArrayStyle{N}}()
function NamedDimsStyle(a::BroadcastStyle, b::BroadcastStyle)
    inner_style = BroadcastStyle(a, b)

    # if the inner_style is Unknown then so is the outer-style
    if inner_style isa Unknown
        return Unknown()
    else
        return NamedDimsStyle(inner_style)
    end
end

function Base.BroadcastStyle(::Type{<:NamedDimsArray{L,T,N,A}}) where {L,T,N,A}
    inner_style = typeof(BroadcastStyle(A))
    return NamedDimsStyle{inner_style}()
end
function Base.BroadcastStyle(::NamedDimsStyle{A}, ::NamedDimsStyle{B}) where {A,B}
    return NamedDimsStyle(A(), B())
end

# Resolve ambiguities
# for all these cases, we define that we win to be the outer style regardless of order
for B in (
    :BroadcastStyle, :DefaultArrayStyle, :AbstractArrayStyle, :(Broadcast.Style{Tuple}),
)
    @eval function Base.BroadcastStyle(::NamedDimsStyle{A}, b::$B) where A
        return NamedDimsStyle(A(), b)
    end
    @eval function Base.BroadcastStyle(b::$B, ::NamedDimsStyle{A}) where A
        return NamedDimsStyle(b, A())
    end
end

"""
    unwrap_broadcasted

Recursively unwraps `NamedDimsArray`s and `NamedDimsStyle`s.
replacing the `NamedDimsArray`s with the wrapped array,
and `NamedDimsStyle` with the wrapped `BroadcastStyle`.
"""
function unwrap_broadcasted(bc::Broadcasted{NamedDimsStyle{S}}) where {S}
    inner_args = map(unwrap_broadcasted, bc.args)
    return Broadcasted{S}(bc.f, inner_args, axes(bc))
end
unwrap_broadcasted(x) = x
unwrap_broadcasted(nda::NamedDimsArray) = parent(nda)

# We need to implement copy because if the wrapper array type does not support setindex
# then the `similar` based default method will not work
function Broadcast.copy(bc::Broadcasted{NamedDimsStyle{S}}) where {S}
    inner_bc = unwrap_broadcasted(bc)
    data = copy(inner_bc)

    L = broadcasted_names(bc)
    return NamedDimsArray{L}(data)
end

function Base.copyto!(dest::AbstractArray, bc::Broadcasted{NamedDimsStyle{S}}) where {S}
    inner_bc = unwrap_broadcasted(bc)
    copyto!(dest, inner_bc)
    L = unify_names(dimnames(dest), broadcasted_names(bc))
    return NamedDimsArray{L}(dest)
end

broadcasted_names(bc::Broadcasted) = broadcasted_names(bc.args...)
function broadcasted_names(a, bs...)
    a_name = broadcasted_names(a)
    b_name = broadcasted_names(bs...)
    return unify_names_longest(a_name, b_name)
end
broadcasted_names(a::AbstractArray) = dimnames(a)
broadcasted_names(a) = tuple()


text:
using Revise

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end
using Revise, Finometrics, Random

try #seems to be a bug with running projects locally
  using Production
catch err
  using Production
end

Random.seed!(11)

const REDEFINE = true

#this is to avoid annoying redefinition warnings when possible
if REDEFINE || !(@isdefined DEFINED)
  const STATES_PER_VAR = 5
  const TOL = 1.0e-10
  const D_TOL = 1e-6
  const MAX_ITER = 10000

  #***************competitive parameters**************************
  const COM_τ = [0.2, 0.3, 0.3]
  const COM_δ = 0.3
  const COM_Θ = (0.8, 0.2, 0.2)
  const COM_α = 1.05
  const COM_β = 0.99
  const COM_γ = Float64(1/3)
  const COM_ζ = (1.0, 0.0)

  const COM_p = [0.81 0.68 0.68; 0.19 0.0 0.0; 0.0 0.32 0.32]
  const COM_M = 3

  #some heuristics for the maximum values
  const COM_KMAX = (COM_Θ[1]/(COM_p[1,1]^(1/COM_α)*COM_δ))^(1/(1-COM_γ))
  const COM_YMAX = COM_KMAX^COM_γ*COM_Θ[1]/COM_p[1,2]^(1/COM_α)

  const COM_YHSTATES = collect(range(TOL, stop=COM_YMAX, length=STATES_PER_VAR))
  const COM_KSTATES = collect(range(TOL, stop=COM_KMAX, length=STATES_PER_VAR))

  #***************monopolistic parameters**************************
  const MON_τ = COM_τ
  const MON_δ = COM_δ
  const MON_Θ = COM_Θ
  const MON_α = COM_α
  const MON_β = COM_β
  const MON_γ = COM_γ
  const MON_ζ = (3.0, 1.0)

  const MON_p = COM_p
  const MON_M = COM_M

  #technically this sould be based on a sepearte
  const MON_YHSTATES = COM_YHSTATES
  const MON_KSTATES = COM_KSTATES

  const DEFINED = true
end

function productionscript(firmtype::Symbol)::Nothing
  local F::Firm
  local Φ::Problem

  #the problem parameters
  local τ::Vector{Float64}
  local δ::Float64
  local Θ::NTuple{3,Float64}
  local α::Float64
  local β::Float64
  local γ::Float64
  local kstates::Vector{Float64}
  local yhstates::Vector{Float64}
  local ζ::NTuple{2,Float64}
  local p::Matrix{Float64}

  #assign values based on the scenario
  if firmtype ∈ (:competitive, :COM, :com)
    τ= COM_τ
    δ= COM_δ
    Θ= COM_Θ
    α= COM_α
    β= COM_β
    γ= COM_γ
    kstates= COM_KSTATES
    yhstates= COM_YHSTATES
    ζ= COM_ζ
    p= COM_p
  elseif firmtype ∈ (:monpolistic, :MON, :mon)
    τ = MON_τ
    δ = MON_δ
    Θ = MON_Θ
    α = MON_α
    β= MON_β
    γ= MON_γ
    kstates = MON_KSTATES
    yhstates = MON_YHSTATES
    ζ= MON_ζ
    p = MON_p
  else
    error("firmtype $firmtype not recognized")
  end

  println("\n#**********RUNNING $firmtype using $STATES_PER_VAR states/var***********#")
  F = Firm(τ=τ, δ=δ, Θ=Θ, α=α, β=β, γ=γ, kstates=kstates, yhstates=yhstates, ζ=ζ, p=p)
  Φ = Problem(F, p, dtol = D_TOL, maxiter = MAX_ITER)

  return nothing
end

@time productionscript(:COM)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
@testset "Pre Whitening filter" begin
    test_signal = complex.(randn(1000, 2), randn(1000, 2)) / sqrt(2)
    Rxx = calc_variance_covariance(test_signal)
    pre_whitening_filter = calc_prewhitening_filter(Rxx)
    @test pre_whitening_filter ≈ I atol = 0.1

    test_signal = complex.(randn(1000, 2), randn(1000, 2)) / sqrt(2) * sqrt(100) # Power: 100
    Rxx = calc_variance_covariance(test_signal)
    pre_whitening_filter = calc_prewhitening_filter(Rxx)
    @test pre_whitening_filter ≈ I atol = 0.1
end

@testset "Amplitude filter" begin
    test_signal = complex.(randn(1000, 2), randn(1000, 2)) / sqrt(2)
    Rxx = calc_variance_covariance(test_signal)
    amplitude_filter = calc_amplitude_filter(Rxx)
    @test amplitude_filter ≈ I atol = 0.1

    test_signal = complex.(randn(1000, 2), randn(1000, 2)) / sqrt(2) * sqrt(100) # Power: 100
    Rxx = calc_variance_covariance(test_signal)
    pre_whitening_filter = calc_amplitude_filter(Rxx)
    @test pre_whitening_filter ≈ I atol = 0.1
end

@testset "Eigen beamformer" begin
    test_signal = complex.(randn(1000, 4), randn(1000, 4)) / sqrt(2) / 5 .+ ones(1000, 4)
    Rxx = calc_variance_covariance(test_signal)
    beamformer = calc_eigen_beamformer(Rxx)
    @test beamformer ≈ -1 * [0.5, 0.5, 0.5, 0.5] atol = 0.03
end


text:
using Test
using WriteVTK

@testset "Extent" begin
    @testset "N = 2" begin
        Ns = (6, 8)
        ext = (1:6, 3:10)
        @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 0"
        @test WriteVTK.extent_attribute(ext) == "0 5 2 9 0 0"
    end

    @testset "N = 3" begin
        Ns = (6, 8, 42)
        ext = (1:6, 1:8, 1:3)
        ext = (1:6, 1:8, 2:43)
        @test WriteVTK.extent_attribute(ext) == "0 5 0 7 1 42"
        @test WriteVTK.extent_attribute(Ns) == "0 5 0 7 0 41"
    end

    @testset "N = 4" begin
        Ns = (3, 4, 5, 6)
        @test_throws ArgumentError WriteVTK.extent_attribute(Ns)
    end
end


text:
"""
$(TYPEDEF)

Implements the Block Urn Design (BUD) of [Zhao and Wang (2011)](https://doi.org/10.1016/j.cct.2011.08.004).

The BUD has conditional allocation probability

```math
P_{im} = \\frac{w_m\\lambda + w_mk_{i-1}-n_{i-1,m}}{W\\lambda + Wk_{i-1} - (i - 1)},
\\quad k_{i-1} = \\min_{m=1,...,M} \\left\\lfloor \\frac{n_{i-1,k}}{w_k} \\right\\rfloor.
```

# Fields
$(TYPEDFIELDS)
"""
struct BlockUrn <: MultiArmRandomisationModel
    target::AbstractVector{<:AbstractFloat}
    blocksize::Integer
    weights::AbstractVector{<:Integer}
    W::Integer
    sequence::AbstractVector{<:Integer}
    dist::AbstractVector{<:Integer}
    BlockUrn(
        target::AbstractVector{<:AbstractFloat},
        blocksize::Integer,
        weights::AbstractVector{<:Integer},
        W::Integer,
        sequence::AbstractVector{<:Integer},
        dist::AbstractVector{<:Integer},
    ) = new(target, blocksize, weights, W, sequence, dist)
end

"""
$(TYPEDSIGNATURES)
"""
function BlockUrn(weights::AbstractVector{<:Integer}, blocksize::Int)
    if !all(weights .≥ zero(eltype(weights)))
        throw(DomainError(weights, "All values must be ≥ 0"))
    end
    target = normalize(weights, 1)
    W = sum(weights)
    if mod(blocksize, W) != 0
        throw(DomainError(blocksize, "blocksize not integer multiple of sum(weights)."))
    end
    BlockUrn(target, blocksize, weights, W, zeros(Int, 0), zeros(Int, length(target)))
end

"""
$(TYPEDSIGNATURES)
"""
function BlockUrn(target::AbstractVector{<:AbstractFloat}, blocksize::Int)
    if !isposvec(target)
        throw(DomainError(target, "All values must be ≥ 0"))
    end
    if !isnormvec(target)
        target = target ./ sum(target)
    end
    weights = round.(Int, target .* blocksize)
    W = sum(weights)
    if mod(blocksize, W) != 0
        errstr = "blocksize not integer multiple of sum(target*blocksize)."
        throw(DomainError(blocksize, errstr))
    end
    BlockUrn(target, blocksize, weights, W, zeros(Int, 0), zeros(Int, length(target)))
end

Base.deepcopy(m::BlockUrn) = BlockUrn(deepcopy(m.weights), deepcopy(m.blocksize))

function show(io::IO, BU::BlockUrn)
    k = length(BU.target)
    n = length(BU.sequence)
    print(io, "BlockUrn(arms = $k, blocksize = $(blocksize(BU)), n = $n)")
end

blocksize(BU::BlockUrn) = BU.blocksize
weights(BU::BlockUrn) = BU.weights
target(BU::BlockUrn) = BU.target

"""
$(TYPEDSIGNATURES)
"""
function prob(BU::BlockUrn)
    b = blocksize(BU)
    w = weights(BU)
    W = sum(w)
    λ = b ÷ W
    n = dist(BU)
    i = sum(n) + 1
    k = minimum(floor.(Int, n ./ w)) # minimal balanced sets in previous assignments
    return (λ * w + w * k - n) / (W * λ + W * k - (i - 1))
end

"""
$(TYPEDSIGNATURES)
"""
function update!(BU::BlockUrn, y::Int)
    push!(BU.sequence, y)
    BU.dist[y] += 1
    return nothing
end


text:
function infoPredio(x_predio, y_predio)

x_min = minimum(x_predio);
y_min = minimum(y_predio);
x_predio = x_predio-x_min*ones(size(x_predio));
y_predio = y_predio-y_min*ones(size(y_predio));
V_predio = [x_predio y_predio];

theta, cr = calculaAnguloRotacion(x_predio, y_predio);
numLados=length(theta);

R  = Array{RotInfo,1}([RotInfo(poly2D.rotationMatrix(theta[1]),cr[1,:],theta[1])]);
for i=2:numLados
    push!(R, RotInfo(poly2D.rotationMatrix(theta[i]),cr[i,:],theta[i]));
end

V_predio, R

end


text:
    @testset "Test IterType" begin
        # const IterType = Union{NTuple{N,Tuple} where N, Vector, Matrix, Base.Iterators.Repeated}
        @test isa([(1,2),(1,2)], NDTools.IterType)
        @test isa(((1,2),(1,2)), NDTools.IterType)
        @test isa([1,2,3,4], NDTools.IterType)
        @test isa([1 2;3 4], NDTools.IterType)
    end

    @testset "Test default_type" begin
        Ts = [Bool, Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128]

        for T in Ts
            @test NDTools.default_type(T, Float16) == Float16
        end

        @test NDTools.default_type(Int32, Float32) == Float32
        @test NDTools.default_type(Float64, Float32) == Float64
        @test NDTools.default_type(ComplexF64, Float32) == ComplexF64
        @test NDTools.default_type(Float32, Float64) == Float32  # Should this really be like this?
    end

    @testset "Test default_type" begin
        @test apply_tuple_list(+, 1, 2) == 3
        @test apply_tuple_list(+, [1,2,3], 2) == (3,4,5)
        @test apply_tuple_list(.+, 1, [(1,2),(3,4)]) == ((2, 3), (4, 5))
    end

    @testset "Test get_complex_datatype" begin
        @test get_complex_datatype(22.2) == ComplexF64
        @test get_complex_datatype(22.2f0) == ComplexF32
        @test get_complex_datatype(22) == Complex{Int64}
        @test get_complex_datatype(1f0 + 1f0 * 1im) == Complex{Float32}
        @test get_complex_datatype([1f0 + 1f0 * 1im]) ==Complex{Float32}
    end


text:
using Test
let ps = Pipe[]
    ulimit_n = tryparse(Int, readchomp(`sh -c 'ulimit -n'`))
    try
        for i = 1:100*something(ulimit_n, 1000)
            p = Pipe()
            Base.link_pipe!(p)
            push!(ps, p)
        end
        if ulimit_n === nothing
            @warn "`ulimit -n` is set to unlimited, fd exhaustion cannot be tested"
            @test_broken false
        else
            @test false
        end
    catch ex
        isa(ex, Base.IOError) || rethrow()
        @test ex.code in (Base.UV_EMFILE, Base.UV_ENFILE)
    finally
        foreach(close, ps)
    end
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
## KafkaProducer

mutable struct KafkaProducer
    client::KafkaClient
    topics::Dict{String, KafkaTopic}
end

function KafkaProducer(conf::Dict; dr_cb=nothing, err_cb=nothing)
    kc = KafkaClient(KAFKA_TYPE_PRODUCER, conf; dr_cb=dr_cb, err_cb=err_cb)
    return KafkaProducer(kc, Dict())
end

function KafkaProducer(bootstrap_servers::String, conf::Dict=Dict(); dr_cb=nothing, err_cb=nothing)
    conf["bootstrap.servers"] = bootstrap_servers
    return KafkaProducer(conf; dr_cb=dr_cb, err_cb=err_cb)
end

function Base.show(io::IO, p::KafkaProducer)
    bootstrap_servers = p.client.conf["bootstrap.servers"]
    print(io, "KafkaProducer($bootstrap_servers)")
end

function produce(kt::KafkaTopic, partition::Integer, key, payload)
    # produce(kt.rkt, partition, convert(Vector{UInt8}, key), convert(Vector{UInt8}, payload))
    produce(kt.rkt, partition, Vector{UInt8}(key), Vector{UInt8}(payload))
end

function produce(p::KafkaProducer, topic::String, partition::Integer, key, payload)
    if !haskey(p.topics, topic)
        p.topics[topic] = KafkaTopic(p.client, topic, Dict())
    end
    produce(p.topics[topic], partition, key, payload)
end

function produce(p::KafkaProducer, topic::String, key, payload)
    partition_unassigned = -1
    produce(p, topic, partition_unassigned, key, payload)
end


text:
export
    ####
    #### mjModel.jl
    ####

    mjModel,
    mjLROpt,
    mjVFS,
    mjOption,
    # TODO export these?
    mjVisual_global,
    mjVisual_quality,
    mjVisual_headlight,
    mjVisual_map,
    mjVisual_scale,
    mjVisual_rgba,
    mjVisual,
    mjStatistic,
    mjModel,

    ####
    #### mjData.jl
    ####

    mjData,
    mjContact,
    mjWarningStat,
    mjTimerStat,
    mjSolverStat,

    ####
    #### mjrender.jl
    ####

    mjrRect,
    mjrContext,

    ####
    #### mjui.jl
    ####

    mjuiState,
    mjuiThemeSpacing,
    mjuiThemeColor,
    mjuiItem,
    mjuiSection,
    mjUI,
    mjuiDef,

    ####
    #### mjvisualize.jl
    ####

    mjvPerturb,
    mjvCamera,
    mjvGLCamera,
    mjvGeom,
    mjvLight,
    mjvOption,
    mjvScene,
    mjvFigure,

    ####
    #### mjmujoco.jl
    ####

    mj_Euler,
    mj_RungeKutta,
    mj_activate,
    mj_addContact,
    mj_addFileVFS,
    mj_addM,
    mj_applyFT,
    mj_camlight,
    mj_certAnswer,
    mj_certCheck,
    mj_certQuestion,
    mj_checkAcc,
    mj_checkPos,
    mj_checkVel,
    mj_collision,
    mj_comPos,
    mj_comVel,
    mj_compareFwdInv,
    mj_constraintUpdate,
    mj_contactForce,
    mj_copyData,
    mj_copyModel,
    mj_crb,
    mj_deactivate,
    mj_defaultLROpt,
    mj_defaultOption,
    mj_defaultSolRefImp,
    mj_defaultVFS,
    mj_defaultVisual,
    mj_deleteData,
    mj_deleteFileVFS,
    mj_deleteModel,
    mj_deleteVFS,
    mj_differentiatePos,
    mj_energyPos,
    mj_energyVel,
    mj_factorM,
    mj_findFileVFS,
    mj_forward,
    mj_forwardSkip,
    mj_freeLastXML,
    mj_fullM,
    mj_fwdAcceleration,
    mj_fwdActuation,
    mj_fwdConstraint,
    mj_fwdPosition,
    mj_fwdVelocity,
    mj_getTotalmass,
    mj_id2name,
    mj_integratePos,
    mj_invConstraint,
    mj_invPosition,
    mj_invVelocity,
    mj_inverse,
    mj_inverseSkip,
    mj_isDual,
    mj_isPyramidal,
    mj_isSparse,
    mj_jac,
    mj_jacBody,
    mj_jacBodyCom,
    mj_jacGeom,
    mj_jacPointAxis,
    mj_jacSite,
    mj_kinematics,
    mj_loadModel,
    mj_loadXML,
    mj_local2Global,
    mj_makeConstraint,
    mj_makeData,
    mj_makeEmptyFileVFS,
    mj_mulJacTVec,
    mj_mulJacVec,
    mj_mulM,
    mj_mulM2,
    mj_name2id,
    mj_normalizeQuat,
    mj_objectAcceleration,
    mj_objectVelocity,
    mj_passive,
    mj_printData,
    mj_printModel,
    mj_projectConstraint,
    mj_ray,
    mj_rayHfield,
    mj_rayMesh,
    mj_referenceConstraint,
    mj_resetCallbacks,
    mj_resetData,
    mj_resetDataDebug,
    mj_resetDataKeyframe,
    mj_rne,
    mj_rnePostConstraint,
    mj_saveLastXML,
    mj_saveModel,
    mj_sensorAcc,
    mj_sensorPos,
    mj_sensorVel,
    mj_setConst,
    mj_setLengthRange,
    mj_setTotalmass,
    mj_sizeModel,
    mj_solveM,
    mj_solveM2,
    mj_step,
    mj_step1,
    mj_step2,
    mj_subtreeVel,
    mj_tendon,
    mj_transmission,
    mj_version,
    mj_warning,
    mjr_addAux,
    mjr_blitAux,
    mjr_blitBuffer,
    mjr_changeFont,
    mjr_defaultContext,
    mjr_drawPixels,
    mjr_figure,
    mjr_findRect,
    mjr_finish,
    mjr_freeContext,
    mjr_getError,
    mjr_makeContext,
    mjr_maxViewport,
    mjr_overlay,
    mjr_readPixels,
    mjr_rectangle,
    mjr_render,
    mjr_restoreBuffer,
    mjr_setAux,
    mjr_setBuffer,
    mjr_text,
    mjr_uploadHField,
    mjr_uploadMesh,
    mjr_uploadTexture,
    mju_Halton,
    mju_clearHandlers,
    mju_d2n,
    mju_decodePyramid,
    mju_encodePyramid,
    mju_error,
    mju_error_i,
    mju_error_s,
    mju_f2n,
    mju_free,
    mju_insertionSort,
    mju_isBad,
    mju_isZero,
    mju_malloc,
    mju_muscleBias,
    mju_muscleDynamics,
    mju_muscleGain,
    mju_n2d,
    mju_n2f,
    mju_printMat,
    mju_printMatSparse,
    mju_rayGeom,
    mju_raySkin,
    mju_springDamper,
    mju_standardNormal,
    mju_str2Type,
    mju_strncpy,
    mju_type2Str,
    mju_user_error_cb,
    mju_user_warning_cb,
    mju_warning,
    mju_warningText,
    mju_warning_i,
    mju_warning_s,
    mju_writeLog,
    mjv_addGeoms,
    mjv_alignToCamera,
    mjv_applyPerturbForce,
    mjv_applyPerturbPose,
    mjv_averageCamera,
    mjv_cameraInModel,
    mjv_cameraInRoom,
    mjv_defaultCamera,
    mjv_defaultFigure,
    mjv_defaultOption,
    mjv_defaultPerturb,
    mjv_defaultScene,
    mjv_freeScene,
    mjv_frustumHeight,
    mjv_initGeom,
    mjv_initPerturb,
    mjv_makeConnector,
    mjv_makeLights,
    mjv_makeScene,
    mjv_model2room,
    mjv_moveCamera,
    mjv_moveModel,
    mjv_movePerturb,
    mjv_room2model,
    mjv_select,
    mjv_updateCamera,
    mjv_updateScene,
    mjv_updateSkin


text:
function plot(track::Track)

    lx = track.boundary_left[:,1]
    ly = track.boundary_left[:,2]
    rx = track.boundary_right[:,1]
    ry = track.boundary_right[:,2]

    lt = [Tuple(track.boundary_left[ii,:]) for ii in 1:length(track.boundary_left[:,1])]
    rt = [Tuple(track.boundary_right[ii,:]) for ii in length(track.boundary_right[:,1]):-1:1]
    tuple_list = [lt; rt]

    trk_shape = Plots.Shape(tuple_list)

    plt = plot(trk_shape,
            fillcolor=plot_color(:grey, 1.0),
            lw=0,
            legend=:none,
            grid=false,
            background_color=palette(:algae)[50],
            size=(1500, 1500),
            xlim=(minimum([lx; rx])-10, maximum([lx; rx])+10),
            ylim=(minimum([ly; ry])-10, maximum([ly; ry])+10),
    )
    plot!(plt, lx, ly, color=:black, lw=3, legend=:none)
    plot!(plt, rx, ry, color=:black, lw=3, legend=:none)
    plot!(plt, aspect_ratio=:equal)
    return plt

end

function plot(
    env::CarRacingEnv;
    car_only::Bool=false,
    car_color::Int=1,
    text_output=false,
    text_xy=(80.0, -60.0),
)
    if mod(car_color, 6) == 0 car_color += 1 end
    x, y, ψ = env.state[1:3]
    l_f = env.params.l_f
    l_r = env.params.l_r
    w = (l_f + l_r)*0.4

    rot_mat = [cos(ψ) -sin(ψ) ;
               sin(ψ)  cos(ψ)
              ]

    fl = rot_mat * [l_f; w/2] + [x; y]
    fr = rot_mat * [l_f; -w/2] + [x; y]
    rl = rot_mat * [-l_r; w/2] + [x; y]
    rr = rot_mat * [-l_r; -w/2] + [x; y]

    ar = rot_mat * [-l_r*0.8; 0] + [x; y]
    r = (l_f+l_r) * 0.8
    u, v = r * cos(ψ), r * sin(ψ)

    car_xy = vcat(fl', fr', rr', rl', fl')

    if !car_only
        plt = plot(env.track)
    end

    plt = plot!(car_xy[:,1], car_xy[:,2],
        linestype = :path,
        linewidth = 2,
        linecolor = palette(:lighttest)[mod1(car_color, 7)],
        legend = false
    )
    arrow0!(ar[1], ar[2] , u, v, as=0.35, lc=:black)

    if text_output && car_color == 1
        out_string = @sprintf("%s\n%-24s:%5.1f",
            "", "Time (s)", env.t*env.dt)
        out_string = @sprintf("%s\n\n%8s%13s%9s",
            out_string, "", "--- CAR 1 ---", "")
        out_string = @sprintf("%s\n%-24s:%5.1f",
            out_string, "Reward", reward(env))
        out_string = @sprintf("%s\n%-24s:%5.1f",
            out_string, "|v| (m/s)", norm(env.state[4:5]))
        out_string = @sprintf("%s\n%-24s:%5.1f",
            out_string, "Dist to Lane Center (m)", within_track(env).dist)
        out_string = @sprintf("%s\n%-24s:%5.1f",
            out_string, "Beta Angle (deg)", rad2deg(calculate_β(env)))
        out_string = @sprintf("%s\n%-24s:%5.1f",
            out_string, "Steering Ang (deg)", rad2deg(env.state[7]))
        out_string = @sprintf("%s\n%-24s:%5.1f",
            out_string, "Gas/Brake (%)", env.state[8]*100)
        annotate!(text_xy[1], text_xy[2], text(out_string, "courier bold", :left, :black, 14))
    end
    return plt
end

function plot(
    env::CarRacingEnv,
    pol::AbstractPathIntegralPolicy,
    perc=1.0;
    text_output=false,
    text_xy=(80.0, -60.0),
)
    K = pol.params.num_samples
    trajs = pol.logger.trajectories
    traj_weights = pol.logger.traj_weights

    mod_fac = 10 / perc
    mod_fac = round(Int, mod_fac)

    order = sortperm(traj_weights, rev=true)

    p = plot(env, text_output=text_output, text_xy=text_xy)

    for (ii, k) ∈ enumerate(order)
        col_idx = (K - ii + 1) / K
        if mod(ii*10, mod_fac) == 0 && ii != 1
            x = trajs[k][:,1]
            y = trajs[k][:,2]
            p = plot!(x,y,
                linewidth=0.05,
                linecolor=cgrad(:RdYlGn_9)[col_idx]
            )
        end
    end
    return p
end

function plot(env::MultiCarRacingEnv; text_output=false, text_xy=(80.0, -60.0))
    p = plot(env.envs[1], car_color=1, text_output=text_output)
    for ii in 2:env.N
        p = plot(env.envs[ii], car_only=true, car_color=ii, text_output=text_output, text_xy=text_xy)
    end
    return p
end

function plot(
    env::MultiCarRacingEnv,
    pol::AbstractPathIntegralPolicy,
    perc=1.0;
    text_output=false,
    text_xy=(80.0, -60.0),
)
    K = pol.params.num_samples
    trajs = pol.logger.trajectories
    traj_weights = pol.logger.traj_weights

    mod_fac = 10 / perc
    mod_fac = round(Int, mod_fac)

    order = sortperm(traj_weights, rev=true)

    p = plot(env, text_output=text_output, text_xy=text_xy)

    for (ii, k) ∈ enumerate(order)
        col_idx = (K - ii + 1) / K
        if mod(ii*10, mod_fac) == 0 && ii != 1
            for jj ∈ 1:env.N
                idx = round(Int, (jj-1)*pol.params.ss/env.N)
                x = trajs[k][:,1+idx]
                y = trajs[k][:,2+idx]
                p = plot!(x,y,
                    linewidth=0.05,
                    linecolor=cgrad(:RdYlGn_9)[col_idx]
                )
            end
        end
    end
    return p
end

function plot(env::DroneEnv; drone_only=false)
    x, y, z = env.state[1:3]
    vx, vy, vz = env.state[4:6]
    L = 2*env.params.L

    ψ = atan(vy,vx)
    rot_mat = [cos(ψ) -sin(ψ) ;
               sin(ψ)  cos(ψ)
              ]

    fl = rot_mat * [L; L] + [x; y]
    fr = rot_mat * [L; -L] + [x; y]
    rl = rot_mat * [-L; L] + [x; y]
    rr = rot_mat * [-L; -L] + [x; y]

    ar = rot_mat * [-L*0.8; 0] + [x; y]
    r = 2*L*0.8
    u, v = r * cos(ψ), r * sin(ψ)

    car_xy = vcat(fl', fr', rr', rl', fl')

    if !drone_only
        plt = plot(env.track)
    end

    plot!(car_xy[:,1], car_xy[:,2], linestype=:path, linewidth=2, linecolor=:red, legend=false)
    arrow0!(ar[1], ar[2] , u, v, as=0.35, lc=:black)
    return plt
end

# as: arrow head size 0-1 (fraction of arrow length)
# la: arrow alpha transparency 0-1
function arrow0!(x, y, u, v; as=0.07, lc=:black, la=1)
    nuv = sqrt(u^2 + v^2)
    v1, v2 = [u;v] / nuv,  [-v;u] / nuv
    v4 = (3*v1 + v2)/3.1623  # sqrt(10) to get unit vector
    v5 = v4 - 2*(v4'*v2)*v2
    v4, v5 = as*nuv*v4, as*nuv*v5
    plot!([x,x+u], [y,y+v], lc=lc,la=la)
    plot!([x+u,x+u-v5[1]], [y+v,y+v-v5[2]], lc=lc, la=la)
    plot!([x+u,x+u-v4[1]], [y+v,y+v-v4[2]], lc=lc, la=la)
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SPRAL"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/ralna/spral.git", "cd2d2e817275f16d586cf72767631b3c7472ce02")
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/spral
if [[ "${target}" == *-freebsd* ]] || [[ "${target}" == *-apple-* ]]; then
    CC=gcc
    CXX=g++
fi
./autogen.sh
mkdir build
cd build
CFLAGS=-fPIC CPPFLAGS=-fPIC CXXFLAGS=-fPIC FFLAGS=-fPIC FCFLAGS=-fPIC \
    ../configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-blas="-L${libdir} -lopenblas" --with-lapack="-L${libdir} -lopenblas" \
    --with-metis="-L${libdir} -lmetis" --with-metis-inc-dir="${prefix}/include"
make
gfortran -fPIC -shared -Wl,$(flagon --whole-archive) libspral.a -Wl,$(flagon --no-whole-archive | cut -d' ' -f1) -lgomp -lopenblas -lhwloc -lmetis -lstdc++ -o ${libdir}/libspral.${dlext}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libspral", :libspral),
    ExecutableProduct("spral_ssids", :spral_ssids)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="METIS_jll", uuid="d00139f3-1899-568f-a2f0-47f597d42d70"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
    Dependency(PackageSpec(name="OpenBLAS32_jll", uuid="656ef2d0-ae68-5445-9ca0-591084a874a2"))
    Dependency(PackageSpec(name="Hwloc_jll", uuid="e33a78d0-f292-5ffc-b300-72abe9b543c8"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9.1.0", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
__precompile__()

module Playground

using Compat
using ArgParse
using Mocking
using FilePaths
using Memento
using YAML

import YAML: load

include("constants.jl")
include("config.jl")
include("env.jl")
include("repl.jl")
include("shell.jl")
include("utils.jl")
include("parsing.jl")
include("install.jl")
include("create.jl")
include("activate.jl")
include("execute.jl")
include("list.jl")
include("clean.jl")

export
    # methods
    main,
    argparse,
    install,
    create,
    activate,
    deactivate,
    execute,
    list,
    clean,

    # Constants
    DEFAULT_CONFIG,

    # Types
    Config,
    Environment

const logger = get_logger(current_module())
const cache = Vector{Dict{Symbol, Any}}()

function main(cmdargs, configargs...)
    args = argparse(cmdargs)

    cmd = args["%COMMAND%"]
    log_level = args["debug"] ? "debug" : "info"
    Memento.config(log_level; fmt="[ {level} ] {msg}")

    args = args[cmd]

    config = Config(configargs...)

    debug(logger, "Arguments: $args")
    if cmd == "install"
        install_cmd = args["%COMMAND%"]
        args = args[install_cmd]

        if install_cmd == "download"
            install(
                config,
                args["version"];
                labels=args["labels"],
            )
        elseif install_cmd == "link"
            install(
                config,
                abs(args["dir"]);
                labels=args["labels"],
            )
        # elseif install_cmd == "build"
        #     error("Building from source isn't supported yet.")
        end
    elseif cmd == "create"
        create(
            config,
            args["dir"],
            args["name"];
            julia=args["julia-version"],
            reqs_file=args["requirements"],
            registry=args["registry"],
            branch=args["branch"],
        )
    elseif cmd == "activate"
        activate(
            config,
            args["dir"],
            args["name"],
        )
    elseif cmd == "exec"
        execute(
            `$(Base.shell_split(args["cmd"]))`,
            config,
            args["dir"],
            args["name"],
        )
    elseif cmd == "list"
        list(
            config;
            show_links=args["show-links"],
        )
    elseif cmd == "clean"
        clean(config)
    elseif cmd == "rm"
        rm(
            config;
            name=args["name"],
            dir=args["dir"],
        )
    end
end

end


text:
"""
    PieceType

`DataType` representing a type of piece.
"""
struct PieceType
    val::UInt8
end

# Define all piecetypes
const PIECETYPES = [:VOID, :PAWN, :KNIGHT, :BISHOP, :ROOK, :QUEEN, :KING]

for (num, pt) in enumerate(PIECETYPES)
    @eval const $pt = PieceType($num - 1)
end

"""
    Piece

`DataType` representing a specific piece with type and color.
"""
struct Piece
    val::UInt8
end

# Define all pieces
const WHITE_PIECES = [:WHITEPAWN, :WHITEKNIGHT, :WHITEBISHOP, :WHITEROOK, :WHITEQUEEN, :WHITEKING]
const BLACK_PIECES = [:BLACKPAWN, :BLACKKNIGHT, :BLACKBISHOP, :BLACKROOK, :BLACKQUEEN, :BLACKKING]
const BLANK = Piece(0)

for (num, p) in enumerate(WHITE_PIECES)
    @eval const $p = Piece($num + 8)
end

for (num, p) in enumerate(BLACK_PIECES)
    @eval const $p = Piece($num + 16)
end

"""
    type(piece::Piece)

Return the type of a given `piece`.
"""
type(piece::Piece) = PieceType(piece.val & 7)

"""
    PieceStack

`DataType` storing a vector of PieceTypes.
"""
mutable struct PieceStack <: AbstractArray{PieceType, 1}
    list::Vector{PieceType}
    idx::Int
end

# Allows a preallocation for PieceStack
PieceStack(size::Int) = PieceStack(repeat([VOID], size), 0)

# define useful array methods for PieceStack
Base.iterate(m::PieceStack, state = 1) = (state > m.idx) ? nothing : (m.list[state], state + 1)
Base.length(m::PieceStack) = m.idx
Base.eltype(::Type{PieceStack}) = PieceType
Base.size(m::PieceStack) = (m.idx, )
Base.IndexStyle(::Type{<:PieceStack}) = IndexLinear()
Base.getindex(m::PieceStack, idx::Int) = m.list[idx]

# add moves to the PieceStack
function push!(m::PieceStack, piece::PieceType)
    m.idx += 1
    @inbounds m.list[m.idx] = piece
end

"""
    Color

`DataType` representing a color.
"""
struct Color
    val::UInt8
end

const COLORS = [:NONE, :WHITE, :BLACK]
for (num, c) in enumerate(COLORS)
    @eval const $c = Color($num - 1)
end

"""
    !(color::Color)

Find the opposite color.
"""
(!)(color::Color) = Color(4 >> color.val)

"""
    color(piece::Piece)

Return the `Color` type of a given `piece`.
"""
color(piece::Piece) = Color(piece.val >> 3)

"""
    makepiece(piece::Piece, color::Color)

Return a `piece` of a given `color`.
"""
makepiece(type::PieceType, color::Color) = Piece((color.val << 3) | type.val)

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
#ConventionalApp: Install interface (deprecated).
#-------------------------------------------------------------------------------

#==Top-level interface
===============================================================================#

"""
`install(appurl::String, destpath::String)`

Download/install specified application (.git repository) to `destpath`.
"""
function install(appurl::String, destpath::String)
	@info("Cloning $appurl...")
	println()
	repo = LibGit2.clone(appurl, destpath)
	proj = Project(LibGit2.path(repo))
	return proj
end

#TODO: Install in pwd() by default??
#install(appurl::String) = install(appurl, pwd())

"""
`update!(repo::LibGit2.GitRepo)`

Update application code (.git repository).
"""
function update!(repo::LibGit2.GitRepo)
	LibGit2.fetch(repo)
	LibGit2.merge!(repo)
end

"""
update_repo!(proj::Project)

Update application code (.git repository).
More explicit name than `update!()` to minimize accidental use.
"""
function update_repo!(proj::Project)
	repo = LibGit2.GitRepo(proj.root)
end

max_stars_repo_name: ma-laforge/JuliaApplications

max_issues_repo_name: ma-laforge/JuliaApplications

max_forks_repo_name: ma-laforge/JuliaApplications


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Hinton diagram generation using text.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-12-08

# TODO: Add a legend?

immutable HintonTextStyle
    bins::Vector{Char}
    pos ::String
    neg ::String
    bg  ::String
    function HintonTextStyle(syms, pos, neg, bg; mirror_syms=true)
        new(mirror_syms ? vcat(reverse(syms[2:end]), syms) : bins, pos, neg, bg)
    end
end
bins(o::HintonTextStyle) = o.bins
pos(o::HintonTextStyle)  = o.pos
neg(o::HintonTextStyle)  = o.neg
bg(o::HintonTextStyle)   = o.bg
const TEXTSTYLE = HintonTextStyle(
    [' ', '·', '▪', '■', '█'],
    "\033[97m",
    "\033[30m",
    "\033[47m",
    )

immutable HintonTextDiagram{T<:Real}
    data ::Matrix{T}
    style::HintonTextStyle
end
data(o::HintonTextDiagram)  = o.data
style(o::HintonTextDiagram) = o.style
import Base: show
function show(io::IO, o::HintonTextDiagram)
    w                = data(o)
    sty              = style(o)
    numbins          = length(bins(sty))
    (minval, maxval) = extrema(w)
    (m, n)           = size(w)
    for j in 1:m
        for i in 1:n
            v      = w[j, i]
            bin    = @compat let
                round(Int, (v - minval)/(maxval - minval)*(numbins - 1)) + 1
            end
            sym    = bins(sty)[bin]
            colour = bin < (numbins + 1)÷2 ? neg(sty) : pos(sty)
            print(io, "\033[1m$(bg(sty))$colour$sym\033[0m")
        end
        j >= m || println(io)
    end
end

function hintontxt{T<:Real}(w::AbstractMatrix{T},
        sty::HintonTextStyle=TEXTSTYLE)
    HintonTextDiagram(copy(w), sty)
end


text:
using EVQueues, Plots, Random

lambda=120.0;
mu=1.0;
gamma=0.5;

C=60.0;

Tfinal=1000.0;

sim = ev_edf(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
compute_statistics!(sim)

sim2 = ev_lifo(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
compute_statistics!(sim2)

att1 = [ev.requestedEnergy - ev.departureWorkload for ev in sim.EVs[30000:end]]
att2 = [ev.requestedEnergy - ev.departureWorkload for ev in sim2.EVs[30000:end]]

scatter(att1,att2)


text:
#=

    This is a port of the R2 model CoinBias.cs

    Output from the R2 model:
    ```
    Mean: 0.421294
    Variance: 0.0162177
    Number of accepted samples = 692
    ```

    This model:
  parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        bias    0.4166    0.1360     0.0014    0.0020   5021.6561    1.0000     2276.3627

=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function coin_bias(x)
    n = length(x)
    # Beta(2,5) has mean about 0.2855
    bias ~ Beta(2,5)
    x ~ filldist(Bernoulli(bias),n)
end

x = parse.(Int,split(readline("coin_bias.txt"),","))
println("x:$x")
model = coin_bias(x)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

show_var_dist_pct(chns,:bias,20)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Random
using STL
using Test
using TestAbstractTypes

include("test-StdMap.jl")
include("test-StdSharedPtr.jl")
include("test-StdString.jl")
include("test-StdVector.jl")


text:
using PDDL
using Revise
using mutex

domain = load_domain("test/pddl/domain.pddl")
problem = load_problem("test/pddl/hard_problem.pddl")

graph = create_graph(domain, problem; max_levels=1000);
# acts = get_all_actions(domain, problem)
1


text:
mutable struct MoschopoulosParameters{T}
        θ₁::T
        C::T
        to_power::AbstractVector{T}
        γ::AbstractVector{T}
        ρ::T
        δ::AbstractVector{T}
end
function MoschopoulosParameters(α,θ)
    T = Base.promote_eltype(α,θ,[1.0])
    α = T.(α)
    θ = T.(θ)
    δ = T.([1])
    θ₁ = Base.minimum(θ)
    C = exp(sum(α .* log.(θ₁ ./ θ)))
    ρ = sum(α)
    to_power = (-θ₁ ./ θ) .+1
    γ = [sum(α .* to_power)] # gamma1
    return MoschopoulosParameters(θ₁,C,to_power,γ,ρ,δ)
end
Base.eltype(::MoschopoulosParameters{T}) where T = T

"""
    UnivariateGammaConvolution(α,θ)

Constructs a distribution that corresponds to the convolutions of Gamma(α[i],θ[i]) distributions.
The distribution can then be used through several methods,
following the Distributions.jl standard, to obtain pdf, cdf, random samples...

The pdf and cdf are handled by the Moshopoulos algorithm, and random samples by simply adding random gammas.
The code is type stable and handles any <:Real types, given by the parameters.

To fit the distribution, a loglikelyhood approach could be used. A more involved approach from Furman might be coded sometimes (but requires tanh-sinh integration and bigfloats...)

# Examples
```julia-repl
julia> dist = UnivariateGammaConvolution([1,0.5, 3.7],[4,2, 10]);
julia> sample = zeros(Float64,10);
julia> Random.rand!(dist,sample);
julia> pdf.((dist,),sample);
```
"""
struct UnivariateGammaConvolution{T<:Real} <: Distributions.ContinuousUnivariateDistribution where T
    α::AbstractVector{T}
    θ::AbstractVector{T}
    P::MoschopoulosParameters{T}
end

n(x::UnivariateGammaConvolution) = length(x.α)

function Base.show(io::IO, m::UnivariateGammaConvolution)
    println("Univariate Gamma Convolutions with parametrisation:")
    display([Text.(["α" "θ"]); [m.α m.θ]])
end

# Constructors :
UnivariateGammaConvolution(α::Real,θ::Real) = Distributions.Gamma(α,θ)
function UnivariateGammaConvolution(α::AbstractVector{T1},θ::AbstractVector{T2}) where {T1 <: Real, T2 <: Real}

    # Start by promoving alpha and theta to the same type :
    T = Base.promote_eltype(α,θ,[1.0]) # At least float.
    α = T.(α)
    θ = T.(θ)
    tol = eps(T) # Arbitrary fixed tolerence !

    order = sortperm(θ)
    n = length(θ)
    θ = θ[order]
    α = α[order]
    for i in 1:n
        if i < length(θ)
            for j in (i+1):n
                if j <= length(θ)
                    if abs(θ[i] - θ[j]) <= tol
                        new_α = α[i]+α[j]
                        new_θ = (α[i]*θ[i]+α[j]*θ[j])/(new_α)
                        α[i] = new_α
                        θ[i] = new_θ
                        deleteat!(α,j)
                        deleteat!(θ,j)
                        j = j-1
                    end
                end
            end
        end
    end
    to_be_kept = (α .> tol) .& (θ .> tol) .* (α.*θ .> tol)
    α = α[to_be_kept]
    θ = θ[to_be_kept]

    if length(α) == 1
        return Distributions.Gamma(α[1],θ[1])
    end
    return UnivariateGammaConvolution(α,θ,MoschopoulosParameters(α,θ))
end

#### Support
Distributions.@distr_support UnivariateGammaConvolution 0.0 Inf

#### Conversions
Base.eltype(d::UnivariateGammaConvolution) = typeof(d.α[1])
Base.convert(::Type{UnivariateGammaConvolution{T}}, d::UnivariateGammaConvolution{S}) where {T <: Real, S <: Real} = UnivariateGammaConvolution(T.(d.α), T.(d.θ))

#### parameters
shapes(d::UnivariateGammaConvolution) = d.α
scales(d::UnivariateGammaConvolution) = d.θ
rates(d::UnivariateGammaConvolution) = 1 / d.θ

#### Statistics
mean(d::UnivariateGammaConvolution) = sum(d.α .* d.θ)
var(d::UnivariateGammaConvolution) = sum(d.α .* d.θ .^ 2)

#### Sampling
function Base.rand(rng::Distributions.AbstractRNG,d::UnivariateGammaConvolution)
    sum(rand.(rng,Distributions.Gamma.(d.α,d.θ)))
end

#### Characteristic functions: pdf, cdf, mgf, cf

# Moshopoulos algorithm for pdf and cdf.
function MoschopoulosAlgo!(d::UnivariateGammaConvolution,x::Real, which)

    @assert(which in ["pdf","cdf"], "which should be etiher pdf or cdf")
    T = Base.eltype(d)
    atol = eps(T)
    rtol = T(0)
    entry_type = typeof(x)
    x = T(x)
    if x < T(0)
        return T(0)
    end
    k = 0
    out = T(0)
    while true
        if length(d.P.δ) < k+1
            # then compute the new deltas:
            n = length(d.P.δ)
            for k in n:(k+1)
                pushfirst!(d.P.γ, sum(d.α .* d.P.to_power .^ (k+1)))
                push!(d.P.δ,sum(d.P.γ[2:end] .* d.P.δ)/(k))
            end
        end
        dist = Distributions.Gamma(T(d.P.ρ + k),T(d.P.θ₁))
        if which == "pdf"
            step = d.P.δ[k+1] * Distributions.pdf(dist,x)
        elseif which == "cdf"
            step = d.P.δ[k+1] * Distributions.cdf(dist,x)
        end
        # if ((!isfinite(step)) & (x > T(0)))
        #     print("x = ", x)
        #     print(d)
        #     error("mince")
        # end
        @assert(!((!isfinite(step)) & (x > T(0))),"inf or nan append, the algorithm did not converge for x = $x")
        out += step
        if isapprox(step,T(0),atol=atol,rtol=rtol)
            break
        end
        k = k+1
    end
    out *= d.P.C
    return entry_type(out)
end

Distributions.pdf(d::UnivariateGammaConvolution,x::Real) = MoschopoulosAlgo!(d,x,"pdf")
Distributions.cdf(d::UnivariateGammaConvolution, x::Real) = MoschopoulosAlgo!(d,x,"cdf")
Distributions.logpdf(d::UnivariateGammaConvolution, x::Real) = log(pdf(d,x))
Distributions.mgf(d::UnivariateGammaConvolution, t::Real) = prod((1 - t .* d.θ) .^ (-d.α))
Distributions.cf(d::UnivariateGammaConvolution, t::Real) = prod((1 - (im * t) .* d.θ) .^ (-d.α))