text:
function constraint_frequency_stab_OPF_sensitivity(pm::_PM.AbstractPowerModel) #gurobi
    if pm.setting["Permanentloss"] == true
        constraint_frequency_stab_OPF_MP_PL_sensitivity(pm)
    elseif pm.setting["FSprotection"] == true || pm.setting["NSprotection"] == true
        constraint_frequency_stab_OPF_MP_FSNS_sensitivity(pm)
    end
end

function constraint_frequency_stab_OPF_MP_PL_sensitivity(pm::_PM.AbstractPowerModel)
    base_nws = pm.setting["base_list"]
    cont_nws = pm.setting["Cont_list"]
    ev_syncarea = pm.setting["syncarea"]

    for (base, nw, br) in cont_nws
        load = _PM.ref(pm, nw, :load, 1) # get index later on
        syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
        Zb1 = _PM.var(pm, nw)[:zb1]

        # for i in syncarea # when there is more than one synchornous areas
        i = ev_syncarea
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
        JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"])
        JuMP.@constraint(pm.model,  Phvdccaux[i]  == 0)

        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]
        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        H = reserves[bi_bp[i]]["H"]
        Tcl = reserves[bi_bp[i]]["Tcl"]
        M = 1 # largest possible Pl=0.15 for now

        z1 = _PM.var(pm, nw, :z1, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))

        ############## frequency minimum in interval 1 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

          k11 = _PM.var(pm, nw, :k11, i)
          k12 = _PM.var(pm, nw, :k12, i)
          k11_dup = _PM.var(pm, nw, :k11_dup, i)
          k12_dup = _PM.var(pm, nw, :k12_dup, i)
          variable_on_off_switch(pm, k11, k11_dup, M, z1)
          variable_on_off_switch(pm, k12, k12_dup, M, z1)

          JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
          JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
          JuMP.@constraint(pm.model,k12_dup^2 <=  k11_dup)

         ############## frequency minimum in interval 2 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  M*(1-z41))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg   - M*(z41))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tf)   <=   Phvdcoaux[i]  + M*(1-z42))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tf)   >= e +  Phvdcoaux[i]  - M*(z42))
         JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
         JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

         k41 = _PM.var(pm, nw, :k41, i)
         k42 = _PM.var(pm, nw, :k42, i)
         k43 = _PM.var(pm, nw, :k43, i)
         k41_dup = _PM.var(pm, nw, :k41_dup, i)
         k42_dup = _PM.var(pm, nw, :k42_dup, i)
         k43_dup = _PM.var(pm, nw, :k43_dup, i)
         variable_on_off_switch(pm, k41, k41_dup, M, z4)
         variable_on_off_switch(pm, k42, k42_dup, M, z4)
         variable_on_off_switch(pm, k43, k43_dup, M, z4)

         JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4) )
         JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
         JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf )/2 )
         JuMP.@constraint(pm.model,[k41_dup/sqrt(2), k42_dup/sqrt(2), k43_dup, 0] in JuMP.RotatedSecondOrderCone() )

         JuMP.@constraint(pm.model, z1 + z4 == 1)
     end
end

function constraint_frequency_stab_OPF_MP_FSNS_sensitivity(pm::_PM.AbstractPowerModel)

    base_nws = pm.setting["base_list"]
    cont_nws = pm.setting["Cont_list"]
    ev_syncarea = pm.setting["syncarea"]

    for (base, nw, br) in cont_nws
        load = _PM.ref(pm, nw, :load, 1)
        syncarea  = _PM.ids(pm, nw, :arcs_reserves_syn)
        Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux]
        Phvdccaux = _PM.var(pm, nw)[:Phvdccaux]
        Pconv1 = _PM.var(pm, base, :pconv_tf_fr)
        Pconv2 = _PM.var(pm, nw, :pconv_tf_fr)
        Zb1 = _PM.var(pm, nw)[:zb1]
        # for i in syncarea
        i = ev_syncarea
        conv_conn = _PM.ref(pm, nw, :bus_arcs_conv, i)

        if pm.setting["FSprotection"] == true
            # display("FSprotection")
            Mmax = 72/ load["pd"]
            display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_convs_dirct, i)) / load["pd"]) )
            display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - 0.5*sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_convs_dirct, i)) /load["pd"]) )
            # display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"] + sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"] ) )
         elseif pm.setting["NSprotection"] == true
            # display("NSprotection")
            # display(load["pd"])
            display(JuMP.@constraint(pm.model,  Phvdcoaux[i]  == - sum(Pconv1[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) /load["pd"]))
            display(JuMP.@constraint(pm.model,  Phvdccaux[i]  == - 0.5*sum(Pconv2[a] for a in _PM.ref(pm, nw, :bus_arcs_conv, i)) / load["pd"]))
        end

        JuMP.set_upper_bound(Phvdcoaux[i], 72/load["pd"])
        JuMP.set_upper_bound(Phvdccaux[i], 72/load["pd"])
        Pg = _PM.var(pm, nw, :Pgg, i)/load["pd"]
        Pf = _PM.var(pm, nw, :Pff, i)/load["pd"]

        reserves = _PM.ref(pm, nw, :reserves)
        bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
        Td = reserves[bi_bp[i]]["Td"]
        Tg= reserves[bi_bp[i]]["Tg"]
        Tf = reserves[bi_bp[i]]["Tf"]
        Tcl = reserves[bi_bp[i]]["Tcl"]
        H = reserves[bi_bp[i]]["H"]

        z1 = _PM.var(pm, nw, :z1, i)
        z2 = _PM.var(pm, nw, :z2, i)
        z3 = _PM.var(pm, nw, :z3, i)
        z4 = _PM.var(pm, nw, :z4, i)
        z11 = _PM.var(pm, nw, :z11, i)
        z21 = _PM.var(pm, nw, :z21, i)
        z31 = _PM.var(pm, nw, :z31, i)
        z41 = _PM.var(pm, nw, :z41, i)
        z12 = _PM.var(pm, nw, :z12, i)
        z22 = _PM.var(pm, nw, :z22, i)
        z32 = _PM.var(pm, nw, :z32, i)
        z42 = _PM.var(pm, nw, :z42, i)
        e = 1e-06
        p2r = (50/(2*H))
        M = 1
        # M = ub(Pf + Pg +  Phvdccaux[i])
       ############## frequency minimum in interval 1 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e <=   Pf + (Pg/Tg)*Tf + M*(1-z11))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tf - M*(z11))
         JuMP.@constraint(pm.model, 0 <=   Phvdcoaux[i]  + M*(1-z12))
         JuMP.@constraint(pm.model,  0  >= e +  Phvdcoaux[i]  - M*(z12))
         JuMP.@constraint(pm.model,  0 <= z11 + z12 - 2*z1)
         JuMP.@constraint(pm.model,  z11 + z12 - 2*z1 <= 1)

         k11 = _PM.var(pm, nw, :k11, i)
         k12 = _PM.var(pm, nw, :k12, i)
         k11_dup = _PM.var(pm, nw, :k11_dup, i)
         k12_dup = _PM.var(pm, nw, :k12_dup, i)
         variable_on_off_switch(pm, k11, k11_dup, M, z1)
         variable_on_off_switch(pm, k12, k12_dup, M, z1)

         JuMP.@constraint(pm.model,k11 ==  (H/50)*(Pg/Tg + Pf/Tf) )
         JuMP.@constraint(pm.model,k12 ==  (Phvdcoaux[i])/2)
         JuMP.@constraint(pm.model,k12_dup^2 <=  k11_dup)

     ############## frequency minimum in interval 2 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e  <=  Pf + (Pg/Tg)*Tcl + M*(1-z21))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*Tcl - M*(z21))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tf <=   Phvdcoaux[i]  + M*(1-z22))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tf  >= e +  Phvdcoaux[i]  - M*(z22))
         JuMP.@constraint(pm.model,  0 <= z21 + z22 - 2*z2)
         JuMP.@constraint(pm.model,  z21 + z22 - 2*z2 <= 1)

         k21 = _PM.var(pm, nw, :k21, i)
         k22 = _PM.var(pm, nw, :k22, i)
         k23 = _PM.var(pm, nw, :k23, i)

         k21_dup = _PM.var(pm, nw, :k21_dup, i)
         k22_dup = _PM.var(pm, nw, :k22_dup, i)
         k23_dup = _PM.var(pm, nw, :k23_dup, i)
         variable_on_off_switch(pm, k21, k21_dup, M, z2)
         variable_on_off_switch(pm, k22, k22_dup, M, z2)
         variable_on_off_switch(pm, k23, k23_dup, M, z2)

         JuMP.@constraint(pm.model,k21 ==  (H/50 - Pf*Tf/4))
         JuMP.@constraint(pm.model,k22 ==  Pg/Tg)
         JuMP.@constraint(pm.model,k23 ==  (Phvdcoaux[i] - Pf)/2)

         JuMP.@constraint(pm.model,[k21_dup/sqrt(2), k22_dup/sqrt(2), k23_dup, 0] in JuMP.RotatedSecondOrderCone() )

         ############## frequency minimum in interval 3 ##########################################
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  +e  <=   Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  + M*(1-z31))
         JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  - M*(z31))
         JuMP.@constraint(pm.model, Pf + (Pg/Tg)*Tcl <=   Phvdcoaux[i]  + M*(1-z32))
         JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*Tcl  >= e+  Phvdcoaux[i]  - M*(z32))
         JuMP.@constraint(pm.model,  0 <= z31 + z32 - 2*z3)
         JuMP.@constraint(pm.model,  z31 + z32 - 2*z3 <= 1)

         k31 = _PM.var(pm, nw, :k31, i)
         k32 = _PM.var(pm, nw, :k32, i)
         k33 = _PM.var(pm, nw, :k33, i)
         k31_dup = _PM.var(pm, nw, :k31_dup, i)
         k32_dup = _PM.var(pm, nw, :k32_dup, i)
         k33_dup = _PM.var(pm, nw, :k33_dup, i)
         variable_on_off_switch(pm, k31, k31_dup, M, z3)
         variable_on_off_switch(pm, k32, k32_dup, M, z3)
         variable_on_off_switch(pm, k33, k33_dup, M, z3)

         JuMP.@constraint(pm.model, k31 ==  (H/50 - Pf*Tf/4 + Phvdccaux[i]*Tcl^2/(Td*4) ) )
         JuMP.@constraint(pm.model, k32 ==  (Pg/Tg + Phvdccaux[i]/Td) )
         JuMP.@constraint(pm.model, k33 ==  (Phvdcoaux[i] - Pf + Phvdccaux[i]*Tcl/Td)/2 )

         display(JuMP.@constraint(pm.model,[k31_dup/sqrt(2), k32_dup/sqrt(2), k33_dup, 0] in JuMP.RotatedSecondOrderCone() ))

        ############## frequency minimum in interval 4 ##########################################
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  + e <=   Pf + Pg +  Phvdccaux[i]  + M*(1-z41))
        JuMP.@constraint(pm.model,  Phvdcoaux[i]  >= Pf + Pg +  Phvdccaux[i]  - M*(z41))
        JuMP.@constraint(pm.model, Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]  <=   Phvdcoaux[i]  + M*(1-z42))
        JuMP.@constraint(pm.model,  Pf + (Pg/Tg)*(Tcl+Td) +  Phvdccaux[i]   >= e +  Phvdcoaux[i]  - M*(z42))
        JuMP.@constraint(pm.model,  0 <= z41 + z42 - 2*z4)
        JuMP.@constraint(pm.model,  z41 + z42 - 2*z4 <= 1)

        k41 = _PM.var(pm, nw, :k41, i)
        k42 = _PM.var(pm, nw, :k42, i)
        k43 = _PM.var(pm, nw, :k43, i)
        k41_dup = _PM.var(pm, nw, :k41_dup, i)
        k42_dup = _PM.var(pm, nw, :k42_dup, i)
        k43_dup = _PM.var(pm, nw, :k43_dup, i)
        variable_on_off_switch(pm, k41, k41_dup, M, z4)
        variable_on_off_switch(pm, k42, k42_dup, M, z4)
        variable_on_off_switch(pm, k43, k43_dup, M, z4)

        JuMP.@constraint(pm.model, k41 ==  (H/50 - Pf*Tf/4 - Phvdccaux[i]*(Td + 2*Tcl)/4) )
        JuMP.@constraint(pm.model, k42 ==  (Pg/Tg) )
        JuMP.@constraint(pm.model, k43 ==  (Phvdcoaux[i] - Pf - Phvdccaux[i])/2 )
        JuMP.@constraint(pm.model,[k41_dup/sqrt(2), k42_dup/sqrt(2), k43_dup, 0] in JuMP.RotatedSecondOrderCone() )

        JuMP.@constraint(pm.model,  z1+z2+z3+z4 == 1)

     end
