text:
# AUTO GENERATED FILE - DO NOT EDIT

export notification

"""
    notification(;kwargs...)
    notification(children::Any;kwargs...)
    notification(children_maker::Function;kwargs...)

A Notification component.
Show dynamic notifications and alerts to user, part of notifications system. For more information, see: https://mantine.dev/core/notification/
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): Notification body, place main text here
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `className` (String; optional): Often used with CSS to style elements with common properties
- `color` (optional): Notification line or icon color
- `loading` (Bool; optional): Replaces colored line or icon with Loader component
- `title` (String; optional): Notification title, displayed before body
"""
function notification(; kwargs...)
        available_props = Symbol[:children, :id, :className, :color, :loading, :title]
        wild_props = Symbol[]
        return Component("notification", "Notification", "dash_mantine_components", available_props, wild_props; kwargs...)
end

notification(children::Any; kwargs...) = notification(;kwargs..., children = children)
notification(children_maker::Function; kwargs...) = notification(children_maker(); kwargs...)

ext: jl


text:
#=
Problems based on:
CEC2021 Real world multi-objective Constrained Optimization Test Suite
See https://github.com/P-N-Suganthan/2021-RW-MOP
=#

const NAME_OF_PROBLEMS_RW_MOP_2021 = [
                       # Mechanical problems
                       "pressure_vessel",
                       "vibrating_platform",
                       "two_bar_Truss_design_problems",
                       "weldan_beam_design",
                       "disc_brake_design",
                       "speed_reducer_design",
                       "gear_train_design",
                       "car_side_impact_design",
                       "four_bar_plane_truss",
                       "two_bar_plane_truss",
                       "water_tesource_management",
                       "simply_supported_I_beam_design",
                       "gear_box_design",
                       "multiple_disk_clutch_brake_design",
                       "spring_design",
                       "cantilever_beam_design",
                       "bulk_carriers_design",
                       "front_rail_design",
                       "multi_product_batch_plant",
                       "hydro_static_thrust_bearing_design",
                       "crash_energy_management_for_high_speed_train",

                       # chemical problems
                       "haverlys_pooling_problem",
                       "reactor_network_design",
                       "heat_exchanger_network_design",

                       # Process Design and Synthesis
                       "process_synthesis",
                       "process_sythesis_and_design",
                       "process_flow_sheeting_problem",
                       "two_reactor_problem",
                       "process_synthesis_problem",

                       # Power electronics Problem
                       "SOPWM_for_3_level_invereters",
                       "SOPWM_for_5_level_Inverters",
                       "SOPWM_for_7_level_inverters",
                       "SOPWM_for_9_level_inverters",
                       "SOPWM_for_11_level_inverters",
                       "SOPWM_for_13_level_inverters",

                       "sizing_of_single_phase_distributed_generation_1",
                       "sizing_of_single_phase_distributed_generation_2",
                       "sizing_of_single_phase_distributed_generation_3",
                       "sizing_of_single_phase_distributed_generation_4",
                       "active_and_reactive_power_loss",
                       "optimal_power_flow_1",
                       "optimal_power_flow_2",
                       "optimal_power_flow_3",
                       "optimal_power_flow_4",
                       "optimal_power_flow_5",
                       "optimal_power_flow_6",
                       "minimization_of_active_power_loss_in_islanded_microgrids_1",
                       "minimization_of_active_power_loss_in_islanded_microgrids_2",
                       "minimization_of_active_power_loss_in_islanded_microgrids_3",
                       "power_distribution_system_planning"
                      ]

function gen_dict()
    # problems_in_issue = [21, 3, 5, 6, 15]

    properties = vcat(
                      # Mechanical problems (21)
                      [Set([ "mechanics", "real-world", "multiobjective"]) for i in 1:21],
                      # chemical problems (3)
                      [Set([ "chemistry", "real-world", "multiobjective"]) for i in 1:3],
                      # Process Design and Synthesis (5)
                      [Set([ "process-design-and-synthesis", "real-world", "multiobjective"]) for i in 1:5],
                      # Power electronics Problem  (6)
                      [Set([ "power-electronics", "real-world", "multiobjective"]) for i in 1:6],
                      # Power System Problems (15)
                      [Set([ "power-systems", "real-world", "multiobjective"]) for i in 1:15]
                     )

    dict = Dict{String, Set{String}}()

    for (name, property) in zip(NAME_OF_PROBLEMS_RW_MOP_2021, properties)
        dict[name] = property
    end

    dict
end

const functionDict_RWMOP_2021 = gen_dict()

function RW_MOP_bounds(prob_k::Int, dimension = 0)

    n = dimension

    ## range
    # bound constraint definitions for all 50 test functions
    xmin  = [ zeros(1) for i in 1:50 ]
    xmax  = [ zeros(1) for i in 1:50 ]
    xmin[1]    = [0.51,0.51,10,10];
    xmax[1]    = [99.49,99.49,200,200];
    xmin[2]    = [0.05,0.2,0.2,0.35,3];
    xmax[2]    = [0.5,0.5,0.6,0.5,6];
    xmin[3]    = [1e-5,1e-5,1];
    xmax[3]    = [100, 100, 3];
    xmin[4]    = [0.125,0.1,0.1,0.125];
    xmax[4]    = [5,10,10,5];
    xmin[5]    = [55,75,1000,11];
    xmax[5]    = [80,110,3000,20];
    xmin[6]    = [2.6,0.7,16.51,7.3,7.3,2.9,5];
    xmax[6]    = [3.6,0.8,28.49,8.3,8.3,3.9,5.5];
    xmin[7]    = [11.51,11.51,11.51,11.51];
    xmax[7]    = [60.49,60.49,60.49,60.49];
    xmin[8]    = [0.5,0.45,0.5,0.5,0.875,0.4,0.4];
    xmax[8]    = [1.5,1.35,1.5,1.5,2.625,1.2,1.2];
    xmin[9]    = [1,sqrt(2),sqrt(2),1];
    xmax[9]    = [3,3,3,3];
    xmin[10]   = [0.1,0.5];
    xmax[10]   = [2,2.5];
    xmin[11]   = [0.01,0.01, 0.01];
    xmax[11]   = [0.45,0.1, 0.1];
    xmin[12]   = [10, 10, 0.9, 0.9];
    xmax[12]   = [80, 50, 5, 5];
    xmin[13]    = [2.6,0.7,16.51,7.3,7.3,2.9,5];
    xmax[13]    = [3.6,0.8,28.49,8.3,8.3,3.9,5.5];
    xmin[14]    =  [60,90, 1, 0, 2];
    xmax[14]    = [80, 110, 3, 1000, 9];
    xmin[15]    = [0.51,0.6,0.51];
    xmax[15]    = [70.49,3,42.49];
    xmin[16]    = [0.01, 0.20];
    xmax[16]    = [0.05, 1];
    xmin[17]    = [150.0, 20.0, 13.0, 10.0, 14.0, 0.63];
    xmax[17]    = [274.32, 32.31, 25.0, 11.71, 18.0, 0.75];
    xmin[18]    = [136,56,1.4];
    xmax[18]    = [146,68,2.2];
    xmin[19]   = [ 0.51,0.51,0.51,250,250,250,6,4,40,10];
    xmax[19]   = [3.49,3.49,3.49,2500,2500,2500,20,16,700,450];
    xmin[20]   = [ 1, 1,  1e-6,1];
    xmax[20]   = [16, 16, 16*1e-6,16];
    xmin[21] = [1.3, 2.5, 1.3, 1.3, 1.3, 1.3];
    xmax[21] = [1.7, 3.5, 1.7, 1.7, 1.7, 1.7];
    xmin[22] = zeros(9);
    xmax[22] = [100,200,100,100,100,100,200,100,200];
    xmin[23] = [0,0,0,0,0.00001,0.00001];
    xmax[23] = [ 1,1,1,1,16,16];
    xmin[24] = [0,0,0,0,1000,0,100,100,100];
    xmax[24] = [10,200,100,200,2000000,600,600,600,900];
    xmin[25] = [0,-0.49];
    xmax[25] = [1.6,1.49];
    xmin[26] = [0.5,0.5,-0.49];
    xmax[26] = [1.4,1.4,1.49];
    xmin[27] = [0.2,-2.22554,-0.49];
    xmax[27] = [1,-1,1.49];
    xmin[28] = [0,0,0,0,-0.49,-0.49,0];
    xmax[28] = [20,20,10,10,1.49,1.49,40];
    xmin[29] = [0,0,0,-0.49,-0.49,-0.49,-0.49];
    xmax[29] = [100,100,100,1.49,1.49,1.49,1.49];
    xmin[30]   = -0*ones(n);
    xmax[30]   = 90*ones(n);
    xmin[31]   = -0*ones(n);
    xmax[31]   = 90*ones(n);
    xmin[32]   = -0*ones(n);
    xmax[32]   = 90*ones(n);
    xmin[33]   = -0*ones(n);
    xmax[33]   = 90*ones(n);
    xmin[34]   = -0*ones(n);
    xmax[34]   = 90*ones(n);
    xmin[35]   = -0*ones(n);
    xmax[35]   = 90*ones(n);

    if 39 >= prob_k >= 36
        xmin[36]   = -1*ones(n);xmin[36][25:28] .= 0;
        xmax[36]   = +1*ones(n);
        xmin[37]   = -1*ones(n);xmin[37][25:28] .= 0;
        xmax[37]   = +1*ones(n);
        xmin[38]   = -1*ones(n);xmin[38][25:28] .= 0;
        xmax[38]   = +1*ones(n);
        xmin[39]   = -1*ones(n);xmin[39][25:28] .= 0;
        xmax[39]   = +1*ones(n);
    end
    if 46 >= prob_k >= 40
        xmin[40]   = -1*ones(n);xmin[40][27:34] .= 0;
        xmax[40]   = +1*ones(n);
        xmin[41]   = -1*ones(n);xmin[41][27:34] .= 0;
        xmax[41]   = +1*ones(n);
        xmin[42]   = -1*ones(n);xmin[42][27:34] .= 0;
        xmax[42]   = +1*ones(n);
        xmin[43]   = -1*ones(n);xmin[43][27:34] .= 0;
        xmax[43]   = +1*ones(n);
        xmin[44]   = -1*ones(n);xmin[44][27:34] .= 0;
        xmax[44]   = +1*ones(n);
        xmin[45]   = -1*ones(n);xmin[45][27:34] .= 0;
        xmax[45]   = +1*ones(n);
        xmin[46]   = -1*ones(n);xmin[46][27:34] .= 0;
        xmax[46]   = +1*ones(n);
    end
    if 49 >= prob_k >= 47
        xmin[47]   = -1*ones(n);xmin[47][11:12] .= 0;xmin[47][13:18] .= 0;
        xmax[47]   = +1*ones(n);xmax[47][11:12] .= 2;xmax[47][13:18] .= 500;
        xmin[48]   = -1*ones(n);xmin[48][11:12] .= 0;xmin[48][13:18] .= 0;
        xmax[48]   = +1*ones(n);xmax[48][11:12] .= 2;xmax[48][13:18] .= 500;
        xmin[49]   = -1*ones(n);xmin[49][11:12] .= 0;xmin[49][13:18] .= 0;
        xmax[49]   = +1*ones(n);xmax[49][11:12] .= 2;xmax[49][13:18] .= 500;
    end

    xmin[50]    = [10, 10, 35, 35, 125, 130];
    xmax[50]  = [125, 150, 210, 225, 315, 325];

    xmin_ = xmin[prob_k]
    xmax_ = xmax[prob_k]

    #nadir = nadir_points[prob_k]

    return xmin_, xmax_
end

function get_RW_MOP_config(function_name::AbstractString)
    fnum = findfirst( s -> function_name == s, NAME_OF_PROBLEMS_RW_MOP_2021)

    isnothing(fnum) && error("$function_name is not implemented yet or does not exist.")

    get_RW_MOP_config(fnum)
end

