text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/using-a-abstractunitrange-inside-a-struct/64463/5

# %%
using BenchmarkTools

struct Instance_RH{P <: AbstractFloat, Q <: AbstractUnitRange}
    cost::Array{P, 2}
    time_by_quantity::Array{P, 2}
    penalty_by_multitask::Array{P, 1}
    penalty_by_know_how::Array{P, 2}
    maximum_budget::P
    I::Q
    J::Q
end

function generate_instance(cost, time_by_quantity,
    penalty_by_multitask, penalty_by_know_how, maximum_budget)

    dimensions_of_instance = size(cost);
    I = Base.OneTo(dimensions_of_instance[1]);
    J = Base.OneTo(dimensions_of_instance[2]);
    return Instance_RH(cost, time_by_quantity,
        penalty_by_multitask, penalty_by_know_how, maximum_budget, I, J);
end

cost = [404.16130257432536 342.6121400760149 334.65400017129275; 609.9826980212956 451.2171757784507 483.83148674153233; 362.93218686863423 414.41083606045623 675.6532058287046; 498.73524270040207 379.8833522296185 372.12075326006203; 507.1394146512986 343.40004673720284 344.87937542652276; 614.5935044087983 343.38652475868105 666.6455577870877; 654.5084251471844 460.7325334872372 343.51651476383165; 514.5598529059743 445.47832539448785 667.5994573488664; 352.0625777996967 357.63523731788064 426.5800063667885; 633.8737992123354 575.4684437771787 526.8544116378721]
time_by_quantity = [551.149318836251 353.5576871795756 544.7850084991596; 540.8786124298001 347.11253457830435 534.2360337745861; 530.9324781581835 340.806606793647 524.2600861079036; 521.2428271420944 334.758141874596 514.891286519538; 512.0603850942899 329.686814227325 506.0319653251741; 503.2315108697694 324.6932839541254 497.68125911747074; 495.06705853446783 320.35799835061465 489.867352509886; 487.4786234404495 316.2110977723632 483.3234557345035; 481.19585646796173 312.4731479370266 477.0680377885712; 476.0273928951479 309.1389856638628 472.01244831930615]
penalty_by_multitask = [3.932768925120902, 2.1629197821510475, 3.747201375134096, 1.2027401788064473, 3.762680990567876, 3.601216190081945, 3.9412748264981996, 1.004961550315905, 3.963908498229274, 1.5845709053894907]
penalty_by_know_how = [0.811984335664391 0.2852506669586609 0.4118216702970811; 0.10635146869948212 0.3497031430472351 0.6630721652898788; 0.13880963807944877 0.5621105391563797 0.22663788295134837; 0.44358375556142055 0.22274423735250545 0.5015070457420686; 0.14826541117764483 0.7966779051319146 0.6769308830979445; 0.07009524283219473 0.08771913744412921 0.5220504363220271; 0.8324525095520683 0.8185583447816812 0.4939598817833735; 0.597792860354574 0.4440247339167606 0.8233927217289345; 0.014589514521600778 0.29122819336562134 0.6274289868926287; 0.21229629326964108 0.2910770390826134 0.8959322834368744]
maximum_budget = 14109.108387884087
instance_0 = generate_instance(cost, time_by_quantity, penalty_by_multitask, penalty_by_know_how, maximum_budget)

function calculate!(y, x, instance) # Quantidade de pessoas i envolvidas no projeto j
    @inbounds begin
        for j in instance.J
            y[j] = 0;
            for i in instance.I
                y[j] += x[i, j];
            end
        end
    end
    return minimum(y) == 0 ? count(i->(i == 0), y) : 0;
end
x = [1 1 0; 0 0 1; 1 0 1; 0 0 0; 0 0 0; 1 1 1; 0 0 0; 0 0 0; 0 0 0; 0 1 1]
y = zeros(Int, size(x, 2))
@benchmark calculate!($y, $x, $instance_0)

# %%
function calculate_2!(y, x, instance) # Quantidade de pessoas i envolvidas no projeto j
    @inbounds begin
        for j in 1:3
            y[j] = 0;
            for i in 1:10
                y[j] += x[i, j];
            end
        end
    end
    return minimum(y) == 0 ? count(i->(i == 0), y) : 0;
end

@benchmark calculate_2!($y, $x, $instance_0)

# %%
function calculate_3!(y, x, instance) # Quantidade de pessoas i envolvidas no projeto j
    m, n = size(x)
    @inbounds begin
        for j in 1:n
            y[j] = 0;
            for i in 1:m
                y[j] += x[i, j];
            end
        end
    end
    return minimum(y) == 0 ? count(i->(i == 0), y) : 0;
end

@benchmark calculate_3!($y, $x, $instance_0)

# %%
calculate!(y, x, instance_0) == calculate_2!(y, x, instance_0) == calculate_3!(y, x, instance_0)

# %%
using BenchmarkTools

struct Instance_RH{P <: AbstractFloat, Q <: AbstractUnitRange}
    cost::Array{P, 2}
    time_by_quantity::Array{P, 2}
    penalty_by_multitask::Array{P, 1}
    penalty_by_know_how::Array{P, 2}
    maximum_budget::P
    I::Q
    J::Q
end

function generate_instance(cost, time_by_quantity,
    penalty_by_multitask, penalty_by_know_how, maximum_budget)

    dimensions_of_instance = size(cost);
    I = Base.OneTo(dimensions_of_instance[1]);
    J = Base.OneTo(dimensions_of_instance[2]);
    return Instance_RH(cost, time_by_quantity,
        penalty_by_multitask, penalty_by_know_how, maximum_budget, I, J);
end

function calculate!(y, x, instance) # Quantidade de pessoas i envolvidas no projeto j
    @inbounds begin
        for j in instance.J
            y[j] = 0;
            for i in instance.I
                y[j] += x[i, j];
            end
        end
    end
    #return minimum(y) == 0 ? count(i->(i == 0), y) : 0;
    y
end

function calculate_2!(y, x, instance) # Quantidade de pessoas i envolvidas no projeto j
    @inbounds begin
        for j in 1:1000
            y[j] = 0;
            for i in 1:3000
                y[j] += x[i, j];
            end
        end
    end
    #return minimum(y) == 0 ? count(i->(i == 0), y) : 0;
    y
end

function calculate_3!(y, x, instance) # Quantidade de pessoas i envolvidas no projeto j
    m, n = size(x)
    @inbounds begin
        for j in 1:n
            y[j] = 0;
            for i in 1:m
                y[j] += x[i, j];
            end
        end
    end
    #return minimum(y) == 0 ? count(i->(i == 0), y) : 0;
    y
end

m, n = 3000, 1000
instance_1 = generate_instance(rand(m, n), rand(m, n), rand(n), rand(m, n), rand())
x_1 = rand([0, 1], m, n)

y_1 = zeros(Int, size(x_1, 2))
y_12 = zeros(Int, size(x_1, 2))
y_13 = zeros(Int, size(x_1, 2))
@show calculate!(y_1, x_1, instance_1) == calculate_2!(y_12, x_1, instance_1) == calculate_3!(y_13, x_1, instance_1)

@btime calculate!($y_1, $x_1, $instance_1)
@btime calculate_2!($y_12, $x_1, $instance_1)
@btime calculate_3!($y_13, $x_1, $instance_1);

# %%

max_stars_repo_path: 0010/Using a AbstractUnitRange inside a struct.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
@testset "display_obj.jl" begin
    label = DS.Label("我是字串", [3, 5])
    DS.render(label)
    @test T.read_stream(T.out_stream) == "\e[32m\e[3;5H我是字串\e[0m"

    label = DS.Label("我是字串", [3, 5], style=Crayon(foreground=:blue, bold=true))
    DS.render(label)
    @test T.read_stream(T.out_stream) == "\e[34;1m\e[3;5H我是字串\e[0m"

    panel = DS.Panel(" 我是標題 ", [3, 20], [3, 5])
    DS.render(panel)
    @test T.read_stream(T.out_stream) ==
        "\e[34;1m\e[3;5H" *
        "\e[s┌ 我是標題 ────────┐\e[u\e[1B" *
        "\e[s│                  │\e[u\e[1B" *
        "\e[s└──────────────────┘\e[u\e[1B" *
        "\e[0m"

    panel = DS.Panel(" 我是標題 ", [3, 20], [3, 5], style=Crayon(foreground=:green))
    DS.render(panel)
    @test T.read_stream(T.out_stream) ==
        "\e[32m\e[3;5H" *
        "\e[s┌ 我是標題 ────────┐\e[u\e[1B" *
        "\e[s│                  │\e[u\e[1B" *
        "\e[s└──────────────────┘\e[u\e[1B" *
        "\e[0m"

    panel = DS.Panel(" 我是標題 ", [3, 10], [3, 5])
    DS.render(panel)
    @test T.read_stream(T.out_stream) ==
        "\e[34;1m\e[3;5H" *
        "\e[s┌────────┐\e[u\e[1B" *
        "\e[s│        │\e[u\e[1B" *
        "\e[s└────────┘\e[u\e[1B" *
        "\e[0m"
end

ext: jl


text:
module Exp
using StatsBase
using Statistics
include("./observable.jl")

const t = [0, 5, 15, 30, 45, 60, 90, 120]

function norm01(;egf1::Vector{Int}, hrg1::Vector{Int}, egf2::Vector{Int},
    hrg2::Vector{Int}, egf3::Vector{Int}, hrg3::Vector{Int})
    data1::Matrix{Float64} = hcat(egf1, hrg1)
    data2::Matrix{Float64} = hcat(egf2, hrg2)
    data3::Matrix{Float64} = hcat(egf3, hrg3)

    data1 .= data1 ./ maximum(data1)
    data2 .= data2 ./ maximum(data2)
    data3 .= data3 ./ maximum(data3)

    egf_ave::Vector{Float64} = zeros(length(t))
    hrg_ave::Vector{Float64} = zeros(length(t))
    egf_sem::Vector{Float64} = zeros(length(t))
    hrg_sem::Vector{Float64} = zeros(length(t))

    for i in eachindex(t)
        egf_ave[i] = mean([data1[i,1], data2[i,1], data3[i,1]])
        hrg_ave[i] = mean([data1[i,2], data2[i,2], data3[i,2]])
    end

    ave_vec::Vector{Float64} = vcat(egf_ave, hrg_ave)
    ave_min::Float64 = minimum(ave_vec)
    ave_max::Float64 = maximum(ave_vec)

    data1 .= (data1 .- ave_min) ./ (ave_max .- ave_min)
    data2 .= (data2 .- ave_min) ./ (ave_max .- ave_min)
    data3 .= (data3 .- ave_min) ./ (ave_max .- ave_min)

    for i in eachindex(t)
        egf_ave[i] = mean([data1[i,1], data2[i,1], data3[i,1]])
        hrg_ave[i] = mean([data1[i,2], data2[i,2], data3[i,2]])
        egf_sem[i] = std([data1[i,1], data2[i,1], data3[i,1]]) ./ sqrt(3)
        hrg_sem[i] = std([data1[i,2], data2[i,2], data3[i,2]]) ./ sqrt(3)
    end

    return egf_ave, hrg_ave, egf_sem, hrg_sem
end

experiments = Array{Dict{String,Array{Float64,1}},1}(undef, length(observables))
error_bars = Array{Dict{String,Array{Float64,1}},1}(undef, length(observables))

mcf7_data_pAkt = norm01(
    egf1=[32101, 156970, 90301, 76709, 63640, 52536, 46414, 57329],
    hrg1=[32101, 565508, 551901, 560064, 489678, 408802, 425323, 451502],
    egf2=[11612, 96189, 43622, 43238, 41007, 29902, 19255, 35079],
    hrg2=[11612, 397931, 432609, 417622, 434519, 509919, 361041, 292523],
    egf3=[66038, 208525, 102689, 117308, 125158, 92086, 68587, 78252],
    hrg3=[66038, 563079, 573540, 521062, 447462, 383774, 434807, 409615],
)
bt474_data_pAkt = norm01(
    egf1=[405198, 356865, 321475, 383445, 346872, 328052, 299123, 316633],
    hrg1=[405198, 357121, 419948, 488508, 495214, 443710, 402765, 451831],
    egf2=[432524, 619289, 581376, 481899, 429541, 399922, 376170, 334923],
    hrg2=[432524, 410919, 413878, 390581, 405359, 408471, 373108, 515120],
    egf3=[150446, 435897, 466378, 443105, 415827, 381441, 398841, 413906],
    hrg3=[150446, 556176, 560385, 539165, 589297, 552227, 540005, 539010],
)
mdamb231_data_pAkt = norm01(
    egf1=[86491, 826975, 575400, 354446, 143728, 107326, 88082, 108892],
    hrg1=[86491, 85990, 114224, 94636, 66139, 92359, 89540, 126099],
    egf2=[56004, 816102, 296684, 165282, 92711, 83090, 55002, 50605],
    hrg2=[56004, 50640, 64646, 55814, 70627, 66657, 68994, 46922],
    egf3=[36061, 729961, 218315, 166109, 68746, 53331, 45191, 21852],
    hrg3=[36061, 75029, 174050, 89546, 88262, 99161, 78673, 58412],
)
skbr3_data_pAkt = norm01(
    egf1=[185525, 566558, 675230, 679968, 669038, 596447, 547016, 439083],
    hrg1=[185525, 530483, 594649, 637933, 733390, 607193, 649356, 709869],
    egf2=[208726, 431186, 439068, 424423, 452676, 414664, 388795, 310282],
    hrg2=[208726, 467749, 509839, 497772, 470593, 429163, 340911, 314740],
    egf3=[290628, 699828, 586893, 512882, 460166, 432428, 442622, 383550],
    hrg3=[290628, 644516, 632090, 591780, 655548, 689898, 632575, 534658],
)
experiments[observables_index("Phosphorylated_Akt")] = Dict(
    "MCF7_EGF" => mcf7_data_pAkt[1],
    "MCF7_HRG" => mcf7_data_pAkt[2],
    "BT474_EGF" => bt474_data_pAkt[1],
    "BT474_HRG" => bt474_data_pAkt[2],
    "MDAMB231_EGF" => mdamb231_data_pAkt[1],
    "MDAMB231_HRG" => mdamb231_data_pAkt[2],
    "SKBR3_EGF" => skbr3_data_pAkt[1],
    "SKBR3_HRG" => skbr3_data_pAkt[2],
)
error_bars[observables_index("Phosphorylated_Akt")] = Dict(
    "MCF7_EGF" => mcf7_data_pAkt[3],
    "MCF7_HRG" => mcf7_data_pAkt[4],
    "BT474_EGF" => bt474_data_pAkt[3],
    "BT474_HRG" => bt474_data_pAkt[4],
    "MDAMB231_EGF" => mdamb231_data_pAkt[3],
    "MDAMB231_HRG" => mdamb231_data_pAkt[4],
    "SKBR3_EGF" => skbr3_data_pAkt[3],
    "SKBR3_HRG" => skbr3_data_pAkt[4],
)

