text:
module OperationsTests

using Test
using Gridap.TensorValues
using Gridap.Arrays
using LinearAlgebra

# Comparison

a = VectorValue(1,2,3)
b = VectorValue(1,3,3)

@test (a < b) == true
@test (a <= b) == true
@test (a == b) == false
@test (a >= b) == false
@test (a > b) == false

a = VectorValue(1,2,3)
b = VectorValue(2,1,6)

@test a==a
@test a ≈ a
@test a!=b
@test [a,a] == [a,a]
@test [a,a] ≈ [a,a]

# Addition / subtraction

c = +a
r = VectorValue(1,2,3)
@test c == r

c = -a
r = VectorValue(-1,-2,-3)
@test c == r

c = a + b
r = VectorValue(3,3,9)
@test c == r

c = a - b
r = VectorValue(-1,1,-3)
@test c == r

a = VectorValue{0,Float64}()
b = VectorValue{0,Int}()

c = +a
r = VectorValue{0,Float64}()
@test c == r

c = -a
r = VectorValue{0,Float64}()
@test c == r

c = a + b
r = VectorValue{0,Float64}()
@test c == r

c = a - b
r = VectorValue{0,Float64}()
@test c == r

c = a*2
r = VectorValue{0,Float64}()
@test c == r

c = 2*a
r = VectorValue{0,Float64}()
@test c == r

a = TensorValue(1,2,3,4)
b = TensorValue(5,6,7,8)

c = +a
r = a
@test c==r

c = -a
r = TensorValue(-1,-2,-3,-4)
@test c==r

c = a - b
r = TensorValue(-4, -4, -4, -4)
@test c==r

a = SymTensorValue(1,2,3)
b = SymTensorValue(5,6,7)

c = -a
r = SymTensorValue(-1,-2,-3)
@test c==r

c = a + b
r = SymTensorValue(6,8,10)
@test c==r

# Matrix Division

a = VectorValue(1,2,3)

t = one(TensorValue{3,3,Int})
c = t\a
@test c == a

st = one(SymTensorValue{3,Int})
c = st\a
@test c == a

# Operations by a scalar

t = TensorValue(1,2,3,4,5,6,7,8,9)
st = SymTensorValue(1,2,3,5,6,9)
s4ot = one(SymFourthOrderTensorValue{2,Int})
a = VectorValue(1,2,3)

c = 2 * a
@test isa(c,VectorValue{3,Int})
r = VectorValue(2,4,6)
@test c == r

c = a * 2
@test isa(c,VectorValue{3,Int})
r = VectorValue(2,4,6)
@test c == r

c = 2 + a
@test isa(c,VectorValue{3,Int})
r = VectorValue(3,4,5)
@test c == r

c = a + 2
@test isa(c,VectorValue{3,Int})
@test c == r

c = a / 2
@test isa(c,VectorValue{3,Float64})
r = VectorValue(1/2,1.0,3/2)
@test c == r

c = 2 * t
@test isa(c,TensorValue{3})
r = TensorValue(2, 4, 6, 8, 10, 12, 14, 16, 18)
@test c == r

c = t * 2
@test isa(c,TensorValue{3})
r = TensorValue(2, 4, 6, 8, 10, 12, 14, 16, 18)
@test c == r

c = t + 2
@test isa(c,TensorValue{3,3,Int})
r = TensorValue(3, 4, 5, 6, 7, 8, 9, 10, 11)
@test c == r

c = 2 * st
@test isa(c,SymTensorValue{3})
r = SymTensorValue(2,4,6,10,12,18)
@test c == r

c = st * 2
@test isa(c,SymTensorValue{3})
r = SymTensorValue(2,4,6,10,12,18)
@test c == r

c = st + 2
@test isa(c,SymTensorValue{3})
r = SymTensorValue(3,4,5,7,8,11)
@test c == r

c = 2 * s4ot
@test isa(c,SymFourthOrderTensorValue{2})
r = SymFourthOrderTensorValue(2,0,0, 0,1,0, 0,0,2)
@test c == r

c = s4ot * 2
@test isa(c,SymFourthOrderTensorValue{2})
r = SymFourthOrderTensorValue(2,0,0, 0,1,0, 0,0,2)
@test c == r

c = c + 0
@test isa(c,SymFourthOrderTensorValue{2})
r = SymFourthOrderTensorValue(2,0,0, 0,1,0, 0,0,2)
@test c == r

# Dot product (simple contraction)

a = VectorValue(1,2,3)
b = VectorValue(2,1,6)

t = TensorValue(1,2,3,4,5,6,7,8,9)
s = TensorValue(9,8,3,4,5,6,7,2,1)
st = SymTensorValue(1,2,3,5,6,9)
st2 = SymTensorValue(9,6,5,3,2,1)

c = a ⋅ b
@test isa(c,Int)
@test c == 2+2+18

c = t ⋅ a
@test isa(c,VectorValue{3,Int})
r = VectorValue(30,36,42)
@test c == r

c = st ⋅ a
@test isa(c,VectorValue{3,Int})
r = VectorValue(14,30,42)
@test c == r

c = s ⋅ t
@test isa(c,TensorValue{3,3,Int})
r = TensorValue(38,24,18,98,69,48,158,114,78)
@test c == r

c = st ⋅ st2
@test isa(c,TensorValue{3,3,Int})
r = TensorValue(36, 78, 108, 18, 39, 54, 12, 26, 36)
@test c == r

c = a ⋅ st
@test isa(c,VectorValue{3,Int})
r = VectorValue(14,30,42)
@test c == r

# Inner product (full contraction)

c = 2 ⊙ 3
@test c == 6

c = a ⊙ b
@test isa(c,Int)
@test c == 2+2+18

c = inner(t,s)
@test isa(c,Int)
@test c == 185

c = inner(st,st2)
c = st ⊙ st2
@test isa(c,Int)
@test c == inner(TensorValue(get_array(st)),TensorValue(get_array(st2)))

# Reductions

a = VectorValue(1,2,3)

c = sum(a)
@test isa(c,Int)
@test c == 1+2+3

c = maximum(a)
@test isa(c,Int)
@test c == 3

c = minimum(a)
@test isa(c,Int)
@test c == 1

# Outer product (aka dyadic product)

a = VectorValue(1,2,3)
e = VectorValue(2,5)

c = outer(2,3)
c = 2 ⊗ 3
@test c == 6

r = VectorValue(2,4,6)
c = outer(2,a)
c = 2 ⊗ a
@test isa(c,VectorValue{3,Int})
@test c == r

c = outer(a,2)
@test isa(c,VectorValue{3,Int})
@test c == r

c = outer(a,e)
c = a ⊗ e
@test isa(c,TensorValue{3,2,Int})
r = TensorValue{3,2,Int}(2,4,6,5,10,15)
@test c == r

e = VectorValue(10,20)
k = TensorValue(1,2,3,4)
c = outer(e,k)
@test c == ThirdOrderTensorValue{2,2,2}(10, 20, 20, 40, 30, 60, 40, 80)

@test tr(c) == VectorValue(50,110)

# Cross product

a = VectorValue(1,2,3)
b = VectorValue(3,1,6)
c = VectorValue(9,3,-5)
@test a × b == c

a = VectorValue(2,3)
b = VectorValue(5,2)
@test a × b == -11

a = VectorValue(1.0,5.0,-4.0)
b = VectorValue(6.0,2.0,3.0)
c = VectorValue(23.0,-27.0,-28.0)
@test cross(a, b) == c

a = VectorValue(4.0,1.0)
b = VectorValue(3.0,-2.0)
@test cross(a, b) == -11.0

# Linear Algebra

t = TensorValue(10,2,30,4,5,6,70,8,9)

c = det(t)
@test c ≈ -8802.0
@test det(t) == det(TensorValue(get_array(t)))
@test inv(t) == inv(TensorValue(get_array(t)))

c = inv(t)
@test isa(c,TensorValue{3})

st = SymTensorValue(9,8,7,5,4,1)
@test det(st) == det(TensorValue(get_array(st)))
@test inv(st) == inv(TensorValue(get_array(st)))

t = TensorValue(10)
@test det(t) == 10
@test inv(t) == TensorValue(1/10)

t = TensorValue(1,4,-1,1)
@test det(t) == det(TensorValue(get_array(t)))
@test inv(t) == inv(TensorValue(get_array(t)))

# Measure

a = VectorValue(1,2,3)
c = meas(a)
@test c ≈ 3.7416573867739413

t = TensorValue(10,2,30,4,5,6,70,8,9)
c = meas(t)
@test c ≈ 8802.0

st = SymTensorValue(1,2,3,5,6,9)
@test meas(st) == meas(TensorValue(get_array(st)))

v = TensorValue{1,2}(10,20)
@test meas(v) == sqrt(500)

v = TensorValue{2,3}(1,0,0,1,0,0)
@test meas(v) ≈ 1.0

v = TensorValue{2,3}(1,0,0,1,1,0)
@test meas(v) ≈ sqrt(2)

# Broadcasted operations

a = VectorValue(1,2,3)
b = VectorValue(2,1,6)

A = Array{VectorValue{3,Int},2}(undef,(4,5))
A .= a
@test A == fill(a,(4,5))

C = .- A

r = VectorValue(-1,-2,-3)
R = fill(r,(4,5))
@test C == R

B = fill(b,(4,1))

C = A .+ B

r = VectorValue(3,3,9)
R = fill(r,(4,5))
@test C == R

C = A .- B

r = VectorValue(-1,1,-3)
R = fill(r,(4,5))
@test C == R

# conj

v = VectorValue(1,0)
@test v == v'
@test conj(v) == v'

t = TensorValue(1,2,3,4)
@test tr(t) == 5

t = TensorValue(1,2,3,4,5,6,7,8,9)
@test tr(t) == 15

st = SymTensorValue(1,2,3,5,6,9)
@test tr(st) == tr(TensorValue(get_array(st)))

@test get_array(symmetric_part(t)) == get_array(TensorValue(1.0, 3.0, 5.0, 3.0, 5.0, 7.0, 5.0, 7.0, 9.0))
@test symmetric_part(st) == symmetric_part(TensorValue(get_array(st)))

a = TensorValue(1,2,3,4)
b = a'
@test adjoint(a) == b
@test b == TensorValue(1,3,2,4)
@test a⋅b == TensorValue(10,14,14,20)

a = TensorValue(1,2,3,4)
b = a'
@test transpose(a) == b
@test b == TensorValue(1,3,2,4)
@test a⋅b == TensorValue(10,14,14,20)

sa = SymTensorValue(1,2,3,5,6,9)
sb = sa'
@test adjoint(sa) == sb
@test sb == SymTensorValue(1,2,3,5,6,9)
@test sa⋅sb == TensorValue(get_array(sa))⋅TensorValue(get_array(sb))

sa = SymTensorValue(1,2,3,5,6,9)
sb = sa'
@test transpose(sa) == sb
@test sb == SymTensorValue(1,2,3,5,6,9)
@test sa⋅sb == TensorValue(get_array(sa))⋅TensorValue(get_array(sb))

u = VectorValue(1.0,2.0)
v = VectorValue(2.0,3.0)
@test dot(u,v) ≈ inner(u,v)
@test norm(u) ≈ sqrt(inner(u,u))

a = VectorValue(1.0,2.0)
b = VectorValue(2.0,3.0)
@test [a,b] ≈ [a,b]

a = VectorValue{0,Int}()
@test a ≈ a

