text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "GoogleSparseHash"
version = v"2.11.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sparsehash/sparsehash-c11.git", "edd6f1180156e76facc1c0449da245208ab39503")
]

script = raw"""
cd sparsehash-*
install_license LICENSE

mkdir ${prefix}/include
mv sparsehash ${prefix}/include/
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

products = [
    FileProduct("include/sparsehash/dense_hash_map", :dense_hash_map),
    FileProduct("include/sparsehash/dense_hash_set", :dense_hash_set),
    FileProduct("include/sparsehash/sparse_hash_map", :sparse_hash_map),
    FileProduct("include/sparsehash/sparse_hash_set", :sparse_hash_set),
    FileProduct("include/sparsehash/sparsetable", :sparsetable),
    FileProduct("include/sparsehash/traits", :traits),
    FileProduct("include/sparsehash/internal/densehashtable.h", :densehashtable_h),
    FileProduct("include/sparsehash/internal/hashtable-common.h", :hashtable_common_h),
    FileProduct("include/sparsehash/internal/libc_allocator_with_realloc.h", :libc_allocator_with_realloc_h),
    FileProduct("include/sparsehash/internal/sparsehashtable.h", :sparsehashtable_h)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version = v"6.1.0")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
function can_form_government(model::Model)
    party_combs = all_party_combinations(model)
    seats = [combined_seats(model, comb) for comb in party_combs]
    party_consensus_dict = Dict(
        party => party_consensus(model, party)
        for party in model.parameter_set.parties
    )
    pairwise_similarities = [
        pair_similarities(pc, party_consensus_dict)
        for pc in party_combs
    ]
    have_consensus_list = [
        have_consensus(model, s)
        for s in pairwise_similarities
    ]
    return sum(have_consensus_list) == 1
end

function have_consensus(model::Model, similarities::AbstractArray)
    entrywise_have_consensus = [
        s .> model.parameter_set.required_consensus
        for s in similarities
    ]
    return sum(entrywise_have_consensus) == length(entrywise_have_consensus)
end

function all_party_combinations(model::Model)
    return collect(Combinatorics.combinations(model.parameter_set.parties))
end

function combined_seats(model::Model, parties::AbstractArray)
     return sum([model.parameter_set.parliament[party] for party in parties])
end

function pair_similarities(
    parties::AbstractArray,
    party_consensus_dict::AbstractDict
)
    all_pairs = Combinatorics.combinations(parties, 2)
    similarities = [
        similarity(
            [party_consensus_dict[pair[1]], party_consensus_dict[pair[2]]]
        )
        for pair in all_pairs
    ]
    return similarities
end

# TODO: refactor
#       setup not ideal because it assumes an array of length 2
function similarity(opinions::AbstractArray)
    entrywise_difference = sum(abs.(opinions[1] .- opinions[2]))
    maximum_difference = 2 * length(opinions[1])
    return 1 - (entrywise_difference / maximum_difference)
end

function party_consensus(model, party)
    agents = filter(a -> a.party == party, model.agents)
    opinions = [a.opinions for a in agents]
    consensus = [
        StatsBase.mode([o[i] for o in opinions])
        for i in 1:length(model.opinions[Symbol(party)])
    ]
    return consensus
end

ext: jl


text:
# pfile.jl

using Test: @test
using MIRTio: header_write
using MIRTio: loadpfile, rdb_hdr_26_002_def
using MIRTio: header_init, header_size, header_string

hd = rdb_hdr_26_002_def()
ht = header_init(hd) # random header values
hs = header_size(hd)

necho = 2
nslice = 7
nview = 64
ndat = 100
ncoil = 5
ht = merge(ht, [:rdbm_rev => Float32(26.002)]) # instead of setindex
ht = merge(ht, [:frame_size => UInt16(ndat)])
ht = merge(ht, [:nslices => UInt16(nslice)])
ht = merge(ht, [:point_size => Int16(2)])
ht = merge(ht, [:nechoes => Int16(necho)])
ht = merge(ht, [:nframes => Int16(nview)])
dab = zeros(Int16,8)
dab[2] = Int16(ncoil-1)
#@show dab
ht = merge(ht, [:dab => dab])
ht = merge(ht, [:off_data => Int32(hs)])

tname = tempname()
open(tname, "w") do fid
    header_write(fid, ht)
end

# Apparent file data order is: [ndat, nview, necho, nslice, ncoil]
fdat = rand(Complex{Int16}, ndat, nview+1, necho, nslice, ncoil) # rand test data
open(tname, "a") do fid
    write(fid, fdat)
end

(odat, hr) = loadpfile(tname ; quiet=true) # load and verify
@test hr == header_string(ht)
pdat = permutedims(odat, [1, 5, 4, 3, 2])
#=
@show size(fdat)
@show size(odat)
@show size(pdat)
@test pdat == fdat[:, 2:end, :, 2:end, :] # skip 1st view and slice
=#

iecho = 2
(dat1, _) = loadpfile(tname, iecho ; quiet=true)
pdat = permutedims(dat1, [1, 5, 4, 3, 2])
gdat = fdat[:, 2:end, [iecho], 2:end, :] # skip 1st view and slice
#=
@show size(fdat)
@show size(dat1)
@show size(pdat)
@show size(gdat)
=#
@test pdat == gdat

nview = 2 # test special case
ht = merge(ht, [:nframes => Int16(nview)])
fdat = zeros(Complex{Int16}, ndat, nview+1, necho, nslice, ncoil)
open(tname, "w") do fid
    header_write(fid, ht)
    write(fid, fdat)
end
loadpfile(tname ; quiet=true)


text:
module IndexNow

using Franklin; const fr = Franklin

const apikey = Ref("")
const search_engines = (
    "bing.com",
    "yandex.com",
)

@doc "Set the indexnow apikey, if KEY is empty use the environment variable INDEXNOW_APIKEY."
function set_apikey!(key="")
	apikey[] = isempty(key) ? get(ENV, "INDEXNOW_APIKEY", "") : key
    try
	    k = read(joinpath(:folder ∈ fr.PATHS ? fr.path(:folder) : pwd(), "$(apikey).txt"), String)
        @assert apikey[] === k
    catch
        throw("Couldn't find the api key text file, ensure it is present in website root directory.")
    end
end

function push_urls(urls=[])
    for u in urls

    end
end

end


text:
struct PlsLda
    fm
    lev::AbstractVector
    ni::AbstractVector
end

"""
    plslda(X, y, weights = ones(size(X, 1)); nlv)
LDA on PLS latent variables.
* `X` : X-data.
* `y` : y-data (class membership).
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to compute.

The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a PLS2 is implemented on the X-data and Ydummy,
returning `nlv` latent variables (LVs).
Finally, a LDA is run on these LVs and `y`.
"""
function plslda(X, y, weights = ones(size(X, 1)); nlv, prior = "unif")
    z = dummy(y)
    fm_pls = plskern(X, z.Y, weights; nlv = nlv)
    fm_da = list(nlv)
    @inbounds for i = 1:nlv
        fm_da[i] = lda(fm_pls.T[:, 1:i], y; prior = prior)
    end
    fm = (fm_pls = fm_pls, fm_da = fm_da)
    PlsLda(fm, z.lev, z.ni)
end

function predict(object::PlsLda, X; nlv = nothing)
    X = ensure_mat(X)
    m = size(X, 1)
    a = size(object.fm.fm_pls.T, 2)
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    le_nlv = length(nlv)
    pred = list(le_nlv, Matrix{Float64})
    posterior = list(le_nlv, Matrix{Float64})
    @inbounds for i = 1:le_nlv
        znlv = nlv[i]
        T = transform(object.fm.fm_pls, X, nlv = znlv)
        zres = predict(object.fm.fm_da[znlv], T)
        z =  mapslices(argmax, zres.posterior; dims = 2)
        pred[i] = reshape(replacebylev(z, object.lev), m, 1)
        posterior[i] = zres.posterior
    end

    if le_nlv == 1
        pred = pred[1]
        posterior = posterior[1]
    end
    (pred = pred, posterior = posterior)
end

"""
    plsqda(X, y, weights = ones(size(X, 1)); nlv)
QDA on PLS latent variables.
* `X` : X-data.
* `y` : y-data (class membership).
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to compute.

The training variable y (univariate class membership) is transformed
to a dummy table (Ydummy) containing nlev columns, where nlev is the number
of classes present in `y`. Each column of Ydummy is a dummy variable (0/1).
Then, a PLS2 is implemented on the X-data and Ydummy, returning `nlv` latent
variables (LVs). Finally, a QDA is run on these LVs and `y`.
"""
function plsqda(X, y, weights = ones(size(X, 1)); nlv, prior = "unif")
    z = dummy(y)
    fm_pls = plskern(X, z.Y, weights; nlv = nlv)
    fm_da = list(nlv)
    for i = 1:nlv
        fm_da[i] = qda(vcol(fm_pls.T, 1:i), y; prior = prior)
    end
    fm = (fm_pls = fm_pls, fm_da = fm_da)
    PlsLda(fm, z.lev, z.ni)
end


text:
export Goal, RectangleGoal, BallGoal, ConvexHullWorkspaceGoal, ConvexHullStateSpaceGoal
export PointGoal, StateGoal, is_goal_pt, sample_goal

abstract type Goal end
abstract type WorkspaceGoal <: Goal end
abstract type StateSpaceGoal <: Goal end

struct RectangleGoal{N,T<:AbstractFloat} <: WorkspaceGoal
    lo::SVector{N,T}
    hi::SVector{N,T}
end
RectangleGoal(lo::AbstractVector, hi::AbstractVector) = RectangleGoal(SVector(lo), SVector(hi))
changeprecision{T<:AbstractFloat}(::Type{T}, G::RectangleGoal) = RectangleGoal(changeprecision(T, G.lo),
                                                                               changeprecision(T, G.hi))

struct BallGoal{N,T<:AbstractFloat} <: WorkspaceGoal
    center::SVector{N,T}
    radius::T
end
BallGoal(center::AbstractVector, radius) = BallGoal(SVector(center), radius)
changeprecision{T<:AbstractFloat}(::Type{T}, G::BallGoal) = BallGoal(changeprecision(T, G.center), T(radius))

struct ConvexHullWorkspaceGoal{N,T<:AbstractFloat} <: WorkspaceGoal
    pts::Vector{SVector{N,T}}
    A::Matrix{Float64}    # for use with the SCS LP solver
    b::Vector{Float64}    # for use with the SCS LP solver
    c::Vector{Float64}    # for use with the SCS LP solver
    m::Int
    n::Int

    function ConvexHullWorkspaceGoal{N,T}(pts::Vector{SVector{N,T}}) where {N,T}
        m, n = N, length(pts)
        A = n > 1 ? convert(Matrix{Float64}, [statevec2mat(pts); ones(1,n); -eye(n)]) : Array{Float64}(0,0)
        b = zeros(m + n + 1); b[m+1] = 1
        c = zeros(n)
        new(pts, A, b, c, m, n)
    end
end
ConvexHullWorkspaceGoal{N,T<:AbstractFloat}(pts::Vector{SVector{N,T}}) = ConvexHullWorkspaceGoal{N,T}(pts)
ConvexHullWorkspaceGoal{V<:AbstractVector}(pts::Vector{V}) = ConvexHullWorkspaceGoal([SVector(p) for p in pts])
PointGoal(pt::AbstractVector) = ConvexHullWorkspaceGoal([pt])
changeprecision{T<:AbstractFloat}(::Type{T}, G::ConvexHullWorkspaceGoal) =
    ConvexHullWorkspaceGoal(changeprecision(T, G.pts))

struct ConvexHullStateSpaceGoal{S<:State} <: StateSpaceGoal
    sts::Vector{S}
    A::Matrix{Float64}    # for use with the SCS LP solver
    b::Vector{Float64}    # for use with the SCS LP solver
    c::Vector{Float64}    # for use with the SCS LP solver
    m::Int
    n::Int

    function ConvexHullStateSpaceGoal{S}(sts::Vector{S}) where S
        m, n = length(sts[1]), length(sts)
        A = n > 1 ? convert(Matrix{Float64}, [statevec2mat(sts); ones(1,n); -eye(n)]) : Array{Float64}(0,0)
        b = zeros(m + n + 1); b[m+1] = 1
        c = zeros(n)
        new(sts, A, b, c, m, n)
    end
end
ConvexHullStateSpaceGoal{S<:State}(sts::Vector{S}) = ConvexHullStateSpaceGoal{S}(sts)
ConvexHullStateSpaceGoal{S<:SVector}(sts::Vector{S}) = ConvexHullStateSpaceGoal{S}(sts)
ConvexHullStateSpaceGoal{S<:FieldVector}(sts::Vector{S}) = ConvexHullStateSpaceGoal{S}(sts)
ConvexHullStateSpaceGoal{V<:AbstractVector}(sts::Vector{V}) = ConvexHullStateSpaceGoal([SVector(s) for s in sts])
StateGoal(s::State) = ConvexHullStateSpaceGoal([s])
changeprecision{T<:AbstractFloat}(::Type{T}, G::ConvexHullStateSpaceGoal) =
    ConvexHullStateSpaceGoal(changeprecision(T, G.sts))

plot(G::RectangleGoal, SS::StateSpace) = plot_rectangle(G.lo, G.hi, color = "green")
plot(G::BallGoal, SS::StateSpace) = plot_circle(G.center, G.radius, color = "green")
function plot(G::ConvexHullWorkspaceGoal, SS::StateSpace)
    if G.n == 1    # PointGoal
        plt.scatter(G.pts[1][1], G.pts[1][2], color = "green", zorder=5)
    elseif G.n == 2
        plot_line_segments([G.pts[1]], [G.pts[2]], color = "green", zorder=5, linewidth=2.5)
    else
        plot_polygon(G.pts, color = "green", zorder=5)
    end
end
function plot(G::ConvexHullStateSpaceGoal, SS::StateSpace)
    pts = [state2workspace(s, SS) for s in G.sts]
    if G.n == 1    # StateGoal
        plt.scatter(pts[1][1], pts[1][2], color = "green", zorder=5)
    elseif G.n == 2
        plot_line_segments([pts[1]], [pts[2]], color = "green", zorder=5, linewidth=2.5)
    else
        plot_polygon(pts, color = "green", zorder=5)
    end
end

sample_goal(G::WorkspaceGoal, SS::StateSpace) = workspace2state(sample_goal(G), SS)
sample_goal(G::StateSpaceGoal, SS::StateSpace) = sample_goal(G)

## Rectangle Goal
is_goal_pt(v::State, G::RectangleGoal, SS::StateSpace) = (reduce(&, G.lo .<= state2workspace(v, SS) .<= G.hi) != 0)
sample_goal(G::RectangleGoal) = G.lo + (G.hi - G.lo).*rand(typeof(G.lo))

## Ball Goal
is_goal_pt(v::State, G::BallGoal, SS::StateSpace) = (norm(state2workspace(v, SS) - G.center) <= G.radius)
function sample_goal(G::BallGoal)
    while true
        v = G.center + 2*G.radius*(rand(typeof(G.center)) - .5)
        if norm(v - G.center) <= G.radius
            return v
        end
    end
end

## Convex Hull Goals
function is_goal_pt(v::State, G::ConvexHullWorkspaceGoal, SS::StateSpace)
    pt = state2workspace(v, SS)
    if G.n == 1
        pt == G.pts[1]
    elseif G.n == 2
        isapprox(norm(G.pts[1] - G.pts[2]), norm(G.pts[1] - pt) + norm(G.pts[2] - pt))    # not the most efficient way
    else
        G.b[1:G.m] = pt
        sol = SCS_solve(G.m+G.n+1, G.n, G.A, G.b, G.c, G.m+1, G.n, Int[], 0, Int[], 0, 0, 0, Float64[], 0, verbose=0)
        sol.status != :Infeasible
    end
end
function sample_goal{N,T}(G::ConvexHullWorkspaceGoal{N,T})
    G.n == 1 ? G.pts[1] : G.n == 2 ? G.pts[1] + rand(T)*(G.pts[2] - G.pts[1]) : sum(G.pts .* rand(Dirichlet(G.n, T(1))))
end

function is_goal_pt(v::State, G::ConvexHullStateSpaceGoal, SS::StateSpace)
    if G.n == 1
        v == G.sts[1]
    elseif G.n == 2
        isapprox(norm(G.sts[1] - G.sts[2]), norm(G.sts[1] - v) + norm(G.sts[2] - v))    # not the most efficient way
    else
        G.b[1:G.m] = v
        sol = SCS_solve(G.m+G.n+1, G.n, G.A, G.b, G.c, G.m+1, G.n, Int[], 0, Int[], 0, 0, 0, Float64[], 0, verbose=0)
        sol.status != :Infeasible
    end
end
function sample_goal{S}(G::ConvexHullStateSpaceGoal{S})
    T = eltype(S)
    G.n == 1 ? G.sts[1] : G.n == 2 ? G.sts[1] + rand(T)*(G.sts[2] - G.sts[1]) : sum(G.sts .* rand(Dirichlet(G.n, T(1))))
end


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

mutable struct ObjectReplicationPolicyProperties <: SwaggerModel
    policyId::Any # spec type: Union{ Nothing, String } # spec name: policyId
    enabledTime::Any # spec type: Union{ Nothing, DateTime } # spec name: enabledTime
    sourceAccount::Any # spec type: Union{ Nothing, String } # spec name: sourceAccount
    destinationAccount::Any # spec type: Union{ Nothing, String } # spec name: destinationAccount
    rules::Any # spec type: Union{ Nothing, Vector{ObjectReplicationPolicyRule} } # spec name: rules

    function ObjectReplicationPolicyProperties(;policyId=nothing, enabledTime=nothing, sourceAccount=nothing, destinationAccount=nothing, rules=nothing)
        o = new()
        validate_property(ObjectReplicationPolicyProperties, Symbol("policyId"), policyId)
        setfield!(o, Symbol("policyId"), policyId)
        validate_property(ObjectReplicationPolicyProperties, Symbol("enabledTime"), enabledTime)
        setfield!(o, Symbol("enabledTime"), enabledTime)
        validate_property(ObjectReplicationPolicyProperties, Symbol("sourceAccount"), sourceAccount)
        setfield!(o, Symbol("sourceAccount"), sourceAccount)
        validate_property(ObjectReplicationPolicyProperties, Symbol("destinationAccount"), destinationAccount)
        setfield!(o, Symbol("destinationAccount"), destinationAccount)
        validate_property(ObjectReplicationPolicyProperties, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type ObjectReplicationPolicyProperties

const _property_map_ObjectReplicationPolicyProperties = Dict{Symbol,Symbol}(Symbol("policyId")=>Symbol("policyId"), Symbol("enabledTime")=>Symbol("enabledTime"), Symbol("sourceAccount")=>Symbol("sourceAccount"), Symbol("destinationAccount")=>Symbol("destinationAccount"), Symbol("rules")=>Symbol("rules"))
const _property_types_ObjectReplicationPolicyProperties = Dict{Symbol,String}(Symbol("policyId")=>"String", Symbol("enabledTime")=>"DateTime", Symbol("sourceAccount")=>"String", Symbol("destinationAccount")=>"String", Symbol("rules")=>"Vector{ObjectReplicationPolicyRule}")
Base.propertynames(::Type{ ObjectReplicationPolicyProperties }) = collect(keys(_property_map_ObjectReplicationPolicyProperties))
Swagger.property_type(::Type{ ObjectReplicationPolicyProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjectReplicationPolicyProperties[name]))}
Swagger.field_name(::Type{ ObjectReplicationPolicyProperties }, property_name::Symbol) =  _property_map_ObjectReplicationPolicyProperties[property_name]

function check_required(o::ObjectReplicationPolicyProperties)
    (getproperty(o, Symbol("sourceAccount")) === nothing) && (return false)
    (getproperty(o, Symbol("destinationAccount")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjectReplicationPolicyProperties }, name::Symbol, val)
end


text:
using Distributions

include("ADVI.jl")

# Dirichlet
d = Dirichlet(collect(1:3)); x = collect(1:3) / sum(1:3)
@assert ADVI.compute_lpdf(d, x) == logpdf(d, x)

# Normal
d = Normal(2, 3); x = -5.0
@assert ADVI.compute_lpdf(d, x) == logpdf(d, x)

# Beta
d = Beta(2, 3); x = .6
@assert ADVI.compute_lpdf(d, x) == logpdf(d, x)

# Gamma
d = Gamma(2,4); x = 5.0
@assert ADVI.compute_lpdf(d, x) == logpdf(d, x)

# LogNormal
d = LogNormal(2, 4); x = 1.5
@assert ADVI.compute_lpdf(d, x) == logpdf(d, x)
# ADVI.compute_lpdf(d, x)
# logpdf(d, x)

# Uniform
d = Uniform(2,4); x = 3.5
@assert ADVI.compute_lpdf(d, x) == logpdf(d, x)


text:
# # Getting Started with Power Graphics

# This uses the plotting package for SIIP [PowerGraphics](http://github.com/nrel-siip/PowerGraphics.jl.git)

# Start by creating a results object using [PowerSimulations](http://github.com/nrel-siip/PowerSimulations.jl.git)

# Set up Power Graphics:
# ```julia
# import PowerGraphics
# const PG = PowerGraphics
# ```

# ## To make a simple GR() backend plot (static plot), simply call

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

# ![this one](plots-01/P__ThermalStandard_Stack.png)
# ![this one](plots-01/Example_GR_Plot.png)

# ## To make an interactive PlotlyJS plot, reset the backend

# ```julia
# Plots.plotlyjs()

# stack_plot(results)
# ```
# ![this one](plots-02/P__ThermalStandard_Stack.png)
# ![this one](plots-02/Example_PlotlyJS_Plot.png)


text:
export fourier_diff

"""
    fourier_diff([T=Float64,] N; order=1)

Create a Fourier differentiation matrix with numerical type T on the domain
`x = range(0, 2π, length=N+1)[1:end-1]`.
"""
function fourier_diff(T::Type{<:Number}, N::Integer; order=1)
    D = zeros(T, N, N)
    n1 = (N - 1) ÷ 2
    n2 = N ÷ 2
    x = LinRange{T}(0, π, N+1)
    if order == 1
        for i in 2:N
            sgn = (one(T)/2 - iseven(i))
            D[i, 1] = iseven(N) ? sgn*cot(x[i]) : sgn*csc(x[i])
        end
    elseif order == 2
        D[1, 1] = iseven(N) ? -N^2*one(T)/12 - one(T)/6 : -N^2*one(T)/12 + one(T)/12
        for i in 2:N
            sgn = -(one(T)/2 - iseven(i))
            D[i, 1] = iseven(N) ? sgn*csc(x[i]).^2 : sgn*cot(x[i])*csc(x[i])
        end
    else
        error("Not implemented")
    end
    for j in 2:N
        D[1, j] = D[N, j-1]
        D[2:N, j] .= D[1:N-1, j-1]
    end
    return D
end
fourier_diff(N::Integer; kwargs...) = fourier_diff(Float64, N; kwargs...)


text:
"""
	loop constraints
"""

struct LoopConstraints <: Constraints
	n
	ineq
	idx_1
	idx_2
end

function loop_constraints(model, idx1, idx2)
	n = model.n
	ineq = (1:0)
	return LoopConstraints(n, ineq, idx1, idx2)
end

function constraints!(c, Z, con::LoopConstraints, model, idx, h, T)
	n = model.n
	a = view(Z, idx.x[con.idx_1])
	b = view(Z, idx.x[con.idx_2])

	c[1:n] = b - a

    return nothing
end

function constraints_jacobian!(∇c, Z, con::LoopConstraints, model, idx, h, T)
	n = model.n
	shift = 0

	r_idx = 1:n
	c_idx = idx.x[con.idx_2]
	len = length(r_idx) * length(c_idx)
	∇c[shift .+ (1:len)] = vec(Diagonal(ones(n)))
	shift += len

	c_idx = idx.x[con.idx_1]
	len = length(r_idx) * length(c_idx)
	∇c[shift .+ (1:len)] = vec(Diagonal(-1.0 * ones(n)))
	shift += len
	return nothing
end

function constraints_sparsity(con::LoopConstraints, model, idx, T;
		shift_row = 0, shift_col = 0)
	row = []
    col = []

	n = model.n

	r_idx = shift_row .+ (1:n)

	c_idx = shift_col .+ idx.x[con.idx_2]
	row_col!(row, col, r_idx, c_idx)

	c_idx = shift_col .+ idx.x[con.idx_1]
	row_col!(row, col, r_idx, c_idx)

    return collect(zip(row, col))
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior.Gherkin: parsescenario!, issuccessful, Given, When, Then, ByLineParser, ScenarioStep

@testset "Scenario Outline     " begin
    @testset "Outline has a Given step; Step is parsed" begin
        text = """
        Scenario Outline: This is one scenario outline
            Given a precondition with field <Foo>

        Examples:
            | Foo |
            | 1   |
            | 2   |
        """
        byline = ByLineParser(text)

        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value
        @test scenario.steps == ScenarioStep[Given("a precondition with field <Foo>")]
    end

    @testset "Scenario Outline has description; Description is parsed" begin
        text = """
        Scenario Outline: This is one scenario outline
            Given a precondition with field <Foo>

        Examples:
            | Foo |
            | 1   |
            | 2   |
        """
        byline = ByLineParser(text)

        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value
        @test scenario.description == "This is one scenario outline"
    end

    @testset "Scenario Outline has tags; Tags are parsed" begin
        text = """
        @tag1 @tag2
        Scenario Outline: This is one scenario outline
            Given a precondition with field <Foo>

        Examples:
            | Foo |
            | 1   |
            | 2   |
        """
        byline = ByLineParser(text)

        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value
        @test scenario.tags == ["@tag1", "@tag2"]
    end

    @testset "Scenario Outline Examples" begin
        @testset "Outline has three placeholders; The placeholders are parsed" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 1   | 2   | 3   |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.placeholders == ["Foo", "Bar", "Baz"]
        end

        @testset "Two examples with three placeholders are provided; Examples array is 3x2" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 1   | 2   | 3   |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test size(scenario.examples) == (3,2)
        end

        @testset "Three examples with four placeholders are provided; Examples array is 4x3" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>, <Quux>

            Examples:
                | Foo | Bar | Baz | Quux |
                | 1   | 2   | 3   | 4    |
                | 1   | 2   | 3   | 4    |
                | 1   | 2   | 3   | 4    |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test size(scenario.examples) == (4,3)
        end

        @testset "Two examples with three placeholders are provided; Examples array has all values" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar>, <Baz>

            Examples:
                | Foo | Bar | Baz |
                | 1   | 2   | 3   |
                | 4   | 5   | 6   |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.examples[:,1] == ["1", "2", "3"]
            @test scenario.examples[:,2] == ["4", "5", "6"]
        end

        @testset "Examples with spaces; Examples are split on | not on spaces" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholder <Foo>

            Examples:
                | Foo       |
                | word      |
                | two words |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.examples[:,1] == ["word"]
            @test scenario.examples[:,2] == ["two words"]
        end

        @testset "Example with an empty element; The empty line results in an empty value" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholder <Foo>

            Examples:
                | Foo       |
                ||
                | two words |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.examples[:,1] == [""]
            @test scenario.examples[:,2] == ["two words"]
        end

        @testset "A placeholder is non-alphanumerical; The placeholders are correct" begin
            text = """
            Scenario Outline: This is one scenario outline
                Given a precondition with placeholders <Foo>, <Bar-baz>, <>

            Examples:
                | Foo | Bar-baz |
                | 1   | 2       |
                | 4   | 5       |
            """
            byline = ByLineParser(text)

            result = parsescenario!(byline)

            @test issuccessful(result)
            scenario = result.value
            @test scenario.placeholders == ["Foo", "Bar-baz"]
        end
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2021/14")
end
using Pkg
Pkg.activate(".")
using StatsBase: countmap
using Base.Iterators

function loadfile(filename::AbstractString)
    file = readchomp(filename)
    target, ps = split(file, r"\r?\n\r?\n")
    rules = Dict(a => first(b) for (a, b) in split.(split(ps, r"\r?\n"), " -> "))
    target, rules
end

"Returns an iterator of each overlapping pair of characters in a string."
charpairs(str::AbstractString) = (view(str, i-1:i) for i in drop(eachindex(str), 1))

function part1(filename::AbstractString)
    target, rules = loadfile(filename)
    for _ = 1:10
        result = ""
        for (a, b) in charpairs(target)
            result *= a
            result *= rules[a*b]
        end
        result *= target[end]
        target = result
    end
    cm = values(countmap(target))
    maximum(cm) - minimum(cm)
end

function part2(filename::AbstractString, n::Integer)
    target, oldrules = loadfile(filename)
    rules = Dict{String,Tuple{String,String}}(
        k => (first(k) * v, v * last(k)) for (k, v) in pairs(oldrules)
    )
    counts = Dict{String,Int}()
    for cs in charpairs(target)
        counts[cs] = get(counts, cs, 0) + 1
    end
    for _ = 1:n
        oldcounts = counts
        counts = empty(counts)
        for (k, v) in pairs(oldcounts)
            if k in keys(rules) && !iszero(v)
                for o in rules[k]
                    counts[o] = get(counts, o, 0) + v
                end
            end
        end
    end
    charcounts = Dict{Char,Int}((
        first(target) => 1,
        last(target) => 1
    ))
    for (lr, v) in pairs(counts)
        l, r = lr
        charcounts[l] = get(charcounts, l, 0) + v
        charcounts[r] = get(charcounts, r, 0) + v
    end
    map!(v -> v ÷ 2, values(charcounts))
    maximum(values(charcounts)) - minimum(values(charcounts))
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
function ensure_domain(m::AbstractMetricParams{T}) where {T}
    min_r = inner_radius(m)
    (u, λ, integrator) -> u[2] ≤ min_r * 1.1 || u[2] > 1200.0
end

function flip_radial_sign!(integrator)
    p = integrator.p
    integrator.p = @set p.r = -p.r
    integrator.sol.prob.p.changes[1] = integrator.t[end]
end

function flip_angular_sign!(integrator)
    p = integrator.p
    integrator.p = @set p.θ = -p.θ
    integrator.sol.prob.p.changes[2] = integrator.t[end]
end

function radial_negative_check(m::CarterMethodBL{T}) where {T}
    (u, λ, integrator) -> begin
        Vr(m, u, integrator.p) < 0
    end
end

function angular_negative_check(m::CarterMethodBL{T}) where {T}
    (u, λ, integrator) -> begin
        Vθ(m, u, integrator.p) < 0
    end
end


text:
using Pkg
Pkg.activate(pwd())

# # Conditional evalutions
# ## `if-elseif-else` statement

function compare(x, y)
    if x < y
        println("x is less than y")
    elseif x > y
        println("x is greater than y")
    else
        println("x is equal to y")
    end
    return
end

#+

compare(1, 2.3)
compare(4.7, 2.3)
compare(2.3, 2.3)

#+

x, y = 2, 1;

if x < y
    println("x is less than y")
elseif x > y
    println("x is greater than y")
end

#+

if x < y
    println("x is less than y")
end

#+

if 1
    println("Hello")
end

#+

x, y = 2, 1;

if x < y
    z = y
else
    z = x
end

z

#+

function compare(x, y)
    if x < y
        z = y
    elseif x > y
        z = x
    end
    return z
end

#+

compare(1, 2.3)
compare(4.7, 2.3)
compare(2.3, 2.3)

#+

function compare(x, y)
    z = if x < y
        y
    else
        x
    end
    return z
end

#+

compare(1, 2.3)
compare(4.7, 2.3)
compare(2.3, 2.3)

# ### Exercise:
# Write the `fact(n)` function that computes the factorial of `n`. Use the following
# function declaration:

function fact(n)
    # some code
end

# Make sure that the input argument is a non-negative integer. For negative input arguments
# and for arguments that can not be represented as an integer, the function should throw an
# error.
#
# **Hint:** use recursion, the `isinteger` function and the `error` function.
# The or operator
# is written by `|`.
#
# ---
# ### Solution:

function fact(n)
    return if n < 0 | !isinteger(n)
        error("argument must be non-negative integer")
    elseif n == 0
        1
    else
        n * fact(n - 1)
    end
end

#+

fact(4)
fact(0)
fact(-5)
fact(1.4)

# ---
#
# ## Ternary operator

x, y = 2, 1;
println(x < y ? "x is less than y" : "x is greater than or equal to y")

# ## Short-circuit evaluation

t(x) = (println(x); true)
f(x) = (println(x); false)

t(1) && println(2) # both expressions are evaluated
f(1) && println(2) # only the first expression is evaluated
t(1) || println(2) # only the first expression is evaluated
f(1) || println(2) # both expressions are evaluated

#+

f(1) & t(2)
f(1) && t(2)

#+

t(1) && t(2) || println(3) # the first two expressions are evaluated
f(1) && t(2) || println(3) # the first and the last expressions are evaluated
f(1) && f(2) || println(3) # the first and the last expressions are evaluated
t(1) && f(2) || println(3) # all expressions are evaluated

#+

t(1) || t(2) && println(3) # the first expression is evaluated
f(1) || t(2) && println(3) # all expressions are evaluated
f(1) || f(2) && println(3) # the first two expressions are evaluated
t(1) || f(2) && println(3) # the first expression is evaluated

# ### Exercise:
# Rewrite the factorial function from the exercises above. Use the short-circuit evaluation
# to check if the given number is a non-negative integer and the ternary operator for
# recursion.
#
# ---
# ### Solution:

function fact(n)
    isinteger(n) || error("argument must be non-negative integer")
    n >= 0 || error("argument must be non-negative integer")
    return n == 0 ? 1 : n * fact(n - 1)
end

#+

function fact(n)
    isinteger(n) && n >= 0 || error("argument must be non-negative integer")
    return n == 0 ? 1 : n * fact(n - 1)
end

#+

fact(4)
fact(0)
fact(-5)
fact(1.4)


text:
using Downloads: download
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

equations = IdealGlmMhdEquations2D(5/3)

function initial_condition_shifted_weak_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
  # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
  # Shift blastwave to center of domain
  inicenter = (sqrt(2)/2, sqrt(2)/2)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)
  phi = atan(y_norm, x_norm)

  # Calculate primitive variables
  rho = r > 0.5 ? 1.0 : 1.1691
  v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi)
  v2 = r > 0.5 ? 0.0 : 0.1882 * sin(phi)
  p = r > 0.5 ? 1.0 : 1.245

  return prim2cons(SVector(rho, v1, v2, 0.0, p, 1.0, 1.0, 1.0, 0.0), equations)
end

initial_condition = initial_condition_shifted_weak_blast_wave

###############################################################################
# Get the DG approximation space

volume_flux = flux_derigs_etal
solver = DGSEM(polydeg=6, surface_flux=FluxRotated(flux_derigs_etal),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

###############################################################################
# Get the unstructured quad mesh from a file (downloads the file if not available locally)

default_mesh_file = joinpath(@__DIR__, "mesh_alfven_wave_with_twist_and_flip.mesh")
isfile(default_mesh_file) || download("https://gist.githubusercontent.com/andrewwinters5000/8f8cd23df27fcd494553f2a89f3c1ba4/raw/85e3c8d976bbe57ca3d559d653087b0889535295/mesh_alfven_wave_with_twist_and_flip.mesh",
                                      default_mesh_file)
mesh_file = default_mesh_file

mesh = UnstructuredMesh2D(mesh_file, periodicity=true)

###############################################################################
# create the semi discretization object

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false,
                                     extra_analysis_integrals=(entropy, energy_total,
                                                               energy_kinetic, energy_internal,
                                                               energy_magnetic, cross_helicity))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)
cfl = 1.0
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        save_solution,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module MemorylessNonlinearities

using AlphaStableDistributions, DelimitedFiles, Interpolations, QuadGK

export
    Arctangent,
    Blanking,
    CauchyNL,
    Clipping,
    HampelThreePart,
    SαSNL,
    SoftClipping,
    TurkeyBiweight,
    InverseHyperbolicSine,

    filt,
    minmaxrescale,

    NLnames

include("utils.jl")

const DATADIR = joinpath(dirname(pathof(@__MODULE__)), "..", "data")

const _xsas = -200:0.1:200
const _αsas = 1.0:0.001:2.0
const _ysas = readdlm(joinpath(DATADIR, "sas.csv"), ',', Float64)

abstract type AbstractMemorylessNonlinearity end

struct Arctangent{T<:Real} <: AbstractMemorylessNonlinearity
    α::T
end

struct Blanking{T<:Real} <: AbstractMemorylessNonlinearity
    k::T
end

struct CauchyNL{T<:Real} <: AbstractMemorylessNonlinearity
    k::T
end

struct Clipping{T<:Real} <: AbstractMemorylessNonlinearity
    k::T
end

struct HampelThreePart{T<:Real} <: AbstractMemorylessNonlinearity
    a::T
    b::T
    c::T
end

struct SαSNL{T<:Real} <: AbstractMemorylessNonlinearity
    α::T
    scale::T
    location::T
    approx::Bool
end
SαSNL(α, scale, location) = SαSNL(α, scale, location, true)

struct SoftClipping{T<:Real} <: AbstractMemorylessNonlinearity
    k::T
end

struct TurkeyBiweight{T<:Real} <: AbstractMemorylessNonlinearity
    k::T
end

struct InverseHyperbolicSine{T<:Real} <: AbstractMemorylessNonlinearity
    θ::T
end

const NLnames = [:Arctangent,
                 :Blanking,
                 :CauchyNL,
                 :Clipping,
                 :HampelThreePart,
                 :SαSNL,
                 :SoftClipping,
                 :TurkeyBiweight,
                 :InverseHyperbolicSine]

"""
Arctangent nonlinearity.
"""
arctangent(x::S, α::T) where {S<:Real,T<:Real} = (2 / π) * atan(α * x)
function filt(f::Arctangent, x::AbstractVector)
    arctangent.(x, f.α)
end

"""
Blanking nonlinearity.
"""
blanking(x::S, k::T) where {S<:Real,T<:Real} = abs(x) > k ? S(0) : x
function filt(f::Blanking, x::AbstractVector)
    blanking.(x, f.k)
end

"""
Cauchy nonlinearity.
"""
cauchynl(x::S, k::T) where {S<:Real,T<:Real} = (2 * (x / k)) / (1 + (x / k)^2)
function filt(f::CauchyNL, x::AbstractVector)
    cauchynl.(x, f.k)
end

"""
Clipping nonlinearity.
"""
clipping(x::S, k::T) where {S<:Real,T<:Real} = abs(x) > k ? sign(x) * k : x
function filt(f::Clipping, x::AbstractVector)
    clipping.(x, f.k)
end

"""
Hampel's Three Part Redescending nonlinearity.
"""
function hampelthreepart(x::S, a::T, b::T, c::T) where {S<:Real,T<:Real}
    absx = abs(x)
    if absx < a
        return x
    elseif absx >= a && absx < b
        return sign(x) * a
    elseif absx >= b && absx < c
        return (absx - c) * (a / (b - c)) * sign(x)
    else
        return 0.
    end
end
function filt(f::HampelThreePart, x::AbstractVector)
    hampelthreepart.(x, f.a, f.b, f.c)
end

"""
Symmetric Alpha Stable nonlinearity.

`x` is standard symmetric alpha stable distributed.
"""
function sαsnl(x::S, α) where {S<:Real}
    num, _ = quadgk(t -> t * sin(t * x) * exp(-t^α), 0, Inf; atol=1e-12, order=3)
    den, _ = quadgk(t -> cos(t * x) * exp(-t^α), 0, Inf; atol=1e-12, order=3)
    num / den
end
function filt(f::SαSNL, x::AbstractVector)
    (abs(f.α - 2.0) < 0.001) && return (x ./ 2) # Gaussian
    (abs(f.α - 1.0) < 0.001) && return (2 .* x) ./ (1 .+ x.^2) # Cauchy
    xstd = (collect(x) .- f.location) ./ f.scale
    if f.approx
        maxval = maximum(_xsas)
        minval = minimum(_xsas)
        xstd[xstd .> maxval] .= (f.α+1) ./ xstd[xstd .> maxval]
        xstd[xstd .< minval] .= (f.α+1) ./ xstd[xstd .< minval]
        itp = interpolate((_αsas, _xsas), _ysas, Gridded(Linear()))
        xstd[(xstd .<= maxval) .& (xstd .>= minval)] .= itp(f.α, xstd[(xstd .<= maxval) .& (xstd .>= minval)])
        xstd
    else
        sαsnl.(xstd, f.α)
    end
end

"""
Soft clipping nonlinearity.
"""
function softclipping(x::S, k::T) where {S<:Real,T<:Real}
    abs(x) > k ? sign(x) * k : k * (3*(x / k) / 2) * (1- ((x/ k)^2) / 3)
end
function filt(f::SoftClipping, x::AbstractVector)
    softclipping.(x, f.k)
end

"""
Turkey's biweight nonlinearity.
"""
turkeybiweight(x::S, k::T) where {S<:Real,T<:Real} = abs(x) > k ? 0. : x * (1 - (x / k)^2)^2
function filt(f::TurkeyBiweight, x::AbstractVector)
    turkeybiweight.(x, f.k)
end

"""
Inverse hyperbolic sine nonlinearity.
"""
inversehyperbolicsine(x::S, θ::T) where {S<:Real,T<:Real} = asinh(x*θ)
function filt(f::InverseHyperbolicSine, x::AbstractVector)
    inversehyperbolicsine.(x, f.θ)
end

end


text:
#!/usr/bin/env julia

import Celeste.AccuracyBenchmark
import Celeste.GalsimBenchmark

test_case_names = String[]
for arg in ARGS
    if ismatch(r"^--test_case_names", arg)
        # test case name: function name in test_case_definitions.py, or CL_DESCR
        # field from FITS header
        test_case_names =
            split(replace(arg, "--test_case_names=", ""), ",")
    else
        println("Usage: run_galsim_benchmark.jl " *
                "[--test_case_names=test1,test2,...]")
        error("Invalid argument: ", arg)
    end
end

srand(12345)

truth_catalog, single_predictions = GalsimBenchmark.run_benchmarks(
    test_case_names=test_case_names,
    joint_inference=false,
)
unused, joint_predictions = GalsimBenchmark.run_benchmarks(
    test_case_names=test_case_names,
    joint_inference=true,
)

score_df = AccuracyBenchmark.score_predictions(
    truth_catalog,
    [single_predictions, joint_predictions],
)
println(repr(score_df))


text:
# Autogenerated wrapper script for battery_cli_jll for armv7l-linux-musleabihf
export battery

JLLWrappers.@generate_wrapper_header("battery_cli")
JLLWrappers.@declare_executable_product(battery)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        battery,
        "bin/battery",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Dates

module TimeSeries

using Dates

export TimeArray, AbstractTimeSeries,
       by, from, to, findwhen, findall, timestamp, values, colnames,
       lag, lead, percentchange, moving, upto,
       basecall,
       merge, collapse,
       readtimearray

###### include ##################

include(".timeseriesrc.jl")
include("timearray.jl")
include("split.jl")
include("apply.jl")
include("combine.jl")
include("readwrite.jl")

end


text:
using AppleAccelerate
using DSP

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

srand(7)
N = 1_000

for T in (Float32, Float64)
    @testset "Element-wise Operators::$T" begin
        X::Vector{T} = randn(N)
        Y::Vector{T} = randn(N)
        @test (X .+ Y) ≈ AppleAccelerate.vadd(X, Y)
        @test (X .- Y) ≈ AppleAccelerate.vsub(X, Y)
        @test (X .* Y) ≈ AppleAccelerate.vmul(X, Y)
        @test (X ./ Y) ≈ AppleAccelerate.vdiv(X, Y)
    end
end

for T in (Float32, Float64)
    @testset "Rounding::$T" begin
        X::Array{T} = 100*randn(N)
        @testset "Testing $f::$T" for f in [:floor,:ceil,:trunc,:round]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end
    end
end

for T in (Float32, Float64)
    @testset "Logarithmic::$T" begin
        X::Array{T} = exp(10*randn(N))
        @testset "Testing $f::$T" for f in [:log,:log2,:log10, :log1p]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end
    end
end

for T in (Float32, Float64)
    @testset "Exponential::$T" begin
        @testset "Testing $f::$T" for f in [:exp,:exp2,:expm1]
            X::Array{T} = 10*randn(N)
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end
    end
end

for T in (Float32, Float64)
    X::Array{T} = 10*randn(N)
    @testset "Trigonometric::$T" begin
        @testset "Testing $f::$T" for f in [:sin,:sinpi,:cos,:cospi,:tan,:atan] # tanpi not defined in Base
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end

        Y::Array{T} = 10*randn(N)
        @testset "Testing $f::$T" for f in [:atan2]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X,Y) ≈ fb(X,Y)
        end

        Z::Array{T} = 2*rand(N)-1
        @testset "Testing $f::$T" for f in [:asin,:acos]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(Z) ≈ fb(Z)
        end
    end
end

for T in (Float32, Float64)
    @testset "Hyperbolic::$T" begin
        X = 10*randn(N)
        @testset "Testing $f::$T" for f in [:sinh,:cosh,:tanh,:asinh]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end

        Y = exp(10*randn(N))+1
        @testset "Testing $f::$T" for f in [:acosh]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(Y) ≈ fb(Y)
        end

        Z = 2*rand(N)-1
        @testset "Testing $f::$T" for f in [:atanh]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(Z) ≈ fb(Z)
        end
    end
end

@testset "DCT::Float32" begin
    r=rand(Float32,2^16)
    d1=dct(r)
    plan_accel = AppleAccelerate.plan_dct(length(r), 2)
    d2=AppleAccelerate.dct(r, plan_accel)
    @test norm(d1[2]/d2[2]*d2[2:end]-d1[2:end])≤1000eps(Float32)
end

for T in (Float32,  Float64)
    @testset "Convolution & Correlation::$T" begin
        X::Vector{T} = randn(N)
        Y::Vector{T} = randn(N)
        @testset "Testing $f::$T" for f in [:conv, :xcorr]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fb(X, Y) ≈ fa(X, Y)
        end

        @testset "Testing $f::$T" for f in [:xcorr]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fb(X, copy(X)) ≈ fa(X)
        end
    end
end

for T in (Float64, )
    @testset "Biquadratic Flitering::$T" begin
        @testset "Single Section::$T" begin
            X::Vector{T} = randn(10)
            d::Vector{T} = zeros(4)
            c::Vector{T} = [x%0.5 for x in randn(5)]
            fdsp = DSP.Biquad(c[1], c[2], c[3], c[4], c[5])
            fa = AppleAccelerate.biquadcreate(c, 1)
            @test DSP.filt(fdsp, X) ≈ AppleAccelerate.biquad(X, d, length(X), fa)
        end
    end
end

for T in (Float32, Float64)
    @testset "Window Functions::$T" begin
        N = 64
        @testset "Testing blackman::$T" begin
            Wa = AppleAccelerate.blackman(N, T)
            Wb = Array(T, N)
            for n in 1:N
                Wb[n] = 0.42-(0.5cos(2pi*(n-1)/N)) + (0.08cos(4pi*(n-1)/N))
            end
            @test Wa ≈ Wb
        end

        @testset "Testing hamming::$T" begin
            Wa = AppleAccelerate.hamming(N, T)
            Wb = Array(T, N)
            for n in 1:N
                Wb[n] = 0.54-0.46cos(2pi*(n-1)/N)
            end
            @test Wa ≈ Wb
        end

        @testset "Testing hanning::$T" begin
            Wa = AppleAccelerate.hanning(N, T)
            Wb = Array(T, N)
            for n in 1:N
                Wb[n] = 0.5(1.0-cos(2pi*(n-1)/N))
            end
            @test Wa ≈ Wb
        end
    end
end

for T in (Float32, Float64)
    @testset "Array Properties::$T" begin
        X::Array{T} = randn(N)
        @testset "Testing $f::$T" for f in [:maximum, :minimum, :mean, :sum]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end

        @testset "Testing $f::$T" for f in [:findmax, :findmin]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X)[1] ≈ fb(X)[1]
            @test fa(X)[2] ≈ fb(X)[2]
        end

        @testset "Testing meansqr::$T" begin
            @test AppleAccelerate.meansqr(X) ≈ mean(X .*X)
        end

        @testset "Testing meanmag::$T" begin
            @test AppleAccelerate.meanmag(X) ≈ mean(abs(X))
        end

    end
end

for T in (Float32, Float64)
    @testset "Misc::$T" begin
        X::Array{T} = exp(10*randn(N))
        @testset "Testing $f::$T" for f in [:sqrt]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X) ≈ fb(X)
        end

        Y::Array{T} = 10*randn(N)
        @testset "Testing $f::$T" for f in [:exponent, :abs]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(Y) ≈ fb(Y)
        end

        Z::Array{T} = 10*randn(N)
        @testset "Testing $f::$T" for f in [:copysign]
            @eval fb = $f
            @eval fa = AppleAccelerate.$f
            @test fa(X,Y) ≈ fb(X,Y)
        end

        @test AppleAccelerate.rem(X,Y) == [rem(X[i], Y[i]) for i=1:length(X)]

    end
end

for T in (Float32, Float64)
    @testset "Extra::$T" begin
        X::Array{T} = randn(N)
        Y::Array{T} = abs(randn(N))

        @test AppleAccelerate.rec(X) ≈ 1./X
        @test AppleAccelerate.rsqrt(Y) ≈ 1./sqrt(Y)
        @test AppleAccelerate.pow(Y,X) ≈ Y.^X
        @test AppleAccelerate.fdiv(X,Y) ≈ X./Y

        @test AppleAccelerate.sincos(X)[1] ≈ sin(X)
        @test AppleAccelerate.sincos(X)[2] ≈ cos(X)
        @test AppleAccelerate.cis(X) ≈ cis(X)

    end
end

AppleAccelerate.@replaceBase(sin, atan2, ./)

@testset "Replace Base::$T" for T in (Float32, Float64)
    X::Array{T} = randn(N)
    Y::Array{T} = abs(randn(N))

    @test Base.sin(X) == AppleAccelerate.sin(X)
    @test Base.atan2(X, Y) == AppleAccelerate.atan2(X, Y)
    @test X ./ Y  == AppleAccelerate.fdiv(X, Y)
end


text:
module ObjectPools
using Base.Threads

export allocate, release

const MAX_POOL_SIZE = 64

const pools = IdDict{Any, Any}()
const poolslock = SpinLock()

function _createpool_unsafe(type::Type{T}) where T
    pool = Vector{T}()
    sizehint!(pool, MAX_POOL_SIZE)
    pools[T] = pool
    return pool
end

# Eventually thread-safe: This solution may recreate a pool when threads race on the same
# type, but subsequent getpool calls will return the last one and the overwitten ones will
# be gc-d.
function _createpool(type)
    lock(poolslock)
    try
        return _createpool_unsafe(type)
    finally
        unlock(poolslock)
    end
end

function _trycreatepool(type)
    if trylock(poolslock)
        try
            return _createpool_unsafe(type)
        finally
            unlock(poolslock)
        end
    else
        return nothing
    end
end

function getpool(@nospecialize type)
    if !haskey(pools, type)
        _createpool(type)
    end
    return pools[type]::Vector{type}
end

@generated function getpool_static(type::Type{T}) where T
    if !haskey(pools, T)
        pool = _trycreatepool(T)
        if isnothing(pool) # Wasn't able to create the pool because pools is locked
            return quote getpool($T) end # Falling back to dynamic
        end
    end
    return quote
        return $(getpool(T))
    end
end

function allocate(type::Type{T}, args...) where T
    pool = getpool_static(T)
    return allocate(pool, T, args...)
end

function allocate(pool, type::Type{T}, args...) where T
    if length(pool) > 0
        return pop!(pool)
    end
    return T(args...)
end

function release(obj::T) where T
    pool = getpool_static(T)
    return release(pool, obj)
end

function release(pool, obj)
    if length(pool) < MAX_POOL_SIZE
        push!(pool, obj)
        return true
    end
    return false
end

end # module


text:
export newTemplate, render2file
################################################################
"""
    newTemplate(name)

Create new destination html file as the template

    newTemplate(name, :function)

Prints a function to be used as a template

# Examples
```julia
# you can create a file and edit the file directly by using
newTemplate("person")

## create person function to store out html template
newTemplate("person", :function)
```
"""
function newTemplate(name, type::Symbol = :file, template::Union{String,Nothing} = nothing)

    if type == :file

        try
            mkpath(joinpath("app", "destination"))
        catch e

        end
        filePath = joinpath("app", "destination", "$(name).html")
        touch(filePath)

        if !isnothing(template)

            file = open(filePath, "w")
            print(file, template)
            close(file)

        else # open file in Atom for editing
            if !isdefined(Main, :Atom)
                Atom.edit(filePath)
            end
        end

    elseif type == :function

        if !isnothing(template) # doesn't work because $ are evaluated outside of the function
            prt = template
        else
            prt = """
            # enter your html code here
            # use \$() to evaluate julia variables or code
            """
        end

    # todo: append this to controller file
        println("""

        function $(name)(; put_variables_here)

        # use your any julia code here to control the outputed html string.

        \"\"\"
        $prt
        \"\"\"

        end
        """)

    end

end
################################################################

"""
    render2file(destination, overwrite, var...)

render variables passed as an input to the destination file.

You should put \$(var) in the destination file/string so var is evaluated there. Pass the variables as keyword arguments with the same name you used in the html string/file. Variables should be string,

If you want to statically overwrite the file pass `true` as the 2nd argument to the function. Useful if you don't want a dynamic website.

# Examples
```julia
# Add the following html code to the generated html file
#=
<person id=\$(id)>
  <age>\$(age)</age>
  <study-field>\$(field)</study-field>
  <GPA>\$(GPA)</GPA>
  <taken-courses>\$(courses[1])</taken-courses>
  <taken-courses>\$(courses[2])</taken-courses>
</person>
=#

# Specify the template (or its path), and also the variables for rendering
out =render2file("person", false,
  id = 1,
  age = 24,
  field = "Mechanical Engineering",
  GPA = 4.5,
  courses = ["Artificial Intelligence", "Robotics"])

# you pass `true` as the 2nd argument to owerwrite person.html statically.

```
"""
function render2file(destination::String, overwrite::Bool=false; var...)

    local file, destinationString, s

    # creates variables from keywords inputed to the function
    for key in keys(var)
        @eval $key = multiString($(var[key]))
    end

    # reading string from html file
    if ispath(destination)

        try
            filePath = destination
            file = open(filePath, "r")
            destinationString = read(file, String)

        catch e
            error("destination file address is not correct")
        end

    else
        filePath = joinpath("app", "destination", "$(destination).html")
        file = open(filePath, "r")
        destinationString = read(file, String)
    end

    try
        s = '"' * '"' * '"' * destinationString * '"' * '"' * '"'
    catch e
        s = destinationString
    end

    # evaluate $() values in the string
    out = eval(Meta.parse(s))

    if overwrite
        print(file, out)
    end
    close(file)

    return out
end

################################################################
# Multi-Stringer

multiString(var::String) = var
multiString(var::Union{Number,Bool}) = string(var)

function multiString(var::Array{T}) where {T<:Union{String, Number, Bool}}

    varNum=length(var)

    # out = NTuple{varNum,String}

    out = Array{String}(undef,varNum)

    for i = 1:varNum
        out[i] = string(var[i])
        # Base.setindex!(out,var[i],i)
    end

    return out
end

function multiString(var::NTuple{N,T}) where {T<:Union{String, Number, Bool}} where{N}

    # out = NTuple{N,String}
    out = Array{String}(undef,N)

    for i = 1:N
        out[i] = string(var[i])
        # Base.setindex!(out,var[i],i)
    end

    return Tuple(out)
end
################################################################


text:
# some definitions used throughout TIcon

"""
  Abstract type for _AssembleElementData.  See that type for details.
"""
abstract type AssembleElementData end

#------------------------------------------------------------------------------
# sparsity pattern enums

"""
  disc_type constant for SparseMatrixCSC constructor, indicating the face nodes
  of an element are connected to the face nodes of the element that shares
  the face.
"""
global const INVISCID=1

"""
  disc_type constant for SparseMatrixCSC constructor, indicating all nodes
  of an element are connected to all face nodes of its neighbours
"""
global const VISCOUS=2

"""
  disc_type constant for SparseMatrixCSC constructor, indicating all nodes
  of an element are connected to all nodes of its neighbors
"""
global const COLORING=3

"""
  disc_type constant for sparsity pattern of viscous terms when the T4 = 0.
  In this case, all nodes in the stencil of the interpolation operator are
  connected to all nodes of the adjacent element.
"""
global const VISCOUSTIGHT=4

#------------------------------------------------------------------------------
# Shock Capturing stuff - see the docs in the Euler module

abstract type AbstractShockSensor end


text:
exp_decay(i::Integer, tot::Integer, a::Real=10., b::Real=0.) = (1-b)*exp(-a*i/tot)+b

pow_decay(i::Integer, tot::Integer, a::Real=1e-3, b::Real=0.) = (1-b)*(a^(i/tot))+b

lin_decay(i::Integer, tot::Integer, a::Real=1e-3, b::Real=0.) = (1-b)*(1/(1+a*i))+b

quad_decay(i::Integer, tot::Integer, a::Real=1e-3, b::Real=0.) = (1-b)*(1/(1+a*abs2(i)))+b


text:
using FunManifolds
using Test
using StaticArrays

include("../utils.jl")

@testset "DiscretizedCurves" begin
    N = 50
    s2 = Sphere(2)
    r2 = Euclidean(2)

    dcss2 = DiscretizedCurves(s2, range(0.0, 1.0, length=N))
    test_manifold(
        dcss2,
        [
            discretize(dcss2, t -> project(s2, @SVector [t + 0.1, t^2, sin(t)])),
            discretize(dcss2, t -> project(s2, @SVector [t - 1, sin(t^2), 1])),
            discretize(dcss2, t -> project(s2, @SVector [t^2, t, cos(t)])),
        ];
        exp_log_atol_multiplier=10.0,
        is_tangent_atol_multiplier=10,
        test_default_vector_transport=true,
        test_project_point=true,
        test_project_tangent=true,
    )

    dcsr2 = DiscretizedCurves(r2, range(0.0, 1.0, length=N))
    test_manifold(
        dcsr2,
        [
            discretize(dcsr2, t -> project(r2, @SVector [t^2, sin(t)])),
            discretize(dcsr2, t -> project(r2, @SVector [sin(t^2), 1])),
            discretize(dcsr2, t -> project(r2, @SVector [t, cos(t)])),
        ];
        is_tangent_atol_multiplier=10,
        test_default_vector_transport=true,
    )

    @testset "SRVF" begin
        r3 = Euclidean(3)
        dcsr1 = DiscretizedCurves(r3, range(0.0, 1.0, length=9))
        f = t -> (@SVector [sin(t), t^2, 2 * t + 1])
        c = discretize(dcsr1, f)
        q = tsrvf(
            dcsr1,
            c,
            (@SVector [0.0, 0.0, 0.0]),
            FunManifolds.ProjectedDifferenceBackend(nothing),
        )
        c_rev = reverse_srvf(dcsr1, q, f(0.0))
        @test c ≈ c_rev atol = 1e-10
    end
end


text:
import ScikitLearnBase
using ScikitLearnBase: @declare_hyperparameters

export SkGLRM, PCA, QPCA, NNMF, KMeans, RPCA

################################################################################
# Shared definitions

# Note: there is redundancy in the hyperparameters. This is
# necessary if we want to offer a simple interface in PCA(), and a full
# interface in SkGLRM(). PCA(abs_tol=0.1, max_iter=200) cannot create
# `ProxGradParams(abs_tol, max_iter)` right away, because abs_tol and
# max_iter are hyperparameters and need to be visible/changeable by
# set_params for grid-search.
# There are other ways of setting it up, but this seems like the simplest.
mutable struct SkGLRM <: ScikitLearnBase.BaseEstimator
    # Hyperparameters: those will be passed to GLRM, so it doesn't matter if
    # they're not typed.
    fit_params # if fit_params != nothing, it has priority over abs_tol, etc.
    loss
    rx
    ry
    # rx/ry_scale can be nothing, in which case they're ignored. This allows
    # ry to be a vector
    rx_scale
    ry_scale
    abs_tol::Float64
    rel_tol::Float64
    max_iter::Int
    inner_iter::Int
    k::Int
    init::Function   # initialization function
    verbose::Bool

    glrm::GLRM     # left undefined by the constructor
end

# This defines `clone`, `get_params` and `set_params!`
@declare_hyperparameters(SkGLRM, [:fit_params, :init, :rx, :ry,
                                 :rx_scale, :ry_scale, :loss,
                                 :abs_tol, :rel_tol, :max_iter, :inner_iter, :k,
                                 :verbose])

function do_fit!(skglrm::SkGLRM, glrm::GLRM)
    fit_params = (skglrm.fit_params === nothing ?
                  ProxGradParams(abs_tol=skglrm.abs_tol,
                                 rel_tol=skglrm.rel_tol,
                                 max_iter=skglrm.max_iter) :
                  skglrm.fit_params)

    fit!(glrm, fit_params; verbose=skglrm.verbose)
end

function ind2sub(a, i)
    i2s[i]
end

function build_glrm(skglrm::SkGLRM, X, missing_values)
    k = skglrm.k == -1 ? size(X, 2) : skglrm.k
    i2s = CartesianIndices(missing_values)
    obs = [i2s[x] for x in  (LinearIndices(.!missing_values))[findall(.!missing_values)] ]
    rx, ry = skglrm.rx, skglrm.ry
    if skglrm.rx_scale !== nothing
        rx = copy(rx)
        mul!(rx, skglrm.rx_scale)
    end
    if skglrm.ry_scale !== nothing
        ry = copy(ry)
        mul!(ry, skglrm.ry_scale)
    end
    GLRM(X, skglrm.loss, rx, ry, k; obs=obs)
end

# The input matrix is called X (instead of A) following ScikitLearn's convention
function ScikitLearnBase.fit_transform!(skglrm::SkGLRM, X, y=nothing;
                                        missing_values=isnan.(X))
    @assert size(X)==size(missing_values)

    # Reuse the standard GLRM constructor and fitting machinery
    skglrm.glrm = build_glrm(skglrm, X, missing_values)
    skglrm.init(skglrm.glrm)
    X, _, _ = do_fit!(skglrm, skglrm.glrm)

    return X'
end

function ScikitLearnBase.fit!(skglrm::SkGLRM, X, y=nothing; kwargs...)
    ScikitLearnBase.fit_transform!(skglrm, X; kwargs...)
    skglrm
end

""" `transform(skglrm::SkGLRM, X)` brings X to low-rank-space """
function ScikitLearnBase.transform(skglrm::SkGLRM, X;
                                   missing_values=isnan.(X))
    glrm = skglrm.glrm
    ry_fixed = [FixedLatentFeaturesConstraint(glrm.Y[:, i])
                for i=1:size(glrm.Y, 2)]
    glrm_fixed = build_glrm(skglrm, X, missing_values)
    X2, _, ch = do_fit!(skglrm, glrm_fixed)
    return X2'
end

""" `transform(skglrm::SkGLRM, X)` brings X from low-rank-space back to the
original input-space """
ScikitLearnBase.inverse_transform(skglrm::SkGLRM, X) = X * skglrm.glrm.Y

# Only makes sense for KMeans
function ScikitLearnBase.predict(km::SkGLRM, X)
    X2 = ScikitLearnBase.transform(km, X)
    # This performs the "argmax" over the columns to get the cluster #
    return mapslices(argmax, X2, 2)[:]
end

################################################################################
# Public constructors

"""
    SkGLRM(; fit_params=nothing, init=glrm->nothing, k::Int=-1,
           loss=QuadLoss(), rx::Regularizer=ZeroReg(), ry=ZeroReg(),
           rx_scale=nothing, ry_scale=nothing,
           # defaults taken from proxgrad.jl
           abs_tol=0.00001, rel_tol=0.0001, max_iter=100, inner_iter=1,
           verbose=false)

Generalized low rank model (GLRM). GLRMs model a data array by a low rank
matrix. GLRM makes it easy to mix and match loss functions and regularizers to
construct a model suitable for a particular data set.

Hyperparameters:

- `fit_params`: algorithm to use in fitting the GLRM. Defaults to
   `ProxGradParams(abs_tol, rel_tol, skglrm.max_iter)`
- `init`: function to initialize the low-rank matrices, before the main gradient
   descent loop.
- `k`: number of components (rank of the latent representation). By default,
   use k=nfeatures (full rank)
- `loss`: loss function. Can be either a single `::Loss` object, or a vector
   of `nfeature` loss objects, allowing for mixed inputs (eg. binary and
   continuous data)
- `rx`: regularization over the hidden coefficient matrix
- `ry`: regularization over the latent features matrix. Can be either a single
   regularizer, or a vector of regularizers of length nfeatures, allowing
   for mixed inputs
- `rx_scale`, `ry_scale`: strength of the regularization (higher is stronger).
   By default, `scale=1`. Cannot be used if `rx/ry` are vectors.
- `abs_tol, rel_tol`: tolerance criteria to stop the gradient descent iteration
- `max_iter, inner_iter`: number of iterations in the gradient descent loops
- `verbose`: print convergence information

All parameters (in particular, `rx/ry_scale`) can be tuned with
`ScikitLearn.GridSearch.GridSearchCV`

For more information on the parameters see [LowRankModels](https://github.com/madeleineudell/LowRankModels.jl)
"""
function SkGLRM(; fit_params=nothing, init=glrm->nothing, k=-1,
                loss=QuadLoss(), rx=ZeroReg(), ry=ZeroReg(),
                rx_scale=nothing, ry_scale=nothing,
                # defaults taken from proxgrad.jl
                abs_tol=0.00001, rel_tol=0.0001, max_iter=100, inner_iter=1,
                verbose=false)
    dummy = pca(zeros(1,1), 1) # it needs an initial value - will be overwritten
    return SkGLRM(fit_params, loss, rx, ry, rx_scale, ry_scale, abs_tol,
                  rel_tol, max_iter,
                  inner_iter, k, init, verbose, dummy)
end

"""    PCA(; k=-1, ...)

Principal Component Analysis with `k` components (defaults to using
`nfeatures`). Equivalent to

    SkGLRM(loss=QuadLoss(), rx=ZeroReg(), ry=ZeroReg(), init=init_svd!)

See ?SkGLRM for more hyperparameters. In particular, increasing `max_iter`
(default 100) may improve convergence. """
function PCA(; kwargs...)
    # principal components analysis
    # minimize ||A - XY||^2
    loss = QuadLoss()
    r = ZeroReg()
    return SkGLRM(; loss=loss, rx=r, ry=r, init=init_svd!, kwargs...)
end

"""    QPCA(k=-1, rx_scale=1, ry_scale=1; ...)

Quadratically Regularized PCA with `k` components
(default: `k = nfeatures`). Equivalent to

    SkGLRM(loss=QuadLoss(), rx=QuadReg(1.0), ry=QuadReg(1.0), init=init_svd!)

Regularization strength is set by `rx_scale` and `ry_scale`. See ?SkGLRM for
more hyperparameters.
"""
function QPCA(; kwargs...)
    # quadratically regularized principal components analysis
    # minimize ||A - XY||^2 + rx_scale*||X||^2 + ry_scale*||Y||^2
    loss = QuadLoss()
    r = QuadReg(1.0) # scale is set in build_glrm
    return SkGLRM(; loss=loss, rx=r, ry=r, init=init_svd!, kwargs...)
end

"""    NNMF(; k=-1, ...)

Non-negative matrix factorization with `k` components (default:
`k=nfeatures`). Equivalent to

    SkGLRM(loss=QuadLoss(), rx=NonNegConstraint(), ry=NonNegConstraint(), init=init_svd!)

See ?SkGLRM for more hyperparameters
"""
function NNMF(; kwargs...)
    # nonnegative matrix factorization
    # minimize_{X>=0, Y>=0} ||A - XY||^2
    loss = QuadLoss()
    r = NonNegConstraint()
    return SkGLRM(; loss=loss,rx=r,ry=r, init=init_svd!, kwargs...)
end

"""    KMeans(; k=2, inner_iter=10, max_iter=100, ...)

K-Means algorithm. Separates the data into `k` clusters. See ?SkGLRM for more
hyperparameters. In particular, increasing `inner_iter` and `max_iter` may
improve convergence.

**IMPORTANT**: This is not the most efficient way of performing K-Means, and
the iteration may not reach convergence.
"""
function KMeans(; k=2, inner_iter=10, kwargs...)
    # minimize_{columns of X are unit vectors} ||A - XY||^2
    loss = QuadLoss()
    rx = UnitOneSparseConstraint()
    ry = ZeroReg()
    return SkGLRM(k=k, loss=loss,rx=rx,ry=ry, inner_iter=inner_iter,
                  init=init_kmeanspp!; kwargs...)
end

"""    RPCA(; k=-1, ...)

Robust PCA with `k` components (default: `k = nfeatures`). Equivalent to

    SkGLRM(loss=HuberLoss(), rx=QuadReg(1.0), ry=QuadReg(1.0), init=init_svd!)

Regularization strength is set by `rx_scale` and `ry_scale`. See ?SkGLRM for
more hyperparameters. In particular, increasing `max_iter` (default 100) may
improve convergence. """
function RPCA(; kwargs...)
    # robust PCA
    # minimize HuberLoss(A - XY) + scale*||X||^2 + scale*||Y||^2
    loss = HuberLoss()
    r = QuadReg(1.0)
    return SkGLRM(; loss=loss,rx=r,ry=r, init=init_svd!, kwargs...)
end


text:
#=Copyright (c) 2017 Gabriel Goh

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.=#

λ = 1
#M(α,β) = [β λ; -α*β (1-α)]
M(α,β) = [β λ; -α*(β^2) 1-α*λ*(1+β)]
#M(α,β) = [1 αλ; -β-1 -α*β*λ]
S(α,β) = max(abs((eig(M(α,β))[1]))...)

n = 200
R1 = linspace(0,2,n);
R2 = linspace(-1,1,2*n);
Z = [S(R1[i],R2[j]) for i = 1:n, j = 1:2*n]

figure(); contour(R1, R2, Z', linspace(0,1,30), cmap = "magma")
contour(R1, R2, Z', [1.,])

# figure(); imshow(flipdim(Z',1), cmap = "magma")

xlabel(L"\alpha\lambda_i")
ylabel(L"\beta")

figure(); plot(R2, Z[6,:]',"k")
xlabel(L"\beta")
ylabel("Convergence Rate")

max_stars_repo_name: sbeleidy/post--momentum

max_issues_repo_name: sbeleidy/post--momentum

max_forks_repo_name: sbeleidy/post--momentum


text:
module Visualization

# user-side functionality
export
    SimulationControls,
    GUI,
    setanimation!

using DocStringExtensions
@template (FUNCTIONS, METHODS, MACROS) =
    """
    $(SIGNATURES)
    $(DOCSTRING)
    """

@template (TYPES,) =
    """
    $(TYPEDEF)
    $(DOCSTRING)
    """

using MeshCatMechanisms

import MeshCat
import WebIO

using Printf: @sprintf
using DiffEqBase: DiscreteCallback, ODESolution, CallbackSet, u_modified!, terminate!
using RigidBodyDynamics: Mechanism, MechanismState, normalize_configuration!, configuration, num_positions
using MeshCat: Animation, atframe, setanimation!
using Observables: Observable
using InteractBase: Widget, button, observe
using WebIO: render, node, Node
using Blink: Window, body!, title

using DataStructures: top

struct SimulationStatus
    terminate::Observable{Bool}
    pause::Observable{Bool}
    time::Observable{Float64}
end

function initialize!(status::SimulationStatus)
    status.terminate[] = false
    status.pause[] = false
    status.time[] = 0.0
end

const TIME_DISPLAY_INTERVAL = 1e-2
const DEFAULT_PAUSE_POLLINT = 0.05

function CommandHandler(status::SimulationStatus; pause_pollint::Float64 = DEFAULT_PAUSE_POLLINT)
    condition = (u, t, integrator) -> true
    action = let status = status
        function (integrator)
            yield()
            while status.pause[] && !status.terminate[]
                sleep(pause_pollint)
            end
            status.terminate[] && (status.terminate[] = false; terminate!(integrator))
            if abs(status.time[] - integrator.t) >= TIME_DISPLAY_INTERVAL
                status.time[] = integrator.t
            end
            u_modified!(integrator, false)
        end
    end
    initialize = (c, t, u, integrator) -> initialize!(status)
    DiscreteCallback(condition, action, initialize = initialize, save_positions=(false, false))
end

function TransformPublisher(vis; max_fps = 60.)
    last_update_time = Ref(-Inf)
    condition = let last_update_time = last_update_time, min_Δt = 1 / max_fps
        function (u, t, integrator)
            last_time_step = length(integrator.opts.tstops) == 1 && t == top(integrator.opts.tstops)
            last_time_step || time() - last_update_time[] >= min_Δt
        end
    end
    action = let vis = vis, last_update_time = last_update_time
        function (integrator)
            last_update_time[] = time()
            copyto!(vis, integrator.u)
            u_modified!(integrator, false)
        end
    end
    DiscreteCallback(condition, action; save_positions = (false, false))
end

struct SimulationControls
    terminate::Widget{:button}
    pause::Widget{:button}
    time::Observable{HTML{String}}
end

"""
Create a new `SimulationControls` object, which may be used to pause and terminate the simulation.

The controls can be displayed in a standalone window using `open(controls, Blink.Window())`.
"""
SimulationControls() = SimulationControls(button(""), button(""), Observable(HTML("0.0")))

# Icons taken from the freely available set at https://icons8.com/color-icons/
function WebIO.render(controls::SimulationControls)
    node(:div,
        node(:style, """
            .rigidbodysim-controls button {
                height: 2.5em;
                width: 2.5em;
                margin: 0.1em;
                padding: 0.2em 0.2em;
                image-rendering: pixelated;
                user-select: none;
            }
            .rigidbodysim-controls button:active {
                filter: brightness(50%);
            }
            .rigidbodysim-controls-pause button {
                background: url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAADAAAAAwCAYAAABXAvmHAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAADUSURBVGhD7Y6xCQMBEMN+lSydRTNAgkGdVR8ELFAlF37GGGOMMcb4T17vz9ckF7aN5HvsTCQXto3ke+xMJBe2jeR77EwkF7aN5HvsTCQXto3ke+xMJBe2jeR77EwkF7aN5HvsTCQXto3ke+xMJBe2jeR77EwkF7aN5HvsTCQXto3ke+xMJBe2jeR77EwkF7aN5HvsTCQXto3ke+xMJBe2jeR77EwkF7aN5HvsTCQXto3ke+xMJBe2jeR77EwkF7aN5HvsTCQXto3kMcYYY4wx/orn+QEf136Pgv8IWQAAAABJRU5ErkJggg==) no-repeat;
                background-size: cover;
            }
            .rigidbodysim-controls-terminate button {
                background: url(data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAADAAAAAwCAYAAABXAvmHAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsMAAA7DAcdvqGQAAACQSURBVGhD7c6xCcBAAMPAXyz7j5U0Kt05GB50oNb4SJIkSfd5n/P+GbM76UQTszvpRBOzO+lEE7M76UQTszvpRBOzO+lEE7M76UQTszvpRBOzO+lEE7M76UQTszvpRBOzO+lEE7M76UQTszvpRBOzO+lEE7M76UQTszvpRBOzO+lEE7M76UQTs5IkSdI1zvkACnV8XnyMZlcAAAAASUVORK5CYII=) no-repeat;
                background-size: cover;
            }
            .rigidbodysim-controls > div {flex-grow: 0; ; }
        """),
        node(:div,
            node(:div,
                node(:div, "Time:"),
                render(controls.time),
                style = Dict(:fontSize => "10pt", :fontFamily => "sans-serif",
                            :userSelect => "none", :cursor => "default",
                            :height => "2.5em", :width => "8em", :marginLeft => "0.5em", :marginRight => "0.5em", :lineHeight => "2.5em",
                            :display => "flex", :justifyContent => "space-between")
            ),
            node(:div, controls.pause, attributes = Dict(:class => "rigidbodysim-controls-pause")),
            node(:div, controls.terminate, attributes = Dict(:class => "rigidbodysim-controls-terminate")),
            attributes = Dict(:class => "rigidbodysim-controls"),
            style = Dict(:display => "flex", :flexWrap => "wrap", :minHeight => "2.5em", :overflow => "hidden")
        )
    )
end

function Base.open(controls::SimulationControls, window::Window)
    size(window, 300, 70)
    title(window, "RigidBodySim controls")
    body!(window, render(controls))
    nothing
end

function SimulationStatus(controls::SimulationControls)
    terminate = map!(!iszero, Observable(false), observe(controls.terminate))
    pause = map!(isodd, Observable(false), observe(controls.pause))
    time = Observable(0.0)
    status = SimulationStatus(terminate, pause, time)
    map!(t -> HTML(@sprintf("%4.2f", t)), controls.time, status.time)
    status
end

CallbackSet(controls::SimulationControls) = CommandHandler(SimulationStatus(controls))

struct GUI{V}
    visualizer::V
    controls::SimulationControls
    usernode::Union{Node, Nothing}
end

"""
Create a new RigidBodySim graphical user interface given a visualizer (e.g. a `MeshCatMechanisms.MechanismVisualizer`).

The visualizer must support:

* `Base.copyto!(visualizer, state::Union{MechanismState, AbstractVector})`
* `Base.wait(visualizer)`
* `MeshCatMechanisms.visualizer(visualizer)`, which should return a `MeshCat.Visualizer`.

Use `open(gui)` to open the GUI in a standalone window.
"""
GUI(visualizer; usernode=nothing) = GUI(visualizer, SimulationControls(), usernode)

"""
Create a new RigidBodySim graphical user interface for the given Mechanism.
All arguments are passed on to the `MeshCatMechanisms.MechanismVisualizer` constructor.

Use `open(gui)` to open the GUI in a standalone window.
"""
GUI(mechanism::Mechanism, args...; usernode=nothing) = GUI(MechanismVisualizer(mechanism, args...); usernode=usernode)

function Base.open(gui::GUI, window::Window)
    title(window, "RigidBodySim")
    body = node(:div,
        WebIO.render(gui.controls),
        WebIO.iframe(MeshCatMechanisms.visualizer(gui.visualizer).core, minHeight="0"),
        WebIO.render(gui.usernode),
        style = Dict(:display => "flex", :flexDirection => "column", :height => "100vh")
    )
    body!(window, body)
    wait(gui)
    nothing
end

Base.open(gui::GUI) = (window = Window(); open(gui, window); window)
Base.wait(gui::GUI) = wait(gui.visualizer)

"""
Create the DifferentialEquations.jl callbacks associated with the [`GUI`](@ref).

`max_fps` is the maximum number of frames per second (in terms of wall time) to draw. Default: `60.0`.
"""
CallbackSet(gui::GUI; max_fps=60) = CallbackSet(CallbackSet(gui.controls), TransformPublisher(gui.visualizer; max_fps = max_fps))

"""
Create the DifferentialEquations.jl callbacks needed for publishing to a
visualizer during simulation.

`max_fps` is the maximum number of frames per second (in terms of wall time) to draw. Default: `60.0`.
"""
CallbackSet(vis::MechanismVisualizer; max_fps = 60.) = CallbackSet(TransformPublisher(vis; max_fps = max_fps))

"""
Create a `MeshCat.Animation` from an `ODESolution` obtained by `solve!`ing an
`ODEProblem` created using this package.

* `vis` is a `MeshCatMechanisms.MechanismVisualizer`
* `sol` is a `DiffEqBase.ODESolution` obtained from a RigidBodySim.jl simulation.

Keyword arguments:

* `fps`: the frame rate to be used for keyframes.
* `realtime_rate`: can be used to slow down or speed up playback compared to wall time. A `realtime_rate` of `2`
  will result in playback that is sped up 2x. Default: `1`.
"""
function MeshCat.Animation(mvis::MechanismVisualizer, sol::ODESolution; fps::Int=30, realtime_rate::Number=1)
    t0, tf = first(sol.t), last(sol.t)
    animation = Animation(fps)
    # MeshCat animations don't support a realtime_rate option
    # (although it can be adjusted in the GUI), so we instead just
    # do some fps scaling.
    num_frames = floor(Int, (tf - t0) * fps / realtime_rate)
    for frame in 0 : num_frames
        time = t0 + frame * realtime_rate / fps
        u = sol(time)
        q = view(u, 1 : num_positions(mvis.state)) # TODO: make nicer
        atframe(animation, MeshCatMechanisms.visualizer(mvis), frame) do frame_vis
            set_configuration!(MechanismVisualizer(mvis.state, frame_vis), q)
        end
    end
    return animation
end

"""
Play back a visualization of a RigidBodySim.jl simulation.

Positional arguments:

* `vis` is a `MeshCatMechanisms.MechanismVisualizer`
* `sol` is a `DiffEqBase.ODESolution` obtained from a RigidBodySim.jl simulation.

`setanimation` accepts the following keyword arguments:

* `max_fps`: the maximum number of frames per second to draw. Default: `60`.
* `realtime_rate`: can be used to slow down or speed up playback compared to wall time. A `realtime_rate` of `2`
  will result in playback that is sped up 2x. Default: `1.0`.

# Examples

Visualizing the result of a simulation of the passive dynamics of an Acrobot (double pendulum) at half speed:

```jldoctest; output = false
using RigidBodySim, RigidBodyDynamics, MeshCatMechanisms, Blink
urdf = joinpath(dirname(pathof(RigidBodySim)), "..", "test", "urdf", "Acrobot.urdf")
mechanism = parse_urdf(Float64, urdf)
state = MechanismState(mechanism)
set_configuration!(state, [0.1; 0.2])
problem = ODEProblem(Dynamics(mechanism), state, (0., 2.))
sol = solve(problem, Vern7())
vis = MechanismVisualizer(mechanism, URDFVisuals(urdf))
# open(vis, Window()) # uncomment to open the visualizer window
setanimation!(vis, sol; realtime_rate = 0.5);

# output

```
"""
function MeshCat.setanimation!(mvis::MechanismVisualizer, sol::ODESolution;
        max_fps::Int = 60, realtime_rate::Number = 1., pause_pollint = nothing)
    # TODO: deprecate.
    if pause_pollint !== nothing
        warn("pause_pollint is no longer used. You can control the animation directly from the visualizer.")
    end
    @assert max_fps > 0
    @assert 0 < realtime_rate < Inf
    q0 = view(sol[1], 1 : num_positions(mvis.state)) # TODO: make nicer
    animation = MeshCat.Animation(mvis, sol; fps=max_fps, realtime_rate=realtime_rate)
    MeshCat.setanimation!(MeshCatMechanisms.visualizer(mvis), animation)
    set_configuration!(mvis.state, q0)
    nothing
end

end # module


text:
using FileType
using Test: @testset, @test, @test_throws

@testset "test_files_types" begin
    @test matcher("files/example.jpeg").extension == "jpg"
    @test matcher("files/example.gif").extension == "gif"
    @test matcher("files/example.webp").extension == "webp"
    @test matcher("files/example.tiff").extension == "tif"
    @test matcher("files/example.mp3").extension == "mp3"
    @test matcher("files/example.ttf").extension == "ttf"
    @test matcher("files/example.m4v").extension == "m4v"
    @test matcher("files/example.tar.gz").extension == "gz"
    @test matcher("files/example.wasm").extension == "wasm"
    @test matcher("files/example.dex").extension == "dex"
    @test matcher("files/example.dey").extension == "dey"
    @test Is(FileType.Images.Image,"files/example.png") == true
    @test is_extension_supported("png") == true
    @test is_mime_supported(MIME("image/jpeg")) == true
    @test_throws ErrorException Is(FileType.Images.Image,"files")
end


text:
"""
This file contains the code for the work by

Y. Chandak, S. Niekum, B. Castro da Silva, E. Learned-Miller, E. Brunskill, and P. S. Thomas.
Universal Off-Policy Evaluation.
In Advances in Neural Information Processing Systems, 2021

Code repo: https://github.com/yashchandak/UnO
"""

using Statistics
using Random
using Bootstrap
using LinearAlgebra

#######################################################################
# General Utilities
#######################################################################??

function decumsum(xs::Array{Float64, 1})
    """
    Create empirical probability distribution from empirical CDF
    """

    N = length(xs)
    temp = fill(0.0, N)
    temp[1] = xs[1]

    # First element stays as-is
    for i in 2:1:N
        temp[i] = xs[i] - xs[i-1]
    end

    return temp
end

function KS_distance(xs1, ys1, xs2, ys2)
    """
    Kolmogorov-Smirnov distance between 2 CDFs

    Inputs
    xs1, ys1, for the first CDF
    xs2, ys2, for the second CDF
    """

    function dist(x1, y1, x2, y2)
        # Check the distance between CDF1 and CDF2,
        # only at the points in CDF1

        n = length(x1)
        n2 = length(x2)
        j = 1
        max_d = -Inf
        for i in 1:1:n
            if j <= n2
                while x2[j] <= x1[i]
                    j = j+1
                    if j > n2; break; end # Prevent overflowing
                end
            end

            # Above while loop stops immediately when x2[j] crosses x1[i]
            # Therefore decrement j by one to get 2nd cdf value at x1[i]
            if j==1
                max_d = max(max_d, abs(y1[i] - 0))
            else
                max_d = max(max_d, abs(y1[i] - y2[j-1]))
            end
        end
        return max_d
    end

    # Search both ways: CDF1 vs CDF2, and CDF2 vs CDF1
    d1 = dist(xs1, ys1, xs2, ys2)
    d2 = dist(xs2, ys2, xs1, ys1)
    return max(d1, d2)
end

function get_equispaced_CDF(xs, ys, props, n_points=100)
    """
    Create a CDF at equispaced points from the sample CDF
    This function is particularly useful for plotting
    and averaging over multiple trials

    Note: This is only an approximation because of discretization
    """

    Gmin, Gmax = props["Gmin"], props["Gmax"]
    steps = (Gmax - Gmin)/(n_points - 1)

    # Initialize equi-spaced CDF
    Y_CDF = fill(0.0, n_points)
    X_points = fill(0.0, n_points)
    n = length(xs)

    j, k = 1, 1                     # Indexing starts at 1 in Julia
    for i in Gmin:steps:Gmax
        X_points[j] = i

        # Search for the smallest observed return greater than i
        if k <= n
            while xs[k] <= X_points[j]
                k = k+1
                if k > n; break; end # handle overflow
            end
        end

        if k == 1
            # If i is less than smallest observed return then set CDF = 0
            Y_CDF[j] = 0
        else
            # Otherwise CDF at i is the
            # estimated CDF of the largest observed return smaller than i
            Y_CDF[j] = ys[k-1]
        end

        j = j + 1
    end

    return X_points, Y_CDF
end

function get_equispaced_CDF2(xs, ys, Gmin, Gmax, n_points=100)
    """
    Create a CDF at equispaced points from the sample CDF
    This function is particularly useful for plotting
    and averaging over multiple trials

    Note: This is only an approximation because of discretization
    """

    steps = (Gmax - Gmin)/(n_points - 1)

    # Initialize equi-spaced CDF
    Y_CDF = fill(0.0, n_points)
    X_points = fill(0.0, n_points)
    n = length(xs)

    j, k = 1, 1                     # Indexing starts at 1 in Julia
    for i in Gmin:steps:Gmax
        X_points[j] = i

        # Search for the smallest observed return greater than i
        if k <= n
            while xs[k] <= X_points[j]
                k = k+1
                if k > n; break; end # handle overflow
            end
        end

        if k == 1
            # If i is less than smallest observed return then set CDF = 0
            Y_CDF[j] = 0
        else
            # Otherwise CDF at i is the
            # estimated CDF of the largest observed return smaller than i
            Y_CDF[j] = ys[k-1]
        end

        j = j + 1
    end

    return X_points, Y_CDF
end

function get_CDF_at_locs(xs, ys, X_points)
    """
    Get CDF values at some alternate points: X_points
    """

    n_points = length(X_points)
    Y_CDF = fill(0.0, n_points)
    n = length(xs)

    j, k = 1, 1                     # Indexing starts at 1 in Julia
    for i in X_points

        # Search for the smallest observed return greater than i
        if k <= n
            while xs[k] <= i
                k = k+1
                if k > n; break; end # handle overflow
            end
        end

        if k == 1
            # If i is less than smallest observed return then set CDF = 0
            Y_CDF[j] = 0
        else
            # Otherwise CDF at i is the
            # estimated CDF of the largest observed return smaller than i
            Y_CDF[j] = ys[k-1]
        end

        j = j + 1
    end

    return Y_CDF
end

function on_policy_CDF(pi_data)
    """
    Empirical CDF computed using on-policy Monte-carlo
    """

    N = length(pi_data)
    xs = sort(pi_data);
    ys =  collect(1.0/N:1.0/N:1.0)

    return xs, ys
end

function behavior_CDF(data, props)
    """
    Empirical CDF for the behavior policy
    """

    n, horizon, d = size(data)
    γ = props["gamma"]
    returns = fill(0.0, n)

    # Extract returns under behavior policy
    for i in 1:1:n
        G = 0
        for j in 1:1:horizon
            G += γ^(j-1) * data[i,j,2]
        end
        returns[i] = G
    end

    return on_policy_CDF(returns)
end

################################################
# Wild bootstrap stuff
# Percentile-based
#
# Y. Chandak, S. Jordan, G. Theocharous, M. White, and P. S. Thomas.
# Towards Safe Policy Improvement for Non-Stationary MDPs.
# In Advances in Neural Information Processing Systems, 2020.
#############################################

function get_coefs(Φ, ϕτ)
    """
    get_coefs(Φ, ϕτ)
    returns least squares coefficients for making predictions
    at observed points Φ and points of interest ϕτ.
    This function is a helper function to be used for predicting
    future points in a time series and for in wild bootstrap.
    """
    H = pinv(Φ' * Φ) * Φ'
    W = Φ * H
    ϕ = ϕτ * H
    return W, ϕ
end

function get_preds_and_residual(Y, W, ϕ)
    """
    get_preds_and_residual(Y, W, ϕ)
    returns the baseline predictions using observed labels Y at points ϕ
    and along with the vector of residuals Y - Ŷ.
    W and ϕ should be the output of get_coefs.
    """
    Ŷ = W*Y
    y = ϕ*Y
    ξ = Y .- Ŷ
    return y, ξ
end

function wildbs_eval(y, ξ, ϕ, σ, f)
        """
        wildbs_eval(y, ξ, ϕ, σ, f)
    evaluates the prediction of a linear timeseries model using
    noise generated for the wild bootstrap.
    y is the prediction at points ϕ using original labels Y
    ξ are the residual between the labels Y and prediction Ŷ
    ϕ are the points (coefficients) at at which the predictions y are made
    σ is the noise sample used to modify the residual, e.g., vector of {-1, 1}
    f is the function to aggregate the result for all points in ϕ, e.g., sum, mean, maximum, minimum
    """
    return f(y .+ ϕ * (ξ .* σ))
end

function wildbs_CI(y, ξ, ϕ, δ, num_boot, aggf, rng, mode=:normal)
    """
    wildbs_CI(y, ξ, ϕ, δ, num_boot, aggf)
    computes the δ percentiles bootstrap using the wild bootstrap method with num_boot bootstrap samples
    for original predictions y, with residuals ξ, at features ϕ, aggregated by aggf.
    """
    # bsamples = [wildbs_eval(y, ξ, ϕ, sign.(randn(rng, length(ξ))), aggf) for i in 1:num_boot]
    # bsamples = [wildbs_eval(y, ξ, ϕ, sign.(randn(length(ξ))), aggf) for i in 1:num_boot]        # Without random seed
    # return quantile(sort(bsamples), δ, sorted=true)
    function get_val()
        return wildbs_eval(y, ξ, ϕ, sign.(randn(length(ξ))), aggf)
    end

    n = length(ξ)
    idxs = collect(1:n)

    if mode == :fast
        bs = bootstrap(boot_id -> get_val(), idxs, BasicSampling(100))
        bci = confint(bs, PercentileConfInt(1 - 2*δ[1]))
        _, lb_param, ub_param = bci[1]
        return lb_param, ub_param
    end

    ##########
    # Use BCa
    ###############

    bs = bootstrap(boot_id -> get_val(), idxs, BasicSampling(num_boot))   # Bootstrap KS distance's distribution

    try
        # Important:
        # In the bootstrap library, due to some unknown bug a variable zalpha becomes NaN
        # Another way to avoid error is by using the following change to the BCa code IN THE BOOTSTRAP LIBRARY:
        # zalpha == zalpha || return (t0, t0, t0)  # Yash: avoid the edge case of zalpha being NaN
        bci = confint(bs, BCaConfInt(1 - 2*δ[1]))              # BCa Bootstrap CI
        _, lb_param, ub_param = bci[1]      # return format: estimate, lower, upper
        return lb_param, ub_param
    catch e
        println("BCa bootstrap error occured. Trying percentile bootstrap...")
        bci = confint(bs, PercentileConfInt(1 - 2*δ[1]))       # Percentile Bootstrap CI
        _, lb_param, ub_param = bci[1]      # return format: estimate, lower, upper
        return lb_param, ub_param
    end

end

function nswildbs_CI(idxs, Y, δ, L, mode,
                     tail=:both,
                     τ=1, num_boot=1000,
                     aggf=mean, rng=MersenneTwister(0))

    # fb = fourierseries(Float64, 8)
    fb = fourierseries(Float64, 0)
    nt = normalize_time(L, τ)
    ϕ(t) = fb(nt(t))

    Φ, ϕτ = create_features(ϕ, idxs, collect(Float64, L+1:L+τ))
    A, B = get_coefs(Φ, ϕτ)

    if tail == :left
        return wildbs_CI(get_preds_and_residual(Y, A, B)..., B, δ, num_boot, aggf, rng)
    elseif tail == :right
        return wildbs_CI(get_preds_and_residual(Y, A, B)..., B, 1-δ, num_boot, aggf, rng)
    elseif tail == :both
        return wildbs_CI(get_preds_and_residual(Y, A, B)..., B, [δ/2.0, 1-(δ/2.0)], num_boot, aggf, rng, mode)
    else
        println("ERROR tail: '$tail' not recognized. Returning NaN")
        return NaN
    end
end

#########################################
# Wild bootstrap-T stuff
# Algortithms 1,2, and 3 in the work by
#
# Y. Chandak, S. Jordan, G. Theocharous, M. White, and P. S. Thomas.
# Towards Safe Policy Improvement for Non-Stationary MDPs.
# In Advances in Neural Information Processing Systems, 2020.
########################################

function get_coefst(Φ, ϕτ)
    H = pinv(Φ' * Φ) * Φ'
    W = Φ * H
    ϕ = ϕτ * H
    C = I - W
    return W, ϕ, C
end

function get_preds_and_residual_t(Y, W, ϕ, C)
    Ŷ = W*Y
    y = ϕ * Ŷ
    x = C * Ŷ
    ξ = Y .- Ŷ
    Σ = ϕ * Diagonal(ξ.^2) * ϕ'
    v = mean(Σ)

    return y, v, x, ξ
end

function wildbst_eval(x, C, ξ, ϕ, σ)
    ξ̂ = ξ .* σ
    ξ̂2 = (x .+ C * ξ̂).^2
    Δŷ = mean(ϕ * ξ̂)
    Σ = ϕ * Diagonal(ξ̂2) * ϕ'
    v̂ = mean(Σ)
    return Δŷ / √v̂
end

function bst_CIs(p, v, bsamples, δ)
    return p - quantile(sort(bsamples), 1.0-δ, sorted=true) * √v
end

function bst_CIs(p, v, bsamples, δ::Array{T,1}) where {T}
    return p .- quantile(sort(bsamples), 1.0 .- δ, sorted=true) .* √v
end

function wildbst_CI(y, v, x, ξ, ϕ, C, δ, num_boot,rng)
    # bsamples = [wildbst_eval(x, C, ξ, ϕ, sign.(randn(rng, length(ξ)))) for i in 1:num_boot]
    bsamples = [wildbst_eval(x, C, ξ, ϕ, sign.(randn(length(ξ)))) for i in 1:num_boot] # Without RNG
    p = mean(y)
    # return p - quantile(sort(bsamples), 1.0-δ, sorted=true) * √v
    return bst_CIs(p, v, bsamples, δ)
end

function nswildbst_CI(idxs, Y, δ, L,
                        tail=:both,
                        τ=1, num_boot=1000,
                        aggf=mean, rng=MersenneTwister(0))

    fb = fourierseries(Float64, 4)
    nt = normalize_time(L, τ)
    ϕ(t) = fb(nt(t))

    Φ, ϕτ = create_features(ϕ, idxs, collect(Float64, L+1:L+τ))
    A, B, C = get_coefst(Φ, ϕτ)

    if tail == :left
        return wildbst_CI(get_preds_and_residual_t(Y, A, B, C)..., B, C, δ, num_boot, rng)
    elseif tail == :right
        return wildbst_CI(get_preds_and_residual_t(Y, A, B, C)..., B, C, 1-δ, num_boot, rng)
    elseif tail == :both
        return wildbst_CI(get_preds_and_residual_t(Y, A, B, C)..., B, C, [δ/2.0, 1-(δ/2.0)], num_boot, rng)
    else
        println("ERROR tail: '$tail' not recognized. Returning NaN")
        return NaN
    end
    return
end

"""
create features for time series using basis function ϕ,
for observed time points x, and future time points tau
"""
function create_features(ϕ, x, τ)
    Φ = hcat(ϕ.(x)...)'
    ϕτ = hcat(ϕ.(τ)...)'
    return Φ, ϕτ
end

function fourierseries(::Type{T}, order::Int) where {T}
    C = collect(T, 0:order) .* π
    return x -> @. cos(C * x)
end

function normalize_time(L, τ)
    return x -> x / (L + τ)
end

max_stars_repo_name: yashchandak/UnO

max_issues_repo_name: yashchandak/UnO

max_forks_repo_name: yashchandak/UnO


text:
"""
    BinaryConvRBM(visible_fields, hidden_fields, w)

Convolutional RBM with binary visible and hidden units,
with fields `a` and `b` and weights `w`.
"""
function BinaryConvRBM(
    visible_fields::AbstractArray,
    hidden_fields::AbstractArray,
    w::AbstractArray;
    kwargs...
)
    return ConvRBM(RBMs.Binary(visible_fields), RBMs.Binary(hidden_fields), w; kwargs...)
end

"""
    BinaryConvRBM(visible_size, hidden_size, kernel_size)
    BinaryConvRBM(T, visible_size, hidden_size, kernel_size)

Convolutional binary RBM with given dimensions, and parameters initialized to zero
of type `T` (= `Float64` by default).
"""
function BinaryConvRBM(
    ::Type{T},
    visible_size::Union{Int,TupleN{Int}},
    hidden_size::Union{Int,TupleN{Int}},
    kernel_size::Union{Int,TupleN{Int}};
    kwargs...
) where {T}
    visible_fields = zeros(T, visible_size...)
    hidden_fields = zeros(T, hidden_size...)
    w = zeros(T, visible_size..., kernel_size..., hidden_size...)
    return BinaryConvRBM(visible_fields, hidden_fields, w; kwargs...)
end

function BinaryConvRBM(
    visible_size::Union{Int,TupleN{Int}},
    hidden_size::Union{Int,TupleN{Int}},
    kernel_size::Union{Int,TupleN{Int}};
    kwargs...
)
    return BinaryConvRBM(Float64, visible_size, hidden_size, kernel_size; kwargs...)
end