text:
export GroupActions,
    SymmetricGroup, UniquePoints, search_in_radius, add!, multiplicities, unique_points

################
# Group actions
################
"""
    GroupActions(actions::Function...)

Store a bunch of group actions `(f1, f2, f3, ...)`.
Each action has to return a tuple.
The actions are applied in the following sense
1) f1 is applied on the original solution `s`
2) f2 is applied on `s` and the results of 1
3) f3 is applied on `s` and the results of 1) and 2)
and so on

## Example
```julia-repl
julia> f1(s) = (s * s,);

julia> f2(s) = (2s, -s, 5s);

julia> f3(s) = (s + 1,);

julia> GroupActions(f1)(3)
(3, 9)

julia> GroupActions(f1, f2)(3)
(3, 9, 6, -3, 15, 18, -9, 45)

julia> GroupActions(f1,f2, f3)(3)
(3, 9, 6, -3, 15, 18, -9, 45, 4, 10, 7, -2, 16, 19, -8, 46)
```
"""
struct GroupActions{T<:Tuple}
    actions::T
end
GroupActions(::Nothing) = GroupActions(())
GroupActions(actions::GroupActions) = actions
GroupActions(actions::Function...) = GroupActions(actions)
GroupActions(actions) = GroupActions(actions...)

function (actions::GroupActions)(s)
    S = [s]
    T = typeof(s)
    apply_actions(actions, s) do sᵢ
        sⱼ = convert(T, sᵢ)
        if sⱼ != s
            push!(S, sⱼ)
        end
        false
    end
    S
end

apply_actions(cb, action::GroupActions, s) = _apply_actions(action.actions, s, cb)
@inline function _apply_actions(actions::Tuple, x, cb::F) where {F}
    f, rest = first(actions), Base.tail(actions)
    y = f(x)
    if isa(x, AbstractVector{<:Number}) && isa(y, AbstractVector{<:Number})
        cb(y) && return true
        if _apply_actions(rest, y, cb)
            return true
        end
    else
        for yᵢ in f(x)
            cb(yᵢ) && return true
            if _apply_actions(rest, yᵢ, cb)
                return true
            end
        end
    end
    _apply_actions(rest, x, cb)
end
@inline _apply_actions(::Tuple{}, s, cb) = false

apply_actions(cb::G, actions::F, s) where {G,F} = actions(cb, s)

# Implemented group actions

"""
    SymmetricGroup(n)

Group action of the symmetric group S(n).
"""
struct SymmetricGroup
    permutations::Vector{Vector{Int}}
end
SymmetricGroup(N::Int) = SymmetricGroup(permutations(N))
function permutations(N::Int)
    s = Vector(1:N)
    perms = [copy(s)]
    while true
        i = N - 1
        while i >= 1 && s[i] >= s[i+1]
            i -= 1
        end
        if i > 0
            j = N
            while j > i && s[i] >= s[j]
                j -= 1
            end
            s[i], s[j] = s[j], s[i]
            reverse!(s, i + 1)
        else
            s[1] = N + 1
        end

        s[1] > N && break
        push!(perms, copy(s))
    end
    perms
end

Base.eltype(::Type{SymmetricGroup}) = Vector{Int}
Base.length(p::SymmetricGroup) = length(p.permutations)
Base.iterate(p::SymmetricGroup) = iterate(p.permutations)
Base.iterate(p::SymmetricGroup, s) = iterate(p.permutations, s)

#############
# UniquePoints
#############

"""
    UniquePoints{T, Id, M}

A data structure for assessing quickly whether a point is close to an indexed point as
determined by the given distances function `M`. The distance function has to be a *metric*.
The indexed points are only stored by their identifiers `Id`.

    UniquePoints(v::AbstractVector{T}, id::Id;
                    metric = EuclideanNorm(),
                    group_actions = nothing)

Initialize the data structure. This *does not* initialize the data structure with the point.

## Example

```julia
x = randn(ComplexF64, 4)
permutation(x) = ([x[2]; x[1]; x[3]; x[4]],)
group_actions = GroupActions(permutation)
X = group_actions(x)

# without group actions
unique_points = UniquePoints(x, 1)
HC.add!.(unique_points, X, 1:length(X), 1e-5)
length(unique_points) # 2

unique_points = UniquePoints(x, 1, group_actions = group_actions)
HC.add!.(unique_points, X, 1:length(X), 1e-5)
length(unique_points) # 1
```
"""
struct UniquePoints{T,Id,M,MaybeGA}
    tree::VoronoiTree{T,Id,M}
    group_actions::MaybeGA
    zero_vec::Vector{T}
end

function UniquePoints(
    v::AbstractVector,
    id;
    metric = EuclideanNorm(),
    group_action = nothing,
    group_actions = isnothing(group_action) ? nothing : GroupActions(group_action),
)
    if (group_actions isa Tuple) || (group_actions isa AbstractVector)
        group_actions = GroupActions(group_actions)
    end

    tree = VoronoiTree(v, id; metric = metric)
    UniquePoints(tree, group_actions, zeros(eltype(v), length(v)))
end

function Base.show(io::IO, UP::UniquePoints)
    print(io, typeof(UP), " with ", length(UP.tree), " points")
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", x::UniquePoints) = x
Base.length(UP::UniquePoints) = length(UP.tree)
Base.collect(UP::UniquePoints) = collect(UP.tree)
Base.broadcastable(UP::UniquePoints) = Ref(UP)
Base.empty!(UP::UniquePoints) = (Base.empty!(UP.tree); UP)

"""
    search_in_radius(unique_points, v, tol)

Search whether `unique_points` contains a point `p` with distances at most `tol` from `v`.
Returns `nothing` if no point exists, otherwise the identifier of `p` is returned.
"""
function search_in_radius(UP::UniquePoints{T,Id,M,GA}, v, tol::Real) where {T,Id,M,GA}
    id = search_in_radius(UP.tree, v, tol)
    if isnothing(id) && !isnothing(UP.group_actions)
        let actions = UP.group_actions::GA
            apply_actions(actions, v) do w
                id′ = search_in_radius(UP.tree, w, tol)
                if !isnothing(id′)
                    id = id′
                    return true
                end
                false
            end
        end
    end
end

"""
    add!(unique_points, v, id; atol = 1e-14, rtol = sqrt(eps()))
    add!(unique_points, v, id, atol)

Search whether `unique_points` contains a point `p` with distances at most
`max(atol, norm(v)rtol)` from `v`. If this is the case the identifier of `p` and `false` is
returned. Otherwise `(id, true)` is returned.
"""
function add!(UP::UniquePoints{T,Id,M,GA}, v, id::Id, tol::Real) where {T,Id,M,GA}
    found_id = search_in_radius(UP.tree, v, tol)
    if isnothing(found_id)
        if isnothing(UP.group_actions)
            insert!(UP.tree, v, id; use_distances = true)
            return (id, true)
        else
            let actions = UP.group_actions::GA
                apply_actions(actions, v) do w
                    found_id′ = search_in_radius(UP.tree, w, tol)
                    if !isnothing(found_id′)
                        found_id = found_id′
                        return true
                    end
                    false
                end
            end
            if isnothing(found_id)
                insert!(UP.tree, v, id)
                return (id, true)
            else
                return (found_id::Id, false)
            end
        end
    else
        return (found_id::Id, false)
    end
end
function add!(
    UP::UniquePoints{T,Id,M,GA},
    v,
    id::Id;
    atol::Float64 = 1e-14,
    rtol::Float64 = sqrt(eps()),
) where {T,Id,M,GA}
    n = UP.tree.metric(v, UP.zero_vec)
    rad = max(atol, rtol * n)
    add!(UP, v, id, rad)
end

####################
## Multiplicities ##
####################
"""
    multiplicities(vectors; metric = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)

Returns a `Vector{Vector{Int}}` `v`. Each vector `w` in 'v' contains all indices `i`,`j`
such that `w[i]` and `w[j]` have `distance` at most `max(atol, rtol * metric(0,w[i]))`.
The remaining `kwargs` are things that can be passed to [`UniquePoints`](@ref).

```julia-repl
julia> multiplicities([[1,0.5], [1,0.5], [1,1]])
[[1,2]]
```
This is the same as
```julia
multiplicities([[1,0.5], [1,0.5], [1,1]]; distance=(x,y) -> LinearAlgebra.norm(x-y))
```
Here is an example for using group actions.
```julia-repl
julia> X = [[1, 2, 3, 4], [2,1,3,4], [1,2,4,3], [2,1,4,3]]
julia> permutation(x) = [x[2], x[1], x[3], x[4]]
julia> m = multiplicities(X, group_action = permutation)
[[1,2], [3,4]]
```
"""
multiplicities(v; kwargs...) = multiplicities(identity, v; kwargs...)
function multiplicities(f::F, v; metric = EuclideanNorm(), kwargs...) where {F<:Function}
    isempty(v) && return Vector{Vector{Int}}()
    _multiplicities(f, v, metric; kwargs...)
end
function _multiplicities(
    f::F,
    V,
    metric;
    atol::Float64 = 1e-14,
    rtol::Float64 = 1e-8,
    kwargs...,
) where {F<:Function}
    unique_points = UniquePoints(f(first(V)), 1; metric = metric, kwargs...)
    mults = Dict{Int,Vector{Int}}()
    for (i, vᵢ) in enumerate(V)
        wᵢ = f(vᵢ)
        k, new_point = add!(unique_points, wᵢ, i; atol = atol, rtol = rtol)
        if !new_point
            if haskey(mults, k)
                push!(mults[k], i)
            else
                mults[k] = [k, i]
            end
        end
    end
    collect(values(mults))
end
"""
    unique_points(vectors; metric = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)

Returns all elements in `vector` for which two elements have `distance` at most `max(atol, rtol * metric(0,w[i]))`.
Note that the output can depend on the order of elements in vectors.
The remaining `kwargs` are things that can be passed to [`UniquePoints`](@ref).
"""
function unique_points(V; metric = EuclideanNorm(), atol = 1e-14, rtol = 1e-8, kwargs...)
    unique_points = UniquePoints(first(V), 1; metric = metric, kwargs...)
    out = Vector{eltype(V)}()
    for (i, vᵢ) in enumerate(V)
        _, new_point = add!(unique_points, vᵢ, i; atol = atol, rtol = rtol)
        if new_point
            push!(out, vᵢ)
        end
    end
    out
end


text:
# DistributionStat objects
# _fit! methods should only update "sufficient statistics"
# value methods should create the distribution

abstract DistributionStat{I<:Input} <: OnlineStat{I}
Base.mean(d::DistributionStat) = mean(value(d))
Base.var(d::DistributionStat) = var(value(d))
Base.std(d::DistributionStat) = std(value(d))
Ds.params(o::DistributionStat) = Ds.params(value(o))
Ds.ncategories(o::DistributionStat) = Ds.ncategories(value(o))
Base.cov(o::DistributionStat) = cov(value(o))

#--------------------------------------------------------------# Beta
type FitBeta{W<:Weight} <: DistributionStat{ScalarInput}
    value::Ds.Beta
    var::Variance{W}
end
FitBeta(wgt::Weight = EqualWeight()) = FitBeta(Ds.Beta(), Variance(wgt))
nobs(o::FitBeta) = nobs(o.var)
_fit!(o::FitBeta, y::Real, γ::Float64) = _fit!(o.var, y, γ)
function value(o::FitBeta)
    if nobs(o) > 1
        m = mean(o.var)
        v = var(o.var)
        α = m * (m * (1 - m) / v - 1)
        β = (1 - m) * (m * (1 - m) / v - 1)
        o.value = Ds.Beta(α, β)
    end
end
updatecounter!(o::FitBeta, n2::Int = 1) = updatecounter!(o.var, n2)
weight(o::FitBeta, n2::Int = 1) = weight(o.var, n2)

#------------------------------------------------------------------# Categorical
# Ignores weight
"""
Find the proportions for each unique input.  Categories are sorted by proportions.
Ignores `Weight`.

```julia
o = FitCategorical(y)
```
"""
type FitCategorical{T<:Any} <: DistributionStat{ScalarInput}
    value::Ds.Categorical
    d::Dict{T, Int}
    nobs::Int
end
FitCategorical() = FitCategorical(Ds.Categorical(1), Dict{Any, Int}(), 0)
function FitCategorical(y)
    o = FitCategorical()
    fit!(o, y)
    o
end
function _fit!(o::FitCategorical, y::Union{Real, AbstractString, Symbol}, γ::Float64)
    if haskey(o.d, y)
        o.d[y] += 1
    else
        o.d[y] = 1
    end
    o
end
# FitCategorical allows more than just Real input, so it needs special fit! methods
function fit!(o::FitCategorical, y::Union{AbstractString, Symbol})
    updatecounter!(o)
    γ = weight(o)
    _fit!(o, y, γ)
    o
end
function fit!{T <: Union{AbstractString, Symbol}}(o::FitCategorical, y::AVec{T})
    for yi in y
        fit!(o, yi)
    end
    o
end

sortpairs(o::FitCategorical) = sort(collect(o.d), by = x -> 1 / x[2])
function value(o::FitCategorical)
    if nobs(o) !== 0
        sortedpairs = sortpairs(o)
        counts = zeros(length(sortedpairs))
        for i in 1:length(sortedpairs)
            counts[i] = sortedpairs[i][2]
        end
        o.value = Ds.Categorical(counts / sum(counts))
    end