mcf7_data_pERK = norm01(
    egf1=[65481, 446949, 221435, 283171, 265152, 266056, 204912, 188972],
    hrg1=[65481, 698717, 766252, 710005, 693622, 691856, 522173, 334410],
    egf2=[41927, 507623, 169918, 193671, 164088, 145916, 110844, 130362],
    hrg2=[41927, 605118, 699511, 654697, 579863, 490649, 299946, 229297],
    egf3=[118995, 807929, 338665, 267160, 253820, 230200, 157620, 153112],
    hrg3=[118995, 710436, 673318, 615206, 612686, 523198, 390301, 257664],
)
bt474_data_pERK = norm01(
    egf1=[358203, 550378, 633802, 632047, 500267, 394009, 339650, 221411],
    hrg1=[358203, 531893, 703437, 663640, 629213, 612525, 613871, 643056],
    egf2=[355065, 1421536, 1307969, 1101679, 939944, 689539, 507787, 468836],
    hrg2=[355065, 929915, 897601, 924274, 865529, 820386, 673456, 788623],
    egf3=[61593, 631017, 754722, 652440, 575812, 432406, 315961, 259708],
    hrg3=[61593, 480140, 487770, 463604, 438917, 452289, 470624, 531531],
)
mdamb231_data_pERK = norm01(
    egf1=[314472, 504819, 607786, 618492, 475195, 376035, 293988, 324600],
    hrg1=[314472, 156705, 183456, 277862, 141450, 199719, 253331, 407923],
    egf2=[458693, 1001334, 875594, 834259, 782639, 815888, 629576, 539187],
    hrg2=[458693, 322542, 403985, 331734, 263578, 262142, 276371, 313541],
    egf3=[365691, 747932, 937413, 945635, 870059, 706306, 510590, 451927],
    hrg3=[365691, 340876, 428510, 303543, 269653, 195660, 215972, 388050],
)
skbr3_data_pERK = norm01(
    egf1=[182424, 757623, 702797, 710750, 704235, 636329, 680400, 637799],
    hrg1=[182424, 662752, 754617, 725132, 637076, 590793, 610035, 541614],
    egf2=[335742, 721244, 699229, 542494, 486792, 449850, 446215, 452925],
    hrg2=[335742, 588327, 607123, 557027, 490253, 485330, 478285, 530225],
    egf3=[153334, 617406, 425659, 408737, 358565, 338492, 350159, 353366],
    hrg3=[153334, 636281, 559026, 436561, 389041, 360622, 333743, 355748],
)
experiments[observables_index("Phosphorylated_ERK")] = Dict(
    "MCF7_EGF" => mcf7_data_pERK[1],
    "MCF7_HRG" => mcf7_data_pERK[2],
    "BT474_EGF" => bt474_data_pERK[1],
    "BT474_HRG" => bt474_data_pERK[2],
    "MDAMB231_EGF" => mdamb231_data_pERK[1],
    "MDAMB231_HRG" => mdamb231_data_pERK[2],
    "SKBR3_EGF" => skbr3_data_pERK[1],
    "SKBR3_HRG" => skbr3_data_pERK[2],
)
error_bars[observables_index("Phosphorylated_ERK")] = Dict(
    "MCF7_EGF" => mcf7_data_pERK[3],
    "MCF7_HRG" => mcf7_data_pERK[4],
    "BT474_EGF" => bt474_data_pERK[3],
    "BT474_HRG" => bt474_data_pERK[4],
    "MDAMB231_EGF" => mdamb231_data_pERK[3],
    "MDAMB231_HRG" => mdamb231_data_pERK[4],
    "SKBR3_EGF" => skbr3_data_pERK[3],
    "SKBR3_HRG" => skbr3_data_pERK[4],
)

mcf7_data_pcMyc = norm01(
    egf1=[115975, 226001, 166894, 194150, 263331, 235172, 126949, 91142],
    hrg1=[115975, 62515, 81364, 155844, 390689, 664641, 848356, 856941],
    egf2=[185069, 276202, 204012, 234391, 290020, 360762, 325531, 242455],
    hrg2=[185069, 234416, 251732, 333993, 550670, 859790, 939956, 769616],
    egf3=[127244, 186118, 163387, 132053, 192949, 220987, 184381, 151547],
    hrg3=[127244, 110676, 152880, 277206, 461217, 637033, 908235, 712427],
)
bt474_data_pcMyc = norm01(
    egf1=[677642, 631653, 761828, 940801, 824185, 1063834, 1248333, 1126031],
    hrg1=[677642, 672273, 647219, 804473, 958061, 1105501, 1088048, 1323065],
    egf2=[382898, 389470, 359455, 414052, 501284, 540238, 621011, 445658],
    hrg2=[382898, 299678, 327189, 491156, 762165, 898141, 835015, 684669],
    egf3=[321190, 366660, 301538, 338201, 395081, 377184, 334796, 277554],
    hrg3=[321190, 208981, 211540, 332880, 434028, 539873, 598716, 603566],
)
mdamb231_data_pcMyc = norm01(
    egf1=[175088, 147000, 158258, 201315, 324002, 362181, 481797, 475243],
    hrg1=[175088, 239248, 242011, 307948, 345390, 408564, 454948, 443210],
    egf2=[398355, 426051, 420120, 513399, 483150, 516812, 529569, 506936],
    hrg2=[398355, 320668, 323452, 315337, 321620, 365778, 434971, 439945],
    egf3=[887432, 906168, 503726, 913939, 1236533, 1305121, 1000144, 955499],
    hrg3=[887432, 803746, 822021, 739253, 833273, 844401, 844497, 1064656],
)
skbr3_data_pcMyc = norm01(
    egf1=[216306, 260513, 398333, 450702, 457053, 431985, 360345, 234876],
    hrg1=[216306, 189928, 337528, 333291, 370834, 298729, 263847, 198003],
    egf2=[1100082, 874469, 1105692, 1121651, 1132333, 984411, 861916, 928244],
    hrg2=[1100082, 620915, 847173, 1070530, 849862, 753615, 516349, 486978],
    egf3=[624420, 636816, 875382, 972561, 1066021, 916545, 669543, 627495],
    hrg3=[624420, 511242, 678598, 826159, 933370, 839597, 600985, 610415],
)
experiments[observables_index("Phosphorylated_c-Myc")] = Dict(
    "MCF7_EGF" => mcf7_data_pcMyc[1],
    "MCF7_HRG" => mcf7_data_pcMyc[2],
    "BT474_EGF" => bt474_data_pcMyc[1],
    "BT474_HRG" => bt474_data_pcMyc[2],
    "MDAMB231_EGF" => mdamb231_data_pcMyc[1],
    "MDAMB231_HRG" => mdamb231_data_pcMyc[2],
    "SKBR3_EGF" => skbr3_data_pcMyc[1],
    "SKBR3_HRG" => skbr3_data_pcMyc[2],
)
error_bars[observables_index("Phosphorylated_c-Myc")] = Dict(
    "MCF7_EGF" => mcf7_data_pcMyc[3],
    "MCF7_HRG" => mcf7_data_pcMyc[4],
    "BT474_EGF" => bt474_data_pcMyc[3],
    "BT474_HRG" => bt474_data_pcMyc[4],
    "MDAMB231_EGF" => mdamb231_data_pcMyc[3],
    "MDAMB231_HRG" => mdamb231_data_pcMyc[4],
    "SKBR3_EGF" => skbr3_data_pcMyc[3],
    "SKBR3_HRG" => skbr3_data_pcMyc[4],
)

function get_timepoint(obs_name::String)::Vector{Float64}
    if obs_name in observables
        return t
    end
end
end # module


text:
include("6.1.10 (c+d x)^m (a+b sinh)^n.jl")
include("6.1.11 (e x)^m (a+b x^n)^p sinh.jl")
include("6.1.12 (e x)^m (a+b sinh(c+d x^n))^p.jl")
include("6.1.13 (d+e x)^m sinh(a+b x+c x^2)^n.jl")

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.1 Hyperbolic sine/.jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AdlsError <: SwaggerModel
    RemoteException::Any # spec type: Union{ Nothing, AdlsRemoteException } # spec name: RemoteException

    function AdlsError(;RemoteException=nothing)
        o = new()
        validate_property(AdlsError, Symbol("RemoteException"), RemoteException)
        setfield!(o, Symbol("RemoteException"), RemoteException)
        o
    end
end # type AdlsError

const _property_map_AdlsError = Dict{Symbol,Symbol}(Symbol("RemoteException")=>Symbol("RemoteException"))
const _property_types_AdlsError = Dict{Symbol,String}(Symbol("RemoteException")=>"AdlsRemoteException")
Base.propertynames(::Type{ AdlsError }) = collect(keys(_property_map_AdlsError))
Swagger.property_type(::Type{ AdlsError }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_AdlsError[name]))}
Swagger.field_name(::Type{ AdlsError }, property_name::Symbol) =  _property_map_AdlsError[property_name]

function check_required(o::AdlsError)
    true
end

function validate_property(::Type{ AdlsError }, name::Symbol, val)
end


text:
include("testsetup.jl")

include("test-suite-preamble.jl")

include("internals.jl")
include("threadingutilities.jl")
if (!parse(Bool, get(ENV, "GITHUB_ACTIONS", "false"))) && Threads.nthreads() > 3
    include("staticarrays.jl")
end
include("threadpool.jl")
include("warnings.jl")

include("aqua.jl") # run the Aqua.jl tests last


text:
# Convert a number of molars and a volume to a number of molecules
# Arguments:  m = the number of molars
#             v = the volume in metres cubed
# Returns:    Int, the number of molecules

function molarsToMolecules(m, v)
  N_av = 6.02214085774e+23
  v_litres = 1000 * v
  num_molecules = try
    floor(Int, N_av * m * v_litres)
  catch
    floor(N_av * m * v_litres)
  end
  return num_molecules
end


text:
function ccArm(k, l)
    #check if we have enough curvature points
    n=length(k)

    #init segments cords
    p=zeros(n+1, 3) #[x, y, θ]
    for i=1:n #each segment, need 6 sets of points
        p[i+1, 3]=p[i, 3]+l*k[i] #update angle
        p[i+1, 1]=p[i, 1]+(sin(p[i+1, 3])-sin(p[i,3]))/k[i] #update x
        p[i+1, 2]=p[i, 2]+(-cos(p[i+1, 3])+cos(p[i,3]))/k[i] #update y
    end
    return p
end

function cantArm(p, l=1, EI=1, res=10)
    n=length(p)
    #init segments cords
    arm=zeros(n*res+1, 3) #[x, y, θ]
    for i=1:n #each segment,
        #@show i
        c=arm[(i-1)*res+1,:] #base of the segment
        for j=1:res
            inx=(i-1)*res+j
            dx=(1.0*j)/res*l
            arm[inx+1,:]=fwdKin(dx, p[i], l, c)

            #arm[i+1, 3]=arm[i, 3]+p[i]*l^3/(6*EI) #update angle
            #arm[i+1, 1]=arm[i, 1]+cos(arm[i, 3])*l-sin(arm[i, 3])*p[i]*l^4/(8*EI) #update x
            #arm[i+1, 2]=arm[i, 2]+sin(arm[i, 3])*l+cos(arm[i, 3])*p[i]*l^4/(8*EI) #update y
        end
    end
    return arm

end
function rotate(t)
    return [cos(t) -sin(t); sin(t) cos(t)]
end

function curv(x, q, l=1, EI=1)
    return q*(3*l^2*x-3*l*x^2+x^3)/(6*EI)
end

function deflec(x, q, l=1, EI=1)
    return q*x^2*(6*l^2-4*l*x+x^2)/(24*EI)
end

function fwdKin(dx, q, l, cord)
    t=curv(dx, q, l) + cord[3]
    w=deflec(dx, q, l)
    c=rotate(t)* [dx;w]+ cord[1:2]
    return vcat(c, t)
end


text:
# type definition
immutable GenoMatrix
    x          :: SharedVector{Int8}
    xt         :: SharedVector{Int8}
    n          :: Int
    p          :: Int
    blocksize  :: Int
    tblocksize :: Int
    snpids     :: Vector{String}

    GenoMatrix(x,xt,n,p,blocksize,tblocksize,snpids) = new(x,xt,n,p,blocksize,tblocksize,snpids)
end

