text:
module Sudoku

using DancingLinks

export create_sudoku, solve_sudoku!, check_solution

const L = 3
const N = L * L

# For each label, there are 9 positions to place. There are 9 labels in each sub-grid. And there
# are 9 sub-grids in total.
const MAX_N_ROWS = N * N * N

# From https:#en.wikipedia.org/wiki/Exact_cover#Sudoku:
#
# - Row-Column: Each intersection of a row and column must contain exactly one number.
# - Row-Number: Each row must contain each number exactly once.
# - Column-Number: Each column must contain each number exactly once.
# - Box-Number: Each box must contain each number exactly once.
const MAX_N_COLS = N * N * 4

# `+ 1` for the root. `* 4` because there are 4 data objects linked together in each row.
const MAX_N = 1 + MAX_N_COLS + MAX_N_ROWS * 4

const Coord = Tuple{Int, Int}

struct SudokuSolver
    d::DLinks
    names::Vector{Tuple{Coord, Int}}  # the metadata, for outputing result
end

function create_sudoku(s::Array{Int})::SudokuSolver
    @assert size(s) == (N, N) "Invalid input"

    # allocates dancing links
    d = let
        L = Vector{Int}(undef, MAX_N)
        R = Vector{Int}(undef, MAX_N)
        U = Vector{Int}(undef, MAX_N)
        D = Vector{Int}(undef, MAX_N)

        h = 0
        C = Vector{Int}(undef, MAX_N)
        S = zeros(Int, MAX_N_COLS)

        DLinks(L, R, U, D, h, C, S)
    end
    names = Vector{Tuple{Coord, Int}}(undef, MAX_N)

    init_headers!(d, MAX_N_COLS)
    i = d.h

    for r in 1:N
        for c in 1:N
            cons = compute_constrait_indices(r, c)

            n_list = if s[r, c] == 0  # empty
                1:N
            else
                n = s[r, c]
                n:n
            end

            for n in n_list
                n_ofs = [n, n, n, 1]
                for ii in 1:4
                    DancingLinks.insert!(d,
                        i + ii, i + mod1(ii - 1, 4), i + mod1(ii + 1, 4),
                        cons[ii] + n_ofs[ii])
                    names[i + ii] = ((r, c), n)
                end
                i += 4
            end
        end
    end

    SudokuSolver(d, names)
end

const N2 = N * N

@inline function compute_constrait_indices(r::Int, c::Int)::Tuple{Int, Int, Int, Int}
    offset = 0
    row_cons = offset + N * (r - 1)

    offset += N2
    col_cons = offset + N * (c - 1)

    offset += N2
    box_cons = offset + N * ((r - 1) ÷ L * L + (c - 1) ÷ L)

    offset += N2
    r_c_cons = offset + ((r - 1) * N + (c - 1))

    (row_cons, col_cons, box_cons, r_c_cons)
end

function solve_sudoku!(sudoku::SudokuSolver)::Union{Array{Int}, Nothing}
    sudoku_solution = Int[]
    solved = dlx!(sudoku.d, sudoku_solution)
    if solved
        sol = zeros(Int, N, N)
        for r in sudoku_solution
            (coord, n) = sudoku.names[r]
            sol[coord...] = n
        end
        sol
    else
        println(stderr, "No solution.")
        nothing
    end
end

const ALL_NUMBERS = Set(1:N)

function check_solution(s::Array{Int}, puzzle::Array{Int})::Bool
    @assert size(s) == (N, N) "Invalid input"

    # Rows
    if !all(map(s -> s == ALL_NUMBERS, mapslices(Set, s, dims=2)))
        return false
    end

    # Columns
    if !all(map(s -> s == ALL_NUMBERS, mapslices(Set, s, dims=1)))
        return false
    end

    # Boxes
    for r in 1:L
        for c in 1:L
            if Set(view(s, (1:L) .+ L * (r - 1), (1:L) .+ L * (c - 1))) != ALL_NUMBERS
                return false
            end
        end
    end

    # Conforms to the completed grids
    all(@. (s == puzzle) | (puzzle == 0))
end

end  # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
import Cairo; const X = Cairo
import Graphics; const G = Graphics
import Colors; const C = Colors

# so, from this we learn that:
# - if the output attribute comes first, we cna handle pdf and png
# - after that, background must come before drawing, but can be mixed
#   with axes.
# - output can set its own axes too, for defaults.
# - heck, it can even set its own background and ink
# - "number of pixels" is defined for png, but not for pdf (duh!)
# in total, then, we need the following stages
# - initial, uninitialized
# - output set (with default coords, default background, default ink)
# - optional paper (background) and axes (these are "absolute" so override)
# - drawing stuff
# where only the first is required.  for the rest, only order matters.
# (ie don't allow baackground or axes once drawwing has started)

function background(ctx)
#    X.save(ctx)
    G.set_coords(ctx, 0, 0, X.width(ctx.surface), X.height(ctx.surface), 0, 1, 0, 1)
    X.rectangle(ctx, 0, 0, 1, 1)
    X.set_source(ctx, parse(C.Colorant, "red"))
    X.fill(ctx)
#    X.restore(ctx)
end

function diagonal(ctx)
    X.set_source(ctx, parse(C.Colorant, "blue"))
    G.set_coords(ctx, 0, 0, X.width(ctx.surface), X.height(ctx.surface), 0, 1, 0, 1)
    X.move_to(ctx, 0, 0)
    X.line_to(ctx, 0.5, 0.5)
    X.stroke(ctx)
end

function rnadom_crap(ctx)
    X.scale(ctx, 0.5, 0.5)
    X.rotate(ctx, pi/2)
end

function middle_bit(ctx)
    # this stuff no longer fucks up the background, because it sets its own
    # coords
    random_crap(ctx)
    # but the order of these is still important if you want to see both
    background(ctx)
    random_crap(ctx)
    diagonal(ctx)
end

ctx = X.CairoContext(X.CairoPDFSurface("formats.pdf", 3*72, 2*72))
middle_bit(ctx)
X.destroy(ctx)

ctx = X.CairoContext(X.CairoRGBSurface(300, 200))
middle_bit(ctx)
X.write_to_png(ctx.surface, "formats.png")
X.destroy(ctx)

ext: jl


text:
# Utilities for tests.

using Distributions
using LinearAlgebra: norm
using Random
using Test: @test

# Add noise on top of one factorisation WH
function add_matrix_noise(;WH::Matrix{T},
                        dist_E::Distribution=Uniform(0,1),
                        rng::AbstractRNG=MersenneTwister(1234)) where T <: Real
    max.(WH + rand(rng, dist_E, size(WH,1), size(WH, 2)), 0)
end

# Simulate one matrix factorisation V = W*H + E
function get_one_simulated_V()
    F, K, N, rng = 10, 5, 25, MersenneTwister(1234)
    W, H = rand(rng, Uniform(0,1), F, K), rand(rng, Uniform(0,1), K, N)

    add_matrix_noise(WH=W*H, dist_E=Uniform(0,0.1), rng=rng)
end

# Simulate different noises on top of one factorisation WH
function get_many_simulated_V()
    F, K, N = 10, 5, 25
    rng = MersenneTwister(1234)
    W, H = rand(rng, Uniform(0,1), F, K), rand(rng, Uniform(0,1), K, N)

    #### simulate different factorisation
    noise_Vs = Matrix{Float64}[]
    noise_names = String[]

    push!(noise_Vs, add_matrix_noise(WH=W*H, dist_E=Uniform(0,0.01), rng=rng))
    push!(noise_names, "uniform_small")
    push!(noise_Vs, add_matrix_noise(WH=W*H, dist_E=Uniform(0,0.2), rng=rng))
    push!(noise_names, "uniform_large")
    push!(noise_Vs, add_matrix_noise(WH=W*H, dist_E=Normal(0,0.01), rng=rng))
    push!(noise_names, "normal_small")
    push!(noise_Vs, add_matrix_noise(WH=W*H, dist_E=Normal(0,0.2), rng=rng))
    push!(noise_names, "normal_large")
    push!(noise_Vs, add_matrix_noise(WH=W*H, dist_E=Gamma(0.05,1), rng=rng))
    push!(noise_names, "gamma_small")
    push!(noise_Vs, add_matrix_noise(WH=W*H, dist_E=Gamma(0.05,4), rng=rng))
    push!(noise_names, "gamma_large")

    noise_Vs, noise_names
end

# Simulate V, W and H.
function simulate_VWH(;F::Integer,
                      N::Integer,
                      K::Integer,
                      dist_W::Distribution=Uniform(0,1),
                      dist_H::Distribution=Uniform(0,1),
                      dist_E::Distribution=Uniform(0,0.1),
                      rng::AbstractRNG=MersenneTwister(1234))
    W = rand(rng, dist_W, F, K)
    H = rand(rng, dist_H, K, N)
    E = rand(rng, dist_E, F, N)

    W*H+E, W, H
end


text:
function dtftau(t,
		fft_u, fft_v,
		A1, A2,
                matr, conjmatr,
		sigma, llambda, Ntaumm)

    # attention ici je n'ai code' que le cas sigma=1

    sigma = 1

    u = ifft(exp(1im * t .* A1) .* fft_u .* matr)
    v = ifft(exp(1im * t .* A1) .* fft_v .* matr)

    du = ifft(exp(1im * t .* A1) * (1im .* A1) .* fft_u) .* matr
    dv = ifft(exp(1im * t .* A1) * (1im .* A1) .* fft_v) .* matr

    z  = (u  + conj(v) ) / 2
    dz = (du + conj(dv)) / 2

    fz1 = 2 * abs(z) .^ 2 * dz + z .^ 2 * conj(dz)
    champu1 = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* fft(conjmatr .* fz1)
    champv1 = -1im * llambda * A2 .* exp.(-1im * t .* A1) .* fft(conjmatr .* conj(fz1))

    fz1 = abs(z) ^ 2 * z
    champu2 = -1im * llambda .* A2 .* exp.(-1im * t * A1) .* (-1im * A1) .* fft(conjmatr * fz1)
    champv2 = -1im * llambda .* A2 .* exp.(-1im * t * A1) .* (-1im * A1) .* fft(conjmatr * conj(fz1))

    champu = champu1 + champu2
    champv = champv1 + champv2

    return champu, champv

end


text:
using DashBootstrapComponents, DashHtmlComponents

loading_spinner = html_div([
    dbc_button("Load", id="loading-button", n_clicks=0),
    dbc_spinner(html_div(id="loading-output")),
]);

callback!(
    app,
    Output("loading-output", "children"),
    Input("loading-button", "n_clicks"),
) do n

    if n > 0
        sleep(1)
        return "Output loaded $n times"
    else
        return "Output not reloaded yet"
    end

end;

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
push!(LOAD_PATH,"../src/")
using Pkg
Pkg.activate(".")

using Documenter, Reporstat

makedocs(
sitename = "Reporstat.jl",
format = Documenter.HTML(),
modules = [Reporstat],
pages=["Home" => "index.md",])

deploydocs(repo = "github.com/mucinoab/Reporstat.git")

max_stars_repo_name: yollotltamayo/Covid

max_issues_repo_name: yollotltamayo/Covid

max_forks_repo_name: yollotltamayo/Covid


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

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

# ╔═╡ 15187690-0403-11eb-2dfd-fd924faa3513
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
    ])
    using Plots, PlutoUI
end

# ╔═╡ 01341648-0403-11eb-2212-db450c299f35
md"_homework 7, version 1_"

# ╔═╡ 06f30b2a-0403-11eb-0f05-8badebe1011d
md"""

# **Homework 7**: _Epidemic modeling_
`18.S191`, Spring 2021

This notebook contains _built-in, live answer checks_! In some exercises you will see a coloured box, which runs a test case on your code, and provides feedback based on the result. Simply edit the code, run it, and the check runs again.

_For MIT students:_ there will also be some additional (secret) test cases that will be run as part of the grading process, and we will look at your notebook and write comments.

Feel free to ask questions!
"""

# ╔═╡ 095cbf46-0403-11eb-0c37-35de9562cebc
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Jazzy Doe", kerberos_id = "jazz")

# you might need to wait until all other cells in this notebook have completed running.
# scroll around the page to see what's up

# ╔═╡ 03a85970-0403-11eb-334a-812b59c0905b
md"""

Submission by: **_$(student.name)_** ($(student.kerberos_id)@mit.edu)
"""

# ╔═╡ 107e65a4-0403-11eb-0c14-37d8d828b469
md"_Let's create a package environment:_"

# ╔═╡ d8797684-0414-11eb-1869-5b1e2c469011
function bernoulli(p::Number)
	rand() < p
end

# ╔═╡ 61789646-0403-11eb-0042-f3b8308f11ba
md"""
## **Exercise 1:** _Agent-based model for an epidemic outbreak -- types_

In this and the following exercises we will develop a simple stochastic model for combined infection and recovery in a population, which may exhibit an **epidemic outbreak** (i.e. a large spike in the number of infectious people).
The population is **well mixed**, i.e. everyone is in contact with everyone else.
[An example of this would be a small school or university in which people are
constantly moving around and interacting with each other.]

The model is an **individual-based** or **agent-based** model:
we explicitly keep track of each individual, or **agent**, in the population and their
infection status. For the moment we will not keep track of their position in space;
we will just assume that there is some mechanism, not included in the model, by which they interact with other individuals.

#### Exercise 1.1

Each agent will have its own **internal state**, modelling its infection status, namely "susceptible", "infectious" or "recovered". We would like to code these as values `S`, `I` and `R`, respectively. One way to do this is using an [**enumerated type**](https://en.wikipedia.org/wiki/Enumerated_type) or **enum**. Variables of this type can take only a pre-defined set of values; the Julia syntax is as follows:
"""

# ╔═╡ 26f84600-041d-11eb-1856-b12a3e5c1dc7
@enum InfectionStatus S I R

# ╔═╡ 271ec5f0-041d-11eb-041b-db46ec1465e0
md"""
We have just defined a new type `InfectionStatus`, as well as names `S`, `I` and `R` that are the (only) possible values that a variable of this type can take.

👉 Define a variable `test_status` whose value is `S`.
"""

# ╔═╡ 7f4e121c-041d-11eb-0dff-cd0cbfdfd606
test_status = missing

# ╔═╡ 7f744644-041d-11eb-08a0-3719cc0adeb7
md"""
👉 Use the `typeof` function to find the type of `test_status`.
"""

# ╔═╡ 88c53208-041d-11eb-3b1e-31b57ba99f05

# ╔═╡ 847d0fc2-041d-11eb-2864-79066e223b45
md"""
👉 Convert `x` to an integer using the `Integer` function. What value does it have? What values do `I` and `R` have?
"""

# ╔═╡ f2792ff5-b0b6-4fcd-94aa-0b6ef048f6ab

# ╔═╡ 860790fc-0403-11eb-2f2e-355f77dcc7af
md"""
#### Exercise 1.2

For each agent we want to keep track of its infection status and the number of *other* agents that it infects during the simulation. A good solution for this is to define a *new type* `Agent` to hold all of the information for one agent, as follows:
"""

# ╔═╡ ae4ac4b4-041f-11eb-14f5-1bcde35d18f2
mutable struct Agent
	status::InfectionStatus
	num_infected::Int64
end

# ╔═╡ ae70625a-041f-11eb-3082-0753419d6d57
md"""
When you define a new type like this, Julia automatically defines one or more **constructors**, which are methods of a generic function with the *same name* as the type. These are used to create objects of that type.

👉 Use the `methods` function to check how many constructors are pre-defined for the `Agent` type.
"""

# ╔═╡ 60a8b708-04c8-11eb-37b1-3daec644ac90

# ╔═╡ 189cae1e-0424-11eb-2666-65bf297d8bdd
md"""
👉 Create an agent `test_agent` with status `S` and `num_infected` equal to 0.
"""

# ╔═╡ 18d308c4-0424-11eb-176d-49feec6889cf
test_agent = missing

# ╔═╡ 190deebc-0424-11eb-19fe-615997093e14
md"""
👉 For convenience, define a new constructor (i.e. a new method for the function) that takes no arguments and creates an `Agent` with status `S` and number infected 0, by calling one of the default constructors that Julia creates. This new method lives *outside* (not inside) the definition of the `struct`. (It is called an **outer constructor**.)

(In Pluto, a struct definition and an outer constructor need to be combined in a single cell using a `begin end` block.)

Let's check that the new method works correctly. How many methods does the constructor have now?

"""

# ╔═╡ 82f2580a-04c8-11eb-1eea-bdb4e50eee3b
Agent()

# ╔═╡ 8631a536-0403-11eb-0379-bb2e56927727
md"""
#### Exercise 1.3
👉 Write functions `set_status!(a)` and `set_num_infected!(a)` which modify the respective fields of an `Agent`. Check that they work. [Note the bang ("`!`") at the end of the function names to signify that these functions *modify* their argument.]

"""

# ╔═╡ 98beb336-0425-11eb-3886-4f8cfd210288
function set_status!(agent::Agent, new_status::InfectionStatus)
	# your code here
end

# ╔═╡ 866299e8-0403-11eb-085d-2b93459cc141
md"""
👉 We will also need functions `is_susceptible` and `is_infected` that check if a given agent is in those respective states.

"""

# ╔═╡ 9a837b52-0425-11eb-231f-a74405ff6e23
function is_susceptible(agent::Agent)

	return missing
end

# ╔═╡ a8dd5cae-0425-11eb-119c-bfcbf832d695
function is_infected(agent::Agent)

	return missing
end

# ╔═╡ 8692bf42-0403-11eb-191f-b7d08895274f
md"""
#### Exericse 1.4
👉 Write a function `generate_agents(N)` that returns a vector of `N` freshly created `Agent`s. They should all be initially susceptible, except one, chosen at random (i.e. uniformly), who is infectious.

"""

# ╔═╡ 7946d83a-04a0-11eb-224b-2b315e87bc84
function generate_agents(N::Integer)

	return missing
end

# ╔═╡ 488771e2-049f-11eb-3b0a-0de260457731
generate_agents(3)

# ╔═╡ 86d98d0a-0403-11eb-215b-c58ad721a90b
md"""
We will also need types representing different infections.

Let's define an (immutable) `struct` called `InfectionRecovery` with parameters `p_infection` and `p_recovery`. We will make it a subtype of an abstract `AbstractInfection` type, because we will define more infection types later.
"""

# ╔═╡ 223933a4-042c-11eb-10d3-852229f25a35
abstract type AbstractInfection end

# ╔═╡ 1a654bdc-0421-11eb-2c38-7d35060e2565
struct InfectionRecovery <: AbstractInfection
	p_infection
	p_recovery
end

# ╔═╡ 2d3bba2a-04a8-11eb-2c40-87794b6aeeac
md"""
#### Exercise 1.5
👉 Write a function `interact!` that takes an affected `agent` of type `Agent`, an `source` of type `Agent` and an `infection` of type `InfectionRecovery`.  It implements a single (one-sided) interaction between two agents:

- If the `agent` is susceptible and the `source` is infectious, then the `source` infects our `agent` with the given infection probability. If the `source` successfully infects the other agent, then its `num_infected` record must be updated.
- If the `agent` is infected then it recovers with the relevant probability.
- Otherwise, nothing happens.

$(html"")
"""

# ╔═╡ 9d0f9564-e393-401f-9dd5-affa4405a9c6
function interact!(agent::Agent, source::Agent, infection::InfectionRecovery)

	# your code here
end

# ╔═╡ b21475c6-04ac-11eb-1366-f3b5e967402d
md"""
Play around with the test case below to test your function! Try changing the definitions of `agent`, `source` and `infection`. Since we are working with randomness, you might want to run the cell multiple times.
"""

# ╔═╡ 9c39974c-04a5-11eb-184d-317eb542452c
let
	agent = Agent(S, 0)
	source = Agent(I, 0)
	infection = InfectionRecovery(0.9, 0.5)

	interact!(agent, source, infection)

	(agent=agent, source=source)
end

# ╔═╡ 619c8a10-0403-11eb-2e89-8b0974fb01d0
md"""
## **Exercise 2:** _Agent-based model for an epidemic outbreak --  Monte Carlo simulation_

In this exercise we will build on Exercise 2 to write a Monte Carlo simulation of how an infection propagates in a population.

Make sure to re-use the functions that we have already written, and introduce new ones if they are helpful! Short functions make it easier to understand what the function does and build up new functionality piece by piece.

You should not use any global variables inside the functions: Each function must accept as arguments all the information it requires to carry out its task. You need to think carefully about what the information each function requires.

#### Exercise 2.1

👉 Write a function `step!` that takes a vector of `Agent`s and an `infection` of type `InfectionRecovery`. It implements a single step of the infection dynamics as follows:

- Choose two random agents: an `agent` and a `source`.
- Apply `interact!(agent, source, infection)`.
- Return `agents`.

"""

# ╔═╡ 2ade2694-0425-11eb-2fb2-390da43d9695
function step!(agents::Vector{Agent}, infection::InfectionRecovery)

end

# ╔═╡ 955321de-0403-11eb-04ce-fb1670dfbb9e
md"""
👉 Write a function `sweep!`. It runs `step!` $N$ times, where $N$ is the number of agents. Thus each agent acts, on average, once per sweep; a sweep is thus the unit of time in our Monte Carlo simulation.
"""

# ╔═╡ 46133a74-04b1-11eb-0b46-0bc74e564680
function sweep!(agents::Vector{Agent}, infection::AbstractInfection)

end

# ╔═╡ 95771ce2-0403-11eb-3056-f1dc3a8b7ec3
md"""
👉 Write a function `simulation` that does the following:

1. Generate the $N$ agents.

2. Run `sweep!` a number $T$ of times. Calculate and store the total number of agents with each status at each step in variables `S_counts`, `I_counts` and `R_counts`.

3. Return the vectors `S_counts`, `I_counts` and `R_counts` in a **named tuple**, with keys `S`, `I` and `R`.

You've seen an example of named tuples before: the `student` variable at the top of the notebook!

_Feel free to store the counts in a different way, as long as the return type is the same._
"""

# ╔═╡ 887d27fc-04bc-11eb-0ab9-eb95ef9607f8
function simulation(N::Integer, T::Integer, infection::AbstractInfection)

	return (S=missing, I=missing, R=missing)
end

# ╔═╡ b92f1cec-04ae-11eb-0072-3535d1118494
simulation(3, 20, InfectionRecovery(0.9, 0.2))

# ╔═╡ 2c62b4ae-04b3-11eb-0080-a1035a7e31a2
simulation(100, 1000, InfectionRecovery(0.005, 0.2))

# ╔═╡ 28db9d98-04ca-11eb-3606-9fb89fa62f36
@bind run_basic_sir Button("Run simulation again!")

# ╔═╡ c5156c72-04af-11eb-1106-b13969b036ca
let
	run_basic_sir

	N = 100
	T = 1000
	sim = simulation(N, T, InfectionRecovery(0.02, 0.002))

	result = plot(1:T, sim.S, ylim=(0, N), label="Susceptible")
	plot!(result, 1:T, sim.I, ylim=(0, N), label="Infectious")
	plot!(result, 1:T, sim.R, ylim=(0, N), label="Recovered")
end

# ╔═╡ bf6fd176-04cc-11eb-008a-2fb6ff70a9cb
md"""
#### Exercise 2.2
Alright! Every time that we run the simulation, we get slightly different results, because it is based on randomness. By running the simulation a number of times, you start to get an idea of the _mean behaviour_ of our model. This is the essence of a Monte Carlo method! You use computer-generated randomness to generate samples.

Instead of pressing the button many times, let's have the computer repeat the simulation. In the next cells, we run your simulation `num_simulations=20` times with $N=100$, $p_\text{infection} = 0.02$, $p_\text{infection} = 0.002$ and $T = 1000$.

Every single simulation returns a named tuple with the status counts, so the result of multiple simulations will be an array of those. Have a look inside the result, `simulations`, and make sure that its structure is clear.
"""

# ╔═╡ 38b1aa5a-04cf-11eb-11a2-930741fc9076
function repeat_simulations(N, T, infection, num_simulations)
	N = 100
	T = 1000

	map(1:num_simulations) do _
		simulation(N, T, infection)
	end
end

# ╔═╡ 80c2cd88-04b1-11eb-326e-0120a39405ea
simulations = repeat_simulations(100, 1000, InfectionRecovery(0.02, 0.002), 20)

# ╔═╡ 80e6f1e0-04b1-11eb-0d4e-475f1d80c2bb
md"""
In the cell below, we plot the evolution of the number of $I$ individuals as a function of time for each of the simulations on the same plot using transparency (`alpha=0.5` inside the plot command).
"""

# ╔═╡ 9cd2bb00-04b1-11eb-1d83-a703907141a7
let
	p = plot()

	for sim in simulations
		plot!(p, 1:1000, sim.I, alpha=.5, label=nothing)
	end

	p
end

# ╔═╡ 95c598d4-0403-11eb-2328-0175ed564915
md"""
👉 Write a function `sir_mean_plot` that returns a plot of the means of $S$, $I$ and $R$ as a function of time on a single graph.
"""

# ╔═╡ 843fd63c-04d0-11eb-0113-c58d346179d6
function sir_mean_plot(simulations::Vector{<:NamedTuple})
	# you might need T for this function, here's a trick to get it:
	T = length(first(simulations).S)

	return missing
end

# ╔═╡ 7f635722-04d0-11eb-3209-4b603c9e843c
sir_mean_plot(simulations)

# ╔═╡ a4c9ccdc-12ca-11eb-072f-e34595520548
let
	T = length(first(simulations).S)

	all_S_counts = map(result -> result.S, simulations)
	all_I_counts = map(result -> result.I, simulations)
	all_R_counts = map(result -> result.R, simulations)

	(S=round.(sum(all_S_counts) ./ length(simulations) ./ 100, digits=4),
	I=round.(sum(all_I_counts) ./ length(simulations) ./ 100, digits=4),
	R=round.(sum(all_R_counts) ./ length(simulations) ./ 100, digits=4))

end

# ╔═╡ dfb99ace-04cf-11eb-0739-7d694c837d59
md"""
👉 Allow $p_\text{infection}$ and $p_\text{recovery}$ to be changed interactively and find parameter values for which you observe an epidemic outbreak.
"""

# ╔═╡ 1c6aa208-04d1-11eb-0b87-cf429e6ff6d0

# ╔═╡ 95eb9f88-0403-11eb-155b-7b2d3a07cff0
md"""
👉 Write a function `sir_mean_error_plot` that does the same as `sir_mean_plot`, which also computes the **standard deviation** $\sigma$ of $S$, $I$, $R$ at each step. Add this to the plot using **error bars**, using the option `yerr=σ` in the plot command; use transparency.

This should confirm that the distribution of $I$ at each step is pretty wide!
"""

# ╔═╡ 287ee7aa-0435-11eb-0ca3-951dbbe69404
function sir_mean_error_plot(simulations::Vector{<:NamedTuple})
	# you might need T for this function, here's a trick to get it:
	T = length(first(simulations).S)

	return missing
end

# ╔═╡ 9611ca24-0403-11eb-3582-b7e3bb243e62
md"""
#### Exercise 2.3

👉 Plot the probability distribution of `num_infected`. Does it have a recognisable shape? (Feel free to increase the number of agents in order to get better statistics.)

"""

# ╔═╡ 26e2978e-0435-11eb-0d61-25f552d2771e

# ╔═╡ 9635c944-0403-11eb-3982-4df509f6a556
md"""
#### Exercse 2.4
👉 What are three *simple* ways in which you could characterise the magnitude (size) of the epidemic outbreak? Find approximate values of these quantities for one of the runs of your simulation.

"""

# ╔═╡ 4ad11052-042c-11eb-3643-8b2b3e1269bc

# ╔═╡ 61c00724-0403-11eb-228d-17c11670e5d1
md"""
## **Exercise 3:** _Reinfection_

In this exercise we will *re-use* our simulation infrastructure to study the dynamics of a different type of infection: there is no immunity, and hence no "recovery" rather, susceptible individuals may now be **re-infected**

#### Exercise 3.1
👉 Make a new infection type `Reinfection`. This has the *same* two fields as `InfectionRecovery` (`p_infection` and `p_recovery`). However, "recovery" now means "becomes susceptible again", instead of "moves to the `R` class.

This new type `Reinfection` should also be a **subtype** of `AbstractInfection`. This allows us to reuse our previous functions, which are defined for the abstract supertype.
"""

# ╔═╡ 8dd97820-04a5-11eb-36c0-8f92d4b859a8

# ╔═╡ 99ef7b2a-0403-11eb-08ef-e1023cd151ae
md"""
👉 Make a *new method* for the `interact!` function that accepts the new infection type as argument, reusing as much functionality as possible from the previous version.

"""

# ╔═╡ bbb103d5-c8f9-485b-9337-40892bb60506

# ╔═╡ 9a13b17c-0403-11eb-024f-9b37e95e211b
md"""
#### Exercise 3.2
👉 Run the simulation 20 times and plot $I$ as a function of time for each one, together with the mean over the 20 simulations (as you did in the previous exercises).

Note that you should be able to re-use the `sweep!` and `simulation` functions , since those should be sufficiently **generic** to work with the new `step!` function! (Modify them if they are not.)

"""

# ╔═╡ 1ac4b33a-0435-11eb-36f8-8f3f81ae7844

# ╔═╡ 9a377b32-0403-11eb-2799-e7e59caa6a45
md"""
👉 Run the new simulation and draw $I$ (averaged over runs) as a function of time. Is the behaviour qualitatively the same or different? Describe what you see.

"""

# ╔═╡ 21c50840-0435-11eb-1307-7138ecde0691

# ╔═╡ 531d13c2-0414-11eb-0acd-4905a684869d
if student.name == "Jazzy Doe"
	md"""
	!!! danger "Before you submit"
	    Remember to fill in your **name** and **Kerberos ID** at the top of this notebook.
	"""
end

# ╔═╡ 4f19e872-0414-11eb-0dfd-e53d2aecc4dc
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 48a16c42-0414-11eb-0e0c-bf52bbb0f618
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ 9cf9080a-04b1-11eb-12a0-17013f2d37f5
md"""
👉 Calculate the **mean number of infectious agents** of our simulations for each time step. Add it to the plot using a heavier line (`lw=3` for "linewidth") by modifying the cell above.

Check the answer yourself: does your curve follow the average trend?

$(hint(md"This exercise requires some creative juggling with arrays, anonymous functions, `map`s, or whatever you see fit!"))
"""

# ╔═╡ 461586dc-0414-11eb-00f3-4984b57bfac5
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ 43e6e856-0414-11eb-19ca-07358aa8b667
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ 41cefa68-0414-11eb-3bad-6530360d6f68
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ 3f5e0af8-0414-11eb-34a7-a71e7aaf6443
yays = [md"Fantastic!", md"Splendid!", md"Great!", md"Yay ❤", md"Great! 🎉", md"Well done!", md"Keep it up!", md"Good job!", md"Awesome!", md"You got the right answer!", md"Let's move on to the next section."]

# ╔═╡ 3d88c056-0414-11eb-0025-05d3aff1588b
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ 3c0528a0-0414-11eb-2f68-a5657ab9e73d
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 7c515a7a-04d5-11eb-0f36-4fcebff709d5
if !@isdefined(set_status!)
	not_defined(:set_status!)
else
	let
		agent = Agent(I,2)

		set_status!(agent, R)

		if agent.status == R
			correct()
		else
			keep_working()
		end
	end
end

# ╔═╡ c4a8694a-04d4-11eb-1eef-c9e037e6b21f
if !@isdefined(is_susceptible)
	not_defined(:is_susceptible)
else
	let
		result1 = is_susceptible(Agent(I,2))
		result2 = is_infected(Agent(I,2))

		if result1 isa Missing || result2 isa Missing
			still_missing()
		elseif !(result1 isa Bool) || !(result2 isa Bool)
			keep_working(md"Make sure that you return either `true` or `false`.")
		elseif result1 === false && result2 === true
			if is_susceptible(Agent(S,3)) && !is_infected(Agent(R,9))
				correct()
			else
				keep_working()
			end
		else
			keep_working()
		end
	end
end

# ╔═╡ 393041ec-049f-11eb-3089-2faf378445f3
if !@isdefined(generate_agents)
	not_defined(:generate_agents)
else
	let
		result = generate_agents(4)

		if result isa Missing
			still_missing()
		elseif result isa Nothing
			keep_working("The function returned `nothing`. Did you forget to return something?")
		elseif !(result isa Vector) || !all(x -> x isa Agent, result)
			keep_working(md"Make sure that you return an array of objects of the type `Agent`.")
		elseif length(result) != 4
			almost(md"Make sure that you return `N` agents.")
		elseif length(Set(result)) != 4
			almost(md"You returned the **same** agent `N` times. You need to call the `Agent` constructor `N` times, not once.")
		else
			if sum(a -> a.status == I, result) != 1
				almost(md"Exactly one of the agents should be infectious.")
			else
				correct()
			end
		end
	end
end

# ╔═╡ 759bc42e-04ab-11eb-0ab1-b12e008c02a9
if !@isdefined(interact!)
	not_defined(:interact!)
else
	let
		agent = Agent(S, 9)
		source = Agent(I, 0)
		interact!(agent, source, InfectionRecovery(0.0, 1.0))

		if source.status != I || source.num_infected != 0
			keep_working(md"The `source` should not be modified if no infection occured.")
		elseif agent.status != S
			keep_working(md"The `agent` should get infected with the right probability.")
		else
			agent = Agent(S, 9)
			source = Agent(S, 0)
			interact!(agent, source, InfectionRecovery(1.0, 1.0))

			if source.status != S || source.num_infected != 0 || agent.status != S
				keep_working(md"The `agent` should get infected with the right probability if the source is infectious.")
			else
				agent = Agent(S, 9)
				source = Agent(I, 3)
				interact!(agent, source, InfectionRecovery(1.0, 1.0))

				if agent.status == R
					almost(md"The agent should not recover immediately after becoming infectious.")
				elseif agent.status == S
					keep_working(md"The `agent` should recover from an infectious state with the right probability.")
				elseif source.status != I || source.num_infected != 4
					almost(md"The `source` did not get updated correctly after infecting the `agent`.")
				else
					correct(md"Your function treats the **susceptible** agent case correctly!")
				end
			end
		end
	end
end

# ╔═╡ 1491a078-04aa-11eb-0106-19a3cf1e94b0
if !@isdefined(interact!)
	not_defined(:interact!)
else
	let
		agent = Agent(I, 9)
		source = Agent(S, 0)

		interact!(agent, source, InfectionRecovery(1.0, 1.0))

		if source.status != S || source.num_infected != 0
			keep_working(md"The `source` should not be modified if `agent` is infectious.")
		elseif agent.status != R
			keep_working(md"The `agent` should recover from an infectious state with the right probability.")
		elseif agent.num_infected != 9
			keep_working(md"`agent.num_infected` should not be modified if `agent` is infectious.")
		else
			let
				agent = Agent(I, 9)
				source = Agent(R, 0)

				interact!(agent, source, InfectionRecovery(1.0, 0.0))

				if agent.status == R
					keep_working(md"The `agent` should recover from an infectious state with the right probability.")
				else
					correct(md"Your function treats the **infectious** agent case correctly!")
				end
			end
		end
	end
end

# ╔═╡ f8e05d94-04ac-11eb-26d4-6f1d2c5ed272
if !@isdefined(interact!)
	not_defined(:interact!)
else
	let
		agent = Agent(R, 9)
		source = Agent(I, 0)
		interact!(agent, source, InfectionRecovery(1.0, 1.0))

		if source.status != I || source.num_infected != 0
			keep_working(md"The `source` should not be modified if no infection occured.")
		elseif agent.status != R || agent.num_infected != 9
			keep_working(md"The `agent` should not be momdified if it is in a recoved state.")
		else
			correct(md"Your function treats the **recovered** agent case correctly!")
		end
	end
end

# ╔═╡ 39dffa3c-0414-11eb-0197-e72b299e9c63
bigbreak = html"

";

# ╔═╡ 2b26dc42-0403-11eb-205f-cd2c23d8cb03
bigbreak

# ╔═╡ 5689841e-0414-11eb-0492-63c77ddbd136
bigbreak

