text:
@testset "Production Problem" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    c = [4.0, 3.0]
    A1 = [2.0, 1.0, 1.0]
    A2 = [1.0, 2.0, 1.0]
    b1 = 4.0
    b2 = 4.0

    x = MOI.add_variables(optimizer, length(c))

    @test typeof(x[1]) == MOI.VariableIndex

    w, cw = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    cons1 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A1, [x[1], x[2], y]),
        0.0,
    )
    MOI.add_constraint(optimizer, cons1, MOI.LessThan(b1))

    cons2 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A2, [x[1], x[2], z]),
        0.0,
    )
    MOI.add_constraint(optimizer, cons2, MOI.LessThan(b2))

    @test cons1.terms[1].coefficient == 2
    @test POI.is_parameter_in_model(optimizer, cons2.terms[3].variable_index)

    obj_func = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([c[1], c[2], 3.0], [x[1], x[2], w]),
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    MOI.get(optimizer, MOI.TerminationStatus())

    MOI.get(optimizer, MOI.PrimalStatus())

    @test isapprox.(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        28 / 3,
        atol = ATOL,
    )
    @test isapprox.(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        4 / 3,
        atol = ATOL,
    )
    @test isapprox.(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        4 / 3,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -5 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        3.0,
        atol = ATOL,
    )

    MOI.get(optimizer, MOI.VariablePrimal(), w)
    MOI.get(optimizer, MOI.VariablePrimal(), y)
    MOI.get(optimizer, MOI.VariablePrimal(), z)

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(2.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 2.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) == 1.0

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 13.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [1.0, 1.0]

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -5 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        3.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(0.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) == 1.0
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [1.0, 1.0]

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-5.0))
    MOI.optimize!(optimizer)
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 12.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [3.0, 0.0]

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 0, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -4,
        atol = ATOL,
    )
end

@testset "Production Problem variation on parameters for duals" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    c = [4.0, 3.0]
    A1 = [2.0, 1.0, 3.0]
    A2 = [1.0, 2.0, 0.5]
    b1 = 4.0
    b2 = 4.0

    x = MOI.add_variables(optimizer, length(c))

    @test typeof(x[1]) == MOI.VariableIndex

    w, cw = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    cons1 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A1, [x[1], x[2], y]),
        0.0,
    )
    ci1 = MOI.add_constraint(optimizer, cons1, MOI.LessThan(b1))

    cons2 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A2, [x[1], x[2], z]),
        0.0,
    )
    ci2 = MOI.add_constraint(optimizer, cons2, MOI.LessThan(b2))

    @test cons1.terms[1].coefficient == 2
    @test POI.is_parameter_in_model(optimizer, cons2.terms[3].variable_index)

    obj_func = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([c[1], c[2], 2.0], [x[1], x[2], w]),
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    MOI.get(optimizer, MOI.TerminationStatus())

    MOI.get(optimizer, MOI.PrimalStatus())

    @test isapprox(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        28 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        4 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        4 / 3,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -5,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2 / 6,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        2.0,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        3 * MOI.get(optimizer, MOI.ConstraintDual(), ci1),
        atol = 1e-4,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        0.5 * MOI.get(optimizer, MOI.ConstraintDual(), ci2),
        atol = 1e-4,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(2.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -9.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -0.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        2.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(0.0))

    MOI.optimize!(optimizer)

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 3.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -9.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -0.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        2.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-5.0))
    MOI.optimize!(optimizer)
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 14.0, atol = ATOL)
    @test isapprox(
        MOI.get.(optimizer, MOI.VariablePrimal(), x),
        [3.5, 0.0],
        atol = ATOL,
    )

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 0, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2,
        atol = ATOL,
    )
end


text:
if VERSION <= v"1.1"
    isnothing(x) = x === nothing
end

function define_docstring()
    path = joinpath(@__DIR__, "..", "README.md")
    include_dependency(path)
    doc = read(path, String)
    # doc = replace(doc, r"^```julia"m => "```jldoctest README")
    @eval TestFunctionRunner $Base.@doc $doc TestFunctionRunner
end

ext: jl


text:
# # Machine Learning in Julia, JuliaCon2020

# A workshop introducing the machine learning toolbox
# [MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)

# ### Set-up

# The following instantiates a package environment and pre-loads some
# packages, to avoid delays later on.

DIR = @__DIR__
include(joinpath(DIR, "setup.jl"))

# Only evaluate the next cell if you're using an IDE (not a notebook)
# and it has issues displaying color/boldface REPL output:

color_off()

# ## General resources

# - [List of methods introduced in this tutorial](methods.md)
# - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
# - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
# - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
# - [Data Science Tutorials in Julia](https://alan-turing-institute.github.io/DataScienceTutorials.jl/)

# ## Contents

# ### Basic

# - [Part 1 - Data Representation](#part-1-data-representation)
# - [Part 2 - Selecting, Training and Evaluating Models](#part-2-selecting-training-and-evaluating-models)
# - [Part 3 - Transformers and Pipelines](#part-3-transformers-and-pipelines)
# - [Part 4 - Tuning Hyper-parameters](#part-4-tuning-hyper-parameters)
# - [Part 5 - Advanced model composition](#part-5-advanced-model-composition)
# - [Solutions to Exercises](#solutions-to-exercises)

#

# ## Part 1 - Data Representation

# > **Goals:**
# > 1. Learn how MLJ specifies it's data requirements using "scientific" types
# > 2. Understand the options for representing tabular data
# > 3. Learn how to inspect and fix the representation of data to meet MLJ requirements

# ### Scientific types

# To help you focus on the intended *purpose* or *interpretation* of
# data, MLJ models specify data requirements using *scientific types*,
# instead of machine types. An example of a scientific type is
# `OrderedFactor`. The other basic "scalar" scientific types are
# illustrated below:

# ![](assets/scitypes.png)

# A scientific type is an ordinary Julia type (so it can be used for
# method dispatch, for example) but it usually has no instances. The
# `scitype` function is used to articulate MLJ's convention about how
# different machine types will be interpreted by MLJ models:

using MLJ
scitype(3.141)

#-

time = [2.3, 4.5, 4.2, 1.8, 7.1]
scitype(time)

# To fix data which MLJ is interpreting incorrectly, we use the
# `coerce` method:

height = [185, 153, 163, 114, 180]
scitype(height)

#-

height = coerce(height, Continuous)

# Here's an example of data we would want interpreted as
# `OrderedFactor` but isn't:

exam_mark = ["rotten", "great", "bla",  missing, "great"]
scitype(exam_mark)

#-

exam_mark = coerce(exam_mark, OrderedFactor)

#-

levels(exam_mark)

# Use `levels!` to put the classes in the right order:

levels!(exam_mark, ["rotten", "bla", "great"])
exam_mark[1] < exam_mark[2]

# When sub-sampling, no levels are lost:

levels(exam_mark[1:2])

# **Note on binary data.** There is no separate scientific type for
# binary data. Binary data is `OrderedFactor{2}` or
# `Multiclass{2}`. If a binary measure like `truepositive` is a
# applied to `OrderedFactor{2}` then the "positive" class is assumed
# to appear *second* in the ordering. If such a measure is applied to
# `Multiclass{2}` data, a warning is issued. A single `OrderedFactor`
# can be coerced to a single `Continuous` variable, for models that
# require this, while a `Multiclass` variable can only be one-hot
# encoded.

# ### Two-dimensional data

# Whenever it makes sense, MLJ Models generally expect two-dimensional
# data to be *tabular*. All the tabular formats implementing the
# [Tables.jl API](https://juliadata.github.io/Tables.jl/stable/) (see
# this
# [list](https://github.com/JuliaData/Tables.jl/blob/master/INTEGRATIONS.md))
# have a scientific type of `Table` and can be used with such models.

# The simplest example of a table is the julia native *column
# table*, which is just a named tuple of equal-length vectors:

column_table = (h=height, e=exam_mark, t=time)

#-

scitype(column_table)

#-

# Notice the `Table{K}` type parameter `K` encodes the scientific
# types of the columns. (This is useful when comparing table scitypes
# with `<:`). To inspect the individual column scitypes, we use the
# `schema` method instead:

schema(column_table)

# Here are four other examples of tables:

row_table = [(a=1, b=3.4),
             (a=2, b=4.5),
             (a=3, b=5.6)]
schema(row_table)

#-

import DataFrames
df = DataFrames.DataFrame(column_table)

#-

schema(df)

#-

using CSV
file = CSV.File(joinpath(DIR, "data", "horse.csv"));
schema(file) # (triggers a file read)

# Most MLJ models do not accept matrix in lieu of a table, but you can
# wrap a matrix as a table:

matrix_table = MLJ.table(rand(2,3))
schema(matrix_table)

# The matrix is *not* copied, only wrapped.

# **Manipulating tabular data.** In this workshop we assume
# familiarity with some kind of tabular data container (although it is
# possible, in principle, to carry out the exercises without this.)
# For a quick start introduction to `DataFrames`, see [this
# tutorial](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/dataframe/)

# ### Fixing scientific types in tabular data

# To show how we can correct the scientific types of data in tables,
# we introduce a cleaned up version of the UCI Horse Colic Data Set
# (the cleaning work-flow is described
# [here](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/horse/#dealing_with_missing_values))

using CSV
file = CSV.File(joinpath(DIR, "data", "horse.csv"));
horse = DataFrames.DataFrame(file); # convert to data frame without copying columns
first(horse, 4)

#-

# From [the UCI
# docs](http://archive.ics.uci.edu/ml/datasets/Horse+Colic) we can
# surmise how each variable ought to be interpreted (a step in our
# work-flow that cannot reliably be left to the computer):

# variable                    | scientific type (interpretation)
# ----------------------------|-----------------------------------
# `:surgery`                  | Multiclass
# `:age`                      | Multiclass
# `:rectal_temperature`       | Continuous
# `:pulse`                    | Continuous
# `:respiratory_rate`         | Continuous
# `:temperature_extremities`  | OrderedFactor
# `:mucous_membranes`         | Multiclass
# `:capillary_refill_time`    | Multiclass
# `:pain`                     | OrderedFactor
# `:peristalsis`              | OrderedFactor
# `:abdominal_distension`     | OrderedFactor
# `:packed_cell_volume`       | Continuous
# `:total_protein`            | Continuous
# `:outcome`                  | Multiclass
# `:surgical_lesion`          | OrderedFactor
# `:cp_data`                  | Multiclass

# Let's see how MLJ will actually interpret the data, as it is
# currently encoded:

schema(horse)

# As a first correction step, we can get MLJ to "guess" the
# appropriate fix, using the `autotype` method:

autotype(horse)

#-

# Okay, this is not perfect, but a step in the right direction, which
# we implement like this:

coerce!(horse, autotype(horse));
schema(horse)

# All remaining `Count` data should be `Continuous`:

coerce!(horse, Count => Continuous);
schema(horse)

# We'll correct the remaining truant entries manually:

coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
schema(horse)

# ### Resources for Part 1
#
# - From the MLJ manual:
#    - [A preview of data type specification in
#   MLJ](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#A-preview-of-data-type-specification-in-MLJ-1)
#    - [Data containers and scientific types](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Data-containers-and-scientific-types-1)
#    - [Working with Categorical Data](https://alan-turing-institute.github.io/MLJ.jl/dev/working_with_categorical_data/)
# - [Summary](https://alan-turing-institute.github.io/MLJScientificTypes.jl/dev/#Summary-of-the-MLJ-convention-1) of the MLJ convention for representing scientific types
# - [MLJScientificTypes.jl](https://alan-turing-institute.github.io/MLJScientificTypes.jl/dev/)
# - From Data Science Tutorials:
#     - [Data interpretation: Scientific Types](https://alan-turing-institute.github.io/DataScienceTutorials.jl/data/scitype/)
#     - [Horse colic data](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/horse/)
# - [UCI Horse Colic Data Set](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)