end
function Base.show(io::IO, o::FitCategorical)
    printheader(io, "FitCategorical ($(length(o.d)) levels)")
    print_item(io, "value", value(o))
    sortedpairs = sortpairs(o)
    print_item(io, "labels", [sortedpairs[i][1] for i in 1:length(sortedpairs)])
    print_item(io, "nobs", nobs(o))
end
updatecounter!(o::FitCategorical, n2::Int = 1) = (o.nobs += n2)
weight(o::FitCategorical, n2::Int = 1) = 0.0
nobs(o::FitCategorical) = o.nobs

#------------------------------------------------------------------# Cauchy
type FitCauchy{W<:Weight} <: DistributionStat{ScalarInput}
    value::Ds.Cauchy
    q::QuantileMM{W}
end
FitCauchy(wgt::Weight = LearningRate()) = FitCauchy(Ds.Cauchy(), QuantileMM(wgt))
_fit!(o::FitCauchy, y::Real, γ::Float64) = _fit!(o.q, y, γ)
nobs(o::FitCauchy) = nobs(o.q)
function value(o::FitCauchy)
    o.value = Ds.Cauchy(o.q.value[2], 0.5 * (o.q.value[3] - o.q.value[1]))
end
updatecounter!(o::FitCauchy, n2::Int = 1) = updatecounter!(o.q, n2)
weight(o::FitCauchy, n2::Int = 1) = weight(o.q, n2)

#------------------------------------------------------------------------# Gamma
# method of moments, TODO: look at Distributions for MLE
type FitGamma{W<:Weight} <: DistributionStat{ScalarInput}
    value::Ds.Gamma
    var::Variance{W}
end
FitGamma(wgt::Weight = EqualWeight()) = FitGamma(Ds.Gamma(), Variance(wgt))
_fit!(o::FitGamma, y::Real, γ::Float64) = _fit!(o.var, y, γ)
nobs(o::FitGamma) = nobs(o.var)
function value(o::FitGamma)
    m = mean(o.var)
    v = var(o.var)
    θ = v / m
    α = m / θ
    if nobs(o) > 1
        o.value = Ds.Gamma(α, θ)
    end
end
updatecounter!(o::FitGamma, n2::Int = 1) = updatecounter!(o.var, n2)
weight(o::FitGamma, n2::Int = 1) = weight(o.var, n2)

#-----------------------------------------------------------------------# LogNormal
type FitLogNormal{W<:Weight} <: DistributionStat{ScalarInput}
    value::Ds.LogNormal
    var::Variance{W}
end
FitLogNormal(wgt::Weight = EqualWeight()) = FitLogNormal(Ds.LogNormal(), Variance(wgt))
nobs(o::FitLogNormal) = nobs(o.var)
_fit!(o::FitLogNormal, y::Real, γ::Float64) = _fit!(o.var, log(y), γ)
function value(o::FitLogNormal)
    if nobs(o) > 1
        o.value = Ds.LogNormal(mean(o.var), std(o.var))
    end
end
updatecounter!(o::FitLogNormal, n2::Int = 1) = updatecounter!(o.var, n2)
weight(o::FitLogNormal, n2::Int = 1) = weight(o.var, n2)

#-----------------------------------------------------------------------# Normal
type FitNormal{W<:Weight} <: DistributionStat{ScalarInput}
    value::Ds.Normal
    var::Variance{W}
end
FitNormal(wgt::Weight = EqualWeight()) = FitNormal(Ds.Normal(), Variance(wgt))
nobs(o::FitNormal) = nobs(o.var)
_fit!(o::FitNormal, y::Real, γ::Float64) = _fit!(o.var, y, γ)
function value(o::FitNormal)
    if nobs(o) > 1
        o.value = Ds.Normal(mean(o.var), std(o.var))
    end
end
updatecounter!(o::FitNormal, n2::Int = 1) = updatecounter!(o.var, n2)
weight(o::FitNormal, n2::Int = 1) = weight(o.var, n2)

#-----------------------------------------------------------------------# Multinomial
type FitMultinomial{W<:Weight} <: DistributionStat{VectorInput}
    value::Ds.Multinomial
    means::Means{W}
end
function FitMultinomial(p::Integer, wgt::Weight = EqualWeight)
    FitMultinomial(Ds.Multinomial(p, ones(p) / p), Means(p, wgt))
end
nobs(o::FitMultinomial) = nobs(o.means)
_fit!{T<:Real}(o::FitMultinomial, y::AVec{T}, γ::Float64) = _fit!(o.means, y, γ)
function value(o::FitMultinomial)
    m = mean(o.means)
    o.value = Ds.Multinomial(round(Int, sum(m)), m / sum(m))
end
updatecounter!(o::FitMultinomial, n2::Int = 1) = updatecounter!(o.means, n2)
weight(o::FitMultinomial, n2::Int = 1) = weight(o.means, n2)

#---------------------------------------------------------------------# MvNormal
type FitMvNormal{W<:Weight} <: DistributionStat{VectorInput}
    value::Ds.MvNormal
    cov::CovMatrix{W}
end
function FitMvNormal(p::Integer, wgt::Weight = EqualWeight)
    FitMvNormal(Ds.MvNormal(zeros(p), eye(p)), CovMatrix(p, wgt))
end
nobs(o::FitMvNormal) = nobs(o.cov)
Base.std(d::FitMvNormal) = sqrt(var(d))  # No std() method from Distributions?
_fit!{T<:Real}(o::FitMvNormal, y::AVec{T}, γ::Float64) = _fit!(o.cov, y, γ)
function value(o::FitMvNormal)
    c = cov(o.cov)
    if isposdef(c)
        o.value = Ds.MvNormal(mean(o.cov), c)
    else
        warn("Covariance not positive definite.  More data needed.")
    end
end
updatecounter!(o::FitMvNormal, n2::Int = 1) = updatecounter!(o.cov, n2)
weight(o::FitMvNormal, n2::Int = 1) = weight(o.cov, n2)

#---------------------------------------------------------------------# convenience constructors
for nm in [:FitBeta, :FitGamma, :FitLogNormal, :FitNormal]
    eval(parse("""
        function $nm{T<:Real}(y::AVec{T}, wgt::Weight = EqualWeight())
            o = $nm(wgt)
            fit!(o, y)
            o
        end
    """))
end

function FitCauchy{T<:Real}(y::AVec{T}, wgt::Weight = LearningRate())
    o = FitCauchy(wgt)
    fit!(o, y)
    o
end

for nm in [:FitMultinomial, :FitMvNormal]
    eval(parse("""
        function $nm{T<:Real}(y::AMat{T}, wgt::Weight = EqualWeight())
            o = $nm(size(y, 2), wgt)
            fit!(o, y)
            o
        end
    """))
end

for nm in[:Beta, :Categorical, :Cauchy, :Gamma, :LogNormal, :Normal, :Multinomial, :MvNormal]
    eval(parse("""
        fitdistribution(::Type{Ds.$nm}, args...) = Fit$nm(args...)
    """))
end

#--------------------------------------------------------------# fitdistribution docs
"""
Estimate the parameters of a distribution.

```julia
using Distributions
# Univariate distributions
o = fitdistribution(Beta, y)
o = fitdistribution(Categorical, y)  # ignores Weight
o = fitdistribution(Cauchy, y)
o = fitdistribution(Gamma, y)
o = fitdistribution(LogNormal, y)
o = fitdistribution(Normal, y)
mean(o)
var(o)
std(o)
params(o)

# Multivariate distributions
o = fitdistribution(Multinomial, x)
o = fitdistribution(MvNormal, x)
mean(o)
var(o)
std(o)
cov(o)
```
"""
fitdistribution

ext: jl


text:
# Unit test functions
function test_bst_empty()
    t = BST{Int64, Int64}()
    @test isbst(t)

    return nothing
end

# should test isbst also

function test_bst_put()
    t = BST{Int64, Int64}()
    put!(t, 5, 1)
    @test isbst(t)
    put!(t, 7, 2)
    @test isbst(t)
    put!(t, 2, 3)
    @test isbst(t)
    @test t.root.key == 5
    @test t.root.right.key == 7
    @test t.root.left.key == 2

    return nothing
end

function test_bst_get()
    res = []

    t = BST{Int64, Int64}()

    put!(t, 5, 1)
    @test isbst(t)
    put!(t, 7, 2)
    @test isbst(t)
    put!(t, 2, 3)
    @test isbst(t)
    @test get!(t, 5) == 1
    @test get!(t, 7) == 2
    @test get!(t, 2) == 3

    return nothing
end

function test_bst_putget()
    n = 10^3
    c = rand((1:n), n)
    t = BST{Int64, Int64}()
    @test isbst(t)

    for i in c
        put!(t, i, 2i)
        @test isbst(t)
    end

    for i in sort(c)
        @test get!(t, i) == 2i
    end

    return nothing
end

function test_bst_delete()
    c = [5, 3, 4, 8, 7, 9]#, 50, 40, 30, 10, 20, 42, 43, 44, 60, 65, 66, 57, 56, 58, 64]
    n = length(c)

    for i in 0:n-1
        t = newt(c)
        for j in i:n+i
            delete!(t, c[mod(j, n)+1])
            @test isbst(t)
        end
    end

    return nothing
end

function newt(c)
    t = BST{Int64, Int64}()

    for i in c
        put!(t, i, 2i)
    end

    return t
end

function test_bst_maximum()
    n = 10^3
    c = rand((1:n), n)
    t = BST{Int64, Int64}()
    @test isbst(t)

    for i in c
        put!(t, i, 2i)
        @test isbst(t)
    end

    @test (maximum(t).key, maximum(t).value) == (maximum(c), 2maximum(c))

    return nothing
end

function test_bst_minimum()
    n = 10^3
    c = rand((1:n), n)
    t = BST{Int64, Int64}()
    @test isbst(t)

    for i in c
        put!(t, i, 2i)
        @test isbst(t)
    end

    @test (minimum(t).key, minimum(t).value) == (minimum(c), 2minimum(c))

    return nothing
end

function test_bst_maximum_node()
    n = 10^3
    c = rand((1:n), n)
    t = BST{Int64, Int64}()
    @test isbst(t)

    for i in c
        put!(t, i, 2i)
        @test isbst(t)
    end

    @test maximum_node(t.root, t.root)[1].key == maximum(c)

    return nothing
end

function test_bst_maximum_collect()
    n = 10^3
    c = shuffle(1:n)
    t = BST{Int64, Int64}()
    @test isbst(t)

    for i in c
        put!(t, i, 2i)
        @test isbst(t)
    end

    @test collect(t) == sort(c)

    return nothing
end

function test_bst_size()
    n = 6
    c = [5, 3, 4, 8, 7, 9]

    for i in 1:n
        t = newt(c)
        cn = setdiff(c, c[i])
        @test size(t) == n

        delete!(t, c[i])
        @test size(t) == n - 1
        @test isbst(t)
    end

    t = newt(c)
    delete!(t, 5)
    @test size(t.root.left) == 1
    @test size(t.root.right) == 3
    @test size(t.root.right.left) == 1
    @test size(t.root.right.right) == 1

    return nothing
end

# Run the unit tests
@testset "Binary search tree" begin
    test_bst_empty()
    test_bst_put()
    test_bst_get()
    test_bst_putget()
    test_bst_delete()
    test_bst_maximum()
    test_bst_minimum()
    test_bst_maximum_node()
    test_bst_maximum_collect()
    test_bst_size()
end


text:
using MultipleViewGeometry, Base.Test
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using MultipleViewGeometry.ModuleConstraints
using MultipleViewGeometry.ModuleConstruct
using MultipleViewGeometry.ModuleDraw
using MultipleViewGeometry.ModuleTriangulation
using BenchmarkTools, Compat
using StaticArrays
using  MAT, Plots

# Load MATLAB matrices that represent a pair of images and that contain
# a set of manually matched corresponding points.
file = matopen("./data/teapot.mat")
X = read(file,"pts3D")
# X is a 4 x N matrix, where N is the number of points expressed in homogeneous
# coordinates.
X = X[:,1:10:end]
close(file)

# Fix random seed.
srand(1234)
plotlyjs()

# Transform the point cloud so that it will be oriented upright and
# in front of two cameras which we will define next.
T1 = eye(4)
T1[:,4] = -mean(X,2)
T2 = [1 0 0 0;
      0 1 0 0;
      0 0 1 -2500;
      0 0 0 1]
R = eye(4)
R[1:3,1:3] = rotx(deg2rad(90))
X = T2*R*T1*X;

# Convert the 4 x N matrix into a list of 3D points expressed in homogeneous coordinates.
𝒳 = [Point3DH(X[1,i],X[2,i],X[3,i],1.0) for i = 1:size(X,2)]

# Intrinsic and extrinsic parameters of camera one.
𝐊₁ = eye(3)
𝐑₁ =  eye(3)
𝐭₁ = [0.0, 0.0, 0.0]
𝐜₁ = [0.0, 0.0, 0.0]
𝐄₁ = [𝐑₁ 𝐭₁]

# Intrinsic and extrinsic parameters of camera two.
𝐊₂ = eye(3)
𝐑₂ = eye(3)
𝐭₂ = [450.0,   0.0, 0.0]
𝐜₂ = [0.0, 0.0, 0.0]
𝐄₂ = [𝐑₂ 𝐭₂]

# Camera projection matrices.
𝐏₁ = construct(ProjectionMatrix(),𝐊₁,𝐑₁,𝐭₁)
𝐏₂ = construct(ProjectionMatrix(),𝐊₂,𝐑₂,𝐭₂)

