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

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

# ╔═╡ 5719ba80-59c7-11eb-09d7-0b79c4fe0b8f
using Plots 	# Plotting backend

# ╔═╡ 2e943e40-59d7-11eb-1b9c-b5ec5d9cf114
html"
Projectile Motion with Euler's Method
"

# ╔═╡ ae334720-59de-11eb-2bb5-81773ca0f2e4
html"
Setting things up for plotting
"

# ╔═╡ 40e2d3e0-59d7-11eb-0a7b-9f108787b759
html"
A struct to store the quantities we care about:
"

# ╔═╡ 89eaa350-59c4-11eb-0187-c5ad24105f84
@Base.kwdef mutable struct Ball
	θ::Float64 				# Angle of launch above neutral, radians
	x::Float64 = 0 			# Initial (x, y) position is (0, 0), in meters
	y::Float64 = 0
	vx::Float64 = 30cos(θ) 	# 30 m/s ≈ 70 mph
	vy::Float64 = 30sin(θ)
	ax::Float64 = 0 		# Ignoring wind resistance, ax is always zero
	ay::Float64 = -9.8 		# Gravity, m/s²
	distance::Float64 = Inf # Distance travelled before landing, recorded later
end

# ╔═╡ 84ce10ae-59d7-11eb-0159-d940a49ea8f3
html"
Time step:
"

# ╔═╡ ccf7c2e0-59c9-11eb-2c30-05b6926c7125
τ = .001 # time step for Euler's method, seconds

# ╔═╡ 9695ea20-59d7-11eb-0295-73f37d0a1922
html"
Euler's Method applied to Ball struct as step function:
"

# ╔═╡ c876563e-59c5-11eb-0960-fd0d317f137e
function step!(b::Ball)
	b.vx += b.ax*τ
	b.vy += b.ay*τ
	b.x += b.vx*τ
	b.y += b.vy*τ
	if b.y <= 0
		b.y = -b.y # let's assume ball bounces perfectly elastically, for fun
		b.vy = -b.vy
		if b.distance == Inf
			b.distance = b.x
		end
	end
end

# ╔═╡ 049cdc20-59df-11eb-0219-0194a1d3ce10
html"
Setting up empty arrays for plotting 3 different balls:
"

# ╔═╡ fd3d4910-59ca-11eb-0533-abe0fa401cc1

begin
	maxsteps::Int64 = 6000
	maxplots::Int64 = 120
	plotstep::Int64 = maxsteps ÷ maxplots
	iplot::Int64 = 1

	xplot = Array{Float64, 2}(undef, 3, maxplots)
	yplot = Array{Float64, 2}(undef, 3, maxplots)

	b1 = Ball(θ=π/4)
	b2 = Ball(θ=π/3)
	b3 = Ball(θ=π/6)
end

# ╔═╡ c5404f30-59de-11eb-1ea5-d5e0a3378e35
html"
Plotting
"

# ╔═╡ e8675532-59de-11eb-20be-a94a386bfd21
html"
Calling main step function:
"

# ╔═╡ c65c4ee2-59d0-11eb-01cd-af3ec58fd3b8
for i ∈ 1:maxsteps
	step!(b1)
	step!(b2)
	step!(b3)
	if i%plotstep == 1
		xplot[:, iplot] = [b1.x, b2.x, b3.x]
		yplot[:, iplot] = [b1.y, b2.y, b3.y]
		iplot += 1
	end
end

# ╔═╡ 2d96f8e0-59df-11eb-118b-a7b8e259d543
html"
Distance travelled before landing for each ball:
"

# ╔═╡ 25419320-59d6-11eb-2d88-19eb6787edbf
b1.distance, b2.distance, b3.distance

# ╔═╡ 51237770-59df-11eb-312f-79544f3b168a
html"
This is what we expect; we launched the first ball at π/4 radians, and the other two equidistant from it:
"

# ╔═╡ 55c72500-59d6-11eb-1274-8b6dfb43cfb9
b1.θ, rad2deg(b1.θ) # first ball launched at 45°

# ╔═╡ eca1d520-59df-11eb-2dc4-d9d2779b6bc3
b2.θ, rad2deg(b2.θ) # launched at 45° + 15°

# ╔═╡ 2863fd40-59e0-11eb-16b5-991137d7867e
b3.θ, rad2deg(b3.θ) # launched at 45° - 15°

# ╔═╡ 1779f3f0-59d0-11eb-0b0d-2503afafab7d
begin
	scatter(xplot[1, :], yplot[1, :], xlim=(0, 100), ylim=(0, 40),
		label="θ=$(round(rad2deg(b1.θ), digits=3))°")
	scatter!(xplot[2, :], yplot[2, :],
		label="θ=$(round(rad2deg(b2.θ), digits=3))°")
	scatter!(xplot[3, :], yplot[3, :],
		label="θ=$(round(rad2deg(b3.θ), digits=3))°")
	scatter!(title="Projectile Motion Over $(τ*maxsteps) Seconds",
		xlabel="Distance [m]", ylabel="Height [m]")
end

# ╔═╡ a34de050-59d8-11eb-0125-553fbef90785
@gif for i ∈ 1:maxplots
	scatter([xplot[1, i]], [yplot[1, i]], ylim=(0, 40), xlim=(0, 100),
	    label="θ=$(round(rad2deg(b1.θ), digits=3))°", markersize=10)
	scatter!([xplot[2, i]], [yplot[2, i]], markersize=10,
		label="θ=$(round(rad2deg(b2.θ), digits=3))°")
	scatter!([xplot[3, i]], [yplot[3, i]], markersize=10,
	    label="θ=$(round(rad2deg(b3.θ), digits=3))°")
	scatter!(title="Projectile Motion Over $(τ*maxsteps) Seconds",
		xlabel="Distance [m]", ylabel="Height [m]")
end

# ╔═╡ 698b9b50-59e2-11eb-2686-13b813de6970
html"
Try it yourself!
"

# ╔═╡ 6b0f3142-59e6-11eb-3450-b30c9151ff64
html"This currently only works in the Julia REPL, not on static webpages!"

# ╔═╡ 74e99c40-59e2-11eb-3ae6-6d9363a2834a
html"
Choose an angle of launch, α, in degrees:
"

# ╔═╡ 1d761150-59e2-11eb-00b4-3b97e23069eb
@bind α html"<input type='range' min='1' max='89' step='1' value='42'>"

# ╔═╡ 9addf310-59e2-11eb-06f7-2d5cce4a9f73
begin
	steps::Int64 = 6400
	plots::Int64 = 200
	sstep::Int64 = steps÷plots

	ball = Ball(θ=deg2rad(α))
	x = Array{Float64, 1}(undef, plots)
	y = Array{Float64, 1}(undef, plots)
	jplot::Int64 = 1
	for i ∈ 1:steps
		step!(ball)
		if i%sstep == 1
			x[jplot] = ball.x
			y[jplot] = ball.y
			jplot += 1
		end
	end
end

# ╔═╡ c87d5f80-59e3-11eb-2fb2-dbaca4ee07f0
begin
	plot(x, y, lw=5, color=:green, ylim=(0, 50), xlim=(0, 100),
		title="Initial Launch at α=$(α)°", xlabel="Distance [m]", ylabel="Height [m]",
		label="")
	scatter!([ball.distance], [.5], markersize=5,
		label="Flight distance: $(round(ball.distance, digits=2)) meters")
end

# ╔═╡ Cell order:
# ╟─2e943e40-59d7-11eb-1b9c-b5ec5d9cf114
# ╟─ae334720-59de-11eb-2bb5-81773ca0f2e4
# ╟─40e2d3e0-59d7-11eb-0a7b-9f108787b759
# ╠═89eaa350-59c4-11eb-0187-c5ad24105f84
# ╟─84ce10ae-59d7-11eb-0159-d940a49ea8f3
# ╠═ccf7c2e0-59c9-11eb-2c30-05b6926c7125
# ╟─9695ea20-59d7-11eb-0295-73f37d0a1922
# ╠═c876563e-59c5-11eb-0960-fd0d317f137e
# ╟─049cdc20-59df-11eb-0219-0194a1d3ce10
# ╠═fd3d4910-59ca-11eb-0533-abe0fa401cc1
# ╟─c5404f30-59de-11eb-1ea5-d5e0a3378e35
# ╠═5719ba80-59c7-11eb-09d7-0b79c4fe0b8f
# ╟─e8675532-59de-11eb-20be-a94a386bfd21
# ╠═c65c4ee2-59d0-11eb-01cd-af3ec58fd3b8
# ╟─2d96f8e0-59df-11eb-118b-a7b8e259d543
# ╠═25419320-59d6-11eb-2d88-19eb6787edbf
# ╟─51237770-59df-11eb-312f-79544f3b168a
# ╠═55c72500-59d6-11eb-1274-8b6dfb43cfb9
# ╠═eca1d520-59df-11eb-2dc4-d9d2779b6bc3
# ╠═2863fd40-59e0-11eb-16b5-991137d7867e
# ╠═1779f3f0-59d0-11eb-0b0d-2503afafab7d
# ╠═a34de050-59d8-11eb-0125-553fbef90785
# ╟─698b9b50-59e2-11eb-2686-13b813de6970
# ╟─6b0f3142-59e6-11eb-3450-b30c9151ff64
# ╟─74e99c40-59e2-11eb-3ae6-6d9363a2834a
# ╟─1d761150-59e2-11eb-00b4-3b97e23069eb
# ╟─9addf310-59e2-11eb-06f7-2d5cce4a9f73
# ╠═c87d5f80-59e3-11eb-2fb2-dbaca4ee07f0

max_stars_repo_name: benharmonics/GeneralPhysics

max_issues_repo_name: benharmonics/GeneralPhysics

max_forks_repo_name: benharmonics/GeneralPhysics


text:
extracttriangle(a::Triangle) = (a.data...,)
function triangleview(::Type{Triangle}, A::AbstractArray)
    if size(img, 1) != 3
        @error "Selected dimension, $i, must be of size 3"
    else
        return mappedarray(Triangle, extracttriangle,
                           A[expand2indices(A,(1,))...],
                           A[expand2indices(A,(2,))...],
                           A[expand2indices(A,(3,))...])
    end
end

ext: jl


text:
# This contains several objectives

##
# TODO make objectives have ! names
##
# Objective: L∞ norm of the disturbance
function max_disturbance!(model::Model, var)
    o = symbolic_infty_norm(var)
    @objective(model, Min, o)
    return o
end

function min_sum!(model::Model, var)
    o = sum(sum.(var))
    @objective(model, Min, o)
    return o
end

function max_sum!(model::Model, var)
    o = sum(sum.(var))
    @objective(model, Max, o)
    return o
end

# This is the default when creating a model. Only used for explicit-ness.
feasibility_problem!(model::Model) = nothing

function linear_objective!(mode::Model, map::HPolytope, var)
    c, d = tosimplehrep(map)
    o = c * var - d
    @objective(model, Min, o)
    return o
end


text:
__precompile__()

"""

A module for solving Cauchy problems for linear-stiff systems of two ordinary differential
equations.

See `Problem`, `Options`, and `solve` for more info.

"""
baremodule A2

using Base

# Load the internal logic
Base.include(A2, "Internal/Internal.jl")

# Export the second-level components
using .Internal

# Export the first-level components
export Problem, Options, solve

end


text:
using Test
using Recommenders: σ, Logloss, grad

@testset "Sigmoid." begin
    @test σ(0.0) == 0.5
    @test σ(100) ≈ 1.0
    @test σ(-100) < 1e-20
    @test σ(1.0) ≈ 1.0 / (1 + exp(-1))
    @test σ(0.1) > 0.5
    @test σ(-0.1) < 0.5
end

@testset "Logloss." begin
    logloss = Logloss()
    @test logloss(100, 1) ≈ 0
    @test logloss(-100, 1) > 20
    @test logloss(3, 0) > 3
    @test logloss(-100, 0) ≈ 0
    @test logloss(0.0, 1) == logloss(0.0, 0)
end

@testset "Grad of Logloss." begin
    logloss = Logloss()

    @test grad(logloss, 0.3243, 0) ≈ σ(0.3243)
    @test grad(logloss, 0.3243, 1) ≈ σ(0.3243) - 1

    # check sign
    logit = 0.4
    label = 1
    loss = logloss(logit, label)
    grad_logloss = grad(logloss, logit, label)
    @test logloss(logit - 0.1 * grad_logloss, label) < loss

    label = 0
    loss = logloss(logit, label)
    grad_logloss = grad(logloss, logit, label)
    @test logloss(logit - 0.1 * grad_logloss, label) < loss

    logit = -0.3
    label = 1
    loss = logloss(logit, label)
    grad_logloss = grad(logloss, logit, label)
    @test logloss(logit - 0.1 * grad_logloss, label) < loss

    label = 0
    loss = logloss(logit, label)
    grad_logloss = grad(logloss, logit, label)
    @test logloss(logit - 0.1 * grad_logloss, label) < loss
end


text:
"""
Provides a collection of utility functions and types that are used in other submodules.
"""
module Utilities

using Base.Meta
import Base: isdeprecated, Docs.Binding
using DocStringExtensions
import Markdown, LibGit2
import Base64: stringmime

# escape characters that has a meaning in regex
regex_escape(str) = sprint(escape_string, str, "\\^\$.|?*+()[{")

# helper to display linerange for error printing
function find_block_in_file(code, file)
    source_file = Base.find_source_file(file)
    source_file === nothing && return nothing
    isfile(source_file) || return nothing
    content = read(source_file, String)
    content = replace(content, "\r\n" => "\n")
    # make a regex of the code that matches leading whitespace
    rcode = "\\h*" * replace(regex_escape(code), "\\n" => "\\n\\h*")
    blockidx = findfirst(Regex(rcode), content)
    blockidx === nothing && return nothing
    startline = countlines(IOBuffer(content[1:prevind(content, first(blockidx))]))
    endline = startline + countlines(IOBuffer(code)) + 1 # +1 to include the closing ```
    return ":$(startline)-$(endline)"
end

# Pretty-printing locations
function locrepr(file, line=nothing)
    str = Base.contractuser(file) # TODO: Maybe print this relative the doc-root??
    line !== nothing && (str = str * "$(line)")
    return str
