text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

# According to file format specified in Matlab script
# written by Mykel Kochenderfer, mykel@stanford.edu
#
# This script converts a log file produced by RLES and converts it to a "scripts" file
# of the following format:
#
#   SCRIPTS FILE:
#   The scripts file contains a set of scripted encounters. Each
#   encounter is defined by a set of encounter scripts associated with a
#   fixed number of aircraft. The file is organized as follows:
#
#   [Header]
#   uint32 (number of encounters)
#   uint32 (number of aircraft)
#       [Encounter 1]
#           [Initial]
#               [Aircraft 1]
#               double (initial airspeed in ft/s)
#               double (initial north position in ft)
#               double (initial east position in ft)
#               double (initial altitude in ft)
#               double (initial heading angle in radians)
#               double (initial flight path angle in radians)
#               double (initial roll angle in radians)
#               double (initial airspeed acceleration in ft/s^2)
#               ...
#               [Aircraft n]
#               double (initial airspeed in ft/s)
#               double (initial north position in ft)
#               double (initial east position in ft)
#               double (initial altitude in ft)
#               double (initial heading angle in radians)
#               double (initial flight path angle in radians)
#               double (initial roll angle in radians)
#               double (initial airspeed acceleration in ft/s^2)
#           [Updates]
#               [Aircraft 1]
#               uint8 (number of updates)
#                   [Update 1]
#                   double (time in s)
#                   double (commanded vertical rate in ft/s)
#                   double (commanded turn rate in rad/s)
#                   double (commanded airspeed acceleration in ft/s^2)
#                   ...
#                   [Update m]
#                   double (time in s)
#                   double (commanded vertical rate in ft/s)
#                   double (commanded turn rate in rad/s)
#                   double (commanded airspeed acceleration in ft/s^2)
#               ...
#               [Aircraft n]
#                   ...
#       ...
#       [Encounter k]
#           ...

module Log_To_Scripted

export log_to_scripted

using ..DefineSave
using ..SaveHelpers

include("corr_aem_save_scripts.jl")

using JSON
using DataFrames

function log_to_scripted{T<:AbstractString}(filenames::Vector{T};
    outfile::AbstractString = "scripted.dat")

    d = trajLoad(filenames[1]) #use the first one as a reference
    num_aircraft = get_num_aircraft(d)
    num_encounters = length(filenames) #one encounter per log file
    encounters = Array(Dict{String, Array{Float64, 2}}, num_aircraft, num_encounters)

    #encounter i
    for (i, file) in enumerate(filenames)
        d = trajLoad(file)

        #make sure all of them have the same number of aircraft
        @assert num_aircraft == get_num_aircraft(d)

        #aircraft j
        for j = 1:num_aircraft
            encounters[j, i] = Dict{String,Array{Float64, 2}}()
            encounters[j, i]["initial"] = j2s_initial(d, j)
            encounters[j, i]["update"] = j2s_update(d, j)'
        end
    end
    save_scripts(outfile, encounters, numupdatetype=UInt8)
    encounters
end

function j2s_initial(d::TrajLog, aircraft_number::Int64)
    #d is the loaded encounter
    #j is aircraft number

    out = Array(Float64, 1, 8)
    initial = get_log(d, "Initial", aircraft_number)
    airspeed            = initial[1, :v]
    pos_north           = initial[1, :y]
    pos_east            = initial[1, :x]
    altitude            = initial[1, :z]
    heading             = initial[1, :psi]   |> deg2rad #to radians
    flight_path_angle   = initial[1, :theta] |> deg2rad #to radians
    roll_angle          = initial[1, :phi]   |> deg2rad #to radians
    airspeed_accel      = initial[1, :v_d]

    out[1, :] = [airspeed, pos_north, pos_east, altitude, heading, flight_path_angle,
              roll_angle, airspeed_accel]
    out
end

function j2s_update(d::TrajLog, aircraft_number::Int64)
    #d is the loaded encounter,
    pr = get_log(d, "Response", aircraft_number)
    t_end = nrow(pr)
    out = Array(Float64, t_end-1, 4)
    for t = 2:t_end #ignore init values
        t_     = pr[t, :t]
        h_d    = pr[t, :h_d]
        psi_d  = pr[t, :psi_d] |> deg2rad #to radians
        v_d    = pr[t, :v_d]
        out[t-1, :] = Float64[t_, h_d, psi_d, v_d]
    end
    out
end

function log_to_scripted(filename::AbstractString)
  log_to_scripted([filename], outfile = string(getSaveFileRoot(filename), "_scripted.dat"))
end

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
import Glob
println("Loading:")
try
    import InsarLOS
catch
    println("insar los import error")
end
try
    import Sario
catch
    println("Sario import error")
end

if ispath("../extra_files/")
    dbfile = Glob.glob("../extra_files/*.db*")[1]
else
    dbfile = Glob.glob("../*.db*")[1]
end

@time InsarLOS.create_los_map(directory = ".", dbfile = dbfile, outfile = "los_map.h5");

ext: jl


text:
"""
# Cooley–Tukey FFT algorithm

## References

- [clrs]: Cormen, Thomas H., Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein.
  Introduction to algorithms. MIT press, 2009.
- [manber]: Manber, Udi. Introduction to algorithms: a creative approach. Addison-Wesley Longman
  Publishing Co., Inc., 1989.
- [wiki]: https://en.wikipedia.org/wiki/Cooley–Tukey_FFT_algorithm
"""
module FastFourierTransform

using OffsetArrays

export
    Direction, Forward, Backward,
    dft, idft,
    fft_recursive_naive, ifft_recursive_naive,
    fft_recursive, ifft_recursive,
    fft_iterative, ifft_iterative,
    eval_poly, poly_mul

#=

Primitive nth root of unity (can be either clockwise or counterclockwise):

    w_n = exp(-im * 2π / n)

Cancellation lemma ([clrs] Lemma 30.3):

    w_{n d}^{d k} = exp(-im * 2π * (k d)/(n d)) = w_n^k

Halving lemma ([clrs] Lemma 30.5) (by Cancellation lemma we have $(w_n^k)^2 = w_{n/2}^k$):

      { (w_n^k)^2 for k in 0..(n - 1) }
    = { w_{n/2}^k for k in 0..(n/2 - 1) } ∪ { w_{n/2}^k for k in (n/2)..(n - 1) }
    = { w_{n/2}^k for k in 0..(n/2 - 1) }

Summation lemma ([clrs] Lemma 30.6):

      \sum_{j = 0}^{n - 1} w_n^{j k}
    = ((w_n^k)^n - 1) / (w_n^k - 1)
    = (1^k - 1) / (w_n^k - 1)
    = 0

Danielson–Lanczos(?) lemma ([manber] Eq. 9.8):

    P(x) = P_e(x^2) + x * P_o(x^2)

So the problem of computing (note it's a list instead of a set)

    [ P(w_n^i) for i in 0..(n - 1) ]

reduces to

    [ P_e((w_n^i)^2) + w_n^i * P_o((w_n^i)^2) for i in 0..(n - 1) ]

which by Halving lemma further reduces to (also uses [clrs] Corollary 30.4)

      [ P_e(w_{n/2}^i) + w_n^i * P_o(w_{n/2}^i) for i in 0..(n - 1) ]
    = [ P_e(w_{n/2}^i) + w_n^i * P_o(w_{n/2}^i) for i in 0..(n/2 - 1) ] ++
      [ P_e(w_{n/2}^i) - w_n^i * P_o(w_{n/2}^i) for i in 0..(n/2 - 1) ]

The inverse of DFT matrix `V^(-1)`:

    V_ij(w) = w^{i j}
    V^{-1}_ij(w) = 1/n V_ij(1/w) = 1/n w^{-i j}
    (V V^{-1})_ij(w) = 1/n w_{i k} w^{-k j}  ...(*)

    (*) = 1/n n = 1,  if i = j
    (*) = 1/n n = 1/n \sum_{k = 0}^{n - 1} (k d) = 0,  if i - j = d ≠ 0  (by Summation lemma)

=#

# All `OffsetArray`s are 0-based.

@inline function zero_based(A::AbstractArray)::OffsetArray
    (r, c) = size(A) .- 1
    OffsetArray(A, 0:r, 0:c)
end

@inline function zero_based(v::AbstractVector)::OffsetVector
    OffsetVector(v, 0:(length(v) - 1))
end

function divide(W::OffsetArray{T})::Tuple{OffsetArray{T}, OffsetArray{T}} where T
    (n, nn) = size(W)
    @assert n == nn && n > 1 && iseven(n) "Invalid input"

    W_e = @view parent(W)[1:(n ÷ 2), 1:2:end]
    W_o = @view parent(W)[1:(n ÷ 2), 2:2:end]
    (W_e, W_o) .|> zero_based
end

function divide_vec(
    x::OffsetVector{T}
)::Tuple{OffsetVector{T}, OffsetVector{T}} where T
    n = length(x)
    # `n` can be 1 -- an edge case
    @assert (n == 1) || (n > 1 && iseven(n)) "Invalid input"

    x_e = @view parent(x)[1:2:end]
    x_o = @view parent(x)[2:2:end]
    (x_e, x_o) .|> zero_based
end

function bit_reverse(x::Int, bit_len::Int)::Int
    y = 0
    for i in 1:bit_len
        y <<= 1
        y |= x & 1
        x >>= 1
    end
    y
end

function bit_reverse_copy(a::OffsetVector{T})::OffsetVector{T} where T
    n = length_checked(a)
    bit_len = trailing_zeros(n)
    b = zero(a)
    for i in 0:(n - 1)
        b[bit_reverse(i, bit_len)] = a[i]
    end
    b
end

@enum Direction Forward Backward

function dft_matrix(T::Type, n::Int, dir::Direction)::OffsetArray{Complex{T}}
    W = OffsetArray(Array{Complex{T}, 2}(undef, n, n), 0:(n - 1), 0:(n - 1))

    # Sets `W_ij = exp(-im * 2π * i*j/n)` (forward). It is already O(n^2).
    for i in 0:(n - 1)
        for j in i:(n - 1)
            W[i, j] = W[j, i] =
                if dir == Forward
                    exp((-2π)im * (i * j / n))
                else
                    # (!) Do not put `/ n` here.
                    exp((2π)im * (i * j / n))
                end
        end
    end

    W
end

function dft_vector(T::Type, n::Int, dir::Direction)::OffsetVector{Complex{T}}
    # Sets `w_i = exp(-im * 2π * i/n)` (forward).
    #
    # Type annotation is required for `n == 1`.
    w::Vector{Complex{T}} =
        if dir == Forward
            [exp((-2 * T(π))im * i / n) for i in 0:(n ÷ 2 - 1)]
        else
            # (!) Do not put `/ n` here as unlike `W`, `V[:, 0]` is not 0
            [exp((2 * T(π))im * i / n) for i in 0:(n ÷ 2 - 1)]
        end
    OffsetVector(w, 0:(n ÷ 2 - 1))
end

function length_checked(a::AbstractVector)::Int
    n = length(a)
    @assert ispow2(n) "Coefficient length is not a power of 2"
    n
end

#===== DFT (by direct matrix-vector multiplication) =====#

# Reference implementation, for verification.
function dft(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    W = dft_matrix(T, n, Forward)
    parent(W) * a
end

function idft(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    W = dft_matrix(T, n, Backward)
    parent(W) * a / n
end

#===== FFT (by recursing on DFT matrix) =====#
#
# This section follows the explanation in [manber] 9.6.

function fft_recursive_naive(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    length_checked(a)
    fft_recursive_naive_impl(a, Forward)
end

function ifft_recursive_naive(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    fft_recursive_naive_impl(a, Backward) / n
end

function fft_recursive_naive_impl(
    a::Vector{Complex{T}},
    dir::Direction,
)::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    a = zero_based(a)

    W = dft_matrix(T, n, dir)

    function fft_recursive_naive_r!(
        W::OffsetArray{Complex{T}},
        a::OffsetVector{Complex{T}},
    )::Vector{Complex{T}} where T<:AbstractFloat
        local n = length(a)
        if n == 1
            return [W[0, 0] * a[0]]
        end

        (W_e, W_o) = divide(W)
        (a_e, a_o) = divide_vec(a)

        # Warning: Tricky code. Should also use `W_e` to compute `p_o`, and scales the result with
        # `w`. This is due to the odd sub matrix (`W_o`) no longer holds the recursing property.
        w = parent(W)[1:(n ÷ 2), 2]
        p_e = fft_recursive_naive_r!(W_e, a_e)
        p_o = fft_recursive_naive_r!(W_e, a_o) .* w

        [
            p_e + p_o
            p_e - p_o
        ]
    end

    fft_recursive_naive_r!(W, a)
end

#===== FFT (the "standard" recursive version) =====#

function fft_recursive(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    w = dft_vector(T, n, Forward)
    fft_recursive_r(zero_based(a), w)
end

function ifft_recursive(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    w = dft_vector(T, n, Backward)
    fft_recursive_r(zero_based(a), w) / n
end

function fft_recursive_r(
    a::OffsetVector{Complex{T}},
    w::OffsetVector{Complex{T}},
)::Vector{Complex{T}} where T<:AbstractFloat
    n = length(a)
    if n == 1
        return [a[0]]
    end

    (a_e, a_o) = divide_vec(a)
    (w_e, _) = divide_vec(w)  # relies on DFT matrix's special structure
    p_e = fft_recursive_r(a_e, w_e)
    p_o = fft_recursive_r(a_o, w_e)

    [
        p_e .+ parent(w) .* p_o
        p_e .- parent(w) .* p_o
    ]
end

#===== FFT (the iterative version) =====#

function fft_iterative(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    fft_iterative_impl(a, Forward)
end

function ifft_iterative(a::Vector{Complex{T}})::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    fft_iterative_impl(a, Backward) / n
end

function fft_iterative_impl(
    a::Vector{Complex{T}},
    dir::Direction,
)::Vector{Complex{T}} where T<:AbstractFloat
    n = length_checked(a)
    log2n = trailing_zeros(n)

    a = zero_based(a)
    a = bit_reverse_copy(a)
    w = dft_vector(T, n, dir)

    # The `m` here is actually `1/2 m` in [clrs].
    for s in 1:log2n
        m = 2 ^ (s - 1)
        k = 0
        while k < n
            # butterfly operation
            for i = 0:(m - 1)
                # Indexes `w_{2 ^ s}^i` into vector `w_n`, so its step is `2 ^ (log2n - s)`.
                t = w[(2 ^ (log2n - s)) * i] * a[k + i + m]
                u = a[k + i]
                a[k + i] = u + t
                a[k + i + m] = u - t
            end

            k += m * 2
        end
    end
    parent(a)
end

#===== Polynomial multiplication =====#

"""
Evaluates polynomial `p(x)` using Horner's rule, where `p` is the coefficient list for `x` from
`x^0` to `x^{n - 1}`.
"""
function eval_poly(p::Vector{T}, x::T)::T where T<:Real
    foldr((a::T, s::T) -> s * x + a, p)
end

function poly_mul(
    a::Vector{T}, b::Vector{T},
    fft_f::Function=fft_iterative, ifft_f::Function=ifft_iterative,
)::Vector{T} where T<:Real
    # zero-padding
    na = length(a)
    nb = length(b)
    # da = na - 1, db = nb - 1; d = da + db = na + nb - 2, n = d + 1 = na + nb - 1
    n = na + nb - 1
    np = nextpow(2, n)
    a = vcat(complex(a), zeros(Complex{T}, np - na))
    b = vcat(complex(b), zeros(Complex{T}, np - nb))

    # By Convolution theorem ([clrs] Theorem 30.8)
    fa = fft_f(a)
    fb = fft_f(b)
    fp = fa .* fb
    real(ifft_f(fp)[1:n])
end

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
# File for generating data for Figure 3

using Pkg
using Random, Distributions
using Statistics
using ProgressMeter
using DataFrames

include("/home/christian/UniMA/EMMA/src/RFR.jl")
include("/home/christian/UniMA/EMMA/src/cross_val.jl")
include("/home/christian/UniMA/EMMA/src/aux_functions.jl")

function get_best_model(n::Int, d::Int, type::String, σ::Int, m_features::Int)

    xtrain1, xtest1, ytrain1, ytest1 = make_data(n, d, type, σ)

    a_list = collect(LinRange(0, 50, 15))

    d1 = Dict{Symbol, Vector{Float64}}(
        :max_features => [m_features],
        :n_trees => [30],
        :α => a_list)

    cv = cross_val(d1, random_state = 0)
    fit!(cv, xtrain1, ytrain1, nfolds=3)

    return(best_model(cv, "mse"))
end

function get_data(n_runs, max_d, σ, opt_α=false)

    # Allocate space for results
    res_mat = zeros(max_d, n_runs)

    # Setting parameters
    n = 4000
    d = 10
    m_features = d

    if opt_α
        # Getting the best model given error term variance
        rf_opt = get_best_model(n, d, "friedman", σ, m_features)
        println(rf_opt.α)
    end

    @showprogress for r in 1:n_runs

        x_train, x_test, y_train, y_test = make_data(n, d, "friedman", σ)

        if opt_α
            # Fitting the forest with optimally chosen α with new data
            rf = rf_opt
        else
            # Else choose CART
            d1 = Dict{Symbol, Vector{Float64}}(
                :max_features => [m_features],
                :n_trees => [30],
                :α => [0.0])
            rf = RFR(param_dict = d1)
        end

        fit!(rf, x_train, y_train)

        result_arr = Array{Float64}(undef, 0, 2)

        for i in 1:length(rf.trees)

            # Save depth and splitting dimension
            res = zeros(length(rf.trees[i].depth_list), 2)
            res[:,1] = rf.trees[i].depth_list
            res[:,2] = rf.trees[i].split_dimensions

            result_arr = vcat(result_arr, res)
        end

        # Only comparing non-terminal nodes
        plot_data = result_arr[result_arr[:,2].!=-2, :]
        # Noisy dimensions are all regressors above 5
        plot_data[:,2] = plot_data[:,2].>5

        # Save figures as dataframe
        plot_data = DataFrame(plot_data, :auto)
        rename!(plot_data, ["depth", "noise"])

        # Grouping them by depth
        gdf = groupby(plot_data, :depth)
        # Getting mean of fraction split on noise
        cdf2 = combine(gdf, :noise => mean)

        res_mat[1:size(cdf2, 1), r] = cdf2[!, :noise_mean]
    end
    return res_mat
end

function data_to_plot(res_mat::Matrix, x_max::Int, σ)
    wide_df = DataFrame(res_mat[1:x_max, :], :auto)
    wide_df_plot = DataFrame()
    wide_df_plot[!, "depth"] = collect(1:x_max)
    wide_df_plot[!, "mean"] = mean.(eachrow(wide_df))

    # Getting minimum and maximum value of fraction in addition to mean
    wide_df_plot[!, "ymin"] = minimum.(eachrow(wide_df))
    wide_df_plot[!, "ymax"] = maximum.(eachrow(wide_df))
    wide_df_plot[!, "approach"] = repeat([String(σ)], size(wide_df_plot, 1))

    return wide_df_plot

end

# Set seed for reproducability
Random.seed!(68151)

n_runs=100
max_d=50 #Max depth allowed - need to set this parameter; max acieved is around 30

# Getting data
res_mat1 = get_data(n_runs, max_d, 1, false)
res_mat3 = get_data(n_runs, max_d, 3, false)
res_mat8 = get_data(n_runs, max_d, 8, false)

# Figure plots until depth of 20
x_max = 20
r1 = data_to_plot(res_mat1, x_max, "1")
r2 = data_to_plot(res_mat3, x_max, "3")
r3 = data_to_plot(res_mat8, x_max, "8")

plot_df = vcat(r1, r2, r3)
jldsave("data/noise_split_dim.jld2"; plot_df)

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
# SPDX-License-Identifier: MIT

#######################################################################
# Model definition

"""
    struct NormalDiag{D} <: AbstractModel
        μ
        λ
    end

Normal distribution with diagonal covariance matrix.
"""
struct NormalDiag{D} <: AbstractModel
    μ::T where T<:BayesParam
    λ::BayesParamList{D,T} where T<:BayesParam
end

function NormalDiag(;datadim, pstrength = 1, init_noise_std = 0, m₀ = nothing,
                    λ₀ = nothing)
    D = datadim
    m₀ = isnothing(m₀) ? zeros(D) : m₀
    λ₀ = isnothing(λ₀) ? ones(D) : λ₀

    μprior = ExpFamilyDistributions.NormalDiag(m₀, λ₀/pstrength)
    μposterior = ExpFamilyDistributions.NormalDiag(m₀ + randn(D)*init_noise_std, λ₀/pstrength)
    μ = BayesParam(μprior, μposterior)

    α₀ = pstrength
    β₀ = pstrength ./ λ₀
    λ = tuple([BayesParam(Gamma{Float64}(α₀, β₀ᵢ), Gamma{Float64}(α₀, β₀ᵢ)) for β₀ᵢ in β₀]...)

    NormalDiag{D}(μ, λ)
end

#######################################################################
# Model interface

basemeasure(::NormalDiag, x::AbstractVector{<:Real}) = -.5*length(x)*log(2π)

function vectorize(m::NormalDiag)
    Tλ = gradlognorm.(getproperty.(m.λ, :posterior), vectorize = false)
    λ = [Tλᵢ[1] for Tλᵢ in Tλ]
    lnλ = [Tλᵢ[2] for Tλᵢ in Tλ]
    μ, μ² = gradlognorm(m.μ.posterior, vectorize = false)
    vcat(λ .* μ, -.5 .* λ, -.5 * (λ' * μ² - sum(lnλ)))
end

statistics(m::NormalDiag, x::AbstractVector{<:Real}) = vcat(x, x.^2, 1)

function loglikelihood(m::NormalDiag, x::AbstractVector{<:Real})
    D = @Zygote.ignore length(x)
    Tη = vectorize(m)
    Tx = @Zygote.ignore statistics(m, x)
    Tη'*Tx + basemeasure(m, x)
end

function loglikelihood(m::NormalDiag, X::AbstractVector{<:AbstractVector})
    D = @Zygote.ignore length(X[1])
    Tη = vectorize(m)
    TX = @Zygote.ignore statistics.([m], X)
    dot.([Tη], TX) .+ basemeasure.([m], X)
end

function getparam_stats(m::NormalDiag, X, resps)
    λ_s = λstats(m, X, resps)
    μ_s = μstats(m, X, resps)
    T = eltype(μ_s)
    retval = Dict{BayesParam, Vector{T}}(m.μ => μ_s)
    for (λᵢ,λᵢ_s) in zip(m.λ, λ_s)
        retval[λᵢ] = λᵢ_s
    end
    retval
end
getparam_stats(m::NormalDiag, X) = getparam_stats(m, X, ones(eltype(X[1]), length(X)))

#######################################################################
# λ statistics

function λstats(m::NormalDiag, X, resps)
    μ, μ² = gradlognorm(m.μ.posterior, vectorize = false)
    s1 = sum(t -> ((x,z) = t ; z * (-.5 * x.^2 + μ .* x -.5*μ²)) , zip(X, resps))
    s2 = .5*ones(eltype(μ), length(μ))*sum(resps)
    map(x -> eltype(μ)[x...], zip(s1, s2))
end

#######################################################################
# μ statistics

function μstats(m::NormalDiag, X, z)
    λ = [gradlognorm(λᵢ.posterior, vectorize = false)[1] for λᵢ in m.λ]
    s = sum(x -> vcat(λ .* x, -λ/2), X)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
@testset "778.swim-in-rising-water.jl" begin
    @test swim_in_water([[0,2],[1,3]]) == 3
    @test swim_in_water([[0,1,2,3,4],[24,23,22,21,5],[12,13,14,15,16],[11,17,18,19,20],[10,9,8,7,6]]) == 16
end


text:
@testset "Traits" begin
    A = rand(5,4,3)
    Aoneto = IndicesArray(A);

    r1 = range(.1, stop = .2, length=5)
    r2 = ["a", "b", "c", "d"]
    r3 = 2:4

    ind1 = asindex(r1)
    ind2 = asindex(r2)
    ind3 = asindex(r3)

    named_ind1 = OneToIndex{:a}(ind1)
    named_ind2 = OneToIndex{:b}(ind2)
    named_ind3 = OneToIndex{:c}(ind3)

    #map(asindex, (r1, r2, r3), axes(A))
    Aindices = IndicesArray(A, r1, r2, r3);
    Anamed = IndicesArray(Aindices, a = r1, b = r2, c = r3);

    # indexnames,
    # unname,
    @testset "filteraxes" begin
        @test filteraxes(i -> isequal(length(i), 4), A) == (Base.OneTo(4),)
        @test filteraxes(i -> isequal(length(i), 4), Aindices) == (ind2,)
        @test filteraxes(i -> isequal(length(i), 4), Anamed) == (OneToIndex{:b}(ind2),)
    end
    @testset "findaxes" begin
        @test findaxes(i -> isequal(length(i), 3), A) == (3,)
        @test findaxes(i -> isequal(length(i), 3), Aindices) == (3,)
        @test findaxes(i -> isequal(length(i), 3), Anamed) == (3,)
    end

    @testset "to_dims" begin
        @test to_dims(A, 2) == 2
        @test to_dims(Aindices, 2) == 2
        @test to_dims(Anamed, 2) == 2
        @test to_dims(Anamed, :a) == 1
        @test to_dims(Anamed, (1, 2)) == (1, 2)
        @test to_dims(Anamed, (:a, :b)) == (1, 2)
        @test to_dims(Anamed, :) == (1, 2, 3)
        @test to_dims(Aindices, 1) == 1
        @test to_dims(Aindices, :a) == 0
        @test to_dims(Anamed, :d) == 0
    end

    @testset "mapaxes" begin
        @test mapaxes(length, Aindices) == (length(named_ind1),
                                            length(named_ind2),
                                            length(named_ind3))
    end

    @testset "dropaxes" begin
        @test dropaxes(Anamed, dims=:b) == (named_ind1, named_ind3)
    end

    @testset "permuteaxes" begin
        @test permuteaxes(Anamed, (:c, :b, :a)) == (named_ind3, named_ind2, named_ind1)
    end

    @testset "reduceaxes" begin
        @test reduceaxes(A, dims=2) == (Base.OneTo(5), 1, Base.OneTo(3))
        @test reduceaxes(Aindices, dims=2) == (ind1, reduceaxis(ind2), ind3)
        @test reduceaxes(Anamed, 2) == (OneToIndex{:a}(ind1), reduceaxis(OneToIndex{:b}(ind2)), OneToIndex{:c}(ind3))
        @test reduceaxes(Anamed, :b) == (OneToIndex{:a}(ind1), reduceaxis(OneToIndex{:b}(ind2)), OneToIndex{:c}(ind3))
        @test reduceaxes(Anamed, (1, 2)) == (reduceaxis(OneToIndex{:a}(ind1)), reduceaxis(OneToIndex{:b}(ind2)), OneToIndex{:c}(ind3))
        @test reduceaxes(Anamed, (:a, :b)) == (reduceaxis(OneToIndex{:a}(ind1)), reduceaxis(OneToIndex{:b}(ind2)), OneToIndex{:c}(ind3))
    end

    # TODO named_axes
    @testset "named_axes" begin
        @test named_axes(Aindices) == (OneToIndex{:dim_1}(ind1),
                                      OneToIndex{:dim_2}(ind2),
                                      OneToIndex{:dim_3}(ind3))
        @test named_axes(Anamed) == (named_ind1, named_ind2, named_ind3)
        @test unnamed_axes(Anamed) == axes(Aindices)
        @test has_indexnames(Anamed) == true
        @test has_indexnames(Aindices) == false
        @test unname((a=1,b=2)) == (1, 2)
        @test indexnames((a=1, b=2)) == (:a, :b)
    end

    @testset "Range interface" begin
        a = [1, 2, 3]
        r = 1:3
        mr = StepMRange(1,1,4)

        @testset "first" begin
            @testset "can_grow_first" begin
                @test @inferred(can_grow_first(a)) == true
                @test @inferred(can_grow_first(r)) == false
                @test @inferred(can_grow_first(mr)) == true
            end

            @testset "set_first!" begin
                @test @inferred(set_first!(a, 2)) == [2, 2, 3]
                @test_throws MethodError set_first!(r, 2)
            end
        end

        @testset "last" begin
            @testset "can_grow_last" begin
                @test @inferred(can_grow_last(a)) == true
                @test @inferred(can_grow_last(r)) == false
                @test @inferred(can_grow_last(mr)) == true
            end

            @testset "set_last!" begin
                @test @inferred(set_last!(a, 2)) == [2, 2, 2]
                @test_throws MethodError set_last!(r, 2)
            end
        end

        @testset "step" begin
            @testset "has_step" begin
                @test @inferred(has_step(a)) == false
                @test @inferred(has_step(r)) == true
                @test @inferred(has_step(mr)) == true
            end

            @testset "can_set_step" begin
                @test @inferred(can_set_step(a)) == false
                @test @inferred(can_set_step(r)) == false
            end

            @testset "set_last!" begin
                @test_throws MethodError set_step!(r, 2)
                set_step!(mr, 2)
                @test @inferred(step(mr)) == 2
            end
        end
    end
end


text:
struct Up   end
struct Down end
struct Rand end

const UP  = Up()
const DWN = Down()
const RND = Rand()

macro defRound(op, eft)
    quote
        function $(esc(op))(::Up, a, b)
            (res, err) = $(esc(eft))(a, b)
            if err > 0
                $(esc(nextfloat))(res)
            else
                res
            end
        end

        function $(esc(op))(::Down, a, b)
            (res, err) = $(esc(eft))(a, b)
            if err < 0
                $(esc(prevfloat))(res)
            else
                res
            end
        end

        function $(esc(op))(::Rand, a, b)
            (res, err) = $(esc(eft))(a, b)
            if err < 0
                Determinism.randBool(a,b) ? $(esc(prevfloat))(res) : res
            elseif err > 0
                Determinism.randBool(a,b) ? $(esc(nextfloat))(res) : res
            else
                res
            end
        end
    end
end

import Base: +, -, *, /
@defRound(+, EFT.twoSum)
@defRound(-, EFT.twoDiff)
@defRound(*, EFT.twoProd)
@defRound(/, EFT.twoDiv)


text:
using InteractiveUtils
using Pope

@testset "wait_for_file_to_exist" begin
  endchannel = Channel{Bool}(1)
  fname = tempname()
  @async begin sleep(1); touch(fname) end
  t1=@elapsed @test Pope.wait_for_file_to_exist(fname, endchannel)
  @test t1>0.5
  @async begin sleep(1); put!(endchannel, true) end
  t2=@elapsed @test !Pope.wait_for_file_to_exist(tempname(), endchannel)
  @test t2>0.5
  #check that it returns true immediatley if file exists even if endchannel is already ready
  t3 = @elapsed @test Pope.wait_for_file_to_exist(fname, endchannel)
  @test t3<0.1
end


text:
# Contiguous on-device arrays

export VEDeviceArray, VEDeviceVector, VEDeviceMatrix

## construction

# NOTE: we can't support the typical `tuple or series of integer` style construction,
#       because we're currently requiring a trailing pointer argument.

struct VEDeviceArray{T,N,A} <: AbstractArray{T,N}
    shape::Dims{N}
    ptr::LLVMPtr{T,A}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    VEDeviceArray{T,N,A}(shape::Dims{N}, ptr::LLVMPtr{T,A}) where {T,A,N} = new(shape,ptr)
end

const VEDeviceVector = VEDeviceArray{T,1,A} where {T,A}
const VEDeviceMatrix = VEDeviceArray{T,2,A} where {T,A}

# outer constructors, non-parameterized
VEDeviceArray(dims::NTuple{N,<:Integer}, p::LLVMPtr{T,A})    where {T,A,N} = VEDeviceArray{T,N,A}(dims, p)
VEDeviceArray(len::Integer,              p::LLVMPtr{T,A})    where {T,A}   = VEDeviceVector{T,A}((len,), p)

# outer constructors, partially parameterized
VEDeviceArray{T}(dims::NTuple{N,<:Integer},   p::LLVMPtr{T,A}) where {T,A,N} = VEDeviceArray{T,N,A}(dims, p)
VEDeviceArray{T}(len::Integer,                p::LLVMPtr{T,A}) where {T,A}   = VEDeviceVector{T,A}((len,), p)
VEDeviceArray{T,N}(dims::NTuple{N,<:Integer}, p::LLVMPtr{T,A}) where {T,A,N} = VEDeviceArray{T,N,A}(dims, p)
VEDeviceVector{T}(len::Integer,               p::LLVMPtr{T,A}) where {T,A}   = VEDeviceVector{T,A}((len,), p)

# outer constructors, fully parameterized
VEDeviceArray{T,N,A}(dims::NTuple{N,<:Integer}, p::LLVMPtr{T,A}) where {T,A,N} = VEDeviceArray{T,N,A}(Int.(dims), p)
VEDeviceVector{T,A}(len::Integer,               p::LLVMPtr{T,A}) where {T,A}   = VEDeviceVector{T,A}((Int(len),), p)

## getters

Base.pointer(a::VEDeviceArray) = a.ptr
Base.pointer(a::VEDeviceArray, i::Integer) =
    pointer(a) + (i - 1) * Base.elsize(a)

Base.elsize(::Type{<:VEDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::VEDeviceArray) = g.shape
Base.length(g::VEDeviceArray) = prod(g.shape)

## conversions

Base.unsafe_convert(::Type{LLVMPtr{T,A}}, a::VEDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

## indexing intrinsics

# NOTE: these intrinsics are now implemented using plain and simple pointer operations;
#       when adding support for isbits union arrays we will need to implement that here.

# TODO: how are allocations aligned by the level zero API? keep track of this
#       because it enables optimizations like Load Store Vectorization
#       (cfr. shared memory and its wider-than-datatype alignment)

@inline function arrayref(A::VEDeviceArray{T}, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_load(pointer(A), index, Val(align))
end

@inline function arrayset(A::VEDeviceArray{T}, x::T, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_store!(pointer(A), x, index, Val(align))
    return A
end

@inline function const_arrayref(A::VEDeviceArray{T}, index::Int) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    unsafe_cached_load(pointer(A), index, Val(align))
end

## indexing

Base.@propagate_inbounds Base.getindex(A::VEDeviceArray{T}, i1::Int) where {T} =
    arrayref(A, i1)
Base.@propagate_inbounds Base.setindex!(A::VEDeviceArray{T}, x, i1::Int) where {T} =
    arrayset(A, convert(T,x)::T, i1)

Base.IndexStyle(::Type{<:VEDeviceArray}) = Base.IndexLinear()

## const indexing

"""
    Const(A::VEDeviceArray)

Mark a VEDeviceArray as constant/read-only. The invariant guaranteed is that you will not
modify an VEDeviceArray for the duration of the current kernel.

This API can only be used on devices with compute capability 3.5 or higher.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct Const{T,N} <: DenseArray{T,N}
    a::VEDeviceArray{T,N}
end
Base.Experimental.Const(A::VEDeviceArray) = Const(A)

Base.IndexStyle(::Type{<:Const}) = IndexLinear()
Base.size(C::Const) = size(C.a)
Base.axes(C::Const) = axes(C.a)
Base.@propagate_inbounds Base.getindex(A::Const, i1::Int) = const_arrayref(A.a, i1)

# deprecated
Base.@propagate_inbounds ldg(A::VEDeviceArray, i1::Integer) = const_arrayref(A, Int(i1))

## other

Base.show(io::IO, a::VEDeviceVector) =
    print(io, "$(length(a))-element device array at $(pointer(a))")
Base.show(io::IO, a::VEDeviceArray) =
    print(io, "$(join(a.shape, '×')) device array at $(pointer(a))")

Base.show(io::IO, mime::MIME"text/plain", a::VEDeviceArray) = show(io, a)

@inline function Base.iterate(A::VEDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end


text:
# Test for GaussLobatto()
n = 2
x,w = gausslobatto(n)
@test_approx_eq x[1] -1.
@test_approx_eq x[2] 1.
@test_approx_eq w[1] 1.
@test_approx_eq w[2] 1.

n = 3
x,w = gausslobatto(n)
@test_approx_eq x[1] -1.
@test abs(x[2])<1e-15
@test_approx_eq x[3] 1.
@test_approx_eq w[1] 1./3
@test_approx_eq w[2] 4./3
@test_approx_eq w[3] 1./3

tol = 1e-14
n = 42
x,w = gausslobatto(n)
@test ( length(x) == n) && ( length(w) == n )
@test abs(x[37] - 0.922259214258616) < tol
@test abs(w[37] - 0.029306411216166) < tol
@test ( x[1] == -1 && x[n] == 1 )

@test_approx_eq dot( w,(x.^2)) 2/3
@test_approx_eq dot( w,exp(x)) exp(1)-exp(-1)


text:
module TestPhylo
using Test

using Phylo

@testset "Deprecations" begin
    t = NamedTree()
    @test_deprecated name = addnode!(t)
    @test_deprecated addnode!(t, "Second")
    @test_deprecated addbranch!(t, name, "Second")
    @test_deprecated branch!(t, name)
    @test nleaves(t) == 2
    @test nroots(t) == 1
    @test nbranches(t) == 2
    @test_deprecated treeiter(t)
    @test_deprecated treenameiter(t)
end

end


text:
#= DSP is removed from test dependencies in this version.
import MechanicalSketch
import MechanicalSketch: empty_figure, PALETTE, O, HE, WI, EM, finish, ∙, Point
import MechanicalSketch: @import_expand, set_scale_sketch, settext, scale_pt_to_unit
import MechanicalSketch: noise, normalize_datarange, place_image, BinLegend, @layer
import MechanicalSketch: poly, dimension_aligned, sethue, arrow, circle, prettypoly, Greys_9
using  MechanicalSketch: MechanicalUnits, upreferred, ComplexQuantity, Quantity
import MechanicalUnits:  dimension, NoDims
import DSP.Periodograms: spectrogram, Spectrogram
import DSP.Util:         nextfastfft
import DSP.Windows:      tukey
import Statistics:       mean
#let
if !@isdefined m²
    @import_expand(~m, s)
end
empty_figure(filename = joinpath(@__DIR__, "test_29.png"));

curpoint = O + (-WI / 2 + EM, -HE / 2 + 3EM)
velocities = range(0.05m/s, stop= 0.5m/s, length = 3)
settext("
The flow field visualization should have a linear relationship between velocity and luminosity.
We need to control the noise image wavelength spectrum to get that.
", curpoint)

curpoint += (0, 4EM)
settext("
For reference, let's first illustrate two equal amplitude wavelengths, λ, as greyscale:
    cos(2π x / 0.5m) + 0.5cos(2π x / 0.05m)
", curpoint, markup = true)

curpoint += (0, 3EM)
n_pixels = round(Int, 0.9WI)
physwidth = 2.0m
set_scale_sketch(physwidth, n_pixels)
length_one_pixel = physwidth / n_pixels
twowave(x) = cos(2π∙x / 0.5m) + 0.5cos(2π∙x / 0.05m)
no = [twowave(x) for x in (1:n_pixels)*length_one_pixel]
nno = normalize_datarange(no)
greylegend = BinLegend(;maxlegend = 1.0, noofbins = 128, colorscheme = reverse(Greys_9))
@layer begin
    place_image(curpoint, greylegend.(nno); centered = false)
    place_image(curpoint + (0,1), greylegend.(nno); centered = false)
    place_image(curpoint + (0,2), greylegend.(nno); centered = false)
    sethue(PALETTE[3])
    dimension_aligned(curpoint, curpoint + (physwidth, 0.0m), offset = -2.5EM)
    dimension_aligned(curpoint, curpoint + (0.5m, 0.0m), offset = -EM, toextension = (0EM, 7EM))
end

curpoint += (0, 3EM)
settext("
The same greyscale values shown as a curve:", curpoint)
curpoint += (0, EM)
points = map(1:n_pixels, nno) do xul, yul
    (xul, -yul * EM)
end
poly(curpoint .+ points, :stroke)

curpoint += (0, 2EM)
settext("
Amplitude-wavelength spectrum from sampling the same values:", curpoint, markup=true)

curpoint += (0, 3EM)
include("test_functions_29.jl")
spectrum = sampledspectrum_29(no .- 0.5, length_one_pixel, maximumwavelength = 0.5m)
draw_spectrum_29(curpoint, spectrum)

curpoint += (0, 2EM)

(λ_1, a_1), (λ_2, a_2) = two_largest_maxima_amplitude_wavelength(spectrum)
settext("
    Local maxima are (λ, amplitude) = [$((λ_1, a_1)), $((λ_2, a_2))]",
    curpoint, markup=true)

curpoint += (0, 2EM)
settext("
We can sample from a longer example and get a finer spectrum resolution:", curpoint)

curpoint += (0, 3EM)
longnoise = [twowave(x) for x in (1:100n_pixels)*length_one_pixel]
spectrum = sampledspectrum_29(longnoise, length_one_pixel, maximumwavelength = 0.5m)
draw_spectrum_29(curpoint, spectrum)

curpoint += (0, 2EM)
(λ_1, a_1), (λ_2, a_2) = two_largest_maxima_amplitude_wavelength(spectrum)
settext("
    Local maxima are (λ, amplitude) = [$((λ_1, a_1)), $((λ_2, a_2))]",
    curpoint, markup=true)

finish()
set_scale_sketch()
#end
=#


text:
using Atum
using Atum.Advection

using Test
using Printf
using StaticArrays: SVector
using LinearAlgebra: norm

if !@isdefined integration_testing
  const integration_testing = parse(
    Bool,
    lowercase(get(ENV, "ATUM_INTEGRATION_TESTING", "false")),
  )
end

function wave(law, x⃗, t)
  ρ = 2 + sin(π * (sum(x⃗) - sum(constants(law).u⃗) * t))
  SVector(ρ)
end

function run(A, FT, N, K; volume_form=WeakForm())
  Nq = N + 1

  law = AdvectionLaw{FT, 1}()

  cell = LobattoCell{FT, A}(Nq)
  v1d = range(FT(-1), stop=FT(1), length=K+1)
  grid = brickgrid(cell, (v1d,); periodic=(true,))

  dg = DGSEM(; law, grid, volume_form,
               surface_numericalflux = RusanovFlux())

  cfl = FT(1 // 2)
  dt = cfl * step(v1d) / N / norm(constants(law).u⃗)
  timeend = FT(0.7)

  q = wave.(Ref(law), points(grid), FT(0))

  @info @sprintf """Starting
  N           = %d
  K           = %d
  volume_form = %s
  norm(q)     = %.16e
  """ N K volume_form weightednorm(dg, q)

  odesolver = LSRK54(dg, q, dt)
  solve!(q, timeend, odesolver)

  qexact = wave.(Ref(law), points(grid), timeend)
  errf = weightednorm(dg, q .- qexact)

  @info @sprintf """Finished
  norm(q)      = %.16e
  norm(q - qe) = %.16e
  """ weightednorm(dg, q) errf
  errf
end

let
  A = Array
  FT = Float64
  N = 4

  expected_error = Dict()

  #form, lev
  expected_error[WeakForm(), 1] = 1.5238393162339734e-04
  expected_error[WeakForm(), 2] = 4.8377348010293410e-06
  expected_error[WeakForm(), 3] = 1.4063471307765498e-07
  expected_error[WeakForm(), 4] = 4.5781066120767573e-09

  expected_error[FluxDifferencingForm(CentralFlux()), 1] = 1.5238393162391597e-04
  expected_error[FluxDifferencingForm(CentralFlux()), 2] = 4.8377348011318472e-06
  expected_error[FluxDifferencingForm(CentralFlux()), 3] = 1.4063471309426943e-07
  expected_error[FluxDifferencingForm(CentralFlux()), 4] = 4.5781064622952347e-09

  nlevels = integration_testing ? 4 : 1

  @testset for volume_form in (WeakForm(), FluxDifferencingForm(CentralFlux()))
    errors = zeros(FT, nlevels)
    for l in 1:nlevels
      K = 5 * 2 ^ (l - 1)
      errors[l] = run(A, FT, N, K; volume_form)
      @test errors[l] ≈ expected_error[volume_form, l]
    end
    if nlevels > 1
      rates = log2.(errors[1:(nlevels-1)] ./ errors[2:nlevels])
      @info "Convergence rates\n" *
        join(["rate for levels $l → $(l + 1) = $(rates[l])" for l in 1:(nlevels - 1)], "\n")
      @test rates[end] ≈ N + 1 atol=0.1
    end
  end
end


text:
@testset "weighted_mean" begin
    for i in 1:10
        len = rand(1:80)
        x = rand(1:100.0,len)
        y = rand(1:100.0,len)
        @test length(NighttimeLights.weighted_mean(x, y)) == 1
    end
end


text:
@with_kw struct Analysis{T, N}
    data::SelectedData{T, N}
    compute_error = nothing
    x = nothing
    y = nothing
    z = nothing
    xfunc = nothing
    yfunc = nothing
    zfunc = nothing
    axis_type = nothing
    smoother = nothing
    plot = nothing
    plot_kwargs = []
end

function Analysis(a::Analysis; kwargs...)
    d = Dict(kwargs)
    Analysis((get(d, f, getfield(a, f)) for f in fieldnames(a))...)
end

struct StatPlotsRecipe; end
struct GroupedError; end

function analysistype(a)
    a.compute_error !== nothing && return GroupedError
    (a.y in colnames(a.data.table) || a.y === nothing) ? StatPlotsRecipe : GroupedError
end

process(a::Analysis) = process(analysistype(a), a)


text:
# using Random
### This doesn't error!
# @testset "does the expand crash?" begin

#     rstructs = readmaster("data/master.json")
#     rids = readids("data/taxon_reactions/1234567892.json")
#     sids = readkeyedids("data/seeds/seeds.json")["Enceladus_20-SAFR-032"]
#     tids = readkeyedids("data/seeds/seeds.json")["targets_Freilich09"]

#     S = System(rstructs,
#                rids,
#                sids,
#                tids)

#     expand(S,"results/simple_test.json")

# end

### This doesn't error!
using SparseArrays

@testset "minimal seed sets" begin

    rstructs = readmaster("data/master.json")
    rids = readids(["R09735","R08827"])
    sids = readids(["C19848","C06232","C18237","C00020","C00001","C00355","C08538","C08543","C00001"])
    tids = readids(["C06232"])
    # sids = readkeyedids("data/seeds/seeds.json")["Enceladus_20-SAFR-032"]
    # tids = readkeyedids("data/seeds/seeds.json")["targets_Freilich09"]

    S = System(rstructs=rstructs,
               rids=rids,
               sids=sids,
               tids=tids)

    ## Below tests just make sure that nothing errors
    @test isa(find_minimal_seed_set(S,write_path="results/simple_minimalseed_test.json"), Any)
    @test isa(find_minimal_seed_set(S), Any)

    allowed_forward = sparsevec([true,false])
    allowed_backward = sparsevec([false,true])

    @test isa(find_minimal_seed_set(S, allowed_forward = allowed_forward, allowed_backward = allowed_backward), Any)

    # myrng = MersenneTwister(1234)
    # [rand(myrng,[0,1]) for i in 1:length(rids)]

end

@testset "test BioXP 'private' functions" begin

    rstructs = readmaster("data/master-subset-3rxns.json")
    rids = readids(["R09735","R08827"])

    @testset "identify_biosystem_compounds" begin

        expected_cids = ["C19848","C06232","C18237","C00020","C00001","C00355","C08538","C08543","C00001"]

        @test Set(expected_cids) == Set(BioXP.identify_biosystem_compounds(rstructs,rids))
        @test 8 == length(BioXP.identify_biosystem_compounds(rstructs,rids))
    end

    @testset "matrixify_compounds" begin

        cids = BioXP.identify_biosystem_compounds(rstructs,rids)
        expected_cid_order = ["C19848", "C06232", "C18237", "C00020", "C00001", "C00355", "C08538", "C08543"]
        # leftids = ["C19848","C06232"] ["C00355","C08538"]
        # rightids = ["C18237","C00020","C00001"] ["C08543","C00001"]
        expected_R = sparse([1 0;1 0;0 0;0 0;0 0;0 1;0 1;0 0])
        expected_P = sparse([0 0;0 0;1 0;1 0;1 1;0 0;0 0;0 1])

        R,P = BioXP.matrixify_compounds(rstructs,cids,rids)
        @test expected_R == R
        @test expected_P == P
        @test isa(R,SparseMatrixCSC)
        @test isa(P,SparseMatrixCSC)

    end

    @testset "matrixify_seeds" begin

        sids = ["C19848","C06232","C00001"]
        cids = BioXP.identify_biosystem_compounds(rstructs,rids)
        expected_cid_order = ["C19848", "C06232", "C18237", "C00020", "C00001", "C00355", "C08538", "C08543"]

        expected_x = sparse([1,1,0,0,1,0,0,0])
        x = BioXP.matrixify_seeds(sids, cids)
        @test expected_x == x
        @test isa(x, SparseVector)

    end

    @testset "matrixify_targets" begin

        tids = ["C18237","C00020","C00001"]
        cids = BioXP.identify_biosystem_compounds(rstructs,rids)
        expected_cid_order = ["C19848", "C06232", "C18237", "C00020", "C00001", "C00355", "C08538", "C08543"]

        expected_t = [0,0,1,1,1,0,0,0]
        t = BioXP.matrixify_targets(cids,tids)
        @test expected_t == t
        @test isa(t, SparseVector)

    end

    ### test BioXP simple expansion"

    @testset "expandmatrices" begin

        sids = ["C19848","C06232","C00001"]
        cids = BioXP.identify_biosystem_compounds(rstructs,rids)
        R,P = BioXP.matrixify_compounds(rstructs,cids,rids)
        # leftids = ["C19848","C06232"] ["C00355","C08538"]
        # rightids = ["C18237","C00020","C00001"] ["C08543","C00001"]
        # expected_cid_order = ["C19848", "C06232", "C18237", "C00020", "C00001", "C00355", "C08538", "C08543"]
        x = BioXP.matrixify_seeds(sids, cids)

        ## 1st X column should be same as seeds
        ## 2nd X column should also include new compounds
        ##     which can be made from first column
        expected_X = [[1,1,0,0,1,0,0,0], # seeds
                      [1,1,1,1,1,0,0,0], # seeds+new compounds
                      [1,1,1,1,1,0,0,0]] # can't make anymore new compounds

        expected_Y = [[1,0],
                      [1,0]]

        X,Y = BioXP.expandmatrices(R,P,x)
        @test expected_X == X
        @test expected_Y == Y
        @test isa(X[1],SparseVector{Int64,Int64})
        @test isa(Y[1],SparseVector{Int64,Int64})

    end

    @testset "expand" begin

    end

end

## Next steps:
## - check old results in `test/data/submission_results` against results from running new code
##    - the formatting should be identical. but perhaps there's no guarantee that the ordering of the compounds will be the same.
##      So I will have to figure out some way to test.
@testset "test expanding original master prokaryotes" begin
    newformat_dir = "../test/data/submission_results/archaea_jgi_newformat/"
    archaea_simpleresults_P_dir = "../test/data/submission_results/archaea_simpleresults_P/"
    seeds_and_targets = readkeyedids("../test/data/seeds/seeds.json")

    sids = seeds_and_targets["Enceladus_20-SAFR-032_P"]
    tids = seeds_and_targets["targets_Freilich09"]
    rstructs = readmaster("../test/data/master_from_redges-og-submission.json")

    for fname in readdir(newformat_dir)
        rids = readids(newformat_dir*fname)
        x, t, cids, X, Y = expand(rstructs,rids,sids,tids=tids) ## new results

        e = JSON.parsefile(archaea_simpleresults_P_dir*fname) ## expected results

        @test length(x) == length(e["x"])
        @test length(t) == length(e["t"])
        @test length(cids) == length(e["compounds"])
        @test length(X) == length(e["X"])
        @test length(Y) == length(e["Y"])

        # Compare calculated against expected results
        for gen in 1:length(e["X"])

            expected_cpds = Set(zip(e["compounds"],e["X"][gen]))
            calcualted_cpds = Set(zip(cids,X[gen]))

            @test expected_cpds == calcualted_cpds

            if gen != length(e["X"])
                expected_rxns = Set(zip(e["reactions"],e["Y"][gen]))
                calculated_rxns = Set(zip(rids,Y[gen]))

                @test expected_rxns == calculated_rxns
            end
        end
    end
end

## This test takes about 60 seconds
# @testset "test expanding original master KEGG w/P" begin
#     expected_results_fname = "../test/data/submission_results/kegg_edge_json_P/reaction_edges_P.json"
#     seeds_and_targets = readkeyedids("../test/data/seeds/seeds.json")
#     sids = seeds_and_targets["Enceladus_20-SAFR-032_P"]
#     tids = seeds_and_targets["targets_Freilich09"]
#     rstructs = readmaster("../test/data/master_from_redges-og-submission.json")

#     rids = readids("../test/data/taxon_reactions/kegg.json")
#     x, t, cids, X, Y = expand(rstructs,rids,sids,tids) ## new results

#     e = JSON.parsefile(expected_results_fname) ## expected results

#     @test length(x) == length(e["x"])
#     @test length(t) == length(e["t"])
#     @test length(cids) == length(e["compounds"])
#     @test length(X) == length(e["X"])
#     @test length(Y) == length(e["Y"])

#     # Compare calculated against expected results
#     for gen in 1:length(e["X"])

#         expected_cpds = Set(zip(e["compounds"],e["X"][gen]))
#         calcualted_cpds = Set(zip(cids,X[gen]))

#         @test expected_cpds == calcualted_cpds

#         if gen != length(e["X"])
#             expected_rxns = Set(zip(e["reactions"],e["Y"][gen]))
#             calculated_rxns = Set(zip(rids,Y[gen]))

#             @test expected_rxns == calculated_rxns
#         end
#     end
# end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
using JuMP, Gurobi
include("function.jl")

SolverBrendon.SolveMIP(m)


text:
export SVHN2

"""
The Street View House Numbers (SVHN) Dataset

- Authors: Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, Andrew Y. Ng
- Website: http://ufldl.stanford.edu/housenumbers

SVHN was obtained from house numbers in Google Street View
images. As such they are quite diverse in terms of orientation
and image background. Similar to MNIST, SVHN has 10 classes (the
digits 0-9), but unlike MNIST there is more data and the images
are a little bigger (32x32 instead of 28x28) with an additional
RGB color channel. The dataset is split up into three subsets:
73257 digits for training, 26032 digits for testing, and 531131
additional to use as extra training data.

## Interface

- [`SVHN2.traintensor`](@ref), [`SVHN2.trainlabels`](@ref), [`SVHN2.traindata`](@ref)
- [`SVHN2.testtensor`](@ref), [`SVHN2.testlabels`](@ref), [`SVHN2.testdata`](@ref)
- [`SVHN2.extratensor`](@ref), [`SVHN2.extralabels`](@ref), [`SVHN2.extradata`](@ref)

## Utilities

- [`SVHN2.download`](@ref)
- [`SVHN2.classnames`](@ref)
- [`SVHN2.convert2image`](@ref)
"""
module SVHN2
    using DataDeps
    using ColorTypes
    using FixedPointNumbers
    using MAT: matopen, matread
    using ..MLDatasets: bytes_to_type, datafile, download_dep, download_docstring,
                        _colorview, _channelview

    export

        traintensor,
        testtensor,
        extratensor,

        trainlabels,
        testlabels,
        extralabels,

        traindata,
        testdata,
        extradata,

        convert2image,

        download

    @deprecate convert2features reshape

    const DEPNAME = "SVHN2"
    const TRAINDATA = "train_32x32.mat"
    const TESTDATA  = "test_32x32.mat"
    const EXTRADATA = "extra_32x32.mat"
    const CLASSES = [1, 2, 3, 4, 5, 6, 7, 8, 9, 0]

    """
        download([dir]; [i_accept_the_terms_of_use])

    Trigger the (interactive) download of the full dataset into
    "`dir`". If no `dir` is provided the dataset will be
    downloaded into "~/.julia/datadeps/$DEPNAME".

    This function will display an interactive dialog unless
    either the keyword parameter `i_accept_the_terms_of_use` or
    the environment variable `DATADEPS_ALWAY_ACCEPT` is set to
    `true`. Note that using the data responsibly and respecting
    copyright/terms-of-use remains your responsibility.
    """
    download(args...; kw...) = download_dep(DEPNAME, args...; kw...)

    include("interface.jl")
    include("utils.jl")

    function __init__()
        register(DataDep(
            DEPNAME,
            """
            Dataset: The Street View House Numbers (SVHN) Dataset
            Authors: Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, Andrew Y. Ng
            Website: http://ufldl.stanford.edu/housenumbers
            Format: Cropped Digits (Format 2 on the website)
            Note: for non-commercial use only

            [Netzer et al., 2011]
                Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu, Andrew Y. Ng
                "Reading Digits in Natural Images with Unsupervised Feature Learning"
                NIPS Workshop on Deep Learning and Unsupervised Feature Learning 2011

            The dataset is split up into three subsets: 73257
            digits for training, 26032 digits for testing, and
            531131 additional to use as extra training data.

            The files are available for download at the official
            website linked above. Note that using the data
            responsibly and respecting copyright remains your
            responsibility. For example the website mentions that
            the data is for non-commercial use only. Please read
            the website to make sure you want to download the
            dataset.
            """,
            "http://ufldl.stanford.edu/housenumbers/" .* [TRAINDATA, TESTDATA, EXTRADATA],
            "2fa3b0b79baf39de36ed7579e6947760e6241f4c52b6b406cabc44d654c13a50"
        ))
    end
end


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

mutable struct Encryption <: SwaggerModel
    services::Any # spec type: Union{ Nothing, EncryptionServices } # spec name: services
    keySource::Any # spec type: Union{ Nothing, String } # spec name: keySource
    requireInfrastructureEncryption::Any # spec type: Union{ Nothing, Bool } # spec name: requireInfrastructureEncryption
    keyvaultproperties::Any # spec type: Union{ Nothing, KeyVaultProperties } # spec name: keyvaultproperties

    function Encryption(;services=nothing, keySource="Microsoft.Storage", requireInfrastructureEncryption=nothing, keyvaultproperties=nothing)
        o = new()
        validate_property(Encryption, Symbol("services"), services)
        setfield!(o, Symbol("services"), services)
        validate_property(Encryption, Symbol("keySource"), keySource)
        setfield!(o, Symbol("keySource"), keySource)
        validate_property(Encryption, Symbol("requireInfrastructureEncryption"), requireInfrastructureEncryption)
        setfield!(o, Symbol("requireInfrastructureEncryption"), requireInfrastructureEncryption)
        validate_property(Encryption, Symbol("keyvaultproperties"), keyvaultproperties)
        setfield!(o, Symbol("keyvaultproperties"), keyvaultproperties)
        o
    end
end # type Encryption

const _property_map_Encryption = Dict{Symbol,Symbol}(Symbol("services")=>Symbol("services"), Symbol("keySource")=>Symbol("keySource"), Symbol("requireInfrastructureEncryption")=>Symbol("requireInfrastructureEncryption"), Symbol("keyvaultproperties")=>Symbol("keyvaultproperties"))
const _property_types_Encryption = Dict{Symbol,String}(Symbol("services")=>"EncryptionServices", Symbol("keySource")=>"String", Symbol("requireInfrastructureEncryption")=>"Bool", Symbol("keyvaultproperties")=>"KeyVaultProperties")
Base.propertynames(::Type{ Encryption }) = collect(keys(_property_map_Encryption))
Swagger.property_type(::Type{ Encryption }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Encryption[name]))}
Swagger.field_name(::Type{ Encryption }, property_name::Symbol) =  _property_map_Encryption[property_name]

const _allowed_Encryption_keySource = ["Microsoft.Storage", "Microsoft.Keyvault"]

function check_required(o::Encryption)
    true
end

function validate_property(::Type{ Encryption }, name::Symbol, val)
    if name === Symbol("keySource")
        Swagger.validate_param(name, "Encryption", :enum, val, _allowed_Encryption_keySource)
    end
end


text:
using SkewDist, Distributions, Base.Test

println("Testing operators...")

srand(2)

N = 1000000

d = 3
ξ = rand(d)
α = rand(d)
Ω = rand_pd(d)/(2*d)
ν = 4.0
A = rand(d-1,d)
x = rand(d)
dist = MvSkewTDist(ξ, Ω, α, 4.0)

println("\tMatrix transformation...")
dist2 = A * dist
y = A * rand(dist, N)
@test_approx_eq_eps mean(y, 2) mean(dist2) 1e-1
@test_approx_eq_eps var(y,2) var(dist2) 1e-1
@test_approx_eq_eps Base.covm(y, mean(y, 2), 2) cov(dist2) 1e-1

println("\tLinear combination...")
dist3 = x * dist
z = x'rand(dist, N)
@test_approx_eq_eps mean(z) mean(dist3) 1e-2
@test_approx_eq_eps var(z) var(dist3) 1e-2

println("\tLinear and matrix transformations consistent...")
dist4 = x' * dist
@test_approx_eq dist4.ξ[1] dist3.ξ
@test_approx_eq sqrt(dist4.Ω.mat[1,1]) dist3.ω
@test_approx_eq dist4.α[1] dist3.α


text:
export Interval

using ..Points

struct Interval <: AbstractShape{1,2}
    start::Float64
    stop::Float64
    origin::Point{1}
    a::Float64
end
(i::Interval)(p::Point) = i.start ≤ p.x ≤ i.stop
(i::Interval)(x::Real,y...) = i(Point(x,y...))

Interval(start,stop;ref::Symbol=:center) = Interval(start,stop,generate_origin(start,stop,ref),abs(stop-start))

function generate_origin(start::Real,stop::Real,ref::Symbol)
    left,right = minmax(start,stop)
    if Base.sym_in(ref,(:left,:l,:Left,:L))
        origin = Point(left)
    elseif Base.sym_in(ref,(:right,:r,:Right,:R))
        origin = Point(right)
    else
        origin = Point((left+right)/2)
    end
    return origin
end

# struct BraggStack{N,T}
#     origin::Float64
#     depths::NTuple{N,Float64}
#     indices::NTuple{N,T}
#     periods::Int
#     integrated_depths::NTuple{N,Float64}
#
#     function BraggStack(origin,depths,indices,periods)
#         new{length(depths),eltype(indices)}(origin,depths,indices,periods,tuple(cumsum(depths...)))
#     end
# end
# function BraggStack(p::Point)
#     0 ≤ (p.x-origin)
# end


text:
"""
    solve_analytical(setup, experiment[, volumes])

Compute the signal in a multilayered cylinder or sphere using an analytical
matrix formalism solution.

This function is based on the following articles and corresponding code:
    [1] D. S. Grebenkov, NMR Survey of Reflected Brownian Motion,
        Rev. Mod.Phys. 79, 1077 (2007)
    [2] D. S. Grebenkov, Pulsed-gradient spin-echo monitoring of restricted diffusion in
        multilayered structures, J. Magn. Reson. 205, 181-195 (2010).
"""
function solve_analytical(setup, experiment, volumes = nothing)

    # Sizes
    namplitude = length(experiment.gradient.values)
    nsequence = length(experiment.gradient.sequences)

    # Extract parameters
    sequences = experiment.gradient.sequences
    rmean = (setup.rmin + setup.rmax) / 2

    # Q-values and b-values
    if experiment.gradient.values_type == "q"
        qvalues = repeat(experiment.gradient.values, 1, nsequence)
        bvalues =
            experiment.gradient.values .^ 2 .* bvalue_no_q.(experiment.gradient.sequences)'
    else
        bvalues = repeat(experiment.gradient.values, 1, nsequence)
        qvalues =
            .√(experiment.gradient.values ./ bvalue_no_q.(experiment.gradient.sequences)')
    end

    if isa(setup, CylinderSetup)
        dim = 2
        get_vol = r -> π * r^2 * setup.height
    else
        dim = 3
        get_vol = r -> 4 * π / 3 * r^3
    end

    # Get OUT parameters
    ρ_out = setup.ρ_out
    r_out = (setup.rmin + setup.rmax) / 2
    D_out = tr(setup.D_out) / 3
    T_out = setup.T₂_out
    W_out = 0

    # Get IN parameters
    if setup.include_in
        ρ_in = setup.ρ_in
        r_in = setup.in_ratio * r_out
        D_in = tr(setup.D_in) / 3
        W_in = setup.κ_in_out
        T_in = setup.T₂_in
    else
        # Empty IN
        ρ_in = []
        r_in = []
        D_in = []
        W_in = []
        T_in = []
    end

    # Get ECS parameters
    if setup.ecs_shape != "no_ecs"
        # Include ECS
        ρ_ecs = setup.ρ_ecs
        r_ecs = r_out + setup.ecs_ratio * rmean
        D_ecs = tr(setup.D_ecs) / 3
        W_ecs = setup.κ_ecs
        T_ecs = setup.T₂_ecs

        # Add interface permeability between out and ecs
        W_out = setup.κ_out_ecs
    else
        # Empty ECS
        ρ_ecs = []
        r_ecs = []
        D_ecs = []
        W_ecs = []
        T_ecs = []

        if dim == 3
            # Add surface relaxivity for the outermost sphere ("out")
            W_out = setup.κ_out
        end
    end

    # Create parameters structure
    ρ = vcat(ρ_in, ρ_out, ρ_ecs)
    r = vcat(r_in, r_out, r_ecs) * 1e-6
    D = vcat(D_in, D_out, D_ecs) * 1e-6
    W = vcat(W_in, W_out, W_ecs) * 1
    T = vcat(T_in, T_out, T_ecs) * 1e-6
    params = (; ρ, r, D, W, T, dim)

    L = length(D)

    # Get spherical volumes
    volumes_exact = get_vol.(r * 1e6)

    # Subtract inner volumes (each volume is contained within the next)
    volumes_exact[2:end] = volumes_exact[2:end] - volumes_exact[1:end-1]

    # Volumes
    if isnothing(volumes)
        display("Using exact volumes")
        volumes = volumes_exact
    else
        display("Exact volumes:")
        display(volumes_exact)
        display("Finite element volumes:")
        display(volumes)
        display("Relative error:")
        display(abs.(volumes - volumes_exact) ./ volumes_exact)
        display("Using finite element volumes")
    end

    # Volume weighted quantities
    D_mean = D' * volumes / sum(volumes)
    S0 = volumes' * ρ

    # Upper Laplace eigenvalue limit
    eiglim = π^2 * D_mean / experiment.analytical.length_scale^2

    # Compute radial and angular eigenvalues (λ and ν), with
    # α^2 = λ and n^2 = ν for cylinders and n(n+1) = ν for spheres
    dα = √experiment.analytical.eigstep * 1e3
    α_min = dα / 100
    α_max = √eiglim * 1e6
    α, n = find_α(params, α_min, α_max, dα)

    N = length(α)
    λ = α .^ 2
    Λ = diagm(λ)

    β = zeros(N)
    J = zeros(N)
    for m = 1:N
        β[m] = compute_β(params, α[m], n[m])
        if n[m] == 0
            J[m] = compute_int_J(params, α[m])[1]
        end
    end

    if dim == 2
        coeff = 2
    else
        coeff = √6
    end
    U = coeff * J .* β

    B = zeros(N, N)
    for m1 = 1:N
        for m2 = m1:N
            if abs(n[m1] - n[m2]) == 1
                K = compute_int_K(params, α[m1], n[m1], α[m2], n[m2])
                if dim == 2
                    ε = √(1 + (n[m1] == 0) + (n[m2] == 0))
                else
                    ε = (n[m1] + n[m2] + 1) / √((2 * n[m1] + 1) * (2 * n[m2] + 1))
                end
                B[m1, m2] = ε * β[m1] * β[m2] * K
                B[m2, m1] = B[m1, m2]
            end
        end
    end

    Bri = zeros(N, N, L)
    for m1 = 1:N
        for m2 = m1:N
            if n[m1] == n[m2]
                int_I = compute_int_I(params, α[m1], α[m2], n[m1])
                Bri[m1, m2, :] = 2 * β[m1] * β[m2] * int_I
                Bri[m2, m1, :] = Bri[m1, m2, :]
            end
        end
    end

    Br = zeros(N, N)
    for i = 1:L
        Br .+= Bri[:, :, i] / T[i]
    end

    # Iterate over experiments
    signal = zeros(namplitude, nsequence)
    for iseq = 1:nsequence, iamp = 1:namplitude
        seq = sequences[iseq]
        w = qvalues[iamp, iseq] * 1e12 * r[end]

        δ = seq.δ * 1e-6
        Δ = seq.Δ * 1e-6

        # Compute signal attenuation
        E = (
            U' *
            exp(-δ * (Λ + Br + im * w * B)) *
            exp(-(Δ - δ) * (Λ + Br)) *
            exp(-δ * (Λ + Br - im * w * B)) *
            U
        )
        signal[iamp, iseq] = E * S0
    end

    signal
end


text:
# ---
# title: 865. Smallest Subtree with all the Deepest Nodes
# id: problem865
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/smallest-subtree-with-all-the-deepest-nodes/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, the depth of each node is **the shortest
# distance to the root**.
#
# Return _the smallest subtree_ such that it contains **all the deepest nodes**
# in the original tree.
#
# A node is called **the  deepest** if it has the largest depth possible among
# any node in the entire tree.
#
# The **subtree** of a node is tree consisting of that node, plus the set of all
# descendants of that node.
#
# **Note:** This question is the same as 1123:
# <https://leetcode.comhttps://leetcode.com/problems/lowest-common-ancestor-of-deepest-leaves/>
#
#
#
# **Example 1:**
#
# ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/07/01/sketch1.png)
#
#
#
#     Input: root = [3,5,1,6,2,0,8,null,null,7,4]
#     Output: [2,7,4]
#     Explanation: We return the node with value 2, colored in yellow in the diagram.
#     The nodes coloured in blue are the deepest nodes of the tree.
#     Notice that nodes 5, 3 and 2 contain the deepest nodes in the tree but node 2 is the smallest subtree among them, so we return it.
#
#
# **Example 2:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#     Explanation: The root is the deepest node in the tree.
#
#
# **Example 3:**
#
#
#
#     Input: root = [0,1,3,null,2]
#     Output: [2]
#     Explanation: The deepest node in the tree is 2, the valid subtrees are the subtrees of nodes 2, 1 and 0 but the subtree of node 2 is the smallest.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree will be in the range `[1, 500]`.
#   * `0 <= Node.val <= 500`
#   * The values of the nodes in the tree are **unique**.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
import Base.@deprecate

import Base:
    dirname,
    ispath,
    realpath,
    normpath,
    abspath,
    relpath,
    filemode,
    isabspath,
    mkpath,
    mv,
    rm

@deprecate dirname(fp::AbstractPath) parent(fp)
@deprecate ispath(fp::AbstractPath) exists(fp)
@deprecate realpath(fp::AbstractPath) real(fp)
@deprecate normpath(fp::AbstractPath) norm(fp)
@deprecate abspath(fp::AbstractPath) abs(fp)
@deprecate relpath(fp::AbstractPath) relative(fp)
@deprecate filemode(fp::AbstractPath) mode(fp)
@deprecate isabspath(fp::AbstractPath) isabs(fp)
@deprecate mkpath(fp::AbstractPath) mkdir(fp; recursive=true, exist_ok=true)
@deprecate parts(fp::AbstractPath) fp.segments
@deprecate drive(fp::Abstractpath) fp.drive
@deprecate root(fp::AbstractPath) fp.root
@deprecate anchor(fp::AbstractPath) fp.anchor


text:
using Douglass

tests = ["Douglass.jl"
		 ]

println("Running tests:")

for test in tests
	try
		include(test)
		println("\t\033[1m\033[32mPASSED\033[0m: $(test)")
	 catch e
	 	println("\t\033[1m\033[31mFAILED\033[0m: $(test)")
	 	showerror(stdout, e, backtrace())
	 	rethrow(e)
	 end
end

# don't do repl test on unix-- TODO improve this
if !Sys.isunix()
	try
		include("repl.jl")
		println("\t\033[1m\033[32mPASSED\033[0m: REPL.jl")
	 catch e
	 	println("\t\033[1m\033[31mFAILED\033[0m: REPL.jl")
	 	showerror(stdout, e, backtrace())
	 	rethrow(e)
	 end
end


text:
module Utils

# Units
const second = float(1)    # second
const minute = 60second    # minute
const ms = 1e-3second      # millisecond
const Hz = inv(second)     # Herz
const kHz = 1e3Hz          # kilohertz
const metre = float(1)     # meter
const cm = 0.01metre       # centimeter
const cm² = cm^2           # square centimeter
const mL = cm^3            # milliliter = cubic centimeter
const Liter = 1e3mL        # liter
const μL = 1E-6Liter
const pL = 1E-12Liter
const mM = float(1)
const Molar = 1000mM       # molar (1000 since the SI units is mM)
const μM = 1E-3mM          # micromolar
const nM = 1E-6mM          # nanomolar
const Amp = float(1)       # ampere
const mA = 1E-3Amp         # milliampere
const μA = 1E-6Amp         # micrpampere
const Volt = float(1)      # volt
const mV = 1E-3Volt        # millivolt
const T0 = 310.0            # Default temp (37C)
const F = 96485.0           # Faraday constant (c / mol)
const R = 8.314            # Ideal gas constant (K/mol)
const VT = R * T0 / F      # Thermal voltage (@37C) (Volts)
const iVT = inv(VT)        # Reciprocal of thermal voltage (@37C)

# Model Constants
const C_MIT = 1.812μM / mV        # Mitochondrial membrane capacitance
const F_M = 3E-4                # Frasction of free Ca in mitochondria
const F_I = 0.01                # Fraction of free Ca in cytosol
const V_I = 0.53                # Relative cytoplasmic volume
const V_MT = 0.06               # Relative mitochondrial volume
const V_MTX = 0.0144            # Relative mitochondrial matrix volume
const iVi = inv(V_I)
const iVmtx = inv(V_MTX)
const iVimtx = inv(V_MTX + V_I)
const iCmt = inv(C_MIT)

##################################
### Commonly-used functions
##################################

"Return `1+x` with the same type as x"
p_one(x) = one(x) + x

"Return `1-x` with the same type as x"
one_m(x) = one(x) - x

"""
Regular Hill function
"""
hill(x, k=one(x)) = x / (x + k)
hill(x, k, n) = hill(x^n, k^n)

"""
Repressive Hill function
"""
hillr(x, k=one(x)) = hill(k, x)
hillr(x, k, n) = hill(k, x, n)

"""
Logistic sigmoid function.
See scipy example https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.expit.html
"""
expit(x) = hillr(exp(-x))

"""
    exprel(x, em1 = expm1(x))

Returns `x / (exp(x)-1)` accurately when x is near zero.
See scipy example https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.exprel.html
Note the fraction is the inverse of `scipy.exprel()`
"""
exprel(x, em1=expm1(x)) = x / em1

"""
Signed sqare root
"""
sqrt_s(x) = flipsign(sqrt(abs(x)), x)

"""
Signed power
"""
pow_s(x, n) = flipsign(abs(x)^n, x)

"""
Signed Hill function
"""
hill_s(x, k, n) = hill(pow_s(x, n), pow_s(k, n))

end # module

max_stars_repo_name: ntumitolab/mitodyn-ode

max_issues_repo_name: ntumitolab/mitodyn-ode

max_forks_repo_name: NTUMitoLab/mitodyn-ode


text:
type EulerDiscretization <: AbstractDiscretization end

drift(euler::EulerDiscretization, process::StochasticProcess1D, t0::Float64, x0::Float64, dt::Float64) = drift(process, t0, x0) * dt

diffusion(euler::EulerDiscretization, process::StochasticProcess1D, t0::Float64, x0::Float64, dt::Float64) = diffusion(process, t0, x0) * sqrt(dt)

function variance(euler::EulerDiscretization, process::StochasticProcess1D, t0::Float64, x0::Float64, dt::Float64)
  sigma = diffusion(process, t0, x0)

  return sigma * sigma * dt
end


text:
module ConvInsight

using HDF5

export

    Logfile,
    log_metadata!,
    log_timeseries!

const ATTR_NETWORKNAME = "name"

const GROUP_TIMESERIES = "timeseries"
const DATA_TIMESERIES_TIME = "time"
const DATA_TIMESERIES_STEPS = "steps"
const DATA_TIMESERIES_VALUE = "value"

include("hdf5.jl")
include("logfile.jl")
include("metadata.jl")
include("timeseries.jl")

end # module


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

name = "opusfile"
version = v"0.11.0"

# Collection of sources required to complete build
sources = [
    "https://downloads.xiph.org/releases/opus/opusfile-0.11.tar.gz" =>
    "74ce9b6cf4da103133e7b5c95df810ceb7195471e1162ed57af415fabf5603bf",
    "./bundled"
]

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

if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 "${WORKSPACE}/srcdir/patches/no-openssl-wincert.patch"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --with-pic
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(!Sys.isfreebsd, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libopusurl", :libopusurl),
    LibraryProduct("libopusfile", :libopusfile)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    PackageSpec(name="Ogg_jll", uuid="e7412a2a-1a6e-54c0-be00-318e2571c051")
    PackageSpec(name="Opus_jll", uuid="91d4177d-7536-5919-b921-800302f37372")
    PackageSpec(name="OpenSSL_jll", uuid="458c3c95-2e84-50aa-8efc-19380b2a3a95")
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
export Space, domainspace, rangespace, maxspace,Space,conversion_type

##
# "enum" types used for whether the basis is
# Real or Complex.  AnyBasis is used when the answer
# is unknown.
##

immutable RealBasis end
immutable ComplexBasis end
immutable AnyBasis end

Base.promote_rule(::Type{RealBasis},::Type{ComplexBasis})=ComplexBasis
Base.promote_rule(::Type{ComplexBasis},::Type{AnyBasis})=AnyBasis
Base.promote_rule(::Type{RealBasis},::Type{AnyBasis})=AnyBasis

# coefficient_type(basis,valuetype) gives the type for coefficients
# for basis of type RealBasis/ComplexBasis and valuetype
# giving the type of function values
coefficient_type{T<:Complex}(::Type{ComplexBasis},::Type{T})=T
coefficient_type{T<:Real}(::Type{ComplexBasis},::Type{T})=Complex{T}
coefficient_type{T}(::Type{RealBasis},::Type{T})=T

#
# eltype for RealBasis/ComplexBasis gives the
# default type.  Maybe should be defaulteltype?
#

Base.eltype(::RealBasis)=Float64
Base.eltype(::ComplexBasis)=Complex{Float64}
Base.eltype(::AnyBasis)=Number

Base.eltype(::Type{RealBasis})=Float64
Base.eltype(::Type{ComplexBasis})=Complex{Float64}
Base.eltype(::Type{AnyBasis})=Number

# T is either RealBasis (cos/sin/polynomial) or ComplexBasis (laurent)
# D is the domain
# d is the dimension
abstract Space{T,D,d}

typealias RealSpace{D,d} Space{RealBasis,D,d}
typealias ComplexSpace{D,d} Space{ComplexBasis,D,d}
typealias UnivariateSpace{T,D} Space{T,D,1}
typealias BivariateSpace{T} Space{T,AnyDomain,2}
typealias RealUnivariateSpace{D} RealSpace{D,1}

Base.eltype{T}(::Space{T}) = eltype(T)
Base.eltype{T,D,d}(::Type{Space{T,D,d}}) = eltype(T)
basistype{T}(::Space{T}) = T
basistype{T,D,d}(::Type{Space{T,D,d}}) = T
basistype{FT<:Space}(::Type{FT}) = basistype(supertype(FT))

domaintype{T,D}(::Space{T,D}) = D
domaintype{T,D,d}(::Type{Space{T,D,d}}) = D
domaintype{FT<:Space}(::Type{FT}) = domaintype(supertype(FT))

coefficient_type{S}(::Space{S},T) = coefficient_type(S,T)

domaindimension{S,D,d}(::Space{S,D,d}) = d
dimension(::Space) = ∞  # We assume infinite-dimensional spaces

Space{D<:Number}(d::AbstractVector{D}) = Space(convert(Domain,d))

abstract AmbiguousSpace <: Space{RealBasis,AnyDomain,1}

domain(::AmbiguousSpace) = AnyDomain()

function setdomain{T,D<:Domain}(sp::Space{T,D},d::D)
    S=typeof(sp)
    @assert length(fieldnames(S))==1
    S(d)
end

# function setdomain(sp::Space,d::Domain)
#     S=typeof(sp)
#     @assert length(fieldnames(S))==1
#     # the domain is not compatible, but maybe we c
#     # can drop the space depence.  For example,
#     # CosSpace{Circle{Float64}} -> CosSpace
#     eval(parse(string(S.name.module)*"."*string(S.name)))(d)
# end

setcanonicaldomain(s) = setdomain(s,canonicaldomain(s))
reverseorientation(S::Space) = setdomain(S,reverse(domain(S)))

# UnsetSpace dictates that an operator is not defined until
#   its domainspace is promoted
# NoSpace is used to indicate no space exists for, e.g.,
# conversion_type

immutable UnsetSpace <: AmbiguousSpace end
immutable NoSpace <: AmbiguousSpace end
immutable ZeroSpace <: AmbiguousSpace end   # ZeroSpace is compatible with all spaces

dimension(::ZeroSpace) = 0

isambiguous(::)=false
isambiguous(::AmbiguousSpace) = true

#TODO: should it default to canonicalspace?
points(d::Space,n) = points(domain(d),n)

canonicalspace(T) = T
canonicaldomain(S::Space) = canonicaldomain(domain(S))

# Check whether spaces are the same, override when you need to check parameters
# This is used in place of == to support AnyDomain
spacescompatible{D<:Space}(f::D,g::D) = error("Override spacescompatible for "*string(D))
spacescompatible(::UnsetSpace,::UnsetSpace) = true
spacescompatible(::NoSpace,::NoSpace) = true
spacescompatible(::ZeroSpace,::ZeroSpace) = true
spacescompatible(f,g) = false
==(A::Space,B::Space) = spacescompatible(A,B)&&domain(A)==domain(B)
spacesequal(A::Space,B::Space) = A==B

# check a list of spaces for compatibility
for OP in (:spacescompatible,:domainscompatible,:spacesequal)
    @eval begin
        function $OP{T<:Space}(v::Vector{T})
            for k=1:length(v)-1
                if !$OP(v[k],v[k+1])
                    return false
                end
            end
            true
        end

        $OP{T<:Space}(v::Array{T}) = $OP(vec(v))
    end
end

domain(A::Space) = A.domain # assume it has a field domain

for op in (:tocanonical,:fromcanonical,:tocanonicalD,:fromcanonicalD,:invfromcanonicalD)
    @eval ($op)(sp::Space,x...)=$op(domain(sp),x...)
end

mappoint(a::Space,b::Space,x)=mappoint(domain(a),domain(b),x)
mappoint(a::Space,b::Domain,x)=mappoint(domain(a),b,x)
mappoint(a::Domain,b::Space,x)=mappoint(a,domain(b),x)

for FUNC in (:conversion_rule,:maxspace_rule,:union_rule)
    @eval begin
        function $FUNC(a,b)
            if spacescompatible(a,b)
                a
            else
                NoSpace()
            end
        end
    end
end

for FUNC in (:conversion_type,:maxspace)
    @eval begin
        $FUNC(::ZeroSpace,::UnsetSpace) = UnsetSpace()
        $FUNC(::UnsetSpace,::ZeroSpace) = UnsetSpace()
        $FUNC(a::UnsetSpace,b::UnsetSpace) = a
        $FUNC(a::UnsetSpace,b::Space) = b
        $FUNC(a::Space,b::UnsetSpace) = a
    end
end

# gives a space c that has a banded conversion operator TO a and b
function conversion_type(a,b)
    if spacescompatible(a,b)
        a
    elseif !domainscompatible(a,b)
        NoSpace()  # this avoids having to check eachtime
    else
        cr=conversion_rule(a,b)
        cr==NoSpace()?conversion_rule(b,a):cr
    end
end

# gives a space c that has a banded conversion operator FROM a and b
maxspace(a,b) = NoSpace()  # TODO: this fixes weird bug with Nothing
function maxspace(a::Space,b::Space)
    if spacescompatible(a,b)
        return a
    elseif !domainscompatible(a,b)
        return NoSpace()  # this avoids having to check eachtime
    end

    cr=maxspace_rule(a,b)
    if !isa(cr,NoSpace)
        return cr
    end

    cr=maxspace_rule(b,a)
    if !isa(cr,NoSpace)
        return cr
    end

    cr=conversion_type(a,b)
    if cr==a
        return b
    elseif cr ==b
        return a
    end

    # check if its banded through canonicalspace
    cspa=canonicalspace(a)
    if spacescompatible(cspa,b)
        # we can't call maxspace(cspa,a)
        # maxspace/conversion_type should be implemented for canonicalspace
        error("Override conversion_type or maxspace for "*string(a)*" and "*string(b))
    end
    if cspa != a && maxspace(cspa,a)==cspa
        return maxspace(b,cspa)
    end

    cspb=canonicalspace(b)
    if spacescompatible(cspb,a)
        # we can't call maxspace(cspb,b)
        error("Override conversion_type or maxspace for "*string(a)*" and "*string(b))
    end
    if cspb !=b && maxspace(cspb,b)==cspb
        return maxspace(a,cspb)
    end

    NoSpace()
end

# union combines two spaces
# this is used primarily for addition of two funs
# that may be incompatible
Base.union(a::AmbiguousSpace,b::AmbiguousSpace)=b
Base.union(a::AmbiguousSpace,b::Space)=b
Base.union(a::Space,b::AmbiguousSpace)=a
function Base.union(a::Space,b::Space)
    if spacescompatible(a,b)
        if isambiguous(domain(a))
            return b
        else
            return a
        end
    end

    cr=union_rule(a,b)
    if !isa(cr,NoSpace)
        return cr
    end

    cr=union_rule(b,a)
    if !isa(cr,NoSpace)
        return cr
    end

    cspa=canonicalspace(a)
    cspb=canonicalspace(b)
    if cspa!=a || cspb!=b
        cr=union(cspa,cspb)
    end
    if !isa(cr,NoSpace)
        return cr
    end

    #TODO: Uncomment when Julia bug is fixed
    # cr=maxspace(a,b)  #Max space since we can convert both to it
    # if !isa(cr,NoSpace)
    #     return cr
    # end

    a⊕b
end

Base.union(a::Space) = a

Base.union(a::Space,b::Space,c::Space) = union(union(a,b),c)
Base.union(a::Space,b::Space,c::Space,d::Space...) =
    union(union(a,b),c,d...)

# tests whether a Conversion operator exists
hasconversion(a,b) = maxspace(a,b)==b

# tests whether a coefficients can be converted to b
isconvertible(a,b) = hasconversion(a,b)

## Conversion routines
#       coefficients(v::Vector,a,b)
# converts from space a to space b
#       coefficients(v::Fun,a)
# is equivalent to coefficients(v.coefficients,v.space,a)
#       coefficients(v::Vector,a,b,c)
# uses an intermediate space b

coefficients(f,sp1,sp2,sp3) = coefficients(coefficients(f,sp1,sp2),sp2,sp3)

coefficients{T1<:Space,T2<:Space}(f::Vector,::Type{T1},::Type{T2}) =
    coefficients(f,T1(),T2())
coefficients{T1<:Space}(f::Vector,::Type{T1},sp2::Space) = coefficients(f,T1(),sp2)
coefficients{T2<:Space}(f::Vector,sp1::Space,::Type{T2}) = coefficients(f,sp1,T2())

## coefficients defaults to calling Conversion, otherwise it tries to pipe through Chebyshev

function defaultcoefficients(f,a,b)
    ct=conversion_type(a,b) # gives a space that has a banded conversion to both a and b

    if spacescompatible(a,b)
        f
    elseif spacescompatible(ct,a)
        (Conversion(a,b)*f).coefficients
    elseif spacescompatible(ct,b)
        (Conversion(b,a)\f).coefficients
    else
        csp=canonicalspace(a)

        if spacescompatible(a,csp)# a is csp, so try b
            csp=canonicalspace(b)
        end
        if spacescompatible(a,csp)||spacescompatible(b,csp)
            # b is csp too, so we are stuck, try Fun constructor
            if domain(b)⊆domain(a)
                coefficients(Fun(x->Fun(f,a)(x),b))
            else
                # we set the value to be zero off the domain of definition
                # but first ensure that domain(b) has a jump
                # TODO: this is disabled as it breaks the case of splitting
                #       one interval into two
                d=domain(a)
#                 if !issubcomponent(d,domain(b))
#                     error("$(d) is not a subcomponent of $(domain(b))")
#                 end

                coefficients(Fun(x->x∈d?Fun(f,a)(x):zero(Fun(f,a)(x)),b))
            end
        else
            coefficients(f,a,csp,b)
        end
    end
end

coefficients(f,a,b) = defaultcoefficients(f,a,b)

## TODO: remove zeros
Base.zero(S::Space) = zeros(S)
Base.zero{T<:Number}(::Type{T},S::Space) = zeros(T,S)
Base.zeros{T<:Number}(::Type{T},S::Space) = Fun(zeros(T,1),S)
Base.zeros(S::Space) = Fun(zeros(1),S)

# catch all
Base.ones(S::Space) = Fun(x->1.0,S)
Base.ones{T<:Number}(::Type{T},S::Space) = Fun(x->one(T),S)

identity_fun(S::Space) = identity_fun(domain(S))

function identity_fun(d::Domain)
    cd=canonicaldomain(d)
    if typeof(d)==typeof(cd)
        Fun(x->x,d) # fall back to constructor
    else
        # this allows support for singularities, that the constructor doesn't
        sf=fromcanonical(d,Fun(identity,cd))
        Fun(coefficients(sf),setdomain(space(sf),d))
    end
end

## rand
# checkpoints is used to give a list of points to double check
# the expansion
Base.rand(d::Space,k...)=rand(domain(d),k...)
checkpoints(d::Space)=checkpoints(domain(d))

## default transforms

# transform converts from values at points(S,n) to coefficients
# itransform converts from coefficients to values at points(S,n)

transform(S::Space,vals)=transform(S,vals,plan_transform(S,vals))
itransform(S::Space,cfs)=itransform(S,cfs,plan_itransform(S,cfs))

function transform(S::Space,vals,plan)
    csp=canonicalspace(S)
    if S==csp
        error("Override transform(::"*string(typeof(S))*",vals,plan)")
    end

    coefficients(transform(csp,vals,plan),csp,S)
end

function itransform(S::Space,cfs,plan)
    csp=canonicalspace(S)
    if S==csp
        error("Override itransform(::"*string(typeof(S))*",cfs,plan)")
    end

    itransform(csp,coefficients(cfs,S,csp),plan)
end

for OP in (:plan_transform,:plan_itransform)
    # plan transform expects a vector
    # this passes an empty Float64 array
    @eval $OP(S::Space,n::Integer)=$OP(S,Array(Float64,n))
end

function plan_transform(S::Space,vals)
    csp=canonicalspace(S)
    if S==csp
        identity #TODO: why identity?
    else
        plan_transform(csp,vals)
    end
end

function plan_itransform(S::Space,cfs)
    csp=canonicalspace(S)
    if S==csp
        identity #TODO: why identity?
    else
        plan_itransform(csp,cfs)
    end
end

## sorting
# we sort spaces lexigraphically by default

for OP in (:<,:(<=),:>,:(>=),:(Base.isless))
    @eval $OP(a::Space,b::Space)=$OP(string(a),string(b))
end

## Important special spaces

"""
`ConstantSpace` Represents a single number.  The remaining
coefficients are ignored.
"""

immutable ConstantSpace{DD} <: UnivariateSpace{RealBasis,DD}
    domain::DD
    ConstantSpace(d::DD)=new(d)
    ConstantSpace(d::AnyDomain)=new(DD(d))
end

ConstantSpace(d::Domain) = ConstantSpace{typeof(d)}(d)
ConstantSpace() = ConstantSpace(AnyDomain())

isconstspace(::ConstantSpace) = true

# we override maxspace instead of maxspace_rule to avoid
# domainscompatible check.
for OP in (:maxspace,:(Base.union))
    @eval begin
        $OP(A::ConstantSpace{AnyDomain},B::ConstantSpace{AnyDomain})=A
        $OP(A::ConstantSpace{AnyDomain},B::ConstantSpace)=B
        $OP(A::ConstantSpace,B::ConstantSpace{AnyDomain})=A
        $OP(A::ConstantSpace,B::ConstantSpace)=ConstantSpace(domain(A) ∪ domain(B))
    end
end

"""
`SequenceSpace` is the space of all sequences, i.e., infinite vectors
"""
immutable SequenceSpace <: Space{RealBasis,PositiveIntegers,0} end
const ℓ⁰ = SequenceSpace()
dimension(::SequenceSpace) = ∞
domain(::SequenceSpace) = ℕ
spacescompatible(::SequenceSpace,::SequenceSpace) = true


text:
decompose_trace(opt::BlackBoxOptim.OptRunController) = BlackBoxOptim.best_candidate(opt)

export BBO

struct BBO
    method::Symbol
    BBO(method) = new(method)
end

BBO() = BBO(:adaptive_de_rand_1_bin_radiuslimited) # the recommended optimizer as default

function __solve(prob::OptimizationProblem, opt::BBO, data = DEFAULT_DATA;
                 cb = (args...) -> (false), maxiters = nothing,
                 progress = false, kwargs...)

    local x, cur, state

    if data != DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
      cb_call = cb(decompose_trace(trace),x...)
      if !(typeof(cb_call) <: Bool)
        error("The callback should return a boolean `halt` for whether to stop the optimization process.")
      end
      if cb_call == true
        BlackBoxOptim.shutdown_optimizer!(trace) #doesn't work
      end
      cur, state = iterate(data, state)
      cb_call
    end

    if !(isnothing(maxiters)) && maxiters <= 0.0
        error("The number of maxiters has to be a non-negative and non-zero number.")
    elseif !(isnothing(maxiters))
        maxiters = convert(Int, maxiters)
    end

    _loss = function(θ)
        x = prob.f(θ, prob.p, cur...)
        return first(x)
    end

    bboptre = !(isnothing(maxiters)) ? BlackBoxOptim.bboptimize(_loss;Method = opt.method, SearchRange = [(prob.lb[i], prob.ub[i]) for i in 1:length(prob.lb)], MaxSteps = maxiters, CallbackFunction = _cb, CallbackInterval = 0.0, kwargs...) : BlackBoxOptim.bboptimize(_loss;Method = opt.method, SearchRange = [(prob.lb[i], prob.ub[i]) for i in 1:length(prob.lb)], CallbackFunction = _cb, CallbackInterval = 0.0, kwargs...)

    SciMLBase.build_solution(prob, opt, BlackBoxOptim.best_candidate(bboptre),
                             BlackBoxOptim.best_fitness(bboptre); original=bboptre)
end


text:
#####################
# WebhookEvent Type #
#####################

type WebhookEvent
    kind::GitHubString
    payload::Dict
    repository::Repo
    sender::Owner
end

function event_from_payload!(kind, data::Dict)
    if haskey(data, "repository")
        repository = Repo(data["repository"])
    elseif kind == "membership"
        repository = Repo("")
    else
        error("event payload is missing repository field")
    end

    if haskey(data, "sender")
        sender = Owner(data["sender"])
    else
        error("event payload is missing sender")
    end

    return WebhookEvent(kind, data, repository, sender)
end

########################
# Validation Functions #
########################

has_event_header(request::HttpCommon.Request) = haskey(request.headers, "X-GitHub-Event")
event_header(request::HttpCommon.Request) = request.headers["X-GitHub-Event"]

has_sig_header(request::HttpCommon.Request) = haskey(request.headers, "X-Hub-Signature")
sig_header(request::HttpCommon.Request) = request.headers["X-Hub-Signature"]

function has_valid_secret(request::HttpCommon.Request, secret)
    if has_sig_header(request)
        secret_sha = "sha1="*bytes2hex(MbedTLS.digest(MbedTLS.MD_SHA1, request.data, secret))
        return sig_header(request) == secret_sha
    end
    return false
end

function is_valid_event(request::HttpCommon.Request, events)
    return (has_event_header(request) && in(event_header(request), events))
end

function from_valid_repo(event, repos)
    return (name(event.repository) == "" || in(name(event.repository), repos))
end

#################
# EventListener #
#################

immutable EventListener
    server::HttpServer.Server
    function EventListener(handle; auth::Authorization = AnonymousAuth(),
                           secret = nothing, events = nothing,
                           repos = nothing, forwards = nothing)
        if !(isa(forwards, Void))
            forwards = map(HttpCommon.URI, forwards)
        end

        if !(isa(repos, Void))
            repos = map(name, repos)
        end

        server = HttpServer.Server() do request, response
            try
                handle_event_request(request, handle; auth = auth,
                                     secret = secret, events = events,
                                     repos = repos, forwards = forwards)
            catch err
                println("SERVER ERROR: $err")
                return HttpCommon.Response(500)
            end
        end

        server.http.events["listen"] = port -> begin
            println("Listening for GitHub events sent to $port;")
            println("Whitelisted events: $(isa(events, Void) ? "All" : events)")
            println("Whitelisted repos: $(isa(repos, Void) ? "All" : repos)")
        end

        return new(server)
    end
end

function handle_event_request(request, handle;
                              auth::Authorization = AnonymousAuth(),
                              secret = nothing, events = nothing,
                              repos = nothing, forwards = nothing)
    if !(isa(secret, Void)) && !(has_valid_secret(request, secret))
        return HttpCommon.Response(400, "invalid signature")
    end

    if !(isa(events, Void)) && !(is_valid_event(request, events))
        return HttpCommon.Response(400, "invalid event")
    end

    event = event_from_payload!(event_header(request), Requests.json(request))

    if !(isa(repos, Void)) && !(from_valid_repo(event, repos))
        return HttpCommon.Response(400, "invalid repo")
    end

    if !(isa(forwards, Void))
        for address in forwards
            Requests.post(address, request)
        end
    end

    return handle(event)
end

function Base.run(listener::EventListener, args...; kwargs...)
    return HttpServer.run(listener.server, args...; kwargs...)
end

###################
# CommentListener #
###################

const COMMENT_EVENTS = ["commit_comment",
                        "pull_request",
                        "pull_request_review_comment",
                        "issues",
                        "issue_comment"]

immutable CommentListener
    listener::EventListener
    function CommentListener(handle, trigger::Regex;
                             auth::Authorization = AnonymousAuth(),
                             check_collab::Bool = true,
                             secret = nothing,
                             repos = nothing,
                             forwards = nothing)
        listener = EventListener(auth=auth, secret=secret,
                                 events=COMMENT_EVENTS, repos=repos,
                                 forwards=forwards) do event
            return handle_comment(handle, event, auth, trigger, check_collab)
        end
        return new(listener)
    end
end

function Base.run(listener::CommentListener, args...; kwargs...)
    return run(listener.listener, args...; kwargs...)
end

function handle_comment(handle, event::WebhookEvent, auth::Authorization,
                        trigger::Regex, check_collab::Bool)
    kind, payload = event.kind, event.payload

    if (kind == "pull_request" || kind == "issues") && payload["action"] == "opened"
        body_container = kind == "issues" ? payload["issue"] : payload["pull_request"]
    elseif haskey(payload, "comment")
        body_container = payload["comment"]
    else
        return HttpCommon.Response(204, "payload does not contain comment")
    end

    if check_collab
        repo = event.repository
        user = body_container["user"]["login"]
        if !(iscollaborator(repo, user; auth = auth))
            return HttpCommon.Response(204, "commenter is not collaborator")
        end
    end

    trigger_match = match(trigger, body_container["body"])

    if trigger_match == nothing
        return HttpCommon.Response(204, "trigger match not found")
    end

    return handle(event, trigger_match)
end