text:
include("models/data_month.jl")
include("models/data_kw.jl")

Ntotal = @combine(groupby(mordat,[:jahr,:monat]), :Ntotal=sum(:N))
Ntotalkw = @combine(groupby(mordatkw,[:jahr,:kw]), :Ntotal=sum(:N))
@transform! groupby(mordat, :date) :P_AG = :N ./ sum(:N)
@transform! groupby(mordatkw, :date) :P_AG = :N ./ sum(:N)

function referenceP_AG(mordat)
    @combine(groupby(mordat,:group),
             :referenceP_AG = mean(:P_AG)
             )
end

function referenceP_AG(mordat,date)
    @select(filter(x->x.date==date, mordat),
            :group,
            :referenceP_AG = :P_AG
            )
end

# (@transform referenceP_AG(md, Date(2000,1,1)) :marginalN = 80e6 .* :referenceP_AG)[:,end] |> sum
function adjustierte_sterbefälle(md, rP, time = [:jahr,:monat], Ntotal = @combine(groupby(md,time), :Ntotal=sum(:N)))
    @info "data" md[time..., :group, :P_AG, :PD, :D]
    mmd = leftjoin(md[:,[time..., :group, :P_AG, :PD, :D]], rP, on=:group)
    mmd = leftjoin(mmd, Ntotal, on=time)
    @transform! mmd :PDreferenceP_AG = :PD .* :referenceP_AG
    @combine(groupby(mmd, time),
             :adjD = sum(:PD .* :referenceP_AG) * :Ntotal[1],
             :D = sum(:D),
             :adjPD = sum(:PD .* :referenceP_AG))
    #)
end

#@transform! groupby(mordat,:date) :stackedN = cumsum(:N)

ag_col=let na = length(DeutschlandDaten.Altersgruppen_Monate)
    m = cgrad([:lightblue, :gray],na+5; scale=:log)
    w = cgrad([:lightpink, :gray],na+5; scale=:log)
    vcat([ [ m[x], w[x] ] for x in 1:na ]...)'
end

let title="Demographische Anteile (kumulative %)\nder Bundesrepublik Deutschland der Jahre",
    file = "Demographie_Deutschland_Prozent",
    p = @df mordat yearsareaplot(:date,  100 * :P_AG, :group;
                                 title=title,
                                 seriescolor=ag_col,
                                 ylim=(00,100), yticks=0:10:100, linecolor=:black, rightmargin=40px, bottommargin=20px)
    plot!(twinx(),yticks = 0:10:100,
          ylim=(00,100), xlim=(Date(2000,1,1),Date(2021,12,1)),
          ymirror=true, xticks=nothing,
          rightmargin=100px
          )
    DeutschlandDaten.savefigs(p, file, data = "Datenquelle: Statistisches Bundesamt Code: 12411-0006")
end

let title="Demographie (Bevölkerungszahl)\nder Bundesrepublik Deutschland der Jahre",
    file = "Demographie_Deutschland_interpoliert",
    yt = 0:8,
    yl=(00,85e6),
    plotattrs = (yticks = (1e7*yt, convert.(Int,10*yt)), ylim=yl, rightmargin=40px, bottommargin=20px)
    p = @df filter(x->x.date <= Date(2020,1,1), mordat) yearsareaplot(:date, :N, :group;
                                                                      title=title,
                                                                      seriescolor=ag_col,
                                                                      ylabel="Millionen",
                                                                      linecolor=:black,
                                                                      plotattrs...)
    plot!(twinx(); xlim=(Date(2000,1,1),Date(2021,1,1)),
          ymirror=true, xticks=nothing,
          plotattrs...
              )
    DeutschlandDaten.savefigs(p, file, data = "Datenquelle: Statistisches Bundesamt Code: 12411-0006")
end

let title="Sterbezahlen\n Bundesrepublik Deutschland",
    file = "Sterbezahlen_Deutschland_monatlich"
    p = @df mordat yearsareaplot(:date, :D, :group; legend=nothing, seriescolor=ag_col, linecolor=:black, yticks=(0000:10000:100000, 0:10:100),
                             title=title,
                             ylabel="Tausend / Monat",
                             annotation_margin =-500, xlim=(Date(1998,5,1),Date(2023,7,1)))
    DeutschlandDaten.savefigs(p, file)
    p
end

if false
    @df @combine(groupby(mordat,:date), :D = sum(:D)) plot(:date, :D, seriestype=:bar, legend=nothing, fillcolor=:black, linecolor=:darkgray, linewidth=.001)
end

let title = "Verstorbenenzahl, Demographie-adjustiert\nBundesrepublik Deutschland",
    file = "adjustierte_Sterbezahlen_Deutschland_monatlich"
    p= @df @combine(groupby(mordat,:date), :D = sum(:D)) plot(; legend=nothing, fillcolor=:black, linecolor=:black, linewidth=2, ylim=(0,130000/1000),
                                                              title =title,
                                                              ylabel="Tausend / Monat",
                                                              xticks=([Date(y,1,1) for y in 2000:2021], 2000:2021),
                                                              xrotation=90, xmirror=true,
                                                              topmargin=30px)
    for j in (2021,2000)
        rP = referenceP_AG(filter(x->x.jahr == j,mordat))
        sp = filter(x->x.jahr in 2000:2021, adjustierte_sterbefälle(mordat, rP, Ntotal))
        faktor = mean(filter(x->x.jahr==j, Ntotal).Ntotal)#lastD/sp.adjPD[end]
        @df sp plot!(p, Date.(:jahr,:monat,1), :adjPD .* faktor/1000, label="auf Demographie $j adjustiert", width=1)
    end
    @df @combine(groupby(mordat,:date), :D = sum(:D)) plot!(p,:date, :D/1000, legend=:bottomright, label="offizielle Sterbefallzahl", fillcolor=:black, linecolor=:black, linewidth=2)
    DeutschlandDaten.savefigs(p, file)
    p
end
@df @combine(groupby(mordat,:date), :D = sum(:D)) plot!(p,:date, :D/1000, legend=:bottomright, label="offizielle Sterbefallzahl", fillcolor=:black, linecolor=:black, linewidth=2)
p
DeutschlandDaten.savefigs(p, "adjustierte_Sterbezahlen_Deutschland_monatlich")

if false
j=2021
d = adjustierte_sterbefälle(mordat, referenceP_AG(filter(x->x.jahr == j,mordat)), Ntotal)
@transform! d :adjD_24m = gliding_mean(:adjD,12,0)

@df d plot(Date.(:jahr,:monat,1), :adjD )

@df d plot(Date.(:jahr,:monat,1), :adjD ./ :adjD_24m )

@df d plot(:monat, :adjD; group=:jahr, linecolor=[ getindex(cgrad(:reds,22),j-2000+1) for j in :jahr], legend=nothing )

p = @df filter(x->x.jahr in 2010:2019,d) boxplot(:monat, :adjD, :jahr, alpha=.25, label="2010-2019",
                                             title="Adjustierte Sterbezahlen\nBundesrepublik Deutschland",
                                             topmargin=10px,
                                             ylabel="Tausend",
                                             xticks=(1:12,DeutschlandDaten.monattick),
                                             ylim=10000 .* (7,13),
                                             yticks=let tsd=7:13; (10000*tsd, 10*tsd); end)
@df filter(x->x.jahr in 2000:2009,d) boxplot!(:monat, :adjD, :jahr, alpha=.1, label="2000-2009",
                                              xticks=(1:12,DeutschlandDaten.monattick),
                                              ylim=10000 .* (7,13),
                                              yticks=let tsd=7:13; (10000*tsd, 10*tsd); end)
@df filter(x->x.jahr in 2020:2021,d) plot!(:monat, :adjD; group=:jahr, legend=:topright, width=2)
    # @df filter(x->x.jahr in 2020:2021,d) plot!(:monat, :D; group=:jahr, legend=:topright, width=2)
DeutschlandDaten.savefigs(p, "adjustierte_Sterbezahlen_2020_2021_monatlich")

@df filter(x->x.jahr in 2010:2019,mordat) plot(; title="Adjustierte Sterbezahlen\nBundesrepublik Deutschland",
                                               topmargin=10px,
                                               ylim=(-10,1),
                                               xlim=(1,12),
                                               xticks=(1:12,DeutschlandDaten.monattick))

for a in DeutschlandDaten.Altersgruppen_Monate
    d = filter(x->x.alter==a && x.geschlecht=="Männlich",mordat)
    p = @df filter(x -> x.jahr in 2010:2019,d) boxplot(:monat, :PD; group=:group, alpha=.25, label="2010-2019", bar_width=[1,1],
                                                   title="Sterbewahrscheinlichkeit $(alterstring(a))",
                                                   xlim=(0,13),
                                                   xticks=(1:12,DeutschlandDaten.monattick))
    @df filter(x -> x.jahr in 2000:2009,d) boxplot!(:monat, :PD; group=:group, alpha=.25, label="2000-2009", bar_width=[1,1],
                                                    xlim=(0,13),
                                                    xticks=(1:12,DeutschlandDaten.monattick))
    @df filter(x -> x.jahr in 2020:2021,d) plot!(:monat, :PD; group=:jahr, width=1)
    DeutschlandDaten.savefigs(p, "Sterbewahrscheinlichkeit_$a")
end

@df filter(x->x.geschlecht=="Männlich" && x.jahr in 2000:2019,mordat) groupedboxplot(:monat, :logitPD; group=:group, alpha=.25, label="2010-2019", bar_width=[1,1],
                                               xlim=(0,13),
                                               xticks=(1:12,DeutschlandDaten.monattick))

plot([
begin
    p = @df filter(x->x.geschlecht==g && x.jahr in 2000:2019,mordat) boxplot(:monat, :logitPD; group=:group, alpha=.25, label="2010-2019", bar_width=[1,1],
                                                                             xrotation=90,
                                                                             color=:gray, linecolor=:darkgray,
                                                                             markersize=.5,
                                                                             xlim=(0,13),
                                                                             xticks=(1:12,DeutschlandDaten.monattick))
    # p = @df filter(x->x.geschlecht==g && x.jahr in 2000:2009,mordat) boxplot!(:monat, :logitPD; group=:group, alpha=.15, label="2010-2019", bar_width=[1,1],
    #                                                                          xrotation=90,
    #                                                                          color=:green, linecolor=:darkgreen,
    #                                                                          markersize=.5,
    #                                                                          xlim=(0,13),
    #                                                                          xticks=(1:12,DeutschlandDaten.monattick))
    @df filter(x->x.geschlecht==g && x.jahr in 2020:2021,mordat) plot!(:monat, :logitPD; group=(:group,:jahr), color=:jahr, legend=nothing, width=2)
    p
end
    for g in ["Männlich", "Weiblich"]
]...)

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
function add_re!(sys)
    re = RenewableDispatch(
        "WindBusA",
        true,
        get_component(Bus, sys, "bus5"),
        0.0,
        0.0,
        1.200,
        PrimeMovers.WT,
        (min = 0.0, max = 0.0),
        1.0,
        TwoPartCost(0.220, 0.0),
        100.0,
    )
    add_component!(sys, re)
    copy_time_series!(re, get_component(PowerLoad, sys, "bus2"))

    fx = RenewableFix(
        "RoofTopSolar",
        true,
        get_component(Bus, sys, "bus5"),
        0.0,
        0.0,
        1.100,
        PrimeMovers.PVe,
        1.0,
        10.0,
    )
    add_component!(sys, fx)
    copy_time_series!(fx, get_component(PowerLoad, sys, "bus2"))

    for g in get_components(HydroEnergyReservoir, sys)
        tpc = get_operation_cost(g)
        smc = StorageManagementCost(
            variable = get_variable(tpc),
            fixed = get_fixed(tpc),
            start_up = 0.0,
            shut_down = 0.0,
            energy_shortage_cost = 10.0,
            energy_surplus_cost = 10.0,
        )
        set_operation_cost!(g, smc)
    end

    batt = GenericBattery(
        "test_batt",
        true,
        get_component(Bus, sys, "bus4"),
        PrimeMovers.BA,
        0.0,
        (min = 0.0, max = 1.0),
        1.0,
        0.0,
        (min = 0.0, max = 1.0),
        (min = 0.0, max = 1.0),
        (in = 1.0, out = 1.0),
        0.0,
        nothing,
        10.0,
    )
    add_component!(sys, batt)
end

