text:
"""
DR_LagrangeDual

Lagrangian dual method for dual decomposition. This `mutable struct` constains:
- `block_model::BlockModel` object
- `var_to_index` mapping coupling variable to the index wrt the master problem
- `masiter::Int` sets the maximum number of iterations
- `tol::Float64` sets the relative tolerance for termination
- `tree` keeps Tree information
"""

mutable struct DR_LagrangeDual <: DD.AbstractLagrangeDual
    block_model::DR_BlockModel
    var_to_index::Dict{Tuple{Int,Any},Int} # maps coupling variable to the index wrt the master problem
    heuristics::Vector{Type}
    subsolve_time::Vector{Dict{Int,Float64}}
    subcomm_time::Vector{Float64}
    subobj_value::Vector{Float64}
    master_time::Vector{Float64}

    dh::Union{Nothing, DD.DataHelper}

    tree::Union{Nothing,DD.Tree{DR_TreeNode}}
    subtrees::Union{Nothing,Dict{Int,DD.SubTree}}
    P_model::Union{Nothing, JuMP.Model}

    function DR_LagrangeDual()
        LD = new()
        LD.block_model = DR_BlockModel()
        LD.var_to_index = Dict()
        LD.heuristics = []
        LD.subsolve_time = []
        LD.subcomm_time = []
        LD.subobj_value = []
        LD.master_time = []

        LD.dh = nothing

        LD.tree = nothing
        LD.subtrees = nothing
        LD.P_model = nothing

        return LD
    end
end

function DR_LagrangeDual(tree::DD.Tree{DR_TreeNode})
    LD = DR_LagrangeDual()
    DD.add_tree!(LD, tree)
    return LD
end

function DD.get_solution!(LD::DR_LagrangeDual, method::BM.AbstractMethod)
    LD.block_model.dual_solution = copy(BM.get_solution(method))
    bundle = BM.get_model(method)
    model = BM.get_model(bundle)
    P = model[:P]
    for id in axes(P)[1]
        LD.block_model.P_solution[id] = JuMP.value(P[id])
    end
end

function DD.add_constraints!(LD::DR_LagrangeDual, method::DD.BundleMaster)
    node_to_couple = sort_couple_by_label(LD.tree, LD.block_model.variables_by_couple)

    model = BM.get_jump_model(method.inner)
    @variable(model, P[2:length(LD.tree.nodes)] >= 0)

    for (id, node) in LD.tree.nodes
        add_non_anticipativity!(LD, model, node, node_to_couple[id])
        add_ambiguity!(LD.tree, model, node, node.set)
    end
    for (id, node) in LD.tree.nodes
        add_ambiguity_link!(model, node, node.set)
    end

    #JuMP.print(model)
end

function sort_couple_by_label(tree::DD.Tree, variables_by_couple::Dict{Any,Vector{DD.CouplingVariableKey}})::Dict{Int,Vector{Any}}
    node_to_couple = Dict{Int, Vector{Any}}()
    for (id, nodes) in tree.nodes
        node_to_couple[id] = Vector{Any}()
    end
    for (couple_id, keys) in variables_by_couple
        loc = findfirst(x -> x=='_', couple_id)
        node_id = parse(Int, couple_id[2:loc-1])
        push!(node_to_couple[node_id], couple_id)
    end
    return node_to_couple
end

function add_non_anticipativity!(LD::DR_LagrangeDual, m::JuMP.Model, node::DR_TreeNode, couple_ids::Vector{Any})
    λ = m[:x]
    P = m[:P]

    if DD.check_root(node)
        for couple_id in couple_ids
            vars = LD.block_model.variables_by_couple[couple_id]
            @constraint(m, sum(λ[DD.index_of_λ(LD, v)] for v in vars) == get_cost(node, couple_id), base_name = couple_id)
        end
    else
        for couple_id in couple_ids
            vars = LD.block_model.variables_by_couple[couple_id]
            @constraint(m, sum(λ[DD.index_of_λ(LD, v)] for v in vars) == get_cost(node, couple_id) * P[DD.get_id(node)], base_name = couple_id)
        end
    end
end

function add_ambiguity!(tree::DD.Tree{DR_TreeNode}, m::JuMP.Model, node::DR_TreeNode, set::AbstractAmbiguitySet) end
function add_ambiguity!(tree::DD.Tree{DR_TreeNode}, m::JuMP.Model, node::DR_TreeNode, set::Nothing) end

function add_ambiguity!(tree::DD.Tree{DR_TreeNode}, m::JuMP.Model, node::DR_TreeNode, set::WassersteinSet)
    P = m[:P]
    if DD.check_root(node)
        @variable(m, w[id=DD.get_children(node),s=1:set.N] >= 0, base_name = "n1_w")
        @constraint(m, sum( sum( w[id, s] * set.norm_func(tree.nodes[id].ξ, set.samples[s].ξ) for s in 1:set.N) for id in DD.get_children(node)) <= set.ϵ)
        @constraint(m, [s=1:set.N], sum( w[id, s] for id in DD.get_children(node)) == set.samples[s].p )
        @constraint(m, sum( P[child] for child in DD.get_children(node)) == 1)
        JuMP.unregister(m, :w)
    elseif !DD.check_leaf(node)
        node_id = DD.get_id(node)
        @variable(m, w[id=DD.get_children(node),s=1:set.N] >= 0, base_name = "n$(node_id)_w")
        @constraint(m, sum( sum( w[id, s] * set.norm_func(tree.nodes[id].ξ, set.samples[s].ξ) for s in 1:set.N) for id in DD.get_children(node)) <= set.ϵ * P[node_id])
        @constraint(m, [s=1:set.N], sum( w[id, s] for id in DD.get_children(node)) == set.samples[s].p * P[node_id])
        @constraint(m, sum( P[child] for child in DD.get_children(node)) == P[node_id])
        JuMP.unregister(m, :w)
    end
end

function add_ambiguity_link!(m::JuMP.Model, node::DR_TreeNode, set::Nothing) end
function add_ambiguity_link!(m::JuMP.Model, node::DR_TreeNode, set::WassersteinSet)
    if !DD.check_leaf(node)
        id = DD.get_id(node)
        P = m[:P]
        w = Dict{Tuple{Int,Int},JuMP.VariableRef}( (child, s) => JuMP.variable_by_name(m, "n$(id)_w[$(child),$(s)]") for child=DD.get_children(node), s in 1:set.N)
        @constraint(m, [child=DD.get_children(node)], sum( w[child, s] for s in 1:set.N) == P[child])
    end
end

function initialize_bundle(tree::DD.Tree{DR_TreeNode}, LD::DR_LagrangeDual, Optimizer=nothing)::Array{Float64,1}
    n = DD.parallel.sum(DD.num_coupling_variables(LD.block_model))
    bundle_init = Array{Float64,1}(undef, n)
    variable_keys = [v.key for v in LD.block_model.coupling_variables]
    all_variable_keys = DD.parallel.allcollect(variable_keys)
    if DD.parallel.is_root()
        P = get_feasible_P(tree, LD, Optimizer)
        #println(P)
        for key in all_variable_keys
            i = LD.var_to_index[(key.block_id,key.coupling_id)]
            N = length(LD.block_model.variables_by_couple[key.coupling_id])
            couple_id = key.coupling_id

            loc = findfirst(x -> x=='_', couple_id)
            node_id = parse(Int, couple_id[2:loc-1])

            if DD.check_root(tree.nodes[node_id])
                bundle_init[i] =  get_cost(tree.nodes[node_id], couple_id) / N
            else
                bundle_init[i] =  get_cost(tree.nodes[node_id], couple_id) / N * P[node_id] #P[corresponding leaf node] ?
            end
            """
            for (id, node) in tree.nodes
                if DD.check_leaf(node)
                    bundle_init[i] =  get_cost(tree.nodes[node_id], couple_id)* P[DD.get_id(node)]
                end
            end
            """
        end
        DD.parallel.bcast(bundle_init)
    else
        bundle_init = DD.parallel.bcast(nothing)
    end
    return bundle_init
end

function get_feasible_P(tree::DD.Tree{DR_TreeNode}, LD::DR_LagrangeDual, Optimizer)::Dict{Int,Float64}
    if isnothing(Optimizer)
        model = JuMP.Model(GLPK.Optimizer)
    else
        model = JuMP.Model(Optimizer)
    end

    @variable(model, P[2:length(tree.nodes)] >= 0)
    for (id, node) in tree.nodes
        add_ambiguity!(tree, model, node, node.set)
    end
    for (id, node) in tree.nodes
        add_ambiguity_link!(model, node, node.set)
    end

    @objective(model, Min, 0)

    LD.P_model = model
    JuMP.optimize!(model)
    #JuMP.print(model)

    Pref = Dict()
    P = model[:P]

    for id in 2:length(tree.nodes)
        Pref[id] = JuMP.value(P[id])
    end
    return Pref
end


text:
"""
Check the correctness of the filtering implementations vs. basic readable math code
"""

using Test
using ProbNumDiffEq
using LinearAlgebra

@testset "PREDICT" begin
    # Setup
    d = 5
    m = rand(d)
    L_p = Matrix(LowerTriangular(rand(d, d)))
    P = L_p * L_p'

    A = rand(d, d)
    L_Q = Matrix(LowerTriangular(rand(d, d)))
    Q = L_Q * L_Q'

    # PREDICT
    m_p = A * m
    P_p = A * P * A' + Q

    x_curr = Gaussian(m, P)
    x_out = copy(x_curr)

    @testset "predict" begin
        x_out = ProbNumDiffEq.predict(x_curr, A, Q)
        @test m_p == x_out.μ
        @test P_p == x_out.Σ
    end

    @testset "predict! with SRMatrix" begin
        x_curr = Gaussian(m, SRMatrix(L_p))
        x_out = copy(x_curr)
        Q_SR = SRMatrix(L_Q)
        cache = SRMatrix(zeros(d, 2d))
        ProbNumDiffEq.predict!(x_out, x_curr, A, Q_SR, cache)
        @test m_p == x_out.μ
        @test P_p ≈ x_out.Σ
    end
end

