text:
@testset "rational.jl" begin
    rng = MersenneTwister(123456)
    x = rand(rng) * 2
    v1 = rand(rng, 3)
    v2 = rand(rng, 3)

    @testset "RationalKernel" begin
        α = rand()
        k = RationalKernel(; α=α)

        @testset "RationalKernel ≈ Exponential for large α" begin
            @test isapprox(
                RationalKernel(; α=1e9)(v1, v2),
                ExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @test metric(RationalKernel()) == Euclidean()
        @test metric(RationalKernel(; α=α)) == Euclidean()
        @test repr(k) == "Rational Kernel (α = $(α), metric = Euclidean(0.0))"

        k2 = RationalKernel(; α=α, metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(x -> RationalKernel(; alpha=exp(x[1])), [α])
        test_params(k, ([α],))
    end

    @testset "RationalQuadraticKernel" begin
        α = rand()
        k = RationalQuadraticKernel(; α=α)

        @testset "RQ ≈ EQ for large α" begin
            @test isapprox(
                RationalQuadraticKernel(; α=1e9)(v1, v2),
                SqExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @test metric(RationalQuadraticKernel()) == SqEuclidean()
        @test metric(RationalQuadraticKernel(; α=α)) == SqEuclidean()
        @test repr(k) == "Rational Quadratic Kernel (α = $(α), metric = Euclidean(0.0))"

        k2 = RationalQuadraticKernel(; α=α, metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        test_ADs(x -> RationalQuadraticKernel(; alpha=exp(x[1])), [α])
        test_params(k, ([α],))
    end

    @testset "GammaRationalKernel" begin
        k = GammaRationalKernel()

        @test repr(k) == "Gamma Rational Kernel (α = 2.0, γ = 1.0, metric = Euclidean(0.0))"

        @testset "GammaRational (γ=2) ≈ RQ with rescaled inputs" begin
            @test isapprox(
                GammaRationalKernel(; γ=2)(v1 ./ sqrt(2), v2 ./ sqrt(2)),
                RationalQuadraticKernel()(v1, v2),
            )
            a = 1 + rand()
            @test isapprox(
                GammaRationalKernel(; α=a, γ=2)(v1 ./ sqrt(2), v2 ./ sqrt(2)),
                RationalQuadraticKernel(; α=a)(v1, v2),
            )
        end

        @testset "GammaRational (γ=2) ≈ EQ for large α with rescaled inputs" begin
            v1 = randn(2)
            v2 = randn(2)
            @test isapprox(
                GammaRationalKernel(; α=1e9, γ=2)(v1 ./ sqrt(2), v2 ./ sqrt(2)),
                SqExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @testset "Default GammaRational ≈ Rational" begin
            @test isapprox(GammaRationalKernel()(v1, v2), RationalKernel()(v1, v2))
            a = 1 + rand()
            @test isapprox(
                GammaRationalKernel(; α=a)(v1, v2), RationalKernel(; α=a)(v1, v2)
            )
        end

        @testset "Default GammaRational ≈ Exponential for large α" begin
            v1 = randn(4)
            v2 = randn(4)
            @test isapprox(
                GammaRationalKernel(; α=1e9)(v1, v2),
                ExponentialKernel()(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @testset "GammaRational ≈ GammaExponential for same γ and large α" begin
            v1 = randn(3)
            v2 = randn(3)
            γ = rand() + 0.5
            @test isapprox(
                GammaRationalKernel(; α=1e9, γ=γ)(v1, v2),
                GammaExponentialKernel(; γ=γ)(v1, v2);
                atol=1e-6,
                rtol=1e-6,
            )
        end

        @test metric(GammaRationalKernel()) == Euclidean()
        @test metric(GammaRationalKernel(; γ=2.0)) == Euclidean()
        @test metric(GammaRationalKernel(; γ=2.0, α=3.0)) == Euclidean()

        k2 = GammaRationalKernel(; metric=WeightedEuclidean(ones(3)))
        @test metric(k2) isa WeightedEuclidean
        @test k2(v1, v2) ≈ k(v1, v2)

        # Standardised tests.
        TestUtils.test_interface(k, Float64)
        a = 1.0 + rand()
        test_ADs(x -> GammaRationalKernel(; α=x[1], γ=x[2]), [a, 1 + 0.5 * rand()])
        test_params(GammaRationalKernel(; α=a, γ=x), ([a], [x]))
    end
end


text:
"""
Unpenalized

`g(θ) = 0`
"""
struct NoPenalty <: ProxableElementPenalty end
value(p::NoPenalty, θ::Number) = zero(θ)
deriv(p::NoPenalty, θ::Number) = zero(θ)
prox(p::NoPenalty,  θ::Number, s::Number) = θ

"""
L1Penalty aka LASSO

`g(θ) = abs(θ)`
"""
struct L1Penalty <: ProxableElementPenalty end
value(p::L1Penalty, θ::Number) = abs(θ)
deriv(p::L1Penalty, θ::Number) = sign(θ)
prox(p::L1Penalty,  θ::Number, s::Number) = soft_thresh(θ, s)

"""
L2Penalty aka Ridge

`g(θ) = .5 * θ ^ 2`
"""
struct L2Penalty <: ProxableElementPenalty end
value(p::L2Penalty, θ::T) where {T <: Number} = inv(T(2)) * θ * θ
deriv(p::L2Penalty, θ::Number) = θ
prox(p::L2Penalty, θ::Number, s::Number) = θ / (one(θ) + s)

"""
ElasticNetPenalty, weighted average of L1Penalty and L2Penalty

`g(θ) = α * abs(θ) + (1 - α) * .5 * θ ^ 2`
"""
struct ElasticNetPenalty{T <: Number} <: ProxableElementPenalty
    α::T
    function ElasticNetPenalty(α::T = 0.5) where {T <: Number}
        0 <= α <= 1 || throw(ArgumentError("α must be in [0, 1]"))
        new{T}(α)
    end
end
for f in (:value, :deriv)
    @eval function ($f)(p::ElasticNetPenalty{T}, θ::Number) where {T <: Number}
        p.α * ($f)(L1Penalty(), θ) + (one(T) - p.α) * ($f)(L2Penalty(), θ)
    end
end
function prox(p::ElasticNetPenalty{T}, θ::Number, s::Number) where {T <: Number}
    αs = p.α * s
    soft_thresh(θ, αs) / (one(T) + s - αs)
end

"""
LogPenalty(η)

`g(θ) = log(1 + η * θ)`
"""
struct LogPenalty{T <: Number} <: ElementPenalty
    η::T
    function LogPenalty(η::T = 1.0) where {T <: Number}
        η > 0 || throw(ArgumentError("η must be > 0"))
        new{T}(η)
    end
end
value(p::LogPenalty, θ::Number) = log1p(p.η * abs(θ))
deriv(p::LogPenalty{T}, θ::Number) where {T <: Number} = p.η * sign(θ) / (one(T) + p.η * abs(θ))

# http://www.pstat.ucsb.edu/student%20seminar%20doc/SCAD%20Jian%20Shi.pdf
"""
Smoothly Clipped Absolute Deviation Penalty
"""
struct SCADPenalty{T <: Number} <: ElementPenalty
    a::T
    γ::T
    function SCADPenalty{T}(a::T, γ::T) where T<:Number
        a > 2 || throw(ArgumentError("First parameter must be > 2"))
        γ > 0 || throw(ArgumentError("Second parameter must be > 0"))
        new{T}(a, γ)
    end
end
SCADPenalty(a::T = 3.7, γ::T = T(1)) where {T <: Number} = SCADPenalty{T}(a, γ)
SCADPenalty(a::Number, γ::Number) = SCADPenalty(promote(a, γ)...)

function value(p::SCADPenalty{T}, θ::S) where {T, S <: Number}
    absθ = abs(θ)
    R = float(promote_type(T, S))
    if absθ < p.γ
        R(p.γ * absθ)
    elseif absθ <= p.γ * p.a
        -R(0.5) * (absθ^2 - R(2) * p.a * p.γ * absθ + p.γ^2) / (p.a - one(R))
    else
        R(0.5) * (p.a + one(R)) * p.γ * p.γ
    end
end
function deriv(p::SCADPenalty{T}, θ::S) where {T, S <: Number}
    absθ = abs(θ)
    R = float(promote_type(T, S))
    if absθ < p.γ
        R(p.γ * sign(θ))
    elseif absθ <= p.γ * p.a
        R(-(absθ - p.a * p.γ * sign(θ)) / (p.a - one(R)))
    else
        zero(R)
    end
end
# function prox{T}(p::SCADPenalty{T}, θ::T, λ::T)
#     absθ = abs(θ)
#     if absθ < 2λ
#         return prox(L1Penalty(), θ, λ)
#     elseif absθ <= λ * p.a
#         return prox(L1Penalty(), θ, λ * p.a / (p.a - 1)) * (p.a - 1) / (p.a - 2)
#     else
#         return θ
#     end
# end

# https://arxiv.org/abs/1002.4734
"""
MCPPenalty(γ) (MC+)
"""
struct MCPPenalty{T <: Number} <: ElementPenalty
    γ::T  # In paper, this is λ * γ
    function MCPPenalty(γ::T = 2.0) where T
        γ > 0 || throw(ArgumentError("γ must be > 0"))
        new{T}(γ)
    end
end
MCPPenalty(γ::Integer) = MCPPenalty(Float64(γ))
function value(p::MCPPenalty{T}, θ::Number) where {T <: Number}
    t = abs(θ)
    t < p.γ ? t - t^2 / (T(2) * p.γ) : (T(1)/T(2)) * p.γ
end
function deriv(p::MCPPenalty{T}, θ::S) where {T <: Number, S <: Number}
    t = abs(θ)
    t < p.γ ? sign(θ) * (T(1) - t / p.γ) : zero(float(promote_type(S,T)))
end

function _scale_check(λ)
    isa(λ, Number) || throw(ArgumentError("Scale factor λ must be a Number"))
    λ >= 0 || throw(ArgumentError("Scale factor λ has to be strictly positive."))
end

struct ScaledElementPenalty{T <: Number, P <: ElementPenalty} <: ElementPenalty
    penalty::P
    λ::T
end
scaled(p::ElementPenalty, λ::Number) = (_scale_check(λ); ScaledElementPenalty(p, λ))
Base.show(io::IO, sp::ScaledElementPenalty) = print(io, "$(sp.λ) * ($(sp.penalty))")

Base.:(*)(λ::Number, p::ElementPenalty) = scaled(p, λ)

value(p::ScaledElementPenalty{<:Number}, θ::Number) = p.λ * value(p.penalty, θ)
deriv(p::ScaledElementPenalty{<:Number}, θ::Number) = p.λ * deriv(p.penalty, θ)
prox(p::ScaledElementPenalty{<:Number},  θ::Number) = prox(p.penalty, θ, p.λ)
prox(p::ScaledElementPenalty{<:Number},  θ::AbstractArray{<:Number}) = prox(p.penalty, θ, p.λ)

# SCAD is special
for f in (:value, :deriv)
    @eval function ($f)(p::ScaledElementPenalty{<:Number, <:SCADPenalty}, θ::Number)
        ($f)(p.penalty, θ, p.λ)
    end
end

ext: jl


text:
struct TelegramLogger <: AbstractLogger
    tg::TelegramClient
    fmt::Function
    min_level::LogLevel
    async::Bool
end

"""
    TelegramLogger(tg::TelegramClient; fmt = tg_formatter,
                   min_level = Logging.Info, async = true)

Creates logger, which output log messages to Telegram channel.

# Arguments
- `tg`: telegram client, should have valid `chat_id` which will be used for outcome messages.
- `fmt`: function which accepts (`level`, `_module`, `group`, `id`, `file`, `line`) arguments and outputs message prefix. More details can be found in [Logging](https://docs.julialang.org/en/v1/stdlib/Logging/#Logging.handle_message) module. By default each messages is prepended with uppercase level, e.g. "INFO: " and the like.
- `min_level`: minimum level of log messages which is going to be processed.
- `async`: send messages in `sync` or `async` mode. Since telegram messages can take some time to be send, it make sense to set this parameter to `true`, so messages is send in the background with little impact on main program. But since in one run scripts main julia process can be finished before async message is sent, it make sense to set this parameter to `false`
"""
function TelegramLogger(tg::TelegramClient; fmt = tg_formatter, min_level = Info, async = true)
    return TelegramLogger(tg, fmt, min_level, async)
end

function tg_formatter(level, _module, group, id, file, line)
   return uppercase(string(level)) * ": "
end

function handle_message(tglogger::TelegramLogger, level, message, _module, group, id,
                        filepath, line; kwargs...)
    iob = IOBuffer()

    prefix = tg_formatter(level, _module, group, id, filepath, line)
    print(iob, prefix, string(message))
    for (key, val) in kwargs
        print(iob, key, " = ", val)
    end
    text = String(take!(iob))
    if tglogger.async
        @async sendMessage(tglogger.tg; text = text)
    else
        sendMessage(tglogger.tg; text = text)
    end
end

shouldlog(tglogger::TelegramLogger, arg...) = true
min_enabled_level(tglogger::TelegramLogger) = tglogger.min_level


text:
# f(t,u,du,res) = 0
"""
$(TYPEDEF)
"""
struct DAEProblem{uType,duType,tType,isinplace,P,F,K,D} <: AbstractDAEProblem{uType,duType,tType,isinplace}
  f::F
  du0::duType
  u0::uType
  tspan::tType
  p::P
  kwargs::K
  differential_vars::D
  @add_kwonly function DAEProblem{iip}(f::AbstractDAEFunction{iip},
                      du0,u0,tspan,p=NullParameters();
                      differential_vars = nothing,
                      kwargs...) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(du0),typeof(_tspan),
               isinplace(f),typeof(p),
               typeof(f),typeof(kwargs),
               typeof(differential_vars)}(
               f,du0,u0,_tspan,p,
               kwargs,differential_vars)
  end

  function DAEProblem{iip}(f,du0,u0,tspan,p=NullParameters();kwargs...) where {iip}
    DAEProblem(convert(DAEFunction{iip},f),du0,u0,tspan,p;kwargs...)
  end
end

function DAEProblem(f::AbstractDAEFunction,du0,u0,tspan,p=NullParameters();kwargs...)
  DAEProblem{isinplace(f)}(f,du0,u0,tspan,p;kwargs...)
end

function DAEProblem(f,du0,u0,tspan,p=NullParameters();kwargs...)
  DAEProblem(convert(DAEFunction,f),du0,u0,tspan,p;kwargs...)
end


text:
using heuristicopt

function main()
    # 50 individual , 100 gene
    population = rand(0:1, 50, 100)
    crossoverRate = 0.5
    mutationRate = 0.1
    iterations = 10000
    gaparam = GA.GAParameter(
                         crossoverRate,
                         mutationRate,
                         iterations
                         )

    select = GATool.selectroulette
    crossover = GATool.crossoverone
    mutate = GATool.mutaterandom
    gafunc = GA.GAFunction(select, crossover, mutate)
    model = GA.GAModel(
        gaparam,
        gafunc,
        objfunc
    )
    predictPopulation, bestFitness = GA.predict(model, population)
    print("bestFitness", bestFitness)
    print(predictPopulation)
end

function objfunc(individual)
    sum(individual)
end

main()

max_stars_repo_name: ariririri/heuristicopt

max_issues_repo_name: ariririri/heuristicopt


text:
using JLD2
using Statistics

num_bins = 12

mutable struct TTMazeStateVisitation <: LoggerKeyData
    state_visitations::Array{Int32,3}
    temp_state_visitations::Array{Int32, 2}
    log_interval::Int

    function TTMazeStateVisitation(logger_init_info)
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1
        new(zeros(7, 9, num_logged_steps), zeros(7, 9), logger_init_info[LoggerInitKey.INTERVAL])
    end
end

function lg_start!(self::TTMazeStateVisitation, env, agent)
end

function lg_step!(self::TTMazeStateVisitation, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)

    mask = valid_state_mask()
    env_state = convert(Integer, s[1])
    self.temp_state_visitations[mask[env_state][2], mask[env_state][1]] += 1

    if (rem(cur_step_total, self.log_interval) == 0)
        ind = fld(cur_step_total, self.log_interval) + 1
        self.state_visitations[:, :, ind] = self.temp_state_visitations

        # Reset temp state Visitation
        # println(sum(self.state_visitations))
        self.temp_state_visitations .*= 0
    end
end

function lg_episode_end!(self::TTMazeStateVisitation, cur_step_in_episode, cur_step_total)
end

function save_log(self::TTMazeStateVisitation, save_dict::Dict)
    save_dict[:ttmaze_state_visitation] = self.state_visitations
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
####################### kernel for SARA's outer loop ###########################
# At SLAC in Nov. 2019, used l_T = .03, l_τ = .1 in normalized domain [0, 1]²
# without normalizing, that corresponds to
# l_T = 1400 * 0.03 = 42, and
# l_τ = 2 * 0.1 = .2 (the first number is the range of the true domain of T and log10_τ)
function oSARA(l::AbstractVector)
    k = Kernel.MaternP(2)
    E = Diagonal(inv.(l).^2)
    Kernel.Energetic(k, E)
end

# length scales of temperature, dwelltime, and composition
oSARA(l_T::Real, l_τ::Real, l_c::Real) = oSARA([l_T, l_τ, l_c])
oSARA(l_T::Real, l_τ::Real) = oSARA([l_T, l_τ])

############################### outer loop #####################################
# returns next stripe conditions and conditional GP
function next_stripe(G::Gaussian, x::AbstractArray{<:NTuple}, y::AbstractArray{<:Real},
                     σ²::AbstractArray{<:Real},
                     stripe_conditions::AbstractVector{<:NTuple}, policy)
    return next_stripe(vec(x), vec(y), vec(σ²), stripe_conditions, policy)
end

function next_stripe(G::Gaussian, x::AbstractVector{<:NTuple}, y::AbstractVector{<:Real},
                     σ²::AbstractVector{<:Real},
                     stripe_conditions::AbstractVector{<:NTuple}, policy)
    length(x) == length(y) == length(σ²) || throw(DimensionMismatch("input data has different lengths: $(length(x)), $(length(y)), $(length(σ²))"))
    k_scale = 1/4 # this parameter is the scaling parameter of the kernel, but isn't incorporated in k due to performance bug
    σ² /= k_scale^2
    C = conditional(G, x, Gaussian(y, σ²), tol = 1e-6) # condition GP on data
    return policy(C, stripe_conditions), C
end


text:
# # How to Make a stack plot

# See [How to set up plots](3.0_set_up_plots.md) to get started

# ```julia
# using PowerGraphics
# const PG = PowerGraphics
# ```

# ### Make stack plots of results

# ```julia
# stack_plot(results)
# ```

# ![this one](plots-10/P__ThermalStandard_Stack.png)
# ![this one](plots-10/Example_Stack_Plot.png)

# ### Save the stack plots to a folder

# ```julia
# folder_path = joinpath(file_path, "plots_1")
# if !isdir(folder_path)
#    mkdir(folder_path)
# end
# stack_plot(results; save = folder_path)
# ```

# ![this one](plots-11/P__ThermalStandard_Stack.png)
# ![this one](plots-11/Example_saved_Stack_Plot.png)

# ### Show reserves in the stack plots

# ```julia
# stack_plot(results; reserves = true)
# ```
# ![this one](plots-12/P__ThermalStandard_Stack.png)
# ![this one](plots-12/Example_Stack_Plot_with_Reserves.png)

# ### Set different colors for the plots

# ```julia
# colors = [:pink :green :blue :magenta :black]
# stack_plot(results; seriescolor = colors)
# ```
# ![this one](plots-13/P__ThermalStandard_Stack.png)
# ![this one](plots-13/Example_Stack_Plot_with_Other_Colors.png)

# ### Create a Stair Plot, instead of interpolating between values

# ```julia
# stair_plot(results)
# ```
# or
# ```julia
# stack_plot(results; stair = true)
# ```
# ![this one](plots-14/P__ThermalStandard_Stair.png)
# ![this one](plots-14/Stair_Plot.png)

# ### Set a title

# ```julia
# title = "Example of a Title"
# stack_plot(results; title = title)
# ```
# ![this one](plots-15/P__ThermalStandard_Stack.png)
# ![this one](plots-15/Example_of_a_Title.png)

# ### For saving the plot with the PlotlyJS backend, you can set a different format for saving
# ```julia
# stack_plot(results; save = path, format = "html")
# ```
# Default format for saving is png.
# Optional formats for saving include png, html, and svg.


text:
using Colors, Interpolations
import Base: getindex, length, setindex!, push!, lastindex, eachindex, size
import Sockets: IPAddr
import FixedPointNumbers

"""
    AbstractChannel

    This channel abstraction allows me to use circular referential types between
    LEDStrip, LEDChannel, and LEDController because Julia doesn't yet support the
    forward lookahead necessary to use this feature with only the child types.
"""
abstract type AbstractChannel end
"""
    AbstractController

    This does the same for Controllers that AbstractChannel did for Channels:
    it allows me to use typed circular references without forward lookahead.
"""
abstract type AbstractController end

"""
    LEDStrip{T<:AbstractChannel, N<:AbstractController}

    Summary:
        Represents a string of LEDs in an Array-like form

    Details:
        This encapsulates the core representation of a string of LED's.
        It contains a view (a reference to a section of an array) of a portion of
        its controller's memory addresses.  It also has a reference to its controller
        which happends to be circular.

    Fields:
        name::String - A string representing the name of the strip of LEDs
            like "living room wall"
        controller::AbstractController - A reference to the LEDController
            which this particular LEDStrip is attached to
        subArray::SubArray - A reference to the relevant portion of the
            controller's LED memory
        idxRange::UnitRange - A range reflecting which addresses this
            LEDStrip refers to on its attached controller

    Inner Constructors:
        LEDStrip{T,N}(name::String, channel::T, controller::N, startAddr::Int, endAddr::Int) where {T<:AbstractChannel, N<:AbstractController}
            returns: LEDStrip{T,N} initialized with provided values
"""
mutable struct LEDStrip{T<:AbstractChannel, N<:AbstractController}
    name::String
    controller::N
    subArray::SubArray
    idxRange::UnitRange
    function LEDStrip{T,N}(name::String, channel::T, controller::N, startAddr::Int, endAddr::Int) where {T<:AbstractChannel, N<:AbstractController}
        x = new(name, controller, view(controller.addrs, startAddr:endAddr, :), startAddr:endAddr)
        push!(channel, x)
        push!(controller, x)
        return x
    end
end

#=
    Array-like interface implementation (Julia doesn't have interfaces,
    but there are collections of functions informally expected to exist)

    Most of these use some form of redirection to the inner view of the
    LEDController memory, but some have significant customizations on that.
=#
"Implements an override for the base length function"
length(strip::LEDStrip{T, N}) where {T<:AbstractChannel, N<:AbstractController}= length(strip.subArray)

"Implements getting an index from the strip"
getindex(strip::LEDStrip{T, N}, i::Any...) where {T<:AbstractChannel, N<:AbstractController}= getindex(strip.subArray, i...)

"Implements the equivalent of Python's enumerate() instance method, returns a range of indices to iterate over"
eachindex(s::LEDStrip{T, N}) where {T<:AbstractChannel, N<:AbstractController}=eachindex(s.subArray)

"Implements a lastindex base function that returns the index of the last element"
lastindex(s::LEDStrip{T, N}) where {T<:AbstractChannel, N<:AbstractController} = lastindex(s.subArray)

"Implements a size function for the LEDStrip, overrides the base function"
size(s::LEDStrip{T, N}, vargs...) where {T<:AbstractChannel, N<:AbstractController} = size(s.subArray, vargs...)

"Implements a function to set a particular index to a value for LEDStrip"
setindex!(strip::LEDStrip{T, N}, val::Any, idx::Any...) where {T<:AbstractChannel, N<:AbstractController} = setindex!(strip.subArray, val, idx...)

"Implements a vectorized form of setindex! to ensure maximum efficiency"
function setindex!(strip::LEDStrip{T, N}, val::Union{Array{UInt8},UInt8}, idx::Union{AbstractArray, Number, Colon}) where {T<:AbstractChannel, N<:AbstractController}
    return setindex!(strip.subArray, val, idx)
end

"""
    Implements another vectorized form to cover slightly different input
    params, but still want to type specialize for maximum performance
"""
function setindex!(
                   strip::LEDStrip{T, N},
                   val::Union{Array{UInt8},UInt8},
                   idx::Union{AbstractArray, Number, Colon},
                   idx2::Union{AbstractArray, Number, Colon}
                  ) where {T<:AbstractChannel, N<:AbstractController}
    return setindex!(strip.subArray, val, idx, idx2)
end

"""
    LEDController <: AbstractController

    Summary:
        Represents a particular networked LED controller.  Contains most of the memory for "real" LEDs

    Details:
        Contains the core representation of a LEDController for Irradiance.
        This class contains a large portion of the actual memory allocations for the system
        because it holds the memory for physical LED addresses.

    Fields:
        addrs::Array{UInt8, 2} - The address space for this controller, one dimension is for different LEDs
            and the other channel is for the different color channels
        strips::Array{LEDStrip} - An array containing all LEDStrips that are controlled by this controller
            and thus have their memory mapped into this controller
        location::Tuple{IPAddr,Int} - A tuple of an IPAddr, which is the IP address of the controller, and an
            Int, which is the port to contact the controller on.
        raw_data::Array{UInt8, 1} - A flattened form of the addrs field, for sending over the wire to the
            controller.  In this array each LED gets three adjacent indexes and then the next LED gets the
            next three indexes, and so on.
    Inner Constructors:
        LEDController(sz::Int, location::Tuple{IPAddr,Int})
            returns: LEDController initialized with no LEDStrips attached and zeroed data arrays

"""
mutable struct LEDController <: AbstractController
    addrs::Array{UInt8, 2}
    strips::Array{LEDStrip}
    location::Tuple{IPAddr, Int}
    raw_data::Array{UInt8, 1}
    function LEDController(sz::Int, location::Tuple{IPAddr, Int})
        return new(zeros(UInt8, sz, 3), Array{LEDStrip, 1}(0), location, Array{UInt8, 1}(sz*3))
    end
end

#=
    The below section implements the Array-like interface for LEDControllers
    Once again, Julia doesn't have proper interfaces, but there are informal ones
    made of collections of methods
=#
"Implements the length function for LEDControllers"
length(controller::LEDController) = length(controller.addrs)

"Implements a push function for LEDControllers"
push!(controller::LEDController, val::LEDStrip{T, N}) where {T<:AbstractChannel, N<:AbstractController}= push!(controller.strips, val)

"Implments a getindex function for LEDControllers"
getindex(controller::LEDController, idx...) = getindex(controller.addrs, idx...)

"Implements the eachindex function for LEDControllers, returns a range of all valid indices for the LEDController"
eachindex(c::LEDController) = 1:length(c)

"Implements the lastindex function for LEDControllers, returns the last valid index of the LEDController"
lastindex(c::LEDController) = length(c)

"Implements the setindex! function for LEDControllers, sets an index of the LEDController to a given value"
function setindex!(controller::LEDController, val::Any, idx::Any...)
    setindex!(controller.addrs, val, idx...)
end

"""
    LEDChannel <: AbstractChannel

    Summary:
        Abstracts LEDStrips into amorpheous channels to which effects are applied.

    Details:
        This contains the implementation of a AbstractChannel for Irradiance. This is the main abstraction used
        in the rest of the program for a collection of LED's, all effects are applied to LEDChannels, not
        LEDStrips, so that effects can be independent of the physical layout of the LED's

    Fields:
        map::Array{Tuple{LEDStrip{LEDChannel,LEDController}, Float64, Float64, AFArray{Float32,1}}, 1} -
            An array of large tuples that contain four things: the LEDStrip, its start on the channel, its
            end on the channel, and the array on the GPU that contains the pseudo-indexes for that LEDStrip
        virtualmem::Array{UInt8} - The abstraction of LEDs provided to effects so that they don't have to
            worry about the details of LED density, position, etc.  This way effects are much more portable
            across setups given that those setups have a reasonable LEDChannel map setup.
        gpu_virtualmem::AFArray{Float32,2} - The memory allocation on the GPU used for the virtual memory
            mentioned above.  Since all of the interpolation is executed on the GPU, we allocate arrays on
            the GPU semi-permanently so that there is less allocation overhead.  All data going through the
            channel passes through here at least once per frame (unless one is using a function eval effect)
        precision::Float64 - A somewhat arbitrary variable that determines how much memory to allocate to the
            virtual memory abstraction.  Allows higher fidelity rendering of effects assuming one has as a dense
            enough LED mapping that it exceeds the default 100 per channel. Higher precisions obviously cause more
            computation workload per frame.

    Inner Constructors:
        LEDChannel(
            map::Array{Tuple{LEDStrip{LEDChannel,LEDController}, Float64, Float64, AFArray{Float32,1}}, 1},
            precision::Int64
        )
            returns: A LEDChannel initialized with zeroed out virtual memory, uses given map and precision.

"""
mutable struct LEDChannel <: AbstractChannel
    map::Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64, AFArray{Float32, 1}},1}
    virtualmem::Array{UInt8,2}
    gpu_virtualmem::AFArray{Float32,2}
    precision::Float64
    function LEDChannel(map::Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64, AFArray{Float32, 1}},1}, precision::Int64)
        init_virtualmem = zeros(UInt8, precision, 3)
        gpu_virtualmem = AFArray(convert(Array{Float32}, init_virtualmem))
        return new(map, init_virtualmem, gpu_virtualmem, precision)
    end
end

"""
    getGPUIdxArray(strip::LEDStrip, x::Float64, y::Float64, precision::Real)::AFArray{Float32}

    Summary:
        Translates the 1:end indexes of the strip into pseudo-indexes for interpolation on the GPU

    Details:
        This function gets the virtual memory addresses corresponding to the real indices of the given
        strip.  Essentially it returns an array of equal size to the input strip, but each element in the
        returned array is the index that the LED at that position should reference in the GPU.  Note that
        these can be Float32s, this is intentional, since this data is then fed into in interpolation algorithm
        that doesn't mind non-integer indices.

    Arguments:
        strip::LEDStrip - The strip for which to get indexes
        x::Float64 - The start point of the strip in the LEDChannel
        y::Float64 - The end point of the strip in the LEDChannel
        precision - The precision of the LEDChannel being indexed

    Returns:
        AFArray{Float32, 1} (Array on GPU) of pseudo-indexes mapping to the LEDs in the LEDStrip
"""
function getGPUIdxArray(strip::LEDStrip, x::Float64, y::Float64, precision::Real)
    return AFArray(convert.(Float32, collect(linspace(x/100*precision, y/100*precision, size(strip,1)))))
end

"""
    LEDChannel(map::Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64},1}, precision::Float64)

    Summary:
        Convenience constructor for an LEDChannel that starts with the given (non-GPU) map and precision

    Details:
        This function is a convenience constructor for LEDChannel that first creates a map with the proper
        gpuIdxArrays initialized on the GPU, then initializes an LEDChannel with that new, valid map.  Prevents
        us from having to do the GPU setup manually with every initialization of an LEDChannel.  This is the
        primary constructor used to build other more convenient constructors.

    Returns:
        LEDChannel initialized with zeroed out virtual memory and the given map and precision.
"""
function LEDChannel(map::Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64},1}, precision::Float64)
    new_map = map(map) do x
        return (x[1:3]..., getGPUIdxArray(x[1:3]..., precision))
    end
    return LEDChannel(new_map, precision)
