text:
if get(ARGS, 1, "Array") == "CuArray"
    using CUDA
    const ArrayType = CUDA.CuArray
    CUDA.allowscalar(false)
else
    const ArrayType = Array
end

using Test
using NonlinearSolvers

@testset "NonlinearSolvers tests" begin
    function f!(F, x)
        t = Tuple(x)
        F_nt = ntuple(Val(length(F))) do i
            i == 1 && (F_i = (t[1] + 3) * (t[2]^3 - 7) + 18)
            i == 2 && (F_i = sin(t[2] * exp(t[1]) - 1))
            F_i
        end
        F .= F_nt
    end

    function j!(J, x)
        t = Tuple(x)
        u = exp(t[1]) * cos(t[2] * exp(t[1]) - 1)
        J_nt = ntuple(Val(length(J))) do i
            i == 1 && (J_ij = t[2]^3 - 7)
            i == 2 && (J_ij = t[2] * u)
            i == 3 && (J_ij = 3 * t[2]^2 * (t[1] + 3))
            i == 4 && (J_ij = u)
            return J_ij
        end
        J[:] .= J_nt
    end
    FT = Float64
    N = 2
    x_initial = FT[0.1; 1.2]
    correct_sol = FT[0, 1]

    @show ArrayType
    x_initial = ArrayType(x_initial)

    @testset "Compact solution converging" begin
        # Jacobian provided
        nls = NewtonsMethod(f!, j!, x_initial)
        sol = solve!(nls)
        @test all(isapprox.(correct_sol, Array(sol.root); atol = 10eps(FT)))
        @test sol.converged == true

        # Jacobian-free
        nls = NewtonsMethodAD(f!, x_initial)
        sol = solve!(nls)
        @test all(isapprox.(correct_sol, Array(sol.root); atol = 10eps(FT)))
        @test sol.converged == true
    end

    @testset "Compact solution non-converging" begin
        # Jacobian provided
        nls = NewtonsMethod(f!, j!, x_initial)
        sol = solve!(nls, CompactSolution(), SolutionTolerance(eps(FT)), 1)
        @test sol.converged == false

        # Jacobian-free
        nls = NewtonsMethodAD(f!, x_initial)
        sol = solve!(nls, CompactSolution(), SolutionTolerance(eps(FT)), 1)
        @test sol.converged == false
    end

    if x_initial isa Array
        @testset "Verbose solution" begin
            # Jacobian provided
            nls = NewtonsMethod(f!, j!, x_initial)
            sol = solve!(nls, VerboseSolution())
            @test all(isapprox.(correct_sol, Array(sol.root); atol = 10eps(FT)))
            @test sol.converged
            @test sol.err[end] < 10eps(FT)
            @test hasproperty(sol, :iter_performed)
            @test sol.root_history isa Vector{typeof(nls.x1)}
            @test sol.err_history isa Vector{typeof(nls.x1)}

            # Jacobian-free
            nls = NewtonsMethodAD(f!, x_initial)
            sol = solve!(nls, VerboseSolution())
            @test all(isapprox.(correct_sol, Array(sol.root); atol = 10eps(FT)))

            @test sol.converged
            @test sol.err[end] < 10eps(FT)
            @test hasproperty(sol, :iter_performed)
            @test sol.root_history isa Vector{typeof(nls.x1)}
            @test sol.err_history isa Vector{typeof(nls.x1)}
        end
    end
end


text:
using Test

function naive!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:col
            for k in 1:colA
                C[i, j] = C[i, j] + A[i, k] * B[k, j]
            end
        end
    end
    C
end

"""
inner-product-form
"""
function ipform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for i in 1:row
        for j in 1:col
            dataC = 0.0
            for k in 1:colA
                dataC = dataC + A[i, k] * B[k, j]
            end
            C[i, j] = dataC
        end
    end
end

"""
outer-product-form
"""
function opform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for k in 1:colA
        for j in 1:col
            dataB = B[k, j]
            for i in 1:row
                C[i, j] = C[i, j] + A[i, k] * dataB
            end
        end
    end
end

"""
middle-product-form
"""
function mpform!(C, A, B)
    row, col  = size(C)
    colA = size(A)[2]
    for j in 1:col
        for k in 1:colA
            dataB = B[k, j]
            for i in 1:row
                C[i, j] = C[i, j] + A[i, k] * dataB
            end
        end
    end
end

function main(N)
    A = reshape(collect(1:1.0:N^2), N, N)
    B = reshape(collect(1:1.0:N^2), N, N)
    C = zeros(N, N)
    @time naive!(C, A, B)
    # copy result
    refC = copy(C)

    C = zeros(N, N)
    @time ipform!(C, A, B)
    @test C ≈ refC

    C = zeros(N, N)
    @time opform!(C, A, B)
    @test C ≈ refC

    C = zeros(N, N)
    @time mpform!(C, A, B)
    @test C ≈ refC
end

ext: jl

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
const curandGenerator_t = Ptr{Cvoid}

mutable struct RNG <: Random.AbstractRNG
    ptr::curandGenerator_t
    typ::Int
end

Base.unsafe_convert(::Type{curandGenerator_t}, rng::RNG) = rng.ptr

const curandDiscreteDistribution_t = Ptr{Cvoid}

mutable struct DiscreteDistribution
    ptr::curandDiscreteDistribution_t
end

Base.unsafe_convert(::Type{curandDiscreteDistribution_t}, dist::DiscreteDistribution) = dist.ptr

# CURAND status codes
const curandStatus_t = UInt32
const CURAND_STATUS_SUCCESS = 0
const CURAND_STATUS_VERSION_MISMATCH = 100
const CURAND_STATUS_NOT_INITIALIZED = 101
const CURAND_STATUS_ALLOCATION_FAILED = 102
const CURAND_STATUS_TYPE_ERROR = 103
const CURAND_STATUS_OUT_OF_RANGE = 104
const CURAND_STATUS_LENGTH_NOT_MULTIPLE = 105
const CURAND_STATUS_DOUBLE_PRECISION_REQUIRED = 106
const CURAND_STATUS_LAUNCH_FAILURE = 201
const CURAND_STATUS_PREEXISTING_FAILURE = 202
const CURAND_STATUS_INITIALIZATION_FAILED = 203
const CURAND_STATUS_ARCH_MISMATCH = 204
const CURAND_STATUS_INTERNAL_ERROR = 999

# CURAND RNG types (curandRngType)
const CURAND_RNG_TEST = 0
const CURAND_RNG_PSEUDO_DEFAULT = 100
const CURAND_RNG_PSEUDO_XORWOW = 101
const CURAND_RNG_PSEUDO_MRG32K3A = 121
const CURAND_RNG_PSEUDO_MTGP32 = 141
const CURAND_RNG_PSEUDO_MT19937 = 142
const CURAND_RNG_PSEUDO_PHILOX4_32_10 = 161
const CURAND_RNG_QUASI_DEFAULT = 200
const CURAND_RNG_QUASI_SOBOL32 = 201
const CURAND_RNG_QUASI_SCRAMBLED_SOBOL32 = 202
const CURAND_RNG_QUASI_SOBOL64 = 203
const CURAND_RNG_QUASI_SCRAMBLED_SOBOL64 = 204

# CURAND ordering of results in memory
const CURAND_ORDERING_PSEUDO_BEST = 100
const CURAND_ORDERING_PSEUDO_DEFAULT = 101
const CURAND_ORDERING_PSEUDO_SEEDED = 102
const CURAND_ORDERING_QUASI_DEFAULT = 201

# CURAND choice of direction vector set
const CURAND_DIRECTION_VECTORS_32_JOEKUO6 = 101
const CURAND_SCRAMBLED_DIRECTION_VECTORS_32_JOEKUO6 = 102
const CURAND_DIRECTION_VECTORS_64_JOEKUO6 = 103
const CURAND_SCRAMBLED_DIRECTION_VECTORS_64_JOEKUO6 = 104

# CURAND method
const CURAND_CHOOSE_BEST = 0
const CURAND_ITR = 1
const CURAND_KNUTH = 2
const CURAND_HITR = 3
const CURAND_M1 = 4
const CURAND_M2 = 5
const CURAND_BINARY_SEARCH = 6
const CURAND_DISCRETE_GAUSS = 7
const CURAND_REJECTION = 8
const CURAND_DEVICE_API = 9
const CURAND_FAST_REJECTION = 10
const CURAND_3RD = 11
const CURAND_DEFINITION = 12
const CURAND_POISSON = 13


text:
#
# This script was used to generate results for the camera-ready AABI submission.
#

println("pseudo_point_synthetic_experiments_irregular_obs_sum_separable")

using Revise

using Pkg
Pkg.activate(".");
Pkg.resolve();
Pkg.instantiate();

using ArgParse
using CSV

# Make it possible to parse strings from the command line into Vector{Int}s.
function ArgParse.parse_item(::Type{Vector{Int}}, x::String)
    @assert x[1] == '[' && x[end] == ']'
    return CSV.read(IOBuffer(x[2:prevind(x, end)]); header=false, transpose=true).Column1
end

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--run-experiments"
            help = "true / false - run the experiments or not"
            arg_type = Bool
            default = false
        "--N-space-per-time"
            help = "number of points to consider in space per time point"
            arg_type = Int
            default = 10
        "--N-time"
            help = "number of time points to use -- expects multiple settings"
            arg_type = Vector{Int}
            default = [
                10, 25, 50,
                100, 250, 500,
                1_000, 2_500, 5_000,
                10_000, 25_000, 50_000,
                100_000,
            ]
        "--M-per-time"
            help = "number of pseudo-points per time point -- expects multiple settings. 0 means exact"
            arg_type = Vector{Int}
            default = [0, 5, 10, 20]
        "--BLAS-threads"
            help = "number of threads to let BLAS use"
            arg_type = Int
            default = 1
    end
    return parse_args(s)
end

args = parse_commandline()

using BenchmarkTools
using DrWatson
using LinearAlgebra
using Random
using Stheno
using TemporalGPs

BLAS.set_num_threads(args["BLAS-threads"])

const exp_dir_name = "pseudo_point_scaling_irregular_obs_sum_separable"

#
# Specify some functionality required to run the experiments.
#

function separable_kernel(σ::Real, λ_space::Real, λ_time::Real)
    k_space = σ^2 * stretch(EQ(), λ_space)
    k_time = stretch(Matern52(), λ_time)
    return TemporalGPs.Separable(k_space, k_time)
end

function sum_separable_kernel(σ::Real, λ_space::Real, λ_time::Real)
    k1 = separable_kernel(σ, λ_space, λ_time)
    k2 = separable_kernel(1.1 * σ, 0.3 * λ_space, 1.05 * λ_time)
    return k1 + k2
end

#
# Specify experiments to run.
#

# Ensure that output directory for results exists.
const data_dir = joinpath(datadir(), exp_dir_name)
mkpath(data_dir)

# Ensure that directory for experiment specification exists.
mkpath(joinpath(data_dir, "meta"))

wsave(
    joinpath(data_dir, "meta", "settings.bson"),
    :settings => dict_list(
        Dict(
            # The number of kernels to include in the summation.
            :model => [
                (
                    name = "sum-separable",
                    θ = (
                        σ = 0.92,
                        λ_space = 1.1,
                        λ_time = 0.83,
                        S_min = 0.1,
                    ),
                    ctor = sum_separable_kernel,
                ),
            ],

            # The total number of points in space to consider.
            :N_space_per_time => args["N-space-per-time"],

            # The total number of points in time to consider.
            :Ns_time => args["N-time"],

            # The total number of pseudo-points per time point.
            :Ms_per_time => args["M-per-time"],
        ),
    ),
)

#
# Run experiments.
#

if args["run-experiments"]

println("Running experiments")

