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

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

# ╔═╡ 2b37ca3a-0970-11eb-3c3d-4f788b411d1a
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ 2dcb18d0-0970-11eb-048a-c1734c6db842
begin
	Pkg.add(["PlutoUI", "Plots"])

	using Plots
	gr()
	using PlutoUI
end

# ╔═╡ 19fe1ee8-0970-11eb-2a0d-7d25e7d773c6
md"_homework 5, version 0_"

# ╔═╡ 49567f8e-09a2-11eb-34c1-bb5c0b642fe8
# WARNING FOR OLD PLUTO VERSIONS, DONT DELETE ME

html"""

"""

# ╔═╡ 181e156c-0970-11eb-0b77-49b143cc0fc0
md"""

# **Homework 5**: _Epidemic modeling II_
`18.S191`, fall 2020

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!
"""

# ╔═╡ 1f299cc6-0970-11eb-195b-3f951f92ceeb
# edit the code below to set your name and kerberos ID (i.e. email without @mit.edu)

student = (name = "Daniel Fryer", 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

# ╔═╡ 1bba5552-0970-11eb-1b9a-87eeee0ecc36
md"""

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

# ╔═╡ 2848996c-0970-11eb-19eb-c719d797c322
md"_Let's create a package environment:_"

# ╔═╡ 69d12414-0952-11eb-213d-2f9e13e4b418
md"""
In this problem set, we will look at a simple **spatial** agent-based epidemic model: agents can interact only with other agents that are *nearby*.  (In the previous homework any agent could interact with any other, which is not realistic.)

A simple approach is to use **discrete space**: each agent lives
in one cell of a square grid. For simplicity we will allow no more than
one agent in each cell, but this requires some care to
design the rules of the model to respect this.

We will adapt some functionality from the previous homework. You should copy and paste your code from that homework into this notebook.
"""

# ╔═╡ 3e54848a-0954-11eb-3948-f9d7f07f5e23
md"""
## **Exercise 1:** _Wandering at random in 2D_

In this exercise we will implement a **random walk** on a 2D lattice (grid). At each time step, a walker jumps to a neighbouring position at random (i.e. chosen with uniform probability from the available adjacent positions).

"""

# ╔═╡ 3e623454-0954-11eb-03f9-79c873d069a0
md"""
#### Exercise 1.1
We define a struct type `Coordinate` that contains integers `x` and `y`.
"""

# ╔═╡ 0ebd35c8-0972-11eb-2e67-698fd2d311d2
begin
struct Coordinate
	x::Int64
	y::Int64
end
Coordinate() = Coordinate(0,0)
get_tuple(c::Coordinate) = (c.x,c.y)
end

# ╔═╡ 027a5f48-0a44-11eb-1fbf-a94d02d0b8e3
md"""
👉 Construct a `Coordinate` located at the origin.
"""

# ╔═╡ b2f90634-0a68-11eb-1618-0b42f956b5a7
origin = Coordinate()

# ╔═╡ 3e858990-0954-11eb-3d10-d10175d8ca1c
md"""
👉 Write a function `make_tuple` that takes an object of type `Coordinate` and returns the corresponding tuple `(x, y)`. Boring, but useful later!
"""

# ╔═╡ 189bafac-0972-11eb-1893-094691b2073c
function make_tuple(c)
 	(c.x,c.y)
end

# ╔═╡ 73ed1384-0a29-11eb-06bd-d3c441b8a5fc
md"""
#### Exercise 1.2
In Julia, operations like `+` and `*` are just functions, and they are treated like any other function in the language. The only special property you can use the _infix notation_: you can write
```julia
1 + 2
```
instead of
```julia
+(1, 2)
```
_(There are [lots of special 'infixable' function names](https://github.com/JuliaLang/julia/blob/master/src/julia-parser.scm#L23-L24) that you can use for your own functions!)_

When you call it with the prefix notation, it becomes clear that it really is 'just another function', with lots of predefined methods.
"""

# ╔═╡ 96707ef0-0a29-11eb-1a3e-6bcdfb7897eb
+(1, 2)

# ╔═╡ b0337d24-0a29-11eb-1fab-876a87c0973f
+

# ╔═╡ 9c9f53b2-09ea-11eb-0cda-639764250cee
md"""
> #### Extending + in the wild
> Because it is a function, we can add our own methods to it! This feature is super useful in general languages like Julia and Python, because it lets you use familiar syntax (`a + b*c`) on objects that are not necessarily numbers!
>
> One example we've see before is the `RGB` type in Homework 1. You are able to do:
> ```julia
> 0.5 * RGB(0.1, 0.7, 0.6)
> ```
> to multiply each color channel by $0.5$. This is possible because `Images.jl` [wrote a method](https://github.com/JuliaGraphics/ColorVectorSpace.jl/blob/master/src/ColorVectorSpace.jl#L131):
> ```julia
> *(::Real, ::AbstractRGB)::AbstractRGB
> ```

👉 Implement addition on two `Coordinate` structs by adding a method to `Base.:+`
"""

# ╔═╡ e24d5796-0a68-11eb-23bb-d55d206f3c40
function Base.:+(a::Coordinate, b::Coordinate)::Coordinate
 	return Coordinate(a.x + b.x, a.y + b.y)
end

# ╔═╡ ec8e4daa-0a2c-11eb-20e1-c5957e1feba3
Coordinate(3,4) + Coordinate(10,10) # uncomment to check + works

# ╔═╡ e144e9d0-0a2d-11eb-016e-0b79eba4b2bb
md"""
_Pluto has some trouble here, you need to manually re-run the cell above!_
"""

# ╔═╡ 71c358d8-0a2f-11eb-29e1-57ff1915e84a
md"""
#### Exercise 1.3
In our model, agents will be able to walk in 4 directions: up, down, left and right. We can define these directions as `Coordinate`s.
"""

# ╔═╡ 5278e232-0972-11eb-19ff-a1a195127297
# uncomment this:

possible_moves = [
 	Coordinate( 1, 0),
 	Coordinate( 0, 1),
 	Coordinate(-1, 0),
 	Coordinate( 0,-1),
]

# ╔═╡ 71c9788c-0aeb-11eb-28d2-8dcc3f6abacd
md"""
👉 `rand(possible_moves)` gives a random possible move. Add this to the coordinate `Coordinate(4,5)` and see that it moves to a valid neighbor.
"""

# ╔═╡ 69151ce6-0aeb-11eb-3a53-290ba46add96
rand(possible_moves) + Coordinate(4,5)

# ╔═╡ 3eb46664-0954-11eb-31d8-d9c0b74cf62b
md"""
We are able to make a `Coordinate` perform one random step, by adding a move to it. Great!

👉 Write a function `trajectory` that calculates a trajectory of a `Wanderer` `w` when performing `n` steps., i.e. the sequence of positions that the walker finds itself in.

Possible steps:
- Use `rand(possible_moves, n)` to generate a vector of `n` random moves. Each possible move will be equally likely.
- To compute the trajectory you can use either of the following two approaches:
  1. 🆒 Use the function `accumulate` (see the live docs for `accumulate`). Use `+` as the function passed to `accumulate` and the `w` as the starting value (`init` keyword argument).
  1. Use a `for` loop calling `+`.

"""

# ╔═╡ edf86a0e-0a68-11eb-2ad3-dbf020037019
function trajectory(w::Coordinate, n::Int)
	return accumulate(+, vcat(w, rand(possible_moves, n)))
end

# ╔═╡ 00f01a64-2993-11eb-2431-47731cdf2d92
function lame_trajectory(w::Coordinate, n::Int)
	moves = rand(possible_moves, n)
	traj = copy(moves)
	for i in 1:length(moves)
		@views traj[i] = sum(moves[1:i])
	end
	return traj
end

# ╔═╡ 478309f4-0a31-11eb-08ea-ade1755f53e0
function plot_trajectory!(p::Plots.Plot, trajectory::Vector; kwargs...)
	plot!(p, make_tuple.(trajectory);
		label=nothing,
		linewidth=2,
		linealpha=LinRange(1.0, 0.2, length(trajectory)),
		kwargs...)
end

# ╔═╡ 3ebd436c-0954-11eb-170d-1d468e2c7a37
md"""
#### Exercise 1.4
👉 Plot 10 trajectories of length 1000 on a single figure, all starting at the origin. Use the function `plot_trajectory!` as demonstrated above.

Remember from last week that you can compose plots like this:

```julia
let
	# Create a new plot with aspect ratio 1:1
	p = plot(ratio=1)

	plot_trajectory!(p, test_trajectory)      # plot one trajectory
	plot_trajectory!(p, another_trajectory)   # plot the second one
	...

	p
end
```
"""

# ╔═╡ b4d5da4a-09a0-11eb-1949-a5807c11c76c
md"""
#### Exercise 1.5
Agents live in a box of side length $2L$, centered at the origin. We need to decide (i.e. model) what happens when they reach the walls of the box (boundaries), in other words what kind of **boundary conditions** to use.

One relatively simple boundary condition is a **collision boundary**:

> Each wall of the box is a wall, modelled using "collision": if the walker tries to jump beyond the wall, it ends up at the position inside the box that is closest to the goal.

👉 Write a function `collide_boundary` which takes a `Coordinate` `c` and a size $L$, and returns a new coordinate that lies inside the box (i.e. ``[-L,L]\times [-L,L]``), but is closest to `c`. This is similar to `extend_mat` from Homework 1.
"""

# ╔═╡ 0237ebac-0a69-11eb-2272-35ea4e845d84
begin
function collide_boundary(x::Int64, L::Number)
	if (abs(x) > L) x = sign(x)*L end
	return x
end
function collide_boundary(c::Coordinate, L::Number)
	return Coordinate(collide_boundary(c.x, L), collide_boundary(c.y, L))
end
end

# ╔═╡ ad832360-0a40-11eb-2857-e7f0350f3b12
collide_boundary(Coordinate(12,4), 10) # uncomment to test

# ╔═╡ b4ed2362-09a0-11eb-0be9-99c91623b28f
md"""
#### Exercise 1.6
👉  Implement a 3-argument method  of `trajectory` where the third argument is a size. The trajectory returned should be within the boundary (use `collide_boundary` from above). You can still use `accumulate` with an anonymous function that makes a move and then reflects the resulting coordinate, or use a for loop.

"""

# ╔═╡ 0665aa3e-0a69-11eb-2b5d-cd718e3c7432
function trajectory(c::Coordinate, n::Int, L::Number)
	return accumulate((c₁,c₂)->collide_boundary(c₁+c₂,L),
		vcat(c,rand(possible_moves, n)))
end

# ╔═╡ 44107808-096c-11eb-013f-7b79a90aaac8
test_trajectory = trajectory(Coordinate(4,4), 30) # uncomment to test

# ╔═╡ 87ea0868-0a35-11eb-0ea8-63e27d8eda6e
try
	p = plot(ratio=1, size=(650,200))
	plot_trajectory!(p, test_trajectory; color="black", showaxis=false, axis=nothing, linewidth=4)
	p
catch
end

# ╔═╡ 51788e8e-0a31-11eb-027e-fd9b0dc716b5
 	let
 		long_trajectory = trajectory(Coordinate(4,4), 1000)

 		p = plot(ratio=1)
 		plot_trajectory!(p, long_trajectory)
 		p
 	end

# ^ uncomment to visualize a trajectory

# ╔═╡ dcefc6fe-0a3f-11eb-2a96-ddf9c0891873
let
	p = plot(ratio=1)
	for _ in 1:10
		long_trajectory = trajectory(Coordinate(), 1000)
		plot_trajectory!(p, long_trajectory)
	end
	p
end

# ╔═╡ 475f01b2-2999-11eb-3993-b5666bd60bfe
let
	p = plot(ratio=1)
	for _ in 1:10
		long_trajectory = trajectory(Coordinate(), 1000, 25)
		plot_trajectory!(p, long_trajectory)
	end
	p
end

# ╔═╡ 3ed06c80-0954-11eb-3aee-69e4ccdc4f9d
md"""
## **Exercise 2:** _Wanderering Agents_

In this exercise we will create Agents which have a location as well as some infection state information.

Let's define a type `Agent`. `Agent` contains a `position` (of type `Coordinate`), as well as a `status` of type `InfectionStatus` (as in Homework 4).)

(For simplicity we will not use a `num_infected` field, but feel free to do so!)
"""

# ╔═╡ 35537320-0a47-11eb-12b3-931310f18dec
@enum InfectionStatus S I R

# ╔═╡ cf2f3b98-09a0-11eb-032a-49cc8c15e89c
begin
abstract type AbstractAgent end
mutable struct Agent <: AbstractAgent
	position::Coordinate
	status::InfectionStatus
end
end

# ╔═╡ 814e888a-0954-11eb-02e5-0964c7410d30
md"""
#### Exercise 2.1
👉 Write a function `initialize` that takes parameters $N$ and $L$, where $N$ is the number of agents abd $2L$ is the side length of the square box where the agents live.

It returns a `Vector` of `N` randomly generated `Agent`s. Their coordinates are randomly sampled in the ``[-L,L] \times [-L,L]`` box, and the agents are all susceptible, except one, chosen at random, which is infectious.
"""

# ╔═╡ 0cfae7ba-0a69-11eb-3690-d973d70e47f4
function initialize(N::Number, L::Number)
	L = floor(Int, L)
	xcoords = rand(-L:L, N)
	ycoords = rand(-L:L, N)
	coords = Coordinate.(xcoords, ycoords)
	agents = Agent.(coords, S)
	random_agent = agents[rand(1:length(agents))]
	random_agent.status = I
	return agents
end

# ╔═╡ 1d0f8eb4-0a46-11eb-38e7-63ecbadbfa20
initialize(3, 10)

# ╔═╡ e0b0880c-0a47-11eb-0db2-f760bbbf9c11
# Color based on infection status
color(s::InfectionStatus) = if s == S
	"blue"
elseif s == I
	"red"
else
	"green"
end

# ╔═╡ b5a88504-0a47-11eb-0eda-f125d419e909
position(a::Agent) = a.position # uncomment this line

# ╔═╡ 87a4cdaa-0a5a-11eb-2a5e-cfaf30e942ca
color(a::Agent) = color(a.status) # uncomment this line

# ╔═╡ 49fa8092-0a43-11eb-0ba9-65785ac6a42f
md"""
#### Exercise 2.2
👉 Write a function `visualize` that takes in a collection of agents as argument, and the box size `L`. It should plot a point for each agent at its location, coloured according to its status.

You can use the keyword argument `c=color.(agents)` inside your call to the plotting function make the point colors correspond to the infection statuses. Don't forget to use `ratio=1`.
"""

# ╔═╡ f953e06e-099f-11eb-3549-73f59fed8132
md"""

### Exercise 3: Spatial epidemic model -- Dynamics

Last week we wrote a function `interact!` that takes two agents, `agent` and `source`, and an infection of type `InfectionRecovery`, which models the interaction between two agent, and possibly modifies `agent` with a new status.

This week, we define a new infection type, `CollisionInfectionRecovery`, and a new method that is the same as last week, except it **only infects `agent` if `agents.position==source.position`**.
"""

# ╔═╡ e6dd8258-0a4b-11eb-24cb-fd5b3554381b
abstract type AbstractInfection end

# ╔═╡ de88b530-0a4b-11eb-05f7-85171594a8e8
struct CollisionInfectionRecovery <: AbstractInfection
	p_infection::Float64
	p_recovery::Float64
end

# ╔═╡ 80f39140-0aef-11eb-21f7-b788c5eab5c9
md"""

Write a function `interact!` that takes two `Agent`s and a `CollisionInfectionRecovery`, and:

- If the agents are at the same spot, causes a susceptible agent to communicate the desease from an infectious one with the correct probability.
- if the first agent is infectious, it recovers with some probability
"""

# ╔═╡ ee5c548e-3d11-11eb-0367-6555ca14c8a7
begin
	function set_status!(agent, new_status)
		agent.status = new_status
	end
end

# ╔═╡ d1bcd5c4-0a4b-11eb-1218-7531e367a7ff
function interact!(agent::Agent, source::Agent, infection::CollisionInfectionRecovery)
	if agent.status == I
		if rand() < infection.p_recovery
			set_status!(agent, R)
	    end
	elseif agent.status == S && source.status == I
		if agent.position == source.position
			set_status!(agent, I)
		end
	end
end

# ╔═╡ 34778744-0a5f-11eb-22b6-abe8b8fc34fd
md"""
#### Exercise 3.1
Your turn!

👉 Write a function `step!` that takes a vector of `Agent`s, a box size `L` and an `infection`. This that does one step of the dynamics on a vector of agents.

- Choose an Agent `source` at random.

- Move the `source` one step, and use `collide_boundary` to ensure that our agent stays within the box.

- For all _other_ agents, call `interact!(other_agent, source, infection)`.

- return the array `agents` again.
"""

# ╔═╡ 2d5f2b1e-3d14-11eb-2d35-7d35f97e9412
function agent_step!(source::AbstractAgent, L::Number)
	source.position = collide_boundary(source.position + rand(possible_moves),   L)
end

# ╔═╡ 1fc3271e-0a45-11eb-0e8d-0fd355f5846b
md"""
#### Exercise 3.2
If we call `step!` `N` times, then every agent will have made one step, on average. Let's call this one _sweep_ of the simulation.

👉 Create a before-and-after plot of ``k_{sweeps}=1000`` sweeps.

- Initialize a new vector of agents (`N=50`, `L=40`, `infection` is given as `pandemic` below).
- Plot the state using `visualize`, and save the plot as a variable `plot_before`.
- Run `k_sweeps` sweeps.
- Plot the state again, and store as `plot_after`.
- Combine the two plots into a single figure using
```julia
plot(plot_before, plot_after)
```
"""

# ╔═╡ 18552c36-0a4d-11eb-19a0-d7d26897af36
pandemic = CollisionInfectionRecovery(0.5, 0.00001)

# ╔═╡ 4e7fd58a-0a62-11eb-1596-c717e0845bd5
@bind k_sweeps Slider(1:10000, default=1000)

# ╔═╡ e964c7f0-0a61-11eb-1782-0b728fab1db0
md"""
#### Exercise 3.3

Every time that you move the slider, a completely new simulation is created an run. This makes it hard to view the progress of a single simulation over time. So in this exercise, we we look at a single simulation, and plot the S, I and R curves.

👉 Plot the SIR curves of a single simulation, with the same parameters as in the previous exercise. Use `k_sweep_max = 10000` as the total number of sweeps.
"""

# ╔═╡ 4d83dbd0-0a63-11eb-0bdc-757f0e721221
k_sweep_max = 10_000

# ╔═╡ 62230726-3d19-11eb-17d3-e7aab3048862
get_status(a::AbstractAgent) = a.status

# ╔═╡ 201a3810-0a45-11eb-0ac9-a90419d0b723
md"""
#### Exercise 3.4 (optional)
Let's make our plot come alive! There are two options to make our visualization dynamic:

👉1⃣ Precompute one simulation run and save its intermediate states using `deepcopy`. You can then write an interactive visualization that shows both the state at time $t$ (using `visualize`) and the history of $S$, $I$ and $R$ from time $0$ up to time $t$. $t$ is controlled by a slider.

👉2⃣ Use `@gif` from Plots.jl to turn a sequence of plots into an animation. Be careful to skip about 50 sweeps between each animation frame, otherwise the GIF becomes too large.

This an optional exercise, and our solution to 2⃣ is given below.
"""

# ╔═╡ cd931cf2-3d33-11eb-0565-abfc49c2e4c0

# ╔═╡ a3f180f4-3d34-11eb-3c3e-0b191891eeaa
@bind t Slider(1:100, default=10)

# ╔═╡ 2031246c-0a45-11eb-18d3-573f336044bf
md"""
#### Exercise 3.5
👉  Using $L=20$ and $N=100$, experiment with the infection and recovery probabilities until you find an epidemic outbreak. (Take the recovery probability quite small.) Modify the two infections below to match your observations.
"""

# ╔═╡ 63dd9478-0a45-11eb-2340-6d3d00f9bb5f
causes_outbreak = CollisionInfectionRecovery(0.5, 0.0001)

# ╔═╡ 269955e4-0a46-11eb-02cc-1946dc918bfa
does_not_cause_outbreak = CollisionInfectionRecovery(0.5, 0.00095)

# ╔═╡ 4d4548fe-0a66-11eb-375a-9313dc6c423d

# ╔═╡ 20477a78-0a45-11eb-39d7-93918212a8bc
md"""
#### Exercise 3.6
👉 With the parameters of Exercise 3.2, run 50 simulations. Plot $S$, $I$ and $R$ as a function of time for each of them (with transparency!). This should look qualitatively similar to what you saw in the previous homework. You probably need different `p_infection` and `p_recovery` values from last week. Why?
"""

# ╔═╡ b9bc85c6-3e0d-11eb-3b70-2bf813cbbb01
[1 2; 3 4]'

# ╔═╡ b1b1afda-0a66-11eb-2988-752405815f95
need_different_parameters_because = md"""
i say so
"""

# ╔═╡ 05c80a0c-09a0-11eb-04dc-f97e306f1603
md"""
## **Exercise 4:** _Effect of socialization_

In this exercise we'll modify the simple mixing model. Instead of a constant mixing probability, i.e. a constant probability that any pair of people interact on a given day, we will have a variable probability associated with each agent, modelling the fact that some people are more or less social or contagious than others.
"""

# ╔═╡ b53d5608-0a41-11eb-2325-016636a22f71
md"""
#### Exercise 4.1
We create a new agent type `SocialAgent` with fields `position`, `status`, `num_infected`, and `social_score`. The attribute `social_score` represents an agent's probability of interacting with any other agent in the population.
"""

# ╔═╡ 1b5e72c6-0a42-11eb-3884-a377c72270c7
mutable struct SocialAgent <: AbstractAgent
	position::Coordinate
	status::InfectionStatus
	num_infected::Int64
	social_score::Float64
end

# ╔═╡ c704ea4c-0aec-11eb-2f2c-859c954aa520
md"""define the `position` and `color` methods for `SocialAgent` as we did for `Agent`. This will allow the `visualize` function to work. on both kinds of Agents"""

# ╔═╡ e97e39aa-0a5d-11eb-3d5f-f90a0acfe5a2
begin
	position(a::SocialAgent) = a.position
	color(a::SocialAgent) = color(a.status)
end

# ╔═╡ 1ccc961e-0a69-11eb-392b-915be07ef38d
function visualize(agents::Vector, L)
	positions = get_tuple.(position.(agents))
	p = plot(positions, seriestype = :scatter, ratio = 1, c = color.(agents),
		legend = false)
	return p
end

# ╔═╡ 1f96c80a-0a46-11eb-0690-f51c60e57c3f
let
	N = 20
	L = 10
	visualize(initialize(N, L), L) # uncomment this line!
end

# ╔═╡ b554b654-0a41-11eb-0e0d-e57ff68ced33
md"""
👉 Create a function `initialize_social` that takes `N` and `L`, and creates N agents  within a 2L x 2L box, with `social_score`s chosen from 10 equally-spaced between 0.1 and 0.5. (see LinRange)
"""

# ╔═╡ 40c1c1d6-0a69-11eb-3913-59e9b9ec4332
function initialize_social(N, L)
	L = floor(Int, L)
	xcoords = rand(-L:L, N)
	ycoords = rand(-L:L, N)
	coords = Coordinate.(xcoords, ycoords)
	scores = rand(LinRange(0.1, 0.5, 10), N)
	n_infect = zeros(Int, N)
	agents = SocialAgent.(coords, S, n_infect, scores)
    agents[rand(1:N)].status = I
	return agents
end

# ╔═╡ a47b5dd8-3eba-11eb-1ce3-f574ecf41690
initialize_social(5, 2)

# ╔═╡ 18ac9926-0aed-11eb-034f-e9849b71c9ac
md"""
Now that we have 2 agent types

1. let's create an AbstractAgent type
2. Go back in the notebook and make the agent types a subtype of AbstractAgent.

"""

# ╔═╡ b56ba420-0a41-11eb-266c-719d39580fa9
md"""
#### Exercise 4.2
Not all two agents who end up in the same grid point may actually interact in an infectious way -- they may just be passing by and do not create enough exposure for communicating the disease.

👉 Write a new `interact!` method on `SocialAgent` which adds together the social_scores for two agents and uses that as the probability that they interact in a risky way. Only if they interact in a risky way, the infection is transmitted with the usual probability.
"""

# ╔═╡ 465e918a-0a69-11eb-1b59-01150b4b0f36
function interact!(agent::SocialAgent, source::SocialAgent, infection::CollisionInfectionRecovery)
	if agent.status == I
		if rand() < infection.p_recovery
			set_status!(agent, R)
	    end
	elseif agent.status == S && source.status == I
		if agent.position == source.position
			score = agent.social_score + source.social_score
			if rand() < score
				set_status!(agent, I)
				source.num_infected += 1
			end
		end
	end
end

# ╔═╡ 24fe0f1a-0a69-11eb-29fe-5fb6cbf281b8
function step!(agents::Vector, L::Number, infection::AbstractInfection)
	indices = 1:length(agents)
	source_index = rand(indices)
	agent_step!(agents[source_index], L)
	for i in indices
		if i != source_index
			interact!(agents[i], agents[source_index], infection)
		end
	end
return agents
end

# ╔═╡ 778c2490-0a62-11eb-2a6c-e7fab01c6822
let
	N = 50
	L = 40
	agents = initialize(N, L)

	plot_before = visualize(agents, L)

	# sweeps
    for _ in 1:k_sweeps*N
		step!(agents, L, pandemic)
	end

	plot_after = visualize(agents, L)

	plot(plot_before, plot_after)
end

# ╔═╡ ef27de84-0a63-11eb-177f-2197439374c5
let
	N = 50
	L = 30

	S_n, I_n, R_n = [], [], []

	agents = initialize(N, L)
	# compute k_sweep_max number of sweeps and plot the SIR

	for _ in 1:k_sweep_max
		for _ in 1:N
			step!(agents, L, pandemic)
		end
		push!(S_n, sum(get_status.(agents) .== S) )
		push!(I_n, sum(get_status.(agents) .== I) )
		push!(R_n, sum(get_status.(agents) .== R) )
	end

	p = plot(S_n, label = "susceptible")
	plot!(p, I_n, label = "infectious")
	plot!(p, R_n, label = "recovered")
end

# ╔═╡ e5040c9e-0a65-11eb-0f45-270ab8161871
counts, states = let
	N = 50
	L = 30

	states = []
	counts = []
	S_n, I_n, R_n = [], [], []
	agents = initialize(N, L)

	for sweep in 1:k_sweep_max
		for _ in 1:N
			step!(agents, L, pandemic)
		end
		push!(S_n, sum(get_status.(agents) .== S) )
		push!(I_n, sum(get_status.(agents) .== I) )
		push!(R_n, sum(get_status.(agents) .== R) )
		if sweep % 50 == 0
			push!(counts, deepcopy((S_n, I_n, R_n)))
			push!(states, deepcopy(agents))
		end
	end
	counts, states
end

# ╔═╡ 90794bba-3d34-11eb-2f70-0da9e08ea54d
begin
p = plot(counts[t][1], color = color(S), label = "susceptible")
plot!(p, counts[t][2], color = color(I), label = "infected")
plot!(p, counts[t][3], color = color(R), label = "recovered")
plot(p, visualize(states[t], 30))
end

# ╔═╡ f55659d2-3d35-11eb-3675-bfc8d6d0d687
let

	N = 100
	L = 20
	pandemic = CollisionInfectionRecovery(0.5, 0.000099)

	agents = initialize(N, L)
	S_n, I_n, R_n = [],[],[]

	Tmax = 100
	@gif for t in 1:Tmax
		for _ in 1:50N
			step!(agents, L, pandemic)
		end

		push!(S_n, sum(get_status.(agents) .== S))
		push!(I_n, sum(get_status.(agents) .== I))
		push!(R_n, sum(get_status.(agents) .== R))

		p = plot(S_n, color = color(S), label = "susceptible",
		         xlim = (1,Tmax), ylim = (1,N))
        plot!(p, I_n, color = color(I), label = "infected")
        plot!(p, R_n, color = color(R), label = "recovered")

		plot(p, visualize(agents, L))
	end
end

# ╔═╡ 601f4f54-0a45-11eb-3d6c-6b9ec75c6d4a
S_m, I_m, R_m = let
	sims = 50
	N = 100
	L = 20
	k_sweep_max= 5000
	pandemic = CollisionInfectionRecovery(0.5, 0.00001)

	S_m = zeros(sims,k_sweep_max)
	I_m = zeros(sims,k_sweep_max)
	R_m = zeros(sims,k_sweep_max)
	for i in 1:sims
		agents = initialize(N, L)
		for k in 1:k_sweep_max
			for j in 1:N
				step!(agents, L, pandemic)
			end
			S_m[i,k] = sum(get_status.(agents) .== S)
			I_m[i,k] = sum(get_status.(agents) .== I)
			R_m[i,k] = sum(get_status.(agents) .== R)
		end
	end

	S_m, I_m, R_m
end

# ╔═╡ 83c7d300-3eb7-11eb-22d3-0dbc0dc0bfa4
begin
	p2 = plot(S_m', legend = false, colour = "blue", alpha = 0.2)
	plot!(p2, I_m', legend = false, colour = "red", alpha = 0.2)
	plot!(p2, R_m', legend = false, colour = "green", alpha = 0.2)
end

# ╔═╡ a0f5cbf6-3ebf-11eb-3200-8f8286079f4d
begin
  	test_agents = initialize_social(5, 2)
 	position(test_agents[1])
	color(test_agents[1])
	step!(test_agents, 2, pandemic)
end

# ╔═╡ a885bf78-0a5c-11eb-2383-9d74c8765847
md"""
Make sure `step!`, `position`, `color`, work on the type `SocialAgent`. If `step!` takes an untyped first argument, it should work for both Agent and SocialAgent types without any changes. We actually only need to specialize `interact!` on SocialAgent.

#### Exercise 4.3
👉 Plot the SIR curves of the resulting simulation.

N = 50;
L = 40;
number of steps = 200

In each step call `step!` 50N times.
"""

# ╔═╡ 1f172700-0a42-11eb-353b-87c0039788bd
let
	N = 50
	L = 40

	global social_agents = initialize_social(N, L)
	Ss, Is, Rs = [], [], []

	pandemic = CollisionInfectionRecovery(0.5, 0.00001)
	Tmax = 200

	@gif for t in 1:Tmax
		for _ in 1:50N
			step!(social_agents, L, pandemic)
		end

		push!(Ss, count(==(S), get_status.(social_agents)))
		push!(Is, count(==(I), get_status.(social_agents)))
		push!(Rs, count(==(R), get_status.(social_agents)))

		v = visualize(social_agents, L)
		p = plot(Ss, color = color(S), label = "susceptible",
		         xlim = (1,Tmax), ylim = (1,N))
        plot!(p, Is, color = color(I), label = "infected")
        plot!(p, Rs, color = color(R), label = "recovered")
		plot(p,v, size = (600,300))
	end
end

# ╔═╡ b59de26c-0a41-11eb-2c67-b5f3c7780c91
md"""
#### Exercise 4.4
👉 Make a scatter plot showing each agent's `social_score` on one axis, and the `num_infected` from the simulation in the other axis. Run this simulation several times and comment on the results.
"""

# ╔═╡ faec52a8-0a60-11eb-082a-f5787b09d88c
begin
	get_score(a::SocialAgent) = a.social_score
	get_num_infected(a::SocialAgent) = a.num_infected
	plot(get_score.(social_agents), get_num_infected.(social_agents),
	seriestype = :scatter, legend = false, xlim = (0.09,0.51))
end

# ╔═╡ b5b4d834-0a41-11eb-1b18-1bd626d18934
md"""
👉 Run a simulation for 100 steps, and then apply a "lockdown" where every agent's social score gets multiplied by 0.25, and then run a second simulation which runs on that same population from there.  What do you notice?  How does changing this factor form 0.25 to other numbers affect things?
"""

# ╔═╡ 64782262-3ec3-11eb-3b62-d5481d23438b
begin
function lockdown_agent!(a::AbstractAgent)
	a.social_score *= 0.1
end
function open_economy!(a::AbstractAgent)
	a.social_score /= 0.1
end
end

# ╔═╡ a83c96e2-0a5a-11eb-0e58-15b5dda7d2d2
let
	N = 100
	L = 40

	global social_agents2 = initialize_social(N, L)
	Ss, Is, Rs = [], [], []

	pandemic = CollisionInfectionRecovery(0.7, 0.00001)
	Tmax = 200

	@gif for t in 1:Tmax
		if (t == 25) lockdown_agent!.(social_agents2) end
		if (t == 50) open_economy!.(social_agents2) end
		for _ in 1:50N
			step!(social_agents2, L, pandemic)
		end

		push!(Ss, count(==(S), get_status.(social_agents2)))
		push!(Is, count(==(I), get_status.(social_agents2)))
		push!(Rs, count(==(R), get_status.(social_agents2)))

		v = visualize(social_agents2, L)
		p = plot(Ss, color = color(S), label = "susceptible",
		         xlim = (1,Tmax), ylim = (1,N))
        plot!(p, Is, color = color(I), label = "infected")
        plot!(p, Rs, color = color(R), label = "recovered")
		plot(p,v, size = (600,300))
	end
end

# ╔═╡ 05fc5634-09a0-11eb-038e-53d63c3edaf2
md"""
## **Exercise 5:** (Optional) _Effect of distancing_

We can use a variant of the above model to investigate the effect of the
mis-named "social distancing"
(we want people to be *socially* close, but *physically* distant).

In this variant, we separate out the two effects "infection" and
"movement": an infected agent chooses a
neighbouring site, and if it finds a susceptible there then it infects it
with probability $p_I$. For simplicity we can ignore recovery.

Separately, an agent chooses a neighbouring site to move to,
and moves there with probability $p_M$ if the site is vacant. (Otherwise it
stays where it is.)

When $p_M = 0$, the agents cannot move, and hence are
completely quarantined in their original locations.

👉 How does the disease spread in this case?

"""

# ╔═╡ 24c2fb0c-0a42-11eb-1a1a-f1246f3420ff

# ╔═╡ c7649966-0a41-11eb-3a3a-57363cea7b06
md"""
👉 Run the dynamics repeatedly, and plot the sites which become infected.
"""

# ╔═╡ 2635b574-0a42-11eb-1daa-971b2596ce44

# ╔═╡ c77b085e-0a41-11eb-2fcb-534238cd3c49
md"""
👉 How does this change as you increase the *density*
    $\rho = N / (L^2)$ of agents?  Start with a small density.

This is basically the [**site percolation**](https://en.wikipedia.org/wiki/Percolation_theory) model.

When we increase $p_M$, we allow some local motion via random walks.
"""

# ╔═╡ 274fe006-0a42-11eb-1869-29193bb84957

# ╔═╡ c792374a-0a41-11eb-1e5b-89d9de2cf1f9
md"""
👉 Investigate how this leaky quarantine affects the infection dynamics with
different densities.

"""

# ╔═╡ d147f7f0-0a66-11eb-2877-2bc6680e396d

# ╔═╡ 0e6b60f6-0970-11eb-0485-636624a0f9d7
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

# ╔═╡ 0a82a274-0970-11eb-20a2-1f590be0e576
md"## Function library

Just some helper functions used in the notebook."

# ╔═╡ 0aa666dc-0970-11eb-2568-99a6340c5ebd
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]))

# ╔═╡ 8475baf0-0a63-11eb-1207-23f789d00802
hint(md"""
After every sweep, count the values $S$, $I$ and $R$ and push! them to 3 arrays.
""")

# ╔═╡ f9b9e242-0a53-11eb-0c6a-4d9985ef1687
hint(md"""
```julia
let
	N = 50
	L = 40

	x = initialize(N, L)

	# initialize to empty arrays
	Ss, Is, Rs = Int[], Int[], Int[]

	Tmax = 200

	@gif for t in 1:Tmax
		for i in 1:50N
			step!(x, L, pandemic)
		end

		#... track S, I, R in Ss Is and Rs

		left = visualize(x, L)

		right = plot(xlim=(1,Tmax), ylim=(1,N), size=(600,300))
		plot!(right, 1:t, Ss, color=color(S), label="S")
		plot!(right, 1:t, Is, color=color(I), label="I")
		plot!(right, 1:t, Rs, color=color(R), label="R")

		plot(left, right)
	end
end
```
""")

# ╔═╡ 0acaf3b2-0970-11eb-1d98-bf9a718deaee
almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))

# ╔═╡ 0afab53c-0970-11eb-3e43-834513e4632e
still_missing(text=md"Replace `missing` with your answer.") = Markdown.MD(Markdown.Admonition("warning", "Here we go!", [text]))

# ╔═╡ 0b21c93a-0970-11eb-33b0-550a39ba0843
keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]))

# ╔═╡ 0b470eb6-0970-11eb-182f-7dfb4662f827
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."]

# ╔═╡ 0b6b27ec-0970-11eb-20c2-89515ee3ab88
correct(text=rand(yays)) = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]))

# ╔═╡ ec576da8-0a2c-11eb-1f7b-43dec5f6e4e7
let
	# we need to call Base.:+ instead of + to make Pluto understand what's going on
	# oops
	if @isdefined(Coordinate)
		result = Base.:+(Coordinate(3,4), Coordinate(10,10))

		if result isa Missing
			still_missing()
		elseif !(result isa Coordinate)
			keep_working(md"Make sure that your return a `Coordinate`. 🧭")
		elseif result.x != 13 || result.y != 14
			keep_working()
		else
			correct()
		end
	end
end

# ╔═╡ 0b901714-0970-11eb-0b6a-ebe739db8037
not_defined(variable_name) = Markdown.MD(Markdown.Admonition("danger", "Oopsie!", [md"Make sure that you define a variable called **$(Markdown.Code(string(variable_name)))**"]))

# ╔═╡ 66663fcc-0a58-11eb-3568-c1f990c75bf2
if !@isdefined(origin)
	not_defined(:origin)
else
	let
		if origin isa Missing
			still_missing()
		elseif !(origin isa Coordinate)
			keep_working(md"Make sure that `origin` is a `Coordinate`.")
		else
			if origin == Coordinate(0,0)
				correct()
			else
				keep_working()
			end
		end
	end
end

# ╔═╡ ad1253f8-0a34-11eb-265e-fffda9b6473f
if !@isdefined(make_tuple)
	not_defined(:make_tuple)
else
	let
		result = make_tuple(Coordinate(2,1))
		if result isa Missing
			still_missing()
		elseif !(result isa Tuple)
			keep_working(md"Make sure that you return a `Tuple`, like so: `return (1, 2)`.")
		else
			if result == (2,1)
				correct()
			else
				keep_working()
			end
		end
	end
end

# ╔═╡ 058e3f84-0a34-11eb-3f87-7118f14e107b
if !@isdefined(trajectory)
	not_defined(:trajectory)
else
	let
		c = Coordinate(8,8)
		t = trajectory(c, 100)

		if t isa Missing
			still_missing()
		elseif !(t isa Vector)
			keep_working(md"Make sure that you return a `Vector`.")
		elseif !(all(x -> isa(x, Coordinate), t))
			keep_working(md"Make sure that you return a `Vector` of `Coordinate`s.")
		else
			if length(t) != 100
				almost(md"Make sure that you return `n` elements.")
			elseif 1 < length(Set(t)) < 90
				correct()
			else
				keep_working(md"Are you sure that you chose each step randomly?")
			end
		end
	end
end

# ╔═╡ 4fac0f36-0a59-11eb-03d0-632dc9db063a
if !@isdefined(initialize)
	not_defined(:initialize)
else
	let
		N = 200
		result = initialize(N, 1)

		if result isa Missing
			still_missing()
		elseif !(result isa Vector) || length(result) != N
			keep_working(md"Make sure that you return a `Vector` of length `N`.")
		elseif any(e -> !(e isa Agent), result)
			keep_working(md"Make sure that you return a `Vector` of `Agent`s.")
		elseif length(Set(result)) != N
			keep_working(md"Make sure that you create `N` **new** `Agent`s. Do not repeat the same agent multiple times.")
		elseif sum(a -> a.status == S, result) == N-1 && sum(a -> a.status == I, result) == 1
			if 8 <= length(Set(a.position for a in result)) <= 9
				correct()
			else
				keep_working(md"The coordinates are not correctly sampled within the box.")
			end
		else
			keep_working(md"`N-1` agents should be Susceptible, 1 should be Infectious.")
		end
	end
end

# ╔═╡ d5cb6b2c-0a66-11eb-1aff-41d0e502d5e5
bigbreak = html"

";

# ╔═╡ fcafe15a-0a66-11eb-3ed7-3f8bbb8f5809
bigbreak

# ╔═╡ ed2d616c-0a66-11eb-1839-edf8d15cf82a
bigbreak

# ╔═╡ e84e0944-0a66-11eb-12d3-e12ae10f39a6
bigbreak

# ╔═╡ e0baf75a-0a66-11eb-0562-938b64a473ac
bigbreak

# ╔═╡ Cell order:
# ╟─19fe1ee8-0970-11eb-2a0d-7d25e7d773c6
# ╟─1bba5552-0970-11eb-1b9a-87eeee0ecc36
# ╟─49567f8e-09a2-11eb-34c1-bb5c0b642fe8
# ╟─181e156c-0970-11eb-0b77-49b143cc0fc0
# ╠═1f299cc6-0970-11eb-195b-3f951f92ceeb
# ╟─2848996c-0970-11eb-19eb-c719d797c322
# ╠═2b37ca3a-0970-11eb-3c3d-4f788b411d1a
# ╠═2dcb18d0-0970-11eb-048a-c1734c6db842
# ╟─69d12414-0952-11eb-213d-2f9e13e4b418
# ╠═fcafe15a-0a66-11eb-3ed7-3f8bbb8f5809
# ╟─3e54848a-0954-11eb-3948-f9d7f07f5e23
# ╟─3e623454-0954-11eb-03f9-79c873d069a0
# ╠═0ebd35c8-0972-11eb-2e67-698fd2d311d2
# ╟─027a5f48-0a44-11eb-1fbf-a94d02d0b8e3
# ╠═b2f90634-0a68-11eb-1618-0b42f956b5a7
# ╟─66663fcc-0a58-11eb-3568-c1f990c75bf2
# ╟─3e858990-0954-11eb-3d10-d10175d8ca1c
# ╠═189bafac-0972-11eb-1893-094691b2073c
# ╟─ad1253f8-0a34-11eb-265e-fffda9b6473f
# ╟─73ed1384-0a29-11eb-06bd-d3c441b8a5fc
# ╠═96707ef0-0a29-11eb-1a3e-6bcdfb7897eb
# ╠═b0337d24-0a29-11eb-1fab-876a87c0973f
# ╟─9c9f53b2-09ea-11eb-0cda-639764250cee
# ╠═e24d5796-0a68-11eb-23bb-d55d206f3c40
# ╠═ec8e4daa-0a2c-11eb-20e1-c5957e1feba3
# ╟─e144e9d0-0a2d-11eb-016e-0b79eba4b2bb
# ╠═ec576da8-0a2c-11eb-1f7b-43dec5f6e4e7
# ╟─71c358d8-0a2f-11eb-29e1-57ff1915e84a
# ╠═5278e232-0972-11eb-19ff-a1a195127297
# ╟─71c9788c-0aeb-11eb-28d2-8dcc3f6abacd
# ╠═69151ce6-0aeb-11eb-3a53-290ba46add96
# ╟─3eb46664-0954-11eb-31d8-d9c0b74cf62b
# ╠═edf86a0e-0a68-11eb-2ad3-dbf020037019
# ╠═00f01a64-2993-11eb-2431-47731cdf2d92
# ╠═44107808-096c-11eb-013f-7b79a90aaac8
# ╟─87ea0868-0a35-11eb-0ea8-63e27d8eda6e
# ╟─058e3f84-0a34-11eb-3f87-7118f14e107b
# ╠═478309f4-0a31-11eb-08ea-ade1755f53e0
# ╠═51788e8e-0a31-11eb-027e-fd9b0dc716b5
# ╠═3ebd436c-0954-11eb-170d-1d468e2c7a37
# ╠═dcefc6fe-0a3f-11eb-2a96-ddf9c0891873
# ╟─b4d5da4a-09a0-11eb-1949-a5807c11c76c
# ╠═0237ebac-0a69-11eb-2272-35ea4e845d84
# ╠═ad832360-0a40-11eb-2857-e7f0350f3b12
# ╟─b4ed2362-09a0-11eb-0be9-99c91623b28f
# ╠═0665aa3e-0a69-11eb-2b5d-cd718e3c7432
# ╠═475f01b2-2999-11eb-3993-b5666bd60bfe
# ╟─ed2d616c-0a66-11eb-1839-edf8d15cf82a
# ╟─3ed06c80-0954-11eb-3aee-69e4ccdc4f9d
# ╠═35537320-0a47-11eb-12b3-931310f18dec
# ╠═cf2f3b98-09a0-11eb-032a-49cc8c15e89c
# ╟─814e888a-0954-11eb-02e5-0964c7410d30
# ╠═0cfae7ba-0a69-11eb-3690-d973d70e47f4
# ╠═1d0f8eb4-0a46-11eb-38e7-63ecbadbfa20
# ╟─4fac0f36-0a59-11eb-03d0-632dc9db063a
# ╠═e0b0880c-0a47-11eb-0db2-f760bbbf9c11
# ╠═b5a88504-0a47-11eb-0eda-f125d419e909
# ╠═87a4cdaa-0a5a-11eb-2a5e-cfaf30e942ca
# ╠═49fa8092-0a43-11eb-0ba9-65785ac6a42f
# ╠═1ccc961e-0a69-11eb-392b-915be07ef38d
# ╠═1f96c80a-0a46-11eb-0690-f51c60e57c3f
# ╟─f953e06e-099f-11eb-3549-73f59fed8132
# ╠═e6dd8258-0a4b-11eb-24cb-fd5b3554381b
# ╠═de88b530-0a4b-11eb-05f7-85171594a8e8
# ╠═80f39140-0aef-11eb-21f7-b788c5eab5c9
# ╠═ee5c548e-3d11-11eb-0367-6555ca14c8a7
# ╠═d1bcd5c4-0a4b-11eb-1218-7531e367a7ff
# ╠═34778744-0a5f-11eb-22b6-abe8b8fc34fd
# ╠═2d5f2b1e-3d14-11eb-2d35-7d35f97e9412
# ╠═24fe0f1a-0a69-11eb-29fe-5fb6cbf281b8
# ╟─1fc3271e-0a45-11eb-0e8d-0fd355f5846b
# ╠═18552c36-0a4d-11eb-19a0-d7d26897af36
# ╠═4e7fd58a-0a62-11eb-1596-c717e0845bd5
# ╠═778c2490-0a62-11eb-2a6c-e7fab01c6822
# ╟─e964c7f0-0a61-11eb-1782-0b728fab1db0
# ╠═4d83dbd0-0a63-11eb-0bdc-757f0e721221
# ╠═62230726-3d19-11eb-17d3-e7aab3048862
# ╠═ef27de84-0a63-11eb-177f-2197439374c5
# ╟─8475baf0-0a63-11eb-1207-23f789d00802
# ╟─201a3810-0a45-11eb-0ac9-a90419d0b723
# ╠═cd931cf2-3d33-11eb-0565-abfc49c2e4c0
# ╠═e5040c9e-0a65-11eb-0f45-270ab8161871
# ╠═a3f180f4-3d34-11eb-3c3e-0b191891eeaa
# ╠═90794bba-3d34-11eb-2f70-0da9e08ea54d
# ╠═f55659d2-3d35-11eb-3675-bfc8d6d0d687
# ╟─f9b9e242-0a53-11eb-0c6a-4d9985ef1687
# ╟─2031246c-0a45-11eb-18d3-573f336044bf
# ╠═63dd9478-0a45-11eb-2340-6d3d00f9bb5f
# ╠═269955e4-0a46-11eb-02cc-1946dc918bfa
# ╠═4d4548fe-0a66-11eb-375a-9313dc6c423d
# ╟─20477a78-0a45-11eb-39d7-93918212a8bc
# ╠═601f4f54-0a45-11eb-3d6c-6b9ec75c6d4a
# ╠═83c7d300-3eb7-11eb-22d3-0dbc0dc0bfa4
# ╠═b9bc85c6-3e0d-11eb-3b70-2bf813cbbb01
# ╠═b1b1afda-0a66-11eb-2988-752405815f95
# ╟─e84e0944-0a66-11eb-12d3-e12ae10f39a6
# ╟─05c80a0c-09a0-11eb-04dc-f97e306f1603
# ╟─b53d5608-0a41-11eb-2325-016636a22f71
# ╠═1b5e72c6-0a42-11eb-3884-a377c72270c7
# ╟─c704ea4c-0aec-11eb-2f2c-859c954aa520
# ╠═e97e39aa-0a5d-11eb-3d5f-f90a0acfe5a2
# ╟─b554b654-0a41-11eb-0e0d-e57ff68ced33
# ╠═40c1c1d6-0a69-11eb-3913-59e9b9ec4332
# ╠═a47b5dd8-3eba-11eb-1ce3-f574ecf41690
# ╟─18ac9926-0aed-11eb-034f-e9849b71c9ac
# ╟─b56ba420-0a41-11eb-266c-719d39580fa9
# ╠═465e918a-0a69-11eb-1b59-01150b4b0f36
# ╠═a0f5cbf6-3ebf-11eb-3200-8f8286079f4d
# ╟─a885bf78-0a5c-11eb-2383-9d74c8765847
# ╠═1f172700-0a42-11eb-353b-87c0039788bd
# ╟─b59de26c-0a41-11eb-2c67-b5f3c7780c91
# ╠═faec52a8-0a60-11eb-082a-f5787b09d88c
# ╟─b5b4d834-0a41-11eb-1b18-1bd626d18934
# ╠═64782262-3ec3-11eb-3b62-d5481d23438b
# ╠═a83c96e2-0a5a-11eb-0e58-15b5dda7d2d2
# ╟─05fc5634-09a0-11eb-038e-53d63c3edaf2
# ╠═24c2fb0c-0a42-11eb-1a1a-f1246f3420ff
# ╟─c7649966-0a41-11eb-3a3a-57363cea7b06
# ╠═2635b574-0a42-11eb-1daa-971b2596ce44
# ╟─c77b085e-0a41-11eb-2fcb-534238cd3c49
# ╠═274fe006-0a42-11eb-1869-29193bb84957
# ╟─c792374a-0a41-11eb-1e5b-89d9de2cf1f9
# ╠═d147f7f0-0a66-11eb-2877-2bc6680e396d
# ╟─e0baf75a-0a66-11eb-0562-938b64a473ac
# ╟─0e6b60f6-0970-11eb-0485-636624a0f9d7
# ╟─0a82a274-0970-11eb-20a2-1f590be0e576
# ╟─0aa666dc-0970-11eb-2568-99a6340c5ebd
# ╟─0acaf3b2-0970-11eb-1d98-bf9a718deaee
# ╟─0afab53c-0970-11eb-3e43-834513e4632e
# ╟─0b21c93a-0970-11eb-33b0-550a39ba0843
# ╟─0b470eb6-0970-11eb-182f-7dfb4662f827
# ╟─0b6b27ec-0970-11eb-20c2-89515ee3ab88
# ╟─0b901714-0970-11eb-0b6a-ebe739db8037
# ╟─d5cb6b2c-0a66-11eb-1aff-41d0e502d5e5


text:
PyABC_Register(s, t) = begin
    r = PyObject_GetAttrString(t, "register")
    isnull(r) && return Cint(-1)
    u = PyObject_CallNice(r, PyObjectRef(s))
    Py_DecRef(r)
    isnull(u) && return Cint(-1)
    Py_DecRef(u)
    Cint(0)
end

for n in [
    :Container,
    :Hashable,
    :Iterable,
    :Iterator,
    :Reversible,
    :Generator,
    :Sized,
    :Callable,
    :Collection,
    :Sequence,
    :MutableSequence,
    :ByteString,
    :Set,
    :MutableSet,
    :Mapping,
    :MutableMapping,
    :MappingView,
    :ItemsView,
    :KeysView,
    :ValuesView,
    :Awaitable,
    :Coroutine,
    :AsyncIterable,
    :AsyncIterator,
    :AsyncGenerator,
]
    p = Symbol(:Py, n, :ABC)
    t = Symbol(p, :_Type)
    tr = Symbol(p, :__ref)
    c = Symbol(p, :_Check)
    @eval const $tr = Ref(PyNULL)
    @eval $t(doimport::Bool = true) = begin
        ptr = $tr[]
        isnull(ptr) || return ptr
        a =
            doimport ? PyImport_ImportModule("collections.abc") :
            PyImport_GetModule("collections.abc")
        isnull(a) && return a
        b = PyObject_GetAttrString(a, $(string(n)))
        Py_DecRef(a)
        isnull(b) && return b
        $tr[] = b
    end
    @eval $c(o) = begin
        t = $t(false)
        isnull(t) && return (PyErr_IsSet() ? Cint(-1) : Cint(0))
        PyObject_IsInstance(o, t)
    end
end

"""
    PyIterable_Map(f, xs::PyPtr) :: Int

Calls `f(x)` for each element `x::PyPtr` of `xs`.

The function `f` must return an integer, where -1 signals a Python error,
0 signals to stop iterating, 1 to continue.

If `f` throws a Julia error, it is caught and converted to a Python error.
In this situation, `x` is automatically decref'd, but it is up to the user
to ensure that other references created by `f` are cleared.

Return -1 on error, 0 if iteration was stopped, 1 if it reached the end.
"""
PyIterable_Map(f, xso::PyPtr) = begin
    it = PyObject_GetIter(xso)
    isnull(it) && return PyNULL
    xo = PyNULL
    try
        while true
            xo = PyIter_Next(it)
            if !isnull(xo)
                r = f(xo)
                Py_DecRef(xo)
                xo = PyNULL
                if r == -1
                    return -1
                elseif r == 0
                    return 0
                end
            elseif PyErr_IsSet()
                return -1
            else
                return 1
            end
        end
    catch err
        PyErr_SetJuliaError(err)
        -1
    finally
        Py_DecRef(it)
        Py_DecRef(xo)
    end
end

"""
    PyIterable_Collect(xs::PyPtr, xs, [skip=false]) :: Vector{T}

Convert the elements of `xs` to type `T` and collect them into a vector.
On error, an empty vector is returned.

If `skip` then elements which cannot be converted to a `T` are skipped over,
instead of raising an error. Other errors are still propagated.
"""
PyIterable_Collect(xso::PyPtr, ::Type{T}, skip::Bool = false) where {T} = begin
    xs = T[]
    r = if skip
        PyIterable_Map(xso) do xo
            r = PyObject_TryConvert(xo, eltype(xs))
            r == -1 && return -1
            r == 0 && return 1
            x = takeresult(eltype(xs))
            push!(xs, x)
            return 1
        end
    else
        PyIterable_Map(xso) do xo
            r = PyObject_Convert(xo, eltype(xs))
            r == -1 && return -1
            x = takeresult(eltype(xs))
            push!(xs, x)
            return 1
        end
    end
    r == -1 && empty!(xs)
    xs
end

PyIterable_ConvertRule_vecorset(o, ::Type{S}) where {S<:Union{Vector,Set}} = begin
    xs = S()
    r = PyIterable_Map(o) do xo
        r = PyObject_TryConvert(xo, eltype(xs))
        r == -1 && return -1
        r == 0 && return 0
        x = takeresult(eltype(xs))
        push!(xs, x)
        return 1
    end
    r == 1 && putresult(xs)
    r
end
PyIterable_ConvertRule_vecorset(o, ::Type{Vector}) =
    PyIterable_ConvertRule_vecorset(o, Vector{Python.PyObject})
PyIterable_ConvertRule_vecorset(o, ::Type{Set}) =
    PyIterable_ConvertRule_vecorset(o, Set{Python.PyObject})

PyIterable_ConvertRule_tuple(o, ::Type{S}) where {S<:Tuple} = begin
    S isa DataType || return 0
    ts = S.parameters
    if !isempty(ts) && Base.isvarargtype(ts[end])
        isvararg = true
        vartype = ts[end].body.parameters[1]
        ts = ts[1:end-1]
    else
        isvararg = false
    end
    xs = Union{ts...,isvararg ? vartype : Union{}}[]
    r = PyIterable_Map(o) do xo
        if length(xs) < length(ts)
            t = ts[length(xs)+1]
        elseif isvararg
            t = vartype
        else
            return 0
        end
        r = PyObject_TryConvert(xo, t)
        r == -1 && return -1
        r == 0 && return 0
        x = takeresult(t)
        push!(xs, x)
        return 1
    end
    r == -1 ? -1 : r == 0 ? 0 : length(xs) ≥ length(ts) ? putresult(S(xs)) : 0
end

PyIterable_ConvertRule_namedtuple(o, ::Type{NamedTuple{names, types}}) where {names, types<:Tuple} = begin
    r = PyIterable_ConvertRule_tuple(o, types)
    r == 1 ? putresult(NamedTuple{names, types}(takeresult(types))) : r
end
PyIterable_ConvertRule_namedtuple(o, ::Type{NamedTuple{names}}) where {names} = begin
    types = NTuple{length(names), Python.PyObject}
    r = PyIterable_ConvertRule_tuple(o, types)
    r == 1 ? putresult(NamedTuple{names}(takeresult(types))) : r
end
PyIterable_ConvertRule_namedtuple(o, ::Type{S}) where {S<:NamedTuple} = 0

PyIterable_ConvertRule_pair(o, ::Type{Pair{K,V}}) where {K,V} = begin
    k = Ref{K}()
    v = Ref{V}()
    i = Ref(0)
    r = PyIterable_Map(o) do xo
        if i[] == 0
            # key
            r = PyObject_TryConvert(xo, eltype(k))
            r == -1 && return -1
            r == 0 && return 0
            i[] = 1
            k[] = takeresult(eltype(k))
            return 1
        elseif i[] == 1
            # value
            r == PyObject_TryConvert(xo, eltype(v))
            r == -1 && return -1
            r == 0 && return 0
            i[] = 2
            v[] = takeresult(eltype(v))
            return 1
        else
            return 0
        end
    end
    r == -1 ? -1 : r == 0 ? 0 : i[] == 2 ? putresult(Pair{K,V}(k[], v[])) : 0
end
PyIterable_ConvertRule_pair(o, ::Type{Pair{K}}) where {K} =
    PyIterable_ConvertRule_pair(o, Pair{K,Python.PyObject})
PyIterable_ConvertRule_pair(o, ::Type{Pair{K,V} where K}) where {V} =
    PyIterable_ConvertRule_pair(o, Pair{Python.PyObject,V})
PyIterable_ConvertRule_pair(o, ::Type{Pair}) =
    PyIterable_ConvertRule_pair(o, Pair{Python.PyObject,Python.PyObject})
PyIterable_ConvertRule_pair(o, ::Type{S}) where {S<:Pair} = begin
    PyErr_SetString(
        PyExc_Exception(),
        "When converting Python iterable to Julia 'Pair', the destination type cannot be too complicated: the two types must either be fully specified or left unspecified. Got '$S'.",
    )
    return -1
end

PyMapping_ConvertRule_dict(o, ::Type{S}) where {S<:Dict} = begin
    xs = S()
    r = PyIterable_Map(o) do ko
        # get the key
        r = PyObject_TryConvert(ko, keytype(xs))
        r == -1 && return -1
        r == 0 && return 0
        k = takeresult(keytype(xs))
        # get the value
        vo = PyObject_GetItem(o, ko)
        isnull(vo) && return -1
        r = PyObject_TryConvert(vo, valtype(xs))
        Py_DecRef(vo)
        r == -1 && return -1
        r == 0 && return 0
        v = takeresult(valtype(xs))
        # done
        xs[k] = v
    end
    r == -1 ? -1 : r == 0 ? 0 : putresult(xs)
end
PyMapping_ConvertRule_dict(o, ::Type{Dict{K}}) where {K} =
    PyMapping_ConvertRule_dict(o, Dict{K,Python.PyObject})
PyMapping_ConvertRule_dict(o, ::Type{Dict{K,V} where K}) where {V} =
    PyMapping_ConvertRule_dict(o, Dict{Python.PyObject,V})
PyMapping_ConvertRule_dict(o, ::Type{Dict}) =
    PyMapping_ConvertRule_dict(o, Dict{Python.PyObject,Python.PyObject})

ext: jl


text:
#!/usr/bin/env python
# coding: utf-8

# # Propagation on the path graph

# ## Loading modules

using QSWalk
using PyPlot # for plot
using LightGraphs # for PathGraph
using LinearAlgebra # for diag method

# ## Evolution on path graph for fixed time, global vs local interaction

# Basic parameters. Note that *dim* should be an odd number. Otherwise the middle point *midpoint* is not unique. Matrix
# *adjacency* can be alternatively generated as *adjacency = spdiagm((ones(dim-1), ones(dim-1)), (-1, 1))*.

dim = 151 # odd for unique middle point
midpoint = ceil(Int, dim/2)
w = 0.5
timepoint = 50.
adjacency = adjacency_matrix(PathGraph(dim));

# Quantum stochastic walk: global interaction case. Note that *evolve_generator* requires a list of Lindblad operators, hence we need extra brackets around the second argument.

op_global = evolve_generator(adjacency, [adjacency], w)
rho_global = evolve(op_global, proj(midpoint, dim), timepoint);

# Quantum stochastic walk: local interaction case. Note that here *lind_local* is already a list of Lindbladian operators, hence we do not need extra brackets, contrary to the global interaction example.

lind_local = local_lind(adjacency)
op_local = evolve_generator(adjacency, lind_local, w)
rho_local = evolve(op_local, proj(midpoint, dim), timepoint);

# To plot the result of the cannonical measurement, we take a diagonal of the states. Note that both *rhoglobal* and *rholocal* are complex matrices, hence we need to take the real part only. Argument *positions* is generated assuming that *midpoint* (the initial position) is at 0. Note that we have very heavy tails in the global interaction case, which confirms fast propagation in this model.

positions = (collect(1:dim) .- midpoint)
plot(positions, real.(diag(rho_local)), "k")
plot(positions, real.(diag(rho_global)), "b")
xlabel("position")
ylabel("probability")
show()

# ## Evolution on path graph for fixed time, propagation analysis

# Basic parameters. Note that *dim* should be an odd number. Otherwise the middle point *midpoint* is not unique. Matrix
# *adjacency* can be alternatively generated as *adjacency = spdiagm((ones(dim-1), ones(dim-1)), (-1, 1))*.

dim = 151 #odd for unique middle point
midpoint = ceil(Int, dim/2)
w = 0.5
timepoints = collect(0:2.5:50)
adjacency = adjacency_matrix(PathGraph(dim));

# Preparation of the global and the local interaction operators.

op_global = evolve_generator(adjacency, [adjacency], w)
lind_local = local_lind(adjacency)
op_local = evolve_generator(adjacency, lind_local, w);

# Calculating resulting states. Note that comparing to global operators, density states usually are small matrices. Hence we can store those. Furthermore, note that here we use variant of *evolve* function, whcih enables us to compute the same evolutions for different time points.

global_states = evolve(op_global, proj(midpoint, dim), timepoints)
local_states = evolve(op_local, proj(midpoint, dim), timepoints);

# Calculation of the second moment.

secmoment_global = Float64[]
secmoment_local = Float64[]
positions = (collect(1:dim) .- midpoint)

for (rho_global, rho_local) = zip(global_states, local_states)
  push!(secmoment_global, sum(positions.^2 .* diag(rho_global)))
  push!(secmoment_local, sum(positions.^2 .* diag(rho_local)))
end

# Plots of the results. Note that second moment for global interaction case grows ballistically.

plot(timepoints, secmoment_local, "k")
plot(timepoints, secmoment_global, "b")
xlabel("t")
ylabel(L"\mu_2")
show()


text:
group = addgroup!(SUITE, "array")

group["construct"] = @benchmarkable CuArray{Int}(undef, 1)

cpu_mat = rand(Float32, 512, 1000)
gpu_mat = CuArray{Float32}(undef, size(cpu_mat))
group["copy"] = @benchmarkable CUDA.@sync copy($gpu_mat)

gpu_mat2 = copy(gpu_mat)
let group = addgroup!(SUITE, "copyto!")
    group["cpu_to_gpu"] = @benchmarkable CUDA.@sync copyto!($gpu_mat, $cpu_mat)
    group["gpu_to_cpu"] = @benchmarkable CUDA.@sync copyto!($cpu_mat, $gpu_mat)
    group["gpu_to_gpu"] = @benchmarkable CUDA.@sync copyto!($gpu_mat2, $gpu_mat)
end

group["fill!"] = @benchmarkable CUDA.@sync fill!($gpu_mat, 0f0)

gpu_vec = reshape(gpu_mat, length(gpu_mat))
let group = addgroup!(SUITE, "reverse")
    group["1d"] = @benchmarkable CUDA.@sync reverse($gpu_vec)
    group["2d"] = @benchmarkable CUDA.@sync reverse($gpu_mat2; dims=1)
    group["1d_inplace"] = @benchmarkable CUDA.@sync reverse!($gpu_vec)
    group["2d_inplace"] = @benchmarkable CUDA.@sync reverse!($gpu_mat2; dims=1)
end

group["broadcast"] = @benchmarkable CUDA.@sync $gpu_mat .= 0f0


text:
# Julia binding for IP Tool.
#
# Uses the IP Tool FFI to provide all of the
# IP Tool functions in a Julia module.
#
# Author:
# Tim Monfette

module Iptool
export valid_ipv4, ipv4_ipv6, ipv4_bin, valid_ipv6, ipv6_ipv4, ipv6_bin

# Call free_string
function free_string(toFree::Cstring)
  ccall((:free_string, "./lib/libiptool"),
        Void, (Cstring,), toFree)
end

# Wrapper for valid_ipv4
function valid_ipv4(userInput::AbstractString)
  result = ccall((:valid_ipv4, "./lib/libiptool"),
                 Cstring, (Cstring,), userInput)

  jstring = unsafe_string(result)
  free_string(result)
  jstring
end

# Wrapper for ipv4_ipv6
function ipv4_ipv6(userInput::AbstractString)
  result = ccall((:ipv4_ipv6, "./lib/libiptool"),
                 Cstring, (Cstring,), userInput)

  jstring = unsafe_string(result)
  free_string(result)
  jstring
end

# Wrapper for ipv4_bin
function ipv4_bin(userInput::AbstractString)
  result = ccall((:ipv4_bin, "./lib/libiptool"),
                 Cstring, (Cstring,), userInput)

  jstring = unsafe_string(result)
  free_string(result)
  jstring
end

# Wrapper for valid_ipv6
function valid_ipv6(userInput::AbstractString)
  result = ccall((:valid_ipv6, "./lib/libiptool"),
                 Cstring, (Cstring,), userInput)

  jstring = unsafe_string(result)
  free_string(result)
  jstring
end

# Wrapper for ipv6_ipv4
function ipv6_ipv4(userInput::AbstractString)
  result = ccall((:ipv6_ipv4, "./lib/libiptool"),
                 Cstring, (Cstring,), userInput)

  jstring = unsafe_string(result)
  free_string(result)
  jstring
end

# Wrapper for ipv6_bin
function ipv6_bin(userInput::AbstractString)
  result = ccall((:ipv6_bin, "./lib/libiptool"),
                 Cstring, (Cstring,), userInput)

  jstring = unsafe_string(result)
  free_string(result)
  jstring
end

end


text:
# see correct result with slanted (narrow) L1 modes at y=+-5.  L1 should not be round blobs
# https://github.com/JuliaRobotics/RoME.jl/pull/434#issuecomment-817246038

using RoME
using Test

##

@testset "Test correlation induced by partials" begin

##

# start with an empty factor graph object
N = 200
fg = initfg()
getSolverParams(fg).N = N

getSolverParams(fg).useMsgLikelihoods = true

addVariable!(fg, :x0, Pose2)
addVariable!(fg, :l1, Point2)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0; 0; 0], diagm([0.3;0.3;0.3].^2)) ))
ppr = Pose2Point2Range(MvNormal([7.3], diagm([0.3].^2)))
addFactor!(fg, [:x0; :l1], ppr )

addVariable!(fg, :x1, Pose2)
pp = Pose2Pose2(MvNormal([9.8;0;0.8], diagm([0.3;0.3;0.05].^2)))
ppr = Pose2Point2Range(MvNormal([6.78], diagm([0.3].^2)))
addFactor!(fg, [:x0; :x1], pp )
addFactor!(fg, [:x1; :l1], ppr )

##

tree, smt, hist = solveTree!(fg)

## check that stuff is where it should be

L1_ = getPoints(getBelief(fg, :l1))

@test 0.3*N < sum(L1_[2,:] .< 0)
@test 0.3*N < sum(0 .< L1_[2,:])

L1_n = L1_[:, L1_[2,:] .< 0]
L1_p = L1_[:, 0 .< L1_[2,:]]

mvn = fit(MvNormal, L1_n)
mvp = fit(MvNormal, L1_p)

# check diagonal structure for correlation
@test isapprox(mvn.Σ.mat[1,1], 1.7, atol=1.2)
@test isapprox(mvn.Σ.mat[2,2], 1.7, atol=1.2)
@test isapprox(mvn.Σ.mat[1,2], 1.1, atol=1.2)

@test isapprox(mvp.Σ.mat[1,1], 1.7, atol=1.2)
@test isapprox(mvp.Σ.mat[2,2], 1.7, atol=1.2)
@test isapprox(mvp.Σ.mat[1,2], -1.1, atol=1.2)

# sanity check for symmetry
@test mvn.Σ.mat - mvn.Σ.mat' |> norm < 0.01

# test means in the right location
@test isapprox(mvn.μ[1], 5.4, atol=1.0)
@test isapprox(mvn.μ[2], -4.8, atol=0.75)

@test isapprox(mvp.μ[1], 5.4, atol=1.0)
@test isapprox(mvp.μ[2], 4.8, atol=0.75)

##

end

##

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)

##

# plotSLAM2D(fg)

# plotKDE(fg, :l1, levels=3)

##


text:
using Mux
using Base.Test
using Lazy
import Requests

@test Mux.notfound()(d())[:status] == 404

# Test basic server
@app test = (
  Mux.defaults,
  page(respond("
Hello World!
")),
  page("/about", respond("
Boo!
")),
  page("/user/:user", req -> "
Hello, $(req[:params][:user])!
"),
  Mux.notfound())
serve(test)
@test Requests.text(Requests.get("http://localhost:8000")) ==
            "
Hello World!
"
@test Requests.text(Requests.get("http://localhost:8000/about")) ==
            "
Boo!
"
@test Requests.text(Requests.get("http://localhost:8000/user/julia")) ==
            "
Hello, julia!
"


text:
# CUDA discovery and dependency management

using Pkg, Pkg.Artifacts
using Libdl

## global state

const __version = Ref{VersionNumber}()

"""
    version()

Returns the version of the CUDA toolkit in use.
"""
version() = @after_init(__version[])

"""
    release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
release() = @after_init(VersionNumber(__version[].major, __version[].minor))

const __nvdisasm = Ref{String}()
const __libdevice = Ref{String}()
const __libcudadevrt = Ref{String}()
const __libcupti = Ref{Union{Nothing,String}}()
const __libnvtx = Ref{Union{Nothing,String}}()

nvdisasm() = @after_init(__nvdisasm[])
libdevice() = @after_init(__libdevice[])
libcudadevrt() = @after_init(__libcudadevrt[])
function libcupti()
    @after_init begin
        @assert has_cupti() "This functionality is unavailable as CUPTI is missing."
        __libcupti[]
    end
end
function libnvtx()
    @after_init begin
        @assert has_nvtx() "This functionality is unavailable as NVTX is missing."
        __libnvtx[]
    end
end

export has_cupti, has_nvtx
has_cupti() = @after_init(__libcupti[]) !== nothing
has_nvtx() = @after_init(__libnvtx[]) !== nothing

## discovery

# NOTE: we don't use autogenerated JLLs, because we have multiple artifacts and need to
#       decide at run time (i.e. not via package dependencies) which one to use.
const cuda_artifacts = Dict(
    v"10.2" => ()->artifact"CUDA10.2",
    v"10.1" => ()->artifact"CUDA10.1",
    v"10.0" => ()->artifact"CUDA10.0",
    v"9.2"  => ()->artifact"CUDA9.2",
    v"9.0"  => ()->artifact"CUDA9.0",
)

# try use CUDA from an artifact
function use_artifact_cuda()
    @debug "Trying to use artifacts..."

    # select compatible artifacts
    if haskey(ENV, "JULIA_CUDA_VERSION")
        wanted_version = VersionNumber(ENV["JULIA_CUDA_VERSION"])
        filter!(((version,artifact),) -> version == wanted_version, cuda_artifacts)
    else
        driver_version = CUDAdrv.release()
        filter!(((version,artifact),) -> version <= driver_version, cuda_artifacts)
    end

    # download and install
    artifact = nothing
    for release in sort(collect(keys(cuda_artifacts)); rev=true)
        try
            artifact = (release=release, dir=cuda_artifacts[release]())
            break
        catch
        end
    end
    if artifact == nothing
        @debug "Could not find a compatible artifact."
        return false
    end

    # utilities to look up stuff in the artifact (at known locations, so not using CUDAapi)
    get_binary(name) = joinpath(artifact.dir, "bin", Sys.iswindows() ? "$name.exe" : name)
    function get_library(name)
        filename = if Sys.iswindows()
            "$name.dll"
        elseif Sys.isapple()
            "lib$name.dylib"
        else
            "lib$name.so"
        end
        joinpath(artifact.dir, Sys.iswindows() ? "bin" : "lib", filename)
    end
    get_static_library(name) = joinpath(artifact.dir, "lib", Sys.iswindows() ? "$name.lib" : "lib$name.a")
    get_file(path) = joinpath(artifact.dir, path)

    __nvdisasm[] = get_binary("nvdisasm")
    @assert isfile(__nvdisasm[])
    __version[] = parse_toolkit_version(__nvdisasm[])

    # Windows libraries are tagged with the CUDA release
    long = "$(artifact.release.major)$(artifact.release.minor)"
    short = artifact.release >= v"10.1" ? string(artifact.release.major) : long

    __libcupti[] = get_library(Sys.iswindows() ? "cupti64_$long" : "cupti")
    @assert isfile(__libcupti[])
    __libnvtx[] = get_library(Sys.iswindows() ? "nvToolsExt64_1" : "nvToolsExt")
    @assert isfile(__libnvtx[])

    __libcudadevrt[] = get_static_library("cudadevrt")
    @assert isfile(__libcudadevrt[])
    __libdevice[] = get_file(joinpath("share", "libdevice", "libdevice.10.bc"))
    @assert isfile(__libdevice[])

    @debug "Using CUDA $(__version[]) from an artifact at $(artifact.dir)"
    return true
end

# try to use CUDA from a local installation
function use_local_cuda()
    @debug "Trying to use local installation..."

    cuda_dirs = find_toolkit()

    let path = find_cuda_binary("nvdisasm", cuda_dirs)
        if path === nothing
            @debug "Could not find nvdisasm"
            return false
        end
        __nvdisasm[] = path
    end

    cuda_version = parse_toolkit_version(__nvdisasm[])
    __version[] = cuda_version

    cupti_dirs = map(dir->joinpath(dir, "extras", "CUPTI"), cuda_dirs) |> x->filter(isdir,x)
    __libcupti[] = find_cuda_library("cupti", [cuda_dirs; cupti_dirs], [cuda_version])
    __libnvtx[] = find_cuda_library("nvtx", cuda_dirs, [v"1"])

    let path = find_libcudadevrt(cuda_dirs)
        if path === nothing
            @debug "Could not find libcudadevrt"
            return false
        end
        __libcudadevrt[] = path
    end
    let path = find_libdevice(cuda_dirs)
        if path === nothing
            @debug "Could not find libdevice"
            return false
        end
        __libdevice[] = path
    end

    @debug "Found local CUDA $(cuda_version) at $(join(cuda_dirs, ", "))"
    return true
end

function __configure_dependencies__(show_reason::Bool)
    found = false

    if parse(Bool, get(ENV, "JULIA_CUDA_USE_BINARYBUILDER", "true"))
        found = use_artifact_cuda()
    end

    if !found
        found = use_local_cuda()
    end

    if !found
        show_reason && @error "Could not find a suitable CUDA installation"
    end

    return found
end


text:
module AntiIdClosedGroup

using ClosedGroupFunctions
using Zd_Arithmetics
using SparseArrays
using LinearAlgebra
using Combinatorics
using Bijections

import Zd_Arithmetics: ℤ₂, ℤ₄
import Base: inv

export

    id,
    anti_id,

    po2ai,
    ai2po,
    po2matrix,
    po2latex,
    all_pos2latex,
    inv,

    first_po,
    next_po,
    construct_generators,

    perm_rows,
    perm_cols,
    print_permutation_cycles,
    apply_permutation,

    ord,
    tr_of_power,
    rank_of_power_plus_id,
    ℤ₂ⁿ,
    ker_g_plus_id,
    quadratic_form,
    GArf,
    generalized_arf

include("anti-identity.jl")
include("notation_converters.jl")
include("generators.jl")
include("permutations.jl")
include("conjugacy_invariants.jl")

end # module


text:
module PackmolInputCreator

using PDBTools

export find_x
export convert_c
export write_input
export interpolate

# conversion factor from mL/mol to Å^3/molecule
const CMV = 1e24 / 6.02e23

# conversion factor from mol/L to molecules/Å^3
const CMC = 6.02e23 / 1e27

"""

Interpolate a value from a (x,y) data matrix

"""
function interpolate(x,ρ)
  i = findfirst(d -> d > x, ρ[:,1])
  dρdx = (ρ[i,2]-ρ[i-1,2])/(ρ[i,1]-ρ[i-1,1])
  d = ρ[i-1,2] + dρdx * (x - ρ[i-1,1])
  return d
end

"""

Find by bisection which is the molar fraction that is consistent
with a desired volume fraction (%), given a data table of densities
as a function of the molar fractions.

"""
function find_x(
  vv, cossolvent_mass, density_pure, densities::Matrix;
  tol = 1e-5, maxit=1000
)

  xl = densities[1,1]
  xr = densities[end,1]
  if xl < xr
    increasing = true
  else
    increasing = false
  end
  x = (xr-xl)/2
  xnew = convert_c(
    vv,"%vv" => "x", density=interpolate(x,densities),
    molar_mass=cossolvent_mass, density_pure=density_pure
  )
  it = 0
  while abs(xnew-x) > tol
    if xnew > x
      increasing ? xl = x : xr = x
    else
      increasing ? xr = x : xl = x
    end
    x = (xr+xl)/2
    xnew = convert_c(
      vv,"%vv" => "x", density=interpolate(x,densities),
      molar_mass=cossolvent_mass, density_pure=density_pure
    )
    it += 1
    it > maxit && error("Maximum number of iterations achieved.")
  end
  return x
end

"""

Convert concentrations.

"""
function convert_c(
  cin,units;
  density=nothing, # solution
  density_pure=nothing, # cossolvent
  molar_mass=18., # cossolvent
  molar_mass_water=18.,
  density_water=1.0
)

  # Why not allow this
  units[1] == units[2] && return cin

  ρ = density
  ρc = density_pure
  ρw = density_water
  Mw = molar_mass_water
  Mc = molar_mass

  # nc and nw are the molar concentrations

  if units[1] == "%vv"
    if ρc == nothing
      error("Density of pure solvent is required to convert from %vv.")
    end
    vv = cin/100
    nc = ( ρc*vv / Mc )
    if units[2] == "x"
      if ρ == nothing
        error("Density of solution is required to convert to molar fraction.")
      end
      nw = ( ρ - nc*Mc ) / Mw
      return nc/(nc + nw)
    end
    if units[2] == "mol/L"
      return 1000*nc
    end
  end

  if units[1] == "x"
    if ρ == nothing
      error("Density of solution is required to convert from molar fraction.")
    end
    x = cin
    nc = ρ / ( Mc + Mw*(1-x)/x )
    if units[2] == "mol/L"
      return 1000*nc
    end
    if units[2] == "%vv"
      if ρc == nothing
        error("Density of pure solvent is required to convert to %vv.")
      end
      nw = nc*(1-x)/x
      Vc = nc * Mc / ρc
      vv = 100*Vc
      return vv
    end
  end

  if units[1] == "mol/L"
    if ρ == nothing
      error("Density of solution is required to convert from molarity.")
    end
    nc = cin/1000
    nw = ( ρ - nc*Mc ) / Mw
    if units[2] == "x"
      return nc/(nc + nw)
    end
    if units[2] == "%vv"
      if ρc == nothing
        error("Density of pure solvent is required to convert to %vv.")
        return nothing
      end
      Vc = nc * Mc / ρc
      vv = 100*Vc
      return vv
    end
  end

end

"""

Function that generates an input file for Packmol. By default, the concentrations is given in mol/L, but it can also be given in molar fraction "x" or volume percentage "%vv", using `cunit="x"` or `cunit="%vv"`.

"""
function write_input(
  pdbfile::String, solvent_file::String, concentration::Real, box_side::Real;
  water_file="tip4p2005.pdb",
  density=1.0,
  density_pure_solvent=nothing,
  cunit="mol/L",
  packmol_input="box.inp",
  packmol_output="system.pdb"

)

  protein = readPDB(pdbfile)
  cossolvent = readPDB(solvent_file)

  # molar masses (g/mol)
  Mp = mass(protein)
  Mc = mass(cossolvent)
  Mw = 18.01528

  # aliases for clearer formulas
  ρ = density
  ρc = density_pure_solvent

  # Convert concentration to mol/L
  cc_mol = convert_c(
    concentration, cunit => "mol/L",
    density=ρ,
    density_pure=ρc,
    molar_mass=Mc
  )
  c_vv = convert_c(
    concentration, cunit => "%vv",
    density=ρ,
    density_pure=ρc,
    molar_mass=Mc
  )
  c_x = convert_c(
    concentration, cunit => "x",
    density=ρ,
    density_pure=ρc,
    molar_mass=Mc
  )

  # Convert cossolvent concentration in molecules/Å³
  cc = CMC*cc_mol

  # Box volume (Å³)
  vbox = box_side^3

  # Solution volume (vbox - vprotein)
  vs = vbox - CMV*Mp/ρ

  # number of cossolvent molecules
  nc = round(Int,cc*vs)

  #number of water molecules
  nw = round(Int,(ρ*vs/CMV - nc*Mc)/Mw)

  # Final density of the solution
  ρ = CMV*(Mc*nc + Mw*nw)/vs

  # Final cossolvent concentration (mol/L)
  cc_f = 1000*(nc/vs)*CMV

  # Final water concentration (mol/L)
  cw_f = 1000*(nw/vs)*CMV

  # Final recovered concentration in %vv
  vv = 100*CMV*(nc*Mc/ρc)/vs

  println(
    """
    Summary:
    ========
    Target concentration = $cc_mol mol/L
                         = $c_vv %vv
                         = $c_x x
                         = $cc molecules/Å³

    Box volume = $vbox Å³
    Solution volume = $vs Å³

    Density = $ρ g/mL
    Protein molar mass = $Mp g/mol
    Cossolvent molar mass = $Mc g/mol

    Number of cossolvent molecules = $nc molecules
    Number of water molecules = $nw molecules

    Final cossolvent concentration = $cc_f mol/L
    Final water concentration = $cw_f mol/L
                              = $(CMC*cw_f) molecules/Å³
    Final solvent density = $ρ g/mL
    Final %vv concentration = $vv %
    Final molar fraction = $(nc/(nc+nw))
    """)

  l = box_side/2
  open(packmol_input,"w") do io
    println(io,
      """
      tolerance 2.0
      output $packmol_output
      add_box_sides 1.0
      filetype pdb
      seed -1

      structure $pdbfile
        number 1
        center
        fixed 0. 0. 0. 0. 0. 0.
      end structure

      structure $water_file
        number $nw
        inside box -$l -$l -$l $l $l $l
      end structure
      """)
    if nc > 0
      println(io,
        """
        structure $solvent_file
          number $nc
          inside box -$l -$l -$l $l $l $l
        end structure
        """)
   end
 end
 println("Wrote file: $packmol_input")

end # function write_input

end # module


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using FoodCachingExperiments, FoodCachingModels, Statistics, DataFrames, PGFPlotsX,
      FoodCachingPlotting
import FoodCachingExperiments: Experiment, _abc
import PlanningVsCompensatoryCaching
const PC = PlanningVsCompensatoryCaching
import StatsBase: sample

postprocess(::Any, data) = data
function postprocess(::Experiment{:Amodio20_exp2}, data)
    cached = []
    for row in 1:nrow(data)
        push!(cached, circshift([data.A[row], data.B[row], data.C[row]],
                                -(_abc(row, 1)-1)))
    end
    data, cached
end
function postprocess(::Experiment{:Amodio20_exp1}, data)
    cached = []
    for row in 1:nrow(data) ÷ 2
        push!(cached, hcat([circshift([data.A[2*row-i],
                                       data.B[2*row-i],
                                       data.C[2*row-i]],
                                      -(_abc(row, 1)-1))
                            for i in (isodd(row) ? (1, 0) : (0, 1))]...))
    end
    data, cached
end
function predict(fp; N = 500, experiment = :Amodio20_exp2)
    ms = rand(fp, N)
    expe = Experiment{experiment}(nothing, nothing, "", "", DataFrame(), Dict(), 0, [])
    data = run!(expe, ms)
    postprocess(expe, data)
end

function posterior(cached, K, ms; only_data = false)
    N = length(cached)
    d_fn = cached[sample(2:2:N, K, replace = false)]
    only_data || (lp_fn = PC.compute_posterior(d_fn, ms, group = :FN, unnormalizedlog = true))
    d_nf = cached[sample(1:2:N, K, replace = false)]
    only_data ||  (lp_nf = PC.compute_posterior(d_nf, ms, group = :NF, unnormalizedlog = true))
    only_data ? 0 : (x -> x ./ sum(x))(exp.(lp_fn .+ lp_nf)), mean(d_fn), mean(d_nf)
end
function posterior2(cached, K, ms)
    N = length(cached)
    d = cached[sample(1:N, K, replace = false)]
    p = PC.compute_posterior(d, ms)
    p, mean(d)
end

function plotcachingpattern(res)
    @pgf Axis({xticklabels = [raw"$K_1$", raw"$K_2$", raw"$K_3$"],
                xtick = 1:3, height = "3cm", width = "3cm",
                font = raw"\scriptsize",
                ylabel = "items cached",
#                 legend_entries = ["group 1", "group 2"]
},
               Plot({"mblue", mark = "square"}, Table(1:3, mean(x[2] for x in res[3]))),
               Plot({"mbrown", mark = "square"}, Table(1:3, mean(x[3] for x in res[3]))))
end
function plot_bayesratioh(bf; colors = ["red", "green!60!black", "blue"],
                              data = nothing,
                              legend = ["6 birds", "12 birds", "18 birds"],
                              opt = @pgf({}),
    )
    p1 = boxplot([bf],
            options = [PGFPlotsX.Options(:draw => c) for c in colors],
            labels = [],
            axisopt = merge(@pgf({ylabel = "Bayes Factor",
#                             ymin = 1e-8, ymax = 1e3,
#                             ytick = [1e-6, 1e-3, 1, 1e3],
                            height = "4.5cm", xtick = [.75, 1, 1.25],
                            font = raw"\scriptsize",
                            width = "3cm",
                            "legend image code/.code = {\\draw [#1,only marks] plot coordinates { (0cm, 0cm) (0.1cm, 0cm) (0.2cm, 0cm)};}",
                            legend_entries = legend
                           }), opt)
           );
    push!(p1.contents, raw"\draw[dotted] (0, 1) -- (4, 1);");
    data === nothing || push!(p1.contents, "\\draw[] (0, $data) -- (4, $data);");
    p1
end
function plotcachingpattern1(res)
    @pgf Axis({xticklabels = [raw"$K_1$", raw"$K_2$", raw"$K_3$"],
                xtick = 1:3, height = "3cm", width = "3cm",
                font = raw"\scriptsize",
                ylabel = "items cached",
#                 legend_entries = ["group 1", "group 2"]
},
              Plot({"mblue", mark = "square"}, Table(1:3, mean(res)[:, 1])),
              Plot({"mbrown", mark = "square"}, Table(1:3, mean(res)[:, 2])))
end
function plotcachingpattern_hyp(res)
    @pgf Axis({xticklabels = [raw"$K_1$", raw"$K_2$", raw"$K_3$"],
                xtick = 1:3, height = "3cm", width = "3cm",
                ymin = 0, ymax = 7,
                font = raw"\scriptsize",
                ylabel = "items cached",
#                 legend_entries = ["group 1", "group 2"]
},
               Plot({"mblue", mark = "square"}, Table(1:2, mean(x[2] for x in res[3]))),
              )
end
function plot_bayesratio(bf;
        colors = ["red", "green!60!black", "blue"],
        labels = [raw"$M_0/M_1$", raw"$M_0/M_2$", raw"$M_1/M_2$"],
        legend = ["6 birds", "12 birds", "18 birds"],
        opt = @pgf({}),
    )
    p1 = boxplot(bf,
                 options = [PGFPlotsX.Options(:draw => c) for c in colors],
                 labels = [],
            axisopt = merge(@pgf({ylabel = "Bayes Factor",
#                             ymin = 1e-6, ymax = 1e5,
#                             xmax = 2.5,
#                             xmin = 1.5,
                            height = "4.5cm", xtick = [.75, 1, 1.25],
                            font = raw"\scriptsize",
                            width = "3cm",
#                             ytick = [1e-12, 1e-3, 1e-0, 1e3],
                            "legend image code/.code = {\\draw [#1,only marks] plot coordinates { (0cm, 0cm) (0.1cm, 0cm) (0.2cm, 0cm)};}",
                            legend_entries = legend
                           }), opt));
    push!(p1.contents, raw"\draw[dotted] (0, 1) -- (4, 1);");
    push!(p1.contents, raw"\draw[] (0, 80) -- (4, 80);");
    p1
end

include(joinpath(@__DIR__, "paths.jl"))
model1 = load(joinpath(datapath, "ReplayAndPlan_Raby07_breakfastchoice"))
model2 = load(joinpath(datapath, "PlasticCaching_Raby07_breakfastchoice"))

data11, cached11 = predict(model1, N = 5000, experiment = :Amodio20_exp1);
data21, cached21 = predict(model2, N = 5000, experiment = :Amodio20_exp1);
data1, cached1 = predict(model1, N = 5000, experiment = :Amodio20_exp2);
data2, cached2 = predict(model2, N = 5000, experiment = :Amodio20_exp2);
datah1 = predict(model1, N = 5000, experiment = :Hypothetical);
datah2 = predict(model2, N = 5000, experiment = :Hypothetical);

ms = [
      (bird_indep = true, compartement_indep = true,
       food_indep = true, hyp = nothing),
      (bird_indep = true, compartement_indep = false,
       food_indep = true, hyp = :cch),
      (bird_indep = true, compartement_indep = false,
       food_indep = true, hyp = :fph1),
      (bird_indep = true, compartement_indep = false,
       food_indep = true, hyp = :fph2),
     ]
ms_hypothetical = [(bird_indep = true, compartement_indep = true,
                    food_indep = true),
                   (bird_indep = true, compartement_indep = false,
                    food_indep = true, constraints = ((1, >, 2),))]
ms_1 = [(bird_indep = true, compartement_indep = false, food_indep = false,
         hyp = :cch),
        (bird_indep = true, compartement_indep = false, food_indep = false,
         hyp = :fph1)]

N = 500
res1 = [[posterior(cached1, n, ms, only_data = false) for _ in 1:N] for n in [3, 6, 9]]
res2 = [[posterior(cached2, n, ms, only_data = false) for _ in 1:N] for n in [3, 6, 9]]
res11 = [[posterior2(cached11, n, ms_1) for _ in 1:N] for n in [6, 12, 18]]
res21 = [[posterior2(cached21, n, ms_1) for _ in 1:N] for n in [6, 12, 18]]
resh1 = [[posterior2([Array(datah1[i, 3:4]) for i in 1:nrow(datah1)],
                                 n, ms_hypothetical) for _ in 1:N] for n in [6, 12, 18]]
resh2 = [[posterior2([Array(datah2[i, 3:4]) for i in 1:nrow(datah1)],
                                 n, ms_hypothetical) for _ in 1:N] for n in [6, 12, 18]]

p1 = plotcachingpattern(res1)
p2 = plotcachingpattern(res2)
pgfsave(joinpath(figpath, "prediction_Amodio_exp2_1.tikz"), p1)
pgfsave(joinpath(figpath, "prediction_Amodio_exp2_2.tikz"), p2)

bayesfactors(res, models = [(1, 2), (1, 4), (2, 4), (2, 3)]) = [[r[1][i] ./ r[1][j] for r in res] for (i, j) in models]
cchconsistency(res) = [r[2][2] > r[2][1] && r[2][2] > r[2][3] &&
                       r[3][2] < r[3][1] && r[3][2] < r[3][3] for r in res]
fph1consistency(res) = [r[2][1] > r[2][2] && r[2][1] > r[2][3] &&
                        r[3][2] > r[3][1] && r[3][2] > r[3][3] for r in res]
fph2consistency(res) = [r[2][1] > r[2][2] && r[2][3] > r[2][2] &&
                        r[3][2] > r[3][1] && r[3][2] > r[3][3] for r in res]

bayesfactors2(res) = [r[1][1] / r[1][2] for r in res]
bfh1 = bayesfactors2.(resh1)
bfh2 = bayesfactors2.(resh2)

p1 = plot_bayesratioh(bfh1, legend = [], opt = @pgf({ymin = 1e-25, ymax = 1e2}))
p2 = plot_bayesratioh(bfh2, opt = @pgf({ymin = 1e-25, ymax = 1e2, legend_pos = "outer north east"}))
f1 = joinpath(figpath, "prediction_bayesfactor_hyp1.tikz")
pgfsave(f1, p1)
f2 = joinpath(figpath, "prediction_bayesfactor_hyp2.tikz")
pgfsave(f2, p2)

p1 = plotcachingpattern_hyp(resh1)
p2 = plotcachingpattern_hyp(resh2)
f1 = joinpath(figpath, "prediction_hyp_1.tikz")
pgfsave(f1, p1)
f2 = joinpath(figpath, "prediction_hyp_2.tikz")
pgfsave(f2, p2)

replnan(x) = isnan(x) ? rand() : x
replnan(x::AbstractVector) = replnan.(x)
bf11 = bayesfactors2.(res11)
bf21 = bayesfactors2.(res21)
p1 = plot_bayesratioh(replnan.(bf11), data = 0.00001/0.0014,
                      legend = [], opt = @pgf({ymin = 1e-25, ymax = 1e4}))
p2 = plot_bayesratioh(replnan.(bf21), data = 0.00001/0.0014,
                      opt = @pgf({ymin = 1e-25, ymax = 1e4, legend_pos = "outer north east"})
                     )
f1 = joinpath(figpath, "prediction_bayesfactor11.tikz")
pgfsave(f1, p1)
f2 = joinpath(figpath, "prediction_bayesfactor12.tikz")
pgfsave(f2, p2)

p1 = plotcachingpattern1(cached11)
p2 = plotcachingpattern1(cached21)
f1 = joinpath(figpath, "prediction_Amodio_exp1_1.tikz")
pgfsave(f1, p1)
f2 = joinpath(figpath, "prediction_Amodio_exp1_2.tikz")
pgfsave(f2, p2)


text:
using CuYao
using CuArrays, GPUArrays
using Test
using CUDAnative

@testset "isapprox-complex" begin
    ca = CuArray(randn(ComplexF64,3,3))
    cb = copy(ca)
    #@test ca ≈ cb still error!
    cb[1:1, 1:1] .+= 1e-7im
    @test isapprox(ca, cb, atol=1e-5)
    @test !isapprox(ca, cb, atol=1e-9)
end

@testset "view general" begin
    a = randn(5,6,8) |> CuArray
    @test view(a, 2:4, 4, [1,4,3]) |> size == (3, 3)
end

@testset "permutedims vector" begin
    ca = randn(ComplexF64,3,4,5,1)
    @test permutedims(CuArray(ca), [2,1,4,3]) |> Array ≈ permutedims(ca, [2,1,4,3])
end

@testset "Complex pow" begin
    for T in [ComplexF64, ComplexF32]
        a = CuArray(randn(T, 4, 4))
        @test Array(CUDAnative.pow.(a, Int32(3))) ≈ Array(a).^3
        @test Array(CUDAnative.pow.(a, real(T)(3))) ≈ Array(a).^3
    end
end


text:
#!/usr/bin/julia
#
# Discrete grids to run simulations on.
export Grid, Grid1D, UniformGrid;

"General abstract type for the implementation of discretised spaces."
abstract type Grid{S<:Space,N,T<:Number} <: AbstractArray{T,N} end;
Base.size(g::Grid) = size(g.grid);
Base.getindex(g::Grid{S,N,T}, I::Vararg{Int,N}) where {S,N,T} =
  get(g.grid, I, zero(T));

"Define Grid1D alias to one dimensional grids for convenience."
Grid1D{S,T} = Grid{S,1,T};
# Use linear indexing for 1D grids and directly access underlying data.
Base.IndexStyle(::Type{<:Grid1D}) = IndexLinear();
Base.getindex(G::Grid1D, I::Int) = G.grid[I][1];
Base.setindex!(G::Grid1D, v, I::Int) = (G.grid[I][1] = v);

"Grid with equally spaced grid points."
struct UniformGrid{S,N,T} <: Grid{S,N,T}
  "Maximum spatial value in each dimension (the minimum is -max)."
  max::Vector{T};

  "Grid spacing in each dimension."
  Δ::Vector{T};

  "The grid points."
  grid::Array{Vector{T},N};

  """
      UniformGrid{S,N,T}(max::Vector{T}, Δ::Vector{T}) where {S,N,T}

  Construct uniform grid of two or more dimensions.

  TODO: Finish implementing.
  """
  function UniformGrid{S,N,T}(max::Vector{J}, Δ::Vector{J}) where
            {S,N,T,J<:Number,K<:Number}
    if (length(max) != length(Δ))
      throw(DimensionMismatch((length(max),length(Δ)),
                              "must have same number of maxium values and grid spacings"));
    elseif (length(max) != N)
      throw(DimensionMismatch(length(max), "size of vectors must match dimension of grid"));
    end

    values = LinRange.(-max,Δ,max);
    #grid = zeros(T, Tuple(length.(values)))
    grid = Array{Vector{T},N}(undef, Tuple(length.(values)));

    # Scope through each dimension.
    for i ∈ CartesianIndices(grid)
      for j ∈ 1:length(values) end
    end

    #=grid = zeros(T, Tuple(2floor.(max./Δ).+1));
    for i ∈ 1:N
      grid[i,:] = -max[i]:Δ[i]:max[i];
    end=#
    new(max, Δ, LinRange.(-max,Δ,max));
  end

  """
      UniformGrid{S,1,T}(max::Number, Δ::Number) where {S,T}

  Construct one-dimensional uniform grid.
  """
  function UniformGrid{S,1,T}(max::Number, Δ::Number) where {S,T}
    values = collect(-max:Δ:max);

    len = length(values);

    # Always have an odd number of elements.
    if (len%2 == 0) push!(values, values[len]+Δ); len+=1; end

    grid = Vector{Vector{T}}(undef, len);

    for i ∈ 1:len grid[i] = Vector{T}([values[i]]); end

    new([ max ], [ Δ ], grid);
  end

  """
      UniformGrid{S,T,1}(max::Number, Δ::Number) where {S<:MomentumSpace,T}

  Construct one-dimensional uniform grid in momentum space.

  The space is not stored sequentially to allow it to be used with Fast Fourier
  transform functions (such as those implemented by the FFTW package).
  """
  function UniformGrid{S,1,T}(max::Number, Δ::Number) where {S<:MomentumSpace,T}
    k = collect(-max:Δ:max);
    len = length(k);

    # Always have an odd number of elements.
    if (len%2 == 0) push!(k, k[len]+Δ); len+=1; end

    # Compute the index of the last negative entry in k.
    lastnegidx = floor(Int,len/2);

    # Split and re-order the array around the zero/sign-change point.
    k1 = [k[i] for i in 1:lastnegidx];
    k2 = [k[i] for i in lastnegidx+1:length(k)]
    k = vcat(k2, k1);

    grid = Vector{Vector{T}}(undef, len);

    for i ∈ 1:len grid[i] = Vector{T}([k[i]]); end

    new([ max ], [ Δ ], grid);
  end
end

"""
    getproperty(g::UniformGrid{S,T,1}, s::Symbol) where {S,T}

Overwrite `getproperty()` for one-dimensional `UniformGrid` objects to allow
direct access to the single value stored within the vector.
"""
function Base.getproperty(g::UniformGrid{S,1,T}, s::Symbol) where {S,T}
  if (s === :max || s === :Δ) return Base.getfield(g, s)[1];
  else return Base.getfield(g, s); end
end

"""
    Base.:*(G::UniformGrid{S,1}, m::Real) where {S<:Space}

Multiplying a `UniformGrid` object by a real scalar returns an appropriately
scaled `UniformGrid` object.
"""
Base.:*(G::UniformGrid{S,1,T}, m::Number) where {S,T} =
  UniformGrid{S,1,T}(G.max * m, G.Δ * m);
Base.:*(m::Number, G::UniformGrid) = Base.:*(G, m);

"""
    UniformGrid{S}(xgrid::UniformGrid{C,N,T}) where
      {S<:MomentumSpace,C<:ConfigurationSpace,N,T}

Construct discrete and uniform momentum space associated with the corresponding
position/configuration space given by `xgrid`.
"""
function UniformGrid{S}(xgrid::UniformGrid{C,N}) where
            {S<:MomentumSpace,C<:ConfigurationSpace,N}
  # Use Δk = (2πħ/x) where ħ = 1.
  Δk = π/xgrid.max;
  kgrid = UniformGrid{MomentumSpace,N,typeof(Δk*1.0)}(floor(Int, length(xgrid)/2) * Δk, Δk);

  # Sanity check.
  check_same_size(xgrid, kgrid);

  return kgrid;
end

"""
    Uniform1DGrid(xmax::Number, Δx::Number; debug::Bool=false)

Construct appropriate discrete, uniform, one dimensional configuration and
momentum spaces using the parameters passed.

The `xmax` and `Δx` parameters refer to the properties of the configuration
space. An optional keyword boolean parameter `debug` may be passed to print
information about the constructed spaces.
"""
function UniformGrid(xmax::Number, Δx::Number; debug::Bool=false)

  xgrid = UniformGrid{ConfigurationSpace}(xmax, Δx);

  if debug println(xgrid, '\n', '\n'); end

  kgrid = UniformGrid{MomentumSpace}(xgrid);

  # Debugging for grids.
  if debug
    println(k, '\n', '\n');
    println("Grid size: x ", length(x), " k ", length(k), '\n');
  end

  return (xgrid, kgrid);
end

"Default to one dimensional grid with `Float64` values."
UniformGrid{S}(max::Number, Δ::Number) where S<:Space = UniformGrid{S,1,typeof(Δ*1.0)}(max,Δ);
"Implement using `Float64` type by default."
UniformGrid{S,1}(max::Number, Δ::Number) where {S<:Space,N} =
  UniformGrid{S,1,typeof(Δ*1.0)}(max,Δ);
"For convinence, specifying the grid space to convert to is optional."
UniformGrid(grid::UniformGrid{S}) where S<:ConfigurationSpace =
  UniformGrid{MomentumSpace}(grid);

"Default grid is one dimensional, uniform, configuration space grid."
Grid(max::Number, Δ::Number) = UniformGrid{ConfigurationSpace,1}(max, Δ);
Grid(grid::UniformGrid) = UniformGrid(grid);

"""
    check_same_size(g1::Grid1D, g2::Grid1D)

Check that the passed discrete, one dimensional spaces are the same size.
"""
function check_same_size(g1::Grid1D, g2::Grid1D)
  if (length(g1) != length(g2))
    throw(DimensionMismatch(string("grids are different sizes: ", length(g1), " and ", length(g2))));
  end
end