end

function variable_frequency_stab_OPF(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    t1 = _PM.var(pm, nw)[:t1] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_t1",  lower_bound = 0,  upper_bound = 10, start = 0)
    t2 = _PM.var(pm, nw)[:t2] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_t2",  lower_bound = 0,  upper_bound = 10, start = 0)
    t3 = _PM.var(pm, nw)[:t3] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_t3",  lower_bound = 0,  upper_bound = 10, start = 0)
    t4 = _PM.var(pm, nw)[:t4] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_t4",  lower_bound = 0,  upper_bound = 10, start = 0)

    # f1 = _PM.var(pm)[:f1] = JuMP.@variable(pm.model, base_name="df1",  lower_bound = -10,  upper_bound = 10, start = 0)
    # f2 = _PM.var(pm)[:f2] = JuMP.@variable(pm.model, base_name="df2",  lower_bound = -10,  upper_bound = 10, start = 0)
    # f3 = _PM.var(pm)[:f3] = JuMP.@variable(pm.model, base_name="df3",  lower_bound = -10,  upper_bound = 10, start = 0)
    # f4 = _PM.var(pm)[:f4] = JuMP.@variable(pm.model, base_name="df4",  lower_bound = -10,  upper_bound = 10, start = 0)

    f1 = _PM.var(pm, nw)[:f1] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_df1",  lower_bound = -10,  upper_bound = 10, start = 0)
    f2 = _PM.var(pm, nw)[:f2] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_df2",  lower_bound = -10,  upper_bound = 10, start = 0)
    f3 = _PM.var(pm, nw)[:f3] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_df3",  lower_bound = -10,  upper_bound = 10, start = 0)
    f4 = _PM.var(pm, nw)[:f4] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_df4",  lower_bound = -10,  upper_bound = 10, start = 0)

    report && _IM.sol_component_value(pm, nw, :reserves, :df, _PM.ids(pm, nw, :arcs_reserves_syn), df)
    report && _IM.sol_component_value(pm, nw, :reserves, :dt, _PM.ids(pm, nw, :arcs_reserves_syn), dt)

    report && _IM.sol_component_value(pm, nw, :reserves, :f1, _PM.ids(pm, nw, :arcs_reserves_syn), f1)
    report && _IM.sol_component_value(pm, nw, :reserves, :f2, _PM.ids(pm, nw, :arcs_reserves_syn), f2)
    report && _IM.sol_component_value(pm, nw, :reserves, :f3, _PM.ids(pm, nw, :arcs_reserves_syn), f3)
    report && _IM.sol_component_value(pm, nw, :reserves, :f4, _PM.ids(pm, nw, :arcs_reserves_syn), f4)
    report && _IM.sol_component_value(pm, nw, :reserves, :t1, _PM.ids(pm, nw, :arcs_reserves_syn), t1)
    report && _IM.sol_component_value(pm, nw, :reserves, :t2, _PM.ids(pm, nw, :arcs_reserves_syn), t2)
    report && _IM.sol_component_value(pm, nw, :reserves, :t3, _PM.ids(pm, nw, :arcs_reserves_syn), t3)
    report && _IM.sol_component_value(pm, nw, :reserves, :t4, _PM.ids(pm, nw, :arcs_reserves_syn), t4)
end

function variable_on_off_switch(pm, k, k_dup, M , z)
    M = 100 # M=1 was binding for nonactivated K
    JuMP.@constraint(pm.model,  k_dup  >=  k - M*(1-z) )
    JuMP.@constraint(pm.model,  k_dup  <=  k + M*(1-z) )
    JuMP.@constraint(pm.model,  k_dup  >=  -M*z )
    JuMP.@constraint(pm.model,  k_dup  <=  M*z )
end


text:
include("poweriteration.jl")

# This code generates a dictionnary where index -> page_name

function get_index_to_page(filename)
    io = open(filename)
    lines = readlines(io)
    dictionnary = Dict()
    println("Construction index to page dictionnary")
    for line in ProgressBar(lines)
        if !startswith(line, "#") && line != ""
            index = ""
            page = ""
            is_index = true
            for element in line
                if element == '\t'
                    is_index = false
                else
                    if is_index
                        index = string(index, element)
                    else
                        page = string(page, element)
                    end
                end
            end
            dictionnary[index] = page
        end
    end

    return dictionnary
end

ext: jl

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
# ---
# title: 1415. The k-th Lexicographical String of All Happy Strings of Length n
# id: problem1415
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/the-k-th-lexicographical-string-of-all-happy-strings-of-length-n/description/>
# hidden: true
# ---
#
# A **happy string** is a string that:
#
#   * consists only of letters of the set `['a', 'b', 'c']`.
#   * `s[i] != s[i + 1]` for all values of `i` from `1` to `s.length - 1` (string is 1-indexed).
#
# For example, strings **" abc", "ac", "b"** and **" abcbabcbcb"** are all happy
# strings and strings **" aa", "baa"** and  **" ababbc"** are not happy strings.
#
# Given two integers `n` and `k`, consider a list of all happy strings of length
# `n` sorted in lexicographical order.
#
# Return _the kth string_ of this list or return an **empty string**  if there
# are less than `k` happy strings of length `n`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1, k = 3
#     Output: "c"
#     Explanation: The list ["a", "b", "c"] contains all happy strings of length 1. The third string is "c".
#
#
# **Example 2:**
#
#
#
#     Input: n = 1, k = 4
#     Output: ""
#     Explanation: There are only 3 happy strings of length 1.
#
#
# **Example 3:**
#
#
#
#     Input: n = 3, k = 9
#     Output: "cab"
#     Explanation: There are 12 different happy string of length 3 ["aba", "abc", "aca", "acb", "bab", "bac", "bca", "bcb", "cab", "cac", "cba", "cbc"]. You will find the 9th string = "cab"
#
#
# **Example 4:**
#
#
#
#     Input: n = 2, k = 7
#     Output: ""
#
#
# **Example 5:**
#
#
#
#     Input: n = 10, k = 100
#     Output: "abacbabacb"
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10`
#   * `1 <= k <= 100`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
"""
ctmc type
"""

mutable struct CTMCModel{Tv}
    Q::AbstractMatrix{Tv}
    initv::Vector{Tv}
    reward::Dict{Symbol,Vector{Tv}}
    states::Vector{Symbol}
end

function Base.show(io::IO, x::CTMCModel{Tv}) where Tv
    Base.show(io, "CTMCModel($(objectid(x)))")
end

function ctmc(Q::AbstractMatrix{Tv}, initv::Vector{Tv}, r::Dict{Symbol,Vector{Tv}}, states::Vector{Symbol}) where Tv
    CTMCModel{Tv}(Q, initv, r, states)
end

function ctmc(m::Markov, modeltype=:SparseCTMC)
    Q, initv, r, states = generate(m, modeltype=modeltype)
    ctmc(Q, initv, r, states)
end

# """
# exrss
# """

# function exrss(m::CTMCModel{Tv}; reward, maxiter=5000, steps=20, rtol=1.0e-6) where {Tv}
#     _exrss(m.Q, m.reward[reward], maxiter, steps, rtol)
# end

# function _exrss(Q::Matrix{Tv}, r::Vector{Tv}, maxiter, steps, rtol) where Tv
#     dot(gth(Q), r)
# end

# function _exrss(Q::SparseCSC{Tv,Ti}, r::Vector{Tv}, maxiter, steps, rtol) where {Tv<:Number,Ti}
#     x, conv, iter, rerror = stgs(Q, maxiter=maxiter, steps=steps, rtol=rtol)
#     if conv == false
#         @warn "GS did not converge", iter, rerror
#     end
#     dot(x, r)
# end

# function _exrss(Q::SparseMatrixCSC{Tv,Ti}, r::Vector{Tv}, maxiter, steps, rtol) where {Tv<:Number,Ti}
#     x, conv, iter, rerror = stgs(Q, maxiter=maxiter, steps=steps, rtol=rtol)
#     if conv == false
#         @warn "GS did not converge", iter, rerror
#     end
#     dot(x, r)
# end

# function _exrss(Q::SparseCSC{<:AbstractSymbolic{Tv},Ti}, r::Vector{<:AbstractSymbolic{Tv}}, maxiter, steps, rtol) where {Tv,Ti}
#     dot(stgs(Q, maxiter=maxiter, steps=steps, rtol=rtol), r)
# end

# function _exrss(Q::SparseMatrixCSC{<:AbstractSymbolic{Tv},Ti}, r::Vector{<:AbstractSymbolic{Tv}}, maxiter, steps, rtol) where {Tv,Ti}
#     dot(stgs(Q, maxiter=maxiter, steps=steps, rtol=rtol), r)
# end

# """
# exrt, cexrt
# """

# function exrt(ts::AbstractVector{<:Number}, m::CTMCModel{Tv}; reward, forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv}
#     _exrt(m.Q, m.initv, m.reward[reward], ts, forward, ufact, eps, rmax)
# end

# function cexrt(ts::AbstractVector{<:Number}, m::CTMCModel{Tv}; reward, forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv}
#     _cexrt(m.Q, m.initv, m.reward[reward], ts, forward, ufact, eps, rmax)
# end

# function exrt(ts::Number, m::CTMCModel{Tv}; reward, forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv}
#     _exrt(m.Q, m.initv, m.reward[reward], [ts], forward, ufact, eps, rmax)[1]
# end

# function cexrt(ts::Number, m::CTMCModel{Tv}; reward, forward=:T, ufact=1.01, eps=1.0e-8, rmax=500) where {Tv}
#     _cexrt(m.Q, m.initv, m.reward[reward], [ts], forward, ufact, eps, rmax)[1]
# end

# function _exrt(Q::AbstractMatrix{Tv}, x::AbstractVector{Tv}, r::AbstractVector{Tv}, ts::AbstractVector{<:Number}, forward, ufact, eps, rmax) where {Tv<:Number}
#     result, _, = tran(Q, x, r, ts, forward=forward, ufact=ufact, eps=eps, rmax=rmax)
#     result
# end

# function _cexrt(Q::AbstractMatrix{Tv}, x::AbstractVector{Tv}, r::AbstractVector{Tv}, ts::AbstractVector{<:Number}, forward, ufact, eps, rmax) where {Tv<:Number}
#     _, cresult, = tran(Q, x, r, ts, forward=forward, ufact=ufact, eps=eps, rmax=rmax)
#     cresult
# end

# function _exrt(Q::AbstractMatrix{<:AbstractSymbolic{Tv}}, x::AbstractVector{<:AbstractSymbolic{Tv}}, r::AbstractVector{<:AbstractSymbolic{Tv}}, ts::AbstractVector{<:Number}, forward, ufact, eps, rmax) where {Tv}
#     tran(Q, x, r, ts, forward=forward, ufact=ufact, eps=eps, rmax=rmax, cumulative=false)
# end

# function _cexrt(Q::AbstractMatrix{<:AbstractSymbolic{Tv}}, x::AbstractVector{<:AbstractSymbolic{Tv}}, r::AbstractVector{<:AbstractSymbolic{Tv}}, ts::AbstractVector{<:Number}, forward, ufact, eps, rmax) where {Tv}
#     tran(Q, x, r, ts, forward=forward, ufact=ufact, eps=eps, rmax=rmax, cumulative=true)
# end


text:
using Printf

import Memento: debug

export info, warn, debug, @timing,
       configure_logger, add_file_logger

global LOGGER = Memento.getlogger(@__MODULE__)

const DEFAULT_LOG_LEVEL = "warn"

"""
    info(msg)

Prints a message on info level.

### Input

- `msg` - message string
"""
@inline function info(msg::String)
    Memento.info(LOGGER, msg)
end

"""
    warn(msg)

Prints a message on warn level.

### Input

- `msg` - message string
"""
@inline function warn(msg::String)
    Memento.warn(LOGGER, msg)
end

"""
    debug(msg)

Prints a message on debug level.

### Input

- `msg` - message string
"""
@inline function debug(msg::String)
    Memento.debug(LOGGER, msg)
end

"""
    @timing(expr, [func]=info)

Executes the expression `expr` and prints the elapsed time using the given log
function.

### Input

- `expr`   -- expression
- `func`   -- (optional, default: `info`) log function

### Output

The macro returns the result of evaluating the expression `expr`.
The timing information is printed to the logger.

### Notes

This function is taken from the
[Julia documentation](https://docs.julialang.org/en/v1/manual/metaprogramming/#Hygiene-1).

### Examples

```julia
julia> @timing(1+1)
[info | Reachability]: elapsed time: 1.269e-06 seconds
2
```
"""
macro timing(expr, func=info)
    return quote
        local t0 = time()
        local val = $(esc(expr))
        local t1 = time()
        $func(@sprintf "elapsed time: %1.3e seconds" t1-t0)
        val
    end
end

"""
    configure_logger(level)

Configures the global log level. If no log level is passed, we use the log level
that is defined by the constant `DEFAULT_LOG_LEVEL`.

### Input

- `level` - (optional) the log level; can be either an integer between 0 and 2
            or a string that is supported by the
            [Memento.jl](https://invenia.github.io/Memento.jl/latest/man/intro.html#Logging-levels-1)
            package.
"""
function configure_logger(level::Union{String, Int, Nothing}=DEFAULT_LOG_LEVEL)
    if level isa String
        level_string = level
    elseif level isa Int
        if level == 0
            level_string = "warn"
        elseif level == 1
            level_string = "info"
        elseif level == 2
            level_string = "debug"
        else
            error("Illegal verbosity input $level.")
        end
    else
        error("Illegal verbosity input $level.")
    end
    return Memento.config!(level_string; fmt="[{level}] {msg}")
end

"""
    add_file_logger(filename)

Sets up an additional logger to a file.

### Input

- `filename` -- (optional, default: `tempname()`) the log file name
"""
function add_file_logger(filename::String=tempname())::Nothing
    formatter = Memento.DefaultFormatter("[{date}|{level}] {msg}")
    handler = Memento.DefaultHandler(filename, formatter)
    push!(LOGGER, handler)
    return nothing
end


text:
"""
    DeepEquilibriumNetwork(model, solver; jacobian_regularization::Bool=false,
                           p=nothing, sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10),
                           kwargs...)

Deep Equilibrium Network as proposed in [baideep2019](@cite)

## Arguments

* `model`: Explicit Neural Network which takes 2 inputs
* `solver`: Solver for the optimization problem (See: [`ContinuousDEQSolver`](@ref) & [`DiscreteDEQSolver`](@ref))
* `jacobian_regularization`: If true, Jacobian Loss is computed and stored in the [`DeepEquilibriumSolution`](@ref)
* `p`: Optional parameters for the `model`
* `sensealg`: See [`SteadyStateAdjoint`](@ref)
* `kwargs`: Additional Parameters that are directly passed to `solve`

## Example

```julia
model = DeepEquilibriumNetwork(
    Parallel(
        +,
        Dense(2, 2; bias=false),
        Dense(2, 2; bias=false)
    ),
    ContinuousDEQSolver(VCABM3(); abstol=0.01f0, reltol=0.01f0)
)

model(rand(Float32, 2, 1))
```

See also: [`SkipDeepEquilibriumNetwork`](@ref), [`MultiScaleDeepEquilibriumNetwork`](@ref), [`MultiScaleSkipDeepEquilibriumNetwork`](@ref)
"""
struct DeepEquilibriumNetwork{J,M,P,RE,A,S,K} <: AbstractDeepEquilibriumNetwork
    jacobian_regularization::Bool
    model::M
    p::P
    re::RE
    solver::A
    kwargs::K
    sensealg::S
    stats::DEQTrainingStats

    function DeepEquilibriumNetwork(jacobian_regularization, model, p, re, solver, kwargs, sensealg, stats)
        _p, re = destructure_parameters(model)
        p = p === nothing ? _p : convert(typeof(_p), p)

        return new{jacobian_regularization,typeof(model),typeof(p),typeof(re),typeof(solver),
                   typeof(sensealg),typeof(kwargs)}(jacobian_regularization, model, p, re,
                                                    solver, kwargs, sensealg, stats)
    end
end

Flux.@functor DeepEquilibriumNetwork

function Base.show(io::IO, l::DeepEquilibriumNetwork{J}) where {J}
    return print(io, "DeepEquilibriumNetwork(jacobian_regularization = $J) ",
                 string(length(l.p)), " Trainable Parameters")
end

function DeepEquilibriumNetwork(model, solver; jacobian_regularization::Bool=false,
                                p=nothing, sensealg=SteadyStateAdjoint(0.1f0, 0.1f0, 10), kwargs...)
    return DeepEquilibriumNetwork(jacobian_regularization, model, p, nothing, solver,
                                  kwargs, sensealg, DEQTrainingStats(0))
end

function (deq::DeepEquilibriumNetwork)(x::AbstractArray{T}) where {T}
    z = zero(x)
    Zygote.@ignore deq.re(deq.p)(z, x)

    current_nfe = deq.stats.nfe

    z_star = solve_steady_state_problem(deq.re, deq.p, x, z, deq.sensealg, deq.solver; dudt=nothing,
                                        update_nfe=() -> (deq.stats.nfe += 1), deq.kwargs...)

    jac_loss = (deq.jacobian_regularization ? compute_deq_jacobian_loss(deq.re, deq.p, z_star, x) : T(0))::T

    residual = Zygote.@ignore z_star .- deq.re(deq.p)(z_star, x)

    return z_star, DeepEquilibriumSolution(z_star, z, residual, jac_loss, deq.stats.nfe - current_nfe)
end


text:
import Base: ∈, split
using Base: product

export AbstractHyperrectangle,
       radius_hyperrectangle,
       constraints_list,
       low, high,
       isflat

"""
    AbstractHyperrectangle{N<:Real} <: AbstractZonotope{N}

Abstract type for hyperrectangular sets.

### Notes

See [`Hyperrectangle`](@ref) for a standard implementation of this interface.

Every concrete `AbstractHyperrectangle` must define the following functions:
- `radius_hyperrectangle(::AbstractHyperrectangle{N})::Vector{N}` -- return the
    hyperrectangle's radius, which is a full-dimensional vector
- `radius_hyperrectangle(::AbstractHyperrectangle{N}, i::Int)::N` -- return the
    hyperrectangle's radius in the `i`-th dimension
- `isflat(::AbstractHyperrectangle{N})::Bool` -- determine whether the
    hyperrectangle's radius is zero in some dimension

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractHyperrectangle)
5-element Array{Any,1}:
 AbstractSingleton
 BallInf
 Hyperrectangle
 Interval
 SymmetricIntervalHull
```
"""
abstract type AbstractHyperrectangle{N<:Real} <: AbstractZonotope{N}
end

# --- AbstractZonotope interface functions ---

"""
   genmat(H::AbstractHyperrectangle)

Return the generator matrix of a hyperrectangular set.

### Input

- `H` -- hyperrectangular set

### Output

A matrix where each column represents one generator of `H`.
"""
function genmat(H::AbstractHyperrectangle)
    return genmat_fallback(H)
end

# iterator that wraps the generator matrix
struct HyperrectangleGeneratorIterator{AH<:AbstractHyperrectangle}
    H::AH
    nonflats::Vector{Int}  # dimensions along which `H` is not flat
    dim::Int  # total number of dimensions of `H` (stored for efficiency)

    function HyperrectangleGeneratorIterator(H::AH) where {N<:Real,
            AH<:AbstractHyperrectangle{N}}
        n = dim(H)
        nonflats = Vector{Int}()
        sizehint!(nonflats, n)
        @inbounds for i in 1:n
            if radius_hyperrectangle(H, i) != zero(N)
                push!(nonflats, i)
            end
        end
        return new{AH}(H, nonflats, n)
    end
end

Base.length(it::HyperrectangleGeneratorIterator) = length(it.nonflats)

Base.eltype(::Type{<:HyperrectangleGeneratorIterator{<:AbstractHyperrectangle{N}}}) where {N} =
    SingleEntryVector{N}

function Base.iterate(it::HyperrectangleGeneratorIterator{<:AH},
                      state::Int=1) where {N, AH<:AbstractHyperrectangle{N}}
    if state > length(it.nonflats)
        return nothing
    end
    i = it.nonflats[state]
    r = radius_hyperrectangle(it.H, i)
    g = SingleEntryVector(i, it.dim, r)
    state += 1
    return (g, state)
end

"""
    generators(H::AbstractHyperrectangle)

Return an iterator over the generators of a hyperrectangular set.

### Input

- `H` -- hyperrectangular set

### Output

An iterator over the generators of `H`.
"""
function generators(H::AbstractHyperrectangle)
    return HyperrectangleGeneratorIterator(H)
end

# --- AbstractPolytope interface functions ---

"""
    vertices_list(H::AbstractHyperrectangle{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a hyperrectangular set.

### Input

- `H` -- hyperrectangular set

### Output

A list of vertices.

### Notes

For high dimensions, it is preferable to develop a `vertex_iterator` approach.
"""
function vertices_list(H::AbstractHyperrectangle{N}
                      )::Vector{Vector{N}} where {N<:Real}
    # fast evaluation if H has radius 0
    if radius_hyperrectangle(H) == zeros(N, dim(H))
        return [center(H)]
    end
    return [center(H) .+ si .* radius_hyperrectangle(H)
        for si in Iterators.product([[1, -1] for i = 1:dim(H)]...)][:]
end

"""
    constraints_list(H::AbstractHyperrectangle{N}) where {N<:Real}

Return the list of constraints of an axis-aligned hyperrectangular set.

### Input

- `H` -- hyperrectangular set

### Output

A list of linear constraints.
"""
function constraints_list(H::AbstractHyperrectangle{N}) where {N<:Real}
    n = dim(H)
    constraints = Vector{LinearConstraint{N}}(undef, 2*n)
    b, c = high(H), -low(H)
    one_N = one(N)
    for i in 1:n
        ei = SingleEntryVector(i, n, one_N)
        constraints[i] = HalfSpace(ei, b[i])
        constraints[i+n] = HalfSpace(-ei, c[i])
    end
    return constraints
end

# --- LazySet interface functions ---

"""
    σ(d::AbstractVector{N}, H::AbstractHyperrectangle{N}) where {N<:Real}

Return the support vector of a hyperrectangular set in a given direction.

### Input

- `d` -- direction
- `H` -- hyperrectangular set

### Output

The support vector in the given direction.
If the direction has norm zero, the vertex with biggest values is returned.
"""
function σ(d::AbstractVector{N}, H::AbstractHyperrectangle{N}) where {N<:Real}
    @assert length(d) == dim(H) "a $(length(d))-dimensional vector is " *
                                "incompatible with a $(dim(H))-dimensional set"
    return center(H) .+ sign_cadlag.(d) .* radius_hyperrectangle(H)
end

"""
    ρ(d::AbstractVector{N}, H::AbstractHyperrectangle{N}) where {N<:Real}

Evaluate the support function of a hyperrectangular set in a given direction.

### Input

- `d` -- direction
- `H` -- hyperrectangular set

### Output

Evaluation of the support function in the given direction.
"""
function ρ(d::AbstractVector{N}, H::AbstractHyperrectangle{N}) where {N<:Real}
    @assert length(d) == dim(H) "a $(length(d))-dimensional vector is " *
                                "incompatible with a $(dim(H))-dimensional set"
    c = center(H)
    res = zero(N)
    @inbounds for (i, di) in enumerate(d)
        if di < zero(N)
            res += di * (c[i] - radius_hyperrectangle(H, i))
        elseif di > zero(N)
            res += di * (c[i] + radius_hyperrectangle(H, i))
        end
    end
    return res
end

"""
    norm(H::AbstractHyperrectangle, [p]::Real=Inf)::Real

Return the norm of a hyperrectangular set.

The norm of a hyperrectangular set is defined as the norm of the enclosing ball,
of the given ``p``-norm, of minimal volume that is centered in the origin.

### Input

- `H` -- hyperrectangular set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the norm.

### Algorithm

Recall that the norm is defined as

```math
‖ X ‖ = \\max_{x ∈ X} ‖ x ‖_p = max_{x ∈ \\text{vertices}(X)} ‖ x ‖_p.
```
The last equality holds because the optimum of a convex function over a polytope
is attained at one of its vertices.

This implementation uses the fact that the maximum is achieved in the vertex
``c + \\text{diag}(\\text{sign}(c)) r``, for any ``p``-norm, hence it suffices to
take the ``p``-norm of this particular vertex. This statement is proved below.
Note that, in particular, there is no need to compute the ``p``-norm for *each*
vertex, which can be very expensive.

If ``X`` is an axis-aligned hyperrectangle and the ``n``-dimensional vectors center
and radius of the hyperrectangle are denoted ``c`` and ``r`` respectively, then
reasoning on the ``2^n`` vertices we have that:

```math
\\max_{x ∈ \\text{vertices}(X)} ‖ x ‖_p = \\max_{α_1, …, α_n ∈ \\{-1, 1\\}} (|c_1 + α_1 r_1|^p + ... + |c_n + α_n r_n|^p)^{1/p}.
```

The function ``x ↦ x^p``, ``p > 0``, is monotonically increasing and thus the
maximum of each term ``|c_i + α_i r_i|^p`` is given by ``|c_i + \\text{sign}(c_i) r_i|^p``
for each ``i``. Hence, ``x^* := \\text{argmax}_{x ∈ X} ‖ x ‖_p`` is the vertex
``c + \\text{diag}(\\text{sign}(c)) r``.
"""
function norm(H::AbstractHyperrectangle, p::Real=Inf)::Real
    c, r = center(H), radius_hyperrectangle(H)
    return norm((@. c + sign_cadlag(c) * r), p)
end

"""
    radius(H::AbstractHyperrectangle, [p]::Real=Inf)::Real

Return the radius of a hyperrectangular set.

### Input

- `H` -- hyperrectangular set
- `p` -- (optional, default: `Inf`) norm

### Output

A real number representing the radius.

### Notes

The radius is defined as the radius of the enclosing ball of the given
``p``-norm of minimal volume with the same center.
It is the same for all corners of a hyperrectangular set.
"""
function radius(H::AbstractHyperrectangle, p::Real=Inf)::Real
    return norm(radius_hyperrectangle(H), p)
end

"""
    ∈(x::AbstractVector{N}, H::AbstractHyperrectangle{N})::Bool where {N<:Real}

Check whether a given point is contained in a hyperrectangular set.

### Input

- `x` -- point/vector
- `H` -- hyperrectangular set

### Output

`true` iff ``x ∈ H``.

### Algorithm

Let ``H`` be an ``n``-dimensional hyperrectangular set, ``c_i`` and ``r_i`` be
the box's center and radius and ``x_i`` be the vector ``x`` in dimension ``i``,
respectively.
Then ``x ∈ H`` iff ``|c_i - x_i| ≤ r_i`` for all ``i=1,…,n``.
"""
function ∈(x::AbstractVector{N},
           H::AbstractHyperrectangle{N})::Bool where {N<:Real}
    @assert length(x) == dim(H)
    for i in eachindex(x)
        if abs(center(H)[i] - x[i]) > radius_hyperrectangle(H, i)
            return false
        end
    end
    return true
end

# --- common AbstractHyperrectangle functions ---

"""
    high(H::AbstractHyperrectangle{N})::Vector{N} where {N<:Real}

Return the higher coordinates of a hyperrectangular set.

### Input

- `H` -- hyperrectangular set

### Output

A vector with the higher coordinates of the hyperrectangular set.
"""
function high(H::AbstractHyperrectangle{N})::Vector{N} where {N<:Real}
    return center(H) .+ radius_hyperrectangle(H)
end

"""
    high(H::AbstractHyperrectangle{N}, i::Int)::N where {N<:Real}

Return the higher coordinate of a hyperrectangular set in a given dimension.

### Input

- `H` -- hyperrectangular set
- `i` -- dimension of interest

### Output

The higher coordinate of the hyperrectangular set in the given dimension.
"""
function high(H::AbstractHyperrectangle{N}, i::Int)::N where {N<:Real}
    return center(H)[i] + radius_hyperrectangle(H, i)
end

"""
    low(H::AbstractHyperrectangle{N})::Vector{N} where {N<:Real}

Return the lower coordinates of a hyperrectangular set.

### Input

- `H` -- hyperrectangular set

### Output

A vector with the lower coordinates of the hyperrectangular set.
"""
function low(H::AbstractHyperrectangle{N})::Vector{N} where {N<:Real}
    return center(H) .- radius_hyperrectangle(H)
end

"""
    low(H::AbstractHyperrectangle{N}, i::Int)::N where {N<:Real}

Return the lower coordinate of a hyperrectangular set in a given dimension.

### Input

- `H` -- hyperrectangular set
- `i` -- dimension of interest

### Output

The lower coordinate of the hyperrectangular set in the given dimension.
"""
function low(H::AbstractHyperrectangle{N}, i::Int)::N where {N<:Real}
    return center(H)[i] - radius_hyperrectangle(H, i)
end

"""
    isflat(H::AbstractHyperrectangle)::Bool

Determine whether a hyperrectangular set is flat, i.e. whether its radius
is zero in some dimension.

### Input

- `H` -- hyperrectangular set

### Output

`true` iff the hyperrectangular set is flat.

### Notes

For robustness with respect to floating-point inputs, this function relies on
the result of `isapproxzero` when applied to the radius in some dimension.
Hence, this function depends on the absolute zero tolerance `ABSZTOL`.
"""
function isflat(H::AbstractHyperrectangle)::Bool
    return any(i -> isapproxzero(radius_hyperrectangle(H, i)), 1:dim(H))
end

"""
    split(H::AbstractHyperrectangle{N}, num_blocks::AbstractVector{Int}
         ) where {N<:Real}

Partition a hyperrectangular set into uniform sub-hyperrectangles.

### Input

- `H`          -- hyperrectangular set
- `num_blocks` -- number of blocks in the partition for each dimension

### Output

A list of `Hyperrectangle`s.
"""
function split(H::AbstractHyperrectangle{N}, num_blocks::AbstractVector{Int}
              )::Vector{Hyperrectangle{N}} where {N<:Real}
    @assert length(num_blocks) == dim(H) "need number of blocks in each dimension"
    radius = copy(radius_hyperrectangle(H))
    total_number = 1
    lo = low(H)
    hi = high(H)

    # precompute center points in each dimension
    centers = Vector{StepRangeLen{N}}(undef, dim(H))
    for (i, m) in enumerate(num_blocks)
        if m <= 0
            throw(ArgumentError(m, "each dimension needs at least one block"))
        elseif m == one(N)
            centers[i] = range(lo[i] + radius[i], length=1)
        else
            radius[i] /= m
            centers[i] = range(lo[i] + radius[i], step=(2 * radius[i]),
                               length=m)
            total_number *= m
        end
    end

    # create hyperrectangles for every combination of the center points
    result = Vector{Hyperrectangle{N}}(undef, total_number)
    for (i, center) in enumerate(product(centers...))
        result[i] = Hyperrectangle(collect(center), copy(radius))
    end
    return result
end


text:
# ZEROFPR
#
# PART OF BAZINGA.jl

using Base.Iterators
using ProximalAlgorithms.IterationTools
using LinearAlgebra: norm
using OptiMo
using Printf

###########################################################################
# Options
###########################################################################
struct ZEROFPR_options{R<:Real}
    alpha::R                        # ∈ (0, 1), e.g.: 0.95
    beta::R                         # ∈ (0, 1), e.g.: 0.5
    sigma::R                        # = beta * (1 - alpha) / 2 ∈ (0, 0.5)
    adaptive::Bool                  # enforce adaptive stepsize even if L is provided
    tau_min::R                      # minimum linesearch stepsize
    gamma_min::R                    # minimum forward-backward stepsize
    bigeps::R                       # small value used for tolerances
    grad_scaling::Bool              # perform gradient-based scaling?
    scaling_grad_max::R             # maximum gradient after scaling
    scaling_grad_target::Maybe{R}   # target value of gradient after scaling
    scaling_min_value::R            # minimum value of gradient-based scaling factor
    warnings::Bool                  # can print warnings?
end

###########################################################################
# Iterable
###########################################################################
struct ZEROFPR_iterable{R<:Real,Tx<:AbstractArray{R},TP,TH,TO<:ZEROFPR_options{R}}
    prob::TP                        # problem
    x0::Tx                          # initial point
    gamma::Maybe{R}                 # forward-backward stepsize
    H::TH                           # quasi-Newton object
    opts::TO                        # options
end

Base.IteratorSize(::Type{<:ZEROFPR_iterable}) = Base.IsInfinite()

###########################################################################
# State
###########################################################################
mutable struct ZEROFPR_state{R<:Real,Tx,TH}
    x::Tx                   # iterate
    f_x::R                  # (unscaled) value of smooth term
    grad_f_x::Tx            # (unscaled) gradient of smooth term
    gamma::R                # stepsize parameter of forward and backward steps
    y::Tx                   # forward point
    xbar::Tx                # forward-backward point
    g_xbar::R               # (unscaled) value of nonsmooth term (at xbar)
    res::Tx                 # fixed-point residual at iterate (= x - xbar)
    optim::R                # optimality measure
    H::TH                   # variable metric
    tau::Maybe{R}           # linesearch stepsize
    objec::R                # (unscaled) objective function at `xbar`
    scalobj::R              # scaling factor (> 0)
    fbe_x::R                # (scaled) forward-backward envelope at `x`
    # additional allocation
    grad_f_xbar::Tx
    xbarbar::Tx
    res_xbar::Tx
    d::Tx                   # search direction at `xbar`
    xbar_old::Maybe{Tx}
    res_xbar_old::Maybe{Tx}
end

###########################################################################
# Quadratic model
###########################################################################
f_model(state::ZEROFPR_state) = f_model(
    state.scalobj * state.f_x,
    state.scalobj .* state.grad_f_x,
    state.res,
    state.gamma,
)

###########################################################################
# Initialize
###########################################################################
function Base.iterate(iter::ZEROFPR_iterable{R}) where {R}
    # initial point
    x = copy(iter.x0)
    g_xbar, x = objprox(iter.prob, x, iter.opts.bigeps)
    f_x, grad_f_x = objgrad(iter.prob, x)
    # scaling
    if iter.opts.grad_scaling
        scalobj = objgradscaling(
            grad_f_x,
            iter.opts.scaling_grad_target,
            iter.opts.scaling_grad_max,
            iter.opts.scaling_min_value,
        )
    else
        scalobj = R(1)
    end
    # forward-backward stepsize
    gamma = iter.gamma
    if gamma === nothing
        # compute lower bound to Lipschitz constant of x ↦ ∇f(x)
        xh = x .+ 0.01 .* (x .+ R(1))
        grad_f_xh = grad(iter.prob, xh)
        L = norm(scalobj .* (grad_f_xh - grad_f_x)) / R(sqrt(length(x)))
        L = max(1.0, L) # this avoids zero-division
        gamma = max(iter.opts.alpha / L, sqrt(iter.opts.gamma_min))
    end

    # initial forward-backward step
    y = x - gamma .* (scalobj .* grad_f_x)
    g_xbar, xbar = objprox(iter.prob, y, gamma * scalobj)

    # initial fixed-point residual and optimality
    res = x - xbar
    #optim = norm( res, Inf ) / gamma
    f_xbar, tmp = objgrad(iter.prob, xbar) # grad_f_xbar
    tmp .+= (res ./ gamma) .- grad_f_x
    optim = norm(tmp, Inf)

    # initial objective
    objec = f_xbar + g_xbar

    # forward-backward envelope
    f_xbar_upp = f_model(scalobj * f_x, scalobj .* grad_f_x, res, gamma)
    fbe_x = f_xbar_upp + scalobj * g_xbar

    state = ZEROFPR_state(
        x,
        f_x,
        grad_f_x,
        gamma,
        y,
        xbar,
        g_xbar,
        res,
        optim,
        iter.H,
        nothing,
        objec,
        scalobj,
        fbe_x,
        zero(x),
        zero(x),
        zero(x),
        zero(x),
        nothing,
        nothing,
    )

    return state, state
end

###########################################################################
# Step
###########################################################################
function Base.iterate(iter::ZEROFPR_iterable{R}, state::ZEROFPR_state{R,Tx}) where {R,Tx}

    # (scaled) quadratic model of `f` around `x` at `xbar`
    f_xbar_upp = f_model(state)

    # `f` at `xbar`
    f_xbar = objgrad!(iter.prob, state.xbar, state.grad_f_xbar)
    f_xbar *= state.scalobj

    # check and possibly backtrack stepsize `gamma`
    while iter.gamma === nothing || iter.opts.adaptive == true
        if state.gamma < iter.opts.gamma_min
            if iter.opts.warnings
                @warn "stepsize `gamma` too small ($(state.gamma))"
            end
            return nothing
        end
        tol = iter.opts.bigeps * (1 + abs(f_xbar))
        if f_xbar <= f_xbar_upp + tol
            break
        end
        state.gamma *= 0.5
        state.y .= state.x .- state.gamma .* (state.scalobj .* state.grad_f_x)
        state.g_xbar = objprox!(iter.prob, state.y, state.gamma * state.scalobj, state.xbar)
        state.res .= state.x .- state.xbar
        reset!(state.H)
        f_xbar_upp = f_model(state)
        f_xbar = objgrad!(iter.prob, state.xbar, state.grad_f_xbar)
        f_xbar *= state.scalobj
    end

    # (scaled) forward-backward envelope
    state.fbe_x = f_xbar_upp + state.scalobj * state.g_xbar

    # residual at `xbar`
    state.y .= state.xbar .- state.gamma .* (state.scalobj .* state.grad_f_xbar)
    g_xbarbar = objprox!(iter.prob, state.y, state.gamma * state.scalobj, state.xbarbar)
    state.res_xbar .= state.xbar .- state.xbarbar

    if state.xbar_old === nothing
        # store vectors for next update
        state.xbar_old = copy(state.xbar)
        state.res_xbar_old = copy(state.res_xbar)
    else
        # update variable metric
        update!(state.H, state.xbar - state.xbar_old, state.res_xbar - state.res_xbar_old)
        # store vectors for next update
        copyto!(state.xbar_old, state.xbar)
        copyto!(state.res_xbar_old, state.res_xbar)
    end

    # search direction
    mul!(state.d, state.H, -state.res_xbar)

    # linesearch over the forward-backward envelope
    tau = R(1)
    sigma = iter.opts.sigma / state.gamma
    tol = iter.opts.bigeps * (1 + abs(state.fbe_x))
    rhs = state.fbe_x - sigma * norm(state.res)^2 + tol

    while true
        state.x .= state.xbar_old .+ tau .* state.d
        state.f_x = objgrad!(iter.prob, state.x, state.grad_f_x)
        state.y .= state.x .- state.gamma .* (state.scalobj .* state.grad_f_x)
        state.g_xbar = objprox!(iter.prob, state.y, state.gamma * state.scalobj, state.xbar)
        state.res .= state.x .- state.xbar
        f_xbar_upp = f_model(state)
        state.fbe_x = f_xbar_upp + state.scalobj * state.g_xbar

        f_xbar = obj(iter.prob, state.xbar)
        f_xbar *= state.scalobj
        tol = iter.opts.bigeps * (1 + abs(f_xbar))
        is_gamma_ok = f_xbar <= f_xbar_upp + tol

        if (state.fbe_x <= rhs) && is_gamma_ok
            break

        elseif tau < iter.opts.tau_min
            if !is_gamma_ok
                state.gamma *= 0.5
                reset!(state.H)
            end
            objec_t = f_xbar + state.g_xbar
            tol = iter.opts.bigeps * (1 + abs(state.objec))
            if objec_t <= state.objec + tol
                # small linesearch stepsize : taking tau > 0, which gives decrease
                # in the objective function.
                break
            end
            # small linesearch stepsize : taking `tau = 0`, which gives decrease in
            # the FBE.
            tau = 0.0
            copyto!(state.x, state.xbar_old)
            state.f_x = objgrad!(iter.prob, state.x, state.grad_f_x)
            state.y .= state.x .- state.gamma .* (state.scalobj .* state.grad_f_x)
            state.g_xbar =
                objprox!(iter.prob, state.y, state.gamma * state.scalobj, state.xbar)
            state.res .= state.x .- state.xbar
            state.fbe_x = f_model(state) + state.scalobj * state.g_xbar
            break

        else
            tau *= 0.5
        end
    end

    state.tau = tau
    f_xbar, tmp = objgrad(iter.prob, state.xbar) # grad_f_xbar
    tmp .+= (state.res ./ state.gamma) .- state.grad_f_x
    state.optim = norm(tmp, Inf)
    state.objec = f_xbar + state.g_xbar
    return state, state
end

###########################################################################
# Solver
###########################################################################
mutable struct ZEROFPR{R<:Real}
    memory::Int
    max_iter::Int
    tol_optim::R
    verbose::Bool
    freq::Int
    tol_diverging::R
    opts::ZEROFPR_options

    function ZEROFPR{R}(;
        alpha::R = R(0.95),
        beta::R = R(0.5),
        adaptive::Bool = true,
        memory::Int = 8,
        max_iter::Int = 10000,
        tol_optim::R = R(1e-8),
        verbose::Bool = false,
        freq::Int = 10,
        tau_min::R = R(1e-8),
        gamma_min::R = R(1e-12),
        grad_scaling::Bool = true,
        scaling_grad_max::R = R(100),
        scaling_grad_target::Maybe{R} = nothing,
        scaling_min_value::R = R(1e-8),
        tol_diverging::R = R(1e20),
        bigeps::R = R(1e-12),
        warnings::Bool = true,
    ) where {R}
        @assert 0 < alpha < 1
        @assert 0 < beta < 1
        @assert 0 <= memory
        @assert 0 < max_iter
        @assert 0 < freq
        @assert 0 < bigeps < 1
        @assert 0 < tol_optim < 1 < tol_diverging
        @assert 0 < scaling_grad_max
        @assert scaling_grad_target === nothing || 0 < scaling_grad_target
        @assert 0 < scaling_min_value < 1
        @assert 0 < gamma_min < 1
        @assert 0 < tau_min < 1
        sigma = beta * 0.5 * (1.0 - alpha)
        opts = ZEROFPR_options(
            alpha,
            beta,
            sigma,
            adaptive,
            tau_min,
            gamma_min,
            bigeps,
            grad_scaling,
            scaling_grad_max,
            scaling_grad_target,
            scaling_min_value,
            warnings,
        )
        new(memory, max_iter, tol_optim, verbose, freq, tol_diverging, opts)
    end
end

###########################################################################
# Solver call
###########################################################################
function (solver::ZEROFPR{R})(
    prob::TP;
    x0::AbstractArray{R} = prob.meta.x0,
    L::Maybe{R} = nothing,
) where {R,TP<:AbstractOptiModel}

    tstart = time()

    stop_solved(state::ZEROFPR_state) = state.optim <= solver.tol_optim
    stop_unbound(state::ZEROFPR_state) = state.objec < -solver.tol_diverging
    is_illegal(x) = x === nothing ? false : (any(isnan.(x)) || any(isinf.(x)))
    stop_illegal(state::ZEROFPR_state) = is_illegal(state.x) || is_illegal(state.xbar)
    stop(state::ZEROFPR_state) =
        stop_solved(state) || stop_unbound(state) || stop_illegal(state)

    if solver.verbose
        @printf(
            "%6s | %8s | %8s  %9s | %9s  %9s\n",
            "iter",
            "objec",
            "fbe",
            "optim",
            "gamma",
            "tau"
        )
        disp((i, state)) = @printf(
            "%6d | %+.1e | %+.1e  %.3e | %.3e  %.3e \n",
            i,
            state.objec,
            state.fbe_x,
            state.optim,
            state.gamma,
            (state.tau === nothing ? 0.0 : state.tau)
        )
    end

    gamma = (L === nothing ? nothing : solver.opts.alpha / L)

    iter = ZEROFPR_iterable(prob, x0, gamma, LBFGS(x0, solver.memory), solver.opts)
    iter = halt(iter, stop)
    iter = take(iter, solver.max_iter)
    iter = enumerate(iter)
    if solver.verbose
        iter = sample(iter, solver.freq)
        iter = tee(iter, disp)
    end

    num_iters, state = loop(iter)

    status = if stop_solved(state)
        :first_order
    elseif stop_unbound(state)
        :unbounded
    elseif stop_illegal(state)
        :ieee_nan_inf
    elseif (num_iters >= solver.max_iter)
        :max_iter
    elseif (state.gamma <= solver.opts.gamma_min)
        :stalled
    else
        :unknown
    end

    return OptiOutput(
        status,
        x = state.xbar,
        objective = state.objec,
        optimality = state.optim,
        iterations = num_iters,
        time = time() - tstart,
        solver_name = "ZEROFPR",
        solver = Dict(:gx => state.g_xbar),
    )
end

###########################################################################
# Outer constructors
###########################################################################
"""
    ZEROFPR([gamma, adaptive, memory, max_iter, tol_optim, verbose, freq, ...])

Instantiate the ZeroFPR algorithm [TSP18] for solving optimization problems
of the form

    minimize f(x) + g(x),

where `f` is smooth and `g` is proximable. If `solver = ZEROFPR(args...)`, then the
above problem is solved with

    solver( prob, [x0, L] )

Optional keyword arguments:

* `gamma::Real` (default: `nothing`), the stepsize to use; defaults to `alpha/L` if not set (but `L` is).
* `adaptive::Bool` (default: `true`), if true, forces the method stepsize to be adaptively adjusted even
    if `L` is provided (this behaviour is always enforced if `L` is not provided).
* `memory::Integer` (default: `5`), memory parameter for L-BFGS.
* `max_iter::Integer` (default: `1000`), maximum number of iterations to perform.
* `tol_optim::Real` (default: `1e-8`), absolute tolerance on the fixed-point residual.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `10`), frequency of verbosity.
* `alpha::Real` (default: `0.95`), stepsize to inverse-Lipschitz-constant ratio; should be in (0, 1).
* `beta::Real` (default: `0.5`), sufficient decrease parameter; should be in (0, 1).
* `grad_scaling::Bool`
* `scaling_grad_max::Real`
* `scaling_grad_target::Maybe{Real}`
* `scaling_min_value::Real`

If `gamma` is not specified at construction time, the following keyword
argument can be used to set the stepsize parameter:

* `L::Real` (default: `nothing`), the Lipschitz constant of the gradient of x ↦ f(Ax).

References:

[TSP18] Themelis, Stella, Patrinos, ``Forward-backward envelope for the sum of two
        nonconvex functions: Further properties and nonmonotone line-search algorithms``,
        SIAM Journal on Optimization, vol. 28, no. 3, pp. 2274–2303 (2018).
        DOI: 10.1137/16M1080240
"""
ZEROFPR(::Type{R}; kwargs...) where {R} = ZEROFPR{R}(; kwargs...)
ZEROFPR(; kwargs...) = ZEROFPR(Float64; kwargs...)


text:
using BinaryBuilder

name = "UserNSSandbox"
version = v"2021.04.22"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/staticfloat/Sandbox.jl",
              "0fc3812cd74a355e5a2c4d45d68fabb9fb059f55"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/Sandbox.jl/deps
mkdir -p ${bindir}
$CC -std=c99 -O2 -static -static-libgcc -g -o ${bindir}/sandbox ./userns_sandbox.c
install_license /usr/share/licenses/MIT
"""

# We only build for Linux
platforms = filter(p -> Sys.islinux(p), supported_platforms())

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

# Dependencies that must be installed before this package can be built
build_tarballs(ARGS, name, version, sources, script, platforms, products, Dependency[])

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(LoopVectorization, Symbol("#1#4")) && precompile(Tuple{getfield(LoopVectorization, Symbol("#1#4")),Expr})
    isdefined(LoopVectorization, Symbol("#1#4")) && precompile(Tuple{getfield(LoopVectorization, Symbol("#1#4")),Expr})
    isdefined(MacroTools, Symbol("#19#20")) && precompile(Tuple{getfield(MacroTools, Symbol("#19#20")),Int64})
    isdefined(MacroTools, Symbol("#19#20")) && precompile(Tuple{getfield(MacroTools, Symbol("#19#20")),Symbol})
    isdefined(MacroTools, Symbol("#19#20")) && precompile(Tuple{getfield(MacroTools, Symbol("#19#20")),Symbol})
    isdefined(MacroTools, Symbol("#19#20")) && precompile(Tuple{getfield(MacroTools, Symbol("#19#20")),Symbol})
    isdefined(MacroTools, Symbol("#19#20")) && precompile(Tuple{getfield(MacroTools, Symbol("#19#20")),Symbol})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Float64})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Int64})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Int64})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),LineNumberNode})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),LineNumberNode})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Module})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Module})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),QuoteNode})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),QuoteNode})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Symbol})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Symbol})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Symbol})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Symbol})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Type{T} where T})
    isdefined(MacroTools, Symbol("#21#22")) && precompile(Tuple{getfield(MacroTools, Symbol("#21#22")),Type{T} where T})
    precompile(Tuple{Core.kwftype(typeof(LoopVectorization._vectorloads!)),NamedTuple{(:itersym, :declared_iter_sym, :VectorizationDict, :mod),Tuple{Symbol,Symbol,Dict{Symbol,Tuple{Symbol,Symbol}},Module}},typeof(LoopVectorization._vectorloads!),Expr,Expr,Tuple{Dict{Symbol,Symbol},Dict{Tuple{Symbol,Symbol},Symbol},Dict{Expr,Symbol},Dict{Expr,Symbol}},Type{NTuple{8,VecElement{Float64}}},Expr,Expr})
    precompile(Tuple{Core.kwftype(typeof(LoopVectorization._vectorloads!)),NamedTuple{(:itersym, :declared_iter_sym, :VectorizationDict, :mod),Tuple{Symbol,Symbol,Dict{Symbol,Tuple{Symbol,Symbol}},Symbol}},typeof(LoopVectorization._vectorloads!),Expr,Expr,Tuple{Dict{Symbol,Symbol},Dict{Tuple{Symbol,Symbol},Symbol},Dict{Expr,Symbol},Dict{Expr,Symbol}},Type{NTuple{8,VecElement{Float64}}},Expr,Expr})
    precompile(Tuple{typeof(LoopVectorization.add_masks),Expr,Symbol,Dict{Tuple{Symbol,Symbol},Symbol},Module})
    precompile(Tuple{typeof(LoopVectorization.add_masks),Expr,Symbol,Dict{Tuple{Symbol,Symbol},Symbol},Symbol})
    precompile(Tuple{typeof(LoopVectorization.vectorize_assign_linear_index),Symbol,Expr,Expr,Dict{Symbol,Symbol},Symbol,Symbol,Symbol})
    precompile(Tuple{typeof(LoopVectorization.vectorize_assign_linear_index),Symbol,Expr,Symbol,Dict{Symbol,Symbol},Symbol,Symbol,Module})
    precompile(Tuple{typeof(LoopVectorization.vectorize_assign_linear_index),Symbol,Expr,Symbol,Dict{Symbol,Symbol},Symbol,Symbol,Symbol})
    precompile(Tuple{typeof(LoopVectorization.vectorize_assign_linear_index),Symbol,Symbol,Symbol,Dict{Symbol,Symbol},Symbol,Symbol,Symbol})
    precompile(Tuple{typeof(LoopVectorization.vectorize_body),Int64,Int64,Type{NTuple{8,VecElement{Float64}}},Int64,Symbol,Array{Any,1},Dict{Symbol,Tuple{Symbol,Symbol}},Bool,Int64,Int64,Module})
    precompile(Tuple{typeof(LoopVectorization.vectorize_body),Int64,Int64,Type{NTuple{8,VecElement{Float64}}},Int64,Symbol,Array{Any,1},Dict{Symbol,Tuple{Symbol,Symbol}},Bool,Int64,Int64,Symbol})
    precompile(Tuple{typeof(LoopVectorization.vectorize_body),Int64,Type{Float64},Int64,Symbol,Array{Any,1},Dict{Symbol,Tuple{Symbol,Symbol}},Any,Bool,Module})
    precompile(Tuple{typeof(LoopVectorization.vectorize_body),Int64,Type{Float64},Int64,Symbol,Array{Any,1},Dict{Symbol,Tuple{Symbol,Symbol}},Any,Bool})
    precompile(Tuple{typeof(LoopVectorization.vectorize_body),Symbol,Symbol,Type{NTuple{8,VecElement{Float64}}},Int64,Symbol,Array{Any,1},Dict{Symbol,Tuple{Symbol,Symbol}},Bool,Int64,Int64,Symbol})
    precompile(Tuple{typeof(LoopVectorization.vectorize_body),Symbol,Type{Float64},Int64,Symbol,Array{Any,1},Dict{Symbol,Tuple{Symbol,Symbol}},Any,Bool})
    precompile(Tuple{typeof(LoopVectorization.vectorize_linear_index!),Expr,Dict{Expr,Symbol},Dict{Symbol,Symbol},Symbol,Expr,Symbol,Symbol,Symbol,Type{T} where T})
    precompile(Tuple{typeof(LoopVectorization.vectorize_linear_index!),Expr,Dict{Expr,Symbol},Dict{Symbol,Symbol},Symbol,Symbol,Symbol,Symbol,Module,Type{T} where T})
    precompile(Tuple{typeof(LoopVectorization.vectorize_linear_index!),Expr,Dict{Expr,Symbol},Dict{Symbol,Symbol},Symbol,Symbol,Symbol,Symbol,Symbol,Type{T} where T})
