text:
module CUSPARSE

using ..CuArrays
using ..CuArrays: active_context, unsafe_free!, @argout, @workspace

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

import CUDAnative

using CEnum

const SparseChar = Char

const libcusparse = Ref("libcusparse")

# core library
include("libcusparse_common.jl")
include("error.jl")
include("libcusparse.jl")

# low-level wrappers
include("array.jl")
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

const _handles = Dict{CuContext,cusparseHandle_t}()
const _handle = Ref{cusparseHandle_t}()

function handle()
    if _handle[] == C_NULL
        CUDAnative.maybe_initialize("CUSPARSE")
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cusparseCreate()
            atexit(()->CUDAdrv.isvalid(context) && cusparseDestroy(handle))
            handle
        end
    end

    return _handle[]
end

end


text:
using Test

@testset "soilstorage_run" begin
    for i in 1:1
        #parameters that don't change
        beta = rand(0.5:0.1:4)[1]
        Ce =  rand(0.4:0.01:0.8)[1]
        #Percolationcapacity = rand(0.1:0.1:4)[1]
        Ratio_Pref = rand(0.1:0.1:0.9)[1]
        Soilstoragecapacity = rand(50.:5:400.)[1]

        #parameters that change
        #Effective_Precipitation = rand(0.1:0.5:15)[1]
        #Interception_Evaporation = rand(0.1:0.1:4)[1]
        #Potential_Evaporation = rand(0.1:0.2:10)[1]
        #Storage = rand(1:5:Soilstoragecapacity)[1]

        # if no precipitation and storage zero, all outputs are zero
        Storage = 0.
        for i in 1:1000000
            Effective_Precipitation = rand(0.1:0.5:15)[1]
            Potential_Evaporation = rand(0.1:0.2:10)[1]
            Interception_Evaporation = rand(0:0.1:Potential_Evaporation * 0.5)[1]
            Overlandflow, Preferentialflow, Soil_Evaporation, Soilstorage = soilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Storage, beta, Ce, Ratio_Pref, Soilstoragecapacity)
            Storage = Soilstorage
            @test Storage <= Soilstoragecapacity
        end
    end
end

@testset "soilstorage_rip_run" begin
    for i in 1:1
        #parameters that don't change
        beta = rand(0.5:0.1:4)[1]
        Ce =  rand(0.4:0.01:0.8)[1]
        Drainagecapacity = rand(0.1:0.1:4)[1]
        #Ratio_Pref = rand(0.1:0.1:0.9)[1]
        Soilstoragecapacity = rand(50.:5:400.)[1]

        #parameters that change

        # if no precipitation and storage zero, all outputs are zero
        Storage = 0.
        for i in 1:1000000
            Effective_Precipitation = rand(0.1:0.5:15)[1]
            Potential_Evaporation = rand(0.1:0.2:5)[1]
            Interception_Evaporation = rand(0.:0.1:Potential_Evaporation * 0.5)[1]
            Riparian_Discharge = rand(0.1:0.2:10)[1]
            Overlandflow, Soil_Evaporation, Soilstorage = ripariansoilstorage(Effective_Precipitation, Interception_Evaporation, Potential_Evaporation, Riparian_Discharge, Storage, beta, Ce, Drainagecapacity, Soilstoragecapacity)
            Storage = Soilstorage
            @test Storage <= Soilstoragecapacity
        end
    end
end

ext: jl

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
function getPawnAttack(attack::UInt64, b::UInt64, white::Bool)
    white ? lx_clear = CLEAR_FILE_A : lx_clear = CLEAR_FILE_H
    white ? lx_shift = 9 : lx_shift = -9
    attack |= (b & lx_clear) << lx_shift
    white ? rx_shift = 7 : rx_shift = -7
    white ? rx_clear = CLEAR_FILE_H : rx_clear = CLEAR_FILE_A
    return attack |= (b & rx_clear) << rx_shift
end

function getPawnMoves!(moves::Moves, bitboard::UInt64, taken::UInt64, friends::UInt64,
                       enemy::ChessSet, white::Bool, enpassant::UInt64)
    white ? shift = 8 : shift = -8
    white ? home = MASK_RANKS[2] : home = MASK_RANKS[7]
    white ? prom = MASK_RANKS[8] : prom = MASK_RANKS[1]

    one_step = (bitboard << shift) & ~taken
    if one_step != EMPTY
        b0, d = bitshift(one_step)
        for j = 0:d
            if b0<<j & one_step != EMPTY
                target = b0<<j
                if target & prom != EMPTY
                    src = target>>shift
                    for p in [:knight, :bishop, :rook, :queen]
                        push!(moves, Move(:pawn, src, target, NONE, EMPTY, p, NOCASTLING))
                    end
                else
                    push!(moves, Move(:pawn, target>>shift, target, NONE, EMPTY, :none, NOCASTLING))
                end
            end
        end
    end

    two_steps = ((((bitboard & home) << shift) & ~taken) << shift) & ~taken
    if two_steps != EMPTY
        b0, d = bitshift(two_steps)
        for j = 0:d
            if b0<<j & two_steps != EMPTY
                target = b0<<j
                push!(moves, Move(:pawn, target>>(shift*2), target, NONE, target>>shift, :none, NOCASTLING))
            end
        end
    end

    white ? lx_clear = CLEAR_FILE_A : lx_clear = CLEAR_FILE_H
    white ? lx_shift = 9 : lx_shift = -9
    x_lx = (((bitboard & lx_clear) << lx_shift) & taken) & ~friends
    if x_lx != EMPTY
        b0, d = bitshift(x_lx)
        for j = 0:d
            if b0<<j & x_lx != EMPTY
                target = b0<<j
                if target & enemy.friends != EMPTY
                    take = Piece(getTypeAt(enemy, target), target)
                    if target & prom != EMPTY
                        src = target>>lx_shift
                        for p in  [:knight, :bishop, :rook, :queen]
                            push!(moves, Move(:pawn, src, target, take, EMPTY, p, NOCASTLING))
                        end
                    else
                        push!(moves, Move(:pawn, target>>lx_shift, target, take, EMPTY, :none, NOCASTLING))
                    end
                end
            end
        end
    end

    en_lx = (((bitboard & lx_clear) << lx_shift) & ~taken) & enpassant
    if en_lx != EMPTY
        b0, d = bitshift(en_lx)
        take = Piece(:pawn, enpassant<<-shift)
        push!(moves, Move(:pawn, b0>>lx_shift, enpassant, take, EMPTY, :none, NOCASTLING))
    end

    white ? rx_shift = 7 : rx_shift = -7
    white ? rx_clear = CLEAR_FILE_H : rx_clear = CLEAR_FILE_A
    x_rx = (((bitboard & rx_clear) << rx_shift) & taken) & ~friends
    if x_rx != EMPTY
        b0, d = bitshift(x_rx)
        for j = 0:d
            if b0<<j & x_rx != EMPTY
                target = b0<<j
                if target & enemy.friends != EMPTY
                    take = Piece(getTypeAt(enemy, target), target)
                    if target & prom != EMPTY
                        src = target>>rx_shift
                        for  p in  [:knight, :bishop, :rook, :queen]
                            push!(moves, Move(:pawn, src, target, take, EMPTY, p, NOCASTLING))
                        end
                    else
                        push!(moves, Move(:pawn, target>>rx_shift, target, take, EMPTY, :none, NOCASTLING))
                    end
                end
            end
        end
    end

    en_rx = (((bitboard & rx_clear) << rx_shift) & ~taken) & enpassant
    if en_rx != EMPTY
        b0, d = bitshift(en_rx)
        take = Piece(:pawn, enpassant<<-shift)
        push!(moves, Move(:pawn, b0>>rx_shift, enpassant, take, EMPTY, :none, NOCASTLING))
    end
end

function pawnMovesGen(white::Bool)
    pawn_moves = Dict{UInt64,UInt64}()
    pawn_attacks = Dict{UInt64,UInt64}()

    white ? home = MASK_RANKS[2] : home = MASK_RANKS[7]
    white ? lx_clear = CLEAR_FILE_A : lx_clear = CLEAR_FILE_H
    white ? rx_clear = CLEAR_FILE_H : rx_clear = CLEAR_FILE_A

    for square in values(PGN2UINT)
        targets = EMPTY
        attacks = EMPTY
        white ? targets |= square << 8 : targets |= square >> 8
        if square & lx_clear != EMPTY
            attacks |= square << (white ? 9 : -9)
        end
        if square & rx_clear != EMPTY
            attacks |= square << (white ? 7 : -7)
        end
        if square & home != EMPTY
            white ? targets |= square << 16 : targets |= square >> 16
        end
        push!(pawn_moves, square=>targets)
        push!(pawn_attacks, square=>attacks)
    end

    return pawn_moves, pawn_attacks
end
const PAWN_WHITE, PAWN_X_WHITE = pawnMovesGen(true)
const PAWN_BLACK, PAWN_X_BLACK = pawnMovesGen(false)


text:
# An example of graph classification

using Flux
using Flux:onecold, onehotbatch
using Flux.Losses: logitbinarycrossentropy
using Flux.Data: DataLoader
using GraphNeuralNetworks
using MLDatasets: TUDataset
using Statistics, Random
using CUDA
CUDA.allowscalar(false)

function eval_loss_accuracy(model, data_loader, device)
    loss = 0.
    acc = 0.
    ntot = 0
    for g in data_loader
        g = g |> device
        n = g.num_graphs
        y = g.gdata.y
        ŷ = model(g, g.ndata.x) |> vec
        loss += logitbinarycrossentropy(ŷ, y) * n
        acc += mean((ŷ .> 0) .== y) * n
        ntot += n
    end
    return (loss = round(loss/ntot, digits=4), acc = round(acc*100/ntot, digits=2))
end

function getdataset()
    data = TUDataset("MUTAG")

    x = Array{Float32}(onehotbatch(data.node_labels, 0:6))
    y = (1 .+ Array{Float32}(data.graph_labels)) ./ 2
    @assert all(∈([0,1]), y) # binary classification
    # The dataset also has edge features but we won't be using them
    e = Array{Float32}(onehotbatch(data.edge_labels, sort(unique(data.edge_labels))))

    return GNNGraph(data.source, data.target,
                num_nodes=data.num_nodes,
                graph_indicator=data.graph_indicator,
                ndata=(; x), edata=(; e), gdata=(; y))
end

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 1f-3             # learning rate
    batchsize = 64      # batch size (number of graphs in each batch)
    epochs = 200         # number of epochs
    seed = 17             # set seed > 0 for reproducibility
    usecuda = true      # if true use cuda (if available)
    nhidden = 128        # dimension of hidden features
    infotime = 10 	     # report every `infotime` epochs
end

function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)

    if args.usecuda && CUDA.functional()
        device = gpu
        args.seed > 0 && CUDA.seed!(args.seed)
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    # LOAD DATA

    NUM_TRAIN = 150

    gfull = getdataset()

    @info gfull

    perm = randperm(gfull.num_graphs)
    gtrain, _ = getgraph(gfull, perm[1:NUM_TRAIN])
    gtest, _ = getgraph(gfull, perm[NUM_TRAIN+1:end])
    train_loader = DataLoader(gtrain, batchsize=args.batchsize, shuffle=true)
    test_loader = DataLoader(gtest, batchsize=args.batchsize, shuffle=false)

    # DEFINE MODEL

    nin = size(gtrain.ndata.x, 1)
    nhidden = args.nhidden

    model = GNNChain(GraphConv(nin => nhidden, relu),
                     Dropout(0.5),
                     GraphConv(nhidden => nhidden, relu),
                     GlobalPool(mean),
                     Dense(nhidden, 1))  |> device

    ps = Flux.params(model)
    opt = ADAM(args.η)

    # LOGGING FUNCTION

    function report(epoch)
        train = eval_loss_accuracy(model, train_loader, device)
        test = eval_loss_accuracy(model, test_loader, device)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
    end

    # TRAIN

    report(0)
    for epoch in 1:args.epochs
        for g in train_loader
            g = g |> device
            gs = Flux.gradient(ps) do
                ŷ = model(g, g.ndata.x) |> vec
                logitbinarycrossentropy(ŷ, g.gdata.y)
            end
            Flux.Optimise.update!(opt, ps, gs)
        end

        epoch % args.infotime == 0 && report(epoch)
    end
end

train()


text:
@unix_only begin
    const path_separator    = "/"
    const path_separator_re = r"/+"
    const path_absolute_re  = r"^/"
    const path_directory_re = r"(?:^|/)\.{0,2}$"
    const path_dir_splitter = r"^(.*?)(/+)([^/]*)$"
    const path_ext_splitter = r"^((?:.*/)?(?:\.|[^/\.])[^/]*?)(\.[^/\.]*|)$"

    splitdrive(path::String) = ("",path)
    homedir() = ENV["HOME"]
end
@windows_only begin
    const path_separator    = "\\"
    const path_separator_re = r"[/\\]+"
    const path_absolute_re  = r"^(?:\w+:)?[/\\]"
    const path_directory_re = r"(?:^|[/\\])\.{0,2}$"
    const path_dir_splitter = r"^(.*?)([/\\]+)([^/\\]*)$"
    const path_ext_splitter = r"^((?:.*[/\\])?(?:\.|[^/\\\.])[^/\\]*?)(\.[^/\\\.]*|)$"

    function splitdrive(path::String)
        m = match(r"^(\w+:|\\\\\w+\\\w+|\\\\\?\\UNC\\\w+\\\w+|\\\\\?\\\w+:|)(.*)$", path)
        bytestring(m.captures[1]), bytestring(m.captures[2])
    end
    homedir() = get(ENV,"HOME",string(ENV["HOMEDRIVE"],ENV["HOMEPATH"]))
end

isabspath(path::String) = ismatch(path_absolute_re, path)
isdirpath(path::String) = ismatch(path_directory_re, splitdrive(path)[2])

function splitdir(path::ByteString)
    a, b = splitdrive(path)
    m = match(path_dir_splitter,b)
    m == nothing && return (a,b)
    a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
    a, bytestring(m.captures[3])
end
splitdir(path::String) = splitdir(bytestring(path))

 dirname(path::String) = splitdir(path)[1]
basename(path::String) = splitdir(path)[2]

function splitext(path::String)
    a, b = splitdrive(path)
    m = match(path_ext_splitter, b)
    m == nothing && return (path,"")
    a*m.captures[1], bytestring(m.captures[2])
end

function pathsep(paths::String...)
    for path in paths
        m = match(path_separator_re, path)
        m != nothing && return m.match[1]
    end
    return path_separator
end

joinpath(a::String) = a
joinpath(a::String, b::String, c::String...) = joinpath(joinpath(a,b), c...)

function joinpath(a::String, b::String)
    isabspath(b) && return b
    A, a = splitdrive(a)
    B, b = splitdrive(b)
    !isempty(B) && A != B && error("drive mismatch: $A$a $B$b")
    C = isempty(B) ? A : B
    isempty(a)                             ? string(C,b) :
    ismatch(path_separator_re, a[end:end]) ? string(C,a,b) :
                                             string(C,a,pathsep(a,b),b)
end

function normpath(path::String)
    isabs = isabspath(path)
    isdir = isdirpath(path)
    drive, path = splitdrive(path)
    parts = split(path, path_separator_re)
    filter!(x->!isempty(x) && x!=".", parts)
    while true
        clean = true
        for j = 1:length(parts)-1
            if parts[j] != ".." && parts[j+1] == ".."
                splice!(parts, j:j+1)
                clean = false
                break
            end
        end
        clean && break
    end
    if isabs
        while !isempty(parts) && parts[1] == ".."
            shift!(parts)
        end
    elseif isempty(parts)
        push!(parts, ".")
    end
    path = join(parts, path_separator)
    if isabs
        path = path_separator*path
    end
    if isdir && !isdirpath(path)
        path *= path_separator
    end
    string(drive,path)
end
normpath(a::String, b::String...) = normpath(joinpath(a,b...))

abspath(a::String) = normpath(isabspath(a) ? a : joinpath(pwd(),a))
abspath(a::String, b::String...) = abspath(joinpath(a,b...))

@windows_only function realpath(path::String)
    buflength = length(path)+1
    buf = zeros(Uint8,buflength)
    p = ccall((:GetFullPathNameA, "Kernel32"), stdcall,
        Uint32, (Ptr{Uint8}, Uint32, Ptr{Uint8}, Ptr{Void}),
        path, buflength, buf, C_NULL)
    if p > buflength
        buf = zeros(Uint8,p)
        p = ccall((:GetFullPathNameA, "Kernel32"), stdcall,
            Uint32, (Ptr{Uint8}, Uint32, Ptr{Uint8}, Ptr{Void}),
            path, p, buf, C_NULL)
    end
    systemerror(:realpath, p == 0)
    return bytestring(buf)[1:end-1]
end

@unix_only function realpath(path::String)
    p = ccall(:realpath, Ptr{Uint8}, (Ptr{Uint8}, Ptr{Uint8}), path, C_NULL)
    systemerror(:realpath, p == C_NULL)
    s = bytestring(p)
    c_free(p)
    return s
end

@windows_only expanduser(path::String) = path # on windows, ~ means "temporary file"
@unix_only function expanduser(path::String)
    i = start(path)
    c, i = next(path,i)
    if c != '~' return path end
    if done(path,i) return homedir() end
    c, j = next(path,i)
    if c == '/' return homedir()*path[i:end] end
    error("~user tilde expansion not yet implemented")
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
using Documenter, Mathieu, PyPlot

DocMeta.setdocmeta!(Mathieu, :DocTestSetup, :(using Mathieu); recursive=true)

makedocs(
    sitename = "Mathieu.jl",
    authors = "Jérémy Béjanin",
    modules = [Mathieu],
    format = Documenter.HTML(prettyurls = false),
    linkcheck = true,
    clean = false,
    pages = [
        "Home" => "index.md",
        "Examples" => [
            "examples/mathieu_funs.md",
            "examples/mathieu_surf.md",
        ],
        "API" => [
            "api/values.md",
            "api/functions.md"
        ],
    ]
)

deploydocs(
    repo = "github.com/jebej/Mathieu.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
type BackTraceElement
    func::Symbol
    file::Symbol
    line::Int
end
BackTraceElement() = BackTraceElement(symbol(), symbol(), symbol(), -1)
typealias BACKTRACE Nullable{BackTraceElement}

function getbacktrace()
    bt = backtrace()
    btout = BackTraceElement[]
    for b in bt
        code = ccall(:jl_lookup_code_address, Any, (Ptr{Void}, Cint), b, true)
        if !code[6]
            push!(btout, BackTraceElement(code[1],code[2],code[3]))
        end
    end
    return btout
end

function moduledir(m::Module)
    mname = string(m)
    mpath = Base.find_in_path(mname)

    # could not find by name, try `eval` function
    if mpath === nothing
        try
            m = eval(:(first(methods($m.eval))))
            mpath = string(m.func.code.file)
        catch
        end
    end

    # no luck - return empty
    mpath === nothing && return ""

    # path does not have module name - no package, remove only file from path
    # otherwise go to th root folder of the package
    rmpaths = contains(mpath, mname) ? 2 : 1

    return join(split(mpath, Base.path_separator)[1:end-rmpaths], Base.path_separator)
end

function isfqmn(fqmn::AbstractString, m::Module=Main)
    submdls = split(fqmn, '.')
    for i in 1:length(submdls)
        ms = symbol(submdls[i])
        !isdefined(m, ms) && return false
        m = getfield(m, ms)
    end
    return true
end

function getmodule(fqmn::AbstractString, m::Module=Main)
    submdls = split(fqmn, '.')
    for i in 1:length(submdls)
        ms = symbol(submdls[i])
        !isdefined(m, ms) && return nothing
        m = getfield(m, ms)
    end
    return m
end

const LOG4JL_CTX_PREFIX = "__Log4jl_"
mlcvar(fqmn::AbstractString) = symbol("$(LOG4JL_CTX_PREFIX)LC$(string(hash(fqmn)))__")

function submodules(m::Module, mdls::Vector{Symbol}= Symbol[])
    for vs in names(m, true)
        v = getfield(m, vs)
        fqmn = symbol(v)
        if isa(v, Module) && vs ∉ [:Main, :Core, :Base] && fqmn ∉ mdls
            push!(mdls, fqmn)
            submodules(v, mdls)
        end
    end
    return mdls
end

function evaltype(stype::AbstractString, smodule::AbstractString="")
    stype = isempty(smodule) ? stype : "$(smodule).$(stype)"
    try
        eval(parse(stype))
    catch
        Void
    end
end

function getlevel(strlevel::AbstractString)
    lt = evaltype(strlevel, "Level")
    return LEVEL(lt === Void ? nothing : lt)
end

function configlevel(conf::Dict)
    return haskey(conf, "level") ? getlevel((conf["level"] |> uppercase)) : LEVEL()
end


text:
g = IntAdjacencyList(4,3)
add_edge!(g, 1, 2, 1)
add_edge!(g, 2, 3, 1)
add_edge!(g, 1, 4, 2)

@test num_edges(g) == 3

@test forward_neighbors(g, 1, 1) == [(2,1), (1,2)]


text:
#=
    Examples file
=#

using Nemo

include("../src/ExactODEReduction.jl")
using .ExactODEReduction: from_dense, invariant_subspace, fullspace

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

gen_matrices = [
    from_dense([
        1  0 0 0;
        0 -1 0 0;
        0  0 0 0;
        0  0 0 0
    ], QQ),
    from_dense([
        0 -1 0 0;
        1  0 0 0;
        0  0 0 0;
        0  0 0 0
    ], QQ),
    from_dense([
        0 0 0  0;
        0 0 0  0;
        0 0 1  0;
        0 0 0 -1
    ], QQ),
    from_dense([
        0 0  0 0;
        0 0  0 0;
        0 0  0 1;
        0 0 -1 0
    ], QQ)
]

invariants = invariant_subspace(gen_matrices)

e1, e2, e3, e4 = fullspace(4, QQ)

# may fail due to randomness of the used algorithm
@assert invariants in [ [e1, e2], [e3, e4] ]


text:
"""
    SumGreaterThan(x<:AbstractIntVar, v::Int)

Summing constraint, states that `x[1] + x[2] + ... + x[length(x)] >= lower`
"""
struct SumGreaterThan <: Constraint
    x                   ::Array{<:AbstractIntVar}
    lower               ::Int
    active              ::StateObject{Bool}
    numberOfFreeVars    ::StateObject{Int}
    sumOfFixedVars      ::StateObject{Int}
    freeIds             ::Array{Int}
    function SumGreaterThan(x::Array{<:AbstractIntVar}, lower, trailer)
        @assert !isempty(x)

        freeIds = zeros(length(x))
        for i in 1:length(x)
            freeIds[i] = i
        end

        constraint = new(x, lower, StateObject{Bool}(true, trailer), StateObject{Int}(length(x), trailer), StateObject{Int}(0, trailer), freeIds)
        for xi in x
            addOnDomainChange!(xi, constraint)
        end
        return constraint
    end
end

"""
    propagate!(constraint::SumGreaterThan, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`SumGreaterThan` propagation function. The pruning is quite superficial.
"""
function propagate!(constraint::SumGreaterThan, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    # Stop propagation if constraint not active
    if !constraint.active.value
        return false
    end

    # Computing maxSum, minSum, and refreshing other variables
    newNumberOfFreeVars = constraint.numberOfFreeVars.value
    sumOfMax, sumOfMin = constraint.sumOfFixedVars.value, constraint.sumOfFixedVars.value
    for i in newNumberOfFreeVars:-1:1
        currentId = constraint.freeIds[i]
        sumOfMax += maximum(constraint.x[currentId].domain)

        if isbound(constraint.x[currentId])
            setValue!(constraint.sumOfFixedVars, constraint.sumOfFixedVars.value + assignedValue(constraint.x[currentId]))
            constraint.freeIds[i] = constraint.freeIds[newNumberOfFreeVars]
            constraint.freeIds[newNumberOfFreeVars] = currentId
            newNumberOfFreeVars -= 1
        end
    end
    setValue!(constraint.numberOfFreeVars, newNumberOfFreeVars)

    # Checking feasibility
    if sumOfMax < constraint.lower
        return false
    end

    ### Filtering ###
    # Here we must have: x_i = - sum(x_j for j != i)
    # But we know that: - sum(x_j for j != i) >= - sum(max(x_j) for j != i)
    # And: - sum(max(x_j) for j != i) = max(x_i) - sumOfMax
    # Hence x_i >= lower + max(x_i) - sumOfMax
    # Hence we remove everything below that last value from the domain of x_i
    # The reasoning is equivalent for the minimization
    for i in newNumberOfFreeVars:-1:1
        currentId = constraint.freeIds[i]
        currentMax = maximum(constraint.x[currentId].domain)
        pruned = vcat(removeBelow!(constraint.x[currentId].domain, constraint.lower + currentMax - sumOfMax))
        if !isempty(pruned)
            addToPrunedDomains!(prunedDomains, constraint.x[currentId], pruned)
            triggerDomainChange!(toPropagate, constraint.x[currentId])

        end
    end

    if newNumberOfFreeVars == 0
        setValue!(constraint.active, false)
    end

    return true
end

function Base.show(io::IO, ::MIME"text/plain", con::SumGreaterThan)
    ids = [var.id for var in con.x]
    print(io, typeof(con), ": ", join(ids, " + "), " ≥ ", con.lower, ", active = ", con.active)
    for var in con.x
        print(io, "\n   ", var)
    end
end

function Base.show(io::IO, con::SumGreaterThan)
    ids = [var.id for var in con.x]
    print(io, typeof(con), ": ", join(ids, " + "), " ≥ ", con.lower)
end


text:
module CMA_ES

using LinearAlgebra

export CMAES_evolution

mutable struct CMAES
    n::Int64
    λ::Int64
    μ::Int64

    x̄::Array{Float64, 1}

    σ::Float64
    ω::Array{Float64, 1}
    μω::Float64
    cσ::Float64
    cc::Float64
    c1::Float64
    cμ::Float64

    pσ::Array{Float64, 1}
    pc::Array{Float64, 1}

    C::Array{Float64, 2}

    CMAES(dim) = begin
        n = dim
        λ = 4 + trunc(Int, 3 * log(n))
        μ = trunc(Int, λ / 2)

        x̄ = rand(n)

        σ = 1
        ω = (-log.(1:μ) .+ log.(μ+0.5)) ./ (sum(-log.(1:μ) .+ log.(μ+0.5)))
        μω = 1 / sum(ω.^2)
        cσ = (μω + 2) / (n+μω+5)
        cc = 4/(n+4)
        c1 = 2/((n+1.3)^2+μω)
        cμ = min(1-c1, 2*(μω-2+1/μω)/((n+2)^2+μω))
        pσ = zeros(n)
        pc = zeros(n)
        C = Array{Float64, 2}(I, (n,n))
        new(n, λ, μ, x̄, σ, ω, μω, cσ, cc, c1, cμ, pσ, pc, C)
    end
end

mutable struct Population
    z::Array{Float64, 2}
    d::Array{Float64, 2}
    x::Array{Float64, 2}
    fit::Array{Float64, 1}
    Population(num, dim) = begin
        z = zeros(num, dim)
        d = zeros(num, dim)
        x = zeros(num, dim)
        fit = zeros(num)
        new(z, d, x, fit)
    end
end

function sample(es::CMAES)
    D, B = eigen(es.C)
    D, B = real(D), real(B)
    M = B * sqrt.(D)
    np = Population(es.λ, es.n)
    for i = 1:es.λ
        np.z[i, :] = randn(es.n)
        np.d[i, :] = M .* np.z[i, :]
        np.x[i, :] = es.x̄ + es.σ * np.d[i, :]
    end
    np
end

function selection(es::CMAES, np::Population, f)
    f(np)
    argx = sortperm(np.fit)
    argx[1:es.μ]
end

function update(es::CMAES, np::Population, argx)
    # update mean value
    es.x̄ += vec(es.σ .* sum(es.ω .* np.d[argx, :], dims=1))

    # update evolution path
    zz = vec(sum(es.ω .* np.z[argx, :], dims=1))
    c = sqrt(es.cσ * (2 - es.cσ) * es.μω)
    es.pσ -= es.cσ .* es.pσ
    es.pσ += c * zz

    dd = vec(sum(es.ω .* np.d[argx, :], dims=1))
    c = sqrt(es.cc * (2 - es.cc) * es.μω)
    es.pc -= es.cc * es.pc
    es.pc += c * dd

    es.C = Array{Float64, 2}(I, (es.n,es.n))
    # update Covariance matrix
    es.C = (1 - es.c1 - es.cμ) * es.C + es.c1 * es.pc * es.pc'
    part3 = zeros(es.n, es.n)
    for i in 1:es.μ
        es.C += es.cμ * es.ω[i] * np.d[argx[i], :] * np.d[argx[i], :]'
    end

    # update step size
    es.σ *= exp(es.cσ / 2 * (sum(es.pσ .^ 2) / es.n - 1))
end

function CMAES_evolution(f, dim::Int64, gen_max::Int64)
    es = CMAES(dim)
    opti_x = zeros(es.n)
    opti_f = 1e10
    for i in 1:gen_max
        np = sample(es)
        argx = selection(es, np, f)

        if opti_f > np.fit[argx[1]]
            opti_x, opti_f = np.x[argx[1]], np.fit[argx[1]]
        end
        println("gen: ", i, " fitness: ", opti_f)

        update(es, np, argx)
    end
    opti_x, opti_f
end

end  # modue EAFO


text:
module TestTrial

import Wordle
import Test: @test, @test_throws

trial = Wordle.Trial(Wordle.Candidate("alone"), Wordle.Result("⬜🟨🟩⬜🟨"))

@test trial.candidate.chars === ('a', 'l', 'o', 'n', 'e')

@test trial.result.codes === (
    Wordle.PRESENT_NOWHERE,
    Wordle.PRESENT_ELSEWHERE,
    Wordle.PRESENT_HERE,
    Wordle.PRESENT_NOWHERE,
    Wordle.PRESENT_ELSEWHERE,
)

end # module


text:
function exponential_1d(; invariant=Universe(1))
    s = @system(x' = -x, x ∈ invariant)
    x0 = Interval(0.4, 0.5)
    prob = InitialValueProblem(s, x0)
    tspan = (0.0, 1.0)
    return prob, tspan
end


text:
"""

The weighted rich club coefficient, Rw, at level k is the fraction of
edge weights that connect nodes of degree k or higher out of the
maximum edge weights that such nodes might share.

ndependent of each other, i.e. each hub has a specialised function, or they perform
some collaborative function. A graph measure to characterise the relation between the
hubs of a net- work is the rich-club phenomenon. The k-density φ(k), is defined
as the internal density of links between the nodes with degree larger than k′:

```math
\\psi(k) = \\frac{E_{>k}}{N_{>k}(N_{>k} - 1)}
```
where ``k`` is a threshold, ``E_{k}`` is the number of edges between nodes with
a degree greater than ``k`` and ``N_{k}`` is the number of nodes with a degree
greater than ``k``.

The weighted corrollary is:

```math
\\psi^{w}(k) = \\frac{W_{>k}}{\\sum_{i=1}^{E_{>k}} w_{i}} w_{i}^{ranked}}
```
Where  ``w_{i}^{ranked}`` is the ith index of the ranked graph weights and
``W_{>k}`` is the sum of weights from those nodes with a degree greater than
``k``.
CIJ = [0 1 1 0 0 0
        1 0 1 1 1 0
        1 1 0 0 0 0
        0 1 0 0 0 0
        0 1 0 0 0 1
        0 0 0 0 1 0]
d = [2
 4
 2
 1
 2
 1]
# Examples
```jldoctest
g = Graph([0 1 1 0 0 0
           1 0 1 1 1 0
           1 1 0 0 0 0
           0 1 0 0 0 0
           0 1 0 0 0 1
           0 0 0 0 1 0])
add_edge!(g, 1, 2)
add_edge!(g, 1, 3)
add_edge!(g, 2, 3)
add_edge!(g, 2, 4)
add_edge!(g, 2, 5)
add_edge!(g, 5, 6)

```

# Citations

@cite(vandenheuvelRichClubOrganizationHuman2011)

Zamora-López G, Zhou C, Kurths J. Cortical hubs form a module for multisensory
integration on top of the hierarchy of cortical networks. Front Neuroinform. 2010
Mar 19;4:1. doi: 10.3389/neuro.11.001.2010. eCollection 2010. PubMed [citation]
PMID: 20428515, PMCID: PMC2859882

"""

# rich club coefficient curve
function rich_club_coefficient(
    g::AbstractGraph,
    krange::AbstractRange{Int},
    dstmx::AbstractMatrix
   )

    E_r = sort(dstmx[:])
    dgs = degree(g)
    [__weighted_rich_club_coefficient(k, dstmx, E_r, dgs) for k in krange]
end

# rich club coefficient at single threshold
function rich_club_coefficient(g::AbstractGraph, k::Int, dstmx::AbstractMatrix)
    __weighted_rich_club_coefficient(k, dstmx, sort(dstmx[:]), degree(g))
end

function __weighted_rich_club_coefficient(
    k::Int,
    dstmx::AbstractMatrix{T},
    E_r::AbstractVector{T},  # ranked edges
    dgs::AbstractVector{<:Integer}  # degrees
   )
    __weighted_rich_club_coefficient(dstmx, dgs .>= k, E_r)
end

function __weighted_rich_club_coefficient(
    dstmx::AbstractMatrix{T},
    E_k::AbstractVector{Bool},
    E_r::AbstractVector{T}) where T
    E_n = 0
    W_k = zero(T)
    for (x,ek_x) in enumerate(E_k)
        for (y,ek_y) in enumerate(E_k)
            if ek_x && ek_y
                W_k = dstmx[x,y]
                E_n += 1
            end
        end
    end
    return W_k / sum(E_r[1:E_n])
end

function rich_club_coefficient(g::AbstractGraph, krange::AbstractRange{T}) where {T<:Integer}
    dstmx = weights(g)
    dgs = degree(g)
    [__rich_club_coefficient(k, dstmx, dgs) for k in krange]
end

function rich_club_coefficient(g::AbstractGraph, k::Int)
    d = outdegree(g)
    d = d[d .> k]
    Nk = sum(d)
    return 2* length(d) / (Nk * (Nk - 1))
end

#    __rich_club_coefficient(k, weights(g), degree(g))
end

function rich_club_coefficient(k::Int, dstmx::AbstractMatrix{<:Integer}, dgs::AbstractVector{<:Integer})
    __rich_club_coefficient(dstmx, dgs .>= k)
end

# unweighted version
"""
"""
function rich_club_coefficient(g::AbstractGraph, kmax::Int)
    if is_directed(g)
        dgs = sort(indegree(g) .+ outdegree(g))
        rcc = zeros(Float64, kmax)
        for k in 1:min(kmax) - 1
            dgs = dgs[dgs .>= k]
            N_k = length(dgs)
            if N_k < 2
                rcc[k:end] .= Inf
                break
            end
            rcc[k] = sum(dgs) / (N_k * (N_k - 1))
        end
        return rcc
    else
        dgs = sort(outdegree(g))
        rcc = zeros(Float64, kmax)
        for k in 1:min(kmax) - 1
            dgs = dgs[dgs .>= k]
            N_k = length(dgs)
            if N_k < 2
                rcc[k:end] .= Inf
                break
            end
            rcc[k] = 2 * sum(dgs) / (N_k * (N_k - 1))
        end
        return rcc
    end
end

# weighted
function rich_club_coefficient(g::AbstractGraph, kmax::Int, w::AbstractMatrix)
    dstmx = copy(w)
    dgs = indegree(g) .+ outdegree(g)
    rcc = zeros(Float64, kmax)
    wrank = sort(dstmx)
    for k in 1:min(maximum(dgs), kmax)
        E_k = dgs .>= k
        dgs = dgs[E_k]
        dstmx = [E_k, E_k]
        N_k = length(dgs)
        if N_k < 2
            rcc[k:end] .= Inf
            break
        end
        wrank = wrank[1:length(dstmx .!= 0)]
        rcc[k] = sum(dstmx) / sum(wrank)
    end
    return rcc
end

function strength(
    g::AbstractGraph,
    v::Integer,
    dstmx::AbstractMatrix=weights(g);
    normalize::Bool=true
   )
    if is_directed(g)
        return outstrength(g, v, dstmx, normalize=normalize)
    else
        return outstrength(g, v, dstmx, normalize=normalize) +
                instrength(g, v, dstmx, normalize=normalize)
    end
end

function strength(
    g::AbstractGraph,
    vs::AbstractVector=vertices(g),
    dstmx::AbstractMatrix=weights(g);
    normalize::Bool=true
   )
    map(v->strength(g, v, dstmx, normalize=normalize), vs)
end

function instrength(
    g::AbstractGraph,
    vs::AbstractVector=vertices(g),
    dstmx::AbstractMatrix=weights(g);
    normalize::Bool=true,
   )
    map(i -> instrength(g, i, dstmx, normalize=normalize), vs)
end

function instrength(
    g::AbstractGraph,
    v::Integer,
    dstmx::AbstractMatrix=weights(g);
    normalize::Bool=true
   )
    maybe_strength(dstmx[inneighbors(g, v), v], normalize)
end

function outstrength(
    g::AbstractGraph,
    vs::AbstractVector=vertices(g),
    dstmx::AbstractMatrix=weights(g);
    normalize::Bool=true
   )
    map(i -> outstrength(g, i, dstmx, normalize=normalize), vs)
end

function outstrength(
    g::AbstractGraph,
    v::Integer,
    dstmx::AbstractMatrix=weights(g);
    normalize::Bool=true
   )
    maybe_maybe(dstmx[v, outneighbors(g, v)], normalize)
end

function _maybe_normalize(x::AbstractVector, normalize::Bool)
    if normalize
        sum(x) / (length(x) - 1)
    else
        sum(x)
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: route53_recovery_control_config
using AWS.Compat
using AWS.UUIDs

"""
    create_cluster(cluster_name)
    create_cluster(cluster_name, params::Dict{String,<:Any})

Create a new cluster. A cluster is a set of redundant Regional endpoints against which you
can run API calls to update or get the state of one or more routing controls. Each cluster
has a name, status, Amazon Resource Name (ARN), and an array of the five cluster endpoints
(one for each supported Amazon Web Services Region) that you can use with API calls to the
Amazon Route 53 Application Recovery Controller cluster data plane.

# Arguments
- `cluster_name`: The name of the cluster.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientToken"`: Unique client idempotency token.
"""
function create_cluster(ClusterName; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "POST",
        "/cluster",
        Dict{String,Any}("ClusterName" => ClusterName, "ClientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_cluster(
    ClusterName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "POST",
        "/cluster",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClusterName" => ClusterName, "ClientToken" => string(uuid4())
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_control_panel(cluster_arn, control_panel_name)
    create_control_panel(cluster_arn, control_panel_name, params::Dict{String,<:Any})

Creates a new control panel. A control panel represents a group of routing controls that
can be changed together in a single transaction. You can use a control panel to centrally
view the operational status of applications across your organization, and trigger multi-app
failovers in a single transaction, for example, to fail over an Availability Zone or AWS
Region.

# Arguments
- `cluster_arn`: The Amazon Resource Name (ARN) of the cluster for the control panel.
- `control_panel_name`: The name of the control panel.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientToken"`: Unique client idempotency token.
"""
function create_control_panel(
    ClusterArn, ControlPanelName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "POST",
        "/controlpanel",
        Dict{String,Any}(
            "ClusterArn" => ClusterArn,
            "ControlPanelName" => ControlPanelName,
            "ClientToken" => string(uuid4()),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_control_panel(
    ClusterArn,
    ControlPanelName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "POST",
        "/controlpanel",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClusterArn" => ClusterArn,
                    "ControlPanelName" => ControlPanelName,
                    "ClientToken" => string(uuid4()),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_routing_control(cluster_arn, routing_control_name)
    create_routing_control(cluster_arn, routing_control_name, params::Dict{String,<:Any})

Creates a new routing control. A routing control has one of two states: ON and OFF. You can
map the routing control state to the state of an Amazon Route 53 health check, which can be
used to control traffic routing. To get or update the routing control state, see the
Recovery Cluster (data plane) API actions for Amazon Route 53 Application Recovery
Controller.

# Arguments
- `cluster_arn`: The Amazon Resource Name (ARN) of the cluster that includes the routing
  control.
- `routing_control_name`: The name of the routing control.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientToken"`: Unique client idempotency token.
- `"ControlPanelArn"`: The Amazon Resource Name (ARN) of the control panel that includes
  the routing control.
"""
function create_routing_control(
    ClusterArn, RoutingControlName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "POST",
        "/routingcontrol",
        Dict{String,Any}(
            "ClusterArn" => ClusterArn,
            "RoutingControlName" => RoutingControlName,
            "ClientToken" => string(uuid4()),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_routing_control(
    ClusterArn,
    RoutingControlName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "POST",
        "/routingcontrol",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ClusterArn" => ClusterArn,
                    "RoutingControlName" => RoutingControlName,
                    "ClientToken" => string(uuid4()),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_safety_rule()
    create_safety_rule(params::Dict{String,<:Any})

Creates a safety rule in a control panel. Safety rules let you add safeguards around
enabling and disabling routing controls, to help prevent unexpected outcomes. There are two
types of safety rules: assertion rules and gating rules. Assertion rule: An assertion rule
enforces that, when a routing control state is changed, the criteria set by the rule
configuration is met. Otherwise, the change to the routing control is not accepted. Gating
rule: A gating rule verifies that a set of gating controls evaluates as true, based on a
rule configuration that you specify. If the gating rule evaluates to true, Amazon Route 53
Application Recovery Controller allows a set of routing control state changes to run and
complete against the set of target controls.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AssertionRule"`:
- `"ClientToken"`: Unique client idempotency token.
- `"GatingRule"`:
"""
function create_safety_rule(; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "POST",
        "/safetyrule",
        Dict{String,Any}("ClientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_safety_rule(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "POST",
        "/safetyrule",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ClientToken" => string(uuid4())), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_cluster(cluster_arn)
    delete_cluster(cluster_arn, params::Dict{String,<:Any})

Delete a cluster.

# Arguments
- `cluster_arn`: The Amazon Resource Name (ARN) of the cluster that you're deleting.

"""
function delete_cluster(ClusterArn; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "DELETE",
        "/cluster/$(ClusterArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_cluster(
    ClusterArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "DELETE",
        "/cluster/$(ClusterArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_control_panel(control_panel_arn)
    delete_control_panel(control_panel_arn, params::Dict{String,<:Any})

Deletes a control panel.

# Arguments
- `control_panel_arn`: The Amazon Resource Name (ARN) of the control panel that you're
  deleting.

"""
function delete_control_panel(
    ControlPanelArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "DELETE",
        "/controlpanel/$(ControlPanelArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_control_panel(
    ControlPanelArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "DELETE",
        "/controlpanel/$(ControlPanelArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_routing_control(routing_control_arn)
    delete_routing_control(routing_control_arn, params::Dict{String,<:Any})

Deletes a routing control.

# Arguments
- `routing_control_arn`: The Amazon Resource Name (ARN) of the routing control that you're
  deleting.

"""
function delete_routing_control(
    RoutingControlArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "DELETE",
        "/routingcontrol/$(RoutingControlArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_routing_control(
    RoutingControlArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "DELETE",
        "/routingcontrol/$(RoutingControlArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_safety_rule(safety_rule_arn)
    delete_safety_rule(safety_rule_arn, params::Dict{String,<:Any})

Deletes a safety rule./>

# Arguments
- `safety_rule_arn`: The request body that you include when you update a safety rule.

"""
function delete_safety_rule(
    SafetyRuleArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "DELETE",
        "/safetyrule/$(SafetyRuleArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_safety_rule(
    SafetyRuleArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "DELETE",
        "/safetyrule/$(SafetyRuleArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_cluster(cluster_arn)
    describe_cluster(cluster_arn, params::Dict{String,<:Any})

Display the details about a cluster. The response includes the cluster name, endpoints,
status, and Amazon Resource Name (ARN).

# Arguments
- `cluster_arn`: The Amazon Resource Name (ARN) of the cluster that you're getting details
  for.

"""
function describe_cluster(ClusterArn; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "GET",
        "/cluster/$(ClusterArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_cluster(
    ClusterArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/cluster/$(ClusterArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_control_panel(control_panel_arn)
    describe_control_panel(control_panel_arn, params::Dict{String,<:Any})

Displays details about a control panel.

# Arguments
- `control_panel_arn`: The Amazon Resource Name (ARN) of the control panel that you're
  getting details for.

"""
function describe_control_panel(
    ControlPanelArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanel/$(ControlPanelArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_control_panel(
    ControlPanelArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanel/$(ControlPanelArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_routing_control(routing_control_arn)
    describe_routing_control(routing_control_arn, params::Dict{String,<:Any})

Displays details about a routing control. A routing control has one of two states: ON and
OFF. You can map the routing control state to the state of an Amazon Route 53 health check,
which can be used to control routing. To get or update the routing control state, see the
Recovery Cluster (data plane) API actions for Amazon Route 53 Application Recovery
Controller.

# Arguments
- `routing_control_arn`: The Amazon Resource Name (ARN) of the routing control that you're
  getting details for.

"""
function describe_routing_control(
    RoutingControlArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/routingcontrol/$(RoutingControlArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_routing_control(
    RoutingControlArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/routingcontrol/$(RoutingControlArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_safety_rule(safety_rule_arn)
    describe_safety_rule(safety_rule_arn, params::Dict{String,<:Any})

Describes the safety rules (that is, the assertion rules and gating rules) for the routing
controls in a control panel.

# Arguments
- `safety_rule_arn`: The request body that you include when you update a safety rule.

"""
function describe_safety_rule(
    SafetyRuleArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/safetyrule/$(SafetyRuleArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_safety_rule(
    SafetyRuleArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/safetyrule/$(SafetyRuleArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_associated_route53_health_checks(routing_control_arn)
    list_associated_route53_health_checks(routing_control_arn, params::Dict{String,<:Any})

Returns an array of all Amazon Route 53 health checks associated with a specific routing
control.

# Arguments
- `routing_control_arn`: The Amazon Resource Name (ARN) of the routing control that you're
  getting details for.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The number of objects that you want to return with this call.
- `"NextToken"`: The token that identifies which batch of results you want to see.
"""
function list_associated_route53_health_checks(
    RoutingControlArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/routingcontrol/$(RoutingControlArn)/associatedRoute53HealthChecks";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_associated_route53_health_checks(
    RoutingControlArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/routingcontrol/$(RoutingControlArn)/associatedRoute53HealthChecks",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_clusters()
    list_clusters(params::Dict{String,<:Any})

Returns an array of all the clusters in an account.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The number of objects that you want to return with this call.
- `"NextToken"`: The token that identifies which batch of results you want to see.
"""
function list_clusters(; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "GET", "/cluster"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_clusters(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET", "/cluster", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_control_panels()
    list_control_panels(params::Dict{String,<:Any})

Returns an array of control panels for a cluster.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClusterArn"`: The Amazon Resource Name (ARN) of a cluster.
- `"MaxResults"`: The number of objects that you want to return with this call.
- `"NextToken"`: The token that identifies which batch of results you want to see.
"""
function list_control_panels(; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "GET", "/controlpanels"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_control_panels(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanels",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_routing_controls(control_panel_arn)
    list_routing_controls(control_panel_arn, params::Dict{String,<:Any})

Returns an array of routing controls for a control panel. A routing control is an Amazon
Route 53 Application Recovery Controller construct that has one of two states: ON and OFF.
You can map the routing control state to the state of an Amazon Route 53 health check,
which can be used to control routing.

# Arguments
- `control_panel_arn`: The Amazon Resource Name (ARN) of the control panel that you're
  getting routing control details for.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The number of objects that you want to return with this call.
- `"NextToken"`: The token that identifies which batch of results you want to see.
"""
function list_routing_controls(
    ControlPanelArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanel/$(ControlPanelArn)/routingcontrols";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_routing_controls(
    ControlPanelArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanel/$(ControlPanelArn)/routingcontrols",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_safety_rules(control_panel_arn)
    list_safety_rules(control_panel_arn, params::Dict{String,<:Any})

List the safety rules (the assertion rules and gating rules) that you've defined for the
routing controls in a control panel.

# Arguments
- `control_panel_arn`: The Amazon Resource Name (ARN) of the control panel that you're
  getting details for.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The number of objects that you want to return with this call.
- `"NextToken"`: The token that identifies which batch of results you want to see.
"""
function list_safety_rules(
    ControlPanelArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanel/$(ControlPanelArn)/safetyrules";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_safety_rules(
    ControlPanelArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "GET",
        "/controlpanel/$(ControlPanelArn)/safetyrules",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_control_panel(control_panel_arn, control_panel_name)
    update_control_panel(control_panel_arn, control_panel_name, params::Dict{String,<:Any})

Updates a control panel. The only update you can make to a control panel is to change the
name of the control panel.

# Arguments
- `control_panel_arn`: The Amazon Resource Name (ARN) of the control panel.
- `control_panel_name`: The name of the control panel.

"""
function update_control_panel(
    ControlPanelArn, ControlPanelName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "PUT",
        "/controlpanel",
        Dict{String,Any}(
            "ControlPanelArn" => ControlPanelArn, "ControlPanelName" => ControlPanelName
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_control_panel(
    ControlPanelArn,
    ControlPanelName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "PUT",
        "/controlpanel",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ControlPanelArn" => ControlPanelArn,
                    "ControlPanelName" => ControlPanelName,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_routing_control(routing_control_arn, routing_control_name)
    update_routing_control(routing_control_arn, routing_control_name, params::Dict{String,<:Any})

Updates a routing control. You can only update the name of the routing control. To get or
update the routing control state, see the Recovery Cluster (data plane) API actions for
Amazon Route 53 Application Recovery Controller.

# Arguments
- `routing_control_arn`: The Amazon Resource Name (ARN) of the routing control.
- `routing_control_name`: The name of the routing control.

"""
function update_routing_control(
    RoutingControlArn, RoutingControlName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "PUT",
        "/routingcontrol",
        Dict{String,Any}(
            "RoutingControlArn" => RoutingControlArn,
            "RoutingControlName" => RoutingControlName,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_routing_control(
    RoutingControlArn,
    RoutingControlName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return route53_recovery_control_config(
        "PUT",
        "/routingcontrol",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "RoutingControlArn" => RoutingControlArn,
                    "RoutingControlName" => RoutingControlName,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_safety_rule()
    update_safety_rule(params::Dict{String,<:Any})

Update a safety rule (an assertion rule or gating rule) for the routing controls in a
control panel. You can only update the name and the waiting period for a safety rule. To
make other updates, delete the safety rule and create a new safety rule.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"AssertionRuleUpdate"`:
- `"GatingRuleUpdate"`:
"""
function update_safety_rule(; aws_config::AbstractAWSConfig=global_aws_config())
    return route53_recovery_control_config(
        "PUT", "/safetyrule"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function update_safety_rule(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return route53_recovery_control_config(
        "PUT", "/safetyrule", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end


text:
"""
Atomic solution for an SDE.

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `ΔW`: Wiener process driving the stochastic process q
* `ΔZ`: Wiener process driving the stochastic process q
* `K`:  integer parameter defining the truncation of the increments of the Wiener process (for strong solutions)
"""
mutable struct AtomicSolutionSDE{DT,TT} <: AtomicSolution{DT,TT}
    t::TT
    t̅::TT

    q::Vector{DT}
    q̅::Vector{DT}
    q̃::Vector{DT}

    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    K::Int

    function AtomicSolutionSDE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real}
        new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nm), zeros(DT, nm))
    end
end

AtomicSolutionSDE(DT, TT, nd, nm) = AtomicSolutionSDE{DT, TT}(nd, nm)

function set_solution!(asol::AtomicSolutionSDE, sol)
    t, q = sol
    asol.t  = t
    asol.q .= q
end

function set_increments!(asol::AtomicSolutionSDE, incs)
    ΔW, ΔZ = incs
    asol.ΔW .= ΔW
    asol.ΔZ .= ΔZ
end

function get_solution(asol::AtomicSolutionSDE)
    (asol.t, asol.q)
end

function get_increments(asol::AtomicSolutionSDE)
    (asol.ΔW, asol.ΔZ)
end

function get_increments!(asol::AtomicSolutionSDE, ΔW, ΔZ)
    ΔW .= asol.ΔW
    ΔZ .= asol.ΔZ
end

function CommonFunctions.reset!(asol::AtomicSolutionSDE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.t += Δt
end

function update!(asol::AtomicSolutionSDE{DT}, y::Vector{DT}) where {DT}
    for k in eachindex(y)
        update!(asol, y[k], k)
    end
end

function update!(asol::AtomicSolutionSDE{DT}, y::DT, k::Int) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(y, asol.q[k], asol.q̃[k])
end


text:
@doc raw"""
    MaternKernel([ν=1 [,θ=1]])

The Matern kernel is a Mercer kernel given by:

```math
\kappa(\mathbf{x},\mathbf{y}) =
\frac{1}{2^{\nu-1}\Gamma(\nu)}
\left(\frac{2\sqrt{\nu}||\mathbf{x}-\mathbf{y}||}{\theta}\right)^{\nu}
K_{\nu}\left(\frac{2\sqrt{\nu}||\mathbf{x}-\mathbf{y}||}{\theta}\right)
```

where ``\Gamma`` is the gamma function, ``K_{\nu}`` is the modified Bessel function of the
second kind, ``\nu > 0`` and ``\theta > 0``.
"""
struct MaternKernel{T<:AbstractFloat} <: MercerKernel{T}
    nu::HyperParameter{T}
    rho::HyperParameter{T}
    MaternKernel{T}(ν::Real, ρ::Real) where {T<:AbstractFloat}  = new{T}(
        HyperParameter(convert(T,ν), interval(OpenBound(zero(T)), nothing)),
        HyperParameter(convert(T,ρ), interval(OpenBound(zero(T)), nothing))
    )
end
function MaternKernel(ν::T1=1.0, ρ::T2=one(T1)) where {T1<:Real,T2<:Real}
    MaternKernel{floattype(T1,T2)}(ν,ρ)
end

@inline function maternkernel(z::T, ν::T, ρ::T) where {T}
    v1 = sqrt(2ν) * z / ρ
    v1 = v1 < eps(T) ? eps(T) : v1  # Overflow risk as z -> Inf
    2 * (v1/2)^(ν) * besselk(ν, v1) / gamma(ν)
end

@inline basefunction(::MaternKernel) = SquaredEuclidean()
@inline function kappa(κ::MaternKernel{T}, z::T) where {T}
    maternkernel(z, getvalue(κ.nu), getvalue(κ.rho))
end


text:
"""
`SpatialPooler` is a learning algorithm that decorrelates the features of an input space,
producing a Sparse Distributed Representation (SDR) of the input space.
If defines the proximal connections of an HTM layer.

# Examples
```julia
sp= SpatialPooler(SPParams(szᵢₙ=(600,), szₛₚ=(2048,)))
z= Random.bitrand(600)  # input
activation= sp(z)
# or, to adapt the SP to the input as well:
activation= step!(sp,z)
```

# Properties

It's called `SpatialPooler` because input patterns that share a large number of co-active neurons
(i.e., that are spatially similar) are grouped together into a common output representation.
It is designed to achieve a set of _computational properties_ that support further downstream computations with SDRs, including:

- preserving topology of the input space by mapping similar inputs to similar outputs
- continuously adapting to changing statistics of the input stream
- forming fixed sparsity representations
- being robust to noise
- being fault tolerant

[Source](https://www.frontiersin.org/articles/10.3389/fncom.2017.00111/full)

# Algorithm overview

## Mapping I/O spaces

The spatial pooler maps an input space `x` to an output space `y` of minicolumns through a matrix of proximal synapses.
The input space can optionally have a topology, which the spatial pooler will preserve by mapping
output minicolumn `yᵢ` to a subset of the input space, a [`Hypercube`](@ref) around center `xᶜ`.
```julia
xᶜ(yᵢ)= floor.(Int, (yᵢ.-1) .* (szᵢₙ./szₛₚ)) .+1
````

## Output activation

1. Calculate the overlap `o(yᵢ)` by propagating the input activations through the proximal synapses
   and adjusting by boosting factors `b` (control mechanism that spreads out the activation pattern
   across understimulated neurons)
2. Inhibition `Z` between `yᵢ` (local/global): competition where only the top-K `yᵢ` with the highest `o(yᵢ)` win; ensures sparsity
3. Activate winning `yᵢ` > activation threshold (`θ_stimulus_activate`)

See also: [`sp_activate`](@ref)

# State variables

- `synapses`: includes the synapse permanence matrix Dₚ
- `åₜ`: [boosting] moving (in time) average (in time) activation of each minicolumn
- `åₙ`: [boosting] moving (in time) average (in neighborhood) activation of each minicolumn
- `φ`: the adaptible radius of local inhibition

---
See also: [`ProximalSynapses`](@ref), [`TemporalMemory`](@ref)
"""
struct SpatialPooler
  params::SPParams
  # includes the permanence matrix Dₚ
  synapses::ProximalSynapses
  φ::Vector{Float32}  # vector instead of Float32 for mutating
  # Boosting
  "boosting factors"
  b::Vector{Float32}
  "[boosting] average in time activation of each minicolumn"
  åₜ::Array{Float32}
  "[boosting] average in neighborhood activation of each minicolumn"
  åₙ::Array{Float32}
end
function SpatialPooler(params::SPParams)
  params= normalize_SPparams(params)
  @unpack szᵢₙ,szₛₚ,prob_synapse,θ_permanence,γ,
          enable_local_inhibit = params

  synapseSparsity= prob_synapse * (enable_local_inhibit ?
                      (α(γ)^length(szᵢₙ))/prod(szᵢₙ) : 1)
  SpatialPooler(params,
      ProximalSynapses(szᵢₙ,szₛₚ,synapseSparsity,γ,
          prob_synapse,θ_permanence, topology= enable_local_inhibit),
      init_φ(γ,prob_synapse,szᵢₙ,szₛₚ),
      ones(prod(szₛₚ)), zeros(szₛₚ), zeros(szₛₚ)
  )
end
init_φ(γ,prob_synapse,szᵢₙ,szₛₚ)= [((2γ+0.5)*prob_synapse*mean(szₛₚ./szᵢₙ) - 1)/2]
b(sp::SpatialPooler)= sp.params.enable_boosting ? ones(length(sp.b)) : sp.b
åₙ(sp::SpatialPooler)= sp.åₙ
Wₚ(sp::SpatialPooler)= Wₚ(sp.synapses)
φ(sp::SpatialPooler)= sp.φ[1]

reset!(sp::SpatialPooler)= begin
  @unpack szᵢₙ,szₛₚ,prob_synapse,γ = sp.params
  reset!(sp.synapses)
  sp.φ.= init_φ(γ,prob_synapse,szᵢₙ,szₛₚ)
  sp.b.= 1.0
  sp.åₜ.= 0.0
  sp.åₙ.= 0.0
end

# SP activation

"""
Return the activation pattern of the [`SpatialPooler`](@ref) for the given input activation.

`z`: [`CellActivity`](@ref)

# Example
```julia
sp= SpatialPooler(SPParams(szᵢₙ=(600,), szₛₚ=(2048,)))
z= Random.bitrand(600)  # input
activation= sp(z)
# or, to adapt the SP to the input as well:
activation= step!(sp,z)
```

For details see: [`sp_activate`](@ref)
"""
(sp::SpatialPooler)(z)= sp_activate(sp,z)

"""
`sp_activate(sp::SpatialPooler, z)` calculates the SP's output activation for given input activation `z`.

# Algorithm

1. Overlap `o(yᵢ)`
    - Propagate the input activations through the proximal synapses (matrix multiply)
    - Apply boosting factors `b`: control mechanism that spreads out the activation pattern across understimulated
      neurons (homeostatic excitability control)
2. Inhibition `Z` between `yᵢ` (local/global): competition where only the top-k `yᵢ` with the highest `o(yᵢ)` win; ensures sparsity
   The competition happens within an area around each neuron.
    - `k`: number of winners depends on desired sparsity (`s`, see [`SPParams`](@ref)) and area size
    - `θ_inhibit`: inhibition threshold per neighborhood ``:= o(k-th y)``
    - `Z(o(yᵢ))`: convolution of `o` with `θ_inhibit`
3. Activate winning `yᵢ` > activation threshold (`θ_stimulus_activate`)

!!! info
    The local sparsity `s` tends to the sparsity of the entire layer as it grows larger, but for small values or small inhibition radius
    it diverges, because of the limited & integral number of neurons winning in each neighborhood.
    This could be addressed by *tie breaking*, but it doesn't seem to have much practical importance.
"""
function sp_activate(sp::SpatialPooler, z)
  @unpack szₛₚ,s,θ_permanence,θ_stimulus_activate,enable_local_inhibit = sp.params

  # overlap
  overlap(z)= @chain vec(z) begin
    Wₚ(sp)' * _
    θ_activation.(_)
    b(sp) .* _
    reshape(szₛₚ)
  end
  θ_activation(x)= x .> θ_stimulus_activate ? x : 0

  # inhibition
  area()= enable_local_inhibit ? α(φ(sp))^length(szₛₚ) : prod(szₛₚ)
  k()=    ceil(Int, s*area())
  # inhibition threshold per area
  # OPTIMIZE: local inhibition is the SP's bottleneck. "mapwindow" is suboptimal;
  #   https://github.com/JuliaImages/Images.jl/issues/751
  θ_inhibit!(X)= @> X vec partialsort!(k(),rev=true)
  # Z: k-th largest overlap in neighborhood
  Z(y)= _Z(Val(enable_local_inhibit),y)
  # we increase the threshold a bit to account stochastically for the error in approximating k and for tiebreaking
  _Z(loc_inhibit::Val{true}, y)= @> mapwindow(θ_inhibit!, y, neighborhood(φ(sp),length(szₛₚ)), border=Fill(0))
  _Z(loc_inhibit::Val{false},y)= @> θ_inhibit!(copy(y))
  Z_perm(y)= partialsortperm(y, 1:k(), rev=true)

  # Tiebreaker
  sErr()= 1 - k() + s*area()
  tieEst(o,Z)= count(o .== Z) * area()/prod(szₛₚ)
  tiebreaker(o,Z)= rand(Float32,szₛₚ) .- (1 - sErr()/tieEst(o,Z))
  activate_tiebreaker(o)= begin
    Z= Z(o)
    (o .+ tiebreaker(o,Z) .>= Z)
  end

  if enable_local_inhibit
    @chain z begin
      overlap
      activate_tiebreaker
    end
  else
    @chain z begin
      overlap
      Z_perm
      bitarray(_, szₛₚ)
    end
  end
end

# SP adaptation

"""
`step!(sp::SpatialPooler, z::CellActivity)` evolves the Spatial Pooler to the next timestep
by evolving each of its constituents (synapses, boosting, inhibition radius)
and returns the output activation.

See also: [`sp_activate`](@ref)
"""
function step!(sp::SpatialPooler, z)
  a= sp_activate(sp, z)
  if sp.params.enable_learning
    step!(sp.synapses, z,a, sp.params)
    step_boost!(sp,a)
    # φ kept static for the moment (no obvious benefit from dynamic φ)
    #nupic_update_φ!(sp)
  end
  return a
end

# # boosting
"""
`step_boost!(sp::SpatialPooler,a)` evolves the boosting factors `b` (see [`sp_activate`](@ref)).
They depend on:
- `åₜ`: moving average in time activation of each minicolumn
- `åₙ`: moving average in neighborhood activation of each minicolum
"""
step_boost!(sp::SpatialPooler,a)= begin
  @unpack Tboost,β, szₛₚ= sp.params
  step_åₙ!(sp)
  # exponential moving average with period Tboost
  sp.åₜ.= (sp.åₜ.*(Tboost-1) .+ reshape(a,szₛₚ))./Tboost
  sp.b.= boostfun(sp.åₜ,sp.åₙ,β)
end
boostfun(åₜ,åₙ,β)= @> exp.(-β .* (åₜ .- åₙ)) vec
step_åₙ!(sp::SpatialPooler)= begin
  # local mean filter
  _step_åₙ!(sp::SpatialPooler,local_inhibit::Val{true})=
      imfilter!(sp.åₙ,sp.åₜ, mean_kernel(φ(sp),length(sp.params.szₛₚ)), "symmetric")
  _step_åₙ!(sp::SpatialPooler,local_inhibit::Val{false})= (sp.åₙ.= mean(sp.åₜ))

  _step_åₙ!(sp,Val(sp.params.enable_local_inhibit))
end

# # Inhibition radius

# This implementation follows the SP paper description and NUPIC, but seems too complex
#   for no reason. Replace with a static inhibition radius instead
nupic_update_φ!(sp::SpatialPooler)= begin
  @unpack szₛₚ, szᵢₙ, prob_synapse = sp.params
  # sample 10% of minicolumns
  mean_receptiveFieldSpan()= (i-> receptivefieldSpan(szᵢₙ,Wₚ(sp)[:,i]) ).( randperm(size(Wₚ(sp),2))[1:100] ) |> mean
  #mean_receptiveFieldSpan()= mapslices(receptivefieldSpan, W, dims=1)|> mean
  # Simplified version:
  #mean_receptiveFieldSpan()= (2γ+0.5) * prob_synapse
  receptiveFieldSpan_yspace()= ( mean_receptiveFieldSpan()*mean(szₛₚ./szᵢₙ) - 1)/2
  sp.φ[1]= max(receptiveFieldSpan_yspace(), 1)
end

# Foreach minicolumn, find the geometric span of its projection to the input space through the connected synapses
# foreach dimension, then return the mean
function receptivefieldSpan(sz_in, Wᵢ)
  c= CartesianIndices(sz_in)
  connectedInputXY= c[findall(Wᵢ)]
  maxc= [mapreduce(c->c.I[d], max, connectedInputXY) for d in 1:length(sz_in)]
  minc= [mapreduce(c->c.I[d], min, connectedInputXY) for d in 1:length(sz_in)]
  mean(maxc .- minc .+ 1)
end

normalize_SPparams(params)= begin
  params.szᵢₙ|> typeof <: Int ? (params = @set params.szᵢₙ = (params.szᵢₙ,)) : nothing
  params.szₛₚ|> typeof <: Int ? (params = @set params.szₛₚ = (params.szₛₚ,)) : nothing
  !params.enable_local_inhibit ? (params = @set params.γ = maximum(params.szᵢₙ)) : nothing
  params
end


text:
# This file was generated, do not modify it. # hide
v = categorical([1, 2, 3, 1, 2, 3, 1, 2, 3], ordered=true)

levels(v)


text:
module LinearAlgebra

#println("You have decided to learn Linear Algebra using Julia. Please refer the Julia_LinAlg.pdf")

end


text:
"""
Trains SeqModel model given a config dict.

Returns dict containing results.
"""
function easy_sample!(
        model::Union{SeqModel,DistributedSeqModel},
        spikes::Vector{Spike},
        initial_assignments::Vector{Int64},
        config::Dict
    )

    # Save copy of initial assignments.
    _inits = copy(initial_assignments)

    # Draw annealed Gibbs samples.
    (
        assignments,
        anneal_assignment_hist,
        anneal_log_p_hist,
        anneal_latent_event_hist,
        anneal_globals_hist
    ) =
    annealed_gibbs!(
        model,
        spikes,
        initial_assignments,
        config[:num_anneals],
        config[:samples_per_anneal],
        config[:max_temperature],
        config[:split_merge_moves_during_anneal],
        config[:split_merge_window],
        config[:save_every_during_anneal];
        verbose=true
    )

    # # Sanity check.
    # for k in model.sequence_events.indices
    #     event = model.sequence_events[k]
    #     @assert event.spike_count == sum(assignments .== k)
    # end

    # Draw regular Gibbs samples.
    (
        assignments,
        assignment_hist,
        log_p_hist,
        latent_event_hist,
        globals_hist
    ) =
    gibbs_sample!(
        model,
        spikes,
        assignments,
        config[:samples_after_anneal],
        config[:split_merge_moves_after_anneal],
        config[:split_merge_window],
        config[:save_every_after_anneal];
        verbose=true
    )

    # # Sanity check.
    # for k in model.sequence_events.indices
    #     event = model.sequence_events[k]
    #     @assert event.spike_count == sum(assignments .== k)
    # end

    return Dict(

        # Initial assignment variables.
        :initial_assignments => _inits,

        # Results during annealing.
        :anneal_assignment_hist => anneal_assignment_hist,
        :anneal_log_p_hist => anneal_log_p_hist,
        :anneal_latent_event_hist => anneal_latent_event_hist,
        :anneal_globals_hist => anneal_globals_hist,

        # Results after annealing.
        :final_assignments => assignments,
        :assignment_hist => assignment_hist,
        :latent_event_hist => latent_event_hist,
        :globals_hist => globals_hist,
        :log_p_hist => log_p_hist
    )
end