text:
# HistoryRecorder
# maintained by @zsunberg

"""
A simulator that records the history for later examination

The simulation will be terminated when either
1) a terminal state is reached (as determined by `isterminal()` or
2) the discount factor is as small as `eps` or
3) max_steps have been executed

Keyword Arguments:
    - `rng`: The random number generator for the simulation
    - `capture_exception::Bool`: whether to capture an exception and store it in the history, or let it go uncaught, potentially killing the script
    - `show_progress::Bool`: show a progress bar for the simulation
    - `eps`
    - `max_steps`

Usage (optional arguments in brackets):

    hr = HistoryRecorder()
    history = simulate(hr, pomdp, policy, [updater [, init_belief [, init_state]]])
"""
mutable struct HistoryRecorder <: Simulator
    rng::AbstractRNG

    # options
    capture_exception::Bool
    show_progress::Bool

    # optional: if these are null, they will be ignored
    max_steps::Union{Nothing,Any}
    eps::Union{Nothing,Any}
end

# This is the only stable constructor
function HistoryRecorder(;rng=Random.GLOBAL_RNG,
                          eps=nothing,
                          max_steps=nothing,
                          capture_exception=false,
                          show_progress=false)
    return HistoryRecorder(rng, capture_exception, show_progress, max_steps, eps)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy) begin
    @req updater(::typeof(policy))
    up = updater(policy)
    @subreq simulate(sim, pomdp, policy, up)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater) begin
    @req initialstate(::typeof(pomdp))
    dist = initialstate(pomdp)
    @subreq simulate(sim, pomdp, policy, bu, dist)
end

function simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater=updater(policy))
    dist = initialstate(pomdp)
    return simulate(sim, pomdp, policy, bu, dist)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater, dist::Any) begin
    P = typeof(pomdp)
    S = statetype(pomdp)
    A = actiontype(pomdp)
    O = obstype(pomdp)
    @req initialize_belief(::typeof(bu), ::typeof(dist))
    @req isterminal(::P, ::S)
    @req discount(::P)
    @req gen(::P, ::S, ::A, ::typeof(sim.rng))
    b = initialize_belief(bu, dist)
    B = typeof(b)
    @req action(::typeof(policy), ::B)
    @req update(::typeof(bu), ::B, ::A, ::O)
end

function simulate(sim::HistoryRecorder,
                           pomdp::POMDP{S,A,O},
                           policy::Policy,
                           bu::Updater,
                           initialstate_dist::Any,
                           is::Any=rand(sim.rng, initialstate(pomdp))
                  ) where {S,A,O}

    initial_belief = initialize_belief(bu, initialstate_dist)
    max_steps = something(sim.max_steps, typemax(Int))
    if sim.eps != nothing
        max_steps = min(max_steps, ceil(Int,log(sim.eps)/log(discount(pomdp))))
    end

    if sim.show_progress
        if (sim.max_steps == nothing) && (sim.eps == nothing)
            error("If show_progress=true in a HistoryRecorder, you must also specify max_steps or eps.")
        end
        prog = Progress(max_steps, "Simulating..." )
    else
        prog = nothing
    end

    it = POMDPSimIterator(default_spec(pomdp),
                          pomdp,
                          policy,
                          bu,
                          sim.rng,
                          initial_belief,
                          is,
                          max_steps)

    history, exception, backtrace = collect_history(it, Val(sim.capture_exception), prog)

    if sim.show_progress
        finish!(prog)
    end

    return SimHistory(promote_history(history), discount(pomdp), exception, backtrace)
end

@POMDP_require simulate(sim::HistoryRecorder, mdp::MDP, policy::Policy) begin
    init_state = rand(sim.rng, initialstate(mdp))
    @subreq simulate(sim, mdp, policy, init_state)
end

@POMDP_require simulate(sim::HistoryRecorder, mdp::MDP, policy::Policy, initialstate::Any) begin
    P = typeof(mdp)
    S = statetype(mdp)
    A = actiontype(mdp)
    @req isterminal(::P, ::S)
    @req action(::typeof(policy), ::S)
    @req gen(::P, ::S, ::A, ::typeof(sim.rng))
    @req discount(::P)
end

function simulate(sim::HistoryRecorder,
                  mdp::MDP{S,A}, policy::Policy,
                  init_state::S=rand(sim.rng, initialstate(mdp))) where {S,A}

    max_steps = something(sim.max_steps, typemax(Int))
    if sim.eps != nothing
        max_steps = min(max_steps, ceil(Int,log(sim.eps)/log(discount(mdp))))
    end

    it = MDPSimIterator(default_spec(mdp),
                        mdp,
                        policy,
                        sim.rng,
                        init_state,
                        max_steps)

    if sim.show_progress
        if (sim.max_steps == nothing) && (sim.eps == nothing)
            error("If show_progress=true in a HistoryRecorder, you must also specify max_steps or eps.")
        end
        prog = Progress(max_steps, "Simulating..." )
    else
        prog = nothing
    end

    history, exception, backtrace = collect_history(it, Val(sim.capture_exception), prog)

    if sim.show_progress
        finish!(prog)
    end

    return SimHistory(promote_history(history), discount(mdp), exception, backtrace)
end

function collect_history(it, cap_ex::Val{true}, prog::Union{Progress,Nothing})
    exception = nothing
    backtrace = nothing
    history = NamedTuple[] # capturing part of the history is more important than this having a concrete type
    try
        for step in it
            push!(history, step)
            if prog !== nothing
                next!(prog)
            end
        end
    catch ex
        exception = ex
        backtrace = catch_backtrace()
    end
    return history, exception, backtrace
end

collect_history(it, cap_ex::Val{false}, prog::Nothing) = collect(it), nothing, nothing
function collect_history(it, cap_ex::Val{false}, prog::Progress)
    h = collect(begin
                    next!(prog)
                    step
                end for step in it)
    return h, nothing, nothing
end

"""
Promotes all NamedTuples in the history to the same type.
"""
function promote_history(hist::AbstractVector)
    if isconcretetype(eltype(hist))
        return hist
    elseif isempty(hist) # note, from above, also does not have concrete type
        return NamedTuple{(), Tuple{}}[]
    else
        # it would really astound me if this branch was type stable
        names = fieldnames(first(hist))
        types = fieldtypes(first(hist))
        for step in hist
            @assert fieldnames(step) == names
            types = map(promote_type, types, fieldtypes(step))
        end
        newtype = NamedTuple{names, Tuple{types...}}
        return convert(Vector{newtype}, hist)
    end
end


text:
using StochasticDiffEq, DiffEqJump, DiffEqBase, Statistics
using Test, LinearAlgebra

function regular_rate(out,u,p,t)
    out[1] = (0.1/1000.0)*u[1]*u[2]
    out[2] = 0.01u[2]
end

const dc = zeros(3, 2)
dc[1,1] = -1
dc[2,1] = 1
dc[2,2] = -1
dc[3,2] = 1

function regular_c(du,u,p,t,counts,mark)
    mul!(du,dc,counts)
end

rj = RegularJump(regular_rate,regular_c,2)
jumps = JumpSet(rj)
iip_prob = DiscreteProblem([999.0,1,0],(0.0,250.0))
jump_iipprob = JumpProblem(iip_prob,Direct(),rj)
@time sol = solve(jump_iipprob,TauLeaping())
@time sol = solve(jump_iipprob,SimpleTauLeaping();dt=1.0)
@time sol = solve(jump_iipprob,TauLeaping();dt=1.0,adaptive=false)
@time sol = solve(jump_iipprob,CaoTauLeaping();dt=1.0)
@time sol = solve(jump_iipprob,CaoTauLeaping())

N = 40_000
sol1 = solve(EnsembleProblem(jump_iipprob),SimpleTauLeaping();dt=1.0,trajectories = N)
sol2 = solve(EnsembleProblem(jump_iipprob),TauLeaping();dt=1.0,adaptive=false,save_everystep=false,trajectories = N)

mean1 = mean([sol1[i][end,end] for i in 1:N])
mean2 = mean([sol2[i][end,end] for i in 1:N])
@test mean1 ≈ mean2 rtol=1e-2

f(du,u,p,t) = (du .= 0)
g(du,u,p,t) = (du .= 0)
iip_sdeprob = SDEProblem(f,g,[999.0,1,0],(0.0,250.0))
jumpdiff_iipprob = JumpProblem(iip_sdeprob,Direct(),rj)
@time sol = solve(jumpdiff_iipprob,EM();dt=1.0)
@time sol = solve(jumpdiff_iipprob,ImplicitEM();dt=1.0,adaptive=false)

sol = solve(EnsembleProblem(jumpdiff_iipprob),EM();dt=1.0,trajectories = N)
meanX = mean([sol[i][end,end] for i in 1:N])
@test mean1 ≈ meanX rtol=1e-2

sol = solve(EnsembleProblem(jumpdiff_iipprob),ImplicitEM();dt=1.0,trajectories = N)
meanX = mean([sol[i][end,end] for i in 1:N])
@test mean1 ≈ meanX rtol=1e-2

iip_prob = DiscreteProblem([999,1,0],(0.0,250.0))
jump_iipprob = JumpProblem(iip_prob,Direct(),rj)
sol = solve(jump_iipprob,TauLeaping())

function rate_oop(u,p,t)
    [(0.1/1000.0)*u[1]*u[2],0.01u[2]]
end

function regular_c(u,p,t,counts,mark)
    dc*counts
end

rj = RegularJump(rate_oop,regular_c,2)
jumps = JumpSet(rj)
prob = DiscreteProblem([999.0,1,0],(0.0,250.0))
jump_prob = JumpProblem(prob,Direct(),rj)
sol = solve(jump_prob,TauLeaping(),reltol=5e-2)

sol2 = solve(EnsembleProblem(jump_prob),TauLeaping();dt=1.0,adaptive=false,save_everystep=false,trajectories = N)
mean2 = mean([sol2[i][end,end] for i in 1:N])
@test mean1 ≈ mean2 rtol=1e-2

foop(u,p,t) = [0.0,0.0,0.0]
goop(u,p,t) = [0.0,0.0,0.0]
oop_sdeprob = SDEProblem(foop,goop,[999.0,1,0],(0.0,250.0))
jumpdiff_prob = JumpProblem(oop_sdeprob,Direct(),rj)
@time sol = solve(jumpdiff_prob,EM();dt=1.0)
@time sol = solve(jumpdiff_prob,ImplicitEM();dt=1.0)

sol = solve(EnsembleProblem(jumpdiff_prob),EM();dt=1.0,trajectories = 10_000)
meanX = mean([sol[i][end,end] for i in 1:10_000])
@test mean1 ≈ meanX rtol=1e-2

sol = solve(EnsembleProblem(jumpdiff_prob),ImplicitEM();dt=1.0,trajectories = 1_000)
meanX = mean([sol[i][end,end] for i in 1:1_000])
@test mean1 ≈ meanX rtol=1e-1

ext: jl


text:
# ------------------------------------------------------------------------------------------
# ## How can we calculate $\pi$?
#
# Given a square of length $2r$, the square's area is
#
# $$A_{square} = (2r)^2 = 4r^2$$
#
# whereas the area of a circle with radius $r$ is
# $$A_{circle} = \pi r^2$$
#
#
#
# Therefore the ratio of the area of the circle to that of the square above is
#
# $$\frac{A_{circle}}{A_{square}} = \frac{\pi r^2}{4r^2} = \frac{\pi}{4}$$
#
# and we can define $\pi$ as
#
# $$\pi = 4\frac{A_{circle}}{A_{square}}$$
#
# This suggests a way to calculate $\pi$: if we have a square and the largest circle that
# fits inside that square, we can determine the ratio of areas of a circle and a square. We
# can calculate this ratio using a monte carlo simulation. We select random points inside a
# square, and we keep track of how often those points also fall inside the circle that fits
# perfectly inside that square.
#
# Given a large enough sampling points, $\frac{A_{circle}}{A_{square}}$ will be equal to the
# fraction of randomly chosen points inside the square that also fall inside the circle.
# Then we can figure out $\pi$!
#
# #### Pseudo-code
#
# Given the above, our algorithm for determining $\pi$ looks like this:
#
# 1. For each of $N$ iterations,
#     1. Select a random point inside a square of area $4r^2$ as Cartesian, $(x, y)$,
# coordinates.
#     1. Determine if the point also falls inside the circle embedded within this square of
# area $\pi r^2$.
#     1. Keep track of whether or not this point fell inside the circle. At the end of $N$
# iterations, you want to know $M$ -- the number of the $N$ random points that fell inside
# the circle!
# 1. Calculate $\pi$ as $4\frac{M}{N}$
#
# #### Exercise
#
# Write a function that calculates $\pi$ using Julia.
#
# The algorithm above should work for any value of $r$ that you choose to use. Make sure you
# make $N$ big enough that the value of $\pi$ is correct to at least a couple numbers after
# the decimal point!
#
# *Hint*:
#
# This will probably be easier if you center your circle and square at the coordinate (0, 0)
# and use a radius of 1. For example, to choose random coordinates within your square at
# position (x, y), you may want to choose x and y so that they are each a value between -1
# and +1. Then any point within a distance of 1 from (0, 0) will fall inside the circle!
#
#
#
#
# #### Solution
#
# In what follows are two functions, `calculate_pi` and `calculate_pi_2`. They use the same
# algorithm to determine `pi` with a couple minor differences in execution. Each of these
# functions is then called to show that as the number of samples, `N`, increases, our value
# for `pi` becomes more precise.
# ------------------------------------------------------------------------------------------

"""
calculate_pi(N = 1000)

Return the value of pi, calculated using a Monte Carlo simulation with N samples. N defaults to 1000.
The radius of the circle, r, defaults to 1.
"""
function calculate_pi(N = 1000)
    how_often_in_circle = 0
    for i in 1:N
        # Generate a random point (x, y) inside the square centered at (0, 0)
        # that has corners at (-1, 1), (-1, 1), (1, -1), and (1, 1)
        x, y = rand([-1, 1])*rand(), rand([-1, 1])*rand()
        # Check if the distance to (x, y) from (0, 0) is less than the radius of 1
        if sqrt(x^2 + y^2) < 1
            how_often_in_circle += 1
        end
    end
    return 4 * how_often_in_circle / N
end

calculate_pi.([10, 100, 1000, 10000, 100_000, 1_000_000_000])

"""
calculate_pi_2(N = 1000)

Return the value of pi, calculated using a Monte Carlo simulation with N samples. N defaults to 1000.
The radius of the circle, r, defaults to 1.
"""
function calculate_pi_2(N = 1000)
    how_often_in_circle = 0
    for i in 1:N
        # Generate a random point (x, y) inside the square centered at (0, 0)
        # that has corners at (-1, 1), (-1, 1), (1, -1), and (1, 1)
        x = 1-2*rand()
        y = 1-2*rand()
        # Check if the distance to (x, y) from (0, 0) is less than the radius of 1
        if x^2 + y^2 < 1
            how_often_in_circle += 1
        end
    end
    return 4 * how_often_in_circle / N
end

calculate_pi.([10, 100, 1000, 10000, 100_000, 1_000_000_000])

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8c1e6dae-b339-11eb-0e86-f1d18858689a
html"""

Section 3.8

 Advection and diffusion in 1D

Lecture Video

"""

# ╔═╡ e0c0dc94-277e-11eb-379e-83d064a93413
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra
end

# ╔═╡ 5648fa26-da0b-41d9-b13f-debd4e0485af
TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)

# ╔═╡ 00877a4a-277c-11eb-3ec0-e71e4094b404
md"""
# Evolution in time and space: Advection and diffusion in 1D
"""

# ╔═╡ 1b25b916-277c-11eb-0274-4b4fb946258d
md"""
So far we have been looking at dynamics in time, for example how does the temperature of the Earth change over time. But the Earth does not have a single, uniform temperature; rather, at a particular moment in time, different places on Earth are at different temperatures, and those different temperatures change over time due to several mechanisms.

In this notebook we will look at two fundamental mechanisms: **advection** and **diffusion**. Let's think about the temperature in the ocean. Since the ocean is a fluid that is in motion, a warm "parcel" of water can flow (move) to a new location due to the physical motion of the water itself; this is **advection**.  And even if the water doesn't move, temperature or a higher concentration of some substance dissolved in the fluid can spread out due to molecular mechanisms; this is **diffusion**.
"""

# ╔═╡ 956f5104-277d-11eb-291d-1faef485a5aa
md"""
In this notebook we will restrict ourselves to one spatial dimension (1D).
So we will think about the temperature $T$, for example, being a function

$$T(t, x)$$

of two independent variables:

- time, $t$
- space, $x$

We want to calculate a value of the temperature $T$ for each possible pair of values $(t, x)$, i.e. for all time ($>0$) and all positions.

The temperature at a given point will change due to different physical processes. We need to model this by writing down equations describing each physical process and how it affects the temperature. Since there are now *two* independent variables, $t$ and $x$, we can expect to end up with derivatives with respect to *both* of these variables, so that the rate of change of temperature in time at a given point depends also on gradients of the temperature in *space*. This will lead to a **partial differential equation** that relates *partial* derivatives of $T$.

In the context of climate modelling, we can think of $x$ as being the **latitude**, supposing that the temperature is the same at all points with the same latitude. In this way we could model the fact that the poles are cold and the equator is warm, and we could model how heat flows from hot to cold.

However, we clearly cannot model actual ocean currents like this, which would require two, or even three, spatial dimensions.

"""

# ╔═╡ b12e76db-1a18-465a-8955-dab29dfde611
md"""
## Visualising advection--diffusion
"""

# ╔═╡ c14470f2-d8a4-4d34-8470-09842b2576a3
md"""
Here is a visualisation of the physical processes of advection and diffusion in one dimension that we will discuss and build during this notebook.
"""

# ╔═╡ 30006c82-695d-40b1-8ded-22d03c3bff41
tt, results = evolve(advection_diffusion, xs, δt, (UU, DD))

# ╔═╡ b04a6f81-3ece-4521-b141-a2e416718948
md"""
U = $(@bind UU Slider(-1:0.01:1, show_value=true, default=0))

D = $(@bind DD Slider(-0.2:0.001:0.2, show_value=true, default=0))
"""

# ╔═╡ 6b2bfc73-d0a9-4a36-970d-c89149238284
md"""
time step = $(@bind n6 Slider(1:length(results), show_value=true))
"""

# ╔═╡ 21eb19f7-467b-4995-be65-8dede4eb7ac1
let
	p1 = plot(xs, results[n6], m=:o, xlim=(0, 1), ylim=(-3.1, 3.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, results[n6])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ 36328b70-277d-11eb-02c7-2f854c1466cc
md"""
# Temperature profiles and discretization
"""

# ╔═╡ 42190984-277d-11eb-1ac2-7d84516c3269
md"""
An ordinary differential equation needs an initial value for each variable. Similarly, we will need an initial *function* $T_0(x)$ that gives us the temperature at each position $x$. Let's suppose that the position is restricted to the interval $[0, L_x]$.

As usual, to represent the continuous function $T_0(x)$ on a computer, we will need to **discretise** it in some way, i.e. *approximate* the continuous function by a *finite* set of numbers in the computer.

The simplest (but far from the only!) discretisation method is to **sample** the function at discrete **grid points** (or **nodes**) $x_i$, for $i = 1, \dots, N_x$. For simplicity we will take these equally spaced, with spacing $x_{i+1} - x_i =: \delta x := L_x / N_x$.
"""

# ╔═╡ d2bed768-277e-11eb-32cf-41f1fedec3cb
md"""
For example, let's consider the following initial temperature profile:
"""

# ╔═╡ e6493da0-277e-11eb-22ff-29752652b576
# T₀(x) = sin(2π * x)^2 + 0.5

T₀(x) = sin(2π * x) + 2*cos(4π * x) + 0.2

# ╔═╡ 0d6271c0-286c-11eb-1c9c-3ba039b49d24
md"""
and define the grid points as follows:
"""

# ╔═╡ f17f7734-277e-11eb-25cf-5f2ba2db5aa3
begin
	Nₓ = 20
	Lₓ = 1.0

	δx = Lₓ / Nₓ

	xs = δx/2:δx:Lₓ
end

# ╔═╡ fa327c08-286b-11eb-0032-2384998a42a8
xs

# ╔═╡ 0db43be2-284c-11eb-2740-4379437fd70c
md"""
It turns out to be a good idea to take the grid points at the *centre* of each interval, so that we have $N_x$ intervals and $N_x$ grid points, starting at $x_1 = \delta x/2$ and finishing at $x_N = L_x - \delta x / 2$.
"""

# ╔═╡ 468a0590-2780-11eb-045c-d1f468fc4e50
md"""
We call such a function of $x$ at a given time a **temperature profile**. Let's draw it both as a function and as a heatmap:
"""

# ╔═╡ 6de1859c-277f-11eb-1ead-8b4794832d59
begin
	p1 = plot(0:0.001:Lₓ, T₀, label="T₀", lw=3)
	scatter!(xs, T₀.(xs), label="sampled")
	scatter!(xs, zero.(xs), label="x nodes", alpha=0.5, ms=3)

	xlabel!("x")
	ylabel!("T₀")

	for x in xs
		plot!([ (x, 0), (x, T₀(x)) ], ls=:dash, c=:black, label="", alpha=0.5)
	end

	hline!([0], ls=:dash, lab=false)

	p2 = temperature_heatmap(xs, T₀.(xs))

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]))

end

# ╔═╡ af30a0d0-2781-11eb-0274-ab423205facb
md"""
We will denote by $T^0_i$ the initial temperature at grid point number $i$.
"""

# ╔═╡ 646bc32e-284c-11eb-2ce8-5f64b1a49534
md"""
A useful way to think about $T^n_i$ is as some kind of (spatial) average of $T(t_n, x)$ over the interval of positions between neighbouring grid points, so $T_i$ is the average over the interval between $x_i - \frac{\delta x}{2}$ and $x_i + \frac{\delta x}{2}$. We can thus think of the following **piecewise constant** approximation to the original continuous function:
"""

# ╔═╡ 79ce4b10-284c-11eb-2258-2155f850171d
let
	δx = xs[2] - xs[1]

	p = plot(0:0.001:Lₓ, T₀, label="T₀", lw=1, ls=:dash)
	scatter!(xs, T₀.(xs), label="sampled")
	scatter!(xs, zero.(xs), label="x nodes", alpha=0.5, ms=3, lw=2)

	for i in 1:length(xs)
		plot!([ (xs[i] - δx/2, T₀(xs[i])), (xs[i] + δx/2, T₀(xs[i])) ], c=:green, lw=4, lab=false)

		plot!([ (xs[i] - δx/2, 0), (xs[i] - δx/2, T₀(xs[i])), (xs[i] + δx/2, T₀(xs[i])), (xs[i] + δx/2, 0)], c=:green, lw=1, lab=false, ls=:dash, alpha=0.3
		)
	end

	xlabel!("x")
	ylabel!("T₀(x)")
end

# ╔═╡ 2494daaa-2780-11eb-3084-2317924048ea
md"""
# Advection
"""

# ╔═╡ 29444ffe-2780-11eb-0875-095302b5d486
md"""
Now let's think of this profile as representing the temperature in each small volume, or "parcel", of fluid. Let's suppose that the fluid is moving to the right with a constant, uniform speed $U$. (**Uniform** here means that the speed is the same in all parts of the fluid.) Then the temperature profile should also *move with the fluid*! We call a quantity, such as the temperature, that is carried along with the fluid a **tracer**.

If we fix our attention at a single, fixed point in space, say the grid point $x_i$, the temperature there will vary over time, due to the fact that the fluid is moving past it. How it varies in time depends on the values at neighbouring grid points, since they determine how much heat will be transported *into* and *out of* the current cell.

[The point of view where we fix our attention at one point in space is called **Eulerian**. The alternative is to follow a parcel of fluid as it moves along in space; this is called **Lagrangian**.]
"""

# ╔═╡ 1dcb9690-6436-49f0-880f-23490fe28ea4
md"""
## Visualising fluxes in a fluid
"""

# ╔═╡ b63bb2e8-1d23-48fb-94b5-60d947465830
md"""
Let's visualise what happens as the fluid moves past a grid point, or rather the cell centered at a grid point. We will visualise tracer particles moving inside the fluid:
"""

# ╔═╡ e94a90c5-f2c1-4b5b-9946-7869ef7775a6
N = 5000

# ╔═╡ dd87fc01-4bf0-44f6-a9f6-560e433754a0
begin
	xx = ( abs.(-2 .+ 4 .* rand(N)) .^ 2) .- 1.5
	yy = rand(N)
end

# ╔═╡ 7ae9f5b8-10ea-42a7-aa01-0e04a7287c77
δ = 0.8

# ╔═╡ addab3e6-f189-41d6-badb-92f0323b6192
# assign colours to particles:

cs = map(xx) do x
	if -U * δ < x < 0
		1
	elseif 1 - (U * δ) < x < 1
		2
	else
		0
	end
end

# ╔═╡ 2f24e0c7-b05c-4f89-835a-081f8e6107e5
md"""
show particles entering and leaving in $\delta t$: $(@bind show_particles CheckBox())
"""

# ╔═╡ 75bc87be-2b66-46b5-8de8-428a63655815
md"""
t = $(@bind t Slider(0:0.001:2, show_value=true, default=0))
"""

# ╔═╡ f684dd94-f1c7-4f79-9776-3a06b8eec39b
begin
	plot([0, 1, 1, 0, 0], [0, 0, 1, 1, 0], series=:shape, alpha=0.5, fill=true, ratio=1, label=false, leg=false)

	new_xx = xx .+ U .* t

	scatter!(xx .+ U .* t, yy, ms=1.5, alpha=0.1, c=:gray)

	if show_particles
		scatter!(new_xx[cs .!= 0], yy[cs .!= 0], ms=1.5, alpha=0.5, c=cs[cs .!= 0])
	end

	plot!([-1.5, 2], [0, 0], c=:black)
	plot!([-1.5, 2], [1, 1], c=:black)

	xlims!(-2, 2)
	ylims!(-0.1, 1.1)

	as_svg(plot!(axis=true, yticks=[0, 1]))
end

# ╔═╡ 3437e53b-9dd0-4afe-a1bd-a556871d1799
md"""
## Time stepping
"""

# ╔═╡ 65df7158-60dc-4809-82a3-913a79bcfc75
md"""
We want to model how the temperature profile changes in time due to the flow of the fluid. We'll do so by looking at each cell and asking how much heat enters and leaves the cell in a given time step, of duration $\delta t$.
"""

# ╔═╡ 7256778a-2785-11eb-0369-f3b43d5dd203
md"""
Let's call $T^n_i$ the approximate (unknown) average value of $T$ in the cell at position $x_i$ and at the $n$th time step $t_n$, i.e. an approximation of $T(t_n, x_i)$, where $t_n = n \, \delta t$.

Then $T^{n+1}_i \simeq T(t_n + \delta t, x_i)$ and $T^{n}_{i+1} \simeq T(t_n, x_i + \delta x)$.

[Note that the superscript $n$ in these algorithms does not mean a power; it's just a label for the time step. We could write $T_i^{(n)}$ instead, but that is annoying to both write and read, so we omit the parentheses.]

"""

