text:
@test HPack.huffman_encode_bytes(b"www.example.com") ==
    [0xf1; 0xe3; 0xc2; 0xe5; 0xf2; 0x3a; 0x6b; 0xa0; 0xab; 0x90; 0xf4; 0xff]

@test HPack.huffman_decode_bytes(HPack.huffman_encode_bytes(b"www.example.com")) == b"www.example.com"

@test HPack.decode(HPack.new_dynamic_table(), IOBuffer([0x82; 0x84])) ==
    Headers(":method" => "GET",
            ":path" => "/")

## Below are examples copied from Appendix C

### C.2.1

headers = Headers("custom-key" => "custom-header")
headers_raw =
    [0x00; 0x0a; 0x63; 0x75; 0x73; 0x74; 0x6f; 0x6d; 0x2d; 0x6b; 0x65; 0x79;
     0x0d; 0x63; 0x75; 0x73; 0x74; 0x6f; 0x6d; 0x2d; 0x68; 0x65; 0x61; 0x64;
     0x65; 0x72]

@test HPack.encode(HPack.new_dynamic_table(), headers; huffman=false) ==
    headers_raw

@test HPack.decode(HPack.new_dynamic_table(), IOBuffer(headers_raw)) == headers

### C.2.2

headers = Headers(":path" => "/sample/path")
headers_raw =
    [0x04; 0x0c; 0x2f; 0x73; 0x61; 0x6d; 0x70; 0x6c; 0x65; 0x2f; 0x70; 0x61;
     0x74; 0x68]

@test HPack.decode(HPack.new_dynamic_table(), IOBuffer(HPack.encode(HPack.new_dynamic_table(), headers; huffman=true))) == headers

@test HPack.decode(HPack.new_dynamic_table(), IOBuffer(headers_raw)) == headers

### Server and client examples

request_headers = Headers(":method" => "GET",
                          ":path" => "/",
                          ":scheme" => "http",
                          ":authority" => "127.0.0.1:9000",
                          "accept" => "*/*",
                          "accept-encoding" => "gzip, deflate",
                          "user-agent" => "HTTP2.jl")

response_headers = Headers(":status" => "404",
                           "server" => "nghttpd nghttpd2/1.10.0",
                           "date" => "Thu, 02 Jun 2016 19:00:13 GMT",
                           "content-type" => "text/html; charset=UTF-8")

client_dt = HPack.new_dynamic_table()
server_dt = HPack.new_dynamic_table()

request_headers_raw = HPack.encode(client_dt, request_headers; huffman=true)
@test HPack.decode(server_dt, IOBuffer(request_headers_raw)) == request_headers

response_headers_raw = HPack.encode(server_dt, response_headers; huffman=true)
@test HPack.decode(client_dt, IOBuffer(response_headers_raw)) == response_headers

#### Recorded nghttp server response

response_headers_raw = [0x8d,0x76,0x90,0xaa,0x69,0xd2,0x9a,0xe4,0x52,0xa9,0xa7,0x4a,0x6b,0x13,0x01,0x5c,0x20,0x5c,0x1f,0x61,0x96,0xdf,0x3d,0xbf,0x4a,0x00,0x4a,0x65,0xb6,0xa5,0x04,0x00,0xb8,0xa0,0x5f,0xb8,0x26,0xee,0x32,0xda,0x98,0xb4,0x6f,0x5f,0x92,0x49,0x7c,0xa5,0x89,0xd3,0x4d,0x1f,0x6a,0x12,0x71,0xd8,0x82,0xa6,0x0e,0x1b,0xf0,0xac,0xf7]

client_dt = HPack.new_dynamic_table()
HPack.encode(client_dt, request_headers; huffman=false)
response_headers = HPack.decode(client_dt, IOBuffer(response_headers_raw))
println("Headers")
println("======================")
for header in response_headers
    println(header)
end


text:
abstract type AbstractPeriodicLattice{T,N_ARR,N_CDT} <: AbstractLattice{T,N_ARR,N_CDT} end

@doc """
A Lattice of points with `extent` describing the length along each dimension and `n_points` describing the number of points representing each dimension.
"""
struct PeriodicLattice{T,N_ARR} <: AbstractPeriodicLattice{T,N_ARR,N_ARR}
    arr::Array{NTuple{N_ARR,T},N_ARR}
    extent::NTuple{N_ARR,T}
end
(t::Type{<:PeriodicLattice})(; extent, n_points) =_periodiclattice(t, extent, n_points)
function _periodiclattice(t, extent::Number, n_points::Number)
    _periodiclattice(t, (extent,), (n_points,))
end
function _periodiclattice(t, extent::Tuple, n_points::Tuple)
    t(discrete_lattice(.-extent ./ 2, extent ./ 2 .- (extent ./ n_points), n_points), extent)
end
extent(p::PeriodicLattice) = p.extent

difference(p_lattice::PeriodicLattice, edge) = abs_difference_periodic(edge, extent(p_lattice))

# @recipe function f(lattice::PeriodicLattice{T,1}, values; val_lim=(0.0, 1.0)) where T
#     θ = coordinate_axes(lattice)[1] .* (2π / extent(lattice)[1])
#     #y = values .* sin.(θ)
#     #x = values .* cos.(θ)
#     seriestype := :scatter
#     projection := :polar
#     min_val, max_val = val_lim
#     min_r = max(abs(min_val), abs(max_val))
#     inner_r = min_r * 2
#     r_values = values .+ inner_r
#     #lims := (0, inner_r + max_val)

#     aspect_ratio := :equal

#     (θ |> collect, r_values)
# end

# @recipe function f(lattice::PeriodicLattice{T,2}, values::Array{T,2}; val_lim=nothing) where T
#     (x, y) = coordinate_axes(lattice) .|> collect
#     seriestype := :heatmap
#     if val_lim != nothing
#
#     (x,y,values)
# end

# TODO: Should be completely uniform in periodic case?
export simpson_weights
function simpson_weights(lattice::PeriodicLattice{T}) where T
    w = ones(T, size(lattice)...)
    return w
end
# function simpson_weights(lattice::PeriodicLattice{T,2}) where T
#     # http://mathfaculty.fullerton.edu/mathews/n2003/SimpsonsRule2DMod.html
#     @assert all(size(lattice) .% 2 .== 0)
#     w = ones(T, size(lattice)...)
#     w[1:2:end-1,:] .*= 4.0
#     w[2:2:end,:] .*= 2.0
#     w[:,1:2:end-1] .*= 4.0
#     w[:,2:2:end] .*= 2.0
#     return w
# end

ext: jl


text:
module FEMMaterials

using Materials, JuliaFEM, FEMBase, LinearAlgebra, SparseArrays, Tensors

# Material simulator to solve global system and run standard one element tests
include("mecamatso.jl")
export get_one_element_material_analysis, AxialStrainLoading, ShearStrainLoading
export update_bc_elements!, Continuum3D, MecaMatSo

# Material point simulator to study material behavior in single integration point
include("simulator.jl")
export Simulator

# tovoigt method for arrays
Tensors.tovoigt(x::Array{T, 2}) where T = [x[1,1], x[2,2], x[3,3], x[2,3], x[1,3], x[1,2]]

material_preprocess_increment!(material::M, element, ip, time) where {M<:AbstractMaterial} = nothing
material_postprocess_analysis!(material::M, element, ip, time) where {M<:AbstractMaterial} = nothing
material_postprocess_increment!(material::M, element, ip, time) where {M<:AbstractMaterial} = nothing
material_postprocess_iteration!(material::M, element, ip, time) where {M<:AbstractMaterial} = nothing

material_preprocess_analysis!(material::Nothing, element, ip, time) = nothing
material_preprocess_increment!(material::Nothing, element, ip, time) = nothing
material_preprocess_iteration!(material::Nothing, element, ip, time) = nothing
material_postprocess_analysis!(material::Nothing, element, ip, time) = nothing
material_postprocess_increment!(material::Nothing, element, ip, time) = nothing
material_postprocess_iteration!(material::Nothing, element, ip, time) = nothing

function update_ip!(material::M, ip, time) where {M<:AbstractMaterial}
    variables = fieldnames(typeof(material.variables))
    for variable in variables
        update!(ip, String(variable), time => copy(getfield(material.variables, variable)))
    end

    drivers = fieldnames(typeof(material.drivers))
    for driver in drivers
        update!(ip, String(driver), time => copy(getfield(material.drivers, driver)))
    end
end

# Copying to ip's the duplicate data?
function material_preprocess_analysis!(material::M, element, ip, time) where {M<:AbstractMaterial}
    update_ip!(material, ip, time)
    # Read parameter values
    expr = [element(String(p), ip, time) for p in fieldnames(typeof(material.parameters))]
    material.parameters = typeof(material.parameters)(expr...)
end

function material_preprocess_iteration!(material::M, element, ip, time) where {M<:AbstractMaterial}
    gradu = element("displacement", ip, time, Val{:Grad})
    strain = SymmetricTensor{2,3, Float64}((i,j) -> 0.5*(gradu[i,j]+gradu[j,i]))
    # Check for vector or tensor form needs to be implemented later
    dstrain = strain - material.drivers.strain
    #@info("time = $time, dstrain = $dstrain")
    material.ddrivers.strain = dstrain
    return nothing
end

material_preprocess_analysis!(material::M, element::Element{Poi1}, ip, time) where {M<:AbstractMaterial} = nothing
material_postprocess_analysis!(material::M, element::Element{Poi1}, ip, time) where {M<:AbstractMaterial} = nothing
material_preprocess_increment!(material::M, element::Element{Poi1}, ip, time) where {M<:AbstractMaterial} = nothing
material_postprocess_increment!(material::M, element::Element{Poi1}, ip, time) where {M<:AbstractMaterial} = nothing
material_preprocess_iteration!(material::M, element::Element{Poi1}, ip, time) where {M<:AbstractMaterial} = nothing
material_postprocess_iteration!(material::M, element::Element{Poi1}, ip, time) where {M<:AbstractMaterial} = nothing

function material_preprocess_increment!(material::M, element, ip, time) where {M<:AbstractMaterial}
    # Update time increment
    dtime = time - material.drivers.time
    material.ddrivers.time = dtime
    # Update parameters
    expr = [element(String(p), ip, time) - getfield(material.parameters, p) for p in fieldnames(typeof(material.parameters))]
    material.dparameters = typeof(material.dparameters)(expr...)
    return nothing
end

function material_postprocess_increment!(material::M, element, ip, time) where {M<:AbstractMaterial}
    update_material!(material)
    # Store history data to integration points
    update_ip!(material, ip, time)
end

end


text:
# Implements a simple versino of the Kalman smoother

using DataStructures

function filter1Dl2!(
    β::AbstractVector{Float64},
    y::AbstractVector{Float64},
    λ::AbstractVector{Float64},
    from::Int,
    to::Int
)
    # initialise with improper prior
    β[from] = y[from]
    ψstack = Stack{Float64}()
    push!(ψstack, 1.0)

    @simd for t in from:(to - 1)
        ψ = top(ψstack)
        ω = λ[t] * ψ / (λ[t] + ψ)
        β[t + 1] = (y[t + 1] + ω * β[t]) / (1.0 + ω)
        push!(ψstack, 1.0 + ω)
    end
    pop!(ψstack)

    ϕ = 1.0
    η = y[to]
    @simd for t in reverse(from:(to - 1))
        ψ = pop!(ψstack)
        ω = ϕ * λ[t] / (ϕ + λ[t])
        β[t] = (ψ * β[t] + ω * η) / (ψ + ω)
        ϕ =  1.0 + ω
        η = (y[t] + ω * η) / ϕ
    end
end

function filter1Dl2(
    y::AbstractVector{Float64},
    λ::AbstractVector{Float64},
    from::Int,
    to::Int
)
    β = zeros(to - from + 1)
    filter1Dl2!(β, y, λ, from, to)
    return β
end

function filter1Dl2!(
    β::AbstractVector{Float64},
    y::AbstractVector{Float64},
    λ::AbstractVector{Float64}
)
    filter1Dl2!(β, y, λ, 1, length(y))
end

function filter1Dl2(
    y::AbstractVector{Float64},
    λ::AbstractVector{Float64}
)
    filter1Dl2(y, λ, 1, length(y))
end


text:
using DocumentDistances: calc_marginal, frequencies_from_counts
using DocumentDistances: make_vocabulary, update_vocabulary!, calculate
using SparseArrays: spzeros

@testset "SinkhornDocumentDistance" begin

@testset "frequencies_from_counts" begin
    fs = frequencies_from_counts(Int[1, 3, 0])
    @test length(fs) == 3
    @test fs[1] == 0.25
    @test fs[2] == 0.75
    @test fs[3] == 0.0
end # @testset "frequencies_from_counts"

@testset "frequencies_from_counts w sparse input" begin
    v = spzeros(3)
    v[1] = 4
    v[2] = 1
    fs = frequencies_from_counts(v)
    @test length(fs) == 3
    @test fs[1] == 0.80
    @test fs[2] == 0.20
    @test fs[3] == 0.00
end # @testset "frequencies_from_counts"

@testset "make_vocabulary and update_vocabulary!" begin
    # @repeat(Vector{String}, Vector{String}) do d1, d2
    for _ in 1:NumRepeatedTests
        d1 = String[randstring(rand(0:5)) for _ in 1:rand(1:20)]
        vocab = make_vocabulary(d1)
        @test length(vocab) == length(unique(d1))
        for w in d1
            @test haskey(vocab, w)
        end
        d2 = String[randstring(rand(0:10)) for _ in 1:rand(1:20)]
        update_vocabulary!(vocab, d2)
        @test length(vocab) == length(unique(vcat(d1, d2)))
        for w in vcat(d1, d2)
            @test haskey(vocab, w)
        end
    end
end # @testset "make_vocabulary and update_vocabulary!"

@testset "calc_marginal" begin
    m1 = calc_marginal(["a"], ["a", "b"])
    @test isa(m1, AbstractVector{Float64})
    @test m1[1] == 1.0
    @test m1[2] == 0.0

    m2 = calc_marginal(["a", "b"], ["a", "b"])
    @test m2[1] == 0.5
    @test m2[2] == 0.5

    m3 = calc_marginal(["c", "b", "b", "c", "c"], ["a", "b", "c"])
    @test m3[1] == 0.0
    @test m3[2] == (2.0/5)
    @test m3[3] == (3.0/5)

    m4 = calc_marginal(String[], ["a", "b", "c"])
    @test length(m4) == 3
end # @testset "calc_marginal"

@testset "calculate and evaluate on strings and docs" begin
    D1 = ["obama", "speaks", "media", "illinois"]
    D2 = ["president", "greets", "press", "chicago"]

    sdd = SinkhornDocumentDistance()
    cost12, plan12 = calculate(sdd, D1, D2)
    d12 = evaluate(sdd, D1, D2)
    @test d12 == cost12
    @test size(plan12, 1) == 8
    @test size(plan12, 2) == 8

    s1 = join(D1, " ")
    s2 = join(D2, " ")
    d12b = evaluate(sdd, s1, s2)
    @test d12b == d12
end

end # @testset "SinkhornDocumentDistance"


text:
using MatrixProfile
using Test, Statistics, LinearAlgebra, Plots

using SlidingDistancesBase
using MatrixProfile: Profile

normdist = (x,y)->norm(znorm(x)-znorm(y))

function naive_matrix_profile(A,B,m)
    res = map(1:length(B)-m+1) do i
        findmin(distance_profile(ZEuclidean(), getwindow(B,m,i), A))
    end
    Profile(B,getindex.(res, 1), getindex.(res, 2), m, A)
end

function profeq(p1,p2)
    cond = p1.P ≈ p2.P
    cond && return true
    @show norm(p1.P - p2.P)
    false
end

@testset "MatrixProfile.jl" begin

   t = range(0, stop=1, step=1/10)
   y0 = sin.(2pi .* t)
   T = [randn(50); y0; randn(50); y0; randn(50)]
   A = sign.([randn(50); y0; randn(50); y0; randn(50)])[1:end÷2] .+ 0.01.*randn.()

   profile = @inferred matrix_profile(T, length(y0))
   P,I = profile.P, profile.I
   @test_nowarn plot(profile)
   # plot(T, layout=2)
   # plot!(P, sp=2)

   m = findmin(P)
   @test m[1] < 1e-6
   @test m[2] == 51 || m[2] == 112

   # Test Euclidean between two series
   profile3 = @inferred matrix_profile(T, T, length(y0))
   @test all(profile3.P .< 1e-6)

   profile5 = matrix_profile(T[1:end÷2], T, length(y0))
   @test sum(profile5.P[1:length(T)÷2-length(y0)]) < 1e-5

   profile6 = matrix_profile(T, T[1:end÷2], length(y0))
   @test sum(profile6.P) < 1e-5

   profile7n = naive_matrix_profile(T, A, length(y0))
   profile7n2 = @inferred matrix_profile(T, A, length(y0), normdist)
   profile7 = matrix_profile(T, A, length(y0))
   @test profeq(profile7, profile7n)
   @test profeq(profile7n, profile7n2)

   profile8n = naive_matrix_profile(A, T, length(y0))
   profile8n2 = @inferred matrix_profile(A, T, length(y0), normdist)
   profile8 = @inferred matrix_profile(A, T, length(y0))
   @test profeq(profile8, profile8n)
   @test profeq(profile8n, profile8n2)

   @testset "Generic matrix profile" begin
       @info "Testing Generic matrix profile"

       # Test generic one serie
       profile2 = @inferred matrix_profile(T, length(y0), normdist)
       @test profile2.P ≈  profile.P
       @test all(eachindex(profile.I)) do i
           profile2.I[i] ∈ (51,112) || profile2.I[i] == profile.I[i]
       end

       # Test generic between two series
       profile4 = @inferred matrix_profile(T, T, length(y0), normdist)
       @test all(profile4.P .< 1e-6)
   end

   @testset "Motifs" begin
       @info "Testing Motifs"
       t = range(0, stop=1, step=1/10)
       y0 = sin.(2pi .* t)
       T = [randn(50); y0; randn(50); y0; randn(50)]
       A = sign.([randn(50); y0; randn(50); y0; randn(50)])[1:end÷2] .+ 0.01.*randn.()
       profile = matrix_profile(T, length(y0))

       @time mot = motifs(profile, 2, r=2, th=5)
       @test MatrixProfile.subseqtype(mot) == "Motif"
       @test onsets(mot[1]) == [51, 112]
       @test_nowarn plot(mot)
       @test_nowarn plot(profile, mot)

       an = anomalies(profile, 3)
       @test length(an) == 3
       @test length(seqs(an)) == 3
       @test MatrixProfile.subseqtype(an) == "Anomaly"
       @test onsets(an)[1] == argmax(profile.P)
       @test_nowarn plot(an)
       @test MatrixProfile.seqlength(an) == length(y0)

   end

   @time matrix_profile(randn(Float32, 2^15), 256)

   @testset "mpdist" begin
       @info "Testing mpdist"

       t = 1:0.01:3
       A = @. sin(2pi * t) + 0.1 * randn()
       B = @. sign(sin(2pi * t)) + 0.1 * randn()
       m = 4
       # p1 = matrix_profile(A,B,m)
       # p2 = matrix_profile(B,A,m)
       # partialsort([p1.P; p2.P], 20)
       # plot(p1)
       # plot!(p2)
       @test @inferred(mpdist(A,B, m)) > 0
       @test mpdist(A,A, m) < 10sqrt(eps())
       @test mpdist(B,B, m) < 10sqrt(eps())
       T = [A; B]
       p = mpdist_profile(T, 50, 5)
       @test_nowarn plot(p)

       snips = snippets(T, 2, 50, m=5)
       @test_nowarn plot(snips)

   end

 end

# Benchmark
# @btime matrix_profile($(randn(Float32, 10_000)), 20)
#
# @time matrix_profile((randn(Float32, 100_000)),  20)
#
# @time matrix_profile((randn(Float32, 10_000)), (randn(Float32, 1_000)), 20, (x,y)->norm(x-y))
#
# @profiler matrix_profile(randn(10_000), 20)
# Q = randn(5)
# T = randn(10)
# d1 = DynamicAxisWarping.window_dot(Q,T)
# d2 = DynamicAxisWarping.window_dot2(Q,T)
# d3 = DynamicAxisWarping.window_dot3(Q,T)
# [d1 d2 d3]

# Eval precision
# a = randn(Float64, 50_000)
# p1 = matrix_profile(a,  20)
# p2 = matrix_profile(Float32.(a),  20)
# @show mean(abs.(p1.P - p2.P)./p1.P)
# plot(abs.(p1.P - p2.P)./p1.P)

# path = "/tmp/MixedBag/01911m_02019m_III_7680_200.txt"
# T = parse.(Int, split(join(Char.(read(path))), ','))
# profile, snips, Cfracs = snippets(T, 3, 20, 10)
# plot(plot(snips, layout = (1, 3), size = (800, 200)), plot(profile, snips, legend = false), layout=(2,1))


text:
############################################
################ Structures ################
############################################

function glb(types...)
    # If a lower bound is in the types then it is greatest
    # else we just return Unknown for now
    for type in types
        all(type <: t_ for t_ in types) && return type
    end
    return Unknown
end

function input_target_scitypes(models, metalearner)
    # The target scitype is defined as the greatest lower bound of the
    # metalearner and the base models in the library
    all_tg_scitypes = [target_scitype(m) for m in models]
    tg_scitype = glb(target_scitype(metalearner), all_tg_scitypes...)
    # The input scitype is defined as the greatest lower bound of the
    # base models in the library
    inp_scitype = glb([input_scitype(m) for m in models]...)

    return inp_scitype, tg_scitype
end

mutable struct DeterministicStack{modelnames, inp_scitype, tg_scitype} <: DeterministicComposite
   models::Vector{Supervised}
   metalearner::Deterministic
   resampling
   measures::Union{Nothing,AbstractVector}
   function DeterministicStack(modelnames, models, metalearner, resampling, measures)
        inp_scitype, tg_scitype = input_target_scitypes(models, metalearner)
        return new{modelnames, inp_scitype, tg_scitype}(models, metalearner, resampling, measures)
   end
end

mutable struct ProbabilisticStack{modelnames, inp_scitype, tg_scitype} <: ProbabilisticComposite
    models::Vector{Supervised}
    metalearner::Probabilistic
    resampling
    measures::Union{Nothing,AbstractVector}
    function ProbabilisticStack(modelnames, models, metalearner, resampling, measures)
        inp_scitype, tg_scitype = input_target_scitypes(models, metalearner)
        return new{modelnames, inp_scitype, tg_scitype}(models, metalearner, resampling, measures)
    end
 end

const Stack{modelnames, inp_scitype, tg_scitype} =
    Union{DeterministicStack{modelnames, inp_scitype, tg_scitype},
            ProbabilisticStack{modelnames, inp_scitype, tg_scitype}}

"""
    Stack(;metalearner=nothing, resampling=CV(), name1=model1, name2=model2, ...)

Implements the two-layer generalized stack algorithm introduced by
[Wolpert
(1992)](https://www.sciencedirect.com/science/article/abs/pii/S0893608005800231)
and generalized by [Van der Laan et al
(2007)](https://biostats.bepress.com/ucbbiostat/paper222/). Returns an
instance of type `ProbablisiticStack` or `DeterministicStack`,
depending on the prediction type of `metalearner`.

When training a machine bound to such an instance:

- The data is split into training/validation sets according to the
  specified `resampling` strategy.

- Each base model `model1`, `model2`, ... is trained on each training
  subset and outputs predictions on the corresponding validation
  sets. The multi-fold predictions are spliced together into a
  so-called out-of-sample prediction for each model.

- The adjudicating model, `metalearner`, is subsequently trained on
  the out-of-sample predictions to learn the best combination of base
  model predictions.

- Each base model is retrained on all supplied data for purposes of
  passing on new production data onto the adjudicator for making new
  predictions

### Arguments

- `metalearner::Supervised`: The model that will optimize the desired
  criterion based on its internals.  For instance, a LinearRegression
  model will optimize the squared error.

- `resampling`: The resampling strategy used
  to prepare out-of-sample predictions of the base learners.
  It can be a user-defined strategy, the only
  caveat being that it should have a `nfolds` attribute.

- `measures`: A measure or iterable over measures, to perform an internal
  evaluation of the learners in the Stack while training. This is not for the
  evaluation of the Stack itself.

- `name1=model1, name2=model2, ...`: the `Supervised` model instances
  to be used as base learners.  The provided names become properties
  of the instance created to allow hyper-parameter access

### Example

The following code defines a `DeterministicStack` instance for
learning a `Continuous` target, and demonstrates that:

- Base models can be `Probabilistic` models even if the stack
  itself is `Deterministic` (`predict_mean` is applied in such cases).

- As an alternative to hyperparameter optimization, one can stack
  multiple copies of given model, mutating the hyper-parameter used in
  each copy.

```julia
using MLJ

DecisionTreeRegressor = @load DecisionTreeRegressor pkg=DecisionTree
EvoTreeRegressor = @load EvoTreeRegressor
XGBoostRegressor = @load XGBoostRegressor
KNNRegressor = @load KNNRegressor pkg=NearestNeighborModels
LinearRegressor = @load LinearRegressor pkg=MLJLinearModels

X, y = make_regression(500, 5)

stack = Stack(;metalearner=LinearRegressor(),
                resampling=CV(),
                measures=rmse,
                constant=ConstantRegressor(),
                tree_2=DecisionTreeRegressor(max_depth=2),
                tree_3=DecisionTreeRegressor(max_depth=3),
                evo=EvoTreeRegressor(),
                knn=KNNRegressor(),
                xgb=XGBoostRegressor())

mach = machine(stack, X, y)
evaluate!(mach; resampling=Holdout(), measure=rmse)

```

The internal evaluation report can be accessed like this
and provides a PerformanceEvaluation object for each model:

```julia
report(mach).cv_report
```

"""
function Stack(;metalearner=nothing, resampling=CV(), measure=nothing, measures=measure, named_models...)
    metalearner === nothing &&
        throw(ArgumentError("No metalearner specified. Use Stack(metalearner=...)"))

    nt = NamedTuple(named_models)
    modelnames = keys(nt)
    models = collect(nt)
    if (measures !== nothing) && !(measures isa AbstractVector)
        measures = [measures, ]
    end

    if metalearner isa Deterministic
        stack =  DeterministicStack(modelnames, models, metalearner, resampling, measures)
    elseif metalearner isa Probabilistic
        stack = ProbabilisticStack(modelnames, models, metalearner, resampling, measures)
    else
        throw(ArgumentError("The metalearner should be a subtype
                    of $(Union{Deterministic, Probabilistic})"))
    end
    # Issuing clean! statement
    message = MMI.clean!(stack)
    isempty(message) || @warn message

    # Warning if either input_scitype/target_scitype is
    # Unknown at construction time
    params = typeof(stack).parameters
    params[end-1] == Unknown && @warn "Could not infer input_scitype of the stack"
    params[end] == Unknown && @warn "Could not infer target_scitype of the stack"

    return stack
end

function MMI.clean!(stack::Stack{modelnames, inp_scitype, tg_scitype}) where {modelnames,inp_scitype,tg_scitype}
    # We only carry checks and don't try to correct the arguments here
    message = ""
    # Checking target_scitype and input_scitype have not been changed from the original stack
    glb_inp_scitype, glb_tg_scitype = input_target_scitypes(getfield(stack, :models), stack.metalearner)
    glb_inp_scitype == inp_scitype ||
            throw(DomainError(inp_scitype, "The newly inferred input_scitype of the stack doesn't
            match its original one. You have probably changed one of the base models or the metalearner
            to a non compatible type."))
    glb_tg_scitype == tg_scitype ||
            throw(DomainError(tg_scitype, "The newly inferred target_scitype of the stack doesn't
            match its original one. You have probably changed one of the base model or the metalearner
            to a non compatible type."))
    # Checking the target scitype is consistent with either Probabilistic/Deterministic Stack
    target_scitype(stack.metalearner) <: Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}} ||
        throw(ArgumentError("The metalearner should have target_scitype:
                $(Union{AbstractArray{<:Continuous}, AbstractArray{<:Finite}})"))

    return message
end

Base.propertynames(::Stack{modelnames}) where modelnames = tuple(:resampling, :metalearner, modelnames...)

function Base.getproperty(stack::Stack{modelnames}, name::Symbol) where modelnames
    name === :metalearner && return getfield(stack, :metalearner)
    name === :resampling && return getfield(stack, :resampling)
    name == :measures && return getfield(stack, :measures)
    models = getfield(stack, :models)
    for j in eachindex(modelnames)
        name === modelnames[j] && return models[j]
    end
    error("type Stack has no property $name")
end

function Base.setproperty!(stack::Stack{modelnames}, _name::Symbol, val) where modelnames
    _name === :metalearner && return setfield!(stack, :metalearner, val)
    _name === :resampling && return setfield!(stack, :resampling, val)
    _name === :measures && return setfield!(stack, :measures, val)
    idx = findfirst(==(_name), modelnames)
    idx isa Nothing || return getfield(stack, :models)[idx] = val
    error("type Stack has no property $name")
end

MMI.target_scitype(::Type{<:Stack{modelnames, input_scitype, target_scitype}}) where
    {modelnames, input_scitype, target_scitype} = target_scitype

MMI.input_scitype(::Type{<:Stack{modelnames, input_scitype, target_scitype}}) where
    {modelnames, input_scitype, target_scitype} = input_scitype

MLJBase.load_path(::Type{<:ProbabilisticStack}) = "MLJBase.ProbabilisticStack"
MLJBase.load_path(::Type{<:DeterministicStack}) = "MLJBase.DeterministicStack"
MLJBase.package_name(::Type{<:Stack}) = "MLJBase"
MLJBase.package_uuid(::Type{<:Stack}) = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
MLJBase.package_url(::Type{<:Stack}) = "https://github.com/alan-turing-institute/MLJBase.jl"
MLJBase.package_license(::Type{<:Stack}) = "MIT"

###########################################################
################# Node operations Methods #################
###########################################################

getfolds(y::AbstractNode, cv::CV, n::Int) =
    source(train_test_pairs(cv, 1:n))

getfolds(y::AbstractNode, cv::StratifiedCV, n::Int) =
    node(YY->train_test_pairs(cv, 1:n, YY), y)

trainrows(X::AbstractNode, folds::AbstractNode, nfold) =
    node((XX, ff) -> selectrows(XX, ff[nfold][1]), X, folds)

testrows(X::AbstractNode, folds::AbstractNode, nfold) =
    node((XX, ff) -> selectrows(XX, ff[nfold][2]), X, folds)

pre_judge_transform(ŷ::Node, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Finite}}) =
    node(ŷ -> pdf(ŷ, levels(first(ŷ))), ŷ)

pre_judge_transform(ŷ::Node, ::Type{<:Probabilistic}, ::Type{<:AbstractArray{<:Continuous}}) =
    node(ŷ->mean.(ŷ), ŷ)

pre_judge_transform(ŷ::Node, ::Type{<:Deterministic}, ::Type{<:AbstractArray{<:Continuous}}) =
    ŷ

store_for_evaluation(mach::Machine, Xtest::AbstractNode, ytest::AbstractNode, measures::Nothing) = nothing
function store_for_evaluation(mach::Machine, Xtest::AbstractNode, ytest::AbstractNode, measures)
    node((ytest, Xtest) -> [mach, Xtest, ytest], ytest, Xtest)
end

"""
    internal_stack_report(m::Stack, verbosity::Int, y::AbstractNode, folds_evaluations::Vararg{Nothing})

When measure/measures is a Nothing, the folds_evaluation won't have been filled by `store_for_evaluation`
and we thus return an empty NamedTuple.
"""
internal_stack_report(m::Stack, verbosity::Int, y::AbstractNode, folds_evaluations::Vararg{Nothing}) = NamedTuple{}()

"""
    internal_stack_report(m::Stack, verbosity::Int, y::AbstractNode, folds_evaluations::Vararg{AbstractNode})

When measure/measures is provided, the folds_evaluation will have been filled by `store_for_evaluation`. This function is
not doing any heavy work (not constructing nodes corresponding to measures) but just unpacking all the folds_evaluations in a single node that
can be evaluated later.
"""
function internal_stack_report(m::Stack, verbosity::Int, y::AbstractNode, folds_evaluations::Vararg{AbstractNode})
    _internal_stack_report(y, folds_evaluations...) =
        internal_stack_report(m, verbosity, y, folds_evaluations...)
    return (report=(cv_report=node(_internal_stack_report, y, folds_evaluations...),),)
end

"""
    internal_stack_report(stack::Stack{modelnames,}, verbosity::Int, y, folds_evaluations...) where modelnames

Returns a `NamedTuple` of `PerformanceEvaluation` objects, one for each model. The folds_evaluations
are  built in a flatten array respecting the order given by:
(fold_1:(model_1:[mach, Xtest, ytest], model_2:[mach, Xtest, ytest], ...), fold_2:(model_1, model_2, ...), ...)
"""
function internal_stack_report(stack::Stack{modelnames,}, verbosity::Int, y, folds_evaluations...) where modelnames

    n_measures = length(stack.measures)
    nfolds = stack.resampling.nfolds

    # For each model we record the results mimicking the fields PerformanceEvaluation
    results = NamedTuple{modelnames}([
            (measure=stack.measures,
            measurement=Vector{Any}(undef, n_measures),
            operation=_actual_operations(nothing, stack.measures, model, verbosity),
            per_fold=[Vector{Any}(undef, nfolds) for _ in 1:n_measures],
            per_observation=Vector{Union{Missing, Vector{Any}}}(missing, n_measures),
            fitted_params_per_fold=[],
            report_per_fold=[],
            train_test_pairs=train_test_pairs(stack.resampling, 1:nrows(y), y)
            )
        for model in getfield(stack, :models)]
    )

    # Update the results
    index = 1
    for foldid in 1:nfolds
        for modelname in modelnames
            model_results = results[modelname]
            mach, Xtest, ytest = folds_evaluations[index]
            # Update report and fitted_params per fold
            push!(model_results.fitted_params_per_fold, fitted_params(mach))
            push!(model_results.report_per_fold, report(mach))
            # Loop over measures to update per_observation and per_fold
            for (i, (measure, operation)) in enumerate(zip(stack.measures, model_results.operation))
                ypred = operation(mach, Xtest)
                loss = measure(ypred, ytest)
                # Update per_observation
                if reports_each_observation(measure)
                    if model_results.per_observation[i] === missing
                        model_results.per_observation[i] = Vector{Any}(undef, nfolds)
                    end
                    model_results.per_observation[i][foldid] = loss
                end

                # Update per_fold
                model_results.per_fold[i][foldid] =
                    reports_each_observation(measure) ? MLJBase.aggregate(loss, measure) : loss
            end
            index += 1
        end
    end

    # Update measurement field by aggregation
    for modelname in modelnames
        for (i, measure) in enumerate(stack.measures)
            model_results = results[modelname]
            model_results.measurement[i] = MLJBase.aggregate(model_results.per_fold[i], measure)
        end
    end

    return NamedTuple{modelnames}([PerformanceEvaluation(r...) for r in results])
end

check_stack_measures(stack, verbosity::Int, measures::Nothing, y) = nothing
"""
    check_stack_measures(stack, measures, y)

Check the measures compatibility for each model in the Stack.
"""
function check_stack_measures(stack, verbosity::Int, measures, y)
    for model in getfield(stack, :models)
        operations = _actual_operations(nothing, measures, model, verbosity)
        _check_measures(measures, operations, model, y)
    end
end

"""
    oos_set(m::Stack, folds::AbstractNode, Xs::Source, ys::Source)

This function is building the out-of-sample dataset that is later used by the `judge`
for its own training. It also returns the folds_evaluations object if internal
cross-validation results are requested.
"""
function oos_set(m::Stack, folds::AbstractNode, Xs::Source, ys::Source)
    Zval = []
    yval = []
    folds_evaluations = []
    # Loop over the cross validation folds to build a training set for the metalearner.
    for nfold in 1:m.resampling.nfolds
        Xtrain = trainrows(Xs, folds, nfold)
        ytrain = trainrows(ys, folds, nfold)
        Xtest = testrows(Xs, folds, nfold)
        ytest = testrows(ys, folds, nfold)

        # Train each model on the train fold and predict on the validation fold
        # predictions are subsequently used as an input to the metalearner
        Zfold = []
        for model in getfield(m, :models)
            mach = machine(model, Xtrain, ytrain)
            ypred = predict(mach, Xtest)
            # Internal evaluation on the fold if required
            push!(folds_evaluations, store_for_evaluation(mach, Xtest, ytest, m.measures))
            # Dispatch the computation of the expected mean based on
            # the model type and target_scytype
            ypred = pre_judge_transform(ypred, typeof(model), target_scitype(model))
            push!(Zfold, ypred)
        end
        Zfold = hcat(Zfold...)

        push!(Zval, Zfold)
        push!(yval, ytest)
    end

    Zval = MLJBase.table(vcat(Zval...))
    yval = vcat(yval...)

    Zval, yval, folds_evaluations
end

#######################################
################# Fit #################
#######################################
"""
    fit(m::Stack, verbosity::Int, X, y)
"""
function fit(m::Stack, verbosity::Int, X, y)
    check_stack_measures(m, verbosity, m.measures, y)

    n = nrows(y)

    Xs = source(X)
    ys = source(y)

    folds = getfolds(ys, m.resampling, n)

    Zval, yval, folds_evaluations = oos_set(m, folds, Xs, ys)

    metamach = machine(m.metalearner, Zval, yval)

    # Each model is retrained on the original full training set
    Zpred = []
    for model in getfield(m, :models)
        mach = machine(model, Xs, ys)
        ypred = predict(mach, Xs)
        ypred = pre_judge_transform(ypred, typeof(model), target_scitype(model))
        push!(Zpred, ypred)
    end

    Zpred = MLJBase.table(hcat(Zpred...))
    ŷ = predict(metamach, Zpred)

    internal_report = internal_stack_report(m, verbosity, ys, folds_evaluations...)

    # We can infer the Surrogate by two calls to supertype
    mach = machine(supertype(supertype(typeof(m)))(), Xs, ys; predict=ŷ, internal_report...)

    return!(mach, m, verbosity)
end


text:
function load_offsets(obj::UTCTAIOffsetsLoader)
    return jcall(obj, "loadOffsets", List, ())
end


text:
# 1-dimensional Total Variation (times a constant)

export TotalVariation1D

"""
**1-dimensional Total Variation**

    TotalVariation1D(λ=1)

With a nonnegative scalar parameter λ, returns the function
```math
f(x) = λ ∑_{i=2}^{n} |x_i - x_{i-1}|.
```
"""
struct TotalVariation1D{T <: Real} <: ProximableFunction
    lambda::T
    function TotalVariation1D{T}(lambda::T) where {T <: Real}
        if lambda < 0
            error("parameter λ must be nonnegative")
        else
            new(lambda)
        end
    end
end

is_separable(f::TotalVariation1D) = false
is_convex(f::TotalVariation1D) = true
is_positively_homogeneous(f::TotalVariation1D) = true

TotalVariation1D(lambda::R=1) where {R <: Real} = TotalVariation1D{R}(lambda)

function (f::TotalVariation1D)(x::AbstractArray)
    return f.lambda * norm(x[2:end] - x[1:end-1], 1)
end

# Condat algorithm
# https://lcondat.github.io/publis/Condat-fast_TV-SPL-2013.pdf
function tvnorm_prox_condat(y::AbstractArray, x::AbstractArray, lambda::Real)
    # solves y = arg min_z lambda*sum_k |z_{k+1}-z_k| + 1/2 * ||z-x||^2
    N = length(x)

    k = k0 = kmin = kplus = 1
    vmin = x[1] - lambda
    vmax = x[1] + lambda
    umin = lambda
    umax = -lambda

    while 0 < 1
        while k == N
            if umin < 0
                y[k0:kmin] .= vmin
                kmin += 1
                k = k0 = kmin
                vmin = x[k]
                umin = lambda
                umax = x[k] + lambda - vmax
            elseif umax > 0
                y[k0:kplus] .= vmax
                kplus +=1
                k = k0 = kplus
                vmax = x[k]
                umax = -lambda
                umin = x[k] - lambda - vmin
            else
                y[k0:N] .= vmin + umin/(k-k0+1)
                return
            end

            if k==N
                y[N] = vmin + umin
                return
            end
        end

        if x[k+1] + umin < vmin - lambda
            y[k0:kmin] .= vmin
            kmin += 1
            k = k0 = kplus = kmin
            vmin = x[k]
            vmax = x[k] + 2*lambda
            umin = lambda
            umax = -lambda
        elseif x[k+1] + umax > vmax + lambda
            y[k0:kplus] .= vmax
            kplus += 1
            k = k0 =kmin = kplus
            vmin = x[k] - 2*lambda
            vmax = x[k]
            umin = lambda
            umax = -lambda
        else
            k += 1
            umin = umin + x[k] - vmin
            umax = umax + x[k] - vmax
            if umin >= lambda
                vmin = vmin + (umin-lambda)/(k-k0+1)
                umin = lambda
                kmin = k
            end

            if umax <= -lambda
                vmax += (umax+lambda)/(k-k0+1)
                umax = -lambda
                kplus = k
            end
        end
    end
end

function prox!(y::AbstractArray{T}, f::TotalVariation1D, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
    a = gamma * f.lambda
    tvnorm_prox_condat(y, x, a)
    return f.lambda * norm(y[2:end] - y[1:end-1], 1)
end

fun_name(f::TotalVariation1D) = "1D Total Variation"
fun_dom(f::TotalVariation1D) = "AbstractArray{Real}"
fun_expr(f::TotalVariation1D) = "x ↦ λ ∑_{i=2}^{n} |x_i - x_{i-1}|"
fun_params(f::TotalVariation1D) = "λ = $(f.lambda)"


text:
# to control chain initialization:
myinit(n, m) = reshape(convert(Vector{Float32}, (1:n*m)), n , m)

mutable struct TESTBuilder <: MLJFlux.Builder end
MLJFlux.build(builder::TESTBuilder, rng, n_in, n_out) =
    Flux.Chain(Flux.Dense(n_in, n_out, init=myinit))

@testset_accelerated "issue #152" accel begin

    # data:
    n = 100
    d = 5
    Xmat = rand(Float32, n, d)
#    Xmat = fill(one(Float32), n, d)
    X = MLJBase.table(Xmat);
    y = X.x1 .^2 + X.x2 .* X.x3 - 4 * X.x4

    # train a model on all the data using batch size > 1:
    model = MLJFlux.NeuralNetworkRegressor(builder=TESTBuilder(),
                                           batch_size=25,
                                           epochs=1,
                                           loss=Flux.mse,
                                           acceleration=accel)
    mach = machine(model, X, y)
    fit!(mach, verbosity=0)

    # extract the pre-training loss computed in the `fit!(chain, ...)` method:
    pretraining_loss = report(mach).training_losses[1]

    # compute by hand:
    chain0 = myinit(1, d)
    pretraining_yhat = Xmat*chain0' |> vec
    @test y isa Vector && pretraining_yhat isa Vector
    pretraining_loss_by_hand =  MLJBase.l2(pretraining_yhat, y) |> mean
    mean(((pretraining_yhat - y).^2)[1:2])

    # compare:
    @test pretraining_loss ≈ pretraining_loss_by_hand

end

@testset_accelerated "Short" accel begin
    builder = MLJFlux.Short(n_hidden=4, σ=Flux.relu, dropout=0)
    chain = MLJFlux.build(builder, StableRNGs.StableRNG(123), 5, 3)
    ps = Flux.params(chain)
    @test size.(ps) == [(4, 5), (4,), (3, 4), (3,)]

    # reproducibility (without dropout):
    chain2 = MLJFlux.build(builder, StableRNGs.StableRNG(123), 5, 3)
    x = rand(5)
    @test chain(x) ≈ chain2(x)
end

@testset_accelerated "@builder" accel begin
    builder = MLJFlux.@builder(Flux.Chain(Flux.Dense(n_in, 4,
                                                     init = (out, in) -> randn(rng, out, in)),
                                     Flux.Dense(4, n_out)))
    rng = StableRNGs.StableRNG(123)
    chain = MLJFlux.build(builder, rng, 5, 3)
    ps = Flux.params(chain)
    @test size.(ps) == [(4, 5), (4,), (3, 4), (3,)]

    chain2 = MLJFlux.build(builder, StableRNGs.StableRNG(1), 5, 3)
    @test chain.layers[1].weight != chain2.layers[1].weight

    chain3 = MLJFlux.build(builder, rng, 5, 3)
    @test chain.layers[1].weight != chain3.layers[1].weight

    conv_builder = MLJFlux.@builder begin
        front = Flux.Chain(Flux.Conv((3, 3), n_channels => 16), Flux.flatten)
        d = Flux.outputsize(front, (n_in..., n_channels, 1)) |> first
        Flux.Chain(front, Flux.Dense(d, n_out));
    end

    chain4 = MLJFlux.build(conv_builder, nothing, (5, 5), 3, 2)
    ps4 = Flux.params(chain4)
    @test size.(ps4) == [(3, 3, 2, 16), (16,), (3, 144), (3,)]
end


text:
using WorldVocoder
using Test

@testset "WorldVocoder.jl" begin
    fs = 44100
    x = [sinpi(2t * 200) for t in 0:1/fs:5]
    f0, timestamps = dio(x, fs)
    refined_f0 = stonemask(x, fs, timestamps, f0)
    spectrogram = cheaptrick(x, fs, timestamps, refined_f0)
    fft_size = (size(spectrogram, 1) - 1) * 2
    aperiodicity = d4c(x, fs, timestamps, refined_f0, fft_size)
    y = synthesis(refined_f0, spectrogram, aperiodicity, fs)

    f0, sp, ap = world(x, fs)
end

@testset "Harvest" begin
    fs = 44100
    x = [sinpi(2t * 200) for t in 0:1/fs:5]

    f0, temporal_positions = harvest(x, fs)
end


text:
export Invite

"""
An invite to a [`Guild`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/invite#invite-object).
"""
struct Invite
    code::String
    guild::Optional{Guild}
    channel::DiscordChannel
    approximate_presence_cound::Optional{Int}
    approximate_member_count::Optional{Int}
end
@boilerplate Invite :constructors :docs :lower :merge :mock


text:
"""
The `HeatmapCanvas` is also Unicode-based.
It has a half the resolution of the `BlockCanvas`.
This canvas effectively turns every character
into two pixels (top and bottom).
"""
struct HeatmapCanvas <: LookupCanvas
    grid::Array{UInt8,2}
    colors::Array{ColorType,2}
    blend::Bool
    visible::Bool
    pixel_width::Int
    pixel_height::Int
    origin_x::Float64
    origin_y::Float64
    width::Float64
    height::Float64
    xscale::Union{Symbol,Function}
    yscale::Union{Symbol,Function}
end

const HALF_BLOCK = '▄'
@inline x_pixel_per_char(::Type{HeatmapCanvas}) = 1
@inline y_pixel_per_char(::Type{HeatmapCanvas}) = 2

@inline lookup_encode(c::HeatmapCanvas) = [0 0; 1 1]
@inline lookup_decode(c::HeatmapCanvas) = [HALF_BLOCK; HALF_BLOCK]

@inline nrows(c::HeatmapCanvas) = div(size(grid(c), 2) + 1, 2)

HeatmapCanvas(args...; kwargs...) = CreateLookupCanvas(
    HeatmapCanvas,
    args...;
    min_char_width = 1,
    min_char_height = 1,
    kwargs...,
)

_toCrayon(c) = c === nothing ? 0 : (c isa Unsigned ? Int(c) : c)

function printrow(io::IO, c::HeatmapCanvas, row::Int)
    0 < row <= nrows(c) || throw(ArgumentError("Argument row out of bounds: $row"))
    y = 2row

    # extend the plot upwards by half a row
    isodd(size(grid(c), 2)) && (y -= 1)

    iscolor = get(io, :color, false)
    for x in 1:ncols(c)
        if iscolor
            fgcol = _toCrayon(c.colors[x, y])
            if y > 1
                bgcol = _toCrayon(c.colors[x, y - 1])
                print(io, Crayon(foreground = fgcol, background = bgcol), HALF_BLOCK)
                # for odd numbers of rows, only print the foreground for the top row
            else
                print(io, Crayon(foreground = fgcol), HALF_BLOCK)
            end
        else
            print(io, HALF_BLOCK)
        end
    end
    iscolor && print(io, Crayon(reset = true))

    nothing
end


text:
__precompile__()

module GenomicFeatures

export
    Strand,
    STRAND_NA,
    STRAND_POS,
    STRAND_NEG,
    STRAND_BOTH,

    Interval,
    seqname,
    leftposition,
    rightposition,
    strand,
    metadata,
    isoverlapping,

    IntervalCollection,
    eachoverlap,
    coverage,

    isfilled,
    hasseqname,
    hasleftposition,
    hasrightposition

import BioGenerics: BioGenerics, seqname, leftposition, rightposition, isoverlapping, isfilled, hasseqname, hasleftposition, hasrightposition, metadata
import DataStructures
import IntervalTrees
import Base.@propagate_inbounds

include("strand.jl")
include("interval.jl")
include("intervalcollection.jl")
include("queue.jl")
include("overlap.jl")
include("coverage.jl")

end # module


text:
# psuedo code for running steps on multiple channels in a single thread in fault tolerant way

function stepsloop!(channels::Vector{Channel}, steps::Vector{AbstractStep})
	while true
		try
			dosteps!(channels, steps)
		catch
			if isa(ex, InterruptException) # break on Ctrl-C for clean exit
				println("stepsloop ended by Ctrl-C")
				break
			end
		end
		if all([isdone(c,steps) for c in channels])
			println("stepsloop ended because all channels isdone=true")
			break
		end
	end
end

function dosteps!(channels::Vector{Channel})
	for c in channels
		try
			dosteps!(c)
		catch
			if isa(ex, InterruptException) # rethrow Ctrl-C
				rethrow(ex)
			end
		end
	end
end

function dosteps!(c::Channel, steps::Vector{AbstractStep})
	haskey(c, :exception_info)  && return # stop doing steps on this channel afte error
	workdone_cumulative = c[:workdone_cumulative]
	stepelapsed_cumulative = c[:stepelapsed_cumulative]
	for (i,s) in enumerate(steps)
		try
			tstart = time()
			wu = workunits(dostep!(c,s))
			stepelapsed = time()-tstart
			workdone_cumulative[j] += wu
			stepelapsed_cumulative[j] += stepelapsed
		catch ex
			if isa(ex, InterruptException) # rethrow Ctrl-C
				rethrow(ex)
			else # store a backtrace of the exception
				ex_str = sprint(io->showerror(io, ex))) do io
					showerror(io, ex)
					print(io,"\n")
					Base.show_backtrace(io, catch_backtrace())
					end
				c[:exception_info] = ex_str
				return # stop doing steps on this channel afte error
			end
		end
	end
end
dosteps!(c::Channel) = dosteps(c, c[:steps])

#### TODO ####
 # switch channel to Dict{AbstractString, Any}
 # put steps, workdone, stepelased, exception_info into channel


text:
"""

```
icell3D
```

Returns the i,j,k coordinates of the cell from the coordinates of an atom and box properties

"""
function icell3D(x::AbstractVector, box::Box)
  i = trunc(Int,(x[1]+box.sides[1]/2)/box.l[1])+1
  j = trunc(Int,(x[2]+box.sides[2]/2)/box.l[2])+1
  k = trunc(Int,(x[3]+box.sides[3]/2)/box.l[3])+1
  return i, j, k
end

#
# given the 1D representation of the index of a 3D cell, returns the 3 indexes of the
# coordinates in 3 dimensions
#

# if the number of cells in each dimension is different

function icell3D(n::AbstractVector{Int}, i1D)
  jk = i1D%(n[2]*n[3])
  if jk == 0
    jk = n[2]*n[3]
  end
  j, k = icell2D(n[2:3],jk)
  i = round(Int,(i1D-jk)/(n[2]*n[3]))+1
  return i, j, k
end

# If the number of cells in each dimension is the same

function icell3D(n::Int, i1D::Int)
  n2 = n^2
  jk = i1D%n2
  if jk == 0
    jk = n2
  end
  j, k = icell2D(n,jk)
  i = round(Int,(i1D-jk)/(n2))+1
  return i, j, k
end

# Given the 1D represetation of the index of a 2D cell, returns the 2
# indexes of the 2D representation

# If the two dimensions are different

function icell2D(n::AbstractVector{Int}, i1D::Int)
  j = i1D%n[2]
  if j == 0
    j = n[2]
  end
  i = round(Int,(i1D - j)/n[2])+1
  return i, j
end

# If the two dimensions are the same

function icell2D(n::Int, i1D::Int)
  j = i1D%n
  if j == 0
    j = n
  end
  i = round(Int,(i1D - j)/n)+1
  return i, j
end


text:
export mk_tree

struct Node
    id::Integer
    parent::Union{Node,Nothing}
    sorg::Integer
    tail::Integer
    cplen::Int
    pflen::Int
    extype::Char
    excep::Union{Any,Nothing}
    children::Array{Node,1}
    at
end

function mk_key(node, cplen, extype, excep)
    key = (node, cplen, extype, excep)
    return key
end

"""
    lcp(text,i1,i2)

    Compute the length of the common prefix of text[i1:end] and text[i2:end]

"""
function lcp(text, i1, i2)

    l = 0

    @views for (a, b) in zip(text[i1:end], text[i2:end])

        if a == b
            l += 1
        else
            break
        end
    end

    return l
end

"""
    prefixparts(node)

    return list of prefixparts and list of positions
"""
function prefixparts(node)

    pfx = []

    cur = node

    while !isnothing(cur.parent)
        push!(pfx, text[cur.parent.tail:cur.parent.tail*cur.cplen])
        cur = cur.parent
    end

    return reverse(pfx)
end

"""
    mk_tree(text)

    compute a tree with unique prefixes in text
"""
function mk_tree(text)

    function addNode(newnode, key)

        push!(nodes, newnode)

        if !isnothing(parent)
            push!(newnode.parent.children, newnode)
        end

        if !isnothing(key)
            exceptions[key] = newnode
        end

    end

    function get_next(node, cplen, extype, excep)
        key = mk_key(node, cplen, extype, excep)
        get(exceptions, key, nothing)
    end

    """
        add(s,sref)

        add the string starting at position s to the tree. ref is e.g. id van string
    """
    function add(s, sref)

        sorg = s

        next = root  # start with root
        pflen = 0  # length of prefix

        while !isnothing(next)
            cur = next

            if s == cur.tail
                push!(cur.at, sref)
                return
            end

            cplen = lcp(text, cur.tail, s)

            if cplen == 0

                excep = text[s]
                extype = '!'

                next = get_next(cur, cplen, extype, excep)
                if !isnothing(next)
                    continue
                else
                    key = mk_key(cur, cplen, extype, excep)

                    addNode(mk_node(parent = cur, tail = s, sorg = sorg, pflen = (cur.pflen + cplen), cplen = cplen, extype = '!', excep = excep,
                            ref = sref), key)
                    return
                end
            end

            if cplen == (lengt - cur.tail) && (lengt - s) > cplen && cur != root

                s_before = cur.tail

                newnode = mk_node(parent = cur, tail = cur.tail + cplen, sorg = cur.sorg, pflen = (cur.pflen + lengt - s_before),
                    cplen = lengt - s_before, extype = '$', excep = cur.excep)
                newnode.at = cur.at
                cur.tail = s
                cur.ref = sref
                cur.at = [sref]
                cur.sorg = sorg

                extype = '$'
                key = make_key(cur, lengt - s_before, extype, None)
                addNode(newnode, key)
                return
            end

            if (lengt - s) >= cplen
                s += cplen
                excep = text[s]
                extype = '!'
            else
                extype = '$'
                excep = nothing
                s = length(text) + 1
            end

            next = get_next(cur, cplen, extype, excep)

            if !isnothing(next)
                continue
            else
                key = mk_key(cur, cplen, extype, excep)
                addNode(
                    mk_node(parent = cur, tail = s, sorg = sorg, pflen = cur.pflen + cplen, cplen = cplen, extype = extype, excep = excep,
                        ref = sref), key)
                return
            end
        end

    end

    nodes = []

    function mk_node(; parent, tail, sorg = nothing, cplen = -1, pflen = -1, extype, excep = nothing, ref = nothing)

        id = 1 + length(nodes)
        node = Node(id, parent, sorg, tail, cplen, pflen, extype, excep, [], isnothing(ref) ? [] : [ref])
        return node
    end

    lengt = length(text)
    root = mk_node(parent = nothing, tail = length(text) + 1, pflen = 0, cplen = 0, sorg = 1, extype = '$')

    nodes = [root]

    exceptions = Dict()

    for (i, c) in enumerate(text)
        add(i, i)
    end

    return (nodes = nodes, root = root)
end


text:
function generaSombraTeor_v3(ps_volteor, matConexionVertices, vecVertices, ps_publico, ps_calles)

    V_ = ps_volteor.Vertices[1]
    numVerticesTotales, numDim = size(V_);

    conjuntoAlturas = unique(V_[:,3])
    numAlturas = length(conjuntoAlturas)

    numVerticesBase = Int(round(numVerticesTotales / numAlturas));
    ps_SombraVolTeor_p = PolyShape([],0)
    ps_SombraVolTeor_o = PolyShape([],0)
    ps_SombraVolTeor_s = PolyShape([],0)

    for k = 1:numAlturas-1
        for j = 1:length(vecVertices[k])
            if j == length(vecVertices[k])
                jl0 = Int.(vecVertices[k][j])
                jl1 = Int.(vecVertices[k][1])
            else
                jl0 = Int.(vecVertices[k][j])
                jl1 = Int.(vecVertices[k][j + 1])
            end
            x1 = V_[jl0,1]; y1 = V_[jl0,2]; z1 = V_[jl0,3];
            x2 = V_[jl1,1]; y2 = V_[jl1,2]; z2 = V_[jl1,3];
            ju0 = Int.(matConexionVertices[matConexionVertices[:,1] .== jl0,2])[1]
            ju1 = Int.(matConexionVertices[matConexionVertices[:,1] .== jl1,2])[1]

            x3 = V_[ju1,1]; y3 = V_[ju1,2]; z3 = V_[ju1,3];
            x4 = V_[ju0,1]; y4 = V_[ju0,2]; z4 = V_[ju0,3];

            verts_p = [[x1 - z1/0.49  y1];
                       [x2 - z2/0.49  y2];
                       [x3 - z3/0.49  y3];
                       [x4 - z4/0.49  y4]]

            verts_o = [[x1 + z1/0.49  y1];
                       [x2 + z2/0.49  y2];
                       [x3 + z3/0.49  y3];
                       [x4 + z4/0.49  y4]]

            verts_s = [[x1  y1 - z1/1.54];
                       [x2  y2 - z2/1.54];
                       [x3  y3 - z3/1.54];
                       [x4  y4 - z4/1.54]]

            if ps_SombraVolTeor_p.NumRegions == 0
                ps_SombraVolTeor_p.Vertices = push!(ps_SombraVolTeor_p.Vertices, verts_p)
                ps_SombraVolTeor_p.NumRegions = 1
                ps_SombraVolTeor_o.Vertices = push!(ps_SombraVolTeor_o.Vertices, verts_o)
                ps_SombraVolTeor_o.NumRegions = 1
                ps_SombraVolTeor_s.Vertices = push!(ps_SombraVolTeor_s.Vertices, verts_s)
                ps_SombraVolTeor_s.NumRegions = 1
            else
                ps_SombraVolTeor_p = polyShape.polyUnion_v2(ps_SombraVolTeor_p, PolyShape([verts_p],1))
                ps_SombraVolTeor_o = polyShape.polyUnion_v2(ps_SombraVolTeor_o, PolyShape([verts_o],1))
                ps_SombraVolTeor_s = polyShape.polyUnion_v2(ps_SombraVolTeor_s, PolyShape([verts_s],1))
            end
        end

    end

    p_p = polyShape.polyDifference_v3(ps_SombraVolTeor_p, ps_publico)
    if length(p_p.Vertices) > 0
        ps_SombraVolTeor_p = PolyShape(p_p.Vertices, length(p_p.Vertices))
    else
        ps_SombraVolTeor_p = PolyShape([], 0)
    end
    p_o = polyShape.polyDifference_v3(ps_SombraVolTeor_o, ps_publico)
    if length(p_o.Vertices) > 0
        ps_SombraVolTeor_o = PolyShape(p_o.Vertices, length(p_o.Vertices))
    else
        ps_SombraVolTeor_o = PolyShape([], 0)
    end
    p_s = polyShape.polyDifference_v3(ps_SombraVolTeor_s, ps_publico)
    if length(p_s.Vertices) > 0
        ps_SombraVolTeor_s = PolyShape(p_s.Vertices, length(p_s.Vertices))
    else
        ps_SombraVolTeor_s = PolyShape([], 0)
    end

    ps_SombraVolTeor_p = polyShape.polyDifference_v3(ps_SombraVolTeor_p, ps_calles)
    ps_SombraVolTeor_o = polyShape.polyDifference_v3(ps_SombraVolTeor_o, ps_calles)
    ps_SombraVolTeor_s = polyShape.polyDifference_v3(ps_SombraVolTeor_s, ps_calles)

    return ps_SombraVolTeor_p, ps_SombraVolTeor_o, ps_SombraVolTeor_s
end


text:
function appearsmorethanonce(var::SymPy.Sym,system::Array{SymPy.Sym})
    appeared = false
    for expression in system
        if var in SymPy.free_symbols(expression)
            if appeared
                return true
            else
                appeared = true
            end
        end
    end
    return false
end

function findstaticexpressions(expressions_sym::Array{SymPy.Sym},leads_sym::Array{SymPy.Sym},contemps_sym::Array{SymPy.Sym},states_sym::Array{SymPy.Sym},var_isstate::BitArray,eq_isexo::BitArray,var_isexo::BitArray)
    numvars = length(contemps_sym)
    var_isstatic = copy(.!var_isstate)
    expression_isstatic = trues(length(expressions_sym))
    for (i,expression) in enumerate(expressions_sym)
        for (j,lead) in enumerate(leads_sym)
            if lead in SymPy.free_symbols(expression)
                var_isstatic[j] = false
                expression_isstatic[i] = false
            end
        end
        # Only check this if haven't already flipped it
        if expression_isstatic[i]
            for state in states_sym
                if state in SymPy.free_symbols(expression)
                    expression_isstatic[i] = false
                end
            end
        end
    end

    # Finally, if a static variable appears only _once_ in the whole system, the equation it appears in has to be a static
    tocheck = BitArray(numvars)
    fill!(tocheck,false)
    for (i,var) in enumerate(contemps_sym)
        if var_isstatic[i] && .!appearsmorethanonce(var,expressions_sym)
            tocheck[i] = true
        end
    end
    for var in contemps_sym[tocheck]
        for (i,expression) in enumerate(expressions_sym)
            if var in SymPy.free_symbols(expression)
                expression_isstatic[i] = true
            end
        end
    end

    # Now negate out any exos, since they'll be solved separately
    expression_isstatic = expression_isstatic .& .!eq_isexo
    var_isstatic = var_isstatic .& .!var_isexo

    return expression_isstatic,var_isstatic
end

function expressioncontains(expression::SymPy.Sym,symbols::Array{SymPy.Sym},shortcircuit::Bool)
    # This is a short-circuit version of expression contains that is equivalent but quicker than any(expressioncontains(expression,symbols))
    for symbol in symbols
        if symbol in SymPy.free_symbols(expression)
            return true
        end
    end
    return false
end

function expressioncontains(expression::SymPy.Sym,symbols::Array{SymPy.Sym})
    count = 0
    isin = BitArray(length(symbols))
    fill!(isin,false)
    for (i,symbol) in enumerate(symbols)
        if symbol in SymPy.free_symbols(expression)
            isin[i] = true
        end
    end
    return isin
end

function checkifcontainsforward(expressions::Array{SymPy.Sym},leads::Array{SymPy.Sym})
    contains = BitArray(length(expressions))
    fill!(contains,false)
    for (i,expression) in enumerate(expressions)
        contains[i] = expressioncontains(expression,leads,true)
    end
    return contains
end

function findexoequations(expressions::Array{SymPy.Sym},leads::Array{SymPy.Sym},contemps::Array{SymPy.Sym})
    contains_forward = checkifcontainsforward(expressions,leads)
    restrictedindices = collect(1:length(expressions))
    restrictedsystem = expressions[.!contains_forward]
    restrictedindices = restrictedindices[.!contains_forward]
    contempindices = collect(1:length(contemps))

    endocontemps = copy(contemps)
    exocontemps = Array{Int}(0)
    exoequations = Array{Int}(0)

    somethingchanged = true
    while somethingchanged
        somethingchanged = false
        keepexpression = trues(length(restrictedsystem))
        for (i,(expression,expression_index)) in enumerate(zip(restrictedsystem,restrictedindices))
            indices = expressioncontains(expression,endocontemps)
            if countnz(indices) == 1
                # Now figure out which one
                newexo = contempindices[indices]
                # Remove it from the list of endos
                endocontemps = endocontemps[.!indices]
                contempindices = contempindices[.!indices]
                push!(exocontemps,newexo[1])
                push!(exoequations,expression_index)
                keepexpression[i] = false
                somethingchanged = true
            end
        end
        # Now replace restricted system with just the equation we want to keep
        restrictedsystem = restrictedsystem[keepexpression]
        restrictedindices = restrictedindices[keepexpression]
    end

    # Convert to boolean arrays
    expression_isexo = BitArray(length(expressions))
    fill!(expression_isexo,false)
    expression_isexo[exoequations] = true
    var_isexo = BitArray(length(contemps))
    fill!(var_isexo,false)
    var_isexo[exocontemps] = true
    if countnz(expression_isexo) != countnz(var_isexo)
        error("Something has done wrong. Found $(countnz(expression_isexo)) exogenous equations, but $(countnz(var_isexo)) exogenous vars.")
    end

    return expression_isexo,var_isexo
end


text:
"""
    clusters(support::PP, cluster::Function)

generate clusters with centers from support
"""
clusters(support, cluster) = reduce(hcat, [cluster() .+ c for c in eachcol(support)])

function surface_points(n, support, metric)
    all_points = rand(size(support, 1), 1000 * n)

    a = [metric(x, support) for x in eachcol(all_points)]

    ix = sortperm(a)[1:n]
    all_points[:, ix]
end

voronoi(n, support) = surface_points(n, support, voronoi_metric)

circles(n, support, r) = surface_points(n, support, (x, s) -> circle_metric(x, s, r))


text:
parse_date(::Missing) = missing
parse_time(::Missing) = missing
parse_datetime(::Missing) = missing
parse_int(::Missing) = missing
parse_date(d) = Date(d)
parse_time(d) = Time(d)
parse_datetime(d) = DateTime(d)
parse_int(x) = Int(x)

check_if_date(d::Missing) = true
check_if_time(d::Missing) = true
check_if_datetime(d::Missing) = true
function check_if_date(d::AbstractString)
    x = match(r"\d{4}-\d{1,2}-\d{1,2}", d)
    x === nothing && return false
    x.match == d && tryparse(Date, d) !== nothing
end
function check_if_time(d::AbstractString)
    x = match(r"\d{2}:\d{2}:\d{2}\.\d{1,4}", d)
    x === nothing && return false
    x.match == d && tryparse(Time, d) !== nothing
end
function check_if_datetime(d::AbstractString)
    x = match(r"\d{4}-\d{1,2}-\d{1,2}T\d{2}:\d{2}:\d{2}\.\d{1,4}", d)
    x === nothing && return false
    x.match == d && tryparse(DateTime, d) !== nothing
end
function check_if_date(d::AbstractVector)
    for x in d
        check_if_date(x) && continue
        return false
    end
    return true
end
function check_if_time(d::AbstractVector)
    for x in d
        check_if_time(x) && continue
        return false
    end
    return true
end
function check_if_datetime(d::AbstractVector)
    for x in d
        check_if_date(x) && continue
        return false
    end
    return true
end

function integer_or_missing(x::AbstractVector)
    for i in x
        ismissing(i) && continue
        isinteger(i) && continue
        return false
    end
    return true
end

"""
`modify_col!` tries to identify the real type of a column, especially for strings
that are actually dates or floats that are actually integers. Almost all data downloaded
from WRDS correctly specifies dates, but all numbers are stored as float8, even items like
`year` which should be an integer. This uses multiple dispatch to check if all elements
in a given column are compatible with changing type and then changes the type of the
column. Note that for strings this function does not by default try to convert integer
like strings to integer (such as GVKey), it only converts strings that look like a date,
datetime, or time.
"""
function modify_col!(df::AbstractDataFrame, col::String, ::Type{<:AbstractString})
    if check_if_date(df[:, col])
        df[!, col] = parse_date.(df[:, col])
    elseif check_if_time(df[:, col])
        df[!, col] = parse_time.(df[:, col])
    elseif check_if_datetime(df[:, col])
        df[!, col] = parse_datetime.(df[:, col])
    end
end

function modify_col!(df::AbstractDataFrame, col::String, ::Type{<:Real})
    if integer_or_missing(df[:, col])
        df[!, col] = parse_int.(df[:, col])
    end
end

function modify_col!(df::AbstractDataFrame, col::String, ::Type{Union{Missing, A}}) where {A <: Any}
    all_missing = all(ismissing.(df[:, col]))
    if !all_missing # if all items are missing, do not modify the column since it is uncertain what
        # it should be modified to
        modify_col!(df, col, A)
        if !any(ismissing.(df[:, col]))
            disallowmissing!(df, col)
        end
    end
end

function modify_col!(df::AbstractDataFrame, col::String, ::Type{<:Any})
end

function run_sql_query(
    conn::LibPQ.Connection,
    q::AbstractString
)
    df = LibPQ.execute(conn, q) |> DataFrame
    for col in names(df)
        modify_col!(df, col, eltype(df[:, col]))
    end
    df
end

function run_sql_query(
    conn::DBInterface.Connection,
    q::AbstractString
)
    df = DBInterface.execute(conn, q) |> DataFrame
    for col in names(df)
        modify_col!(df, col, eltype(df[:, col]))
    end
    df
end

struct Conditions
    fun::Function
    l::Union{Symbol, String}
    r::Union{Symbol, String}
end

function Conditions(
    l::Union{Symbol, String},
    fun::Function,
    r::Union{Symbol, String}
)
    Conditions(
        fun,
        l,
        r
    )
end

function Conditions(
    x::Tuple{Function, Symbol, Symbol}
)
    Conditions(
        x[1],
        x[2],
        x[3]
    )
end

function range_join(
    df1::DataFrame,
    df2::DataFrame,
    on,
    conditions::Array{Tuple{Function, Symbol, Symbol}};
    minimize=nothing,
    join_conditions::Union{Array{Symbol}, Symbol}=:and,
    validate::Tuple{Bool, Bool}=(false, false),
    joinfun::Function=leftjoin
)
    jointype = if joinfun == leftjoin
        :left
    elseif joinfun == rightjoin
        :right
    elseif joinfun == innerjoin
        :inner
    else
        :outer
    end

    range_join(
        df1,
        df2,
        on,
        Conditions.(conditions);
        minimize,
        join_conditions,
        validate=(validate[2], validate[1]),
        jointype
    )
end

function change_function(fun)
    if fun == <=
        return >=
    elseif fun == >=
        return <=
    elseif fun == >
        return <
    elseif fun == <
        return >
    else
        return !fun
    end
end

function filter_data(
    df_row::DataFrameRow,
    df_partial::AbstractDataFrame,
    conditions::Array{Conditions};
    join_conditions::Union{Array{Symbol}, Symbol}=:and,
)

    fil = ones(Bool, nrow(df_partial))

    for (j, condition) in enumerate(conditions)
        temp_join = if typeof(join_conditions) <: Symbol
            join_conditions
        elseif j == 1 # if it is the first time through, all the current values are
            :and
        else
            join_conditions[j-1]
        end
        if temp_join == :and
            fil = fil .& broadcast(
                condition.fun,
                df_row[condition.l],
                df_partial[:, condition.r]
            )
        else
            fil = fil .| broadcast(
                condition.fun,
                df_row[condition.l],
                df_partial[:, condition.r]
            )
        end
    end
    return fil
end

"""
I need the keymap returned instead of the full
SubDataFrame
"""
function special_get(gdf, key)
    if haskey(gdf.keymap, key)
        x = gdf.keymap[key]
        return (gdf.starts[x], gdf.ends[x])
    else
        return (0, 0)
    end
end

r_col_names(x::Conditions) = x.r
l_col_names(x::Conditions) = x.l

function validate_error(df)
    if nrow(df) > 0
        if nrow(df) == 1
            error_message = "df1 contains 1 duplicate key: " *
                            "$(NamedTuple(df[1, :])). "
        elseif nrow(df) == 2
            error_message = "df1 contains 2 duplicate keys: " *
                            "$(NamedTuple(df[1, :])) and " *
                            "$(NamedTuple(df[2, :])). "
        else
            error_message = "df1 contains $(nrow(df)) duplicate keys: " *
                            "$(NamedTuple(df[1, :])), ..., " *
                            "$(NamedTuple(df[end, :])). "
        end
        error(error_message)
    end
end

"""

    function range_join(
        df1::DataFrame,
        df2::DataFrame,
        on,
        conditions::Array{Conditions};
        minimize=nothing,
        join_conditions::Union{Array{Symbol}, Symbol}=:and,
        validate::Tuple{Bool, Bool}=(false, false),
        jointype::Symbol=:inner
    )
Joins the dataframes based on a series of conditions, designed
to work with ranges

### Arguments
- `df1::DataFrame`: left DataFrame
- `df2::DataFrame`: right DataFrame
- `on`: either array of column names or matched pairs
- `conditions::Array{Conditions}`: array of `Conditions`, which specifies the function (<=, >, etc.), left and right columns
- `jointype::Symbol=:inner`: type of join, options are :inner, :outer, :left, and :right
- `minimize`: either `nothing` or an array of column names or matched pairs, minimization will take place in order
- `join_conditions::Union{Array{Symbol}, Symbol}`: defaults to `:and`, otherwise an array of symbols that is 1 less than the length of conditions that the joins will happen in (:or or :and)
- `validate::Tuple{Bool, Bool}`: Whether to validate, this works differently than the equivalent in DataFrames joins,
  here, validate insures that a single row from the dataframe is not duplicated. So validate=(true, false) means that
  there are no duplicated rows from the left dataframe in the result.

### Example

```
df1 = DataFrame(
    firm=1:10,
    date=Date.(2013, 1:10, 1:10)
)

df2 = df1[:, :]
df2[!, :date_low] = df2.date .- Day(2)
df2[!, :date_high] = df2.date .+ Day(2)
select!(df2, Not(:date))

range_join(
    df1,
    df2,
    [:firm],
    [
        Conditions(<, :date, :date_high),
        Conditions(>, :date, :date_low)
    ],
    join_conditions=[:and]
)
```

"""
function range_join(
    df1::DataFrame,
    df2::DataFrame,
    on,
    conditions::Array{Conditions};
    minimize=nothing,
    join_conditions::Union{Array{Symbol}, Symbol}=:and,
    validate::Tuple{Bool, Bool}=(false, false),
    jointype::Symbol=:inner,
)

    if nrow(df1) > nrow(df2) && minimize === nothing
        # the fewer main loops this goes through
        # the faster it is overall, if they are about equal
        # this likely slows things down, but it is
        # significantly faster if it is flipped for large sets
        # I added minimize === nothing since, under the current setup
        # the minimization is done as minimize assumes a groupby on
        # the left, so it does not work properly if this goes through
        # one thought on how to fix that is to create a minimize object
        # that would allow for setting the condition under which
        # minimizaiton makes sense
        new_cond = Conditions[]
        for con in conditions
            push!(new_cond, Conditions(change_function(con.fun), con.r, con.l))
        end
        on1, on2 = parse_ons(on)
        new_join = if jointype == :left
            new_join = :right
        elseif jointype == :right
            new_join = :left
        else
            new_join = jointype
        end
        return range_join(
            df2,
            df1,
            on2 .=> on1,
            new_cond;
            join_conditions,
            validate=(validate[2], validate[1]),
            jointype=new_join,
        )
    end

    on1, on2 = WRDSMerger.parse_ons(on)

    res_left = Int[]
    res_right = Int[]

    sizehint!(res_left, max(nrow(df1), nrow(df2)))
    sizehint!(res_right, max(nrow(df1), nrow(df2)))

    gdf = groupby(df2, on2)

    idx = gdf.idx

    # under the current setup, I do not think threading will work here
    for (i, key) in enumerate(Tuple.(copy.(eachrow(df1[:, on1]))))

        s, e = special_get(gdf, key)

        e == 0 && continue

        fil = if join_conditions == :and
            broadcast(
                &,
                [
                    broadcast(
                        condition.fun,
                        df1[i, condition.l],
                        df2[idx[s:e], condition.r]
                    )
                    for condition in conditions
                ]...
            )
        elseif join_conditions == :or
            broadcast(
                |,
                [
                    broadcast(
                        condition.fun,
                        df1[i, condition.l],
                        df2[idx[s:e], condition.r]
                    )
                    for condition in conditions
                ]...
            )
        else
            WRDSMerger.filter_data(
                df1[i, :],
                df2[idx[s:e], :],
                conditions;
                join_conditions
            )
        end

        sum(fil) == 0 && continue
        cur = length(res_left)
        to_grow = sum(fil) + cur
        resize!(res_left, to_grow)
        resize!(res_right, to_grow)
        res_left[cur+1:end] .= i
        @inbounds res_right[cur+1:end] = idx[s:e][fil]

    end
    if minimize !== nothing
        min1, min2 = parse_ons(minimize)
        df_temp = hcat(
            df1[res_left, string.(min1) |> unique],
            df2[res_right, string.(min2) |> unique]
        )
        df_temp[!, :_idx_left] = res_left
        df_temp[!, :_idx_right] = res_right
        for (l, r) in zip(min1, min2)
            #group_col = x.group_col_left ? :_idx_left : :_idx_right
            group_col = :_idx_left
            df_temp = subset(
                groupby(df_temp, group_col),
                [l, r] =>
                (l, r) -> abs.(l .- r) .== minimum(abs.(l .- r))
            )
        end
        res_left = df_temp[:, :_idx_left]
        res_right = df_temp[:, :_idx_right]
    end
    df = if jointype == :right
        hcat(df2[res_right, :], select(df1[res_left, :], Not(on1)))
    else
        hcat(df1[res_left, :], select(df2[res_right, :], Not(on2)))
    end

    if any(validate)
        df[!, :_idx_left] = res_left
        df[!, :_idx_right] = res_right
        if validate[1]
            cols = string.(vcat(on1, l_col_names.(conditions))) |> unique
            temp = df[nonunique(df[:, [:_idx_left]]), cols]
            validate_error(temp)
        end

        if validate[2]
            cols = string.(vcat(on1, r_col_names.(conditions))) |> unique
            temp = df[nonunique(df[:, [:_idx_right]]), cols]
            validate_error(temp)
        end
        select!(df, Not([:_idx_left, :_idx_right]))
    end

    if jointype == :left || jointype == :outer
        idx_l_add = DataFrames.find_missing_idxs(res_left, nrow(df1))
        temp = df1[idx_l_add, :]
        insertcols!(temp, [col => missing for col in names(df) if col ∉ names(temp)]...)
        df = vcat(df, temp)
    end
    if jointype == :right || jointype == :outer
        idx_r_add = DataFrames.find_missing_idxs(res_right, nrow(df2))
        temp = df2[idx_r_add, :]
        insertcols!(temp, [col => missing for col in names(df) if col ∉ names(temp)]...)
        df = vcat(df, temp)
    end

    return df
end

function parse_ons(on)
    on1 = String[]
    on2 = String[]
    for x in on
        if typeof(x) <: Pair
            push!(on1, String(x[1]))
            push!(on2, String(x[2]))
        else
            push!(on1, String(x))
            push!(on2, String(x))
        end
    end
    return on1, on2
end

# I wrote this macro a while ago, do not know if it currently works
# function parse_expression(
#     expression::Expr
# )
#     out = Expr[]
#     if expression.head == :call || expression.head == :comparison
#         push!(out, expression)
#         return out
#     end
#     for a in expression.args
#         if a.head == :&& || a.head == :||
#             out = vcat(out, parse_expression(a))
#         else a.head == :call
#             push!(out, a)
#         end
#     end
#     return out
# end

# function return_function(
#     val::Symbol
# )
#     if val == :<
#         return <
#     elseif val == :>
#         return >
#     elseif val == :<=
#         return <=
#     elseif val == :>=
#         return >=
#     else
#         error("Function Symbol must be a comparison")
#     end
# end

# function reverse_return_function(
#     val::Symbol
# )
#     if val == :<
#         return return_function(:>)
#     elseif val == :>
#         return return_function(:<)
#     elseif val == :<=
#         return return_function(:>=)
#     elseif val == :>=
#         return return_function(:<=)
#     else
#         error("Function Symbol must be a comparison")
#     end
# end

# function push_condition!(
#     conditions::Array{Tuple{Function, Symbol, Symbol}},
#     f::Symbol,
#     first::Expr,
#     second::Expr
# )
#     if first.args[1] == :left && second.args[1] == :right
#         push!(
#             conditions,
#             (
#                 return_function(f),
#                 eval(first.args[2]),
#                 eval(second.args[2])
#             )
#         )
#     elseif first.args[1] == :right && second.args[1] == :left
#         push!(
#             conditions,
#             (
#                 reverse_return_function(f),
#                 eval(second.args[2]),
#                 eval(first.args[2])
#             )
#         )
#     else
#         error("Comparison must have right and left as labels")
#     end
# end

# function expressions_to_conditions(
#     expressions::Array{Expr}
# )
#     out = Tuple{Function, Symbol, Symbol}[]
#     for x in expressions
#         if x.head == :call
#             push_condition!(
#                 out,
#                 x.args[1],
#                 x.args[2],
#                 x.args[3]
#             )
#         elseif x.head == :comparison
#             for i in 1:2:length(x.args)-1
#                 push_condition!(
#                     out,
#                     x.args[i+1],
#                     x.args[i],
#                     x.args[i+2]
#                 )
#             end
#         end
#     end
#     return out
# end

# function parse_expr(fil)
#     fil = string(fil)
#     for (s, r) in [("||", ") .| ("), ("&&", ") .& ("), ("<", ".<"), (">", ".>"), (r"left\.([^\s]*)", s"df1[i, :\1]"), (r"right\.([^\s]*)", s"temp[:, :\1]")]
#         fil = replace(fil, s => r)
#     end
#     fil = "($fil)"
#     Meta.parse(fil)
# end

# function join_helper(
#     df1,
#     df2,
#     on,
#     conditions,
#     args...
# )
#     #new_conditions = conditions |> parse_expression |> expressions_to_conditions
#     quote
#         $range_join(
#             $df1,
#             $df2,
#             $on,
#             $conditions;
#             $(args...)
#         )
#     end
# end

# macro join(
#     df1,
#     df2,
#     on,
#     conditions,
#     args...
# )
#     local new_conditions = conditions |> parse_expression |> expressions_to_conditions
#     #local aakws = [esc(a) for a in args]
#     esc(join_helper(df1, df2, on, new_conditions, args...))
# end


text:
# f(x,y)=x-y

function Fxmy(x::Float64,y::Float64)
    return ((x-y)*ones(1,1),nothing)
end

function Fxmy(x::Float64,y)
    return (x*ones(size(y))-y,nothing)
end

function Fxmy(x,y::Float64)
    return (x-y*ones(size(y)),nothing)
end

function Fxmy(x,y)
    if size(x)==(1,1)
        return (x[1]*ones(size(y))-y,nothing)
    elseif size(y)==(1,1)
        return (x-y[1]*ones(size(x)),nothing)
    else
        return (x.-y,nothing)
    end
end

function Fxmy_inplace(value,auxvalue,x,y)
    if size(x)==(1,1)
        copy!(value,x[1]*ones(size(y))); axpy!(-1.0,y,value)
    elseif size(y)==(1,1)
        copy!(value,x); axpy!(-1.0,y[1]*ones(size(x)),value)
    else
        copy!(value,x); axpy!(-1.0,y,value)
    end
end

function Dxmy(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if derivativeIDX==1
        if size(x)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c,grad_n)
        end
    elseif derivativeIDX==2
        if size(y)==(1,1)
            axpy!(-1.0,[sum(grad_c)],grad_n)
        else
            axpy!(-1.0,grad_c,grad_n)
        end
    end
end

if PROC=="GPU"
function Fxmy(x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            tmp=CudaArray(Float64,size(y))
            gfill!(tmp,x)
            axpy!(-1.0,y,tmp)
        elseif size(y)==(1,1)
            tmp=CudaArray(Float64,size(x))
            gfill!(tmp,y); scale!(-1.0,tmp)
            axpy!(1.0,x,tmp)
        else
            tmp=CudaArray(zeros(size(x))); copy!(tmp,x); axpy!(1.0,y,tmp)
        end
        return (tmp,nothing)
    end

    function Fxmy_inplace(value::CudaArray,auxvalue,x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            gfill!(value,x); axpy!(-1.0,y,value)
        elseif size(y)==(1,1)
            gfill!(value,y); scale!(-1.0,value); axpy!(1.0,x,value)
        else
            copy!(value,x); axpy!(1.0,y,value)
        end
    end

    function Dxmy(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)

        if derivativeIDX==1
            if size(x)==(1,1)
                axpy!(1.0,sum(grad_c),grad_n)
            else
                axpy!(1.0,grad_c,grad_n)
            end
        elseif derivativeIDX==2
            if size(y)==(1,1)
                axpy!(-1.0,sum(grad_c),grad_n)
            else
                axpy!(-1.0,grad_c,grad_n)
            end
        end
    end

end

Derivative[Fxmy]=Dxmy
Inplace[Fxmy]=Fxmy_inplace

import Base.-
-(A::ADnode,B::ADnode)=ADnode(Fxmy,[A B])
-(A::Real,B::ADnode)=ADnode(Fxmy,[ADconst(A) B])
-(A::ADnode,B::Real)=ADnode(Fxmy,[A ADconst(B)])

export Fxmy, -

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
## Copyright (c) 2015, Júlio Hoffimann Mendes <juliohm@stanford.edu>
##
## Permission to use, copy, modify, and/or distribute this software for any
## purpose with or without fee is hereby granted, provided that the above
## copyright notice and this permission notice appear in all copies.
##
## THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
## WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
## MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
## ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
## WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
## ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
## OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

"""
    voxelreuse(training_image::AbstractArray,
               tplsizex::Integer, tplsizey::Integer, tplsizez::Integer;
               overlapx::Real=1/6, overlapy::Real=1/6, overlapz::Real=1/6,
               cut::Symbol=:boykov, simplex::Bool=false, nreal::Integer=10,
               threads::Integer=CPU_PHYSICAL_CORES, gpu::Bool=false,
			   soft::AbstractVector=[], hard::HardData=HardData(), tol::Real=.1)

Returns the mean voxel reuse in `[0,1]` and its standard deviation.

### Notes

The approximation gets better as `nreal` is made larger.
"""
function voxelreuse(training_image::AbstractArray,
                    tplsizex::Integer, tplsizey::Integer, tplsizez::Integer;
                    overlapx::Real=1/6, overlapy::Real=1/6, overlapz::Real=1/6,
                    cut::Symbol=:boykov, simplex::Bool=false, nreal::Integer=10,
                    threads::Integer=CPU_PHYSICAL_CORES, gpu::Bool=false,
					soft::AbstractVector=[], hard::HardData=HardData(), tol::Real=.1)

  # calculate the overlap from given percentage
  ovx = ceil(Int, overlapx * tplsizex)
  ovy = ceil(Int, overlapy * tplsizey)
  ovz = ceil(Int, overlapz * tplsizez)

  # elementary raster path
  ntilex = ovx > 1 ? 2 : 1
  ntiley = ovy > 1 ? 2 : 1
  ntilez = ovz > 1 ? 2 : 1

  # simulation grid dimensions
  gridsizex = ntilex * (tplsizex - ovx) + ovx
  gridsizey = ntiley * (tplsizey - ovy) + ovy
  gridsizez = ntilez * (tplsizez - ovz) + ovz

  _, _, voxs = iqsim(training_image,
                     tplsizex, tplsizey, tplsizez,
                     gridsizex, gridsizey, gridsizez,
                     overlapx=overlapx, overlapy=overlapy, overlapz=overlapz,
                     cut=cut, simplex=simplex, nreal=nreal,
                     threads=threads, gpu=gpu, debug=true,
					 soft=soft,hard=hard,tol=tol)

  μ = mean(voxs)
  σ = std(voxs, mean=μ)

  μ, σ
end


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.12.7/agnostic/x64/nccl_2.12.7-1+cuda10.2_x86_64.txz",
                      "922b875850e25d209e0c73da7f2f123e6c88a2667002af7212295797401da7ee")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.12.7/agnostic/ppc64le/nccl_2.12.7-1+cuda10.2_ppc64le.txz",
                      "870185ba54e6e1aee0909ea2f43cea6b643a8da15d033d62b5e0867c8ad35d8e")],
)

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using Documenter, ParticleDA

# Load the `Model` module from the tests to show the docstring of the parameters
test_dir = joinpath(dirname(pathof(ParticleDA)), "..", "test")
module_src = joinpath(test_dir, "model", "model.jl")
include(module_src)
using .Model

makedocs(
    modules = [ParticleDA, Model],
    sitename = "ParticleDA",
)

deploydocs(
    repo = "github.com/Team-RADDISH/ParticleDA.jl",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using Dates
tags = split(read(`git tag`,String))
version_base = "$(Dates.year(Dates.today())-2000).$(Dates.dayofyear(Dates.today()))"
i = 1
version = version_base
while version in tags
  global i, version, version_base
  i = i + 1
  version = version_base*".$i"
end
println("Will create version: $version")

# Update version file:

version_func = "
function version()
  version = \"$version\"
  println(\"Version: \$version\")
end
"

file = open("./version.jl","w")
write(file,version_func)
close(file)
run(`git add -A`)
run(`git commit -m "updated version file to $version"`)
run(`git tag -a $version -m "Release $version"`)
run(`git push origin master tag $version`)

range = "$(tags[length(tags)])...$version"
tagdiff = split(read(`git log --pretty=oneline $range`,String),'\n')

println("----------------------")
println("CHANGE LOG:")
println("----------------------")
for line in tagdiff
  if length(line) < 1
    continue
  end
  num = split(line)[1]
  line = replace(line, num => "-")
  println(line)
end
println("----------------------")


text:
# This file was generated, do not modify it. # hide
m11 = machine(model1, corestrict(X, f, 1), corestrict(y, f, 1))
m12 = machine(model1, corestrict(X, f, 2), corestrict(y, f, 2))
m13 = machine(model1, corestrict(X, f, 3), corestrict(y, f, 3))