text:
### A Pluto.jl notebook ###
# v0.17.7

using Markdown
using InteractiveUtils

# ╔═╡ 64bac944-8008-498d-a89b-b9f8ee54aa98
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using PlutoUI
    PlutoUI.TableOfContents()
end

# ╔═╡ a022be72-b0b0-43b7-8f39-c02a53875ff6
using OpenML, MLJ, MLJXGBoostInterface, DataFrames, MLJLinearModels, MLJDecisionTreeInterface

# ╔═╡ cf267836-0ac5-494d-9b5e-12eb314ac170
using MLJLIBSVMInterface, Random, Plots

# ╔═╡ 01a467a5-7389-44d1-984d-244dfb1ea39f
begin
    using MLCourse
    MLCourse.list_notebooks(@__FILE__)
end

# ╔═╡ 7cfd274d-75d8-4cd9-b38b-4176466c9e26
Markdown.parse("# Fitting MNIST with a Convolutional Neural Networks

In the following we fit a convolutional neural network to the MNIST data set.

Running this example takes more than 10 minutes.

We load the data
```julia
using MLJ, Plots, MLJFlux, Flux, OpenML, DataFrames
mnist_x, mnist_y = let df = OpenML.load(554) |> DataFrame
    Float32.(df[:, 1:end-1] ./ 255), # scaling the input
    df.class
end
```
and transform the input to an image representation
```julia
images = coerce(PermutedDimsArray(reshape(Array(mnist_x), :, 28, 28), (3, 2, 1)),
                GrayImage);
plot(images[1])
```
$(MLCourse.embed_figure("mnist_example1.png"))

Now we define our network builder
```julia
builder = MLJFlux.@builder(
              begin
                  front = Chain(Conv((8, 8), n_channels => 16, relu),
                                Conv((4, 4), 16 => 32, relu),
                                Flux.flatten)
                  d = first(Flux.outputsize(front, (n_in..., n_channels, 1)))
                  Chain(front, Dense(d, n_out))
              end)
```
The `front` consists of two convolutional layers. The volume of the second
convolutional layer is flattened and taken as input to a single dense layer.
```julia
m = machine(ImageClassifier(builder = builder,
                            batch_size = 32,
                            epochs = 5),
            images, mnist_y)
fit!(m, rows = 1:60000, verbosity = 2)
```
Finally, let us compute the test error rate
```julia
mean(predict_mode(m, images[60001:70000]) .!= mnist_y[60001:70000])
```
The result is a misclassification rate of 1.4% which is again a strong improvement
over the MLP result. Note that we achieved this value without tuning of the hyper-parameters. Well-tuned convolutional networks achieve performances below 1% error.

If you are interested in learning more about convolutional neural networks I can highly recommend to have a look at one of the excellent tutorials on the internet, like [this](https://towardsdatascience.com/intuitively-understanding-convolutions-for-deep-learning-1f6f42faee1), [this](https://towardsdatascience.com/a-comprehensive-guide-to-convolutional-neural-networks-the-eli5-way-3bd2b1164a53) or [this](https://distill.pub/2017/feature-visualization/) one.
")

# ╔═╡ dd12ea16-f087-4507-a571-c95bf38aa976
Markdown.parse("# Language Detection with a Recurrent Neural Network

The following code (adapted from [here](https://github.com/FluxML/model-zoo/blob/master/text/lang-detection/model.jl)) uses wikipedia articles to detect the language of a sentence
by simply reading each character of the sentence with a recurrent neural network
and applying a classifier to the last hidden state.

Here are some example sentences in English, French, Italian, Spanish and Danish
after preprocessing by the `get_processed_data` function below. The preprocessing replaces most special characters by a single special character (we use the special character ˌ below). Ultimately, the recurrent neural network gets these sentences in one-hot encoded form; for better readability for humans we print them here in the standard form:

\"jewish philosophy refers to the conjunction between serious study of philosophy and jewish theology\"
 \"during the early years of the second templeˌ the highest religious authority was a council known as the great assemblyˌ led by ezra of the book of ezra\"
 \"en 756ˌ abderraman i ˌun omeya superviviente del exterminio de la familia califal destronada por los abbasiesˌ fue acogido por sus partidarios en alˌandalus y se impuso como emir\"
 \"they also saw an elite population convert to judaism ˌthe khazarsˌˌ only to disappear as the centers of power in the lands once occupied by that elite fell to the people of rus and then the mongols\"
 \"ˌ167ˌ\\ncritics argue that wikipediaˌs open nature and a lack of proper sources for most of the information makes it unreliable\"
 \"tambien existia una clase media formada por comerciantesˌ pequenos propietariosˌ funcionarios y profesionales liberalesˌ que en general toleraba el franquismo pero no comulgaba totalmente con sus ideas\"
 \"wikipedia ˌpronunciaˌ vedi sottoˌ e unˌenciclopediaonline a contenuto liberoˌ collaborativaˌ multilingue e gratuitaˌ nata nel 2001ˌ sostenuta e ospitata dalla wikimedia foundationˌ unˌorganizzazione non a scopo di lucrostatunitense\"
 \"dˌautre partˌ des jalons sont poses pour une integration de toutes les economies nationales dans un grand espace economique domine par lˌallemagneˌ65ˌ\"
 \"ˌ423ˌˌ el proyecto preveia grandes infraestructuras viarias ˌbulevaresˌ grandes plazasˌ paseos de rondaˌ diagonalesˌ paseos maritimosˌˌ parques y jardinesˌ enlaces ferroviarios ˌcon las lineas interiores soterradasˌˌ edificios publicos y colectivos en los puntos centrales de las trazas viariasˌ equipamientos y areas de servicios\"
 \"ˌ24ˌ\\nif the stoppage is earlier in the gameˌ or if it is a playoff or grey cup gameˌ play may be stopped for up to 3 hours and then resume\"
 \"ˌ195ˌˌ el estudio se realizo comparando 42 articulos de ambas obras por un comite de expertos sin que estos supieran de cual de las dos enciclopedias provenian\"
 \"las evidencias arqueologicas de establecimientos fenicios ˌebusus ˌibizaˌˌ sexi ˌalmunecarˌˌ malaka ˌmalagaˌˌ permiten hablar de un monopolio fenicio de las rutas comerciales en torno al estrecho de gibraltar ˌincluyendo las del atlanticoˌ como la ruta del estanoˌˌ que limito la colonizacion griega al norte mediterraneo ˌemporionˌ la actual ampuriasˌ\"
 \"ˌ158ˌˌ en cuanto al sistema monetarioˌ al dinero carolingio sucedio tras la independizacion del condado el dinero de vellonˌ con un 75ˌˌ menos de plata\"
 \"es responsable ante las cortes generales\"
 \"ˌ818ˌˌ\\n\\nlas elecciones del 21 de diciembre fueron ganadas en escanos y votos por ciudadanosˌ cuya candidata era ines arrimadas\"
 \"ˌ102ˌ\\ndopo pochi giorni il marito mori\"
 \"ˌ116ˌ\\na kittel ˌyiddishˌ ˌˌˌˌˌˌ a white kneeˌlength overgarmentˌ is worn by prayer leaders and some observant traditional jews on the high holidays\"
 \"esta es una lista de las grafias alternativas usadas para nombrar a wikipedia en las diferentes edicionesˌˌ53ˌˌ\\n\\nnotaˌ cada asterisco ˌˌˌ denota una referencia respectiva por la multiple adopcion linguistica\"
 \"afin dˌassurer la perennite des liens vers ces sourcesˌ internet archive a repare en 2018 pres de 9 millions de liens brisesˌ37ˌ\"
 \"de hechoˌ 18 de los 25 mayores hoteles del mundo se encuentran en ella\"

The code in this section is somewhat more advanced than in the rest of this course.
You do not need to understand everything. But if you are interested in learning
about recurrent neural networks, this may give you an idea for how to implement
them.

WARNING: This code takes more than 30 minutes to run.

```julia
using Flux
using Flux: onehot, onehotbatch, logitcrossentropy, reset!, throttle
using Statistics: mean
using Random
using Unicode

Base.@kwdef mutable struct Args
    lr::Float64 = 1e-3     # learning rate
    N::Int = 32            # Number of perceptrons in hidden layer
    test_len::Int = 100    # length of test data
    langs_len::Int = 0     # Number of different languages in Corpora
    alphabet_len::Int = 0  # Total number of characters possible, in corpora
    throttle::Int = 120    # throttle timeout
end

function get_processed_data(args)
    corpora = Dict()

    dir = joinpath(@__DIR__, \"..\", \"data\", \"corpus\")
    for file in readdir(dir)
        lang = Symbol(match(r\"(.*)\\.txt\", file).captures[1])
        corpus = split(String(read(joinpath(dir, file))), \".\")
        corpus = strip.(Unicode.normalize.(corpus, casefold=true, stripmark=true))
        corpus = filter(!isempty, corpus)
        corpora[lang] = corpus
    end

    langs = collect(keys(corpora))
    args.langs_len = length(langs)
    alphabet = ['a':'z'; '0':'9'; ' '; '\\n'; '_']
    args.alphabet_len = length(alphabet)

    # See which chars will be represented as \"unknown\"
    unique(filter(x -> x ∉ alphabet, join(vcat(values(corpora)...))))

    dataset = [(onehotbatch(s, alphabet, '_'), onehot(l, langs))
               for l in langs for s in corpora[l]] |> shuffle

    train, test = dataset[1:end-args.test_len], dataset[end-args.test_len+1:end]
    return train, test
end

function build_model(args)
    # The scanner is the recurrent neural network; it uses a special kind of
    # neuron in the recurrent layer, called LSTM.
    scanner = Chain(Dense(args.alphabet_len, args.N, relu), LSTM(args.N, args.N))
    # The encoder takes the last element of the sequence and computes the
    # log-probability of each language for the given sequence using a standard
    # multi-layer perceptron (MLP) with one hidden layer.
    encoder = Chain(Dense(args.N, args.N, relu), Dense(args.N, args.langs_len))
    return scanner, encoder
end

function model(x, scanner, encoder)
    # scan all characters of sentence x
    state = mapslices(scanner, x, dims = 1)[:, end]
    reset!(scanner)
    # take the last state of the neural network
    # to compute the probability of each language.
    encoder(state)
end

args = Args()
train_data, test_data = get_processed_data(args)

@info(\"Constructing Model...\")
scanner, encoder = build_model(args)

# we use the negative loglikelihood loss (here called logitcrossentropy)
loss(x, y) = logitcrossentropy(model(x, scanner, encoder), y)
function evaluate(scanner, encoder, dataset)
    () -> begin
        tmp = [(model(x, scanner, encoder), y) for (x, y) in dataset]
        (loss = mean(logitcrossentropy(ŷ, y) for (ŷ, y) in tmp),
         accuracy = mean(argmax(ŷ) == findfirst(y) for (ŷ, y) in tmp))
    end
end
traineval = evaluate(scanner, encoder, train_data)
testeval = evaluate(scanner, encoder, test_data)

opt = ADAM(args.lr)
ps = params(scanner, encoder)
evalcb = throttle(() -> @show(traineval(), testeval()), args.throttle)
@info(\"Training...\")
for epoch in 1:100
    @show epoch
    Flux.train!(loss, ps, train_data, opt, cb = evalcb)
end
```
The result is around 60% accuracy on the test set (which is much better than the
chance level of 20%). With longer training, more neurons in the scanner and the
encoder and potentially well chosen regularization, one could improve this
result a lot.
")

# ╔═╡ e68603b3-6f95-4c55-aded-14eaa7615176
md"# Tree-Based Methods

In this section we will analyse the cars data set and try to predict the miles per gallon a car drives based on different predictors. We will use a one-hot encoding for all categorical predictors.

You should observe the following:
1. It is fast to tune a gradient boosting machine in the most important hyper-parameters learning rate `:eta`, number of rounds `:num_round` and maximal tree depth `:max_depth`.
2. A tuned gradient boosting machine is clearly better than a linear regressor.
3. Even without any tuning, a random forest regressor with 500 trees reaches almost the same performance as the the tuned gradient boosting machine.
4. It takes time to tune a neural network regressor that works well.
"

# ╔═╡ 987264db-aca7-41c8-aa3a-b9437fcbd7cd
cars = dropmissing(DataFrame(OpenML.load(455)));

# ╔═╡ 5be3a1e7-de72-49e6-8170-fc5aac38e3bb
cars_input = MLJ.transform(machine(OneHotEncoder(drop_last = true),
	                               select(cars, Not([:name, :mpg]))) |> fit!,
	                       select(cars, Not([:name, :mpg])))

# ╔═╡ ef75fad9-41ca-4960-8c86-3f71c9581582
m1 = machine(LinearRegressor(), cars_input, cars.mpg) |> fit!;

# ╔═╡ 89a6c574-aa3d-4eec-8bf3-50c0d64bee40
evaluate!(m1, measure = rmse)

# ╔═╡ ba5563e7-844e-46b9-a4fd-c2d277d3ed5c
begin
	xgb = XGBoostRegressor()
    m2 = machine(TunedModel(model = xgb,
                            resampling = CV(nfolds = 4),
                            tuning = Grid(goal = 25),
                            range = [range(xgb, :eta,
                                           lower = 1e-2, upper = .1, scale = :log),
                                     range(xgb, :num_round, lower = 50, upper = 500),
                                     range(xgb, :max_depth, lower = 2, upper = 6)]),
                  cars_input, cars.mpg) |> fit!
end;

# ╔═╡ 56b1d5da-eb6d-4a80-92ac-1de73dc76e9f
evaluate!(machine(report(m2).best_model, cars_input, cars.mpg), measure = rmse)

# ╔═╡ df56c35c-713a-494c-b59f-3f54644bac71
m3 = machine(RandomForestRegressor(n_trees = 500),
	         cars_input, cars.mpg);

# ╔═╡ 2942f37c-b848-49c9-a599-0deb5040045c
evaluate!(m3, measure = rmse)

# ╔═╡ c4b21c4a-f8ec-468e-b86b-7e4c3e0203df
md"The following code takes a moment to run. It is certainly not the best neural network regressor, but it would take quite a bit of time to tune it and find one that outperforms gradient boosting.
```julia
using MLJFlux
m4 = machine(NeuralNetworkRegressor(builder = MLJFlux.Short(n_hidden = 64,
	                                                        dropout = 0.1,
	                                                        σ = MLJFlux.Flux.relu),
                                    epochs = 5e3, batch_size = 20),
	         cars_input, cars.mpg);
evaluate!(m4, measure = rmse)
```
The result is a cross-validation estimate of the test error of approximately 3.13,
which is better than linear regression, but worse than the results we found quickly with the tree based approaches.
"

# ╔═╡ 9fd69edd-9d2d-4341-ac34-efe9b4c86f68
md"# Support Vector Machines

In this section we fit a support vector machine to some binary classification task, where the data is rather well, but non-linearly, separated. It is charateristic of support vector machines to find a flexible decision boundary that matches well the decision boundary that humans would intuitively draw in such situations.
"

# ╔═╡ 2d9cf0a4-6ba7-438b-8e8f-4a911a340f06
Random.seed!(23); X, y = make_blobs(centers = 3); # generate some random data

# ╔═╡ 87a48194-fc87-414f-8d44-8bf00a87b3ec
y[y .== 3] .= 1; # we want only two classes here

# ╔═╡ 3a7fa61a-69e0-4082-9bb2-13df40089bb2
m5 = machine(SVC(), X, y) |> fit!;

# ╔═╡ a48c53e8-99f7-4b54-9a80-0be243b4cddb
let grid = MLCourse.grid(-5:.5:11.5, -9:.5:12, names = (:x1, :x2))
	pred = predict(m5, grid)
	scatter(grid.x1, grid.x2, c = Int.(int.(pred)),
		    markersize = 2, markerstrokewidth = 0, legend = false)
	scatter!(X.x1, X.x2, c = Int.(int.(y)), xlabel = "x₁", ylabel = "x₂")
end

# ╔═╡ 2bbf9876-0676-11ec-3985-73f4dcaea02f
Markdown.parse("# Exercises
## Conceptual

1. Here below is an image (with padding 1 already applied). We would like to process it with a convolutional network with one convolution layer with two ``3 \\times 3`` filters (depicted below the image), stride 1 and relu non-linearity.
    - Determine the width, height and depth of the volume after the convolutional layer.
    - Compute the output of the convolutional layer assuming the two biases to be zero.
$(MLCourse.embed_figure("conv_exercise.png"))
2. Given a volume of width ``n``, height ``n`` and depth ``c`` and a convolutional layer with ``k`` filters of size ``f\\times f\\times c`` with stride ``s`` and padding ``p``.
    - Convince yourself that this convolution layer has an output volume of size ``\\left( \\frac{n + 2p - f}{s} + 1\\right)\\times \\left( \\frac{n + 2p - f}{s} + 1\\right)\\times k\\, .``
    - Flux knows the padding option `pad = SamePad()`. It means that the output volume should have the same x and y dimension as the input volume. In this exercise we compute the padding needed this option.  What padding ``p`` do you have to choose for a given ``n`` and ``f`` such that the input and output volumes have the same width and depth for stride ``s=1``. Check you result for the special case of ``n=4`` and ``f=3``.
    - Compute the number of weights and biases of this convolutional layer.
    - Compute the number of weights and biases for a dense layer (standard MLP layer) that has the same number of activations as this convolutional layer. Is the number of parameters larger for the dense layer or the convolutional layer?

## Applied
1. In this exercise our goal is to find a good machine learning model to classify images of Zalando's articles. You can load a description of the so-called Fashion-MNIST data set with `OpenML.describe_dataset(40996)` and load the data set with `OpenML.load(40996)`. Take our recipe for supervised learning (last slide of the presentation on \"Model Assessment and Hyperparameter Tuning\") as a guideline. Hints: cleaning is not necessary, but plotting some examples is advisable; linear classification is a good starting point for a first benchmark, but you should also explore other models like random forests (`RandomForestClassifier`), multilayer perceptrons (`NeuralNetworkClassifier`) and convolutional neural networks (`ImageClassifer`) and play manually a bit with their hyper-parameters (proper tuning with `TunedModel` may be too time-intensive). To reduce the computation time, we will only use the first 5000 images of the full dataset for training and we will not do cross-validation here but instead use samples 5001 to 10000 as a validation set to estimate the classification accuracy. *Hint:* you can use the resampling strategy `Holdout` for model tuning and evaluation.
")

# ╔═╡ 3bd0ed27-7d3b-4036-87a4-2dd0e2874b58
MLCourse.footer()

# ╔═╡ Cell order:
# ╟─7cfd274d-75d8-4cd9-b38b-4176466c9e26
# ╟─dd12ea16-f087-4507-a571-c95bf38aa976
# ╟─e68603b3-6f95-4c55-aded-14eaa7615176
# ╠═a022be72-b0b0-43b7-8f39-c02a53875ff6
# ╠═987264db-aca7-41c8-aa3a-b9437fcbd7cd
# ╠═5be3a1e7-de72-49e6-8170-fc5aac38e3bb
# ╠═ef75fad9-41ca-4960-8c86-3f71c9581582
# ╠═89a6c574-aa3d-4eec-8bf3-50c0d64bee40
# ╠═ba5563e7-844e-46b9-a4fd-c2d277d3ed5c
# ╠═56b1d5da-eb6d-4a80-92ac-1de73dc76e9f
# ╠═df56c35c-713a-494c-b59f-3f54644bac71
# ╠═2942f37c-b848-49c9-a599-0deb5040045c
# ╟─c4b21c4a-f8ec-468e-b86b-7e4c3e0203df
# ╟─9fd69edd-9d2d-4341-ac34-efe9b4c86f68
# ╠═cf267836-0ac5-494d-9b5e-12eb314ac170
# ╠═2d9cf0a4-6ba7-438b-8e8f-4a911a340f06
# ╠═87a48194-fc87-414f-8d44-8bf00a87b3ec
# ╠═3a7fa61a-69e0-4082-9bb2-13df40089bb2
# ╠═a48c53e8-99f7-4b54-9a80-0be243b4cddb
# ╟─2bbf9876-0676-11ec-3985-73f4dcaea02f
# ╟─01a467a5-7389-44d1-984d-244dfb1ea39f
# ╟─64bac944-8008-498d-a89b-b9f8ee54aa98
# ╟─3bd0ed27-7d3b-4036-87a4-2dd0e2874b58

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
module day18

function evalnext(s; addfirst=false)
    while true
        m = match(r"\((\d+[\+\*])+\d+\)", s)
        if isnothing(m) break end
        val = evalnext(m.match[2:end-1]; addfirst)
        s = replace(s, m.match => "$val")
    end

    if addfirst
        while true
            m = match(r"\d+\+\d+", s)
            if isnothing(m) break end
            val = +(parse.(Int64, split(m.match, '+'))...)
            s = replace(s, m.match => "$val")
        end
    end

    lhs = match(r"\d+", s).match
    current = parse(Int64, lhs)
    if length(lhs) == length(s) return current end
    s = s[length(lhs)+1:end]
    while true
        op = Dict('+' => +, '*' => *)[s[1]]
        rhs = match(r"\d+", s[2:end]).match
        current = op(current, parse(Int64, rhs))
        if length(rhs) + 1 == length(s) return current end
        s = s[length(rhs)+2:end]
    end
end

equations = map(x -> replace(x, " "=>""), readlines("day18_input.txt"))
println("Part 1: sum = $(sum(evalnext.(equations)))")

# Part 2
println("Part 2: sum = $(sum(evalnext.(equations; addfirst=true)))")

end

ext: jl

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
#==============================================================================#
# Athena.jl
#
# This file is generated from:
# https://github.com/aws/aws-sdk-js/blob/master/apis/athena-2017-05-18.normal.json
#==============================================================================#

__precompile__()

module Athena

using AWSCore

"""
    using AWSSDK.Athena.batch_get_named_query
    batch_get_named_query([::AWSConfig], arguments::Dict)
    batch_get_named_query([::AWSConfig]; NamedQueryIds=)

    using AWSCore.Services.athena
    athena([::AWSConfig], "BatchGetNamedQuery", arguments::Dict)
    athena([::AWSConfig], "BatchGetNamedQuery", NamedQueryIds=)

# BatchGetNamedQuery Operation

Returns the details of a single named query or a list of up to 50 queries, which you provide as an array of query ID strings. Use [ListNamedQueries](@ref) to get the list of named query IDs. If information could not be retrieved for a submitted query ID, information about the query ID submitted is listed under [UnprocessedNamedQueryId](@ref). Named queries are different from executed queries. Use [BatchGetQueryExecution](@ref) to get details about each unique query execution, and [ListQueryExecutions](@ref) to get a list of query execution IDs.

# Arguments

## `NamedQueryIds = [::String, ...]` -- *Required*
An array of query IDs.

# Returns

`BatchGetNamedQueryOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/BatchGetNamedQuery)
"""
@inline batch_get_named_query(aws::AWSConfig=default_aws_config(); args...) = batch_get_named_query(aws, args)

@inline batch_get_named_query(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "BatchGetNamedQuery", args)

@inline batch_get_named_query(args) = batch_get_named_query(default_aws_config(), args)

"""
    using AWSSDK.Athena.batch_get_query_execution
    batch_get_query_execution([::AWSConfig], arguments::Dict)
    batch_get_query_execution([::AWSConfig]; QueryExecutionIds=)

    using AWSCore.Services.athena
    athena([::AWSConfig], "BatchGetQueryExecution", arguments::Dict)
    athena([::AWSConfig], "BatchGetQueryExecution", QueryExecutionIds=)

# BatchGetQueryExecution Operation

Returns the details of a single query execution or a list of up to 50 query executions, which you provide as an array of query execution ID strings. To get a list of query execution IDs, use [ListQueryExecutions](@ref). Query executions are different from named (saved) queries. Use [BatchGetNamedQuery](@ref) to get details about named queries.

# Arguments

## `QueryExecutionIds = [::String, ...]` -- *Required*
An array of query execution IDs.

# Returns

`BatchGetQueryExecutionOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/BatchGetQueryExecution)
"""
@inline batch_get_query_execution(aws::AWSConfig=default_aws_config(); args...) = batch_get_query_execution(aws, args)

@inline batch_get_query_execution(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "BatchGetQueryExecution", args)

@inline batch_get_query_execution(args) = batch_get_query_execution(default_aws_config(), args)

"""
    using AWSSDK.Athena.create_named_query
    create_named_query([::AWSConfig], arguments::Dict)
    create_named_query([::AWSConfig]; Name=, Database=, QueryString=, <keyword arguments>)

    using AWSCore.Services.athena
    athena([::AWSConfig], "CreateNamedQuery", arguments::Dict)
    athena([::AWSConfig], "CreateNamedQuery", Name=, Database=, QueryString=, <keyword arguments>)

# CreateNamedQuery Operation

Creates a named query.

For code samples using the AWS SDK for Java, see [Examples and Code Samples](http://docs.aws.amazon.com/athena/latest/ug/code-samples.html) in the *Amazon Athena User Guide*.

# Arguments

## `Name = ::String` -- *Required*
The plain language name for the query.

## `Description = ::String`
A brief explanation of the query.

## `Database = ::String` -- *Required*
The database to which the query belongs.

## `QueryString = ::String` -- *Required*
The text of the query itself. In other words, all query statements.

## `ClientRequestToken = ::String`
A unique case-sensitive string used to ensure the request to create the query is idempotent (executes only once). If another `CreateNamedQuery` request is received, the same response is returned and another query is not created. If a parameter has changed, for example, the `QueryString`, an error is returned.

**Important**
> This token is listed as not required because AWS SDKs (for example the AWS SDK for Java) auto-generate the token for users. If you are not using the AWS SDK or the AWS CLI, you must provide this token or the action will fail.

# Returns

`CreateNamedQueryOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/CreateNamedQuery)
"""
@inline create_named_query(aws::AWSConfig=default_aws_config(); args...) = create_named_query(aws, args)

@inline create_named_query(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "CreateNamedQuery", args)

@inline create_named_query(args) = create_named_query(default_aws_config(), args)

"""
    using AWSSDK.Athena.delete_named_query
    delete_named_query([::AWSConfig], arguments::Dict)
    delete_named_query([::AWSConfig]; NamedQueryId=)

    using AWSCore.Services.athena
    athena([::AWSConfig], "DeleteNamedQuery", arguments::Dict)
    athena([::AWSConfig], "DeleteNamedQuery", NamedQueryId=)

# DeleteNamedQuery Operation

Deletes a named query.

For code samples using the AWS SDK for Java, see [Examples and Code Samples](http://docs.aws.amazon.com/athena/latest/ug/code-samples.html) in the *Amazon Athena User Guide*.

# Arguments

## `NamedQueryId = ::String` -- *Required*
The unique ID of the query to delete.

# Returns

`DeleteNamedQueryOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/DeleteNamedQuery)
"""
@inline delete_named_query(aws::AWSConfig=default_aws_config(); args...) = delete_named_query(aws, args)

@inline delete_named_query(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "DeleteNamedQuery", args)

@inline delete_named_query(args) = delete_named_query(default_aws_config(), args)

"""
    using AWSSDK.Athena.get_named_query
    get_named_query([::AWSConfig], arguments::Dict)
    get_named_query([::AWSConfig]; NamedQueryId=)

    using AWSCore.Services.athena
    athena([::AWSConfig], "GetNamedQuery", arguments::Dict)
    athena([::AWSConfig], "GetNamedQuery", NamedQueryId=)

# GetNamedQuery Operation

Returns information about a single query.

# Arguments

## `NamedQueryId = ::String` -- *Required*
The unique ID of the query. Use [ListNamedQueries](@ref) to get query IDs.

# Returns

`GetNamedQueryOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/GetNamedQuery)
"""
@inline get_named_query(aws::AWSConfig=default_aws_config(); args...) = get_named_query(aws, args)

@inline get_named_query(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "GetNamedQuery", args)

@inline get_named_query(args) = get_named_query(default_aws_config(), args)

"""
    using AWSSDK.Athena.get_query_execution
    get_query_execution([::AWSConfig], arguments::Dict)
    get_query_execution([::AWSConfig]; QueryExecutionId=)

    using AWSCore.Services.athena
    athena([::AWSConfig], "GetQueryExecution", arguments::Dict)
    athena([::AWSConfig], "GetQueryExecution", QueryExecutionId=)

# GetQueryExecution Operation

Returns information about a single execution of a query. Each time a query executes, information about the query execution is saved with a unique ID.

# Arguments

## `QueryExecutionId = ::String` -- *Required*
The unique ID of the query execution.

# Returns

`GetQueryExecutionOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/GetQueryExecution)
"""
@inline get_query_execution(aws::AWSConfig=default_aws_config(); args...) = get_query_execution(aws, args)

@inline get_query_execution(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "GetQueryExecution", args)

@inline get_query_execution(args) = get_query_execution(default_aws_config(), args)

"""
    using AWSSDK.Athena.get_query_results
    get_query_results([::AWSConfig], arguments::Dict)
    get_query_results([::AWSConfig]; QueryExecutionId=, <keyword arguments>)

    using AWSCore.Services.athena
    athena([::AWSConfig], "GetQueryResults", arguments::Dict)
    athena([::AWSConfig], "GetQueryResults", QueryExecutionId=, <keyword arguments>)

# GetQueryResults Operation

Returns the results of a single query execution specified by `QueryExecutionId`. This request does not execute the query but returns results. Use [StartQueryExecution](@ref) to run a query.

# Arguments

## `QueryExecutionId = ::String` -- *Required*
The unique ID of the query execution.

## `NextToken = ::String`
The token that specifies where to start pagination if a previous request was truncated.

## `MaxResults = ::Int`
The maximum number of results (rows) to return in this request.

# Returns

`GetQueryResultsOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/GetQueryResults)
"""
@inline get_query_results(aws::AWSConfig=default_aws_config(); args...) = get_query_results(aws, args)

@inline get_query_results(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "GetQueryResults", args)

@inline get_query_results(args) = get_query_results(default_aws_config(), args)

"""
    using AWSSDK.Athena.list_named_queries
    list_named_queries([::AWSConfig], arguments::Dict)
    list_named_queries([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.athena
    athena([::AWSConfig], "ListNamedQueries", arguments::Dict)
    athena([::AWSConfig], "ListNamedQueries", <keyword arguments>)

# ListNamedQueries Operation

Provides a list of all available query IDs.

For code samples using the AWS SDK for Java, see [Examples and Code Samples](http://docs.aws.amazon.com/athena/latest/ug/code-samples.html) in the *Amazon Athena User Guide*.

# Arguments

## `NextToken = ::String`
The token that specifies where to start pagination if a previous request was truncated.

## `MaxResults = ::Int`
The maximum number of queries to return in this request.

# Returns

`ListNamedQueriesOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/ListNamedQueries)
"""
@inline list_named_queries(aws::AWSConfig=default_aws_config(); args...) = list_named_queries(aws, args)

@inline list_named_queries(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "ListNamedQueries", args)

@inline list_named_queries(args) = list_named_queries(default_aws_config(), args)

"""
    using AWSSDK.Athena.list_query_executions
    list_query_executions([::AWSConfig], arguments::Dict)
    list_query_executions([::AWSConfig]; <keyword arguments>)

    using AWSCore.Services.athena
    athena([::AWSConfig], "ListQueryExecutions", arguments::Dict)
    athena([::AWSConfig], "ListQueryExecutions", <keyword arguments>)

# ListQueryExecutions Operation

Provides a list of all available query execution IDs.

For code samples using the AWS SDK for Java, see [Examples and Code Samples](http://docs.aws.amazon.com/athena/latest/ug/code-samples.html) in the *Amazon Athena User Guide*.

# Arguments

## `NextToken = ::String`
The token that specifies where to start pagination if a previous request was truncated.

## `MaxResults = ::Int`
The maximum number of query executions to return in this request.

# Returns

`ListQueryExecutionsOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/ListQueryExecutions)
"""
@inline list_query_executions(aws::AWSConfig=default_aws_config(); args...) = list_query_executions(aws, args)

@inline list_query_executions(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "ListQueryExecutions", args)

@inline list_query_executions(args) = list_query_executions(default_aws_config(), args)

"""
    using AWSSDK.Athena.start_query_execution
    start_query_execution([::AWSConfig], arguments::Dict)
    start_query_execution([::AWSConfig]; QueryString=, ResultConfiguration=, <keyword arguments>)

    using AWSCore.Services.athena
    athena([::AWSConfig], "StartQueryExecution", arguments::Dict)
    athena([::AWSConfig], "StartQueryExecution", QueryString=, ResultConfiguration=, <keyword arguments>)

# StartQueryExecution Operation

Runs (executes) the SQL query statements contained in the `Query` string.

For code samples using the AWS SDK for Java, see [Examples and Code Samples](http://docs.aws.amazon.com/athena/latest/ug/code-samples.html) in the *Amazon Athena User Guide*.

# Arguments

## `QueryString = ::String` -- *Required*
The SQL query statements to be executed.

## `ClientRequestToken = ::String`
A unique case-sensitive string used to ensure the request to create the query is idempotent (executes only once). If another `StartQueryExecution` request is received, the same response is returned and another query is not created. If a parameter has changed, for example, the `QueryString`, an error is returned.

**Important**
> This token is listed as not required because AWS SDKs (for example the AWS SDK for Java) auto-generate the token for users. If you are not using the AWS SDK or the AWS CLI, you must provide this token or the action will fail.

## `QueryExecutionContext = ["Database" =>  ::String]`
The database within which the query executes.

## `ResultConfiguration = [ ... ]` -- *Required*
Specifies information about where and how to save the results of the query execution.
```
 ResultConfiguration = [
        "OutputLocation" => <required> ::String,
        "EncryptionConfiguration" =>  [
            "EncryptionOption" => <required> "SSE_S3", "SSE_KMS" or "CSE_KMS",
            "KmsKey" =>  ::String
        ]
    ]
```

# Returns

`StartQueryExecutionOutput`

# Exceptions

`InternalServerException`, `InvalidRequestException` or `TooManyRequestsException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/StartQueryExecution)
"""
@inline start_query_execution(aws::AWSConfig=default_aws_config(); args...) = start_query_execution(aws, args)

@inline start_query_execution(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "StartQueryExecution", args)

@inline start_query_execution(args) = start_query_execution(default_aws_config(), args)

"""
    using AWSSDK.Athena.stop_query_execution
    stop_query_execution([::AWSConfig], arguments::Dict)
    stop_query_execution([::AWSConfig]; QueryExecutionId=)

    using AWSCore.Services.athena
    athena([::AWSConfig], "StopQueryExecution", arguments::Dict)
    athena([::AWSConfig], "StopQueryExecution", QueryExecutionId=)

# StopQueryExecution Operation

Stops a query execution.

For code samples using the AWS SDK for Java, see [Examples and Code Samples](http://docs.aws.amazon.com/athena/latest/ug/code-samples.html) in the *Amazon Athena User Guide*.

# Arguments

## `QueryExecutionId = ::String` -- *Required*
The unique ID of the query execution to stop.

# Returns

`StopQueryExecutionOutput`

# Exceptions

`InternalServerException` or `InvalidRequestException`.

See also: [AWS API Documentation](https://docs.aws.amazon.com/goto/WebAPI/athena-2017-05-18/StopQueryExecution)
"""
@inline stop_query_execution(aws::AWSConfig=default_aws_config(); args...) = stop_query_execution(aws, args)

@inline stop_query_execution(aws::AWSConfig, args) = AWSCore.Services.athena(aws, "StopQueryExecution", args)

@inline stop_query_execution(args) = stop_query_execution(default_aws_config(), args)

end # module Athena

#==============================================================================#
# End of file
#==============================================================================#


text:
# Here are little snippets of code that aren't necessarily used anywhere,
# but are worth saving :)

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Define tuple of inputs:
regions = ["co","md","nd","wi"]
sectors = ["agr","fof"]
other = ["a","b"]

inp = tuple(regions, sectors, other)

# Define dictionary missing some keys:
k = vcat(collect(Base.Iterators.product(inp...))...)
d = Dict(k => 1.0 for k in k[1:end-2])

function fill_zero(source::Tuple, tofill::Dict)

    # Assume all possible permutations of keys should be present
    # and determine which are missing.
    allkeys = vcat(collect(Base.Iterators.product(source...))...)
    missingkeys = setdiff(allkeys, collect(keys(tofill)))

    # Add
    [push!(tofill, fill=>0) for fill in missingkeys]
    return tofill
end

fill_zero(inp, d)

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Get dictionary values with some condition.
d = Dict(k => v > 0.0 for (k,v) in d)
[k => v > 0.0 for (k,v) in d]

# * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
# Regular expressions
# GSP FULL_CODE -> STATE_CODE, COUNTY_CODE
str = "01001"
match(r"(?<state_code>\d{2})(?<county_code>\d*)", str)

# NASS NAICS_CODE -> NAICS_CODE
str = "NAICS CLASSIFICATION: (111)"
match(r"\((?<naics_code>.*)\)", str)

str = "1114 Mushrooms, Nursery & Related Products"
match(r"(?<naics_industry_group_code>\d*) (?<naics_industry_group_desc>.*)", str)

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Test, Mocking, EzXML

Mocking.activate()

gpx_sample_file2=joinpath(@__DIR__, gpx_sample_file)
tcx_sample_file2=joinpath(@__DIR__, tcx_sample_file)
tcx_centrypark_run_file2=joinpath(@__DIR__, tcx_centrypark_run_file)
tcx_treadmill_run_file2=joinpath(@__DIR__, tcx_treadmill_run_file)

@testset "TESTSETS: Test read TCX file" begin
    @testset "CASE: Test file not XML" begin
        err, data =  TCX.parse_tcx_file(joinpath(@__DIR__, "runtests.jl"))
        @test err == TCX.CLIENT_ERROR
    end

    @testset "CASE: Test file not exists" begin
        err, data =  TCX.parse_tcx_file("non-existing-file.tcx")
        @test err == TCX.NOT_FOUND
    end

    @testset "CASE: Test file exist but not a TCX file" begin
        err, _ =  TCX.parse_tcx_file(gpx_sample_file2)
        @test err == TCX.CLIENT_TCX_ERROR
    end

    @testset "CASE: Test a TCX records run." begin
        err, data =  TCX.parse_tcx_file(tcx_sample_file2)
        @test (err == TCX.OK) & (getActivityType(data) == "Running")
    end

    @testset "CASE: Test a TCX contains track points." begin
        err, data =  TCX.parse_tcx_file(tcx_sample_file2)
        @test (err == TCX.OK) & (size(getDataFrame(data), 1) == 3427)
    end

    @testset "CASE: Test a TCX treadmill run, no logitude/latitude." begin
        err, data =  TCX.parse_tcx_file(tcx_treadmill_run_file2)
        @test (err == TCX.OK) & (getDistance2(data) == 0)
    end

    @testset "CASE: Test a TCX static running distance." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == TCX.OK) & (getDistance(data) == 27020)
    end

    @testset "CASE: Test a TCX running distance from Geodesy." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == TCX.OK) & (getDistance2(data) > 0)
    end

    @testset "CASE: Test a TCX static average pace." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == TCX.OK) & (getAveragePace(data) > 0 )
    end

    @testset "CASE: Test to get TCX duration." begin
        err, data =  TCX.parse_tcx_file(tcx_centrypark_run_file2)
        @test (err == TCX.OK) & (getDuration(data) > 0 )
    end

    @testset "CASE: Test throws error on unexpected condition" begin
        patch = @patch EzXML.readxml(path::String) = fn(throw(ReadOnlyMemoryError()))
        apply(patch) do
            @test_throws ReadOnlyMemoryError TCX.parse_tcx_file(tcx_centrypark_run_file2)
        end
    end
end


text:
function vtk_connectivity_map_lagrange(Nq1, Nq2 = 1, Nq3 = 1)
    connectivity = Array{Int, 1}(undef, Nq1 * Nq2 * Nq3)
    L = LinearIndices((1:Nq1, 1:Nq2, 1:Nq3))

    corners = (
        (1, 1, 1),
        (Nq1, 1, 1),
        (Nq1, Nq2, 1),
        (1, Nq2, 1),
        (1, 1, Nq3),
        (Nq1, 1, Nq3),
        (Nq1, Nq2, Nq3),
        (1, Nq2, Nq3),
    )
    edges = (
        (2:(Nq1 - 1), 1:1, 1:1),
        (Nq1:Nq1, 2:(Nq2 - 1), 1:1),
        (2:(Nq1 - 1), Nq2:Nq2, 1:1),
        (1:1, 2:(Nq2 - 1), 1:1, 1:1),
        (2:(Nq1 - 1), 1:1, Nq3:Nq3),
        (Nq1:Nq1, 2:(Nq2 - 1), Nq3:Nq3),
        (2:(Nq1 - 1), Nq2:Nq2, Nq3:Nq3),
        (1:1, 2:(Nq2 - 1), Nq3:Nq3),
        (1:1, 1:1, 2:(Nq3 - 1)),
        (Nq1:Nq1, 1:1, 2:(Nq3 - 1)),
        (1:1, Nq2:Nq2, 2:(Nq3 - 1)),
        (Nq1:Nq1, Nq2:Nq2, 2:(Nq3 - 1)),
    )
    faces = (
        (1:1, 2:(Nq2 - 1), 2:(Nq3 - 1)),
        (Nq1:Nq1, 2:(Nq2 - 1), 2:(Nq3 - 1)),
        (2:(Nq1 - 1), 1:1, 2:(Nq3 - 1)),
        (2:(Nq1 - 1), Nq2:Nq2, 2:(Nq3 - 1)),
        (2:(Nq1 - 1), 2:(Nq2 - 1), 1:1),
        (2:(Nq1 - 1), 2:(Nq2 - 1), Nq3:Nq3),
    )
    if Nq2 == Nq3 == 1
        @assert Nq1 > 1
        corners = (corners[1:2]...,)
        edges = (edges[1],)
        faces = ()
    elseif Nq3 == 1
        @assert Nq1 > 1
        @assert Nq2 > 1
        corners = (corners[1:4]...,)
        edges = (edges[1:4]...,)
        faces = (faces[5],)
    end

    # corners
    iter = 1
    for (i, j, k) in corners
        connectivity[iter] = L[i, j, k]
        iter += 1
    end
    # edges
    for (is, js, ks) in edges
        for k in ks, j in js, i in is
            connectivity[iter] = L[i, j, k]
            iter += 1
        end
    end
    # faces
    for (is, js, ks) in faces
        for k in ks, j in js, i in is
            connectivity[iter] = L[i, j, k]
            iter += 1
        end
    end
    # interior
    for k in 2:(Nq3 - 1), j in 2:(Nq2 - 1), i in 2:(Nq1 - 1)
        connectivity[iter] = L[i, j, k]
        iter += 1
    end

    return connectivity
end

"""
  vtk_cells_lagrange(space)

Construct a vector of `MeshCell` objects representing the elements of `space` as
an unstuctured mesh of Lagrange polynomial cells, suitable for passing to
`vtk_grid`.
"""
function vtk_cells_lagrange(
    gspace::Spaces.SpectralElementSpace2D{
        T,
        Spaces.Quadratures.ClosedUniform{Nq},
    },
) where {T, Nq}
    # TODO: this should depend on the backing DataLayouts (e.g. IJFH)
    con_map = vtk_connectivity_map_lagrange(Nq, Nq)
    [
        MeshCell(
            VTKCellTypes.VTK_LAGRANGE_QUADRILATERAL,
            ((e - 1) * Nq * Nq) .+ con_map,
        ) for e in 1:Topologies.nlocalelems(gspace)
    ]
end

"""
  vtk_cells_linear(space)

Construct a vector of `MeshCell` objects representing the elements of `space` as
an unstuctured mesh of linear cells, suitable for passing to
`vtk_grid`.
"""
function vtk_cells_linear(gridspace::Spaces.SpectralElementSpace2D)
    Nq = Spaces.Quadratures.degrees_of_freedom(gridspace.quadrature_style)
    Nh = Topologies.nlocalelems(gridspace)
    ind = LinearIndices((1:Nq, 1:Nq, 1:Nh))
    cells = [
        MeshCell(
            VTKCellTypes.VTK_QUAD,
            (
                ind[i, j, e],
                ind[i + 1, j, e],
                ind[i + 1, j + 1, e],
                ind[i, j + 1, e],
            ),
        ) for i in 1:(Nq - 1), j in 1:(Nq - 1), e in 1:Nh
    ]
    return vec(cells)
end

function vtk_cells_linear(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace2D)
    Nq = Spaces.Quadratures.degrees_of_freedom(
        gridspace.horizontal_space.quadrature_style,
    )
    Nh = Topologies.nlocalelems(gridspace.horizontal_space)
    Nv = Spaces.nlevels(gridspace)
    ind = LinearIndices((1:Nv, 1:Nq, 1:Nh)) # assume VIFH
    cells = [
        MeshCell(
            VTKCellTypes.VTK_QUAD,
            (
                ind[v, i, e],
                ind[v + 1, i, e],
                ind[v + 1, i + 1, e],
                ind[v, i + 1, e],
            ),
        ) for v in 1:(Nv - 1), i in 1:(Nq - 1), e in 1:Nh
    ]
    return vec(cells)
end
function vtk_cells_linear(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace3D)
    Nq = Spaces.Quadratures.degrees_of_freedom(
        gridspace.horizontal_space.quadrature_style,
    )
    Nh = Topologies.nlocalelems(gridspace.horizontal_space)
    Nv = Spaces.nlevels(gridspace)
    ind = LinearIndices((1:Nv, 1:Nq, 1:Nq, 1:Nh)) # assumes VIJFH
    cells = [
        MeshCell(
            VTKCellTypes.VTK_HEXAHEDRON,
            (
                ind[v, i, j, e],
                ind[v + 1, i, j, e],
                ind[v + 1, i + 1, j, e],
                ind[v, i + 1, j, e],
                ind[v, i, j + 1, e],
                ind[v + 1, i, j + 1, e],
                ind[v + 1, i + 1, j + 1, e],
                ind[v, i + 1, j + 1, e],
            ),
        ) for v in 1:(Nv - 1), i in 1:(Nq - 1), j in 1:(Nq - 1), e in 1:Nh
    ]
    return vec(cells)
end

"""
    vtk_grid_space(space::ClimaCore.Spaces.AbstractSpace)

The space for the grid used by VTK, for any field on `space`.

This generally does two things:
 - Modifies the horizontal space to use a `ClosedUniform` quadrature rule, which
   will use equispaced nodal points in the reference element. This is required
   for using [VTK Lagrange elements](https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/)
   (see [1](https://discourse.paraview.org/t/node-positions-of-high-order-lagrange-quadrilateral-cells/7012/3)).
 - Modifies the vertical space to be on the faces.
"""
function vtk_grid_space(space::Spaces.SpectralElementSpace1D)
    if space.quadrature_style isa Spaces.Quadratures.ClosedUniform
        return space
    end
    Nq = Spaces.Quadratures.degrees_of_freedom(space.quadrature_style)
    lagrange_quad = Spaces.Quadratures.ClosedUniform{Nq}()
    return Spaces.SpectralElementSpace1D(space.topology, lagrange_quad)
end
function vtk_grid_space(space::Spaces.SpectralElementSpace2D)
    if space.quadrature_style isa Spaces.Quadratures.ClosedUniform
        return space
    end
    Nq = Spaces.Quadratures.degrees_of_freedom(space.quadrature_style)
    lagrange_quad = Spaces.Quadratures.ClosedUniform{Nq}()
    return Spaces.SpectralElementSpace2D(space.topology, lagrange_quad)
end
function vtk_grid_space(space::Spaces.FaceExtrudedFiniteDifferenceSpace)
    # this will need to be updated for warped meshes
    horizontal_space = vtk_grid_space(space.horizontal_space)
    vertical_space = Spaces.FaceFiniteDifferenceSpace(space.vertical_topology)
    return Spaces.ExtrudedFiniteDifferenceSpace(
        horizontal_space,
        vertical_space,
    )
end
function vtk_grid_space(space::Spaces.CenterExtrudedFiniteDifferenceSpace)
    return vtk_grid_space(Spaces.FaceExtrudedFiniteDifferenceSpace(space))
end

vtk_grid_space(field::Fields.Field) = vtk_grid_space(axes(field))
vtk_grid_space(fields::NamedTuple) = vtk_grid_space(first(fields))
vtk_grid_space(fieldvec::Fields.FieldVector) =
    vtk_grid_space(Fields._values(fieldvec))

"""
    vtk_cell_space(gridspace::ClimaCore.Spaces.AbstractSpace)

Construct a space for outputting cell data, when using outputting a grid `gridspace`.
be stored.

This generally does two things:
 - Modifies the horizontal space to use a `Uniform` quadrature rule, which
   will use equispaced nodal points in the reference element (excluding the boundary).
 - Modifies the vertical space to be on the centers.
"""
function vtk_cell_space(gridspace::Spaces.SpectralElementSpace1D)
    @assert gridspace.quadrature_style isa Spaces.Quadratures.ClosedUniform
    Nq = Spaces.Quadratures.degrees_of_freedom(gridspace.quadrature_style)
    quad = Spaces.Quadratures.Uniform{Nq - 1}()
    return Spaces.SpectralElementSpace1D(gridspace.topology, quad)
end
function vtk_cell_space(gridspace::Spaces.SpectralElementSpace2D)
    @assert gridspace.quadrature_style isa Spaces.Quadratures.ClosedUniform
    Nq = Spaces.Quadratures.degrees_of_freedom(gridspace.quadrature_style)
    quad = Spaces.Quadratures.Uniform{Nq - 1}()
    return Spaces.SpectralElementSpace2D(gridspace.topology, quad)
end
function vtk_cell_space(gridspace::Spaces.FaceExtrudedFiniteDifferenceSpace)
    # this will need to be updated for warped meshes
    horizontal_space = vtk_cell_space(gridspace.horizontal_space)
    vertical_space =
        Spaces.CenterFiniteDifferenceSpace(gridspace.vertical_topology)
    return Spaces.ExtrudedFiniteDifferenceSpace(
        horizontal_space,
        vertical_space,
    )
end


text:
# These are the additions to the types_and_consts.jl file not from clang.jl

function Base.convert(::Type{Matrix}, J::DlsMat)
    _dlsmat = unsafe_load(J)
    # own is false as memory is allocated by sundials
    unsafe_wrap(Array, _dlsmat.data, (_dlsmat.M, _dlsmat.N), own=false)
end

CVRhsFn_wrapper(fp::CVRhsFn) = fp
CVRhsFn_wrapper(f) = @cfunction($f,Cint,(realtype,N_Vector,N_Vector,Ptr{Cvoid})).ptr

ARKRhsFn_wrapper(fp::ARKRhsFn) = fp
ARKRhsFn_wrapper(f) = @cfunction($f,Cint,(realtype,N_Vector,N_Vector,Ptr{Cvoid})).ptr

CVRootFn_wrapper(fp::CVRootFn) = fp
CVRootFn_wrapper(f) = @cfunction($f,Cint,(realtype,N_Vector,Ptr{realtype},Ptr{Cvoid})).ptr

CVQuadRhsFn_wrapper(fp::CVQuadRhsFn) = fp
CVQuadRhsFn_wrapper(f) = @cfunction($f,Cint,(realtype,N_Vector,N_Vector,Ptr{Cvoid})).ptr

IDAResFn_wrapper(fp::IDAResFn) = fp
IDAResFn_wrapper(f) = @cfunction($f,Cint,(realtype,N_Vector,N_Vector,N_Vector,Ptr{Cvoid})).ptr

IDARootFn_wrapper(fp::IDARootFn) = fp
IDARootFn_wrapper(f) = @cfunction($f,Cint,(realtype,N_Vector,N_Vector,Ptr{realtype},Ptr{Cvoid})).ptr

KINSysFn_wrapper(fp::KINSysFn) = fp
KINSysFn_wrapper(f) = @cfunction($f,Cint,(N_Vector,N_Vector,Ptr{Cvoid})).ptr

function Base.convert(::Type{Matrix}, J::SUNMatrix)
    _sunmat = unsafe_load(J)
    _mat = convert(SUNMatrixContent_Dense, _sunmat.content)
    mat = unsafe_load(_mat)
    # own is false as memory is allocated by sundials
    unsafe_wrap(Array, mat.data, (mat.M, mat.N), own=false)
end

function Base.convert(::Type{SparseMatrixCSC}, J::SUNMatrix)
    _sunmat = unsafe_load(J)
    _mat = convert(SUNMatrixContent_Sparse, _sunmat.content)
    mat = unsafe_load(_mat)
    # own is false as memory is allocated by sundials
    # TODO: Get rid of allocation for 1-based index change
    rowval = unsafe_wrap(Array, mat.indexvals, (mat.NNZ), own=false)
    colptr = unsafe_wrap(Array, mat.indexptrs, (mat.NP+1), own=false)
    m = mat.M
    n = mat.N
    nzval = unsafe_wrap(Array,mat.data, (mat.NNZ), own=false)
    SparseMatrixCSC(m,n,colptr,rowval,nzval)
end

abstract type SundialsMatrix end
struct DenseMatrix <: SundialsMatrix end
struct BandMatrix <: SundialsMatrix end
struct SparseMatrix <: SundialsMatrix end

abstract type SundialsLinSol end
struct Dense <: SundialsLinSol end
struct Band <: SundialsLinSol end
struct SPGMR <: SundialsLinSol end
struct SPFGMR <: SundialsLinSol end
struct SPBCGS <: SundialsLinSol end
struct PCG <: SundialsLinSol end
struct PTFQMR <: SundialsLinSol end
struct KLU <: SundialsLinSol end

abstract type StiffnessChoice end
struct Explicit <: StiffnessChoice end
struct Implicit <: StiffnessChoice end


text:
function test_evaluate_explicit()
    println("-- Testing evaluate function for explicit feedback")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)
    recommender = MF(data, 2)
    build!(recommender)

    truth_m = [1 3 4 1 2 1 2 4
               1 2 4 1 3 2 2 3
               5 2 3 3 4 5 2 1]
    truth_data = DataAccessor(truth_m)

    # average error should be less than 2.5 in the 1-to-5 rating
    @test evaluate(recommender, truth_data, RMSE()) < 2.5
end

function test_evaluate_implicit()
    println("-- Testing evaluate function for implicit feedback")

    m = [NaN 1 NaN 0 0 0 NaN 1
         0 0 NaN NaN 1 0 NaN 1
         NaN 0 1 1 NaN 1 NaN 0]
    data = DataAccessor(m)
    recommender = MF(data, 2)
    build!(recommender)

    truth_m = [0 1 1 0 0 0 0 1
               0 0 1 0 1 0 0 1
               1 0 1 1 1 1 0 0]
    truth_data = DataAccessor(truth_m)

    # average error should be less than or equal to 0.5 in the [0, 1] scale
    @test evaluate(recommender, truth_data, Recall(), 4) <= 0.5
end

test_evaluate_explicit()
test_evaluate_implicit()


text:
using AbstractGPs
using BayesianLinearRegressors
using Distributions
using FiniteDifferences
using LinearAlgebra
using Random
using Test
using Zygote
using PDMats

using BayesianLinearRegressors: BayesianLinearRegressor, posterior, marginals, cov, mean
using FiniteDifferences: j′vp

include("test_utils.jl")

@testset "BayesianLinearRegressors" begin
    include("bayesian_linear_regression.jl")
    include("basis_function_regression.jl")
    include("sampling_functions.jl")
end


text:
using IterativeSolvers
using FactCheck
using Base.Test

facts("lsqr") do

context("Small dense matrix") do
    A = rand(10, 5)
    b = rand(10)
    x, = lsqr(A, b)
    @fact norm(x - A\b) --> less_than(√eps())
end

context("SOL test") do
    # Test adapted from the BSD-licensed Matlab implementation at
    #    http://www.stanford.edu/group/SOL/software/lsqr.html
    #              Michael Saunders, Systems Optimization Laboratory,
    #              Dept of MS&E, Stanford University.
    #-----------------------------------------------------------------------
    function lsqrSOLtest( m, n, damp )
        fmul  = (out, b) -> Aprodxxx!(out,b,1,m,n)
        fcmul = (out, b) -> Aprodxxx!(out,b,2,m,n)
        A = MatrixCFcn{Int}(m, n, fmul, fcmul)
        xtrue = n:-1:1
        b = float(A*xtrue)
        x, = lsqr(A, b, atol = 1e-6, btol = 1e-6, conlim = 1e10, maxiter = 10n)
        r = b - A*x
        @fact norm(r) --> less_than_or_equal(1e-4)
        x
    end

    function Aprodxxx!(y, x, mode, m, n )
        # if mode = 1, computes y = A*x
        # if mode = 2, computes y = A'*x
        # for some matrix  A.
        #
        # This is a simple example for testing  LSQR.
        # It uses the leading m*n submatrix from
        # A = [ 1
        #       1 2
        #         2 3
        #           3 4
        #             ...
        #               n ]
        # suitably padded by zeros.
        #
        # 11 Apr 1996: First version for distribution with lsqr.m.
        #              Michael Saunders, Dept of EESOR, Stanford University.

        if mode == 1
            y[1] = x[1]
            for i = 2:n
                y[i] = i*x[i] + (i-1)*x[i-1]
            end
            for i = n+1:m
                y[i] = 0
            end
        else
            mn = min(m, n)
            for i = 1:mn-1
                y[i] = i*(x[i]+x[i+1])
            end
            y[mn] = mn*x[mn]
            for i = m+1:n
                y[i] = 0
            end
        end
        y
    end

    x = lsqrSOLtest( 10,10, 0   )
    x = lsqrSOLtest( 20,10, 0   )
    x = lsqrSOLtest( 20,10, 0.1 )
end

context("Issue 64") do
    srand(224)
    n=10
    A=sprand(n,n,.5)
    b=rand(n)

    x, ch = lsqr(A,b,maxiter=100)
    resnorm = norm(A*x - b)
    @fact resnorm --> less_than(√eps())
    @fact ch.isconverged --> true
    @fact last(ch[:resnorm]) --> roughly(resnorm, atol=√eps())
end

end


text:
## Defining and showing some examples
cheers_1 = NVModel(truncated(Normal(90, 30), 0, 180), 1, 5, 0)
@test_nowarn show(cheers_1)

cheers_2 = NVModel(Uniform(0, 60), 4.0, 5.5)
@test_nowarn show(cheers_2)

cheers_4 = NVModel(cost=10, price=42, demand=DiscreteNonParametric([0,1,2],[.3,.5,.2]))
@test_nowarn show(cheers_4)

bagel = NVModel(Normal(150, 50), 1, 4)
@test_nowarn show(bagel)

xs = vec([5_000	10_000	15_000	20_000	25_000	30_000	35_000	40_000	45_000	50_000	55_000	60_000	65_000	70_000	75_000])
ps = vec([0.0181	0.0733	0.1467	0.1954	0.1954	0.1563	0.1042	0.0595	0.0298	0.0132	0.0053	0.0019	0.0006	0.0002	0.0001])
elvis = NVModel(DiscreteNonParametric(xs,ps), 6, 12, 2.5)
@test_nowarn show(elvis)

ksweet = NVModel(Binomial(2000, 0.02), 28, 45, 27)
@test_nowarn show(ksweet)

beer = NVModel(Uniform(0, 300), holding = 0.5, backorder = 1)
@test_nowarn show(beer)

maxtest = NVModel(Normal(90, 30), 0, 3, fixcost = 100, q_min=90, q_max = 120)
@test_nowarn show(maxtest)

mintest = NVModel(Normal(90, 30), 1, 0, q_min=90, q_max = 120)
@test_nowarn show(mintest)

unboundtest = NVModel(Normal(90, 30), 0, 3, q_min=90)
@test_nowarn show(unboundtest)

thesaurus_nvm = NVModel(cost = 1.15, price = 2.75, demand = Normal(18, 6),
                        substitute = 2.75-1.15, backorder = 0.5, salvage = 1.15,
                        holding = 1.15 * 0.2 / 12)
@test_nowarn show(thesaurus_nvm)

## Testing solve and showing result
res = solve(cheers_1)
@test_nowarn show(res)
nvm = NewsvendorModel.nvmodel(res)
@test_nowarn show(nvm)

res = solve(cheers_1, rounded=false)
@test q_opt(res) == q_opt(solve(cheers_1, rounded=false))

@test_nowarn show(solve(cheers_1) )
@test_nowarn show(solve(cheers_2) )
@test_nowarn show(solve(cheers_4) )
@test_nowarn show(solve(bagel) )
@test_nowarn show(solve(elvis) )
@test_nowarn show(solve(elvis, rounded=false) )
@test_nowarn show(solve(ksweet) )
@test_nowarn show(solve(beer) )
@test_nowarn show(solve(maxtest) )
@test_nowarn show(solve(mintest) )
@test_nowarn show(solve(unboundtest) )
@test_nowarn show(solve(thesaurus_nvm) )


text:
module OpenSSLCrypto

include("Crypto.jl")

end


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 86f770fe-74a1-11eb-01f7-5b3ecf057124
begin
	import Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(name="ImageIO", version="0.5"),
		Pkg.PackageSpec(name="ImageShow", version="0.2"),
		Pkg.PackageSpec(name="FileIO", version="1.6"),
		Pkg.PackageSpec(name="PNGFiles", version="0.3.6"),
		Pkg.PackageSpec(name="Colors", version="0.12"),
		Pkg.PackageSpec(name="ColorVectorSpace", version="0.8"),

		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
		Pkg.PackageSpec(name="Unitful", version="1.6"),
		Pkg.PackageSpec(name="ImageFiltering", version="0.6"),
		Pkg.PackageSpec(name="OffsetArrays", version="1.6"),
		Pkg.PackageSpec(name="Plots", version="1.10")
	])

	using PlutoUI
	using Colors, ColorVectorSpace, ImageShow, FileIO
	using Unitful
	using ImageFiltering
	using OffsetArrays
	using Plots
end

# ╔═╡ b310756a-af08-48b0-ae10-ee2e8dd0c968
filter!(LOAD_PATH) do path
	path != "@v#.#"
end;

# ╔═╡ febfa62a-74fa-11eb-2fe6-df7de43ef4b6
html"""

Section 1.4

 Transformations with Images

Lecture Video

"""

# ╔═╡ 8d389d80-74a1-11eb-3452-f38eff03483b
PlutoUI.TableOfContents(aside=true)

# ╔═╡ 9f1a72da-7532-11eb-079c-b7baccc6614a
md"""
#### Intializing packages

_When running this notebook for the first time, this could take up to 15 minutes. Hang in there!_
"""

# ╔═╡ 4d332c7e-74f8-11eb-1f49-a518246d1db8
md"""
# Announcement: Lectures will be nearly an hour
"""

# ╔═╡ f7689472-74a8-11eb-32a1-8379ae5c88e1
rotabook = PlutoUI.Resource("https://i.gr-assets.com/images/S/compressed.photo.goodreads.com/books/1348902666l/1646354.jpg")

# ╔═╡ 0f2f9004-74a8-11eb-01a2-973dbe80f166
md"""
##  **Never run overtime** (a microcentury with UnitFul)

Running overtime is the one unforgivable error a lecturer can make.
After fifty minutes (one microcentury as von Neumann used to say)
everybody's attention will turn elsewhere even if we are trying to prove
the Riemann hypothesis. One minute overtime can destroy the best of
lectures. (from "Indiscrete Thoughts" by Rota, Chpt 18, 10 Lessons I Wish I Had Been Taught)
"""