function run_test_sim(result_dir::String)
    sim_name = "results_sim"
    sim_path = joinpath(result_dir, sim_name)

    if ispath(sim_path)
        @info "Reading UC system from" sim_path
        c_sys5_hy_uc = System(
            joinpath(sim_path, "..", "c_sys5_hy_uc.json"),
            time_series_read_only = true,
        )
        @info "Reading ED system from" sim_path
        c_sys5_hy_ed = System(
            joinpath(sim_path, "..", "c_sys5_hy_ed.json"),
            time_series_read_only = true,
        )
        executions = tryparse.(Int, readdir(sim_path))
        sim = joinpath(sim_path, string(maximum(executions)))
        @info "Reading results from last execution" sim
    else
        @info "Building UC system from"
        c_sys5_hy_uc = PSB.build_system(PSB.SIIPExampleSystems, "5_bus_hydro_uc_sys")
        @info "Building ED system from"
        c_sys5_hy_ed = PSB.build_system(PSB.SIIPExampleSystems, "5_bus_hydro_ed_sys")

        @info "Adding extra RE"
        add_re!(c_sys5_hy_uc)
        add_re!(c_sys5_hy_ed)
        to_json(c_sys5_hy_uc, joinpath(sim_path, "..", "c_sys5_hy_uc.json"), force = true)
        to_json(c_sys5_hy_ed, joinpath(sim_path, "..", "c_sys5_hy_ed.json"), force = true)

        mkpath(result_dir)
        GLPK_optimizer =
            optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)

        template_hydro_st_uc = template_unit_commitment()
        set_device_model!(template_hydro_st_uc, HydroDispatch, FixedOutput)
        set_device_model!(
            template_hydro_st_uc,
            HydroEnergyReservoir,
            HydroDispatchReservoirStorage,
        )
        set_device_model!(template_hydro_st_uc, GenericBattery, BookKeeping)

        template_hydro_st_ed = template_economic_dispatch(
            network = CopperPlatePowerModel,
            use_slacks = true,
            duals = [CopperPlateBalanceConstraint],
        )
        set_device_model!(template_hydro_st_ed, HydroDispatch, FixedOutput)
        set_device_model!(
            template_hydro_st_ed,
            HydroEnergyReservoir,
            HydroDispatchReservoirStorage,
        )
        set_device_model!(template_hydro_st_ed, GenericBattery, BookKeeping)
        template_hydro_st_ed.services = Dict() #remove ed services
        models = SimulationModels(
            decision_models = [
                DecisionModel(
                    template_hydro_st_uc,
                    c_sys5_hy_uc,
                    optimizer = GLPK_optimizer,
                    name = "UC",
                    system_to_file = false,
                ),
                DecisionModel(
                    template_hydro_st_ed,
                    c_sys5_hy_ed,
                    optimizer = GLPK_optimizer,
                    name = "ED",
                    system_to_file = false,
                ),
            ],
        )

        sequence = SimulationSequence(
            models = models,
            feedforwards = feedforward = Dict(
                "ED" => [
                    SemiContinuousFeedforward(
                        component_type = ThermalStandard,
                        source = OnVariable,
                        affected_values = [ActivePowerVariable],
                    ),
                    EnergyLimitFeedforward(
                        component_type = HydroEnergyReservoir,
                        source = ActivePowerVariable,
                        affected_values = [ActivePowerVariable],
                        number_of_periods = 12,
                    ),
                ],
            ),
            ini_cond_chronology = InterProblemChronology(),
        )
        sim = Simulation(
            name = "results_sim",
            steps = 2,
            models = models,
            sequence = sequence,
            simulation_folder = result_dir,
        )
        build!(sim)
        execute!(sim)
    end

    results = SimulationResults(sim)
    results_uc = get_decision_problem_results(results, "UC")
    set_system!(results_uc, c_sys5_hy_uc)
    results_ed = get_decision_problem_results(results, "ED")
    set_system!(results_ed, c_sys5_hy_ed)

    return results_uc, results_ed
end

function run_test_prob()
    c_sys5_hy_uc = PSB.build_system(PSB.SIIPExampleSystems, "5_bus_hydro_uc_sys")
    add_re!(c_sys5_hy_uc)
    GLPK_optimizer =
        optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => GLPK.GLP_MSG_OFF)

    template_hydro_st_uc = template_unit_commitment()
    set_device_model!(template_hydro_st_uc, HydroDispatch, FixedOutput)
    set_device_model!(
        template_hydro_st_uc,
        HydroEnergyReservoir,
        HydroDispatchReservoirStorage,
    )

    prob = DecisionModel(
        template_hydro_st_uc,
        c_sys5_hy_uc,
        optimizer = GLPK_optimizer,
        horizon = 12,
    )
    build!(prob, output_dir = mktempdir())
    solve!(prob)
    res = ProblemResults(prob)
    return res
end

ext: jl

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
module EquationModule

import ..ProgramConstantsModule: CONST_TYPE
import ..OptionsStructModule: Options

# Define a serialization format for the symbolic equations:
mutable struct Node
    #Holds operators, variables, constants in a tree
    degree::Int #0 for constant/variable, 1 for cos/sin, 2 for +/* etc.
    constant::Bool #false if variable
    val::CONST_TYPE
    # ------------------- (possibly undefined below)
    feature::Int #Either const value, or enumerates variable.
    op::Int #enumerates operator (separately for degree=1,2)
    l::Node
    r::Node

    Node(d::Int, c::Bool, v::CONST_TYPE) = new(d, c, v)
    Node(d::Int, c::Bool, v::CONST_TYPE, f::Int) = new(d, c, v, f)
    Node(d::Int, c::Bool, v::CONST_TYPE, f::Int, o::Int, l::Node) = new(d, c, v, f, o, l)
    function Node(d::Int, c::Bool, v::CONST_TYPE, f::Int, o::Int, l::Node, r::Node)
        return new(d, c, v, f, o, l, r)
    end
end

Node(val::CONST_TYPE) = Node(0, true, val) #Leave other values undefined
"""
    Node(feature::Int)

Create a variable node using feature `feature::Int`
"""
Node(feature::Int) = Node(0, false, convert(CONST_TYPE, 0.0f0), feature)
"""
    Node(op::Int, l::Node)

Apply unary operator `op` (enumerating over the order given) to `Node` `l`
"""
Node(op::Int, l::Node) = Node(1, false, convert(CONST_TYPE, 0.0f0), 0, op, l)
"""
    Node(op::Int, l::Union{AbstractFloat, Int})

Short-form for creating a scalar/variable node, and applying a unary operator
"""
function Node(op::Int, l::Union{AbstractFloat,Int})
    return Node(1, false, convert(CONST_TYPE, 0.0f0), 0, op, Node(l))
end
"""
    Node(op::Int, l::Node, r::Node)

Apply binary operator `op` (enumerating over the order given) to `Node`s `l` and `r`
"""
Node(op::Int, l::Node, r::Node) = Node(2, false, convert(CONST_TYPE, 0.0f0), 0, op, l, r)
"""
    Node(op::Int, l::Union{AbstractFloat, Int}, r::Node)

Short-form to create a scalar/variable node, and apply a binary operator
"""
function Node(op::Int, l::Union{AbstractFloat,Int}, r::Node)
    return Node(2, false, convert(CONST_TYPE, 0.0f0), 0, op, Node(l), r)
end
"""
    Node(op::Int, l::Node, r::Union{AbstractFloat, Int})

Short-form to create a scalar/variable node, and apply a binary operator
"""
function Node(op::Int, l::Node, r::Union{AbstractFloat,Int})
    return Node(2, false, convert(CONST_TYPE, 0.0f0), 0, op, l, Node(r))
end
"""
    Node(op::Int, l::Union{AbstractFloat, Int}, r::Union{AbstractFloat, Int})

Short-form for creating two scalar/variable node, and applying a binary operator
"""
function Node(op::Int, l::Union{AbstractFloat,Int}, r::Union{AbstractFloat,Int})
    return Node(2, false, convert(CONST_TYPE, 0.0f0), 0, op, Node(l), Node(r))
end
"""
    Node(val::AbstractFloat)

Create a scalar constant node
"""
Node(val::AbstractFloat) = Node(convert(CONST_TYPE, val))
"""
    Node(var_string::String)

Create a variable node, using the format `"x1"` to mean feature 1
"""
Node(var_string::String) = Node(parse(Int, var_string[2:end]))
"""
    Node(var_string::String, varMap::Array{String, 1})

Create a variable node, using a user-passed format
"""
function Node(var_string::String, varMap::Array{String,1})
    return Node(
        [i for (i, _variable) in enumerate(varMap) if _variable == var_string][1]::Int
    )
end

# Copy an equation (faster than deepcopy)
function copy_node(tree::Node)::Node
    if tree.degree == 0
        if tree.constant
            return Node(copy(tree.val))
        else
            return Node(copy(tree.feature))
        end
    elseif tree.degree == 1
        return Node(copy(tree.op), copy_node(tree.l))
    else
        return Node(copy(tree.op), copy_node(tree.l), copy_node(tree.r))
    end
end

function string_op(
    op::F,
    tree::Node,
    options::Options;
    bracketed::Bool=false,
    varMap::Union{Array{String,1},Nothing}=nothing,
)::String where {F}
    if op in [+, -, *, /, ^]
        l = string_tree(tree.l, options; bracketed=false, varMap=varMap)
        r = string_tree(tree.r, options; bracketed=false, varMap=varMap)
        if bracketed
            return "$l $(string(op)) $r"
        else
            return "($l $(string(op)) $r)"
        end
    else
        l = string_tree(tree.l, options; bracketed=true, varMap=varMap)
        r = string_tree(tree.r, options; bracketed=true, varMap=varMap)
        return "$(string(op))($l, $r)"
    end
end

"""
    string_tree(tree::Node, options::Options; kws...)

Convert an equation to a string.

# Arguments

- `varMap::Union{Array{String, 1}, Nothing}=nothing`: what variables
    to print for each feature.
"""
function string_tree(
    tree::Node,
    options::Options;
    bracketed::Bool=false,
    varMap::Union{Array{String,1},Nothing}=nothing,
)::String
    if tree.degree == 0
        if tree.constant
            return string(tree.val)
        else
            if varMap === nothing
                return "x$(tree.feature)"
            else
                return varMap[tree.feature]
            end
        end
    elseif tree.degree == 1
        return "$(options.unaops[tree.op])($(string_tree(tree.l, options, bracketed=true, varMap=varMap)))"
    else
        return string_op(
            options.binops[tree.op], tree, options; bracketed=bracketed, varMap=varMap
        )
    end
end

# Print an equation
function print_tree(
    io::IO, tree::Node, options::Options; varMap::Union{Array{String,1},Nothing}=nothing
)
    return println(io, string_tree(tree, options; varMap=varMap))
end

function print_tree(
    tree::Node, options::Options; varMap::Union{Array{String,1},Nothing}=nothing
)
    return println(string_tree(tree, options; varMap=varMap))
end

end


text:
struct Linear
    W
    b
end

foo(x) = x + 2

function linear_cost(m::Linear, x, y)
    sum((m.W * x .+ m.b .- y) .^ 2)
end

function linear_cost_flat(W, b, x, y)
    sum((W * x .+ b .- y) .^ 2)
end

let
    W = randn(500, 784)
    b = randn(500)
    m = Linear(W, b)
    x = randn(784, 100)
    y = randn(500)
    inputs = [:m => m, :x => x, :y => y]
    flat_inputs = [:W => W, :b => b, :x => x, :y => y]

    # first check that df with flat model is correct
    test_compare(linear_cost_flat; flat_inputs...)

    # now check that results from flat and structured models are the same
    df_flat = xdiff(linear_cost_flat; flat_inputs...)
    results_flat = df_flat(W, b, x, y)
    df = xdiff(linear_cost; inputs...)
    results_flat = df(W, b, x, y)
    results_struct = df(m, x, y)
    @test results_flat == results_struct
end


text:
# promotion.jl
#

# Promotion
promote(a::TaylorModelN{N,Interval{T},S}, b::S) where {N, T<:Real, S<:Real} =
    (a, TaylorModelN(IntervalArithmetic.atomic(Interval{T}, b), get_order(a), a.x0, a.I))
promote(b::S, a::TaylorModelN{N,Interval{T},S}) where {N, T<:Real, S<:Real} =
    reverse( promote(a, b) )
#
promote(a::TaylorModelN{N,Interval{T},S}, b::R) where {N, T<:Real, S<:Real, R<:Real} =
    promote(a, IntervalArithmetic.atomic(Interval{T}, b))
