text:
# default node types in IIF

SamplableBelief = Union{Distributions.Distribution, KernelDensityEstimate.BallTreeDensity, AliasingScalarSampler}

"""
$(TYPEDEF)

Most basic continuous scalar variable a `::FactorGraph` object.
"""
struct ContinuousScalar <: InferenceVariable
  dims::Int
  labels::Vector{String}
  manifolds::Tuple{Symbol}
  ContinuousScalar(;labels::Vector{<:AbstractString}=String[], manifolds::Tuple{Symbol}=(:Euclid,)) = new(1, labels, manifolds)
end

"""
$(TYPEDEF)

Continuous variable of dimension `.dims` on manifold `.manifolds`.
"""
struct ContinuousMultivariate{T1 <: Tuple} <: InferenceVariable
  dims::Int
  labels::Vector{String}
  manifolds::T1
  ContinuousMultivariate{T}() where {T} = new()
  ContinuousMultivariate{T}(x::Int;labels::Vector{<:AbstractString}=String[], manifolds::T=(:Euclid,)) where {T <: Tuple} = new(x, labels, manifolds)
end

function ContinuousMultivariate(x::Int;
                                labels::Vector{<:AbstractString}=String[],
                                manifolds::T1=(:Euclid,)  )  where {T1 <: Tuple}
  #
  @show maniT = length(manifolds) < x ? ([manifolds[1] for i in 1:x]...,) : manifolds
  ContinuousMultivariate{typeof(maniT)}(x, labels=labels, manifolds=maniT)
end

"""
$(TYPEDEF)

Default prior on all dimensions of a variable node in the factor graph.  `Prior` is
not recommended when non-Euclidean dimensions are used in variables.
"""
struct Prior{T} <: IncrementalInference.FunctorSingleton where T <: SamplableBelief
  Z::T
end
getSample(s::Prior, N::Int=1) = (reshape(rand(s.Z,N),:,N), )

"""
$(TYPEDEF)

Message prior on all dimensions of a variable node in the factor graph.

Notes
- Only temporary existance during CSM operations.
"""
struct MsgPrior{T} <: IncrementalInference.FunctorSingleton where T <: SamplableBelief
  Z::T
  inferdim::Float64
  MsgPrior{T}() where {T} = new{T}()
  MsgPrior{T}(z::T, infd::R) where {T <: SamplableBelief, R <: Real} = new{T}(z, infd)
end
function MsgPrior(z::T, infd::R) where {T <: SamplableBelief, R <: Real}
    MsgPrior{T}(z, infd)
end
getSample(s::MsgPrior, N::Int=1) = (reshape(rand(s.Z,N),:,N), )

struct PackedMsgPrior <: PackedInferenceType where T
  Z::String
  inferdim::Float64
  PackedMsgPrior() = new()
  PackedMsgPrior(z::S, infd::R) where {S <: AbstractString, R <: Real} = new(string(z), infd)
end

function convert(::Type{PackedMsgPrior}, d::MsgPrior)
  PackedMsgPrior(string(d.Z), d.inferdim)
end
function convert(::Type{MsgPrior}, d::PackedMsgPrior)
  MsgPrior(extractdistribution(d.Z), d.inferdim)
end

"""
Converter: Prior -> PackedPrior::Dict{String, Any}
"""
function convert(::Type{Dict{String, Any}}, prior::IncrementalInference.Prior)
    z = convert(Type{Dict{String, Any}}, prior.Z)
    return Packed_Factor([z], "Prior")
end

"""
Converter: PackedPrior::Dict{String, Any} -> Prior
"""
function convert(::Type{Prior}, prior::Dict{String, Any})
    # Genericize to any packed type next.
    z = prior["measurement"][1]
    z = convert(_evalType(z["distType"]), z)
    return Prior(z)
end

"""
$(TYPEDEF)

Partial prior belief (absolute data) on any variable, given `<:SamplableBelief` and which dimensions of the intended variable.
"""
struct PartialPrior{T,P} <: IncrementalInference.FunctorSingleton where {T <: SamplableBelief, P <: Tuple}
  Z::T
  partial::P
end
PartialPrior(z::T, par::P) where {T <: SamplableBelief, P<:Tuple} = PartialPrior{T, P}(z, par)
getSample(s::PartialPrior, N::Int=1) = (reshape(rand(s.Z,N),:,N), )

"""
$(TYPEDEF)

Default linear offset between two scalar variables.
"""
struct LinearConditional{T} <: IncrementalInference.FunctorPairwise where T <: SamplableBelief
  Z::T
end
getSample(s::LinearConditional, N::Int=1) = (reshape(rand(s.Z,N),:,N), )
function (s::LinearConditional)(res::Array{Float64},
                                userdata::FactorMetadata,
                                idx::Int,
                                meas::Tuple{Array{Float64, 2}},
                                X1::Array{Float64,2},
                                X2::Array{Float64,2}  )
  #
  res[1] = meas[1][idx] - (X2[1,idx] - X1[1,idx])
  nothing
end

"""
$(TYPEDEF)

Define a categorical mixture of prior beliefs on a variable.
"""
struct MixturePrior{T} <: IncrementalInference.FunctorSingleton where {T <: SamplableBelief}
  Z::Vector{T}
  C::Distributions.Categorical
  MixturePrior{T}() where T  = new{T}()
  MixturePrior{T}(z::Vector{T}, c::Distributions.Categorical) where {T <: SamplableBelief} = new{T}(z, c)
  MixturePrior{T}(z::Vector{T}, p::Vector{Float64}) where {T <: SamplableBelief} = MixturePrior{T}(z, Distributions.Categorical(p))
end
MixturePrior(z::Vector{T}, c::Union{Distributions.Categorical, Vector{Float64}}) where {T <: SamplableBelief} = MixturePrior{T}(z, c)

getSample(s::MixturePrior, N::Int=1) = (reshape.(rand.(s.Z, N),1,:)..., rand(s.C, N))

"""
$(TYPEDEF)

Define a categorical mixture of (relative) likelihood beliefs between any two variables.
"""
struct MixtureLinearConditional{T} <: IncrementalInference.FunctorPairwise
  Z::Vector{T}
  C::Distributions.Categorical
  MixtureLinearConditional{T}() where T  = new{T}()
  MixtureLinearConditional{T}(z::Vector{T}, c::Distributions.Categorical) where {T <: SamplableBelief} = new{T}(z, c)
  MixtureLinearConditional{T}(z::Vector{T}, p::Vector{Float64}) where {T <: SamplableBelief} = MixtureLinearConditional{T}(z, Distributions.Categorical(p))
end
MixtureLinearConditional(z::Vector{T}, c::Union{Distributions.Categorical, Vector{Float64}}) where {T <: SamplableBelief} = MixtureLinearConditional{T}(z, c)

getSample(s::MixtureLinearConditional, N::Int=1) = (reshape.(rand.(s.Z, N),1,:)..., rand(s.C, N))
function (s::MixtureLinearConditional)(res::Array{Float64},
                               userdata::FactorMetadata,
                               idx::Int,
                               meas::Tuple,
                               X1::Array{Float64,2},
                               X2::Array{Float64,2}  )
  #
  res[1] = meas[meas[end][idx]][idx] - (X2[1,idx] - X1[1,idx])
  nothing
end

## packed types are still developed by hand.  Future versions would likely use a @packable macro to write Protobuf safe versions of factors

"""
$(TYPEDEF)

Serialization type for Prior.
"""
mutable struct PackedPrior <: PackedInferenceType
  Z::String
  PackedPrior() = new()
  PackedPrior(z::AS) where {AS <: AbstractString} = new(z)
end
function convert(::Type{PackedPrior}, d::Prior)
  PackedPrior(string(d.Z))
end
function convert(::Type{Prior}, d::PackedPrior)
  Prior(extractdistribution(d.Z))
end

"""
$(TYPEDEF)

Serialization type for `PartialPrior`.
"""
mutable struct PackedPartialPrior <: PackedInferenceType
  Z::String
  partials::Vector{Int}
  PackedPartialPrior() = new()
  PackedPartialPrior(z::AS) where {AS <: AbstractString} = new(z)
end
function convert(::Type{PackedPartialPrior}, d::PartialPrior)
  PackedPartialPrior(string(d.Z), [d.partial...;])
end
function convert(::Type{PartialPrior}, d::PackedPartialPrior)
  PartialPrior(extractdistribution(d.Z),(d.partials...,))
end

"""
$(TYPEDEF)

Serialization type for `LinearConditional` binary factor.
"""
mutable struct PackedLinearConditional <: PackedInferenceType
  Z::String
  PackedLinearConditional() = new()
  PackedLinearConditional(z::AS) where {AS <: AbstractString} = new(z)
end
function convert(::Type{PackedLinearConditional}, d::LinearConditional)
  PackedLinearConditional(string(d.Z))
end
function convert(::Type{LinearConditional}, d::PackedLinearConditional)
  LinearConditional(extractdistribution(d.Z))
end

"""
$(TYPEDEF)

Serialization type for `MixtureLinearConditional`.
"""
mutable struct PackedMixtureLinearConditional <: PackedInferenceType
  strs::Vector{String}
  cat::String
  PackedMixtureLinearConditional() = new()
  PackedMixtureLinearConditional(z::Vector{<:AbstractString}, cstr::AS) where {AS <: AbstractString} = new(z, cstr)
end
function convert(::Type{PackedMixtureLinearConditional}, d::MixtureLinearConditional)
  PackedMixtureLinearConditional(string.(d.Z), string(d.C))
end
function convert(::Type{MixtureLinearConditional}, d::PackedMixtureLinearConditional)
  MixtureLinearConditional(extractdistribution.(d.strs), extractdistribution(d.cat))
end

"""
$(TYPEDEF)

Serialization type for `MixedPrior`.
"""
mutable struct PackedMixturePrior <: PackedInferenceType
  strs::Vector{String}
  cat::String
  PackedMixturePrior() = new()
  PackedMixturePrior(z::Vector{<:AbstractString}, cstr::AS) where {AS <: AbstractString} = new(z, cstr)
end
function convert(::Type{PackedMixturePrior}, d::MixturePrior)
  PackedMixturePrior(string.(d.Z), string(d.C))
end
function convert(::Type{MixturePrior}, d::PackedMixturePrior)
  MixturePrior(extractdistribution.(d.strs), extractdistribution(d.cat))
end


text:
using Documenter, DiffEqFlux

include("pages.jl")

makedocs(
    sitename = "DiffEqFlux.jl",
    authors="Chris Rackauckas et al.",
    clean = true,
    doctest = false,
    modules = [DiffEqFlux],

    format = Documenter.HTML(analytics = "UA-90474609-3",
                             assets = ["assets/favicon.ico"],
                             canonical="https://diffeqflux.sciml.ai/stable/"),
    pages=pages
)

deploydocs(
   repo = "github.com/SciML/DiffEqFlux.jl.git";
   push_preview = true
)

ext: jl


text:
# statefuliterator
it = ShapleyValues.statefuliterator(1:3)
@test length(it) == 3
@test done(it) == false
@test next(it) == 1
@test next(it) == 2
@test done(it) == false
tmp = next(it)
@test tmp == 3
@test done(it) == true
@test typeof(tmp) == Int64

# allocate_samples
@assert ShapleyValues.allocate_samples([1,1,1,1], 4, ones(Int64, 4)) == ones(4)
@assert ShapleyValues.allocate_samples([1,1,1,1], 8, ones(Int64, 4)) == ones(4)
@assert ShapleyValues.allocate_samples([1,2,1,1], 10, 4*ones(Int64, 4)) == [2,4,2,2]
@assert ShapleyValues.allocate_samples([1,2,1,1], 11, [4,4,4,1]) == [3,4,3,1]
@assert sum(ShapleyValues.allocate_samples([1,2,1,1], 10, 3*ones(Int64, 4))) == 10
@assert maximum(ShapleyValues.allocate_samples([1,2,1,1], 10, 3*ones(Int64, 4))) == 3
@assert minimum(ShapleyValues.allocate_samples([1,2,1,1], 10, 3*ones(Int64, 4))) == 2
@assert sum(ShapleyValues.allocate_samples([1,2,1,1], 21, 10*ones(Int64, 4))) == 21

# MeanVarianceAccumulator
v = [1,2,3.0]
mv = ShapleyValues.MeanVarianceAccumulator()
@test isnan(std(mv))
ShapleyValues.observe!(mv, v[1], 1)
@test mean(mv) == v[1]
@test isnan(std(mv))
ShapleyValues.observe!(mv, v[2], 1)
@test mean(mv) == mean(v[1:2])
@test std(mv) ≈ std(v[1:2])
ShapleyValues.observe!(mv, v[3], 1)
@test mean(mv) == mean(v)
@test std(mv) ≈ std(v)
@test length(mv) == 3

w = [0.1,0.4,0.5]
v = [1,2,3.0]
mv = ShapleyValues.MeanVarianceAccumulator()
@test isnan(std(mv))
ShapleyValues.observe!(mv, v[1], w[1])
@test mean(mv) == v[1]
@test isnan(std(mv))
ShapleyValues.observe!(mv, v[2], w[2])
@test mean(mv) == dot(v[1:2], w[1:2])/sum(w[1:2])
@test sum(((v[1:2] - mean(mv)).^2).*w[1:2])/sum(w[1:2]) * (2/1) ≈ var(mv)
ShapleyValues.observe!(mv, v[3], w[3])
@test mean(mv) ==  dot(v, w)/sum(w)
@test sum(((v - mean(mv)).^2).*w)/sum(w) * (3/2) ≈ var(mv)