λ = 1
μ = 1
ε = SymTensorValue(1,2,3)
σ = λ*tr(ε)*one(ε) + 2*μ*ε
@test isa(σ,SymTensorValue)
@test (σ ⊙ ε) == 52
#@test σ:ε == 52

I = one(SymFourthOrderTensorValue{2,Int})
@test I[1,1,1,1] == 1
@test I[1,2,1,2] == 0.5
@test I[2,1,1,2] == 0.5
@test I[2,2,2,2] == 1

@test I ⊙ ε == ε
#@test I : ε == ε

a = TensorValue(1,2,3,4)
b = I ⊙ a
@test b == symmetric_part(a)
#b = I : a
#@test b == symmetric_part(a)

σ1 = λ*tr(ε)*one(ε) + 2*μ*ε
C = 2*μ*one(ε⊗ε) + λ*one(ε)⊗one(ε)
σ2 = C ⊙ ε
@test σ1 == σ2
#σ2 = C : ε
#@test σ1 == σ2

end # module OperationsTests


text:
using juliaset
using Base.Test

# non-divergence check;
@test juliaset.createjuliaset(ones(1),zeros(1),10,100,0) == 11*ones(1,1)

ext: jl

max_stars_repo_name: asgoldin/juliaset

max_issues_repo_name: asgoldin/juliaset

max_forks_repo_name: asgoldin/juliaset


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("logdval.jl")
include("abstract_density.jl")
include("shaped_densities.jl")
include("generic_density.jl")
include("const_density.jl")
include("density_product.jl")
include("distribution_density.jl")
include("posterior_density.jl")
include("renormalize_density.jl")
include("truncate_density.jl")
include("transformed_density.jl")
include("external_density.jl")
include("testDensities.jl")


text:
################################################################################
# Davis-Yin splitting iterator
#
# See:
# [1] Davis, Yin "A Three-Operator Splitting Scheme and its Optimization Applications",
# Set-Valued and Variational Analysis, vol. 25, no. 4, pp 829–858 (2017).
#

struct DavisYinIterator{I <: Integer, R <: Real, T1 <: BlockArray{R}, T2 <: BlockArray{R}} <: ProximalAlgorithm{I, Tuple{T1, T2}}
    z::T1
    f
    g
    h
    L
    gamma::R
    lambda::R
    maxit::I
    tol::R
    verbose::I
    verbose_freq::I
    x::T1
    y::T2
    z_half::T1
    w::T1
    FPR_x::T1
end

################################################################################
# Constructor(s)

function DavisYinIterator(z0::BlockArray{R}; f=Zero(), g=Zero(), h=Zero(), L=Identity(blocksize(z0)), gamma::R=1.0, maxit::I=10000, tol::R=1e-4, verbose=1, verbose_freq = 100) where {I, R}
    # TODO
end

################################################################################
# Utility methods

maxit(sol::DavisYinIterator) = sol.maxit

converged(sol::DavisYinIterator, it) = it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol

verbose(sol::DavisYinIterator) = sol.verbose > 0
verbose(sol::DavisYinIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0))

function display(sol::DavisYinIterator)
    # TODO
end

function display(sol::DavisYinIterator, it)
    # TODO
end

function Base.show(io::IO, sol::DavisYinIterator)
    # TODO
end

################################################################################
# Initialization

function initialize!(sol::DavisYinIterator)
    return
end

################################################################################
# Iteration

function iterate!(sol::DavisYinIterator{I, T}, it::I) where {I, T}
    # x = prox(g, z)
    # y = Lx
    # w = ∇h(y)
    # z_half = 2x - z - γL'w
    # z = z + λ(prox(f, z_half) - x)
end

################################################################################
# Solver interface

function DavisYin(x0; kwargs...)
    sol = DavisYinIterator(x0; kwargs...)
    it, point = run!(sol)
    return (it, point, sol)
end


text:
module OutputWriters

export
    write_output,
    JLD2OutputWriter, FieldOutput, FieldOutputs,
    NetCDFOutputWriter, write_grid_and_attributes,
    Checkpointer, restore_from_checkpoint

using CUDA

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Fields
using Oceananigans.Architectures

using Oceananigans: AbstractOutputWriter, @hascuda
using Oceananigans.Fields: OffsetArray

Base.open(ow::AbstractOutputWriter) = nothing
Base.close(ow::AbstractOutputWriter) = nothing

include("output_writer_utils.jl")
include("jld2_output_writer.jl")
include("netcdf_output_writer.jl")
include("checkpointer.jl")

end


text:
# fisher.jl
# Fisher's exact test
#
# Copyright (C) 2013   Simon Kornblith
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export FisherExactTest

immutable FisherExactTest <: HypothesisTest
    # Format:
    # X1  X2
    # Y1  a  b
    # Y2  c  d
    a::Int
    b::Int
    c::Int
    d::Int

    # conditional maximum likehood estimate of odd ratio
    ω::Float64

    function FisherExactTest(a::Int, b::Int, c::Int, d::Int)
        ω = cond_mle_odds_ratio(a, b, c, d)
        new(a, b, c, d, ω)
    end
end

testname(::FisherExactTest) = "Fisher's exact test"
population_param_of_interest(x::FisherExactTest) = ("Odds ratio", 1.0, x.ω) # parameter of interest: name, value under h0, point estimate

# The sizing argument to print_matrix was removed during the 0.5 dev period
if VERSION < v"0.5.0-dev+1936"
    function _print_matrix(io::IO, X::AbstractVecOrMat, pre::AbstractString)
        m = typemax(Int)
        Base.print_matrix(io, X, (m, m), pre)
    end
else
    function _print_matrix(io::IO, X::AbstractVecOrMat, pre::AbstractString)
        Base.print_matrix(io, X, pre)
    end
end

function show_params(io::IO, x::FisherExactTest, ident="")
    println(io, ident, "contingency table:")
    _print_matrix(io, [x.a x.b; x.c x.d], repeat(ident, 2))
    println(io)
end

# DOC: for tail=:both there exist multiple ``method``s for computing a pvalue and the corresponding ci.
function pvalue(x::FisherExactTest; tail=:both, method=:central)
    if tail == :both && method != :central
        if method == :minlike
            p = pvalue_both_minlike(x)
        else
            throw(ArgumentError("method=$(method) is not implemented yet"))
        end
    else
        p = pvalue(Hypergeometric(x.a + x.b, x.c + x.d, x.a + x.c), x.a, tail=tail)
    end
    p = max(min(p, 1.0), 0.0)

    return p
end

function pvalue_both_minlike(x::FisherExactTest, ω::Float64=1.0)
    a, b, c, d = reorder(x.a, x.b, x.c, x.d)

    dist = FisherNoncentralHypergeometric(a+b, c+d, a+c, ω)

    p = pdf(dist, a)
    v = nextfloat(p)
    if a != 0
        p += cdf(dist, a-1)
    end

    # Add p-values of all tables in other tail equally or less probable
    for i = a+c:-1:a+1
        curp = pdf(dist, i)
        if curp > v
            break
        end
        p += curp
    end
    p
end

# confidence interval by inversion of p-value
function ci(x::FisherExactTest, alpha::Float64=0.05; tail=:both, method=:central)
    check_alpha(alpha)
    dist(ω) = FisherNoncentralHypergeometric(x.a+x.b, x.c+x.d, x.a+x.c, ω)
    obj(ω) = pvalue(dist(ω), x.a, tail=tail) - alpha

    if tail == :left # upper bound
        if (x.a == maximum(dist(1.0)))
            (0.0, Inf)
        else
            (0.0, fzero(obj, find_brackets(obj)...))
        end
    elseif tail == :right # lower bound
        if (x.a == minimum(dist(1.0)))
            (0.0, Inf)
        else
            (fzero(obj, find_brackets(obj)...), Inf)
        end
    elseif tail == :both
        if method == :central
            (ci(x, alpha/2; tail=:right)[1], ci(x, alpha/2; tail=:left)[2])
        else
            throw(ArgumentError("method=$(method) is not implemented yet"))
        end
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

## helpers

function reorder(a,b,c,d)
    if a + c > b + d
        a, b, c, d = b, a, d, c
    end
    if a/c > b/d
        a, b, c, d = c, d, a, b
    end
    (a, b, c, d)
end

# find values x_lower and x_upper, s.t. f(x_lower) < 0 and f(x_upper) > 0 or vice versa
function find_brackets(f::Function, x_init::Float64=1.0)
    f_init = f(x_init)

    if f_init > f(x_init + 1.0)
        find_brackets(x -> -f(x), x_init)
    else
        x_upper = x_lower = x_init
        if f_init > 0.0
            while x_lower > eps(0.0) && f(x_lower) > 0.0
                x_lower /= 2
            end
        else
            while f(x_upper) < 0.0
                x_upper *= 2
            end
        end
        (x_lower, x_upper)
    end
end

# find odds ratio ω that maximizes the (conditional) likelihood;
# since the mode and mean of Fisher's Noncentral Hypergeometric distribution
# coincide, this is equivalent to find ω, s.t., mean(dist(ω)) = a
function cond_mle_odds_ratio(a::Int, b::Int, c::Int, d::Int)
    dist(ω) = FisherNoncentralHypergeometric(a+b, c+d, a+c, ω)

    if (a == minimum(dist(1.0)))
        0.0
    elseif a == maximum(dist(1.0))
        Inf
    else
        obj(ω) = mean(dist(ω))-a
        fzero(obj, find_brackets(obj)...)
    end
end


text:
using Kinetic

# initialization
set, ctr, xface, yface, t = initialize("config.toml")

# solution algorithm
t = solve!(set, ctr, xface, yface, t)

# visualization
plot_contour(set, ctr)


text:
# generate examples
import Literate

# TODO: Remove items from `SKIPFILE` as soon as they run on the latest
# stable `Optim` (or other dependency)
#ONLYSTATIC = ["ipnewton_basics.jl",]
ONLYSTATIC = []

EXAMPLEDIR = joinpath(@__DIR__, "src", "examples")
GENERATEDDIR = joinpath(@__DIR__, "src", "examples", "generated")
for example in filter!(x -> endswith(x, ".jl"), readdir(EXAMPLEDIR))
    input = abspath(joinpath(EXAMPLEDIR, example))
    script = Literate.script(input, GENERATEDDIR)
    code = strip(read(script, String))
    mdpost(str) = replace(str, "@__CODE__" => code)
    Literate.markdown(input, GENERATEDDIR, postprocess = mdpost,
                      documenter = !(example in ONLYSTATIC))
    Literate.notebook(input, GENERATEDDIR, execute = !(example in ONLYSTATIC))
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# run simulation

"""
	function simulate!(sim::Simulation;
		time::Float = Inf, duration::Float = Inf, numEvents::Int = -1,
		doPrint::Bool = false, printingInterval::Float = 1.0)
Run the simulation to completion, or to any specified stopping point (according to `time`, `duration`, or `numEvents`), whichever comes first. Returns `true` if simulation is complete; `false` otherwise.

# Keyword arguments
- `time` is the maximum time to simulate to; will stop if time of next event is greater
- `duration` is the maximum duration to simulate for; equivalent to setting `time = sim.time + duration`
- `numEvents` is the maximum number of events to simulate
- `doPrint` controls whether any lines are printed while simulating
- `printingInterval` is the time interval (seconds) between printing while simulating, if `doPrint = true`
"""
function simulate!(sim::Simulation;
	time::Real = Inf, duration::Real = Inf, numEvents::Real = Inf,
	doPrint::Bool = false, printingInterval::Real = 1.0)

	@assert(time == Inf || duration == Inf, "can only set one of: time, duration")
	@assert(time >= sim.time)
	@assert(duration >= 0)
	@assert(numEvents >= 0)
	@assert(printingInterval >= 1e-2, "printing too often can slow the simulation")

	duration != Inf && (time = sim.time + duration) # set end time based on duration

	# for printing progress
	startTime = Base.time()
	doPrint || (printingInterval = Inf)
	nextPrintTime = startTime + printingInterval
	eventCount = 0
	printProgress() = doPrint && print(@sprintf("\rsim time: %-9.2f sim duration: %-9.2f events simulated: %-9d real duration: %.2f seconds", sim.time, sim.time - sim.startTime, eventCount, Base.time() - startTime))

	# simulate
	stats = sim.stats # shorthand
	while !sim.complete && sim.eventList[end].time <= time && eventCount < numEvents
		if stats.doCapture && stats.nextCaptureTime <= sim.eventList[end].time
			captureSimStats!(sim, stats.nextCaptureTime)
		end

		simulateNextEvent!(sim)
		eventCount += 1

		if doPrint && Base.time() >= nextPrintTime
			printProgress()
			nextPrintTime += printingInterval
		end
	end

	printProgress()
	doPrint && println()

	if stats.doCapture && sim.complete
		captureSimStats!(sim, sim.endTime)
		populateSimStats!(sim)
	end

	return sim.complete
end

# simulate up to last event with time <= given time
# similar to calling simulate!(sim, time = time), but does not record statistics
# kept this for animation, to avoid problems with animating as an asynchronous task
function simulateToTime!(sim::Simulation, time::Float)
	while !sim.complete && sim.eventList[end].time <= time # eventList must remain sorted by non-increasing time
		simulateNextEvent!(sim)
	end
end

# run simulation to end
# similar to calling simulate!(sim), but does not record statistics
function simulateToEnd!(sim::Simulation)
	while !sim.complete
		simulateNextEvent!(sim)
	end
end

# set sim.backup to copy of sim (before running)
# note that for reducing memory usage, sim.backup does not contain backups of all fields
function backup!(sim::Simulation)
	@assert(!sim.used)
	if !isdefined(sim, :backup) sim.backup = Simulation() end
	fnamesDontCopy = (:backup, :net, :travel, :grid, :resim, :calls, :demand, :demandCoverage, :reps)
	for fname in setdiff(fieldnames(Simulation), fnamesDontCopy)
		setfield!(sim.backup, fname, deepcopy(getfield(sim, fname)))
	end
end
backupSim!(sim::Simulation) = backup!(sim) # compat

# reset sim from sim.backup
function reset!(sim::Simulation)
	@assert(isdefined(sim, :backup))
	@assert(!sim.backup.used)

	if sim.used
		fnamesDontCopy = (:backup, :net, :travel, :grid, :resim, :calls, :demand, :demandCoverage, :reps, :animating)
		# fnamesDontCopy lists fields that are not backed up in sim.backup, or will be reset by another method
		for fname in setdiff(fieldnames(Simulation), fnamesDontCopy)
			setfield!(sim, fname, deepcopy(getfield(sim.backup, fname)))
		end

		reset!(sim.calls)

		# reset resimulation state
		sim.resim.prevEventIndex = 0

		# reset travel and demand state
		sim.travel.recentSetsStartTimesIndex = 1
		sim.demand.recentSetsStartTimesIndex = 1
	end

	return nothing
end
resetSim!(sim::Simulation) = reset!(sim) # compat

# simulate next event in list
function simulateNextEvent!(sim::Simulation)
	# get next event, update event index and sim time
	event = getNextEvent!(sim.eventList)
	if event.form == nullEvent
		error()
	end
	sim.used = true
	sim.eventIndex += 1
	event.index = sim.eventIndex
	sim.time = event.time

	if sim.resim.use
		resimCheckCurrentEvent!(sim, event)
	elseif sim.writeOutput
		writeEventToFile!(sim, event)
	end

	simulateEvent!(sim, event)

	if length(sim.eventList) == 0
		# simulation complete
		@assert(sim.endTime == nullTime)
		@assert(sim.complete == false)
		sim.endTime = sim.time
		sim.complete = true
		for amb in sim.ambulances
			setAmbStatus!(sim, amb, amb.status, sim.time) # to account for duration spent with final status
		end
	end
end

# simulate event
function simulateEvent!(sim::Simulation, event::Event)
	@assert(sim.time == event.time)

	# Find event simulation function to match event form.
	# Use of if/else here is faster than dict lookup, and using multiple dispatch was very slow (at least doubled simulation run time).
	form = event.form # shorthand
	if form == ambGoesToSleep simulateEventAmbGoesToSleep!(sim, event)
	elseif form == ambWakesUp simulateEventAmbWakesUp!(sim, event)
	elseif form == callArrives simulateEventCallArrives!(sim, event)
	elseif form == considerDispatch simulateEventConsiderDispatch!(sim, event)
	elseif form == ambDispatched simulateEventAmbDispatched!(sim, event)
	elseif form == ambMobilised simulateEventAmbMobilised!(sim, event)
	elseif form == ambReachesCall simulateEventAmbReachesCall!(sim, event)
	elseif form == ambGoesToHospital simulateEventAmbGoesToHospital!(sim, event)
	elseif form == ambReachesHospital simulateEventAmbReachesHospital!(sim, event)
	elseif form == ambBecomesFree simulateEventAmbBecomesFree!(sim, event)
	elseif form == ambReturnsToStation simulateEventAmbReturnsToStation!(sim, event)
	elseif form == ambReachesStation simulateEventAmbReachesStation!(sim, event)
	elseif form == considerMoveUp simulateEventConsiderMoveUp!(sim, event)
	elseif form == ambMoveUpToStation simulateEventAmbMoveUpToStation!(sim, event)
	else error("Unknown event: ", form, ".")
	end
	# Format of event simulation functions:
	# - event may change relevant ambulance / call fields at event.time
	# - collect statistics
	# - event may then trigger future events / cancel scheduled events
end

function simulateEventAmbGoesToSleep!(sim::Simulation, event::Event)
	@assert(event.form == ambGoesToSleep)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambIdleAtStation) # ambulance should be at station before sleeping
	@assert(event.callIndex == nullIndex)

	setAmbStatus!(sim, ambulance, ambSleeping, sim.time)

	station = sim.stations[ambulance.stationIndex]
	updateStationStats!(station; numIdleAmbsChange = -1, time = sim.time)

	addEvent!(sim.eventList; parentEvent = event, form = ambWakesUp, time = sim.time + sleepDuration, ambulance = ambulance, station = station)
end

function simulateEventAmbWakesUp!(sim::Simulation, event::Event)
	@assert(event.form == ambWakesUp)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambSleeping) # ambulance should have been sleeping
	@assert(event.callIndex == nullIndex)

	setAmbStatus!(sim, ambulance, ambIdleAtStation, sim.time)
	ambulance.event = Event() # no event currently

	updateStationStats!(sim.stations[ambulance.stationIndex]; numIdleAmbsChange = +1, time = sim.time)
end

function simulateEventCallArrives!(sim::Simulation, event::Event)
	@assert(event.form == callArrives)
	@assert(event.ambIndex == nullIndex) # no ambulance should be assigned yet
	call = sim.calls[event.callIndex]
	@assert(call.status == callNullStatus)

	push!(sim.currentCalls, sim.calls[event.callIndex])
	setCallStatus!(call, callScreening, sim.time)

	addEvent!(sim.eventList; parentEvent = event, form = considerDispatch, time = sim.time + call.dispatchDelay, call = call)

	# add next call arrival to event queue
	if call.index < sim.numCalls
		addEvent!(sim.eventList, sim.calls[call.index + 1])
	end
end

function simulateEventConsiderDispatch!(sim::Simulation, event::Event)
	@assert(event.form == considerDispatch)
	@assert(event.ambIndex == nullIndex) # no ambulance should be assigned yet
	call = sim.calls[event.callIndex]
	@assert(in(call.status, (callScreening, callWaitingForAmb))) # callWaitingForAmb if call bumped

	# find ambulance to respond
	if sim.resim.use && sim.resim.doDispatch
		ambIndex = resimFindAmbToDispatch(sim, call)
	else
		ambIndex = sim.findAmbToDispatch!(sim, call)
	end
	# if no ambulance found, add call to queue,
	# will respond when an ambulance becomes free
	if ambIndex == nullIndex
		sim.addCallToQueue!(sim.queuedCallList, call)
		setCallStatus!(call, callQueued, sim.time)
	else
		ambulance = sim.ambulances[ambIndex]

		# remove any previously scheduled events
		if isGoingToStation(ambulance.status)
			deleteEvent!(sim.eventList, ambulance.event)
		elseif in(ambulance.status, (ambMobilising, ambGoingToCall))
			# if ambulance was going to call, redirect ambulance
			deleteEvent!(sim.eventList, ambulance.event)

			# bumped call has no ambulance assigned yet
			bumpedCall = sim.calls[ambulance.callIndex]
			bumpedCall.ambIndex = nullIndex
			bumpedCall.numBumps += 1 # stats

			# consider new dispatch for bumped call
			addEvent!(sim.eventList; parentEvent = event, form = considerDispatch, time = sim.time, call = bumpedCall)
		end

		# dispatch chosen ambulance
		addEvent!(sim.eventList; parentEvent = event, form = ambDispatched, time = sim.time, ambulance = ambulance, call = call)
	end
end

function simulateEventAmbDispatched!(sim::Simulation, event::Event)
	@assert(event.form == ambDispatched)
	ambulance = sim.ambulances[event.ambIndex]
	status = ambulance.status # shorthand
	@assert(isFree(status) || in(status, (ambMobilising, ambGoingToCall)))
	call = sim.calls[event.callIndex]
	@assert(in(call.status, (callScreening, callQueued, callWaitingForAmb))) # callWaitingForAmb if call bumped

	ambulance.callIndex = call.index
	ambulance.dispatchTime = sim.time
	mobiliseAmbulance = sim.mobilisationDelay.use && in(ambulance.status, (ambIdleAtStation, ambMobilising))
	if mobiliseAmbulance
		setAmbStatus!(sim, ambulance, ambMobilising, sim.time)
	else
		setAmbStatus!(sim, ambulance, ambGoingToCall, sim.time)
		changeRoute!(sim, ambulance.route, sim.responseTravelPriorities[call.priority], sim.time, call.location, call.nearestNodeIndex)
	end

	setCallStatus!(call, callWaitingForAmb, sim.time)
	call.ambIndex = event.ambIndex
	ambLoc = getRouteCurrentLocation!(sim.net, ambulance.route, sim.time)
	copy!(call.ambDispatchLoc, ambLoc)
	call.ambStatusBeforeDispatch = ambulance.prevStatus

	# stats
	if ambulance.prevStatus == ambIdleAtStation
		updateStationStats!(sim.stations[ambulance.stationIndex]; numIdleAmbsChange = -1, time = sim.time)
	end
	if sim.stats.recordDispatchStartLocCounts
		ambulance.dispatchStartLocCounts[ambLoc] = get(ambulance.dispatchStartLocCounts, ambLoc, 0) + 1
	end

	if mobiliseAmbulance
		mobilisationDelay = 0.0 # init
		if ambulance.prevStatus == ambMobilising
			# ambulance.mobilisation time unchanged
			@assert(ambulance.mobilisationTime != nullTime)
		else
			@assert(ambulance.prevStatus == ambIdleAtStation)
			ambulance.mobilisationTime = sim.time + rand(sim.mobilisationDelay.distrRng)
		end
		@assert(ambulance.mobilisationTime >= sim.time)
		addEvent!(sim.eventList; parentEvent = event, form = ambMobilised, time = ambulance.mobilisationTime, ambulance = ambulance, call = call)
	else
		addEvent!(sim.eventList; parentEvent = event, form = ambReachesCall, time = ambulance.route.endTime, ambulance = ambulance, call = call)
	end

	if sim.moveUpData.useMoveUp
		if in(sim.moveUpData.moveUpModule, (compTableModule, ddsmModule, multiCompTableModule, zhangIpModule, temp0Module, temp1Module, temp2Module))
			addEvent!(sim.eventList; parentEvent = event, form = considerMoveUp, time = sim.time, ambulance = ambulance, addEventToAmb = false)
		end
	end