end

"LEDChannel overload to pass only the map and use a default precision"
LEDChannel(map::Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64, AFArray{Float32, 1}},1}) = LEDChannel(map, 100)

"LEDChannel overload to pass only the precision and use a default (empty) map"
LEDChannel(precision::Int64) = LEDChannel(Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64, AFArray{Float32, 1}},1}(0), precision)

"LEDChannel overload to pass no arguments at all and use a default map, precision, etc."
LEDChannel() = LEDChannel(Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64, AFArray{Float32, 1}},1}(0), 100)

"Implements getindex for LEDChannel in a maximally generic way"
getindex(channel::LEDChannel, idx::Any...) = getindex(channel.virtualmem, idx...)

"Implements setindex! for LEDChannel in a very generic way"
setindex!(channel::LEDChannel, val, idx::Any...) = setindex!(channel.virtualmem, val, idx...)

"Implements setindex! for LEDChannel for a specialized value type of RGB for convenience in writing effects"
function setindex!(channel::LEDChannel, val::ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}, idx...)
    channel.virtualmem[idx...,1] = val.r.i
    channel.virtualmem[idx...,2] = val.g.i
    channel.virtualmem[idx...,3] = val.b.i
    return channel.virtualmem[idx...,:]
end

"Implements setindex! for LEDChannel for an array type of RGB for convenience and performance in writing effects"
function setindex!(channel::LEDChannel, val::Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}}}, idx...)
    setindex!(channel, [getfield(v, f).i for v in val, f in fieldnames(eltype(val))], idx..., :)
end

"Implements a length function for LEDChannels"
length(channel::LEDChannel) = length(channel.virtualmem)

"Implements a size function for LEDChannels, the most basic overload"
size(channel::LEDChannel, x::Int) = size(channel.virtualmem, x)

"""
    Implements a more generic size function for LEDChannels that catches more
    at the expense of being less specialized (and thus slower, usually)
"""
size(channel::LEDChannel, x...) = size(channel.virtualmem, x...)

"""
    Implements an eachindex function for LEDChannels that allows easy iteration
    by returning a range of valid indicies
"""
eachindex(channel::LEDChannel) = 1:length(channel)

"Implements the lastindex function for LEDChannel, returns the last valid index"
lastindex(c::LEDChannel) = length(c)

"""
    update!(f:Function, channel::LEDChannel)

    Summary:
        Updates the representations of physical LEDs according to a continuous function over [0.0,100.0]

    Details:
        The update! functions are probably two of the most carefully implemented functions in the
        whole project because there is a potentially a huge amount of processing to be done and
        this function has to be performant because it is called at least once every frame of
        execution.

        This version takes in a function and evaluates that function at all of the virtual indexes
        that the LEDChannel mapping specifies.  It writes the results of that function evaluation
        to the LEDStrip memory (and by reference the LEDController memory), which then can be sent
        out to the real controller and rendered onto the real LEDs.

        Note: Since the function is the first argument, this supports Julia's "do syntax".

    Arguments:
        f::Function - A continuous function from 0.0 to 100.0 that outputs a color for the LED at
            virtual index
        channel::LEDChannel - The channel on which to render the functional effect.

    Returns:
        None
"""
function update!(f::Function, channel::LEDChannel)
    for m in channel.map
        indicies = linspace(m[2]/100*precision, m[3]/100*channel.precision, size(m[1],1))
        m[1] .= f.(indices)
    end