# ╔═╡ 44433a34-2782-11eb-0079-837c9306c5bd
md"""
Suppose the fluid is moving to the right with speed $U$. During a time step of duration $\delta t$, the temperature $T^n_i$ at cell $i$ changes for two reasons:

- some heat enters cell $i$
- some heat leaves cell $i$

Note that most of the fluid that starts within cell $i$ remains within that cell during the time step (if the time step is short enough), as we see from the visualisation above.

To calculate how much heat enters and leaves, note that only heat in the region of fluid within a distance $U \, \delta t$ from the boundary of the cell will cross into that cell. So a *proportion* $(U \, \delta t) / \delta x$ of the amount in cell $i$ crosses the boundary.

[We will blur the distinction between "amount of heat" and temperature.]

Hence, roughly an amount $T^n_i (U \delta t) / \delta x$ will leave cell number $i$ and cross into cell $i+1$ (the cell to the right). Similarly, an amount $T^n_{i-1} (U \delta t) / \delta x$ will *enter* cell $i$ from the neighbouring cell $i-1$ on the left.

Hence we arrive at the following:

$$T^{n+1}_i = T^{n}_i + (T^n_{i-1} - T^n_{i})\, U \, \delta t / \delta x.$$

Note that on the right-hand side we have quantities at the time step $n$, and on the left at time step $n+1$. So this tells us how to *update* our quantities from slice $n$ to slice $n+1$.

"""

# ╔═╡ 87e2be25-227c-498c-94fa-6e404c8918f1
md"""
## Continuous limit: Advection equation PDE
"""

# ╔═╡ 72c0ab0c-2781-11eb-1f59-9b22a52b0be0
md"""
Rearranging the previous equation we get

$$\frac{T^{n+1}_i - T^{n}_i}{\delta t} =  \frac{T^n_{i-1} - T^n_{i}}{\delta x}\,  U.$$
"""

# ╔═╡ e5761990-278b-11eb-134e-7954b577b1ac
md"""
Taking the continuum limit when $\delta t \to 0$ and $\delta x \to 0$, we recognise the definition of **partial derivatives** with respect to time and space variables from multivariable calculus. (Note the different indices that change on the two sides of the equation.)

Denoting these partial derivatives using $\partial$, we arrive at the **advection equation**:

$$\frac{\partial T(t, x)}{\partial t} = -U \frac{\partial T(t, x)}{\partial x},$$

or for short

$$\frac{\partial T}{\partial t} = -U \frac{\partial T}{\partial x}.$$

Since $T$ is a function of both $x$ and $t$, and this equation involves partial derivatives with respect to both of the independent variables, this is a **partial differential equation** (PDE). It describes how the function $T(t, x)$ changes continuously as a function both of time and space.

Although there are some analytical methods to solve PDEs, often it's necessary to use numerical methods. Here we'll look at simple numerical methods to solve such equations.
"""

# ╔═╡ 2033364e-278c-11eb-2936-17598ce14a41
md"""
## Numerics for the advection equation
"""

# ╔═╡ e9a37908-278c-11eb-278e-9bd155f0cae6
md"""
Let's return to the version of the equation in which the value at the *following* time step is isolated:

$$T^{n+1}_i = T^{n}_i - \left( U \frac{\delta t}{\delta x} \right) (T^n_{i} - T^n_{i-1}).$$

In the last term on the right-hand side, we see that we require combinations of values of $T$ at the *same* time step from *different* places, with certain coefficients.
"""

# ╔═╡ bcf1ceca-f557-4d75-9058-bbaa58665fb7
md"""
There are many approaches to implementing this numerically. The simplest is to directly transcribe the equation for the $i$th entry of the vector.

Calling `T` the current vector, i.e. $\mathbf{T}^n := (T^n_i)_{i=1, \ldots, N_x}$, and `T′` the new vector at the next time step, we have the following basic expression:

	T′[i] = T[i] + δt * U * (T[i-1] - T[i]) / δx

But now we realise a problem: What should we do when $i=1$? This will try to access the index 0 of the vector `T`, which does not exist!
"""

# ╔═╡ 3736a25e-4dec-46ac-9bf6-9712e3d00e7a
md"""

### Boundary conditions

This illustrates the necessity of choosing **boundary conditions** that specify what happens at the edge of the domain.

For simplicity we will choose to use **periodic boundary conditions**. This is a convenient mathematical fiction that allows us to treat all cells as being on the same footing, by wrapping the system around a torus, so that cells $i=1$ and $i=N_x$ are neighbours.
"""

# ╔═╡ e542a8da-284e-11eb-3297-6bbbf052284b
md"""
We can then write this as follows, where we separate out the case $i=1$:
"""

# ╔═╡ b15f4f44-284b-11eb-37c5-ab0153f7fe92
function advection(T, δt, δx, U)
	N = length(T)
	T′ = similar(T)  # create new vector of the same length

	# bulk cells:
	for i in 2:N
		T′[i] = T[i] - δt * U * (T[i] - T[i-1]) / δx
	end

	# boundary cells:
	T′[1] = T[1] - δt * U * (T[1] - T[N]) / δx   # periodic

	return T′
end

# ╔═╡ fcbec610-d9fc-4e41-8e76-729dbbc61d92
md"""
This performs a single time step of the advection equation; it takes in the current vector of $T$s and returns the new $T$s after the step.

Note that this is just like a step of the Euler method for solving ODEs, but where many spatial coordinates are updated at the same time. Effectively we are solving a system of coupled ODEs!
"""

# ╔═╡ af79e360-286e-11eb-2a4d-3d6d7564088c
δt = 0.001;

# ╔═╡ dce9e53a-28f4-11eb-070b-17e10779a38b
U = 0.2;

# ╔═╡ 02a893e4-2852-11eb-358a-371459191da7
ts, evolution = evolve(advection, xs, δt, U)

# ╔═╡ e6ae447e-2851-11eb-3fe1-096459167f2b
@bind n Slider(1:length(evolution), show_value=true)

# ╔═╡ 014e2530-2852-11eb-103f-1d647cb999b0
let
	p1 = plot(xs, evolution[n], m=:o, xlim=(0, 1), ylim=(-1.1, 1.1), title="t = $(round(ts[n], digits=2))", leg=false)

	p2 = temperature_heatmap(xs, evolution[n])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]))
end

# ╔═╡ 8c05e3cc-2858-11eb-1e1c-9781c30738c3
md"""
Unfortunately this does *not* behave as we expect: instead of preserving the shape of the profile over time, it is decaying. This is due to the way we are approximating.

A better way to discretize the spatial derivative is using the following **centered difference**:

$$\frac{\partial T(t_n, x_i)}{\partial x} \simeq \frac{T^n_{i+1} - T^n_{i-1}}{2 \delta x}$$

"""

# ╔═╡ a29fecac-285a-11eb-14b0-9313f8994fbb
function advection2(T, δt, δx, U)
	N = length(T)
	T′ = similar(T)  # create new vector of the same length

	for i in 2:N-1
		T′[i] = T[i] - δt * U * (T[i+1] - T[i-1]) / (2δx)
	end

	# periodic boundary:
	T′[1] = T[1] - δt * U * (T[2] - T[N]) / (2δx)
	T′[N] = T[N] - δt * U * (T[1] - T[N-1]) / (2δx)

	return T′
end

# ╔═╡ e42ec13e-285a-11eb-3cc0-7dc41ed5495b
ts2, evolution2 = evolve(advection2, xs, δt, 0.1)

# ╔═╡ f60a8b5e-285a-11eb-0d35-8daf23cf92ae
n2_slider = @bind n2 Slider(1:length(evolution2), show_value=true)

# ╔═╡ f1b5d130-285a-11eb-001c-67035925f43d
let
	p1 = plot(xs, evolution2[n2], m=:o, xlim=(0, 1), ylim=(-3.1, 3.1), title="t = $(round(ts2[n2], digits=2))", leg=false)

	p2 = temperature_heatmap(xs, evolution2[n2])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]))
end

# ╔═╡ c59388ea-286e-11eb-0f21-eb18e5ba516f
md"""
# Diffusion
"""

# ╔═╡ 3c944998-2888-11eb-087d-492b9d0ee32e
md"""
Another key physical process is **diffusion**. This models how temperature or mass spreads out from hot or high concentration regions towards regions where it is cold or where there is a low concentration.

## Physical mechanism: Random walks

The physical mechanism behind this is **random motion**: this is the continuous limit of equations describing the evolution of the probability distribution in space and time of a cloud of random walkers.

This is the same process that we studied in lecture 2.6.
Using our current notation, there we showed that the probability distribution of a cloud of random walkers satisfies the following time evolution:

$$p^{n+1}_i = \frac{1}{2}(p^n_{i-1} + p^n_{i+1})$$

If now we say that the walkers jump only with a certain probability, with a large probability to stay in the same place, and that these random walkers are the carriers of heat, then we get

$$T^{n+1}_i = \kappa (T^n_{i-1} - 2 T^n_i + T^n_{i+1}).$$

Watch [this video](https://www.youtube.com/watch?v=a3V0BJLIo_c) from last semester's class to see Grant Sanderson explaining this.
"""

# ╔═╡ 6ac74e34-ed58-4903-8c53-82be13b6c21f
md"""
## Continuous limit: Heat equation PDE
"""

# ╔═╡ de42149c-85ce-4e73-8503-84f64a173cbb
md"""
Introducing $\delta x$ as the spatial discretisation, and $\delta t$ as the time step, we get

$$T^{n+1}_i = \kappa \frac{\delta t}{\delta x^2}  (T^n_{i-1} - 2 T^n_i + T^n_{i+1}).$$

"""

# ╔═╡ ef42d541-74a1-433a-9773-5e6cca525350
md"""
The continuous limit is the following **heat equation** or **diffusion equation**:
"""

# ╔═╡ 6b7cea44-2888-11eb-0208-990860d6a152
md"""
$$\frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2}.$$
"""

# ╔═╡ 83a1e1f5-0946-422c-83f4-d7a19e9c0789
md"""
Here, $\kappa$ is the **heat diffusivity**, which says how quickly heat spreads out. In the context of diffusion of mass the equivalent is the **diffusion coefficient**, $D$.
"""

# ╔═╡ 68db3372-2888-11eb-1b03-b5ebca4c2bd5
md"""
To obtain a numerical method to solve this equation, we again need to discretise this, in particular the second derivative. One possible discretisation is

$$\frac{\partial^2 T}{\partial x^2}(t_n, x_i) \simeq \frac{T^n_{i+1} - 2 T^n_i + T^n_{i-1}}{\delta x^2}.$$
"""

# ╔═╡ d6131ad0-2889-11eb-3085-15d17e33ee7a
md"""
This may again be transcribed directly into code:
"""

# ╔═╡ 630314bc-2868-11eb-1b93-b7b08a4b2887
function diffusion(T, δt, δx, D)
	N = length(T)
	T′ = similar(T)  # create new vector of the same length

	for i in 2:N-1
		T′[i] = T[i] + δt * D * (T[i+1] -2T[i] + T[i-1]) / (δx^2)
	end

	# periodic boundary:
	T′[1] = T[1] + δt * D * (T[2] - 2T[1] + T[N]) / (δx^2)
	T′[N] = T[N] + δt * D * (T[1] - 2T[N] + T[N-1]) / (δx^2)

	return T′
end

# ╔═╡ 121255d2-288a-11eb-1fa5-9db68af8c232
ts3, evolution3 = evolve(diffusion, xs, δt, 0.01)

# ╔═╡ 09bc3c40-288a-11eb-0339-59f0b70e03a3
@bind n3 Slider(1:length(evolution2), show_value=true)

# ╔═╡ 175d9902-288a-11eb-3700-390ccd1caa5b
let
	p1 = plot(xs, evolution3[n3], m=:o, xlim=(0, 1), ylim=(-3.1, 3.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, evolution3[n3])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ e74d920a-28fa-11eb-3c91-9133a01effc5
@bind n4 Slider(1:length(evolution4), show_value=true)

# ╔═╡ dc7b6328-28fa-11eb-38b7-71a6f8d0a751
let
	p1 = plot(xs, evolution4[n4], m=:o, xlim=(0, 1), ylim=(-0.1, 1.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, evolution4[n4])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ e63cfa84-2889-11eb-1ea2-51726645ddd9
md"""
# The advection--diffusion PDE
"""

# ╔═╡ eee3008e-2889-11eb-088a-73aff304e736
md"""
Finally we can combine both mechanisms, to describe a tracer that is both being advected at a constant speed and diffusing. This basically utilises the composition of the advection and diffusion functions:
"""

# ╔═╡ ffd2a838-2889-11eb-1a7c-b35992543b8a
function advection_diffusion(T, δt, δx, (U, D))
	temp = advection2(T, δt, δx, U)
	return diffusion(temp, δt, δx, D)
end

# ╔═╡ f6fa3770-288d-11eb-32de-f95e03705791
ts5, evolution5 = evolve(advection_diffusion, xs, δt, (1.0, 0.01))

# ╔═╡ 6eb00a02-288d-11eb-354b-b56cf5a8380e
@bind n5 Slider(1:length(evolution5), show_value=true)

# ╔═╡ 65126bfc-288d-11eb-2bfc-493588365164
let
	p1 = plot(xs, evolution5[n5], m=:o, xlim=(0, 1), ylim=(-1.1, 1.1), title="t = $(round(ts3[n3], digits=2))", leg=false)
	p2 = temperature_heatmap(xs, evolution5[n5])

	plot(p1, p2, layout = grid(2, 1, heights=[0.9, 0.1]), clim=(-1, 1))
end

# ╔═╡ 575a5f3c-2780-11eb-2119-27a4114ceac5
md"""
# Function library
"""

# ╔═╡ 5a3eec86-2780-11eb-0341-39a5c343fc52
function temperature_heatmap(x, T)

	p = heatmap(x, [0.], collect(T'),
			   clims=(-1., 1.), cbar=false, xticks=nothing, yticks=nothing)

	return p
end

# ╔═╡ 9187350a-2851-11eb-05f0-d3a6eef190fe
function evolve(method, xs, δt, U, t_final=10.0, f₀=T₀)

	T = f₀.(xs)
	δx = xs[2] - xs[1]

	t = 0.0
	ts = [t]

	results = [T]

	while t < t_final
		T′ = method(T, δt, δx, U)  # new
		push!(results, T′)

		t += δt
		push!(ts, t)

		T = copy(T′)

	end

	return ts, results
end

# ╔═╡ Cell order:
# ╟─8c1e6dae-b339-11eb-0e86-f1d18858689a
# ╠═e0c0dc94-277e-11eb-379e-83d064a93413
# ╠═5648fa26-da0b-41d9-b13f-debd4e0485af
# ╟─00877a4a-277c-11eb-3ec0-e71e4094b404
# ╟─1b25b916-277c-11eb-0274-4b4fb946258d
# ╟─956f5104-277d-11eb-291d-1faef485a5aa
# ╟─b12e76db-1a18-465a-8955-dab29dfde611
# ╟─c14470f2-d8a4-4d34-8470-09842b2576a3
# ╟─30006c82-695d-40b1-8ded-22d03c3bff41
# ╟─b04a6f81-3ece-4521-b141-a2e416718948
# ╟─6b2bfc73-d0a9-4a36-970d-c89149238284
# ╟─21eb19f7-467b-4995-be65-8dede4eb7ac1
# ╟─36328b70-277d-11eb-02c7-2f854c1466cc
# ╟─42190984-277d-11eb-1ac2-7d84516c3269
# ╟─d2bed768-277e-11eb-32cf-41f1fedec3cb
# ╠═e6493da0-277e-11eb-22ff-29752652b576
# ╟─0d6271c0-286c-11eb-1c9c-3ba039b49d24
# ╠═f17f7734-277e-11eb-25cf-5f2ba2db5aa3
# ╠═fa327c08-286b-11eb-0032-2384998a42a8
# ╟─0db43be2-284c-11eb-2740-4379437fd70c
# ╟─468a0590-2780-11eb-045c-d1f468fc4e50
# ╠═6de1859c-277f-11eb-1ead-8b4794832d59
# ╟─af30a0d0-2781-11eb-0274-ab423205facb
# ╟─646bc32e-284c-11eb-2ce8-5f64b1a49534
# ╟─79ce4b10-284c-11eb-2258-2155f850171d
# ╟─2494daaa-2780-11eb-3084-2317924048ea
# ╟─29444ffe-2780-11eb-0875-095302b5d486
# ╟─1dcb9690-6436-49f0-880f-23490fe28ea4
# ╟─b63bb2e8-1d23-48fb-94b5-60d947465830
# ╠═e94a90c5-f2c1-4b5b-9946-7869ef7775a6
# ╠═dd87fc01-4bf0-44f6-a9f6-560e433754a0
# ╠═7ae9f5b8-10ea-42a7-aa01-0e04a7287c77
# ╟─addab3e6-f189-41d6-badb-92f0323b6192
# ╟─2f24e0c7-b05c-4f89-835a-081f8e6107e5
# ╟─75bc87be-2b66-46b5-8de8-428a63655815
# ╟─f684dd94-f1c7-4f79-9776-3a06b8eec39b
# ╟─3437e53b-9dd0-4afe-a1bd-a556871d1799
# ╟─65df7158-60dc-4809-82a3-913a79bcfc75
# ╟─7256778a-2785-11eb-0369-f3b43d5dd203
# ╟─44433a34-2782-11eb-0079-837c9306c5bd
# ╟─87e2be25-227c-498c-94fa-6e404c8918f1
# ╟─72c0ab0c-2781-11eb-1f59-9b22a52b0be0
# ╟─e5761990-278b-11eb-134e-7954b577b1ac
# ╟─2033364e-278c-11eb-2936-17598ce14a41
# ╟─e9a37908-278c-11eb-278e-9bd155f0cae6
# ╟─bcf1ceca-f557-4d75-9058-bbaa58665fb7
# ╟─3736a25e-4dec-46ac-9bf6-9712e3d00e7a
# ╟─e542a8da-284e-11eb-3297-6bbbf052284b
# ╠═b15f4f44-284b-11eb-37c5-ab0153f7fe92
# ╟─fcbec610-d9fc-4e41-8e76-729dbbc61d92
# ╠═af79e360-286e-11eb-2a4d-3d6d7564088c
# ╠═dce9e53a-28f4-11eb-070b-17e10779a38b
# ╠═02a893e4-2852-11eb-358a-371459191da7
# ╠═e6ae447e-2851-11eb-3fe1-096459167f2b
# ╟─014e2530-2852-11eb-103f-1d647cb999b0
# ╟─8c05e3cc-2858-11eb-1e1c-9781c30738c3
# ╠═a29fecac-285a-11eb-14b0-9313f8994fbb
# ╠═e42ec13e-285a-11eb-3cc0-7dc41ed5495b
# ╠═f60a8b5e-285a-11eb-0d35-8daf23cf92ae
# ╠═f1b5d130-285a-11eb-001c-67035925f43d
# ╟─c59388ea-286e-11eb-0f21-eb18e5ba516f
# ╟─3c944998-2888-11eb-087d-492b9d0ee32e
# ╟─6ac74e34-ed58-4903-8c53-82be13b6c21f
# ╟─de42149c-85ce-4e73-8503-84f64a173cbb
# ╟─ef42d541-74a1-433a-9773-5e6cca525350
# ╟─6b7cea44-2888-11eb-0208-990860d6a152
# ╟─83a1e1f5-0946-422c-83f4-d7a19e9c0789
# ╟─68db3372-2888-11eb-1b03-b5ebca4c2bd5
# ╟─d6131ad0-2889-11eb-3085-15d17e33ee7a
# ╠═630314bc-2868-11eb-1b93-b7b08a4b2887
# ╠═121255d2-288a-11eb-1fa5-9db68af8c232
# ╠═09bc3c40-288a-11eb-0339-59f0b70e03a3
# ╠═175d9902-288a-11eb-3700-390ccd1caa5b
# ╠═e74d920a-28fa-11eb-3c91-9133a01effc5
# ╟─dc7b6328-28fa-11eb-38b7-71a6f8d0a751
# ╟─e63cfa84-2889-11eb-1ea2-51726645ddd9
# ╟─eee3008e-2889-11eb-088a-73aff304e736
# ╠═ffd2a838-2889-11eb-1a7c-b35992543b8a
# ╠═f6fa3770-288d-11eb-32de-f95e03705791
# ╠═6eb00a02-288d-11eb-354b-b56cf5a8380e
# ╠═65126bfc-288d-11eb-2bfc-493588365164
# ╟─575a5f3c-2780-11eb-2119-27a4114ceac5
# ╟─5a3eec86-2780-11eb-0341-39a5c343fc52
# ╠═9187350a-2851-11eb-05f0-d3a6eef190fe


text:
should_precompile = true

# Don't edit the following! Instead change the script for `snoop_bot`.
ismultios = true
ismultiversion = true
# precompile_enclosure
@static if !should_precompile
    # nothing
elseif !ismultios && !ismultiversion
    include("../deps/SnoopCompile/precompile/precompile_AcuteML.jl")
    _precompile_()
else
    @static if Sys.islinux()
        @static if v"1.3.0-DEV" <= VERSION <= v"1.3.9"
            include("../deps/SnoopCompile/precompile/linux/1.3/precompile_AcuteML.jl")
            _precompile_()
        elseif v"1.4.0-DEV" <= VERSION <= v"1.4.9"
            include("../deps/SnoopCompile/precompile/linux/1.4/precompile_AcuteML.jl")
            _precompile_()
        else
        end

    elseif Sys.iswindows()
        @static if v"1.3.0-DEV" <= VERSION <= v"1.3.9"
            include("../deps/SnoopCompile/precompile/windows/1.3/precompile_AcuteML.jl")
            _precompile_()
        elseif v"1.4.0-DEV" <= VERSION <= v"1.4.9"
            include("../deps/SnoopCompile/precompile/windows/1.4/precompile_AcuteML.jl")
            _precompile_()
        else
        end

    elseif Sys.isapple()
        @static if v"1.3.0-DEV" <= VERSION <= v"1.3.9"
            include("../deps/SnoopCompile/precompile/apple/1.3/precompile_AcuteML.jl")
            _precompile_()
        elseif v"1.4.0-DEV" <= VERSION <= v"1.4.9"
            include("../deps/SnoopCompile/precompile/apple/1.4/precompile_AcuteML.jl")
            _precompile_()
        else
        end

    else
    end

end # precompile_enclosure


text:
@testset "classification" begin

    generate_caos_rules("data/S10593.nex", "data/output")
    tree, character_labels, taxa_labels = load_tree("data/output")

    @testset "classify_sequence" begin
        classification = classify_new_sequence(tree, character_labels, taxa_labels, "data/dna_seq.txt", "data/output")
        @test classification == "Tubificoides_amplivasatus_II_CE1247"
    end

    generate_caos_rules("data/E1E2L1.nex", "data/output2", protein=true)
    tree2, character_labels2, taxa_labels2 = load_tree("data/output2")

    @testset "classify_sequence" begin

        classification2 = classify_new_sequence(tree2, character_labels2, taxa_labels2, "data/protein_seq.txt", "data/output2", protein=true)

        @test classification2 == "HPV200"

        # @test classification == "Tubificoides_amplivasatus_II_CE1247"
    end
end


text:
using AutoHashEquals

@auto_hash_equals struct RDateCompound <: RDate
    parts::Vector{RDate}
end

apply(rdate::RDateCompound, date::Dates.Date, cal_mgr::CalendarManager) = Base.foldl((x,y) -> apply(y, x, cal_mgr), rdate.parts, init=date)
multiply_no_roll(rdate::RDateCompound, count::Integer) = RDateCompound(map(x -> multiply_no_roll(x, count), rdate.parts))
Base.:-(rdate::RDateCompound) = RDateCompound(map(-, rdate.parts))
Base.:+(x::RDateCompound, y::RDate) = RDateCompound(vcat(x.parts, y))
Base.:+(x::RDate, y::RDateCompound) = RDateCompound(vcat(x, y.parts))
combine(left::RDate, right::RDate) = RDateCompound([left,right])
function Base.show(io::IO, rdate::RDateCompound)
    for (i,part) in enumerate(rdate.parts)
        if i > 1 print(io, "+") end
        show(io, part)
    end
end

@auto_hash_equals struct RDateRepeat <: RDate
    count::Int64
    part::RDate
end

apply(rdate::RDateRepeat, date::Dates.Date, cal_mgr::CalendarManager) = Base.foldl((x,y) -> apply(y, x, cal_mgr), fill(rdate.part, rdate.count), init=date)
multiply_no_roll(rdate::RDateRepeat, count::Integer) = RDateRepeat(rdate.count, multiply_no_roll(rdate.part, count))
Base.:-(rdate::RDateRepeat) = RDateRepeat(rdate.count, -rdate.part)
Base.show(io::IO, rdate::RDateRepeat) = (print(io, "$(rdate.count)*roll("), show(io, rdate.part), print(io,")"))

Base.:+(left::RDate, right::RDate) = combine(left, right)
Base.:-(left::RDate, right::RDate) = combine(left, -right)
multiply_roll(rdate::RDate, count::Integer) = count >= 0 ? RDateRepeat(count, rdate) : RDateRepeat(-count, -rdate)

@auto_hash_equals struct CalendarAdj{R <: RDate, S <: HolidayRoundingConvention} <: RDate
    calendar_names::Vector{String}
    part::R
    rounding::S

    CalendarAdj(calendar_names, part::R, rounding::S) where {R <: RDate, S <: HolidayRoundingConvention} = new{R, S}(calendar_names, part, rounding)
end

function apply(rdate::CalendarAdj, date::Dates.Date, cal_mgr::CalendarManager)
    base_date = apply(rdate.part, date, cal_mgr)
    cal = calendar(cal_mgr, rdate.calendar_names)
    apply(rdate.rounding, base_date, cal)
end
multiply_no_roll(rdate::CalendarAdj, count::Integer) = CalendarAdj(rdate.calendar_names, multiply_no_roll(rdate.part, count), rdate.rounding)
Base.:-(x::CalendarAdj) = CalendarAdj(x.calendar_names, -x.part, x.rounding)
Base.show(io::IO, rdate::CalendarAdj) = print(io, "$(rdate.part)@$(join(rdate.calendar_names, "|"))[$(rdate.rounding)]")


text:
"""
prob = create_problem(background_σ, background_ϵ, desired_σ, desired_ϵ, incident_electric_field, frequency, resolution).\n
supplying the background and desired dielectrics, the incident electric field, frequency and resolution\n
this will create a problem tuple that can be solved by solve_problem().\n
"""
function create_problem(
    background_σ,
    background_ϵ,
    desired_σ,
    desired_ϵ,
    incident_electric_field,
    frequency,
    resolution
)

    constants = get_constants(frequency)

    p, q, r = size(background_σ);
    number_of_cells = [p-2, q-2, r-2];
    a = minimum(resolution)/2;
    scaling = -1im*constants.μ₀*constants.ω/(constants.kb^2);

    x_axis = 0f0:resolution[1]:p*resolution[1];
    y_axis = 0f0:resolution[2]:q*resolution[2];
    z_axis = 0f0:resolution[3]:r*resolution[3];

    background_dielectric = cell_to_yee_dielectric(background_σ, background_ϵ, x_axis, y_axis, z_axis);
    desired_dielectric = cell_to_yee_dielectric(desired_σ, desired_ϵ, x_axis, y_axis, z_axis);
    difference_in_dielectric, locations_of_change = get_update_maps(background_dielectric, desired_dielectric);

    #get the number of edges that need to be updated
    number_of_changes = _get_number_of_updated_edges(locations_of_change);

    #operators
    AToE = _electric_vector_potential_to_electric_field_operator(number_of_cells, resolution, constants.kb);
    Ig   = _create_Greens_functions_restriction_operators(number_of_cells);
    G    = create_Greens_function(number_of_cells, resolution, constants.kb);

    dielectric = _set_dielectric(background_dielectric, constants, number_of_cells);

    #allocate memory for inplace operations
    sources, vacuum_electric_field, vector_potential_1d, vector_potential_3d, fft_tmp, pfft, pifft = _allocate_memory_for_volume_integral_equation(number_of_cells, number_of_changes);
    x, p, r, rt, u, v, q, uq = _allocate_memory_cgs(number_of_cells);
    #get the sparse S matrix for each field component
    S = get_S_matrix(locations_of_change);
    C = get_C_matrix(difference_in_dielectric, constants);

    b, nb = get_b_vector(C, S, incident_electric_field)

    output_volume_integral_equation = zeros(ComplexF32, sum(number_of_changes));

    volume_integral_equation_map =  LinearMap{ComplexF32}(vecxyz-> _volume_integral_equation_wrapper_explicit_return(output_volume_integral_equation, sources, vacuum_electric_field, C, vecxyz, S, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq), sum(number_of_changes))
    scattered_field_map = LinearMap{ComplexF32}(current_density -> library_matrix_times_current_density!(current_density, sources, vacuum_electric_field, S, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq), sum(number_of_changes))

    return (g_map = volume_integral_equation_map, b_vector = b, norm_b = nb, scattered_electric_field_map = scattered_field_map)
end

"""
prob = create_problem(background_σ, background_ϵ, desired_σ, desired_ϵ, incident_electric_field, frequency, resolution).\n
supplying the background and desired dielectrics, the incident electric field, frequency and resolution\n
this will create a problem tuple that can be solved using a GPU by solve_problem().\n
"""
function create_problem_gpu(
    background_σ,
    background_ϵ,
    desired_σ,
    desired_ϵ,
    incident_electric_field,
    frequency,
    resolution
)

    constants = get_constants(frequency)

    p, q, r = size(background_σ);
    number_of_cells = [p-2, q-2, r-2];
    a = minimum(resolution)/2;

    scaling = -1im*constants.μ₀*constants.ω/(constants.kb^2);

    x_axis = 0f0:resolution[1]:p*resolution[1];
    y_axis = 0f0:resolution[2]:q*resolution[2];
    z_axis = 0f0:resolution[3]:r*resolution[3];

    background_dielectric = cell_to_yee_dielectric(background_σ, background_ϵ, x_axis, y_axis, z_axis);
    desired_dielectric = cell_to_yee_dielectric(desired_σ, desired_ϵ, x_axis, y_axis, z_axis);
    difference_in_dielectric, locations_of_change = get_update_maps(background_dielectric, desired_dielectric);

    #get the number of edges that need to be updated
    number_of_changes = _get_number_of_updated_edges(locations_of_change);

    #operators
    AToE = _electric_vector_potential_to_electric_field_operator(number_of_cells, resolution, constants.kb) |> CUDA.CUSPARSE.CuSparseMatrixCSC;
    Ig   = _create_Greens_functions_restriction_operators_gpu(number_of_cells);
    G    = create_Greens_function(number_of_cells, resolution, constants.kb) |> cu;

    dielectric = _set_dielectric_gpu(background_dielectric, constants, number_of_cells);

    #allocate memory for inplace operations
    sources, vacuum_electric_field, vector_potential_1d, vector_potential_3d, fft_tmp, pfft, pifft = _allocate_memory_for_volume_integral_equation_gpu(number_of_cells, number_of_changes);
    x, p, r, rt, u, v, q, uq = _allocate_memory_cgs_gpu(number_of_cells);
    #get the sparse S matrix for each field component
    S = get_S_matrix_gpu(locations_of_change);
    C = get_C_matrix(difference_in_dielectric, constants) |> cu;

    b, nb = get_b_vector(C, S, incident_electric_field)

    outputVIE = CUDA.zeros(ComplexF32, sum(number_of_changes));

    volume_integral_equation_map =  LinearMap{ComplexF32}(vecxyz-> _volume_integral_equation_wrapper_explicit_return_gpu(outputVIE, sources, vacuum_electric_field, C, vecxyz, S, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq), sum(number_of_changes))
    scattered_field_map = LinearMap{ComplexF32}(current_density -> library_matrix_times_current_density_gpu!(current_density, sources, vacuum_electric_field, S, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, scaling, pfft, pifft, x, p, r, rt, u, v, q, uq), sum(number_of_changes))

    return (g_map = volume_integral_equation_map, b_vector = b, norm_b = nb, scattered_electric_field_map = scattered_field_map)
end

"""
scattered_current_density, scattered_electric_field = solve_problem(problem; verbose = true, max_iterations = 100, restart = 50, tol = 1f-4).\n
the problem is created by calling the create_problem function.\n
the problem is solved using IterativeSolvers gmres implementation.\n
it can also be solved by other solvers if you believe they are more accurate or faster.\n
personally, I found this solver to be both fast, stable, and finding the correct solution.\n
Some standard arguments are supplied, but can be tweaked according to preference/complexity of the problem.\n
"""
function solve_problem(
    problem;
    verbose = true,
    max_iterations = 100,
    restart = 50,
    tol = 1f-4
)
    # solve the matrix inverse to obtain the scattered current density
    scattered_current_density = IterativeSolvers.gmres(problem.g_map, problem.b_vector, verbose=verbose, initially_zero=true, reltol=tol, restart = restart, maxiter=max_iterations);
    # the vector b is normalized to a norm of 1, thus we rescale the solution since the problem is linear.
    scattered_current_density .*= problem.norm_b
    # we calculate the scattered electric field that is created by the dielectric update that is performed
    scattered_electric_field  = problem.scattered_electric_field_map*scattered_current_density
    return scattered_current_density, scattered_electric_field
end

"""
electric_field = calculate_electric_field_vacuum_to_dielectric(σ, ϵr, resolution, source, frequency)\n
this will calculate an electric field distribution using the VIE method.\n
Warning: this is not stable for very high conductivities or high permittivities.\n
"""
function calculate_electric_field_vacuum_to_dielectric(
    σ,
    ϵr,
    resolution,
    source,
    frequency
)
    constants = get_constants(frequency);
    np, nq, nr   = size(σ);
    number_of_cells  = [np-2, nq-2, nr-2];
    a       = minimum(resolution)/2f0;
    divωϵim = 1/(1im*constants.ω*constants.ϵ₀);

    x_axis = 0:resolution[1]:np*resolution[1]
    y_axis = 0:resolution[2]:nq*resolution[2]
    z_axis = 0:resolution[3]:nr*resolution[3]

    yee_dielectric = cell_to_yee_dielectric(σ, ϵr, x_axis, y_axis, z_axis);

    dielectric = _set_dielectric(yee_dielectric, constants, number_of_cells);
    # operators
    AToE = _electric_vector_potential_to_electric_field_operator(number_of_cells, resolution, constants.kb);
    Ig   = _create_Greens_functions_restriction_operators(number_of_cells);
    G = create_Greens_function(number_of_cells, resolution, constants.kb);
    #allocate memory for the VIE method
    _, _, vector_potential_1d, vector_potential_3d, fft_tmp, pfft, pifft = _allocate_memory_for_volume_integral_equation(number_of_cells, [1,0,0]);

    incident_electric_field = copy(source);
    #from source in E incident to the actual incident electric field
    source_to_incident_field!(incident_electric_field, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, divωϵim);

    #create map for VIE
    v_out = zeros(ComplexF32, length(incident_electric_field))
    volume_integral_equation_map = LinearMap{ComplexF32}(x_in -> total_minus_scattered_electric_field_map!(v_out, x_in, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, one(ComplexF32)), length(v_out))

    #solve for the total electric field
    total_electric_field = IterativeSolvers.gmres(volume_integral_equation_map, incident_electric_field, verbose=true, restart = 50, maxiter=50)

    return total_electric_field
end


text:
using Cropbox
using CropRootBox
using Test

root_maize = @config(
    :RootArchitecture => :maxB => 5,
    :BaseRoot => :T => [
        # P F S    T
          0 1 0 #= 0 =#; # P
          0 0 1 #= 0 =#; # F
          0 0 0 #= 0 =#; # S
       #= 0 0 0    0   ; # T =#
    ],
    :PrimaryRoot => (;
        lb = 0.1 ± 0.01,
        la = 18.0 ± 1.8,
        ln = 0.6 ± 0.06,
        lmax = 89.7 ± 7.4,
        r = 6.0 ± 0.6,
        Δx = 0.5,
        σ = 10,
        θ = 80 ± 8,
        N = 1.5,
        a = 0.04 ± 0.004,
        color = CropRootBox.RGBA(1, 0, 0, 1),
    ),
    :FirstOrderLateralRoot => (;
        lb = 0.2 ± 0.04,
        la = 0.4 ± 0.04,
        ln = 0.4 ± 0.03,
        lmax = 0.6 ± 1.6,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 15,
        N = 1,
        a = 0.03 ± 0.003,
        color = CropRootBox.RGBA(0, 1, 0, 1),
    ),
    :SecondOrderLateralRoot => (;
        lb = 0,
        la = 0.4 ± 0.02,
        ln = 0,
        lmax = 0.4,
        r = 2.0 ± 0.2,
        Δx = 0.1,
        σ = 20,
        θ = 70 ± 10,
        N = 2,
        a = 0.02 ± 0.002,
        color = CropRootBox.RGBA(0, 0, 1, 1),
    )
)

@testset "root" begin
    b = instance(CropRootBox.Pot)
    s = instance(CropRootBox.RootArchitecture; config = root_maize, options = (; box = b), seed = 0)
    r = simulate!(s, stop = 100u"d")
    @test r.time[end] == 100u"d"
    # using GLMakie
    # scn = CropRootBox.render(s)
    # GLMakie.save("root_maize.png", scn)
    CropRootBox.writevtk(tempname(), s)
    # CropRootBox.writepvd(tempname(), CropRootBox.RootArchitecture, config = root_maize, stop = 50)
    CropRootBox.writestl(tempname(), s)
end


text:
"""
    gevfitbayes(..., niter::Int=5000, warmup::Int=2000)

Generate a sample from the GEV parameters' posterior distribution.

# Arguments
- `niter::Int = 5000`: The total number of MCMC iterations.
- `warmup::Int = 2000`: The number of warmup iterations (burn-in).

# Implementation

The function uses the No-U-Turn Sampler (NUTS; [Hoffman and Gelman, 2014](http://jmlr.org/papers/v15/hoffman14a.html))
implemented in the [Mamba.jl](https://mambajl.readthedocs.io/en/latest/index.html)
package to generate a random sample from the posterior distribution.

Currently, only the improper uniform prior is implemented, *i.e.*
```math
f_{(β₁,β₂,β₃)}(β₁,β₂,β₃) ∝ 1,
```
where
```math
μ = X₁ × β₁,
```
```math
ϕ = X₂ × β₂,
```
```math
ξ = X₃ × β₃.
```
In the stationary case, this improper prior yields to a proper posterior if the
sample size is larger than 3 ([Northrop and Attalides, 2016](https://www.jstor.org/stable/24721296?seq=1)).

The covariates are standardized before estimating the parameters to help fit the
 model. They are transformed back on their original scales before returning the
 fitted model.

See also [`gevfitbayes`](@ref) for the other methods, [`gevfitpwm`](@ref), [`gevfit`](@ref) and [`BlockMaxima`](@ref).

# References

Hoffman M. D. and Gelman A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. *Journal of Machine Learning Research*, 15:1593–1623.

Paul J. Northrop P. J. and Attalides N. (2016). Posterior propriety in Bayesian extreme value analyses using reference priors. *Statistica Sinica*, 26:721-743.

"""
function gevfitbayes end

"""
    gevfitbayes(y,
        locationcov = Vector{Variable}(),
        logscalecov = Vector{Variable}(),
        shapecov = Vector{Variable}(),
        niter::Int=5000,
        warmup::Int=2000
        )

Generate a sample from the GEV parameters' posterior distribution.

# Arguments

- `y::Vector{<:Real}`: The vector of block maxima.
- `locationcov::Vector{<:DataItem} = Vector{Variable}()`: The covariates of the location parameter.
- `logscalecov::Vector{<:DataItem} = Vector{Variable}()`: The covariates of the log-scale parameter.
- `shapecov::Vector{<:DataItem} = Vector{Variable}()`: The covariates of the shape parameter.

"""
function gevfitbayes(y::Vector{<:Real};
    locationcov::Vector{<:DataItem} = Vector{Variable}(),
    logscalecov::Vector{<:DataItem} = Vector{Variable}(),
    shapecov::Vector{<:DataItem} = Vector{Variable}(),
    niter::Int=5000, warmup::Int=2000)::BayesianEVA

    locationcovstd = standardize.(locationcov)
    logscalecovstd = standardize.(logscalecov)
    shapecovstd = standardize.(shapecov)

    model = BlockMaxima(Variable("y", y), locationcov = locationcovstd, logscalecov = logscalecovstd, shapecov = shapecovstd)

    fittedmodel = fitbayes(model, niter=niter, warmup=warmup)

    return transform(fittedmodel)

end

"""
    gevfitbayes(df::DataFrame,
        datacol::Symbol,
        locationcovid = Vector{Symbol}(),
        logscalecovid = Vector{Symbol}(),
        shapecovid = Vector{Symbol}(),
        niter::Int=5000,
        warmup::Int=2000)

Generate a sample from the GEV parameters' posterior distribution.

# Arguments
- `df::DataFrame`: The dataframe containing the data.
- `datacol::Symbol`: The symbol of the column of `df` containing the block maxima data.
- `locationcovid::Vector{Symbol} = Vector{Symbol}()`: The symbols of the columns of `df` containing the covariates of the location parameter.
- `logscalecovid::Vector{Symbol} = Vector{Symbol}()`: The symbols of the columns of `df` containing the covariates of the log-scale parameter.
- `shapecovid::Vector{Symbol} = Vector{Symbol}()`: The symbols of the columns of `df` containing the covariates of the shape parameter.
"""
function gevfitbayes(df::DataFrame, datacol::Symbol;
    locationcovid::Vector{Symbol}=Symbol[],
    logscalecovid::Vector{Symbol}=Symbol[],
    shapecovid::Vector{Symbol}=Symbol[],
    niter::Int=5000, warmup::Int=2000)::BayesianEVA

    locationcovstd = standardize.(buildVariables(df, locationcovid))
    logscalecovstd = standardize.(buildVariables(df, logscalecovid))
    shapecovstd = standardize.(buildVariables(df, shapecovid))

    model = BlockMaxima(Variable(string(datacol), df[:, datacol]), locationcov = locationcovstd, logscalecov = logscalecovstd, shapecov = shapecovstd)

    fittedmodel = fitbayes(model, niter=niter, warmup=warmup)

    return transform(fittedmodel)

end

"""
    gevfitbayes(model::BlockMaxima;
        niter::Int=5000,
        warmup::Int=2000)

Generate a sample from the `BlockMaxima` model parameters' posterior distribution.
"""
function gevfitbayes(model::BlockMaxima; niter::Int=5000, warmup::Int=2000)::BayesianEVA

    return fitbayes(model, niter=niter, warmup=warmup)

end


text:
using LinearCovarianceModels
using LinearAlgebra, Test
using HomotopyContinuation

const HC = HomotopyContinuation
const LC = LinearCovarianceModels

@testset "LinearCovariance.jl" begin
    @testset "vec to sym and back" begin
        v = [1,2,3, 4, 5, 6]
        @test vec_to_sym(v) == [1 2 3; 2 4 5; 3 5 6]
        @test sym_to_vec(vec_to_sym(v)) == v
    end

    @testset "LCModels" begin
        A = toeplitz(3)
        @test A isa LCModel
        x = variables(vec(A.Σ))
        @test A.Σ == [x[1] x[2] x[3]
                    x[2] x[1] x[2]
                    x[3] x[2] x[1]]
        @test dim(A) == 3
        T = tree(4, "{1,2},{1,2,3}")
        @test dim(T) == 7

        Ts = trees(5)
        @test all(isa.(last.(Ts), LCModel))

        D = generic_diagonal(6, 3)
        @test dim(D) == 3
        @test_throws ArgumentError generic_diagonal(6, 0)
        @test_throws ArgumentError generic_diagonal(6, 7)

        D = generic_subspace(6, 4)
        @test dim(D) == 4
        @test_throws ArgumentError generic_subspace(6, binomial(6+1,2)+1)
        @test_throws ArgumentError generic_subspace(6, 0)

        # @test_throws ArgumentError LCModel(toeplitz(3).Σ .^2)

        # throw for non-symmetric input
        @test_throws ArgumentError LCModel([x[1] x[1]; x[2] x[1]])

        # handle special matrix types
        @var θ[1:7]
        M = LCModel(SymTridiagonal(θ[1:4],θ[5:7]))
        @test dim(M) == 7
    end

    @testset "mle system" begin
        @var x y z
        Σ = [x y; y z]
        F = mle_system(LCModel(Σ))
        @test F isa HC.System
    end

    @testset "ml_degree_witness" begin
        Σ = toeplitz(3)
        W = ml_degree_witness(Σ)
        @test model(W) isa LCModel
        @test is_dual(W) == false
        @test ml_degree(W) == 3
        @test length(solutions(W)) == 3
        @test dim(model(W)) == 3
        @test parameters(W) isa AbstractVector
        @test verify(W)
        @test verify(W; trace_tol=1e-6)
    end

    @testset "solve" begin
        S = [4/5 -9/5 -1/25
            -9/5 79/16 25/24
            -1/25 25/24 17/16]
        Σ = toeplitz(3)
        W = ml_degree_witness(Σ)
        crits = critical_points(W, S)
        @test length(crits) == 3
        sort!(crits; by=s -> s[2], rev=true)
        p1,p2,p3 = last.(crits)
        @test p1 == :global_maximum
        @test p2 == :local_maximum
        @test p3 == :saddle_point
        @test crits[1][1] ≈ mle(W, S) atol=1e-8
    end
end


text:
using Base.Meta: quot, show_sexpr, isexpr

export @sexpr, @methods, @methodswith, @dump, @macroexpand, @esc

macro sexpr(expr)
    :((show_sexpr($(quot(expr))); println()))
end

macro methods(func::Symbol)
    :(methods($func))
end

# Simplify life when calling methodswith. I always do it wrong.
Base.methodswith(x, showparents::Bool=false) = methodswith(typeof(x), showparents)

macro methodswith(expr)
    :(methodswith($expr))
end

macro dump(expr)
    :(dump($(quot(expr))))
end

macro macroexpand(expr)
    :(macroexpand($(quot(expr))))
end

macro esc(sym)
    isexpr(sym, :quote) || throw(ArgumentError("expected quoted symbol"))
    esc(sym.args[1])
end


text:
#################### Model Simulation ####################

function gettune(m::Model, block::Integer)
  block == 0 && return gettune(m)
  m.samplers[block].tune
end

function gettune(m::Model)
  Any[gettune(m, i) for i in 1:length(m.samplers)]
end

function settune!(m::Model, tune, block::Integer)
  block == 0 && return settune!(m, tune)
  m.samplers[block].tune = tune
end

function settune!(m::Model, tune::Vector{Any})
  nsamplers = length(m.samplers)
  ntune = length(tune)
  nsamplers == ntune ||
    throw(DimensionMismatch(
      "tried to assign $ntune tune elements to $nsamplers samplers"
    ))

  for i in 1:nsamplers
    settune!(m, tune[i], i)
  end
end

function gradlogpdf(m::Model, block::Integer=0, transform::Bool=false;
                    dtype::Symbol=:forward)
  x0 = unlist(m, block, transform)
  value = gradlogpdf!(m, x0, block, transform, dtype=dtype)
  relist!(m, x0, block, transform)
  value
end

function gradlogpdf{T<:Real}(m::Model, x::AbstractVector{T}, block::Integer=0,
                             transform::Bool=false; dtype::Symbol=:forward)
  x0 = unlist(m, block)
  value = gradlogpdf!(m, x, block, transform, dtype=dtype)
  relist!(m, x0, block)
  value
end

function gradlogpdf!{T<:Real}(m::Model, x::AbstractVector{T}, block::Integer=0,
                              transform::Bool=false; dtype::Symbol=:forward)
  f = x -> logpdf!(m, x, block, transform)
  gradient(f, convert(Vector{T}, x), dtype)
end

function logpdf(m::Model, block::Integer=0, transform::Bool=false)
  params = keys(m, :block, block)
  targets = keys(m, :target, block)
  logpdf(m, params, transform) + logpdf(m, setdiff(targets, params))
end

function logpdf(m::Model, nodekeys::Vector{Symbol}, transform::Bool=false)
  lp = 0.0
  for key in nodekeys
    lp += logpdf(m[key], transform)
    isfinite(lp) || break
  end
  lp
end

function logpdf{T<:Real}(m::Model, x::AbstractArray{T}, block::Integer=0,
                         transform::Bool=false)
  x0 = unlist(m, block)
  lp = logpdf!(m, x, block, transform)
  relist!(m, x0, block)
  lp
end

function logpdf!{T<:Real}(m::Model, x::AbstractArray{T}, block::Integer=0,
                          transform::Bool=false)
  params = keys(m, :block, block)
  targets = keys(m, :target, block)
  m[params] = relist(m, x, params, transform)
  lp = logpdf(m, setdiff(params, targets), transform)
  for key in targets
    isfinite(lp) || break
    node = m[key]
    update!(node, m)
    lp += key in params ? logpdf(node, transform) : logpdf(node)
  end
  lp
end

function sample!(m::Model, block::Integer=0)
  m.iter += 1
  isoneblock = block != 0
  blocks = isoneblock ? block : 1:length(m.samplers)
  for b in blocks
    sampler = m.samplers[b]
    value = sampler.eval(m, b)
    if value != nothing
      m[sampler.params] = value
      update!(m, b)
    end
  end
  m.iter -= isoneblock
  m
end

function unlist(m::Model, block::Integer=0, transform::Bool=false)
  unlist(m, keys(m, :block, block), transform)
end

function unlist(m::Model, monitoronly::Bool)
  f = function(key)
    node = m[key]
    lvalue = unlist(node)
    monitoronly ? lvalue[node.monitor] : lvalue
  end
  vcat(map(f, keys(m, :dependent))...)
end

function unlist(m::Model, nodekeys::Vector{Symbol}, transform::Bool=false)
  vcat(map(key -> unlist(m[key], transform), nodekeys)...)
end

function relist{T<:Real}(m::Model, x::AbstractArray{T}, block::Integer=0,
                         transform::Bool=false)
  relist(m, x, keys(m, :block, block), transform)
end

function relist{T<:Real}(m::Model, x::AbstractArray{T},
                         nodekeys::Vector{Symbol}, transform::Bool=false)
  values = Dict{Symbol,Any}()
  N = length(x)
  offset = 0
  for key in nodekeys
    value, n = relistlength(m[key], sub(x, (offset + 1):N), transform)
    values[key] = value
    offset += n
  end
  offset == length(x) ||
    throw(ArgumentError("incompatible number of values to put in nodes"))
  values
end

function relist!{T<:Real}(m::Model, x::AbstractArray{T}, block::Integer=0,
                 transform::Bool=false)
  nodekeys = keys(m, :block, block)
  values = relist(m, x, nodekeys, transform)
  for key in nodekeys
    m[key].value = values[key]
  end
  update!(m, block)
end

function relist!{T<:Real}(m::Model, x::AbstractArray{T}, nodekey::Symbol,
                          transform::Bool=false)
  node = m[nodekey]
  m[nodekey] = relist(node, x, transform)
  update!(m, node.targets)
end

function update!(m::Model, block::Integer=0)
  nodekeys = block == 0 ? keys(m, :dependent) : m.samplers[block].targets
  update!(m, nodekeys)
end

function update!(m::Model, nodekeys::Vector{Symbol})
  for key in nodekeys
    update!(m[key], m)
  end
  m
end


text:
export exactpruning!

# Exact pruning
"""
    exactpruning!(man::AbstractCutPruner, solver::MathProgBase.AbstractMathProgSolver;
                  ub=Inf, lb=-Inf, epsilon=1e-5)

Remove dominated cuts in CutPruner `man`.

We use a LP solver to determine whether a cut is dominated or not.

# Arguments
* `man::AbstractCutPruner`
    Cut pruner where to remove cuts
* `solver`
    Solver used to solve LP
* `ub::Union{Float64, Vector{Float64}}`
    State x upper bound
* `lb::Union{Float64, Vector{Float64}}`
    State x lower bound
* `epsilon::Float64`
    Pruning's tolerance
"""
function exactpruning!(man::AbstractCutPruner, solver::MathProgBase.AbstractMathProgSolver; ub=Inf, lb=-Inf, epsilon=1e-5)
    K = getdominated(man.A, man.b, man.islb, man.isfun, solver, lb, ub, epsilon)
    removecuts!(man, K)
end

"""Return dominated cuts."""
function getdominated(A, b, islb, isfun, solver, lb, ub, epsilon)
    red = Int[]
    if size(A, 1) == 1
        return red
    end
    for i in 1:size(A, 1)
        if isdominated(A, b, islb, isfun, i, solver, lb, ub, epsilon)
            push!(red, i)
        end
    end
    red
end

"""State whether a cut is dominated with a tolerance epsilon."""
function isdominated(A, b, islb, isfun, k, solver, lb, ub, epsilon)
    # we use MathProgBase to solve the test
    # For instance, if islb & isfun, the LP becomes:
    # min - y
    #     x ∈ R^n, y ∈ R
    #     y <= (A[k, :] - A[i, :])*x + b[k] - b[i]     ∀ i != k

    # we formulate an equivalent problem as
    # min c'*z
    # s.t H*z <= h

    # get problem's dimension
    ncuts, nx = size(A)
    # allocate arrays
    h = zeros(ncuts - 1)
    H = zeros(ncuts - 1, nx + 1)
    c = zeros(nx + 1)
    c[1] = (islb) ? -1 : 1

    λk = @view A[k, :]

    ic = 0
    for ix in 1:ncuts
        if ix != k
            ic += 1
            δb = b[k] - b[ix]
            @inbounds h[ic] = ((isfun&islb) || (~isfun&~islb)) ? δb : -δb
            @inbounds H[ic, 1] = (islb) ? 1 : -1

            for jx in 1:nx
                dl = A[ix, jx] - λk[jx]
                @inbounds H[ic, jx+1] = (islb) ? dl : -dl
            end
        end
    end

    # update lower and upper bound if Vector
    lbx = isa(lb, Vector) ? vcat(-Inf, lb) : lb
    ubx = isa(ub, Vector) ? vcat( Inf, ub) : ub

    # solve the LP with MathProgBase
    res = linprog(c, H, -Inf, h, lbx, ubx, solver)
    if res.status == :Optimal
        res = res.objval
        return (islb) ? -res < epsilon : res > -epsilon
    else
        return false
    end
end


text:
using Test

@testset "Quite partial to it really" begin
  add3 = @p(3 + _)
  @test add3(7) == 10
  add4r = @p(_ + 4)
  @test add4r(7) == 11

  lt10 = @partial _ < 10
  @test lt10(7) == true
  @test lt10(27) == false

  gt10 = @p(10 < _)
  @test gt10(17) == true
  @test gt10(7) == false

  fn_of7_and_10 = @p _(7,10)
  @test fn_of7_and_10(*) == 70
  @test fn_of7_and_10(+) == 17
  @test fn_of7_and_10(-) == -3
  @test fn_of7_and_10(/) == 0.7
  @test fn_of7_and_10(%) == 7

  @test map(@p(_(7,10)),[*,+,-,/,%]) == [70,17,-3,0.7,7]

  v1 = 3
  v2 = 5
  fn_of_v1_and_v2 = @p _(v1,v2)
  @test fn_of_v1_and_v2(*) == 15
  @test fn_of_v1_and_v2(+) == 8
  @test fn_of_v1_and_v2(-) == -2
  @test fn_of_v1_and_v2(/) == 0.6
  @test fn_of_v1_and_v2(%) == 3

end

@testset "partial multi sub" begin
  d = 10
  zeros_dxN = @p zeros(_,d,_)
  @test zeros_dxN(Int,3) == zeros(Int,10,3)
  T = Float64
  d2 = 12
  zerosF64 = @p zeros(T,_,_)
  @test zerosF64(d2,4) == zeros(Float64,d2,4)
end

@testset "partial can handle nested syms" begin
  len_gt_10 = @partial(length(_) > 10)
  @test len_gt_10("fred") == false
  @test len_gt_10("freddy mercury") == true

  fred = Array{Real}[]
  push!(fred, [1,2,3])
  push!(fred, [4.0,5.0,6.0])
  fred_sub_ind = @p(fred[_][_])
  @test fred_sub_ind(1,2) == 2
  @test fred_sub_ind(1,3) == 3
  @test fred_sub_ind(2,1) == 4.0
  @test fred_sub_ind(2,3) == 6.0

  len_gt_x = @partial(length(_) > _)
  @test len_gt_x("fred",2) == true
  @test len_gt_x("fred",10) == false
end

@testset "partial indexed sub" begin
  d = 10
  zeros_dxN = @p zeros(Int,_1,_1)
  @test zeros_dxN(3) == zeros(Int,3,3)
  zeros_dxN = @p zeros(Int,_1,_)
  @test zeros_dxN(5) == zeros(Int,5,5)
  T = Float64
  d2 = 12
  zerosF64 = @p zeros(T,_,_1)
  @test zerosF64(d2) == zeros(Float64,12,12)
end

@testset "partial chained" begin
  # unary fn
  x10plus3 = @p _*10 _+3
  for (inp, outp) in zip((2,3,4), (23,33,43))
    @test x10plus3(inp) == outp
  end
  # binary fn
  b_minus_a_x3 = @p _2-_1 _*3
  for (inp, outp) in zip([(2,4), (7,3), (9,2)], (6, -12, -21))
    @test b_minus_a_x3(inp...) == outp
  end
end


text:
function endforces!(node::Int64, params::Parameters)
  if node == params.noofnodes
    l = params.nodes[node].z - params.nodes[node-1].z
    # Average rate of change in y direction
    θ = -(params.nodes[node].y - params.nodes[node-1].y) / l
    # Average rate of change in x direction
    γ = (params.nodes[node].x - params.nodes[node-1].x) / l

    params.f[1, 1, node] += params.fn[params.noofnodes] * γ
    params.f[2, 1, node] -= params.fn[params.noofnodes] * θ
    params.f[3, 1, node] += params.fn[params.noofnodes]
  end

  if node <= 3
    params.del[4, 3, node] = sqrt(params.del[4, 2, node]^2 + params.del[4, 3, node]^2)
    params.del[4, 2, node] = 0.0
  end
end


text:
# This file was generated, do not modify it. # hide
x = 3.141592
y = sin(x)
@show y


text:
# Base.null_safe_op is not defined on v0.5
if VERSION < v"0.6.0-dev.0"
    # TODO: Find parsimonious way to support null_safe_op optimizations on 0.5
    null_safe_op(f, xs...) = false
else
    using Base: null_safe_op
end

lift(f, x) = f(x)

abstract MaybeNullable
immutable YesNullable <: MaybeNullable end
immutable NoNullable <: MaybeNullable end
Base.@pure any_nullable(x,y) = ifelse((x<:Nullable) | (y<:Nullable), YesNullable(), NoNullable())
Base.@pure any_nullable(xs) = ifelse(mapreduce(x->x<:Nullable, |, xs.parameters), YesNullable(), NoNullable())

# # back-up in case @pure gives us grief
# any_nullable(x::Nullable, xs...) = true
# any_nullable(x, xs...) = any_nullables(xs...)
# any_nullable(x) = false
# any_nullable() = false

@inline lift(f, x, y) = lift(f, x, y, any_nullable(typeof(x), typeof(y)))
@inline lift(f, xs...) = lift(f, xs, any_nullable(typeof(xs)))
@inline lift(f, x, y, ::NoNullable) = f(x, y)
@inline lift(f, xs, ::NoNullable) = f(xs...)

##############################################################################
##
## Standard lifting semantics
##
## For a function call f(xs...), return null if any x in xs is null;
## otherwise, return f applied to values of xs.
##
##############################################################################

@inline function lift(f, x::Nullable)
    if null_safe_op(f, typeof(x))
        return @compat Nullable(f(x.value), !isnull(x))
    else
        U = Core.Inference.return_type(f, Tuple{eltype(typeof(x))})
        if isnull(x)
            return Nullable{U}()
        else
            return Nullable(f(unsafe_get(x)))
        end
    end
end

@inline function lift(f, x1, x2, ::YesNullable)
    if null_safe_op(f, typeof(x1), typeof(x2))
        return @compat Nullable(
            f(x1.value, x2.value), !(isnull(x1) | isnull(x2)))
    else
        U = Core.Inference.return_type(
            f, Tuple{eltype(typeof(x1)), eltype(typeof(x2))})
        if isnull(x1) | isnull(x2)
            return Nullable{U}()
        else
            return Nullable(f(unsafe_get(x1), unsafe_get(x2)))
        end
    end
end

@inline function lift(f, xs, ::YesNullable)
    if null_safe_op(f, map(typeof, xs)...)
        return @compat Nullable(
            f(map(unsafe_get, xs)...), !(mapreduce(isnull, |, xs)))
    else
        U = Core.Inference.return_type(
            f, Tuple{map(x->eltype(typeof(x)), xs)...})
        if hasnulls(xs)
            return Nullable{U}()
        else
            return Nullable(f(map(unsafe_get, xs)...))
        end
    end
end

##############################################################################
##
## Non-standard lifting semantics
##
##############################################################################

# three-valued logic implementation
@inline function lift(::typeof(&), x, y, ::YesNullable)::Nullable{Bool}
    return ifelse( isnull(x),
        ifelse( isnull(y),
            Nullable{Bool}(),                       # x, y null
            ifelse( unsafe_get(y),
                Nullable{Bool}(),                   # x null, y == true
                Nullable(false)                     # x null, y == false
            )
        ),
        ifelse( isnull(y),
            ifelse( unsafe_get(x),
                Nullable{Bool}(),                   # x == true, y null
                Nullable(false)                     # x == false, y null
            ),
            Nullable(unsafe_get(x) & unsafe_get(y)) # x, y not null
        )
    )
end

# three-valued logic implementation
@inline function lift(::typeof(|), x, y, ::YesNullable)::Nullable{Bool}
    return ifelse( isnull(x),
        ifelse( isnull(y),
            Nullable{Bool}(),                       # x, y null
            ifelse( unsafe_get(y),
                Nullable(true),                     # x null, y == true
                Nullable{Bool}()                    # x null, y == false
            )
        ),
        ifelse( isnull(y),
            ifelse( unsafe_get(x),
                Nullable(true),                     # x == true, y null
                Nullable{Bool}()                    # x == false, y null
            ),
            Nullable(unsafe_get(x) | unsafe_get(y)) # x, y not null
        )
    )
end

# TODO: Decide on semantics for isequal and uncomment the following
# @inline function lift(::typeof(isequal), x, y)
#     return ifelse( isnull(x),
#         ifelse( isnull(y),
#             true,                                   # x, y null
#             false                                   # x null, y not null
#         ),
#         ifelse( isnull(y),
#             false,                                   # x not null, y null
#             isequal(unsafe_get(x), unsafe_get(y))    # x, y not null
#         )
#     )
# end

@inline function lift(::typeof(isless), x, y, ::YesNullable)::Bool
    if null_safe_op(isless, typeof(x), typeof(y))
        return ifelse( isnull(x),
            false,                                      # x null
            ifelse( isnull(y),
                true,                                   # x not null, y null
                isless(unsafe_get(x), unsafe_get(y))    # x, y not null
            )
        )
    else
        return  isnull(x) ? false :
                isnull(y) ? true  : isless(unsafe_get(x), unsafe_get(y))
    end
end

@inline lift(::typeof(isnull), x::Nullable) = isnull(x)
@inline lift(::typeof(isnull), x, ::MaybeNullable) = isnull(x)
@inline lift(::typeof(get), x::Nullable) = get(x)
@inline lift(::typeof(get), x::Nullable, y) = get(x, y)


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

include("AppDemo.jl")
using .AppDemo

include("../Renderer.jl")
using .Renderer

## create store
store = create_store(AppDemo.reducer, AppDemo.IMGUI_DEMO_STATE)

## draw UI
Renderer.render(
    () -> AppDemo.ImGui_Demo(store),
    width = 1280,
    height = 720,
    title = "IMGUI Demo",
)


text:
# using Dates

@testset "Gar.jl" begin
    filename = "data\\test.gar"
    dc = DDR2import.Gar.read(filename)
    @test dc["000EG"][5,1] == 57.0
    @test dc["000EG"][6,2] ≈ -13.8213888889 atol = 0.00001
    #LAST
    @test dc["000EU"][5,1] ≈ 39.6677777778 atol = 0.00001
end


text:
using Visualisations
using Documenter

DocMeta.setdocmeta!(Visualisations, :DocTestSetup, :(using Visualisations); recursive=true)

makedocs(;
    modules=[Visualisations],
    authors="Graham Stark",
    repo="https://github.com/grahamstark/Visualisations.jl/blob/{commit}{path}#{line}",
    sitename="Visualisations.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://grahamstark.github.io/Visualisations.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/grahamstark/Visualisations.jl",
)


text:
function wk3odes(t::Float64,y::Vector{Float64},dy::Vector{Float64},mparams::WK3.ModelParams,k::Float64)
    # ventricular elastance and pressure

    E,dE = elastancefn(t*ts,mparams,k);
    Pv = E.*(y[2].*Vs.-mparams.V0)./Ps;

    # derivatives
    # println("ODEs, first loop:")
    if y[3] != 0
        dy[1] = ts/Ps*((1+mparams.Zc/mparams.R)*mparams.C)^-1*((1+mparams.Zc/mparams.R)*Qs*y[3]
            .+ mparams.C*mparams.Zc/mparams.R*(-E*Qs*y[3].+dE*(Vs*y[2]-mparams.V0)) .- Ps*y[1]/mparams.R);
    else
        dy[1] = ts/Ps*((1+mparams.Zc/mparams.R)*mparams.C)^-1.*(-Ps.*y[1]./mparams.R);
    end
    # println("ODEs, second loop:")
    if Pv > y[1] || (Pv < y[1] && y[3] > 0)
        dy[2] = -ts/Vs*Qs*y[3];
        dy[3] = ts/(Qs*(mparams.Rav+mparams.Zc))*(-E*Qs*y[3] .+ dE*(Vs*y[2].-mparams.V0) .-
            (1+mparams.Zc/mparams.R)*Qs/mparams.C*y[3] .+ Ps/(mparams.R*mparams.C)*y[1]);
    elseif Pv < y[1] && y[3] <= 0
        dy[2] = 0.;
        dy[3] = 0.;
    end

    # for i = 1:length(dy)
    #     if dy[i] > 1e6
    #         warn("dy/dt above tolerance at t = $t. dE/dt = $dE.")
    #     end
    # end
end


text:
#ENV["GKSwstype"]="svg"

function get_dicts(time_db, feed_db)
    allids = union(unique(time_db[1][:,1]), unique(feed_db[1][:,1]))
    tdict = Dict(map(p -> p => [], allids))
    fdict = Dict(map(p -> p => [], allids))
    for i in 1:size(time_db[1])[1]
        r = time_db[1][i,:]
        push!(tdict[r[1]], r[5])
    end
    for i in 1:size(feed_db[1])[1]
        r = feed_db[1][i,:]
        push!(fdict[r[1]], (r[2:end]))
    end
    return tdict, fdict
end

function send_cal_msg(id::AbstractString)
    domain = Ft2k.conf["Watchman"]["domain"]
    run(`mosquitto_pub --cafile /etc/mosquitto/ca_certificates/ca.crt -h $domain -t deltacal -m $id,no,88 -p 8883`)
    println("sent mqtt msg to deltacal for $id")
end

function send_deltaorder_msg(id, scheme)
    domain = Ft2k.conf["Watchman"]["domain"]
    println("sending deltaorder for $id for $scheme")
    feedingtimes = []
    times = []
    if scheme == "7times"
        feedingtimes = [7.2 9 11 13 15 17 18.8]
    elseif scheme == "12times"
        feedingtimes = collect(7.5:1:18.5)
    elseif scheme == "manual"
        feedingtimes = [7.4 8.4 13.4 14.4]
    elseif scheme == "DR"
        feedingtimes = [7.4 7.9 8.4]
    elseif scheme == "comp7"
        feedingtimes = [7.2 7.5 7.8 8.0 8.2 8.4 8.6]
    end
    for t in feedingtimes
        push!(times, Int(round(t*60*60+2)))
    end
    times = [x+rand(-240:240) for x in times]
    joinedtimes = join(times, ",")
    # actual publication of feedings
    try
        run(`mosquitto_pub --cafile /etc/mosquitto/ca_certificates/ca.crt -h $domain -p 8883 -t deltaorder -m $id:$joinedtimes`)
        println("sent mqtt msg to deltaorder for $id")
        println("feedings now: $joinedtimes")
    catch
        println("ERROR: could not update commands")
    end
end
# map(p -> send_deltaorder_msg(p, "7times"), lifespancohort)

function pullfeedingcountsbyday(fdict)
    daysoi = unique(map(p -> (p[5]), fdict))
    countdict = Dict(map(p -> p[1] => p[2], [(i,0) for i in daysoi]))
    [countdict[(i[5])] += 1 for i in fdict]
    lcd = length(countdict)
    feed_day = Array{Date,1}(undef, lcd)
    feed_num =  Array{Int64,1}(undef, lcd)
    i = 1
    for (k,v) in countdict
        feed_day[i] = k
        feed_num[i] = v
        i += 1
    end
    return feed_day, feed_num
end

function setup_app_comps(tdict, fdict, allids, sqlite_db)
    charray = []
    cha_dict = Dict()
    cohort_dict = Dict("cohort1" => "cohort1", "cohort2"=> "cohort2")
    for i in 1:size(sqlite_db[1])[1]
        r = sqlite_db[1][i,:]
        id = r[1]
        if haskey(tdict, id) && length(tdict[id])>0
            feedings = r[2]
            calibrated =r[3]
            num_feedings = (length(split(feedings, ",")))
            nowtime = now(Dates.UTC)-Hour(8)
            last_checkedin = (Int(round((Dates.datetime2unix(nowtime) - maximum(tdict[id]) )/60)))
            if length(fdict[id]) > 0
                fdlookup = sort(fdict[id], by = p -> p[5], rev = true)[1]
                last_fed_diff = (fdlookup[3] - fdlookup[2])
            else
                last_fed_diff = 2000
            end
            cha_dict[id] = feedings
            push!(charray, [id feedings calibrated num_feedings last_checkedin last_fed_diff] )
        end
    end

    charray = sort(charray, by = p -> p[6])
    return charray, cha_dict
end

function getcheckdict(charray)
    domain = Ft2k.conf["Watchman"]["domain"]
    checkdict =OrderedDict()
    for r in charray
        id = r[1]
        checkdict[id] = hbox(
                             pad(1em,node(:div, HTML("$id"))),
                             pad(1em,node(:div, string(r[4]))), #
        pad(1em,node(:div, string(r[5]))), #
        pad(1em, node(:div, string(r[6]))), #
        pad(1em, node(:div, r[3])),
        pad(1em, node(:div, r[2])) #
       )
    end
    return checkdict
end

function getfinalapp(tdict, fdict, ids, sqlite_db,apps, checkdict)
    apps = values(apps)
    ui = vbox( # put things one on top of the other
              node(:div, "Dashboard 3000"),
              vbox(values(checkdict)...)
             )
    base_app = page("/", req -> ui)
    base_app = Mux.mux(base_app)

    appfinal = Mux.stack(base_app, apps...)
    return appfinal
end

function get_dbvalues()
    time_db, feed_db, sqlite_db = pull_configure_data(Ft2k.conf)
    tdict, fdict = get_dicts(time_db, feed_db)
    ids= union(collect(keys(tdict)),collect(keys(fdict)), unique(sqlite_db[1][:,1]))
    return tdict, fdict, ids, sqlite_db
end

# goal here is to load labels from csv file
function load_labels()
    labelsbackup = joinpath(pwd(), "data", "labelsbackup.csv")
    # load csv file, create if not there already
    if !isfile(labelsbackup)
        writedlm(labelsbackup, ["feeder" "label"],',')
    else
        labels = readdlm(labelsbackup, ',')
    end

    # parse into dict
    labels = Dict(zip(labels[:,1], labels[:,2]))

    # return Dict
    return labels
end

# goal here is to backup labels to csv file, overwriting
function backup_labels(labels)
    labelsbackup = joinpath(pwd(), "data", "labelsbackup.csv")
    writedlm(labelsbackup, labels, ',')
end

#=
function set_label(id, label, labels)

println("Label for $id set to $label")
end
=#

function update_plots(fdicto, ido)
    id = ido[]
    p11 = Array{Any,1}()
    p12 = Array{Any,1}()
    p21 = Array{Any,1}()
    p22 = Array{Any,1}()
    p31 = Array{Any,1}()
    p32 = Array{Any,1}()
    if !haskey(fdicto[], id)
        println("No $id in fdict")
    end
    feed_day, feed_num= pullfeedingcountsbyday(fdicto[][ido[]])
    sortedfeed = sortslices(hcat(feed_day, feed_num), dims = 1)
    println(sortedfeed)

    if length(fdicto[][id])> 50
        p11 = collect(1:50)
        p12 = map(p-> p[3]-p[2], fdicto[][id][end-49:end])
        p21 = map(p -> p[6],fdicto[][id][end-49:end])
        p22 = map(p -> p[3]-p[2], fdicto[][id][end-49:end])
        p31 = sortedfeed[:,1] # hcat(sortedfeed[:,1], sortedfeed[:,2])
        p32 = sortedfeed[:,2] # hcat(sortedfeed[:,1], sortedfeed[:,2])
    elseif length(fdicto[][id])< 2
        p11 = "No feedings"
        p21 = "No feedings"
        p31 = "No feedings"
        p12 = "No feedings"
        p22 = "No feedings"
        p32 = "No feedings"
    else
        ll = length(fdicto[][id])
        p11 = collect(1:ll)
        p21 = map(p -> p[6],fdicto[][id])
        p31 = sortedfeed[:,1] #hcat(sortedfeed[:,1], sortedfeed[:,2])
        p12 = map(p-> p[3]-p[2], fdicto[][id])
        p22 =map(p -> p[3]-p[2], fdicto[][id])
        p32 = sortedfeed[:,2] #hcat(sortedfeed[:,1], sortedfeed[:,2])
    end
    # strictly debugging
    #map(p -> checkp(p), (p11, p12, p21, p22, p31, p32))
    return p11, p12, p21, p22, p31, p32
end

function reroute(p, domain, ido)
    println("Testing")
    println(p)
    ido[] = p
    # labelchooser[] = labelso[ido[]]
    #    labelmaker = textbox(hint= "Create New Label"; value = "")
    println(domain)
    node(:div, HTML(" <meta http-equiv=\"refresh\" content=\"1; URL=http://$domain:8500/\"/> "))
    #  node(:div,"
Hi!
")
    #   respond("
Hi!
")
    #,respond(" <meta http-equiv=\"refresh\" content=\"0; URL=http://localhost:8500/\"/> "
end

# this will not update, need an immutable app structure that allows new units to be added
# This is the main loop for restarting the static pages every 10 min
function ml()
    ENV["GKSwstype"]="svg"

    tdict, fdict, ids, sqlite_db = get_dbvalues()
    #host = Ft2k.conf["Watchman"]["domain"]
    ids = convert(Array{String}, ids)

    # TODO remove for debug
    #ids = ids[1:1]
    ido = Observable(ids[1])

    # this loads any previously saved labels/metadata
    labels = load_labels()

    # make values observables
    tdicto = Observable(tdict)
    fdicto = Observable(fdict)
    idso = Observable(ids)
    sqlite_dbo = Observable(sqlite_db)
    labelso = Observable(labels)

    # setup feeding scheme:
    feedingscheme = sqlite_dbo[][1][findfirst(p -> p == ido[], sqlite_dbo[][1][:,1]), 2]
    feedingscheme = Observable(feedingscheme)
    on(sqlite_dbo) do p
        feedingscheme[] = sqlite_dbo[][1][findfirst(p -> p == ido[], sqlite_dbo[][1][:,1]), 2]
    end
    feedschemetext = pad(1em, feedingscheme)

    # check boxes
    toggledfeeders = Observable(Dict(map(i -> i => true, ids)))
    toggledfeeder = checkboxes(Dict("Enrolled?" => true))
    on(p -> toggledfeeders[ido[]] = !toggledfeeders[ido[]], toggledfeeder)
    #
    # this loads any previously saved labels/metadata
    labels =load_labels()
    labelso = Observable(labels)
    labelslist = Observable(unique(values(labelso[])))# Observable(["", "LifespanCohort1"])

    labelmaker = textbox(value = "")
    submit_button = button("Make New Label")

    on(submit_button) do p
        push!(labelslist[], labelmaker[])
        # x = unique(values(labelso[]))
        # labelslist[] = x
        labelslist[] = labelslist[]
    end

    # setup labels
    if !haskey(labelso[], ido[])
        labelso[][ido[]] = "no label"
    end

    # dropdown input
    textlabel = Observable(labelso[][ido[]])
    labelchooser = dropdown(labelslist)
    on(labelchooser) do p
        println("labelchoooser chooses $p")
        labelso[][ido[]] = p
        textlabel[] = p
    end

    # need to get calibration value here TODO
    calibrated = Observable("no")
    calibratedtext = pad(2em, calibrated)

    # calibrate button
    calibrate_button = button("Calibrate")
    on(calibrate_button) do p
        calibrated[] = "yes"
        send_cal_msg(ido[])
    end

    # debug 7
    debug_box_7 = button("Set orders for 7 times")
    on(n -> send_deltaorder_msg(ido[], "7times"), debug_box_7)
    # debug 12
    debug_box_12 = button("Set orders for 12 times")
    on(n -> send_deltaorder_msg(ido[], "12times"), debug_box_12)
    # debug manual
    debug_box_manual = button("Set orders for manual")
    on(n -> send_deltaorder_msg(ido[], "manual"), debug_box_manual)
    # debug DR
    debug_box_DR = button("Set orders for DR")
    on(n -> send_deltaorder_msg(ido[], "DR"), debug_box_DR)
    # debug compressed 7
    debug_box_comp7 =button("Set orders for Compressed7")
    on(n -> send_deltaorder_msg(ido[], "comp7"), debug_box_comp7)

    # plots
    #    p11, p12, p21, p22, p31, p32 = Ft2k.setupplots(fdicto, ido[])
    p11 =""
    p12 = ""
    p21 = ""
    p22 =""
    p31 =""
    p32 =""
    p11, p12, p21, p22, p31, p32= update_plots(fdicto, ido) # , p31, p32

    println("plotting plots")
    plt1 = Observable(plot());
    plt2 = Observable(scatter()); # scatter
    plt3 = Observable(bar());

    #plt1 = Observable(plot(p11, p12, legend = false, ylim = (-10, 1024)));
    #plt2 = Observable(scatter(p21, p22, legend = false, ylim = (-10, 1024))); # scatter
    #plt1 = Interact.@map plot(&p11, &p12, legend = false, ylim=(-10,1024));
    #plt2 = Interact.@map scatter(&p21, &p22, legend = false, ylim = (-10, 1024)); # scatter
    #plt3 = Interact.@map plot(&p31, &p32, legend = false, ylim = (0, 15));

    on(ido) do p
        println("trying to update plots")
        # strictly debugging
        #map(p -> checkpo(p), (p11, p12, p21, p22, p31, p32))
        p1, p2, p3, p4, p5, p6 = update_plots(fdicto, ido) #, p5, p6
        plt1[] = plot(p1, p2, legend = false, ylim = (-10, 1024));
        plt2[] = scatter(p3, p4, legend = false, ylim = (-10, 1024));
        plt3[] = bar(p5, p6, ylim = (0, 16));
        plt3[] = vline!([p5[end]])
        feedingscheme[] = sqlite_dbo[][1][findfirst(p -> p == ido[], sqlite_dbo[][1][:,1]), 2]
    #    labelmaker = textbox(hint= "Create New Label"; value = "")
    end

    textout = pad(2em, textlabel)
    plt1out = pad(1em, plt1)
    plt2out = pad(1em, plt2)
    plt3out = pad(1em, plt3)
    labelmakerpad = pad(1em, labelmaker)
    submitpad = pad(1em, submit_button)
    labelchooserpad = pad(1em, labelchooser)

    # would be good to have colors on these debug boxes to distinguish them
    pageui = hbox(
                  vbox(pad(1em, ido),calibratedtext, feedschemetext,
                       hbox(toggledfeeder,
                            textout,
                            labelchooserpad, labelmakerpad, submitpad),
                       plt1out, plt2out, plt3out),#, plt3
                  vbox(pad(2em, calibrate_button),
                       pad(2em, debug_box_7),
                       pad(2em, debug_box_12),
                       pad(2em, debug_box_manual),
                       pad(2em, debug_box_DR),
                       pad(2em, debug_box_comp7),
                      ))

    rerouting_apps  = []
    for p in ids
        test=page("/$p", req->reroute(p, "www.ft2k.com", ido))#respond("
Hi!
")) #
        push!(rerouting_apps, test)
    end

    p2 = pad(1em, ido)

    base_app = page("/", req -> pageui)
    overview_app = page("/all", req -> p2)
    appfinal = Mux.stack(base_app, overview_app, rerouting_apps...)

    # this updates feeding values, but only every 60 seconds, need to make a function?
    logged = []
   #@async while true
   @async while true
        sleep(120)
        tdict, fdict, ids, sqlite_db = get_dbvalues()
        println("just updated fdict and co")
        ids = convert(Array{String}, ids)

        # update all observables somehow, this seems off
        fdicto[] = fdict
        tdicto[] = tdict
        idso[] = ids
        sqlite_dbo[] = sqlite_db

        # email check
        feedersDown = check_man_down(tdicto[])

        # check whether any feedings have been missed outside given window
        missedfeedings =check_missed_feedings(fdicto[], sqlite_dbo[])

        # Email out as needed
        logged = Ft2k.checkAndEmail(logged, feedersDown, missedfeedings)

        # backup labels
        println("backing up labels")
        println(labels)
        backup_labels(labels)
        GC.gc()
        println("restarting async loop")
    end

    host = "0.0.0.0"
    serving = WebIO.webio_serve(appfinal, host, 8500)
end

function debug_deltaorder(id, secondsfromnow, interval, numberoffeedings)
    domain = Ft2k.conf["Watchman"]["domain"]
    nowtime = Dates.Time(Dates.now() - Hour(7))
    # currently set for PST day light savings
    nowtime = Hour(nowtime).value * 60*60 + Minute(nowtime).value * 60 + Second(nowtime).value #- 7*60*60
    println(nowtime)

    feedingtimes = [t for t in (nowtime + 60):interval:(nowtime+120*numberoffeedings)]
    joinedtimes = join(feedingtimes, ",")
    try
        run(`mosquitto_pub --cafile /etc/mosquitto/ca_certificates/ca.crt -h $domain -p 8883 -t deltaorder -m $id:$joinedtimes`)
        println("sent mqtt msg to deltaorder for $id")
        println("feedings now: $joinedtimes")
    catch
        println("ERROR: could not update commands")
    end
end