text:
# Kernel implementations.

function evalkernel(p, q::Vector{T},
                    θ::AdditiveVarianceKernelType{KT,T2})::T where {T,KT,T2}
    #
    out = θ.base_gain * evalkernel(p, q, θ.base_kernel_parameters)

    if norm(p-q) < θ.zero_tol
        return out + θ.additive_function(p)
    end

    return out
end

# SOS kernel evaluation.
function evalkernel(p, q::Vector{T}, θ::SOSKernelType{KT})::T where {T,KT}

    out = evalkernel(p, q, θ.base_params)
    return out^2
end

# KR kernel
function evalkernel(p, q::Vector{T}, θ::ElementaryKRKernelType{KT,T})::T where {T,KT}

    out_v = evalkernel(p[1:end-1], q[1:end-1], θ.θ_a)
    out_x = evalkernel(p[end], q[end], θ.θ_canonical)

    return out_v*out_x
end

# Adaptive kernel evaluation.
function evalkernel(p, q::Vector{T},
            θ::AdaptiveKernelType{KT} )::T where {T,KT}

    # if norm(p-q) < 1e-15 #zero_tol
    #     return one(T)
    # end

    warpfunc = θ.warpfunc
    p3 = warpfunc(p)
    q3 = warpfunc(q)
    r_new = p3-q3

    r = p-q
    #τ² = LinearAlgebra.dot(r,r) + LinearAlgebra.dot(r_new,r_new)#(p3-q3)*(p3-q3)
    #τ = sqrt(τ²)
    τ = sqrt(sum( r[d]^2 for d = 1:length(r) ) + r_new^2)

    #println("a: τ = ", τ)
    return evalkernel(τ, θ.canonical_params)
end

function evalkernel(p, q::Vector{T},
            θ::FastAdaptiveKernelType{KT,T,D} )::T where {T,KT,D}

    sum_r_warp_sq = sum( (θ.s[i]*(θ.warpfuncs[i](p) - θ.warpfuncs[i](q)))^2
                            for i = 1:length(θ.warpfuncs))

    # sb = one(T) - sum(θ.s)
    # @assert one(T) >= sb > zero(T)
    sb = one(T)

    sum_r_sq = sum( (sb*(p[d]-q[d]))^2 for d = 1:length(q) )

    τ = sqrt( sum_r_sq + sum_r_warp_sq )

    return evalkernel(τ, θ.canonical_kernel)
end

function evalkernel(p, q::Vector{T},
            θ::AdaptiveKernelDPPType{KT},
            zero_tol::T = eps(T)*2 )::T where {T,KT}

    if norm(p-q) < zero_tol
        return one(T) + (θ.warpfunc(p))^2
    end

    warpfunc = θ.warpfunc
    p3 = warpfunc(p)
    q3 = warpfunc(q)
    r_new = p3-q3

    r = p-q
    τ² = LinearAlgebra.dot(r,r) + LinearAlgebra.dot(r_new,r_new)#(p3-q3)*(p3-q3)
    τ = sqrt(τ²)

    return evalkernel(τ, θ.canonical_params)
end

# Adaptive kernel evaluation.
function evalkernel(p, q::Vector{T},
            θ::AdaptiveKernelMultiWarpType{KT},
            zero_tol::T = eps(T)*2 )::T where {T,KT}

    if norm(p-q) < zero_tol
        return one(T)
    end

    return evalkernel(θ.warpfuncs, θ.a, p, q, θ.canonical_params)
end

function evalkernel(p, q::Vector{T},
            θ::AdaptiveKernelMultiWarpDPPType{KT},
            zero_tol::T = eps(T)*2 )::T where {T,KT}

    if norm(p-q) < zero_tol
        #return one(T) + sum( θ.a[m]*(θ.warpfuncs[m](p))^2 for m = 1:length(θ.warpfuncs) )
        return one(T) + θ.self_gain*sum( θ.a[m]*abs(θ.warpfuncs[m](p)) for m = 1:length(θ.warpfuncs) )
    end

    return evalkernel(θ.warpfuncs, θ.a, p, q,
                            θ.canonical_params)
end

function evalkernel(warpfuncs::Vector{Function},
                    weights::Vector{T},
                    p, q::Vector{T},
                    θ_canonical::KT)::T where {T,KT}

    warp_terms = zero(T)
    for m = 1:length(warpfuncs)
        p3 = warpfuncs[m](p)
        q3 = warpfuncs[m](q)

        r_new = p3-q3
        warp_terms += weights[m]*LinearAlgebra.dot(r_new,r_new)
    end

    r = p-q
    τ² = LinearAlgebra.dot(r,r) + warp_terms
    τ = sqrt(τ²)

    return evalkernel(τ, θ_canonical)
end

# Adaptive kernel evaluation.
function evalkernel(p, q::Vector{T},
            θ::AdaptiveKernelΨType{KT})::T where {T,KT}

    if p == q
        return one(T)
    end

    r_new = θ.ψ(p,q)

    r = p-q
    τ² = LinearAlgebra.dot(r,r) + LinearAlgebra.dot(r_new,r_new)#(p3-q3)*(p3-q3)
    τ = sqrt(τ²)
#println("a: τ = ", τ)
    return evalkernel(τ, θ.canonical_params)
end

####### Brownian Bridge kernels. Use tensor product for multivariate input. Defined only on [0,1].

function evalkernel(p, q::T, θ::BrownianBridge10)::T where T<: Real
    return min(p,q)-p*q
end

# function evalkernel(x, z::T, θ::BrownianBridge20)::T where T<: Real
#     if x < z
#         return x*(1-z)*(x^2 +z^2 -2*z)/6
#     end
#
#     return (1-x)*z*(x^2 +z^2 -2*x)/6
# end

function evalkernel(x, z::T, θ::BrownianBridge1ϵ)::T where T<: Real
    ϵ = θ.ϵ
    denominator = ϵ*sinh(ϵ)
    numerator = sinh(ϵ*min(x,z))*sinh(ϵ*(1-max(x,z)))

    return numerator/denominator
end

function evalkernel(x, z::T, θ::BrownianBridge2ϵ)::T where T<: Real
    ϵ = θ.ϵ

    numerator = exp(-ϵ*(x+z))
    denominator = 4*ϵ^3*(exp(2*ϵ)-1)^2
    multiplier = numerator/denominator

    term1 = exp(2*ϵ)*(2*ϵ-ϵ*(x+z)-1)
    term2 = exp(4*ϵ)*(ϵ*(x+z)+1)
    term3 = exp(2*ϵ*(1+x+z))*(2*ϵ-ϵ*(x+z)+1)
    term4 = exp(2*ϵ*(x+z))*(ϵ*(x+z)-1)
    term5 = exp(2*ϵ*(2+min(x,z)))*(-ϵ*abs(x-z)-1)
    term6 = exp(2*ϵ*max(x,z))*(-ϵ*abs(x-z)+1)
    term7 = exp(2*ϵ*(1+min(x,z)))*(1-2*ϵ+ϵ*abs(x-z))
    term8 = exp(2*ϵ*(1+max(x,z)))*(1+2*ϵ-ϵ*abs(x-z))

    return multiplier*(term1+term2+term3+term4+term5+term6+term7+term8)
end

# tensor product kernels.
function evalkernel(p, q::Vector{T}, θ::BrownianBridgeKernelType)::T where T<: Real
    return prod( evalkernel(p[d],q[d],θ) for d = 1:length(q) )
end

function evalkernel(p, q::Vector{T}, θ::BrownianBridgeSemiInfDomain)::T where T<: Real
    return prod( evalkernel(p[d],q[d],θ) for d = 1:length(q) )
end

function evalkernel(p, q::Vector{T}, θ::BrownianBridgeCompactDomain)::T where T<: Real
    return prod( evalkernel(p[d],q[d],θ) for d = 1:length(q) )
end

# equation 4.2 from McCourt's "An Introduction to the Hilbert-Schmidt SVD using iterated Brownian bridge kernels".
# function evalkernel(x, z::T, θ::BrownianBridge20)::T where T<: Real
#     β::Int = 2
#
#     factor1 = (-1)^(β-1)*2^(2*β-1)/factorial(2*β)
#     factor2 = Bernoullipolynomial4(abs(x-z)/2) - Bernoullipolynomial4((x+z)/2)
#     return factor1*factor2
# end

# an un-numbered equation shortly after McCourt's equation 4.3.
function evalkernel(x, z::T, θ::BrownianBridge20)::T where T<: Real

    if z < x
        return -1/6*z*(1-x)*(x^2+z^2-2*x)
    end

    return -1/6*x*(1-z)*(x^2+z^2-2*z)
end

function evalkernel(x_in, z_in::T, θ::BrownianBridge20CompactDomain, d::Int = 1)::T where T<: Real

    # Map input, with ϵ as buffer from the end points.
    ϵ = 1e-6
    x_tmp = interval2itpindex(x_in, θ.a[d], θ.b[d], 10)
    z_tmp = interval2itpindex(z_in, θ.a[d], θ.b[d], 10)

    # assign a value of zero for inputs that are not within [1,10].
    if !(1 < x_tmp < 10)
        return zero(T)
    end

    if !(1 < z_tmp < 10)
        return zero(T)
    end

    # Map inputs to (0,1)
    x = (x_tmp -1)/(9 + ϵ)
    z = (z_tmp -1)/(9 + ϵ)

    # evaluate the Brownian Bridge 20 kernel.
    if z < x
        return -1/6*z*(1-x)*(x^2+z^2-2*x)
    end

    return -1/6*x*(1-z)*(x^2+z^2-2*z)
end

function evalkernel(x_in, z_in::T, θ::BrownianBridgeSemiInfDomain{KT})::T where {T<: Real, KT}

    # map [0,Inf) to [0,0.5).
    x = x_in/(2*(1+x_in))
    z = z_in/(2*(1+z_in))

    return evalkernel(x, z, θ.θ_base)
end

function evalkernel(x_in, z_in::T, θ::BrownianBridgeCompactDomain{KT,T})::T where {T<: Real, KT}

    # map [0,Inf) to [0,0.5).
    x = Utilities.convertcompactdomain(x_in, θ.a, θ.b, zero(T), one(T))
    z = Utilities.convertcompactdomain(z_in, θ.a, θ.b, zero(T), one(T))

    return evalkernel(x, z, θ.θ_base)
end

######## Conventional kernel function implementations
# All kernels are normalized to have maximum value of 1
# Support is r <= 1.
function evalkernel(x1, x2::Vector{T}, θ::StationaryKernelType)::T where T

    # norm(x1-x2)
    #r = x1-x2
    #τ² = LinearAlgebra.dot(r,r)
    #τ = sqrt(τ²)
    τ = norm(x1-x2)

#println("s: τ = ", τ)
    return evalkernel(τ,θ)
end

# 1D version. For speed.
function evalkernel(x1, x2::T, θ::StationaryKernelType)::T where T
    τ = abs(x1-x2)
    #τ = x1-x2

    return evalkernel(τ,θ)
end

# This is the canonical kernel: spline 34.
# outputs 1 at τ = 0, outputs 0 at τ >= 1.
function evalkernel(τ::T, θ::Spline34KernelType) where T

    r::T = τ*θ.a[1]

    tmp = one(T)-r
    if sign(tmp) < zero(T)
        return zero(T)
    end

    tmp = tmp^6
    out::T = ( (35.0*r^2 + 18.0*r + 3.0)*tmp )/3.0

    #@assert isfinite(out)
    return out
end

# This is the canonical kernel: spline 12.
function evalkernel(τ::T, θ::Spline12KernelType) where T

    r::T = τ*θ.a[1]

    tmp = one(T)-r
    if sign(tmp) < zero(T)
        return zero(T)
    end

    tmp = tmp^3
    out::T = (3.0*r + 1.0)*tmp

    #@assert isfinite(out)
    return out
end

# This is the canonical kernel: spline 32.
function evalkernel(τ::T, θ::Spline32KernelType) where T

    r::T = τ*θ.a[1]

    tmp = one(T)-r
    if sign(tmp) < zero(T)
        return zero(T)
    end

    tmp = tmp^4
    out::T = (4.0*r + 1.0)*tmp

    #@assert isfinite(out)
    return out
end

# This is the canonical kernel: 1D Gaussian.
function evalkernel(τ::T, θ::GaussianKernel1DType)::T where T <: Real

    # a is ϵ^2 in "Stable eval of Gaussian radial basis func interpolants", 2012.
    out::T = exp(-θ.ϵ_sq[1]*τ^2)

    #@assert isfinite(out)
    return out
end

# algebraic sigmoid's univariate rational quadratic kernel.
function evalkernel(τ::T, θ::RationalQuadraticKernelType)::T where T <: Real
    denominator = sqrt(θ.a[1]+τ^2)^3

    #return θ.a/denominator
    #return sqrt(θ.a)/denominator
    return sqrt(θ.a[1])^3/denominator
end

function evalkernel(τ::T, θ::TunableRationalQuadraticKernelType)::T where T <: Real
    denominator = sqrt(θ.a[1]+τ^2)^3

    #return θ.a/denominator
    #return sqrt(θ.a)/denominator
    return θ.w[1]*sqrt(θ.a[1])^3/denominator
end

function evalkernel(τ::T, θ::ModulatedSqExpKernelType)::T where T <: Real

    out::T = exp(-θ.ϵ_sq[1]*τ^2)*cos(θ.ν*τ)

    return out
end

function evalkernel(x, z::Vector{T}, θ::ModulatedSqExpKernelType)::T where T <: Real

    r = x-z
    τ = norm(r)

    out::T = exp(-θ.ϵ_sq[1]*τ^2)*cos(dot(r,θ.ν[1]))

    return out
end

function evalkernel(x, z::Vector{T}, θ::AdaptiveModulatedSqExpKernelType)::T where T <: Real
    ϕ = θ.ϕ
    ψ = θ.ψ

    r1 = [x; ϕ(x)] - [z; ϕ(z)]
    τ1 = norm(r1)

    r2 = [x; ψ(x)] - [z; ψ(z)]
    τ2 = norm(r2)

    #out::T = exp(-θ.ϵ_sq*τ1^2)*cos(dot(r2,θ.ν))
    out::T = exp(-θ.ϵ_sq[1]*τ1^2)*cos(r2[end]*θ.ν[end])

    return out
end

# Adaptive kernel for conditional univariate pdf evaluation:
# p(x[end] | x[1:end-1]).
function evalkernel(p, q::Vector{T},
                    θ::CPDF1DKernelType{T,KT})::T where {T,KT}

    d = length(p)

    # k(p,q)
    factor1 = evalkernel(p[1:d-1], q[1:d-1], θ.B)

    # h(p,q)
    factor2 = evalkernel(p[d], q[d], θ.A)

    return factor1*factor2
end

function evalkernel(x, z::Vector{T}, θ::KRWarpKernelType{KTv,KTx})::T where {T,KTv,KTx}

    out::T = evalkernel(x[1:end-1], z[1:end-1], θ.θ_v)
    return out*evalkernel(x[end], z[end], θ.θ_x)
end

# M is such that 1 <= M <= length(x)-1 == length(z)-1
function evalkernel(x, z::Vector{T}, θ::DoubleProductKernelType{KT})::T where {T,KT}
    M = θ.M
    return evalkernel(x[1:M], z[1:M], θ.θ_base)*evalkernel(x[M+1:end], z[M+1:end], θ.θ_base)
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
using Test
using Ripserer
using Suppressor
using StaticArrays

using Ripserer: circumcenter_radius2

include("interfacetest.jl")

@static if Sys.iswindows()
    @test_broken begin
        Alpha(points)
        true
    end
else
    @testset "circumcenter_radius2" begin
        c1, r2 = circumcenter_radius2([SA[1.0, 2.0, 3.0, 4.0]])
        @test c1 == [1.0, 2.0, 3.0, 4.0]
        @test r2 == 0

        c2, r2 = circumcenter_radius2([SA[1, 1, 1], SA[-1, -1, -1]])
        @test c2 == [0, 0, 0]
        @test r2 == 3.0

        c3, r3 = circumcenter_radius2([SA[1, 0], SA[-1, 0], SA[0, √3]])
        @test c3 ≈ [0, 1 / √3]
        @test r3 ≈ (2 / √3)^2

        c4, r4 = circumcenter_radius2([SA[1, 0], SA[0, 0], SA[-1, 0]])
        c4 = [NaN, NaN]
        r4 = Inf
    end

    @testset "Compare to Cechmate" begin
        # This test compares the results to cechmate: https://github.com/scikit-tda/cechmate/
        points = [(t * sinpi(t), t * cospi(t)) for t in range(0, 2; length=20)]

        d0_cechmate_deaths =
            2 .* [
                0.05263158,
                0.05805554,
                0.06761044,
                0.07982644,
                0.09366810,
                0.10851510,
                0.12400688,
                0.13992943,
                0.15615105,
                0.17258742,
                0.18918257,
                0.20589813,
                0.22270697,
                0.23958947,
                0.25653108,
                0.27352082,
                0.29055025,
                0.30761277,
                0.32470317,
                Inf,
            ]
        d1_cechmate_births = 2 .* [0.83075319, 0.86707034, 0.91920916, 0.89414800]
        d1_cechmate_deaths = 2 .* [0.83143060, 0.86803555, 0.92415392, 0.89985412]

        d0, d1, d2 = ripserer(Alpha(points); dim_max=2)
        @test all(iszero, birth.(d0))
        @test death.(d0) ≈ d0_cechmate_deaths atol = 1e-7
        @test birth.(d1) ≈ d1_cechmate_births atol = 1e-7
        @test death.(d1) ≈ d1_cechmate_deaths atol = 1e-7
        @test isempty(d2)
    end

    @testset "Interface" begin
        points = [(t * sinpi(t), t * cospi(t)) for t in range(0, 2; length=20)]
        test_filtration(Alpha, points)
    end
end

ext: jl


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""
    setdiff(x::Interval, y::Interval)

Calculate the set difference `x ∖ y`, i.e. the set of values
that are inside the interval `x` but not inside `y`.

Returns an array of intervals.
The array may:

- be empty if `x ⊆ y`;
- contain a single interval, if `y` overlaps `x`
- contain two intervals, if `y` is strictly contained within `x`.
"""
function setdiff(x::Interval{T}, y::Interval{T}) where T
    intersection = x ∩ y

    isempty(intersection) && return [x]
    intersection == x && return [emptyinterval(T)]

    if x.lo == intersection.lo
        if isfinite(intersection.hi)
            return [unsafe_interval(intersection.hi, x.hi)]
        else
            return [emptyinterval(T)]
        end
    end
    if x.hi == intersection.hi
        if isfinite(intersection.lo)
            return [unsafe_interval(x.lo, intersection.lo)]
        else
            return [emptyinterval(T)]
        end
    end

    return [unsafe_interval(x.lo, y.lo),
            unsafe_interval(y.hi, x.hi)]
end
setdiff(x::ValidInterval, ::ValidInterval) = [x]

"""
    distance(a::Interval, b::Interval)

Hausdorff distance between two intervals.
"""
distance(a::Interval, b::Interval) = max(abs(a.lo-b.lo), abs(a.hi-b.hi))

const where_bisect = 0.49609375

"""
    bisect(X::Interval, α=0.49609375)

Split the interval `X` at position α; α=0.5 corresponds to the midpoint.
Returns a tuple of the new intervals.
"""
function bisect(X::Interval, α=where_bisect)
    @assert 0 ≤ α ≤ 1

    m = mid(X, α)

    return (unsafe_interval(X.lo, m), unsafe_interval(m, X.hi))
end

"""
    mid(a::Interval, α=0.5)

Find an intermediate point at a relative position `α`` in the interval `a`.
The default is the true midpoint at `α = 0.5`.

Assumes 0 ≤ α ≤ 1.

Warning: if the parameter `α = 0.5` is explicitly set, the behavior differs
from the default case if the provided `Interval` is not finite, since when
`α` is provided `mid` simply replaces `+Inf` (respectively `-Inf`) by `prevfloat(+Inf)`
(respecively `nextfloat(-Inf)`) for the computation of the intermediate point.
"""
function mid(a::Interval{T}, α) where T

    lo = (a.lo == -Inf ? nextfloat(T(-Inf)) : a.lo)
    hi = (a.hi == +Inf ? prevfloat(T(+Inf)) : a.hi)

    β = T(α)

    midpoint = β * (hi - lo) + lo
    isfinite(midpoint) && return midpoint
    #= Fallback in case of overflow: hi - lo == +Inf.
       This case can not be the default one as it does not pass several
       IEEE1788-2015 tests for small floats.
    =#
    return (1 - β) * lo + β * hi
end

mid(::EmptyInterval{T}, α) where T = T(NaN)


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Sphere

"""
    Sphere <: SimpleShape

A [`SimpleShape`](@ref) representing a sphere.

This is a unitary sphere centered in the origin. A generic sphere can be specified by applying a [`Transformation`](@ref).

# Members

- `transformation::Transformation`: the `Transformation` associated with the sphere.
- `material::Material`: the [`Material`](@ref) of the spere.
"""
Base.@kwdef struct Sphere <: SimpleShape
    transformation::Transformation = Transformation()
    material::Material = Material()
end

@doc """
    Sphere(transformation::Transformation, material::Material)

Constructor for a [`Sphere`](@ref) instance.
""" Sphere(::Transformation, ::Material)

@doc """
    Sphere(transformation::Transformation = Transformation(),
           material::Material = Material())

Constructor for a [`Sphere`](@ref) instance.
""" Sphere(; ::Transformation, ::Material)

function get_all_ts(::Type{Sphere}, ray::Ray)
    # compute intersection
    origin_vec = convert(Vec, ray.origin)
    a = norm²(ray.dir)
    b = 2f0 * origin_vec ⋅ ray.dir
    c = norm²(origin_vec) - 1f0
    Δ = b^2 - 4f0 * a * c
    Δ < 0 && return Vector{Float32}()
    sqrt_Δ = sqrt(Δ)
    t_1 = (-b - sqrt_Δ) / (2f0 * a)
    t_2 = (-b + sqrt_Δ) / (2f0 * a)

    # @assert isfinite(t_1)
    # @assert isfinite(t_2)

    return [t_1, t_2]
end

function get_t(::Type{Sphere}, ray::Ray)
    # compute intersection
    origin_vec = convert(Vec, ray.origin)
    a = norm²(ray.dir)
    b = 2f0 * origin_vec ⋅ ray.dir
    c = norm²(origin_vec) - 1f0
    Δ = b^2 - 4f0 * a * c
    Δ < 0 && return Inf32
    sqrt_Δ = sqrt(Δ)
    t_1 = (-b - sqrt_Δ) / (2f0 * a)
    t_2 = (-b + sqrt_Δ) / (2f0 * a)
    # nearest point
    if ray.tmin < t_1 < ray.tmax
        return t_1
    elseif ray.tmin < t_2 < ray.tmax
        return t_2
    else
        return Inf32
    end
end

function get_uv(::Type{Sphere}, point::Point)
    u = atan(point[2], point[1]) / (2f0 * π)
    u = u >= 0 ? u : u+1f0
    v = acos(clamp(point[3], -1f0, 1f0)) / π
    Vec2D(u, v)
end

function get_normal(::Type{Sphere}, point::Point, ray::Ray)
    normal = convert(Normal, point) |> normalize
    (normal ⋅ ray.dir < 0) ? normal : -normal
end


text:
using Random, Plots
include("conditional-value-at-risk.jl")

function scale(x::Float64, low::Float64, high::Float64)
    return x * (high - low) + low
end

function random_probability_distribution(rng::AbstractRNG, n_neg::Int, n_pos::Int, low::Real, high::Real)
    low < high || throw(DomainError(""))
    x_neg = rand(rng, n_neg)
    x_neg = scale.(x_neg, low, 0.0)
    x_pos = rand(rng, n_pos)
    x_pos = scale.(x_pos, 0.0, high)
    x = [x_neg; x_pos]
    f = rand(rng, n_neg + n_pos)
    f = f / sum(f)
    i = sortperm(x)
    return x[i], f[i]
end

function hair(plt, x, y)
    plot!(plt, [x, x], [0, y], linewidth=3, color=:grey, alpha=0.5, label=false)
end

function plot_distribution(x, f, x_α; plt=plot())
    for (x2, f2) in zip(x, f)
        hair(plt, x2, f2)
    end
    tail = x .≤ x_α
    head = x .> x_α
    plot!(plt, ylims = (0, 1.15 * maximum(f)))
    plot!(plt, x[tail], f[tail],
        linewidth=0, markershape=:circle,
        markercolor=:darkred, label="Tail")
    plot!(plt, x[head], f[head],
        linewidth=0, markershape=:circle,
        label="")
    return plt
end

function plot_distributions(x, f, α)
    mean = sum(x.*f)
    VaR = value_at_risk(x, f, α)
    CVaR = conditional_value_at_risk(x, f, α)

    # Probability distribution
    plt1 = plot(
        title="Probability distribution",
        xlabel="x", ylabel="f(x)",
        legend=false
    )
    plot_distribution(x, f, VaR; plt=plt1)
    plot!(plt1, [mean], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:green,
        label="E")
    plot!(plt1, [VaR], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:darkblue,
        label="VaR")
    plot!(plt1, [CVaR], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:darkorange,
        label="CVaR")

    # Cumulative probability distribution
    plt2 = plot(
        title="Cumulative distribution",
        xlabel="x", ylabel="F(x)",
        legend=false
    )
    plot_distribution(x, cumsum(f), VaR; plt=plt2)
    plot!(plt2, x, [α for _ in x],
        linewidth=2, label="α",
        linecolor=:darkorange)
    plot!(plt2, [VaR], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:darkblue,
        label="VaR")

    # Stacked plots
    plt = plot(plt1, plt2,
        layout=(2, 1), legend=:outerright, size=(720, 480))
end

function plot_VaR(x, f, α; plt=plot(size=(720, 480), legend=false, title="Value-at-Risk"))
    # Cumulative distribution
    F2 = cumsum(f)
    F1 = [0; F2[1:end-1]]

    plt = plot!(plt, xlabel="p")

    for i in LinearIndices(x)
        # Integral from zero to probability level α
        if F2[i] <= α
            plot!(plt, [F1[i], F2[i]], [x[i], x[i]], color=false, fill=(0, :darkgreen), fillalpha=0.9)
        end
        if F1[i] <= α <= F2[i]
            plot!(plt, [F1[i], α], [x[i], x[i]], color=false, fill=(0, :darkgreen), fillalpha=0.9)
            plot!(plt, [α, F2[i]], [x[i], x[i]], color=false, fill=(0, :darkblue), fillalpha=0.5)
        end
        # Value-at-risk
        plot!(plt, [F1[i], F2[i]], [x[i], x[i]], color=:black)
    end

    # Inverse of cumlative distribution function
    plot!(plt, F2, x, linewidth=0, markershape=:circle, markercolor=:black, markersize=3)
    plot!(plt, F1, x, linewidth=0, markershape=:circle, markercolor=:white, markersize=3)

    # Probability level α
    plot!(plt, [α], [0], marker=:circle, color=:darkorange)

    return plt
end

function plot_CVaR(x, f; plt=plot(size=(720, 480), legend=false, title="Conditional Value-at-Risk"))
    αs = sort([0:0.01:1; f])
    plot!(plt, xlabel="α")
    plot!(plt, [0, 1], repeat([minimum(x)], 2), linewidth=2, linestyle=:dash, alpha=0.5)
    plot!(plt, [0, 1], repeat([sum(x .* f)], 2), linewidth=2, linestyle=:dash, alpha=0.5)
    plot!(plt, αs, [conditional_value_at_risk(x, f, α) for α in αs], linewidth=2)
    return plt
end

max_stars_repo_name: jaantollander/ConditionalValueAtRisk

max_issues_repo_name: jaantollander/ConditionalValueAtRisk

max_forks_repo_name: jaantollander/ConditionalValueAtRisk


text:
# This file makes plots for investigating the Hellinger metric
using PyPlot, PyCall
using LaTeXStrings
@pyimport mpl_toolkits.axes_grid1 as mpl_kit
ImageGrid = mpl_kit.ImageGrid
@pyimport matplotlib.patches as patches
@pyimport matplotlib.ticker as ticker

include("self_confidence.jl")
pygui(false) # turn off gui, instead just write to files

function hellinger_surf(m_rng::Float64,s1::Float64;m2::Float64=0.,s2::Float64=1.)
    return hellinger_normal(m_rng,s1,m2,s2)
end
function sq_surf(m_rng::Float64,f_rng::Float64;s1::Float64=1.,m2::Float64=0.,s2::Float64=1.,return_raw_sq::Bool=false,alpha=1/2)
    SQ = X3(Normal(m_rng,s1),Normal(m2,s2),global_rwd_range=[f_rng],return_raw_sq=return_raw_sq,alpha=alpha)
    if return_raw_sq
        return SQ[2][:H_scaled]
    else
        return SQ
    end
end

sq_rng_m = [-3.;3.]
sq_rng_f = [0.1;2.]
ar = diff(sq_rng_m)./diff(sq_rng_f)[1]

m_rng_sq = linspace(sq_rng_m[1],sq_rng_m[2],100)
f_rng = linspace(sq_rng_f[1],sq_rng_f[2],100)

# Solver Quality Surface
if true
    fig = plt[:figure](figsize=(6.0,8.0))

    grid = ImageGrid(fig, 111,          # as in plt.subplot(111)
                         nrows_ncols=(3,2),
                         direction="row",
                         aspect = true,
                         axes_pad=0.65,
                         share_all=true,
                         cbar_location="right",
                         cbar_mode="each",
                         cbar_size="3%",
                         cbar_pad=0.10,
                         )

    as = [1.,1/2,1/8]
    hgs_a1 = [sq_surf(x,y,return_raw_sq=true,alpha=as[1]) for x in m_rng_sq, y in f_rng]'
    hgs_a2 = [sq_surf(x,y,return_raw_sq=true,alpha=as[2]) for x in m_rng_sq, y in f_rng]'
    hgs_a3 = [sq_surf(x,y,return_raw_sq=true,alpha=as[3]) for x in m_rng_sq, y in f_rng]'

    sq_cmap = :viridis
    cont_cmap = :binary
    im1 = grid[1][:imshow](hgs_a1,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[sq_rng_m;sq_rng_f])
    grid[1][:contour](m_rng_sq,f_rng,hgs_a1,cmap=cont_cmap,linewidths=0.5)
    grid[1][:set_ylabel](L"f, \alpha=1")
    grid[1][:set_title](L"q vs. $\Delta\mu$, and $f$")
    grid[1][:cax][:colorbar](im1)
    grid[1][:set_aspect](ar)

    im3 = grid[3][:imshow](hgs_a2,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[sq_rng_m;sq_rng_f])
    grid[3][:contour](m_rng_sq,f_rng,hgs_a2,cmap=cont_cmap,linewidths=0.5)
    grid[3][:set_ylabel](L"f,\alpha=1/2")
    grid[3][:cax][:colorbar](im3)
    grid[3][:set_aspect](ar)

    im5 = grid[5][:imshow](hgs_a3,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[sq_rng_m;sq_rng_f])
    grid[5][:contour](m_rng_sq,f_rng,hgs_a3,cmap=cont_cmap,linewidths=0.5)
    grid[5][:set_xlabel](L"\Delta\mu")
    grid[5][:set_ylabel](L"f,\alpha=1/8")
    grid[5][:cax][:colorbar](im5)
    grid[5][:set_aspect](ar)

    sqs_a1 = [sq_surf(x,y,return_raw_sq=false,alpha=as[1]) for x in m_rng_sq, y in f_rng]'
    sqs_a2 = [sq_surf(x,y,return_raw_sq=false,alpha=as[2]) for x in m_rng_sq, y in f_rng]'
    sqs_a3 = [sq_surf(x,y,return_raw_sq=false,alpha=as[3]) for x in m_rng_sq, y in f_rng]'

    im2 = grid[2][:imshow](sqs_a1,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[sq_rng_m;sq_rng_f],vmin=0.,vmax=2)
    grid[2][:contour](m_rng_sq,f_rng,sqs_a1,cmap=cont_cmap,linewidths=0.5)
    grid[2][:set_title](L"$x_Q$ vs. $\Delta\mu$, and $f$")
    grid[2][:cax][:colorbar](im2)
    grid[2][:set_aspect](ar)
    im4 = grid[4][:imshow](sqs_a2,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[sq_rng_m;sq_rng_f],vmin=0.,vmax=2)
    grid[4][:contour](m_rng_sq,f_rng,sqs_a2,cmap=cont_cmap,linewidths=0.5)
    grid[4][:cax][:colorbar](im4)
    grid[4][:set_aspect](ar)
    im6 = grid[6][:imshow](sqs_a3,origin="lower",cmap=sq_cmap,interpolation="bilinear",extent=[sq_rng_m;sq_rng_f],vmin=0.,vmax=2)
    grid[6][:contour](m_rng_sq,f_rng,sqs_a3,cmap=cont_cmap,linewidths=0.5)
    grid[6][:set_xlabel](L"$\Delta\mu$")
    grid[6][:cax][:colorbar](im6)
    grid[6][:set_aspect](ar)

    savefig("figs/sq_surf.pdf",transparent=true,dpi=300)
end

if true
    res = 100
    m_max = 5.
    m_min = -5.
    m_rng = linspace(m_min,m_max,res)
    s_max = 5.0
    s_min = 0.1
    s_rng = linspace(s_min,s_max,res)

    ar = diff([m_min,m_max])./diff([s_min,s_max])[1]

    # Hellinger Distance Surface
    fig2, ax2 = PyPlot.subplots()
    fig2[:set_size_inches](4.0,4.0)

    hs = [hellinger_surf(x,y) for x in m_rng, y in s_rng]'

    cx2 = ax2[:imshow](hs,origin="lower",interpolation="bilinear",extent=[m_min,m_max,s_min,s_max],clim=[0.,1.])
    ax2[:set_aspect](ar)
    fig2[:colorbar](cx2,shrink=0.75)
    ax2[:set_xlabel](L"$\Delta\mu$")
    ax2[:set_ylabel](L"$\sigma_P$")
    ax2[:set_title](L"H vs. $\Delta\mu$, and $\sigma_P$")
    PyPlot.tight_layout()
    #
    savefig("figs/hellinger_surf.pdf",transparent=true,dpi=300)
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
function is_sorted(data::Array)
	sorted = true
	for i in 1 : length(data) - 1
		sorted &= data[i + 1] >= data[i]
	end
	return sorted
end

function insertionsort!(data::Array)
	for i in 2 : length(data)
		current_value = data[i]
		for j in i - 1 : -1 : 1
			if current_value < data[j]
				data[j + 1] = data[j]
				if j == 1
					data[j] = current_value
				end
			else
				data[j + 1] = current_value
				break
			end
		end
	end
end

function bubblesort!(data::Array)
	for i in length(data) - 1 : -1 : 1
		for j in i : length(data) - 1
			if data[j] > data[j + 1]
				data[j], data[j + 1] = data[j + 1], data[j]
			end
		end
	end
end

function builtinsort!(data::Array)
	sort!(data)
end

function quicksort!(data::Array, start=1, stop=-1)
	if stop == -1
		stop = length(data)
	end

	if start == stop
		return

	p = rand(start:stop, 1)

	for i in start:stop
		if i < p
			if data[i] > data[p]
				data[i], data[p] = data[p], data[i]
				p = i
			end
		elseif i > p
			if data[i] < data[p]
				data[i], data[p] = data[p], data[i]
				p = i
			end
		end
	end
	quicksort!(data, start, p - 1)
	quicksort!(data, p + 1, stop)
end

function quicksort_dp!(data::Array, start=-1, stop=-1, lowest=-1, highest=-1)
	if start == -1
		start = 1
		stop = length(data)
		pivot = 0
		lowest = data[1]
		highest = data[1]
		for i in 2 : stop
			if data[i] < lowest
				lowest = data[i]
			end
			if data[i] > highest
				highest = data[i]
			end
		end
	end

	pivot = (highest + lowest) / 2

#=	print(start)
	print(", ")
	print(stop)
	print(", ")
	print(pivot)
	print(", ")
	print(lowest)
	print(", ")
	println(highest)
=#
	if highest == lowest
		return
	end

	if stop - start <= 2
		if stop - start == 2
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
			if data[start + 1] > data[start + 2]
				buffer = data[start + 1]
				data[start + 1] = data[start + 2]
				data[start + 2] = buffer
			end
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
#			println(data)
			return
		end

		if stop - start == 1
			if data[start] > data[start + 1]
				buffer = data[start]
				data[start] = data[start + 1]
				data[start + 1] = buffer
			end
#			println(data)
			return
		end

		if stop - start == 0
			return
		end
	end

	leftPointer = start
	rightPointer = stop

	leftHighest = lowest
	rightLowest = highest

	done = false

	while !done
#=		print(leftPointer)
		print(", ")
		println(rightPointer)
=#
		checkForLeftExtreme = false
		checkForRightExtreme = false

		left = data[leftPointer] <= pivot
		right = data[rightPointer] > pivot

		if rightPointer - leftPointer <= 0
			done = true
		elseif !left && !right
			data[leftPointer], data[rightPointer] = data[rightPointer], data[leftPointer]
			if rightPointer - leftPointer > 1
				checkForLeftExtreme = true
				checkForRightExtreme = true
				leftPointer += 1
				rightPointer -= 1
			else
				done = true
				if data[leftPointer] > leftHighest
					leftHighest = data[leftPointer]
				end
				if data[rightPointer] < rightLowest
					rightLowest = data[rightPointer]
				end
			end
		elseif !left && right
			rightPointer -= 1
			checkForRightExtreme = true
		elseif left && !right
			leftPointer += 1
			checkForLeftExtreme = true
		elseif rightPointer - leftPointer > 1
			leftPointer += 1
			rightPointer -= 1
			checkForLeftExtreme = true
			checkForRightExtreme = true
		else
			done = true
		end

		if checkForLeftExtreme
			if data[leftPointer - 1] > leftHighest
				leftHighest = data[leftPointer - 1]
			end
		end
		if checkForRightExtreme
			if data[rightPointer + 1] < rightLowest
				rightLowest = data[rightPointer + 1]
			end
		end

#		println(data)
	end

	pivotPointer = leftPointer

	if leftPointer == rightPointer
		if data[leftPointer] > pivot
			pivotPointer -= 1
			if data[leftPointer] < rightLowest
				rightLowest = data[leftPointer]
			end
		else
			if data[leftPointer] > leftHighest
				leftHighest = data[leftPointer]
			end
		end
	else
		if data[leftPointer] > leftHighest
			leftHighest = data[leftPointer]
		end
		if data[rightPointer] < rightLowest
			rightLowest = data[rightPointer]
		end
	end

	quicksort!(data, start, pivotPointer, lowest, leftHighest)
	quicksort!(data, pivotPointer + 1, stop, rightLowest, highest)
end


text:
using GNSSTools
using Statistics
using ProgressMeter
using JLD
using PyPlot
pygui(true)

h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
CN0 = 45.  # dB⋅Hz
f_d = 0
nADC = 16
code_start_idx = 1000
iterations = 100
Ts = Array(1e-3:1e-3:100e-3)
# t_length = maximum(Ts)

# L1 C/A
f_s = 2.048e6
sigtype = define_l1ca_code_type(; sig_freq=1*L1_freq)
# generatesignal!(signal)
l1ca_SNRs = zeros(length(Ts))
p = Progress(length(Ts)*iterations, 1, "Processing L1 C/A...")
for i in 1:length(Ts)
    SNR = 0
    T = Ts[i]
    signal = definesignal(sigtype, f_s, T;
                          prn=1, code_start_idx=code_start_idx, f_d=f_d, fd_rate=0,
                          CN0=CN0, include_phase_noise=false, f_if=1.25e6,
                          include_thermal_noise=true, include_adc=false, nADC=nADC,
                          include_carrier=true, receiver_h_parms=h_parms,
                          skip_noise_generation=false)
    replica = definesignal(sigtype, f_s, T; f_d=f_d, allocate_noise_vectors=false,
                           skip_noise_generation=true)
    generatereplica!(replica)
    for j in 1:iterations
        # definesignal!(signal; new_thermal_noise=true, new_phase_noise=false)
        # generatesignal!(signal)
        generatereplica!(signal)
        A = sqrt(2*GNSSTools.k*535*10^(CN0/10))
        σₙ = sqrt(GNSSTools.k*2.046e6*535)
        signal.data .= A.*signal.data .+ σₙ.*randn(signal.sample_num)
        corr_result = fft_correlate(signal.data[1:replica.sample_num], replica.data)
        corr_result_fft = fft(corr_result)
        PS = 2*maximum(abs2.(corr_result_fft))
        # PS = abs2(corr_result[code_start_idx])
        PN = (sum(abs2.(corr_result_fft)) - PS) / (length(corr_result_fft) - 2)
        SNR += PS/PN
        next!(p)
    end
    l1ca_SNRs[i] = 10log10(SNR/iterations)
end

fig = figure()
plot(Ts.*1000, l1ca_SNRs, label="L1 C/A")
xlabel("Integration Times (ms)")
ylabel("Acquired SNR (dB)")
legend()
xscale("log")

save("ch2_snr_vs_integration_time.jld", "l1ca_SNRs", l1ca_SNRs)
savefig("figures/ch2_snr_vs_integration_time.pdf", dpi=300)


text:
import Base

# Base.show(io::IO, data::Data) = print(io, data.label)

Base.show(io::IO, data::Data) = print(io, _show(data)...)
Base.show(io::IO, cascade::Cascade{Data}) = print(io, _show(cascade; maxlen=3)...)

function _show(cascade::Cascade{T}; kwargs...) where T<:Sampling
    fields = fieldnames(Cascade)
    return vcat(
        Printf.@sprintf("%s\n",typeof(cascade)),
        indent.(vcat([_show(cascade, field; kwargs...) for field in fields]...),1),
    )[1:end-1]
end

function _show(plot::Cascade{T}; kwargs...) where T<:Sampling
    fields = fieldnames(Cascade)
    return vcat(
        Printf.@sprintf("%s\n",typeof(cascade)),
        indent.(vcat([_show(cascade, field; kwargs...) for field in fields]...),1),
    )[1:end-1]
end

function _show(data::Data; show_type=true, kwargs...)
    # fields = collect(fieldnames(Data))
    fields = setdiff(collect(fieldnames(Data)),[:label])
    str = vcat(
        Printf.@sprintf("%s\"%s\"\n", (show_type ? "Data: " : ""), data.label),
        [Printf.@sprintf("%s%-11s %s\n", indent(), "$field:", _show(getproperty(data,field)))
            for field in fields]
    )
    return str
end

function _show(axis::Axis; kwargs...)
    fields = setdiff(collect(fieldnames(Axis)),[:label])
    str = vcat(
        Printf.@sprintf("\"%s\"\n", axis.label),
        [Printf.@sprintf("%s%-11s %s\n", indent(), "$field:", _show(getproperty(axis,field)))
            for field in fields]
    )
    return str
end

function _show(cascade::Cascade, field::Symbol; kwargs...)
    return vcat(_show(field, getproperty(cascade,field); kwargs...), "\n")
end

# Should be able to merge these two.....
function _show(field::Symbol, property::Vector; kwargs...)
    return vcat(
        Printf.@sprintf("%s: %d-element %s\n", field, length(property), typeof(property)),
        indent.(_show(property; kwargs...), 1),
    )
end

function _show(field::Symbol, property; kwargs...)
    return vcat(
        Printf.@sprintf("%s: %s\n", field, typeof(property)),
        indent.(_show(property; kwargs...), 1),
    )
end

function _show(lst::Vector{Data}; maxlen=missing)
    N = length(lst)
    if N==1
        str = _show(lst[1])
    else
        !ismissing(maxlen) && (maxlen=max(3,maxlen))
        shorten = !ismissing(maxlen) && N>maxlen
        N1 = shorten ? maxlen-2 : N-1

        idx = [1:N1;N]
        str = itemize.(_show.(lst[idx]), idx)

        str = vcat(str[1:N1]..., " ⋮\n"^shorten, str[end])
        str = str[.!isempty.(str)]
    end
    return str
end

function _show(mat::Matrix; kwargs...)
    return if size(mat,2)==1
        _show(mat[:,1]; kwargs...)
    end
end

function _show(lst::Vector; maxlen=5)
    all(typeof.(pdata.xaxis.ticklabels).==String) && (maxlen=4)
    if length(lst)==1
        str = _show(lst[1])
    else
        shorten = length(lst)>maxlen

        N1 = shorten ? maxlen-2 : length(lst)-1

        str1 = string(string.(_show.(lst[1:N1]),", ")...)
        str2 = _show(lst[end])

        str = Printf.@sprintf("%s%s%s", str1, (shorten ? "… " : ""), str2)
    end
    return str
end

_show(x::String) = Printf.@sprintf("\"%s\"", x)
_show(x::Float64) = Printf.@sprintf("%.5g", x)
_show(x::Integer) = Printf.@sprintf("%d", x)
_show(x::Tuple) = _show(collect(x))

_show(x::Symbol, T::DataType) = Printf.@sprintf("%s: %s", x, T)

function indent(str::String, args...; kwargs...)
    if str!=="\n"
        m = match.(r"^(\n)?(.*\n?)", str)
        str = string((isnothing(m[1]) ? "" : m[1]), indent(args...; kwargs...), m[2])
    end
    return str
end

indent(level::Integer=1; spaces::Integer=2) = " "^(spaces*level)

function itemize!(lst::Vector{String}, ii::Integer)
    lst[1] = Printf.@sprintf("%2d) %s", ii, lst[1])
    lst[2:end] .= indent.(lst[2:end], 2)
    return lst
end

itemize(lst, ii) = itemize!(copy(lst),ii)

# print(_show(data; maxlen=4)...)


text:
#= Custom currencies =#

@testset "Custom" begin

@usingcustomcurrency xbt "Bitcoin (100 satoshi unit)" 2
custom = newcurrency!(:custom, "Custom Currency", 6)

@testset "@usingcustomcurrency" begin
    @test occursin("xbt", stringmime("text/plain", 10xbt))
    @test occursin("10.00", stringmime("text/plain", 10xbt))
    @test 10xbt - 5xbt == 5xbt
    @test format(-1111.11xbt, styles=[:us, :finance]) == "xbt (1,111.11)"
    @test Basket([10xbt, 10USD]) - 10USD == 10xbt

    @test currencyinfo(Monetary{:xbt}) == "Bitcoin (100 satoshi unit)"

    @test iso4217alpha(Monetary{:xbt}) == "xbt"
    @test iso4217alpha(100xbt) == "xbt"

    @test currency(10xbt) == :xbt
end

@testset "newcurrency!()" begin
    @test stringmime("text/plain", 20custom) == "20.000000 custom"
    @test string(20custom) == "20.0custom"
    @test format(-20custom) == "(20.000000) custom"
    @test 10custom / 10000000 == 0.000001custom

    @test currencyinfo(:custom) == "Custom Currency"
    @test currencyinfo(10custom) == "Custom Currency"

    @test iso4217num(:custom) == 0
    @test iso4217num(10custom) == 0

    @test currency(custom) == :custom
end

end  # testset Custom


text:
function banded_chol!(::BandedColumns{DenseColumnMajor},
                       A::AbstractMatrix{T}, ::Type{UpperTriangular}) where T<:BlasFloat
    _, info = pbtrf!('U', size(A, 1), bandwidth(A,2), bandeddata(A))
    UpperTriangular(A), info
end

function banded_chol!(::BandedColumns{DenseColumnMajor},
                      A::AbstractMatrix{T}, ::Type{LowerTriangular}) where T<:BlasFloat
    _, info = pbtrf!('L', size(A, 1), bandwidth(A,1), bandeddata(A))
    LowerTriangular(A), info
end

## Non BLAS/LAPACK element types (generic)
function banded_chol!(_, A::AbstractMatrix, ::Type{UpperTriangular})
    require_one_based_indexing(A)
    n = checksquare(A)
    u = bandwidth(A,2)
    @inbounds begin
        for k = 1:n
            for i = max(1,k-u):k - 1
                A[k,k] -= A[i,k]'A[i,k]
            end
            Akk, info = _chol!(A[k,k], UpperTriangular)
            if info != 0
                return UpperTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(copy(Akk'))
            for j = k + 1:min(k+u,n)
                for i = max(1,j-u):k - 1
                    A[k,j] -= A[i,k]'A[i,j]
                end
                A[k,j] = AkkInv*A[k,j]
            end
        end
    end
    return UpperTriangular(A), convert(BlasInt, 0)
end
function banded_chol!(_, A::AbstractMatrix, ::Type{LowerTriangular})
    require_one_based_indexing(A)
    n = checksquare(A)
    l = bandwidth(A,1)
    @inbounds begin
        for k = 1:n
            for i = max(1,k-l):k - 1
                A[k,k] -= A[k,i]*A[k,i]'
            end
            Akk, info = _chol!(A[k,k], LowerTriangular)
            if info != 0
                return LowerTriangular(A), info
            end
            A[k,k] = Akk
            AkkInv = inv(Akk)
            for j = max(1,k-l):k - 1
                @simd for i = k + 1:min(n,j+l)
                    A[i,k] -= A[i,j]*A[k,j]'
                end
            end
            for i = k + 1:min(n,k+l)
                A[i,k] *= AkkInv'
            end
        end
     end
    return LowerTriangular(A), convert(BlasInt, 0)
end

banded_chol!(A, ::Type{T}) where T = banded_chol!(MemoryLayout(typeof(A)), A, T)
_chol!(A::AbstractBandedMatrix, ::Type{UpperTriangular}) = banded_chol!(A, UpperTriangular)
_chol!(A::AbstractBandedMatrix, ::Type{LowerTriangular}) = banded_chol!(A, LowerTriangular)

function _ldiv!(::BandedColumns{DenseColumnMajor}, ::AbstractStridedLayout,
                A::Cholesky{T}, B::AbstractVecOrMat{T}) where T<:BlasFloat
    if A.uplo == 'U'
        pbtrs!('U', size(A, 1), bandwidth(A.factors,2), bandeddata(A.factors), B)
    else
        pbtrs!('L', size(A, 1), bandwidth(A.factors,1), bandeddata(A.factors), B)
    end
    B
end

function _ldiv!(::AbstractBandedLayout, _, C::Cholesky, B::AbstractMatrix)
    if C.uplo == 'L'
        return ldiv!(adjoint(LowerTriangular(C.factors)), ldiv!(LowerTriangular(C.factors), B))
    else
        return ldiv!(UpperTriangular(C.factors), ldiv!(adjoint(UpperTriangular(C.factors)), B))
    end
end

ldiv!(A::Cholesky{T,<:AbstractBandedMatrix}, B::StridedVecOrMat{T}) where T<:BlasFloat =
    _ldiv!(MemoryLayout(typeof(A.factors)), MemoryLayout(typeof(B)), A, B)

# For some bizarre reason this isnt in LinearAlgebra
ArrayLayouts._cholesky(::Union{SymmetricLayout{<:AbstractBandedLayout},HermitianLayout{<:AbstractBandedLayout}}, ax, A, ::Val{false}=Val(false); check::Bool = true) = cholesky!(cholcopy(A); check = check)

# cholesky(A::Hermitian{T,<:BandedMatrix{T}},
#     ::Val{false}=Val(false); check::Bool = true) where T = cholesky!(cholcopy(A); check = check)


text:
@testset "autodiff" begin

    @testset "fill" begin
        tp = Gen.new_tape()
        x = Gen.track(2.3, tp)
        y = fill(x, 3)
        grad = rand(3)
        Gen.deriv!(y, grad)
        Gen.reverse_pass!(tp)
        @test isapprox(Gen.deriv(x), sum(grad))

        tp = Gen.new_tape()
        x = Gen.track(2.3, tp)
        y = fill(x, 3, 2)
        grad = rand(3, 2)
        Gen.deriv!(y, grad)
        Gen.reverse_pass!(tp)
        @test isapprox(Gen.deriv(x), sum(grad))
    end

end


text:
immutable MersenneFunctor{T} <: Functor{1}
    mt::MersenneTwister
end
@compat( @inline (rf::MersenneFunctor{T}){T}(i...) = rand(rf.mt, T))

immutable RandFunctor{T} <: Functor{1}
    valuerange::T
end
@compat( @inline (rf::RandFunctor)(i...) = rand(rf.valuerange) )
immutable RandnFunctor{T} <: Functor{1}
    mt::MersenneTwister
end
@compat( @inline (rf::RandnFunctor{Float64})(i...) = randn(rf.mt) )
@compat( @inline (rf::RandnFunctor{Complex{Float64}})(i...) = randn(rf.mt) + im*randn(rf.mt) )

immutable ConstFunctor{T} <: Functor{1}
    args::T
end
@compat (f::ConstFunctor)(i...) = f.args

immutable EyeFunc{T} <: Functor{1} end
@compat( @inline (::Type{EyeFunc{T}}){T}(i::Int, j::Int) = (i==j ? one(T) : zero(T)) )

immutable UnitFunctor <: Functor{1}
    i::Int
    eltype::DataType
end
@compat (ef::UnitFunctor)(i) = ef.i==i ? one(ef.eltype) : zero(ef.eltype)

immutable ConversionIndexFunctor{T, T2} <: Functor{1}
    args1::T
    target::Type{T2}
end
@compat (f::ConversionIndexFunctor)(i...) = f.target(f.args1[i...])

immutable IndexFunctor{T} <: Functor{1}
    args1::T
end
@compat (f::IndexFunctor)(i) = f.args1[i]

immutable IndexFunc{T} <: Functor{1}
    arg::T
    i::Int
end
@compat (a::IndexFunc{T}){T}(j) = a.arg[a.i,j]

immutable CIndexFunc{T} <: Functor{1}
    arg::T
    i::Int
end
@compat (a::CIndexFunc{T}){T}(j) = a.arg[a.i,j]'

immutable CRowFunctor{M}
    mat::M
end
@compat (r::CRowFunctor)(i::Int) = crow(r.mat, i)

immutable SetindexFunctor{T <: FixedArray, V, N} <: Functor{1}
    target::T
    value::V
    index::NTuple{N, Int}
end

@compat function (sf::SetindexFunctor)(i::Int...)
    sf.index == i && return eltype(sf.target)(sf.value)
    sf.target[i...]
end

immutable RowFunctor{M}
    mat::M
end
@compat (r::RowFunctor)(i::Int) = row(r.mat, i)


text:
using LogClustering.KATE
using LogClustering.KATE: KCompetetive
using Base.Test
using Flux
using Flux: throttle
using Flux: crossentropy #, binarycrossentropy
using Juno: @progress

cd(@__DIR__)

isfile("input.txt") ||
  download("http://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt",
           "input.txt")

input, wc = KATE.transform_text_to_input("input.txt", 1000)
vocabulary = collect(keys(wc))

# Xs, Ys = deepcopy(input), deepcopy(input)
# Xs, Ys = deepcopy([input]), deepcopy([input])
Xs, Ys = [deepcopy(input)], [deepcopy(input)]
@show typeof(Xs), typeof(Ys)
@show length(Xs), length(Ys), length(Xs[1]), length(Ys[1])
# @show Xs, Ys

N = length(input)

k = 64

m = Chain(
    KCompetetive(N, k, tanh),
    Dense(k, N, sigmoid)
)

function loss(xs, ys)
    # @show typeof(xs), typeof(ys)
    l = crossentropy(m(xs), ys)
    # @show typeof(l), l
    return l
end

function evaluation_callback()
    @show loss(Xs[1], Ys[1])
end

opt = Flux.ADADelta(params(m))

@progress for i = 1:10
    info("Epoch $i")
    Flux.train!(
        loss,
        zip(Xs, Ys),
        opt,
        cb = throttle(evaluation_callback, 30))
end

# @show KATE.get_similar_words(m, 1, vocabulary)


text:
include("importance_sampling.jl")
include("train.jl")
include("variational.jl")
include("particle_filter.jl")
include("mala.jl")
include("hmc.jl")
include("map_optimize.jl")
include("elliptical_slice.jl")
include("mh.jl")
include("trace_translators.jl")


text:
@testset "Colormaps" begin
    @test length(colormap("RdBu", 100)) == 100
end


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = [1.0, 1.0, 1.0, 1.0];
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    gam = mu[1];
    gam1 = gam - 1.0;
    r = u[1];
    ru = u[2];
    rv = u[3];
    rE = u[4];
    r1 = 1/r;
    uv = ru*r1;
    vv = rv*r1;
    E = rE*r1;
    p = gam1*(rE-r*0.5*(uv*uv+vv*vv));
    h = E+p*r1;
    f = [ru, ru*uv+p, rv*uv, ru*h, rv, ru*vv, rv*vv+p, rv*h];
    f = reshape(f,(4,2));
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = [0.0, 0.0, 0.0, 0.0];
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = zeros((4,2));
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);

    gam = mu[1];
    gam1 = gam - 1.0;
    r = u[1];
    ru = u[2];
    rv = u[3];
    rE = u[4];
    nx = n[1];
    ny = n[2];

    r1 = 1/r;
    uv = ru*r1;
    vv = rv*r1;
    E = rE*r1;
    p = gam1*(rE-r*0.5*(uv*uv+vv*vv));
    h = E+p*r1;
    a = sqrt(gam*p*r1);

    run = ru*nx + rv*ny;
    rut = -ru*ny + rv*nx;
    un = run/r;
    ut = rut/r;

    K = [ 1  1  0  1 ; un-a  un  0  un+a ; ut  ut  1  ut ; h - un*a  (1/2)*(un^2 + ut^2)  ut  h+un*a];
    Kinv = (gam1/(2*a^2))*[ h + (a/gam1)*(un-a)  -(un+a/gam1)  -ut  1; -2*h + (4/gam1)*a^2  2*un  2*ut  -2 ; -2*(ut*a^2)/gam1  0  2*(a^2)/gam1  0 ; h - a*(un+a)/gam1  -un+a/gam1  -ut  1];
    T = [ 1 0 0 0; 0 nx ny 0; 0 -ny nx 0; 0 0 0 1];
    Tinv = [ 1 0 0 0; 0 nx -ny 0; 0 ny nx 0; 0 0 0 1];
    Lambda = [ tanh(1e3*(un-a)) 0 0 0 ; 0 tanh(1e3*(un)) 0 0 ; 0 0 tanh(1e3*(un)) 0 ; 0 0 0 tanh(1e3*(un+a)) ];
    # E = simplify(K * Lambda * Kinv);
    # An = simplify(Tinv * E * T);
    E = (K * Lambda * Kinv);
    An = (Tinv * E * T);

    # freestream boundary condition
    uinf = mu[3:6]; # freestream flow
    uinf = uinf[:];
    u = u[:];          # state variables
    ui = 0.5*((u.+uinf) + An*(u.-uinf));  # Riemann solution
    fi = f[:,1]*n[1] + f[:,2]*n[2] + tau[1]*(u-ui);

    uw = u[:];
    uw[2] = ru - run*nx;
    uw[3] = rv - run*ny;
    uw = 0.5*((u.+uw) + An*(u.-uw));  # Riemann solution
    fw = f[:,1]*n[1] + f[:,2]*n[2] + tau[1]*(u-uw); # numerical flux at freestream boundary
    fw[1] = 0;
    fw[end] = 0;
    fb = hcat(fw[:], fi[:]); # wall and freestream boundary conditions

    return fb;
end
function initu(x, mu, eta)

    u0 = mu[3:6]; # freestream flow
    return u0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
"""シミュレーション"""

using LinearAlgebra  # 逆行列とノルムの計算に使用
using Plots  # グラフ作成
#using ProgressBars
#using Parameters
using DifferentialEquations  # 微分方程式のソルバ
using LaTeXStrings

include("kinematics.jl")
include("dynamics.jl")
include("controller.jl")
include("plot_utils.jl")

using .Kinematics: Phi0, Arm
using .Dynamics: H_dot, calc_q_dot_dot, uncertain_K, uncertain_D
using .Controller

const R = 0.015

"""適当な目標アクチュエータ変位"""
function calc_qd(t::T) where T
    [
        R * sin(3t),
        R * sin(3t + π/2),
        R * sin(3t + π),
    ]
end

"""適当な目標アクチュエータ速度"""
function calc_qd_dot(t::T) where T
    [
        R * 3cos(3t),
        R * 3cos(3t + π/2),
        R * 3cos(3t + π),
    ]
end

"""適当な目標アクチュエータ加速度"""
function calc_qd_dot_dot(t::T) where T
    [
        R * -9sin(3t),
        R * -9sin(3t + π/2),
        R * -9sin(3t + π),
    ]
end

"""ソルバで使うやつ"""
function state_eq!(
    X_dot::Vector{T}, X::Vector{T},
    p,
    t::T
    ) where T
    q = X[1:3]
    q_dot = X[4:6]
    H = X[7:9]
    #println(t)
    τ = calc_torque(
    p, q, q_dot,
    calc_qd(t), calc_qd_dot(t), calc_qd_dot_dot(t),
    )
    X_dot[1:3] = q_dot
    X_dot[4:6] = calc_q_dot_dot(τ, q, q_dot, H)
    X_dot[7:9] = H_dot(H, q, q_dot)
end

"""ソルバで使うやつ

受動性に基づく適応制御
"""
function state_eq!(
    X_dot::Vector{T}, X::Vector{T},
    p::PassivityBasedAdaptiveController{T},
    t::T
    ) where T
    q = X[1:3]
    q_dot = X[4:6]
    H = X[7:9]
    θp = X[10:15]  # パラメータの推定値

    #println(t)
    #println(θp)
    τ = calc_torque(
    p, q, q_dot,
    calc_qd(t), calc_qd_dot(t),
    θp
    )
    X_dot[1:3] = q_dot
    X_dot[4:6] = calc_q_dot_dot(τ, q, q_dot, H)
    X_dot[7:9] = H_dot(H, q, q_dot)
    X_dot[10:15] = θp_dot(p, q, q_dot, calc_qd(t), calc_qd_dot(t))
end

"""ソルバで使うやつ"""
function state_eq!(
    X_dot::Vector{T}, X::Vector{T},
    p::MPCController,
    t::T
    ) where T
    q = X[1:3]
    q_dot = X[4:6]
    H = X[7:9]
    #println(t)
    τ = calc_torque(
    p, q, q_dot, H,
    calc_qd, calc_qd_dot, t,
    )
    X_dot[1:3] = q_dot
    X_dot[4:6] = calc_q_dot_dot(τ, q, q_dot, H)
    X_dot[7:9] = H_dot(H, q, q_dot)
end

"""微分方程式を解く"""
function run_simulation(
    TIME_SPAN::T, method_name::String, controller_param
    ) where T

    X₀ = zeros(T, 9)
    t_span = (0.0, TIME_SPAN)
    println(method_name * " now...")
    prob = ODEProblem(state_eq!, X₀, t_span, controller_param)
    solve(prob)
end

"""微分方程式を解く"""
function run_simulation(
    TIME_SPAN::T, method_name::String,
    controller_param::PassivityBasedAdaptiveController{T}
    ) where T

    X₀ = [
        zeros(T, 9); diag(uncertain_K); diag(uncertain_D)
    ]
    t_span = (0.0, TIME_SPAN)
    println(method_name * " now...")
    prob = ODEProblem(state_eq!, X₀, t_span, controller_param)

    solve(prob)
end

# """微分方程式を解く"""
# function run_simulation(
#     TIME_SPAN::T, method_name::String,
#     controller_param::MPCController{T}
#     ) where T

#     X₀ = zeros(T, 9)
#     t_span = (0.0, TIME_SPAN)
#     println(method_name * " now...")
#     prob = ODEProblem(state_eq!, X₀, t_span, controller_param)

#     solve(prob)
# end

function reproduce_τ(p::PassivityBasedAdaptiveController{T}, t::T, u::Vector{T}) where T
    calc_torque(
        p, u[1:3], u[4:6],
        calc_qd(t), calc_qd_dot(t),
        u[10:15]
    )
end

function reproduce_τ(p, t::T, u::Vector{T}) where T
    calc_torque(
        p, u[1:3], u[4:6],
        calc_qd(t), calc_qd_dot(t), calc_qd_dot_dot(t),
    )
end

"""実行"""
function exmample()

    TIME_SPAN = 0.05

    hutashikasa = false  # 剛性行列と減衰行列に不確かさがあるかないか

    params = [
        # (
        #     name = "kine",
        #     p = KinematicController(
        #         K_kin = Dynamics.K,
        #         isUncertainty = hutashikasa,
        #     ),
        #     color = :red
        # ),
        # (
        #     name = "pdfb",
        #     p = PDandFBController(
        #         Kd = Matrix{Float64}(I, 3, 3)*200,
        #         Kp = Matrix{Float64}(I, 3, 3)*10000,
        #         isUncertainty = hutashikasa,
        #     ),
        #     color = :blue,
        # ),
        # (
        #     name = "psiv",
        #     p = PassivityBasedController(
        #         Λ = Matrix{Float64}(I, 3, 3)*10,
        #         KG = Matrix{Float64}(I, 3, 3)*10,
        #         isUncertainty = hutashikasa,
        #     ),
        #     color = :magenta
        # ),
        # (
        #     name = "psad",
        #     p = PassivityBasedAdaptiveController(
        #         invΓ = Matrix{Float64}(I, 6, 6)*1e+5,
        #         Λ = Matrix{Float64}(I, 3, 3)*10,
        #         KG = Matrix{Float64}(I, 3, 3)*10,
        #         isUncertainty = hutashikasa,
        #     ),
        #     color = :black
        # ),
        (
            name = "sdre",
            p = SDREController(
                Q = Matrix{Float64}(I, 6, 6)*2000,
                R = Matrix{Float64}(I, 3, 3)*0.1,
                isUncertainty = hutashikasa,
            ),
            color = :green
        ),
        # (
        #     name = "dmpc",
        #     p = MPCController(
        #         Q = diagm([1, 1, 1, 1, 1, 1, 0.0, 0.0, 0.0]),
        #         R = diagm([1.0, 1.0, 1.0,]),
        #         n = 10,
        #         Δt = 1.0e-15,
        #         isUncertainty = hutashikasa
        #     ),
        #     color = :cyan
        # )
    ]
    sols = []
    fig0 = plot(xlims=(0.0, TIME_SPAN), ylims=(0,0.02))
    fig_tau = plot(xlims=(0.0, TIME_SPAN), ylims=(-150.0, 150.0))
    fig1 = plot(xlims=(0.0, TIME_SPAN))
    fig2 = plot(xlims=(0.0, TIME_SPAN))
    fig3 = plot(xlims=(0.0, TIME_SPAN))

    for param in params
        @time sol = run_simulation(TIME_SPAN, param.name, param.p)
        push!(sols, sol)

        error = Vector{Vector{Float64}}(undef, length(sol.t))
        τ = Vector{Vector{Float64}}(undef, length(sol.t))
        for (i, t) in enumerate(sol.t)
            error[i] = calc_qd(t) .- sol.u[i][1:3]
            τ[i] = reproduce_τ(param.p, t, sol.u[i])
        end
        L2_error = norm.(error)
        plot!(
            fig0,
            sol.t, L2_error,
            label=(param.name * "-") * "e",
            legend=:outerright,
            c=param.color,
            ylabel="error [m]"
        )

        x, y, z = split_vec_of_arrays(τ)
        plot!(fig_tau, sol.t, x, label=param.name*"-"*"1", c=param.color, linestyl=:solid)
        plot!(fig_tau, sol.t, y, label=param.name*"-"*"2", c=param.color, linestyle=:dash)
        plot!(fig_tau, sol.t, z, label=param.name*"-"*"3", c=param.color, linestyle=:dot)
        plot!(
            fig_tau, legend=:outerright,
            ylabel="τ [N/rad]"
        )

        plot!(
            fig1,
            sol, vars=[(0,1), (0,2), (0,3)],
            label=(param.name * "-") .* ["1" "2" "3"],
            legend=:outerright,
            c=param.color,
            linestyle=[:solid :dash :dot],
            ylabel="l [m]"
        )
        plot!(
            fig2,
            sol, vars=[(0,4), (0,5), (0,6)],
            label=(param.name * "-") .* ["1" "2" "3"],
            legend=:outerright,
            c=param.color,
            linestyle=[:solid :dash :dot],
            ylabel="l_dot  [m/s]"
        )
        plot!(
            fig3,
            sol, vars=[(0,7), (0,8), (0,9)],
            label=(param.name * "-") .* ["1" "2" "3"],
            legend=:outerright,
            c=param.color,
            linestyle=[:solid :dash :dot],
            ylabel="h  [N/rad]"
        )

    end

    fig_I = plot(
        fig0, fig_tau, fig1, fig2, fig3,
        layout=(5,1), size=(800, 1600),
        margin = 1.2Plots.cm,
    )
    savefig(fig_I, "solve.png")

    # fig_param_K = plot(xlims=(0.0, TIME_SPAN))
    # plot!(
    #     fig_param_K,
    #     sols[2], vars=[(0,10), (0,11), (0,12)],
    #     label=["K1" "K2" "K3"],
    #     legend=:outerright,
    # )
    # # plot!(
    # #     fig_param_K,
    # #     [sols[2].t[1], sols[2].t[end]], [Dynamics.K[1], Dynamics.K[1],],
    # #     label="real_K"
    # # )

    # fig_param_D = plot(xlims=(0.0, TIME_SPAN))
    # plot!(
    #     fig_param_D,
    #     sols[2], vars=[(0,13), (0,14), (0,15)],
    #     label=["D1" "D2" "D3"],
    #     legend=:outerright,
    # )
    # # plot!(
    # #     fig_param_D,
    # #     [sols[2].t[1], sols[2].t[end]], [Dynamics.D[1], Dynamics.D[1],],
    # #     label="real_D"
    # # )
    # fig_hog = plot(
    #     fig_param_K, fig_param_D,
    #     layout=(2,1),
    # )
    # savefig(fig_hog, "param.png")

    sols
    println("done!")
end

@time _ = exmample()


text:
using KrylovKit, ArnoldiMethod, LinearMaps, LinearAlgebra, RecursiveArrayTools

struct PrecPartialSchur{Ts, Tu, Tsm1, Teigen}
	S::Ts
	U::Tu
	Sm1::Tsm1
	eigenvalues::Teigen
end

function LinearAlgebra.ldiv!(out, Pl::PrecPartialSchur, rhs::AbstractArray)
	########################################################
	# U * Sm1 * Ut + (I - U Ut)
	# U(Sm1 - I)Ut + I
	########################################################
	# y = (Pl.U * Pl.Sm1 * transpose(Pl.U) .+ (I - Pl.U * transpose(Pl.U))) * rhs
	y = transpose(Pl.U) * rhs
	y .= Pl.Sm1 * y .- y
	out .= rhs .+ Pl.U * y
end

function LinearAlgebra.ldiv!(Pl::PrecPartialSchur, rhs::AbstractArray)
	########################################################
	# U * Sm1 * Ut + (I - U Ut)
	# U(Sm1 - I)Ut + I
	########################################################
	# y = (Pl.U * Pl.Sm1 * transpose(Pl.U) .+ (I - Pl.U * transpose(Pl.U))) * rhs
	y = transpose(Pl.U) * rhs
	y .= Pl.Sm1 * y .- y
	rhs .= rhs .+ Pl.U * y
end

function LinearAlgebra.:\(Pl::PrecPartialSchur, rhs::AbstractArray)
	out = similar(rhs)
	ldiv!(out, Pl, rhs)
	return out
end

"""
	PrecPartialSchurKrylovKit(J, x0, nev, which = :LM; krylovdim = max(2nev, 20), verbosity = 0)

Builds a preconditioner based on deflation of `nev` eigenvalues chosen according to `which`. A partial Schur decomposition is computed (Matrix-Free), using the package `KrylovKit.jl`, from which a projection is built. The options are similar to the ones of `EigKrylovKit()`.
"""
function PrecPartialSchurKrylovKit(J, x0, nev, which = :LM; krylovdim = max(2nev, 20), verbosity = 0)
	H, V, vals, info = KrylovKit.schursolve(J, x0, nev, which, Arnoldi(krylovdim = krylovdim, verbosity = verbosity))
	Q, S = qr(H)
	U = VectorOfArray(V) * Q
	return PrecPartialSchur(S, U, inv(S), vals)
end

"""
	PrecPartialSchurArnoldiMethod(J, N, nev, which = LM(); tol = 1e-9, kwargs...)

Builds a preconditioner based on deflation of `nev` eigenvalues chosen according to `which`. A partial Schur decomposition is computed (Matrix-Free), using the package `ArnoldiMethod.jl`, from which a projection is built. See the package `ArnoldiMethod.jl` for how to pass the proper options.
"""
function PrecPartialSchurArnoldiMethod(J, N, nev, which = LM(); tol = 1e-9, kwargs...)
	if J isa AbstractMatrix
		decomp, history = ArnoldiMethod.partialschur(J; nev = nev, tol = tol, which = which, kwargs...)
	else
		Jmap = LinearMap{Float64}(J, N, N ; ismutating = false)
		decomp, history = ArnoldiMethod.partialschur(Jmap; nev = nev, tol = tol, which = which, kwargs...)
	end
	return PrecPartialSchur(decomp.R, decomp.Q, inv(decomp.R), decomp.eigenvalues)
end

PrecPartialSchurArnoldiMethod(J::AbstractArray, nev, which = LM(); tol = 1e-9, kwargs...) = PrecPartialSchurArnoldiMethod(J, size(J)[1], nev, which ; tol = tol, kwargs...)


text:
"""
`PSSFSS` is a software package for the analysis of polarization and frequency selective surfaces (PSSs and FSSs).
The user specifies the geometry to be analyzed as a `Vector` containing two or more dielectric [`Layer`](@ref)s
and zero or more [`Sheet`](@ref) objects denoting the PSS/FSS surfaces.  After also specifying the scan angles or
unit cell incremental phasings, frequencies to be considered, and optionally some output parameters to be written
to CSV file(s),
the user then invokes the [`analyze`](@refs) function to perform the analysis.  FSS/PSS triangulations
can be conveniently visualized using the `plot` command of the `Plots` package.

"""
module PSSFSS

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 3
end

using Reexport
using PkgVersion
using InteractiveUtils: versioninfo
using Dates: now
using DelimitedFiles: writedlm
using Printf: @sprintf
using LinearAlgebra: ×, norm, ⋅, factorize, lu!, ldiv!
using StaticArrays: StaticArrays, SVector, SArray, @SVector
using Unitful: ustrip, @u_str
using Logging: with_logger
using ProgressMeter

include("Constants.jl")
include("Log.jl")
include("PSSFSSLen.jl")
include("Rings.jl")
include("Layers.jl")
include("Sheets.jl")
include("Meshsub.jl")
include("Elements.jl")
include("RWG.jl")
include("PGF.jl")
include("Zint.jl")
include("FillZY.jl")
include("GSMs.jl")
include("Modes.jl")
include("Outputs.jl")

using .Rings
@reexport using .Sheets: Sheet, RWGSheet, read_sheet_data
using .RWG: setup_rwg, rwgbfft!, RWGData
using .GSMs: GSM, cascade, cascade!, gsm_electric_gblock, gsm_magnetic_gblock,
             gsm_slab_interface, translate_gsm!, choose_gblocks, Gblock, pecgsm, pmcgsm
using .FillZY: fillz, filly
using .Modes: zhatcross, choose_layer_modes!, setup_modes!
using .Constants: twopi, c₀, tdigits, dbmin
using .Log: pssfss_logger, @logfile
@reexport using .PSSFSSLen
@reexport using .Layers: Layer
@reexport using .Elements: rectstrip, diagstrip, polyring, meander, loadedcross,
                           jerusalemcross, pecsheet, pmcsheet
@reexport using .Outputs: @outputs, extract_result_file, extract_result
using .Outputs: Result, append_result_data

export analyze

Base.isfile(f::Base.DevNull) = false
Base.open(f::Base.DevNull,::AbstractString) = f

"""
    analyze(strata::Vector, flist, steering; outlist=[], logfile="pssfss.log", resultfile="pssfss.res", showprogress::Bool=true)

Analyze a full FSS/PSS structure over a range of frequencies and steering angles/phasings.
Generate output files as specified in `outlist`.

## Positional Arguments
- `strata`:  A vector of `Layer` and `Sheet` objects. The first and last entries must be of type `Layer`.

- `flist`: An iterable containing the analysis frequencies in GHz.

- `steering`: A length 2 `NamedTuple` containing as keys the steering parameter labels and as values
  the iterables that define the values of steering parameters to be analyzed.

    - one of {`:phi` ,`:ϕ`} and one of {`:theta`, `:θ`}, or

    - one of {`:psi1` ,`:ψ₁`} and one of {`:psi2`, `:ψ₂`}.

  All steering parameters are input in degrees.

  The program will analyze while iterating over a triple loop over the two steering
  parameters and frequency, with frequency in the innermost loop (i.e. varying the fastest).
  The steering parameter listed first will be in the outermost loop and will therefore
  vary most slowly.

## Keyword Arguments

- `outlist`:  A matrix with 2 columns.  The first column in each row is a string
  containing the name of the CSV file to write the output to.  The second entry in
  each row is a tuple generated by the `@outputs` macro of the `Outputs` module. The
  contents of the specified file(s) will be updated as the program completes each analysis
  frequency.

- `logfile`:  A string containing the name of the log file to which timing and other
  information about the run is written. Defaults to `"pssfss.log"`.
  If this file already exists, it will be overwritten.

- `resultfile`:  A string containing the name of the results file. Defaults to `pssfss.res`.
  If this file already exists, it will be overwritten.  It is a binary
  file that contains information (including the generalized scattering matrix) from
  the analysis performed for each scan condition and frequency. The result file can be
  post-processed to produce similar or additional outputs that were requested at run time
  via the `outlist` argument.

- `showprogress`: If true (default), then show progress bar during execution.

"""
function analyze(strata::Vector, flist, steering; outlist=[], logfile="pssfss.log",
    resultfile="pssfss.res", showprogress::Bool=true)
    layers = Layer[deepcopy(s) for s in strata if s isa Layer]
    sheets = RWGSheet[s for s in strata if s isa Sheet]
    islayer = map(x -> x isa Layer, strata)
    issheet = map(x -> x isa RWGSheet, strata)
    nl = length(layers)
    nj = nl - 1
    ns = length(sheets)
    sint = cumsum(islayer)[issheet] # sint[k] contains dielectric interface number of k'th sheet
    junc = zeros(Int, nj)
    junc[sint] = 1:ns #  junc[i] is the sheet number present at interface i, or 0 if no sheet is there
    freqstemp = float.(collect(flist))
    if length(freqstemp) < 2
        freqs = Float64[freqstemp]
    else
        freqs::Vector{Float64} = freqstemp
    end

    stkeys::Tuple{Symbol, Symbol} = keys(steering)
    stvaluestemp = [float.(collect(s)) for s in steering]
    stvalues = Vector{Float64}[]
    for (i,s) in pairs(stvaluestemp)
        if length(s) < 2
            push!(stvalues, Float64[s])
        else
            push!(stvalues, s)
        end
    end

    with_logger(pssfss_logger(logfile)) do
        _analyze(layers, sheets, junc, freqs, stkeys, stvalues;
                 outlist, resultfile, showprogress)
    end
end # function

"""
function _analyze(layers, sheets, junc, freqs, stkeys, stvalues;
    outlist=Any[], resultfile="pssfss.res", showprogress::Bool=true)

## Positional Arguments
- `layers`:  A vector of `Layer` objects.

- `sheets`:  A vector of `RWGSheet` objects.

- `junc`: `junc[k]` contains the sheet number present at layer interface `k` or `0` if no sheet
  is present there.

- `freqs`: A vector containing the analysis frequencies in GHz.

- `stkeys`: A length 2 `Tuple` containing as the steering angle labels as `Symbols` Either

    - one of {`:phi` ,`:ϕ`} and one of {`:theta`, `:θ`}, or

    - one of {`:psi1` ,`:ψ₁`} and one of {`:psi2`, `:ψ₂`}.

  The angular steering parameters are input in degrees, while the incremental phase shift
  parameters are input in radians.

  The program will analyze while iterating over a triple loop over the two steering
  parameters and frequency, with frequency in the innermost loop (i.e. varying the fastest).
  The steering parameter listed first will be in the outermost loop and will therefore
  vary most slowly.

## Keyword Arguments

- `outlist`:  A matrix with 2 columns.  The first column in each row is a string
  containing the name of the CSV file to write the output to.  The second entry in
  each row is a tuple generated by the `@outputs` macro of the `Outputs` module. The
  contents of the specified file(s) will be updated as the program completes each analysis
  frequency.

- `resultfile`:  A string containing the name of the results file. Defaults to `pssfss.res`.
  If this file already exists, it will be overwritten.  It is a binary
  file that contains information (including the generalized scattering matrix) from
  the analysis performed for each scan condition and frequency. The result file can be
  post-processed to produce similar or additional outputs that were requested at run time
  via the `outlist` argument.

- `showprogress`: If true, use ProgressMeter to show execution progress
"""
function _analyze(layers, sheets, junc, freqs, stkeys, stvalues;
    outlist=[], resultfile="pssfss.res", showprogress::Bool=true)
    showprogress && println("Beginning PSSFSS Analysis")
    ncount = 0 # Number of analyses performed
    ntotal = length(freqs) * length(stvalues[1]) * length(stvalues[2])
    showprogress && (progress = Progress(ntotal,1))
    showprogress && update!(progress, 0)
    isfile(resultfile) && rm(resultfile)
    date, clock = split(string(now()),'T')
    pssfssv = PkgVersion.Version(PSSFSS)
    ss = "Environment"
    io=IOBuffer(); versioninfo(io); juliainfo = String(take!(io)) * ss
    i = findfirst(ss, juliainfo)
    juliainfo = juliainfo[1:first(i)-1]
    juliainfo = juliainfo * "  Threads.nthreads() = $(Threads.nthreads())\n"
    @logfile "\n\nStarting PSSFSS $(pssfssv) analysis on $(date) at $(clock)\n$(juliainfo)\n\n"
    check_inputs(layers, sheets, junc, freqs, stkeys, stvalues, outlist)
    k0min, k0max = twopi*1e9/c₀ .* extrema(freqs)
    gbls = choose_gblocks(layers, sheets, junc, k0min)
    gsm_save = Vector{GSM}(undef, length(gbls)) # Storage for reusable GSMs
    choose_layer_modes!(layers, sheets, junc, gbls, k0max, dbmin)
    gbldup = get_gbldup(gbls, layers, sheets, junc)
    usi = unique_indices(sheets)
    rwgdat = RWGData[setup_rwg(sheet) for sheet in sheets]
    report_layers_sheets(layers, sheets, junc, rwgdat, usi)
    uvec::Vector{Float64} = map(sheets) do sh  # Green's function smoothing factors
        sh.style == "NULL" && (return 0.0)
        ufactor = 0.5 * ustrip(Float64, sh.units, 1u"m")
        ufactor * max(norm(sh.β₁), norm(sh.β₂))
    end
    # Begin analysis loops over steering angles and frequency
    firstoutput = true
    results = Result[]
    for stout in stvalues[1], stin in stvalues[2]
        steer = getsttuple(stkeys, stout, stin)
        if keys(steer)[1] == :ψ₁
            ψ₁, ψ₂ = deg2rad.([steer...]) # radians
            upm::Float64 = ustrip(Float64, sheets[1].units, 1u"m")
            β₁, β₂ = sheets[1].β₁*upm, sheets[1].β₂*upm
            β⃗₀₀ = (ψ₁*β₁ + ψ₂*β₂) / twopi # Eq. (2.13b)
        else
            θ, ϕ = steer # degrees
            st = sind(θ)
            sp, cp = sincosd(ϕ)
        end
        @logfile "Beginning $(steer)"
        for fghz in freqs
            @logfile "  $(fghz) GHz"
            t_freq = time()
            k0 = twopi*fghz*1e9/c₀
            if keys(steer)[1] == :θ
                k1 = k0 * sqrt(real(layers[1].ϵᵣ * layers[1].μᵣ))
                β⃗₀₀ = @SVector([k1*st*cp, k1*st*sp])
            end
            setup_modes!.(layers, k0, Ref(β⃗₀₀))
            if !(angle(layers[begin].γ[1]) ≈ angle(layers[end].γ[1]) ≈ π/2)
                @logfile "  Skipping $(fghz) GHz due to cutoff principal modes in ambient medium"
                continue
            end
            # Initialize overall GSM and propagate it through layer 1's width:
            n1 = length(layers[1].P)
            gsma = GSM(n1,n1)
            gsmc = GSM(n1,n1)
            cascade!(gsma, layers[1])
            for (ig,gbl) in pairs(gbls) # Walk through the Gblocks
                i1 = first(gbl.rng) # Index of layer to left of Gblock
                i2 = 1 + last(gbl.rng) # Index of layer to right of Gblock
                i_junc = gbl.j # junction where FSS is located, or 0 if no sheet
                i_sheet = i_junc == 0 ? 0 : junc[i_junc]
                if i_sheet ≠ 0
                    if gbldup[ig] > 0
                        gsmb = deepcopy(gsm_save[gbldup[ig]]) # Use previously calculated GSM
                    else
                        region = @view layers[i1:i2]
                        sheet = sheets[i_sheet]
                        s = gbl.j - i1 + 1 # sheet interface location within `region`
                        if sheet.class == 'J'
                            gsmb = calculate_jtype_gsm(region, sheet, uvec[i_sheet],
                                                       rwgdat[i_sheet], s, k0, β⃗₀₀, usi[i_sheet])
                        elseif sheet.class == 'M'
                            gsmb = calculate_mtype_gsm(region, sheet, uvec[i_sheet],
                                                       rwgdat[i_sheet], s, k0, β⃗₀₀, usi[i_sheet])
                        elseif sheet.class == 'E'
                            gsmb = pecgsm(length(region[1].P), length(region[end].P))
                        elseif sheet.class == 'H'
                            gsmb = pmcgsm(length(region[1].P), length(region[end].P))
                        else
                            error("Illegal sheet class: $(sheet.class)")
                        end

                        gbldup[ig] < 0 && (gsm_save[ig] = deepcopy(gsmb))
                    end
                    # Apply translations if requested:
                    sheet = sheets[i_sheet]
                    if sheet.dx ≠ 0 || sheet.dy ≠ 0
                        upm = ustrip(Float64, sheet.units, 1u"m")
                        dx = sheet.dx / upm
                        dy = sheet.dy / upm
                        translate_gsm!(gsmb, dx, dy, first(region), last(region))
                    end

                else # no sheet
                    @assert i2 - i1 == 1
                    gsmb = gsm_slab_interface(layers[i1], layers[i2], k0)
                end
                gsmc = cascade(gsma, gsmb)
                cascade!(gsmc, layers[i2])
                gsma = gsmc
            end # Gblock loop
            t_freq = round(time()-t_freq, digits=tdigits)
            @logfile "  $(t_freq) seconds total at $(fghz) GHz"

            result = Result(gsmc, steer, β⃗₀₀, fghz, layers[1].ϵᵣ, layers[1].μᵣ,
            layers[1].β₁, layers[1].β₂, layers[end].ϵᵣ, layers[end].μᵣ,
            layers[end].β₁, layers[end].β₂)
            push!(results, result)
            ncount += 1
            # Write to output files
            append_result_data(resultfile,string(ncount),result)
            for row in eachrow(outlist)
                if firstoutput
                    open(row[1], "w") do io
                        writedlm(io, permutedims([r.label for r in row[2]]), ',')
                    end
                end
                open(row[1], "a") do io
                    writedlm(io, permutedims([r(result) for r in row[2]]), ',')
                end
            end
            firstoutput = false
            showprogress && next!(progress) # Bump progress meter
        end # Frequency loop
    end # steering angle loop

    date, clock = split(string(now()),'T')
    @logfile "\n\n PSSFSS analysis exiting on $(date) at $(clock)\n\n"
    return results
end # function

"""
    calculate_jtype_gsm(layers, sheet::RWGSheet, u::Real, rwgdat::RWGData, s::Int, k0, k⃗inc, is::Int) -> gsm

Compute the generalized scattering matrix for a sheet of class `'J'`.

### Input Arguments

- `layers`: An iterable of `Layer` instances containing the layers for the `Gblock`
    associated with the `sheet` under consideration.
- `sheet`:  A sheet of class `'J'` for which the GSM is desired.
- `u`: The Green's function smoothing parameter for the `sheet`.
- `rwgdat`: The `RWGData` object associated with the `sheet` argument.
- `s`: The interface number within `layers` at which the sheet is located.
- `k0`: The free-space wavenumber in radians/meter.
- `k⃗inc`: A 2-vector containing the incident field wave vector x and y components. Note
    that by the phase match condition this vector is the same for all layers in the entire FSS
    structure.
- `is`: The sheet index for the `sheet` argument within the global list of sheets.

### Return Value

- `gsm::GSM`  The full GSM for the GBlock including incident fields and scattered fields due
    to currents induced on the sheet surface.
"""
function calculate_jtype_gsm(layers::AbstractVector{Layer}, sheet::RWGSheet, u::Float64,
                                 rwgdat::RWGData, s::Int, k0::Float64, k⃗inc::SVector{2, Float64}, is::Int)
    one_meter = ustrip(Float64, sheet.units, 1u"m")
    area = norm(sheet.s₁ × sheet.s₂) / one_meter^2 # Unit cell area (m^2).
    nmodesmax = max(length(layers[begin].P), length(layers[end].P))
    nbf = size(rwgdat.bfe,2) # Number of basis functions
    if isempty(rwgdat.bfftstore)
        rwgdat.bfftstore = zeros(SArray{Tuple{2},ComplexF64,1,2}, (nbf, 2, nmodesmax))
    end
    bfftstore = rwgdat.bfftstore
    # Compute area correction factors for the mode normalization constants of
    # the two end regions:
    tempvec = zeros(2)
    for (i,l) in enumerate(@view layers[[begin,end]])
      area_i = twopi * twopi / norm(l.β₁ × l.β₂)
      tempvec[i] = √(area_i / area)
    end
    acf = SVector(tempvec[1], tempvec[2])

    # Set up the partial GSM due to incident field
    (gsm, tlgfvi, vincs) = gsm_electric_gblock(layers, s, k0)
    if sheet.style == "NULL"
        return gsm
    end

    # Calculate the scattered field partial scattering matrix of the
    # FSS sheet at this junction. Then add it to GSM already computed
    # for the dielectric discontinuity...

    #  Fill the interaction matrix for the current sheet:
    t_temp = time()
    ψ₁ = k⃗inc ⋅ sheet.s₁ / one_meter
    ψ₂ = k⃗inc ⋅ sheet.s₂ / one_meter
    @logfile "    Beginning matrix fill for sheet $(is)"
    zmat = fillz(k0,u,layers,s,ψ₁,ψ₂,sheet,rwgdat)
    t_fill = round(time() - t_temp, digits=tdigits)
    @logfile "      $(t_fill) seconds total matrix fill time for sheet $(is)"
    # Factor the matrix:
    t_temp = time()
    zmatf = lu!(zmat)
    t_factor = round(time() - t_temp, digits=tdigits)
    @logfile "      $(t_factor) seconds to factor matrix for sheet $(is)"
    t_temp = time()
    # Compute and store the basis function Fourier transforms:
    i_ft = 0
    for (sr,l) in enumerate(@view layers[[begin,end]]) # loop over possible source regions
        for qp = 1:length(l.P)
            kvec = l.β[qp]
            # If desired F.T. has already been computed, then copy it.
            if qp > 1 && kvec ≈ l.β[qp-1]
                bfftstore[:,sr,qp] = bfftstore[:,sr,qp-1]
                continue
            end
            if sr == 2 && length(layers[1].β) ≥ qp && kvec ≈ layers[1].β[qp]
                bfftstore[:,2,qp] = bfftstore[:,1,qp]
                continue
            end
            bfft = @view bfftstore[:,sr,qp]
            rwgbfft!(bfft, rwgdat, sheet, kvec, ψ₁, ψ₂) # Otherwise, compute from scratch
            i_ft += 1
        end
    end
    t_fft = round(time() - t_temp, digits=tdigits)
    #@logfile "      $(t_fft) seconds for basis function Fourier transforms at $(i_ft) points"
    nsolve = 0
    t_extract = 0.0
    i_extract = 0
    t_solve = 0.0
    imat = rwgdat.rhs
    for (sr,ls) in enumerate(@view layers[[begin,end]]) # Loop over source regions
        for qp in 1:length(ls.P) # Loop over srce reg modes
            # Incident field for source layer in absence of the FSS sheet:
            sourcevec = vincs[qp,sr] * ls.c[qp] * acf[sr] * ls.tvec[qp]
            # Compute generalized voltage vector:
            imat .= (b ⋅ sourcevec for b in bfftstore[:,sr,qp]) # Eq. (7.39)
            # Solve the matrix equation
            t_solve1 = time()
            ldiv!(zmatf, imat)
            t_solve2 = time()
            t_solve += t_solve2 - t_solve1
            nsolve += 1
            t_extract1 = time()
            for (or,lo) in enumerate(@view layers[[begin,end]]) # Loop over obs. regions
                smat = gsm[or,sr]
                for q in 1:length(lo.P)  # Loop obs. regn modes
                    # Extract partial scattering parameter due to scattered fields...
                    FTJ = sum((imat[n] * bfftstore[n,or,q] for n in 1:nbf)) # FT of total current
                    smat[q,qp] -= (lo.tvec[q] ⋅ FTJ) * (tlgfvi[q,or] /
                                                (lo.c[q] * acf[or] * area)) # Eq (6.18)
                    i_extract += 1
                end
            end
            t_extract2 = time()
            t_extract = t_extract + (t_extract2 - t_extract1)
        end
    end
    #@logfile "      $(round(t_extract,digits=tdigits)) seconds to extract $(i_extract) GSM entries"
    return gsm
end

"""
    calculate_mtype_gsm(layers, sheet::RWGSheet, u::Real, rwgdat::RWGData, s::Int, k⃗inc, is::Int) -> gsm

Compute the generalized scattering matrix for a sheet of class `'M'`.

### Input Arguments

- `layers`: An iterable of `Layer` instances containing the layers for the `Gblock`
    associated with the `sheet` under consideration.
- `sheet`:  A sheet of class `'M'` for which the GSM is desired.
- `u`: The Green's function smoothing parameter for the `sheet`.
- `rwgdat`: The `RWGData` object associated with the `sheet` argument.
- `s`: The interface number within `layers` at which the sheet is located.
- `k0`: The free-space wavenumber in radians/meter.
- `k⃗inc`: A 2-vector containing the incident field wave vector x and y components. Note
    that by the phase match condition this vector is the same for all layers in the entire FSS
    structure.
- `is`: The sheet index for the `sheet` argument within the global list of sheets.

### Return Value

- `gsm::GSM`  The full GSM for the GBlock including incident fields and scattered fields due
    to magnetic currents induced in the gaps on the sheet surface.
"""
function calculate_mtype_gsm(layers::AbstractVector{Layer}, sheet::RWGSheet, u::Float64,
                                 rwgdat::RWGData, s::Int, k0::Float64, k⃗inc::SVector{2, Float64}, is::Int)
    one_meter = ustrip(Float64, sheet.units, 1u"m")
    area = norm(sheet.s₁ × sheet.s₂) / one_meter^2 # Unit cell area (m^2).
    nmodesmax = max(length(layers[begin].P), length(layers[end].P))
    nbf = size(rwgdat.bfe,2) # Number of basis functions
    if isempty(rwgdat.bfftstore)
        rwgdat.bfftstore = zeros(SArray{Tuple{2},ComplexF64,1,2}, (nbf, 2, nmodesmax))
    end
    bfftstore = rwgdat.bfftstore
    # Compute area correction factors for the mode normalization constants of
    # the two end regions:
    tempvec = zeros(2)
    for (i,l) in enumerate(@view layers[[begin,end]])
      area_i = twopi * twopi / norm(l.β₁ × l.β₂)
      tempvec[i] = √(area_i / area)
    end
    acf = SVector(tempvec[1], tempvec[2])

    # Set up the partial GSM due to incident field
    (gsm, tlgfiv, iincs) = gsm_magnetic_gblock(layers, s, k0)
    if sheet.style == "NULL"
        return gsm
    end

    # Calculate the scattered field partial scattering matrix of the
    # FSS sheet at this junction. Then add it to GSM already computed
    # for the dielectric discontinuity...

    #  Fill the interaction matrix for the current sheet:
    t_temp = time()
    ψ₁ = k⃗inc ⋅ sheet.s₁ / one_meter
    ψ₂ = k⃗inc ⋅ sheet.s₂ / one_meter
    @logfile "    Beginning matrix fill for sheet $(is)"
    ymat = filly(k0,u,layers,s,ψ₁,ψ₂,sheet,rwgdat)
    t_fill = round(time() - t_temp, digits=tdigits)
    @logfile "      $(t_fill) seconds total matrix fill time for sheet $(is)"
    # Factor the matrix:
    t_temp = time()
    ymatf = lu!(ymat)
    t_factor = round(time() - t_temp, digits=tdigits)
    @logfile "      $(t_factor) seconds to factor matrix for sheet $(is)"
    t_temp = time()
    # Compute and store the basis function Fourier transforms:
    i_ft = 0
    for (sr,l) in enumerate(@view layers[[begin,end]]) # loop over possible source regions
        for qp = 1:length(l.P)
            kvec = l.β[qp]
            # If desired F.T. has already been computed, then copy it.
            if qp > 1 && kvec ≈ l.β[qp-1]
                bfftstore[:,sr,qp] = bfftstore[:,sr,qp-1]
                continue
            end
            if sr == 2 && length(layers[1].β) ≥ qp && kvec ≈ layers[1].β[qp]
                bfftstore[:,2,qp] = bfftstore[:,1,qp]
                continue
            end
            bfft = @view bfftstore[:,sr,qp]
            rwgbfft!(bfft, rwgdat, sheet, kvec, ψ₁, ψ₂) # Otherwise, compute from scratch
            i_ft += 1
        end
    end
    t_fft = round(time() - t_temp, digits=tdigits)
    @logfile "      $(t_fft) seconds for basis function Fourier transforms at $(i_ft) points"
    nsolve = 0
    t_extract = 0.0
    t_solve = 0.0
    vmat = rwgdat.rhs
    i_extract = 0
    σ = -1
    for (sr,ls) in enumerate(@view layers[[begin,end]]) # Loop over source regions
        σ *= -1 # 1 for sr == 1, and -1 for sr == 2
        for qp in 1:length(ls.P) # Loop over srce reg modes
            # Incident field for Region sr (Eq. (7.64))
            sourcevec = iincs[qp,sr] * ls.c[qp] * ls.Y[qp] * zhatcross(ls.tvec[qp])
            # Compute generalized current vector:
            vmat .= (b ⋅ sourcevec for b in bfftstore[:,sr,qp]) # Eq. (7.64)
            # Solve the matrix equation
            t_solve1 = time()
            ldiv!(ymatf, vmat)
            t_solve2 = time()
            t_solve += t_solve2 - t_solve1
            nsolve += 1
            t_extract1 = time()
            for (or,lo) in enumerate(@view layers[[begin,end]]) # Loop over obs. regions
                smat = gsm[or,sr]
                for q in 1:length(lo.P)  # Loop obs. regn. modes
                    # Extract partial scattering parameter due to scattered fields...
                    FTM = sum((vmat[n] * bfftstore[n,or,q] for n in 1:nbf)) # FT of total mag. current
                    smat[q,qp] += (zhatcross(lo.tvec[q]) ⋅ FTM) *
                                        (σ * tlgfiv[q,or] * lo.c[q]) # Eq. (6.37)
                    i_extract += 1
                end
            end
            t_extract2 = time()
            t_extract += t_extract2 - t_extract1
        end
    end
    @logfile "      $(t_extract) seconds to extract $(i_extract) GSM entries"
    return gsm
end

"""
    getsttuple(stkeys::Tuple{Symbol,Symbol}, stout::Float64, stin::Float64) -> NamedTuple

Return a named tuple either of the form `(θ = θ, ϕ = ϕ)` or `(ψ₁ = ψ₁, ψ₂ = ψ₂)` that serves
to define the current steering situation. Actually, the input field names can be spelled out in
English as `:theta`, `:phi`, `:psi1`, and `:psi2`.

### Arguments

- `stkeys`: A 2-tuple containing steering parameters as `Symbols`, either (`:ψ₁` and `ψ₂`) or (`:θ` and `:ϕ`)
  (or their spelled-out English versions as detailed above), either of which
  could be listed in either order.  The order is significant in that the first member of the pair
  defines the outer steering loop.

- `stout` and `stin`: These are the current values of the outer and inner steering variables,
  respectively.
"""
function getsttuple(stkeys::Tuple{Symbol, Symbol}, stout::Float64, stin::Float64)
    if stkeys[1] ∈ (:phi, :ϕ, :Phi, :PHI, :Φ)
        return (θ = stin, ϕ = stout)
    elseif stkeys[2] ∈ (:phi, :ϕ, :Phi, :PHI, :Φ)
        return (θ = stout, ϕ = stin)
    elseif keys(steering)[1] ∈ (:psi1, :Psi1, :PSI1, :ψ₁, :Ψ₁, :ψ1, :Ψ1)
        return (ψ₁ = stout, ψ₂ = stin)
    else
        return (ψ₁ = stin, ψ₂ = stout)
    end
end

function check_inputs(layers, sheets, junc, freqs, stkeys, stvalues, outlist)
    # Check that input and output media are lossless and the same
    imag(first(layers).ϵᵣ) == imag(last(layers).ϵᵣ) ==
    imag(first(layers).μᵣ) == imag(last(layers).μᵣ) == 0 ||
    error("First and last layers must be lossless")
    (first(layers).ϵᵣ == last(layers).ϵᵣ) &&
    (first(layers).μᵣ == last(layers).μᵣ) ||
    error("First and last layers must have identical electrical parameters")

    # todo:

    # Check that ψ₁ and ψ₂ are not specified when there are no nonnull sheets

    return
end

"""
    unique_indices(v::Vector)

Return a vector `ui` of the same length as `v`.
`ui[k]` contains the smallest `i` such that 'i ≤ k` and `ui[i] === ui[k]`
"""
function unique_indices(v::Vector)
    n = length(v)
    ui = collect(1:n)
    for io in 1:n-1
        ui[io] ≠ io && continue
        for it in io+1:n
            ui[it] ≠ it && continue
            v[io] === v[it]  && (ui[it] = io)
        end
    end
    return ui
end

"""
    get_gbldup(gbls::Vector{Gblock}, layers::Vector{Layer}, sheets::Vector{RWGSheet}, junc::Vector{Int})
    -> (gbldup, junc)

Return `gbldup::Vector{Int}` of the same length as `gbls`.
`gbldup[k]` contains `0` for an ordinary Gblock.  `gbldup[k] == -1` means that
`gbls[k]` is the first occurence of a repeated Gblock and that its GSM should be
saved for reuse.  `gbldup[k] == i` where `0<i<k` means that `gbls[k]` is identical
to `gbls[i]` and they can both use the same GSM.
Two `Gblock`s are considered identical if they

1. Contain identical (`===`) `Sheet` objects at the same location within the block.
2. Comprise the same number of dielectric layers with identical widths analyzed
   electrical characteristics.
3. Are embedded within similar adjacent dielectric layers, having identical electrical
   properties and numbers of modes.

`junc::Vector{Int}` is has length `length(Layers)-1`. `junc[i]`` is the sheet number
present at dielectric interface `i`, or `0` if no sheet is present there.
"""
function get_gbldup(gbls::Vector{Gblock}, layers::Vector{Layer}, sheets::Vector{RWGSheet}, junc::Vector{Int})
    gbldup = zeros(Int, length(gbls))
    for (g1,gbl1) in pairs(gbls)
        (gbldup[g1] ≠ 0 || gbl1.j == 0) && continue
        j1, rng1 = gbl1.j, gbl1.rng
        n1 = length(layers[first(rng1)].P) # modes at left side
        n2 = length(layers[last(rng1)+1].P) # modes at right side
        # Examine succeeding Gblocks to see if they match:
        for g2 in g1+1:length(gbls)
            gbl2 = gbls[g2]
            rng2 = gbl2.rng
            j2 = gbl2.j
            (gbldup[g2] ≠ 0 || j2 == 0) && continue
            sheets[junc[j1]] === sheets[junc[j2]] || continue
            last(rng1)-j1 ≠ last(rng2)-j2 && continue
            # Check that layers within blocks gbl1 and gbl2 are identical:
            length(rng1) ≠ length(rng2) && continue
            n1 ≠ length(layers[first(rng2)].P) && continue
            n2 ≠ length(layers[1+last(rng2)].P) && continue
            all(zip(rng1,rng2)) do (i1,i2)
                i1 == first(rng1) && (return true)
                layers[i1] == layers[i2]
            end || continue
            # If we made it to here, the two Gblocks are identical:
            gbldup[g1] = -1  # Indicate that GSM of Gblock g1 is to be saved
            gbldup[g2] = g1  # GSM of Gblock g2 is obtained from saved GSM of block g1
        end
    end
    return gbldup
end # function

function report_layers_sheets(layers, sheets, junc, rwgdat, usi)
    @logfile "Dielectric layer information... \n"
    @logfile " Layer  Width  units  epsr   tandel   mur  mtandel modes  beta1x  beta1y  beta2x  beta2y"
    @logfile " ----- ------------- ------- ------ ------- ------ ----- ------- ------- ------- -------"
    for (jl,l) in pairs(layers)
        eps = real(l.ϵᵣ)
        tandel = -imag(l.ϵᵣ) / eps
        mu = real(l.μᵣ)
        mtandel = -imag(l.μᵣ) / mu
        nmode = length(l.P)
        units = string(unit(l.user_width))
        units == "inch" && (units = "in")
        uw_unitless = ustrip(l.user_width)
        str = @sprintf(" %5i %9.4f %3s %7.2f %6.4f %7.2f %6.4f %5i %7.1f %7.1f %7.1f %7.1f",
              jl,uw_unitless, units, eps, tandel, mu, mtandel, nmode, l.β₁[1], l.β₁[2],
              l.β₂[1], l.β₂[2])
        @logfile "$str"
        if jl < length(layers) && junc[jl] ≠ 0
            js = junc[jl]
            s = sheets[js]
            om = ustrip(Float64, s.units, 1.0u"m")
            str = @sprintf(
                " ==================  Sheet %3i  ======================== %7.1f %7.1f %7.1f %7.1f",
                usi[js], s.β₁[1] * om, s.β₁[2] * om, s.β₂[1] * om, s.β₂[2] * om)
            @logfile "$str"
        end
    end
    sint = findall(junc .≠ 0) # sint[k] contains dielectric interface number of k'th sheet

    @logfile "\n\n\nPSS/FSS sheet information...\n"
    @logfile "Sheet  Loc         Style      Rot  J/M Faces Edges Nodes Unknowns  NUFP"
    @logfile "-----  ---  ---------------- ----- --- ----- ----- ----- -------- ------"
    for (js,s) in pairs(sheets)
        str = @sprintf("%4i   %3i  %16s %5.1f  %1s  %5i %5i %5i  %6i %7i",
             usi[js], sint[js], s.style, s.rot, s.class, size(s.fe,2), length(s.e1),
             length(s.ρ), size(rwgdat[js].bfe,2), length(rwgdat[js].ufp2fp))
        @logfile "$str"
    end
    @logfile "\n\n"
    nothing
end

end # module


text:
path = dirname(@__FILE__)

# Set up arguments
m = AnSchorfheide(testing = true)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:forecast_horizons, 1)

system = JLD2.jldopen("$path/../reference/forecast_args.jld2","r") do file
    read(file, "system")
end
z0 = zeros(n_states_augmented(m))

# Read expected output
exp_states, exp_obs, exp_pseudo, exp_shocks =
    JLD2.jldopen("$path/../reference/forecast_out.jld2", "r") do file
        read(file, "exp_states"),
        read(file, "exp_obs"),
        read(file, "exp_pseudo"),
        read(file, "exp_shocks")
    end

# Without shocks
global states, obs, pseudo, shocks = forecast(m, system, z0; draw_shocks = false)

@testset "Testing forecasting without drawing shocks" begin
    @test @test_matrix_approx_eq exp_states states
    @test @test_matrix_approx_eq exp_obs    obs
    @test @test_matrix_approx_eq exp_pseudo pseudo
    @test @test_matrix_approx_eq exp_shocks shocks
end

# Supplying shocks
global states, obs, pseudo, shocks = forecast(m, system, z0; shocks = shocks)

@testset "Testing forecasting with pre-supplied shocks" begin
    @test @test_matrix_approx_eq exp_states states
    @test @test_matrix_approx_eq exp_obs    obs
    @test @test_matrix_approx_eq exp_pseudo pseudo
    @test @test_matrix_approx_eq exp_shocks shocks
end

# Draw normally distributed shocks
global states, obs, pseudo, shocks = forecast(m, system, z0; draw_shocks = true)

# Draw t-distributed shocks
m <= Setting(:forecast_tdist_shocks, true)
global states, obs, pseudo, shocks = forecast(m, system, z0; draw_shocks = true)
m <= Setting(:forecast_tdist_shocks, false)

# Enforce ZLB
ind_r = m.observables[:obs_nominalrate]
ind_r_sh = m.exogenous_shocks[:rm_sh]
zlb_value = forecast_zlb_value(m)
shocks = zeros(n_shocks_exogenous(m), forecast_horizons(m))
shocks[ind_r_sh, :] .= -10.

@testset "Ensure valid forecasting at the ZLB" begin
    global states, obs, pseudo, shocks = forecast(m, system, z0; shocks = shocks)
    @test all(x -> x < zlb_value, obs[ind_r, :])

    global states, obs, pseudo, shocks = forecast(m, system, z0; shocks = shocks, enforce_zlb = true)
    @test all(x -> abs(x - zlb_value) < 0.01, obs[ind_r, :])
    @test all(x -> x != -10.,                 shocks[ind_r_sh, :])
end

@testset "Test that forecasting and populating shocks under alternative policy works" begin
    m <= Setting(:alternative_policy, AltPolicy(:taylor93, eqcond, solve))
    @test typeof(forecast(m, system, z0; draw_shocks = true)) == NTuple{4, Array{Float64, 2}}
end

nothing