promote(b::R, a::TaylorModelN{N,Interval{T},S}) where {N, T<:Real, S<:Real, R<:Real} =
    reverse( promote(a, b) )
#
promote(a::TaylorModelN{N,Interval{T},S}, b::Interval{T}) where {N, T<:Real, S<:Real} =
    (a, TaylorModelN(b, get_order(a), a.x0, a.I))
promote(b::Interval{T}, a::TaylorModelN{N,Interval{T},S}) where {N, T<:Real, S<:Real} =
    reverse( promote(a, b) )
#
function promote(a::TaylorModelN{N,T,S}, b::TaylorN{R}) where {N, T, S, R}
    RR = promote_type(T,R)
    aa = TaylorModelN(convert(TaylorN{RR},a.pol), a.rem, a.x0, a.I)
    bb = TaylorModelN(convert(TaylorN{RR},b), 0..0, a.x0, a.I)
    (aa, bb)
end
promote(b::TaylorN{R}, a::TaylorModelN{N,T,S}) where {N, T, S, R} = reverse( promote(a, b) )

function promote(a::Taylor1{TaylorModelN{N,Interval{T},S}}, b::Taylor1{Interval{T}}) where
        {N, T<:Real, S<:Real}

    orderTMN = get_order(a[0])
    bTN = Array{TaylorModelN{N,Interval{T},S}}(get_order(a)+1)
    for ord = 1:length(bTN)
        bTN[ord] = b.coeffs[ord]*TaylorModelN(1..1, orderTMN, a[0].x0, a[0].I)
    end
    (a, Taylor1(bTN))
end
promote(b::Taylor1{Interval{T}}, a::Taylor1{TaylorModelN{N,Interval{T},S}}) where
    {N, T<:Real, S<:Real} = reverse( promote(a, b) )

# function promote(a::Taylor1{TaylorModelN{N,T,S}}, b::Taylor1{S}) where {N,T<:Real,S<:Real}
#     (a, Taylor1(TaylorModelN(interval(b), get_order(a), a.x0, a.I), b.order))
# end
# promote(b::Taylor1{Interval{T}}, a::Taylor1{TaylorModelN{N,Interval{T},S}}) where
#     {N, T<:Real, S<:Real} = promote(a, b)

for TM in tupleTMs
    @eval promote(a::$TM{T,S}, b::T) where {T,S} =
        (a, $(TM)(Taylor1([b], get_order(a)), zero(a.rem), a.x0, a.I))
    @eval promote(b::T, a::$TM{T,S}) where {T,S} = reverse( promote(a,b) )
    #
    @eval promote(a::$TM{T,S}, b::S) where {T,S} =
        (a, $TM(Taylor1([convert(T, b)], get_order(a)), zero(a.rem), a.x0, a.I))
    @eval promote(b::S, a::$TM{T,S}) where {T,S} = reverse( promote(a,b) )
    #
    @eval promote(a::$TM{T,S}, b::R) where {T,S,R} = promote(a, convert(S, b))
    @eval promote(b::R, a::$TM{T,S}) where {T,S,R} = reverse( promote(a,b) )
    #
    @eval function promote(a::$TM{TaylorModelN{N,T,S},S}, b::T) where {N,T,S}
        tmN = TaylorModelN(b, get_order(a.pol[0]), a.pol[0].x0, a.pol[0].I)
        (a, $TM(Taylor1([tmN], get_order(a)), zero(a.rem), a.x0, a.I))
    end
    @eval promote(b::T, a::$TM{TaylorModelN{N,T,S},S}) where {N,T,S} = reverse( promote(a,b) )
end


text:
"""
    struct LinearCombinationDofVector{T} <: AbstractVector{Dof}
      change_of_basis::Matrix{T}
      dof_basis::AbstractVector{<:Dof}
    end

Type that implements a dof basis (a) as the linear combination of a dof basis
(b). The dofs are first evaluated at dof basis (b) (field `dof_basis`) and the
dof values are next mapped to dof basis (a) applying a change of basis (field
`change_of_basis`).

Fields:

- `change_of_basis::Matrix{T}` the matrix of the change from dof basis (b) to (a)
- `dof_basis::AbstractVector{<:Dof}` A type representing dof basis (b)
"""
struct LinearCombinationDofVector{T} <: AbstractVector{Dof}
  change_of_basis::Matrix{T}
  dof_basis::AbstractVector{<:Dof}
end

@inline Base.size(a::LinearCombinationDofVector) = size(a.dof_basis)
@inline Base.axes(a::LinearCombinationDofVector) = axes(a.dof_basis)
@inline Base.getindex(a::LinearCombinationDofVector,i::Integer) = getindex(a.dof_basis,i)
@inline Base.IndexStyle(::LinearCombinationDofVector) = IndexLinear()

function linear_combination(a::AbstractMatrix{<:Number},
                            b::AbstractVector{<:Dof})
  LinearCombinationDofVector(a,b)
end

function linear_combination(a::LinearCombinationDofVector{T},
                            b::AbstractVector{<:Dof}) where T
  linear_combination(a.change_of_basis,b)
end

function return_cache(b::LinearCombinationDofVector,field)
  c, cf = return_cache(b.dof_basis,field)
  c, cf, return_cache(*,b.change_of_basis,c)
end

@inline function evaluate!(cache,b::LinearCombinationDofVector,field)
  c, cf, cc = cache
  vals = evaluate!(cache,b.dof_basis,field)
  evaluate!(cc,*,b.change_of_basis,vals)
end


text:
"""
    (gc, x) = gc_content(region::String, refseq)

Calculate the base composition for a DNA region.
Return values: (GC, [G, C, A, T, missing])
"""
function gc_content(region::String, refseq::NamedTuple)
    chr,range = split(region, ":")
    start,stop = split(range, "-")
    start = parse(Int64, start)
    stop = parse(Int64, stop)

    target = uppercase(refseq[Symbol(chr)][start:stop])

    g = sum( [char for char in target] .== 'G' )
    c = sum( [char for char in target] .== 'C' )
    a = sum( [char for char in target] .== 'A' )
    t = sum( [char for char in target] .== 'T' )
    n = sum( [char for char in target] .== 'N' )

    d = (stop - start + 1) - n
    bases = [g, c, a, t] ./ d
    gc_content = (g + c) / d
    miss = n / (stop - start + 1)
    push!(bases, miss)

    return (gc_content, bases)

end


text:
function parseinput(lines::Vector{String})::Dict{CartesianIndex{3},Int}
    grid = Dict{CartesianIndex{3},Int}()
    for (i, line) in enumerate(lines)
        for (j, char) in enumerate(line)
            grid[CartesianIndex(i, j, 0)] = char == '#'
        end
    end
    return grid
end

function nextstate(
    grid::Dict{CartesianIndex{N},Int},
    pos::CartesianIndex{N},
    neighbors::Vector{CartesianIndex{N}}) where N
    state = get(grid, pos, 0)
    if state == 1
        return sum(get(grid, pos + n, 0) for n in neighbors) ∈ [2, 3] ? 1 : 0
    else
        return sum(get(grid, pos + n, 0) for n in neighbors) == 3 ? 1 : 0
    end
end

function oneround!(
    grid::Dict{CartesianIndex{N},Int},
    neighbors::Vector{CartesianIndex{N}}) where N
    prev_grid = copy(grid)
    for k in keys(prev_grid)
        for n in neighbors
            grid[k + n] = nextstate(prev_grid, k + n, neighbors)
        end
    end
end

const neighbors_p1 = [
    CartesianIndex(i, j, k)
    for i = -1:1 for j = -1:1 for k = -1:1
    if !(i == j == k == 0)
]

const data = readlines("day17.txt")
const grid_p1 = parseinput(data)

function run!(grid::Dict{CartesianIndex{N},Int}, neighbors::Vector{CartesianIndex{N}}) where N
    for i in 1:6
        oneround!(grid, neighbors)
    end
    return sum(values(grid))
end

println("Part 1: $(run!(grid_p1, neighbors_p1))")

function parseinput_p2(lines::Vector{String})::Dict{CartesianIndex{4},Int}
    grid = Dict{CartesianIndex{4},Int}()
    for (i, line) in enumerate(lines)
        for (j, char) in enumerate(line)
            grid[CartesianIndex(i, j, 0, 0)] = char == '#'
        end
    end
    return grid
end

const neighbors_p2 = [
    CartesianIndex(i, j, k, l)
    for i = -1:1 for j = -1:1 for k = -1:1 for l = -1:1
    if !(i == j == k == l == 0)
]

const grid_p2 = parseinput_p2(data)

println("Part 2: $(run!(grid_p2, neighbors_p2))")


text:
using Meshing
using GeometryTypes
using LinearAlgebra: dot, norm
using FileIO

f(v) = sqrt(sum(dot(v,v))) - 1
sdf = SignedDistanceField(f,HyperRectangle(Vec(-1,-1,-1.), Vec(2,2,2.)))

mc = HomogenousMesh(sdf, MarchingCubes())
mt = HomogenousMesh(sdf, MarchingTetrahedra())
ns = HomogenousMesh(sdf, NaiveSurfaceNets())

# save the Sphere as a PLY file
save("sphere_mc.ply",mc)
save("sphere_mt.ply",mt)
save("sphere_ns.ply",ns)


text:
#!/usr/bin/julia

# Simulate the toppling of sandpiles.

# See also:
#   https://en.wikipedia.org/wiki/Abelian_sandpile_model
#   https://www.youtube.com/watch?v=1MtEUErz7Gg -- Sandpiles - Numberphile
#   https://www.youtube.com/watch?v=diGjw5tghYU -- Coding Challenge #107: Sandpiles (by Daniel Shiffman)

using Images

function topple(plane, nextplane, width, height)

    for y in 1:height, x in 1:width
        pile = plane[y,x]

        if (pile < 4)
            nextplane[y,x] = pile
        end
    end

    for y in 2:height-1, x in 2:width-1
        pile = plane[y,x]

        if (pile >= 4)
            nextplane[y, x] += pile - 4
            nextplane[y - 1, x] += 1
            nextplane[y + 1, x] += 1
            nextplane[y, x - 1] += 1
            nextplane[y, x + 1] += 1
        end
    end

    return nextplane
end

function contains_four(plane, width, height)

    for y in 1:height, x in 1:width
        if (plane[y,x] >= 4)
            return true
        end
    end

    return false
end

function sandpiles()

    w, h = 240, 240

    plane = zeros(Int64, h, w)
    plane[div(h,2), div(w,2)] = 10^5

    while (contains_four(plane, w, h))
        nextplane = zeros(Int64, h, w)
        plane = topple(plane, nextplane, w, h)
    end

    colors = [
                RGB{Float64}(0, 0, 0),
                RGB{Float64}(0, 0, 255),
                RGB{Float64}(0, 255, 0),
                RGB{Float64}(255, 255, 255),
            ]

    background = RGB{Float64}(0, 0, 0)
    img = zeros(RGB{Float64}, h, w)

    for y in 1:h, x in 1:w
        pile = plane[y,x]
        col = background

        if (pile <= 3)
            col = colors[pile+1]
        end

        img[y,x] = col
    end

    img = map(clamp01nan, img)

    return img
end

img = sandpiles()
println("Saving image...")
save("sandpiles.png", img)

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
@doc doc"""
Minimise le problème : ``min_{||s||< \delta_{k}} q_k(s) = s^{t}g + (1/2)s^{t}Hs``
                        pour la ``k^{ème}`` itération de l'algorithme des régions de confiance

# Syntaxe
```julia
sk = Gradient_Conjugue_Tronque(fk,gradfk,hessfk,option)
```

# Entrées :
   * **gradfk**           : (Array{Float,1}) le gradient de la fonction f appliqué au point xk
   * **hessfk**           : (Array{Float,2}) la Hessienne de la fonction f appliqué au point xk
   * **options**          : (Array{Float,1})
      - **delta**    : le rayon de la région de confiance
      - **max_iter** : le nombre maximal d'iterations
      - **tol**      : la tolérance pour la condition d'arrêt sur le gradient

# Sorties:
   * **s** : (Array{Float,1}) le pas s qui approche la solution du problème : ``min_{||s||< \delta_{k}} q(s)``

# Exemple d'appel:
```julia
gradf(x)=[-400*x[1]*(x[2]-x[1]^2)-2*(1-x[1]) ; 200*(x[2]-x[1]^2)]
hessf(x)=[-400*(x[2]-3*x[1]^2)+2  -400*x[1];-400*x[1]  200]
xk = [1; 0]
options = []
s = Gradient_Conjugue_Tronque(gradf(xk),hessf(xk),options)
```
"""
function Gradient_Conjugue_Tronque(gradfk,hessfk,options)

    "# Si option est vide on initialise les 3 paramètres par défaut"
    if options == []
        deltak = 2
        max_iter = 50
        tol = 1e-6
    else
        deltak = options[1]
        max_iter = options[2]
        tol = options[3]
    end

   n = length(gradfk)
   s = zeros(n)
   sortir = false

#Initisalisation
j = 0
sj = zeros(n)
s = zeros(n)
gj = gradfk
pj = -gj
q(s) = transpose(gradfk)*s + (1/2) * transpose(s)*hessfk*s

while j <= max_iter && sortir == false

 if norm(gj) != 0
    kj = transpose(pj)*hessfk*pj

    discriminant = (transpose(pj)*sj)^2 + (norm(pj)^2)*(deltak^2 - norm(sj)^2)
    sigma1 = ((-transpose(sj)*pj)+sqrt(discriminant))/(norm(pj)^2)
    sigma2 = ((-transpose(sj)*pj)-sqrt(discriminant))/(norm(pj)^2)

    if kj <= 0
        if q(sj + sigma1*pj) <= q(sj + sigma2*pj)
            s = sj + sigma1 * pj
            #s = sj
        else
            s = sj + sigma2 * pj
            #s = sj
        end
        sortir = true
    end

  if sortir == false
    alphaj = (transpose(gj)*gj)/kj
    if norm(sj + alphaj * pj) >= deltak
        if sigma1 >= 0
            s = sj + sigma1 * pj
        elseif sigma2 >= 0
            s = sj + sigma2 * pj
        end
        sortir = true
    end
  end

  if sortir == false

        sj += alphaj * pj
        gj0 =  gj
        gj += alphaj * hessfk * pj
        betaj = (transpose(gj)*gj)/(transpose(gj0)*gj0)
        pj  = -gj + betaj * pj

        if norm(gj) < tol * norm(gradfk)
            s = sj
            sortir = true
        end
  end

 end
  j = j + 1
end
   return s
end

max_stars_repo_name: faicaltoubali/ENSEEIHT

max_issues_repo_name: faicaltoubali/ENSEEIHT

max_forks_repo_name: faicaltoubali/ENSEEIHT


text:
# ---
# title: 1442. Count Triplets That Can Form Two Arrays of Equal XOR
# id: problem1442
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Math, Bit Manipulation
# link: <https://leetcode.com/problems/count-triplets-that-can-form-two-arrays-of-equal-xor/description/>
# hidden: true
# ---
#
# Given an array of integers `arr`.
#
# We want to select three indices `i`, `j` and `k` where `(0 <= i < j <= k <
# arr.length)`.
#
# Let's define `a` and `b` as follows:
#
#   * `a = arr[i] ^ arr[i + 1] ^ ... ^ arr[j - 1]`
#   * `b = arr[j] ^ arr[j + 1] ^ ... ^ arr[k]`
#
# Note that **^** denotes the **bitwise-xor** operation.
#
# Return _the number of triplets_ (`i`, `j` and `k`) Where `a == b`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [2,3,1,6,7]
#     Output: 4
#     Explanation: The triplets are (0,1,2), (0,2,2), (2,3,4) and (2,4,4)
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,1,1,1,1]
#     Output: 10
#
#
# **Example 3:**
#
#
#
#     Input: arr = [2,3]
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: arr = [1,3,5,7,9]
#     Output: 3
#
#
# **Example 5:**
#
#
#
#     Input: arr = [7,11,12,9,5,2,7,17,22]
#     Output: 8
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 300`
#   * `1 <= arr[i] <= 10^8`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function Base.:^(a::UncertainIndexDataset, b::UncertainIndexDataset; n = 30000)
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainIndexDataset([a[i] ^ b[i] for i = 1:N])
end

function Base.:^(a::T, b::UncertainIndexDataset; n = 30000) where {T <: Number}
    N = length(a)
    UncertainIndexDataset([a ^ b[i] for i = 1:N])
end

function Base.:^(a::Vector{T}, b::UncertainIndexDataset; n = 30000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainIndexDataset([a[i] ^ b[i] for i = 1:N])
end

function Base.:^(a::UncertainIndexDataset, b::T; n = 30000) where {T <: Number}
    N = length(a)

    UncertainIndexDataset([a[i] ^ b for i = 1:N])
end

function Base.:^(a::UncertainIndexDataset, b::Vector{T}; n = 30000) where {T <: Number}
    N = length(a)
    n_vals_b = length(b)

    if N != n_vals_b
        throw(ArgumentError("Dataset lengths do not match ($N, $n_vals_b)"))
    end

    UncertainIndexDataset([a[i] ^ b[i] for i = 1:N])
end


text:
using Test
@test ENV["A"] == ENV["B"]


text:
export isIndependent, isRelevant
export independencyIndex, toroidalIndex, divisorialIndex, minMaxDivisorial
export remove!, getIndex, coneContains
export BerghC

"""
    remove!(::Array{Int64},::Int64)

    In-place removes a given item from a vector.

# Examples
```jldoctest
julia> A=[1,2,3,4];

julia> remove!(A,1);

julia> A
[2,3,4]
```
"""
function remove!(a::Array{Int64,1}, item::Int64)
    return deleteat!(a, findall(x->x==item, a))
end

"""
    getIndex(::Array{Int64,1},::Array{Int64,2})

    Returns the first index at which a vector appears as a row of a matrix.

# Examples
```jldoctest
getIndex([0,1,0],[1 0 0; 0 1 0; 0 0 1])
2
```
"""
function getIndex(ray::Array{Int64,1},rayMatrix::Array{Int64,2})
    slice=slicematrix(rayMatrix)
    index=findfirst(x->x==ray,slice)
    return index
end

"""
    isIndependent(::Int64,::Array{Int64,1},::Array{Int64,2})

    Takes a ray matrix, a list of indices representing a cone, and an index represeting a ray of that cone. Determines whether the given ray is independent in the cone (i.e. does not contribute to the multiplicity of the cone).

# Examples
```jldoctest
julia> isIndependent(3,[1,2,3],[1 0 0; 0 1 0; 1 2 3])
false

julia> isIndependent(3,[1,2,3],[1 0 0; 0 1 0; 1 1 1])
true
```
"""
function isIndependent(rayIndex::Int64,cone::Array{Int64,1},rayMatrix::Array{Int64,2})
    if size(cone,1)==1
        return true
    end
    scone=copy(cone)
    subcone=remove!(scone,rayIndex)
    mult=getMultiplicity(cone,rayMatrix)
    submult=getMultiplicity(subcone,rayMatrix)
    return mult==submult
end

"""
    independencyIndex(::Array{Int64,1},::Array{Int64,2})

    Returns the number of non-independent rays in a cone. Input in indices-ray matrix format.

# Examples
```jldoctest
julia> independencyIndex([1,2,3],[1 0 0 ; 1 2 0; 2 0 3; 0 0 5])
2
```
"""
function independencyIndex(cone::Array{Int64,1},rayMatrix::Array{Int64,2})
    return count(elt -> !isIndependent(elt, cone, rayMatrix), cone)
end

"""
    isRelevant(::Array{Int64,1},::Array{Int64,1},::StackyFan)

    Determines if the given ray, relative to the given cone, either has a stacky value greater than 1 or is not independent.

# Examples
```jldoctest
julia> F=makeStackyFan([1 0 0; 1 2 0; 0 0 1],[[0,1,2]],[1,1,2]);

julia> isRelevant([1,2,0],[1,2,3],F)
true

julia> F=makeStackyFan([1 0 0; 0 1 0; 0 0 1],[[0,1,2]],[1,1,2]);

julia> isRelevant([0,1,0],[1,2,3],F)
false
```
"""
function isRelevant(ray::Array{Int64,1},cone::Array{Int64,1},F::StackyFan)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(F.fan.RAYS)))
    rayStack=F.stacks[encode(ray)]
    rayIndex=getIndex(ray,rayMatrix)
    rayIndependent=isIndependent(rayIndex,cone,rayMatrix)
    return rayStack != 1 || !rayIndependent
end

"""
    toroidalIndex(::Array{Int64,1},::StackyFan,::Dict)

    Calculates the toroidal index of the given cone of a divisorial stacky fan, or the number of relevant non-divisorial rays. Compare to divisorialIndex().

# Examples
```jldoctest
julia> F=makeStackyFan([1 0 0; 0 1 0; 1 0 2],[[0,1,2]],[1,2,1]);

julia> div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>0);

julia> toroidalIndex([1,2,3],F,div)
3
julia> div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>1);

julia> toroidalIndex([1,2,3],F,div)
2
```
"""
function toroidalIndex(cone::Array{Int64,1},F::StackyFan,div::Dict)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(F.fan.RAYS)))
    slice=slicematrix(rayMatrix)
    # Find number of divisorial rays
    # div is a dictionary that represents which rays are divisorial, 0 represents a non-divisorial ray and 1 represents a divisorial ray
    s=count(x->div[slice[x]]==1,cone)
    #flipt counts the number of non-divisorial and irrelevant rays
    flipt = count(i -> div[slice[i]]==0 && !isRelevant(slice[i], cone, F), cone)
    # number of relevant cones
    t=size(cone,1)-flipt
    # return the toroidal index, which is the number of relevant residual (non-divisorial) rays
    return t-s
end

"""
    divisorialIndex(::Array{Int64,1},::StackyFan,::Dict)

    Calculates the divisorial index (defined by Daniel Bergh) of a given cone in a fan with divisorial rays.
    Specifically, takes the subcone consisting of all relevant non-divisorial rays in a cone, and counts the number of rays that are relevant in that subcone.

# Examples
```jldoctest
julia> F=makeStackyFan([1 0 0; 0 1 0; 1 0 2],[[0,1,2]],[1,2,1]);

julia> div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>0);

julia> divisorialIndex([1,2,3],F,div)
3
julia> div=Dict([1,0,0]=>0,[0,1,0]=>0,[1,0,2]=>1);

julia> divisorialIndex([1,2,3],F,div)
1
```
"""
function divisorialIndex(cone::Array{Int64,1},F::StackyFan,div::Dict)
    slicedRayMatrix=slicematrix(convert(Array{Int64,2},Array(Polymake.common.primitive(F.fan.RAYS))))
    # relevant residual cones
    relRes=Array{Int64,1}[]
    relResStack=Int64[]
    # c is the number of non-divisorial relevant cones
    c=0
    for i in cone
        ray=slicedRayMatrix[i]
        stack=F.stacks[encode(ray)]
        # If the ray is non-divisorial and relevant, increment c by one, add ray to relRes, and stack to relResStack. The rays in relRes are used to build a new cone.
        if div[ray]==0 && isRelevant(ray,cone,F)
            c+=1
            push!(relRes,ray)
            push!(relResStack,stack)
        end
    end
    # If there are no relevant residual cones in F, the divisorial index is 0
    if c==0
        return 0
    else
        # convert to 0-indexing
        relResIndz=[[i-1 for i in 1:c]]
        relResInd=[i for i in 1:c]
        relResCat=Array{Int64}(transpose(hcat(relRes...)))
        # Construct a subfan consisting of all relevant residual rays
        subfan=makeStackyFan(relResCat,relResIndz,relResStack)
        divInd=0
        for ray in relRes
            # Iterate over relevant residual cones to count the number of relevant rays in the subfan
            # This count represents the divisorial index
            if isRelevant(ray,relResInd,subfan)
                divInd+=1
            end
        end
        return divInd
    end
end

"""
    coneContains(::Array{Int64,1},::Array{Int64,1})

    Checks whether every index in the first input is also contained in the second input.

# Examples
```jldoctest
julia> coneContains([1,2,3],[1,2,3,4])
true
julia> coneContains([1,2,5],[1,2,3,4])
false
```
"""
function coneContains(A::Array{Int64,1},B::Array{Int64,1})
    return issubset(A, B)
end

"""
    minMaxDivisorial(::StackyFan,::Dict)

    Calculates the maximal divisorial index of all cones in a stacky fan.
    Each maximal cone of the fan will contain at most one minimal subcone of maximal divisorial index; a list of such cones is returned.

# Examples
```jldoctest
julia> F=makeStackyFan([1 2 0;1 3 0; 3 0 1],[[0,1,2]],[1,1,5]);

julia> div=Dict([1,2,0]=>0,[1,3,0]=>0,[3,0,1]=>0);

julia> minMaxDivisorial(F,div)
[[3]]

julia> F=makeStackyFan([1 1 0;1 3 0; 3 0 1],[[0,1,2]],[1,1,5]);

julia> div=Dict([1,1,0]=>0,[1,3,0]=>0,[3,0,1]=>0);

julia> minMaxDivisorial(F,div)
[[1,2,3]]
```
"""
function minMaxDivisorial(F::StackyFan,div::Dict)
    # Calculates the maximal divisorial index of any cone in the fan
    divMax=0
    coneList=getCones(F.fan)
    # dictionary that represents each cone with its divisorial index
    divisorialDict=Dict()
    for cone in coneList
        d=divisorialIndex(cone,F,div)
        divisorialDict[cone]=d
        if d>divMax
            divMax=d
        end
    end
    if divMax==0
        return nothing
    end

    # cones with maximal divisorial index
    divMaxCones=Array{Int64,1}[]
    for cone in coneList
        if divisorialDict[cone]==divMax
             # if the cone's divisorial index is the fan's maximal divisorial index, add the cone to divMaxCones
            push!(divMaxCones,cone)
        end
    end

    #divMaxConesRefined stores the cones in divMaxCones that are minimal with respect to inclusion
    divMaxConesRefined=Array{Int64,1}[]
    # List of maximal cones in F
    maxconeList=convertIncidenceMatrix(F.fan.MAXIMAL_CONES)
    for maxcone in maxconeList
        # if the div index of the current maxcone is the fan's max div index, its minimal subcone with maximal divisorial index is calculated
        if divisorialDict[maxcone]==divMax
            maxconeContains=Array{Int64,1}[]
            mincone=maxcone
            for cone in divMaxCones
                if coneContains(cone,maxcone) && size(cone,1)<size(mincone,1)
                    mincone=cone
                end
            end
            if !(mincone in divMaxConesRefined)
                push!(divMaxConesRefined,mincone)
            end
        end
    end
    return divMaxConesRefined
end

"""
    BerghC(::StackyFan,::Array{Int64,1})

    Takes a stacky fan and a binary array indicating which rays are divisorial, and runs Daniel Bergh's algorithm C.
    This algorithm performs a series of stacky blowups to reduce the maximal divisorial index of the fan, and returns a fan with a maximal divisorial index of 0.

# Examples
```jldoctest
julia> F=makeStackyFan([1 1 0;1 3 0; 0 0 1],[[0,1,2]],[1,1,5]);

julia> H, div = BerghC(F,[0,0,0]);

julia> convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))
5×3 Matrix{Int64}:
 1  1  0
 0  0  1
 2  4  5
 1  3  0
 1  2  0

julia> div
Dict{Any, Any} with 5 entries:
  [1, 1, 0] => 0
  [0, 0, 1] => 1
  [2, 4, 5] => 1
  [1, 2, 0] => 1
  [1, 3, 0] => 0

julia> F=makeStackyFan([1 0;1 3; 5 17],[[0,1],[1,2]],[1,1,5]);

julia> H, div = BerghC(F,[0,0,0]);

julia> convert(Array{Int64,2},Polymake.common.primitive(H.fan.RAYS))
5×2 Matrix{Int64}:
  1   0
  2   3
  1   3
 13  44
  5  17
julia> div
Dict{Any, Any} with 5 entries:
  [1, 0]   => 0
  [13, 44] => 1
  [1, 3]   => 0
  [2, 3]   => 1
  [5, 17]  => 1