# additional constructors, all of which will infer snpids field from BIM file
# consequently the snpids field never appers as an argument
# first constructor uses all remaining fields of GenoMatrix as arguments
function GenoMatrix(
    filename   :: String,
    tfilename  :: String,
    n          :: Int,
    p          :: Int,
    blocksize  :: Int,
    tblocksize :: Int;
    pids       :: Vector{Int} = procs(),
)
    # find SNP ids from the corresponding BIM file
    # first create filepath to BIM file
    bimfile = filename[1:(endof(filename)-3)] * "bim"

    # specify column element types of BIM file
    #eltypes = [Int, String, Int, Int, String, String]

    # load BIM
    #df = readtable(bimfile, header = false, separator = '\t', eltypes = eltypes)

    # second column of BIM file contains the SNP ids
    #snpids = convert(Vector{String}, df[:x2]) :: Vector{String}
    snpids = convert(Vector{String}, split(chomp(readstring(`awk -F " " '{print $2}' $bimfile`)), '\n', keep=false))

    x = GenoMatrix(
        read_bedfile(filename, pids=pids),
        read_bedfile(tfilename, transpose = true, pids=pids),
        n, p, blocksize, tblocksize, snpids)

    return x
end

# another constructor for when blocksizes are not precomputed
function GenoMatrix(
    filename  :: String,
    tfilename :: String,
    n         :: Int,
    p         :: Int;
    pids      :: Vector{Int} = procs()
)
    blocksize  = ( (n-1) >>> 2) + 1
    tblocksize = ( (p-1) >>> 2) + 1
    return GenoMatrix(filename, tfilename, n, p, blocksize, tblocksize, pids=pids) :: PLINK.GenoMatrix
end

# constructor to load entirely from file
function GenoMatrix(
    filename  :: String,
    tfilename :: String;
    pids      :: Vector{Int} = procs()
)
    # find n from the corresponding FAM file
    famfile = filename[1:(endof(filename)-3)] * "fam"
    n       = countlines(famfile)

    # find p from the corresponding BIM file
    bimfile = filename[1:(endof(filename)-3)] * "bim"
    p       = countlines(bimfile)

    # blocksizes are easy to calculate
    blocksize  = ((n-1) >>> 2) + 1
    tblocksize = ((p-1) >>> 2) + 1

    # can now create GenoMatrix
    return GenoMatrix(filename, tfilename, n, p, blocksize, tblocksize, pids=pids) :: PLINK.GenoMatrix
end

# subroutines
Base.size(x::GenoMatrix) = (x.n, x.p)
Base.size(x::GenoMatrix, d::Int) = d == 1 ? x.n : x.p

Base.length(x::GenoMatrix) = x.n * x.p

Base.eltype(x::GenoMatrix) = Int8

Base.ndims(x::GenoMatrix) = 2

Base.copy(x::GenoMatrix) = GenoMatrix(x.x,x.xt,x.n,x.p,x.blocksize,x.tblocksize)

Base.IndexStyle(x::Type{GenoMatrix}) = Base.IndexLinear()

Base.procs(x::Type{GenoMatrix}) = procs(x.x)

function ==(x::GenoMatrix, y::GenoMatrix)
    x.x  == y.x  &&
    x.xt == y.xt &&
    x.n  == y.n  &&
    x.p  == y.p  &&
    x.blocksize  == y.blocksize &&
    x.tblocksize == y.tblocksize
end

isequal(x::GenoMatrix, y::GenoMatrix) = x == y

# important matrix indexing! returns the Int8
function getindex(
    X   :: GenoMatrix,
    x   :: SharedVector{Int8},
    row :: Int,
    col :: Int
)
    genotype_block = x[(col-1)*X.blocksize + ((row - 1) >>> 2) + 1]
    k = 2*((row-1) & 3)
    genotype = (genotype_block >>> k) & THREE8
    return genotype
end

# internal routine to index the transposed matrix
# used for dott()
function getindex_t(
    X   :: GenoMatrix,
    x   :: SharedVector{Int8},
    row :: Int,
    col :: Int
)
    genotype_block = x[(col-1)*X.tblocksize + ((row - 1) >>> 2) + 1]
    k = 2*((row-1) & 3)
    genotype = (genotype_block >>> k) & THREE8
    return genotype
end

# default indexing is to column-major matrix
# can use x.xt to index the transpose
getindex(x::GenoMatrix, row::Int, col::Int) = getindex(x, x.x, row, col)

"""
    read_bedfile(filename [, transpose=false, pids=procs]) -> SharedVector{Int8}

This function reads a PLINK binary file (BED) and returns a `SharedArray` of `Int8` numbers.
It discards the first three bytes ("magic numbers") since they are not needed here.

Arguments:

- `filename` is the path to the BED file

Optional Arguments:

- `transpose` indicates if the compressed matrix is transposed. Defaults to `false`.
- `pids` indicates the processes over which to distribute the output `SharedArray`. Defaults to `procs()` (all available processes).

Output:

- A vector of `Int8` numbers. For a BED file encoding `n` cases and `p` SNPs,
  there should be *at least* `(n*p/4)` numbers. The scaling factor of `4` comes from the
  compression of four genotypes into each byte. But PLINK stores each column in blocks
  of bytes instead of a continuous bitstream, which sometimes entails extra unused bits
  at the end of each block.
"""
function read_bedfile(
    filename  :: String;
    transpose :: Bool = false,
    pids      :: Vector{Int} = procs()
)

    # check that file is BED file
    @assert contains(filename, ".bed") "Filename must explicitly point to a PLINK \".bed\" file."

#    # check magic numbers and mode
#   isequal(read(xstream, Int8), MNUM1) || throw(error("Problem with first byte of magic number, is this a true BED file?"))
#   isequal(read(xstream, Int8), MNUM2) || throw(error("Problem with second byte of magic number, is this a true BED file?"))
#   (transpose && isequal(read(xstream, Int8), ONE8)) && throw(error("For transposed matrix, third byte of BED file must indicate individual-major format."))

    # how many bytes do we have?
    nbytes = filesize(filename)

    # file seems to be a true BED file
    # will close filestream and slurp entire file into SharedArray
    x = SharedArray{Int8}(abspath(filename), (nbytes-3,), 3, pids=pids) :: SharedVector{Int8}

    # return the genotypes
    return x
end

function display(x::GenoMatrix)
    println("A compressed GenoMatrix object with the following features:")
    println("\tnumber of cases        = $(x.n)")
    println("\tgenetic covariates     = $(x.p)")
end

max_stars_repo_name: klkeys/PLINK

max_issues_repo_name: klkeys/PLINK


text:
module RefShapes

using LinearAlgebra

"""
    AbstractRefShape{MANIFDIM}

Abstract type of a reference shape.
"""
abstract type AbstractRefShape{MANIFDIM} end

"""
    RefShapePoint <: AbstractRefShape{0}

Type of a reference shape for a zero-dimensional manifold (point).
"""
struct RefShapePoint <: AbstractRefShape{0} end

"""
    RefShapeInterval <: AbstractRefShape{1}

Type of a reference shape for a 1-dimensional manifold (curve).
"""
struct RefShapeInterval <: AbstractRefShape{1} end

"""
    RefShapeSquare <: AbstractRefShape{2}

Type of a logically rectangular reference shape for a 2-dimensional manifold
(surface).
"""
struct RefShapeSquare <: AbstractRefShape{2} end

"""
    RefShapeCube <: AbstractRefShape{3}

Type of a reference shape for a 3-dimensional manifold (solid) bounded by six quadrilaterals.
"""
struct RefShapeCube <: AbstractRefShape{3} end

"""
    RefShapeTriangle <: AbstractRefShape{2}

Type of a logically triangular reference shape for a 2-dimensional manifold
(surface).
"""
struct RefShapeTriangle <: AbstractRefShape{2} end

"""
    RefShapeTetrahedron <: AbstractRefShape{3}

Type of a reference shape for a 3-dimensional manifold (solid) bounded by 4 triangles.
"""
struct RefShapeTetrahedron <: AbstractRefShape{3} end

"""
    manifdim(rs)

Get the manifold dimension of the reference shape.
"""
manifdim(::Type{T}) where {T<:AbstractRefShape{MANIFDIM}} where {MANIFDIM} = MANIFDIM

"""
    manifdimv(::Type{T}) where {T<:AbstractRefShape{MANIFDIM}} where {MANIFDIM}

Get the manifold dimension of the reference shape as `Val`.
"""
manifdimv(::Type{T}) where {T<:AbstractRefShape{MANIFDIM}} where {MANIFDIM} = Val(MANIFDIM)

"""
    IntegRule

Type for integration rule.
"""
struct IntegRule
    npts::Int64
    param_coords::Array{Float64, 2}
    weights::Array{Float64, 1}
end

npts(qr) = qr.npts
param_coords(qr) = qr.param_coords
weights(qr) = qr.weights

function _gauss1(order)
    function opengauss(N, a = -1.0, b = +1.0)
        F = eigen(SymTridiagonal(zeros(N), [n/sqrt(4n^2 - 1) for n = 1:N-1]))
        return [(F.values[i]+1)*(b-a)/2 + a for i = 1:N], [2*F.vectors[1, i]^2 for i = 1:N]*(b-a)/2
    end

    if     (order==1)
        param_coords = vec([ 0.0 ]);
        weights  = vec([ 2.0 ]);
    elseif (order==2)
        param_coords = vec([ -0.577350269189626 0.577350269189626 ]);
        weights  = vec([ 1.0 1.0 ]);
    elseif (order==3)
        param_coords = vec([ -0.774596669241483  0.0  0.774596669241483 ]);
        weights  = vec([ 0.5555555555555556 0.8888888888888889 0.5555555555555556 ]);
    elseif (order==4)
        param_coords = vec([ -0.86113631159405  -0.33998104358486   0.33998104358486   0.86113631159405]);
        weights  = vec([ 0.34785484513745   0.65214515486255   0.65214515486255   0.34785484513745]);
    elseif (order==5)
        param_coords = vec([ -0.906179845938664        -0.538469310105683         0.000000000000000         0.538469310105683         0.906179845938664])
        weights  = vec([0.236926885056189        0.478628670499367        0.568888888888889        0.478628670499367        0.236926885056189])
    else
        param_coords, weights  = opengauss(order)
    end
    return length(param_coords), param_coords, weights
end

function _triangle(npts=1)
    if npts == 1 # integrates exactly linear polynomials
        param_coords = [1.0/3. 1.0/3.];
        weights = reshape([1.0]/2.0,1,1);
    elseif npts == 3 # integrates exactly quadratic polynomials
        param_coords = [ 2.0/3 1.0/6; 1.0/6 2.0/3; 1.0/6 1.0/6 ];
        weights = [1.0/3 1.0/3 1.0/3]/2;
    elseif npts == 4 # integrates exactly quadratic polynomials
        param_coords = [0.333333333333333   0.333333333333333
                        0.200000000000000   0.200000000000000
                        0.600000000000000   0.200000000000000
                        0.200000000000000   0.600000000000000];
        weights = [-0.281250000000000
                    0.260416666666667
                    0.260416666666667
                    0.260416666666667];
    elseif npts == 6 # integrates exactly quartic polynomials
        param_coords = [ 0.816847572980459 0.091576213509771;
                        0.091576213509771 0.816847572980459;
                        0.091576213509771 0.091576213509771;
                        0.108103018168070 0.445948490915965;
                        0.445948490915965 0.108103018168070;
                        0.445948490915965 0.445948490915965];
        weights = [0.109951743655322*[1, 1, 1] 0.223381589678011*[1, 1, 1] ]/2;
    elseif npts == 7 # integrates exactly ? polynomials
        param_coords = [0.101286507323456   0.101286507323456
                        0.797426958353087   0.101286507323456
                        0.101286507323456   0.797426958353087
                        0.470142064105115   0.470142064105115
                        0.059715871789770   0.470142064105115
                        0.470142064105115   0.059715871789770
                        0.333333333333333   0.333333333333333];
        weights = [0.062969590272414
                    0.062969590272414
                    0.062969590272414
                    0.066197076394253
                    0.066197076394253
                    0.066197076394253
                    0.112500000000000];
    elseif npts == 9 # integrates exactly ? polynomials
        param_coords = [   0.437525248383384   0.437525248383384
            0.124949503233232   0.437525248383384
            0.437525248383384   0.124949503233232
            0.165409927389841   0.037477420750088
            0.037477420750088   0.165409927389841
            0.797112651860071   0.165409927389841
            0.165409927389841   0.797112651860071
            0.037477420750088   0.797112651860071
            0.797112651860071   0.037477420750088];
        weights = [   0.205950504760887
            0.205950504760887
            0.205950504760887
            0.063691414286223
            0.063691414286223
            0.063691414286223
            0.063691414286223
            0.063691414286223
            0.063691414286223] ./ 2;
    elseif npts == 12 # integrates exactly ? polynomials
        param_coords = [   0.063089014491502   0.063089014491502
            0.873821971016996   0.063089014491502
            0.063089014491502   0.873821971016996
            0.249286745170910   0.249286745170910
            0.501426509658179   0.249286745170910
            0.249286745170910   0.501426509658179
            0.310352451033785   0.053145049844816
            0.053145049844816   0.310352451033785
            0.636502499121399   0.310352451033785
            0.310352451033785   0.636502499121399
            0.053145049844816   0.636502499121399
            0.636502499121399   0.053145049844816];
        weights = [0.050844906370207
            0.050844906370207
            0.050844906370207
            0.116786275726379
            0.116786275726379
            0.116786275726379
            0.082851075618374
            0.082851075618374
            0.082851075618374
            0.082851075618374
            0.082851075618374
            0.082851075618374] ./ 2;
    elseif npts == 13 # integrates exactly ? polynomials
        param_coords = [0.333333333333333  0.333333333333333
                        0.479308067841923  0.260345966079038
                        0.260345966079038  0.479308067841923
                        0.260345966079038  0.260345966079038
                        0.869739794195568  0.065130102902216
                        0.065130102902216  0.869739794195568
                        0.065130102902216  0.065130102902216
                        0.638444188569809  0.312865496004875
                        0.638444188569809  0.048690315425316
                        0.312865496004875  0.638444188569809
                        0.312865496004875  0.048690315425316
                        0.048690315425316  0.638444188569809
                        0.048690315425316  0.312865496004875
        ];
        weights = [ -0.149570044467670
                    0.175615257433204
                    0.175615257433204
                    0.175615257433204
                    0.053347235608839
                    0.053347235608839
                    0.053347235608839
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
                    0.077113760890257
        ]'/2;
    else
        #nothing doing: this input is wrong
        error( "Unknown number of integration points $(npts)" )
    end
    return npts, param_coords, weights