end


text:
# Moved from ModiaBase/src/Symbolic.jl, because dependent on FloatType that is only available in ModiaLang

prependPar(ex, prefix, parameters=[], inputs=[]) = ex

function prependPar(ex::Symbol, prefix, parameters=[], inputs=[])
    if prefix == nothing; ex elseif ex in [:time, :instantiatedModel, :_leq_mode, :_x]; ex else Expr(:ref, prefix, QuoteNode(ex)) end
end

function prependPar(ex::Expr, prefix, parameters=[], inputs=[])
    if isexpr(ex, :.)
        Expr(:ref, prependPar(ex.args[1], prefix, parameters, inputs), QuoteNode(ex.args[2].value))
    else
        Expr(ex.head, [prependPar(arg, prefix, parameters, inputs) for arg in ex.args]...)
    end
end

function addCastAndTypeCheckBasicValueType(ex, value, FloatType)
    # typeof(FloatType) is neither Measurements nor MonteCarloMeasurements
    T = baseType(FloatType)
    if isQuantity(typeof(value))
        ustr = unitAsString(unit(value))
        :( Modia.quantity($T, @u_str($ustr))($ex)::Modia.quantity($T, @u_str($ustr)) )
    else
        :( $T($ex)::$T )
    end
end

"""
    e = addCastAndTypeCheck(ex,value,FloatType)

If `eltype(value) <: Number` add either

1. RequiredType(ex)::RequiredType or
2. typeof(value)(ex)::typeof(value)

If it is manageable, variant (1) is used and the type cast is made with respect to _FloatType.
This is done for scalars that are no Measurements and no MonteCarloMeasurements.
Otherwise, variant (2) is used.

As a result, the generated code is type stable and does not work with Any types
(although the values are stored in an Any dictionary).
So the code will be more efficient and memory allocated at run-time will be reduced.
Furthermore, if value is a scalar that has a unit, the conversion includes also the conversion
to the unit used during code generation.

Note, this function should only be used on parameter values.
"""
function addCastAndTypeCheck(ex,value,FloatType)
    if eltype(value) <: Number
        # Numeric scalar or array
        valueType = typeof(value)

        if valueType <: Number
            # Numeric scalar
            T = baseType(FloatType)
            if FloatType <: Measurements.Measurement
                if isMeasurements(valueType)
                    if isQuantity(valueType)
                        ustr = unitAsString(unit(value))
                        :( Modia.quantity(Modia.Measurements.Measurement{$T}, @u_str($ustr))($ex)::Modia.quantity(Modia.Measurements.Measurement{$T}, @u_str($ustr)) )
                    else
                        :( Modia.Measurements.Measurement{$T}($ex)::Modia.Measurements.Measurement{$T} )
                    end
                else
                    addCastAndTypeCheckBasicValueType(ex,value,T)
                end

            elseif FloatType <: MonteCarloMeasurements.AbstractParticles
                if isMonteCarloMeasurements(valueType)
                    if isQuantity(valueType)
                        error("Error in Modia (MonteCarloMeasurements with Quantities appears - which is not supported.)")
                    else
                        :( ($ex)::$valueType )
                    end
                else
                    addCastAndTypeCheckBasicValueType(ex,value,T)
                end

            elseif valueType <: AbstractFloat || isQuantity(valueType)
                # FloatType is neither Measurements nor MonteCarloMeasurements and
                # is either an AbstractFloat or has a unit (is treated as Float)
                addCastAndTypeCheckBasicValueType(ex,value,T)

            else
                # For example Bool or Int values: No conversion, just define the type in the code
                :( ($ex)::$valueType )
            end
        else
            # Numeric array
            :( ($ex)::$valueType )
        end
    else
        # No numeric scalar or array
        ex
    end
end

"""
    e = makeDerVar(ex, parameters, inputs, FloatType, evaluateParameters=false)

Recursively converts der(x) to Symbol(:(der(x))) in expression `ex`

* `ex`: Expression or array of expressions
* `return `e`: ex with der(x) converted
"""
function makeDerVar(ex, parameters, inputs, FloatType, evaluateParameters=false)
    if typeof(ex) in [Symbol, Expr]
        if ex in keys(parameters)
            addCastAndTypeCheck( prependPar(ex, :(_p), parameters, inputs),  parameters[ex], FloatType )
        elseif ex in keys(inputs)
            prependPar(ex, :(_p), parameters, inputs)
        else
            ex
        end
    else
        ex
    end
end

function makeDerVar(ex::Expr, parameters, inputs, FloatType, evaluateParameters=false)
    if ex.head == :call && ex.args[1] == :der
        Symbol(ex)
	elseif isexpr(ex, :.) && ex in keys(parameters)
        if evaluateParameters
            parameters[ex]
        else
            addCastAndTypeCheck( prependPar(ex, :(_p), parameters, inputs), parameters[ex], FloatType )
        end
	elseif isexpr(ex, :.) && ex in keys(inputs)
        if evaluateParameters
            inputs[ex]
        else
            prependPar(ex, :(_p), parameters, inputs)
        end
    elseif ex.head == :.
        Symbol(ex)
    elseif ex.head == :call # Don't change dot-notation for function calls
        Expr(ex.head, ex.args[1], [makeDerVar(arg, parameters, inputs, FloatType, evaluateParameters) for arg in ex.args[2:end]]...)
    else
        Expr(ex.head, [makeDerVar(arg, parameters, inputs, FloatType, evaluateParameters) for arg in ex.args]...)
    end