# ### Exercises for Part 1

# #### Exercise 1

# Try to guess how each code snippet below will evaluate:

scitype(42)

#-

questions = ["who", "why", "what", "when"]
scitype(questions)

#-

elscitype(questions)

#-

t = (3.141, 42, "how")
scitype(t)

#-

A = rand(2, 3)

# -

scitype(A)

#-

elscitype(A)

#-

using SparseArrays
Asparse = sparse(A)

#-

scitype(Asparse)

#-

using CategoricalArrays
C1 = categorical(A)

#-

scitype(C1)

#-

elscitype(C1)

#-

C2 = categorical(A, ordered=true)
scitype(C2)

#-

v = [1, 2, missing, 4]
scitype(v)

#-

elscitype(v)

#-

scitype(v[1:2])

# Can you guess at the general behavior of
# `scitype` with respect to tuples, abstract arrays and missing
# values? The answers are
# [here](https://github.com/alan-turing-institute/ScientificTypes.jl#2-the-scitype-and-scitype-methods)
# (ignore "Property 1").

# #### Exercise 2

# Coerce the following vector to make MLJ recognize it as a vector of
# ordered factors (with an appropriate ordering):

quality = ["good", "poor", "poor", "excellent", missing, "good", "excellent"]

#-

# #### Exercise 3 (fixing scitypes in a table)

# Fix the scitypes for the [House Prices in King
# County](https://mlr3gallery.mlr-org.com/posts/2020-01-30-house-prices-in-king-county/)
# dataset:

file = CSV.File(joinpath(DIR, "data", "house.csv"));
house = DataFrames.DataFrame(file); # convert to data frame without copying columns
first(house, 4)

# (Two features in the original data set have been deemed uninformative
# and dropped, namely `:id` and `:date`. The original feature
# `:yr_renovated` has been replaced by the `Bool` feature `is_renovated`.)

#

# ## Part 2 - Selecting, Training and Evaluating Models

# > **Goals:**
# > 1. Search MLJ's database of model metadata to identify model candidates for a supervised learning task.
# > 2. Evaluate the performance of a model on a holdout set using basic `fit!`/`predict` work-flow.
# > 3. Inspect the outcomes of training and save these to a file.
# > 3. Evaluate performance using other resampling strategies, such as cross-validation, in one line, using `evaluate!`
# > 4. Plot a "learning curve", to inspect performance as a function of some model hyper-parameter, such as an iteration parameter

# The "Hello World!" of machine learning is to classify Fisher's
# famous iris data set. This time, we'll grab the data from
# [OpenML](https://www.openml.org):

iris = OpenML.load(61); # a row table
iris = DataFrames.DataFrame(iris);
first(iris, 4)

# **Main goal.** To build and evaluate models for predicting the
# `:class` variable, given the four remaining measurement variables.

# ### Step 1. Inspect and fix scientific types

schema(iris)

#-

coerce!(iris, :class => Multiclass);
schema(iris)

# ### Step 2. Split data into input and target parts

# Here's how we split the data into target and input features, which
# is needed for MLJ supervised models. We randomize the data at the
# same time:

y, X = unpack(iris, ==(:class), name->true; rng=123);
scitype(y)

# Do `?unpack` to learn more:

@doc unpack

# ### On searching for a model

# Here's how to see *all* models (not immediately useful):

all_models = models()

# Each entry contains metadata for a model whose defining code is not yet loaded:

meta = all_models[3]

#-

targetscitype = meta.target_scitype

#-

scitype(y) <: targetscitype

# So this model won't do. Let's  find all pure julia classifiers:

filter_julia_classifiers(meta) =
    AbstractVector{Finite} <: meta.target_scitype &&
    meta.is_pure_julia

models(filter_julia_classifiers)

# Find all models with "Classifier" in `name` (or `docstring`):

models("Classifier")

# Find all (supervised) models that match my data!

models(matching(X, y))

# ### Step 3. Select and instantiate a model

# To load the code defining a new model type we use the `@load` macro,
# which returns an *instance* of the type, with default
# hyperparameters:

model = @load NeuralNetworkClassifier

#-

info(model)

# In MLJ a *model* is just a struct containing hyper-parameters, and
# that's all. A model does not store *learned* parameters. Models are
# mutable:

model.epochs = 12

# And all models have a key-word constructor that works once `@load`
# has been performed:

NeuralNetworkClassifier(epochs=12) == model

# ### On fitting, predicting, and inspecting models

# In MLJ a model and training/validation data are typically bound
# together in a machine:

mach = machine(model, X, y)

# A machine stores *learned* parameters, among other things. We'll
# train this machine on 70% of the data and evaluate on a 30% holdout
# set. Let's start by dividing all row indices into `train` and `test`
# subsets:

train, test = partition(eachindex(y), 0.7)

# Now we can `fit!`...

fit!(mach, rows=train, verbosity=2)

# ... and `predict`:

predict(mach, rows=test)  # or `predict(mach, Xnew)`

# After training, one can inspect the learned parameters:

fitted_params(mach)

#-

# Everything else the user might be interested in is accessed from the
# training *report*:

report(mach)

# You save a machine like this:

MLJ.save("neural_net.jlso", mach)

# And retrieve it like this:

mach2 = machine("neural_net.jlso")
predict(mach2, X)[1:3]

# If you want to fit a retrieved model, you will need to bind some data to it:

mach3 = machine("neural_net.jlso", X, y)
fit!(mach3)

# Machines remember the last set of hyper-parameters used during fit,
# which, in the case of iterative models, allows for a warm restart of
# computations in the case that only the iteration parameter is
# increased:

model.epochs = model.epochs + 4
fit!(mach, rows=train, verbosity=2)

# By default (for this particular model) we can also increase
# `:learning_rate` without triggering a cold restart:

model.epochs = model.epochs + 4
model.optimiser.eta = 10*model.optimiser.eta
fit!(mach, rows=train, verbosity=2)

# However, change any other parameter and training will restart from
# scratch:

model.lambda = 0.001
fit!(mach, rows=train, verbosity=2)

# Let's train silently for a total of 50 epochs, and look at a prediction:

model.epochs = 50
fit!(mach, rows=train)
yhat = predict(mach, X[test,:]); # or predict(mach, rows=test)
yhat[1]

# What's going on here?

info(model).prediction_type

# **Important**:
# - In MLJ, a model that can predict probabilities (and not just point values) will do so by default. (These models have supertype `Probabilistic`, while point-estimate predictors have supertype `Deterministic`.)
# - For most probabilistic predictors, the predicted object is a `Distributions.Distribution` object, supporting the `Distributions.jl` [API](https://juliastats.org/Distributions.jl/latest/extends/#Create-a-Distribution-1) for such objects. In particular, the methods `rand`,  `pdf`, `mode`, `median` and `mean` will apply, where appropriate.

# So, to obtain the probability of "Iris-virginica" in the first test
# prediction, we do

pdf(yhat[1], "Iris-virginica")

# To get the most likely observation, we do

mode(yhat[1])

# These can be broadcast over multiple predictions in the usual way:

broadcast(pdf, yhat[1:4], "Iris-versicolor")

#-

mode.(yhat[1:4])

# Or, alternatively, you can use the `predict_mode` operation instead
# of `predict`:

predict_mode(mach, X[test,:])[1:4] # or predict_mode(mach, rows=test)[1:4]

# For a more conventional matrix of probabilities you can do this:

L = levels(y)
pdf(yhat, L)[1:4, :]

# However, in a typical MLJ work-flow, this is not as useful as you
# might imagine. In particular, all probabilistic performance measures
# in MLJ expect distribution objects in their first slot:

cross_entropy(yhat, y[test]) |> mean

# To apply a deterministic measure, we first need to obtain point-estimates:

misclassification_rate(mode.(yhat), y[test])

# We note in passing that there is also a search tool for measures
# analogous to `models`:

measures()

#-

measures(matching(y)) # experimental

# ### Step 4. Evaluate the model performance

# Naturally, MLJ provides boilerplate code for carrying out a model
# evaluation with a lot less fuss. Let's repeat the performance
# evaluation above and add an extra measure, `brier_score`:

evaluate!(mach, resampling=Holdout(fraction_train=0.7),
          measures=[cross_entropy, brier_score])

# Or applying cross-validation instead:

evaluate!(mach, resampling=CV(nfolds=6),
          measures=[cross_entropy, brier_score])

# Or, Monte Carlo cross-validation (cross-validation repeated
# randomized folds)

e = evaluate!(mach, resampling=CV(nfolds=6, rng=123),
              repeats=3,
              measures=[cross_entropy, brier_score])

# One can access the following properties of the output `e` of an
# evaluation: `measure`, `measurement`, `per_fold` (measurement for
# each fold) and `per_observation` (measurement per observation, if
# reported).

# We finally note that you can restrict the rows of observations from
# which train and test folds are drawn, by specifying `rows=...`. For
# example, imagining the last 30% of target observations are `missing`
# you might have a work-flow like this:

train, test = partition(eachindex(y), 0.7)
mach = machine(model, X, y)
evaluate!(mach, resampling=CV(nfolds=6),
          measures=[cross_entropy, brier_score],
          rows=train)     # cv estimate, resampling from `train`
fit!(mach, rows=train)    # re-train using all of `train` observations
predict(mach, rows=test); # and predict missing targets

# ### On learning curves

# Since our model is an iterative one, we might want to inspect the
# out-of-sample performance as a function of the iteration
# parameter. For this we can use the `learning_curve` function (which,
# incidentally can be applied to any model hyper-parameter). This
# starts by defining a one-dimensional range object for the parameter
# (more on this when we discuss tuning in Part 4):

r = range(model, :epochs, lower=1, upper=50, scale=:log)
curve = learning_curve(mach,
                       range=r,
                       resampling=Holdout(fraction_train=0.7), # (default)
                       measure=cross_entropy)

using Plots
pyplot(size=(490,300))
plt=plot(curve.parameter_values, curve.measurements)
xlabel!(plt, "epochs")
ylabel!(plt, "cross entropy on holdout set")
plt

#@vlplot(curve.parameter_values,curve.measurements)

# We will return to learning curves when we look at tuning in Part 4.

# ### Resources for Part 2

# - From the MLJ manual:
#     - [Getting Started](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/)
#     - [Model Search](https://alan-turing-institute.github.io/MLJ.jl/dev/model_search/)
#     - [Evaluating Performance](https://alan-turing-institute.github.io/MLJ.jl/dev/evaluating_model_performance/) (using `evaluate!`)
#     - [Learning Curves](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_curves/)
#     - [Performance Measures](https://alan-turing-institute.github.io/MLJ.jl/dev/performance_measures/) (loss functions, scores, etc)
# - From Data Science Tutorials:
#     - [Choosing and evaluating a model](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/choosing-a-model/)
#     - [Fit, predict, transform](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/fit-and-predict/)

# ### Exercises for Part 2

# #### Exercise 4

# (a) Identify all supervised MLJ models that can be applied (without
# type coercion or one-hot encoding) to a supervised learning problem
# with input features `X4` and target `y4` defined below:

import Distributions
poisson = Distributions.Poisson

age = 18 .+ 60*rand(10);
salary = coerce(rand([:small, :big, :huge], 10), OrderedFactor);
levels!(salary, [:small, :big, :huge]);
X4 = DataFrames.DataFrame(age=age, salary=salary)

n_devices(salary) = salary > :small ? rand(poisson(1.3)) : rand(poisson(2.9))
y4 = [n_devices(row.salary) for row in eachrow(X4)]

# (b) What models can be applied if you coerce the salary to a
# `Continuous` scitype?

# #### Exercise 5 (unpack)