end

function simulateEventAmbMobilised!(sim::Simulation, event::Event)
	@assert(event.form == ambMobilised)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambMobilising)
	call = sim.calls[event.callIndex]
	@assert(call.status == callWaitingForAmb)
	@assert(sim.mobilisationDelay.use)

	setAmbStatus!(sim, ambulance, ambGoingToCall, sim.time)
	changeRoute!(sim, ambulance.route, sim.responseTravelPriorities[call.priority], sim.time, call.location, call.nearestNodeIndex)
	ambulance.mobilisationTime = nullTime # reset

	# call status unchanged

	addEvent!(sim.eventList; parentEvent = event, form = ambReachesCall, time = ambulance.route.endTime, ambulance = ambulance, call = call)
end

function simulateEventAmbReachesCall!(sim::Simulation, event::Event)
	@assert(event.form == ambReachesCall)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambGoingToCall)
	call = sim.calls[event.callIndex]
	@assert(call.status == callWaitingForAmb)

	setAmbStatus!(sim, ambulance, ambAtCall, sim.time)

	setCallStatus!(call, callOnSceneTreatment, sim.time)

	# transport call to hospital if needed, otherwise amb becomes free
	if call.transport
		# transport to hospital
		addEvent!(sim.eventList; parentEvent = event, form = ambGoesToHospital, time = sim.time + call.onSceneDuration, ambulance = ambulance, call = call)
	else
		# amb becomes free
		addEvent!(sim.eventList; parentEvent = event, form = ambBecomesFree, time = sim.time + call.onSceneDuration, ambulance = ambulance, call = call)
	end
end

function simulateEventAmbGoesToHospital!(sim::Simulation, event::Event)
	@assert(event.form == ambGoesToHospital)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambAtCall)
	call = sim.calls[event.callIndex]
	@assert(call.status == callOnSceneTreatment)
	@assert(call.transport)

	# if hospital not specified for call, find closest hospital
	hospitalIndex = call.hospitalIndex
	if hospitalIndex == nullIndex
		hospitalIndex = nearestHospitalToCall!(sim, call, lowPriority)
	end
	@assert(hospitalIndex != nullIndex)
	hospital = sim.hospitals[hospitalIndex]

	setAmbStatus!(sim, ambulance, ambGoingToHospital, sim.time)
	changeRoute!(sim, ambulance.route, lowPriority, sim.time, hospital.location, hospital.nearestNodeIndex)

	setCallStatus!(call, callGoingToHospital, sim.time)
	call.chosenHospitalIndex = hospitalIndex

	addEvent!(sim.eventList; parentEvent = event, form = ambReachesHospital, time = ambulance.route.endTime, ambulance = ambulance, call = call)
end

function simulateEventAmbReachesHospital!(sim::Simulation, event::Event)
	@assert(event.form == ambReachesHospital)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambGoingToHospital)
	call = sim.calls[event.callIndex]
	@assert(call.status == callGoingToHospital)
	@assert(call.transport)

	setAmbStatus!(sim, ambulance, ambAtHospital, sim.time)

	setCallStatus!(call, callAtHospital, sim.time)
	sim.hospitals[call.chosenHospitalIndex].numCalls += 1 # stats

	addEvent!(sim.eventList; parentEvent = event, form = ambBecomesFree, time = sim.time + call.handoverDuration, ambulance = ambulance, call = call)
end

function simulateEventAmbBecomesFree!(sim::Simulation, event::Event)
	@assert(event.form == ambBecomesFree)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(in(ambulance.status, (ambAtCall, ambAtHospital)))
	call = sim.calls[event.callIndex]
	@assert(in(call.status, (callOnSceneTreatment, callAtHospital)))

	# remove call, processing is finished
	delete!(sim.currentCalls, call)
	setCallStatus!(call, callProcessed, sim.time)

	setAmbStatus!(sim, ambulance, ambFreeAfterCall, sim.time)
	ambulance.callIndex = nullIndex

	# if queued call exists, respond
	# otherwise return to station
	if length(sim.queuedCallList) > 0
		call = getNextCall!(sim.queuedCallList)
		@assert(call != nothing)

		# dispatch ambulance
		addEvent!(sim.eventList; parentEvent = event, form = ambDispatched, time = sim.time, ambulance = ambulance, call = call)
	else
		# return to station (can be cancelled by move up)
		addEvent!(sim.eventList; parentEvent = event, form = ambReturnsToStation, time = sim.time, ambulance = ambulance, station = sim.stations[ambulance.stationIndex])

		# consider move up
		# note that this event is created after above ambReturnsToStation event, so that it will happen first and can delete the ambReturnsToStation event if needed
		if sim.moveUpData.useMoveUp
			if in(sim.moveUpData.moveUpModule, (compTableModule, dmexclpModule, multiCompTableModule, priorityListModule, zhangIpModule, temp0Module, temp1Module, temp2Module))
				addEvent!(sim.eventList; parentEvent = event, form = considerMoveUp, time = sim.time, ambulance = ambulance, addEventToAmb = false)
			end
		end
	end
end

function simulateEventAmbReturnsToStation!(sim::Simulation, event::Event)
	@assert(event.form == ambReturnsToStation)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(ambulance.status == ambFreeAfterCall)
	@assert(event.callIndex == nullIndex)

	setAmbStatus!(sim, ambulance, ambReturningToStation, sim.time)
	station = sim.stations[ambulance.stationIndex]
	changeRoute!(sim, ambulance.route, lowPriority, sim.time, station.location, station.nearestNodeIndex)

	addEvent!(sim.eventList; parentEvent = event, form = ambReachesStation, time = ambulance.route.endTime, ambulance = ambulance, station = station)
end

function simulateEventAmbReachesStation!(sim::Simulation, event::Event)
	@assert(event.form == ambReachesStation)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(isGoingToStation(ambulance.status))
	@assert(event.callIndex == nullIndex)

	setAmbStatus!(sim, ambulance, ambIdleAtStation, sim.time)
	ambulance.event = Event() # no event currently

	updateStationStats!(sim.stations[ambulance.stationIndex]; numIdleAmbsChange = +1, time = sim.time)
end

function simulateEventConsiderMoveUp!(sim::Simulation, event::Event)
	@assert(event.form == considerMoveUp)
	ambulance = sim.ambulances[event.ambIndex] # ambulance that triggered consideration of move up
	@assert(event.callIndex == nullIndex)
	mud = sim.moveUpData # shorthand
	@assert(mud.useMoveUp)

	# call move up function
	movableAmbs = [] # init
	ambStations = [] # init

	if sim.resim.use && sim.resim.doMoveUp
		(movableAmbs, ambStations) = resimMoveUp(sim)
	else
		if mud.moveUpModule == compTableModule
			(movableAmbs, ambStations) = compTableMoveUp(sim)
		elseif mud.moveUpModule == ddsmModule
			(movableAmbs, ambStations) = ddsmMoveUp(sim)
		elseif mud.moveUpModule == dmexclpModule
			(movableAmbs, ambStations) = dmexclpMoveUp(sim, ambulance)
		elseif mud.moveUpModule == multiCompTableModule
			(movableAmbs, ambStations) = multiCompTableMoveUp(sim)
		elseif mud.moveUpModule == priorityListModule
			(movableAmbs, ambStations) = priorityListMoveUp(sim, ambulance)
		elseif mud.moveUpModule == zhangIpModule
			(movableAmbs, ambStations) = zhangIpMoveUp(sim)
		elseif mud.moveUpModule == temp0Module
			(movableAmbs, ambStations) = temp0MoveUp(sim)
		elseif mud.moveUpModule == temp1Module
			(movableAmbs, ambStations) = temp1MoveUp(sim)
		elseif mud.moveUpModule == temp2Module
			(movableAmbs, ambStations) = temp2MoveUp(sim)
		else
			error()
		end
	end

	for i = 1:length(movableAmbs)
		ambulance = movableAmbs[i]
		station = ambStations[i]

		# move up ambulance if ambulance station has changed
		if ambulance.stationIndex != station.index

			if isGoingToStation(ambulance.status) # ambulance.event.form == ambReachesStation
				# delete station arrival event for this ambulance
				deleteEvent!(sim.eventList, ambulance.event)
			elseif ambulance.status == ambFreeAfterCall && ambulance.event.form == ambReturnsToStation
				deleteEvent!(sim.eventList, ambulance.event)
			end

			addEvent!(sim.eventList; parentEvent = event, form = ambMoveUpToStation, time = sim.time, ambulance = ambulance, station = station)
		end
	end
end

function simulateEventAmbMoveUpToStation!(sim::Simulation, event::Event)
	@assert(event.form == ambMoveUpToStation)
	ambulance = sim.ambulances[event.ambIndex]
	@assert(event.callIndex == nullIndex)
	station = sim.stations[event.stationIndex] # station to move up to
	@assert(ambulance.stationIndex != station.index) # otherwise it is not move up

	setAmbStatus!(sim, ambulance, ambMovingUpToStation, sim.time)
	changeRoute!(sim, ambulance.route, lowPriority, sim.time, station.location, station.nearestNodeIndex)
	prevStationIndex = ambulance.stationIndex
	ambulance.stationIndex = station.index

	# stats
	if ambulance.moveUpFromStationIndex == station.index
		@assert(ambulance.prevStatus == ambMovingUpToStation)
		ambulance.numMoveUpsReturnToPrevStation += 1
	end
	if ambulance.prevStatus == ambIdleAtStation
		updateStationStats!(sim.stations[prevStationIndex]; numIdleAmbsChange = -1, time = sim.time)
		ambulance.moveUpFromStationIndex = prevStationIndex
	end
	if sim.stats.recordMoveUpStartLocCounts
		startLoc = getRouteCurrentLocation!(sim.net, ambulance.route, sim.time) # startLoc should be same as ambulance.route.startLoc
		station.moveUpStartLocCounts[startLoc] = get(station.moveUpStartLocCounts, startLoc, 0) + 1
	end

	addEvent!(sim.eventList; parentEvent = event, form = ambReachesStation, time = ambulance.route.endTime, ambulance = ambulance, station = station)
end


text:
using GPUCompiler

if !@isdefined(TestRuntime)
    include("../util.jl")
end

# create a native test compiler, and generate reflection methods for it

function bpf_job(@nospecialize(func), @nospecialize(types); kernel::Bool=false, kwargs...)
    source = FunctionSpec(func, Base.to_tuple_type(types), kernel)
    target = BPFCompilerTarget()
    params = TestCompilerParams()
    CompilerJob(target, source, params), kwargs
end