@testset "UPDATE" begin
    # Setup
    d = 5
    m_p = rand(d)
    L_P_p = Matrix(LowerTriangular(rand(d, d)))
    P_p = L_P_p * L_P_p'

    # Measure
    o = 3
    H = rand(o, d)
    L_R = 0.0 * rand(o, o)
    R = L_R * L_R'

    z_data = zeros(o)
    z = H * m_p
    S = Symmetric(H * P_p * H' + R)

    # UPDATE
    S_inv = inv(S)
    K = P_p * H' * S_inv
    m = m_p + K * (z_data .- z)
    P = P_p - K * S * K'

    x_pred = Gaussian(m_p, P_p)
    x_out = copy(x_pred)
    measurement = Gaussian(z, S)

    @testset "update" begin
        x_out = ProbNumDiffEq.update(x_pred, measurement, H)
        @test m ≈ x_out.μ
        @test P ≈ x_out.Σ
    end

    @testset "update!" begin
        K_cache = copy(K)
        M_cache = zeros(d, d)
        m_tmp = copy(measurement)
        ProbNumDiffEq.update!(x_out, x_pred, measurement, H, K_cache, M_cache, m_tmp)
        @test m ≈ x_out.μ
        @test P ≈ x_out.Σ
    end
end

@testset "SMOOTH" begin
    # Setup
    d = 5
    m, m_s = rand(d), rand(d)
    L_P, L_P_s = Matrix(LowerTriangular(rand(d, d))), Matrix(LowerTriangular(rand(d, d)))
    P, P_s = L_P * L_P', L_P_s * L_P_s'

    A = rand(d, d)
    L_Q = Matrix(LowerTriangular(rand(d, d)))
    Q = L_Q * L_Q'
    Q_SR = SRMatrix(L_Q)

    # PREDICT first
    m_p = A * m
    P_p = A * P * A' + Q

    # SMOOTH
    G = P * A' * inv(P_p)
    m_smoothed = m + G * (m_s - m_p)
    P_smoothed = P + G * (P_s - P_p) * G'

    x_curr = Gaussian(m, SRMatrix(L_P))
    x_smoothed = Gaussian(m_s, SRMatrix(L_P_s))

    @testset "smooth" begin
        x_out, _ = ProbNumDiffEq.smooth(x_curr, x_smoothed, A, Q)
        @test m_smoothed ≈ x_out.μ
        @test P_smoothed ≈ x_out.Σ
    end
    @testset "smooth with SRMatrix" begin
        x_curr = Gaussian(m, SRMatrix(L_P))
        x_smoothed = Gaussian(m_s, SRMatrix(L_P_s))
        x_out, _ = ProbNumDiffEq.smooth(x_curr, x_smoothed, A, Q_SR)
        @test m_smoothed ≈ x_out.μ
        @test P_smoothed ≈ x_out.Σ
    end
end

ext: jl


text:
using Pkg
packages = ["ArgParse", "MriResearchTools"]
Pkg.add(packages)
using ArgParse

max_stars_repo_name: aswinnarayanan/neurocontainers

max_issues_repo_name: aswinnarayanan/neurocontainers

max_forks_repo_name: aswinnarayanan/neurocontainers


text:
using Test
using SimpleTools

sq = Dict{Int,Int}()
for k=1:10
    sq[k] = k*k
end

rt = Dict{Int,Int}()
for k=1:10
    rt[k*k] = k
end

f = rt*sq
for k=1:10
    @test f[k] == k
end

y1 = [x^2 for x=0:.01:3]
y2 = [mod(y,1) for y in y1 ]
make_continuous!(y2,1)
@test y1==y2

@test flush_print(23,5) == "   23"
@test flush_print(23,5,false) == "23   "

A = dcat(1,1,1)
M = [1 0 0 ; 0 1 0 ; 0 0 1]
@test A==M

@test sage(A)==nothing

@test sum(eye(5)) == 5
@test cofactor_det(2eye(Int,3))==8

using Polynomials
x = Poly([0,1])
A = [1 1; 3 2]
@test char_poly(A) == (x-1)*(x-2)-3


text:
#TODO sort treatmentplan by default
function generateTreatmentplan!(visits::IndexedTable,row,bestord,columns::Dict{Symbol,String},patientID)
        for col in columns
            if !ismissing(row[col[1]])
                if row[col[1]] == 1 || rand() < row[col[1]]
                    push!(rows(visits),(intID=length(visits)+1,patientID=patientID,bestord=bestord[1],bestord_date=bestord[2],req_type=col[2]))
                end
            end
        end

end

"""
    readPatientTable(path::String,sheet::String,columns::Dict{Symbol,String},mastercalendar::Dict{Int64,Date})

Read aggregated table data of the format found in 'test/Sample data/PatientOverview_test.xlsx'.
Produces a table of patients and a table of visists with the patients best.ord./deadline spread uniformly across everyday of the master calendar.
columns is a Dict{Symbol,String} making it possible to translate the column headers into types of visits

# Examples
```julia-repl
julia> columns = Dict(:Consultation => "Consultation" , :Telefon =>"Telephone" , :TTE => "TTE", :AEKG => "AEKG",  :MR => "MR", :Holter=>"Holter")
Dict{Symbol,String} with 6 entries:
  :MR           => "MR"
  :Telefon      => "Telephone"
  :AEKG         => "AEKG"
  :Holter       => "Holter"
  :Consultation => "Consultation"
  :TTE          => "TTE"

julia> patients, visits = readPatientTable(path,sheet,columns,mastercalendar)
(Table with 40 rows, 4 columns:
intID  diagnosis  bestord  bestord_date
───────────────────────────────────────
1      "test"     31       2019-01-31
2      "test"     22       2019-01-22
3      "test"     16       2019-01-16
4      "test"     28       2019-01-28
5      "test"     8        2019-01-08
6      "test"     2        2019-01-02
⋮
37     "test"     8        2019-01-08
38     "test"     2        2019-01-02
39     "test"     25       2019-01-25
40     "test"     17       2019-01-17, Table with 130 rows, 5 columns:
intID  patientID  bestord  bestord_date  req_type
───────────────────────────────────────────────────────
1      1          31       2019-01-31    "AEKG"
2      1          31       2019-01-31    "Consultation"
3      1          31       2019-01-31    "TTE"
4      2          22       2019-01-22    "AEKG"
5      2          22       2019-01-22    "Consultation"
6      2          22       2019-01-22    "TTE"
⋮
127    39         25       2019-01-25    "TTE"
128    40         17       2019-01-17    "MR"
129    40         17       2019-01-17    "Consultation"
130    40         17       2019-01-17    "TTE")
```
"""
function readPatientTable(path,sheet,columns,mastercalendar)
    patientOverview = DataFrame(XLSX.readtable(path,sheet)...)
    patients = JuliaDB.table((intID=Int64[],diagnosis=String[],bestord=Int64[],bestord_date=Date[]);pkey=[:intID])
    visits =  JuliaDB.table((intID=Int64[],patientID=Int64[],bestord=Int64[],bestord_date=Date[],req_type=String[]);pkey=[:intID])
    for row in eachrow(patientOverview)
        if row.Kategori === missing
            continue
        end
        for i = 1:Int(row.Antal)
            bestord = rand(mastercalendar)
            intID=length(patients)+1
            generateTreatmentplan!(visits,row,bestord,columns,intID)
            push!(rows(patients),(intID=intID,diagnosis="test",bestord=bestord[1],bestord_date=bestord[2]))
        end
    end
    patients, visits
end
"""
    readWorkPattern(path_resourceOverview,sheet_amb)

Creates a list of resources, a workpattern which is a list of timeslots in even and odd weeks, and lastly an empty list of timeslots. The data is take from sheet 'sheet_amb' in the excel file given in path_resourceOverview. The format for the files can be seen in 'test/Sample data/GUCHamb_Timeslot_test.xlsx'

# Examples
```julia-repl
julia> readWorkPattern(path_resourceOverview,sheet_amb)
(Table with 3 rows, 4 columns:
intID  id                     type            name
──────────────────────────────────────────────────────
1      "Consultation_GUCH L"  "Consultation"  "GUCH L"
2      "Consultation_MRS"     "Consultation"  "MRS"
3      "Consultation_NV"      "Consultation"  "NV",
Table with 0 rows, 7 columns:
resourceID  type  dayID  timeslotID  startTime  endTime  booked
───────────────────────────────────────────────────────────────,
Table with 82 rows, 7 columns:
resourceID  type            weekdayID  oddWeek  timeslotID  startTime  endTime
───────────────────────────────────────────────────────────────────────────────
1           "Consultation"  5          false    1           09:00:00   09:30:00
1           "Consultation"  5          true     2           09:00:00   09:30:00
⋮
3           "Consultation"  1          false    81          11:15:00   12:00:00
3           "Consultation"  1          false    82          11:15:00   12:00:00)
```
"""
function readWorkPattern(path::String,sheet::String)
    resources = JuliaDB.table((intID=Int64[],id=String[],type=String[],name=String[],);pkey=[:intID])
    #resources = JuliaDB.table((intID=Int64[],id=String[],type=String[],name=String[],qualifications=Dict{String,String}[],);pkey=[:intID])
    workpattern = JuliaDB.table((resourceID=Int64[],type = String[],weekdayID=Int64[],oddWeek=Bool[],timeslotID=Int64[],startTime=Time[],endTime=Time[]);pkey=[:resourceID,:weekdayID,:timeslotID])
    timeslots = JuliaDB.table((resourceID=Int64[],type = String[], dayID=Int64[],timeslotID=Int64[],startTime=Time[],endTime=Time[],booked=Bool[]);pkey=[:resourceID,:dayID,:timeslotID])
    readWorkPattern!(resources,timeslots,workpattern,path::String,sheet::String,)
end

"""
    readWorkPattern!(resources,timeslots,workpattern,path::String,sheet::String,)

Adds additional resources to the list of resources and creates workpattern for these added to workpattern

"""
function readWorkPattern!(resources,timeslots,workpattern,path::String,sheet::String,)
        columns = Dict{String,Int}("Monday"=>1,"Tuesday"=> 2, "Wednesday" => 3, "Thursday" => 4, "Friday"=> 5)
        wp_df = DataFrame(XLSX.readtable(path,sheet)...)
        for resource_df in DataFrames.groupby(wp_df,[:Type,:Resource],skipmissing = true)
            cur_resourceid = string(resource_df[1,:Type],"_",resource_df[1,:Resource])
            tempID = filter(i -> i.id == cur_resourceid,resources)
            if length(tempID) == 0
                type = string(resource_df[1,:Type])
                name = string(resource_df[1,:Resource])
                resourceid = string(resource_df[1,:Type],"_",resource_df[1,:Resource])
                intID = length(resources)+1

                push!(rows(resources),(intID=intID,id=resourceid,type=type,name=name,qualifications=Dict("name" => name,"type"=>type)))

            elseif length(tempID) == 1
                intID = first(tempID)
            else
                throw("Error: Multiple resources with same ID")
            end
            for day in columns
                if ismissing(resource_df[1,Symbol(string(day[1],"_start"))]) continue end;
                for row in eachrow(resource_df)
                    if ismissing(row[Symbol(string(day[1],"_start"))]) ||row[Symbol(string(day[1],"_start"))]=="PAUSE"
                         continue
                    end
                    push!(rows(workpattern),(resourceID = intID,type = string(resource_df[1,:Type]),weekdayID = day[2],oddWeek = !(row.Weeks in ("All","Even")),timeslotID = length(workpattern)+1, startTime = row[Symbol(string(day[1],"_start"))],endTime = row[Symbol(string(day[1],"_end"))]))
                    push!(rows(workpattern),(resourceID = intID,type = string(resource_df[1,:Type]),weekdayID = day[2],oddWeek = (row.Weeks in ("All","Odd")),timeslotID = length(workpattern)+1, startTime = row[Symbol(string(day[1],"_start"))],endTime = row[Symbol(string(day[1],"_end"))]))
                end
            end
        end
        resources, timeslots, workpattern
end

"""
    generateCalendarFromPattern!(timeslots,resources::IndexedTable,workpattern,masterCalendar::Dict{Int64,Date},occupancyrate = 0)

Populates the variable timeslots with appropriate timeslots from workpattern list for everyday in the mastercalendar. occupancyrate is the part of the timeslots that should be randomly booked (meaning that they are completely taken out of the planning problem). 0 means no booked, 0.5 means half of all slots are booked and 1 means all slots are booked.

"""
function generateCalendarFromPattern!(timeslots,resources::IndexedTable,workpattern,masterCalendar::Dict{Int64,Date},occupancyrate = 0)#
    for cur_resource in rows(resources)
        for _day in sort(collect(masterCalendar))
            xtemp  =  @from i in workpattern begin
                    @where i.resourceID == cur_resource.intID && i.weekdayID == dayofweek(_day[2]) && week(_day[2])%2 == i.oddWeek
                    @select {resourceID=i.resourceID,type= i.type,startTime= i.startTime,endTime = i.endTime}
                    @collect
                end
            if length(xtemp) > 0
                x = table(xtemp)
                y = (:dayID => fill(_day[1],length(x)),:timeslotID => [i for i in (length(rows(timeslots))+1):(length(rows(timeslots))+length(x))],:booked => [y <= (1-occupancyrate)  for y in rand(length(x))] )
                x = transform(x, y)
                append!(rows(timeslots),rows(x))
            end
        end
    end
end


text:
"This script implement a threads version for the Laguerre-Kravchuk function, so we can distribute the calculation over a multithreaded processor"

using PixelsGT # Call the parent module

# Splits the calculations on lower [0, 2j] and upper [2j, 4j] Laguerre-Kravchuk modes

"`rd(Q₁, Q₂, q₁, q₂, n, m, θ, j)` gives a tetradimensional array of lower modes"
rd(Q₁, Q₂, q₁, q₂, n, m, θ, j) = real(conj(LKd(Q₁, Q₂, n, m, j)) * exp(-1im * θ * m) * LKd(q₁, q₂, n, m, j))

"`ru(Q₁, Q₂, q₁, q₂, n, m, θ, j)` gives a tetradimensional array of upper modes"
ru(Q₁, Q₂, q₁, q₂, n, m, θ, j) = real(conj(LKu(Q₁, Q₂, n, m, j)) * exp(-1im * θ * m) * LKu(q₁, q₂, n, m, j))

export RdM, RuM

"`RdM(θ, j)` gives the twodimensional array of lower modes at angle θ for a j-dimensional representation"
function RdM(θ, j)
    m = zeros((Integer(2 * j + 1), Integer(2 * j + 1), Integer(2 * j + 1), Integer(2 * j + 1)))
    Threads.@threads  for Q₁ in 1:Integer(2*j+1)
        for Q₂ in 1:Integer(2*j+1)
            for q₁ in 1:Integer(2*j+1)
                for q₂ in 1:Integer(2*j+1)
                    m[Q₁, Q₂, q₁, q₂] = sum([rd(Q₁-j-1, Q₂-j-1, q₁-j-1, q₂-j-1, n, m, θ, j) for n in 0:Integer(2*j) for m in -Integer(n):2:Integer(n)])
                end
            end
        end
    end
    return m
end

"`RuM(θ, j)` gives the twodimensional array of lower modes at angle θ for a j-dimensional representation"
function RuM(θ, j)
    m = zeros((Integer(2 * j + 1), Integer(2 * j + 1), Integer(2 * j + 1), Integer(2 * j + 1)))
    Threads.@threads  for Q₁ in 1:Integer(2*j+1)
        for Q₂ in 1:Integer(2*j+1)
            for q₁ in 1:Integer(2*j+1)
                for q₂ in 1:Integer(2*j+1)
                    m[Q₁, Q₂, q₁, q₂] = sum([ru(Q₁-j-1, Q₂-j-1, q₁-j-1, q₂-j-1, n, m, θ, j) for n in Integer(2*j + 1):Integer(4*j) for m in -Integer(4 * j - n):2:Integer(4 * j - n)])
                end
            end
        end
    end
    return m
end


text:
# map-reduce

# GPUArrays' mapreduce methods build on `Base.mapreducedim!`, but with an additional
# argument `init` value to avoid eager initialization of `R` (if set to something).
mapreducedim!(f, op, R::AbstractGPUArray, A::AbstractArray, init=nothing) = error("Not implemented") # COV_EXCL_LINE
Base.mapreducedim!(f, op, R::AbstractGPUArray, A::AbstractArray) = mapreducedim!(f, op, R, A)

neutral_element(op, T) =
    error("""GPUArrays.jl needs to know the neutral element for your operator `$op`.
             Please pass it as an explicit argument to (if possible), or register it
             globally your operator by defining `GPUArrays.neutral_element(::typeof($op), T)`.""")
neutral_element(::typeof(Base.:(|)), T) = zero(T)
neutral_element(::typeof(Base.:(+)), T) = zero(T)
neutral_element(::typeof(Base.add_sum), T) = zero(T)
neutral_element(::typeof(Base.:(&)), T) = one(T)
neutral_element(::typeof(Base.:(*)), T) = one(T)
neutral_element(::typeof(Base.mul_prod), T) = one(T)
neutral_element(::typeof(Base.min), T) = typemax(T)
neutral_element(::typeof(Base.max), T) = typemin(T)

function Base.mapreduce(f, op, A::AbstractGPUArray; dims=:, init=nothing)
    # figure out the destination container type by looking at the initializer element,
    # or by relying on inference to reason through the map and reduce functions.
    if init === nothing
        ET = Base.promote_op(f, eltype(A))
        ET = Base.promote_op(op, ET, ET)
        (ET === Union{} || ET === Any) &&
            error("mapreduce cannot figure the output element type, please pass an explicit init value")

        init = neutral_element(op, ET)
    else
        ET = typeof(init)
    end

    sz = size(A)
    red = ntuple(i->(dims==Colon() || i in dims) ? 1 : sz[i], ndims(A))
    R = similar(A, ET, red)
    mapreducedim!(f, op, R, A, init)

    if dims==Colon()
        @allowscalar R[]
    else
        R
    end
end

Base.any(A::AbstractGPUArray{Bool}) = mapreduce(identity, |, A)
Base.all(A::AbstractGPUArray{Bool}) = mapreduce(identity, &, A)

Base.any(f::Function, A::AbstractGPUArray) = mapreduce(f, |, A)
Base.all(f::Function, A::AbstractGPUArray) = mapreduce(f, &, A)
Base.count(pred::Function, A::AbstractGPUArray) = mapreduce(pred, +, A; init = 0)

Base.:(==)(A::AbstractGPUArray, B::AbstractGPUArray) = Bool(mapreduce(==, &, A, B))

# avoid calling into `initarray!``
Base.sum!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(Base.add_sum, R, A)
Base.prod!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(Base.mul_prod, R, A)
Base.maximum!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(max, R, A)
Base.minimum!(R::AbstractGPUArray, A::AbstractGPUArray) = Base.reducedim!(min, R, A)

LinearAlgebra.ishermitian(A::AbstractGPUMatrix) = mapreduce(==, &, A, adjoint(A))


text:
using KdotP, Crystalline, LinearAlgebra, GellMannMatrices

timereversal=true
@profview for sgnum in 1:230
    lgirsd = lgirreps(sgnum, Val(3))
    printed_sg = false
    for (klab, lgirs) in lgirsd
        timereversal && (lgirs = realify(lgirs))
        for lgir in lgirs
            try
                kdotp(lgir; timereversal)
            catch err
                printed_sg || (print("SG ", sgnum, ":\n"); printed_sg=true)
                print("   ", Crystalline.formatirreplabel(label(lgir)), ": ")
                println(err)
            end
        end
    end
    printed_sg && println()
end


text:
##
## rand for Flint-Finite fields
##
## fq_nmod has no base ring
function rand(R::FqNmodFiniteField)
  #gen is not a generator for the group!
  r = R(0)
  for i=0:degree(R)
    r *= gen(R)
    r += rand(1:characteristic(R))
  end

  return r
end

function rand(R::FqFiniteField)
  #gen is not a generator for the group!
  r = R(0)
  for i=0:degree(R)
    r *= gen(R)
    r += rand(1:characteristic(R))
  end

  return r
end


text:
using JacobiDavidson
using Documenter

makedocs(;
    modules=[JacobiDavidson],
    authors="Harmen Stoppels",
    repo="https://github.com/haampie/JacobiDavidson.jl/blob/{commit}{path}#{line}",
    sitename="JacobiDavidson.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://haampie.github.io/JacobiDavidson.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Correction equation" => "solvers.md"
    ],
    doctest = false
)

deploydocs(;
    repo="github.com/haampie/JacobiDavidson.jl",
)

max_stars_repo_name: haampie/jacobi-davidson-julia

max_forks_repo_name: haampie/jacobi-davidson-julia


text:
using ShuHua
using ShuHua.Parse
using ShuHua.Write

pr = parse_slide(ShuHua.project("example", "simple.md"))
pr.meta
html(pr)
write("index.html", pr)
pr = parse_slide(ShuHua.project("example", "huan"))
html(pr)
using ShuHua.BuildTools

pr.meta
build(pr)


text:
using Animated_CMA_ES, BenchmarkTools, Profile

N=50
function h(x)
    h=0.0
    for i in 1:N
        h=h+(x[i]+Float64(i))^2
    end
    return h
end

cma_es(h,zeros(N),tol=1e-5)
new_cma_es(h,zeros(N),tol=1e-5)

@btime cma_es(h,zeros(N),tol=1e-5)
@btime new_cma_es(h,zeros(N),tol=1e-5)

function prof()
    cma_es(h,zeros(100),tol=1e-7)
end

function new_prof()
    new_cma_es(h,zeros(100),tol=1e-7)
end

Profile.clear()
Profile.init(n=10^8,delay=0.000000001)
@profile prof()
Profile.print(noisefloor=2.0)

Profile.clear()
Profile.init(n=10^8,delay=0.000000001)
@profile new_prof()
Profile.print(noisefloor=2.0)


text:
#LBER_F( void )
#ber_free LDAP_P((
#    BerElement *ber,
#    int freebuf ));
function ber_free(ber_element_handle::Ptr{Cvoid}, freebuf::Integer)
    ccall((:ber_free, liblber), Cvoid, (Ptr{Cvoid}, Cint), ber_element_handle, freebuf)
end


text:
"""
"""
abstract type AssemblyStrategy end

"""
"""
function row_map(a::AssemblyStrategy,row)
  @abstractmethod
end

"""
"""
function col_map(a::AssemblyStrategy,col)
  @abstractmethod
end

"""
"""
function row_mask(a::AssemblyStrategy,row)
  @abstractmethod
end

"""
"""
function col_mask(a::AssemblyStrategy,col)
  @abstractmethod
end

@inline function map_rows!(gids,a::AssemblyStrategy,rows)
  u = -one(eltype(gids))
  for i in eachindex(rows)
    ri = rows[i]
    if ri>0 && row_mask(a,ri)
      gids[i] = row_map(a,ri)
    else
      gids[i] = u
    end
  end
  nothing
end

@inline function map_cols!(gids,a::AssemblyStrategy,cols)
  u = -one(eltype(gids))
  for i in eachindex(cols)
    ri = cols[i]
    if ri>0 && col_mask(a,ri)
      gids[i] = col_map(a,ri)
    else
      gids[i] = u
    end
  end
  nothing
end

function map_cell_rows(strategy::AssemblyStrategy,cell_ids)
  k = AssemblyStrategyMap{:rows}(strategy)
  lazy_map(k,cell_ids)
end

function map_cell_cols(strategy::AssemblyStrategy,cell_ids)
  k = AssemblyStrategyMap{:cols}(strategy)
  lazy_map(k,cell_ids)
end

struct AssemblyStrategyMap{S,T} <: Map
  strategy::T
  function AssemblyStrategyMap{S}(strategy::T) where {S,T}
    new{S,T}(strategy)
  end
end

function Arrays.return_cache(a::AssemblyStrategyMap,ids::AbstractArray)
  gids = similar(ids,Int,axes(ids))
  CachedArray(gids)
end

function Arrays.evaluate!(cache,a::AssemblyStrategyMap{:cols},ids::AbstractArray)
  setsize!(cache,size(ids))
  gids = cache.array
  map_cols!(gids,a.strategy,ids)
  gids
end

function Arrays.evaluate!(cache,a::AssemblyStrategyMap{:rows},ids::AbstractArray)
  setsize!(cache,size(ids))
  gids = cache.array
  map_rows!(gids,a.strategy,ids)
  gids
end

function Arrays.return_cache(k::AssemblyStrategyMap,ids::ArrayBlock)
  fi = testitem(ids)
  ci = return_cache(k,fi)
  gi = evaluate!(ci,k,fi)
  b = Array{typeof(ci),ndims(ids)}(undef,size(ids))
  for i in eachindex(ids.array)
    if ids.touched[i]
      ki = return_cache(k,ids.array[i])
      b[i] = return_cache(k,ids.array[i])
    end
  end
  array = Array{typeof(gi),ndims(ids)}(undef,size(ids))
  ArrayBlock(array,ids.touched), b
end

function Arrays.evaluate!(cache,k::AssemblyStrategyMap,ids::ArrayBlock)
  a,b = cache
  for i in eachindex(ids.array)
    if ids.touched[i]
      a.array[i] = evaluate!(b[i],k,ids.array[i])
    end
  end
  a
end

struct DefaultAssemblyStrategy <: AssemblyStrategy end

row_map(a::DefaultAssemblyStrategy,row) = row

col_map(a::DefaultAssemblyStrategy,col) = col

row_mask(a::DefaultAssemblyStrategy,row) = true

col_mask(a::DefaultAssemblyStrategy,col) = true

map_cell_rows(strategy::DefaultAssemblyStrategy,cell_ids) = cell_ids

map_cell_cols(strategy::DefaultAssemblyStrategy,cell_ids) = cell_ids

struct GenericAssemblyStrategy{A,B,C,D} <: AssemblyStrategy
  row_map::A
  col_map::B
  row_mask::C
  col_mask::D
end

row_map(a::GenericAssemblyStrategy,row) = a.row_map(row)

col_map(a::GenericAssemblyStrategy,col) = a.col_map(col)

row_mask(a::GenericAssemblyStrategy,row) = a.row_mask(row)

col_mask(a::GenericAssemblyStrategy,col) = a.col_mask(col)

"""
"""
abstract type Assembler <: GridapType end

"""
"""
function get_rows(a::Assembler)
  @abstractmethod
end

"""
"""
function get_cols(a::Assembler)
  @abstractmethod
end

num_rows(a::Assembler) = length(get_rows(a))
num_cols(a::Assembler) = length(get_cols(a))

Base.axes(a::Assembler) = (get_rows(a),get_cols(a))
Base.size(a::Assembler) = (num_rows(a),num_cols(a))

"""
"""
function get_assembly_strategy(a::Assembler)
  @abstractmethod
end

"""
"""
function allocate_matrix(a::Assembler,matdata)
  @abstractmethod
end

"""
"""
function allocate_vector(a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function allocate_matrix_and_vector(a::Assembler,data)
  @abstractmethod
end

"""
"""
function assemble_matrix!(A,a::Assembler,matdata)
  @abstractmethod
end

"""
"""
function assemble_matrix_add!(mat,a::Assembler, matdata)
  @abstractmethod
end

"""
"""
function assemble_vector!(b,a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function assemble_vector_add!(b,a::Assembler,vecdata)
  @abstractmethod
end

"""
"""
function assemble_matrix_and_vector!(A,b,a::Assembler, data)
  @abstractmethod
end

function assemble_matrix_and_vector_add!(A,b,a::Assembler, data)
  @abstractmethod
end

"""
"""
function assemble_matrix(a::Assembler,matdata)
  A = allocate_matrix(a,matdata)
  assemble_matrix!(A,a,matdata)
  A
end

"""
"""
function assemble_vector(a::Assembler,vecdata)
  b = allocate_vector(a,vecdata)
  assemble_vector!(b,a,vecdata)
  b
end

"""
"""
function assemble_matrix_and_vector(a::Assembler,data)
  A, b = allocate_matrix_and_vector(a,data)
  assemble_matrix_and_vector!(A,b,a,data)
  (A, b)
end

"""
"""
function test_assembler(a::Assembler,matdata,vecdata,data)
  A = allocate_matrix(a,matdata)
  @test num_cols(a) == size(A,2)
  @test num_rows(a) == size(A,1)
  assemble_matrix!(A,a,matdata)
  assemble_matrix_add!(A,a,matdata)
  A = assemble_matrix(a,matdata)
  @test num_cols(a) == size(A,2)
  @test num_rows(a) == size(A,1)
  b = allocate_vector(a,vecdata)
  @test num_rows(a) == length(b)
  assemble_vector!(b,a,vecdata)
  assemble_vector_add!(b,a,vecdata)
  b = assemble_vector(a,vecdata)
  @test num_rows(a) == length(b)
  A, b = allocate_matrix_and_vector(a,data)
  assemble_matrix_and_vector!(A,b,a,data)
  assemble_matrix_and_vector_add!(A,b,a,data)
  @test num_cols(a) == size(A,2)
  @test num_rows(a) == size(A,1)
  @test num_rows(a) == length(b)
  A, b = assemble_matrix_and_vector(a,data)
  @test num_cols(a) == size(A,2)
  @test num_rows(a) == size(A,1)
  @test num_rows(a) == length(b)
end

# Some syntactic sugar for assembling from anonymous functions
# and objects from which one can collect cell matrices/vectors

function assemble_matrix(f::Function,a::Assembler,U::FESpace,V::FESpace)
  v = get_fe_basis(V)
  u = get_trial_fe_basis(U)
  assemble_matrix(a,collect_cell_matrix(U,V,f(u,v)))
end

function assemble_matrix!(f::Function,A::AbstractMatrix,a::Assembler,U::FESpace,V::FESpace)
  v = get_fe_basis(V)
  u = get_trial_fe_basis(U)
  assemble_matrix!(A,a,collect_cell_matrix(U,V,f(u,v)))
end

function assemble_vector(f::Function,a::Assembler,V::FESpace)
  v = get_fe_basis(V)
  assemble_vector(a,collect_cell_vector(V,f(v)))
end

function assemble_vector!(f::Function,b::AbstractVector,a::Assembler,V::FESpace)
  v = get_fe_basis(V)
  assemble_vector!(b,a,collect_cell_vector(V,f(v)))
end

function assemble_matrix_and_vector(f::Function,b::Function,a::Assembler,U::FESpace,V::FESpace)
  v = get_fe_basis(V)
  u = get_trial_fe_basis(U)
  assemble_matrix_and_vector(a,collect_cell_matrix_and_vector(U,V,f(u,v),b(v)))
end

function assemble_matrix_and_vector!(f::Function,b::Function,M::AbstractMatrix,r::AbstractVector,a::Assembler,U::FESpace,V::FESpace)
  v = get_fe_basis(V)
  u = get_trial_fe_basis(U)
  assemble_matrix_and_vector!(M,r,a,collect_cell_matrix_and_vector(U,V,f(u,v),b(v)))
end

function assemble_matrix(f,a::Assembler,U::FESpace,V::FESpace)
  assemble_matrix(a,collect_cell_matrix(U,V,f))
end

function assemble_matrix!(f,A::AbstractMatrix,a::Assembler,U::FESpace,V::FESpace)
  assemble_matrix!(A,a,collect_cell_matrix(U,V,f))
end

function assemble_vector(f,a::Assembler,V::FESpace)
  assemble_vector(a,collect_cell_vector(V,f))
end

function assemble_vector!(f,b::AbstractVector,a::Assembler,V::FESpace)
  assemble_vector!(b,a,collect_cell_vector(V,f))
end

function assemble_matrix_and_vector(f,b,a::Assembler,U::FESpace,V::FESpace)
  assemble_matrix_and_vector(a,collect_cell_matrix_and_vector(U,V,f,b))
end

function assemble_matrix_and_vector!(f,b,M::AbstractMatrix,r::AbstractVector,a::Assembler,U::FESpace,V::FESpace)
  assemble_matrix_and_vector!(M,r,a,collect_cell_matrix_and_vector(U,V,f,b))
end

function assemble_matrix(f,U::FESpace,V::FESpace)
  a = SparseMatrixAssembler(U,V)
  assemble_matrix(f,a,U,V)
end

function assemble_matrix!(f,A::AbstractMatrix,U::FESpace,V::FESpace)
  a = SparseMatrixAssembler(U,V)
  assemble_matrix!(f,A,a,U,V)
end

function assemble_vector(f,V::FESpace)
  a = SparseMatrixAssembler(V,V)
  assemble_vector(f,a,V)
end

function assemble_vector!(f,b::AbstractVector,V::FESpace)
  a = SparseMatrixAssembler(V,V)
  assemble_vector!(f,b,a,V)
end

function assemble_matrix_and_vector(f,b,U::FESpace,V::FESpace)
  a = SparseMatrixAssembler(U,V)
  assemble_matrix_and_vector(f,b,a,U,V)
end

function assemble_matrix_and_vector!(f,b,M::AbstractMatrix,r::AbstractVector,U::FESpace,V::FESpace)
  a = SparseMatrixAssembler(U,V)
  assemble_matrix_and_vector!(f,b,M,r,a,U,V)
end

# Abstract interface for computing the data to be sent to the assembler

function collect_cell_matrix(trial::FESpace,test::FESpace,mat_contributions)
  @abstractmethod
end

function collect_cell_vector(test::FESpace,vec_contributions)
  @abstractmethod
end

function collect_cell_matrix_and_vector(
  trial::FESpace,test::FESpace,mat_contributions,vec_contributions)
  @abstractmethod
end

function collect_cell_matrix_and_vector(
  trial::FESpace,test::FESpace,mat_contributions,vec_contributions,uhd::FEFunction)
  @abstractmethod
end

# Implementation of this interface for DomainContribution

function collect_cell_matrix(trial::FESpace,test::FESpace,a::DomainContribution)
  w = []
  r = []
  c = []
  for strian in get_domains(a)
    scell_mat = get_contribution(a,strian)
    cell_mat, trian = move_contributions(scell_mat,strian)
    @assert ndims(eltype(cell_mat)) == 2
    cell_mat_c = attach_constraints_cols(trial,cell_mat,trian)
    cell_mat_rc = attach_constraints_rows(test,cell_mat_c,trian)
    rows = get_cell_dof_ids(test,trian)
    cols = get_cell_dof_ids(trial,trian)
    #push!(w,compress_contributions(cell_mat_rc,trian))
    #push!(r,compress_ids(rows,trian))
    #push!(c,compress_ids(cols,trian))
    push!(w,cell_mat_rc)
    push!(r,rows)
    push!(c,cols)
  end
  (w,r,c)
end

function collect_cell_vector(test::FESpace,a::DomainContribution)
  w = []
  r = []
  for strian in get_domains(a)
    scell_vec = get_contribution(a,strian)
    cell_vec, trian = move_contributions(scell_vec,strian)
    @assert ndims(eltype(cell_vec)) == 1
    cell_vec_r = attach_constraints_rows(test,cell_vec,trian)
    rows = get_cell_dof_ids(test,trian)
    #push!(w,compress_contributions(cell_vec_r,trian))
    #push!(r,compress_ids(rows,trian))
    push!(w,cell_vec_r)
    push!(r,rows)
  end
  (w,r)
end

function _collect_cell_matvec(trial::FESpace,test::FESpace,a::DomainContribution)
  w = []
  r = []
  c = []
  for strian in get_domains(a)
    scell_mat = get_contribution(a,strian)
    cell_mat, trian = move_contributions(scell_mat,strian)
    @assert eltype(cell_mat) <: Tuple
    cell_mat_c = attach_constraints_cols(trial,cell_mat,trian)
    cell_mat_rc = attach_constraints_rows(test,cell_mat_c,trian)
    rows = get_cell_dof_ids(test,trian)
    cols = get_cell_dof_ids(trial,trian)
    #push!(w,compress_contributions(cell_mat_rc,trian))
    #push!(r,compress_ids(rows,trian))
    #push!(c,compress_ids(cols,trian))
    push!(w,cell_mat_rc)
    push!(r,rows)
    push!(c,cols)
  end
  (w,r,c)
end

function collect_cell_matrix_and_vector(
  trial::FESpace,test::FESpace,
  biform::DomainContribution,liform::DomainContribution)

  matvec, mat, vec = _pair_contribution_when_possible(biform,liform)
  matvecdata = _collect_cell_matvec(trial,test,matvec)
  matdata = collect_cell_matrix(trial,test,mat)
  vecdata = collect_cell_vector(test,vec)
  (matvecdata, matdata, vecdata)
end

function collect_cell_matrix_and_vector(
  trial::FESpace,test::FESpace,
  biform::DomainContribution,liform::DomainContribution,uhd::FEFunction)

  matvec, mat, vec = _pair_contribution_when_possible(biform,liform,uhd)

  matvecdata = _collect_cell_matvec(trial,test,matvec)
  matdata = collect_cell_matrix(trial,test,mat)
  vecdata = collect_cell_vector(test,vec)
  (matvecdata, matdata, vecdata)
end

function _pair_contribution_when_possible(biform,liform)
  matvec = DomainContribution()
  mat = DomainContribution()
  vec = DomainContribution()
  for (trian,t) in biform.dict
    if haskey(liform.dict,trian)
      matvec.dict[trian] = pair_arrays(t,liform.dict[trian])
    else
      mat.dict[trian] = t
    end
  end
  for (trian,t) in liform.dict
    if ! haskey(biform.dict,trian)
      vec.dict[trian] = t
    end
  end
  matvec, mat, vec
end

function _pair_contribution_when_possible(biform,liform,uhd)
  _matvec, _mat, _vec = _pair_contribution_when_possible(biform,liform)
  matvec = DomainContribution()
  mat = DomainContribution()
  for (trian,t) in _matvec.dict
    cellvals = get_cell_dof_values(uhd,trian)
    cellmask = get_cell_is_dirichlet(uhd,trian)
    matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
  end
  for (trian,t) in _mat.dict
    cellvals = get_cell_dof_values(uhd,trian)
    cellmask = get_cell_is_dirichlet(uhd,trian)
    matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)
  end
  matvec, mat, _vec
end

# allow linear forms like `l(v) = 0

function collect_cell_vector(test::FESpace,l::Number)
  @notimplementedif l != 0
  w = []
  r = []
  (w,r)
end

function collect_cell_matrix_and_vector(
  trial::FESpace,test::FESpace,mat_contributions::DomainContribution,l::Number)
  @notimplementedif l != 0
  vec_contributions = DomainContribution()
  collect_cell_matrix_and_vector(test,trial,mat_contributions,vec_contributions)
end

function collect_cell_matrix_and_vector(
  trial::FESpace,test::FESpace,biform::DomainContribution,l::Number,uhd::FEFunction)
  liform = DomainContribution()
  collect_cell_matrix_and_vector(test,trial,biform,liform,uhd)
end


text:
using Epidemia
using Dates: DateTime, Day
using Test

using ReverseDiff
import Zygote
using Optim
using StatsPlots
using Distributions
using ParameterHandling
using ParameterHandling: value
using UnPack
using TransformVariables
using LogDensityProblems
using AdvancedHMC
using DynamicHMC

# @testset "Epidemia.jl" begin
struct RandomWalkProblem2{Ty, Tx, Trt, Tσ}
    y::Ty
    x::Tx
    epirt::Trt
    σ::Tσ
end

p = let
    d = dates = DateTime("2020-04-01"):Day(1):DateTime("2021-01-01") |> collect
    σ = 1e-6
    T = 100
    x = 1:length(d)
    y = sin.(2π/T*x) + rand(Normal(0,σ), length(x))

    epirt = RandomWalkModel(
        d,
        r0_prior = Normal(0.,0.1),
        σ_rw_prior = InverseGamma(3., .01),
        link = identity,
        invlink = identity
    )

    RandomWalkProblem2( y, d, epirt, σ )
end

plot(p.x, p.y)

## ===========================================================================

function (p::RandomWalkProblem2)(θ)
    @unpack y, epirt, σ = p
    @unpack r0, rw = θ

    ℓ  = 0.
    ℓ += logpdf(epirt, θ)
    ŷ = vcat(r0, rw)
    ℓ += logpdf(MvNormal(ŷ, σ), y)
    return ℓ
end
Distributions.logpdf(p::RandomWalkProblem2, θ) = p(θ)
## ===========================================================================
# Estimate parameters with TransformVariables + LogDensityProblems + Optim

function problem_transformation(p::RandomWalkProblem2)
    as(( σ_rw = asℝ₊, r0 = asℝ, rw = as(Array, length(p.y)-1)))
end

t = problem_transformation(p)
θ = rand(epirt)
x = TransformVariables.inverse(t, θ)
P = TransformedLogDensity(t, p)

∇P = ADgradient(Val(:Zygote), P)
loss(x) = -LogDensityProblems.logdensity(P, x)
grad(x) = -LogDensityProblems.logdensity_and_gradient(∇P, x)[2]

opt = optimize(loss, grad, x0; method=LBFGS(), inplace=false)

x = opt.minimizer
θ = TransformVariables.transform(t, x)
ŷ = vcat(vcat(θ.r0, θ.rw))
plot!(p.x, ŷ)

## ===========================================================================
# Estimate parameters with ParameterHandling + DynamicHMC
using Random
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000; reporter = ProgressMeterReport())

θs = transform.(t, results.chain)
plot( mean( getindex.(θs, :rw) ) )
@show DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics)
## ===========================================================================
# Estimate parameters with ParameterHandling + Optim

epirt_params = rand(epirt)
t = (
    σ_rw = positive(epirt_params.σ_rw),
    r0   = epirt_params.r0,
    rw   = epirt_params.rw
)
x0, unflatten = flatten( t )
unpack(x) = value(unflatten(x))
θ0 = unpack(x0)

loss(x) = -logpdf(p, unpack(x))
grad(x) = only(Zygote.gradient( loss, x ))

opt = optimize(loss, grad, x0; method=LBFGS(), inplace=false)

x = opt.minimizer
θ = unpack(x)
ŷ = vcat(vcat(θ.r0, θ.rw))
plot!(p.x, ŷ)

## ===========================================================================
# Estimate parameters with ParameterHandling + AdvancedHMC

Distributions.logpdf(p::RandomWalkProblem2, x::AbstractArray) =
    Distributions.logpdf(p, unpack(x))

function _logdensity_getter(p, x)
    logdensity(x) = logpdf(p, x)
end
logdensity = _logdensity_getter(p, x)

function _logdensity_and_gradient_getter(p, x)
    function logdensity_and_gradient(x)
        logdensity, back = Zygote.pullback(x->logpdf(p, x), x)
        ∂θ = first(back(1.0))
        return logdensity, ∂θ
    end
end

logdensity_and_gradient = _logdensity_and_gradient_getter(p, x)

# Sampling parameter settings
n_samples, n_adapts = 2000, 1000

# Define metric space, Hamiltonian, sampling method and adaptor
metric      = DiagEuclideanMetric( length(x) )
hamiltonian = Hamiltonian(metric, logdensity, logdensity_and_gradient)
initial_ϵ   = find_good_stepsize(hamiltonian, x0)
integrator  = Leapfrog(initial_ϵ)
proposal    = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
adaptor     = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

# Perform inference.
samples, stats = sample(hamiltonian, proposal, x0, n_samples, adaptor, n_adapts; progress=true)
θs = unpack.(samples)
plot( mean( getindex.(θs, :rw) ) )

## Plotting / inspection utilities from MCMCChains
using MCMCChains
using OrderedCollections

(t,_) = iterate(ts)
"Converts a vector of named tuples to a tuple of vectors."
function _vectup2tupvec(ts::AbstractVector{<:NamedTuple})
    ks = keys(first(ts))
    t  = first(ts)
    dic= OrderedDict()
    for t in ts
        for (k,v) in zip(keys(t), values(t))
            if isa(v, Number)
                if k in keys(dic)
                    push!(dic[k], v)
                else
                    dic[k] = [v]
                end
            else
                for (i, x) in enumerate(v)
                    kx = Symbol("$k[$i]")
                    if kx in keys(dic)
                        push!(dic[kx], x)
                    else
                        dic[kx] = [x]
                    end
                end
            end
        end
    end
    return (; dic...)
end

function _tupvec2chain(tv)
    nsamples = length(first(tv))
    nparams  = length(tv)
    ks = keys(tv)
    d = Array{Float64}(undef, nsamples, nparams, 1)
    n = Vector{Symbol}(undef, nparams)
    for (i, k) in enumerate(ks)
        d[:, i, 1] = tv[k]
        n[i] = k
    end
    MCMCChains.Chains(d, n)
end

function samples2chain(s)
    tv = vectup2tupvec(s)
    return tupvec2chain(tv)
end

n = filter( x->!occursin(r"\[", x), String.(names(ch)))
p = plot(ch[n])

## =========================================================================
# integrate into Turing

# using Turing
# @model function turingmodel(p)
#     θ ~ p
#     Turing.@addlogprob! logpdf(p, θ)
# end
#
# m = turingmodel(p)
# m()


text:
using StructuresKit

#Confirm that VX is the first moment area about the x axis in the CrossSection.CUFSMProperties function.

#Use Beer and Johnston Mechanics of Materials, Example 6.06

#Cee section without lips

CorZ = 1
t = 0.15
H = 6.0
B = 4.0
D = 0.0
θ = 0.0

r = 0.0
kipin = 0
center = 1
n = 4

prop,node,elem,lengths,springs,constraints,geom,cz = CrossSection.CUFSMtemplate(CorZ,H,B,B,D,D,r,r,r,r,θ,θ,t,n,n,n,n,n,n,n,n,n,kipin,center)

using Plots
plot(node[:,2], node[:,3], seriestype=:scatter)

A, xc, yc, Ixc, Iyc, Ixyc, Imax, Imin, Th_p, Cw, J, Xs, Ys, w, Bx, By, B1, B2, VX, VY = CrossSection.CUFSMproperties(node, elem)

#Stress at bottom flange-web intersection is 2.22 ksi (from left to right) for V=2.5 kips down.

#VX is the same as the first moment area taken about the x axis, call it Qx

τ_solution = 2.22

Q = VX[5]
V = -2.5
τ =  V*Q/(Ixc*t)


text:
"""
Abstract type for all kinds of gene product associations
"""
abstract type GeneProductAssociation end

"""
Gene product reference in the association expression
"""
struct GPARef <: GeneProductAssociation
    gene_product::String
end

"""
Boolean binary "and" in the association expression
"""
struct GPAAnd <: GeneProductAssociation
    terms::Vector{GeneProductAssociation}
end

"""
Boolean binary "or" in the association expression
"""
struct GPAOr <: GeneProductAssociation
    terms::Vector{GeneProductAssociation}
end

"""
A simplified representation of MathML-specified math AST
"""
abstract type Math end

"""
A literal value (usually a numeric constant) in mathematical expression
"""
struct MathVal{T} <: Math where {T}
    val::T
end

"""
An identifier (usually a variable name) in mathematical expression
"""
struct MathIdent <: Math
    id::String
end

"""
A constant identified by name (usually something like `pi`, `e` or `true`) in
mathematical expression
"""
struct MathConst <: Math
    id::String
end

"""
A special value representing the current time of the simulation, with a special
name.
"""
struct MathTime <: Math
    id::String
end

"""
Function application ("call by name", no tricks allowed) in mathematical expression
"""
struct MathApply <: Math
    fn::String
    args::Vector{Math}
end

"""
Function definition (aka "lambda") in mathematical expression
"""
struct MathLambda <: Math
    args::Vector{String}
    body::Math
end

"""
SBML Compartment with sizing information.
"""
struct Compartment
    name::Maybe{String}
    constant::Maybe{Bool}
    spatial_dimensions::Maybe{Int}
    size::Maybe{Float64}
    units::Maybe{String}
    notes::Maybe{String}
    annotation::Maybe{String}
    Compartment(na, c, sd, s, u, no = nothing, an = nothing) = new(na, c, sd, s, u, no, an)
end

"""
Reaction with stoichiometry that assigns reactants and products their relative
consumption/production rates (accessible in field `stoichiometry`), lower/upper
bounds (in tuples `lb` and `ub`, with unit names), and objective coefficient
(`oc`). Also may contains `notes` and `annotation`.
"""
struct Reaction
    stoichiometry::Dict{String,Float64}
    lb::Tuple{Float64,String}
    ub::Tuple{Float64,String}
    oc::Float64
    gene_product_association::Maybe{GeneProductAssociation}
    kinetic_math::Maybe{Math}
    reversible::Bool
    notes::Maybe{String}
    annotation::Maybe{String}
    Reaction(s, l, u, o, as, km, r, n = nothing, an = nothing) =
        new(s, l, u, o, as, km, r, n, an)
end

"""
Species metadata -- contains a human-readable `name`, a `compartment`
identifier, `formula`, `charge`, and additional `notes` and `annotation`.
"""
struct Species
    name::Maybe{String}
    compartment::String
    boundary_condition::Maybe{Bool}
    formula::Maybe{String}
    charge::Maybe{Int}
    initial_amount::Maybe{Tuple{Float64,Maybe{String}}}
    initial_concentration::Maybe{Tuple{Float64,Maybe{String}}}
    only_substance_units::Maybe{Bool}
    notes::Maybe{String}
    annotation::Maybe{String}
    Species(na, co, b, f, ch, ia, ic, osu, no = nothing, a = nothing) =
        new(na, co, b, f, ch, ia, ic, osu, no, a)
end

"""
Gene product metadata.
"""
struct GeneProduct
    name::Maybe{String}
    label::Maybe{String}
    notes::Maybe{String}
    annotation::Maybe{String}
    GeneProduct(na, l, no = nothing, a = nothing) = new(na, l, no, a)
end

"""
Custom function definition.
"""
struct FunctionDefinition
    name::Maybe{String}
    body::Maybe{Math}
    notes::Maybe{String}
    annotation::Maybe{String}
    FunctionDefinition(na, b, no = nothing, a = nothing) = new(na, b, no, a)
end

"""
Structure that collects the model-related data. Contains `parameters`, `units`,
`compartments`, `species` and `reactions` and `gene_products`, and additional
`notes` and `annotation` (also present internally in some of the data fields).
The contained dictionaries are indexed by identifiers of the corresponding
objects.
"""
struct Model
    parameters::Dict{String,Float64}
    units::Dict{String,Number}
    compartments::Dict{String,Compartment}
    species::Dict{String,Species}
    reactions::Dict{String,Reaction}
    gene_products::Dict{String,GeneProduct}
    function_definitions::Dict{String,FunctionDefinition}
    notes::Maybe{String}
    annotation::Maybe{String}
    Model(p, u, c, s, r, g, f, n = nothing, a = nothing) = new(p, u, c, s, r, g, f, n, a)
end


text:
"""
    qrdelete!(Q, R, k)
Delete the left-most column of F = Q[:, 1:k] * R[1:k, 1:k] by updating Q and R.
Only Q[:, 1:(k-1)] and R[1:(k-1), 1:(k-1)] are valid on exit.
"""
function qrdelete!(Q::AbstractMatrix, R::AbstractMatrix, k::Int)
  n, m = size(Q)
  m == checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  # apply Givens rotations
  for i in 2:k
      g = first(givens(R, i - 1, i, i))
      lmul!(g, R)
      rmul!(Q, g')
  end

  # move columns of R
  @inbounds for j in 1:(k-1)
    for i in 1:(k-1)
      R[i, j] = R[i, j + 1]
    end
  end

  Q, R
end

"""
    qradd!(Q, R, v, k)
Replace the right-most column of F = Q[:, 1:k] * R[1:k, 1:k] with v by updating Q and R.
This implementation modifies vector v as well. Only Q[:, 1:k] and R[1:k, 1:k] are valid on
exit.
"""
function qradd!(Q::AbstractMatrix, R::AbstractMatrix, v::AbstractVector, k::Int)
  n, m = size(Q)
  n == length(v) || throw(DimensionMismatch())
  m == checksquare(R) || throw(DimensionMismatch())
  1 ≤ k ≤ m || throw(ArgumentError())

  @inbounds for i in 1:(k-1)
    q = view(Q, :, i)
    r = dot(q, v)

    R[i, k] = r
    axpy!(-r, q, v)
  end

  @inbounds begin
    d = norm(v)
    R[k, k] = d
    @. Q[:, k] = v / d
  end

  Q, R
end


text:
include("../../Modules/Gridcells/HexGrid.jl")

function dist(p1::Array{Float64},p2::Array{Float64})
    d = sqrt((p1[1]-p2[1])^2+(p1[2]-p2[2])^2)
    return d
end

                                # for the intensity to index coding (never equals one)
gauss(d::Float64, node_size::Float64) = 0.999999 * exp(-2*(d/node_size)) # considered 2sigma as d in the gaussian

function input_spikes(obs::Array{Float64}, cells_pos::Array{HexModule}, cells_vel::Array{HexModule})
    pos = [obs[1],obs[3]%pi]     # obs pos and theta
    vel = [obs[2],obs[4]]     # obs pos_v and theta_v
    n_module = length(cells_pos)
    spikes = Tuple{Int64,Float64}[]
    neuron_tot = sum([cells_pos[i].n_grid for i in 1:n_module])
    neuron_passed = 0
    for i in 1:n_module
        d_pos = cells_pos[i].grids[1].node_size
        d_vel = cells_vel[i].grids[1].node_size
        for j in 1:cells_pos[i].n_grid           # parameter
            close_pos = filter(p -> dist(p,pos) < 1.5*d_pos, cells_pos[i].grids[j].grid) # 1.5 is to reach up to 3sigmas
            close_vel =  filter(p -> dist(p,vel) < 1.5*d_vel, cells_vel[i].grids[j].grid)
            if length(close_pos) != 0
                neuron = neuron_passed + j
                d = dist(pos,close_pos[1])
                intensity = gauss(d,d_pos)
                push!(spikes,(neuron,intensity))
            end
            if length(close_vel) != 0
                neuron = neuron_tot + neuron_passed + j
                d = dist(vel,close_vel[1])
                intensity = gauss(d,d_vel)
                push!(spikes,(neuron,intensity))
            end
        end
        neuron_passed += cells_pos[i].n_grid
    end
    return spikes
end

function latency_spikes(obs::Array{Float64}, cells_pos::Array{HexModule}, cells_vel::Array{HexModule}, win_size::Int64)
    pos = [obs[1],obs[3]%pi]     # obs pos and theta
    vel = [obs[2],obs[4]]     # obs pos_v and theta_v
    n_module = length(cells_pos)
    spikes = [Int64[] for _ in 1:win_size]
    neuron_tot = sum([cells_pos[i].n_grid for i in 1:n_module])
    neuron_passed = 0
    for i in 1:n_module
        d_pos = cells_pos[i].grids[1].node_size
        d_vel = cells_vel[i].grids[1].node_size
        for j in 1:cells_pos[i].n_grid           # parameter
            close_pos = filter(p -> dist(p,pos) < 1.25*d_pos, cells_pos[i].grids[j].grid) # 1.5 is to reach up to 3sigmas
            close_vel =  filter(p -> dist(p,vel) < 1.25*d_vel, cells_vel[i].grids[j].grid)
            if length(close_pos) != 0
                neuron = neuron_passed + j
                d = dist(pos,close_pos[1])
                intensity = gauss(d,d_pos)
                ind = win_size - floor(Int64, win_size * intensity)
                push!(spikes[ind],neuron)
            end
            if length(close_vel) != 0
                neuron = neuron_tot + neuron_passed + j
                d = dist(vel,close_vel[1])
                intensity = gauss(d,d_vel)
                ind = win_size - floor(Int64, win_size * intensity)
                push!(spikes[ind],neuron)
            end
        end
        neuron_passed += cells_pos[i].n_grid
    end
    return spikes
end

function input_spikes(obs::Array{Float64}, grid_cells::Array{HexGrid})
    pos = [obs[1],obs[3]%pi]     # obs pos and theta
    vel = [obs[2],obs[4]]     # obs pos_v and theta_v
    bound = [-15.0, 15.0]
    pos_scaled = (pos .- bound[1]) / (bound[2]-bound[1])
    vel_scaled = (vel .- bound[1]) / (bound[2]-bound[1])
    n_grid = length(grid_cells)
    spikes = Tuple{Int64,Float64}[]
    for i in 1:n_grid
        d = grid_cells[i].node_size
        close_pos = filter(p -> dist(p,pos_scaled) < 1.5*d, grid_cells[i].grid) # 1.5 is to reach up to 3sigmas
        close_vel =  filter(p -> dist(p,vel_scaled) < 1.5*d, grid_cells[i].grid)
        if length(close_pos) != 0
            neuron = i
            d_pos = dist(pos_scaled,close_pos[1])
            intensity = gauss(d_pos,d)
            push!(spikes,(neuron,intensity))
        end
        if length(close_vel) != 0
            neuron = n_grid + i
            d_vel = dist(vel_scaled,close_vel[1])
            intensity = gauss(d_vel,d)
            push!(spikes,(neuron,intensity))
        end
    end
    return spikes
end

function latency_spikes(obs::Array{Float64}, grid_cells::Array{HexGrid}, win_size::Int64)
    pos = [obs[1],obs[3]%pi]     # obs pos and theta
    vel = [obs[2],obs[4]]     # obs pos_v and theta_v
    bound = [-15.0, 15.0]
    pos_scaled = (pos .- bound[1]) / (bound[2]-bound[1])
    vel_scaled = (vel .- bound[1]) / (bound[2]-bound[1])
    n_grid = length(grid_cells)
    spikes = [Int64[] for _ in 1:win_size]
    for i in 1:n_grid
        d = grid_cells[i].node_size
        close_pos = filter(p -> dist(p,pos_scaled) < 1.25*d, grid_cells[i].grid) # 1.5 is to reach up to 3sigmas
        close_vel =  filter(p -> dist(p,vel_scaled) < 1.25*d, grid_cells[i].grid)
        if length(close_pos) != 0
            neuron = i
            d_pos = dist(pos_scaled,close_pos[1])
            intensity = gauss(d_pos,d)
            ind = win_size - floor(Int64, win_size * intensity)
            push!(spikes[ind],neuron)
        end
        if length(close_vel) != 0
            neuron = n_grid + i
            d_vel = dist(vel_scaled,close_vel[1])
            intensity = gauss(d_vel,d)
            ind = win_size - floor(Int64, win_size * intensity)
            push!(spikes[ind],neuron)
        end
    end
    return spikes
end

# # %% test input_spikes
#
# cells_tmp_pos = HexModule[]
# cells_tmp_vel = HexModule[]
#
#
# borpos = [-4.1,4.1]
# borvpos = [-50.0,50.0]
# bortheta = [-pi,pi]
# borvtheta = [-50.0,50.0]
#
# for neuron in 1:16 # parameter
#     res = 0.07 # parameter
#     dilat = rand()*(1.4-0.7)+0.7 # parameter
#     theta = rand()*(2pi/6)-pi/6
#     mod_pos = HexModule(res,dilat = dilat, theta = theta, bor_x = borpos, bor_y = bortheta)
#     mod_vel = HexModule(res,dilat = dilat, theta = theta, bor_x = borvpos, bor_y = borvtheta)
#     push!(cells_tmp_pos,mod_pos)
#     push!(cells_tmp_vel,mod_vel)
# end
#
# const cells_pos = cells_tmp_pos
# const cells_vel = cells_tmp_vel
#
# # %%
#
# obs = [0.523324424,0.323242424243,0.4242424225,0.224253525253]
#
#
# @time spikes = input_spikes(obs,cells_pos,cells_vel)

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
module TestExpressionsTVL

import Test: @testset, @test
import Volcanito: tvl, @tvl

@testset "TVL || truth table" begin
    @test @tvl(true || true) === true
    @test @tvl(true || false) === true
    @test @tvl(true || missing) === true
    @test @tvl(false || true) === true
    @test @tvl(false || false) === false
    @test @tvl(false || missing) === missing
    @test @tvl(missing || true) === true
    @test @tvl(missing || false) === missing
    @test @tvl(missing || missing) === missing
end

@testset "TVL && truth table" begin
    @test @tvl(true && true) === true
    @test @tvl(true && false) === false
    @test @tvl(true && missing) === missing
    @test @tvl(false && true) === false
    @test @tvl(false && false) === false
    @test @tvl(false && missing) === false
    @test @tvl(missing && true) === missing
    @test @tvl(missing && false) === false
    @test @tvl(missing && missing) === missing
end

# We define a function that prints out a unique ID for each argument
# to a Boolean operator. By wrapping all Boolean values in calls to this
# function, we're able to check that the order of evaluation and
# side-effects of the short-circuiting operators are retained by our
# macro rewrites.

function f(io, i, x)
    print(io, i)
    x
end

@testset "Order of evaluation for tvl" begin
    for x in (true, false)
        for y in (true, false)
            for z in (true, false)
                io = IOBuffer()

                a = f(io, 1, x) && f(io, 2, y) || f(io, 3, z)
                order_a = String(take!(io))
                b = @tvl f(io, 1, x) && f(io, 2, y) || f(io, 3, z)
                order_b = String(take!(io))

                @test a === b
                @test order_a === order_b

                a = f(io, 1, x) || f(io, 2, y) && f(io, 3, z)
                order_a = String(take!(io))
                b = @tvl f(io, 1, x) || f(io, 2, y) && f(io, 3, z)
                order_b = String(take!(io))

                @test a === b
                @test order_a === order_b
            end
        end
    end
end

@testset "@tvl with local variables" begin
    let x = 1, y = 2
        @test (@tvl x > 1 || y < 2) === false
    end
end

end


text:
using SpecialFunctions
using DiffRules: diffrule, hasdiffrule

@testset "Functional" begin
    let rng = MersenneTwister(123456)
        import Nabla.fmad

        # Check that `broadcastsum!` works as intended.
        let
            Z, X, Y = [1 2; 3 4], [1 2], [5 6; 7 8]
            @test Nabla.broadcastsum!(x->2x, false, copy(X), Z) == [2 + 6 4 + 8]
            @test Nabla.broadcastsum!(x->2x, true, copy(X), Z) == X + [2 + 6 4 + 8]
            @test Nabla.broadcastsum!(x->2x, false, copy(Y), Z) == 2Z
            @test Nabla.broadcastsum!(x->2x, true, copy(Y), Z) == Y + 2Z
        end

        # Check that `broadcast` returns the correct gradient under the defined unary functions.
        function check_unary_broadcast(f, x)
            x_ = Leaf(Tape(), x)
            s = broadcast(f, x_)
            return ∇(s, oneslike(unbox(s)))[x_] ≈ ∇.(f, Arg{1}, x)
        end
        for (package, f) in Nabla.unary_sensitivities
            domain = domain1(eval(f))
            domain === nothing && error("Could not determine domain for $f.")
            x_dist = Uniform(domain...)
            x = rand(rng, x_dist, 100)
            @test check_unary_broadcast(eval(f), x)
        end

        # Check that `broadcast` returns the correct gradient under each implemented binary
        # function.
        function check_binary_broadcast(f, x, y)
            tape = Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            s = broadcast(f, x_, y_)
            o = oneslike(unbox(s))
            ∇s = ∇(s, o)
            ∇x = broadcast((z, z̄, x, y)->∇(f, Arg{1}, nothing, z, z̄, x, y),
                           unbox(s), o, x, y)
            ∇y = broadcast((z, z̄, x, y)->∇(f, Arg{2}, nothing, z, z̄, x, y),
                           unbox(s), o, x, y)
            @test broadcast(f, x, y) == unbox(s)
            @test ∇s[x_] ≈ ∇x
            @test ∇s[y_] ≈ ∇y
        end
        function check_binary_broadcast(f, x::Real, y)
            tape = Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            s = broadcast(f, x_, y_)
            o = oneslike(unbox(s))
            ∇s = ∇(s, o)
            ∇x = sum(broadcast((z, z̄, x, y)->∇(f, Arg{1}, nothing, z, z̄, x, y),
                               unbox(s), o, x, y))
            ∇y = broadcast((z, z̄, x, y)->∇(f, Arg{2}, nothing, z, z̄, x, y),
                           unbox(s), o, x, y)
            @test broadcast(f, x, y) == unbox(s)
            @test ∇s[x_] ≈ ∇x
            @test ∇s[y_] ≈ ∇y
        end
        function check_binary_broadcast(f, x, y::Real)
            tape = Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            s = broadcast(f, x_, y_)
            o = oneslike(unbox(s))
            ∇s = ∇(s, o)
            ∇x = broadcast((z, z̄, x, y)->∇(f, Arg{1}, nothing, z, z̄, x, y),
                           unbox(s), o, x, y)
            ∇y = sum(broadcast((z, z̄, x, y)->∇(f, Arg{2}, nothing, z, z̄, x, y),
                               unbox(s), o, x, y))
            @test broadcast(f, x, y) == unbox(s)
            @test ∇s[x_] ≈ ∇x
            @test ∇s[y_] ≈ ∇y
        end
        for (package, f) in Nabla.binary_sensitivities

            # TODO: More care needs to be taken to test the following.
            if hasdiffrule(package, f, 2)
                ∂f∂x, ∂f∂y = diffrule(package, f, :x, :y)
            else
                ∂f∂x, ∂f∂y = :∂f∂x, :∂f∂y
            end

            # TODO: Implement the edge cases for functions differentiable in only either
            # argument.
            (∂f∂x == :NaN || ∂f∂y == :NaN) && continue
            domain = domain2(eval(f))
            domain === nothing && error("Could not determine domain for $f.")
            (x_lb, x_ub), (y_lb, y_ub) = domain
            x_dist, y_dist = Uniform(x_lb, x_ub), Uniform(y_lb, y_ub)
            x, y = rand(rng, x_dist, 100), rand(rng, y_dist, 100)
            check_binary_broadcast(eval(f), x, y)
            check_binary_broadcast(eval(f), rand(rng, x_dist), y)
            check_binary_broadcast(eval(f), x, rand(rng, y_dist))
        end
        #
        let # Ternary functions (because it's useful to check I guess.)
            f = (x, y, z)->x * y + y * z + x * z
            x, y, z = randn(rng, 5), randn(rng, 5), randn(rng, 5)
            x_, y_, z_ = Leaf.(Tape(), (x, y, z))
            s_ = broadcast(f, x_, y_, z_)
            ∇s = ∇(s_, oneslike(unbox(s_)))
            @test unbox(s_) == broadcast(f, x, y, z)
            @test ∇s[x_] == getindex.(broadcast((x, y, z)->fmad(f, (x, y, z)), x, y, z), 1)
            @test ∇s[y_] == getindex.(broadcast((x, y, z)->fmad(f, (x, y, z)), x, y, z), 2)
            @test ∇s[z_] == getindex.(broadcast((x, y, z)->fmad(f, (x, y, z)), x, y, z), 3)
        end

        let
            x, y, tape = 5.0, randn(rng, 5), Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            z_ = x_ .+ y_
            z2_ = broadcast(+, x_, y_)
            @test unbox(z_) == x .+ y
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
            @test ∇(z_, oneslike(unbox(z_)))[y_] == ∇(z2_, oneslike(unbox(z2_)))[y_]
        end
        let
            x, y, tape = randn(rng, 5), 5.0, Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            z_ = x_ * y_
            z2_ = broadcast(*, x_, y_)
            @test unbox(z_) == x .* y
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
            @test ∇(z_, oneslike(unbox(z_)))[y_] == ∇(z2_, oneslike(unbox(z2_)))[y_]
        end
        let
            x, y, tape = randn(rng, 5), 5.0, Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            z_ = x_ .- y_
            z2_ = broadcast(-, x_, y_)
            @test unbox(z_) == x .- y
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
            @test ∇(z_, oneslike(unbox(z_)))[y_] == ∇(z2_, oneslike(unbox(z2_)))[y_]
        end
        let
            x, y, tape = randn(rng, 5), 5.0, Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            z_ = x_ / y_
            z2_ = broadcast(/, x_, y_)
            @test unbox(z_) == x ./ y
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
            @test ∇(z_, oneslike(unbox(z_)))[y_] == ∇(z2_, oneslike(unbox(z2_)))[y_]
        end
        let
            x, y, tape = 5.0, randn(rng, 5), Tape()
            x_, y_ = Leaf(tape, x), Leaf(tape, y)
            z_ = x_ \ y_
            z2_ = broadcast(\, x_, y_)
            @test unbox(z_) == x .\ y
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
            @test ∇(z_, oneslike(unbox(z_)))[y_] == ∇(z2_, oneslike(unbox(z2_)))[y_]
        end

        # Check that dot notation works as expected for all unary function in Nabla for both
        # scalars and arrays.
        function check_unary_dot(f, x)
            x_ = Leaf(Tape(), x)
            z_ = f.(x_)
            z2_ = broadcast(f, x_)
            @test unbox(z_) == f.(x)
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
        end
        function check_unary_dot(f, x::∇Scalar)
            x_ = Leaf(Tape(), x)
            z_ = f.(x_)
            @test unbox(z_) == f.(x)
            @test ∇(z_)[x_] == ∇(broadcast(f, x_))[x_]
        end
        for (package, f) in Nabla.unary_sensitivities
            domain = domain1(eval(f))
            domain === nothing && error("Could not determine domain for $f.")
            x_dist = Uniform(domain...)
            check_unary_dot(eval(f), rand(rng, x_dist))
            check_unary_dot(eval(f), rand(rng, x_dist, 100))
        end

        # Check that the dot notation works as expected for all of the binary functions in
        # Nabla for each permutation of scalar / array input.
        function check_binary_dot(f, x, y)
            x_, y_ = Leaf.(Tape(), (x, y))
            z_ = f.(x_, y_)
            z2_ = broadcast(f, x_, y_)
            @test unbox(z_) == f.(x, y)
            @test ∇(z_, oneslike(unbox(z_)))[x_] == ∇(z2_, oneslike(unbox(z2_)))[x_]
            @test ∇(z_, oneslike(unbox(z_)))[y_] == ∇(z2_, oneslike(unbox(z2_)))[y_]
        end
        function check_binary_dot(f, x::∇Scalar, y::∇Scalar)
            x_, y_ = Leaf.(Tape(), (x, y))
            z_ = f.(x_, y_)
            @test ∇(z_)[x_] == ∇(broadcast(f, x_, y_))[x_]
            @test ∇(z_)[y_] == ∇(broadcast(f, x_, y_))[y_]
        end
        for (package, f) in Nabla.binary_sensitivities
            # TODO: More care needs to be taken to test the following.
            f in [:atan, :mod, :rem] && continue
            if hasdiffrule(package, f, 2)
                ∂f∂x, ∂f∂y = diffrule(package, f, :x, :y)
            else
                ∂f∂x, ∂f∂y = :∂f∂x, :∂f∂y
            end
            # TODO: Implement the edge cases for functions differentiable in only either
            # argument.
            (∂f∂x == :NaN || ∂f∂y == :NaN) && continue
            domain = domain2(eval(f))
            domain === nothing && error("Could not determine domain for $f.")
            (x_lb, x_ub), (y_lb, y_ub) = domain
            x_distr = Uniform(x_lb, x_ub)
            y_distr = Uniform(y_lb, y_ub)
            x = rand(rng, x_distr, 100)
            y = rand(rng, y_distr, 100)
            check_binary_dot(eval(f), x, y)
            check_binary_dot(eval(f), rand(rng, x_distr), y)
            check_binary_dot(eval(f), x, rand(rng, y_distr))
            check_binary_dot(eval(f), Ref(rand(rng, x_distr)), y)
            check_binary_dot(eval(f), x, Ref(rand(rng, y_distr)))
            check_binary_dot(eval(f), rand(rng, x_distr), rand(rng, y_distr))
        end

        # test with other broadcast styles
        let
            a = Diagonal(ones(3))
            b = ones(3, 3)
            check_binary_dot(+, a, b)
            check_binary_dot(+, b, a)
            check_binary_dot(+, a, a)
        end
    end

    # Check that the number of allocations which happen in the reverse pass of `map` and
    # `broadcast` is invariant of the size of the data structure of which we are
    # mapping / broadcasting.
    for f in [:map, :broadcast]
        let
            @eval foo_small() = sum($f(tanh, Leaf(Tape(), randn(10, 10))))
            @eval foo_large() = sum($f(tanh, Leaf(Tape(), randn(10, 100))))
            @test allocs(@benchmark foo_small()) == allocs(@benchmark foo_large())
            @test allocs(@benchmark ∇(foo_small())) == allocs(@benchmark ∇(foo_large()))
        end
    end

    # #111
    let
        f(x) = sum(Float64[1,2,3] .* (x .+ Float64[3,2,1]))
        @test ∇(f)(Float64[1,2,3]) isa Tuple{Vector{Float64}}
        @test unbox(∇(f; get_output=true)(Float64[1,2,3])[1]) == f(Float64[1,2,3])
    end

    # #117
    let
        f(x) = sum(x .+ [1,2,4] .* [4,2,1])
        @test ∇(f)(Float64[1,2,3])[1] == ones(Float64, 3)
        @test unbox(∇(f; get_output=true)(Float64[1,2,3])[1]) == f(Float64[1,2,3])
    end

    # broadcasting literal_pow
    let
        f(x) = sum(x .^ 2)
        @test ∇(f)(Float64[1,2,3])[1] == Float64[2,4,6]
        @test unbox(∇(f; get_output=true)(Float64[1,2,3])[1]) == f(Float64[1,2,3])
    end

    # fused broadcasting with different styles
    let
        f(x) = sum(Symmetric(x) .+ 0.0001 .* Diagonal(ones(size(x, 1))))
        a = rand(3, 3)
        a += transpose(a)
        @test ∇(f)(a)[1] == Float64[1 2 2; 0 1 2; 0 0 1]
        @test unbox(∇(f; get_output=true)(a)[1]) == f(a)
    end
end

struct CoolArray{T} <: AbstractVector{T}
    x::Vector{T}
end
Base.map(f, x::CoolArray, y::AbstractArray) = map(f, x.x, y)
@testset "Issue #136" begin
    # This doesn't even involve Nabla at all; we're just making sure that we aren't
    # introducing method ambiguities with our `map` overloading that can interfere with
    # other packages that extend `map`, e.g. StaticArrays
    @test map(-, CoolArray([1,2,3]), [1,2,3]) == [0,0,0]

    # `map` with lots of inputs
    for i = 1:10
        args = Any[1:3 for _ = 1:20]
        args[i] = Leaf(Tape(), 1:3)
        x = map(+, args...)
        @test x isa Branch{Vector{Int}}
        @test unbox(x) == [20,40,60]
    end
end


text:
function _file_content_hash(file)
    hash_ = hash("")
    isfile(file) && for line in eachline(file)
        hash_ = hash(line, hash_)
    end
    return hash_
end

function _file_content_mtime_event(e::CustomEvent)
    get!(extras(e), "_MTIME_") do
        FileMTimeEvent()
    end
end

function _get_file_content_hash_cache(e::CustomEvent)
    get!(extras(e), "_HASH_CACHE_") do
        _file_content_hash("")
    end
end

function _up_file_content_hash_cache!(e::CustomEvent, file)
    hash_ = _file_content_hash(file)
    setindex!(extras(e), hash_, "_HASH_CACHE_")
    return hash_
end

function _dft_file_content_old_state(e::CustomEvent, file::String)
    get!(e, file) do
        _file_content_hash(file)
    end
end

function _dft_file_content_new_state(e::CustomEvent, file::String)
    mtime_event = _file_content_mtime_event(e)
    if has_event!(mtime_event, file)
        # Recompute only is file was modify (better performance)
        # And cached
        hash_ = _up_file_content_hash_cache!(e, file)
        return hash_
    end
    return _get_file_content_hash_cache(e)
end

function _dft_file_content_update!(e::CustomEvent, file::String)
    new_hash = new_state(e, file)
    setindex!(e, new_hash, file)
    # Update also the mtime event
    mtime_event = _file_content_mtime_event(e)
    update!(mtime_event, file)
    return e
end

function _dft_file_content_trigger(::CustomEvent, old_hash, new_hash)
    return (old_hash != new_hash)
end

function FileContentEvent()
    T = Union{CustomEvent{String, Float64}, UInt64}
    return CustomEvent{String, T}(;
        old_state = _dft_file_content_old_state,
        new_state = _dft_file_content_new_state,
        update! = _dft_file_content_update!,
        trigger = _dft_file_content_trigger
    )
end


text:
#=
Tasks aka
* symmetric coroutines
* lightweight threads
* cooperative multitasking
* one-shot continuations
NOTE: tasks are not scheduled to run on separate CPU cores

Channel is a waitable first-in first-out queue
A special Channel constructor which accepts a 1-arg function as an argument
can be used to run a task bound to a channel.  Channel object is bound to the
task and automatically closes when the Task finishes.
=#
function producer(c::Channel)
    put!(c, "start")
    for n=1:4
        put!(c, 2n)
    end
    put!(c, "stop")
end

chan1 = Channel(producer)
println(take!(chan1))
println(take!(chan1))
println(take!(chan1))
println(take!(chan1))
println(take!(chan1))
println(take!(chan1))

# example using a for loop
for x in Channel(producer)
   println(x)
end

function mytask(myarg)
    println("mytask received $(myarg)")
end

# tasks expect a 0 arg function so can use an anonymous function
taskHdl = Task(() -> mytask(7))
# or, equivalently using a convenience macro
taskHdl = @task mytask(7)

max_stars_repo_name: suicidejack/julia-practice

max_issues_repo_name: suicidejack/julia-practice

max_forks_repo_name: suicidejack/julia-practice


text:
# Day 10: Syntax Scoring
# I'm out of my comfort zone here, so everything is naively implemented!

load(::Day{10}, path) = [first.(l) for l in split.(readlines(path), "")]

# Part 1 – find corrupted codes: lines where a chunk closes with the wrong character
char_score = Dict(')' => 3, ']' => 57, '}' => 1197, '>' => 25137)
match_closing = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
is_opening(c::Char) = c in "([{<"

function error_score!(line)
    expected_closing = Char[]
    while !isempty(line)
        c = popfirst!(line)
        if is_opening(c)
            pushfirst!(expected_closing, match_closing[c])
        else
            if c == first(expected_closing)
                popfirst!(expected_closing)
            else
                return char_score[c]
            end
        end
    end
    return 0 # zero error score if no error was found
end

function solve1(::Day{10}, data)
    data = deepcopy(data) # in-place `solve1!` could comment out this line for performance
    return sum(error_score!.(data))
end

# Part 2 – autocomplete correct chunks
char_score2 = Dict(')' => 1, ']' => 2, '}' => 3, '>' => 4) # new scoring rules

# we only need to modify the outputs of `error_score` from part 1
function autocomplete!(line)
    expected_closing = Char[]
    while !isempty(line)
        c = popfirst!(line)
        if is_opening(c)
            pushfirst!(expected_closing, match_closing[c])
        else
            if c == first(expected_closing)
                popfirst!(expected_closing)
            else
                return nothing # error was found
            end
        end
    end
    return expected_closing
end

function score_completion(completion)
    score = 0
    for c in completion
        score *= 5
        score += char_score2[c]
    end
    return score
end

function solve2(::Day{10}, data)
    data = deepcopy(data) # in-place `solve2!` could comment out this line for performance
    scores = score_completion.(filter(!isnothing, autocomplete!.(data))) # filter errors
    return partialsort!(scores, div(length(scores), 2, RoundUp)) # same as `median!`
end

testresult1(::Day{10}) = 26397
testresult2(::Day{10}) = 288957


text:
# powm1 is based on the Boost implementation
# https://github.com/boostorg/math/blob/242e9d669198415d165bb0213f730c00985bee9a/include/boost/math/special_functions/powm1.hpp
#
# (C) Copyright John Maddock 2006.
#
# Boost Software License - Version 1.0 - August 17th, 2003
#
# Permission is hereby granted, free of charge, to any person or organization
# obtaining a copy of the software and accompanying documentation covered by
# this license (the "Software") to use, reproduce, display, distribute,
# execute, and transmit the Software, and to prepare derivative works of the
# Software, and to permit third-parties to whom the Software is furnished to
# do so, all subject to the following:
#
# The copyright notices in the Software and this entire statement, including
# the above license grant, this restriction and the following disclaimer,
# must be included in all copies of the Software, in whole or in part, and
# all derivative works of the Software, unless such copies or derivative
# works are solely in the form of machine-executable object code generated by
# a source language processor.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
# SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
# FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
"""
    powm1(x,y)

Compute `x^y-1` accurately when `x` is near 1. Based on Boost's implementation

"""
function powm1(x::Real,y::Real)
    if x < 0 && isinteger(y/2) # y is an even integer
        x = -x
    end
    if x > 0 && (abs(y * (x - 1)) < 0.5 || abs(y) < 0.2)
        l = y * log(x)
        if l < 0.5
            return expm1(l)
        end
    end
    x^y-1
end

"""
    poscubicroot(p3,p2,p1,0)

Find the smallest positive real root of the polynomial

    p3*x^3 + p2*x^2 + p1*x + p0

"""
function poscubicroot(p3,p2,p1,p0)
    a = p2/p3
    b = p1/p3
    C = p0/p3

    Q = (a^2 - 3*b) / 9
    r = (2*a^3 - 9*a*b + 27*C) / 54

    if r^2 - Q^3 >= 0
        # 1 real root
        A = -copysign(cbrt(abs(r) + sqrt(r^2 - Q^3)), r)
        B = A==0 ? A : Q/A
        return A + B - a/3

    else
        # 3 real roots
        θ = acos(r/(Q*sqrt(Q)))

        r1 = -2 * sqrt(Q) * cos(θ / 3) - a / 3
        r2 = -2 * sqrt(Q) * cos((θ + 2π)/3) - a / 3
        r3 = -2 * sqrt(Q) * cos((θ - 2π)/3) - a / 3

        r = Inf
        if 0 < r1
            r = r1
        end
        if 0 < r2 < r
            r = r2
        end
        if 0 < r3 < r
            r = r3
        end
        if r < Inf
            return r
        else
            throw(DomainError())
        end
    end
end


text:
@testset "Pattern" begin
    manifold = IdealManifold(1575420e3, 0.1904 / 4 * SVector(SVector(1, 1, 0), SVector(-1, 1, 0), SVector(1, -1, 0), SVector(-1, -1, 0)))
    pattern = Pattern(manifold, norm, num_az = 360, num_el = 91, max_el = π / 2)
    @test first(pattern.azs) == 0.0
    @test step(pattern.azs) == 2π/360
    @test length(pattern.azs) == 360
    @test isa(pattern.azs, AbstractRange)
    @test first(pattern.els) == 0.0
    @test step(pattern.els) == π/2 / 90
    @test length(pattern.els) == 91
    @test isa(pattern.els, AbstractRange)
    @test pattern.values ≈ ones(91, 360) * 2
    @test pattern.max_el == π / 2
end

@testset "Pattern3D" begin
    manifold = IdealManifold(1575420e3, 0.1904 / 4 * SVector(SVector(1, 1, 0), SVector(-1, 1, 0), SVector(1, -1, 0), SVector(-1, -1, 0)))
    pattern = Pattern3D(manifold, norm, num_az = 360, num_el = 181, max_el = π)
    @test pattern.gains ≈ ones(181, 360) * 2
    @test pattern.max_el == 1π
end


text:
module Absolute
    using LoopVectorization

    function func(output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        loss = 0.0
        @avxt for i in eachindex(output_matrix)
            loss += abs(output_matrix[i]-sample_matrix[i])
        end
        return loss
    end
end


text:
include("utils.jl")

height_ = (xs) -> 45f-2Float32(sin(3*xs)) + 55f-2
convert_x = (xs, parameters) -> (xs - -12f-1)*parameters.scale  # min_position = -12f-1
convert_y = (ys, parameters) -> (parameters.screen_height - ys*parameters.scale)

struct MountainCarDrawParams <: AbstractDrawParams
    screen_height::UInt32
    screen_width::UInt32
    scale::Float32
    world_width::Float32
    car_width::Float32
    car_height::Float32
end

function MountainCarDrawParams()

    MountainCarDrawParams(
        400,        # screen_height
        600,        # screen_width
        600/18f-1,  # scale
        18f-1,      # world_width
        40f0,       # car_width
        20f0,       # car_height
    )
end

function Ctx(env::MountainCarEnv, mode::Symbol=:no_render)
    if mode == :human_pane
        draw_parameters = MountainCarDrawParams()
        viewer = CairoRGBSurface(draw_parameters.screen_width, draw_parameters.screen_height)

        CairoCtx(draw_parameters, viewer)
    elseif mode == :human_window
        draw_parameters = MountainCarDrawParams()
        viewer = CairoRGBSurface(draw_parameters.screen_width, draw_parameters.screen_height)

        canvas = @GtkCanvas()
        canvas.backcc = CairoContext(viewer)
        win = GtkWindow(canvas, "MountainCar",
                draw_parameters.screen_width, draw_parameters.screen_height; resizable=false)
        show(canvas)
        visible(win, false)
        signal_connect(win, "delete-event") do widget, event
            ccall((:gtk_widget_hide_on_delete, Gtk.libgtk), Bool, (Ptr{GObject},), win)
        end

        GtkCtx(draw_parameters, canvas, win)
    elseif mode == :rgb
        draw_parameters = MountainCarDrawParams()
        viewer = CairoRGBSurface(draw_parameters.screen_width, draw_parameters.screen_height)

        RGBCtx(draw_parameters, viewer)
    elseif mode == :no_render
        NoCtx()
    else
        error("Unrecognized mode in Ctx(): $(mode)")
    end
end

function drawcanvas!(env::MountainCarEnv, viewer::CairoContext, parameters::MountainCarDrawParams)
    xs = collect(Float32.(range(-12f-1, 6f-1, length=100)))
    ys = convert_y.(height_.(xs), [parameters])
    xs = convert_x.(xs, [parameters])
    set_source_rgb(viewer, 1, 1, 1)
    rectangle(viewer, 0, 0, parameters.screen_width, parameters.screen_height)
    fill(viewer)

    # Track
    set_source_rgb(viewer, 0, 0, 0)
    move_to(viewer, xs[1], ys[1])
    for i=2:length(xs)
        line_to(viewer, xs[i], ys[i])
    end
    stroke(viewer)
    close_path(viewer)

    #Goalpost
    goal_coordinates = Pair(convert_x(env.goal_position, parameters),           # x_cordinate
                            convert_y(height_(env.goal_position), parameters))  # y_coordinate
    set_source_rgb(viewer, 1, 0, 0)
    move_to(viewer, goal_coordinates.first, goal_coordinates.second)
    line_to(viewer, goal_coordinates.first, goal_coordinates.second-50)
    set_line_width(viewer, 8)
    stroke(viewer)
    set_line_width(viewer, 1)  # Resetting line width
    close_path(viewer)

    # Car
    cartx = convert_x(env.state[1].data, parameters)
    carty = convert_y(height_(env.state[1].data), parameters)
    set_source_rgb(viewer, 0, 0, 0)
    translate(viewer, cartx, carty)
    #rotate(viewer, cos(3state[1]))
    rectangle(viewer, parameters.car_width/2, 0, -parameters.car_width, -parameters.car_height)  # cartx-car_width/2
    fill(viewer)

    # Wheels
    set_source_rgb(viewer, 0.5, 0.5, 0.5)
    circle(viewer, parameters.car_width/4, -5, 5)
    circle(viewer, -parameters.car_width/4, -5, 5)
    fill(viewer)
    translate(viewer, -cartx, -carty)  # Resetting the cursor

    #=
    # Calibration circle
    set_source_rgb(viewer, 0, 0, 0)
    calibration_radius = 10
    circle(viewer, calibration_radius, calibration_radius, calibration_radius)
    set_line_width(viewer, 1)
    stroke(viewer)
    =#
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

bat_default(
    ::typeof(bat_eff_sample_size),
    ::Val{:algorithm},
    samples::AbstractVectorOfSimilarVectors{<:Real}
) = EffSampleSizeFromAC()

bat_default(
    ::typeof(bat_eff_sample_size),
    ::Val{:algorithm},
    samples::DensitySampleVector
) = EffSampleSizeFromAC()


text:
module BenchForeach

import ThreadsX
using BenchmarkTools
using Referenceables: referenceable

@inline function f(a, b, c)
    a[] = b + c
end

n = 2000
B = randn(n, n)
C = randn(n, n)
A = similar(B)
suite = BenchmarkGroup()

for (label, _foreach) in [("base", foreach), ("tx", ThreadsX.foreach)]
    s = suite[label] = BenchmarkGroup()
    s["A .= B .+ C"] = @benchmarkable($_foreach(f, referenceable($A), $B, $C))
    s["A .= B .+ B'"] = @benchmarkable($_foreach(f, referenceable($A), $B, $B'))
end

let s = suite["broadcast"] = BenchmarkGroup()
    s["A .= B .+ C"] = @benchmarkable $A .= $B .+ $C
    s["A .= B .+ B'"] = @benchmarkable $A .= $B .+ $B'
end

end  # module
BenchForeach.suite


text:
using CSV
using DataFrames

function bom(filename::String, c::Circuit, query::Bool=false)
    df = DataFrame(ref=String[], manf_num=Union{String,Missing}[])
    for e in c.elements
        push!(df, (string(e), get(e.meta, :manf_num, missing)))
    end
    count_parts = dropmissing(by(df, :manf_num, :ref => length),
                              disallowmissing=true)
    rename!(count_parts, :manf_num => Symbol("manf#"), :ref_length => :qty)
    CSV.write(filename, count_parts)
    if query
        xl = rsplit(filename, ".", limit=2)[1]*".xlsx"
        isfile(xl) && rm(xl)
        run(`kicost -i $(filename) --eda csv`)
        xl
    else
        filename
    end
end

bom(c::Circuit, query::Bool=false) = bom("bom.csv", c, query)


text:
import IntertemporalChoiceHeuristics:
    load,
    update_weights!,
    fit,
    models,
    canonical_parameters

import Optim: minimizer

using ProgressMeter
using ArgParse

function main(demo::Bool, output_path::String)
    io = open(output_path, "w")

    inputs, weights = load()

    new_weights = copy(weights)

    if demo
        n_replicates = 100
    else
        n_replicates = 10_000
    end

    @showprogress 1 "Bootstrapping:" for replicate in 1:n_replicates
        update_weights!(new_weights, weights, 1.0)

        for model in models
            try
                res = fit(model, inputs, new_weights)
                parameters = minimizer(res)
                m = model(parameters)
                params = canonical_parameters(m)
                K = length(params)
                @printf(
                    io,
                    "%d\t%s\t%s\n",
                    replicate,
                    model,
                    join([i > K ? "NA" : string(params[i]) for i in 1:5], '\t')
                )
            end
        end
    end

    close(io)
end

s = ArgParseSettings(description = "generate_replicates.jl")

@add_arg_table s begin
    "--demo"
        action = :store_true
    "output_path"
        required = true
end

args = parse_args(s)

main(args["demo"], args["output_path"])


text:
var_idx(x::JuMP.VariableRef) = JuMP.optimizer_index(x).value
var_idx(x::MOI.VariableIndex) = x.value
# support for @variable(m, x, Set)
function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    set::T,
) where {T<:MOI.AbstractScalarSet}
    return JuMP.VariableConstrainedOnCreation(JuMP.ScalarVariable(info), set)
end

"""
Optimizer struct
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::CS.ConstraintSolverModel
    variable_info::Vector{Variable}
    # which variable index, (:leq,:geq,:eq,:Int,:Bin), and lower and upper bound
    var_constraints::Vector{Tuple{Int,Symbol,Int,Int}}
    status::MOI.TerminationStatusCode
    options::SolverOptions
end

include("util.jl")
include("variables.jl")
include("Bridges/util.jl")
include("Bridges/complement.jl")
include("Bridges/indicator.jl")
include("Bridges/reified.jl")
include("Bridges/strictly_greater_than.jl")
include("Bridges/bool.jl")
include("constraints.jl")
include("objective.jl")
include("results.jl")

MOI.get(::Optimizer, ::MOI.SolverName) = "ConstraintSolver"

"""
Optimizer struct constructor
"""
function Optimizer(; options...)
    options = combine_options(options)
    com = CS.ConstraintSolverModel(options.solution_type)
    optimizer = Optimizer(com, [], [], MOI.OPTIMIZE_NOT_CALLED, options)
    lbo = MOIB.full_bridge_optimizer(optimizer, options.solution_type)
    greater2less_bridges = [
        MOIBC.GreaterToLessBridge{options.solution_type},
        CS.StrictlyGreaterToStrictlyLessBridge{options.solution_type}
    ]
    inner_bridges = greater2less_bridges
    for i in 1:length(inner_bridges)
        push!(inner_bridges, CS.ComplementBridge{options.solution_type, inner_bridges[i]})
    end
    # the bool bridge supports all inner bridges
    for i in 1:length(inner_bridges)
        push!(inner_bridges, CS.BoolBridge{options.solution_type, inner_bridges[i]})
    end
    for i in 1:length(inner_bridges)
        push!(inner_bridges, CS.ComplementBridge{options.solution_type, inner_bridges[i]})
    end

    for inner_bridge in inner_bridges
        MOIB.add_bridge(lbo, inner_bridge)
        MOIB.add_bridge(lbo, CS.IndicatorBridge{options.solution_type, inner_bridge})
        MOIB.add_bridge(lbo, CS.ReifiedBridge{options.solution_type, inner_bridge})
    end
    return lbo
end

"""
    MOI.is_empty(model::Optimizer)
"""
function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) && isempty(model.var_constraints)
end

"""
    MOI.empty!(model::Optimizer)
"""
function MOI.empty!(model::Optimizer)
    model.inner = CS.ConstraintSolverModel(model.options.solution_type)
    empty!(model.variable_info)
    empty!(model.var_constraints)
    model.status = MOI.OPTIMIZE_NOT_CALLED
    # !important => don't remove the options
end

"""
Copy constructor for the optimizer
"""
MOIU.supports_default_copy_to(model::Optimizer, copy_names::Bool) = !copy_names
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(model, src; kws...)
end

MOI.supports(::Optimizer, ::MOI.RawParameter) = true
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

"""
    MOI.set(model::Optimizer, p::MOI.RawParameter, value)

Set a RawParameter to `value`
"""
function MOI.set(model::Optimizer, p::MOI.RawParameter, value)
    current_options_type = SolverOptions
    current_options_obj = model.options
    p_symbol = Symbol(p.name)

    # being able to parse simple options like "traverse_strategy"
    # as well as sub category options like "activity.decay"
    num_subcat = 0
    parts = split(p.name, ".")
    for pname in parts
        num_subcat += 1
        cp_symbol = Symbol(pname)
        if in(cp_symbol, fieldnames(current_options_type))
            type_of_param = fieldtype(current_options_type, cp_symbol)
            if num_subcat == length(parts)
                if hasmethod(convert, (Type{type_of_param}, typeof(value)))
                    if is_possible_option_value(p, value)
                        setfield!(
                            current_options_obj,
                            cp_symbol,
                            convert(type_of_param, value),
                        )
                    else
                        @error "The option $(cp_symbol) doesn't have $(value) as a possible value. Possible values are: $(POSSIBLE_OPTIONS[p_symbol])"
                        break
                    end
                else
                    @error "The option $(p.name) has a different type ($(type_of_param))"
                    break
                end
            else
                current_options_type = type_of_param
                current_options_obj = getfield(current_options_obj, cp_symbol)
            end
        else
            @error "The option $(p.name) doesn't exist."
            break
        end
    end
    return
end

"""
    MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)

Set the time limit
"""
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Union{Nothing,Float64})
    if value === nothing
        model.options.time_limit = Inf
    else
        model.options.time_limit = value
    end
    return
end

"""
    MOI.optimize!(model::Optimizer)
"""
function MOI.optimize!(model::Optimizer)
    model.inner.options = model.options
    # check if every variable has bounds and is an Integer
    check_var_bounds(model)

    set_pvals!(model)
    set_var_in_all_different!(model)

    create_lp_model!(model)

    status = solve!(model.inner)
    set_status!(model, status)

    if status == :Solved
        com = model.inner
        com.solutions = unique!(sol -> sol.hash, com.solutions)
        if !com.options.all_solutions
            filter!(sol -> sol.incumbent == com.best_sol, com.solutions)
        end
        sort_solutions!(com)
    end
end


text:
F84(κ::F,
    πA::F, πC::F, πG::F, πT::F,
    safe::Bool=true) where F <: Float64 =
  F84rel(κ, πA, πC, πG, πT, safe)

F84(κ::F, β::F,
    πA::F, πC::F, πG::F, πT::F,
    safe::Bool=true) where F <: Float64 =
  F84abs(κ, β, πA, πC, πG, πT, safe)

function F84(θ_vec::A,
             π_vec::A,
             safe::Bool=true) where A <: AbstractArray
  if safe && length(π_vec) != 4
    error("Incorrect base frequency vector length")
  end
  if length(θ_vec) == 1
    return F84rel(θ_vec[1],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe)
  elseif length(θ_vec) == 2
    return F84abs(θ_vec[1], θ_vec[2],
                  π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T],
                  safe)
  else
    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
  end
end

function F84rel(θ_vec::A,
                π_vec::A,
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 1
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return F84rel(θ_vec[1], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
end

function F84abs(θ_vec::A,
                π_vec::A,
                safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 2
      error("Incorrect parameter vector length")
    elseif length(π_vec) != 4
      error("Incorrect base frequency vector length")
    end
  end
  return F84abs(θ_vec[1], θ_vec[2], π_vec[DNA_A], π_vec[DNA_C], π_vec[DNA_G], π_vec[DNA_T], safe)
end