# varying_groups
@test ShapleyValues.varying_groups(zeros(3), [0 0 1; 0 1 1; 0 0 1]', Array{Int64,1}[[1,2],[3]]) == [1,2]
@test ShapleyValues.varying_groups(zeros(3), [0 0 1; 0 1 1; 0 0 1]', Array{Int64,1}[[1],[2,3]]) == [2]
@test ShapleyValues.varying_groups(zeros(3), zeros(3,10), Array{Int64,1}[[1],[2,3]]) == Int64[]


text:
#changing the field
Base.convert(::Type{DiffCheck{T, S, N}}, z::DiffCheck{T, S, N}) where {T<:ReComp, S, N} = z
Base.convert(::Type{DiffCheck{T, S, N}}, z::DiffCheck{R, S, N}) where {T, R<:ReComp, S, N} =
	DiffCheck{T, S}(convert(T, value(z)), StaticTuple{T,S,N}(convert.(T, values(derivatives(z)))))

#conversion for scalars
Base.convert(::Type{DiffCheck{T, (), 0}}, x::T) where {T<:ReComp} = DiffCheck{T}(x)
Base.convert(::Type{DiffCheck{T, S, N}}, x::T) where {T<:ReComp, S, N} = convert(DiffCheck{T, S, N}, DiffCheck{T}(x))

Base.convert(::Type{DiffCheck}, x::T) where {T<:ReComp} = Base.convert(DiffCheck{T, (), 0}, x)
Base.convert(::Type{T}, z::DiffCheck) where {T<:ReComp} = (isempty(derivatives(z)) ? convert(T, value(z)) : throw(InexactError(Symbol(T), T, z)))

#This convert is only used to convert upward, i.e., S₂ ⊂ S₁
@generated function Base.convert(::Type{DiffCheck{T₁, S₁, N₁}}, x::DiffCheck{T₂, S₂, N₂}) where {T₁, T₂<:ReComp, S₁, S₂, N₁, N₂}
    if issubset(Set(S₂), Set(S₁))
        str = "("
        for z in S₁
            if z in S₂ # Since Z is the union both cannot be false
               str*="$T₁(x.der[:$z]), "
            else
               str*="zero(T₁) ,"
            end
        end
        str*= ")"
        expr = Meta.parse(str)
        return :(DiffCheck{$T₁, $S₁}($T₁(value(x)), NamedTuple{$S₁, NTuple{$N₁, $T₁}}($expr)))
    else
        return :(throw(InexactError(:convert, Type{DiffCheck{T₁, S₁, N₁}}, x))) #check
    end
end

# Promote rule for
Base.promote_rule(::Type{DiffCheck{T, S, N}}, ::Type{DiffCheck{R, S, N}}) where {T<:ReComp,R<:ReComp, S, N} = DiffCheck{promote_type(T, R), S, N}

# Unsure whether to use this
#@generated function Base.promote_rule(::Type{DiffCheck{T₁, S₁, N₁}}, ::Type{DiffCheck{T₂, S₂, N₂}}) where {T₁, T₂<:ReComp, S₁, S₂, N₁, N₂}
#    Z = Tuple(union(Set(S₁), Set(S₂)))
#    N = length(Z)
#    return :(DiffCheck{promote_type($T₁, $T₂), $Z, $N})
#end

Base.promote_rule(::Type{DiffCheck{T, S, N}}, ::Type{R}) where {T, R<:ReComp, S, N} = DiffCheck{promote_type(T, R), S, N}
Base.promote_rule(::Type{DiffCheck{T, S, N}}, ::Type{T}) where {T<:ReComp, S, N} = DiffCheck{T, S, N}
Base.promote_rule(::Type{DiffCheck{T, (), 0}}, ::Type{T}) where {T<:ReComp, S, N} = DiffCheck{T, (), 0}

Base.widen(::Type{DiffCheck{T, S, N}}) where {T, S, N} = DiffCheck{widen(T),S, N}


text:
"""
    LinearEquationOfState{FT} <: AbstractEquationOfState

Linear equation of state for seawater.
"""
struct LinearEquationOfState{FT} <: AbstractEquationOfState
    thermal_expansion :: FT
    haline_contraction :: FT
end

Base.summary(eos::LinearEquationOfState) =
    string("LinearEquationOfState(thermal_expansion=", scalar_summary(eos.thermal_expansion),
                               ", haline_contraction=", scalar_summary(eos.haline_contraction), ")")

Base.show(io, eos::LinearEquationOfState) = print(io, summary(eos))

"""
    LinearEquationOfState([FT=Float64;] thermal_expansion=1.67e-4, haline_contraction=7.80e-4)

Return `LinearEquationOfState` for `SeawaterBuoyancy` with
`thermal_expansion` coefficient and `haline_contraction` coefficient.
The buoyancy perturbation ``b`` for `LinearEquationOfState` is

```math
    b = g (α T - β S),
```

where ``g`` is gravitational acceleration, ``α`` is `thermal_expansion`, ``β`` is
`haline_contraction`, ``T`` is temperature, and ``S`` is practical salinity units.

Default constants in units inverse Kelvin and practical salinity units
for `thermal_expansion` and `haline_contraction`, respectively,
are taken from Table 1.2 (page 33) of Vallis, "Atmospheric and Oceanic Fluid
Dynamics: Fundamentals and Large-Scale Circulation" (2nd ed, 2017).
"""
LinearEquationOfState(FT=Float64; thermal_expansion=1.67e-4, haline_contraction=7.80e-4) =
    LinearEquationOfState{FT}(thermal_expansion, haline_contraction)

#####
##### Thermal expansion and haline contraction coefficients
#####

@inline  thermal_expansion(Θ, sᴬ, D, eos::LinearEquationOfState) = eos.thermal_expansion
@inline haline_contraction(Θ, sᴬ, D, eos::LinearEquationOfState) = eos.haline_contraction

# Shortcuts
@inline  thermal_expansionᶜᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.thermal_expansion
@inline  thermal_expansionᶠᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.thermal_expansion
@inline  thermal_expansionᶜᶠᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.thermal_expansion
@inline  thermal_expansionᶜᶜᶠ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.thermal_expansion

@inline haline_contractionᶜᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.haline_contraction
@inline haline_contractionᶠᶜᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.haline_contraction
@inline haline_contractionᶜᶠᶜ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.haline_contraction
@inline haline_contractionᶜᶜᶠ(i, j, k, grid, eos::LinearEquationOfState, C) = eos.haline_contraction

#####
##### Convinient aliases to dispatch on
#####

const LinearSeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState} where FT
const LinearTemperatureSeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState, <:Nothing, <:Number} where FT
const LinearSalinitySeawaterBuoyancy = SeawaterBuoyancy{FT, <:LinearEquationOfState, <:Number, <:Nothing} where FT

#####
##### BuoyancyModels perturbation
#####

@inline buoyancy_perturbation(i, j, k, grid, b::LinearSeawaterBuoyancy, C) =
    @inbounds b.gravitational_acceleration * (b.equation_of_state.thermal_expansion * C.T[i, j, k] -
                                              b.equation_of_state.haline_contraction * C.S[i, j, k])

@inline buoyancy_perturbation(i, j, k, grid, b::LinearTemperatureSeawaterBuoyancy, C) =
    @inbounds b.gravitational_acceleration * b.equation_of_state.thermal_expansion * C.T[i, j, k]

@inline buoyancy_perturbation(i, j, k, grid, b::LinearSalinitySeawaterBuoyancy, C) =
    @inbounds - b.gravitational_acceleration * b.equation_of_state.haline_contraction * C.S[i, j, k]


text:
############################################################################################
# Proposal steps and post processing
#=
iter=1
_resample = resamplemethods[1]
_reference = referencemethods[1]
generated = generating[2]
=#
for iter in eachindex(objectives)
    @testset "Kernel construction and propose, all models" begin
        ## Resampling methods
        for _resample in resamplemethods
            ## Referencing methods
            for _reference in referencemethods
                for generated in generating
                    _obj = deepcopy(objectives[iter])
                    ## Define PF default tuning parameter
                    pfdefault = ParticleFilterDefault(;
                        referencing = _reference,
                        resampling = _resample,
                        generated = generated
                    )
                    ## Check if Default options work
                    ParticleFilter(_rng, _obj, pfdefault)
                    ParticleFilter(_rng, _obj, pfdefault, SampleDefault())
                    ## Check if we can initiate from Constructor
                    constructor = ParticleFilterConstructor(:latent, pfdefault)
                    constructor(_rng, _obj.model, _obj.data, 1., SampleDefault())
                    ParticleFilter((:latent,))
                    ## Initialize kernel and check if it can be run
                    pfkernel = ParticleFilter(
                        _rng,
                        _obj,
                        pfdefault
                    )
                    _val, _diag = propose(_rng, pfkernel, _obj)
                    @test size(pfkernel.particles.val, 2) == length(_obj.data)
                    ## Check if all particles are correct ~ Easy to check for Semi-Markov particles
                    if _obj.model.id isa HSMM
                        @test check_correctness(pfkernel.particles.kernel, pfkernel.particles.val) == 0
                    end
                    ## Check if all ancestors are correct
                    @test check_ancestors(pfkernel.particles.ancestor)
                    ## Postprocessing
                    BaytesFilters.generate_showvalues(_diag)()
                    _type = infer(_rng, BaytesFilters.AbstractDiagnostics, pfkernel, _obj.model, _obj.data)
                    @test _diag isa _type
                    _type = infer(_rng, pfkernel.tune, pfkernel.particles.kernel, _obj.model, _obj.data)
                    @test _diag.base.prediction isa _type
                    _type = infer(_rng, pfkernel, _obj.model, _obj.data)
                    @test _diag.base.prediction isa _type
                    _type = BaytesFilters.infer_generated(_rng, pfkernel, _obj.model, _obj.data)
                    @test _diag.generated isa _type
                    @test predict(_rng, pfkernel, _obj) isa typeof(_diag.base.prediction)
                end
            end
        end
    end
end

############################################################################################
# Propagation steps
for iter in eachindex(objectives)
    Ndata_extension = [1, 10]
    @testset "Kernel propagation, all models" begin
        for newdata in Ndata_extension
        # First extend 1 data point, then 10
            # -> once PF will only extend Ndata, then Nparticles and Ndata
            data2 = vcat(objectives[iter].data, objectives[iter].data[1:newdata]) #randn(length(objectives[iter].data)+10)
            ## Resampling methods
            for _resample in resamplemethods
                ## Referencing methods
                _reference = Marginal()
                _obj = deepcopy(objectives[iter])
                ## Define PF default tuning parameter
                pfdefault = ParticleFilterDefault(;
                    #!NOTE: So Nparticles does not change if 1 more data point is added
                    coverage = .5,
                    referencing = _reference,
                    resampling = _resample
                )
                ## Initialize kernel and check if it can be run
                pfkernel = ParticleFilter(
                    _rng,
                    _obj,
                    pfdefault
                )
                propose(_rng, pfkernel, _obj)
                ## Propagate forward
                propagate!(_rng, pfkernel, _obj.model, data2)
                @test size(pfkernel.particles.val, 2) == length(data2)
                propose!(_rng, pfkernel, _obj.model, data2)
                ## Check if Nparticles change for more data accordingly even if not propagated
                pfkernel2 = ParticleFilter(
                    _rng,
                    _obj,
                    pfdefault
                )
                ChainsInit = pfkernel2.tune.chains.Nchains
                @test ChainsInit == size(pfkernel2.particles.val,1)
                propose!(_rng, pfkernel2, deepcopy(objectives[iter].model), data2)
                @test size(pfkernel2.particles.val,1) <= size(data2,1)*pfkernel2.tune.chains.coverage
                ## Check if Nparticles change for less data accordingly even if not propagated
                propose!(_rng, pfkernel2, deepcopy(objectives[iter].model), data2[1:Int(round(length(data2)/2))])
                @test size(pfkernel2.particles.val,1) <= size(data2[1:Int(round(length(data2)/2))],1)*pfkernel2.tune.chains.coverage
            end
        end
    end
end


text:
# Computing some distributions for discrete values
using Distributions

function make_distributions(c_min,c_max,number,N)

values = linspace(c_min,c_max,number)

    return sample(values,N)

end


text:
using Test

include("reduce.jl")

if capability(device()) < v"3.0"
    @warn("this example requires a newer GPU")
    exit(0)
end

len = 10^7
input = ones(Int32, len)
output = similar(input)

# CPU
cpu_val = reduce(+, input)

# CUDAnative
let
    gpu_input = CuTestArray(input)
    gpu_output = CuTestArray(output)
    gpu_reduce(+, gpu_input, gpu_output)
    gpu_val = Array(gpu_output)[1]
    @assert cpu_val == gpu_val
end


text:
include("_tile_coder.jl")
"""
    TileCoder(num_tilings, num_tiles, num_features, num_ints)
Tile coder for coding all features together.
"""
mutable struct TileCoder <: AbstractFeatureConstructor
    # Main Arguments
    tilings::Int64
    tiles::Int64
    dims::Int64
    ints::Int64

    # Optional Arguments
    wrap::Bool
    wrapwidths::Float64

    iht::TileCoding.IHT
    TileCoder(num_tilings, num_tiles, num_features, num_ints=1; wrap=false, wrapwidths=0.0) =
        new(num_tilings, num_tiles, num_features, num_ints, wrap, 0.0, TileCoding.IHT(num_tilings*(num_tiles+1)^num_features * num_ints))
end

function create_features(fc::TileCoder, s; ints=[], readonly=false)
    if fc.wrap
        return 1 .+ TileCoding.tileswrap!(fc.iht, fc.tilings, s.*fc.tiles, fc.wrapwidths, ints, readonly)
    else
        return 1 .+ TileCoding.tiles!(fc.iht, fc.tilings, s.*fc.tiles, ints, readonly)
    end
end

feature_size(fc::TileCoder) = fc.tilings*(fc.tiles+1)^fc.dims * num_ints


text:
export isPowerOfTwo, isMersenneNumber, getMersenneNumber

function isPowerOfTwo(n::Integer)::Bool
    return (n != 0) && (n & (n - 1)) == 0
end

function isMersenneNumber(n::Integer)::Bool
    return isPowerOfTwo(n + 1)
end

function getMersenneNumber(n::Integer)::Integer
    return 2^n - 1
end

max_stars_repo_name: EllianCarlos/number-theory-but-in-julia

max_issues_repo_name: EllianCarlos/number-theory-but-in-julia

max_forks_repo_name: EllianCarlos/number-theory-but-in-julia


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
hist(res, bins=30, density=true, color="green")

xx = range(-1100, 1100, length=100)
ndfit = D.fit_mle(D.Normal, res)
lfit  = D.fit_mle(D.Laplace, res)

plot(xx, pdf.(ndfit, xx), lw=3, color="orange", label="Normal fit")
plot(xx, pdf.(lfit, xx), lw=3, color="magenta", label="Laplace fit")

legend(fontsize=12)

xticks(fontsize=12); yticks(fontsize=12)
xlabel("Residual (ŷ - y)", fontsize=14); ylabel("Density", fontsize=14)
xlim([-1100, 1100])

savefig(joinpath(@OUTPUT, "ISL-lab-6-g4.svg")) # hide


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

using Markdown
using InteractiveUtils

# ╔═╡ 5057c590-c5a4-11eb-2621-4b67f12ca95d
md"""
# K Nearest Neighbors
"""

# ╔═╡ 74f4afa1-dd5a-4931-ab48-a1d4fc77bbe5
md"""
[Reference 1](https://scikit-learn.org/stable/modules/neighbors.html)\
[Reference 2](https://en.wikipedia.org/wiki/Nearest_neighbor_search)\
[Reference 3](https://booking.ai/k-nearest-neighbours-from-slow-to-fast-thanks-to-maths-bec682357ccd)

### Naive Approach

No transforming on original dataset (No training)\
For prediction, iter through the original dataset, and find the nearest K data, by a given metrics

Can use euclidean distance:

$$\text{dist}(X_1, X_2) = \|X_1 - X_2\|$$

A better approach is cosine similarity:

$$\text{sim}(X_1, X_2) = \frac{X_1 \cdot X_2}{\|X_1\| \|X_2\|}$$

which computes the cos value between two vectors\
1 for 0 degree, and less than 1 for $(0, \pi]$
"""

