text:
#---------------------------------------------------------
# # [copyto! test](@id 2-copyto)
#---------------------------------------------------------

# This page examines `copyto!` speed of the lazy grids in the Julia package
# [`LazyGrids`](https://github.com/JuliaArrays/LazyGrids.jl).

# ### Setup

# Packages needed here.

using LazyGrids: ndgrid, ndgrid_array
using LazyGrids: btime, @timeo # not exported; just for timing tests here
using BenchmarkTools: @benchmark

# ### Overview

# There are 4 sub-types of `AbstractGrids`.
# Here we focus on the simplest (using `OneTo`)
# and most general (`AbstractVector`).

# #### `OneTo`

dims = (2^7,2^8,2^9)

(xl, _, _) = ndgrid(dims...) # lazy version
xa = Array(xl) # regular dense Array
out = zeros(Int, dims)
sizeof.((xl, xa))

#
ta = @benchmark copyto!(out, xa) # 12.6ms
btime(ta)

#
tl = @benchmark copyto!(out, xl) # 27.3ms
btime(tl)

# #### `AbstractVector`

x,y,z = map(rand, dims)

(xl, _, _) = ndgrid(dims...)
xa = Array(xl)
out = zeros(eltype(x), dims)
sizeof.((xl, xa))

#
ta = @benchmark copyto!(out, xa) # 15.7ms
btime(ta)

#
tl = @benchmark copyto!(out, xl) # 21.7ms
btime(tl)

# These results suggest that `copyto!` is somewhat slower
# for a lazy grid than for an `Array`.
# This drawback could be reduced or possibly even eliminated
# by adding a dedicated `copyto!` method
# for lazy grids.
# Submit an issue or PR if there is a use case
# that needs faster `copyto!`.
#
# See
# [broadcasting](https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting).


text:
# # Machine Learning in Julia, JuliaCon2020

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

# The following loads a Julia environment and forces pre-compilation
# of some packages.

# However, if this the **binder** notebook version of the tutorial, we
# recommend you *skip* evaluation of this first cell.

DIR = @__DIR__
include(joinpath(DIR, "setup.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)

# ### Advanced

# - [Part 4 - Tuning Hyper-parameters](#part-4-tuning-hyper-parameters)

#

# ## 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 subsampling, no levels are not lost:

levels(exam_mark[1:2])

# **Note on binary data.** There is no separate scientific type for binary
# data. Binary data is `OrderedFactor{2}` if it has an intrinsic
# "true" class (eg, "pass"/"fail") and `Multiclass{2}` otherwise (eg,
# "male"/"female").

# ### 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 a 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)

# Under the hood many algorithms convert tabular data to matrices. If
# your table is a wrapped matrix like the above, then the compiler
# will generally collapse the conversions to a no-op.

# **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 workflow 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
# workflow 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 behaviour 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` workflow.
# > 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):

kitchen_sink = models()

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

meta = kitchen_sink[3]

#-

targetscitype = meta.target_scitype

#-

scitype(y) <: targetscitype

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

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

# 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

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)

#-

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

# 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 `Proababilistic`, 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 workflow, this is not as useful as you
# might imagine. In particular, all probablistic 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(matching(y))

# ### 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
# randomizied 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 workflow 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

# 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 modelling 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. (To convert a vector of pairs `v =
#   [x1 => y1, x2 => y2, ...]` into a dictionary, do `Dict(v)`.)

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

# - ☆(iii) In how many `test` observations does the predicted
#   probablility 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
#   substatially 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 tranformers, 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 `ContinousEncoder` 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 elscityp `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 workflow 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!;
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) |> fit!
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 standarize 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 centred 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 alread 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, "epochs")
ylabel!(plt, "cross entropy on holdout set")

#-

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 `ContiuousEncoder` 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, 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.

# In `RandomSearch` the `scale` attribute of a one-dimensional range
# only plays a role if we specify a *function*, which means we'll need
# to apply the corresponding inverse transform to our bounds, like
# this:

r = range(model, :(logistic_classifier.lambda), lower = -2, upper=2, scale=x->10^x)

# By default, a *bounded* range is sampled uniformly (before the
# `:scale` function is applied). We can see what this means like this:

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

#-

# Alternatively, we can replace `r` with a positive *unbounded* range
# which, by default, is sampled using a `Gamma` distribution (which
# has an infinite decaying tail). An positive unbounded range is specified in
# this way:

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

# And we then get this kind of distribution:

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])
#-

# 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` workflow 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 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.

# ☆(b) Define the one-dimensional range

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

# and try to guess the outcome of evaluating the following two code blocks:

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

#-

sort(unique(samples))

# (c) Optimize these two parameters over the ranges `r1` and `r2`
# using a random search with uniform priors. 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 a
# total 40 models (combinations of hyper-parameters).  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.

# ## 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

ext: jl


text:
#=
Create wrappers for cores that will interact with DES to create the simulation
environment.
=#
mutable struct SimWrapper{T}
    obj :: T
    callback :: Callback
    visits :: Int
end

SimWrapper(obj) = SimWrapper(obj, Callback(), 0)

unwrap(s::SimWrapper) = s.obj
callback(s::SimWrapper) = s.callback
setcallback!(s::SimWrapper, cb) = s.callback = cb

function update!(sim, wrapper :: SimWrapper)
    update!(unwrap(wrapper))
    schedule!(sim, callback(wrapper), clockperiod(unwrap(wrapper)))
    wrapper.visits += 1
    return nothing
end

function wrap!(sim, obj :: T) where T
    wrapper = SimWrapper(obj)
    callback = Callback(sim -> update!(sim, wrapper))
    setcallback!(wrapper, callback)

    return wrapper
end

LightDES.schedule!(sim::Simulation, wrapper::SimWrapper) = schedule!(sim, callback(wrapper), 0)


text:
using DiffEqBase.InternalEuler, DiffEqBase, Test

# Here's the problem to solve

struct LorenzFunction <: Function
    syms::Vector{Symbol}
end

function (::LorenzFunction)(u,p,t)
 [10.0(u[2]-u[1]),u[1]*(28.0-u[3]) - u[2],u[1]*u[2] - (8/3)*u[3]]
end
lorenz = LorenzFunction([:x,:y,:z])

u0 = [1., 5., 10.]
tspan = (0., 100.)
prob = ODEProblem(lorenz, u0, tspan)
dt = 0.1
sol = solve(prob,InternalEuler.FwdEulerAlg(), tstops=0:dt:1)
syms = [:x,:y,:z]

@test DiffEqBase.has_syms(prob.f) == true
@test DiffEqBase.interpret_vars([(0,1), (1,3), (4,5)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([1, (1,3), (4,5)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([1, 3, 4],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,3), (DiffEqBase.DEFAULT_PLOT_FUNC,0,4)]
@test DiffEqBase.interpret_vars(([1,2,3], [4,5,6]),sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,4), (DiffEqBase.DEFAULT_PLOT_FUNC,2,5), (DiffEqBase.DEFAULT_PLOT_FUNC,3,6)]
@test DiffEqBase.interpret_vars((1, [2,3,4]),sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,1,4)]

@test DiffEqBase.interpret_vars([(:t,:x),(:t,:y)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,2)]
@test DiffEqBase.interpret_vars([:x, (0,:x), (:x,:y)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,2)]
@test DiffEqBase.interpret_vars([:x, :y, :z],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,2), (DiffEqBase.DEFAULT_PLOT_FUNC,0,3)]
@test DiffEqBase.interpret_vars(([:x,:x], [:y,:z]),sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3)]
@test DiffEqBase.interpret_vars((:x, [:y,:z]),sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3)]

f(x,y) = (x+y,y)
@test DiffEqBase.interpret_vars([(f,0,1), (1,3), (4,5)],sol,syms) == [(f,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([1, (f,1,3), (4,5)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (f,1,3), (DiffEqBase.DEFAULT_PLOT_FUNC,4,5)]
@test DiffEqBase.interpret_vars([(f,:t,:x),(:t,:y)],sol,syms) == [(f,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,0,2)]
@test DiffEqBase.interpret_vars([:x, (f,0,:x), (:x,:y)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,0,1), (f,0,1), (DiffEqBase.DEFAULT_PLOT_FUNC,1,2)]
@test DiffEqBase.interpret_vars([(:x,:y)],sol,syms) == [(DiffEqBase.DEFAULT_PLOT_FUNC,1,2)]
@test DiffEqBase.interpret_vars((f,:x,:y),sol,syms) == [(f,1,2)]


text:
# Equality

# equals() now depends on == instead
# of the other way round.
function ==(x::Decimal, y::Decimal)
    # return early on zero
    x_is_zero = iszero(x)
    y_is_zero = iszero(y)
    if x_is_zero || y_is_zero
        return x_is_zero === y_is_zero
    end

    a = normalize(x)
    b = normalize(y)
    a.c == b.c && a.q == b.q && a.s == b.s
end

Base.iszero(x::Decimal) = iszero(x.c)

function <(x::Decimal, y::Decimal)
    # return early on zero
    if iszero(x) && iszero(y)
        return false
    end

    # avoid normalization if possible
    if x.q == y.q
        return isless(x.s == 0 ? x.c : -x.c, y.s == 0 ? y.c : -y.c)
    end

    diff = y - x

    farther_from_0 = diff.c > 0 || (iszero(diff.c) && diff.q > 0)

    if diff.s == 1
        return !farther_from_0
    else
        return farther_from_0
    end
end

function <=(x::Decimal, y::Decimal)
    return x < y || x == y
end

# Special case equality with AbstractFloat to allow comparison against Inf/Nan
# which are not representable in Decimal

==(a::AbstractFloat, b::Decimal) = b == a
function ==(a::Decimal, b::AbstractFloat)
    # Decimal does not represent NaN/Inf
    (isinf(b) || isnan(b)) && return false
    ==(promote(a, b)...)
end

function Base.min(a::AbstractFloat, b::Decimal)
    !signbit(a) && isinf(a) && return convert(promote_type(typeof(a), typeof(b)), b)
    min(promote(a, b)...)
end
function Base.min(a::Decimal, b::AbstractFloat)
    !signbit(b) && isinf(b) && return convert(promote_type(typeof(a), typeof(b)), a)
    min(promote(a, b)...)
end
Base.min(a::Decimal, b::Decimal) = invoke(min, Tuple{T, T} where T<:AbstractFloat, a, b)

function Base.max(a::AbstractFloat, b::Decimal)
    signbit(a) && isinf(a) && return convert(promote_type(typeof(a), typeof(b)), b)
    max(promote(a, b)...)
end
function Base.max(a::Decimal, b::AbstractFloat)
    signbit(b) && isinf(b) && return convert(promote_type(typeof(a), typeof(b)), a)
    max(promote(a, b)...)
end
Base.max(a::Decimal, b::Decimal) = invoke(max, Tuple{T, T} where T<:AbstractFloat, a, b)


text:
using Flight
using OrdinaryDiffEq
using SciMLBase
using LinearAlgebra
using BenchmarkTools
using Base.Iterators

################################## AC2 ######################################

function demo_rt()

    h_trn = AltOrth(609.5);

    trn = HorizontalTerrain(altitude = h_trn);
    atm = System(AtmosphereDescriptor());
    ac = System(C172Aircraft());
    ac_mdl = Model(ac, (trn, atm); dt = 0.02, t_end = 20, adaptive = false, solver = RK4(), y_saveat = 0.02);

    kin_init = KinInit(v_eOb_b = [30, 0, 0],
                        ω_lb_b = [0, 0, 0],
                        q_nb = REuler(ψ = π, θ = 0.1, φ = 0.05),
                        Ob = Geographic(LatLon(ϕ = deg2rad(40.531818), λ = deg2rad(-3.574862)),
                                        h_trn + 0.85 + 0.4));
    ac.u.throttle = 0
    ac.u.pedals = 0.1
    ac.u.brake_left = 1
    ac.u.brake_right = 1
    atm.u.wind.v_ew_n[1] = 0

    # kin_init = KinInit(v_eOb_b = [0, 0, 0],
    #                     ω_lb_b = [0, 0, 0],
    #                     q_nb = REuler(ψ = π, θ = 0, φ = 0),
    #                     Ob = Geographic(LatLon(ϕ = deg2rad(40.531818), λ = deg2rad(-3.574862)),
    #                                     h_trn + 0.85));

    # ac.u.throttle = 1
    # ac.u.brake_left = 1
    # ac.u.brake_right = 1

    init!(ac, kin_init)
    #if the model was instantiated in advance, we need this to change its internal state vector!
    reinit!(ac_mdl, ac.x)

    xp = XPInterface()
    disable_physics(xp)
    set_position(xp, kin_init)

    output_div = 1

    t_wall = time()
    t_wall_0 = t_wall

    println("Generalize this for non-rt")

    # for i in take(ac_mdl.integrator, 5)
    for i in ac_mdl.integrator

        #when using the integrator as an iterator, we don't need to step, it is done
        #automatically at the beginning of each iteration. therefore, the dt we need
        #for pacing the simulation is not the proposed dt for the next step, but the
        #one in the step the integrator has already taken.

        #retrieve the dt just taken by the integrator
        dt = ac_mdl.dt

        #compute the wall time corresponding to the newly updated simulation
        t_wall_next = t_wall + dt

        #busy wait until the wall time catches up
        while (time() < t_wall_next) end
        # println(time()-t_wall_next)

        t_wall = t_wall_next
        # println(t_wall - t_wall_0)

        if ac_mdl.success_iter % output_div == 0
            set_position(xp, ac_mdl.sys.y.kin.pos)
        end

        #when GLFW is implemented, build the model with t_end = Inf and use this to break
        # if !GLFW.WindowShouldClose(window)
        #     break
        # end

    end

    plot_settings = (linewidth=2, margin = 10mm, guidefontsize = 12)
    plots(ac_mdl; save_path = joinpath("tmp", "plots_demo_rt"), plot_settings...)

end

# plots(ac_mdl; save_path = joinpath("tmp", "plots_ac_rt"), plot_settings...)


text:
# ******************************************************************************************
# Notices:
#
# Copyright © 2022 United States Government as represented by the Administrator of the
# National Aeronautics and Space Administration.  All Rights Reserved.
#
# Disclaimers
#
# No Warranty: THE SUBJECT SOFTWARE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY OF ANY KIND,
# EITHER EXPRESSED, IMPLIED, OR STATUTORY, INCLUDING, BUT NOT LIMITED TO, ANY WARRANTY
# THAT THE SUBJECT SOFTWARE WILL CONFORM TO SPECIFICATIONS, ANY IMPLIED WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR FREEDOM FROM INFRINGEMENT, ANY
# WARRANTY THAT THE SUBJECT SOFTWARE WILL BE ERROR FREE, OR ANY WARRANTY THAT DOCUMENTATION,
# IF PROVIDED, WILL CONFORM TO THE SUBJECT SOFTWARE. THIS AGREEMENT DOES NOT, IN ANY MANNER,
# CONSTITUTE AN ENDORSEMENT BY GOVERNMENT AGENCY OR ANY PRIOR RECIPIENT OF ANY RESULTS,
# RESULTING DESIGNS, HARDWARE, SOFTWARE PRODUCTS OR ANY OTHER APPLICATIONS RESULTING FROM
# USE OF THE SUBJECT SOFTWARE.  FURTHER, GOVERNMENT AGENCY DISCLAIMS ALL WARRANTIES AND
# LIABILITIES REGARDING THIRD-PARTY SOFTWARE, IF PRESENT IN THE ORIGINAL SOFTWARE, AND
# DISTRIBUTES IT "AS IS."
#
# Waiver and Indemnity:  RECIPIENT AGREES TO WAIVE ANY AND ALL CLAIMS AGAINST THE UNITED
# STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL AS ANY PRIOR RECIPIENT.
# IF RECIPIENT'S USE OF THE SUBJECT SOFTWARE RESULTS IN ANY LIABILITIES, DEMANDS, DAMAGES,
# EXPENSES OR LOSSES ARISING FROM SUCH USE, INCLUDING ANY DAMAGES FROM PRODUCTS BASED ON,
# OR RESULTING FROM, RECIPIENT'S USE OF THE SUBJECT SOFTWARE, RECIPIENT SHALL INDEMNIFY AND
# HOLD HARMLESS THE UNITED STATES GOVERNMENT, ITS CONTRACTORS AND SUBCONTRACTORS, AS WELL
# AS ANY PRIOR RECIPIENT, TO THE EXTENT PERMITTED BY LAW.  RECIPIENT'S SOLE REMEDY FOR ANY
# SUCH MATTER SHALL BE THE IMMEDIATE, UNILATERAL TERMINATION OF THIS AGREEMENT.
# ******************************************************************************************

"""
Abstract base type for solvers.
"""
abstract type Solver end

"""
Abstract type for solvers that return an adversarial policy.
"""
abstract type GlobalSolver <: Solver end

"""
Abstract type for solvers that return failure examples.
"""
abstract type LocalSolver <: Solver end

"""
Abstract type for solver output.
"""
abstract type Result end

"""
Abstract type for global result (adversarial policy).
"""
abstract type GlobalResult <: Result end

"""
Abstract type for local result (failure examples).
"""
abstract type LocalResult <: Result end

function (solver::Solver)(env::Any)
    @warn """Solver was not passed an environment generator.
    Suboptimal performance may occur if environment cannot be auto-replicated.
    """
    return solver(() -> env)
end

"""
    solve(::Solver, ::Function)

Apply solver to MDP generator function.
"""
function solve(::Solver, ::Function) end
(solver::Solver)(env_fn::Function) = solve(solver, env_fn)


text:
#
# types.jl -
#
# Type definitions for Julia interface to ActiveSilicon Phoenix (PHX) library.
#
#------------------------------------------------------------------------------
#
# This file is part of the `Phoenix.jl` package which is licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2017-2019, Éric Thiébaut (https://github.com/emmt/Phoenix.jl).
# Copyright (C) 2016, Éric Thiébaut & Jonathan Léger.
#

# Custom exception to report errors.
struct PHXError <: Exception
   status::Status
end

struct TimeVal
    sec::_typeof_tv_sec    # seconds
    usec::_typeof_tv_usec  # microseconds
end

struct TimeSpec
    sec::_typeof_tv_sec    # seconds
    nsec::_typeof_tv_nsec  # nanoseconds
end

"""

`AcquisitionContext` is used to share information between the acquisition
callback and the main Julia thread.  It has the following fields:

- `mutex` is a lock to protect the contents of the structure;

- `cond` is a condition variable to signal events to the processing thread;

- `events` is the mask of events to be signaled to the processing thread;

- `number` is the total number frames so far (it tries to take into account
  all the frames: captured frames but also overflows and synchronization losts);

- `overflows`: is the number of overflows (or dropped frames) so far;

- `synclosts` is the number of synchronization losts so far

- `pending` is the number of image buffers ready for being processed;

See also: [`read`](@ref), [`start`](@ref), [`wait`](@ref).

"""
mutable struct AcquisitionContext
    mutex::Ptr{Cvoid}     # Lock to protect this structure
    cond::Ptr{Cvoid}      # Condition to signal events
    # The 2 following fields should exactly match `ImageBuff` structure
    imgbuf::Ptr{Cvoid}    # Address of last captured image buffer
    imgctx::Ptr{Cvoid}    # Address of context associated with last ...
    index::Int            # Index of last captured image buffer
    # The 2 following fields should exactly match `TimeVal` structure
    sec::_typeof_tv_sec   # Time stamp (seconds) of last captured image
    usec::_typeof_tv_usec # Time stamp (microseconds) of last captured image
    number::Int           # Number of image buffers so far
    overflows::Int        # Number of overflows so far
    synclosts::Int        # Number of synchronization losts so far
    pending::Int          # Number of pending image buffers
    events::UInt          # Mask of events to be signaled
    function AcquisitionContext()
        ctx = new(C_NULL, C_NULL, C_NULL, C_NULL,
                  0, 0, 0, 0, 0, 0, 0, 0)
        ctx.mutex = Libc.malloc(_sizeof_pthread_mutex_t)
        if ctx.mutex == C_NULL
            throw(OutOfMemoryError())
        end
        if ccall(:pthread_mutex_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                 ctx.mutex, C_NULL) != SUCCESS
            _destroy(ctx)
            error("pthread_mutex_init failed")
        end
        ctx.cond = Libc.malloc(_sizeof_pthread_cond_t)
        if ctx.cond == C_NULL
            _destroy(ctx)
            throw(OutOfMemoryError())
        end
        if ccall(:pthread_cond_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}),
                 ctx.cond, C_NULL) != SUCCESS
            _destroy(ctx)
            error("pthread_cond_init failed")
        end
        return finalizer(_destroy, ctx)
    end
end

# Beware must not be mutable!
struct FrameData
    sec::_typeof_tv_sec   # Time stamp (seconds)
    usec::_typeof_tv_usec # Time stamp (microseconds)
    index::Int
end

"""

`_destroy(obj)` is the finalizer method for a Phoenix camera instance or an
acquisition context.  It *must not* be called directly.

See also: [`Phoenix.Camera`](@ref), [`Phoenix.AcquisitionContext`](@ref).

"""
function _destroy(ctx::AcquisitionContext)
    # This is the destructor of an AcquisitionContext instance.
    if ctx.cond != C_NULL
        cond = ctx.cond
        ctx.cond = C_NULL
        ccall(:pthread_cond_destroy, Cint, (Ptr{Cvoid},), cond)
        Libc.free(cond)
    end
    if ctx.mutex != C_NULL
        mutex = ctx.mutex
        ctx.mutex = C_NULL
        ccall(:pthread_mutex_destroy, Cint, (Ptr{Cvoid},), mutex)
        Libc.free(mutex)
    end
end

"""

Concrete types derived from abstract type `CameraModel` are used to uniquely
identify the different camera models.

"""
abstract type CameraModel end

abstract type Configuration{M<:CameraModel} end

mutable struct Camera{M<:CameraModel} <: ScientificCamera
    state::Int # 0 initially, 1 when camera open, 2 while acquiring
    handle::Handle # 0 before calling PHX_Create or after calling PHX_Destroy
    imgs::Vector{Array{T,2}} where {T} # images for acquisition
    ctxs::Vector{FrameData} # metadata for captured images
    bufs::Vector{ImageBuff} # virtual image buffers currently used
    context::AcquisitionContext # context shared with acquisition callback
    timeout::UInt32 # time out (in ms) for reading/writing registers
    swap::Bool # swap bytes for read/write control connection?
    coaxpress::Bool # is it a CoaXPress camera?
    debug::Bool # print debug messages?

    function Camera{M}(errorhandler::Ptr{Cvoid} = _ERRORHANDLER_REF[]) where {M}
        # Create a new PHX handle structure.
        handle = Ref{Handle}(0)
        status = ccall(_PHX_Create[], Status, (Ptr{Handle}, Ptr{Cvoid}),
                       handle, errorhandler)
        status == PHX_OK || throw(PHXError(status))

        # Create the instance and attach the destroy callback.
        cam = new{M}(0, handle[],
                     Vector{Array{UInt8,2}}(undef, 0),
                     Vector{FrameData}(undef, 0),
                     Vector{ImageBuff}(undef, 0),
                     AcquisitionContext(),
                     500, false, false, false)
        return finalizer(_destroy, cam)
    end
end

function _destroy(cam::Camera)
    if cam.handle != 0
        if cam.state == 2
            # Abort acquisition (using the private routine which does not throw
            # exceptions).
            _readstream(cam, PHX_ABORT, C_NULL)
            _readstream(cam, PHX_UNLOCK, C_NULL)
            _stophook(cam)
            cam.state = 1
        end
        if cam.state == 1
            # Close the camera.
            cam.state = 0
            _close(cam)
        end
        # Release other ressources.
        ref = Ref(cam.handle)
        cam.handle = 0 # to avoid doing this more than once
        ccall(_PHX_Destroy[], Status, (Ptr{Handle},), ref)
    end
end


text:
using Luxor, Random

function juliatables1(cpos::Point, radius, twidth, theight, nrows, ncols)
    # clockwise, from bottom LEFT...
    gsave()
    points3 = ngon(cpos, radius, 3, pi/6, vertices=true)
    theta = 0
    for (n, centerp) in enumerate(points3)
        gsave()
        translate(centerp)
        tiles = Tiler(twidth, theight, nrows, ncols)
        for (pos, p) in tiles
            sethue("black")
            if tiles.currentrow == 1 && tiles.currentcol == 1

            elseif tiles.currentrow == 1
                sethue("grey70")
                box(pos, tiles.tilewidth, tiles.tileheight, :stroke)
                sethue(cols[mod1(tiles.currentcol, 3)])
                circle(pos, 3, :fill)
            elseif tiles.currentcol == 1
                sethue("grey70")
                box(pos, tiles.tilewidth, tiles.tileheight, :stroke)
                box(pos, 6, 4, :fill)
            elseif tiles.currentrow > 1
                box(pos, tiles.tilewidth - 2, tiles.tileheight - 2, :stroke)
                sethue(cols[mod1(tiles.currentcol, 3)])
                box(pos, tiles.tilewidth - 6, tiles.tileheight - 6, :fill)
            end
        end
        grestore()
    end
    grestore()
end

function drawflask(pos;
    scalefactor=1)
    flask = Point[Point(100, 0), Point(25, -140), Point(25, -220),
        Point(-25, -220), Point(-25, -140), Point(-100, 0)]
    @layer begin
        translate(pos + (0, 100))
        scale(scalefactor)
        sethue("black")
        setline(20)
        #polysmooth(flask, 10, :fill)
        polysmooth(flask, 10, :clip)
        positions = between.(
            boxtopcenter(BoundingBox(flask)),
            boxbottomcenter(BoundingBox(flask)),
            range(0.5, 1.0, length=3))
        @layer begin
            setopacity(0.9)
            for i in 1:3
                sethue(cols[mod1(i, end)])
                box(positions[i], 200, boxheight(BoundingBox(flask))/2, :fill)
            end
            setopacity(0.3)
            for i in 1:15
                sethue("white")
                circle(rand(BoundingBox(flask) * 0.5), 4, :fill)
            end
        end
        clipreset()
        sethue("black")
        polysmooth(flask, 10, :stroke)
    end
end

const cols = [Luxor.julia_green, Luxor.julia_red,Luxor.julia_purple, Luxor.julia_blue]

function datascienceicon()
    Random.seed!(1)
    sethue("antiquewhite")
    squircle(O, 240, 240, :fill, rt=0.1)

    sethue("plum1")
    circle(O, 200, :fill)

    @layer begin
        setopacity(0.98)
        scale(0.8)
        setline(0.5)
        # cpos, radius, twidth, theight, nrows, ncols)
        juliatables1(Point(-120, -100), 100, 150, 150, 7, 7)
    end
    drawflask(O + (30, 100), scalefactor=1.7)
end

function drawicon()
    Drawing(500, 500, "/tmp/juliadatascienceicon.svg")
    origin()
    datascienceicon()
    finish()
    preview()
end

function drawlogo()
    Drawing(960, 540, "/tmp/juliadatascienceandtitle.svg")
    Luxor.background("ivory")
    origin()
    table = Table([540], [540, 960-540])
    translate(table[1])
    datascienceicon()
    translate(table[2])
    sethue("black")
    fontface("MacklinVariable_Slab-Black")
    fontsize(120)
    Luxor.textbox(["Data", "Science"], O - (0, 80), leading=90)
    finish()
    preview()
end

drawicon()
drawlogo()

max_stars_repo_name: radovankavicky/graphics

max_issues_repo_name: radovankavicky/graphics

max_forks_repo_name: cormullion/graphics


text:
# ---
# title: 523. Continuous Subarray Sum
# id: problem523
# author: Indigo
# date: 2021-06-02
# difficulty: Medium
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/continuous-subarray-sum/description/>
# hidden: true
# ---
#
# Given a list of **non-negative** numbers and a target **integer** k, write a
# function to check if the array has a continuous subarray of size at least 2
# that sums up to a multiple of **k** , that is, sums up to n*k where n is also
# an **integer**.
#
#
#
# **Example 1:**
#
#
#
#     Input: [23, 2, 4, 6, 7],  k=6
#     Output: True
#     Explanation: Because [2, 4] is a continuous subarray of size 2 and sums up to 6.
#
#
# **Example 2:**
#
#
#
#     Input: [23, 2, 6, 4, 7],  k=6
#     Output: True
#     Explanation: Because [23, 2, 6, 4, 7] is an continuous subarray of size 5 and sums up to 42.
#
#
#
#
# **Constraints:**
#
#   * The length of the array won't exceed 10,000.
#   * You may assume the sum of all the numbers is in the range of a signed 32-bit integer.
#
#
## @lc code=start
using LeetCode

function check_subarray_sum(nums::Vector{Int}, k::Int)
    pushfirst!(nums, 0)
    cumsum!(nums, nums)
    nums .%= k
    st = Set{Int}()
    len = length(nums)
    for i in 3:len
        nums[i] ∈ st && return true
        push!(st, nums[i - 1])
    end
    return false
end
## @lc code=end


text:
function triangle_vbroadcast_quote(V::DataType, Alength::Int)
    q = quote end
    outtup = Expr(:tuple,)
    for m ∈ 1:Alength
        Asym = Symbol(:A_, m)
        push!(q.args, :($Asym = vbroadcast($V, L[$m])))
        push!(outtup.args, Asym)
    end
    q, outtup
end

# @generated
@generated function SIMDPirates.vbroadcast(::Type{Vec{W,T}}, L::AbstractUpperTriangularMatrix{M,T}) where {W,T,M}
    Alength = (M*(M+1)) >> 1
    V = Vec{W,T}
    q, outtup = triangle_vbroadcast_quote(V, Alength)
    # push!(q.args, :($(A.name){$M,$V,$Alength}($outtup) ) )
    quote
        @inbounds begin
            $q
        end
        $(Expr(:call, Expr(:curly, UpperTriangularMatrix, M, V, Alength), outtup))
        # $(A.name){$M,$V,$Alength}($outtup)
    end
end
@generated function SIMDPirates.vbroadcast(::Type{Vec{W,T}}, L::AbstractLowerTriangularMatrix{M,T}) where {W,T,M}
    Alength = (M*(M+1)) >> 1
    V = Vec{W,T}
    q, outtup = triangle_vbroadcast_quote(V, Alength)
    # push!(q.args, :($(A.name){$M,$V,$Alength}($outtup) ) )
    quote
        @inbounds begin
            $q
        end
        $(Expr(:call, Expr(:curly, LowerTriangularMatrix, M, V, Alength), outtup))
        # $(A.name){$M,$V,$Alength}($outtup)
    end
end

@generated function Base.zero(::Type{LT}) where {N,W,T, LT <: AbstractMutableDiagMatrix{N,Vec{W,T}}}
    quote
        $(Expr(:meta,:inline))
        out = $LT(undef)
        z = vbroadcast(Vec{$W,$T}, zero($T))
        @inbounds for i ∈ 1:$(binomial2(N+1))
            out[i] = z
        end
        out
    end
end

function triangle_vsum_quote(T::DataType,M::Int,L::Int,Lfull::Int)
    quote
        @inbounds for l ∈ 1:$L
            out[l] = SIMDPirates.vsum(triangle[l])
        end
        @inbounds for l ∈ $(L+1):$Lfull
            out[l] = zero($T)
        end
    end
end
@generated function SIMDPirates.vsum(triangle::AbstractLowerTriangularMatrix{M,Vec{W,T},L}) where {T,W,M,L}
    Wm1 = VectorizationBase.pick_vector_width(L, T) - 1
    Lfull = (L + Wm1) & ~Wm1
    q = triangle_vsum_quote(T,M,L,Lfull)
    quote
        $(Expr(:meta,:inline))
        out = MutableLowerTriangularMatrix{$M,$T,$Lfull}(undef)
        $q
        out
    end
end
@generated function SIMDPirates.vsum(triangle::AbstractUpperTriangularMatrix{M,Vec{W,T},L}) where {T,W,M,L}
    Wm1 = VectorizationBase.pick_vector_width(L, T) - 1
    Lfull = (L + Wm1) & ~Wm1
    q = triangle_vsum_quote(T,M,L,Lfull)
    quote
        $(Expr(:meta,:inline))
        out = MutableUpperTriangularMatrix{$M,$T,$Lfull}(undef)
        $q
        out
    end
end
@inline function vsumvec(triangle::AbstractDiagTriangularMatrix{M,Vec{W,T},L}) where {W,T,M,L}
    out = FixedSizeVector{L,T}(undef)
    @inbounds for l ∈ 1:L
        out[l] = SIMDPirates.vsum(triangle[l])
    end
    out
end

# function triangle_vsum_quote(T::DataType,M::Int,L::Int,Lfull::Int)
#     q = quote end
#     outtup = Expr(:tuple,)
#     for l ∈ 1:L
#         sym = Symbol(:triangle_, l)
#         push!(q.args, :($sym = SIMDPirates.vsum(triangle[$l])))
#         push!(outtup.args, sym)
#     end
#     for l ∈ L+1:Lfull
#         push!(outtup.args, zero(T))
#     end
#     q, outtup
# end
#
# @generated function SIMDPirates.vsum(triangle::AbstractLowerTriangularMatrix{M,Vec{W,T},L}) where {W,T,M,L}
#     Lfull = VectorizationBase.pick_vector_width(L, T)
#     q, outtup = triangle_vsum_quote(T,M,L,Lfull)
#     quote
#         @inbounds begin
#             $q
#         end
#         $(Expr(:call, Expr(:curly, LowerTriangularMatrix, M, T, Lfull), outtup))
#     end
# end
# @generated function SIMDPirates.vsum(triangle::AbstractUpperTriangularMatrix{M,Vec{W,T},L}) where {W,T,M,L}
#     Lfull = VectorizationBase.pick_vector_width(L, T)
#     q, outtup = triangle_vsum_quote(T,M,L,Lfull)
#     quote
#         @inbounds begin
#             $q
#         end
#         $(Expr(:call, Expr(:curly, UpperTriangularMatrix, M, T, Lfull), outtup))
#     end
# end

function mkernel_triangle_AU_quote(kernel_size_m::Int, kernel_size_n::Int, mindexpr, nindexpr, add::Bool = true, W::Int = 0, T::DataType = Float64)
    muladdfunc = add ? :vmuladd : :vfnmadd

    # kernel_size = kernel_size_n #min(kernel_size_m, kernel_size_n)
    # @assert kernel_size == kernel_size_m
    # Assumes A_m_n block is already loaded

    q = quote end
    for nouter in 1:kernel_size_n
        for m in 1:kernel_size_m
            push!(q.args, :( $(Symbol(:AU_,m)) = AU[$m + $mindexpr, $nindexpr + $nouter ] ) )
        end
        for ninner in 1:nouter
            if ninner == nouter
                Uind = :(U[ $nindexpr + $ninner ])
            else
                Uind = :(U[ triangle_ind_tuple[$nouter + $nindexpr] + $ninner + $nindexpr ])
            end
            W == 0 ? push!(q.args, :(vU = $Uind )) : push!(q.args, :(vU = vbroadcast(Vec{$W,$T}, $Uind) ))
            for m in 1:kernel_size_m
                AUsym = Symbol(:AU_,m)
                push!(q.args, :( $AUsym = SIMDPirates.$muladdfunc($(Symbol(:A_,m,:_,ninner)), vU, $AUsym) ) )
            end
        end
        for m in 1:kernel_size_m
            push!(q.args, :( AU[$m + $mindexpr, $nindexpr + $nouter ] =  $(Symbol(:AU_,m)) ) )
        end
    end
    q
end
function square_iterations_AU_quote(kernel_size_m::Int, kernel_size_n::Int, mindexpr, nindexpr, N::Int, add::Bool, W::Int, T::DataType)
    muladdfunc = add ? :vmuladd : :vfnmadd
    vUind = :(U[ $nindexpr + n + triangle_ind_tuple[ncol_iter] ])
    vUexpr = W == 0 ? vUind : :(vbroadcast(Vec{$W,$T}, $vUind))
    quote
        # Base.Cartesian.@nexprs $kernel_size_n n -> Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[$mindexpr+m,$nindexpr+n]
        for ncol_iter in $(nindexpr)+$(kernel_size_n+1):$N
            Base.Cartesian.@nexprs $kernel_size_m m -> AU_m = AU[m + $mindexpr, ncol_iter]
            Base.Cartesian.@nexprs $kernel_size_n n -> begin
                vU = $vUexpr
                Base.Cartesian.@nexprs $kernel_size_m m -> begin
                    AU_m = SIMDPirates.$muladdfunc(A_m_n, vU, AU_m)
                end
            end
            Base.Cartesian.@nexprs $kernel_size_m m -> AU[m + $mindexpr, ncol_iter] = AU_m
        end
    end
end
function mkern_iteration_AU_quote(kernel_size_m::Int, kernel_size_n::Int, nk::Int, nr::Int, mindexpr, N::Int, add::Bool, W::Int, T::DataType)
    nindexpr = :(nkern * $kernel_size_n)
    if nr == 0
        q = quote
            for nkern in 0:$(nk-1)
                Base.Cartesian.@nexprs $kernel_size_n n -> begin
                    Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[m+$mindexpr,n+$nindexpr]
                end
                $(mkernel_triangle_AU_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, add, W, T))
                $(square_iterations_AU_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, N, add, W, T))
            end
        end
    else
        # nindexpr = :($nr + nkern * $kernel_size_n)
        q = quote
            for nkern in 0:$(nk-1)
                Base.Cartesian.@nexprs $kernel_size_n n -> begin
                    Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[m+$mindexpr,n+$nindexpr]
                end
                $(mkernel_triangle_AU_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, add, W, T))
                $(square_iterations_AU_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, N, add, W, T))
            end
            Base.Cartesian.@nexprs $nr n -> begin
                Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[m+$mindexpr,n+$(nk*kernel_size_n)]
            end
            $(mkernel_triangle_AU_quote(kernel_size_m, nr, mindexpr, nk * kernel_size_n, add, W, T))
        end
    end
    q
end

function matrix_upper_triangle_mul_quote(M::Int, N::Int, T::DataType, add::Bool = true, W::Int = 0)
    mk, mr, nk, nr, kernel_size_m, kernel_size_n = PaddedMatrices.determine_pattern(M,N)

    if nr > 0
        # @assert nk > 0
        # nk -= 1
        # nr = kernel_size_n
        if (N % (nk+1)) == 0
            nr = 0
            nk += 1
            kernel_size_n = N ÷ nk
        end
    end
    if mr > 0
        if (M % (mk+1)) == 0
            mr = 0
            mk += 1
            kernel_size_m = M ÷ mk
        end
    end
    q_triangle_ind_defs = Expr(:tuple, Int32(N))
    triangle_ind = N
    triangle_increment = 0
    for n in 2:N
        push!(q_triangle_ind_defs.args, Int32(triangle_ind))
        triangle_increment += 1
        triangle_ind += triangle_increment
    end

    q = quote
        triangle_ind_tuple = $q_triangle_ind_defs
        for mkern in 0:$(mk-1)
            $(mkern_iteration_AU_quote(kernel_size_m, kernel_size_n, nk, nr, :(mkern * $kernel_size_m), N, add, W, T))
        end
    end
    if mr > 0
        push!(q.args, mkern_iteration_AU_quote(mr, kernel_size_n, nk, nr, mk * kernel_size_m, N, add, W, T))
    end
    q