end

# Directory paths.

"""
Returns the current directory.
"""
function currentdir()
    d = Base.source_dir()
    d === nothing ? pwd() : d
end

"""
Returns the path to the Documenter `assets` directory.
"""
assetsdir() = normpath(joinpath(dirname(@__FILE__), "..", "..", "assets"))

cleandir(d::AbstractString) = (isdir(d) && rm(d, recursive = true); mkdir(d))

"""
Find the path of a file relative to the `source` directory. `root` is the path
to the directory containing the file `file`.

It is meant to be used with `walkdir(source)`.
"""
srcpath(source, root, file) = normpath(joinpath(relpath(root, source), file))

# Slugify text.

"""
Slugify a string into a suitable URL.
"""
function slugify(s::AbstractString)
    s = replace(s, r"\s+" => "-")
    s = replace(s, r"^\d+" => "")
    s = replace(s, r"&" => "-and-")
    s = replace(s, r"[^\p{L}\p{P}\d\-]+" => "")
    s = strip(replace(s, r"\-\-+" => "-"), '-')
end
slugify(object) = string(object) # Non-string slugifying doesn't do anything.

# Parse code blocks.

"""
Returns a vector of parsed expressions and their corresponding raw strings.

Returns a `Vector` of tuples `(expr, code)`, where `expr` is the corresponding expression
(e.g. a `Expr` or `Symbol` object) and `code` is the string of code the expression was
parsed from.

The keyword argument `skip = N` drops the leading `N` lines from the input string.

If `raise=false` is passed, the `Meta.parse` does not raise an exception on parse errors,
but instead returns an expression that will raise an error when evaluated. `parseblock`
returns this expression normally and it must be handled appropriately by the caller.
"""
function parseblock(code::AbstractString, doc, file; skip = 0, keywords = true, raise=true)
    # Drop `skip` leading lines from the code block. Needed for deprecated `{docs}` syntax.
    code = string(code, '\n')
    code = last(split(code, '\n', limit = skip + 1))
    endofstr = lastindex(code)
    results = []
    cursor = 1
    while cursor < endofstr
        # Check for keywords first since they will throw parse errors if we `parse` them.
        line = match(r"^(.*)\r?\n"m, SubString(code, cursor)).match
        keyword = Symbol(strip(line))
        (ex, ncursor) =
            # TODO: On 0.7 Symbol("") is in Docs.keywords, remove that check when dropping 0.6
            if keywords && (haskey(Docs.keywords, keyword) || keyword == Symbol(""))
                (QuoteNode(keyword), cursor + lastindex(line))
            else
                try
                    Meta.parse(code, cursor; raise=raise)
                catch err
                    push!(doc.internal.errors, :parse_error)
                    @warn "failed to parse exception in $(Utilities.locrepr(file))" exception = err
                    break
                end
            end
        str = SubString(code, cursor, prevind(code, ncursor))
        if !isempty(strip(str))
            push!(results, (ex, str))
        end
        cursor = ncursor
    end
    results
end
isassign(x) = isexpr(x, :(=), 2) && isa(x.args[1], Symbol)

# Checking arguments.

"""
Prints a formatted warning to the user listing unrecognised keyword arguments.
"""
function check_kwargs(kws)
    isempty(kws) && return
    out = IOBuffer()
    println(out, "Unknown keywords:\n")
    for (k, v) in kws
        println(out, "  ", k, " = ", v)
    end
    @warn(String(take!(out)))
end

# Finding submodules.

const ModVec = Union{Module, Vector{Module}}

"""
Returns the set of submodules of a given root module/s.
"""
function submodules(modules::Vector{Module})
    out = Set{Module}()
    for each in modules
        submodules(each, out)
    end
    out
end
function submodules(root::Module, seen = Set{Module}())
    push!(seen, root)
    for name in names(root, all=true)
        if Base.isidentifier(name) && isdefined(root, name) && !isdeprecated(root, name)
            object = getfield(root, name)
            if isa(object, Module) && !(object in seen) && parentmodule(object::Module) == root
                submodules(object, seen)
            end
        end
    end
    return seen
end

## objects
## =======

"""
Represents an object stored in the docsystem by its binding and signature.
"""
struct Object
    binding   :: Binding
    signature :: Type

    function Object(b::Binding, signature::Type)
        m = nameof(b.mod) === b.var ? parentmodule(b.mod) : b.mod
        new(Binding(m, b.var), signature)
    end
end

function splitexpr(x::Expr)
    isexpr(x, :macrocall) ? splitexpr(x.args[1]) :
    isexpr(x, :.)         ? (x.args[1], x.args[2]) :
    error("Invalid @var syntax `$x`.")
end
splitexpr(s::Symbol) = :(Main), quot(s)
splitexpr(other)     = error("Invalid @var syntax `$other`.")

"""
    object(ex, str)

Returns a expression that, when evaluated, returns an [`Object`](@ref) representing `ex`.
"""
function object(ex::Union{Symbol, Expr}, str::AbstractString)
    binding   = Expr(:call, Binding, splitexpr(Docs.namify(ex))...)
    signature = Base.Docs.signature(ex)
    isexpr(ex, :macrocall, 2) && !endswith(str, "()") && (signature = :(Union{}))
    Expr(:call, Object, binding, signature)
end

function object(qn::QuoteNode, str::AbstractString)
    if haskey(Base.Docs.keywords, qn.value)
        binding = Expr(:call, Binding, Main, qn)
        Expr(:call, Object, binding, Union{})
    else
        error("'$(qn.value)' is not a documented keyword.")
    end
end

function Base.print(io::IO, obj::Object)
    print(io, obj.binding)
    print_signature(io, obj.signature)
end
print_signature(io::IO, signature::Union{Union, Type{Union{}}}) = nothing
print_signature(io::IO, signature)        = print(io, '-', signature)

## docs
## ====

"""
    docs(ex, str)

Returns an expression that, when evaluated, returns the docstrings associated with `ex`.
"""
function docs end

# Macro representation changed between 0.4 and 0.5.
function docs(ex::Union{Symbol, Expr}, str::AbstractString)
    isexpr(ex, :macrocall, 2) && !endswith(rstrip(str), "()") && (ex = quot(ex))
    :(Base.Docs.@doc $ex)
end
docs(qn::QuoteNode, str::AbstractString) = :(Base.Docs.@doc $(qn.value))

"""
Returns the category name of the provided [`Object`](@ref).
"""
doccat(obj::Object) = startswith(string(obj.binding.var), '@') ?
    "Macro" : doccat(obj.binding, obj.signature)

function doccat(b::Binding, ::Union{Union, Type{Union{}}})
    if b.mod === Main && haskey(Base.Docs.keywords, b.var)
        "Keyword"
    elseif startswith(string(b.var), '@')
        "Macro"
    else
        doccat(getfield(b.mod, b.var))
    end
end

doccat(b::Binding, ::Type)  = "Method"

doccat(::Function) = "Function"
doccat(::DataType) = "Type"
doccat(x::UnionAll) = doccat(Base.unwrap_unionall(x))
doccat(::Module)   = "Module"
doccat(::Any)      = "Constant"

"""
    filterdocs(doc, modules)

Remove docstrings from the markdown object, `doc`, that are not from one of `modules`.
"""
function filterdocs(doc::Markdown.MD, modules::Set{Module})
    if isempty(modules)
        # When no modules are specified in `makedocs` then don't filter anything.
        doc
    else
        if haskey(doc.meta, :module)
            doc.meta[:module] ∈ modules ? doc : nothing
        else
            if haskey(doc.meta, :results)
                out = []
                results = []
                for (each, result) in zip(doc.content, doc.meta[:results])
                    r = filterdocs(each, modules)
                    if r !== nothing
                        push!(out, r)
                        push!(results, result)
                    end
                end
                if isempty(out)
                    nothing
                else
                    md = Markdown.MD(out)
                    md.meta[:results] = results
                    md
                end
            else
                out = []
                for each in doc.content
                    r = filterdocs(each, modules)
                    r === nothing || push!(out, r)
                end
                isempty(out) ? nothing : Markdown.MD(out)
            end
        end
    end
end
# Non-markdown docs won't have a `.meta` field so always just accept those.
filterdocs(other, modules::Set{Module}) = other

"""
Does the given docstring represent actual documentation or a no docs error message?
"""
nodocs(x) = occursin("No documentation found.", stringmime("text/plain", x))
nodocs(::Nothing) = false

header_level(::Markdown.Header{N}) where {N} = N

"""
    repo_root(file; dbdir=".git")

Tries to determine the root directory of the repository containing `file`. If the file is
not in a repository, the function returns `nothing`.

The `dbdir` keyword argument specifies the name of the directory we are searching for to
determine if this is a repostory or not. If there is a file called `dbdir`, then it's
contents is checked under the assumption that it is a Git worktree or a submodule.
"""
function repo_root(file; dbdir=".git")
    parent_dir, parent_dir_last = dirname(abspath(file)), ""
    while parent_dir != parent_dir_last
        dbdir_path = joinpath(parent_dir, dbdir)
        isdir(dbdir_path) && return parent_dir
        # Let's see if this is a worktree checkout
        if isfile(dbdir_path)
            contents = chomp(read(dbdir_path, String))
            if startswith(contents, "gitdir: ")
                if isdir(joinpath(parent_dir, contents[9:end]))
                    return parent_dir
                end
            end
        end
        parent_dir, parent_dir_last = dirname(parent_dir), parent_dir
    end
    return nothing
end

"""
    $(SIGNATURES)

Returns the path of `file`, relative to the root of the Git repository, or `nothing` if the
file is not in a Git repository.
"""
function relpath_from_repo_root(file)
    cd(dirname(file)) do
        root = repo_root(file)
        root !== nothing && startswith(file, root) ? relpath(file, root) : nothing
    end
end

function repo_commit(file)
    cd(dirname(file)) do
        readchomp(`git rev-parse HEAD`)
    end
end

function url(repo, file; commit=nothing)
    file = abspath(file)
    if !isfile(file)
        @warn "couldn't find file \"$file\" when generating URL"
        return nothing
    end
    file = realpath(file)
    remote = getremote(dirname(file))
    isempty(repo) && (repo = "https://github.com/$remote/blob/{commit}{path}")
    path = relpath_from_repo_root(file)
    if path === nothing
        nothing
    else
        repo = replace(repo, "{commit}" => commit === nothing ? repo_commit(file) : commit)
        # Note: replacing any backslashes in path (e.g. if building the docs on Windows)
        repo = replace(repo, "{path}" => string("/", replace(path, '\\' => '/')))
        repo = replace(repo, "{line}" => "")
        repo
    end
end

url(remote, repo, doc) = url(remote, repo, doc.data[:module], doc.data[:path], linerange(doc))

function url(remote, repo, mod, file, linerange)
    file === nothing && return nothing # needed on julia v0.6, see #689
    remote = getremote(dirname(file))
    isabspath(file) && isempty(remote) && isempty(repo) && return nothing

    # make sure we get the true path, as otherwise we will get different paths when we compute `root` below
    if isfile(file)
        file = realpath(abspath(file))
    end

    # Format the line range.
    line = format_line(linerange, LineRangeFormatting(repo_host_from_url(repo)))
    # Macro-generated methods such as those produced by `@deprecate` list their file as
    # `deprecated.jl` since that is where the macro is defined. Use that to help
    # determine the correct URL.
    if inbase(mod) || !isabspath(file)
        file = replace(file, '\\' => '/')
        base = "https://github.com/JuliaLang/julia/blob"
        dest = "base/$file#$line"
        if isempty(Base.GIT_VERSION_INFO.commit)
            "$base/v$VERSION/$dest"
        else
            commit = Base.GIT_VERSION_INFO.commit
            "$base/$commit/$dest"
        end
    else
        path = relpath_from_repo_root(file)
        if isempty(repo)
            repo = "https://github.com/$remote/blob/{commit}{path}#{line}"
        end
        if path === nothing
            nothing
        else
            repo = replace(repo, "{commit}" => repo_commit(file))
            # Note: replacing any backslashes in path (e.g. if building the docs on Windows)
            repo = replace(repo, "{path}" => string("/", replace(path, '\\' => '/')))
            repo = replace(repo, "{line}" => line)
            repo
        end
    end
end

function getremote(dir::AbstractString)
    remote =
        try
            cd(() -> readchomp(`git config --get remote.origin.url`), dir)
        catch err
            ""
        end
    m = match(LibGit2.GITHUB_REGEX, remote)
    if m === nothing
        travis = get(ENV, "TRAVIS_REPO_SLUG", "")
        isempty(travis) ? "" : travis
    else
        m[1]
    end
end

"""
$(SIGNATURES)

Returns the first 5 characters of the current git commit hash of the directory `dir`.
"""
function get_commit_short(dir)
    commit = cd(dir) do
        readchomp(`git rev-parse HEAD`)
    end
    (length(commit) > 5) ? commit[1:5] : commit
end

function inbase(m::Module)
    if m ≡ Base
        true
    else
        parent = parentmodule(m)
        parent ≡ m ? false : inbase(parent)
    end
end

# Repository hosts
#   RepoUnknown denotes that the repository type could not be determined automatically
@enum RepoHost RepoGithub RepoBitbucket RepoGitlab RepoUnknown

# Repository host from repository url
# i.e. "https://github.com/something" => RepoGithub
#      "https://bitbucket.org/xxx" => RepoBitbucket
# If no match, returns RepoUnknown
function repo_host_from_url(repoURL::String)
    if occursin("bitbucket", repoURL)
        return RepoBitbucket
    elseif occursin("github", repoURL) || isempty(repoURL)
        return RepoGithub
    elseif occursin("gitlab", repoURL)
        return RepoGitlab
    else
        return RepoUnknown
    end
end

# Find line numbers.
# ------------------

linerange(doc) = linerange(doc.text, doc.data[:linenumber])

function linerange(text, from)
    # text is assumed to be a Core.SimpleVector (svec) from the .text field of a Docs.DocStr object.
    # Hence, we need to be careful when summing over an empty svec below.
    #
    # Also, the isodd logic _appears_ to be there to handle variable interpolation into docstrings. In that case,
    # the .text field seems to become longer than just 1 element and every even element is the interpolated object,
    # and only the odd ones actually contain the docstring text as a string.
    lines = sum(Int[isodd(n) ? newlines(s) : 0 for (n, s) in enumerate(text)])
    return lines > 0 ? (from:(from + lines + 1)) : (from:from)
end

struct LineRangeFormatting
    prefix::String
    separator::String

    function LineRangeFormatting(host::RepoHost)
        if host == RepoBitbucket
            new("", ":")
        elseif host == RepoGitlab
            new("L", "-")
        else
            # default is github-style
            new("L", "-L")
        end
    end
end

function format_line(range::AbstractRange, format::LineRangeFormatting)
    if length(range) <= 1
        string(format.prefix, first(range))
    else
        string(format.prefix, first(range), format.separator, last(range))
    end
end

newlines(s::AbstractString) = count(c -> c === '\n', s)
newlines(other) = 0

# Output redirection.
# -------------------
using Logging

"""
Call a function and capture all `stdout` and `stderr` output.

    withoutput(f) --> (result, success, backtrace, output)

where

  * `result` is the value returned from calling function `f`.
  * `success` signals whether `f` has thrown an error, in which case `result` stores the
    `Exception` that was raised.
  * `backtrace` a `Vector{Ptr{Cvoid}}` produced by `catch_backtrace()` if an error is thrown.
  * `output` is the combined output of `stdout` and `stderr` during execution of `f`.

"""
function withoutput(f)
    # Save the default output streams.
    default_stdout = stdout
    default_stderr = stderr

    # Redirect both the `stdout` and `stderr` streams to a single `Pipe` object.
    pipe = Pipe()
    Base.link_pipe!(pipe; reader_supports_async = true, writer_supports_async = true)
    redirect_stdout(pipe.in)
    redirect_stderr(pipe.in)
    # Also redirect logging stream to the same pipe
    logger = ConsoleLogger(pipe.in)

    # Bytes written to the `pipe` are captured in `output` and converted to a `String`.
    output = UInt8[]

    # Run the function `f`, capturing all output that it might have generated.
    # Success signals whether the function `f` did or did not throw an exception.
    result, success, backtrace = with_logger(logger) do
        try
            f(), true, Vector{Ptr{Cvoid}}()
        catch err
            # InterruptException should never happen during normal doc-testing
            # and not being able to abort the doc-build is annoying (#687).
            isa(err, InterruptException) && rethrow(err)

            err, false, catch_backtrace()
        finally
            # Force at least a single write to `pipe`, otherwise `readavailable` blocks.
            println()
            # Restore the original output streams.
            redirect_stdout(default_stdout)
            redirect_stderr(default_stderr)
            # NOTE: `close` must always be called *after* `readavailable`.
            append!(output, readavailable(pipe))
            close(pipe)
        end
    end
    return result, success, backtrace, chomp(String(output))
end

"""
    issubmodule(sub, mod)

Checks whether `sub` is a submodule of `mod`. A module is also considered to be
its own submodule.

E.g. `A.B.C` is a submodule of `A`, `A.B` and `A.B.C`, but it is not a submodule
of `D`, `A.D` nor `A.B.C.D`.
"""
function issubmodule(sub, mod)
    if (sub === Main) && (mod !== Main)
        return false
    end
    (sub === mod) || issubmodule(parentmodule(sub), mod)
end

"""
    isabsurl(url)

Checks whether `url` is an absolute URL (as opposed to a relative one).
"""
isabsurl(url) = occursin(ABSURL_REGEX, url)
const ABSURL_REGEX = r"^[[:alpha:]+-.]+://"

"""
    mdparse(s::AbstractString; mode=:single)

Parses the given string as Markdown using `Markdown.parse`, but strips away the surrounding
layers, such as the outermost `Markdown.MD`. What exactly is returned depends on the `mode`
keyword.

The `mode` keyword argument can be one of the following:

* `:single` (default) -- returns a single block-level object (e.g. `Markdown.Paragraph` or
  `Markdown.Admonition`) and errors if the string parses into multiple blocks.
* `:blocks` -- the function returns a `Vector{Any}` of Markdown blocks.
* `:span` -- Returns a `Vector{Any}` of span-level items, stripping away the outer block.
  This requires the string to parse into a single `Markdown.Paragraph`, the contents of
  which gets returned.
"""
function mdparse(s::AbstractString; mode=:single)
    mode in [:single, :blocks, :span] || throw(ArgumentError("Invalid mode keyword $(mode)"))
    md = Markdown.parse(s)
    if mode == :blocks
        md.content
    elseif length(md.content) == 0
        # case where s == "". We'll just return an empty string / paragraph.
        (mode == :single) ? Markdown.Paragraph(Any[""]) : Any[""]
    elseif (mode == :single || mode == :span) && length(md.content) > 1
        @error "mode == :$(mode) requires the Markdown string to parse into a single block" s md.content
        throw(ArgumentError("Unsuitable string for mode=:$(mode)"))
    else
        @assert length(md.content) == 1
        @assert mode == :span || mode == :single
        if mode == :span && !isa(md.content[1], Markdown.Paragraph)
            @error "mode == :$(mode) requires the Markdown string to parse into a Markdown.Paragraph" s md.content
            throw(ArgumentError("Unsuitable string for mode=:$(mode)"))
        end
        (mode == :single) ? md.content[1] : md.content[1].content
    end
end

# Capturing output in different representations similar to IJulia.jl
function limitstringmime(m::MIME"text/plain", x)
    io = IOBuffer()
    show(IOContext(io, :limit=> true), m, x)
    return String(take!(io))
end
function display_dict(x)
    out = Dict{MIME,Any}()
    x === nothing && return out
    # Always generate text/plain
    out[MIME"text/plain"()] = limitstringmime(MIME"text/plain"(), x)
    for m in [MIME"text/html"(), MIME"image/svg+xml"(), MIME"image/png"(),
              MIME"image/webp"(), MIME"image/gif"(), MIME"image/jpeg"(),
              MIME"text/latex"(), MIME"text/markdown"()]
        showable(m, x) && (out[m] = stringmime(m, x))
    end
    return out
end

"""
    struct Default{T}

Internal wrapper type that is meant to be used in situations where it is necessary to
distinguish whether the user explicitly passed the same value as the default value to a
keyword argument, or whether the keyword argument was not passed at all.

```julia
function foo(; kwarg = Default("default value"))
    if isa(kwarg, Default)
        # User did not explicitly pass a value for kwarg
    else kwarg === "default value"
        # User passed "default value" explicitly
    end
end
```
"""
struct Default{T}
    value :: T
end
Base.getindex(default::Default) = default.value

include("DOM.jl")
include("MDFlatten.jl")
include("TextDiff.jl")
include("Selectors.jl")
include("Markdown2.jl")
include("JSDependencies.jl")

end


text:
mutable struct Num <: Integer
	val::Int
end
Num() = Num(0)

function get_next_word!(b::Array{UInt8, 1}, pos::Num) :: UInt32
	w = Int32(0)
	j = pos.val
	for i in 1:4
		w = w<<8 | UInt32(b[j+1])
		j += 1
		if j >= length(b)
			j = 0
		end
	end
	pos.val = j
	return w
end

function ExpandKey!(key::Array{UInt8, 1}, c::Cipher)
	j = Int(0)
	for i in 1:18
		d = UInt32(0)
		for k in 1:4
			d = d<<8 | UInt32(key[j+1])
			j += 1
			if (j >= length(key))
				j = 0
			end
		end
		c.p[i] ⊻= d
	end

	l, r = UInt32(0), UInt32(0)
	for i in 1:2:18
		l, r = encrypt_block(l, r, c)
		c.p[i], c.p[i+1] = l, r
	end
	for i in 1:2:256
		l, r = encrypt_block(l, r, c)
		c.s0[i], c.s0[i+1] = l, r
	end
	for i in 1:2:256
		l, r = encrypt_block(l, r, c)
		c.s1[i], c.s1[i+1] = l, r
	end
	for i in 1:2:256
		l, r = encrypt_block(l, r, c)
		c.s2[i], c.s2[i+1] = l, r
	end
	for i in 1:2:256
		l, r = encrypt_block(l, r, c)
		c.s3[i], c.s3[i+1] = l, r
	end
end

function expand_key_with_salt!(key::Array{UInt8, 1}, salt::Array{UInt8, 1}, c::Cipher)
	j = Num()
	for i in 1:18
		c.p[i] ⊻= get_next_word!(key, j)
	end

	j = Num()
	l, r = UInt32(0), UInt32(0)
	for i in 1:2:18
		l ⊻= get_next_word!(salt, j)
		r ⊻= get_next_word!(salt, j)
		l, r = encrypt_block(l, r, c)
		c.p[i], c.p[i+1] = l, r
	end
	for i in 1:2:256
		l ⊻= get_next_word!(salt, j)
		r ⊻= get_next_word!(salt, j)
		l, r = encrypt_block(l, r, c)
		c.s0[i], c.s0[i+1] = l, r
	end
	for i in 1:2:256
		l ⊻= get_next_word!(salt, j)
		r ⊻= get_next_word!(salt, j)
		l, r = encrypt_block(l, r, c)
		c.s1[i], c.s1[i+1] = l, r
	end
	for i in 1:2:256
		l ⊻= get_next_word!(salt, j)
		r ⊻= get_next_word!(salt, j)
		l, r = encrypt_block(l, r, c)
		c.s2[i], c.s2[i+1] = l, r
	end
	for i in 1:2:256
		l ⊻= get_next_word!(salt, j)
		r ⊻= get_next_word!(salt, j)
		l, r = encrypt_block(l, r, c)
		c.s3[i], c.s3[i+1] = l, r
	end
end

function encrypt_block(l::UInt32, r::UInt32, c::Cipher) :: Tuple{UInt32, UInt32}
	xl, xr = l, r
	xl ⊻= c.p[1]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[2]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[3]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[4]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[5]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[6]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[7]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[8]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[9]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[10]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[11]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[12]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[13]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[14]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[15]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[16]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[17]
	xr ⊻= c.p[18]
	return xr, xl
end

function decrypt_block(l::UInt32, r::UInt32, c::Cipher) :: Tuple{UInt32, UInt32}
	xl, xr = l, r
	xl ⊻= c.p[18]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[17]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[16]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[15]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[14]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[13]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[12]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[11]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[10]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[9]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[8]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[7]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[6]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[5]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[4]
	xr ⊻= ((c.s0[xl>>24 & 0xff + 1] + c.s1[xl>>16 & 0xff + 1]) ⊻ c.s2[xl>>8 & 0xff + 1]) + c.s3[xl & 0xff + 1] ⊻ c.p[3]
	xl ⊻= ((c.s0[xr>>24 & 0xff + 1] + c.s1[xr>>16 & 0xff + 1]) ⊻ c.s2[xr>>8 & 0xff + 1]) + c.s3[xr & 0xff + 1] ⊻ c.p[2]
	xr ⊻= c.p[1]
	return xr, xl
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 06755ffe-8e78-11eb-0492-154107f31423
using Pkg; Pkg.activate("MLJ_env", shared=true)

# ╔═╡ ac463e7a-8b59-11eb-229e-db560e17c5f5
begin
	using Test
	using Random
	using Distributions
	using PlutoUI
	using Plots
end

# ╔═╡ 8c80e072-8b59-11eb-3c21-a18fe43c4536
md"""
## Deep Learning (DL)

ref. from book **"Data Science from Scratch"**, Chap 19

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ e7373726-8b59-11eb-2a2b-b5138e4f5268
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ f5ee64b2-8b59-11eb-2751-0778efd589cd
html"""

"""

# ╔═╡ 164b4054-8b5a-11eb-03cc-9fc52eed5937
html"""

"""

# ╔═╡ 31d3fa52-8b5d-11eb-017f-4308a2a06930
md"""
Julia implements natively multi-dimensional arrays (namely tensors in other languages/DL frameworks) which we are going to use here.
"""

# ╔═╡ d783876e-8b59-11eb-3ec1-39663403c77d
begin
  const Tensor = Union{BitArray, Array{T}} where {T <: Real}
  const F = Float64
end

# ╔═╡ 25828140-8b5a-11eb-375c-61061d409d46
function mean(t::Tensor; dims::Integer=0, keep_shape=false)
  """
  Compute mean of tensor up to order 3
  FIXME: how to generalize this?
  """
  @assert 0 ≤ dims ≤ length(size(t))
  p = size(t) |> length

  if dims == 0
    return sum(t)/*(size(t)...)

  elseif p == 2
    d, n = dims, size(t)[dims]
    xr, yc = size(t)
    sum.([view(t, ((d == 1 ? i : 1:xr),
            (d == 2 ? i : 1:yc))...) for i ∈ 1:n]) / n

  elseif dims ≤ length(size(t))
    d, n = dims, size(t)[dims]
    xr, yc, zd = size(t)
    sum.([view(t, ((d == 1 ? i : 1:xr),
            (d == 2 ? i : 1:yc), (d == 3 ? i : 1:zd))...) for i ∈ 1:n]) / n
    else
      throws(ArgumentError("Not Implemented"))
    end
end

# ╔═╡ 59828562-8b5a-11eb-32df-c5eb8302912f
begin
  m = [0.8705316232370359 0.809449166245158;
    0.401883439532462 0.3131013871278703;
    0.3040771908917186 0.007551710346251239]

  @test mean(m) ≈ 0.4510990862300826

  ## mean over rows
  @test mean(m; dims=1) ≈ [0.5599935964940647, 0.23832827555344407, 0.10387630041265661]

  ## mean over cols
  @test mean(m; dims=2) ≈ [0.7882461268306082, 0.5650511318596397]
end

# ╔═╡ a9572022-8b5b-11eb-371d-f13966c4d526
begin
  t = zeros(F, (3, 2, 4))
  t[:, :, 1] = [0.164626  0.0597695; 0.604697  0.286739; 0.836432  0.302445]
  t[:, :, 2] = [0.592242 0.505178; 0.923708 0.473227; 0.0661196 0.787647]
  t[:, :, 3] = [0.16899 0.163547; 0.297063 0.385749; 0.326159 0.345408]
  t[:, :, 4] = [0.0881336 0.266533; 0.54628 0.798231; 0.760609 0.464476];

  @test mean(t) ≈ 0.4255836958333334

  @test mean(t, dims=1) ≈ [0.6696730333333333, 1.438564666666667, 1.2964318666666665]
  @test mean(t, dims=2) ≈ [2.6875296, 2.4194747500000005]
  @test mean(t, dims=3) ≈ [0.563677125, 0.8370304, 0.421729, 0.7310656499999999]
end

# ╔═╡ b57b356a-8b5c-11eb-037e-9125388bd979
md"""
Let us implement a (meta-) function which apply a function to each element of a multi-array (aka tensor). In `Julia` this is easy as we just need to use broadcasting
"""

# ╔═╡ 1be547ba-8b5c-11eb-1ae4-cdaa28a20a2b
tensor_apply(fn::Function, t::Tensor) = fn.(t)

# ╔═╡ 20c0eabc-8b5e-11eb-32dd-cd738a6e8f68
double_fn = x -> eltype(x)(2) * x

# ╔═╡ 3c5a9d7c-8b5e-11eb-287f-91f80252e3ee
begin
  _rt = tensor_apply(double_fn, t)

  @test _rt[:, :, 1] ≈ [0.329252 0.119539; 1.209394 0.573478; 1.672864 0.60489]
  @test _rt[:, :, 2] ≈ [1.184484 1.010356; 1.847416 0.946454; 0.1322392 1.575294]
  @test _rt[:, :, 3] ≈ [0.33798 0.327094; 0.594126 0.771498; 0.652318 0.690816]
  @test _rt[:, :, 4] ≈ [0.1762672 0.533066; 1.09256 1.596462; 1.521218 0.928952]
end

# ╔═╡ ddc961a2-8b5e-11eb-06bb-d7335d0d6eaa
md"""
Which is equivalent to:
"""

# ╔═╡ ea2b2c5a-8b5e-11eb-12fb-33885f709e3a
begin
  @test _rt ≈ 2. * t   # ≡ scalar × tensor (implicit broadcast)
  @test _rt ≈ 2. .* t  # explicit broadcast
end

# ╔═╡ 65c8011c-8b5f-11eb-2d14-29cdd068bbf7
begin
  ## Addition of 2 tensors
  @test t + _rt ≈ 3. * t

  ## Difference of 2 tensors
  @test t - _rt ≈ -1. * t
  @test _rt - t ≈ t

  ## ...
end

# ╔═╡ ab6ef588-8b62-11eb-04fa-ab3051ddd7db
html"""

"""

# ╔═╡ ab54224e-8b62-11eb-3f76-8506247531dc
md"""
### The Layer Abstraction

Our fundamental abstraction is the layer, which is something that knows how to apply some function over its input and how to backpropagate gradients.

*Layer Interface*:
  - forward(...)
  - backward(...)
  - params(...)
  - gradients(...)
"""

# ╔═╡ 97901c1a-8b79-11eb-397e-dda6619ddb03
begin
  abstract type AbstractLayer end

  ## FIXME: review
  forward(::AbstractLayer, ::Tensor) = throws(ArgumentError("Not Implemented"))
  backward(::AbstractLayer, ::Tensor) = throws(ArgumentError("Not Implemented"))
  parms(::AbstractLayer) = throws(ArgumentError("Not Implemented"))
  ∇parms(::AbstractLayer) = throws(ArgumentError("Not Implemented"));
end

# ╔═╡ 96ef5f36-8b63-11eb-3503-617be0c0415b
md"""
##### Initialization of a Layer

We will implement three different shemes for initializing our weigths parameters:
  - pick values from random uniform distribution
  - pick values from standerd normal distribution
  - *Xavier* initialisation

"""

# ╔═╡ 86ca6718-8b75-11eb-1d18-c16694e4482d
function init_rand_normal(shape::Tuple;
		seed=42, DT::DataType=Float64,
		kwargs...)
	Random.seed!(seed)
	nd = Normal(haskey(kwargs, :μ) ? kwargs[:μ] : 0.,
		haskey(kwargs, :σ) ? kwargs[:σ] : 1.)
	DT == Float64 ? rand(nd, shape) : DT[rand(nd, shape)...]
end

# ╔═╡ 43d5f21a-8b77-11eb-397e-dda6619ddb03
mₒ = init_rand_normal((3, 2); DT=Float32)

# ╔═╡ b2b89bc8-8b76-11eb-397e-dda6619ddb03
m₁ = init_rand_normal((3, 2); μ=10., σ=2.)

# ╔═╡ a66e2e68-8b63-11eb-0153-3bb8586dd8b2
function init_rand_uniform(shape::Tuple;
    seed=42, DT::DataType=Float64, _kwargs...)
  Random.seed!(seed)
  ud = Uniform(0., 1.)
  DT == Float64 ? rand(ud, shape) : DT[rand(ud, shape)...]
end

# ╔═╡ 3c685ea4-8b76-11eb-0c25-b917a390c8ba
m₂ = init_rand_uniform((3, 2))

# ╔═╡ a654fbbe-8b63-11eb-10a6-5f34a3164a8e
function init_xavier(shape::Tuple;
    seed=42, DT::DataType=Float64, _kwargs...)
  σ = length(shape) / sum(shape)
  init_rand_normal(shape; σ, seed, DT)
end

# ╔═╡ 2c779860-8b77-11eb-397e-dda6619ddb03
m₃ = init_xavier((3, 2); DT=Float32)

# ╔═╡ 6345dd34-8b77-11eb-397e-dda6619ddb03
const INIT_FNs = [init_rand_normal, init_rand_uniform, init_xavier]

# ╔═╡ 3276d728-8b75-11eb-3bcb-95a228498748
md"""
##### Activation functions

We will implement the follwoing activation funtions:
  - Sigmoid
  - Tanh
  - ReLU
"""

# ╔═╡ d2d26a08-8b83-11eb-3975-57c2a8e1fd58
struct Activation <: AbstractLayer
  name::Symbol
  fn::Function
  der_fn::Function
  #
  _type::Symbol
  store::Dict

  function Activation(name, fn, der_fn)
    new(name, fn, der_fn, :activation, Dict())
  end
end

# ╔═╡ 18a07872-8b84-11eb-3046-331119a0dc9b
begin
  function forward(self::Activation, input::Tensor)::Tensor
    self.store[self.name] = self.fn.(input)
    self.store[self.name]
  end

  function backward(self::Activation, ∇p::Tensor)::Tensor
    self.der_fn(self.store[self.name]) .* ∇p
  end

  parms(self::Activation) = []

  ∇parms(self::Activation) = []
end

# ╔═╡ 879291b6-8b7a-11eb-3340-337957fd81b7
md"""
###### Sigmoid function
"""

# ╔═╡ 325b5476-8b75-11eb-0867-739680374b56
Sigmoid = let
	σ = z -> 1. ./ (1. .+ exp.(-z))
    der_σ = z -> σ(z) .* (1 .- σ(z))
	Activation(:sigmoid, σ, der_σ)
end

# ╔═╡ 9337d830-8c17-11eb-397c-c56dc9ffeea4
begin
  @test Sigmoid.name == :sigmoid

  @test Sigmoid.fn(0.) ≈ 0.5
  @test Sigmoid.fn(1.) ≈ 0.73105857863
  @test Sigmoid.fn(5.) ≈ 0.9933071490757
  @test Sigmoid.fn([5. 1. 0.]) ≈ [0.9933071490757 0.73105857863 0.5]

  @test Sigmoid.der_fn(0.) ≈ 0.25
  @test Sigmoid.der_fn(1.) ≈ 0.1966119332414841
  @test Sigmoid.der_fn(5.) ≈ 0.0066480566708051
end

# ╔═╡ 323ad930-8b75-11eb-315c-5398adedfb78
Tanh = let
	tanₕ = z -> (x = exp.(z); y = exp.(-z); (x .- y) ./ (x .+ y))
	der_tanₕ = z -> 1 .- tanₕ.(z) .^ 2
	Activation(:tanh, tanₕ, der_tanₕ)
end

# ╔═╡ 2e1816d2-8c19-11eb-3046-331119a0dc9b
begin
  @test Tanh.name == :tanh

  @test Tanh.fn(0.) ≈ 0.
  @test Tanh.fn(1.) ≈ 0.7615941559557649
  @test Tanh.fn(1.) ≈ -Tanh.fn(-1.)
  @test Tanh.fn(5.) ≈ 0.999909204262595
  @test Tanh.fn(5.) ≈ -Tanh.fn(-5.)

  @test Tanh.der_fn(0.) ≈ 1.0
  @test Tanh.der_fn(1.) ≈ 0.41997434161402614
  @test Tanh.der_fn(5.) ≈ 0.00018158323094408235
end

# ╔═╡ 41080436-8b86-11eb-2d5c-bd298437f953
md"""
###### ReLU function
"""

# ╔═╡ 321c84bc-8b75-11eb-3869-5b3dfd0fcba7
ReLU = let
	relu_fn = z -> max(zero(eltype(z)), z)
	der_relu = z -> (z₀ =	zero(eltype(z)); z .< z₀ ? z₀ : one(eltype(z)))
 	Activation(:relu, relu_fn, der_relu)
end

# ╔═╡ 04a96d56-8c1a-11eb-1de6-fd84daec8930
begin
  @test ReLU.name == :relu

  @test ReLU.fn(0.) ≈ 0.
  @test ReLU.fn(1.) ≈ 1.
  @test ReLU.fn(-1.) ≈ 0.
  @test ReLU.fn.([-1. 1]) ≈ [0. 1.]

  @test ReLU.der_fn(0.1) ≈ 1.0
  @test ReLU.der_fn(1.) ≈ 1.0
  @test ReLU.der_fn(5.) ≈ 1.0
  @test ReLU.der_fn(-5.) ≈ 0.0
end

# ╔═╡ a6305b04-8b63-11eb-25ea-fd23665f9606
md"""
##### Linear (or Dense) Layer
"""

# ╔═╡ 0a8f0468-8b63-11eb-2e53-9fa288538f73
struct Linear <: AbstractLayer
  idim::Integer   ## input dim
  odim::Integer   ## output dim
  ## parameters
  w::Tensor
  b::Tensor
  #
  _type::Symbol
  store::Dict # {Any. Any}

  ## Assume type is DT == Float64
  function Linear(idim::Integer, odim::Integer;
      DT=Float64, init_fn=init_xavier, seed=42, kwargs...)
    ##
    @assert 1 ≤ idim
    @assert 1 ≤ odim
    @assert init_fn ∈ INIT_FNs
    ##
    w = init_fn((odim, idim); seed, DT, kwargs...)
    b = zeros(DT, (odim, 1))
    new(idim, odim, w, b, :dense, Dict())
  end
end

# ╔═╡ 6657fdee-8bb1-11eb-397c-c56dc9ffeea4
# gen_key(self::Linear, prefix::Symbol) = Symbol(prefix, self._id[1])

# ╔═╡ 86d2aac4-8b64-11eb-2ce3-554d89f33f35
function forward(self::Linear, input::Tensor)::Tensor
  self.store[:input] = input  ## Storing the input for backward pass
  self.w * input + self.b
end

# ╔═╡ 964d2ac2-8b72-11eb-28b5-ddc35ed07aa9
function backward(self::Linear, ∇p::Tensor)::Tensor
  ## as each bᵢ is added to output oᵢ, ∇b is the same as output ∇
  self.store[:∇b] = ∇p
  self.store[:∇w] = ∇p * self.store[:input]'
  r = sum.(self.w' * ∇p)
  @assert size(self.w) == size(self.store[:∇w]) "w and ∇w should have same shape: $(size(self.w)) == $(size(self.store[:∇w]))"
  @assert size(self.b) == size(self.store[:∇b]) "b and ∇b should have same shape: $(size(self.b)) == $(size(self.store[:∇b]))"
  r
end

# ╔═╡ a928209e-8b74-11eb-192c-8d53b41fb0a0
parms(self::Linear) = [self.w, self.b]

# ╔═╡ c03c06ce-8b74-11eb-3158-c758e41c696d
∇parms(self::Linear) = [self.store[:∇w], self.store[:∇b]]

# ╔═╡ c01aade4-8b74-11eb-0af8-73ad690afddb
html"""

"""

# ╔═╡ c0024696-8b74-11eb-321f-0dd3d94a8d48
md"""
### Neural Network as a sequence of Layers

We can think of a neural network as a sequence of layers (Linear/Dense, Activation,...).
"""

# ╔═╡ bfe7710e-8b74-11eb-0cc5-4946d943c398
struct Sequential <: AbstractLayer
  layers::Vector{AbstractLayer}
  _type::Symbol

  function Sequential(layers::Vector{AbstractLayer})
    ## output of prev. layer == input of curr layer
    @assert length(layers) > 0
    pl = layers[1]
    for cl ∈ layers[2:end]
      cl._type == :activation && continue
      @assert pl.odim == cl.idim
      pl = cl
    end

    new(layers, :sequential)
  end
end

# ╔═╡ bfc98a2a-8b74-11eb-3458-153202bf0b91
function forward(self::Sequential, input::Tensor)::Tensor
  for (ix, l) ∈ enumerate(self.layers)
    input = forward(l, input)
  end
  input
end

# ╔═╡ bfb091d4-8b74-11eb-0850-ab7a2e10ff52
function backward(self::Sequential, ∇::Tensor)::Tensor
  for (ix, l) ∈ enumerate(reverse(self.layers))
    ∇ = backward(l, ∇)
  end
  ∇
end

# ╔═╡ bf95bfd0-8b74-11eb-115e-0fc0dd1bc84f
parms(self::Sequential) = [p for l ∈ self.layers for p ∈ parms(l)]

# ╔═╡ a9082e56-8b74-11eb-2359-d91cbb8c5e23
∇parms(self::Sequential) = [∇p for l ∈ self.layers for ∇p ∈ ∇parms(l)]

# ╔═╡ 5b0074ac-8b88-11eb-397a-4971fb64f900
md"""
And now we can define a simple network for calculating the XOR function.
"""

# ╔═╡ dc8498c0-8b81-11eb-397e-dda6619ddb03
begin
  xor_seed = 19999
  xor_init = init_rand_normal

  xor_net = Sequential([
      Linear(2, 2; seed=xor_seed, init_fn=xor_init, μ=1., σ=10.),
      Sigmoid,
      Linear(2, 1; seed=xor_seed, init_fn=xor_init, μ=1, σ=10.),
  ])
end

# ╔═╡ 7c12872a-8b88-11eb-1de6-fd84daec8930
html"""

"""

# ╔═╡ 7bf6740e-8b88-11eb-34ee-4db60fb0db8a
md"""
#### Loss and Optimization

"""

# ╔═╡ be7ae9c0-8b8a-11eb-1de6-fd84daec8930
md"""
##### Loss
"""

# ╔═╡ 7bdcd058-8b88-11eb-3340-337957fd81b7
begin
  abstract type AbstractLoss end

  loss(::AbstractLoss, _ŷ::Tensor, _y::Tensor) =
    throws(ArgumentError("Not Implemented"))

  ∇loss(::AbstractLoss, _ŷ::Tensor, _y::Tensor) =
    throws(ArgumentError("Not Implemented"))
end

# ╔═╡ 7bbf92cc-8b88-11eb-3046-331119a0dc9b
md"""
###### Sum squared Error (SSE)
"""

# ╔═╡ 051a9e80-8b8a-11eb-34ee-4db60fb0db8a
begin
  struct SSE <: AbstractLoss
  end

  function loss(_self::SSE, ŷ::Tensor, y::Tensor)::F
    sum((ŷ .- y) .^ 2)
  end

  function ∇loss(_self::AbstractLoss, ŷ::Tensor, y::Tensor)::Tensor
    2 * (ŷ .- y)
  end
end

# ╔═╡ 04dcc506-8b8a-11eb-3046-331119a0dc9b
md"""
##### Optimization
"""

# ╔═╡ 04c3d762-8b8a-11eb-397a-4971fb64f900
begin
  abstract type AbstractOptimizer end

  astep(::AbstractOptimizer, _l::AbstractLayer) =
    throws(ArgumentError("Not Implemented"))
end

# ╔═╡ 247b100e-8b8b-11eb-2d5c-bd298437f953
md"""
###### (Vanilla) Gradient Descent (GD)
"""

# ╔═╡ 33ad413e-8b8b-11eb-170e-0f17904c9f2c
begin
  struct GD <: AbstractOptimizer
    η::F
  end

  function astep(self::GD, sl::AbstractLayer)
    ## sl ≡ seq. layer
    for (_parms, _∇parms) ∈ zip(parms(sl), ∇parms(sl))
      _parms[:] = _parms - self.η .* _∇parms
    end
  end
end

# ╔═╡ e19d82b4-8b8a-11eb-3c22-e788d025e8b4
md"""
###### GD with momentum
"""

# ╔═╡ 0c1df400-8b8c-11eb-1ed2-f78b522d0ae9
begin
	mutable struct MomentumGD <: AbstractOptimizer
		η::F  ## learning rate
		α::F  ## momentum rate
		updates::Vector{Tensor}
		#
		function MomentumGD(;η=0.1, α=0.9, updates=Vector{Tensor}[])
			new(η, α, updates)
		end
	end

	function astep(self::MomentumGD, sl::AbstractLayer)
		## sl ≡ seq. layer
		length(self.updates) == 0 &&
			(self.updates = [zeros(eltype(∇p), size(∇p)) for ∇p ∈ ∇parms(sl)])
		#
		for (_upd, _parms, _∇parms) ∈ zip(self.updates, parms(sl), ∇parms(sl))
			## apply momentum
			_upd[:] = self.α * _upd .+ (1. - self.α) * _∇parms
			## take step
			_parms[:] = _parms - self.η .* _upd
		end
	end
end

# ╔═╡ 21614b5a-8b8c-11eb-01e9-bd2517364992
html"""

"""

# ╔═╡ 212a4800-8b8c-11eb-0b29-659a3be01512
md"""
#### XOR example

We will use xor_net as defined above. But first le us defien a generic training loop.
"""

# ╔═╡ 3e135428-8c39-11eb-0307-d932886e4b10
begin
	const AL = AbstractLayer
	const AA = AbstractArray
	const ALoss = AbstractLoss
	const AOpt = AbstractOptimizer
end

# ╔═╡ 2111ca94-8b8c-11eb-0307-d932886e4b10
begin
  ## Training data
  xs = [0. 0.; 0. 1.; 1. 0.; 1. 1.]
  ys = [0.; 1.; 1.; 0.]
end

# ╔═╡ 580afffa-8c3a-11eb-0b29-659a3be01512
function zipper(x::AA, y::AA)
  @assert size(x)[1] == size(y)[1]

  ([x[ix, :], y[ix, :]] for ix ∈ 1:size(x)[1]) ## returns a gen.
end

# ╔═╡ 254daee6-8db7-11eb-231c-f331a607203c
function xor_accuracy(net::AbstractLayer, xs, ys, ϵ=1e-9):: F
  ncorrect = 0
  for (x, y) ∈ zipper(xs, ys)
    ŷ = forward(net, x)
    @show ŷ, y, abs(ŷ[1] - y[1]) ≤ ϵ
     abs(ŷ[1] - y[1]) ≤ ϵ && (ncorrect += 1)
  end
  ncorrect / length(ys)
end

# ╔═╡ 61904ed8-8b96-11eb-3046-331119a0dc9b
size(xs[1, :]), size(reshape(xs[1, :], 2, :)), size(ys[1, :]), size(reshape(ys[1, :], 1, 1))

# ╔═╡ ac9fbad6-8c3a-11eb-231c-f331a607203c
with_terminal() do
  for (x, y) in zip(xs, ys)
    println(x, " / ", y)
  end
end

# ╔═╡ 9577b10a-8db9-11eb-3046-331119a0dc9b
with_terminal() do
  for (x, y) in zipper(xs, ys)
    println(x, " / ", y)
  end
end

# ╔═╡ 276897b0-8c27-11eb-3046-331119a0dc9b
with_terminal() do
	for p in parms(xor_net); println(p); end
end

# ╔═╡ 8880560a-8db7-11eb-01e9-bd2517364992
xor_accuracy(xor_net, xs, ys)

# ╔═╡ 593338dc-8c3f-11eb-3978-75a4a2df9116
html"""

"""

# ╔═╡ 319822fc-8c2a-11eb-397a-4971fb64f900
md"""
#### FizzBuzz example
"""

# ╔═╡ 9ef5592a-8c37-11eb-3046-331119a0dc9b
md"""
First some helper functions
"""

# ╔═╡ 47eaec34-8c45-11eb-3046-331119a0dc9b
function binary_encoder(x::Integer; n::Integer=10)::BitArray
  vbin::BitArray = fill(false, n)
  for ix ∈ 1:n
    vbin[ix] = x % 2
    x ÷= 2
  end
  vbin
end

# ╔═╡ 4baf8304-8c45-11eb-3340-337957fd81b7
begin
  @test binary_encoder(0) == [0, 0, 0, 0, 0, 0, 0, 0,  0,  0]
  @test binary_encoder(1)  == [1, 0, 0, 0, 0, 0, 0, 0,  0,  0]
  @test binary_encoder(10)  == [0, 1, 0, 1, 0, 0, 0, 0,  0,  0]
  @test binary_encoder(101) == [1, 0, 1, 0, 0, 1, 1, 0,  0,  0]
  @test binary_encoder(999) == [1, 1, 1, 0, 0, 1, 1, 1,  1,  1]
end

# ╔═╡ 5c2d2b14-8c45-11eb-2d5c-bd298437f953
function fizzbuzz_encoder(x::Integer)::BitArray
    if x % 15 == 0
        [0, 0, 0, 1]
    elseif x % 5 == 0
        [0, 0, 1, 0]
    elseif x % 3 == 0
        [0, 1, 0, 0]
    else
        [1, 0, 0, 0]
  end
end

# ╔═╡ 5c168c88-8c45-11eb-3c22-e788d025e8b4
begin
  @test fizzbuzz_encoder(2) == [1, 0, 0, 0]
  @test fizzbuzz_encoder(6) == [0, 1, 0, 0]
  @test fizzbuzz_encoder(10) == [0, 0, 1, 0]
  @test fizzbuzz_encoder(30) == [0, 0, 0, 1]
end

# ╔═╡ 5bfa74da-8c45-11eb-1de6-fd84daec8930
md"""
Second, let us define the accuracy. We will use thsi during the training of our fizzbuzz model.
"""

# ╔═╡ 5bdf8ee8-8c45-11eb-34ee-4db60fb0db8a
function fb_accuracy(net::AbstractLayer, low::Integer, high::Integer):: F
	ncorrect = 0
	for n ∈ low:high
		x = binary_encoder(n)
		r₁ = forward(net, x)
		r₂ = fizzbuzz_encoder(n)
		## NOTE: argmax(r₁) => CartesianIndex
		ŷ, y = argmax(r₁)[1], argmax(r₂)
		ŷ == y && (ncorrect += 1)
	end
	ncorrect / (high - low)
end

# ╔═╡ aa26a570-8c45-11eb-170e-0f17904c9f2c
md"""
Third, let us prepare some taining data (using our helper functions).
"""

# ╔═╡ b324dd40-8c45-11eb-231c-f331a607203c
begin
	## Training data
	rₛ = 101:1023
	xₛ = binary_encoder.(collect(rₛ))
	yₛ = fizzbuzz_encoder.(collect(rₛ))
	(xₛ, yₛ)
end

# ╔═╡ b3092534-8c45-11eb-0b29-659a3be01512
md"""
Fourth, let us define our model and optimizer.
"""

# ╔═╡ b2ecba78-8c45-11eb-0307-d932886e4b10
begin
  const Num_Hidden = 25
  #
  # random seed
  #
  fb_net = Sequential([
     Linear(10, Num_Hidden; init_fn=init_rand_uniform),
     Tanh,
     Linear(Num_Hidden, 4; init_fn=init_rand_uniform),
     Sigmoid
  ])
  #
  fb_opt = MomentumGD(;η=.1, α=.9)
  fb_loss = SSE()
end

# ╔═╡ b2d39ee4-8c45-11eb-1ed2-f78b522d0ae9
md"""
Fifth, let us define our training loop.
"""

# ╔═╡ 3e34451e-8c51-11eb-3978-75a4a2df9116
html"""

"""

# ╔═╡ 48ba66fa-8c51-11eb-34ee-4db60fb0db8a
md"""
#### Softmax and Cross-Entropy Loss

The softmax function converts a vector of real numbers into a vector of probabilities.
"""

# ╔═╡ 489fe8d2-8c51-11eb-3340-337957fd81b7
function softmax(tensor::Tensor)::Tensor
	"""Softmax along ths last dimension"""
	## for numerical stability - subtract largest value
	exps = exp.(tensor .- maximum(tensor))
	sum_exps = sum(exps)
	exps ./ sum_exps
end

# ╔═╡ 4882bcee-8c51-11eb-3046-331119a0dc9b
begin
	@test softmax([1. 3. 2.]) ≈ [0.09003057 0.66524096 0.24472847]
	@test softmax([0. 1. 0.]) ≈ [0.21194155761708 0.5761168847658 0.21194155761708]
	@test softmax([8. 5. 0. 0.]) ≈ [0.9519657197813 0.04739558237461 0.0003193489220309 0.0003193489220309]
end

# ╔═╡ bd860472-8c53-11eb-0307-d932886e4b10
begin
	struct SoftmaxXEntropy <: AbstractLoss; end

	function loss(_self::SoftmaxXEntropy, ŷ::Tensor, y::Tensor; ϵ=1e-20)::F
		probs = softmax(ŷ)
		-sum(log.(probs .+ ϵ) .* y)
	end

	function ∇loss(_self::SoftmaxXEntropy, ŷ::Tensor, y::Tensor)::Tensor
		softmax(ŷ) .- y
	end
end

# ╔═╡ 6474f75a-8c38-11eb-1ed2-f78b522d0ae9
function train_loop(model::AL, xs::AA, ys::AA, loss_fn::ALoss, opt_fn::AOpt;
    epochs=3000, zip_fn=zip, acc_fn=nothing, verbose=false)
	#
	rec_loss = []
	for epoch ∈ 1:epochs
		epoch_loss = 0.
		for (x, y) ∈ zip_fn(xs, ys)
			ŷ = forward(model, x)
			epoch_loss += loss(loss_fn, ŷ, y)
			∇p = ∇loss(loss_fn, ŷ, y)
			backward(model, ∇p)
            astep(opt_fn, model)
		end
		if !isnothing(acc_fn)
			acc = acc_fn(model, 101, 1024)
			epoch % 100 == 0 && (@show epoch, epoch_loss, acc)
		else
			verbose && epoch % 100 == 0 && (@show epoch, epoch_loss)
		end
		push!(rec_loss, epoch_loss)
	end
	rec_loss
end

# ╔═╡ 7d69f7f8-8b8c-11eb-1de6-fd84daec8930
begin
	rec_loss = nothing
	xor_epochs = 1_000 # 3_000
	with_terminal() do
		xor_loss = SSE()
		xor_opt = GD(0.1)
		global rec_loss = train_loop(xor_net, xs, ys, xor_loss, xor_opt;
			epochs=xor_epochs, verbose=true, zip_fn=zipper)
	end
end

# ╔═╡ 70569cb6-8c44-11eb-397a-4971fb64f900
plot(1:xor_epochs, rec_loss, title="XOR training loss", legend=false)

# ╔═╡ ce39053e-8c45-11eb-01e9-bd2517364992
begin
	fb_rec_loss = nothing
	n_epo = 500

	with_terminal() do
		global fb_rec_loss = train_loop(fb_net, xₛ, yₛ, fb_loss, fb_opt;
			verbose=true, epochs=n_epo, acc_fn=fb_accuracy)
	end
end

# ╔═╡ 8090be8e-8c49-11eb-3978-75a4a2df9116
plot(1:n_epo, fb_rec_loss, title="FizzBuzz training loss", legend=false)

# ╔═╡ b0a16012-8c49-11eb-3046-331119a0dc9b
("Results (after training for $(n_epo) epochs): ", fb_accuracy(fb_net, 1, 101))

# ╔═╡ bd69d86a-8c53-11eb-1ed2-f78b522d0ae9
md"""
Now let us redo FizzBuzz, using our `SoftmaxXEntropy`.
"""

# ╔═╡ a3a76e8a-8c56-11eb-3046-331119a0dc9b
begin
	fb_netₓ = Sequential([
		Linear(10, Num_Hidden; init_fn=init_rand_uniform),
		Tanh,
		Linear(Num_Hidden, 4; init_fn=init_rand_uniform),
		# No final Sigmoid
	])
	fb_lossₓ = SoftmaxXEntropy()
	fb_optₓ = MomentumGD(;η=.05, α=.95)
end

# ╔═╡ bd4dc13e-8c53-11eb-170e-0f17904c9f2c
begin
	fb_rec_lossₓ = nothing
	n_epoₓ = 500

	with_terminal() do
		global fb_rec_lossₓ = train_loop(fb_netₓ, xₛ, yₛ, fb_lossₓ, fb_optₓ;
			verbose=true, epochs=n_epoₓ, acc_fn=fb_accuracy)
	end
end

# ╔═╡ e723003c-8c58-11eb-3046-331119a0dc9b
plot(1:n_epoₓ, fb_rec_lossₓ, title="FizzBuzz training loss (Xentropy)", legend=false)

# ╔═╡ 49b61c8e-8c56-11eb-3978-75a4a2df9116
("Results (after training for $(n_epoₓ) epochs): ", fb_accuracy(fb_netₓ, 1, 101))

# ╔═╡ Cell order:
# ╟─8c80e072-8b59-11eb-3c21-a18fe43c4536
# ╠═06755ffe-8e78-11eb-0492-154107f31423
# ╠═ac463e7a-8b59-11eb-229e-db560e17c5f5
# ╟─e7373726-8b59-11eb-2a2b-b5138e4f5268
# ╟─f5ee64b2-8b59-11eb-2751-0778efd589cd
# ╟─164b4054-8b5a-11eb-03cc-9fc52eed5937
# ╟─31d3fa52-8b5d-11eb-017f-4308a2a06930
# ╠═d783876e-8b59-11eb-3ec1-39663403c77d
# ╠═25828140-8b5a-11eb-375c-61061d409d46
# ╠═59828562-8b5a-11eb-32df-c5eb8302912f
# ╠═a9572022-8b5b-11eb-371d-f13966c4d526
# ╟─b57b356a-8b5c-11eb-037e-9125388bd979
# ╠═1be547ba-8b5c-11eb-1ae4-cdaa28a20a2b
# ╠═20c0eabc-8b5e-11eb-32dd-cd738a6e8f68
# ╠═3c5a9d7c-8b5e-11eb-287f-91f80252e3ee
# ╟─ddc961a2-8b5e-11eb-06bb-d7335d0d6eaa
# ╠═ea2b2c5a-8b5e-11eb-12fb-33885f709e3a
# ╠═65c8011c-8b5f-11eb-2d14-29cdd068bbf7
# ╟─ab6ef588-8b62-11eb-04fa-ab3051ddd7db
# ╟─ab54224e-8b62-11eb-3f76-8506247531dc
# ╠═97901c1a-8b79-11eb-397e-dda6619ddb03
# ╟─96ef5f36-8b63-11eb-3503-617be0c0415b
# ╠═86ca6718-8b75-11eb-1d18-c16694e4482d
# ╠═43d5f21a-8b77-11eb-397e-dda6619ddb03
# ╠═b2b89bc8-8b76-11eb-397e-dda6619ddb03
# ╠═a66e2e68-8b63-11eb-0153-3bb8586dd8b2
# ╠═3c685ea4-8b76-11eb-0c25-b917a390c8ba
# ╠═a654fbbe-8b63-11eb-10a6-5f34a3164a8e
# ╠═2c779860-8b77-11eb-397e-dda6619ddb03
# ╠═6345dd34-8b77-11eb-397e-dda6619ddb03
# ╟─3276d728-8b75-11eb-3bcb-95a228498748
# ╠═d2d26a08-8b83-11eb-3975-57c2a8e1fd58
# ╠═18a07872-8b84-11eb-3046-331119a0dc9b
# ╟─879291b6-8b7a-11eb-3340-337957fd81b7
# ╠═325b5476-8b75-11eb-0867-739680374b56
# ╠═9337d830-8c17-11eb-397c-c56dc9ffeea4
# ╠═323ad930-8b75-11eb-315c-5398adedfb78
# ╠═2e1816d2-8c19-11eb-3046-331119a0dc9b
# ╟─41080436-8b86-11eb-2d5c-bd298437f953
# ╠═321c84bc-8b75-11eb-3869-5b3dfd0fcba7
# ╠═04a96d56-8c1a-11eb-1de6-fd84daec8930
# ╟─a6305b04-8b63-11eb-25ea-fd23665f9606
# ╠═0a8f0468-8b63-11eb-2e53-9fa288538f73
# ╠═6657fdee-8bb1-11eb-397c-c56dc9ffeea4
# ╠═86d2aac4-8b64-11eb-2ce3-554d89f33f35
# ╠═964d2ac2-8b72-11eb-28b5-ddc35ed07aa9
# ╠═a928209e-8b74-11eb-192c-8d53b41fb0a0
# ╠═c03c06ce-8b74-11eb-3158-c758e41c696d
# ╟─c01aade4-8b74-11eb-0af8-73ad690afddb
# ╟─c0024696-8b74-11eb-321f-0dd3d94a8d48
# ╠═bfe7710e-8b74-11eb-0cc5-4946d943c398
# ╠═bfc98a2a-8b74-11eb-3458-153202bf0b91
# ╠═bfb091d4-8b74-11eb-0850-ab7a2e10ff52
# ╠═bf95bfd0-8b74-11eb-115e-0fc0dd1bc84f
# ╠═a9082e56-8b74-11eb-2359-d91cbb8c5e23
# ╟─5b0074ac-8b88-11eb-397a-4971fb64f900
# ╠═dc8498c0-8b81-11eb-397e-dda6619ddb03
# ╟─7c12872a-8b88-11eb-1de6-fd84daec8930
# ╟─7bf6740e-8b88-11eb-34ee-4db60fb0db8a
# ╟─be7ae9c0-8b8a-11eb-1de6-fd84daec8930
# ╠═7bdcd058-8b88-11eb-3340-337957fd81b7
# ╟─7bbf92cc-8b88-11eb-3046-331119a0dc9b
# ╠═051a9e80-8b8a-11eb-34ee-4db60fb0db8a
# ╟─04dcc506-8b8a-11eb-3046-331119a0dc9b
# ╠═04c3d762-8b8a-11eb-397a-4971fb64f900
# ╟─247b100e-8b8b-11eb-2d5c-bd298437f953
# ╠═33ad413e-8b8b-11eb-170e-0f17904c9f2c
# ╟─e19d82b4-8b8a-11eb-3c22-e788d025e8b4
# ╠═0c1df400-8b8c-11eb-1ed2-f78b522d0ae9
# ╟─21614b5a-8b8c-11eb-01e9-bd2517364992
# ╟─212a4800-8b8c-11eb-0b29-659a3be01512
# ╠═3e135428-8c39-11eb-0307-d932886e4b10
# ╠═2111ca94-8b8c-11eb-0307-d932886e4b10
# ╠═580afffa-8c3a-11eb-0b29-659a3be01512
# ╠═6474f75a-8c38-11eb-1ed2-f78b522d0ae9
# ╠═254daee6-8db7-11eb-231c-f331a607203c
# ╠═61904ed8-8b96-11eb-3046-331119a0dc9b
# ╠═ac9fbad6-8c3a-11eb-231c-f331a607203c
# ╠═9577b10a-8db9-11eb-3046-331119a0dc9b
# ╠═7d69f7f8-8b8c-11eb-1de6-fd84daec8930
# ╠═276897b0-8c27-11eb-3046-331119a0dc9b
# ╠═8880560a-8db7-11eb-01e9-bd2517364992
# ╠═70569cb6-8c44-11eb-397a-4971fb64f900
# ╟─593338dc-8c3f-11eb-3978-75a4a2df9116
# ╟─319822fc-8c2a-11eb-397a-4971fb64f900
# ╟─9ef5592a-8c37-11eb-3046-331119a0dc9b
# ╠═47eaec34-8c45-11eb-3046-331119a0dc9b
# ╠═4baf8304-8c45-11eb-3340-337957fd81b7
# ╠═5c2d2b14-8c45-11eb-2d5c-bd298437f953
# ╠═5c168c88-8c45-11eb-3c22-e788d025e8b4
# ╟─5bfa74da-8c45-11eb-1de6-fd84daec8930
# ╠═5bdf8ee8-8c45-11eb-34ee-4db60fb0db8a
# ╟─aa26a570-8c45-11eb-170e-0f17904c9f2c
# ╠═b324dd40-8c45-11eb-231c-f331a607203c
# ╟─b3092534-8c45-11eb-0b29-659a3be01512
# ╠═b2ecba78-8c45-11eb-0307-d932886e4b10
# ╟─b2d39ee4-8c45-11eb-1ed2-f78b522d0ae9
# ╠═ce39053e-8c45-11eb-01e9-bd2517364992
# ╠═8090be8e-8c49-11eb-3978-75a4a2df9116
# ╠═b0a16012-8c49-11eb-3046-331119a0dc9b
# ╟─3e34451e-8c51-11eb-3978-75a4a2df9116
# ╟─48ba66fa-8c51-11eb-34ee-4db60fb0db8a
# ╠═489fe8d2-8c51-11eb-3340-337957fd81b7
# ╠═4882bcee-8c51-11eb-3046-331119a0dc9b
# ╠═bd860472-8c53-11eb-0307-d932886e4b10
# ╟─bd69d86a-8c53-11eb-1ed2-f78b522d0ae9
# ╠═a3a76e8a-8c56-11eb-3046-331119a0dc9b
# ╠═bd4dc13e-8c53-11eb-170e-0f17904c9f2c
# ╠═e723003c-8c58-11eb-3046-331119a0dc9b
# ╠═49b61c8e-8c56-11eb-3978-75a4a2df9116

max_stars_repo_path: Data Science From Scratch/19-DL.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/19-DL.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/19-DL.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
using Revise
using ZonalFlow

tspan = (0.0,1000.0);
tsargs  = ( dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            save_noise=false,
            saveat=10
           );

eqs     = [GQL(Λ=0),CE2(),CE2()];

"""
   Single mode excitation
"""
dn      = "validate/stochastic/nonlinear/1.0/"

coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.0,ν₄=1.0,linear=false);
forcing = Stochastic(kf=3,dk=0,ε=1.0); #note dk = 0 => only kf is excited

labels  = ["ql_6x6_m3","ce2_6x6_m3","ce2_6x6_qlic_m3"];

domain  = Domain(extent=(2π,2π),res=(6,6));
prob    = BetaPlane(domain,coeffs,forcing);

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

labels  = ["ql_8x8_m3","ce2_8x8_m3","ce2_8x8_qlic_m3"];

domain  = Domain(extent=(2π,2π),res=(8,8));
prob    = BetaPlane(domain,coeffs,forcing);

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

labels  = ["ql_10x10_m3","ce2_10x10_m3","ce2_10x10_qlic_m3"];

domain  = Domain(extent=(2π,2π),res=(10,10));
prob    = BetaPlane(domain,coeffs,forcing);

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

"""
   Three-mode excitations
"""
dn      = "validate/stochastic/nonlinear/e0.1m5d1/"

coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.001,ν₄=1.0,linear=false);
forcing = Stochastic(kf=5,dk=1,ε=0.1);

# 6x6

labels  = ["ql_6x6","ce2_6x6","ce2_6x6_qlic"];

domain  = Domain(extent=(2π,2π),res=(6,6));
prob    = BetaPlane(domain,coeffs,forcing);

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 8x8

domain  = Domain(extent=(2π,2π),res=(8,8));
prob    = BetaPlane(domain,coeffs,forcing);

labels  = ["ql_8x8","ce2_8x8","ce2_8x8_qlic"];

ql     = integrate(prob,eqs[1],(0.0,10000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 10x10

domain  = Domain(extent=(2π,2π),res=(10,10));
prob    = BetaPlane(domain,coeffs,forcing);

labels  = ["ql_10x10","ce2_10x10","ce2_10x10_qlic"];

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,CE2(poscheckat=10),tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# 12x12

domain  = Domain(extent=(2π,2π),res=(12,12));
prob    = BetaPlane(domain,coeffs,forcing);

labels  = ["ql_12x12","ce2_12x12","ce2_12x12_qlic"];

ql     = integrate(prob,eqs[1],(0.0,5000.0);tsargs...);
write(prob,eqs[1],ql,dn=dn,fn=labels[1])

ce2     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2,dn=dn,fn=labels[2])

ce2_qlic     = integrate(prob,eqs[3],tspan;u0=ql.u[end],tsargs...);
write(prob,eqs[3],ce2_qlic,dn=dn,fn=labels[3])

# GCE2
# eqs     = [NL(),GQL(Λ=prob.d.nx-1),GCE2(Λ=prob.d.nx-1)];
# labels  = ["nl_linear","gql_5_linear","gce2_5_linear"];
# sols    = integrate(prob,eqs,tspan;tsargs...);
# write(prob,eqs,sols,dn="validate/stochastic/6x6/",labels=labels)


text:
using MOBC, vOptGeneric, CPLEX, vOptGeneric, PyPlot, Suppressor
using Base.Test

function random_instance_min(range, n)
	m = vModel(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
	p1 = rand(range, n)
	p2 = rand(range, n)
	w = [rand(range, n) for i = 1:3]
	c = 3*sum.(w).÷5

	@variable(m, x[1:length(p1)], Bin)
	@addobjective(m, Min, dot(x, p1))
	@addobjective(m, Min, dot(x, p2))
	@constraint(m, dot(x, w[1]) >= c[1])
	@constraint(m, dot(x, w[2]) >= c[2])
	@constraint(m, dot(x, w[3]) >= c[3])
	m, p1, p2
end

# srand(14)
function test_min()
	m, p1, p2 = random_instance_min(50:100, 20)
	solve(m, method=:dicho)
	ub = getY_N(m);
	ub = MOBC.toPointList(ub);
	nadir = (maximum(first.(ub)), maximum(last.(ub)));
	@constraint(m, dot(m[:x], p1) <= nadir[1])
	@constraint(m, dot(m[:x], p2) <= nadir[2])
	@constraint(m, m[:x][rand(1:7)] == rand(0:1))
	@constraint(m, m[:x][rand(8:14)] == rand(0:1))
	@constraint(m, m[:x][rand(15:20)] == rand(0:1))
	solve(m, method=:dicho, relax=true)
	lb = getY_N(m);
	lb = MOBC.toSegmentList(lb, nadir, MOBC.Min);
	clf() ; MOBC.plotdualbound(lb) ; plot(first.(ub), last.(ub), "ko")
	MOBC.filterDualBound(lb, ub, true)
end

function random_instance_max(range, n)
	m = vModel(solver = CplexSolver(CPX_PARAM_SCRIND = 0))
	p1 = rand(range, n)
	p2 = rand(range, n)
	w = [rand(range, n) for i = 1:3]
	c = sum.(w).÷2

	@variable(m, x[1:length(p1)], Bin)
	@addobjective(m, Max, dot(x, p1))
	@addobjective(m, Max, dot(x, p2))
	@constraint(m, dot(x, w[1]) <= c[1])
	@constraint(m, dot(x, w[2]) <= c[2])
	@constraint(m, dot(x, w[3]) <= c[3])
	m, p1, p2
end

# srand(14)
function test_max()
	m, p1, p2 = random_instance_max(50:100, 20)
	solve(m, method=:dicho)
	ub = getY_N(m);
	ub = MOBC.toPointList(ub);
	nadir = (minimum(first.(ub)), minimum(last.(ub)));
	@constraint(m, dot(m[:x], p1) >= nadir[1])
	@constraint(m, dot(m[:x], p2) >= nadir[2])
	@constraint(m, m[:x][rand(1:7)] == rand(0:1))
	@constraint(m, m[:x][rand(8:14)] == rand(0:1))
	@constraint(m, m[:x][rand(15:20)] == rand(0:1))
	solve(m, method=:dicho, relax=true)
	lb = getY_N(m);
	lb = MOBC.toSegmentList(lb, nadir, MOBC.Max);
	clf() ; MOBC.plotdualbound(lb) ; plot(first.(ub), last.(ub), "ko")
	MOBC.filterDualBound(lb, ub, true)
end

# ub = MOBC.toPointList([[2,13], [4,10], [6,6], [10,4], [14,1]]);
# lb = MOBC.toSegmentList([[1,15], [3,11], [5,8], [9,6], [14,5]], (14,13));


text:
module Rings
export Ring

"""
    Ring(r::Integer)