# ╔═╡ Cell order:
# ╟─01341648-0403-11eb-2212-db450c299f35
# ╟─03a85970-0403-11eb-334a-812b59c0905b
# ╟─06f30b2a-0403-11eb-0f05-8badebe1011d
# ╠═095cbf46-0403-11eb-0c37-35de9562cebc
# ╟─107e65a4-0403-11eb-0c14-37d8d828b469
# ╠═15187690-0403-11eb-2dfd-fd924faa3513
# ╟─2b26dc42-0403-11eb-205f-cd2c23d8cb03
# ╠═d8797684-0414-11eb-1869-5b1e2c469011
# ╟─61789646-0403-11eb-0042-f3b8308f11ba
# ╠═26f84600-041d-11eb-1856-b12a3e5c1dc7
# ╟─271ec5f0-041d-11eb-041b-db46ec1465e0
# ╠═7f4e121c-041d-11eb-0dff-cd0cbfdfd606
# ╟─7f744644-041d-11eb-08a0-3719cc0adeb7
# ╠═88c53208-041d-11eb-3b1e-31b57ba99f05
# ╟─847d0fc2-041d-11eb-2864-79066e223b45
# ╠═f2792ff5-b0b6-4fcd-94aa-0b6ef048f6ab
# ╟─860790fc-0403-11eb-2f2e-355f77dcc7af
# ╠═ae4ac4b4-041f-11eb-14f5-1bcde35d18f2
# ╟─ae70625a-041f-11eb-3082-0753419d6d57
# ╠═60a8b708-04c8-11eb-37b1-3daec644ac90
# ╟─189cae1e-0424-11eb-2666-65bf297d8bdd
# ╠═18d308c4-0424-11eb-176d-49feec6889cf
# ╟─190deebc-0424-11eb-19fe-615997093e14
# ╠═82f2580a-04c8-11eb-1eea-bdb4e50eee3b
# ╟─8631a536-0403-11eb-0379-bb2e56927727
# ╠═98beb336-0425-11eb-3886-4f8cfd210288
# ╟─7c515a7a-04d5-11eb-0f36-4fcebff709d5
# ╟─866299e8-0403-11eb-085d-2b93459cc141
# ╠═9a837b52-0425-11eb-231f-a74405ff6e23
# ╠═a8dd5cae-0425-11eb-119c-bfcbf832d695
# ╟─c4a8694a-04d4-11eb-1eef-c9e037e6b21f
# ╟─8692bf42-0403-11eb-191f-b7d08895274f
# ╠═7946d83a-04a0-11eb-224b-2b315e87bc84
# ╠═488771e2-049f-11eb-3b0a-0de260457731
# ╟─393041ec-049f-11eb-3089-2faf378445f3
# ╟─86d98d0a-0403-11eb-215b-c58ad721a90b
# ╠═223933a4-042c-11eb-10d3-852229f25a35
# ╠═1a654bdc-0421-11eb-2c38-7d35060e2565
# ╟─2d3bba2a-04a8-11eb-2c40-87794b6aeeac
# ╠═9d0f9564-e393-401f-9dd5-affa4405a9c6
# ╟─b21475c6-04ac-11eb-1366-f3b5e967402d
# ╠═9c39974c-04a5-11eb-184d-317eb542452c
# ╟─759bc42e-04ab-11eb-0ab1-b12e008c02a9
# ╟─1491a078-04aa-11eb-0106-19a3cf1e94b0
# ╟─f8e05d94-04ac-11eb-26d4-6f1d2c5ed272
# ╟─619c8a10-0403-11eb-2e89-8b0974fb01d0
# ╠═2ade2694-0425-11eb-2fb2-390da43d9695
# ╟─955321de-0403-11eb-04ce-fb1670dfbb9e
# ╠═46133a74-04b1-11eb-0b46-0bc74e564680
# ╟─95771ce2-0403-11eb-3056-f1dc3a8b7ec3
# ╠═887d27fc-04bc-11eb-0ab9-eb95ef9607f8
# ╠═b92f1cec-04ae-11eb-0072-3535d1118494
# ╠═2c62b4ae-04b3-11eb-0080-a1035a7e31a2
# ╠═c5156c72-04af-11eb-1106-b13969b036ca
# ╟─28db9d98-04ca-11eb-3606-9fb89fa62f36
# ╟─bf6fd176-04cc-11eb-008a-2fb6ff70a9cb
# ╠═38b1aa5a-04cf-11eb-11a2-930741fc9076
# ╠═80c2cd88-04b1-11eb-326e-0120a39405ea
# ╟─80e6f1e0-04b1-11eb-0d4e-475f1d80c2bb
# ╠═9cd2bb00-04b1-11eb-1d83-a703907141a7
# ╟─9cf9080a-04b1-11eb-12a0-17013f2d37f5
# ╟─95c598d4-0403-11eb-2328-0175ed564915
# ╠═843fd63c-04d0-11eb-0113-c58d346179d6
# ╠═7f635722-04d0-11eb-3209-4b603c9e843c
# ╠═a4c9ccdc-12ca-11eb-072f-e34595520548
# ╟─dfb99ace-04cf-11eb-0739-7d694c837d59
# ╠═1c6aa208-04d1-11eb-0b87-cf429e6ff6d0
# ╟─95eb9f88-0403-11eb-155b-7b2d3a07cff0
# ╠═287ee7aa-0435-11eb-0ca3-951dbbe69404
# ╟─9611ca24-0403-11eb-3582-b7e3bb243e62
# ╠═26e2978e-0435-11eb-0d61-25f552d2771e
# ╟─9635c944-0403-11eb-3982-4df509f6a556
# ╠═4ad11052-042c-11eb-3643-8b2b3e1269bc
# ╟─61c00724-0403-11eb-228d-17c11670e5d1
# ╠═8dd97820-04a5-11eb-36c0-8f92d4b859a8
# ╟─99ef7b2a-0403-11eb-08ef-e1023cd151ae
# ╠═bbb103d5-c8f9-485b-9337-40892bb60506
# ╟─9a13b17c-0403-11eb-024f-9b37e95e211b
# ╠═1ac4b33a-0435-11eb-36f8-8f3f81ae7844
# ╟─9a377b32-0403-11eb-2799-e7e59caa6a45
# ╠═21c50840-0435-11eb-1307-7138ecde0691
# ╟─5689841e-0414-11eb-0492-63c77ddbd136
# ╟─531d13c2-0414-11eb-0acd-4905a684869d
# ╟─4f19e872-0414-11eb-0dfd-e53d2aecc4dc
# ╟─48a16c42-0414-11eb-0e0c-bf52bbb0f618
# ╟─461586dc-0414-11eb-00f3-4984b57bfac5
# ╟─43e6e856-0414-11eb-19ca-07358aa8b667
# ╟─41cefa68-0414-11eb-3bad-6530360d6f68
# ╟─3f5e0af8-0414-11eb-34a7-a71e7aaf6443
# ╟─3d88c056-0414-11eb-0025-05d3aff1588b
# ╟─3c0528a0-0414-11eb-2f68-a5657ab9e73d
# ╟─39dffa3c-0414-11eb-0197-e72b299e9c63


text:
using FactCheck, Images, Color, FixedPointNumbers
import Images.ColorTypes, Images.Gray, Images.GrayAlpha

macro test_colortype_approx_eq(a, b)
    :(test_colortype_approx_eq($(esc(a)), $(esc(b)), $(string(a)), $(string(b))))
end