```
"""
function BerghC(F::StackyFan,divlist::Array{Int64,1})
    X=deepcopy(F)
    rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS)))
    slicedRayMatrix=slicematrix(rayMatrix)
    div=Dict()
    # Populate dictionary of divisorial rays (div) from divlist
    for i in 1:size(slicedRayMatrix,1)
        div[slicedRayMatrix[i]]=divlist[i]
    end
    while(true)
        rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS)))
        slicedRayMatrix=slicematrix(rayMatrix)
        # Find the cones with maximal divisorial index in X
        subdivTargetCones=minMaxDivisorial(X,div)
        # If there are no such cones, the algorithm terminates
        if subdivTargetCones==nothing
            break
        end
        blowupList=Array{Array{Int64,1},1}[]
        # Iterate through cones with maximal divisorial index
        for cone in subdivTargetCones
            # Add each cone's ray representation to blowupList
            push!(blowupList,slicematrix(rowMinors(rayMatrix,cone)))
        end
        for raycone in blowupList
            indices=Int64[]
            slicedRayMatrix=slicematrix(convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS))))
            for ray in raycone
                push!(indices,findall(x->x==ray,slicedRayMatrix)[1])
            end
            # Sort indices in ascending order
            cone=sort(indices)
            if size(cone,1)==1
                div[slicedRayMatrix[cone[1]]]=1
            else
                exceptional=findStackyBarycenter(cone,X)
                # perform the blowup
                X=stackyBlowup(X,[x-1 for x in cone], exceptional)
                # convert exceptional ray to its primitive form
                primExcep=Array{Int64,1}(Polymake.common.primitive(exceptional))
                div[primExcep]=1
            end
        end
    end
    return X, div
end

"""
    BerghC(::StackyFan,::Dict)

    Functions identically to the previous BerghC method, but takes a dictionary that assigns 0 or 1 to the rays of the input fan to indicate divisoriality. In the previous BerghC method, this dictionary is computed.

"""
function BerghC(F::StackyFan,div::Dict)
    X=deepcopy(F)
    while(true)
        rayMatrix=convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS)))
        slicedRayMatrix=slicematrix(rayMatrix)
        subdivTargetCones=minMaxDivisorial(F,div)
        if subdivTargetCones==nothing
            break
        end
        blowupList=Set[]
        for cone in subdivTargetCones
            push!(blowupList,slicematrix(rowMinors(rayMatrix,cone)))
        end
        for raycone in blowupList
            indices=Int64[]
            slicedRayMatrix=slicematrix(convert(Array{Int64,2},Array(Polymake.common.primitive(X.fan.RAYS))))
            for ray in raycone
                push!(indices,findall(x->x==ray,slicedRayMatrix)[1])
            end
            cone=sort(indices)
            if size(cone,1)==1
                div[slicedRayMatrix[cone[1]]]=1
            else
                exceptional=findStackyBarycenter(cone,X)
                X=stackyBlowup(X,[x-1 for x in cone], exceptional)
                div[exceptional]=1
            end
        end
    end
    return X, div
end

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
using LinearAlgebra
using Printf

include("back_substitution.jl")
include("qr_factorization.jl")

"""
Finds UV decomposition of rank k that minimizes f(U_i,V_i) = ||A - U_iV_i||_f by alternating minimization

# Arguments
- `A::Array{Float64,2}`: matrix of size m x n.
- `k::Integer`: the rank of the UV decomposition. k must be <= rank(A).
- `initV::Array{Float64,2}=rand(k, size(A)[2])`: initial V_0 from which the algorithm starts. V_0 must be full rank.
- `min_error::Integer=1e-8`: the algorithm will stop if f(U_i,V_i) - f(U_i+1,V_i+1) <= min_error.
- `max_iterations::Integer=100`: the algorithm will perform at most max iteration steps.
- `print_error::Bool=true`: if it is true, the algorithm will print the error at every iteration.
- `save errors::Bool=false`: if it is true, the algorithm will return an array containing the history of the errors, besides U and V
- `opt::Array{Float64,2}=A`: optimal factorization of rank k used to compute the residuals if save_errors is true. (Used only for tests)

# Return
- `U::Array{Float64,2}`: matrix of size m x k.
- `V::Array{Float64,2}`: matrix of size k x n.
- `errors::Array{Float64,1}`: only if print_error is true.
- `residuals::Array{Float64,1}`: only if print_error is true.
"""
function cresetqr(A, k; initV = rand(k, size(A)[2]), min_error=1e-8, max_iterations=100, print_error=true, save_errors=false, opt=A)
    (m, n) = size(A)

    # allocating matrices U and V in which the program will allocate the results
    U = zeros(m, k)
    V = zeros(k,n)
    V .= initV

    num_step = 0

    # number of iterations the algorithm have perfomed
    num_steps = 1

    # allocating matrices in which store QR factorization and back_substitution results
    QR_V, v_V, u_V = allocate_matrices(V')
    QR_U, v_U, u_U = allocate_matrices(U)
    W_V = zeros(n, m)
    Q_U = zeros(m, k)

    # computing the initial function value
    f = norm(A - U * V, 2)

    # allocating vectors in which store the error and the residuals during the algorithm execution
    if save_errors
        errors = zeros(max_iterations)
        residuals = zeros(max_iterations)
    end

    if print_error
        println("step\t ||A - UV||_f\t ||f_old-f||")
    end

    @views @inbounds for num_step = 1:max_iterations
        # gradient = [U*V*V' - A*V', U'*U*V - U'*A]

        # Q,R = qr(V')
        qr_factorization!(V', QR_V, v_V, u_V)

        # taking the R factor
        R = triu(QR_V)[1:k,:]

        # W = Q' * A'
        Q_t_times_A!(QR_V, A', W_V)

        # R * U' = W_V by back substitution
        back_substitution!(R, U', W_V[1:k,1:m])

        # Q,R = qr(U)
        qr_factorization!(U, QR_U, v_U, u_U)

        # taking the R factor
        R = triu(QR_U)[1:k,:]

        # taking the Q factor
        get_Q!(QR_U, Q_U)

        # R * V = Q_U' * A by back substitution
        back_substitution!(R, V, Q_U' * A)

        #additional steps for CRESETQR (w.r.t. RESETQR)
        U .= Q_U[:,1:k]
        V .= R * V

        #computing new function value
        f_old = f
        f = norm(A - U * V, 2)
        delta = norm(f_old - f, 2)

        if print_error
            @printf("%d)\t%.13f\t%.13f\n",num_step,f,delta)
        end

        #saving function value and residual
        if save_errors
            errors[num_step] = f
            residuals[num_step] = norm(U * V - opt)
            num_steps = num_step
        end

        #Stopping criteria check
        if delta <= min_error
            if print_error
                println("Reached good approximation with ", num_step, " iterations")
            end
            break
        end
    end

    if save_errors
        return U, V, errors[1:num_steps], residuals[1:num_steps]
    else
        return U, V
    end
end


text:
using Documenter#, FeedbackConvNets

makedocs(
    sitename = "FeedbackConvNets Documentation",
    #modules = [FeedbackConvNets],
    pages = [
        "Home" => "index.md"
        "Guide" => [
            "Getting Started" => "guide/getting_started.md",
            "Chains vs Trees" => "guide/chains_vs_trees.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/cJarvers/FeedbackNets.jl.git"
)


text:
# %% Using Cao's method to estimate embedding
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

lo = Systems.lorenz([0, 10, 0.0])

tr = trajectory(lo, 1000; Ttr=10)
x, y, z = columns(tr)
fig = figure(figsize=(figx/2, figy));
w = x
ψ = z .- y

for (i, (s, l)) in enumerate(zip((w, ψ), ("\$w=x\$", "\$w=z-y\$")))
    τ = estimate_delay(s, "mi_min")
    Ed = delay_afnn(s, τ, 2:7)
    plot(2:7, Ed, marker = "o", label = l, color = "C$([0, 2][i])")
end

xlabel("\$d\$"; labelpad = -10)
ylabel("\$E_{d}\$")
legend(title="measurement")
fig.tight_layout(pad = 0.4)
fig.subplots_adjust(left = 0.15, bottom = 0.15)
wsave(plotsdir("6", "caodemonstration"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


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

# %%
Fact = (u -> u(u))(u -> n -> n == 0 ? 1 : n * u(u)(n - 1))
@show [Fact(n) for n in 0:10];

# %%
f(u, n) = n == 0 ? 1 : n * u(u, n - 1)
@show [f(f, n) for n in 0:10];

# %%
@code_warntype f(f, 5)

# %%
@code_typed f(f, 5)

# %%
@code_llvm debuginfo=:none f(f, 5)

# %%
@code_native debuginfo=:none f(f, 5)

# %%
g(u, n, a, b) = n == 0 ? a : u(u, n-1, b, a+b)
@show [g(g, n, 0, 1) for n in 0:10];

# %%
U(u) = u(u)

# %%
U(U)

# %%

max_stars_repo_path: 0013/recursive computation by U-combinator.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# File copied from SMC.jl
"""
```
function scalar_reduce(args...)
```
Each individual iteration returns n scalars. The output is reduced to n vectors,
where the i-th vector contains all of the i-th scalars from each iteration.

The return type of reduce functions must be the same type as the tuple of
arguments passed in. If args is a tuple of Vector{Float64}, then the return
argument will be a Vector{Float64}.

e.g.
a, b = @parallel (scalar_reduce) for i in 1:10000
           [[1], [2]]
       end
a = [1, 1, 1, ...]
b = [2, 2, 2, ...]

Input/Output type: Vector{Vector{Float64}}
"""
function scalar_reduce(args...)
    return_arg = args[1]
    for (i, arg) in enumerate(args[2:end])
        for (j, el) in enumerate(arg)
            append!(return_arg[j], el)
        end
    end
    return return_arg
end

"""
```
function vector_reduce(args...)
```
Each individual iteration returns n Vector types; we vector-reduce to n matrices, where
the i-th column of that matrix corresponds to the i-th vector from an individual iteration.

Input/Output type: Vector{Matrix{Float64}}
"""
function vector_reduce(args...)
    nargs1 = length(args)    # The number of times the loop is run
    nargs2 = length(args[1]) # The number of variables output by a single run

    return_arg = args[1]
    for i in 1:nargs2
        for j in 2:nargs1
            return_arg[i] = hcat(return_arg[i], args[j][i])
        end
    end
    return return_arg
end

"""
```
function vec_red_scal(args...)
```
vector_reduce but scalar_reduce for last argument
Hard coded for 3 arguments
"""
function vec_scal_reduce(args...)
    nargs1 = length(args)    # The number of times the loop is run
    nargs2 = length(args[1]) # The number of variables output by a single run

    return_arg = args[1]
    for i in 1:nargs2
        for j in 2:nargs1
            if i == nargs2
                append!(return_arg[i], args[j][i][1])
                #return_arg[i] = vcat(return_arg[i], args[j][i])
            else
                return_arg[i] = hcat(return_arg[i], args[j][i])
            end
        end
    end
    return return_arg
#=
    nargs1 = length(args)    # The number of times the loop is run
    nargs2 = length(args[1]) # The number of variables output by a single run

    if nargs2 == 1
        return args
    end

    arg1, arg2, arg3 = args[1]
    #arg2 = args[1][2]
    #arg3 = args[1][3]

    for j in 2:nargs1
        arg1 = hcat(arg1, args[j][1])
        arg2 = hcat(arg2, args[j][2])
        arg3 = vcat(arg3, args[j][3])
    end

    return_arg = (arg1, arg2, arg3)
    return return_arg
=#
end

"""
```
function scalar_reshape(args...)
```
Function ensures type conformity of the return arguments.
"""
function scalar_reshape(args...)
    n_args = length(args)
    return_arg = Vector{Vector{Float64}}(undef, n_args)
    for i in 1:n_args
        arg = typeof(args[i]) <: Vector ? args[i] : [args[i]]
        return_arg[i] = arg
    end
    return return_arg
end

"""
```
function vector_reshape(args...)
```
Function ensures type conformity of the return arguments.
"""
function vector_reshape(args...)
    n_args = length(args)
    return_arg = Vector{Matrix{Float64}}(undef, n_args)
    for i in 1:n_args
        arg = typeof(args[i]) <: Vector ? args[i] : [args[i]]
        return_arg[i] = reshape(arg, length(arg), 1)
    end
    return return_arg
end

"""
```
sendto(p::Int; args...)
```
Function to send data from master process to particular worker, p.
Code from ChrisRackauckas, avavailable at:
 https://github.com/ChrisRackauckas/ParallelDataTransfer.jl/blob/master/src/ParallelDataTransfer.jl.
"""
function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, Core.eval(Main, Expr(:(=), nm, val)))
    end
end

"""
```
sendto(ps::AbstractVector{Int}; args...)
```
Function to send data from master process to list of workers.
Code from ChrisRackauckas, available at:
https://github.com/ChrisRackauckas/ParallelDataTransfer.jl/blob/master/src/ParallelDataTransfer.jl.
"""
function sendto(ps::AbstractVector{Int}; args...)
    for p in ps
        sendto(p; args...)
    end
end


text:
@everywhere pri = Factored(Normal(1, 0.5), DiscreteUniform(1, 10))
@everywhere sim((n, du)) = (n * n + du) * (n + randn() * 0.01)
@everywhere cost(x) = abs(sim(x) - 5.5)
@everywhere model_abc = ApproxPosterior(pri, cost, 0.01)
@everywhere res = sample(model_abc, AIS(100), 1000, discard_initial = 10000, progress = false)
@show sim(Tuple(res)) ≈ 5.5
res = smc(pri, cost; verbose = true)
display(res)
@show res.P[2] ≈ 5
@everywhere function brownianrms((μ, σ), N, samples = 200)
    t = 0:N
    #    rand()<1/20 && sleep(0.001)
    @.(sqrt(μ * μ * t * t + σ * σ * t)) .* (0.95 + 0.1 * rand())
end

@everywhere params = (0.5, 2.0)
@everywhere tdata = brownianrms(params, 30, 10000)
@everywhere prior = Factored(Uniform(0, 1), Uniform(0, 4))
@everywhere cost(x) = sum(abs, brownianrms(x, 30) .- tdata) / length(tdata)
@everywhere modelabc = ApproxPosterior(prior, cost, 0.1)
all(smc(prior, cost, verbose=false,min_r_ess=0.55).P .≈ params)


text:
using Triangulate
using WriteVTK
using CSV, DataFrames
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../../test/useful_routines.jl")
include("../../compute-stress.jl")

function lame_lambda(k, m)
    return k - 2m / 3
end

function bulk_modulus(l, m)
    return l + 2m / 3
end

function onleftboundary(x, L, W)
    return x[1] ≈ 0.0
end

function onbottomboundary(x, L, W)
    return x[2] ≈ 0.0
end

function onrightboundary(x, L, W)
    return x[1] ≈ L
end

function ontopboundary(x, L, W)
    return x[2] ≈ W
end

function onboundary(x, L, W)
    return onleftboundary(x, L, W) ||
           onbottomboundary(x, L, W) ||
           onrightboundary(x, L, W) ||
           ontopboundary(x, L, W)
end

function solve_incoherent_cube_and_compute_stress(
    width,
    corner,
    stiffness,
    theta0,
    initial_pressure,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

    lambda1, mu1 = CutCell.lame_coefficients(stiffness, +1)
    lambda2, mu2 = CutCell.lame_coefficients(stiffness, -1)
    transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

    dx = width / nelmts
    meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
    penalty = penaltyfactor / dx * meanmoduli

    basis = TensorProductBasis(2, polyorder)
    mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
    levelset = InterpolatingPolynomial(1, basis)
    levelsetcoeffs =
        CutCell.levelset_coefficients(x -> -corner_distance_function(x, corner), mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.incoherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    leftdisplacementbc = CutCell.DisplacementComponentCondition(
        x -> 0.0,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onleftboundary(x, width, width),
        [1.0, 0.0],
        penalty,
    )
    bottomdisplacementbc = CutCell.DisplacementComponentCondition(
        x -> 0.0,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onbottomboundary(x, width, width),
        [0.0, 1.0],
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        leftdisplacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(
        sysmatrix,
        sysrhs,
        bottomdisplacementbc,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_component_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onleftboundary(x, width, width),
        [1.0, 0.0],
    )
    CutCell.assemble_penalty_displacement_component_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onbottomboundary(x, width, width),
        [0.0, 1.0],
    )

    CutCell.assemble_traction_force_component_linear_form!(
        sysrhs,
        x -> -initial_pressure,
        basis,
        facequads,
        cutmesh,
        x -> onrightboundary(x, width, width),
        [1.0, 0.0],
    )
    CutCell.assemble_traction_force_component_linear_form!(
        sysrhs,
        x -> -initial_pressure,
        basis,
        facequads,
        cutmesh,
        x -> ontopboundary(x, width, width),
        [0.0, 1.0],
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    displacement = matrix \ rhs

    qpstress = compute_stress_at_quadrature_points(
        displacement,
        basis,
        stiffness,
        transfstress,
        theta0,
        cellquads,
        cutmesh,
    )
    qpcoords = compute_quadrature_points(cellquads, cutmesh)
    nodal_displacement = reshape(displacement, 2, :)

    return nodal_displacement, qpstress, qpcoords
end

K1, K2 = 247.0, 192.0
mu1, mu2 = 126.0, 87.0
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)

theta0 = -0.067
initial_pressure = 13.5
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

width = 1.0
corner = [0.8, 0.8]
penaltyfactor = 1e2

nelmts = 37
polyorder = 3
numqp = required_quadrature_order(polyorder) + 2

nodal_displacement, qpstress, qpcoords = solve_incoherent_cube_and_compute_stress(
    width,
    corner,
    stiffness,
    theta0,
    initial_pressure,
    nelmts,
    polyorder,
    numqp,
    penaltyfactor,
)

pressure = -(qpstress[1, :] + qpstress[2, :] + qpstress[4, :]) / 3

triin = Triangulate.TriangulateIO()
triin.pointlist = qpcoords
(triout, vorout) = triangulate("", triin)
connectivity = triout.trianglelist
cells = [
    MeshCell(VTKCellTypes.VTK_TRIANGLE, connectivity[:, i]) for i = 1:size(connectivity)[2]
]
filename = "stress-poly-" * string(polyorder) * "-nelmts-" * string(nelmts)
vtkfile = vtk_grid(
    "examples/transformation-strain/kubo-cube/incoherent/" * filename,
    qpcoords[1, :],
    qpcoords[2, :],
    cells,
)
vtkfile["pressure"] = pressure
vtkfile["s11"] = qpstress[1, :]
vtkfile["s22"] = qpstress[2, :]
vtkfile["s12"] = qpstress[3, :]
vtkfile["s33"] = qpstress[4, :]
outfiles = vtk_save(vtkfile)


text:
# D: the dimension of the manifold
# E: the type of the singularity
abstract type Singularity end
abstract type Singularity6D <: Singularity end
abstract type Singularity5D <: Singularity end
abstract type Singularity4D <: Singularity end

struct Singularity6DPoint  <: Singularity6D T::SVector{1,Int64}; S::SVector{1,Int64} end
struct Singularity6DEdge   <: Singularity6D T::SVector{2,Int64}; S::SVector{2,Int64} end
struct Singularity6DFace   <: Singularity6D T::SVector{3,Int64}; S::SVector{3,Int64} end
struct Singularity6DVolume <: Singularity6D T::SVector{4,Int64}; S::SVector{4,Int64} end

struct Singularity5DPositiveDistance  <: Singularity5D end
struct Singularity5DPoint  <: Singularity5D T::SVector{1,Int64}; S::SVector{1,Int64} end
struct Singularity5DEdge   <: Singularity5D T::SVector{2,Int64}; S::SVector{2,Int64} end
struct Singularity5DFace   <: Singularity5D T::SVector{3,Int64}; S::SVector{3,Int64} end

struct Singularity4DPoint  <: Singularity4D T::SVector{1,Int64}; S::SVector{1,Int64} end
struct Singularity4DEdge   <: Singularity4D T::SVector{2,Int64}; S::SVector{2,Int64} end
struct Singularity4DFace   <: Singularity4D T::SVector{3,Int64}; S::SVector{3,Int64} end

#=
struct Singularity{D,E}
    T::Array{Int64}
    S::Array{Int64}
end
=#

function singularity_detection(t,s)

    sing = 0;

    D = dimension(t)+dimension(s)
    idx_t = []
    idx_s = []
    for i in 1:length(t)
        v = t[i]
        for j in 1:length(s)
            w = s[j]
            if norm(w-v) < eps(eltype(v)) * 1.0e3
                sing += 1
                push!(idx_t,i)
                push!(idx_s,j)
                break
            end
        end
    end

    if D == 4
        sing == 1 && return Singularity4DPoint(idx_t,idx_t)
        sing == 2 && return Singularity4DEdge(idx_t,idx_t)
        sing == 3 && return Singularity4DFace(idx_t,idx_t)
    elseif D == 5
        sing == 1 && return Singularity5DPoint(idx_t,idx_t)
        sing == 2 && return Singularity5DEdge(idx_t,idx_t)
        sing == 3 && return Singularity5DFace(idx_t,idx_t)
    elseif D == 6
        sing == 1 && return Singularity6DPoint(idx_t,idx_t)
        sing == 2 && return Singularity6DEdge(idx_t,idx_t)
        sing == 3 && return Singularity6DFace(idx_t,idx_t)
        sing == 4 && return Singularity6DVolume(idx_t,idx_t)
    end
    #return Singularity{D,sing}(idx_t, idx_s )
end

#=
function tetrahedron_face_order(I)
    n = length(I)
    K = zeros(Int64,4)
    for i in 1:n
        for j in 1:n
            if I[j] == i
                K[i] = j
                break
            end
        end
    end

    return SVector{4}(K)
end

function tetrahedron_edge_order(I)

    n = length(I)
    J = zeros(Int64,4)
    for i in 1:n
        for j in 1:n
            if I[j] == i
                J[i] = j
                break
            end
        end
    end

    edges = collect(combinations(J,2))
    ref_edges  = collect(combinations([1,2,3,4],2))
    n = length(edges)
    K = zeros(Int64,6)
    O = zeros(Int64,6)

    for i in 1:n
        for j in 1:n
            if edges[i] == ref_edges[j]
                K[i] = j
                O[i] = 1
                break
            elseif reverse(edges[i]) == ref_edges[j]
                K[i] = j
                O[i] = -1
                break
            end
        end
    end

    return SVector{6}(K),SVector{6}(O)
end
=#

""""
Tetrahedron-Tetrahedron reordering
"""
#Common Vertex
function reorder(sing::Singularity6DPoint)

    # Find the permutation P of t and s that make
    # Pt = [P, A1, A2, A3]
    # Ps = [P, B1, B2, B3]

    i = sing.T[1]-1

    I = circshift([1,2,3,4],-i)
    for k in 1:i
        reverse!(I,2,4)
    end

    j = sing.S[1]-1
    J = circshift([1,2,3,4],-j)
    for k in 1:j

        reverse!(J,2,4)
    end

    return SVector{4}(I), SVector{4}(J)
end

#Common Edge
function reorder(sing::Singularity6DEdge)

    # Find the permutation P of t and s that make
    # Pt = [P1, A1, A2, P2]
    # Ps = [P1, B1, B2, P2]

    I = sort(sing.T)

    I2 = setdiff([1,2,3,4],I)

    for i in 1:I2[2]-I2[1]-1
        reverse!(I2)
    end
    I3 = hcat(I,I2)

    J = sort(sing.S)

    J2 = setdiff([1,2,3,4],J)

    for i in 1:J2[2]-J2[1]-1
        reverse!(J2)
    end
    J3 = hcat(J,J2)

    return SVector{4}(I3), SVector{4}(J3)
end

#Common Face
function reorder(sing::Singularity6DFace)
    # Find the permutation P of t and s that make
    # Pt = [P1, P2, A1, P3]
    # Ps = [P1, P2, A1, P3]

    i = setdiff([1,2,3,4],sing.T)[1]-1

    I = circshift([1,2,3,4],-i)
    for k in 1:i

        reverse!(I,2,4)
    end
    I = circshift(I,2)

    j = setdiff([1,2,3,4],sing.S)[1]-1
    J = circshift([1,2,3,4],-j)
    for k in 1:j

        reverse!(J,2,4)
    end
    J = circshift(J,2)

    return SVector{4}(I), SVector{4}(J)
end

#Common Volume
function reorder(sing::Singularity6DVolume)

    I = SVector{4,Int64}([1,2,3,4])
    J = SVector{4,Int64}([1,2,3,4])

    return I, J
end

"""
Tetrahedron-Triangle reordering
"""

#Common Volume
function reorder(sing::Singularity5DPositiveDistance)

    I = SVector{4,Int64}([1,2,3,4])
    J = SVector{3,Int64}([1,2,3])

    return I, J
end

#Common Vertex
function reorder(sing::Singularity5DPoint)

    # Find the permutation P of t and s that make
    # Pt = [P, A1, A2, A3]
    # Ps = [P, B1, B2]

    i = sing.T[1]-1
    I = circshift([1,2,3,4],-i)
    for k in 1:i
        reverse!(I,2,4)
    end

    J = circshift([1,2,3],-sing.S[1]+1)

    return SVector{4}(I), SVector{3}(J)
end

#Common Edge
function reorder(sing::Singularity5DEdge)

    # Find the permutation P of t and s that make
    # Pt = [P1, A1, A2, P2]
    # Ps = [P1, B1, P2]

    #Tetrahedron
    I = sort(sing.T)

    I2 = setdiff([1,2,3,4],I)

    for i in 1:I2[2]-I2[1]-1
        reverse!(I2)
    end
    I3 = hcat(I,I2)

    #Triangle
    s = setdiff([1,2,3],sing.S)
    J = circshift([1,2,3],-s[1]+2)

    return SVector{4}(I3), SVector{3}(J)
end

#Common Face
function reorder(sing::Singularity5DFace)
    # Find the permutation P of t and s that make
    # Pt = [P1, P2, A1, P3]
    # Ps = [P1, P2, P3]

    #Tetrahedron
    i = setdiff([1,2,3,4],sing.T)[1]-1

    I = circshift([1,2,3,4],-i)
    for k in 1:i

        reverse!(I,2,4)
    end
    I = circshift(I,2)

    #Triangle
    J = SVector{3,Int64}([1,2,3])

    return SVector{4}(I), J
end

"""
Triangle-Trangle reordering
"""

#Commmon Vertex
function reorder(sing::Singularity4DPoint)

    # Find the permutation P of t and s that make
    # Pt = [P, A1, A2]
    # Ps = [P, B1, B2]

    I = circshift([1,2,3],-sing.T[1]+1)

    J = circshift([1,2,3],-sing.S[1]+1)

    return I, J
end

#Commmon Edge
function reorder(sing::Singularity4DEdge)

    # Find the permutation P of t and s that make
    # Pt = [P1, A1, P2]
    # Ps = [P1, B1, P2]

    t = setdiff([1,2,3],sing.T)
    I = circshift([1,2,3],-t[1]+2)
    s = setdiff([1,2,3],sing.S)
    J = circshift([1,2,3],-s[1]+2)

    return I, J
end

#Commmon Face
function reorder(sing::Singularity4DFace)

    # Find the permutation P of t and s that make
    # Pt = [P1, P2, P3]
    # Ps = [P1, P2, P3]

    I = [1,2,3]
    J = [1,2,3]

    return I, J
end


text:
"""
$(SIGNATURES)

Track the signal `signal` based on the current tracking `state`, the sampling frequency
`sampling_frequency` and PRN `prn`. Optional configurations are:
- Post correlation filter `post_corr_filter` defaults to `get_default_post_corr_filter(...)`
- Intermediate frequency `intermediate_frequency` defaults to 0Hz
- Maximal integration time `max_integration_time` defaults to 1ms. The actual integration
  time can be less in order to find the secondary code or the bit shift. Once this is found
  the integration time is identical to the maximal integration time.
- Minimal integration time `min_integration_time` defaults to 0.75ms. It's the minimal
  integration time, that leads to a valid correlation result. It is only used for the first
  integration period.
- Sample shift of the correlator replica `correlator_sample_shifts` defaults to
  `get_correlator_sample_shifts(...)`
- Bandwidth of the carrier loop `carrier_loop_filter_bandwidth` defaults to 18Hz
- Bandwidth of the code loop `code_loop_filter_bandwidth` defaults to 1Hz
- Velocity aiding `velocity_aiding` defaults to 0Hz
"""
function track(
        signal,
        state::TrackingState{S, C, CALF, COLF, CN, DS, CAR, COR},
        sampling_frequency;
        post_corr_filter = get_default_post_corr_filter(get_correlator(state)),
        intermediate_frequency = 0.0Hz,
        max_integration_time::typeof(1ms) = 1ms,
        min_integration_time::typeof(1.0ms) = 0.75ms,
        correlator_sample_shifts = get_correlator_sample_shifts(get_system(state),
            get_correlator(state), sampling_frequency, 0.5),
        early_late_index_shift = get_early_late_index_shift(get_system(state),
            correlator_sample_shifts, get_correlator(state), sampling_frequency, 0.5),
        carrier_loop_filter_bandwidth = 18Hz,
        code_loop_filter_bandwidth = 1Hz,
        velocity_aiding = 0Hz,
) where {
    S <: AbstractGNSS,
    C <: AbstractCorrelator,
    CALF <: AbstractLoopFilter,
    COLF <: AbstractLoopFilter,
    CN <: AbstractCN0Estimator,
    DS,
    CAR,
    COR
}
    prn = get_prn(state)
    system = get_system(state)
    if get_data_frequency(system) != 0Hz
        @assert rem(1 / get_data_frequency(system), max_integration_time) == 0ms
    end
    correlator = get_correlator(state)
    num_ants = get_num_ants(correlator)
    size(signal, 2) == num_ants || throw(ArgumentError("The second dimension of the signal should be equal to the number of antennas specified by num_ants = NumAnts(N) in the TrackingState."))
    if typeof(system.codes) <: CuMatrix
        typeof(signal) <: StructArray || throw(ArgumentError("Signal is not a StructArray, initialize the signal properly and try again."))
        typeof(signal.re) <: CuArray && typeof(state.system.codes) <: CuArray || throw(ArgumentError("Signal and GNSS codes are not of the same type. Please check if CPU or GPU is used."))
    end
    downconverted_signal_temp = get_downconverted_signal(state)
    downconverted_signal = resize!(downconverted_signal_temp, size(signal, 1), signal)
    carrier_replica = get_carrier(state)
    resize!(choose(carrier_replica, signal), size(signal, 1))
    code_replica = get_code(state)
    resize!(code_replica, size(signal, 1) + correlator_sample_shifts[end]-correlator_sample_shifts[1])
    init_carrier_doppler = get_init_carrier_doppler(state)
    init_code_doppler = get_init_code_doppler(state)
    carrier_doppler = get_carrier_doppler(state)
    code_doppler = get_code_doppler(state)
    carrier_phase = get_carrier_phase(state)
    code_phase = get_code_phase(state)
    sc_bit_detector = get_sc_bit_detector(state)
    carrier_loop_filter = get_carrier_loop_filter(state)
    code_loop_filter = get_code_loop_filter(state)
    prompt_accumulator = get_prompt_accumulator(state)
    integrated_samples = get_integrated_samples(state)
    cn0_estimator = get_cn0_estimator(state)
    signal_start_sample = 1
    bit_buffer = BitBuffer()
    valid_correlator = zero(correlator)
    valid_correlator_carrier_phase = 0.0
    valid_correlator_carrier_frequency = 0.0Hz
    got_correlator = false
    while true
        num_samples_left_to_integrate = get_num_samples_left_to_integrate(
            system,
            max_integration_time,
            sampling_frequency,
            code_doppler,
            code_phase,
            found(sc_bit_detector)
        )
        signal_samples_left = get_num_samples(signal) - signal_start_sample + 1
        num_samples_left = min(num_samples_left_to_integrate, signal_samples_left)
        carrier_frequency = get_current_carrier_frequency(
            intermediate_frequency,
            carrier_doppler
        )
        code_frequency = get_current_code_frequency(system, code_doppler)
        correlator = downconvert_and_correlate!(
            system,
            signal,
            correlator,
            code_replica,
            code_phase,
            carrier_replica,
            carrier_phase,
            downconverted_signal,
            code_frequency,
            correlator_sample_shifts,
            carrier_frequency,
            sampling_frequency,
            signal_start_sample,
            num_samples_left,
            prn
        )
        integrated_samples += num_samples_left
        carrier_phase = update_carrier_phase(
            num_samples_left,
            carrier_frequency,
            sampling_frequency,
            carrier_phase,
        )
        prev_code_phase = code_phase
        code_phase = update_code_phase(
            system,
            num_samples_left,
            code_frequency,
            sampling_frequency,
            code_phase,
            found(sc_bit_detector)
        )
        integration_time = integrated_samples / sampling_frequency
        if num_samples_left == num_samples_left_to_integrate &&
                integration_time >= min_integration_time
            got_correlator = true
            correlator = normalize(correlator, integrated_samples)
            valid_correlator = correlator
            valid_correlator_carrier_phase = carrier_phase
            valid_correlator_carrier_frequency = carrier_frequency
            filtered_correlator = filter(post_corr_filter, correlator)
            pll_discriminator = pll_disc(
                system,
                filtered_correlator,
                correlator_sample_shifts
            )
            dll_discriminator = dll_disc(
                system,
                filtered_correlator,
                correlator_sample_shifts,
                early_late_index_shift,
                code_frequency / sampling_frequency
            )
            carrier_freq_update, carrier_loop_filter = filter_loop(
                carrier_loop_filter,
                pll_discriminator,
                integration_time,
                carrier_loop_filter_bandwidth
            )
            code_freq_update, code_loop_filter = filter_loop(
                code_loop_filter,
                dll_discriminator,
                integration_time,
                code_loop_filter_bandwidth
            )
            carrier_doppler, code_doppler = aid_dopplers(
                system,
                init_carrier_doppler,
                init_code_doppler,
                carrier_freq_update,
                code_freq_update,
                velocity_aiding
            )
            cn0_estimator = update(
                cn0_estimator,
                get_prompt(filtered_correlator, correlator_sample_shifts)
            )
            bit_buffer, prompt_accumulator = buffer(
                system,
                bit_buffer,
                prompt_accumulator,
                found(sc_bit_detector),
                prev_code_phase,
                code_phase,
                max_integration_time,
                get_prompt(filtered_correlator, correlator_sample_shifts)
            )
            sc_bit_detector = find(
                system,
                sc_bit_detector,
                get_prompt(filtered_correlator, correlator_sample_shifts)
            )
            correlator = zero(correlator)
            integrated_samples = 0
        end

        num_samples_left == signal_samples_left && break
        signal_start_sample += num_samples_left
    end
    next_state = TrackingState{S, C, CALF, COLF, CN, DS, CAR, COR}(
        prn,
        system,
        init_carrier_doppler,
        init_code_doppler,
        carrier_doppler,
        code_doppler,
        carrier_phase,
        code_phase,
        correlator,
        carrier_loop_filter,
        code_loop_filter,
        sc_bit_detector,
        integrated_samples,
        prompt_accumulator,
        cn0_estimator,
        downconverted_signal,
        carrier_replica,
        code_replica
    )
    estimated_cn0 = estimate_cn0(cn0_estimator, max_integration_time)
    TrackingResults(
        next_state,
        valid_correlator,
        correlator_sample_shifts,
        early_late_index_shift,
        valid_correlator_carrier_frequency,
        valid_correlator_carrier_phase,
        got_correlator,
        bit_buffer,
        estimated_cn0
    )
end

function downconvert_and_correlate!(
    system,
    signal,
    correlator,
    code_replica,
    code_phase,
    carrier_replica,
    carrier_phase,
    downconverted_signal,
    code_frequency,
    correlator_sample_shifts,
    carrier_frequency,
    sampling_frequency,
    signal_start_sample,
    num_samples_left,
    prn
)
    gen_code_replica!(
        code_replica,
        system,
        code_frequency,
        sampling_frequency,
        code_phase,
        signal_start_sample,
        num_samples_left,
        correlator_sample_shifts,
        prn
    )
    gen_carrier_replica!(
        choose(carrier_replica, signal),
        carrier_frequency,
        sampling_frequency,
        carrier_phase,
        signal_start_sample,
        num_samples_left
    )
    downconvert!(
        choose(downconverted_signal, signal),
        signal,
        choose(carrier_replica, signal),
        signal_start_sample,
        num_samples_left
    )
    correlate(
        correlator,
        choose(downconverted_signal, signal),
        code_replica,
        correlator_sample_shifts,
        signal_start_sample,
        num_samples_left
    )
end

# CUDA downconvert_and_correlate for num_ants > 1
function downconvert_and_correlate!(
    system::AbstractGNSS{C},
    signal::AbstractMatrix,
    correlator::T,
    code_replica,
    code_phase,
    carrier_replica,
    carrier_phase,
    downconverted_signal,
    code_frequency,
    correlator_sample_shifts,
    carrier_frequency,
    sampling_frequency,
    signal_start_sample,
    num_samples_left,
    prn
) where {C <: CuMatrix, T <: AbstractCorrelator}
    accumulator_result = downconvert_and_correlate_kernel_wrapper(
        system,
        view(signal, signal_start_sample:signal_start_sample - 1 + num_samples_left,:),
        correlator,
        code_phase,
        carrier_phase,
        code_frequency,
        correlator_sample_shifts,
        carrier_frequency,
        sampling_frequency,
        signal_start_sample,
        num_samples_left,
        prn
    )
    return T(map(+, get_accumulators(correlator), eachcol(Array(accumulator_result[1,:,:]))))
end

# CUDA downconvert_and_correlate for num_ants = 1
function downconvert_and_correlate!(
    system::AbstractGNSS{C},
    signal::AbstractVector,
    correlator::T,
    code_replica,
    code_phase,
    carrier_replica,
    carrier_phase,
    downconverted_signal,
    code_frequency,
    correlator_sample_shifts,
    carrier_frequency,
    sampling_frequency,
    signal_start_sample,
    num_samples_left,
    prn
) where {C <: CuMatrix, T <: AbstractCorrelator}
    accumulator_result = downconvert_and_correlate_kernel_wrapper(
        system,
        view(signal, signal_start_sample:signal_start_sample - 1 + num_samples_left),
        correlator,
        code_phase,
        carrier_phase,
        code_frequency,
        correlator_sample_shifts,
        carrier_frequency,
        sampling_frequency,
        signal_start_sample,
        num_samples_left,
        prn
    )
    addition(a,b) = a + first(b)
    return T(map(addition, get_accumulators(correlator), eachcol(Array(accumulator_result[1,:,:]))))
end

function choose(replica::CarrierReplicaCPU, signal::AbstractArray{Complex{Float64}})
    replica.carrier_f64
end
function choose(replica::CarrierReplicaCPU, signal::AbstractArray{Complex{T}}) where T <: Number
    replica.carrier_f32
end
function choose(replica::DownconvertedSignalCPU, signal::AbstractArray{Complex{Float64}})
    replica.downconverted_signal_f64
end
function choose(replica::DownconvertedSignalCPU, signal::AbstractArray{Complex{T}}) where T <: Number
    replica.downconverted_signal_f32
end
function choose(replica::Nothing, signal::AbstractArray)
    nothing
end

"""
$(SIGNATURES)

Default post corr filter
"""
function get_default_post_corr_filter(correlator::AbstractCorrelator{T}) where T <: SVector
    x -> x[1]
end

function get_default_post_corr_filter(correlator)
    x -> x
end

"""
$(SIGNATURES)

Returns the appropiate integration time. It will be the maximum integration time once the
secondary code or the bit shift has been found.
"""
function get_integration_time(
    system::AbstractGNSS,
    max_integration_time,
    secondary_code_or_bit_found::Bool
)
    ifelse(
        secondary_code_or_bit_found,
        max_integration_time,
        min(
            ceil(typeof(1ms), get_code_length(system) / get_code_frequency(system)),
            max_integration_time
        )
    )
end

"""
$(SIGNATURES)

Calculates the number of chips to integrate.
"""
function get_num_chips_to_integrate(
    system::AbstractGNSS,
    max_integration_time,
    current_code_phase,
    secondary_code_or_bit_found
)
    max_phase = Int(upreferred(get_code_frequency(system) *
        get_integration_time(system, max_integration_time, secondary_code_or_bit_found)))
    current_phase_mod_max_phase = mod(current_code_phase, max_phase)
    return max_phase - current_phase_mod_max_phase
end

"""
$(SIGNATURES)

Calculates the number of samples to integrate.
"""
function get_num_samples_left_to_integrate(
    system::AbstractGNSS,
    max_integration_time,
    sampling_frequency,
    current_code_doppler,
    current_code_phase,
    secondary_code_or_bit_found
)
    phase_to_integrate = get_num_chips_to_integrate(
        system,
        max_integration_time,
        current_code_phase,
        secondary_code_or_bit_found
    )
    code_frequency = get_code_frequency(system) + current_code_doppler
    ceil(Int, phase_to_integrate * sampling_frequency / code_frequency)
end

"""
$(SIGNATURES)

Aid dopplers. That is velocity aiding for the carrier doppler and carrier aiding
for the code doppler.
"""
function aid_dopplers(
    system::AbstractGNSS,
    init_carrier_doppler,
    init_code_doppler,
    carrier_freq_update,
    code_freq_update,
    velocity_aiding
)
    carrier_doppler = carrier_freq_update + velocity_aiding
    code_doppler = code_freq_update + carrier_doppler * get_code_center_frequency_ratio(system)
    init_carrier_doppler + carrier_doppler, init_code_doppler + code_doppler
end

"""
$(SIGNATURES)

Get the number of samples in the signal.
"""
@inline function get_num_samples(signal)
    length(signal)
end

@inline function get_num_samples(signal::AbstractMatrix)
    size(signal, 1)
end

function resize!(ds::DownconvertedSignalCPU, b::Integer, signal::AbstractVector)
    resize!(choose(ds, signal), b)
    ds
end

function resize!(ds::DownconvertedSignalCPU, b::Integer, signal::AbstractMatrix{Complex{Float64}})
    num_ants = size(signal, 2)
    DownconvertedSignalCPU(
        ds.downconverted_signal_f32,
        size(ds.downconverted_signal_f64, 1) == b ?
            ds.downconverted_signal_f64 :
            StructArray{Complex{Float64}}((Matrix{Float64}(undef, b, num_ants), Matrix{Float64}(undef, b, num_ants)))
    )
end

function resize!(ds::DownconvertedSignalCPU, b::Integer, signal::AbstractMatrix{Complex{T}}) where T <: Number
    num_ants = size(signal, 2)
    DownconvertedSignalCPU(
        size(ds.downconverted_signal_f32, 1) == b ?
            ds.downconverted_signal_f32 :
            StructArray{Complex{Float32}}((Matrix{Float32}(undef, b, num_ants), Matrix{Float32}(undef, b, num_ants))),
        ds.downconverted_signal_f64
    )
end

# No need for resizing when dealing with GPU signals
function resize!(ds::Nothing, b::Integer, signal::AbstractArray)
    return ds
end
# No need for resizing the GPU GNSS codes
function resize!(codes::Nothing, b::Integer)
    return codes
end


text:
"""
    update_wavelets!(model::AbstractModel)

Compute wavelet transform of velocities to define the coarse grid used in multigrid preconditioner.
"""
function update_wavelets!(model::AbstractModel)
    @unpack wu,wv,gu,gv=model.fields
    @unpack params,solver_params=model

    wu.wavelets[:] .= gu.dwt*(gu.crop*gu.u[:])
    wv.wavelets[:] .= gv.dwt*(gv.crop*gv.v[:])

    wu.mask .= (abs.(wu.wavelets) .>= solver_params.wavelet_threshold)
    wv.mask .= (abs.(wv.wavelets) .>= solver_params.wavelet_threshold)

    wu.n[] = count(wu.mask)
    wv.n[] = count(wv.mask)

    wu.crop[] .= Diagonal(float(wu.mask[:]))
    wv.crop[] .= Diagonal(float(wv.mask[:]))

    wu.samp[] = wu.crop[][wu.mask[:],:]
    wv.samp[] = wv.crop[][wv.mask[:],:]

    wu.spread[] = sparse(wu.samp[]')
    wv.spread[] = sparse(wv.samp[]')

    return model
end
"""
    wavelet_matrix(n,levels,direction)

Compute matrix used to apply wavelet transform in one dimension.
"""
function wavelet_matrix(n,levels,direction);

    identity=spdiagm(n,n, 0 => ones(n))

    isforward = lowercase(direction) == "forward"
    isreverse = lowercase(direction) == "reverse"

    if isforward
        fwav=identity
    elseif isreverse
        rwav=identity
    else
        error("Direction of wavelet transform not defined properly")
    end

    step=1

    for j=levels:-1:0

        step=step*2

        ixodd=1:step:(n-div(step,2))
        ixeven=ixodd.+div(step,2)
        ixothers=setdiff(1:n, union(ixodd,ixeven))

        odds=identity[ixodd,:]
        evens=identity[ixeven,:]
        others=identity[ixothers,:]

        if isforward
            #Forward Haar wavelet lifting scheme
            fwav = (others'*others+0.5*(evens'*(evens - odds) + odds'*(odds + evens)))*fwav
        elseif isreverse
            #Reverse Haar wavelet lifting scheme
            rwav = rwav*(others'*others + evens'*(evens + odds) + odds'*(odds - evens))
        else
            error("Direction of wavelet transform not defined properly")
        end

    end

    if isforward
            #Forward Haar wavelet lifting scheme
            return fwav
    elseif isreverse
            #Reverse Haar wavelet lifting scheme
            return rwav
    else
            error("Direction of wavelet transform not defined properly")
    end

end