text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct IEEETurbineGov1 <: TurbineGov
        K::Float64
        T1::Float64
        T2::Float64
        T3::Float64
        U0::Float64
        U_c::Float64
        valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        T4::Float64
        K1::Float64
        K2::Float64
        T5::Float64
        K3::Float64
        K4::Float64
        T6::Float64
        K5::Float64
        K6::Float64
        T7::Float64
        K7::Float64
        K8::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

IEEE Type 1 Speed-Governing Model

# Arguments
- `K::Float64`: Governor Gain, validation range: `(5, 30)`, action if invalid: `warn`
- `T1::Float64`: Input Filter Lag, validation range: `(0, 5)`, action if invalid: `warn`
- `T2::Float64`: Input Filter Lead, validation range: `(0, 10)`, action if invalid: `warn`
- `T3::Float64`: Valve position Time Constant, validation range: `(eps(), 1)`, action if invalid: `error`
- `U0::Float64`: Maximum Valve Opening Rate, validation range: `(0.01, 0.03)`, action if invalid: `warn`
- `U_c::Float64`: Maximum Valve closing rate, validation range: `(-0.3, 0)`, action if invalid: `warn`
- `valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Valve position limits in MW
- `T4::Float64`: Time Constant inlet steam, validation range: `(0, 1)`, action if invalid: `warn`
- `K1::Float64`: Fraction of high presure shaft power, validation range: `(-2, 1)`, action if invalid: `warn`
- `K2::Float64`: Fraction of low presure shaft power, validation range: `(0, nothing)`, action if invalid: `warn`
- `T5::Float64`: Time constant for second boiler pass, validation range: `(0, 10)`, action if invalid: `warn`
- `K3::Float64`: Fraction of high presure shaft power second boiler pass, validation range: `(0, 0.5)`, action if invalid: `warn`
- `K4::Float64`: Fraction of low presure shaft power second boiler pass, validation range: `(0, 0.5)`, action if invalid: `warn`
- `T6::Float64`: Time constant for third boiler pass, validation range: `(0, 10)`, action if invalid: `warn`
- `K5::Float64`: Fraction of high presure shaft power third boiler pass, validation range: `(0, 0.35)`, action if invalid: `warn`
- `K6::Float64`: Fraction of low presure shaft power third boiler pass, validation range: `(0, 0.55)`, action if invalid: `warn`
- `T7::Float64`: Time constant for fourth boiler pass, validation range: `(0, 10)`, action if invalid: `warn`
- `K7::Float64`: Fraction of high presure shaft power fourth boiler pass, validation range: `(0, 0.3)`, action if invalid: `warn`
- `K8::Float64`: Fraction of low presure shaft power fourth boiler pass, validation range: `(0, 0.3)`, action if invalid: `warn`
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the IEEETurbineGov model are:
	x_g1: First Governor integrator,
	x_g2: Governor output,
	x_g3: First Turbine integrator,
	x_g4: Second Turbine Integrator,
	x_g5: Third Turbine Integrator,
	x_g6: Fourth Turbine Integrator,
- `n_states::Int`: IEEEG1 has 6 states
- `states_types::Vector{StateTypes}`: IEEEG1 has 6 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct IEEETurbineGov1 <: TurbineGov
    "Governor Gain"
    K::Float64
    "Input Filter Lag"
    T1::Float64
    "Input Filter Lead"
    T2::Float64
    "Valve position Time Constant"
    T3::Float64
    "Maximum Valve Opening Rate"
    U0::Float64
    "Maximum Valve closing rate"
    U_c::Float64
    "Valve position limits in MW"
    valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Time Constant inlet steam"
    T4::Float64
    "Fraction of high presure shaft power"
    K1::Float64
    "Fraction of low presure shaft power"
    K2::Float64
    "Time constant for second boiler pass"
    T5::Float64
    "Fraction of high presure shaft power second boiler pass"
    K3::Float64
    "Fraction of low presure shaft power second boiler pass"
    K4::Float64
    "Time constant for third boiler pass"
    T6::Float64
    "Fraction of high presure shaft power third boiler pass"
    K5::Float64
    "Fraction of low presure shaft power third boiler pass"
    K6::Float64
    "Time constant for fourth boiler pass"
    T7::Float64
    "Fraction of high presure shaft power fourth boiler pass"
    K7::Float64
    "Fraction of low presure shaft power fourth boiler pass"
    K8::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the IEEETurbineGov model are:
	x_g1: First Governor integrator,
	x_g2: Governor output,
	x_g3: First Turbine integrator,
	x_g4: Second Turbine Integrator,
	x_g5: Third Turbine Integrator,
	x_g6: Fourth Turbine Integrator, "
    states::Vector{Symbol}
    "IEEEG1 has 6 states"
    n_states::Int
    "IEEEG1 has 6 differential states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function IEEETurbineGov1(K, T1, T2, T3, U0, U_c, valve_position_limits, T4, K1, K2, T5, K3, K4, T6, K5, K6, T7, K7, K8, P_ref=1.0, ext=Dict{String, Any}(), )
    IEEETurbineGov1(K, T1, T2, T3, U0, U_c, valve_position_limits, T4, K1, K2, T5, K3, K4, T6, K5, K6, T7, K7, K8, P_ref, ext, [:x_g1, :x_g2, :x_g3, :x_g4, :x_g5, :x_g6], 6, [StateTypes.Differential, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid], InfrastructureSystemsInternal(), )
end

function IEEETurbineGov1(; K, T1, T2, T3, U0, U_c, valve_position_limits, T4, K1, K2, T5, K3, K4, T6, K5, K6, T7, K7, K8, P_ref=1.0, ext=Dict{String, Any}(), states=[:x_g1, :x_g2, :x_g3, :x_g4, :x_g5, :x_g6], n_states=6, states_types=[StateTypes.Differential, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid], internal=InfrastructureSystemsInternal(), )
    IEEETurbineGov1(K, T1, T2, T3, U0, U_c, valve_position_limits, T4, K1, K2, T5, K3, K4, T6, K5, K6, T7, K7, K8, P_ref, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function IEEETurbineGov1(::Nothing)
    IEEETurbineGov1(;
        K=0,
        T1=0,
        T2=0,
        T3=0,
        U0=0,
        U_c=0,
        valve_position_limits=(min=0.0, max=0.0),
        T4=0,
        K1=0,
        K2=0,
        T5=0,
        K3=0,
        K4=0,
        T6=0,
        K5=0,
        K6=0,
        T7=0,
        K7=0,
        K8=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`IEEETurbineGov1`](@ref) `K`."""
get_K(value::IEEETurbineGov1) = value.K
"""Get [`IEEETurbineGov1`](@ref) `T1`."""
get_T1(value::IEEETurbineGov1) = value.T1
"""Get [`IEEETurbineGov1`](@ref) `T2`."""
get_T2(value::IEEETurbineGov1) = value.T2
"""Get [`IEEETurbineGov1`](@ref) `T3`."""
get_T3(value::IEEETurbineGov1) = value.T3
"""Get [`IEEETurbineGov1`](@ref) `U0`."""
get_U0(value::IEEETurbineGov1) = value.U0
"""Get [`IEEETurbineGov1`](@ref) `U_c`."""
get_U_c(value::IEEETurbineGov1) = value.U_c
"""Get [`IEEETurbineGov1`](@ref) `valve_position_limits`."""
get_valve_position_limits(value::IEEETurbineGov1) = value.valve_position_limits
"""Get [`IEEETurbineGov1`](@ref) `T4`."""
get_T4(value::IEEETurbineGov1) = value.T4
"""Get [`IEEETurbineGov1`](@ref) `K1`."""
get_K1(value::IEEETurbineGov1) = value.K1
"""Get [`IEEETurbineGov1`](@ref) `K2`."""
get_K2(value::IEEETurbineGov1) = value.K2
"""Get [`IEEETurbineGov1`](@ref) `T5`."""
get_T5(value::IEEETurbineGov1) = value.T5
"""Get [`IEEETurbineGov1`](@ref) `K3`."""
get_K3(value::IEEETurbineGov1) = value.K3
"""Get [`IEEETurbineGov1`](@ref) `K4`."""
get_K4(value::IEEETurbineGov1) = value.K4
"""Get [`IEEETurbineGov1`](@ref) `T6`."""
get_T6(value::IEEETurbineGov1) = value.T6
"""Get [`IEEETurbineGov1`](@ref) `K5`."""
get_K5(value::IEEETurbineGov1) = value.K5
"""Get [`IEEETurbineGov1`](@ref) `K6`."""
get_K6(value::IEEETurbineGov1) = value.K6
"""Get [`IEEETurbineGov1`](@ref) `T7`."""
get_T7(value::IEEETurbineGov1) = value.T7
"""Get [`IEEETurbineGov1`](@ref) `K7`."""
get_K7(value::IEEETurbineGov1) = value.K7
"""Get [`IEEETurbineGov1`](@ref) `K8`."""
get_K8(value::IEEETurbineGov1) = value.K8
"""Get [`IEEETurbineGov1`](@ref) `P_ref`."""
get_P_ref(value::IEEETurbineGov1) = value.P_ref
"""Get [`IEEETurbineGov1`](@ref) `ext`."""
get_ext(value::IEEETurbineGov1) = value.ext
"""Get [`IEEETurbineGov1`](@ref) `states`."""
get_states(value::IEEETurbineGov1) = value.states
"""Get [`IEEETurbineGov1`](@ref) `n_states`."""
get_n_states(value::IEEETurbineGov1) = value.n_states
"""Get [`IEEETurbineGov1`](@ref) `states_types`."""
get_states_types(value::IEEETurbineGov1) = value.states_types
"""Get [`IEEETurbineGov1`](@ref) `internal`."""
get_internal(value::IEEETurbineGov1) = value.internal

"""Set [`IEEETurbineGov1`](@ref) `K`."""
set_K!(value::IEEETurbineGov1, val) = value.K = val
"""Set [`IEEETurbineGov1`](@ref) `T1`."""
set_T1!(value::IEEETurbineGov1, val) = value.T1 = val
"""Set [`IEEETurbineGov1`](@ref) `T2`."""
set_T2!(value::IEEETurbineGov1, val) = value.T2 = val
"""Set [`IEEETurbineGov1`](@ref) `T3`."""
set_T3!(value::IEEETurbineGov1, val) = value.T3 = val
"""Set [`IEEETurbineGov1`](@ref) `U0`."""
set_U0!(value::IEEETurbineGov1, val) = value.U0 = val
"""Set [`IEEETurbineGov1`](@ref) `U_c`."""
set_U_c!(value::IEEETurbineGov1, val) = value.U_c = val
"""Set [`IEEETurbineGov1`](@ref) `valve_position_limits`."""
set_valve_position_limits!(value::IEEETurbineGov1, val) = value.valve_position_limits = val
"""Set [`IEEETurbineGov1`](@ref) `T4`."""
set_T4!(value::IEEETurbineGov1, val) = value.T4 = val
"""Set [`IEEETurbineGov1`](@ref) `K1`."""
set_K1!(value::IEEETurbineGov1, val) = value.K1 = val
"""Set [`IEEETurbineGov1`](@ref) `K2`."""
set_K2!(value::IEEETurbineGov1, val) = value.K2 = val
"""Set [`IEEETurbineGov1`](@ref) `T5`."""
set_T5!(value::IEEETurbineGov1, val) = value.T5 = val
"""Set [`IEEETurbineGov1`](@ref) `K3`."""
set_K3!(value::IEEETurbineGov1, val) = value.K3 = val
"""Set [`IEEETurbineGov1`](@ref) `K4`."""
set_K4!(value::IEEETurbineGov1, val) = value.K4 = val
"""Set [`IEEETurbineGov1`](@ref) `T6`."""
set_T6!(value::IEEETurbineGov1, val) = value.T6 = val
"""Set [`IEEETurbineGov1`](@ref) `K5`."""
set_K5!(value::IEEETurbineGov1, val) = value.K5 = val
"""Set [`IEEETurbineGov1`](@ref) `K6`."""
set_K6!(value::IEEETurbineGov1, val) = value.K6 = val
"""Set [`IEEETurbineGov1`](@ref) `T7`."""
set_T7!(value::IEEETurbineGov1, val) = value.T7 = val
"""Set [`IEEETurbineGov1`](@ref) `K7`."""
set_K7!(value::IEEETurbineGov1, val) = value.K7 = val
"""Set [`IEEETurbineGov1`](@ref) `K8`."""
set_K8!(value::IEEETurbineGov1, val) = value.K8 = val
"""Set [`IEEETurbineGov1`](@ref) `P_ref`."""
set_P_ref!(value::IEEETurbineGov1, val) = value.P_ref = val
"""Set [`IEEETurbineGov1`](@ref) `ext`."""
set_ext!(value::IEEETurbineGov1, val) = value.ext = val
"""Set [`IEEETurbineGov1`](@ref) `states_types`."""
set_states_types!(value::IEEETurbineGov1, val) = value.states_types = val


text:
#------------------------------------------------------------------------------
# Author: Edward Tate <edward.tate@erlang-solutions.com>
#------------------------------------------------------------------------------

export decode, decode_message
export decode_integer

#------------------------------------------------------------------------------
# External API
#------------------------------------------------------------------------------

function decode_message (bytes::Array{Uint8})
    arity = decode_integer(bytes[1:4])
    # Skip byte 5 as this is the magic number 131
    bytes = bytes[6:end]
    pkt, nbytes = decode(bytes)
    pkt, nbytes + 5
end

#------------------------------------------------------------------------------
# Decoding Functions
#------------------------------------------------------------------------------

function decode (bytes::Array{Uint8})
    decoded, nbytes = decode(bytes[1], bytes[2:end])
    decoded, nbytes + 1
end

function decode (tag::Uint8, rest::Array{Uint8})
    tag == 70 ? decode(ErlFloat, rest) :
    tag == 97 ? decode(ErlSmallInteger, rest) :
    tag == 98 ? decode(ErlInteger, rest) :
    tag == 100 ? decode(ErlAtom, rest) :
    tag == 104 ? decode(ErlSmallTuple, rest) :
    tag == 105 ? decode(ErlLargeTuple, rest) :
    tag == 106 ? decode(ErlNil, rest) :
    tag == 107 ? decode(ErlString, rest) :
    tag == 108 ? decode(ErlList, rest) :
    tag == 109 ? decode(ErlBinary, rest) :
    tag == 110 ? decode(ErlSmallBigInt, rest) :
    tag == 111 ? decode(ErlLargeBigInt, rest) :
    error("Decoding tag not recognized")
end

function decode (::Type{ErlFloat}, bytes::Array{Uint8})
    hex2num(bytes2hex(bytes[1:8])), 8
end

function decode (::Type{ErlSmallInteger}, bytes::Array{Uint8})
    bytes[1], 1
end

function decode (::Type{ErlInteger}, bytes::Array{Uint8})
    decode_integer(bytes[1:4]), 4
end

function decode (::Type{ErlAtom}, bytes::Array{Uint8})
    arity = decode_integer(bytes[1:2])
    atom = bytes[3:arity+2]
    symbol(ascii(atom)), arity + 2
end

function decode (::Type{ErlSmallTuple}, bytes::Array{Uint8})
    arity = bytes[1]
    seq, nbytes, rest = decode_sequence(bytes[2:end], arity)
    tuple(seq...), nbytes + 1
end

function decode (::Type{ErlLargeTuple}, bytes::Array{Uint8})
    arity = decode_integer(bytes[1:4])
    seq, nbytes, rest = decode_sequence(bytes[5:end], arity)
    tuple(seq...), nbytes + 4
end

function decode (::Type{ErlNil}, bytes::Array{Uint8})
    [], 1
end

function decode (::Type{ErlString}, bytes::Array{Uint8})
    arity = decode_integer(bytes[1:2])
    ascii(bytes[3:arity+2]), arity + 2
end

function decode (::Type{ErlList}, bytes::Array{Uint8})
    arity = decode_integer(bytes[1:4])
    seq, nbytes, rest = decode_sequence(bytes[5:end], arity)
    seq, nbytes + 4
end

function decode (::Type{ErlBinary}, bytes::Array{Uint8})
    arity = decode_integer(bytes[1:4])
    bytes[5:arity+4], arity + 4
end

function decode (::Type{ErlSmallBigInt}, bytes::Array{Uint8})
    const b = 256
    arity = bytes[1]
    sign = bytes[2]
    digits = bytes[3:arity+2]

    n = 0
    for i = 1:arity
        n += digits[i] * (b ^ (i - 1))
    end

    BigInt(n), arity + 2
end

function decode (::Type{ErlLargeBigInt}, bytes::Array{Uint8})
    error("Not Implemented")
end

#------------------------------------------------------------------------------
# Utility
#------------------------------------------------------------------------------

function decode_next (bytes::Array{Uint8})
    decoded, nbytes = decode(bytes)
    decoded, bytes[nbytes+1:end], nbytes
end

function decode_sequence (bytes::Array{Uint8}, arity::Integer)
    seq = cell(arity)
    byte_count = 0
    for i = 1:arity
        decoded, next, nbytes = decode_next(bytes)
        bytes = next
        seq[i] = decoded
        byte_count += nbytes
    end
    seq, byte_count, bytes
end

function decode_integer (bytes::Array{Uint8})
    arity = length(bytes)
    bytes = reverse(bytes)
    arity == 1 ? reinterpret(Uint8, bytes)[1] :
    arity == 2 ? reinterpret(Uint16, bytes)[1] :
    arity == 4 ? reinterpret(Uint32, bytes)[1] :
    arity == 8 ? reinterpret(Uint64, bytes)[1] :
    reinterpret(Uint32, bytes)[1]
end

ext: jl

max_stars_repo_name: esl/FsBert

max_issues_repo_name: esl/FsBert

max_forks_repo_name: esl/FsBert


text:
export NALU

"""
    NALU(in::Int, out::Int; initNAC=glorot_uniform, initG=glorot_uniform, initb=glorot_uniform)

Neural Arithmetic Logic Unit. Layer that is capable of learing multiplication,
division, power functions, addition, and subtraction.
As proposed in: https://arxiv.org/abs/1808.00508
"""
struct NALU{Tg<:AbstractMatrix,Tb<:AbstractVector,Tϵ<:Real}
    nac::NAC
    G::Tg
    b::Tb
    ϵ::Tϵ
end

NALU(nac::NAC, G::AbstractMatrix, b::AbstractVector) = NALU(nac, G, b, 1e-8)

function NALU(in::Int, out::Int;
              initNAC=glorot_uniform, initG=glorot_uniform, initb=glorot_uniform)
    nac = NAC(in, out, initW=initNAC, initM=initNAC)
    G = initG(out, in)
    b = initb(out)
    NALU(nac, G, b)
end

add(nalu::NALU, x) = nalu.nac(x)
mult(nalu::NALU, x) = exp.(nalu.nac(log.(abs.(x) .+ nalu.ϵ)))
gate(nalu::NALU, x) = σ.(nalu.G*x .+ nalu.b)

function (nalu::NALU)(x)
    a = add(nalu, x)
    m = mult(nalu, x)
    g = gate(nalu, x)
    g .* a .+ (1 .- g) .* m
end

function Base.show(io::IO, l::NALU)
    in = size(l.G, 2)
    out = size(l.G, 1)
    print(io, "NALU(in=$in, out=$out)")
end

Flux.@functor NALU


text:
loaddict(filename::AbstractString) = Dict(lval => rval for (rval, lval) in split.(eachline(filename), " -> "))
function prettyprint(dict)
    cache = Dict{AbstractString, UInt16}()
    foreach(sort!(collect(keys(dict)))) do key
        println(key, ": ", bobbycalculate(dict, key, cache))
    end
end

bobbycalculate(filename::AbstractString, var) = bobbycalculate(loaddict(filename), var)
function bobbycalculate(dict::Dict, var::AbstractString, cache = Dict{AbstractString, UInt16}())
    get!(cache, var) do
        str = dict[var]
        if all(isdigit, str)
            parse(UInt16, str)
        elseif startswith(str, "NOT")
            ~bobbycalculate(dict, str[5:end], cache)
        elseif str in keys(dict)
            bobbycalculate(dict, str, cache)
        else
            l, i, r = match(r"(.+) (AND|OR|[LR]SHIFT) (.+)", str)
            lval = all(isdigit, l) ? parse(UInt16, l) : bobbycalculate(dict, l, cache)
            instr = if i == "AND"
                (&)
            elseif i == "OR"
                (|)
            elseif i == "LSHIFT"
                (<<)
            elseif i == "RSHIFT"
                (>>)
            end
            rval = all(isdigit, r) ? parse(UInt16, r) : bobbycalculate(dict, r, cache)
            instr(lval, rval)
        end
    end
end

if basename(pwd()) == "aoc"
    cd("2015/7")
end
part1() = bobbycalculate(loaddict("input.txt"), "a")
part2() = bobbycalculate(loaddict("input.txt"), "a", Dict("b"=>0xb3f1))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
@named s = sesam3CN()
@named pl = plant_const()

# @named s = Sesam.sesam3CNC()
# sp = compose(ODESystem([
#     s.i_L ~ pl.i_L,
#     # colimitation
#     s.syn_B ~ s.C_synBC,
#   ], t; name=:sp), s, pl)
#sp = structural_simplify(sp)

@named sp = plant_sesam_system(s,pl)
states(sp)

p = pC = Dict(
    s.ϵ_tvr => 0.45,   # carbon use efficiency of microbial tvr (part by predators
        #which respire and corresponding amount of N must be mineralized)
    s.κ_E => 0.8,     ##<< amount of recycling enzyme turnover by biomass (
        # added to assimilable, i.e. uptake instead of R)
    s.a_E => 0.001*365,   ##<< C biomass allocated to enzymes 1/day /microbial biomass
    s.m => 0.005*365,    ##<< maintenance respiration rate   1/day /microbial biomass,
    s.τ => 1/60*365,  ##<< biomass turnover rate (12 days)
    s.k_L => 5.0,       ##<< 1/(x years)   # formerly 1 year
    s.k_R => 1/(20.0),        ##<< 1/(x years) # to demonstrate changes on short time scale
    s.k_mN_L => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
        # /yr enzyme turnover 60 times a year
    s.k_mN_R => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
    # /yr enzyme turnover 60 times a year
    s.ϵ => 0.5,      ##<< carbon use efficiency for growth respiration
    #i_L => t -> 1 - exp(-t),  # litter input
    #s.ρ_CBtvr => 0.0,    # no carbon resorption
    pl.i_L0 => 400.0,         # g/m2 input per year (half NPP)
    pl.β_Ni0 => 25,
    pl.i_IN0 => 0,   ##<< input of mineral N,
    #
    # P from plant model parameters not used in CN-Sesam soil model
    pl.β_Pi0 => Inf, #25*20, ## leaf litter N:P ~20(massratio Kang10)
    pl.i_IP0 => Inf, #0.65,   ##<< input of mineral P, weathering: Table3 mixed sedimentary rocks 0.65g/m2/yr Hartmann14 10.1016/j.chemgeo.2013.10.025
    pl.s_EP0 => Inf, # 0.5, # plant 1/20 of typical total microbial enzyme synthesis flux
    pl.u_PlantPmax0 => Inf,
    pl.k_PlantP0 => Inf,
)
pN = Dict(
    s.i_BN => 0.4, ##<< potential immobilization flux rate
    s.β_NEnz => 3.1,     # Sterner02: Protein (Fig. 2.2.), high N investment (low P) need
    s.β_NB => 11.0,
    #s.ρ_NBtvr => 0.0,    # no nitrogen resorption
    #s.l_N => 0.96,       #0.00262647*365     ##<< leaching rate of mineralN lN IN
    s.l_N => 0.0,
    s.ν_N =>  0.9,     # microbial N use efficiency accounting for apparent
    ## minceralization of N during uptake in heterogeneous soils
    # for N uptake: take the defaults which take as much N as supplied by litter
    # pl.u_PlantNmax0 => Inf32, # only constrained by k_PlantN in min function
    # pl.k_PlantN0 => 10.57, #0.0289652*365     ##<< plant uptake rate first order of IN
)
p = merge(pC, pN)

u0 = u0C = Dict(
    s.B => 17,
    s.L => 100,
    s.R => 1100,
    #s.cumresp => 0.0,
    s.α_R => 0.1, # enzyme synthesis into L # TODO model by optimality
)
#u0C[s.α_L] = 1.0 - u0C[s.α_R]
u0N = Dict(
    s.I_N => 0.04, ##<< inorganic pool gN/m2
    s.L_N => u0[s.L]/p[pl.β_Ni0],
    s.R_N => u0[s.R]/calculate_β_NR_sesam3(p,s) #p[s.β_NB],
    )
u0 = merge(u0C, u0N)
#u0[s.R]/u0[s.R_N] # smaller p[s.β_NB]

tspan = (0.0,100.0)

#prob = ODEProblem(sp,[t for t in u0], tspan, [t for t in p])
#prob = ODEProblem(sp, remove_units(u0), tspan, remove_units(p))
prob = ODEProblem(sp, u0, tspan, p)
#prob = ODEProblem(sp,u0, tspan, p, jac=true)
sol = sol_sesam3CN = solve(prob);

i_plot = () -> begin
    #using Plots
    plot(sol)
    plot(sol, vars=[s.R])
    Plots.plot(sol, vars=[s.β_NBtvr, p[s.β_NB] * (1-p[s.ρ_CBtvr])/(1-p[s.ρ_NBtvr])], tspan=ts)
    calculate_β_NR_sesam3(p,s), sol[s.β_NR,end]
end

@testset "no N recycling" begin
    #include("test/test_sesam3CN_sol.jl")
    include("test_sesam3CN_sol.jl")
end;

@testset "N recycled" begin
    p2 = copy(p0)
    p2[s.ρ_NBtvr] = 0.1  # 10% N recycled during turnover
    prob2 = ODEProblem(sp, u0, tspan, p2)
    sol = sol2 = solve(prob2, Rodas4());
    #include("test/test_sesam3CN_sol.jl")
    include("test_sesam3CN_sol.jl")
end


text:
@testset "1038.binary-search-tree-to-greater-sum-tree.jl" begin
    @test bst_to_gst(
        TreeNode{Int}([
            4, 1, 6, 0, 2, 5, 7, nothing, nothing, nothing, 3, nothing, nothing, nothing, 8
        ]),
    ) == TreeNode{Int}([
        30,
        36,
        21,
        36,
        35,
        26,
        15,
        nothing,
        nothing,
        nothing,
        33,
        nothing,
        nothing,
        nothing,
        8,
    ])
    @test bst_to_gst(TreeNode{Int}([0, nothing, 1])) == TreeNode{Int}([1, nothing, 1])
    @test bst_to_gst(TreeNode{Int}([1, 0, 2])) == TreeNode{Int}([3, 3, 2])
    @test bst_to_gst(TreeNode{Int}([3, 2, 4, 1])) == TreeNode{Int}([7, 9, 4, 10])
    @test bst_to_gst(nothing) == nothing
end


text:
decode_target(target::Nothing, lost_levels::Bool) = nothing, nothing, nothing
decode_target(target::AbstractArray{Int, 1}, lost_levels::Bool) = decode_target(target..., lost_levels)
function decode_target(target_world::Int, target_stage::Int, lost_levels::Bool)
    if lost_levels
        !(1 ≤ target_world ≤ 12) && error("target_world must be in {1, ..., 12}")
    elseif !(1 ≤ target_world ≤ 12)
        error("target_world must be in {1, ..., 8}")
    end

    !(1 ≤ target_stage ≤ 4) && error("target_stage must be in {1, ..., 4}")

    target_area = target_stage

    if lost_levels
        if target_world ∈ (1, 3)
            if target_stage ≥ 2
                target_area = target_area + 1
            end
        elseif target_world ≥ 5
            worlds = (5, 6, 7, 8, 9, 10, 11, 12)
            error("Lost levels [$worlds] not supported")
        end
    else
        target_world ∈ (1, 2, 4, 7) && target_stage ≥ 2 && (target_area += 1)
    end

    target_world, target_stage, target_area
end


text:
using Pkg

Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate();

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter
using SkipLists

makedocs(
    sitename = "SkipLists",
    format = Documenter.HTML(),
    modules = [SkipLists]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/kernelmethod/SkipLists.jl.git"
)


text:
##############################################################
### copying
##############################################################

import Base.copy
export copy, copy_estimate, GLRM

for T in :[Loss, Regularizer, AbstractGLRM].args
  @eval function copy(r::$T)
    fieldvals = [getfield(r, f) for f in fieldnames(r)]
    return typeof(r)(fieldvals...)
  end
end
# points to all the same problem data as the original input GLRM,
# but copies the estimate of the model parameters
function copy_estimate(g::GLRM)
  return GLRM(g.A,g.losses,g.rx,g.ry,g.k,
              g.observed_features,g.observed_examples,
              copy(g.X),copy(g.Y))
end
# domains are immutable, so this is ok
copy(d::Domain) = d

##############################################################
### fill singleton losses and regularizers to the right shapes
##############################################################

# fill an array of length n with copies of the object foo
fillcopies(foo, n::Int; arraytype=typeof(foo)) = arraytype[copy(foo) for i=1:n]

# singleton loss:
GLRM(A, loss::Loss, rx::Array, ry::Regularizer, k::Int; kwargs...) =
    GLRM(A, fillcopies(loss, size(A, 2), arraytype=Loss), rx, fillcopies(ry, size(A, 2), arraytype=Regularizer), k; kwargs...)
GLRM(A, loss::Loss, rx::Regularizer, ry::Array, k::Int; kwargs...) =
    GLRM(A, fillcopies(loss, size(A, 2), arraytype=Loss), fillcopies(rx, size(A, 1), arraytype=Regularizer), ry, k; kwargs...)
GLRM(A, loss::Loss, rx::Array, ry::Array, k::Int; kwargs...) =
    GLRM(A, fillcopies(loss, size(A, 2), arraytype=Loss), rx, ry, k; kwargs...)

# singleton regularizer on x and/or y:
GLRM(A, losses::Array, rx::Regularizer, ry::Array, k::Int; kwargs...) =
    GLRM(A, losses, fillcopies(rx, size(A, 1), arraytype=Regularizer), ry, k::Int; kwargs...)
GLRM(A, losses::Array, rx::Array, ry::Regularizer, k::Int; kwargs...) =
    GLRM(A, losses, rx, fillcopies(ry, size(A, 2), arraytype=Regularizer), k::Int; kwargs...)
GLRM(A, losses::Array, rx::Regularizer, ry::Regularizer, k::Int; kwargs...) =
    GLRM(A, losses, fillcopies(rx, size(A, 1), arraytype=Regularizer), fillcopies(ry, size(A, 2), arraytype=Regularizer), k::Int; kwargs...)

# singleton everything
GLRM(A, loss::Loss, rx::Regularizer, ry::Regularizer, k::Int; kwargs...) =
    GLRM(A, fillcopies(loss, size(A, 2), arraytype=Loss), fillcopies(rx, size(A, 1), arraytype=Regularizer), fillcopies(ry, size(A, 2), arraytype=Regularizer), k::Int; kwargs...)


text:
__precompile__(true)

module MathKernel

using DocStringExtensions
using StaticArrays
using LinearAlgebra
using Rotations
using ForwardDiff: Dual, value

include("geometry_kernel.jl")
include("poly_approx.jl")
include("rotations.jl")
include("vector_projections.jl")
include("utility.jl")
include("basic_dh.jl")

export
  # geometry_kernel.jl
  area,
  centroid,
  triangleNormal,
  triangle_area,
  vector_area,
  volume,

  # rotations.jl
  quatErr,
  cheapRV,
  components,

  # vector_projections.jl
  vec_proj,
  vec_sub_vec_proj,
  a_dot_one_pad_b,

  # utility.jl
  first_3_of_6,
  last_3_of_6,
  getTop,
  unPad,
  onePad,
  zeroPad,
  mul_then_un_pad,
  one_pad_then_mul,
  weightPoly,
  make_pd_gains,
  rand_pd,

  # basic_dh.jl
  basic_dh,
  dh_R_t,
  povray_12,
  dh_vector_mul

end


text:
"""Cassette Based Counter-factual

```julia
g(x) = @show 2x + 1
h(x) = @show 2x / 3
x_(rng) = g(rand(rng))
Cassette.@Context G2HCtx
Cassette.overdub(::G2HCtx, ::typeof(g), x) = h(x)
x = ciid(x)
x_ = cf(rv, G2HCtx())
y = x_ + x
rand(y)
```
"""
cf(rv::RandVar, ctx::Cassette.Context) = ciid(ω -> Cassette.overdub(ctx, rv, ω))

Cassette.@context SubCtx

function dispatch(ctx, (f_, g_)::Pair{FT, FG}, f::FT, args...) where {FT, FG}
  Cassette.recurse(ctx, g_, args...)
end

dispatch(ctx, metadata, f, args...) = Cassette.recurse(ctx, f, args...)

Cassette.overdub(ctx::SubCtx, f, args...) = dispatch(ctx, ctx.metadata, f, args...)

"""
Cassette based function replacement

`replace(rv::RandVar, function => function)`

```julia
x = normal(0, 1)
noisy(f) = x -> (println("applying f"); f(x))
noisysin = noisy(sin)
y = lift(noisysin)(x)
y_ = replace(y, noisysin => noisy(cos))
rand(y_)
rand(y)
xsampl, ysampl, y_sampl = rand((x, y, y_))
@show cos(xsampl) == y_sampl
```
"""
function Base.replace(rv::RandVar,  replmap::Pair{<:Function, <:Function})
  let ctx = SubCtx(metadata = replmap)
    ciid(ω -> Cassette.overdub(ctx, rv, ω))
  end
end


text:
# ---
# title: 190. Reverse Bits
# id: problem190
# author: Qling
# date: 2021-03-29
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/reverse-bits/description/>
# hidden: true
# ---
#
# Reverse bits of a given 32 bits unsigned integer.
#
# **Note:**
#
#   * Note that in some languages such as Java, there is no unsigned integer type. In this case, both input and output will be given as a signed integer type. They should not affect your implementation, as the integer's internal binary representation is the same, whether it is signed or unsigned.
#   * In Java, the compiler represents the signed integers using [2's complement notation](https://en.wikipedia.org/wiki/Two%27s_complement). Therefore, in **Example 2**  above, the input represents the signed integer `-3` and the output represents the signed integer `-1073741825`.
#
# **Follow up** :
#
# If this function is called many times, how would you optimize it?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 00000010100101000001111010011100
#     Output:    964176192 (00111001011110000010100101000000)
#     Explanation: The input binary string **00000010100101000001111010011100** represents the unsigned integer 43261596, so return 964176192 which its binary representation is **00111001011110000010100101000000**.
#
#
# **Example 2:**
#
#
#
#     Input: n = 11111111111111111111111111111101
#     Output:   3221225471 (10111111111111111111111111111111)
#     Explanation: The input binary string **11111111111111111111111111111101** represents the unsigned integer 4294967293, so return 3221225471 which its binary representation is **10111111111111111111111111111111**.
#
#
#
#
# **Constraints:**
#
#   * The input must be a **binary string** of length `32`
#
#
## @lc code=start
using LeetCode

function reverse_bit(n::UInt32)::UInt32
    ret, power = 0, 31
    while n != 0
        ret += (n & 1) << power
        power -= 1
        n = n >> 1
    end

    return ret
end
## @lc code=end


text:
include("./util.jl")
using Test: @test

instructions = util.load_input(3)

function deliveries(instructions::String)::Set{Tuple{Int,Int}}
    visited = Set()
    curr_pos = (0, 0)
    push!(visited, curr_pos)

    for ins in instructions
        cx, cy = curr_pos
        if ins == '^'
            curr_pos = (cx, cy - 1)
        elseif ins == 'v'
            curr_pos = (cx, cy + 1)
        elseif ins == '>'
            curr_pos = (cx + 1, cy)
        elseif ins == '<'
            curr_pos = (cx - 1, cy)
        else
            throw(DomainError("Unknown instruction $ins"))
        end

        push!(visited, curr_pos)
    end
    visited
end

num_deliveries(instructions) = length(deliveries(instructions))

@test num_deliveries(">") == 2
@test num_deliveries("^>v<") == 4
@test num_deliveries("^v^v^v^v^v") == 2

println("Part 1: $(num_deliveries(instructions))")

normal = ""
robo = ""

for (i, ins) in enumerate(instructions)
    if i % 2 == 0
        global normal *= ins
    else
        global robo *= ins
    end
end

houses = union(deliveries(normal), deliveries(robo))
println("Part 2: $(length(houses))")


text:
using Plots
using SpikingNeuralNetworks
SNN.@load_units

E = SNN.HH(;N = 3200)
I = SNN.HH(;N = 800)
EE = SNN.SpikingSynapse(E, E, :ge; σ = 6nS, p = 0.02)
EI = SNN.SpikingSynapse(E, I, :ge; σ = 6nS, p = 0.02)
IE = SNN.SpikingSynapse(I, E, :gi; σ = 67nS, p = 0.02)
II = SNN.SpikingSynapse(I, I, :gi; σ = 67nS, p = 0.02)
P = [E, I]
C = [EE, EI, IE, II]

SNN.monitor(E, [(:v, [1, 10, 100])])
SNN.sim!(P, C; dt = 0.01ms, duration = 100ms)
SNN.vecplot(E, :v) |> display

using Plots, SNN

RS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -65, d = 8))
IB = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -55, d = 4))
CH = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.2, c = -50, d = 2))
FS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.2, c = -65, d = 2))
TC1 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
TC2 = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.02, b = 0.25, c = -65, d = 0.05))
RZ = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.26, c = -65, d = 2))
LTS = SNN.IZ(;N = 1, param = SNN.IZParameter(;a = 0.1, b = 0.25, c = -65, d = 2))
P = [RS, IB, CH, FS, TC1, TC2, RZ, LTS]

SNN.monitor(P, [:v])
T = 2second
for t = 0:T
    for p in [RS, IB, CH, FS, LTS]
        p.I = [10]
    end
    TC1.I = [(t < 0.2T) ? 0mV : 2mV]
    TC2.I = [(t < 0.2T) ? -30mV : 0mV]
    RZ.I =  [(0.5T < t < 0.6T) ? 10mV : 0mV]
    SNN.sim!(P, [], 0.1ms)
end
SNN.vecplot(P, :v) |> display


text:
mutable struct Ant
    route:: Array{Int64,1}
    currentIndexNode:: Int64
    fitness:: Int64
end

function generateAntPopulation(antCount::Int64,antsStartNode::Int64)
    ants = []
    for i in 1:antCount
        push!(ants,Ant([],antsStartNode,0))
    end
    return ants
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

import Base: trailingsize
import LinearAlgebra.BLAS

const SMALL_FILT_CUTOFF = 58

_zerosi(b,a,T) = zeros(promote_type(eltype(b), eltype(a), T), max(length(a), length(b))-1)

"""
    filt(b, a, x, [si])

Apply filter described by vectors `a` and `b` to vector `x`, with an optional initial filter
state vector `si` (defaults to zeros).
"""
function filt(b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number},
              x::AbstractArray{T}, si::AbstractArray{S} = _zerosi(b,a,T)) where {T,S}
    filt!(Array{promote_type(eltype(b), eltype(a), T, S)}(undef, size(x)), b, a, x, si)
end

# in-place filtering: returns results in the out argument, which may shadow x
# (and does so by default)

"""
    filt!(out, b, a, x, [si])

Same as [`filt`](@ref) but writes the result into the `out` argument, which may
alias the input `x` to modify it in-place.
"""
function filt!(out::AbstractArray, b::Union{AbstractVector, Number}, a::Union{AbstractVector, Number},
               x::AbstractArray{T}, si::AbstractArray{S,N} = _zerosi(b,a,T)) where {T,S,N}
    isempty(b) && throw(ArgumentError("filter vector b must be non-empty"))
    isempty(a) && throw(ArgumentError("filter vector a must be non-empty"))
    a[1] == 0  && throw(ArgumentError("filter vector a[1] must be nonzero"))
    if size(x) != size(out)
        throw(ArgumentError("output size $(size(out)) must match input size $(size(x))"))
    end

    as = length(a)
    bs = length(b)
    sz = max(as, bs)
    silen = sz - 1
    ncols = trailingsize(x,2)

    if size(si, 1) != silen
        throw(ArgumentError("initial state vector si must have max(length(a),length(b))-1 rows"))
    end
    if N > 1 && trailingsize(si,2) != ncols
        throw(ArgumentError("initial state vector si must be a vector or have the same number of columns as x"))
    end

    size(x,1) == 0 && return out
    sz == 1 && return mul!(out, x, b[1]/a[1]) # Simple scaling without memory

    # Filter coefficient normalization
    if a[1] != 1
        norml = a[1]
        a = a ./ norml
        b = b ./ norml
    end
    # Pad the coefficients with zeros if needed
    bs<sz   && (b = copyto!(zeros(eltype(b), sz), b))
    1<as<sz && (a = copyto!(zeros(eltype(a), sz), a))

    initial_si = si
    for col = 1:ncols
        # Reset the filter state
        si = initial_si[:, N > 1 ? col : 1]
        if as > 1
            _filt_iir!(out, b, a, x, si, col)
        elseif bs <= SMALL_FILT_CUTOFF
            _small_filt_fir!(out, b, x, si, col)
        else
            _filt_fir!(out, b, x, si, col)
        end
    end
    return out
end

# Transposed direct form II
function _filt_iir!(out, b, a, x, si, col)
    silen = length(si)
    @inbounds for i=1:size(x, 1)
        xi = x[i,col]
        val = muladd(xi, b[1], si[1])
        for j=1:(silen-1)
            si[j] = muladd(val, -a[j+1], muladd(xi, b[j+1], si[j+1]))
        end
        si[silen] = muladd(xi, b[silen+1], -a[silen+1]*val)
        out[i,col] = val
    end
end

# Transposed direct form II
function _filt_fir!(out, b, x, si, col)
    silen = length(si)
    @inbounds for i=1:size(x, 1)
        xi = x[i,col]
        val = muladd(xi, b[1], si[1])
        for j=1:(silen-1)
            si[j] = muladd(xi, b[j+1], si[j+1])
        end
        si[silen] = b[silen+1]*xi
        out[i,col] = val
    end
end

#
# filt implementation for FIR filters (faster than Base)
#

for n = 2:SMALL_FILT_CUTOFF
    silen = n-1
    si = [Symbol("si$i") for i = 1:silen]
    # Transposed direct form II
    @eval function _filt_fir!(out, b::NTuple{$n,T}, x, siarr, col) where {T}
        offset = (col - 1) * size(x, 1)

        $(Expr(:block, [:(@inbounds $(si[i]) = siarr[$i]) for i = 1:silen]...))
        @inbounds for i=1:size(x, 1)
            xi = x[i+offset]
            val = muladd(xi, b[1], $(si[1]))
            $(Expr(:block, [:($(si[j]) = muladd(xi, b[$(j+1)], $(si[j+1]))) for j = 1:(silen-1)]...))
            $(si[silen]) = b[$(silen+1)]*xi
            out[i+offset] = val
        end
    end
end

# Convert array filter tap input to tuple for small-filtering
let chain = :(throw(ArgumentError("invalid tuple size")))
    for n = SMALL_FILT_CUTOFF:-1:2
        chain = quote
            if length(h) == $n
                _filt_fir!(
                    out,
                    ($([:(@inbounds(h[$i])) for i = 1:n]...),),
                    x,
                    si,
                    col
                )
            else
                $chain
            end
        end
    end

    @eval function _small_filt_fir!(
        out::AbstractArray, h::AbstractVector{T}, x::AbstractArray, si, col
    ) where T
        $chain
    end
end

"""
    deconv(b,a) -> c

Construct vector `c` such that `b = conv(a,c) + r`.
Equivalent to polynomial division.
"""
function deconv(b::StridedVector{T}, a::StridedVector{T}) where T
    lb = size(b,1)
    la = size(a,1)
    if lb < la
        return [zero(T)]
    end
    lx = lb-la+1
    x = zeros(T, lx)
    x[1] = 1
    filt(b, a, x)
end

"""
    _zeropad!(padded::AbstractVector,
              u::AbstractVector,
              padded_axes = axes(padded),
              data_dest::Tuple = (1,),
              data_region = CartesianIndices(u))

Place the portion of `u` specified by `data_region` into `padded`, starting at
location `data_dest`. Sets the rest of `padded` to zero. This will mutate
`padded`. `padded_axes` must correspond to the axes of `padded`.

"""
@inline function _zeropad!(
    padded::AbstractVector,
    u::AbstractVector,
    padded_axes = axes(padded),
    data_dest::Tuple = (first(padded_axes[1]),),
    data_region = CartesianIndices(u),
)
    datasize = length(data_region)
    # Use axes to accommodate arrays that do not start at index 1
    data_first_i = first(data_region)[1]
    dest_first_i = data_dest[1]
    copyto!(padded, dest_first_i, u, data_first_i, datasize)
    padded[first(padded_axes[1]):dest_first_i - 1] .= 0
    padded[dest_first_i + datasize : end] .= 0

    padded
end

@inline function _zeropad!(
    padded::AbstractArray,
    u::AbstractArray,
    padded_axes = axes(padded),
    data_dest::Tuple = first.(padded_axes),
    data_region = CartesianIndices(u),
)
    fill!(padded, zero(eltype(padded)))
    dest_axes = UnitRange.(data_dest, data_dest .+ size(data_region) .- 1)
    dest_region = CartesianIndices(dest_axes)
    copyto!(padded, dest_region, u, data_region)

    padded
end

"""
    _zeropad(u, padded_size, [data_dest, data_region])

Creates and returns a new base-1 index array of size `padded_size`, with the
section of `u` specified by `data_region` copied into the region of the new
 array as specified by `data_dest`. All other values will be initialized to
 zero.

If either `data_dest` or `data_region` is not specified, then the defaults
described in [`_zeropad!`](@ref) will be used.
"""
function _zeropad(u, padded_size, args...)
    padded = similar(u, padded_size)
    _zeropad!(padded, u, axes(padded), args...)
end

function _zeropad_keep_offset(u, padded_size, u_axes, args...)
    ax_starts = first.(u_axes)
    new_axes = UnitRange.(ax_starts, ax_starts .+ padded_size .- 1)
    _zeropad!(similar(u, new_axes), u, args...)
end

function _zeropad_keep_offset(
    u, padded_size, ::NTuple{<:Any, Base.OneTo{Int}}, args...
)
    _zeropad(u, padded_size, args...)
end

"""
    _zeropad_keep_offset(u, padded_size, [data_dest, dat_region])

Like [`_zeropad`](@ref), but retains the first index of `u` when creating a new
array.
"""
function _zeropad_keep_offset(u, padded_size, args...)
    _zeropad_keep_offset(u, padded_size, axes(u), args...)
end

"""
Estimate the number of floating point multiplications per output sample for an
overlap-save algorithm with fft size `nfft`, and filter size `nb`.
"""
os_fft_complexity(nfft, nb) =  (nfft * log2(nfft) + nfft) / (nfft - nb + 1)

"""
Determine the length of FFT that minimizes the number of multiplications per
output sample for an overlap-save convolution of vectors of size `nb` and `nx`.
"""
function optimalfftfiltlength(nb, nx)
    nfull = nb + nx - 1

    # Step through possible nffts and find the nfft that minimizes complexity
    # Assumes that complexity is convex
    first_pow2 = ceil(Int, log2(nb))
    max_pow2 = ceil(Int, log2(nfull))
    prev_complexity = os_fft_complexity(2 ^ first_pow2, nb)
    pow2 = first_pow2 + 1
    while pow2 <= max_pow2
        new_complexity = os_fft_complexity(2 ^ pow2, nb)
        new_complexity > prev_complexity && break
        prev_complexity = new_complexity
        pow2 += 1
    end
    nfft = pow2 > max_pow2 ? 2 ^ max_pow2 : 2 ^ (pow2 - 1)

    if nfft > nfull
        # If nfft > nfull, then it's better to find next fast power
        nfft = nextfastfft(nfull)
    end

    nfft
end

"""
Prepare buffers and FFTW plans for convolution. The two buffers, tdbuff and
fdbuff may be an alias of each other, because complex convolution only requires
one buffer. The plans are mutating where possible, and the inverse plan is
unnormalized.
"""
@inline function os_prepare_conv(u::AbstractArray{T, N},
                                 nffts) where {T<:Real, N}
    tdbuff = similar(u, nffts)
    bufsize = ntuple(i -> i == 1 ? nffts[i] >> 1 + 1 : nffts[i], N)
    fdbuff = similar(u, Complex{T}, NTuple{N, Int}(bufsize))

    p = plan_rfft(tdbuff)
    ip = plan_brfft(fdbuff, nffts[1])

    tdbuff, fdbuff, p, ip
end

@inline function os_prepare_conv(u::AbstractArray{<:Complex}, nffts)
    buff = similar(u, nffts)

    p = plan_fft!(buff)
    ip = plan_bfft!(buff)

    buff, buff, p, ip # Only one buffer for complex
end

"""
Transform the smaller convolution input to frequency domain, and return it in a
new array. However, the contents of `buff` may be modified.
"""
@inline function os_filter_transform!(buff::AbstractArray{<:Real}, p)
    p * buff
end

@inline function os_filter_transform!(buff::AbstractArray{<:Complex}, p!)
    copy(p! * buff) # p operates in place on buff
end

"""
Take a block of data, and convolve it with the smaller convolution input. This
may modify the contents of `tdbuff` and `fdbuff`, and the result will be in
`tdbuff`.
"""
@inline function os_conv_block!(tdbuff::AbstractArray{<:Real},
                                fdbuff::AbstractArray,
                                filter_fd,
                                p,
                                ip)
    mul!(fdbuff, p, tdbuff)
    fdbuff .*= filter_fd
    mul!(tdbuff, ip, fdbuff)
end

"Like the real version, but only operates on one buffer"
@inline function os_conv_block!(buff::AbstractArray{<:Complex},
                                ::AbstractArray, # Only one buffer for complex
                                filter_fd,
                                p!,
                                ip!)
    p! * buff # p! operates in place on buff
    buff .*= filter_fd
    ip! * buff # ip! operates in place on buff
end

# Used by `unsafe_conv_kern_os!` to handle blocks of input data that need to be padded.
#
# For a given number of edge dimensions, convolve all regions along the
# perimeter that have that number of edge dimensions
#
# For a 3d cube, if n_edges = 1, this means all faces. If n_edges = 2, then
# all edges. Finally, if n_edges = 3, then all corners.
#
# This needs to be a separate function for subsets to generate tuple elements,
# which is only the case if the number of edges is a Val{n} type. Iterating over
# the number of edges with Val{n} is inherently type unstable, so this function
# boundary allows dispatch to make efficient code for each number of edge
# dimensions.
function unsafe_conv_kern_os_edge!(
    # These arrays and buffers will be mutated
    out::AbstractArray{<:Any, N},
    tdbuff,
    fdbuff,
    perimeter_range,
    # Number of edge dimensions to pad and convolve
    n_edges::Val,
    # Data to be convolved
    u,
    filter_fd,
    # FFTW plans
    p,
    ip,
    # Sizes, ranges, and other pre-calculated constants
    #
    ## ranges listing center and edge blocks
    edge_ranges,
    center_block_ranges,
    ## size and axis information
    all_dims, # 1:N
    su,
    u_start,
    sv,
    nffts,
    out_start,
    out_stop,
    save_blocksize,
    sout_deficit, # How many output samples are missing if nffts > sout
    tdbuff_axes,
) where N
    # Iterate over all possible combinations of edge dimensions for a number of
    # edges
    #
    # For a 3d cube with n_edges = 1, this will specify the top and bottom faces
    # (edge_dims = (1,)), then the left and right faces (edge_dims = (2,)), then
    # the front and back faces (edge_dims = (3,))
    for edge_dims in subsets(all_dims, n_edges)
        # Specify a region on the perimeter by combining an edge block index for
        # each dimension on an edge, and the central blocks for dimensions not
        # on an edge.
        #
        # First make all entries equal to the center blocks:
        @inbounds copyto!(perimeter_range, 1, center_block_ranges, 1, N)

        # For the dimensions chosen to be on an edge (edge_dims), get the
        # ranges of the blocks that would need to be padded (lie on an edge)
        # in that dimension.
        #
        # There can be one to two such ranges for each dimension, because with
        # some inputs sizes the whole convolution is just one range
        # (one edge block), or the padding will be needed at both the leading
        # and trailing side of that dimension
        selected_edge_ranges = getindex.(Ref(edge_ranges), edge_dims)

        # Visit each combination of edge ranges for the edge dimensions chosen.
        # For a 3d cube with n_edges = 1 and edge_dims = (1,), this will visit
        # the top face, and then the bottom face.
        for perimeter_edge_ranges in Iterators.ProductIterator(selected_edge_ranges)
            # The center region for non-edge dimensions has been specified above,
            # so finish specifying the region of the perimeter for this edge
            # block
            @inbounds for (i, dim) in enumerate(edge_dims)
                perimeter_range[dim] = perimeter_edge_ranges[i]
            end

            # Region of block indices, not data indices!
            block_region = CartesianIndices(
                NTuple{N, UnitRange{Int}}(perimeter_range)
            )
            @inbounds for block_pos in block_region
                # Figure out which portion of the input data should be transformed

                block_idx = convert(NTuple{N, Int}, block_pos)
                ## data_offset is NOT the beginning of the region that will be
                ## convolved, but is instead the beginning of the unaliased data.
                data_offset = save_blocksize .* (block_idx .- 1)
                ## How many zeros will need to be added before the data
                pad_before = max.(0, sv .- data_offset .- 1)
                data_ideal_stop = data_offset .+ save_blocksize
                ## How many zeros will need to be added after the data
                pad_after = max.(0, data_ideal_stop .- su)

                ## Data indices, not block indices
                data_region = CartesianIndices(
                    UnitRange.(
                        u_start .+ data_offset .- sv .+ pad_before .+ 1,
                        u_start .+ data_ideal_stop .- pad_after .- 1
                    )
                )

                # Convolve portion of input

                _zeropad!(tdbuff, u, tdbuff_axes, pad_before .+ 1, data_region)
                os_conv_block!(tdbuff, fdbuff, filter_fd, p, ip)

                # Save convolved result to output

                block_out_stop = min.(
                    out_start .+ data_offset .+ save_blocksize .- 1,
                    out_stop
                )
                block_out_region = CartesianIndices(
                    UnitRange.(out_start .+ data_offset, block_out_stop)
                )
                ## If the input could not fill tdbuff, account for that before
                ## copying the convolution result to the output
                u_deficit = max.(0, pad_after .- sv .+ 1)
                valid_buff_region = CartesianIndices(
                    UnitRange.(sv, nffts .- u_deficit .- sout_deficit)
                )
                copyto!(out, block_out_region, tdbuff, valid_buff_region)
            end
        end
    end
end

# Assumes u is larger than, or the same size as, v
# nfft should be greater than or equal to 2*sv-1
function unsafe_conv_kern_os!(out,
                        u::AbstractArray{<:Any, N},
                        v,
                        su,
                        sv,
                        sout,
                        nffts) where N
    u_start = first.(axes(u))
    out_axes = axes(out)
    out_start = first.(out_axes)
    out_stop = last.(out_axes)
    ideal_save_blocksize = nffts .- sv .+ 1
    # Number of samples that are "missing" if the output is smaller than the
    # valid portion of the convolution
    sout_deficit = max.(0, ideal_save_blocksize .- sout)
    # Size of the valid portion of the convolution result
    save_blocksize = ideal_save_blocksize .- sout_deficit
    nblocks = cld.(sout, save_blocksize)

    # Pre-allocation
    tdbuff, fdbuff, p, ip = os_prepare_conv(u, nffts)
    tdbuff_axes = axes(tdbuff)

    # Transform the smaller filter
    _zeropad!(tdbuff, v)
    filter_fd = os_filter_transform!(tdbuff, p)
    filter_fd .*= 1 / prod(nffts) # Normalize once for brfft

    # block indices for center blocks, which need no padding
    first_center_blocks = cld.(sv .- 1, save_blocksize) .+ 1
    last_center_blocks = fld.(su, save_blocksize)
    center_block_ranges = UnitRange.(first_center_blocks, last_center_blocks)

    # block index ranges for blocks that need to be padded
    # Corresponds to the leading and trailing side of a dimension, or if there
    # are no center blocks, corresponds to the whole dimension
    edge_ranges = map(nblocks, first_center_blocks, last_center_blocks) do nblock, firstfull, lastfull
        lastfull > 1 ? [1:firstfull - 1, lastfull + 1 : nblock] : [1:nblock]
    end
    all_dims = 1:N
    # Buffer to store ranges of indices for a single region of the perimeter
    perimeter_range = Vector{UnitRange{Int}}(undef, N)

    # Convolve all blocks that require padding.
    #
    # This is accomplished by dividing the perimeter of the volume into
    # subsections, where the division is done by the number of edge dimensions.
    # For a 3d cube, this convolves the perimeter in the following order:
    #
    # Number of Edge Dimensions | Convolved Region
    # --------------------------+-----------------
    #                         1 | Faces of Cube
    #                         2 | Edges of Cube
    #                         3 | Corners of Cube
    #
    for n_edges in all_dims
        unsafe_conv_kern_os_edge!(
            # These arrays and buffers will be mutated
            out,
            tdbuff,
            fdbuff,
            perimeter_range,
            # Number of edge dimensions to pad and convolve
            Val(n_edges),
            # Data to be convolved
            u,
            filter_fd,
            # FFTW plans
            p,
            ip,
            # Sizes, ranges, and other pre-calculated constants
            #
            ## ranges listing center and edge blocks
            edge_ranges,
            center_block_ranges,
            ## size and axis information
            all_dims, # 1:N
            su,
            u_start,
            sv,
            nffts,
            out_start,
            out_stop,
            save_blocksize,
            sout_deficit,
            tdbuff_axes) # How many output samples are missing if nffts > sout
    end

    tdbuff_region = CartesianIndices(tdbuff)
    # Portion of buffer with valid result of convolution
    valid_buff_region = CartesianIndices(UnitRange.(sv, nffts))
    # Iterate over block indices (not data indices) that do not need to be padded
    @inbounds for block_pos in CartesianIndices(center_block_ranges)
        # Calculate portion of data to transform

        block_idx = convert(NTuple{N, Int}, block_pos)
        ## data_offset is NOT the beginning of the region that will be
        ## convolved, but is instead the beginning of the unaliased data.
        data_offset = save_blocksize .* (block_idx .- 1)
        data_stop = data_offset .+ save_blocksize
        data_region = CartesianIndices(
            UnitRange.(u_start .+ data_offset .- sv .+ 1, u_start .+ data_stop .- 1)
        )

        # Convolve this portion of the data

        copyto!(tdbuff, tdbuff_region, u, data_region)
        os_conv_block!(tdbuff, fdbuff, filter_fd, p, ip)

        # Save convolved result to output

        block_out_region = CartesianIndices(
            UnitRange.(data_offset .+ out_start, data_stop .+ out_start .- 1)
        )
        copyto!(out, block_out_region, tdbuff, valid_buff_region)
    end

    out
end

function _conv_kern_fft!(out,
                         u::AbstractArray{T, N},
                         v::AbstractArray{T, N},
                         su,
                         sv,
                         outsize,
                         nffts) where {T<:Real, N}
    padded = _zeropad(u, nffts)
    p = plan_rfft(padded)
    uf = p * padded
    _zeropad!(padded, v)
    vf = p * padded
    uf .*= vf
    raw_out = irfft(uf, nffts[1])
    copyto!(out,
            CartesianIndices(out),
            raw_out,
            CartesianIndices(UnitRange.(1, outsize)))
end
function _conv_kern_fft!(out, u, v, su, sv, outsize, nffts)
    upad = _zeropad(u, nffts)
    vpad = _zeropad(v, nffts)
    p! = plan_fft!(upad)
    p! * upad # Operates in place on upad
    p! * vpad
    upad .*= vpad
    ifft!(upad)
    copyto!(out,
            CartesianIndices(out),
            upad,
            CartesianIndices(UnitRange.(1, outsize)))
end

# v should be smaller than u for good performance
function _conv_fft!(out, u, v, su, sv, outsize)
    os_nffts = map(optimalfftfiltlength, sv, su)
    if any(os_nffts .< outsize)
        unsafe_conv_kern_os!(out, u, v, su, sv, outsize, os_nffts)
    else
        nffts = nextfastfft(outsize)
        _conv_kern_fft!(out, u, v, su, sv, outsize, nffts)
    end
end

# For arrays with weird offsets
function _conv_similar(u, outsize, axesu, axesv)
    out_offsets = first.(axesu) .+ first.(axesv)
    out_axes = UnitRange.(out_offsets, out_offsets .+ outsize .- 1)
    similar(u, out_axes)
end
function _conv_similar(
    u, outsize, ::NTuple{<:Any, Base.OneTo{Int}}, ::NTuple{<:Any, Base.OneTo{Int}}
)
    similar(u, outsize)
end
_conv_similar(u, v, outsize) = _conv_similar(u, outsize, axes(u), axes(v))

# Does convolution, will not switch argument order
function _conv!(out, u, v, su, sv, outsize)
    # TODO: Add spatial / time domain algorithm
    _conv_fft!(out, u, v, su, sv, outsize)
end

# Does convolution, will not switch argument order
function _conv(u, v, su, sv)
    outsize = su .+ sv .- 1
    out = _conv_similar(u, v, outsize)
    _conv!(out, u, v, su, sv, outsize)
end

# May switch argument order
"""
    conv(u,v)

Convolution of two arrays. Uses either FFT convolution or overlap-save,
depending on the size of the input. `u` and `v` can be  N-dimensional arrays,
with arbitrary indexing offsets, but their axes must be a `UnitRange`.
"""
function conv(u::AbstractArray{T, N},
              v::AbstractArray{T, N}) where {T<:BLAS.BlasFloat, N}
    su = size(u)
    sv = size(v)
    if prod(su) >= prod(sv)
        _conv(u, v, su, sv)
    else
        _conv(v, u, sv, su)
    end
end

function conv(u::AbstractArray{<:BLAS.BlasFloat, N},
              v::AbstractArray{<:BLAS.BlasFloat, N}) where N
    fu, fv = promote(u, v)
    conv(fu, fv)
end

conv(u::AbstractArray{<:Integer, N}, v::AbstractArray{<:Integer, N}) where {N} =
    round.(Int, conv(float(u), float(v)))

conv(u::AbstractArray{<:Number, N}, v::AbstractArray{<:Number, N}) where {N} =
    conv(float(u), float(v))

function conv(u::AbstractArray{<:Number, N},
              v::AbstractArray{<:BLAS.BlasFloat, N}) where N
    conv(float(u), v)
end

function conv(u::AbstractArray{<:BLAS.BlasFloat, N},
              v::AbstractArray{<:Number, N}) where N
    conv(u, float(v))
end

function conv(A::AbstractArray{<:Number, M},
              B::AbstractArray{<:Number, N}) where {M, N}
    if (M < N)
        conv(cat(A, dims=N)::AbstractArray{eltype(A), N}, B)
    else
        @assert M > N
        conv(A, cat(B, dims=M)::AbstractArray{eltype(B), M})
    end
end

"""
    conv(u,v,A)

2-D convolution of the matrix `A` with the 2-D separable kernel generated by
the vectors `u` and `v`.
Uses 2-D FFT algorithm.
"""
function conv(u::AbstractVector{T}, v::AbstractVector{T}, A::AbstractMatrix{T}) where T
    # Arbitrary indexing offsets not implemented
    @assert !Base.has_offset_axes(u, v, A)
    m = length(u)+size(A,1)-1
    n = length(v)+size(A,2)-1
    B = zeros(T, m, n)
    B[1:size(A,1),1:size(A,2)] = A
    u = fft([u;zeros(T,m-length(u))])
    v = fft([v;zeros(T,n-length(v))])
    C = ifft(fft(B) .* (u * transpose(v)))
    if T <: Real
        return real(C)
    end
    return C
end

dsp_reverse(v, ::NTuple{<:Any, Base.OneTo{Int}}) = reverse(v, dims = 1)
function dsp_reverse(v, vaxes)
    vsize = length(v)
    reflected_start = - first(vaxes[1]) - vsize + 1
    reflected_axes = (reflected_start : reflected_start + vsize - 1,)
    out = similar(v, reflected_axes)
    copyto!(out, reflected_start, Iterators.reverse(v), 1, vsize)
end

"""
    xcorr(u,v; padmode = :longest)

Compute the cross-correlation of two vectors, by calculating the similarity
between `u` and `v` with various offsets of `v`. Delaying `u` relative to `v`
will shift the result to the right.

The size of the output depends on the padmode keyword argument: with padmode =
:none the length of the result will be length(u) + length(v) - 1, as with conv.
With padmode = :longest the shorter of the arguments will be padded so they are
equal length. This gives a result with length 2*max(length(u), length(v))-1,
with the zero-lag condition at the center.
"""
function xcorr(
    u::AbstractVector, v::AbstractVector; padmode::Symbol = :none
)
    su = size(u,1); sv = size(v,1)
    if padmode == :longest
        if su < sv
            u = _zeropad_keep_offset(u, sv)
        elseif sv < su
            v = _zeropad_keep_offset(v, su)
        end
        conv(u, dsp_reverse(conj(v), axes(v)))
    elseif padmode == :none
        conv(u, dsp_reverse(conj(v), axes(v)))
    else
        throw(ArgumentError("padmode keyword argument must be either :none or :longest"))
    end
end


text:
_vsl_errors = Dict()
macro vsl_error(code)
    _vsl_errors[code.args[2]] = code.args[1]
    return esc(Meta.parse("const $code"))
end
struct VSLError <: Exception
    errno :: Int32
    message :: AbstractString
end
show(io::IO, E::VSLError) = print(io, string("VSL Error: ", strerror(E.errno), " -- ", E.message))
VSL_ERROR(errno::Integer) = VSLError(errno, _vsl_errors[errno])

@vsl_error VSL_STATUS_OK = 0

@vsl_error VSL_ERROR_FEATURE_NOT_IMPLEMENTED = -1
@vsl_error VSL_ERROR_UNKNOWN                 = -2
@vsl_error VSL_ERROR_BADARGS                 = -3
@vsl_error VSL_ERROR_MEM_FAILURE             = -4
@vsl_error VSL_ERROR_NULL_PTR                = -5
@vsl_error VSL_ERROR_CPU_NOT_SUPPORTED       = -6

@vsl_error VSL_RNG_ERROR_INVALID_BRNG_INDEX        = -1000
@vsl_error VSL_RNG_ERROR_LEAPFROG_UNSUPPORTED      = -1002
@vsl_error VSL_RNG_ERROR_SKIPAHEAD_UNSUPPORTED     = -1003
@vsl_error VSL_RNG_ERROR_BRNGS_INCOMPATIBLE        = -1005
@vsl_error VSL_RNG_ERROR_BAD_STREAM                = -1006
@vsl_error VSL_RNG_ERROR_BRNG_TABLE_FULL           = -1007
@vsl_error VSL_RNG_ERROR_BAD_STREAM_STATE_SIZE     = -1008
@vsl_error VSL_RNG_ERROR_BAD_WORD_SIZE             = -1009
@vsl_error VSL_RNG_ERROR_BAD_NSEEDS                = -1010
@vsl_error VSL_RNG_ERROR_BAD_NBITS                 = -1011
@vsl_error VSL_RNG_ERROR_QRNG_PERIOD_ELAPSED       = -1012
@vsl_error VSL_RNG_ERROR_LEAPFROG_NSTREAMS_TOO_BIG = -1013
@vsl_error VSL_RNG_ERROR_BRNG_NOT_SUPPORTED        = -1014

@vsl_error VSL_RNG_ERROR_BAD_UPDATE              = -1120
@vsl_error VSL_RNG_ERROR_NO_NUMBERS              = -1121
@vsl_error VSL_RNG_ERROR_INVALID_ABSTRACT_STREAM = -1122

@vsl_error VSL_RNG_ERROR_NONDETERM_NOT_SUPPORTED     = -1130
@vsl_error VSL_RNG_ERROR_NONDETERM_NRETRIES_EXCEEDED = -1131

@vsl_error VSL_RNG_ERROR_ARS5_NOT_SUPPORTED = -1140

@vsl_error VSL_RNG_ERROR_FILE_CLOSE = -1100
@vsl_error VSL_RNG_ERROR_FILE_OPEN  = -1101
@vsl_error VSL_RNG_ERROR_FILE_WRITE = -1102
@vsl_error VSL_RNG_ERROR_FILE_READ  = -1103

@vsl_error VSL_RNG_ERROR_BAD_FILE_FORMAT      = -1110
@vsl_error VSL_RNG_ERROR_UNSUPPORTED_FILE_VER = -1111

@vsl_error VSL_RNG_ERROR_BAD_MEM_FORMAT = -1200


text:
""" Conversion pipe dreams... """
function bbr_to_bbc!(gm::GlobalModel, idx::Array{Int64} = collect(1:length(gm.bbls)))
    for i in idx
        if gm.bbls[i] isa BlackBoxRegressor
            bbr = gm.bbls[i]
            if bbr.constraint isa JuMP.ConstraintRef
                throw(OCTException("Cannot convert BBR $(bbr.name) to a BBC " *
                    "since the bbr.constraint isa JuMP.ConstraintRef."))
            end
            if !isempty(bbr.mi_constraints)
                clear_tree_constraints!(gm, bbr)
            end
            clear_data!(bbr)
            new_bbc = BlackBoxClassifier(constraint = bbr.constraint,
                                         vars = [bbr.vars..., bbr.dependent_var],
                                         expr_vars = [bbr.expr_vars..., bbr.dependent_var],
                                         equality = bbr.equality, name = bbr.name)
            new_bbc.f(new_bbc.expr_vars...) = bbr.f()
            f = functionify(expr)
            gradable_fn = x -> Base.invokelatest(new_bbc.f, [x[i] for i in var_ranges]...)
            return x -> ForwardDiff.gradient(gradable_fn, x)

            if get_param(bbr, :gradients)
                var_ranges = get_var_ranges(new_bbc.expr_vars)
            f = functionify(expr)
            gradable_fn = x -> Base.invokelatest(new_bbc.f, [x[i] for i in var_ranges]...)
            return x -> ForwardDiff.gradient(gradable_fn, x)

            expr_vars = vars
            if isnothing(dependent_var)
                new_bbl = BlackBoxClassifier(constraint = constraint, vars = vars, expr_vars = expr_vars,
                                                equality = equality, name = name)
                add_data!(new_bbl, X, Y)
                push!(gm.bbls, new_bbl)
                return
            end
        end
    end
end


text:
using VectorEngine, Test

mutable struct xm
      x::Int32
      m::Int64
end

function pass_struct!(r::xm)
    r.x -= 2
    r.m += 1
    return
end

a = xm(1, 100)

@veda pass_struct!(a)
synchronize()
@test a.x == -1 && a.m == 101


text:
include("runtests_preamble.jl")

@testset "LESbrary" begin
    include("test_fit_cubic.jl")
    include("test_interpolated_profiles.jl")
    include("test_diagnose_buoyancy_flux.jl")
    include("test_turbulence_statistics.jl")
    include("test_examples.jl")
    include("test_experiments.jl")
end

max_stars_repo_name: ali-ramadhan/LESbrary

max_issues_repo_name: ali-ramadhan/LESbrary

max_forks_repo_name: ali-ramadhan/LESbrary


text:
arg_dict = Dict{Any, Any}(
  # "operator_type"   => "SBPGamma",
  "operator_type"   => "SBPOmega",
  # "operator_type"   => "SBPDiagE",
  "physics"         => "Elliptic",
  "use_DG"          => true,
  "delta_t"         => 1.0e-3,
  "t_max"           => 5.0,
  "res_tol"         => 1e-12,
  "step_tol"        => 1e-12,
  "itermax"         => 1,
  "order"           => 2,
  "Relfunc_name"    => true,
  "smb_name"        => "SRCMESHES/square16x16.smb",
  "dmg_name"        => ".null",
  "numBC"           => 1,
  "BC1"             => [0, 1, 2, 3],
  "exactSolution"   =>"ExactExpTrig",
  "BC1_name"        => "DirichletTrig",
  "Diffusion"       => "poly0th",
  "SRC_name"        => "SrcExpTrigPoly0thDiffn",
  "Functional"      => "volumeAverage",
  "exactFunctional" => 1.846230857168755189823348538541e-02,
  "IC_name"         => "ICRandom",

  # "exactSolution"   =>"ExactTrig",
  # "BC1_name"        => "DirichletTrig",
  # "Diffusion"       => "poly2nd", # lambda = [x^2+1 & xy \\ xy & y^2+1]
  # "SRC_name"        => "SrcTrigPoly2ndDiffn",

  # "exactSolution"   =>"ExactTrig",
  # "BC1_name"        => "DirichletTrig",
  # "Diffusion"       => "poly6th",
  # "SRC_name"        => "SrcTrigPoly6thDiffn",

  # "exactSolution"   =>"ExactPoly2nd",
  # "Diffusion"       => "DiffnPoly0th"
  # "BC1_name"        => "DirichletPolynial2nd",
  # "SRC_name"        => "SrcPoly2nd",

  # "SRC_name"        => "SRC0",		# for energy stability test

  # "Flux_name"       => "SAT",			# SAT-BR2
  "Flux_name"       => "SAT0",		# SAT-SIPG
  # "Flux_name"       => "Shahbazi",	# Shahbazi
  "Cip"             => -3.0,
  "run_type"        => 5,
  "jac_type"        => 2,
  "jac_method"      => 2,
  "write_eigs"      => false,
  "write_jac"       => false,
  "write_conditioning" => false,
  "write_energy"    => true,
  "solve"           => true,
  "real_time"       => false
)


text:
function cost_no_Etot(L::Array{Float64, 1}, N::Array{Float64, 1}, Ltot::Float64, γ::Float64)
    #### Returns a Float64. It is the cost function evaluated in γ (in the case where there is no constrain on the energy) ####
    # getting the constants we will need
    n_pop = length(L)

    # Let's start the integration
    int = 0.0
    for a = 1:n_pop
        Lₐ = L[a]
        nₐ = N[a]

        x = γ * Lₐ

        int += nₐ * (Lₐ * coth(x) - 1 / γ)
    end

    return(Ltot - int)
end

function Dichotomy_step(L::Array{Float64, 1}, N::Array{Float64, 1}, Ltot::Float64, γ_min::Float64, γ_max::Float64)
    #### Returns a tuple of 2 floats. They are the new upper and lower bounds at the end of the dichotomy step ####
    # Let's get the current middle of the segment and the value of the cost function there
    γ_mid = (γ_max + γ_min) / 2
    c_mid = cost_no_Etot(L, N, Ltot, γ_mid)

    # Let's get the new bounds
    if (c_mid > 0) # We use the fact that the cost function is monotonous and decreasing
        return (γ_mid, γ_max)
    elseif (c_mid < 0)
        return (γ_min, γ_mid)
    else
        return (γ_mid, γ_mid)
    end
end

function Solution_no_Etot_dichotomy(L::Array{Float64, 1}, N::Array{Float64, 1}, Ltot::Float64)
    #### Returns a Float64. It is the optimal γ found by the dichotomy method ####
    #### This function performs the whole dichotomy                           ####

    # first let's make define the starting bounds of the dichotomy
    min_L, useless = findmin(L)
    max_L, useless = findmax(L)

    γ_min = 10.0^(-7) / max_L # these bounds come from looking at the shape of coth(x) - 1/x
    γ_max = 10.0^(2) / min_L
    γ_mid = (γ_max + γ_min) / 2

    # then let's define an acceptable relative error and maximum number of steps
    ϵ = 10^(-6)
    n_step = 100

    # then let's do the dichotomy method itself
    i = 0
    cost = Inf
    while ((cost > ϵ) && (i <= n_step))
        γ_min, γ_max = Dichotomy_step(L, N, Ltot, γ_min, γ_max) # update the bounds

        γ_mid = (γ_max + γ_min) / 2
        cost = abs(cost_no_Etot(L, N, Ltot, γ_mid)) # update the cost and test wether it is acceptable

        i += 1
    end

    return(γ_mid)
end

function Magnetisations_no_Etot(L::Array{Float64, 1}, N::Array{Float64, 1}, γ::Float64,
                                lmax::Int64,
                                List_x::Array{Float64,1}, List_w::Array{Float64,1}, List_Y::Array{Float64, 2})
    #### Returns an array of size n_pop $ lmax_half. Each cell contains one (axisymmetric, even) magnetisation                               ####
    #### This function computes the magnetization from a solution given by the dichotomy method, but doesn't run the dichotomy method itself ####
    # getting the constants we will need
    n_pop = length(L)
    lmax_half = div(lmax, 2)
    n_quad = length(List_x)

    # Preallocating the arrays we will need
    M = zeros(n_pop * lmax_half) # Will hold the magnetisations

    # Let's start to compute the magnetisations
    for star = 1:n_pop
        Lₛ = L[star]
        nₛ = N[star]

        # first let's compute the renrmalizing factor (lower integral)
        int_down = 0.0
        for i = 1:n_quad
            weight = List_w[i]
            x = List_x[i]

            int_down += weight * exp(γ * Lₛ * x)
        end

        # then let's compute the upper integral
        for l_half = 1:lmax_half
            index = (star - 1) * lmax_half + l_half

            for i = 1:n_quad
                spherical_harmonic = List_Y[l_half, i]
                weight = List_w[i]
                x = List_x[i]

                M[index] += weight * nₛ * spherical_harmonic * exp(γ * Lₛ * x)
            end

            M[index] /= int_down
        end
    end

    return(M)
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
# ======================= SHOWING ====================

function Base.show(io::IO, x::Pattern)
    expr = to_expr(x)
    print(io, expr)
end

to_expr(x::PatVar) = x.name

to_expr(x::PatLiteral) =
    if x.val isa Symbol
        QuoteNode(x.val)
    else
        x.val
    end

function to_expr(x::PatTypeAssertion)
    Expr(Symbol("::"), to_expr(x.var), x.type)
end

function to_expr(x::PatSplatVar)
    Expr(Symbol("..."), to_expr(x.var))
end

function to_expr(x::PatEquiv)
    Expr(:call, Symbol("≡ₙ"), to_expr(x.left), to_expr(x.right))
end

function to_expr(x::PatTerm)
    pl = operation(x)
    similarterm(Expr, pl, arguments(x); exprhead=exprhead(x))
end

function to_expr(x::PatAllTerm)
    # TODO change me ?
    head = Symbol("~", x.head.name)
    Expr(:call, head, to_expr.(x.args)...)
end


text:
"""
	$(SIGNATURES)

Stacked bar graphs for a matrix by (x,y,z).
Dimension `subPlotIdx` is used for subplots. Dimension `groupIdx` is used to group bars. The remaining dimension makes the bars in each group.
"""
function grouped_bar_xyz(xyzM :: AbstractArray{F1, 3}, xyzLabelV;
    yLabel :: AbstractString, legendPos = :best,
    groupIdx :: Integer = 2, subPlotIdx :: Integer = 3) where F1 <: Number

    nSub = size(xyzM, subPlotIdx);
    legendIdx = legend_index(groupIdx, subPlotIdx);
    subLabelV = xyzLabelV[subPlotIdx];
    pV = Vector{Any}(undef, nSub);
    for iSub = 1 : nSub
        dataM = selectdim(xyzM, subPlotIdx, iSub);
        if legendIdx > groupIdx
            dataM = dataM';
        end
        # Legend only for first panel
        if iSub == 1
            subLegPos = legendPos;
        else
            subLegPos = :none;
        end
        pV[iSub] = grouped_bar_xy(dataM, xyzLabelV[groupIdx];
            legendPos = subLegPos,  legendV = make_legend(xyzLabelV[legendIdx]),
            xLabel = subLabelV[iSub], yLabel = yLabel)
    end
    p = plot(pV..., link = :all);
    return p
end

function legend_index(groupIdx, subPlotIdx)
    @assert groupIdx != subPlotIdx
    availV = trues(3);
    availV[groupIdx] = false;
    availV[subPlotIdx] = false;
    legIdx = findfirst(availV);
    return legIdx
end

function group_indices(subPlotIdx)
    if subPlotIdx == 3
        groupIdx = 2;
        legendIdx = 1;
    else
        groupIdx = 3;
        if subPlotIdx == 2
            legendIdx = 1;
        else
            legendIdx = 2;
        end
    end
    @assert sort([subPlotIdx, groupIdx, legendIdx]) == [1,2,3]
    return groupIdx, legendIdx
end

# -----------------


text:
function getMul(a,n)
  ##a is the array, n is the number upto which we need the sum ##
  i,j=0,0
  li=[]
  for i in a
    for j in 1:n
      if j%i == 0
        push!(li,j)
      end end end
  return li
end

function getSum(li)
  sum=0
  for i in li
    sum+=i
  end
  return sum
end
function main(a,n)
  li=getMul(a,n)
  @show(li)
end
print("::the Project euler problem 1::")
l = main([3,5],1000)
getSum(l)

max_stars_repo_name: vbsteja/code

max_issues_repo_name: vbsteja/code

max_forks_repo_name: vbsteja/code


text:
using JuliaParser
using JuliaParser.Diagnostics: AbstractDiagnostic, Incomplete
using Base.Test

function do_diag_test(text)
    ts = Lexer.TokenStream{Lexer.SourceLocToken}(text)
    ts.filename = "test"
    try
        Parser.parse(ts)
        error("Should have failed")
    catch err
        !isa(err, AbstractDiagnostic) && rethrow(err)
        return err
    end
end

let diags = do_diag_test("""
      foo(a,b
      bar(
  """)
#  none:1:6 error: Expected ')' or ','
#  abc(d
#       ^
#  none:1:4 note: to match '(' here
#  abc(d
#     ^
    @test length(diags.elements) == 2
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("a && && b")
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("print x")
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("a ? b c")
#  none:1:7 error: colon expected in "?" expression
#  a ? b c
#        ^
#  none:1:3 note: "?" was here
#  a ? b c
#    ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("1:)")
#  none:1:3 error: missing last argument in range expression
#  1:)
#    ^
#  none:1:1 note: start of range was here
#  1:)
#  ^~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("1:\n")
#   none:1:3 error: line break in ":" expression
#   1:
#     ^
#   none:1:1 note: start of range was here
#   1:
#   ^~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("1:<Int")
#   none:1:2 error: Invalid "<" in range expression. Did you mean "<:"?
#   1:<Int
#    ^~
#    <:
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :fixit
end

let diags = do_diag_test("function f() 1")
#    none:1:15 error: incomplete: "function" requires end
#    function f() 1
#                  ^
#    none:1:1 note: "function" began here
#    function f() 1
#    ^~~~~~~~
    @test isa(diags, Incomplete)
    diags = diags.d
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("if\nend")
#    none:1:1 error: missing condition in "if"
#    if
#    ^~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("if true\nelseif\nend")
#   none:2:1 error: missing condition in "elseif"
#   elseif
#   ^~~~~~
#   none:1:1 note: "if" was here
#   if true
#   ^~")
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("if true\nelseif\nend")
#   none:2:1 error: missing condition in "elseif"
#   elseif
#   ^~~~~~
#   none:1:1 note: "if" was here
#   if true
#   ^~")
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("if true\nelse if false\nend")
#   none:2:1 error: use elseif instead of else if
#   else if false
#   ^~~~~~~
#   elseif
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :fixit
end

let diags = do_diag_test("if true 1\nfinally\nend")
#   none:2:1 error: Unexpected "finally" in if expression
#   finally
#   ^~~~~~~
#   none:1:1 note: previous "if" was here
#   if true 1
#   ^~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("let a b")
#   none:1:7 error: let variables should end in ";" or newline
#   let a b
#         ^
#   none:1:1 note: "let" was here
#   let a b
#   ^~~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("function foo{T} end")
# none:1:16 error: expected "(" in function definition
# function foo{T} end
#                ^
# none:1:1 note: "function" was here
# function foo{T} end
# ^~~~~~~~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("type type\nend")
# none:1:6 error: invalid type name "type"
# type type
#      ^~~~
# none:1:1 note: type was here
# type type
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("try 1\nelse\nend")
# none:2:1 error: Unexpected "else" in try expression
# else
# ^~~~
# none:1:1 note: "try" was here
# try 1
# ^~~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("const sin(2)")
# none:1:7 error: expected assignment after "const"
# const sin(2)
#       ^~~~~~
# none:1:1 note: "const" was here
# const sin(2)
# ^~~~~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("export true")
# none:1:8 error: invalid literal "true" in "export" statement
# export true
#        ^~~~
# none:1:1 note: "export" was here
# export true
# ^~~~~~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("ccall|")
# none:1:6 error: Expected '('
# ccall|
#      ^
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("do")
# none:1:1 error: invalid "do" syntax
# do
# ^~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("[a,a}")
# none:1:5 error: Expected "]", got "}"
# [a,a}
#     ^
# none:1:1 note: Expression began here
# [a,a}
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("[a,a b")
# none:1:6 error: missing separator in array expression
# [a,a b
#      ^
# none:1:1 note: Expression began here
# [a,a b
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("export ]")
# none:1:8 error: invalid syntax: "]"
# export ]
#        ^
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("/= = 1")
# none:1:1 error: invalid identifier name "/="
# /= = 1
# ^~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("@[1]")
# none:1:3 error: invalid macro use "@([1])"
# @[1]
#   ^~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("[x for (x ? true : false)]")
# none:1:9 error: invalid iteration spec
# [x for (x ? true : false)]
#         ^~~~~~~~~~~~~~~~
# none:1:4 note: for this `for`
# [x for (x ? true : false)]
#    ^~~
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("[ _ for i = 1:10}")
# none:1:17 error: expected ']' not "}"
# [ _ for i = 1:10}
#                 ^
# none:1:1 note: in comprehension that began here
# [ _ for i = 1:10}
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
    @test diags.elements[2].location.offset == 0
end

let diags = do_diag_test("[a b;\nc d}")
# none:2:5 error: unexpected "}"
#  c d}
#     ^
# none:1:1 note: In matrix expression that began here
# [a b;
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("[a b;\nc d for a = 1:10]")
# none:2:6 error: invalid comprehension syntax
#  c d for a = 1:10]
#      ^~~
# none:1:1 note: In matrix expression that began here
# [a b;
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("(a,b c)")
# none:1:5 error: missing separator in tuple
# (a,b c)
#     ^
# none:1:3 note: tuple began here
# (a,b c)
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

for expr in ["'","''","'\\","'abc'","`\\","\"","\"\"\"","\"\"\"\"","\"\"\"\"\""]
    let diags = do_diag_test(expr)
        isa(diags, Incomplete) && (diags = diags.d)
        @test diags.elements[1].severity == :error
    end
end

let diags = do_diag_test("=")
# none:1:1 error: unexpected "="
# =
# ^
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test(":(= a)")
# none:1:4 error: Expected ')'
# :(= a)
#    ^
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("(x]")
# none:1:3 error: expected ')' not ']'
# (x]
#   ^
# none:1:1 note: to match '(' here
# (x]
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("(a b)")
# none:1:3 error: missing spearator (did you mean ',' or ';')
# (a b)
#   ^
# none:1:1 note: in expression that began here
# (a b)
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("\"\$(a,)\"")
# none:1:5 error: Expected ')'
# "$(a,)"
#     ^
# none:1:2 note: In interpolation syntax starting here
# "$(a,)"
#  ^
# none:1:1 note: In string expression beginning here
# "$(a,)"
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
    @test diags.elements[3].severity == :note
end

let diags = do_diag_test("\"\$[1]\"")
# none:1:2 error: invalid interpolation syntax: "["
# "$[1]"
#  ^
# none:1:2 note: In interpolation syntax starting here
# "$[1]"
#  ^
# none:1:1 note: In string expression beginning here
# "$[1]"
# ^
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
    @test diags.elements[3].severity == :note
end

let diags = do_diag_test("\\xff")
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("0f1.2")
# none:1:1 error: invalid numeric constant "0f1."
# 0f1.2
# ^~~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("1._a")
# none:1:1 error: invalid use of '_' in numeric constant ""
# 1._a
# ^~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("0b12")
# none:1:1 error: invalid numeric constant "0b12"
# 0b12
# ^~~
    @test diags.elements[1].severity == :error
end

let diags = do_diag_test("#=")
# none:1:3 error: incomplete: unterminated multi-line comment #= ... =#
# #=
#   ^
# none:1:1 note: starting here
# #=
# ^
    @test isa(diags, Incomplete)
    diags = diags.d
    @test diags.elements[1].severity == :error
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("for i = 1:")
# none:1:11 error: missing last argument in range expression
# for i = 1:
#           ^
# none:1:9 note: start of range was here
# for i = 1:
#         ^~
    @test diags.elements[1].severity == :error
    @test diags.elements[1].location.offset == endof("for i = 1:")
    @test diags.elements[2].severity == :note
end

let diags = do_diag_test("ccall((:SetVarDeclInit, libcxxffi), Void, (Ptr{Void},) VD)")
# none:1:56 error: Expected ')' or ','
# ccall((:SetVarDeclInit, libcxxffi), Void, (Ptr{Void},) VD)
#                                                        ^
# none:1:6 note: to match '(' here
# ccall((:SetVarDeclInit, libcxxffi), Void, (Ptr{Void},) VD)
#      ^
    @test diags.elements[1].severity == :error
    @test diags.elements[1].location.offset == 55
    @test diags.elements[2].severity == :note
end

# fix issue #47

# test string with charwidth < 1
let code = "(foo̧ + b̂ar))",
    diag = do_diag_test(code)

    io = IOBuffer()
    JuliaParser.Diagnostics.display_diagnostic(io, code, diag, filename="test")
    res = String(take!(io))
    @test startswith(res, "test:1:12")
end

# test string with charwidth > 1
let code = "(foo̧ + 🔨))",
    diag = do_diag_test(code)
    io = IOBuffer()
    JuliaParser.Diagnostics.display_diagnostic(io, code, diag, filename="test")
    res = String(take!(io))
    @test startswith(res, "test:1:11")
end


text:
export buildModel, copyModel, cutModel, flip, rotate

"""
    buildModel(m::Array{Block, 3}, pos::Tuple{Real, Real, Real})

Build a model from a 3-D block array at `pos`.
"""
function buildModel(m::Array{Block,3}, pos::Tuple{Real, Real, Real})
    for X in 1:size(m, 1)
        for Y in 1:size(m, 2)
            for Z in 1:size(m, 3)
                setBlock(pos .+ (X, Y - 1, Z) , m[X, Y, Z])
            end
        end
    end
end

"""
    copyModel(p1::Tuple{Real, Real, Real}, p2::Tuple{Real, Real, Real})

Copy the Blocks between the diagonally opposite blocks and store in a 3-D Block Array.
"""
function copyModel(p1::Tuple{Real, Real, Real}, p2::Tuple{Real, Real, Real})
    p1 = round.(Int, p1)
    p2 = round.(Int, p2)
    dims = (abs.(p2 .- p1) .+ (1, 1, 1))
    m = Array{Block, 3}(undef, dims)
    for X in 1:dims[1]
        for Y in 1:dims[2]
            for Z in 1:dims[3]
                m[X, Y, Z] = getBlock(min.(p1, p2) .+ (X, Y, Z) .- (1, 1, 1))
            end
        end
    end
    return m
end

"""
    cutModel(p1::Tuple{Real, Real, Real}, p2::Tuple{Real, Real, Real}, block::Block)

Copy blocks between diagonally opposite blocks and fill the space with `block`.
"""
function cutModel(p1::Tuple{Real, Real, Real}, p2::Tuple{Real, Real, Real}, block::Block = Block(0, 0))
    m = copyModel(p1, p2)
    setBlocks(p1, p2, block)
    return m
end

"""
    flip(A::Array{Block, 3}, dim::Symbol)

Flip a 3-D Block Array in a particular dimension. Dimensions are `x`, `y` and `z`

# Example
```julia-repl
julia> A = reshape(Block.(collect(1:8)), (2,2,2))
2×2×2 Array{PiCraft.Block,3}:
[:, :, 1] =
 PiCraft.Block(1, 0)  PiCraft.Block(3, 0)
 PiCraft.Block(2, 0)  PiCraft.Block(4, 0)

[:, :, 2] =
 PiCraft.Block(5, 0)  PiCraft.Block(7, 0)
 PiCraft.Block(6, 0)  PiCraft.Block(8, 0)

julia> flip(A, :z)
2×2×2 Array{PiCraft.Block,3}:
[:, :, 1] =
 PiCraft.Block(5, 0)  PiCraft.Block(7, 0)
 PiCraft.Block(6, 0)  PiCraft.Block(8, 0)

[:, :, 2] =
 PiCraft.Block(1, 0)  PiCraft.Block(3, 0)
 PiCraft.Block(2, 0)  PiCraft.Block(4, 0)
```
"""
flip(A::Array{Block, 3}, dim::Symbol) = flip(A, Val{dim})

function flip(A, ::Type{Val{:x}})
    B = similar(A)
    for i = 1:size(A, 1)
        B[i, :, :] = A[size(A, 1) - i + 1, :, :]
    end
    return B
end

function flip(A, ::Type{Val{:y}})
    B = similar(A)
    for i = 1:size(A, 2)
        B[:, i, :] = A[:, size(A, 2) - i + 1, :]
    end
    return B
end

function flip(A, ::Type{Val{:z}})
    B = similar(A)
    for i = 1:size(A, 3)
        B[:, :, i] = A[:, :, size(A, 3) - i + 1]
    end
    return B
end

"""
    rotate(A::Array{Block, 3}; dir::Symbol = :l, dim::Symbol = :x)

Rotate a 3-D Block Array 90 degrees in a dimension and direction.
Dimensions are `x`, `y` and `z`. Directions are `l` & `r`.

# Example
```julia-repl
julia> A = reshape(Block.(collect(1:8)), (2,2,2))
2×2×2 Array{PiCraft.Block,3}:
[:, :, 1] =
 PiCraft.Block(1, 0)  PiCraft.Block(3, 0)
 PiCraft.Block(2, 0)  PiCraft.Block(4, 0)

[:, :, 2] =
 PiCraft.Block(5, 0)  PiCraft.Block(7, 0)
 PiCraft.Block(6, 0)  PiCraft.Block(8, 0)

julia> rotate(A, dir = :r, dim = :y)
 2×2×2 Array{PiCraft.Block,3}:
 [:, :, 1] =
  PiCraft.Block(3, 0)  PiCraft.Block(7, 0)
  PiCraft.Block(4, 0)  PiCraft.Block(8, 0)

 [:, :, 2] =
  PiCraft.Block(1, 0)  PiCraft.Block(5, 0)
  PiCraft.Block(2, 0)  PiCraft.Block(6, 0)
```
"""
rotate(A::Array{Block, 3}, dir::Symbol, dim::Symbol) = rotate(A, Val{dir}, Val{dim})
rotate(A::Array{Block, 3}; dir::Symbol = :l, dim::Symbol = :x) = rotate(A, Val{dir}, Val{dim})

function rotate(A, ::Type{Val{:l}}, ::Type{Val{:x}})
    B = similar(A, (size(A, 1), size(A, 3), size(A, 2)))
    for i in 1:size(A, 1)
        B[i,:,:] = rotl90(A[i,:,:])
    end
    return B
end

function rotate(A, ::Type{Val{:l}}, ::Type{Val{:y}})
    B = Array{Block, 3}(undef, (size(A, 3), size(A, 2), size(A, 1)))
    for i in 1:size(A, 2)
        B[:, i, :] = rotl90(A[:, i, :])
    end
    return B
end

function rotate(A, ::Type{Val{:l}}, ::Type{Val{:z}})
    B = Array{Block, 3}(undef, (size(A, 2), size(A, 1), size(A, 3)))
    for i in 1:size(A, 1)
        B[:, :, i] = rotl90(A[:, :, i])
    end
    return B
end

function rotate(A, ::Type{Val{:r}}, ::Type{Val{:y}})
    B = Array{Block, 3}(undef, (size(A, 1), size(A, 3), size(A, 2)))
    for i in 1:size(A, 1)
        B[i,:,:] = rotr90(A[i,:,:])
    end
    return B
end

function rotate(A, ::Type{Val{:r}}, ::Type{Val{:z}})
    B = Array{Block, 3}(undef, (size(A, 3), size(A, 2), size(A, 1)))
    for i in 1:size(A, 1)
        B[:, i, :] = rotr90(A[:, i, :])
    end
    return B
end

function rotate(A, ::Type{Val{:r}}, ::Type{Val{:x}})
    B = Array{Block, 3}(undef, (size(A, 2), size(A, 1), size(A, 3)))
    for i in 1:size(A, 1)
        B[:, :, i] = rotr90(A[:, :, i])
    end
    return B
end


text:
include(joinpath(pwd(), "examples/implicit_dynamics/solver/interior_point.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/solver/lu.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/solver/cones.jl"))

include(joinpath(pwd(), "examples/implicit_dynamics/utils.jl"))
include(joinpath(pwd(), "examples/implicit_dynamics/dynamics2.jl"))

max_issues_repo_name: thowell/DirectMotionPlanning


text:
# RGI-based runs for ESA

using BSON, VAWTools
region = [3, 4, 14, 11, 7][3]
runtyp = ["test", "testmid", "prodlow", "prod"][3]
runtyp == :test && println("\n\nTEST RUN !!!!!!!!\n\n")
parallel = true
all_glaciers = parallel
repeat_mode_1 = 1 # number of mode-1 repetitions to be able to calculate Rhat
repeat_mode_2 = 1 # number of mode-2 repetitions to be able to calculate Rhat

fit_sigma = true
use_synthetic_dhdt = true

today = Dates.today()
if runtyp in ["test", "testmid"]
    dir_output = "output/rasters-test/$(lowercase(Dates.monthabbr(today)))-$(VAWTools.int2str2(Dates.day(today)))"
else
    if fit_sigma
        dir_output = "output/rasters-fit_sigma/$(lowercase(Dates.monthabbr(today)))-$(VAWTools.int2str2(Dates.day(today)))"
    else
        dir_output = "output/rasters/$(lowercase(Dates.monthabbr(today)))-$(VAWTools.int2str2(Dates.day(today)))"
    end
end
if use_synthetic_dhdt
    dir_output *= "-without-dhdt/"
else
    dir_output *= "/"
end

!isdir(dir_output) && mkpath(dir_output)

region_str = VAWTools.int2str2(region)

@show region, runtyp, parallel, dir_output

include("rgi-setup.jl") # makes sure all is pre-compiles

tmp = BM.get_git_shas(error_on_dirty=false)
description = """
    Now looking into getting Rhat convergence stats after all.

    BITEModel commit $(tmp["BITEModel"][1])
              branch $(tmp["BITEModel"][2])
    """

if parallel && nprocs()<2
    # if Sys.CPU_CORES>=17

    #     addprocs(17) # /2 to get to physical cores
    # else
    #     addprocs(Sys.CPU_CORES÷2 + 1) # /2 to get to physical cores
    # end
    # addprocs(Sys.CPU_CORES÷5) # five regions
    addprocs(Sys.CPU_CORES)
end

@everywhere begin
    const debugging = true
    include("rgi-setup.jl")
    using NamedTuples
    const progress_bar = true
end

update_cache = true
glnrs1, glnrs2, glnrs_all = get_glaciers(region)

# To reduce to just a few glacier both of mode1 & 2
if !all_glaciers
    glnrs1 = length(glnrs1)>1 ? glnrs1[1:2] : glnrs1
    glnrs2 = glnrs2[1:5]
end

## MCMC

# parameters to play with
sigma_of_model = @NT(sigma_h_model = 10.0,
                     sigma_iv_model = 10.0)

# number of divisions in elevation band variables
fit_vars = @NT(n_1d_btilde=3,
               n_1d_fsl=3,
               n_1d_temp=1)

fit_target = [BM.FitTarget.h, BM.FitTarget.h_iv, BM.FitTarget.length][1]

@show region, debugging, runtyp, length(glnrs1)

# Full fitting Mode 1
#####################
start_tic = time()

# Method errors contain their arguments, which in our case can be huge.
# Remove them in case of error in pmap
@everywhere strip_MethodErrors(err) = err isa MethodError ? MethodError(0, err.f, err.world) : err
for rm1 in 1:repeat_mode_1
    tic = time()

    # runtyp=:prodlow
    # nr = 2819
    # out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(1), dir_output="", retsol=true, fit_sigma=fit_sigma)
    # error()

    # out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(1), dir_output="", retsol=true, fit_sigma=fit_sigma)
    global out1 = pmap(nr -> fit_it(nr, region, sigma_of_model, fit_vars,
                                    fit_target, runtyp, Val(1), dir_output=dir_output,
                                    fit_sigma=fit_sigma, store_thetas=true, use_synthetic_dhdt=use_synthetic_dhdt),
                       glnrs1, on_error=strip_MethodErrors, retry_delays = zeros(3))

    toc1 = time() - tic
    println("Iteration over mode-1 glaciers completed in $(toc1/60)min")

    ext = repeat_mode_1>1 ? "-$rm1" : ""
    println("Storing $dir_output/$(region_str)_out1$ext.bson ...")
    write(dir_output *"/$(region_str)_out$ext.txt", description)
    @time bson(dir_output *"/$(region_str)_out1$ext.bson", Dict(:out1 => deepcopy(out1),
                                                                :description => description,
                                                                :commit => BM.get_git_shas(error_on_dirty=false),
                                                                :glnrs1 => glnrs1,
                                                                :glnrs2 => glnrs2,
                                                                :runtyp => runtyp,
                                                                ))
    #out1 = BSON.load(dir_output *"/$(region_str)_out1.bson")[:out1];
    println("... done.")

    if length(out1)>0 && length([o[:thetas_expect] for o in out1 if o isa Dict])==0
        error("All Mode-1 runs errored")
    end
end

# Make the priors of the fitted parameters
##########################################
#
# I'm just doing the easy thing here and using the expectation and std of
# the thetas as component-wise prior uppdate.

function calculate_prior_update_reg14()
    # specify files
    fls = ["output/rasters-fit_sigma/production-with-dhdt/07_out1.bson",
           "output/rasters-fit_sigma/production-with-dhdt/11_out1.bson",
           "output/rasters-fit_sigma/production-without-dhdt/03_out1.bson",
           "output/rasters-fit_sigma/production-without-dhdt/04_out1.bson"]

    ## use a mean
    # exs = []
    # sts = []
    # for f in fls
    #     out1 = BSON.load(f)[:out1];
    #     expect, expect_median, expect_mean, stdev = calculate_prior_update(out1, true)
    #     push!(exs, expect_median)
    #     push!(sts, stdev)
    # end
    # return exs, sts, mean(hcat(exs...),2)[:], mean(hcat(sts...),2)[:]

    # merge all runs into one.  I.e. treat it as one big region.
    out1 = []
    for f in fls
        o = BSON.load(f)[:out1]
        # @show calculate_prior_update(o, true)[1]
        append!(out1, o)
    end
    expect, expect_median, expect_mean, stdev = calculate_prior_update(out1, true)
    return expect, stdev

end
if region==14
    expect, stdev = calculate_prior_update_reg14()
    expect_mean = nothing
    expect_median = expect
else
    expect, expect_median, expect_mean, stdev = calculate_prior_update(out1, true)
end

# Priors-only fitting Mode 2
############################
for rm1 in 1:repeat_mode_2
    tic2 = time()
    #out,sol = fit_it(nr, region, sigma_of_model, fit_vars, fit_target, runtyp, Val(2), expect, stdev, dir_output="", retsol=true, fit_sigma=false)
    global out2 = pmap(nr -> fit_it(nr, region, sigma_of_model, fit_vars,
                             fit_target, runtyp, Val(2),
                             expect, stdev, dir_output=dir_output,
                             fit_sigma=false, store_thetas=true, use_synthetic_dhdt=use_synthetic_dhdt),
                glnrs2, on_error=strip_MethodErrors, retry_delays = zeros(3))

    toc2 = time() - tic2

    println("Iteration over mode-2 glacier completed in $(toc2/60)min")
    ext = repeat_mode_2>1 ? "-$rm1" : ""
    println("Storing $dir_output/$(region_str)_out2$ext.bson ...")
    @time bson(dir_output *"/$(region_str)_out2$ext.bson", Dict(:out2 => deepcopy(out2),
                                                                :description => description,
                                                                :commit => BM.get_git_shas(error_on_dirty=false),
                                                                :expect => expect,
                                                                :expect_mean => expect_mean,
                                                                :expect_median => expect_median,
                                                                :stdev => stdev,
                                                                ))
    println("... done.")
end

toc = time() - start_tic
println("Total time $(toc/60/60) hours")

if length(out2)>0 && length([o[:thetas_expect] for o in out2 if o isa Dict])==0
    error("All Mode-2 runs errored")
end


text:
module DataGenerator

defaultβ(p) = collect(range(-1, stop=1, length=p))

function linregdata(n, p; β = defaultβ(p), σ = 1.0)
	x = randn(n, p)
	y = x * β + σ * randn(n)
	return x, y, β
end

function logregdata(n, p, sgn = true; β = defaultβ(p))
	x = randn(n, p)
	y = Float64[rand() < 1/(1 + exp(-η)) for η in x*β]
	if sgn
		y .= 2 .* y .- 1
	end
	return x, y, β
end

# function poissonregdata(n, p, β = defaultβ(p), V = eye(p), μ = zeros(p))
# 	x = rand(D.MvNormal(μ, V), n)'
# 	y = Float64[rand(D.Poisson(exp(η))) for η in x * β]
# 	return x, y, β
# end
end #module