# ╔═╡ 962143a8-74a7-11eb-26c3-c10548f326ee
century = 100u"yr" #  a u"yr" is a special kind of string denoting a unit of a year

# ╔═╡ c2964c80-74f8-11eb-3a74-b1bdd9e4ae02
century * 2

# ╔═╡ caf488d8-74f8-11eb-0075-0586d66c23c1
century/200

# ╔═╡ 02dd4a02-74f9-11eb-3d1e-53d83cee8062
century^2

# ╔═╡ 10ef13d2-74f9-11eb-2849-fb9f83db6ae9
g = 9.8u"m"/u"s"^2

# ╔═╡ b76a56f4-74a9-11eb-1739-fbfc5e4958e8
uconvert(u"minute", century / 10e3) # convert into minutes the value of a microcentury

# ╔═╡ 77fbf18a-74f9-11eb-1d9e-3f9d2097388f
PotentialEnergy = (10u"kg") * g * (50u"m")

# ╔═╡ bcb69db6-74f9-11eb-100a-29d1d23963ab
uconvert( u"J", PotentialEnergy)

# ╔═╡ fc70c4d2-74f8-11eb-33f5-539c278ed6b6
md"""
Adding units to numbers **just works** in Julia, and furthermore, does not slow down execution.  We are sneaking in an example of the power of generic programming and Julia's type system, some of the underlying technology that makes us love working with Julia.  More on this later in the book.  Meanwhile if this helps you do your problem sets in some other class, go for it.
"""

# ╔═╡ 2f7cde78-74a2-11eb-1e2f-81b5b2465819
md"""
# Reminder

**Try your own pictures everywhere!**
"""

# ╔═╡ e099815e-74a1-11eb-1541-033f6abe9f8e
md"""
# Transforming Images
"""

# ╔═╡ e82a4dd8-74b0-11eb-1108-6b09e67a80c1
md"""
## 2.1. Downsampling / Upsampling
"""

# ╔═╡ 39552b7a-74fb-11eb-04e0-3981ada52c92
md"""
How can we pixelate a corgi? Found this cute picture online, but we'll pixelate
a real corgi.
"""

# ╔═╡ 14f2b85e-74ad-11eb-2682-d9de646aedf3
pixelated_corgi = load("pixelated_corgi.png")

# ╔═╡ 516e73e2-74fb-11eb-213e-9dbd9472e0db
philip = load(download("https://user-images.githubusercontent.com/6933510/107239146-dcc3fd00-6a28-11eb-8c7b-41aaf6618935.png"))

# ╔═╡ b5d0ef90-74fb-11eb-3126-792f954c7be7
@bind r Slider(1:40, show_value=true, default=40)

# ╔═╡ 754c3704-74fb-11eb-1199-2b9798d7251f
downsample_philip = philip[1:r:end, 1:r:end]

# ╔═╡ 9eb917ba-74fb-11eb-0527-15e981ce9c6a
upsample_philip = kron(downsample_philip, fill(1, r, r))

# ╔═╡ 486d3022-74ff-11eb-1865-e15436bd9aad
md"""
  Note the use of kron and fill. See [Wikipedia Kron](https://en.wikipedia.org/wiki/Kronecker_product)
"""

# ╔═╡ b9da7332-74ff-11eb-241b-fb87e77d646a
md"""
Exercise: Use the nose selection tool from Section 1.1 to pixelate a rectangle of an image.  Warning: you'll have to worry about sizes if not exact multiples.
"""

# ╔═╡ 339ccfca-74b1-11eb-0c35-774da6b189ed
md"""
## 2.2 Linear Combinations (Combining Images)
"""

# ╔═╡ 8711c698-7500-11eb-2505-d35a4de169b4
md"""
One big idea in mathematics is the [linear combination](https://en.wikipedia.org/wiki/Linear_combination).
The idea combines
- scaling an object
- combining two or more objects
by combining scaled versions of multiple objects.
"""

# ╔═╡ 84350cb8-7501-11eb-095e-8f1a7e015f25
md"""
Let's scale some corgis.
"""

# ╔═╡ 91a1bca4-74aa-11eb-3917-1dfd73d0ad9c
corgis = load(download("https://user-images.githubusercontent.com/6933510/108605549-fb28e180-73b4-11eb-8520-7e29db0cc965.png"))

# ╔═╡ 8e698bdc-7501-11eb-1d2e-c336ccbde0b0
@bind c Slider(0:.1:3, show_value=true, default=1)

# ╔═╡ ab2bc924-7501-11eb-03ba-8dfc1ffe3f36
c .* corgis  # scaling the corgis changes intensity

# ╔═╡ e11d6300-7501-11eb-239a-135596309d20
md"""
 You might wonder about the **dot times** or **pointwise times**. You can
delete the dot, but it is recommended for clarity
and performance.  The dot emphasizes that the multiplication by c is happening
pixel by pixel or that the scalar is being "broadcast" to every pixel.
"""

# ╔═╡ 9a66c07e-7503-11eb-3127-7fce91b3a24a
md"""
Scaling too far saturates the image.  (Any r,g,b ≥ 1, saturates at 1.)
"""

# ╔═╡ 47d40406-7502-11eb-2f43-cd5c848f25a6
md"""
We need another image.  We could grab one from somewhere or we can just transform the one we have.  Let's do the latter and turn the corgis upsidedown.
"""

