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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ b97724e4-d7b0-4085-b88e-eb3c5bcbe441
begin
    using Pkg
	Pkg.activate(joinpath(Pkg.devdir(), "MLCourse"))
    using StatsBase, DataFrames, MLCourse, Plots
end

# ╔═╡ b3793299-c916-40d3-bd87-31153fc3781a
using PlutoUI, PlutoUI.BuiltinsNotebook.HypertextLiteral; PlutoUI.TableOfContents()

# ╔═╡ 1a17e9b2-a439-4521-842c-96ebe0378919
using ReinforcementLearning

# ╔═╡ ce405f97-6d60-4ae4-b183-79e6c88d9811
md"# Chasse au trésor"

# ╔═╡ b0b89b30-38f1-45c1-8c7a-796ea2f41e8d
md"## Learning Q-Values with Monte Carlo Estimation
In Monte Carlo estimation, the Q-values ``Q(s, a)`` of state-action pair ``(s, a)`` are updated every time action ``a`` is taken in state ``s`` such that the current estimate is given by the average of all cumulative returns obtained so far when taking action ``a`` in state ``s``.

Additionally to an implementation of Monte Carlo Estimation of Q-values you see in the following cell a way to organize code in Julia using custom types.
We define the structure `MCLearner` to have the fields `N` and `Q` (with some default values adapted to our chasse au trésor) and some functions to update objects of type `MCLearner`. We also write doc-strings for our new functions. Open the live docs and click on the `update!` function to see the rendered doc-string.
"

# ╔═╡ 7c756af7-08a0-4878-8f53-9eaa79c16c1a
begin
    Base.@kwdef struct MCLearner # this defines a type MCLearner
        ns = 7                 # default initial number of states
        na = 2                 # default initial number of actions
        N = zeros(Int, na, ns) # default initial counts
        Q = zeros(na, ns)      # default initial Q-values
    end
	"""
	    update!(l::MCLearner, a, s, r)

	Updates an `MCLearner` object `l` for a given
	action `a`, state `s` and reward `r`.
	"""
    function update!(learner::MCLearner, a, s, r)
        N, Q = learner.N, learner.Q
        n = N[a, s]
        Q[a, s] = 1/(n+1) * (n * Q[a, s] + r) # efficient averaging
        N[a, s] += 1
        learner
    end
	"""
	    update!(l::MCLearner, episode)

	Updates an `MCLearner` object `l` for a full episode.
	"""
    function update!(learner::MCLearner, episode)
		G = 0 # initialise the cumulative reward
        for (s, a, r, ) in reverse(episode)
			G += r
            update!(learner, a, s, G)
        end
        learner
    end
    reset!(learner::MCLearner) = learner.Q .= learner.N .= 0
end;

# ╔═╡ 6dced47e-cc0f-4ae5-bdd0-a551c6d5a5b3
md"## An Epsilon-Greedy Agent

Up until now we took the actions. Here we introduce a very simple agent that takes action currently believed to be the best one with probability ``1 - \epsilon`` and a random action with probability ``\epsilon``.
"

# ╔═╡ 64863a37-1c25-4d1c-9f2d-33d87b4039a3
function epsilon_greedy_policy(q; ϵ = .25)
    if rand() < ϵ # this happens with probability ϵ
        rand(1:length(q))  # random exploration
    else
        rand(findall(==(maximum(q)), q)) # greedy exploitation with ties broken
    end
end

# ╔═╡ 0537554c-76e4-4827-b92a-bf5b944685c7
md"Now we let this agent play the game for 10^6 episodes and learn the Q-values with Monte Carlo estimation."

# ╔═╡ dc7c5cb0-30b2-4427-8aeb-f312a88effd1
md"# Q-Learning

Q-learning is an alternative method for updating Q-values that relies on the idea that ``Q(s_t, a_t)`` for the optimal policy can be estimated by knowing ``r_t`` and ``\max_a Q(s_{t+1}, a)``.
"

# ╔═╡ e39227c8-8148-43cb-9351-774682b65646
begin
    Base.@kwdef struct QLearner # this defines a type QLearner
        ns = 7                 # default initial number of states
        na = 2                 # default initial number of actions
        Q = zeros(na, ns)      # default initial Q-values
        λ = 0.01
    end
	"""
	    update!(l::QLearner, a, s, r)

	Updates a `QLearner` object `l` for a given
	state `s`, action `a`, reward `r` and next state `s′`.
	"""
    function update!(learner::QLearner, s, a, r, s′)
        Q, λ = learner.Q, learner.λ
        δ = (r + maximum(Q[:, s′]) - Q[a, s])
        Q[a, s] = Q[a, s] + λ * δ
        learner
    end
    reset!(learner::QLearner) = learner.Q .= 0
end;

# ╔═╡ aff5b0d2-c5e5-4fda-a93c-54a8bca5187f
md"# Learning to Play Tic-Tac-Toe

