text:
@testset "Filter" begin
    # Defining our missing datasets
    a = allowmissing(1.0:1.0:20.0)
    a[[2, 3, 7]] .= missing
    mask = map(!ismissing, a)

    # We call collect to not have a wrapper type that references the same data.
    m = collect(reshape(a, 5, 4))
    m_no_missing = reshape(1.0:1.0:20.0, 5, 4)

    aa = AxisArray(
        deepcopy(m),
        Axis{:time}(DateTime(2017, 6, 5, 5):Hour(1):DateTime(2017, 6, 5, 9)),
        Axis{:id}(1:4)
    )

    table = DataFrame(
        :sin => allowmissing(sin.(1.0:1.0:20.0)),
        :cos => allowmissing(sin.(1.0:1.0:20.0)),
    )

    table.sin[[2, 3, 7, 12, 19]] .= missing

    @test Filter() == Filter()

    @testset "Vector" begin
        result = apply(a, Filter())
        expected = deleteat!(deepcopy(a), [2, 3, 7])

        @test result == expected
        @test result == Impute.filter(a)

        a2 = deepcopy(a)
        Impute.filter!(a2)
        @test a2 == expected
    end

    @testset "Matrix" begin
        @test_throws UndefKeywordError apply(m, Filter())
        @test_throws UndefKeywordError Impute.filter(m)
        @test_throws MethodError Impute.filter!(m)

        @testset "no missing" begin
            result = apply(m_no_missing, Filter())
            @test isequal(result, m_no_missing)

            result = apply(m_no_missing, Filter(); kw=:foo)  # kwargs unused
            @test isequal(result, m_no_missing)
        end

        @testset "rows" begin
            result = apply(m, Filter(); dims=:rows)
            expected = m[[1, 4, 5], :]

            @test isequal(result, expected)
            @test isequal(result, (Impute.filter(m; dims=:rows)))
            @test isequal(collect(result'), Impute.filter(collect(m'); dims=:cols))
        end

        @testset "cols" begin
            result = apply(m, Filter(); dims=:cols)
            expected = copy(m)[:, 3:4]

            @test isequal(result, expected)
            @test isequal(result, Impute.filter(m; dims=:cols))
            @test isequal(collect(result'), Impute.filter(collect(m'); dims=:rows))
        end
    end

    @testset "Tables" begin
        @testset "DataFrame" begin
            df = deepcopy(table)

            @test_throws UndefKeywordError apply(df, Filter())
            @test_throws UndefKeywordError Impute.filter(df)
            @test_throws MethodError Impute.filter!(df)

            @testset "rows" begin
                result = apply(df, Filter(); dims=:rows)
                expected = dropmissing(df)

                @test isequal(result, expected)
                @test isequal(result, Impute.filter(df; dims=:rows))
            end

            @testset "cols" begin
                result = apply(df, Filter(); dims=:cols)
                expected = select(df, :cos)

                @test isequal(result, expected)
                @test isequal(result, Impute.filter(df; dims=:cols))
            end
        end

        @testset "Column Table" begin
            coltab = Tables.columntable(table)

            @test_throws UndefKeywordError apply(coltab, Filter())
            @test_throws UndefKeywordError Impute.filter(coltab)
            @test_throws MethodError Impute.filter!(coltab)

            @testset "rows" begin
                result = apply(coltab, Filter(); dims=:rows)
                expected = Tables.columntable(dropmissing(table))

                @test isequal(result, expected)
                @test isequal(result, Impute.filter(coltab; dims=:rows))
            end

            @testset "cols" begin
                result = apply(coltab, Filter(); dims=:cols)
                expected = Tables.columntable(TableOperations.select(coltab, :cos))

                @test isequal(result, expected)
                @test isequal(result, Impute.filter(coltab; dims=:cols))
            end
        end

        @testset "Row Table" begin
            @testset "rows" begin
                rowtab = Tables.rowtable(table)

                result = apply(rowtab, Filter(); dims=:rows)
                expected = Tables.rowtable(dropmissing(table))

                @test isequal(result, expected)
                @test isequal(result, Impute.filter(rowtab; dims=:rows))
                @test isequal(result, Impute.filter(rowtab))

                rowtab_ = Impute.filter!(rowtab)
                @test isequal(rowtab, expected)
                @test isequal(rowtab_, expected)
            end

            @testset "cols" begin
                rowtab = Tables.rowtable(table)

                @test_throws ArgumentError Impute.filter!(rowtab; dims=:cols)
                result = apply(rowtab, Filter(); dims=:cols)
                expected = Tables.rowtable(TableOperations.select(rowtab, :cos))

                @test isequal(result, expected)
                @test isequal(result, Impute.filter(rowtab; dims=:cols))
            end

        end
    end

    @testset "AxisArray" begin
        @test_throws UndefKeywordError apply(aa, Filter())
        @test_throws UndefKeywordError Impute.filter(aa)
        @test_throws MethodError Impute.filter!(aa)

        @testset "rows" begin
            result = apply(aa, Filter(); dims=:rows)
            expected = aa[[1, 4, 5], :]

            @test isa(result, AxisArray)
            @test isequal(result, expected)
            @test isequal(result, Impute.filter(aa; dims=:rows))
            @test isequal(collect(result'), Impute.filter(collect(aa'); dims=:cols))
        end

        @testset "cols" begin
            result = apply(aa, Filter(); dims=:cols)
            expected = copy(aa)[:, 3:4]

            @test isa(result, AxisArray)
            @test isequal(result, expected)
            @test isequal(result, Impute.filter(aa; dims=:cols))
            @test isequal(collect(result'), Impute.filter(collect(aa'); dims=:rows))
        end
    end

    @testset "functional" begin
        expected = deleteat!(deepcopy(a), [2, 3, 7])
        @test Impute.filter(!ismissing, a) == expected

        b = deepcopy(a)
        @test Impute.filter!(!ismissing, b) == expected
        @test b == expected
    end
end


text:
using Documenter, Cuas
push!(LOAD_PATH,"../src/")
makedocs(sitename="My Documentation")

makedocs(
         sitename = "Cuas.jl",
         modules  = [Cuas],
         pages=[
                "Home" => "index.md"
               ])
deploydocs(;
    repo="github.com/YayeIrene/Cuas.jl",
)

ext: jl


text:
using Test
using ScottishTaxBenefitModel
using .Utils
using .Definitions

@enum A a b c d e

@testset "tests of miscellaneous functions" begin
    gr = Dict{A,Real}( a=>1.0, b=>2, e=> 4.1 )
    ne = Dict{A,Real}( c=>2, b=>9 )
    inc = Dict{A,Real}( a=>1, b=>1.0, c=>1, d =>1 )
    m = mult( data=gr, calculated=ne, included=inc )
    @test m == 12

    pid = 120181906404
    p = from_pid( pid )
    @test p.datasource == FRS
    @test p.year == 2018
    @test p.hid == 19064
    @test p.pno == 4

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module RocketLambdaActorTest

using Test
using Rocket

@testset "LambdaActor" begin

    println("Testing: actor LambdaActor")

    @testset begin
        values = []
        source = from(1:3)

        n = (d) -> push!(values, d)
        e = (e) -> push!(values, e)
        c = ()  -> push!(values, "completed")

        actor  = LambdaActor{Int, typeof(n), typeof(e), typeof(c)}(n, e, c)

        subscribe!(source, actor)

        @test values == [ 1, 2, 3, "completed" ]
    end

    @testset begin
        values = []
        source = from(1:3)
        actor  = lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        )

        subscribe!(source, actor)

        @test values == [ 1, 2, 3, "completed" ]
    end

    @testset begin
        values = []
        actor  = lambda(Any,
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        )

        next!(actor, 1)

        @test values == [ 1 ]
    end

    @testset begin
        values = []
        actor  = lambda(Any,
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        )

        error!(actor, 'e')

        @test values == [ 'e' ]
    end

    @testset begin
        values = []
        actor  = lambda(Any,
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        )

        complete!(actor)

        @test values == [ "completed" ]
    end

    @testset begin
        n = (d) -> println(d)
        e = (e) -> println(e)
        c = ( ) -> println(0)

        @test lambda(Int; on_next     = n) isa LambdaActor{Int, typeof(n), Nothing, Nothing}
        @test lambda(Int; on_error    = e) isa LambdaActor{Int, Nothing, typeof(e), Nothing}
        @test lambda(Int; on_complete = c) isa LambdaActor{Int, Nothing, Nothing, typeof(c)}
        @test lambda(Int; on_next     = n, on_error    = e) isa LambdaActor{Int, typeof(n), typeof(e), Nothing  }
        @test lambda(Int; on_error    = e, on_complete = c) isa LambdaActor{Int, Nothing,   typeof(e), typeof(c)}
        @test lambda(Int; on_complete = c, on_next     = n) isa LambdaActor{Int, typeof(n), Nothing,   typeof(c)}
        @test lambda(Int; on_next     = n, on_error    = e, on_complete = c) isa LambdaActor{Int, typeof(n), typeof(e), typeof(c)}

        @test lambda(on_next     = n) isa Rocket.LambdaActorFactory{typeof(n), Nothing, Nothing}
        @test lambda(on_error    = e) isa Rocket.LambdaActorFactory{Nothing, typeof(e), Nothing}
        @test lambda(on_complete = c) isa Rocket.LambdaActorFactory{Nothing, Nothing, typeof(c)}
        @test lambda(on_next     = n, on_error    = e) isa Rocket.LambdaActorFactory{typeof(n), typeof(e), Nothing}
        @test lambda(on_error    = e, on_complete = c) isa Rocket.LambdaActorFactory{Nothing, typeof(e), typeof(c)}
        @test lambda(on_complete = c, on_next     = n) isa Rocket.LambdaActorFactory{typeof(n), Nothing, typeof(c)}
        @test lambda(on_next     = n, on_error    = e, on_complete = c) isa Rocket.LambdaActorFactory{typeof(n), typeof(e), typeof(c)}
    end
end

end


text:
type RODEProblem{uType,tType,isinplace,NP,F,C,MM,ND} <: AbstractRODEProblem{uType,tType,isinplace,ND}
  f::F
  u0::uType
  tspan::Tuple{tType,tType}
  noise::NP
  callback::C
  mass_matrix::MM
  rand_prototype::ND
end