# Set of corresponding points.
ℳ = project(Pinhole(),𝐏₁,𝒳)
ℳʹ = project(Pinhole(),𝐏₂,𝒳)

# Visualise the data points
p1 = Plots.plot(X[1,:],X[2,:],X[3,:],seriestype = :scatter, ms=1,grid = false, box = :none, legend = false)
draw!(WorldCoordinateSystem3D(), 1, p1)
draw!(Camera3D(), 𝐊₁, 𝐑₁, 𝐭₁, 1, p1)
draw!(Camera3D(), 𝐊₂, 𝐑₂, 𝐭₂, 1, p1)

# Plot the projections of the point cloud in the image pair.
p2 =  Plots.plot();
for n = 1:length(ℳ)
    m = ℳ[n]
    Plots.plot!([m[1]],[m[2]], grid = false, box = :none, legend = false,
                    seriestype = :scatter, ms = 2, markercolor=:Black,
                    aspect_ratio = :equal, yflip = true)
end

p3 =  Plots.plot();
for n = 1:length(ℳʹ)
    mʹ = ℳʹ[n]
    Plots.plot!([mʹ[1]],[mʹ[2]], grid = false, box = :none, legend = false,
                    seriestype = :scatter, ms = 2, markercolor=:Black,
                    aspect_ratio = :equal, yflip = true)
end

# Visualise the 3D point cloud, as well as the projected images.
l = @layout [ a; [b c] ]
p4 = Plots.plot(p1,p2, p3, layout = l)

# Estimate of the fundamental matrix and the true fundamental matrix.
𝐅 = estimate(FundamentalMatrix(), DirectLinearTransform(), (ℳ, ℳʹ))
𝐄 = construct(EssentialMatrix(), 𝐅,  𝐊₁, 𝐊₂)

𝐏₁, 𝐏₂ = construct(ProjectionMatrix(), 𝐄, (ℳ, ℳʹ))

𝒴 = triangulate(DirectLinearTransform(),𝐏₁,𝐏₂,(ℳ,ℳʹ))
Y = reinterpret(Float64,map(SVector{4,Float64},𝒴),(4,length(𝒴)))

# The point cloud can be reconsructed only up to an unknown scale factor.
p5 = Plots.plot(Y[1,:],Y[2,:],Y[3,:],seriestype = :scatter, ms=1,grid = false, box = :none, legend = false)
draw!(WorldCoordinateSystem3D(), 0.002, p5)
draw!(Camera3D(), 𝐊₁, 𝐏₁[1:3,1:3], 𝐏₁[:,4], 0.002, p5)
draw!(Camera3D(), 𝐊₂, 𝐏₂[1:3,1:3], 𝐏₂[:,4], 0.002, p5)

display(p4)
#display(p5)


text:
module test_libgit2_fetch

using Test
using LibGit2

@test hasmethod(LibGit2.GitRepo, Tuple{String})
@test hasmethod(LibGit2.fetch, Tuple{LibGit2.GitRepo})

end # module test_libgit2_fetch


text:
using KernelFunctions
using Zygote, ForwardDiff
using Test, LinearAlgebra
using FiniteDifferences
include("utils_AD.jl")
dims = [10,5]

A = rand(dims...)
B = rand(dims...)
K = [zeros(dims[1],dims[1]),zeros(dims[2],dims[2])]
kernels_noparams = [:SqExponentialKernel,:ExponentialKernel,:GammaExponentialKernel,
 :MaternKernel,:Matern32Kernel,:Matern52Kernel,
 :LinearKernel,:PolynomialKernel,
 :RationalQuadraticKernel,:GammaRationalQuadraticKernel,
 :ExponentiatedKernel]
l = 2.0
ds = [0.0,3.0]
vl = l*ones(dims[1])
testfunction(k,A,B) = det(kernelmatrix(k,A,B))
testfunction(k,A) = det(kernelmatrix(k,A))
ADs = [:Zygote,:ForwardDiff]

## Test kappa functions
@testset "Kappa functions" begin
    for AD in ADs
        @testset "$AD" begin
            for k in kernels_noparams
                for d in ds
                    @eval begin @test kappa_AD(Val(Symbol($AD)),$k(),$d) ≈ kappa_fdm($k(),$d) atol=1e-8 end
                end
            end
            # Linear -> C
            # Polynomial -> C,D
            # Gamma (etc) -> gamma
            #
        end
    end
end

@testset "Transform Operations" begin
    for AD in ADs
        @testset "$AD" begin
            @eval begin
            # Scale Transform
            transform_AD(Val(Symbol($AD)),ScaleTransform(l),A)
            # ARD Transform
            transform_AD(Val(Symbol($AD)),ARDTransform(vl),A)
            # LowRankTransform
            transform_AD(Val(Symbol($AD)),LowRankTransform(rand(2,10)),A)
            # Chain Transform
            # transform_AD(Val(Symbol($AD)),LowRankTransform,A)
            end
        end
    end
end

##TODO Eventually store real results in file
@testset "Zygote Automatic Differentiation test" begin
    @testset "ARD" begin
        for k in kernels
            @testset "$k" begin
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A,B),vl)[1], ForwardDiff.gradient(x->testfunction(k(x),A,B),vl)))
                @test  all(isapprox.(Zygote.gradient(x->testfunction(k(vl),x,B),A)[1],ForwardDiff.gradient(x->testfunction(k(vl),x,B),A)))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A),vl)[1],ForwardDiff.gradient(x->testfunction(k(x),A),vl)))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(vl),x),A)[1],ForwardDiff.gradient(x->testfunction(k(vl),x),A)))
            end
        end
    end
    @testset "ISO" begin
        for k in kernels
            @testset "$k" begin
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A,B),l)[1],ForwardDiff.gradient(x->testfunction(k(x[1]),A,B),[l])[1]))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(l),x,B),A)[1],ForwardDiff.gradient(x->testfunction(k(l),x,B),A)))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(x),A),l)[1],ForwardDiff.gradient(x->testfunction(k(x[1]),A),[l])))
                @test all(isapprox.(Zygote.gradient(x->testfunction(k(l),x),A)[1],ForwardDiff.gradient(x->testfunction(k(l[1]),x),A)))
            end
        end
    end
end

@testset "ForwardDiff AutomaticDifferentation test" begin
    @testset "ARD" begin
        for k in kernels
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x),A,B),vl)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(vl),x,B),A)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x),A),vl)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(vl),x),A)
        end
    end
    @testset "ISO" begin
        for k in kernels
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x[1]),A,B),[l])
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(l),x,B),A)
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(x[1]),A),[l])
            @test_nowarn ForwardDiff.gradient(x->testfunction(k(l[1]),x),A)
        end
    end
end

@testset "Tracker AutomaticDifferentation test" begin
    @testset "ARD" begin
        for k in kernels
            @test_broken all(Tracker.gradient(x->testfunction(k(x),A,B),vl)[1] .≈ ForwardDiff.gradient(x->testfunction(k(x),A,B),vl))
            @test_broken all(Tracker.gradient(x->testfunction(k(vl),x,B),A)[1] .≈ ForwardDiff.gradient(x->testfunction(k(vl),x,B),A))
            @test_broken all(Tracker.gradient(x->testfunction(k(x),A),vl)[1] .≈  ForwardDiff.gradient(x->testfunction(k(x),A),vl))
            @test_broken all.(Tracker.gradient(x->testfunction(k(vl),x),A) .≈ ForwardDiff.gradient(x->testfunction(k(vl),x),A))
        end
    end
    @testset "ISO" begin
        for k in kernels
            @test_broken Tracker.gradient(x->testfunction(k(x[1]),A,B),[l])
            @test_broken Tracker.gradient(x->testfunction(k(l),x,B),A)
            @test_broken Tracker.gradient(x->testfunction(k(x[1]),A),[l])
            @test_broken Tracker.gradient(x->testfunction(k(l),x),A)

        end
    end
end


text:
using Parameters
using Printf
using LinearAlgebra

export calc_radius_factor, sort_and_cut, get_channels

"""
Object containing the scattering channels and the asympotic energy
"""
struct ScatteringChannels
    "reciprical lattice points"
    g
    "reciprical lattice index"
    gi
    "kz2 after scattering"
    kz2::Array{Real}
end

function Base.show(ch::ScatteringChannels)
    @unpack g,gi,kz2 = ch
    println("g","\t","index","\t","kz²")
    for i in 1:length(g)
        @printf("%8.2f %8.2f %3d %3d %8.2f\n", g[i][1],g[i][2],gi[i][1],gi[i][1],kz2[i])
        end
end

"""
Calculate a hexagonal grid with a radius n and a lattice constant of 1
angle between a1 and a2 = 120 deg

# Arguments
- `n:Int64`
- `a₁=1`: base vector length 1
- `a₂=1`: base vector length 2

pattern (n = 3):
```
  + + +
 + + + +
+ + + + +
 + + + +
  + + +
```

# returns (g,gi)
 - g .. 2d array of lattice point coordinates
 - gi ... index of lattice point
"""
function hexgrid(n::Int64, a₁ = 1,a₂ = 1)
    a = [1.0 -0.5;
         0.0 sqrt(3.0)/2.0]

    max_index = 2*n+1

    nelements = 1+6*n+3*(n)*(n-1)
    g = zeros(Float64,nelements,2)
    gi = zeros(Int32,nelements,2)
    index = 1
    for i = 1:max_index
        maxj = i+n
        if maxj > max_index
            maxj = max_index
        end
        minj = i-n
        if minj < 1
            minj = 1
        end
        for j = minj:maxj
            t = a*[i-n-1 ; j-n-1]
            g[index,:]=t
            gi[index,:]=[i-n-1; j-n-1]
            index+=1
        end
    end
    return g,gi
end

"""
# Arguments
* lattice Lattice object which is used for the reciprical lattice vectors
* ki wavevector of incident beam
* theta angle of incident beam and the surface normal in radians
* n range of generated channels
"""
function get_channels(b1,b2,ki, theta, n::T = 10) where (T <: Int)
    g = []
    gi = []
    kz2= []
    for i in -n:n, j in -n:n
        push!(gi,[i,j])
        cur_g = b1.*i+b2.*j
        push!(g,cur_g)
        tmp1 = ki.*[sin(theta),0].+cur_g
        push!(kz2,ki^2-tmp1'tmp1)
    end
    ScatteringChannels(g,gi,kz2)
end

function filter_channels_radius(ch::ScatteringChannels,
                                max_norm)
    @unpack kz2,g,gi = ch
    new_g = []
    new_gi = []
    new_kz²= []
    for (i,cur_g) in enumerate(g)
        if norm(cur_g) <= max_norm
            push!(new_g,cur_g)
            push!(new_gi,gi[i])
            push!(new_kz²,kz²[i])
        end
    end
    ScatteringChannels(new_g,new_gi,new_kz²)
end

"""
Calculate radius of a circle fitting in 4 cells
cells are drawn by g = i*b1+j+b2 for i and j in -1:1
"""
function calc_radius_factor(lattice::Generic2DLattice)
    @unpack b1,b2 = lattice
    d1 = [b1[2],-b1[1]]./norm(b1)
    d2 = [b2[2],-b2[1]]./norm(b2)
    r1 = abs(d1'b2)
    r2 = abs(d2'b1)
    min(r1,r2)
end

"""
sort by kz² and cut away channels
# Arguments
* ch channels
* n_open number of open channels to keep
* n_closed number of closed channels to keep

"""
function sort_and_cut(ch::ScatteringChannels, n_open::T, n_closed::U) where {T <: Int, U <:Int}

    @unpack g,gi,kz2 = ch

    index = sortperm(collect(zip(kz2,getindex.(gi,1),getindex.(gi,2))),rev=true)
    sorted_kz2 = kz2[index]

    #find lowest open channel
    first_closed = findfirst(x -> x < 0, sorted_kz2)

    open_index = 1:first_closed-1
    if n_open >= 0 && n_open < length(open_index)
        open_index = (first_closed-n_open:first_closed-1)
    end

    close_index = first_closed:length(kz2)
    if n_closed >= 0 &&  n_closed < length(close_index)
        close_index = first_closed:first_closed+n_closed-1
    end

    selected_index = vcat(open_index,close_index)

    ScatteringChannels(g[index,:][selected_index],
                       gi[index,:][selected_index],
                       sorted_kz2[selected_index])

end


text:
"""
    eigTSD_Distance(P::Matrix{Float64}, 𝚽::Matrix{Float64}, 𝛌::Vector{Float64},
                    Q::SparseMatrixCSC{Int64,Int64}; length::Any = 1,
                    T::Any = :Inf, tol::Float64 = 1e-5)

computes the TSD distance matrix of P's column vectors on a graph.

# Input Argument
- `P::Matrix{Float64}`: vector measures with the same total mass 0.
- `𝚽::Matrix{Float64}`: matrix of the unweighted graph Laplacian eigenvectors.
- `𝛌::Vector{Float64}`: vector of eigenvalues.
- `Q::SparseMatrixCSC{Int64,Int64}`: the unweighted incidence matrix.
- `length::Any`: vector of edge lengths (default: `1` represents unweighted graphs)
- `T::Any`: the stopping time T in K_functional (default: `:Inf`)
- `tol::Float64`: tolerance for integral convergence (default: `1e-5`)

# Output Argument
- `dis::Matrix{Float64}`: distance matrix, d_{TSD}(φᵢ, φⱼ; T).

"""
function eigTSD_Distance(P::Matrix{Float64}, 𝚽::Matrix{Float64}, 𝛌::Vector{Float64},
                         Q::SparseMatrixCSC{Int64,Int64}; length::Any = 1,
                         T::Any = :Inf, tol::Float64 = 1e-5)
    N, ncols = Base.size(P)
    total_mass = sum(P, dims = 1)[:]
    if norm(total_mass - total_mass[1] * ones(ncols), Inf) > 10^4 * eps()
        @error("P's column measures do not share the same total mass.")
        return
    end
    # initialize the distance matrix
    dis = zeros(ncols, ncols)
    # store gradient of 𝚽 to avoid repeated computation
    ∇𝚽 = Q' * 𝚽

    for i = 1:(ncols - 1), j = (i + 1):ncols
        dis[i, j] = K_functional(P[:, i], P[:, j], 𝚽, ∇𝚽, 𝛌; length = length,
                                    T = T, tol = tol)[1]
    end
    return dis + dis'
end

"""
    K_functional(𝐩::Vector{Float64}, 𝐪::Vector{Float64}, 𝚽::Matrix{Float64},
                 ∇𝚽::Matrix{Float64}, 𝛌::Vector{Float64}; length::Any = 1,
                 T::Any = :Inf, dt::Float64 = 0.5/maximum(𝛌), tol::Float64 = 1e-5)

computes the K_functional between two vector meassures 𝐩 and 𝐪 on a graph.

# Input Argument
- `𝐩::Vector{Float64}`: the source vector measure.
- `𝐪::Vector{Float64}`: the destination vector measure.
- `𝚽::Matrix{Float64}`: matrix of the unweighted graph Laplacian eigenvectors.
- `∇𝚽::Matrix{Float64}`: gradient of unweighted graph Laplacian eigenvectors.
- `𝛌::Vector{Float64}`: vector of eigenvalues.
- `length::Any`: vector of edge lengths (default: 1 represents unweighted graphs)
- `T::Any`: the stopping time T in K_functional (default: :Inf)
- `tol::Float64`: tolerance for convergence (default: 1e-5)

# Output Argument
- `K::Float64`: TSD distance d_{TSD}(p, q; T).
- `E::Float64`: an estimated upper bound on the absolute error. In general,
    `E <= tol * norm(K)`.

"""
function K_functional(𝐩::Vector{Float64}, 𝐪::Vector{Float64}, 𝚽::Matrix{Float64},
                        ∇𝚽::Matrix{Float64}, 𝛌::Vector{Float64}; length::Any = 1,
                        T::Any = :Inf, tol::Float64 = 1e-5)
    if abs(sum(𝐩 - 𝐪)) > 10^4 * eps()
        @error("𝐩 and 𝐪 do not have the same total mass.")
    end
    f₀ = 𝐪 - 𝐩
    # store b to avoid repeated computation
    b = 𝚽' * f₀

    if T == :Inf
        # choose stopping time T such that exp(-λ₁T) = tol
        T = -log(tol) / 𝛌[2]
    end
    # use QuadGK.jl to numerically evaluate the integral
    K, E = quadgk(t -> weighted_1norm(∇f(t, ∇𝚽, b, 𝛌), length), 0, T, rtol=tol)
    return K, E
end

function ∇f(t, ∇𝚽, b, 𝛌)
    gu = ∇𝚽 * (exp.(-t * 𝛌) .* b)
    return gu
end

function weighted_1norm(x, length)
    if length == 1
        return norm(x, 1)
    else
        return sum(abs.(x) .* length)
    end
end


text:
@testset "SequentialQueryStrategies" begin

    dummy_data, _ = load_data(TEST_DATA_FILE)

    @testset "multi kde" begin
        data = hcat(fill([1, 2, 3, 4,], 3)...)
        @test_throws KDEException multi_kde(data)
        data2 = hcat(data, hcat(fill([1, 2, 3, 4,], 2)...))
        @test_throws KDEException multi_kde(data2)
        data3 = hcat(data2, rand(4,20) * 4)
        @test isa(OneClassActiveLearning.multi_kde(data3), PyObject)
    end

    params = Dict{Symbol, Any}()
    @testset "initialize_qs" begin
        @testset "Standard" begin
            qs_types = [RandomPQs]
            qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
            for qs in qs_objs
                for labels in [labelmap(fill(:U, 10)), labelmap(fill(:Lin, 10)), labelmap(fill(:Lout, 10))]
                    scores = qs_score(qs, dummy_data, labels)
                    @test length(scores) == size(dummy_data, 2)
                end
            end
        end

        qs_types = [MinimumMarginPQs, MinimumLossPQs]
        for x in qs_types
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict{Symbol, Any}())
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => nothing))
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => 1.5))
            @test_throws ArgumentError initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => -1))
        end
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, Dict(:p_inlier => 0.05)), qs_types)
        for qs in qs_objs
            @testset "DataBasedPQs $(typeof(qs))" begin
                @testset "only :U" begin
                        @test isdefined(qs, :bw_method)
                        @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:U, 10)))
                end

                @testset "only :Lin" begin
                    # number of :Lin must be sufficient that leave-one-out does not result in singular covariance matrix
                    scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 15)))
                    @test length(scores) == size(dummy_data, 2)
                end
                @testset "only :Lout" begin
                    @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:Lout, 10)))
                end
            end
        end
        @testset "DataBasedPQs GaussianKernel" begin
            classifier = SVDD.VanillaSVDD(dummy_data)
            init_strategy = SVDD.SimpleCombinedStrategy(FixedGammaStrategy(MLKernels.GaussianKernel(0.5)), FixedCStrategy(1))
            SVDD.initialize!(classifier, init_strategy)
            qs = initialize_qs(MinimumMarginPQs, classifier, dummy_data, Dict(:p_inlier => 0.05))
            scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 15)))
            @test length(scores) == size(dummy_data, 2)
        end

        qs_types = [ExpectedMinimumMarginPQs, ExpectedMaximumEntropyPQs]
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
        for qs in qs_objs
            @testset "DataBasedPQs $(typeof(qs))" begin
                @testset "only :U" begin
                        @test isdefined(qs, :bw_method)
                        @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:U, 10)))
                end

                @testset "only :Lin" begin
                    # number of :Lin must be sufficient that leave-one-out does not result in singular covariance matrix
                    scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 15)))
                    @test length(scores) == size(dummy_data, 2)
                end
                @testset "only :Lout" begin
                    @test_throws OneClassActiveLearning.MissingLabelTypeException qs_score(qs, dummy_data, labelmap(fill(:Lout, 10)))
                end
            end
        end

        qs_types = [RandomOutlierPQs, HighConfidencePQs, DecisionBoundaryPQs]
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
        for qs in qs_objs
            @testset "ModelBasedPQs $(typeof(qs))" begin
                scores = qs_score(qs, dummy_data, labelmap(fill(:Lin, 10)))
                @test length(scores) == size(dummy_data, 2)
            end
        end

        qs_types = [NeighborhoodBasedPQs, BoundaryNeighborCombinationPQs{SVDD.RandomOCClassifier}]
        qs_objs = map(x -> initialize_qs(x, SVDD.RandomOCClassifier(dummy_data), dummy_data, params), qs_types)
        for qs in qs_objs
            @testset "HybridPQs $(typeof(qs))" begin
                scores = qs_score(qs, dummy_data, labelmap(fill(:U, 10)))
                @test length(scores) == size(dummy_data, 2)
            end
        end

        @testset "Unknown Strategy" begin
            @test_throws ErrorException initialize_qs(Vector{Int}, SVDD.RandomOCClassifier(dummy_data), dummy_data, params)
        end

        @testset "multiple classifiers" begin

            pools = fill(:U, size(dummy_data, 2))
            init_strategy = SVDD.FixedParameterInitialization(MLKernels.GaussianKernel(0.5), 0.5)

            models = []
            model = SVDD.SVDDneg(dummy_data, pools)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            model = SVDD.VanillaSVDD(dummy_data)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            model = SVDD.SSAD(dummy_data, pools)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            model = SVDD.RandomOCClassifier(dummy_data)
            SVDD.initialize!(model, init_strategy)
            SVDD.fit!(model, TEST_SOLVER)
            push!(models, model)

            qs_types = [MinimumMarginPQs, MinimumLossPQs]
            for qst in qs_types
                for m in models
                    @testset "initialize_qs $qst, $(typeof(m))" begin
                        # check if qs is instantiated
                        qs = OneClassActiveLearning.initialize_qs(qst, model, dummy_data, Dict(:p_inlier => 0.05))
                    end
                end
            end

            qs_types = [ExpectedMinimumMarginPQs, ExpectedMaximumEntropyPQs] ∪ subtypes(ModelBasedPQs) ∪ [NeighborhoodBasedPQs, BoundaryNeighborCombinationPQs{SVDD.RandomOCClassifier}]
            for qst in qs_types
                for m in models
                    @testset "initialize_qs $qst, $(typeof(m))" begin
                        # check if qs is instantiated
                        qs = OneClassActiveLearning.initialize_qs(qst, model, dummy_data, params)
                    end
                end
            end
        end
    end

    @testset "get_query_objects" begin
        qs = OneClassActiveLearning.TestPQs()
        @testset "a" begin
            data = rand(2, 6)
            pools = fill(:U, 6)
            pools[5] = :Lin
            indices = [1, 3, 5, 7, 9, 11]
            history = [[7]]
            @test scores = OneClassActiveLearning.qs_score(qs, data, MLLabelUtils.labelmap(pools)) == collect(1:6)
            @test_throws ArgumentError OneClassActiveLearning.get_query_objects(qs, data, fill(:Lin, 5), indices, history)
            @test OneClassActiveLearning.get_query_objects(qs, data, pools, indices, history) == [11]
        end

        @testset "b" begin
            data = rand(2, 5)
            pools = fill(:U, 5)
            pools[5] = :Lin
            indices = [1, 2, 4, 7, 9]
            history = [[7]]
            @test scores = OneClassActiveLearning.qs_score(qs, data, MLLabelUtils.labelmap(pools)) == collect(1:5)
            @test_throws ArgumentError OneClassActiveLearning.get_query_objects(qs, data, fill(:Lin, 5), indices, history)
            @test OneClassActiveLearning.get_query_objects(qs, data, pools, indices, history) == [4]
        end
    end

    @testset "HybridQuerySynthesisPQs" begin
        params = Dict{Symbol, Any}(:qss_type => :RandomQss)
        qs = OneClassActiveLearning.initialize_qs(HybridQuerySynthesisPQs, SVDD.RandomOCClassifier(dummy_data), dummy_data, params)
        data = rand(2, 5)
        pools = fill(:U, 5)
        pools[5] = :Lin
        indices = [1, 2, 4, 7, 9]
        history = [[7]]
        @test_throws ArgumentError OneClassActiveLearning.get_query_objects(qs, data, fill(:Lin, 5), indices, history)
        query_index = OneClassActiveLearning.get_query_objects(qs, data, pools, indices, history)
        @test length(query_index) == 1
        @test query_index[1] ∈ indices
        @test query_index[1] ∉ history
    end
end


text:
using JuMP, EAGO

m = Model(EAGO.Optimizer)
set_optimizer_attribute(m, "mul_relax_style", 1)
set_optimizer_attribute(m, "verbosity", 4)
set_optimizer_attribute(m, "output_iterations", 1)
set_optimizer_attribute(m, "iteration_limit", 1000)
set_optimizer_attribute(m, "cut_max_iterations", 2)
set_optimizer_attribute(m, "subgrad_tighten", false)

# OBBT depth 0 -> 20... increases number of iterations...
set_optimizer_attribute(m, "obbt_depth", 8)
set_optimizer_attribute(m, "obbt_repetitions", 2)

# ----- Variables ----- #
x_Idx = Any[1,2]
@variable(m, x[x_Idx])
set_lower_bound(x[1], 0.1)
set_upper_bound(x[1], 0.9)

set_lower_bound(x[2], 0.1)
set_upper_bound(x[2], 0.9)

# ----- Constraints ----- #

@NLobjective(m, Min, (x[1]^2 - x[1])*(x[2]^2 - x[2]))
s = time()
optimize!(m)
@show termination_status(m)
sout = s - time()
@show solve_time(m)
@show sout

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from Euler angle to Euler angle and
#   axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angle_to_angleaxis

"""
    angle_to_angleaxis(θ₁::Number, θ₂::Number, θ₃::Number, rot_seq::Symbol = :ZYX)
    angle_to_angleaxis(Θ::EulerAngles)

Convert the Euler angles `θ₁`, `θ₂`, and `θ₃` [rad] with the rotation sequence
`rot_seq` to an Euler angle and axis representation.

Those values can also be passed inside the structure `Θ` (see
[`EulerAngles`](@ref)).

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```jldoctest
julia> angle_to_angleaxis(1, 0, 0, :XYZ)
EulerAngleAxis{Float64}:
  Euler angle : 1.0 rad  (57.2958°)
  Euler axis  : [1.0, 0.0, 0.0]

julia> Θ = EulerAngles(1, 1, 1, :XYZ);

julia> angle_to_angleaxis(Θ)
EulerAngleAxis{Float64}:
  Euler angle : 1.93909 rad  (111.102°)
  Euler axis  : [0.692363, 0.203145, 0.692363]
```
"""
@inline function angle_to_angleaxis(
    θ₁::Number,
    θ₂::Number,
    θ₃::Number,
    rot_seq::Symbol = :ZYX
)
    # First convert to DCM and then to Euler angle and axis.
    return dcm_to_angleaxis(angle_to_dcm(θ₁, θ₂, θ₃, rot_seq))
end

@inline function angle_to_angleaxis(Θ::EulerAngles)
    return angle_to_angleaxis(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)
end


text:
include("center.jl")
include("cross_validation.jl")
include("mse.jl")

function ridge( X, y, lambda::Number )

    colmeans_X = center!( X )
    mean_y = center!( y )

    #A = X'*X + diagm( fill( lambda, size(X,2) ) )
    #w = A\(X'*y)

    A = [ X ; fill( sqrt( lambda ), ( 1, size(X,2) ) ) ]
    B = [ y ; fill( 0, ( 1, size(y,2) ) ) ]
    w = A\B

    broadcast!( +, X, X, colmeans_X )
    broadcast!( +, y, y, mean_y )

    bias = mean_y - colmeans_X * w

    return w, bias

end

function ridge( X, y, ls, gen = Kfold( size(X,1), 5 ) )
    if Base.length( ls ) == 1
        return ridge( X, y, ls[1] )
    end
    l_errs = Array( Float64, size(ls,1) )
    for i in 1:size(ls,1)
        l = ls[i]
        train_f(is) = ridge( X[is,:], y[is,:], l )
        error_f(model,is) = mse( predict( model, X[is,:] ), y[is,:] )
        l_errs[i] = mean( cross_validate( train_f, error_f, gen ) )
    end
    l = ls[ indmin( l_errs ) ]
    return ridge( X, y, l ), l, l_errs
end

function predict( w, bias, X )
    return broadcast( +, X*w, bias )
end

function ridgeKernel( K, y, lambda )
    m = size( K, 1 )

    colmeans_K = center!( K )
    mean_y = center!( y )

    for i = 1:m
        K[i,i] += lambda
    end
    w = K \ y
    for i = 1:m
        K[i,i] -= lambda
    end

    broadcast!( +, y, y, mean_y )
    broadcast!( +, K, K, colmeans_K )

    bias = mean_y - colmeans_K * w

    return w, bias
end

function ridgeKernel!( K, y, lambda )
    m = size( K, 1 )

    colmeans_K = center!( K )
    mean_y = center!( y )

    for i = 1:m
        K[i,i] += lambda
    end
    w = K \ y

    broadcast!( +, y, y, mean_y )

    bias = mean_y - colmeans_K * w

    return w, bias
end

function predictKernel( w, bias, K )
    return broadcast( +, K*w, bias )
end

function ridge( kernel::Function, X, y, lambda::Number )

    m = size( X, 1 )
    K = Array( eltype(y), m, m )
    for i = 1:m
        for j = i:m
            K[i,j] = kernel( X[i,:], X[j,:] )
            K[j,i] = K[i,j]
        end
    end

    w,bias = ridgeKernel!( K, y, lambda )

    return w, bias, kernel, X

end

function ridge( kernel::Function, X, y, ls, gen = Kfold( size(X,1), 5 ) )

    if Base.length( ls ) == 1
        return ridge( X, y, ls[1], kernel )
    end

    m = size( X, 1 )
    K = Array( eltype(y), m, m )
    for i = 1:m
        for j = i:m
            K[i,j] = kernel( X[i,:], X[j,:] )
            K[j,i] = K[i,j]
        end
    end

    l_errs = Array( Float64, size(ls,1) )
    for i in 1:size(ls,1)
        l = ls[i]
        train_f(is) = ridgeKernel( K[is,is], y[is,:], l ), is

        function error_f(v,is)
            model, train_is = v
            yh = predictKernel( model..., K[is,train_is] )
            return mse( yh, y[is,:] )
        end

        l_errs[i] = mean( cross_validate( train_f, error_f, gen ) )
    end
    l = ls[ indmin( l_errs ) ]
    w,bias = ridgeKernel!( K, y, l )
    return (w, bias, kernel, X), l, l_errs
end

function predict( w, bias, kernel, X_model, X )
    m = size( X, 1 )
    n = size( X_model, 1 )
    K = Array( eltype(w), m, n )
    for i = 1:m
        for j = 1:n
            K[i,j] = kernel( X[i,:], X_model[j,:] )
        end
    end
    return predictKernel( w, bias, K )
end

function predict( model, X )
    predict( model..., X )
end

max_stars_repo_name: tbepler/julia-regression

max_issues_repo_name: tbepler/julia-regression

max_forks_repo_name: tbepler/julia-regression


text:
using FuzzySets, Test

@time begin
@time @testset "Fuzzy Sets" begin include("FuzzySets.jl") end
@time @testset "Fuzzy Numbers" begin include("FuzzyNumbers.jl") end
@time @testset "Fuzzy Vectors" begin include("FuzzyVectors.jl") end
@time @testset "Intervals" begin include("Intervals.jl") end
@time @testset "Fuzzy Weighted Averages" begin include("FuzzyWeightedAverages.jl") end
@time @testset "IT2-FCM" begin include("IT2FCM.jl") end
@time @testset "LFCM" begin include("LFCM.jl") end
@time @testset "Utilities" begin include("Utilities.jl") end
end


text:
expected_coskew = [ 0.147575   0.154872   0.126323   0.154872   0.414994  -0.033756   0.126323  -0.033756   0.013083
                    0.154872   0.414994  -0.033756   0.414994   0.861588  -0.301254  -0.033756  -0.301254  -0.287069
                    0.126323  -0.033756   0.013083  -0.033756  -0.301254  -0.287069   0.013083  -0.287069   1.017824]

expected_cokurt = [ 2.126762   1.118844   0.474779   1.118844   1.122935   0.187329   0.474779   0.187329   1.155241   1.118844   1.122935   0.187329   1.122935   1.404624  -0.026636   0.187329  -0.026636   0.276558   0.474779   0.187329  1.155241   0.187329  -0.026636   0.276558   1.155241   0.276558   0.178076
                    1.118844   1.122935   0.187329   1.122935   1.404624  -0.026636   0.187329  -0.026636   0.276558   1.122935   1.404624  -0.026636   1.404624   3.102880  -0.517198  -0.026636  -0.517198   0.779222   0.187329  -0.026636  0.276558  -0.026636  -0.517198   0.779222   0.276558   0.779222   0.218735
                    0.474779   0.187329   1.155241   0.187329  -0.026636   0.276558   1.155241   0.276558   0.178076   0.187329  -0.026636   0.276558  -0.026636  -0.517198   0.779222   0.276558   0.779222   0.218735   1.155241   0.276558  0.178076   0.276558   0.779222   0.218735   0.178076   0.218735   5.989492]

# Expected skewness, kurtosis, and central moments are calculated
# using StatsBase.
# (Note: + 3 has been added to the kurtosis values, to change the
# cumulant into a standardized central moment.)
expected_skewness = [0.2036496291131231, 1.055142530854644, 0.5466699866166419]
expected_kurtosis = [3.2674593510418695, 4.065496961897709, 2.6149221488093146]

expected_third_moment = [0.14757650128987837, 0.8615880265444518, 1.0178235320342266]
expected_fourth_moment = [2.1267754676390918, 3.1028802049533724, 5.989470624146371]


text:
"""
_determine_types(meshtype, fieldtype=Float64, facelen=3)

Given a subtype of AbstractMesh, determine the
type of vertex/point and face to use for internal computations.

Preference is given to the types specified by the Mesh call,
and will default to the `FieldType` for `SignedDistanceField`,
and Point{3,Float64}/Face{3,Int} for direct function sampling.
"""
function _determine_types(meshtype, fieldtype=Float64, facelen=3)
    # determine the point and face types
    # preference is given to the Mesh types
    # followed by SDF if unspecified
    if GT.vertextype(meshtype) !== Any
        VertType = GT.vertextype(meshtype)
    else
        VertType = GT.Point{3, fieldtype}
    end
    if GT.facetype(meshtype) !== Any
        FaceType = GT.facetype(meshtype)
    else
        FaceType = GT.Face{facelen, Int}
    end
    VertType, FaceType
end

function (::Type{MT})(df::GT.SignedDistanceField{3,ST,FT},
                      method::AbstractMeshingAlgorithm) where {MT <: GT.AbstractMesh, ST, FT}
    vertex_eltype = promote_type(FT, typeof(method.iso), typeof(method.eps))
    VertType, FaceType = _determine_types(MT, vertex_eltype, default_face_length(method))
    h = df.bounds
    vts, fcs = isosurface(df.data, method, VertType, FaceType, origin=VertType(GT.origin(h)), widths=VertType(GT.widths(h)))
    MT(vts, fcs)::MT
end

function (::Type{MT})(f::Function,
                      h::GT.HyperRectangle,
                      samples::NTuple{3,T},
                      method::AbstractMeshingAlgorithm)where {MT <: GT.AbstractMesh, T <: Integer}
    vertex_eltype = promote_type(T, typeof(method.iso), typeof(method.eps))
    VertType, FaceType = _determine_types(MT,vertex_eltype, default_face_length(method))
    vts, fcs = isosurface(f, method, VertType, FaceType, samples=samples, origin=VertType(GT.origin(h)), widths=VertType(GT.widths(h)))
    MT(vts, fcs)::MT
end

function (::Type{MT})(f::Function,
                      h::GT.HyperRectangle,
                      method::AbstractMeshingAlgorithm;
                      samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {MT <: GT.AbstractMesh, T <: Integer}
    MT(f, h, samples, method)
end

function (::Type{MT})(volume::AbstractArray{T, 3},
                      method::AbstractMeshingAlgorithm;
                      kwargs...) where {MT <: GT.AbstractMesh, T}
    vertex_eltype = promote_type(T, typeof(method.iso), typeof(method.eps))
    VertType, FaceType = _determine_types(MT,vertex_eltype, default_face_length(method))
    vts, fcs = isosurface(volume, method, VertType, FaceType; kwargs...)
    MT(vts, fcs)::MT
end


text:
module TestSetup

using WhiskerTracking

end


text:
using JLD
using DiffPrecTest
using DataFrames, CSV
using StatsPlots
using Distributions, Random, LinearAlgebra

function coverage(res::Vector{BootstrapEstimates}, α::Float64=0.95)
    NUM_REP = length(res)
    coverage = 0
    for rep=1:NUM_REP
        CI = simulCI(res[rep], α)
        coverage += all(0 .<= CI[:, 2]) * all(0 .>= CI[:, 1]) ? 1 : 0
    end
    coverage / NUM_REP
end

function coverage_studentized(res::Vector{BootstrapEstimates}, α::Float64=0.95)
    NUM_REP = length(res)
    coverage = 0
    for rep=1:NUM_REP
        CI = simulCIstudentized(res[rep], α)
        coverage += all(0 .<= CI[:, 2]) * all(0 .>= CI[:, 1]) ? 1 : 0
    end
    coverage / NUM_REP
end

function coverage_file_studentized(folderName, α::Float64=0.95;  NUM_REP::Int64=1000)
    coverage = 0
    for rep=1:NUM_REP
        res = load("$(folderName)_$(rep).jld", "boot_res")
        CI = simulCIstudentized(res, α)
        coverage += all(0 .<= CI[:, 2]) * all(0 .>= CI[:, 1]) ? 1 : 0
    end
    coverage / NUM_REP
end

function coverage_file(folderName, α::Float64=0.95;  NUM_REP::Int64=1000)
    coverage = 0
    for rep=1:NUM_REP
        res = load("$(folderName)_$(rep).jld", "boot_res")
        CI = simulCI(res, α)
        coverage += all(0 .<= CI[:, 2]) * all(0 .>= CI[:, 1]) ? 1 : 0
    end
    coverage / NUM_REP
end

coverageTable = zeros(3, 4)
sim = [8,9,10,11]
setting = [1,2,3]

for iSim = 1:4
    for iSetting = 1:3
        @show iSim, iSetting
        global sim
        global setting
        global coverageTable
        s1 = sim[iSim]
        s2 = setting[iSetting]
        folderName = "/scratch/midway2/mkolar/diffTest/sim$(s1)_o/res_$(s2)"
        coverageTable[iSetting, iSim] = coverage_file(folderName)
    end
end

@show coverageTable


text:
# Load packages
using Genesys, JLD, Dates

#=
First, you have to define your own designer and controller type. Make sure that the genesys "supertype" is added.
=#

# Define your own designer...
mutable struct bar <: Genesys.AbstractDesigner
    decisions::NamedTuple
    bar() = new()
end

# ...and controller
mutable struct foo <: Genesys.AbstractController
    decisions::NamedTuple
    foo() = new()
end

# Define the "offline" functions for both the designer and controller
function Genesys.initialize_controller!(mg::Microgrid, controller::foo, ω::Scenarios)
    # Preallocate
    Genesys.preallocate!(mg, controller)
    return controller
end

function Genesys.initialize_designer!(mg::Microgrid, designer::bar, ω::Scenarios)
    # Preallocate
    Genesys.preallocate!(mg, designer)
    return designer
end

# Define the "online" functions for both the designer and controller
function Genesys.compute_operation_decisions!(h::Int64, y::Int64, s::Int64, mg::Microgrid, controller::foo)
    return controller
end

function Genesys.compute_investment_decisions!(y::Int64, s::Int64, mg::Microgrid, designer::bar)
    return designer
end

#=
Let's simulate the microgrid with the dummies controller and designer...
=#

# Parameters of the simulation
const nh, ny, ns = 8760, 2, 1

# Load input data
data = load(joinpath("data","ausgrid_5_twostage.jld"))

# Initialize the microgrid
microgrid = Microgrid(parameters = GlobalParameters(nh, ny, ns))

# Add the equipment to the microgrid
add!(microgrid, Demand(carrier = Electricity()), Solar(), Liion(), Grid(carrier = Electricity()))

# Initialize scenarios
ω_d, ω_a = Scenarios(microgrid, data["ω_optim"]), Scenarios(microgrid, data["ω_simu"])

# Initialize designer
designer = initialize_designer!(microgrid, bar(), ω_d)

# Initialize controller
controller = initialize_controller!(microgrid, foo(), ω_d)

# Assessment
simulate!(microgrid, controller, designer, ω_a)


text:
# an implementation of Readout must implement the liquidState method which returns a VecF
abstract Readout

# default is to do nothing on update
OnlineStats.fit!{T <: SpikingNeuron}(readout::Readout, neurons::AbstractVector{T}) = nothing

# ----------------------------------------------------

# simple readout of the current firings in the output neurons
type FireReadout <: Readout end

function liquidState(lsm, readout::FireReadout)
  Float64[float(neuron.fired) for neuron in lsm.liquid.outputNeurons]
end

# ----------------------------------------------------

# takes the state of the neuron, whatever that may be
type StateReadout <: Readout end

function liquidState(lsm, readout::StateReadout)
  Float64[float(state(neuron)) for neuron in lsm.liquid.outputNeurons]
end

# ----------------------------------------------------

# keep a history of firings for each neuron
type FireWindowReadout <: Readout
  lookbackWindow::Int
  firings::Vector{CircularBuffer{Float64}}  # one CB per readout neuron
  FireWindowReadout(lookbackWindow::Integer) = new(lookbackWindow, CircularBuffer{Float64}[])
end

function OnlineStats.fit!{T <: SpikingNeuron}(readout::Readout, neurons::AbstractVector{T})

  # initialize on first pass
  if isempty(readout.firings)
    readout.firings = [CircularBuffer(Float64, readout.lookbackWindow) for i in 1:length(neurons)]
  end

  # add to window
  for (i,neuron) in enumerate(neurons)
    push!(readout.firings[i], Float64(neuron.fired))
  end
end

function liquidState(lsm, readout::FireWindowReadout)
  Float64[sum(cb) for cb in readout.firings]
end


text:
using QuadGK

function BlackScholes(S, T, K, r, σ)
    f(x) = exp(-0.5*x^2)
    N(d) = 1/sqrt(2*pi) * quadgk(f, -Inf, d)[1]
    d1 = ( log(S/K) + (r + 0.5 * σ^2)*(T))/(σ*sqrt(T))
    d2 = d1 - σ*sqrt(T)
    function call()
        return  S * N(d1) - K * exp(-r*(T))*N(d2)
    end
    function put()
        return K*exp(-r*(T))*N(-d2) - S*N(-d1)
    end
    return call(), put()
end

function BlackScholesDividend(S, T, K, r, σ, q)
    f(x) = exp(-0.5*x^2)
    N(d) = 1/sqrt(2*pi) * quadgk(f, -Inf, d)[1]
    d1 = ( log(S/K) + (r - q + 0.5 * σ^2)*(T))/(σ*sqrt(T))
    d2 = d1 - σ*sqrt(T)
    function call()
        return S*exp(-q*(T))*N(d1) - K*exp(-r*(T))*N(d2)
    end
    function put()
        return K*exp(-r*(T))*N(-d2) - S*exp(-q*(T))*N(d1)
    end
    return call(), put()
end


text:
using Plots

function PD(;
    data,
    complex
    )

    if complex=="Vietoris-Rips"
    end

    if complex=="Cech"
    end
end


text:
"""
    MRCExtendedHeader{T}

Extended header of an MRC file.

    MRCExtendedHeader(data)

Store `data` in an extended header. `data` must be directly writeable to a file with
`write(::IO, data)`.

    MRCExtendedHeader()

Create an extended header whose data is an empty `Vector{UInt8}`.
"""
mutable struct MRCExtendedHeader{T}
    data::T
end
MRCExtendedHeader() = MRCExtendedHeader(UInt8[])

Base.sizeof(eh::MRCExtendedHeader) = sizeof(eh.data)


text:
@testset "expected_loglik" begin
    # Test that the various methods of computing expectations return the same
    # result.
    rng = MersenneTwister(123456)
    q_f = Normal.(zeros(10), ones(10))

    likelihoods_to_test = [
        ExponentialLikelihood(),
        GammaLikelihood(),
        PoissonLikelihood(),
        GaussianLikelihood(),
    ]

    @testset "testing all analytic implementations" begin
        # Test that we're not missing any analytic implementation in `likelihoods_to_test`!
        implementation_types = [
            (; quadrature=m.sig.types[2], lik=m.sig.types[5]) for
            m in methods(ApproximateGPs.expected_loglik)
        ]
        analytic_likelihoods = [
            m.lik for m in implementation_types if
            m.quadrature == ApproximateGPs.Analytic && m.lik != Any
        ]
        for lik_type in analytic_likelihoods
            @test any(lik isa lik_type for lik in likelihoods_to_test)
        end
    end

    @testset "$(nameof(typeof(lik)))" for lik in likelihoods_to_test
        methods = [GaussHermite(100), MonteCarlo(1e7)]
        def = ApproximateGPs._default_quadrature(lik)
        if def isa Analytic
            push!(methods, def)
        end
        y = rand.(rng, lik.(zeros(10)))

        results = map(m -> ApproximateGPs.expected_loglik(m, y, q_f, lik), methods)
        @test all(x -> isapprox(x, results[end]; atol=1e-6, rtol=1e-3), results)
    end

    @test ApproximateGPs.expected_loglik(
        MonteCarlo(), zeros(10), q_f, GaussianLikelihood()
    ) isa Real
    @test ApproximateGPs.expected_loglik(
        GaussHermite(), zeros(10), q_f, GaussianLikelihood()
    ) isa Real
    @test ApproximateGPs._default_quadrature(θ -> Normal(0, θ)) isa GaussHermite

    @testset "testing Zygote compatibility with GaussHermite" begin # see issue #82
        N = 10
        gh = GaussHermite(12)
        μs = randn(rng, N)
        σs = rand(rng, N)
        # Test differentiation with variational parameters
        for lik in likelihoods_to_test
            y = rand.(rng, lik.(rand.(Normal.(μs, σs))))
            gμ, glogσ = Zygote.gradient(μs, log.(σs)) do μ, logσ
                ApproximateGPs.expected_loglik(gh, y, Normal.(μ, exp.(logσ)), lik)
            end
            @test all(isfinite, gμ)
            @test all(isfinite, glogσ)
        end
        # Test differentiation with likelihood parameters
        # Test GaussianLikelihood parameter
        σ = 1.0
        y = randn(rng, N)
        glogσ = only(
            Zygote.gradient(log(σ)) do x
                ApproximateGPs.expected_loglik(
                    gh, y, Normal.(μs, σs), GaussianLikelihood(exp(x))
                )
            end,
        )
        @test isfinite(glogσ)
        # Test GammaLikelihood parameter
        α = 2.0
        y = rand.(rng, Gamma.(α, rand(N)))
        glogα = only(
            Zygote.gradient(log(α)) do x
                ApproximateGPs.expected_loglik(
                    gh, y, Normal.(μs, σs), GammaLikelihood(exp(x))
                )
            end,
        )
        @test isfinite(glogα)
    end
end


text:
"""
    VelocityFields(arch, grid; [u, v, w])

Return a NamedTuple with fields `u`, `v`, `w` initialized on the architecture `arch`
and `grid`. Fields may be passed via the optional keyword arguments `u`, `v`, and `w`.
"""
function VelocityFields(arch, grid; u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)),
                                    v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)),
                                    w = ZFaceField(arch, grid, WVelocityBoundaryConditions(grid)))
    return (u=u, v=v, w=w)
end

"""
    VelocityFields(arch, grid, bcs::NamedTuple)

Return a NamedTuple with fields `u`, `v`, `w` initialized on the architecture `arch`
and `grid`. Boundary conditions `bcs` may be specified via a named tuple of
`FieldBoundaryCondition`s.
"""
function VelocityFields(arch, grid, bcs::NamedTuple)
    u_bcs = :u ∈ keys(bcs) ? bcs[:u] : UVelocityBoundaryConditions(grid)
    v_bcs = :v ∈ keys(bcs) ? bcs[:v] : VVelocityBoundaryConditions(grid)
    w_bcs = :w ∈ keys(bcs) ? bcs[:w] : WVelocityBoundaryConditions(grid)

    u = XFaceField(arch, grid, u_bcs)
    v = YFaceField(arch, grid, v_bcs)
    w = ZFaceField(arch, grid, w_bcs)

    return (u=u, v=v, w=w)
end

"""
    TracerFields(arch, grid, tracer_names; kwargs...)

Return a NamedTuple with tracer fields specified by `tracer_names` initialized as
`CellField`s on the architecture `arch` and `grid`. Fields may be passed via optional
keyword arguments `kwargs` for each field.

# Examples
```julia
arch = CPU()
topology = (Periodic, Periodic, Bounded)
grid = RegularCartesianGrid(topology=topology, size=(16, 16, 16), size=(1, 2, 3))
tracers = (:T, :S, :random)
noisy_field = CellField(arch, grid, TracerBoundaryConditions(grid), randn(16, 16))
tracer_fields = TracerFields(arch, grid, tracers, random=noisy_field)
```
"""
function TracerFields(arch, grid, names; kwargs...)
    tracer_names = tracernames(names) # filter `names` if it contains velocity fields
    tracer_fields =
        Tuple(c ∈ keys(kwargs) ?
              kwargs[c] :
              CellField(arch, grid, TracerBoundaryConditions(grid))
              for c in tracer_names)
    return NamedTuple{tracer_names}(tracer_fields)
end

"""
    TracerFields(arch, grid, tracer_names, bcs)

Return a NamedTuple with tracer fields specified by `tracer_names` initialized as
`CellField`s on the architecture `arch` and `grid`. Boundary conditions `bcs` may
be specified via a named tuple of `FieldBoundaryCondition`s.
"""
function TracerFields(arch, grid, names, bcs)
    tracer_names = tracernames(names) # filter `names` if it contains velocity fields
    tracer_fields =
        Tuple(c ∈ keys(bcs) ?
              CellField(arch, grid, bcs[c]) :
              CellField(arch, grid, TracerBoundaryConditions(grid))
              for c in tracer_names)
    return NamedTuple{tracer_names}(tracer_fields)
end

TracerFields(arch, grid, ::Union{Tuple{}, Nothing}, args...; kwargs...) = NamedTuple()
TracerFields(arch, grid, tracer::Symbol; kwargs...) = TracerFields(arch, grid, tuple(tracer); kwargs...)

"""
    TracerFields(arch, grid, tracer_fields::NamedTuple; kwargs...)

Convenience method for restoring checkpointed models that returns the already-instantiated
`tracer_fields` with non-default boundary conditions.
"""
function TracerFields(arch, grid, proposed_tracer_fields::NamedTuple, bcs; kwargs...)
    grid = proposed_tracer_fields[1].grid

    tracer_fields =
        Tuple(c ∈ keys(bcs) ?
              Field{Cell, Cell, Cell}(proposed_tracer_fields[c].data, grid, bcs[c]) :
              proposed_tracer_fields[c]
              for c in tracernames(proposed_tracer_fields))

    return NamedTuple{tracernames(proposed_tracer_fields)}(tracer_fields)
end

"Shortcut constructor for empty tracer fields."
TracerFields(arch, grid, empty_tracer_fields::NamedTuple{(),Tuple{}}, args...; kwargs...) = NamedTuple()

"Returns true if the first three elements of `names` are `(:u, :v, :w)`."
has_velocities(names) = :u == names[1] && :v == names[2] && :w == names[3]

# Tuples of length 0-2 cannot contain velocity fields
has_velocities(::Tuple{}) = false
has_velocities(::Tuple{X}) where X = false
has_velocities(::Tuple{X, Y}) where {X, Y} = false

tracernames(::Nothing) = ()
tracernames(name::Symbol) = tuple(name)
tracernames(names::NTuple{N, Symbol}) where N = has_velocities(names) ? names[4:end] : names
tracernames(::NamedTuple{names}) where names = tracernames(names)

"""
    PressureFields(arch, grid; [pHY′, pNHS])

Return a NamedTuple with pressure fields `pHY′` and `pNHS` initialized as
`CellField`s on the architecture `arch` and `grid`. Fields may be passed via the
optional keyword arguments `pHY′` and `pNHS`.
"""
function PressureFields(arch, grid; pHY′ = CellField(arch, grid, PressureBoundaryConditions(grid)),
                                    pNHS = CellField(arch, grid, PressureBoundaryConditions(grid)))

    return (pHY′=pHY′, pNHS=pNHS)
end

"""
    PressureFields(arch, grid, bcs::NamedTuple)

Return a NamedTuple with pressure fields `pHY′` and `pNHS` initialized as
`CellField`s on the architecture `arch` and `grid`.  Boundary conditions `bcs` may
be specified via a named tuple of `FieldBoundaryCondition`s.
"""
function PressureFields(arch, grid, bcs::NamedTuple)
    pHY′_bcs = :pHY′ ∈ keys(bcs) ? bcs[:pHY′] : PressureBoundaryConditions(grid)
    pNHS_bcs = :pNHS ∈ keys(bcs) ? bcs[:pNHS] : PressureBoundaryConditions(grid)

    pHY′ = CellField(arch, grid, pHY′_bcs)
    pNHS = CellField(arch, grid, pNHS_bcs)

    return (pHY′=pHY′, pNHS=pNHS)
end

"""
    TendencyFields(arch, grid, tracer_names; kwargs...)

Return a NamedTuple with tendencies for all solution fields (velocity fields and
tracer fields), initialized on the architecture `arch` and `grid`. Optional `kwargs`
can be specified to assign data arrays to each tendency field.
"""
function TendencyFields(arch, grid, tracer_names;
                        u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid)),
                        v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid)),
                        w = ZFaceField(arch, grid, WVelocityBoundaryConditions(grid)),
                        kwargs...)

    velocities = (u=u, v=v, w=w)
    tracers = TracerFields(arch, grid, tracer_names; kwargs...)
    return merge(velocities, tracers)
end


text:
import Base: show, ==, push!, append!

## nice-valued ranges for histograms
function histrange{T<:FloatingPoint}(v::AbstractArray{T}, n::Integer, closed::Symbol)
    if length(v) == 0
        return 0.0:1.0:0.0
    end
    lo, hi = extrema(v)
    if hi == lo
        step = 1.0
    else
        bw = (hi - lo) / n
        e = 10.0^floor(log10(bw))
        r = bw / e
        if r <= 2
            step = 2*e
        elseif r <= 5
            step = 5*e
        else
            step = 10*e
        end
    end
    if closed == :right
        start = step*(ceil(lo/step)-1)
        nm1 = iceil((hi - start)/step)
    else
        start = step*floor(lo/step)
        nm1 = ifloor((hi - start)/step)+1
    end
    start:step:(start + nm1*step)
end

function histrange{T<:Integer}(v::AbstractArray{T}, n::Integer, closed::Symbol)
    if length(v) == 0
        return 0:1:0
    end
    lo, hi = extrema(v)
    if hi == lo
        step = 1
    else
        bw = (hi - lo) / n
        e = int(10^max(0,ifloor(log10(bw))))
        r = bw / e
        if r <= 1
            step = e
        elseif r <= 2
            step = 2*e
        elseif r <= 5
            step = 5*e
        else
            step = 10*e
        end
    end
    if closed == :right
        start = step*(iceil(lo/step)-1)
        nm1 = iceil((hi - start)/step)
    else
        start = step*ifloor(lo/step)
        nm1 = ifloor((hi - start)/step)+1
    end
    start:step:(start + nm1*step)
end

histrange{N}(vs::NTuple{N,AbstractVector},nbins::NTuple{N,Integer},closed::Symbol) = map((v,n) -> histrange(v,n,closed),vs,nbins)
histrange{N}(vs::NTuple{N,AbstractVector},nbins::Integer,closed::Symbol) = map(v -> histrange(v,nbins,closed),vs)

## midpoints of intervals
midpoints(r::Range) = r[1:length(r)-1] + 0.5*step(r)
midpoints(v::AbstractVector) = [0.5*(v[i] + v[i+1]) for i in 1:length(v)-1]

## histograms ##
function sturges(n)  # Sturges' formula
    n==0 && return one(n)
    iceil(log2(n))+1
end

# N-dimensional histogram object
type Histogram{T<:Real,N,E}
    edges::E
    weights::Array{T,N}
    closed::Symbol
    function Histogram(edges::NTuple{N,AbstractArray},weights::Array{T,N},closed::Symbol)
        closed == :right || closed == :left || error("closed must :left or :right")
        map(x -> length(x)-1,edges) == size(weights) || error("Histogram edge vectors must be 1 longer than corresponding weight dimensions")
        new(edges,weights,closed)
    end
end
Histogram{T,N}(edges::NTuple{N,AbstractVector},weights::AbstractArray{T,N},closed::Symbol=:right) = Histogram{T,N,typeof(edges)}(edges,weights,closed)
Histogram{T,N}(edges::NTuple{N,AbstractVector},::Type{T},closed::Symbol=:right) = Histogram(edges,zeros(T,map(x -> length(x)-1,edges)...),closed)
Histogram{N}(edges::NTuple{N,AbstractVector},closed::Symbol=:right) = Histogram(edges,Int,closed)

function show(io::IO, h::Histogram)
    println(io, typeof(h))
    println(io,"edges:")
    for e in h.edges
        println(io,"  ",e)
    end
    println(io,"weights: ",h.weights)
    print(io,"closed: ",h.closed)
end

(==)(h1::Histogram,h2::Histogram) = (==)(h1.edges,h2.edges) && (==)(h1.weights,h2.weights) && (==)(h1.closed,h2.closed)

# 1-dimensional
Histogram{T}(edge::AbstractVector,weights::AbstractVector{T},closed::Symbol=:right) = Histogram{T,1,(typeof(edge),)}((edge,),weights,closed)
Histogram{T}(edge::AbstractVector,::Type{T},closed::Symbol=:right) = Histogram(edge,zeros(T,length(edge)-1),closed)
Histogram(edge::AbstractVector,closed::Symbol=:right) = Histogram(edge,Int,closed)

function push!{T,E}(h::Histogram{T,1,E}, x::Real,w::Real)
    i = if h.closed == :right
        searchsortedfirst(h.edges[1], x) - 1
    else
        searchsortedlast(h.edges[1], x)
    end
    if 1 <= i <= length(h.weights)
        @inbounds h.weights[i] += w
    end
    h
end
push!{T,E}(h::Histogram{T,1,E}, x::Real) = push!(h,x,one(T))

function append!{T}(h::Histogram{T,1}, v::AbstractVector)
    for x in v
        push!(h,x)
    end
    h
end
function append!{T}(h::Histogram{T,1}, v::AbstractVector,wv::WeightVec)
    for (x,w) in zip(v,wv.values)
        push!(h,x,w)
    end
    h
end

fit(::Type{Histogram},v::AbstractVector, edg::AbstractVector; closed::Symbol=:right) =
    append!(Histogram(edg,closed), v)
fit(::Type{Histogram},v::AbstractVector; closed::Symbol=:right, nbins=sturges(length(v))) =
    fit(Histogram, v, histrange(v,nbins,closed); closed=closed)

fit{W}(::Type{Histogram},v::AbstractVector, wv::WeightVec{W}, edg::AbstractVector; closed::Symbol=:right) =
    append!(Histogram(edg,W,closed), v, wv)
fit(::Type{Histogram},v::AbstractVector, wv::WeightVec; closed::Symbol=:right, nbins=sturges(length(v))) =
    fit(Histogram, v, wv, histrange(v,nbins,closed); closed=closed)

# N-dimensional
function push!{T,N}(h::Histogram{T,N},xs::NTuple{N,Real},w::Real)
    is = if h.closed == :right
        map((edge, x) -> searchsortedfirst(edge,x) - 1, h.edges, xs)
    else
        map(searchsortedlast, h.edges, xs)
    end
    try
        h.weights[is...] += w
    catch e
        isa(e,BoundsError) || rethrow(e)
    end
    h
end
push!{T,N}(h::Histogram{T,N},xs::NTuple{N,Real}) = push!(h,xs,one(T))

function append!{T,N}(h::Histogram{T,N}, vs::NTuple{N,AbstractVector})
    for xs in zip(vs...)
        push!(h,xs)
    end
    h
end
function append!{T,N}(h::Histogram{T,N}, vs::NTuple{N,AbstractVector},wv::WeightVec)
    for (xs,w) in zip(zip(vs...),wv.values)
        push!(h,xs,w)
    end
    h
end

fit{N}(::Type{Histogram}, vs::NTuple{N,AbstractVector}, edges::NTuple{N,AbstractVector}; closed::Symbol=:right) =
    append!(Histogram(edges,closed), vs)
fit{N}(::Type{Histogram}, vs::NTuple{N,AbstractVector}; closed::Symbol=:right, nbins=sturges(length(vs[1]))) =
    fit(Histogram, vs, histrange(vs,nbins,closed); closed=closed)

fit{N,W}(::Type{Histogram}, vs::NTuple{N,AbstractVector}, wv::WeightVec{W}, edges::NTuple{N,AbstractVector}; closed::Symbol=:right) =
    append!(Histogram(edges,W,closed), vs, wv)
fit{N}(::Type{Histogram},vs::NTuple{N,AbstractVector}, wv::WeightVec; closed::Symbol=:right, nbins=sturges(length(vs[1]))) =
    fit(Histogram, vs, wv, histrange(vs,nbins,closed); closed=closed)


text:
using VizCLIMA, Documenter

pages = Any[
    "Home" => "index.md",
    "VizCLIMA Scripts" => "vizclima_script_examples.md",
    "End-to-end Modeling (SLURM)" => "slurm_users_end_to_end.md",
    "Nimbus" => "datavis.md",
    "Gallery" => "gallery.md",
    "Third party software" => "third_party.md"
]

makedocs(
    sitename = "VizCLIMA.jl",
    doctest = false,
    strict = true,
    format = Documenter.HTML(
          prettyurls = get(ENV, "CI", nothing) == "true",
          mathengine = MathJax(Dict(
              :TeX => Dict(
                  :equationNumbers => Dict(:autoNumber => "AMS"),
                  :Macros => Dict()
              )
          ))
    ),
    clean = false,
    modules = [VizCLIMA],
    pages = pages,
)

deploydocs(
    repo = "github.com/CliMA/VizCLIMA.jl.git",
    target = "build",
    push_preview = true,
    forcepush = true,
)


text:
@testset "mixture sampler" begin

	cc = dummyChoiceContext()

	@testset "getparams, numparams, paramranges, sample" begin

		subsA = DataGenerators.BernoulliSampler()
		subsB = DataGenerators.DiscreteUniformSampler([7.0,9.0])
		s = DataGenerators.MixtureSampler(subsA,subsB)

	    @testset "numparams and paramranges" begin
			@test DataGenerators.numparams(s) == 2 + DataGenerators.numparams(subsA) + DataGenerators.numparams(subsB)
			prs = DataGenerators.paramranges(s)
	        @test typeof(prs) <: Vector{Tuple{Float64,Float64}}
			@test prs == [(0.0,1.0); (0.0,1.0); DataGenerators.paramranges(subsA); DataGenerators.paramranges(subsB)]
	    end

	    @testset "default params" begin
			@test DataGenerators.getparams(s) == [0.5; 0.5; DataGenerators.getparams(subsA); DataGenerators.getparams(subsB)]
	    end

	    @mtestset "default sampling" reps=Main.REPS alpha=Main.ALPHA begin
	        x, trace = DataGenerators.sample(s, (0,1), cc)
	        @test typeof(x) <: Int
 			@mtest_values_are [0,1,7,8,9] x
			if x <= 1
				@mtest_distributed_as Bernoulli(0.5) x
			else
				@mtest_distributed_as DiscreteUniform(7,9) x
			end
	    end

	end

	@testset "setparams" begin

		subsA = DataGenerators.BernoulliSampler()
		subsB = DataGenerators.DiscreteUniformSampler()
		subsC = DataGenerators.GeometricSampler()
		s = DataGenerators.MixtureSampler(subsA,subsB,subsC)
		prs = DataGenerators.paramranges(s)
		midparams = map(pr->robustmidpoint(pr[1],pr[2]), prs)

		@testset "parameters of subsamplers" begin

			params = [0.1, 0.2, 0.7, 0.4, 1.0, 6.0, 0.3]
			DataGenerators.setparams!(s, params)
			@test DataGenerators.getparams(subsA) == [0.4]
			@test DataGenerators.getparams(subsB) == [1.0, 6.0]
			@test DataGenerators.getparams(subsC) == [0.3]

			@mtestset "consistent with mixture" reps=Main.REPS alpha=Main.ALPHA begin
            	x, trace = DataGenerators.sample(s, (0,1), cc)
		        @mtest_distributed_as MixtureModel([Bernoulli(0.4),DiscreteUniform(1,6),Geometric(0.3)],[0.1,0.2,0.7]) x
			end

		end

		@testset "boundary values choice params $choiceparams" for choiceparams in [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0],]

            DataGenerators.setparams!(s, [choiceparams; [0.4, 1.0, 6.0, 0.3]])

			@mtestset "consistent with mixture" reps=Main.REPS alpha=Main.ALPHA begin
            	x, trace = DataGenerators.sample(s, (0,1), cc)
				if choiceparams[1] == 1.0
					@mtest_distributed_as Bernoulli(0.4) x
				elseif choiceparams[2] == 1.0
					@mtest_distributed_as DiscreteUniform(1,6) x
				elseif choiceparams[3] == 1.0
					@mtest_distributed_as Geometric(0.3) x
				end
			end

		end

		@testset "range check exception choice params $pidx bound $bidx" for pidx in 1:3, bidx in 1:2
			choiceparams = [0.5, 0.4, 0.1]
            choiceparams[pidx] = bidx == 1 ? prevfloat(prs[pidx][bidx]) : nextfloat(prs[pidx][bidx])
            @test_throws ErrorException DataGenerators.setparams!(s, [choiceparams, [0.4, 1.0, 6.0, 0.3]])
		end

	    @testset "wrong number of parameters" begin
	        @test_throws ErrorException DataGenerators.setparams!(s, midparams[1:end-1])
	        @test_throws ErrorException DataGenerators.setparams!(s, [midparams; 0.5])
	    end

	end

	@testset "estimateparams" begin

		subs1A = DataGenerators.BernoulliSampler()
		subs1B = DataGenerators.DiscreteUniformSampler()
		s1 = DataGenerators.MixtureSampler(subs1A,subs1B)

		@testset "from parameters $params" for params in [[0.6, 0.4, 0.7, 10.0, 13.0],]

			DataGenerators.setparams!(s1, params)

			subs2A = DataGenerators.BernoulliSampler()
			subs2B = DataGenerators.DiscreteUniformSampler()
			s2 = DataGenerators.MixtureSampler(subs2A,subs2B)
			otherparams = [0.3, 0.7, 0.3, -40.0, -27.0]
			DataGenerators.setparams!(s2, otherparams)

			traces = map(1:100) do i
				x, trace = DataGenerators.sample(s1, (0,1), cc)
				trace
			end
			estimateparams!(s2, traces)

			@mtestset "consistent with mixture" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as MixtureModel([Bernoulli(params[3]),DiscreteUniform(Int(params[4]),Int(params[5]))],params[1:2]) x
			end

		end

	    @testset "too few traces" begin

	        params = [0.2, 0.8, 0.5, 1.0, 6.0]
			DataGenerators.setparams!(s1, params)

			subs2A = DataGenerators.BernoulliSampler()
			subs2B = DataGenerators.DiscreteUniformSampler()
			s2 = DataGenerators.MixtureSampler(subs2A,subs2B)
			otherparams = [0.3, 0.7, 0.3, -40.0, -27.0]
			DataGenerators.setparams!(s2, otherparams)

	        traces = map(1:0) do i
	            x, trace = DataGenerators.sample(s1, (0,1), cc)
	            trace
	        end
	        estimateparams!(s2, traces)

			@mtestset "consistent with mixture" reps=Main.REPS alpha=Main.ALPHA begin
	        	x, trace = DataGenerators.sample(s2, (0,1), cc)
		        @mtest_distributed_as MixtureModel([Bernoulli(otherparams[3]),DiscreteUniform(Int(otherparams[4]),Int(otherparams[5]))],otherparams[1:2]) x
			end

	    end

	end

end


text:
using CSV, DataFrames, DecisionTree, Distributions, Gurobi, JuMP, LightGraphs, Parameters, SparseArrays, Statistics, ArgParse, Random

include("oracles/shortest_path_oracle.jl")
include("solver/util.jl")
include("solver/sgd.jl")
include("solver/reformulation.jl")
include("solver/random_forests_po.jl")
include("solver/validation_set.jl")
include("experiments/replication_functions.jl")

"""
Test to run one instance of the shortest path experiment.
Using the reformulation approach (spo+ formulation). (reformulation.jl)
replication_functions.jl : Functions for generating a single replication of a particular experiment
		(Dependencies: util.jl, sgd.jl, reformulation.jl, random_forests_po.jl, validation_set.jl)
		To change validation loss , will need to adjust: util.jl
		To change data generation method , will need to adjust: util.jl
		To change training loss algorithm (currently spo+), will need to adjust: util.jl, sgd.jl, validation_set.jl

Above dependencies from SPO repository:
https://github.com/paulgrigas/SmartPredictThenOptimize

"""

#environment
envOracle = setup_gurobi_env(method_type = :default, use_time_limit = false)
envReform = setup_gurobi_env(method_type = :method3, use_time_limit = false)

# parameter settings
p_features = 5
grid_dim = 5

num_lambda = 10
lambda_max = 100
lambda_min_ratio = 10.0^(-6)
holdout_percent = 0.25
n_holdout = round(Int, holdout_percent*n_train)

#For regularization used :lasso and :ridge

#setseed
rng_seed = 2128
Random.seed!(rng_seed)

#one instance of the problem:
n_train = 100
n_test = 1000
polykernel_degree = 2  #used 1, 2, 4, 8
polykernel_noise_half_width = 0 #and 0.5

sp_results = shortest_path_replication(grid_dim,
                    	n_train, n_holdout, n_test,
                        p_features, polykernel_degree, polykernel_noise_half_width;
                        num_lambda = 10, lambda_max = missing, lambda_min_ratio = 0.0001, regularization = :ridge,
                        gurobiEnvOracle = envOracle, gurobiEnvReform = envReform,
                        different_validation_losses = false,
                        include_rf = true)

max_stars_repo_name: katerinagiannoutsou/dso-coding-project

max_issues_repo_name: katerinagiannoutsou/dso-coding-project

max_forks_repo_name: katerinagiannoutsou/dso-coding-project


text:
__precompile__() # this module is safe to precompile
"""
   DifferentialExtinction
   Original Python code (solar_diffex.py) by Andrea Lin
   Ported to Julia by Eric Ford
   Based on Collier-Cameron et al. (2019) MNRAS 487, 1082.  doi:10.1093/mnras/stz1215
   See also Astronomical Algorithms (Meeus, 1991) p151 & p178
"""
module Horizons

using PyCall
using ..SunAsAStar

using Dates
using DataFrames, CSV

# Import Python package for querying JPL Horizons database
const JplHorizons = PyNULL()

function __init__()
	copy!(JplHorizons , pyimport("astroquery.jplhorizons") )
	@assert JplHorizons != PyNULL()
end

function query_horizons(;epochs::Union{Vector{Float64},Dict}, obs::Symbol, max_airmass::Float64 = 10.0, verbose::Bool = false)
	if typeof(epochs) <: Dict
		@assert haskey(epochs,"start")
		@assert haskey(epochs,"stop")
		@assert haskey(epochs,"step")
	elseif typeof(epochs) <: Vector{Float64}
		@assert 1 <= length(epochs) <= 10000
	end
	@assert 1 < max_airmass <= 10.0
	loc = get_obs_loc(obs)
	# Sometimes this fails.  Why?
	obj = JplHorizons.Horizons(id="10", location=loc, epochs=epochs, id_type="majorbody")
	if verbose println("# obj = ", obj) end
	#eph = obj.ephemerides(quantities="1,4,8,17,20,42", airmass_lessthan=max_airmass, refraction=true)  # WARNING: refraction=true stopped working! Why?
	eph = obj.ephemerides(quantities="1,4,8,17,20,42", airmass_lessthan=max_airmass )  # WARNING: refraction=true stopped working! Why?
	if verbose  println("# eph = ", eph)  end
	return eph
end

function cache_horizons_output_year(year::Int64, obs::Symbol; verbose::Bool = false)
	fn = "horizons_" * string(obs) * "_" * string(year) * ".csv"
	filename = joinpath(dirname(pathof(SunAsAStar)),"..","data","horizons",fn)
	if isfile(filename) && fielsize(filename) > 0
		return
	end
	df = DataFrame(:jd=>Float64[], :sol_npole_ang=>Float64[], :sol_dist_au=>Float64[])
	epoch_start = string(year) * "-01-01 00:00:00"
	epoch_stop  = string(year+1) * "-01-01 00:00:00"
	step = "6h"
	epochs = Dict("start"=>epoch_start, "stop"=>epoch_stop, "step"=>step)
	eph = query_horizons(epochs=epochs, obs=obs, max_airmass=10.0, verbose=verbose)
	df = DataFrame(:jd=>eph.columns["datetime_jd"], :sol_npole_ang=>eph.columns["NPole_ang"], :sol_dist_au=>eph.columns["delta"])
	CSV.write(filename,df)
	return df
end

function read_horizons_output(jd::Union{Float64,Vector{Float64}}, obs::Symbol; verbose::Bool = false)
	#@assert obs ∈ valid_obs
	years = unique(year.(Dates.julian2datetime.(jd)))
	@assert length(years) >= 1
	cache = Dict{Int64,Any}()
	for year in years
		cache[year] = read_horizons_output_year(year, obs=obs, verbose=verbose)
		#=
		fn = "horizons_" * string(obs) * "_" * string(year) * ".csv"
		filename = joinpath(dirname(pathof(SunAsAStar)),"..","data","horizons",fn)
		if isfile(filename)
			cache[year] = CSV.read(filename, DataFrame)
		else
			cache[year] = cache_horizons_output_year(year, obs)
		end
		cache[year].sol_npole_ang[cache[year].sol_npole_ang .> 180] .-= 360
		=#
	end
	if length(years) == 1
		return cache[first(keys(cache))]
	else
		sorted_keys = sort(collect(keys(cache)))
		output = cache[first(sorted_keys)]
		for k in sorted_keys[2:end]
			append!(output,cache[k])
		end
		return output
	end
end

function read_horizons_output_year(year::Integer ; obs::Symbol, verbose::Bool = false)
	fn = "horizons_" * string(obs) * "_" * string(year) * ".csv"
	filename = joinpath(dirname(pathof(SunAsAStar)),"..","data","horizons",fn)
	if isfile(filename) && filesize(filename) > 0
		cache = CSV.read(filename, DataFrame)
	else
		cache = cache_horizons_output_year(year, obs, verbose=verbose)
	end
	cache.sol_npole_ang[cache.sol_npole_ang .> 180] .-= 360
	return cache
	#first(keys(cache))
end

function get_horizons_output(jd::Float64; obs::Symbol, verbose::Bool = false)
  df = read_horizons_output(jd, obs, verbose=verbose)
  idxlo = searchsortedlast(df.jd, jd)
  @assert 1 <= idxlo < size(df,1)
  idxhi = idxlo + 1
  result = Dict{Symbol,Float64}()
  num = jd-df.jd[idxlo]
  denom = df.jd[idxhi]-df.jd[idxlo]
  whi = num/denom
  wlo = 1-whi
  for k in Symbol.(names(df))
	if k == :jd continue end
  	result[k] = df[idxlo,k] * wlo +  df[idxhi,k] * whi
  end
  result[:jd] = jd
  return result
end

#function get_horizons_output(jd::Vector{Float64}; obs::Symbol)
#  get_horizons_output.(jd,obs)
  #=
  df = read_horizons_output(jd, obs)
  idxlo = map(t->searchsortedlast(df.jd, t),jd)
  mask = idxlo .>= size(df,1)
  if length(mask) >= 1
  	println("# Maximum Δt = ", maximum(jd[mask] .- df.jd[end]) )
  	idxlo[mask] .= size(df,1)-1
end
  @assert all(1 .<= idxlo .<= size(df,1))
  result = DataFrame(:jd=>jd)
  whi = (jd.-df.jd[idxlo]) ./ (df.jd[idxlo.+1].-df.jd[idxlo])
  #wlo = 1.-whi
  for k in  Symbol.(names(df))
	if k == :jd continue end
  	result[!,k] = df[idxlo,k] .* (1.0.-whi) .+  df[idxlo.+1,k] .* whi
  end
  return result
  =#
#end

export query_horizons, get_horizons_output

end # module Horizons


text:
# This file calculates the solution to Laplace's equation via the adaptive spectral method.
# Δu = 0,
# u|Γ = u_0,
# u^i = log|z-z_s|,
# u = u^i + u^s.
# The normal derivative ∂u/∂n of the solution is calculated on the constant-charge boundary Γ.
# The reflected solution is calculated by convolving ∂u/∂n with the fundamental solution.
# Then, the total solution is obtained by summing the incident and the reflected solutions.

using ApproxFun, SingularIntegralEquations

z_s = 2.0
ui = (x,y) -> logabs(complex(x,y)-z_s)
g1 = (x,y) -> 1/2

# Set the domains.
N = 10
r = 1e-1
cr = exp.(im*2π*(0:N-1)/N)
crl = (1-2im*r)cr
crr = (1+2im*r)cr
dom = ∪(Segment.(crl,crr)) # All infinitesimal plates
#dom = ∪(Circle,cr,ones(length(cr))r) # All wires
#dom = ∪(Segment.(crl[1:2:end],crr[1:2:end])) ∪ ∪(Circle,cr[2:2:end],ones(length(cr[2:2:end]))r) # Interlaced wires and plates

sp = Space(dom)
cwsp = CauchyWeight(sp⊗sp,0)
uiΓ,⨍ = Fun(t->ui(real(t),imag(t))+0im,sp),DefiniteLineIntegral(dom)

@time G = GreensFun(g1,cwsp;method=:Cholesky)

@time φ0,∂u∂n=[0 ⨍;1 ⨍[G]]\[0.;uiΓ]

println("The length of ∂u∂n is: ",ncoefficients(∂u∂n))

us = (x,y) -> -logkernel(∂u∂n,complex(x,y))/2
ut = (x,y) -> ui(x,y) + us(x,y)
println("This is the approximate gradient: ",((ut(1e-5,0.)-ut(-1e-5,0.))/2e-5))


text:
module ExponentialUtilities
using LinearAlgebra, SparseArrays
using LinearAlgebra: exp!

include("phi.jl")
include("arnoldi.jl")
include("krylov_phiv.jl")
include("krylov_phiv_adaptive.jl")
include("krylov_expv.jl")

export phi, phi!, KrylovSubspace, arnoldi, arnoldi!, lanczos!, ExpvCache, PhivCache,
    expv, expv!, phiv, phiv!, expv_timestep, expv_timestep!, phiv_timestep, phiv_timestep!,
    StegrCache, get_subspace_cache
end


text:
using GeneticAlgorithms
include("equalityga.jl")

model = runga(equalityga; initial_pop_size = 16)

population(model)