function get_RW_MOP_config(fnum::Int)
    @assert 1 <= fnum <= 50

    problems_dimension = [4,5,3,4,4,7,4,7,4,2,3,4,7,5,3,2,6,3,10,4,6,9,6,9,2,3,3,7,7,25,
                          25,25,30,30,30,28,28,28,28,34,34,34,34,34,34,34,18,18,18,6];

    problems_objectives = [2,2,2,2,2,2,2,3,2,2,5,2,3,2,2,2,3,2,3,2,2,2,2,3,2,
                           2,2,2,2,2,2,2,2,2,2,2,2,2,3,2,3,2,2,3,3,4,2,2,3,2];

    inequility_constrains = [2,5,3,4,4,11,1,9,1,2,7,1,11,8,8,2,9,3,10,7,4,2,1,0,2,1,3,4,9,
                            24,24,24,29,29,29,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];

    equility_constrains   = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,6,0,1,0,4,0,0,
                            0,0,0,0,0,24,24,24,24,26,26,26,26,26,26,26,12,12,12,1];

    nadir_points = [
                    [359648.85, -7330.3829],
                    [-0.0012746083, 318.25489],
                    [0.1, 100000.0],
                    [36.679325, 0.013066667],
                    [5.3067, 3.0281682],
                    [5969.8221, 1300.0],
                    [3.4655, 45.417795],
                    [92.596587, 4.0, 12.699733],
                    [3048.5281, 0.04],
                    [187.04862, 6.7710178e-5],
                    [73450.511, 1350.0, 2.853469e6, 6.620032e6, 25000.0],
                    [437.95472, 0.061459097],
                    [6103.6015, 1300.0, 1004.6759],
                    [1.3967521, 0.014920745],
                    [27.942662, 187991.19],
                    [3.0630528, 0.0020408763],
                    [-3151.4157, 8260.6298, 812.60004],
                    [0.93366407, 1.196596],
                    [244351.82, 48572.551, 6000.0],
                    [267.25846, -2.7672651e-5],
                    [1.3157336, 26.297736],
                    [-105.90871, 2000.0],
                    [-0.00040194083, 4.0],
                    [6.6395236, -3.6323008e-5, -1.9999995e6],
                    [3.2, -1.25],
                    [2.9263226, -1.5793956],
                    [-0.243, 0.0],
                    [133.43214, 13.1082],
                    [2.999999, 14.0],
                    [0.28549866, 0.10239999],
                    [0.73429846, 0.10239997],
                    [0.58881362, 0.12959996],
                    [0.52171456, 0.10239985],
                    [0.61063821, 0.11108847],
                    [1.9499692, 11.076986],
                    [0.69938497, 0.0040801029],
                    [0.64793246, 0.0057752416],
                    [0.0050223653, 0.0036060355],
                    [0.70181358, 0.0040803309, 0.0057776563],
                    [6.1887155, -0.29920474],
                    [6.1954982, -0.18135402, 4.060485],
                    [0.0068788651, 0.0029695767],
                    [5.7139516, 1.2579325],
                    [5.720257, 1.0231384, 4.0471196],
                    [8.7353948, 2.2861141, 4.6710252],
                    [5.7381081, 2.3578333, 8.1909894, 3.9032935],
                    [0.0093764479, 0.0024737407],
                    [0.14296272, 0.00086761855],
                    [0.14253454, 0.020364141, 0.00058302635],
                    [64895.619, 1253.8936],
                   ]

    nadir = nadir_points[fnum]
    #problem_name = keys(functionDict_RWMOP_2021)[fnum]

    xmin, xmax= RW_MOP_bounds(fnum, problems_dimension[fnum])

    return Dict(:function => NAME_OF_PROBLEMS_RW_MOP_2021[fnum],
                :fn    => problems_objectives[fnum],
                :n     => problems_dimension[fnum],
                :nadir => nadir,
                :gn    => inequility_constrains[fnum],
                :hn    => equility_constrains[fnum],
                :xmin  => xmin,
                :xmax  => xmax)

end

function get_RW_MOP_problem(function_name::AbstractString)
    fnum = findfirst( s -> function_name == s, NAME_OF_PROBLEMS_RW_MOP_2021)

    isnothing(fnum) && error("$function_name is not implemented yet or does not exist.")

    return get_RW_MOP_problem(fnum)

end

function get_RW_MOP_problem(fnum::Int)

    !(1 <= fnum <= 50) && error("Problem $fnum is not implemented yet or does not exist.")

    # objetive function x -> f(x)
    f = eval(Symbol(NAME_OF_PROBLEMS_RW_MOP_2021[fnum]))

    # problem configuration
    conf = get_RW_MOP_config(fnum)

    return f, conf
end


text:
using Test, OffsetArrays
using PooledArrays
using DataAPI: refarray, refvalue, refpool, invrefpool
using PooledArrays: refcount
import Future.copy!

if Threads.nthreads() < 2
    @warn("Running with only one thread: correctness of parallel operations is not tested")
else
    @show Threads.nthreads()
end

@testset "PooledArrays" begin
    @test eltype(PooledArray(Int[])) === Int

    a = rand(10)
    b = rand(10,10)
    c = rand(1:10, 1000)

    @test PooledArray(a) == a
    @test PooledArray(b) == b
    pc = PooledArray(c)
    @test pc == c
    @test convert(Array, pc) == c
    #@test issorted(pc.invpool)

    @test copy(pc) == pc

    @test isa(similar(pc), typeof(pc))

    @test reverse(pc) == reverse(c)

    @test pc[34:250] == c[34:250]

    @test sort(pc) == sort(c)
    @test sortperm(pc) == sortperm(c)

    push!(pc, -10)
    push!(c, -10)
    @test pc == c
    @test maximum(values(pc.invpool)) == 11
    @test length(pc.pool) == 11

    pc2 = copy(pc)
    @test append!(pc2, [3, 1]) === pc2
    @test pc2 == [pc; 3; 1]

    pc2 = copy(pc)
    @test deleteat!(pc2, 2) === pc2
    @test pc2 == pc[[1; 3:end]]
    @test deleteat!(pc2, [2, 4]) === pc2
    @test pc2 == pc[[1; 4; 6:end]]

    #@test issorted(pc.invpool)

    @test map(identity, pc) == pc
    @test map(x->2x, pc) == map(x->2x, c)

    # case where the outputs are one-to-many
    pa = PooledArray([1,2,3,4,5,6])
    @test map(isodd, pa) == [true,false,true,false,true,false]

    px = PooledArray(rand(128))
    py = PooledArray(rand(200))

    @test isa(vcat(px, py), PooledArray{Float64, UInt16})
    @test vcat(px, Array(py)) == vcat(px, py)

    px2 = copy(px)
    resize!(px2, 100)
    @test px2 == px[1:100]

    px3 = PooledArray(px2)
    px3[1] = 0
    @test px2 !== px3
    @test px2 != px3

    @test findall(PooledArray([true,false,true])) == [1,3]

    @test PooledArray{Union{Int,Missing}}([1, 2]) isa PooledArray{Union{Int,Missing}}

    @test eltype(PooledArray(rand(128)).refs) == UInt32
    @test eltype(PooledArray(rand(300)).refs) == UInt32
    @test eltype(PooledArray(rand(128), UInt8).refs) == UInt8
    @test eltype(PooledArray(rand(300), UInt16).refs) == UInt16
    @test PooledVector == PooledArray{T, R, 1} where {T, R}
    @test PooledMatrix == PooledArray{T, R, 2} where {T, R}

    s = PooledArray(["a", "a", "b"])
    @test eltype(PooledArray(s).refs) == UInt32
    @test eltype(PooledArray(s, signed=true).refs) == Int32
    @test eltype(PooledArray(s, compress=true).refs) == UInt8
    @test eltype(PooledArray(s, signed=true, compress=true).refs) == Int8
    @test eltype(PooledArray(rand(300), signed=true, compress=true).refs) == Int16
    @test all(refarray(s) .== [1, 1, 2])
    for i in 1:3
        @test refvalue(s, refarray(s)[i]) == s[i]
    end
    @test refpool(s) == ["a", "b"]
    @test invrefpool(s) == Dict("a" => 1, "b" => 2)

    @testset "push!" begin
        xs = PooledArray([10, 20, 30])
        @test xs === push!(xs, -100)
        @test xs == [10, 20, 30, -100]
    end

    @testset "pushfirst!" begin
        ys = PooledArray([10, 20, 30])
        @test ys === pushfirst!(ys, -100)
        @test ys == [-100, 10, 20, 30]
    end

    v1 = PooledArray([1, 3, 2, 4])
    v2 = PooledArray(BigInt.([1, 3, 2, 4]))
    v3 = PooledArray(["a", "c", "b", "d"])

    @test PooledArrays.fast_sortable(v1) === v1
    @test isbitstype(eltype(PooledArrays.fast_sortable(v1)))
    Base.Order.Perm(Base.Order.Forward, v1).data === v1

    @test PooledArrays.fast_sortable(v2) == PooledArray([1, 3, 2, 4])
    @test isbitstype(eltype(PooledArrays.fast_sortable(v2)))
    Base.Order.Perm(Base.Order.Forward, v2).data == PooledArray([1, 3, 2, 4])

    @test PooledArrays.fast_sortable(v3) == PooledArray([1, 3, 2, 4])
    @test isbitstype(eltype(PooledArrays.fast_sortable(v3)))
    Base.Order.Perm(Base.Order.Forward, v3).data == PooledArray([1, 3, 2, 4])

    for T in (Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64)
        @inferred PooledArray([1, 2, 3], T)
        @inferred PooledArray([1 2 3], T)
    end
    @test typeof(PooledArray([1 2; 3 4], Int8)) === PooledMatrix{Int, Int8, Matrix{Int8}}

    for signed in (true, false), compress in (true, false)
        @test_throws ErrorException @inferred PooledArray([1, 2, 3],
                                                          signed=signed,
                                                          compress=compress)
    end
end

@testset "pool non-copying constructor and copy tests" begin
    pa = PooledArray([1, 2, 3])
    @test pa.refcount[] == 1
    pa2 = PooledArray(pa)
    @test pa.refcount[] == 2
    @test pa.refs == pa2.refs
    @test pa.refs !== pa2.refs
    @test pa.pool === pa2.pool
    @test pa.invpool === pa2.invpool
    @test pa.refcount === pa2.refcount

    pav = @view pa[[3, 1]]

    @test pav == [3, 1]
    @test refarray(pav) isa SubArray{UInt32,1,Array{UInt32,1},Tuple{Array{Int,1}},false}
    @test refpool(pav) === pa.pool
    @test invrefpool(pav) === pa.invpool
    @test_throws BoundsError refvalue(pav, 0)
    @test refvalue(pav, 1) === 1
    @test refvalue(pav, 2) === 2
    @test refvalue(pav, 3) === 3
    @test_throws BoundsError refvalue(pav, 4)

    @test pa.refcount[] == 2
    pa3 = PooledArray(pav)
    @test pa.refcount[] == 3
    @test pa.refs[[3, 1]] == pa3.refs
    @test pa.refs !== pa3.refs
    @test pa.pool === pa3.pool
    @test pa.invpool === pa3.invpool
    @test pa.refcount === pa3.refcount
    pa2 = pa3
    # try to force GC to check finalizer
    GC.gc(); GC.gc(); GC.gc(); GC.gc()
    if pa.refcount[] != 2
        @warn "finalizer of PooledArray not triggered; excess refs: $(pa.refcount[] - 2)"
    end

    pa = PooledArray([1, 2, 3])
    @test pa.refcount[] == 1
    pa2 = copy(pa)
    @test pa.refcount[] == 2
    @test pa.refs == pa2.refs
    @test pa.refs !== pa2.refs
    @test pa.pool === pa2.pool
    @test pa.invpool === pa2.invpool
    @test pa.refcount === pa2.refcount

    pav = @view pa[1]
    pa3 = copy(pav)
    @test pa3 isa PooledArray{Int, UInt32, 0}
    @test pa.refcount[] == 3
    @test refarray(pav) == pa3.refs
    @test pa.refs !== pa3.refs
    @test pa.pool === pa3.pool
    @test pa.invpool === pa3.invpool
    @test pa.refcount === pa3.refcount
end

@testset "test de-referencing on setindex!" begin
    pa = PooledArray([1, 2, 3])
    @test pa.refcount[] == 1
    pa2 = copy(pa)
    @test pa.refcount[] == 2
    old_pool = pa.pool
    old_invpool = pa.invpool
    old_refcount = pa.refcount

    # within pool
    pa2[1] = 3
    @test pa == 1:3
    @test pa.pool === old_pool
    @test pa.invpool === old_invpool
    @test pa.refcount === old_refcount
    @test pa.refcount[] == 2
    @test pa2 == [3, 2, 3]
    @test pa2.pool === old_pool
    @test pa2.invpool === old_invpool
    @test pa2.refcount === old_refcount

    # new value
    pa2[1] = 4
    @test pa == 1:3
    @test pa.pool == old_pool
    @test pa.invpool == old_invpool
    @test pa.refcount == old_refcount
    @test pa.refcount[] == 1
    @test pa2 == [4, 2, 3]
    @test pa2.pool !== old_pool
    @test pa2.invpool !== old_invpool
    @test pa2.refcount !== old_refcount
    @test pa2.refcount[] == 1
end

@testset "copy! tests" begin
    pa = PooledArray(1:4)
    pa1 = pa[1:2]
    pav1 = @view pa[2:3]
    pa2 = PooledArray(fill(2))
    pav2 = @view pa[3]

    pat1 = PooledArray([0, 0])
    copy!(pat1, pa1)
    @test pat1 == pa1
    @test refpool(pat1) === refpool(pa1)
    @test invrefpool(pat1) === invrefpool(pa1)
    @test refcount(pat1) === refcount(pa1)
    @test refcount(pat1)[] == 3

    copy!(pat1, pav1)
    @test pat1 == pav1
    @test refpool(pat1) === refpool(pav1)
    @test invrefpool(pat1) === invrefpool(pav1)
    @test refcount(pat1) === refcount(pav1)
    @test refcount(pat1)[] == 3

    pat2 = PooledArray(fill(0))
    copy!(pat2, pa2)
    @test pat2 == pa2
    @test refpool(pat2) === refpool(pa2)
    @test invrefpool(pat2) === invrefpool(pa2)
    @test refcount(pat2) === refcount(pa2)
    @test refcount(pat2)[] == 2

    copy!(pat2, pav2)
    @test pat2 == pav2
    @test refpool(pat2) === refpool(pav2)
    @test invrefpool(pat2) === invrefpool(pav2)
    @test refcount(pat2) === refcount(pav2)
    @test refcount(pat2)[] == 4
    @test refcount(pa2)[] == 1
end

@testset "correct refcount when treading" begin
    pa = PooledArray([1 2; 3 4])
    x = Vector{Any}(undef, 120)
    Threads.@threads for i in 1:120
        x[i] = copy(pa)
    end
    @test pa.refcount[] == 121
    Threads.@threads for i in 1:61
        @test x[i].refcount === pa.refcount
        x[i][1] = 2
        @test x[i].refcount === pa.refcount
        x[i][1] = 5
        @test x[i].refcount[] == 1
    end
    @test pa.refcount[] == 60
    x = nothing
    # try to force GC to check finalizer
    GC.gc(); GC.gc(); GC.gc(); GC.gc()
    if pa.refcount[] != 1
        @warn "finalizer of PooledArray not triggered; excess refs: $(pa.refcount[] - 1)"
    end