An iterator of (m,n) pairs comprising the r'th summation ring.
"""
struct Ring
    r::Int
    Ring(r) = r < 0 ? error("ring must be non-negative") : new(r)
end

function Base.iterate(ring::Ring)
    r = ring.r
    #topbot = ((m,n) for m in (r,-r) for n in -r:r)
    #leftright = ((m,n) for m in 1-r:r-1 for n in (-r,r))
    m = r
    n = -r
    topdone = botdone = leftdone = rightdone = (r == 0)
    item = (m,n)
    state = (m, n, topdone, botdone, leftdone, rightdone)
    return (item, state)
end

function Base.iterate(ring::Ring, state)
    r = ring.r
    (m, n, topdone, botdone, leftdone, rightdone) = state
    rightdone && return nothing
    if leftdone
        # Doing right
        (m,n) == (r-1,-r) && ((m,n) = (-r,r)) # Just finished left
        m == r-2 && (rightdone = true) # last iteration
        m += 1
        return (m,n), (m, n, topdone, botdone, leftdone, rightdone)
    end
    if botdone
        # Doing left
        (m,n) == (-r,r) && (n = -r) # Just finished bottom
        m == r-2 && (leftdone = true)
        m += 1
        return (m,n), (m, n, topdone, botdone, leftdone, rightdone)
    end
    if topdone
        # Doing bottom
        (m,n) == (r,r) && ((m,n) = (-r, -(r+1))) # Just finished top
        n == r-1 && (botdone = true)
        n += 1
        return (m,n), (m, n, topdone, botdone, leftdone, rightdone)
    end
    # Doing top
    n == r-1 && (topdone = true)
    n += 1
    return (m,n), (m, n, topdone, botdone, leftdone, rightdone)
end

Base.IteratorSize(::Type{Ring}) = Base.HasLength()
Base.IteratorEltype(::Type{Ring}) = Base.HasEltype()
Base.eltype(::Type{Ring}) = Tuple{Int,Int}
Base.length(ring::Ring) = ring.r == 0 ? 1 : 8*ring.r

end # module


text:
#####
##### Convenient definitions
#####

"""
    second

A `Float64` constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Δt = 1second`.
"""
const second = 1.0

"""
    minute

A `Float64` constant equal to 60`second`. Useful for increasing the clarity of scripts, e.g. `Δt = 15minute`.
"""
const minute = 60second

"""
    hour

A `Float64` constant equal to 60`minute`. Useful for increasing the clarity of scripts, e.g. `Δt = 3hour`.
"""
const hour   = 60minute

"""
    day

A `Float64` constant equal to 24`hour`. Useful for increasing the clarity of scripts, e.g. `Δt = 0.5day`.
"""
const day    = 24hour

"""
    meter

A `Float64` constant equal to 1.0. Useful for increasing the clarity of scripts, e.g. `Lx = 100meter`.
"""
const meter = 1.0

"""
    kilometer

A `Float64` constant equal to 1000`meter`. Useful for increasing the clarity of scripts, e.g. `Lx = 250kilometer`.
"""
const kilometer = 1000meter

"""
    KiB

A `Float64` constant equal to 1024.0. Useful for increasing the clarity of scripts, e.g. `max_filesize = 250KiB`.
"""
const KiB = 1024.0

"""
    MiB

A `Float64` constant equal to 1024`KiB`. Useful for increasing the clarity of scripts, e.g. `max_filesize = 100MiB`.
"""
const MiB = 1024KiB

"""
    GiB

A `Float64` constant equal to 1024`MiB`. Useful for increasing the clarity of scripts, e.g. `max_filesize = 50GiB`.
"""
const GiB = 1024MiB

"""
    TiB

A `Float64` constant equal to 1024`GiB`. Useful for increasing the clarity of scripts, e.g. `max_filesize = 2TiB`.
"""
const TiB = 1024GiB


text:
using RealNeuralNetworks
using Test
using HDF5
using BigArrays
using BigArrays.GSDicts
using RealNeuralNetworks.NodeNets
using RealNeuralNetworks.SWCs

using LinearAlgebra

const CELL_ID = UInt32(76880)
const EXPANSION= (UInt32(80), UInt32(80), UInt32(40))
const MIP = 4

function get_seg_from_h5()
    # read seg data
    f = h5open("/usr/people/jingpeng/seungmount/research/Ashwin/Scripts/NG_scripts/77605.h5")
    seg = f["main"][1,:,:,101:500]
    close(f)
    seg = reshape(seg, size(seg)[2:4])
    @assert ndims(seg) == 3
    for i in eachindex(seg)
        if seg[i] == 77605
            seg[i] = convert(UInt32,1)
        else
            seg[i] = convert(UInt32, 0)
        end
    end
    return seg
end

function get_seg_from_gs()
    ba = BigArray(GSDict( SEG_PATH ))
    seg = ba[2457:2968, 1777:2288, 16401:16912]
end

function create_fake_seg()
    seg = zeros(UInt32,(100,100,100))
    seg[50,50,:] .= one(UInt32)
    seg[49:52, 49:52, 48:52] .= one(UInt32)
    seg[47:54, 47:54, 71:78] .= one(UInt32)
    return seg
end

@testset "test NodeNets module" begin
    seg = create_fake_seg()
    # @time seg = get_seg_from_h5()
    # @time seg = get_seg_from_gs()
    println("building nodeNet ...")
    #@time nodeNet = NodeNet( seg; obj_id = CELL_ID )
    @time nodeNet = NodeNet( seg; obj_id = one(UInt32) )
    NodeNets.add_offset!(nodeNet, (-one(Float32),-one(Float32),-one(Float32)))
    bin = NodeNets.get_neuroglancer_precomputed(nodeNet)
    # open("/tmp/fake.bin", "w") do f write(f, bin)  end
    @time swc = SWCs.SWC( nodeNet )
    SWCs.stretch_coordinates!(swc, MIP)
    SWCs.save(swc, "/tmp/test.swc")
    @test length(swc) > 1
end


text:
@testset "Steady" begin

    submodel = Submodel(Steady())

    # test provided jacobians
    submodel_jacobian_tests(submodel)

    # test consistency of input/output functions
    submodel_io_tests(submodel)

end


text:
import AxisArrays: AxisArray
using PaddedViews
using OffsetArrays

@testset "Non-Standard Arrays" begin
    function _load_and_save(io, img)
        PNGFiles.save(io, img)
        seek(io, 0)
        PNGFiles.load(io)
    end

    @testset "OffsetArrays" begin
        for CT in (Gray{N0f8}, RGB{N0f8}, Gray{Float32}, RGB{Float32})
            img = OffsetArray(rand(CT, 4, 4), -1, -1)

            ref = _load_and_save(IOBuffer(), parent(img))
            actual = _load_and_save(IOBuffer(), img)
            @test ref == actual
            open(joinpath(PNG_TEST_PATH, "test_offsetarray.png"), "w+") do fio
                ref = _load_and_save(fio, parent(img))
                actual = _load_and_save(fio, img)
                @test ref == actual
            end
        end
    end

    @testset "AxisArrays" begin
        for CT in (Gray{N0f8}, RGB{N0f8}, Gray{Float32}, RGB{Float32})
            img = AxisArray(rand(CT, 4, 4))

            ref = _load_and_save(IOBuffer(), convert(Array, img))
            actual = _load_and_save(IOBuffer(), img)
            @test ref == actual
            open(joinpath(PNG_TEST_PATH, "test_axisarray.png"), "w+") do fio
                ref = _load_and_save(fio, convert(Array, img))
                actual = _load_and_save(fio, img)
                @test ref == actual
            end
        end
    end

    @testset "PaddedViews" begin
        for CT in (Gray{N0f8}, RGB{N0f8}, Gray{Float32}, RGB{Float32})
            img = PaddedView(zero(CT), rand(CT, 4, 4), (-1:5, -1:5))
            ref = _load_and_save(IOBuffer(), collect(img))
            actual = _load_and_save(IOBuffer(), img)
            @test ref == actual
            open(joinpath(PNG_TEST_PATH, "test_paddedviews.png"), "w+") do fio
                ref = _load_and_save(fio, collect(img))
                actual = _load_and_save(fio, img)
                @test ref == actual
            end
        end
    end

    @testset "indexed image with special palatte type" begin
        index = rand(1:5, 3, 4)
        palatte = colorview(RGB{N0f8}, rand(UInt8, 3, 5))
        img = IndirectArray(index, palatte)
        ref = _load_and_save(IOBuffer(), collect(img))
        actual = _load_and_save(IOBuffer(), img)
        @test ref == actual
        @test actual isa IndirectArray{RGB{N0f8}, 2}
    end
end


text:
using Test
using Mmap
using OrderedCollections

include("../../src/io.jl")
include("../../src/strpack.jl")

function packb(values)
    stream_out = IOBuffer()
    packer = Packer(stream_out)
    for value in values
        pack(packer, value)
    end
    return take!(stream_out)
end

function assert_packable(value, packed_value)
    stream_out = IOBuffer()
    packer = Packer(stream_out)
    packer.supports_bytes = true
    pack(packer, value)
    packed = take!(stream_out)

    @test packed == packed_value

    unpacker = Unpacker()
    unpacker.source = MessageFrame(Base.view(packed, :), Tuple{Integer, Integer}[(0, length(packed))])
    unpacked = unpack(unpacker)

    @test unpacked == value
end

@testset "PackStream Tests" begin

    @testset "test_none" begin
        assert_packable(nothing, b"\xC0")
    end

    @testset "test_boolean" begin
        assert_packable(true, b"\xC3")
        assert_packable(false, b"\xC2")
    end

    @testset "test_negative_tiny_int" begin
        for z in -16:-1
            assert_packable(z, UInt8[z + 0x100])
        end
    end

    @testset "test_positive_tiny_int" begin
        for z in 0:127
            assert_packable(z, UInt8[z])
        end
    end

    @testset "test numbers using known byte strings" begin
        @testset "test_positive_int16" begin
            z = 32767
            expected = b"\xc9\x7f\xff"
            assert_packable(z, expected)
        end

        @testset "test_negative_int16" begin
            z = -32768
            expected = b"\xc9\x80\x00"
            assert_packable(z, expected)
        end

        @testset "test_positive_int32" begin
            e = 30
            z = 2 ^ e
            expected = b"\xca@\x00\x00\x00"
            assert_packable(z, expected)
        end

        @testset "test_negative_int32" begin
            e = 30
            z = -(2 ^ e + 1)
            expected = b"\xca\xbf\xff\xff\xff"
            assert_packable(z, expected)
        end

        @testset "test_positive_int64" begin
            e = 62
            z = 2^e
            expected = b"\xcb@\x00\x00\x00\x00\x00\x00\x00"
            assert_packable(z, expected)
        end

        @testset "test_negative_int64" begin
            e = 62
            z = -(2 ^ e + 1)
            expected = b"\xcb\xbf\xff\xff\xff\xff\xff\xff\xff"
            assert_packable(z, expected)
        end

        @testset "test_zero_float64" begin
            zero = 0.0
            expected = b"\xc1\x00\x00\x00\x00\x00\x00\x00\x00"
            assert_packable(zero, expected)
        end

        @testset "test_tau_float64" begin
            tau = 2 * pi
            expected = b"\xc1@\x19!\xfbTD-\x18"
            assert_packable(tau, expected)
        end

        @testset "test_positive_float64" begin
            e = 99.0
            r = 2.0 ^ e + 0.5
            expected = b"\xc1F \x00\x00\x00\x00\x00\x00"
            assert_packable(r, expected)
        end

        @testset "test_negative_float64" begin
            e = 99.0
            r = -(2.0 ^ e + 0.5)
            expected = b"\xc1\xc6 \x00\x00\x00\x00\x00\x00"
            assert_packable(r, expected)
        end
    end

    @testset "test numbers using struct_pack" begin
        @testset "test_positive_int16" begin
            for z in 128:32767
                expected = vcat(b"\xC9", struct_pack(Int16(z)))
                assert_packable(z, expected)
            end
        end

        @testset "test_negative_int16" begin
            for z in -32768:-129
                expected = vcat(b"\xC9", struct_pack(Int16(z)))
                assert_packable(z, expected)
            end
        end

        @testset "test_positive_int32" begin
            for e in 15:30
                z = 2 ^ e
                expected = vcat(b"\xCA", struct_pack(Int32(z)))
                assert_packable(z, expected)
            end
        end

        @testset "test_negative_int32" begin
            for e in 15:30
                z = -(2 ^ e + 1)
                expected = vcat(b"\xCA", struct_pack(Int32(z)))
                assert_packable(z, expected)
            end
        end

        @testset "test_positive_int64" begin
            for e in 31:62
                z = 2^e
                expected = vcat(b"\xCB", struct_pack(Int64(z)))
                assert_packable(z, expected)
            end
        end

        @testset "test_negative_int64" begin
            for e in 31:62
                z = -(2^e + 1)
                expected = vcat(b"\xCB", struct_pack(Int64(z)))
                assert_packable(z, expected)
            end
        end

        @testset "test_zero_float64" begin
            zero = 0.0
            expected = vcat(b"\xC1", struct_pack(Float64(zero)))
            assert_packable(zero, expected)
        end

        @testset "test_tau_float64" begin
            tau = 2 * pi
            expected = vcat(b"\xC1", struct_pack(Float64(tau)))
            assert_packable(tau, expected)
        end

        @testset "test_positive_float64" begin
            for e in 0:99
                r = 2.0 ^ e + 0.5
                expected = vcat(b"\xC1", struct_pack(Float64(r)))
                assert_packable(r, expected)
            end
        end

        @testset "test_negative_float64" begin
            for e in 0:99
                r = -(2.0 ^ e + 0.5)
                expected = vcat(b"\xC1", struct_pack(Float64(r)))
                assert_packable(r, expected)
            end
        end
    end

    @testset "test_empty_bytes" begin
        assert_packable(UInt8[], b"\xCC\x00")
    end

    @testset "test_bytes_8" begin
        assert_packable(Array{UInt8}(b"hello"), b"\xCC\x05hello")
    end

    @testset "test_bytes_16" begin
        b = zeros(UInt8, 40000)
        assert_packable(b, vcat(b"\xCD\x9C\x40", b))
    end

    @testset "test_bytes_32" begin
        b = zeros(UInt8, 80000)
        assert_packable(b, vcat(b"\xCE\x00\x01\x38\x80", b))
    end

    @testset "test_empty_string" begin
        assert_packable("", b"\x80")
    end

    @testset "test_tiny_string" begin
        assert_packable("hello", b"\x85hello")
    end

    @testset "test_string_8" begin
        t = repeat("A", 40)
        assert_packable(t, vcat(b"\xD0\x28", Array{UInt8}(t)))
    end

    @testset "test_string_16" begin
        t = repeat("A", 40000)
        assert_packable(t, vcat(b"\xD1\x9C\x40", Array{UInt8}(t)))
    end

    @testset "test_string_32" begin
        t = repeat("A", 80000)
        assert_packable(t, vcat(b"\xD2\x00\x01\x38\x80", Array{UInt8}(t)))
    end

    @testset "test_unicode_string" begin
        t = "héllö"
        b = Array{UInt8}("héllö")
        assert_packable(t, vcat(UInt8[0x80 + length(b)], b))
    end

    @testset "test_empty_list" begin
        assert_packable([], b"\x90")
    end

    @testset "test_list_8" begin
        l = fill(1, 40)
        assert_packable(l, vcat(b"\xD4\x28", fill(0x01, 40)))
    end

    @testset "test_list_16" begin
        l = fill(1, 40000)
        assert_packable(l, vcat(b"\xD5\x9C\x40", fill(0x01, 40000)))
    end

    @testset "test_list_32" begin
        l = fill(1, 80000)
        assert_packable(l, vcat(b"\xD6\x00\x01\x38\x80", fill(0x01, 80000)))
    end

    @testset "test_nested_lists" begin
        assert_packable([[[]]], b"\x91\x91\x90")
    end

    @testset "test_list_stream" begin
        packed_value = b"\xD7\x01\x02\x03\xDF"
        unpacked_value = [1, 2, 3]
        stream_out = IOBuffer()
        packer = Packer(stream_out)
        pack_list_stream_header(packer)
        pack(packer, 1)
        pack(packer, 2)
        pack(packer, 3)
        pack_end_of_stream(packer)
        packed = take!(stream_out)

        @test packed == packed_value

        unpacker = Unpacker()
        unpacker.source = MessageFrame(Base.view(packed, :), Tuple{Integer, Integer}[(0, length(packed))])
        unpacked = unpack(unpacker)

        @test unpacked == unpacked_value
    end

    @testset "test_empty_map" begin
        assert_packable(Dict(), b"\xA0")
    end

    @testset "test_tiny_map" begin
        d = OrderedDict([("A",1),("B",2)])
        assert_packable(d, b"\xA2\x81A\x01\x81B\x02")
    end

    @testset "test_map_8" begin
        d = OrderedDict([("A$(i)", 1) for i in 0:39])
        b = [packb(("A$(i)", 1)) for i in 0:39]
        assert_packable(d, vcat(b"\xD8\x28", vcat(b...)))
    end

    @testset "test_map_16" begin
        d = OrderedDict([("A$(i)", 1) for i in 0:39999])
        b = [packb(("A$(i)", 1)) for i in 0:39999]
        assert_packable(d, vcat(b"\xD9\x9C\x40", vcat(b...)))
    end

    @testset "test_map_32" begin
        d = OrderedDict([("A$(i)", 1) for i in 0:79999])
        b = [packb(("A$(i)", 1)) for i in 0:79999]
        assert_packable(d, vcat(b"\xDA\x00\x01\x38\x80", vcat(b...)))
    end

    @testset "test_map_stream" begin
        packed_value = b"\xDB\x81A\x01\x81B\x02\xDF"
        unpacked_value = Dict("A"=>1,"B"=>2)
        stream_out = IOBuffer()
        packer = Packer(stream_out)
        pack_map_stream_header(packer)
        pack(packer, "A")
        pack(packer, 1)
        pack(packer, "B")
        pack(packer, 2)
        pack_end_of_stream(packer)
        packed = take!(stream_out)

        @test packed == packed_value

        unpacker = Unpacker()
        unpacker.source = MessageFrame(Base.view(packed, :), Tuple{Integer, Integer}[(0, length(packed))])
        unpacked = unpack(unpacker)

        @test unpacked == unpacked_value
    end

    @testset "test_illegal_signature" begin
        @test_throws MethodError assert_packable(Structure(b"XXX", []), b"\xB0XXX")
    end

    @testset "test_empty_struct" begin
        assert_packable(Structure(UInt8('X')), b"\xB0X")
    end

    @testset "test_tiny_struct" begin
        assert_packable(Structure(UInt8('Z'), ["A", 1]), b"\xB2Z\x81A\x01")
    end

    @testset "test_illegal_value_type" begin
        @test_throws MethodError assert_packable(v"1.0.0", b"\xB0XXX")
    end
end

max_stars_repo_name: virtualgraham/JuliaBolt

max_issues_repo_name: virtualgraham/JuliaBolt

max_forks_repo_name: virtualgraham/JuliaBolt


text:
# This application is licensed under the MIT License.
#
# Copyright (c) 2020 Gustavo A. Hoffmann
#
# https://github.com/gusthoff/chaos-music
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

# using Pkg; Pkg.add("DynamicalSystems")
# using Pkg; Pkg.add("MusicXML")
# using Pkg; Pkg.add("Plots")

using DynamicalSystems
using MusicXML
@importMX # imports all the MusicXML types
# using Plots

##########################################################################
# Create logistic map
##########################################################################

lg = Systems.logistic()
tr = trajectory(lg, 99)

# tr = rand(Float64, 100)

# plot(1:101, tr)

##########################################################################
# MusicXML file: general definition for piano part
#
# Based on example code from MusicXML file
##########################################################################

midiinstrument_piano = MidiInstrument(channel= 1,
                                      program =1,
                                      volume = 100,
                                      pan = 0,
                                      id = "P1-I1")
scorepart_piano = ScorePart(name = "Piano",
                            midiinstrument = midiinstrument_piano,
                            id = "P1")

partlist = PartList(scoreparts = [scorepart_piano])

attributes1_piano = Attributes(
   time = Time(beats = 4, beattype = 4), # 4/4
   divisions = 4, # we want to use 16th notes at minimum
   clef = [Clef(number = 1, sign = "G", line = 2),
           Clef(number = 2, sign = "F", line = 4)], # Piano clefs
   staves = 2, # Piano staves
   key = Key(fifths = 0, mode = "major"), # no accidentals, major key
)

##########################################################################
# Generate notes based on logistic map
##########################################################################

notes = ["C", "D", "E", "F", "G", "A", "B", "C"]

measures_piano = Measure[]
measure1_notes_piano = Note[]

first = true

for (i, j) in enumerate(tr)
    curr_note = notes[round(Int, j*length(notes) + 0.5)]
    push!(measure1_notes_piano,
          Note(pitch = Pitch(step = curr_note, alter = 0, octave = 4),
               duration = 4))

    # Add measure if complete
    if i % 4 == 0
        if first
            push!(measures_piano,
                  Measure(attributes = attributes1_piano,
                          notes = measure1_notes_piano))
            global first = false
        else
            push!(measures_piano,
                  Measure(notes = measure1_notes_piano))
        end
        global measure1_notes_piano = Note[]
    end
end

# Add remaining notes to measure

if length(measure1_notes_piano) > 0
    if first
        push!(measures_piano,
              Measure(attributes = attributes1_piano,
                      notes = measure1_notes_piano))
        global first = false
    else
        push!(measures_piano,
              Measure(notes = measure1_notes_piano))
    end
end

##########################################################################
# MusicXML file: write output file
##########################################################################

part_piano = Part(measures = measures_piano, id = "P1")

score = ScorePartwise(partlist = partlist, parts = [part_piano],)

writemusicxml("myscore.musicxml", score)

max_stars_repo_name: gusthoff/chaos-music

max_issues_repo_name: gusthoff/chaos-music

max_forks_repo_name: gusthoff/chaos-music


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# multi-compliance table for move up

# initialise data relevant to move up
function initMultiCompTable!(sim::Simulation, multiCompTable::MultiCompTable)
	# shorthand names:
	@unpack numAmbs, numStations = sim
	mctd = sim.moveUpData.multiCompTableData

	# set multi-compliance table
	checkMultiCompTable(multiCompTable, sim) # not checking station capacities
	mctd.multiCompTable = multiCompTable
end
function initMultiCompTable!(sim::Simulation, multiCompTableFilename::String)
	multiCompTable = readMultiCompTableFile(multiCompTableFilename)
	initMultiCompTable!(sim, multiCompTable)
end

function multiCompTableMoveUp(sim::Simulation)
	@assert(sim.moveUpData.useMoveUp)

	# shorthand
	@unpack ambulances, stations, numAmbs, numStations = sim
	mctd = sim.moveUpData.multiCompTableData
	multiCompTable = mctd.multiCompTable

	# get movable ambulances (movableAmbs)
	ambMovable = isAmbMovable.(ambulances)
	movableAmbs = ambulances[ambMovable]
	numMovableAmbs = length(movableAmbs)

	if numMovableAmbs == 0
		return moveUpNull()
	end

	# calculate travel time for each movable ambulance
	ambToStationTimes = zeros(Float, numMovableAmbs, numStations) # ambToStationTimes[i,j] = time for movable ambulance i to travel to station j
	for i = 1:numMovableAmbs
		ambToStationTimes[i,:] = ambMoveUpTravelTimes!(sim, movableAmbs[i])
	end

	# for all applicable comp table states, solve assignment problem and calculate total travel time
	ambStationIndicesList = []
	totalTravelTimes = []
	compTableStates = multiCompTable[numMovableAmbs]
	for compTableState in compTableStates
		stationSlots = vcat([[stationIndex for j = 1:stationSlotCount] for (stationIndex, stationSlotCount) in enumerate(compTableState)]...)
		results = Dict()
		ambStationIndices = solveCompTableAssignmentProblem(stationSlots, ambToStationTimes, results = results)
		push!(ambStationIndicesList, ambStationIndices)
		push!(totalTravelTimes, results[:objVal])
	end

	# use comp table state with smallest total travel time
	i = argmin(totalTravelTimes)
	ambStationIndices = ambStationIndicesList[i]
	compTableState = compTableStates[i]

	ambStations = stations[ambStationIndices]

	if checkMode
		# check that comp table state is followed
		compTableStateCopy = deepcopy(compTableState)
		for j in ambStationIndices
			compTableStateCopy[j] -= 1
		end
		@assert(iszero(compTableStateCopy))
	end

	return movableAmbs, ambStations
end

# check that multi-compliance table is valid
function checkMultiCompTable(multiCompTable::MultiCompTable;
	numAmbs::Int = nullIndex, numStations::Int = nullIndex,
	stationCapacities::Union{Vector{Int},Nothing} = nothing)

	@assert(numAmbs == nullIndex || length(multiCompTable) == numAmbs)
	if numStations == nullIndex
		numStations = multiCompTable[1][1]
	end
	for i = 1:numAmbs
		states = multiCompTable[i]
		@assert(length(states) >= 1)
		for state in states
			# state[k] gives number of ambulances to assign to station k
			@assert(length(state) == numStations)
			@assert(all(x -> x >= 0, state))
			@assert(sum(state) == i) # match number of ambs
			if stationCapacities !== nothing
				@assert(all(state .<= stationCapacities))
			end
		end
	end
	return true
end
function checkMultiCompTable(multiCompTable::MultiCompTable, sim::Simulation)
	return checkMultiCompTable(multiCompTable; numAmbs = sim.numAmbs, numStations = sim.numStations, stationCapacities = [s.capacity for s in sim.stations])
end