function bpf_code_llvm(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = bpf_job(func, types; kwargs...)
    GPUCompiler.compile(:llvm, job; kwargs...)
end

function bpf_code_native(@nospecialize(func), @nospecialize(types); kwargs...)
    job, kwargs = bpf_job(func, types; kwargs...)
    GPUCompiler.compile(:asm, job; kwargs...)
end


text:
include("utils.jl")
using Lazy
using Plots
using StatsBase
using Match
using LaTeXStrings

gr()
theme(:dark)

# ### Day 1 ###
s0 = [1 0]
plot(layout=@layout([a b c; e f g]))
p_samples = [(0.5, 0.5), (0.2, 0.1), (0.95, 0.7)]
for i in 1:length(p_samples)
    p, q = p_samples[i]
    P = [1-p p; q 1-q]
    states = reductions((x,y) -> x*P, s0, @lazy range(1, 20))
    plot!([x[1] for x in states], ylim=(0,1),subplot=i, label="p=$p\nq=$q")
    plot!([x[2] for x in states], ylim=(0,1), subplot=i+3, label="p=$p\nq=$q")
end

INIT_STATE = 4
REWARD = [3, 0, 0, 0, 0, 0, 0, 0, 5]
ACTIONS = [-1, 1]
STATES = 1:length(REWARD)

function env(init_state, policy, reward_calculator, state_transformer)
    s = init_state
    function play()
        sₜ  = s
        aₜ = policy(sₜ)
        rₜ = reward_calculator(sₜ, aₜ)
        s = state_transformer(sₜ, aₜ)
        (sₜ, aₜ, rₜ)  # TODO: should return a NamedTuple instead in Julia-0.7
    end
end

# get intermediate State-Action-Reward sequence
get_SAR_seq(play) = @>> begin
    play
    repeatedly
    takeuntil(res -> ((s, a, r) = res; r > 0))  # TODO: use NamedTuple here
end

random_policy(s) = rand(ACTIONS)
get_reward(s, a) = REWARD[s]
state_transform(s, a) = @match (s, a) begin
    (1, _) || (9, _) => s
    _ => s + a
end

begin
    init_play() = env(INIT_STATE, random_policy, get_reward, state_transform)
    get_steps_count() = @> init_play() get_SAR_seq length
    get_final_reward() = @> init_play() get_SAR_seq last last  # TODO: use NamedTuple here
    N = 10000
    steps_samples = convert(Array{Int}, collect(repeatedly(N, get_steps_count)))
    rewards_samples = convert(Array{Int}, collect(repeatedly(N, get_final_reward)))

    bar(counts(steps_samples, maximum(steps_samples)))
    savefig("eve_steps_count.png")
    bar(counts(rewards_samples, maximum(rewards_samples)))
    savefig("eve_reward_count.png")
    println("""Mean of steps: $(mean(steps_samples));
            Mean of rewards: $(mean(rewards_samples))""")
end
# Mean of steps: 15.9168;
# Mean of rewards: 3.78

# ### Day 2 ###
calc_gain(s, γ) = begin
    init_play() = env(s, random_policy, get_reward, state_transform)
    rewards = @>> init_play() get_SAR_seq map(res -> last(res))
    r, n = reduce((temp, r) -> ((g, n) = temp; (g + γ^n * r, n + 1)),  # TODO: use Argument destructuring here
                  (0, 0),
                  rewards)
    r
end

begin
    params = [0.95, 0.98, 1.0]
    plot(layout=@layout([a;b;c]))
    for i in 1:length(params)
        γ = params[i]
        V = zeros(length(REWARD));
        for s in STATES
            V[s] = mean(repeatedly(1000, () -> calc_gain(s, γ)))
        end
        println(V)
        bar!(V, subplot=i, label=latexstring("\\gamma = $(γ)"))
    end
    savefig("eve_v_estimate.png")
end
# [3.0, 2.40616, 2.07831, 2.00239, 2.03833, 2.32305, 2.91446, 3.70394, 5.0]
# [3.0, 2.83945, 2.75356, 2.76483, 2.93804, 3.24011, 3.64397, 4.20483, 5.0]
# [3.0, 3.254, 3.498, 3.736, 3.99, 4.252, 4.558, 4.772, 5.0]

# Value Iteration

function gen_value_iteration(γ)
    function value_iteration(values)
        function calc_V(s)
            maximum(map(a -> get_reward(s, a) + γ * values[state_transform(s, a)],
                        ACTIONS))
        end
        map(calc_V, STATES)
    end
end

function find_stable_value(γ)
    is_terminate(Vₜ, Vₜ₊₁) = maximum(abs.(Vₜ₊₁ - Vₜ)) < 0.1;
    value_iter_fn = gen_value_iteration(γ);
    gen_pair(seq) = zip(seq, tail(seq));
    init_values = zeros(Int, length(STATES))

    @>> iterate(value_iter_fn, init_values) begin
        gen_pair
        filter(pair -> is_terminate(pair...))
        first
        first
    end
end

begin
    params = [0.90, 0.95, 0.99]
    plot(layout=@layout([a;b;c]))
    for i in 1:length(params)
        γ = params[i]
        V = find_stable_value(γ)
        println(V)
        bar!(V, subplot=i, label=latexstring("\\gamma = $γ"))
    end
    savefig("eve_v_iteration.png")
end

# Policy Iteration

function policy_evaluation(V, π, γ)
    get_value(s, a) = @match (s, a) begin
        (1, -1)  => V[1]
        (9, 1) => V[9]
        _ => V[s + a]
    end

    reward = [get_reward(sa...) for sa in zip(STATES, π)]
    future_reward = γ * map(s -> get_value(s, π[s]), STATES)
    reward + future_reward
end

function policy_update(V)
    π = []
    push!(π, V[1] < V[2] ? 1 : -1)
    for i in 2:8
        push!(π, (V[i - 1] < V[i + 1]) ? 1 : -1)
    end
    push!(π, V[end] < V[end - 1] ? -1 : 1)  # ensure go right at position 5
end

function policy_iter(PV, γ)  # TODO: use argument dispatch here
    πₜ, Vₜ = PV
    println(PV)
    Vₜ₊₁ = policy_evaluation(Vₜ, πₜ, γ)
    πₜ₊₁ = policy_update(Vₜ₊₁)
    πₜ₊₁, Vₜ₊₁
end

begin
    V₀ = zeros(Float32, length(STATES))
    π₀ = [random_policy(s) for s in STATES]
    γ = 0.9  # try to change this value
    temp_PV = @>> iterate(PV -> policy_iter(PV, γ), (π₀, V₀)) take(20)
    plot()
    @gif for PV in temp_PV
        π, V = PV
        bar(V, label="V")
        for pair in zip(π, 1:9)
            dir, i = pair
            if dir == -1
                plot!([i, i - 0.5], [0,0], arrow=10, linewidth=5, color = :orange, label = "")
            else
                plot!([i, i + 0.5], [0,0], arrow=10, linewidth=5, color = :blue, label = "")
            end
        end
    end
end


text:
import NNConv1d
import NNlib
using Test
using Random
using ChainRulesTestUtils
import Flux

function to_whcn(x)
    nw,nc,nb = size(x)
    nh = 1
    reshape(x, (nw, nh, nc, nb))
end

# inp = randn(Float32, 1000, 16, 10)
# ker = randn(Float32, 5, 16, 8)
# @btime NNConv1d.conv(inp, ker)
# @btime NNlib.conv(to_whcn(inp), to_whcn(ker))

@testset "Against NNlib" begin
    for _ in 1:10
        nwin = rand(10:100)
        ncin = rand(1:8)
        nb = rand(1:5)
        ncout = rand(1:8)
        nwker = rand(1:5)
        T = rand([Float32, Float64])
        img = randn(T, nwin, ncin, nb)
        ker = randn(T, nwker, ncin, ncout)
        out_nnlib = NNlib.conv(to_whcn(img), to_whcn(ker))
        out2 = @inferred NNConv1d.conv(img, ker)
        @test eltype(out2) === T
        @test to_whcn(out2) ≈ out_nnlib
    end
end

@testset "rrule_test" begin
    for i in 1:10
        nwin = rand(10:20)
        ncin = rand(1:4)
        nb = rand(1:3)
        ncout = rand(1:4)
        nwker = rand(1:3)
        T = rand([Float32, Float64])

        x,x̄ = [randn(T, nwin, ncin, nb) for _ in 1:2]
        ker, ker̄ = [randn(T, nwker, ncin, ncout) for _ in 1:2]
        y = NNConv1d.conv(x, ker)
        NNConv1d.wcn_conv1d(nothing, x, ker)
        NNConv1d.wcn_conv1d(y, nothing, ker)
        NNConv1d.wcn_conv1d(y, x, nothing)
        ȳ = rand!(y)
        rrule_test(NNConv1d.conv, ȳ, (x,x̄), (ker, ker̄), rtol=sqrt(eps(T)))
    end
end

@testset "test Flux training" begin
    net = NNConv1d.Conv1d(3, 2 => 1, Flux.relu)
    f(x) = Flux.relu.(x[2:end-1,1:1, :])
    loss(x,y) = Flux.mse(net(x), y)
    x = randn(10, 2, 4)
    y = f(x)
    data = [(x,y) for _ in 1:10000]
    loss_before = loss(x, y)
    Flux.train!(loss, Flux.params(net), data, Flux.ADAM(1e-3))
    loss_after = loss(x, y)
    @test loss_after <= 1e-3*loss_before
end


text:
using UnitfulAtomic: austrip
using LinearAlgebra: norm
using ComponentArrays: ComponentVector

using ..InitialConditions: QuantisedDiatomic
using NQCDynamics: Estimators
using NQCModels: nstates
using NQCDistributions: Diabatic, Adiabatic

function output_template(output, u0)
    return zero(output(ComponentVector(u=[u0]), 1)[1])
end

"""
Output the end point of each trajectory.
"""
struct OutputFinal end

(::OutputFinal)(sol, i) = (last(sol.u), false)

"""
Output a 1 if the molecule has dissociated, 0 otherwise.
"""
struct OutputDissociation{T}
    "The maximum distance at which the two atoms can be considered bonded."
    distance::T
    "The indices of the two atoms in the molecule of interest."
    atom_indices::Tuple{Int,Int}
    OutputDissociation(distance, atom_indices) = new{typeof(distance)}(austrip(distance), atom_indices)
end

function (output::OutputDissociation)(sol, i)
    R = DynamicsUtils.get_positions(last(sol.u))
    dissociated = norm(R[:,output.atom_indices[1]] .- R[:,output.atom_indices[2]]) > output.distance
    return dissociated ? (1, false) : (0, false)
end

"""
Output the vibrational and rotational quantum numbers of the final image.
"""
struct OutputQuantisedDiatomic{S,H,V}
    sim::S
    height::H
    normal_vector::V
end
OutputQuantisedDiatomic(sim; height=10, normal_vector=[0, 0, 1]) = OutputQuantisedDiatomic(sim, height, normal_vector)

function (output::OutputQuantisedDiatomic)(sol, i)
    final = last(sol.u)
    ν, J = QuantisedDiatomic.quantise_diatomic(output.sim,
        DynamicsUtils.get_velocities(final), DynamicsUtils.get_positions(final);
        height=output.height, normal_vector=output.normal_vector)
    return ((ν, J), false)
end
output_template(::OutputQuantisedDiatomic, u0) = (0.0, 0.0)

"""
Output a `ComponentVector` with fields `reflection` and `transmission` containing
the probability of the outcome.
Each index in the arrays refers to the adiabatic state.
"""
struct OutputStateResolvedScattering1D{S}
    sim::S
    type::Symbol
end
function (output::OutputStateResolvedScattering1D)(sol, i)
    final = last(sol.u) # get final configuration from trajectory
    if output.type == :adiabatic
        populations = Estimators.adiabatic_population(output.sim, final)
    elseif output.type == :diabatic
        populations = Estimators.diabatic_population(output.sim, final)
    else
        throw(ArgumentError("$(output.type) not recognised.
            Only `:diabatic` or `:adiabatic` accepted."))
    end
    output = ComponentVector(
        reflection=zeros(nstates(output.sim)),
        transmission=zeros(nstates(output.sim))
    )
    x = DynamicsUtils.get_positions(final)[1]
    if x > 0 # If final position past 0 then we count as transmission
        output.transmission .= populations
    else # If final position left of 0 then we count as reflection
        output.reflection .= populations
    end
    return (output, false)
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase, StatsPlots, Dates

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# Notes:
#   - figure 1: evolution of total imports/exports from BE to UK vis-a-vis FR, NL, DE and EU
#       + HP filter available, what do chose for λ?
#   - figure 2: evolution of imports/exports from BE to UK at 2digit level
#       + HP filter available, what do chose for λ?
#   - figure 3: evolution of share of imports/exports from BE to UK of total imports/exports from BE
#       + also add FR, DE, NL, EU27
#   - table 1: evolution of most important products at 2digit level
#       + shares per product available for bilteral trade (i.e. sum to 100%) or for total trade (i.e. sum to figure 3)
#       + could compute figure 3 from data of table 1 as well!
#
#   - all figures in values and unit values (unit prices):
#       + compute a price index as using a weighted sum
#   - include vertical lines for different Brexit events

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Brexit events
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# vertical line
# vline!([Date(2016,6,23)], label="", color=:black, lw=2) # referendum
# vline!([Date(2016,12,07)], label="") # vote to trigger Article 50
# vline!([Date(2017,03,29)], label="") # Article 50
# vline!([Date(2017,06,19)], label="") # negotiations commence
# vline!([Date(2019,03,21)], label="") # first extension of Article 50
# vline!([Date(2019,04,10)], label="") # second extension of Article 50
# vline!([Date(2019,10,28)], label="") # third extension of Article 50
# vline!([Date(2019,12,20)], label="") # withdrawl bill passes UK
# vline!([Date(2020,01,31)], label="") # exit of UK from EU
# vline!([Date(2020,12,20)], label="") # UK passes Trade and Cooperation Agreement
# vline!([Date(2020,12,31)], label="") # transition period ends
# vline!([Date(2021,04,27)], label="") # EU passes Trade and Cooperation Agreement

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df_comext = CSV.read(dir_dropbox * "rawdata/" * "df_comext_BE" * ".csv", DataFrame)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 1
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = transform(df_comext, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false) # need to remove missing also from EU
subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))

# remove outliers based on medians abs dev.
# Notes:
#   - where to put MAD (need to drop EU/WORLD for sure)
outlier_cutoff = 3
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
gdf = groupby(df, "PRODUCT_NC")
df = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
subset!(df, :MAD => ByRow(x -> x < outlier_cutoff)) # lose further 1-42083391/50055769 ~16% observations

# produce data for plotting
df = data_fig1(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"])

# -----------------------

# formatting
transform!(df, [:DECLARANT_ISO, :PARTNER_ISO, :FLOW] .=> ByRow(string), renamecols=false)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
transform!(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX] .=> (x -> HP(x, λ)) .=> [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP, :PRICE_INDEX_HP])

# create unweighted unit prices
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
transform!(df, [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP] => ByRow((v,q) -> v/q) => :UNIT_PRICE_HP)

# create 3MMA as observations are very volatile
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :PRICE_INDEX => (x -> movingaverage(x,3)) => :PRICE_INDEX_3MMA)

# plotting
for flow in ["imports", "exports"]

    # values
    # Notes:
    #   - take out EU for scale
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_IN_EUROS,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="euros", title="Belgian "* flow*": values")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/" * "fig1_" * flow * "_values" * ".png") # export image dropbox

    # price index (weighted unit prices)
    # Notes:
    #   - unit prices are extremely volatile, 3 MMA would be useful
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :PRICE_INDEX,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (weighted)", title="Belgian "*flow*": price index (weighted)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/" * "fig1_" * flow * "_prices_weighted" * ".png")

    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :PRICE_INDEX_3MMA,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (weighted, 3MMA)", title="Belgian "*flow*": price index (weighted, 3MMA)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/" * "fig1_" * flow * "_prices_weighted_3MMA" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :UNIT_PRICE,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (unweighted)", title="Belgian "*flow*": price index (unweighted)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/" * "fig1_" * flow * "_prices_unweighted" * ".png")

end

# plotting of HP
for flow in ["imports", "exports"]

    # values
    # Notes:
    #   - take out EU for scale
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :VALUE_IN_EUROS_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="euros", title="Belgian "*flow*": values (HP, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/HP/" * "fig1_" * flow * "_values" * "_HP" * ".png") # export image dropbox

    # price index (weighted unit prices)
    # Notes:
    #   - unit prices are extremely volatile, 3 MMA would be useful
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :PRICE_INDEX_HP,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (weighted)", title="Belgian "*flow*": price index (weighted, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/HP/" * "fig1_" * flow * "_prices_weighted" * "_HP" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :UNIT_PRICE_HP,
        group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="price index (unweighted)", title="Belgian "*flow*": price index (unweighted, λ=$λ)")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig1/HP/" * "fig1_" * flow * "_prices_unweighted" * "_HP" * ".png")

end

# BE vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "GB")) plot(:DATE, :UNIT_PRICE_HP,
        group=:FLOW, lw=2, legend=:topleft, ylabel="euros", title="Price index (unweighted, λ=$λ): \n Belgium vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/comext/fig1/HP/" * "fig1_" * "BE_GB" * "_prices_unweighted" * "_HP" * ".png") # export image

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 4
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# rolling standard deviation
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX, :UNIT_PRICE, :VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP, :PRICE_INDEX_HP, :UNIT_PRICE_HP]
                .=> (x -> rolling_std(x, 6)) .=>
                [:STD_VALUE, :STD_QUANTITY, :STD_PRICE, :STD_UNIT, :STD_VALUE_HP, :STD_QUANTITY_HP, :STD_PRICE_HP, :STD_UNIT_HP])

# plotting HP
for flow in ["imports", "exports"]

    # values
    # Notes:
    #   - take out EU for scale
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> x != "EU")) plot(:DATE, :STD_UNIT_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="euros", title="Belgian "* flow*": std unit prices (unweighted, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig4/HP/" * "fig4_" * flow * "_std" * "_prices_unweighted" * "_HP" * ".png") # export image dropbox

end

# BE vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "GB")) plot(:DATE, :STD_UNIT_HP,
        group=:FLOW, lw=2, legend=:topleft, ylabel="euros", title="6 months STD price index (unweighted, λ=$λ): \n Belgium vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/comext/fig4/HP/" * "fig4_" * "BE_GB" * "_STD_prices_unweighted" * "_HP" * ".png") # export image

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 2
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = transform(df_comext, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false) # need to remove missing also from EU
subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))
subset!(df, :PARTNER_ISO => ByRow(x -> !(x in ["EU", "WORLD"])))

# remove outliers based on medians abs dev.
outlier_cutoff = 3
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
gdf = groupby(df, "PRODUCT_NC")
df = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
subset!(df, :MAD => ByRow(x -> x < outlier_cutoff))

# produce data for plotting
df = data_fig2(df, ["BE"], ["GB"], 2)

# ------------------

# formatting
transform!(df, [:DECLARANT_ISO, :PARTNER_ISO, :FLOW] .=> ByRow(string), renamecols=false)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

df.PRODUCT_NC_digits = parse.(Int64, String.(df.PRODUCT_NC_digits)) # transform into Int64

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PRODUCT_NC_digits"]
gdf = groupby(df, cols_grouping)
transform!(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG, :PRICE_INDEX] .=> (x -> HP(x, λ)) .=> [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP, :PRICE_INDEX_HP])

# create unweighted unit prices
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
transform!(df, [:VALUE_IN_EUROS_HP, :QUANTITY_IN_KG_HP] => ByRow((v,q) -> v/q) => :UNIT_PRICE_HP)

# plotting
# subsetting products for figures
sort!(df, :DATE) # for some reason needed to have correct DATE axis
products = 1:8:89

for flow in ["imports", "exports"]
    for i in products

    product_range = i:i+7
    end_range = string(product_range[end])

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :VALUE_IN_EUROS,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="euros", title="Belgian "*flow*": values")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig2/values/" * "fig2_" * flow * "_" * end_range * "_values" * ".png") # export image

    # price index (weighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :PRICE_INDEX,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (weighted)", title="Belgian "*flow*": unit prices (weighted)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig2/prices_weighted/" * "fig2_" * flow * "_" * end_range * "_prices_weighted" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :UNIT_PRICE,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (unweighted)", title="Belgian "*flow*": unit prices (unweighted)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig2/prices_unweighted/" * "fig2_" * flow * "_" * end_range * "_prices_unweighted" * ".png")

    end
end

for flow in ["imports", "exports"]
    for i in products

    product_range = i:i+7
    end_range = string(product_range[end])

    # values
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :VALUE_IN_EUROS_HP,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="euros", title="Belgian "*flow*": values (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig2/HP/values/" * "fig2_" * flow * "_" * end_range * "_values" * "_HP" * ".png") # export image

    # price index (weighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :PRICE_INDEX_HP,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (weighted)", title="Belgian "*flow*": unit prices (weighted, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig2/HP/prices_weighted/" * "fig2_" * flow * "_" * end_range * "_prices_weighted" * "_HP" * ".png")

    # price index (unweighted unit prices)
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PRODUCT_NC_digits => ByRow(x -> x in product_range)) plot(:DATE, :UNIT_PRICE_HP,
            group=:PRODUCT_NC_digits, lw=2, legend=:bottomleft, ylabel="unit prices (unweighted)", title="Belgian "*flow*": unit prices (unweighted, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig2/HP/prices_unweighted/" * "fig2_" * flow * "_" * end_range * "_prices_unweighted" * "_HP" * ".png")

    end
end

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 3
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = transform(df_comext, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false) # need to remove missing also from EU
subset!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> !ismissing(x)))
#subset!(df, :PARTNER_ISO => ByRow(x -> !(x in ["EU", "WORLD"])))

# remove outliers based on medians abs dev.
outlier_cutoff = 3
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)
gdf = groupby(df, "PRODUCT_NC")
df = transform(gdf, :UNIT_PRICE => MAD_method => :MAD)
subset!(df, :MAD => ByRow(x -> x < outlier_cutoff))

# produce data for plotting
df = data_fig3(df, ["BE"], ["GB", "EU", "DE", "FR", "NL"])

# ------------------

# formatting
transform!(df, [:DECLARANT_ISO, :PARTNER_ISO, :FLOW] .=> ByRow(string), renamecols=false)
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
transform!(gdf, [:VALUE_SHARE, :QUANTITY_SHARE] .=> (x -> HP(x, λ)) .=> [:VALUE_SHARE_HP, :QUANTITY_SHARE_HP])

# plotting
sort!(df, :DATE)

for flow in ["imports", "exports"]

    # value share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :VALUE_SHARE,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Belgian "*flow*": value share")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig3/" * "fig3_" * flow * "_value_share" * ".png") # export image

    # quantity share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :QUANTITY_SHARE,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="kg", title="Belgian "*flow*": quantity share")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig3/" * "fig3_" * flow * "_quantity_share" * ".png") # export image

end

# plotting HP
for flow in ["imports", "exports"]

    # value share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :VALUE_SHARE_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="percentages", title="Belgian "*flow*": value share (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig3/HP/" * "fig3_" * flow * "_value_share" * "_HP" * ".png") # export image

    # quantity share
    p = @df subset(df, :FLOW => ByRow(x -> x == flow), :PARTNER_ISO => ByRow(x -> !(x == "EU"))) plot(:DATE, :QUANTITY_SHARE_HP,
            group=:PARTNER_ISO, lw=2, legend=:bottomleft, ylabel="kg", title="Belgian "*flow*": quantity share (HP, λ=$λ)")
        vline!([Date(2016,6,23)], label="vote", color=:black, lw=1, ls=:solid) # refer
        vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
        vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/comext/fig3/HP/" * "fig3_" * flow * "_quantity_share" * "_HP" * ".png") # export image

end

# BE vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "GB")) plot(:DATE, :VALUE_SHARE_HP,
        group=:FLOW, lw=2, legend=:bottomleft, ylabel="percentages", title="Total trade share (values, λ=$λ): \n Belgium vis-a-vis Great Britain")
    vline!([Date(2016,6,23)], label="vote", color=:black, lw=1) # refer
    vline!([Date(2020,01,31)], label="exit", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="trans end", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/comext/fig3/HP/" * "fig3_" * "BE_GB" * "_value_share" * "_HP" * ".png") # export image


text:
using Term
import Term.Consoles: clear
import Term: chars
install_term_logger()

clear()

pprint(pan) = begin
    print(" " * hLine(pan.measure.w; style = "red"))
    print(vLine(pan.measure.h; style = "red") * pan)
    println(pan.measure, "  ", length(pan.segments))
    # print(pan)
end

# ---------------------------------------------------------------------------- #
#                                  no content                                  #
# ---------------------------------------------------------------------------- #
# ---------------------------------- fitted ---------------------------------- #
@time pprint(Panel(; fit = true))

# # --------------------------------- unfitted --------------------------------- #
@time pprint(Panel())

@time pprint(Panel(; width = 12, height = 4))

# ---------------------------------------------------------------------------- #
#                                   text only                                  #
# ---------------------------------------------------------------------------- #
# ---------------------------------- fitted ---------------------------------- #
@time pprint(Panel("t"; fit = true))

@time pprint(Panel("test"; fit = true))

@time pprint(Panel("1234\n123456789012"; fit = true))

@time pprint(Panel("나랏말싸미 듕귁에 달아"; fit = true))

@time pprint(Panel("나랏말싸미 듕귁에 달아\n1234567890123456789012"; fit = true))

@time pprint(Panel("."^500; fit = true))

# --------------------------------- unfitted --------------------------------- #
@time pprint(Panel("t"))

@time pprint(Panel("test"))

@time pprint(Panel("1234\n123456789012"))

@time pprint(Panel("나랏말싸미 듕귁에 달아\n1234567890123456789012"))

@time pprint(Panel("."^500))

# ---------------------------------------------------------------------------- #
#                                 nested panels                                #
# ---------------------------------------------------------------------------- #
# ---------------------------------- fitted ---------------------------------- #
@time pprint(Panel(Panel("test"; fit = true); fit = true))

@time pprint(Panel(Panel(Panel("."; fit = true); fit = true); fit = true))

@time pprint(Panel(Panel("."^500; fit = true); fit = true))

# --------------------------------- unfitted --------------------------------- #
@time pprint(Panel(Panel("test"); fit = true))

@time pprint(Panel(Panel(Panel("."); fit = true); fit = true))

@time pprint(Panel(Panel("."^250); fit = true))

@time pprint(Panel(Panel("test");))

@time pprint(Panel(Panel(Panel("."););))

@time pprint(Panel(Panel("."^250);))

@time pprint(Panel(Panel("t1"), Panel("t2")))

@time pprint(Panel(Panel("test"; width = 22); width = 30, height = 8))

@time pprint(Panel(Panel("test"; width = 42); width = 30, height = 8))

@time pprint(Panel(Panel("test"; width = 42, height = 12); width = 30, height = 8))

@time pprint(Panel(Panel("test"; width = 42, height = 12); width = 30, height = 14))

# -------------------------- with other renderables -------------------------- #
@time pprint(Panel(RenderableText("x" .^ 5)))

@time pprint(Panel(RenderableText("x" .^ 500)))

@time pprint(Panel(RenderableText("x" .^ 5); fit = true))

@time pprint(Panel(RenderableText("x" .^ 500); fit = true))

# ---------------------------------------------------------------------------- #
#                                    titles                                    #
# ---------------------------------------------------------------------------- #

@time pprint(Panel(Panel("[red].[/red]"^50; title = "test", subtitle = "subtest")))

@time pprint(
    Panel(
        "."^50;
        title = "test",
        subtitle = "subtest",
        subtitle_style = "red",
        fit = false,
    ),
)

# ---------------------------------------------------------------------------- #
#                                    PADDING                                   #
# ---------------------------------------------------------------------------- #

@time pprint(Panel("."^24; padding = (4, 4, 2, 2)))

@time pprint(Panel("."^24; padding = (4, 4, 2, 2), fit = true))

@time pprint(Panel("."^24; padding = (0, 0, 0, 0), fit = true))


text:
using Test
using LibASICamera

devices = []
try
    devices = get_connected_devices()
catch err
    println("\nConnect the camera for testing and make sure you can access the camera
    without being root by calling \"sudo install asi.rules /lib/udev/rules.d\"
    from the lib subdir of the SDK and then relogging / rebooting.\n")
    rethrow(err)
end

@testset "$(cam.info.Name)" for cam in devices
    control_caps = get_control_caps(cam)
    @testset "set/get $(cap.Name)" for cap in control_caps
        default = cap.DefaultValue
        control_type = cap.ControlType
        if cap.IsWritable == ASI_TRUE
            set_control_value(cam, control_type, default, auto = false)
            @test default == get_control_value(cam, control_type)[1]
        end
    end

    @testset "Capture Still" begin
        @test isa(capture_still(cam), Array)
    end

    close_camera(cam)
end


text:
using Statistics: mean

@testset "sparsemax integer, float input" begin
    @test sparsemax(Int[0 0; 0 0]) == [0.5 0.5; 0.5 0.5]
    @test sparsemax(Int[1 0; 0 1]) == [1.0 0.0; 0.0 1.0]
    @test sparsemax(Int[1 2; 1 1]) == [0.5 1.0; 0.5 0.0]

    @test sparsemax(Float32[0 0; 0 0]) == Float32[0.5 0.5; 0.5 0.5]
    @test sparsemax(Float64[0 0; 0 0]) == Float64[0.5 0.5; 0.5 0.5]
end

@testset "sparsemax vector input" begin
    @test sparsemax(Float32[1, 2, 2]; dims=1) == Float32[0.0, 0.5, 0.5]
    @test sparsemax(Float32[1, 2, 2]; dims=2) == Float32[1.0, 1.0, 1.0]
end

@testset "sparsemax inferred return type" begin
    for T in [Float16, Float32, Float64]
        val = @inferred sparsemax(T[1,2,3])
        @test typeof(val) == Vector{T}
    end
end

@testset "AutoDiff" begin
    gradtest(x -> sparsemax(x) .* (1:3), 3)
    gradtest(x -> sparsemax(x) .* (1:3), (3,5), atol = 1e-4)
    gradtest(x -> sparsemax(x, dims = 2) .* (1:3), (3,5), atol = 1e-4)

    for d in (1, 2)
        gradtest(sparsemax, (3,4), fkwargs = (dims = d,))
    end
end


text:
# Here we implement the banded matrix interface for some key examples
isbanded(::Zeros) = true
bandwidths(::Zeros) = (0,0)
inbands_getindex(::Zeros{T}, k::Integer, j::Integer) where T = zero(T)

isbanded(::Eye) = true
bandwidths(::Eye) = (0,0)
inbands_getindex(::Eye{T}, k::Integer, j::Integer) where T = one(T)

isbanded(::Diagonal) = true
bandwidths(::Diagonal) = (0,0)
inbands_getindex(D::Diagonal, k::Integer, j::Integer) = D.diag[k]
inbands_setindex!(D::Diagonal, v, k::Integer, j::Integer) = (D.diag[k] = v)

isbanded(::SymTridiagonal) = true
bandwidths(::SymTridiagonal) = (1,1)
inbands_getindex(J::SymTridiagonal, k::Integer, j::Integer) =
    k == j ? J.dv[k] : J.ev[k]
inbands_setindex!(J::SymTridiagonal, v, k::Integer, j::Integer) =
    k == j ? (J.dv[k] = v) : (J.ev[k] = v)

isbanded(K::Kron{<:Any,2}) = all(isbanded, K.arrays)
function bandwidths(K::Kron{<:Any,2})
    A,B = K.arrays
    (size(B,1)*bandwidth(A,1) + max(0,size(B,1)-size(B,2))*size(A,1)   + bandwidth(B,1),
        size(B,2)*bandwidth(A,2) + max(0,size(B,2)-size(B,1))*size(A,2) + bandwidth(B,2))
end

const BandedMatrixTypes = (:AbstractBandedMatrix, :(AdjOrTrans{<:Any,<:AbstractBandedMatrix}),
                                    :(AbstractTriangular{<:Any, <:AbstractBandedMatrix}),
                                    :(Symmetric{<:Any, <:AbstractBandedMatrix}),
                                    :Zeros, :Eye, :Diagonal, :SymTridiagonal)

for T1 in BandedMatrixTypes, T2 in BandedMatrixTypes
    @eval kron(A::$T1, B::$T2) = BandedMatrix(Kron(A,B))
end


text:
using Pluto
Pluto.run("0.0.0.0.0", 1234; configuration=Pluto.ServerConfiguration(launch_browser=false))


text:
module OpticalFlowUtils

export colorflow, angular_error, endpoint_error

include("io.jl")
include("colorflow.jl")
include("measures.jl")

end # module


text:
"""
    finddatasets(name=nothing, blocks=Any)

Find preconfigured dataset recipes for datasets that match block
types `blocks` in all data sources (if `name == nothing`) or dataset source
with `name`. `blocks` can be given as a type or a nested tuple of block types.

!!! warning "Deprecated"

    This function is deprecated and will be removed in a future version
    of FastAI.jl. Use `filter(datarecipes(); blocks, id = name)` instead.

Return a vector of `Pair`s `datasetname => recipe`

#### Examples

Loading a result

```julia
datasetname, recipe = finddatasets(blocks=(Image, Label))[1]
data, blocks = loadrecipe(recipe, datasetpath(datasetname))
```

Example searches

```
# Single-label image classification
finddatasets(blocks=(Image, Label))

# Single-label classification from any data
finddatasets(blocks=(Any, Label))

# Datasets with images as input data
finddatasets(blocks=(Image, Any))

# All ways to load `pascal2007`
finddatasets(name="pascal2007")
```
"""
finddatasets(; name = "", blocks = Any) = datarecipes(; blocks, id = name)

"""
    loaddataset(name[, blocks = Any]) -> (data, blocks)

Load dataset `name` with a recipe that matches block types
`blocks`. The first matching recipe is selected and loaded.

!!! warning "Deprecated"

    This function is deprecated and will be removed in a future version
    of FastAI.jl. Use `findfirst(datarecipes(); name, blocks)` instead.

## Examples

Load a data container suitable for single-label image classification:
```julia
data, blocks = loaddataset("imagenette2-160", (Image, Label))
```

Load dataset with any recipe:
```julia
data, blocks = loaddataset(name)
```
"""
loaddataset(name::String, blocks = Any) = load(findfirst(datarecipes(); id = name, blocks))

"""
    listdatasources()

!!! warning "Deprecated"

    This function is deprecated and will be removed in a future version
    of FastAI.jl.

List the dataset sources registered in `registry` (defaults to
`FastAI.defaultdataregistry()`).
"""
listdatasources() = getfield(datasets(), :data).id

"""
    datasetpath(name)

(Down)load registered dataset source named `name` from the dataset registry.
Use [`listdatasources`](#) for a list of all dataset sources.
"""
datasetpath(name::String) = load(datasets()[name])


text:
# for pomcp solver
function sparse_actions(pomdp::CommunityAffiliationPOMDP, s::CommunityNetwork, b::RevealedGraph, n::Int)
    @assert n <= 0 # doesn't work with a limited number of actions
    if s.budget > 0
        # println("ProbeNode actions (budget=$(s.budget))")
        return [ProbeNode(i) for i in setdiff(b.revealed_nodes,b.probed)]
    else
        # println("GuessAffiliation actions (budget=$(s.budget))")
        return [GuessAffiliation(i) for i in 1:b.n_comms]
    end
end


text:
using LOLTools.MatchV4

api_key = get(ENV, "RIOT_TOKEN", "")
platform = "euw1"
match = MatchV4.match_by_id(api_key, platform, 4249110676)
@info match.participantIdentities[10].player.summonerName == "FPX Doinb"


text:
using BinDeps

@BinDeps.setup

libpppack = library_dependency("libpppack")

prefix = joinpath(BinDeps.depsdir(libpppack), "usr")
src_dir = joinpath(BinDeps.depsdir(libpppack), "src", "pppack")
build_dir = joinpath(BinDeps.depsdir(libpppack), "builds", "pppack")

provides(SimpleBuild,
         (@build_steps begin
	  `git submodule update --init --recursive`
          CreateDirectory(build_dir)
          @build_steps begin
          ChangeDirectory(build_dir)
          @build_steps begin
          `cmake -DCMAKE_INSTALL_PREFIX="$prefix" $src_dir`
          `make`
          `make install`
          end
          end
          end),
         libpppack)
@BinDeps.install Dict([(:libpppack, :_jl_libpppack)])


text:
# This file is a part of GERDADeepLearning.jl, licensed under the MIT License (MIT).

export mapvalues
function mapvalues(d::Dict, f, post_params...)
  result = copy(d)
  for (key,value) in d
    result[key] = f(value, post_params...)
  end
  return result
end

function mapvalues(data::DLData, f, post_params...)
  return DLData([f(lib, post_params...) for lib in data])
end


text:
const thirty_days = 30days
const Nmonths = 12

#@inline current_time_index(time, interval=2592000, length=12) = mod(unsafe_trunc(Int32, time / interval),     length) + 1
#@inline next_time_index(time, interval=2592000, length=12)    = mod(unsafe_trunc(Int32, time / interval) + 1, length) + 1
#@inline cyclic_interpolate(u₁::Number, u₂, time, interval=2592000) = u₁ + mod(time / interval, 1) * (u₂ - u₁)

@inline current_time_index(time) = mod(unsafe_trunc(Int32, time / 2592000), 12) + 1
@inline next_time_index(time) = mod(unsafe_trunc(Int32, time / 2592000) + 1, 12) + 1
@inline cyclic_interpolate(u₁::Number, u₂, time) = u₁ + mod(time / 2592000, 1) * (u₂ - u₁)

@inline function cyclic_interpolate(τ::AbstractArray, time, interval=2592000, length=12)
    n₁ = current_time_index(time, interval, length)
    n₂ = next_time_index(time, interval, length)
    return cyclic_interpolate.(view(τ, :, :, n₁), view(τ, :, :, n₂), time, interval)
end


text:
cd(@__DIR__); include("setups/simpletree.jl");
gr(dpi=200)

for l = 97:100
    plot(size = (200, 500), framestyle = :none)
    gplot!(W, X, width = 1)
    scatter_gplot!(X; marker = 𝚽[:, l])
    plt = plot!(cbar = false, clims = (-0.3, 0.3))
    savefig(plt, "../figs/SimpleTree_eigenvector$(l-1)_red.png")
end


text:
function readcell(cell_filename)
    cellstring = open(cell_filename,"r") do f
        cellstring = read(f,String)
    end

    cell_object = JSON2.read(cellstring,struct_cell_general)

    return cell_object
end

function writecell(cell_filename,cell)
    cellstring = JSON2.write(cell)

    return_code = open(cell_filename,"w") do f
        write(f,cellstring)
    end

    return return_code
end

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
# Aesthetics is a set of bindings of typed values to symbols (Wilkinson calls
# this a Varset). Each variable controls how geometries are realized.
type Aesthetics
    x::Union(Nothing, Vector{Float64}, Vector{Int64})
    y::Union(Nothing, Vector{Float64}, Vector{Int64})
    size::Maybe(Vector{Measure})
    color::Maybe(AbstractDataVector{ColorValue})
    label::Maybe(PooledDataVector)

    x_min::Union(Nothing, Vector{Float64}, Vector{Int64})
    x_max::Union(Nothing, Vector{Float64}, Vector{Int64})
    y_min::Union(Nothing, Vector{Float64}, Vector{Int64})
    y_max::Union(Nothing, Vector{Float64}, Vector{Int64})

    # Boxplot aesthetics
    middle::Maybe(Vector{Float64})
    lower_hinge::Maybe(Vector{Float64})
    upper_hinge::Maybe(Vector{Float64})
    lower_fence::Maybe(Vector{Float64})
    upper_fence::Maybe(Vector{Float64})
    outliers::Maybe(Vector{Vector{Float64}})

    # Aesthetics pertaining to guides
    xtick::Maybe(Vector{Float64})
    ytick::Maybe(Vector{Float64})
    xgrid::Maybe(Vector{Float64})
    ygrid::Maybe(Vector{Float64})

    # Pesudo-aesthetics used to indicate that drawing might
    # occur beyond any x/y value.
    x_drawmin::Maybe(Float64)
    x_drawmax::Maybe(Float64)
    y_drawmin::Maybe(Float64)
    y_drawmax::Maybe(Float64)

    # Plot viewport extents
    x_viewmin::Maybe(Float64)
    x_viewmax::Maybe(Float64)
    y_viewmin::Maybe(Float64)
    y_viewmax::Maybe(Float64)

    color_key_colors::Maybe(Vector{ColorValue})
    color_key_title::Maybe(String)
    color_key_continuous::Maybe(Bool)

    # Labels. These are not aesthetics per se, but functions that assign lables
    # to values taken by aesthetics. Often this means simply inverting the
    # application of a scale to arrive at the original value.
    x_label::Function
    y_label::Function
    xtick_label::Function
    ytick_label::Function
    color_label::Function

    function Aesthetics()
        new(nothing, nothing, nothing, nothing,
            nothing, nothing, nothing, nothing,
            nothing, nothing, nothing, nothing,
            nothing, nothing, nothing, nothing,
            nothing, nothing, nothing, nothing,
            nothing, nothing, nothing, nothing,
            nothing, nothing, nothing, nothing,
            nothing, nothing,
            fmt_float, fmt_float,
            string, string, string)
    end

    # shallow copy constructor
    function Aesthetics(a::Aesthetics)
        b = new()
        for name in Aesthetics.names
            setfield(b, name, getfield(a, name))
        end
        b
    end
end

# Index as if this were a data frame
function getindex(aes::Aesthetics, i::Integer, j::String)
    getfield(aes, symbol(j))[i]
end

# Return the set of variables that are non-nothing.
function defined_aesthetics(aes::Aesthetics)
    vars = Set{Symbol}()
    for name in Aesthetics.names
        if !is(getfield(aes, name), nothing)
            add!(vars, name)
        end
    end
    vars
end

# Checking aesthetics and giving reasonable error messages.

# Raise an error if any of thu given aesthetics are not defined.
#
# Args:
#   who: A string naming the caller which is printed in the error message.
#   aes: An Aesthetics object.
#   vars: Symbol that must be defined in the aesthetics.
#
function assert_aesthetics_defined(who::String, aes::Aesthetics, vars::Symbol...)
    undefined_vars = setdiff(Set(vars...), defined_aesthetics(aes))
    if !isempty(undefined_vars)
        error(@sprintf("The following aesthetics are required by %s but are not defined: %s\n",
                       who, join(undefined_vars, ", ")))
    end
end

function assert_aesthetics_equal_length(who::String, aes::Aesthetics, vars::Symbol...)
    defined_vars = Symbol[]
    for var in filter(var -> !(getfield(aes, var) === nothing), vars)
        push!(defined_vars, var)
    end

    n = length(getfield(aes, vars[1]))
    for i in 2:length(vars)
        if length(getfield(aes, vars[1])) != n
            error(@sprintf("The following aesthetics are required by %s to be of equal length: %s\n",
                           who, join(vars, ", ")))
        end
    end
end

# Create a shallow copy of an Aesthetics instance.
#
# Args:
#   a: aesthetics to copy
#
# Returns:
#   Copied aesthetics.
#
copy(a::Aesthetics) = Aesthetics(a)

# Replace values in a with non-nothing values in b.
#
# Args:
#   a: Destination.
#   b: Source.
#
# Returns: nothing
#
# Modifies: a
#
function update!(a::Aesthetics, b::Aesthetics)
    for name in Aesthetics.names
        if issomething(getfield(b, name))
            setfield(a, name, getfield(b, name))
        end
    end

    nothing
end

# Serialize aesthetics to JSON.

# Args:
#  a: aesthetics to serialize.
#
# Returns:
#   JSON data as a string.
#
function json(a::Aesthetics)
    join([string(a, ":", json(getfield(a, var))) for var in aes_vars], ",\n")
end

# Concatenate aesthetics.
#
# A new Aesthetics instance is produced with data vectors in each of the given
# Aesthetics concatenated, nothing being treated as an empty vector.
#
# Args:
#   aess: One or more aesthetics.
#
# Returns:
#   A new Aesthetics instance with vectors concatenated.
#
function cat(aess::Aesthetics...)
    cataes = Aesthetics()
    for aes in aess
        for var in Aesthetics.names
            setfield(cataes, var,
                     cat_aes_var!(getfield(cataes, var), getfield(aes, var)))
        end
    end
    cataes
end

cat_aes_var!(a::Nothing, b::Nothing) = a
cat_aes_var!(a::Nothing, b) = copy(b)
cat_aes_var!(a, b::Nothing) = a
cat_aes_var!(a::Function, b::Function) = a === string || a === fmt_float ? b : a
function cat_aes_var!(a, b)
    append!(a, b)
    a
end

function cat_aes_var!{T}(xs::PooledDataVector{T}, ys::PooledDataVector{T})
    newpool = T[x for x in union(Set(xs.pool...), Set(ys.pool...))]
    newdata = vcat(T[x for x in xs], T[y for y in ys])
    PooledDataArray(newdata, newpool, [false for _ in newdata])
end