end

@testset "copyto! tests" begin
    pa1 = PooledArray([1, 2, 3])
    pa2 = similar(pa1, 4)
    @test_throws BoundsError copyto!(pa1, pa2)
    copyto!(pa2, pa1)
    @test refpool(pa2) === refpool(pa1)
    @test invrefpool(pa2) === invrefpool(pa1)
    @test refcount(pa2) === refcount(pa1)
    @test refcount(pa2)[] == 2

    pa1 = view(PooledArray([1, 2, 3]), :)
    pa2 = similar(pa1, 4)
    @test_throws BoundsError copyto!(pa1, pa2)
    copyto!(pa2, pa1)
    @test refpool(pa2) === refpool(pa1)
    @test invrefpool(pa2) === invrefpool(pa1)
    @test refcount(pa2) === refcount(pa1)
    @test refcount(pa2)[] == 2

    pa1 = PooledArray([1, 2, 3])
    pa2 = view(similar(pa1, 4), :)
    @test_throws BoundsError copyto!(pa1, pa2)
    copyto!(pa2, pa1)
    @test refpool(pa2) === refpool(pa1)
    @test invrefpool(pa2) === invrefpool(pa1)
    @test refcount(pa2) === refcount(pa1)
    @test refcount(pa2)[] == 2

    pa1 = view(PooledArray([1, 2, 3]), :)
    pa2 = view(similar(pa1, 4), :)
    @test_throws BoundsError copyto!(pa1, pa2)
    copyto!(pa2, pa1)
    @test refpool(pa2) === refpool(pa1)
    @test invrefpool(pa2) === invrefpool(pa1)
    @test refcount(pa2) === refcount(pa1)
    @test refcount(pa2)[] == 2

    pa1 = PooledArray([1, 2, 3])
    pa2 = similar(pa1, 4)
    copyto!(pa2, 1, view(pa1, [1, 1]), 1, 2)
    @test refpool(pa2) === refpool(pa1)
    @test invrefpool(pa2) === invrefpool(pa1)
    @test refcount(pa2) === refcount(pa1)
    @test refcount(pa2)[] == 2
    copyto!(pa2, 3, view(pa1, [1, 1]), 1, 2)
    @test refpool(pa2) === refpool(pa1)
    @test invrefpool(pa2) === invrefpool(pa1)
    @test refcount(pa2) === refcount(pa1)
    @test refcount(pa2)[] == 2
    @test pa2 == [1, 1, 1, 1]

    pa1 = PooledArray([1, 2, 3])
    pa2 = similar(pa1, Float64, 3)
    copyto!(pa2, 1, pa1, 1, 3)
    @test refpool(pa2) !== refpool(pa1)
    @test invrefpool(pa2) !== invrefpool(pa1)
    @test refcount(pa2) !== refcount(pa1)
    @test refcount(pa1)[] == 1
    @test refcount(pa2)[] == 1
    @test pa2 == [1, 2, 3]

    pa1 = PooledArray([1, 2, 3])
    pa2 = similar(pa1, 3)
    pa2 .= 1
    copyto!(pa2, 1, pa1, 1, 3)
    @test refpool(pa2) !== refpool(pa1)
    @test invrefpool(pa2) !== invrefpool(pa1)
    @test refcount(pa2) !== refcount(pa1)
    @test refcount(pa1)[] == 1
    @test refcount(pa2)[] == 1
    @test pa2 == [1, 2, 3]
end

@testset "reverse" begin
    pa1 = PooledArray([1, 2, 3])
    pa2 = reverse(pa1)
    pa3 = reverse(pa2)
    @test pa2 == [3, 2, 1]
    @test pa3 == pa1
    @test refpool(pa1) === refpool(pa2) === refpool(pa3)
    @test invrefpool(pa1) === invrefpool(pa2) === invrefpool(pa3)
    @test refcount(pa1) === refcount(pa2) === refcount(pa3)
    @test refcount(pa1)[] == 3
end

@testset "convert" begin
    pa1 = PooledArray([1, 2, 3])
    @test convert(PooledArray, pa1) === pa1
    @test eltype(convert(PooledArray{Float64}, pa1)) === Float64
    pa1c = convert(PooledArray{Int, UInt64, 1}, pa1)
    @test pa1c isa PooledArray{Int,UInt64,1,Array{UInt64,1}}
    @test pa1c == pa1
    @test !(pa1c isa typeof(pa1))
end

@testset "indexing" begin
    pa = PooledArray([1 2; 3 4])
    @test pa[2, 2, 1] == 4
    @test pa[2, 2] == 4
    @test pa[big(2), 2] == 4
    pav = view(pa, :, :)
    @test pav[2, 2, 1] == 4
    @test pav[2, 2] == 4
    @test pav[big(2), 2] == 4

    @test refcount(pa)[] == 1
    pa2 = pa[[true, false, false, true]]
    @test pa2 == [1, 4]
    @test refpool(pa) === refpool(pa2)
    @test invrefpool(pa) === invrefpool(pa2)
    @test refcount(pa) === refcount(pa2)
    @test refcount(pa)[] == 2
    pa3 = pav[[true, false, false, true]]
    @test pa3 == [1, 4]
    @test refpool(pa) === refpool(pa3)
    @test invrefpool(pa) === invrefpool(pa3)
    @test refcount(pa) === refcount(pa3)
    @test refcount(pa)[] == 3

    # these checks are mostly needed to check for dispatch ambiguities
    @test pa[1] == 1
    @test pa[1, 1] == 1
    @test pa[1, 1, 1] == 1
    @test pa[:] == [1, 3, 2, 4]
    @test pa[1:4] == [1, 3, 2, 4]
    @test pa[collect(1:4)] == [1, 3, 2, 4]
    @test pa[1, 1:2] == [1, 2]
    @test pa[1, [1, 2]] == [1, 2]
    @test pa[1:1, 1:2] == [1 2]
    @test pa[1:1, [1, 2]] == [1 2]
    @test pa[[1], 1:2] == [1 2]
    @test pa[[1], [1, 2]] == [1 2]
    @test pav[1] == 1
    @test pav[1, 1] == 1
    @test pav[1, 1, 1] == 1
    @test pav[:] == [1, 3, 2, 4]
    @test pav[1:4] == [1, 3, 2, 4]
    @test pav[collect(1:4)] == [1, 3, 2, 4]
    @test pav[1, 1:2] == [1, 2]
    @test pav[1, [1, 2]] == [1, 2]
    @test pav[1:1, 1:2] == [1 2]
    @test pav[1:1, [1, 2]] == [1 2]
    @test pav[[1], 1:2] == [1 2]
    @test pav[[1], [1, 2]] == [1 2]

    pav2 = view(PooledArray([1]), 1)
    pa2 = similar(pav2)
    pa2[] = 10

    @test pav2[] == 1
    @test pa2[] == 10
end

@testset "isassigned" begin
    pa1 = PooledArray(["a"])
    pa2 = similar(pa1, 2)
    pa2v = view(pa2, 1)
    @test !isassigned(pa2, 1)
    @test !isassigned(pa2v)
end

@testset "setindex!" begin
    pa = PooledArray([1 2; 3 4])
    pa[1, 1] = 10
    @test pa == [10 2; 3 4]
    @test [pa pa] == [10 2 10 2; 3 4 3 4]
    pa[2] = 1000
    @test pa == [10 2; 1000 4]
    pa[1, :] = [11, 12]
    @test pa == [11 12; 1000 4]
    pa[1:2, 1:1] = [111, 222]
    @test pa == [111 12; 222 4]
    pa[1, 1, 1] = 0
    @test pa == [0 12; 222 4]
    pa[:] = [1 2; 3 4]
    @test pa == [1 2; 3 4]
end

@testset "repeat" begin
    pa1 = PooledArray([1, 2, 3])
    pa2 = repeat(pa1)
    pa3 = repeat(pa1, 2)
    pa4 = repeat(pa1, inner = 2)
    @test pa2 == pa1
    @test pa3 == [1, 2, 3, 1, 2, 3]
    @test pa4 == [1, 1, 2, 2, 3, 3]
    @test refpool(pa1) === refpool(pa2) === refpool(pa3) === refpool(pa4)
    @test invrefpool(pa1) === invrefpool(pa2) === invrefpool(pa3) === invrefpool(pa4)
    @test refcount(pa1) === refcount(pa2) === refcount(pa3) === refcount(pa4)
    @test refcount(pa1)[] == 4

    pa1 = PooledArray(["one", "two"])
    pa2 = repeat(pa1, outer = 3)
    pa3 = repeat(pa1, inner = 3, outer = 2)
    @test pa2 == ["one", "two", "one", "two", "one", "two"]
    @test pa3 == ["one", "one", "one", "two", "two", "two", "one", "one", "one", "two", "two", "two"]

    # missing shouldn't be a problem

    pa1 = PooledArray(["one", missing, "two"])
    pa2 = repeat(pa1, 2)
    pa3 = repeat(pa1, 0)
    @test isequal(pa2, ["one", missing, "two", "one", missing, "two"])
    @test isequal(pa2.pool, ["one", missing, "two"])
    @test size(pa3) == (0,)
    @test isempty(pa3.refs)

    # two dimensional
    pa1 = PooledArray([true false; false true; true true])
    pa2 = repeat(pa1, 2)
    @test pa2 == Bool[1 0; 0 1; 1 1; 1 0; 0 1; 1 1]

    pa1 = PooledArray([1 2; 3 4])
    pa2 = repeat(pa1, inner = (2, 1))
    @test pa2 == [1 2; 1 2; 3 4; 3 4]
end

@testset "map pure tests" begin
    x = PooledArray([1, 2, 3])
    x[3] = 1
    y = map(-, x, pure=true)
    @test refpool(y) == [-1, -2, -3]
    @test y == [-1, -2, -1]

    y = map(-, x)
    @test refpool(y) == [-1, -2]
    @test y == [-1, -2, -1]

    function f()
        i = Ref(0)
        return x -> (i[] -= 1; i[])
    end

    # the order is strange as we iterate invpool which is a Dict
    # and it depends on the version of Julia
    y = map(f(), x, pure=true)
    d = Dict(Set(1:3) .=> -1:-1:-3)
    @test refpool(y) == [d[i] for i in 1:3]
    @test y == [d[v] for v in x]

    y = map(f(), x)
    @test refpool(y) == [-1, -2, -3]
    @test y == [-1, -2, -3]

    x = PooledArray([1, missing, 2])
    y = map(identity, x)
    @test isequal(y, [1, missing, 2])
    @test typeof(y) === PooledVector{Union{Missing, Int}, UInt32, Vector{UInt32}}

    x = PooledArray([1, missing, 2], signed=true, compress=true)
    y = map(identity, x)
    @test isequal(y, [1, missing, 2])
    @test typeof(y) === PooledVector{Union{Missing, Int}, Int8, Vector{Int8}}

    x = PooledArray(fill(1, 200), signed=true, compress=true)
    y = map(f(), x)
    @test y == -1:-1:-200
    @test typeof(y) === PooledVector{Int, Int, Vector{Int}}

    x = PooledArray(reshape(fill(1, 200), 2, :), signed=true, compress=true)
    y = map(f(), x)
    @test y == reshape(-1:-1:-200, 2, :)
    @test typeof(y) === PooledMatrix{Int, Int, Matrix{Int}}

    x = PooledArray(fill("a"), signed=true, compress=true)
    y = map(f(), x)
    @test y == fill(-1)
    @test typeof(y) === PooledArray{Int, Int8, 0, Array{Int8, 0}}

    @static if VERSION >= v"1.6"
        for signed in (true, false), compress in (true, false), len in (1, 100, 1000)
            x = PooledArray(fill(1, len), signed=signed, compress=compress)
            @inferred PooledVector{Int, Int, Vector{Int}} map(identity, x)
        end
    end
end

@testset "insert! test" begin
    x = PooledArray([1, 2, 3])
    @test insert!(x, 2, 10) == [1, 10, 2, 3]
    @test_throws ArgumentError insert!(x, true, 10)
    @test_throws BoundsError insert!(x, 0, 10)
    @test x == [1, 10, 2, 3]
    @test x.pool == [1, 2, 3, 10]
    @test insert!(x, 3, 'c') == [1, 10, 99, 2, 3]
    @test x.pool == [1, 2, 3, 10, 99]
    @test insert!(x, 1, true) == [1, 1, 10, 99, 2, 3]
    @test x.pool == [1, 2, 3, 10, 99]
    @test insert!(x, 7, true) == [1, 1, 10, 99, 2, 3, 1]
    @test_throws BoundsError insert!(x, 9, true)
    @test x == [1, 1, 10, 99, 2, 3, 1]
end

@testset "pop! and popfirst!" begin
    x = PooledArray([1, 2, 3])
    @test pop!(x) == 3
    @test x == [1, 2]
    @test popfirst!(x) == 1
    @test x == [2]
    x = PooledArray(["1", "2", "3"])
    @test pop!(x) == "3"
    @test x == ["1", "2"]
    @test popfirst!(x) == "1"
    @test x == ["2"]
end

@testset "constructor corner cases" begin
    x = Vector{Any}(undef, 3)
    y = PooledArray(x)
    @test y isa PooledArray{Any}
    @test !any(i -> isassigned(y, i), eachindex(y))
    @test all(iszero, y.refs)
    @test isempty(y.pool)
    @test isempty(y.invpool)

    x[2] = "a"
    for v in (x, OffsetVector(x, -5))
        y = PooledArray(v)
        @test y isa PooledArray{Any}
        @test !isassigned(x, 1)
        @test x[2] == "a"
        @test !isassigned(x, 3)
        @test y.refs == [0, 1, 0]
        @test y.pool == ["a"]
        @test y.invpool == Dict("a" => 1)
    end
end


text:
using OrderedCollections

abstract type PlottingRecipe end

Plots.plot(pr::PlottingRecipe, args...; kwargs...) = plot!(plot(), pr, args...; kwargs...)

# =============================================================================
# for a single region

struct ObservationPlottingRecipe1{Tsd, Ted, To, Te, Tp, Tl} <: PlottingRecipe
    startdate ::Tsd
    enddate   ::Ted
    observed  ::To
    expected  ::Te
    predicted ::Tp
    label     ::Tl
end

struct ObservationWithPredictorsPlottingRecipe1{Tsd, Ted, To, Te, Tp, Tl, Psd} <: PlottingRecipe
    startdate ::Tsd
    enddate   ::Ted
    observed  ::To
    expected  ::Te
    predicted ::Tp
    label     ::Tl
    pstartdate::Psd
end

function expected(data, gp, region, label)
    @assert label in ["cases", "hospitalizations", "deaths", "infections"]
    i  = region isa Integer ? region : findfirst(==(region), data.regions)
    dates = data.dates[i]
    if label == "cases"
        values = gp.expected_daily_cases[i]
        # start  = data.turing_data.casemodel.starts[i]
        # stop   = data.turing_data.casemodel.stops[i]
        # dates  =
        return (; dates, values )
    elseif label == "hospitalizations"
        values = gp.expected_daily_hospits[i]
        return (; dates, values )
    elseif label == "deaths"
        values = gp.expected_daily_deaths[i]
        return (; dates, values )
    else
        values = gp.newly_infecteds[i]
        return (; dates, values )
    end
end

function observed(data, region, label)
    @assert label in ["cases", "hospitalizations", "deaths", "infections"]
    i = region isa Integer ? region : findfirst(==(region), data.regions)
    r = region isa Integer ? Regional.regions[region] : region
    dates, values = if label == "cases"
        data.cases.date, data.cases[:,r]
    elseif label == "hospitalizations"
        data.hospit.date, data.hospit[:,r]
    else
        data.deaths.date, data.deaths[:,r]
    end
    return (; dates, values )
end

function predicted(data, gp, region, label)
    @assert label in ["cases", "hospitalizations", "deaths", "infections"]
    i  = region isa Integer ? region : findfirst(==(region), data.regions)
    dates = data.dates[i]
    if label == "cases"
        values = gp.predicted_daily_cases[i]
        # start  = data.turing_data.casemodel.starts[i]
        # stop   = data.turing_data.casemodel.stops[i]
        # dates  =
        return (; dates, values )
    elseif label == "hospitalizations"
        values = gp.predicted_daily_hospits[i]
        return (; dates, values )
    else
        values = gp.predicted_daily_deaths[i]
        return (; dates, values )
    end
end

function startdate(data::CovidSurvey.Data, region, label::String)
    @assert label in ["cases", "hospitalizations", "deaths"]
    i = region isa Integer ? region : findfirst(==(region), data.regions)
    s = if label == "cases"
        data.turing_data.casemodel.starts[i]
    elseif label == "hospitalizations"
        data.turing_data.hospitmodel.starts[i]
    else
        data.turing_data.deathmodel.starts[i]
    end
    return data.dates[i][s]
end

function enddate(data::CovidSurvey.Data, region, label::String)
    @assert label in ["cases", "hospitalizations", "deaths"]
    i = region isa Integer ? region : findfirst(==(region), data.regions)
    s = if label == "cases"
        data.turing_data.casemodel.stops[i]
    elseif label == "hospitalizations"
        data.turing_data.hospitmodel.stops[i]
    else
        data.turing_data.deathmodel.stops[i]
    end
    return data.dates[i][s]
end

function ObservationPlottingRecipe1(data::CovidSurvey.Data, gp, region, label)
    e = expected(data, gp, region, label)
    o = observed(data, region, label)
    p = predicted(data, gp, region, label)
    sd = startdate(data, region, label)
    ed = enddate(data, region, label)
    isnothing(data.predictors) && return ObservationPlottingRecipe1( sd, ed, o, e, p, label)
    pstartdate = Date("2020-11-10")
    return ObservationWithPredictorsPlottingRecipe1( sd, ed, o, e, p, label, pstartdate)
end

function Plots.plot!(p::Plots.Plot, r::ObservationPlottingRecipe1; plot_only_fit=false)
    o  = r.observed
    e  = r.expected
    pr = r.predicted
    ed = r.enddate
    sd = r.startdate

    if !plot_only_fit
        !isnothing(ed) && vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
        !isnothing(ed) && vline!(p, [sd], lab="start observations", lw=2, lc=:black, hover="$sd", ls=:dash)
    end
    plot!(p, o.dates, o.values, α=0.5, lc=:match, lab="observed $(r.label)", c=:midnightblue, lw=4, ylab="cases")
    plot_confidence_timeseries!(p, pr.dates, pr.values; label = "predicted $(r.label)", c=:midnightblue) #Dict(hover=>strdates)
    plot_confidence_timeseries!(p, e.dates, e.values; label = "expected $(r.label)") #Dict(hover=>strdates)
    if plot_only_fit
        xlo = Dates.value.(sd)
        xup = Dates.value.(ed)
        yup = maximum( o.values[o.dates .< ed] ) * 1.1
        xlims!(p, xlo, xup)
        ylims!(p, 0, yup)
    end
    return p
end

function Plots.plot!(p::Plots.Plot, r::ObservationWithPredictorsPlottingRecipe1; plot_only_fit = false)
    o  = r.observed
    e  = r.expected
    pr = r.predicted
    ed = r.enddate
    sd = r.startdate
    psd= r.pstartdate

    if !plot_only_fit
        !isnothing(ed) && vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
        !isnothing(ed) && vline!(p, [sd], lab="start observations", lw=2, lc=:black, hover="$sd", ls=:dash)
    end
    vspan!(p, [psd, ed], label="fit to survey data", α=0.2, fc=:midnightblue)
    bar!(p, o.dates, o.values, α=0.5, lc=:match, lab="observed $(r.label)", c=:midnightblue, lw=4, ylab="cases")
    plot_confidence_timeseries!(p, pr.dates, pr.values; label = "predicted $(r.label)", c=:royalblue ) #Dict(hover=>strdates)
    plot_confidence_timeseries!(p, e.dates, e.values; label = "expected $(r.label)", c=:peru) #Dict(hover=>strdates)
    if plot_only_fit
        xlo = Dates.value.(sd)
        xup = Dates.value.(ed)
        yup = maximum( o.values[o.dates .< ed] ) * 1.1
        xlims!(p, xlo, xup)
        ylims!(p, 0, yup)
    end
    return p
end

# ============================================================================
# latent infections
struct LatentInfectionsPlottingRecipe3{Ted, Te, Tl} <: PlottingRecipe
    enddate   ::Ted
    expected  ::Te
    label     ::Tl
end

function LatentInfectionsPlottingRecipe3(data::CovidSurvey.Data, gp, region, label)
    e = expected(data, gp, region, label)
    ed = data.observations_end
    #isnothing(data.predictors) &&
    return LatentInfectionsPlottingRecipe3( ed, e, label)
    #pstartdate = Date("2020-11-10")
    #return LatentInfectionsPlottingRecipe1( ed, e, p, label, pstartdate)
end

function Plots.plot!(p::Plots.Plot, r::Regional.LatentInfectionsPlottingRecipe3; plot_only_fit=false)
    e  = r.expected
    ed = r.enddate

    if !plot_only_fit
        !isnothing(ed) && vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
    end
    ylabel!("cases")
    plot_confidence_timeseries!(p, e.dates, e.values; label = "latent $(r.label)") #Dict(hover=>strdates)
    if plot_only_fit
        xlo = Dates.value(e.dates[1])
        xup = Dates.value.(ed)
        # yup = maximum( e.values[e.dates .< ed] ) * 1.1
        xlims!(p, xlo, xup)
        # ylims!(p, 0, yup)
    end
    return p
end
# ============================================================================
# seropos

struct SeroPlottingRecipe1{SD, ED, E, O, L} <: Regional.PlottingRecipe
    startdate  ::SD
    enddate    ::ED
    expecteds  ::E
    observeds  ::O
    label     ::L
end

function SeroPlottingRecipe1(data::CovidSurvey.Data, gp, region, label)
    @unpack dates, means, CIs, delay, populations = data.turing_data.seromodel
    i  = region isa Integer ? region : findfirst(==(region), data.regions)
    sd = first(data.dates[i])# startdate(data, region, "hospitalizations")
    ed = Date(data.observations_end)

    observed   = let
        mean = means[i] * 100
        error = [(m-l*100, u*100-m) for (m,(l,u)) in zip(mean,CIs[i])]
        (; dates, mean, error )
    end
    expected = let
        values = gp.expected_seropos[i]  ./ populations[i] .* 100
        dates  = data.dates[i] + Day(delay)
        (; dates, values)
    end
    SeroPlottingRecipe1( sd, ed, expected, observed, label)
end

function Plots.plot!(p::Plots.Plot, r::SeroPlottingRecipe1; plot_only_fit = false)
    sd= r.startdate
    ed= r.enddate
    e = r.expecteds
    o = r.observeds

    plot_confidence_timeseries!(p, e.dates, e.values; r.label)
    y  = o.mean
    ye = o.error
    scatter!(p, o.dates, y, yerror=ye, hover=["$d: $val" for (d, val) in zip(o.dates, y)], lab="serological survey")
    if plot_only_fit
        xlo = Dates.value.(sd)
        xup = Dates.value.(ed)
        yup = maximum([ maximum( v[e.dates .< ed] ) for v in e.values ]) * 1.1
        xlims!(p, xlo, xup)
        ylims!(p, 0, yup)
    end
    return p
end

# ============================================================================
# reproduction number for a single region

struct RtPlottingRecipe1{Tlo, Tst, Ted, Te, Tl} <: Regional.PlottingRecipe
    lockdown  ::Tlo
    startdate ::Tst
    enddate  ::Ted
    expected  ::Te
    label     ::Tl
end

struct RtWithPredictorsPlottingRecipe1{Tlo, Tst, Ted, Te, Tl, Psd} <: Regional.PlottingRecipe
    lockdown  ::Tlo
    startdate ::Tst
    enddate  ::Ted
    expected  ::Te
    label     ::Tl
    pstartdate::Psd
end

function RtPlottingRecipe1(data::CovidSurvey.Data, generated_posterior, region, label)
    i  = region isa Integer ? region : findfirst(==(region), data.regions)
    dates      = data.dates[i]
    values     = generated_posterior.Rts[i]
    expected   = (; dates, values )
    lo = isnothing(data.lockdown) ? nothing : Date(data.lockdown)
    sd = first(data.dates[i])
    ed = Date(data.observations_end)

    isnothing(data.predictors) && return RtPlottingRecipe1( lo, sd, ed, expected, label)
    psd = Date("2020-11-10")
    return RtWithPredictorsPlottingRecipe1( lo, sd, ed, expected, label, psd)
end

function Plots.plot!(p::Plots.Plot, r::RtPlottingRecipe1; plot_only_fit = false)
    e  = r.expected
    sd = r.startdate
    ed = r.enddate
    lo = r.lockdown

    if !plot_only_fit
        !isnothing(ed) && vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
        !isnothing(lo) && vline!(p, [lo], lab="lockdown", lw=2, lc=:black, hover="$lo", ls=:dash)
    end
    plot_confidence_timeseries!(p, e.dates, e.values; label = r.label) #Dict(hover=>strdates)
    if plot_only_fit
        xlo = Dates.value.(sd)
        xup = Dates.value.(ed)
        # yup = maximum( o.values[o.dates .< ed] ) * 1.1
        xlims!(p, xlo, xup)
        # ylims!(p, 0, yup)
    end
    return p
end

function Plots.plot!(p::Plots.Plot, r::RtWithPredictorsPlottingRecipe1; plot_only_fit = false)
    e  = r.expected
    sd = r.startdate
    ed = r.enddate
    lo = r.lockdown
    psd= r.pstartdate

    if !plot_only_fit
        vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
        vline!(p, [lo], lab="lockdown", lw=2, lc=:black, hover="$lo", ls=:dash)
    end
    vspan!(p, [psd, ed], label="fit to survey data", α=0.2, fc=:midnightblue)
    plot_confidence_timeseries!(p, e.dates, e.values; label = r.label) #Dict(hover=>strdates)
    if plot_only_fit
        xlo = Dates.value.(sd)
        xup = Dates.value.(ed)
        # yup = maximum( e.values[e.dates .< ed] ) * 1.1
        xlims!(p, xlo, xup)
        # ylims!(p, 0, yup)
    end
    return p
end
# ============================================================================
# region plot

struct RegionPlottingRecipe{Tr,Tt,T} <: Regional.PlottingRecipe
    recipes::Tr
    titles ::Tt
    region ::T
end

posterior2recipe = OrderedDict(
    :expected_daily_cases   => ObservationPlottingRecipe1,
    :expected_daily_hospits => ObservationPlottingRecipe1,
    :expected_daily_deaths  => ObservationPlottingRecipe1,
    :Rts                    => RtPlottingRecipe1,
    # :iar                    => IARPlottingRecipe,
    # :effective_Rt           => RtPlottingRecipe1,
    :expected_seropos       => SeroPlottingRecipe1,
)

posterior2label = OrderedDict(
    :expected_daily_cases   => "cases",
    :expected_daily_hospits => "hospitalizations",
    :expected_daily_deaths  => "deaths",
    :Rts                    => "reproduction number",
    # :iar                   => "infections ascertainment rate",
    # :effective_Rt          => "effective reproduction number",
    :expected_seropos      => "total infected"
)

function RegionPlottingRecipe(data::CovidSurvey.Data, generated_posterior, region)
    ks = keys(generated_posterior)
    recipes = Vector{PlottingRecipe}()
    titles  = Vector{String}()
    for (k,recipe) in posterior2recipe
        if k in ks
            label = posterior2label[k]
            r     = recipe(data, generated_posterior, region, label)
            title = label * " for the $region"
            push!(recipes, r)
            push!(titles, title)
        end
    end
    RegionPlottingRecipe(recipes, titles, region)
end

function Plots.plot(r::RegionPlottingRecipe; kwargs...)
    plots  = Vector{Plots.Plot}()
    nplots = length(r.recipes)
    for (recipe, title) in zip(r.recipes, r.titles)
        p = plot(; xaxis = true, legend = :outertopright, title)
        plot!(p, recipe; kwargs...)
        push!(plots, p)
    end
    plot(plots..., layout=(nplots,1), size=(1000, nplots*250), sharex=true, link=:x)
end

# Plots.plot(data::CovidSurvey.Data, generated_posterior) =
#     plot( RegionPlottingRecipe(data, generated_posterior) )
# Plots.plot(data::CovidSurvey.Data, generated_posterior) = plot!(plot(), data, generated_posterior)

# ============================================================================
# reproduction number for all regions

struct RtsPlottingRecipe{Tlo, Tsd, Ted, Tre, Te} <: Regional.PlottingRecipe
    lockdown  ::Tlo
    startdates::Tsd
    enddate   ::Ted
    regions   ::Tre
    expecteds ::Te
end

function RtsPlottingRecipe(data::CovidSurvey.Data, generated_posterior)
    dates      = data.dates
    samples    = generated_posterior.Rts
    expecteds  = [(dates = d, values = g) for (d,g) in zip(dates, samples)]
    lockdown   = isnothing(data.lockdown) ? nothing : Date(data.lockdown) # NOTE use data["lockdown"] instead
    startdates = first.(data.dates)
    enddate    = Date(data.observations_end) # NOTE use data["observations_end"] instead
    RtsPlottingRecipe( lockdown, startdates, enddate, data.regions, expecteds)
end

function Plots.plot(r::RtsPlottingRecipe; kwargs...)
    plots  = Vector{Plots.Plot}()
    nplots = length(r.regions)
    for (expected, region, startdate) in zip(r.expecteds, r.regions, r.startdates)
        recipe = RtPlottingRecipe1(r.lockdown, startdate, r.enddate, expected, "reproduction number")
        p = plot(; xaxis = true, legend = :outertopright, title="$region")
        plot!(p, recipe; kwargs...)
        push!(plots, p)
    end
    plot(plots..., layout=(nplots,1), size=(1000, nplots*200), sharex=true, link=:x)
end

# ============================================================================
# observations for all regions

struct RegionsOverviewPlottingRecipe{Tr,Tt,T} <: Regional.PlottingRecipe
    recipes::Tr
    titles ::Tt
    region ::T
end

function RegionsOverviewPlottingRecipe(data::CovidSurvey.Data, generated_posterior, posterior)
    # i  = region isa Integer ? region : findfirst(==(region), data.regions)
    recipes = Vector{PlottingRecipe}()
    titles  = Vector{String}()
    for (i,region) in enumerate(data.regions)
        label  = posterior2label[posterior]
        recipe = posterior2recipe[posterior]
        r      = recipe(data, generated_posterior, region, label)
        title  = i == 1 ? label*" for the $region" : "$region"
        push!(recipes, r)
        push!(titles, title)
    end
    return RegionsOverviewPlottingRecipe(recipes, titles, data.regions)
end

function Plots.plot(r::RegionsOverviewPlottingRecipe; kwargs...)
    # plots  = Vector{Plots.Plot}()
    nplots = length(r.recipes)

    title = r.titles[1]
    p = plot(; xaxis = true, legend = :outertopright, title)
    plot!(p, r.recipes[1]; kwargs...)
    plots = [p]

    for (recipe, title) in Iterators.drop(zip(r.recipes, r.titles),1)
        p = plot(; xaxis = true, legend = nothing, title)
        plot!(p, recipe; kwargs...)
        push!(plots, p)
    end
    plot(plots..., layout=(nplots,1), size=(1000, nplots*250), sharex=true, link=:x)
end

# # ============================================================================
# # reproduction number for a single region
#
# struct SingleRtPlottingRecipe1{Tlo, Ted, Te} <: Regional.PlottingRecipe
#     lockdown  ::Tlo
#     enddates  ::Ted
#     expected  ::Te
# end
#
# function SingleRtPlottingRecipe1(data::CovidSurvey.Data, generated_posterior, region)
#     i  = region isa Integer ? region : findfirst(==(region), data.regions)
#     dates      = data.dates[i]
#     values     = generated_posterior.Rts[i]
#     expected   = (; dates, values )
#     lockdown   = Date(data.lockdown)
#     enddates   = Date(data.observations_end)
#     SingleRtPlottingRecipe1( lockdown, enddates, expected)
# end
#
# function Plots.plot!(p::Plots.Plot, rt::SingleRtPlottingRecipe1)
#     e  = rt.expected
#     ed = rt.enddates
#     lo = rt.lockdown
#
#     vline!(p, [ed], lab="end observations", lw=2, lc=:black, hover="$ed")
#     vline!(p, [lo], lab="lockdown", lw=2, lc=:black, hover="$lo", ls=:dash)
#     plot_confidence_timeseries!(p, e.dates, e.values; label = "reproduction number") #Dict(hover=>strdates)
# end
#
# Plots.plot(pr::SingleRtPlottingRecipe1) = plot!(plot(), pr)

# # ============================================================================
# # region plot
#
# struct RegionPlottingRecipe{Tr,Tt,T} <: Regional.PlottingRecipe
#     recipes::Tr
#     titles ::Tt
#     region ::T
# end
#
# posterior2recipe = OrderedDict(
#     :expected_daily_hospits => Regional.SingleHospitPlottingRecipe,
#     :expected_daily_deaths  => Regional.SingleDeathsPlottingRecipe,
#     :Rts                    => Regional.SingleRtPlottingRecipe1
# )
#
# posterior2title = OrderedDict(
#     :expected_daily_hospits => "daily hospitalizations",
#     :expected_daily_deaths  => "daily deaths",
#     :Rts                    => "effective reproduction number"
# )
#
# function RegionPlottingRecipe(data::CovidSurvey.Data, generated_posterior, region)
#     ks = keys(generated_posterior)
#     recipes = Vector{Regional.PlottingRecipe}()
#     titles  = Vector{String}()
#     for (k,v) in posterior2recipe
#         if k in ks
#             recipe = v(data, generated_posterior, region)
#             title  = posterior2title[k] * " for the $region"
#             push!(recipes, recipe)
#             push!(titles, title)
#         end
#     end
#     RegionPlottingRecipe(recipes, titles, region)
# end
#
# function Plots.plot(r::RegionPlottingRecipe)
#     plots  = Vector{Plots.Plot}()
#     nplots = length(r.recipes)
#     for (recipe, title) in zip(r.recipes, r.titles)
#         p = plot(; xaxis = true, legend = :outertopright, title)
#         plot!(p, recipe)
#         push!(plots, p)
#     end
#     plot(plots..., layout=(nplots,1), size=(1000, nplots*250), sharex=true, link=:x)
# end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
module PythonIO

include("utils.jl")
include("PyObjects/PyObjects.jl")
include("PyExprIO/PyExprIO.jl")
include("PickleIO/PickleIO.jl")
include("NumpyIO/NumpyIO.jl")

import .PyExprIO: readpyexpr
import .PickleIO: readpkl
import .NumpyIO: readnpy, writenpy, readnpz, writenpz

export readpyexpr, readnpy, writenpy, readnpz, writenpz, readpkl

end


text:
function WriteCℓ!(Probes::String, Cℓ::AbstractCℓ, Filename::String)
    h5write(Filename*".h5", "cls/"*Probes*"/c_lij",
    Cℓ.CℓArray)
end

function WriteWeightFunctions!(Probes::String, WA::AbstractWeightFunction, Filename::String)
    h5write(Filename*".h5", "weight_functions/"*Probes*"/w_bin_z",
    WA.WeightFunctionArray)
end

#TODO Probably those read-write dictionaries functions can be rewritten as a couple of
#macros!

function WriteCosmology!(Cosmology::w0waCDMCosmology, Filename::String)
    CosmoDict = Dict{String,Float64}()
    CosmoDict["w0"]  = Cosmology.w0
    CosmoDict["wa"]  = Cosmology.wa
    CosmoDict["Mν"]  = Cosmology.Mν
    CosmoDict["H0"]  = Cosmology.H0
    CosmoDict["ΩM"]  = Cosmology.ΩM
    CosmoDict["ΩB"]  = Cosmology.ΩB
    CosmoDict["ΩDE"] = Cosmology.ΩDE
    CosmoDict["Ωk"]  = Cosmology.Ωk
    CosmoDict["Ωr"]  = Cosmology.Ωr
    CosmoDict["ns"]  = Cosmology.ns
    CosmoDict["σ8"]  = Cosmology.σ8
    JSON3.write(Filename*"/cosmology.json", CosmoDict)
end

function WriteCosmology!(Cosmology::Flatw0waCDMCosmology, Filename::String)
    CosmoDict = Dict{String,Float64}()
    CosmoDict["w0"]  = Cosmology.w0
    CosmoDict["wa"]  = Cosmology.wa
    CosmoDict["Mν"]  = Cosmology.Mν
    CosmoDict["H0"]  = Cosmology.H0
    CosmoDict["ΩM"]  = Cosmology.ΩM
    CosmoDict["ΩB"]  = Cosmology.ΩB
    CosmoDict["ns"]  = Cosmology.ns
    CosmoDict["σ8"]  = Cosmology.σ8
    JSON3.write(Filename*"/cosmology.json", CosmoDict)
end

function WriteParameters!(CosmoDict::Dict, Filename::String)
    JSON3.write(Filename*"/parameters.json", CosmoDict)
end

function ReadCosmology(CosmoDict::Dict)
    if "ΩDE" in keys(CosmoDict)
        Cosmology = w0waCDMCosmology(
        w0 = CosmoDict["w0"],
        wa = CosmoDict["wa"],
        Mν = CosmoDict["Mν"],
        H0 = CosmoDict["H0"],
        ΩM = CosmoDict["ΩM"],
        ΩB = CosmoDict["ΩB"],
        ΩDE = CosmoDict["ΩDE"],
        Ωk = CosmoDict["Ωk"],
        Ωr = CosmoDict["Ωr"],
        ns = CosmoDict["ns"],
        σ8 = CosmoDict["σ8"])
    else
        Cosmology = Flatw0waCDMCosmology(
        w0 = CosmoDict["w0"],
        wa = CosmoDict["wa"],
        Mν = CosmoDict["Mν"],
        H0 = CosmoDict["H0"],
        ΩM = CosmoDict["ΩM"],
        ΩB = CosmoDict["ΩB"],
        ns = CosmoDict["ns"],
        σ8 = CosmoDict["σ8"])
    end
    return Cosmology
end

function ReadCosmologyForecast(CosmoDict::Dict, CosmoModel::String)
    if "ΩDE" in keys(CosmoDict)
        Cosmology = w0waCDMCosmology(
        w0 = CosmoDict["w0"][1],
        wa = CosmoDict["wa"][1],
        Mν = CosmoDict["Mν"][1],
        H0 = CosmoDict["H0"][1],
        ΩM = CosmoDict["ΩM"][1],
        ΩB = CosmoDict["ΩB"][1],
        ΩDE = CosmoDict["ΩDE"][1],
        Ωk = CosmoDict["Ωk"][1],
        Ωr = CosmoDict["Ωr"][1],
        ns = CosmoDict["ns"][1],
        σ8 = CosmoDict["σ8"][1])
    else
        Cosmology = Flatw0waCDMCosmology(
            w0 = CosmoDict["w0"][1],
            wa = CosmoDict["wa"][1],
            Mν = CosmoDict["Mν"][1],
            H0 = CosmoDict["H0"][1],
            ΩM = CosmoDict["ΩM"][1],
            ΩB = CosmoDict["ΩB"][1],
            ns = CosmoDict["ns"][1],
            σ8 = CosmoDict["σ8"][1])
    end
    return Cosmology
end

function ReadCosmology(CosmoDict::JSON3.Object)
    if "ΩDE" in keys(CosmoDict)
        Cosmology = w0waCDMCosmology(
        w0 = CosmoDict["w0"],
        wa = CosmoDict["wa"],
        Mν = CosmoDict["Mν"],
        H0 = CosmoDict["H0"],
        ΩM = CosmoDict["ΩM"],
        ΩB = CosmoDict["ΩB"],
        ΩDE = CosmoDict["ΩDE"],
        Ωk = CosmoDict["Ωk"],
        Ωr = CosmoDict["Ωr"],
        ns = CosmoDict["ns"],
        σ8 = CosmoDict["σ8"])
    else
        Cosmology = Flatw0waCDMCosmology(
        w0 = CosmoDict["w0"],
        wa = CosmoDict["wa"],
        Mν = CosmoDict["Mν"],
        H0 = CosmoDict["H0"],
        ΩM = CosmoDict["ΩM"],
        ΩB = CosmoDict["ΩB"],
        ns = CosmoDict["ns"],
        σ8 = CosmoDict["σ8"])
    end
    return Cosmology
end

function ReadIntrinsicAlignment(CosmoDict::Dict)
    intrinsicalignment = ExtendedNLIA()
    intrinsicalignment.𝓐IA = CosmoDict["𝓐IA"]
    intrinsicalignment.βIA = CosmoDict["βIA"]
    intrinsicalignment.𝓒IA = CosmoDict["𝓒IA"]
    intrinsicalignment.ηIA = CosmoDict["ηIA"]
    return intrinsicalignment
end

function ReadIntrinsicAlignmentForecast(IADict::Dict, IAModel::String)
    if IAModel == "ExtendedNLIA"
        intrinsicalignment = ExtendedNLIA()
        intrinsicalignment.𝓐IA = IADict["𝓐IA"][1]
        intrinsicalignment.βIA = IADict["βIA"][1]
        intrinsicalignment.𝓒IA = IADict["𝓒IA"][1]
        intrinsicalignment.ηIA = IADict["ηIA"][1]
    elseif IAModel == "None"
        intrinsicalignment = AbsentIA()
    else
        ErrorException("Intrinsc Alignment model not defined correctly.")
    end
    return intrinsicalignment
end

function ReadBias(CosmoDict::Dict)
    bias = EuclidBias()
    bias.A = CosmoDict["A"]
    bias.B = CosmoDict["B"]
    bias.C = CosmoDict["C"]
    bias.D = CosmoDict["D"]
    return bias
end

function ReadBiasForecast(BiasDict::Dict, BiasModel::String)
    if BiasModel == "EuclidBias"
        bias = EuclidBias()
        bias.A = BiasDict["A"][1]
        bias.B = BiasDict["B"][1]
        bias.C = BiasDict["C"][1]
        bias.D = BiasDict["D"][1]
    elseif BiasModel == "PiecewiseBias"
        bias = PiecewiseBias()
        bias.BiasMultiplier = BiasDict["BiasMultiplier"]
    else
        ErrorException("Bias model not defined correctly.")
    end
    return bias
end

function ReadCℓ(Filename::String, Probes::String)
    Filename *= ".h5"
    file = HDF5.h5open(Filename, "r")
    c_lij = HDF5.read(file["cls"][Probes]["c_lij"])
    return Cℓ(CℓArray = c_lij)
end

function Read∂Cℓ(Filename::String, Probe::String)
    Filename *= ".h5"
    file = HDF5.h5open(Filename, "r")
    dc_lij = HDF5.read(file["dcls"][Probe]["dc_lij"])
    return ∂Cℓ(∂CℓArray = dc_lij)
end

function Write∂Cℓ!(DerivativeArray::AbstractArray{Float64, 3},
    Filename::String, Probes::String)
    h5write(Filename*".h5", "dcls/"*Probes*"/dc_lij",
    DerivativeArray)
end

"""
    WritePowerSpectrumBackground(PowerSpectrum::PowerSpectrum,
    BackgroundQuantities::BackgroundQuantities,
    CosmologicalGrid::CosmologicalGrid, Filename::String)

This function writes the Power Spectrum, the Background quantities and the
Cosmological Grid in a HDF5 file.
"""
function WritePowerSpectrumBackground(PowerSpectrum::AbstractPowerSpectrum,
    BackgroundQuantities::AbstractBackgroundQuantities,
    CosmologicalGrid::AbstractCosmologicalGrid, Filename::String)
    h5write(Filename*".h5",
    "cosmology/comoving_distance_array",
    BackgroundQuantities.rZArray)
    h5write(Filename*".h5",
    "cosmology/hubble_array",
    BackgroundQuantities.HZArray)
    h5write(Filename*".h5",
    "cosmology/z_grid",
    CosmologicalGrid.ZArray)
    h5write(Filename*".h5",
    "power_spectrum/z_grid",
    CosmologicalGrid.ZArray)
    h5write(Filename*".h5",
    "power_spectrum/k_grid",
    CosmologicalGrid.KArray)
    h5write(Filename*".h5",
    "power_spectrum/lin_p_mm_k_z",
    PowerSpectrum.PowerSpectrumLinArray)
    h5write(Filename*".h5",
    "power_spectrum/nonlin_p_mm_k_z",
    PowerSpectrum.PowerSpectrumNonlinArray)
end

"""
    ReadPowerSpectrumBackground(Filename::String,
    ℓBinCenters::Vector{Float64})

This function reads the Power Spectrum, the Background quantities and the
Cosmological Grid from a HDF5 file.
"""
function ReadPowerSpectrumBackground(Filename::String, ℓBinCenters::Vector{Float64},
    ℓBinWidths::Vector{Float64})
    Filename *= ".h5"
    file = HDF5.h5open(Filename, "r")
    nonlin_p_mm_k_z = HDF5.read(file["power_spectrum"]["nonlin_p_mm_k_z"])
    lin_p_mm_k_z    = HDF5.read(file["power_spectrum"]["lin_p_mm_k_z"])
    z_grid = HDF5.read(file["power_spectrum"]["z_grid"])
    k_grid = HDF5.read(file["power_spectrum"]["k_grid"])
    comoving_distance_array = HDF5.read(file["cosmology"]["comoving_distance_array"])
    hubble_array = HDF5.read(file["cosmology"]["hubble_array"])
    BackgroundQuantitiesRead = BackgroundQuantities(HZArray = hubble_array,
    rZArray = comoving_distance_array)
    CosmologicalGridRead = CosmologicalGrid(ZArray = z_grid, KArray = k_grid,
    ℓBinCenters = ℓBinCenters, ℓBinWidths = ℓBinWidths)
    PowerSpectrumRead = PowerSpectrum(PowerSpectrumLinArray = lin_p_mm_k_z,
    PowerSpectrumNonlinArray = nonlin_p_mm_k_z,
    InterpolatedPowerSpectrum = zeros(length(CosmologicalGridRead.ℓBinCenters),
    length(CosmologicalGridRead.ZArray)))
    return PowerSpectrumRead, BackgroundQuantitiesRead, CosmologicalGridRead
end

function ReadPowerSpectrumBackgroundSeyfert(Filename::String,
    ℓBinCenters::Vector{Float64})
    c_0 = 2.99792458e5
    Filename *= ".h5"
    file = HDF5.h5open(Filename, "r")
    nonlin_p_mm_k_z = HDF5.read(file["power_spectrum"]["nonlin_p_mm_z_k"])
    lin_p_mm_k_z    = HDF5.read(file["power_spectrum"]["lin_p_mm_z_k"])
    z_grid = HDF5.read(file["power_spectrum"]["z_grid"])
    k_grid = HDF5.read(file["power_spectrum"]["k_grid"])
    dimensionless_comoving_distance_array = HDF5.read(file["cosmology"]["dimensionless_comoving_distance_array"])
    dimensionless_hubble_array = HDF5.read(file["cosmology"]["dimensionless_hubble_array"])
    BackgroundQuantities = BackgroundQuantities(HZArray =
    dimensionless_hubble_array*67,
    rZArray = dimensionless_comoving_distance_array*c_0/67)
    CosmologicalGrid = CosmologicalGrid(ZArray = z_grid, KArray = k_grid,
    ℓBinCenters = ℓBinCenters)
    PowerSpectrum = PowerSpectrum(PowerSpectrumLinArray = lin_p_mm_k_z,
    PowerSpectrumNonlinArray = nonlin_p_mm_k_z,
    InterpolatedPowerSpectrum = zeros(length(CosmologicalGrid.ℓBinCenters),
    length(CosmologicalGrid.ZArray)))
    return PowerSpectrum, BackgroundQuantities, CosmologicalGrid
end

function WriteCosmologicalGrid!(Filename::String, cosmogrid::CosmologicalGrid)
    h5write(Filename*".h5", "cosmologicalgrid/l_bin_centers", cosmogrid.ℓBinCenters)
    h5write(Filename*".h5", "cosmologicalgrid/l_bin_widths", cosmogrid.ℓBinWidths)
end

function WriteFisher!(Fisher::CosmoCentral.Fisherαβ, Filename::String, Probes::String)
    h5write(Filename*".h5", "FisherMatrices/"*Probes*"/FisherMatrix",Fisher.FisherMatrix)
    h5write(Filename*".h5", "FisherMatrices/"*Probes*"/FisherParameters",Fisher.SelectedParametersList)
end

function ReadFisher(Filename::String, Probes::String)
    Filename *= ".h5"
    file = HDF5.h5open(Filename, "r")
    fishermatrix = HDF5.read(file["FisherMatrices"][Probes]["FisherMatrix"])
    parameters = HDF5.read(file["FisherMatrices"][Probes]["FisherParameters"])
    return Fisherαβ(FisherMatrix=fishermatrix, SelectedParametersList=parameters, ParametersList = parameters)
end


text:
function HatanakaCompressFilter()
    return HatanakaCompressFilter(())
end

function filter(obj::HatanakaCompressFilter, arg0::NamedData)
    return jcall(obj, "filter", NamedData, (NamedData,), arg0)
end


text:
using Distributions
using IncrementalInference

import Unmarshal: unmarshal
import Base: convert

# Already in IIF DefaultNodeTypes.jl Ln151
# """
# Converter: Prior -> PackedPrior::Dict{String, Any}
# """
# function convert(::Type{Dict{String, Any}}, prior::IncrementalInference.Prior)
#     z = convert(Type{Dict{String, Any}}, prior.Z)
#     return Packed_Factor([z], "Prior")
# end

# """
# Converter: PackedPrior::Dict{String, Any} -> Prior
# """
# function convert(::Type{Prior}, prior::Dict{String, Any})
#     # Genericize to any packed type next.
#     z = prior["measurement"][1]
#     z = convert(_evalType(z["distType"]), z)
#     return Prior(z)
# end

"""
Converter: PackedPrior::Dict{String, Any} -> PriorPose2
"""
function convert(::Type{PriorPose2}, prior::Dict{String, Any})
    # Genericize to any packed type next.
    z = prior["measurement"][1]
    z = convert(_evalType(z["distType"]), z)
    return PriorPose2(z)
end

"""
Converter: PackedPrior::Dict{String, Any} -> PriorPose3
"""
function convert(::Type{PriorPose3}, prior::Dict{String, Any})
    # Genericize to any packed type next.
    z = prior["measurement"][1]
    z = convert(_evalType(z["distType"]), z)
    return PriorPose3(z)
end

"""
Converter: Pose2Pose2::Dict{String, Any} -> Pose2Pose2
"""
function convert(::Type{Pose2Pose2}, prior::Dict{String, Any})
    # Genericize to any packed type next.
    z = prior["measurement"][1]
    z = convert(_evalType(z["distType"]), z)
    return Pose2Pose2(z)
end

"""
Converter: Pose2Point2BearingRange::Dict{String, Any} -> Pose2Point2BearingRange
"""
function convert(::Type{Pose2Point2BearingRange}, fact::Dict{String, Any})
    # Genericize to any packed type next.
    b = fact["measurement"][1]
    r = fact["measurement"][2]
    b = convert(_evalType(b["distType"]), b)
    r = convert(_evalType(r["distType"]), r)
    return Pose2Point2BearingRange(b, r)
end

"""
Converter: Pose3Pose3::Dict{String, Any} -> Pose3Pose3
"""
function convert(::Type{Pose3Pose3}, fact::Dict{String, Any})
    # Genericize to any packed type next.
    z = fact["measurement"][1]
    zdist = convert(_evalType(z["distType"]), z)
    return Pose3Pose3(zdist)
end

# already added in RoME PartialPose3.jl Ln230
# """
# Converter: PartialPose3XYYaw::Dict{String, Any} -> PartialPose3XYYaw
# """
# function convert(::Type{PartialPose3XYYaw}, fact::Dict{String, Any})
#     # Genericize to any packed type next.
#     xy = fact["measurement"][1]
#     yaw = fact["measurement"][2]
#     xydistr = convert(_evalType(xy["distType"]), xy)
#     yawdistr = convert(_evalType(yaw["distType"]), yaw)
#     return PartialPose3XYYaw(xydistr, yawdistr)
# end

# already added in RoME PartialPose3.jl Ln89
# """
# Converter: PartialPriorRollPitchZ::Dict{String, Any} -> PartialPriorRollPitchZ
# """
# function convert(::Type{PartialPriorRollPitchZ}, fact::Dict{String, Any})
#     # Genericize to any packed type next.
#     z = fact["measurement"][1]
#     rp = fact["measurement"][2]
#     rpdistr = convert(_evalType(rp["distType"]), rp)
#     zdistr = convert(_evalType(z["distType"]), z)
#     return PartialPriorRollPitchZ(rpdistr, zdistr)
# end

# TODO: Write the other side of the converters for above before moving.


text:
#! /Applications/Julia-0.4.5.app/Contents/Resources/julia/bin/julia
# Geoffrey Hannigan
# University of Michigan

# Setting dependencies
using Bio.Align
using Bio.Seq
using ArgParse

# Setup ArgParse
argx = ArgParseSettings()
@add_arg_table argx begin
    "--input", "-i"
        help = "Input fasta file for identifying circular contigs."
    "--output", "-o"
        help = "Output file to write reults."
    "--length", "-l"
        help = "Length of sequence window for detecting repeats."
        default = 100
    "--thresh", "-t"
        help = "Threshold (percent of length) for circular contig."
        default = 0.97
end

parsed_args = parse_args(ARGS, argx)

problem = LocalAlignment()

scoremodel = AffineGapScoreModel(
    match=5,
    mismatch=-4,
    gap_open=-4,
    gap_extend=-1
)

filepath = parsed_args["input"]
fileout = parsed_args["output"]
outputfile = open(fileout, "w")
hitlength = parsed_args["length"]
threshold = hitlength * 5 * parsed_args["thresh"]
maxscore = hitlength * 5
print("Perfect alignment score is $maxscore\n")
inputfasta = FASTA.Reader(open( filepath))

for s in inputfasta
    alignment_result = pairalign(problem, FASTA.sequence(s , 1:hitlength), FASTA.sequence(s , hitlength:length(FASTA.sequence(s))), scoremodel)
    scorestring = score(alignment_result)
    seqname = FASTA.identifier(s)
    if scorestring > threshold
        write(outputfile, "$seqname\t$scorestring\n")
    end
end

print("Complete\n")

max_stars_repo_name: Microbiology/ccontigs

max_issues_repo_name: Microbiology/ccontigs

max_forks_repo_name: Microbiology/ccontigs


text:
## standalone script to extract recursively the stellar clusters in a votable with or w/o
## optimization of the DBSCAN parameters.
##

using DataFrames, Query
using CSV, Glob, Dates
using Statistics, Random, UUIDs
using Printf, ArgParse

rootdir =  ENV["GAIA_ROOT"]

push!(LOAD_PATH,"$rootdir/run/src")
using GaiaClustering

##
## parse extra options/flags
##
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "-m"
            help = "configuration file"
            arg_type = String
        "--ncycle", "-n"
            help = "maximum number of cycles"
            arg_type = Int
        "-o"
            help = "optimization of the weightings/DBSCAN"
            action = :store_true
        "--maxdist" , "-d"
            help = "maximum distance for stars in pc"
            arg_type = Float64
        "--mindist"
            help = "miniimum distance for stars in pc"
            arg_type = Float64
        "--qmetric", "-q"
            help = "Q metric method to choose best solution in final DBSCAN clusters. The options are: Qc, Qn, QcQn, QcQnhigh"
            arg_type = String
        "--w3d"
            help = "XYZ weighting"
            arg_type = Float64
        "--wvel"
            help = "Velocity weighting"
            arg_type = Float64
        "--whrd"
            help = "Magnitud/color weighting"
            arg_type = Float64
        "--eps"
            help = "Epsilon parameter of DBSCAN"
            arg_type = Float64
        "--mcl"
            help = "Min_cluster parameter of DBSCAN"
            arg_type = Int
        "--mnei"
            help = "Min_neighbor parameter of DBSCAN"
            arg_type = Int
        "votable"
            help = "votable"
            required = true
    end

    return parse_args(s)
end
##
## get the data. The weightings are fake. Should be applied later
##
function getdata(m::GaiaClustering.meta)
    println("## Distance cut : $(m.mindist) $(m.maxdist) pc")

    data       = read_votable(m.votname)
    df         = filter_data(data, [m.mindist, m.maxdist])
    dfcart     = add_cartesian(df)
    blck       = [[1,2,3],[4,5], [6,7,8]]
    wghtblck   = [4.0,5.0,1.0]
    norm       = "identity"

    dfcartnorm , scale8 = normalization_PerBlock(dfcart, blck, wghtblck , norm, false)

    return(df, dfcart , dfcartnorm)
end
##
## Main function
##
function main(m::GaiaClustering.meta, optim)
    tstart= now()
    rng = MersenneTwister()
    uuid=uuid4(rng)
    m.uuid= uuid

    println("###########################")
    println("## Starting with $(m.votname)")
    println("## Starting at $tstart")
    println("## Id $uuid")

    df , dfcart , dfcartnorm = getdata(m)

    cycle, flag= cycle_extraction_optim(df, dfcart, m, optim)

    tend= now()
    println("## Ending at $tend")
    println("## number of cycle: $cycle , flag:$flag ")
    println("##########################")
    println("##")

    duration= Dates.value(tend-tstart) / (1000*3600)
    durationstr= @sprintf("%3.3f", duration)
    @printf("## %s \n",specialstr("Duration: $durationstr hours","YELLOW"))
    @printf("## %s \n",specialstr("Votable done: $(m.votname)","YELLOW"))
    println("##\n##")
end
##################################### MAIN
let
    parsed_args = parse_commandline()

    metafile= parsed_args["m"]
    votable= parsed_args["votable"]
    ncycle= parsed_args["ncycle"]
    qmetric= parsed_args["qmetric"]
    maxdist= parsed_args["maxdist"]
    mindist= parsed_args["mindist"]
    eps= parsed_args["eps"]
    mcl= parsed_args["mcl"]
    mnei= parsed_args["mnei"]
    w3d= parsed_args["w3d"]
    wvel= parsed_args["wvel"]
    whrd= parsed_args["whrd"]

    if parsed_args["o"] opt="yes" else opt= "no" end

    ############
    header_extract()

    if metafile != nothing
        m= read_params(metafile, false)
        opt= m.optim
    else
        println("## The default options are used.")
        m= set_default_params()
    end

    m.optim= opt
    m.votname= votable

    if m.optim == "yes"
        isoptimize= true
    else
        isoptimize= false
    end

    if qmetric != nothing m.labels= qmetric end
    if ncycle != nothing m.cyclemax= ncycle end
    if maxdist != nothing m.maxdist= maxdist end
    if mindist != nothing m.mindist= mindist end
    if eps != nothing m.eps= eps end
    if mcl != nothing m.mcl= mcl end
    if mnei != nothing m.mnei= mnei end
    if w3d != nothing m.w3d= w3d end
    if wvel != nothing m.wvel= wvel end
    if whrd != nothing m.whrd= whrd end

    main(m, isoptimize)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
# Simulate user inputs
function eventbutton(c, event_type, btn, x=DeviceUnit(0), y=DeviceUnit(0), state=0)
    xd, yd = GtkObservables.convertunits(DeviceUnit, c, x, y)
    Gtk.GdkEventButton(event_type,
                       Gtk.gdk_window(widget(c)),
                       Int8(0),
                       UInt32(0),
                       convert(Float64, xd), convert(Float64, yd),
                       convert(Ptr{Float64},C_NULL),
                       UInt32(state),
                       UInt32(btn),
                       C_NULL,
                       0.0, 0.0)
end
function eventscroll(c, direction, x=DeviceUnit(0), y=DeviceUnit(0), state=0)
    xd, yd = GtkObservables.convertunits(DeviceUnit, c, x, y)
    Gtk.GdkEventScroll(Gtk.GdkEventType.SCROLL,
                       Gtk.gdk_window(widget(c)),
                       Int8(0),
                       UInt32(0),
                       convert(Float64, xd), convert(Float64, yd),
                       UInt32(state),
                       direction,
                       convert(Ptr{Float64},C_NULL),
                       0.0, 0.0,
                       0.0, 0.0)
end
function eventmotion(c, btn, x, y)
    xd, yd = GtkObservables.convertunits(DeviceUnit, c, x, y)
    Gtk.GdkEventMotion(Gtk.GdkEventType.MOTION_NOTIFY,
                       Gtk.gdk_window(widget(c)),
                       Int8(0),
                       UInt32(0),
                       convert(Float64, xd), convert(Float64, yd),
                       convert(Ptr{Float64},C_NULL),
                       UInt32(btn),
                       Int16(0),
                       C_NULL,
                       0.0, 0.0)
end

const ModType = Gtk.GConstants.GdkModifierType
mask(btn) =
    btn == 1 ? ModType.GDK_BUTTON1_MASK :
    btn == 2 ? ModType.GDK_BUTTON2_MASK :
    btn == 3 ? ModType.GDK_BUTTON3_MASK :
    error(btn, " not recognized")


text:
########################### Thermal Generation Models ######################################

abstract type AbstractThermalFormulation <: AbstractDeviceFormulation end

abstract type AbstractThermalDispatchFormulation <: AbstractThermalFormulation end

abstract type AbstractThermalUnitCommitment <: AbstractThermalFormulation end

struct ThermalBasicUnitCommitment <: AbstractThermalUnitCommitment end

struct ThermalStandardUnitCommitment <: AbstractThermalUnitCommitment end

struct ThermalDispatch <: AbstractThermalDispatchFormulation end

struct ThermalRampLimited <: AbstractThermalDispatchFormulation end

struct ThermalDispatchNoMin <: AbstractThermalDispatchFormulation end

########################### Active Dispatch Variables ######################################

"""
This function add the variables for power generation output to the model
"""
function activepower_variables!(canonical_model::CanonicalModel,
                           devices::IS.FlattenIteratorWrapper{T}) where {T<:PSY.ThermalGen}

    add_variable(canonical_model,
                 devices,
                 Symbol("P_$(T)"),
                 false,
                 :nodal_balance_active;
                 ub_value = d -> d.tech.activepowerlimits.max,
                 lb_value = d -> d.tech.activepowerlimits.min,
                 init_value = d -> PSY.get_activepower(PSY.get_tech(d)))

    return

end

"""
This function add the variables for power generation output to the model
"""
function reactivepower_variables!(canonical_model::CanonicalModel,
                           devices::IS.FlattenIteratorWrapper{T}) where {T<:PSY.ThermalGen}

    add_variable(canonical_model,
                 devices,
                 Symbol("Q_$(T)"),
                 false,
                 :nodal_balance_reactive;
                 ub_value = d -> d.tech.reactivepowerlimits.max,
                 lb_value = d -> d.tech.reactivepowerlimits.min,
                 init_value = d -> d.tech.reactivepower)

    return

end

"""
This function add the variables for power generation commitment to the model
"""
function commitment_variables!(canonical_model::CanonicalModel,
                           devices::IS.FlattenIteratorWrapper{T}) where {T<:PSY.ThermalGen}

    time_steps = model_time_steps(canonical_model)
    var_names = [Symbol("ON_$(T)"), Symbol("START_$(T)"), Symbol("STOP_$(T)")]

    for v in var_names
        add_variable(canonical_model, devices, v, true)
    end

    return

end

"""
This function adds the active power limits of generators when there are no CommitmentVariables
"""
function activepower_constraints!(canonical_model::CanonicalModel,
                                 devices::IS.FlattenIteratorWrapper{T},
                                 device_formulation::Type{D},
                                 system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                                     D<:AbstractThermalDispatchFormulation,
                                                                     S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(g),  PSY.get_activepowerlimits(PSY.get_tech(g))) for g in devices]

    device_range(canonical_model,
                 range_data,
                 Symbol("activerange_$(T)"),
                 Symbol("P_$(T)"))
    return

end

"""
This function adds the active power limits of generators when there are CommitmentVariables
"""
function activepower_constraints!(canonical_model::CanonicalModel,
                                 devices::IS.FlattenIteratorWrapper{T},
                                 device_formulation::Type{D},
                                 system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                                      D<:AbstractThermalFormulation,
                                                                      S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(g),  PSY.get_activepowerlimits(PSY.get_tech(g))) for g in devices]
    device_semicontinuousrange(canonical_model,
                               range_data,
                               Symbol("activerange_$(T)"),
                               Symbol("P_$(T)"),
                               Symbol("ON_$(T)"))

    return

end

"""
This function adds the active power limits of generators when there are
    no CommitmentVariables
"""
function activepower_constraints!(canonical_model::CanonicalModel,
                                  devices::IS.FlattenIteratorWrapper{T},
                                  device_formulation::Type{ThermalDispatchNoMin},
                                  system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                                     S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(g), (min = 0.0, max=( PSY.get_activepowerlimits(PSY.get_tech(g))).max)) for g in devices]
    var_key = Symbol("P_$(T)")
    variable = var(canonical_model, var_key)

    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            JuMP.set_lower_bound(v, 0.0)
        end
    end

    device_range(canonical_model,
                range_data,
                Symbol("activerange_$(T)"),
                Symbol("P_$(T)")
                )

    return

end

"""
This function adds the reactive  power limits of generators when there are CommitmentVariables
"""
function reactivepower_constraints!(canonical_model::CanonicalModel,
                                   devices::IS.FlattenIteratorWrapper{T},
                                   device_formulation::Type{D},
                                   system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                                       D<:AbstractThermalDispatchFormulation,
                                                                       S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(g),  PSY.get_reactivepowerlimits(PSY.get_tech(g))) for g in devices]

    device_range(canonical_model,
                 range_data ,
                 Symbol("reactiverange_$(T)"),
                 Symbol("Q_$(T)"))

    return

end

"""
This function adds the reactive power limits of generators when there CommitmentVariables
"""
function reactivepower_constraints!(canonical_model::CanonicalModel,
                                   devices::IS.FlattenIteratorWrapper{T},
                                   device_formulation::Type{D},
                                   system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                                        D<:AbstractThermalFormulation,
                                                                        S<:PM.AbstractPowerModel}

    range_data = [(PSY.get_name(g),  PSY.get_reactivepowerlimits(PSY.get_tech(g))) for g in devices]

    device_semicontinuousrange(canonical_model,
                               range_data,
                               Symbol("reactiverange_$(T)"),
                               Symbol("Q_$(T)"),
                               Symbol("ON_$(T)"))

    return

end

### Constraints for Thermal Generation without commitment variables ####
"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitment_constraints!(canonical_model::CanonicalModel,
                                 devices::IS.FlattenIteratorWrapper{T},
                                 device_formulation::Type{D},
                                 system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                                     D<:AbstractThermalFormulation,
                                                                     S<:PM.AbstractPowerModel}

    key = ICKey(DeviceStatus, T)

    if !(key in keys(canonical_model.initial_conditions))
        error("Initial status conditions not provided. This can lead to unwanted results")
    end

    device_commitment(canonical_model,
                     canonical_model.initial_conditions[key],
                     Symbol("commitment_$(T)"),
                     (Symbol("START_$(T)"),
                      Symbol("STOP_$(T)"),
                      Symbol("ON_$(T)"))
                      )

    return

end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(canonical_model::CanonicalModel,
                            devices::IS.FlattenIteratorWrapper{T},
                            device_formulation::Type{D}) where {T<:PSY.ThermalGen,
                                                                D<:AbstractThermalFormulation}

    status_init(canonical_model, devices)
    output_init(canonical_model, devices)
    duration_init(canonical_model, devices)

    return

end

function initial_conditions!(canonical_model::CanonicalModel,
                            devices::IS.FlattenIteratorWrapper{T},
                            device_formulation::Type{D}) where {T<:PSY.ThermalGen,
                                                                D<:AbstractThermalDispatchFormulation}

    output_init(canonical_model, devices)

    return

end

########################### Ramp/Rate of Change constraints ################################
"""
This function gets the data for the generators
"""
function _get_data_for_rocc(initial_conditions::Vector{InitialCondition},
                            resolution::Dates.Period)

    if resolution > Dates.Minute(1)
        minutes_per_period = Dates.value(Dates.Minute(resolution))
    else
        minutes_per_period = Dates.value(Dates.Second(resolution))/60
    end

    lenght_devices = length(initial_conditions)
    ini_conds = Vector{InitialCondition}(undef, lenght_devices)
    ramp_params = Vector{UpDown}(undef, lenght_devices)
    minmax_params = Vector{MinMax}(undef, lenght_devices)

    idx = 0
    for ic in initial_conditions
        g = ic.device
        gen_tech = PSY.get_tech(g)
        name = PSY.get_name(g)
        non_binding_up = false
        non_binding_down = false
        ramplimits =  PSY.get_ramplimits(gen_tech)
        rating = PSY.get_rating(gen_tech)
        if !isnothing(ramplimits)
            p_lims = PSY.get_activepowerlimits(gen_tech)
            max_rate = abs(p_lims.min - p_lims.max)/minutes_per_period
            if (ramplimits.up*rating >= max_rate) & (ramplimits.down*rating >= max_rate)
                @info "Generator $(name) has a nonbinding ramp limits. Constraints Skipped"
                continue
            else
                idx += 1
            end
            ini_conds[idx] = ic
            ramp_params[idx] = (up = ramplimits.up*rating*minutes_per_period,
                                down = ramplimits.down*rating*minutes_per_period)
            minmax_params[idx] = p_lims
        end
    end

    if idx < lenght_devices
        deleteat!(ini_conds, idx+1:lenght_devices)
        deleteat!(ramp_params, idx+1:lenght_devices)
        deleteat!(minmax_params, idx+1:lenght_devices)
    end

    return ini_conds, ramp_params, minmax_params

end

"""
This function adds the ramping limits of generators when there are CommitmentVariables
"""
function ramp_constraints!(canonical_model::CanonicalModel,
                           devices::IS.FlattenIteratorWrapper{T},
                           device_formulation::Type{D},
                           system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                    D<:AbstractThermalFormulation,
                                                    S<:PM.AbstractPowerModel}
    key = ICKey(DevicePower, T)

    if !(key in keys(canonical_model.initial_conditions))
        error("Initial Conditions for $(T) Rate of Change Constraints not in the model")
    end

    time_steps = model_time_steps(canonical_model)
    resolution = model_resolution(canonical_model)
    initial_conditions = get_ini_cond(canonical_model, key)
    rate_data = _get_data_for_rocc(initial_conditions, resolution)
    ini_conds, ramp_params, minmax_params = _get_data_for_rocc(initial_conditions, resolution)

    if !isempty(ini_conds)
        # Here goes the reactive power ramp limits when versions for AC and DC are added
        device_mixedinteger_rateofchange(canonical_model,
                                         (ramp_params, minmax_params),
                                         ini_conds,
                                         Symbol("ramp_$(T)"),
                                        (Symbol("P_$(T)"),
                                         Symbol("START_$(T)"),
                                         Symbol("STOP_$(T)"))
                                        )
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end

    return

end

function ramp_constraints!(canonical_model::CanonicalModel,
                          devices::IS.FlattenIteratorWrapper{T},
                          device_formulation::Type{D},
                          system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                   D<:AbstractThermalDispatchFormulation,
                                                   S<:PM.AbstractPowerModel}

    key = ICKey(DevicePower, T)

    if !(key in keys(canonical_model.initial_conditions))
        error("Initial Conditions for $(T) Rate of Change Constraints not in the model")
    end

    time_steps = model_time_steps(canonical_model)
    resolution = model_resolution(canonical_model)
    initial_conditions = get_ini_cond(canonical_model, key)
    rate_data = _get_data_for_rocc(initial_conditions, resolution)
    ini_conds, ramp_params, minmax_params = _get_data_for_rocc(initial_conditions, resolution)

    if !isempty(ini_conds)
        # Here goes the reactive power ramp limits when versions for AC and DC are added
        device_linear_rateofchange(canonical_model,
                                  ramp_params,
                                  ini_conds,
                                   Symbol("ramp_$(T)"),
                                   Symbol("P_$(T)"))
    else
        @warn "Data doesn't contain generators with ramp limits, consider adjusting your formulation"
    end

    return

end

########################### time duration constraints ######################################
"""
If the fraction of hours that a generator has a duration constraint is less than
the fraction of hours that a single time_step represents then it is not binding.
"""
function _get_data_for_tdc(initial_conditions_on::Vector{InitialCondition},
                           initial_conditions_off::Vector{InitialCondition},
                           resolution::Dates.Period)

    steps_per_hour = 60/Dates.value(Dates.Minute(resolution))
    fraction_of_hour = 1/steps_per_hour
    lenght_devices_on = length(initial_conditions_on)
    lenght_devices_off = length(initial_conditions_off)
    @assert lenght_devices_off == lenght_devices_on
    time_params = Vector{UpDown}(undef, lenght_devices_on)
    ini_conds = Matrix{InitialCondition}(undef, lenght_devices_on, 2)

    idx = 0
    for (ix, ic) in enumerate(initial_conditions_on)
        g = ic.device
        @assert g == initial_conditions_off[ix].device
        tech = PSY.get_tech(g)
        non_binding_up = false
        non_binding_down = false
        timelimits =  PSY.get_timelimits(tech)
        name = PSY.get_name(g)
        if !isnothing(timelimits)
            if (timelimits.up <= fraction_of_hour) & (timelimits.down <= fraction_of_hour)
                @info "Generator $(name) has a nonbinding time limits. Constraints Skipped"
            else
                idx += 1
            end
            ini_conds[idx, 1] = ic
            ini_conds[idx, 2] = initial_conditions_off[ix]
            up_val = round(timelimits.up * steps_per_hour, RoundUp)
            down_val = round(timelimits.down * steps_per_hour, RoundUp)
            time_params[idx] = time_params[idx] = (up = up_val, down = down_val)
        end
    end

    if idx < lenght_devices_on
        ini_conds = ini_conds[1:idx,:]
        deleteat!(time_params, idx+1:lenght_devices_on)
    end

    return ini_conds, time_params

end

function time_constraints!(canonical_model::CanonicalModel,
                          devices::IS.FlattenIteratorWrapper{T},
                          device_formulation::Type{D},
                          system_formulation::Type{S}) where {T<:PSY.ThermalGen,
                                                   D<:AbstractThermalFormulation,
                                                   S<:PM.AbstractPowerModel}

    ic_keys = [ICKey(TimeDurationON, T), ICKey(TimeDurationOFF, T)]
    for key in ic_keys
        if !(key in keys(canonical_model.initial_conditions))
            error("Initial Conditions for $(T) Time Constraint not in the model")
        end
    end

    parameters = model_has_parameters(canonical_model)
    resolution = model_resolution(canonical_model)
    initial_conditions_on  = get_ini_cond(canonical_model, ic_keys[1])
    initial_conditions_off = get_ini_cond(canonical_model, ic_keys[2])
    ini_conds, time_params = _get_data_for_tdc(initial_conditions_on,
                                               initial_conditions_off,
                                               resolution)

    if !(isempty(ini_conds))
       if parameters
            device_duration_param(canonical_model,
                                time_params,
                                ini_conds,
                                Symbol("duration_$(T)"),
                                (Symbol("ON_$(T)"),
                                Symbol("START_$(T)"),
                                Symbol("STOP_$(T)"))
                                      )
        else
            device_duration_retrospective(canonical_model,
                                        time_params,
                                        ini_conds,
                                        Symbol("duration_$(T)"),
                                        (Symbol("ON_$(T)"),
                                        Symbol("START_$(T)"),
                                        Symbol("STOP_$(T)"))
                                        )
        end
    else
        @warn "Data doesn't contain generators with time-up/down limits, consider adjusting your formulation"
    end

    return

end

########################### Cost Function Calls#############################################

function cost_function(canonical_model::CanonicalModel,
                       devices::IS.FlattenIteratorWrapper{T},
                       ::Type{D},
                       ::Type{S}) where {T<:PSY.ThermalGen,
                                         D<:AbstractThermalDispatchFormulation,
                                         S<:PM.AbstractPowerModel}

    add_to_cost(canonical_model,
                devices,
                Symbol("P_$(T)"),
                :variable)

    return

end

function cost_function(canonical_model::CanonicalModel,
                       devices::IS.FlattenIteratorWrapper{T},
                       ::Type{D},
                       ::Type{S}) where {T<:PSY.ThermalGen,
                                         D<:AbstractThermalFormulation,
                                         S<:PM.AbstractPowerModel}

    #Variable Cost component
    add_to_cost(canonical_model, devices, Symbol("P_$(T)"), :variable)

    #Commitment Cost Components
    add_to_cost(canonical_model, devices, Symbol("START_$(T)"), :startup)
    add_to_cost(canonical_model, devices, Symbol("STOP_$(T)"), :shutdn)
    add_to_cost(canonical_model, devices, Symbol("ON_$(T)"), :fixed)

    return

end


text:
using Documenter
using SatelliteToolbox
using SatelliteToolbox.SGP4
using SatelliteToolbox.SatelliteToolboxTLE

makedocs(
    modules = [SatelliteToolbox,
               SatelliteToolbox.SGP4,
               SatelliteToolbox.SatelliteToolboxTLE],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://juliaspace.github.io/SatelliteToolbox.jl/stable/",
    ),
    sitename = "Satellite Toolbox",
    authors = "Ronan Arraes Jardim Chagas",
    pages = [
        "Home" => "index.md",
        "Earth" => Any[
            "Earth atmospheric models" => "man/earth/atmospheric_models.md",
            "Earth geomagnetic field models" => "man/earth/geomagnetic_field_models.md",
            "Space indices" => "man/earth/space_indices.md",
           ],
        "Orbit" => Any[
            "Anomalies" => "man/orbit/anomalies.md",
            "General analysis" => "man/orbit/general.md",
            "Orbit propagators" => "man/orbit/propagators.md",
            "TLE" => "man/orbit/tle.md",
        ],
        "Transformations" => Any[
            "ECEF and ECI" => "man/transformations/ecef_eci.md",
            "Geodetic and Geocentric" => "man/transformations/geodetic_geocentric.md",
        ],
        "Library" => "lib/library.md",
    ],
)

deploydocs(
    repo = "github.com/JuliaSpace/SatelliteToolbox.jl.git",
    target = "build",
)


text:
using Test

@testset "MultiThreadedCaches.jl" begin
    include("MultiThreadedCaches.jl")
end


text:
using Random;
using LinearAlgebra;
using Statistics;
using ImageView, Images, Gtk.ShortNames;

#===============================================#
#                                               #
#                Tools functions                #
#                                               #
#===============================================#
function load_mnist(path, kind)
    labels_path = joinpath(path, "$(kind)-labels-idx1-ubyte")
    images_path = joinpath(path, "$(kind)-images-idx3-ubyte")

    labels = []
    open(labels_path, "r") do labels_file
        num, magic = reinterpret(Int32, reverse(read(labels_file,8)))
        labels = read(labels_file, num)
    end

    images = []
    open(images_path, "r") do images_file
        cols, rows, num, magic = reinterpret(Int32, reverse(read(images_file,16)))
        images = read(images_file, num * rows * cols)
        images = reshape(images, (rows * cols, num))
        # images->reshape(num, rows, cols)
    end

    return labels, images
end

#===============================================#
#                                               #
#                      PCA                      #
#                                               #
#===============================================#
function PCA(Data, Top_n_feat=10)
    # 特征标准化，平衡各个特征尺度
    μ = mean(Data, dims=2)
    mean_removed = Data .- μ

    # 计算协方差矩阵（Find covariance matrix）
    Σ = cov(mean_removed, dims=2)

    # 计算特征值(Find eigenvalues and eigenvectors)
    U, S, Vt = svd(Σ)

    # 降维
    Ur = U[:,1:Top_n_feat]
    DataReduce = Ur' * mean_removed

    # 恢复原数据
    Rebuild = (Ur * DataReduce) .+ μ
    return DataReduce, Rebuild
end

labels, images = load_mnist("data/mnist","train")

grid, frames, canvases = canvasgrid((2,2))  # 2 row, 2 columns
imshow(canvases[1,1], reshape(images[:,610],28,28)')

DataReduce, Rebuild = PCA(images, 1)
imshow(canvases[1,2], reshape(Rebuild[:,610],28,28)')

DataReduce, Rebuild = PCA(images, 50)
imshow(canvases[2,1], reshape(Rebuild[:,610],28,28)')

DataReduce, Rebuild = PCA(images, 100)
imshow(canvases[2,2], reshape(Rebuild[:,610],28,28)')

win = Window(grid)
Gtk.showall(win)

MSE = sum((Rebuild .- images) .^ 2 , dims = 1) ./ 784
MSE_means = sum(MSE) ./ 60000

max_stars_repo_name: AlongWY/HIT-WorkShop

max_issues_repo_name: AlongWY/HIT-WorkShop

max_forks_repo_name: AlongWY/HIT-WorkShop