# After evaluating the following ...

data = (a = [1, 2, 3, 4],
     b = rand(4),
     c = rand(4),
     d = coerce(["male", "female", "female", "male"], OrderedFactor));
pretty(data)

using Tables
y, X, w = unpack(data, ==(:a),
                 name -> elscitype(Tables.getcolumn(data, name)) == Continuous,
                 name -> true);

# ...attempt to guess the evaluations of the following:

y

#-

pretty(X)

#-

w

# #### Exercise 6 (first steps in modeling Horse Colic)

# (a) Suppose we want to use predict the `:outcome` variable in the
# Horse Colic study introduced in Part 1, based on the remaining
# variables that are `Continuous` (one-hot encoding categorical
# variables is discussed later in Part 3) *while ignoring the others*.
# Extract from the `horse` data set (defined in Part 1) appropriate
# input features `X` and target variable `y`. (Do not, however,
# randomize the observations.)

# (b) Create a 70:30 `train`/`test` split of the data and train a
# `LogisticClassifier` model, from the `MLJLinearModels` package, on
# the `train` rows. Use `lambda=100` and default values for the
# other hyper-parameters. (Although one would normally standardize
# (whiten) the continuous features for this model, do not do so here.)
# After training:

# - (i) Recalling that a logistic classifier (aka logistic regressor) is
#   a linear-based model learning a *vector* of coefficients for each
#   feature (one coefficient for each target class), use the
#   `fitted_params` method to find this vector of coefficients in the
#   case of the `:pulse` feature. (You can convert a vector of pairs `v =
#   [x1 => y1, x2 => y2, ...]` into a dictionary with `Dict(v)`.)

# - (ii) Evaluate the `cross_entropy` performance on the `test`
#   observations.

# - ☆(iii) In how many `test` observations does the predicted
#   probability of the observed class exceed 50%?

# - (iv) Find the `misclassification_rate` in the `test`
#   set. (*Hint.* As this measure is deterministic, you will either
#   need to broadcast `mode` or use `predict_mode` instead of
#   `predict`.)

# (c) Instead use a `RandomForestClassifier` model from the
#     `DecisionTree` package and:
#
# - (i) Generate an appropriate learning curve to convince yourself
#   that out-of-sample estimates of the `cross_entropy` loss do not
#   substantially improve for `n_trees > 50`. Use default values for
#   all other hyper-parameters, and feel free to use all available
#   data to generate the curve.

# - (ii) Fix `n_trees=90` and use `evaluate!` to obtain a 9-fold
#   cross-validation estimate of the `cross_entropy`, restricting
#   sub-sampling to the `train` observations.

# - (iii) Now use *all* available data but set
#   `resampling=Holdout(fraction_train=0.7)` to obtain a score you can
#   compare with the `KNNClassifier` in part (b)(iii). Which model is
#   better?

#

# ## Part 3 - Transformers and Pipelines

# ### Transformers

# Unsupervised models, which receive no target `y` during training,
# always have a `transform` operation. They sometimes also support an
# `inverse_transform` operation, with obvious meaning, and sometimes
# support a `predict` operation (see the clustering example discussed
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/#Transformers-that-also-predict-1)).
# Otherwise, they are handled much like supervised models.

# Here's a simple standardization example:

x = rand(100);
@show mean(x) std(x);

#-

model = Standardizer() # a built-in model
mach = machine(model, x)
fit!(mach)
x̂ = transform(mach, x);
@show mean(x̂) std(x̂);

# This particular model has an `inverse_transform`:

inverse_transform(mach, x̂) ≈ x

# ### Re-encoding the King County House data as continuous

# For further illustrations of transformers, let's re-encode *all* of the
# King County House input features (see [Ex
# 3](#ex-3-fixing-scitypes-in-a-table)) into a set of `Continuous`
# features. We do this with the `ContinuousEncoder` model, which, by
# default, will:

# - one-hot encode all `Multiclass` features
# - coerce all `OrderedFactor` features to `Continuous` ones
# - coerce all `Count` features to `Continuous` ones (there aren't any)
# - drop any remaining non-Continuous features (none of these either)

# First, we reload the data and fix the scitypes (Exercise 3):

file = CSV.File(joinpath(DIR, "data", "house.csv"));
house = DataFrames.DataFrame(file)
coerce!(house, autotype(file))
coerce!(house, Count => Continuous, :zipcode => Multiclass);
schema(house)

#-

y, X = unpack(house, ==(:price), name -> true, rng=123);

# Instantiate the unsupervised model (transformer):

encoder = ContinuousEncoder() # a built-in model; no need to @load it

# Bind the model to the data and fit!

mach = machine(encoder, X) |> fit!;

# Transform and inspect the result:

Xcont = transform(mach, X);
schema(Xcont)

# ### More transformers

# Here's how to list all of MLJ's unsupervised models:

models(m->!m.is_supervised)

# Some commonly used ones are built-in (do not require `@load`ing):

# model type                  | does what?
# ----------------------------|----------------------------------------------
# ContinuousEncoder | transform input table to a table of `Continuous` features (see above)
# FeatureSelector | retain or dump selected features
# FillImputer | impute missing values
# OneHotEncoder | one-hot encoder `Multiclass` (and optionally `OrderedFactor`) features
# Standardizer | standardize (whiten) a vector or all `Continuous` features of a table
# UnivariateBoxCoxTransformer | apply a learned Box-Cox transformation to a vector
# UnivariateDiscretizer | discretize a `Continuous` vector, and hence render its elscitypw `OrderedFactor`

# In addition to "dynamic" transformers (ones that learn something
# from the data and must be `fit!`) users can wrap ordinary functions
# as transformers, and such *static* transformers can depend on
# parameters, like the dynamic ones. See
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/#Static-transformers-1)
# for how to define your own static transformers.

# ### Pipelines

length(schema(Xcont).names)

# Let's suppose that additionally we'd like to reduce the dimension of
# our data.  A model that will do this is `PCA` from
# `MultivariateStats`:

reducer = @load PCA

# Now, rather simply repeating the work-flow above, applying the new
# transformation to `Xcont`, we can combine both the encoding and the
# dimension-reducing models into a single model, known as a
# *pipeline*. While MLJ offers a powerful interface for composing
# models in a variety of ways, we'll stick to these simplest class of
# composite models for now. The easiest way to construct them is using
# the `@pipeline` macro:

pipe = @pipeline encoder reducer

# Notice that `pipe` is an *instance* of an automatically generated
# type (called `Pipeline<some digits>`).

# The new model behaves like any other transformer:

mach = machine(pipe, X)
fit!(mach)
Xsmall = transform(mach, X)
schema(Xsmall)

# Want to combine this pre-processing with ridge regression?

rgs = @load RidgeRegressor pkg=MLJLinearModels
pipe2 = @pipeline encoder reducer rgs

# Now our pipeline is a supervised model, instead of a transformer,
# whose performance we can evaluate:

mach = machine(pipe2, X, y)
evaluate!(mach, measure=mae, resampling=Holdout()) # CV(nfolds=6) is default

# ### Training of composite models is "smart"

# Now notice what happens if we train on all the data, then change a
# regressor hyper-parameter and retrain:

fit!(mach)

#-

pipe2.ridge_regressor.lambda = 0.1
fit!(mach)

# Second time only the ridge regressor is retrained!

# Mutate a hyper-parameter of the `PCA` model and every model except
# the `ContinuousEncoder` (which comes before it will be retrained):

pipe2.pca.pratio = 0.9999
fit!(mach)

# ### Inspecting composite models

# The dot syntax used above to change the values of *nested*
# hyper-parameters is also useful when inspecting the learned
# parameters and report generated when training a composite model:

fitted_params(mach).ridge_regressor

#-

report(mach).pca

# ### Incorporating target transformations

# Next, suppose that instead of using the raw `:price` as the
# training target, we want to use the log-price (a common practice in
# dealing with house price data). However, suppose that we still want
# to report final *predictions* on the original linear scale (and use
# these for evaluation purposes). Then we supply appropriate functions
# to key-word arguments `target` and `inverse`.

# First we'll overload `log` and `exp` for broadcasting:
Base.log(v::AbstractArray) = log.(v)
Base.exp(v::AbstractArray) = exp.(v)

# Now for the new pipeline:

pipe3 = @pipeline encoder reducer rgs target=log inverse=exp
mach = machine(pipe3, X, y)
evaluate!(mach, measure=mae)

# MLJ will also allow you to insert *learned* target
# transformations. For example, we might want to apply
# `Standardizer()` to the target, to standardize it, or
# `UnivariateBoxCoxTransformer()` to make it look Gaussian. Then
# instead of specifying a *function* for `target`, we specify a
# unsupervised *model* (or model type). One does not specify `inverse`
# because only models implementing `inverse_transform` are
# allowed.

# Let's see which of these two options results in a better outcome:

box = UnivariateBoxCoxTransformer(n=20)
stand = Standardizer()

pipe4 = @pipeline encoder reducer rgs target=box
mach = machine(pipe4, X, y)
evaluate!(mach, measure=mae)

#-

pipe4.target = stand
evaluate!(mach, measure=mae)

# ### Resources for Part 3

# - From the MLJ manual:
#     - [Transformers and other unsupervised models](https://alan-turing-institute.github.io/MLJ.jl/dev/transformers/)
#     - [Linear pipelines](https://alan-turing-institute.github.io/MLJ.jl/dev/composing_models/#Linear-pipelines-1)
# - From Data Science Tutorials:
#     - [Composing models](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/composing-models/)

# ### Exercises for Part 3

# #### Exercise 7

# Consider again the Horse Colic classification problem considered in
# Exercise 6, but with all features, `Finite` and `Infinite`:

y, X = unpack(horse, ==(:outcome), name -> true);
schema(X)

# (a) Define a pipeline that:
# - uses `Standardizer` to ensure that features that are already
#   continuous are centered at zero and have unit variance
# - re-encodes the full set of features as `Continuous`, using
#   `ContinuousEncoder`
# - uses the `KMeans` clustering model from `Clustering.jl`
#   to reduce the dimension of the feature space to `k=10`.
# - trains a `EvoTreeClassifier` (a gradient tree boosting
#   algorithm in `EvoTrees.jl`) on the reduced data, using
#   `nrounds=50` and default values for the other
#    hyper-parameters

# (b) Evaluate the pipeline on all data, using 6-fold cross-validation
# and `cross_entropy` loss.

# ☆(c) Plot a learning curve which examines the effect on this loss
# as the tree booster parameter `max_depth` varies from 2 to 10.

#

# ## Part 4 - Tuning Hyper-parameters

# ### Naive tuning of a single parameter

# The most naive way to tune a single hyper-parameter is to use
# `learning_curve`, which we already saw in Part 2. Let's see this in
# the Horse Colic classification problem, in a case where the parameter
# to be tuned is *nested* (because the model is a pipeline):

y, X = unpack(horse, ==(:outcome), name -> true);

@load LogisticClassifier pkg=MLJLinearModels
model = @pipeline Standardizer ContinuousEncoder LogisticClassifier
mach = machine(model, X, y)

#-

r = range(model, :(logistic_classifier.lambda), lower = 1e-2, upper=100, scale=:log10)

# If you're curious, you can see what `lambda` values this range will
# generate for a given resolution:

iterator(r, 5)

#-

_, _, lambdas, losses = learning_curve(mach,
                                       range=r,
                                       resampling=CV(nfolds=6),
                                       resolution=30, # default
                                       measure=cross_entropy)
plt=plot(lambdas, losses, xscale=:log10)
xlabel!(plt, "lambda")
ylabel!(plt, "cross entropy using 6-fold CV")

#-

best_lambda = lambdas[argmin(losses)]

# ### Self tuning models

# A more sophisticated way to view hyper-parameter tuning (inspired by
# MLR) is as a model *wrapper*. The wrapped model is a new model in
# its own right and when you fit it, it tunes specified
# hyper-parameters of the model being wrapped, before training on all
# supplied data. Calling `predict` on the wrapped model is like
# calling `predict` on the original model, but with the
# hyper-parameters already optimized.

# In other words, we can think of the wrapped model as a "self-tuning"
# version of the original.

# We now create a self-tuning version of the pipeline above, adding a
# parameter from the `ContinuousEncoder` to the parameters we want
# optimized.

# First, let's choose a tuning strategy (from [these
# options](https://github.com/alan-turing-institute/MLJTuning.jl#what-is-provided-here)). MLJ
# supports ordinary `Grid` search (query `?Grid` for
# details). However, as the utility of `Grid` search is limited to a
# small number of parameters, and as `Grid` searches are demonstrated
# elsewhere (see the [resources below](#resources-for-part-4)) we'll
# demonstrate `RandomSearch` here:

tuning = RandomSearch(rng=123)

# In this strategy each parameter is sampled according to a
# pre-specified prior distribution that is fit to the one-dimensional
# range object constructed using `range` as before. While one has a
# lot of control over the specification of the priors (run
# `?RandomSearch` for details) we'll let the algorithm generate these
# priors automatically.

# #### Unbounded ranges and sampling

# In MLJ a range does not have to be bounded. Furthermore, in
# `RandomSearch` a A positive unbounded range is specified will be
# sampled using a `Gamma` distribution, by default:

r = range(model, :(logistic_classifier.lambda), lower=0, origin=6, unit=5)

# Let's see what sampling using a Gamma distribution is going to mean
# for this range:

import Distributions
sampler_r = sampler(r, Distributions.Gamma)
histogram(rand(sampler_r, 10000), nbins=50)

# The second parameter we'll add to this is *nominal* (finite) and, by
# default, will be sampled uniformly. Since it is nominal, we specify
# `values` instead of `upper` and `lower` bounds:

s  = range(model, :(continuous_encoder.one_hot_ordered_factors),
           values = [true, false])
#-

# #### The tuning wrapper

# Now for the wrapper, which is an instance of `TunedModel`:

tuned_model = TunedModel(model=model,
                         ranges=[r, s],
                         resampling=CV(nfolds=6),
                         measures=cross_entropy,
                         tuning=tuning,
                         n=15)

# We can apply the `fit!/predict` work-flow to `tuned_model` just as
# for any other model:

tuned_mach = machine(tuned_model, X, y);
fit!(tuned_mach);
predict(tuned_mach, rows=1:3)

# The outcomes of the tuning can be inspected from a detailed
# report. For example, we have:

rep = report(tuned_mach);
rep.best_model

# By default, sampling of a bounded range is uniform. Lets

# In the special case of two-parameters, you can also plot the results:

plot(tuned_mach)

# Finally, let's compare cross-validation estimate of the performance
# of the self-tuning model with that of the original model (an example
# of [*nested
# resampling*](https://mlr3book.mlr-org.com/nested-resampling.html)
# here):

err = evaluate!(mach, resampling=CV(nfolds=3), measure=cross_entropy)

#-

tuned_err = evaluate!(tuned_mach, resampling=CV(nfolds=3), measure=cross_entropy)

#

# ### Resources for Part 4
#
# - From the MLJ manual:
#    - [Learning Curves](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_curves/)
#    - [Tuning Models](https://alan-turing-institute.github.io/MLJ.jl/dev/tuning_models/)
# - The [MLJTuning repo](https://github.com/alan-turing-institute/MLJTuning.jl#who-is-this-repo-for) - mostly for developers
#
# - From Data Science Tutorials:
#     - [Tuning a model](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/model-tuning/)
#     - [Crabs with XGBoost](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/crabs-xgb/) `Grid` tuning in stages for a tree-boosting model with many parameters
#     - [Boston with LightGBM](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/boston-lgbm/) -  `Grid` tuning for another popular tree-booster
#     - [Boston with Flux](https://alan-turing-institute.github.io/DataScienceTutorials.jl/end-to-end/boston-flux/) - optimizing batch size in a simple neural network regressor
# - [UCI Horse Colic Data Set](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)

# ### Exercises for Part 4

# #### Exercise 8

# This exercise continues our analysis of the King County House price
# prediction problem:

y, X = unpack(house, ==(:price), name -> true, rng=123);

# Your task will be to tune the following pipeline regression model:

@load(EvoTreeRegressor)
tree_booster = EvoTreeRegressor(nrounds = 70)
model = @pipeline ContinuousEncoder tree_booster

# (a) Construct a bounded range `r1` for the `evo_tree_booster`
# parameter `max_depth`, varying between 1 and 12.

# \star&(b) For the `nbins` parameter of the `EvoTreeRegressor`, define the range

r2 = range(model,
           :(evo_tree_regressor.nbins),
           lower = 2.5,
           upper= 7.5, scale=x->2^round(Int, x))

# Notice that in this case we've specified a *function* instead of a
# canned scale, like `:log10`. In this case the `scale` function is
# applied after sampling (uniformly) between the limits of `lower` and
# `upper`. Perhaps you can guess the outputs of the following lines of
# code?

r2_sampler = sampler(r2, Distributions.Uniform)
samples = rand(r2_sampler, 1000);
histogram(samples, nbins=50)

#-

sort(unique(samples))

# (c) Optimize `model` over these the parameter ranges `r1` and `r2`
# using a random search with uniform priors (the default). Use
# `Holdout()` resampling, and implement your search by first
# constructing a "self-tuning" wrap of `model`, as described
# above. Make `mae` (mean absolute error) the loss function that you
# optimize, and search over a total of 40 combinations of
# hyper-parameters.  If you have time, plot the results of your
# search. Feel free to use all available data.

# (d) Evaluate the best model found in the search using 3-fold
# cross-validation and compare with that of the self-tuning model
# (which is different!). Setting data hygiene concerns aside, feel
# free to use all available data.

#

# ## Part 5 - Advanced Model Composition

# > **Goals:**
# > 1. Learn how to build a prototypes of a composite model, called a *learning network*
# > 2. Learn how to use the `@from_network` macro to export a learning network as a new stand-alone model type

# While `@pipeline` is great for composing models in an unbranching
# sequence, for more complicated model composition you'll want to use
# MLJ's generic model composition syntax. There are two main steps:

# - **Prototype** the composite model by building a *learning
#   network*, which can be tested on some (dummy) data as you build
#   it.

# - **Export** the learning network as a new stand-alone model type.

# Like pipeline models, instances of the exported model type behave
# like any other model (and are not bound to any data, until you wrap
# them in a machine).

# ### Building a pipeline using the generic composition syntax

# To warm up, we'll do the equivalent of

pipe = @pipeline Standardizer LogisticClassifier;

# using the generic syntax.

# Here's some dummy data we'll be using to test our learning network:

X, y = make_blobs(5, 3)
pretty(X)

# **Step 0** - Proceed as if you were combining the models "by hand",
# using all the data available for training, transforming and
# prediction:

stand = Standardizer();
linear = LogisticClassifier();

mach1 = machine(stand, X);
fit!(mach1);
Xstand = transform(mach1, X);

mach2 = machine(linear, Xstand, y);
fit!(mach2);
yhat = predict(mach2, Xstand)

# **Step 1** - Edit your code as follows:

# - pre-wrap the data in `Source` nodes

# - delete the `fit!` calls

X = source(X)  # or X = source() if not testing
y = source(y)  # or y = source()

stand = Standardizer();
linear = LogisticClassifier();

mach1 = machine(stand, X);
Xstand = transform(mach1, X);

mach2 = machine(linear, Xstand, y);
yhat = predict(mach2, Xstand)

# Now `X`, `y`, `Xstand` and `yhat` are *nodes* ("variables" or
# "dynammic data") instead of data. All training, predicting and
# transforming is now executed lazily, whenever we `fit!` one of these
# nodes. We *call* a node to retrieve the data it represents in the
# original manual workflow.

fit!(Xstand)
Xstand() |> pretty

#-

fit!(yhat);
yhat()

# The node `yhat` is the "descendant" (in an associated DAG we have
# defined) of a unique source node:

sources(yhat)

#-

# The data at the source node is replaced by `Xnew` to obtain a
# new prediction when we call `yhat` like this:

Xnew, _ = make_blobs(2, 3);
yhat(Xnew)

# **Step 2** - Export the learning network as a new stand-alone model type

# Now, somewhat paradoxically, we can wrap the whole network in a
# special machine - called a *learning network machine* - before have
# defined the new model type. Indeed doing so is a necessary step in
# the export process, for this machine will tell the export macro:

# - what kind of model the composite will be (`Deterministic`,
#   `Probabilistic` or `Unsupervised`)a

# - which source nodes are input nodes and which are for the target

# - which nodes correspond to each operation (`predict`, `transform`,
#   etc) that we might want to define

surrogate = Probabilistic()     # a model with no fields!
mach = machine(surrogate, X, y; predict=yhat)

# Although we have no real need to use it, this machine behaves like
# you'd expect it to:

Xnew, _ = make_blobs(2, 3)
fit!(mach)
predict(mach, Xnew)

#-

# Now we create a new model type using a Julia `struct` definition
# appropriately decorated:

@from_network mach begin
    mutable struct YourPipe
        standardizer = stand
        classifier = linear::Probabilistic
    end
end

# Instantiating and evaluating on some new data:

pipe = YourPipe()
X, y = @load_iris;   # built-in data set
mach = machine(pipe, X, y)
evaluate!(mach, measure=misclassification_rate, operation=predict_mode)

# ### A composite model to average two regressor predictors

# The following is condensed version of
# [this](https://github.com/alan-turing-institute/MLJ.jl/blob/master/binder/MLJ_demo.ipynb)
# tutorial. We will define a composite model that:

# - standardizes the input data

# - learns and applies a Box-Cox transformation to the target variable

# - blends the predictions of two supervised learning models - a ridge
#  regressor and a random forest regressor; we'll blend using a simple
#  average (for a more sophisticated stacking example, see
#  [here](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/stacking/))

# - applies the *inverse* Box-Cox transformation to this blended prediction

@load RandomForestRegressor pkg=DecisionTree

# **Input layer**

X = source()
y = source()

# **First layer and target transformation**

std_model = Standardizer()
stand = machine(std_model, X)
W = MLJ.transform(stand, X)

box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, y)
z = MLJ.transform(box, y)

# **Second layer**

ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)

forest_model = RandomForestRegressor(n_trees=50)
forest = machine(forest_model, W, z)

ẑ = 0.5*predict(ridge, W) + 0.5*predict(forest, W)

# **Output**

ŷ = inverse_transform(box, ẑ)

# With the learning network defined, we're ready to export:

@from_network machine(Deterministic(), X, y, predict=ŷ) begin
    mutable struct CompositeModel
        rgs1 = ridge_model
        rgs2 = forest_model
    end
end

# Let's instantiate the new model type and try it out on some data:

composite = CompositeModel()

#-

X, y = @load_boston;
mach = machine(composite, X, y);
evaluate!(mach,
          resampling=CV(nfolds=6, shuffle=true),
          measures=[rms, mae])

# ### Resources for Part 5
#
# - From the MLJ manual:
#    - [Learning Networks](https://alan-turing-institute.github.io/MLJ.jl/stable/composing_models/#Learning-Networks-1)
# - From Data Science Tutorials:
#     - [Learning Networks](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/learning-networks/)
#     - [Learning Networks 2](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/learning-networks-2/)

#     - [Stacking](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/learning-networks-2/)
#        an advanced example of model compostion

#     - [Finer Control](https://alan-turing-institute.github.io/MLJ.jl/dev/composing_models/#Method-II:-Finer-control-(advanced)-1)
#       exporting learning networks without a macro for finer control

#

# ## Solutions to exercises

# #### Exercise 2 solution

quality = coerce(quality, OrderedFactor);
levels!(quality, ["poor", "good", "excellent"]);
elscitype(quality)

# #### Exercise 3 solution

# First pass:

coerce!(house, autotype(house));
schema(house)

#-

# All the "sqft" fields refer to "square feet" so are
# really `Continuous`. We'll regard `:yr_built` (the other `Count`
# variable above) as `Continuous` as well. So:

coerce!(house, Count => Continuous);

# And `:zipcode` should not be ordered:

coerce!(house, :zipcode => Multiclass);
schema(house)

# `:bathrooms` looks like it has a lot of levels, but on further
# inspection we see why, and `OrderedFactor` remains appropriate:

import StatsBase.countmap
countmap(house.bathrooms)

# #### Exercise 4 solution

# 4(a)

# There are *no* models that apply immediately:

models(matching(X4, y4))

# 4(b)

y4 = coerce(y4, Continuous);
models(matching(X4, y4))

# #### Exercise 6 solution

# 6(a)

y, X = unpack(horse,
              ==(:outcome),
              name -> elscitype(Tables.getcolumn(horse, name)) == Continuous);

# 6(b)(i)

model = @load LogisticClassifier pkg=MLJLinearModels;
model.lambda = 100
mach = machine(model, X, y)
fit!(mach, rows=train)
fitted_params(mach)

#-

coefs_given_feature = Dict(fitted_params(mach).coefs)
coefs_given_feature[:pulse]

#6(b)(ii)

yhat = predict(mach, rows=test); # or predict(mach, X[test,:])
err = cross_entropy(yhat, y[test]) |> mean

# 6(b)(iii)

# The predicted probabilities of the actual observations in the test
# are given by

p = broadcast(pdf, yhat, y[test]);

# The number of times this probability exceeds 50% is:
n50 = filter(x -> x > 0.5, p) |> length

# Or, as a proportion:

n50/length(test)

# 6(b)(iv)

misclassification_rate(mode.(yhat), y[test])

# 6(c)(i)

model = @load RandomForestClassifier pkg=DecisionTree
mach = machine(model, X, y)
evaluate!(mach, resampling=CV(nfolds=6), measure=cross_entropy)

r = range(model, :n_trees, lower=10, upper=70, scale=:log)

# Since random forests are inherently randomized, we generate multiple
# curves:

plt = plot()
for i in 1:4
    curve = learning_curve(mach,
                           range=r,
                           resampling=Holdout(),
                           measure=cross_entropy)
    plt=plot!(curve.parameter_values, curve.measurements)
end
xlabel!(plt, "n_trees")
ylabel!(plt, "cross entropy")
plt

# 6(c)(ii)

evaluate!(mach, resampling=CV(nfolds=9),
                measure=cross_entropy,
                rows=train).measurement[1]

model.n_trees = 90

# 6(c)(iii)

err_forest = evaluate!(mach, resampling=Holdout(),
                       measure=cross_entropy).measurement[1]

# #### Exercise 7

# (a)

@load KMeans pkg=Clustering
@load EvoTreeClassifier
pipe = @pipeline(Standardizer,
                 ContinuousEncoder,
                 KMeans(k=10),
                 EvoTreeClassifier(nrounds=50))

# (b)

mach = machine(pipe, X, y)
evaluate!(mach, resampling=CV(nfolds=6), measure=cross_entropy)

# (c)

r = range(pipe, :(evo_tree_classifier.max_depth), lower=1, upper=10)

curve = learning_curve(mach,
                       range=r,
                       resampling=CV(nfolds=6),
                       measure=cross_entropy)

plt = plot(curve.parameter_values, curve.measurements)
xlabel!(plt, "max_depth")
ylabel!(plt, "CV estimate of cross entropy")
plt

# Here's a second curve using a different random seed for the booster:

pipe.evo_tree_classifier.seed = 123
curve = learning_curve(mach,
                       range=r,
                       resampling=CV(nfolds=6),
                       measure=cross_entropy)
plot!(curve.parameter_values, curve.measurements)

# #### Exercise 8

y, X = unpack(house, ==(:price), name -> true, rng=123);

@load(EvoTreeRegressor)
tree_booster = EvoTreeRegressor(nrounds = 70)
model = @pipeline ContinuousEncoder tree_booster

# (a)

r1 = range(model, :(evo_tree_regressor.max_depth), lower=1, upper=12)

# (c)

tuned_model = TunedModel(model=model,
                         ranges=[r1, r2],
                         resampling=Holdout(),
                         measures=mae,
                         tuning=RandomSearch(rng=123),
                         n=40)

tuned_mach = machine(tuned_model, X, y) |> fit!
plot(tuned_mach)

# (d)

best_model = report(tuned_mach).best_model;
best_mach = machine(best_model, X, y);
best_err = evaluate!(best_mach, resampling=CV(nfolds=3), measure=mae)

#-

tuned_err = evaluate!(tuned_mach, resampling=CV(nfolds=3), measure=mae)

using Literate #src
Literate.markdown(@__FILE__, DIR, execute=true) #src
Literate.notebook(@__FILE__, DIR, execute=false) #src


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#
using DiscreteEvents, .Threads

println("... check if threading is available ...")
clk = Clock()

if VERSION ≥ v"1.3"
    if nthreads() > 1
        include("test_threads1.jl")
        include("test_threads2.jl")
    else
        @test_warn "no parallel threads available!" fork!(clk)
    end
else
    @test_warn "threading is available for ≥ 1.3!" fork!(clk)
end


text:
using DocStringExtensions
using ProgressMeter
using Gadfly
using Cairo
using RoMEPlotting
Gadfly.set_default_plot_size(35cm,25cm)

function layerBeamPatternRose(bear::BallTreeDensity;
                              scale::Float64=1.0,
                              c=colorant"magenta",
                              wRr=TU.R(0.0),
                              wTRr=zeros(2) )
  #
  tp = reshape([0:0.01:(2pi);], 1, :)
  belp = scale*bear(tp)
  Gadfly.plot(x=tp, y=belp, Geom.path)
  belRose = zeros(2, length(tp))
  belRose[1,:] = belp

  idx = 0
  for rRc in TU.R.(tp)
    idx += 1
    belRose[:,idx] = wRr*rRc*belRose[:,idx] + wTRr
  end

  Gadfly.layer(x=belRose[1,:], y=belRose[2,:], Geom.path, Theme(default_color=c))
end

function plotSandsharkFromDFG(dfg::AbstractDFG;
                              radix::Float64=1.8,
                              scale::Float64=0.4,
                              spscale::Float64=1.0,
                              drawTriads::Bool=true,
                              lbls::Bool=true  )
  # bearing range factors
  brf = ls(dfg, Pose2Point2BearingRange)
  #variables and poses connected to br factors
  brvars = 0 < length(brf) ? union(map(x->ls(dfg, x), brf)...) : Symbol[]
  brposes = setdiff(brvars, ls(dfg, r"l"))
  brposes .= sortDFG(brposes)

  # get body frame beam estimate
  brfacts = Dict{Symbol, BallTreeDensity}()

  PL =  drawPosesLandms(dfg, spscale=spscale, contour=false, meanmax=:mean, drawTriads=drawTriads, lbls=lbls)

  # TODO: not matched properly
  for i in 1:length(brposes)
    @show brfi = intersect(ls(dfg, brposes[i]), brf)[1]
    fct = getFactorType(getFactor(dfg, brfi))
    brfacts[brposes[i]] = fct.bearing
    rVo = getKDEMax(getKDE(dfg, brposes[i]))
    pll = AMP.plotKDECircular([fct.bearing;], rVo=rVo, radix=radix, scale=scale)
    for lyr in pll.layers
      push!(PL.layers, lyr)
    end
  end

  PL
end

function drawPosesLandmarksAndOdo(fg::G,
                                  ppbrDict::Dict,
                                  navkeys::Array,
                                  X::Array,
                                  Y::Array,
                                  lblkeys::Array,
                                  lblX::Array,
                                  lblY::Array,
                                  gTitle::String="",
                                  drawhist=true  ) where G <: AbstractDFG
    #
    PLL = []
    # xx = sortVarNested(ls(fg, r"x"))
    # idx = 0
    # for idx in 1:length(xx)
    #   ep = epochs[idx]
    #   theta = getKDEMax(getBelief(fg, xx[idx]))
    #   wRr = TU.R(theta[3])
    #   pll = layerBeamPatternRose(ppbrDict[ep].bearing, wTRr=theta[1:2], scale=5.0 , wRr=wRr)
    #   push!(PLL, pll)
    # end
    #
    # pllandmarks = drawPosesLandms(fg, spscale=2.5, drawhist=drawhist)
    pllandmarks = plotSandsharkFromDFG(fg, scale=1.8)

    # Add odo
    navdf = DataFrame(
      ts = navkeys,
      x = X,
      y = Y
    )
    # pl = Gadfly.layer(navdf, x=:x, y=:y, Geom.path())
    push!(pllandmarks.layers, Gadfly.layer(navdf, x=:x, y=:y, Geom.path(), Theme(default_color=colorant"gray80"))[1])
    lbldf = DataFrame(
      ts = lblkeys .- lblkeys[1],
      x = lblX,
      y = lblY
    )
    push!(pllandmarks.layers, Gadfly.layer(lbldf, x=:x, y=:y, Geom.path(), Theme(default_color=colorant"magenta"))[1])
    # Gadfly.plot(pl.layers)

    # Make X/Y range same so no distorted
    # push!(PLL, Coord.Cartesian(xmin=-160.0,xmax=10.0,ymin=-135.0,ymax=35.0))
    pla = Gadfly.plot([PLL...;pllandmarks.layers; ]...,
        Guide.manual_color_key("Legend", ["Dead reckoning", "LBL Path", "Non-Gaussian SLAM"], ["gray80", "green", "light blue"]),
        Guide.title(gTitle))
    return pla
end

function plotRawBeamFormer(fg1::AbstractDFG, ppbrDict, NAV, interp_x, interp_y, interp_yaw; show::Bool=true, all::Bool=true)

  meta = readdlm(joinpath(getSolverParams(fg1).logpath, "metadata.txt"),',')
  mkpath(joinpath(getSolverParams(fg1).logpath, "rawplots"))

  files = String[]
  @showprogress "draw raw" for i in 1:(all ? size(meta,1) : 1)
    ep = Int(meta[i,2]) # 1531152901000000000
    bfspls = readdlm("/home/dehann/data/sandshark/full_wombat_2018_07_09/extracted/beamformer/particles/$(ep).csv",',')

    azipts = bfspls[:,1]
    aziptsw = TU.wrapRad.(azipts)
    aziprobl = kde!(azipts)
    npts = rand(aziprobl, 200)
    aziprob = manikde!(npts, Sphere1)

    plraw = Gadfly.plot(x=azipts, Geom.histogram(density=true), Guide.title("$ep, $(meta[i,3])"))
    pc = plotKDECircular(aziprob, title="recalc from raw (b-frame)")
    pbr= plotKDECircular(ppbrDict[ep].bearing, title="Pose2Point2BearingRange for fg (body)")
    pll = Gadfly.plot(x=[17;interp_x(ep)],y=[1.8;interp_y(ep)], Geom.line, Coord.Cartesian(xmin=-2.0+interp_x(ep),xmax=2.0+interp_x(ep),ymin=-  2.0+interp_y(ep),ymax=2.0+interp_y(ep)), Guide.title("n-frame, NAV yaw=$(round(interp_yaw(ep), digits=3))"))
    pbR= plotKDECircular(ppbrDict[ep].bearing, rVo=[interp_x(ep);interp_y(ep);interp_yaw(ep)], title="n-frame, NAV yaw=$(round(interp_yaw(ep), digits=3))")

    union!(pll.layers, pbR.layers)

    RoMEPlotting.addXYLineLayers!(pll, [interp_x(ep);], [interp_y(ep);], [interp_yaw(ep);], l=0.5)

    file = joinpath(getSolverParams(fg1).logpath, "rawplots", "bf_$(meta[i,3])_$ep.pdf")
    hstack(plraw, vstack(pc, pbr), pll) |> PDF(file)
    push!(files, file)
  end
  # compile into one PDF
  outfile = joinpath(getSolverParams(fg1).logpath, "bf_out.pdf")
  push!(files, outfile)
  run(`pdfunite $files`)

  !show ? nothing : (@async run(`evince $outfile`))
  return outfile
end

"""
    $SIGNATURES

Plot the timing for Sandshark AUV frontend timing analysis.

"""
function plotFrontendTiming(wtdsh)
  @show wtdsh[1]
  # line = [  now()
            # getLastPoses(fg,filterLabel=r"x\d",number=1)[1]
            # dashboard[:poseStride]
            # dashboard[:canTakePoses]
            # dashboard[:solveInProgress]
            # dashboard[:realTimeSlack].value
            # length(dashboard[:poseSolveToken].data) ]
  #
  ttt = (x->datetime2unix(DateTime(x[1]))).(wtdsh)
  ttt .-= ttt[1]
  wtv = (x->string(x[2])).(wtdsh)
  wtr = (x->parse(Int,string(x[3]))).(wtdsh)
  wth = (x->getfield(Main,Symbol(string(x[4])))).(wtdsh)
  wts = (x->getfield(Main,Symbol(string(x[5])))).(wtdsh)
  wtt = (x->parse(Float64,string(x[6]))).(wtdsh)
  wttn = wtt./wtt[end]
  wtl = (x->parse(Int,string(x[7]))).(wtdsh)

  lpn = wtv .|> x->parse(Int, string(x)[2:end])

  # Gadfly.set_default_plot_size(45cm,25cm)

  pl = Gadfly.plot(
    Gadfly.layer(x=ttt,y=(wts .|> Int)*0.5.+1, Geom.line, Theme(default_color=colorant"red")),
    Gadfly.layer(x=ttt,y=(wtl .|> Int)*0.5, Geom.line, Theme(default_color=colorant"khaki4")),
    Gadfly.layer(x=ttt,y=(wth .|> Int)*0.333.-1, Geom.line),
    Gadfly.layer(x=ttt,y=(lpn.%10)*0.1.-2, Geom.line, Theme(default_color=colorant"magenta")),
    # Gadfly.layer(x=ttt,y=(wtr.%10)*0.1.-2, Geom.line, Theme(default_color=colorant"green")),
    Gadfly.layer(x=ttt,y=wttn.-3, Geom.line, Theme(default_color=colorant"green")),
  )
end
# Gadfly.plot(y=wtt, Geom.line)
# Gadfly.plot(y=wttn, Geom.line)

#


text:
encode_init(vi::VorbisInfo,channels,rate,max_bitrate,nominal,min_bitrate) =
    ccall((:vorbis_encode_init,vorbisenc),Cint,
          (Ref{VorbisInfo},Clong,Clong,Clong,Clong,Clong),
          vi,channels,rate,max_bitrate,nominal,min_bitrate)

encode_setup_managed(vi::VorbisInfo,channels,rate,max_bitrate,nominal,min_bitrate) =
    ccall((:vorbis_encode_setup_managed,vorbisenc),Cint,
          (Ref{VorbisInfo},Clong,Clong,Clong,Clong,Clong),
          vi,channels,rate,max_bitrate,nominal,min_bitrate)

encode_setup_vbr(vi::VorbisInfo,channels,rate,quality) =
    ccall((:vorbis_encode_setup_vbr,vorbisenc),Cint,
          (Ref{VorbisInfo},Clong,Clong,Cfloat),vi,channels,rate,quality)

encode_init_vbr(vi::VorbisInfo,channels,rate,quality) =
    ccall((:vorbis_encode_init_vbr,vorbisenc),Cint,
          (Ref{VorbisInfo},Clong,Clong,Cfloat),vi,channels,rate,quality)

encode_setup_init(vi::VorbisInfo) =
    ccall((:vorbis_encode_setup_init,vorbisenc),Cint,(Ref{VorbisInfo},),vi)

immutable RateManage
    mangement_active::Cint
    bitrate_hard_min::Clong
    bitrate_hard_max::Clong
    bitrate_hard_window::Cdouble
    bitrate_av_lo::Clong
    bitrate_av_hi::Clong
    bitrate_av_window::Cdouble
    bitrate_av_window_center::Cdouble
end

immutable RateManage2
    management_active::Cint
    bitrate_limit_min_kbps::Clong
    bitrate_limit_max_kbps::Clong
    bitrate_limit_reservoir_bits::Clong
    bitrate_limit_reservoir_bias::Cdouble
    bitrate_average_kbps::Clong
    bitrate_average_damping::Cdouble
end


text:
## Imports

using BenchmarkTools
using Graphs
using MultiAgentPathFinding
using PythonCall
using ProgressMeter
using UnicodePlots

## Test

rail_generators = pyimport("flatland.envs.rail_generators")
line_generators = pyimport("flatland.envs.line_generators")
rail_env = pyimport("flatland.envs.rail_env")

rail_generator = rail_generators.sparse_rail_generator(; max_num_cities=6)
line_generator = line_generators.sparse_line_generator()

pyenv = rail_env.RailEnv(;
    width=40,
    height=40,
    number_of_agents=50,
    rail_generator=rail_generator,
    line_generator=line_generator,
    random_seed=11,
)

pyenv.reset();
mapf = flatland_mapf(pyenv);

## Local search

solution_indep = independent_astar(mapf);
is_feasible(solution_indep, mapf)
flowtime(solution_indep, mapf)

solution_indep2 = independent_dijkstra(mapf);
is_feasible(solution_indep2, mapf)
flowtime(solution_indep2, mapf)

solution_lns2 = feasibility_search(
    mapf; neighborhood_size=5, conflict_price=1, conflict_price_increase=1e-2
);
is_feasible(solution_lns2, mapf)
flowtime(solution_lns2, mapf)

solution_coop = cooperative_astar(mapf);
is_feasible(solution_coop, mapf)
flowtime(solution_coop, mapf)

solution_coop_soft = cooperative_astar(mapf; conflict_price=100);
is_feasible(solution_coop_soft, mapf)
flowtime(solution_coop_soft, mapf)

solution_lns = large_neighborhood_search(mapf; neighborhood_size=5, steps=1000);
is_feasible(solution_lns, mapf)
flowtime(solution_lns, mapf)

## (I)LP

# tmax = maximum(t for path in solution_lns for (t, v) in path)

# _, _, solution_lp = solve_lp(mapf, T=tmax+10, integer=false, capacity=true);
# is_feasible(solution_lp, mapf)
# flowtime(solution_lp, mapf)

# _, _, solution_lp_indep = solve_lp(mapf, T=tmax+10, integer=true, capacity=false);
# is_feasible(solution_lp_indep, mapf)
# flowtime(solution_lp_indep, mapf)

# _, _, solution_ilp = solve_lp(mapf, T=tmax+10, integer=true, capacity=true);
# is_feasible(solution_ilp, mapf)
# flowtime(solution_ilp, mapf)

## Animation

# using GLMakie

# fig, (A, XY, M) = plot_flatland_graph(mapf);
# fig
# solution = copy(solution_coop);
# framerate = 5
# tmax = maximum(t for path in solution for (t, v) in path)
# @showprogress for t in 1:tmax
#     A[], XY[], M[] = flatland_agent_coords(mapf, solution, t)
#     sleep(1 / framerate)
# end


text:
using Klara

# Precision matrix
C = Hyperparameter(:C)

p = BasicContMuvParameter(
  :p,
  logtarget=(p::Vector{Float64}, v::Vector) -> -dot(p, v[1]*p),
  gradlogtarget=(p::Vector{Float64}, v::Vector) -> -2*v[1]*p,
  tensorlogtarget=(p::Vector{Float64}, v::Vector) -> -2*v[1],
  nkeys=2
)

model = GenericModel([C, p], isindexed=false)

mcsampler = SMMALA(1.25, H -> softabs(H, 1000.))

tuner = VanillaMCTuner()

mcrange = BasicMCRange(nsteps=10000, burnin=1000)

v0 = Dict(:C=>inv([1. 0.8; 0.8 1.]), :p=>Float64[1.25, 3.11])

outopts = Dict(:monitor=>[:value, :logtarget, :gradlogtarget], :diagnostics=>[:accept])

job = BasicMCJob(model, mcsampler, mcrange, v0, outopts=outopts)

@time run(job)

chain = output(job)

mean(chain)

acceptance(chain)


text:
struct X
    X
end


text:
using LinearAlgebra
using Plots
using CSV
using KernelFunctions
using Distributions
using DataFrames

include("../algo/kernels.jl")
include("../algo/regularized_Picard.jl")
include("../algo/utils.jl")

"""
    convergence_exact_solution()

Visualize the convergence of regularized Picard to exact solution (see Proposition 10)
"""
function convergence_exact_solution()

    # width
    sigma = .01
    # regularizer
    lambda =  .05
    # regularizer for positive definiteness
    epsilon = 1e-10;
    # relative objective tolerance
    tol = 1e-6

    ############# load 1 DPP samples #############
    # intensity of the generating DPP
    intensity = 100; # 50 or 100
    strg_intensity = string(intensity)

    # create an array of arrays
    indices_DPP_samples = Array{Int64,1}[];

    print("Loading DPP samples from: \n")
    i = 1; # id of the DPP sample
    file_name = "data/statspats/samples/GaussDPPsample_alpha0_00p5_rho0_"*strg_intensity*"_nb_"*string(i)*".csv"
    print(file_name, " \n")
    temp = CSV.File(file_name; header=true) |> Tables.matrix
    temp = temp[:,2:3];
    id_temp = collect(1:size(temp,1));
    push!(indices_DPP_samples,id_temp);
    total_samples = temp;

    # for approximating Fredholm determinant
    indices_Fredholm_sample = indices_DPP_samples[1]

    ############# exact solution #############

    # create full kernel matrix
    x = (total_samples)'/sigma;
    k = SqExponentialKernel();
    K = kernelmatrix(k, x) + epsilon *I ;

    # Chol decomposition
    R = cholesky(K).U;
    Rinv = inv(R);

    # Exact solution
    m = size(K,1);
    X_exact = 0.5*(real(sqrt(m^2*I + 4*m*K/lambda))-m*I);
    B_exact = Rinv'*X_exact*Rinv;

    # objective exact solution
    obj_det_exact,obj_reg_exact = Picard_objective(B_exact, indices_DPP_samples, indices_Fredholm_sample, R,lambda)
    obj_exact = obj_det_exact + obj_reg_exact;

    ############# approximate solution #############

    # initial positive definite iterate
    X = randn(size(K));
    B = X*X'+ UniformScaling(epsilon);

    use_inverse = false

    n_steps = 50;
    it_max = 1;
    err = ones(n_steps,1);
    obj = ones(n_steps,1);

    for j=1:n_steps
        B, R, obj_j, i_stop = regularized_Picard(B, K, indices_DPP_samples, indices_Fredholm_sample, lambda, it_max ,tol,use_inverse)
        err[j] = norm(B-B_exact)/norm(B_exact);
        obj[j] = obj_j[end];
    end

    #############  plotting #############

    plt_error = plot(err, yaxis=:log,legend=false,ylabel = "relative error",xlabel = "iteration")
    display(plt_error)

    plt_obj = plot(obj,legend=false,ylabel = "objective",xlabel = "iteration")
    plot!(obj_exact*ones(size(obj)))
    display(plt_obj)
end


text:
"""
A [`User`](@ref) connection to an external service (Twitch, YouTube, etc.).
More details [here](https://discordapp.com/developers/docs/resources/user#connection-object).
"""
struct Connection
    id::String
    name::String
    type::String
    revoked::Bool
    integrations::Vector{Integration}
end
@boilerplate Connection :constructors :docs :lower :merge :mock


text:
"""
    ESS(chain, k)

Calculate the Effective Sample Size.

# Arguments
- `chain`: matrix of draws from the target distribution
- `k`: (integer) number of lags

# Output
- `ESS`: vector of effective sample sizes
"""

function ESS(chain::Matrix{T} where T<:AbstractFloat, k = 10::Integer)
    M, J = size(chain)
    @assert M > k "The chain of β must be longer than k."
    ρ = sum(pacf(chain, 1:k); dims = 1)
    τ = 1 .+ 2 .* sum(ρ; dims = 1)
    ESS = round.(Int, vec(M ./ τ))
    return ESS
end

"""
    probability(x, chain)

Calculate probabilities based on the probit link.ß

# Arguments
- `x`: matrix of inputs
- `chain`: matrix of draws from the target distribution

# Output
- `p`: matrix of probabilities
"""

function probability(x::Matrix{T} where T<:Real, chain::Matrix{T} where T<:AbstractFloat)
    z = x * chain'
    p = cdf(Normal(0, 1), z)
    return p
end
probability(x::Vector{T} where T<:Real, chain::Matrix{T} where T<:AbstractFloat) = probability(repeat(x, 1, 1), chain)


text:
# This is a simple test file that checks (on generated data) whether the
# Kalman.jl functions for AR(1), CARMA(p,q), and multi-epoch CARMA(p,q) can
# accurately reconstruct the parameters from a pre-generated data set.

using DelimitedFiles
using Ensemble
using HDF5
using Statistics

@testset "Kalman.jl Test Suite" begin
    @testset "AR(1) Fit Test" begin
        ar1_parameters = Dict("sigma" => Float64(pi), "mu" => 123.32, "tau" => sqrt(2), "nu" => 0.77)
        data = readdlm("ar1_data.txt")
        ts = data[:,1]
        ys = data[:,2]
        dys = data[:,3]

        post = Kalman.AR1KalmanPosterior(ts, ys, dys)

        pts = Kalman.init(post, 128)
        ns = EnsembleNest.NestState(x->Kalman.log_likelihood(post,x), x->Kalman.log_prior(post,x), pts, 128)
        EnsembleNest.run!(ns, 0.1, verbose=true)
        ps, lls = EnsembleNest.postsample(ns)

        mus = [Kalman.to_params(post, ps[:,j]).mu for j in 1:size(ps,2)]
        sigmas = [Kalman.to_params(post, ps[:,j]).sigma for j in 1:size(ps,2)]
        taus = [Kalman.to_params(post, ps[:,j]).tau for j in 1:size(ps,2)]
        nus = [Kalman.to_params(post, ps[:,j]).nu for j in 1:size(ps,2)]

        @test abs(mean(mus) - ar1_parameters["mu"]) < 2*std(mus)
        @test abs(mean(sigmas) - ar1_parameters["sigma"]) < 2*std(sigmas)
        @test abs(mean(taus) - ar1_parameters["tau"]) < 2*std(taus)
        @test abs(mean(nus) - ar1_parameters["nu"]) < 2*std(nus)
    end

    @testset "CARMA Fit Test" begin
        mu = 123.32
        sigma = Float64(pi)
        nu = 0.77
        f0 = 0.866
        r0 = 1.0/1.33
        b1 = 1.0/0.235
        b2 = 1.0/1.57
        Q = 10.0
        arroots = [-r0, -f0/Q+f0*1im, -f0/Q-f0*1im]
        maroots = [-b1, -b2]

        data = readdlm("carma_data.txt")
        ts = data[:,1]
        ys = data[:,2]
        dys = data[:,3]

        post = Kalman.CARMAKalmanPosterior(ts, ys, dys, 3, 2)
        pts = Kalman.init(post, 128)
        ns = EnsembleNest.NestState(x->Kalman.log_likelihood(post,x), x->Kalman.log_prior(post,x), pts, 128)
        EnsembleNest.run!(ns, 0.1, verbose=true)
        ps, lls = EnsembleNest.postsample(ns)

        mus = [Kalman.to_params(post, ps[:,j]).mu for j in 1:size(ps,2)]
        sigmas = [Kalman.to_params(post, ps[:,j]).sigma for j in 1:size(ps,2)]
        nus = [Kalman.to_params(post, ps[:,j]).nu for j in 1:size(ps,2)]
        fs = vec(maximum(Kalman.frequencies(post, ps), dims=1))

        @test abs(mean(mus) - mu) < 2*std(mus)
        @test abs(mean(sigmas) - sigma) < 2*std(sigmas)
        @test abs(mean(nus) - nu) < 2*std(nus)
        @test abs(mean(fs) - f0) < 2*std(fs)
    end

    @testset "Multi-Segment CARMA Fit" begin
        mu = 123.32
        sigma = Float64(pi)
        nu = 0.77
        f0 = 0.866
        r0 = 1.0/1.33
        b1 = 1.0/0.235
        b2 = 1.0/1.57
        Q = 10.0
        arroots = [-r0, -f0/Q+f0*1im, -f0/Q-f0*1im]
        maroots = [-b1, -b2]

        data = readdlm("carma_data.txt")
        ts = data[:,1]
        ys = data[:,2]
        dys = data[:,3]

        dmu2 = -25.3234
        dmu3 = +136.0

        nufac2 = 1.33
        nufac3 = 0.756

        ts = [ts[1:3:end], ts[2:3:end], ts[3:3:end]]
        ys = [ys[1:3:end], ys[2:3:end].+dmu2, ys[3:3:end].+dmu3]
        dys = [dys[1:3:end], dys[2:3:end].*nufac2, dys[3:3:end].*nufac3]

        post = Kalman.MultiSegmentCARMAKalmanPosterior(ts, ys, dys, 3, 2)
        pts = Kalman.init(post, 128)
        ns = EnsembleNest.NestState(x->Kalman.log_likelihood(post,x), x->Kalman.log_prior(post,x), pts, 128)
        EnsembleNest.run!(ns, 0.1, verbose=true)
        ps, lls = EnsembleNest.postsample(ns)

        mus = hcat([Kalman.to_params(post, ps[:,j]).mu for j in 1:size(ps,2)]...)
        nus = hcat([Kalman.to_params(post, ps[:,j]).nu for j in 1:size(ps,2)]...)
        sigmas = [Kalman.to_params(post, ps[:,j]).sigma for j in 1:size(ps, 2)]
        fs = vec(maximum(Kalman.frequencies(post, ps), dims=1))

        mu0 = [mu, mu+dmu2, mu+dmu3]
        nu0 = [nu, nu/nufac2, nu/nufac3]

        for i in 1:3
            @test abs(mean(mus[i,:]) - mu0[i]) < 2*std(mus[i,:])
            @test abs(mean(nus[i,:]) - nu0[i]) < 2*std(nus[i,:])
        end
        @test abs(mean(sigmas) - sigma) < 2*std(sigmas)
        @test abs(mean(fs) - f0) < 2*std(fs)
    end

    @testset "Read/write to HDF5" begin
        N = 128
        p = 3
        q = 2
        ts = sort(randn(N))
        ys = randn(N)
        dys = abs.(randn(N))

        post = Kalman.CARMAKalmanPosterior(ts, ys, dys, p, q)

        try
            h5open("test-CARMA-Kepler.h5", "w") do f
                write(f, post)
            end

            h5open("test-CARMA-Kepler.h5", "r") do f
                global post2 = Kalman.CARMAKalmanPosterior(f)
            end

            @test all(post.ts .== post2.ts)
            @test all(post.ys .== post2.ys)
            @test all(post.dys .== post2.dys)
            @test post.p == post2.p
            @test post.q == post2.q
        finally
            rm("test-CARMA-Kepler.h5")
        end
    end
end


text:
using InstantiateFromURL
activate_github("QuantEcon/QuantEconLecturePackages", tag = "v0.9.5");

max_stars_repo_name: aderdzyan/QuantEcon-lecture-source-jl

max_issues_repo_name: aderdzyan/QuantEcon-lecture-source-jl

max_forks_repo_name: aderdzyan/QuantEcon-lecture-source-jl


text:
#
# This file is part of ThreadingExperiments.jl
#
# Paul Bayer, 2020
# license: MIT
#
"""
    onthread(f::F, id::Int) where {F<:Function}

Execute a function f on thread id.

# Examples, usage

```julia
julia> using ThreadingExperiments, .Threads

julia> onthread(threadid, 2)
2

julia> onthread(3) do; threadid(); end
3

julia> onthread(4) do
           threadid()
       end
4
```
"""
function onthread(f::F, id::Int) where {F<:Function}
    t = Task(nothing)
    @assert id in 1:nthreads() "thread $id not available!"
    @threads for i in 1:nthreads()
        if i == id
            t = @async f()
        end
    end
    fetch(t)
end


text:
module TypeInstantiationTests

using Interpolations, Test

# NO DIMSPECS
# tests that we forward types correctly to the instance constructors
const A1 = rand(15)
const A2 = rand(15, 10)

# B-splines
for GT in (OnCell, OnGrid)
    # ...without boundary conditions
    for D in (Constant, Linear)
        @inferred(interpolate(A1, BSpline(D), GT))
        @inferred(interpolate(A1, BSpline(D), GT()))
        @inferred(interpolate(A1, BSpline(D()), GT))
        @inferred(interpolate(A2, BSpline(D), GT))
        @inferred(interpolate(A2, BSpline(D), GT()))
        @inferred(interpolate(A2, BSpline(D()), GT))

        @inferred(interpolate!(copy(A1), BSpline(D), GT))
        @inferred(interpolate!(copy(A1), BSpline(D), GT()))
        @inferred(interpolate!(copy(A1), BSpline(D()), GT))
        @inferred(interpolate!(copy(A2), BSpline(D), GT))
        @inferred(interpolate!(copy(A2), BSpline(D), GT()))
        @inferred(interpolate!(copy(A2), BSpline(D()), GT))
    end

    # Quadratic
    for BC in (Flat, Line, Periodic, Reflect, Free)
        if BC != InPlace
            @inferred(interpolate(A1, BSpline(Quadratic(BC)), GT))
            @inferred(interpolate(A1, BSpline(Quadratic(BC())), GT))
            @inferred(interpolate(A1, BSpline(Quadratic(BC)), GT()))
            @inferred(interpolate(A2, BSpline(Quadratic(BC)), GT))
            @inferred(interpolate(A2, BSpline(Quadratic(BC())), GT))
            @inferred(interpolate(A2, BSpline(Quadratic(BC)), GT()))
        end

        @inferred(interpolate!(copy(A1), BSpline(Quadratic(BC)), GT))
        @inferred(interpolate!(copy(A1), BSpline(Quadratic(BC())), GT))
        @inferred(interpolate!(copy(A1), BSpline(Quadratic(BC)), GT()))
        @inferred(interpolate!(copy(A2), BSpline(Quadratic(BC)), GT))
        @inferred(interpolate!(copy(A2), BSpline(Quadratic(BC())), GT))
        @inferred(interpolate!(copy(A2), BSpline(Quadratic(BC)), GT()))
    end
end

# Gridded
const knots1 = (sort(rand(15)),)
const knots2 = (sort(rand(15)), sort(rand(10)))
for D in (Constant, Linear)
    @inferred(interpolate(knots1, A1, Gridded(D)))
    @inferred(interpolate(knots2, A2, Gridded(D)))

    @inferred(interpolate!(knots1, copy(A1), Gridded(D)))
    @inferred(interpolate!(knots2, copy(A2), Gridded(D)))
end

# DIMSPECS

# test that constructing dimspecs work
for T in (
    Tuple{OnGrid,OnCell},
    Tuple{BSpline{Linear},BSpline{Quadratic{Free}}}
    )
    @test isa(@inferred(Interpolations.construct_instance(T)), T)
end
# sample one dimspec for each interpolation constructor to see that it
# calls the construct_instance() function correctly
@inferred(interpolate(A2, Tuple{BSpline{Linear},BSpline{Constant}}, Tuple{OnGrid,OnCell}))
@inferred(interpolate!(copy(A2), Tuple{BSpline{Linear},BSpline{Constant}}, Tuple{OnGrid,OnCell}))
@inferred(interpolate(A2, Tuple{BSpline{Linear},BSpline{Constant}}, (OnGrid(),OnCell())))
@inferred(interpolate!(copy(A2), Tuple{BSpline{Linear},BSpline{Constant}}, (OnGrid(),OnCell())))
@inferred(interpolate(A2, (BSpline(Linear()),BSpline(Constant())), Tuple{OnGrid,OnCell}))
@inferred(interpolate!(copy(A2), (BSpline(Linear()),BSpline(Constant())), Tuple{OnGrid,OnCell}))

end


text:
# This file has been generated automatically by 'generate.jl'

__precompile__()

module Traits

export VertexTrait, EdgeTrait, refines

@enum VertexTrait begin
  CustomScheduler
  InstrumentedCommunicationInterconnect
  Goal
  TimeDivisionMultiplexer
  AbstractOrdering
  FixedPriorityScheduler
  BufferSignal
  ExtraFunctional
  SporadicTask
  TimeTriggeredScheduler
  SYComb
  ReactorTimer
  WCET
  LabelSignal
  HardRequirement
  Signal
  RoundRobinScheduler
  InstrumentedSignal
  AbstractPhysicalComponent
  AbstractProcessingComponent
  ForSyDeFunction
  ReactorElement
  WCCT
  AbstractGrouping
  TriggeredTask
  Requirement
  InstrumentedFunction
  SDFComb
  AbsractInterfaceComponent
  ReactorActor
  SYPrefix
  InstrumentedProcessorTile
  LocationRequirement
  MinimumThroughput
  SDFPrefix
  AbstractStorageComponent
  Instrumented
  AbstractCommunicationComponent
end

function refines(t1::VertexTrait, t2::VertexTrait)
  if t1 == CustomScheduler
    if t2 == CustomScheduler
      return true
    end
  end
  if t1 == InstrumentedCommunicationInterconnect
    if t2 == InstrumentedCommunicationInterconnect
      return true
    end
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbstractCommunicationComponent
      return true
    end
  end
  if t1 == Goal
    if t2 == Goal
      return true
    end
  end
  if t1 == TimeDivisionMultiplexer
    if t2 == TimeDivisionMultiplexer
      return true
    end
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbstractCommunicationComponent
      return true
    end
  end
  if t1 == AbstractOrdering
    if t2 == AbstractOrdering
      return true
    end
    if t2 == AbstractGrouping
      return true
    end
  end
  if t1 == FixedPriorityScheduler
    if t2 == FixedPriorityScheduler
      return true
    end
  end
  if t1 == BufferSignal
    if t2 == BufferSignal
      return true
    end
    if t2 == Signal
      return true
    end
  end
  if t1 == ExtraFunctional
    if t2 == ExtraFunctional
      return true
    end
  end
  if t1 == SporadicTask
    if t2 == AbstractOrdering
      return true
    end
    if t2 == SporadicTask
      return true
    end
    if t2 == AbstractGrouping
      return true
    end
  end
  if t1 == TimeTriggeredScheduler
    if t2 == TimeTriggeredScheduler
      return true
    end
    if t2 == AbstractGrouping
      return true
    end
  end
  if t1 == SYComb
    if t2 == SYComb
      return true
    end
    if t2 == ForSyDeFunction
      return true
    end
  end
  if t1 == ReactorTimer
    if t2 == ReactorTimer
      return true
    end
    if t2 == ReactorElement
      return true
    end
  end
  if t1 == WCET
    if t2 == ExtraFunctional
      return true
    end
    if t2 == WCET
      return true
    end
  end
  if t1 == LabelSignal
    if t2 == LabelSignal
      return true
    end
    if t2 == Signal
      return true
    end
  end
  if t1 == HardRequirement
    if t2 == HardRequirement
      return true
    end
    if t2 == Requirement
      return true
    end
  end
  if t1 == Signal
    if t2 == Signal
      return true
    end
  end
  if t1 == RoundRobinScheduler
    if t2 == RoundRobinScheduler
      return true
    end
  end
  if t1 == InstrumentedSignal
    if t2 == Signal
      return true
    end
    if t2 == InstrumentedSignal
      return true
    end
    if t2 == Instrumented
      return true
    end
  end
  if t1 == AbstractPhysicalComponent
    if t2 == AbstractPhysicalComponent
      return true
    end
  end
  if t1 == AbstractProcessingComponent
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbstractProcessingComponent
      return true
    end
  end
  if t1 == ForSyDeFunction
    if t2 == ForSyDeFunction
      return true
    end
  end
  if t1 == ReactorElement
    if t2 == ReactorElement
      return true
    end
  end
  if t1 == WCCT
    if t2 == ExtraFunctional
      return true
    end
    if t2 == WCCT
      return true
    end
  end
  if t1 == AbstractGrouping
    if t2 == AbstractGrouping
      return true
    end
  end
  if t1 == TriggeredTask
    if t2 == AbstractOrdering
      return true
    end
    if t2 == AbstractGrouping
      return true
    end
    if t2 == TriggeredTask
      return true
    end
  end
  if t1 == Requirement
    if t2 == Requirement
      return true
    end
  end
  if t1 == InstrumentedFunction
    if t2 == ForSyDeFunction
      return true
    end
    if t2 == InstrumentedFunction
      return true
    end
  end
  if t1 == SDFComb
    if t2 == ForSyDeFunction
      return true
    end
    if t2 == SDFComb
      return true
    end
  end
  if t1 == AbsractInterfaceComponent
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbsractInterfaceComponent
      return true
    end
  end
  if t1 == ReactorActor
    if t2 == ForSyDeFunction
      return true
    end
    if t2 == ReactorElement
      return true
    end
    if t2 == ReactorActor
      return true
    end
  end
  if t1 == SYPrefix
    if t2 == ForSyDeFunction
      return true
    end
    if t2 == SYPrefix
      return true
    end
  end
  if t1 == InstrumentedProcessorTile
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbstractProcessingComponent
      return true
    end
    if t2 == InstrumentedProcessorTile
      return true
    end
    if t2 == Instrumented
      return true
    end
  end
  if t1 == LocationRequirement
    if t2 == Requirement
      return true
    end
    if t2 == LocationRequirement
      return true
    end
  end
  if t1 == MinimumThroughput
    if t2 == Goal
      return true
    end
    if t2 == MinimumThroughput
      return true
    end
  end
  if t1 == SDFPrefix
    if t2 == ForSyDeFunction
      return true
    end
    if t2 == SDFPrefix
      return true
    end
  end
  if t1 == AbstractStorageComponent
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbstractStorageComponent
      return true
    end
  end
  if t1 == Instrumented
    if t2 == Instrumented
      return true
    end
  end
  if t1 == AbstractCommunicationComponent
    if t2 == AbstractPhysicalComponent
      return true
    end
    if t2 == AbstractCommunicationComponent
      return true
    end
  end
  return false
end

@enum EdgeTrait begin
  Output
  Composition
  AbstractScheduling
  AbstractMapping
  AbstractPhysicalConnection
  AbstractDecision
  Annotation
  Input
end

function refines(t1::EdgeTrait, t2::EdgeTrait)
  if t1 == Output
    if t2 == Output
      return true
    end
  end
  if t1 == Composition
    if t2 == Composition
      return true
    end
  end
  if t1 == AbstractScheduling
    if t2 == AbstractScheduling
      return true
    end
    if t2 == AbstractDecision
      return true
    end
  end
  if t1 == AbstractMapping
    if t2 == AbstractMapping
      return true
    end
    if t2 == AbstractDecision
      return true
    end
  end
  if t1 == AbstractPhysicalConnection
    if t2 == AbstractPhysicalConnection
      return true
    end
  end
  if t1 == AbstractDecision
    if t2 == AbstractDecision
      return true
    end
  end
  if t1 == Annotation
    if t2 == Annotation
      return true
    end
  end
  if t1 == Input
    if t2 == Input
      return true
    end
  end
  return false
end

const vertex_trait_lut = Dict(
  "CustomScheduler" => CustomScheduler,
  "InstrumentedCommunicationInterconnect" => InstrumentedCommunicationInterconnect,
  "Goal" => Goal,
  "TimeDivisionMultiplexer" => TimeDivisionMultiplexer,
  "AbstractOrdering" => AbstractOrdering,
  "FixedPriorityScheduler" => FixedPriorityScheduler,
  "BufferSignal" => BufferSignal,
  "ExtraFunctional" => ExtraFunctional,
  "SporadicTask" => SporadicTask,
  "TimeTriggeredScheduler" => TimeTriggeredScheduler,
  "SYComb" => SYComb,
  "ReactorTimer" => ReactorTimer,
  "WCET" => WCET,
  "LabelSignal" => LabelSignal,
  "HardRequirement" => HardRequirement,
  "Signal" => Signal,
  "RoundRobinScheduler" => RoundRobinScheduler,
  "InstrumentedSignal" => InstrumentedSignal,
  "AbstractPhysicalComponent" => AbstractPhysicalComponent,
  "AbstractProcessingComponent" => AbstractProcessingComponent,
  "ForSyDeFunction" => ForSyDeFunction,
  "ReactorElement" => ReactorElement,
  "WCCT" => WCCT,
  "AbstractGrouping" => AbstractGrouping,
  "TriggeredTask" => TriggeredTask,
  "Requirement" => Requirement,
  "InstrumentedFunction" => InstrumentedFunction,
  "SDFComb" => SDFComb,
  "AbsractInterfaceComponent" => AbsractInterfaceComponent,
  "ReactorActor" => ReactorActor,
  "SYPrefix" => SYPrefix,
  "InstrumentedProcessorTile" => InstrumentedProcessorTile,
  "LocationRequirement" => LocationRequirement,
  "MinimumThroughput" => MinimumThroughput,
  "SDFPrefix" => SDFPrefix,
  "AbstractStorageComponent" => AbstractStorageComponent,
  "Instrumented" => Instrumented,
  "AbstractCommunicationComponent" => AbstractCommunicationComponent)

const edge_trait_lut = Dict(
  "Output" => Output,
  "Composition" => Composition,
  "AbstractScheduling" => AbstractScheduling,
  "AbstractMapping" => AbstractMapping,
  "AbstractPhysicalConnection" => AbstractPhysicalConnection,
  "AbstractDecision" => AbstractDecision,
  "Annotation" => Annotation,
  "Input" => Input)

make_trait_vertex(label::AbstractString) = vertex_trait_lut[label]
make_trait_edge(label::AbstractString) = edge_trait_lut[label]
is_trait(label::AbstractString)::Bool = haskey(vertex_trait_lut, label) || haskey(edge_trait_lut, label)
end # module

max_stars_repo_name: forsyde/forsyde-io

max_issues_repo_name: forsyde/forsyde-io

max_forks_repo_name: forsyde/forsyde-io