end

"""
    update!(channel::LEDChannel)

    Summary:
        Updates the representations of physical LEDs by interpolating the virtual memory

    Details:
        This function is the bread and butter of most of the existing effects because it
        most closely mimics programming real LEDs, just in this case they aren't actually real.
        It works by taking the current state of the virtual memory and transferring that to a
        pre-allocated array on the GPU.  Then it considers a strip and within that strip a single
        color channel.  Using the LEDChannel's map's gpuIdxArray to provide the pseudo-indexes of
        the strip, it interpolates for each channel.  Then it transfers the data back to main
        memory, after which it applies bounds that prevent the data from overflowing a UInt8.
        Lastly it rounds the whole array to UInt8 and assigns those values to the strip.

        This can be a rather large amount of work when there are many channels (i.e. 55) and large
        numbers of LEDs per channel (i.e. 397).  Those numbers mentioned just so happen to be the
        maximum currently possible with only one controller because UDP packet size limits me to
        21835 LEDs per controller.  Thus a decent bit of care was taken to make this function efficient.

        Note: Presently this function is not thread safe as it does not have any sort of lock on
        the gpu_virtualmem as it transfers that to the GPU, nor does it lock the output path, like
        the controller memory or the strip memory view.

    Arguments:
        channel::LEDChannel - The channel to interpolate to concrete LED values.

    Returns:
        None

"""
function update!(channel::LEDChannel)
    channel.gpu_virtualmem::AFArray{Float32,2} = AFArray(convert.(Float32, channel.virtualmem))
    for m in channel.map
        strip = m[1]
        for k in 1:size(strip, 2)
            strip.subArray[:,k] .= round.(UInt8, min.( max.( Array( approx1(channel.gpu_virtualmem[:,k], m[4], AF_INTERP_CUBIC_SPLINE, 0f0) ), 0.0f0), 255.0f0) )
            #setindex!(strip.subArray, tmp, 1:size(strip.subArray,1), k)
        end
    end
end

"""
    push!(channel::LEDChannel, val::LEDStrip{T,N}) where {T<:AbstractChannel, N<:AbstractChannel}

    Summary:
        This function adds an LEDStrip to the LEDChannel's map in the default configuration

    Details:
        Adds an LEDStrip to the LEDChannel's map in the default configuration which is that the
        LEDStrip spans the whole channel range, [0.0,100.0].  If it already exists in the map,
        nothing is changed.

    Arguments:
        channel::LEDChannel - The channel to which to add the strip
        val::LEDStrip - The strip to add to the channel's map.

    Returns:
        LEDStrip which was added to the channel

"""
function push!(channel::LEDChannel, val::LEDStrip{T,N}) where {T<:AbstractChannel, N<:AbstractController}
    return push!(channel, val, 0.0, 100.0)
end

"""
    push!(channel::LEDChannel, val::LEDStrip{T,N}, startLoc::Real, endLoc::Real) where {T<:AbstractChannel, N:<AbstractController}

    Summary:
        This function adds an LEDStrip to the LEDChannel's map in the specified mapping

    Details:
        Adds an LEDStrip to the LEDChannel's map with the specified start and end location.
        This function does nothing if the value is already present in the LEDChannel's map.

    Arguments:
        channel::LEDChannel - The channel to which to add the strip
        val::LEDStrip - The strip to add to the channel's map
        startLoc::Real - A number representing the beginning of the strip's pseudo-indexes on the channel
        endLoc::Real - A number representing the end of the strip's pseudo-indexes on the channel

    Returns:
        LEDStrip which was added or already present on the channel

"""
function push!(channel::LEDChannel, val::LEDStrip{T,N}, startLoc::Real, endLoc::Real) where {T<:AbstractChannel, N<:AbstractController}
    if !((startLoc >= 0.0 && startLoc <= 100.0) && (endLoc >= 0.0 && endLoc <= 100.0))
        error("Strips must be assigned to locations between 0.0 and 100.0 on a channel")
    end
    if !(val in getindex.(channel.map, 3))
        return push!(channel.map, (val, startLoc, endLoc, getGPUIdxArray(val, startLoc, endLoc, channel.precision)))
    else
        return val
    end
end

"""
    LEDChannel(channel1::LEDChannel, channel2::LEDChannel, offset::Float64)

    Summary:
        Another constructor overload for LEDChannel that allows merging and concatenating channels

    Details:
        Another constructor overload for LEDChannel that creates a new channel by merging
        two existing channels with a given offset.  So if you wanted to set them up to mirror
        each other, the offset would be zero, if you wanted to concatenate them, you would specify
        an offset of 100.0.  Negative offsets are also supported in case one wants to overlay the
        second channel to the left of the first channel.

    Arguments:
        channel1::LEDChannel - the first LEDChannel to merge or concatenate
        channel2::LEDChannel - the second LEDChannel to merge or concatenate
        offset::Float64 - A number representing how much to offset the channels when remapping,
            the resultant channel will be renomalized against the new overall width so
            values from ~(-1000.0) to ~(+1000.0) are allowable.

    Returns:
        LEDChannel the resultant channel from the merge.
"""
function LEDChannel(channel1::LEDChannel, channel2::LEDChannel, offset::Float64)
    new_virtual_space = 100 + abs(offset)
    multiplier = 100/new_virtual_space
    old_map1 = deepcopy(channel1.map)
    old_map2 = deepcopy(channel2.map)
    new_map = Array{Tuple{LEDStrip{LEDChannel, LEDController}, Float64, Float64, AFArray{Float32, 1}},1}(length(channel1.map)+length(channel2.map))
    new_precision = max(channel1.precision, channel2.precision)
    indexed = 0
    offset1 = 0
    offset2 = 0
    if offset < 0
        offset1 = abs(offset)
    else
        offset2 = abs(offset)
    end
    for i in eachindex(channel1.map)
        indexed = i
        new_map[i] = old_map1[i]
        @. new_map[i][2:3] = (new_map[i][2:3]+offset1)*multiplier
        new_map[i][4] = getGPUIdxArray()
    end
    for i in eachindex(channel2.map)
        new_map[i+indexed] = old_map2[i]
        @. new_map[i+indexed][2:3] = (new_map[i+indexed][2:3]+offset2)*multiplier
        new_map[i+indexed][4] = getGPUIdxArray(new_map[i+indexed][1:3]..., new_precision)
    end
    return LEDChannel(new_map, new_precision)
end

"""
    LEDChannel(channels::Array{LEDChannels}, offsets::Array{Float64})

    Summary:
        Takes an Array of LEDChannels and merges them according to an array of offsets

    Details:
        A constructor overload that creates a new channel from an array of channels
        and an array of offsets (the first of which is the offset of the first channel
        and so on).  The merging process normalizes everything against the channel offset
        furthest to the left (furthest negative, even!).

    Arguments:
        channels::Array{LEDChannel} - An array of channels to be merged together
        offsets::Array{Float64} - An array (equal length as channels) of the offsets
            to be used to merge the channels

    Returns:
        LEDChannel the resultant channel from the merge.
"""
function LEDChannel(channels::Array{LEDChannel}, offsets::Array{Float64})
    root_offset = 0
    max_offset = maximum(offsets)
    min_offset = minimum(offsets)
    if min_offset < 0
        root_offset -= min_offset
        offsets .-= min_offset
    end
    new_virtual_space = (max_offset+100) - min_offset
    multiplier = 100/new_virtual_space
    new_precision = maximum(getfield.(maps, :precision))
    maps = deepcopy.(getfield.(channels, :map))
    for i in eachindex(maps)
        if i == 1
            map!(maps[i]) do x
                tmp = (x[1], ((x[2:3]+root_offset)*multiplier)..., x[4])
                tmp[4] = getGPUIdxArray(tmp[1:3]..., new_precision)
                return tmp
            end
        else
            map!(maps[i]) do x
                tmp = (x[1], ((x[2:3]+offsets[i-1])*multiplier)..., x[4])
                tmp[4] = getGPUIdxArray(tmp[1:3]..., new_precision)
                return tmp
            end
        end
    end
    return LEDChannel(vcat(maps...), new_precision)
end

"""
    LEDArray

    Summary:
        A convenience bundle that contains all of the relavent data to update
        the LEDs.

    Details:
        This structure bundles all of the important structures used by Irradiance
        into one large structure for convenience.

    Fields:
        controllers::Array{LEDController,1} - All of the active LEDControllers
        inactive_channels::Array{LEDChannel,1} - All of the inactive LEDChannels that
            are worth keeping around, such as "root channels" which only hold one
            LEDStrip mapped across the whole channel.
        channels::Array{LEDChannel,1} - All of the active channels that are being rendered to
        strips::Array{LEDStrip{LEDController,LEDChannel},1} - All of the currently known strips

    Inner Constructors:
        (default) LEDArray(
            controllers::Array{LEDController},
            inactive_channels::Array{LEDChannel,1},
            channels::Array{LEDChannel,1},
            strips::Array{LEDStrip{LEDController, LEDChannel},1}
            )
            Returns:
                LEDArray with specified values exactly as passed

"""
mutable struct LEDArray
    controllers::Array{LEDController,1}
    inactive_channels::Array{LEDChannel,1}
    channels::Array{LEDChannel,1}
    strips::Array{LEDStrip{LEDChannel, LEDController},1}
end

"""
    LEDStrip(name::String, channel::LEDChannel, controller::LEDController, startAddr::Int, endAddr::Int)

    Summary:
        Concrete constructor free of the generic types previously needed to get around lookahead issues.

    Details:
        This is the concrete constructor for LEDStrip that uses the concrete implementations of LEDController
        and LEDChannel. Using as concrete of types as possible is good performance in Julia, so this was worth
        the parameterization etc. to get here.

    Arguments:
        name::String - A name for the strip
        channel::LEDChannel - The channel to attach this strip to
        controller::LEDController - The controller to attach this strip to
        startAddr::Int - The first index on the controller that this strip represents
        endAddr::Int - The last index on the controller that this strip represents

    Returns:
        LEDStrip{LEDChannel, LEDController} with the specified arguments set, and a reference to the relevant
            controller memory.
"""
function LEDStrip(name::String, channel::LEDChannel, controller::LEDController, startAddr::Int, endAddr::Int)
    return LEDStrip{LEDChannel, LEDController}(name, channel, controller, startAddr, endAddr)
end


text:
Base.@kwdef struct CHSettings
    host::String = "localhost"
    port::Int = 9000
    database::String = ""
    username::String = "default"
    password::String = ""
    connection_timeout::Int = DBMS_DEFAULT_CONNECT_TIMEOUT
    max_insert_block_size::Int = DBMS_DEFAULT_MAX_INSERT_BLOCK
    send_buffer_size::Int = DBMS_DEFAULT_BUFFER_SIZE
end