function RODEProblem(f,u0,tspan; iip = isinplace(f,4),
                      rand_prototype = nothing,
                      noise= nothing,
                     callback=nothing,mass_matrix=I)
  RODEProblem{typeof(u0),eltype(tspan),iip,typeof(noise),typeof(f),typeof(callback),typeof(mass_matrix),typeof(rand_prototype)}(f,u0,tspan,noise,callback,mass_matrix,rand_prototype)
end


text:
module RocketFindOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: find()" begin

    println("Testing: operator find()")

    run_testset([
        (
            source      = from(0:5) |> find(d -> d !== 0 && d % 2 == 0),
            values      = @ts([ 2, c ]),
            source_type = Int
        ),
        (
            source      = of(1) |> find(d -> d !== 0 && d % 2 == 0),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = of("Something") |> find(d -> true),
            values      = @ts([ "Something", c ]),
            source_type = String
        ),
        (
            source      = completed() |> find(d -> d !== 0 && d % 2 == 0),
            values      = @ts(c)
        ),
        (
            source      = completed() |> find(d -> true),
            values      = @ts(c)
        ),
        (
            source      = never() |> find(d -> true),
            values      = @ts()
        )
    ])

end

end


text:
using Base.Test

using DataFrames

import Celeste: GalsimBenchmark

@testset "GalSim benchmark tests" begin
    results = GalsimBenchmark.main(
        test_case_names=["simple_star", "angle_and_axis_ratio_1"],
        print_fn=x -> 0,
    )
    for row_index in 1:size(results, 1)
        row = results[row_index, :]
        if isna(row[1, :expected])
            continue
        end
        for result_column in [:single_infer_actual, :joint_infer_actual]
            @test_approx_eq_eps(
                results[1, :expected],
                results[1, result_column],
                0.05 * results[1, :expected], # maximum permissible error
            )
        end
    end
end


text:
#
# Correctness Tests
#

using Base.Test
# Stealing your sweet testrunner, thanks DataFrames.
# using DataFrames
using Parsimonious

my_tests = [
    "test/test_benchmarks.jl",
    "test/test_expressions.jl",
    # test_firmware.jl
    "test/test_grammars.jl",
    "test/test_nodes.jl"
]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
using StatsModelComparisons, StanSample
using StatsFuns, RDatasets

df = RDatasets.dataset("datasets", "cars")

cars_stan = "
data {
    int N;
    vector[N] speed;
    vector[N] dist;
}
parameters {
    real a;
    real b;
    real sigma;
}
transformed parameters{
    vector[N] mu;
    mu = a + b * speed;
}
model {
    a ~ normal(0, 100);
    b ~ normal(0, 10);
    sigma ~ exponential(1);
    dist ~ normal(mu, sigma)    ;
}
generated quantities {
    vector[N] log_lik;
    for (i in 1:N)
        log_lik[i] = normal_lpdf(dist[i] | mu[i], sigma);
}
"

cars_stan_model = SampleModel("cars.model", cars_stan)
data = (N = size(df, 1), speed = df.Speed, dist = df.Dist)
rc = stan_sample(cars_stan_model; data)

if success(rc)
    stan_summary(cars_stan_model, true)
    nt_cars = read_samples(cars_stan_model);
end

log_lik = nt_cars.log_lik'
n_sam, n_obs = size(log_lik)
lppd = reshape(logsumexp(log_lik .- log(n_sam); dims=1), n_obs)
pwaic = [var(log_lik[:, i]) for i in 1:n_obs]
-2(sum(lppd) - sum(pwaic)) |> display

waic(log_lik) |> display


text:
# AbstractSimulator

"Simulator on random/fixed trees"
abstract type AbstractSimulator end

"""
    print_simulation(sim::AbstractSimulator, nsample)

Print `nsample` samples in simulation `sim`.
"""
function print_simulation(sim::AbstractSimulator, nsample)
    printfmtln("# {} -- {} sample(s)", sim, nsample)
    #printfmtln("# Seed of random number generator {}", Random.GLOBAL_RNG.seed)

    @sync @distributed for i = 1:nsample
        #host = gethostname()

        #print("$host ")
        for r = simulation(sim)
            printfmt("{} ", r)
        end
        println()
        if i % 10 == 0
            flush(stdout)
        end
    end
end

"""
    simulation(sim::AbstractSimulator, nsample)

Return `nsample` samples in simulation `sim`.
"""
function simulation(sim::AbstractSimulator, nsample)
    results = Vector{Any}(undef, nsample)
    for i in 1:nsample
        results[i] = simulation(sim)
    end
    results
end

# KcutSimulator

"""
Simulation of the [k-cut number](https://arxiv.org/abs/1804.03069) of trees.
"""
struct KcutSimulator <: AbstractSimulator
    tree::FiniteTree
    k::Int
end

function simulation(sim::KcutSimulator)
    walker = KcutWalker(sim.k)
    walk(sim.tree, walker)
    result(walker)
end

show(io::IO, sim::KcutSimulator) = printfmt(io, "{}-cut simulation on {}", sim.k, sim.tree)

# SubtreeSizeSimulator

"""
Simulation of the additive function in the form of sum of func(subtree size) applied to random trees.
"""
struct SubtreeSizeSimulator{F<:Function} <: AbstractSimulator
    tree::FiniteTree
    func::F
    funcname::String
end

"""
    simulation(sim::SubtreeSizeSimulator)

Return one sample in simulation `sim`.
"""
function simulation(sim::SubtreeSizeSimulator)
    # simulate the random tree and get subtree sizes
    walker = SubtreeSizeWalker(sim.tree)
    walk(sim.tree, walker)
    subtree_sizes = result(walker)

    # sum over log of subtree sizes
    ret = 0
    for size in subtree_sizes
        @fastmath ret += sim.func(size)
        #@fastmath ret += log(size)
    end
    ret
end

show(io::IO, sim::SubtreeSizeSimulator) = printfmt(io, "sum of {} simulation on {}",
                                                   sim.funcname, sim.tree)
# LogProductSimulator

"""
    LogProductSimulator(tree)

Construct a `SubtreeSizeSimulator` that simulate the sum of log(subtree size) of trees.
"""
LogProductSimulator(tree) = LogProductSimulator(tree, 1)

"""
    LogProductSimulator(tree, power)

Construct a `SubtreeSizeSimulator` that simulate the sum of log(subtree size)^`power` of trees.
"""
function LogProductSimulator(tree::FiniteTree, power::Int)
    function func(size::Int)
        @fastmath log(size)^power
    end
    return SubtreeSizeSimulator(tree, func, format("log(subtree size)^{}", power))
end

# SizePowerSimulator

"""
    SizePowerSimulator(tree)

Construct a `SubtreeSizeSimulator` that simulate the sum of sqrt(subtree size) of trees.
"""
SizePowerSimulator(tree) = SizePowerSimulator(tree, 1/2)

"""
    SizePowerSimulator(tree, power)

Construct a `SubtreeSizeSimulator` that simulate the sum of (subtree size)^`power` of trees.
"""
function SizePowerSimulator(tree::FiniteTree, power)
    function func(size::Int)
        @fastmath size^power
    end
    return SubtreeSizeSimulator(tree, func, format("(subtree size)^{}", power))
end

"""
Simulation of the [total path length](https://arxiv.org/abs/1102.2541) of trees.
"""
struct TotalPathSimulator <: AbstractSimulator
    tree::FiniteTree
end

function simulation(sim::TotalPathSimulator)
    # simulate the tree and count subtree sizes
    walker = SubtreeSizeWalker(sim.tree)
    walk(sim.tree, walker)

    sum(result(walker)[2:end])
end

show(io::IO, sim::TotalPathSimulator) = printfmt(io, "total path length simulation of {}", sim.tree)

# HeightSimulator

"""
Simulation of height of trees.
"""
struct HeightSimulator <: AbstractSimulator
    tree::FiniteTree
end

function simulation(sim::HeightSimulator)
    walker = DepthWalker(sim.tree)
    walk(sim.tree, walker)
    maximum(result(walker))
end

show(io::IO, sim::HeightSimulator) = printfmt(io, "height simulation of {}", sim.tree)

# LeafSimulator

"""
Simulation of the number of leafs.
"""
struct LeafSimulator <: AbstractSimulator
    tree::FiniteTree
end

function simulation(sim::LeafSimulator )
    degseq = degrees(sim.tree)
    leafnum = 0
    for degree in degseq
        if degree == 0
            leafnum += 1
        end
    end
    leafnum
end

show(io::IO, sim::LeafSimulator) = printfmt(io, "leaf number simulation of {}", sim.tree)


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

module RPI

include("../extimports.jl")

include("../aceimports.jl")

import ACE.SphericalHarmonics: SHBasis, index_y

export RPIBasis, SparsePSHDegree, BasicPSH1pBasis,
       diagonal_regulariser

# specify the `Rotations3D` submodule for CLebsch-Gordan and related
include("rotations3d.jl")
using ACE.RPI.Rotations3D

# some basic degree types useful for RPI type constructions
# (this file also specifies the PSH1pBasisFcn
include("rpi_degrees.jl")

# the basic RPI type 1-particle basis
include("rpi_basic1pbasis.jl")
include("rpi_1pbasis.jl")

# RPI basis and RPI potential
# (the RPI potential is specified through the combine function in this file)
include("rpi_basis.jl")

end


text:
"""
    plot(C::ClimGrid, titlefig::String, gridfig::Bool, label::String, color, lw, linestyle)

Plots the spatial average timeserie of ClimGrid `C`.
"""
function PyPlot.plot(C::ClimGrid; level=1, titlestr::String="", gridfig::Bool=true, label::String="", lw=1.5, linestyle="-", xlimit=[], ylimit=[], filename::String="")

    data = C[1].data
    timevec = get_timevec(C)

    if typeof(timevec[1]) != Date
        if typeof(timevec[1]) == Dates.Year
            timevec = Date.(timevec)
        end
    end

    average = fill(NaN, length(timevec))

    # Spatial mean for each timestep
    for t in 1:length(timevec)
        if ndims(data) == 3
            datatmp = data[:, :, t]
            # average[t] = Images.meanfinite(datatmp, 1:2)[1]
            average[t] = Statistics.mean(datatmp[.!isnan.(datatmp)])
        elseif ndims(data) == 4
            datatmp = data[:, :, level, t]
            average[t] = Statistics.mean(datatmp[.!isnan.(datatmp)])
        end
    end

    # figh, ax = subplots()

    if isempty(label)
        label = C.model
    end

    # Convert timevec to an array of string
    timevec_str = string.(timevec)
    if length(timevec) >= 20
        nb_interval_tmp = length(timevec)/8
        nb_int = ClimatePlots.roundup(nb_interval_tmp, 5)
    else
        nb_int = 1
    end

    # PLOTTING
    figh = PyPlot.plot(1:length(timevec_str), average, lw=lw, label=label, linestyle=linestyle)
    xticks(1:nb_int:length(timevec_str), timevec_str[1:nb_int:end], rotation=10)
    if !isempty(xlimit)
        xlim(xlimit[1], xlimit[2])
    end
    if !isempty(ylimit)
        ylim(ylimit[1], ylimit[2])
    end
    # xlabel("Time")
    ylabel(C.dataunits)
    legend()
    # xticks(timevec[1:20:end])
    if isempty(titlestr)
        titlestr = C.variable
    end
    title(titlestr)
    if gridfig
        grid(true)
    end

    # Save to "filename" if not empty
    if !isempty(filename)
        PyPlot.savefig(filename, dpi=300)
    end

    return true, figh

end

"""
    hist(C::ClimGrid; bins::Int=10, level=1, range_x=[], poly=[], start_date::Tuple=(Inf,), end_date::Tuple=(Inf,), titlestr::String="", gridfig::Bool=true, label::String="", ylimit=[])
"""
function PyPlot.hist(C::ClimGrid; bins::Int=10, level=1, range_x=[], titlestr::String="", gridfig::Bool=true, label::String="", ylimit=[])

    if isempty(label)
        label = C.model
    end

    # PLOTTING
    if gridfig
        grid(true)
    end
    if isempty(range_x)
        range_x = (findmin(C, skipnan=true)[1], findmax(C, skipnan=true)[1])
    end

    figh = hist(C[1].data[:], bins=bins, label=label, range=range_x)

    if !isempty(ylimit)
        ylim(ylimit[1], ylimit[2])
    end
    # xlabel("Time")
    ylabel("Frequency")
    legend()
    # xticks(timevec[1:20:end])
    if isempty(titlestr)
        titlestr = C.variable
    end
    title(titlestr)

    return true, figh

end


text:
struct OrdinalFactor{T<:Unsigned} <: AbstractFactor{T}
    name::String
    levels::AbstractVector{<:AbstractString}
    basefactor::AbstractFactor{T}
    newindex::Vector{T}
    islessfun::Union{Function, Nothing}
end

Base.length(factor::OrdinalFactor{T}) where {T<:Unsigned} = length(factor.basefactor)

function OrdinalFactor(name::String, basefactor::AbstractFactor{T}, islessfun::Function) where {T<:Unsigned}
    baselevels = getlevels(basefactor)
    levelcount = length(baselevels)
    perm = sortperm(baselevels, lt = islessfun)
    sortlevels = baselevels[perm]
    newindex = Vector{T}(undef, levelcount + 1)
    newindex[1] = 0
    for i in 1:levelcount
        newindex[i + 1] = perm[i]
    end
    OrdinalFactor{T}(name, sortlevels, basefactor, newindex, islessfun)
end

function OrdinalFactor(factor::AbstractFactor{T}) where {T<:Unsigned}
    OrdinalFactor{T}(getname(factor), getlevels(factor), factor, Vector{T}(), nothing)
end

function slice(factor::OrdinalFactor{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned}
    if factor.islessfun === nothing
        slice(factor.basefactor, fromobs, toobs, slicelength)
    else
        newindex = factor.newindex
        f = (i::T -> newindex[i + one(T)])
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        slices = slice(factor.basefactor, fromobs, toobs, slicelength)
        mapslice(f, slices, slicelength, T)
    end
end

function isordinal(factor::OrdinalFactor{T}) where {T<:Unsigned}
   true
end

function Base.map(factor::OrdinalFactor{T}, dataframe::AbstractDataFrame) where {T<:Unsigned}
    if factor.islessfun === nothing
        OrdinalFactor(map(factor.basefactor, dataframe))
    else
        OrdinalFactor(factor.name, map(factor.basefactor, dataframe), factor.islessfun)
    end
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# The quaternions form an associative normed division algebra

module FundamentalTests
    using Test: @test
    using Quaternionic

    # Algebra
    ## Vector space
    function test_vector_promotion(a, v::Quaternion)
        @test a + v == Quaternion(a) + v
        @test v + a == v + Quaternion(a)
        @test a - v == Quaternion(a) - v
        @test v - a == v - Quaternion(a)
    end
    test_vector_associativity(u::Quaternion, v::Quaternion, w::Quaternion) = @test u + (v + w) ≈ (u + v) + w rtol=eps(v)
    test_vector_commutativity(u::Quaternion, v::Quaternion) = @test u + v ≈ v + u rtol=eps(v)
    function test_vector_identity(v::Quaternion)
        @test v + zero(v) == v
        @test zero(v) + v == v
    end
    test_vector_inverse(v::Quaternion) = @test v + (-v) == zero(v)
    function test_vector_scalar_identity(v::Quaternion)
        @test one(eltype(v)) * v == v
        @test v * one(eltype(v)) == v
    end
    function test_vector_scalar_multiplication(a, v::Quaternion)
        @test a * v == Quaternion(a*v.w, a*v.x, a*v.y, a*v.z)
        @test v * a == Quaternion(a*v.w, a*v.x, a*v.y, a*v.z)
    end
    function test_vector_scalar_division(a, v::Quaternion)
        if !iszero(a)
            @test v / a == Quaternion(v.w / a, v.x / a, v.y / a, v.z / a)
        end
        if !iszero(v)
            @test a / v ≈ a * conj(v) / abs2(v) rtol=eps(v)
        end
    end
    test_vector_scalar_associativity(a, b, v::Quaternion) = @test a*(b*v) ≈ (a*b)*v rtol=eps(v)
    test_vector_scalar_distributivity1(a, u::Quaternion, v::Quaternion) = @test a*(u+v) ≈ (a*u) + (a*v) rtol=eps(v)
    test_vector_scalar_distributivity2(a, b, v::Quaternion) = @test (a+b)*v ≈ (a*v) + (b*v) rtol=eps(v)
    ## Product
    test_right_distributivity(x::Quaternion, y::Quaternion, z::Quaternion) = @test (x+y)*z ≈ (x*z) + (y*z) rtol=eps(x)
    test_left_distributivity(x::Quaternion, y::Quaternion, z::Quaternion) = @test z*(x+y) ≈ (z*x) + (z*y) rtol=eps(x)
    test_scalar_compatibility(a, b, x::Quaternion, y::Quaternion) = @test (a*x) * (b*y) ≈ (a*b) * (x*y) rtol=eps(x)

    # Division
    function test_identity(v::Quaternion)
        @test one(v) * v == v
        @test v * one(v) == v
    end
    function test_inverse(v::Quaternion)
        if !iszero(v)
            @test v * inv(v) == one(v)
            @test inv(v) * v == one(v)
        end
    end
    test_involution(x::Quaternion) = @test conj(conj(x)) == x
    test_involution_norm_real(x::Quaternion) = @test real(x * conj(x)) ≈ abs2(x) rtol=eps(x)
    test_involution_norm_imag(x::Quaternion) = @test absvec(x * conj(x)) ≈ 0 atol=eps(abs(x))

    # Normed
    test_norm_maps_to_field(q::Quaternion) = @test typeof(abs2(q)) == eltype(q)
    test_norm_nondegenerate(v::Quaternion) = @test iszero(v) ⊻ !iszero(abs2(v))
    test_norm_quadratic(v::Quaternion) = @test abs2(2*v) ≈ 4*abs2(v) rtol=eps(v)
    test_norm_quadratic(a, v::Quaternion) = @test abs2(a*v) ≈ a^2*abs2(v) rtol=eps(v)
    test_norm_composition(x::Quaternion, y::Quaternion) = @test abs2(x*y) ≈ abs2(x)*abs2(y) rtol=eps(x)

    # Associative
    test_associativity(u::Quaternion, v::Quaternion, w::Quaternion) = @test (u * v) * w ≈ u * (v * w) rtol=eps(v)
end

@testset verbose=true "Fundamentals" begin
    @testset "$T" for T in [FloatTypes...; IntTypes...]
        println("    Testing fundamentals with type $T")

        # Construct a variety of arguments
        scalars = [zero(T), one(T), -one(T)]#, 2*one(T), -2*one(T)]#, eps(T), -eps(T)]
        quaternions = [Quaternion(a, b, c, d) for a in scalars for b in scalars for c in scalars for d in scalars]

        # Iterate over all tests above
        for n in names(FundamentalTests, all=true)
            if T <: Integer && (n == :test_inverse)
                continue  # Don't try to invert integer quaternions
            end
            f = getproperty(FundamentalTests, n)
            if !isempty(methods(f)) && startswith(string(f), "test_")
                types = methods(f).ms[1].sig.parameters[2:end]
                args = Base.Iterators.product([type===Any ? scalars : quaternions for type in types]...)
                for arg in args
                    f(arg...)
                end
            end
        end
    end
    @testset "$T" for T in SymbolicTypes
        println("    Testing fundamentals with type $T")
        chars = Iterators.Stateful(Iterators.cycle("abcdefghijkl"))
        function next_scalar!(chars)
            x = Symbol(popfirst!(chars))
            xvar = @variables $x
            xvar[1]
        end
        function next_quaternion!(chars)
            x = Symbol(popfirst!(chars))
            # Work around <https://github.com/JuliaSymbolics/Symbolics.jl/issues/379>:
            # xvar = @variables $x[1:4]
            xvar = [[Symbolics.variable(:x, i) for i in 1:4]]
            Quaternion(xvar[1]...)
        end

        # Iterate over all tests above
        for n in names(FundamentalTests, all=true)
            f = getproperty(FundamentalTests, n)
            if !isempty(methods(f)) && startswith(string(f), "test_")
                types = methods(f).ms[1].sig.parameters[2:end]
                args = [type===Any ? next_scalar!(chars) : next_quaternion!(chars) for type in types]
                f(args...)
            end
        end
    end
end


text:
# # A time-varying PDE with forcing

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
In a previous example of the Poisson equation we demonstrated the use of point forcing. In this example
we will explore other ways to apply forcing to a PDE. Our target
problem will be similar to the most recent example -- transient heat conduction --
but now it will be in an unbounded domain. Our forcing will comprise local
*area heating* and *line heating*. We will also include convection by an
externally-imposed velocity field.

The governing equations are

$$\dfrac{\partial T}{\partial t} + \mathbf{v}\cdot\nabla T = \kappa \nabla^2 T + q$$

where $\mathbf{v}$ is a known convection velocity field and $q$ is a known heating field (scaled
by the density and specific heat)

In our discrete formulation, the problem takes the form

$$\mathcal{L}_C^\kappa T = -N(\mathbf{v},T) + q$$

where $N$ is a discrete approximation for the convection term $\mathbf{v}\cdot\nabla T$,
and $\mathcal{L}_C^\kappa = \mathrm{d}/\mathrm{d}t - \kappa L_C$ is the same
diffusion operator as in the previous example. Note that these equations are
no longer constrained. We only need to supply the right-hand side function.
Our job here is to provide everything needed to compute this right-hand side.

We will apply the heating and cooling inside of local regions.
As we discussed before, for each forcing, we supply information about the **shape**
of the forcing region and the **model function**, supplying the forcing strength
whenever we need it. We package these together into a forcing model and
provide it to the problem via the `forcing` keyword.

The convection velocity also must be provided. We also make use of the `forcing`
keyword to supply this to the problem, as well, in the form of a function
that returns the current convection velocity at a given time `t`. The
convective derivative term `N` requires a special cache, which we also
generate in the `prob_cache`.

We will highlight these aspects in the example that follows. In this example,
we will create
* a square-shaped line forcing region that delivers a certain amount of heat per unit length of line
* a circular-shaped area forcing region that supplies heat based on the heat-transfer model $h\cdot(T_0\cos(2\pi f t)-T)$, where $h$ is the heat transfer coefficient, $T$ is the current local temperature,
and $T_0$ is the amplitude of an oscillatory heater temperature, oscillating
at frequency $f$
* a background rotational velocity field, with angular velocity $\Omega$.
=#

using ImmersedLayers
#!jl using Plots
using UnPack

#=
### Set up the grid
Let's set up the grid first before we go any further
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx);

#=
### Specify the physical parameters
We supply here the diffusivity ($\kappa$) and the uniform convection velocity.
We also supply the grid Fourier number and a CFL number, which we will make
use of later to calculate the time step size.
=#
phys_params = Dict("diffusivity" => 0.005,
                    "angular velocity" => 0.5,
                    "length scale" => Lx/2,
                     "Fourier" => 1.0,
                     "CFL" => 0.5,
                     "lineheater_flux" => -2.0,
                     "areaheater_freq" => 1.0,
                     "areaheater_temp" => 2.0,
                     "areaheater_coeff" => 100.0
                     )

#=
## Specifying the forcing regions and models and convection velocity model
We first create the square line heater and place it at $(0,1)$. Its
associated model function is very simple, since it just sets the
strength uniformly. But note the function signature, which must always take this
form. It can make use of the current temperature, time, and physical parameters,
to return the strength of the line forcing. We bundle these together
using [`LineForcingModel`](@ref).
=#
fregion1 = Square(0.5,1.4*Δx)
T = RigidTransform((0.0,1.0),0.0)
T(fregion1)

function model1!(σ,T,t,fr::LineRegionCache,phys_params)
    σ .= phys_params["lineheater_flux"]
end
lfm = LineForcingModel(fregion1,model1!);

#=
Now the oscillatory heater, which we place at $(0,-0.5)$. This one has a few more parameters,
since we provide the heat transfer coefficient and the amplitude and frequency
of the target temperature. These are bundled with [`AreaForcingModel`](@ref).
=#

fregion2 = Circle(0.2,1.4*Δx)
T = RigidTransform((0.0,-0.5),0.0)
T(fregion2)

function model2!(σ,T,t,fr::AreaRegionCache,phys_params)
    f = phys_params["areaheater_freq"]
    T0 = phys_params["areaheater_temp"]
    h = phys_params["areaheater_coeff"]
    σ .= h*(T0*cos(2π*f*t) - T)
end
afm = AreaForcingModel(fregion2,model2!);

#=
Finally, the convection velocity model. Here, we make use of the
coordinate function [`x_gridgrad`](@ref) and [`y_gridgrad`](@ref)
to supply the coordinates of the velocity grid points. Since this
is a staggered grid, the velocity components live at different places.
For example, `yg.u` denotes the $y$ coordinates for the horizontal velocity
components, and `xg.v` denotes the $x$ coordinates for the vertical
velocity components.
=#
function my_velocity!(vel,t,cache,phys_params)
    xg, yg = x_gridgrad(cache), y_gridgrad(cache)
    Ω = phys_params["angular velocity"]
    vel.u .= -Ω*yg.u
    vel.v .= Ω*xg.v
    return vel
end

#=
We pack the forcing and convection together into the *forcing* `Dict`:
=#
forcing_dict = Dict("heating models" => [lfm,afm],
                    "convection velocity model" => my_velocity!)

#=
## Construct the ODE function and the extra cache
For the RHS of the heat conduction equation, we calculate the convective
derivative and the external heating.
=#
function heatconduction_rhs!(dT,T,sys::ILMSystem,t)
    @unpack forcing, phys_params, extra_cache, base_cache = sys
    @unpack cdcache, fcache, v, dT_tmp = extra_cache

    ## This provides the convection velocity at time `t`
    forcing["convection velocity model"](v,t,base_cache,phys_params)

    ## Compute the convective derivative term `N(v,T)`
    fill!(dT_tmp,0.0)
    convective_derivative!(dT_tmp,v,T,base_cache,cdcache)
    dT .= -dT_tmp

    ## Compute the contribution from the forcing models to the right-hand side
    fill!(dT_tmp,0.0)
    apply_forcing!(dT_tmp,T,t,fcache,phys_params)
    dT .+= dT_tmp

    return dT
end

#=
We create a problem type for this, define the extra cache, and extend `prob_cache`
=#
@ilmproblem UnboundedHeatConduction scalar

struct UnboundedHeatConductionCache{VT,CDT,FRT,DTT,FT} <: AbstractExtraILMCache
   v :: VT
   cdcache :: CDT
   fcache :: FRT
   dT_tmp :: DTT
   f :: FT
end

function ImmersedLayers.prob_cache(prob::UnboundedHeatConductionProblem,
                                   base_cache::BasicILMCache{N,scaling}) where {N,scaling}
    @unpack phys_params, forcing = prob
    @unpack gdata_cache, g = base_cache

    ## Construct a Lapacian outfitted with the diffusivity
    κ = phys_params["diffusivity"]
    heat_L = Laplacian(base_cache,gdata_cache,κ)

    ## Create cache for the convective derivative
    v = zeros_gridgrad(base_cache)
    cdcache = ConvectiveDerivativeCache(base_cache)

    ## Create cache for the forcing regions
    fcache = ForcingModelAndRegion(forcing["heating models"],base_cache);

    dT_tmp = zeros_grid(base_cache)

    ## The state here is temperature, and we just supply the RHS function
    f = ODEFunctionList(state = zeros_grid(base_cache),
                        ode_rhs=heatconduction_rhs!,
                        lin_op=heat_L)

    UnboundedHeatConductionCache(v,cdcache,fcache,dT_tmp,f)
end

#=
The last definition we need is for a timestep function. This time,
we take into account both the Fourier and the CFL conditions:
=#
function timestep_fourier_cfl(g,phys_params)
    κ = phys_params["diffusivity"]
    Ω = phys_params["angular velocity"]
    L = phys_params["length scale"]
    Fo = phys_params["Fourier"]
    Co = phys_params["CFL"]
    Δt = min(Fo*cellsize(g)^2/κ,Co*cellsize(g)/Ω/L)
    return Δt
end

#=
As in the last example, we also define the temperature output function for output
=#
temperature(T,σ,sys::ILMSystem,t) = T
@snapshotoutput temperature

#=
## Set up the problem and system
This is similar to previous problems.
=#
prob = UnboundedHeatConductionProblem(g,scaling=GridScaling,
                                        phys_params=phys_params,
                                        forcing=forcing_dict,
                                        timestep_func=timestep_fourier_cfl)

sys = construct_system(prob);

#=
## Solve the problem
As before, we first initialize the state, then we create an integrator,
and finally, advance the solution in time
=#
u0 = init_sol(sys)
tspan = (0.0,2.0)
integrator = init(u0,tspan,sys)

#=
Run the problem for one time unit
=#
step!(integrator,1.0)

#=
Let's see what this looks like. We will plot a set of snapshots in an array.
=#
sol = integrator.sol

#!jl plt = plot(layout = (2,3), size = (700, 500), legend=:false)
#!jl for (i,t) in enumerate(0:0.2:1.0)
#!jl    plot!(plt[i],temperature(sol,sys,t),sys,levels=range(-10,2,length=30),clim=(-10,2))
#!jl end
#!jl plt


text:
module NumericalMethodsforEngineers

using IterativeSolvers
using Reexport
@reexport using Symata, SparseArrays, LinearAlgebra, Compat

#using BandedMatrices

# package code goes here
### Imports ###
#import Symata: @sym, setsymata, getsymata, symeval

### Includes ###
# Type definitions in nmlib
include("nmlib/SkylineMatrices.jl")
include("nmlib/SymmetricBandedMatrices.jl")
# Direct methods in nmlib
include("nmlib/checkit.jl")
include("nmlib/nmex.jl")
include("nmlib/chobac.jl")
include("nmlib/cholin.jl")
include("nmlib/ldlfor.jl")
include("nmlib/ldlt.jl")
include("nmlib/lufac.jl")
include("nmlib/sparin.jl")
include("nmlib/spabac.jl")
include("nmlib/subbac.jl")
include("nmlib/subfor.jl")
include("nmlib/bicgstab.jl")
# Iterative methods in ch02
#include("ch02", "bicgstabl.jl"))
# Initial Value methods
include("ch05/lagrangianpolynomial.jl")
# Initial Value methods
include("ch07/ivp.jl")
# Boundary Value methods
include("ch07/bvp.jl")
# Symata functions
include("ch07/sym_fcns.jl")

### Exports ###

export
  # nmlib (types)
  SkylineMatrix,
  fromskyline,
  SymmetricBandedMatrix,
  tosymmetricbandedmatrix,
  # nmlib (functions)
  nmex,
  checkit,
  bicgstab,
  chobac!,
  cholin!,
  ldlt!,
  ldltfac!,
  ldlfor!,
  lufac,
  sparin!,
  spabac!,
  subbac!,
  subfor!,
  # ch05
  lagrangianpolynomial,
  # ch07
  euler,
  modified_euler,
  mid_point_euler,
  runga_kutta_4,
  shootingmethod,
  bigcstab,

  # Symata
  @sym,
  setsymata,
  getsymata,
  symeval

### Deprecated ###

  include("deprecated.jl")

end # module


text:
module ExcelWriter

using DataFrames
using PyCall

global const pyModnumpy = PyCall.PyNULL()
global const pyModPandas = PyCall.PyNULL()
global const pyModxlsxwriter = PyCall.PyNULL()

function __init__()
	#the following lines may trigger the installation of the respective python packages
	copy!(pyModnumpy, PyCall.pyimport_conda("numpy","numpy"))
	copy!(pyModPandas, PyCall.pyimport_conda("pandas","pandas"))
	copy!(pyModxlsxwriter, PyCall.pyimport_conda("xlsxwriter","xlsxwriter"))
end

include("excel.jl")

end # module


text:
using CorpusLoaders
using MultiResolutionIterators
using TextAnalysis
using TextAnalysis: CRF, crf_loss
using WordTokenizers
using Embeddings
using Flux
using Flux: onehot, param, onehotbatch, Conv, LSTM, flip, Momentum, reset!
using Tracker
using BSON: @save

CHAR_EMBED_DIMS = 25
WORD_EMBED_DIMS = 50
DESIRED_TAG_SCHEME = "BIOES"

train_set = load(CoNLL(), "train") # training set
# test_set = load(CoNLL(), "test") # test set
# dev_set = load(CoNLL(), "dev") # dev set

dataset = flatten_levels(train_set, lvls(CoNLL, :document)) |> full_consolidate
# dev_dataset = flatten_levels(dev_set, lvls(CoNLL, :document)) |> full_consolidate

X_train = [CorpusLoaders.word.(sent) for sent in dataset]
Y_train = [CorpusLoaders.named_entity.(sent) for sent in dataset]
tag_scheme!.(Y_train, "BIO2", DESIRED_TAG_SCHEME)

@assert length.(X_train) == length.(Y_train)

# X_dev = [CorpusLoaders.word.(sent) for sent in test_dataset]
# Y_dev = [CorpusLoaders.named_entity.(sent) for sent in test_dataset]
# tag_scheme!.(Y_dev, "BIO2", DESIRED_TAG_SCHEME)
#
# @assert length.(X_dev) == length.(Y_dev)

words_vocab = unique(vcat(X_train...))
alphabets = unique(vcat(collect.(words_vocab)...))
labels = unique(vcat(Y_train...))
DENSE_OUT_SIZE = num_labels = length(labels)

# TODO: Preprocessing: Change n't => not and Shuffle and all numbers to 0

#### Word Embeddings
# using these indices to generate the flux one hot input word embeddings vectors.

# TODO: EMBEDDINGS OF the input size.
embtable = load_embeddings(GloVe)
get_word_index = Dict(word => ii for (ii, word) in enumerate(embtable.vocab))
get_word_from_index = Dict(value => key for (key, value) in get_word_index)
W_word_Embed = (embtable.embeddings)

# One Vec for unknown chars
UNK_Word = "<UNK>"
UNK_Word_Idx = length(get_word_index) + 1
embedding_vocab_length = length(get_word_index) + 1
@assert UNK_Word ∉ collect(keys(get_word_index))

get_word_index[UNK_Word] = UNK_Word_Idx
get_word_from_index[UNK_Word_Idx] = UNK_Word
W_word_Embed = hcat(W_word_Embed, rand(WORD_EMBED_DIMS))
W_word_Embed = Float32.(W_word_Embed)
@assert size(W_word_Embed, 2) == embedding_vocab_length

##### Char Embeddings
UNK_char = '¿'
@assert UNK_char ∉ alphabets
push!(alphabets, UNK_char)
sort!(unique!(alphabets))

# using these indices to generate the flux one hot input Char embeddings vectors.

W_Char_Embed = rand(CHAR_EMBED_DIMS, length(alphabets))
W_Char_Embed = Float32.(W_Char_Embed)
W_Char_Embed = param(W_Char_Embed)
get_char_index = Dict(char => ii for (ii, char) in enumerate(alphabets))
get_char_from_index = Dict(value => key for (key, value) in get_char_index)

############# Creating onehoot vectors (useful for embeddings lookup), convenient
# TODO: Minibatches
onehotword(word) = onehot(get(get_word_index, lowercase(word), get_word_index[UNK_Word]), 1:embedding_vocab_length)
onehotchars(word) = onehotbatch([get(get_char_index, c, get_char_index[UNK_char]) for c in word], 1:length(alphabets))
onehotlabel(label) = onehot(label, labels)
onehotinput(word) = (onehot(get(get_word_index, lowercase(word), get_word_index[UNK_Word]), 1:embedding_vocab_length),
                onehotbatch([get(get_char_index, c, get_char_index[UNK_char]) for c in word], 1:length(alphabets)))

oh_seq(arr, f) = [f(element) for element in arr]

# X_words_train = [oh_seq(sentence, onehotword) for sentence in X_train] # A Bunch of Sequences of words, i.e. sentences
# X_chars_train = [oh_seq(sentence, onehotchars) for sentence in X_train] # A Bunch Sequences of Array of Chars, done to prevent repeated computations.
X_input_train = [oh_seq(sentence, onehotinput) for sentence in X_train]
Y_oh_train = [oh_seq(tags_sequence, onehotlabel) for tags_sequence in Y_train]

# X_input_dev = [oh_seq(sentence, onehotinput) for sentence in X_dev]
# Y_oh_dev = [oh_seq(tags_sequence, onehotlabel) for tags_sequence in Y_dev]

@assert length.(X_train) == length.(X_input_train) == length.(Y_oh_train)

#################### MODEL ######################
# 1. Character Embeddings with CNNs
# 2. Obtain Word Embeddings, concatenate Char & Word Embeddings
# 3. Bi-LSTM
# 4. Softmax or CRF on each unit in sequence

CNN_OUTPUT_SIZE = 53
CONV_WINDOW_LENGTH = 3
LSTM_STATE_SIZE = 253
DROPOUT_RATE_CNN = DROPOUT_INPUT_EMBEDDING = DROPOUT_OUT_LAYER = 0.68

# 1. Dropout before conv, Max poll layer, PADDING on both sides, hyperparams = window size, output vector size.
char_features = Chain(x -> W_Char_Embed * x,
                      x -> reshape(x, size(x)..., 1,1),
                      Dropout(DROPOUT_RATE_CNN),
                      Conv((CHAR_EMBED_DIMS, CONV_WINDOW_LENGTH), 1=>CNN_OUTPUT_SIZE, pad=(0,2)),
                      x -> maximum(x, dims=2),
                      x -> reshape(x, length(x),1))

# 2. Input embeddings:
# Maybe could use only the embeddings for the words needed

W_word_Embed = param(W_word_Embed) # For trainable
get_word_embedding(w) = W_word_Embed * w # works coz - onehot

# Dropout before LSTM
dropout_embed = Dropout(DROPOUT_INPUT_EMBEDDING)
input_embeddings((w, cs)) = dropout_embed(vcat(get_word_embedding(w), char_features(cs)))

# 3. Bi-LSTM

forward_lstm = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE)
backward = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE)
backward_lstm(x) = flip(backward, x)

bilstm_layer(x) = vcat.(forward_lstm.(x), backward_lstm(x))

# TODO: Bias vectors are initialized to zero, except the bias bf for the forget gate in LSTM , which is initialized to 1.0

# 4. Softmax / CRF Layer
# Dropout after LSTM

dropout = Dropout(DROPOUT_OUT_LAYER)
m(w_cs) = dropout.(bilstm_layer(input_embeddings.(w_cs)))

using TextAnalysis: crf_loss, CRF
Flux.@treelike TextAnalysis.CRF
c = TextAnalysis.CRF(num_labels, LSTM_STATE_SIZE * 2) |> gpu

loss(w_cs, y) =  crf_loss(c, m(w_cs), y)

η = 0.01 #
β = 0.05 # rate decay
ρ = 0.9 # momentum

# TODO: rate decay
opt = Momentum(η, ρ)
data = zip(X_input_train, Y_oh_train)

ps = params(params(char_features)..., params(W_word_Embed)..., params(W_Char_Embed)..., params(forward_lstm, backward)..., params(c)...)

NUM_EPOCHS = 5

function save_weights(char_features, W_word_Embed, W_Char_Embed, forward_lstm, backward_lstm, c)
    char_f_cpu = char_features |> cpu
    W_word_cpu = W_word_Embed |> cpu
    W_char_cpu = W_Char_Embed |> cpu
    forward_lstm_cpu = forward_lstm |> cpu
    backward_lstm_cpu = backward_lstm |> cpu
    crf_cpu = c |> cpu

    @save "./weights/char_f_cpu.bson" char_f_cpu
    @save "./weights/W_word_cpu.bson" W_word_cpu
    @save "./weights/W_char_cpu.bson" W_char_cpu
    @save "./weights/forward_lstm.bson" forward_lstm_cpu
    @save "./weights/backward_lstm.bson" backward_lstm_cpu
    @save "./weights/crf.bson" crf_cpu
end

function train()
    i = 0
    reset!(forward_lstm)
    reset!(backward)

    for epoch in 1:NUM_EPOCHS
        println("----------------------- EPOCH : $epoch ----------------------")
        for d in data
            grads = Tracker.gradient(() -> loss(d[1], d[2]), ps)
            reset!(forward_lstm)
            reset!(backward)

            Flux.Optimise.update!(opt, ps, grads)

            if i % 1000 == 0
                save_weights(char_features, W_word_Embed, W_Char_Embed, forward_lstm, backward_lstm, c)
                println(sum([loss(dd[1], dd[2]) for dd in data])/length(data))
            end
            i += 1
        end
    end
end

# TODO: gradient clipping and rate decay
# batch size = 10

# function test_raw_sentence # Convert unknown to UNK and to lowercase
# TODO: Try with and without lowercased chars in char embedding
# TODO: Try with and without trainable embeddings.


text:
__precompile__(true)

module GaussianFreeFields

import Contour,
       Interpolations,
       FFTW,
       SparseArrays,
       LinearAlgebra

export DGFF,
       fix_boundary_values,
       harmonicextension,
       inside,
       flowline,
       interpolate

function DGFF(n::Int)
    h = complex(zeros(n,n))
    for j=1:n
        for k=1:n
            h[j,k] = 1/sqrt(2.0)*(j+k == 2 ? 0 : (randn() + im*randn()) * 1/sqrt(sin((j-1)*pi/n)^2+ sin((k-1)*pi/n)^2))
        end
    end
    return real(n*1/sqrt(2)*FFTW.ifft(h))
end

function DGFF(m::Int,n::Int)
    h = complex(zeros(m,n))
    for j=1:m
        for k=1:n
            h[j,k] = 1/sqrt(2.0)*(j+k == 2 ? 0 :
                                  randn() + im*randn() *
                                  1/sqrt(sin((j-1)*pi/m)^2+(sin((k-1)*pi/n))^2))
        end
    end
    return real(n/sqrt(2)*FFTW.ifft(h))
end

"Return the array obtained by subtracting from `h` the
harmonic extension of the values around the edge of
the array `boundary_values`"
function fix_boundary_values(h::Array{Float64,2},
                             boundary_values::Array{Float64,2}=zeros(h))
    m,n = size(h)
    #sparse matrix with ones in top left, lower right positions
    cornersonly(m) = SparseArrays.sparse([1,m],[1,m],[1.0,1.0])
    #sparse matrix with ones along diagonal except top-left, lower-right
    eyenocorners(n) = SparseArrays.sparse(LinearAlgebra.I,n,n) - cornersonly(n)
    A = cornersonly(m)
    B = SparseArrays.sparse(LinearAlgebra.I,n,n)
    D_dense = zeros(n,n);
    for i=1:n-1
        D_dense[i,i+1] = -1;
        D_dense[i+1,i] = -1;
        D_dense[i,i] = 4;
    end
    D_dense[1,1] = 1
    D_dense[1,2] = 0
    D_dense[n,n-1] = 0
    D_dense[n,n] = 1
    D = SparseArrays.sparse(D_dense);
    E_dense = zeros(m,m)
    for i=1:m-1
        E_dense[i,i+1]=-1
        E_dense[i+1,i]=-1
    end
    E_dense[m,m-1] = 0
    E_dense[1,2] = 0
    E = SparseArrays.sparse(E_dense)
    Δ = kron(A,B) + kron(eyenocorners(m),D) + kron(E,eyenocorners(n))
    boundary = zeros(m*n)
    for i in 1:m-1
        boundary[n*i+1] = h[i+1,1] - boundary_values[i+1,1]
        boundary[n*i] = h[i,n] - boundary_values[i,n]
    end
    for i in 1:n-1
        boundary[i] = h[1,i] - boundary_values[1,i]
        boundary[end-n+i] = h[m,i] - boundary_values[m,i]
    end
    boundary[end] = h[m,n]

    return h - transpose(reshape(Array(LinearAlgebra.lu(Δ) \ boundary),n,m))
end

function harmonicextension(h::Array{T,2},vertices::Set) where T<:Number
    m,n = size(h)
    function neighbors(i::Int,j::Int)
        return [(i-1,j),(i+1,j),(i,j-1),(i,j+1)]
    end
    vertexmap = Dict(map(reverse,enumerate([(i,j) for i=1:m,j=1:n][:])));
    I = Int64[]
    J = Int64[]
    V = Float64[]
    b = zeros(m*n)
    for v in keys(vertexmap)
        if v in vertices
            k = vertexmap[v]
            push!(I,k); push!(J,k); push!(V,1)
            b[k] = h[k]
        else
            k = vertexmap[v]
            for nb in neighbors(v...)
                try
                    j = vertexmap[nb]
                    push!(I,k); push!(J,j); push!(V,1)
                    push!(I,k); push!(J,k); push!(V,-1)
                catch
                end
            end
        end
    end
    A = SparseArrays.sparse(I,J,V)
    return reshape(LinearAlgebra.lu(A) \ b,m,n)
end

# `inside` determines whether a point p is inside a curve γ
# by checking the number of points of intersection between γ
# and a ray emanating from p
function inside(p::Array{Float64,1},γ::Array{Float64,2})
    if γ[1,:] != γ[length(γ[:,1]),:]
        return false
    end
    cntr = 0; m = sqrt(2); # the slope is an arbitrary irrational number
    for i=1:length(γ[:,1])-1
        (x1,y1,x2,y2) = (γ[i,1],γ[i,2],γ[i+1,1],γ[i+1,2])
        if ((y2 - p[2] - m*(x2-p[1]))*(y1 - p[2] - m*(x1-p[1])) < 0)
            if (m*p[1]*x1 - p[2]*x1 - m*p[1]*x2 + p[2]*x2 - x2*y1 + x1*y2) /
                (m*x1 - m*x2 - y1 + y2)  - p[1] > 0
                cntr += 1
            end
        end
    end
    return isodd(cntr)
end

import Interpolations.interpolate
function interpolate(h::Array{T,2}) where T <: Real
    return interpolate(h,
                       Interpolations.BSpline(Interpolations.Linear()),
                       Interpolations.OnGrid())
end

import Base.getindex
getindex(h::Interpolations.BSplineInterpolation,
         z0::Complex) = h[real(z0),imag(z0)]

"""
Find the flow line of angle θ in the field h.
In other words, solve the DE η'(t) = exp(ih(η(t))/χ + iθ)".
"""
function flowline(h::Interpolations.BSplineInterpolation,
                  z0::Complex,
                  χ::Real,
                  θ::Real;
                  δ::Real=0.01,
                  S::Set{Complex}=Set{Complex}())
    (a,b) = size(h.coefs)
    η = [z0]
    while 1.0 ≤ real(η[end]) ≤ a && 1.0 ≤ imag(η[end]) ≤ b
        push!(η, η[end] + δ * exp(im*h[η[end]]/χ + im*θ))
        w = η[end]
        for z in S
            if abs2(z-w) < 1e-3
                return η
            end
        end
    end
    return η
end

flowline(h::Array{Float64,2},
         z0::Complex{Float64},
         χ::AbstractFloat,
         θ::AbstractFloat;
         δ::AbstractFloat=0.01,
         S::Set{Complex{Float64}}=Set()) = flowline(interpolate(h),z0,χ,θ,δ,S)

end # module


text:
#=
Julia implementation of a Region Proposal Network (RPN)
=#

struct rpn
    anchors_per_loc
    anchor_stride
    depth
end


text:
using Statistics
using BenchmarkTools
using LinearAlgebra

# Looking into pointers from http://fa.bianp.net/blog/2019/evaluate_logistic/
#

σ(x) = 1 / (1 + exp(-x))

## Order 0
function f_naive(x, A, b)
    z = A*x
    tmp = σ.(z)
    @. tmp = - b * log(tmp) - (1 - b) * log(1 - tmp)
    return sum(tmp)/length(tmp)
end

## blog logistic
function logsig_naive(t)
    return log(1 / (1 + exp(-t)))
end

function logsig_log1pexp(t)
    if t < -33.3
        return t
    elseif t <= -18
        return t - exp(t)
    elseif t <= 37
        return -log1p(exp(-t))
    else
        return -exp(-t)
    end
end

function f_logsumexp(x, A, b, logsig)
    m = size(A, 1)

    Ax = A * x
    fval = 0.0
    @inbounds @simd for i in 1:m
        fval += -logsig(b[i] * Ax[i])
    end
    fval /= m

    return fval
end

### Order 1
function ∇f_current(x, A, b)
    m = size(A, 1)

    σbAx = σ.(-b .* (A * x))
    res = transpose(A) * (σbAx .* -b)

    res /= m
    return res
end

function expit_b(x, b)
    out = zeros(size(x))

    @inbounds for (i, xi) in enumerate(x)
        if xi < 0
            out[i] = ((1 - b[i]) * exp(xi) - b[i]) / (1 + exp(xi))
        else
            out[i] = ((1 - b[i]) - b[i] * exp(-xi)) / (1 + exp(-xi))
        end
    end
    return out
end

function ∇f_prop(x, A, b)
    """Computes the gradient of the logistic loss.

    Parameters
    ----------
    x: array-like, shape (n_features,)
        Coefficients

    A: array-like, shape (n_samples, n_features)
        Data matrix

    b: array-like, shape (n_samples,)
        Labels

    Returns
    -------
    grad: array-like, shape (n_features,)
    """
    z = A*x
    s = expit_b(z, b)
    return A'*s / size(A, 1)
end

function main()
    A = [1 1;]
    b = [1]
    x = [20, 20]

    # mode = [:fvals]
    mode = [:∇f]

    if :fvals in mode
        @show f_naive(x, A, b)
        @show f_logsumexp(x, A, b, logsig_naive)
        @show f_logsumexp(x, A, b, logsig_log1pexp)
    elseif :∇f in mode
        @show ∇f_current(x, A, b)
        @show ∇f_prop(x, A, b)
    end

    println("")
    println("")

    println("Larger test")
    n, m = 500, 400
    A = rand(m, n)
    x = rand(n)
    b = rand(m)

    if :fvals in mode
        @show f_naive(x, A, b)
        @show f_logsumexp(x, A, b, logsig_naive)
        @show f_logsumexp(x, A, b, logsig_log1pexp)

        @btime f_naive($x, $A, $b)
        @btime f_logsumexp($x, $A, $b, $logsig_naive)
        @btime f_logsumexp($x, $A, $b, $logsig_log1pexp)
    elseif :∇f in mode
        ∇f_cu = ∇f_current(x, A, b)
        ∇f_pr = ∇f_prop(x, A, b)

        @show norm(∇f_cu - ∇f_pr)

        @btime ∇f_current($x, $A, $b)
        @btime ∇f_prop($x, $A, $b)
    end

    return
end

main()


text:
export setUpWindow, PushToRow, FinalizeRow
# ======================================================================================
### Types exported from Graphics
# ======================================================================================
function setUpWindow(document::Page, w::Float64,h::Float64)

    n = document.children[1]
    n.shape = NBox()
     n.shape.color = [1,1,1]
     n.shape.padding = BoxOutline(0,0,0,0,0,0)
     n.shape.border = Border(0,0, 0,0, 0,0, "solid",[.5,.5,.5,1],[0,0,0,0])
     n.shape.left    = 0
     n.shape.top     = 0
     n.shape.width   = w # S hould be window width  - controls
     n.shape.height  = h # S hould be window height - controls

    document.fixed.rows = [] #Row(0, 0, w)
    fn = document.fixed
    fn.shape = NBox()
      fn.shape.color = [1,1,1]
      fn.shape.padding = BoxOutline(0,0,0,0,0,0)
      fn.shape.border = Border(0,0, 0,0, 0,0, "solid",[.5,.5,.5,1],[0,0,0,0])
      fn.shape.left    = 0
      fn.shape.top     = 0
      fn.shape.width   = w # S hould be window width  - controls
      fn.shape.height  = h # S hould be window height - controls
end
#======================================================================================#
#======================================================================================#
# Types exported from Graphics
#
# move objects up or down (within row space) depending on layout options.
# Be sure that the heights of all objects have been set.
# Float the floats
# The contents of a container may effect the container itself. Such as height; set that!
# Set any node heights that are % values.

getFlags(node) = isa(node, TextElement) ? node.parent.font.flags : node.shape.flags
#======================================================================================#
function FinalizeRow(row::Row)
    if !row.flags[RowFinalized]

        setVertHorizDisplacement(row)
        #..........................................
        # float LEFT: MoveNodeToLeft(row, index)
        #..........................................
        for i in 2:length(row.nodes)
            n = row.nodes[i]
            flags = getFlags(n)
            if flags[FloatLeft] && !flags[Absolute]
                w, h = getSize(n.shape) # we really only need width here.
                MoveAll(n, -(n.shape.left - row.nodes[1].shape.left) ,0.0) # move to row-start (left)
                for j in 1:(i-1) # move everything before N to the right by w
                  MoveAll(row.nodes[j],w,0.0)
                end
                MoveNodeToLeft(row, i)
            end
        end
        #..........................................
        # float RIGHT: MoveNodeToRight(row, index)
        #    Begin with this...
        #--|-1-| |-2-| |-FR-| |-4-|--empty-space--#
        #    Desired result...
        #--|-1-| |-2-| |-4-|--empty-space-- |-FR-|#
        #..........................................
        for i in length(row.nodes):-1:1
            n = row.nodes[i]
            flags = getFlags(n)
            if flags[FloatRight] && !flags[Absolute]
                w = getWidthfromContentLeft(n.shape)
                MoveAll(n, row.width-n.shape.left-w, 0.0) # move to row-start (left)
                for j in (i+1):length(row.nodes)
                    w, h = getSize(n.shape)
                    MoveAll(row.nodes[j], -w, 0.0)
                end
                wide,high = getSize(row.nodes[end].shape)
                n.shape.left = row.nodes[end].shape.left + wide + row.space
            end

        end # for loop

        # Mark row as finalized!
        row.flags[RowFinalized] = true
    end
    return row.top + row.height
end
# ======================================================================================
# vertical/hirizintal shift
# ======================================================================================
function setVertHorizDisplacement(row::Row)
    shiftAll, x_space = 0.0, 0.0

    for i in 1:length(row.nodes)
        width, height = getSize(row.nodes[i].shape)
        row.height = max(row.height, height)
    end

  for i in 1:length(row.nodes)
      node = row.nodes[i]
      flags = getFlags(node) #node.parent.font.flags
      if !flags[Absolute]   # flags[Absolute] && return row.top + row.height
              shape = node.shape #.re ference
              y_space = 0.0
              # set row height: this may need changed to take into account
              width, height = getSize(shape)
              #row.height = max(row.height, height)

              if flags[TextCenter]
                      x_space = row.space * .5
              elseif flags[TextRight]
                      x_space = row.space
              end
              row.space = 0

              if flags[AlignBase]
                  println("row.height: ", row.height)
                  println("Shape height: ", height)
                  println("row.height - height: ", row.height - height)
                  y_space = (row.height - height)
              end

              if flags[AlignMiddle]  &&  row.height > height
                      y_space = (row.height - height) *.5
              end

              if x_space > 0 || y_space > 0
                  MoveAll(row.nodes[i], x_space, y_space)
              end
      end # ! Absolute
  end
end
# ======================================================================================
# function rowSpaceToWidth(row, shape, width::Float64)
#     shape.width = width
#     row.space -= width
# end
# rowSpaceToWidth(row, shape) = rowSpaceToWidth(row, shape, row.space)

function PushToRow(document::Page, node, l::Float64,t::Float64,w::Float64) # height not needed
    parent = node.parent
    shape = node.shape
    width, height = getSize(shape)
    #rows = parent.rows
    row = getCreateLastRow(parent.rows, l, t, w) # return the last row (creatig one if nesesary).
    flags = node.shape.flags

    # Fixed /////////////////////////////////////////////////////
    # push to document instead of current parent! Fixed nodes' positions are
    # measured relative to the viewport, not any specific node.
    # TODO: I think Fixed nodes should be left in the tree where they belong but
    #       the document width and height should be fixed first (...not being updated).
    if flags[Fixed]
        row = getCreateLastRow(document.fixed.rows, l, t, w) # Not redundant: (attached to doc).
        setDistanceFromBounds(0,0,document.width, document.height,  shape )
        push!(row.nodes, node)
        return
    end
    # Absolute /////////////////////////////////////////////////////
    # The only way to set an Absolute position is to wait until
    # the parent's extents are calculated and set.
    if flags[Absolute]
        push!(row.nodes, node) # push to row but don't change any Row metrics
        return
    end
     # Inline /////////////////////////////////////////////////////
    if flags[DisplayInline]
        if flags[Marked]
            println("node.shape.width: ", node.shape.width)
            println("typeof(node.shape): ", typeof(node.shape))
            println("node.parent.shape.width: ", node.parent.shape.width)
        end
        shape.width = row.space
    end
    # InlineBlock /////////////////////////////////////////////////////
    if flags[DisplayInlineBlock] && shape.width < 1
        shape.width = parent.rows[end].space-width # parent.rows[end].width
    end
    # Block /////////////////////////////////////////////////////
    if flags[DisplayBlock]
        if length(row.nodes) > 0
            row = Row(parent)  #row = Row(rows, l, FinalizeRow(row), w)
        end
        shape.width < 1    &&    (shape.width = w - getPackingWidth(shape)) #- getWidth(shape) #(border.width + padding.width + margin.width)
        row.space = 0 # Make sure nothing else gets put on this row.
    elseif row.space < shape.width  # not enough space.. new row!
        row = Row(parent) # row = Row(rows, l, FinalizeRow(row), w)
    end

    feedRow(parent, node)
end
#======================================================================================#
#
#======================================================================================#

#======================================================================================#
#
#======================================================================================#

#======================================================================================#


text:
using DelayEmbeddings
using Test
using Random
using DelimitedFiles

println("\nTesting uzal_cost.jl...")
@testset "Uzal cost" begin
@testset "Random vectors" begin
    ## Check on random vector
    Random.seed!(1516578735)
    tr1 = randn(10000)
    tr2 = randn(5000)
    tr1 = Dataset(tr1)
    tr2 = Dataset(tr2)
    L = uzal_cost(tr1;
        Tw = 60, K= 3, w = 1, samplesize = 1.0,
        metric = Euclidean()
    )
    L2 = uzal_cost(tr2;
        Tw = 60, K= 3, w = 1, samplesize = 1.0,
        metric = Euclidean()
    )

    L_max = -0.06
    L_min = -0.08
    @test L_min < L < L_max
    @test L_min < L2 < L_max

    ## check on clean step function (0-distances)
    tr = zeros(10000)
    tr[1:5000] = zeros(5000)
    tr[5001:end] = ones(5000)
    tr = Dataset(tr)
    L = uzal_cost(tr;
        Tw = 60, K= 3, w = 1, samplesize = 1.0,
        metric = Euclidean()
    )
    @test isnan(L)

    ## check on noisy step function (should yield a proper value)
    Random.seed!(1516578735)
    tr = zeros(10000)
    tr[1:5000] = zeros(5000) .+ 0.001 .* randn(5000)
    tr[5001:end] = ones(5000) .+ 0.001 .* randn(5000)
    tr = Dataset(tr)
    L = uzal_cost(tr;
        Tw = 60, K= 3, w = 1, samplesize = 1.0,
        metric = Euclidean()
    )
    @test -1.6 < L < -1.5
end

@testset "Uzal local cost (Lorenz)" begin
    ## Test Lorenz example
    # For comparison reasons using Travis CI we carry out the integration on a UNIX
    # OS and save the resulting time series
    # See https://github.com/JuliaDynamics/JuliaDynamics for the storage of
    # the time series used for testing
    #
    # u0 = [0, 10.0, 0.0]
    # lo = Systems.lorenz(u0; σ=10, ρ=28, β=8/3)
    # tr = trajectory(lo, 100; Δt = 0.01, Ttr = 100)
    tr = readdlm(joinpath(tsfolder, "test_time_series_lorenz_standard_N_10000_multivariate.csv"))
    tr = Dataset(tr)

    # check local cost function output
    Tw = 60
    L = uzal_cost(tr;
        Tw = Tw, K= 3, w = 12, samplesize = 1.0,
        metric = Euclidean()
    )
    L_local= uzal_cost_local(tr; Tw = Tw, K= 3, samplesize = 1.0,w = 12, metric = Euclidean())
    @test length(L_local) == length(tr)-Tw
    @test maximum(L_local)>L
    @test minimum(L_local)<L
    @test -0.3 < L < -0.28

end

@testset "Roessler system" begin
    ## Test Roessler example as in Fig. 7 in the paper with internal data
    # For comparison reasons using Travis CI we carry out the integration on a UNIX
    # OS and save the resulting time series
    # See https://github.com/JuliaDynamics/JuliaDynamics for the storage of
    # the time series used for testing
    #
    # u0 = [1.0, 1.0, 1.0]
    # ro = Systems.roessler(u0; a=0.15, b = 0.2, c=10)
    # tr = trajectory(ro, 1000; Δt = 0.1, Ttr = 100)
    tr = readdlm(joinpath(tsfolder, "test_time_series_roessler_N_10000_multivariate.csv"))
    tr = Dataset(tr)

    x = tr[:, 1]

    # theiler window
    w  = 12
    # Time horizon
    Tw = 50
    # sample size
    SampleSize = .5
    # metric
    metric = Euclidean()
    # embedding dimension
    m = 3
    # maximum neighbours
    k_max = 5
    # number of considered tw-values
    tw_max = 15 # corresponds to a maximum of tw = 30

    # preallocation
    L = zeros(k_max,tw_max)
    tws = zeros(tw_max)

    for K = 1:k_max
        for i = 1:tw_max
            tws[i] = i*(m-1)
            Y = embed(x,m,i)
            L[K,i] = uzal_cost(Y; Tw=Tw, K=K, w=w, samplesize=SampleSize, metric=metric)
        end
    end

    tau_min = 7
    tau_max = 12

    min1_idx = sortperm(L[1,:])
    min1 = tws[min1_idx[1]]
    @test tau_min < min1 < tau_max
    L_min = -1.2
    L_max = -0.8
    @test L_min < L[1,min1_idx[1]] < L_max

    min2_idx = sortperm(L[2,:])
    min2 = tws[min2_idx[1]]
    @test tau_min < min2 < tau_max
    L_min = -0.79
    L_max = -0.75
    @test L_min < L[2,min2_idx[1]] < L_max

    min3_idx = sortperm(L[3,:])
    min3 = tws[min3_idx[1]]
    @test tau_min < min3 < tau_max
    L_min = -0.62
    L_max = -0.58
    @test L_min < L[3,min3_idx[1]] < L_max

    min4_idx = sortperm(L[4,:])
    min4 = tws[min4_idx[1]]
    @test tau_min < min4 < tau_max
    L_min = -0.51
    L_max = -0.47
    @test L_min < L[4,min4_idx[1]] < L_max

    min5_idx = sortperm(L[5,:])
    min5 = tws[min5_idx[1]]
    @test tau_min < min5 < tau_max
    L_min = -0.43
    L_max = -0.39
    @test L_min < L[5,min5_idx[1]] < L_max

    L_local = uzal_cost_local(tr;
        Tw = Tw, K= 3, w = 12, samplesize=1.0,
        metric = Chebyshev()
    )
    L = uzal_cost(tr;
        Tw = Tw, K= 3, w = 12,
        metric = Chebyshev(), samplesize=1.0
    )
    @test length(L_local) == length(tr)-Tw
    @test maximum(L_local) > L
    @test minimum(L_local) < L

    # # plot results using PyPlot
    # using PyPlot
    # pygui(true)
    # labels = ["k=1", "k=2", "k=3", "k=4", "k=5"]
    # figure()
    # plot(tws,L[1,:], linewidth = 2)
    # plot(tws,L[2,:], linewidth = 2)
    # plot(tws,L[3,:], linewidth = 2)
    # plot(tws,L[4,:], linewidth = 2)
    # plot(tws,L[5,:], linewidth = 2)
    # xlabel(L"$t_w$")
    # ylabel(L"$L_k$")
    # legend(labels)
    # #yscale("symlog")
    # title("Roessler System as in Fig. 7(b) in the Uzal Paper")
    # grid()

    # # Display local cost function for two different embeddings
    # tau_value_1 = 8
    # tau_value_2 = 20
    # m = 2
    #
    # # embedding in two dimensions
    # Y_1 = embed(x,m,tau_value_1)
    # Y_2 = embed(x,m,tau_value_2)
    #
    # # compute local cost functions
    # L_local_1= uzal_cost_local(Y_1;
    #     Tw = Tw, K= 3, w = 12, metric = Euclidean()
    # )
    # L1= uzal_cost(Y_1;
    #     Tw = Tw, K= 3, w = 12, metric = Euclidean(), samplesize=1.0
    # )
    #
    # L_local_2= uzal_cost_local(Y_2;
    #     Tw = Tw, K= 3, w = 12, metric = Euclidean()
    # )
    # L2= uzal_cost(Y_2;
    #     Tw = Tw, K= 3, w = 12, metric = Euclidean(), samplesize=1.0
    # )
    # # plot results using PyPlot
    # using PyPlot
    # pygui(true)
    #
    # x_val1 = Y_1[1:length(L_local_1),1]
    # y_val1 = Y_1[1:length(L_local_1),2]
    #
    # x_val2 = Y_2[1:length(L_local_2),1]
    # y_val2 = Y_2[1:length(L_local_2),2]
    #
    # markersize = 10
    #
    # figure()
    # subplot(1,2,1)
    # scatter(x_val1,y_val1,c=L_local_1,s=markersize)
    # plot(x_val1,y_val1,linewidth=0.1)
    # xlabel("x(t)")
    # ylabel("x(t+$tau_value_1)")
    # title("Local L for Roessler System first embedding cycle (τ=$tau_value_1)")
    # grid()
    #
    # subplot(1,2,2)
    # scatter(x_val2,y_val2,c=L_local_2,s=markersize)
    # plot(x_val2,y_val2,linewidth=0.1)
    # xlabel("x(t)")
    # ylabel("x(t+$tau_value_2)")
    # title("Local L for Roessler System first embedding cycle (τ=$tau_value_2)")
    # grid()
    # cbar = colorbar()
    # cbar.set_label("Local cost function")

end

end # Uzal cost testset


text:
using Documenter, PowerSystems
import DataStructures: OrderedDict
using Literate

# This is commented out because the output is not user-friendly. Deliberation on how to best
# communicate this information to users is ongoing.
#include(joinpath(@__DIR__, "src", "generate_validation_table.jl"))
include(joinpath(@__DIR__, "make_model_library.jl"))

pages = OrderedDict(
        "Welcome Page" => "index.md",
        "User Guide" => Any[
            "user_guide/installation.md",
            "user_guide/quick_start_guide.md",
            "user_guide/system.md",
        ],
        "Modeler" => Any[
            "modeler/parsing.md",
            "modeler/data.md",
            "modeler/example_dynamic_data.md"
        ],
        "Model Developer" => Any[
            "Extending Parsing" => "model_developer/extending_parsing.md",
            "Adding Types" => "model_developer/adding_types.md",
        ],
        "Developer" => Any[
            "Style Guide" => "developer/style.md",
            "Testing" => "developer/testing.md",
            "Logging" => "developer/logging.md",
        ],
        "Model Library" => Any[],
        "Internal API" => "api/internal.md"
)

pages["Model Library"] = make_model_library(
     categories = [
        Topology,
        StaticInjection,
        Service,
        Branch,
        PSY.DeviceParameter,
    ],
    exceptions = [PSY.DynamicComponent,
                  PSY.ActivePowerControl,
                  PSY.ReactivePowerControl,
                  PSY.DynamicBranch
                  ],
    manual_additions =
        Dict("Service" => ["Reserves" => "model_library/reserves.md"],
        "StaticInjection" => ["Regulation Device" => "model_library/regulation_device.md"],
        "DynamicInjection" => ["Dynamic Inverter" => "model_library/dynamic_inverter.md",
        "Dynamic Generator" => "model_library/dynamic_generator.md",
        ],
        "StaticInjection" => ["HybridSystem" => "model_library/hybrid_device.md"],
        "Branch" => ["Dynamic Lines" => "model_library/dynamic_branch.md"]
        )
)

# This code performs the automated addition of Literate - Generated Markdowns. The desired
# section name should be the name of the file for instance network_matrices.jl -> Network Matrices
julia_file_filter = x -> occursin(".jl", x)
folders = Dict(
    "User Guide" => filter(julia_file_filter, readdir("docs/src/user_guide")),
    "Model Library" => filter(julia_file_filter, readdir("docs/src/model_library")),
    "Modeler" => filter(julia_file_filter, readdir("docs/src/modeler")),
    "Model Developer" => filter(julia_file_filter, readdir("docs/src/model_developer")),
    "Developer" => filter(julia_file_filter, readdir("docs/src/developer")),
)

for (section, folder) in folders
    for file in folder
        section_folder_name = lowercase(replace(section, " " => "_"))
        outputdir = joinpath(pwd(), "docs", "src", "$section_folder_name")
        inputfile = joinpath("$section_folder_name", "$file")
        outputfile = string("generated_", replace("$file", ".jl" => ""))
        Literate.markdown(joinpath(pwd(), "docs", "src", inputfile),
                          outputdir;
                          name = outputfile,
                          credit = false,
                          execute=true)
        subsection = titlecase(replace(split(file, ".")[1], "_" => " "))
        push!(pages[section], ("$subsection" =>  joinpath("$section_folder_name", "$(outputfile).md")))
    end
end

makedocs(
    modules = [PowerSystems],
    format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_ACTIONS"),),
    sitename = "PowerSystems.jl",
    authors = "Jose Daniel Lara, Daniel Thom and Clayton Barrows",
    pages = Any[p for p in pages]
)

deploydocs(
    repo = "github.com/NREL-SIIP/PowerSystems.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    push_preview=true,
    versions = ["stable" => "v^", "v#.#"],
)


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# 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.

using Test
using Distributed
include(joinpath(Sys.BINDIR, "..", "share", "julia", "test", "testenv.jl"))
addprocs_with_testenv(3)
@test nworkers() == 3
@everywhere using Logging
for w in workers()
    # Do not log on worker nodes
    remotecall(()->global_logger(NullLogger()),w)
end

@everywhere using StochasticPrograms
@everywhere using JuMP
@everywhere using LinearAlgebra
@everywhere using GLPK
@everywhere using Ipopt

@everywhere using Distributions

@everywhere import StochasticPrograms: probability, expected

include("../decisions/decisions.jl")
TestDecisionVariable.run_dtests()
TestDecisionConstraint.run_dtests()
TestDecisionObjective.run_dtests()
include("../problems/problem_load.jl")
include("dfunctional_tests.jl")
include("dsolver_tests.jl")


text:
using PkgBenchmark: judge, target_result, baseline_result,
    writeresults, export_markdown

mkpath("tmp")
JUDGE_BASELINE = get(ENV, "JUDGE_BASELINE", "^HEAD")
results = judge("LyapunovExponents", JUDGE_BASELINE)
writeresults("tmp/benchmarks_target.json", target_result(results))
writeresults("tmp/benchmarks_baseline.json", baseline_result(results))
export_markdown("tmp/benchmarks_judge.md", results)
showall(results)
println()