end


text:
using NetworkEpidemics
using LightGraphs
using Statistics
using Plots
using ColorSchemes
using Random

Random.seed!(2021)

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SI}) = Metapopulation(mp.h, mp.D, SI(χ*mp.dynamics.β))

N = 1000
M = 10

h = complete_graph(M)
g = watts_strogatz(N, 100, 0.0)

k = degree(g)

χ = mean(k)/N # doesn't change so only compute it once

β = 0.1
D = 0.1

mp = Metapopulation(h, D, SI(β))
mpc = CorrectedMetapopulation(χ, mp)

mpx = Metaplex(g, h, D, SI(β))

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = rand(1:M, N)
x0_mp = zeros(Int, M, 2)
for i in 1:N
    x0_mp[x0_μ[i], x0_i[i]] += 1
end

tmax = 20.0
nmax = 10000
nsims = 200
nbins = 200

ts_mp, output_mp = gillespie(mp, x0_mp, tmax=tmax, nmax=nmax)
ts_mpx, output_mpx = gillespie(mpx, [x0_i,x0_μ], tmax=tmax, nmax=nmax)

plot(ts_mp, output_mp[2], label="")
plot(ts_mpx, output_mpx[2], label="")

ts = LinRange(0.0, tmax, nbins)
ps = LinRange(0.0, 0.4, 5)

