text:
# Autoregressive-Normal
#

#######################################################################
# Model structure and constructor

"""
    struct ARNormal
        hprior::Normal
        hposterior::Normal
        λprior::Gamma
        λposterior::Gamma
    end

AutoRegressive Normal distribution of order `K`. When `K` is `0`, the
distribution is a univariate Normal distribution. Note that `hprior`
and `hposterior` are `K+1` multivariate Normal distribution.

Constructor:

    ARNormal1D(μ₀, Σ₀, a₀, b₀[, μ = missing, Σ = missing, a = missing, b = missing])

The filter order `K` is set to `length(μ₀) - 1`.
"""
struct ARNormal1D{K}
    # Prior/posterior over `h`, the IIR filter
    hprior::Normal
    hposterior::Normal

    # Prior/posterior over `λ`, noise precision
    λprior::Gamma
    λposterior::Gamma
end

function ARNormal1D(
    μ₀::Vector{T},
    Σ₀::Matrix{T},
    a₀::Real,
    b₀::Real;
    μ::Union{Vector{T}, Missing} = missing,
    Σ::Union{Matrix{T}, Missing} = missing,
    a::Union{Real, Missing} = missing,
    b::Union{Real, Missing} = missing
) where T <: AbstractFloat
    K = length(μ₀) - 1
    μ = ismissing(μ) ? copy(μ₀) : μ
    Σ = ismissing(Σ) ? copy(Σ₀) : Σ
    a = ismissing(a) ? a₀ : a
    b = ismissing(b) ? b₀ : b

    a₀ = T(a₀)
    b₀ = T(b₀)
    a = T(a)
    b = T(b)
    ARNormal1D{K}(Normal(μ₀, Σ₀), Normal(μ, Σ), Gamma([a₀], [b₀]), Gamma([a], [b]))
end

"""
    const ARNormal1DSet{K} = Vector{ARNormal1D{K}} where K

Set of AutoRegressive Normal distributions.
"""
const ARNormal1DSet{K} = Vector{ARNormal1D{K}} where K

#######################################################################
# Expected log-likelihood

"""
    loglikelihood(model::ARNormal1D{K}, x::Vector{T})
    loglikelihood(models::ARNormal1DSet{K}, x::Vector{T})
    loglikelihood(models::Vector{ARNormal1DSet{K}}, X::Matrix{T})

Evaluate the expectation of the log-likelihood w.r.t. the variational
posterior. When the model is a (vector of) `ARNomralSet`, the
log-likelihood per component is returned.
"""
function loglikelihood(
    model::ARNormal1D{K},
    x::Vector{T}
) where {K, T <: AbstractFloat}
    r = Regressors1D(K, x)

    # Expectation of the natural parameters
    η_h = gradlognorm(model.hposterior)
    η_λ = gradlognorm(model.λposterior)

    # Sufficient statistics for h
    s1 = hcat([xₜ * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_h = vcat(s1, -.5 * s2)

    # Sufficient statistics for λ
    stats_λ = hcat(-.5 * x.^2 .+ stats_h' * η_h, .5 * ones(T, length(x)))

    stats_λ * η_λ .- .5 * log(2π)
end

function loglikelihood(
    models::ARNormal1DSet{K},
    x::Vector{T}
) where {K, T<:AbstractFloat}
    r = Regressors1D(K, x)

    η_hs = vcat([gradlognorm(model.hposterior)' for model in models]...)
    η_λs = vcat([gradlognorm(model.λposterior)' for model in models]...)

    # Sufficient statistics for h
    s1 = hcat([xₜ * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_h = vcat(s1, -.5 * s2)

    iter = zip(eachrow(η_hs * stats_h), eachrow(η_λs))
    vcat([(hcat(-.5 * x.^2 .+ row, .5 * ones(T, length(x))) * η_λ)'
          for (row, η_λ) in iter]...) .- .5 * log(2π)
end

function loglikelihood(
    models::Vector{ARNormal1DSet{K}},
    X::Matrix{T}
) where {K, T<:AbstractFloat}
    reduce(+, [loglikelihood(models[d], X[d, :]) for d in 1:size(X, 1)])
end

#######################################################################
# Posterior predictive
#
# NOTE: we use the MAP of precision parameters and we integrate over the
#       the filters (h_1, h_2, ...)

"""
    predict(model::ARNormal1D{K}, x::Vector{T})
    predict(model::ARNormal1DSet{K}, x::Vector{T})
    predict(model::Vector{ARNormal1DSet{K}}, x::Vector{T})

Return the logarithm of the posterior predictive distribution for each
frame `x[t]`.
"""
function predict(
    model::ARNormal1D{K},
    x::Vector{T}
) where {K, T<:AbstractFloat}

    r = Regressors1D(K, x)
    λ_map_inv = 1 / mean(model.λposterior)[1]
    m = model.hposterior.μ
    Σ = model.hposterior.Σ

    # Sufficient statistics for h
    s1 = hcat([x̂ₜ for x̂ₜ in r]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)

    ψ = s1' * m
    ϑ_inv = s2' * vec(Σ)

    prec = 1 ./ (λ_map_inv .+ ϑ_inv)

    -.5 .* prec .* (x .- ψ).^2 .+ .5 .* log.(prec) .- .5 * log(2π)
end

function predict(
    models::ARNormal1DSet{K},
    x::Vector{T}
) where {K, T<:AbstractFloat}
    vcat([predict(models[s], x)' for s in 1:length(models)]...)
end

function predict(
    models::Vector{Vector{ARNormal1D{K}}},
    X::Matrix{T}
) where {K, T<:AbstractFloat}
    reduce(+, [predict(models[d], X[d, :]) for d in 1:size(X, 1)])
end

########################################################################
# Accumulate statistics

function accstats_h(
    models::Vector{ARNormal1D{K}},
    x::Vector{T},
    resps::Matrix{T}
) where {K, T<:AbstractFloat}
    r = Regressors1D(K, x)

    # Sufficient statistics
    s1 = hcat([xₜ * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_h = vcat(s1, -.5 * s2)

    accstats = (resps * stats_h')

    [mean(m.λposterior) .* s for (m, s) in zip(models, eachrow(accstats))]
end

function accstats_λ(
    models::Vector{ARNormal1D{K}},
    x::Vector{T},
    resps::Matrix{T}
) where {K, T<:AbstractFloat}
    r = Regressors1D(K, x)

    # Sufficient statistics
    s1 = hcat([xₜ * x̂ₜ for (xₜ, x̂ₜ) in zip(x, r)]...)
    s2 = hcat([vec(x̂ₜ* x̂ₜ') for x̂ₜ in r]...)
    stats_h = vcat(s1, -.5 * s2)

    η_hs = vcat([gradlognorm(model.hposterior)' for model in models]...)

    accstats = Vector{Vector{T}}()
    for (hs, γ)  in zip(eachrow(η_hs * stats_h), eachrow(resps))
        λs = hcat(-.5 * x.^2 .+ hs, .5 * ones(T, length(x)))'
        s = λs * γ
        push!(accstats, s)
    end
    accstats
end

#######################################################################
# Utilities

function save(
    fname::AbstractString,
    model::Vector{Vector{ARNormal1D{K}}}
) where K
    bson(fname, Dict(:model => model, :order=>K, :type=>typeof(model)))
end

function load(
    fname::AbstractString,
)
    data = BSON.load(fname)
    m = data[:model]
    K = data[:order]
    T = data[:type]
    convert(T, m)
end

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
module Auxiliary

using LinearAlgebra

export Indices, FloatRC, Float2, Float3, Int2, Int3
export indvec, threshold, compfloateps
export qraddcols!, qraddcols, lqaddrows

const Indices = Union{Vector{Int},Int,UnitRange{Int},StepRange{Int,K} where K,NTuple{M,Int} where M,Vector{Any},Nothing,Colon}
const FloatRC = Union{T,Complex{T}} where T<:AbstractFloat
const Float2 = Union{Float64,Vector{Float64}}
const Float3 = Union{Float64,Vector{Float64},Vector{Vector{Float64}}}
const Int2 = Union{Int,Vector{Int}}
const Int3 = Union{Int,Vector{Int},Vector{Vector{Int}}}

compfloateps(::Type{S}) where {T<:AbstractFloat,S<:FloatRC{T}} = eps(T)

function indvec(σ::Indices; min::Int=1, max::Int=0)
	isa(σ, Colon) && return collect(min:max)
	isa(σ, Vector{Int}) && return copy(σ)
	isa(σ, Int) && return [σ]
	isa(σ, UnitRange{Int}) && return collect(σ)
	isa(σ, StepRange{Int,S} where S) && return collect(σ)
	isa(σ, Tuple{}) && return Vector{Int}(undef, 0)
	isa(σ, NTuple{N,Int} where N) && return collect(σ)
	isa(σ, Vector{Any}) && return Vector{Int}(undef, 0)
	isa(σ, Nothing) && return Vector{Int}(undef, 0)
end

function threshold(δ::Vector{T}, ε::T, ρ::Int) where {T<:Real}
	r = length(δ)
	if r == 0
		throw(ArgumentError("δ is required to be nonempty"))
	end
	if ε < 0
		throw(ArgumentError("ε is required to be nonnegative"))
	end
	if ρ < 0
		throw(ArgumentError("ρ is required to be nonnegative"))
	end
	if ρ == 0
		ρ = r
	end
	ρ = min(ρ, r)
	push!(δ, convert(T, 0))
	δ = reverse(δ)
	δ = cumsum(δ)
	k = findlast(x -> x ≤ ε, δ)
	k = max(k, r+1-ρ)
	ε = δ[k]
	ρ = r+1-k
	return ε, ρ
end

function qraddcols(Q::Matrix{T}, R::Matrix{T}, U::Matrix{T}) where T<:FloatRC
	p,r = size(Q); rR,q = size(R); pU,s = size(U)
	if rR ≠ r
		throw(ArgumentError("Q and R are incompatible in rank"))
	end
	if pU ≠ p
		throw(ArgumentError("Q and U are incompatible in the first rank"))
	end
	Q = [Q U]
	Q,Z = qr!(Q); Q = Matrix(Q); ρ = size(Z, 1)
	d = diag(view(Z, 1:r, 1:r))
	Q[:,1:r] = Q[:,1:r]*Diagonal(d)
	d = conj(d)
	R = [R Diagonal(d)*Z[1:r,r+1:r+s]; zeros(T, ρ-r, q) Z[r+1:ρ,r+1:r+s]]
	return Q,R
end

function lqaddrows(L::Matrix{T}, Q::Matrix{T}, U::Matrix{T}) where T<:FloatRC
	p,rr = size(L); r,q = size(Q); s,qq = size(U)
	if rr ≠ r
		throw(ArgumentError("L and Q are incompatible in size"))
	end
	if qq ≠ q
		throw(ArgumentError("Q and U are incompatible in size"))
	end
	throw(ErrorException("NOT IMPLEMENTED"))
	return Q,R
end

end

ext: jl


text:
# This file holds the definition of the PersonnelDatabase type. This type holds
#   the database of all personnel members in the military organisation.

# The PersonnelDatabase type requires the Personnel type.
requiredTypes = [ "personnel" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export PersonnelDatabase
type PersonnelDatabase
    # The unique id key field.
    idKey::Symbol

    # The list of attributes.
    attrs::Vector{Symbol}

    # The database.
    dbase::Vector{Personnel}

    # The number of entries.
    persSize::Int

    # This constructor creates a new database.
    function PersonnelDatabase( id::Symbol = :id,
        attrList::Vector{Symbol} = Vector{Symbol}() )
        # Make sure the ID field is in the list of fields.
        push!( attrList, id )

        newBase = new()
        newBase.idKey = id
        newBase.attrs = unique( attrList )
        newBase.dbase = Vector{Personnel}()
        newBase.persSize = 0

        return newBase
    end  # PersonnelDatabase( id, attrdList )

    function PersonnelDatabase( id::String,
        attrList::Vector{Symbol} = Vector{Symbol}() )
        return PersonnelDatabase( Symbol( id ), attrlist )
    end  # PersonnelDatabase( id, fieldList )
end  # PersonnelDatabase


text:
using ClimaCore.Geometry: ⊗
"""
calculate_surface_fluxes_atmos_grid!(integrator)

- TODO: generalize interface for regridding and take land state out of atmos's integrator.p
"""
function calculate_surface_fluxes_atmos_grid!(integrator, T_sfc)
    p = integrator.p
    (; ᶜts, dif_flux_energy, dif_flux_ρq_tot, dif_flux_uₕ, params, Cd, Ch) = p

    Y = integrator.u

    z0m, z0b = (FT(1e-3), FT(1e-5)) # roughness lengths (TODO: import from land)
    # Turbulent surface flux calculation
    tsf =
        constant_T_saturated_surface_coefs_coupled.(
            Spaces.level(ᶜts, 1),
            Geometry.UVVector.(Spaces.level(Y.c.uₕ, 1)),
            Spaces.level(Fields.coordinate_field(Y.c).z, 1),
            FT(0), # TODO: get actual value of z_sfc
            Fields.Field(Fields.field_values(T_sfc), axes(Spaces.level(Y.c, 1))), # remove when same instance issue is resolved
            params,
            z0m, # TODO: get these roughness lengths from land
            z0b,
            Cd,
            Ch,
        )

    # Radiation fluxes
    t = integrator.t
    Rn = FT(10 * sin(t / 2π * 1000)) # TODO: link to the SW and LW fluxes from atmos

    # Total energy flux
    if :ρe in propertynames(Y.c)
        @. dif_flux_energy = #Geometry.WVector(FT(10))
            Geometry.WVector(tsf.shf + tsf.shf + Rn)
    end

    # Moisture mass flux
    if :ρq_tot in propertynames(Y.c)
        @. dif_flux_ρq_tot = Geometry.WVector(tsf.E)
    end

    # Momentum flux

    u_space = axes(tsf.ρτxz) # TODO: delete when "space not the same instance" error is dealt with
    normal = Geometry.WVector.(ones(u_space)) # TODO: this will need to change for topography
    ρ_1 = Fields.Field(Fields.field_values(Fields.level(Y.c.ρ, 1)), u_space) # TODO: delete when "space not the same instance" error is dealt with
    if :uₕ in propertynames(Y.c)
        parent(dif_flux_uₕ) .=  # TODO: remove parent when "space not the same instance" error is dealt with
            parent(
                Geometry.Contravariant3Vector.(normal) .⊗
                Geometry.Covariant12Vector.(Geometry.UVVector.(tsf.ρτxz ./ ρ_1, tsf.ρτyz ./ ρ_1)),
            )
    end

    return nothing
end

function variable_T_saturated_surface_coefs(ts_int, uₕ_int, z_int, z_sfc, T_sfc, params, z0m, z0b) # NB: this currently breaks in SF.jl for neutral conditions

    # get the near-surface thermal state
    T_int = TD.air_temperature(params, ts_int)
    Rm_int = TD.gas_constant_air(params, ts_int)
    ρ_sfc = TD.air_density(params, ts_int) * (T_sfc / T_int)^(TD.cv_m(params, ts_int) / Rm_int) # use ideal gas law and hydrostatic balance to extrapolate for surface density

    q_sfc = TD.q_vap_saturation_generic(params, T_sfc, ρ_sfc, TD.Liquid()) # TODO: assumes all surface is water covered. Generalize!
    ts_sfc = TD.PhaseEquil_ρTq(params, ρ_sfc, T_sfc, q_sfc)

    # wrap state values
    sc = SF.ValuesOnly{FT}(;
        state_in = SF.InteriorValues(z_int, (uₕ_int.u, uₕ_int.v), ts_int),
        state_sfc = SF.SurfaceValues(z_sfc, (FT(0), FT(0)), ts_sfc),
        z0m = z0m,
        z0b = z0b,
    )

    # calculate all fluxes
    tsf = SF.surface_conditions(params, sc, SF.UniversalFunctions.Businger())

    _ρ_liq = FT(1e3)# TODO: use CLIMAParameters
    E = SF.evaporation(sc, params, tsf.Ch) / _ρ_liq
    return (; shf = tsf.shf, lhf = tsf.lhf, E = E, ρτxz = tsf.ρτxz, ρτyz = tsf.ρτyz)
end

function constant_T_saturated_surface_coefs_coupled(ts_int, uₕ_int, z_int, z_sfc, T_sfc, params, z0m, z0b, Cd, Ch)

    # get the near-surface thermal state
    T_int = TD.air_temperature(params, ts_int)
    Rm_int = TD.gas_constant_air(params, ts_int)
    ρ_sfc = TD.air_density(params, ts_int) * (T_sfc / T_int)^(TD.cv_m(params, ts_int) / Rm_int) # use ideal gas law and hydrostatic balance to extrapolate for surface density

    q_sfc = TD.q_vap_saturation_generic(params, T_sfc, ρ_sfc, TD.Liquid()) # TODO: assumes all surface is water covered. Generalize!
    ts_sfc = TD.PhaseEquil_ρTq(params, ρ_sfc, T_sfc, q_sfc)

    # wrap state values
    sc = SF.Coefficients{FT}(;
        state_in = SF.InteriorValues(z_int, (uₕ_int.u, uₕ_int.v), ts_int),
        state_sfc = SF.SurfaceValues(z_sfc, (FT(0), FT(0)), ts_sfc),
        z0m = z0m, #FT(1e-3),
        z0b = z0b, #FT(1e-5),
        Cd = Cd, #FT(0.001),
        Ch = Ch, #FT(0.0001),
    )

    # calculate all fluxes
    tsf = SF.surface_conditions(params, sc, SF.UniversalFunctions.Businger()) # here can specify tol, maxiter

    _ρ_liq = FT(1e3)# TODO: use CLIMAParameters
    E = SF.evaporation(sc, params, tsf.Ch) / _ρ_liq
    return (; shf = tsf.shf, lhf = tsf.lhf, E = E, ρτxz = tsf.ρτxz, ρτyz = tsf.ρτyz)
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
#levels.jl
#Algorithm 5.7 of Random Eigenvalues by Alan Edelman

function levels(n)
    # Plot exact semicircle formula for GUE
    # This is given in formula (5.2.16) in Mehta as the sum of phi_j^2
    # but Christoffel-Darboux seems smarter
    ## Parameter
    # n: Dimension of matrix for which to determine exact semicircle
    x = (-1:0.001:1) * sqrt(2n) * 1.3
    pold = zeros(size(x))   # -1st Hermite Polynomial
    p = ones(size(x))       #  0th Hermite Polynomial
    k = p
    for j = 1:n # Compute the three-term recurrence
        pnew = (sqrt(2) * x .* p - sqrt(j - 1) * pold) / sqrt(j)
        pold = p
        p = pnew
    end
    pnew = (sqrt(2) * x .* p - sqrt(n) * pold)/sqrt(n + 1)

    k = n * p .^ 2 - sqrt(n * (n+1)) * pnew .* pold # Use p.420 of Mehta
    k = k .* exp(-x .^ 2) / sqrt(pi) # Normalize

    # Rescale so that "semicircle" is on [-1, 1] and area is pi/2
    (x / sqrt(2n), k * pi / sqrt(2n))
end


text:
# arithmetic.jl

# Addition, substraction and other functions
for TM in tupleTMs
    @eval begin
        tmdata(f::$TM) = (f.x0, f.dom)

        zero(a::$TM) = $TM(zero(a.pol), zero(a.rem), a.x0, a.dom)
        one(a::$TM) = $TM(one(a.pol), zero(a.rem), a.x0, a.dom)

        # iszero(a::$TM) = iszero(a.pol) && iszero(zero(a.rem))

        findfirst(a::$TM) = findfirst(a.pol)

        ==(a::$TM, b::$TM) =
            a.pol == b.pol && a.rem == b.rem && a.x0 == b.x0 && a.dom == b.dom

        # Addition
        +(a::$TM) = $TM(a.pol, a.rem, a.x0, a.dom)

        function +(a::$TM, b::$TM)
            @assert tmdata(a) == tmdata(b)
            return $TM(a.pol+b.pol, a.rem+b.rem, a.x0, a.dom)
        end

        +(a::$TM, b::T) where {T<:NumberNotSeries} = $TM(a.pol+b, a.rem, a.x0, a.dom)

        +(b::T, a::$TM) where {T<:NumberNotSeries} = $TM(b+a.pol, a.rem, a.x0, a.dom)

        # Substraction
        -(a::$TM) = $TM(-a.pol, -a.rem, a.x0, a.dom)

        function -(a::$TM, b::$TM)
            @assert tmdata(a) == tmdata(b)
            return $TM(a.pol-b.pol, a.rem-b.rem, a.x0, a.dom)
        end

        -(a::$TM, b::T) where {T<:NumberNotSeries} = $TM(a.pol-b, a.rem, a.x0, a.dom)

        -(b::T, a::$TM) where {T<:NumberNotSeries} = $TM(b-a.pol, -a.rem, a.x0, a.dom)

        # Basic division
        function basediv(a::$TM, b::$TM)
            invb = rpa(x->inv(x), b)
            return a * invb
        end

        # Multiplication by numbers
        *(a::$TM, b::T) where {T<:NumberNotSeries} = $TM(a.pol*b, b*a.rem, a.x0, a.dom)

        *(b::T, a::$TM) where {T<:NumberNotSeries} = $TM(a.pol*b, b*a.rem, a.x0, a.dom)

        # Division by numbers
        /(a::$TM, b::T) where {T<:NumberNotSeries} = a * inv(b)

        /(b::T, a::$TM) where {T<:NumberNotSeries} = b * inv(a)

        # Power
        ^(a::$TM, r) = rpa(x->x^r, a)

        ^(a::$TM, n::Integer) = rpa(x->x^n, a)
    end
end

# Multiplication
function *(a::TaylorModel1, b::TaylorModel1)
    @assert tmdata(a) == tmdata(b)

    # Polynomial product with extended order
    order = max(get_order(a), get_order(b))
    aext = Taylor1(copy(a.pol.coeffs), 2*order)
    bext = Taylor1(copy(b.pol.coeffs), 2*order)
    res = aext * bext

    # Returned polynomial
    bext = Taylor1( copy(res.coeffs[1:order+1]) )

    # Bound for the neglected part of the product of polynomials
    res[0:order] .= zero(res[0])
    aux = a.dom - a.x0
    Δnegl = res(aux)

    # Remainder of the product
    Δa = a.pol(aux)
    Δb = b.pol(aux)
    Δ = Δnegl + Δb * a.rem + Δa * b.rem + a.rem * b.rem

    return TaylorModel1(bext, Δ, a.x0, a.dom)
end

function *(a::RTaylorModel1, b::RTaylorModel1)
    @assert tmdata(a) == tmdata(b)

    # Polynomial product with extended order
    order = max(get_order(a), get_order(b))
    aext = Taylor1(copy(a.pol.coeffs), 2*order)
    bext = Taylor1(copy(b.pol.coeffs), 2*order)
    res = aext * bext

    # Returned polynomial
    bext = Taylor1( copy(res.coeffs[1:order+1]) )

    # Bound for the neglected part of the product (properly factorized)
    res = Taylor1(copy(res.coeffs[order+2:2*order+1]), order-1)
    aux = a.dom - a.x0
    Δnegl = res(aux)

    # Remainder of the product
    Δa = a.pol(aux)
    Δb = b.pol(aux)
    V = aux^(order+1)
    Δ = Δnegl + Δb * a.rem + Δa * b.rem + a.rem * b.rem * V

    return RTaylorModel1(bext, Δ, a.x0, a.dom)
end

# Division
function /(a::TaylorModel1, b::TaylorModel1)
    @assert a.x0 == b.x0 && a.dom == b.dom
    return basediv(a, b)
end

function /(a::RTaylorModel1, b::RTaylorModel1)
    @assert tmdata(a) == tmdata(b)

    # DetermineRootOrderUpperBound seems equivalent (optimized?) to `findfirst`
    bk = findfirst(b)
    bk ≤ 0 && return basediv(a, b)

    # In 2.3.12, `a` and `b` are now extended in order by `bk`,
    # but how can we do this without knowing the explicit function
    # that created them?
    #
    # Below we reduce the original order by bk.
    #
    order = get_order(a)
    ared = reducetoorder(
        RTaylorModel1(Taylor1(a.pol.coeffs[bk+1:order+1]), a.rem, a.x0, a.dom), order-bk)
    order = get_order(b)
    bred = reducetoorder(
        RTaylorModel1(Taylor1(b.pol.coeffs[bk+1:order+1]), b.rem, b.x0, b.dom), order-bk)

    return basediv( ared, bred )
end

"""
    reducetoorder(a::RTaylorModel1, m::Integer)

From `a::RTaylorModel1`, it returns a RTaylorModel1 of order `m`.

"""
function reducetoorder(a::RTaylorModel1, m::Integer)
    order = get_order(a)
    @assert order ≥ m ≥ 0

    bpol = Taylor1(copy(a.pol.coeffs))
    zz = zero(bpol[0])
    for ind in 0:m
        bpol[ind] = zz
    end
    bf = bpol(a.dom-a.x0)
    Δ = bf + a.rem * (a.dom-a.x0)^(order-m)
    return RTaylorModel1( Taylor1(copy(a.pol.coeffs[1:m+1])), Δ, a.x0, a.dom )
end

# Same as above, for TaylorModelN
zero(a::TaylorModelN) = TaylorModelN(zero(a.pol), zero(a.rem), a.x0, a.dom)
one(a::TaylorModelN) = TaylorModelN(one(a.pol), zero(a.rem), a.x0, a.dom)

# iszero(a::TaylorModelN) = iszero(a.pol) && iszero(zero(a.rem))

findfirst(a::TaylorModelN) = findfirst(a.pol)

==(a::TaylorModelN, b::TaylorModelN) =
    a.pol == b.pol && a.rem == b.rem && a.x0 == b.x0 && a.dom == b.dom

# Addition and substraction
for op in (:+, :-)
    @eval begin
        $(op)(a::TaylorModelN) = TaylorModelN($(op)(a.pol), $(op)(a.rem), a.x0, a.dom)

        function $(op)(a::TaylorModelN, b::TaylorModelN)
            @assert a.x0 == b.x0 && a.dom == b.dom
            return TaylorModelN($(op)(a.pol,b.pol), $(op)(a.rem,b.rem), a.x0, a.dom)
        end

        $(op)(a::TaylorModelN, b::T) where {T<:NumberNotSeries} =
            TaylorModelN($(op)(a.pol, b), a.rem, a.x0, a.dom)

        $(op)(b::T, a::TaylorModelN) where {T<:NumberNotSeries} =
            TaylorModelN($(op)(b, a.pol), $(op)(a.rem), a.x0, a.dom)
    end
end

# Multiplication
function *(a::TaylorModelN, b::TaylorModelN)
    @assert a.x0 == b.x0 && a.dom == b.dom

    # Polynomial product with extended order
    order = max(get_order(a), get_order(b))
    @assert 2*order ≤ get_order()
    aext = TaylorN(copy(a.pol.coeffs), 2*order)
    bext = TaylorN(copy(b.pol.coeffs), 2*order)
    res = aext * bext

    # Returned polynomial
    bext = TaylorN( copy(res.coeffs[1:order+1]) )

    # Bound for the neglected part of the product of polynomials
    res[0:order] .= zero(res[0])
    aux = a.dom - a.x0
    Δnegl = res(aux)

    # Remainder of the product
    Δa = a.pol(aux)
    Δb = b.pol(aux)
    Δ = Δnegl + Δb * a.rem + Δa * b.rem + a.rem * b.rem

    return TaylorModelN(bext, Δ, a.x0, a.dom)
end

# Multiplication by numbers
function *(b::T, a::TaylorModelN) where {T<:NumberNotSeries}
    pol = a.pol * b
    rem = b * a.rem
    return TaylorModelN(pol, rem, a.x0, a.dom)
end

*(a::TaylorModelN, b::T) where {T<:NumberNotSeries} = b * a

# Basic division
function basediv(a::TaylorModelN, b::TaylorModelN)
    invb = rpa(x->inv(x), b)
    return a * invb
end

# Division by numbers
/(a::TaylorModelN, b::T) where {T<:Number} = a * inv(b)
/(b::T, a::TaylorModelN) where {T<:NumberNotSeries} = b * inv(a)

# Power
^(a::TaylorModelN, r) = rpa(x->x^r, a)
^(a::TaylorModelN, n::Integer) = rpa(x->x^n, a)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PolicySetDefinitionProperties <: SwaggerModel
    policyType::Any # spec type: Union{ Nothing, String } # spec name: policyType
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    metadata::Any # spec type: Union{ Nothing, Any } # spec name: metadata
    parameters::Any # spec type: Union{ Nothing, ParameterDefinitions } # spec name: parameters
    policyDefinitions::Any # spec type: Union{ Nothing, Vector{PolicyDefinitionReference} } # spec name: policyDefinitions
    policyDefinitionGroups::Any # spec type: Union{ Nothing, Vector{PolicyDefinitionGroup} } # spec name: policyDefinitionGroups

    function PolicySetDefinitionProperties(;policyType=nothing, displayName=nothing, description=nothing, metadata=nothing, parameters=nothing, policyDefinitions=nothing, policyDefinitionGroups=nothing)
        o = new()
        validate_property(PolicySetDefinitionProperties, Symbol("policyType"), policyType)
        setfield!(o, Symbol("policyType"), policyType)
        validate_property(PolicySetDefinitionProperties, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(PolicySetDefinitionProperties, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(PolicySetDefinitionProperties, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(PolicySetDefinitionProperties, Symbol("parameters"), parameters)
        setfield!(o, Symbol("parameters"), parameters)
        validate_property(PolicySetDefinitionProperties, Symbol("policyDefinitions"), policyDefinitions)
        setfield!(o, Symbol("policyDefinitions"), policyDefinitions)
        validate_property(PolicySetDefinitionProperties, Symbol("policyDefinitionGroups"), policyDefinitionGroups)
        setfield!(o, Symbol("policyDefinitionGroups"), policyDefinitionGroups)
        o
    end
end # type PolicySetDefinitionProperties

const _property_map_PolicySetDefinitionProperties = Dict{Symbol,Symbol}(Symbol("policyType")=>Symbol("policyType"), Symbol("displayName")=>Symbol("displayName"), Symbol("description")=>Symbol("description"), Symbol("metadata")=>Symbol("metadata"), Symbol("parameters")=>Symbol("parameters"), Symbol("policyDefinitions")=>Symbol("policyDefinitions"), Symbol("policyDefinitionGroups")=>Symbol("policyDefinitionGroups"))
const _property_types_PolicySetDefinitionProperties = Dict{Symbol,String}(Symbol("policyType")=>"String", Symbol("displayName")=>"String", Symbol("description")=>"String", Symbol("metadata")=>"Any", Symbol("parameters")=>"ParameterDefinitions", Symbol("policyDefinitions")=>"Vector{PolicyDefinitionReference}", Symbol("policyDefinitionGroups")=>"Vector{PolicyDefinitionGroup}")
Base.propertynames(::Type{ PolicySetDefinitionProperties }) = collect(keys(_property_map_PolicySetDefinitionProperties))
Swagger.property_type(::Type{ PolicySetDefinitionProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PolicySetDefinitionProperties[name]))}
Swagger.field_name(::Type{ PolicySetDefinitionProperties }, property_name::Symbol) =  _property_map_PolicySetDefinitionProperties[property_name]

const _allowed_PolicySetDefinitionProperties_policyType = ["NotSpecified", "BuiltIn", "Custom", "Static"]

function check_required(o::PolicySetDefinitionProperties)
    (getproperty(o, Symbol("policyDefinitions")) === nothing) && (return false)
    true
end

function validate_property(::Type{ PolicySetDefinitionProperties }, name::Symbol, val)
    if name === Symbol("policyType")
        Swagger.validate_param(name, "PolicySetDefinitionProperties", :enum, val, _allowed_PolicySetDefinitionProperties_policyType)
    end
end


text:
using SafeTestsets, Test

@safetestset "Model" begin include("model.jl") end
@testset "MMA" begin
    @safetestset "Approximation" begin include("approximation.jl") end
    @safetestset "Algorithm" begin include("mma.jl") end
end
@safetestset "AugLag/Percival" begin include("percival.jl") end
@safetestset "AugLag2" begin include("auglag.jl") end
@safetestset "Ipopt" begin include("ipopt.jl") end
@safetestset "NLopt" begin include("nlopt.jl") end


text:
using Test
import ArchGDAL as AG

@testset "test_rasterio.jl" begin
    AG.read("ospy/data4/aster.img") do ds
        @testset "version 1" begin
            band = AG.getband(ds, 1)
            count = 0
            total = 0
            buffer = Array{AG.pixeltype(band)}(undef, AG.blocksize(band)..., 1)
            for (cols, rows) in AG.windows(band)
                AG.rasterio!(ds, buffer, [1], rows, cols)
                data = buffer[1:length(cols), 1:length(rows)]
                count += sum(data .> 0)
                total += sum(data)
            end
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 2" begin
            band = AG.getband(ds, 1)
            count = 0
            total = 0
            buffer = Array{AG.pixeltype(band)}(undef, AG.blocksize(band)..., 1)
            for (cols, rows) in AG.windows(band)
                AG.read!(ds, buffer, [1], rows, cols)
                data = buffer[1:length(cols), 1:length(rows)]
                count += sum(data .> 0)
                total += sum(data)
            end
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 3" begin
            band = AG.getband(ds, 1)
            count = 0
            total = 0
            buffer = Matrix{AG.pixeltype(band)}(undef, AG.blocksize(band)...)
            for (cols, rows) in AG.windows(band)
                AG.read!(ds, buffer, 1, rows, cols)
                data = buffer[1:length(cols), 1:length(rows)]
                count += sum(data .> 0)
                total += sum(data)
            end
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 4" begin
            band = AG.getband(ds, 1)
            count = 0
            total = 0
            buffer = Matrix{AG.pixeltype(band)}(undef, AG.blocksize(band)...)
            for (cols, rows) in AG.windows(band)
                AG.read!(band, buffer, rows, cols)
                data = buffer[1:length(cols), 1:length(rows)]
                count += sum(data .> 0)
                total += sum(data)
            end
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 5" begin
            band = AG.getband(ds, 1)
            count = 0
            total = 0
            xbsize, ybsize = AG.blocksize(band)
            buffer = Matrix{AG.pixeltype(band)}(undef, ybsize, xbsize)
            for ((i, j), (nrows, ncols)) in AG.blocks(band)
                # AG.rasterio!(ds,buffer,[1],i,j,nrows,ncols)
                # AG.read!(band, buffer, j, i, ncols, nrows)
                AG.readblock!(band, j, i, buffer)
                data = buffer[1:nrows, 1:ncols]
                count += sum(data .> 0)
                total += sum(data)
            end
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 6" begin
            band = AG.getband(ds, 1)
            buffer =
                Array{AG.pixeltype(band)}(undef, AG.width(ds), AG.height(ds), 1)
            AG.rasterio!(ds, buffer, [1])
            count = sum(buffer .> 0)
            total = sum(buffer)
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 7" begin
            band = AG.getband(ds, 1)
            buffer =
                Matrix{AG.pixeltype(band)}(undef, AG.width(ds), AG.height(ds))
            AG.read!(band, buffer)
            count = sum(buffer .> 0)
            total = sum(buffer)
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 8" begin
            band = AG.getband(ds, 1)
            buffer =
                Matrix{AG.pixeltype(band)}(undef, AG.width(ds), AG.height(ds))
            AG.read!(ds, buffer, 1)
            count = sum(buffer .> 0)
            total = sum(buffer)
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 9" begin
            band = AG.getband(ds, 1)
            buffer =
                Array{AG.pixeltype(band)}(undef, AG.width(ds), AG.height(ds), 1)
            AG.read!(ds, buffer, [1])
            count = sum(buffer .> 0)
            total = sum(buffer)
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "version 10" begin
            band = AG.getband(ds, 1)
            buffer =
                Array{AG.pixeltype(band)}(undef, AG.width(ds), AG.height(ds), 3)
            AG.read!(ds, buffer)
            count = sum(buffer[:, :, 1] .> 0)
            total = sum(buffer[:, :, 1])
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        # check for calling with Tuple
        @testset "version 11" begin
            band = AG.getband(ds, 1)
            count = 0
            total = 0
            buffer = Array{AG.pixeltype(band)}(undef, AG.blocksize(band)..., 1)
            for (cols, rows) in AG.windows(band)
                AG.rasterio!(ds, buffer, (1,), rows, cols)
                data = buffer[1:length(cols), 1:length(rows)]
                count += sum(data .> 0)
                total += sum(data)
            end
            @test total / count ≈ 76.33891347095299
            @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
        end

        @testset "buffer size" begin
            @test size(AG.read(ds, 1, 0, 0, 20, 10)) === (20, 10)
            @test size(AG.read(ds, [1, 3], 0, 0, 20, 10)) === (20, 10, 2)
            @test size(AG.read(ds, 1, 1:10, 31:50)) === (20, 10)
            @test size(AG.read(ds, [1, 3], 1:10, 31:50)) === (20, 10, 2)
            @test size(AG.read(ds, 1:2)) == (5665, 5033, 2)
            band = AG.getband(ds, 1)
            @test size(AG.read(band, 0, 0, 20, 10)) === (20, 10)
        end

        @testset "Writing to buffers" begin
            band = AG.getband(ds, 1)
            @test AG.pixeltype(band) == UInt8
            xbsize, ybsize = AG.blocksize(band)
            AG.create(
                AG.getdriver("MEM"),
                width = AG.width(band),
                height = AG.height(band),
                nbands = 2,
                dtype = AG.pixeltype(band),
            ) do dsout
                bandout = AG.getband(dsout, 1)

                @testset "writeblock!(rb::RasterBand, xoffset, yoffset, buffer)" begin
                    # We write everything to typemax(UInt8)
                    for ((i, j), (nrows, ncols)) in AG.blocks(bandout)
                        AG.writeblock!(
                            bandout,
                            j,
                            i,
                            fill(typemax(UInt8), ncols, nrows),
                        )
                    end
                    buffer = AG.read(bandout)
                    nnzero = sum(buffer .> 0)
                    @test nnzero == AG.height(bandout) * AG.width(bandout)
                    @test sum(buffer) / nnzero ≈ Float64(typemax(UInt8))

                    # Now we write everything to 0
                    for ((i, j), (nrows, ncols)) in AG.blocks(bandout)
                        AG.writeblock!(bandout, j, i, fill(0x00, ncols, nrows))
                    end
                    @test sum(AG.read(bandout) .> 0) == 0
                end

                @testset "write!(rb::RasterBand, buffer::Matrix[, rows, cols])" begin
                    # We write everything to typemax(UInt8)
                    AG.write!(
                        bandout,
                        fill(
                            typemax(UInt8),
                            AG.height(bandout),
                            AG.width(bandout),
                        ),
                    )
                    buffer = AG.read(bandout)
                    nnzero = sum(buffer .> 0)
                    @test nnzero == AG.height(bandout) * AG.width(bandout)
                    @test sum(buffer) / nnzero ≈ Float64(typemax(UInt8))

                    # Now we write everything to 0
                    AG.write!(
                        bandout,
                        fill(0x00, AG.height(bandout), AG.width(bandout)),
                    )
                    @test sum(AG.read(bandout) .> 0) == 0
                end

                @testset "write!(dataset::Dataset, buffer::Matrix, i::Integer[, rows, cols])" begin
                    # We write everything to typemax(UInt8)
                    AG.write!(
                        dsout,
                        fill(
                            typemax(UInt8),
                            AG.height(bandout),
                            AG.width(bandout),
                        ),
                        1,
                    )
                    buffer = AG.read(bandout)
                    nnzero = sum(buffer .> 0)
                    @test nnzero == AG.height(bandout) * AG.width(bandout)
                    @test sum(buffer) / nnzero ≈ Float64(typemax(UInt8))

                    # Now we write everything to 0
                    AG.write!(
                        dsout,
                        fill(0x00, AG.height(bandout), AG.width(bandout)),
                        1,
                    )
                    @test sum(AG.read(bandout) .> 0) == 0
                end

                @testset "write!(dataset::Dataset, buffer::Array, indices[, rows, cols])" begin
                    # We write everything to typemax(UInt8)
                    AG.write!(
                        dsout,
                        fill(
                            typemax(UInt8),
                            AG.height(dsout),
                            AG.width(dsout),
                            2,
                        ),
                        1:2,
                    )
                    buffer = AG.read(dsout)
                    nnzero = sum(buffer .> 0)
                    @test nnzero == 2 * AG.height(dsout) * AG.width(dsout)
                    @test sum(buffer) / nnzero ≈ Float64(typemax(UInt8))

                    # Now we write everything to 0
                    AG.write!(
                        dsout,
                        fill(0x00, AG.height(dsout), AG.width(dsout), 2),
                        1:2,
                    )
                    @test sum(AG.read(dsout) .> 0) == 0
                end
            end
        end
    end
end

@testset "Complex IO" begin
    a = rand(ComplexF32, 10, 10)

    outpath = "/vsimem/" * "testcomplex1.tif"
    AG.create(
        outpath,
        driver = AG.getdriver("GTiff"),
        width = 10,
        height = 10,
        nbands = 1,
        dtype = ComplexF32,
    ) do ds
        return AG.write!(ds, a, 1)
    end

    AG.read(outpath) do ds_copy
        @test AG.getband(ds_copy, 1) == a
    end
end


text:
module Tutorial
using Revise #useful when revising code frequently

#these are all fairly well established packages
using Distributions, LinearAlgebra, Random, Revise, CUDA, GLM, BenchmarkTools, UnPack, CSV
sleep(0.1)

#a simple ols wrapper
#use parametric types to show high performance code
struct OLS{TX, Ty, Tcoef, T}
  X::TX
  y::Ty
  coef::Tcoef

  #inner construction method
  function OLS(X::TX, y::Ty, #main inputs
    ::Type{T}=eltype(Ty); #types for generic code- not necessary
      addintercept::Bool = true) where {TX,Ty, Tcoef ,T}
    K = size(X,2) + addintercept #set the number of coefficients
    N = size(X,1)
    @assert N ≥ K #checks like these are easy in Julia- if this is not true, it will throw an error

    coef = Vector{Union{Missing,T}}(undef, K) #preallocate
    if addintercept
      #augment a vector of ones to X
      some1s = ones(T, N) |> Ty
      Xout = [some1s X]
      @assert typeof(Xout) == TX
      return new{TX,Ty,typeof(coef), T}([X some1s],y,coef) #create the object
    else
      return new{TX,Ty,typeof(coef), T}(X,y,coef) #create the object
    end
  end
end

function simulatedata(N,K)
  ε = rand(Normal(), N)
  β = 1:K |> collect #true coefficient values
  X = rand(Normal(), N, K-1)
  y = X * β[1:(end-1)] .+ β[end] .+ ε #note the '.' implies a vectorized broadcast operation

  return (X=X,y=y,β=β)
end

function estimatebasic!(ols::OLS)
  ols.coef .= inv(ols.X'*ols.X)*(ols.X'*ols.y) #elementwise assign
  return nothing
end

rms(test,ver) = ((test .- ver).^2 |> mean)^0.5

#Basic regression
function benchbasic(N,K)
  d = simulatedata(N,K)
  ols = OLS(d.X,d.y)

  @info "Running regresion using basic estimation methedology"
  @btime estimatebasic!($ols)
  println("Coef error: $(rms(ols.coef, d.β))")
end

benchbasic(10^6,250)

#can we do better?
function estimatecholesky!(ols::OLS)
  ols.coef .= cholesky!(ols.X'*ols.X)\(ols.X'*ols.y) #elementwise assign
  return nothing
end

#cholesky
function benchcholesky(N,K)
  d = simulatedata(N,K)
  ols = OLS(d.X,d.y)

  @info "Running regresion estimates with cholesky methedology"
  @btime estimatecholesky!($ols)
  println("Coef error: $(rms(ols.coef, d.β))")
end

benchcholesky(10^6,250)
#use the GPU!
function estimatecuda!(ols::OLS{TX,<:Any, <:Any, <:Any}) where TX
  ols.coef .= cholesky!(ols.X'*ols.X)\(ols.X'*ols.y) #elementwise assign
  return nothing
end

function benchcuda(N,K)
  d = simulatedata(N,K)

  ols = OLS(d.X |> CuMatrix{Float64},d.y |> CuVector{Float64})

  @info "Running regresion estimates with CUDA methedology"
  @btime estimatecuda!($ols)
  println("Coef error: $(rms(ols.coef, d.β))")
end

benchcuda(10^6,250)

end #module end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#==============================================================================#
# AWSSES.jl
#
# SES API. See http://aws.amazon.com/documentation/ses/
#
# Copyright OC Technology Pty Ltd 2016 - All rights reserved
#==============================================================================#

__precompile__()

module AWSSES

export ses_send, ses_send_raw, ses_send_attachments

using AWSCore
using SymDict
using Retry

ses(aws::AWSConfig, action, args) = AWSCore.Services.email(aws, action, args)

"""
    ses_send([::AWSConfig]; to=, from=, subject=, body=)

Send plain text email.

```
ses_send(
    to="sam@octech.com.au"
    from="sam@octech.com.au"
    subject="Hi Sam!"
    body="Hello!"
)
```
"""

function ses_send(aws::AWSConfig; to="", from="", subject="", body="")

    ses(aws, "SendEmail", Dict(
        "Source" => from,
        "Destination.ToAddresses.member.1" => to,
        "Message.Subject.Data" => subject,
        "Message.Body.Text.Data" => body))
end

ses_send(; a...) = ses_send(default_aws_config(); a...)

"""
    ses_send_raw([::AWSConfig]; to=, from=, raw=)

Send a raw email.

`raw` must contain both headers and message body.
"""

function ses_send_raw(aws::AWSConfig; to="", from="", raw="")

    ses(aws, "SendRawEmail", Dict(
        "Source" => from,
        "Destination.ToAddresses.member.1" => to,
        "RawMessage.Data" => base64encode(raw)))
end

ses_send_raw(;a...) = ses_send_raw(default_aws_config(); a...)

"""
    ses_send_attachments([::AWSConfig],
                         to=, from=, subject=, body=,
                         attachments=[])

Send an email with attachments.
```
ses_send_attachments(to="sam@octech.com.au", from="sam@octech.com.au"
                     subject="Attachments attached",
                     body="See attached attachments...",
                     [("bar.txt", "text/plain", "bar\\n"),
                      ("foo.txt", "text/plain", "foo\\n")])
```
"""

function  ses_send_attachments(aws::AWSConfig;
                               to="", from="", subject="", body="",
                               attachments=[])

    ses_send_raw(
        aws,
        to=to,
        from=from,
        raw = AWSCore.mime_multipart(
            """
            To: $to
            From: $from
            Subject: $subject
            """,
            [("", "text/plain", body), attachments...]))
end

ses_send_attachments(; a...) = ses_send_attachments(default_aws_config(); a...)

end # module AWSSES

#==============================================================================#
# End of file.
#==============================================================================#


text:
struct SSparseVector{N,Ti,Tv}
    nzind::NTuple{N,Ti}
    nzval::NTuple{N,Tv}
    len::Int
end

Base.length(x::SSparseVector) = x.len

+(x::SSparseVector, y::SSparseVector) = SSparseVector((x.nzind..., y.nzind...), (x.nzval..., y.nzval...), x.len)

Base.map(::typeof(*), x, y::SSparseVector{1}) = SSparseVector(y.nzind, (x[y.nzind[1]] * y.nzval[1],), y.len)

Base.broadcasted(::typeof(*), x, y::SSparseVector{1}) = SSparseVector(y.nzind, (x[y.nzind[1]] * y.nzval[1],), y.len)

Base.broadcasted(::typeof(*), x::SSparseVector{1}, y) = SSparseVector(x.nzind, (x.nzval[1] * y[x.nzind[1]],), x.len)

Base.broadcasted(::typeof(*), x::SSparseVector{1}, y::SSparseVector{1}) = SSparseVector(x.nzind, (x.nzval[1] * y[x.nzind[1]],), x.len)

Base.sum(x::SSparseVector) = sum(x.nzval)

Base.eltype(x::SSparseVector) = eltype(x.nzval)

Base.getindex(x::SSparseVector{1}, i) = x.nzind[1] == i ? x.nzval[1] : zero(eltype(x))

function Base.setindex!(x::Matrix, y::SSparseVector, ::Colon, column::Int)
    map((ind, val) -> x[ind, column] += val, y.nzind, y.nzval)
    return x
end

function Base.circshift(x::SSparseVector{1}, i)
    j = x.nzind[1]
    newj = mod(j+i-1, length(x)) + 1

    return SSparseVector((newj,), x.nzval, length(x))
end

function materialize_hessian(f´´, x)
    n = length(x)
    f´´_x = f´´(x)

    H = zeros(n, n)

    for i in 1:n
        fi = f´´_x(SSparseVector((i,), (1.0,), n))

        for j in 1:n
            H[j,i] = fi(SSparseVector((j,), (1.0,), n))
        end
    end

    return H
end


text:
# Hierarchical Gaussian Filter

module HGF1Benchmark

using Rocket
using ReactiveMP
using GraphPPL
using Distributions

@model function hgf()

    real_k     = constvar(1.0)
    real_w     = constvar(-5.0)
    y_variance = constvar(0.01)
    z_variance = constvar(1.0)

    xt_min_mean = datavar(Float64)
    xt_min_var  = datavar(Float64)

    zt_min_mean = datavar(Float64)
    zt_min_var  = datavar(Float64)

    xt_min ~ NormalMeanVariance(xt_min_mean, xt_min_var)
    zt_min ~ NormalMeanVariance(zt_min_mean, zt_min_var)

    zt ~ NormalMeanVariance(zt_min, z_variance) where { q = q(zt_min)q(z_variance)q(zt) }

    gcv_node, xt ~ GCV(xt_min, zt, real_k, real_w) where { q = q(xt, xt_min)q(zt)q(κ)q(ω) }

    y = datavar(Float64)

    y ~ NormalMeanVariance(xt, y_variance)

    return zt, xt, y, gcv_node, xt_min_mean, xt_min_var, zt_min_mean, zt_min_var
end

function generate_input(rng, n)
    real_k = 1.0
    real_w = -5.0

    z_prev = 0.0
    x_prev = 0.0

    z = Vector{Float64}(undef, n)
    v = Vector{Float64}(undef, n)
    x = Vector{Float64}(undef, n)
    y = Vector{Float64}(undef, n)

    y_variance = 0.01
    z_variance = 1.0

    for i in 1:n
        z[i] = rand(rng, Normal(z_prev, sqrt(z_variance)))
        v[i] = exp(real_k * z[i] + real_w)
        x[i] = rand(rng, Normal(x_prev, sqrt(v[i])))
        y[i] = rand(rng, Normal(x[i], sqrt(y_variance)))

        z_prev = z[i]
        x_prev = x[i]
    end

    return y
end

function benchmark(input)
    iters      = 10
    real_k     = 1.0
    real_w     = -5.0
    y_variance = 0.01
    z_variance = 1.0

    n = length(input)

    ms_scheduler = PendingScheduler()
    fe_scheduler = PendingScheduler()

    mz = keep(Marginal)
    mx = keep(Marginal)
    fe = keep(Float64)

    model, (zt, xt, y, gcv_node, xt_min_mean, xt_min_var, zt_min_mean, zt_min_var) = hgf()

    s_mz = subscribe!(getmarginal(zt) |> schedule_on(ms_scheduler), mz)
    s_mx = subscribe!(getmarginal(xt) |> schedule_on(ms_scheduler), mx)
    s_fe = subscribe!(score(Float64, BetheFreeEnergy(), model, fe_scheduler), fe)

    # Initial prior messages
    current_zt = NormalMeanVariance(0.0, 10.0)
    current_xt = NormalMeanVariance(0.0, 10.0)

    # Prior marginals
    setmarginal!(gcv_node, :y_x, MvNormalMeanCovariance([ 0.0, 0.0 ], [ 5.0, 5.0 ]))
    setmarginal!(gcv_node, :z, NormalMeanVariance(0.0, 5.0))

    for i in 1:n

        for _ in 1:iters
            update!(y, input[i])
            update!(zt_min_mean, mean(current_zt))
            update!(zt_min_var, var(current_zt))
            update!(xt_min_mean, mean(current_xt))
            update!(xt_min_var, var(current_xt))

            release!(fe_scheduler)
        end

        release!(ms_scheduler)

        current_zt = getvalues(mz)[end]
        current_xt = getvalues(mx)[end]
    end

    unsubscribe!(s_mz)
    unsubscribe!(s_mx)
    unsubscribe!(s_fe)

    return map(getvalues, (mz, mx, fe))
end

end


text:
using GLTF
using Test

const TEST_ASSETS_URL = "https://github.com/KhronosGroup/glTF-Asset-Generator/releases/download/v0.6.1/GeneratedAssets-0.6.1.zip"
const TEST_LOCAL_PATH = joinpath(@__DIR__, "download.zip")
const ASSETS_LOCAL_PATH = joinpath(@__DIR__, "assets")

if !isfile(TEST_LOCAL_PATH)
    @info "downloading test assets..."
    download(TEST_ASSETS_URL, TEST_LOCAL_PATH)
end

@static if Sys.iswindows()
    const exe7z = joinpath(Sys.BINDIR, "7z.exe")
    isdir(ASSETS_LOCAL_PATH) || run(`$exe7z x $TEST_LOCAL_PATH -o$ASSETS_LOCAL_PATH`)
else
    isdir(ASSETS_LOCAL_PATH) || run(`unzip -x $TEST_LOCAL_PATH -d $ASSETS_LOCAL_PATH`)
end

@testset "glTF-Asset-Generator" begin
    for (root, dirs, files) in walkdir(ASSETS_LOCAL_PATH)
        for file in files
            if endswith(file, ".gltf")
                path = joinpath(root, file)
                @info "loading $path..."
                @test_nowarn load(path)
            end
        end
    end
end


text:
# ~/zjul/corona/wrld/jmh/src/jR0.jl   Sensitivity
# 200504 006
# 006 : differential data
#      : gtx = dgamma(8.4, 3.8) # Default Discrete Gamma Distribution
#      : Need to change sensitivity settings accordingly !
#      : Identify 8.4/3.8 on plots and files
#
# Adapted to joint processing, i.e. after up.all src/jcor.jl and src/jdfs.jl
# Based on       ~/zjul/corona/R0/jR0td.jl
# Again Based on ~/zr/R0/jR0.R and jtd.R
#
# https://github.com/JuliaStats/Distributions.jl/blob/8bdcda68f903cf377ce37ee4540a3e4d936877e1/src/samplers/multinomial.jl#L53
# https://discourse.julialang.org/t/distributions-jl-using-a-custom-random-number-generator/5106
#

#

using CSV, Missings, DataFrames, Distributions, RandomNumbers, Printf, Plots
println("running jR0: reading data ...")
df_cfm            = CSV.read("dfs/df_cfm.csv")
df_cfm_100        = CSV.read("dfs/df_cfm_100.csv")
df_dcfm           = CSV.read("dfs/df_dcfm.csv")
df_dcfm_100       = CSV.read("dfs/df_dcfm_100.csv")
df_dcfm_100mv     = CSV.read("dfs/df_dcfm_100mv.csv")
#df = df_cfm_100    # 005
#df = df_dcfm_100   # 006 Differential
df = df_dcfm_100mv  # 006 Moving Average
ths = "100"
updated  = @sprintf("%s", last(df_cfm.Date))
plotstring = "Confirmed Cases"

function dgamma(mid, std) # gt Discretized Gamma distribution
    midv = mid  #3.0
    stdv = std  #1.5
    tmax = ceil(midv+10*stdv)
    a = midv*midv/(stdv*stdv)
    s = stdv*stdv / midv
    gd=Gamma(a,s)
    tscale = range(0.5, stop=13.5, step=1)
    gtg = zeros(size(tscale)[1])
    for i in 2:size(gtg)[1] gtg[i] = cdf(gd,tscale[i])-cdf(gd,tscale[i-1]) end
    return gtg
end # function dgamma

function R0( gt, ep)
    Tmax = length(ep)
    nsim = 100
    q = [0.025, 0.975]

    gtpad = gt
    if length(gtpad) < Tmax gtpad = vcat(gtpad, zeros(Tmax-length(gtpad))) end

    P = zeros(Tmax, Tmax)
    p = zeros(Tmax, Tmax)
    ms  = Matrix{Float64}[]
    ms0 = zeros(Tmax, nsim)
    for s in 1:Tmax push!(ms, ms0) end

    # Loop on epidemic duration:
    for s in 2:Tmax
        ms[s] = ms0
        if ep[s] > 0
            gttmp=convert(Array, view(gtpad, s:-1:1))
            wcs = ep[1:s].*gttmp[1:s]
            wcs = wcs./sum(wcs)
            pcs = wcs .* (ep[s]./ep[1:s])
            for i in 1:(s-1) if ep[i]==0 pcs[i] = 0.0 end end
            if ep[s] == 1 pcs[s] = 0.0 end
            P[1:s,s] = pcs
            p[1:s,s] = wcs
            nsims=zeros(s, nsim)

            for i in 1:nsim  nsims[1:s,i] = rand(Multinomial(convert(Int64,ep[s]), p[1:s,s])) end
            ms[s][1:s,1:nsim] = ms[s-1][1:s,1:nsim] + nsims
        else
            P[1:s,s] .= 0.0
            p[1:s,s] .= 0.0
            ms[s][1:s,1:nsim] = ms[s-1][1:s,1:nsim]
        end
    end

    #We now have enough data to compute R (from infection network P)
    #along with its 5% and 95% quantiles (from multiple simulations and p)
    # R.WT<-apply(P,1,sum)  # Hvorfor gir dim=2 korrekt resultat?
    rwt = mapslices(sum, P, dims = [2])
    rcr = rwt./(cumsum(gtpad[1:Tmax]))[Tmax:-1:1]
    if isnan(rcr[length(ep)])  rcr[length(ep)] = 0.0 end

    #Simulated incidence at each time unit is the sum of all cases,
    #stored in the last element of multinom.sim list
    tnf   = ms[length(ep)]
    rsm   = tnf./ep            # Belgium s06 ->NaN divisjon på 0 -> Glatte input data nødvendig.
    rsmcr = rsm./(cumsum(gtpad[1:Tmax]))[Tmax:-1:1]
    rsmcr[length(ep),:] = zeros(nsim)

    # Quantiles:
    qsm    = zeros(Tmax,2)
    qsmcr  = zeros(Tmax,2)
    for s in 1:Tmax
        if (ep[s] == 0)
          rwt[s] <- 0
          rsm[s] <- 0
          rcr[s] <- 0
          rsmcr[s] <- 0
        end
        qsm[s,:]   = quantile(rsm[s,:], q)
        qsmcr[s,:] = quantile(rsmcr[s,:], q)
    end

    # Confidence Interval:
    # multinomial simulations at each time step with the expected value of R.
    correct=true
    cint = zeros(Tmax,2) # colnames?
    if correct == true
        R = rcr
        cint = qsmcr
    else
        R =rwt
        cint = qsm
    end

    # Theoretical Predictions:
    pred = vcat(ep, gt)
    for i in 2:length(ep)+length(gt) pred[i]=0.0 end
    for s in 1:Tmax
        pred[s:s+length(gt)-1] = pred[s:s+length(gt)-1] + R[s] .* ep[s] * gt
    end
    pred = pred[1:Tmax]

return R, pred, cint
end # function R0

# MAIN:

df_R0    = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
df_Pred  = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
df_Clow  = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
df_Chgh  = DataFrame(Array{Union{Missing, Float64}}(missing,size(df)), names(df))
mtxCint =Array{Tuple{Union{Missing,Float64},Union{Missing,Float64}},2}(undef, size(df))
for i in 1:size(mtxCint)[1] for j in 1:size(mtxCint)[2] mtxCint[i,j] = (missing, missing) end end
df_Cint = DataFrame(mtxCint, names(df))

#gtx = dgamma(3.0, 1.5) # Default Discrete Gamma Distribution
gtx = dgamma(8.4, 3.8) # Default Discrete Gamma Distribution
# COUNTRY LOOP:
for j in 1:size(df)[2]
    cname=names(df)[j]
    epx = dropmissing(df, cname)[:,j]
    println(j," ", cname," ", sum(epx), " ",   length(epx))
    R, pred, cint = R0(gtx,epx)
    for i in 1:length(R)-1      df_R0[i,j] = R[i] end
    for i in 1:length(pred)   df_Pred[i,j] = pred[i] end
    for i in 1:size(cint)[1]-1  df_Clow[i,j] = cint[i,1] end
    for i in 1:size(cint)[1]-1  df_Chgh[i,j] = cint[i,2] end
    for i in 1:size(cint)[1]-1  df_Cint[i,j] = (cint[i,1], cint[i,2]) end
end

CSV.write("dfs/df_R0.csv", df_R0)
CSV.write("dfs/df_Pred.csv", df_Pred)
CSV.write("dfs/df_Clow.csv", df_Clow)
CSV.write("dfs/df_Chgh.csv", df_Chgh)
CSV.write("dfs/df_Cint.csv", df_Cint)

# SENSITIVITY
println("jR0: Running Sensitivity Analysis for Norway ...")
j=4 # Country = Norway
cname=names(df)[j]
epx = dropmissing(df, cname)[:,j];
println(j, cname," ", sum(epx), " ", length(epx))
# Mean Value:
s  = 1.5
mlist  = [ 2.5, 3.0, 3.5, 4.0, 5.0, 7.0];
mlistsym = [Symbol(@sprintf("m=%s", mlist[i])) for i in 1:length(mlist)];
df_R0SM = DataFrame(Array{Union{Missing, Float64}}(missing,(size(df)[1],length(mlist))), mlistsym);
for j in 1:length(mlist)
    gtx = dgamma(mlist[j], s)
    R, pred, cint = R0(gtx,epx)
    for i in 1:length(R)-1      df_R0SM[i,j] = R[i] end
end
# Standard Deviation:
m  = 3.5
slist  = [ 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 7.0];
slistsym = [Symbol(@sprintf("s=%s", slist[i])) for i in 1:length(slist)];
df_R0SS = DataFrame(Array{Union{Missing, Float64}}(missing,(size(df)[1],length(slist))), slistsym);
for j in 1:length(slist)
    gtx = dgamma(m, slist[j])
    R, pred, cint = R0(gtx,epx)
    for i in 1:length(R)-1      df_R0SS[i,j] = R[i] end
end
# Write:
CSV.write("dfs/df_R0SM.csv", df_R0SM)
CSV.write("dfs/df_R0SS.csv", df_R0SS)

# PLOTS:
println("jR0: Generating Plots ...")
mtxR0 = Matrix(df_R0)
labels1 = String[]
for i in 1:length(names(df_R0)) push!(labels1, String(names(df_R0)[i])) end
labels = reshape(labels1,1,:)
nmlast = [findfirst( x-> ismissing(x), mtxR0[:,i])-1 for i in 1:size(mtxR0)[2]]

plot(mtxR0, label = labels,
     legend=:outertopright,
     ylabel = "Effective Reproduction Number (R)",
     xlabel = "Days since first "*ths*" cases\nEnd of data: "*updated, xguidefontsize=8,
     title  = "Corona 2020\n(4)"*plotstring*" -Selected Countries")
savefig("figs/cfm_04.png")

plot(10:size(mtxR0)[1],mtxR0[10:end,:], label = labels,
     legend=:outertopright,
     ylabel = "Effective Reproduction Number (R)",
     xlabel = "Days since first "*ths*" cases\nEnd of data: "*updated, xguidefontsize=8,
     annotations=([(nmlast[i], mtxR0[nmlast[i],i],text(labels[i][1:2], :left, 8)) for i in 1:length(nmlast)]),
     title  = "Corona 2020\n(5)"*plotstring*" -Selected Countries")

savefig("figs/cfm_05.png")

plot(30:size(mtxR0)[1],mtxR0[30:end,:], label = labels,
     legend=:outertopright,
     ylabel = "Effective Reproduction Number (R)",
     xlabel = "Days since first "*ths*" cases\nEnd of data: "*updated, xguidefontsize=8,
     annotations=([(nmlast[i], mtxR0[nmlast[i],i],text(labels[i][1:2], :left, 8)) for i in 1:length(nmlast)]),
     title  = "Corona 2020\n(6)"*plotstring*" -Selected Countries")
savefig("figs/cfm_06.png")

# Plot Sensitivity:
mtxR0SM = Matrix(df_R0SM)
labels1 = String[]
for i in 1:length(names(df_R0SM)) push!(labels1, String(names(df_R0SM)[i])) end
labels = reshape(labels1,1,:)
nmlast = [findfirst( x-> ismissing(x), mtxR0SM[:,i])-1 for i in 1:size(mtxR0SM)[2]]
plot(30:size(mtxR0SM)[1],mtxR0SM[30:end,:],
     label = labels,
     legend=:outertopright,
     ylabel = "Effective Reproduction Number (R)",
     xlabel = "Days since first "*ths*" cases\nEnd of data: "*updated, xguidefontsize=8,
     # annotations=([(nmlast[i], mtxR0SM[nmlast[i],i],text(labels[i][1:2], :left, 6)) for i in 1:length(nmlast)]),
     title  = "Corona 2020\n(1)"*plotstring*" -Sensitivity R0 Norway")
savefig("figs5/R0SM.png")

mtxR0SS = Matrix(df_R0SS)
labels1 = String[]
for i in 1:length(names(df_R0SS)) push!(labels1, String(names(df_R0SS)[i])) end
labels = reshape(labels1,1,:)
nmlast = [findfirst( x-> ismissing(x), mtxR0SS[:,i])-1 for i in 1:size(mtxR0SS)[2]]
plot(30:size(mtxR0SS)[1],mtxR0SS[30:end,:],
     label = labels,
     legend=:outertopright,
     ylabel = "Effective Reproduction Number (R)",
     xlabel = "Days since first "*ths*" cases\nEnd of data: "*updated, xguidefontsize=8,
     # annotations=([(nmlast[i], mtxR0SS[nmlast[i],i],text(labels[i][1:2], :left, 8)) for i in 1:length(nmlast)]),
     title  = "Corona 2020\n(6)"*plotstring*" -Sensitity Norway")
savefig("figs5/R0SS.png")
println("jR0: Finished")

# EOF

max_stars_repo_name: dyrehaugen/jdt

max_issues_repo_name: dyrehaugen/jdt

max_forks_repo_name: dyrehaugen/jdt


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

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

This operation cancels a job. Jobs can be cancelled only when they are in the WAITING state.

# Arguments
- `job_id`: The unique identifier for a job.

"""
function cancel_job(JobId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "DELETE",
        "/v1/jobs/$(JobId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function cancel_job(
    JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "DELETE",
        "/v1/jobs/$(JobId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_data_set(asset_type, description, name)
    create_data_set(asset_type, description, name, params::Dict{String,<:Any})

This operation creates a data set.

# Arguments
- `asset_type`: The type of asset that is added to a data set.
- `description`: A description for the data set. This value can be up to 16,348 characters
  long.
- `name`: The name of the data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Tags"`: A data set tag is an optional label that you can assign to a data set when you
  create it. Each tag consists of a key and an optional value, both of which you define. When
  you use tagging, you can also use tag-based access control in IAM policies to control
  access to these data sets and revisions.
"""
function create_data_set(
    AssetType, Description, Name; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "POST",
        "/v1/data-sets",
        Dict{String,Any}(
            "AssetType" => AssetType, "Description" => Description, "Name" => Name
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_data_set(
    AssetType,
    Description,
    Name,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "POST",
        "/v1/data-sets",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "AssetType" => AssetType, "Description" => Description, "Name" => Name
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_event_action(action, event)
    create_event_action(action, event, params::Dict{String,<:Any})

This operation creates an event action.

# Arguments
- `action`: What occurs after a certain event.
- `event`: What occurs to start an action.

"""
function create_event_action(
    Action, Event; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "POST",
        "/v1/event-actions",
        Dict{String,Any}("Action" => Action, "Event" => Event);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_event_action(
    Action,
    Event,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "POST",
        "/v1/event-actions",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("Action" => Action, "Event" => Event), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_job(details, type)
    create_job(details, type, params::Dict{String,<:Any})

This operation creates a job.

# Arguments
- `details`: The details for the CreateJob request.
- `type`: The type of job to be created.

"""
function create_job(Details, Type; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "POST",
        "/v1/jobs",
        Dict{String,Any}("Details" => Details, "Type" => Type);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_job(
    Details,
    Type,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "POST",
        "/v1/jobs",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("Details" => Details, "Type" => Type), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation creates a revision for a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Comment"`: An optional comment about the revision.
- `"Tags"`: A revision tag is an optional label that you can assign to a revision when you
  create it. Each tag consists of a key and an optional value, both of which you define. When
  you use tagging, you can also use tag-based access control in IAM policies to control
  access to these data sets and revisions.
"""
function create_revision(DataSetId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "POST",
        "/v1/data-sets/$(DataSetId)/revisions";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_revision(
    DataSetId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "POST",
        "/v1/data-sets/$(DataSetId)/revisions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_asset(asset_id, data_set_id, revision_id)
    delete_asset(asset_id, data_set_id, revision_id, params::Dict{String,<:Any})

This operation deletes an asset.

# Arguments
- `asset_id`: The unique identifier for an asset.
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
function delete_asset(
    AssetId, DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "DELETE",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_asset(
    AssetId,
    DataSetId,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "DELETE",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation deletes a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

"""
function delete_data_set(DataSetId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "DELETE",
        "/v1/data-sets/$(DataSetId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_data_set(
    DataSetId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "DELETE",
        "/v1/data-sets/$(DataSetId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation deletes the event action.

# Arguments
- `event_action_id`: The unique identifier for the event action.

"""
function delete_event_action(
    EventActionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "DELETE",
        "/v1/event-actions/$(EventActionId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_event_action(
    EventActionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "DELETE",
        "/v1/event-actions/$(EventActionId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_revision(data_set_id, revision_id)
    delete_revision(data_set_id, revision_id, params::Dict{String,<:Any})

This operation deletes a revision.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
function delete_revision(
    DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "DELETE",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_revision(
    DataSetId,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "DELETE",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_asset(asset_id, data_set_id, revision_id)
    get_asset(asset_id, data_set_id, revision_id, params::Dict{String,<:Any})

This operation returns information about an asset.

# Arguments
- `asset_id`: The unique identifier for an asset.
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
function get_asset(
    AssetId, DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_asset(
    AssetId,
    DataSetId,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation returns information about a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

"""
function get_data_set(DataSetId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_data_set(
    DataSetId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation retrieves information about an event action.

# Arguments
- `event_action_id`: The unique identifier for the event action.

"""
function get_event_action(EventActionId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "GET",
        "/v1/event-actions/$(EventActionId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_event_action(
    EventActionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/v1/event-actions/$(EventActionId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation returns information about a job.

# Arguments
- `job_id`: The unique identifier for a job.

"""
function get_job(JobId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "GET", "/v1/jobs/$(JobId)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function get_job(
    JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/jobs/$(JobId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_revision(data_set_id, revision_id)
    get_revision(data_set_id, revision_id, params::Dict{String,<:Any})

This operation returns information about a revision.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
function get_revision(
    DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_revision(
    DataSetId,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation lists a data set's revisions sorted by CreatedAt in descending order.

# Arguments
- `data_set_id`: The unique identifier for a data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
"""
function list_data_set_revisions(
    DataSetId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_data_set_revisions(
    DataSetId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation lists your data sets. When listing by origin OWNED, results are sorted by
CreatedAt in descending order. When listing by origin ENTITLED, there is no order and the
maxResults parameter is ignored.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
- `"origin"`: A property that defines the data set as OWNED by the account (for providers)
  or ENTITLED to the account (for subscribers).
"""
function list_data_sets(; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "GET", "/v1/data-sets"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_data_sets(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/data-sets",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation lists your event actions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"eventSourceId"`: The unique identifier for the event source.
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
"""
function list_event_actions(; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "GET", "/v1/event-actions"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_event_actions(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/event-actions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation lists your jobs sorted by CreatedAt in descending order.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"dataSetId"`: The unique identifier for a data set.
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
- `"revisionId"`: The unique identifier for a revision.
"""
function list_jobs(; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "GET", "/v1/jobs"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_jobs(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET", "/v1/jobs", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_revision_assets(data_set_id, revision_id)
    list_revision_assets(data_set_id, revision_id, params::Dict{String,<:Any})

This operation lists a revision's assets sorted alphabetically in descending order.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
"""
function list_revision_assets(
    DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_revision_assets(
    DataSetId,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags_for_resource(resource-arn)
    list_tags_for_resource(resource-arn, params::Dict{String,<:Any})

This operation lists the tags on the resource.

# Arguments
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.

"""
function list_tags_for_resource(
    resource_arn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "GET",
        "/tags/$(resource-arn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    resource_arn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "GET",
        "/tags/$(resource-arn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation starts a job.

# Arguments
- `job_id`: The unique identifier for a job.

"""
function start_job(JobId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "PATCH", "/v1/jobs/$(JobId)"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function start_job(
    JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "PATCH",
        "/v1/jobs/$(JobId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(resource-arn, tags)
    tag_resource(resource-arn, tags, params::Dict{String,<:Any})

This operation tags a resource.

# Arguments
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.
- `tags`: A label that consists of a customer-defined key and an optional value.

"""
function tag_resource(resource_arn, tags; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "POST",
        "/tags/$(resource-arn)",
        Dict{String,Any}("tags" => tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    resource_arn,
    tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "POST",
        "/tags/$(resource-arn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tags" => tags), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource-arn, tag_keys)
    untag_resource(resource-arn, tag_keys, params::Dict{String,<:Any})

This operation removes one or more tags from a resource.

# Arguments
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.
- `tag_keys`: The key tags.

"""
function untag_resource(
    resource_arn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "DELETE",
        "/tags/$(resource-arn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    resource_arn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "DELETE",
        "/tags/$(resource-arn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_asset(asset_id, data_set_id, name, revision_id)
    update_asset(asset_id, data_set_id, name, revision_id, params::Dict{String,<:Any})

This operation updates an asset.

# Arguments
- `asset_id`: The unique identifier for an asset.
- `data_set_id`: The unique identifier for a data set.
- `name`: The name of the asset. When importing from Amazon S3, the S3 object key is used
  as the asset name. When exporting to Amazon S3, the asset name is used as default target S3
  object key.
- `revision_id`: The unique identifier for a revision.

"""
function update_asset(
    AssetId, DataSetId, Name, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "PATCH",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)",
        Dict{String,Any}("Name" => Name);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_asset(
    AssetId,
    DataSetId,
    Name,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "PATCH",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("Name" => Name), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation updates a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: The description for the data set.
- `"Name"`: The name of the data set.
"""
function update_data_set(DataSetId; aws_config::AbstractAWSConfig=global_aws_config())
    return dataexchange(
        "PATCH",
        "/v1/data-sets/$(DataSetId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_data_set(
    DataSetId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "PATCH",
        "/v1/data-sets/$(DataSetId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

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

This operation updates the event action.

# Arguments
- `event_action_id`: The unique identifier for the event action.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Action"`: What occurs after a certain event.
"""
function update_event_action(
    EventActionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "PATCH",
        "/v1/event-actions/$(EventActionId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_event_action(
    EventActionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "PATCH",
        "/v1/event-actions/$(EventActionId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_revision(data_set_id, revision_id)
    update_revision(data_set_id, revision_id, params::Dict{String,<:Any})

This operation updates a revision.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Comment"`: An optional comment about the revision.
- `"Finalized"`: Finalizing a revision tells AWS Data Exchange that your changes to the
  assets in the revision are complete. After it's in this read-only state, you can publish
  the revision to your products.
"""
function update_revision(
    DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()
)
    return dataexchange(
        "PATCH",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_revision(
    DataSetId,
    RevisionId,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return dataexchange(
        "PATCH",
        "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
struct SimpleNormal{Tm, Tσ}
    m::Tm
    σ::Tσ
end

function AbstractGPs.logpdf(d::SimpleNormal, y::Real)
    m = Zygote.literal_getfield(d, Val(:m))
    σ = Zygote.literal_getfield(d, Val(:σ))
    return -(log(2π) + 2log(σ) + ((y - m) / σ)^2) / 2
end

function ChainRulesCore.rrule(::typeof(AbstractGPs.mean), d::T) where {T<:Normal}
    function mean_pullback(Δ::Real)
        return NoTangent, Tangent{T}(μ=Δ)
    end
    return d.μ, mean_pullback
end

function ChainRulesCore.rrule(::typeof(AbstractGPs.var), d::T) where {T<:Normal}
    function var_pullback(Δ::Real)
        NoTangent, Tangent{T}(σ=2Δ * d.σ)
    end
    return d.σ^2, var_pullback
end

"""
    gaussian_reconstruction_term(
        y::AbstractVector{<:Real},
        σ²::AbstractVector{<:Real},
        m̃::AbstractVector{<:Real},
        σ̃²::AbstractVector{<:Real},
    )

Computes
```julia
∫ logpdf(Normal(f, σ²[n]), y[n]) pdf(Normal(m̃[n], σ̃²[n]), f) df
```
for each `n ∈ eachindex(y)`.
"""
function gaussian_reconstruction_term(
    y::AbstractVector{<:Real},
    Σy::Diagonal{<:Real},
    mq::AbstractVector{<:Real},
    Σq::Diagonal{<:Real},
)
    return sum(
        map(
            (y, σ², m̃, σ̃²) -> logpdf(SimpleNormal(m̃, sqrt(σ²)), y) - σ̃² / (2σ²),
            y, diag(Σy), mq, diag(Σq),
        ),
    )
end

"""
    approx_posterior(
        f::AbstractGP,
        x::AbstractVector,
        η1::AbstractVector{<:Real},
        η2::AbstractVector{<:Real},
    )

Compute the approximate posterior. This is just a posterior GP some kind, which is another
AbstractGP.
"""
function approx_posterior(
    f::AbstractGP,
    x::AbstractVector,
    η1::AbstractVector{<:Real},
    η2::Diagonal{<:Real},
)
    y, Σ = canonical_from_natural(η1, η2)
    return posterior(f(x, Σ + 1e-6 * I), y)
end

"""
    update_approx_posterior(
        f::AbstractGP,
        x::AbstractVector,
        η1::AbstractVector{<:Real},
        η2::AbstractVector{<:Real},
        r,
        ρ::Real,
    )
"""
function update_approx_posterior(
    f::AbstractGP,
    x::AbstractVector,
    η1::AbstractVector{<:Real},
    η2::Diagonal{<:Real},
    r,
    ρ::Real,
)
    # Check that the step size makes sense.
    (ρ > 1 || ρ <= 0) && throw(error("Bad step size"))

    # Compute the approximate posterior marginals.
    mq, σ²q = mean_and_var(approx_posterior(f, x, η1, η2)(x))

    # Compute both of the expectation parameters.
    m1, m2 = expectation_from_canonical(mq, Diagonal(σ²q))

    # Compute the gradient w.r.t. both of the expectation parameters. This is equivalent to
    # the natural gradient w.r.t. the natural parameters.
    g1, g2 = Zygote.gradient(
        (m1, m2) -> begin
            y, σ² = canonical_from_expectation(m1, m2)
            return r(y, σ²)
        end,
        m1, m2,
    )

    # Perform a step of NGA in the first natural pseudo-observation vector.
    η1_new = (1 - ρ) .* η1 .+ ρ .* g1
    η2_new = (1 - ρ) .* η2 .+ ρ .* g2

    return η1_new, η2_new, sqrt(sum(abs2, g1 - η1) + sum(abs2, g2 - η2))
end

"""
    optimise_approx_posterior(
        f::AbstractGP,
        x::AbstractVector,
        η1::AbstractVector{<:Real},
        η2::AbstractVector{<:Real},
        ∇r,
        ρ::Real;
        max_iterations=1_000,
        gtol=1e-8,
    )
"""
function optimise_approx_posterior(
    f::AbstractGP,
    x::AbstractVector,
    η1::AbstractVector{<:Real},
    η2::Diagonal{<:Real},
    r,
    ρ::Real;
    max_iterations=1_000,
    tol=1e-8,
)
    # Perform initial iteration.
    η1, η2, delta_norm = update_approx_posterior(f, x, η1, η2, r, ρ)
    iteration = 1

    # Iterate further until convergence met or max iterations exceeded.
    while delta_norm > tol && iteration < max_iterations
        η1, η2, delta_norm = update_approx_posterior(f, x, η1, η2, r, ρ)
        iteration += 1
    end
    return η1, η2, iteration, delta_norm
end

"""
    elbo(
        f::AbstractGP,
        x::AbstractVector,
        η1::AbstractVector{<:Real},
        η2::AbstractVector{<:Real},
        r,
    )

Compute the evidence lower bound associated with the GP `f`, with reconstruction term `r`,
variational parameters `η1` and `η2`, and collection of inputs `x`.
"""
function AbstractGPs.elbo(
    f::AbstractGP,
    x::AbstractVector,
    η1::AbstractVector{<:Real},
    η2::Diagonal{<:Real},
    r,
)
    ỹ, σ̃² = canonical_from_natural(η1, η2)
    mq, σ²q = mean_and_var(approx_posterior(f, x, η1, η2)(x))
    return logpdf(f(x, σ̃²), ỹ) + r(mq, Diagonal(σ²q)) - gaussian_reconstruction_term(ỹ, σ̃², mq, Diagonal(σ²q))
end

abstract type AbstractIntegrater end

"""
    GaussHermiteQuadrature(num_points::Int)

Specify that the expectation required for the reconstruction term is to be computed using
Gauss-Hermite quadrature with `num_points` quadrature points.
"""
struct GaussHermiteQuadrature <: AbstractIntegrater
    num_points::Int
end

"""
    UnivariateFactorisedLikelihood(build_lik)

A likelihood function which factorises across inputs.
Applying this likelihood to any vector `f::AbstractVector{<:Real}` yields a `Product`
distribution.
Exposing this structure is useful when working with quadrature methods.
"""
struct UnivariateFactorisedLikelihood{Tf}
    build_lik::Tf
end

function (l::UnivariateFactorisedLikelihood)(x::AbstractVector{<:Real})
    return AbstractGPs.Product(map(l.build_lik, x))
end

function build_reconstruction_term(
    integrater::GaussHermiteQuadrature,
    latent_gp::LatentGP,
    y::AbstractVector,
)
    integrands = build_integrands(latent_gp, y)
    function reconstruction_term(mq::AbstractVector{<:Real}, Σq::Diagonal{<:Real})
        return sum(batch_quadrature(integrands, mq, map(sqrt, diag(Σq)), integrater.num_points))
    end
    return reconstruction_term
end

function build_integrands(
    latent_gp::LatentGP{<:AbstractGP, <:UnivariateFactorisedLikelihood},
    y::AbstractVector,
)
    lik = latent_gp.lik.build_lik
    return map(y_ -> (f -> logpdf(lik(f), y_)), y)
end

"""
    optimize_elbo(
        build_latent_gp,
        integrater::AbstractIntegrater,
        x::AbstractVector,
        y::AbstractVector,
        θ0::AbstractVector{<:Real},
        optimiser,
        options,
    )

Optimise the elbo w.r.t. model parameters.

# Arguments
- `build_latent_gp`: a unary function accepting an `AbstractVector{<:Real}` returning a
    LatentGP
- `integrater`::AbstractIntegrater: an object specifying how to perform inference
- `x::AbstractVector`: a collection of inputs
- `y::AbstractVector`: a collection of outputs
- `θ0::AbstractVector{<:Real}`: initial model parameter values
- `optimiser`: an optimiser that can be passed to `Optim.optimize`
- `options`: options to be passed to `Optim.optimize`

# Returns
- `AbstractGP`: the approximate posterior at the optimum
- `NamedTuple`: results summary containing various useful bits of information
"""
function optimize_elbo(
    build_latent_gp,
    integrater::AbstractIntegrater,
    x::AbstractVector,
    y::AbstractVector,
    θ0::AbstractVector{<:Real},
    optimiser,
    options,
)
    # Initialise variational parameters. We'll be mutating these to enable warm-starts
    # during each outer iteration of the algorithm.
    __η1 = zeros(length(x))
    __η2 = Diagonal(-ones(length(x)))

    function objective(θ::AbstractVector{<:Real})

        # Unflatten θ and build the model at the current hyperparameters.
        l = build_latent_gp(θ)
        r = build_reconstruction_term(integrater, l, y)

        # Optimise the approximate posterior. Drop the gradient because we're
        # differentiating through the optimum.
        η1_opt, η2_opt = Zygote.ignore() do
            η1, η2, iters, delta = optimise_approx_posterior(
                l.f, x, __η1, __η2, r, 1; tol=1e-4,
            )
            __η1 .= η1
            __η2 .= η2
            return η1, η2
        end

        # Compute the negation of the elbo.
        return -elbo(l.f, x, η1_opt, η2_opt, r)
    end

    training_results = Optim.optimize(
        objective, θ -> only(Zygote.gradient(objective, θ)), θ0, optimiser, options;
        inplace=false,
    )

    l = build_latent_gp(training_results.minimizer)
    r = build_reconstruction_term(integrater, l, y)
    η1, η2, iters, delta = optimise_approx_posterior(l.f, x, __η1, __η2, r, 1; tol=1e-4)
    approx_post = approx_posterior(l.f, x, η1, η2)

    results_summary = (
        training_results=training_results,
        η1=__η1,
        η2=__η2,
        iters=iters,
        delta=delta,
    )

    return approx_post, results_summary
end