end

function _tetrahedron(npts=1)
    if npts == 1 # integrates exactly linear polynomials
        param_coords = reshape([0.25,0.25,0.25],1,3);
        weights = reshape([1.0]/6.0,1,1);
    elseif npts == 4 # integrates exactly quadratic polynomials
        param_coords = [[0.13819660 0.13819660 0.13819660];
        [0.58541020 0.13819660 0.13819660];
        [0.13819660 0.58541020 0.13819660];
        [0.13819660 0.13819660 0.58541020]];;
        weights = [ 0.041666666666666666667   0.041666666666666666667   0.041666666666666666667   0.041666666666666666667];
    elseif npts == 5 #  Zienkiewicz #3.
        a =   1.0 / 6.0;
        b =   0.25;
        c =   0.5;
        d = - 0.8;
        e =   0.45;
        param_coords = [[b b b];
        [c a a];
        [a c a];
        [a a c];
        [a a a]];
        weights = [d  e  e  e  e]/6;
    else
        #nothing doing: this input is wrong
        error( "Unknown number of integration points $(npts)" )
    end
    return npts, param_coords, weights
end

"""
    quadrature(::Type{RefShapeInterval}, quadraturesettings = (kind = :default,))

Create a quadrature rule for the reference shape of an interval.

The default is Gauss integration rule, where the order is set with the keyword
`order`.
"""
function quadrature(::Type{RefShapeInterval}, quadraturesettings = (kind = :default,))
	kind = :default
	for apair in pairs(quadraturesettings)
	    sy, val = apair
	    if sy == :kind
	        kind = val
	    end
	end
    if (kind == :Gauss) || (kind == :default)
        # Extract arguments
        order = 1;
        for apair in pairs(quadraturesettings)
            sy, val = apair
            if sy == :order
                order = val
            end
        end
        npts, param_coords, weights = _gauss1(order)
        return IntegRule(npts, reshape(param_coords, size(param_coords, 1), 1), vec(weights))
    else
        error("Integration rule $(kind) not available")
    end
end

"""
    quadrature(::Type{RefShapeTriangle}, quadraturesettings = (kind = :default,))

Create a quadrature rule for the reference shape of an triangle.

The default is a triangle rule, distinguished by the number of points set with
the keyword `npts`.
"""
function quadrature(::Type{RefShapeTriangle}, quadraturesettings = (kind = :default,))
	kind = :default
	for apair in pairs(quadraturesettings)
	    sy, val = apair
	    if sy == :kind
	        kind = val
	    end
	end
    if kind == :default
        # Extract arguments
        npts = 1;
        for apair in pairs(quadraturesettings)
            sy, val = apair
            if sy == :npts
                npts = val
            end
        end
        npts, param_coords, weights = _triangle(npts)
        return IntegRule(npts, reshape(param_coords, size(param_coords, 1), 2), vec(weights))
    else
        error("Integration rule $(kind) not available")
    end
end

"""
    quadrature(::Type{RefShapeSquare}, quadraturesettings = (kind = :default,))

Create a quadrature rule for the reference shape of a square.

The default is Gauss integration rule, where the order is set with the keyword
`order`.
"""
function quadrature(::Type{RefShapeSquare}, quadraturesettings = (kind = :default,))
	kind = :default
	for apair in pairs(quadraturesettings)
	    sy, val = apair
	    if sy == :kind
	        kind = val
	    end
	end
    if (kind == :Gauss) || (kind == :default)
        # Extract arguments
        order = 1;
        for apair in pairs(quadraturesettings)
            sy, val = apair
            if sy == :order
                order = val
            end
        end
        np, pc, w = _gauss1(order)
        param_coords = zeros(eltype(w),np^2,2);
        weights = zeros(eltype(w),np^2);
        r=1
        for i in 1:order
        	for j in 1:order
        		param_coords[r,:] = [pc[i] pc[j]];
        		weights[r] = w[i]*w[j];
        		r=r+1
        	end
        end
        npts = r - 1 # back off one to reach the total number of points
        return IntegRule(npts, param_coords, weights)
    else
        error("Integration rule $(kind) not available")
    end
end

"""
    quadrature(::Type{RefShapeTetrahedron}, quadraturesettings = (kind = :default,))

Create a quadrature rule for the reference shape of a tetrahedron.

The default is a one-point tetrahedron rule; other rules may be chosen based on
the number of points set with the keyword `npts`.
"""
function quadrature(::Type{RefShapeTetrahedron}, quadraturesettings = (kind = :default,))
    kind = :default
    for apair in pairs(quadraturesettings)
        sy, val = apair
        if sy == :kind
            kind = val
        end
    end
    if kind == :default
        # Extract arguments
        npts = 1;
        for apair in pairs(quadraturesettings)
            sy, val = apair
            if sy == :npts
                npts = val
            end
        end
        npts, param_coords, weights = _tetrahedron(npts)
        return IntegRule(npts, reshape(param_coords, size(param_coords, 1), 3), vec(weights))
    else
        error("Integration rule $(kind) not available")
    end
end

end # module


text:
@testset "type inference by use" begin
cst = parse_and_pass("""
struct T
end

struct S
end

f(x::T) = 1
g(x::S) = 1

function ex1(x)
    f(x)
end

function ex2(x)
    f(x)
    g(x)
end

function ex3(x)
    if 1
        f(x)
    else
        f(x)
    end
end

function ex4(x)
    x
    if 1
        f(x)
    end
end

function ex5(x)
    if 1
        f(x)
    else
        g(x)
    end
end

function ex6(x)
    if 1
        y = x
        f(y)
    else
        g(x)
    end
end
""");

T = cst.meta.scope.names["T"]
S = cst.meta.scope.names["S"]

@test cst.meta.scope.names["ex1"].val.meta.scope.names["x"].type == T
@test cst.meta.scope.names["ex2"].val.meta.scope.names["x"].type === nothing
@test cst.meta.scope.names["ex3"].val.meta.scope.names["x"].type === nothing
@test cst.meta.scope.names["ex4"].val.meta.scope.names["x"].type === nothing
@test cst.meta.scope.names["ex5"].val.meta.scope.names["x"].type === nothing
@test cst.meta.scope.names["ex6"].val.meta.scope.names["y"].type === T
end

@testset "loop iterator inference" begin
cst = parse_and_pass("""
begin
abstract type T end
X = Int[]
Y = T[]
end

for x in 1 end
for x in "abc" end
for x in 1:10 end
for x in 1.0:10.0 end
for x in Int[1,2,3] end
for x in X end
for y in Y end
""");

@test cst.args[2].meta.scope.names["x"].type === nothing
@test StaticLint.CoreTypes.ischar(cst.args[3].meta.scope.names["x"].type)
@test StaticLint.CoreTypes.isint(cst.args[4].meta.scope.names["x"].type)
@test StaticLint.CoreTypes.isfloat(cst.args[5].meta.scope.names["x"].type)
@test StaticLint.CoreTypes.isint(cst.args[6].meta.scope.names["x"].type)
@test StaticLint.CoreTypes.isint(cst.args[7].meta.scope.names["x"].type)
@test StaticLint.CoreTypes.isint(cst.args[7].meta.scope.names["x"].type)
@test cst.args[8].meta.scope.names["y"].type == cst.meta.scope.names["T"]
end

@testset "Vector{T} infer" begin
cst = parse_and_pass("""
struct T
    t1
end
struct S
    s1::Vector{T}
end

function f(s::S)
    t = s.s1[1]
    t # This should be inferred as T
end
""")

@test cst[3].meta.scope.names["t"].type == cst.meta.scope.names["T"]
end


text:
# Zero mean function

"""
    MeanZero <: Mean

Zero mean function
```math
m(x) = 0.
```
"""
struct MeanZero <: Mean end

num_params(::MeanZero) = 0
grad_mean(::MeanZero, ::VecF64) = Float64[]
mean(::MeanZero, ::VecF64) = 0.0
mean(mZero::MeanZero, X::MatF64) =  fill(0.0, size(X, 2))
get_params(::MeanZero) = Float64[]
get_param_names(::MeanZero) = Symbol[]

function set_params!(::MeanZero, hyp::VecF64)
    length(hyp) == 0 || throw(ArgumentError("Zero mean function has no parameters"))
end


text:
using FredData
using Base.Test

f = Fred()
s1 = get_data(f, "GDPC1")
s2 = get_data(f, "GDPC1"; units="chg", vintage_dates="2015-01-01")

# Bad requests
@test_throws Exception get_data(f, "GDPC1"; limit="foo")
@test_throws Exception get_data(f, "GDPC1"; vintage_dates="bar")

nothing


text:
include("../src/geom_module.jl")

module dtest
using Random
using BenchmarkTools
using ..geom
using PyPlot
PyPlot.pygui(true)

# x_ = vec([0 1 0.3 1.5 1.7 -0.9 -1])
# y_ = vec([ 0 0.2 1 1.5 3 2.55 0.65])

N = 1000
R = 10.0
# # R_x = 30
# # R_y = 3
# #
# # phi = LinRange(0, 2*π, N)
# # x_ = R_x .* map(cos, phi)
# # y_ = R_y .* map(sin, phi)
#
r_ = R .* Random.rand(N)
phi_ = (2*π) .* Random.rand(N)
x_ = r_ .* cos.(phi_)
y_ = r_ .* sin.(phi_)

point_array = map(geom.MutablePoint, x_, y_)

dmesh = geom.build_delaunay_mesh!(point_array)

point_array2 = copy(point_array)
rnb_idx = Random.randperm(length(point_array2))

(p_x, p_y) = geom.discretize_for_plot(point_array)
(init_tri_x, init_tri_y) = geom.discretize_for_plot(
  [dmesh.vertex_array[1], dmesh.vertex_array[2], dmesh.vertex_array[3]], true)

tri_vec = geom.triangle_vector_for_plot(dmesh)

PyPlot.figure()
PyPlot.grid()
PyPlot.plot(p_x, p_y,
    color = "k",
    alpha = 0.7,
    linestyle = "none",
    marker = ".",
    markersize = 10)
PyPlot.plot(init_tri_x, init_tri_y,
  marker = "o", linewidth = 2)

for kk = 1:length(tri_vec)
    (p_x, p_y) = geom.discretize_for_plot(tri_vec[kk])
    PyPlot.plot(p_x, p_y,
      marker = "o",
      linewidth = 1,
      color = "r",
      alpha = 0.7,)
end

end # module dtest

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# # Particle Set
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/random_flow_field.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/random_flow_field.ipynb)
#
# Simulate trajectories of a particle cloud in a two-dimensional flow field.
# A doubly-periodic domain and randomly-generated flow fields are initially used.
# For additional documentation e.g. see :
# [1](https://JuliaClimate.github.io/IndividualDisplacements.jl/dev/),
# [2](https://JuliaClimate.github.io/MeshArrays.jl/dev/),
# [3](https://docs.juliadiffeq.org/latest/solvers/ode_solve.html),
# [4](https://en.wikipedia.org/wiki/Displacement_(vector))
#
# Exercises:
# - change the initial distribution of partices
# - increase the duration of the trajectories simulation
# - treat the non-periodic domain case by padding `u,v` with zeros
# - replace `u,v` with your own two-dimensional flow fields
#
# ![particles in random flow](https://github.com/JuliaClimate/IndividualDisplacements.jl/raw/master/examples/figs/RandomFlow.gif)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 1. Import Software

using IndividualDisplacements, DataFrames
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/helper_functions.jl"))

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 2. Flow Field
#
# The `u,v` arrays below can be replaced with any other pair provided by the user.
# A couple of important considerations, however:

# - `u,v` are staggered on a C-grid; by `-0.5` grid point in direction `1` for `u` (`2` for `v`)
#  from the grid cell center (0.5,0.5)
# - `u,v` here derive from streamfunction `ϕ`, defined at the corner point, which ensures that
#  the resulting `u,v` is non-divergent, purely rotational, over the C-grid domain.
# In brief:
#
# ```
# u=-(circshift(ϕ, (0,-1))-ϕ)
# v=(circshift(ϕ, (-1,0))-ϕ)
# ```

u,v,ϕ=setup_random_flow()

# If user were to start with collocated velocity (`uC,vC` at the grid cell center) then
# one can easily obtain the staggered velocity (`u,v`) as follows. These may contain both
# [rotational and divergent](https://en.wikipedia.org/wiki/Helmholtz_decomposition) components.
#
# ```
# u=0.5*(circshift(uC, (0,1))+uC)
# v=0.5*(circshift(vC, (1,0))+vC)
# ```

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 3. Initialize Individuals
#
# For example, we can initialize 100 particles within a central subdomain as follows.

np,nq=size(u)
x=np*(0. .+ 1.0*rand(1000))
y=nq*(0. .+ 1.0*rand(1000));