end

@generated function PaddedMatrices.zero!(U::AbstractDiagTriangularMatrix{M,NTuple{W,Core.VecElement{T}}}) where {M,W,T}
    L = binomial2(M+1)
    quote
        $(Expr(:meta,:inline))
        @inbounds for l ∈ 1:$L
            U[l] = SIMDPirates.vbroadcast(Vec{$W,$T}, zero($T))
        end
        nothing
    end
end

@generated function Base.:*(
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}
        ) where {M,N,T,W,L}

    quote
        $(Expr(:meta,:inline))
        AU = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        @inbounds for i in 1:$(M*N)
            AU[i] = vbroadcast(Vec{$W,$T}, zero($T))
        end
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, true, 0))
        end
        AU
    end
end
@generated function Base.:*(
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,T,L}
        ) where {M,N,T,W,L}

    quote
        $(Expr(:meta,:inline))
        AU = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        @inbounds for i in 1:$(M*N)
            AU[i] = vbroadcast(Vec{$W,$T}, zero($T))
        end
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, true, W))
        end
        AU
    end
end

@generated function LinearAlgebra.mul!(
            AU::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        PaddedMatrices.zero!(AU)
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, true, 0))
        end
        AU
    end
end
@generated function LinearAlgebra.mul!(
            AU::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,T,L}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        PaddedMatrices.zero!(AU)
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, true, W))
        end
        AU
    end
end
@generated function addmul!(
            AU::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, true, 0))
        end
        AU
    end
end
@generated function addmul!(
            AU::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,T,L}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, true, W))
        end
        AU
    end
end

@generated function submul!(
            AU::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}
        ) where {M,N,W,T,L}

    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, false, 0))
        end
        AU
    end
end
@generated function submul!(
            AU::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            U::AbstractUpperTriangularMatrix{N,T,L}
        ) where {M,N,W,T,L}

    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_upper_triangle_mul_quote(M, N, T, false, W))
        end
        AU
    end
end

function mkernel_end_AUt_quote(kernel_size_m::Int, kernel_size_n::Int, mindexpr, nindexpr, add::Bool = true, W::Int = 0, T::DataType = Float64)
    muladdfunc = add ? :vmuladd : :vfnmadd

    # kernel_size = min(kernel_size_m, kernel_size_n)
    # @assert kernel_size == kernel_size_m
    # Assumes A_m_n block is already loaded
    q = quote end
    for nouter in 1:kernel_size_n
        for m in 1:kernel_size_m
            push!(q.args, :( $(Symbol(:AUt_,m)) = AUt[$m + $mindexpr, $nindexpr + $nouter ] ) )
        end
        for ninner in nouter:kernel_size_n
            if ninner == nouter
                Utind = :(Ut[ $nindexpr + $nouter ])
            else
                Utind = :(Ut[ $(Symbol(:triangle_ind_,  ninner)) + $nindexpr + $nouter ] )
            end
            W == 0 ? push!(q.args, :(vUt = $Utind)) : push!(q.args, :(vUt = vbroadcast(Vec{$W,$T}, $Utind)))
            for m in 1:kernel_size_m
                AUtsym = Symbol(:AUt_,m)
                push!(q.args, :( $AUtsym = SIMDPirates.$muladdfunc($(Symbol(:A_,m,:_,ninner)), vUt, $AUtsym) ) )
            end
        end
        for m in 1:kernel_size_m
            push!(q.args, :( AUt[$m + $mindexpr, $nindexpr + $nouter ] =  $(Symbol(:AUt_,m)) ) )
        end
    end
    q
end
function square_iterations_AUt_quote(kernel_size_m::Int, kernel_size_n::Int, mindexpr, nindexpr, add::Bool, W::Int, T::DataType)
    muladdfunc = add ? :vmuladd : :vfnmadd
    vUtexpr = W == 0 ? :(Ut[ triangle_ind_n + ncol_iter ]) : :(vbroadcast(Vec{$W,$T}, Ut[ triangle_ind_n + ncol_iter ]))
    quote
        # Base.Cartesian.@nexprs $kernel_size_n n -> Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[$mindexpr+m,$nindexpr+n]
        for ncol_iter in 1:$(nindexpr)
            Base.Cartesian.@nexprs $kernel_size_m m -> AUt_m = AUt[m + $mindexpr, ncol_iter]
            Base.Cartesian.@nexprs $kernel_size_n n -> begin
                vUt = $vUtexpr
                Base.Cartesian.@nexprs $kernel_size_m m -> begin
                    AUt_m = SIMDPirates.$muladdfunc(A_m_n, vUt, AUt_m)
                end
            end
            Base.Cartesian.@nexprs $kernel_size_m m -> AUt[m + $mindexpr, ncol_iter] = AUt_m
        end
    end
end
function mkern_iteration_AUt_quote(kernel_size_m::Int, kernel_size_n::Int, nk::Int, nr::Int, mindexpr, add::Bool, W::Int, T::DataType)
    if nr == 0
        nindexpr = :(nkern * $kernel_size_n)
        q = quote
            for nkern in 0:$(nk-1)
                Base.Cartesian.@nexprs $kernel_size_n n -> begin
                    triangle_ind_n = triangle_ind_tuple[n + $nindexpr]
                    Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[m+$mindexpr,n+$nindexpr]
                end
                $(square_iterations_AUt_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, add, W, T))
                $(mkernel_end_AUt_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, add, W, T))
            end
        end
    else
        nindexpr = :($nr + nkern * $kernel_size_n)
        q = quote
            Base.Cartesian.@nexprs $nr n -> begin
                triangle_ind_n = triangle_ind_tuple[n]
                Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[m+$mindexpr,n]
            end
            $(mkernel_end_AUt_quote(kernel_size_m, nr, mindexpr, 0, add, W, T))
            for nkern in 0:$(nk-1)
                Base.Cartesian.@nexprs $kernel_size_n n -> begin
                    triangle_ind_n = triangle_ind_tuple[n + $nindexpr]
                    Base.Cartesian.@nexprs $kernel_size_m m -> A_m_n = A[m+$mindexpr,n+$nindexpr]
                end
                $(square_iterations_AUt_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, add, W, T))
                $(mkernel_end_AUt_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr, add, W, T))
            end
        end
    end
    q
end

function matrix_adjoint_upper_triangle_mul_quote(M::Int, N::Int, T::DataType, add::Bool = true, W::Int = 0)
    mk, mr, nk, nr, kernel_size_m, kernel_size_n = PaddedMatrices.determine_pattern(M,N)

    if nr > 0
        # @assert nk > 0
        # nk -= 1
        # nr = kernel_size_n
        if (N % (nk+1)) == 0
            nr = 0
            nk += 1
            kernel_size_n = N ÷ nk
        end
    end
    if mr > 0
        if (M % (mk+1)) == 0
            mr = 0
            mk += 1
            kernel_size_m = M ÷ mk
        end
    end
    q_triangle_ind_defs = Expr(:tuple, Int32(N))
    triangle_ind = N
    triangle_increment = 0
    for n in 2:N
        push!(q_triangle_ind_defs.args, Int32(triangle_ind))
        triangle_increment += 1
        triangle_ind += triangle_increment
    end

    q = quote
        triangle_ind_tuple = $q_triangle_ind_defs
        for mkern in 0:$(mk-1)
            $(mkern_iteration_AUt_quote(kernel_size_m, kernel_size_n, nk, nr, :(mkern * $kernel_size_m), add, W, T))
        end
    end
    if mr > 0
        push!(q.args, mkern_iteration_AUt_quote(mr, kernel_size_n, nk, nr, mk * kernel_size_m, add, W, T))
    end
    q
end

@generated function Base.:*(
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            Ut::LinearAlgebra.Adjoint{Union{},<: AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        AUt = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        # AU = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        @inbounds for i in 1:$(M*N)
            AUt[i] = vbroadcast(Vec{$W,$T}, zero($T))
        end
        @inbounds begin
            $(matrix_adjoint_upper_triangle_mul_quote(M, N, T, true, 0))
        end
        AUt
    end
end
@generated function Base.:*(
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            Ut::LinearAlgebra.Adjoint{T,<: AbstractUpperTriangularMatrix{N,T,L}}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        AUt = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        # AU = FixedSizeMatrix{$M,$N,NTuple{$W,Core.VecElement{$T}}}(undef)
        @inbounds for i in 1:$(M*N)
            AUt[i] = vbroadcast(Vec{$W,$T}, zero($T))
        end
        @inbounds begin
            $(matrix_adjoint_upper_triangle_mul_quote(M, N, T, true, W))
        end
        AUt
    end
end

@generated function addmul!(
            AUt::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            Ut::LinearAlgebra.Adjoint{Union{},<: AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_adjoint_upper_triangle_mul_quote(M, N, T, true, 0))
        end
        AUt
    end
end
@generated function addmul!(
            AUt::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            Ut::LinearAlgebra.Adjoint{T,<: AbstractUpperTriangularMatrix{N,T,L}}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_adjoint_upper_triangle_mul_quote(M, N, T, true, W))
        end
        AUt
    end
end

@generated function submul!(
            AUt::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            Ut::LinearAlgebra.Adjoint{Union{},<: AbstractUpperTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}}
        ) where {M,N,W,T,L}
    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_adjoint_upper_triangle_mul_quote(M, N, T, false, 0))
        end
        AUt
    end
end
@generated function submul!(
            AUt::PaddedMatrices.AbstractMutableFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            A::AbstractFixedSizeMatrix{M,N,NTuple{W,Core.VecElement{T}}},
            Ut::LinearAlgebra.Adjoint{T,<: AbstractUpperTriangularMatrix{N,T,L}}
        ) where {M,N,W,L,T}
    quote
        $(Expr(:meta,:inline))
        @inbounds begin
            $(matrix_adjoint_upper_triangle_mul_quote(M, N, T, false, W))
        end
        AUt
    end
end

function mkernal_start_triangle_view_quote(kernel_size_m::Int, kernel_size_n::Int, mindexpr, nindexpr, add::Bool)#, W::Int, T::DataType)
    muladdfunc = add ? :vmuladd : :vfnmadd

    load_As = quote
        $([:($(Symbol(:A_,m,:_,n)) = A[$mindexpr + $m, $nindexpr + $n]) for m in 1:kernel_size_m, n in 1:kernel_size_n]...)
    end

    if kernel_size_m == 1
        mkern_start_nloop = quote
            U_1_1 = SIMDPirates.$muladdfunc(A_1_n, B[n+$nindexpr,m_col], U_1_1)
        end

        return quote
            $load_As
            tri_ind = tri_ind_base
            tri_increment = tri_increment_base
            m_col = 1+$mindexpr
            U_1_1 = U[m_col]
            Base.Cartesian.@nexprs $kernel_size_n n -> begin
                $mkern_start_nloop
            end
            U[m_col]   = U_1_1
            tri_ind += tri_increment
            tri_ind_col = tri_ind - $kernel_size_m
            tri_increment_col = tri_increment
            tri_increment += 1
        end

    end

    mkern_start_nloop = quote
        vB = B[n+$nindexpr,m_col+1]
        U_1_1 = SIMDPirates.$muladdfunc(A_1_n, B[n+$nindexpr,m_col], U_1_1)
        U_1_2 = SIMDPirates.$muladdfunc(A_1_n, vB, U_1_2)
        U_2_2 = SIMDPirates.$muladdfunc(A_2_n, vB, U_2_2)
    end

    mkern_start = quote
        $load_As
        tri_ind = tri_ind_base
        tri_increment = tri_increment_base
        m_col = 1+$mindexpr
        U_1_1 = U[m_col]
        U_1_2 = U[tri_ind]
        U_2_2 = U[m_col+1]
        Base.Cartesian.@nexprs $kernel_size_n n -> begin
            $mkern_start_nloop
        end
        U[m_col]   = U_1_1
        U[tri_ind] = U_1_2
        tri_ind += tri_increment
        U[m_col+1] = U_2_2
    end
    if kernel_size_m == 2
        push!(mkern_start.args, quote
            tri_ind_col = tri_ind - $kernel_size_m
            tri_increment_col = tri_increment
            tri_increment += 1
        end)
        return mkern_start
    else
        push!(mkern_start.args, :(tri_increment += 1))
    end

    push!(mkern_start.args, quote
    # println("tri_ind = $tri_ind; column: $(m_col + 2) ")
        U_1_3 = U[tri_ind-1]
        U_2_3 = U[tri_ind]
        U_3_3 = U[m_col+2]
    end)
    for n in 1:kernel_size_n
        push!(mkern_start.args, quote
            vB = B[$n+$nindexpr,m_col+2]
            U_1_3 = SIMDPirates.$muladdfunc($(Symbol(:A_1_,n)), vB, U_1_3)
            U_2_3 = SIMDPirates.$muladdfunc($(Symbol(:A_2_,n)), vB, U_2_3)
            U_3_3 = SIMDPirates.$muladdfunc($(Symbol(:A_3_,n)), vB, U_3_3)
        end)
    end
    push!(mkern_start.args, quote
        U[tri_ind-1] = U_1_3
        U[tri_ind] = U_2_3
        tri_ind += tri_increment
        U[m_col+2] = U_3_3
    end)
    if kernel_size_m == 4
        push!(mkern_start.args, quote
            tri_increment += 1
            # println("tri_ind = $tri_ind; column: $(m_col +3) ")
            U_1_4 = U[tri_ind-2]
            U_2_4 = U[tri_ind-1]
            U_3_4 = U[tri_ind]
            U_4_4 = U[m_col+3]
        end)
        for n in 1:kernel_size_n
            push!(mkern_start.args, quote
                vB = B[$n+$nindexpr,m_col+3]
                U_1_4 = SIMDPirates.$muladdfunc($(Symbol(:A_1_,n)), vB, U_1_4)
                U_2_4 = SIMDPirates.$muladdfunc($(Symbol(:A_2_,n)), vB, U_2_4)
                U_3_4 = SIMDPirates.$muladdfunc($(Symbol(:A_3_,n)), vB, U_3_4)
                U_4_4 = SIMDPirates.$muladdfunc($(Symbol(:A_4_,n)), vB, U_4_4)
            end)
        end
        push!(mkern_start.args, quote
            U[tri_ind-2] = U_1_4
            U[tri_ind-1] = U_2_4
            U[tri_ind] = U_3_4
            tri_ind += tri_increment
            # tri_increment += 1
            tri_increment_col = tri_increment
            tri_increment += 1
            U[m_col+3] = U_4_4
            tri_ind_col = tri_ind - $kernel_size_m
        end)
    else
        push!(mkern_start.args, quote
            tri_ind_col = tri_ind - $kernel_size_m
            # tri_increment += 1
            tri_increment_col = tri_increment
            tri_increment += 1
        end)
    end
    mkern_start
end

function square_iterations_triangle_view_quote(kernel_size_m::Int, kernel_size_n::Int, nindexpr, M::Int, add::Bool)
    muladdfunc = add ? :vmuladd : :vfnmadd
    ## The operation is UpperTriangle(MxN * NxM) matrices A and B
    ## this function iterates across the columns of the upper triangular product
    ## The idea of this function is to iterate from
    quote
        for m_coliter in m_col+$kernel_size_m:$M
            Base.Cartesian.@nexprs $kernel_size_m m -> begin
                # println("tri_ind_col + m = $tri_ind_col + $m; column: $m_coliter")
                U_m = U[tri_ind_col + m]
            end
            Base.Cartesian.@nexprs $kernel_size_n n -> begin
            # for n ∈ 1:$kernel_size_n
                vB = B[n+$nindexpr,m_coliter]
                Base.Cartesian.@nexprs $kernel_size_m m -> begin
                    U_m = SIMDPirates.$muladdfunc(A_m_n, vB, U_m)
                end
            end
            Base.Cartesian.@nexprs $kernel_size_m m -> U[tri_ind_col + m] = U_m
            tri_ind_col += tri_increment_col
            tri_increment_col += 1
        end
    end
end

function mul_mat_of_vecs_upper_triangle_view_quote(M::Int,N::Int,W::Int,T::DataType,add::Bool = true)

    mk, mr, nk, nr, kernel_size_m, kernel_size_n = PaddedMatrices.determine_pattern(M,N)
    P = M

    if mr > 0
        if (M % (mk+1)) == 0
            mr = 0
            mk += 1
            kernel_size_m = div(M, mk)
        end
    end
    if nk == 1 && nr > 0
        # More evenly split the kernel.
        kernel_size_n += nr
        nr = (kernel_size_n >> 1)
        kernel_size_n -= nr
    end
    mindexpr = :(mkern * $kernel_size_m)
    nindexpr = :(nkern * $kernel_size_n)

    mindexpr_complete = mk * kernel_size_m
    nindexpr_complete = nk * kernel_size_n
    # At start, A block is already loaded. Now calc U upper.

    mkern_loop_body = quote
        for nkern in 0:$(nk-1)
            $(mkernal_start_triangle_view_quote(kernel_size_m, kernel_size_n, mindexpr, nindexpr,add))
            $(square_iterations_triangle_view_quote(kernel_size_m, kernel_size_n, nindexpr,M, add))
        end
    end
    if nr > 0
        push!(mkern_loop_body.args, mkernal_start_triangle_view_quote(kernel_size_m, nr, mindexpr, nindexpr_complete,add))
        push!(mkern_loop_body.args, square_iterations_triangle_view_quote(kernel_size_m, nr, nindexpr_complete,M,add))
    end

    q = quote
        tri_ind = 0
        tri_increment = 0
        tri_ind_base = $(M + 1)
        tri_increment_base = 2
        @inbounds for mkern in 0:$(mk-1)
            $mkern_loop_body

            tri_ind_base = tri_ind + tri_increment
            tri_increment_base = tri_increment + 1
        end
    end
    if mr > 0
        push!(q.args,
            quote
                @inbounds begin
                    for nkern in 0:$(nk-1)
                        $(mkernal_start_triangle_view_quote(mr, kernel_size_n, mindexpr_complete, nindexpr,add))
                    end
                    $(nr == 0 ? nothing : mkernal_start_triangle_view_quote(mr, kernel_size_n, mindexpr_complete, nindexpr_complete,add))
                end
            end
        )
    end
    push!(q.args, :U)
    q
end

# @generated function mul_upper_triangle_view!(
@generated function addmul!(
            U::MutableUpperTriangularMatrix{M,NTuple{W,Core.VecElement{T}}},
            A::LinearAlgebra.Adjoint{Union{},<: AbstractFixedSizeMatrix{N,M,NTuple{W,Core.VecElement{T}}}},
            B::AbstractFixedSizeMatrix{N,M,NTuple{W,Core.VecElement{T}}}
        ) where {M,N,W,T}

    quote
        $(Expr(:meta,:inline))
        $(mul_mat_of_vecs_upper_triangle_view_quote(M,N,W,T,true))
    end
end
@generated function submul!(
            U::MutableUpperTriangularMatrix{M,NTuple{W,Core.VecElement{T}}},
            A::LinearAlgebra.Adjoint{Union{},<: AbstractFixedSizeMatrix{N,M,NTuple{W,Core.VecElement{T}}}},
            B::AbstractFixedSizeMatrix{N,M,NTuple{W,Core.VecElement{T}}}
        ) where {M,N,T,W}
        # ) where {M,N,W,T}

    quote
        $(Expr(:meta,:inline))
        $(mul_mat_of_vecs_upper_triangle_view_quote(M,N,W,T,false))
    end
end

@generated function mul_upper_triangle_view(
            A::LinearAlgebra.Adjoint{Union{},<: AbstractFixedSizeMatrix{N,M,NTuple{W,Core.VecElement{T}}}},
            B::AbstractFixedSizeMatrix{N,M,NTuple{W,Core.VecElement{T}}}
        ) where {M,N,T,W}
        # ) where {M,N,W,T}

    L = binomial2(M+1)
    quote
        $(Expr(:meta,:inline))
        U = MutableUpperTriangularMatrix{$M,NTuple{$W,Core.VecElement{$T}},$L}(undef)
        vzero = vbroadcast(NTuple{$W,Core.VecElement{$T}}, zero($T))
        @inbounds for l in 1:$L
            U[l] = vzero
        end
        $(mul_mat_of_vecs_upper_triangle_view_quote(M,N,W,T,true))
    end
end

function all_finite(U::AbstractDiagTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}) where {N,W,T,L}
    @inbounds begin
        reduction = SIMDPirates.visfinite(U[1])
        for l ∈ 2:L
            reduction = SIMDPirates.vand(reduction, SIMDPirates.visfinite(U[l]))
        end
    end
    SIMDPirates.vall(reduction)
end

function PaddedMatrices.zero!(A::AbstractDiagTriangularMatrix{N,T,L}) where {N,T,L}
    @inbounds for i ∈ 1:L
        A[i] = zero(T)
    end
end
# function PaddedMatrices.zero!(A::AbstractDiagTriangularMatrix{N,NTuple{W,Core.VecElement{T}},L}) where {N,W,T,L}
#     @inbounds for i ∈ 1:L
#         A[i] = SIMDPirates.vbroadcast(NTuple{W,Core.VecElement{T}}, zero(T))
#     end
# end