# ╔═╡ cb1dd9ac-ca78-4169-b073-943f64160dbc
module tools include("../tools.jl") end

# ╔═╡ f1fa76da-f656-4f12-97d7-a8b847ca0733
JuTools = tools.JuTools

# ╔═╡ c3ad6021-8080-4870-8350-1337d548f23a
import Statistics

# ╔═╡ 83864e6e-4b62-45d5-962e-39e36bbcc43d
import Random

# ╔═╡ 0d7f31a9-7ae2-4be5-81b2-313ac58af62d
import LinearAlgebra

# ╔═╡ 73b4e7e0-7294-46dc-9545-9728277df652
function cosine_sim(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number)::AbstractFloat
    @assert size(X1) == size(X2)
    @assert ndims(X1) == ndims(X2) == 1
    product = LinearAlgebra.dot(X1, X2)
    X1_norm = LinearAlgebra.norm(X1, 2)
    X2_norm = LinearAlgebra.norm(X2, 2)
    return product / (X1_norm * X2_norm)
end

# ╔═╡ 19b1c000-6fa3-42b7-9fc7-983c7caf44b9
begin
X_data, Y_data = JuTools.data_generate_linear_2d()
size(X_data), size(Y_data)
end

# ╔═╡ 063389c7-07b6-4050-9241-ea953c4f4338
X_data[1:2, :]

# ╔═╡ cf291226-c0a1-4206-a8cb-53bedb67b139
cosine_sim(X_data[1, :], X_data[2, :])

# ╔═╡ fbf7a012-dcf8-498f-ad29-ec8c4854406c
cosine_sim(X_data[1, :], X_data[3, :])

# ╔═╡ 2cf4e990-f788-44e3-94cc-9e210e47c762
# define majority vote function
function majority_vote(y::Array{T} where T<:Number)::Number
    @assert ndims(y) == 1
    unique_votes = Dict{Number, Integer}()
    for y_val in y
        if !haskey(unique_votes, y_val)
            push!(unique_votes, y_val => 1)
        else
            unique_votes[y_val] += 1
        end
    end
    result = sort(collect(unique_votes), by=m->m[2])
    return result[end][1]
end

# ╔═╡ c9d27318-5e3d-4c9c-9ba5-556c649baaf9
majority_vote([1,1,0])

# ╔═╡ ba24b3c5-f350-4774-a018-9acf1f4597ad
md"""
Output ordering is affected by input ordering
"""

# ╔═╡ 5e564f3b-1bda-41b1-8b22-bdf6d91314d2
majority_vote([1,1,0,0])

# ╔═╡ ed575b52-7509-4140-b5e7-295920f73e3d
majority_vote([1,1,0,0,0])

# ╔═╡ 8406faec-66f4-46aa-9582-964ebdbecfba
begin
X_train, X_test, Y_train, Y_test = JuTools.split_data(X_data, Y_data, shuffle=true, ratio=0.3)
size(X_train), size(X_test), size(Y_train), size(Y_test)
end

# ╔═╡ 3001da4a-3288-4168-b514-c0b12f552fdd
# define predict function, naive approach
function predict_naive(X_predict::Array{T} where T<:Number, K::Integer, X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number)::Array
    @assert ndims(X_data) == 2
    @assert ndims(Y_data) == 1
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert 0 < ndims(X_predict) <= 2
    @assert 0 < K < size(X_data)[1]
    if ndims(X_predict) < 2
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert size(X_predict)[2] == size(X_data)[2]
    result = Array{Number}(undef, size(X_predict)[1])
    sim = Array{Tuple{Integer, AbstractFloat}}(undef, size(X_data)[1])
    for i in 1:size(X_predict)[1]
        vec_predict = X_predict[i, :]
        for j in 1:size(X_data)[1]
            vec_data = X_data[j, :]
            vec_similarity = cosine_sim(vec_predict, vec_data)
            sim[j] = (j, vec_similarity)
        end
        sort!(sim, by=m->m[2], rev=true)
        K_nearest_votes = Y_data[[m[1] for m in sim[1:K]]]
        result[i] = majority_vote(K_nearest_votes)
    end
    return result
end

# ╔═╡ 30320c6e-7ada-41b6-a938-f1fd19cc120a
Y_predict = predict_naive(X_test, 5, X_train, Y_train)

# ╔═╡ 27b26afb-a91c-47c4-94ef-1d64761f291d
JuTools.compute_accuracy(Y_predict, Y_test)

# ╔═╡ 0e1b5c16-f444-41dd-9e44-999969ba647c
# what about dist similarity?
function dist_sim(X1::Array{T} where T<:Number, X2::Array{T} where T<:Number)::AbstractFloat
    @assert size(X1) == size(X2)
    @assert ndims(X1) == ndims(X2) == 1
    return sqrt(sum((X1 .- X2).^2))
end

# ╔═╡ c42f2b8a-a88d-4e8c-a1c0-1a371ae6a669
dist_sim(X_data[1, :], X_data[2, :])

# ╔═╡ 158f77bf-2eb8-490a-a062-6364a6ec70e9
dist_sim(X_data[2, :], X_data[3, :])

# ╔═╡ 323e10a0-c68f-433a-afc5-0ac2514a90d1
md"""
It's greatly affected by the scale of data!
"""

# ╔═╡ 52010e3f-2e4f-4778-b53e-4fad04dd75bf
function predict_naive_fun(X_predict::Array{T} where T<:Number, K::Integer, X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number)::Array
    @assert ndims(X_data) == 2
    @assert ndims(Y_data) == 1
    @assert size(X_data)[1] == size(Y_data)[1]
    @assert 0 < ndims(X_predict) <= 2
    @assert 0 < K < size(X_data)[1]
    if ndims(X_predict) < 2
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert size(X_predict)[2] == size(X_data)[2]
    result = Array{Number}(undef, size(X_predict)[1])
    sim = Array{Tuple{Integer, AbstractFloat}}(undef, size(X_data)[1])
    for i in 1:size(X_predict)[1]
        vec_predict = X_predict[i, :]
        for j in 1:size(X_data)[1]
            vec_data = X_data[j, :]
            vec_similarity = dist_sim(vec_predict, vec_data)
            sim[j] = (j, vec_similarity)
        end
        sort!(sim, by=m->m[2])
        K_nearest_votes = Y_data[[m[1] for m in sim[1:K]]]
        result[i] = majority_vote(K_nearest_votes)
    end
    return result
end

# ╔═╡ 36ceae2a-4c04-4f67-9a6d-19a44cb2fa5b
JuTools.compute_accuracy(predict_naive_fun(X_test, 5, X_train, Y_train), Y_test)

# ╔═╡ 2c7f7263-5ba0-4786-8b68-8ca8fb3ffc26
md"""
Although it (`predict_naive`) may be slow on large dataset, it is easy to implement and it works as expected
"""

# ╔═╡ b2aac2a9-cf86-44f0-ae97-b2a52327b0f9
md"""
### K-Dimensional Tree (K-d tree) Approach

A space partitioning technique\
Treat each data row as a point in `k`-dimensional space\
[Wikipedia](https://en.wikipedia.org/wiki/K-d_tree)
"""

# ╔═╡ b127297e-cb0c-46f5-8551-3ba2c4dae4c6
mutable struct KdTree
    X_data::Array{T} where T<:Number # 1d vector
    Y_data::Number                   # number
    child_l::Union{KdTree,Nothing}
    child_r::Union{KdTree,Nothing}
end

# ╔═╡ c3c4fe2c-6e92-4c98-bce3-e77f394c405d
# K-d tree generator function
function create_kdtree(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number)::KdTree
    @assert ndims(X_data) == 2
    @assert ndims(Y_data) == 1
    @assert size(X_data)[1] == size(Y_data)[1]
    function kdtree_recursive_generate(X_data::Array, Y_data::Array, depth::Integer, n_axes::Integer)::KdTree
        curr_axis = mod(depth, n_axes) + 1 # array starts from 1
        data_combined = hcat(X_data, Y_data)
        data_combined = sortslices(data_combined, by=m->m[curr_axis], dims=1)
        X_data = data_combined[:, 1:end-1]
        Y_data = data_combined[:, end]
        i_mid = div(size(X_data)[1], 2) + 1
        node_X_data = X_data[i_mid, :]
        node_Y_data = Y_data[i_mid]
        node = KdTree(node_X_data, node_Y_data, nothing, nothing)
        if i_mid > 1
            node.child_l = kdtree_recursive_generate(X_data[1:i_mid-1,:], Y_data[1:i_mid-1], depth+1, n_axes)
        end
        if i_mid < size(X_data)[1]
            node.child_r = kdtree_recursive_generate(X_data[i_mid+1:end,:], Y_data[i_mid+1:end], depth+1, n_axes)
        end
        return node
    end
    return kdtree_recursive_generate(X_data, Y_data, 0, size(X_data)[2])
end

# ╔═╡ a76dcf16-060e-4ece-9b8d-76c45deb38c4
kdtree = create_kdtree(X_data, Y_data)

# ╔═╡ 4c3444dd-f2cb-4954-8dad-d2e9466252e8
kdtree_test = create_kdtree(reshape([30,5,10,70,50,35], (6, 1)), [1,1,1,1,1,1])

# ╔═╡ 94c43601-cff7-4c7f-a116-8a82aefed5b0
# inspired from https://stackoverflow.com/questions/1627305/nearest-neighbor-k-d-tree-wikipedia-proof/37107030#37107030
# note that for kdtree search, we use euclidean distance
function predict_kdtree(X_predict::Array{T} where T<:Number, kdtree::KdTree; K::Integer=5)::Array
    @assert K > 0
    @assert 0 < ndims(X_predict) <= 2
    if ndims(X_predict) == 1
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert size(X_predict)[2] == size(kdtree.X_data)[1]

    function kdtree_closest_max(kdtree_closest::Array{Union{KdTree, Nothing}},
                kdtree_closest_val::Array{AbstractFloat})::Tuple{Integer, AbstractFloat}
        default = (0, 0.0)
        for i in 1:size(kdtree_closest)[1]
            if kdtree_closest[i] === nothing
                break
            elseif default[1] == 0 || (kdtree_closest_val[i] > default[2])
                default = (i, kdtree_closest_val[i])
            end
        end
        return default
    end

    function kdtree_update_nearest!(X_vec::Array, kdtree::KdTree, kdtree_closest::Array{Union{KdTree, Nothing}},
                kdtree_closest_val::Array{AbstractFloat})
        @assert size(kdtree_closest) == size(kdtree_closest_val)
        distance = dist_sim(kdtree.X_data, X_vec)
        if nothing in kdtree_closest
            for i in 1:size(kdtree_closest)[1]
                if kdtree_closest[i] === nothing
                    kdtree_closest[i] = KdTree(kdtree.X_data, kdtree.Y_data, nothing, nothing)
                    kdtree_closest_val[i] = distance
                    break
                end
            end
        else
            curr_max = kdtree_closest_max(kdtree_closest, kdtree_closest_val)
            if distance < curr_max[2]
                kdtree_closest[curr_max[1]] = KdTree(kdtree.X_data, kdtree.Y_data, nothing, nothing)
                kdtree_closest_val[curr_max[1]] = distance
            end
        end
    end

    function kdtree_recursive_search!(X_vec::Array, kdtree::KdTree, depth::Integer, n_axes::Integer,
                kdtree_closest::Array{Union{KdTree, Nothing}}, kdtree_closest_val::Array{AbstractFloat})
        @assert size(kdtree_closest) == size(kdtree_closest_val)
        # check current node
        kdtree_update_nearest!(X_vec, kdtree, kdtree_closest, kdtree_closest_val)
        # run on children
        curr_axis = mod(depth, n_axes) + 1 # array starts from 1
        if X_vec[curr_axis] < kdtree.X_data[curr_axis]
            if kdtree.child_l !== nothing
                kdtree_recursive_search!(X_vec, kdtree.child_l, depth+1, n_axes, kdtree_closest, kdtree_closest_val)
            end
            if (X_vec[curr_axis] + kdtree_closest_max(kdtree_closest, kdtree_closest_val)[2] >= kdtree.X_data[curr_axis]) && kdtree.child_r !== nothing
                kdtree_recursive_search!(X_vec, kdtree.child_r, depth+1, n_axes, kdtree_closest, kdtree_closest_val)
            end
        else
            if kdtree.child_r !== nothing
                kdtree_recursive_search!(X_vec, kdtree.child_r, depth+1, n_axes, kdtree_closest, kdtree_closest_val)
            end
            if (X_vec[curr_axis] - kdtree_closest_max(kdtree_closest, kdtree_closest_val)[2] <= kdtree.X_data[curr_axis]) && kdtree.child_l !== nothing
                kdtree_recursive_search!(X_vec, kdtree.child_l, depth+1, n_axes, kdtree_closest, kdtree_closest_val)
            end
        end
    end

    result = Array{Number}(undef, size(X_predict)[1])
    for i in 1:size(X_predict)[1]
        kdtree_closest = Array{Union{KdTree, Nothing}}(nothing, K)
        kdtree_closest_val = Array{AbstractFloat}(undef, K)
        kdtree_recursive_search!(X_predict[i, :], kdtree, 0, size(X_predict)[2], kdtree_closest, kdtree_closest_val)
        K_nearest_votes = Number[]
        for i in 1:K
            if kdtree_closest[i] === nothing
                break
            else
                push!(K_nearest_votes, kdtree_closest[i].Y_data)
            end
        end
        result[i] = majority_vote(K_nearest_votes)
    end
    return result
end

# ╔═╡ 569da607-dea8-4994-9ffe-05959071a855
kdtree_train = create_kdtree(X_train, Y_train)

# ╔═╡ 0ea4cb86-6be5-49f1-b9b8-56f6a8cc8ed0
JuTools.compute_accuracy(predict_kdtree(X_test, kdtree_train, K=10), Y_test)

# ╔═╡ 98be952f-0591-45c6-be83-3c0f6cc5d26b
md"""
### Ball Tree Approach

A better space partition approach\
More efficient than K-d Tree when searching\
[Wikipedia](https://en.wikipedia.org/wiki/Ball_tree)
"""

# ╔═╡ dc218b5e-73f5-4dc8-b523-62bed53a71e7
mutable struct BallTree
    X_data::Array{T} where T<:Number
    Y_data::Number
    pivot::Union{Array{T},Nothing} where T<:Number # defines pivot point of hypersphere
    radius::AbstractFloat                            # defines radius of hypersphere
    child_l::Union{BallTree,Nothing}
    child_r::Union{BallTree,Nothing}
end

# ╔═╡ d7cc4b1c-f953-4f64-95bb-09e1d51a38f5
# inspired from https://gist.github.com/jakevdp/5216193
# ball tree generator function
function create_balltree(X_data::Array{T} where T<:Number, Y_data::Array{T} where T<:Number)::BallTree
    @assert ndims(X_data) == 2
    @assert ndims(Y_data) == 1
    @assert size(X_data)[1] == size(Y_data)[1]
    balltree = nothing
    if size(X_data)[1] == 1
        balltree = BallTree(X_data[1, :], Y_data[1], nothing, 0.0, nothing, nothing)
    else
        # find pivot
        pivot = vec(sum(X_data, dims=1)) ./ size(X_data)[1]
        # find radius
        radius = 0.0
        for i in 1:size(X_data)[1]
            X_vec = X_data[i, :]
            dist = dist_sim(pivot, X_vec)
            if dist > radius
                radius = dist
            end
        end
        # find greatest spread dimension
        d_greatest_spread = 1
        n_spread = 0.0
        for i in 1:size(X_data)[2]
            X_vec = X_data[:, i]
            current_spread = abs(maximum(X_vec) - minimum(X_vec))
            if current_spread > n_spread
                d_greatest_spread = i
                n_spread = current_spread
            end
        end
        data_combined = hcat(X_data, Y_data)
        data_combined = sortslices(data_combined, by=m->m[d_greatest_spread], dims=1)
        X_data = data_combined[:, 1:end-1]
        Y_data = data_combined[:, end]
        i_mid = div(size(X_data)[1], 2) + 1
        node_X_data = X_data[i_mid, :]
        node_Y_data = Y_data[i_mid]
        balltree = BallTree(node_X_data, node_Y_data, pivot, radius, nothing, nothing)
        if i_mid > 1
            balltree.child_l = create_balltree(X_data[1:i_mid-1,:], Y_data[1:i_mid-1])
        end
        if i_mid < size(X_data)[1]
            balltree.child_r = create_balltree(X_data[i_mid+1:end,:], Y_data[i_mid+1:end])
        end
    end
    return balltree
end

# ╔═╡ ffb15055-c887-47a9-88d7-71a33a3c2dba
balltree = create_balltree(X_data, Y_data)

# ╔═╡ 77ed169b-d426-4f3f-8362-809cb0c8dd22
balltree.X_data

# ╔═╡ 79ddb157-8ee4-4c44-b8bf-b4f36f88acce
balltree.Y_data

# ╔═╡ 6c589940-e00a-40c9-b2d9-c10384e4ac62
balltree.pivot

# ╔═╡ 5b5c5f78-75c2-4a85-b918-e7485f492d9e
balltree.radius

# ╔═╡ 75368c67-8774-4f67-9163-d31475f4e289
balltree_test = create_balltree(reshape([30,5,10,70,50,35], (6, 1)), [1,1,1,1,1,1])

# ╔═╡ a14bd64d-4188-4682-8a47-c3999afb0d09
# search function for ball tree
# similar to kdtree search
function predict_balltree(X_predict::Array{T} where T<:Number, balltree::BallTree; K::Integer=5)::Array
    @assert K > 0
    @assert 0 < ndims(X_predict) <= 2
    if ndims(X_predict) == 1
        X_predict = reshape(X_predict, (1, size(X_predict)[1]))
    end
    @assert size(X_predict)[2] == size(balltree.X_data)[1]

    function balltree_closest_max(balltree_closest::Array{Union{BallTree, Nothing}},
                balltree_closest_val::Array{AbstractFloat})::Tuple{Integer, AbstractFloat}
        # find the maximum value in balltree_closest_val
        # return its index and value
        default = (0, 0.0)
        for i in 1:size(balltree_closest)[1]
            if balltree_closest[i] === nothing
                break
            elseif default[1] == 0 || (balltree_closest_val[i] > default[2])
                default = (i, balltree_closest_val[i])
            end
        end
        return default
    end

    function balltree_update_nearest!(X_vec::Array, balltree::BallTree, balltree_closest::Array{Union{BallTree, Nothing}},
                balltree_closest_val::Array{AbstractFloat})
        # update current node by distance
        @assert size(balltree_closest) == size(balltree_closest_val)
        distance = dist_sim(balltree.X_data, X_vec)
        if nothing in balltree_closest
            for i in 1:size(balltree_closest)[1]
                if balltree_closest[i] === nothing
                    balltree_closest[i] = balltree
                    balltree_closest_val[i] = distance
                    break
                end
            end
        else
            curr_max = balltree_closest_max(balltree_closest, balltree_closest_val)
            if distance < curr_max[2]
                balltree_closest[curr_max[1]] = balltree
                balltree_closest_val[curr_max[1]] = distance
            end
        end
    end

    function balltree_recursive_search!(X_vec::Array, balltree::BallTree, balltree_closest::Array{Union{BallTree, Nothing}},
            balltree_closest_val::Array{AbstractFloat})
        # recursively search a balltree for K nearest neighbors
        @assert size(balltree_closest) == size(balltree_closest_val)
        if (!(nothing in balltree_closest) && (balltree.pivot !== nothing)
                && (dist_sim(X_vec, balltree.pivot) - balltree.radius >=
                    balltree_closest_max(balltree_closest, balltree_closest_val)[2]))
            return nothing
        end
        # check current node
        balltree_update_nearest!(X_vec, balltree, balltree_closest, balltree_closest_val)
        # run on children
        if balltree.child_l === nothing || balltree.child_r === nothing
            if balltree.child_l !== nothing
                balltree_recursive_search!(X_vec, balltree.child_l, balltree_closest, balltree_closest_val)
            end
            if balltree.child_r !== nothing
                balltree_recursive_search!(X_vec, balltree.child_r, balltree_closest, balltree_closest_val)
            end
        else
            dist_left = dist_sim(X_vec, balltree.child_l.X_data)
            dist_right = dist_sim(X_vec, balltree.child_r.X_data)
            if dist_left < dist_right
                balltree_recursive_search!(X_vec, balltree.child_l, balltree_closest, balltree_closest_val)
                balltree_recursive_search!(X_vec, balltree.child_r, balltree_closest, balltree_closest_val)
            else
                balltree_recursive_search!(X_vec, balltree.child_r, balltree_closest, balltree_closest_val)
                balltree_recursive_search!(X_vec, balltree.child_l, balltree_closest, balltree_closest_val)
            end
        end
    end

    result = Array{Number}(undef, size(X_predict)[1])
    for i in 1:size(X_predict)[1]
        balltree_closest = Array{Union{BallTree, Nothing}}(nothing, K)
        balltree_closest_val = Array{AbstractFloat}(undef, K)
        balltree_recursive_search!(X_predict[i, :], balltree, balltree_closest, balltree_closest_val)
        K_nearest_votes = Number[]
        for i in 1:K
            if balltree_closest[i] === nothing
                break
            else
                push!(K_nearest_votes, balltree_closest[i].Y_data)
            end
        end
        result[i] = majority_vote(K_nearest_votes)
    end
    return result
end

# ╔═╡ 06e6d38d-b08e-4583-863a-e6e06726143b
balltree_train = create_balltree(X_train, Y_train)

# ╔═╡ 41e118e2-2172-4910-ba52-12e199737c6a
JuTools.compute_accuracy(predict_balltree(X_test, balltree_train, K=20), Y_test)

# ╔═╡ Cell order:
# ╟─5057c590-c5a4-11eb-2621-4b67f12ca95d
# ╟─74f4afa1-dd5a-4931-ab48-a1d4fc77bbe5
# ╠═cb1dd9ac-ca78-4169-b073-943f64160dbc
# ╠═f1fa76da-f656-4f12-97d7-a8b847ca0733
# ╠═c3ad6021-8080-4870-8350-1337d548f23a
# ╠═83864e6e-4b62-45d5-962e-39e36bbcc43d
# ╠═0d7f31a9-7ae2-4be5-81b2-313ac58af62d
# ╠═73b4e7e0-7294-46dc-9545-9728277df652
# ╠═19b1c000-6fa3-42b7-9fc7-983c7caf44b9
# ╠═063389c7-07b6-4050-9241-ea953c4f4338
# ╠═cf291226-c0a1-4206-a8cb-53bedb67b139
# ╠═fbf7a012-dcf8-498f-ad29-ec8c4854406c
# ╠═2cf4e990-f788-44e3-94cc-9e210e47c762
# ╠═c9d27318-5e3d-4c9c-9ba5-556c649baaf9
# ╟─ba24b3c5-f350-4774-a018-9acf1f4597ad
# ╠═5e564f3b-1bda-41b1-8b22-bdf6d91314d2
# ╠═ed575b52-7509-4140-b5e7-295920f73e3d
# ╠═8406faec-66f4-46aa-9582-964ebdbecfba
# ╠═3001da4a-3288-4168-b514-c0b12f552fdd
# ╠═30320c6e-7ada-41b6-a938-f1fd19cc120a
# ╠═27b26afb-a91c-47c4-94ef-1d64761f291d
# ╠═0e1b5c16-f444-41dd-9e44-999969ba647c
# ╠═c42f2b8a-a88d-4e8c-a1c0-1a371ae6a669
# ╠═158f77bf-2eb8-490a-a062-6364a6ec70e9
# ╟─323e10a0-c68f-433a-afc5-0ac2514a90d1
# ╠═52010e3f-2e4f-4778-b53e-4fad04dd75bf
# ╠═36ceae2a-4c04-4f67-9a6d-19a44cb2fa5b
# ╟─2c7f7263-5ba0-4786-8b68-8ca8fb3ffc26
# ╟─b2aac2a9-cf86-44f0-ae97-b2a52327b0f9
# ╠═b127297e-cb0c-46f5-8551-3ba2c4dae4c6
# ╠═c3c4fe2c-6e92-4c98-bce3-e77f394c405d
# ╠═a76dcf16-060e-4ece-9b8d-76c45deb38c4
# ╠═4c3444dd-f2cb-4954-8dad-d2e9466252e8
# ╠═94c43601-cff7-4c7f-a116-8a82aefed5b0
# ╠═569da607-dea8-4994-9ffe-05959071a855
# ╠═0ea4cb86-6be5-49f1-b9b8-56f6a8cc8ed0
# ╟─98be952f-0591-45c6-be83-3c0f6cc5d26b
# ╠═dc218b5e-73f5-4dc8-b523-62bed53a71e7
# ╠═d7cc4b1c-f953-4f64-95bb-09e1d51a38f5
# ╠═ffb15055-c887-47a9-88d7-71a33a3c2dba
# ╠═77ed169b-d426-4f3f-8362-809cb0c8dd22
# ╠═79ddb157-8ee4-4c44-b8bf-b4f36f88acce
# ╠═6c589940-e00a-40c9-b2d9-c10384e4ac62
# ╠═5b5c5f78-75c2-4a85-b918-e7485f492d9e
# ╠═75368c67-8774-4f67-9163-d31475f4e289
# ╠═a14bd64d-4188-4682-8a47-c3999afb0d09
# ╠═06e6d38d-b08e-4583-863a-e6e06726143b
# ╠═41e118e2-2172-4910-ba52-12e199737c6a

max_stars_repo_name: teamclouday/JuliaLearn

max_issues_repo_name: teamclouday/JuliaLearn

max_forks_repo_name: teamclouday/JuliaLearn


text:
using DiffEqBase.InternalEuler, SciMLBase, DiffEqBase, Test

# Here's the problem to solve

struct LorenzFunction <: Function
    syms::Vector{Symbol}
end

function (::LorenzFunction)(u,p,t)
 [10.0(u[2]-u[1]),u[1]*(28.0-u[3]) - u[2],u[1]*u[2] - (8/3)*u[3]]
end
lorenz = LorenzFunction([:x,:y,:z])

u0 = [1., 5., 10.]
tspan = (0., 100.)
prob = ODEProblem(lorenz, u0, tspan)
dt = 0.1
sol = solve(prob,InternalEuler.FwdEulerAlg(), tstops=0:dt:1)
syms = [:x,:y,:z]

@test SciMLBase.has_syms(prob.f) == true
@test SciMLBase.interpret_vars([(0,1), (1,3), (4,5)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,1,3), (SciMLBase.DEFAULT_PLOT_FUNC,4,5)]
@test SciMLBase.interpret_vars([1, (1,3), (4,5)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,1,3), (SciMLBase.DEFAULT_PLOT_FUNC,4,5)]
@test SciMLBase.interpret_vars([1, 3, 4],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,0,3), (SciMLBase.DEFAULT_PLOT_FUNC,0,4)]
@test SciMLBase.interpret_vars(([1,2,3], [4,5,6]),sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,1,4), (SciMLBase.DEFAULT_PLOT_FUNC,2,5), (SciMLBase.DEFAULT_PLOT_FUNC,3,6)]
@test SciMLBase.interpret_vars((1, [2,3,4]),sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,1,2), (SciMLBase.DEFAULT_PLOT_FUNC,1,3), (SciMLBase.DEFAULT_PLOT_FUNC,1,4)]

@test SciMLBase.interpret_vars([(:t,:x),(:t,:y)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,0,2)]
@test SciMLBase.interpret_vars([:x, (0,:x), (:x,:y)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,1,2)]
@test SciMLBase.interpret_vars([:x, :y, :z],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,0,2), (SciMLBase.DEFAULT_PLOT_FUNC,0,3)]
@test SciMLBase.interpret_vars(([:x,:x], [:y,:z]),sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,1,2), (SciMLBase.DEFAULT_PLOT_FUNC,1,3)]
@test SciMLBase.interpret_vars((:x, [:y,:z]),sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,1,2), (SciMLBase.DEFAULT_PLOT_FUNC,1,3)]

f(x,y) = (x+y,y)
@test SciMLBase.interpret_vars([(f,0,1), (1,3), (4,5)],sol,syms) == [(f,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,1,3), (SciMLBase.DEFAULT_PLOT_FUNC,4,5)]
@test SciMLBase.interpret_vars([1, (f,1,3), (4,5)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (f,1,3), (SciMLBase.DEFAULT_PLOT_FUNC,4,5)]
@test SciMLBase.interpret_vars([(f,:t,:x),(:t,:y)],sol,syms) == [(f,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,0,2)]
@test SciMLBase.interpret_vars([:x, (f,0,:x), (:x,:y)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,0,1), (f,0,1), (SciMLBase.DEFAULT_PLOT_FUNC,1,2)]
@test SciMLBase.interpret_vars([(:x,:y)],sol,syms) == [(SciMLBase.DEFAULT_PLOT_FUNC,1,2)]
@test SciMLBase.interpret_vars((f,:x,:y),sol,syms) == [(f,1,2)]


text:
using TestPackage
using Test

@testset "TestPackage.jl" begin
    # Write your tests here.
    @test addition(1, 2) == 3
    @test eigvals(ones(2, 2)) == [0., 2.]
end


text:
# Dispatchable and non-dispatchable generators

## Expressions

"Curtailed power of a non-dispatchable generator as the difference between its reference power and the generated power."
function expression_gen_curtailment(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, report::Bool=true)
    pgcurt = _PM.var(pm, nw)[:pgcurt] = Dict{Int,Any}(
        i => ndgen["pmax"] - _PM.var(pm,nw,:pg,i) for (i,ndgen) in _PM.ref(pm,nw,:ndgen)
    )
    if report
        _IM.sol_component_fixed(pm, _PM.pm_it_sym, nw, :gen, :pgcurt, _PM.ids(pm, nw, :dgen), 0.0)
        _PM.sol_component_value(pm, nw, :gen, :pgcurt, _PM.ids(pm, nw, :ndgen), pgcurt)
    end
end


text:
function addone!(a::Vector{Int})
    for i = reverse(1:length(a))
        if a[i] == 0
            a[i] = 1
            return
        else
            a[i] = 0
        end
    end
    # error("Overflow adding 1 to $a".)
end


text:
using SafeTestsets

@safetestset "Bandit Tests" begin include("bandittests.jl") end
@safetestset "MDP Tests" begin include("mdptest.jl") end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 1b737805-a411-4585-b215-d0f99eafac0c
using HypertextLiteral

# ╔═╡ 881b75d3-cebe-4a53-8cf5-beaeeddafd35
function skip_as_script(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ dadf2f40-1764-47a4-b560-683b6479d77f
if skip_as_script(@__MODULE__)
	import Pkg
	Pkg.activate(Base.current_project(@__DIR__))
	Text("Project env active")
end

# ╔═╡ 38b7eeb9-80bb-4a3a-a2d2-809fc423625c
if skip_as_script(@__MODULE__)
	using PlutoUI
end

# ╔═╡ ba94c9b3-55bf-486a-8bf3-eb0b3e65e536
"""
	@skip_as_script expression

Marks a expression as Pluto-only, which means that it won't be executed when running outside Pluto. Do not use this for your own projects.
"""
macro skip_as_script(ex) skip_as_script(__module__) ? esc(ex) : nothing end

# ╔═╡ 9e56dbb9-2e93-4136-b10b-7b1ecb07e6f1
macro only_as_script(ex) skip_as_script(__module__) ? nothing : esc(ex) end

# ╔═╡ 5028fc01-4e14-4550-8ae6-48703b86dc21
import AbstractPlutoDingetjes.Bonds

# ╔═╡ aa3b170e-86fd-42e2-a82e-99e25436c65e
import AbstractPlutoDingetjes

# ╔═╡ ccbe4423-801b-4c20-ad2d-ee89d5cfa859
md"""
# Combining bonds
"""

# ╔═╡ 85918609-5e1f-4040-99be-61c2dd8ff654
md"""
## The magic
"""

# ╔═╡ 957858f8-4ace-4cae-bb16-49569daa9869
const set_input_value_compat = HypertextLiteral.JavaScript("""
(() => {
	let result = null
	try {
	result = setBoundElementValueLikePluto
} catch (e) {
	result = ((input, new_value) => {
	// fallback in case https://github.com/fonsp/Pluto.jl/pull/1755 is not available
    if (new_value == null) {
        //@ts-ignore
        input.value = new_value
        return
    }
    if (input instanceof HTMLInputElement) {
        switch (input.type) {
            case "range":
            case "number": {
                if (input.valueAsNumber !== new_value) {
                    input.valueAsNumber = new_value
                }
                return
            }
            case "date": {
                if (input.valueAsDate == null || Number(input.valueAsDate) !== Number(new_value)) {
                    input.valueAsDate = new_value
                }
                return
            }
            case "checkbox": {
                if (input.checked !== new_value) {
                    input.checked = new_value
                }
                return
            }
            case "file": {
                // Can't set files :(
                return
            }
        }
    } else if (input instanceof HTMLSelectElement && input.multiple) {
        for (let option of Array.from(input.options)) {
            option.selected = new_value.includes(option.value)
        }
        return
    }
    //@ts-ignore
    if (input.value !== new_value) {
        //@ts-ignore
        input.value = new_value
    }
})
}
return result
})()
""")

# ╔═╡ 3f6a3cbc-632d-413d-990e-0f06730bd27c
begin
	local output = begin
	"""
	```julia
	RenderCallback(callback::Function, x::Any)
	```
	An HTML display passthrough of `x` (displays the same content), but when it is displayed, a callback function is invoked. `disable_callback!` can remove a callback.
	"""
	struct RenderCallback
		callback_ref::Ref{Union{Function,Nothing}}
		content::Any
	end
	end

	function disable_callback!(rc::RenderCallback)
		rc.callback_ref[] = nothing
	end
	function Base.show(io::IO, m::MIME"text/html", rc::RenderCallback)
		if rc.callback_ref[] !== nothing
			rc.callback_ref[]()
		end
		Base.show(io, m, rc.content)
	end
	output
end

# ╔═╡ 19613f3f-5825-45a4-8951-8ff1043d0867
begin
	Base.@kwdef struct CombinedBonds
		display_content::Any
		captured_names::Union{Nothing,NTuple{N,Symbol} where N}
		captured_bonds::Vector{Any}
		secret_key::String
	end

	function Base.show(io::IO, m::MIME"text/html", cb::CombinedBonds)
		if !AbstractPlutoDingetjes.is_supported_by_display(io, Bonds.transform_value)
			return Base.show(io, m, HTML("❌ You need to update Pluto to use this PlutoUI element."))
		end
		output = @htl(
			"""$(
				cb.display_content
			)""")
		Base.show(io, m, output)
	end

	function Bonds.initial_value(cb::CombinedBonds)
		vals = tuple((Bonds.initial_value(b) for b in cb.captured_bonds)...)

		if cb.captured_names === nothing
			vals
		else
			NamedTuple{cb.captured_names}(vals)
		end
	end
	function Bonds.validate_value(cb::CombinedBonds, from_js)
		if from_js isa Vector && length(from_js) == length(cb.captured_bonds)
			all((
				Bonds.validate_value(bond, val_js)
				for (bond, val_js) in zip(cb.captured_bonds, from_js)
			))
		else
			false
		end
	end
	function Bonds.transform_value(cb::CombinedBonds, from_js)
		@assert from_js isa Vector
		@assert length(from_js) == length(cb.captured_bonds)

		vals = tuple((
			Bonds.transform_value(bond, val_js)
			for (bond, val_js) in zip(cb.captured_bonds, from_js)
		)...)

		if cb.captured_names === nothing
			vals
		else
			NamedTuple{cb.captured_names}(vals)
		end

		# [
		# 	Bonds.transform_value(bond, val_js)
		# 	for (bond, val_js) in zip(cb.captured_bonds, from_js)
		# ]
	end

	# TODO:
	# function Bonds.possible_values (cb::CombinedBonds, from_js)
	# end

end

# ╔═╡ c8e01722-3372-4692-bc8a-037846c950c2
"""
Same as `repr(MIME"text/html"(), x)` but the `IOContext` is set up to match the one used by Pluto to render. This means that if the object being rendered wants to use `AbstractPlutoDingetjes.is_supported_by_display`, they can.
"""
render_hmtl_with_pluto_display_features(x) = sprint() do io
	Base.show(
		IOContext(io,
			:pluto_supported_integration_features => [
				AbstractPlutoDingetjes,
				AbstractPlutoDingetjes.Bonds,
				AbstractPlutoDingetjes.Bonds.initial_value,
				AbstractPlutoDingetjes.Bonds.transform_value,
				AbstractPlutoDingetjes.Bonds.possible_values,
			],
			:color => false, :limit => true,
			:displaysize => (18, 88), :is_pluto => true,
		),
		MIME"text/html"(),
		x
	)
end

# ╔═╡ bafa80c7-bb2b-4c08-a9ae-5e2b7d9abe07
function _combine(f::Function)
	key = String(rand('a':'z', 10))

	captured_names = Symbol[]
	captured_bonds = []

	function combined_child_element(x)
		@htl("""<pl-combined_child key=$(key)>$(x)</pl-combined_child>""")
	end

	created_callbacks = RenderCallback[]

	#=

	This is the function that will wrap around contained input elements.

	Besides wrapping the input inside a special HTML element
	(with `combined_child_element`), it also secretly adds the element to
	`captured_bonds`.

	This allows us to know exactly which elements are contained in the combine,
	which we use for `Bonds.initial_value`, `Bonds.transform_value`, etc.

	This function could be a lot simpler:

	```
	function Child(x)
		push!(captured_bonds, x)
		combined_child_element(x)
	end
	```

	But instead, it's complicated!

	The reason is that we want to capture the combined bonds not in the order
	that they are wrapped (i.e. when `Child` is called), but in the order that they
	are **displayed** in, because this matches the order that our JavaScript code
	will find the special elements in.

	=#
	function Child(x)
		rc = RenderCallback(combined_child_element(x)) do
			# @info "Rendering" x stacktrace()
			push!(captured_bonds, x)
		end
		push!(created_callbacks, rc)
		rc
	end
	# the same, but with `name`
	function Child(name::Union{String,Symbol}, x)
		rc = RenderCallback(combined_child_element(x)) do
			push!(captured_bonds, x)
			push!(captured_names, Symbol(name))
		end
		push!(created_callbacks, rc)
		rc
	end

	# call the user's render function
	result = f(Child)

	# Trigger HTML display, which will also render the Child elements, and fire our callbacks.
	# We store the result because we don't want to re-render the contents every time the combine is rendered: if the displayed content contains lazy generators, then blablablalbl difficult but fixed now -fons
	display_html = render_hmtl_with_pluto_display_features(result)

	# @info "Captured" captured_bonds length(created_callbacks)

	# disable callbacks
	disable_callback!.(created_callbacks)

	# lets see what we got
	@assert isempty(captured_names) || length(captured_names) == length(captured_bonds) "Some children do not have a name. Make sure that all calls of `Child` provide two arguments."

	CombinedBonds(;
		secret_key = key,
		captured_names = isempty(captured_names) ? nothing : tuple(captured_names...),
		captured_bonds = captured_bonds,
		display_content = HTML(display_html),
	)
end

# ╔═╡ 157a2e04-4ccd-4de6-b998-ef94fcdd962b
"""
```julia
PlutoUI.combine(render_function::Function)
```

Combine multiple input elements into one. The combined values are sent to `@bind` as a single tuple.

`render_function` is a function that you write yourself, take a look at the examples below.

# Examples

## 🐶 & 🐱

We use the [`do` syntax](https://docs.julialang.org/en/v1/manual/functions/#Do-Block-Syntax-for-Function-Arguments) to write our `render_function`. The `Child` function is wrapped around each input, to let `combine` know which values to combine.

```julia
@bind values PlutoUI.combine() do Child
	md""\"
	# Hi there!

	I have \$(
		Child(Slider(1:10))
	) dogs and \$(
		Child(Slider(5:100))
	) cats.

	Would you like to see them? \$(Child(CheckBox(true)))
	""\"
end

values == (1, 5, true) # (initially)
```

> The output looks like:
>
> ![screenshot of running the code above inside Pluto](https://user-images.githubusercontent.com/6933510/145589918-25a3c732-c02e-482b-831b-06131b283597.png)

## 🎏

The `combine` function is most useful when you want to generate your input elements _dynamically_. This example uses [HypertextLiteral.jl](https://github.com/MechanicalRabbit/HypertextLiteral.jl) for the `@htl` macro:

```julia
function wind_speeds(directions)
	PlutoUI.combine() do Child
		@htl(""\"
		<h6>Wind speeds</h6>
		<ul>
		\$([
			@htl("<li>\$(name): \$(Child(Slider(1:100)))</li>")
			for name in directions
		])
		</ul>
		""\")
	end
end

@bind speeds wind_speeds(["North", "East", "South", "West"])

speeds == (1, 1, 1, 1) # (initially)

# after moving the sliders:
speeds == (100, 36, 73, 60)
```

> The output looks like:
>
> ![screenshot of running the code above inside Pluto](https://user-images.githubusercontent.com/6933510/145588612-14824654-5c73-45f8-983c-8913c7101a78.png)

# Named variant

In the last example, we used `Child` to wrap around contained inputs:
```julia
Child(Slider(1:100))
```
We can also use the **named variant**, which looks like:
```julia
Child("East", Slider(1:100))
```

When you use the *named variant* for all children, **the bound value will be `NamedTuple`, instead of a `Tuple`**.

Let's rewrite our example to use the *named variant*:

```julia
function wind_speeds(directions)
	PlutoUI.combine() do Child
		@htl(""\"
		<h6>Wind speeds</h6>
		<ul>
		\$([
			@htl("<li>\$(name): \$(Child(name, Slider(1:100)))</li>")
			for name in directions
		])
		</ul>
		""\")
	end
end

@bind speeds wind_speeds(["North", "East", "South", "West"])

speeds == (North=1, East=1, South=1, West=1) # (initially)

# after moving the sliders:
speeds == (North=100, East=36, South=73, West=60)

md"The Eastern wind speed is \$(speeds.East)."
```

> The output looks like:
>
> ![screenshot of running the code above inside Pluto](https://user-images.githubusercontent.com/6933510/145615489-b3fb910d-0dc1-408b-882f-b05ca0129b18.gif)

# Why?
## You can make a new widget!
You can use `combine` to **create a brand new widget** yourself, by combining existing ones!

In the example above, we created a new widget called `wind_speeds`. You could **put this function in a package** (e.g. `WeatherUI.jl`) and people could use it like so:

```julia
import WeatherUI: wind_speeds

@bind speeds wind_speeds(["North", "East"])
```

In other words: you can use `combine` to create application-specific UI elements, and you can put those in your package!

## Difference with repeated `@bind`
The standard way to combine multiple inputs into one output is to use `@bind` multiple times. Our initial example could more easily be written as:
```julia
md""\"
# Hi there!

I have \$(@bind num_dogs Slider(1:10)) dogs and \$(@bind num_cats Slider(5:10)) cats.

Would you like to see them? \$(@bind want_to_see CheckBox(true))
""\"
```

The `combine` function is useful when you are generating inputs **dynamically**, like in our second example. This is useful when:
- The number of parameters is very large, and you don't want to write `@bind parameter1 ...`, `@bind parameter2 ...`, etc.
- The number of parameters is dynamic! For example, you can load in a table in one cell, and then use `combine` in another cell to select which rows you want to use.

"""
combine(f::Function) = _combine(f)

# ╔═╡ ad5cffa5-313c-4de9-9360-005365b40780
export combine

# ╔═╡ 5f8af778-d25b-462f-afb9-295c0e4cb12e
RenderCallback(@htl("hello")) do
	nothing
end

# ╔═╡ 0d6800f3-9cdc-4531-982d-d7607748d8f5
let
	values = []
	rc = RenderCallback(@htl("asdf")) do
		push!(values, 123)
	end
	repr(MIME"text/html"(), rc)
	disable_callback!(rc)
	repr(MIME"text/html"(), rc)
	repr(MIME"text/html"(), rc)
	values
end

# ╔═╡ 6c8a03e4-7d8e-4aa4-a750-7b815622147d
md"""
## Examples
"""

# ╔═╡ f08680b2-ed28-4fac-838c-2eca7e75c6dc
@skip_as_script @bind values combine() do Child
	md"""
	# Hi there!

	I have $(Child(Slider(1:10))) dogs and $(Child(Slider(5:10))) cats.

	Would you like to see them? $(Child(CheckBox(true)))
	"""
end

# ╔═╡ 8dcb2498-00cf-49a8-8074-301fe88b76ea
@skip_as_script values

# ╔═╡ 801fb021-73a0-4114-a36a-328e84f00b51
@skip_as_script @bind speeds combine() do Child
	@htl("""

Wind speeds

	$([
		@htl("
$(name): $(Child(Slider(1:100)))")
		for name in ["North", "East", "South", "West"]
	])

	""")
end

# ╔═╡ d7985844-5944-42b9-ad41-599cd72eea82
@skip_as_script speeds

# ╔═╡ 39be652f-bcd7-4a4f-8b03-5167285dc7a2
@skip_as_script @bind speeds_named combine() do Child
	@htl("""

Wind speeds

	$([
		@htl("
$(name): $(Child(name, Slider(1:100)))")
		for name in ["North", "East", "South", "West"]
	])

	""")
end

# ╔═╡ c8e09bf5-eae4-492f-946e-4bb55ff7e161
@skip_as_script speeds_named

# ╔═╡ c3176549-7c46-4d89-9c76-817c6bfcbb71
@skip_as_script rb = @bind together combine() do Child
	@htl("""

Hello world!

	$(Child(TextField()))
	$(Child(Slider([sin, cos, tan])))
	$([
		Child(Scrubbable(7))
		for _ in 1:3
	])

	""")
end

# ╔═╡ 5300b80d-4097-4b85-ad2a-3a54b01559bf
@skip_as_script let
	sleep(.5)
	together
end

# ╔═╡ 6d4f62a6-e3d9-416f-a1e3-694c31cfbb31
@skip_as_script rb

# ╔═╡ 5a0196d0-e19f-4202-b36d-18ab9be839b3
md"""
## Tests
"""

# ╔═╡ 3d3a6abb-bea7-41e2-862d-9536a9687ea5
md"""
### Initial value & transform
"""

# ╔═╡ 3cd1f92a-a088-4695-8eea-ada08e01a7c2
@skip_as_script begin
	it2vs = []
	it2b = @bind it2v combine() do Child

		z = Child(Slider(1:10))

		@htl("""

Hello world!

		$(z)
		$(Child(Slider([sin, cos, tan])))
		$(z)

		""")
	end
end

# ╔═╡ 62368acd-d035-4795-804b-0f024d7333a7
@skip_as_script push!(it2vs, it2v)

# ╔═╡ 0b643ed9-914a-475b-a1ce-d840df1fc223
@skip_as_script begin
	itvs = []
	itb = @bind itv combine() do Child
		@htl("""

Hello world!

		$(Child(@htl("<input type=range>")))
		$(Child(Slider([sin, cos, tan])))

		""")
	end
end

# ╔═╡ 59a199f4-4ccc-4319-9d32-03da3adbb5db
# @skip_as_script itb

# ╔═╡ cc3e475e-d2c2-4b40-a1ba-033576aefdae
md"""
Should be `[[missing, sin], [50, sin]]`
"""

# ╔═╡ 127dd58b-9ba0-41d8-91e6-357d9ec63e6b
@skip_as_script push!(itvs, itv)

# ╔═╡ 55133b9c-ea1e-4f19-bb52-d541bc461fcd
@skip_as_script begin
	☎s = []
	☎c = combine() do Child
		@htl("""
		$((
			Child(Slider(LinRange(rand(), rand()+1, 10))) for x in 1:5
		))
		""")
	end
	☎b = @bind ☎ ☎c
end

# ╔═╡ 1fa7b945-3102-4a05-b7a1-5610c65aac0f
@skip_as_script ☎c.captured_bonds

# ╔═╡ 0a87a4d3-299e-4ec6-82be-5b43a1e29077
@skip_as_script push!(☎s, ☎)

# ╔═╡ 7de50435-e962-4405-b06d-83fbe9436cff
md"""
### Combine inside combine
"""

# ╔═╡ 3396a5ba-533e-4e0d-ab4f-c633459bd81a
@skip_as_script cb1 = combine() do Child
	md"""
	Left: $(Child(:left, Slider(1:10))), right: $(Child(:right, Scrubbable(5)))
	"""
end

# ╔═╡ 829f90ac-c4d9-4af9-984a-b4e52a075460
@skip_as_script cb2b = @bind wowz combine() do Child
	md"""
	Do a thing: $(Child(CheckBox()))

	#### Up
	$(Child(cb1))

	#### Down
	$(Child(cb1))

	"""
end

# ╔═╡ a1a9b22f-9df7-4c5f-a73d-38e184da7c35
@skip_as_script wowz

# ╔═╡ c7a32ab2-2a5b-485b-9795-4a743604de82
@skip_as_script cb2b

# ╔═╡ Cell order:
# ╟─881b75d3-cebe-4a53-8cf5-beaeeddafd35
# ╟─ba94c9b3-55bf-486a-8bf3-eb0b3e65e536
# ╟─9e56dbb9-2e93-4136-b10b-7b1ecb07e6f1
# ╟─dadf2f40-1764-47a4-b560-683b6479d77f
# ╠═1b737805-a411-4585-b215-d0f99eafac0c
# ╠═5028fc01-4e14-4550-8ae6-48703b86dc21
# ╠═aa3b170e-86fd-42e2-a82e-99e25436c65e
# ╟─ccbe4423-801b-4c20-ad2d-ee89d5cfa859
# ╠═ad5cffa5-313c-4de9-9360-005365b40780
# ╠═157a2e04-4ccd-4de6-b998-ef94fcdd962b
# ╠═bafa80c7-bb2b-4c08-a9ae-5e2b7d9abe07
# ╟─85918609-5e1f-4040-99be-61c2dd8ff654
# ╟─c8e01722-3372-4692-bc8a-037846c950c2
# ╠═19613f3f-5825-45a4-8951-8ff1043d0867
# ╟─957858f8-4ace-4cae-bb16-49569daa9869
# ╠═3f6a3cbc-632d-413d-990e-0f06730bd27c
# ╠═5f8af778-d25b-462f-afb9-295c0e4cb12e
# ╠═0d6800f3-9cdc-4531-982d-d7607748d8f5
# ╟─6c8a03e4-7d8e-4aa4-a750-7b815622147d
# ╠═38b7eeb9-80bb-4a3a-a2d2-809fc423625c
# ╠═f08680b2-ed28-4fac-838c-2eca7e75c6dc
# ╠═8dcb2498-00cf-49a8-8074-301fe88b76ea
# ╠═801fb021-73a0-4114-a36a-328e84f00b51
# ╠═d7985844-5944-42b9-ad41-599cd72eea82
# ╠═39be652f-bcd7-4a4f-8b03-5167285dc7a2
# ╠═c8e09bf5-eae4-492f-946e-4bb55ff7e161
# ╠═c3176549-7c46-4d89-9c76-817c6bfcbb71
# ╠═5300b80d-4097-4b85-ad2a-3a54b01559bf
# ╠═6d4f62a6-e3d9-416f-a1e3-694c31cfbb31
# ╟─5a0196d0-e19f-4202-b36d-18ab9be839b3
# ╟─3d3a6abb-bea7-41e2-862d-9536a9687ea5
# ╠═3cd1f92a-a088-4695-8eea-ada08e01a7c2
# ╠═62368acd-d035-4795-804b-0f024d7333a7
# ╠═0b643ed9-914a-475b-a1ce-d840df1fc223
# ╠═59a199f4-4ccc-4319-9d32-03da3adbb5db
# ╟─cc3e475e-d2c2-4b40-a1ba-033576aefdae
# ╠═127dd58b-9ba0-41d8-91e6-357d9ec63e6b
# ╠═55133b9c-ea1e-4f19-bb52-d541bc461fcd
# ╠═1fa7b945-3102-4a05-b7a1-5610c65aac0f
# ╠═0a87a4d3-299e-4ec6-82be-5b43a1e29077
# ╟─7de50435-e962-4405-b06d-83fbe9436cff
# ╠═3396a5ba-533e-4e0d-ab4f-c633459bd81a
# ╠═829f90ac-c4d9-4af9-984a-b4e52a075460
# ╠═a1a9b22f-9df7-4c5f-a73d-38e184da7c35
# ╠═c7a32ab2-2a5b-485b-9795-4a743604de82


text:
include("base.jl")

function analyze(widget::QStack, ctx::PlotContext)
    local dates = Array{DateTime,1}()
    local srcCtx = Dict{String,DataAttributeContext}()
    local mini = typemax(Float64)
    local maxi = typemin(Float64)
    local titles = Array{String,1}()

    for attr in widget.stacked
        local _dates, _srcCtx, _mini, _maxi, _titles = analyze(attr, ctx)

        if length(dates) > 0
            @assert length(dates) == length(_dates) "$(length(dates)) == $(length(_dates))"
        end
        dates = _dates

        if widget.ctxs !== nothing
            if widget.ctxs |> length === 0
                _srcCtx = Dict()
            end
            for c in widget.ctxs
                if "$c" in _srcCtx
                    mini = min(mini, _mini)
                    maxi = max(maxi, _maxi)
                    delete!(_srcCtx, "$c")
                end
            end
        end

        for (k, v) in _srcCtx
            srcCtx[k] = v
        end
        mini = min(mini, _mini)
        maxi = Inf # We can't calculate the maximum. Inf is "automatic" in plot.jl ylims

        titles = [titles..., _titles...]
    end

    return dates, srcCtx, mini, maxi, titles
end

function renderWidget!(widget::QStack, ctx::PlotContext)

    local myStyle = ctx.styles
    local i = 1
    local accum = 0.0

    # TODO: What about stacking lines with different resolution?
    # One would have to interpolate values and sort the dates...

    for attr in widget.stacked
        ctx.styles = deepcopy(myStyle)
        delete!(ctx.styles, "source context")
        local _, srcCtx, _, _, _ = analyze(QStack(widget.stacked), ctx)

        local ctxs = []
        if widget.ctxs === nothing
            ctxs = [nothing]
        else
            ctxs = collect(values(srcCtx))
            if length(widget.ctxs) > 0
                ctxs = filter(y -> "$y" in map(x -> "$x", widget.ctxs), ctxs)
            end
        end
        # TODO: Does only support single source contexts!

        for c in ctxs
            try
                ctx.styles = deepcopy(myStyle)

                if c !== nothing
                    logger(c, "Plotting where stack context is", true)
                    ctx.styles["source context"] = c
                    ctx.styles["append to label"] = ""
                    for cc in c.ctx
                        ctx.styles["append to label"] *= "$(cc[2])=$(cc[3])"
                    end
                end

                ctx.styles["fillrange"] = accum
                ctx.styles["offset"] = accum
                ctx.styles["stacked"] = true
                local dates, numbers = renderWidget!(attr, ctx)
                accum = numbers
            catch exception
                @show "exception in stacks"
                bt = catch_backtrace()
                msg = sprint(showerror, exception, bt)
                logger(msg, "failed", true, failed = true)
                # TODO: only catch context-not-found exceptions
            end
        end
    end

    return dates, accum
end

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
# Test of reduction on scalar backend

# Unit testing for element-wise expressions on scalar backend

import Devectorize
import Devectorize.@devec
import Devectorize.@inspect_devec
import Devectorize.dump_devec
import Devectorize.sqr

using Base.Test

# data

a = [3., 4., 5., 6., 8., 7., 6., 5.]
b = [9., 8., 7., 6., 4., 2., 3., 1.]
c = [1., 2., 4., 3., 5., 7., 6., 8.]
abc = [a b c]

#################################################
#
#   full reduction
#
#################################################

r = zeros(1)
@devec r = sum(a)
@test isequal(r, sum(a))

@devec r = sum(a[:,:])
@test isequal(r, sum(a))

@devec r = sum(abc[:,:])
@test isequal(r, sum(abc))

@devec r = max(a)
@test isequal(r, max(a))

@devec r = max(c)
@test isequal(r, max(c))

@devec r = min(a)
@test isequal(r, min(a))

@devec r = min(c)
@test isequal(r, min(c))

@devec r = mean(a)
@test isequal(r, mean(a))

@devec r = mean(abc[:,:])
@test isequal(r, mean(abc))

@devec r = dot(a, b)
@test isequal(r, dot(a, b))

@devec r = dot(a[:,:], b[:,:])
@test isequal(r, dot(a, b))

@devec r = dot(abc[:,:], abc)
@test isequal(r, dot(abc[:], abc[:]))

#################################################
#
#   partial reduction
#
#################################################

@devec r = sum(abc, 1)
@test isequal(r, sum(abc, 1))

@devec r = sum(abc, 2)
@test isequal(r, sum(abc, 2))

r = zeros(size(abc, 2))
r0 = r
@devec r[:] = sum(abc, 1)
@test r === r0
@test isequal(r, vec(sum(abc, 1)))

r = zeros(size(abc, 1))
r0 = r
@devec r[:] = sum(abc, 2)
@test r === r0
@test isequal(r, vec(sum(abc, 2)))

@devec r = mean(abc, 1)
@test isequal(r, sum(abc, 1) / size(abc, 1))

@devec r = mean(abc, 2)
@test isequal(r, sum(abc, 2) / size(abc, 2))

@devec r = max(abc, (), 1)
@test isequal(r, max(abc, (), 1))

@devec r = max(abc, (), 2)
@test isequal(r, max(abc, (), 2))

@devec r = min(abc, (), 1)
@test isequal(r, min(abc, (), 1))

@devec r = min(abc, (), 2)
@test isequal(r, min(abc, (), 2))

@devec r = sum(sqr(abc), 1)
@test isequal(r, sum(abc .* abc, 1))

@devec r = sum(sqr(abc), 2)
@test isequal(r, sum(abc .* abc, 2))


text:
const DimSetters = Union{LookupArraySetters,Type,UnionAll,Dimension,Symbol}

set(dim::Dimension, x::DimSetters) = _set(dim, x)
# Convert args/kw to dims and set
_set(dims_::DimTuple, args::Dimension...; kw...) = _set(dims_, (args..., kwdims(kw)...))
# Convert pairs to wrapped dims and set
_set(dims_::DimTuple, p::Pair, ps::Vararg{<:Pair}) = _set(dims_, (p, ps...))
_set(dims_::DimTuple, ps::Tuple{Vararg{<:Pair}}) = _set(dims_, pairdims(ps...))
# Set dims with (possibly unsorted) wrapper vals
_set(dims::DimTuple, wrappers::DimTuple) = begin
    # Check the dimension types match
    map(wrappers) do w
        hasdim(dims, w) || _wrongdimserr(dims, w)
    end
    # Missing dims return `nothing` from sortdims
    newdims = map(_set, dims, sortdims(wrappers, dims))
    # Swaps existing dims with non-nothing new dims
    swapdims(dims, newdims)
end

# Set things wrapped in dims
_set(dim::Dimension, wrapper::Dimension{<:DimSetters}) = _set(dim::Dimension, val(wrapper))
# Set the dim, checking the lookup
_set(dim::Dimension, newdim::Dimension) = _set(newdim, _set(val(dim), val(newdim)))
# Construct types
_set(dim::Dimension, ::Type{T}) where T = _set(dim, T())
_set(dim::Dimension, key::Symbol) = _set(dim, key2dim(key))
_set(dim::Dimension, dt::DimType) = basetypeof(dt)(val(dim))
_set(dim::Dimension, x) = rebuild(dim; val=_set(val(dim), x))
# Set the lookup
# Otherwise pass this on to set fields on the lookup
_set(dim::Dimension, x::LookupArrayTrait) = rebuild(dim, _set(lookup(dim), x))

# Metadata
_set(dim::Dimension, newmetadata::AllMetadata) = rebuild(dim, _set(lookup(dim), newmetadata))

_set(x::Dimension, ::Nothing) = x
_set(::Nothing, x::Dimension) = x
_set(::Nothing, ::Nothing) = nothing
_set(x, ::Nothing) = x
_set(::Nothing, x) = x

@noinline _wrongdimserr(dims, w) = throw(ArgumentError("dim $(basetypeof(w))) not in $(map(basetypeof, dims))"))


text:
using RadiativeTransfer.vSmartMOM
using RadiativeTransfer
using RadiativeTransfer.Scattering
using CUDA
using Test
# using TensorOperations
using KernelAbstractions
using KernelAbstractions.Extras
using StaticArrays
using BenchmarkTools
using TensorOperations
using LinearAlgebra
using NNlib

if has_cuda_gpu()
    CUDA.allowscalar(false)
end

# Needs some warning if memory is getting too large !
FT = Float64
n = 16
nSpec = 100

# Create CPU Matrices:
r⁻⁺_ = (randn(FT, n, n, nSpec));
t⁺⁺_ = (randn(FT, n, n, nSpec));
r⁺⁻_ = (randn(FT, n, n, nSpec));
t⁻⁻_ = (randn(FT, n, n, nSpec));
R⁻⁺_ = (randn(FT, n, n, nSpec));
T⁺⁺_ = (randn(FT, n, n, nSpec));
R⁺⁻_ = (randn(FT, n, n, nSpec));
T⁻⁻_ = (randn(FT, n, n, nSpec));

# And move to GPU as CuArray
r⁻⁺ = CuArray(r⁻⁺_);
t⁺⁺ = CuArray(t⁺⁺_);
r⁺⁻ = CuArray(r⁺⁻_);
t⁻⁻ = CuArray(t⁻⁻_);
R⁻⁺ = CuArray(R⁻⁺_);
T⁺⁺ = CuArray(T⁺⁺_);
R⁺⁻ = CuArray(R⁺⁻_);
T⁻⁻ = CuArray(T⁻⁻_);

# This is still a bit stupid but noy sure how to make this better...
# I_static_ = Diagonal{FT}(ones(n))
# I_static = CuArray(repeat(I_static, 1, 1))
I_static_ = Diagonal{FT}(ones(n))
I_static  = Diagonal(CuArray(I_static_));

function rt_interaction!(R⁻⁺::AbstractArray{FT,3},
                         T⁺⁺::AbstractArray{FT,3},
                         R⁺⁻::AbstractArray{FT,3},
                         T⁻⁻::AbstractArray{FT,3},
                         r⁻⁺::AbstractArray{FT,3},
                         t⁺⁺::AbstractArray{FT,3},
                         r⁺⁻::AbstractArray{FT,3},
                         t⁻⁻::AbstractArray{FT,3},
                         I_static::AbstractArray) where {FT}
    # Used to store `(I - R⁺⁻ * r⁻⁺)⁻¹ * T⁺⁺`
    tmp_inv = similar(t⁺⁺)

    # Compute and store `(I - R⁺⁻ * r⁻⁺)⁻¹ * T⁺⁺`
    vSmartMOM.batch_solve!(tmp_inv, I_static .- R⁺⁻ ⊠ r⁻⁺, T⁺⁺)
    # synchronize()
    R⁻⁺ += (T⁻⁻ ⊠ r⁻⁺ ⊠ tmp_inv)
    # synchronize()
    T⁺⁺[:] = t⁺⁺ ⊠ tmp_inv

    # Repeating for mirror-reflected directions

    # Compute and store `(I - r⁻⁺ * R⁺⁻)⁻¹ * t⁻⁻`
    # synchronize()
    vSmartMOM.batch_solve!(tmp_inv, I_static .- r⁻⁺ ⊠ R⁺⁻, t⁻⁻)
    synchronize()
    R⁺⁻ .= r⁺⁻ + t⁺⁺ ⊠ R⁺⁻ ⊠ tmp_inv
    T⁻⁻[:] = tmp_inv
    # synchronize()
    T⁻⁻[:] = T⁺⁺ ⊠ tmp_inv

    return nothing
end

@show n, nSpec
function run_rt_interaction(R⁻⁺, T⁺⁺, R⁺⁻, T⁻⁻, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻,  I_static)
    rt_interaction!(R⁻⁺, T⁺⁺, R⁺⁻, T⁻⁻, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻,  I_static)
    synchronize()
end

# Testing RT Interaction time:
println("RT Interaction GPU time:")
@time run_rt_interaction(R⁻⁺, T⁺⁺, R⁺⁻, T⁻⁻, r⁻⁺, t⁺⁺, r⁺⁻, t⁻⁻,  I_static)
println("RT Interaction CPU time (uses multi-threading in LAPACK/BLAS!):")
@time run_rt_interaction(R⁻⁺_, T⁺⁺_, R⁺⁻_, T⁻⁻_, r⁻⁺_, t⁺⁺_, r⁺⁻_, t⁻⁻_, I_static_)

@testset "GPU-CPU consistency" begin
    for (matGPU, matCPU) in ((R⁻⁺, R⁻⁺_),
                             (T⁺⁺, T⁺⁺_),
                             (R⁺⁻, R⁺⁻_),
                             (T⁻⁻, T⁻⁻_),
                             (r⁻⁺, r⁻⁺_),
                             (t⁺⁺, t⁺⁺_),
                             (r⁺⁻, r⁺⁻_),
                             (t⁻⁻, t⁻⁻_))
        @show Array(matGPU) ≈ matCPU
        @test Array(matGPU) ≈ matCPU
    end
end


text:
using Flux, Optim, OrdinaryDiffEq, RecursiveArrayTools, DifferentialEquations
include("./DiffEqFlux.jl")

function main()
  n = 5  # number of ODEs
  tspan = (0.0, 1.0)

  x_data = rand(n, 5)
  y_data = rand(n, 5)

  training_data = Flux.Data.DataLoader(x_data, y_data, batchsize=1)

  NN = Chain(Dense(n, 10n, tanh),
             Dense(10n, n))

  function cb(p, l, x, pred)
    println("This is the size of the output of nODE ", size(pred))
    false
  end

  nODE = DiffEqFlux.NeuralODEBatchConv(NN, tspan, ROCK4(),
                              reltol=1e-4, saveat=[tspan[end]])

  function loss_function(θ, x, y)
    pred = nODE(x, θ)
    Flux.mse(y, pred), x, pred
  end
  res = DiffEqFlux.sciml_train(loss_function, nODE.p, ADAM(0.01), training_data,
                               cb=cb, save_best=false, progress=true)
end
@timev main()


text:
# Harmonic "water" model.

"""
    HarmonicWater

Purely harmonic "water" model.
"""
struct HarmonicWater <: WaterModel
    "Potential terms (returning kJ / mol)."
    potential::TermContainer
    "Force terms (returning kJ / nm mol)."
    force::TermContainer
end

function HarmonicWater(N::Int)
    build_model(HarmonicWater,
                harmonic_potential_intra, harmonic_potential_inter,
                harmonic_force_intra, harmonic_force_inter,
                N)
end

const harmonic_k_intra = 1.0e5 # kJ / nm^2 mol
const harmonic_k_inter = 1.0e2 # kJ / nm^2 mol

function harmonic_potential_intra(n::Int)
    function term_f(j::Int, s::State{P,D,A,N}) where {P,D,A,N}
        result = 0.0 # kJ / mol

        for a1 in 1:A
            for a2 in (a1+1):A
                R = norm(s.qs[j, :, a1, n] .- s.qs[j, :, a2, n]) # nm
                result += harmonic_k_intra * R^2
            end
        end

        0.5 * result
    end
end

function harmonic_potential_inter(n1::Int, n2::Int)
    function term_f(j::Int, s::State{P,D,A,N}) where {P,D,A,N}
        result = 0.0 # kJ / mol

        for a1 in 1:A
            for a2 in 1:A
                R = norm(s.qs[j, :, a1, n1] .- s.qs[j, :, a2, n2]) # nm
                result += harmonic_k_inter * R^2
            end
        end

        0.5 * result
    end
end

function harmonic_force_intra(n::Int)
    function term_f!(result::ForceContainer, j::Int,
                     s::State{P,D,A,N}) where {P,D,A,N}
        Rvec = Array{Float64}(undef, 3) # nm
        F = Array{Float64}(undef, 3) # kJ / nm mol

        for a1 in 1:A
            for a2 in (a1+1):A
                @. Rvec = s.qs[j, :, a1, n] - s.qs[j, :, a2, n]
                @. F = -harmonic_k_intra * Rvec
                @. result.f[j, :, a1, n] += F
                @. result.f[j, :, a2, n] -= F
            end
        end

        nothing
    end
end

function harmonic_force_inter(n1, n2)
    function term_f!(result::ForceContainer, j::Int,
                     s::State{P,D,A,N}) where {P,D,A,N}
        Rvec = Array{Float64}(undef, 3) # nm
        F = Array{Float64}(undef, 3) # kJ / nm mol

        for a1 in 1:A
            for a2 in 1:A
                @. Rvec = s.qs[j, :, a1, n1] - s.qs[j, :, a2, n2]
                @. F = -harmonic_k_inter * Rvec
                @. result.f[j, :, a1, n1] += F
                @. result.f[j, :, a2, n2] -= F
            end
        end

        nothing
    end
end

registered_models["harmonic"] = HarmonicWater

function harmonic_exact_slope(N::Int, T::Float64)
    # Force constant matrix (kJ / nm^2 mol).
    K = zeros(Float64, 3N, 3N)

    idx1 = 0

    for n1 in 1:N
        for a1 in 1:3
            idx1 += 1
            idx2 = 0

            for n2 in 1:N
                for a2 in 1:3
                    idx2 += 1

                    idx1 == idx2 && continue

                    if n1 == n2
                        K[idx2, idx1] += harmonic_k_intra
                    else
                        K[idx2, idx1] += harmonic_k_inter
                    end
                end
            end
        end
    end

    A = diagm(0 => sum(K; dims=2)[:, 1]) - K # kJ / nm^2 mol
    sqrtM = diagm(0 => repeat(sqrt.(MS), N)) # g^1/2 / mol^1/2
    sqrtMinv = diagm(0 => repeat(1 ./ sqrt.(MS), N)) # mol^1/2 / g^1/2

    eigvals, eigvecs = eigen(sqrtMinv * A * sqrtMinv) # 1 / ps^2, 1
    A_norm = norm(sqrtMinv * A * sqrtMinv) # 1 / ps^2
    # Indices of zero eigenvalues.
    zidxs = Set{Int}()

    for (i, eigval) in enumerate(eigvals)
        if abs(eigval) / A_norm < 1e-10
            push!(zidxs, i)
        end
    end

    # Rank of A.
    R = 3N - length(zidxs)

    omega = zeros(Float64, R, R) # mol / g nm^2
    omega_inv = zeros(Float64, R, R) # g nm^2 / mol
    idx = 0

    for i in 1:R
        idx += 1

        while idx in zidxs
            idx += 1
        end

        freq = sqrt(eigvals[idx]) # 1 / ps
        omega[i, i] += 2 * freq * tanh(0.5 * beta(T) * hbar * freq) / hbar
        omega_inv[i, i] += 1 / omega[i, i]
    end

    # g^1/2 / mol^1/2
    G = (eigvecs' * sqrtM)[setdiff(begin:end, zidxs), :]
    C_r = zeros(Float64, R) # g^1/2 / mol^1/2
    C_Q = zeros(Float64, R, 3N-1) # g^1/2 / mol^1/2

    for i in 1:3
        @. C_r += 0.5 * G[:, i]
        @. C_r -= 0.5 * G[:, i+3]
        @. C_Q[:, 1] += G[:, i]
        @. C_Q[:, 1] += G[:, i+3]
    end

    # Total mass of molecule (g / mol).
    M_mol = sum(MS)

    @. C_Q[:, 2] += (MS[2]/(MS[1] + MS[2]) * G[:, 1]
                     - MS[1]/(MS[1] + MS[2]) * G[:, 2])
    @. C_Q[:, 3] += (MS[3]/M_mol * G[:, 1]
                     + MS[3]/M_mol * G[:, 2]
                     - (MS[1] + MS[2])/M_mol * G[:, 3])
    @. C_Q[:, 4] += (MS[2]/(MS[1] + MS[2]) * G[:, 4]
                     - MS[1]/(MS[1] + MS[2]) * G[:, 5])
    @. C_Q[:, 5] += (MS[3]/M_mol * G[:, 4]
                     + MS[3]/M_mol * G[:, 5]
                     - (MS[1] + MS[2])/M_mol * G[:, 6])

    for i in 1:(3N-6)
        @. C_Q[:, i+5] += G[:, i+6]
    end

    v = omega * C_r # mol^1/2 / g^1/2 nm^2
    # 1 / nm^2
    slope = v' * (omega_inv - C_Q * pinv(C_Q' * omega * C_Q) * C_Q') * v

    beta(T), slope
end


text:
#=
Hackerrank problem - Interquartile Range
https://www.hackerrank.com/challenges/s10-interquartile-range/problem

Calculate the range between Q3 and Q1 (Q3 - Q1)

Date: 2020-06-30

Note: The code expands the array with the frequency, which is inefficient.
=#

using Printf
# Example values
# array = [6, 12, 8, 10, 20, 16]
# freq = [5, 4, 3, 2, 1, 5]
n = parse(Int, readline(stdin))
array = parse.(Int, split(readline(stdin)))
freq = parse.(Int, split(readline(stdin)))
array_with_freq = sort(collect(zip(array, freq)), by=first)
array = map((x) -> fill(x[1], x[2]), array_with_freq) |> Iterators.flatten |> collect

"""
  median(arr)

Calculate a median value in the way given in the problem.

# Arguments
- `arr`: An array of numbers
"""
function median(arr)
  n = size(arr)[1]
  n % 2 == 0 ? (arr[n ÷ 2] + arr[n ÷ 2 + 1]) / 2.0 : arr[n ÷ 2 + 1]
end

n = size(array)[1]
q1 = median(array[1:(n ÷ 2)])
q3 = median(array[(n % 2 == 0 ? n ÷ 2 + 1 : n ÷ 2 + 2):end])
println("$(@sprintf("%.1f", q3 - q1))")


text:
"""
Makes a list of actions.

`make_action_list(max_step, n_actions, headings; stay=true)`

`max_step` and `n_actions` determine the lateral movements

`headings::Vector{T<:Real}`
"""
function make_action_list(max_step, n_actions, headings; stay=true)

    # what we return
    actions = Action[]

    # loop over lateral actions
    angles = linspace(0.0, 360.0 - 360/n_actions, n_actions)
    for ang in angles
        ax = max_step * sind(ang)
        ay = max_step * cosd(ang)

        for h in headings
            push!(actions, (ax, ay, h*1.0))
        end
    end

    # include stay action
    if stay
        for h in headings
            push!(actions, (0.0, 0.0, h*1.0) )
        end
    end

    return actions
end


text:
# ----------
# Parameters and data structure for HD and MHD problem
# ----------

struct MVars{Aphys, Atrans, usr_var} <: MHDVars
    "x-component of velocity"
        ux :: Aphys
    "y-component of velocity"
        uy :: Aphys
    "z-component of velocity"
        uz :: Aphys
    "x-component of B-field"
        bx :: Aphys
    "y-component of B-field"
        by :: Aphys
    "z-component of B-field"
        bz :: Aphys
    "Fourier transform of x-component of velocity"
       uxh :: Atrans
    "Fourier transform of y-component of velocity"
       uyh :: Atrans
    "Fourier transform of z-component of velocity"
       uzh :: Atrans
    "Fourier transform of x-component of B-field"
       bxh :: Atrans
    "Fourier transform of y-component of B-field"
       byh :: Atrans
    "Fourier transform of z-component of B-field"
       bzh :: Atrans

    # Temperatory Cache
    "Non-linear term 1"
     nonlin1 :: Aphys
    "Fourier transform of Non-linear term"
    nonlinh1 :: Atrans

    # User Defined Vars
    "User Defined Vars"
    usr_vars :: usr_var
end

struct HVars{Aphys, Atrans, usr_var} <: MHDVars
  "x-component of velocity"
    ux :: Aphys
  "y-component of velocity"
    uy :: Aphys
  "z-component of velocity"
    uz :: Aphys
  "Fourier transform of x-component of velocity"
   uxh :: Atrans
  "Fourier transform of y-component of velocity"
   uyh :: Atrans
  "Fourier transform of z-component of velocity"
   uzh :: Atrans
  # Temperatory Cache
  "Non-linear term 1"
   nonlin1 :: Aphys
  "Fourier transform of Non-linear term"
   nonlinh1 :: Atrans

  # User Defined Vars
  "User Defined Vars"
    usr_vars :: usr_var
end

struct MHDParams_VP{Aphys,usr_param} <: AbstractParams

  "small-scale (hyper)-viscosity coefficient for v"
    ν :: Number
  "small-scale (hyper)-viscosity coefficient for b"
    η :: Number
  "(hyper)-viscosity order, `nν```≥ 1``"
    nν :: Int

  "Array Indexing for velocity"
    ux_ind :: Int
    uy_ind :: Int
    uz_ind :: Int

  "Array Indexing for B-field"
    bx_ind :: Int
    by_ind :: Int
    bz_ind :: Int

  "function that calculates the Fourier transform of the forcing, ``F̂``"
    calcF! :: Function

  "Volume penzlization method paramter"
    χ   :: Aphys
    U₀x :: Aphys
    U₀y :: Aphys
    U₀z :: Aphys
    B₀x :: Aphys
    B₀y :: Aphys
    B₀z :: Aphys

  "User defined params"
  usr_params :: usr_param

end

struct HDParams_VP{Aphys,usr_param} <: AbstractParams

  "small-scale (hyper)-viscosity coefficient for v"
    ν :: Number
  "(hyper)-viscosity order, `nν```≥ 1``"
    nν :: Int

  "Array Indexing for velocity"
    ux_ind :: Int
    uy_ind :: Int
    uz_ind :: Int

  "function that calculates the Fourier transform of the forcing, ``F̂``"
    calcF! :: Function

  "Volume penzlization method paramter"
    χ   :: Aphys
    U₀x :: Aphys
    U₀y :: Aphys
    U₀z :: Aphys

  "User defined params"
  usr_params :: usr_param

end

struct HDParams{usr_param} <: AbstractParams

  "small-scale (hyper)-viscosity coefficient for v"
   ν :: Number
  "(hyper)-viscosity order, `nν```≥ 1``"
   nν :: Int

  "Array Indexing for velocity"
   ux_ind :: Int
   uy_ind :: Int
   uz_ind :: Int

  "function that calculates the Fourier transform of the forcing, ``F̂``"
   calcF! :: Function

  "User defined params"
 usr_params :: usr_param
end

struct MHDParams{usr_param} <: AbstractParams

  "small-scale (hyper)-viscosity coefficient for v"
    ν :: Number
  "small-scale (hyper)-viscosity coefficient for b"
    η :: Number
  "(hyper)-viscosity order, `nν```≥ 1``"
    nν :: Int

  "Array Indexing for velocity"
    ux_ind :: Int
    uy_ind :: Int
    uz_ind :: Int

  "Array Indexing for B-field"
    bx_ind :: Int
    by_ind :: Int
    bz_ind :: Int
  "function that calculates the Fourier transform of the forcing, ``F̂``"
    calcF! :: Function

  "User defined params"
 usr_params :: usr_param

end

function SetMHDVars(::Dev, grid::AbstractGrid, usr_vars) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny, grid.nz) ux  uy  uz  bx  by bz nonlin1
  @devzeros Dev Complex{T} (grid.nkr, grid.nl, grid.nm) uxh uyh uzh bxh byh bzh nonlinh1

  return MVars( ux,  uy,  uz,  bx,  by,  bz,
              uxh, uyh, uzh, bxh, byh, bzh,
              nonlin1, nonlinh1, usr_vars);
end

function SetHDVars(::Dev, grid::AbstractGrid, usr_vars) where Dev
  T = eltype(grid)

  @devzeros Dev T (grid.nx, grid.ny, grid.nz) ux  uy  uz nonlin1
  @devzeros Dev Complex{T} (grid.nkr, grid.nl, grid.nm) uxh uyh uzh nonlinh1

  return HVars( ux,  uy,  uz,  uxh, uyh, uzh,
              nonlin1, nonlinh1, usr_vars);
end

# Functions of setting up the Vars and Params struct
function SetVars(dev, grid, usr_vars; B = false, VP = false)
  setvars = ifelse(B,SetMHDVars,SetHDVars);
  return setvars(dev, grid, usr_vars);
end

 function SetParams(::Dev, grid::AbstractGrid, calcF::Function, usr_params;
                     B = false, VP = false, ν = 0, η = 0, nν = 0) where Dev
  T = eltype(grid);
  usr_param = typeof(usr_params);
  if (B)
    if (VP)
      @devzeros Dev T (grid.nx, grid.ny, grid.nz) χ U₀x U₀y U₀z B₀x B₀y B₀z
      params = MHDParams_VP(ν, η, nν, 1, 2, 3, 4, 5, 6, calcF, χ, U₀x, U₀y, U₀z, B₀x, B₀y, B₀z, usr_params)
    else
      params = MHDParams(ν, η, nν, 1, 2, 3, 4, 5, 6, calcF, usr_params);
    end
  else
    if (VP)
      @devzeros Dev T (grid.nx, grid.ny, grid.nz) χ U₀x U₀y U₀z
      params = HDParams_VP(ν, nν, 1, 2, 3, calcF, χ, U₀x, U₀y, U₀z, usr_params);
    else
      params = HDParams(ν, nν, 1, 2, 3, calcF, usr_params);
    end
  end

  return params

end

max_stars_repo_name: MHDFlows/MHDFlows

max_issues_repo_name: MHDFlows/MHDFlows

max_forks_repo_name: MHDFlows/MHDFlows


text:
A = [true false false
     false true false
     false false true]
B = copy(A); B[1,2] = true;
sz_img = (3,3)
sz_img_3 = (3, 3, 3)
m, n = 3, 5

type_list = generate_test_types([Bool], [Gray])

@testset "GenericHausdorff" begin
    @info "test: GenericHausdorff"

    dist = Hausdorff()
    A_zero = zeros(Bool, 3, 3)
    @testset "Hausdorff" begin
        for T in type_list
            @testset "$T" begin
                test_Metric(dist, sz_img, T)
                test_ndarray(dist, sz_img_3, T)
                test_colwise(dist, n, sz_img, T)
                test_pairwise(dist, m, n, sz_img, T)

                a = A .|> T
                b = B .|> T
                test_numeric(dist, a, b, T; filename="references/Hausdorff_2d_Gray")
            end
        end
        @test dist(A_zero, A_zero) == 0.0
        @test dist(A_zero, A) == Inf
        @test dist(A, A_zero) == Inf
    end
    test_cross_type(dist, A, B, type_list)

    dist = ModifiedHausdorff()
    @testset "ModifiedHausdorff" begin
        for T in type_list
            @testset "$T" begin
                test_SemiMetric(dist, sz_img, T)
                test_ndarray(dist, sz_img_3, T)
                test_colwise(dist, n, sz_img, T)
                test_pairwise(dist, m, n, sz_img, T)

                a = A .|> T
                b = B .|> T
                test_numeric(dist, a, b, T; filename="references/ModifiedHausdorff_2d_Gray")
            end
        end
    end
    test_cross_type(dist, A, B, type_list)

    dist_2 = AxisWeightedHausdorff((1.0, 1.5))
    dist_3 = AxisWeightedHausdorff((1.0, 1.5, 2.0))
    @testset "AxisWeightedHausdorff" begin
        Aw = [false false false
            false true false
            false false true]
        Bw = [false true false
            false true false
            false false true]
        Cw = [false false false
            true true false
            false false true]
        @test dist_2(Aw, Bw) == 1.0
        @test dist_2(Aw, Cw) == 1.5

        for T in type_list
            @testset "$T" begin
                test_SemiMetric(dist_2, sz_img, T)
                test_ndarray(dist_3, sz_img_3, T)
                test_colwise(dist_2, n, sz_img, T)
                test_pairwise(dist_2, m, n, sz_img, T)
            end
        end
    end
    test_cross_type(dist_2, A, B, type_list)
end


text:
struct MockAbstractPotential <: AbstractPotential
    seed::Float64
end
function InteratomicPotentials.energy_and_force(s::AbstractSystem, p::MockAbstractPotential)
    (
        e=(p.seed)u"hartree",
        f=fill((@SVector[p.seed + 1, p.seed + 2, p.seed + 3])u"hartree/bohr", length(s))
    )
end
function InteratomicPotentials.virial_stress(s::AbstractSystem, p::MockAbstractPotential)
    (@SVector [p.seed - 1, p.seed - 2, p.seed - 3, p.seed - 4, p.seed - 5, p.seed - 6])u"hartree"
end

struct MockNonTrainablePotential <: NonTrainablePotential end

Base.@kwdef struct MockTrainablePotential <: TrainablePotential{NamedTuple{(:A, :B)},NamedTuple{(:C,)}}
    A::Float64
    B::Float64
    C::Float64
end

struct MockEmpiricalPotential <: EmpiricalPotential{NamedTuple{},NamedTuple{}}
    seed::Float64
    rcutoff::Float64
    species::Tuple
end
InteratomicPotentials.get_rcutoff(p::MockEmpiricalPotential) = p.rcutoff
InteratomicPotentials.get_species(p::MockEmpiricalPotential) = p.species
function InteratomicPotentials.potential_energy(R::AbstractFloat, p::MockEmpiricalPotential)
    R^2 * p.seed
end
function InteratomicPotentials.force(R::AbstractFloat, r::SVector{3}, p::MockEmpiricalPotential)
    R^2 * p.seed * r
end


text:
# """
# get source spray weights
# """
# function get_spray_weights!(
#     weights::AbstractArray, #denomI::AbstractArray,
#     ix1::AbstractArray{Int64},
#     ix2::AbstractArray{Int64},
#     iz1::AbstractArray{Int64},
#     iz2::AbstractArray{Int64},
#     mesh_x::AbstractArray,
#     mesh_z::AbstractArray,
#     xval::Float64,
#     zval::Float64,
# )

#     ix1[1], ix2[1] = indminn(mesh_x, xval, 2)
#     iz1[1], iz2[1] = indminn(mesh_z, zval, 2)

#     denomI = inv((mesh_x[ix2[1]] - mesh_x[ix1[1]]) * (mesh_z[iz2[1]] - mesh_z[iz1[1]]))
#     "for iz1, ix1"
#     weights[1] = (mesh_x[ix2[1]] - xval) * (mesh_z[iz2[1]] - zval) * denomI
#     "for iz1, ix2"
#     weights[2] = (xval - mesh_x[ix1[1]]) * (mesh_z[iz2[1]] - zval) * denomI
#     "for iz2, ix1"
#     weights[3] = (mesh_x[ix2[1]] - xval) * (zval - mesh_z[iz1[1]]) * denomI
#     "for iz2, ix2"
#     weights[4] = (xval - mesh_x[ix1[1]]) * (zval - mesh_z[iz1[1]]) * denomI

#     return weights, denomI, ix1, ix2, iz1, iz2
# end