# The `setup_point_cloud` function then wraps everything in the `Individuals` data structure.

𝐼=setup_point_cloud(u,v,X=x,Y=y)
𝐼.𝑃.𝑇[2]=10.
𝐼.🔴

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 3. Compute Trajectories
#
# The time period is `𝐼.𝑃.𝑇` by default, unless `∫!(𝐼,𝑇)` is called instead.

∫!(𝐼)

#nb # %% {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ## 4. Plot Results
#
# For example, generate a simple animation:

#md p=dirname(pathof(IndividualDisplacements))
#md include(joinpath(p,"../examples/recipes_plots.jl"));

#md 🔴_by_t = groupby(𝐼.🔴, :t)
#md anim = @animate for t in eachindex(🔴_by_t)
#md    phi_scatter(ϕ,🔴_by_t[t])
#md end

#md pth=tempdir()*"/"
#md gif(anim, pth*"RandomFlow.gif", fps = 10)


text:
#numerator of asteroid radial velocity wrt Earth
function rvelea(dx, x, params, t)
    jd0 = params[4] # Julian date of start time
    dsj2k = t+(jd0-JD_J2000) # days since J2000.0 = 2.451545e6
    ss16asteph_t = evaleph(params[1], dsj2k, x[1]) # params[2](t)*one(q[1]) #ss16asteph(t)
    N = params[6]
    xe = ss16asteph_t[nbodyind(N-1,ea)]
    return true, (x[1]-xe[1])*(x[4]-xe[4]) + (x[2]-xe[2])*(x[5]-xe[5]) + (x[3]-xe[3])*(x[6]-xe[6])
end

function loadeph(ss16asteph_::TaylorInterpolant, μ::Vector)
    # read Solar System ephemeris (Sun+8 planets+Moon+Pluto+16 main belt asteroids)
    ss16asteph_t0 = (ss16asteph_.t0 ./ daysec)
    ss16asteph_t = (ss16asteph_.t ./ daysec)
    ephord = ss16asteph_.x[1].order
    ss16asteph_x = map(x->x(Taylor1(ephord)*daysec), ss16asteph_.x)
    ss16asteph = TaylorInterpolant(ss16asteph_t0, ss16asteph_t, ss16asteph_x)
    #compute point-mass Newtonian accelerations from ephemeris: all bodies except asteroid
    # accelerations of "everybody else" are needed when evaluating asteroid post-Newtonian acceleration
    Nm1 = (size(ss16asteph_x)[2]-13) ÷ 6
    acc_eph = TaylorInterpolant(ss16asteph.t0, ss16asteph.t, Matrix{eltype(ss16asteph.x)}(undef, length(ss16asteph.t)-1, 3Nm1))
    newtonianNb_Potential = TaylorInterpolant(ss16asteph.t0, ss16asteph.t, Matrix{eltype(ss16asteph.x)}(undef, length(ss16asteph.t)-1, Nm1))
    fill!(acc_eph.x, zero(ss16asteph.x[1]))
    fill!(newtonianNb_Potential.x, zero(ss16asteph.x[1]))
    _1_to_Nm1 = Base.OneTo(Nm1) # iterator over all bodies except asteroid
    for j in _1_to_Nm1
        for i in _1_to_Nm1
            # i == j && continue
            if i == j
            else
                X_ij = ss16asteph.x[:,3i-2] .- ss16asteph.x[:,3j-2]
                Y_ij = ss16asteph.x[:,3i-1] .- ss16asteph.x[:,3j-1]
                Z_ij = ss16asteph.x[:,3i  ] .- ss16asteph.x[:,3j  ]
                r_p2_ij = ( (X_ij.^2) .+ (Y_ij.^2) ) .+ (Z_ij.^2)
                r_ij = sqrt.(r_p2_ij)
                newtonianNb_Potential.x[:,j] .+= (μ[i]./r_ij)
            end #if i != j
        end #for, i
        acc_eph.x[:,3j-2] .= differentiate.(ss16asteph.x[:,3(Nm1+j)-2])
        acc_eph.x[:,3j-1] .= differentiate.(ss16asteph.x[:,3(Nm1+j)-1])
        acc_eph.x[:,3j  ] .= differentiate.(ss16asteph.x[:,3(Nm1+j)  ])
    end #for, j
    return ss16asteph, acc_eph, newtonianNb_Potential
end

function save2jldandcheck(objname, sol)
    outfilename = string(objname, "_jt.jld")
    return __save2jldandcheck(outfilename, sol)
end

function __save2jldandcheck(outfilename, sol)
    println("Saving solution to file: $outfilename")
    jldopen(outfilename, "w") do file
        #loop over solution variables
        for ind in eachindex(sol)
            varname = string(ind)
            println("Saving variable: ", varname)
            write(file, varname, sol[ind])
        end
    end
    println("Checking that all variables were saved correctly...")
    #loop over solution variables
    for ind in eachindex(sol)
        varname = string(ind)
        #read varname from files and assign recovered variable to recovered_sol_i
        recovered_sol_i = JLD.load(outfilename, varname)
        #check that varname was recovered succesfully
        @show recovered_sol_i == sol[ind]
    end
    println("Saved solution")
    return outfilename
end

function taylor_minimum(pol::Taylor1{T}, x0::T; niters::Int=10) where {T<:Real}
    dpol = differentiate(pol)
    dpol2 = differentiate(dpol)
    xnewton::T = x0
    #@show xnewton
    for i in 1:niters
        xnewton -= dpol(xnewton)/dpol2(xnewton)
        #@show xnewton, dpol(xnewton)
    end
    return xnewton
end

function taylor_roots(pol::Taylor1{T}, x0::T; niters::Int=10) where {T<:Real}
    dpol = differentiate(pol)
    xnewton::T = x0
    #@show xnewton
    for i in 1:niters
        xnewton -= pol(xnewton)/dpol(xnewton)
        #@show xnewton, pol(xnewton)
    end
    return xnewton
end

function scaling(a::Taylor1{Taylor1{T}}, c::T) where {T<:Real}
    x = c*Taylor1( Taylor1(a.order).coeffs*one(a[0]) )
    return a(x)
end

function propagate(objname::String, dynamics::Function, maxsteps::Int, jd0::T,
        tspan::T, ephfile::String; kwargs...) where {T<:Real}
    # load ephemeris
    ss16asteph_et = JLD.load(ephfile, "ss16ast_eph")
    propagate(objname::String, dynamics::Function, maxsteps::Int, jd0::T, tspan::T, ss16asteph_et; kwargs...)
end

function propagate(objname::String, dynamics::Function, maxsteps::Int, jd0::T,
        tspan::T, ss16asteph_et::TaylorInterpolant; output::Bool=true, newtoniter::Int=10,
        dense::Bool=false, q0::Vector=initialcond(), radarobsfile::String="",
        opticalobsfile::String="", quadmath::Bool=false,
        debias_table::String="2018", μ_ast::Vector=μ_ast343_DE430[1:end],
        lyap::Bool=false, order::Int=order, abstol::T=abstol, tord::Int=10,
        niter::Int=5) where {T<:Real}
    # get asteroid initial conditions
    @assert length(q0) == 8
    @show q0
    @show jd0, jd0-JD_J2000
    # Number of bodies
    Nm1 = (size(ss16asteph_et.x)[2]-13) ÷ 6 # number of massive bodies
    @show Nm1
    N = Nm1 + 1 # number of bodies, including NEA
    # vector of G*m values
    μ = vcat(μ_DE430[1:11], μ_ast[1:Nm1-11], zero(μ_DE430[1]))

    # check: number of SS bodies (N) in ephemeris must be equal to length of GM vector (μ)
    @assert N == length(μ) "Total number of bodies in ephemeris must be equal to length of GM vector μ"

    # process ephemeris (switch from km,km/s units to au,au/day)
    # compute Newtonian accelerations and potentials (used in post-Newtonian accelerations)
    ss16asteph_auday, acc_eph, newtonianNb_Potential = loadeph(ss16asteph_et, μ)

    # interaction matrix with flattened bodies
    UJ_interaction = fill(false, N)
    # UJ_interaction[su] = true
    UJ_interaction[ea] = true
    # UJ_interaction[mo] = true
    params = (ss16asteph_auday, acc_eph, newtonianNb_Potential, jd0, UJ_interaction, N, μ)

    if quadmath
        _q0 = one(Float128)*q0
        _t0 = zero(Float128)
        _abstol = Float128(abstol)
        _ss16asteph = TaylorInterpolant(Float128(ss16asteph_auday.t0), Float128.(ss16asteph_auday.t), map(x->Taylor1(Float128.(x.coeffs)), ss16asteph_auday.x))
        _acc_eph = TaylorInterpolant(Float128(acc_eph.t0), Float128.(acc_eph.t), map(x->Taylor1(Float128.(x.coeffs)), acc_eph.x))
        _newtonianNb_Potential = TaylorInterpolant(Float128(newtonianNb_Potential.t0), Float128.(newtonianNb_Potential.t), map(x->Taylor1(Float128.(x.coeffs)), newtonianNb_Potential.x))
        _params = (_ss16asteph, _acc_eph, _newtonianNb_Potential, Float128(jd0), UJ_interaction, N, μ)
    else
        _q0 = q0
        _t0 = zero(Float64)
        _abstol = abstol
        _params = params
    end

    @show _tmax = _t0+tspan*yr #final time of integration

    # propagate orbit
    if dense
        @time interp = apophisinteg(dynamics, _q0, _t0, _tmax, order, _abstol, _params; maxsteps=maxsteps, dense=dense)
        if quadmath
            asteph_t0 = Float64(jd0-JD_J2000) # days since J2000 until initial integration time
            asteph_t = Float64.(interp.t[:])
            asteph_x = convert(Array{Taylor1{eltype(q0)}}, interp.x[:,:])
            asteph = TaylorInterpolant(asteph_t0, asteph_t, asteph_x)
            sol = (asteph=asteph,
            )
        else
            asteph_t0 = (jd0-JD_J2000) # days since J2000 until initial integration time
            asteph_t = interp.t[:]
            asteph_x = interp.x[:,:]
            asteph = TaylorInterpolant(asteph_t0, asteph_t, asteph_x)
            sol = (asteph=asteph,
            )
        end
    elseif lyap
        @time sol_objs = lyap_apophisinteg(dynamics, _q0, _t0, _tmax, order, _abstol, _params; maxsteps=maxsteps)
        sol = (
            tv=convert(Array{eltype(q0)}, sol_objs[1][:]),
            xv=convert(Array{eltype(q0)}, sol_objs[2][:,:]),
            λv=convert(Array{eltype(q0)}, sol_objs[3][:,:])
        )
    else
        @time sol_objs = apophisinteg(dynamics, rvelea, _q0, _t0, _tmax, order, _abstol, _params; maxsteps=maxsteps, newtoniter=newtoniter, dense=true)
        asteph_t0 = (_params[4]-JD_J2000) # days since J2000 until initial integration time
        asteph_t = sol_objs[1].t[:]
        asteph_x = sol_objs[1].x[:,:]
        asteph = TaylorInterpolant(asteph_t0, asteph_t, asteph_x)
        sol = (
            asteph=asteph,
            tv = asteph_t0 .+ asteph_t,
            xv = asteph_x(),
            tvS1=convert(Array{eltype(q0)}, sol_objs[2][:]),
            xvS1=convert(Array{eltype(q0)}, sol_objs[3][:,:]),
            gvS1=convert(Array{eltype(q0)}, sol_objs[4][:])
        )
    end

    #write solution and predicted values of observations (if requested) to .jld files
    if output
        outfilename = save2jldandcheck(objname, sol)
        # if requested by user, calculate computed (i.e., predicted) values of observations
        if !lyap # don't compute observation ephemeris when "in Lyapunov spectra mode"
            try
                furnsh(
                    joinpath(artifact"naif0012", "naif0012.tls"), # load leapseconds kernel
                    joinpath(artifact"de430", "de430_1850-2150.bsp"), # at least one SPK file must be loaded to read .tls file
                )
                if radarobsfile != ""
                    println("compute_radar_obs")
                    @time compute_radar_obs("deldop_"*basename(radarobsfile)*".jld", radarobsfile, asteph, ss16asteph_et, tord=tord, niter=niter)
                end
                if opticalobsfile != ""
                    println("compute_optical_obs")
                    @time compute_optical_obs("radec_"*basename(opticalobsfile)*".jld", opticalobsfile, asteph, ss16asteph_et, debias_table=debias_table, niter=niter)
                end
            catch e
                @error "Unable to compute observation residuals" exception=(e, catch_backtrace())
            end
        end
    end

    return sol
end

function compute_radar_obs(outfilename::String, radarobsfile::String,
        asteph::TaylorInterpolant, ss16asteph::TaylorInterpolant; tc::Real=1.0,
        autodiff::Bool=true, tord::Int=10, niter::Int=5)
    @assert isfile(radarobsfile) "Cannot open file: $radarobsfile"
    Nm1 = (size(ss16asteph.x)[2]-13) ÷ 6
    N = Nm1 + 1
    # TODO: check that first and last observation times are within interpolation interval
    # asteroid/small-body
    # Change t, x, v units, resp., from days, au, au/day to sec, km, km/sec
    asteph_ord = asteph.x[1].order
    asteph_t0 = asteph.t0*daysec
    asteph_t = asteph.t*daysec
    asteph_r = au*map(x->x(Taylor1(asteph_ord)/daysec), asteph.x[:,1:3])
    asteph_v = (au/daysec)*map(x->x(Taylor1(asteph_ord)/daysec), asteph.x[:,4:6])
    asteph_x = hcat(asteph_r, asteph_v)
    asteph_et = TaylorInterpolant(asteph_t0, asteph_t, asteph_x)
    # Sun (su=1)
    # Change x, v units, resp., from au, au/day to km, km/sec
    sseph_t0 = ss16asteph.t0
    sseph_t = ss16asteph.t
    sseph_x = ss16asteph.x
    sun_r = au*sseph_x[:,nbodyind(Nm1,su)[1:3]]
    sun_v = (au/daysec)*sseph_x[:,nbodyind(Nm1,su)[4:6]]
    sun_x = hcat(sun_r, sun_v)
    sun_et = TaylorInterpolant(sseph_t0, sseph_t, sun_x)
    # Earth (ea=4)
    # Change x, v units, resp., from au, au/day to km, km/sec
    earth_r = au*sseph_x[:,nbodyind(Nm1,ea)[1:3]]
    earth_v = (au/daysec)*sseph_x[:,nbodyind(Nm1,ea)[4:6]]
    earth_x = hcat(earth_r, earth_v)
    earth_et = TaylorInterpolant(sseph_t0, sseph_t, earth_x)
    # construct DataFrame with delay/doppler data from JPL radar obs file, including delay/doppler ephemeris (i.e., predicted values)
    deldop_table_jld = delay_doppler(radarobsfile, niter, xve=earth_et, xvs=sun_et, xva=asteph_et, tc=tc, tord=tord, autodiff=autodiff)
    #save data to file
    println("Saving data to file: $outfilename")
    jldopen(outfilename, "w") do file
        addrequire(file, DataFrames)
        addrequire(file, TaylorSeries)
        # write variables to jld file
        JLD.write(file, "deldop_table", deldop_table_jld)
    end
    return nothing
end

function compute_optical_obs(outfilename::String, opticalobsfile::String,
        asteph, ss16asteph; debias_table::String="2018", niter::Int=5)
    @assert isfile(opticalobsfile) "Cannot open file: $opticalobsfile"
    Nm1 = (size(ss16asteph.x)[2]-13) ÷ 6
    N = Nm1 + 1
    # TODO: check that first and last observation times are within interpolation interval
    function asteph_et(et)
        return auday2kmsec(asteph(et/daysec)[1:6])
    end
    function earth_et(et)
        return auday2kmsec(ss16asteph(et)[union(3*4-2:3*4,3*(N-1+4)-2:3*(N-1+4))])
    end
    function sun_et(et)
        return auday2kmsec(ss16asteph(et)[union(3*1-2:3*1,3*(N-1+1)-2:3*(N-1+1))])
    end
    # construct DataFrame with ra/dec data from MPC optical obs file, including ra/dec ephemeris (i.e., predicted values)
    radec_table_jdb = radec_table(opticalobsfile, niter, xve=earth_et, xvs=sun_et, xva=asteph_et, debias_table=debias_table)
    #save data to file
    println("Saving data to file: $outfilename")
    jldopen(outfilename, "w") do file
        addrequire(file, DataFrames)
        addrequire(file, TaylorSeries)
        # write variables to jld file
        JLD.write(file, "radec_table", radec_table_jdb)
    end
    return nothing
end


text:
using HCIToolbox
using PlotUtils: zscale
using PyCall
using SiriusB

fits = pyimport("astropy.io.fits")
pro = pyimport("proplot")

# rcParams
pro.rc["image.origin"] = "lower"
pro.rc["image.cmap"] = "inferno"
pro.rc["grid"] = false

cube= fits.getdata(datadir("epoch_2020nov21", "processed", "2020nov21_sirius-b_cube_calib_registered_crop.fits"))

flat = collapse(cube)

function norm_psf(frame)
    tmp = frame .- minimum(frame) .+ 1
    return tmp ./ maximum(tmp)
end
psf = norm_psf(flat)

tick_locs = range(25, 174, length=5)
tick_labs = @. string(round(SiriusB.auscale * (tick_locs - 99.5), digits=1))

py"""
import proplot as pro
from matplotlib.colors import LogNorm

fig, axs = pro.subplots(figwidth="4in")

axs.imshow($psf,
    norm=LogNorm(vmin=3e-3, vmax=5e-2),
    colorbar="r",
)
# axs.format(title="Epoch 2020-11-21")
axs.format(
    xticks=$tick_locs,
    xticklabels=$tick_labs,
    yticks=$tick_locs,
    yticklabels=$tick_labs,
    xlabel="x [AU]",
    ylabel="y [AU]",
)
fig.save($(figuredir("psf.pdf")))
"""

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
# NOTE: Bolt's background functions are in terms of x ≡ ln(a), the log scale factor

# derived quantities (I've chosen natural units, possibly the wrong choice)
const km_s_Mpc_100 = ustrip(natural(100.0u"km/s/Mpc"))  # [eV]
const G_natural = ustrip(natural(float(NewtonianConstantOfGravitation)))

H₀(par::AbstractCosmoParams) = par.h * km_s_Mpc_100
ρ_crit(par::AbstractCosmoParams) = (3 / 8π) * H₀(par)^2 / G_natural  # [eV⁴]
Ω_Λ(par::AbstractCosmoParams) = 1 - (par.Ω_r*(1+(7par.N_ν/8)*(4/11)^4/3)
                                     + par.Ω_b + par.Ω_m)  # dark energy density

#need to account for neutrinos

# Hubble parameter ȧ/a in Friedmann background
H_a(a, par::AbstractCosmoParams) = H₀(par) * √((par.Ω_m + par.Ω_b) * a^(-3)
                                                + par.Ω_r*(1+(7par.N_ν/8)*(4/11)^4/3) * a^(-4)
                                                + Ω_Λ(par))
# conformal time Hubble parameter, aH
ℋ_a(a, par::AbstractCosmoParams) = a * H_a(a, par)

# functions in terms of x
H(x, par::AbstractCosmoParams) = H_a(x2a(x),par)
ℋ(x, par::AbstractCosmoParams) = ℋ_a(x2a(x), par)

# conformal time
function η(x, par::AbstractCosmoParams)
    return quadgk(a -> 1.0 / (a * ℋ_a(a, par)), 0.0, x2a(x))[1]
end

# now build a Background with these functions

# a background is parametrized on the scalar type T, the interpolator type IT,
# and a type for the grid GT
abstract type AbstractBackground{T, IT<:AbstractInterpolation{T,1}, GT} end

struct Background{T, IT, GT} <: AbstractBackground{T, IT, GT}
    H₀::T
    η₀::T
    ρ_crit::T
    Ω_Λ::T

    x_grid::GT
    ℋ::IT
    ℋ′::IT
    ℋ′′::IT
    η::IT
    η′::IT
    η′′::IT
end

function Background(par::AbstractCosmoParams{T}; x_grid=-20.0:0.01:0.0) where T
    ℋ_ = spline(x_grid, [ℋ(x, par) for x in x_grid])
    η_ = spline(x_grid, [η(x, par) for x in x_grid])

    return Background(
        T(H₀(par)),
        T(η(0.0, par)),
        T(ρ_crit(par)),
        T(Ω_Λ(par)),
        x_grid,

        ℋ_,
        spline_∂ₓ(ℋ_, x_grid),
        spline_∂ₓ²(ℋ_, x_grid),

        η_,
        spline_∂ₓ(η_, x_grid),
        spline_∂ₓ²(η_, x_grid),
    )
end


text:
module TestParametricBootstrap

using Bootstrap
using Test

using DataFrames
using Distributions
using GLM
using RDatasets
using StatsBase

const city = dataset("boot", "city")
const aircondit = convert(Vector{Float64}, dataset("boot", "aircondit")[!, :Hours])

@testset "Parametric bootstraps" begin

    function test_bootsample(bs, ref, raw_data, n)

        show(IOBuffer(), bs)
        @test typeof(bs) <: ParametricBootstrapSample
        t0 = original(bs)
        @test length(t0) == length(ref)
        [@test t ≈ r for (t, r) in zip(t0, ref)]

        t1 = straps(bs)
        @test length(t1) == length(t0)
        [@test length(t) == n for t in t1]
        [@test (minimum(t) <= tr && maximum(t) >= tr) for (t, tr) in zip(t1, t0)]
        [@test eltype(t) == eltype(tr) for (t, tr) in zip(t1, t0)]

        @test Bootstrap.data(bs) == raw_data

        @test nrun(sampling(bs)) == n
        @test nrun(bs) == n

        @test length(bias(bs)) == length(ref)
        [@test (b[1] > -Inf && b[1] < Inf) for b in bias(bs)]
        @test length(stderror(bs)) == length(ref)
        [@test s >= 0 for s in stderror(bs)]

        [@test original(bs, i) == original(bs)[i]  for i in 1:nvar(bs)]
        [@test straps(bs, i) == straps(bs)[i]  for i in 1:nvar(bs)]
        [@test bias(bs, i) == bias(bs)[i]  for i in 1:nvar(bs)]
        [@test stderror(bs, i) == stderror(bs)[i]  for i in 1:nvar(bs)]

        m = model(bs)
        @test typeof(m) <: Bootstrap.Model

        return Nothing
    end

    function test_confint(bs)

        cim_all = (BasicConfInt(), PercentileConfInt(), NormalConfInt(), BCaConfInt())
        for cim in cim_all
            c = confint(bs, cim)
            [@test (x[1] >= x[2]  && x[1] <= x[3]) for x in c]
            [@test x[1] ≈ t0 for (x, t0) in zip(c, original(bs))]
        end

        return Nothing
    end

    n = 100

    @testset "Parametric distributions" begin

        r = randn(20)
        ref = mean(r)

        @testset "Basic resampling: Normal distribution" begin
            bs = bootstrap(mean, r, Model(Normal), BasicSampling(n))
            test_bootsample(bs, ref, r, n)
            test_confint(bs)
        end

        @testset "Balanced resampling: Normal distribution" begin
            bs = bootstrap(mean, r, Model(Normal), BalancedSampling(n))
            test_bootsample(bs, ref, r, n)
            test_confint(bs)
        end

        ref = mean(fit(Exponential, aircondit))

        @testset "Basic resampling: Exponential distribution" begin
            bs = bootstrap(mean, aircondit, Model(Exponential), BasicSampling(n))
            test_bootsample(bs, ref, aircondit, n)
            test_confint(bs)
        end

        @testset "Balanced resampling: Exponential distribution" begin
            bs = bootstrap(mean, aircondit, Model(Exponential), BalancedSampling(n))
            test_bootsample(bs, ref, aircondit, n)
            test_confint(bs)
        end

    end

    @testset "Linear regression models" begin

        ref = coef(fit(LinearModel, @formula(log10(X) ~ log10(U)), city))

        @testset "Residual resampling" begin
            bs = bootstrap(coef, city, Model(LinearModel, @formula(log10(X) ~ log10(U))), ResidualSampling(n))
            test_bootsample(bs, ref, city, n)
        end

        @testset "Wild resampling: Rademacher" begin
            bs = bootstrap(coef, city, Model(LinearModel, @formula(log10(X) ~ log10(U))), WildSampling(n, rademacher))
            test_bootsample(bs, ref, city, n)
            @test isa(noise(sampling(bs)), Function)
        end

        @testset "Wild resampling: Mammen" begin
            bs = bootstrap(coef, city, Model(LinearModel, @formula(log10(X) ~ log10(U))), WildSampling(n, mammen))
            test_bootsample(bs, ref, city, n)
            @test isa(noise(sampling(bs)), Function)
        end

    end

    @testset "Generalized linear regression models" begin

        ref = coef(fit(GeneralizedLinearModel, @formula(log10(X) ~ log10(U)), city, Normal()))
        maxiter = 200
        atol = 1e-3

        @testset "Residual resampling" begin
            bs = bootstrap(coef, city,
                           Model(GeneralizedLinearModel, @formula(log10(X) ~ log10(U)), Normal(), maxiter = maxiter, atol = atol),
                           ResidualSampling(n))
            test_bootsample(bs, ref, city, n)
        end

        @testset "Residual resampling with link function" begin
            bs = bootstrap(coef, city,
                           Model(GeneralizedLinearModel, @formula(log10(X) ~ log10(U)), Normal(), IdentityLink(), maxiter = maxiter, atol = atol),
                           ResidualSampling(n))
            test_bootsample(bs, ref, city, n)
        end

        @testset "Wild resampling: Rademacher" begin
            bs = bootstrap(coef, city,
                           Model(GeneralizedLinearModel, @formula(log10(X) ~ log10(U)), Normal(), maxiter = maxiter, atol = atol),
                           WildSampling(n, rademacher))
            test_bootsample(bs, ref, city, n)
            @test isa(noise(sampling(bs)), Function)
        end

        @testset "Wild resampling with link function: Mammen" begin
            bs = bootstrap(coef, city,
                           Model(GeneralizedLinearModel, @formula(log10(X) ~ log10(U)), Normal(), IdentityLink(), maxiter = maxiter, atol = atol),
                           WildSampling(n, mammen))
            test_bootsample(bs, ref, city, n)
            @test isa(noise(sampling(bs)), Function)
        end

    end

end

end


text:
mutable struct Tape
    tape
    pos
    relbase
end

function Tape(tape)
    Tape(tape, 0, 0)
end

function Base.copy(tape)
    Tape(copy(tape.tape), tape.pos, tape.relbase)
end

function parse_data(program_str)
    [parse(Int64, v) for v in split(program_str, ',')]
end

function read(tape::Tape, pos)
    return tape.tape[pos + 1]
end

function write!(tape::Tape, pos, value)
    tape.tape[pos + 1] = value
end

POSITION_MODE = 0
IMMEDIATE_MODE = 1
RELATIVE_MODE = 2

function parse_opcode(instruction)
    instruction, opcode = divrem(instruction, 100)
    param_modes = []  # TODO memory assignment...
    while instruction > 0
        instruction, mode = divrem(instruction, 10)
        push!(param_modes, mode)
    end
    param_modes, opcode
end

function get_op_mode(tape::Tape, offset)
    instruction = read(tape, tape.pos)
    param_modes, opcode = parse_opcode(instruction)
    if offset <= length(param_modes)
        param_modes[offset]
    else
        POSITION_MODE
    end
end

function get_arg(tape::Tape, offset)
    mode = get_op_mode(tape, offset)
    if mode == POSITION_MODE
        pos_arg = read(tape, tape.pos + offset)
        arg = read(tape, pos_arg)
    elseif mode == IMMEDIATE_MODE
        arg = read(tape, tape.pos + offset)
    elseif mode == RELATIVE_MODE
        pos_arg = read(tape, tape.pos + offset)
        arg = read(tape, tape.relbase + pos_arg)
    else
        error("invalid mode: " * mode * " !")
    end
    arg
end

function set_arg!(tape::Tape, offset, value)
    mode = get_op_mode(tape, offset)
    pos = read(tape, tape.pos + offset)
    if mode == RELATIVE_MODE
        pos += tape.relbase
    end
    write!(tape, pos, value)
end

function op_add!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    result = arg1 + arg2
    set_arg!(tape, 3, result)
    tape.pos += 4
end

function op_mul!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    result = arg1 * arg2
    set_arg!(tape, 3, result)
    tape.pos += 4
end

function op_input!(tape)
    print("INPUT: ")
    input = parse(Int64, readline())
    set_arg!(tape, 1, input)
    tape.pos += 2
end

function op_output!(tape)
    output = get_arg(tape, 1)
    println("OUTPUT: ", output)
    tape.pos += 2
end

function op_jumptrue!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    if arg1 != 0
        tape.pos = arg2
    else
        tape.pos += 3
    end
end

function op_jumpfalse!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    if arg1 == 0
        tape.pos = arg2
    else
        tape.pos += 3
    end
end

function op_lessthan!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    set_arg!(tape, 3, arg1 < arg2 ? 1 : 0)
    tape.pos += 4
end

function op_equal!(tape)
    arg1 = get_arg(tape, 1)
    arg2 = get_arg(tape, 2)
    set_arg!(tape, 3, arg1 == arg2 ? 1 : 0)
    tape.pos += 4
end

function op_adjustrelbase!(tape)
    arg = get_arg(tape, 1)
    tape.relbase += arg
    tape.pos += 2
end

function run_opcode!(tape)
    instruction = read(tape, tape.pos)
    param_modes, opcode = parse_opcode(instruction)
    if opcode == 1
        op_add!(tape)
    elseif opcode == 2
        op_mul!(tape)
    elseif opcode == 3
        op_input!(tape)
    elseif opcode == 4
        op_output!(tape)
    elseif opcode == 5
        op_jumptrue!(tape)
    elseif opcode == 6
        op_jumpfalse!(tape)
    elseif opcode == 7
        op_lessthan!(tape)
    elseif opcode == 8
        op_equal!(tape)
    elseif opcode == 9
        op_adjustrelbase!(tape)
    elseif opcode == 99
        return opcode
    else
        error("unknown opcode " * string(opcode) * "!")
    end
    opcode
end

function run_program!(tape)
    while true
        opcode = run_opcode!(tape)
        if opcode == 99
            return tape
        end
    end
end

function set_memory!(tape, size)
    memory = zeros(Int64, size)
    memory[1:length(tape.tape)] = tape.tape
    Tape(memory, tape.pos, tape.relbase)
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

function _uvbybeta(by::T, m1::T, c1::T, hbeta::T, eby_in::T,
                   n::Integer) where {T<:AbstractFloat}
    # Rm1 = -0.33 & Rc1 = 0.19 & Rub = 1.53
    if n < 1 || n > 8
      error("Input should be an integer in the range 1:8, giving approximate
            stellar classification")
    end
    ub =  c1 + 2 * (m1 + by)
    # For group 1, beta is a luminosity indicator, c0 is a temperature indicator.
    # (u-b) is also a suitable temperature indicator.
    if n == 1
        # For dereddening, linear relation used between the intrinsic (b-y)
        # and (u-b) (Crawford 1978, AJ 83, 48)
        if isnan(eby_in)
           eby_in::T = (13.608 * by - ub + 1.467) / (12.078)
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)
        # When beta is not given, it is estimated using a cubic fit to the c0-beta
        # relation for luminosity class V given in Crawford (1978).
        if isnan(hbeta)
            hbeta::T = @evalpoly c0 2.61033 0.132557 0.161463 -0.027352
        end
        # Calculation of the absolute magnitude by applying the calibration of
        # Balona & Shobbrock (1974, MNRAS 211, 375)
        g = log10(hbeta - 2.515) - 1.6 * log10(c0 + 0.322)
        mv = 3.4994 + 7.2026 * log10(hbeta - 2.515) + @evalpoly(g, 0, -2.3192, 0, 2.9375)
        te = 5040 / (0.2917 * c0 + 0.2)
        # The ZAMS value of m0 is calculated from a fit to the data of Crawford (1978),
        # modified by Hilditch, Hill & Barnes (1983, MNRAS 204, 241)
        delm0 = @evalpoly(c0, 0.07473, 0.109804, -0.139003, 0.0957758) - m0
    elseif n == 2
        # For dereddening the linear relations between c0 and (u-b) determined from
        # Zhang (1983, AJ 88, 825) is used.
        if isnan(eby_in)
            eby_in = ((1.5 * c1 - ub + 0.035) / (1.5/(1.53/0.19) - 1)) / 1.53
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)

        if isnan(hbeta)
            hbeta = 2.542 + 0.037 * c0
        end
        delm0 = NaN
    elseif n == 3
        # For dereddening the linear relations between c0 and (u-b) determined from
        # Zhang (1983, AJ 88, 825) is used.
        if isnan(eby_in)
            eby_in = ((1.36 * c1 - ub + 0.004) / (1.36/(1.53/0.19) - 1)) / 1.53
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)
        # When beta is not given, it is derived from a fit of the c0-beta
        # relation of Zhang (1983).
        if isnan(hbeta)
            hbeta = 2.578 + 0.047 * c0
        end
        delm0 = NaN
    elseif n == 4
        # For dereddening the linear relations between c0 and (u-b) determined from
        # Zhang (1983, AJ 88, 825) is used.
        if isnan(eby_in)
            eby_in = ((1.32 * c1 - ub - 0.056) / (1.32/(1.53/0.19) - 1)) / 1.53
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)
        # When beta is not given, it is derived from a fit of the c0-beta
        # relation of Zhang (1983).
        if isnan(hbeta)
            hbeta = 2.59 + 0.066 * c0
        end
        delm0 = NaN
    elseif n == 5
        # For group 5, the hydrogen Balmer lines are at maximum; hence two new
        # parameters, a0 = f{(b-y),(u-b)} and r = f{beta,[c1]} are defined in
        # order to calculate absolute magnitude and metallicity.
        if isnan(eby_in)
            by0 = @evalpoly(m1 + 0.33 * by, -0.0235, -0.53921, 4.2608)
            while true
                bycorr = by0
                by0 = @evalpoly(m1 + 0.33 * (by - bycorr), 0.175709, -3.36225, 14.0881)

                if abs(bycorr - by0) < 0.001
                    break
                end
            end
            eby_in = by - by0
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)

        if isnan(hbeta)
            hbeta = 2.7905 - 0.6105 * by + 0.5 * m0 + 0.0355 * c0
        end
        a0 = by0 + 0.18 * (ub0 - 1.36)
        # MV is calculated according to Stroemgren (1966, ARA&A 4, 433)
        # with corrections by Moon & Dworetsky (1984, Observatory 104, 273)
        mv = 1.5 + 6 * a0 - 17 * (0.35 * (c1 - 0.19 * by) - (hbeta - 2.565))
        te = 5040 / (0.7536 * a0 + 0.5282)
        delm0 = @evalpoly(by0, 0.1598, 0.86888, -3.95105) - m0
    elseif n == 6

        if isnan(hbeta)
            hbeta = 3.06 - 1.221 * by - 0.104 * c1
        end
        m1zams = @evalpoly(hbeta, -17.209, 12.26, -2.158)

        if hbeta <= 2.74
            c1zams = 3 * hbeta - 7.56
            mvzams = 22.14 - 7 * hbeta
        elseif 2.74 < hbeta <= 2.82
            c1zams = 2 * hbeta - 4.82
            mvzams = 11.16 - 3 * hbeta
        else
            c1zams = 2 * hbeta - 4.83
            mvzams = @evalpoly hbeta -696.41 497.2 -88.4
        end

        if isnan(eby_in)
            delm1 = m1zams - m1
            delc1 = c1 - c1zams

            if delm1 < 0
                by0 = 2.946 - hbeta - 0.1 * delc1 - 0.25 * delm1
            else
                by0 = 2.946 - hbeta - 0.1 * delc1
            end
            eby_in = by - by0
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)
        delm0 =  m1zams - m0
        mv = mvzams - 9 * (c0 - c1zams)
        te = 5040 / (0.771453 * by0 + 0.546544)
    elseif n == 7
        # For group 7 c1 is the luminosity indicator for a particular beta, while
        # beta {or (b-y)0} indicates temperature.

        # Where beta is not available iteration is necessary to evaluate a corrected
        # (b-y) from which beta is then estimated.
        if isnan(hbeta)
            byinit = by
            m1init = m1
            for i in 1:1000
                m1by = @evalpoly byinit 0.345 -1.32 2.5
                bycorr = byinit + (m1by - m1init) / 2

                if abs(bycorr - byinit) <= 0.0001
                    break
                end
                byinit = bycorr
                m1init = m1by
            end
            hbeta = @evalpoly byinit 2.96618 -1.32861 1.01425
        end
        # m1(ZAMS) and mv(ZAMS) are calculated according to Crawford (1975) with
        # corrections suggested by Hilditch, Hill & Barnes (1983, MNRAS 204, 241)
        # and Olson (1984, A&AS 57, 443).
        m1zams = @evalpoly hbeta 46.4167 -34.4538 6.41701
        mvzams = @evalpoly hbeta 324.482 -188.748 11.0494 5.48012
        # c1(ZAMS) calculated according to Crawford (1975)
        if hbeta <= 2.65
            c1zams = 2 * hbeta - 4.91
        else
            c1zams = @evalpoly hbeta 72.879 -56.9164 11.1555
        end

        if isnan(eby_in)
            dbeta = 2.72 - hbeta
            eby_in = by - (0.222 - 0.05 * (c1 - c1zams) +
                    (1.11 - (0.1 + 3.6 * (m1zams - m1))) * dbeta + 2.7 * dbeta^2)
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)
        delm0 = m1zams - m0
        mv = mvzams - (9 + 20 * dbeta) * (c0 - c1zams)
        te = 5040/(0.771453 * by0 + 0.546544)
    elseif n == 8
        # Dereddening is done using color-color relations derived from
        # Olson (1984, A&AS 57, 443)
        if by <= 0.65
            eby_in = (5.8651 * by - ub - 0.8975) / 4.3351
        elseif 0.65 < by < 0.79
            eby_in = (-0.7875 * by - c1 + 0.6585) / -0.9775
            by0 = by - eby_in

            if by0 < 0.65
                eby_in = (5.8651 * by - ub - 0.8975) / 4.3351
            end
        else
            eby_in = (0.5126 * by - c1 - 0.3645) / 0.3226
            by0 = by - eby_in

            if by0 < 0.79
                eby_in = (-0.7875 * by - c1 + 0.6585) / -0.9775
            end
            by0 = by - eby_in
        end
        by0, m0, c0, ub0 = deredd(eby_in, by, m1, c1, ub)
        # m1(ZAMS), c1(ZAMS), and MV(ZAMS) are calculated according to Olson (1984)
        m1zams = @evalpoly by0 7.18436 -49.43695 122.1875 -122.466 42.93678

        if by0 < 0.65
            c1zams = @evalpoly by0 3.78514 -21.278 42.7486 -28.7056
            mvzams = @evalpoly by0 -59.2095 432.156 -1101.257 1272.503 -552.48
        elseif 0.65 <= by0 < 0.79
            c1zams = @evalpoly by0 0.33657 0.116031 -0.631821
            mvzams = @evalpoly by0 3.4305 4.97911 1.37632
        else
            c1zams = @evalpoly by0 -0.37237 0.530426 -0.010028
            mvzams = @evalpoly by0 4.37507 3.92776 1.18298
        end
        delm0 = m1zams - m0
        # Teff and MV calibration of Olson (1984)
        if by0 < 0.505
            f = 10 - 80 * (by0 - 0.38)
            te = exp10(3.924 - 0.416 * by0)
        else
            f = zero(T)
            te = exp10(3.869 -0.341 * by0)
        end
        mv = mvzams + f * (c1zams - c0) + 3.2 * delm0 - 0.07
    end

    if 1 < n < 5
        # c0-beta relation for ZAMS stars according to Crawford (1978,AJ 83, 48),
        # modified by Hilditch, Hill & Barnes (1983, MNRAS 204, 241).
        b_betaza = @evalpoly(c0, 2.62745, 0.228638, -0.099623, 0.277363, -0.160402) - 2.5
        # MV(ZAMS) calculated according to Balona & Shobbrock (1984, MNRAS 211, 375)
        mvzams = @evalpoly b_betaza -9.563 77.18 -206.98 203.704
        # MV is calculated from the d(beta)-d(MV) relation of Zhang (1983)
        dbeta = b_betaza - hbeta + 2.5
        mv = mvzams - @evalpoly(dbeta, 0.08, 61, -121.6)
        # Estimate of Teff by coupling the relations of Boehm-Vitense
        # (1981, ARA&A 19, 295) and Zhang (1983)
        te = 5040 / (0.27346 + 0.35866 * ub0)
    end
    # Transformation according to the FV-(b-y)0 relation of Moon (1984, MNRAS 211, 21P)
    if by0 <= 0.335
        fv = @evalpoly by0 3.981 -1.092 3.731 -6.759
    else
        fv = 3.959 - 0.534 * by0
    end
    radius = exp10(2 * (4.236 - 0.1 * mv - fv))
    eby = eby_in
    return te, mv, eby, delm0, radius