In the following example we will learn to play Tic-Tac-Toe with a Monte Carlo Learner.
The Tic-Tac-Toe environment is loaded from the [`ReinforcementLearning` package](https://juliareinforcementlearning.org/). This package also provides multiple helper functions to interact with this environment, like `state`, `is_terminated`, `legal_action_space`, `ReinforcementLearning.reward`, `ReinforcementLearning.reset!` etc.
"

# ╔═╡ ae72b2c3-fe7b-4ffb-b45c-e589529209c7
tictactoe = TicTacToeEnv();

# ╔═╡ d094b8d7-5cc7-4bc0-a406-2e3c85f9fb60
md"To perform an action we can simply call the environment with an integer from 1 to 9.
1 means placing a cross `x` or nought `o` (depending on whose turn it is) at the top left corner, 3 is at the bottom left corner, 7 top right and 9 bottom right."

# ╔═╡ 4b90f08e-0183-4f8b-a6cd-e66b6938f4c8
legal_action_space(tictactoe) # now action 5 is no longer available

# ╔═╡ 59029032-1c91-4da1-a61b-6a56449dcd2c
md"Let us now actually play the game. You can choose different modes below to play against yourself (or another human player), against the computer (trained in self-play with our `MCLearner`; see below) or computer against computer. To advance the game when the computer plays against itself you have to use the `step` button."

# ╔═╡ d99a218a-56e8-4081-bd1a-ba7729f529cf
md"The computer player is using a greedy policy with Q-values learned in self-play with the code below."

# ╔═╡ e61c6d43-a097-4f44-a343-05371b4932ef
begin
    is_terminated(env::TicTacToeEnv) = RLBase.is_terminated(env)
    act!(env::TicTacToeEnv, a) = env(a)
    reset!(env::TicTacToeEnv) = RLBase.reset!(env)
    reward(env::TicTacToeEnv) = RLBase.reward(env)
    state(env::TicTacToeEnv) = RLBase.state(env, Observation{Int}())
end;

# ╔═╡ b6b835d1-8f84-4148-8d5b-c7aea6b0c312
md"""# Exercises
## Conceptual

Consider an agent that experiences episode 1:

```math
\begin{align*}
S_1 &= s_1, A_1 = a_1, R_2 = 2,\\
S_2 &= s_7, A_2 = a_2, R_3 = 1,\\
S_3 &= s_6, A_2 = a_2, R_4 = 4,\\
S_4 &= s_3, A_4 = a_1, R_5 = 0
\end{align*}
```

and episode 2:

```math
\begin{align*}
S_1 &= s_1, A_1 = a_1, R_2 = 1,\\
S_2 &= s_6, A_2 = a_1, R_3 = -1,\\
S_3 &= s_2, A_3 = a_2, R_4 = 1
\end{align*}
```.

(a) Compute ``Q(s_1, a_1)`` with Monte Carlo Estimation.

(b) Compute ``Q(s_1, a_1)`` with Q-Learning, learning rate ``\lambda = 0.5`` and initial ``Q(s, a) = 0`` for all ``s`` and ``a``.

## Applied

Consider the cliff-walking task shown below. This is a standard episodic task with an agent (red point) moving around in a grid world with start (bottom left) and goal state (green square at bottom right), and the usual actions causing movement up, down, right, and left. Reward is -1 on all transitions except those into the cliff region marked black and the goal state at the bottom right. Stepping into this region incurs a reward of -100. Episodes end when the agent steps into the cliff or reaches the goal at the bottom right; in this case `is_terminated(environment) == true` and `reset!(environment)` sends the agent back to the start state at the bottom left.

(a) Determine the cumulative reward that would be obtained when following the optimal strategy in this task.

(b) Train a Monte Carlo learner with epsilon-greedy exploration (``\epsilon = 0.1``) for ``10^3`` episodes. Keep track of the cumulative reward in each episode. *Hints:* You can use the `environment = CliffWalkingEnv()` and the usual functions like `state(environment)`, `reward(environment)` and `act!(environment, action)` etc. Because we do not know the length of the episodes, you can use the following pattern:
```julia
while true
    ...
    if is_terminated(environment)
        ...
        break
    end
end
```

(c) Plot the cumulative reward per episode, i.e. x-axis: episode, y-axis: cumulative reward.

(d) Use the following functions to plot the learned greedy policy and the Q values of the greedy policy (i.e. the maximal Q values for each state)
```julia
function plot_greedy_policy(environment, Q)
	[x[1] == 1 ? "<" : x[1] == 2 ? ">" : x[1] == 3 ? "∧" : "∨"
		for x in argmax(Q, dims = 1)][environment.params.states] |>
	string |>
	s -> replace(s, r"[\\"\[\]]" => "") |>
	s -> replace(s, "; " => "\n") |>
	Text
end
function plot_q(environment, Q)
    heatmap(maximum(Q, dims = 1)[environment.params.states], yflip = true)
end
```

(e) Explain why the Monte Carlo learner failed to learn a good policy.

(f) With `reset!(environment)` the agent is put to the start state at the end of each episode. You can use instead the following function to reset the agent to a random position.
```julia
function random_reset!(env::CliffWalkingEnv)
	env.position = rand(CartesianIndices((4, 12)))
end
```
Run ``10^5`` episodes with this random resetting of the environment and plot again the learned greedy policy and Q values.

(g) Repeat (b-d) for a Q learner.

(h) Explain, why the cumulative reward of the Q learner is typically lower than the cumulative reward of the optimal policy, even though the Q learner finds the optimal greedy policy.

(i) Explain, why the Monte Carlo learner with random reset learns a greedy policy that does not follow as closely the cliff as the greedy policy learned with Q learning.
"""

# 2. Devise three example tasks of your own that fit into the MDP framework, identifying for each its states, actions, and rewards. Make the three examples as different from each other as possible. The framework is abstract and flexible and can be applied in many different ways. Stretch its limits in some way in at least one of your examples.
# 3. Show for infinite-horizon problems with discount factor ``0 < \\gamma < 1`` the Bellman equations for a determinstic policy ``\\pi`` are given by ``Q_\\pi(s, a) = \\bar r(s, a) + \\gamma \\sum_{s'}P(s'|s, a)Q_\\pi(s', \\pi(s'))``. Start with the definition ``Q_\\pi(S_t, A_t) = \\mathrm{E}\\left[R_{t+1} + \\gamma R_{t+2} + \\gamma^2 R_{t+3} + \\gamma^3R_{t+4} + \\cdots \\right]`` and show all the intermediate steps that lead to the result.

# ╔═╡ 7c4a9aff-c3c1-48ef-8e83-9d5aa7e75b03
begin
    Base.@kwdef struct CliffWalkingParams
        nx::Int = 4
        ny::Int = 12
        states::LinearIndices{2, NTuple{2, Base.OneTo{Int}}} = LinearIndices((4, ny))
        start::CartesianIndex{2} = CartesianIndex(4, 1)
        goal::CartesianIndex{2} = CartesianIndex(4, ny)
        actions::Vector{CartesianIndex{2}} = [CartesianIndex(0, -1),  # left
                                              CartesianIndex(0, 1),   # right
                                              CartesianIndex(-1, 0),  # up
                                              CartesianIndex(1, 0),   # down
                                             ]
    end
    Base.@kwdef mutable struct CliffWalkingEnv <: AbstractEnv
        params::CliffWalkingParams = CliffWalkingParams()
		position::CartesianIndex{2} = params.start
        cumulative_reward::Int = 0
	end
    iscliff(p, nx, ny) = p[1] == nx && p[2] > 1 && p[2] < ny
    iscliff(env::CliffWalkingEnv) = iscliff(env.position, env.params.nx, env.params.ny)
	function act!(env::CliffWalkingEnv, a::Int)
        x, y = Tuple(env.position + env.params.actions[a])
		env.position = CartesianIndex(min(max(x, 1), env.params.nx),
                                      min(max(y, 1), env.params.ny))
	end
    act!(::CliffWalkingEnv, ::Nothing) = nothing
    state(env::CliffWalkingEnv) = env.params.states[env.position]
    reward(env::CliffWalkingEnv) = env.position == env.params.goal ? 0.0 : (iscliff(env) ? -100.0 : -1.0)
    is_terminated(env::CliffWalkingEnv) = env.position == env.params.goal || iscliff(env)
    reset!(env::CliffWalkingEnv) = env.position = env.params.start
    function Plots.plot(env::CliffWalkingEnv)
        m = float.((!iscliff).(CartesianIndices((env.params.nx, env.params.ny)),
                               env.params.nx, env.params.ny))
        m[env.params.nx, env.params.ny] = 2
        heatmap(m;
                yflip = true, legend = false, aspect_ratio = 1,
                c = palette([:black, :yellow, :green]),
            axis = false, ticks = false)
        for i in 0:env.params.nx
            plot!([0.5, env.params.ny + .5], [i + .5, i + .5], c = :black)
        end
        for i in 0:env.params.ny
            plot!([i + .5, i + .5], [0.5, env.params.nx + .5], c = :black)
        end
        scatter!([env.position[2]], [env.position[1]],
                 markersize = 10, c = :red)
    end
end;

# ╔═╡ c692cc6e-dbb5-40e9-aeaa-486b098c3af1
begin
    import ReinforcementLearning.ReinforcementLearningEnvironments: Cross, Nought
    function autoplaying(player1, player2, tictactoe)
        is_terminated(tictactoe) && return false
        player1 == "machine" && player2 == "machine" && return true
        player1 == "machine" && current_player(tictactoe) == Cross() && return true
        player2 == "machine" && current_player(tictactoe) == Nought() && return true
        false
    end
end;

# ╔═╡ 8bd459cb-20bb-483e-a849-e18caae3beef
MLCourse.list_notebooks(@__FILE__)

# ╔═╡ 412d8fcb-8f98-43b6-9235-a4c228317427
MLCourse.footer()

# ╔═╡ d88c1d9b-3396-42a2-8ebd-81851f778602
begin
    struct CounterButtons
        labels
    end

	function Base.show(io::IO, m::MIME"text/html", button::CounterButtons)
		buttons = [@htl("""<input type="button" value="$label"> """)
                   for label in button.labels]
		show(io, m, @htl("""
		<span>
        $buttons
		<script>
		let count = 0
		const span = currentScript.parentElement
		span.value = [count, 0]
		Array.prototype.forEach.call(span.children, function(child, index){
			child.addEventListener("click", (e) => {
				count += 1
				span.value = [count, index + 1]
				span.dispatchEvent(new CustomEvent("input"))
				e.stopPropagation()
			})
		});
		</script>
		</span>
		"""))
	end

	Base.get(button::CounterButtons) = button.labels
	PlutoUI.BuiltinsNotebook.Bonds.initial_value(::CounterButtons) = [0, 0]
	PlutoUI.BuiltinsNotebook.Bonds.possible_values(::CounterButtons) = PlutoUI.BuiltinsNotebook.Bonds.InfinitePossibilities()
	function PlutoUI.BuiltinsNotebook.Bonds.validate_value(::CounterButtons, val)
		val isa Vector{Int} && val[1] >= 0 && val[2] >= 0
	end
end

# ╔═╡ 7cf26d6c-5a67-4bd6-8ef2-34559b53685b
@bind chasse_actions CounterButtons(["open left door", "open right door",
                                     "reset episode", "reset learner"])

# ╔═╡ bf600a9b-484e-4a7b-bf3b-409ebad51cd0
md"""
If you want to see Q-learning in action, choose "qlearner" below and play our
chasse au trésor. The displayed Q-values will then be for the `qlearner`.

Learner $(@bind learner Select(["mclearner", "qlearner"]))

Below you see the Q-values of an epsilon-greedy policy played for ``10^6`` steps.
Note that is comes close to the optimal solution. This is different than what we found with the Monte Carlo Learner. The reason is that the Monte Carlo Learner estimates the Q-values for the epsilon-greedy policy, whereas the Q-Learner estimates the Q-values for the greedy policy, even though the actions are choosen according to the epsilon-greedy policy. The Monte Carlo Learner is called *on-policy* because it evaluates the policy used for action selection. The Q-Learner is called *off-policy*, because it evaluates the greedy policy, while selecting actions with the epsilon-greedy policy.
"""

# ╔═╡ 5a705ba7-dd9c-411c-a910-659fb1ec9f82
md"Below you see the Q-values of the \`$learner\`."

# ╔═╡ 94c0859f-f2d4-4b00-9793-14c823bbc705
@bind tictactoe_action CounterButtons([string.(1:9); "reset game"; "step"])

# ╔═╡ 4a9fb8a0-81fb-4e59-8208-61df1dbd8255
md"""Player Cross: $(@bind player1 Select(["human", "machine"])) Player Nought: $(@bind player2 Select(["human", "machine"]))"""

# ╔═╡ 8c0a5e46-c790-4cec-ac57-1b2813b81358
@bind cw_action CounterButtons(["left", "right", "up", "down", "reset"])

# ╔═╡ 761d690d-5c73-40dd-b38c-5af67ee837c0
begin
    function onehot(i)
        x = zeros(7)
        x[i] = 1
        x
    end
    act(state, action) = wsample(1:7, T[action, state])
    _reward(r::Number) = r
    _reward(r::AbstractArray) = rand(r)
    reward(state, action) = _reward(R[state, action])
    T = reshape([[0, .7, .3, 0, 0, 0, 0], [0, 0, 0, .3, .7, 0, 0],
                  onehot(6), onehot(6),
                  onehot(7), onehot(6),
                  onehot(6), onehot(6),
                  onehot(6), onehot(7),
                  onehot(6), onehot(6),
                  onehot(7), onehot(7)], 2, :)
    R = [-.5 -.3
         [1:4] 1
         -5  1
         1 [3:6]
         1 -5
         0 0
         0 0]
    function showQ(Q; states = ["red room", "blue without guard", "blue with guard", "green without guard", "green with guard", "treasure room", "K.O."], actions = ["left", "right"])
        df = DataFrame(Q, states)
        df.action = actions
        df
    end
    Base.@kwdef mutable struct ChasseAuTresorEnv
        state::Int = 1
        reward::Float64 = 0
        episode_recorder = Tuple{Int, Int, Float64, Int}[]
        action::Union{Int,Nothing} = nothing
    end
    function act!(env::ChasseAuTresorEnv, a)
        length(env.episode_recorder) == 2 && return
        s = env.state
        r = reward(s, a)
        env.reward = r
        s′ = act(env.state, a)
        env.state = s′
        push!(env.episode_recorder, (s, a, r, s′))
    end
    reward(env::ChasseAuTresorEnv) = env.reward
    state(env::ChasseAuTresorEnv) = env.state
    function reset!(env::ChasseAuTresorEnv)
        env.state = 1
        env.reward = 0
        empty!(env.episode_recorder)
    end
    function distributed_representation(s)
        [s == 1; # 1 if in room 1
         s ∈ (2, 3); # 1 if in room 2
         s ∈ (4, 5); # 1 if in room 3
         s == 6; # 1 if in treasure room
         s == 7; # 1 if KO
         s ∈ (3, 5)] # 1 if guard present
    end
#     chasse = ChasseAuTresorEnv()
end;

# ╔═╡ 712c2a9e-4413-4d7a-b729-cfb219723256
let mclearner = MCLearner(na = 2, ns = 7),
    chasse = ChasseAuTresorEnv()
    for _ in 1:10^6
		reset!(chasse)
        episode = []
        for steps in 1:2
            s = state(chasse) # observe the environment
            a = epsilon_greedy_policy(mclearner.Q[:, s]) # use policy to find action
            act!(chasse, a) # act in the environment
            r = reward(chasse) # observe the reward
            push!(episode, (s, a, r)) # store in the episode buffer
        end
        update!(mclearner, episode) # update learner after the end of the episode
    end
    showQ(mclearner.Q)
end

# ╔═╡ bd8557bc-86f2-4ccc-93e9-a6bd843e80be
let qlearner = QLearner(na = 2, ns = 7),
    chasse = ChasseAuTresorEnv()
    for _ in 1:10^6
        reset!(chasse)
        for steps in 1:2
            s = state(chasse) # observe the environment
            a = epsilon_greedy_policy(qlearner.Q[:, s]) # use policy to find action
            act!(chasse, a) # act in the environment
            r = reward(chasse) # observe the reward
            s′ = state(chasse) # observe the next state
            update!(qlearner, s, a, r, s′) # update the learner
        end
    end
    showQ(qlearner.Q)
end

# ╔═╡ ddd02bd9-9577-44a8-8bb6-2b1d11938121
state(tictactoe)

# ╔═╡ aca72848-c488-46d1-a3df-0cf6fb04f4a3
let
	act!(tictactoe, 5) # place a cross in the center
	Text(RLBase.state(tictactoe)) # diplay the state nicely
end

# ╔═╡ 86a78734-2b31-4ee5-8560-8a9388672b45
reset!(tictactoe);

# ╔═╡ 6a96c33a-b6b3-4a0a-83c8-a0df113887d0
mcl = let tictactoe = TicTacToeEnv()
    mcl = MCLearner(na = 9, ns = 5478) # total number of actions and states
    episode_cross = [] # buffer where we store the positions and actions of cross
    episode_nought = [] # buffer where we store the positions and actions of nought
    for game in 1:10^5
        reset!(tictactoe)
        move = 0
        while true
            move += 1
            legal_a = legal_action_space(tictactoe)
            s = state(tictactoe)
            i = epsilon_greedy_policy(mcl.Q[legal_a, s])
            a = legal_a[i]
            act!(tictactoe, a)
            if is_terminated(tictactoe)
                r = reward(tictactoe) # reward of next player
                if isodd(move)
                    r *= -1 # cross finished
                end
                push!(episode_cross, (s, a, r))
                push!(episode_nought, (s, a, -r))
                break
            else
                push!(isodd(move) ? episode_cross : episode_nought,
                      (s, a, 0.))
            end
        end
        update!(mcl, episode_cross)
        update!(mcl, episode_nought)
        empty!(episode_cross)
        empty!(episode_nought)
    end
    mcl
end;

# ╔═╡ 37d0f305-3aea-4e22-ba0b-184fe880f381
all_states = Dict();

# ╔═╡ 36e0ad57-08bd-4ece-81ba-df180d9c476f
begin
    for (k, v) in all_states # housholding
        if v.time < time() + 86400
            pop!(all_states, k)
        end
    end
    @bind(user_id, html"""<script>
            currentScript.value = Math.random()
        </script>""")
end

# ╔═╡ 402d6e7b-fdd1-4082-a995-d78fc7d4cb69
create_initial_state() = (time = time(),
                          chasse = ChasseAuTresorEnv(),
                          tictactoe = TicTacToeEnv(),
                          qlearner = QLearner(),
                          mclearner = MCLearner(),
                          cwenv = CliffWalkingEnv());

# ╔═╡ e98b3e4f-d17e-4fdd-af1c-a8744ce7ecc3
let
    states = get!(all_states, user_id, create_initial_state())
    chasse = states.chasse
    chasse.action = chasse_actions[1] == 0 ? 3 : chasse_actions[2]
    if chasse.state ≤ 5 || (isa(chasse.action, Int) && chasse.action > 2)
        _learner = if learner == "mclearner"
            if length(chasse.episode_recorder) > 0 && chasse.action == 3
                update!(states.mclearner, chasse.episode_recorder)
            end
            states.mclearner
        else
            if length(chasse.episode_recorder) > 0
                update!(states.qlearner, last(chasse.episode_recorder)...)
            end
            states.qlearner
        end
        if chasse.action == 3
            reset!(chasse)
        elseif chasse.action == 4
            reset!(_learner)
        elseif isa(chasse.action, Int)
            act!(chasse, chasse.action)
        end
    end
    d = distributed_representation(chasse.state)
    room = findfirst(d)
    guard = d[6]
    guard_door = room > 2
    gold = d[4]
    room_colors = [:red, :blue, :green, :orange, :black]
    plot(xlim = (0, 1), ylim = (0, 1), size = (600, 400),
         bg = room_colors[room], framestyle = :none, legend = false)
    if room < 4
        plot!([.1, .1, .4, .4], [0, .7, .7, 0], w = 5, c = :black)
        plot!([.14, .12, .17], [.37, .35, .35], w = 4, c = :black)
        plot!(.5 .+ [.1, .1, .4, .4], [0, .7, .7, 0], w = 5, c = :black)
        plot!(.5 .+ [.14, .12, .17], [.37, .35, .35], w = 4, c = :black)
    end
    if guard
        xshift = guard_door * .5
        plot!(xshift .+ [.14, .17, .25, .33, .37], [0, 0, .3, 0, 0], w = 5, c = :black)
        plot!(xshift .+ [.25, .25], [.3, .52], w = 5, c = :black)
        scatter!(xshift .+ [.25], [.58], markersize = 25, markerstrokewidth = 0, c = :black)
        plot!(xshift .+ [.33, .33, .18, .18], [.23, .47, .47, .23],
              w = 5, c = :black)
    end
    if gold
        x = [.25, .15, .2, .3, .35, .25] .+ .25
        y = [.3, .1, .2, .2, .1, .1]
        r = floor(Int, chasse.reward)
        x = x[1:r]
        y = y[1:r]
        scatter!(x, y, markerstrokewidth = 3, c = :yellow, markersize = 28)
    end
    if room == 5
        scatter!([0], [0], c = :black, legend = false, markerstrokewidth = 0) # dummy
        annotate!([(.5, .5, "K.O.", :red)])
    end
    rs = length(chasse.episode_recorder) == 0 ? [0] : getindex.(chasse.episode_recorder, 3)
    annotate!([(.5, .9, "reward = $(chasse.reward)", :white),
               (.5, .8, "cumulative reward = $(join(rs, " + ")) = $(sum(rs))", :white)
              ])
end

# ╔═╡ e876c526-30f9-458d-abf5-e20e6aa0268e
let
    chasse_actions
    states = get!(all_states, user_id, create_initial_state())
    if learner == "mclearner"
        showQ(states.mclearner.Q)
    else
        showQ(states.qlearner.Q)
    end
end

# ╔═╡ ebebd97a-9dc2-4b39-a998-9279d52c57e5
let
    states = get!(all_states, user_id, create_initial_state())
    tictactoe = states.tictactoe
    a = last(tictactoe_action)
    if a in legal_action_space(tictactoe)
        act!(tictactoe, a)
    elseif a == 10
        reset!(tictactoe)
    end
    if autoplaying(player1, player2, tictactoe)
        legal_a = legal_action_space(tictactoe)
        a2 = legal_a[argmax(mcl.Q[legal_a, state(tictactoe)])]
        act!(tictactoe, a2)
    end
    cpl = split(string(current_player(tictactoe)), ".")[end][1:end-2]
    s = if is_terminated(tictactoe)
        if ReinforcementLearning.reward(tictactoe) == -1
            "$cpl lost the game."
        else
            "The game ended in a draw."
        end
    else
        "\nCurrent player: $cpl"
    end
    legal_a = legal_action_space(tictactoe)
    s *= "\n\n" * RLBase.state(tictactoe)
    s *= "\n\n Legal actions:\n\n$(join([join([k in legal_a ? string(k) : " " for k in 3 .* (0:2) .+ j], " ") for j in 1:3], "\n"))"
    s *= "\n\n Game state: $(state(tictactoe))"
    Text(s)
end

# ╔═╡ 61b3c6ca-c680-41d1-9eb5-6ec2f799f0d1
let
    states = get!(all_states, user_id, create_initial_state())
    cwenv = states.cwenv
    a = last(cw_action)
    if a == 5
        r = 0
        cwenv.cumulative_reward = 0
        reset!(cwenv)
    elseif a == 0
        r = 0
    elseif !is_terminated(cwenv)
        act!(cwenv, a)
        r::Int = reward(cwenv)
        cwenv.cumulative_reward += r
    end;
    plot(cwenv)
    plot!(size = (700, 400))
    annotate!([(cwenv.params.ny ÷ 2, -1, "reward = $r"),
               (cwenv.params.ny ÷ 2, -.2, "cumulative reward = $(cwenv.cumulative_reward)")])
end

# ╔═╡ Cell order:
# ╠═b97724e4-d7b0-4085-b88e-eb3c5bcbe441
# ╟─b3793299-c916-40d3-bd87-31153fc3781a
# ╟─ce405f97-6d60-4ae4-b183-79e6c88d9811
# ╟─e98b3e4f-d17e-4fdd-af1c-a8744ce7ecc3
# ╟─7cf26d6c-5a67-4bd6-8ef2-34559b53685b
# ╟─5a705ba7-dd9c-411c-a910-659fb1ec9f82
# ╟─e876c526-30f9-458d-abf5-e20e6aa0268e
# ╟─b0b89b30-38f1-45c1-8c7a-796ea2f41e8d
# ╠═7c756af7-08a0-4878-8f53-9eaa79c16c1a
# ╟─6dced47e-cc0f-4ae5-bdd0-a551c6d5a5b3
# ╠═64863a37-1c25-4d1c-9f2d-33d87b4039a3
# ╟─0537554c-76e4-4827-b92a-bf5b944685c7
# ╠═712c2a9e-4413-4d7a-b729-cfb219723256
# ╟─dc7c5cb0-30b2-4427-8aeb-f312a88effd1
# ╠═e39227c8-8148-43cb-9351-774682b65646
# ╟─bf600a9b-484e-4a7b-bf3b-409ebad51cd0
# ╠═bd8557bc-86f2-4ccc-93e9-a6bd843e80be
# ╟─aff5b0d2-c5e5-4fda-a93c-54a8bca5187f
# ╠═1a17e9b2-a439-4521-842c-96ebe0378919
# ╠═ae72b2c3-fe7b-4ffb-b45c-e589529209c7
# ╠═ddd02bd9-9577-44a8-8bb6-2b1d11938121
# ╟─d094b8d7-5cc7-4bc0-a406-2e3c85f9fb60
# ╠═aca72848-c488-46d1-a3df-0cf6fb04f4a3
# ╠═4b90f08e-0183-4f8b-a6cd-e66b6938f4c8
# ╠═86a78734-2b31-4ee5-8560-8a9388672b45
# ╟─59029032-1c91-4da1-a61b-6a56449dcd2c
# ╟─94c0859f-f2d4-4b00-9793-14c823bbc705
# ╟─4a9fb8a0-81fb-4e59-8208-61df1dbd8255
# ╟─ebebd97a-9dc2-4b39-a998-9279d52c57e5
# ╟─d99a218a-56e8-4081-bd1a-ba7729f529cf
# ╠═6a96c33a-b6b3-4a0a-83c8-a0df113887d0
# ╟─e61c6d43-a097-4f44-a343-05371b4932ef
# ╟─c692cc6e-dbb5-40e9-aeaa-486b098c3af1
# ╟─b6b835d1-8f84-4148-8d5b-c7aea6b0c312
# ╟─8c0a5e46-c790-4cec-ac57-1b2813b81358
# ╟─61b3c6ca-c680-41d1-9eb5-6ec2f799f0d1
# ╟─7c4a9aff-c3c1-48ef-8e83-9d5aa7e75b03
# ╟─8bd459cb-20bb-483e-a849-e18caae3beef
# ╟─412d8fcb-8f98-43b6-9235-a4c228317427
# ╟─d88c1d9b-3396-42a2-8ebd-81851f778602
# ╟─761d690d-5c73-40dd-b38c-5af67ee837c0
# ╟─37d0f305-3aea-4e22-ba0b-184fe880f381
# ╟─36e0ad57-08bd-4ece-81ba-df180d9c476f
# ╟─402d6e7b-fdd1-4082-a995-d78fc7d4cb69

max_stars_repo_name: jbrea/MLCourse

max_issues_repo_name: davibarreira/MLCourse

max_forks_repo_name: davibarreira/MLCourse


text:
function reparam_CIR(δ, γ, σ)
    a = 2*γ
    b = δ*σ^2/(2*γ)
    σ_prime = 2*σ
    return a, b, σ_prime
end

function inverse_reparam_CIR(a, b, σ_prime)
    γ = a/2
    δ = 4*a*b/σ_prime^2
    σ = σ_prime/2
    return δ, γ, σ
end

function joint_sampler_CIR_reparam_pruning_precompute(data, λ, prior_logpdf, niter, do_the_pruning::Function; final_chain_length = 1000, silence = false, jump_size = 0.5*Matrix(I,3,3), θ_init = [1.,1.,1.])

    joint_sampler_CIR_pruning_precompute(data, λ, prior_logpdf, θ_init, niter, do_the_pruning::Function; final_chain_length = final_chain_length, silence = silence, jump_size = jump_size, reparam = true)
end

function joint_sampler_CIR_reparam_keep_fixed_number_precompute(data, λ, prior_logpdf, niter, fixed_number::Int64; final_chain_length = 1000, silence = false, jump_size = 0.5*Matrix(I,3,3), θ_init = [1.,1.,1.])
    # println("filter_WF_mem2")

    function prune_keeping_fixed_number(Λ_of_t, wms_of_t)
        Λ_of_t_kept, wms_of_t_kept = keep_fixed_number_of_weights(Λ_of_t, wms_of_t, fixed_number)
        return Λ_of_t_kept, normalise(wms_of_t_kept)
    end

    joint_sampler_CIR_reparam_pruning_precompute(data, λ, prior_logpdf, niter, prune_keeping_fixed_number; final_chain_length = 1000, silence = silence, jump_size = jump_size, θ_init = θ_init)

end

function reparam_joint_loglikelihood_CIR(data, trajectory, times, a, b, σ_prime, λ)
    δ, γ, σ = inverse_reparam_CIR(a, b, σ_prime)
    return joint_loglikelihood_CIR(data, trajectory, times, δ, γ, σ, λ)
end

ext: jl

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
# permutations(a) = Combinatorics.permutations(a)

Combinatorics.permutations(a, mode::Type{Val{:unique}}) = permutations(a, length(a), mode)
Combinatorics.permutations(a, mode::Type{Val{:repeated}}) = permutations(a, length(a), mode)
Combinatorics.permutations(a::T, k::Integer, ::Type{Val{:unique}}) where {T} = Permutations{T, Val{:unique}}(a, k)
Combinatorics.permutations(a::T, k::Integer, ::Type{Val{:repeated}}) where {T} = Permutations{T, Val{:repeated}}(a, k)

struct Permutations{T,M}
    a::T
    k::Int
end

Base.start(v::Permutations{T, Val{:unique}}) where {T} = v.k <= length(v.a) ? collect(1:v.k) : [length(v.a) + 1]
Base.start(v::Permutations{T, Val{:repeated}}) where {T} = ones(Int, v.k)

Base.next(v::Permutations{T,M}, s) where {T,M} = begin
    permutation = [v.a[si] for si in s]
    v.k > 0 || return permutation, [length(v.a) + 1]

    s = copy(s)
    for i in length(s):-1:1
        next!(v, s, i)

        if s[i] <= length(v.a)
            s[i+1:end] = tail(v, s, i)
            break
        end
    end

    permutation, s
end

next!(::Permutations{T, Val{:unique}}, s, i) where {T} = begin
    s[i] += 1
    while s[i] ∈ s[1:i-1]
        s[i] += 1
    end
end
next!(::Permutations{T, Val{:repeated}}, s, i) where {T} = s[i] += 1

tail(v::Permutations{T, Val{:unique}}, s, i) where {T} = setdiff(collect(1:v.k), s[1:i])[1:(length(s)-i)]
tail(::Permutations{T, Val{:repeated}}, s, i) where {T} = ones(Int, length(s) - i)

Base.done(v::Permutations, s) = !isempty(s) && s[1] > length(v.a)

Base.eltype(v::Permutations) = eltype(typeof(v))
Base.eltype(::Type{Permutations{T,M}}) where {T,M} = Array{eltype(T),1}

Base.length(v::Permutations{T, Val{:unique}}) where {T} = begin
    v.k <= length(v.a) || return big(0)
    factorial(big(length(v.a)), length(v.a) - v.k)
end
Base.length{T}(v::Permutations{T, Val{:repeated}}) = big(length(v.a))^v.k


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Hanabi"
version = v"0.1.1"

# Collection of sources required to build HanabiBuilder
sources = [
    "https://github.com/findmyway/hanabi-learning-environment.git" =>
    "4751515fe07404af5119927acdabb790e35147c2",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd hanabi-learning-environment/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain .
make
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:i686, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpyhanabi", :libpyhanabi)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


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

using Markdown
using InteractiveUtils

# ╔═╡ eb807d0e-1dfd-11eb-155d-df6f207266b9
using Pkg, DrWatson

# ╔═╡ fa048db0-1de8-11eb-20cd-9b8d6a6fd1f1
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, LinearAlgebra, Plots
	using LaTeXStrings, Measures
end

# ╔═╡ f9c7005e-1dfd-11eb-08de-4d71184b2337
md"## Listing 1.14"

# ╔═╡ 9da9d70a-1dfe-11eb-054f-ffe619ac258d
Random.seed!(1974); rand()

# ╔═╡ 9daa179c-1dfe-11eb-206d-4f94676a8104
Random.seed!(1975); rand()

# ╔═╡ 9daaa4f0-1dfe-11eb-1bb8-4f632207e49b
Random.seed!(1974); rand()

# ╔═╡ 2b564d56-1de9-11eb-120b-7137b81725dd
md"## Listing 1.15"

# ╔═╡ e9553c4a-1dfd-11eb-28d8-f3f4c0297d2a
Random.seed!()

# ╔═╡ b1edef30-1dfe-11eb-26fd-9782f1f70d10
begin
	N = 10^4
	data1     = [[rand(),rand()] for _ in 1:N]
	indata   = filter((x)-> (norm(x) <= 1), data1)
	outdata  = filter((x)-> (norm(x) > 1), data1)
	piApprox = 4*length(indata)/N
	md"Pi Estimate: $(piApprox)"
end

# ╔═╡ b1ee272a-1dfe-11eb-1956-91da67abf143
begin
	scatter(first.(indata),last.(indata), c=:blue, ms=1, msw=0)
	scatter!(first.(outdata),last.(outdata), c=:red, ms=1, msw=0,
		xlims=(0,1), ylims=(0,1), legend=:none, ratio=:equal)
end

# ╔═╡ 2e3fb9ba-1de9-11eb-0835-db03a7547b7c
md"## Listing 1.16"

# ╔═╡ 3c2312fe-1dfd-11eb-1afe-11077862ebb0
begin
	a, c, m = 69069, 1, 2^32
	next(z) = (a*z + c) % m

	#N = 10^6
	data2 = Array{Float64,1}(undef, N)

	x = 808

	for i in 1:N
    	data2[i] = x/m
    	global x = next(x)
	end
end

# ╔═╡ e729992c-1e00-11eb-13dc-0b9f2d96ef0d
begin
	p1 = scatter(1:1000, data2[1:1000], xlim=(0, 1000), ylim=(0,1),
		c=:blue, m=4, msw=0, xlabel=L"n", ylabel=L"x_n")
	p2 = histogram(data2, xlim=(0,1.05), bins=50, normed=:true,
		ylims=(0,1.1), xlabel="Support", ylabel="Density")
	plot(p1, p2, size=(800, 400), legend=:none, margin = 5mm)
end

# ╔═╡ 082010d8-1de9-11eb-23df-3fd3c98c4c1d
md"## Listing 1.17"

# ╔═╡ 4441a4ba-1e01-11eb-16d8-67b2c22f7138
function path(rng, alpha, n=5000)
    x, y = 0.0, 0.0
    xDat, yDat = [], []
    for _ in 1:n
        flip = rand(rng,1:4)
        if flip == 1
            x += 1
        elseif flip == 2
            y += 1
        elseif flip == 3
            x -= (2+alpha)*rand(rng)
        elseif flip == 4
            y -= (2+alpha)*rand(rng)
        end
        push!(xDat,x)
        push!(yDat,y)
    end
    return xDat, yDat
end

# ╔═╡ 4441d4e4-1e01-11eb-201e-91694cb132aa
alphaRange = [0.2, 0.21, 0.22]

# ╔═╡ 44428c8e-1e01-11eb-0cb0-fb6fd51c9e04
default(xlabel = "x", ylabel = "y", xlims=(-150,50), ylims=(-250,50))

# ╔═╡ 4454e5fc-1e01-11eb-3d1f-bd6396d87471
begin
	p3 = plot(path(MersenneTwister(27), alphaRange[1]), c=:blue)
	plot!(path(MersenneTwister(27), alphaRange[2]), c=:red)
	plot!(path(MersenneTwister(27), alphaRange[3]), c=:green)
end;

# ╔═╡ 447a4342-1e01-11eb-22ea-a3ceff1f0df1
rng = MersenneTwister(27)

# ╔═╡ eaf93996-1e04-11eb-00d7-95f753122789
begin
	p4 = plot(path(rng, alphaRange[1]), c=:blue)
	plot!(path(rng, alphaRange[2]), c=:red)
	plot!(path(rng, alphaRange[3]), c=:green)
end;

# ╔═╡ 44a7eb76-1e01-11eb-134f-1d67bf11c36e
plot(p3, p4, size=(500, 400), legend=:none, margin=5mm)

# ╔═╡ ded51f3e-1e01-11eb-0717-af179845759e
md"## End of listing 1.14"

# ╔═╡ Cell order:
# ╟─f9c7005e-1dfd-11eb-08de-4d71184b2337
# ╠═eb807d0e-1dfd-11eb-155d-df6f207266b9
# ╠═fa048db0-1de8-11eb-20cd-9b8d6a6fd1f1
# ╠═9da9d70a-1dfe-11eb-054f-ffe619ac258d
# ╠═9daa179c-1dfe-11eb-206d-4f94676a8104
# ╠═9daaa4f0-1dfe-11eb-1bb8-4f632207e49b
# ╟─2b564d56-1de9-11eb-120b-7137b81725dd
# ╠═e9553c4a-1dfd-11eb-28d8-f3f4c0297d2a
# ╠═b1edef30-1dfe-11eb-26fd-9782f1f70d10
# ╠═b1ee272a-1dfe-11eb-1956-91da67abf143
# ╟─2e3fb9ba-1de9-11eb-0835-db03a7547b7c
# ╠═3c2312fe-1dfd-11eb-1afe-11077862ebb0
# ╠═e729992c-1e00-11eb-13dc-0b9f2d96ef0d
# ╟─082010d8-1de9-11eb-23df-3fd3c98c4c1d
# ╠═4441a4ba-1e01-11eb-16d8-67b2c22f7138
# ╠═4441d4e4-1e01-11eb-201e-91694cb132aa
# ╠═44428c8e-1e01-11eb-0cb0-fb6fd51c9e04
# ╠═4454e5fc-1e01-11eb-3d1f-bd6396d87471
# ╠═447a4342-1e01-11eb-22ea-a3ceff1f0df1
# ╠═eaf93996-1e04-11eb-00d7-95f753122789
# ╠═44a7eb76-1e01-11eb-134f-1d67bf11c36e
# ╟─ded51f3e-1e01-11eb-0717-af179845759e


text:
using Distributions
using ConjugatePriors: posterior_canon, NormalInverseChisq
using Particles: marginal_lhood, marginal_log_lhood

@testset "Normal-χ^-2" begin
    x = rand(Normal(3, 5), 100)
    ss = suffstats(Normal, x)
    prior = NormalInverseChisq(3.0, 5.0, 2., 3.)
    posterior = posterior_canon(prior, ss)

    @testset "Posterior predictive" begin
        pp = posterior_predictive(posterior)

        xpp = rand(pp, 1_000_000)
        xpp2 = [rand(Normal(μ, sqrt(σ2))) for (μ, σ2) in (rand(posterior) for _ in 1:1_000_000)]

        @test isapprox(mean(xpp), mean(xpp2), rtol=0.01)
        @test isapprox(std(xpp), std(xpp2), rtol=0.01)
    end

    @testset "marginal likelihood" begin
        @test marginal_log_lhood(prior, ss) ≈ log(marginal_lhood(prior, ss))
        mean_logpdf = log(mean(exp(sum(logpdf.(Normal(μ, sqrt(σ2)), x)))
                               for (μ, σ2)
                               in (rand(prior) for _ in 1:1000)))#0_000)]
        @test isapprox(marginal_log_lhood(prior, ss), mean_logpdf, rtol=0.01)
    end

end


text:
push!(LOAD_PATH, "..")

using
    Documenter,
    Literate,
    Oceananigans,
    Oceananigans.Operators,
    Oceananigans.TurbulenceClosures

####
#### Generate examples
####

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

examples = [
    "simple_diffusion.jl",
    "two_dimensional_turbulence.jl",
    "ocean_wind_mixing_and_convection.jl",
    "ocean_convection_with_plankton.jl",
    "internal_wave.jl"
]

for example in examples
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR, documenter=true)
end

####
#### Build docs
####

makedocs(
   modules = [Oceananigans, Oceananigans.Operators, Oceananigans.TurbulenceClosures],
   doctest = true,
   clean   = true,
 checkdocs = :all,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                             canonical = "https://climate-machine.github.io/Oceananigans.jl/latest/"),
   authors = "Ali Ramadhan, Gregory Wagner, John Marshall, Jean-Michel Campin, Chris Hill",
  sitename = "Oceananigans.jl",
     pages = [
         "Home"   => "index.md",
         "Manual" => [
             #"Model setup" => "manual/model_setup.md",
             "Examples" => [
                 "One-dimensional diffusion"        => "generated/simple_diffusion.md",
                 "Two-dimensional turbulence"       => "generated/two_dimensional_turbulence.md",
                 "Ocean wind mixing and convection" => "generated/ocean_wind_mixing_and_convection.md",
                 "Ocean convection with plankton"   => "generated/ocean_convection_with_plankton.md",
                 "Internal wave"                    => "generated/internal_wave.md"
             ],
             "Physics" => "manual/physics.md",
             "Numerical implementation" => [
                 #"Overview"               => "manual/overview.md",
                 "Pressure decomposition" => "manual/pressure_decomposition.md",
                 "Time stepping"          => "manual/time_stepping.md",
                 "Finite volume method"   => "manual/finite_volume.md",
                 "Spatial operators"      => "manual/spatial_operators.md",
                 "Boundary conditions"    => "manual/boundary_conditions.md",
                 "Poisson solvers"        => "manual/poisson_solvers.md",
                 "Large eddy simulation"  => "manual/large_eddy_simulation.md"
             ],
             "Verification" => [
                 "Taylor-Green vortex"     => "verification/taylor_green_vortex.md",
                 "Stratified Couette flow" => "verification/stratified_couette_flow.md"
             ],
             #"Appendix" => [
                 #"Staggered grid"         => "manual/staggered_grid.md",
                 #"Fractional step method" => "manual/fractional_step.md",
             #],
         ],
         "Gallery"    => "gallery.md",
         "Benchmarks" => "benchmarks.md",
         "Library"    => "library.md",
         "Index"      => "subject_index.md"
     ]
)

deploydocs(repo = "github.com/climate-machine/Oceananigans.jl.git")

####
#### Delete leftover JLD2 files.
#### See: https://github.com/climate-machine/Oceananigans.jl/issues/509
####

const GENERATED_DIR = joinpath(@__DIR__, "build/generated")

@info "Deleting leftover JLD2 files..."
leftovers = filter(x -> occursin(".jld2", x), readdir(GENERATED_DIR))
for fname in leftovers
    fpath = joinpath(GENERATED_DIR, fname)
    rm(fpath, force=true)
    @info "Deleted: $fpath"
end


text:
using RetroSignalModel
using Test

@testset "Data files" begin
    include("load_data.jl")
end

@testset "Build models" begin
    include("build_model.jl")
end

@testset "Parameter searching" begin
    include("params.jl")
end

@testset "Steady state searching" begin
    include("steadystates.jl")
end


text:
# --- L2 Cost

type L2Cost <: Layer
    output
    grad_input
    L2Cost() = new()
end
export L2Cost

# Defines forward & backward
forward!(cost::L2Cost, inputs) = cost.output = .5 * norm(inputs[1] - inputs[2])
backward!(cost::L2Cost, inputs, scale=1.) = cost.grad_input = scale * (inputs[1] - inputs[2])


text:
# Primer8.2
include("labyrinth.jl")
include("primer82BFS.jl")

nbRows = 8;
nbColumns = 8;
# generisanje lavirinta
G = emptyLabyrinth(nbRows, nbColumns)

# definisanje tipa za odredjeni podskup cvorova iz grafa
defineNodesType!(G, 1:nbRows*nbColumns, 'P'); #put
defineNodesType!(G, [5 9 10 11 13 15 21 23 24 26 28 29 31 34 44 45 46 47 50 55 58 60 61 63], 'Z'); #zid
defineNodesType!(G, [62], 'C'); #cilj

println("Izgled grafa pre obilaska")
printLabyrinthType(G)
println()

path = BFS!(G, 57)
println("Ispis grafa posle obilaska")
printLabyrinthColor(G)
println()

println("Ispis putanje i izgled lavirinta posle obilaska")
println("PATH = $path")
printLabyrinthPath(G, path)


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

using Markdown
using InteractiveUtils

# ╔═╡ 8c925150-4c23-409e-a5a6-b275cd8b1e2c
using Pkg, DrWatson

# ╔═╡ ac8218e1-f8a9-476e-8622-59f4fd2748c2
begin
	@quickactivate "NumericalMethodsforEngineers"
	using Symbolics, SymbolicUtils, ModelingToolkit
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ a52a571a-d590-4540-af19-b57039203d65
begin
	@parameters t
	@variables x y
	@variables γ[1:2]
end;

# ╔═╡ 11d9e598-990e-484f-ba9a-fbe015a4050c
function f(x, y)
    @variables z
    z = x^2 + y
    [x - z, x^2 - y, z + y]
end

# ╔═╡ c3d7eccd-2437-4c4c-b3e2-7be634c68ce2
L.(f(x, y))

# ╔═╡ cd32f519-6a7b-4213-bb77-2e86a51cf55b
v = [4, 4]

# ╔═╡ 534d69c7-7f10-4b02-a589-511a13e8867a
f(v...)

# ╔═╡ 2192bb67-7e70-42aa-bea2-cf1350d4e417
sum(f(v...))

# ╔═╡ 76e59dd6-0d53-4633-aa9d-0536518f6ef1
L.(f(γ...))

# ╔═╡ 10b84a16-9b8f-42d3-988a-7af9f27e0697
SymbolicUtils.istree(f)

# ╔═╡ 86458fb9-3d59-410b-b5b7-2bd7e33e9ccb
begin
	@variables z1
	z1 = x^2 + y
	L(z1)
end

# ╔═╡ 3969d188-b5a0-49ce-8745-0c6c4aa45994
begin
	ex2 = [x - z1, x^2 - y, z1 + y]
	L(ex2)
end

# ╔═╡ bd8e2cab-63f1-4c62-8cdb-f35d1242f789
testf = build_function(ex2, [x, y, z1], expression=Val{false});

# ╔═╡ a74b2190-63e3-4ab8-a8bb-7503b300f597
L(testf[1]([4, 4, z1]))

# ╔═╡ Cell order:
# ╠═8c925150-4c23-409e-a5a6-b275cd8b1e2c
# ╠═ac8218e1-f8a9-476e-8622-59f4fd2748c2
# ╠═a52a571a-d590-4540-af19-b57039203d65
# ╠═11d9e598-990e-484f-ba9a-fbe015a4050c
# ╠═c3d7eccd-2437-4c4c-b3e2-7be634c68ce2
# ╠═cd32f519-6a7b-4213-bb77-2e86a51cf55b
# ╠═534d69c7-7f10-4b02-a589-511a13e8867a
# ╠═2192bb67-7e70-42aa-bea2-cf1350d4e417
# ╠═76e59dd6-0d53-4633-aa9d-0536518f6ef1
# ╠═10b84a16-9b8f-42d3-988a-7af9f27e0697
# ╠═86458fb9-3d59-410b-b5b7-2bd7e33e9ccb
# ╠═3969d188-b5a0-49ce-8745-0c6c4aa45994
# ╠═bd8e2cab-63f1-4c62-8cdb-f35d1242f789
# ╠═a74b2190-63e3-4ab8-a8bb-7503b300f597


text:
using Photodynamics
using Test

include("common.jl")

@testset "Photodynamics.jl" begin
    print("Impact Parameter... ")
    @testset "Impact Parameter" begin
        include("test_impact_parameter.jl")
    end
    println("Done.")

    print("Lightcurve... ")
    @testset "Lightcurve" begin
        include("test_simpson.jl")
        include("test_integrate_timestep.jl")
        include("test_compute_flux.jl")
        include("test_compute_lightcurve.jl")
    end
end


text:
using Test
using CSV
using DataFrames
using Random

include("../src/pions.jl")
include("../src/proton.jl")

using .ProtonFunctions
using .PionFunctions

data = CSV.read("../src/data/p-gamma_cross_section.csv", DataFrame)
ϵ = data[!,1]*1e3 # photon energy in MeV
σ = data[!,2]*1e-13 # cross-section in nanometer^2

mp = 938e6 # eV

@testset "Testing interaction_rate function" begin
    # interaction_rate below the threshold (5e19 eV) should be zero
    @test interaction_rate((1e16/mp), 0, ϵ, σ) ≈ 0
    @test interaction_rate((1e17/mp), 0, ϵ, σ) ≈ 0
    @test interaction_rate((1e18/mp), 0, ϵ, σ) ≈ 0
    @test interaction_rate((1e19/mp), 0, ϵ, σ) ≈ 0
end;

@testset "Testing adiabatic_energy_loss function" begin
    # after propagating to Earth, E = E/(1+z)
    @test adiabatic_energy_loss(1e10, 1, 0) ≈ (1e10)/2
    @test adiabatic_energy_loss(1e10, 4, 0) ≈ (1e10)/5
end;

@testset "Testing survival_probability function" begin
    # if there is no interaction the survival probability should be 1
    @test survival_probability(1, 0, 1, 0, 0) ≈ 1
    # for a high rate the survival probability should be 0
    @test survival_probability(1, 1, 1, 1, 0) ≈ 0
end;

# Testing assertions:

rng = MersenneTwister()

@testset "Testing neutral_pion_production function" begin
    # low energy protons cannot create pions
    @test_throws AssertionError neutral_pion_production((1e10/mp), 0, 0, rng)
    @test_throws AssertionError neutral_pion_production((1e11/mp), 0, 0, rng)
    @test_throws AssertionError neutral_pion_production((1e12/mp), 0, 0, rng)
end;

@testset "Testing charged_pion_production function" begin
    # low energy protons cannot create pions
    @test_throws AssertionError charged_pion_production((1e10/mp), 0, 0, rng)
    @test_throws AssertionError charged_pion_production((1e11/mp), 0, 0, rng)
    @test_throws AssertionError charged_pion_production((1e12/mp), 0, 0, rng)
end;


text:
using Distributed

# Apply a sum of matrices to a vector in
# parallel using multi-processing.
# The vector exists on every process, but
# each matrix in the sum only exists on
# a single process.

# Set the number of processors Julia will use
# @everywhere denotes that it will be defined
# on every worker.
@everywhere nprocs = 2
addprocs(nprocs)

# Define the vector dimension.
# @everywhere denotes that it will be defined
# on every worker.
@everywhere n = 4

# Define the vector we will apply the set
# of matrices to.
# @everywhere denotes that it will be defined
# on every worker.
@everywhere v = randn(n)

# Define the matrices in the matrix sum
# that will be applied in parallel to the
# vector `v`.
# The first matrix is defined on worker 1,
# the second one on worker 2.
A1 = @spawnat workers()[1] randn(n, n)
A2 = @spawnat workers()[2] randn(n, n)

# Define the individual matrix-vector multiplication
# operations on the appropriate worker.
# Alternatively, you can use:
#
# A1v = @spawnat :any fetch(A1) * v
#
# for Julia to choose the worker automatically.
A1v = @spawnat workers()[1] fetch(A1) * v
A2v = @spawnat workers()[2] fetch(A2) * v

# Actually perform the matrix-vector multiplications,
# fetching the results to the main worker for
# local usage.
Av = fetch(A1v) + fetch(A2v)

@show norm((fetch(A1) + fetch(A2)) * v - Av)

max_stars_repo_name: mtfishman/ITensorParallel


text:
using Test
using CompScienceMeshes

sp = CompScienceMeshes.SphereChart(point(1,2,3), 1.0)
for L in 0:10
    @eval @test sum(w for (p,w) in quadpoints(sp, $L)) ≈ 4π
end

function Ylm(K, r)

    @assert K >= 2 "Multipole order needs to be equal or larger than 2."

    r = r / norm(r)
    x = r[1]; y = r[2]; z = r[3]

    T = eltype(r)
    Y = zeros(Complex{T}, K*K)

    Y[1] = 1/√(4π)

    facphip =  x + im*y;
    facphim = -x + im*y;
    for m ∈ 1:(K-1)
        Y[(m+1)*(m+1)] = -Y[m*m] * facphip * sqrt((m+0.5)/m)
        Y[m*m+1] = -Y[(m-1)*(m-1)+1] * facphim * sqrt((m+0.5)/m)
    end

    Y[3] = √3 * z * Y[1]
    for l ∈ 2:(K-1)
        Y[(l+1)*(l+1)-1] = sqrt(2l+1) * z * Y[l*l]
        Y[l*l+2] = sqrt(2*l+1) * z * Y[(l-1)*(l-1)+1]
    end

    for l ∈ 2:(K-1)
        for m ∈ (2-l):(l-2)
            Y[l*(l+1)+m+1] =
                sqrt((4*l^2-1)/(l^2-m^2)) * z * Y[l*(l-1)+m+1] -
                sqrt(((2*l+1)*((l-1)^2-m^2))/((l^2-m^2)*(2*l-3))) * Y[(l-1)*(l-2)+m+1]
        end
    end

    Y
end

K = 3
I = zeros(K^2)
for (_p,w) in quadpoints(sp, 20)
    Y = Ylm(K, cartesian(_p))
    for l in 0:K-1
        for _m in -l:l
            i = l*(l+1)+_m+1
            I[i] += w * real(Y[i] * conj(Y[i]))
        end
    end
end

@test all(isapprox.(I, 1.0))


text:
VERSION >= v"0.4" && __precompile__()

module JPLEphemeris

using Compat
import Compat: ASCIIString

include("naif.jl")
include("daf.jl")
include("spk.jl")

end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module Symbolic_CE

export symbolic_ce, configure

using ExprSearch, SymbolicRegression, DerivTreeVis
using ExprSearch.CE
using RLESUtils, LogSystems, Loggers, Configure
using DataFrames

import Configure.configure

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

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

"""
Example usage:
symbolic_ce(; seed=1)
"""
function symbolic_ce(;outdir::AbstractString=joinpath(RESULTDIR, "Symbolic_CE"),
                     seed=1,
                     logfileroot::AbstractString="symbolic_ce_log",

                     num_samples::Int64=2000,
                     iterations::Int64=50,
                     elite_frac::Float64=0.6,
                     w_new::Float64=0.4,
                     w_prior::Float64=0.1,
                     maxdepth::Int64=10,
                     randchannel_width::Int64=60,
                     default_code::Any=0.0,

                     ver::Symbol=:easy,

                     vis::Bool=true,
                     vis_type::Symbol=:TEX,

                     deepdive_log::Bool=false
                     )
    srand(seed)
    mkpath(outdir)

    logsys = CE.logsystem()
    empty_listeners!(logsys)
    send_to!(STDOUT, logsys, ["verbose1", "current_best_print", "result"])
    logs = TaggedDFLogger()
    send_to!(logs, logsys, ["code", "computeinfo", "current_best", "elapsed_cpu_s", "fitness",
        "parameters", "result"])

    if deepdive_log
        dlog = TaggedDFLogger()
        send_to!(logsys, ["pcfg_probs"]) do x
            iter, probs = x
            for (s, ps) in probs
                s = string(s)
                if !haskey(dlog, s)
                    add_folder!(dlog, s, fill(Float64,length(ps)))
                end
                push!(dlog, s, ps)
            end
        end
        add_folder!(dlog, "elite_exprs", [Int64, String], ["iteration", "expr"])
        send_to!(logsys, ["elite_exprs"]) do x
            iter, exprs = x
            for ex in exprs
                push!(dlog, "elite_exprs", [iter, ex])
            end
        end
    end

    problem = Symbolic(ver)

    ce_params = CEESParams(num_samples, iterations, elite_frac, w_new,
        w_prior, maxdepth, default_code, randchannel_width, logsys)

    result = exprsearch(ce_params, problem)

    #manually push! extra info to log
    push!(logs, "parameters", ["seed", seed])
    push!(logs, "parameters", ["version", ver])

    outfile = joinpath(outdir, logfileroot)
    save_log(LogFile(outfile), logs)

    if deepdive_log
        outfile = joinpath(outdir, "$(logfileroot)_deepdive")
        save_log(LogFile(outfile), dlog)
    end

    if vis
        derivtreevis(get_derivtree(result), joinpath(outdir, "$(logfileroot)_derivtreevis");
            format=vis_type)
    end
    @show result.expr
    return result
end

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
module Output

using Sockets

using Flight.Attitude
using Flight.Geodesy
using Flight.Kinematics

export XPInterface, disable_physics, set_position

Base.@kwdef struct XPInterface
    socket::UDPSocket = UDPSocket()
    host::IPv4 = IPv4("127.0.0.1")
    port::Integer = 49009
end

function set_dref(xp::XPInterface, dref_id::AbstractString, dref_value::Real)

    #ascii() ensures ASCII data, codeunits returns a CodeUnits object, which
    #behaves similarly to a byte array. this is equivalent to b"text".
    #Vector{UInt8}(dref_id) would also work
    buffer = IOBuffer()
    write(buffer,
        b"DREF\0",
        dref_id |> length |> UInt8,
        dref_id |> ascii |> codeunits,
        UInt8(1),
        Float32(dref_value))

    send(xp.socket, xp.host, xp.port, buffer.data)
end

function set_dref(xp::XPInterface, dref_id::AbstractString, dref_value::AbstractVector{<:Real})

    buffer = IOBuffer()
    write(buffer,
        b"DREF\0",
        dref_id |> length |> UInt8,
        dref_id |> ascii |> codeunits,
        dref_value |> length |> UInt8,
        Vector{Float32}(dref_value))

    send(xp.socket, xp.host, xp.port, buffer.data)
end

disable_physics(xp::XPInterface) = set_dref(xp, "sim/operation/override/override_planepath", 1)

function set_position(xp::XPInterface, init::KinInit, aircraft::Integer = 0)
    set_position(xp, KinData(init).pos, aircraft)
end

function set_position(xp::XPInterface, pos::PosData, aircraft::Integer = 0)

    llh = Geographic(pos.ϕ_λ, pos.h_o)
    euler = REuler(pos.q_nb)

    lat = rad2deg(llh.l2d.ϕ)
    lon = rad2deg(llh.l2d.λ)
    alt = llh.alt
    psi = rad2deg(euler.ψ)
    theta = rad2deg(euler.θ)
    phi = rad2deg(euler.φ)
    set_position(xp; lat, lon, alt, psi, theta, phi, aircraft)

end

function set_position(xp::XPInterface; lat = -998, lon = -998, alt = -998,
                        psi = -998, theta = -998, phi = -998,
                        aircraft::Integer = 0)

    #all angles must be in degrees
    buffer = IOBuffer()
    write(buffer,
        b"POSI\0", UInt8(aircraft),
        Float64(lat), Float64(lon), Float64(alt),
        Float32(theta), Float32(phi), Float32(psi),
        Float32(-998)) #last one is landing gear (?!)

    send(xp.socket, xp.host, xp.port, buffer.data)

end

end #module


text:
using BenchmarkTools

# Originally from Teo ShaoWei via Zulip

const INT = UInt32
const cups = parse.(INT, collect("389125467"))

@inline function peek(next, at, n::Int, result=similar(next,n))
  for i in 1:n
    result[i] = next[at]
    at = next[at]
  end
  result
end

@inline function peek(next, at, ::Val{3})
    a = next[at]
    b = next[ a ]
    c = next[ b ]
    ( a, b, c )
end

function next1(cups)
  N = length(cups)
  next = similar(cups)
  for i in 1:N
    next[cups[i]] = cups[mod1(i+1,N)]
  end
  next
end
function next2(cups)
  next = similar(cups)
  next[cups] = circshift(cups,-1)
  next
end

const cups_next = next2(cups)

function run(current::INT, next::Vector{INT}, steps::Int=1)
  N = INT(length(next))

  @inbounds for i in 1:steps
    #pickups = peek(next, current, Val(3))
    a = next[ current ]
    b = next[ a ]
    c = next[ b ]

    #dst = mod1(current-1, N)
    #while dst in pickups
    #  dst = mod1(dst-1, N)
    #end
    dst = current - INT(1)
    if dst == INT(0)
        dst = N
    end
    while dst == a || dst == b || dst == c
        dst -= INT(1)
        if dst == INT(0)
            dst = N
        end
    end

    next[current] = next[c]
    next[c] = next[dst]
    next[dst] = a
    current = next[current]
  end

  return next
end

part1() = join(peek(run(cups[1], cups_next, 10), 1, 8))
#part2() = prod(peek(run(vcat(cups, 10:1_000_000), 10_000_000), 1, 2))
function part2()
    #cups2 = vcat(cups, Int32(10):Int32(1_000_000))
    next = vcat(cups_next, INT(11):INT(1_000_000), cups[1])
    next[cups[end]] = 10
    # cups2_next == next2(cups2) == next1(cups2)
    next = run(cups[1], next, 10_000_000)
    prod( peek( next , 1 , 2) )
end
println( part2() )
@btime part2() #   171.104 ms (30 allocations: 3.82 MiB)
# 149245887792


text:
function xyz_to_geometry(xyzfile)
  input_file = []
  open(joinpath(@__DIR__, xyzfile)) do file
    input_file = readlines(file)
  end

  natoms = parse(Int, input_file[1])
  #println(natoms)
  geometry = input_file[3:end]
  #display(geometry); println()

  coords_vector = Vector{Float64}([])
  symbols = Vector{String}([])
  for geometry_line in geometry
    coords_match = eachmatch(r"([-]?[\d]{1,}\.[\d]{1,})",geometry_line)
    atom_coords = collect(coords_match)
    atom_coords_match = map(x -> x.match, atom_coords)
    append!(coords_vector, parse.(Float64,atom_coords_match))

    symbol_match = match(r"[A-Z]{1,2}",geometry_line)
    push!(symbols, symbol_match.match)
  end
  #coords = reshape(coords_vector, (3,length(geometry)))
  #coords = transpose(coords)
  coords = coords_vector

  #@assert natoms == size(coords)[1]
  #@assert natoms == length(symbols)

  #display(coords); println()
  #display(symbols); println()

  charge = 0

  return coords, symbols, charge
end

function create_input(input)
  #== create json file if it does not exist ==#
  filename = input[1:(end-4)] * ".json"

  #== write input json file ==#
  output = open(joinpath(@__DIR__, filename), write=true) do file
    write(file, "{\n")
    #== write in molecule information ==#
    write(file, "  \"molecule\": {\n")
    coords, symbols, charge = xyz_to_geometry(input)
    #== write fragment geometry ==#
    write(file, "    \"geometry\" : [\n")
    for icoord in 1:length(coords)
      write(file, "       ")

      value = coords[icoord]
      #println(value)
      if icoord == length(coords)
        write(file, "$value\n")
      elseif icoord%3 == 0
        write(file, "$value,\n")
      else
        write(file, "$value,")
      end
    end
    write(file, "    ],\n")

    #== write fragment symbols ==#
    write(file, "    \"symbols\" : [\n")
    for iatom in 1:length(symbols)
      write(file, "      ")
      symbol = symbols[iatom]
      if iatom == length(symbols)
        write(file, "\"$symbol\"\n")
      elseif iatom%5 == 0
        write(file, "\"$symbol\",\n")
      else
        write(file, "\"$symbol\",")
      end
    end
    write(file, "    ],\n")
    #== write fragment charge ==#
    write(file, "    \"molecular_charge\": $charge\n")
    write(file, "  },\n")

    #== write calculation driver and model information ==#
    write(file,"  \"driver\": \"energy\",\n")
    write(file,"  \"model\": {\n")
    write(file,"    \"method\": \"RHF\",\n")
    write(file,"    \"basis\": \"6-31G*\"\n")
    write(file,"  },\n")

    #== write keywords ==#
    write(file,"  \"keywords\": {\n")
    write(file,"    \"scf\": {\n")
    write(file,"      \"dele\":1.0,\n")
    write(file,"      \"rmsd\":2E-5\n")
    write(file,"    }\n")
    write(file,"  }\n")
    write(file,"}")
  end
end


text:
using SpikingNN
using Plots

dt = 0.1
spikes = [3, 20, 23, 30]
N = 100
synapse = Synapse.EPSP(ϵ₀ = 1, τm = 2, τs = 0.5)

# excite the synapse
excite!(synapse, spikes)

# get response
y = [synapse(t; dt = dt) for t in 1:N]

# plot results
scatter(dt .* spikes .- dt, zeros(length(spikes)), label = "Spikes")
plot!(dt .* collect(0:(N - 1)), y, label = "Output", minorticks = true)


text:
function colortype(bg::AbstractArray{T}) where {S, T <: Color{S}}
    RGB{S}
end
# function addcolor(A, bg::AbstractArray{T}) where {S, T <: Color{S}}
#     TC = RGB{S}
#     colors = Dict{eltype(A), TC}(s => rand(TC) for s in unique(A))
#     delete!(colors, 0)
#     bg1 = view(bg, axes(A)...)
#     ((x, y) -> if iszero(x) TC(y) else colors[x] end).(A, bg1);
# end

function addcolor(::Type{T}, A) where T
    colors = Dict{eltype(A), T}(s => rand(T) for s in unique(A))
    delete!(colors, 0)
    (x, y) -> if iszero(x) T(y) else colors[x] end
end
function addcolor(::Type{T}, A::AbstractArray{Bool}) where T
    (x, y) -> if iszero(x) T(y) else zero(T) end
end
function addcolor(::Type{T}, A::AbstractArray{<:Color}) where T
    (x, y) -> x
end
# function addcolor(A::AbstractArray{Bool}, bg::AbstractArray{T}) where {S, T <: Color{S}}
#     bg1 = view(bg, axes(A)[1:3]...)
#     ((x, y) -> if iszero(x) T(y) else zero(T) end).(A, bg1);
# end
function catcolor(bg::AbstractArray{T}, AA...) where {S, T <: Color{S}}
    TC = RGB{S}
    size3 = min(size(bg, 3), minimum(size(A, 3) for A in AA))
    size4 = 1 + sum(size(A, 4) for A in AA)
    size3, size4
    ret = Array{TC, 4}(undef, size(bg)[1:2]..., size3, size4)
    ret[:,:,:,1] .= TC.(view(bg, :, :, 1:size3))
    i4 = 2
    for A in AA
        col = addcolor(TC, view(A, :, :, 1:size3, :))
        ret[:,:,:, i4:i4+size(A,4)-1] .= col.(view(A, :, :, 1:size3, :), view(bg, :, :, 1:size3, :))
        i4 += size(A, 4)
    end
    ret
end
# function addbackground(A, bg::AbstractArray{T}, col=zero(T)) where T
#     Aoi = broadcast((x, y) -> ifelse(x, col, y), A, @view bg[axes(A)[1:3]...])
#     Aoi
# end
#
# function applycolors(f, bg, col::AbstractArray{T}) where T
#     ((x, y) -> if iszero(x) T(y) else col[x] end).(f, @view bg[axes(f)...]);
# end


text:
# This is a working script to plot an MTG using GraphMakie
using NetworkLayout, GraphMakie, GeometryBasics, LightGraphs
node_names = String[traverse(mtg, x -> join([string("(", x.id, ")"), x.MTG.symbol, x.MTG.index], " "))...]
links = String[traverse(mtg, x -> x.MTG.link)[2:end]...]
f, ax, p = graphplot(MetaGraph(mtg); layout = Buchheim(), nlabels = node_names, elabels = links)

# We can also use a custom Layout (like we use with other plots):
function mtg_layout(g::MetaGraph)
    x = zeros(nv(g))
    y = zeros(nv(g))
    z = zeros(nv(g))

    for i = 1:nv(g)
        node_string = "node_" * string(i)
        x[i] = meta_mtg[node_string][:XX]
        y[i] = meta_mtg[node_string][:YY]
        z[i] = meta_mtg[node_string][:ZZ]
    end

    return Point.(zip(x, y, z))
end

f, ax, p = graphplot(MetaGraph(mtg); layout = mtg_layout, nlabels = node_names, elabels = links)


text: push!( LOAD_PATH, "src/")

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Run tests on the block SS contour integral method

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@testset "block SS" begin
    nep=nep_gallery("dep0",3)

    # Circle
    λ,V=contour_block_SS(nep,radius=1.0,N=1000,σ=0.1,k=3,K=3);
    @test norm(compute_Mlincomb(nep,λ[1],V[:,1])) < sqrt(eps())

    # Ellipse
    λ,V=contour_block_SS(nep,radius=[1.0,2.0],N=1000,σ=0.1,k=3,K=3);
    @test norm(compute_Mlincomb(nep,λ[1],V[:,1])) < sqrt(eps())

    # JSIAM Mode
    λ,V=contour_block_SS(nep,radius=1.0,N=1000,σ=0.1,k=3,K=4,Shat_mode=:JSIAM);
    @test norm(compute_Mlincomb(nep,λ[1],V[:,1])) < sqrt(eps())

end


text:
using Statistics
using Random

################################################################

obj(Y, X, b) = mean((Y .- X * b) .^ 2.0) # objective to minimize

# modified example from GradDescent.jl documentation
@testset "linreg" begin

    Random.seed!(1) # set seed
    n = 1000 # number of observations
    d = 10   # number of covariates
    X = randn(n, d) # simulated covariates
    b = randn(d)    # generated coefficients
    ϵ = randn(n) * 0.1 # noise
    Y = X * b + ϵ # observed outcome

    epochs = 100 # number of epochs

    θ = randn(d) # initialize model parameters

    loss_first, g = grad(obj, Y, X, θ)
    loss_last = loss_first
    for i in 1:epochs
        loss_last, g = grad(obj, Y, X, θ)
        # println("Epoch: $i; loss = $loss_last")
        δ = 0.01 * g[4]
        θ -= δ
    end
    @test loss_last < loss_first / 10   # loss reduced at least 10 times

end

mutable struct Linear2{T}
    W::AbstractMatrix{T}
    b::AbstractVector{T}
end

linear_loss(m::Linear2, X, y) = mean((m.W * X .+ m.b .- y) .^ 2.0)

@testset "linreg 2" begin

    n_vars = 10
    n_out = 2
    n_obs = 1000
    X = randn(n_vars, n_obs)
    W_true = rand(n_out, n_vars)
    b_true = rand(n_out)
    ϵ = randn(n_out, n_obs) * 0.1 # noise
    Y = W_true * X .+ b_true .+ ϵ

    m = Linear2(rand(n_out, n_vars), rand(n_out))

    epochs = 100

    for i in 1:epochs
        val, g = grad(linear_loss, m, X, Y)
        # println("Epoch: $i; loss = $val")
        update!(m, g[2], (x, gx) -> x - 0.1 * gx)
    end

    @test isapprox(m.W, W_true; atol=0.1)
    @test isapprox(m.b, b_true; atol=0.1)

end


text:
using HTTP, JSON

function getBody(resp)
    jsonb=String(resp.body)
    bodyData = JSON.parse(jsonb)
end

function histoday()
    base_url="https://min-api.cryptocompare.com"
    path= "data/histoday?fsym=BTC&tsym=USD&limit=100"
    url="$base_url/$path"
    bd = getBody(HTTP.get(url))["Data"]
end

x = histoday()
# vcat(DataFrame.(x)...)


text:
export
    COC_ACTION,
    TERMINAL_POLAR_STATE,
    UAVState,
    convert,
    get_speed,
    State,
    to_polar,
    norm_angle,
    get_horizontal_distance,
    get_min_horizontal_distance

#=
This file defines the State type, which maintains the state of the
systems and is provided to the other components in order to allow
then to either step the system forward or make decisions based on
the state.
=#

# const values
const COC_ACTION = -1.

# this assumes that we do not include
# responding / not responding
# but this does not influence simulator
# because it is used as a flag, not for
# indexing into any tables
const TERMINAL_POLAR_STATE = ones(5) * 1e5

#=
State associated with a single uav.
=#
type UAVState
    x::Float64
    xdot::Float64
    y::Float64
    ydot::Float64
    psi::Float64 # heading
    psidot::Float64
    phi::Float64 # bank angle
    phidot::Float64
    function UAVState(x::Float64, xdot::Float64, y::Float64, ydot::Float64,
                        psi::Float64, psidot::Float64, phi::Float64, phidot::Float64)
        return new(x, xdot, y, ydot, psi, psidot, phi, phidot)
    end
end

# convert a UAVState type to a vector
function Base.convert(::Type{Vector{Float64}}, uav_state::UAVState)
    return [getfield(uav_state, f) for f in fieldnames(uav_state)]
end

# convert a vector to a UAVState type
function Base.convert(::Type{UAVState}, uav_state_vector::Vector{Float64})
    @assert length(uav_state_vector) == length(fieldnames(UAVState)) "vector must have the \
        same number of elements as there are fieldnames in the UAVState type"
    return UAVState(uav_state_vector...)
end

#=
Description:
Return the speed of the uav, calculated from state variables

Parameters:
uav_state_vector: state of uav as a vector

Return Value:
speed: float value for the speed of the uav
=#
function get_speed(uav_state_vector::Vector{Float64})
    return sqrt(uav_state_vector[2]^2 + uav_state_vector[4]^2)
end

#=
Description:
Convenience method for determining if two uav states are equal

Parameters:
- u1: first uav state
- u2: second uav state

Return Values:
- whether or not the two uav states are the same
=#
function Base.:(==)(u1::UAVState, u2::UAVState)
    return convert(Vector{Float64}, u1) == convert(Vector{Float64}, u2)
end

#=
The State type contains all the non-constant information needed by the uavs
to decide what action to take and all the non-constant information needed by
the dynamics type to step the entire environment forward in time.
=#
type State
    uav_states::Vector{UAVState}
    num_uavs::Int64
    # these are the coordination signals
    # emitted the previous timestep
    signals::Array{Symbol}
    responding_states::Vector{Int64}
    coordinator::Int64
    function State(uav_states::Vector{UAVState})
        num_uavs = length(uav_states)
        signals = [:no_signal, :no_signal]
        responding_states = zeros(Int64, num_uavs)
        coordinator = 1 # default is master
        return new(uav_states, num_uavs, signals, responding_states, coordinator)
    end
end

#=
Description:
Convert a state object from cartesian to polar coordinates.
Note that this function does not do bounds checking on the result.
This is performed after calling this function, in the uav
get_polar_state call.

Parameters:
- state: state object to be converted

Return Value:
- polar_state: the state in polar format
=#
function to_polar(state::State)
    @assert length(state.uav_states) == 2 "invalid number of uav_states in state: $state"

    sown = state.uav_states[1]
    sint = state.uav_states[2]

    dx = sint.x - sown.x
    dy = sint.y - sown.y
    xr = dx * cos(sown.psi) + dy * sin(sown.psi)
    yr = -dx * sin(sown.psi) + dy * cos(sown.psi)

    # heading is relative to ownship
    # and is caluated by subtracting the
    # heading of the intruder and normalizing
    # the resulting angle to between
    # 0 and 2pi
    pr = norm_angle(sint.psi - sown.psi)

    vown = sqrt(sown.xdot^2 + sown.ydot^2)
    vint = sqrt(sint.xdot^2 + sint.ydot^2)

    r = sqrt(xr^2 + yr^2)
    # norm theta to be between -pi and pi
    theta = norm_angle(atan2(yr, xr), neg_pi_to_pi = true)

    own_responding = state.responding_states[1]
    int_responding = state.responding_states[2]

    polar_state = [r, theta, pr, vown, vint, own_responding, int_responding]
    return polar_state

end

function to_polar_non_joint(state::State)
    @assert length(state.uav_states) == 2 "invalid number of uav_states in state: $state"

    sown = state.uav_states[1]
    sint = state.uav_states[2]

    dx = sint.x - sown.x
    dy = sint.y - sown.y
    xr = dx * cos(sown.psi) + dy * sin(sown.psi)
    yr = -dx * sin(sown.psi) + dy * cos(sown.psi)

    # heading is relative to ownship
    # and is caluated by subtracting the
    # heading of the intruder and normalizing
    # the resulting angle to between
    # 0 and 2pi
    pr = norm_angle(sint.psi - sown.psi)

    vown = sqrt(sown.xdot^2 + sown.ydot^2)
    vint = sqrt(sint.xdot^2 + sint.ydot^2)

    r = sqrt(xr^2 + yr^2)
    # norm theta to be between -pi and pi
    theta = norm_angle(atan2(yr, xr), neg_pi_to_pi = true)

    own_responding = state.responding_states[1]

    polar_state = [r, theta, pr, vown, vint, own_responding]
    return polar_state

end

#=
Description:
Normalize an angle between 0 and 360
If neg_pi_to_pi flag is true, then between -pi and pi

Parameters:
- angle: angle to Normalize
- neg_pi_to_pi: normalize to between -pi and pi

Return Value:
- noramlized angle
=#
function norm_angle(angle::Float64; neg_pi_to_pi::Bool = false)
    new_angle = ((angle % (2 * pi)) + 2 * pi) % (2 * pi)
    if neg_pi_to_pi
        if new_angle > pi
            new_angle -= 2 * pi
        end
    end
    return new_angle
end

#=
Description:
Returns the horizontal distance between the two uavs

Parameters:
- state: current state of the system

Return Values:
- horizontal_distance: the horizontal distance between the two uavs
=#
function get_horizontal_distance(state::State)
    sown = state.uav_states[1]
    sint = state.uav_states[2]
    horizontal_distance = sqrt((sown.x - sint.x)^2 + (sown.y - sint.y)^2)
    return horizontal_distance
end

#=
Description:
Get the minimum horizontal distance that occurred in a set of
intermediate states.

Parameters:
- state: the intermediate uav states (shape = (num_uavs, num states))

Return Value:
- min horizontal distance: min horizontal dist in meters, float
=#
function get_min_horizontal_distance(states::Array{Array{UAVState}})
    min_horizontal_dist = Inf
    for (sown, sint) in zip(states[1], states[2])
        horizontal_distance = sqrt((sown.x - sint.x)^2 + (sown.y - sint.y)^2)
        min_horizontal_dist = min(min_horizontal_dist, horizontal_distance)
    end
    return min_horizontal_dist
end

#=
Description:
Convenience method for determining if two states are equal

Parameters:
- s1: first state
- s2: second state

Return Values:
- whether or not the two states are the same
=#
function Base.:(==)(s1::State, s2::State)
    if s1.num_uavs != s2.num_uavs
        return false
    end
    for (u1, u2) in zip(s1.uav_states, s2.uav_states)
        if u1 != u2
            return false
        end
    end
    return true
end


text:
"""
    count_triplet_elements(array::Array{Int64,1}, sum::Int64)

Counting the number of three elements in the array, which sum is equal to the reference sum.

# Arguments
- `array::Array{Int64,1}`: Unsorted array
- `array::Array{Int64,1}`: Refernece sum

# Examples
```julia-repl
julia> import ClassicAlgorithmsCollections
julia> arr = [5, 1, 3, 4, 7]
julia> sum = 12
julia> ClassicAlgorithmsCollections.count_triplet_elements(arr, sum)
2

# Notes:
---

This should be implemented recursively instead of using three for-loops because it will
allow using any conditions.

"""
function count_triplet_elements(array::Array{Int64,1}, sum::Int64)

    # Initialize result
    ans = 0
    n = length(array)

    for i in 1:n-2
        for j in i+1:n-1
            for k in j+1:n
                if (array[i] + array[j] + array[k] == sum)
                    ans += 1
                end
            end
        end
    end
    return ans
end

max_stars_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_issues_repo_name: Anselmoo/ClassicAlgorithmsCollections

max_forks_repo_name: Anselmoo/ClassicAlgorithmsCollections


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using CompilerTools
using CompilerTools.LambdaHandling
using Base.Test

## Tests for CompilerTools.AliasAnalysis

function test_alias_2(x::Int, y::Int, z::Int, s::Int)
     A = rand(x, y)
     D_arr = rand(x, y)
     C = rand(x,y)
     B = D_arr

     for i = 1:x
          if( i%2 == 0 )
              for j = 1:y
              D_arr[i,j] /= 2
          end
          end
      end

      E = B
      if( s > 0 )
          E = A
      end

      return E

end

ast = CompilerTools.Helper.LambdaInfo(test_alias_2, (Int,Int,Int,Int), code_typed(test_alias_2, (Int,Int,Int,Int))[1])
linfo, body = lambdaToLambdaVarInfo(ast)
#cfg_2 = CompilerTools.CFGs.from_lambda(ast) :: CompilerTools.CFGs.CFG

#CompilerTools.AliasAnalysis.set_debug_level(3)
lives = CompilerTools.LivenessAnalysis.from_lambda(linfo, body)
handled = CompilerTools.AliasAnalysis.from_lambda(linfo, body, lives)
handled = map(x -> lookupVariableName(x, linfo), handled)
#println(handled)
#@test (in(:C, handled))
@test (in(:E, handled) == false)
@test (in(:B, handled) == false)
#@test (in(:D_arr, handled) == false)
#@test (in(:A, handled) == false)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiNetworkingV1NetworkPolicyIngressRule <: SwaggerModel
    from::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1NetworkPolicyPeer} } # spec name: from
    ports::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1NetworkPolicyPort} } # spec name: ports

    function IoK8sApiNetworkingV1NetworkPolicyIngressRule(;from=nothing, ports=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1NetworkPolicyIngressRule, Symbol("from"), from)
        setfield!(o, Symbol("from"), from)
        validate_property(IoK8sApiNetworkingV1NetworkPolicyIngressRule, Symbol("ports"), ports)
        setfield!(o, Symbol("ports"), ports)
        o
    end
end # type IoK8sApiNetworkingV1NetworkPolicyIngressRule

const _property_map_IoK8sApiNetworkingV1NetworkPolicyIngressRule = Dict{Symbol,Symbol}(Symbol("from")=>Symbol("from"), Symbol("ports")=>Symbol("ports"))
const _property_types_IoK8sApiNetworkingV1NetworkPolicyIngressRule = Dict{Symbol,String}(Symbol("from")=>"Vector{IoK8sApiNetworkingV1NetworkPolicyPeer}", Symbol("ports")=>"Vector{IoK8sApiNetworkingV1NetworkPolicyPort}")
Base.propertynames(::Type{ IoK8sApiNetworkingV1NetworkPolicyIngressRule }) = collect(keys(_property_map_IoK8sApiNetworkingV1NetworkPolicyIngressRule))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1NetworkPolicyIngressRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1NetworkPolicyIngressRule[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1NetworkPolicyIngressRule }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1NetworkPolicyIngressRule[property_name]

function check_required(o::IoK8sApiNetworkingV1NetworkPolicyIngressRule)
    true
end

function validate_property(::Type{ IoK8sApiNetworkingV1NetworkPolicyIngressRule }, name::Symbol, val)
end


text:
###############
# Commit Type #
###############

@ghdef mutable struct Commit
    sha::Union{String, Nothing}
    message::Union{String, Nothing}
    author::Union{Owner, Nothing}
    committer::Union{Owner, Nothing}
    commit::Union{Commit, Nothing}
    url::Union{HTTP.URI, Nothing}
    html_url::Union{HTTP.URI, Nothing}
    comments_url::Union{HTTP.URI, Nothing}
    parents::Union{Vector{Commit}, Nothing}
    stats::Union{Dict, Nothing}
    files::Union{Vector{Content}, Nothing}
    comment_count::Union{Int, Nothing}
end

Commit(sha::AbstractString) = Commit(Dict("sha" => sha))

namefield(commit::Commit) = commit.sha

###############
# API Methods #
###############

@api_default function commits(api::GitHubAPI, repo; options...)
    results, page_data = gh_get_paged_json(api, "/repos/$(name(repo))/commits"; options...)
    return map(Commit, results), page_data
end

@api_default function commit(api::GitHubAPI, repo, sha; options...)
    result = gh_get_json(api, "/repos/$(name(repo))/commits/$(name(sha))"; options...)
    return Commit(result)
end