# ╔═╡ 9ce0b980-74aa-11eb-0678-01209451fb65
upsidedown_corgis = corgis[end: -1: 1, :]

# ╔═╡ 68821bf4-7502-11eb-0d3c-03d7a00fdba4
md"""
Now let's scaled version of the two images to see what that does.
"""

# ╔═╡ 447e7c9e-74b1-11eb-27ea-71aa4338b11a
(.5 * upsidedown_corgis .+ .5 * corgis)

# ╔═╡ c9dff6f4-7503-11eb-2715-0bf9d3ece9e1
md"""
### Convex Combinations
"""

# ╔═╡ d834103c-7503-11eb-1a94-1fbad43801ff
md"""
If all the coefficients are positive and add to 1, we say we have a **convex combination**.  Let's take α and (1-α) as the two coefficients adding to 1, and
scale the two corgi pictures with different α's, thereby giving different weights to the rightside-up and upside-down corgis.
"""

# ╔═╡ aa541288-74aa-11eb-1edc-ab6d7786f271
@bind α Slider(0:.01:1, show_value=true, default = 1.0)

# ╔═╡ c9dcac48-74aa-11eb-31a6-23357180c1c8
α .* corgis .+ (1-α) .* upsidedown_corgis

# ╔═╡ 30b1c1f0-7504-11eb-1be7-a9463caea809
md"""
The moment I did this with α = .5, I noticed my brain's tendency to see the
rightsisde-up corgis even though both have equal weight.  For me maybe
around α = .39 which gives weight .61 to the upside-down corgis "feels" balanced
to me.  I think this is what the field of psychology called psychometrics
tries to measure -- perhaps someone can tell me if there are studies of the
brain's tendency to use world experience to prefer rightside-up corgis,
and in particular to put a numerical value to this tendency.
"""

# ╔═╡ 1fe70e38-751b-11eb-25b8-c741e1726613
md"""
10 seconds with google and I found there is a thing about faces:
[The Face Inversion effect](https://en.wikipedia.org/wiki/Face_inversion_effect#:~:text=The%20face%20inversion%20effect%20is,same%20for%20non%2Dfacial%20objects.&text=The%20most%20supported%20explanation%20for,is%20the%20configural%20information%20hypothesis)
and also the [Thatcher Effect](https://en.wikipedia.org/wiki/Thatcher_effect#:~:text=The%20Thatcher%20effect%20or%20Thatcher,obvious%20in%20an%20upright%20face) seems related.

... the article suggests objects don't suffer in the same way as faces,
so I put forth that the phenomenon applies to corgi faces as much as human faces,
suggesting maybe that corgi faces are processed in the **face processing** part of the brain,not the **object processing** part of the brain.

Corgis are human, after all, right?

(Note, this is 5 minutes of armchair science, not a professional opinion.)
"""

# ╔═╡ 215291ec-74a2-11eb-3476-0dab43fd5a5e
md"""
## 2.3 Fun with Photoshop (What does "filter" mean in this context?)
"""

# ╔═╡ 61db42c6-7505-11eb-1ddf-05e906234572
md"""
[Photshop Filter Reference](https://helpx.adobe.com/photoshop/using/filter-effects-reference.html)
"""

# ╔═╡ cdd4cffc-74b1-11eb-1aa4-e333cb8601d1
md"""
Let's play with photoshop if for no other reason, let's see what image transformations are available considered useful by the pros.

Some worth emphasizing are
1. Blur
2. Sharpen
3. Stylize -> Find Edges
3. Pixelate
4. Distort

Some of these transformations (e.g. Blur, Sharpen, Find Edges) are examples of
convolutions which are very efficient, and show up these days in
machine learning particularly in image recognition.
"""

# ╔═╡ 7489a570-74a3-11eb-1d0b-09d41604ffe1
md"""
## 2.4 Image Filtering (convolutions)
"""

# ╔═╡ 8a8e3f5e-74b2-11eb-3eed-e5468e573e45
md"""
Last semester Grant Sanderson (3Blue1Brown) lectured in this course.  This lecture on convolutions in image processing was popular.  Let's watch an excerpt (from 1:04 to 2:48).  (We pick a few exercepts, but we wouldn't blame you if you just wanted to
watch the whole video.)
"""

# ╔═╡ 5864294a-74a5-11eb-23ef-f38a582f2c2d
html"""

"""

# ╔═╡ fa9c465e-74b2-11eb-2f3c-4be0e7f93bb5
md"""
### Definition of convolutions and kernels
"""

# ╔═╡ 4fab4616-74b0-11eb-0088-6b50237d7d54
md"""
[Wikipedia Page on Kernels]
(https://en.wikipedia.org/wiki/Kernel_(image_processing)#Details)
"""

# ╔═╡ 275bf7ac-74b3-11eb-32c3-cda1e4f1f8c2
md"""
### Computer Science: Complexity

The number of multiplications = (Number of Pixels in the Image) * (Number of Cells in the kernel)
"""

# ╔═╡ 537c54e4-74b3-11eb-341f-951b4a1e0b40
md"""
Thought Problem: Why are small kernels better than large kernels from a complexity viewpoint?
"""

# ╔═╡ c6e340ee-751e-11eb-3ca7-69595b3693b7
md"""
### Computer Science: Architectures: GPUs or Graphical Processing Units

Some important computations can be greatly accelerated through the use of specialized hardware such as the GPU processors that were originally designed as image renderers but it has turned out that these processors can be quite fast at other very regular computations.  Convolutions is a very GPU friendly operation due to its regular structure.
"""

# ╔═╡ 54448d18-7528-11eb-209a-9717affa0d02
 kernelize(M) = OffsetArray(M, -1: 1, -1: 1)

# ╔═╡ acbc563a-7528-11eb-3c38-75a5b66c9241
begin
	identity = [0 0 0 ; 0 1 0 ; 0 0 0]
	edge_detect = [0 -1 0; -1 4 -1; 0 -1 0]
	sharpen = identity .+ edge_detect  # Superposition!
	box_blur = [1 1 1;1 1 1;1 1 1]/9
	∇x = [-1 0 1;-1 0 1;-1 0 1]/2 # centered deriv in x
	∇y = ∇x'

	kernels = [identity, edge_detect, sharpen, box_blur, ∇x, ∇y]
	kernel_keys =["identity", "edge_detect", "sharpen", "box_blur", "∇x", "∇y"]
	selections = kernel_keys .=> kernel_keys
	kernel_matrix = Dict(kernel_keys .=> kernels)
	md"$(@bind kernel_name Select(selections))"
end

# ╔═╡ 995392ee-752a-11eb-3394-0de331e24f40
kernel_matrix[kernel_name]

# ╔═╡ d22903d6-7529-11eb-2dcd-132cd27104c2
[imfilter(corgis, kernelize(kernel_matrix[kernel_name])) Gray.(1.5 .* abs.(imfilter( corgis, kernelize(kernel_matrix[kernel_name]))))]

# ╔═╡ 844ed844-74b3-11eb-2ee1-2de664b26bc6
md"""
  ### Gaussian Filter
"""

# ╔═╡ 4ffe927c-74b4-11eb-23a7-a18d7e51c75b
md"""
In our next Grant Sanderson segment from Fall 2020 (4:35 to 7:00), we hear about
convolving images with a Gaussian kernel.
"""

# ╔═╡ 91109e5c-74b3-11eb-1f31-c50e436bc6e0
html"""

"""

# ╔═╡ 34109062-7525-11eb-10b3-d59d3a6dfda6
round.(Kernel.gaussian(1), digits=3)

# ╔═╡ 9ab89a3a-7525-11eb-186d-29e4b61deb7f
md"""
We could have defined this ourselves with calls to the exponential function.
"""

# ╔═╡ 50034058-7525-11eb-345b-3334e71ac50e
begin
	G = [exp( -(i^2+j^2)/2) for i=-2:2, j=-2:2]
	round.(G ./ sum(G), digits=3)
end

# ╔═╡ c0aec7ae-7505-11eb-2822-a151aad48fc9
md"""
This is often known as Gausssian blur to emphasize the result of this operation.
[Adobe on Gaussian blur](https://www.adobe.com/creativecloud/photography/discover/gaussian-blur.html).
"""

# ╔═╡ 628aea22-7521-11eb-2edc-39ac62683aea
md"""
Focus around 5:23
"""

# ╔═╡ 99eeb11c-7524-11eb-2154-df7d84976445
@bind gparam Slider(0:9, show_value=true, default=1)

# ╔═╡ 2ddcfb90-7520-11eb-2df7-172d07118b7e
kernel = Kernel.gaussian(gparam)

# ╔═╡ d62496b0-7524-11eb-3410-7177e7c7f8eb
plotly()

# ╔═╡ 6aa8a76e-7524-11eb-22b5-015aab4191b0
surface([kernel;])

# ╔═╡ ee93eeb2-7524-11eb-342d-0343d8aebf59
md"""
Note: black lines are contours
"""

# ╔═╡ 662d73b6-74b3-11eb-333d-f1323a001000
md"""
### Computer Science: Data Structure: Offset Arrays
"""

# ╔═╡ d127303a-7521-11eb-3507-7341a416211f
kernel[0, 0]

# ╔═╡ d4581b56-7522-11eb-2c15-991c0c790e67
kernel[-2, 2]

# ╔═╡ 40c15c3a-7523-11eb-1f2a-bd90b127dad2
M = [ 1  2  3  4  5
	  6  7  8  9 10
	 11 12 13 14 15]

# ╔═╡ 08642690-7523-11eb-00dd-63d4cf6513dc
Z = OffsetArray(M, -1:1, -2:2)

# ╔═╡ deac4cf2-7523-11eb-2832-7b9d31389b08
the_indices = [c.I for c ∈ CartesianIndices(Z)]

# ╔═╡ 32887dfa-7524-11eb-35cd-051eff594fa9
Z[1, -2]

# ╔═╡ 0f765670-7506-11eb-2a37-931b15bb387f
md"""
## 2.5. Discrete vs Continuous
"""

# ╔═╡ 82737d28-7507-11eb-1e39-c7dc12e18882
md"""
Some folks only like discrete objects, others continuous.  The computer makes clear what many mathematicians already know, that while different language has evolved to describe discrete objects vs continuous objects, often the underlying conceptual idea is similar or the same.  Here is one analogy:
"""

# ╔═╡ 40d538b2-7506-11eb-116b-efeb16b3478d
md"""
### Blurring Kernels :: Integrals  ≡ Sharpening Kernels :: Derivatives
"""

# ╔═╡ df060a88-7507-11eb-034b-5346d67a0e0d
md"""
Think about integrals vs derivatives in one dimension.
If you replace f(x) with g(x) = ∫ f(t) dt for x-r ≤ t ≤ x+r, that will blur or smooth out the features of f.  However if you take the derivative,you will emphasize the changes, i.e., you will sharpen or "edge-detect."
"""

# ╔═╡ 60c8db60-7506-11eb-1468-c989809c933a
md"""
## 2.6 Respect my Boundaries
"""

# ╔═╡ 8ed0be60-7506-11eb-2769-5f7da1c66243
md"""
Applying the convolution on a boundary requires special thought because it is literally an **edge case**.  Once again Grant said this so very well: (2:53-4:19)
"""

# ╔═╡ b9d636da-7506-11eb-37a6-3116d47b2787
html"""

"""

# ╔═╡ Cell order:
# ╟─febfa62a-74fa-11eb-2fe6-df7de43ef4b6
# ╟─8d389d80-74a1-11eb-3452-f38eff03483b
# ╟─9f1a72da-7532-11eb-079c-b7baccc6614a
# ╠═86f770fe-74a1-11eb-01f7-5b3ecf057124
# ╟─b310756a-af08-48b0-ae10-ee2e8dd0c968
# ╟─4d332c7e-74f8-11eb-1f49-a518246d1db8
# ╟─f7689472-74a8-11eb-32a1-8379ae5c88e1
# ╟─0f2f9004-74a8-11eb-01a2-973dbe80f166
# ╠═962143a8-74a7-11eb-26c3-c10548f326ee
# ╠═c2964c80-74f8-11eb-3a74-b1bdd9e4ae02
# ╠═caf488d8-74f8-11eb-0075-0586d66c23c1
# ╠═02dd4a02-74f9-11eb-3d1e-53d83cee8062
# ╠═10ef13d2-74f9-11eb-2849-fb9f83db6ae9
# ╠═b76a56f4-74a9-11eb-1739-fbfc5e4958e8
# ╠═77fbf18a-74f9-11eb-1d9e-3f9d2097388f
# ╠═bcb69db6-74f9-11eb-100a-29d1d23963ab
# ╟─fc70c4d2-74f8-11eb-33f5-539c278ed6b6
# ╟─2f7cde78-74a2-11eb-1e2f-81b5b2465819
# ╟─e099815e-74a1-11eb-1541-033f6abe9f8e
# ╟─e82a4dd8-74b0-11eb-1108-6b09e67a80c1
# ╟─39552b7a-74fb-11eb-04e0-3981ada52c92
# ╠═14f2b85e-74ad-11eb-2682-d9de646aedf3
# ╠═516e73e2-74fb-11eb-213e-9dbd9472e0db
# ╠═b5d0ef90-74fb-11eb-3126-792f954c7be7
# ╠═754c3704-74fb-11eb-1199-2b9798d7251f
# ╠═9eb917ba-74fb-11eb-0527-15e981ce9c6a
# ╟─486d3022-74ff-11eb-1865-e15436bd9aad
# ╟─b9da7332-74ff-11eb-241b-fb87e77d646a
# ╟─339ccfca-74b1-11eb-0c35-774da6b189ed
# ╟─8711c698-7500-11eb-2505-d35a4de169b4
# ╟─84350cb8-7501-11eb-095e-8f1a7e015f25
# ╠═91a1bca4-74aa-11eb-3917-1dfd73d0ad9c
# ╠═8e698bdc-7501-11eb-1d2e-c336ccbde0b0
# ╠═ab2bc924-7501-11eb-03ba-8dfc1ffe3f36
# ╟─e11d6300-7501-11eb-239a-135596309d20
# ╟─9a66c07e-7503-11eb-3127-7fce91b3a24a
# ╟─47d40406-7502-11eb-2f43-cd5c848f25a6
# ╠═9ce0b980-74aa-11eb-0678-01209451fb65
# ╟─68821bf4-7502-11eb-0d3c-03d7a00fdba4
# ╠═447e7c9e-74b1-11eb-27ea-71aa4338b11a
# ╟─c9dff6f4-7503-11eb-2715-0bf9d3ece9e1
# ╟─d834103c-7503-11eb-1a94-1fbad43801ff
# ╠═aa541288-74aa-11eb-1edc-ab6d7786f271
# ╠═c9dcac48-74aa-11eb-31a6-23357180c1c8
# ╟─30b1c1f0-7504-11eb-1be7-a9463caea809
# ╟─1fe70e38-751b-11eb-25b8-c741e1726613
# ╟─215291ec-74a2-11eb-3476-0dab43fd5a5e
# ╟─61db42c6-7505-11eb-1ddf-05e906234572
# ╟─cdd4cffc-74b1-11eb-1aa4-e333cb8601d1
# ╟─7489a570-74a3-11eb-1d0b-09d41604ffe1
# ╟─8a8e3f5e-74b2-11eb-3eed-e5468e573e45
# ╟─5864294a-74a5-11eb-23ef-f38a582f2c2d
# ╟─fa9c465e-74b2-11eb-2f3c-4be0e7f93bb5
# ╟─4fab4616-74b0-11eb-0088-6b50237d7d54
# ╟─275bf7ac-74b3-11eb-32c3-cda1e4f1f8c2
# ╟─537c54e4-74b3-11eb-341f-951b4a1e0b40
# ╟─c6e340ee-751e-11eb-3ca7-69595b3693b7
# ╠═54448d18-7528-11eb-209a-9717affa0d02
# ╠═acbc563a-7528-11eb-3c38-75a5b66c9241
# ╠═995392ee-752a-11eb-3394-0de331e24f40
# ╠═d22903d6-7529-11eb-2dcd-132cd27104c2
# ╟─844ed844-74b3-11eb-2ee1-2de664b26bc6
# ╟─4ffe927c-74b4-11eb-23a7-a18d7e51c75b
# ╟─91109e5c-74b3-11eb-1f31-c50e436bc6e0
# ╠═34109062-7525-11eb-10b3-d59d3a6dfda6
# ╟─9ab89a3a-7525-11eb-186d-29e4b61deb7f
# ╠═50034058-7525-11eb-345b-3334e71ac50e
# ╟─c0aec7ae-7505-11eb-2822-a151aad48fc9
# ╟─628aea22-7521-11eb-2edc-39ac62683aea
# ╠═99eeb11c-7524-11eb-2154-df7d84976445
# ╠═2ddcfb90-7520-11eb-2df7-172d07118b7e
# ╠═d62496b0-7524-11eb-3410-7177e7c7f8eb
# ╠═6aa8a76e-7524-11eb-22b5-015aab4191b0
# ╟─ee93eeb2-7524-11eb-342d-0343d8aebf59
# ╟─662d73b6-74b3-11eb-333d-f1323a001000
# ╠═d127303a-7521-11eb-3507-7341a416211f
# ╠═d4581b56-7522-11eb-2c15-991c0c790e67
# ╠═40c15c3a-7523-11eb-1f2a-bd90b127dad2
# ╠═08642690-7523-11eb-00dd-63d4cf6513dc
# ╠═deac4cf2-7523-11eb-2832-7b9d31389b08
# ╠═32887dfa-7524-11eb-35cd-051eff594fa9
# ╟─0f765670-7506-11eb-2a37-931b15bb387f
# ╟─82737d28-7507-11eb-1e39-c7dc12e18882
# ╟─40d538b2-7506-11eb-116b-efeb16b3478d
# ╟─df060a88-7507-11eb-034b-5346d67a0e0d
# ╟─60c8db60-7506-11eb-1468-c989809c933a
# ╟─8ed0be60-7506-11eb-2769-5f7da1c66243
# ╟─b9d636da-7506-11eb-37a6-3116d47b2787


text:
using Viruses, Test

ddf = 5.9
@test meltrate(-1, ddf)==0
@test meltrate(1, ddf)==ddf

# issue 1&2
@test meltrate(0, 1.0)==0


text:
@testset "fourier transform" begin
    𝐱 = rand(30, 40, 50, 6, 7) # where ch == 6 and batch == 7

    ft = FourierTransform((3, 4, 5))

    @test size(transform(ft, 𝐱)) == (30, 40, 50, 6, 7)
    @test size(truncate_modes(ft, transform(ft, 𝐱))) == (3, 4, 5, 6, 7)
    @test size(inverse(ft, truncate_modes(ft, transform(ft, 𝐱)))) == (3, 4, 5, 6, 7)
end


text:
function compute(file::String)::Int
    chars = Dict(Char(Int('A') + i) => i + 1 for i ∈ 0:25)
    return sum(i * chars[c] for (i, name) ∈ enumerate(sort(split(replace(file, "\"" => ""), ","))) for c ∈ name)
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
module ComradeAdaptMCMC

using AdaptiveMCMC
using Comrade
using TupleVectors
using AbstractMCMC

export AdaptMCMC

Base.@kwdef struct AdaptMCMC
    ntemp::Int
    swap::Symbol = :nonrev
    algorithm::Symbol = :ram
    fulladapt::Bool = true
    acc_sw::Float64 = 0.234
    all_levels::Bool = false
end

Comrade.samplertype(::Type{<:AdaptMCMC}) = Comrade.IsCube()

function AbstractMCMC.sample(post::Comrade.TransformedPosterior, sampler::AdaptMCMC, nsamples, burnin=nsamples÷2, args...; init_params=nothing, kwargs...)
    ℓ = logdensityof(post)
    function lpr(xx)
        for x in xx
            (x > 1.0 || x < 0.0) && return -Inf
        end
        return 0.0
    end

    θ0 = init_params
    if isnothing(init_params)
        @warn "No starting location chosen, picking start from random"
        θ0 = prior_sample(post)
    end

    apt = adaptive_rwm(θ0, ℓ, nsamples;
                       algorithm = sampler.algorithm,
                       b = burnin,
                       fulladapt=sampler.fulladapt,
                       L = sampler.ntemp,
                       acc_sw = sampler.acc_sw,
                       log_pr = lpr,
                       all_levels=sampler.all_levels,
                       swaps = sampler.swap,
                       kwargs...
                       )

    stats = (logl = apt.D, state = apt.R, accexp = apt.accRWM, accswp=apt.accSW)
    if sampler.all_levels
        chains = Tuple(TupleVector(transform.(Ref(post), eachcol(apt.allX[i]))) for i in eachindex(apt.allX))
    else
        chains = transform.(Ref(post), eachcol(apt.X)) |> TupleVector
    end
    return chains, stats
end

end


text:
"""
    getEnergySpectrum( CR_N::Vector{T}, CR_S::Vector{T},
                            CR_Cut::T, ρ::T;
                            par::CRMomentumDistributionConfig,
                            mode::Integer=3) where T

Compute the energy spectrum in ´[10^10 erg/g]´.
"""
function get_energy_spectrum(CR_N::Vector{T}, CR_S::Vector{T},
                            CR_Cut::T, ρ::T;
                            par::CRMomentumDistributionConfig,
                            mode::Integer = 3) where {T<:Real}

    # transform the norm dependent on IO mode
    norm = transform_norm(CR_N, mode)

    bounds = 10.0 .^ collect(log10(par.pmin):par.bin_width:log10(par.pmax))
    bound_up = Vector{Float64}(undef, par.Nbins)

    for i = 1:par.Nbins-1
        bound_up[i] = (bounds[i+1] > CR_Cut) ? CR_Cut : bounds[i+1]
    end

    bound_up[end] = (bounds[end] > CR_Cut) ? CR_Cut : bounds[end]

    bin_centers = Vector{Float64}(undef, par.Nbins)
    energy = Vector{Float64}(undef, par.Nbins)

    @inbounds for i = 1:par.Nbins
        bin_centers[i] = bounds[i] # 10.0^(0.5* (log10(bounds[i]) + log10(bound_up[i])))
        energy[i] = energy_integral(bounds[i], bound_up[i], norm[i], CR_S[i], ρ)

    end

    return CR_EnergySpectrum(bin_centers, energy)
end

"""
    getEnergySpectrum( spec::CRMomentumDistribution, ρ::T ) where T

Compute the energy spectrum in ´[10^10 erg/g]´.
"""
function getEnergySpectrum( spec::CRMomentumDistribution, ρ::T ) where T

    Nbins = Int64( 0.5* size(spec.norm,1) )
    bin_centers = Vector{Float64}(undef, Nbins)
    energy      = Vector{Float64}(undef, Nbins)

    k = 1
    @inbounds for i = 1:Nbins
        bin_centers[i] = spec.bound[k] #10.0^(0.5* (log10(spec.bound[k]) + log10(spec.bound[k+1])))
        slope          = log10( spec.bound[k] / spec.bound[k+1] ) /
                         log10( spec.norm[k]  / spec.norm[k+1] )
        energy[i]      = energy_integral(spec.bound[k], spec.bound[k+1], spec.norm[k], slope, ρ)
        k += 2
    end

    return CR_EnergySpectrum(bin_centers, energy)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/unsafe_utilities.jl
# Inbounds and non-allocating versions of simple Julia functions.
#############################################################################

"""
$(TYPEDSIGNATURES)

Performs `map!(f, y, x)` in an unsafe manner if y[i] is true, else no-op.
Assumes `n == length(x) == length(y)`. About 2x faster for small arrays (n < 1000).
"""
function unsafe_check_fill!(f, y::Vector{T}, x::T, n::Int) where T
    i = 1
    m = n + 1
    while i < m
        flag = f(@inbounds y[i])
        if flag
            @inbounds y[i] = x
        end
        i += 1
    end
    nothing
end

"""
Performs foreach and returns false if op returns false
"""
function foreach_until(op, x, init)
    if !init
        return false
    end
    for d in x
        if !op(d)
            return false
        end
    end
    return true
end


text:
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __julia-version__ = 1.7.2
# __author__        = "Max Kannenberg"
# __copyright__     = "2020-2022"
# __license__       = "ISC"

function createOutput(train::Train, settings::Settings, path::Path, movingSection::Dict, drivingCourse::Vector{Dict})
    if settings.outputDetail == :running_time
        output = movingSection[:t]  # TODO: or use drivingCourse[end][:t]

    elseif settings.outputDetail == :points_of_interest
        # add points of interest
        if !isempty(path.poi)

            # for elem in 1:length(driving_course)
            # end

            output = Dict[]
            POI = 1
            i = 1
            while POI <= length(path.poi) && i <= drivingCourse[end][:i]
                if path.poi[POI][:station] == drivingCourse[i][:s]
                    push!(output, drivingCourse[i])
                    POI = POI+1
                end
                i = i+1
            end
        end

    elseif settings.outputDetail == :driving_course
            output = drivingCourse

    elseif settings.outputDetail == :everything
        output = Dict{Symbol,Any}()
        merge!(output, Dict(:train => train, :path => path, :settings => settings))

        # add moving section and driving courses
        if settings[:operationModeMinimumRunningTime] == true
            merge!(output, Dict(:movingSectionMinimumRunningTime => movingSection,
                                    :drivingCourseMinimumRunningTime => drivingCourse))
        elseif settings[:operationModeMinimumEnergyConsumption] == true
            merge!(output, Dict(:movingSectionMinimumEnergyConsumption => movingSection,
                                    :drivingCourseMinimumEnergyConsumption => drivingCourse))
        end

        # add points of interest
        if !isempty(path.poi)
            pointsOfInterest = Vector{Dict}()
            POI = 1
            i = 1
            while POI <= length(path.poi) && i <= drivingCourse[end][:i]
                if path.poi[POI] == drivingCourse[i][:s]
                    push!(pointsOfInterest, drivingCourse[i])
                    POI = POI+1
                end
                i = i+1
            end

            if settings[:operationModeMinimumRunningTime] == true
                merge!(output, Dict(:pointsOfInterestMinimumRunningTime => pointsOfInterest))
            elseif settings[:operationModeMinimumEnergyConsumption] == true
                merge!(output, Dict(:pointsOfInterestMinimumEnergyConsumption => pointsOfInterest))
            end
        end
    else
        output = nothing
    end
    return output
end

#=
function createOutputDict(train::Train, settings::Settings, path::Path, movingSection::Dict, drivingCourse::Vector{Dict})
    outputDict = Dict{Symbol,Any}()
    merge!(outputDict, Dict(:train => train, :path => path, :settings => settings))

    # add moving section and driving courses
    if settings[:operationModeMinimumRunningTime] == true
        merge!(outputDict, Dict(:movingSectionMinimumRunningTime => movingSection,
                                :drivingCourseMinimumRunningTime => drivingCourse))
    elseif settings[:operationModeMinimumEnergyConsumption] == true
        merge!(outputDict, Dict(:movingSectionMinimumEnergyConsumption => movingSection,
                                :drivingCourseMinimumEnergyConsumption => drivingCourse))
    end

    # add points of interest
    if !isempty(path.poi)
        pointsOfInterest = Vector{Dict}()
        POI = 1
        i = 1
        while POI <= length(path.poi) && i <= drivingCourse[end][:i]
            if path.poi[POI] == drivingCourse[i][:s]
                push!(pointsOfInterest, drivingCourse[i])
                POI = POI+1
            end
            i = i+1
        end

        if settings[:operationModeMinimumRunningTime] == true
            merge!(outputDict, Dict(:pointsOfInterestMinimumRunningTime => pointsOfInterest))
        elseif settings[:operationModeMinimumEnergyConsumption] == true
            merge!(outputDict, Dict(:pointsOfInterestMinimumEnergyConsumption => pointsOfInterest))
        end
    end

    return outputDict
end # function createOutputDict
=#


text:
S!(metnet::MetNet, metider::IDER_TYPE, rxnider::IDER_TYPE, s::Real) =
    metnet.S[metindex(metnet, metider), rxnindex(metnet, rxnider)] = s

b!(metnet::MetNet, metider::IDER_TYPE, b::Real) = metnet.b[metindex(metnet, metider)] = b

ub!(metnet::MetNet, rxnider::IDER_TYPE, ub::Real) = metnet.ub[rxnindex(metnet, rxnider)] = Float64(ub)

lb!(metnet::MetNet, rxnider::IDER_TYPE, lb::Real) = metnet.lb[rxnindex(metnet, rxnider)] = Float64(lb)

bounds!(metnet::MetNet, ider::IDER_TYPE, lb::Real, ub::Real) =
    (idx = rxnindex(metnet, ider); metnet.lb[idx] = lb; metnet.ub[idx] = ub; nothing)


text:
function xbook(x)
  2x
end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
include("../src/LinearAlgebraicRepresentation.jl")

function generate_perpendicular_lines(steps::Int, minlen, maxlen)
    V = zeros(0,2)

    function rec(o, d, s)
        if s == 0 return end

        a = (maxlen-minlen)*rand() + minlen
        p = o + a*d
        V = [V; o; p]

        b = (a-minlen)*rand() + minlen
        p = o + b*d
        rec(p, d, s-1)

        b = (a-minlen)*rand() + minlen
        p = o + b*d
        rec(p, perpendicular(d), s-1)
    end

    function perpendicular(vec)
        v = zeros(size(vec))
        v[1] = vec[2]
        v[2] = vec[1]
        return v
    end

    rec([0 0], [1 0], steps)
    rec([0 0], [0 1], steps)
    vnum = size(V, 1)
    enum = vnum >> 1
    EV = spzeros(Int8, enum, vnum)
    for i in 1:enum
        EV[i, i*2-1:i*2] = 1
    end
    V, EV
end

function generate_random_lines(n, points_range, alphas_range)
    origins = points_range[1] + (points_range[2]-points_range[1])*rand(n, 2)
    directions = mapslices(normalize, rand(n, 2) - .5*ones(n, 2), dims=2)
    alphas = alphas_range[1] + (alphas_range[2]-alphas_range[1])*rand(n)
    new_points = Array{Float64, 2}(n, 2)
    for i in 1:n
        new_points[i, :] = origins[i, :] + alphas[i]*directions[i, :]
    end
    V = [origins; new_points]
    EV = spzeros(Int8, n, n*2)
    for i in 1:n
        EV[i, i] = 1
        EV[i, n+i] = 1
    end
    V, EV
end

function rubiks_example(ncubes = 3)
    V = Float64[
        0 0 0; 0 1 0;
        1 1 0; 1 0 0;
        0 0 1; 0 1 1;
        1 1 1; 1 0 1
    ]

    EV = sparse(Int8[
        -1  1  0  0  0  0  0  0;
        0 -1  1  0  0  0  0  0;
        0  0 -1  1  0  0  0  0;
        -1  0  0  1  0  0  0  0;
        -1  0  0  0  1  0  0  0;
        0 -1  0  0  0  1  0  0;
        0  0 -1  0  0  0  1  0;
        0  0  0 -1  0  0  0  1;
        0  0  0  0 -1  1  0  0;
        0  0  0  0  0 -1  1  0;
        0  0  0  0  0  0 -1  1;
        0  0  0  0 -1  0  0  1;
    ])

    FE = sparse(Int8[
        1  1  1 -1  0  0  0  0  0  0  0  0;
        0  0  0  0  0  0  0  0 -1 -1 -1  1;
        -1  0  0  0  1 -1  0  0  1  0  0  0;
        0 -1  0  0  0  1 -1  0  0  1  0  0;
        0  0 -1  0  0  0  1 -1  0  0  1  0;
        0  0  0  1 -1  0  0  1  0  0  0 -1;
    ])

    cube = [V, EV, FE]
    cubesRow = (zeros(0,3),spzeros(Int8,0,0),spzeros(Int8,0,0))

    for i in 1:ncubes
        cubesRow = Lar.skel_merge(cubesRow..., cube...)
        cube[1] = cube[1] + [zeros(8) zeros(8) ones(8)]
    end

    cubesRow = collect(cubesRow)
    cubesPlane = cubesRow
    num = size(cubesRow[1], 1)
    for i in 1:ncubes
        cubesPlane = Lar.skel_merge(cubesPlane..., cubesRow...)
        cubesRow[1] = cubesRow[1] + [zeros(num) ones(num) zeros(num)]
    end

    cubesPlane = collect(cubesPlane)
    cubesCube = cubesPlane
    num = size(cubesPlane[1], 1)
    for i in 1:ncubes
        cubesCube = Lar.skel_merge(cubesCube..., cubesPlane...)
        cubesPlane[1] = cubesPlane[1] + [ones(num) zeros(num) zeros(num)]
    end

    println("Arranging a cube of ", ncubes^3," cubes...")
    rubik = Lar.spatial_arrangement(cubesCube...)
    println("DONE")

    rubik = rubik[1] - 1.5, rubik[2:3]...
    c = cos(pi/6); s = sin(pi/6)
    M1 = [1  0 0; 0 c -s; 0 s c]
    M2 = [c -s 0; s c  0; 0 0 1]
    rot_rubik = rubik[1]*M1*M2, rubik[2:3]...

    println("Arranging two rubik cubes...")
    two_rubiks = Lar.skel_merge(rubik..., rot_rubik...)
    println("DONE")

    arranged_rubiks = Lar.spatial_arrangement(two_rubiks...)
end


text:
#!/usr/bin/env julia

using ColorSchemes, DataFrames, Luxor,  Colors

"""
change an RA value from 24 -> 0 to left edge -> right edge allowing for margin
"""

function rascale(n)
    rescale(n, 0, 24, (width/2) - margin, -(width/2) + margin)
end

"""
change a Declination value from -90 to 90
"""

function decscale(n)
    rescale(n, -90, 90, (height/2) - margin, (-height/2) + margin)
end

function magscale(n, highest_mag)
    # change mag value
    # lowest is about -2 (Sirius)
    # biggest circle 3, smallest 0.1
    rescale((2 + n), 0, highest_mag, 1.5, 0.01)
end

function drawgrid(theme)
    gsave()
    setline(theme[:gridline_width])
    sethue(theme[:gridline_color])
    fontsize(theme[:grid_label_font_size])
    fontface(theme[:grid_label_font_face])
    setopacity(theme[:grid_line_opacity])
    # vertical rules
    if  true # originally a test for rectangular or polar grids
        # rectangular
        # vertical lines for each RA value
        for ra in 0:24
            move(rascale(ra), -height/2 + margin)
            line(rascale(ra),  height/2 - margin)
            stroke()
            gsave()
            sethue(theme[:grid_label_font_color])
            setopacity(theme[:grid_label_font_opacity])
            text(string(ra), rascale(ra),  -(height/2) + 40)
            text(string(ra), rascale(ra),   (height/2) - 40)
            grestore()
        end
        # horizontal lines for each declination tick
        for dec in -90:10:90
            move(-width/2 + margin, decscale(dec))
            line(width/2 - margin, decscale(dec))
            stroke()
            gsave()
            sethue(theme[:grid_label_font_color])
            setopacity(theme[:grid_label_font_opacity])
            text(string(dec), width/2 - margin/2, decscale(dec))
            text(string(dec), -(width/2) + margin/2, decscale(dec))
            grestore()
        end
    end
    grestore()
end

function drawstars(table, theme)
  global starcounter
  for star in eachrow(table)
    starcounter += 1
    proper, ra, dec, mag, ci = map(last, star)
    if mag < -2 # don't draw the Sun ! :)
        continue
    end
    if !isna(ci)
        # set color of star
        #=
        1.41	M0	Red
        0.82	K0	Orange
        0.59	G0	Yellow
        0.31	F0	Yellowish
        0.00	A0	White
        -0.29	B0	Blue
        =#
        # rescale: blue is 0, red is 1
        sethue(colorscheme(theme[:starcolorscheme], rescale(ci, -0.3, 1.5, 0, 1)))
    else
        sethue("grey50")
    end
    circle(rascale(ra), decscale(dec), magscale(mag, highest_mag), :fill)
    # text labels for stars if mag is big enough
    fontface(theme[:starname_font_face])
    fontsize(theme[:starname_font_size])
    if isna(proper)
        # this star has no name
    else
        if mag <=  theme[:star_min_mag_label_size]
            fontsize(theme[:starname_font_size] - mag/10)
            sethue(theme[:starname_font_color])
            rect(rascale(ra), decscale(dec), 0.2, -5, :fill)
            text(proper, rascale(ra), decscale(dec) - 5)
        end
    end
  end
end

function drawconstellation(constell, theme)
  sethue(theme[:constellation_boundary_stroke_color])
  setline(theme[:constellation_boundary_stroke_width])
  v = constellation_boundaries[constell]
  realconstellationname = first(constellation_names[constell])
  gsave()
  setdash(theme[:constellation_boundary_stroke_pattern])
  for pgon in v
      pgon_scaled = map(pt -> Point(rascale((pt.x)/3600), decscale((pt.y)/3600)), pgon)
      pgon_scaled = map(pt -> Point(round(Int, pt.x), round(Int, pt.y)),  pgon_scaled)
      poly(pgon_scaled, close=true, :stroke)
      sethue(theme[:constellation_name_font_color])
      poly(pgon_scaled, close=true, :stroke)
      centrep = polycentroid(pgon_scaled)
      setopacity(theme[:constellation_name_opacity])
      fontsize(theme[:constellation_name_font_size])
      fontface(theme[:constellation_name_font_face])
      sethue(theme[:constellation_name_font_color])
      # draw text
      textcentred(string(realconstellationname), centrep)
  end
  grestore()
end

function main(d, w, h, margin, highest_mag)
  d1 = d[d[:mag] .< highest_mag, [:proper, :ra, :dec, :mag, :ci]]

  # all rows which have constellations
  # d[complete_cases(d[:,[:con]]), :]
  ## d[complete_cases(d[:,[:con, :proper]]), :]

  # a theme controls all the fills and strokes
  # this is the default theme
  default_theme = Dict(
      :name                                   => "bright stars on dark background",
      :highest_magnitude                      => 15,
      :starcolorscheme                        => loadcolorscheme("thermometercolors"),
      :starname_font_face                     => "Times-Italic",
      :starname_font_size                     => 12,
      :starname_font_color                    => "grey90",
      :star_min_mag_label_size                => 6,
      :background_color                       => "black",
      :constellation_name_opacity             => 0.8,
      :constellation_boundary_stroke_color    => "blue",
      :constellation_boundary_fill_color      => "",
      :constellation_boundary_stroke_pattern  => "solid",
      :constellation_boundary_stroke_width    => 0.2,
      :constellation_name_font_face           => "BodoniMT-UltraBold",
      :constellation_name_font_size           => 26,
      :constellation_name_font_color          => "grey90",
      :gridline_color                         => "red",
      :gridline_width                         => .2,
      :gridline_pattern                       => "dot",
      :grid_label_font_face                   => "BodoniMT-UltraBold",
      :grid_label_font_size                   => 12,
      :grid_label_font_color                  => "grey90",
      :grid_label_font_opacity                => 1.0,
      :grid_line_opacity                      => 0.3
  )

  antique_theme = Dict(
      :name                                   => "antique theme",
      :highest_magnitude                      => 10,
      :starcolorscheme                        => loadcolorscheme("thermometercolors"),
      :starname_font_face                     => "InputMono-Medium",
      :starname_font_size                     => 9,
      :starname_font_color                    => "grey30",
      :star_min_mag_label_size                => 10,
      :background_color                       => "ivory",
      :constellation_name_opacity             => 0.8,
      :constellation_boundary_stroke_color    => "orange",
      :constellation_boundary_fill_color      => "",
      :constellation_boundary_stroke_pattern  => "solid",
      :constellation_boundary_stroke_width    => 0.4,
      :constellation_name_font_face           => "BodoniMT-UltraBold",
      :constellation_name_font_size           => 25,
      :constellation_name_font_color          => "Grey70",
      :gridline_color                         => "blue",
      :gridline_width                         => .2,
      :gridline_pattern                       => "dotted",
      :grid_label_font_face                   => "BodoniMT-UltraBold",
      :grid_label_font_size                   => 14,
      :grid_label_font_color                  => "grey40",
      :grid_label_font_opacity                => 1.0,
      :grid_line_opacity                      => 0.2
  )

  Drawing(width, height, "/tmp/starmap.pdf")
  origin()
  global starcounter = 0
  theme=default_theme
  background(theme[:background_color])
  drawgrid(theme)
  x = -width/2 + 350
  y = -height/2 + 350
  for k in collect(keys(constellation_boundaries))[1:end]
      gsave()
      drawconstellation(k, theme)
      grestore()
  end

  drawstars(d1, theme)

  finish()
  preview()
end

# constallation boundaries
# some constellations are in two parts
# coordinates are in RA seconds float and Declination seconds float

info("loading constellations")
tic()
include("constellations.jl")
toc()

info("loading star data")

# for convenience I copied this file from https://github.com/astronexus/HYG-Database

astrodbfile = "docs/examples/hygdata_v3.csv"

tic()
db = readtable(astrodbfile, header = true)
toc()

width, height = 4000, 2000
margin = 80
highest_mag = 15

info("drawing chart")
tic()
main(db, width, height, margin, highest_mag)
toc()
@show starcounter