text:
function variable_frequency_stab(pm::_PM.AbstractPowerModel; kwargs...)
    variable_frequency_reserves(pm; kwargs...)
    variable_frequency_reserves_bin(pm; kwargs...)
    variable_frequency_power_dev(pm; kwargs...)
    variable_frequency_power_dev_aux(pm; kwargs...)
end

function variable_frequency_stab_ne(pm::_PM.AbstractPowerModel; kwargs...)
    variable_frequency_reserves_ne(pm; kwargs...)
    variable_frequency_reserves_bin_ne(pm; kwargs...)
    variable_frequency_power_dev_ne(pm; kwargs...)
    variable_frequency_power_dev_aux_ne(pm; kwargs...)
end

function variable_frequency_reserves(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
    reserves = _PM.ref(pm, nw, :reserves)
    Pg = _PM.var(pm, nw)[:Pgg] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_pgg",
    lower_bound = 0, upper_bound = reserves[bi_bp[i]]["Pgmax"], start = 0)

    Pf = _PM.var(pm, nw)[:Pff] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_pff",
    lower_bound = 0, upper_bound = reserves[bi_bp[i]]["Pfmax"], start = 0)

    report && _IM.sol_component_value(pm, nw, :reserves, :Pgg, _PM.ids(pm, nw, :arcs_reserves_syn), Pg)
    report && _IM.sol_component_value(pm, nw, :reserves, :Pff, _PM.ids(pm, nw, :arcs_reserves_syn), Pf)
end

function variable_frequency_reserves_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    bi_bp = Dict([((reserves["syncarea"]),i ) for (i,reserves) in _PM.ref(pm, nw, :reserves)])
    reserves = _PM.ref(pm, nw, :reserves)
    Pg = _PM.var(pm, nw)[:Pgg] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_pgg",
    lower_bound = 0, upper_bound = reserves[bi_bp[i]]["Pgmax"], start = 0)

    Pf = _PM.var(pm, nw)[:Pff] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_pff",
    lower_bound = 0, upper_bound = reserves[bi_bp[i]]["Pfmax"], start = 0)

    report && _IM.sol_component_value(pm, nw, :reserves, :Pgg, _PM.ids(pm, nw, :arcs_reserves_syn_ne), Pg)
    report && _IM.sol_component_value(pm, nw, :reserves, :Pff, _PM.ids(pm, nw, :arcs_reserves_syn_ne), Pf)
end

function variable_frequency_reserves_bin(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    z1 = _PM.var(pm, nw)[:z1] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z1", binary = true, start = 0 )
    z2 = _PM.var(pm, nw)[:z2] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z2", binary = true, start = 0 )
    z3 = _PM.var(pm, nw)[:z3] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z3", binary = true, start = 0 )
    z4 = _PM.var(pm, nw)[:z4] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z4", binary = true, start = 0 )
    z5 = _PM.var(pm, nw)[:z5] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z5", binary = true, start = 0 )

    ev_syncarea = pm.setting["syncarea"]
    zb1 = _PM.var(pm, nw)[:zb1] = JuMP.@variable(pm.model,
    [i in _PM.ref(pm, nw, :bus_arcs_conv, ev_syncarea)], base_name="$(nw)_zb1", binary = true, start = 0 )

    report && _IM.sol_component_value(pm, nw, :reserves, :z1, _PM.ids(pm, nw, :arcs_reserves_syn), z1)
    report && _IM.sol_component_value(pm, nw, :reserves, :z2, _PM.ids(pm, nw, :arcs_reserves_syn), z2)
    report && _IM.sol_component_value(pm, nw, :reserves, :z3, _PM.ids(pm, nw, :arcs_reserves_syn), z3)
    report && _IM.sol_component_value(pm, nw, :reserves, :z4, _PM.ids(pm, nw, :arcs_reserves_syn), z4)
    report && _IM.sol_component_value(pm, nw, :reserves, :z5, _PM.ids(pm, nw, :arcs_reserves_syn), z5)
    report && _IM.sol_component_value(pm, nw, :convdc, :zb1, _PM.ref(pm, nw, :bus_arcs_conv, ev_syncarea), zb1)

    z11 = _PM.var(pm, nw)[:z11] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z11", binary = true, start = 0 )
    z21 = _PM.var(pm, nw)[:z21] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z21", binary = true, start = 0 )
    z31 = _PM.var(pm, nw)[:z31] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z31", binary = true, start = 0 )
    z41 = _PM.var(pm, nw)[:z41] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z41", binary = true, start = 0 )
    z51 = _PM.var(pm, nw)[:z51] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z51", binary = true, start = 0 )
    z12 = _PM.var(pm, nw)[:z12] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z12", binary = true, start = 0 )
    z22 = _PM.var(pm, nw)[:z22] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z22", binary = true, start = 0 )
    z32 = _PM.var(pm, nw)[:z32] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z32", binary = true, start = 0 )
    z42 = _PM.var(pm, nw)[:z42] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z42", binary = true, start = 0 )
    z52 = _PM.var(pm, nw)[:z52] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_z52", binary = true, start = 0 )

    report && _IM.sol_component_value(pm, nw, :reserves, :z11, _PM.ids(pm, nw, :arcs_reserves_syn), z11)
    report && _IM.sol_component_value(pm, nw, :reserves, :z21, _PM.ids(pm, nw, :arcs_reserves_syn), z21)
    report && _IM.sol_component_value(pm, nw, :reserves, :z31, _PM.ids(pm, nw, :arcs_reserves_syn), z31)
    report && _IM.sol_component_value(pm, nw, :reserves, :z41, _PM.ids(pm, nw, :arcs_reserves_syn), z41)
    report && _IM.sol_component_value(pm, nw, :reserves, :z51, _PM.ids(pm, nw, :arcs_reserves_syn), z51)

    report && _IM.sol_component_value(pm, nw, :reserves, :z12, _PM.ids(pm, nw, :arcs_reserves_syn), z12)
    report && _IM.sol_component_value(pm, nw, :reserves, :z22, _PM.ids(pm, nw, :arcs_reserves_syn), z22)
    report && _IM.sol_component_value(pm, nw, :reserves, :z32, _PM.ids(pm, nw, :arcs_reserves_syn), z32)
    report && _IM.sol_component_value(pm, nw, :reserves, :z42, _PM.ids(pm, nw, :arcs_reserves_syn), z42)
    report && _IM.sol_component_value(pm, nw, :reserves, :z52, _PM.ids(pm, nw, :arcs_reserves_syn), z52)

end

function variable_frequency_reserves_bin_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    z1 = _PM.var(pm, nw)[:z1] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z1", binary = true, start = 0 )
    z2 = _PM.var(pm, nw)[:z2] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z2", binary = true, start = 0 )
    z3 = _PM.var(pm, nw)[:z3] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z3", binary = true, start = 0 )
    z4 = _PM.var(pm, nw)[:z4] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z4", binary = true, start = 0 )
    z5 = _PM.var(pm, nw)[:z5] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z5", binary = true, start = 0 )
    ev_syncarea = pm.setting["syncarea"]
    zb1 = _PM.var(pm, nw)[:zb1] = JuMP.@variable(pm.model,
    [i in _PM.ref(pm, nw, :bus_arcs_conv_ne, ev_syncarea)], base_name="$(nw)_zb1", binary = true, start = 0 )

    report && _IM.sol_component_value(pm, nw, :reserves, :z1, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z1)
    report && _IM.sol_component_value(pm, nw, :reserves, :z2, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z2)
    report && _IM.sol_component_value(pm, nw, :reserves, :z3, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z3)
    report && _IM.sol_component_value(pm, nw, :reserves, :z4, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z4)
    report && _IM.sol_component_value(pm, nw, :convdc_ne, :zb1, _PM.ref(pm, nw, :bus_arcs_conv_ne, ev_syncarea), zb1)

    z11 = _PM.var(pm, nw)[:z11] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z11", binary = true, start = 0 )
    z21 = _PM.var(pm, nw)[:z21] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z21", binary = true, start = 0 )
    z31 = _PM.var(pm, nw)[:z31] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z31", binary = true, start = 0 )
    z41 = _PM.var(pm, nw)[:z41] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z41", binary = true, start = 0 )

    z12 = _PM.var(pm, nw)[:z12] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z12", binary = true, start = 0 )
    z22 = _PM.var(pm, nw)[:z22] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z22", binary = true, start = 0 )
    z32 = _PM.var(pm, nw)[:z32] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z32", binary = true, start = 0 )
    z42 = _PM.var(pm, nw)[:z42] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_z42", binary = true, start = 0 )

    report && _IM.sol_component_value(pm, nw, :reserves, :z11, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z11)
    report && _IM.sol_component_value(pm, nw, :reserves, :z21, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z21)
    report && _IM.sol_component_value(pm, nw, :reserves, :z31, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z31)
    report && _IM.sol_component_value(pm, nw, :reserves, :z41, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z41)

    report && _IM.sol_component_value(pm, nw, :reserves, :z12, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z12)
    report && _IM.sol_component_value(pm, nw, :reserves, :z22, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z22)
    report && _IM.sol_component_value(pm, nw, :reserves, :z32, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z32)
    report && _IM.sol_component_value(pm, nw, :reserves, :z42, _PM.ids(pm, nw, :arcs_reserves_syn_ne), z42)
end

function variable_frequency_power_dev(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)

    Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_Phvdcoaux", lower_bound = 0, start = 0)

    Phvdccaux = _PM.var(pm, nw)[:Phvdccaux] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_Phvdccaux", lower_bound = 0, start = 0)

    report && _IM.sol_component_value(pm, nw, :reserves, :Phvdcoaux, _PM.ids(pm, nw, :arcs_reserves_syn), Phvdcoaux)
    report && _IM.sol_component_value(pm, nw, :reserves, :Phvdccaux, _PM.ids(pm, nw, :arcs_reserves_syn), Phvdccaux)

end

function variable_frequency_power_dev_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)

    Phvdcoaux = _PM.var(pm, nw)[:Phvdcoaux] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_Phvdcoaux", lower_bound = 0, start = 0)

    Phvdccaux = _PM.var(pm, nw)[:Phvdccaux] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_Phvdccaux", lower_bound = 0, start = 0)

    Phvdcoaux_dup = _PM.var(pm, nw)[:Phvdcoaux_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_Phvdcoaux_dup", lower_bound = 0, start = 0)

    Phvdccaux_dup = _PM.var(pm, nw)[:Phvdccaux_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_Phvdccaux_dup", lower_bound = 0, start = 0)

    report && _IM.sol_component_value(pm, nw, :reserves, :Phvdcoaux, _PM.ids(pm, nw, :arcs_reserves_syn_ne), Phvdcoaux)
    report && _IM.sol_component_value(pm, nw, :reserves, :Phvdccaux, _PM.ids(pm, nw, :arcs_reserves_syn_ne), Phvdccaux)
    report && _IM.sol_component_value(pm, nw, :reserves, :Phvdcoaux_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), Phvdcoaux_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :Phvdccaux_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), Phvdccaux_dup)

end

function variable_frequency_power_dev_aux(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    k11 = _PM.var(pm, nw)[:k11] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k11", start = 0)
    k12 = _PM.var(pm, nw)[:k12] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k12", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k11, _PM.ids(pm, nw, :arcs_reserves_syn), k11)
    report && _IM.sol_component_value(pm, nw, :reserves, :k12, _PM.ids(pm, nw, :arcs_reserves_syn), k12)

    k21 = _PM.var(pm, nw)[:k21] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k21", start = 0)
    k22 = _PM.var(pm, nw)[:k22] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k22", start = 0)
    k23 = _PM.var(pm, nw)[:k23] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k23", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k21, _PM.ids(pm, nw, :arcs_reserves_syn), k21)
    report && _IM.sol_component_value(pm, nw, :reserves, :k22, _PM.ids(pm, nw, :arcs_reserves_syn), k22)
    report && _IM.sol_component_value(pm, nw, :reserves, :k23, _PM.ids(pm, nw, :arcs_reserves_syn), k23)

    k31 = _PM.var(pm, nw)[:k31] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k31", start = 0)
    k32 = _PM.var(pm, nw)[:k32] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k32", start = 0)
    k33 = _PM.var(pm, nw)[:k33] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k33", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k31, _PM.ids(pm, nw, :arcs_reserves_syn), k31)
    report && _IM.sol_component_value(pm, nw, :reserves, :k32, _PM.ids(pm, nw, :arcs_reserves_syn), k32)
    report && _IM.sol_component_value(pm, nw, :reserves, :k33, _PM.ids(pm, nw, :arcs_reserves_syn), k33)

    k41 = _PM.var(pm, nw)[:k41] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k41", start = 0)
    k42 = _PM.var(pm, nw)[:k42] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k42", start = 0)
    k43 = _PM.var(pm, nw)[:k43] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k43", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k41, _PM.ids(pm, nw, :arcs_reserves_syn), k41)
    report && _IM.sol_component_value(pm, nw, :reserves, :k42, _PM.ids(pm, nw, :arcs_reserves_syn), k42)
    report && _IM.sol_component_value(pm, nw, :reserves, :k43, _PM.ids(pm, nw, :arcs_reserves_syn), k43)

    k11_dup = _PM.var(pm, nw)[:k11_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k11_dup", start = 0)
    k12_dup = _PM.var(pm, nw)[:k12_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k12_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k11_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k11_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k12_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k12_dup)

    k21_dup = _PM.var(pm, nw)[:k21_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k21_dup", start = 0)
    k22_dup = _PM.var(pm, nw)[:k22_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k22_dup", start = 0)
    k23_dup = _PM.var(pm, nw)[:k23_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k23_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k21_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k21_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k22_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k22_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k23_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k23_dup)

    k31_dup = _PM.var(pm, nw)[:k31_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k31_dup", start = 0)
    k32_dup = _PM.var(pm, nw)[:k32_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k32_dup", start = 0)
    k33_dup = _PM.var(pm, nw)[:k33_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k33_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k31_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k31_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k32_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k32_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k33_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k33_dup)

    k41_dup = _PM.var(pm, nw)[:k41_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k41_dup", start = 0)
    k42_dup = _PM.var(pm, nw)[:k42_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k42_dup", start = 0)
    k43_dup = _PM.var(pm, nw)[:k43_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn)], base_name="$(nw)_k43_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k41_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k41_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k42_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k42_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k43_dup, _PM.ids(pm, nw, :arcs_reserves_syn), k43_dup)

end

function variable_frequency_power_dev_aux_ne(pm::_PM.AbstractPowerModel; nw::Int=pm.cnw, bounded::Bool = true, report::Bool=true)
    k11 = _PM.var(pm, nw)[:k11] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k11", start = 0)
    k12 = _PM.var(pm, nw)[:k12] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k12", start = 0)

    k21 = _PM.var(pm, nw)[:k21] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k21", start = 0)
    k22 = _PM.var(pm, nw)[:k22] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k22", start = 0)
    k23 = _PM.var(pm, nw)[:k23] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k23", start = 0)

    k31 = _PM.var(pm, nw)[:k31] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k31", start = 0)
    k32 = _PM.var(pm, nw)[:k32] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k32", start = 0)
    k33 = _PM.var(pm, nw)[:k33] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k33", start = 0)

    k41 = _PM.var(pm, nw)[:k41] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k41", start = 0)
    k42 = _PM.var(pm, nw)[:k42] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k42", start = 0)
    k43 = _PM.var(pm, nw)[:k43] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k43", start = 0)

    k11_dup = _PM.var(pm, nw)[:k11_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k11_dup", start = 0)
    k12_dup = _PM.var(pm, nw)[:k12_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k12_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k11_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k11_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k12_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k12_dup)

    k21_dup = _PM.var(pm, nw)[:k21_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k21_dup", start = 0)
    k22_dup = _PM.var(pm, nw)[:k22_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k22_dup", start = 0)
    k23_dup = _PM.var(pm, nw)[:k23_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k23_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k21_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k21_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k22_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k22_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k23_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k23_dup)

    k31_dup = _PM.var(pm, nw)[:k31_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k31_dup", start = 0)
    k32_dup = _PM.var(pm, nw)[:k32_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k32_dup", start = 0)
    k33_dup = _PM.var(pm, nw)[:k33_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k33_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k31_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k31_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k32_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k32_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k33_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k33_dup)

    k41_dup = _PM.var(pm, nw)[:k41_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k41_dup", start = 0)
    k42_dup = _PM.var(pm, nw)[:k42_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k42_dup", start = 0)
    k43_dup = _PM.var(pm, nw)[:k43_dup] = JuMP.@variable(pm.model,
    [i in _PM.ids(pm, nw, :arcs_reserves_syn_ne)], base_name="$(nw)_k43_dup", start = 0)
    report && _IM.sol_component_value(pm, nw, :reserves, :k41_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k41_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k42_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k42_dup)
    report && _IM.sol_component_value(pm, nw, :reserves, :k43_dup, _PM.ids(pm, nw, :arcs_reserves_syn_ne), k43_dup)

end

function sol_component_value_mod(pm::_PM.AbstractPowerModel, n::Int, comp_name::Symbol, variables)
    _PM.sol(pm, n)[comp_name] = variables
end


text:
"""
    abstract type updateInfo
- abstract supertype for all types of weight update
- holds algorithm-dependent information about training procedure
(e.g. am I using gradient descent or Reinforce etc, and what do I need for them)
"""
abstract type updateInfo  end

"""
    gradientDependentInfo <: updateInfo
abstract supertype for training procedures that require gradient calculation
"""
abstract type gradientDependentInfo <: updateInfo  end

"""
    gradientFreeInfo <: updateInfo
abstract supertype for training procedures that are gradient free
"""
abstract type gradientFreeInfo <: updateInfo end
"""
    fixedNormFluctuatingNoisyGD(sn_ratio::Float32, norm::Float32, intrinsic_norm::Float32) <: gradientDependentInfo
update = k1*normalised_gradient + k2*normalised_systematic_noise + k3*normalised_intrinsic_noise
- sn_ratio is k1/k2
- norm is √(k1^2 + k2^2)
- intrinsic norm is norm(k3)
"""
struct fixedNormFluctuatingNoisyGD <: gradientDependentInfo
    sn_ratio::Float32
    norm::Float32
    intrinsic_norm::Float32
end

"""
    fixedNormNoisyGD <: fixedNormFluctuatingNoisyGD
special case of fixedNormFluctuatingNoisyGD with intrinsic_norm = 0
"""
fixedNormNoisyGD(sn_ratio, norm) = fixedNormFluctuatingNoisyGD(sn_ratio,norm,0)
# fixedNormNoisyGD(;sn_ratio, norm) = fixedNormFluctuatingNoisyGD(sn_ratio,norm,0) #allow keyword arguments
# fixedNormNoisyGD(;sn_ratio, norm, intrinsic_norm) = fixedNormFluctuatingNoisyGD(sn_ratio,norm, intrinsic_norm)

non_updateable(upd::updateInfo) = (:lossInfo, :trainingInfo, :updateInfo, :updated)

"""
return things that are required in the store by :name, and have their own update method (i.e. things with potential conflicts wrt double updating)
"""
function get_storage_dependencies(names::Array{Symbol,1})
    to_return = []
    for name in names
        push!(to_return, name)
        if name == :local_task_difficulty
            push!(to_return,[:test_loss, :test_grad, :test_grad_norm, :update_norm]...)
        elseif name == :test_grad
            push!(to_return,[:test_loss]...)
        elseif name == :training_test_grad_correlation
            push!(to_return, [:full_training_grad_norm, :test_grad_norm, :test_grad, :full_training_grad, :full_training_loss]...)
        elseif name == :full_training_grad
            push!(to_return, [:full_training_loss]...)
        end
    end
    to_return = unique(to_return)
    return (to_return...,)
end

"""
    pnas_update(γ, T, N)
- Holds information about update from PNAS paper: Fundamental bounds on learning performance.
- update = T(γ1*normalised_gradient + γ2*normalised_systematic_noise + γ3*√N/√T*normalised_intrinsic_noise)
- not used explicitly: converts to 'equivalent' fixedNormNoisyGD for training
"""
pnas_update(γ::Array{NN,1}, T::TT, N::UU) where NN<:Number where TT <: Number where UU <: Number = fixedNormFluctuatingNoisyGD(T*γ[1], T*γ[2], sqrt(N*T)*γ[3])
pnas_update(; γ::Array{NN,1}, T::TT, N::UU) where NN<:Number where TT <: Number where UU <: Number = fixedNormFluctuatingNoisyGD(T*γ[1], T*γ[2], sqrt(N*T)*γ[3])

# function pnas_update_iterable(;  γ1::Union{AbstractRange, Number}, γ2::Union{AbstractRange, Number}, γ3::Union{AbstractRange,Number},
#                         T::Union{AbstractRange, Number},
#                         N::Union{AbstractRange, Number,Array{TT,1}}, repeats::Integer=1
#                     ) where TT
#     return [pnas_update( [_γ1,_γ2,_γ3], _T, _N) for (_γ1,_γ2,_γ3,_T,_N,r) in Base.product(γ1, γ2, γ3, T, N,1:repeats)]
# end

cdw_update(c::Number, normalised_corr::Number) = fixedNormFluctuatingNoisyGD(
    (normalised_corr/(c*(1-normalised_corr^2))), c, 0.)

################################################################################################################################################

"""
Struct for vanilla gradient descent
"""

"""
    vanillaGD(learning_rate)
implements gradient descent on the training set
"""
struct vanillaGD <: gradientDependentInfo
    learning_rate::Float32
end

"""
    beelineLearning(step_size, intrinsic_noise_norm, beehive)
moves the weights directly in the direction of 'beehive' (a predefined set of weights)
meanwhile intrinsic gaussian noise with a specified total norm makes the weights fluctuates
"""
struct beelineLearning <: gradientFreeInfo
    step_size::Float32
    intrinsic_norm::Float32
end
# beelineLearning(s, i, teacher) = beelineLearning(s,i

ext: jl


text:
module FileSessionAdapter

using Sessions, Genie, Logger, Genie.Configuration, App

const SESSION_FOLDER = IS_IN_APP ? App.config.session_folder : tempdir()

"""
    write(session::Sessions.Session) :: Sessions.Session

Persists the `Session` object to the file system, using the configured sessions folder and returns it.
"""
function write(session::Sessions.Session) :: Sessions.Session
  try
    open(joinpath(SESSION_FOLDER, session.id), "w") do (io)
      serialize(io, session)
    end
  catch ex
    Logger.log("Error when serializing session $(escape_string(session)) in $(@__FILE__):$(@__LINE__)", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    rethrow(ex)
  end

  session
end

"""
    read(session_id::Union{String,Symbol}) :: Nullable{Sessions.Session}
    read(session::Sessions.Session) :: Nullable{Sessions.Session}

Attempts to read from file the session object serialized as `session_id`.
"""
function read(session_id::Union{String,Symbol}) :: Nullable{Sessions.Session}
  try
    isfile(joinpath(SESSION_FOLDER, session_id)) || return Nullable{Sessions.Session}(write(Session(session_id)))
  catch ex
    Logger.log("Can't check session file", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    Nullable{Sessions.Session}(write(Session(session_id)))
  end

  try
    session = open(joinpath(SESSION_FOLDER, session_id), "r") do (io)
      deserialize(io)
    end

    Nullable{Sessions.Session}(session)
  catch ex
    Logger.log("Can't read session", :err)
    Logger.log(string(ex), :err)
    Logger.log("$(@__FILE__):$(@__LINE__)", :err)

    Nullable{Sessions.Session}(write(Session(session_id)))
  end
end
function read(session::Sessions.Session) :: Nullable{Sessions.Session}
  read(session.id)
end

end


text:
@testset "value access" begin
    @test year(nd) == yr
    @test month(nd) == mn
    @test day(nd) == dy
    @test hour(nd) == hr
    @test minute(nd) == mi
    @test second(nd) == sc
    @test millisecond(nd) == ms
    @test microsecond(nd) == cs
    @test nanosecond(nd) == ns
end

@testset "period access" begin
    @test Year(nd) == Yr
    @test Month(nd) == Mn
    @test Day(nd) == Dy
    @test Hour(nd) == Hr
    @test Minute(nd) == Mi
    @test Second(nd) == Sc
    @test Millisecond(nd) == Ms
    @test Microsecond(nd) == Cs
    @test Nanosecond(nd) == Ns
end


text:
module epilepsy

using Statistics
using LinearAlgebra
using Distances

export sample_entropy, multiscale_entropy

# Write your package code here.
include("entropy/sample_entropy.jl")

end


text:
using IRTracker
using DynamicPPL

export trackmodel, strip_model_layers

# F is always the type of the model evaluator closure

# dispatching on values
ismodelcall(::Type{F}, f, args...) where {F} = ismodelcall(F, typeof(f), typeof.(args)...)
ismodelcall(::Type{F}, ::NestedCallNode{<:Any, G, TArgs}) where {F, G, TArgs} = ismodelcall(F, G, TArgs.parameters...)

# dispatching on types
ismodelcall(::Type{F}, ::Type{<:Model{F}}, ::Vararg{<:Type}) where {F} = true
ismodelcall(::Type{F}, ::Type{typeof(Core._apply)}, ::Type{<:Model{F}}, ::Vararg{<:Type}) where {F} = true
ismodelcall(::Type{F}, ::Type{typeof(DynamicPPL.evaluate_threadsafe)}, ::Type, ::Type{<:Model{F}}, ::Vararg{<:Type}) where {F} = true
ismodelcall(::Type{F}, ::Type{typeof(DynamicPPL.evaluate_threadunsafe)}, ::Type, ::Type{<:Model{F}}, ::Vararg{<:Type}) where {F} = true
ismodelcall(::Type{F}, ::Type{F}, ::Vararg{<:Type}) where {F} = true
ismodelcall(::Type{F}, ::Type{G}, ::Vararg{<:Type}) where {F, G} = false

struct AutoGibbsContext{F} <: AbstractTrackingContext end

IRTracker.canrecur(ctx::AutoGibbsContext{F}, f, args...) where {F} = ismodelcall(F, f, args...)

trackmodel(model::Model{F}, args...) where {F} = track(AutoGibbsContext{F}(), model, args...)

function strip_model_layers(::Type{F}, node::NestedCallNode{<:Any, F}) where {F}
    # base case: call of the model closure, F
    return node
end

function strip_model_layers(::Type{F}, node::NestedCallNode) where {F}
    if ismodelcall(F, node)
        children = getchildren(node)
        ix = findfirst(Base.Fix1(ismodelcall, F), children)

        if ix !== nothing
            return strip_model_layers(F, children[ix])
        else
            throw(ArgumentError("$node does not contain any model call, you should not have reached this point!"))
        end
    else
        throw(ArgumentError("$node is not a model call, you should not have reached this point!"))
    end
end


text:
# We need these as some functions have moved from Base to Stdlib
#
using Printf, SpecialFunctions, LinearAlgebra

x = 2;   typeof(x)

x = 2.0;  typeof(x)

x0 = 2^65
x1 = big(2)^65
@assert x0 == x1

for T = Any[Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128]
    println("$(lpad(T,7)): [$(typemin(T)),$(typemax(T))]")
end

x = 0xbb31; y = 0xaa5f;  xor(x,y)

x = 0xbb31;  x << 8

x = 0xbb31; p = (2 < 3); x + p

# v1.0 requires the Statistics module for mean(), std() etc.
#
using Statistics

# Mean of 15 random numbers in range 0:100
#
A = rand(0:100,15)
mean(A)

# Create an empty array, note new syntax
# 'undef' will not initialise the elements
#
A = Array{Int64,1}(undef, 15)

# Verify: Tuple of the element type and the dimension sizes
#
(eltype(A),size(A))

# Fill array A with the first 15 Fibonnaci series
#
A[1] = 1
A[2] = 1
[A[i] = A[i-1] + A[i-2] for i = 3:length(A)]

# The 'recursive' definition of factorial function
# A simple loop is much quicker
#
function fac(n::Integer)
  @assert n > 0
  (n == 1) ? 1 : n*fac(n-1)
end

# This has difficulties with integer overflow
# We now need the Printf module to use the @printf macro
#
using Printf
for i = 1:30
  @printf "%3d : %d \n" i fac(i)
end

# But since a BigInt <: Integer if we pass a BigInt the reoutine returns one
#
fac(big(30))

# Find stdlib, location is O/S dependent

cd(Sys.BINDIR)
cd("../share/julia/stdlib/v1.0")
pwd()

# We can check this using the gamma function
# Again we need a module (SpecialFunctions)
#
using SpecialFunctions
gamma(31)     # Γ(n+1)  <=>  n!

# This non-recursive one liner works!
# Note that this returns a BigInt regardless of the input
#
fac(N::Integer) =
  (N < 1) ? throw(ArgumentError("N must be positive")) : reduce(*,big.(collect(1:N)))

@time(fac(402))

gamma(big(403.0))

# The 'standard' recursive definition

function fib(k::Integer)
  @assert k > 0
  (k < 3) ? 1 : fib(k-1) + fib(k-2)
end

@time fib(15)

# A better version

function fib(n::Integer)
  @assert n > 0
  a = Array{typeof(n),1}(undef,n)
  a[1] = 0
  a[2] = 1
  for i = 3:n
    a[i] = a[i-1] + a[i-2]
  end
  return a[n]
end

@time(fib(big(402)))

# A still better version
# This requires no array storage

function fib(n::Integer)
  @assert n > 0
  (a, b) = (big(0), big(1))
  while n > 0
    (a, b) = (b, a+b)
    n -= 1
  end
  return a
end

# Golden ratio

@printf "%.15f" fib(101)/fib(100)

# Check with the actual value

γ = (1.0 + sqrt(5.0))/2.0

using Random # stdlib module needed for srand() => seed!()
#
tm = round(time());
seed = convert(Int64,tm);
Random.seed!(seed);

function bacs()
         bulls = cows = turns = 0
         a = Any[]
         while length(unique(a)) < 4
           push!(a,rand('0':'9'))
         end
         my_guess = unique(a)
         println("Bulls and Cows")
         while (bulls != 4)
           print("Guess? > ")
           s = chomp(readline(stdin))
           if (s == "q")
             print("My guess was "); [print(my_guess[i]) for i=1:4]
             return
           end
           guess = collect(s)
           if !(length(unique(guess)) == length(guess) == 4 && all(isdigit,guess))
             print("\nEnter four distinct digits or q to quit: ")
             continue
           end
           bulls = sum(map(==, guess, my_guess))
           cows = length(intersect(guess,my_guess)) - bulls
           println("$bulls bulls and $cows cows!")
           turns += 1
         end
         println("\nYou guessed my number in $turns turns.")
       end

bacs()

# The matrix file is in Files subdirectory
# Check we are at the correct location
#
cd() # Change as needed
pwd()

#http://en.wikipedia.org/wiki/Stochastic_matrix
#
# Create stochastic matrix and write to disk
#
open("./cm3.txt","w") do f
  write(f,"0.0,0.0,0.5,0.0\n")
  write(f,"0.0,0.0,1.0,0.0\n")
  write(f,"0.25,0.25,0.0,0.25\n")
  write(f,"0.0,0.0,0.5,0.0\n")
end

using DelimitedFiles

I = zeros(4,4);
[I[i,i] = 1 for i in 1:4];

f = open("./cm3.txt","r")
T = readdlm(f,',');
close(f);

T

Ep = [0 1 0 0]*inv(I - T)*[1,1,1,1];

println("Expected lifetime for the mouse is $(Ep[1]) hops.")

# Look at different definitions of the norm function
# For a Gaussian distribution of size N we should expect the answer ~= √N
# The first call f1(1) is to run in the function and not affect the timing
# This version uses the function in the stdlib LinearAlgebra module

using LinearAlgebra

f1(n) = norm(randn(n))
f1(10);
@time f1(100_000_000)

# We can get the same result using a mapreduce procedure
# Note that it is a new set of random number, so the answer is slightly different
# The time is about the same

f2(n) = sqrt(mapreduce(x -> x*x, +, randn(n)))
f2(10);
@time f2(100_000_000)

# Using a conventional mapping we need to pipe the result to sum it
# and then take the square root
# This takes a little longer tha the previous 2

f3(n) = map(x -> x*x,randn(n)) |> sum |> sqrt
f3(10);
@time f3(100_000_000)

# Finally we can non-vectorize the code, which is much quicker,
# In Julia non-vectorized (i.e loopy) code is invariably faster
# than the vectorized equivalent.

function f4(n)
    t = 0.0
    for i = 1:n
        t += randn()^2
    end
    return sqrt(t)
end

f4(10);
@time f4(100_000_000)

# Define a very simple type to represent a 2-D point

struct Point
  x::Real
  y::Real
end

# We can define how to add and scale points
# This needs importing the + and * function from Base

import Base: +,*
+(u::Point,v::Point) = Point(u.x + v.x, u.y + v.y )
*(a::Real,u::Point)  = Point(a*u.x, a*u.y)
*(u::Point,a::Real)  = a*u

# Just test the type structure

u1 = Point(1.0,2.3)
u1*(17//13)

# Using Julia's aggregate object model this type "knows" all about arrays
# Note: I'll deal with the object model in greater detail the next two chapters
#
aa = [Point(randn(),randn()) for i = 1:100];
ab = reshape(aa,10,10)

# It is useful to define a zero function (not sure about a one())

zero(Point) = Point(0.0,0.0)

# The dot product is the sum of the product of the x and y coordinates

dot(u::Point)::Real = u.x^2 * u.y^2
dot(u::Point,v::Point)::Real = u.x*v.x * u.y*v.y
dot(ab[1])

# The distance between two points is determined by Pythagoras's rule

dist(u::Point,v::Point)::Real = sqrt((u.x - v.x)^2 + (u.y - v.y)^2)
dist(ab[4,1],ab[2,7])

# The distance of the point from the origin is equivalent to it's norm
# We need to import the 'norm' function

import LinearAlgebra.norm

norm(u::Point)::Real = sqrt(u.x^2 + u.y^2)
norm(aa[17])

# It is also possible to define this as:

dist(u::Point)::Real = dist(u::Point,zero(Point))

# Although this requires slightly more work to compute and it may be
# better just to define it as dist(u::Point) = norm(u)

@assert(dist(aa[17]) == norm(aa[17]))  ## Should produce NO output if it is TRUE

typeof(aa)  # Note that this is an array of points
            # So norm will not work on this yet

# One way is to overload the norm function using the mapreduce function above
# (or write a non-vectorized one)

norm(a::Array{Point,1}) = sqrt(mapreduce(x -> dist(x), +, a))

# And now we can estimate PI (as in the previous chapter)

N = 1000000; ac = [Point(rand(),rand()) for i = 1:N];

# We need the let/end because of the scoping rules
let
  count = 0
  for i = 1:N
    (dist(ac[i]) < 1.0) && (count += 1)
  end
  4.0*(count/N)
end

# Note as yet norm(ab) will not work so we need to be a little more imaginative
# with the function definition.

# In the next chapter we will generalise this type definition to 3-Vector and N-vector
# and revisit the quiestion then

norm(ab)

function juliaset(z, z0, nmax::Int64)
    for n = 1:nmax
        if abs(z) > 2 (return n-1) end
        z = z^2 + z0
    end
    return nmax
end

function create_pgmfile(img, outf::String)
    s = open(outf, "w")
    write(s, "P5\n")
    n, m = size(img)
    write(s, "$m $n 255\n")
    for i=1:n, j=1:m
        p = min(img[i,j],255)
        write(s, UInt8(p))
    end
    close(s)
end

h = 400;
w = 800;
m = Array{Int64,2}(undef,h,w);

c0 = -0.8 + 0.16im;
pgm_name = "jset.pgm";

t0 = time();
for y=1:h, x=1:w
    c = complex((x-w/2)/(w/2), (y-h/2)/(w/2))
    m[y,x] = juliaset(c, c0, 256)
end
t1 = time();

# You should find the file in the same chapter as the notebook

create_pgmfile(m, pgm_name);
@printf "Written %s\nFinished in %.4f seconds.\n" pgm_name (t1-t0);

# Display the image using Imagemagick's display command
# Clicking the file may display it (OSX/XQuartz and Centos/Gnome certainly does)
# On OSX use:  brew reinstall imagemagick --with-x11
#
# Otherwise you make be able to click onit (OSX or Linux)
# or on WIndows use an image processing program such as Irfanview
# We will be looking at how to do this in Julia later.

here = pwd()
run(`display $here/$pgm_name`)


text:
#
# Manager of jevo plugins. Plugins - are additional modules, which
# listens different events and do different stuff like statistics,
# phylogenetic tree creation and so on. They compiled only if there
# is an option in Config module (see Config.ConfigData.plugIncluded).
# They also work in listen mode. It means that they shouldn't have
# public interface. Only event handlers. So we may disable them
# by changing configuration from command line. The meaning of configs
# plugIncluded and plugExcluded is very simple: plugIncluded contains
# default list of active plugins. Using plugExcluded you may disable
# (exclude) some of them. You may use "all" value in both settings to
# include or exclude all plugins. Example of excluding of Phylogen
# plugin:
#
# >julia src\AppManager.jl plugExcluded:["Phylogen"]
#
# @author DeadbraiN
#
import Config
#
# Means all plugins in include or exclude lists
#
const ALL_PLUGINS = "all"
#
# Is used for right padding with "." symbol
#
const PLUGIN_NAME_MAX_LEN = 20
#
# Initializes all available plugins according to configuration.
# Included and excluded plugins are set in Config.DataConfig.plugIncluded,
# plugExcluded
# @param cfg Configuration data object
# @param obs Observer for events
# @return {Array{String,1}} Activated plugins array
#
function initPlugins(man::ManagerTypes.ManagerData)
  local plugin::String
  local plugins::Array{String, 1}
  local included::Array{String, 1} = man.cfg.plugIncluded
  local excluded::Array{String, 1} = man.cfg.plugExcluded
  #
  # "all" setting is in exclude list. We have to disable all plugins
  # by returning nothing. require() will not be called for them
  #
  if findfirst(excluded, ALL_PLUGINS) > 0 return false end
  #
  # if "all" setting is in include list, we have to include all plugins.
  # otherwise excluded plugins should be excluded from included list
  #
  plugins = findfirst(included, ALL_PLUGINS) > 0 ? included : setdiff(included, excluded)
  if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info("Initializing plugins: ") end
  for plugin in plugins
    if plugin === ALL_PLUGINS continue end
    if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.info(string("  ", rpad(plugin, PLUGIN_NAME_MAX_LEN, ".")), false) end
    Base.require(Symbol(plugin))
    getfield(Main.eval(Symbol(plugin)), :init)(man)
    if man.cfg.modeQuiet < Config.MODE_QUIET_NO Helper.done() end
  end

  true
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
module CategoricalTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable, prod!, unsafeMean, unsafeVar, vague, dims, logPdf, naturalParams, standardDistribution
using ForneyLab: SPCategoricalOutNP, VBCategoricalOut, VBCategoricalIn1
using SparseArrays: SparseVector, spzeros

@testset "Categorical ProbabilityDistribution and Message construction" begin
    @test ProbabilityDistribution(Univariate, Categorical, p=[0.1, 0.8, 0.1]) == ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>[0.1, 0.8, 0.1]))
    @test_throws Exception ProbabilityDistribution(Multivariate, Categorical)
    @test ProbabilityDistribution(Categorical, p=[0.1, 0.8, 0.1]) == ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>[0.1, 0.8, 0.1]))
    @test ProbabilityDistribution(Categorical) == ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>[1/3, 1/3, 1/3]))
    @test Message(Categorical) == Message{Categorical, Univariate}(ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>[1/3, 1/3, 1/3])))
    @test Message(Univariate, Categorical) == Message{Categorical, Univariate}(ProbabilityDistribution{Univariate, Categorical}(Dict(:p=>[1/3, 1/3, 1/3])))
    @test_throws Exception Message(Multivariate, Categorical)
end

@testset "dims" begin
    @test dims(ProbabilityDistribution(Categorical, p=[0.1, 0.8, 0.1])) == ()
end

@testset "vague" begin
    @test vague(Categorical) == ProbabilityDistribution(Categorical, p=[1/3, 1/3, 1/3])
    @test vague(Categorical, 4) == ProbabilityDistribution(Categorical, p=[1/4, 1/4, 1/4, 1/4])
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(Categorical, p=[0.1, 0.8, 0.1])) == [0.1, 0.8, 0.1]
end

@testset "sample" begin
    s = sample(ProbabilityDistribution(Categorical, p=[0.1, 0.8, 0.1]))
    @test isa(s, SparseVector)
    @test sum(s.==1.0) == 1
    @test sum(s.==0.0) == 2
end

@testset "log pdf" begin
    x = spzeros(Float64, 3)
    x[3] = 1.0
    @test isapprox(logPdf(ProbabilityDistribution(Categorical, p=[0.2, 0.3, 0.5]),x), -0.6931471805599453)
end

@testset "prod!" begin
    @test ProbabilityDistribution(Categorical, p=[0.2, 0.8])*ProbabilityDistribution(Categorical, p=[0.8, 0.2]) == ProbabilityDistribution(Categorical, p=[0.5, 0.5])
    @test ProbabilityDistribution(Categorical, p=[0.25, 0.5, 0.25]) * ProbabilityDistribution(Categorical, p=[1/3, 1/3, 1/3]) == ProbabilityDistribution(Categorical, p=[0.25, 0.5, 0.25])
    @test ProbabilityDistribution(Categorical, p=[0.0, 0.5, 0.5]) * ProbabilityDistribution(Categorical, p=[1.0, 0.0, 0.0]) == ProbabilityDistribution(Categorical, p=ones(3)/3)
end

@testset "natural parameters" begin
    d = ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8])
    η = naturalParams(d)
    s = standardDistribution(Univariate, Categorical, η=η)
    @test d.params[:p] == s.params[:p] # Test conversion consistency

    x = [[0.0, 1.0], [1.0, 0.0]]
    d_x = logPdf.([d], x)
    η_x = logPdf.(Univariate, Categorical, x; η=η)
    @test isapprox(d_x, η_x) # Test pdf consistency
end

#-------------
# Update rules
#-------------

@testset "SPCategoricalOutNP" begin
    @test SPCategoricalOutNP <: SumProductRule{Categorical}
    @test outboundType(SPCategoricalOutNP) == Message{Categorical}
    @test isApplicable(SPCategoricalOutNP, [Nothing, Message{PointMass}])

    @test ruleSPCategoricalOutNP(nothing, Message(Multivariate, PointMass, m=[0.1, 0.8, 0.1])) == Message(Univariate, Categorical, p=[0.1, 0.8, 0.1])
end

@testset "VBCategoricalOut" begin
    @test VBCategoricalOut <: NaiveVariationalRule{Categorical}
    @test outboundType(VBCategoricalOut) == Message{Categorical}
    @test isApplicable(VBCategoricalOut, [Nothing, ProbabilityDistribution])
    @test !isApplicable(VBCategoricalOut, [ProbabilityDistribution, Nothing])

    @test ruleVBCategoricalOut(nothing, ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1])) == Message(Univariate, Categorical, p=[0.10000000000000002, 0.8, 0.10000000000000002])
    @test ruleVBCategoricalOut(nothing, ProbabilityDistribution(Multivariate, Dirichlet, a=[0.1, 0.8, 0.1])) == Message(Univariate, Categorical, p=[7.799215056092699e-5, 0.999844015698878, 7.799215056092699e-5])
end

@testset "VBCategoricalIn1" begin
    @test VBCategoricalIn1 <: NaiveVariationalRule{Categorical}
    @test outboundType(VBCategoricalIn1) == Message{Dirichlet}
    @test isApplicable(VBCategoricalIn1, [ProbabilityDistribution, Nothing])

    @test ruleVBCategoricalIn1(ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1]), nothing) == Message(Multivariate, Dirichlet, a=[1.1, 1.8, 1.1])
    @test ruleVBCategoricalIn1(ProbabilityDistribution(Univariate, Categorical, p=[0.1, 0.8, 0.1]), nothing) == Message(Multivariate, Dirichlet, a=[1.1, 1.8, 1.1])
    @test ruleVBCategoricalIn1(ProbabilityDistribution(Univariate, Bernoulli, p=0.1), nothing) == Message(Multivariate, Dirichlet, a=[1.1, 1.9])
end

@testset "averageEnergy and differentialEntropy" begin
    @test differentialEntropy(ProbabilityDistribution(Univariate, Categorical, p=[0.1, 0.8, 0.1])) == averageEnergy(Categorical, ProbabilityDistribution(Univariate, Categorical, p=[0.1, 0.8, 0.1]), ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.8, 0.1]))
    @test differentialEntropy(ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8])) == differentialEntropy(ProbabilityDistribution(Univariate, Bernoulli, p=0.2))
    @test averageEnergy(Categorical, ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), ProbabilityDistribution(Multivariate, PointMass, m=[0.1, 0.9])) == averageEnergy(Bernoulli, ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Univariate, PointMass, m=0.1))
end

end # module


text:
function macroexpand(o::OrgDoc, m::Macro)
    if m.name in ("title", "author", "email") &&
        haskey(o.keywords, m.name)
        join(o.keywords[m.name], ' ')
    elseif m.name == "keyword" && length(m.arguments) == 1 &&
        haskey(o.keywords, m.arguments[1])
        join(o.keywords[m.arguments[1]], ' ')
    elseif m.name == "date" && haskey(o.keywords, "date")
        if length(m.arguments) == 0
            o.keywords["date"][1]
        else
            date_timestamp = consume(Timestamp, o.keywords["date"][1])
            if !isnothing(date_timestamp)
                seconds = if !isnothing(date_timestamp[2].time)
                    DateTime(date_timestamp[2].date, date_timestamp[2].time)
                else
                    DateTime(date_timestamp[2].date)
                end |> datetime2unix
                Base.Libc.strftime(m.arguments[1], seconds)
            end
        end
    elseif m.name == "time" && length(m.arguments) > 0
        Base.Libc.strftime(m.arguments[1], time())
    # TODO support n macro ... somehow
    elseif m.name in keys(o.macros)
        o.macros[m.name](m.arguments)
    elseif m.name == "results" && length(m.arguments) == 1
        "=$(m.arguments[1])="
    end
end


text:
export Opt_SpectralRatio_SA
include("Spectral_GapR.jl")
using Statistics,  Random, Distributions
# return SG Ratio

###### neighbor perturb
function sample_neighbour(nR::Int64, rec::Vector, r::Float64, p::Int64)
    #rec: the subsampled locations
    #r: the ratio of subsampled points to be perturbed
    #p: subsampling factor

    #the number of all
    n = length(rec);

    #random choose the perturb locations
    Num = Int(round(length(rec)*r));
    A = sortperm(randn(length(rec)));
    Perturb_rec = rec[A[1:Num]];

    #constructed perturbed area based on selected locations and neightbor range
    test_range = Int(floor((Perturb_rec[1] -1) / p))
    Perturb_range = Array(collect(test_range*p+1:1:(test_range+1)*p)')
    for i = 2: Num
    	test_range = Int(floor((Perturb_rec[i] -1) / p))
    	Perturb_range = vcat(Perturb_range, Array(collect(test_range*p+1:1:(test_range+1)*p)'))
    end

    # range constraint especially for the edge
    Perturb_range[findall(x-> x.<=0, Perturb_range)] .= 1
    Perturb_range[findall(x-> x.>=nR, Perturb_range)] .= nR

    # replace range and remove repeated locations
    new_rec = copy(rec);
    new_rec[A[1:Num]] .= Perturb_range[CartesianIndex.(Tuple.(eachrow(hcat(1:Num,rand(1:5,Num)))))];

    return sort(unique(new_rec))
end

function Opt_SpectralRatio_SA(nR::Int64, locs1::Vector, pprob::Int64, r::Float64)
    #Input: nR: number of sources
    #       locs1: given locations (Uniform random or optimal jittered)
    #       p: subsampling factor
    #       r: the percentage of data to be removed and here we set r=20%
    #Output: best_locs: the subsampled locations with smaller spectral ratio
    #        best_sgr[1]: the initial SG ratio
    #        best_sgr[end]: the output SG ratio

	best_locs = locs1;

	# calculate the evaluation (Spectral ratio)
	cur_sgr,_= Spectral_GapR(nR, locs1)
	best_sgr = [cur_sgr];

    # set the probability used to control the moves
	function prob(n, res)
	    return exp(-res/(95f-2 * (70f-2)^n))
	end

	max_itr = 4000 #iteration numbers
	prob_hist = [];
	for itr in range(1, max_itr, step=1)
        # sample neighbour to get the new location and calculate the SG ratio
	    locs =  sample_neighbour(nR, best_locs, r, pprob);
	    cur_sgr, _ = Spectral_GapR(nR, locs);

        # if better, move to the new locations
        if cur_sgr < best_sgr[end]
            push!(best_sgr, cur_sgr)
            push!(prob_hist, 1)
            best_locs = locs
        else #else use the probability to decide move or not
            p = prob(itr, cur_sgr-best_sgr[end])
            if rand() < p
                best_locs = locs;
                push!(best_sgr, cur_sgr)
            end
	    end
	end

    # return the new locations, the input SG ratio and the output SG ratio
	return best_locs, best_sgr[1], best_sgr[end];
end


text:
abstract type AbstractConstruct{T} end
######################################
# Base functionalities for AbstractConstruct
######################################

"""
    eltype(::AbstractConstruct)

Returns the type that is generated by the iterator on an `AbstractConstruct` structure
"""
Base.eltype(::AbstractConstruct{T}) where {T} = T

"""
    getindex(Construct::AbstractConstruct,i::Int)

Returns the module on the i-th possition in the construct
Output of type `Mod`
"""
Base.getindex(c::AbstractConstruct,i::Int) = c.c[i]

"""
    getindex(c::AbstractConstruct,i::UnitRange)

Allows to index an `AbstractConstruct` using an `UnitRange` input.
Return a sub-construct of `c`,
base on the indexes defined on in the `UnitRange` object.
"""
Base.getindex(c::AbstractConstruct,i::UnitRange) = c.c[i] |> typeof(c)

"""
    length(Construct::AbstractConstruct)

Returns the length of the input construct, i.e., the number of modules that are
used.
"""
Base.length(c::AbstractConstruct) = length(c.c)

"""
    lastindex(c::AbstractConstruct)

Returns the last module of a construct, allow the uses of `end` when indexing on a `c`.
"""
Base.lastindex(c::AbstractConstruct) =  last(eachindex(IndexLinear(), c.c))

"""
    iterate(Construct::AbstractConstruct,state=1)

Iterator to loop over all modules in a given construct.
"""
function Base.iterate(Construct::AbstractConstruct,state=1)
    if  state <= length(Construct)
        mod = Construct[state]
        state += 1
        return (mod,state)
    else
        return
    end
end

Base.convert(::Type{Array},c::AbstractConstruct) = [m.m for m in c]
Base.convert(::Type{Array{Mod}},c::AbstractConstruct) = [m for m in c]
##################
#New structs
##################

"""
    OrderedConstruct{T}  <: AbstractConstruct{T}

Structure that contains an ordered construct. In an ordered construct, the order
of the different modules is important.

For example: `OrderedConstruct([a,b,c])` is different compared to
`OrderedConstruct([c,b,a])`.
"""
struct OrderedConstruct{T<:Mod}  <: AbstractConstruct{T}
    c::Array{T}
end

"""
    isequal(c1::OrderedConstruct,c2::OrderedConstruct)

Evaluates if two OrderedConstruct constructs are equal.
"""
Base.isequal(c1::OrderedConstruct,c2::OrderedConstruct) = isequal(c1.c,c2.c)
==(c1::OrderedConstruct,c2::OrderedConstruct) = c1.c == c2.c

function Base.show(io::IO,a::OrderedConstruct)
    construct = convert(Array,a)
    print(io,"{OrderedConstruct}",construct)
end

"""
    *(m1::Mod, m2::Mod)

Returns a `OrderedConstruct` containing both modules `m1` and `m2`. The order of
the input argument is maintained in the returned construct.

To allow efficient on the fly calculation using Kronecker product, the `Base.*`
should be defined between two `Mod`.
"""
*(m1::Mod, m2::Mod) = OrderedConstruct([m1, m2])

"""
    *(c::OrderedConstruct , m::Mod)

Returns a new `OrderedConstruct` where the module `m` is concatenated to the
end of the input construct `c`.
"""
*(c::OrderedConstruct, m::Mod) = push!(c.c,m) |> OrderedConstruct

"""
    *(m::Mod, c::OrderedConstruct)

Returns a new `OrderedConstruct` where the module `m` is concatenated in front
of the input construct `c`.

Probably less efficient than `*(c::OrderedConstruct ,m::Mod)` because of a
`firstpush!`
"""
*(m::Mod, c::OrderedConstruct) = pushfirst!(c.c,m) |> OrderedConstruct

"""
    *(c1::OrderedConstruct,c2::OrderedConstruct)

Concatenated two ordered constructs `c1` and `c2`,
where the construct that is in the second argument position will be appended to the end of the construct in the first argument position

```jldoctest
m1 = Mod("a")
m2 = Mod("b")
m3 = Mod("c")
c1 = OrderedConstruct([m1,m2])
c2 = OrderedConstruct([m1,m3])
c1*c2

# output
OrderedConstruct([m1,m2,m1,m4])
```
"""

*(c1::OrderedConstruct ,c2::OrderedConstruct) = vcat(c1.c,c2.c) |> OrderedConstruct

*(m::Mod, i::Int) = repeat([m],i) |> OrderedConstruct
*(i::Int, m::Mod) = repeat([m],i) |> OrderedConstruct
"""
    UnorderedConstruct{T} <: AbstractConstruct{T}

Structure that contains an unordered construct. In unordered construct, the
order of the different modules has no meaning.

For example: `UnorderedConstruct([a,b,c])` is considered equal to
`UnorderedConstruct([c,b,a])`.
"""
# currently the underlying structure is an array, a Set is can be a more logical alternative
# still for the first implementation and test, an array structure was preferred. An array can always be transformed Set where needed.
# The users shouldn't worry too much about this more to run everything smoothly internally

struct UnorderedConstruct{T<:Mod} <: AbstractConstruct{T}
    c::Array{T}
end

function Base.show(io::IO,a::UnorderedConstruct)
    construct = convert(Array,a)
    s = replace(string(construct), "[" => "{") |> x->  replace(x, "]" => "}")
    print(io,"{UnorderedConstruct}",s)
end

"""
    +(m1::Mod,m2::Mod)

Returns a `UnorderedConstruct` containing both modules `m1` and `m2`. Order of input argument is maintained in the returned construct.
The order of the modules doesn't have useful meaning, but it has some advantages internally to keep the order.
"""

+(m1::Mod,m2::Mod) = UnorderedConstruct([m1,m2])

"""
    +(c::UnorderedConstruct ,m::Mod)

Returns a new `UnorderedConstruct` where the module `m` is concatenated to the
end of the input construct `c`.
"""
+(c::UnorderedConstruct, m::Mod) = push!(c.c,m) |> UnorderedConstruct

"""
    +(c::OrderedConstruct, m::Mod)

Returns a new `UnorderedConstruct` where the module `m` is concatenated in front
of the input construct `c`
"""
+(m::Mod, c::UnorderedConstruct) = pushfirst!(c.c,m) |> UnorderedConstruct

"""
    +(c1::UnorderedConstruct, c2::UnorderedConstruct)

Concatenated two `UnorderedConstruct` `c1` and `c2`,
where the construct that is in the second argument position will be appended to
the end of the construct in the first argument position
"""
+(c1::UnorderedConstruct, c2::UnorderedConstruct) = vcat(c1.c, c2.c) |> UnorderedConstruct

"""
    Base.isequal(c1::UnorderedConstruct, c2::UnorderedConstruct)

Evaluates if two `UnorderedConstruct` are equal.
"""
Base.isequal(c1::UnorderedConstruct, c2::UnorderedConstruct) = isequal(Set(c1.c),Set(c2.c))
==(c1::UnorderedConstruct, c2::UnorderedConstruct) = Set(c1.c) == Set(c2.c)

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
function synth(
    spatial_modes,
    temporal_modes,
    center,
    xvals,
    coefficients,
    cormat,
    vars;
    corthresh = 0.6,
)

    ntrials = size(cormat, 2)

    abscormat = abs.(cormat)
    maxcor = maximum(abscormat, dims = 1)[1, :]
    maxcoridx = mapslices(argmax, abscormat, dims = 1)[1, :]

    proctrials = findall(x -> x >= corthresh, maxcor)

    np = length(proctrials)
    ntime = size(xvals,2)[1]
    spaceshp = size(center)
    nspace = prod(spaceshp)
    nmodes = size(spatial_modes,3)[1]

    # tempout = zeros(
    #     Union{Missing,Float64},
    #     size(spatial_modes)[1:2]...,
    #     size(xvals, 2),
    #     length(proctrials),
    # )

    yvals = zeros(ntime,np)

    for (i, idx) in enumerate(proctrials)

        bestidx = maxcoridx[idx]

        model = Polynomial(coefficients[:, idx])

        yvals[:,i] = model.(xvals[bestidx, :])
    end

    # vars = [0.308162,0.227067,0.195983,0.174960,0.093828]
    varscale = vars[proctrials] / sum(vars[proctrials])

    sm = reshape(spatial_modes,(nspace,nmodes))

    out_ = (yvals .* reshape(varscale,(1,np))) * sm[:,proctrials]'
    # out_ = yvals * sm[:,proctrials]'

    out = reshape(transpose(out_),(spaceshp...,ntime)) .+ reshape(center, (spaceshp..., 1))

    return out
end

function extractwater(img::Array{<:Any,2}; init_thresh = -16)
    binary = img .< init_thresh

    invalidx = findall(x->ismissing(x),binary)
    binary[invalidx] .= 0

    diff1, diff2 = Kernel.sobel()

    xedges = imfilter(binary,diff1)
    yedges = imfilter(binary,diff2)
    edgeimg = @. sqrt(xedges^2 + yedges^2)

    # imgm = (mapwindow(maximum, edgeimg .> 0.5, (11,11))) .& (.!ismissing.(img))
    imgm = (mapwindow(maximum, edgeimg .> 0.5, (11,11)))

    fill_img =img[:,:]
    fill_img[invalidx] = rand(init_thresh*2:0.005:init_thresh,length(invalidx))

    sampleidx = findall(imgm)
    # samples = img[sampleidx]
    samples = fill_img[sampleidx]

    edges, counts = HistogramThresholding.build_histogram(Float32.(samples),500)

    t = find_threshold(counts[1:end], edges, Otsu())

    return img .< t

end

# diff = est .+ ( obs .* 2)
# 0 = both no water
# 1 = est water, obs no water
# 2 = est no water, obs water
# 3 = both water

max_stars_repo_name: SERVIR/fier-cli

max_issues_repo_name: SERVIR/fier-cli

max_forks_repo_name: SERVIR/fier-cli


text:
@enum VaoKind SIMPLE ELEMENTS ELEMENTS_INSTANCED EMPTY

# buffers which are not instanced have divisor = -1
const GEOMETRY_DIVISOR = GLint(-1)
const UNIFORM_DIVISOR = GLint(1)

struct BufferAttachmentInfo{T}
	name     ::Symbol
    location ::GLint
    buffer   ::Buffer{T}
    divisor  ::GLint
end

Base.eltype(b::BufferAttachmentInfo{T}) where T = T

function generate_buffers(program::Program, divisor::GLint; name_buffers...)
	buflen  = 0
    buffers = BufferAttachmentInfo[]
    for (name, val) in pairs(name_buffers)
        loc = attribute_location(program, name)
        if loc != INVALID_ATTRIBUTE
	        buflen = buflen == 0 ? length(val) : buflen
            vallen = length(val)
            if vallen == buflen
                push!(buffers, BufferAttachmentInfo(name, loc, Buffer(val, usage=GL_DYNAMIC_DRAW), divisor))
            elseif !isa(val, Vector)
                push!(buffers, BufferAttachmentInfo(name, loc, Buffer(fill(val, buflen), usage=GL_DYNAMIC_DRAW), divisor))
            end
        end
    end
    return buffers
end

mutable struct VertexArray{Vertex, Kind}
    id         ::GLuint
    bufferinfos::Vector{<:BufferAttachmentInfo}
    indices    ::Union{Buffer, Nothing}
    nverts     ::GLint #total vertices to be drawn in drawcall
    ninst      ::GLint
    face       ::GLenum
    context    ::Context
    function VertexArray(kind::VaoKind, bufferinfos::Vector{<:BufferAttachmentInfo}, indices, ninst, face)
        id = glGenVertexArrays()
        glBindVertexArray(id)

        nverts = 0
        for b in bufferinfos
            nverts_ = length(b.buffer)
            if nverts != 0 && nverts != nverts_ && b.divisor == GEOMETRY_DIVISOR
                error("Amount of vertices is not equal.")
            end
            nverts = nverts_
        end

        if indices != nothing
            bind(indices)
            nverts = length(indices) * cardinality(indices)
        end

        attach2vao.(bufferinfos)

        glBindVertexArray(0)
        if length(bufferinfos) == 1
            if !is_glsl_primitive(eltype(bufferinfos[1].buffer))
                vert_type = eltype(bufferinfos[1].buffer)
            else
                vert_type = Tuple{eltype(bufferinfos[1].buffer)}
            end
        else
            vert_type = Tuple{eltype.(bufferinfos)...}
        end
        obj = new{vert_type, kind}(id, bufferinfos, indices, nverts, ninst, face, current_context())
        finalizer(free!, obj)
        obj
    end
    function VertexArray()
	    return new{Nothing, EMPTY}(GLuint(0),
	                               Buffer[],
	                               nothing,
	                               GLint(0),
	                               GLint(0),
	                               GL_POINTS,
	                               current_context())
	end
end

VertexArray(bufferinfos::Vector{<:BufferAttachmentInfo}, face::GLenum=GL_TRIANGLES) =
    VertexArray(SIMPLE, bufferinfos, nothing, 1, face)

VertexArray(bufferinfos::Vector{<:BufferAttachmentInfo}, indices::Vector{<:Integer}, face::GLenum) =
    VertexArray(ELEMENTS, bufferinfos, indexbuffer(GLint.(indices)), 1, face)

#TODO this is messy
VertexArray(bufferinfos::Vector{<:BufferAttachmentInfo}, indices::Vector{F}; facelength=F) where F =
    VertexArray(ELEMENTS, bufferinfos, indexbuffer(indices), 1, face2glenum(facelength))

VertexArray(bufferinfos::Vector{<:BufferAttachmentInfo}, indices::Vector{<:Integer}, facelength::Union{GLenum, Int}, ninst::Int) =
    VertexArray(ELEMENTS_INSTANCED, bufferinfos, indexbuffer(GLint.(indices)), ninst, face2glenum(facelength))

VertexArray(bufferinfos::Vector{<:BufferAttachmentInfo}, indices::Vector{F}, ninst::Int) where F =
    VertexArray(ELEMENTS_INSTANCED, bufferinfos, indexbuffer(indices), ninst, face2glenum(F))

free!(x::VertexArray) =
    context_command(x.context, () -> glDeleteVertexArrays(1, [x.id]))

is_null(vao::VertexArray{Nothing, EMPTY}) = true
is_null(vao::VertexArray)                 = false

bufferinfo(vao::VertexArray, name::Symbol) = (id =findfirst(x->x.name == name, vao.bufferinfos); id !== nothing ? vao.bufferinfos[id] : id)

# the instanced ones assume that there is at least one buffer with the vertextype (=has fields, bit whishy washy) and the others are the instanced things
# It is assumed that when an attribute is longer than 4 bytes, the rest is stored in consecutive locations
function attach2vao(bufferinfo::BufferAttachmentInfo{T}) where T

    function enable_attrib(loc)
        glEnableVertexAttribArray(loc)
        if bufferinfo.divisor != GEOMETRY_DIVISOR
            glVertexAttribDivisor(loc, bufferinfo.divisor)
        end
    end

    bind(bufferinfo.buffer)
    if !is_glsl_primitive(T)
        # This is for a buffer that holds all the attributes in a OpenGL defined way.
        # This requires us to find the fieldoffset
        # TODO this is not tested
        for i = 1:nfields(T)
            loc = bufferinfo.location + i - 1
            FT = fieldtype(T, i); ET = eltype(FT)
            glVertexAttribPointer(loc,
                                  cardinality(FT), julia2glenum(ET),
                                  GL_FALSE, sizeof(T), Ptr{Nothing}(fieldoffset(T, i)))
            enable_attrib(loc)
        end
    else
        # This is for when the buffer holds a single attribute, no need to
        # calculate fieldoffsets and stuff like that.
        # TODO Assumes everything larger than vec4 is a matrix, is this ok?
        FT = T; ET = eltype(FT)
        cardi = cardinality(FT)
        gltype = julia2glenum(ET)
        if cardi > 4
            s = size(FT)
            loc_size = s[2]
            for li=0:s[1]-1
                loc = bufferinfo.location + li
                offset = sizeof(gltype) * loc_size * li
                glVertexAttribPointer(loc, loc_size, gltype, GL_FALSE, cardi*sizeof(gltype), Ptr{Nothing}(offset))
                enable_attrib(loc)
            end
        else
            glVertexAttribPointer(bufferinfo.location, cardi, gltype, GL_FALSE, 0, C_NULL)
            enable_attrib(bufferinfo.location)
        end
    end
end

# TODO
Base.convert(::Type{VertexArray}, x) = VertexArray(x)
Base.convert(::Type{VertexArray}, x::VertexArray) = x

function face2glenum(face)
    facelength = typeof(face) <: Integer ? face : (face <: Integer ? 1 : length(face))
    facelength == 1  && return GL_POINTS
    facelength == 2  && return GL_LINES
    facelength == 3  && return GL_TRIANGLES
    facelength == 4  && return GL_QUADS
    facelength == 5  && return GL_TRIANGLE_STRIP
    facelength == 10 && return GL_LINES_ADJACENCY
    facelength == 11 && return GL_LINE_STRIP_ADJACENCY
    return GL_TRIANGLES
end

function glenum2face(glenum)
    glenum == GL_POINTS               && return 1
    glenum == GL_LINES                && return 2
    glenum == GL_TRIANGLES            && return 3
    glenum == GL_QUADS                && return 4
    glenum == GL_TRIANGLE_STRIP       && return 5
    glenum == GL_LINES_ADJACENCY      && return 10
    glenum == GL_LINE_STRIP_ADJACENCY && return 11
    return 1
end

is_struct(::Type{T}) where T = !(sizeof(T) != 0 && nfields(T) == 0)
# is_glsl_primitive{T <: StaticVector}(::Type{T}) = true
is_glsl_primitive(::Type{T}) where {T <: Union{Float32, Int32}}= true
function is_glsl_primitive(T)
    glasserteltype(T)
    true
end

_typeof(::Type{T}) where T = Type{T}
_typeof(::T) where T = T

glitype(vao::VertexArray) = julia2glenum(eltype(vao.indices))

Base.length(vao::VertexArray) = vao.nverts

bind(vao::VertexArray) = glBindVertexArray(vao.id)
unbind(vao::VertexArray) = glBindVertexArray(0)

#does this ever work with anything aside from an unsigned int??
draw(vao::VertexArray{V, ELEMENTS} where V) = glDrawElements(vao.face, vao.nverts, GL_UNSIGNED_INT, C_NULL)

draw(vao::VertexArray{V, ELEMENTS_INSTANCED} where V) = glDrawElementsInstanced(vao.face, vao.nverts, GL_UNSIGNED_INT, C_NULL, vao.ninst)

draw(vao::VertexArray{V, SIMPLE} where V) = glDrawArrays(vao.face, 0, vao.nverts)

function Base.show(io::IO, vao::T) where T<:VertexArray
    fields = filter(x->x != :buffers && x!=:indices, [fieldnames(T)...])
    for field in fields
        show(io, getfield(vao, field))
        println(io,"")
    end
end

Base.eltype(::Type{VertexArray{ElTypes, Kind}}) where {ElTypes, Kind} = (ElTypes, Kind)

#This could also be done just with the symbols, not needing new bufferattachment infos
function upload!(vao::VertexArray; name_buffers...)
    for (n, b) in name_buffers
        existing_b = findfirst(x -> x.name == n, vao.bufferinfos)
        upload_data!(vao.bufferinfos[existing_b].buffer, b)
    end
end


text:
"""
    ConstantLinearAdvection{T,Dim}

Linear advection equation with constant coefficients
``
  \\partial_t u + \\sum_{i=1}^{\\mathrm{Dim}} \\partial_{x_i} u = 0
``
in `Dim` space dimensions using `T` as scalar type.
"""
struct ConstantLinearAdvection{T,Dim} <: ScalarBalanceLaw{T,Dim} end

function ConstantLinearAdvection(T=Float64, Dim=1)
  ConstantLinearAdvection{T,Dim}()
end

function Base.show(io::IO, model::ConstantLinearAdvection{T,Dim}) where {T,Dim}
  print(io, "Linear advection equation with constant coefficient {T=", T, ", Dim=", Dim, "}",
            " with flux f(u) = u * (1,...,1)")
end

"""
    flux{T}(u, model::ConstantLinearAdvection{T,1})

Compute the flux of `u` for `model`.
"""
@inline flux(u, model::ConstantLinearAdvection{T,1}) where {T} = u

"""
    speed(u::Real, model::ConstantLinearAdvection)

Compute the speed f'(`u`) for `model`.
"""
@inline speed(u, model::ConstantLinearAdvection) = one(u)

@inline flux(u, model::ConstantLinearAdvection{T,2}, direction) where {T} = u

################################################################################

"""
    (::GodunovFlux){T}(uₗ::T, uᵣ::T, model::ConstantLinearAdvection{T,1})

Compute Godunov's flux between `uₗ` and `uᵣ` for `model`.
"""
function (::GodunovFlux)(uₗ::T, uᵣ::T, model::ConstantLinearAdvection{T,1}) where {T}
    uₗ
end

function (flux::EnergyConservativeFlux)(uₗ::T, uᵣ::T, model::ConstantLinearAdvection{T,1}) where T<:Real
    (uₗ + uᵣ) / 2
end

################################################################################

@inline function add_numerical_fluxes_inner_loop2!(du, fluxes, u, balance_law::ConstantLinearAdvection,
                                                    Nx, basis, boundaries_included::Val{false},
                                                    jacx, parallel)
    Pp1 = length(basis.nodes)
    Rl = basis.interp_left
    Rr = basis.interp_right

    # add numerical fluxes
    @inbounds for ix in Base.OneTo(Nx)
        Rul = zero(eltype(u))
        Rur = zero(eltype(u))
        for nx in Base.OneTo(Pp1)
            Rul += Rl[nx]*u[nx,ix]
            Rur += Rr[nx]*u[nx,ix]
        end

        for nx in Base.OneTo(Pp1)
            du[nx,ix] += ((fluxes[ix] - Rul) * Rl[nx] - (fluxes[ix+1] - Rur) * Rr[nx]
                            ) * jacx / basis.weights[nx]
        end
    end

    nothing
end

################################################################################

"""
    ConstantLinearAdvectionRiemannSolution{T,T1}

The solution of a Riemann problem `prob` for the linear advection equation with
constant coefficients.
"""
struct ConstantLinearAdvectionRiemannSolution{T,T1} <: ScalarRiemannSolution
    prob::RiemannProblem{ConstantLinearAdvection{T,1},T,T1}
    σ::T
end

"""
    minmax_speeds(sol::ConstantLinearAdvectionRiemannSolution)

Return the minimal and maximal speeds `σ⁻, σ⁺` that appear in the solution `sol`.
"""
function minmax_speeds(sol::ConstantLinearAdvectionRiemannSolution)
    sol.σ, sol.σ
end

"""
    (sol::ConstantLinearAdvectionRiemannSolution)(ξ::Real)

Evaluate the solution `sol` at the value `ξ` of the self-similarity variable
`ξ = (x - x₀) / (t - t₀)`.
"""
function (sol::ConstantLinearAdvectionRiemannSolution)(ξ::Real)
    @unpack σ = sol
    @unpack uₗ, uᵣ = sol.prob

    if ξ < σ
        uₗ
    else
        uᵣ
    end
end

"""
    (sol::ConstantLinearAdvectionRiemannSolution)(t::Real, x::Real)

Evaluate the solution `sol` at the time and space coordinates `t` and `x`.
"""
function (sol::ConstantLinearAdvectionRiemannSolution)(t::Real, x::Real)
    @unpack x₀, t₀ = sol.prob

    sol((x-x₀)/(t-t₀))
end

"""
    solve{T,T1}(prob::RiemannProblem{ConstantLinearAdvection{T,1},T,T1})

Compute the solution of the Riemann prolem `prob`.
"""
function solve(prob::RiemannProblem{ConstantLinearAdvection{T,1},T,T1}) where {T,T1}
    ConstantLinearAdvectionRiemannSolution(prob, one(T))
end


text:
abstract type AbstractModelSet end
# using CSV
# using DataFrames

"""
    WindFarmModelSet(wakedeficitmodel, wake_deflection_model, wake_combination_model, local_ti_model)

Container for objects defining models to use in wind farm calculations

# Arguments
- `wake_defiict_model::AbstractWakeDeficitModel`: contains a struct defining the desired wake deficit model
- `wake_deflection_model::AbstractWakeDeflectionModel`: contains a struct defining the desired wake deflection model
- `wake_combination_model::AbstractWakeCombinationModel`: contains a struct defining the desired wake combination model
- `local_ti_model::AbstractTurbulenceIntensityModel`: contains a struct defining the desired turbulence intensity model
"""
struct WindFarmModelSet{DTM,DNM,CM,TIM} <: AbstractModelSet

    wake_deficit_model::DTM
    wake_deflection_model::DNM
    wake_combination_model::CM
    local_ti_model::TIM

end

"""
    point_velocity(loc, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
    wind_resource, model_set::AbstractModelSet;
    wind_farm_state_id=1, downwind_turbine_id=0)

Calculates the wind speed at a given point for a given state

# Arguments
- `loc::Array{TF,3}`: Location of interest
- `turbine_x::Array{TF,nTurbines}`: turbine east-west locations in the state
    reference frame
- `turbine_y::Array{TF,nTurbines}`: turbine north-south locations in the state
    reference frame
- `turbine_z::Array{TF,nTurbines}`: turbine base height in the state reference frame
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `turbine_ct::Array{TF,nTurbines}`: turbine thrust coefficients for the given state
- `turbine_ai::Array{TF,nTurbines}`: turbine axial induction for the given state
- `rotor_diameter::Array{TF,nTurbines}`: turbine rotor diameters
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_local_ti::Array{TF,nTurbines}`: turbine local turbulence intensity for
    the given state
- `sorted_turbine_index::Array{TF,nTurbines}`: array containing indices of wind turbines
    from most upwind to most downwind turbine in the given state
- `wtvelocities::Array{TF,nTurbines}`: effective inflow wind speed for given state
- `wind_resource::DiscretizedWindResource`: contains wind resource discreption (directions,
    speeds, frequencies, etc)
- `wind_farm_state_id::Int`: index to correct state to use from wind resource provided.
    Defaults to 1
- `downwind_turbine_id::Int`: index of wind turbine of interest (if any). If not a point for
    calculating effective wind speed of a turbine, then provide 0 (default)
"""
function point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
                    wind_resource, model_set::AbstractModelSet;
                    wind_farm_state_id=1, downwind_turbine_id=0)

    wakedeficitmodel = model_set.wake_deficit_model
    wakedeflectionmodel = model_set.wake_deflection_model
    wakecombinationmodel = model_set.wake_combination_model

    # extract flow information
    wind_speed = wind_resource.wind_speeds[wind_farm_state_id]
    reference_height = wind_resource.measurement_heights[wind_farm_state_id]

    # set ground height
    ground_height = wind_resource.wind_shear_model.ground_height    # TODO: allow topology to be given

    # find order for wind shear and deficit calculations
    shear_order = wind_resource.wind_shear_model.shear_order

    # adjust wind speed for wind shear
    if shear_order == "nothing"
        wind_speed_internal = wind_speed
    elseif shear_order == "first"
        wind_speed_internal = adjust_for_wind_shear(locz, wind_speed, reference_height, ground_height, wind_resource.wind_shear_model)
    else
        wind_speed_internal = wind_speed
    end

    # get number of turbines
    nturbines = length(turbine_x)

    # initialize deficit summation term to zero
    deficit_sum = 0.0

    # loop through all turbines
    for u=1:nturbines

        # get index of upstream turbine
        upwind_turb_id = Int(sorted_turbine_index[u])

        # don't allow turbine to impact itself
        if upwind_turb_id == downwind_turbine_id; continue; end

        # downstream distance between upstream turbine and point
        x = locx - turbine_x[upwind_turb_id]

        # check turbine relative locations
        if x > 1E-6
            # skip this loop if it would include a turbine's impact on itself)
            if upwind_turb_id==downwind_turbine_id; continue; end

            # calculate wake deflection of the current wake at the point of interest
            horizontal_deflection = wake_deflection_model(locx, locy, locz, turbine_x, turbine_yaw, turbine_ct,
                            upwind_turb_id, rotor_diameter, turbine_local_ti, wakedeflectionmodel)

            vertical_deflection = 0.0

            # velocity difference in the wake
            deltav = wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, horizontal_deflection, vertical_deflection,
                            upwind_turb_id, downwind_turbine_id, hub_height, rotor_diameter, turbine_ai,
                            turbine_local_ti, turbine_ct, turbine_yaw, wakedeficitmodel)

            # combine deficits according to selected wake combination method
            deficit_sum = wake_combination_model(deltav, wind_speed_internal, wtvelocities[upwind_turb_id], deficit_sum, wakecombinationmodel)
            # println(deficit_sum, " ", downwind_turbine_id, " ", upwind_turb_id)

        end
    end

    # find velocity at point without shear
    point_velocity = wind_speed_internal - deficit_sum

    if shear_order == "nothing"
        point_velocity_out = point_velocity
    elseif shear_order == "first"
        point_velocity_out = point_velocity
    else
        point_velocity_out = adjust_for_wind_shear(locz, point_velocity, reference_height, ground_height, wind_resource.wind_shear_model)
    end

    return point_velocity_out

end

"""
    point_velocity(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
    sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
    model_set::AbstractModelSet; wind_farm_state_id=1)

Calculates the wind speed at a given point for a given state

# Arguments
- `turbine_x::Array{TF,nTurbines}`: turbine east-west locations in the state
    reference frame
- `turbine_y::Array{TF,nTurbines}`: turbine north-south locations in the state
    reference frame
- `turbine_z::Array{TF,nTurbines}`: turbine base height in the state reference frame
- `rotor_diameter::Array{TF,nTurbines}`: turbine rotor diameters
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `sorted_turbine_index::Array{TF,nTurbines}`: turbine sorted order upstream to downstream
    for given state
- `ct_model::AbstractThrustCoefficientModel`: defines how the thrust coefficient changes
    with state etc
- rotor_sample_points_y::Array{TF,N}`: horizontal wind location of points to sample across
    the rotor swept area when calculating the effective wind speed for the wind turbine.
    Points are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- rotor_sample_points_z::Array{TF,N}`: vertical wind location of points to sample across the
    rotor swept area when calculating the effective wind speed for the wind turbine. Points
    are centered at the hub (0,0) and scaled by the radius (1=tip of blades)
- `wind_resource::DiscretizedWindResource`: wind resource discreption (directions, speeds,
    frequencies, etc)
- `model_set::AbstractModelSet`: defines wake-realated models to be used in analysis
- `wind_farm_state_id::Int`: index to correct state to use from wind resource provided.
    Defaults to 1
"""
function turbine_velocities_one_direction(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
                    sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
                    model_set::AbstractModelSet; wind_farm_state_id=1, velocity_only=true)

    # get number of turbines and rotor sample point
    n_turbines = length(turbine_x)
    n_rotor_sample_points = length(rotor_sample_points_y)

    # initialize correct array types
    arr_type = promote_type(typeof(turbine_x[1]),typeof(turbine_y[1]),typeof(turbine_z[1]),typeof(rotor_diameter[1]),
                            typeof(hub_height[1]),typeof(turbine_yaw[1]))

    # initialize arrays
    turbine_velocities = zeros(arr_type, n_turbines)
    turbine_ct = zeros(arr_type, n_turbines)
    turbine_ai = zeros(arr_type, n_turbines)
    turbine_local_ti = zeros(arr_type, n_turbines)

    # loop over all turbines
    for d=1:n_turbines

        # get index of downstream turbine
        downwind_turbine_id = Int(sorted_turbine_index[d])

        # initialize downstream wind turbine velocity to zero
        # println("start array: ", turbine_velocities[downwind_turbine_id])
        # wind_turbine_velocity = typeof(turbine_velocities[downwind_turbine_id])(0.0)
        wind_turbine_velocity = 0.0
        # turbine_velocities[downwind_turbine_id] = 0.0

        # loop over all rotor sample points to approximate the effective inflow velocity
        for p=1:n_rotor_sample_points

            # scale rotor sample point coordinate by rotor diameter (in rotor hub ref. frame)
            local_rotor_sample_point_y = rotor_sample_points_y[p]*0.5*rotor_diameter[downwind_turbine_id]
            local_rotor_sample_point_z = rotor_sample_points_z[p]*0.5*rotor_diameter[downwind_turbine_id]

            # put rotor sample points in wind direction coordinate system, and account for yaw
            locx = turbine_x[downwind_turbine_id] .+ local_rotor_sample_point_y*sin(turbine_yaw[downwind_turbine_id])
            locy = turbine_y[downwind_turbine_id] .+ local_rotor_sample_point_y*cos(turbine_yaw[downwind_turbine_id])
            locz = turbine_z[downwind_turbine_id] .+ hub_height[downwind_turbine_id] + local_rotor_sample_point_z

            # calculate the velocity at given point
            point_velocity_with_shear = point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, turbine_velocities,
                                    wind_resource, model_set,
                                    wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=downwind_turbine_id)

            # add sample point velocity to turbine velocity to be averaged later
            wind_turbine_velocity += point_velocity_with_shear

        end

        # final velocity calculation for downstream turbine (average equally across all points)
        wind_turbine_velocity /= n_rotor_sample_points

        turbine_velocities[downwind_turbine_id] = deepcopy(wind_turbine_velocity)

        # update thrust coefficient for downstream turbine
        turbine_ct[downwind_turbine_id] = calculate_ct(turbine_velocities[downwind_turbine_id], ct_model[downwind_turbine_id])

        # update axial induction for downstream turbine
        turbine_ai[downwind_turbine_id] = _ct_to_axial_ind_func(turbine_ct[downwind_turbine_id])

        # get local turbulence intensity for this wind state
        ambient_ti = wind_resource.ambient_tis[wind_farm_state_id]

        # update local turbulence intensity for downstream turbine
        turbine_local_ti[downwind_turbine_id] = calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
                            turbine_velocities, turbine_ct, model_set.local_ti_model; turbine_id=downwind_turbine_id, tol=1E-6)

        # println("local ti turb 9: ", turbine_local_ti[downwind_turbine_id])
    end

    # df = DataFrame(ID=1:n_turbines, V=turbine_velocities, TI=turbine_local_ti, CT=turbine_ct)
    # CSV.write("internaldata.txt", df)

    if velocity_only
        return turbine_velocities
    else
        return turbine_velocities, turbine_ct, turbine_ai, turbine_local_ti
    end
end

function turbine_velocities_one_direction(x, turbine_z, rotor_diameter, hub_height, turbine_yaw,
    sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
    model_set::AbstractModelSet; wind_farm_state_id=1, velocity_only=true)

    n_turbines = Int(length(x)/2)
    # println(typeof(x), n_turbines)
    turbine_x = x[1:n_turbines]
    turbine_y = x[n_turbines+1:end]
    # println(turbine_x)
    # println("turbine_x type ", typeof(turbine_x))
    # println("type of x ", typeof(x))

    # get number of turbines and rotor sample point
    # n_turbines = length(turbine_x)
    n_rotor_sample_points = length(rotor_sample_points_y)

    arr_type = promote_type(typeof(turbine_x[1]),typeof(turbine_y[1]),typeof(turbine_z[1]),typeof(rotor_diameter[1]),
                typeof(hub_height[1]),typeof(turbine_yaw[1]))
    turbine_velocities = zeros(arr_type, n_turbines)
    turbine_ct = zeros(arr_type, n_turbines)
    turbine_ai = zeros(arr_type, n_turbines)
    turbine_local_ti = zeros(arr_type, n_turbines)

    for d=1:n_turbines

        # get index of downstream turbine
        downwind_turbine_id = Int(sorted_turbine_index[d])

        # initialize downstream wind turbine velocity to zero
        # println("start array: ", turbine_velocities[downwind_turbine_id])
        # wind_turbine_velocity = typeof(turbine_velocities[downwind_turbine_id])(0.0)
        wind_turbine_velocity = 0.0
        # turbine_velocities[downwind_turbine_id] = 0.0

        for p=1:n_rotor_sample_points

            # scale rotor sample point coordinate by rotor diameter (in rotor hub ref. frame)
            local_rotor_sample_point_y = rotor_sample_points_y[p]*0.5*rotor_diameter[downwind_turbine_id]
            local_rotor_sample_point_z = rotor_sample_points_z[p]*0.5*rotor_diameter[downwind_turbine_id]

            locx = turbine_x[downwind_turbine_id] .+ local_rotor_sample_point_y*sin(turbine_yaw[downwind_turbine_id])
            locy = turbine_y[downwind_turbine_id] .+ local_rotor_sample_point_y*cos(turbine_yaw[downwind_turbine_id])
            locz = turbine_z[downwind_turbine_id] .+ hub_height[downwind_turbine_id] + local_rotor_sample_point_z

            # calculate the velocity at given point
            point_velocity_with_shear = point_velocity(locx, locy, locz, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                                rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, turbine_velocities,
                                wind_resource, model_set,
                                wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=downwind_turbine_id)

            # add sample point velocity to turbine velocity to be averaged later
            wind_turbine_velocity += point_velocity_with_shear

        end

        # final velocity calculation for downstream turbine (average equally across all points)
        wind_turbine_velocity /= n_rotor_sample_points

        turbine_velocities[downwind_turbine_id] = deepcopy(wind_turbine_velocity)

        # update thrust coefficient for downstream turbine
        turbine_ct[downwind_turbine_id] = calculate_ct(turbine_velocities[downwind_turbine_id], ct_model[downwind_turbine_id])

        # update axial induction for downstream turbine
        turbine_ai[downwind_turbine_id] = _ct_to_axial_ind_func(turbine_ct[downwind_turbine_id])

        # update local turbulence intensity for downstream turbine
        ambient_ti = wind_resource.ambient_tis[wind_farm_state_id]
        turbine_local_ti[downwind_turbine_id] = calculate_local_ti(turbine_x, turbine_y, ambient_ti, rotor_diameter, hub_height, turbine_yaw, turbine_local_ti, sorted_turbine_index,
                    turbine_velocities, turbine_ct, model_set.local_ti_model; turbine_id=downwind_turbine_id, tol=1E-6)

    end

    if velocity_only
        return turbine_velocities
    else
        return turbine_velocities, turbine_ct, turbine_ai, turbine_local_ti
    end

end

# turbine_velocities_one_direction!(model_set::AbstractModelSet, problem_description::AbstractWindFarmProblem; wind_farm_state_id=1) = turbine_velocities_one_direction!([0.0], [0.0],
# model_set::AbstractModelSet, problem_description::AbstractWindFarmProblem; wind_farm_state_id=1)

"""
calculate_flow_field(xrange, yrange, zrange, model_set::AbstractModelSet, turbine_x,
    turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai, rotor_diameter, hub_height,
    turbine_local_ti, sorted_turbine_index, wtvelocities, wind_resource; wind_farm_state_id=1)

Generates a flow field for a given state and cross section

# Arguments
- `xrange::Range`: range defining east-west locations to sample in global reference frame
- `yrange::Range`: range defining north-west locations to sample in global reference frame
- `zrange::Range`: range defining vertical locations to sample in global reference frame
- `model_set::AbstractModelSet`: defines wake-realated models to be used in analysis
- `turbine_x::Array{TF,nTurbines}`: turbine east-west locations in the global
    reference frame
- `turbine_y::Array{TF,nTurbines}`: turbine north-south locations in the global
    reference frame
- `turbine_z::Array{TF,nTurbines}`: turbine base height in the global reference frame
- `turbine_yaw::Array{TF,nTurbines}`: turbine yaw for the given wind direction in
    radians
- `turbine_ct::Array{TF,nTurbines}`: thrust coefficient of each turbine for the given state
- `turbine_ai::Array{TF,nTurbines}`: turbine axial induction for the given state
- `rotor_diameter::Array{TF,nTurbines}`: turbine rotor diameters
- `hub_height::Array{TF,nTurbines}`: turbine hub heights
- `turbine_local_ti::Array{TF,nTurbines}`: turbine local turbulence intensity for
    the given state
- `sorted_turbine_index::Array{TF,nTurbines}`: turbine north-south locations in the
    global reference frame
- `wtvelocities::Array{TF,nTurbines}`: effective inflow wind speed for given state
- `wind_resource::DiscretizedWindResource`: wind resource discreption (directions, speeds,
    frequencies, etc)
- `wind_farm_state_id::Int`: index to correct state to use from wind resource provided.
    Defaults to 1
"""
function calculate_flow_field(xrange, yrange, zrange,
    model_set::AbstractModelSet, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
    wind_resource; wind_farm_state_id=1)

    xlen = length(xrange)
    ylen = length(yrange)
    zlen = length(zrange)
    npoints = xlen*ylen*zlen
    point_velocities = zeros(npoints)
    point_velocities = reshape(point_velocities, (zlen, ylen, xlen))

    # rotate to direction frame for velocity calculations
    rot_tx, rot_ty = rotate_to_wind_direction(turbine_x, turbine_y, wind_resource.wind_directions[wind_farm_state_id])

    # sort the turbines
    sorted_turbine_index = sortperm(rot_tx)

    for zi in 1:zlen
        for yi in 1:ylen
            for xi in 1:xlen
                locx = xrange[xi]
                locy = yrange[yi]
                locz = zrange[zi]
                locx, locy = rotate_to_wind_direction(locx, locy, wind_resource.wind_directions[wind_farm_state_id])

                point_velocities[zi, yi, xi] = point_velocity(locx, locy, locz, rot_tx, rot_ty, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
                    rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, wtvelocities,
                    wind_resource, model_set,
                    wind_farm_state_id=wind_farm_state_id, downwind_turbine_id=0)

            end
        end
    end

    # if zlen == 1
    #     return point_velocities[1,1:ylen,1:xlen]
    # elseif ylen == 1
    #     return point_velocities[1:zlen,1,1:xlen]
    # elseif xlen == 1
    #     return point_velocities[1:zlen,1:ylen,1]
    # else
    return point_velocities[1:zlen,1:ylen,1:xlen]
    # end

end

function calculate_flow_field(xrange, yrange, zrange,
    model_set::AbstractModelSet, turbine_x, turbine_y, turbine_z, turbine_yaw,
    rotor_diameter, hub_height, ct_models, rotor_sample_points_y, rotor_sample_points_z,
    wind_resource; wind_farm_state_id=1)

    # rotate to direction frame for velocity calculations
    rot_tx, rot_ty = rotate_to_wind_direction(turbine_x, turbine_y, wind_resource.wind_directions[wind_farm_state_id])

    # sort the turbines
    sorted_turbine_index = sortperm(rot_tx)

    turbine_velocities, turbine_ct, turbine_ai, turbine_local_ti = turbine_velocities_one_direction(rot_tx, rot_ty, turbine_z, rotor_diameter, hub_height, turbine_yaw,
    sorted_turbine_index, ct_models, rotor_sample_points_y, rotor_sample_points_z, wind_resource,
    model_set, wind_farm_state_id=wind_farm_state_id, velocity_only=false)

    return calculate_flow_field(xrange, yrange, zrange,
        model_set, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
        rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, turbine_velocities,
        wind_resource, wind_farm_state_id=wind_farm_state_id)

end


text:
module Simulation73

using DrWatson
using DiffEqBase, OrdinaryDiffEq, DiffEqCallbacks#, DiffEqParamEstim
import OrdinaryDiffEq: OrdinaryDiffEqAlgorithm, ODEProblem, Euler
import DiffEqBase: DESolution, solve, ODEProblem, DECallback
using OrdinaryDiffEq
using Lazy
using AxisIndices
using Parameters

# "variables.jl"
export AbstractVariable, UnboundedVariable, BoundedVariable,
	default_value, bounds, pops, MaybeVariable,
	AbstractParameter, AbstractAction, AbstractSpaceAction

# space.jl
export AbstractSpace, AbstractLattice, AbstractPeriodicLattice, AbstractCompactLattice,
	AbstractEmbeddedLattice

export CompactLattice, PeriodicLattice

export RandomlyEmbeddedLattice, unembed_values

export coordinates, origin_idx, differences, coordinate_axes, timepoints, space,
    extent, abs_difference, abs_difference_periodic, discrete_segment, fft_center_dx

# "action.jl"
export AbstractAction, AbstractInteraction, AbstractSpaceAction, AbstractSpaceInteraction
export NullifyParameter, NullAction

# "population.jl"
export AbstractPopulationActionsParameters, AbstractPopulationInteractionsParameters,
    AbstractPopulationActions, AbstractPopulationInteractions,
    population, population_coordinates, population_timepoint, population_repeat,
    PopulationActionsParameters, PopulationInteractionsParameters,
    PopulationActions2, PopActParam, PopInteractParam, PopAct, PopInteract

# "subsampling.jl" (note: should probably be meshed with meshes)
export scalar_to_idx_window, subsampling_Δidx, subsampling_idxs,
	subsampling_time_idxs, subsampling_space_idxs, AbstractSubsampler,
    IndexSubsampler, ValueSubsampler, ValueWindower, RadialSlice, StrideToEnd

# "stimulus.jl"
export AbstractStimulusParameter, AbstractStimulusAction,
    NoStimulusParameter, NoStimulusAction

# "analysing.jl"
export AbstractPlotSpecification, AbstractSpaceTimePlotSpecification, Analyses,
	output_name, plot_and_save, analyse, subsample, subsampling_idxs

# "targets.jl"
export AbstractTarget, target_loss

export execute

# "simulating.jl"
export AbstractModel, AbstractODEModel,
     AbstractModelwithDelay, Solver, Simulation,
    AbstractExecution, Execution, FailedExecution,
	initial_value, history, time_span, saved_dt, saved_dx,
	generate_problem, solve, run_simulation,
	make_mutators, make_system_mutator,
    BareSolution, reduced_space, frame_xs

# # "exploring.jl"
# export Search, SearchExecution, make_problem_generator, search, run_search

abstract type AbstractParameter{T} end
include("helpers.jl")
include("deconstructing.jl")
include("variables.jl")
include("space.jl")
include("parameter.jl")
include("action.jl")
include("population.jl")
include("subsampling.jl") # depends on space.jl
include("stimulus.jl") # depends on variables.jl, population.jl, and space.jl
include("solutions.jl")
#include("solvers.jl")
include("simulating.jl")
include("targets.jl")
# include("exploring.jl")
include("analysing.jl")

end