mutable struct ClickHouseSock
    io ::Union{IO, Nothing}
    settings ::CHSettings
    cond ::ReentrantLock
    busy ::Bool
    server_name ::String
    server_rev ::Int
    server_timezone ::Union{String, Nothing}

    function ClickHouseSock(io, settings = CHSettings())
        return new(
            io,
            settings,
            ReentrantLock(),
            false,
            "", 0, nothing
        )
    end
end

"""
    @guarded(sock::ClickHouseSock, expr)

    Run `expr` thread-safe under lock of `sock.cond`.
"""
macro guarded(sock, expr)
    quote
        lock($(esc(sock)).cond)
        local res = try
            $(esc(expr))
        catch e
            unlock($(esc(sock)).cond)
            rethrow(e)
        end
        unlock($(esc(sock)).cond)
        res
    end
end

is_connected(sock::ClickHouseSock) = @guarded sock !isnothing(sock.io)
is_busy(sock::ClickHouseSock) = @guarded sock sock.busy

set_busy!(sock::ClickHouseSock, value::Bool) =
                        @guarded sock sock.busy = value

"""
    @using_socket(sock::ClickHouseSock, expr)

    Set `sock.busy` status and run `expr`.
    Raises an exception if `sock` is not connected or already busy.
    If an exception occurs during the execution of an expression and this exception
    is not a server exception, the socket will disconnect.
"""
macro using_socket(sock, expr)
    quote
        @guarded $(esc(sock)) begin
            isnothing($(esc(sock)).io) && error("ClickHouseSock not connected")
            $(esc(sock)).busy && error("ClickHouseSock is busy")
            $(esc(sock)).busy = true
        end
        local res = try
            $(esc(expr))
        catch e
            typeof(e) == ClickHouseServerException ?
                    set_busy!($(esc(sock)), false) :
                    close($(esc(sock)))
            rethrow(e)
        end
        set_busy!($(esc(sock)), false)
        res
    end
end


text:
# WARNING: almost certainly outdated.
# https://github.com/JuliaGaussianProcesses/Stheno.jl/issues/223

# Set up the environment to run this example. Make sure you're within the folder that this
# file lives in.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Stheno
using Plots; pyplot();
using Random; Random.seed!(4);
using Flux
using Zygote
using DelimitedFiles
using Statistics

# read AirPass data
data = readdlm("AirPassengers.csv", ',')
year = data[2:end,2]; passengers = data[2:end,3];
# Split the data into training and testing data
oxtrain = year[year.<1958]; oytrain = passengers[year.<1958];
oxtest = year[year.>=1958]; oytest = passengers[year.>=1958];

#data preprocessing
## standardize X and y
xtrain_mean = mean(oxtrain)
ytrain_mean = mean(oytrain)
xtrain_std = std(oxtrain)
ytrain_std = std(oytrain)
xtrain = @. (oxtrain-xtrain_mean)/xtrain_std
ytrain = @. (oytrain-ytrain_mean)/ytrain_std

xtest = @. (oxtest-xtrain_mean)/xtrain_std
ytest = @. (oytest-ytrain_mean)/ytrain_std

## input data
Xtrain = reshape(xtrain, 1, length(xtrain));
Xtest = reshape(xtest, 1, length(xtest));
Year = hcat(Xtrain, Xtest);
Passengers = vcat(ytrain, ytest)

# kernel parameter initialization
function median_distance_local(x)
    n = length(x)
    dist = []
    for i in 1:n
        for j in i:n
            push!(dist, abs(x[j]-x[i]))
        end
    end
    median(dist)
end
l = median_distance_local(xtrain)

# NKN
# NOTE: the numerical values are in log-scale
# therefore, stretch(Linear(), 0.0, exp) ≡ stretch(Linear(), 1.0)
# `exp` is a constraint added to the kernel parameters to maintain
# their positiveness during the optimisation.
primitives = Primitive(
    stretch(Linear(), 0.0, x->exp(-x)),
    stretch(PerEQ(log(l), exp), log(l), x->exp(-x)),
    stretch(SEKernel(), log(l/4.0), x->exp(-x)),
    stretch(RQ(log(0.2), exp), log(2.0*l), x->exp(-x)),
    stretch(Linear(), 0.0, x->exp(-x)),
    stretch(RQ(log(0.1), exp), log(l), x->exp(-x)),
    stretch(SEKernel(), log(l), x->exp(-x)),
    stretch(PerEQ(log(l/4.0), exp), log(l/4.0), x->exp(-x))
)
nn = Chain(
    LinearLayer(8, 8),
    product,
    LinearLayer(4, 4),
    product,
    LinearLayer(2, 1),
)
nkn = NeuralKernelNetwork(primitives, nn)

# build GP model
σ²_n = 0.1
gp = GP(nkn, GPC())
gp_Xtrain = gp(ColVecs(Xtrain), σ²_n)
ps = params(nkn)

# optimize
using Flux.Optimise: update!

optimizer = ADAM(0.001)
loss = []
for i in 1:5000
    ll = .-logpdf(gp_Xtrain, ytrain)
    push!(loss, ll)
    if i==1 || i%100 == 0
        @info "step=$i, loss=$ll"
    end
    gs = gradient(()->.-logpdf(gp_Xtrain, ytrain), ps)
    for p in ps
        update!(optimizer, p, gs[p])
    end
end

png(plot(loss), "loss.png")

# predict
function predict(X, Xtrain, ytrain)
    noisy_prior = gp(ColVecs(Xtrain), σ²_n)
    posterior = gp | Obs(noisy_prior, ytrain)
    posterior(ColVecs(X))
end

posterior = predict(Year, Xtrain, ytrain)
post_marginals = marginals(posterior)
pred_y = mean.(post_marginals)
var_y = std.(post_marginals)
pred_oy = @. pred_y*ytrain_std+ytrain_mean
pred_oσ = @. var_y*ytrain_std

plt = plot(xlabel="Year", ylabel="Airline Passenger number", legend=true)
plot!(plt, year, pred_oy;
    ribbons=3*pred_oσ,
    title="Time series prediction",
    label="95% predictive confidence region",
)
scatter!(plt, oxtest, oytest, label="Observations(test)", color=:red)
scatter!(plt, oxtrain, oytrain, label="Observations(train)", color=:black)
png(plt, "predict.png")


text:
using Pkg
Pkg.activate(".")
pkg"dev .. ../CairoMakie ../GLMakie ../WGLMakie"
Pkg.instantiate()
Pkg.precompile()

using NodeJS
run(`$(npm_cmd()) install highlight.js`)
run(`$(npm_cmd()) install lunr`)
run(`$(npm_cmd()) install cheerio`)

using Franklin
using Documenter: deploydocs, deploy_folder, GitHubActions
using Gumbo
using AbstractTrees

cfg = GitHubActions() # this should pick up all details via GHA environment variables

repo = "github.com/JuliaPlots/Makie.jl.git"
push_preview = true

deploydecision = deploy_folder(cfg; repo, push_preview,
    devbranch = "master",
    devurl = "dev",
)

@info "Setting PREVIEW_FRANKLIN_WEBSITE_URL to $repo"
ENV["PREVIEW_FRANKLIN_WEBSITE_URL"] = repo

"""
Converts the string `s` which might be an absolute path,
to a relative path, relative to the current location `here`
"""
function make_relative(s, here)
    if !startswith(s, "/")
        return s
    end

    there = s[2:end]

    if here == ""
        return there
    end

    hereparts = split(here, "/")
    thereparts = split(there, "/")

    # tutorials/layout-tutorial   tutorials/
    closest_common = 0
    for i in 1:min(length(thereparts), length(hereparts))
        if hereparts[i] == thereparts[i]
            closest_common = i
        else
            break
        end
    end
    n_up = length(hereparts) - closest_common
    therepart = join(thereparts[closest_common+1:end], "/")
    if n_up == 0
        therepart
    else
        up_part = join((".." for i in 1:n_up), "/")
        if therepart == ""
            up_part
        else
            up_part * "/" * therepart
        end
    end
end

"""
Replaces all absolute links in all html files in the __site folder with
relative links.
"""
function make_links_relative()

    old = pwd()
    try
        cd("__site")
        for (root, dirs, files) in walkdir(".")
            path = join(splitpath(root)[2:end], "/")

            html_files = filter(endswith(".html"), files)
            for file in html_files
                s = read(joinpath(root, file), String)

                html = parsehtml(s)

                for e in PreOrderDFS(html.root)
                    if (e isa HTMLElement{:script} ||
                        e isa HTMLElement{:img} ||
                        e isa HTMLElement{:video}) && haskey(e.attributes, "src")

                        link = e.attributes["src"]
                        e.attributes["src"] = make_relative(link, path)

                    elseif (e isa HTMLElement{:link} ||
                            e isa HTMLElement{:a}) && haskey(e.attributes, "href")

                        link = e.attributes["href"]
                        e.attributes["href"] = make_relative(link, path)

                    elseif e isa HTMLElement{:form} && haskey(e.attributes, "action")

                        link = e.attributes["action"]
                        e.attributes["action"] = make_relative(link, path)
                    end
                end

                open(joinpath(root, file), "w") do f
                    print(f, html)
                end
            end
        end
    finally
        cd(old)
    end
end

serve(single=true, cleanup=false)
lunr()
optimize(minify=false, prerender=false)

# by making all links relative, we can forgo the `prepath` setting of Franklin
# which means that files in some `vX.Y.Z` subfolder which happens to be `stable`
# at the time, link relatively within `stable` so that users don't accidentally
# copy & paste versioned links if they started out on `stable`
@info "Rewriting all absolute links as relative"
make_links_relative()

deploydocs(;
    repo,
    push_preview,
    target = "__site",
)


text:
using LinearAlgebra, Test
using TensorToolbox

include("doctest.jl") #check for erros in README examples

include("tensorTest.jl")
include("ttensorTest.jl")
include("ktensorTest.jl")
include("dimtreeTest.jl")
include("htensorTest.jl")
include("TTtensorTest.jl")

println("\n\n")


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_antdchord

"""
    ''_antdchord(;kwargs...)

An AntdChord component.

Keyword arguments:
- `id` (String; optional)
- `appendPadding` (Real | Array of Reals; optional)
- `autoFit` (Bool; optional)
- `className` (String; optional)
- `data` (Array of Dicts; optional)
- `edgeStyle` (optional): . edgeStyle has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `height` (Real; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-CN', 'en-US'; optional)
- `meta` (optional)
- `nodePaddingRatio` (Real; optional)
- `nodeStyle` (optional): . nodeStyle has the following type: lists containing elements 'func'.
Those elements have the following types:
  - `func` (String; optional)
- `nodeWidthRatio` (Real; optional)
- `padding` (Real | Array of Reals | a value equal to: 'auto'; optional)
- `renderer` (a value equal to: 'canvas', 'svg'; optional)
- `sourceField` (String; optional)
- `style` (Dict; optional)
- `targetField` (String; optional)
- `weightField` (String; optional)
- `width` (Real; optional)
"""
function ''_antdchord(; kwargs...)
        available_props = Symbol[:id, :appendPadding, :autoFit, :className, :data, :edgeStyle, :height, :loading_state, :locale, :meta, :nodePaddingRatio, :nodeStyle, :nodeWidthRatio, :padding, :renderer, :sourceField, :style, :targetField, :weightField, :width]
        wild_props = Symbol[]
        return Component("''_antdchord", "AntdChord", "feffery_antd_charts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-antd-charts

max_issues_repo_name: CNFeffery/feffery-antd-charts

max_forks_repo_name: CNFeffery/feffery-antd-charts


text:
###ODE for the FitzHugh-Nagumo model: http://www.scholarpedia.org/article/FitzHugh-Nagumo_model
### t = timepoints to evaluate the ODE at
### y = the model's state variables (Membrane Potential,Refractory Variable)
### ydot = the derivate values
### paras = the equation parameters (a,b,c)
function fitzhughnagumoode(t,y,ydot,paras)

    a = paras[1]
    b = paras[2]
    c = paras[3]

    ydot[1] = c*(y[1]-(y[1]^3)/3+y[2])
    ydot[2] = -(y[1]-a+b*y[2])/c

    nothing
end

###Convenience function to create a FitzHugh-Nagumo ODEModel
### y0: initial conditions of the dynamic system
### variance: the variance of the noise model
### paraminit: the default and boundary values of the parameters
function fitzhughnagumomodel(y0::AbstractVector,variance::AbstractVector,paraminit...)
    modelparameters = parameters([:a,:b,:c],paraminit...)
    noisemodel = noise(:gaussian,variance)
    measurementdata = data(:array,fitzhughnagumodataset1()...)
    model(:ode,modelparameters,measurementdata,noisemodel,fitzhughnagumoode,y0,2,[1,2];name = "FitzHugh-Nagumo")
end


text:
# This is a list of transitive Frobenius groups G with abelian Frobenius
# kernel K and solvable H = G/K and degree(G) <= 30.
#
# Frob := []l;
# for i in [1..30] do
#   for j in [1..NumberOfTransitiveGroups(i)] do
#     G := TransitiveGroup(i, j);
#     fl := IsFrobenius(G);
#     if not fl then
#       continue;
#     end if;
#     i, j, GroupName(G);
#     Append(~Frob, <i, jj>);
#   end for;
# end for
#
# for ID in Frob do
#   G := TransitiveGroup(ID[1], ID[2]);
#   K := FittingSubgroup(G);
#   sid := IdentifyGroup(G);
#   H := Complements(G, K);
#   assert #H eq 1; H := H[1];
#   if not IsSolvable(H) or not IsAbelian(K) then
#     continue;
#   end if;
#   <sid, IdentifyGroup(H), InvariantFactors(AbelianGroup(K))>;
#   end for;
#
# For each Frobenius group we store we have a pair
# SID -> (KID, H)
# where
#  - SID is the small group ID
#  - KID is the small group ID of G/K
#  - H is the list of invariant factors of K.

global _frobenius_groups = Dict{Tuple{Int, Int}, Tuple{Tuple{Int, Int}, Vector{Int}}}(
  (600, 150) => ((24, 3), [5, 5]),
  (342, 7) => ((18, 2), [19]),
  (39, 1) => ((3, 1), [13]),
  (300, 24) => ((12, 2), [5, 5]),
  (600, 149) => ((24, 2), [5, 5]),
  (406, 1) => ((14, 2), [29]),
  (203, 1) => ((7, 1), [29]),
  (702, 47) => ((26, 2), [3, 3, 3]),
  (56, 11) => ((7, 1), [2, 2, 2]),
  (34, 1) => ((2, 1), [17]),
  (50, 1) => ((2, 1), [25]),
  (150, 6) => ((6, 2), [5, 5]),
  (54, 7) => ((2, 1), [3, 9]),
  (10, 1) => ((2, 1), [5]),
  (18, 1) => ((2, 1), [9]),
  (26, 1) => ((2, 1), [13]),
  (42, 1) => ((6, 2), [7]),
  (80, 49) => ((5, 1), [2, 2, 2, 2]),
  (136, 12) => ((8, 1), [17]),
  (46, 1) => ((2, 1), [23]),
  (156, 7) => ((12, 2), [13]),
  (57, 1) => ((3, 1), [19]),
  (21, 1) => ((3, 1), [7]),
  (38, 1) => ((2, 1), [19]),
  (78, 1) => ((6, 2), [13]),
  (48, 3) => ((3, 1), [4, 4]),
  (351, 12) => ((13, 1), [3, 3, 3]),
  (84, 11) => ((3, 1), [2, 14]),
  (22, 1) => ((2, 1), [11]),
  (300, 23) => ((12, 1), [5, 5]),
  (14, 1) => ((2, 1), [7]),
  (30, 3) => ((2, 1), [15]),
  (42, 5) => ((2, 1), [21]),
  (100, 11) => ((4, 1), [5, 5]),
  (36, 9) => ((4, 1), [3, 3]),
  (55, 1) => ((5, 1), [11]),
  (114, 1) => ((6, 2), [19]),
  (110, 1) => ((10, 2), [11]),
  (812, 7) => ((28, 2), [29]),
  (48, 50) => ((3, 1), [2, 2, 2, 2]),
  (54, 14) => ((2, 1), [3, 3, 3]),
  (100, 12) => ((4, 1), [5, 5]),
  (72, 39) => ((8, 1), [3, 3]),
  (6, 1) => ((2, 1), [3]),
  (52, 3) => ((4, 1), [13]),
  (12, 3) => ((3, 1), [2, 2]),
  (253, 1) => ((11, 1), [23]),
  (200, 40) => ((8, 1), [5, 5]),
  (171, 3) => ((9, 1), [19]),
  (68, 3) => ((4, 1), [17]),
  (200, 44) => ((8, 4), [5, 5]),
  (600, 148) => ((24, 1), [5, 5]),
  (272, 50) => ((16, 1), [17]),
  (54, 1) => ((2, 1), [27]),
  (506, 1) => ((22, 2), [23]),
  (72, 41) => ((8, 4), [3, 3]),
  (240, 191) => ((15, 1), [2, 2, 2, 2]),
  (58, 1) => ((2, 1), [29]),
  (50, 4) => ((2, 1), [5, 5]),
  (18, 4) => ((2, 1), [3, 3]),
  (116, 3) => ((4, 1), [29]),
  (100, 3) => ((4, 1), [25]),
  (20, 3) => ((4, 1), [5]),
  (75, 2) => ((3, 1), [5, 5])
)

function primitive_frobenius_extensions(::FlintRationalField, id::Tuple{Int, Int}, B::fmpz; only_real::Bool = false, only_non_real::Bool = false)
  @req haskey(_frobenius_groups, id) "id ($id) must be small group id of a transitive Frobenius group " *
                                     "with abelian Frobenius kernel and solvable quotient and degree " *
                                     "bounded by 30."
  @req !(only_real && only_non_real) "Either \"only_real\" or \"only_non_real\" must be set"
  sid, invfac = _frobenius_groups[id]
  R = ArbField(64, cached = false)
  res = AnticNumberField[]
  maxB = fmpz(1)
  reldeg = prod(invfac)
  k = 1
  s = (reldeg - 1)//sid[1] # (|F| - 1)/|H|
  #Ms = abelian_fields(QQ, [l], max(fmpz(1), upper_bound(fmpz, R(B)^(1//s))))
  Ms = fields(sid[1], sid[2], max(fmpz(1), upper_bound(fmpz, R(B)^(1//s))))
  ll = length(Ms)
  target_id = id

  for M in Ms
    k += 1
    Mnf = number_field(M)
    Mabs = Mnf
    dM = abs(discriminant(maximal_order(Mabs)))

    newB = upper_bound(fmpz, (B//R(dM)^s)^sid[1] * dM^reldeg)

    # If I want only real degree l fields, then the normal closure N can be anything
    #
    # If I want only non-real degree l fields, then the normal closure must not be real,
    # hence must be totally complex.
    kk = 1
    Ns = abelian_normal_extensions(Mabs, invfac, newB, only_complex = only_non_real)

    for N in Ns
      kk += 1
      Nabs, = absolute_simple_field(number_field(N, using_norm_relation = true))
      _A = absolute_automorphism_group(N)
      _A = closure(_A, *)
      id, = find_small_group(generic_group(_A, *)[1])
      if id != target_id
        continue
      end

      A = automorphisms(Nabs)
      N, AtoN, NtoA = generic_group(A, *)

      cur = nothing

      for (K, KtoN) in subgroups(N, normal = true)
        _,KK, = find_small_group(K)
        if KK.isnilpotent == 1 && (cur === nothing || order(KK) > order(cur[1]))
          cur = K, KtoN
        end
        if cur !== nothing && order(cur[1]) == reldeg
          break
        end
      end
      K, KtoN = cur
      @assert order(K) == reldeg
      # I found a kernel
      # Now find a complement
      comp = (K, KtoN)
      KinN = [ KtoN(k) for k in K]
      for (H, HtoN) in subgroups(N, order = sid[1])
        HinN = [HtoN(h) for h in H]
        if length(unique!([ k * h for k in KinN for h in HinN ])) == order(N)
          comp = (H, HtoN)
          break
        end
      end

      H, HtoN = comp
      @assert order(H) == sid[1]
      # trigger better automorphism group computation
      lll(maximal_order(Nabs))

      target_field,  = fixed_field1(Nabs, [ NtoA[HtoN(k)] for k in gens(H) ])
      @assert abs(discriminant(maximal_order(target_field))) <= B
      # In the only_non_real case, we already have the right fields,
      if (!only_real && !only_non_real) || (only_non_real && !istotally_real(target_field)) || (only_real && istotally_real(target_field))
        push!(res, target_field)
      end
    end
  end
  return res
end


text:
# Aliases (types)
const FloatVector  = Array{Float64,1};
const FloatArray   = Array{Float64};
const SymMatrix    = Symmetric{Float64,Array{Float64,2}};
const DiagMatrix   = Diagonal{Float64,Array{Float64,1}};
const JVector{T}   = Array{Union{Missing, T}, 1};
const JArray{T, N} = Array{Union{Missing, T}, N};

# Kalman structures

abstract type KalmanSettings end

"""
    ImmutableKalmanSettings(...)

Define an immutable structure that includes all the Kalman filter and smoother inputs.

# Model
The state space model used below is,

``Y_{t} = B*X_{t} + e_{t}``

``X_{t} = C*X_{t-1} + v_{t}``

Where ``e_{t} ~ N(0, R)`` and ``v_{t} ~ N(0, V)``.

# Arguments
- `Y`: observed measurements (`nxT`)
- `B`: Measurement equations' coefficients
- `R`: Covariance matrix of the measurement equations' error terms
- `C`: Transition equations' coefficients
- `V`: Covariance matrix of the transition equations' error terms
- `X0`: Mean vector for the states at time t=0
- `P0`: Covariance matrix for the states at time t=0
- `n`: Number of series
- `T`: Number of observations
- `m`: Number of latent states
- `compute_loglik`: Boolean (true for computing the loglikelihood in the Kalman filter)
- `store_history`: Boolean (true to store the history of the filter and smoother)
"""
struct ImmutableKalmanSettings <: KalmanSettings
    Y::Union{FloatArray, JArray{Float64}}
    B::FloatArray
    R::SymMatrix
    C::FloatArray
    V::SymMatrix
    X0::FloatVector
    P0::SymMatrix
    n::Int64
    T::Int64
    m::Int64
    compute_loglik::Bool
    store_history::Bool
end

# ImmutableKalmanSettings constructor
function ImmutableKalmanSettings(Y::Union{FloatArray, JArray{Float64}}, B::FloatArray, R::SymMatrix, C::FloatArray, V::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B,2);
    X0 = zeros(m);
    P0 = Symmetric(reshape((I-kron(C, C))\V[:], m, m));

    # Return ImmutableKalmanSettings
    return ImmutableKalmanSettings(Y, B, R, C, V, X0, P0, n, T, m, compute_loglik, store_history);
end

# ImmutableKalmanSettings constructor
function ImmutableKalmanSettings(Y::Union{FloatArray, JArray{Float64}}, B::FloatArray, R::SymMatrix, C::FloatArray, V::SymMatrix, X0::FloatVector, P0::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B,2);

    # Return ImmutableKalmanSettings
    return ImmutableKalmanSettings(Y, B, R, C, V, X0, P0, n, T, m, compute_loglik, store_history);
end

"""
    MutableKalmanSettings(...)

Define a mutable structure identical in shape to ImmutableKalmanSettings.

See the docstring of ImmutableKalmanSettings for more information.
"""
mutable struct MutableKalmanSettings <: KalmanSettings
    Y::Union{FloatArray, JArray{Float64}}
    B::FloatArray
    R::SymMatrix
    C::FloatArray
    V::SymMatrix
    X0::FloatVector
    P0::SymMatrix
    n::Int64
    T::Int64
    m::Int64
    compute_loglik::Bool
    store_history::Bool
end

# MutableKalmanSettings constructor
function MutableKalmanSettings(Y::Union{FloatArray, JArray{Float64}}, B::FloatArray, R::SymMatrix, C::FloatArray, V::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B,2);
    X0 = zeros(m);
    P0 = Symmetric(reshape((I-kron(C, C))\V[:], m, m));

    # Return MutableKalmanSettings
    return MutableKalmanSettings(Y, B, R, C, V, X0, P0, n, T, m, compute_loglik, store_history);
end

# MutableKalmanSettings constructor
function MutableKalmanSettings(Y::Union{FloatArray, JArray{Float64}}, B::FloatArray, R::SymMatrix, C::FloatArray, V::SymMatrix, X0::FloatVector, P0::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B,2);

    # Return MutableKalmanSettings
    return MutableKalmanSettings(Y, B, R, C, V, X0, P0, n, T, m, compute_loglik, store_history);
end

"""
    KalmanStatus(...)

Define an mutable structure to manage the status of the Kalman filter and smoother.

# Arguments
- `t`: Current point in time
- `loglik`: Loglikelihood
- `X_prior`: Latest a-priori X
- `X_post`: Latest a-posteriori X
- `P_prior`: Latest a-priori P
- `P_post`: Latest a-posteriori P
- `history_X_prior`: History of a-priori X
- `history_X_post`: History of a-posteriori X
- `history_P_prior`: History of a-priori P
- `history_P_post`: History of a-posteriori P
"""
mutable struct KalmanStatus
    t::Int64
    loglik::Union{Float64, Nothing}
    X_prior::Union{FloatVector, Nothing}
    X_post::Union{FloatVector, Nothing}
    P_prior::Union{SymMatrix, Nothing}
    P_post::Union{SymMatrix, Nothing}
    history_X_prior::Union{Array{FloatVector,1}, Nothing}
    history_X_post::Union{Array{FloatVector,1}, Nothing}
    history_P_prior::Union{Array{SymMatrix,1}, Nothing}
    history_P_post::Union{Array{SymMatrix,1}, Nothing}
end

# KalmanStatus constructors
KalmanStatus() = KalmanStatus(0, [nothing for i=1:9]...);


text:
function outputsizes(fs::Tuple, x)
    res = Flux.outputsize(first(fs), x)
    return (res, outputsizes(Base.tail(fs), res)...)
end
outputsizes(::Tuple{}, x) = ()
outputsizes(c::Chain, x) = outputsizes(c.layers, x)

# trim off unused network parts
get_encoding_network_part(m::ResNet) = Chain(Chain(identity), m.layers[1][1:end-1]...)
get_encoding_network_part(m::VGG) = m.layers[1][1:end-1]

function encoder(m)
    encoder = get_encoding_network_part(m)
    os = outputsizes(encoder, (512,512,3,1))  # snoop intermediate output sizes
    os = map(x->x[1], os)  # take only the width
    os = (os..., 0)  # the below loop checks for differences in size, so we have to set an endpoint to also catch the last layers

    layers = []
    last_size = os[1]
    last = 1
    for i in 1:length(os)
        current_size = os[i]
        if current_size != last_size
            push!(layers, deepcopy(encoder[last:i-1]))  # actually I don't know if the deepcopy is needed, but just to be sure...
            last = i
            last_size = current_size
        end
    end
    return Chain(layers...)
end

function encoder_channels(m, input_channels=3)
    os = outputsizes(m, (512, 512, input_channels, 1))
    return map(x->x[3], os)
end


text:
# Copyright (c) 2021, Rolfe Power

"""
    GRAVITATIONAL_CONSTANT

Universal constant of gravitation. [km^3 / (kg . s^2)]

Source: [NIST](https://physics.nist.gov/cgi-bin/cuu/Value?bg|search_for=universal_in!)
"""
const GRAVITATIONAL_CONSTANT = 6.6743e-20

"""
    ASTRONOMICAL_UNIT

Astronomical unit as defined by the IAU. [km]

Source: [IAU](https://www.iau.org/static/resolutions/IAU2012_English.pdf)
"""
const ASTRONOMICAL_UNIT = 149_597_870.700

"""
    EARTH_STANDARD_GRAVITY

Standard gravitational acceleration at surface of Earth. [km/s^2]

Source: [NIST](https://physics.nist.gov/cgi-bin/cuu/Value?gn|search_for=adopted_in!)
"""
const EARTH_STANDARD_GRAVITY = 9.80665e-03

"""
    PROTON_MASS

Mass of the proton. [kg]

Source: [NIST](https://physics.nist.gov/cgi-bin/cuu/Value?mp|search_for=atomnuc!)
"""
const PROTON_MASS = 1.67262192369e-27

"""
    NEUTRON_MASS

Mass of the neutron. [kg]

Source: [NIST](https://physics.nist.gov/cuu/Constants/Table/allascii.txt)
"""
const NEUTRON_MASS = 1.67492749804e-27

"""
    ELECTRON_MASS

Mass of the electron. [kg]

Source: [NIST](https://physics.nist.gov/cgi-bin/cuu/Value?mp|search_for=atomnuc!)
"""
const ELECTRON_MASS = 9.1093837015e-31

"""
    SPEED_OF_LIGHT

Speed of light in a vacuum. [km/s]

Source: [NIST](https://physics.nist.gov/cgi-bin/cuu/Value?c|search_for=adopted_in!)
"""
const SPEED_OF_LIGHT = 299792458e-3

"""
    BOLTZMANN_CONSTANT

Boltzmann constant. [J/K]
"""
const BOLTZMANN_CONSTANT = 1.380649e-23


text:
function Rect(geometry::AbstractArray{<: Point{N, T}}) where {N,T}
    return Rect{N,T}(geometry)
end

"""
Construct a HyperRectangle enclosing all points.
"""
function Rect{N1, T1}(geometry::AbstractArray{PT}) where {N1, T1, PT <: AbstractPoint}
    N2, T2 = length(PT), eltype(PT)
    @assert N1 >= N2
    vmin = Point{N2, T2}(typemax(T2))
    vmax = Point{N2, T2}(typemin(T2))
    for p in geometry
        vmin, vmax = minmax(p, vmin, vmax)
    end
    o = vmin
    w = vmax - vmin
    if N1 > N2
        z = zero(Vec{N1 - N2, T1})
        return Rect{N1, T1}(vcat(o, z), vcat(w, z))
    else
        return Rect{N1, T1}(o, w)
   end
end

function Rect(primitive::GeometryPrimitive{N, T}) where {N, T}
    Rect{N, T}(primitive)
end

function Rect{T}(primitive::GeometryPrimitive{N, T}) where {N, T}
    Rect{N, T}(primitive)
end

function Rect{T}(a::Pyramid) where T
    w,h = a.width / T(2), a.length
    m = Vec{3,T}(a.middle)
    Rect{T}(m .- Vec{3,T}(w,w,0), m .+ Vec{3,T}(w, w, h))
end

function Rect{T}(a::Sphere) where T
    mini, maxi = extrema(a)
    Rect{T}(mini, maxi .- mini)
end

Rect{T}(a) where T = Rect{T}(coordinates(a))
Rect{N, T}(a) where {N, T} = Rect{N, T}(coordinates(a))


text:
struct Newton
end
struct NewtonCache{Tx} <: AbstractSolverCache
    x::Tx
    xold::Tx
    p::Tx
    g::Tx
end
function NewtonCache(df)
    x = copy(df.x_f)
    xold = copy(x)
    p = copy(x)
    g = copy(x)
    NewtonCache(x, xold, p, g)
end

function newtontrace(stepnorm, tracing, extended_trace, cache, df, it, tr, store_trace, show_trace)
    if tracing
        dt = Dict()
        if extended_trace
            dt["x"] = copy(cache.x)
            dt["f(x)"] = copy(value(df))
            dt["g(x)"] = copy(jacobian(df))
        end
        update!(tr,
                it,
                maximum(abs, value(df)),
                stepnorm,
                dt,
                store_trace,
                show_trace)
    end
end

function newton_(df::OnceDifferentiable,
                 initial_x::AbstractArray{T},
                 xtol::Real,
                 ftol::Real,
                 iterations::Integer,
                 store_trace::Bool,
                 show_trace::Bool,
                 extended_trace::Bool,
                 linesearch,
                 linsolve,
                 cache = NewtonCache(df)) where T
    n = length(initial_x)
    copyto!(cache.x, initial_x)
    value_jacobian!!(df, cache.x)
    check_isfinite(value(df))
    vecvalue = vec(value(df))
    it = 0
    x_converged, f_converged = assess_convergence(initial_x, cache.xold, value(df), NaN, ftol)
    stopped = any(isnan, cache.x) || any(isnan, value(df)) ? true : false

    converged = x_converged || f_converged
    x_ls = copy(cache.x)
    tr = SolverTrace()
    tracing = store_trace || show_trace || extended_trace
    newtontrace(convert(real(T), NaN), tracing, extended_trace, cache, df, it, tr, store_trace, show_trace)

    # Create objective function for the linesearch.
    # This function is defined as fo(x) = 0.5 * f(x) ⋅ f(x) and thus
    # has the gradient ∇fo(x) = ∇f(x) ⋅ f(x)
    function fo(xlin)
        value!(df, xlin)
        dot(value(df), value(df)) / 2
    end

    # The line search algorithm will want to first compute ∇fo(xₖ).
    # We have already computed ∇f(xₖ) and it is possible that it
    # is expensive to recompute.
    # We solve this using the already computed ∇f(xₖ)
    # in case of the line search asking us for the gradient at xₖ.
    function go!(storage, xlin)
        value_jacobian!(df, xlin)
        mul!(vec(storage), transpose(jacobian(df)), vecvalue)
    end
    function fgo!(storage, xlin)
        value_jacobian!(df, xlin)
        mul!(vec(storage), transpose(jacobian(df)), vecvalue)
        dot(value(df), value(df)) / 2
    end
    dfo = OnceDifferentiable(fo, go!, fgo!, cache.x, zero(real(T)))

    while !stopped && !converged && it < iterations

        it += 1

        if it > 1
            value_jacobian!(df, cache.x)
        end

        try
            mul!(vec(cache.g), transpose(jacobian(df)), vec(value(df)))
            linsolve(cache.p, jacobian(df), vec(value(df)))
            rmul!(cache.p, -1)
        catch e
            if isa(e, LAPACKException) || isa(e, SingularException)
                # Modify the search direction if the jacobian is singular
                # FIXME: better selection for lambda, see Nocedal & Wright p. 289
                fjac2 = jacobian(df)'*jacobian(df)
                lambda = convert(real(T),1e6)*sqrt(n*eps())*norm(fjac2, 1)
                linsolve(cache.p, -(fjac2 + lambda * I), vec(value(df)))
            else
                throw(e)
            end
        end

        copyto!(cache.xold, cache.x)

        if linesearch isa Static
            x_ls .= cache.x .+ cache.p
            value_jacobian!(df, x_ls)
            alpha, ϕalpha = one(real(T)), value(dfo)
        else
            value_gradient!(dfo, cache.x)
            alpha, ϕalpha = linesearch(dfo, cache.x, cache.p, one(real(T)), x_ls, value(dfo), dot(cache.g, cache.p))
        end
        # fvec is here also updated in the linesearch so no need to call f again.
        copyto!(cache.x, x_ls)
        x_converged, f_converged = assess_convergence(cache.x, cache.xold, value(df), xtol, ftol)
        stopped = any(isnan, cache.x) || any(isnan, value(df)) ? true : false

        converged = x_converged || f_converged
        newtontrace(sqeuclidean(cache.x, cache.xold), tracing, extended_trace, cache, df, it, tr, store_trace, show_trace)
    end
    return SolverResults("Newton with line-search",
                         initial_x, copy(cache.x), norm(value(df), Inf),
                         it, x_converged, xtol, f_converged, ftol, tr,
                         first(df.f_calls), first(df.df_calls))
end

function newton(df::OnceDifferentiable,
                initial_x::AbstractArray{T},
                xtol::Real,
                ftol::Real,
                iterations::Integer,
                store_trace::Bool,
                show_trace::Bool,
                extended_trace::Bool,
                linesearch,
                cache = NewtonCache(df);
                linsolve=(x, A, b) -> copyto!(x, A\b)) where T
    newton_(df, initial_x, convert(real(T), xtol), convert(real(T), ftol), iterations, store_trace, show_trace, extended_trace, linesearch, linsolve, cache)
end


text:
# type
abstract type Shape{T} end

struct Ellipse{T} <: Shape{T}
    a::T
    b::T
end

# parametrics
struct Rectangle{T1,T2} <: Shape{T1}
    len::T1
    wid::T2
end

struct Square{T} <: Shape{T}
    side::T
end

function area(s::Shape)

    if isa(s, Rectangle)
        return s.len * s.wid
    else
        return s.side^2
    end

end

function main()

    r = Rectangle{Float64, Float64}(5.0, 4.0)
    s = Square{Float64}(5.0)

    println(area(r))

end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
# TODO: Is there a way to generate docs without needing to
# generate them on a gpu node
using RAPIDS
using Documenter

Documenter.makedocs(
    modules = [RAPIDS],
    sitename = "RAPIDS.jl",
)


text:
using JLD, JLD2, Test, Distributions, Random, StateSpaceRoutines, BenchmarkTools, ParallelDataTransfer, FLoops

# Read in from JLD
tpf_main_input = load("reference/tpf_main_inputs.jld2")
data = tpf_main_input["data"]
TTT = tpf_main_input["TTT"]
RRR = tpf_main_input["RRR"]
CCC = tpf_main_input["CCC"]
ZZ = tpf_main_input["ZZ"]
DD = tpf_main_input["DD"]
F_ϵ = tpf_main_input["F_epsilon"]
F_u = tpf_main_input["F_u"]
s_init = tpf_main_input["s_init"]

# s_init = repeat(s_init, outer=(1,1000))

# Tune algorithm
tuning = Dict(:r_star => 2., :c_init => 0.3, :target_accept_rate => 0.4,
              :resampling_method => :systematic, :n_mh_steps => 1,
              :n_particles => 1000,#1000000,
              :n_presample_periods => 0,
              :allout => true, :parallel => true)

# Define Φ and Ψ (can't be saved to JLD)
Φ(s_t::AbstractVector{Float64}, ϵ_t::AbstractVector{Float64}) = TTT*s_t + RRR*ϵ_t + CCC
Ψ(s_t::AbstractVector{Float64}) = ZZ*s_t + DD

# Load in test inputs and outputs
test_file_inputs = load("reference/tpf_aux_inputs.jld2")
test_file_outputs = load("reference/tpf_aux_outputs.jld2")

φ_old = test_file_inputs["phi_old"]
norm_weights = test_file_inputs["norm_weights"]
coeff_terms = test_file_inputs["coeff_terms"]
log_e_1_terms = test_file_inputs["log_e_1_terms"]
log_e_2_terms = test_file_inputs["log_e_2_terms"]
inc_weights = test_file_inputs["inc_weights"]
HH = cov(F_u)
s_t_nontemp = test_file_inputs["s_t_nontemp"]

#=s_t_nontemp = repeat(s_t_nontemp, outer=(1,1000))
coeff_terms = repeat(coeff_terms,1000)
log_e_1_terms = repeat(log_e_1_terms, 1000)
log_e_2_terms = repeat(log_e_2_terms, 1000)
inc_weights = repeat(inc_weights, 1000)
norm_weights = repeat(norm_weights, 1000)=#

ENV["frbnyjuliamemory"] = "5G"
myprocs = addprocs_frbny(48)
@everywhere using JLD, JLD2, Test, Distributions, Random, StateSpaceRoutines, BenchmarkTools, ParallelDataTransfer

#=
x = randn(10,10)
sendto(workers(), x=x)
@show fetch(@spawnat 5 x)
x[1:10,1:10] .= 1.0
@show fetch(@spawnat 5 x)
=#

@everywhere tpf_main_input = load("reference/tpf_main_inputs.jld2")
@everywhere TTT = tpf_main_input["TTT"]
@everywhere RRR = tpf_main_input["RRR"]
@everywhere CCC = tpf_main_input["CCC"]
@everywhere ZZ = tpf_main_input["ZZ"]
@everywhere DD = tpf_main_input["DD"]
@everywhere Φ(s_t::AbstractVector{Float64}, ϵ_t::AbstractVector{Float64}) = TTT*s_t + RRR*ϵ_t + CCC
@everywhere Ψ(s_t::AbstractVector{Float64}) = ZZ*s_t + DD

StateSpaceRoutines.sendto(workers(), Φ = Φ)
StateSpaceRoutines.sendto(workers(), Ψ = Ψ)#=
StateSpaceRoutines.sendto(workers(), coeff_terms = coeff_terms)
StateSpaceRoutines.sendto(workers(), log_e_1_terms = log_e_1_terms)
StateSpaceRoutines.sendto(workers(), log_e_2_terms = log_e_2_terms)
StateSpaceRoutines.sendto(workers(), data = data)
StateSpaceRoutines.sendto(workers(), s_t_nontemp = s_t_nontemp)
StateSpaceRoutines.sendto(workers(), HH = HH)=#

sendto(workers(), HH=HH)
sendto(workers(), s_t_nontemp=s_t_nontemp)
@everywhere inv_HH = inv(HH)

#=weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ, data[:, 47], s_t_nontemp, det(HH), inv_HH;
               initialize = false, parallel = false)=#
@btime weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ, data[:, 47], s_t_nontemp, det(HH), inv_HH;
               initialize = false, parallel = true)
weight_kernel!(coeff_terms, log_e_1_terms, log_e_2_terms, φ_old, Ψ, data[:, 47], s_t_nontemp, det(HH), inv_HH;
               initialize = false, parallel = true)
φ_new = next_φ(φ_old, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]), tuning[:r_star], 2)
correction!(inc_weights, norm_weights, φ_new, coeff_terms, log_e_1_terms, log_e_2_terms, length(data[:,47]))

@testset "Corection and Auxiliary Tests" begin
    @test coeff_terms[1] ≈ test_file_outputs["coeff_terms"][1]
    @test log_e_1_terms[1] ≈ test_file_outputs["log_e_1_terms"][1]
    @test log_e_2_terms[1] ≈ test_file_outputs["log_e_2_terms"][1]
    @test φ_new ≈ test_file_outputs["phi_new"]
    @test inc_weights[1] ≈ test_file_outputs["inc_weights"][1]
end

## Selection Tests
s_t1_temp = test_file_inputs["s_t1_temp"]
ϵ_t = test_file_inputs["eps_t"]

#=s_t1_temp = repeat(s_t1_temp, outer=(1,1000))
ϵ_t = repeat(ϵ_t, outer=(1,1000))=#

#=
ENV["frbnyjuliamemory"] = "5G"
myprocs = addprocs_frbny(48)
@everywhere using JLD2, Test, Distributions, Random, StateSpaceRoutines, BenchmarkTools, ParallelDataTransfer
#@everywhere include("../src/filters/tempered_particle_filter/correction.jl")
#@everywhere include("../src/util.jl")
=#
Random.seed!(47)
selection!(norm_weights, s_t1_temp, s_t_nontemp,ϵ_t, resampling_method = tuning[:resampling_method])
@testset "Selection Tests" begin
    @test s_t1_temp[1] ≈ test_file_outputs["s_t1_temp"][1]
    @test s_t_nontemp[1] ≈ test_file_outputs["s_t_nontemp"][1]
    @test ϵ_t[1] ≈ test_file_outputs["eps_t"][1]
end

## Mutation Tests
QQ = cov(F_ϵ)
accept_rate = test_file_inputs["accept_rate"]
c = test_file_inputs["c"]

c = update_c(c, accept_rate, tuning[:target_accept_rate])
Random.seed!(47)

sendto(workers(), s_t1_temp = s_t1_temp)
sendto(workers(), ϵ_t = ϵ_t)
sendto(workers(), φ_new = φ_new)

StateSpaceRoutines.mutation!(Φ, Ψ, QQ, det(HH), inv_HH, φ_new, data[:,47], s_t_nontemp, s_t1_temp, ϵ_t, c, tuning[:n_mh_steps],
                             parallel = true)

@testset "Mutation Tests" begin
    @test s_t_nontemp[1] ≈ test_file_outputs["s_t_nontemp_mutation"][1]
    @test ϵ_t[1] ≈ test_file_outputs["eps_t_mutation"][1]
end

sendto(workers(), y_t = data[:,47])

@btime StateSpaceRoutines.mutation!(Φ, Ψ, QQ, det(HH), inv_HH, φ_new, data[:,47], s_t_nontemp, s_t1_temp, ϵ_t, c,
                                    tuning[:n_mh_steps], parallel = true)
sendto(workers(), data=data)
sendto(workers(), s_init = s_init)
sendto(workers(), F_ϵ = F_ϵ)
sendto(workers(), F_u = F_u)
@btime out_parallel_one_worker = tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; tuning..., verbose = :none, parallel = true)

## Whole TPF Tests
Random.seed!(47)
out_no_parallel = tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; tuning..., verbose = :none, parallel = false)
Random.seed!(47)
out_parallel_one_worker = tempered_particle_filter(data, Φ, Ψ, F_ϵ, F_u, s_init; tuning..., verbose = :none, parallel = true)
@testset "TPF tests" begin
    @test out_no_parallel[1] ≈ -302.99967306704133
    @test out_parallel_one_worker[1] ≈ -302.99967306704133
    ## Note when using more than 1 worker, this is not true because there is still a random step in mh_steps
#=
    # In Julia 1.5, seeding appears to be different
    if VERSION >=  v"1.5"
        @test out_parallel_one_worker[1] ≈ -307.9285106003482
    elseif VERSION >= v"1.0" # This is different than in Julia6 because @distributed seeds differently than @parallel
        @test out_parallel_one_worker[1] ≈ -303.64727963725073 # Julia 6 was: -306.8211172094595
    end
=#
end

for i in workers()
    rmprocs(i)
end