let
    settings = load(joinpath(data_dir, "meta", "settings.bson"))[:settings]

    for (n, setting) in enumerate(settings)
        println("Starting experiment $n / $(length(settings))")

        # Set seed for reproducibility.
        rng = MersenneTwister(123456)

        # Extract experiment parameters.
        model_name = setting[:model].name
        @show model_name
        θ = setting[:model].θ
        ctor = setting[:model].ctor
        T = setting[:Ns_time]
        N_space_per_time = setting[:N_space_per_time]
        M_per_time = setting[:Ms_per_time]

        # Specify input locations. We place these irregularly through time.
        t = sort(rand(rng, T)) .* 50
        r = [rand(rng, N_space_per_time) * 10 for _ in 1:T]
        x = TemporalGPs.RegularInTime(t, r)
        x_naive = collect(x)

        # Specify inducing inputs. They are regularly-spaced between 0 and 10.
        z_r = collect(range(0.0, 10.0; length=M_per_time))
        z_naive = collect(TemporalGPs.RectilinearGrid(z_r, t))

        # Specify observation noise.
        S = [Diagonal(rand(rng, N_space_per_time) .+ θ.S_min) for _ in 1:T]
        S_naive = vcat(diag.(S)...)

        # Generate data exactly / naively.
        k = ctor(θ.σ, θ.λ_space, θ.λ_time)
        f = GP(k, GPC())

        # Hack together generation.
        N_per_slice = 5_000
        ys_naive = map(1:div(length(x_naive), N_per_slice)) do j
            idx = N_per_slice * (j - 1) + 1:min(N_per_slice * j, length(x_naive))
            return rand(rng, f(x_naive[idx], S_naive[idx]))
        end
        ys_naive = if rem(length(x_naive), N_per_slice) == 0
            ys_naive
        else
            idx = length(x_naive) - rem(length(x_naive), N_per_slice) + 1:length(x_naive)
            vcat(
                ys_naive,
                rand(rng, f(x_naive[idx], S_naive[idx]))
            )
        end
        y_naive = vcat(ys_naive...)
        y = TemporalGPs.match_to(y_naive, x)

        display((length(x_naive), length(y_naive)))
        println()
        if M_per_time == 0 && length(y_naive) <= 10_000
            fx_naive = f(x_naive, S_naive)
            println("Performing exact inference.")
            wsave(
                joinpath(
                    data_dir,
                    savename(
                        Dict(
                            :k => model_name,
                            :N_space_per_time => N_space_per_time,
                            :N_time => T,
                            :M_per_time => M_per_time,
                            :method => :exact,
                        ),
                        "bson",
                    ),
                ),
                Dict(
                    :setting => setting,
                    :lml_result => logpdf(fx_naive, y_naive),
                    :timing => (@benchmark logpdf($fx_naive, $y_naive)),
                    :method => :exact,
                ),
            )
        end

        # State-Space ELBO
        println("Using the state-space ELBO")
        if M_per_time > 0
            wsave(
                joinpath(
                    data_dir,
                    savename(
                        Dict(
                            :k => model_name,
                            :N_space_per_time => N_space_per_time,
                            :N_time => T,
                            :M_per_time => M_per_time,
                            :method => :state_space_elbo,
                        ),
                        "bson",
                    ),
                ),
                Dict(
                    :setting => setting,
                    :lml_result => elbo(to_sde(f)(x, S_naive), y_naive, z_r),
                    :timing => (@benchmark elbo(to_sde($f)($x, $S_naive), $y_naive, $z_r)),
                    :method => :state_space_elbo,
                ),
            )
        end
    end
end

end # run-experiments

#
# Summarise results. Generate output (plots and tables).
#

using ColorTypes, DataFrames, LaTeXStrings, PGFPlotsX, Plots, TexTables
pgfplotsx();

Plots.default(
    tickfontsize=nothing,
    legendfontsize=nothing,
    guidefontsize=nothing,
);

# Use default document font-size if nothing provided.
function Plots.pgfx_font(fontsize::Nothing, thickness_scaling = 1, font = "\\selectfont")
    return string("{", font, "}")
end

# Ensure that results directory exists.
const results_dir = joinpath(projectdir(), "results", exp_dir_name)
mkpath(results_dir)

let
    # Load the results and augment columns with extra info from `settings` column.
    results_unsorted = DataFrame(
        map(eachrow(collect_results(data_dir))) do row
            merge(
                copy(row),
                (
                    N_time = row.setting[:Ns_time],
                    M_per_time = row.setting[:Ms_per_time],
                ),
            )
        end
    )

    results = sort(results_unsorted, [:method, :M_per_time]; rev=true)

    # Each inference method gets its own marker shape.
    marker_shapes = Dict(
        :exact => :square,
        :naive_elbo => :circle,
        :state_space_elbo => :circle,
    )

    # Each inference method gets its own colour.
    colours = Dict(
        :exact => :black,
        :naive_elbo => :red,
        :state_space_elbo => :blue,
    )

    method_name_map = Dict(
        :exact => "exact",
        :naive_elbo => "elbo (naive)",
        :state_space_elbo => "elbo"
    )

    # Each M gets its own line style. This means that the maximum value of `M` is 5.
    # We really shouldn't be using that many anyway though...
    line_styles = [:solid, :solid, :dash, :dot, :dashdot, :dashdotdot]
    Ms_per_time = sort(unique(results.M_per_time))
    line_styles_map = Dict(
        M => style for (M, style) in zip(Ms_per_time, line_styles[eachindex(Ms_per_time)])
    )

    # Locations to save the results to.
    save_dirs = [
        results_dir,
        "/Users/willtebbutt/Dropbox/University/PhD/first_year_report/chapter_3/pseudo_points"
    ]

    # Ensure that save_dirs exist.
    foreach(mkpath, save_dirs)

    let

        # width_pts = 246.0  # or any other value
        # inches_per_points = 1.0 / 72.27
        # width_inches = width_pts *inches_per_points
        width_inches = 2.8
        width_px = width_inches * Plots.DPI  # or  width_inches*DPI
        height_px = width_px

        plt_lml = plot(
            xlabel="T",
            ylabel="ELBO / T",
            legend=:bottomright,
            legendfontsize=9,
            background_color_legend=RGBA(1, 1, 1, 0.6),
            size=(width_px, height_px),
            grid = (1, 0.4),
        );
        plt_comp_time = plot(
            xlabel="T",
            ylabel="Compute Time (s)",
            legend=:bottomright,
            legendfontsize=9,
            background_color_legend=RGBA(1, 1, 1, 0.6),
            size=(width_px, height_px),
            grid = (1, 0.4),
        )
        for group in groupby(results, [:method, :M_per_time])
            method_name = only(unique(group.method))
            M_per_time = only(unique(group.M_per_time))
            if method_name == :exact && M_per_time != first(Ms_per_time)
                continue
            end
            if method_name == :naive_elbo
                continue
            end

            idx = sortperm(group.N_time)
            label = if method_name == :exact
                method_name_map[method_name]
            else
                L"M_\tau=%$M_per_time"
            end

            Ts = group.N_time[idx]

            plot!(plt_lml, Ts, group.lml_result[idx] ./Ts ;
                markershape=marker_shapes[method_name],
                color=colours[method_name],
                linestyle=line_styles_map[M_per_time],
                markerstrokewidth=0,
                markersize=3,
                label=label,
                xscale=:log10,
            )

            min_time = getfield.(minimum.(group.timing), :time)[idx] ./ 1e9
            plot!(plt_comp_time, Ts, min_time;
                markershape=marker_shapes[method_name],
                color=colours[method_name],
                linestyle=line_styles_map[M_per_time],
                markersize=3,
                markerstrokewidth=0,
                label=label,
                yscale=:log10,
                xscale=:log10,
            )
        end
        foreach(
            save_dir -> begin
                savefig(
                    plt_lml,
                    joinpath(save_dir, "irregular_lml_plot_sum_separable.tikz"),
                )
                savefig(
                    plt_comp_time,
                    joinpath(save_dir, "irregular_timing_plot_sum_separable.tikz"),
                )
                savefig(
                    plot(plt_lml, plt_comp_time; layout=(2, 1)),
                    joinpath(save_dir, "irregular_lml_timing_plot_sum_separable.tikz"),
                )
            end,
            save_dirs,
        )
    end
end


text:
const PInv{Style, Typ} = Applied{Style, typeof(pinv), <:Tuple{Typ}}
const Inv{Style, Typ} = Applied{Style, typeof(inv), <:Tuple{Typ}}

Inv(A) = applied(inv, A)
PInv(A) = applied(pinv, A)

const InvOrPInv = Union{PInv, Inv}

parent(A::InvOrPInv) = first(A.args)

pinv(A::PInv) = parent(A)
function inv(A::PInv)
    checksquare(parent(A))
    parent(A)
end

inv(A::Inv) = parent(A)
pinv(A::Inv) = inv(A)

ndims(A::InvOrPInv) = ndims(parent(A))

size(A::InvOrPInv) = reverse(size(parent(A)))
axes(A::InvOrPInv) = reverse(axes(parent(A)))
size(A::InvOrPInv, k) = size(A)[k]
axes(A::InvOrPInv, k) = axes(A)[k]
eltype(A::InvOrPInv) = Base.promote_op(inv, eltype(parent(A)))

struct LdivApplyStyle <: ApplyStyle end

ldivapplystyle(_, _) = LdivApplyStyle()
ldivapplystyle(::LazyLayout, ::LazyLayout) = LazyArrayApplyStyle()
ldivapplystyle(::LazyLayout, _) = LazyArrayApplyStyle()
ldivapplystyle(_, ::LazyLayout) = LazyArrayApplyStyle()
ApplyStyle(::typeof(\), ::Type{A}, ::Type{B}) where {A<:AbstractArray,B<:AbstractArray} =
    ldivapplystyle(MemoryLayout(A), MemoryLayout(B))

axes(M::Applied{Style,typeof(\)}) where Style = ldivaxes(M.args...)
axes(M::Applied{Style,typeof(\)}, p::Int)  where Style = axes(M)[p]
size(M::Applied{Style,typeof(\)}) where Style = length.(axes(M))

@inline eltype(M::Applied{Style,typeof(\)}) where Style = eltype(Ldiv(M.args...))
@inline ndims(M::Applied{Style,typeof(\)}) where Style = ndims(last(M.args))

check_applied_axes(A::Applied{<:Any,typeof(\)}) = check_ldiv_axes(A.args...)

######
# PInv/Inv
########

const PInvMatrix{T,Arg} = ApplyMatrix{T,typeof(pinv),<:Tuple{Arg}}
const InvMatrix{T,Arg} = ApplyMatrix{T,typeof(inv),<:Tuple{Arg}}

PInvMatrix(A) = ApplyMatrix(pinv, A)
function InvMatrix(A)
    checksquare(A)
    ApplyMatrix(inv, A)
end

parent(A::PInvMatrix) = first(A.args)
parent(A::InvMatrix) = first(A.args)
axes(A::PInvMatrix) = reverse(axes(parent(A)))
size(A::PInvMatrix) = map(length, axes(A))

@propagate_inbounds getindex(A::PInvMatrix{T}, k::Int, j::Int) where T =
    (parent(A)\[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k]

@propagate_inbounds getindex(A::InvMatrix{T}, k::Int, j::Int) where T =
    (parent(A)\[Zeros(j-1); one(T); Zeros(size(A,2) - j)])[k]

struct InvLayout{L} <: MemoryLayout end
struct PInvLayout{L} <: MemoryLayout end

applylayout(::Type{typeof(inv)}, ::A) where A = InvLayout{A}()
applylayout(::Type{typeof(pinv)}, ::A) where A = PInvLayout{A}()

mulapplystyle(::InvLayout{A}, B) where A = ldivapplystyle(A, B)
mulapplystyle(::PInvLayout{A}, B) where A = ldivapplystyle(A, B)

@inline function Ldiv(M::Mul)
    Ai,b = M.args
    Ldiv(parent(Ai), b)
end
Ldiv(A::Applied{<:Any,typeof(\)}) = Ldiv(A.args...)

similar(M::Applied{LdivApplyStyle}, ::Type{T}) where T = similar(Ldiv(M), T)
copy(M::Applied{LdivApplyStyle}) = copy(Ldiv(M))
@inline copyto!(dest::AbstractArray, M::Applied{LdivApplyStyle}) = copyto!(dest, Ldiv(M))
@inline materialize!(M::Applied{LdivApplyStyle}) = materialize!(Ldiv(M))

@propagate_inbounds getindex(A::Applied{LazyArrayApplyStyle,typeof(\)}, kj...) =
    materialize(Ldiv(A))[kj...]

###
# * layout
###

function copy(L::Ldiv{<:Any,ApplyLayout{typeof(*)}})
    args = arguments(L.B)
    apply(*, L.A \  first(args),  tail(args)...)
end


text:
using Test

using iLQGames:
    ProductSystem,
    Unicycle4D,
    n_states,
    n_controls,
    dx,
    linearize

using StaticArrays

# create two separate dynamics systems
unicycle1 = Unicycle4D{0.1}()
unicycle2 = Unicycle4D{0.1}()
xid1 = SVector{4}(1:4)
xid2 = SVector{4}(5:8)
uid1 = SVector{2}(1:2)
uid2 = SVector{2}(3:4)

# and combine them in the product system
dyn = ProductSystem((unicycle1, unicycle2))

# test some operations on the joint system
x = @SVector rand(n_states(dyn))
u = @SVector rand(n_controls(dyn))
t = 0.0;

@testset "dx" begin
    dx_dyn = dx(dyn, x, u, t)
    dx_unicycle1 = dx(unicycle1, x[xid1], u[uid1], t)
    dx_unicycle2 = dx(unicycle2, x[xid2], u[uid2], t)

    @test isapprox(dx_dyn[xid1], dx_unicycle1)
    @test isapprox(dx_dyn[xid2], dx_unicycle2)
end

@testset "linearize" begin
    lin_dyn = linearize(dyn, x, u, t)
    lin_unicycle1 = linearize(unicycle1, x[xid1], u[uid1], t)
    lin_unicycle2 = linearize(unicycle2, x[xid2], u[uid2], t)

    # linearizing the combined system should be the same as linearizing the
    # subsystems separately
    @test isapprox(lin_dyn.A[xid1, xid1], lin_unicycle1.A)
    @test isapprox(lin_dyn.A[xid2, xid2], lin_unicycle2.A)
    @test isapprox(lin_dyn.B[xid1, uid1], lin_unicycle1.B)
    @test isapprox(lin_dyn.B[xid2, uid2], lin_unicycle2.B)
end


text:
using SplittablesBase: halve, amount
using Base.Threads: @spawn, nthreads, threadid

"""
    abstract type ThreadingStrategy

Controls how threading is performed in [`lomc!`](@ref).

# Interface

* [`fciqmc_step!`](@ref)
* [`working_memory`](@ref)

# Implemented Strategies

* [`NoThreading`](@ref)
* [`ThreadsThreading`](@ref)
* [`SplittablesThreading`](@ref)
* [`ThreadsXSumThreading`](@ref)
* [`ThreadsXForeachThreading`](@ref)

"""
abstract type ThreadingStrategy end

"""
    ntuple_working_memory(dv::AbstractDVec)

Create a `NTuple{N}` of vectors that are `similar` to `dv`, where `N = Threads.nthreads()`.
"""
ntuple_working_memory(dv) = ntuple_working_memory(localpart(dv))
function ntuple_working_memory(v::AbstractDVec)
    return Tuple(similar(v) for _ in 1:nthreads())
end
function ntuple_working_memory(v::AbstractVector)
    return Tuple(similar(v) for _ in 1:nthreads())
end

"""
    fciqmc_step!(t_strat::ThreadingStrategy, w, ham, v, shift, dτ) -> stat_names, stats

Perform a single matrix(/operator)-vector multiplication:

```math
v^{(n + 1)} = [1 - dτ(\\hat{H} - S)]⋅v^n ,
```

where `Ĥ == ham` and `S == shift`.

Whether the operation is performed in stochastic, semistochastic, or determistic way is
controlled by the trait `StochasticStyle(w)`. See [`StochasticStyle`](@ref).

Whether the multiplication is performed on multiple threads is controlled by `t_strat` (see
[`ThreadingStrategy`](@ref)). `w` is the working memory corresponding to `t_strat` (see
[`working_memory`](@ref)).

Returns the step stats generated by the `StochasticStyle`.
"""
fciqmc_step!

"""
    working_memory(t_strat::ThreadingStrategy, dv)

Create a working memory instance compatible with `t_strat`. The working memory must be
compatible with [`sort_into_targets!`](@ref).
"""
working_memory

"""
    NoThreading <: ThreadingStrategy

[`ThreadingStrategy`](@ref) that disables threading.
"""
struct NoThreading <: ThreadingStrategy end

working_memory(::NoThreading, dv) = similar(localpart(dv))

function fciqmc_step!(::NoThreading, w, ham, dv, shift, dτ)
    # single-threaded version suitable for MPI
    v = localpart(dv)
    @assert w ≢ v "`w` and `v` must not be the same object"
    zero!(w) # clear working memory

    stat_names, stats = step_stats(v, Val(1))
    for (add, val) in pairs(v)
        stats += fciqmc_col!(w, ham, add, val, shift, dτ)
    end
    return stat_names, stats
end

"""
    ThreadsThreading <: ThreadingStrategy

[`ThreadingStrategy`](@ref) based on `Base.Threads`.
"""
struct ThreadsThreading <: ThreadingStrategy end

working_memory(::ThreadsThreading, dv) = ntuple_working_memory(dv)

function fciqmc_step!(::ThreadsThreading, ws::NTuple{N}, ham, dv, shift, dτ) where {N}
    nthreads = Threads.nthreads()
    batchsize = max(100, min(length(dv) ÷ nthreads, round(Int, sqrt(length(dv)) * 10)))
    @assert N == nthreads "`nthreads()` not matching dimension of `ws`"
    v = localpart(dv)
    stat_names, stats = step_stats(v, Val(N))
    zero!.(ws)
    @sync for btr in Iterators.partition(pairs(v), batchsize)
        Threads.@spawn for (add, num) in btr
            ss = fciqmc_col!(ws[Threads.threadid()], ham, add, num, shift, dτ)
            stats[Threads.threadid()] += ss
        end
    end
    return stat_names, stats
end

"""
    SplittablesThreading <: ThreadingStrategy

[`ThreadingStrategy`](@ref) based on [SplittablesBase](https://github.com/JuliaFolds/SplittablesBase.jl).
"""
struct SplittablesThreading <: ThreadingStrategy end

working_memory(::SplittablesThreading, dv) = ntuple_working_memory(dv)

@inline function _loop_configs!(ws, stats, ham, pairs, shift, dτ, batchsize)
    if amount(pairs) > batchsize
        two_halves = halve(pairs)
        first_half = @spawn _loop_configs!(
            ws, stats, ham, two_halves[1], shift, dτ, batchsize
        )
        _loop_configs!(
            ws, stats, ham, two_halves[2], shift, dτ, batchsize
        )
        wait(first_half)
    else
        for (add, num) in pairs
            stats[threadid()] += fciqmc_col!(ws[threadid()], ham, add, num, shift, dτ)
        end
    end
    return nothing
end

function fciqmc_step!(::SplittablesThreading, ws::NTuple{N}, ham, dv, shift, dτ) where {N}
    @assert N > 1 "attempted to run threaded code with one thread"
    v = localpart(dv)
    stat_names, stats = step_stats(v, Val(N))
    batchsize = max(100.0, min(amount(pairs(v))/N, sqrt(amount(pairs(v))) * 10))

    zero!.(ws)
    _loop_configs!(ws, stats, ham, pairs(v), shift, dτ, batchsize)

    return stat_names, stats
end

"""
    ThreadsXSumThreading <: ThreadingStrategy

[`ThreadingStrategy`](@ref) based on [`ThreadsX`](https://github.com/tkf/ThreadsX.jl)`.sum`.
"""
struct ThreadsXSumThreading <: ThreadingStrategy end

working_memory(::ThreadsXSumThreading, dv) = ntuple_working_memory(dv)

function fciqmc_step!(::ThreadsXSumThreading, ws::NTuple{N}, ham, dv, shift, dτ) where {N}
    @assert N == Threads.nthreads() "`nthreads()` not matching dimension of `ws`"
    v = localpart(dv)
    zero!.(ws)

    stat_names, _ = step_stats(v, Val(1))

    stats = ThreadsX.sum(pairs(v)) do (add, val)
        MultiScalar(fciqmc_col!(ws[Threads.threadid()], ham, add, val, shift, dτ))
    end
    return stat_names, stats
end

"""
    ThreadsXForeachThreading <: ThreadingStrategy

[`ThreadingStrategy`](@ref) based on [`ThreadsX`](https://github.com/tkf/ThreadsX.jl)`.foreach`.
"""
struct ThreadsXForeachThreading <: ThreadingStrategy end

working_memory(::ThreadsXForeachThreading, dv) = ntuple_working_memory(dv)

function fciqmc_step!(::ThreadsXForeachThreading, ws::NTuple{N}, ham, dv, shift, dτ) where {N}
    # multithreaded version; should also work with MPI
    @assert N == Threads.nthreads() "`nthreads()` not matching dimension of `ws`"
    v = localpart(dv)

    stat_names, stats = step_stats(v, Val(N))
    zero!.(ws) # clear working memory

    # parallel execution happens here:
    ThreadsX.foreach(pairs(v)) do (add, val)
        tid = Threads.threadid()
        stats[tid] += fciqmc_col!(ws[tid], ham, add, val, shift, dτ)
    end
    return stat_names, stats
end

"""
    select_threading_strategy(threading, targetwalkers)

Select a [`ThreadingStrategy`](@ref) to control threading in [`lomc!`](@ref).

`threading` can be:

* `:auto`: decide whether threading should be done or not based on `targetwalkers` and
  whether threads are available.
* `true` or `false`: use the default [`ThreadingStrategy`](@ref) or [`NoThreading`](@ref).
* Any [`ThreadingStrategy`](@ref).

The default [`ThreadingStrategy`](@ref) is currently [`SplittablesThreading`](@ref).
"""
select_threading_strategy(t::ThreadingStrategy, _) = t

function select_threading_strategy(threading::Symbol, targetwalkers)
    if threading == :auto
        t = targetwalkers > 500 && Threads.nthreads() > 1
        return select_threading_strategy(t, targetwalkers)
    else
        error("invalid threading strategy `$threading`")
    end
end
function select_threading_strategy(threading::Bool, _)
    if threading
        if Threads.nthreads() == 1
            @warn "threading was requested, but only one thread is available"
            return NoThreading()
        else
            return SplittablesThreading()
        end
    else
        return NoThreading()
    end
end


text:
import Pkg
import UUIDs

@inline function _convert_to_package_spec(versioned_package::VersionedPackage)
    name = versioned_package.name::String
    uuid = versioned_package.uuid::UUIDs.UUID
    version = versioned_package.version::String
    tree = versioned_package.tree::String
    package_spec = Pkg.API.Package(
        name = name,
        uuid = uuid,
        version = version,
    )::Pkg.Types.PackageSpec
    return package_spec
end

@inline function Pkg.Types.PackageSpec(versioned_package::VersionedPackage)
    package_spec = _convert_to_package_spec(versioned_package)::Pkg.Types.PackageSpec
    return package_spec
end

# Note: Pkg.API.Package == Pkg.PackageSpec

@inline function Pkg.API.Package(versioned_package::VersionedPackage)
    package_spec = _convert_to_package_spec(versioned_package)::Pkg.Types.PackageSpec
    return package_spec
end

@inline function Base.convert(::Type{P}, versioned_package::VersionedPackage) where P <: Pkg.Types.PackageSpec
    package_spec = _convert_to_package_spec(versioned_package)::Pkg.Types.PackageSpec
    return package_spec
end

# Note: Pkg.API.add == Pkg.add

@inline function Pkg.API.add(versioned_package::VersionedPackage)
    package_spec = _convert_to_package_spec(versioned_package)::Pkg.Types.PackageSpec
    return Pkg.API.add(package_spec)
end

@inline function Pkg.API.add(versioned_packages::AbstractVector{<:VersionedPackage})
    package_specs = _convert_to_package_spec.(versioned_packages)
    return Pkg.API.add(package_specs)
end

# Note: Pkg.API.pin == Pkg.pin

@inline function Pkg.API.pin(versioned_package::VersionedPackage)
    package_spec = _convert_to_package_spec(versioned_package)::Pkg.Types.PackageSpec
    return Pkg.API.pin(package_spec)
end

@inline function Pkg.API.pin(versioned_packages::AbstractVector{<:VersionedPackage})
    package_specs = _convert_to_package_spec.(versioned_packages)
    return Pkg.API.pin(package_specs)
end


text:
using BFREPL
using REPL
using Test

struct MockREPL
    t::REPL.Terminals.TTYTerminal
end
(::Type{MockREPL})(inio::IO, outio::IO) = MockREPL(inio, outio, outio)
function (::Type{MockREPL})(inio::IO, outio::IO, errio::IO)
    MockREPL(REPL.Terminals.TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), inio, outio, errio))
end
REPL.outstream(m::MockREPL) = m.t

@testset "Hello, world!" begin
    inio = IOBuffer()
    outio = IOBuffer()
    repl = MockREPL(inio, outio)

    hello_world_code = """
    +++++++++[>++++++++>+++++++++++>+++++<<<-]>.>++.+++++++..+++.>-.
    ------------.<++++++++.--------.+++.------.--------.>+.
    """
    BFREPL.do_cmd(repl, hello_world_code)
    result = String(take!(outio))
    @test result == "Hello, world!"

    I_♥_Julia_code = """
    >++[>++[>+++[>+++>++++++>+<<<-]<-]<-]>>>>+.<----.[<-<---<---<+>>>>-]<++.++++
    <-------.<+++++.<.>>>>>+.<<<<<[>-->->>>+>+<<<<<<-]>>----.>>>++.---.<<<<----.>>>>>++.
    """
    BFREPL.do_cmd(repl, I_♥_Julia_code)
    result = String(take!(outio))
    @test result == "I ♥ Julia."
end

@testset "lc_echo" begin
    lc_echo_code = """
    >>,+[-[-<+<+>>]++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    +<[->->+<[>-]>[>]<[-<<[->+<]++++++++++++++++++++++++++[->>+<[>-]>[>]<[-<<<+++++
    +++++++++++++++++++++++++++>[-]>>]<-<]>>]<<]<.[-]>>,+]
    """
    outio = IOBuffer()

    A2Z_repl = MockREPL(IOBuffer("ABCDEFGHIJKLMNOPQRSTUVWXYZ"), outio)
    BFREPL.do_cmd(A2Z_repl, lc_echo_code)
    result = String(take!(outio))
    @test result == "abcdefghijklmnopqrstuvwxyz"

    I_♥_Julia_repl = MockREPL(IOBuffer("I ♥ Julia."), outio)
    BFREPL.do_cmd(I_♥_Julia_repl, lc_echo_code)
    result = String(take!(outio))
    @test result == "i ♥ julia."
end

@testset "error handlings" begin
    inio = IOBuffer()
    outio = IOBuffer()
    errio = IOBuffer()
    repl = MockREPL(inio, outio, errio)

    BFREPL.do_cmd(repl, "[]]")
    result = String(take!(errio))
    @test startswith(result, "ERROR: syntax: ']' mismatch")

    BFREPL.do_cmd(repl, "[[]")
    result = String(take!(errio))
    @test startswith(result, "ERROR: syntax: '[' mismatch")

    BFREPL.do_cmd(repl, "<+")
    result = String(take!(errio))
    @test startswith(result, "ERROR: BoundsError:")

    BFREPL.do_cmd(repl, "<-")
    result = String(take!(errio))
    @test startswith(result, "ERROR: BoundsError:")

    BFREPL.do_cmd(repl, "<.")
    result = String(take!(errio))
    @test startswith(result, "ERROR: BoundsError:")

    BFREPL.do_cmd(MockREPL(IOBuffer("x"), outio, errio), "<,")
    result = String(take!(errio))
    @test startswith(result, "ERROR: BoundsError:")

    BFREPL.do_cmd(repl, "+[<]")
    result = String(take!(errio))
    @test startswith(result, "ERROR: BoundsError:")

    BFREPL.do_cmd(repl, "+[>+]")
    result = String(take!(errio))
    @test startswith(result, "ERROR: BoundsError:")
end


text:
Base.@propagate_inbounds function setindex(args...)
    Base.setindex(args...)
end

@inline setindex(::Base.RefValue, val) = Ref(val)

Base.@propagate_inbounds function setindex(xs::AbstractArray, v, I...)
    T = promote_type(eltype(xs), typeof(v))
    ys = similar(xs, T)
    if eltype(xs) !== Union{}
        copy!(ys, xs)
    end
    ys[I...] = v
    return ys
end

Base.@propagate_inbounds function setindex(d0::AbstractDict, v, k)
    K = promote_type(keytype(d0), typeof(k))
    V = promote_type(valtype(d0), typeof(v))
    d = empty(d0, K, V)
    copy!(d, d0)
    d[k] = v
    return d
end


text:
# Statistics code

export ShiftedExponential, ExponentialOrder, OrderStatistic, NonIDOrderStatistic

"""
    ExponentialOrder(scale::Real, total::Int, order::Int)

Random variable representing the order-th largest value out of total
realizations of an exponential random variable with given scale.

"""
function ExponentialOrder(scale::Real, total::Int, order::Int)
    scale > 0 || throw(DomainError(scale, "scale must be positive"))
    total > 0 || throw(DomainError((total, order), "total must be positive"))
    order > 0 || throw(DomainError((total, order), "order must be positive"))
    order <= total || throw(DomainError((total, order), "order must be <= total"))
    var = sum(1/(i^2)  for i=(total-order+1):total) * scale^2
    mean = sum(1/i for i=(total-order+1):total) * scale
    alpha = mean^2 / var # shape parameter
    theta = var / mean # scale
    return Gamma(alpha, theta)
end

### Shifted exponential distribution

struct ShiftedExponential{T<:Real} <: ContinuousUnivariateDistribution
    shift::T
    exponential::Exponential{T}
    ShiftedExponential{T}(s::T, θ::T) where T = new{T}(s, Exponential(θ))
end

function ShiftedExponential(s::T, θ::T) where {T<:Real}
    θ > zero(θ) || throw(DomainError(θ), "θ must be positive")
    ShiftedExponential{T}(s, θ)
end

function ShiftedExponential(s::T, θ::T) where {T<:Integer}
    ShiftedExponential(float(s), float(θ))
end

struct ShiftedExponentialStats <: SufficientStats
    vx::Float64     # variance of x
    sx::Float64     # (weighted) sum of x
    sw::Float64     # sum of sample weights
end

Distributions.suffstats(::Type{<:ShiftedExponential}, x::AbstractArray{T}) where {T<:Real} = ShiftedExponentialStats(var(x), sum(x), length(x))

function Distributions.fit_mle(::Type{<:ShiftedExponential}, ss::ShiftedExponentialStats)
    θ = sqrt(ss.vx)
    s = ss.sx / ss.sw - θ
    ShiftedExponential(s, θ)
end

Distributions.quantile(d::ShiftedExponential, p::Real) = Distributions.quantile(d.exponential, p) + d.shift

Distributions.pdf(d::ShiftedExponential, p::Real) = pdf(d.exponential, p-d.shift)
Distributions.cdf(d::ShiftedExponential, p::Real) = cdf(d.exponential, p-d.shift)
Base.rand(rng::AbstractRNG, d::ShiftedExponential) = rand(rng, d.exponential) + d.shift
Base.minimum(d::ShiftedExponential) = d.shift
Base.maximum(::ShiftedExponential) = Inf
Statistics.mean(d::ShiftedExponential) = mean(d.exponential) + d.shift
Statistics.var(d::ShiftedExponential) = var(d.exponential)
Distributions.params(d::ShiftedExponential) = (d.shift, params(d.exponential)...)

### Fitting a distribution from its mean and variance

function distribution_from_mean_variance(::Type{ShiftedExponential}, m, v)
    θ = sqrt(v)
    s = m .- θ
    ShiftedExponential(s, θ)
end

function distribution_from_mean_variance(::Type{Gamma}, m, v)
    θ = v / m
    α = m / θ
    Gamma(α, θ)
end

### Order statistics sampling

struct OrderStatistic{S<:Union{Discrete,Continuous},Spl<:Sampleable{Univariate,S},T} <: Sampleable{Univariate,S}
    spl::Spl
    k::Int
    buffer::Vector{T}
end

Base.show(io::IO, s::OrderStatistic) = print(io, "OrderStatistic($(s.spl), k=$(s.k), n=$(length(s.buffer)))")

OrderStatistic(s::Sampleable, k::Integer, n::Integer) = OrderStatistic(s, k, Vector{eltype(s)}(undef, n))

function Random.rand(rng::AbstractRNG, s::OrderStatistic)
    Distributions.rand!(rng, s.spl, s.buffer)
    partialsort!(s.buffer, s.k)
    s.buffer[s.k]
end

struct NonIDOrderStatistic{S<:Union{Discrete,Continuous},Spl<:Sampleable{Univariate,S},T} <: Sampleable{Univariate,S}
    spls::Vector{Spl}
    k::Int
    buffer::Vector{T}
end

Base.show(io::IO, s::NonIDOrderStatistic) = print(io, "NonIDOrderStatistic($(eltype(s.spls)), k=$(s.k), n=$(length(s.buffer)))")

NonIDOrderStatistic(spls::AbstractVector{<:Sampleable}, k::Integer) = NonIDOrderStatistic(spls, k, Vector{promote_type((eltype(spl) for spl in spls)...)}(undef, length(spls)))

function Random.rand(rng::AbstractRNG, s::NonIDOrderStatistic)
    for (i, spl) in enumerate(s.spls)
        s.buffer[i] = Distributions.rand(rng, spl)
    end
    partialsort!(s.buffer, s.k)
    s.buffer[s.k]
end


text:
p_cstdout = unsafe_load(cglobal(:stdout, Ptr{Void}))

function Matrix_Read_Array{T<:Integer}(arr::Array{T, 2})
	nrow, ncol = size(arr)
	pmat = Matrix_Alloc(Cuint(nrow), Cuint(ncol))
	if nrow > 0 && ncol > 0
		m = unsafe_load(pmat)
		a = pointer_to_array(m.p_Init, length(arr))
		a[:] = arr'[:]
	end
	return pmat
end


text:
### utils.jl
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file defines some utility functions.
#
### Code:

export stdscore, weightedmean
using LinearAlgebra

# Standard Score
"""
    stdscore(measure::Measurement, expected_value::Real) -> standard_score

Gives the standard score between a measure, with uncertainty, and its expected
value:

    (measure.val - expected_value)/measure.err
"""
stdscore(a::Measurement, b) = (a.val - b)/(a.err)

"""
    stdscore(measure_1::Measurement, measure_2::Measurement) -> standard_score

Gives the standard score between two measurements (both with uncertainty)
computed as the standard score between their difference and 0:

    stdscore(measure_1 - measure_2, 0)
"""
stdscore(a::Measurement{S}, b::Measurement{T}) where {S<:AbstractFloat,T<:AbstractFloat} =
    stdscore(a - b, zero(promote_type(S, T)))

# Weighted Average with Inverse-Variance Weighting
"""
    weightedmean(iterable) -> measurement(weighted_mean, standard_deviation)

Return the weighted mean of measurements listed in `iterable`, using
inverse-variance weighting.  NOTA BENE: correlation is not taken into account.
"""
function weightedmean(iterable)
    v = [el.val for el in iterable]
    w = [inv(el.err)^2 for el in iterable]
    invsumw = inv(sum(w))
    return measurement(dot(v, w)*invsumw, sqrt(invsumw))
end

# Derivative and Gradient
derivative(a::Measurement{F},
           tag::Tuple{T, T, UInt64}) where {F<:AbstractFloat, T<:AbstractFloat} =
               get(a.der, tag, zero(F))

"""
    derivative(x::Measurement, y::Measurement)

Return the value of the partial derivative of `x` with respect to the
independent measurement `y`, calculated on the nominal value of `y`.  Return
zero if `x` does not depend on `y`.

Use `Measurements.derivative.(x, array)` to calculate the gradient of `x` with respect to an
array of independent measurements.
"""
derivative(a::Measurement, b::Measurement) =
    derivative(a, (b.val, b.err, b.tag))

# value and uncertainty
for (f, field) in ((:value, :val), (:uncertainty, :err))
    @eval begin
        ($f)(a::Measurement) = a.$field
        ($f)(a::Complex{<:Real}) = complex(($f)(a.re), ($f)(a.im))
    end
end
value(a::Real) = a
uncertainty(a::Real) = zero(a)

"""
    Measurements.value(x::Measurement)
    Measurements.value(x::Complex{Measurement})

Return the nominal value of measurement `x`.
"""
value

"""
    Measurements.uncertainty(x::Measurement)
    Measurements.uncertainty(x::Complex{Measurement})

Return the uncertainty of measurement `x`.
"""
uncertainty

"""
    Measurements.uncertainty_components(x::Measurement)

Return the components to the uncertainty of the dependent quantity `x` in the
form of a `Dict`.
"""
function uncertainty_components(x::Measurement{T}) where {T<:AbstractFloat}
    out = Dict{Tuple{T, T, UInt64}, T}()
    for var in keys(x.der)
        out[var] = abs(var[2] * Measurements.derivative(x, var))
    end
    return out
end


text:
input = joinpath(@__DIR__, "input")
lines = readlines(input)

m = hcat(map(line -> map(n -> parse(Int, n), collect(line)), lines)...)

function neighbors(m, idx)
    x, y = Tuple(idx)
    return filter(c -> checkbounds(Bool, m, c), [
        CartesianIndex(x + 1, y),
        CartesianIndex(x - 1, y),
        CartesianIndex(x, y + 1),
        CartesianIndex(x, y - 1),
        CartesianIndex(x + 1, y + 1),
        CartesianIndex(x - 1, y - 1),
        CartesianIndex(x - 1, y + 1),
        CartesianIndex(x + 1, y - 1)
    ])
end

function step(m, limit)
    counter = 0
    for i in 1:limit
        m = m .+ 1
        octopii = Set()
        flashed = Set()
        for idx in CartesianIndices(m)
            push!(octopii, idx)
            while !isempty(octopii)
                o = pop!(octopii)
                if m[o] > 9 && !(o in flashed)
                    push!(flashed, o)
                    ns = neighbors(m, o)
                    for n in ns
                        m[n] = m[n] + 1
                        push!(octopii, n)
                    end
                end
            end
        end
        for f in flashed
            m[f] = 0
        end
        if length(flashed) == 100
            return (i, counter)
        end
        counter += length(flashed)
    end
    return (limit, counter)
end

(_, p1) = step(m, 100)
(p2, _) = step(m, 1000)

println("-----------------------------------------------------------------------")
println("dumbo octopus -- part one :: $p1")
println("dumbo octopus -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 1632)
@assert(p2 == 303)


text:
module Day17

struct TargetArea
    x::UnitRange{Int}
    y::UnitRange{Int}
end

function parse_target(input::String)
    input = split(input, ": ")[2]
    xpart, ypart = split(input, ", ")
    xpart = split(split(xpart, "=")[2], "..")
    ypart = split(split(ypart, "=")[2], "..")

    x = parse(Int, xpart[1]):parse(Int, xpart[2])
    y = parse(Int, ypart[1]):parse(Int, ypart[2])

    TargetArea(x, y)
end

function part1(target::TargetArea)
    S = abs(target.y[1])
    Int((S - 1) * S / 2)
end

function match_x_in_steps(range::UnitRange{Int}, n)
    # for 2 steps,
    # l <= v + v-1 <= r
    # for 3 steps,
    # l <= v + v-1 + v-2 <= r --> l+3 <= 3v <= r+3 --> l/3+1 <= v <= r/3+1
    l = range[1]
    r = range[end]
    triangle = n * (n-1) / 2
    lower = (l+triangle)/n |> ceil |> Int
    upper = (r + triangle)/n |> floor |> Int
    return lower:upper
end

function match_y_in_steps(range::UnitRange{Int}, n)
    # the problem is symmetric, so only need to consider negative vy.
    # for 1 step,
    # t <= v <= b
    # for 2 steps,
    # t <= v + v+1 <= b
    # for 3 steps,
    # t <= v + v+1 + v+2 <= b
    # ...
    # t-tri <= nv <= b-tri
    t = range[end] |> abs
    b = range[1] |> abs
    triangle = n * (n-1) / 2
    lower = (t - triangle) / n |> ceil |> Int
    upper = (b - triangle) / n |> floor |> Int
    return -(lower:upper)
end

function match(target::TargetArea)
    vxs = Dict()
    stalled_out_vxs = Int[]
    n = 1
    while n*(n+1)/2 <= target.x[end]
        matches = match_x_in_steps(target.x, n)
        push!(vxs, n => Array(matches))
        # For all of the x velocities which arrive in vx steps,
        # they are "stalled out".
        if n in matches
            push!(stalled_out_vxs, n)
        end
        n += 1
    end

    vys = Dict()
    n = 1
    while n*(n+1)/2 <= abs(target.y[1])
        push!(vys, n => Array(match_y_in_steps(target.y, n)))
        n += 1
    end
    i = 0
    k = 0
    for (n, vy) in vys
        vy = vy[vy .< 0]
        parabola_steps = 2abs.(vy) .- 1
        original_vy = -1 .- vy
        for (pn, vy) in zip(parabola_steps, original_vy)
            if (pn+n) ∈ keys(vys)
                push!(vys[pn+n], vy)
            else
                push!(vys, (pn+n) => [vy])
            end
        end
        i += 1
    end

    velocities = Set()
    # Need to find a solution for each vy
    for (n, vy) in vys
        if n ∈ keys(vxs)
            n_step_velocities = tuple.(vxs[n], Array(vys[n])') |> collect
            union!(velocities, vec(n_step_velocities))
        else
            stalled_out_matches = stalled_out_vxs[stalled_out_vxs .<= n]
            vs = tuple.(stalled_out_vxs, Array(vys[n])') |> collect
            union!(velocities, vs)
        end
    end
    velocities
end

function run()
    test = "target area: x=20..30, y=-10..-5"
    testarea = parse_target(test)

    @assert part1(testarea) == 45
    @assert length(match(testarea)) == 112

    input = "target area: x=155..182, y=-117..-67"
    target = parse_target(input)
    display(part1(target))
    display(length(match(target)))
end

function test_part2()
    expected = """23,-10  25,-9   27,-5   29,-6   22,-6   21,-7   9,0     27,-7   24,-5
25,-7   26,-6   25,-5   6,8     11,-2   20,-5   29,-10  6,3     28,-7
8,0     30,-6   29,-8   20,-10  6,7     6,4     6,1     14,-4   21,-6
26,-10  7,-1    7,7     8,-1    21,-9   6,2     20,-7   30,-10  14,-3
20,-8   13,-2   7,3     28,-8   29,-9   15,-3   22,-5   26,-8   25,-8
25,-6   15,-4   9,-2    15,-2   12,-2   28,-9   12,-3   24,-6   23,-7
25,-10  7,8     11,-3   26,-7   7,1     23,-9   6,0     22,-10  27,-6
8,1     22,-8   13,-4   7,6     28,-6   11,-4   12,-4   26,-9   7,4
24,-10  23,-8   30,-8   7,0     9,-1    10,-1   26,-5   22,-9   6,5
7,5     23,-6   28,-10  10,-2   11,-1   20,-9   14,-2   29,-7   13,-3
23,-5   24,-8   27,-9   30,-7   28,-5   21,-10  7,9     6,6     21,-5
27,-10  7,2     30,-9   21,-8   22,-7   24,-9   20,-6   6,9     29,-5
8,-2    27,-8   30,-5   24,-7"""
    expected = split(expected) |> Set

    test = "target area: x=20..30, y=-10..-5"
    testarea = parse_target(test)
    velocities = match(testarea)
    actual = [join(v, ",") for v in velocities] |> Set

    setdiff(actual, expected) |> display
    setdiff(expected, actual) |> display
end

end


text:
using Test

include("test_env.jl")
include("test_agent.jl")
include("test_episode.jl")
# include("test_hashlesstilecoder.jl")

function runtests()
    test_env()
    test_agent()
    test_episode()
    # test_hashlesstilecoder()
end

runtests()


text:
#<----Mean---->
"""
    ## Mean
    ### Description
      Returns the mean of an array.\n
      --------------------\n
    ### Input
      mean(x)\n
      --------------------\n
      #### Positional Arguments
      Array{Any} - x:: Array to obtain the mean of.
      --------------------\n
     ### Output
     mu:: The mean of the provided array.
       """
mean(x) = sum(x) / length(x)
#<----Median---->
"""
      Calculates the median (numerical center) of a given array.\n
      --------------------\n
      array = [5,10,15]\n
      median = Lathe.stats.median(array)\n
      println(median)\n
        10
       """
function median(array)
    n = length(array)
    half = n / 2
    current = 1
    sorted = sort!(array,rev = false)
    median = 0
    for i in sorted
        if current >= half
            median = i
        else
        current = current + 1
    end
    end
    return(median)
end
#<----Mode---->
"""
      Gives the digit most common in a given array\n
      --------------------\n
      array = [5,10,15,15,10,5,10]\n
      mode = Lathe.stats.mode(array)\n
      println(mode)\n
        10
       """
function mode(array)
    m = findmax(array)
    return(m)
end
#<----Variance---->
"""
      Gives the variance of an array..\n
      --------------------\n
      array = [5,10,15]\n
      variance = Lathe.stats.variance(array)\n
       """
function variance(array)
    me = mean(array)
    sq = sum(array) - me
    squared_mean = sq ^ 2
    return(squared_mean)
end
#<----Confidence Intervals---->
"""
      Returns the confidence intervals of given data.\n
      --------------------\n
      array = [5,10,15]\n
      confidence = .98\n
      low, high = Lathe.stats.confints(array,confidence)\n
      --------------------\n
      PARAMETERS\n
      confidence: Confidence is a float percentage representing the level of
      confidence required in your test. The confidence metric is used
      exclusively for calculating the interval.
       """
function confints(data, confidence=.95)
    mean = stats.mean(data)
    std = std(data)
    stderr = ste(data)
#    interval = stderr * scs.t.ppf((1 + confidence) / 2.0, n-1)
#    return (mean-interval, mean+interval)
end
#<----Standard Error---->
@doc """
      Calculates the Standard Error of an array.\n
      --------------------\n
      array = [5,10,15]\n
      ste = Lathe.stats.standarderror(array)"""
function ste(data)
    std = std(data)
    sample = length(data)
    ste = (std/sqrt(sample))
    return(ste)
end
#<----Standard Deviation---->
"""
      Calculates the Standard Deviation of a given array.\n
      --------------------\n
      array = [5,10,15]\n
      std = Lathe.stats.std(array)\n
       """
function std(array3)
    m = mean(array3)
    [i = (i-m) ^ 2 for i in array3]
    m = mean(array3)
    try
        m = sqrt(m)
    catch
        m = sqrt(Complex(m))
    end
    return(m)
end
#<----Quartiles---->
# - First
"""
      Returns the point in an array located at 25 percent of the sorted data.\n
      --------------------\n
      array = [5,10,15]\n
      q1 = Lathe.stats.firstquar(array)\n
       """
function q1(array)
    m = median(array)
    q1 = array * .5
    return(q1)
end
# - Third
"""
      Returns the point in an array located at 75 percent of the sorted data.\n
      --------------------\n
      array = [5,10,15]\n
      q2 = Lathe.stats.secondquar(array)\n
       """
function q3(array)
    q = median(array)
    q3 = q * 1.5
    return(q3)
end
# <---- Rank ---->
"""
      Ranks indices in an array based on quantitative weight (count of the
      numbers) and returns a new array of the ranks of each column. This
      function is made primarily for the Wilcox Rank-Sum test.\n
      --------------------\n
      array = [5,10,15]\n
      q2 = Lathe.stats.secondquar(array)\n
       """
function getranks(array,rev = false)
    sortedar = sort(array,rev=rev)
    num = 1
    list = []
    for i in sortedar
        append!(list,i)
        num = num + 1
    end
    return(list)
end
"""
    ## Factorials
    ### Description
      Calculates the factorial of a number.\n
      --------------------\n
    ### Input
      fact(n)\n
      --------------------\n
      #### Positional Arguments
      Int64 - n:: The number for the factorial.\n
      --------------------\n
     ### Output
     f:: Factorial of n
       """
function fact(n)
    if n == 1
        return(1)
    else
        return n * fact(n-1)
    end
end

is_prime(n) = φ(n) == n - 1


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Variables, constraints, and objective functions

# **Originally Contributed by**: Arpit Bhatia

# While the last tutorial introduced you to basics of of JuMP code, this
# tutorial will go in depth focusing on how to work with different parts of a
# JuMP program.

using JuMP
model = Model();

# ## Variable basics

# ### Variable bounds

# All of the variables we have created till now have had a bound. We can also
# create a free variable.

@variable(model, free_x)

# While creating a variable, instead of using the <= and >= syntax, we can also
# use the `lower_bound` and `upper_bound` keyword arguments.

@variable(model, keyword_x, lower_bound = 1, upper_bound = 2)

# We can query whether a variable has a bound using the `has_lower_bound` and
# `has_upper_bound` functions. The values of the bound can be obtained using the
# `lower_bound` and `upper_bound` functions.

has_upper_bound(keyword_x)

#-

upper_bound(keyword_x)

# Note querying the value of a bound that does not exist will result in an error.

try #hide
lower_bound(free_x)
catch err; showerror(stderr, err); end  #hide

# JuMP also allows us to change the bounds on variable. We will learn this in
# the problem modification tutorial.

# ### [Containers](@id tutorial_variable_container)

# We have already seen how to add a single variable to a model using the
# [`@variable`](@ref) macro. Let's now look at more ways to add variables to a
# JuMP model.

# JuMP provides data structures for adding collections of variables to a model.
# These data structures are referred to as Containers and are of three types:
# `Arrays`, `DenseAxisArrays`, and `SparseAxisArrays`.

# #### Arrays

# JuMP arrays are created in a similar syntax to Julia arrays with the addition
# of specifying that the indices start with 1. If we do not tell JuMP that the
# indices start at 1, it will create a `DenseAxisArray` instead.

@variable(model, a[1:2, 1:2])

# An n-dimensional variable $x \in {R}^n$ having a bound $l \preceq x \preceq u$
# ($l, u \in {R}^n$) is added in the following manner.

n = 10
l = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10]
u = [10; 11; 12; 13; 14; 15; 16; 17; 18; 19]

@variable(model, l[i] <= x[i = 1:n] <= u[i])

# Note that while working with Containers, we can also create variable bounds
# depending upon the indices:

@variable(model, y[i = 1:2, j = 1:2] >= 2i + j)

# #### DenseAxisArrays

# `DenseAxisArrays` are used when the required indices are not one-based integer
# ranges. The syntax is similar except with an arbitrary vector as an index as
# opposed to a one-based range.

# An example where the indices are integers but do not start with one.

@variable(model, z[i = 2:3, j = 1:2:3] >= 0)

# Another example where the indices are an arbitrary vector.

@variable(model, w[1:5,["red", "blue"]] <= 1)

# #### SparseAxisArrays

# `SparseAxisArrays` are created when the indices do not form a rectangular set.
# For example, this applies when indices have a dependence upon previous indices
# (called triangular indexing).

@variable(model, u[i = 1:3, j = i:5])

# We can also conditionally create variables by adding a comparison check that
# depends upon the named indices and is separated from the indices by a
# semi-colon (;).

@variable(model, v[i = 1:9; mod(i, 3) == 0])

# ### Variable types

# The last arguement to the `@variable` macro is usually the variable type. Here
# we'll look at how to specifiy the variable type.

# #### Integer variables

# Integer optimization variables are constrained to the set $x \in {Z}$

@variable(model, integer_x, Int)

# or

@variable(model, integer_z, integer = true)

# #### Binary variables

# Binary optimization variables are constrained to the set $x \in \{0, 1\}$.

@variable(model, binary_x, Bin)

# or

@variable(model, binary_z, binary = true)

# #### Semidefinite variables

# JuMP also supports modeling with semidefinite variables. A square symmetric
# matrix X is positive semidefinite if all eigenvalues are nonnegative.

@variable(model, psd_x[1:2, 1:2], PSD)

# We can also impose a weaker constraint that the square matrix is only
# symmetric (instead of positive semidefinite) as follows:

@variable(model, sym_x[1:2, 1:2], Symmetric)

# ## Constraint basics

model = Model()
@variable(model, x)
@variable(model, y)
@variable(model, z[1:10]);

# ### Constraint references

# While calling the `@constraint` macro, we can also set up a constraint
# reference. Such a referrence is useful for obtaining additional information
# about the constraint, such as its dual solution.

@constraint(model, con, x <= 4)

# ### [Containers](@id tutorial_constraint_container)

# Just as we had containers for variables, JuMP also provides `Arrays`,
# `DenseAxisArrays`, and `SparseAxisArrays` for storing collections of
# constraints. Examples for each container type are given below.

# #### Arrays

@constraint(model, [i = 1:3], i * x <= i + 1)

# #### DenseAxisArrays

@constraint(model, [i = 1:2, j = 2:3], i * x <= j + 1)

# #### SparseAxisArrays

@constraint(model, [i = 1:2, j = 1:2; i != j], i * x <= j + 1)

# ### Constraints in a loop

# We can add constraints using regular Julia loops

for i in 1:3
    @constraint(model, 6x + 4y >= 5i)
end

# or use for each loops inside the `@constraint` macro.

@constraint(model, [i in 1:3], 6x + 4y >= 5i)

# We can also create constraints such as $\sum _{i = 1}^{10} z_i \leq 1$

@constraint(model, sum(z[i] for i in 1:10) <= 1)

# ## Objective functions

# While the recommended way to set the objective is with the [`@objective`](@ref)
# macro, the functions [`set_objective_sense`](@ref) and [`set_objective_function`](@ref)
# provide an equivalent lower-level interface.

using GLPK

model = Model(GLPK.Optimizer)
@variable(model, x >= 0)
@variable(model, y >= 0)
set_objective_sense(model, MOI.MIN_SENSE)
set_objective_function(model, x + y)

optimize!(model)

#-

objective_value(model)

# To query the objective function from a model, we use the [`objective_sense`](@ref),
# [`objective_function`](@ref), and [`objective_function_type`](@ref) functions.

objective_sense(model)

#-

objective_function(model)

#-

objective_function_type(model)

# # Vectorized syntax

# We can also add constraints and an objective to JuMP using vectorized linear
# algebra. We'll illustrate this by solving an LP in standard form i.e.

# ```math
# \begin{aligned}
# & \min & c^T x \\
# & \;\;\text{s.t.} & A x = b \\
# & & x \succeq 0 \\
# & & x \in \mathbb{R}^n
# \end{aligned}
# ```

vector_model = Model(GLPK.Optimizer)

A= [ 1 1 9  5;
     3 5 0  8;
     2 0 6 13]

b = [7; 3; 5]

c = [1; 3; 5; 2]

@variable(vector_model, x[1:4] >= 0)
@constraint(vector_model, A * x .== b)
@objective(vector_model, Min, c' * x)

optimize!(vector_model)

#-

objective_value(vector_model)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export generateSignal, stretchData, addNoise, addNoise!

using Random, DSP, Dierckx, FileIO, JLD2, SeisIO, FFTW

include("waves.jl")

# This script generates or reads a synthetic time series, stretches it by some
# value, and (optionally) adds noise to the signal.

"""

    generateSignal(type::String; params::Dict{String,Real}=Dict(), sparse::Int64=0, seed::Int64=66743)

Generate a signal of type "ricker" or "dampedSinusoid".

# Arguments
- `type::String,`    : Type of signal to generate: "ricker" or "dampedSinusoid"
- `params::Dict{String,Real}`    : parameters used to define the Ricker wavelet or damped sinusoid
- `sparse::Int64`    : Interval over which to keep elements in array, e.g. sparse=4 corresponds to keeping every 4th element
- `seed::Int64`    : Seed for random number generator

# Output
- `u0`::Array{Float64,1}    : Computed signal
- `t`::Array{Float64,1}    : Time axis

"""
function generateSignal(type::String, params::Dict{String,Real}; sparse::Int64=0, seed::Int64=66473, stretchSource::Float64=0.0)
    if type=="ricker"
        # unpack parameters
        f    = params["f"]    # peak frequency
        dt   = params["dt"]   # sampling interval
        npr  = params["npr"]  # number of points in the Ricker wavelet
        npts = params["npts"] # number of points in the generated signal
        m    = params["m"]    # power for reflectivity series

        # generate a ricker wavelet
        (w, t) = ricker(f=f, n=npr, dt=dt)

        if stretchSource != 0.0
            spectrum = rfft(w)
            stAmpSpec, st = stretchData(real(spectrum), dt, stretchSource)
            spectrum = stAmpSpec .+ imag(spectrum).*im
            w = irfft(spectrum, length(w))
        end

        # generate a random reflectivity series
        rng = MersenneTwister(seed)
        f = randn(rng, Float64, npts)

        if sparse!=0
            r=collect(1:npts)
            f[r .% sparse .!= 0] .= 0
        end
        # raise output from randn to an integral power. See <<Numerical Methods of
        # Explorational Seismology>> page 77.
        f = sign.(f) .* abs.(f).^m

        # convolve reflectivity series with ricker wavelet
        u0 = conv(f,w)[convert(Int64, npts-floor(npts/2)):convert(Int64, npts+floor(npts/2))]

    elseif type=="dampedSinusoid"
        # unpack parameters
        A    = params["A"]    # amplitude
        ω    = params["ω"]    # angular frequency of the sinusoid
        ϕ    = params["ϕ"]    # phase shift for the sinusoid
        η    = params["η"]    # shift in function maximum
        λ    = params["λ"]    # decay constant for the exponential decay
        dt   = params["dt"]   # sampling interval
        t0   = params["t0"]   # start time
        npts = params["npts"] # number of points in the generated signal

        u0, t = dampedSinusoid(A=A, ω=ω, ϕ=ϕ, η=η, n=npts, dt=dt, t0=t0, λ=λ)

    elseif type=="sinc"
        # unpack parameters
        A    = params["A"]
        ω    = params["ω"]    # angular frequency of the sinusoid
        ϕ    = params["ϕ"]    # phase shift for the sinusoid
        dt   = params["dt"]   # sampling interval
        t0   = params["t0"]   # start time
        npts = params["npts"] # number of points in the generated signal

        u0, t = sinc(A=A, ω=ω, ϕ=ϕ, dt=dt, t0=t0, n=npts)

    elseif type=="chirp"
        # unpack parameters
        c     = params["c"]     # phase velocity
        tp    = params["tp"]    # period
        mintp = params["mintp"] # minimum period
        maxtp = params["maxtp"] # maximum period
        dist  = params["dist"]  # distance
        npts  = params["npts"]  # number of points in the generated signal
        dt    = params["dt"]    # sampling interval
        t0    = params["t0"]    # start time

        u0, t = chirp(c=c, tp=tp, mintp=mintp, maxtp=maxtp, dist=dist, n=npts, dt=dt, t0=t0)

    elseif type=="spectSynth"
        # unpack parameters
        A    = params["A"]
        dt   = params["dt"]
        t0   = params["t0"]
        npts = params["npts"]

        u0, t = spectSynth(A=A, dt=dt, t0=t0, n=npts)
    end

    return u0, t
end

"""

    stretchData(u0::Array{Float64,1}, dt::Float64, dvV::Float64; n::Float64=0.0, seed::Int64=66743)

Linearly stretch a given signal, u0, by some factor given by a homogenous relative velocity change dvV.

# Arguments
- `u0::Array{Float64,1},`    : Input signal to be stretched
- `dt::Float64`    : Sampling interval
- `dvV::Float64`    : Stretching factor's causal relative velocity change
- `n::Float64`    : Noise level as a scalar of maximum time
- `seed::Int64`    : Seed for random number generator

# Output
- `u1`::Array{Float64,1}    : Stretched signal
- `tvec`::Array{Float64,1}    : Unstretched time axis

"""
function stretchData(u0::Union{Array{Float32,1},Array{Float64,1}}, dt::Float64, dvV::Float64; starttime::Float64=0.0, stloc::Float64=0.0, n::Float64=0.0, seed::Int64=66473)
    tvec = collect((0:length(u0)-1) .* dt) .+ starttime
    st = (tvec.-stloc) .* dvV

    if (n != 0.0) addNoise!(st, n, seed=seed) end

    tvec2 = tvec .+ st

    # interpolation
    spl = Spline1D(tvec, u0)
    u1 = spl(tvec2)

    return u1, tvec
end

"""

    addNoise!(signal::Array{Float64,1}, level::Float64; seed::Int64=66743)

Add noise to an array by some given mulitple of the maximum value.

# Arguments
- `signal::Array{Float64,1},`    : Input signal to add noise to
- `level::Float64`    : Multiple of maximum value of signal to set noise amplitude
- `seed::Int64`    : Seed for random number generator

# Output
- `signal`::Array{Float64,1}    : Signal with added noise

"""
function addNoise!(signal::Union{Array{Float32,1},Array{Float64,1}}, level::Float64; percent::Bool=true, seed::Int64=66473, freqmin::Float64=0.1, freqmax::Float64=9.99, fs=20.0, corners::Int64=4, zerophase::Bool=false)
    rng = MersenneTwister(seed)
    noise = randn(rng, Float64, length(signal))
    bandpass!(noise, freqmin, freqmax, fs, corners=corners, zerophase=zerophase)
    noise ./= maximum(abs.(noise)) # normalize noise
    noise .*= level # scale noise to given level
    if percent noise .*= maximum(abs.(signal)) end
    signal .+= noise

    return nothing
end
addNoise(signal::Union{Array{Float32,1},Array{Float64,1}}, level::Float64;
         seed::Int64=66473, freqmin::Float64=0.1, freqmax::Float64=9.99, fs=20.0,
         corners::Int64=4, zerophase::Bool=false) = (U = deepcopy(signal);
         addNoise!(U, level, seed=seed, freqmin=freqmin, freqmax=freqmax, fs=fs,
         corners=corners, zerophase=zerophase); return U)

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
const Idx = Union{Real,Colon,AbstractArray{Int}}

using Base: ViewIndex, @propagate_inbounds, tail

# Defer IndexStyle to the wrapped array
@compat Base.IndexStyle{T,N,D,Ax}(::Type{AxisArray{T,N,D,Ax}}) = IndexStyle(D)

# Simple scalar indexing where we just set or return scalars
@propagate_inbounds Base.getindex(A::AxisArray, idxs::Int...) = A.data[idxs...]
@propagate_inbounds Base.setindex!(A::AxisArray, v, idxs::Int...) = (A.data[idxs...] = v)

# Cartesian iteration
Base.eachindex(A::AxisArray) = eachindex(A.data)

"""
    reaxis(A::AxisArray, I...)

This internal function determines the new set of axes that are constructed upon
indexing with I.
"""
reaxis(A::AxisArray, I::Idx...) = _reaxis(make_axes_match(axes(A), I), I)
# Ensure the number of axes matches the number of indexing dimensions
@inline make_axes_match(axs, idxs) = _make_axes_match((), axs, Base.index_ndims(idxs...))
# Move the axes into newaxes, until we run out of both simultaneously
@inline _make_axes_match(newaxes, axs::Tuple, nidxs::Tuple) =
    _make_axes_match((newaxes..., axs[1]), tail(axs), tail(nidxs))
@inline _make_axes_match(newaxes, axs::Tuple{}, nidxs::Tuple{}) = newaxes
# Drop trailing axes, replacing it with a default name for the linear span
@inline _make_axes_match(newaxes, axs::Tuple, nidxs::Tuple{}) =
    (maybefront(newaxes)..., _nextaxistype(newaxes)(Base.OneTo(length(newaxes[end]) * prod(map(length, axs)))))
# Insert phony singleton trailing axes
@inline _make_axes_match(newaxes, axs::Tuple{}, nidxs::Tuple) =
    _make_axes_match((newaxes..., _nextaxistype(newaxes)(Base.OneTo(1))), (), tail(nidxs))

@inline maybefront(::Tuple{}) = ()
@inline maybefront(t::Tuple) = Base.front(t)

# Now we can reaxis without worrying about mismatched axes/indices
@inline _reaxis(axs::Tuple{}, idxs::Tuple{}) = ()
# Scalars are dropped
const ScalarIndex = @compat Union{Real, AbstractArray{<:Any, 0}}
@inline _reaxis(axs::Tuple, idxs::Tuple{ScalarIndex, Vararg{Any}}) = _reaxis(tail(axs), tail(idxs))
# Colon passes straight through
@inline _reaxis(axs::Tuple, idxs::Tuple{Colon, Vararg{Any}}) = (axs[1], _reaxis(tail(axs), tail(idxs))...)
# But arrays can add or change dimensions and accompanying axis names
@inline _reaxis(axs::Tuple, idxs::Tuple{AbstractArray, Vararg{Any}}) =
    (_new_axes(axs[1], idxs[1])..., _reaxis(tail(axs), tail(idxs))...)

# Vectors simply create new axes with the same name; just subsetted by their value
@inline _new_axes{name}(ax::Axis{name}, idx::AbstractVector) = (Axis{name}(ax.val[idx]),)
# Arrays create multiple axes with _N appended to the axis name containing their indices
@generated function _new_axes{name, N}(ax::Axis{name}, idx::@compat(AbstractArray{<:Any,N}))
    newaxes = Expr(:tuple)
    for i=1:N
        push!(newaxes.args, :($(Axis{Symbol(name, "_", i)})(indices(idx, $i))))
    end
    newaxes
end
# And indexing with an AxisArray joins the name and overrides the values
@generated function _new_axes{name, N}(ax::Axis{name}, idx::@compat(AxisArray{<:Any, N}))
    newaxes = Expr(:tuple)
    idxnames = axisnames(idx)
    for i=1:N
        push!(newaxes.args, :($(Axis{Symbol(name, "_", idxnames[i])})(idx.axes[$i].val)))
    end
    newaxes
end

@propagate_inbounds function Base.getindex(A::AxisArray, idxs::Idx...)
    AxisArray(A.data[idxs...], reaxis(A, idxs...))
end

# To resolve ambiguities, we need several definitions
if VERSION >= v"0.6.0-dev.672"
    using Base.AbstractCartesianIndex
    @propagate_inbounds Base.view(A::AxisArray, idxs::Idx...) = AxisArray(view(A.data, idxs...), reaxis(A, idxs...))
else
    @propagate_inbounds function Base.view{T,N}(A::AxisArray{T,N}, idxs::Vararg{Idx,N})
        AxisArray(view(A.data, idxs...), reaxis(A, idxs...))
    end
    @propagate_inbounds function Base.view(A::AxisArray, idx::Idx)
        AxisArray(view(A.data, idx), reaxis(A, idx))
    end
    @propagate_inbounds function Base.view{N}(A::AxisArray, idxs::Vararg{Idx,N})
        # this should eventually be deleted, see julia #14770
        AxisArray(view(A.data, idxs...), reaxis(A, idxs...))
    end
end

# Setindex is so much simpler. Just assign it to the data:
@propagate_inbounds Base.setindex!(A::AxisArray, v, idxs::Idx...) = (A.data[idxs...] = v)

### Fancier indexing capabilities provided only by AxisArrays ###
@propagate_inbounds Base.getindex(A::AxisArray, idxs...) = A[to_index(A,idxs...)...]
@propagate_inbounds Base.setindex!(A::AxisArray, v, idxs...) = (A[to_index(A,idxs...)...] = v)
# Deal with lots of ambiguities here
if VERSION >= v"0.6.0-dev.672"
    @propagate_inbounds Base.view(A::AxisArray, idxs::ViewIndex...) = view(A, to_index(A,idxs...)...)
    @propagate_inbounds Base.view(A::AxisArray, idxs::Union{ViewIndex,AbstractCartesianIndex}...) = view(A, to_index(A,Base.IteratorsMD.flatten(idxs)...)...)
    @propagate_inbounds Base.view(A::AxisArray, idxs...) = view(A, to_index(A,idxs...)...)
else
    for T in (:ViewIndex, :Any)
        @eval begin
            @propagate_inbounds function Base.view{T,N}(A::AxisArray{T,N}, idxs::Vararg{$T,N})
                view(A, to_index(A,idxs...)...)
            end
            @propagate_inbounds function Base.view(A::AxisArray, idx::$T)
                view(A, to_index(A,idx)...)
            end
            @propagate_inbounds function Base.view{N}(A::AxisArray, idsx::Vararg{$T,N})
                # this should eventually be deleted, see julia #14770
                view(A, to_index(A,idxs...)...)
            end
        end
    end
end

# First is indexing by named axis. We simply sort the axes and re-dispatch.
# When indexing by named axis the shapes of omitted dimensions are preserved
# TODO: should we handle multidimensional Axis indexes? It could be interpreted
#       as adding dimensions in the middle of an AxisArray.
# TODO: should we allow repeated axes? As a union of indices of the duplicates?
@generated function to_index{T,N,D,Ax}(A::AxisArray{T,N,D,Ax}, I::Axis...)
    dims = Int[axisdim(A, ax) for ax in I]
    idxs = Expr[:(Colon()) for d = 1:N]
    names = axisnames(A)
    for i=1:length(dims)
        idxs[dims[i]] == :(Colon()) ||
            return :(throw(ArgumentError(string("multiple indices provided ",
                "on axis ", $(string(names[dims[i]]))))))
        idxs[dims[i]] = :(I[$i].val)
    end

    meta = Expr(:meta, :inline)
    return :($meta; to_index(A, $(idxs...)))
end

function Base.reshape{N}(A::AxisArray, ::Type{Val{N}})
    # axN, _ = Base.IteratorsMD.split(axes(A), Val{N})
    # AxisArray(reshape(A.data, Val{N}), reaxis(A, Base.fill_to_length(axN, :, Val{N})...))
    AxisArray(reshape(A.data, Val{N}), reaxis(A, ntuple(d->Colon(), Val{N})...))
end

### Indexing along values of the axes ###

# Default axes indexing throws an error
axisindexes(ax, idx) = axisindexes(axistrait(ax.val), ax.val, idx)
axisindexes(::Type{Unsupported}, ax, idx) = error("elementwise indexing is not supported for axes of type $(typeof(ax))")
axisindexes(t, ax, idx) = error("cannot index $(typeof(ax)) with $(typeof(idx)); expected $(eltype(ax)) axis value or Integer index")

# Dimensional axes may be indexed directy by their elements if Non-Real and unique
# Maybe extend error message to all <: Numbers if Base allows it?
axisindexes{T<:Real}(::Type{Dimensional}, ax::AbstractVector{T}, idx::T) = error("indexing by axis value is not supported for axes with $(eltype(ax)) elements; use an ClosedInterval instead")
function axisindexes(::Type{Dimensional}, ax::AbstractVector, idx)
    idxs = searchsorted(ax, ClosedInterval(idx,idx))
    length(idxs) > 1 && error("more than one datapoint lies on axis value $idx; use an interval to return all values")
    idxs[1]
end

# Dimensional axes may be indexed by intervals to select a range
axisindexes{T}(::Type{Dimensional}, ax::AbstractVector{T}, idx::ClosedInterval) = searchsorted(ax, idx)

# Or repeated intervals, which only work if the axis is a range since otherwise
# there will be a non-constant number of indices in each repetition.
# Two small tricks are used here:
# * Compute the resulting interval axis with unsafe indexing without any offset
#   - Since it's a range, we can do this, and it makes the resulting axis useful
# * Snap the offsets to the nearest datapoint to avoid fencepost problems
# Adds a dimension to the result; rows represent the interval and columns are offsets.
axisindexes(::Type{Dimensional}, ax::AbstractVector, idx::RepeatedInterval) = error("repeated intervals might select a varying number of elements for non-range axes; use a repeated Range of indices instead")
function axisindexes(::Type{Dimensional}, ax::Range, idx::RepeatedInterval)
    n = length(idx.offsets)
    idxs = unsafe_searchsorted(ax, idx.window)
    offsets = [searchsortednearest(ax, idx.offsets[i]) for i=1:n]
    AxisArray(RepeatedRangeMatrix(idxs, offsets), Axis{:sub}(inbounds_getindex(ax, idxs)), Axis{:rep}(ax[offsets]))
end

# We also have special datatypes to represent intervals about indices
axisindexes(::Type{Dimensional}, ax::AbstractVector, idx::IntervalAtIndex) = searchsorted(ax, idx.window + ax[idx.index])
function axisindexes(::Type{Dimensional}, ax::Range, idx::IntervalAtIndex)
    idxs = unsafe_searchsorted(ax, idx.window)
    AxisArray(idxs + idx.index, Axis{:sub}(inbounds_getindex(ax, idxs)))
end
axisindexes(::Type{Dimensional}, ax::AbstractVector, idx::RepeatedIntervalAtIndexes) = error("repeated intervals might select a varying number of elements for non-range axes; use a repeated Range of indices instead")
function axisindexes(::Type{Dimensional}, ax::Range, idx::RepeatedIntervalAtIndexes)
    n = length(idx.indexes)
    idxs = unsafe_searchsorted(ax, idx.window)
    AxisArray(RepeatedRangeMatrix(idxs, idx.indexes), Axis{:sub}(inbounds_getindex(ax, idxs)), Axis{:rep}(ax[idx.indexes]))
end

# Categorical axes may be indexed by their elements
function axisindexes(::Type{Categorical}, ax::AbstractVector, idx)
    i = findfirst(ax, idx)
    i == 0 && throw(ArgumentError("index $idx not found"))
    i
end
# Categorical axes may be indexed by a vector of their elements
function axisindexes(::Type{Categorical}, ax::AbstractVector, idx::AbstractVector)
    res = findin(ax, idx)
    length(res) == length(idx) || throw(ArgumentError("index $(setdiff(idx,ax)) not found"))
    res
end

# This catch-all method attempts to convert any axis-specific non-standard
# indexing types to their integer or integer range equivalents using axisindexes
# It is separate from the `Base.getindex` function to allow reuse between
# set- and get- index.
@generated function to_index{T,N,D,Ax}(A::AxisArray{T,N,D,Ax}, I...)
    ex = Expr(:tuple)
    n = 0
    for i=1:length(I)
        if I[i] <: Idx
            push!(ex.args, :(I[$i]))
            n += 1
        elseif I[i] <: AbstractArray{Bool}
            push!(ex.args, :(find(I[$i])))
            n += 1
        elseif I[i] <: CartesianIndex
            for j = 1:length(I[i])
                push!(ex.args, :(I[$i][$j]))
            end
            n += length(I[i])
        elseif i <= length(Ax.parameters)
            push!(ex.args, :(axisindexes(A.axes[$i], I[$i])))
            n += 1
        else
            push!(ex.args, :(error("dimension ", $i, " does not have an axis to index")))
        end
    end
    for _=n+1:N
        push!(ex.args, :(Colon()))
    end
    meta = Expr(:meta, :inline)
    return :($meta; $ex)
end

## Extracting the full axis (name + values) from the Axis{:name} type
@inline Base.getindex{Ax<:Axis}(A::AxisArray, ::Type{Ax}) = getaxis(Ax, axes(A)...)
@inline getaxis{Ax<:Axis}(::Type{Ax}, ax::Ax, axs...) = ax
@inline getaxis{Ax<:Axis}(::Type{Ax}, ax::Axis, axs...) = getaxis(Ax, axs...)
@noinline getaxis{Ax<:Axis}(::Type{Ax}) = throw(ArgumentError("no axis of type $Ax was found"))

# Boundschecking specialization: defer to the data array.
# Note that we could unwrap AxisArrays when they are used as indices into other
# arrays within Base's to_index/to_indices methods, but that requires a bigger
# refactor to merge our to_index method with Base's.
@inline Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, A::AxisArray) = Base.checkindex(Bool, inds, A.data)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/constraints.jl
# Defines constraints supported by optimizer and how to store them.
#############################################################################

const INEQ_SETS = Union{LT, GT, ET}
const VAR_SETS = Union{LT, GT, ET, ZO, MOI.Integer}

##### Utilities for checking that JuMP model contains variables used in expression
function check_inbounds!(m::Optimizer, vi::VI)
    if !(1 <= vi.value <= m._input_problem._variable_count)
        error("Invalid variable index $vi. ($(m._input_problem._variable_count) variables in the model.)")
    end
    return nothing
end
check_inbounds!(m::Optimizer, f::SAF) = foreach(x -> check_inbounds!(m, x.variable), f.terms)
check_inbounds!(m::Optimizer, f::VECOFVAR) = foreach(x -> check_inbounds!(m, x), f.variables)
function check_inbounds!(m::Optimizer, f::SQF)
    foreach(x -> check_inbounds!(m, x.variable), f.affine_terms)
    for term in f.quadratic_terms
        check_inbounds!(m, term.variable_1)
        check_inbounds!(m, term.variable_2)
    end
    return nothing
end

MOI.supports_constraint(::Optimizer, ::Type{VI}, ::Type{S}) where S <: VAR_SETS = true
MOI.supports_constraint(::Optimizer,::Type{T},::Type{S}) where {T<:Union{SAF,SQF},S<:INEQ_SETS} = true

MOI.is_valid(m::Optimizer, v::VI) = (1 <= v.value <= m._input_problem._variable_count)
MOI.is_valid(m::Optimizer, c::CI{VI,S}) where S <: VAR_SETS = (1 <= c.value <= m._input_problem._variable_count)
MOI.is_valid(m::Optimizer, c::CI{F,S}) where {F<:Union{SAF,SQF}, S<:INEQ_SETS} = (1 <= c.value <= length(_constraints(m,F,S)))

MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{VI,S}) where S<:VAR_SETS = length(_constraints(m,VI,S))
MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F<:Union{SAF,SQF},S<:INEQ_SETS} = length(_constraints(m,F,S))

MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{VI,S}) where S<:VAR_SETS = collect(keys(_constraints(m,VI,S)))
MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{F,S}) where {F<:Union{SAF,SQF},S<:INEQ_SETS} = collect(keys(_constraints(m,F,S)))

MOI.add_variable(m::Optimizer) = VI(m._input_problem._variable_count += 1)

function MOI.add_constraint(m::Optimizer, f::F, s::S) where {F<:Union{SAF,SQF},S<:INEQ_SETS}
    check_inbounds!(m, f)
    ci = CI{F, S}(m._input_problem._constraint_count += 1)
    _constraints(m, F, S)[ci] = (f, s)
    return ci
end
function MOI.add_constraint(m::Optimizer, f::VI, s::S) where S<:VAR_SETS
    check_inbounds!(m, f)
    ci = CI{VI,S}(f.value)
    _constraints(m,VI,S)[ci] = (f, s)
    return ci
end

result_index_1_error(v::T) where T = throw(MOI.ResultIndexBoundsError{T}(v, 1))
function MOI.get(model::Optimizer, v::MOI.VariablePrimal, vi::MOI.VariableIndex)
    check_inbounds!(model, vi)
    (v.result_index != 1) && result_index_1_error(v)
    return model._global_optimizer._continuous_solution[vi.value]
end
function MOI.get(m::Optimizer{R,S,T}, v::MOI.ConstraintPrimal, c::CI{VI,<:Any}) where {R,S,T}
    (v.result_index != 1) && result_index_1_error(v)
    return MOI.get(m, MOI.VariablePrimal(), MOI.VariableIndex(c.value))
end
function MOI.get(m::Optimizer{R,Q,T}, v::MOI.ConstraintPrimal, c::CI{F,S}) where {R,Q,T,F,S}
    (v.result_index != 1) && result_index_1_error(v)
    return _constraint_primal(m._global_optimizer, F, S)[c]
end

MOI.get(m::Optimizer, ::MOI.ConstraintFunction, c::CI{F,S}) where {F,S} = _constraints(m,F,S)[c][1]
MOI.get(m::Optimizer, ::MOI.ConstraintSet, c::CI{F,S}) where {F,S} = _constraints(m,F,S)[c][2]

function MOI.empty!(m::Optimizer{R,S,T}) where {R,S,T}

    MOI.empty!(m.subsolver_block)
    MOI.empty!(m._global_optimizer)
    m._input_problem = InputProblem()
    m._working_problem = ParsedProblem()

    m._termination_status_code = MOI.OPTIMIZE_NOT_CALLED
    m._result_status_code = MOI.OTHER_RESULT_STATUS
    m._run_time = 0.0
    m._objective_value  = -Inf
    m._objective_bound  =  Inf
    m. _relative_gap     = Inf
    m._iteration_count  = 0
    m._node_count       = 0

    return nothing
end

function MOI.is_empty(m::Optimizer{R,S,T}) where {R,S,T}

    flag = true
    flag &= MOI.is_empty(m._global_optimizer)
    flag &= isempty(m._input_problem)
    flag &= isempty(m._working_problem)
    flag &= isempty(m.subsolver_block)
    flag &= m._termination_status_code == MOI.OPTIMIZE_NOT_CALLED
    flag &= m._result_status_code == MOI.OTHER_RESULT_STATUS

    # set constructor reset on empty! and  to zero in initial parse! in parse.jl
    flag &= iszero(m._run_time)
    flag &= iszero(m._iteration_count)
    flag &= iszero(m._node_count)
    flag &= m._objective_value == -Inf
    flag &= m._objective_bound ==  Inf
    flag &= m. _relative_gap   == Inf

    return flag
end

MOI.supports_incremental_interface(m::Optimizer) = true
MOI.copy_to(model::Optimizer, src::MOI.ModelLike) = MOIU.default_copy_to(model, src)

#####
##### Set & get attributes of model
#####
MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(m::Optimizer, ::MOI.Silent, value)
    if value
        m._parameters.verbosity = 0
        m._parameters.log_on = false
    else
        m._parameters.verbosity = 1
    end
    return
end

MOI.set(m::Optimizer, ::MOI.TimeLimitSec, ::Nothing) = m._parameters.time_limit = Inf
MOI.set(m::Optimizer, ::MOI.TimeLimitSec, v::Float64) = m._parameters.time_limit = v

function MOI.get(m::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    constraint_types = []
    for S in (ZO, MOI.Integer)
        if MOI.get(m, MOI.NumberOfConstraints{VI,S}()) > 0
            push!(constraint_types, (VI,S))
        end
    end
    for S in (LT, GT, ET), T in (VI, SAF, SQF)
        if MOI.get(m, MOI.NumberOfConstraints{T,S}()) > 0
            push!(constraint_types, (VI,S))
        end
    end
    return constraint_types
end

MOI.get(m::Optimizer, v::MOI.ObjectiveValue) = !isone(v.result_index) ? result_index_1_error(v) : m._objective_value
MOI.get(m::Optimizer, v::MOI.PrimalStatus) = !isone(v.result_index) ? MOI.NO_SOLUTION : m._result_status_code
MOI.get(m::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION
MOI.get(m::Optimizer, ::MOI.ObjectiveBound) = m._objective_bound
MOI.get(m::Optimizer, ::MOI.NumberOfVariables) = m._input_problem._variable_count
MOI.get(m::Optimizer, ::MOI.SolverName) = "EAGO: Easy Advanced Global Optimization"
MOI.get(m::Optimizer, ::MOI.SolverVersion) = "0.7.0"
MOI.get(m::Optimizer, ::MOI.TerminationStatus) = m._termination_status_code
MOI.get(m::Optimizer, ::MOI.SolveTimeSec) = m._run_time
MOI.get(m::Optimizer, ::MOI.NodeCount) = m._node_count
MOI.get(m::Optimizer, ::MOI.ResultCount) = (m._result_status_code === MOI.FEASIBLE_POINT) ? 1 : 0
MOI.get(m::Optimizer, ::MOI.TimeLimitSec) = m._parameters.time_limit
MOI.get(m::Optimizer, ::MOI.Silent) = m._parameters.verbosity == 0
MOI.get(m::Optimizer, ::MOI.ListOfVariableIndices) = [VI(i) for i = 1:m._input_problem._variable_count]

function MOI.get(m::Optimizer, ::MOI.RelativeGap)
    b = MOI.get(m, MOI.ObjectiveBound())
    v = MOI.get(m, MOI.ObjectiveValue())
    return relative_gap(b,v)
end

_to_sym(d) = error("EAGO only supports raw parameters with Symbol or String names.")
_to_sym(d::String) = Symbol(d)
_to_sym(d::Symbol) = d
function MOI.get(m::Optimizer, p::MOI.RawOptimizerAttribute)
    s = _to_sym(p.name)
    s in EAGO_PARAMETERS ? getfield(m._parameters, s) : getfield(m, s)
end
function MOI.set(m::Optimizer, p::MOI.RawOptimizerAttribute , x)
    s = _to_sym(p.name)
    if (s == :relaxed_optimizer) || (s == :upper_optimizer)
        setfield!(m, s, Incremental(x))
    else
        s in EAGO_PARAMETERS ? setfield!(m._parameters, s, x) : setfield!(m, s, x)
    end
end
MOI.get(m::Optimizer, p::MOI.RawStatusString) = string(m._global_optimizer._end_state)

#####
##### Support, set, and evaluate objective functions
#####
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F <: Union{VI, SAF, SQF}} = true

function MOI.set(m::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    if nlp_data.has_objective
        m._input_problem._objective = nothing
    end
    m._input_problem._nlp_data = nlp_data
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{T}, f::T) where T <: Union{VI,SAF,SQF}
    check_inbounds!(m, f)
    m._input_problem._objective = f
end
MOI.get(m::Optimizer, ::MOI.ObjectiveFunction{T}) where T <: Union{VI,SAF,SQF} = m._input_problem._objective
MOI.get(m::Optimizer, ::MOI.ObjectiveFunctionType) = typeof(m._input_problem._objective)

MOI.set(m::Optimizer, ::MOI.ObjectiveSense, s::MOI.OptimizationSense) = m._input_problem._optimization_sense = s
MOI.get(m::Optimizer, ::MOI.ObjectiveSense) = m._input_problem._optimization_sense

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions