text:
function simplextest{Lib<:PolyhedraLibrary}(lib::Lib)
  A = [1 1; -1 0; 0 -1]
  b = [1, 0, 0]
  linset = IntSet([1])
  V = [0 1; 1 0]

  ine = SimpleHRepresentation(A, b, linset)
  poly1 = polyhedron(ine, lib)
  @test !isempty(poly1)
  inequality_fulltest(poly1, A, b, linset)
  generator_fulltest(poly1, V)

  ext = SimpleVRepresentation(V)
  poly2 = polyhedron(ext, lib)
  @test !isempty(poly2)
  inequality_fulltest(poly2, A, b, linset)
  generator_fulltest(poly2, V)

  # x_1 cannot be 2
  @test isempty(polyhedron(SimpleHRepresentation([A; 1 0], [b; 2], union(linset, IntSet([4]))), lib))

  # We now add the vertex (0, 0)
  V0 = [0 0]
  ext0 = SimpleVRepresentation(V0)

  push!(poly1, ext0)
  inequality_fulltest(poly1, A, b, IntSet([]))
  generator_fulltest(poly1, [V0; V])
  push!(poly2, ext0)
  inequality_fulltest(poly2, A, b, IntSet([]))
  generator_fulltest(poly2, [V0; V])

  # nonnegative orthant cut by x_1 + x_2 = 1
  Vray = [1 0; 0 1]
  extray = SimpleVRepresentation(Matrix{Int}(0,2), Vray)
  poly3 = polyhedron(extray, lib)
  Acut = [1 1]
  bcut = [1]
  linsetcut = IntSet([1])
  inecut = SimpleHRepresentation(Acut, bcut, linsetcut)
  push!(poly3, inecut)
  inequality_fulltest(poly3, A, b, linset)
  generator_fulltest(poly3, V)

  poly4 = project(poly1, [1; 0])
  inequality_fulltest(poly4, [-1; 1], [0, 1], IntSet())
  generator_fulltest(poly4, [0; 1], [])

  #\
  # \
  # |\
  # |_\
  #    \
  Alin = [1 1]
  blin = [1]
  linsetlin = IntSet(1)
  Vlin = [1 0]
  Rlin = [1 -1]
  inelin = SimpleHRepresentation([1 1; -1 -1], [1, -1], IntSet())
  plin = polyhedron(inelin, lib)
  inequality_fulltest(plin, Alin, blin, linsetlin)
  generator_fulltest(plin, Vlin, Rlin, IntSet(), IntSet(1))
  ineout = getinequalities(plin)
  @test ineout.linset == IntSet(1)
  Vlin = [1 0]
  Rlin = [1 -1]
  extlin = SimpleVRepresentation(Vlin, [1 -1; -1 1])
  plin = polyhedron(extlin, lib)
  inequality_fulltest(plin, Alin, blin, linsetlin)
  generator_fulltest(plin, Vlin, Rlin, IntSet(), IntSet(1))
  extout = SimpleVRepresentation(getgenerators(plin))
  @test extout.Rlinset == IntSet(1)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1bc961f0-66f7-11eb-17c8-2b54d4c4328d
begin
	push!(LOAD_PATH, ".")
	using CellularBase
end

# ╔═╡ f8149ae0-66f6-11eb-0335-39baafdb24c8
begin
	using Plots
	using Images
	using PlutoUI
	using Statistics
end

# ╔═╡ 9a74ef90-67ad-11eb-1fc3-e7080fbfc7e1
html"""
Term Paper 1

This work is a submission for the course Modelling Complex Systems(IDC-631) for the year 2020-21 Spring (Even) Semester by Dhruva Sambrani, MS18163
"""

# ╔═╡ eb2aad32-687f-11eb-0c6b-91b443cdb3d3
md"""# Abstract
This paper looks at the effect of a Lockdown on the spread of a pathogen. This setup is modelled via a Probabilisitic/Non Deterministic 2D cellular automaton.

This paper first [introduces](#7cb22830-67c9-11eb-0e0e-41f2316ed468) some working definitions along with the same definitions in code.
Then features are added, and parallelly introduced to the reader, [layer](#c0318a00-67ac-11eb-39b7-e3c13a5aff57) upon [layer](#1cf29850-66fa-11eb-3f76-db755c0a4c50) with logical reasoning for their implementations.
Following this, the actual topic of the paper is introduced with some [results](#5b1e6e50-6877-11eb-1ca4-3df17b85575e) that show the efficiency of [Lockdowns](#9e71a7c0-67f1-11eb-0282-672d358d2292).
Finally, [future scope of this work](#c7ea4fa0-687b-11eb-0dbc-7f316bd4aefd) is also discussed.

!!! note "Note to Evaluator"
    The author believes that papers which involve a lot of simulation and code should also be accompanied with the code they use so that readers can gain deeper insight into the processes described in the paper. Hence, the author has attempted to create a paper which interlaces the information with the corresponding code.

The coding language used is `Julia`[^julia] and following document is made with `Pluto.jl` [^pluto].

"""

# ╔═╡ 6f7bce82-67ad-11eb-1ec1-79ce2620fe9b
PlutoUI.TableOfContents(title="Term Paper 1 - Index")

# ╔═╡ 3459c040-6858-11eb-09c1-95c7bdb02722
html""

# ╔═╡ 7cb22830-67c9-11eb-0e0e-41f2316ed468
md"# Code Setup"

# ╔═╡ 65ff1d80-66f9-11eb-11c2-69c47b8458e4
md"## `Include packages`"

# ╔═╡ 330799c0-6884-11eb-230d-71c512e4a447
md"`CellularBase`[^cellularbase] is a package written by the author themselves to abstract out certain repetitive funcitions for all Cellular Automata."

# ╔═╡ db4fd9e0-6883-11eb-384e-6df2c1d7a52b
md"The `Plots.jl`[^plots] and `Images.jl`[^images] packages are used for plotting and the `Statistics`[^stats] Package for some statistical analysis. `PlutoUI.jl`[^plutoui] is used for some UI elements such as the ToC"

# ╔═╡ a95c8c70-67ad-11eb-124b-7b66d5547321
md"## Define `Structures` and `Functions`"

# ╔═╡ 73ea0ae0-66f9-11eb-19cf-b7ab6151ef63
md"### Define `LifeGrid`"

# ╔═╡ 038cbdf2-6880-11eb-1b2b-2b4c4974be49
md"""This is the basic structure that defines all the parameters of a particular "grid". The name `LifeGrid` was given by mistake but was not changed as it would entail too many other changes"""

# ╔═╡ 91afcec0-66f9-11eb-37b1-f78de5e0ba8a
md"### Define method of evolution"

# ╔═╡ fcc23790-686e-11eb-03d7-1f0946446f49
md"The method of evolution follows the following algorithm -
1. If the agent is susceptible, infect the agent depending on neighborhood and return. The exact dependacy is discussed later in [Non Deterministic runs](#2563e210-67d5-11eb-089e-f376da973765)
2. Otherwise, the state is increased by 1
3. If the state is recovered, set it to 0

Mathematically this can be written as

$s_{i+1} = \begin{cases}
    1, & s_i=0 \text{ and } \rho(\text{neighbours}) \\
    0, & s_i=0 \text{ and } 1-\rho(\text{neighbours}) \\
	(s_i+1) \mod (\tau_i+\tau_r+1) & \text{otherwise}
\end{cases}$

"

# ╔═╡ b8e7c920-66f9-11eb-397f-3760662a648c
md"### Define States"

# ╔═╡ 8b6833e0-6870-11eb-2593-69edaec0a7b8
md"States are defined in the following way
1. Susceptible $\implies 0$
2. Infected $\implies i \in [1, \tau_i]$
3. Refactory $\implies i \in [\tau_i + 1, \tau_i+\tau_r]$
4. Recovered $\implies i = (\tau_i+\tau_r + 1)$
"

# ╔═╡ de865ed2-66f9-11eb-31dc-c782d4338af3
md"### Define helper `functions`"

# ╔═╡ c3e6cc10-66f7-11eb-039f-df82addf1a5e
begin
	CellularBase.possible_states(l::LifeGrid) = 0:(l.τi + l.τr)
	CellularBase.newstate(grid::LifeGrid) = grid.otherstate;
	function CellularBase.state!(grid::LifeGrid, newstate)
	    grid.otherstate = grid.state
	    grid.state = newstate
	end
end

# ╔═╡ 4d4cdcc0-67b0-11eb-2278-8de64128d6f1
function Base.show(io::IO, ::MIME"text/html", grid::LifeGrid)
	final = Markdown.parse(
		"---" *
		"\n\n**Grid Type:** " * repr(typeof(grid), context=io) *
		"\n\n**\$\\tau_i:\$** " * string(grid.τi) *
		"\n\n**\$\\tau_r:\$** " * string(grid.τr) *
		"\n\n**Max neighbors for probabilistic spread:** " * string(grid.max) *
		"\n\n**Possible States:** " * string(possible_states(grid)) *
		"\n\n**Boundary Condition:** " * string(boundary_conditions(grid)) *
		"\n\n**Size:** " * string(size(grid)) *
		"\n\n**Neighborhood:** "
	)
	show(io, MIME("text/html"), final)
	show(io, MIME("text/html"), grid.neighborhood)
	show(io, MIME("text/html"), md"---")
end

# ╔═╡ 7cdebc60-684b-11eb-38a2-45873e4b51d0
function makegraph(ldgrid, ldresults, ldradii)
    p = plot(1:length(ldresults), ldresults .|> frame -> count(frame .|> p -> susceptible(p, ldgrid)) / length(state(ldgrid)), color=:green, label="Sus")
    plot!(1:length(ldresults), ldresults .|> frame -> count(frame .|> p -> infectious(p, ldgrid)) / length(state(ldgrid)), color=:red, label="Inf")
    plot!(1:length(ldresults), ldresults .|> frame -> count(frame .|> p -> refactory(p, ldgrid)) / length(state(ldgrid)), color=:blue, label="Ref")
    q = plot(ldradii, label="Neighbourhood Radius", legend=:bottomright)
    plot(p, q, size=(600, 500), layout=(2, 1))
end;

# ╔═╡ 5031733e-66f9-11eb-107f-6d69089d9384
function color(i, grid)
    if susceptible(i, grid)
        RGB(0.84, 0.95, 0.80)
    elseif infectious(i, grid)
        RGB(0.83, 0.11, 0.31)
    else
        RGB(0.16, 0.11, 0.42)
    end
end;

# ╔═╡ c0318a00-67ac-11eb-39b7-e3c13a5aff57
md"# Deterministic Case - A First Run"

# ╔═╡ 47b81340-67cb-11eb-372e-972e4d4ff746
md"## Introduction
As a first trial, we can try to run the simulation for the deterministic case, where the infection can spread even if 1 neighbor is infected. We also verify that the results observed are very similar to results obtained in [^Moitra2019]
"

# ╔═╡ 02c24a70-67ae-11eb-2944-61a7963eb947
LifeGrid((100, 100), 4, 9, max=1, nh=VonNeumannNeighborhood(1, 2))

# ╔═╡ fbdba862-67cf-11eb-39bc-f37eb3662785
md"## `state(center patient) = 1`"

# ╔═╡ fe0dda70-66fa-11eb-0e02-73c05666aace
dgrid, dresults = let
	grid = LifeGrid(
		(100, 100),
		4, 9,
		max=1,
		nh=VonNeumannNeighborhood(1, 2));
	results = simulate!(grid, 150)
	grid, results
end;

# ╔═╡ 8315d810-677a-11eb-3953-8faad42ada40
@gif for frame ∈ dresults[5:end]
	plot(frame .|> p -> color(p, dgrid))
end

# ╔═╡ 7d24e300-67d0-11eb-0864-41c6166f13ac
md"## Random states from 0 to $\tau_i$"

# ╔═╡ b8523010-67ce-11eb-2396-4780284dfa01
randgrid, randresults = let
	arr = rand(0:4, 100, 100)
	grid = LifeGrid(
		arr,
		4, 5,
		max=1,
		nh=VonNeumannNeighborhood(1, 2));
	results = simulate!(grid, 50)
	grid, results
end;

# ╔═╡ cf22a450-67ce-11eb-2141-fd6f3efc689a
@gif for frame ∈ randresults
	plot(frame .|> p -> color(p, randgrid))
end

# ╔═╡ 07c440a0-67d1-11eb-3a8b-d559cac6276e
md"## `state(center agent) = τᵢ; τᵢ=τᵣ`"

# ╔═╡ 30ce7a10-67bd-11eb-027d-8dc583201bc1
ϕgrid, ϕresults = let
	arr = zeros(Int, 100, 100)
	arr[50,50] = 4
	grid = LifeGrid(
		arr,
		4, 4,
		max=1,
		nh=VonNeumannNeighborhood(1, 2));
	results = simulate!(grid, 50)
	grid, results
end;

# ╔═╡ 4ca7fe50-67bd-11eb-382c-1304f3c1c3c4
@gif for frame ∈ ϕresults
	plot(frame .|> p -> color(p, ϕgrid))
end

# ╔═╡ 643ffb90-67ad-11eb-20c6-7744050e5448
md"""## Observations in the deterministic case

We can see that the infection spreads outward and then dies off.

The refactory agents act as a buffer between the leading edge of infected agents and the recovered(susceptible) agents inside the wave. This is similar to how a batch of burnt material doesn't allow the fire to move back in.

!!! tip "Spread Back"
    The infection can spread back into the wave only if either -
    -
      1. the radius of the neighborhood is larger than the refactory period
      2. if $\tau_i \ge \tau_r$ and $\text{state(patient zero)}\ge\tau_r$

"""

# ╔═╡ 1cf29850-66fa-11eb-3f76-db755c0a4c50
md"# Non Deterministic Case"

# ╔═╡ 393a5bc0-67cb-11eb-3be1-a9ec7943b5f3
md"""## Why?
Mainly because determinism is so 1900s.

Seriously though, the reason why we don't see too many interesting things in the deterministic case is because of the refactory buffer separating the infected leading edge and the recovered members.

Allowing some agents to escape the infection wave allows breaks in the refactory buffer allowing the infection to come back into the recovered agents. It seems very sadistic to play for the pathogen, but what's the fun otherwise.
"""

# ╔═╡ 2563e210-67d5-11eb-089e-f376da973765
md"""## How is it implemented?
In code, obviously. Oh you mean what's the logic behind it?

If the spread is probabilistic, one way to do it would be to set a non zero probability if the number of infected neighbors is non-negative. While this does make the spread probabilistic, it doesn't model _real life_ very well. In real life, one's chances of infection increases with the increase in infected neighbors. While the exact distribution is unknown, we can model the spread by assuming that the chance of spreading follows

$\rho_\text{infect}{}(\text{count}) = \begin{cases}
  \frac{\text{count}}{\text{max}} & 0\leq \text{count}\leq \text{max} \\
  1 & \text{max}\leq \text{count}
\end{cases}$

where $\text{max}$ is a parameter of the simulation. This basically means that for every infected neighbour, the chance of infection is linear, until a limit, after which the agent will always get infected.

This definition has the added benefit of being independent on the definition of the neighbourhood, and can be attributed only to the pathogen in question.

Every `LifeGrid` has a critical limit of infected neighbors, denoted by `max`. For every agent, $\rho$ is calculated and infection is done by a Bernoulli trial.

!!! info "High values of max"

    Hence, it also makes sense to set `max` higher than the size of the neighbourhood itself, for that would allow us to set the highest possible probability to be lesser than 1.

    After all, what if the neighborhood changes during the run? How? Keep reading
"""

# ╔═╡ dafac300-67e9-11eb-0972-f9070193399e
md"## Max = 4"

# ╔═╡ 46055e90-67e9-11eb-2b6f-dffa58138005
ρgrid, ρresults = let
	grid = LifeGrid(
		(100, 100),
		4, 5,
		max=4,
		nh=VonNeumannNeighborhood(1, 2));
	results = simulate!(grid, 400)
	grid, results
end;

# ╔═╡ 5626db00-67e9-11eb-0e12-a77b6b9c58be
@gif for frame ∈ ρresults
	plot(frame .|> p -> color(p, ρgrid))
end every 2

# ╔═╡ d4c2cff0-67e9-11eb-1f5b-03c7bab91782
md"## Max = 9"

# ╔═╡ 3056cdae-67ec-11eb-0b35-b3d6cde9b0fa
ρ7grid, ρ7results = let
	grid = LifeGrid(
		(100, 100),
		4, 5,
		max=9,
		init=5,
		nh=VonNeumannNeighborhood(1, 2));
	results = simulate!(grid, 100)
	grid, results
end;

# ╔═╡ 3b0024a0-67ec-11eb-3398-eb07250758a9
@gif for frame ∈ ρ7results
	plot(frame .|> p -> color(p, ρ7grid))
end

# ╔═╡ 6d822d52-67ed-11eb-1b71-679d491bf133
md"## Max = 2"

# ╔═╡ 76e43730-67ed-11eb-201f-53bf8d3e3bfb
ρ2grid, ρ2results = let
	grid = LifeGrid(
		(100, 100),
		4, 5,
		max=2,
		nh=VonNeumannNeighborhood(1, 2));
	results = simulate!(grid, 150)
	grid, results
end;

# ╔═╡ 8445d050-67ed-11eb-2411-f776a43c7f40
@gif for frame ∈ ρ2results
	plot(frame .|> p -> color(p, ρ2grid))
end every 2

# ╔═╡ 0f5ba100-67ef-11eb-15f5-3fbaad70dde2
md"""## Observations

Evidently, higher max values imply smaller probability of spreading. This is evident in the run with `max = 9`, where very few infected cases are seen at any point in time, and agents often recover before they spread it on, leading to the pathogen dying out.

On the flip side, if `max` is too close to 1, then there it acts very similar to the deterministic case _most of the time_. This is also evident in the fact that the case where `max` is 1, it is exactly the deterministic case which we've seen also dies out.

!!! tip "Relevant Quote"
    It is possible to have too much of a good thing

    \- **Aesop**
"""

# ╔═╡ 9e71a7c0-67f1-11eb-0282-672d358d2292
md"""# Lockdown
## Introduction
The word _lockdown_ seems to have taken over our lives and completely changed the way that we go about our daily lives. But why are we doing this? And how effective are these lockdowns? We explore this question by looking at how we can simulate a lockdown period for our agents and see how this affects the infection spread.
"""

# ╔═╡ ef666b40-67f4-11eb-2f31-05145a6dfcba
md"""## The logic behind lockdowns

!!! info "Definition"
    A lockdown is a restriction policy for people or community to stay where they are, usually due to specific risks to themselves or to others if they can move and interact freely - [Wikipedia](https://en.wikipedia.org/wiki/Lockdown)

From the definition, it is clear that the intent is to prevent people from interacting and hence hopefully reducing the spread of the disease. If an agent interacts with fewer people, they have a lower chance of getting infected.
"""

# ╔═╡ 5586a260-67f4-11eb-11e2-b706fd0d4dfe
md"""## Neighborhoods, Virulence and Infection Rate
Neighbourhood size and intrinsic virulence of a pathogen together decide the pathogen's spreading capabilities.

Neighbourhood size is a representation of the amount of inter-agent interactions. Due to the way the model is defined, it also shows how _far_ the agent can influence. It defines the number of agents one agent can meet.

Intrinsic virulence of a pathogen is how well it spreads given that two agents meet. There is sufficient discussion regarding this in the above section.

If both these parameters define, in some sense, the same thing, why not just club them together into a single parameter which gives the net spread?

The reason is 3 fold-

1. The **neighborhood is a function of the agents' behavior**, the **virulence is a function of the pathogen's behaviour**. It makes sense to keep them separate to deal with different cases of spread changes.
2. The model is such that there are **geometric differences** that arise by changing the neighborhood radius. The _buffering_ action of the refactory line is reduced by a larger neighborhood. It is not possible to increase the infection rate without also reducing this phenomenon if both were clubbed together.
3. The code would need too many changes.

!!! tips "Relevant quote"
    Any sufficiently advanced bug is indistinguishable from a feature.

    \- **Rich Kulawiec**

"""

# ╔═╡ 21b80dc0-67f9-11eb-0e87-6db59c88c296
md"""## Lockdowns and Neighborhoods

From the above discussion, it must be clear that a lockdown can be simulated pretty easily by varying the size of the neighborhood. The lockdown must also start at some predetermined point during the simulation. The size of the neighbourhood is changed by changing `grid.neighborhood` and the in-simulation change of the neighborhood can be executed by using the `postrunhook` of the `simulate!` function

"""

# ╔═╡ ead37ab2-66f9-11eb-2e56-07f105753991
function lockdown(grid, step; step_dict, radii)
    push!(radii, radius(grid.neighborhood))
    if step in keys(step_dict)
        grid.neighborhood = VonNeumannNeighborhood(step_dict[step], 2);
    end
end

# ╔═╡ 5641a270-684b-11eb-29b4-ef8c59ef9679
md"## 200 → 1"

# ╔═╡ 0f20b710-6840-11eb-1712-9f4614e9ed87
ld1 = let
	step_dict = Dict(200 => 1) # Step -> radius
	radii = []
	grid = LifeGrid(
		(150, 150),
		4, 5,
		max=9,
		init=5,
		nh=VonNeumannNeighborhood(3, 2));
	results = simulate!(grid, 400;
		postrunhook=lockdown, step_dict=step_dict, radii=radii);
	grid, results, radii
end;

# ╔═╡ 6f9dc770-6847-11eb-14ff-2b2e8562076f
@gif for frame ∈ ld1[2]
	plot(frame .|> p -> color(p, ld1[1]))
end every 2

# ╔═╡ 2bf1c112-6848-11eb-3b4a-a3ce569da696
makegraph(ld1...)

# ╔═╡ a9588a00-684b-11eb-2585-09a2cdf9ad80
md"## 50 → 2; 150 → 1; 300 → 3"

# ╔═╡ 4f4fdfd0-684c-11eb-3769-49d1677a1d52
ld2 = let
	step_dict = Dict(50 => 2, 150 => 1, 300 => 3) # Step -> radius
	radii = []
	grid = LifeGrid(
		(150, 150),
		4, 5,
		max=9,
		init=5,
		nh=VonNeumannNeighborhood(3, 2));
	results = simulate!(grid, 400;
		postrunhook=lockdown, step_dict=step_dict, radii=radii);
	grid, results, radii
end;

# ╔═╡ d97dddfe-684c-11eb-1d11-5734830507a2
@gif for frame ∈ ld2[2]
	plot(frame .|> p -> color(p, ld2[1]))
end every 2

# ╔═╡ e1baa1c0-684c-11eb-292c-3336d6d715b1
makegraph(ld2...)

# ╔═╡ 52123eb0-684d-11eb-0179-895599ae7a2f
md"## 100 → 2; 300 → 3"

# ╔═╡ 4e182db0-684d-11eb-0bab-b7c22d3f3ff5
ld3 = let
	step_dict = Dict(100 => 2, 300 => 3) # Step -> radius
	radii = []
	grid = LifeGrid(
		(150, 150),
		4, 5,
		max=9,
		init=5,
		nh=VonNeumannNeighborhood(3, 2));
	results = simulate!(grid, 400;
		postrunhook=lockdown, step_dict=step_dict, radii=radii);
	grid, results, radii
end;

# ╔═╡ 53bb3e50-684e-11eb-1bd3-d582200be8b5
@gif for frame ∈ ld3[2]
	plot(frame .|> p -> color(p, ld3[1]))
end every 2

# ╔═╡ 660d5660-684e-11eb-27ad-cbde3f616410
makegraph(ld3...)

# ╔═╡ f0c45ab2-684e-11eb-24ca-c5942b62f396
md"""## Observations

Apart from being very flashy, we see an interesting phenomenon. When the radius of the neighbourhood is changed, it takes about 25-30 steps for the grid to find a new equilibrium. During this transient phase, the distribution of states can vary wildly.

Further, the equilibrium fractions of the states averaged over a few cycles is nearly equal, but the exact variation may vary on the basis of initial conditions ([See 100 → 2; 300 → 3](#52123eb0-684d-11eb-0179-895599ae7a2f)). This can be attributed to the geometry of the neighborhood, the refactory buffer and the distribution of infected agents.

Because most large scale structures don't last for long times, we can assume that the final state distribution is memoryless or has only short-term memory. Thus our study of the effectiveness of a lockdown can be reduced to a simpler study of equilibrium infection as a function of neighborhood radius and [`grid.max`](#2563e210-67d5-11eb-089e-f376da973765), which, as you probably guessed it, is our next section.
"""

# ╔═╡ 61e28250-6857-11eb-39f9-2ff5565519ab
md"""
# Equilibrium infection distribution
**As a function of Neighbourhood Radius and Virulence of pathogen**
"""

# ╔═╡ 2734ec90-6859-11eb-2430-e3bd97302fd8
md"""## Introduction
This section is mainly computational. With the functionalities we have built up in the previous sections, we simply simulate a random 200X200 grid for 200 steps for different values of neighbourhood radius and virulence, find the average(over 5 cycle-lenghts) infected fraction then plot a scatter plot to see the distribution.
"""

# ╔═╡ f7bf5830-6865-11eb-29d7-2be3d482ab20
function ss(grid, step; kwargs...)
	if count(i -> infectious(i, grid), state(grid)) == 0
		return :shortcurcuit
	end
end

# ╔═╡ 2e5c48b0-6859-11eb-2b69-0b78043bfc13
function equilibrium_infection(nhradius, p_max)
	arr = rand(1:9, 100, 100)
	grid = LifeGrid(arr, 4, 5, max=p_max, nh=VonNeumannNeighborhood(nhradius, 2))
	results = simulate!(grid, 200, postrunhook=ss)
	mean(count.(p -> infectious(p, grid), results[end - 45:end])) / 10000
end

# ╔═╡ 16c3dea0-685b-11eb-1b89-1da8ccf00941
md" ### Verification
We verify that the function is working correctly by checking for the case `nhradius = 3, p_max = 3` against the run [100 → 2; 300 → 3](#52123eb0-684d-11eb-0179-895599ae7a2f)"

# ╔═╡ ba6cd7a0-685b-11eb-3375-1b00d96dc710
begin
	res = equilibrium_infection(3, 9)
	if 0.35 < res < 0.4
		Markdown.parse("""
!!! tips "Success!"
    \$0.35< \\texttt{equilibrium\\_infection(3, 9)} = $(round(res, sigdigits=2)) <0.4\$
""")
	else
		Markdown.parse(md"""
!!! warning "Oh No!"
    \$0.35 \\nless \\texttt{equilibrium\\_infection(3, 9)} = $(round(res, sigdigits=2)) \\nless 0.4\$
    **This should never HAPPEN!!!**
""")
	end
end

# ╔═╡ 977a4062-6861-11eb-17eb-03874dfb3c79
md"""
!!! danger "Long running Code!"
    The following cell takes a long time (≈ 400 secs) to run. Be careful!
    You may want to reduce `nhs` and `maxs`, or skip some elements.
"""

# ╔═╡ 5b1e6e50-6877-11eb-1ca4-3df17b85575e
md"## Results"

# ╔═╡ 104a9570-685e-11eb-06c2-f9df4e63d1ba
begin
	nhs = 1:4
	maxs = 1:2:25
	points = [equilibrium_infection(nhradius, p_max) for nhradius in nhs, p_max in maxs]
end;

# ╔═╡ c6754200-685e-11eb-319f-f1d06650a620

# ╔═╡ 3f6db120-685e-11eb-3fa6-fd5ec5dfe981
let
	plotlyjs()
	p = plot(maxs, nhs, points, st=:surface, ylabel="NH Radius", xlabel="Virulence(max)", zlabel="Frac infected", title="Surface Plot")
	gr()
	p
end

# ╔═╡ 8af3d922-686e-11eb-03bb-d96c3ac24780
plot(maxs, nhs, points, ylabel="NH Radius", xlabel="Virulence(max)", title="Contour Plot", width=3)

# ╔═╡ 97bff6c2-6873-11eb-2f3d-c108ce4accde
begin
	plotlyjs()
	p = plot(maxs, nhs, points, ylabel="NH Radius", xlabel="Virulence(max)", title="HeatMap", st=:heatmap, hover = string.(reshape(points, :)))
	gr()
	p
end

# ╔═╡ 44881ce0-6877-11eb-0cef-f35b60baf357
md"
## Observations
1. As expected, runs with low virulence and low radius (bottom-right) have very low equilibrium infection levels.
2. For very virulent pathogens, even very strong lockdowns do not have much of an effect.
3. Even for mid-level virulence pathogens, weak lockdowns are only able to reduce the infection probability, but are unable to completely kill off the pathogen.
4. As we have seen before, even very low initial viral load can quickly explode in high interaction neighborhoods, so when the lockdown is removed, the infection quickly increases to pre-lockdown levels

## Discussion on the effectiveness of lockdowns

Given that lockdowns don't help in killing off the pathogen, was the whole last year a waste? Was it all just a government hoax to install 5G towers?

While the author can't say much about government plans to install 5G towers, it is clear that lockdowns have a dampening effect in the number of active cases. This decrease in cases allow the society to focus on combating the disease by spending more efforts into research for vaccines, treatments and preventive steps, along with building infrastructure and materials which are necessary to combat the recurrence of the pathogen. It also gives epidemologists and people working in Complex Systems time to make pretty plots and animations.
"

# ╔═╡ c7ea4fa0-687b-11eb-0dbc-7f316bd4aefd
md"# Further scope of this work
While this work has largely focused on the case of pathogen spread, the same work can be extended within the same topic and to other related topics.

## Things within the same topic
1. Fitting a mathematical model of the effect of lockdown and pathogen spread to data obtained by this work, then fitting it for other real life pathogens
2. Verification of model reflecting real pathogens via case studies
3. More extended statistical analysis to quantify equilibrium
4. Search for multi-equilibrium states
5. Contact tracing and quarantine simulations and comparing effectiveness of blanket lockdowns with more surgical methods.

## Other topics which can gain from this work
While we have called agents as people and the spread as a pathogen, any system which relaxes after some time can be modelled with a similar setup. Accordingly the meaning of lockdown can change.
Some examples include -
1. Spread of radioactivity with lockdown as shielding.
2. Complex web of processors which heat up after a while. There, the lockdown would be the tradeoff between a less complex web but a less efficient system of workload balancing.
3. Spread of algae and lockdown being the speed of flow of water.
"

# ╔═╡ 01866aa0-6854-11eb-3369-7deb6c18c3fe
md"""# References
[^Moitra2019]: Moitra, P., Sinha, S. Localized spatial distributions of disease phases yield long-term persistence of infection. Sci Rep 9, 20309 (2019). [https://doi.org/10.1038/s41598-019-56616-3](https://doi.org/10.1038/s41598-019-56616-3)
[^julia]: Julia: A Fresh Approach to Numerical Computing. Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah. (2017) SIAM Review, 59: 65–98. doi: [10.1137/141000671](https://dx.doi.org/10.1137/141000671)
[^plots]: Plots - powerful convenience for visualization in Julia. Thomas Breloff. [docs.juliaplots.org](docs.juliaplots.org)
[^images]: JuliaImages: image processing and machine vision for Julia. [https://juliaimages.org/stable/](https://juliaimages.org/stable/)
[^pluto]: Pluto - Reactive Notebooks for Julia. Fons van der Plas and Mikołaj Bochenski. [https://github.com/fonsp/Pluto.jl](https://github.com/fonsp/Pluto.jl)
[^plutoui]: PlutoUI. Fons van der Plas. [https://github.com/fonsp/PlutoUI.jl](https://github.com/fonsp/PlutoUI.jl)
[^stats]: Statistics - Julia stdlib for Statistics. [https://docs.julialang.org/en/v1/stdlib/Statistics/](https://docs.julialang.org/en/v1/stdlib/Statistics/)
[^cellularbase]: CellularBase.jl: Abstract definitions for Cellular Automaton. Dhruva Sambrani. [https://dhruvasambrani.github.io/idc621/CellularBase](https://dhruvasambrani.github.io/idc621/CellularBase)
"""

# ╔═╡ 263db440-66f9-11eb-03b5-1fce784fad2d
function (_g::LifeGrid)(neighbors::Array{Int}; kwargs...)
	current = neighbors[1]
	if susceptible(current, _g)
		if infect(count(i -> infectious(i, _g), neighbors), _g)
			return 1
		else
			return 0
		end
	else
		return (current + 1) % (_g.τi + _g.τr)
	end
end

# ╔═╡ 35d0be40-66f7-11eb-2f60-01982265ccb0
mutable struct LifeGrid <: AbstractGrid{Int}
	state::Array{Int}
	otherstate::Array{Int}
	bc::BoundaryCondition
	neighborhood::AbstractNeighborhood
	τi::Int64
	τr::Int64
	max::Int64
	function LifeGrid(
			size::Tuple{Int,Int}, τi::Int, τr::Int;
			max::Int64, init=nothing,
			nh::AbstractNeighborhood=VonNeumannNeighborhood(1, 2))
		z = zeros(size)
		if isnothing(init) z[(size .÷ 2)...] = 1
		else z[rand(CartesianIndices(z), init)] .= 1
		end
		new(z, copy(z), FixedMax, nh, τi, τr, max)
	end
	function LifeGrid(
			arr::Array{Int}, τi::Int, τr::Int;
			max::Int64,
			nh::AbstractNeighborhood=VonNeumannNeighborhood(1, 2)
		)
		new(arr, copy(arr), FixedMax, nh, τi, τr, max)
	end
end

# ╔═╡ 21d1eb20-66f8-11eb-26d3-4fb244ca8646
begin
	infect(_count::Int64, g::LifeGrid) = if (_count == 0) false else rand() < 1 / g.max * _count end
	susceptible(i::Int, g::LifeGrid) = i == 0
	infectious(i::Int, g::LifeGrid) = i ∈ 1:g.τi
	refactory(i::Int, g::LifeGrid) = i ∈ (g.τi + 1):(g.τr + g.τi)
	recovered(i::Int, g::LifeGrid) = i == (g.τr + g.τi) + 1
end

# ╔═╡ Cell order:
# ╟─9a74ef90-67ad-11eb-1fc3-e7080fbfc7e1
# ╟─eb2aad32-687f-11eb-0c6b-91b443cdb3d3
# ╟─6f7bce82-67ad-11eb-1ec1-79ce2620fe9b
# ╟─3459c040-6858-11eb-09c1-95c7bdb02722
# ╟─7cb22830-67c9-11eb-0e0e-41f2316ed468
# ╟─65ff1d80-66f9-11eb-11c2-69c47b8458e4
# ╟─330799c0-6884-11eb-230d-71c512e4a447
# ╠═1bc961f0-66f7-11eb-17c8-2b54d4c4328d
# ╟─db4fd9e0-6883-11eb-384e-6df2c1d7a52b
# ╠═f8149ae0-66f6-11eb-0335-39baafdb24c8
# ╟─a95c8c70-67ad-11eb-124b-7b66d5547321
# ╟─73ea0ae0-66f9-11eb-19cf-b7ab6151ef63
# ╟─038cbdf2-6880-11eb-1b2b-2b4c4974be49
# ╠═35d0be40-66f7-11eb-2f60-01982265ccb0
# ╟─91afcec0-66f9-11eb-37b1-f78de5e0ba8a
# ╟─fcc23790-686e-11eb-03d7-1f0946446f49
# ╠═263db440-66f9-11eb-03b5-1fce784fad2d
# ╟─b8e7c920-66f9-11eb-397f-3760662a648c
# ╟─8b6833e0-6870-11eb-2593-69edaec0a7b8
# ╠═21d1eb20-66f8-11eb-26d3-4fb244ca8646
# ╟─de865ed2-66f9-11eb-31dc-c782d4338af3
# ╠═c3e6cc10-66f7-11eb-039f-df82addf1a5e
# ╠═4d4cdcc0-67b0-11eb-2278-8de64128d6f1
# ╠═7cdebc60-684b-11eb-38a2-45873e4b51d0
# ╠═5031733e-66f9-11eb-107f-6d69089d9384
# ╟─c0318a00-67ac-11eb-39b7-e3c13a5aff57
# ╟─47b81340-67cb-11eb-372e-972e4d4ff746
# ╠═02c24a70-67ae-11eb-2944-61a7963eb947
# ╟─fbdba862-67cf-11eb-39bc-f37eb3662785
# ╠═fe0dda70-66fa-11eb-0e02-73c05666aace
# ╠═8315d810-677a-11eb-3953-8faad42ada40
# ╟─7d24e300-67d0-11eb-0864-41c6166f13ac
# ╠═b8523010-67ce-11eb-2396-4780284dfa01
# ╠═cf22a450-67ce-11eb-2141-fd6f3efc689a
# ╟─07c440a0-67d1-11eb-3a8b-d559cac6276e
# ╠═30ce7a10-67bd-11eb-027d-8dc583201bc1
# ╠═4ca7fe50-67bd-11eb-382c-1304f3c1c3c4
# ╟─643ffb90-67ad-11eb-20c6-7744050e5448
# ╟─1cf29850-66fa-11eb-3f76-db755c0a4c50
# ╟─393a5bc0-67cb-11eb-3be1-a9ec7943b5f3
# ╟─2563e210-67d5-11eb-089e-f376da973765
# ╟─dafac300-67e9-11eb-0972-f9070193399e
# ╠═46055e90-67e9-11eb-2b6f-dffa58138005
# ╠═5626db00-67e9-11eb-0e12-a77b6b9c58be
# ╟─d4c2cff0-67e9-11eb-1f5b-03c7bab91782
# ╠═3056cdae-67ec-11eb-0b35-b3d6cde9b0fa
# ╠═3b0024a0-67ec-11eb-3398-eb07250758a9
# ╟─6d822d52-67ed-11eb-1b71-679d491bf133
# ╠═76e43730-67ed-11eb-201f-53bf8d3e3bfb
# ╠═8445d050-67ed-11eb-2411-f776a43c7f40
# ╟─0f5ba100-67ef-11eb-15f5-3fbaad70dde2
# ╟─9e71a7c0-67f1-11eb-0282-672d358d2292
# ╟─ef666b40-67f4-11eb-2f31-05145a6dfcba
# ╟─5586a260-67f4-11eb-11e2-b706fd0d4dfe
# ╟─21b80dc0-67f9-11eb-0e87-6db59c88c296
# ╠═ead37ab2-66f9-11eb-2e56-07f105753991
# ╟─5641a270-684b-11eb-29b4-ef8c59ef9679
# ╠═0f20b710-6840-11eb-1712-9f4614e9ed87
# ╠═6f9dc770-6847-11eb-14ff-2b2e8562076f
# ╠═2bf1c112-6848-11eb-3b4a-a3ce569da696
# ╟─a9588a00-684b-11eb-2585-09a2cdf9ad80
# ╠═4f4fdfd0-684c-11eb-3769-49d1677a1d52
# ╠═d97dddfe-684c-11eb-1d11-5734830507a2
# ╠═e1baa1c0-684c-11eb-292c-3336d6d715b1
# ╟─52123eb0-684d-11eb-0179-895599ae7a2f
# ╠═4e182db0-684d-11eb-0bab-b7c22d3f3ff5
# ╠═53bb3e50-684e-11eb-1bd3-d582200be8b5
# ╠═660d5660-684e-11eb-27ad-cbde3f616410
# ╟─f0c45ab2-684e-11eb-24ca-c5942b62f396
# ╟─61e28250-6857-11eb-39f9-2ff5565519ab
# ╟─2734ec90-6859-11eb-2430-e3bd97302fd8
# ╠═f7bf5830-6865-11eb-29d7-2be3d482ab20
# ╠═2e5c48b0-6859-11eb-2b69-0b78043bfc13
# ╟─16c3dea0-685b-11eb-1b89-1da8ccf00941
# ╟─ba6cd7a0-685b-11eb-3375-1b00d96dc710
# ╟─977a4062-6861-11eb-17eb-03874dfb3c79
# ╟─5b1e6e50-6877-11eb-1ca4-3df17b85575e
# ╠═104a9570-685e-11eb-06c2-f9df4e63d1ba
# ╟─c6754200-685e-11eb-319f-f1d06650a620
# ╠═3f6db120-685e-11eb-3fa6-fd5ec5dfe981
# ╠═8af3d922-686e-11eb-03bb-d96c3ac24780
# ╠═97bff6c2-6873-11eb-2f3d-c108ce4accde
# ╟─44881ce0-6877-11eb-0cef-f35b60baf357
# ╟─c7ea4fa0-687b-11eb-0dbc-7f316bd4aefd
# ╟─01866aa0-6854-11eb-3369-7deb6c18c3fe

ext: jl


text:
using Test
using NeXLCore

@testset "Standardize" begin
    @testset "Standardize KRatio" begin
        fe1 = KRatio(
            characteristic(n"Fe", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.6517, 0.02),
        )
        @test isstandard(fe1)
        fe2 = KRatio(
            characteristic(n"Fe", ltransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.3217, 0.02),
        )
        @test isstandard(fe2)
        ca1 = KRatio(
            characteristic(n"Ca", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Ca5(PO4)3F"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"CaF2",
            uv(0.7280, 0.02),
        )
        @test isstandard(ca1)
        ca2 = KRatio(
            characteristic(n"Ca", kbeta),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Ca5(PO4)3F"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"CaF2",
            uv(0.7326, 0.02),
        )
        @test isstandard(ca2)
        si1 = KRatio(
            characteristic(n"Si", ktransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"SiO2"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Si",
            uv(0.3801, 0.02),
        )
        @test isstandard(si1)
        @test NeXLCore.matches(KRatio(
            characteristic(n"Fe", ltransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.6517, 0.02),
            ), fe2)
        @test !NeXLCore.matches(KRatio(
            characteristic(n"Fe", ltransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.6517, 0.02),
            ), fe1)
        @test !NeXLCore.matches(KRatio(
            characteristic(n"Fe", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.45, 0.02),
        ), ca1)
        stdized = standardize([ KRatio(
                characteristic(n"Fe", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"Fe",
                uv(0.45, 0.02)),
            KRatio(
                characteristic(n"Ca", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"CaF2",
                uv(0.32, 0.02)),
            KRatio(
                characteristic(n"F", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"LiF",
                uv(0.333, 0.02)) ],
            [ fe1, fe2, ca1, ca2, si1 ])
        @test value(stdized[1]) ==0.45/0.6517
        @test isequal(stdized[1].standard, mat"Fe2O3")
        @test value(stdized[2])==0.32/0.7280
        @test isequal(stdized[2].standard, mat"Ca5(PO4)3F")
        @test isequal(stdized[3].standard, mat"LiF")
        @test value(stdized[3])==0.333
    end
    @testset "Standardize KRatios" begin
        fe1 = KRatio(
            characteristic(n"Fe", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.6517, 0.02),
        )
        fe2 = KRatio(
            characteristic(n"Fe", ltransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.3217, 0.02),
        )
        ca1 = KRatio(
            characteristic(n"Ca", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Ca5(PO4)3F"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"CaF2",
            uv(0.7280, 0.02),
        )
        ca2 = KRatio(
            characteristic(n"Ca", kbeta),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Ca5(PO4)3F"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"CaF2",
            uv(0.7326, 0.02),
        )
        si1 = KRatio(
            characteristic(n"Si", ktransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"SiO2"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Si",
            uv(0.3801, 0.02),
        )
        @test NeXLCore.matches(KRatios(
            characteristic(n"Fe", ltransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            [ uv(0.6517, 0.02), uv(0.6517, 0.02) ]
            ), fe2)
        @test !NeXLCore.matches(KRatios(
            characteristic(n"Fe", ltransitions),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            [ uv(0.6517, 0.02), uv(0.6517, 0.02) ]
            ), fe1)
        @test !NeXLCore.matches(KRatios(
            characteristic(n"Fe", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            [ uv(0.45, 0.02), uv(0.45, 0.02) ]
        ), ca1)
        stdized = standardize([ KRatios(
                characteristic(n"Fe", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"Fe",
                [ uv(0.45, 0.02), uv(0.49, 0.03)]),
            KRatios(
                characteristic(n"Ca", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"CaF2",
                [ uv(0.32, 0.02), uv(0.32, 0.02)]),
            KRatios(
                characteristic(n"F", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"LiF",
                [ uv(0.333, 0.02), uv(0.333, 0.02) ]) ],
            [ fe1, fe2, ca1, ca2, si1 ])
        @test value(stdized[1].kratios[1]) == 0.45/0.6517
        @test isapprox(σ(stdized[1].kratios[1]),0.037,atol=0.001)
        @test value(stdized[1].kratios[2]) == 0.49/0.6517
        @test isapprox(σ(stdized[1].kratios[2]),0.051,atol=0.001)
        @test isequal(stdized[1].standard, mat"Fe2O3")
        @test value(stdized[2].kratios[1])==0.32/0.7280
        @test isequal(stdized[2].standard, mat"Ca5(PO4)3F")
        @test isequal(stdized[3].standard, mat"LiF")
        @test value(stdized[3].kratios[1])==0.333
    end
    @testset "Standardize KRatios - mixed" begin
        fe1 = KRatio(
            characteristic(n"Fe", kalpha),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0), :Composition => mat"Fe2O3"),
            Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
            mat"Fe",
            uv(0.6517, 0.02),
        )
        stdized = standardize([
            KRatios(
                characteristic(n"Fe", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"Fe",
                [ 0.45, 0.49 ]),  # Float64
            KRatios(
                characteristic(n"Fe", kalpha),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                Dict(:BeamEnergy => 20.0e3, :TakeOffAngle => deg2rad(40.0)),
                mat"Fe",
                [ uv(0.45, 0.02), uv(0.49, 0.03) ]) ],  # UncertainValue
            [ fe1 ])
        # Check mixed UncertainValue and Float64
        @test stdized[1].kratios[1] == 0.45/0.6517
        @test σ(stdized[1].kratios[1]) == 0.0
        @test stdized[1].kratios[2] == 0.49/0.6517
        @test value(stdized[2].kratios[1]) == 0.45/0.6517
        @test isapprox(σ(stdized[2].kratios[1]),0.037,atol=0.001)
        @test isequal(stdized[1].standard, mat"Fe2O3")
    end
end


text:
function loadSweeps(path)
    f = h5open(path, "r")
    dfDat = read(f["df"]["block0_values"])'
    headers = read(f["df"]["block0_items"])
    temp = DataFrame(dfDat, :auto)
    rename!(temp, headers)
    nSweeps = div(length(temp.index), length(unique(temp.index)))
    insertcols!(temp, 1,
        :sweep => repeat(1:nSweeps, inner=length(unique(temp.index))))

    return temp
end

function plotSweeps(df, sweep;window=[650, 1250], title=false;
    force_ylims=[-10,300], work_ylims=[-10,400], current_ylims=[-2000, 100],
     kwargs...)
    sweepGrps = groupby(df, :sweep)
    plotRegion = subset(sweepGrps[sweep],
        :ti => ByRow(>=(window[1])), :ti => ByRow(<=(window[2])))

    l = @layout [a{0.1h}; b{0.1h}; c{0.1h}; d{0.7h}]

    p1 = plot(plotRegion.tz, -1 .* plotRegion.position, ylim=(-6000,3000);
            kwargs...)
    p2 = plot(plotRegion.tin0, plotRegion.force, ylim=(-10,300); kwargs...)
    p3 = plot(plotRegion.tin0, plotRegion.work, ylim=(-10,400); kwargs...)
    p4 = plot(plotRegion.ti, plotRegion.i, ylim=(-1800,100); kwargs...)

    plt = plot(p1, p2, p3, p4, layout=l)
    if title == true
        annotate!(450, 2800, subplot=1, text("Position", :black,:left, 10))
        annotate!(450, 300, subplot=2, text("Force", :black,:left, 10))
        annotate!(450, 300, subplot=3, text("Work", :black,:left, 10))
        annotate!(450, 100, subplot=4, text("Current", :black,:left, 10))
    end

    return plt
end


text:
##################################################
# Test CoD from Zhao et al 2011
##################################################
bnd = MCBN.BayesNetDAI(3)
MCBN.add_edge!(bnd, 1, 2)
MCBN.add_edge!(bnd, 1, 3)
DAI.clearBackups!(bnd)
c10 = 0.4
c11 = 1. - c10
n13 = n12 = 0.2
d12 = d13 = 0.1
MCBN.set_factor!(bnd, 1, [c10, c11])
MCBN.set_factor!(bnd, 2, [1-n12, d12, n12, 1-d12])
MCBN.set_factor!(bnd, 3, [1-n13, d13, n13, 1-d13])

MCBN.check_bnd(bnd)
@show true_cod1 = MCBN.cod(bnd, 1, [2,3])
@show true_cod2 = MCBN.cod(bnd, 1, [2])
@show true_cod3 = MCBN.cod(bnd, 2, [1])

codbnd = bnd
bn = MCBN.BayesNetSampler(codbnd, MCBN.draw_data(codbnd, 20))

@show MCBN.ting_cod(bn, 1, [2,3])
@show MCBN.ting_cod(bn, 1, [2])
@show MCBN.ting_cod(bn, 2, [1])

mh = SAMC.set_energy_limits(bn)
mh.stepscale = 10
mh.thin = 100
mh.burn = 1000
SAMC.sample!(mh, 100000)

char_cod1(bnd::MCBN.BayesNetDAI) = MCBN.cod(bnd, 1, [2,3])
char_cod2(bnd::MCBN.BayesNetDAI) = MCBN.cod(bnd, 1, [2])
char_cod3(bnd::MCBN.BayesNetDAI) = MCBN.cod(bnd, 2, [1])

@show SAMC.posterior_e(char_cod1,mh)
@show SAMC.posterior_e(char_cod2,mh)
@show SAMC.posterior_e(char_cod3,mh)

@show map_cod = char_cod1(SAMC.mapvalue(mh))
@show map_cod = char_cod2(SAMC.mapvalue(mh))
@show map_cod = char_cod3(SAMC.mapvalue(mh))

chains = 10
burn = 500
thin = 100
iters = 100000
mh = SAMC.MHRecord[SAMC.MHRecord(deepcopy(bn), burn=burn, thin=thin) for i=1:chains]
SAMC.sample!(mh, div(iters,chains))

@show SAMC.posterior_e(char_cod1,mh)
@show SAMC.posterior_e(char_cod2,mh)
@show SAMC.posterior_e(char_cod3,mh)

@show map_cod = char_cod1(SAMC.mapvalue(mh))
@show map_cod = char_cod2(SAMC.mapvalue(mh))
@show map_cod = char_cod3(SAMC.mapvalue(mh))


text:
# Functions that implement formulae for the information measures

# Information measures are reviewed in:
# Timme, Nicholas; Alford, Wesley; Flecker, Benjamin; Beggs, John M. (2013-07-03).
# "Synergy, redundancy, and multivariate information measures: an experimentalist's perspective".
# Journal of Computational Neuroscience. 36 (2): 119–140.
# http://link.springer.com/article/10.1007%2Fs10827-013-0458-4

# TODO: Add further measures (dual total correlation, delta i)

export apply_entropy_formula, apply_conditional_entropy_formula, apply_mutual_information_formula,
	apply_conditional_mutual_information_formula, apply_interaction_information_formula,
	apply_total_correlation_formula, apply_specific_information_formula, apply_redundancy_formula,
	apply_unique_information_formula, apply_synergy_formula, apply_cross_entropy_formula

function remove_non_finite(x)
	return isfinite(x) ? x : zero(x)
end

function remove_negative(x)
	return x < 0 ? zero(x) : x
end

# Parameters:
# 	- probabilities, array of floats
# 	- base, number
function apply_entropy_formula(p::AbstractArray{T}, base::R) where {T <: AbstractFloat, R <: Real}
	return -sum(remove_non_finite.(p .* log.(base, p)))
end

# Parameters:
# 	- joint entropy of all variables, number
#	- entropy of conditioned variables, number
function apply_conditional_entropy_formula(entropy_xy, entropy_y)
	return entropy_xy - entropy_y
end

# Parameters:
# 	- joint probabilities, array of floats
# 	- probabilities (first variable), array of floats
# 	- probabilities (second variable), array of floats
# 	- base, number
function apply_mutual_information_formula(p_xy::AbstractArray{T}, p_x::AbstractArray{T}, p_y::AbstractArray{T}, base::R) where {T <: AbstractFloat, R <: Real}
	return sum(remove_non_finite.(p_xy .* log.(base, p_xy ./ (p_x .* p_y))))
end
# Parameters:
# 	- entropy of first variable, number
# 	- entropy of second variable, number
#	- joint entropy of both variables, number
function apply_mutual_information_formula(entropy_x, entropy_y, entropy_xy)
	return entropy_x + entropy_y - entropy_xy
end

# Parameters:
# 	- joint entropy of first and conditioned variables, number
# 	- joint entropy of second and conditioned variables, number
#	- joint entropy of all three variables, number
function apply_conditional_mutual_information_formula(entropy_xz, entropy_yz, entropy_xyz, entropy_z)
	return entropy_xz + entropy_yz - entropy_xyz - entropy_z
end

# Parameters:
# 	- mutual information of first two variables, number
# 	- conditional mutual information of first two variables on the third, number
function apply_interaction_information_formula(conditional_mutual_information, mutual_information)
	return conditional_mutual_information - mutual_information
end

# Parameters:
# 	- entropy of first variable, number
# 	- entropy of second variable, number
# 	- entropy of third variable, number
#	- joint entropy of all three variables, number
function apply_total_correlation_formula(entropy_x, entropy_y, entropy_z, entropy_xyz)
	return entropy_x + entropy_y + entropy_z - entropy_xyz
end

# Parameters:
# 	- joint probabilities, array of floats
# 	- probabilities (source), array of floats
# 	- probabilities (target), array of floats
# 	- dimension along which to sum, integer
# 	- base, number
function apply_specific_information_formula(p_xz, p_x, p_z, dim_sum, base)
	return vec(sum(remove_non_finite.((p_xz ./ p_z) .* log.(base, p_xz ./ (p_x .* p_z))), dims = dim_sum))
end

# Parameters:
# 	- joint probabilities (source 1 and target), array of floats
# 	- joint probabilities (source 2 and target), array of floats
# 	- probabilities (source 1), array of floats
# 	- probabilities (source 2), array of floats
# 	- probabilities (target), array of floats
# 	- dimensions along which to sum, tuple of integers
# 	- base, number
function apply_redundancy_formula(p_xz::AbstractArray{T}, p_yz::AbstractArray{T}, p_x::AbstractArray{T}, p_y::AbstractArray{T}, p_z::AbstractArray{T}, dim_sum::Tuple{I, I}, base::R) where {T <: AbstractFloat, R <: Real, I <: Integer}
	minimum_specific_information = min.(
		apply_specific_information_formula(p_xz, p_x, p_z, dim_sum[1], base),
		apply_specific_information_formula(p_yz, p_y, p_z, dim_sum[2], base)
	)
	return sum(vec(p_z) .* vec(minimum_specific_information))
end
# Parameters:
# 	- probabilities (target), array of floats
# 	- specific information of source 1 and target, array of floats
# 	- specific information of source 2 and target, array of floats
# 	- base, number
function apply_redundancy_formula(p_z::AbstractArray{T}, specific_information_1::AbstractArray{T}, specific_information_2::AbstractArray{T}, base::R) where {T <: AbstractFloat, R <: Real}
	minimum_specific_information = min.(specific_information_1, specific_information_2)
	return sum(vec(p_z) .* vec(minimum_specific_information))
end

# Parameters:
# 	- mutual information of source 1 and target, number
# 	- redundancy of both sources and target, number
function apply_unique_information_formula(mutual_information, redundancy)
	# Rounding errors may lead to slightly negative results
	return remove_negative.(mutual_information - redundancy)
end

# Parameters:
# 	- interaction information of both sources and target, number
# 	- redundancy of both sources and target, number
function apply_synergy_formula(interaction_information, redundancy)
	return interaction_information + redundancy
end

# Parameters:
#	- probabilities of first variable
#	- probabilities of second variable
#	- base of exponent
function apply_cross_entropy_formula(p_x::AbstractArray{T}, p_y::AbstractArray{T}, base::R) where {T <: AbstractFloat, R <: Real}
    return -sum(remove_non_finite.(p_x .* log.(base, p_y)))
end


text:
#!/usr/bin/env julia

# Converts data in MAT files (created using raw2mat.jl) to png
# thumbnails.

# e.g. thumbnail Sv38 ./
# e.g. thumbnail.jl ntheta120 *.mat

using EchogramImages
using FileIO
using MAT

function readvar(filename, varname)
    file = matopen(filename)
    v = read(file, varname)
    close(file)
    return v
end

function main(args)

    if isempty(args)
        info("thumbnail.jl <varname> <matfiles>")
        exit(1)
    end

    varname = args[1] # The variable name, e.g. Sv38

    filenames = args[2:end]

    if length(filenames) == 1
        dir = filenames[1]
        if isdir(dir)
            filenames =  filter(x->endswith(x,".mat"), readdir(dir))
            filenames = [joinpath(dir,x) for x in filenames]
        end
    end

    for filename in filenames

        out ="$(basename(filename))-$varname.png"
        if isfile(out)
            info("Skipping $out")
        else
            try
                info("Processing $filename ...")

                d = readvar(filename, varname)

                if startswith(varname,"Sv")
                    img = imagesc(d, vmin=-95, vmax = -40)
                else
                    img = imagesc(d)
                end

                info("Writing $out ...")
                save(out, img)
            catch ex
                warn(ex)
            end
        end
    end
end

main(ARGS)


text:
# Classical Multidimensional Scaling

## convert Gram matrix to Distance matrix

function gram2dmat!{DT}(D::AbstractMatrix{DT}, G::AbstractMatrix)
    # argument checking
    m = size(G, 1)
    n = size(G, 2)
    m == n || error("D should be a square matrix.")
    size(D) == (m, n) ||
        throw(DimensionMismatch("Sizes of D and G do not match."))

    # implementation
    for j = 1:n
        for i = 1:j-1
            @inbounds D[i,j] = D[j,i]
        end
        D[j,j] = zero(DT)
        for i = j+1:n
            @inbounds D[i,j] = sqrt(G[i,i] + G[j,j] - 2 * G[i,j])
        end
    end
    return D
end

gram2dmat{T<:Real}(G::AbstractMatrix{T}) = gram2dmat!(similar(G, Base.momenttype(T)), G)

## convert Distance matrix to Gram matrix

function dmat2gram!{GT}(G::AbstractMatrix{GT}, D::AbstractMatrix)
    # argument checking
    n = Compat.LinAlg.checksquare(D)
    size(G) == (n, n) ||
        throw(DimensionMismatch("Sizes of G and D do not match."))

    # implementation
    u = zeros(GT, n)
    s = 0.0
    for j = 1:n
        s += (u[j] = sum(abs2, view(D,:,j)) / n)
    end
    s /= n

    for j = 1:n
        for i = 1:j-1
            @inbounds G[i,j] = G[j,i]
        end
        for i = j:n
            @inbounds G[i,j] = (u[i] + u[j] - abs2(D[i,j]) - s) / 2
        end
    end
    return G
end

dmat2gram{T<:Real}(D::AbstractMatrix{T}) = dmat2gram!(similar(D, Base.momenttype(T)), D)

## classical MDS

function classical_mds{T<:Real}(D::AbstractMatrix{T}, p::Int;
        dowarn::Bool=true)

    n = size(D, 1)
    m = min(p, n) #Actual number of eigenpairs wanted

    G = dmat2gram(D)

    #Get m largest eigenpairs
    E = eigfact!(Symmetric(G))

    #Sometimes dmat2gram produces a negative definite matrix, and the sign just
    #needs to be flipped. The heuristic to check for this robustly is to check
    #if there is a negative eigenvalue of magnitude larger than the largest
    #positive eigenvalue, and flip the sign of eigenvalues if necessary.

    mineig, maxeig = extrema(E[:values])
    if mineig < 0 && abs(mineig) > abs(maxeig)
        #do flip
        ord = sortperm(E.values)
        v = -E[:values][ord[1:m]]
    else
        ord = sortperm(E.values; rev=true)
        v = E[:values][ord[1:m]]
    end

    for i = 1:m
        if v[i] > 0
            v[i] = √v[i]
        else #Keeping all remaining eigenpairs would not change solution (if 0)
             #or make the answer _worse_ (if <0).
             #The least squares solution would want to throw all these away.
            dowarn && warn("Gramian has only $(i-1) positive eigenvalue(s)")
            m = i-1
            ord = ord[1:m]
            v = v[1:m]
            break
        end
    end

    #Check if the last considered eigenvalue is degenerate
    if m>0
        nevalsmore = sum(abs.(E[:values][ord[m+1:end]] .- v[m]^2) .< n*eps())
        nevals = sum(abs.(E[:values] .- v[m]^2) .< n*eps())
        if nevalsmore > 1
            dowarn && warn("The last eigenpair is degenerate with $(nevals-1) others; $nevalsmore were ignored. Answer is not unique")
        end
    end
    U = E[:vectors][:, ord[1:m]]
    scale!(U, v)

    #Add trailing zero coordinates if dimension of embedding space (p) exceeds
    #number of eigenpairs used (m)
    if m < p
        U = [U zeros(n, p-m)]
    end

    U' #Return each coordinate in a column
end


text:
using Libz, BufferedStreams, Random, Test

Random.seed!(0x123456)

@testset "Source" begin
    function test_round_trip(data)
        return data == read(data |> ZlibDeflateInputStream |> ZlibInflateInputStream)
    end

    @test test_round_trip(UInt8[])
    @test test_round_trip(rand(UInt8, 1))
    @test test_round_trip(rand(UInt8, 1000000))
    @test test_round_trip(zeros(UInt8, 1000000))

    function test_round_trip2(data, bufsize, raw, gzip, reset_on_end)
        return data == read(
            ZlibInflateInputStream(
                ZlibDeflateInputStream(data, bufsize=bufsize, raw=raw, gzip=gzip),
                bufsize=bufsize, raw=raw, gzip=gzip, reset_on_end=reset_on_end))
    end

    for bufsize in [1, 5, 10, 50, 128*2^10], gzip in [false, true], reset_on_end in [false, true]
        @test test_round_trip2(UInt8[], bufsize, false, gzip, reset_on_end)
        @test test_round_trip2(rand(UInt8, 1), bufsize, false, gzip, reset_on_end)
        @test test_round_trip2(rand(UInt8, 1000000), bufsize, false, gzip, reset_on_end)
        @test test_round_trip2(zeros(UInt8, 1000000), bufsize, false, gzip, reset_on_end)
    end

    # raw deflate
    for bufsize in [1, 5, 10, 50, 128*2^10],  reset_on_end in [false, true]
        @test test_round_trip2(UInt8[], bufsize, true, false, reset_on_end)
        @test test_round_trip2(rand(UInt8, 1), bufsize, true, false, reset_on_end)
        @test test_round_trip2(rand(UInt8, 1000000), bufsize, true, false, reset_on_end)
        @test test_round_trip2(zeros(UInt8, 1000000), bufsize, true, false, reset_on_end)
    end

    @test_throws ArgumentError ZlibDeflateInputStream(UInt8[], bufsize=0)
    @test_throws ArgumentError ZlibDeflateInputStream(UInt8[], level=10)
    @test_throws ArgumentError ZlibInflateInputStream(UInt8[], bufsize=0)
    @test_throws ErrorException read(ZlibInflateInputStream([0x00, 0x01]))

    # check state transition
    stream = ZlibInflateInputStream(ZlibDeflateInputStream([0x00,0x00], bufsize=1), bufsize=1)
    @test stream.source.state === Libz.initialized
    # read 1 byte
    read(stream, UInt8)
    @test stream.source.state === Libz.inprogress
    # read the rest
    read(stream)
    @test stream.source.state === Libz.finished
    # close and release resources
    close(stream)
    @test stream.source.state === Libz.finalized
    # close again
    try
        close(stream)
        @test true
    catch
        @test false
    end
    @test stream.source.state === Libz.finalized
end

@testset "Sink" begin
    function test_round_trip(data)
        outbuf = IOBuffer()
        stream = outbuf |> ZlibInflateOutputStream |> ZlibDeflateOutputStream
        write(stream, data)
        flush(stream)
        actual = take!(outbuf)
        return actual == data
    end

    @test test_round_trip(UInt8[])
    @test test_round_trip(rand(UInt8, 1))
    @test test_round_trip(rand(UInt8, 1000000))

    function test_round_trip2(data, bufsize, raw, gzip)
        outbuf = IOBuffer()
        stream = ZlibDeflateOutputStream(
            ZlibInflateOutputStream(outbuf, bufsize=bufsize, raw=raw, gzip=gzip),
            bufsize=bufsize, raw=raw, gzip=gzip)
        write(stream, data)
        flush(stream)
        actual = take!(outbuf)
        return actual == data
    end

    for bufsize in [1, 5, 10, 50, 128*2^10], gzip in [false, true]
        @test test_round_trip2(UInt8[], bufsize, false, gzip)
        @test test_round_trip2(rand(UInt8, 1), bufsize, false, gzip)
        @test test_round_trip2(rand(UInt8, 1000000), bufsize, false, gzip)
        @test test_round_trip2(zeros(UInt8, 1000000), bufsize, false, gzip)
    end

    for bufsize in [1, 5, 10, 50, 128*2^10]
        @test test_round_trip2(UInt8[], bufsize, true, false)
        @test test_round_trip2(rand(UInt8, 1), bufsize, true, false)
        @test test_round_trip2(rand(UInt8, 1000000), bufsize, true, false)
        @test test_round_trip2(zeros(UInt8, 1000000), bufsize, true, false)
    end

    @test_throws ArgumentError ZlibDeflateOutputStream(UInt8[], bufsize=0)
    @test_throws ArgumentError ZlibDeflateOutputStream(UInt8[], level=10)
    @test_throws ArgumentError ZlibInflateOutputStream(UInt8[], bufsize=0)

    deflated = read(ZlibDeflateInputStream(Vector{UInt8}(codeunits("foo"))))
    buf = IOBuffer()
    out = ZlibInflateOutputStream(buf)
    BufferedStreams.writebytes(out, deflated, length(deflated), true)
    flush(out)
    @test String(take!((buf))) == "foo"

    @testset "sink vector" begin
        # Python: zlib.compress(bytearray([0x40, 0x41, 0x42]))
        x = [0x40, 0x41, 0x42]
        y = b"x\x9cspt\x02\x00\x01\x87\x00\xc4"

        out = ZlibDeflateOutputStream(UInt8[], gzip=false)
        write(out, x)
        flush(out)
        @test take!(out.sink.output) == y

        out = ZlibInflateOutputStream(UInt8[], gzip=false)
        write(out, y)
        flush(out)
        @test take!(out.sink.output) == x
    end

    # check state transition
    buf = IOBuffer()
    stream = ZlibDeflateOutputStream(ZlibInflateOutputStream(buf, bufsize=1), bufsize=1)
    @test stream.sink.state === Libz.initialized
    # write 2 bytes
    write(stream, [0x00, 0x01])
    @test stream.sink.state === Libz.inprogress
    # flush data
    flush(stream)
    @test stream.sink.state === Libz.finished
    # close and release resources
    close(stream)
    @test stream.sink.state === Libz.finalized
    # close again
    try
        close(stream)
        @test true
    catch
        @test false
    end
    @test stream.sink.state === Libz.finalized
end

@testset "Inflate/Deflate" begin
    data = rand(UInt8, 100000)
    @test Libz.inflate(Libz.deflate(data)) == data
end

@testset "Checksums" begin
    # checking correctness isn't our job, just make sure they're usable
    data = rand(UInt8, 100000)

    c32 = crc32(data)
    @test isa(c32, UInt32)

    a32 = adler32(data)
    @test isa(a32, UInt32)

    @test crc32(BufferedInputStream(IOBuffer(data), 1024)) === c32
    @test adler32(BufferedInputStream(IOBuffer(data), 1024)) === a32
end

@testset "Concatenated gzip files" begin
    filepath = joinpath(dirname(@__FILE__), "foobar.txt.gz")
    s = read(open(filepath) |> ZlibInflateInputStream, String)
    @test s == "foo\nbar\n"
end

@testset "Error" begin
    @test_throws ErrorException Libz.zerror(Libz.Z_ERRNO)
    @test_throws ErrorException Libz.zerror(Libz.Z_STREAM_ERROR)
    @test_throws ErrorException Libz.zerror(Libz.Z_DATA_ERROR)
    @test_throws ErrorException Libz.zerror(Libz.Z_MEM_ERROR)
    @test_throws ErrorException Libz.zerror(Libz.Z_BUF_ERROR)
    @test_throws ErrorException Libz.zerror(Libz.Z_VERSION_ERROR)
end

@testset "Low-level APIs" begin
    zstream = Libz.ZStream()
    ret = Libz.init_inflate!(zstream, 15)
    @test ret == Libz.Z_OK
    zstream = Libz.ZStream()
    @test_throws ErrorException (Libz.@zcheck Libz.init_inflate!(zstream, 100))
end


text:
const c₂                 = 1.4387769
const cMassMol           = 1.66053873e-30 #kilograms per molecule for unit molec. mass
const cSqrtLn2divSqrtPi  = 0.469718639319144059835 #√(ln2/π)
const cLn2               = 0.6931471805599 #ln2
const cSqrtLn2           = 0.8325546111577 #√(ln2)
const cSqrt2Ln2          = 1.1774100225 #√(2ln2)
const cc_                = 2.99792458e8 #speed of light [m/s]
const cBolts_            = 1.3806503e-23 #Boltzmann const. [J/K]
const p_ref              = 1013.25  # reference pressure [hPa]
const t_ref              = 296.0    # reference temperature [K]
const nm_per_m           = 1.0e7
# Note: ν stands for wavenumber in the following (NOT frequency)
function apply_lineshape!(Δνᵢ, σᵢ,  # discrete transitions
                λ₀,                 # incident  wavelength [nm]
                # Note:  λ₀ can either be used to denote an individual
                # wavelength or as a representative (e.g. central)
                # wavelength for a specific band
                Δν_out,             # Output grid (equidistant)
                σ_out,              # σ at output grid
                pressure::Real,     # actual pressure [hPa]
                temperature::Real,  # Temperature (K)
                molMass;
                wavelength_flag::Bool=false)

    # Notify user of wavelength grid
    if (wavelength_flag)
        @info """
        Note: Rayleigh/Raman Cross-section reported to wavelength grid (nm)
        """  maxlog = 5
    end

    # Max min range (ignoring wing cutoff here)
    grid_max = maximum(Δν_out)
    grid_min = minimum(Δν_out)

    # Interpolators from grid bounds to index values
    grid_idx_interp_low  = LinearInterpolation(Δν_out, 1:1:length(Δν_out), extrapolation_bc=1)
    grid_idx_interp_high = LinearInterpolation(Δν_out, 1:1:length(Δν_out), extrapolation_bc=length(Δν_out))

    # Loop through all transition lines:
    for j in eachindex(Δνᵢ)
        # Test that this ν lies within the grid
        if grid_min < Δνᵢ[j] < grid_max

            ν = Δνᵢ[j] + nm_per_m/λ₀#13500.0 #Dummy for now #Suniti

            # Compute Doppler HWHM, ν still needs to be supplied, @Suniti?:
            γ_d = ((cSqrt2Ln2 / cc_) * sqrt(cBolts_ / cMassMol) * sqrt(temperature) * ν / sqrt(molMass))

            @show γ_d, Δνᵢ[j]

            # line intensity
            S = σᵢ[j] *  ν^4 #Suniti

            wing_cutoff = 4γ_d

            # Calculate index range that this transition impacts
            ind_start = Int64(floor(grid_idx_interp_low(Δνᵢ[j] - wing_cutoff)))
            ind_stop  = Int64(ceil(grid_idx_interp_high(Δνᵢ[j] + wing_cutoff)))

            # Create views from the result and grid arrays
            result_view   = view(σ_out,  ind_start:ind_stop);
            grid_view     = view(Δν_out, ind_start:ind_stop);

            # Just Doppler broadening for now, can be modified with any line-shape later (using a kernel ideally)
            for I in eachindex(grid_view)
                # If we undersample the line-width, we have to make sure the integral is conserved (almost), TBD
                @inbounds result_view[I] += S * cSqrtLn2divSqrtPi * exp(-cLn2 * (((grid_view[I]) - Δνᵢ[j]) / γ_d)^2) / γ_d
            end
        end
    end
end

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
#######################################################################
# Install Jupyter kernel-spec file.

copy_config(src, dest) = cp(src, joinpath(dest, basename(src)), remove_destination=true)

"""
    installkernel(name, options...; specname=replace(lowercase(name), " ", "-")

Install a new Julia kernel, where the given `options` are passed to the `julia`
executable, and the user-visible kernel name is given by `name` followed by the
Julia version.

Internally, the Jupyter name for the kernel (for the `jupyter kernelspec`
command is given by the optional keyword `specname` (which defaults to
`name`, converted to lowercase with spaces replaced by hyphens),
followed by the Julia version number.

Both the `kernelspec` command (a `Cmd` object)
and the new kernel name are returned by `installkernel`.
For example:
```
kernelspec, kernelname = installkernel("Julia O3", "-O3")
```
creates a new Julia kernel in which `julia` is launched with the `-O3`
optimization flag.  The returned `kernelspec` command will be something
like `jupyter kernelspec` (perhaps with different path), and `kernelname`
will be something like `julia-O3-0.6` (in Julia 0.6).   You could
uninstall the kernel by running e.g.
```
run(`\$kernelspec remove -f \$kernelname`)
```
"""
function installkernel(name::AbstractString, julia_options::AbstractString...;
                   specname::AbstractString = replace(lowercase(name), " ", "-"))
    # Is IJulia being built from a debug build? If so, add "debug" to the description.
    debugdesc = ccall(:jl_is_debugbuild,Cint,())==1 ? "-debug" : ""

    # name of the Jupyter kernelspec directory to install
    spec_name = "$specname-$(VERSION.major).$(VERSION.minor)$debugdesc"

    juliakspec = joinpath(tempdir(), spec_name)
    try
        binary_name = Compat.Sys.iswindows() ? "julia.exe" : "julia"
        kernelcmd_array = String[joinpath(Sys.BINDIR,"$binary_name"), "-i",
                                 "--startup-file=yes", "--color=yes"]
        append!(kernelcmd_array, julia_options)
        ijulia_dir = get(ENV, "IJULIA_DIR", dirname(@__DIR__)) # support non-Pkg IJulia installs
        append!(kernelcmd_array, [joinpath(ijulia_dir,"src","kernel.jl"), "{connection_file}"])

        ks = Dict(
            "argv" => kernelcmd_array,
            "display_name" => name * " " * Base.VERSION_STRING * debugdesc,
            "language" => "julia",
        )

        destname = "kernel.json"
        mkpath(juliakspec)
        dest = joinpath(juliakspec, destname)

        open(dest, "w") do f
            # indent by 2 for readability of file
            write(f, JSON.json(ks, 2))
        end

        copy_config(joinpath(ijulia_dir,"deps","logo-32x32.png"), juliakspec)
        copy_config(joinpath(ijulia_dir,"deps","logo-64x64.png"), juliakspec)

        info("Installing $name kernelspec $spec_name")

        # remove these hacks when
        # https://github.com/jupyter/notebook/issues/448 is closed and the fix
        # is widely available -- just run `$jupyter kernelspec ...` then.
        kspec_cmd = String[] # keep track of the kernelspec command used
        try
            run(`$jupyter kernelspec install --replace --user $juliakspec`)
            push!(kspec_cmd, jupyter, "kernelspec")
        catch
            @static if Compat.Sys.isunix()
                run(`$jupyter-kernelspec install --replace --user $juliakspec`)
                push!(kspec_cmd, jupyter * "-kernelspec")
            end

            # issue #363:
            @static if Compat.Sys.iswindows()
                jupyter_dir = dirname(jupyter)
                jks_exe = ""
                if jupyter_dir == abspath(Conda.SCRIPTDIR)
                    jk_path = "$jupyter-kernelspec"
                    if isfile(jk_path * "-script.py")
                        jk_path *= "-script.py"
                    end
                    python = abspath(Conda.PYTHONDIR, "python.exe")
                else
                    jks_exe = joinpath(jupyter_dir, "jupyter-kernelspec.exe")
                    if !isfile(jks_exe)
                        jk_path = readchomp(`where.exe $jupyter-kernelspec`)
                        # jupyter-kernelspec should start with "#!/path/to/python":
                        python = strip(chomp(open(readline, jk_path, "r"))[3:end])
                        # strip quotes, if any
                        if python[1] == python[end] == '"'
                            python = python[2:end-1]
                        end
                    end
                end
                if isfile(jks_exe)
                    run(`$jks_exe install --replace --user $juliakspec`)
                    push!(kspec_cmd, jks_exe)
                else
                    run(`$python $jk_path install --replace --user $juliakspec`)
                    push!(kspec_cmd, python, jk_path)
                end
            end
        end

        return `$kspec_cmd`, spec_name
    finally
        rm(juliakspec, force=true, recursive=true)
    end
end


text:
#dummy file to get test to run with CI
module HolyLabRegistry
end

max_stars_repo_name: HolyLab/HolyLabRegistry

max_issues_repo_name: HolyLab/HolyLabRegistry

max_forks_repo_name: HolyLab/HolyLabRegistry


text:
module FPRS

using FFTW
using LinearAlgebra
using Random         # for densities
using DelimitedFiles # readdlm for loadElser

include("utilities.jl")
include("densities.jl")
include("projections.jl")
include("loadElser.jl")
include("bookkeeping.jl")

export
# Utilities
randBall
randRotationMatrix
# Densities
genGaussWGrid
genGauss
genRandGaussWGrid
genRandGauss
# Projections
pM
pM0
pP
pTa
pTn
eM
eM0
eP
eTa
eTn
# loadElser
loadElser
loadElserPure
# bookkeeping
pri
pri1D
lpri

end # module

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
import Term.Tables: Table
import Term.Layout: PlaceHolder

t = 1:5
data = hcat(t, ones(length(t)), string.(ones(length(t)) .* 100000))

t1 = Table(data)
t2 = Table(
    data;
    header = ["Num", "Const.", "Values"],
    header_style = "bold white",
    columns_style = ["dim", "bold", "red"],
)
t3 = Table(
    data;
    header = ["Num", "Const.", "Values"],
    header_style = "bold white",
    columns_style = ["dim", "bold", "red"],
    hpad = [1, 2, 5],
    columns_justify = [:center, :right, :left],
    vpad = 1,
)
t4 = Table(data; footer = ["get", "a", "footer"], footer_justify = :center)

t5 = Table(data; footer = sum, footer_justify = :center, footer_style = "dim bold")

t7 = Table(data; columns_widths = [25, 7, 7], footer = sum, box = :SIMPLE)

ph1 = PlaceHolder(25, 5)
ph2 = PlaceHolder(23, 9)
ph3 = PlaceHolder(22, 11)

data = Dict(
    "first\ncol." => [ph1, ph2, ph3],
    "second\ncol." => [ph3, ph2, ph3],
    "third\ncol." => [ph2, ph2, ph1],
)

t6 = Table(data)

X = rand(RNG, 5, 3)
t8 = Table(
    X;
    columns_widths = [12, 10, 22],
    hpad = 2,
    columns_justify = [:left, :center, :left],
)

# save tables as strings to files
tbls = [t1, t2, t3, t4, t5, t6, t7, t8]
for (i, t) in enumerate(tbls)
    tofile(string(t), "./txtfiles/table_$i.txt")
end

@testset "TABLE" begin
    for (i, t) in enumerate(tbls)
        @test fromfile("./txtfiles/table_$i.txt") == cleanstring(t)
    end
end


text:
module SimplePool

# simple scan into a list of free buffers

using ..CuArrays
using ..CuArrays: @pool_timeit

using CUDAdrv

using Base: @lock
using Base.Threads: SpinLock

const pool_lock = ReentrantLock()

## tunables

# how much larger a buf can be to fullfil an allocation request.
# lower values improve efficiency, but increase pressure on the underlying GC.
function max_oversize(sz)
    if sz <= 2^20       # 1 MiB
        # small buffers are fine no matter
        return typemax(Int)
    elseif sz <= 2^20   # 32 MiB
        return 2^20
    else
        return 2^22
    end
end

## block of memory

struct Block
    ptr::CuPtr{Nothing}
    sz::Int
end

Base.pointer(block::Block) = block.ptr
Base.sizeof(block::Block) = block.sz

@inline function actual_alloc(sz)
    ptr = CuArrays.actual_alloc(sz)
    block = ptr === nothing ? nothing : Block(ptr, sz)
end

function actual_free(block::Block)
    CuArrays.actual_free(pointer(block))
    return
end

## pooling

const available = Set{Block}()
const allocated = Dict{CuPtr{Nothing},Block}()

const freed = Vector{Block}()
const freed_lock = SpinLock()

function scan(sz)
    @lock pool_lock for block in available
        if sz <= sizeof(block) <= max_oversize(sz)
            delete!(available, block)
            return block
        end
    end
    return
end

function repopulate()
    blocks = @lock freed_lock begin
        isempty(freed) && return
        blocks = Set(freed)
        empty!(freed)
        blocks
    end

    @lock pool_lock begin
        for block in blocks
            @assert !in(block, available)
            push!(available, block)
        end
    end

    return
end

function reclaim(sz::Int=typemax(Int))
    repopulate()

    @lock pool_lock begin
        freed_bytes = 0
        while freed_bytes < sz && !isempty(available)
            block = pop!(available)
            freed_bytes += sizeof(block)
            actual_free(block)
        end
        return freed_bytes
    end
end

function pool_alloc(sz)
    block = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(true)
        end

        @pool_timeit "$phase.1 repopulate" repopulate()

        @pool_timeit "$phase.2 scan" begin
            block = scan(sz)
        end
        block === nothing || break

        @pool_timeit "$phase.3 alloc" begin
            block = actual_alloc(sz)
        end
        block === nothing || break

        @pool_timeit "$phase.4 reclaim + alloc" begin
            reclaim(sz)
            block = actual_alloc(sz)
        end
        block === nothing || break
    end

    return block
end

function pool_free(block)
    # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
    # and to simplify locking (preventing concurrent access during GC interventions)
    @lock freed_lock begin
        push!(freed, block)
    end
end

## interface

init() = return

function alloc(sz)
    block = pool_alloc(sz)
    if block !== nothing
        ptr = pointer(block)
        @lock pool_lock begin
            allocated[ptr] = block
        end
        return ptr
    else
        return nothing
    end
end

function free(ptr)
    block = @lock pool_lock begin
        block = allocated[ptr]
        delete!(allocated, ptr)
        block
    end
    pool_free(block)
    return
end

used_memory() = @lock pool_lock mapreduce(sizeof, +, values(allocated); init=0)

cached_memory() = @lock pool_lock mapreduce(sizeof, +, union(available, freed); init=0)

end


text:
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.

module FAME

include("Types.jl")
include("ChliLibrary.jl")
include("Command.jl")
include("Databases.jl")
include("Objects.jl")
include("Read.jl")
include("Write.jl")
include("Bridge.jl")

export version
"""
    version()

Return the version of FAME in use.
"""
function version()
    version = Ref{Cfloat}(0.0)
    @cfm_call_check(cfmver, (Ref{Cfloat},), version)
    return VersionNumber(Base.string(version.x))
end

"""
    close_chli()

Finalize FAME and unload the CHLI library. Normally, it should not be necessary
to call this function. If you want to restart FAME, e.g., in order to clear the
work database, call [`init_chli`](@ref).
"""
function close_chli()
    if chli.lib !== nothing
        @cfm_call_check(cfmfin)
        Libdl.dlclose(chli.lib)
        global chli = Chli{FameDatabase}(nothing)
    end
    return
end

"""
    init_chli()

Load the CHLI library and initialize FAME. If it is already loaded, we first
close it and then load it fresh. See also [`close_chli`](@ref).
"""
function init_chli()
    close_chli()
    global chli = Chli()

    chli.lib === nothing && return

    @cfm_call_check(cfmini)

    global FAME_INDEX_NC = @cfm_global(FAME_INDEX_NC, FameIndex)
    global FAME_INDEX_NA = @cfm_global(FAME_INDEX_NA, FameIndex)
    global FAME_INDEX_ND = @cfm_global(FAME_INDEX_ND, FameIndex)

    global FPRCNA = @cfm_global(FPRCNA, Cdouble)
    global FPRCNC = @cfm_global(FPRCNC, Cdouble)
    global FPRCND = @cfm_global(FPRCND, Cdouble)

    global FNUMNA = @cfm_global(FNUMNA, Cfloat)
    global FNUMNC = @cfm_global(FNUMNC, Cfloat)
    global FNUMND = @cfm_global(FNUMND, Cfloat)

    global FBOONA = @cfm_global(FBOONA, Cint)
    global FBOONC = @cfm_global(FBOONC, Cint)
    global FBOOND = @cfm_global(FBOOND, Cint)

    global FSTRNA = @cfm_global(FSTRNA, Cstring)
    global FSTRNC = @cfm_global(FSTRNC, Cstring)
    global FSTRND = @cfm_global(FSTRND, Cstring)

    return
end

function __init__()
    init_chli()
    atexit(close_chli)
end

end # module


text:
function groupedbar(; x::AbstractVector = Int[], y::AbstractVector = Int[], group::AbstractVector = Int[], horizontal::Bool = false)

    v = VegaVisualization(height = 300, width = 500, name = "groupedbar")

    #This seems like it should be a method, but so easy to do right here
    if typeof(x) == Vector{Char}
        x = String[string(s) for s in x]
    end

    add_data!(v, x = x, y = y, group = group)
    default_axes!(v)
    legend!(v)

    #Get unique table name
    table = v.data[1].name

    default_scales!(v, typeXaxis = "ordinal")
    v.scales[1].padding = 0.2
    v.scales[2].nice = true

    v.marks = [VegaMark(_type = "group")]
    v.marks[1].from = VegaMarkFrom(data = table,
                                   transform = [VegaTransform(Dict{Any, Any}("groupby" => ["x"], "type" => "facet"))])

    v.marks[1].marks = [VegaMark(_type = "rect")]
    v.marks[1].marks[1].properties = VegaMarkProperties(enter = VegaMarkPropertySet(fill = VegaValueRef(field = "group", scale = "group"),
                                                                                width = VegaValueRef(band = true, offset = -1, scale = "pos"),
                                                                                x = VegaValueRef(field = "group", scale = "pos"),
                                                                                y = VegaValueRef(field = "y", scale = "y"),
                                                                                y2 = VegaValueRef(scale = "y", value = 0)
                                                                                )
                                                        )
    v.marks[1].properties = VegaMarkProperties(enter = VegaMarkPropertySet(width = VegaValueRef(band = true, scale = "x"),
                                                                            x = VegaValueRef(field = "key", scale = "x")
                                                )
                                                )

    v.marks[1].scales = [VegaScale(name = "pos",
                                    _type = "ordinal",
                                    range = "width",
                                    domain = VegaDataRef(field = "group"))]

    if horizontal == true
        coord_flip!(v)
    end

    return v

end


text:
using Documenter, FilePathsBase, FilePathsBase.TestPaths

const SETUP_CODE = quote
    using FilePathsBase
    using FilePathsBase: /, join
end

DocMeta.setdocmeta!(FilePathsBase, :DocTestSetup, SETUP_CODE; recursive=true)

makedocs(
    modules=[FilePathsBase],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    pages=[
        "Home" => "index.md",
        "Design" => "design.md",
        "FAQ" => "faq.md",
        "API" => "api.md",
    ],
    repo="https://github.com/rofinn/FilePathsBase.jl/blob/{commit}{path}#L{line}",
    sitename="FilePathsBase.jl",
    authors="Rory Finnegan",
    # checkdocs = :exports,
    # strict = true,
)

deploydocs(
    repo = "github.com/rofinn/FilePathsBase.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
type XYZ <: AtomFileType end

file_exts(::Type{XYZ}) = (".xyz",)

function Atoms(::Type{XYZ}, num_atoms::Int)
    Atoms(num_atoms,
        :type => Vector{Symbol}(num_atoms),
        :coords => Matrix{Float64}(3, num_atoms),
        :molecule => Vector{Int}(num_atoms)
    )
end

function read_frame_header{T}(::Type{T}, ::Type{XYZ}, chunk::Chunk)
    frame_start = chunk.i

    num_atoms = read_int(chunk)
    read_newlines(chunk, 2)

    atom_data_start = chunk.i
    read_newlines(chunk, num_atoms)
    atom_data_size = chunk.i - atom_data_start

    frame_size = chunk.i - frame_start

    frame_header = FrameHeader(frame_start, frame_size, num_atoms, atom_data_start, atom_data_size)

    frame = Frame()

    return frame_header, frame
end

function read_atoms{T}(::Type{T}, ::Type{XYZ}, chunk::Chunk, atom::Atoms, indices::UnitRange{Int})
    types = atom.props[:type]::Vector{Symbol}
    coords = atom.props[:coords]::Matrix{T}
    molecules = atom.props[:molecule]::Vector{Int}

    for i in indices
        atom_type = read_symbol(chunk)
        read_whitespace(chunk)

        x = read_float(T, chunk)
        read_whitespace(chunk)

        y = read_float(T, chunk)
        read_whitespace(chunk)

        z = read_float(T, chunk)
        read_whitespace(chunk)

        molecule = read_int(chunk)
        read_newline(chunk)

        types[i] = atom_type
        coords[1, i] = x
        coords[2, i] = y
        coords[3, i] = z
        molecules[i] = molecule
    end
end

function write_atoms(::Type{XYZ}, io::IO, frame::Frame)
    haskey(frame, :atoms) || throw(ArgumentError("Need frame atoms"))

    atoms = frame[:atoms]
    haskey(atoms, :type) || throw(ArgumentError("Need atom types"))
    haskey(atoms, :coords) || throw(ArgumentError("Need atom coords"))

    num_atoms = length(atoms)
    println(io, num_atoms)

    haskey(frame, :timestep) && print(io, "t = ", frame[:timestep], ' ')

    if haskey(frame, :box_max)
        xhi, yhi, zhi = frame[:box_max]
        xlo, ylo, zlo = 0, 0, 0

        if haskey(frame, :box_min)
            xlo, ylo, zlo = frame[:box_min]
        end

        print(io, "box = ", xhi - xlo, ' ', yhi - ylo, ' ', zhi - zlo, ' ')
    end

    println(io)

    has_extra = haskey(atoms, :molecule) || haskey(atoms, :energy)

    for (i, atom) in enumerate(atoms)
        print(io, atom[:type], ' ', atom[:coords][1], ' ', atom[:coords][2], ' ', atom[:coords][3])

        haskey(atoms, :charge) && print(io, ' ', atom[:charge])

        if has_extra
            print(io, " #")
            haskey(atoms, :molecule) && print(io, ' ', atom[:molecule])
            haskey(atoms, :energy) && print(io, ' ', atom[:energy])
        end

        println(io)
    end
end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
module BioMedQuery

# --------Database Utilities---------
module DBUtils
export init_mysql_database,
       assemble_vals,
       assemble_cols_and_vals_select,
       assemble_cols_and_vals,
       assemble_cols_and_vals_string,
       insert_row!,
       db_select,
       db_query,
       db_clean_string,
       create_server,
       insert_df!,
       assemble_cols,
       col_match,
       set_innodb_checks!,
       select_all_tables

include("DBUtils/DBUtils.jl")
end

#-----------PubMed----------------
module PubMed

# sql utilities
export  create_tables!,
        create_pmid_table!,
        all_pmids,
        abstracts,
        get_article_mesh,
        get_article_mesh_by_concept,
        db_insert!,
        abstracts_by_year,
        add_mysql_keys!,
        drop_mysql_keys!
include("PubMed/pubmed_sql_utils.jl")

# eutils -> sql
export  save_efetch!,
        save_pmids!
include("PubMed/eutils_sql_save.jl")

# eutils -> dfs -> csv
export parse_articles,
        dfs_to_csv
include("PubMed/pubmed_to_csv.jl")

# citation formats
export  CitationOutput,
        citations_endnote,
        citations_bibtex,
        save_efetch!
include("PubMed/citation_manager.jl")

end

#--------Clinical Trials------------
module CT
export search_ct
include("CT/CT.jl")
end
#-----------------------------------

#--------Processes------------
module Processes
export pubmed_search_and_save!,
       pubmed_pmid_search,
       pubmed_pmid_search_and_save!,
       pubmed_search_and_parse,
       map_mesh_to_umls_async!,
       map_mesh_to_umls_async,
       umls_semantic_occurrences,
       filter_mesh_by_concept,
       export_citation,
       load_medline!
include("Processes/pubmed_search_and_save.jl")
include("Processes/pubmed_mesh_to_umls_map.jl")
include("Processes/pubmed_occurrance_filtering.jl")
include("Processes/pubmed_export_citations.jl")
include("Processes/medline_load.jl")

end
#----------------------------------- d

end #BioMedQuery


text:
module Statifiers

using StatsBase: std, skewness, kurtosis, variation, sem, mad,  entropy, summarystats, autocor, pacf, rle, quantile
using Dates
using DataFrames
using Random
using Statistics

export fit!,transform!
export Statifier

export fullstat,statifierrun

using TSML.TSMLTypes
import TSML.TSMLTypes.fit! # to overload
import TSML.TSMLTypes.transform! # to overload
using TSML.Utils

# Transforms instances with nominal features into one-hot form
# and coerces the instance matrix to be of element type Float64.
mutable struct Statifier <: Transformer
  model
  args

  function Statifier(args=Dict())
    default_args = Dict(
      :processmissing => true
    )
    new(nothing, mergedict(default_args, args))
  end
end

function fit!(st::Statifier, features::T=[], labels::Vector=[]) where {T<:Union{Vector,Matrix,DataFrame}}
  typeof(features) <: DataFrame || error("Statifier.fit!: data should be a dataframe: Date,Val ")
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  st.model = st.args
end

function transform!(st::Statifier, features::T=[]) where {T<:Union{Vector,Matrix,DataFrame}}
  features != [] || return DataFrame()
  typeof(features) <: DataFrame || error("Statifier.fit!: data should be a dataframe: Date,Val ")
  ncol(features) == 2 || error("dataframe must have 2 columns: Date, Val")
  sum(names(features) .== (:Date,:Value))  == 2 || error("wrong column names")
  fstat = fullstat(features[:Value])
  timestat = timevalstat(features)
  if st.args[:processmissing] == true
    # full namedtuple: stat1,stat2,bstat
    hcat(timestat,fstat...)
  else
    hcat(timestat,fstat.stat1,fstat.stat2)
  end
end

function timevalstat(features::DataFrame)
  ldates=features[:Date] |> skipmissing |> collect
  totalhours = sum(diff(ldates)).value/1000/3600 # to hours
  lcount = length(ldates)
  timestart=first(features[:Date])
  timeend=last(features[:Date])
  sfreq = totalhours/lcount
  dftime = DataFrame(tstart=timestart,tend=timeend,sfreq=sfreq)
end

function fullstat(dat::Vector)
  data = skipmissing(dat) |> collect
  lcount = length(data)
  lmax = maximum(data)
  lmin = minimum(data)
  lsm = summarystats(data)
  q1 = quantile(data,0.1)
  q2 = quantile(data,0.2)
  q8 = quantile(data,0.8)
  q9 = quantile(data,0.9)
  lks = kurtosis(data)
  lsk = skewness(data)
  lvar = variation(data)
  lsem = sem(data)
  lentropy = entropy(data)
  # assume 24-hour period
  _autolags = 1:minimum([24,length(data)-1])
  _pacflags = 1:minimum([24,div(length(data),2)-1])
  lautocor = autocor(data,_autolags) .|> abs2 |> sum |> sqrt
  lpacf = pacf(data,_pacflags;method=:yulewalker) .|> abs2 |> sum |> sqrt
  lbrle = rlestatmissingblocks(dat)
  df1=DataFrame(count=lcount,max=lmax,min=lmin,median=lsm.median,
                mean=lsm.mean,q1=q1,q2=q2,q25=lsm.q25,
                q75=lsm.q75,q8=q8,q9=q9)
  df2=DataFrame(kurtosis=lks,skewness=lsk,variation=lvar,entropy=lentropy,
            autocor=lautocor,pacf=lpacf)
  df3=DataFrame(bmedian=lbrle.bmedian,bmean=lbrle.bmean,bq25=lbrle.bq25,bq75=lbrle.bq75,
                bmin=lbrle.bmin,bmax=lbrle.bmax)
  return (stat1=df1,stat2=df2,bstat=df3)
end

function rlestatmissingblocks(dat::Vector{T}) where {T<:Union{AbstractFloat,Integer,Missing}}
  data = deepcopy(dat)
  # replace missing with dummy value for rle processing
  indxmissing=findall(x->ismissing(x),data)
  dummy = -99999
  if eltype(data) <: Union{AbstractFloat,Missing}
    data[indxmissing] .= dummy
  elseif eltype(data) <: Union{Integer,Missing}
    data[indxmissing] .= dummy
  else
    error("not a valid type")
  end
  _rle = rle(data)
  indx=findall(x->x == dummy,_rle[1])
  blocks = _rle[2][indx]
  _sm = summarystats(blocks)
  return (bq25=_sm.q25,bmean=_sm.mean,bmedian=_sm.median,bq75=_sm.q75,bmax=_sm.max,bmin=_sm.min,bmiss=_sm.nmiss)
end

end


text:
struct Peek{Rows}
    rows::Rows
    maximum_length::Int
end

export Peek
"""
    Peek(rows, maximum_length = 4)

Peek an iterator which returns named tuples. Will show no more than `maximum_length` rows.

```jldoctest Peek
julia> using LightQuery

julia> @name Peek(Rows((a = 1:5, b = 5:-1:1)))
Showing 4 of 5 rows
| `a` | `b` |
| ---:| ---:|
|   1 |   5 |
|   2 |   4 |
|   3 |   3 |
|   4 |   2 |
```
"""
function Peek(rows)
    Peek(rows, 4)
end

function make_any(values)
    Any[values...]
end
@inline function justification(column)
    :r
end
function show(output::IO, peek::Peek)
    rows = peek.rows
    maximum_length = peek.maximum_length
    if isa(IteratorSize(rows), Union{HasLength, HasShape})
        if length(rows) > maximum_length
            println(output, "Showing $(maximum_length) of $(length(rows)) rows")
        end
    else
        println(output, "Showing at most $(maximum_length) rows")
    end
    columns = make_columns(take(rows, maximum_length))
    rows = map(make_any, zip(map(value, columns)...))
    pushfirst!(rows, make_any(map_unrolled(key, columns)))
    show(output, MD(Table(
        rows,
        make_any(map_unrolled(justification, columns))
    )))
end

@inline function reduce_values(a_function, columns, row1, row2)
    map_unrolled(tuple,
        columns,
        map_unrolled(
            a_function,
            map_unrolled(value, columns(row1)),
            map_unrolled(value, columns(row2))
        )
    )
end

"""
    function reduce_rows(rows, a_function, columns...)

Reduce a function over each of `columns` in `rows`.

```jldoctest
julia> using LightQuery

julia> using Test: @inferred

julia> @name @inferred reduce_rows(Rows((a = [1, 1], b = [1.0, 1.0])), +, :a, :b)
((`a`, 2), (`b`, 2.0))
```
"""
@inline function reduce_rows(rows, a_function, columns...)
    reduce(
        let a_function = a_function, columns = columns
            @inline function reduce_rows_capture(row1, row2)
                reduce_values(a_function, columns, row1, row2)
            end
        end,
        rows
    )
end
export reduce_rows