text:
# A tuple of (Pix,Spin,Basis) defines the generic behavior of our fields
abstract type Pix end
abstract type Spin end
abstract type Basis end

# All fields are a subtype of this.
abstract type Field{B<:Basis, S<:Spin, P<:Pix, T} <: AbstractVector{T} end

# Spin types, "S0" is spin-0, i.e. a scalar map. "S2" is spin-2 like QU, and S02
# is a tuple of S0 and S2 like TQU.
abstract type S0 <: Spin end
abstract type S2 <: Spin end
abstract type S02 <: Spin end

# the basis type for a FieldTuple can be a BasisTuple which just holds the bases
# of the sub-fields
abstract type BasisTuple{T} <: Basis end
promote_type(::Type{BasisTuple{BT1}}, ::Type{BasisTuple{BT2}}) where {BT1,BT2} = BasisTuple{Tuple{map_tupleargs(promote_type,BT1,BT2)...}}

# Basis types
abstract type Map       <: Basis end
abstract type Fourier   <: Basis end
abstract type QUMap     <: Basis end
abstract type EBMap     <: Basis end
abstract type QUFourier <: Basis end
abstract type EBFourier <: Basis end
const IQUFourier = BasisTuple{Tuple{Fourier,QUFourier}}
const IEBFourier = BasisTuple{Tuple{Fourier,EBFourier}}
const IQUMap     = BasisTuple{Tuple{Map,QUMap}}
const IEBMap     = BasisTuple{Tuple{Map,EBMap}}

basis_promotion_rules = Dict(
    # S0
    (Map,       Fourier)    => Map,
    # S2
    (QUMap,     QUFourier)  => QUMap,
    (EBMap,     EBFourier)  => EBFourier,
    (QUMap,     EBMap)      => QUMap,
    (QUFourier, EBFourier)  => QUFourier,
    (QUMap,     EBFourier)  => QUMap,
    (QUFourier, EBMap)      => QUFourier
)
for ((B1,B2),B′) in basis_promotion_rules
    @eval promote_rule(::Type{$B1}, ::Type{$B2}) = $B′
end

# A "basis-like" object, e.g. the lensing basis Ł or derivative basis Ð. For any
# particular types of fields, these might be different actual bases, e.g. the
# lensing basis is Map for S0 but QUMap for S2.
abstract type Basislike <: Basis end
(::Type{B})(f::F) where {F<:Field,B<:Basislike} = B(F)(f)
(::Type{B})(f′::Field, f::F) where {F<:Field,B<:Basislike} = B(F)(f′,f)
(::Type{B})(a::AbstractArray{<:Field}...) where {B<:Basis} = B.(a...)

# Basis in which derivatives are sparse for a given field
abstract type DerivBasis <: Basislike end
const Ð = DerivBasis
Ð!(args...) = Ð(args...)

# Basis in which lensing is a pixel remapping for a given field
abstract type LenseBasis <: Basislike end
const Ł = LenseBasis
Ł!(args...) = Ł(args...)

# B(f) where B is a basis converts f to that basis. This is the fallback if the
# field is already in the right basis.
(::Type{B})(f::Field{B}) where {B} = f

# The abstract `Basis` type means "any basis", hence this conversion rule:
Basis(f::Field) = f

# B(f′, f) converts f to basis B and stores the result inplace in f′. If f is
# already in basis B, we just return f (but note, we never actually set f′ in
# this case, which is more efficient, but necessitates some care when using this
# construct)
(::Type{B})(f′::Field{B}, f::Field{B}) where {B} = f

### ImplicitOp
# ImplicitOps and ImplicitFields are operators and fields which may not be
# explicitly stored, but can be multiplied, etc... and generally implicitly take
# on whatever size they need. The specification of Float32 as the the Array type
# is only for inference, such that its widened to whatever its needed to be the
# end.
abstract type ImplicitOp{B<:Basis, S<:Spin, P<:Pix} <: AbstractMatrix{Float32} end
abstract type ImplicitField{B<:Basis, S<:Spin, P<:Pix} <: Field{B,S,P,Float32} end

# no size and 0-length represents the fact that ImplicitOps aquire the size of
# the fields they're applied to. it also helps makes the generic printing
# methods for AbstractMatrix work well.
size(::Union{ImplicitOp,ImplicitField}) = ()
size(::Union{ImplicitOp,ImplicitField},i) = nothing
length(::Union{ImplicitOp,ImplicitField}) = 0
checksquare(::ImplicitOp) = nothing

adapt_structure(to, x::Union{ImplicitOp,ImplicitField}) = x
adapt_structure(to, arr::Array{<:Field}) = map(f -> adapt(to,f), arr)

# printing
show(io::IO, ::MIME"text/plain", L::ImplicitOp) = show(io,L)
show(io::IO, ::MIME"text/plain", L::Adjoint{<:Any,<:ImplicitOp}) = show(io,L)
show(io::IO, L::Adjoint{<:Any,<:ImplicitOp}) = (print(io,"Adjoint{"); show(io,parent(L)); print(io,"}"))
# this is the main function ImplicitOps should specialize if this default behavior isn't enough:
show(io::IO, L::ImplicitOp) = showarg(io, L, true)

# All CMBLensing operators are then either Diagonals or ImplicitOps.
# ImplicitOrAdjOp are things for which algebra is done lazily. This used to
# include Diagonal{<:ImplicitField}, but that was leading so some really
# annoying ambiguities, so its removed for now.
const DiagOp{F<:Field,T} = Diagonal{T,F}
const LinOp{B,S,P} = Union{ImplicitOp{B,S,P},DiagOp{<:Field{B,S,P}}}
const ImplicitOrAdjOp{B,S,P} = Union{ImplicitOp{B,S,P}, Adjoint{<:Any,<:ImplicitOp{B,S,P}}}
const LinOrAdjOp{B,S,P} = Union{ImplicitOrAdjOp{B,S,P},DiagOp{<:Field{B,S,P}}}

# assume no dependence on parameters θ unless otherwise specified
(L::LinOrAdjOp)(θ::NamedTuple) = L

### Scalars
# scalars which are allowed in our expressions must be real because we
# implicitly assume our maps are real, and addition/multiplication by a complex
# number, even of the fourier transform, would break this.
const Scalar = Real
const FieldOrOp = Union{Field,LinOp}
const FieldOpScal = Union{Field,LinOp,Scalar}

"""
    logdet(L::LinOp, θ)

If L depends on θ, evaluates `logdet(L(θ))` offset by its fiducial value at
`L()`. Otherwise, returns 0.
"""
logdet(L::LinOp, θ) = depends_on(L,θ) ? logdet(L()\L(θ)) : 0
logdet(L::Int, θ) = 0 # the default returns Float64 which unwantedly poisons the backprop to Float64
logdet(L, θ) = logdet(L)

### Simulation

@doc doc"""
    simulate(Σ; rng=global_rng_for(Σ), seed=nothing)

Draw a simulation from the covariance matrix `Σ`, i.e. draw a random vector
$\xi$ such that the covariance $\langle \xi \xi^\dagger \rangle = \Sigma$.

The random number generator `rng` will be used and advanced in the proccess, and
is by default the appropriate one depending on if `Σ` is backed by `Array` or
`CuArray`.

The `seed` argument can also be used to seed the `rng`.
"""
function simulate(Σ; rng=global_rng_for(Σ), seed=nothing)
    (seed != nothing) && Random.seed!(rng, seed)
    simulate(rng, Σ)
end
function white_noise(Σ; rng=global_rng_for(Σ), seed=nothing)
    (seed != nothing) && Random.seed!(rng, seed)
    white_noise(rng, Σ)
end
function fixed_white_noise(Σ; rng=global_rng_for(Σ), seed=nothing)
    (seed != nothing) && Random.seed!(rng, seed)
    fixed_white_noise(rng, Σ)
end
global_rng_for(x::T) where {T<:AbstractArray} = global_rng_for(T)
global_rng_for(::Type{<:Array}) = Random.GLOBAL_RNG

#
# ### Matrix conversion
#
# # We can build explicit matrix representations of linear operators by applying
# # them to a set of basis vectors which form a complete basis (in practice this
# # is prohibitive except for fairly small maps, e.g. 16x16 pixels, but is quite
# # useful)
#
# """
#     full(::Type{Fin}, ::Type{Fout}, L::LinOp; progress=true)
#     full(::Type{F}, L::LinOp; progress=true)
#
# Construct an explicit matrix representation of the linear operator `L` by
# applying it to a set of vectors which form a complete basis. The `Fin` and
# `Fout` types should be fields which specify the input and output bases for the
# representation (or just `F` if `L` is square and we want the same input/output
# bases)
#
# The name `full` is to be consistent with Julia's SparseArrays where `full`
# builds a full matrix from a sparse one.
# """
# full(::Type{Fin}, ::Type{Fout}, L::LinOp; progress=true) where {Fin<:Field, Fout<:Field} =
#     hcat(@showprogress(progress ? 1 : Inf, [(Fout(L*(x=zeros(length(Fin)); x[i]=1; x)[Tuple{Fin}]))[:] for i=1:length(Fin)])...)
# full(::Type{F}, L::LinOp; kwargs...) where {F<:Field} = full(F,F,L; kwargs...)
#
#

### Other generic stuff

# allows one to pass `1` to something which expect a LenseFlow-like operator
(s::Scalar)(::Field) = s

alloc_cache(x, ::Any) = x
cache(x, ::Any) = x
cache!(x, ::Any) = x