ts_mf_mp, u_mf_mp = meanfield(mp, x0_mp, tmax=tmax, saveat=ts)
ts_mf_mpc, u_mf_mpc = meanfield(mpc, x0_mp, tmax=tmax, saveat=ts)

u_mf_mpx = zeros(length(ps), length(ts))
for i in 1:length(ps), j in 1:100
    g = watts_strogatz(N, 100, ps[i])
    mpx = Metaplex(g, h, D, SI(β))
    t, u = meanfield(mpx, [x0_i,x0_μ], tmax=tmax, saveat=ts)
    u_mf_mpx[i,:] += sum(u[2], dims=2)
end
u_mf_mpx /= 100

colors = ColorSchemes.GnBu_9;
orange = ColorSchemes.Oranges_3[3];

plot(ts, sum(u_mf_mpc[2], dims=2)/N,
    #label="Meanfield (No correction)",
    label = "Mean degree correction",
    xlabel="Time",
    ylabel="Fraction of infected individuals",
    fontfamily="Deja Vu Sans",
    legend=:bottomright,
    color = orange
    )
    plot!(ts, u_mf_mpx[1,:]/N,
    label = "p=0.0",
    color = colors[5]
    )
    plot!(ts, u_mf_mpx[2,:]/N,
    label = "p=0.1",
    color = colors[6]
    )
    plot!(ts, u_mf_mpx[3,:]/N,
    label = "p=0.2",
    color = colors[7]
    )
    plot!(ts, u_mf_mpx[4,:]/N,
    label = "p=0.3",
    color = colors[8]
    )
    plot!(ts, u_mf_mpx[5,:]/N,
    label = "p=0.4",
    color = colors[9]
    )
    plot!(ts, sum(u_mf_mpc[2], dims=2)/N,
    label="",
    color = orange
    )

plot(ts, (sum(u_mf_mpc[2], dims=2) - u_mf_mpx[1,:])/N,
    xlabel="Time",
    ylabel="Residual",
    label = "p=0.0",
    color = colors[5]
    )
    plot!(ts, (sum(u_mf_mpc[2], dims=2) - u_mf_mpx[2,:])/N,
    label = "p=0.1",
    color = colors[6]
    )
    plot!(ts, (sum(u_mf_mpc[2], dims=2) - u_mf_mpx[3,:])/N,
    label = "p=0.2",
    color = colors[7]
    )
    plot!(ts, (sum(u_mf_mpc[2], dims=2) - u_mf_mpx[4,:])/N,
    label = "p=0.3",
    color = colors[8]
    )
    plot!(ts, (sum(u_mf_mpc[2], dims=2) - u_mf_mpx[5,:])/N,
    label = "p=0.4",
    color = colors[9]
    )

cols = map(x->[x.r, x.g, x.b], colors.colors)
or = [orange.r, orange.g, orange.b]

using DelimitedFiles

writedlm("umpc.dat", sum(u_mf_mpc[2], dims=2)/N)
writedlm("umpx.dat", u_mf_mpx/N)
writedlm("ts.dat", ts)
writedlm("colors.dat", [cols..., or])

using MATLAB;

mat"""
figure
xlabel('Time')
ylabel('Fraction of Infected Individuals')
hold on
plot($ts, $(sum(u_mf_mpc[2], dims=2)/N), 'Color', $(or))
plot($ts, $(sum(u_mf_mpx[1], dims=2)/N), 'Color', $(cols[5]))
plot($ts, $(sum(u_mf_mpx[2], dims=2)/N), 'Color', $(cols[6]))
plot($ts, $(sum(u_mf_mpx[3], dims=2)/N), 'Color', $(cols[7]))
plot($ts, $(sum(u_mf_mpx[4], dims=2)/N), 'Color', $(cols[8]))
plot($ts, $(sum(u_mf_mpx[5], dims=2)/N), 'Color', $(cols[9]))
hold off
legend('Mean degree correction', 'p=0.0', 'p=0.1', 'p=0.2', 'p=0.3', 'p=0.4')
""";

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
using Test
using Empirikos

int1 = Interval(nothing, 0.0)
ebs1 = StandardNormalSample(int1)

@test Empirikos.likelihood(ebs1, 0.0) == 0.5

int2 = Interval{Open, Unbounded}(0.0, nothing)
ebs2 = StandardNormalSample(int2)

@test Empirikos.likelihood(ebs2, 0.0) == 0.5

vec_typed = StandardNormalSample.([int1; int1])
vec_not_typed = StandardNormalSample.([int1; int2])
vec_union_typed = StandardNormalSample{EBInterval{Float64}}.([int1;int2])

int3 = Interval{Closed, Closed}(1.0,2.0)
eb3 = NormalSample{typeof(int3),Float64}(int3, 1.0)
@inferred Empirikos.likelihood(eb3, 0.0)

Empirikos.likelihood.(vec_not_typed, 0.0)

@inferred broadcast(Empirikos.likelihood, vec_typed, 0.0)
#@test_throws @inferred broadcast(Empirikos.likelihood, vec_not_typed, 0.0)
@inferred broadcast(Empirikos.likelihood, vec_union_typed, 0.0)

@test isa(vec_not_typed, AbstractVector{EB} where EB <: EBayesSample)
@test isa(vec_union_typed, AbstractVector{EB} where EB <: EBayesSample)

@code_warntype Empirikos.likelihood(ebs2, 0.0)

@test skedasticity(vec_typed) == Empirikos.Homoskedastic()
@test skedasticity(vec_not_typed) == Empirikos.Homoskedastic()

vec_heteroskedastic = [NormalSample(0.0, 2.0); NormalSample(0.0, 3.0)]
@test skedasticity(vec_heteroskedastic) == Empirikos.Heteroskedastic()

int = Interval{Open,Closed}(1.0,2.)

interval_normal = NormalSample{typeof(int), Float64}(int, 1.0)
Base.summarysize(interval_normal)
@inferred likelihood(interval_normal, 0.0)

# What happens with union type?
union_interval_normal = NormalSample(int, 1.0)
Base.summarysize(union_interval_normal)
@inferred likelihood(union_interval_normal, 0.0)

any_normal = NormalSample{Any, Float64}(int, 1.0);
Base.summarysize(any_normal)
#@inferred likelihood(any_normal, 0.0)

ints = [Interval{Open,Closed}(1.0,2.0);Interval{Unbounded,Closed}(nothing,2.0);Interval{Open,Unbounded}(1.0,nothing)]
multinormals = NormalSample.(ints, 2.0)
untyped_multinormals = [NormalSample{typeof(int), Float64}(int, 2.0) for int in ints]
Base.summarysize(multinormals)
Base.summarysize(untyped_multinormals)

# Test some discrete Interval functionality
tmp = BinomialSample(Interval(0,1), 10)
@test cdf(Binomial(10, 0.5), 1) == likelihood(tmp, 0.5)
@test likelihood(BinomialSample(Interval(0,10), 10), 0.6) == 1

@test pdf(Dirac(0.7), BinomialSample(Interval(9,nothing), 10))  ≈ ccdf(Binomial(10,0.7), 8)
@test pdf(Dirac(0.7), BinomialSample(Interval(9,nothing), 10))  ≈  pdf(Dirac(0.7), BinomialSample(Interval(9,10), 10))

#pdf(DiscreteNonParametric([2.0],[1.0]), PoissonSample(Interval(5,nothing)))
#ccdf(Poisson(2.0), 4)


text:
"""
$(TYPEDEF)

A nonlinear system of equations.

# Fields
$(FIELDS)

# Examples

```
@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])
```
"""
struct NonlinearSystem <: AbstractSystem
    """Vector of equations defining the system."""
    eqs::Vector{Equation}
    """Unknown variables."""
    states::Vector
    """Parameters."""
    ps::Vector
    observed::Vector{Equation}
    """
    Name: the name of the system
    """
    name::Symbol
    """
    systems: The internal systems
    """
    systems::Vector{NonlinearSystem}
    """
    default_u0: The default initial conditions to use when initial conditions
    are not supplied in `ODEProblem`.
    """
    default_u0::Dict
    """
    default_p: The default parameters to use when parameters are not supplied
    in `ODEProblem`.
    """
    default_p::Dict
    """
    structure: structural information of the system
    """
    structure::Any
end

function NonlinearSystem(eqs, states, ps;
                         observed = [],
                         name = gensym(:NonlinearSystem),
                         default_u0=Dict(),
                         default_p=Dict(),
                         systems = NonlinearSystem[])
    default_u0 isa Dict || (default_u0 = Dict(default_u0))
    default_p isa Dict || (default_p = Dict(default_p))
    default_u0 = Dict(value(k) => value(default_u0[k]) for k in keys(default_u0))
    default_p = Dict(value(k) => value(default_p[k]) for k in keys(default_p))
    NonlinearSystem(eqs, value.(states), value.(ps), observed, name, systems, default_u0, default_p, nothing)
end

function calculate_jacobian(sys::NonlinearSystem;sparse=false,simplify=false)
    rhs = [eq.rhs for eq ∈ equations(sys)]
    vals = [dv for dv in states(sys)]
    if sparse
        jac = sparsejacobian(rhs, vals, simplify=simplify)
    else
        jac = jacobian(rhs, vals, simplify=simplify)
    end
    return jac
end

function generate_jacobian(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys);
                           sparse = false, simplify=false, kwargs...)
    jac = calculate_jacobian(sys,sparse=sparse, simplify=simplify)
    return build_function(jac, vs, ps;
                          conv = AbstractSysToExpr(sys), kwargs...)
end

function generate_function(sys::NonlinearSystem, dvs = states(sys), ps = parameters(sys); kwargs...)
    #obsvars = map(eq->eq.lhs, observed(sys))
    #fulldvs = [dvs; obsvars]
    fulldvs = dvs
    fulldvs′ = makesym.(value.(fulldvs))

    sub = Dict(fulldvs .=> fulldvs′)
    # substitute x(t) by just x
    rhss = [substitute(deq.rhs, sub) for deq ∈ equations(sys)]
    #obss = [makesym(value(eq.lhs)) ~ substitute(eq.rhs, sub) for eq ∈ observed(sys)]
    #rhss = Let(obss, rhss)

    dvs′ = fulldvs′[1:length(dvs)]
    ps′ = makesym.(value.(ps), states=())
    return build_function(rhss, dvs′, ps′;
                          conv = AbstractSysToExpr(sys), kwargs...)
end

jacobian_sparsity(sys::NonlinearSystem) =
    jacobian_sparsity([eq.rhs for eq ∈ equations(sys)],
                      states(sys))

function DiffEqBase.NonlinearFunction(sys::NonlinearSystem, args...; kwargs...)
    NonlinearFunction{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create an `NonlinearFunction` from the [`NonlinearSystem`](@ref). The arguments
`dvs` and `ps` are used to set the order of the dependent variable and parameter
vectors, respectively.
"""
function DiffEqBase.NonlinearFunction{iip}(sys::NonlinearSystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing,
                                     jac = false,
                                     eval_expression = true,
                                     sparse = false, simplify=false,
                                     kwargs...) where {iip}

    f_gen = generate_function(sys, dvs, ps; expression=Val{eval_expression}, kwargs...)
    f_oop,f_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in f_gen) : f_gen
    f(u,p) = f_oop(u,p)
    f(du,u,p) = f_iip(du,u,p)

    if jac
        jac_gen = generate_jacobian(sys, dvs, ps;
                                    simplify=simplify, sparse = sparse,
                                    expression=Val{eval_expression}, kwargs...)
        jac_oop,jac_iip = eval_expression ? (@RuntimeGeneratedFunction(ex) for ex in jac_gen) : jac_gen
        _jac(u,p) = jac_oop(u,p)
        _jac(J,u,p) = jac_iip(J,u,p)
    else
        _jac = nothing
    end

    NonlinearFunction{iip}(f,
                     jac = _jac === nothing ? nothing : _jac,
                     jac_prototype = sparse ? similar(sys.jac[],Float64) : nothing,
                     syms = Symbol.(states(sys)))
end

"""
```julia
function DiffEqBase.NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = states(sys),
                                     ps = parameters(sys);
                                     version = nothing,
                                     jac = false,
                                     sparse = false,
                                     kwargs...) where {iip}
```

Create a Julia expression for an `ODEFunction` from the [`ODESystem`](@ref).
The arguments `dvs` and `ps` are used to set the order of the dependent
variable and parameter vectors, respectively.
"""
struct NonlinearFunctionExpr{iip} end

function NonlinearFunctionExpr{iip}(sys::NonlinearSystem, dvs = states(sys),
                                     ps = parameters(sys), u0 = nothing;
                                     version = nothing, tgrad=false,
                                     jac = false,
                                     linenumbers = false,
                                     sparse = false, simplify=false,
                                     kwargs...) where {iip}

    idx = iip ? 2 : 1
    f = generate_function(sys, dvs, ps; expression=Val{true}, kwargs...)[idx]

    if jac
        _jac = generate_jacobian(sys, dvs, ps;
                                 sparse=sparse, simplify=simplify,
                                 expression=Val{true}, kwargs...)[idx]
    else
        _jac = :nothing
    end

    jp_expr = sparse ? :(similar($(sys.jac[]),Float64)) : :nothing

    ex = quote
        f = $f
        jac = $_jac
        NonlinearFunction{$iip}(f,
                         jac = jac,
                         jac_prototype = $jp_expr,
                         syms = $(Symbol.(states(sys))))
    end
    !linenumbers ? striplines(ex) : ex
end

function process_NonlinearProblem(constructor, sys::NonlinearSystem,u0map,parammap;
                           version = nothing,
                           jac = false,
                           checkbounds = false, sparse = false,
                           simplify=false,
                           linenumbers = true, parallel=SerialForm(),
                           eval_expression = true,
                           kwargs...)
    dvs = states(sys)
    ps = parameters(sys)
    u0map′ = lower_mapnames(u0map)
    u0 = varmap_to_vars(u0map′,dvs; defaults=default_u0(sys))
    defp = default_p(sys)

    if !(parammap isa DiffEqBase.NullParameters)
        parammap′ = lower_mapnames(parammap)
        p = varmap_to_vars(parammap′,ps; defaults=defp)
    elseif !isempty(defp)
        p = varmap_to_vars(Dict(),ps; defaults=defp)
    else
        p = ps
    end

    f = constructor(sys,dvs,ps,u0;jac=jac,checkbounds=checkbounds,
                    linenumbers=linenumbers,parallel=parallel,simplify=simplify,
                    sparse=sparse,eval_expression=eval_expression,kwargs...)
    return f, u0, p
end

function DiffEqBase.NonlinearProblem(sys::NonlinearSystem, args...; kwargs...)
    NonlinearProblem{true}(sys, args...; kwargs...)
end

"""
```julia
function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,
                                          parammap=DiffEqBase.NullParameters();
                                          jac = false, sparse=false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
```

Generates an NonlinearProblem from a NonlinearSystem and allows for automatically
symbolically calculating numerical enhancements.
"""
function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,
                                    parammap=DiffEqBase.NullParameters();kwargs...) where iip
    f, u0, p = process_NonlinearProblem(NonlinearFunction{iip}, sys, u0map, parammap; kwargs...)
    NonlinearProblem{iip}(f,u0,p;kwargs...)
end

"""
```julia
function DiffEqBase.NonlinearProblemExpr{iip}(sys::NonlinearSystem,u0map,
                                          parammap=DiffEqBase.NullParameters();
                                          jac = false, sparse=false,
                                          checkbounds = false,
                                          linenumbers = true, parallel=SerialForm(),
                                          kwargs...) where iip
```

Generates a Julia expression for a NonlinearProblem from a
NonlinearSystem and allows for automatically symbolically calculating
numerical enhancements.
"""
struct NonlinearProblemExpr{iip} end

function NonlinearProblemExpr(sys::NonlinearSystem, args...; kwargs...)
    NonlinearProblemExpr{true}(sys, args...; kwargs...)
end

function NonlinearProblemExpr{iip}(sys::NonlinearSystem,u0map,
                             parammap=DiffEqBase.NullParameters();
                             kwargs...) where iip

    f, u0, p = process_NonlinearProblem(NonlinearFunctionExpr{iip}, sys, u0map, parammap; kwargs...)
    linenumbers = get(kwargs, :linenumbers, true)

    ex = quote
        f = $f
        u0 = $u0
        p = $p
        NonlinearProblem(f,u0,p;$(kwargs...))
    end
    !linenumbers ? striplines(ex) : ex
end

function flatten(sys::NonlinearSystem)
    systems = get_systems(sys)
    if isempty(systems)
        return sys
    else
        return NonlinearSystem(
                               equations(sys),
                               states(sys),
                               parameters(sys),
                               observed=observed(sys),
                               default_u0=default_u0(sys),
                               default_p=default_p(sys),
                              )
    end
end


text:
using Test
using NQCBase

atoms = Atoms([:H, :C, :O, :N])
cell = PeriodicCell(rand(3, 3) .* 10)
R = rand(3, 4) .* 10

@testset "to/from_extxyz_dict" begin
    dict = NQCBase.to_extxyz_dict(atoms, R, cell)
    @test dict["cell"] ≈ au_to_ang.(permutedims(cell.vectors, (2,1)))
    natoms, nR, ncell = NQCBase.from_extxyz_dict(dict)
    @test ncell.vectors ≈ cell.vectors
    @test ncell.inverse ≈ cell.inverse
end

@testset "Single frame" begin
    write_extxyz("output.xyz", atoms, R, cell)
    new_atoms, new_R, new_cell = read_extxyz("output.xyz")
    @test new_atoms == atoms
    @test new_cell.vectors ≈ cell.vectors
    @test new_cell.inverse ≈ cell.inverse
    @test new_cell.periodicity ≈ cell.periodicity
    @test new_R ≈ [R]
end

@testset "Multiple frames" begin
    atoms = Atoms([:H, :C, :O, :N])
    cell = PeriodicCell(rand(3, 3) .* 10)
    R = [rand(3, 4) .* 10 for _=1:100]
    write_extxyz("output.xyz", atoms, R, cell)
    new_atoms, new_R, new_cell = read_extxyz("output.xyz")
    @test new_atoms == atoms
    @test new_cell.vectors ≈ cell.vectors
    @test new_cell.inverse ≈ cell.inverse
    @test new_cell.periodicity ≈ cell.periodicity
    @test new_R ≈ R
end

rm("output.xyz")


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -7.107954588851326 <= q <= 8.439913276188518)

                     add_NL_constraint(m, :(log(1 + exp(-0.4953615570310643 + -0.5197659840811282*log(1 + exp(0.7941458217392698 + 0.49653685598476827*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.8533482485058403*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + -0.13525586788633204*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))) + -0.07122185861106178*log(1 + exp(-0.06420985369672083 + 0.9131512287471262*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.03263883937508805*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + 0.0015134736925910275*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))) + 0.004105284291951339*log(1 + exp(-0.8124442669202567 + 0.4437121866105662*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.26832550078458484*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + 0.2804882536336004*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))))) + log(1 + exp(0.7344840166166842 + -0.5219551557665656*log(1 + exp(0.7941458217392698 + 0.49653685598476827*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.8533482485058403*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + -0.13525586788633204*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))) + -0.12237894168399821*log(1 + exp(-0.06420985369672083 + 0.9131512287471262*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.03263883937508805*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + 0.0015134736925910275*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))) + 0.5945370485807278*log(1 + exp(-0.8124442669202567 + 0.4437121866105662*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.26832550078458484*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + 0.2804882536336004*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))))) + log(1 + exp(-0.7546740103129919 + -0.5743271753296177*log(1 + exp(0.7941458217392698 + 0.49653685598476827*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.8533482485058403*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + -0.13525586788633204*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))) + 0.959961363428667*log(1 + exp(-0.06420985369672083 + 0.9131512287471262*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.03263883937508805*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + 0.0015134736925910275*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))) + -0.0685664132974031*log(1 + exp(-0.8124442669202567 + 0.4437121866105662*log(1 + exp(-0.2081348932235647 + -0.9819593709419507*$(x[1]) + 0.4754911924134797*$(x[2]) + -0.6006888917392823*$(x[3]))) + 0.26832550078458484*log(1 + exp(0.14812040410194793 + -0.8433474401842953*$(x[1]) + -0.1761112323855083*$(x[2]) + 0.6144836466445156*$(x[3]))) + 0.2804882536336004*log(1 + exp(-0.3937152244783322 + 0.36708048782945024*$(x[1]) + 0.138776295405171*$(x[2]) + -0.8492811105063662*$(x[3]))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions