text:
module RulesNormalMeanVarianceOutTest

using Test
using ReactiveMP
using Random
using Distributions

import ReactiveMP: @test_rules

@testset "rules:NormalMeanVariance:out" begin

    @testset "Belief Propagation: (m_μ::PointMass, m_v::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = PointMass(-1.0), m_v = PointMass(2.0)), output = NormalMeanVariance(-1.0, 2.0)),
            (input = (m_μ = PointMass(1.0), m_v = PointMass(2.0)),  output = NormalMeanVariance(1.0, 2.0)),
            (input = (m_μ = PointMass(2.0), m_v = PointMass(1.0)),  output = NormalMeanVariance(2.0, 1.0))
        ]

    end

    @testset "Belief Propagation: (m_μ::UnivariateNormalDistributionsFamily, m_v::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = NormalMeanVariance(0.0, 1.0),  m_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_μ = NormalMeanVariance(-1.0, 1.0), m_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_μ = NormalMeanVariance(2.0, 0.5),  m_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 1.5)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = NormalMeanPrecision(0.0, 1.0),  m_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_μ = NormalMeanPrecision(-1.0, 1.0), m_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_μ = NormalMeanPrecision(2.0, 0.5),  m_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = NormalWeightedMeanPrecision(0.0, 1.0),  m_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_μ = NormalWeightedMeanPrecision(-1.0, 1.0), m_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_μ = NormalWeightedMeanPrecision(2.0, 0.5),  m_v = PointMass(1.0)), output = NormalMeanVariance(4.0, 3.0)),
        ]

    end

    @testset "Belief Propagation: (m_μ::UnivariateNormalDistributionsFamily, q_v::PointMass)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = NormalMeanVariance(0.0, 1.0),  q_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_μ = NormalMeanVariance(-1.0, 1.0), q_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_μ = NormalMeanVariance(2.0, 0.5),  q_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 1.5)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = NormalMeanPrecision(0.0, 1.0),  q_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_μ = NormalMeanPrecision(-1.0, 1.0), q_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_μ = NormalMeanPrecision(2.0, 0.5),  q_v = PointMass(1.0)), output = NormalMeanVariance(2.0, 3.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (m_μ = NormalWeightedMeanPrecision(0.0, 1.0),  q_v = PointMass(2.0)), output = NormalMeanVariance(0.0, 3.0)),
            (input = (m_μ = NormalWeightedMeanPrecision(-1.0, 1.0), q_v = PointMass(1.5)), output = NormalMeanVariance(-1.0, 2.5)),
            (input = (m_μ = NormalWeightedMeanPrecision(2.0, 0.5),  q_v = PointMass(1.0)), output = NormalMeanVariance(4.0, 3.0)),
        ]

    end

    @testset "Variational: (q_μ::Any, q_v::Any)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (q_μ = PointMass(-1.0), q_v = PointMass(2.0)), output = NormalMeanVariance(-1.0, 2.0)),
            (input = (q_μ = PointMass(1.0), q_v = PointMass(2.0)),  output = NormalMeanVariance(1.0, 2.0)),
            (input = (q_μ = PointMass(2.0), q_v = PointMass(1.0)),  output = NormalMeanVariance(2.0, 1.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (q_μ = NormalMeanVariance(-1.0, 2.0), q_v = PointMass(2.0)), output = NormalMeanVariance(-1.0, 2.0)),
            (input = (q_μ = NormalMeanPrecision(1.0, 4.0), q_v = PointMass(3.0)),  output = NormalMeanVariance(1.0, 3.0)),
            (input = (q_μ = NormalWeightedMeanPrecision(2.0, 4.0), q_v = PointMass(1.0)),  output = NormalMeanVariance(0.5, 1.0))
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanVariance(:out, Marginalisation) [
            (input = (q_μ = PointMass(-1.0), q_v = InverseGamma(2.0, 1.0)), output = NormalMeanVariance(-1.0, 1.0)),
            (input = (q_μ = PointMass(1.0), q_v = InverseGamma(4.0, 2.0)),  output = NormalMeanVariance(1.0, 2.0 / 3.0)),
            (input = (q_μ = PointMass(2.0), q_v = InverseGamma(4.0, 6.0)),  output = NormalMeanVariance(2.0, 2.0))
        ]

    end

end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to including a table in a file.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Include Pretty Table to file" begin

    # Text
    # ==========================================================================

    path = "test.txt"

    orig = """
    This is one line.

    This is "another" line.

    % <PrettyTables Table 1> This should be deleted.
    This should be deleted.
    This should be deleted.
    This should be deleted.
    This should be deleted.

    % </PrettyTables>

    <PrettyTables Table 2></PrettyTables>

    <PrettyTables Table 3>This should not be deleted.
    This should not be deleted.
    This should not be deleted.
    This should not be deleted."""

    open(path,"w") do f
        write(f, orig)
    end

    data_table_1 = [1 2 3
                    4 5 6]

    data_table_2 = [7 8 9
                    1 2 3]

    include_pt_in_file(path, "Table 2", data_table_2, tf = tf_mysql,
                       body_hlines = [1])
    include_pt_in_file(path, "Table 1", data_table_1, alignment = :c,
                       show_row_number = true, backup_file = false)
    include_pt_in_file(path, "Table 3", data_table_2, alignment = :c,
                       show_row_number = true, backup_file = false)

    result = read(path, String)
    backup = read(path * "_backup", String)

    expected = """
    This is one line.

    This is "another" line.

    % <PrettyTables Table 1>
    ┌─────┬────────┬────────┬────────┐
    │ Row │ Col. 1 │ Col. 2 │ Col. 3 │
    ├─────┼────────┼────────┼────────┤
    │   1 │   1    │   2    │   3    │
    │   2 │   4    │   5    │   6    │
    └─────┴────────┴────────┴────────┘
    % </PrettyTables>

    <PrettyTables Table 2>
    +--------+--------+--------+
    | Col. 1 | Col. 2 | Col. 3 |
    +--------+--------+--------+
    |      7 |      8 |      9 |
    +--------+--------+--------+
    |      1 |      2 |      3 |
    +--------+--------+--------+
    </PrettyTables>

    <PrettyTables Table 3>This should not be deleted.
    This should not be deleted.
    This should not be deleted.
    This should not be deleted."""

    @test backup == orig
    @test result == expected

    # HTML (`tag_append` option)
    # ==========================================================================

    data_table_1 = [1 2 3
                    4 5 6]

    path = "test.html"

    orig = """

This is a table:

    """

    open(path,"w") do f
        write(f, orig)
    end

    data_table_1 = [1 2 3
                    4 5 6]

    include_pt_in_file(path, "Table 1", data_table_1, backend = :html,
                       standalone = false, tag_append = " -->")

    result = read(path, String)
    backup = read(path * "_backup", String)

    expected = """

This is a table:

        Col. 1
        Col. 2
        Col. 3

        1
        2
        3

        4
        5
        6

    """

    @test backup == orig
    @test result == expected

    # Markdown (`remove_tags` option)
    # ==========================================================================

    data_table_1 = [1 2 3
                    4 5 6]

    path = "test.md"

    orig = """
    # Markdown

    This is a markdown table.

    <PrettyTables Table 1> This should be removed.
    This should be removed.
    This should be removed.
    This should be removed.
    This should be removed.
    </PrettyTables>
    """

    open(path,"w") do f
        write(f, orig)
    end

    data_table_1 = [1 2 3
                    4 5 6]

    include_pt_in_file(path, "Table 1", data_table_1, tf = tf_markdown,
                       remove_tags = true)
    result = read(path, String)
    backup = read(path * "_backup", String)

    expected = """
    # Markdown

    This is a markdown table.

    | Col. 1 | Col. 2 | Col. 3 |
    |--------|--------|--------|
    |      1 |      2 |      3 |
    |      4 |      5 |      6 |

    """

    @test backup == orig
    @test result == expected
end

ext: jl


text:
module Readables

export Readable,
       readable, readablestring,
       decpoint, setdecpoint,
       intsep, setintsep, intgroup, setintgroup,
       fracsep, setfracsep, fracgroup, setfracgroup

const IMAG_UNIT_STR = ["𝛊"]
const DUAL_UNIT_STR = ["ε"]

const IntGroupSize  = Ref(3)
const FracGroupSize = Ref(5)
const IntSepChar    = Ref(',')
const FracSepChar   = Ref('_')
const DecimalPoint  = Ref('.')

mutable struct Readable
    intgroupsize::Int
    fracgroupsize::Int
    intsepchar::Char
    fracsepchar::Char
    decimalpoint::Char

    function Readable(intgroupsize::Int, fracgroupsize::Int, intsepchar::Char, fracsepchar::Char, decimalpoint::Char)
        if !(0 < intgroupsize && 0 < fracgroupsize)
            throw(ErrorException("groups must be > 0 ($intgroupsize, $fracgroupsize)"))
        end
        return new(intgroupsize, fracgroupsize, intsepchar, fracsepchar, decimalpoint)
    end
end

Readable(;intgroup::Int=IntGroupSize[],
          fracgroup::Int=FracGroupSize[],
          intsep::Char=IntSepChar[],
          fracsep::Char=FracSepChar[],
          decimalpoint::Char=DecimalPoint[]
          ) =
    Readable(intgroup, fracgroup, intsep, fracsep, decimalpoint)

const baseprefixes = Dict(2=>"0b", 8=>"0o", 10=>"", 16=>"0x")

function baseprefix(x::Int)
    res = get(baseprefixes, x, nothing)
    res === nothing && throw(ErrorException("base $x is not supported"))
    return res
end

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

const READABLE = Readable()

decpoint(x::Readable)  = x.decimalpoint
intsep(x::Readable)    = x.intsepchar
intgroup(x::Readable)  = x.intgroupsize
fracsep(x::Readable)   = x.fracsepchar
fracgroup(x::Readable) = x.fracgroupsize

function setdecpoint(x::Readable, decpoint::Char)
    x.decimalpoint = decpoint
    return x
end
function setintsep(x::Readable, intsep::Char)
    x.intsepchar = intsep
    return x
end
function setintgroup(x::Readable, intgroup::Int)
    x.intgroupsize = intgroup
    return x
end
function setfracsep(x::Readable, fracsep::Char)
    x.fracsepchar = fracsep
    return x
end
function setfracgroup(x::Readable, fracgroup::Int)
    x.fracgroupsize = fracgroup
    return x
end

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

function readablestring(x::T; base::Int=10, sepwith::Union{Nothing,Char}=nothing, groupby::Union{Nothing,Int}=nothing) where {T}
    if isnothing(sepwith)
        if isnothing(groupby)
            readablestring(x, base=base)
        else
            readablestring(Readable(groupby=groupby, sepwith = intsep(READABLE)), x, base=base)
        end
    elseif isnothing(groupby)
        readablestring(Readable(sepwith=sepwith, groupby=intgroup(READABLE)), x, base=base)
    else
        readablestring(Readable(groupby=groupby, sepwith=sepwith), x, base=base)
    end
end

readablestring(r::Readable, x::T, base::Int=10) where {T<:Signed} =
    readable_int(r, x, base)

readablestring(x::T, base::Int=10) where {T<:Signed} =
    readablestring(READABLE, x, base)

function readablestring(r::Readable, x::T, base::Int=10) where {T<:AbstractFloat}
    str = string(x)
    return readablestring(r, str, base)
end

readablestring(x::T, base::Int=10) where {T<:AbstractFloat} =
    readablestring(READABLE, x, base)

function readable(io::IO, r::Readable, x::T, base::Int=10) where {T<:Signed}
    str = readablestring(r, x, base)
    print(io, str)
end

readable(io::IO, x::T, base::Int=10) where {T<:Signed} =
    readable(io, READABLE, x, base)

readable(r::Readable, x::T, base::Int=10) where {T<:Signed} =
    readable(Base.stdout, r, x, base)

readable(x::T, base::Int=10) where {F, T<:Signed} =
    readable(Base.stdout, READABLE, x, base)

function readable(io::IO, r::Readable, x::T, base::Int=10) where {T<:AbstractFloat}
    str = readablestring(r, x, base)
    print(io, str)
end

readable(io::IO, x::T, base::Int=10) where {T<:AbstractFloat} =
    readable(io, READABLE, x, base)

readable(r::Readable, x::T, base::Int=10) where {T<:AbstractFloat} =
    readable(Base.stdout, r, x, base)

readable(x::T, base::Int=10) where {F, T<:AbstractFloat} =
    readable(Base.stdout, READABLE, x, base)

function readablestring(r::Readable, x::T, base::Int=10) where {T<:Real}
    str = string(x)
    return readablestring(r, str, base)
end

readablestring(x::T, base::Int=10) where {T<:Real} =
    readablestring(READABLE, x, base)

function readable(io::IO, r::Readable, x::T, base::Int=10) where {T<:Real}
    str = readablestring(r, x, base)
    print(io, str)
end

readable(io::IO, x::T, base::Int=10) where {T<:Real} =
    readable(io, READABLE, x, base)

readable(r::Readable, x::T, base::Int=10) where {T<:Real} =
    readable(Base.stdout, r, x, base)

readable(x::T, base::Int=10) where {F, T<:Real} =
    readable(Base.stdout, READABLE, x, base)

function readablestring(r::Readable, x::T, base::Int=10) where {F, T<:Complex{F}}
    re = real(x)
    im = imag(x)
    sgn = signbit(im) ? " - " : " + "
    im = abs(im)
    re_str = readable(r, string(re), base)
    im_str = readable(r, string(im), base)
    return string(re_str, sgn, im_str, IMAG_UNIT_STR[1])
end

readablestring(x::T, base::Int=10) where {F, T<:Complex{F}} =
    readablestring(READABLE, x, base)

function readable(io::IO, r::Readable, x::T, base::Int=10) where {F, T<:Complex{F}}
    str = readablestring(r, x, base)
    print(io, str)
end

readable(io::IO, x::T, base::Int=10) where {F, T<:Complex{F}} =
    readable(io, READABLE, x, base)

readable(r::Readable, x::T, base::Int=10) where {F, T<:Complex{F}} =
    readable(Base.stdout, r, x, base)

readable(x::T, base::Int=10) where {F, T<:Complex{F}} =
    readable(Base.stdout, READABLE, x, base)

function readablestring(r::Readable, x::T, base::Int=10) where {T<:Number}
     if hasmethod(real, (T,))
        re = real(x)
        if hasmethod(imag, (T,))
            im = imag(x)
            if isa(im, Real)
                readablestring(r, re, im, IMAG_UNIT_STR[1], base)
            else
                throw(DomainError("$T is not supported"))
            end
        elseif hasmethod(dual, (T,))
            du = dual(x)
            if isa(im, Real)
                readablestring(r, re, du, DUAL_UNIT_STR[1], base)
            else
                throw(DomainError("$T is not supported"))
            end
        else
            throw(DomainError("$T is not supported"))
        end
    else
       throw(DomainError("$T is not supported"))
    end
end

readablestring(x::T, base::Int=10) where {T<:Number} =
    readablestring(READABLE, x, base)

function readablestring(r::Readable, x::T, y::T, unitstr::String, base::Int=10)  where {T<:Real}
    sgn = signbit(y) ? " - " : " + "
    y = abs(y)
    xstr = readablestring(r, x, base)
    ystr = readablestring(r, y, base)
    return string(xstr, sgn, ystr, unitstr)
end

readablestring(x::T, y::T, unitstr::String, base::Int=10) where {T<:Real} =
    readablestring(READABLE, x, y, unitstr, base)

function readable(io::IO, r::Readable, x::T, base::Int=10) where {T<:Number}
    str = readablestring(r, x, base)
    print(io, str)
end

readable(io::IO, x::T, base::Int=10) where {T<:Number} =
    readable(io, READABLE, x, base)

readable(r::Readable, x::T, base::Int=10) where {T<:Number} =
    readable(Base.stdout, r, x, base)

readable(x::T, base::Int=10) where {T<:Number} =
    readable(Base.stdout, READABLE, x, base)

splitstr(str::AbstractString, at::Union{String, Char}) = String.(split(str, at))
stripstr(str::AbstractString) = String(strip(str))

function readablestring(r::Readable, str::String, base::Int=10)
    if !occursin(READABLE.decimalpoint, str)
       readable_int(r, BigInt(str), base)
    else
       ipart, fpart = splitstr(str, READABLE.decimalpoint)
       if occursin("e", fpart)
          fpart, epart = splitstr(fpart, "e")
          epart = (epart[1] !== '-' && epart[1] !== '+') ? string("e+", epart) : string("e", epart)
       else
          epart = ""
       end
       ripart = readable_int(r, BigInt(ipart), base)
       rfpart = readable_frac(r, BigInt(fpart), base)
       string(ripart, r.decimalpoint, rfpart, epart)
    end
end

readablestring(x::String, base::Int=10) =
    readablestring(READABLE, x, base)

function readable_int(r::Readable, x::I, base::Int=10) where {I<:Signed}
    numsign = signbit(x) ? "-" : ""
    str = string(abs(x), base=base)
    ndigs = length(str)
    ngroups, firstgroup = divrem(ndigs, r.intgroupsize)

    ngroups == 0 && return str

    idx = firstgroup
    if idx > 0
        res = string(str[1:idx], r.intsepchar)
    else
        res = ""
    end

    while ngroups > 1
        res = string(res, str[idx+1:idx+r.intgroupsize], r.intsepchar)
        idx += r.intgroupsize
        ngroups -= 1
    end
    res = string(res, str[idx+1:idx+r.intgroupsize])

    return string(numsign, baseprefix(base), res)
end

readable_int(x::I, base::Int=10) where {I<:Signed} = readable(READABLE, x, base)

function readable_frac(r::Readable, x::I, base::Int=10) where {I<:Signed}
    signbit(x) && throw(ErrorException("negative fractional parts ($x) are not allowed"))
    str = string(abs(x), base=base)
    ndigs = length(str)
    ngroups, lastgroup = divrem(ndigs, r.fracgroupsize)

    ngroups == 0 && return str

    idx = 0
    res = ""

    while ngroups > 1
        res = string(res, str[idx+1:idx+r.fracgroupsize], r.fracsepchar)
        idx += r.fracgroupsize
        ngroups -= 1
    end
    if lastgroup == 0
        res = string(res, str[idx+1:end])
    else
        res = string(res, str[idx+1:idx+r.fracgroupsize], r.fracsepchar, str[idx+r.fracgroupsize+1:end])
    end

    return res
end

readable_frac(x::I, base::Int=10) where {I<:Signed} = readable_frac(READABLE, x, base)

function Base.BigInt(str::AbstractString)
   s = stripstr(str)
   nchars = length(s)
   prec = ceil(Int, log2(10) * nchars) + 16
   holdprec = precision(BigFloat)
   setprecision(BigFloat, prec)
   res = BigInt(BigFloat(s))
   setprecision(BigFloat, holdprec)
   return res
end

Base.BigInt(str::SubString) = BigInt(String(str))

end # Readables


text:
function ispositive(cumulant::Array{ComplexF64,4},nx::Int,ny::Int,Λ::Int)
    twopoint = reshape(cumulant,(2*ny-1)*(nx-Λ),(2*ny-1)*(nx-Λ))
    D = eigvals(twopoint)
    !any(x->x<0.0,D)
end

function positivity!(cumulant::Array{ComplexF64,4},nx::Int,ny::Int,Λ::Int)
    twopoint = reshape(cumulant,(2*ny-1)*(nx-Λ),(2*ny-1)*(nx-Λ))
    D,V = eigen(twopoint)
    @info "Removing following eignvalues from second cumulant ", D[isless.(D,0)]
    Dpos = max.(D,0.0)
    twopoint = V*diagm(Dpos)*inv(V)
    # optimise further:
    # mul!(twopoint,V,lmul!(diagm(Dpos),inv(V)))
    cumulant = reshape(twopoint,2*ny-1,nx-Λ,2*ny-1,nx-Λ)
end

function rankis(cumulant::Array{ComplexF64,4},nx::Int,ny::Int,Λ::Int)
    twopoint = reshape(cumulant,(2*ny-1)*(nx-Λ),(2*ny-1)*(nx-Λ))
    D = eigvals(twopoint)
    return ((2*ny-1)*(nx-Λ) - length(D[D .< 1e-6]),D)
end


text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

# World Bank Population database

const WORLD_POP_WB_URL = "http://api.worldbank.org/v2/en/indicator/SP.POP.TOTL?downloadformat=csv"
const WORLD_POP_YEAR_0 = 2000

"""
Total population of each country, extracted from https://data.worldbank.org/
"""
struct WorldPopulationWB <: AbstractDataSource
    function WorldPopulationWB()
        new()
    end
end
WorldPopulationWB(::Nothing; kwargs...) = WorldPopulationWB()

DATA_SOURCES[:world_population] = WorldPopulationWB

function import_data(
    source::WorldPopulationWB;
    update_period = Year(1),
    url = WORLD_POP_WB_URL,
    kwargs...,
)
    @debug "Importing population data from World Bank."
    curyear = year(today())
    years = [string(y) for y ∈ WORLD_POP_YEAR_0:curyear]
    keycols = ["Country Name", "Country Code"]
    col_select = [keycols; years]
    metacols = ["Country Code", "Region", "IncomeGroup"]

    file = import_data(DownloadSource(url); update_period = update_period, kwargs...)
    global r = ZipFile.Reader(string(file))
    zipfiles = r.files
    file = zipfiles[findfirst(f -> f.name[1:10] == "API_SP.POP", zipfiles)]
    metafile = zipfiles[findfirst(f -> f.name[1:16] == "Metadata_Country", zipfiles)]

    data = csv_read(file; copycols = true, header = 5, select = col_select)
    metadata = csv_read(metafile; copycols = true, select = metacols)

    close(file)
    close(metafile)
    global r = nothing

    @debug "Joining data and metadata." _debuginfo(data) _debuginfo(metadata)
    data = leftjoin(data, metadata; on = "Country Code")
    data = select!(data, "Country Name", metacols..., intersect(years, names(data))...)

    rename!(
        data,
        "Country Name" => :country,
        "Country Code" => :country_code,
        "Region" => :region,
        "IncomeGroup" => :income_group,
    )
    data[!, :country] = correct_countries!(data.country)
    sort!(data, :country)

    ncols = ncol(data)
    for nr ∈ ncols:-1:1
        if !all(ismissing.(data[!, nr]))
            if nr < ncols
                select!(data, Not((nr+1):ncols))
            end
            break
        end
    end
    years = intersect(Symbol.(years), Symbol.(names(data)))
    fyear = years[1]
    lyear = years[end]

    delete!(data, data[!, fyear] .=== missing)
    myear = years[findlast(y -> all(data[!, y] .!== missing), years)]

    fyear_int = parse(Int, string(fyear))
    lyear_int = parse(Int, string(lyear))
    myear_int = parse(Int, string(myear))

    n = (lyear_int - fyear_int)
    avg_pop_factor = (data[!, lyear] ./ data[!, fyear]) .^ (1 / n)
    navg_factor = Float64[
        factor === missing ? (data[i, myear] / data[i, fyear])^(1 / n) : factor
        for (i, factor) ∈ enumerate(avg_pop_factor)
    ]
    insertcols!(data, ncol(data) + 1, :avg_pop_factor => navg_factor)

    n = (curyear - lyear_int)
    pop_est = (data[!, lyear] .* (data.avg_pop_factor .^ n))
    pop_est = [
        ismissing(val) ? data[i, myear] * (data[i, :avg_pop_factor] .^ n) : val
        for (i, val) ∈ enumerate(pop_est)
    ]
    pop_est = pop_est .|> round .|> Int
    insertcols!(data, ncol(data) + 1, :estimated_population => pop_est)

    select(
        data,
        :country,
        :country_code,
        :region,
        :income_group,
        :avg_pop_factor,
        :estimated_population,
    )
end

# World Bank GDP per capita database

const WORLD_GDP_PER_CAPITA_WB_URL = "http://api.worldbank.org/v2/en/indicator/NY.GDP.PCAP.CD?downloadformat=csv"

const WORLD_GDP_PER_CAPITA_YEAR_0 = 2000

"""
GDP per capita of each country, extracted from https://data.worldbank.org/
"""
struct WorldGDPPerCapitaWB <: AbstractDataSource
    function WorldGDPPerCapitaWB()
        new()
    end
end
WorldGDPPerCapitaWB(::Nothing; kwargs...) = WorldGDPPerCapitaWB()

DATA_SOURCES[:world_gdp_per_capita] = WorldGDPPerCapitaWB

function import_data(
    source::WorldGDPPerCapitaWB;
    update_period = Year(1),
    url = WORLD_GDP_PER_CAPITA_WB_URL,
    kwargs...,
)
    @debug "Importing GDP per capita data from World Bank."
    curyear = year(today())
    years = [string(y) for y ∈ WORLD_GDP_PER_CAPITA_YEAR_0:curyear]
    keycols = ["Country Name", "Country Code"]
    col_select = [keycols; years]
    metacols = ["Country Code", "Region", "IncomeGroup"]

    file = import_data(DownloadSource(url); update_period = update_period, kwargs...)
    global r = ZipFile.Reader(string(file))
    zipfiles = r.files
    file = zipfiles[findfirst(f -> f.name[1:15] == "API_NY.GDP.PCAP", zipfiles)]
    metafile = zipfiles[findfirst(f -> f.name[1:16] == "Metadata_Country", zipfiles)]

    data = csv_read(file; copycols = true, header = 5, select = col_select)
    metadata = csv_read(metafile; copycols = true, select = metacols)

    close(file)
    close(metafile)
    global r = nothing

    @debug "Joining data and metadata." _debuginfo(data) _debuginfo(metadata)
    data = leftjoin(data, metadata; on = "Country Code")
    data = select!(data, "Country Name", metacols..., intersect(years, names(data))...)

    rename!(
        data,
        "Country Name" => :country,
        "Country Code" => :country_code,
        "Region" => :region,
        "IncomeGroup" => :income_group,
    )
    data.country = correct_countries!(data.country)
    sort!(data, :country)

    ncols = ncol(data)
    for nr ∈ ncols:-1:1
        if !all(ismissing.(data[!, nr]))
            if nr < ncols
                select!(data, Not((nr+1):ncols))
            end
            break
        end
    end

    years = intersect(Symbol.(years), Symbol.(names(data)))
    fyear = years[1]
    lyear = years[end]

    delete!(data, data[!, fyear] .=== missing)
    myear = years[findlast(y -> all(data[!, y] .!== missing), years)]

    fyear_int = parse(Int, string(fyear))
    lyear_int = parse(Int, string(lyear))
    myear_int = parse(Int, string(myear))

    n = (lyear_int - fyear_int)
    avg_gdp_pc_factor = (data[!, lyear] ./ data[!, fyear]) .^ (1 / n)
    navg_factor = Float64[
        factor === missing ? (data[i, myear] / data[i, fyear])^(1 / n) : factor
        for (i, factor) ∈ enumerate(avg_gdp_pc_factor)
    ]
    insertcols!(data, ncol(data) + 1, :avg_gdp_pc_factor => navg_factor)

    n = (curyear - lyear_int)
    gdp_pc_est = data[!, lyear]
    gdp_pc_est =
        [ismissing(val) ? data[i, myear] : val for (i, val) ∈ enumerate(gdp_pc_est)]
    gdp_pc_est = gdp_pc_est .|> round .|> Int
    insertcols!(data, ncol(data) + 1, :gdp_per_capita => gdp_pc_est)

    select!(data, :country, :country_code, :avg_gdp_pc_factor, :gdp_per_capita)
end


text:
export VBasedPolicy

using MacroTools: @forward

"""
    VBasedPolicy(;learner, mapping, explorer=GreedyExplorer())

# Key words & Fields

- `learner`::[`AbstractLearner`](@ref), learn how to estimate state values.
- `mapping`, a customized function `(obs, learner) -> action_values`
- `explorer`::[`AbstractExplorer`](@ref), decide which action to take based on action values.
"""
Base.@kwdef struct VBasedPolicy{L<:AbstractLearner,M,E<:AbstractExplorer} <: AbstractPolicy
    learner::L
    mapping::M
    explorer::E = GreedyExplorer()
end

(p::VBasedPolicy)(obs) = p(obs, ActionStyle(obs))

(p::VBasedPolicy)(obs, ::MinimalActionSet) = p.mapping(obs, p.learner) |> p.explorer

function (p::VBasedPolicy)(obs, ::FullActionSet)
    action_values = p.mapping(obs, p.learner)
    p.explorer(action_values, get_legal_actions_mask(obs))
end

RLBase.get_prob(p::VBasedPolicy, obs, action::Integer) =
    get_prob(p, obs, ActionStyle(obs), action)

RLBase.get_prob(p::VBasedPolicy, obs, ::MinimalActionSet) =
    get_prob(p.explorer, p.mapping(obs, p.learner))
RLBase.get_prob(p::VBasedPolicy, obs, ::MinimalActionSet, action) =
    get_prob(p.explorer, p.mapping(obs, p.learner), action)

function RLBase.get_prob(p::VBasedPolicy, obs, ::FullActionSet)
    action_values = p.mapping(obs, p.learner)
    get_prob(p.explorer, action_values, get_legal_actions_mask(obs))
end

function RLBase.get_prob(p::VBasedPolicy, obs, ::FullActionSet, action)
    action_values = p.mapping(obs, p.learner)
    get_prob(p.explorer, action_values, get_legal_actions_mask(obs), action)
end

@forward VBasedPolicy.learner RLBase.get_priority, RLBase.update!


text:
module Atmos

export AtmosModel, AtmosAcousticLinearModel, AtmosAcousticGravityLinearModel

using CLIMAParameters
using CLIMAParameters.Planet: grav, cp_d
using CLIMAParameters.Atmos.SubgridScale: C_smag
using DocStringExtensions
using LinearAlgebra, StaticArrays
using ..ConfigTypes
using ..Orientations
import ..Orientations:
    vertical_unit_vector,
    altitude,
    latitude,
    longitude,
    projection_normal,
    gravitational_potential,
    ∇gravitational_potential,
    projection_tangential

using ..VariableTemplates
using ..Thermodynamics
using ..TemperatureProfiles

using ..TurbulenceClosures
import ..TurbulenceClosures: turbulence_tensors
using ..TurbulenceConvection

import ..Thermodynamics: internal_energy
using ..MPIStateArrays: MPIStateArray
using ..Mesh.Grids:
    VerticalDirection,
    HorizontalDirection,
    min_node_distance,
    EveryDirection,
    Direction

using ClimateMachine.BalanceLaws
using ClimateMachine.Problems

import ClimateMachine.BalanceLaws:
    vars_state,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_state!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    transform_post_gradient_laplacian!,
    init_state_auxiliary!,
    init_state_prognostic!,
    update_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ClimateMachine.DGMethods:
    LocalGeometry, lengthscale, resolutionmetric, DGModel

import ..DGMethods.NumericalFluxes:
    boundary_state!,
    boundary_flux_second_order!,
    normal_boundary_flux_second_order!,
    NumericalFluxFirstOrder,
    NumericalFluxGradient,
    NumericalFluxSecondOrder,
    CentralNumericalFluxHigherOrder,
    CentralNumericalFluxDivergence,
    CentralNumericalFluxFirstOrder,
    numerical_flux_first_order!,
    NumericalFluxFirstOrder
using ..DGMethods.NumericalFluxes:
    RoeNumericalFlux, HLLCNumericalFlux, RusanovNumericalFlux

import ..Courant: advective_courant, nondiffusive_courant, diffusive_courant

"""
    AtmosModel <: BalanceLaw

A `BalanceLaw` for atmosphere modeling. Users may over-ride prescribed
default values for each field.

# Usage

    AtmosModel(
        param_set,
        problem,
        orientation,
        ref_state,
        turbulence,
        hyperdiffusion,
        spongelayer,
        moisture,
        radiation,
        source,
        tracers,
        data_config,
    )

# Fields
$(DocStringExtensions.FIELDS)
"""
struct AtmosModel{FT, PS, PR, O, RS, T, TC, HD, VS, M, P, R, S, TR, DC} <:
       BalanceLaw
    "Parameter Set (type to dispatch on, e.g., planet parameters. See CLIMAParameters.jl package)"
    param_set::PS
    "Problem (initial and boundary conditions)"
    problem::PR
    "An orientation model"
    orientation::O
    "Reference State (For initial conditions, or for linearisation when using implicit solvers)"
    ref_state::RS
    "Turbulence Closure (Equations for dynamics of under-resolved turbulent flows)"
    turbulence::T
    "Turbulence Convection Closure (e.g., EDMF)"
    turbconv::TC
    "Hyperdiffusion Model (Equations for dynamics of high-order spatial wave attenuation)"
    hyperdiffusion::HD
    "Viscous sponge layers"
    viscoussponge::VS
    "Moisture Model (Equations for dynamics of moist variables)"
    moisture::M
    "Precipitation Model (Equations for dynamics of precipitating species)"
    precipitation::P
    "Radiation Model (Equations for radiative fluxes)"
    radiation::R
    "Source Terms (Problem specific source terms)"
    source::S
    "Tracer Terms (Equations for dynamics of active and passive tracers)"
    tracers::TR
    "Data Configuration (Helper field for experiment configuration)"
    data_config::DC
end

"""
    AtmosModel{FT}()

Constructor for `AtmosModel` (where `AtmosModel <: BalanceLaw`) for LES
and single stack configurations.
"""
function AtmosModel{FT}(
    ::Union{Type{AtmosLESConfigType}, Type{SingleStackConfigType}},
    param_set::AbstractParameterSet;
    init_state_prognostic::ISP = nothing,
    problem::PR = AtmosProblem(init_state_prognostic = init_state_prognostic),
    orientation::O = FlatOrientation(),
    ref_state::RS = HydrostaticState(DecayingTemperatureProfile{FT}(param_set),),
    turbulence::T = SmagorinskyLilly{FT}(0.21),
    turbconv::TC = NoTurbConv(),
    hyperdiffusion::HD = NoHyperDiffusion(),
    viscoussponge::VS = NoViscousSponge(),
    moisture::M = EquilMoist{FT}(),
    precipitation::P = NoPrecipitation(),
    radiation::R = NoRadiation(),
    source::S = (
        Gravity(),
        Coriolis(),
        GeostrophicForcing{FT}(7.62e-5, 0, 0),
        turbconv_sources(turbconv)...,
    ),
    tracers::TR = NoTracers(),
    data_config::DC = nothing,
) where {FT <: AbstractFloat, ISP, PR, O, RS, T, TC, HD, VS, M, P, R, S, TR, DC}

    atmos = (
        param_set,
        problem,
        orientation,
        ref_state,
        turbulence,
        turbconv,
        hyperdiffusion,
        viscoussponge,
        moisture,
        precipitation,
        radiation,
        source,
        tracers,
        data_config,
    )

    return AtmosModel{FT, typeof.(atmos)...}(atmos...)
end

"""
    AtmosModel{FT}()

Constructor for `AtmosModel` (where `AtmosModel <: BalanceLaw`) for GCM
configurations.
"""
function AtmosModel{FT}(
    ::Type{AtmosGCMConfigType},
    param_set::AbstractParameterSet;
    init_state_prognostic::ISP = nothing,
    problem::PR = AtmosProblem(init_state_prognostic = init_state_prognostic),
    orientation::O = SphericalOrientation(),
    ref_state::RS = HydrostaticState(DecayingTemperatureProfile{FT}(param_set),),
    turbulence::T = SmagorinskyLilly{FT}(C_smag(param_set)),
    turbconv::TC = NoTurbConv(),
    hyperdiffusion::HD = NoHyperDiffusion(),
    viscoussponge::VS = NoViscousSponge(),
    moisture::M = EquilMoist{FT}(),
    precipitation::P = NoPrecipitation(),
    radiation::R = NoRadiation(),
    source::S = (Gravity(), Coriolis(), turbconv_sources(turbconv)...),
    tracers::TR = NoTracers(),
    data_config::DC = nothing,
) where {FT <: AbstractFloat, ISP, PR, O, RS, T, TC, HD, VS, M, P, R, S, TR, DC}

    atmos = (
        param_set,
        problem,
        orientation,
        ref_state,
        turbulence,
        turbconv,
        hyperdiffusion,
        viscoussponge,
        moisture,
        precipitation,
        radiation,
        source,
        tracers,
        data_config,
    )

    return AtmosModel{FT, typeof.(atmos)...}(atmos...)
end

"""
    vars_state(m::AtmosModel, ::Prognostic, FT)

Conserved state variables (prognostic variables).
"""
function vars_state(m::AtmosModel, st::Prognostic, FT)
    @vars begin
        ρ::FT
        ρu::SVector{3, FT}
        ρe::FT
        turbulence::vars_state(m.turbulence, st, FT)
        turbconv::vars_state(m.turbconv, st, FT)
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
        moisture::vars_state(m.moisture, st, FT)
        radiation::vars_state(m.radiation, st, FT)
        tracers::vars_state(m.tracers, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::Gradient, FT)

Pre-transform gradient variables.
"""
function vars_state(m::AtmosModel, st::Gradient, FT)
    @vars begin
        u::SVector{3, FT}
        h_tot::FT
        turbulence::vars_state(m.turbulence, st, FT)
        turbconv::vars_state(m.turbconv, st, FT)
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
        moisture::vars_state(m.moisture, st, FT)
        tracers::vars_state(m.tracers, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::GradientFlux, FT)

Post-transform gradient variables.
"""
function vars_state(m::AtmosModel, st::GradientFlux, FT)
    @vars begin
        ∇h_tot::SVector{3, FT}
        turbulence::vars_state(m.turbulence, st, FT)
        turbconv::vars_state(m.turbconv, st, FT)
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
        moisture::vars_state(m.moisture, st, FT)
        tracers::vars_state(m.tracers, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::GradientLaplacian, FT)

Pre-transform hyperdiffusive variables.
"""
function vars_state(m::AtmosModel, st::GradientLaplacian, FT)
    @vars begin
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::Hyperdiffusive, FT)

Post-transform hyperdiffusive variables.
"""
function vars_state(m::AtmosModel, st::Hyperdiffusive, FT)
    @vars begin
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::Auxiliary, FT)

Auxiliary variables, such as vertical (stack) integrals, coordinates,
orientation information, reference states, subcomponent auxiliary vars,
debug variables.
"""
function vars_state(m::AtmosModel, st::Auxiliary, FT)
    @vars begin
        ∫dz::vars_state(m, UpwardIntegrals(), FT)
        ∫dnz::vars_state(m, DownwardIntegrals(), FT)
        coord::SVector{3, FT}
        orientation::vars_state(m.orientation, st, FT)
        ref_state::vars_state(m.ref_state, st, FT)
        turbulence::vars_state(m.turbulence, st, FT)
        turbconv::vars_state(m.turbconv, st, FT)
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
        moisture::vars_state(m.moisture, st, FT)
        tracers::vars_state(m.tracers, st, FT)
        radiation::vars_state(m.radiation, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::UpwardIntegrals, FT)
"""
function vars_state(m::AtmosModel, st::UpwardIntegrals, FT)
    @vars begin
        radiation::vars_state(m.radiation, st, FT)
        turbconv::vars_state(m.turbconv, st, FT)
    end
end

"""
    vars_state(m::AtmosModel, ::DownwardIntegrals, FT)
"""
function vars_state(m::AtmosModel, st::DownwardIntegrals, FT)
    @vars begin
        radiation::vars_state(m.radiation, st, FT)
    end
end

####
#### Forward orientation methods
####
projection_normal(bl, aux, u⃗) =
    projection_normal(bl.orientation, bl.param_set, aux, u⃗)
projection_tangential(bl, aux, u⃗) =
    projection_tangential(bl.orientation, bl.param_set, aux, u⃗)
latitude(bl, aux) = latitude(bl.orientation, aux)
longitude(bl, aux) = longitude(bl.orientation, aux)
altitude(bl, aux) = altitude(bl.orientation, bl.param_set, aux)
vertical_unit_vector(bl, aux) =
    vertical_unit_vector(bl.orientation, bl.param_set, aux)
gravitational_potential(bl, aux) = gravitational_potential(bl.orientation, aux)
∇gravitational_potential(bl, aux) =
    ∇gravitational_potential(bl.orientation, aux)

turbulence_tensors(atmos::AtmosModel, args...) =
    turbulence_tensors(atmos.turbulence, atmos, args...)

include("problem.jl")
include("ref_state.jl")
include("moisture.jl")
include("thermo_states.jl")
include("precipitation.jl")
include("radiation.jl")
include("source.jl")
include("tracers.jl")
include("linear.jl")
include("courant.jl")
include("filters.jl")

"""
    flux_first_order!(
        m::AtmosModel,
        flux::Grad,
        state::Vars,
        aux::Vars,
        t::Real
    )

Computes and assembles non-diffusive fluxes in the model
equations.
"""
@inline function flux_first_order!(
    m::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρ = state.ρ
    ρinv = 1 / ρ
    ρu = state.ρu
    u = ρinv * ρu

    # advective terms
    flux.ρ = ρ * u
    flux.ρu = ρ * u .* u'
    flux.ρe = u * state.ρe

    # pressure terms
    ts = recover_thermo_state(m, state, aux)
    p = air_pressure(ts)
    if m.ref_state isa HydrostaticState
        flux.ρu += (p - aux.ref_state.p) * I
    else
        flux.ρu += p * I
    end
    flux.ρe += u * p
    flux_radiation!(m.radiation, m, flux, state, aux, t)
    flux_moisture!(m.moisture, m, flux, state, aux, t)
    flux_tracers!(m.tracers, m, flux, state, aux, t)
    flux_first_order!(m.turbconv, m, flux, state, aux, t)
end

function compute_gradient_argument!(
    atmos::AtmosModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.u = ρinv * state.ρu
    ts = recover_thermo_state(atmos, state, aux)
    e_tot = state.ρe * (1 / state.ρ)
    transform.h_tot = total_specific_enthalpy(ts, e_tot)

    compute_gradient_argument!(atmos.moisture, transform, state, aux, t)
    compute_gradient_argument!(atmos.turbulence, transform, state, aux, t)
    compute_gradient_argument!(
        atmos.hyperdiffusion,
        atmos,
        transform,
        state,
        aux,
        t,
    )
    compute_gradient_argument!(atmos.tracers, transform, state, aux, t)
    compute_gradient_argument!(atmos.turbconv, atmos, transform, state, aux, t)
end

function compute_gradient_flux!(
    atmos::AtmosModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.∇h_tot = ∇transform.h_tot

    # diffusion terms required for SGS turbulence computations
    compute_gradient_flux!(
        atmos.turbulence,
        atmos.orientation,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )
    # diffusivity of moisture components
    compute_gradient_flux!(atmos.moisture, diffusive, ∇transform, state, aux, t)
    compute_gradient_flux!(atmos.tracers, diffusive, ∇transform, state, aux, t)
    compute_gradient_flux!(
        atmos.turbconv,
        atmos,
        diffusive,
        ∇transform,
        state,
        aux,
        t,
    )
end

function transform_post_gradient_laplacian!(
    atmos::AtmosModel,
    hyperdiffusive::Vars,
    hypertransform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform_post_gradient_laplacian!(
        atmos.hyperdiffusion,
        atmos,
        hyperdiffusive,
        hypertransform,
        state,
        aux,
        t,
    )
end

"""
    flux_second_order!(
        atmos::AtmosModel,
        flux::Grad,
        state::Vars,
        diffusive::Vars,
        hyperdiffusive::Vars,
        aux::Vars,
        t::Real
    )
Diffusive fluxes in AtmosModel. Viscosity, diffusivity are calculated
in the turbulence subcomponent and accessed within the diffusive flux
function. Contributions from subcomponents are then assembled (pointwise).
"""
@inline function flux_second_order!(
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    ν, D_t, τ = sponge_viscosity_modifier(atmos, atmos.viscoussponge, ν, D_t, τ, aux)
    d_h_tot = -D_t .* diffusive.∇h_tot
    flux_second_order!(atmos, flux, state, τ, d_h_tot)
    flux_second_order!(atmos.moisture, flux, state, diffusive, aux, t, D_t)
    flux_second_order!(
        atmos.hyperdiffusion,
        flux,
        state,
        diffusive,
        hyperdiffusive,
        aux,
        t,
    )
    flux_second_order!(atmos.tracers, flux, state, diffusive, aux, t, D_t)
    flux_second_order!(atmos.turbconv, atmos, flux, state, diffusive, aux, t)
end

#TODO: Consider whether to not pass ρ and ρu (not state), foc BCs reasons
@inline function flux_second_order!(
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    τ,
    d_h_tot,
)
    flux.ρu += τ * state.ρ
    flux.ρe += τ * state.ρu
    flux.ρe += d_h_tot * state.ρ
end

@inline function wavespeed(
    m::AtmosModel,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρinv = 1 / state.ρ
    u = ρinv * state.ρu
    uN = abs(dot(nM, u))
    ts = recover_thermo_state(m, state, aux)
    ss = soundspeed_air(ts)

    FT = typeof(state.ρ)
    ws = fill(uN + ss, MVector{number_states(m, Prognostic()), FT})
    vars_ws = Vars{vars_state(m, Prognostic(), FT)}(ws)

    wavespeed_tracers!(m.tracers, vars_ws, nM, state, aux, t)

    return ws
end

function update_auxiliary_state!(
    dg::DGModel,
    m::AtmosModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    FT = eltype(Q)
    state_auxiliary = dg.state_auxiliary

    if number_states(m, UpwardIntegrals()) > 0
        indefinite_stack_integral!(dg, m, Q, state_auxiliary, t, elems)
        reverse_indefinite_stack_integral!(dg, m, Q, state_auxiliary, t, elems)
    end

    update_auxiliary_state!(nodal_update_auxiliary_state!, dg, m, Q, t, elems)

    # TODO: Remove this hook. This hook was added for implementing
    # the first draft of EDMF, and should be removed so that we can
    # rely on a single vertical element traversal. This hook allows
    # us to compute globally vertical quantities specific to EDMF
    # until we're able to remove them or somehow incorporate them
    # into a higher level hierarchy.
    update_auxiliary_state!(dg, m.turbconv, m, Q, t, elems)

    return true
end

function nodal_update_auxiliary_state!(
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    atmos_nodal_update_auxiliary_state!(m.moisture, m, state, aux, t)
    atmos_nodal_update_auxiliary_state!(m.radiation, m, state, aux, t)
    atmos_nodal_update_auxiliary_state!(m.tracers, m, state, aux, t)
    turbulence_nodal_update_auxiliary_state!(m.turbulence, m, state, aux, t)
    turbconv_nodal_update_auxiliary_state!(m.turbconv, m, state, aux, t)
end

function integral_load_auxiliary_state!(
    m::AtmosModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
)
    integral_load_auxiliary_state!(m.radiation, integ, state, aux)
    integral_load_auxiliary_state!(m.turbconv, m, integ, state, aux)
end

function integral_set_auxiliary_state!(m::AtmosModel, aux::Vars, integ::Vars)
    integral_set_auxiliary_state!(m.radiation, aux, integ)
    integral_set_auxiliary_state!(m.turbconv, m, aux, integ)
end

function reverse_integral_load_auxiliary_state!(
    m::AtmosModel,
    integ::Vars,
    state::Vars,
    aux::Vars,
)
    reverse_integral_load_auxiliary_state!(m.radiation, integ, state, aux)
end

function reverse_integral_set_auxiliary_state!(
    m::AtmosModel,
    aux::Vars,
    integ::Vars,
)
    reverse_integral_set_auxiliary_state!(m.radiation, aux, integ)
end

function atmos_nodal_init_state_auxiliary!(
    m::AtmosModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    aux.coord = geom.coord
    init_aux_turbulence!(m.turbulence, m, aux, geom)
    atmos_init_aux!(m.ref_state, m, aux, tmp, geom)
    init_aux_hyperdiffusion!(m.hyperdiffusion, m, aux, geom)
    atmos_init_aux!(m.tracers, m, aux, geom)
    init_aux_turbconv!(m.turbconv, m, aux, geom)
    m.problem.init_state_auxiliary(m.problem, m, aux, geom)
end

"""
    init_state_auxiliary!(
        m::AtmosModel,
        aux::Vars,
        grid,
        direction
    )

Initialise auxiliary variables for each AtmosModel subcomponent.
Store Cartesian coordinate information in `aux.coord`.
"""
function init_state_auxiliary!(
    m::AtmosModel,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    init_aux!(m, m.orientation, state_auxiliary, grid, direction)

    init_state_auxiliary!(
        m,
        (m, aux, tmp, geom) ->
            atmos_init_ref_state_pressure!(m.ref_state, m, aux, geom),
        state_auxiliary,
        grid,
        direction,
    )

    ∇p = ∇reference_pressure(m.ref_state, state_auxiliary, grid)

    init_state_auxiliary!(
        m,
        atmos_nodal_init_state_auxiliary!,
        state_auxiliary,
        grid,
        direction;
        state_temporary = ∇p,
    )
end

"""
    source!(
        m::AtmosModel,
        source::Vars,
        state::Vars,
        diffusive::Vars,
        aux::Vars,
        t::Real,
        direction::Direction,
    )

Computes (and assembles) source terms `S(Y)` in:
```
∂Y
-- = - ∇ • F + S(Y)
∂t
```
"""
function source!(
    m::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    atmos_source!(m.source, m, source, state, diffusive, aux, t, direction)
end

"""
    init_state_prognostic!(
        m::AtmosModel,
        state::Vars,
        aux::Vars,
        coords,
        t,
        args...,
    )

Initialise state variables. `args...` provides an option to include
configuration data (current use cases include problem constants,
spline-interpolants).
"""
function init_state_prognostic!(
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    coords,
    t,
    args...,
)
    m.problem.init_state_prognostic(
        m.problem,
        m,
        state,
        aux,
        coords,
        t,
        args...,
    )
end

roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
    (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))

function numerical_flux_first_order!(
    numerical_flux::RoeNumericalFlux,
    balance_law::AtmosModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}
    @assert balance_law.moisture isa DryModel

    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        balance_law,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

    FT = eltype(fluxᵀn)
    param_set = balance_law.param_set
    _cv_d::FT = cv_d(param_set)
    _T_0::FT = T_0(param_set)

    Φ = gravitational_potential(balance_law, state_auxiliary⁻)

    ρ⁻ = state_prognostic⁻.ρ
    ρu⁻ = state_prognostic⁻.ρu
    ρe⁻ = state_prognostic⁻.ρe
    ts⁻ = recover_thermo_state(
        balance_law,
        balance_law.moisture,
        state_prognostic⁻,
        state_auxiliary⁻,
    )

    u⁻ = ρu⁻ / ρ⁻
    uᵀn⁻ = u⁻' * normal_vector
    e⁻ = ρe⁻ / ρ⁻
    h⁻ = total_specific_enthalpy(ts⁻, e⁻)
    p⁻ = air_pressure(ts⁻)
    c⁻ = soundspeed_air(ts⁻)

    ρ⁺ = state_prognostic⁺.ρ
    ρu⁺ = state_prognostic⁺.ρu
    ρe⁺ = state_prognostic⁺.ρe

    # TODO: state_auxiliary⁺ is not up-to-date
    # with state_prognostic⁺ on the boundaries
    ts⁺ = recover_thermo_state(
        balance_law,
        balance_law.moisture,
        state_prognostic⁺,
        state_auxiliary⁺,
    )

    u⁺ = ρu⁺ / ρ⁺
    uᵀn⁺ = u⁺' * normal_vector
    e⁺ = ρe⁺ / ρ⁺
    h⁺ = total_specific_enthalpy(ts⁺, e⁺)
    p⁺ = air_pressure(ts⁺)
    c⁺ = soundspeed_air(ts⁺)

    ρ̃ = sqrt(ρ⁻ * ρ⁺)
    ũ = roe_average(ρ⁻, ρ⁺, u⁻, u⁺)
    h̃ = roe_average(ρ⁻, ρ⁺, h⁻, h⁺)
    c̃ = sqrt(roe_average(ρ⁻, ρ⁺, c⁻^2, c⁺^2))

    ũᵀn = ũ' * normal_vector

    Δρ = ρ⁺ - ρ⁻
    Δp = p⁺ - p⁻
    Δu = u⁺ - u⁻
    Δuᵀn = Δu' * normal_vector

    w1 = abs(ũᵀn - c̃) * (Δp - ρ̃ * c̃ * Δuᵀn) / (2 * c̃^2)
    w2 = abs(ũᵀn + c̃) * (Δp + ρ̃ * c̃ * Δuᵀn) / (2 * c̃^2)
    w3 = abs(ũᵀn) * (Δρ - Δp / c̃^2)
    w4 = abs(ũᵀn) * ρ̃

    fluxᵀn.ρ -= (w1 + w2 + w3) / 2
    fluxᵀn.ρu -=
        (
            w1 * (ũ - c̃ * normal_vector) +
            w2 * (ũ + c̃ * normal_vector) +
            w3 * ũ +
            w4 * (Δu - Δuᵀn * normal_vector)
        ) / 2
    fluxᵀn.ρe -=
        (
            w1 * (h̃ - c̃ * ũᵀn) +
            w2 * (h̃ + c̃ * ũᵀn) +
            w3 * (ũ' * ũ / 2 + Φ - _T_0 * _cv_d) +
            w4 * (ũ' * Δu - ũᵀn * Δuᵀn)
        ) / 2

    if !(balance_law.tracers isa NoTracers)
        ρχ⁻ = state_prognostic⁻.tracers.ρχ
        χ⁻ = ρχ⁻ / ρ⁻

        ρχ⁺ = state_prognostic⁺.tracers.ρχ
        χ⁺ = ρχ⁺ / ρ⁺

        χ̃ = roe_average(ρ⁻, ρ⁺, χ⁻, χ⁺)
        Δρχ = ρχ⁺ - ρχ⁻

        wt = abs(ũᵀn) * (Δρχ - χ̃ * Δp / c̃^2)

        fluxᵀn.tracers.ρχ -= ((w1 + w2) * χ̃ + wt) / 2
    end
end

"""
    NumericalFluxFirstOrder()
        ::HLLCNumericalFlux,
        balance_law::AtmosModel,
        fluxᵀn,
        normal_vector,
        state_prognostic⁻,
        state_auxiliary⁻,
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )

An implementation of the numerical flux based on the HLLC method for
the AtmosModel. For more information on this particular implementation,
see Chapter 10.4 in the provided reference below.

## References
    @book{toro2013riemann,
        title={Riemann solvers and numerical methods for fluid dynamics: a practical introduction},
        author={Toro, Eleuterio F},
        year={2013},
        publisher={Springer Science & Business Media}
    }
"""
function numerical_flux_first_order!(
    ::HLLCNumericalFlux,
    balance_law::AtmosModel,
    fluxᵀn::Vars{S},
    normal_vector::SVector,
    state_prognostic⁻::Vars{S},
    state_auxiliary⁻::Vars{A},
    state_prognostic⁺::Vars{S},
    state_auxiliary⁺::Vars{A},
    t,
    direction,
) where {S, A}
    FT = eltype(fluxᵀn)
    num_state_prognostic = number_states(balance_law, Prognostic())
    param_set = balance_law.param_set

    # Extract the first-order fluxes from the AtmosModel (underlying BalanceLaw)
    # and compute normals on the positive + and negative - sides of the
    # interior facets
    flux⁻ = similar(parent(fluxᵀn), Size(3, num_state_prognostic))
    fill!(flux⁻, -zero(FT))
    flux_first_order!(
        balance_law,
        Grad{S}(flux⁻),
        state_prognostic⁻,
        state_auxiliary⁻,
        t,
        direction,
    )
    fluxᵀn⁻ = flux⁻' * normal_vector

    flux⁺ = similar(flux⁻)
    fill!(flux⁺, -zero(FT))
    flux_first_order!(
        balance_law,
        Grad{S}(flux⁺),
        state_prognostic⁺,
        state_auxiliary⁺,
        t,
        direction,
    )
    fluxᵀn⁺ = flux⁺' * normal_vector

    # Extract relevant fields and thermodynamic variables defined on
    # the positive + and negative - sides of the interior facets
    ρ⁻ = state_prognostic⁻.ρ
    ρu⁻ = state_prognostic⁻.ρu
    ρe⁻ = state_prognostic⁻.ρe
    ts⁻ = recover_thermo_state(
        balance_law,
        balance_law.moisture,
        state_prognostic⁻,
        state_auxiliary⁻,
    )

    u⁻ = ρu⁻ / ρ⁻
    c⁻ = soundspeed_air(ts⁻)

    uᵀn⁻ = u⁻' * normal_vector
    e⁻ = ρe⁻ / ρ⁻
    p⁻ = air_pressure(ts⁻)

    ρ⁺ = state_prognostic⁺.ρ
    ρu⁺ = state_prognostic⁺.ρu
    ρe⁺ = state_prognostic⁺.ρe
    ts⁺ = recover_thermo_state(
        balance_law,
        balance_law.moisture,
        state_prognostic⁺,
        state_auxiliary⁺,
    )

    u⁺ = ρu⁺ / ρ⁺
    uᵀn⁺ = u⁺' * normal_vector
    e⁺ = ρe⁺ / ρ⁺
    p⁺ = air_pressure(ts⁺)
    c⁺ = soundspeed_air(ts⁺)

    # Wave speeds estimates S⁻ and S⁺
    S⁻ = min(uᵀn⁻ - c⁻, uᵀn⁺ - c⁺)
    S⁺ = max(uᵀn⁻ + c⁻, uᵀn⁺ + c⁺)

    # Compute the middle wave speed S⁰ in the contact/star region
    S⁰ =
        (p⁺ - p⁻ + ρ⁻ * uᵀn⁻ * (S⁻ - uᵀn⁻) - ρ⁺ * uᵀn⁺ * (S⁺ - uᵀn⁺)) /
        (ρ⁻ * (S⁻ - uᵀn⁻) - ρ⁺ * (S⁺ - uᵀn⁺))

    p⁰ =
        (
            p⁺ +
            p⁻ +
            ρ⁻ * (S⁻ - uᵀn⁻) * (S⁰ - uᵀn⁻) +
            ρ⁺ * (S⁺ - uᵀn⁺) * (S⁰ - uᵀn⁺)
        ) / 2

    # Compute p * D = p * (0, n₁, n₂, n₃, S⁰)
    pD = @MVector zeros(FT, num_state_prognostic)
    if balance_law.ref_state isa HydrostaticState
        # pressure should be continuous but it doesn't hurt to average
        ref_p⁻ = state_auxiliary⁻.ref_state.p
        ref_p⁺ = state_auxiliary⁺.ref_state.p
        ref_p⁰ = (ref_p⁻ + ref_p⁺) / 2

        momentum_p = p⁰ - ref_p⁰
    else
        momentum_p = p⁰
    end

    pD[2] = momentum_p * normal_vector[1]
    pD[3] = momentum_p * normal_vector[2]
    pD[4] = momentum_p * normal_vector[3]
    pD[5] = p⁰ * S⁰

    # Computes both +/- sides of intermediate flux term flux⁰
    flux⁰⁻ =
        (S⁰ * (S⁻ * parent(state_prognostic⁻) - fluxᵀn⁻) + S⁻ * pD) / (S⁻ - S⁰)
    flux⁰⁺ =
        (S⁰ * (S⁺ * parent(state_prognostic⁺) - fluxᵀn⁺) + S⁺ * pD) / (S⁺ - S⁰)

    if 0 <= S⁻
        parent(fluxᵀn) .= fluxᵀn⁻
    elseif S⁻ < 0 <= S⁰
        parent(fluxᵀn) .= flux⁰⁻
    elseif S⁰ < 0 <= S⁺
        parent(fluxᵀn) .= flux⁰⁺
    else # 0 > S⁺
        parent(fluxᵀn) .= fluxᵀn⁺
    end
end

end # module

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
module Deeplearning

using Flux

include("basic.jl")
include("recurrent.jl")

end


text:
@testset "General" begin
  g = smallgraph(:bull)
  qws = QWSearch(Szegedy(g), [1])

  @test check_qwdynamics(qws) == nothing
  @test_throws AssertionError QWSearch(Szegedy(g), Dict{Symbol, Any}(), [1], 0., check=true)
  @test_no_error QWSearch(Szegedy(g), Dict{Symbol, Any}(), [1], 0., check=false)
end

@testset "Continuous search" begin
  @testset "General" begin
    g = smallgraph(:bull)
    qws = QWSearch(CTQW(g), [1])
    @test check_qwdynamics(qws) == nothing
  end

  @testset "Norm preservation" begin
    g = smallgraph(:bull)
    qws = QWSearch(CTQW(g), [1])
    @test norm(execute_single(qws, 10).state) ≈ 1

    qws = QWSearch(CTQWDense(g), [1])
    @test norm(execute_single(qws, 10).state) ≈ 1
  end

  @testset "CTQW and CTQWDense arch" begin
    n = 5
    g = complete_graph(n)
    qws = QWSearch(CTQW(g), [1], 0, 1/n)
    s0 =  QSearchState(qws, real(initial_state(qws)), 0.0)
    @test sum(execute_single(qws, pi*sqrt(n)/2).probability) ≈ 1
    @test evolve(qws,evolve(qws,s0,0.1),0.1) ≈ evolve(qws,s0,0.2)

    qws = QWSearch(CTQWDense(g), [1], 0, 1/n)
    @test sum(execute_single(qws, pi*sqrt(n)/2).probability) ≈ 1

    n = 5
    g = cycle_graph(5)
    qss_adj = QWSearch(CTQW(g, :adjacency), [1], 0, 1/2)
    qss_lap = QWSearch(CTQW(g, :laplacian), [1], 0, 1/2)
    qss_normlap = QWSearch(CTQW(g, :normalized_laplacian), [1], 0, 1)
    @test sum(execute_single(qss_adj, 2.).probability) ≈ sum(execute_single(qss_lap, 2.).probability)
    @test sum(execute_single(qss_adj, 2.).probability) ≈ sum(execute_single(qss_normlap, 2.).probability)
  end
end

@testset "Discrete search" begin
  g = smallgraph(:bull)
  qws = QWSearch(Szegedy(g), [1])

  @testset "General" begin
    @test check_qwdynamics(qws) == nothing
  end

  @testset "Norm preservation" begin
    @test norm(execute_single(qws, 10).state) ≈ 1
  end

  @testset "Evolution" begin
    steps = 10
    res0 = execute(qws,steps)
    res1 = execute(qws,steps,all=true)
    res2 = execute_all(qws,steps)
    res3 = execute_single(qws,steps)
    # some funky constructions
    @test res1[2].state == QSearchState(qws,evolve(qws,res2[1].state),1).state
    @test res1[10].probability == res2[10].probability
    @test res3.state == res1[11].state
    @test res0.state == res1[11].state
    @test length(res1) == steps+1
    # more straightforward tests
    @test execute_all_measured(qws,steps) == execute(qws,steps,all=true,measure=true)
    @test execute_single_measured(qws,steps) == execute(qws,steps,measure=true)
    @test probability(res3) == res1[steps+1].probability
    @test state(res3) == res1[steps+1].state
    @test state(execute_single(qws,1)) == execute_single(qws,1).state

    s = evolve(qws, QuantumWalk.initial_state(qws))
    @test evolve(qws, s) == evolve(qws, evolve(qws, QuantumWalk.initial_state(qws)))
  end

  @testset "Measurement" begin
    s = QSearchState(qws, initial_state(qws),0)
    @test sum(measure(qws,s)) ≈ 1
    @test sum(measure(qws,s,[1,2,3])) + sum(measure(qws,s,[4,5])) ≈ 1
  end
end


text:
Random.seed!(1234)

p   = 2
P1  = randn(p)
P1 *= P1'
mu  = zeros(p)

mvg            = PDSampler.MvGaussianCanon(mu,P1)
gll(x)         = PDSampler.gradloglik(mvg, x)
nextevent(x,v) = PDSampler.nextevent_bps(mvg, x, v)

listfactors = [ Factor( nextevent, gll, i) for i in 1:3]
chain       = chaingraph(listfactors)
compc       = FactorGraphStruct([[1,2],[2,3],[3,4]])
chain2      = FactorGraph(FactorGraphStruct([[1,2],[2,3],[3,4]]), listfactors)
chain3      = FactorGraph([[1,2],[2,3],[3,4]], listfactors)

i = rand(1:3)
x = randn(p)
v = randn(p)

@test   chain.structure.flist == compc.flist &&
        chain.structure.vlist == compc.vlist

@test   chain2.structure.flist == chain.structure.flist &&
        chain2.structure.vlist == chain.structure.vlist

@test   chain3.structure.flist == chain.structure.flist &&
        chain3.structure.vlist == chain.structure.vlist

@test chain.factors[i].gll(x) == chain2.factors[i].gll(x)

Random.seed!(12); a = chain2.factors[i].nextevent(x,v)
Random.seed!(12); b = chain.factors[i].nextevent(x,v)

@test a==b

fgs = FactorGraphStruct( [[1,2,5],
                          [1,3,6],
                          [3,4,5,6],
                          [4,5]] )

@test assocfactors(fgs,3) == [2,3]
@test assocvariables(fgs,3) == [3,4,5,6]
@test sort(linkedfactors(fgs,3)) == [1,2,4]
@test sort(linkedfactors(fgs,4)) == [1,3]
@test chainstruct(5).flist == [[1,2],[2,3],[3,4],[4,5]]

@test assocfactors(chain,  2) == [1,2]
@test linkedfactors(chain, 2) == [1,3]

xtest = randn(p)
vtest = randn(p)

Random.seed!(12); ne1 = chain.factors[3].nextevent(xtest,vtest)
Random.seed!(12); ne2 = nextevent(xtest,vtest)

@test   chain.structure.flist == [[1,2],[2,3],[3,4]] &&
        chain.factors[1].gll(xtest) == gll(xtest) &&
        ne1.tau == ne2.tau

chain = chainstruct(3)
@test assocvariables(chain, 1) == [1,2] && assocvariables(chain, 2) == [2,3]


text:
export AbstractKernel, RBFKernel, kmat, rbf_kernel, kernel_expect

abstract type AbstractKernel end

"""
    RBFKernel(σ, matrix)

RBF Kernel with dense array as kernel matrix.
"""
struct RBFKernel <: AbstractKernel
    sigma::Float64
    matrix::Matrix{Float64}
end

"""
    kmat(k::AbstractKernel) -> Matrix

Returns Kernel Matrix.
"""
function kmat end
kmat(mbf::RBFKernel) = mbf.matrix

"""
    rbf_kernel(basis, σ::Real) -> RBFKernel

Returns RBF Kernel Matrix.
"""
function rbf_kernel(basis, σ::Real)
    dx2 = (basis .- basis').^2
    RBFKernel(σ, exp.(-1/2σ * dx2))
end

"""
    kernel_expect(kernel::AbstractKernel, px::Vector{Float64}, py::Vector{Float64}) -> Float64

Returns the expectation value of kernel on specific distributions.
"""
kernel_expect(kernel::AbstractKernel, px::Vector, py::Vector=px) = px' * kmat(kernel) * py


text:
using SubHunt
using Base.Test

using POMDPs
using POMDPToolbox
using DiscreteValueIteration
using ProfileView

pomdp = SubHuntPOMDP(size=10)

# show(STDOUT, MIME("text/plain"), SubVis(pomdp))

rng = MersenneTwister(5)
# policy = RandomPolicy(pomdp, rng=rng)

solver = ValueIterationSolver()

policy = solve(solver, pomdp)
Profile.clear()
@time solve(solver, pomdp)
# @profile solve(solver, pomdp, verbose=true)
# ProfileView.view()

for (s, a, r, sp) in stepthrough(pomdp, policy, "s,a,r,sp", max_steps=200, rng=rng)
    v = SubVis(pomdp, s=s, a=a, r=r)
    show(STDOUT, MIME("text/plain"), v)
end


text:
bothat_ref(img, dims, r) = bothat_ref(img, strel_box(img, dims; r))
bothat_ref(img, se) = float.(closing(img, se)) - float.(img)

@testset "bothat" begin
    test_ranges = [Bool, 1:10, Float64, Gray{N0f8}, Gray{Float64}]
    for T in test_ranges
        for N in (1, 2, 3)
            sz = ntuple(_ -> 32, N)
            img = rand(T, sz...)

            out = bothat(img)
            @test out == bothat(img, strel_box(img, ntuple(identity, N); r=1))
            @test out == bothat_ref(img, ntuple(identity, N), 1)

            @test bothat(img; dims=(1,), r=2) == bothat_ref(img, (1,), 2)

            se = centered(rand(Bool, ntuple(_ -> 3, N)))
            @test bothat(img, se) == bothat_ref(img, se)
        end
    end

    img = rand(1:5, 7, 7)
    out = similar(img)
    bothat!(out, img, similar(img))
    @test out == bothat(img)

    img = rand(RGB, 7, 7)
    @test_throws ArgumentError("color image is not supported") bothat(img)
end


text:
using Documenter, .JuliaInterface

makedocs(sitename="JuliaInterface Documentation")

max_stars_repo_name: sonolac/JuliaInterface

max_issues_repo_name: sonolac/JuliaInterface

max_forks_repo_name: sonolac/JuliaInterface


text:
using MetaGraphs
using LightGraphs

L = 2
N = L^2

instance = "$(@__DIR__)/../instances/$(N)_001.txt"

ig = ising_graph(instance)
r = fill(2, N)
set_prop!(ig, :rank, r)
dβ = 0.01
β = 1

sgn = -1.
ϵ = 1E-6
D = prod(r) + 1
var_ϵ = 1E-8
sweeps = 4
control = MPSControl(D, var_ϵ, sweeps, β, dβ)

states = all_states(get_prop(ig, :rank))
ϱ = cu(gibbs_tensor(ig))
@test sum(ϱ) ≈ 1

@testset "MPS from gates" begin

    @testset "Exact Gibbs pure state (MPS)" begin
        L = nv(ig)
        rank = get_prop(ig, :rank)

        @info "Generating Gibbs state - |ρ>" L rank β ϵ

        ψ = ones(rank...)

        for σ ∈ states
            for i ∈ 1:L
                h = get_prop(ig, i, :h)

                nbrs = unique_neighbors(ig, i)
                ψ[idx.(σ)...] *= exp(sgn * 0.5 * β * h * σ[i])

                for j ∈ nbrs
                    J = get_prop(ig, i, j, :J)
                    ψ[idx.(σ)...] *= exp(sgn * 0.5 * β * σ[i] * J * σ[j])
                end
            end
        end
        ψ = cu(ψ)
        ρ = abs.(ψ) .^ 2
        @test ρ / sum(ρ) ≈ ϱ

        @info "Generating MPS from |ρ>"
        rψ = CuMPS(ψ)
        @test dot(rψ, rψ) ≈ 1
        @test_nowarn is_right_normalized(rψ)

        lψ = CuMPS(ψ, :left)
        @test dot(lψ, lψ) ≈ 1

        vlψ = vec(tensor(lψ))
        vrψ = vec(tensor(rψ))

        vψ = vec(ψ)
        vψ /= norm(vψ)

        @test abs(1 - abs(CUDA.dot(vlψ, vrψ))) < ϵ
        @test abs(1 - abs(CUDA.dot(vlψ, vψ))) < ϵ

        @info "Verifying MPS from gates"

        Gψ = CuMPS(ig, control)

        @test_nowarn is_right_normalized(Gψ)
        @test bond_dimension(Gψ) > 1
        @test dot(Gψ, Gψ) ≈ 1
        @test_nowarn verify_bonds(Gψ)

        @test abs(1 - abs(dot(Gψ, rψ))) < ϵ

        @info "Verifying probabilities" L β
        for σ ∈ states
            p = dot(rψ, σ)
            r = dot(rψ, cuproj(σ, rank), rψ)

            @test p ≈ r
            @test Array(ϱ)[idx.(σ)...] ≈ p
        end

        for max_states ∈ [1, N, 2*N, N^2]

            @info "Verifying low energy spectrum" max_states
            @info "Testing spectrum"
            states, prob, pCut = solve(rψ, max_states)
            sp = brute_force(ig, num_states = max_states)

            @info "The largest discarded probability" pCut

            eng = zeros(length(prob))

            for (j, (p, e)) ∈ enumerate(zip(prob, sp.energies))
                σ = states[j, :]
                eng[j] = energy(σ, ig)
                @test log(Array(ϱ)[idx.(σ)...]) ≈ p
            end

            perm = partialsortperm(eng, 1:max_states)
            eng = eng[perm]
            states = states[perm, :]
            prob = prob[perm]
            state = states[1, :]
            @info "The largest discarded probability" pCut
            @test maximum(prob) > pCut
            @info "State with the lowest energy" state
            @info "Probability of the state with the lowest energy" prob[1]
            @info "The lowest energy" eng[1]

        end

    end
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -23.82106505746807 <= q <= 20.87497636142748)

                     add_NL_constraint(m, :(swish(-0.7432570628691102 + 0.9028765283094948*swish(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.9177729976456397*swish(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.6091404444372248*swish(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + 0.7075296204914552*swish(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + -0.9115635329002854*swish(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + swish(-0.2671327260424561 + -0.018499351907995898*swish(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + 0.12529725503678657*swish(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + 0.9839396030143157*swish(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + 0.28500560463058955*swish(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.970334521625404*swish(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + swish(0.40040910594674495 + -0.7454849743358074*swish(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + 0.7108986625228169*swish(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.30349624151173904*swish(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.26057109207128715*swish(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.7858005505148937*swish(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + swish(-0.7120753920031948 + -0.005433419997351141*swish(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.8839514190404207*swish(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + -0.2933551592941672*swish(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.5418193939730438*swish(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.8978172795430157*swish(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) + swish(0.7950909223530527 + -0.6380481687930124*swish(-0.11514299287103702 + 0.5608817277242992*$(x[1]) + 0.6700218006141032*$(x[2]) + -0.9693868545458724*$(x[3])) + -0.8467501205760102*swish(-0.5568784315984976 + 0.281764667612443*$(x[1]) + 0.42150753955170206*$(x[2]) + -0.5829321104114529*$(x[3])) + 0.5001198683847332*swish(-0.04388549286647825 + -0.9010617152865183*$(x[1]) + -0.7540101668519821*$(x[2]) + -0.4758905699629783*$(x[3])) + -0.03855272462979853*swish(-0.511868403698823 + 0.5804725192481515*$(x[1]) + 0.7124288843879394*$(x[2]) + -0.06834010138653612*$(x[3])) + 0.3301524112226937*swish(-0.9281357966735606 + -0.33421670762637223*$(x[1]) + -0.5412233562926474*$(x[2]) + -0.2936515077150612*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module RocketSubjectTest

using Test
using Rocket

@testset "Subject" begin

    struct DummySubjectType end

    struct NotImplementedSubject{T} end
    Rocket.as_subject(::Type{<:NotImplementedSubject{T}}) where T = ValidSubjectTrait{T}()

    struct ActorMissingSubject{T} end
    Rocket.as_subject(::Type{<:ActorMissingSubject{T}})      where T = ValidSubjectTrait{T}()
    Rocket.as_subscribable(::Type{<:ActorMissingSubject{T}}) where T = SimpleSubscribableTrait{T}()

    struct SubscribableMissingSubject{T} end
    Rocket.as_subject(::Type{<:SubscribableMissingSubject{T}})      where T = ValidSubjectTrait{T}()
    Rocket.as_actor(::Type{<:SubscribableMissingSubject{T}})        where T = BaseActorTrait{T}()

    struct ImplementedSubject{T}
        values :: Vector{T}

        ImplementedSubject{T}() where T = new(Vector{T}())
    end

    Rocket.as_subject(::Type{<:ImplementedSubject{T}})      where T = ValidSubjectTrait{T}()
    Rocket.as_actor(::Type{<:ImplementedSubject{T}})        where T = BaseActorTrait{T}()
    Rocket.as_subscribable(::Type{<:ImplementedSubject{T}}) where T = SimpleSubscribableTrait{T}()

    Rocket.on_next!(subject::ImplementedSubject{T}, data::T) where T = push!(subject.values, data)
    Rocket.on_error!(subject::ImplementedSubject, err)               = error(err)
    Rocket.on_complete!(subject::ImplementedSubject)                 = begin end

    function Rocket.on_subscribe!(subject::ImplementedSubject, actor)
        complete!(actor)
        return voidTeardown
    end

    struct ImplementedAutoSubject{T} <: AbstractSubject{T}
        values :: Vector{T}

        ImplementedAutoSubject{T}() where T = new(Vector{T}())
    end

    Rocket.on_next!(subject::ImplementedAutoSubject{T}, data::T) where T = push!(subject.values, data)
    Rocket.on_error!(subject::ImplementedAutoSubject, err)               = error(err)
    Rocket.on_complete!(subject::ImplementedAutoSubject)                 = begin end

    function Rocket.on_subscribe!(subject::ImplementedAutoSubject, actor)
        complete!(actor)
        return voidTeardown
    end

    @testset "as_subject" begin
        # Check if arbitrary dummy type has invalid subject type
        @test as_subject(DummySubjectType) === InvalidSubjectTrait()

        # Check if as_subject returns valid subject type for an implemented subject object
        @test as_subject(ImplementedSubject{Int})     === ValidSubjectTrait{Int}()
        @test as_subject(ImplementedAutoSubject{Int}) === ValidSubjectTrait{Int}()

        @test as_actor(ImplementedSubject{Int})     === BaseActorTrait{Int}()
        @test as_actor(ImplementedAutoSubject{Int}) === BaseActorTrait{Int}()

        @test as_subscribable(ImplementedSubject{Int})     === SimpleSubscribableTrait{Int}()
        @test as_subscribable(ImplementedAutoSubject{Int}) === SimpleSubscribableTrait{Int}()
    end

    @testset "subscribe! as a subscribable" begin
        actor   = void(Int)
        actor_s = void(String)

        # Check if subscribe! throws an error for not valid subject being as a source
        @test_throws InvalidSubscribableTraitUsageError    subscribe!(DummySubjectType(), actor)
        @test_throws InvalidSubscribableTraitUsageError    subscribe!(SubscribableMissingSubject{Int}(), actor)
        @test_throws MissingOnSubscribeImplementationError subscribe!(ActorMissingSubject{Int}(), actor)

        # Check if subscribe! subscribes to a valid subscribable
        @test subscribe!(ImplementedSubject{Int}(), actor)     === voidTeardown
        @test subscribe!(ImplementedAutoSubject{Int}(), actor) === voidTeardown

        # Check if subscribe! throws an error if subscribable and actor data types does not match
        @test_throws InconsistentActorWithSubscribableDataTypesError subscribe!(ImplementedSubject{Int}(), actor_s)
        @test_throws InconsistentActorWithSubscribableDataTypesError subscribe!(ImplementedAutoSubject{Int}(), actor_s)
    end

    @testset "subscribe! as an actor" begin
        source = from(1:5)

        # Check if subscribe! throws an error for not valid subject being as an actor
        @test_throws InvalidActorTraitUsageError           subscribe!(source, DummySubjectType())
        @test_throws InvalidActorTraitUsageError           subscribe!(source, ActorMissingSubject{Int}())
        @test_throws MissingOnNextImplementationError      subscribe!(source, SubscribableMissingSubject{Int}())

        # Check if subscribe! subscribes to a valid subscribable
        subject = ImplementedSubject{Int}()
        @test subscribe!(source, subject) === voidTeardown
        @test subject.values == [ 1, 2, 3, 4, 5 ]

        # Check if subscribe! subscribes to a valid subscribable
        subject = ImplementedAutoSubject{Int}()
        @test subscribe!(source, subject) === voidTeardown
        @test subject.values == [ 1, 2, 3, 4, 5 ]

        # Check if subscribe! throws an error if subscribable and actor data types does not match
        @test_throws InconsistentActorWithSubscribableDataTypesError subscribe!(source, ImplementedSubject{String}())
        @test_throws InconsistentActorWithSubscribableDataTypesError subscribe!(source, ImplementedAutoSubject{String}())
    end

    struct DummyFactory end

    struct NotImplementedSubjectFactory <: AbstractSubjectFactory end

    struct ImplementedSubjectFactory <: AbstractSubjectFactory end

    create_subject(::Type{L}, factory::ImplementedSubjectFactory) where L = ImplementedSubject{L}()

    @testset "test AbstractSubjectFactory" begin
        @test_throws MethodError create_subject(Int, DummyFactory())
        @test_throws MethodError create_subject(Int, NotImplementedSubjectFactory())

        actor = create_subject(String, ImplementedSubjectFactory())
        @test actor isa ImplementedSubject{String}
    end

end

end


text:
using Documenter, DocumenterMarkdown, BioServices

makedocs(
	format	= Markdown(),
	sitename= "BioServices.jl"
)

deploydocs(
    deps 	= Deps.pip("mkdocs", "mkdocs-material", "pygments", "python-markdown-math"),
    repo 	= "github.com/BioJulia/BioServices.jl.git",
    make 	= () -> run(`mkdocs build`),
    target 	= "site"
)


text:
# SD_AIN Functions for Keysight M31XXA/M33XXA Digitizers
export SD_AIN_channelInputConfig, SD_AIN_channelPrescalerConfig,
	SD_AIN_channelPrescalerConfigMultiple, SD_AIN_clockIOconfig,
	SD_AIN_clockGetFrequency, SD_AIN_clockGetSyncFrequency,
	SD_AIN_clockSetFrequency, SD_AIN_clockResetPhase, SD_AIN_triggerIOconfig,
	SD_AIN_triggerIOwrite, SD_AIN_triggerIOread, SD_AIN_DAQconfig,
	SD_AIN_DAQstart, SD_AIN_DAQstop, SD_AIN_DAQflush, SD_AIN_DAQtrigger,
	SD_AIN_DAQstartMultiple, SD_AIN_DAQflushMultiple, SD_AIN_DAQstopMultiple,
	SD_AIN_DAQtriggerMultiple, SD_AIN_DAQtriggerConfig, SD_AIN_DAQanalogTriggerConfig,
	SD_AIN_DAQdigitalTriggerConfig, SD_AIN_channelTriggerConfig,
	SD_AIN_DAQtriggerExternalConfig, SD_AIN_channelPrescaler, SD_AIN_channelFullScale,
	SD_AIN_channelMinFullScale, SD_AIN_channelMaxFullScale, SD_AIN_channelImpedance,
	SD_AIN_channelCoupling, SD_AIN_DAQread, SD_AIN_DAQcounterRead, SD_AIN_DAQnPoints,
	SD_AIN_DAQbufferPoolConfig, SD_AIN_DAQbufferPoolRelease, SD_AIN_DAQbufferAdd,
	SD_AIN_DAQbufferRemove, SD_AIN_DAQbufferGet, SD_AIN_FFT

## int SD_AIN_channelInputConfig(int moduleID, int channel, double fullScale, int impedance, int coupling);
"""
This function configures the input full scale, impedance and coupling as
applicable according to the product.
### Input Impedance
- High Impedance : Input impedance is high --> AIN_IMPEDANCE_HZ = 0
- 50Ω : Input impedance is 50Ω  --> AIN_IMPEDANCE_50 = 1
### Input Coupling
- DC : DC coupling --> AIN_COUPLING_DC = 0
- AC : AC coupling --> AIN_COUPLING_AC = 1
"""
SD_AIN_channelInputConfig(moduleID::Integer, channel::Integer, fullScale::Real,
	impedance::Integer, coupling::Integer) =
	ccall((:SD_AIN_channelInputConfig, klib), Cint, (Cint, Cint, Cdouble, Cint, Cint),
		moduleID, channel, fullScale, impedance, coupling)

## int SD_AIN_channelPrescalerConfig(int moduleID, int channel, int prescaler);
"""
This function configures the input prescaler.
"""
SD_AIN_channelPrescalerConfig(moduleID::Integer, channel::Integer, prescaler::Integer) =
	ccall((:SD_AIN_channelPrescalerConfig, klib), Cint, (Cint, Cint, Cint),
		moduleID, channel, prescaler)

## int SD_AIN_channelPrescalerConfigMultiple(int moduleID, int mask, int prescaler);
SD_AIN_channelPrescalerConfigMultiple(moduleID::Integer, mask::Integer,
	prescaler::Integer) =
	ccall((:SD_AIN_channelPrescalerConfigMultiple, klib), Cint,
		(Cint, Cint, Cint), moduleID, mask, prescaler)

# Clock
## int SD_AIN_clockIOconfig(int moduleID, int clockConfig);
SD_AIN_clockIOconfig(moduleID::Integer, clockConfig::Integer) =
	ccall((:SD_AIN_clockIOconfig, klib), Cint, (Cint, Cint), moduleID, clockConfig)

## double SD_AIN_clockGetFrequency(int moduleID);
"""
This function returns the real hardware clock frequency. It may differ from the
frequency set with the function clockSetFrequency, due to the hardware
frequency resolution.
"""
SD_AIN_clockGetFrequency(moduleID::Integer) =
	ccall((:SD_AIN_clockGetFrequency, klib), Cdouble, (Cint,), moduleID)

## double SD_AIN_clockGetSyncFrequency(int moduleID);
"""
This function returns the frequency of Clock System.
"""
SD_AIN_clockGetSyncFrequency(moduleID::Integer) =
	ccall((:SD_AIN_clockGetSyncFrequency, klib), Cdouble, (Cint,), moduleID)

## double SD_AIN_clockSetFrequency(int moduleID, double frequency, int mode SD_DEFAULT_CLOCK_FREQ_MODE);
"""
This function sets the module clock frequency, see FlexCLK Technology (models
with variable sampling rate).
### CLK Set Frequency Mode
- Low Jitter Mode : The clock system is set to achieve the lowest jitter,
sacrificing tuning speed --> CLK_LOW_JITTER = 0
- Fast Tuning Mode : The clock system is set to achieve the lowest tuning time,
sacrificing jitter performance --> CLK_FAST_TUNE = 1
"""
SD_AIN_clockSetFrequency(moduleID::Integer, frequency::Real, mode::Integer) =
    ccall((:SD_AIN_clockSetFrequency, klib), Cdouble, (Cint, Cdouble, Cint),
        moduleID, frequency, mode)

## int SD_AIN_clockResetPhase(int moduleID, int triggerBehavior, int PXItrigger, double skew SD_DEFAULT_CLOCK_SKEW);
"""
This function sets the module in a sync state, waiting for the first trigger to
reset the phase of the internal clocks CLKsync and CLKsys.
"""
SD_AIN_clockResetPhase(moduleID::Integer, triggerBehavior::Integer, PXItrigger::Integer,
	skew::Real) =
	ccall((:SD_AIN_clockResetPhase, klib), Cint, (Cint, Cint, Cint, Cdouble),
		moduleID, triggerBehavior, PXItrigger, skew)

# Trigger
## int SD_AIN_triggerIOconfig(int moduleID, int direction);
"""
This function configures the trigger connector/line direction and
synchronization/sampling method.
### Trigger I/O Options
- Trigger Output (readable) : TRG operates as a general purpose digital output
signal, which can be written by the user software --> AOU_TRIG_OUT = 0
- Trigger Input : TRG operates as a trigger input, or as general purpose digital
input signal, which can be read by the user software --> AOU_TRIG_IN = 1
"""
SD_AIN_triggerIOconfig(moduleID::Integer, direction::Integer) =
    ccall((:SD_AIN_triggerIOconfig, klib), Cint, (Cint, Cint), moduleID, direction)

## int SD_AIN_triggerIOwrite(int moduleID, int value, int syncMode);
"""
This function sets the trigger output. The trigger must be configured as output
using function triggerIOconfig and I/O Triggers.
### Trigger Synchronization/Sampling Options
- Non-synchronized mode : The trigger is sampled with an internal 100 MHz clock
--> SYNC_NONE = 0
- Synchronized mode : (PXI form factor only) The trigger is sampled using CLK10
--> SYNC_CLK_0 = 1
"""
SD_AIN_triggerIOwrite(moduleID::Integer, value::Integer, syncMode::Integer) =
    ccall((:SD_AIN_triggerIOwrite, klib), Cint, (Cint, Cint, Integer),
		moduleID, value, syncMode)

## int SD_AIN_triggerIOread(int moduleID);
"""
This function reads the trigger input.
"""
SD_AIN_triggerIOread(moduleID::Integer) =
    ccall((:SD_AIN_triggerIOread, klib), Cint, (Cint,), moduleID)

## int SD_AIN_DAQconfig(int moduleID, int nDAQ, int nDAQpointsPerCycle, int nCycles, int triggerDelay, int triggerMode);
"""
This function configures the acquisition of words in two possible reading modes:
- Blocking: Using the function to read the words. DAQread is a blocking function
that is released when the amount of words specified in DAQpoints is acquired or
when timeout elapses. This mode is enabled when a callback function is not
specified (it is set to null).
- Non-blocking: The user specifies a callback function which is called whenever
the DAQeventDataReady event is signaled or when timeout elapses. In the latter
condition, there may be words available, but less than the amount specified in
DAQpoints.
### DAQ Trigger Mode
- Auto (Immediate) : The acquisition starts automatically after a call to
function DAQstart --> AUTOTRIG = 0
- Software / HVI : Software trigger. The acquisition is triggered by the
function DAQtrigger, provided that the DAQ is running. DAQtrigger can be
executed from the user application (VI) or from an HVI --> SWHVITRIG = 1
- Hardware Digital Trigger : Hardware trigger. The DAQ waits for an external
digital trigger --> HWDIGTRIG = 2
- Hardware Analog Trigger : Hardware trigger. The DAQ waits for an external
analog trigger (only products with analog inputs) --> HWANATRIG = 3
"""
SD_AIN_DAQconfig(moduleID::Integer, nDAQ::Integer, nDAQpointsPerCycle::Integer,
	nCycles::Integer, triggerDelay::Integer, triggerMode::Integer) =
	ccall((:SD_AIN_DAQconfig, klib), Cint, (Cint, Cint, Cint, Cint, Cint, Cint),
		moduleID, nDAQ, nDAQpointsPerCycle, nCycles, triggerDelay, triggerMode)

## int SD_AIN_DAQstart(int moduleID, int nDAQ);
"""
This function starts acquisition on the selected DAQ. Acquisition will start
when a trigger is received.
"""
SD_AIN_DAQstart(moduleID::Integer, nDAQ::Integer) =
	ccall((:SD_AIN_DAQstart, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQstop(int moduleID, int nDAQ);
"""
This function stops the words acquisition on the selected DAQ.
"""
SD_AIN_DAQstop(moduleID::Integer, nDAQ::Integer) =
	ccall((:SD_AIN_DAQstop, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQflush(int moduleID, int nDAQ);
"""
This function flushes the acquisition buffers and resets the acquisition counter
included in a Data Acquisition block.
"""
SD_AIN_DAQflush(moduleID::Integer, nDAQ::Integer) =
    ccall((:SD_AIN_DAQflush, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQtrigger(int moduleID, int nDAQ);
"""
This function triggers the acquisition of words in the selected DAQs provided
that they are configured for VI/HVI Trigger.
"""
SD_AIN_DAQtrigger(moduleID::Integer, nDAQ::Integer) =
    ccall((:SD_AIN_DAQtrigger, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQstartMultiple(int moduleID, int DAQmask);
"""
This function starts acquisition on the selected DAQs. Acquisition will start
when a trigger is received.
"""
SD_AIN_DAQstartMultiple(moduleID::Integer, DAQmask::Integer) =
	ccall((:SD_AIN_DAQstartMultiple, klib), Cint, (Cint, Cint),
		moduleID, DAQmask)

## int SD_AIN_DAQflushMultiple(int moduleID, int DAQmask);
"""
This function flushes the acquisition buffers and resets the acquisition counter
included in a Data Acquisition block.
"""
SD_AIN_DAQflushMultiple(moduleID::Integer, DAQmask::Integer) =
    ccall((:SD_AIN_DAQflushMultiple, klib), Cint, (Cint, Cint),
        moduleID, DAQmask)

## int SD_AIN_DAQstopMultiple(int moduleID, int DAQmask);
"""
This function pauses the words acquisition on the selected DAQs.
"""
SD_AIN_DAQstopMultiple(moduleID::Integer, DAQmask::Integer) =
	ccall((:SD_AIN_DAQstopMultiple, klib), Cint, (Cint, Cint),
		moduleID, DAQmask)

## int SD_AIN_DAQtriggerMultiple(int moduleID, int DAQmask);
"""
This function triggers the acquisition of words in the selected DAQs provided
that they are configured for VI/HVI Trigger.
"""
SD_AIN_DAQtriggerMultiple(moduleID::Integer, DAQmask::Integer) =
    ccall((:SD_AIN_DAQtriggerMultiple, klib), Cint, (Cint, Cint),
        moduleID, DAQmask)

## int SD_AIN_DAQtriggerConfig(int moduleID, int nDAQ, int digitalTriggerMode, int digitalTriggerSource, int analogTriggerMask);
SD_AIN_DAQtriggerConfig(moduleID::Integer, nDAQ::Integer, digitalTriggerMode::Integer,
	digitalTriggerSource::Integer, analogTriggerMask::Integer) =
	ccall((:SD_AIN_DAQtriggerConfig, klib), Cint, (Cint, Cint, Cint, Cint, Cint),
		moduleID, nDAQ, digitalTriggerMode, digitalTriggerSource,
		analogTriggerMask)

## int SD_AIN_DAQanalogTriggerConfig(int moduleID, int nDAQ, int analogTriggerMask);
"""
This function configures the analog hardware trigger for the selected DAQ.
"""
SD_AIN_DAQanalogTriggerConfig(moduleID::Integer, nDAQ::Integer, analogTriggerMask::Integer) =
	ccall((:SD_AIN_DAQanalogTriggerConfig, klib), Cint, (Cint, Cint, Cint),
		moduleID, nDAQ, analogTriggerMask)

## int SD_AIN_DAQdigitalTriggerConfig(int moduleID, int nDAQ, int triggerSource, int triggerBehavior);
"""
This function configures the digital hardware triggers for the selected DAQ.
### DAQ Hardware Digital Trigger Source
- External I/O Trigger : The DAQ trigger is a TRG connector/line of the product.
PXI form factor only: this trigger can be synchronized to CLK10
--> TRIG_EXTERNAL = 0
- PXI Trigger : PXI form factor only. The DAQ trigger is a PXI trigger line
and it is synchronized to CLK10 --> TRIG_PXI = 1
### DAQ Hardware Digital Trigger Behavior
- Active High : Trigger is active when it is at level high --> TRIG_HIGH = 1
- Active Low : Trigger is active when it is at level low --> TRIG_LOW = 2
- Rising Edge : Trigger is active on the rising edge --> TRIG_RISE = 3
- Falling Edge : Trigger is active on the falling edge --> TRIG_FALL = 4
"""
SD_AIN_DAQdigitalTriggerConfig(moduleID::Integer, nDAQ::Integer,
	triggerSource::Integer, triggerBehavior::Integer) =
	ccall((:SD_AIN_DAQdigitalTriggerConfig, klib), Cint, (Cint, Cint, Cint, Cint),
		moduleID, nDAQ, triggerSource, triggerBehavior)

## int SD_AIN_channelTriggerConfig(int moduleID, int nChannel, int analogTriggerMode, double threshold);
"""
This function configures the analog trigger block for each channel.
### Analog Trigger Mode
- Rising Edge : Trigger is generated when the input signal is rising and crosses
the threshold --> AIN_RISING_EDGE = 0
- Falling Edge : Trigger is generated when the input signal is falling and
crosses the threshold --> AIN_FALLING_EDGE = 1
- Both Edges : Trigger is generated when the input signal crosses the threshold,
no matter if it is rising or falling --> AIN_BOTH_EDGES = 3
"""
SD_AIN_channelTriggerConfig(moduleID::Integer, nChannel::Integer,
	analogTriggerMode::Integer, threshold::Real) =
	ccall((:SD_AIN_channelTriggerConfig, klib), Cint, (Cint, Cint, Cint,
		Cdouble), moduleID, nChannel, analogTriggerMode, threshold)

## int SD_AIN_DAQtriggerExternalConfig(int moduleID, int nDAQ, int externalSource, int triggerBehavior);
SD_AIN_DAQtriggerExternalConfig(moduleID::Integer, nDAQ::Integer,
	externalSource::Integer, triggerBehavior::Integer) =
	ccall((:SD_AIN_DAQtriggerExternalConfig, klib), Cint, (Cint, Cint, Cint, Cint),
		moduleID, nDAQ, externalSource, triggerBehavior)

## int SD_AIN_channelPrescaler(int moduleID, int channel);
SD_AIN_channelPrescaler(moduleID::Integer, channel::Integer) =
	ccall((:SD_AIN_channelPrescaler, klib), Cint, (Cint, Cint), moduleID, channel)

## double SD_AIN_channelFullScale(int moduleID, int channel);
SD_AIN_channelFullScale(moduleID::Integer, channel::Integer) =
	ccall((:SD_AIN_channelFullScale, klib), Cdouble, (Cint, Cint), moduleID, channel)

## double SD_AIN_channelMinFullScale(int moduleID, int impedance, int coupling);
SD_AIN_channelMinFullScale(moduleID::Integer, impedance::Integer,
	coupling::Integer) =
	ccall((:SD_AIN_channelMinFullScale, klib), Cdouble, (Cint, Cint, Cint),
		moduleID, impedance, coupling)

## double SD_AIN_channelMaxFullScale(int moduleID, int impedance, int coupling);
SD_AIN_channelMaxFullScale(moduleID::Integer, impedance::Integer,
	coupling::Integer) =
	ccall((:SD_AIN_channelMaxFullScale, klib), Cdouble, (Cint, Cint, Cint),
		moduleID, impedance, coupling)

## int SD_AIN_channelImpedance(int moduleID, int channel);
SD_AIN_channelImpedance(moduleID::Integer, channel::Integer) =
	ccall((:SD_AIN_channelImpedance, klib), Cint, (Cint, Cint), moduleID, channel)

## int SD_AIN_channelCoupling(int moduleID, int channel);
SD_AIN_channelCoupling(moduleID::Integer, channel::Integer) =
	ccall((:SD_AIN_channelCoupling, klib), Cint, (Cint, Cint), moduleID, channel)

## int SD_AIN_DAQread(int moduleID, int nDAQ, short *dataBuffer, int nPoints, int timeOut);
"""
This function reads the words acquired with the selected DAQ. It can be used
only after calling the function DAQconfig and when a callback function is not
configured. DAQread is a blocking function released when the configured amount
of words is acquired, or when the configured timeout elapses (if timeout is set
to ”0” , then DAQread waits until DAQpoints are acquired). In the timeout
elapses, there may be words available, but less than the configured amount.
"""
function SD_AIN_DAQread(moduleID::Integer, nDAQ::Integer, nPoints::Integer,
	timeOut::Integer=0)
	dataBuffer = Vector{Cshort}(nPoints) 	# Create an array to retrieve data
	# dataBuffer contains nPoints words. Its size is nPoints × 2bytes/word.
 	val = ccall((:SD_AIN_DAQread, klib), Cint, (Cint, Cint, Ref{Cshort}, Cint,
		Cint), moduleID, nDAQ, dataBuffer, nPoints, timeOut)
	if val > 0
		return dataBuffer
	else
		return val
	end
end

## int SD_AIN_DAQcounterRead(int moduleID, int nDAQ);
"""
This function reads the number of words acquired by the selected DAQ since the
last call to DAQflush or DAQ.
"""
SD_AIN_DAQcounterRead(moduleID::Integer, nDAQ::Integer) =
    ccall((:SD_AIN_DAQcounterRead, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQnPoints(int moduleID, int nDAQ);
SD_AIN_DAQnPoints(moduleID::Integer, nDAQ::Integer) =
	ccall((:SD_AIN_DAQnPoints, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQbufferPoolConfig(int moduleID, int nDAQ, short *dataBuffer, int nPoints, int timeOut SD_DEFAULT_NULL, callbackEventPtr callbackFunction SD_DEFAULT_NULL, void *callbackUserObj SD_DEFAULT_NULL);
"""
This function configures buffer pool that will be filled with the data of the
channel to be transferred to PC.
"""
function SD_AIN_DAQbufferPoolConfig(moduleID::Integer, nDAQ::Integer,
	nPoints::Integer, timeOut::Integer=0)
	dataBuffer = Vector{Cshort}(nPoints)
	val = ccall((:SD_AIN_DAQbufferPoolConfig, klib), Cint,
		(Cint, Cint, Ref{Cshort}, Cint, Cint, Ptr{Void}, Ptr{Void}),
		moduleID, nDAQ, dataBuffer, nPoints, timeOut, C_NULL, C_NULL)
	return val
end

## int SD_AIN_DAQbufferPoolRelease(int moduleID, int nDAQ);
"""
Releases the channel buffer pool and its resources. After this call, user has to
call DAQbufferRemove consecutively to get all buffers back and release them.
"""
SD_AIN_DAQbufferPoolRelease(moduleID::Integer, nDAQ::Integer) =
	ccall((:SD_AIN_DAQbufferPoolRelease, klib), Cint, (Cint, Cint), moduleID, nDAQ)

## int SD_AIN_DAQbufferAdd(int moduleID, int nDAQ, short *dataBuffer, int nPoints);
"""
Adds an additional buffer to the channel’s previously configured pool.
"""
function SD_AIN_DAQbufferAdd(moduleID::Integer, nDAQ::Integer,
	dataBuffer::Vector{Int16}, nPoints::Integer)
	nPoints = length(dataBuffer)
	val = ccall((:SD_AIN_DAQbufferAdd, klib), Cint, (Cint, Cint, Ref{Cshort}, Cint),
		moduleID, nDAQ, dataBuffer, nPoints)
	return val
end

## short* SD_AIN_DAQbufferRemove(int moduleID, int nDAQ);
"""
Ask for a buffer to be removed from the channel buffer pool. If NULL pointer is
returned, no more buffers remains in buffer pool. Returned buffer is a
previously added buffer from user and user has to release/delete it.
"""
function SD_AIN_DAQbufferRemove(moduleID::Integer, nDAQ::Integer)
	ptr = ccall((:SD_AIN_DAQbufferRemove, klib), Ptr{Cshort}, (Cint, Cint),
		moduleID, nDAQ)
	return unsafe_load(ptr)
end

## short* SD_AIN_DAQbufferGet(int moduleID, int nDAQ, int SD_DLL_API_REF_C readPointsOut, int SD_DLL_API_REF_C errorOut);
"""
Gets a filled buffer from the channel buffer pool. User has to call DAQbufferAdd
with this buffer to tell the pool that the buffer can be used again.
"""
function SD_AIN_DAQbufferGet(moduleID::Integer, nDAQ::Integer)
	errorOut = Cint(0)
	readPointsOut = Cint(0)
	ptr = ccall((:SD_AIN_DAQbufferGet, klib), Ptr{Cshort},
		(Cint, Cint, Ref{Cint}, Ref{Cint}), moduleID, nDAQ, readPointsOut, errorOut)
	if errorOut < 0
		return errorOut
	else
		nPoints = readPointsOut
		if nPoints > 0
			return Vector{Cshort}([unsafe_load(ptr, n) for n in 1:nPoints])
		else
			return Vector{Cshort}(0)
		end
	end
end

## int SD_AIN_FFT(int moduleID, int channel, short *data, int size, double *result, int resultSize, double *resultPhase, SD_BOOL dB, int windowType);
"""
Calculates the FFT of data captured by DAQread for the selected channel.
### Window types used in FFT function
- Rectangular : Simplest B-spine window --> WINDOW_RECTANGULAR = 0 (default)
- Bartlett : Hybrid window --> WINDOW_BARTLETT = 1
- Hanning : Side-lobes roll off about 18 dB per octave --> WINDOW_HANNING = 2
- Hamming : Optimized to minimize the maximum nearest side lobe
--> WINDOW_HAMMING = 3
- Blackman : Higher-order generalized cosine windows for applications that
require windowing by the convolution in the frequency-domain
--> WINDOW_BLACKMAN = 4
- Kaiser : Adjustable window maximizing energy concentration in the main lobe
--> WINDOW_KAISER = 5
- Gauss : Adjustable window (can be used for quadratic interpolation in
frequency estimation) --> WINDOW_GAUSS = 6
"""
function SD_AIN_FFT(moduleID::Integer, channel::Integer, data::Vector{Int16},
	dB::Bool=false, windowType::Integer=0)
	c_data = Vector{Cshort}(data) # conversion of data into appropriate Cshort counterpart
	dataSize = length(c_data)
	resultSize = dataSize
	resultMag = Vector{Cdouble}(resultSize)
	resultPhase = Vector{Cdouble}(resultSize)
	val = ccall((:SD_AIN_FFT, klib), Cint, (Cint, Cint, Ref{Cshort}, Cint,
		Ref{Cdouble}, Cint, Ref{Cdouble}, Bool, Cint), moduleID, channel, c_data,
		dataSize, resultMag, resultSize, resultPhase, dB, windowType)
	if val > 0
		return resultMag, resultPhase
	else
		return val
	end
end

# Functions that exist in the manual, but not in the library `SD_AIN.h` file
# ## DAQpause
# """
# This function pauses the words acquisition on the selected DAQ. Acquisition can
# be resumed using DAQresume.
# """
# SD_AIN_DAQpause(moduleID::Integer, nDAQ::Integer) =
#     ccall((:SD_AIN_DAQpause, klib), Cint, (Cint, Cint), moduleID, nDAQ)
#
# ## DAQpauseMultiple
# """
# This function pauses the words acquisition on the selected DAQs. Acquisition can
# be resumed using DAQresume.
# """
# SD_AIN_DAQpauseMultiple(moduleID::Integer, DAQmask::Integer) =
#     ccall((:SD_AIN_DAQpauseMultiple, klib), Cint, (Cint, Cint),
#         moduleID, DAQmask)
#
# ## DAQresume
# """
# This function resumes acquisition on the selected DAQ.
# """
# SD_AIN_DAQresume(moduleID::Integer, nDAQ::Integer) =
# 	ccall((:SD_AIN_DAQresume, klib), Cint, (Cint, Cint), moduleID, nDAQ)
#
# ## DAQresumeMultiple
# """
# This function resumes acquisition on the selected DAQs.
# """
# SD_AIN_DAQresumeMultiple(moduleID::Integer, DAQmask::Integer) =
#     ccall((:SD_AIN_DAQresumeMultiple, klib), Cint, (Cint, Cint),
#         moduleID, DAQmask)