end

"""
    uvbybeta(by, m1, c1, n[, hbeta=NaN, eby_in=NaN]) -> te, mv, eby, delm0, radius

### Purpose ###

Derive dereddened colors, metallicity, and Teff from Stromgren colors.

### Arguments ###

* `by`: Stromgren b-y color
* `m1`: Stromgren line-blanketing parameter
* `c1`: Stromgren Balmer discontinuity parameter
* `n`: Integer which can be any value between 1 to 8, giving approximate stellar
  classification.
  (1) B0 - A0, classes III - V, 2.59 < Hbeta < 2.88,-0.20 <   c0   < 1.00
  (2) B0 - A0, class   Ia     , 2.52 < Hbeta < 2.59,-0.15 <   c0   < 0.40
  (3) B0 - A0, class   Ib     , 2.56 < Hbeta < 2.61,-0.10 <   c0   < 0.50
  (4) B0 - A0, class   II     , 2.58 < Hbeta < 2.63,-0.10 <   c0   < 0.10
  (5) A0 - A3, classes III - V, 2.87 < Hbeta < 2.93,-0.01 < (b-y)o < 0.06
  (6) A3 - F0, classes III - V, 2.72 < Hbeta < 2.88, 0.05 < (b-y)o < 0.22
  (7) F1 - G2, classes III - V, 2.60 < Hbeta < 2.72, 0.22 < (b-y)o < 0.39
  (8) G2 - M2, classes  IV - V, 0.20 < m0    < 0.76, 0.39 < (b-y)o < 1.00
* `hbeta` (optional): H-beta line strength index. If it is not supplied, then by
  default its value will be `NaN` and the code will estimate a value based on by,
  m1,and c1. It is not used for stars in group 8.
* `eby_in` (optional): specifies the E(b-y) color to use. If not supplied, then by
  default its value will be `NaN` and E(b-y) will be estimated from the Stromgren colors.

### Output ###

* `te`: approximate effective temperature
* `mv`: absolute visible magnitude
* `eby`: color excess E(b-y)
* `delm0`: metallicity index, delta m0, may not be calculable for early B stars
  and so returns `NaN`.
* `radius`: stellar radius (R/R(solar))

### Example ###

Suppose 5 stars have the following Stromgren parameters

by = [-0.001 ,0.403, 0.244, 0.216, 0.394]
m1 = [0.105, -0.074, -0.053, 0.167, 0.186]
c1 = [0.647, 0.215, 0.051, 0.785, 0.362]
hbeta = [2.75, 2.552, 2.568, 2.743, 0]
nn = [1,2,3,7,8]

Determine the stellar parameters

```jldoctest
julia> using AstroLib

julia> by = [-0.001 ,0.403, 0.244, 0.216, 0.394];

julia> m1 = [0.105, -0.074, -0.053, 0.167, 0.186];

julia> c1 = [0.647, 0.215, 0.051, 0.785, 0.362];

julia> hbeta = [2.75, 2.552, 2.568, 2.743, 0];

julia> nn = [1,2,3,7,8];

julia> uvbybeta.(by, m1, c1, nn, hbeta)
5-element Array{NTuple{5,Float64},1}:
 (13057.535222326893, -0.27375469585031265, 0.04954396423248884, -0.008292894218734928, 2.7136529525371897)
 (14025.053834219656, -6.907050783073221, 0.4140562248995983, NaN, 73.50771722263974)
 (18423.76405400214, -5.935816553877892, 0.2828247876690783, NaN, 39.84106215808709)
 (7210.507090112837, 2.2180408083364167, 0.018404079180028038, 0.018750927360588615, 2.0459018065648165)
 (5755.671513413262, 3.9449408311022, -0.025062997393370458, 0.03241423718769865, 1.5339239690774464)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
uvbybeta(by::Real, m1::Real, c1::Real, n::Integer, hbeta::Real=NaN, eby_in::Real=NaN) =
    _uvbybeta(promote(float(by), float(m1), float(c1), float(hbeta), float(eby_in))..., n)


text:
"""
    make_optimization_model(
        model::MetabolicModel,
        optimizer;
        sense = MOI.MAX_SENSE,
    )

Convert `MetabolicModel`s to a JuMP model, place objectives and the equality
constraint.
"""
function make_optimization_model(model::MetabolicModel, optimizer; sense = MOI.MAX_SENSE)

    precache!(model)

    m, n = size(stoichiometry(model))
    xl, xu = bounds(model)

    optimization_model = Model(optimizer)
    @variable(optimization_model, x[i = 1:n])
    @objective(optimization_model, sense, objective(model)' * x)
    @constraint(optimization_model, mb, stoichiometry(model) * x .== balance(model)) # mass balance
    @constraint(optimization_model, lbs, xl .<= x) # lower bounds
    @constraint(optimization_model, ubs, x .<= xu) # upper bounds

    C = coupling(model) # empty if no coupling
    cl, cu = coupling_bounds(model)
    isempty(C) || @constraint(optimization_model, c_lbs, cl .<= coupling(model) * x) # coupling lower bounds
    isempty(C) || @constraint(optimization_model, c_ubs, coupling(model) * x .<= cu) # coupling upper bounds

    return optimization_model
end

"""
    is_solved(opt_model)

Return `true` if `opt_model` solved successfully (solution is optimal or locally
optimal).  Return `false` if any other termination status is reached.
Termination status is defined in the documentation of `JuMP`.
"""
is_solved(opt_model) = termination_status(opt_model) in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]

"""
    optimize_objective(opt_model)::Maybe{Float64}

Shortcut for running JuMP `optimize!` on a model and returning the objective
value, if solved.
"""
function optimize_objective(opt_model)::Maybe{Float64}
    optimize!(opt_model)
    solved_objective_value(opt_model)
end

"""
    get_optmodel_bounds(opt_model)

Returns vectors of the lower and upper bounds of `opt_model` constraints, where
`opt_model` is a JuMP model constructed by e.g.
[`make_optimization_model`](@ref) or [`flux_balance_analysis`](@ref).
"""
get_optmodel_bounds(opt_model) = (
    [-normalized_rhs(lb) for lb in opt_model[:lbs]],
    [normalized_rhs(ub) for ub in opt_model[:ubs]],
)

"""
    set_optmodel_bound!(vidx, opt_model;
        ub::Maybe{Real} = nothing,
        lb::Maybe{Real} = nothing,
    )

Helper function to set the bounds of a variable in the model. Internally calls
`set_normalized_rhs` from JuMP. If the bounds are set to `nothing`, they will
not be changed.
"""
function set_optmodel_bound!(
    vidx,
    opt_model;
    lb::Maybe{Real} = nothing,
    ub::Maybe{Real} = nothing,
)
    isnothing(lb) || set_normalized_rhs(opt_model[:lbs][vidx], -lb)
    isnothing(ub) || set_normalized_rhs(opt_model[:ubs][vidx], ub)
end

"""
    solved_objective_value(opt_model)::Maybe{Float64}

Returns the current objective value of a model, if solved.

# Example
```
solved_objective_value(flux_balance_analysis(model, ...))
```
"""
solved_objective_value(opt_model)::Maybe{Float64} =
    is_solved(opt_model) ? objective_value(opt_model) : nothing

"""
    flux_vector(opt_model)::Maybe{Vector{Float64}}

Returns a vector of fluxes of the model, if solved.

# Example
```
flux_vector(flux_balance_analysis(model, ...))
```
"""
flux_vector(opt_model)::Maybe{Vector{Float64}} =
    is_solved(opt_model) ? value.(opt_model[:x]) : nothing

"""
    flux_dict(model::MetabolicModel, opt_model)::Maybe{Dict{String, Float64}, Nothing}

Returns the fluxes of the model as a reaction-keyed dictionary, if solved.

# Example
```
flux_dict(model, flux_balance_analysis(model, ...))
```
"""
flux_dict(model::MetabolicModel, opt_model)::Maybe{Dict{String,Float64}} =
    is_solved(opt_model) ? Dict(reactions(model) .=> value.(opt_model[:x])) : nothing


text:
using Rauschen
using Test

# Tests will be added here
@testset "Basic Tests" begin
    @testset "Test1" begin
        @testset "Part 1" begin
            @assert true
        end
        @testset "Part 2" begin
            @assert true
        end
    end
end


text:
abstract type AbstractOptic; end

struct OpticRepr <: AbstractOptic
  l::Function
  r::Function
end
OpticRepr(o::OpticRepr) = o
OpticRepr(o::AbstractOptic) =
  OpticRepr(a->left(o, a),
    (m, b′)->right(o, m, b′))

left(o::OpticRepr, a) = o.l(a)
right(o::OpticRepr, m, b′) = o.r(m, b′)

function ⨟(o₁::AbstractOptic, o₂::AbstractOptic)
  OpticRepr(
    function (a)
      x = left(o₁, a)
      y = left(o₂, x[2])
      (x[1], y[1]), y[2]
    end,
    ((m₁, m₂), c′)->right(o₁, m₁, right(o₂, m₂, c′))
  )
end


text:
using Parameters

@component struct Test1
    p::Int
end

Test1() = Test1(0)

@component @with_kw struct Test2
    p::Int = 1
end

@shared_component struct Test3
    p::Int
end
Test3() = Test3(1)

@shared_component @with_kw struct Test4
    p::Int = 1
end

@test Overseer.component_type(Test1) == Component
@test Overseer.component_type(Test2) == Component
@test Overseer.component_type(Test3) == SharedComponent
@test Overseer.component_type(Test4) == SharedComponent

c1 = Overseer.component_type(Test1){Test1}()
c2 = Overseer.component_type(Test2){Test2}()
c3 = Overseer.component_type(Test3){Test3}()
c4 = Overseer.component_type(Test4){Test4}()

entities1 = [Entity(i) for i in 2:2:10]
entities2 = [Entity(i) for i in 10:3:20]
entities3 = [Entity(i) for i in 3:10]
entities4 = [Entity(1)]

for (c, es) in zip((c1, c2, c3, c4), (entities1, entities2, entities3, entities4))
    for e in es
        c[e] = eltype(c)()
    end
end

for (c, es) in zip((c1, c2, c3, c4), (entities1, entities2, entities3, entities4))
    for e in es
        @test in(e, c)
    end
end

t = 0
for e in @entities_in(((c1 && c3) || c4) && !c2)
    global t += e.id
end
@test t == 4+6+8+1

t = 0
for e in @entities_in((c1 || c3) && !c2)
    if e in c1
        global t += e.id
    end
    if e in c3
        global t += c3[e].p
    end
end
@test t == 27

t = 0
for e in @entities_in(c1)
    global t += e.id
end
@test t == sum(2:2:10)

@test pop!(c1, Entity(10)) == Test1()

@test length(c1) == length(entities1) - 1

@test pop!(c2, Entity(10)) == Test2()

@test length(c2) == length(entities2) - 1
@test c1[1] == Test1()
@test c3[1] == Test3()

c2[Entity(13)] = Test2(50)
@test c2[Entity(13)] == Test2(50)

c3[Entity(13)] = Test3(50)
@test c3[Entity(13)] == Test3(50)

pop!(c3, Entity(13))
@test !in(Entity(13), c3)

empty!(c1)
@test isempty(c1)

empty!(c3)
@test isempty(c3)

# swap_ordering
c2[Entity(12)] = Test2()

@test_throws BoundsError swap_order!(c2, Entity(14), Entity(15))
@test_throws BoundsError swap_order!(c2, Entity(13), Entity(14))

orig1 = c2[Entity(12)]
orig2 = c2[Entity(13)]

orig_id1 = c2.indices[12]
orig_id2 = c2.indices[13]

swap_order!(c2, Entity(12), Entity(13))
@test c2[Entity(12)] == orig1
@test c2[Entity(13)] == orig2

@test c2.indices[13] == orig_id1
@test c2.indices[12] == orig_id2

c3[Entity(12)] = Test3()
c3[Entity(13)] = Test3(50)

@test_throws BoundsError swap_order!(c3, Entity(14), Entity(15))
@test_throws BoundsError swap_order!(c3, Entity(13), Entity(14))

orig1 = c3[Entity(12)]
orig2 = c3[Entity(13)]

orig_id1 = c3.indices[12]
orig_id2 = c3.indices[13]

swap_order!(c3, Entity(12), Entity(13))
@test c3[Entity(12)] == orig1
@test c3[Entity(13)] == orig2

@test c3.indices[13] == orig_id1
@test c3.indices[12] == orig_id2

# Issue 4: collect() and iterator length
@testset "collect" begin
    e1 = Entity(1)
    e2 = Entity(2)
    e3 = Entity(3)

    comp1 = Overseer.component_type(Test1){Test1}()
    comp2 = Overseer.component_type(Test2){Test2}()
    comp1[e1] = Test1(1)
    comp1[e2] = Test1(1); comp2[e2] = Test2(1)
    comp2[e3] = Test2(1)

    iter = @entities_in(comp1 && comp2)
    es = collect(iter)
    @test es == [e2]
    @test eltype(es) == Entity
end