# caching for adjoints
cache(L::Adjoint, f) = cache(L',f)'
cache!(L::Adjoint, f) = cache!(L',f)'

# convenience "getter" functions for the Basis/Spin/Pix
# basis of UnionAlls like basis(Field) will return Basis (which means any Basis)
basis(::Type{<:Field}) = Basis
basis(::Type{<:Field{B,S,P}}) where {B<:Basis,S<:Spin,P<:Pix} = B
basis(::F) where {F<:Field} = basis(F)
spin(::Type{<:Field}) = Spin
spin(::Type{<:Field{B,S,P}}) where {B<:Basis,S<:Spin,P<:Pix} = S
spin(::F) where {F<:Field} = spin(F)
pix(::Type{<:Field}) = Pix
pix(::Type{<:Field{B,S,P}}) where {B<:Basis,S<:Spin,P<:Pix} = P
pix(::F) where {F<:Field} = pix(F)

# we use Field cat'ing mainly for plotting, e.g. plot([f f; f f]) plots a 2×2
# matrix of maps. the following definitions make it so that Fields aren't
# splatted into a giant matrix when doing [f f; f f] (which they would othewise
# be since they're Arrays)
hvcat(rows::Tuple{Vararg{Int}}, values::Field...) = hvcat(rows, ([x] for x in values)...)
hcat(values::Field...) = hcat(([x] for x in values)...)

### printing
print_array(io::IO, f::Field) = print_array(io, f[:])
show_vector(io::IO, f::Field) = show_vector(io, f[:])
show_vector(io::IO, f::ImplicitField) = print(io, "[…]")

# addition/subtraction works between any fields and scalars, promotion is done
# automatically if fields are in different bases
for op in (:+,:-), (T1,T2,promote) in ((:Field,:Scalar,false),(:Scalar,:Field,false),(:Field,:Field,true))
    @eval ($op)(a::$T1, b::$T2) = broadcast($op, ($promote ? promote(a,b) : (a,b))...)
end

≈(a::Field, b::Field) = ≈(promote(a,b)...)

# multiplication/division is not strictly defined for abstract vectors, but
# make it work anyway if the two fields are exactly the same type, in which case
# its clear we wanted broadcasted multiplication/division.
for op in (:*, :/)
    @eval ($op)(a::Field{B,S,P}, b::Field{B,S,P}) where {B,S,P} = broadcast($op, a, b)
end

# needed unless I convince them to undo the changes here:
# https://github.com/JuliaLang/julia/pull/35257#issuecomment-657901503
if VERSION>v"1.4"
    *(x::Adjoint{<:Number,<:Field}, y::Field) = dot(x.parent,y)
end

(::Type{T})(f::Field{<:Any,<:Any,<:Any,<:Real}) where {T<:Real} = T.(f)
(::Type{T})(f::Field{<:Any,<:Any,<:Any,<:Complex}) where {T<:Real} = Complex{T}.(f)

# misc
one(f::Field) = fill!(similar(f), one(eltype(f)))
norm(f::Field) = sqrt(dot(f,f))
sum_kbn(f::Field) = sum_kbn(f[:])

function invalid_broadcast_error(B1,F1,B2,F2)
    if B1!=B2
        error("""Can't broadcast across fields in $B1 and $B2 bases.
        Try the same operation without broadcasting (which will do an automatic basis conversion).""")
    else
        error("""Broadcasting across fields with the following differing broadcast styles is not implemented:
        * $F1
        * $F2
        """)
    end
end

@init @require PyCall="438e738f-606a-5dbb-bf0a-cddfbfd45ab0" begin
    # never try to auto-convert Fields or LinOps to Python arrays
    PyCall.PyObject(x::Union{LinOp,Field}) = PyCall.pyjlwrap_new(x)
end

gc = () -> GC.gc(true)


text:
using BinDeps

@BinDeps.setup

mpi = library_dependency("juliampi", aliases=["libjuliampi", "libjuliampi.so", "libjuliampi.dylib"])

builddir = joinpath(BinDeps.depsdir(mpi),"build")

provides(SimpleBuild,
    (@build_steps begin
        CreateDirectory(builddir)
        (@build_steps begin
            ChangeDirectory(builddir)
            (@build_steps begin
                `cmake -DCMAKE_INSTALL_PREFIX=$(pwd())/../src -DCMAKE_LIB_INSTALL_PREFIX=$(pwd())/../usr/lib ..`
                `make`
                `make install`
                `make test`
            end)
        end)
    end), [mpi], os = :Darwin)

provides(SimpleBuild,
    (@build_steps begin
        CreateDirectory(builddir)
        (@build_steps begin
            ChangeDirectory(builddir)
            (@build_steps begin
                `cmake -DCMAKE_INSTALL_PREFIX=$(pwd())/../src -DCMAKE_LIB_INSTALL_PREFIX=$(pwd())/../usr/lib ..`
                `make`
                `make install`
                `make test`
            end)
        end)
    end), [mpi], os = :Unix)

@BinDeps.install [:mpi => :mpi]

ext: jl


text:
# Li, Lin, "Accelerated Proximal Gradient Methods for Nonconvex Programming",
# Proceedings of NIPS 2015 (2015).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

struct LiLin_iterable{R <: Real, C <: Union{R, Complex{R}}, Tx <: AbstractArray{C}, Tf, TA, Tg}
    f::Tf             # smooth term
    A::TA             # matrix/linear operator
    g::Tg             # (possibly) nonsmooth, proximable term
    x0::Tx            # initial point
    gamma::Maybe{R}   # stepsize parameter of forward and backward steps
    adaptive::Bool    # enforce adaptive stepsize even if L is provided
    delta::R          #
    eta::R            #
end

mutable struct LiLin_state{R <: Real, Tx, TAx}
    x::Tx             # iterate
    y::Tx             # extrapolated point
    Ay::TAx           # A times y
    f_Ay::R           # value of smooth term at y
    grad_f_Ay::TAx    # gradient of f at Ay
    At_grad_f_Ay::Tx  # gradient of smooth term at y
    # TODO: *two* gammas should be used in general, one for y and one for x
    gamma::R          # stepsize parameter of forward and backward steps
    y_forward::Tx     # forward point at y
    z::Tx             # forward-backward point
    g_z::R            # value of nonsmooth term at z
    res::Tx           # fixed-point-residual (at y)
    theta::R          # auxiliary sequence to compute extrapolated points
    F_average::R      # moving average of objective values
    q::R              # auxiliary sequence to compute moving average
end

function Base.iterate(iter::LiLin_iterable{R}) where R
    y = iter.x0
    Ay = iter.A * y
    grad_f_Ay, f_Ay = gradient(iter.f, Ay)

    # TODO: initialize gamma if not provided
    # TODO: authors suggest Barzilai-Borwein rule?
    # TODO: *two* gammas should be used in general, one for y and one for x

    # compute initial forward-backward step
    At_grad_f_Ay = iter.A' * grad_f_Ay
    y_forward = y - iter.gamma .* At_grad_f_Ay
    z, g_z = prox(iter.g, y_forward, iter.gamma)

    Fy = f_Ay + iter.g(y)

    @assert isfinite(Fy) "initial point must be feasible"

    # compute initial fixed-point residual
    res = y - z

    state = LiLin_state(
        copy(iter.x0), y, Ay, f_Ay, grad_f_Ay, At_grad_f_Ay,
        iter.gamma, y_forward, z, g_z, res, R(1), Fy, R(1)
    )

    return state, state
end

function Base.iterate(iter::LiLin_iterable{R}, state::LiLin_state{R, Tx, TAx}) where {R, Tx, TAx}
    # TODO: backtrack gamma at y

    Fz = iter.f(state.z) + state.g_z

    theta1 = (R(1)+sqrt(R(1)+4*state.theta^2))/R(2)

    if Fz <= state.F_average - iter.delta * norm(state.res)^2
        case = 1
    else
        # TODO: re-use available space in state?
        # TODO: backtrack gamma at x
        Ax = iter.A * state.x
        grad_f_Ax, f_Ax = gradient(iter.f, Ax)
        At_grad_f_Ax = iter.A' * grad_f_Ax
        x_forward = state.x - state.gamma .* At_grad_f_Ax
        v, g_v = prox(iter.g, x_forward, state.gamma)
        Fv = iter.f(v) + g_v
        case = Fz <= Fv ? 1 : 2
    end

    if case == 1
        state.y .= state.z .+ ((state.theta - R(1)) / theta1) .* (state.z .- state.x)
        state.x, state.z = state.z, state.x
        Fx = Fz
    elseif case == 2
        state.y .= state.z .+ (state.theta / theta1) .* (state.z .- v) .+ ((state.theta - R(1)) / theta1) .* (v .- state.x)
        state.x = v
        Fx = Fv
    end

    mul!(state.Ay, iter.A, state.y)
    state.f_Ay = gradient!(state.grad_f_Ay, iter.f, state.Ay)
    mul!(state.At_grad_f_Ay, adjoint(iter.A), state.grad_f_Ay)
    state.y_forward .= state.y .- state.gamma .* state.At_grad_f_Ay
    state.g_z = prox!(state.z, iter.g, state.y_forward, state.gamma)

    state.res .= state.y - state.z

    state.theta = theta1

    # NOTE: the following can be simplified
    q1 = iter.eta * state.q + 1
    state.F_average = (iter.eta * state.q * state.F_average + Fx)/q1
    state.q = q1

    return state, state
end

# Solver

struct LiLin{R <: Real}
    gamma::Maybe{R}
    adaptive::Bool
    delta::R
    eta::R
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int

    function LiLin{R}(; gamma::Maybe{R}=nothing, adaptive::Bool=false,
        delta::R=R(1e-3), eta::R=R(0.8), maxit::Int=10000, tol::R=R(1e-8),
        verbose::Bool=false, freq::Int=100
    ) where R
        @assert gamma === nothing || gamma > 0
        @assert delta > 0
        @assert 0 < eta < 1
        @assert maxit > 0
        @assert tol > 0
        @assert freq > 0
        new(gamma, adaptive, delta, eta, maxit, tol, verbose, freq)
    end
end

function (solver::LiLin{R})(
    x0::AbstractArray{C}; f=Zero(), A=I, g=Zero(), L::Maybe{R}=nothing
) where {R, C <: Union{R, Complex{R}}}

    stop(state::LiLin_state) = norm(state.res, Inf)/state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e\n",
        it, state.gamma, norm(state.res, Inf)/state.gamma
    )

    if solver.gamma === nothing && L !== nothing
        gamma = R(1)/L
    elseif solver.gamma !== nothing
        gamma = solver.gamma
    else
        gamma = nothing
    end

    iter = LiLin_iterable(f, A, g, x0, gamma, solver.adaptive, solver.delta, solver.eta)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose iter = tee(sample(iter, solver.freq), disp) end

    num_iters, state_final = loop(iter)

    return state_final.z, num_iters

end

# Outer constructors

"""
    LiLin([gamma, adaptive, fast, maxit, tol, verbose, freq])

Instantiate the nonconvex accelerated proximal gradient method by Li and Lin
(see Algorithm 2 in [1]) for solving optimization problems of the form

    minimize f(Ax) + g(x),

where `f` is smooth and `A` is a linear mapping (for example, a matrix).
If `solver = LiLin(args...)`, then the above problem is solved with

    solver(x0, [f, A, g, L])

Optional keyword arguments:

* `gamma::Real` (default: `nothing`), the stepsize to use; defaults to `1/L` if not set (but `L` is).
* `adaptive::Bool` (default: `false`), if true, forces the method stepsize to be adaptively adjusted.
* `delta::Real` (default: `1e-3`), parameter determinining when extrapolated steps are to be accepted.
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `tol::Real` (default: `1e-8`), absolute tolerance on the fixed-point residual.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `10`), frequency of verbosity.

If `gamma` is not specified at construction time, the following keyword
argument can be used to set the stepsize parameter:

* `L::Real` (default: `nothing`), the Lipschitz constant of the gradient of x ↦ f(Ax).

References:

[1] Li, Lin, "Accelerated Proximal Gradient Methods for Nonconvex Programming",
Proceedings of NIPS 2015 (2015).
"""
LiLin(::Type{R}; kwargs...) where R = LiLin{R}(; kwargs...)
LiLin(; kwargs...) = LiLin(Float64; kwargs...)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@recipe function f(samples::DensitySampleVector;
                params=collect(1:5),
                mean=false,
                std_dev=false,
                globalmode=false,
                localmode=false,
                diagonal = Dict(),
                upper = Dict(),
                lower = Dict(),
                param_labels = [])

    parsel = params
    parsel = parsel[parsel .<= Base.size(samples.params[1], 1)]

    if Base.size(param_labels, 1) == 0
        param_labels = [latexstring("\\theta_$i") for i in parsel]
        param_labels_y = [latexstring("p(\\theta_$i)") for i in parsel]
    else
        param_labels_y = [latexstring("p("*param_labels[i]*")") for i in 1:length(param_labels)]
        param_labels = [latexstring(param_labels[i]) for i in 1:length(param_labels)]
    end

    nparams = length(parsel)
    layout --> nparams^2
    size --> (1000, 600)

    for i in 1:nparams

        # diagonal
        @series begin
            subplot := i + (i-1)*nparams

            seriestype --> get(diagonal, "seriestype", :smallest_intervals)
            bins --> get(diagonal, "bins", 200)
            colors --> get(diagonal, "colors", standard_colors)
            intervals --> get(diagonal, "intervals", standard_confidence_vals)
            legend --> get(diagonal, "legend", false)
            mean --> get(diagonal, "mean", mean)
            std_dev --> get(diagonal, "std_dev", std_dev)
            globalmode --> get(diagonal, "globalmode", globalmode)
            localmode --> get(diagonal, "localmode", localmode)
            xlims --> get(diagonal, "xlims", :auto)
            ylims --> get(diagonal, "ylims", :auto)
            xguide --> param_labels[i]
            yguide --> param_labels_y[i]

            samples, (parsel[i])
        end

        # upper right plots
        for j in i+1:nparams

            @series begin
                subplot := j + (i-1)*nparams

                seriestype --> get(upper, "seriestype", :histogram)
                bins --> get(upper, "bins", 200)
                colors --> get(upper, "colors", standard_colors)
                intervals --> get(upper, "intervals", standard_confidence_vals)
                legend --> get(upper, "legend", false)
                mean --> get(upper, "mean", mean)
                std_dev --> get(upper, "std_dev", std_dev)
                globalmode --> get(upper, "globalmode", globalmode)
                localmode --> get(upper, "localmode", localmode)
                xlims --> get(upper, "xlims", :auto)
                ylims --> get(upper, "ylims", :auto)
                xguide --> param_labels[i]
                yguide --> param_labels[j]

                samples, (parsel[i], parsel[j])
            end

            # lower left plots
            @series begin
                subplot := i + (j-1)*nparams

                seriestype --> get(lower, "seriestype", :smallest_intervals)
                bins --> get(lower, "bins", 200)
                colors --> get(lower, "colors", standard_colors)
                intervals --> get(lower, "intervals", standard_confidence_vals)
                legend --> get(lower, "legend", false)
                mean --> get(lower, "mean", mean)
                std_dev --> get(lower, "std_dev", std_dev)
                globalmode --> get(lower, "globalmode", globalmode)
                localmode --> get(lower, "localmode", localmode)
                xlims --> get(lower, "xlims", :auto)
                ylims --> get(lower, "ylims", :auto)
                xguide --> param_labels[i]
                yguide --> param_labels[j]

                samples, (parsel[i], parsel[j])
            end

        end
    end

end


text:
using FileIO
using RobotOSData

# Load all of the sample bags:
sample_bags = [joinpath("data", name) for name in readdir("data") if endswith(name, ".bag")]
for bagfile in sample_bags
    bag = load(bagfile)
end

# Load a bag with a custom type:
include("test_types.jl")
bag = load("data/convergent_gen1.bag", TestTypes)
msgs = bag["convergent"]
@assert msgs isa Vector{MessageData{TestTypes.test_rosbag.Convergent}}

# TODO:
# - test message filters
# - test for message content


text:
@testset "593.valid-square.jl" begin
    @test valid_square([[0, 0], [1, 1], [1, 0], [0, 1]]) == true
    @test valid_square([[0, 0], [1, 1], [1, 0], [0, -1]]) == false
    @test valid_square([[0, 1], [2, 0], [3, 2], [1, 3]]) == true
end


text:
function add_modifier(obj::MeasurementBuilder, arg0::EstimationModifier)
    return jcall(obj, "addModifier", void, (EstimationModifier,), arg0)
end

function get_modifiers(obj::MeasurementBuilder)
    return jcall(obj, "getModifiers", List, ())
end

function init(obj::MeasurementBuilder, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end


text:
# this script basially only handles `Base.ARGS`

Base.push!(LOAD_PATH, joinpath(@__DIR__, "..", "packages"))
using VSCodeServer
pop!(LOAD_PATH)

# load Revise ?
if "USE_REVISE" in Base.ARGS
    try
        @eval using Revise
        Revise.async_steal_repl_backend()
    catch err
        @warn "failed to load Revise: $err"
    end
end

atreplinit() do repl
    "USE_PLOTPANE" in Base.ARGS && Base.Multimedia.pushdisplay(VSCodeServer.InlineDisplay())
end

let
    conn_pipeline, telemetry_pipeline = Base.ARGS[1:2]
    VSCodeServer.serve(conn_pipeline; is_dev = "DEBUG_MODE" in Base.ARGS, crashreporting_pipename = telemetry_pipeline)
end


text:
# src/rounding.jl
# Generic Rounding Protocols

"""```
_round_rational_internal(::Type{T}, x::AbstractRational{Tr}, inner_term)
```
Internal implementation of the rational rounding.
Adapted from Julia internals.
"""
function _round_rational_internal(::Type{T}, x::AbstractRational{Tr}, inner_term) where {T,Tr}
    if denominator(x) == zero(Tr) && T <: Integer
        throw(DivideError())
    elseif denominator(x) == zero(Tr)
        return convert(T, copysign(one(Tr) // zero(Tr), numerator(x)))
    end
    q, r = divrem(numden(x)...)
    s = q
    if abs(r) >= abs((denominator(x) - copysign(Tr(4), numerator(x)) + one(Tr) +
                      inner_term(x, q)) >> 1 + copysign(Tr(2), numerator(x)))
        s += copysign(one(Tr), numerator(x))
    end
    return convert(T, s)
end

"""```
_round_rational(::Type{T}, x::AbstractRational{Tr}, ::RoundingMode{:Nearest})
```
Internal implementation of the `Nearest` mode rational rounding.
Adapted from Julia internals.
"""
function _round_rational(
    ::Type{T},
    x::AbstractRational{Tr},
    ::RoundingMode{:Nearest},
) where {T,Tr}
    return _round_rational_internal(T, x, (_, q) -> iseven(q))
end

"""```
_round_rational(::Type{T}, x::AbstractRational{Tr}, ::RoundingMode{:NearestTiesAway})
```
Internal implementation of the `NearestTiesAway` mode rational rounding.
Adapted from Julia internals.
"""
function _round_rational(
    ::Type{T},
    x::AbstractRational{Tr},
    ::RoundingMode{:NearestTiesAway},
) where {T,Tr}
    return _round_rational_internal(T, x, (_, q) -> zero(q))
end

"""```
_round_rational(::Type{T}, x::AbstractRational{Tr}, ::RoundingMode{:NearestTiesUp})
```
Internal implementation of the `NearestTiesUp` mode rational rounding.
Adapted from Julia internals.
"""
function _round_rational(
    ::Type{T},
    x::AbstractRational{Tr},
    ::RoundingMode{:NearestTiesUp},
) where {T,Tr}
    return _round_rational_internal(T, x, (x, q) -> numerator(x) < zero(q))
end


text:
import Base.Docs: Binding, @var
import Base.Markdown: MD, Code, Paragraph

# flat_content(md) = md
# flat_content(xs::Vector) = reduce((xs, x) -> vcat(xs,flat_content(x)), [], xs)
# flat_content(md::MD) = flat_content(md.content)
flatten(md::MD) = MD(flat_content(md))

# Faster version

function flat_content!(xs, out = [])
  for x in xs
    if isa(x, MD)
      flat_content!(x.content, out)
    else
      push!(out, x)
    end
  end
  return out
end

flat_content(md::MD) = flat_content!(md.content)

function hasdoc(b::Binding)
    for m in Docs.modules
        meta = Docs.meta(m)
        if haskey(meta, b) || haskey(meta, b[])
            return true
        end
    end
    false
end

hasdoc(b) = hasdoc(Docs.aliasof(b, typeof(b)))

function fullsignature(b::Binding)
  hasdoc(b) || return
  docs = Docs.doc(b)
  isa(docs, MD) || return
  first = flatten(docs).content[1]
  code =
    isa(first, Code) ? first.code :
    isa(first, Paragraph) && isa(first.content[1], Code) ?
      first.content[1].code :
      ""
  if Base.startswith(code, string(b.var))
    split(code, "\n")[1]
  end
end

function signature(b::Binding)
  sig = fullsignature(b)
  sig == nothing && return
  replace(sig, r" -> .*$", "")
end

function returns(b::Binding)
  r = r" -> (.*)"
  sig = fullsignature(b)
  sig == nothing && return
  if ismatch(r, sig)
    ret = match(r, sig).captures[1]
    if length(ret) < 10
      ret
    end
  end
end

function description(b::Binding)
  hasdoc(b) || return
  docs = Docs.doc(b)
  isa(docs, MD) || return
  md = flatten(docs)
  first = md.content[1]
  if isa(first, Code)
    length(md.content) < 2 && return
    first = md.content[2]
  end
  if isa(first, Paragraph)
    desc = Markdown.plain(first)
    return length(desc) > 100 ?
      desc[1:100]*"..." :
      desc
  end
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 37dc27ba-f78a-11ea-1d7d-2569658183fa
md"# _Introduction_

#### Who shoud read this tutorial?
Anyone who is interested in knowing what Quantum Computing is about or is new to Quantum Computing.

#### Requirements:
##### What should I know?
Knowledge of high school Mathematics should be enough.
Also, this course is in Julia, so knowing basic syntax of Julia might 			  help... but knowing Julia syntax is not necessary. Even if you know the 				basics of Python, R or Ruby, and how to write programs in them, you 				  should be fine.
##### The required software :-
- [Julia] (https://julialang.org/downloads/)
- [Yao.jl] (https://yaoquantum.org/)
Install Julia on your system if you haven't already. This tutorial assumes you already have Julia and you opened this notebook in your system."

# ╔═╡ 51da098a-071f-11eb-368f-9171f0d5a43e
md"As you've opened this tutorial, all the required packages including Yao are being installed currently. It requires an active internet connection. Don't close Julia or this notebook before all the required packages are installed. You can check whether the required packages have been installed or not. If it displays \"Installed\" before the next heading, which is Quantum, it means the required packages have been installed. If it's blank, then wait until you see it."

# ╔═╡ e5004b1c-071e-11eb-1846-5991805d5c2c
begin
	import Pkg
	Pkg.add("Yao")
	Pkg.add("YaoExtensions")
	Pkg.add("YaoPlots")
	Pkg.add("StatsBase")
	Pkg.add("Makie")
	Pkg.add("SymEngine")
	md"Installed"
end

# ╔═╡ b0a72718-f78e-11ea-1ef3-01bd13cfcd17
md"## _*Quantum*_

Now recently, if you're in any way affiliated to Science or Computer Science or even if you're just a tech lover, you must have stumbled across the word \"_quantum_\".
You must have wondered, \" What does it mean?\" .

So, let's just say that historically, some experiments were conducted by some scientists in Physics, and their results happened to defy what laws of Physics stated, before then. Now this was around the year 1900. The scientists came up with theories that explained those results and the most widely accepted (and working) theory was quantum mechanics, or in Einstein's words, \"*Real Black Magic Calculus*\". On the basis of quantum mechanics, we now have Quantum Physics, Quantum Chemistry, Quantum Biology, Quantum Computers and lots of other new stuff. The quantum computers just exploit some weird phenomenons of quantum mechanics.

Quantum Mechanics will not be explained in detail, in this tutorial. Only the part required for computation. If you're interested in Quantum Mechanics, [Modern Quantum Mechanics by J. J. Sakurai] (https://doi.org/10.1017/9781108499996) might be a good read."

# ╔═╡ 8c5d75be-f794-11ea-2b14-d108a10dc9d9
md"## *Bits*
Okay, what makes quantum computers different than a present day or _classical_ one.
To understand that, we must understand a few things about classical computers.

\"*A bit (short for binary digit) is the smallest unit of data in a computer*\". You must have read this in a book or article about computers. Computers use bits, represented by the digits 0 and 1, to store data. We organize these bits to store information and manipulate these bits and perform operations on them to get a variety of things done.
Like storing numbers, or images or videos. Adding numbers, subtracting numbers.
How? Well... Consider storing a number. In this case, the collection of bits 1 0 1 1 can be considered as, $ 2^3 + 0 + 2^1 + 2^0 = 11 $. The genereal idea being, the summation of $ 2^{position of 1s(starting from 0) from right} $.

Then just divide your computer screen into a matrix of 1000s of cells, and every cell containing a number, corresponding to the colour in that cell. Yeah, images are stored that way, those numbers are called pixels. Addition looks a bit more complicated than this. It's done using something called 'gates'. Gates take one or more bits and perform a logical operation on them to give a certain output. There's an \"AND gate\" which, takes in 2 bits and multiplies them to give an output. There're more gates like the OR gate, XOR gate etc. These gates are arranged in a certain manner to perform feats like addition, subtraction, etc."

# ╔═╡ Cell order:
# ╟─37dc27ba-f78a-11ea-1d7d-2569658183fa
# ╟─51da098a-071f-11eb-368f-9171f0d5a43e
# ╟─e5004b1c-071e-11eb-1846-5991805d5c2c
# ╟─b0a72718-f78e-11ea-1ef3-01bd13cfcd17
# ╟─8c5d75be-f794-11ea-2b14-d108a10dc9d9

max_stars_repo_name: logankilpatrick/tutorials

max_issues_repo_name: logankilpatrick/tutorials

max_forks_repo_name: logankilpatrick/tutorials


text:
export AbstractApproximator,
    ApproximatorStyle, Q_APPROXIMATOR, QApproximator, V_APPROXIMATOR, VApproximator

"""
    (app::AbstractApproximator)(env)

An approximator is a functional object for value estimation.
It serves as a black box to provides an abstraction over different
kinds of approximate methods (for example DNN provided by Flux or Knet).
"""
abstract type AbstractApproximator end

"""
    update!(a::AbstractApproximator, correction)

Usually the `correction` is the gradient of inner parameters.
"""
function RLBase.update!(a::AbstractApproximator, correction) end

#####
# traits
#####

abstract type AbstractApproximatorStyle end

"""
Used to detect what an [`AbstractApproximator`](@ref) is approximating.
"""
function ApproximatorStyle(::AbstractApproximator) end

struct QApproximator <: AbstractApproximatorStyle end

const Q_APPROXIMATOR = QApproximator()

struct VApproximator <: AbstractApproximatorStyle end

const V_APPROXIMATOR = VApproximator()


text:
using Base.Test
using SmoothingKernels

for testfile in ["kernels.jl"]
    include(testfile)
end


text:
# Goal of this file is to test unicode support across platforms we'll be using (Atom, Terminal, Jupyter) and across OSs.

# Glyphs
    L̃ = 2.0 # L-tilde
    u₀ = 3.0 # u-naught
    â = 4.0 # a-hat
    z̲ = 5.0 # z-underbar

# Ordinary unicode characters
    println("π, α, ζ, ξ, δ, ρ, γ, β, χ, ω, θ, υ, ν, Γ, Δ")

# Git
    "⋅" # cdot
    "⊗" # otimes
    "∂" # partial
    "ℒ" # scrL
    "𝒟" # scrD
    "∈" # in
    "∉" # notin
    "⊆" # subseteq
    "≠" # ne
    "≈" # approx

# Ones to avoid!
# Combining characters don't show on Windows Atom/vscode/jupyter
    γ̃  = 1 # gamma-tilde
    α̂ = 1 # alpha-hat
    π̲ = 1 # pi-underbar
    π̄ = 1 # pi-bar

max_stars_repo_name: Shunsuke-Hori/lecture-source-jl

max_issues_repo_name: Shunsuke-Hori/lecture-source-jl

max_forks_repo_name: Shunsuke-Hori/lecture-source-jl


text:
using LinearLeastSquares
using Gadfly

temps = readdlm("melbourne_temps.txt", ',')
n = size(temps)[1]
p = plot(
  x=1:1500, y=temps[1:1500], Geom.line,
  Theme(panel_fill=color("white"))
)
# draw(PNG("melbourne.png", 16cm, 12cm), p)

yearly = Variable(n)
eq_constraints = []
for i in 365 + 1 : n
  eq_constraints += yearly[i] == yearly[i - 365]
end

smoothing = 100
smooth_objective = sum_squares(yearly[1 : n - 1] - yearly[2 : n])
optval = minimize!(sum_squares(temps - yearly) + smoothing * smooth_objective, eq_constraints)
residuals = temps - evaluate(yearly)

# Plot smooth fit
p = plot(
  layer(x=1:1500, y=evaluate(yearly)[1:1500], Geom.line, Theme(default_color=color("red"), line_width=2px)),
  layer(x=1:1500, y=temps[1:1500], Geom.line),
  Theme(panel_fill=color("white"))
)
# draw(PNG("yearly_fit.png", 16cm, 12cm), p)

# Plot residuals for a few days
p = plot(
  x=1:100, y=residuals[1:100], Geom.line,
  Theme(default_color=color("green"), panel_fill=color("white"))
)
# draw(PNG("residuals.png", 16cm, 12cm), p)

# Generate the residuals matrix
ar_len = 5
residuals_mat = residuals[ar_len : n - 1]
for i = 1:ar_len - 1
  residuals_mat = [residuals_mat residuals[ar_len - i : n - i - 1]]
end

# Solve autoregressive problem
ar_coef = Variable(ar_len)
optval2 = minimize!(sum_squares(residuals_mat * ar_coef - residuals[ar_len + 1 : end]))

# plot autoregressive fit of daily fluctuations for a few days
ar_range = 1:145
day_range = ar_range + ar_len
p = plot(
  layer(x=day_range, y=residuals[day_range], Geom.line, Theme(default_color=color("green"))),
  layer(x=day_range, y=residuals_mat[ar_range, :] * evaluate(ar_coef), Geom.line, Theme(default_color=color("red"))),
  Theme(panel_fill=color("white"))
)
# draw(PNG("ar_fit.png", 16cm, 12cm), p)

total_estimate = evaluate(yearly)
total_estimate[ar_len + 1 : end] += residuals_mat * evaluate(ar_coef)

# plot final fit of data
p = plot(
  layer(x=1:1500, y=total_estimate[1:1500], Geom.line, Theme(default_color=color("red"))),
  layer(x=1:1500, y=temps[1:1500], Geom.line),
  Theme(panel_fill=color("white"))
)
# draw(PNG("total_fit.png", 16cm, 12cm), p)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""The PlanarShape interface:

`distancefromplane(p::PlanarShape,point)`  returns distance of the point from the plane the planar shape lies within
`normal(p::PlanarShape)` returns normal of plane
`interface(p::PlanarShape)` returns optical interface of plane
`vertices(p::PlanarShape)` returns vertices of shape. For Ellipse this is an approximation.

There are default functions for plane,normal,interface,vertices which assume each PlanarShape type has a field of the same name
`
    plane(a::PlanarShape) = a.plane
    normal(a::PlanaShape) = a.plane.normal
`
etc.

If your type doesn't have these fields then you should define a more specialized method to handle this.
"""
abstract type PlanarShape{T} <: Surface{T} end

"""All planar shapes lie on a plane. This function computes the distance from a point to that plane. This is a signed distance. If the point is on the positive side of the plane (the side the normal points toward) the distance will be positive, otherwise negative or 0 if the point lies in the plane."""
distancefromplane(p::PlanarShape, point::SVector{3}) = distancefromplane(p.plane,point)
normal(p::PlanarShape) = normal(p.plane)
interface(p::PlanarShape) = interface(p.plane)
"""The vertices of planar shapes are defined in a plane so they are two dimensional. In the local coordinate frame this is the x,y plane, so the implied z coordinate is 0"""
vertices(p::PlanarShape) = throw(ErrorException("This function should be defined for any concrete type that is a subtype of PlanarShape")) #don't think this function can ever be called because you can't instantiate something of PlanarShape type; it's an abstract type.
export vertices


text: println(return (39 + 3))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
"""
    Cubic() <: Kernel
"""
struct Cubic <: Kernel
    # TODO
end


text:
const ∇! = grad!
const ∇ = grad
const ∫ = integrate


text:
# Try to approximate squared-exponential kernel:
using PyPlot
using Optim

function gaussian_approx5()

ntime = 1000
time = collect(linspace(0,10,ntime))
gaussian=exp(-time.^2/2.)

function compute_chi(x)
model = zeros(ntime)
c1 = 1.0
c2 = 0.0
for k=1:nexp
  model += x[1+(k-1)*2].*exp(-x[2+(k-1)*2].*time)
  c1 -= x[1+(k-1)*2]
  c2 -= x[1+(k-1)*2]*x[2+(k-1)*2]
end
c2 = c2/c1
model += c1.*exp(-c2.*time)
return sum((gaussian-model).^2)
end

nexp = 6
#nexp = 2
#nexp = 4
npar = 2*nexp
abest=zeros(npar)
abest=[2.7518074963750734,1.5312499407621147,-0.3546256365636662,4.753504292338284,-1.3968344706386124,1.7967431129812566,1.0,0.5,0.1,0.1,-0.2,0.5]

atrial = zeros(npar)
for k=1:npar
  atrial[k] = abest[k]
end
model = zeros(ntime)
c1 = 1.0
c2 = 0.0
for k=1:nexp-1
  model += atrial[1+(k-1)*2].*exp(-atrial[2+(k-1)*2].*time)
  c1 -= atrial[1+(k-1)*2]
  c2 -= atrial[1+(k-1)*2]*atrial[2+(k-1)*2]
end
c2 = c2/c1
model += c1.*exp(-c2.*time)
clf()
plot(time,gaussian)
plot(time,model)
#chibest=sum((gaussian-model).^2)
chibest=compute_chi(atrial)

println("Initial chi-square: ",chibest)
read(STDIN,Char)

result = optimize(compute_chi, atrial, BFGS(), OptimizationOptions(autodiff = true))

println(Optim.minimizer(result),Optim.minimum(result))
abest= Optim.minimizer(result)
clf()
plot(time,gaussian)
plot(time,model)
model = zeros(ntime)
c1 = 1.0
c2 = 0.0
for k=1:nexp
  model += abest[1+(k-1)*2].*exp(-abest[2+(k-1)*2].*time)
  c1 -= abest[1+(k-1)*2]
  c2 -= abest[1+(k-1)*2]*abest[2+(k-1)*2]
end
c2 = c2/c1
model += c1.*exp(-c2.*time)
plot(time,(gaussian-model))
chi=sum((gaussian-model).^2)
println("New minimum: ",chi," ",abest," ",c1," ",c2," ",std(gaussian-model)," ",maximum(abs(gaussian-model)))
read(STDIN,Char)

return
end

max_stars_repo_name: dfm/ess

max_issues_repo_name: dfm/ess

max_forks_repo_name: dfm/ess


text:
push!(LOAD_PATH, "/Users/skohn/Programming/projects/Units.jl/src")
using Base.Test
using units

# Tests of BaseUnit

# Factory constructors
@test m == Meter{1}()
@test m^3 == Meter{3}()
@test m^-2.5 == Meter{-2.5}()

# Existence of other units
@test s == Second{1}()
@test kg == Kilogram{1}()
@test rad == Radian{1}()

# In-unit arithmetic
@test m * m == Meter{2}()
@test m / m^3 == Meter{-2}()
@test reciprocal(m) == m^-1

# Tests of Unit

# Constructor
@test Unit(m^5.2) == Unit{Tuple{Meter{5.2}}}()
@test Unit(m^2) == Unit{Tuple{Meter{2}}}()

# In-unit arithmetic
@test m * s == Unit{Tuple{Meter{1}, Second{1}}}()
@test m * s == s * m
@test (m * s) * m == Unit{Tuple{Meter{2}, Second{1}}}()
@test m * (s * m) == Unit{Tuple{Meter{2}, Second{1}}}()
@test (m * s) * m == m * (s * m)
@test m * s * s == Unit{Tuple{Meter{1}, Second{2}}}()
@test reciprocal(m*s) == Unit{ Tuple{Meter{-1}, Second{-1}}}()
@test m / s == Unit{Tuple{Meter{1}, Second{-1}}}()
@test Unit(m) * Unit(s^-1) == m/s
@test m/(kg * s^-1) == Unit{Tuple{Kilogram{-1}, Meter{1}, Second{1}}}()


text:
"""
    module Datasets

Commonly used datasets and utilities for creating data containers.

In the future, contents will be integrated into packages:

- FastAI datasets and data containers will be moved into MLDatasets.jl
- data container transformations will be moved to MLDataPattern.jl

This submodule will then reexport the same definitions.
"""
module Datasets

using ..FastAI
using ..FastAI: typify, TableRow

using DataDeps
using Glob
using FilePathsBase
import DataAugmentation
using FilePathsBase: filename
import FileIO
using IndirectArrays: IndirectArray
using MLDataPattern
using MLDataPattern: splitobs
import LearnBase
using Colors
using FixedPointNumbers
using DataFrames
using Tables
using CSV
using ShowCases

include("fastaidatasets.jl")

function __init__()
    initdatadeps()
end

include("containers.jl")
include("transformations.jl")

include("load.jl")
include("recipes/recipe.jl")
include("recipes/vision.jl")
include("recipes/tabular.jl")
include("registry.jl")
include("fastairegistry.jl")

export
    # reexports from MLDataPattern
    splitobs,

    # container transformations
    mapobs,
    filterobs,
    groupobs,
    joinobs,
    eachobs,

    # primitive containers
    FileDataset,
    TableDataset,

    # utilities
    isimagefile,
    matches,
    loadfile,
    loadmask,
    pathname,
    pathparent,
    parentname,
    grandparentname,

    # datasets
    DATASETS,
    loadfolderdata,
    datasetpath,

    # recipes
    loadrecipe,
    finddatasets,
    listdatasources,
    loaddataset

end  # module


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

name = "FLAC"
version = v"1.3.3"

# Collection of sources required to build FLAC
sources = [
    ArchiveSource("https://downloads.xiph.org/releases/flac/flac-$(version).tar.xz",
                  "213e82bd716c9de6db2f98bcadbc4c24c7e2efe8c75939a1a84e28539c4e1748"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/flac-*/

# Include patch for finding definition of `AT_HWCAP2` within the Linux
# kernel headers, rather than the glibc headers, sicne our glibc is too old
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/flac_linux_headers.patch"

./configure --prefix=$prefix --host=$target  --build=${MACHTYPE}
make -j${nproc}
make install
install_license COPYING.Xiph
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct("libFLAC", :libflac),
    LibraryProduct("libFLAC++", :libflacpp),
    ExecutableProduct("metaflac", :metaflac),
    ExecutableProduct("flac", :flac)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Ogg_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
using GeoData, Test, Dates
using GeoData.LookupArrays, GeoData.Dimensions

data1 = cumsum(cumsum(ones(10, 11); dims=1); dims=2)
data2 = 2cumsum(cumsum(ones(10, 11, 1); dims=1); dims=2)
dims1 = X(10:10:100), Y(-50:10:50)
dims2 = (dims1..., Ti([DateTime(2019)]))
refdimz = ()
mval = -9999.0
meta = NoMetadata()
nme = :test

ga2 = GeoArray(data2, dims2)
ga1 = GeoArray(data1; dims=dims1, refdims=refdimz, name=nme, metadata=meta, missingval=mval)

@test ga1 == data1
@test ga2 == data2

@testset "arary dims have been formatted" begin
    @test index(ga2) == (10.0:10:100, -50.0:10:50.0, [DateTime(2019)])
    @test dims(ga1)[1:2] == dims(ga2)[1:2]
    @test name(ga1) == :test
    @test missingval(ga1) == -9999.0
end

@testset "constructors" begin
    x = X(Projected(1.0:1.0:2.0; crs=ProjString("+proj=")))
    y = Y(Projected(1.0:1.0:2.0; crs=ProjString("+proj=")))

    A1 = rand(Float32, x, y)
    @test A1 isa GeoArray{Float32,2}
    @test size(A1) == (2, 2)

    A2 = zeros(Int, x, y)
    @test A2 isa GeoArray{Int,2}
    @test A2 == [0 0; 0 0]

    A3 = ones(x, y)
    @test A3 isa GeoArray{Float64,2}
    @test A3 == [1.0 1.0; 1.0 1.0]

    A4 = fill(99.0, x, y)
    @test A4 isa GeoArray{Float64,2}
    @test A4 == [99.0 99.0; 99.0 99.0]

    A5 = trues(x, y)
    @test A5 isa GeoArray{Bool,2}
    @test A5 == [true true; true true]

    A6 = falses(x, y)
    @test A6 isa GeoArray{Bool,2}
    @test A6 == [false false; false false]
end

@testset "show" begin
    sh = sprint(show, MIME("text/plain"), ga1)
    # Test but don't lock this down too much
    @test occursin("GeoArray", sh)
    @test occursin("Y", sh)
    @test occursin("X", sh)
end


text:
import Pkg; Pkg.add("Scanf")
using Scanf

### Part 1
N = 50
S = 2 * N + 1
grid = zeros(Bool, S, S, S)

function updateInit(v, xs)
  cut(r) = intersect(r, -N:N) .+ (N + 1)
  xl, xr, yl, yr, zl, zr = xs
  grid[cut(xl:xr), cut(yl:yr), cut(zl:zr)] .= v
end

### Part 2
concatMap(f, xs) = collect(Iterators.flatten(map(f, xs)))

# from p, remove q
function cut1D(p, q)
  q = intersect(p, q)
  if (length(q) == 0) return [p] end
  if (p == q) return [] end
  l, r = extrema(p)
  ll, rr = extrema(q)
  return filter(s -> length(s) != 0, [l:ll-1, rr+1:r])
end

# from c1, cut out c2
function cutOverlap(c1, c2)
  x1, y1, z1, v = c1
  x2, y2, z2, _ = c2
  res = []
  for xx in cut1D(x1, x2)
    push!(res, (xx, y1, z1, v))
  end
  x1 = intersect(x1, x2)
  if (length(x1) == 0) return res end
  for yy in cut1D(y1, y2)
    push!(res, (x1, yy, z1, v))
  end
  y1 = intersect(y1, y2)
  if (length(y1) == 0) return res end
  for zz in cut1D(z1, z2)
    push!(res, (x1, y1, zz, v))
  end
  return res
end

cubes = []
function update(v, xs)
  xl, xr, yl, yr, zl, zr = xs
  cur = (xl:xr, yl:yr, zl:zr, v)
  global cubes = concatMap(c -> cutOverlap(c, cur), cubes)
  push!(cubes, cur)
end

function volume(c)
  x, y, z, v = c
  return if (v) length(x) * length(y) * length(z) else 0 end
end

### Compute
for line in eachline()
  _, state, xs... = @scanf(line, "%s x=%d..%d,y=%d..%d,z=%d..%d",
                           String, Int, Int, Int, Int, Int, Int)
  updateInit(state == "on", xs)
  update(state == "on", xs)
end

println(sum(grid))
println(sum(map(volume, cubes)))


text:
module AFDM
	include("afd_functions_pile.jl")
	export czt, unitDiskGrid, reconstructTimeDomain, afd
end


text:
export continuation_simple

function continuation_simple(state::State; x_perturbation=1e-2, ops=(:zero, :+, :-), dc_tolerance=1e-3, verbose=true, kwargs...)
    # NOTE: Only designed for DiffCorrectAxisymmetric!
    corrector = DiffCorrectAxisymmetric()

    # Try to do a simple continuation from sols
    sols = Dict{Symbol,Vector{Trajectory}}(:zero => [], :+ => [], :- => [])
    orig_state = deepcopy(state)

    # Properties
    props = R3BPSystemProperties(state.model)

    # Limits
    min_x = uconvert(NoUnits, props.R1[1] / props.L) - props.μ
    max_x = (1 - props.μ) - uconvert(NoUnits, props.R2[1] / props.L)

    # Solutions to return
    new_sols::Array{Trajectory} = []

    # Threads.@threads for id in ops
    try
        for id in ops
            op = getfield(Base, id)
            p = ProgressThresh(id == :+ ? max_x : min_x; desc="Continuing orbits along $(op):", enabled=verbose)
            last_state = orig_state
            while true
                # Build a new solution by perturbing the last solution
                new_u0 = convert_to_frame(last_state, corrector.frame).u0
                new_u0[1] += op(x_perturbation)
                new_state = State(state.model, corrector.frame, new_u0, orig_state.tspan)

                # Check conditions
                state_check = convert_to_frame(new_state, SynodicFrame(true))
                if state_check.u0[1] >= max_x || state_check.u0[1] <= min_x
                    # Initial position as not between primary and secondary body...
                    break
                end

                try
                    new_sol = solve(new_state, corrector; dc_tolerance, kwargs...)
                    if new_sol.u[end][1] < (1 - props.μ)
                        # Solution did not cross x-axis on the right side of secondary body...
                        break
                    end

                    # Append the new solution
                    push!(sols[id], new_sol)
                    ProgressMeter.update!(p, state_check.u0[1])
                    last_state = State(new_sol)
                catch ex
                    if isa(ex, AxisymmetricException)
                        break
                    else
                        rethrow(ex)
                    end
                end

                id == :zero && break
            end
            ProgressMeter.finish!(p)
        end
    catch ex
        if isa(ex, InterruptException)
            @warn "Interrupting continuation! Found $(length(sols)) solutions."
        else
            rethrow()
        end
    end

    new_sols = vcat(values(sols)...)  # Concatenate all solutions
    sort!(new_sols, by=x -> x.prob.u0[1])  # Sort solutions by X value

    return new_sols
end


text:
@testset "Utils" begin
    De, Dx = 3, 2
    g = Flux.batch([GNNGraph(erdos_renyi(10, 30),
                             ndata=rand(Dx, 10),
                             edata=rand(De, 30),
                             graph_type=GRAPH_T) for i=1:5])
    x = g.ndata.x
    e = g.edata.e

    @testset "reduce_nodes" begin
        r = reduce_nodes(mean, g, x)
        @test size(r) == (Dx, g.num_graphs)
        @test r[:,2] ≈ mean(getgraph(g, 2).ndata.x, dims=2)
    end

    @testset "reduce_edges" begin
        r = reduce_edges(mean, g, e)
        @test size(r) == (De, g.num_graphs)
        @test r[:,2] ≈ mean(getgraph(g, 2).edata.e, dims=2)
    end

    @testset "softmax_nodes" begin
        r = softmax_nodes(g, x)
        @test size(r) == size(x)
        @test r[:,1:10] ≈ softmax(getgraph(g, 1).ndata.x, dims=2)
    end

    @testset "softmax_edges" begin
        r = softmax_edges(g, e)
        @test size(r) == size(e)
        @test r[:,1:60] ≈ softmax(getgraph(g, 1).edata.e, dims=2)
    end

    @testset "broadcast_nodes" begin
        z = rand(4, g.num_graphs)
        r = broadcast_nodes(g, z)
        @test size(r) == (4, g.num_nodes)
        @test r[:,1] ≈ z[:,1]
        @test r[:,10] ≈ z[:,1]
        @test r[:,11] ≈ z[:,2]
    end

    @testset "broadcast_edges" begin
        z = rand(4, g.num_graphs)
        r = broadcast_edges(g, z)
        @test size(r) == (4, g.num_edges)
        @test r[:,1] ≈ z[:,1]
        @test r[:,60] ≈ z[:,1]
        @test r[:,61] ≈ z[:,2]
    end

    @testset "softmax_edge_neighbors" begin
        s = [1,2,3,4]
        t = [5,5,6,6]
        g2 = GNNGraph(s, t)
        e2 = randn(Float32, 3, g2.num_edges)
        z = softmax_edge_neighbors(g2, e2)
        @test size(z) == size(e2)
        @test z[:,1:2] ≈ NNlib.softmax(e2[:,1:2], dims=2)
        @test z[:,3:4] ≈ NNlib.softmax(e2[:,3:4], dims=2)
    end
end


text:
module IProfile

export @iprofile

PROFILE_LINES = 1
PROFILE_DESCEND = 2
PROFILE_STATE = PROFILE_LINES | PROFILE_DESCEND    # state is a bitfield

# Record of expressions to parse according to the current PROFILE_STATE
PROFILE_EXPR = []

# To avoid the performance penalty of global variables, we design the
# macro to create local variables through a "let" block that shares
# the timing and count data with "reporting" and "clearing" functions.
PROFILE_REPORTS = []  # list of reporting functions
PROFILE_CLEARS = []   # list of clearing functions
PROFILE_TAGS = []     # line #s for all timing variables

# Profile calibration, to compensate for the overhead of calling time()
PROFILE_CALIB = 0
# Do it inside a let block, just like in real profiling, in case of
# extra overhead
let # tlast::UInt64 = 0x0, tnow::UInt64 = 0x0
global profile_calib
function profile_calib(n_iter)
    trec = Array(UInt64, n_iter)
    for i = 1:n_iter
        tlast = time_ns()
        blast = Base.gc_bytes()
        tnow = time_ns()
        bnow = Base.gc_bytes()
        trec[i] = tnow - tlast
    end
    return trec
end
end
PROFILE_CALIB = minimum(profile_calib(100))

# Utilities
# Generic expression type testing
is_expr_head(ex::Expr, s::Symbol) = ex.head == s
is_expr_head(nonex, s::Symbol) = false

# Test whether an expression is a function declaration
function isfuncexpr(ex::LineNumberNode)
    return false
end
function isfuncexpr(ex::Expr)
    return ex.head == :function || (ex.head == :(=) && typeof(ex.args[1]) == Expr && ex.args[1].head == :call)
end

# Get the "full syntax" of the function call
function funcsyntax(ex::Expr)
    return ex.args[1]
end

# Get the symbol associated with the function call
function funcsym(ex::Expr)
    tmp = funcsyntax(ex)
    tmp = tmp.args[1]
    if is_expr_head(tmp, :curly)
        tmp = tmp.args[1]
    end
    return tmp
end

# Test for control-flow statements
#is_cf_expr(ex::Expr) = contains([:for, :while, :if, :try], ex.head)
is_cf_expr(ex::Expr) = in(ex.head, [:for, :while, :if, :try])
is_cf_expr(ex) = false

# General switchyard function
function insert_profile(ex::Expr, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx::Int, retsym, rettest)
    if ex.head == :block
        insert_profile_block(ex, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym, rettest)
    elseif is_cf_expr(ex)
        insert_profile_cf(ex, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym)
    else
        error("Don't know what to do")
    end
end

# A variant for anything but an expression (a no-op).
insert_profile(notex, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx::Int, retsym, rettest) = notex, indx

# Insert profiling statements into a code block
# rettest is a function with the following syntax:
#    rettest(Expr, Int)
# and evaluates to true if the return value of Expr needs to be saved
# before inserting profiling statements.
function insert_profile_block(fblock::Expr, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx::Int, retsym, rettest)
    global PROFILE_STATE, PROFILE_DESCEND
    if fblock.head != :block
        println(fblock)
        error("expression is not a block")
    end
    descend = PROFILE_STATE & PROFILE_DESCEND > 0
    fblocknewargs = []
    local lasttag
    for i = 1:length(fblock.args)
        if isa(fblock.args[i],LineNumberNode) || is_expr_head(fblock.args[i], :line)
            # This is a line expression, so no counters/timers required
            push!(fblocknewargs,fblock.args[i])
             # ...but keep track of the line # for use during reporting
            lasttag = fblock.args[i]
        elseif descend && is_cf_expr(fblock.args[i])
            # This is a control-flow statement, it requires special
            # handling (recursive)
            cfnew, indx = insert_profile_cf(fblock.args[i], tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym)
            push!(fblocknewargs, cfnew)
        else
            # This is an "ordinary" statement
            saveret = rettest(fblock, i)
            push!(tags, Any[lasttag, fblock.args[i]])
            if saveret
                if is_expr_head(fblock.args[i], :return)
                    push!(fblocknewargs, :($retsym = $(fblock.args[i].args[1])))
                else
                    push!(fblocknewargs, :($retsym = $(fblock.args[i])))
                end
            else
                push!(fblocknewargs, fblock.args[i])
            end
            # This next line inserts timing statements between two
            # lines of code, equivalent to:
            #   timehr(PROFILE_USE_CLOCK, tnow)  # end time of prev
            #   timers[indx] += timehr_diff(tlast, tnow)
            #   counters[indx] += 1
            #   timehr(PROFILE_USE_CLOCK, tlast) # start time for next
            append!(fblocknewargs,Any[:($tnow = time_ns()), :($bnow = Base.gc_bytes()), :(($timers)[($indx)] += $tnow - $tlast),  :(($byters)[($indx)] += $bnow - $blast), :(($counters)[($indx)] += 1), :($tlast = time_ns()), :($blast = Base.gc_bytes())])
            indx += 1
            if saveret
                push!(fblocknewargs, :(return $retsym))
            end
        end
    end
    return Expr(:block, fblocknewargs...), indx
end

# Handling control-flow statements
function insert_profile_cf(ex::Expr, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx::Int, retsym)
    rettest = (ex, i) -> is_expr_head(ex, :return)
    if length(ex.args) == 2
        # This is a for, while, or 2-argument if or try block
        block1, indx = insert_profile(ex.args[2], tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym, rettest)
        return Expr(ex.head, Any[ex.args[1], block1]...), indx
    elseif length(ex.args) == 3
        # This is for a 3-argument if or try block
        block1, indx = insert_profile(ex.args[2], tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym, rettest)
        block2, indx = insert_profile(ex.args[3], tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym, rettest)
        return Expr(ex.head, Any[ex.args[1], block1, block2]...), indx
    else
        error("Wrong number of arguments")
    end
end

# Insert timing and counters into function body.
# Function bodies differ from blocks in these respects:
#   - We need to initialize tlast
#   - The final statement of a function should be returned, even if there
#     is no explicit return statement
#   - Functions can be defined in "short-form" (e.g.,
#     "isempty(x) = numel(x)==0"), and the return value for these
#     needs to be managed, too
function insert_profile_function(ex::Expr, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx::Int, retsym)
    fblock = ex.args[2]
    if fblock.head != :block
        error("Can't parse func expression")
    end
    # Prepare the test for whether we need to save the return value of
    # a given line of code.  We may need to store the return value
    # because we need to run timing operations after computing the
    # output.
    # For a function, this will be true in three cases:
    #   - For a "f1(x) = x+1" type of function declaration
    #   - For an explicit return statement
    #   - For the last line of a function that does not have
    #     an explicit return statement in it.
    if ex.head == :(=)
        # This is a short-form function declaration
        savefunc = (ex, i) -> true
    else
        # Long form, check to see if it's a return or the last line
        savefunc = (ex, i) -> i == length(fblock.args) || is_expr_head(ex, :return)
    end
    # Insert the profiling statements in the function
    fblocknewargs, indx = insert_profile_block(fblock, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym, savefunc)
    # Prepend the initialization of tlast
    fblocknewargs = vcat(Any[:($tlast = time_ns())], Any[:($blast = Base.gc_bytes())], fblocknewargs.args)
    return Expr(:function,Any[funcsyntax(ex),Expr(:block,fblocknewargs...)]...), indx
end

function profile_parse(ex::Expr)
    if PROFILE_STATE & PROFILE_LINES > 0
        # Create the "let" variables for timing and counting
        tlast = gensym()
        tnow = gensym()
        timers = gensym()
        blast = gensym()
        bnow = gensym()
        byters = gensym()
        counters = gensym()
        # Keep track of line numbers
        tags = []
        # Preserve return values
        retsym = gensym()
        # Create the symbols used for reporting and clearing the data
        # for this block
        funcreport = gensym()
        funcclear = gensym()
        # Parse the block and insert instructions
        indx = 1
        coreargs = []
        if ex.head == :block
            # This is a block which may contain many function declarations
            for i = 1:length(ex.args)
                if isfuncexpr(ex.args[i])
                    # Insert "global" statement for each function
                    push!(coreargs,Expr(:global,funcsym(ex.args[i])))
                    # Insert function-call counters
                    newfuncexpr, indx = insert_profile_function(ex.args[i], tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym)
                    push!(coreargs, newfuncexpr)
                else
                    push!(coreargs,ex.args[i])
                end
            end
        elseif isfuncexpr(ex)
            # This is a single function declaration
            push!(coreargs,Expr(:global,funcsym(ex)))
            newfuncexpr, indx = insert_profile_function(ex, tlast, tnow, timers, blast, bnow, byters, counters, tags, indx, retsym)
            push!(coreargs, newfuncexpr)
        else
            error("Could not parse expression")
        end
        n_lines = indx-1
        # Insert reporting function
        # Because we're using a gensym for the function name, we can't
        # quote the whole thing
        push!(coreargs, Expr(:global, funcreport))
        push!(coreargs, Expr(:function, Any[Expr(:call, Any[funcreport]...), Expr(:block,Any[:(return $timers, $byters, $counters)]...)]...))
        # Insert clearing function
        push!(coreargs, Expr(:global, funcclear))
        push!(coreargs, Expr(:function, Any[Expr(:call, Any[funcclear]...), Expr(:block,Any[:(fill!($timers,0)), :(fill!($byters,0)), :(fill!($counters,0))]...)]...))
        # Put all this inside a let block
        excore = Expr(:block,coreargs...)
        exlet = Expr(:let,Any[Expr(:block,excore), :($timers = zeros(UInt64, $n_lines)), :($byters = zeros(UInt64, $n_lines)), :($counters = zeros(UInt64, $n_lines))]...)
        # Export the reporting and clearing functions, in case we're inside a module
        exret = Expr(:toplevel, Any[esc(exlet), Expr(:export, Any[esc(funcclear), esc(funcreport)]...)]...)
        return exret, tags, funcreport, funcclear
    else
        return ex ,[], :funcnoop, :funcnoop
    end
end

function funcnoop()
end

function profile_parse_all()
    empty!(PROFILE_REPORTS)
    empty!(PROFILE_CLEARS)
    empty!(PROFILE_TAGS)
    retargs = []
    for i = 1:length(PROFILE_EXPR)
        newblock, tags, funcreport, funcclear = profile_parse(PROFILE_EXPR[i])
        retargs = vcat(retargs, newblock.args)
        if !isempty(tags)
            push!(PROFILE_TAGS, tags)
            push!(PROFILE_REPORTS, esc(funcreport))
            push!(PROFILE_CLEARS, esc(funcclear))
        end
    end
    push!(retargs,:(return nothing))
    return esc(Expr(:block,retargs...))
end

function profile_report()
    exret = cell(length(PROFILE_REPORTS)+2)
    ret = gensym()
    exret[1] = :($ret = [])
    for i = 1:length(PROFILE_REPORTS)
        exret[i+1] = :(push!($ret,$(Expr(:call,Any[PROFILE_REPORTS[i]]...))))
    end
    exret[end] = :(profile_print($ret))
    return Expr(:block,exret...)
end

compensated_time(t, c) = t >= c*PROFILE_CALIB ? t-c*PROFILE_CALIB : 0
function show_unquoted(arg::Array{Any,1})
  if length(arg[1].args)==1
    ret_val = string("\t# line ", arg[1].args[1])
  else
    ret_val = string("\t# ", arg[1].args[2], ", line ", arg[1].args[1])
  end
  ret_val = string(ret_val, ", ", string(arg[2]))
end
#show_unquoted(linex::Expr) = show_linenumber(linex.args...)
#show_unquoted(lnn::LineNumberNode) = show_linenumber(lnn.line)
#show_linenumber(line)       = string("\t#  line ",line)
#show_linenumber(line, file) = string("\t#  ",file,", line ",line)

function profile_print(tc)
    # Compute total elapsed time
    ttotal = 0.0
    btotal = 0.0
    for i = 1:length(tc)
        timers = tc[i][1]
        byters = tc[i][2]
        counters = tc[i][3]
        for j = 1:length(counters)
            comp_time = compensated_time(timers[j], counters[j])
            ttotal += comp_time
            btotal += byters[j]
        end
    end
    # Display output
    for i = 1:length(tc)
        timers = tc[i][1]
        byters = tc[i][2]
        counters = tc[i][3]
        println("  count  time(%)   time(s) bytes(%) bytes(k)")
        for j = 1:length(counters)
            if counters[j] != 0
                comp_time = compensated_time(timers[j], counters[j])
                @printf("%8d  %5.2f  %9.6f  %5.2f  %8d  %s\n", counters[j],
                        100*(comp_time/ttotal),
                        comp_time*1e-9,
                        btotal > 0.0? 100*(byters[j]/btotal) : 0.0,
                        byters[j] / 1e3,
                        show_unquoted(PROFILE_TAGS[i][j]))
            end
        end
    end
end

function profile_clear()
    exret = cell(length(PROFILE_CLEARS)+1)
    for i = 1:length(PROFILE_CLEARS)
        exret[i] = Expr(:call,Any[PROFILE_CLEARS[i]]...)
    end
    exret[end] = :(return nothing)
    return Expr(:block,exret...)
end

macro iprofile(ex)
    global PROFILE_STATE
    if isa(ex,Symbol)
        # State changes
        if ex == :off
            PROFILE_STATE = 0
        elseif ex == :on
            PROFILE_STATE = PROFILE_LINES
        elseif ex == :reset
        elseif ex == :report
            return profile_report()
        elseif ex == :clear
            return profile_clear()
        else
            error("Profile mode not recognized")
        end
        return profile_parse_all()
    elseif isa(ex,Expr)
        push!(PROFILE_EXPR,ex)
        exret, tags, funcreport, funcclear = profile_parse(ex)
        if !isempty(tags)
            push!(PROFILE_TAGS, tags)
            push!(PROFILE_REPORTS, esc(funcreport))
            push!(PROFILE_CLEARS, esc(funcclear))
        end
        return exret
    end
end

export @iprofile

end # module Profile


text:
# Decision stumps
module DecisionStump

	using StaticArrays, MLDataPattern
	using StatsBase: entropy, trim, sample
	using Distances: sqeuclidean
	using j4pr: countappw, countappw!, gini, misclassification, linearsplit, densitysplit
	export AbstractStump, StumpClassifier, StumpRegressor, stump_train, stump_exec

	# Types
	abstract type AbstractStump end

	struct ClassificationStump{N} <: AbstractStump
		idx::Int			# Index of the threshold variable
		th::SVector{N}			# Thresholds (N=1 for ordinal/real split, N>1 for nominal split)
		prob::Matrix{Float64}		# Class probabilities corresponding to each class (lines) and threshold (columns)

		# Define consistency checks in inner constructor
		function ClassificationStump(idx::Int, th::SVector{N}, prob::Matrix{Float64}) where N
			@assert N == size(prob,2)-1 "The number of thresholds + 1 has to match the number of columns in the probability matrix"
			return new{N}(idx, th, prob)
		end
	end

	struct RegressionStump{N} <: AbstractStump
		idx::Int			# Index of the threshold variable
		th::SVector{N}			# Thresholds (N=1 for ordinal/real split, N>1 for nominal split)
		coeff::Matrix{Float64}		# Regression coefficient matrix (the number of rows is equivalent to the
						# order of the polynomial fitted i.e. 1 row, a constant, 2 rows a line etc.)

		# Define consistency checks in inner constructor
		function RegressionStump(idx::Int, th::SVector{N}, coeff::Matrix{Float64}) where N
			@assert N == size(coeff,2)-1 "The number of thresholds has to match the number of columns in the regression coefficient matrix"
			return new{N}(idx, th, coeff)
		end
	end

	# Aliases
	const ClassificationStumpOR = ClassificationStump{1}					# stump classifier - ordinal or real variable
	const ClassificationStumpN{N} = ClassificationStump{N} where N 				# stump classifier - nominal variable
	const RegressionStumpOR = RegressionStump{1}						# stump regressor - ordinal or real variable
	const RegressionStumpN = RegressionStump{N} where N					# stump regressor - nominal variable

	# Printers
	Base.show(io::IO, m::ClassificationStumpOR) = print(io, "Decision stump, classification, var=$(m.idx), th=$(m.th[1])")
	Base.show(io::IO, m::ClassificationStumpN) = print(io, "Decision stump, classification, var=$(m.idx), $(length(m.th)) thresholds")
	Base.show(io::IO, m::RegressionStumpOR) = print(io, "Decision stump, regression, var=$(m.idx), th=$(m.th[1]), order $(size(m.coeff,2)-1) fit")
	Base.show(io::IO, m::RegressionStumpN) = print(io, "Decision stump, regression, var=$(m.idx), $(length(m.th)) thresholds, order 0 fit")

	# Calculate the information gain for ordinal or real variables, for all thresholds provided
	_gain_real_(x::AbstractVector{Float64}, y::AbstractVector{Int}, yu::AbstractVector{Int}, th::Vector{Float64}, fp::Function) = begin
		# Initialize
		nt = length(th)									# number of thresholds
		n = length(y)									# number of samples
		Imax::Float64 = -Inf;
		I::Float64=0.0
		pos::Int = -1;
		nleft::Float64 = 0.0
		@inbounds for i = 1:nt
			nleft = sum(x.<=th[i])/n
			I = fp(countappw(y,yu)) - nleft*fp(countappw(view(y,x.<=th[i]),yu)) - (1-nleft)*fp(countappw(view(y,x.>th[i]),yu))
			# Check if the information gain was improved and update
			if I >= Imax
				Imax = I							# update maximum gain
				pos = i 							# update threshold index
			end
		end
		if pos == -1
			# No split is better
			return Imax, [th[div(nt/2)]], [countappw(y,yu) countappw(y,yu)]
		else
			# Return split found
			return Imax, th[pos:pos], [countappw(view(y,x.<=th[pos]),yu) countappw(view(y,x.>th[pos]),yu)]
		end
	end

	# Calculate the information gain for nominal variables, for all thresholds provided
	_gain_nominal_(x::AbstractVector{Float64}, y::AbstractVector{Int}, yu::AbstractVector{Int}, th::Vector{Float64}, fp::Function) = begin
		# For nominal splits, N+1 thresholds, if N is the number of distinct values of the variable, the additional one for 'unknown'
		# 'th' is ignored
		ux = unique(x) 									# unique values of the variable (i.e. the thresholds)
		nt = length(ux)									# number of thresholds
		n = length(y)  									# number of samples
		C = length(yu)									# number of classes
		priors = [sum(yi.==y)/(n+C) for yi in yu]					# priors
		probs = zeros(Float64, C, nt+1)							# preallocate probability matrix (last row for unseen values)
		probs[:,end] = priors

		@inbounds for i = 1:nt
			ym = view(y,x.==ux[i])							# the labels for a the leaf ux[i]
			#probs[:,i] = length(ym)/(n+C)*countappw(ym,yu) # with re-normalization to the total number of samples
			probs[:,i] = countappw(ym,yu) 						# count
		end

		# Calculate gain
		Isplit = 0.0
		@inbounds @simd for i in 1:nt
			Isplit += fp(view(probs,:,i))
		end
		I = fp(priors) - Isplit # no actual threshold selection is made
		return I, ux, probs
	end

	# Calculate the fit error for ordinal or real variables, for all thresholds provided
	_regressor_real_(x::AbstractVector{Float64}, y::AbstractVector{Float64}, th::Vector{Float64}, fm::Function, fe::Function) = begin
		# Initialize
		nt = length(th)									# number of thresholds
		n = length(y)									# number of samples
		Emin::Float64 = Inf;
		E::Float64=0.0
		pos::Int = -1;
		nleft::Float64 = 0.0
		@inbounds for i = 1:nt
			nleft = sum(x.<=th[i])
			(nleft == 0 || nleft == n) && continue 					# skip theshold if extreme
			E = -fe(y,fm(x,y)[2]) + fe(view(y,x.<=th[i]), fm(view(x,x.<=th[i]),view(y,x.<=th[i]))[2]) +
			                        fe(view(y,x.>th[i]), fm(view(x,x.>th[i]),view(y,x.>th[i]))[2])
			# Check if the error was improved and update
			if E <= Emin
				Emin = E							# update error difference
				pos = i 							# update threshold index
			end
		end

		if pos == -1
			return Emin, [th[div(nt,2)]], [fm(x,y)[1] fm(x,y)[1]]
		else
			return Emin, th[pos:pos], [fm(view(x,x.<=th[pos]),view(y,x.<=th[pos]))[1] fm(view(x,x.>th[pos]),view(y,x.>th[pos]))[1]]
		end
	end

	# Calculate the fit error for nominal variables, for all thresholds provided
	_regressor_nominal_(x::AbstractVector{Float64}, y::AbstractVector{Float64}, th::Vector{Float64}, fm::Function, fe::Function) = begin
		# For nominal splits, N+1 thresholds, if N is the number of distinct values of the variable, the additional one for 'unknown'
		# 'th' is ignored
		ux = sort(unique(x))								# unique values of the variable (i.e. the thresholds)
		nt = length(ux)									# number of thresholds
		n = length(y)  									# number of samples
		deflt = fm(x,y)[1]								# default coefficients if value is not found
		coeff = zeros(Float64, 2, nt+1)							# preallocate coefficient matrix (last row for unseen values)
		coeff[:,end] = deflt

		Esplit = 0.0
		@inbounds @simd for i = 1:nt
			ym = view(y,x.==ux[i])							# select values corresponding to current unique data value
			c, ŷ = fm(fill(ux[i], length(ym)), ym)					# fit model
			coeff[:,i] = c								# get coefficients
			Esplit += fe(ym,ŷ)							# update split error
		end

		# Calculate the error difference
		E = Esplit - fe(y,fm(x,y)[2]) # no actual threshold selection is made
		return E, ux, coeff
	end

	"""
		stump_train(x,y [;kwargs])

	Trains a decision stump classifier using the data `x` and labels `y`.

	# Arguments
	  * `x::AbstractMatrix{Float64}` the input data; it is assumed that each column of the input matrix represents an observation
	  * `y::AbstractVector{Int}` the labels vector

	# Keyword arguments
	  * `crit::Symbol` is the purity criterion used to assess the split. Supported: `:gini`,`:entropy` and `:misclassification` (default `:gini`)
	  * `nthresh::Int` is the number of thresholds for a ordinal or real variable that will be investigated (default `10`)
	  * `split::Symbol` describes the variable split used in generating the thresholds. Supported: `:linear`, `:density` (default `:linear`)
	  * `count::Int` is the number of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0`)
	  * `prop::Float64` is the proportion of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0.0`)
	  * `vartypes::Union{Dict{Int=>Symbol}, Symbol}` describes the type of variables from the dataset; if it is a `Symbol`, all variables
	are assumed to be of the same type. If it is a `Dict{Int=>:Symbol}`, the variable indicated by the key is assumed to be of the type of the value;
	two values are supported: `:nominal` and any other symbol; if a variable is specified as nominal, the classifier will have `N+1` leaves, where `N` is
	the number of distinct values seen in training, the last leaf corresponding to unseen values. If the variable is real, only two leaves are present.
	"""
	# Training method (classification)
	function stump_train(x::AbstractMatrix{Float64}, y::AbstractVector{Int};
		      		crit::Symbol=:gini,
		      		vartypes::Union{Dict{Int,Symbol},Symbol}=:real,
		      		nthresh::Int=10, split::Symbol=:linear,
				count::Int=0, prop::Float64=0.0)
		m = size(x,1)					# number of variables
		n = size(x,2)   				# number of observations
		yu = sort(unique(y)) 				# unique classes
		C = length(yu)
		@assert C > 1 "[stump_train] Number of classes has to be larger than 1."

		# Parse purity criterion
		if crit == :gini
			purity_function = gini
		elseif crit == :entropy
			purity_function = entropy
		elseif crit == :misclassification
			purity_function = misclassification
		else
			warn("[stump_train] Unrecognized purity criterion, defaulting to ':gini'.")
			purity_function = gini
		end

		# Parse number of thresholds
		nthresh = max(nthresh, 1) # at least one

		# Parse variable split
		if split == :linear
			split_function = linearsplit
		elseif split == :density
			split_function = densitysplit
		else
			warn("[stump_train] Unrecognized variable split, defaulting to ':linear' split.")
			split_function = linearsplit
		end

		# Parse count and prop
		count = clamp(count, 0, div(n,4)) 	# do no remove more than half the values of the variable
		prop = clamp(prop, 0.0, 0.25) 		# do not remove more than half of the values of the variable

		# Parse variable types
		_parse_vartypes_!(gain, vartypes::Symbol) = begin
			for i in 1:length(gain)
				if vartypes == :nominal
					gain[i] = _gain_nominal_
				else
					gain[i] = _gain_real_
				end
			end
			return gain
		end
		_parse_vartypes_!(gain, vartypes::Dict{Int,Symbol}) = begin
			for i in 1:length(gain)
				if i in keys(vartypes)
					if vartypes[i] == :nominal
						gain[i] = _gain_nominal_
					else
						gain[i] = _gain_real_
					end
				else
					gain[i] = _gain_real_
				end
			end
			return gain
		end

		# Parse the variable specifications to obtain the gain functions corresponding to each
		gain_functions = Vector{Function}(m)
		_parse_vartypes_!(gain_functions, vartypes)

		# Loop through variables, looking at the specified variable type
		Imax = -Inf;
		vidx = 0;									# the index of the best variable
		th = -Inf
		prob = Matrix(Float64[])
		@inbounds for i in 1:m
			variable = view(x,i,:)							# variable data
			tv = split_function(variable, nthresh, prop=prop, count=count)		# thresholds to be investivated
			I, t, v = gain_functions[i](variable, y, yu, tv, purity_function) 	# get the information gain, threshold and split for

			if I > Imax
				vidx = i; th = t; prob = v
			end

		end
		return ClassificationStump(vidx, SVector(th...), prob)
	end

	"""
		stump_train(x,y [;kwargs])

	Trains a decision stump regressor using the data `x` and regression targets `y`.

	# Arguments
	  * `x::AbstractMatrix{Float64}` the input data; it is assumed that each column of the input matrix represents an observation
	  * `y::AbstractVector{Float64}` the regression targets vector

	# Keyword arguments
	  * `model::Symbol` is the type of model assumed for the data. Supported: `:mean`,`:median` and `:linear` (default `:mean`); regardless of
	the model, there are always two coefficients / leaf of the stump, corresponding to the linear and constant terms of a linear model.
	  * `errcrit::Function` A function that returns a measure of the difference between two vectors or scalars (default `Distances.sqeuclidean`)
	  * `nthresh::Int` is the number of thresholds for a ordinal or real variable that will be investigated (default `10`)
	  * `split::Symbol` describes the variable split used in generating the thresholds. Supported: `:linear`, `:density` (default `:linear`)
	  * `count::Int` is the number of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0`)
	  * `prop::Int` is the proportion of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0.0`)
	  * `vartypes::Union{Dict{Int=>Symbol}, Symbol}` describes the type of variables from the dataset; if it is a `Symbol`, all variables
	are assumed to be of the same type. If it is a `Dict{Int=>:Symbol}`, the variable indicated by the key is assumed to be of the type of the value;
	two values are supported: `:nominal` and any other symbol; if a variable is specified as nominal, the regressor will have `N+1` sets of coefficients,
	where `N` is the number of distinct variable values seen in training, the last coefficients corresponding to unseen values.
	If the variable is real, only two sets of coefficients are present.

	!!! Note
	  `model=:linear` is still experimental and it might produce bad results for small values or singular data matrices.
	"""
	# Training method (regression)
	function stump_train(x::AbstractMatrix{Float64}, y::AbstractVector{Float64};
				model::Symbol=:mean,
		      		errcrit::Function=sqeuclidean,
		      		vartypes::Union{Dict{Int,Symbol},Symbol}=:real,
		      		nthresh::Int=10, split::Symbol=:linear,
				count::Int=0, prop::Float64=0.0)
		m = size(x,1)					# number of variables
		n = size(x,2)   				# number of observations

		# Define regression functions
		# These have the same signature and same return arguments:
		#	- a 2 coefficient regression vector
		#       - the fit of that model considering the inut data
		_mean_ = (x,y::AbstractVector{Float64})->begin
			m = mean(y)
			return [0.0,m], fill(m,length(y))
		end

		_median_ = (x,y::AbstractVector{Float64})->begin
			m = median(y)
			return [0.0,m], fill(m,length(y))
		end

		_linear_ = (x::AbstractVector{Float64},y::AbstractVector{Float64})->begin
			if length(x)==1
				m=[0.0,y[1]]
			else
				m=[x+1e-12*rand() ones(x)]\y+1e-12*rand() #add just a bit of noise
			end
			return m, m[1].*x.+m[2] # ŷ
		end

		# Parse model type
		if  model == :mean
			model_function = _mean_
		elseif model == :median
			model_function = _median_
		elseif model == :linear
			model_function = _linear_
		else
			warn("[stump_train] Unrecognized model, defaulting to ':mean'.")
			model_function = _mean_
		end

		# Parse number of thresholds
		nthresh = max(nthresh, 1) # at least one

		# Parse variable split
		if split == :linear
			split_function = linearsplit
		elseif split == :density
			split_function = densitysplit
		else
			warn("[stump_train] Unrecognized variable split, defaulting to ':linear' split.")
			split_function = linearsplit
		end

		# Parse count and prop
		count = clamp(count, 0, div(n,4)) 	# do no remove more than half the values of the variable
		prop = clamp(prop, 0.0, 0.25) 		# do not remove more than half of the values of the variable

		# Parse variable types
		_parse_vartypes_!(regressor, vartypes::Symbol) = begin
			for i in 1:length(regressor)
				if vartypes == :nominal
					regressor[i] = _regressor_nominal_
					if model == :linear
						warn("[stump_train] A ':linear' regression model is not supported with ':nominal' variables. Changing model to ':mean'.")
						model_function = _mean_
					end
				else
					regressor[i] = _regressor_real_
				end
			end
			return regressor
		end
		_parse_vartypes_!(regressor, vartypes::Dict{Int,Symbol}) = begin
			for i in 1:length(regressor)
				if i in keys(vartypes)
					if vartypes[i] == :nominal
						regressor[i] = _regressor_nominal_
						if model == :linear
							warn("[stump_train] A ':linear' regression model is not supported with ':nominal' variables. Changing model to ':mean'.")
							model_function = _mean_
						end
					else
						regressor[i] = _regressor_real_
					end
				else
					regressor[i] = _regressor_real_
				end
			end
			return regressor
		end

		# Parse the variable specifications to obtain the corresponding regressor functions
		regressor_functions = Vector{Function}(m)
		_parse_vartypes_!(regressor_functions, vartypes)

		# Loop through variables, looking at the specified variable type
		Emin = Inf;
		vidx = 0;									# the index of the best variable
		th = Inf
		coeff = Matrix(Float64[])
		@inbounds for i in 1:m
			variable = view(x,i,:)							# variable data
			tv = split_function(variable, nthresh, prop=prop, count=count)		# thresholds to be investivated
			E, t, c = regressor_functions[i](variable, y, tv, model_function, errcrit) # get the error, threshold and split for

			if E < Emin
				vidx = i; th = t; coeff = c
			end

		end
		return RegressionStump(vidx, SVector(th...), coeff)
	end

	# Execution methods (classification)
	stump_exec(m::ClassificationStumpOR, x::Matrix{Float64})=begin
		out = zeros(Float64, size(m.prob,1), size(x,2))
		@inbounds for j = 1:size(x,2)
			# Check input value for the variable: take left side of prob if smaller, right if larger
			out[:,j] = ifelse(x[m.idx,j] <= m.th[1], m.prob[:,1], m.prob[:,2])
		end
		return out
	end

	stump_exec(m::ClassificationStumpN, x::Matrix{Float64})=begin
		out = zeros(Float64, size(m.prob,1), size(x,2))
		@inbounds for j = 1:size(x,2)
			# Search for input value of the variable in thresholds:
			# take corresponding prob if found, default if not
			#out[:,j] = ifelse(x[m.idx,j] in m.th, m.prob[:,indexin(x[m.idx,j:j],m.th)], m.prob[:,end])
			if x[m.idx,j] in m.th
				out[:,j] = m.prob[:,indexin(x[m.idx,j:j],m.th)]
			else
				out[:,j] = m.prob[:,end]
			end
		end
		return out
	end

	# Execution methods (regression)
	stump_exec(m::RegressionStumpOR, x::Matrix{Float64})=begin
		out = zeros(Float64, 1, size(x,2))
		@inbounds @simd for j = 1:size(x,2)
			# Check input value for the variable: take left side coefficients if smaller, right if larger
			out[:,j] = ifelse(x[m.idx,j] <= m.th[1], x[m.idx,j]*m.coeff[1,1]+m.coeff[2,1], # left side
		     						 x[m.idx,j]*m.coeff[1,2]+m.coeff[2,2]  # right side
		     )
		end
		return out
	end

	stump_exec(m::RegressionStumpN, x::Matrix{Float64})=begin
		out = zeros(Float64, 1, size(x,2))
		@inbounds @simd  for j = 1:size(x,2)
			# Search for closest position in thresholds, not an exact match
			pos = findlast(x[m.idx,j].>m.th)
			pos = ifelse(pos==0,1,pos)
			out[:,j] = x[m.idx,j]*m.coeff[1,pos] + m.coeff[2,pos]
		end
		return out
	end

end

##########################
# FunctionCell Interface #
##########################
"""
	stump([;kwargs])

Constructs an untrained cell that when piped data inside, returns a decision stump
classifier trained function cell based on the input data and labels.

# Keyword arguments
  * `crit::Symbol` is the purity criterion used to assess the split. Supported: `:gini`,`:entropy` and `:misclassification` (default `:gini`)
  * `nthresh::Int` is the number of thresholds for a ordinal or real variable that will be investigated (default `10`)
  * `split::Symbol` describes the variable split used in generating the thresholds. Supported: `:linear`, `:density` (default `:linear`)
  * `count::Int` is the number of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0`)
  * `prop::Int` is the proportion of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0.0`)
  * `vartypes::Union{Dict{Int=>Symbol}, Symbol}` describes the type of variables from the dataset; if it is a `Symbol`, all variables
are assumed to be of the same type. If it is a `Dict{Int=>:Symbol}`, the variable indicated by the key is assumed to be of the type of the value;
two values are supported: `:nominal` and any other symbol; if a variable is specified as nominal, the classifier will have `N+1` leaves, where `N` is
the number of distinct values seen in training, the last leaf corresponding to unseen values. If the variable is real, only two leaves are present.
"""
stump(;kwargs...) = FunctionCell(stump, (), ModelProperties(), kwtitle("Decision stump classifier", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	stump(x [;kwargs])

Trains a decision stump classifcation model that using the data `x`.
"""
# Training
stump(x::T where T<:CellDataL; kwargs...) = stump((getx!(x), gety(x)); kwargs...)
stump(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector; kwargs...) =
	stump((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]); kwargs...)
stump(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[stump] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	stumpdata = DecisionStump.stump_train(getobs(x[1]), yenc; kwargs...)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(stump, Model(stumpdata, modelprops), kwtitle("Decision stump classifier, $(length(stumpdata.th)) threshold(s)", kwargs))
end

# Execution
stump(x::T where T<:CellData, model::Model{<:DecisionStump.ClassificationStump}) =
	datacell(stump(getx!(x), model), gety(x))

stump(x::T where T<:AbstractVector, model::Model{<:DecisionStump.ClassificationStump}) =
	stump(mat(x, LearnBase.ObsDim.Constant{2}()), model)

stump(x::T where T<:AbstractMatrix, model::Model{<:DecisionStump.ClassificationStump}) =
	DecisionStump.stump_exec(model.data, getobs(x))

##########################
# FunctionCell Interface #
##########################
"""
	stumpr([;kwargs])

Constructs an untrained cell that when piped data inside, returns a decision stump
trained regressor function cell.

# Keyword arguments
  * `model::Symbol` is the type of model assumed for the data. Supported: `:mean`,`:median` and `:linear` (default `:mean`); regardless of
the model, there are always two coefficients / leaf of the stump, corresponding to the linear and constant terms of a linear model.
  * `errcrit::Function` A function that returns a measure of the difference between two vectors or scalars (default `Distances.sqeuclidean`)
  * `nthresh::Int` is the number of thresholds for a ordinal or real variable that will be investigated (default `10`)
  * `split::Symbol` describes the variable split used in generating the thresholds. Supported: `:linear`, `:density` (default `:linear`)
  * `count::Int` is the number of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0`)
  * `prop::Int` is the proportion of samples in a variable that are discarded at both ends of its value interval in training; see `StatsBase.trim` (default `0.0`)
  * `vartypes::Union{Dict{Int=>Symbol}, Symbol}` describes the type of variables from the dataset; if it is a `Symbol`, all variables
are assumed to be of the same type. If it is a `Dict{Int=>:Symbol}`, the variable indicated by the key is assumed to be of the type of the value;
two values are supported: `:nominal` and any other symbol; if a variable is specified as nominal, the regressor will have `N+1` sets of coefficients,
where `N` is the number of distinct variable values seen in training, the last coefficients corresponding to unseen values.
If the variable is real, only two sets of coefficients are present.

!!! Note
  `model=:linear` is still experimental and it might produce bad results for small values or singular data matrices.
"""
stumpr(;kwargs...) = FunctionCell(stumpr, (), ModelProperties(), kwtitle("Decision stump regressor", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	stumpr(x, [;kwargs])

Trains a decision stump regression model using the data `x`.
"""
# Training
stumpr(x::T where T<:CellDataL; kwargs...) = stumpr((getx!(x), gety(x)); kwargs...)
stumpr(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector; kwargs...) =
	stumpr((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]); kwargs...)
stumpr(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[stumpr] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Train model
	stumpdata = DecisionStump.stump_train(getobs(x[1]), getobs(x[2]); kwargs...)

	FunctionCell(stumpr, Model(stumpdata, ModelProperties(nvars(x[1]),1)),
	      kwtitle("Decision stump regressor", kwargs))

end

# Execution
stumpr(x::T where T<:CellData, model::Model{<:DecisionStump.RegressionStump}) =
	datacell(stumpr(getx!(x), model), gety(x))
stumpr(x::T where T<:AbstractVector, model::Model{<:DecisionStump.RegressionStump}) =
	stumpr(mat(x, LearnBase.ObsDim.Constant{2}()), model)
stumpr(x::T where T<:AbstractMatrix, model::Model{<:DecisionStump.RegressionStump}) =
	DecisionStump.stump_exec(model.data, getobs(x))


text:
# module for drawing moving robots
# one of them move normally and the others sometimes are kidnapped

module AnimeMoveKidnap
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function animate_per_time(time_interval, world, kidnapped_robot,
                            no_kidnapped_robot)
    draw(world)
    annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
    draw!(kidnapped_robot)
    draw!(no_kidnapped_robot)
  end

  function main(time_interval=0.1; is_test=false)
    w = World(-5.0, 5.0, -5.0, 5.0)

    init_pose = [0.0, 0.0, 0.0]
    circling = Agent(0.2, 10.0/180*pi)

    kr = DifferentialWheeledRobot(init_pose, 0.2, "gray",
                                  circling, time_interval,
                                  exp_kidnap_time=3.0,
                                  kidnap_rx=[-5.0, 5.0],
                                  kidnap_ry=[-5.0, 5.0])

    nkr = DifferentialWheeledRobot(init_pose, 0.2, "red",
                                   circling, time_interval)

    if is_test == false
      anime = @animate for t in 0:time_interval:30
        animate_per_time(t, w, kr, nkr)
      end
      path = "src/model/uncertainty/movement/kidnap/anime_move_kidnap.gif"
      gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], path))
    else
      for t in 0:time_interval:10
        animate_per_time(t, w, kr, nkr)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
    abstract type NodeOp{T} end

Represent a time-series operation whose output will be a time-series with value type `T`.
"""
abstract type NodeOp{T} end

"""
    Node(parents, op)

A node in the computational graph that combines zero or more `parents` with `op` to produce
a timeseries.

!!! warning
    Note that a `Node` is only declared mutable so that we can attach finalizers to
    instances. This is required for the [`WeakIdentityMap`](@ref) to work.
    Nodes should NEVER actually be mutated!

Due to subgraph elimination, nodes that are equivalent should always be identical objects.
We therefore leave `hash` & `==` defined in terms of the `objectid`.
"""
mutable struct Node{T}
    parents::NTuple{N,Node} where {N}
    op::NodeOp{T}
end

# Node & NodeOps are immutable for duplication purposes.
duplicate_internal(x::Node, ::IdDict) = x
duplicate_internal(x::NodeOp, ::IdDict) = x

Base.show(io::IO, node::Node) = show(io, node.op)

# Enable AbstractTrees to understand the graph.
# TODO It might be nice to elide repeated subtrees. This would require modifying the
#   iteration procedure within AbstractTrees, so ostriching for now.
AbstractTrees.children(node::Node) = parents(node)
AbstractTrees.nodetype(::Node) = Node

"""
    value_type(node::Node{T}) -> T

The type of each value emitted for this node.
"""
value_type(::Node{T}) where {T} = T
value_type(::NodeOp{T}) where {T} = T

"""
    abstract type NodeEvaluationState end

Represents any and all state that a node must retain between evaluating batches.

Instances of subtypes of `NodeEvaluationState` are given to [`run_node!`](@ref).
"""
abstract type NodeEvaluationState end

"""
    EmptyNodeEvaluationState()

A [`NodeEvaluationState`](@ref) which has no content, to be used by non-stateful nodes.

In practice, the common instance [`TimeDag.EMPTY_NODE_STATE`](@ref) can be used.
"""
struct EmptyNodeEvaluationState <: NodeEvaluationState end

# Can have a singleton instance, since it is just a placeholder.
"""
    const EMPTY_NODE_STATE

A singleton instance of [`TimeDag.EmptyNodeEvaluationState`](@ref).
"""
const EMPTY_NODE_STATE = EmptyNodeEvaluationState()

"""
    parents(node::Node) -> NTuple{N, Node} where {N}

Get immediate parents of the given node.
"""
parents(node::Node) = node.parents

"""
    create_evaluation_state(node::Node) -> NodeEvaluationState
    create_evaluation_state(parents, node::NodeOp) -> NodeEvaluationState

Create an empty evaluation state for the given node, when starting evaluation at the
specified time.
"""
create_evaluation_state(node::Node) = create_evaluation_state(node.parents, node.op)

"""
    run_node!(
        op::NodeOp{T},
        state::NodeEvaluationState,
        time_start::DateTime,
        time_end::DateTime,
        input_blocks::Block...
    ) -> Block{T}

Evaluate the given node from `time_start` until `time_end`, with the initial `state`.
Zero or more blocks will be passed as an input; these correspond to the parents of a node,
and are passed in the same order as that returned by `parents(node)`.

We return a new [`Block`](@ref) of output knots from this node.

!!! warning
    The implementer of `run_node!` must ensure:
    * **No future peeking** occurs: i.e. that no output knot is dependent on input knots
        that occur subsequently.
    * **Correct time range**: all output timestamps must be in the interval
        `[time_start, time_end)`.
    * **Consistency**: Calling `run_node!` over a single interval should give the same
        result as calling it multiple times over a decomposition of that same interval.
    * **Determinism**: `run_node!` should *always* be fully deterministic. If a
        pseudo-random number generator is required, it should be held on the evaluation
        state.
"""
function run_node! end

"""
    stateless(node) -> Bool
    stateless(op) -> Bool

Returns true iff `op` can be assumed to be stateless; that is, if the node evaluation state
is [`TimeDag.EMPTY_NODE_STATE`](@ref).
"""
stateless(node::Node) = stateless(node.op)
stateless(::NodeOp) = false

"""
    duplicate(x)

Return an object that is equal to `x`, but fully independent of it.

Note that for any parts of `x` that `TimeDag` considers to be immutable (e.g. `Block`s),
this can return the identical object.

Conceptually this is otherwise very similar to `deepcopy(x)`.
"""
function duplicate(x)
    isbitstype(typeof(x)) && return x
    return duplicate_internal(x, IdDict())::typeof(x)
end

"""
Internal implementation of duplicate.

By default delegates to deepcopy, but this can be avoided where it is known to be
unnecessary.
"""
duplicate_internal(x, stackdict::IdDict) = Base.deepcopy_internal(x, stackdict)