facts("Colortypes") do
    function test_colortype_approx_eq(a::ColorValue, b::ColorValue, astr, bstr)
        @fact typeof(a) --> typeof(b)
        n = length(fieldnames(typeof(a)))
        for i = 1:n
            @fact getfield(a, i) --> roughly(getfield(b,i))
        end
    end

    function test_colortype_approx_eq(a::AbstractAlphaColorValue, b::AbstractAlphaColorValue, astr, bstr)
        @fact typeof(a) --> typeof(b)
        n = length(fieldnames(typeof(a.c)))
        for i = 1:n
            @fact getfield(a.c, i) --> roughly(getfield(b.c, i))
        end
        @fact a.alpha --> roughly(b.alpha)
    end

    context("General") do
        @fact length(RGB) --> 3
        @fact length(ColorTypes.BGR) --> 3
        @fact length(AlphaColorValue{ColorTypes.BGR{Float32},Float32}) --> 4  # TypeConstructor mess, would be length(BGRA) otherwise
        @fact length(ColorTypes.AlphaColor{RGB, Ufixed8}) --> 4  # length(ARGB)
        @fact length(ColorTypes.RGB1) --> 4
        @fact length(ColorTypes.RGB4) --> 4
        @fact length(RGB24) --> 4
        @fact length(ColorTypes.Gray) --> 1
        @fact length(AlphaColorValue{ColorTypes.Gray,Ufixed8}) --> 2
        @fact length(ColorTypes.Gray24) --> 4
        # @fact length(ColorTypes.AGray32) --> 4
        @fact length(ColorTypes.ARGB{Float32}) --> 4
        @fact length(ColorTypes.YIQ) --> 3
        c = ColorTypes.AlphaColor(RGB{Ufixed8}(0.8,0.2,0.4), Ufixed8(0.5))
        @fact length(c) --> 4
        @fact length(zero(GrayAlpha{Float32})) --> 2
        @fact length(one(RGB{Float32})) --> 3
    end

    context("Arithmetic with Gray") do
        cf = Gray{Float32}(0.1)
        ccmp = Gray{Float32}(0.2)
        @fact 2*cf --> ccmp
        @fact cf*2 --> ccmp
        @fact ccmp/2 --> cf
        @fact 2.0f0*cf --> ccmp
        @test_colortype_approx_eq cf*cf Gray{Float32}(0.01)
        @test_colortype_approx_eq cf^2 Gray{Float32}(0.01)
        @test_colortype_approx_eq cf^3.0f0 Gray{Float32}(0.001)
        @fact eltype(2.0*cf) --> Float64
        cu = Gray{Ufixed8}(0.1)
        @fact 2*cu --> Gray(2*cu.val)
        @fact 2.0f0*cu --> Gray(2.0f0*cu.val)
        f = Ufixed8(0.5)
        @fact (f*cu).val --> roughly(f*cu.val)
        @fact 2.*cf --> ccmp
        @fact cf.*2 --> ccmp
        @fact cf/2.0f0 --> Gray{Float32}(0.05)
        @fact cu/2 --> Gray(cu.val/2)
        @fact cu/0.5f0 --> Gray(cu.val/0.5f0)
        @fact cf+cf --> ccmp
        @fact cf --> isfinite
        @fact cf --> not(isinf)
        @fact cf --> not(isnan)
        @fact Gray(NaN) --> not(isfinite)
        @fact Gray(NaN) --> not(isinf)
        @fact Gray(NaN) --> isnan
        @fact Gray(Inf) --> not(isfinite)
        @fact Gray(Inf) --> isinf
        @fact Gray(Inf) --> not(isnan)
        @fact abs(Gray(0.1)) --> roughly(0.1)
        @fact eps(Gray{Ufixed8}) --> Gray(eps(Ufixed8))  # #282

        acu = Gray{Ufixed8}[cu]
        acf = Gray{Float32}[cf]
        @fact typeof(acu+acf) --> Vector{Gray{Float32}}
        @fact typeof(2*acf) --> Vector{Gray{Float32}}
        @fact typeof(0x02*acu) --> Vector{Gray{Float32}}
        @fact typeof(acu/2) --> Vector{Gray{typeof(Ufixed8(0.5)/2)}}
        @fact typeof(acf.^2) --> Vector{Gray{Float32}}
    end

    context("Arithmetic with GrayAlpha") do
        p1 = GrayAlpha{Float32}(Gray(0.8), 0.2)
        p2 = GrayAlpha{Float32}(Gray(0.6), 0.3)
        @test_colortype_approx_eq p1+p2 GrayAlpha{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq (p1+p2)/2 GrayAlpha{Float32}(Gray(0.7),0.25)
        @test_colortype_approx_eq 0.4f0*p1+0.6f0*p2 GrayAlpha{Float32}(Gray(0.68),0.26)
        @test_colortype_approx_eq ([p1]+[p2])[1] GrayAlpha{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1].+[p2])[1] GrayAlpha{Float32}(Gray(1.4),0.5)
        @test_colortype_approx_eq ([p1]/2)[1] GrayAlpha{Float32}(Gray(0.4),0.1)
        @test_colortype_approx_eq (0.4f0*[p1]+0.6f0*[p2])[1] GrayAlpha{Float32}(Gray(0.68),0.26)
    end

    context("Arithemtic with RGB") do
        cf = RGB{Float32}(0.1,0.2,0.3)
        ccmp = RGB{Float32}(0.2,0.4,0.6)
        @fact 2*cf --> ccmp
        @fact cf*2 --> ccmp
        @fact ccmp/2 --> cf
        @fact 2.0f0*cf --> ccmp
        @fact eltype(2.0*cf) --> Float64
        cu = RGB{Ufixed8}(0.1,0.2,0.3)
        @fact 2*cu --> RGB(2*cu.r, 2*cu.g, 2*cu.b)
        @fact 2.0f0*cu --> RGB(2.0f0*cu.r, 2.0f0*cu.g, 2.0f0*cu.b)
        f = Ufixed8(0.5)
        @fact (f*cu).r --> roughly(f*cu.r)
        @fact 2.*cf --> ccmp
        @fact cf.*2 --> ccmp
        @fact cf/2.0f0 --> RGB{Float32}(0.05,0.1,0.15)
        @fact cu/2 --> RGB(cu.r/2,cu.g/2,cu.b/2)
        @fact cu/0.5f0 --> RGB(cu.r/0.5f0, cu.g/0.5f0, cu.b/0.5f0)
        @fact cf+cf --> ccmp
        @test_colortype_approx_eq (cf*[0.8f0])[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq ([0.8f0]*cf)[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq (cf.*[0.8f0])[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @test_colortype_approx_eq ([0.8f0].*cf)[1] RGB{Float32}(0.8*0.1,0.8*0.2,0.8*0.3)
        @fact cf --> isfinite
        @fact cf --> not(isinf)
        @fact cf --> not(isnan)
        @fact RGB(NaN, 1, 0.5) --> not(isfinite)
        @fact RGB(NaN, 1, 0.5) --> not(isinf)
        @fact RGB(NaN, 1, 0.5) --> isnan
        @fact RGB(1, Inf, 0.5) --> not(isfinite)
        @fact RGB(1, Inf, 0.5) --> isinf
        @fact RGB(1, Inf, 0.5) --> not(isnan)
        @fact abs(RGB(0.1,0.2,0.3)) --> roughly(0.6)

        acu = RGB{Ufixed8}[cu]
        acf = RGB{Float32}[cf]
        @fact typeof(acu+acf) --> Vector{RGB{Float32}}
        @fact typeof(2*acf) --> Vector{RGB{Float32}}
        @fact typeof(convert(UInt8, 2)*acu) --> Vector{RGB{Float32}}
        @fact typeof(acu/2) --> Vector{RGB{typeof(Ufixed8(0.5)/2)}}

        c = ColorTypes.Gray{Ufixed16}(0.8)
        @fact convert(RGB, c) --> RGB{Ufixed16}(0.8,0.8,0.8)
        @fact convert(RGB{Float32}, c) --> RGB{Float32}(0.8,0.8,0.8)
        r4 = ColorTypes.RGB4(1,0,0)
        @fact convert(RGB, r4) --> RGB(1,0,0)
        @fact convert(RGB{Ufixed8}, r4) --> RGB{Ufixed8}(1,0,0)
        @fact convert(ColorTypes.RGB4{Ufixed8}, r4) --> ColorTypes.RGB4{Ufixed8}(1,0,0)
        @fact convert(ColorTypes.RGB4{Float32}, r4) --> ColorTypes.RGB4{Float32}(1,0,0)
        @fact convert(ColorTypes.BGR{Float32}, r4) --> ColorTypes.BGR{Float32}(1,0,0)

        iob = IOBuffer()
        c = ColorTypes.BGR{Ufixed8}(0.1,0.2,0.301)
        show(iob, c)
        @fact takebuf_string(iob) --> "BGR{Ufixed8}(0.102,0.2,0.302)"
        c = ColorTypes.Gray{Ufixed16}(0.8)
        show(iob, c)
        @fact takebuf_string(iob) --> "Gray{Ufixed16}(0.8)"
        ca = AlphaColorValue(c, Ufixed16(0.2))
        show(iob, ca)
        @fact takebuf_string(iob) --> "GrayAlpha{Ufixed16}(0.8,0.2)"
    end

    context("YIQ") do
        @fact convert(ColorTypes.YIQ, RGB(1,0,0)) --> ColorTypes.YIQ(0.299, 0.595716, 0.211456)
        @fact convert(ColorTypes.YIQ, RGB(0,1,0)) --> ColorTypes.YIQ(0.587, -0.274453, -0.522591)
        @fact convert(ColorTypes.YIQ, RGB(0,0,1)) --> ColorTypes.YIQ(0.114, -0.321263, 0.311135)
        @fact convert(RGB, ColorTypes.YIQ(1.0,0.0,0.0)) --> RGB(1,1,1)
        v = 0.5957
        @fact convert(RGB, ColorTypes.YIQ(0.0,1.0,0.0)) --> RGB(0.9563*v,-0.2721*v,-1.1070*v)  # will be clamped to v
        v = -0.5226
        @fact convert(RGB, ColorTypes.YIQ(0.0,0.0,-1.0)) --> RGB(0.6210*v,-0.6474*v,1.7046*v)  # will be clamped to v
    end

end


text:
# Some utils for DG issues
export dg_nodes_element_to_face, dg_nodes_face_to_element

# For face integrals, change a vector from full element to face components only
# and the other way around just inserting zeros.
function dg_nodes_element_to_face(v, refel, face)
    nfp = refel.Nfp[face];
    vf = zeros(nfp);

    for i=1:nfp
        vf[i] = v[refel.face2local[face][i]];
    end

    return vf;
end

function dg_nodes_face_to_element(vf, refel, face)
    nfp = refel.Nfp[face];
    v = zeros(refel.Np);

    for i=1:nfp
        v[refel.face2local[face][i]] = vf[i];
    end

    return v;
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Printf
using BenchmarkTools
using FFTW
using CUDA
using Oceananigans
using Benchmarks

# Benchmark functions

function benchmark_fft(::Type{CPU}, N, dims; FT=Float64, planner_flag=FFTW.PATIENT)
    A = zeros(complex(FT), N, N, N)
    FFT! = FFTW.plan_fft!(A, dims, flags=planner_flag)
    trial = @benchmark ($FFT! * $A) samples=10
    return trial
end

function benchmark_fft(::Type{GPU}, N, dims; FT=Float64, planner_flag=FFTW.PATIENT)
    A = zeros(complex(FT), N, N, N) |> CuArray

    # Cannot do CUDA FFTs along non-batched dims so dim=2 must
    # be computed via tranposing.
    if dims == 2
        B = similar(A)
        FFT! = CUDA.CUFFT.plan_fft!(A, 1)

        trial = @benchmark begin
            CUDA.@sync begin
                permutedims!($B, $A, (2, 1, 3))
                $FFT! * $B
                permutedims!($A, $B, (2, 1, 3))
            end
        end samples=10
    else
        FFT! = CUDA.CUFFT.plan_fft!(A, dims)

        trial = @benchmark begin
            CUDA.@sync ($FFT! * $A)
        end samples=10
    end

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Ns = [16, 64, 256]
dims = [1, 2, 3, (1, 2, 3)]

# Run and summarize benchmarks

print_system_info()
suite = run_benchmarks(benchmark_fft; Architectures, Ns, dims)

df = benchmarks_dataframe(suite)
sort!(df, [:Architectures, :dims, :Ns], by=(string, string, identity))
benchmarks_pretty_table(df, title="FFT benchmarks")

println("3D FFT --> 3 × 1D FFTs slowdown:")
for Arch in Architectures, N in Ns
    fft_x_time = median(suite[(Arch, N, 1)]).time
    fft_y_time = median(suite[(Arch, N, 2)]).time
    fft_z_time = median(suite[(Arch, N, 3)]).time
    fft_xyz_time = median(suite[(Arch, N, (1, 2, 3))]).time
    slowdown = (fft_x_time + fft_y_time + fft_z_time) / fft_xyz_time
    @printf("%s, %3d: %.4fx\n", Arch, N, slowdown)
end


text:
struct Duffy <: QuadratureName end

const duffy = Duffy()

function Quadrature(p::Polytope,name::Duffy,degree::Integer)
  @assert is_simplex(p)
  D = num_dims(p)
  x,w = _duffy_quad_data(degree,D)
  msg = "Simplex quadrature of degree $degree obtained with Duffy transformation and tensor product of 1d Gauss-Jacobi and Gauss-Legendre rules."
  GenericQuadrature(x,w,msg)
end

function _duffy_quad_data(order::Integer,D::Int)

  beta = 0
  dim_to_quad_1d = [
    _gauss_jacobi_in_0_to_1(order,(D-1)-(d-1),beta) for d in 1:(D-1) ]

  quad_1d = _gauss_legendre_in_0_to_1(order)
  push!(dim_to_quad_1d,quad_1d)

  x_pos = 1
  w_pos = 2
  dim_to_xs_1d = [quad_1d[x_pos] for quad_1d in dim_to_quad_1d]
  dim_to_ws_1d = [quad_1d[w_pos] for quad_1d in dim_to_quad_1d]

  a = 0.5
  for d in (D-1):-1:1
    ws_1d = dim_to_ws_1d[d]
    ws_1d[:] *= a
    a *= 0.5
  end

  x,w = _tensor_product_duffy(dim_to_xs_1d,dim_to_ws_1d)

  (_duffy_map.(x),w)

end

# Duffy map from the n-cube in [0,1]^d to the n-simplex in [0,1]^d
function _duffy_map(q::Point{D,T}) where {D,T}
  m = zero(Mutable(Point{D,T}))
  m[1] = q[1]
  a = one(T)
  for i in 2:D
    a *= (1-q[i-1])
    m[i] = a*q[i]
  end
  Point{D,T}(m)
end

_duffy_map(q::Point{1,T}) where T = q

function _gauss_jacobi_in_0_to_1(order,alpha,beta)
  n = _npoints_from_order(order)
  x,w = gaussjacobi(n,alpha,beta)
  _map_to(0,1,x,w)
end

function _gauss_legendre_in_0_to_1(order)
  n = _npoints_from_order(order)
  x,w = gausslegendre(n)
  _map_to(0,1,x,w)
end

# Transforms a 1-D quadrature from `[-1,1]` to `[a,b]`, with `a<b`.
function _map_to(a,b,points,weights)
  points_ab = 0.5*(b-a)*points .+ 0.5*(a+b)
  weights_ab = 0.5*(b-a)*weights
  (points_ab, weights_ab)
end

function _npoints_from_order(order)
  ceil(Int, (order + 1.0) / 2.0 )
end

function _tensor_product_duffy(
  dim_to_xs_1d::Vector{Vector{T}},
  dim_to_ws_1d::Vector{Vector{W}}) where {T,W}

  D = length(dim_to_ws_1d)
  @assert D == length(dim_to_xs_1d)
  dim_to_n = [length(ws_1d) for ws_1d in dim_to_ws_1d]
  n = prod(dim_to_n)
  xs = zeros(Point{D,T},n)
  ws = zeros(W,n)
  cis = CartesianIndices(tuple(dim_to_n...))
  m = zero(Mutable(Point{D,T}))
  _tensor_product_duffy!(xs,ws,dim_to_xs_1d,dim_to_ws_1d,cis,m)
  (xs,ws)
end

function _tensor_product_duffy!(
  xs,ws,dim_to_xs_1d,dim_to_ws_1d,cis::CartesianIndices{D},m) where D
  k = 1
  for ci in cis
    w = 1.0
    for d in 1:D
      xs_1d = dim_to_xs_1d[d]
      ws_1d = dim_to_ws_1d[d]
      i = ci[d]
      xi = xs_1d[i]
      wi = ws_1d[i]
      w *= wi
      m[d] = xi
    end
    xs[k] = m
    ws[k] = w
    k += 1
  end
end


text:
module Visuals

    using PlotlyJS

    import ..Utils: x, y, z

    export plot_nodes, plot_curve, plot_surface

    function plot_nodes(
        nodes::AbstractArray;
        color::String="black",
        name="nodes",
        node_size::Number=8,
    )
        N = ndims(nodes)
        if size(nodes, 1) == 3
            _type = "scatter3d"
        else
            _type = "scatter"
        end

        scatter(
            x= (N==2) ? x(nodes) : vec(x(nodes)),
            y= (N==2) ? y(nodes) : vec(y(nodes)),
            z= (N==2) ? z(nodes) : vec(z(nodes)),
            mode="markers",
            type=_type,
            marker=attr(
                    size=node_size,
                    color=color,
                    opacity=1,
                ),
            name=name
            )
    end

    """
        plot_curve(
            curve::Curve;
            as_scatter::Bool=false,
            color::String="black",
            width::Number=10,
            )

    Creates a scatter or line representation of a 1
    dimensional curve in N dimensional space
    """
    function plot_curve(
            points;
            as_scatter::Bool=false,
            color::String="black",
            width::Number=10,
            name::String="curve"
            )

            # prep params
            N = size(points, 1)
            if N == 3
                _type = "scatter3d"
            else
                _type = "scatter"
            end

            if as_scatter
                mode= "marker"
            else
                mode = "lines"
            end

            # plot curve
            scatter(
                x=x(points),
                y=y(points),
                z=z(points),
                mode=mode,
                type=_type,
                line = attr(color=color, width=width),
                name=name
            )
    end

    function plot_surface(points; with_points::Bool=true, name::String="surface", opacity=1.0, kwargs...)
        srf = surface(
            x=x(points),
            y=y(points),
            z=z(points),
            name=name,
            showscale=false,
            opacity=opacity,
            kwargs...
            )

        if with_points == false
            return srf
        end

        sct = scatter(
            x=vec(x(points)),
            y=vec(y(points)),
            z=vec(z(points)),
            mode="markers",
            type="scatter3d",
            marker=attr(
                    size=1,
                    color="black",
                ),
            name=name,
            )

        return [srf, sct]

    end
end

max_stars_repo_name: FedeClaudi/Splines

max_issues_repo_name: FedeClaudi/Splines

max_forks_repo_name: FedeClaudi/Splines


text:
mutable struct DirectJumpAggregation{T,S,F1,F2,RNG} <: AbstractSSAJumpAggregator
  next_jump::Int
  prev_jump::Int
  next_jump_time::T
  end_time::T
  cur_rates::Vector{T}
  sum_rate::T
  ma_jumps::S
  rates::F1
  affects!::F2
  save_positions::Tuple{Bool,Bool}
  rng::RNG
end
DirectJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; kwargs...) where {T,S,F1,F2,RNG} =
  DirectJumpAggregation{T,S,F1,F2,RNG}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng)

############################# Required Functions #############################

# creating the JumpAggregation structure (tuple-based constant jumps)
function aggregate(aggregator::Direct, u, p, t, end_time, constant_jumps,
                    ma_jumps, save_positions, rng; kwargs...)

  # handle constant jumps using tuples
  rates, affects! = get_jump_info_tuples(constant_jumps)

  build_jump_aggregation(DirectJumpAggregation, u, p, t, end_time, ma_jumps,
                          rates, affects!, save_positions, rng; kwargs...)
end

# creating the JumpAggregation structure (function wrapper-based constant jumps)
function aggregate(aggregator::DirectFW, u, p, t, end_time, constant_jumps,
                    ma_jumps, save_positions, rng; kwargs...)

  # handle constant jumps using function wrappers
  rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)

  build_jump_aggregation(DirectJumpAggregation, u, p, t, end_time, ma_jumps,
                          rates, affects!, save_positions, rng; kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::DirectJumpAggregation, integrator, u, params, t)
  generate_jumps!(p, integrator, u, params, t)
  nothing
end

# execute one jump, changing the system state
@inline function execute_jumps!(p::DirectJumpAggregation, integrator, u, params, t)
  update_state!(p, integrator, u)
  nothing
end

# calculate the next jump / jump time
function generate_jumps!(p::DirectJumpAggregation, integrator, u, params, t)
  p.sum_rate, ttnj = time_to_next_jump(p, u, params, t)
  @fastmath p.next_jump_time = t + ttnj
  @inbounds p.next_jump = searchsortedfirst(p.cur_rates, rand(p.rng) * p.sum_rate)
  nothing
end

######################## SSA specific helper routines ########################

# tuple-based constant jumps
@fastmath function time_to_next_jump(p::DirectJumpAggregation{T,S,F1,F2,RNG}, u, params, t) where {T,S,F1 <: Tuple, F2 <: Tuple, RNG}
  prev_rate = zero(t)
  new_rate  = zero(t)
  cur_rates = p.cur_rates

  # mass action rates
  majumps   = p.ma_jumps
  idx       = get_num_majumps(majumps)
  @inbounds for i in 1:idx
    new_rate     = evalrxrate(u, i, majumps)
    cur_rates[i] = new_rate + prev_rate
    prev_rate    = cur_rates[i]
  end

  # constant jump rates
  rates = p.rates
  if !isempty(rates)
    idx  += 1
    fill_cur_rates(u, params, t, cur_rates, idx, rates...)
    @inbounds for i in idx:length(cur_rates)
      cur_rates[i] = cur_rates[i] + prev_rate
      prev_rate    = cur_rates[i]
    end
  end

  @inbounds sum_rate = cur_rates[end]
  sum_rate, randexp(p.rng) / sum_rate
end

@inline function fill_cur_rates(u, p, t, cur_rates, idx, rate, rates...)
  @inbounds cur_rates[idx] = rate(u, p, t)
  idx += 1
  fill_cur_rates(u, p, t, cur_rates, idx, rates...)
end

@inline function fill_cur_rates(u, p, t, cur_rates, idx, rate)
  @inbounds cur_rates[idx] = rate(u, p, t)
  nothing
end

# function wrapper-based constant jumps
@fastmath function time_to_next_jump(p::DirectJumpAggregation{T,S,F1,F2,RNG}, u, params, t) where {T,S,F1 <: AbstractArray,F2 <: AbstractArray, RNG}
  prev_rate = zero(t)
  new_rate  = zero(t)
  cur_rates = p.cur_rates

  # mass action rates
  majumps   = p.ma_jumps
  idx       = get_num_majumps(majumps)
  @inbounds for i in 1:idx
    new_rate     = evalrxrate(u, i, majumps)
    cur_rates[i] = new_rate + prev_rate
    prev_rate    = cur_rates[i]
  end

  # constant jump rates
  idx  += 1
  rates = p.rates
  @inbounds for i in 1:length(p.rates)
    new_rate       = rates[i](u, params, t)
    cur_rates[idx] = new_rate + prev_rate
    prev_rate      = cur_rates[idx]
    idx           += 1
  end

  @inbounds sum_rate = cur_rates[end]
  sum_rate, randexp(p.rng) / sum_rate
end


text:
## Rotators Our rotators have field names c, s where c nad s are
## either T or Complex{T} It proved to be faster to have immutable
## rotators, rather than mutable ones so there are no "setter" methods
## anymore, a new instance (via `Rotator(c,s,i)`) should be used

abstract type CoreTransform{T,S} end
abstract type AbstractRotator{T,S} <: CoreTransform{T,S} end

## our rotators are [c s; -s conj(c)]
## We *could* just directly use `LinearAlgebra.Givens` here, but
## instead implement Rotator below to ensure the types Real/Real and Complex/Real for c and s
## XXX: for twisted, it would have been easier to have Complex/Complex rotators, as then
## `passthrough_phase` could have been dramatically simplified.
##
## Interface
## vals: return c,s
## idx: return i
## adjoint: return adjoint

## get values
@inline vals(r::AbstractRotator{T}) where {T} = (r.c, r.s)
@inline idx(r::AbstractRotator) = r.i

Base.copy(a::AbstractRotator) = AbstractRotator(a.c, a.s, a.i)
is_diagonal(r::AbstractRotator{T,S}) where {T,S} = norm(r.s) <= eps(T)

## Multiplication of a rotator and a matrix
## this uses [c s; -conj(s) conj(c)] for rotator!
function *(a::AbstractRotator, M::AbstractArray)
    N = copy(M)
    lmul!(a, N)
    N
end

function LinearAlgebra.lmul!(a::AbstractRotator, M::AbstractArray)
    c, s = vals(a)
    i = idx(a); j = i+1
    n = size(M)[2]
    for k in 1:n
        mik, mjk = M[i,k], M[j,k]
        M[i,k]  =  c * mik + s * mjk
        M[j,k]  =  -conj(s) * mik + conj(c) * mjk
    end
    M
end

function *(M::AbstractArray, a::AbstractRotator)
    N = copy(M)
    rmul!(N, a)
    N
end

function LinearAlgebra.rmul!(M::AbstractArray, a::AbstractRotator)
    c, s = vals(a)
    i = idx(a); j = i+1
    n = size(M)[1]
    for k in 1:n
        mki, mkj = M[k,i], M[k,j]
        M[k, i] = c * mki - conj(s) * mkj
        M[k, j] = s * mki + conj(c) * mkj
    end
    M
end

*(Qs::Vector{R}, M::Matrix) where {R <: CoreTransform} = foldr(*, Qs, init=M)
*(M::Matrix, Qs::Vector{R}) where {R <: CoreTransform} = foldl(*, Qs, init=M)
function LinearAlgebra.lmul!(Qs::Vector{R}, M::Matrix) where {R <: CoreTransform}
    foldr(lmul!, Qs, init=M)
end
function LinearAlgebra.rmul!(M::Matrix, Qs::Vector{R}) where {R <: CoreTransform}
    foldl(rmul!, Qs, init=M)
end

# A rotator
# c is real or complex
# s is always real. (The Complex/Complex case was dropped)
# i the index, adds storage space, but simplifies some algorithms
struct Rotator{T,S} <: AbstractRotator{T,S}
c::S
s::T
i::Int
end

function LinearAlgebra.adjoint(U::Rotator)
    c,s = vals(U)
    Rotator(conj(c), -conj(s), idx(U))
end

Base.eltype(::Type{Rotator{T,S}}) where {T, S} = S
Base.one(::Type{Rotator{T,S}}) where {T, S} = Rotator(one(S), zero(T), 0)
Base.copy(U::Rotator) = Rotator(vals(U)..., idx(U))

## A diagonal rotator has iszero(s) being true
struct DiagonalRotator{T,S} <: AbstractRotator{S, S}
c::S
i::Int
DiagonalRotator{T,S}(c,i) where {T,S} = new(convert(S,c),i)
DiagonalRotator{T,S}(c,s,i) where {T,S} = new(convert(S,c),i)
DiagonalRotator{S}(c,i) where {S} = DiagonalRotator{real(S),S}(c,i)
DiagonalRotator(c::S, i) where {S} = DiagonalRotator{real(S),S}(c,i)
end

vals(D::DiagonalRotator{T,S}) where {T,S} = D.c, zero(T)
LinearAlgebra.adjoint(U::DiagonalRotator) = DiagaonalRotator(conj(U.c), U.i)

## for real case, we have identity diagonal
struct IdentityRotator{T,S} <: AbstractRotator{T,S}
i::Int
IdentityRotator{T,S}(i) where {T,S} = new(i)
end

vals(D::IdentityRotator{T,S}) where {T,S} = (one(S), zero(T))
idx(D::IdentityRotator{T,S}) where {T,S} = D.i
LinearAlgebra.adjoint(U::IdentityRotator) = U


text:
export End

"""
    End()

`End` objects `show` a FastHenry .end command.
"""
mutable struct End <: Element
end


text:
@test typeof(@inferred(NoOp())) <: NoOp <: Augmentor.AffineOperation

@testset "constructor" begin
    @test str_show(NoOp()) == "Augmentor.NoOp()"
    @test str_showconst(NoOp()) == "NoOp()"
    @test str_showcompact(NoOp()) == "No operation"
end
@testset "eager" begin
    @test_throws MethodError Augmentor.applyeager(NoOp(), nothing)
    @test Augmentor.supports_eager(NoOp) === false
    res1 = rect
    res2 = OffsetArray(rect, -2, -1)
    res3 = OffsetArray(rect, 0, 0)
    imgs = [
        (rect, res1),
        (view(rect, :, :), res1),
        (Augmentor.prepareaffine(rect), res3),
        (OffsetArray(rect, -2, -1), res2),
        (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res3),
    ]
    @testset "single image" begin
        for (img_in, img_out) in imgs
            res = @inferred(Augmentor.applyeager(NoOp(), img_in))
            @test res == img_out
            @test typeof(res) == typeof(img_out)
        end
    end
end
@testset "affine" begin
    @test_throws MethodError Augmentor.toaffinemap(NoOp(), nothing)
    @test Augmentor.supports_affine(NoOp) === true
    @test @inferred(Augmentor.toaffinemap(NoOp(), rect)) == AffineMap(@SMatrix([1. 0; 0 1]), @SVector([0., 0.]))
    wv = @inferred Augmentor.applyaffine(NoOp(), rect)
    @test wv == rect
    @test typeof(wv) <: InvWarpedView{eltype(rect),2}
end
@testset "affineview" begin
    @test Augmentor.supports_affineview(NoOp) === true
    wv = @inferred Augmentor.applyaffineview(NoOp(), rect)
    @test typeof(wv) <: SubArray{eltype(rect),2}
    @test typeof(parent(wv)) <: InvWarpedView
    @test parent(parent(wv)) === rect
    @test wv == rect
end
@testset "lazy" begin
    @test Augmentor.supports_lazy(NoOp) === true
    @test @inferred(Augmentor.applylazy(NoOp(), rect)) === rect
    oa = OffsetArray(rect, -2, -1)
    @test @inferred(Augmentor.applylazy(NoOp(), oa)) === oa
    wv = Augmentor.prepareaffine(rect)
    @test @inferred(Augmentor.applylazy(NoOp(), wv)) === wv
    res1 = rect
    res2 = OffsetArray(rect, -2, -1)
    res3 = OffsetArray(rect, 0, 0)
    imgs = [
        (rect, res1),
        (view(rect, :, :), res1),
        (Augmentor.prepareaffine(rect), res3),
        (OffsetArray(rect, -2, -1), res2),
        (view(rect, IdentityRange(1:2), IdentityRange(1:3)), res3),
    ]
    @testset "single image" begin
        for (img_in, img_out) in imgs
            res = @inferred(Augmentor.applylazy(NoOp(), img_in))
            @test res == img_out
        end
    end
end
@testset "view" begin
    @test Augmentor.supports_view(NoOp) === true
    @test @inferred(Augmentor.applyview(NoOp(), rect)) === view(rect, IdentityRange(1:2), IdentityRange(1:3))
end
@testset "stepview" begin
    @test Augmentor.supports_stepview(NoOp) === true
    @test @inferred(Augmentor.applystepview(NoOp(), rect)) === view(rect, 1:1:2, 1:1:3)
end
@testset "permute" begin
    @test Augmentor.supports_permute(NoOp) === false
end


text:
module Simulation

# Write your package code here.
    # install necessary packages (open source)
    using Statistics, Plots, Distributions, FFTW, LinearAlgebra

    # install necessary packages (original)
    using .FundamentalStats

    # function name
    #OK
    #export SpectraPlot
    export BarcodePlot
    export BarcodePlot!
    export BinningPlot
    export BinningPlot!
    export Correlogram
    export Correlogram!
    export VariogramPlot
    export VariogramPlot!
    export RCDFPlot
    export RCDFPlot!
    export ECDFPlot
    export ECDFPlot!
end

max_stars_repo_name: yukki-jpn/JuliaTSA

max_issues_repo_name: yukki-jpn/JuliaTSA

max_forks_repo_name: yukki-jpn/JuliaTSA


text:
using ToyFHE
using OffsetArrays
using Test
using Primes
using GaloisFields

N = 2^5
q₀ = nextprime(Int128(2)^40 + 1, 1; interval=2N)
q₁ = nextprime(q₀ + 2N, 1; interval=2N)
ps = nextprime(q₁ + 2N, 1; interval=2N)

ℛ = let CT = ToyFHE.CRTEncoded{3, Tuple{GaloisField.((q₀,q₁,ps))...}}
    ζ₂n = GaloisFields.minimal_primitive_root(CT, 2N)
    ToyFHE.NegacyclicRing{CT, N}(ζ₂n)
end

scale = 2^60
Tscale = FixedRational{scale}

plain = CKKSEncoding{Tscale}(zero(ℛ))
plain .= 2

Tscale2 = FixedRational{scale/ps}
@test isapprox(CKKSEncoding{Tscale2}(ToyFHE.modswitch(convert(ToyFHE.NTT.RingElement, plain)))[0], 2.0; atol=10^-5)

# The same but with encryption noise
ckks_params = CKKSParams(ℛ, 1, 3.2)
kp = keygen(ckks_params)

switched = let c = ToyFHE.modswitch(encrypt(kp, plain))
    decrypt(kp, c)
end
@test switched ≈ plain atol=10^-3


text:
@testset "Exterior map" begin
  z1 = -1.0+0.0im
  z2 = exp(im*5π/4)
  sing1 = 2
  beta = [0.5,0.5,0.5,0.5]
  z = [1.0im,-1.0,-1.0im,1.0]
  nqpts = 8
  qdat = SchwarzChristoffel.Exterior.Integration.qdata(beta,nqpts)
  dquad = SchwarzChristoffel.Exterior.DQuad(beta,qdat)
  @test dquad([z1],[z2],[sing1],z)[1] ≈ -0.599070117 - 0.599070117im
  dabsquad = SchwarzChristoffel.Exterior.DabsQuad(beta,qdat)
  @test dabsquad([z1],[z2],[sing1],z)[1] ≈ 0.847213085

  p = SchwarzChristoffel.Polygon([-1.0,0.2,1.0,-1.0],[-1.0,-1.0,0.5,1.0])
  m = SchwarzChristoffel.ExteriorMap(p)
  dm = SchwarzChristoffel.DerivativeMap(m)
  m⁻¹ = SchwarzChristoffel.InverseMap(m)
  prev, C = SchwarzChristoffel.parameters(m)
  @test prev[3] ≈ -0.902383088-0.430934755im

  zeta = [0.1,0.5-0.75im,-0.25-0.3im]
  z = m(zeta;inside=true)
  @test z ≈
        ComplexF64[-6.934402676 - 7.689645145im,
                    0.043977374 - 1.112493651im,
                    2.411813129 - 0.044778980im]
  @test m⁻¹(z;inside=true) ≈ zeta

  (dz,ddz) = dm(zeta;inside=true)
  @test dz ≈
        ComplexF64[67.206798744 + 76.628383501im,
                   -1.116655339 + 0.5445759699im,
                    3.991294815 - 5.3064069537im]
  zeta = 0.5-0.75im
  z = m(zeta;inside=true)
  @test z ≈ 0.043977374 - 1.112493651im
  @test m⁻¹(z;inside=true) ≈ zeta
  (dz,ddz) = dm(zeta;inside=true)
  @test dz ≈ -1.116655339 + 0.5445759699im
  zeta = [1.0+3.0im,-2.0-2.0im,0.0+1.1im]
  z = m(zeta)
  @test z ≈
        ComplexF64[ 0.816139620 + 3.029559043im,
                   -2.252366325 - 2.085230469im,
                   -0.333104076 + 0.975837123im]
  @test m⁻¹(z) ≈ zeta
  (dz,ddz) = dm(zeta)
  @test dz ≈
        ComplexF64[ 1.030528002 + 0.004449924im,
                    1.006955879 - 0.011501136im,
                    1.300781073 - 0.266624652im]
  zeta = 1.0+3.0im
  z = m(zeta)
  @test z ≈ 0.816139620 + 3.029559043im
  @test m⁻¹(z) ≈ zeta
  dz, ddz = dm(zeta)
  @test dz ≈ 1.030528002 + 0.004449924im

  ccoeff, dcoeff = SchwarzChristoffel.coefficients(m)
  @test ccoeff[1] ≈ 1.019789030 - 4.043787137e-8im

end

@testset "Power series map" begin

  a = 1
  b = 1
  c = ComplexF64[0.5(a+b),0,0.5(a-b)]
  m = SchwarzChristoffel.PowerMap(c)
  @test area(m) ≈    3.141592654
  @test centroid(m) ≈ ComplexF64(0)
  @test Jmoment(m) ≈ 1.570796327
  a = 1
  b = 0
  c = ComplexF64[0.5(a+b),0,0.5(a-b)]
  m = SchwarzChristoffel.PowerMap(c)
  @test area(m) ≈    Float64(0)
  @test centroid(m) ≈ ComplexF64(0)
  @test Jmoment(m) ≈ Float64(0)

  m = SchwarzChristoffel.PowerMap(ComplexF64[1,0,0,0,0,0.1])
  @test area(m) ≈ 3.015928947
  @test centroid(m) ≈ ComplexF64(0)
  @test Jmoment(m) ≈ 1.475920229

  a = 1
  b = 0.1
  c = ComplexF64[0.5(a+b),0,0.5(a-b)]
  m = SchwarzChristoffel.PowerMap(c)
  m⁻¹ = SchwarzChristoffel.InverseMap(m)
  dm = SchwarzChristoffel.DerivativeMap(m)
  zeta = [1.0+3.0im,-2.0-2.0im,0.0+1.1im]
  z = m(zeta)
  @test z ≈
      ComplexF64[0.595+1.515im,-1.2125-0.9875im,
                 0.195909090909im]
  @test m⁻¹(z) ≈ zeta
  dz, ddz = dm(zeta)
  @test dz ≈ ComplexF64[0.586 + 0.027im,
                        0.55 + 0.05625im,
                        0.921900826]
  zeta = 1.0+3.0im
  z = m(zeta)
  @test z ≈ 0.595+1.515im
  @test m⁻¹(z) ≈ zeta
  dz, ddz = dm(zeta)
  @test dz ≈ 0.586 + 0.027im

  a = 1
  b = 0
  c = ComplexF64[0.5(a+b),0,0.5(a-b)]
  m = SchwarzChristoffel.PowerMap(c)
  m⁻¹ = SchwarzChristoffel.InverseMap(m)
  dm = SchwarzChristoffel.DerivativeMap(m)
  z = [-1.0+0im,-1.1+0im,1.0+0im,1.1+0im]
  @test m(m⁻¹(z)) ≈ z

end

@testset "Karman-Trefftz map" begin
  ν = 1.94; ϵ = 0.08*sqrt(2); δ = 3π/4; C = 0.25;
  m = SchwarzChristoffel.KarmanTrefftzMap(ν,ϵ,δ,C)
  @test area(m) ≈ 0.08510380847113
  @test centroid(m) ≈ -0.0796138968044 + 0.0271627776676im
  @test Jmoment(m) ≈ 0.00538941353392

  zeta = [1.0+3.0im,-2.0-2.0im,0.0+1.1im]
  z = m(zeta)
  @test all(abs.(z - ComplexF64[0.269784+0.768849im,
                      -0.616452-0.470251im,
                      -0.0310697+0.138951im]) .≤ 1e-6)

  m⁻¹ = SchwarzChristoffel.InverseMap(m)
  @test m⁻¹(m(zeta)) ≈ zeta

end


text:
@assert length(ARGS) == 3
file_loc = joinpath(dirname(@__FILE__), "join_performance.jl")
llen = ARGS[1]
rlen = ARGS[2]
type = ARGS[3]

# re-enable "cat" when
# https://github.com/JuliaData/CategoricalArrays.jl/pull/327
# and
# https://github.com/JuliaData/CategoricalArrays.jl/pull/331
# are merged and tagged
for a3 in ["str", "int", "pool"], # "cat"],
    a4 in ["uniq", "dup", "manydup"],
    a5 in ["sort", "rand"],
    a6 in ["1", "2"]
    a4 == "uniq" && a3 in ["pool", "cat"] && continue
    run(`julia $file_loc $llen $rlen $a3 $a4 $a5 $a6 $type`)
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# 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.
# *****************************************************************************

"""
Simulated Annealing for the ACASX problem.
Example usage: config=configure(ACASX_SA,"normal","nvn_dasc"); acasx_sa(;config...)
"""
module ACASX_SA

export configure, acasx_sa, acasx_sa1, acasx_temp_params

import Compat.ASCIIString
using ExprSearch.SA
using Reexport

using GrammarExpts
using ACASXProblem, DerivTreeVis, SA_Logs
using RLESUtils, Configure
import Configure.configure

const CONFIGDIR = joinpath(dirname(@__FILE__), "..", "config")
const RESULTDIR = joinpath(dirname(@__FILE__), "..", "..", "..", "results")

configure(::Type{Val{:ACASX_SA}}, configs::AbstractString...) = configure_path(CONFIGDIR, configs...)

"""
Example call:
config=configure(ACASX_SA, "nvn_dasc", "normal")
acasx_sa(; config...)
"""
function acasx_sa(;outdir::AbstractString=joinpath(RESULTDIR, "ACASX_SA"),
                  seed=1,
                  logfileroot::AbstractString="acasx_sa_log",

                  runtype::Symbol=:nmacs_vs_nonnmacs,
                  data::AbstractString="dasc",
                  manuals::AbstractString="dasc_manual",
                  clusterdataname::AbstractString="josh1",

                  maxsteps::Int64=20,
                  T1::Float64=12.184,
                  alpha::Float64=0.99976,
                  n_epochs::Int64=50,
                  n_starts::Int64=1,
                  n_threads::Int64=1,

                  loginterval::Int64=100,
                  vis::Bool=true)

  srand(seed)
  mkpath(outdir)

  problem = ACASXClustering(runtype, data, manuals, clusterdataname)

  observer = Observer()
  par_observer = Observer()

  logs = default_logs(par_observer)
  default_console!(observer, loginterval)

  sa_params = SAESParams(maxsteps, T1, alpha, n_epochs, n_starts, observer)
  psa_params = PSAESParams(n_threads, sa_params, par_observer)

  result = exprsearch(psa_params, problem)

  add_members_to_log!(logs, problem, result.expr)
  outfile = joinpath(outdir, "$(logfileroot).txt")
  save_log(outfile, logs)

  if vis
    derivtreevis(get_derivtree(result), joinpath(outdir, "$(logfileroot)_derivtreevis"))
  end

  return result
end

function acasx_sa1(;outdir::AbstractString=joinpath(RESULTDIR, "ACASX_SA1"),
                  seed=1,
                  logfileroot::AbstractString="acasx_sa_log",

                  runtype::Symbol=:nmacs_vs_nonnmacs,
                  data::AbstractString="dasc",
                  manuals::AbstractString="dasc_manual",
                  clusterdataname::AbstractString="josh1",

                  maxsteps::Int64=20,
                  T1::Float64=12.184,
                  alpha::Float64=0.99976,
                  n_epochs::Int64=50,
                  n_starts::Int64=1,

                  loginterval::Int64=100,
                  vis::Bool=true)

  srand(seed)
  mkpath(outdir)

  problem = ACASXClustering(runtype, data, manuals, clusterdataname)

  observer = Observer()

  logs = default_logs1(observer, loginterval)
  default_console!(observer, loginterval)

  sa_params = SAESParams(maxsteps, T1, alpha, n_epochs, n_starts, observer)
  result = exprsearch(sa_params, problem)

  add_members_to_log!(logs, problem, result.expr)
  outfile = joinpath(outdir, "$(logfileroot).txt")
  save_log(outfile, logs)

  if vis
    derivtreevis(get_derivtree(result), joinpath(outdir, "$(logfileroot)_derivtreevis"))
  end

  return result
end

function add_members_to_log!{T}(logs::TaggedDFLogger, problem::ACASXClustering{T}, expr)
  add_folder!(logs, "members", [ASCIIString, ASCIIString], ["members_true", "members_false"])
  members_true, members_false = get_members(problem, expr)
  push!(logs, "members", [join(members_true, ","), join(members_false, ",")])
end

"""
Get recommended temperature commands
"""
function acasx_temp_params(P1::Float64=0.8; seed=1,
                           n_epochs::Int64=1,
                           Tfinal::Float64=1.0,
                           runtype::Symbol=:nmacs_vs_nonnmacs,
                           data::AbstractString="dasc",
                           manuals::AbstractString="dasc_manual",
                           clusterdataname::AbstractString="josh1",
                           maxsteps::Int64=20,
                           N::Int64=1000,
                           ntrials::Int64=10)

  srand(seed)

  problem = ACASXClustering(runtype, data, clusterdataname)
  T1, alpha, n_epochs = estimate_temp_params(problem, P1, n_epochs, Tfinal, maxsteps, N, ntrials)

  return T1, alpha, n_epochs
end

end #module


text:
THArray{T}{T,N}(size::NTuple{N,Integer}) = THArray{T,N}(size)
THArray{T}{T,N}(size::Vararg{Integer,N}) = THArray{T}(size)

THArray(args...) = THArray{Float64}(args...)

THArray{T}{T}(data::AbstractArray) = copy!(THArray{T}(size(data)), data)

THArray(data::AbstractArray) = THArray{eltype(data)}(data)

Base.size(xs::THArray) = ntuple(i -> size(xs, i), ndims(xs))

Base.similar(xs::THArray) = typeof(xs)(size(xs))

# TODO: promotion for THArrays of different types

import Base: +, -

(xs::THArray{T} + ys::THArray{T}) where T = add!(xs, ys, out = similar(xs))
(xs::THArray{T} - ys::THArray{T}) where T = sub!(xs, ys, out = similar(xs))


text:
function FieldLatitudeCrossingDetector(arg0::Field, arg1::OneAxisEllipsoid, arg2::jdouble)
    return FieldLatitudeCrossingDetector((Field, OneAxisEllipsoid, jdouble), arg0, arg1, arg2)
end

function FieldLatitudeCrossingDetector(arg0::RealFieldElement, arg1::RealFieldElement, arg2::OneAxisEllipsoid, arg3::jdouble)
    return FieldLatitudeCrossingDetector((RealFieldElement, RealFieldElement, OneAxisEllipsoid, jdouble), arg0, arg1, arg2, arg3)
end

function g(obj::FieldLatitudeCrossingDetector, arg0::FieldSpacecraftState)
    return jcall(obj, "g", RealFieldElement, (FieldSpacecraftState,), arg0)
end

function get_body(obj::FieldLatitudeCrossingDetector)
    return jcall(obj, "getBody", OneAxisEllipsoid, ())
end

function get_latitude(obj::FieldLatitudeCrossingDetector)
    return jcall(obj, "getLatitude", jdouble, ())
end


text:
function update_automerge_comment!(repo::GitHub.Repo,
                                   pr::GitHub.PullRequest;
                                   body::AbstractString,
                                   auth::GitHub.Authorization,
                                   whoami)::Nothing
    my_comments = my_retry(() -> get_all_my_pull_request_comments(repo,
                                                                  pr;
                                                                  auth = auth,
                                                                  whoami = whoami))
    num_comments = length(my_comments)
    _body::String = convert(String,
                                    strip(string(strip(body),
                                    "\n",
                                    "<!---\n",
                                    "this_is_the_single_automerge_comment\n",
                                    "--->\n")))::String
    if num_comments > 1
        for i = 2:num_comments
            comment_to_delete = my_comments[i]
            try
                my_retry(() -> delete_comment!(repo,
                                               pr,
                                               comment_to_delete;
                                               auth = auth))
            catch ex
                @error("Ignoring error: ", exception=(ex, catch_backtrace()))
            end
        end
        comment_to_update = my_comments[1]
        if strip(comment_to_update.body) != strip(_body)
            my_retry(() -> edit_comment!(repo,
                                         pr,
                                         comment_to_update,
                                         _body;
                                         auth = auth))
        end
    elseif num_comments == 1
        comment_to_update = my_comments[1]
        if strip(comment_to_update.body) != strip(_body)
            my_retry(() -> edit_comment!(repo,
                                         pr,
                                         comment_to_update,
                                         _body;
                                         auth = auth))
        end
    else
        always_assert(num_comments < 1)
        my_retry(() -> post_comment!(repo,
                                     pr,
                                     _body;
                                     auth = auth))
    end
    return nothing
end


text:
struct FilePicker{T}
    files::Observable{T}
end

function jsrender(session::Session, fp::FilePicker)
    text = Observable{String}("no files selected")
    input = DOM.input(
        type="file",
        onchange=js"""
            $(UtilitiesJS).readFiles(this.files, $(fp.files));
            const text = [...this.files].map(file => file.name).join(' ');
            text && JSServe.update_obs($(text), text);
        """,
        multiple=true,
        style="display:none;",
    )
    trigger = js"$(input).click()"
    btn = DOM.button(
        class="text-blue-800 text-xl font-semibold rounded mr-4 p-4 text-left w-full font-semibold hover:bg-gray-200 focus:outline-none",
        onclick=trigger,
        "Upload Files",
    )
    p = DOM.p(
        class="px-4 py-2 bg-white w-full text-gray-700 hover:text-gray-900 cursor-pointer hover:bg-gray-200 max-h-64 overflow-y-hidden",
        onclick=trigger,
        text
    )
    return jsrender(session, DOM.div(input, btn, p))
end

function to_vec(dict)
    list = valtype(dict)[]
    i = 0
    while (el = get(dict, string(i), nothing)) !== nothing
        push!(list, el)
        i += 1
    end
    return list
end

function readfiles(files::AbstractDict)
    datasets = [CSV.File(codeunits(file)) for file in to_vec(files)]
    return SimpleTable(reduce(vcat, datasets))
end


text:
# Make modification to the content clear the cache
function Base.setproperty!(o::OrgDoc, name::Symbol, value)
    if name == :contents
        o.cache = OrgCache(o)
        setfield!(o, name, value)
    elseif name == :settings
        setfield!(o, name, value)
    else # Anything other than .contents or .settings is a cache reference
        setfield!(getfield(o, :cache), name, value)
    end
end

function Base.getproperty(o::OrgDoc, name::Symbol)
    if name == :contents || name == :settings || name == :cache
        getfield(o, name)
    else
        getproperty(getfield(o, :cache), name)
    end
end

function Base.getproperty(c::OrgCache, name::Symbol)
    if isnothing(getfield(c, name))
        setfield!(c, name, gencache(c, Val(name)))
    end
    getfield(c, name)
end

function gencache(c::OrgCache, ::Val{:elements})
    Vector{Union{Heading, Section, GreaterElement, Element}}(
        OrgElementIterator(c.doc::OrgDoc) |> collect)
end

function gencache(c::OrgCache, ::Val{:components})
    Vector{Component}(OrgIterator(c.doc::OrgDoc) |> collect)
end

function gencache(c::OrgCache, ::Val{:headings})
    Vector{Heading}(filter(h -> h isa Heading, c.elements::Vector{Element}))
end

function gencache(c::OrgCache, ::Val{:keywords})
    kwdict = Dict{SubString{String}, Vector{SubString{String}}}()
    for kw::Keyword in filter(k -> k isa Keyword, c.elements::Vector{Element})
        if haskey(kwdict, kw.key)
            push!(kwdict[kw.key], kw.value)
        else
            kwdict[kw.key] = [kw.value]
        end
    end
    kwdict
end

function gencache(c::OrgCache, ::Val{:macros})
    macros = Dict{SubString{String}, Function}()
    if haskey(c.keywords::Dict{SubString{String}, Vector{SubString{String}}}, "macro")
        for (name, replacement) in split.(c.keywords["macro"], r"[ \t]+", limit=2)::Vector{Vector{SubString{String}}}
            n = getproperty.(eachmatch(r"\$\d+", replacement), :match) |> unique |> length
            if startswith(replacement, "(eval ")
                if !isnothing(Sys.which("emacs"))
                    macros[name] = function(arguments)
                        if length(arguments) == n
                            elisp = replace(replacement[7:end-1], r"\$\d+" => function(i)
                                                val = arguments[parse(Int, i[2:end])]
                                                if !isnothing(match(r"^[\d\.]+$", val))
                                                    val
                                                else
                                                    sprint(show, val)
                                                end
                                            end)
                            try
                                # TODO gate this behind some sort of unsafe switch
                                read(`emacs --batch --eval "(princ $elisp)"`, String)
                            catch _
                            end
                        end
                    end
                else
                    macros[name] = _ -> nothing
                end
            else
                macros[name] = arguments -> if length(arguments) == n
                    replace(replacement, r"\$\d+" => i -> arguments[parse(Int, i[2:end])])
                end
            end
        end
    end
    macros
end

function gencache(c::OrgCache, ::Val{:footnotes})
    footnotes = Dict{Union{SubString{String}, FootnoteReference}, Tuple{Int, Union{FootnoteReference, FootnoteDefinition}}}()
    i = 1
    for f in filter(f -> f isa FootnoteDefinition || f isa FootnoteReference, c.components) |> Vector{Union{FootnoteDefinition, FootnoteReference}}
        if !isnothing(f.definition)
            footnotes[something(f.label, f)] = (i, f)
            i += 1
        end
    end
    footnotes
end


text:
#
# Types definition module. Will be included in many Manager
# related files/modules as shared object. This module should
# contain only shared object like types and may be something
# more...
#
# @author DeadbraiN
# TODO: do we need this module at all?!
module ManagerTypes
  import Creature
  import Server
  import Client
  import Config
  import World
  import Event

  export Plugin
  export ManagerData
  export OrganismTask
  export Connections

  export orgAmount
  #
  # Abstraction for plugin type. All plugins should be inherited
  # from this type like this:
  #
  # type MyPlugin <: Plugin ... end
  #
  abstract Plugin
  #
  # One task related to one organism
  #
  type OrganismTask
    #
    # Organism unique id
    #
    id::UInt
    #
    # Task object. With it we may use green
    #
    task::Task
    #
    # One organism
    #
    organism::Creature.Organism
  end
  #
  # Current Manager connection objects. They are: server and
  # all four clients. "frozen" field is used for storing "frozen"
  # organisms (which are transfserring from current Manager to
  # another one by network). "streaming" flag means, that streaming
  # mode is on or off. Here streaming is a world dots streaming.
  #
  type Connections
    streamInit::Bool
    server    ::Server.ServerConnection
    fastServer::Server.ServerConnection
    left      ::Client.ClientConnection
    right     ::Client.ClientConnection
    up        ::Client.ClientConnection
    down      ::Client.ClientConnection
    frozen    ::Dict{UInt, Creature.Organism}
    Connections() = new()
    Connections(
      streamInit::Bool,
      server    ::Server.ServerConnection,
      fastServer::Server.ServerConnection,
      left      ::Client.ClientConnection,
      right     ::Client.ClientConnection,
      up        ::Client.ClientConnection,
      down      ::Client.ClientConnection,
      frozen    ::Dict{UInt, Creature.Organism}
    ) = new(streamInit, server, fastServer, left, right, up, down, frozen)
  end
  #
  # Manager's related type. Contains world, command line parameters,
  # organisms map and so on... If some fields will be changed, don't
  # forget to change them in recover() function.
  #
  type ManagerData
    #
    # Application wide configuration
    #
    cfg::Config.ConfigData
    #
    # Instance of the world
    #
    world::World.Plane
    #
    # Positions map, which stores positions of all organisms. Is used
    # for fast access to the organism by it's coordinates.
    #
    positions::Dict{Int, Creature.Organism}
    #
    # Map of organisms by id
    #
    organisms::Dict{UInt, Creature.Organism}
    #
    # All available organism's tasks
    #
    tasks::Array{OrganismTask, 1}
    #
    # Parameters passed through command line
    #
    params::Dict{String, String}
    #
    # Unique id of organism. It's increased every time, when new
    # organism will be created
    #
    organismId::UInt
    #
    # Total amount of organisms: alive + dead
    #
    totalOrganisms::UInt
    #
    # If true, then minimum terminal messages will be posted
    #
    quiet::Bool
    #
    # Callback, which is called when at least one dot in a
    # world has changed it's color
    #
    dotCallback::Function
    #
    # Callback, which is called when one dot in a
    # world has changed it's position (moves from one
    # position to another).
    #
    moveCallback::Function
    #
    # Manager's task (main task)
    #
    task::Task
    #
    # Energy of maximum energetic organism
    #
    maxEnergy::Int
    #
    # Amount of iterations per second
    #
    ips::Float64
    #
    # Stack of killed organism id's
    #
    killed::Array{Int, 1}
    #
    # All available(compiled) Manager's plugins. Key is a plugin string
    # name (filename). Value - plugin related data type.
    #
    plugins::Dict{String, Plugin}
    #
    # Manager observer
    #
    obs::Event.Observer
    #
    # Manager connections (with other managers, terminals, visualizer etc...)
    #
    cons::Connections
    #
    # Short constructor
    #
    ManagerData(
      cfg::Config.ConfigData,
      world::World.Plane,
      positions::Dict{Int, Creature.Organism},
      organisms::Dict{UInt, Creature.Organism},
      tasks::Array{OrganismTask, 1},
      params::Dict{String, String},
      organismId::UInt,
      totalOrganisms::UInt,
      quiet::Bool,
      dotCallback::Function,
      moveCallback::Function,
      task::Task,
      maxEnergy::Int,
      ips::Float64,
      killed::Array{Int, 1},
      plugins::Dict{String, Plugin},
      obs::Event.Observer
    ) = new(
      cfg,
      world,
      positions,
      organisms,
      tasks,
      params,
      organismId,
      totalOrganisms,
      quiet,
      dotCallback,
      moveCallback,
      task,
      maxEnergy,
      ips,
      killed,
      plugins,
      obs
    )
  end

  #
  # Returns amount of alive organisms at the moment. This function is Here
  # and not in Manager module, because it's used in different other modules.
  # In case if you need some function, which work with manager's data and is
  # called only within Manager module and it's sub modules you have to add
  # it there.
  # @return {Int}
  #
  function orgAmount(man::ManagerTypes.ManagerData)
    length(man.tasks) - length(man.killed)
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo