text:
abstract type AbstractParameter end

"""
    value(x)

Return the "value" of an object.
For `AbstractParameter`s this typically applies some transformation to some data
contained in the parameter, and returns a plain data type.
It might, for example, return a transformation of some internal data, the result of which
is guaranteed to satisfy some contraint.
"""
value(x)

# Various basic `value` definitions.
value(x::Number) = x
value(x::AbstractArray{<:Number}) = x
value(x::AbstractArray) = map(value, x)
value(x::Tuple) = map(value, x)
value(x::NamedTuple) = map(value, x)
value(x::Dict) = Dict(k => value(v) for (k, v) in x)

"""
    positive(val::T, transform::Bijector=Bijectors.Exp(), ε=sqrt(eps(T))) where T<:Real

Returns a `Postive`.
The `value` of a `Positive` is a `Real` number that is constrained to be positive.
This is represented in terms of an a `transform` that maps an `unconstrained_value` to the
positive reals.
Satisfies `val ≈ transform(unconstrained_value)`
"""
function positive(
    val::T, transform::Bijector=Bijectors.Exp(), ε=sqrt(eps(T)),
) where T<:Real
    val > 0 || throw(ArgumentError("Value ($val) is not positive."))
    val > ε || throw(ArgumentError("Value ($val) is too small, relative to ε ($ε)."))
    unconstrained_value = inv(transform)(val - ε)
    return Positive(unconstrained_value, transform, convert(typeof(unconstrained_value), ε))
end

struct Positive{T<:Real, V<:Bijector, Tε<:Real} <: AbstractParameter
    unconstrained_value::T
    transform::V
    ε::Tε
end

value(x::Positive) = x.transform(x.unconstrained_value) + x.ε

function flatten(::Type{T}, x::Positive) where T<:Real
    v, unflatten_to_Real = flatten(T, x.unconstrained_value)

    function unflatten_Positive(v_new::Vector{T})
        return Positive(unflatten_to_Real(v_new), x.transform, x.ε)
    end

    return v, unflatten_Positive
end

"""
    bounded(val::Real, lower_bound::Real, upper_bound::Real)

Constructs a `Bounded`.
The `value` of a `Bounded` is a `Real` number that is constrained to be within the interval
(`lower_bound`, `upper_bound`), and is equal to `val`.
This is represented internally in terms of an `unconstrained_value` and a `transform` that
maps any real to this interval. `unconstrained_value` is `inv(transform)(val)`.
"""
function bounded(val::Real, lower_bound::Real, upper_bound::Real)
    lb = convert(typeof(val), lower_bound)
    ub = convert(typeof(val), upper_bound)
    ε = convert(typeof(val), 1e-12)

    if val > upper_bound || val < lower_bound
        throw(ArgumentError(
            "Value, $val, outside of specified bounds ($lower_bound, $upper_bound).",
        ))
    end

    inv_transform = Bijectors.Logit(lb + ε, ub - ε)
    transform = inv(inv_transform)

    # Bijectors defines only Logit struct so we use Logistic as the inverse of Logit
    return Bounded(inv_transform(val), lb, ub, transform, ε)
end

struct Bounded{T<:Real, V<:Bijector, Tε<:Real} <: AbstractParameter
    unconstrained_value::T
    lower_bound::T
    upper_bound::T
    transform::V
    ε::Tε
end

value(x::Bounded) = x.transform(x.unconstrained_value)

function flatten(::Type{T}, x::Bounded) where T<:Real
    v, unflatten_to_Real = flatten(T, x.unconstrained_value)

    function unflatten_Bounded(v_new::Vector{T})
        return Bounded(
            unflatten_to_Real(v_new), x.lower_bound, x.upper_bound, x.transform, x.ε,
        )
    end

    return v, unflatten_Bounded
end

"""
    fixed(val)

Represents a parameter whose value is required to stay constant. The `value` of a `Fixed` is
simply `val`. Constantness of the parameter is enforced by returning an empty
vector from `flatten`.
"""
fixed(val) = Fixed(val)

struct Fixed{T} <: AbstractParameter
    value::T
end

value(x::Fixed) = x.value

function flatten(::Type{T}, x::Fixed) where T<:Real
    unflatten_Fixed(v_new::Vector{T}) = x
    return T[], unflatten_Fixed
end

"""
    deferred(f, args...)

The `value` of a `deferred` is `f(value(args)...)`. This makes it possible to make the value
of the `args` e.g. `AbstractParameter`s and, therefore, enforce constraints on them even if
`f` knows nothing about `AbstractParameters`.

It can be helpful to use `deferred` recursively when constructing complicated objects.
"""
deferred(f, args...) = Deferred(f, args)

struct Deferred{Tf, Targs} <: AbstractParameter
    f::Tf
    args::Targs
end

Base.:(==)(a::Deferred, b::Deferred) = (a.f == b.f) && (a.args == b.args)

value(x::Deferred) = x.f(value(x.args)...)

function flatten(::Type{T}, x::Deferred) where T<:Real
    v, unflatten = flatten(T, x.args)
    unflatten_Deferred(v_new::Vector{T}) = Deferred(x.f, unflatten(v_new))
    return v, unflatten_Deferred
end

"""
    nearest_orthogonal_matrix(X::StridedMatrix)

Project `X` onto the closest orthogonal matrix in Frobenius norm.

Originally used in varz: https://github.com/wesselb/varz/blob/master/varz/vars.py#L446
"""
@inline function nearest_orthogonal_matrix(X::StridedMatrix{<:Union{Real, Complex}})
    # Inlining necessary for type inference for some reason.
    U, _, V = svd(X)
    return U * V'
end

"""
    orthogonal(X::StridedMatrix{<:Real})

Produce a parameter whose `value` is constrained to be an orthogonal matrix. The argument `X` need not
be orthogonal.

This functionality projects `X` onto the nearest element subspace of orthogonal matrices (in
Frobenius norm) and is overparametrised as a consequence.

Originally used in varz: https://github.com/wesselb/varz/blob/master/varz/vars.py#L446
"""
orthogonal(X::StridedMatrix{<:Real}) = Orthogonal(X)

struct Orthogonal{TX<:StridedMatrix{<:Real}} <: AbstractParameter
    X::TX
end

Base.:(==)(X::Orthogonal, Y::Orthogonal) = X.X == Y.X

value(X::Orthogonal) = nearest_orthogonal_matrix(X.X)

function flatten(::Type{T}, X::Orthogonal) where {T<:Real}
    v, unflatten_to_Array = flatten(T, X.X)
    unflatten_Orthogonal(v_new::Vector{T}) = Orthogonal(unflatten_to_Array(v_new))
    return v, unflatten_Orthogonal
end


text:
# target a univariate measure μ(dx) ∝ (x-\mu)^2/2\sigma^2 + δ_0 (x) dx
using ZigZagBoomerang
const ZZB = ZigZagBoomerang

###Boomerang freexing time with μstar = 0.0
function ZZB.freezing_time(x, θ)
    if θ*x >= 0.0
        return π - atan(x/θ)
    else
        return atan(-x/θ)
    end
end
###Boomerang freexing time with μstar = 0.0
function ZZB.freezing_time(x, θ, F::Boomerang1d)
    if F.μ == 0.0
        if x == 0.0
            return 2π
        else
            return ZZB.freezing_time(x, θ)
        end
    elseif abs(F.μ) - sqrt((x - F.μ)^2 + θ^2) >= 0.0
            return Inf
    end
    a = acos(F.μ/(sqrt((x - F.μ)^2 + θ^2)))
    b = atan(θ/(x - F.μ))
    if (x - F.μ) < 0.0
        if x < 0.0  #  x'<0, x<0
            return a + b
        elseif θ > 0 # x'<0, x>0, θ>0
            return 2π - a + b
        else   # x'<0, x>0, θ<0
            return - a + b
        end
    elseif x > 0   # x' > 0, x > 0
        return π - a + b
    elseif θ < 0   # x' > 0, x < 0, θ < 0
        return π + a + b
    else    # x' > 0, x < 0, θ > 0
        return -π + a + b
    end
end

function freezing_time_at0(v, B::Boomerang1d)
    if B.μ == 0.0
        return pi
    end
    a = acos(B.μ/(sqrt((B.μ)^2 + v^2)))
    b = atan(v/(-B.μ))
    if B.μ > 0
        if v > 0
            return 2π - a + b
        else
            return a + b
        end
    elseif v > 0
        return π - a + b
    else
        return π + a + b
    end
end

# k determines the weight on the Dirac measure. The smaller k, the higher the wegiht
function ss_pdmp(∇ϕ, x, θ, T, c, k, Flow::ZZB.ContinuousDynamics; adapt=false, factor=2.0)
    t = zero(T)
    Ξ = [(t, x, θ)]
    num = acc = 0
    a, b = ZZB.ab(x, θ, c, Flow)
    t_ref = t + ZZB.waiting_time_ref(Flow)
    t′ =  t + poisson_time(a, b, rand())
    tˣ = t + ZZB.freezing_time(x, θ, Flow)
    while t < T
        if  tˣ < min(t_ref, t′) # it s freezing
            t, x, θ = ZZB.move_forward(tˣ - t, t, x, θ, Flow) # go to 0
            @assert -0.0001 < x < 0.0001 #check
            #t, x , θ = tˣ, 0.0, θ #go to 0
            push!(Ξ, (t, x, 0.0))
            t = t - log(rand())/k #wait exponential time
            if t_ref < t
                θ = sqrt(Flow.Σ)*randn()
                t_ref = t +  ZZB.waiting_time_ref(Flow)
            end
            push!(Ξ, (t, x, θ))
            tˣ  = t + freezing_time_at0(θ, Flow)
        elseif t_ref < t′
            t, x, θ = ZZB.move_forward(t_ref - t, t, x, θ, Flow)
            θ = sqrt(Flow.Σ)*randn()
            t_ref = t +  ZZB.waiting_time_ref(Flow)
            tˣ = t + ZZB.freezing_time(x, θ, Flow)
            push!(Ξ, (t, x, θ))
        else
            τ = t′ - t
            t, x, θ = ZZB.move_forward(τ, t, x, θ, Flow)
            l, lb = ZZB.λ(∇ϕ, x, θ, Flow), ZZB.λ_bar(τ, a, b)
            num += 1
            if rand()*lb < l
                acc += 1
                if l >= lb
                    !adapt && error("Tuning parameter `c` too small.")
                    c *= factor
                end
                θ = -θ  # In multi dimensions the change of velocity is different:
                        # reflection symmetric on the normal vector of the contour
                tˣ = t + ZZB.freezing_time(x, θ, Flow)
                push!(Ξ, (t, x, θ))
            end
        end
        a, b = ZZB.ab(x, θ, c, Flow)
        t′ = t + poisson_time(a, b, rand())
    end
    return Ξ, acc/num
end

ϕ(x) = x^2/2 #not used
∇ϕ(x) = x

# Example: ZigZag
x0, θ0 = -randn(), randn()
T = 1000.0
c = 0.0
k = 0.5
μ = 1.0
ref = 0.5
B = Boomerang1d(μ, ref)
c = 0.0
out1, acc = ss_pdmp(∇ϕ, x0, θ0, T, c, k, B)
out1
dt = 0.1
xx = ZigZagBoomerang.discretize(out1, B, dt)
using Makie, CairoMakie
using Plots
p2 = Makie.lines(xx.t, xx.x, linewidth=0.4)

## TEST probability of x being 0 euqal to total time divided by time spent on 0
function prob_of_zero(out1)
    k = 0
    freeze =false
    for event in out1
        if freeze == true
            k += event[1]
        end
        if event[3] == 0.0
            freeze = true
            k -= event[1]
        else
            freeze = false
        end
    end
    k/out1[end][1]
end
prob_of_zero(out1)

A = zeros(1, 200)
C = zeros(1, 200)
μ, k = 0.0, 0.01
    for i in 1:1
        global μ += 0.0
        for j in 1:200
            global k += 0.03
            x0, θ0 = -randn(), randn()
            T = 100000.0
            c = 0.0
            ref = 0.5
            B = Boomerang1d(μ, ref)
            c = 0.0
            out1, acc = ss_pdmp(∇ϕ, x0, θ0, T, c, k, B)
            A[i,j] = prob_of_zero(out1)
            #C[i,j] = 1/(1 + k*sqrt(2*pi)*exp(0.5*(-μ)^2))
            C[i,j] = 1/(1 + k*sqrt(2pi)*exp(0.5*(-μ)^2)/sqrt(2/pi))
        end
    end
b = sum([(A[1,j] - sum(A[1,:])/200)*(C[1,j] - sum(C[1,:])/200) for j in 1:200])/
            sum([(A[1,j] - sum(A[1,:])/200)^2 for j in 1:200])
a = sum(C[1,:])/200 - sum(A[1,:])/200*b
#b1 =  sum([A[1,j]*C[1,j] for j in 1:100])/
#            sum([(A[1,j])^2 for j in 1:100])
pp = Makie.plot([A[1,j] for j in 1:200], [C[1,j] for j in 1:200], markersize=3)
Makie.lines!(0.1:0.1:1.0, [a + b*x for x in 0.1:0.1:1.0])
#Makie.lines!(0.1:0.1:0.7, [b1*x for x in 0.1:0.1:0.7])
save("figures/boomerang_prob0.png", title(pp, "Boomerang 1d, varying k"))

for i in 2:10
    Makie.plot!(a, [A[i,j] for j in 1:10], [C[i,j] for j in 1:10])
end
a

ext: jl


text:
module TestAqua

using Aqua
using GeneratorsX

Aqua.test_all(
    GeneratorsX;
    ambiguities = false,  # upstream
    project_extras = true,
    stale_deps = true,
    deps_compat = true,
    project_toml_formatting = true,
)

end  # module


text:
using UnitfulPlasmas
using Base.Test


text:
function includeall(dir)
    for (root,dirs,files) in walkdir(dir)
        for file in files
            if file[end-2:end] == ".jl" && file[1:4]!="Test"
                include(joinpath(root,file))
            end
        end
    end
end

@inline function evalpoly(x::N,coefs::T) where {T<:AbstractArray,N<:Number}
    out = 0.0
    for i in length(coefs):-1:1
        @inbounds out += coefs[i]
        if i != 1
            @fastmath out *= x
        end
    end
    return out
end

export evalpoly

@inline function getBoundingIndsSorted(el::Q,x::T) where {T<:AbstractArray,Q<:Any}
    if el <= x[1]
        return (1,-1)
    end
    for i in 1:length(x)
        if @inbounds x[i] >= el
            if @inbounds x[i] == el
                return (i,-1)
            else
                return (i-1,i)
            end
        end
    end
    return (length(x),-1)
end
export getBoundingIndsSorted


text:
# Tests on Multivariate LogNormal distributions

using Distributions, Compat
using Base.Test

####### Core testing procedure

function test_mvlognormal(g::MvLogNormal, n_tsamples::Int=10^6)
    d = length(g)
    mn = mean(g)
    md = median(g)
    mo = mode(g)
    S = cov(g)
    s = var(g)
    e = entropy(g)
    @test partype(g) == Float64
    @test isa(mn, Vector{Float64})
    @test isa(md, Vector{Float64})
    @test isa(mo, Vector{Float64})
    @test isa(s, Vector{Float64})
    @test isa(S, Matrix{Float64})
    @test length(mn) == d
    @test length(md) == d
    @test length(mo) == d
    @test length(s) == d
    @test size(S) == (d, d)
    @test_approx_eq s diag(S)
    @test_approx_eq md @compat(exp.(mean(g.normal)))
    @test_approx_eq mn @compat(exp.(mean(g.normal) + var(g.normal)/2))
    @test_approx_eq mo @compat(exp.(mean(g.normal) - var(g.normal)))
    @test_approx_eq entropy(g) d*(1 + Distributions.log2π)/2 + logdetcov(g.normal)/2 + sum(mean(g.normal))
    gg = typeof(g)(MvNormal(params(g)...))
    @test full(g.normal.μ) == full(gg.normal.μ)
    @test full(g.normal.Σ) == full(gg.normal.Σ)
    @test insupport(g,ones(d))
    @test !insupport(g,zeros(d))
    @test !insupport(g,-ones(d))

    # sampling
    X = rand(g, n_tsamples)
    emp_mn = vec(mean(X, 2))
    emp_md = vec(median(X, 2))
    Z = X .- emp_mn
    emp_cov = A_mul_Bt(Z, Z) * (1.0 / n_tsamples)
    for i = 1:d
        @test_approx_eq_eps emp_mn[i] mn[i] (sqrt(s[i] / n_tsamples) * 8.0)
    end
    for i = 1:d
        @test_approx_eq_eps emp_md[i] md[i] (sqrt(s[i] / n_tsamples) * 8.0)
    end
    for i = 1:d, j = 1:d
        @test_approx_eq_eps emp_cov[i,j] S[i,j] (sqrt(s[i] * s[j]) * 20.0) / sqrt(n_tsamples)
    end

    # evaluation of logpdf and pdf
    for i = 1:min(100, n_tsamples)
        @test_approx_eq logpdf(g, X[:,i]) log(pdf(g, X[:,i]))
    end
    @test_approx_eq logpdf(g, X) @compat(log.(pdf(g, X)))
    @test isequal(logpdf(g, zeros(d)),-Inf)
    @test isequal(logpdf(g, -mn),-Inf)
    @test isequal(pdf(g, zeros(d)),0.0)
    @test isequal(pdf(g, -mn),0.0)

    # test the location and scale functions
    @test_approx_eq_eps location(g) location(MvLogNormal,:meancov,mean(g),cov(g)) 1e-8
    @test_approx_eq_eps location(g) location(MvLogNormal,:mean,mean(g),scale(g)) 1e-8
    @test_approx_eq_eps location(g) location(MvLogNormal,:median,median(g),scale(g)) 1e-8
    @test_approx_eq_eps location(g) location(MvLogNormal,:mode,mode(g),scale(g)) 1e-8
    @test_approx_eq_eps scale(g) scale(MvLogNormal,:meancov,mean(g),cov(g)) 1e-8

    @test_approx_eq_eps location(g) location!(MvLogNormal,:meancov,mean(g),cov(g),zeros(mn)) 1e-8
    @test_approx_eq_eps location(g) location!(MvLogNormal,:mean,mean(g),scale(g),zeros(mn)) 1e-8
    @test_approx_eq_eps location(g) location!(MvLogNormal,:median,median(g),scale(g),zeros(mn)) 1e-8
    @test_approx_eq_eps location(g) location!(MvLogNormal,:mode,mode(g),scale(g),zeros(mn)) 1e-8
    @test_approx_eq_eps scale(g) scale!(MvLogNormal,:meancov,mean(g),cov(g),zeros(S)) 1e-8

    lc1,sc1 = params(MvLogNormal,mean(g),cov(g))
    lc2,sc2 = params!(MvLogNormal,mean(g),cov(g),similar(mn),similar(S))
    @test_approx_eq_eps location(g) lc1  1e-8
    @test_approx_eq_eps location(g) lc2  1e-8
    @test_approx_eq_eps scale(g) sc1  1e-8
    @test_approx_eq_eps scale(g) sc2  1e-8
end

####### Validate results for a single-dimension MvLogNormal by comparing with univariate LogNormal
println("    comparing results from MvLogNormal with univariate LogNormal")
l1 = LogNormal(0.1,0.4)
l2 = MvLogNormal(0.1*ones(1),0.4)
@test_approx_eq [mean(l1)] mean(l2)
@test_approx_eq [median(l1)] median(l2)
@test_approx_eq [mode(l1)] mode(l2)
@test_approx_eq [var(l1)] var(l2)
@test_approx_eq [entropy(l1)] entropy(l2)
@test_approx_eq logpdf(l1,5.0) logpdf(l2,[5.0])
@test_approx_eq pdf(l1,5.0) pdf(l2,[5.0])
@test (srand(78393) ; [rand(l1)]) == (srand(78393) ; rand(l2))

###### General Testing

mu = [0.1, 0.2, 0.3]
va = [0.16, 0.25, 0.36]
C = [0.4 -0.2 -0.1; -0.2 0.5 -0.1; -0.1 -0.1 0.6]

for (g, μ, Σ) in [
    (MvLogNormal(mu,PDMats.PDMat(C)), mu, C),
    (MvLogNormal(PDMats.PDiagMat(Vector{Float64}(@compat(sqrt.(va))))), zeros(3), diagm(va)), # Julia 0.4 loses type information so Vector{Float64} can be dropped when we don't support 0.4
    (MvLogNormal(mu, sqrt(0.2)), mu, 0.2 * eye(3)),
    (MvLogNormal(3, sqrt(0.2)), zeros(3), 0.2 * eye(3)),
    (MvLogNormal(mu, Vector{Float64}(@compat(sqrt.(va)))), mu, diagm(va)), # Julia 0.4 loses type information so Vector{Float64} can be dropped when we don't support 0.4
    (MvLogNormal(Vector{Float64}(@compat(sqrt.(va)))), zeros(3), diagm(va)), # Julia 0.4 loses type information so Vector{Float64} can be dropped when we don't support 0.4
    (MvLogNormal(mu, C), mu, C),
    (MvLogNormal(C), zeros(3), C) ]

    println("    testing $(typeof(g)) with normal distribution $(Distributions.distrname(g.normal))")

    m,s = params(g)
    @test_approx_eq full(m) μ
    test_mvlognormal(g, 10^4)
end

##### Constructors and conversions
d = MvLogNormal(Array{Float32}(mu), PDMats.PDMat(Array{Float32}(C)))
@test typeof(convert(MvLogNormal{Float64}, d)) == typeof(MvLogNormal(mu, PDMats.PDMat(C)))
@test typeof(convert(MvLogNormal{Float64}, d.normal.μ, d.normal.Σ)) == typeof(MvLogNormal(mu, PDMats.PDMat(C)))


text:
"""
	grd2kml(cmd0::String="", arg1=[], kwargs...)

Reads a 2-D grid file and makes a quadtree of PNG images and KML wrappers for Google Earth
using the selected tile size [256x256 pixels].

Full option list at [`grd2kml`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html)

Parameters
----------

- $(GMT.opt_C)
- **E** : **url** : -- Str --		Flags = ``url``

    Instead of hosting the files locally, prepend a site URL. The top-level prefix.kml file
    will then use this URL to find the other files it references.``
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#e)
- **F** : **filter** : -- Str --

    Specifies the filter to use for the downsampling of the grid for more distant viewing.
    Choose among boxcar, cosine arch, gaussian, or median [Gaussian].
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#f)
- **H** : **sub_pixel** : -- Int --         Flags = ``factor``

    Improve the quality of rasterization by passing the sub-pixel smoothing factor to psconvert.
    [`-H`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#h)
- **I** : **shade** : **intensity** : **intensfile** : -- Str or GMTgrid --

    Gives the name of a grid file or GMTgrid with intensities in the (-1,+1) range,
    or a grdgradient shading flags.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#i)
- **L** : **tile_size** : -- Number --			Flags = tilesize

    Sets the fixed size of the image building blocks. Must be an integer that is radix 2.
    Typical values are 256 or 512 [256].
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#l)
- **N** : **prefix** -- Str --		            Flags = ``prefix``

    Sets a unique name prefixed used for the top-level KML filename and the directory where all
    referenced KML files and PNG images will be written [GMT_Quadtree].
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#n)
- **Q** : **nan_t** : **nan_alpha** : -- Bool or [] --

    Make grid nodes with z = NaN transparent, using the color-masking feature in PostScript Level 3.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#q)
- **T** : **title** : -- Str --			        Flags = ``title``

    Sets the title of the top-level document (i.e., its description).
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grd2kml.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_f)
"""
function grd2kml(cmd0::String="", arg1=[], arg2=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grd2kml", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_f(cmd, d)

	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', 0, arg1, arg2)

	for sym in [:I :shade :intensity :intensfile]
		if (haskey(d, sym))
			if (!isa(d[sym], GMTgrid))		# Uff, simple. Either a file name or a -A type modifier
				cmd = cmd * " -I" * arg2str(d[sym])
			else
				cmd,N_shade = put_in_slot(cmd, d[sym], 'I', [arg1, arg2])
				if (N_shade == 1)     arg1 = d[sym]
				elseif (N_shade == 2) arg2 = d[sym]
				end
			end
			break
		end
    end

	cmd = add_opt(cmd, 'E', d, [:E :url])
	cmd = add_opt(cmd, 'F', d, [:F :filter])
	cmd = add_opt(cmd, 'H', d, [:H :sub_pixel])
	cmd = add_opt(cmd, 'L', d, [:L :tile_size])
	cmd = add_opt(cmd, 'N', d, [:N :prefix])
	cmd = add_opt(cmd, 'Q', d, [:Q :nan_t :nan_alpha])
	cmd = add_opt(cmd, 'T', d, [:T :title])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end
	return common_grd(d, cmd, got_fname, 1, "grd2kml", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grd2kml(arg1=[], arg2=[], cmd0::String=""; kw...) = grd2kml(cmd0, arg1, arg2; kw...)


text:
struct ODELocalSensitivityFunction{iip,F,A,Tt,J,JP,PJ,TW,TWt,UF,PF,JC,PJC,Alg,fc,JM,pJM,MM,CV} <: DiffEqBase.AbstractODEFunction{iip}
  f::F
  analytic::A
  tgrad::Tt
  jac::J
  jac_prototype::JP
  paramjac::PJ
  Wfact::TW
  Wfact_t::TWt
  uf::UF
  pf::PF
  J::JM
  pJ::pJM
  jac_config::JC
  paramjac_config::PJC
  alg::Alg
  numparams::Int
  numindvar::Int
  f_cache::fc
  mass_matrix::MM
  isautojacvec::Bool
  colorvec::CV
end

function ODELocalSensitivityFunction(f,analytic,tgrad,jac,jac_prototype,paramjac,Wfact,Wfact_t,uf,pf,u0,
                                    jac_config,paramjac_config,alg,p,f_cache,mm,
                                    isautojacvec,colorvec)
  numparams = length(p)
  numindvar = length(u0)
  J = isautojacvec ? nothing : Matrix{eltype(u0)}(undef,numindvar,numindvar)
  pJ = Matrix{eltype(u0)}(undef,numindvar,numparams) # number of funcs size
  ODELocalSensitivityFunction{isinplace(f),typeof(f),typeof(analytic),
                             typeof(tgrad),typeof(jac),typeof(jac_prototype),typeof(paramjac),
                             typeof(Wfact),typeof(Wfact_t),typeof(uf),
                             typeof(pf),typeof(jac_config),
                             typeof(paramjac_config),typeof(alg),
                             typeof(f_cache),
                             typeof(J),typeof(pJ),typeof(mm),typeof(f.colorvec)}(
                             f,analytic,tgrad,jac,jac_prototype,paramjac,Wfact,Wfact_t,uf,pf,J,pJ,
                             jac_config,paramjac_config,alg,
                             numparams,numindvar,f_cache,mm,isautojacvec,colorvec)
end

function (S::ODELocalSensitivityFunction)(du,u,p,t)
  y = @view u[1:S.numindvar] # These are the independent variables
  dy = @view du[1:S.numindvar]
  S.f(dy,y,p,t) # Make the first part be the ODE

  # Now do sensitivities
  # Compute the Jacobian

  if !S.isautojacvec
    if DiffEqBase.has_jac(S.f)
      S.jac(S.J,y,p,t) # Calculate the Jacobian into J
    else
      S.uf.t = t
      jacobian!(S.J, S.uf, y, S.f_cache, S.alg, S.jac_config)
    end
  end

  if DiffEqBase.has_paramjac(S.f)
    S.paramjac(S.pJ,y,p,t) # Calculate the parameter Jacobian into pJ
  else
    S.pf.t = t
    S.pf.u .= y
    jacobian!(S.pJ, S.pf, p, S.f_cache, S.alg, S.paramjac_config)
  end

  # Compute the parameter derivatives
  for i in eachindex(p)
    Sj = @view u[i*S.numindvar+1:(i+1)*S.numindvar]
    dp = @view du[i*S.numindvar+1:(i+1)*S.numindvar]
    if !S.isautojacvec
      mul!(dp,S.J,Sj)
    else
      jacobianvec!(dp, S.uf, y, Sj, S.alg, S.jac_config)
    end
    dp .+= @view S.pJ[:,i]
  end
end

function ODELocalSensitivityProblem(f::DiffEqBase.AbstractODEFunction,u0,
                                    tspan,p=nothing,
                                    alg = SensitivityAlg(autojacvec=true);
                                    kwargs...)
  isinplace = DiffEqBase.isinplace(f)
  # if there is an analytical Jacobian provided, we are not going to do automatic `jac*vec`
  isautojacvec = get_jacvec(alg)
  p == nothing && error("You must have parameters to use parameter sensitivity calculations!")
  uf = DiffEqDiffTools.UJacobianWrapper(f,tspan[1],p)
  pf = DiffEqDiffTools.ParamJacobianWrapper(f,tspan[1],copy(u0))
  if isautojacvec
    if alg_autodiff(alg)
      # if we are using automatic `jac*vec`, then we need to use a `jac_config`
      # that is a tuple in the form of `(seed, buffer)`
      jac_config_seed = ForwardDiff.Dual{:___jac_tag}.(u0,u0)
      jac_config_buffer = similar(jac_config_seed)
      jac_config = jac_config_seed, jac_config_buffer
    else
      jac_config = (similar(u0),similar(u0))
    end
  elseif DiffEqBase.has_jac(f)
    jac_config = nothing
  else
    jac_config = build_jac_config(alg,uf,u0)
  end

  if DiffEqBase.has_paramjac(f)
    paramjac_config = nothing
  else
    paramjac_config = build_param_jac_config(alg,pf,u0,p)
  end

  # TODO: Use user tgrad. iW can be safely ignored here.
  sense = ODELocalSensitivityFunction(f,f.analytic,nothing,f.jac,f.jac_prototype,f.paramjac,nothing,nothing,
                                     uf,pf,u0,jac_config,
                                     paramjac_config,alg,
                                     p,similar(u0),f.mass_matrix,
                                     isautojacvec,f.colorvec)
  sense_u0 = [u0;zeros(sense.numindvar*sense.numparams)]
  ODEProblem(sense,sense_u0,tspan,p;kwargs...)
end

function ODELocalSensitivityProblem(f,args...;kwargs...)
  ODELocalSensitivityProblem(ODEFunction(f),args...;kwargs...)
end

# Get ODE u vector and sensitivity values from all time points
function extract_local_sensitivities(sol)
  ni = sol.prob.f.numindvar
  u = sol[1:ni, :]
  du = [sol[ni*j+1:ni*(j+1),:] for j in 1:sol.prob.f.numparams]
  return u, du
end

extract_local_sensitivities(sol, i::Integer, asmatrix::Val=Val(false)) = _extract(sol, sol[i], asmatrix)
extract_local_sensitivities(sol, i::Integer, asmatrix::Bool) = extract_local_sensitivities(sol, i, Val{asmatrix}())
extract_local_sensitivities(sol, t, asmatrix::Val=Val(false)) = _extract(sol, sol(t), asmatrix)
extract_local_sensitivities(sol, t, asmatrix::Bool) = extract_local_sensitivities(sol, t, Val{asmatrix}())
extract_local_sensitivities(tmp, sol, t, asmatrix::Val=Val(false)) = _extract(sol, sol(tmp, t), asmatrix)
extract_local_sensitivities(tmp, sol, t, asmatrix::Bool) = extract_local_sensitivities(tmp, sol, t, Val{asmatrix}())

# Get ODE u vector and sensitivity values from sensitivity problem u vector
function _extract(sol, su::Vector, asmatrix::Val = Val(false))
  u = view(su, 1:sol.prob.f.numindvar)
  du = _extract_du(sol, su, asmatrix)
  return u, du
end

# Get sensitivity values from sensitivity problem u vector (nested form)
function _extract_du(sol, su::Vector, ::Val{false})
  ni = sol.prob.f.numindvar
  return [view(su, ni*j+1:ni*(j+1)) for j in 1:sol.prob.f.numparams]
end

# Get sensitivity values from sensitivity problem u vector (matrix form)
function _extract_du(sol, su::Vector, ::Val{true})
  ni = sol.prob.f.numindvar
  np = sol.prob.f.numparams
  return view(reshape(su, ni, np+1), :, 2:np+1)
end


text:
@testset "APE" begin

    @testset "Conversions" begin

        @testset "From UInt8 to DNAbin" begin
            @test DNAbin(0x88) === DNAbin_A
            @test DNAbin(0x18) === DNAbin_T
            @test DNAbin(0x28) === DNAbin_C
            @test DNAbin(0x48) === DNAbin_G
            @test DNAbin(0xC0) === DNAbin_R
            @test DNAbin(0x30) === DNAbin_Y
            @test DNAbin(0x60) === DNAbin_S
            @test DNAbin(0x90) === DNAbin_W
            @test DNAbin(0x50) === DNAbin_K
            @test DNAbin(0xA0) === DNAbin_M
            @test DNAbin(0x70) === DNAbin_B
            @test DNAbin(0xD0) === DNAbin_D
            @test DNAbin(0xB0) === DNAbin_H
            @test DNAbin(0xE0) === DNAbin_V
            @test DNAbin(0xF0) === DNAbin_N
            @test DNAbin(0x04) === DNAbin_Gap
        end

        @testset "From DNA to DNAbin" begin
            @test DNAbin(DNA_A) === DNAbin_A
            @test DNAbin(DNA_T) === DNAbin_T
            @test DNAbin(DNA_C) === DNAbin_C
            @test DNAbin(DNA_G) === DNAbin_G
            @test DNAbin(DNA_R) === DNAbin_R
            @test DNAbin(DNA_Y) === DNAbin_Y
            @test DNAbin(DNA_S) === DNAbin_S
            @test DNAbin(DNA_W) === DNAbin_W
            @test DNAbin(DNA_K) === DNAbin_K
            @test DNAbin(DNA_M) === DNAbin_M
            @test DNAbin(DNA_B) === DNAbin_B
            @test DNAbin(DNA_D) === DNAbin_D
            @test DNAbin(DNA_H) === DNAbin_H
            @test DNAbin(DNA_V) === DNAbin_V
            @test DNAbin(DNA_N) === DNAbin_N
            @test DNAbin(DNA_Gap) === DNAbin_Gap
        end

        @testset "From DNAbin to DNA" begin
            @test DNA(DNAbin_A) === DNA_A
            @test DNA(DNAbin_T) === DNA_T
            @test DNA(DNAbin_C) === DNA_C
            @test DNA(DNAbin_G) === DNA_G
            @test DNA(DNAbin_R) === DNA_R
            @test DNA(DNAbin_Y) === DNA_Y
            @test DNA(DNAbin_S) === DNA_S
            @test DNA(DNAbin_W) === DNA_W
            @test DNA(DNAbin_K) === DNA_K
            @test DNA(DNAbin_M) === DNA_M
            @test DNA(DNAbin_B) === DNA_B
            @test DNA(DNAbin_D) === DNA_D
            @test DNA(DNAbin_H) === DNA_H
            @test DNA(DNAbin_V) === DNA_V
            @test DNA(DNAbin_N) === DNA_N
            @test DNA(DNAbin_Gap) === DNA_Gap
        end

        @testset "Char" begin
            @testset "DNAbin conversions from Char" begin
                @test convert(DNAbin, 'A') === DNAbin_A
                @test convert(DNAbin, 'C') === DNAbin_C
                @test convert(DNAbin, 'G') === DNAbin_G
                @test convert(DNAbin, 'T') === DNAbin_T
                @test convert(DNAbin, 'N') === DNAbin_N
                @test_throws InexactError convert(DNAbin, 'Z')
                @test_throws InexactError convert(DNAbin, '核')
            end

            @testset "DNAbin conversions to Char" begin
                @test convert(Char, DNAbin_A) == 'A'
                @test convert(Char, DNAbin_C) == 'C'
                @test convert(Char, DNAbin_G) == 'G'
                @test convert(Char, DNAbin_T) == 'T'
                @test convert(Char, DNAbin_N) == 'N'
            end
        end

    end

    @testset "iscompatible" begin
        @test  iscompatible(DNAbin_A, DNAbin_A)
        @test  iscompatible(DNAbin_A, DNAbin_R)
        @test !iscompatible(DNAbin_C, DNAbin_A)
        @test !iscompatible(DNAbin_C, DNAbin_R)

        for x ∈ alphabet(DNAbin)
            @test iscompatible(x, DNAbin_N) == (x != DNAbin_Gap)
            @test iscompatible(DNAbin_N, x) == (x != DNAbin_Gap)
        end
    end

    @testset "isambiguous" begin
        for nt ∈ alphabet(DNAbin)
            @test isambiguous(nt) == (nt ∉ (DNAbin_A, DNAbin_C, DNAbin_G, DNAbin_T, DNAbin_Gap))
        end
    end

    @testset "isGC" begin
        for nt ∈ alphabet(DNAbin)
            @test isGC(nt) == (nt ∈ (DNAbin_G, DNAbin_C, DNAbin_S))
        end
    end

    @testset "ispurine" begin
        for nt ∈ alphabet(DNAbin)
            @test ispurine(nt) == (nt == DNAbin_A || nt == DNAbin_G || nt == DNAbin_R)
        end
    end

    @testset "ispyrimidine" begin
        for nt ∈ alphabet(DNAbin)
            @test ispyrimidine(nt) == (nt == DNAbin_T || nt == DNAbin_C || nt == DNAbin_Y)
        end
    end

    @testset "iscertain" begin
        for nt ∈ alphabet(DNAbin)
            @test iscertain(nt) == (nt ∈ (DNAbin_A, DNAbin_G, DNAbin_T, DNAbin_C,))
        end
    end

    @testset "isgap" begin
        for nt ∈ alphabet(DNAbin)
            @test isgap(nt) == (nt === gap(DNAbin))
        end
    end

    @testset "complement" begin
        @test complement(DNAbin_A) === DNAbin_T
        @test complement(DNAbin_C) === DNAbin_G
        @test complement(DNAbin_G) === DNAbin_C
        @test complement(DNAbin_T) === DNAbin_A
        @test complement(DNAbin_Gap) === DNAbin_Gap
        @test complement(DNAbin_N) === DNAbin_N
    end

    @testset "show" begin
        buf = IOBuffer()
        for nt ∈ [DNAbin_A, DNAbin_C, DNAbin_G, DNAbin_T, DNAbin_N, DNAbin_Gap]
            show(buf, nt)
            write(buf, ' ')
        end
        @test String(take!(buf)) == "DNAbin_A DNAbin_C DNAbin_G DNAbin_T DNAbin_N DNAbin_Gap "
    end

    @testset "Data transfer" begin
        rseq = R"""
        library(ape)
        a <- as.DNAbin(
            matrix(
                c('a', 't', 'c', 'g', 'a', 't', 'c', 'g', '-', 't', 'c', 'm'),
                nrow = 3
            )
        )
        """
        dba = [DNAbin_A  DNAbin_G  DNAbin_C    DNAbin_T
               DNAbin_T  DNAbin_A  DNAbin_G    DNAbin_C
               DNAbin_C  DNAbin_T  DNAbin_Gap  DNAbin_M]
        bjdba = [DNA_A  DNA_G  DNA_C    DNA_T
                 DNA_T  DNA_A  DNA_G    DNA_C
                 DNA_C  DNA_T  DNA_Gap  DNA_M]
        vs = DNASequence[dna"AGCT", dna"TAGC", dna"CT-M"]
        @test rcopy(Array{DNAbin, 2}, rseq) == dba
        @test rcopy(Vector{DNASequence}, rseq) == vs
        @test rcopy(Array{DNA, 2}, rseq) == bjdba
        @test all(i == j for (i, j) ∈ zip(rseq, RObject(RCall.sexp(RCall.RawSxp, dba))))
        @test all(i == j for (i, j) ∈ zip(rseq, RObject(RCall.sexp(RCall.RawSxp, vs))))
    end

end


text:
using DifferentialEquations, PGFPlotsX, LaTeXStrings, Colors
export lif_pair_get_solution, lnp_pair_get_solution, setup_pgfplotsx

setup_pgfplotsx() = begin
push!(PGFPlotsX.CUSTOM_PREAMBLE,
raw"""
\usepackage{amsmath}
\usepgfplotslibrary{colorbrewer}
\usepgfplotslibrary{fillbetween}
\pgfplotsset{
    cycle list/Set1,
}
"""
)
end

setup_pgfplotsx()

example_trange=(0,10.0)
example_yrange=(0,2.0)
example_input=t-> 0.5*(sin(t)+1) + 0.3*(cos(0.7+1.5*t)+1)
example_input_integral=t-> -0.5*cos(t) + 0.3*sin(0.7+1.5*t)/1.5 + 0.8*t

discretize(x, levels, values) = values[searchsortedlast(levels, x)]

function lif_pair_f′(du, u, p, t)
    du[1] = p.α*p.input(t) - p.α*u[1]
    du[2] = -p.β * u[2]
    du[3] = p.β*(p.input(t) - u[3])
end

lif_pair_spike_condition(u,t,integrator) = u[1]-integrator.p.δ

function lif_pair_spike!(integrator)
    integrator.u[1] = 0.0
    integrator.u[2] += integrator.p.δ
    push!(integrator.p.spiketimes, integrator.t)
end

function lif_pair_get_solution(tspan, input;  α = 0.5, β = 0.5, δ = 0.2, u₀ = [0.0,1.0,1.0])
    spiketimes=Float64[]
    p = (;input,α,β,δ,spiketimes)
    cb = ContinuousCallback(lif_pair_spike_condition, lif_pair_spike!, nothing)
    sol = solve(ODEProblem(lif_pair_f′, u₀, tspan, p); callback=cb)
    return sol, spiketimes
end

function lnp_pair_f′(du, u, p, t)
    du[1] = -p.β * u[1]
    du[2] = p.β*(p.input(t) - u[2])
end

function lnp_pair_spike!(integrator)
    integrator.u[1] += integrator.p.δ
    push!(integrator.p.spiketimes, integrator.t)
end

function lnp_pair_get_solution(tspan, input;  β = 0.5, δ=0.2, u₀ = [1.0,1.0], kwargs...)
    spiketimes = Float64[]
    p = (;input,β, δ, spiketimes)
    rate_wrapped(u,p,t) = input(t)*β/δ
    prob = ODEProblem(lnp_pair_f′, u₀, tspan, p)
    jump = VariableRateJump(rate_wrapped, lnp_pair_spike!)
    jump_prob=JumpProblem(prob, Direct(),jump)
    sol = solve(jump_prob,Tsit5(); kwargs...)

    return sol, spiketimes
end

max_stars_repo_name: jleugeri/phd

max_issues_repo_name: jleugeri/phd

max_forks_repo_name: jleugeri/phd


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module NHosts1Vector
# # Model Specifications

# ## Host Specifications
#
# The host may be in one of 4 states:
#
# - susceptible (S_H)
#
# - incubation period (E_H)
#
# - infectious (I_H)
#
# - recovered (R_H)
#
# ## Vector Specifications
#
# The vector may be in one of 3 states:
#
# - susceptible (S_V)
#
# - incubation period (E_V)
#
# - infectious (I_V)

# # Model Parameters
#
# σ_H, σ_V: the incubation rates for hosts & vectors (units: per time)
#
# μ_H, μ_V: the mortality rates for hosts & vectors (units: per time)
#
# λ: the clearance (or recovery) rate for hosts (units: per time)
#
# β: the infection rate (units: per capita per time)

# $$N_H^i = S_H^i + E_H^i + I_H^i + R_H^i, \dot{N_H^i}=0$$
#
# $$N_H = \sum_{i=1}^N N_H^i$$
#
# $$\dot{S_H}^i = \mu_H^i(E_H^i + I_H^i + R_H^i) - \frac{\beta^i}{N_H} S_H^i I_V$$
#
# $$\dot{E_H}^i = \frac{\beta^i}{N_H} S_H^i I_V - (\sigma_H^i + \mu_H^i) E_H^i$$
#
# $$\dot{I_H}^i = \sigma_H^i E_H^i - (\lambda^i + \mu_H^i) I_H^i$$
#
# $$\dot{R_H}^i = \lambda^i I_H^i - \mu_H^i R_H^i$$
#
# $$\dot{S_V} = \mu_V(E_V + I_V) - \frac{\sum_{i=1}^N\beta^i I_H^i}{N_H} S_V$$
#
# $$\dot{E_V} = \frac{\sum_{i=1}^N\beta^i I_H^i}{N_H} S_V - (\sigma_V + \mu_V) E_V$$
#
# $$\dot{I_V} = \sigma_V E_V - \mu_V I_V$$

# This version considers multiple hosts, which introduces some indexing complications. This approach adopts the indexing scheme of:
#
# - all of host $i$ compartments ($S_H$, $E_H$, etc), for $i∈1…N$
# - vector compartments last
# With this approach, we can re-use the solution for one host, with slight modifications:

# +
H_comps = 4
V_comps = 3

# in sub functions, du / u are the particular relevant slices only
function F1H(du,u,p,t,I_V,N_H)
    S_H, E_H, I_H, R_H = u

    # host dynamics
    host_infection = (p.β*S_H*I_V)/N_H
    host_mortality = p.μ_H .* u # include S_H, so easier to remove mortality
    host_births = sum(host_mortality)
    host_progression = p.σ_H*E_H
    recovery = p.λ*I_H

    du[1] = -host_infection + host_births
    du[2] = host_infection - host_progression
    du[3] = host_progression - recovery
    du[4] = recovery
    du[1:end] -= host_mortality
end

# in sub functions, du / u are the particular relevant slices only
function FV(du,u,p,t,sum_β_I_H)
    S_V, E_V, I_V = u
    vec_infection = sum_β_I_H*S_V/p.N_H
    vec_mortality = p.μ_V .* u # include S_V, so easier to remove mortality
    vec_births = sum(vec_mortality)
    vec_progression = p.σ_V*E_V

    du[1] = -vec_infection + vec_births
    du[2] = vec_infection - vec_progression
    du[3] = vec_progression
    du[1:end] -= vec_mortality
end

function F(du,u,p,t)
    uvec = @view(u[(p.nHosts*H_comps+1):end]) # grab the vector compartments
    S_V, E_V, I_V = uvec
    sum_β_I_H = 0.0
    for host in 0:(p.nHosts-1)
        slice = (1:H_comps).+(H_comps*host)
        F1H(@view(du[slice]), @view(u[slice]), p.host[host+1], t, I_V, p.vec.N_H)
        # must use @view here, so that these arrays can be modified in F1H
        sum_β_I_H += p.host[host+1].β * u[slice[3]] # this host's I compartment
    end
    FV(@view(du[(p.nHosts*4+1):end]), uvec, p.vec,t,sum_β_I_H)
end
# -

# First, state initial conditions. This code generates them randomly for convenience, though they could be assigned based on data, desired parameter space, or algorithmically as part of a fitting process:

# +
using Random
nH = 5
Random.seed!(0)

S_Hs = ones(nH) .* 100.0
E_Hs = zeros(nH)
I_Hs = shuffle(vcat(zeros(nH-1),[1.0]))
R_Hs = zeros(nH)
host0 = reshape(hcat(S_Hs,E_Hs,I_Hs,R_Hs)', nH*4, 1)
vec0 = [10000.0, 0.0, 0.0]
u0 = vcat(host0, vec0);
# -

# Now, generate dynamic parameters. Again: this code generates them randomly for convenience, though they could be assigned based on data, desired parameter space, or algorithmically as part of a fitting process:

# +
Random.seed!(1)

μs = 1 ./ (rand(nH) .* 360)
σs = 1 ./ (rand(nH) .* 6)
λs = 1 ./ (rand(nH) .* 28)
βs = rand(nH) ./ 10.0

p = (
  nHosts = nH,
  vec = (μ_V=1/30, σ_V=1/7, N_H = sum(host0)),
  host = [(μ_H=μs[i], σ_H=σs[i], λ=λs[i], β=βs[i]) for i in 1:nH]
  # just building up a random collection of params for demonstration
);

# +
using DifferentialEquations
using IterableTables, DataFrames

tspan = (0.0, 365.0)
prob = ODEProblem(F, u0, tspan, p)
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8,saveat=collect(range(0,stop = 365,length = 365*10+1)));
# -

df = DataFrame(sol)
namespace = zip(names(df)[2:end],[
    :S_H1,:E_H1,:I_H1,:R_H1,
    :S_H2,:E_H2,:I_H2,:R_H2,
    :S_H3,:E_H3,:I_H3,:R_H3,
    :S_H4,:E_H4,:I_H4,:R_H4,
    :S_H5,:E_H5,:I_H5,:R_H5,
    :S_V,:E_V,:I_V])
namespace = collect(namespace)
rename!(df,namespace);

using Plots
host_1 = plot(df[:timestamp],[df[:,:S_H1],df[:,:E_H1],df[:,:I_H1],df[:,:R_H1]])
host_2 = plot(df[:timestamp],[df[:,:S_H2],df[:,:E_H2],df[:,:I_H2],df[:,:R_H2]])
host_3 = plot(df[:timestamp],[df[:,:S_H3],df[:,:E_H3],df[:,:I_H3],df[:,:R_H3]])
host_4 = plot(df[:timestamp],[df[:,:S_H4],df[:,:E_H4],df[:,:I_H4],df[:,:R_H4]])
plot(host_1,host_2,host_3,host_4,layout=(4,1),label=["S","E","I","R"])

end


text:
# Load modules
using Plots
using NCDatasets
using Statistics: mean

# Setup run-time enviromnent
ENV["GKSwstype"] = "100"

# Specify needed directories and filenames: these will be replaced automatically by the bash run script or can be modified manually
CLIMA_NETCDF = "../netcdf/";

# Get the current and previous  GCM netcdf file names in the CLIMA_NETCDF directory
fnames = filter(x -> occursin(".nc", x), readdir( CLIMA_NETCDF ) );
filename = "$CLIMA_NETCDF"fnames[end] # selects the last file on list

# extract data
ds = NCDataset(filename, "r");
lon = ds["long"][:];
lat = ds["lat"][:];
lev = ds["level"][:] / 1e3; # height in kilometers
time = ds["time"][:];
u = ds["u"][:];
v = ds["v"][:];
T = ds["temp"][:];
close(ds)

# zonal mean:T and u at last diagnostic time
time_index = length(time)
u_zm = mean( u[:,:,:,:], dims=1)[1,:,:,:]; # lon, lat,lev, time
T_zm = mean( T[:,:,:,:], dims=1)[1,:,:,:]; # lon, lat,lev, time

plot1 = contourf( lat, lev, (u_zm[:,:,time_index])', title="u", xlabel="lat (deg N)", ylabel="z (km)");
plot_array = [plot1]
plot2 = contourf( lat, lev, (T_zm[:,:,time_index])', title="T", xlabel="lat (deg N)", ylabel="z (km)");
push!(plot_array,plot2);
fig=plot(plot_array... , layout=(1, 2), size=(800, 400) )
savefig(fig, "zontal_mean.pdf");
display(fig)

# vertical slice of v at lev_index
lev_index = 10
z_in_km = lev[lev_index]
v_plot = contourf( lon, lat, (v[:,:,lev_index,time_index])', title="v @ $z_in_km", xlabel="lon (deg)", ylabel="lat (deg N)");
fig=plot(v_plot, layout=(1, 1), size=(800, 400) )
savefig(fig, "vertical_slice.pdf");
display(fig)


text:
using Unitful, UnitfulOfu, Test

import Unitful: ft, inch, m, g, kg, °C, °, W, kW, oz, Hz, rad, s, N, uconvert
import UnitfulOfu: rpm, lbf, klbf, hp

@testset "OFU conversions" begin
  @test uconvert(N, 1.0lbf) == 4.45N
  @test 1000.0lbf == 1.0klbf
  @test 1.0hp == 735.4990027920646W
  @test 5.0rpm - 1//60*1.0Hz == 0.06666666666666667s^-1
  @test uconvert(Hz, 61.0rpm - 1.0Hz) == 0.016666666666666607Hz
end

#=
@testset "OFU string macro" begin
    @test macroexpand(:(ofu"lbf")) == :(UnitfulOfu.ustrcheck(UnitfulOfu.lbf))
    @test macroexpand(:(ofu"1.0")) == 1.0
    @test macroexpand(:(ofu"rpm")) == :(UnitfulOfu.ustrcheck(UnitfulOfu.rpm))
    @test macroexpand(:(ofu"1.0hp")) == :(1.0 * UnitfulOfu.ustrcheck(UnitfulOfu.hp))
end
=#


text:
#
# Runtime representation of input and output values.
#

#
# Structured input or output type.
#

abstract Kind{T}

eltype{T}(::Type{Kind{T}}) = T

# Internal representation of the structure.
data(X::Kind) = X.data

#
# Structure-free input or output.
#

immutable Iso{T} <: Kind{T}
    data::T

    Iso(val::T) = new(val)
end

Iso{T}(val::T) = Iso{T}(val)

data{T}(::Type{Iso{T}}) = T

Kind{K<:Iso}(::Type{K}, T::Type) = Iso{T}

start(X::Iso) = true
next(X::Iso, state) = (data(X), false)
done(X::Iso, state) = state

#
# Partial output.
#

typealias OptData{T} Nullable{T}

immutable Opt{T} <: Kind{T}
    data::OptData{T}

    Opt(data::OptData{T}) = new(data)
    Opt(val::T) = new(Nullable{T}(val))
    Opt() = new(OptData{T}())
end

Opt{T}(val::T) = Opt{T}(Nullable{T}(val))
Opt() = Opt{Union{}}(Nullable())

isnull(X::Opt) = isnull(X.data)
get(X::Opt) = get(X.data)
get(X::Opt, default) = get(X.data, default)

data{T}(::Type{Opt{T}}) = OptData{T}

Kind{K<:Opt}(::Type{K}, T::Type) = Opt{T}

start(X::Opt) = !isnull(X)
next(X::Opt, state) = (get(X), false)
done(X::Opt, state) = state

#
# Plural output.
#

typealias SeqData{T} Vector{T}

immutable Seq{T} <: Kind{T}
    data::SeqData{T}

    Seq(data::SeqData{T}) = new(data)
    Seq() = new(T[])
end

Seq{T}(vals::AbstractVector{T}) = Seq{T}(vals)

length(X::Seq) = length(X.data)
isempty(X::Seq) = isempty(X.data)
endof(X::Seq) = endof(X.data)
eachindex(X::Seq) = eachindex(X.data)
getindex(X::Seq, keys...) = getindex(X.data, keys...)

data{T}(::Type{Seq{T}}) = SeqData{T}

Kind{K<:Seq}(::Type{K}, T::Type) = Seq{T}

start(X::Seq) = start(data(X))
next(X::Seq, state) = next(data(X), state)
done(X::Seq, state) = done(data(X), state)

#
# Input with surrounding context (that is, preceeding and subsequent values).
#

typealias RelData{T} Pair{Int,Vector{T}}

immutable Rel{T} <: Kind{T}
    data::RelData{T}

    Rel(data::RelData{T}) = new(data)
    Rel(ptr::Int, vals::Vector{T}) = new(RelData{T}(ptr, vals))
end

Rel{T}(ptr::Int, vals::Vector{T}) = Rel{T}(RelData(ptr, vals))

data{T}(::Type{Rel{T}}) = RelData{T}

Kind{K<:Rel}(::Type{K}, T::Type) = Rel{T}

#
# Input with named parameters.
#

# Ns is Val{(:x, :y, :z, ...)}, Vs is Tuple{Int, Int, Float64, ...}.
typealias EnvData{Ns,Vs,T} Pair{Vs,T}

immutable Env{Ns,Vs,T} <: Kind{T}
    data::EnvData{Ns,Vs,T}

    Env(data::EnvData{Ns,Vs,T}) = new(data)
    Env(env::Vs, val::T) = new(EnvData{Ns,Vs,T}(env, val))
end

Env{Ns<:Val,Vs<:Tuple,T}(::Type{Ns}, env::Vs, val::T) =
    Env{Ns,Vs,T}(EnvData(env, val))

data{Ns,Vs,T}(::Type{Env{Ns,Vs,T}}) = EnvData{Ns,Vs,T}

envsig{Ns,Vs,T}(::Type{Env{Ns,Vs,T}}) = Ns
envsig{Ns,Vs,T}(::Env{Ns,Vs,T}) = Ns
envtype{Ns,Vs,T}(::Type{Env{Ns,Vs,T}}) = Vs
envtype{Ns,Vs,T}(::Env{Ns,Vs,T}) = Vs

Kind{K<:Env}(::Type{K}, T::Type) = Env{envsig(K), envtype(K), T}

#
# Input with both named parameters and surrounding context.
#

typealias EnvRelData{Ns,Vs,T} EnvData{Ns,Vs,RelData{T}}

immutable EnvRel{Ns,Vs,T} <: Kind{T}
    data::EnvRelData{Ns,Vs,T}

    EnvRel(data::EnvRelData{Ns,Vs,T}) = new(data)
    EnvRel(env::Vs, ptr::Int, vals::Vector{T}) =
        new(EnvRelData{Ns,Vs,T}(env, RelData{T}(ptr, vals)))
end

EnvRel{Ns<:Val,Vs<:Tuple,T}(::Type{Ns}, env::Vs, ptr::Int, vals::Vector{T}) =
    EnvRel{Ns,Vs,T}(EnvData(env, RelData(ptr, vals)))

data{Ns,Vs,T}(::Type{EnvRel{Ns,Vs,T}}) = EnvRelData{Ns,Vs,T}

envsig{Ns,Vs,T}(::Type{EnvRel{Ns,Vs,T}}) = Ns
envsig{Ns,Vs,T}(::EnvRel{Ns,Vs,T}) = Ns
envtype{Ns,Vs,T}(::Type{EnvRel{Ns,Vs,T}}) = Vs
envtype{Ns,Vs,T}(::EnvRel{Ns,Vs,T}) = Vs

Kind{K<:EnvRel}(::Type{K}, T::Type) = EnvRel{envsig(K), envtype(K), T}

#
# Wrapping plain Julia values.
#

# Output.

convert{T}(K::Type{Iso{T}}, X::Iso{T}) = X
convert{T}(K::Type{Iso{T}}, val) =
    K(convert(T, val))

convert{T}(K::Type{Opt{T}}, X::Opt{T}) = X
convert{T}(K::Type{Opt{T}}, val) =
    K(convert(Nullable{T}, val))
convert{T}(K::Type{Opt{T}}) =
    K(Nullable{T}())

convert{T}(K::Type{Seq{T}}, X::Seq{T}) = X
convert{T}(K::Type{Seq{T}}, vals) =
    K(convert(Vector{T}, vals))

# Input.

convert{T}(K::Type{Iso{T}}) =
    K(nothing)

convert{T}(K::Type{Rel{T}}, X::Rel{T}) = X
convert{T}(K::Type{Rel{T}}, data::RelData{T}) =
    K(data)
convert{T}(K::Type{Rel{T}}, ptr::Int, vals) =
    K(RelData{T}(ptr, vals))
convert{T}(K::Type{Rel{T}}, val) =
    K(RelData{T}(1, T[val]))
convert{T}(K::Type{Rel{T}}) =
    K(RelData{T}(1, T[nothing]))

convert{Ns,Vs,T}(K::Type{Env{Ns,Vs,T}}, X::Env{Ns,Vs,T}) = X
convert{Ns,Vs,T}(K::Type{Env{Ns,Vs,T}}, data::EnvData{Ns,Vs,T}) =
    K(data)
convert{Ns,Vs,T}(K::Type{Env{Ns,Vs,T}}, env::Vs, val=nothing) =
    K(EnvData{Ns,Vs,T}(env, val))
convert{Ns,Vs,T}(K::Type{Env{Ns,Vs,T}}, params::Dict{Symbol,Any}=Dict{Symbol,Any}(), val=nothing) =
    let env = convertparams(Ns, Vs, params)
        K(env, val)
    end
convert{Ns,Vs,T}(K::Type{Env{Ns,Vs,T}}, args...; params...) =
    K(Dict{Symbol,Any}(params), args...)

convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, X::EnvRel{Ns,Vs,T}) = X
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, data::EnvRelData{Ns,Vs,T}) =
    K(data)
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, env::Vs, ptr::Int, vals) =
    K(EnvRelData{Ns,Vs,T}(env, RelData{T}(ptr, vals)))
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, env::Vs, val=nothing) =
    K(env, 1, T[val])
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, params::Dict{Symbol,Any}, ptr::Int, vals) =
    let env = convertparams(Ns,Vs, params)
        K(env, ptr, vals)
    end
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, ptr::Int, vals) =
    K(Dict{Symbol,Any}(), ptr, vals)
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, params::Dict{Symbol,Any}=Dict{Symbol,Any}(), val=nothing) =
    let env = convertparams(Ns, Vs, params)
        K(env, 1, T[val])
    end
convert{Ns,Vs,T}(K::Type{EnvRel{Ns,Vs,T}}, args...; params...) =
    K(Dict{Symbol,Any}(params), args...)

convertparams{names,Vs<:Tuple}(::Type{Val{names}}, ::Type{Vs}, params::Dict{Symbol,Any}) =
        ntuple(n -> convert(Vs.parameters[n], get(params, names[n], nothing)), length(names))

#
# Getting the value out of the input structure.
#

get(X::Iso) = data(X)

get(X::Rel) =
    begin
        ptr, vals = data(X)
        return vals[ptr]
    end

get(X::Env) =
    begin
        env, val = data(X)
        return val
    end

get(X::EnvRel) =
    begin
        env, (ptr, vals) = data(X)
        return vals[ptr]
    end

#
# Converting an input/output structure to a smaller/larger structure.
#

lift{K<:Kind}(::Type{K}, X::K) = X

# Output.

lift{T}(K::Type{Iso{T}}, X::Iso) = K(data(X))

lift{T}(K::Type{Opt{T}}, X::Iso) = K(Nullable{T}(data(X)))
lift{T}(K::Type{Opt{T}}, X::Opt) = K(Nullable{T}(data(X)))

lift{T}(K::Type{Seq{T}}, X::Iso) = K(T[data(X)])
lift{T}(K::Type{Seq{T}}, X::Opt) = K(isnull(X) ? T[] : T[get(X)])
lift{T}(K::Type{Seq{T}}, X::Seq) = K(Vector{T}(data(X)))

# Input.

lift{K<:Iso}(::Type{K}, X::Rel) =
    begin
        ptr, vals = data(X)
        return K(vals[ptr])
    end
lift{K<:Iso}(::Type{K}, X::Env) =
    begin
        env, val = data(X)
        return K(val)
    end
lift{K<:Iso}(::Type{K}, X::EnvRel) =
    begin
        env, (ptr, vals) = data(X)
        return K(vals[ptr])
    end

lift{K<:Rel}(::Type{K}, X::Rel) =
    begin
        ptr, vals = data(X)
        return K(ptr, vals)
    end
lift{K<:Rel}(::Type{K}, X::EnvRel) =
    begin
        env, (ptr, vals) = data(X)
        return K(ptr, vals)
    end

lift{K<:Env}(::Type{K}, X::Env) =
    begin
        Ns, Vs = envsig(K), envtype(K)
        Ns0, Vs0 = envsig(X), envtype(X)
        env, val = data(X)
        return K(liftenv(Ns, Vs, Ns0, Vs0, env), val)
    end
lift{K<:Env}(::Type{K}, X::EnvRel) =
    begin
        Ns, Vs = envsig(K), envtype(K)
        Ns0, Vs0 = envsig(X), envtype(X)
        env, (ptr, vals) = data(X)
        return K(liftenv(Ns, Vs, Ns0, Vs0, env), vals[ptr])
    end

lift{K<:EnvRel}(::Type{K}, X::EnvRel) =
    begin
        Ns, Vs = envsig(K), envtype(K)
        Ns0, Vs0 = envsig(X), envtype(X)
        env, (ptr, vals) = data(X)
        return K(liftenv(Ns, Vs, Ns0, Vs0, env), ptr, vals)
    end

@generated function liftenv{names,Vs,names0,Vs0}(
    ::Type{Val{names}}, ::Type{Vs}, ::Type{Val{names0}}, ::Type{Vs0}, env::Vs0)
    if names == names0 && Vs == Vs0
        return :(env)
    else
        items = []
        indexes = Dict{Symbol,Int}([name => n for (n, name) in enumerate(names0)])
        for (name, V) in zip(names, Vs.parameters)
            index = indexes[name]
            push!(items, :(env[$index]))
        end
        return :(convert(Vs, $(Expr(:tuple, items...))))
    end
end

#
# Adding parameters and context.
#

bind_rel{T}(X::Iso{T}) =
    begin
        val = data(X)
        Rel{T}(1, T[val])
    end

bind_rel{Ns,Vs,T}(X::Env{Ns,Vs,T}) =
    begin
        env, val = data(X)
        EnvRel{Ns,Vs,T}(env, 1, T[val])
    end

bind_env{Ns,Vs,T}(X::Iso{T}, ::Type{Ns}, env::Vs) =
    begin
        val = data(X)
        Env{Ns,Vs,T}(env, val)
    end

@generated function bind_env{Ns1,Vs1,Ns2,Vs2,T}(X::Env{Ns1,Vs1,T}, ::Type{Ns2}, env2::Vs2)
    Ns, Vs, idxmap = _envjoin(Ns1,Vs1,Ns2,Vs2)
    return quote
        env1, val = data(X)
        env = $(Expr(:tuple, [idx > 0 ? :( env1[$idx] ) : :( env2[-$idx] ) for idx in idxmap]...))
        return Env{$Ns,$Vs,$T}(env, val)
    end
end

@generated function bind_env{Ns1,Vs1,Ns2,Vs2,T}(X::EnvRel{Ns1,Vs1,T}, ::Type{Ns2}, env2::Vs2)
    Ns, Vs, idxmap = _envjoin(Ns1,Vs1,Ns2,Vs2)
    return quote
        env1, ptr, vals = data(X)
        env = ($([idx > 0 ? :( env1[$idx] ) : :( env2[-$idx] ) for idx in idxmap]...))
        return EnvRel{$Ns,$Vs,$T}(env, ptr, vals)
    end
end

_envjoin(Ns1,Vs1,Ns2,Vs2) =
    begin
        names1 = Ns1.parameters[1]
        names2 = Ns2.parameters[1]
        idxs = Dict{Symbol,Int}()
        for (j, name) in enumerate(names1)
            idxs[name] = j
        end
        for (j, name) in enumerate(names2)
            idxs[name] = -j
        end
        names = sort(collect(Symbol, keys(idxs)))
        idxmap = Int[idxs[name] for name in names]
        types = Type[idx > 0 ? Vs1.parameters[idx] : Vs2.parameters[-idx] for idx in idxmap]
        Ns = Val{tuple(names...)}
        Vs = Tuple{types...}
        return Ns, Vs, idxmap
end

#
# Duplicate input: W{T} -> W_2{W_1{T}}.
#

dup{K<:Kind,K0<:Kind}(KK::Type{Iso{K}}, X::K0) = KK(lift(K, X))

dup{K<:Kind,K0<:Rel}(KK::Type{Rel{K}}, X::K0) =
    begin
        ptr, vals = data(X)
        KK(ptr, K[lift(K, K0(n, vals)) for n = 1:length(vals)])
    end
dup{K<:Kind,K0<:EnvRel}(KK::Type{Rel{K}}, X::K0) =
    begin
        env, (ptr, vals) = data(X)
        KK(ptr, K[lift(K, K0(env, n, vals)) for n = 1:length(vals)])
    end

dup{Ns,Vs,K<:Kind,K0<:Env}(KK::Type{Env{Ns,Vs,K}}, X::K0) =
    begin
        Ns0, Vs0 = envsig(K0), envtype(K0)
        env, val = data(X)
        return KK(liftenv(Ns, Vs, Ns0, Vs0, env), lift(K, X))
    end
dup{Ns,Vs,K<:Kind,K0<:EnvRel}(KK::Type{Env{Ns,Vs,K}}, X::K0) =
    begin
        Ns0, Vs0 = envsig(K0), envtype(K0)
        env, (ptr, vals) = data(X)
        return KK(liftenv(Ns, Vs, Ns0, Vs0, env), lift(K, X))
    end

dup{Ns,Vs,K<:Kind,K0<:EnvRel}(KK::Type{EnvRel{Ns,Vs,K}}, X::K0) =
    begin
        Ns0, Vs0 = envsig(K0), envtype(K0)
        env, (ptr, vals) = data(X)
        return KK(
            liftenv(Ns, Vs, Ns0, Vs0, env),
            ptr,
            K[lift(K, K0(env, n, vals)) for n = 1:length(vals)])
    end

#
# Distributive law: W_2{M_1{T}} -> M_1{W_2{T}}.
#

dist{T}(XX::Iso{Iso{T}}) = XX
dist{T}(XX::Iso{Opt{T}}) =
    let K = Iso{T},
        X = data(XX)
        isnull(X) ? Opt{K}() : Opt{K}(K(get(X)))
    end
dist{T}(XX::Iso{Seq{T}}) =
    let K = Iso{T}
        Seq{K}(K[K(x) for x in data(XX)])
    end

dist{T}(XX::Rel{Iso{T}}) =
    begin
        K = Rel{T}
        ptr, Xs = data(XX)
        return Iso{K}(K(ptr, T[data(X) for X in Xs]))
    end
dist{T}(XX::Rel{Opt{T}}) =
    begin
        K = Rel{T}
        ptr, Xs = data(XX)
        if isnull(Xs[ptr])
            return Opt{K}()
        else
            LXs = sub(Xs, 1:ptr)
            RXs = sub(Xs, ptr+1:endof(Xs))
            vals = T[]
            for X in LXs
                if !isnull(X)
                    push!(vals, get(X))
                else
                    ptr -= 1
                end
            end
            for X in RXs
                if !isnull(X)
                    push!(vals, get(X))
                end
            end
            return Opt{K}(K(ptr, vals))
        end
    end
dist{T}(XX::Rel{Seq{T}}) =
    begin
        K = Rel{T}
        ptr, Xs = data(XX)
        if isempty(Xs)
            return Seq{K}(K{T}[])
        else
            vals = T[]
            lptr = 0
            rptr = 0
            for n in eachindex(Xs)
                if n == ptr
                    lptr = endof(vals)+1
                end
                append!(vals, data(Xs[n]))
                if n == ptr
                    rptr = endof(vals)
                end
            end
            return Seq{K}(K[K(n, vals) for n=lptr:rptr])
        end
    end

dist{Ns,Vs,T}(XX::Env{Ns,Vs,Iso{T}}) =
    begin
        K = Env{Ns,Vs,T}
        env, X = data(XX)
        return Iso{K}(K(env, data(X)))
    end
dist{Ns,Vs,T}(XX::Env{Ns,Vs,Opt{T}}) =
    begin
        K = Env{Ns,Vs,T}
        env, X = data(XX)
        return isnull(X) ? Opt{K}() : Opt{K}(K(env, get(X)))
    end
dist{Ns,Vs,T}(XX::Env{Ns,Vs,Seq{T}}) =
    begin
        K = Env{Ns,Vs,T}
        env, X = data(XX)
        return Seq{K}(K[K(env, x) for x in X])
    end

dist{Ns,Vs,T}(XX::EnvRel{Ns,Vs,Iso{T}}) =
    begin
        K = EnvRel{Ns,Vs,T}
        env, (ptr, Xs) = data(XX)
        return Iso{K}(K(env, ptr, T[data(X) for X in Xs]))
    end
dist{Ns,Vs,T}(XX::EnvRel{Ns,Vs,Opt{T}}) =
    begin
        K = EnvRel{Ns,Vs,T}
        env, (ptr, Xs) = data(XX)
        if isnull(Xs[ptr])
            return Opt{K}()
        else
            LXs = sub(Xs, 1:ptr)
            RXs = sub(Xs, ptr+1:endof(Xs))
            vals = T[]
            for X in LXs
                if !isnull(X)
                    push!(vals, get(X))
                else
                    ptr -= 1
                end
            end
            for X in RXs
                if !isnull(X)
                    push!(vals, get(X))
                end
            end
            return Opt{K}(K(env, ptr, vals))
        end
    end
dist{Ns,Vs,T}(XX::EnvRel{Ns,Vs,Seq{T}}) =
    begin
        K = EnvRel{Ns,Vs,T}
        env, (ptr, Xs) = data(XX)
        if isempty(Xs)
            return Seq{K}(K{T}[])
        else
            vals = T[]
            lptr = 0
            rptr = 0
            for n in eachindex(Xs)
                if n == ptr
                    lptr = endof(vals)+1
                end
                append!(vals, data(Xs[n]))
                if n == ptr
                    rptr = endof(vals)
                end
            end
            return Seq{K}(K[K(env, n, vals) for n=lptr:rptr])
        end
    end

#
# Flatten output: M_1{M_2{T}} -> M{T}.
#

flat{T}(XX::Iso{Iso{T}}) = data(XX)
flat{T}(XX::Iso{Opt{T}}) = data(XX)
flat{T}(XX::Iso{Seq{T}}) = data(XX)

flat{T}(XX::Opt{Iso{T}}) = isnull(XX) ? Opt{T}() : Opt{T}(data(get(XX)))
flat{T}(XX::Opt{Opt{T}}) = isnull(XX) || isnull(get(XX)) ? Opt{T}() : Opt{T}(get(get(XX)))
flat{T}(XX::Opt{Seq{T}}) = isnull(XX) ? Seq{T}(T[]) : get(XX)

flat{T}(XX::Seq{Iso{T}}) = Seq{T}(T[data(X) for X in XX])
flat{T}(XX::Seq{Opt{T}}) =
    let vals = T[]
        for X in XX
            if !isnull(X)
                push!(vals, get(X))
            end
        end
        Seq{T}(vals)
    end
flat{T}(XX::Seq{Seq{T}}) =
    let vals = T[]
        for X in XX
            append!(vals, data(X))
        end
        Seq{T}(vals)
    end

#
# Evaluation.
#

codegen{I<:Kind,I0<:Kind}(F, ::Type{I0}, X, ::Type{I}) =
    if I <: I0
        codegen(F, X, I0)
    else
        codegen(F, :( lift($I0, $X) ), I0)
    end

codegen(F, X, I) = F(X, I)

#
# Mapping over a structure.
#

codegen_map{I<:Kind,O<:Kind,K<:Kind}(F, ::Type{I}, ::Type{O}, XX, ::Type{Iso{K}}) =
    quote
        Iso{$O}($(codegen(F, I, :( data($XX) ), K)))
    end

codegen_map{I<:Kind,O<:Kind,K<:Kind}(F, ::Type{I}, ::Type{O}, XX, ::Type{Opt{K}}) =
    let
        @gensym xx
        quote
            $xx = data($XX)::Nullable{$K}
            isnull($xx) ? Opt{$O}() : Opt{$O}($(codegen(F, I, :( get($xx) ), K)))
        end
    end

codegen_map{I<:Kind,O<:Kind,K<:Kind}(F, ::Type{I}, ::Type{O}, XX, ::Type{Seq{K}}) =
    let
        @gensym X
        quote
            Seq{$O}($O[$(codegen(F, I, X, K)) for $X in data($XX) ])
        end
    end

codegen_map{I<:Kind,O<:Kind,K<:Kind}(F, ::Type{I}, ::Type{O}, XX, ::Type{Rel{K}}) =
    let
        @gensym ptr X Xs
        quote
            $ptr, $Xs = data($XX)::RelData{$K}
            Rel{$O}($ptr, $O[$(codegen(F, I, X, K)) for $X in $Xs])
        end
    end

codegen_map{I<:Kind,O<:Kind,Ns,Vs,K<:Kind}(F, ::Type{I}, ::Type{O}, XX, ::Type{Env{Ns,Vs,K}}) =
    let
        @gensym env X
        quote
            $env, $X = data($XX)::EnvData{$Ns,$Vs,$K}
            Env{$Ns,$Vs,$O}($env, $(codegen(F, I, X, K)))
        end
    end

codegen_map{I<:Kind,O<:Kind,Ns,Vs,K<:Kind}(F, ::Type{I}, ::Type{O}, XX, ::Type{EnvRel{Ns,Vs,K}}) =
    let
        @gensym env ptr X Xs
        quote
            $env, ($ptr, $Xs) = data($XX)::EnvRelData{$Ns,$Vs,$K}
            EnvRel{$Ns,$Vs,$O}($env, $ptr, $O[$(codegen(F, I, X, K)) for $X in $Xs])
        end
    end

#
# Composition.
#

codegen_compose{I<:Kind,I1<:Kind,O1<:Kind,I2<:Kind,O2<:Kind}(
    G, ::Type{I2}, ::Type{O2}, F, ::Type{I1}, ::Type{O1}, X, ::Type{I}) =
    begin
        I2I1 = Kind(I2, I1)
        O1I2 = Kind(O1, I2)
        @gensym XX XY YX YY Y
        quote
            $XX = dup($I2I1, $X)
            $XY = $( codegen_map(F, I1, O1, XX, I2I1) )
            $YX = dist($XY)
            $YY = $( codegen_map(G, I2, O2, YX, O1I2) )
            $Y = flat($YY)
            $Y
        end
    end

# Faster shortcuts.

codegen_compose{A,B,C}(G, ::Type{Iso{B}}, ::Type{Iso{C}}, F, ::Type{Iso{A}}, ::Type{Iso{B}}, X, ::Type{Iso{A}}) =
    codegen(G, Iso{B}, codegen(F, Iso{A}, X, Iso{A}), Iso{B})

codegen_compose{A,B,C}(G, ::Type{Iso{B}}, ::Type{Iso{C}}, F, ::Type{Iso{A}}, ::Type{Seq{B}}, X, ::Type{Iso{A}}) =
    begin
        @gensym Y
        quote
            Seq{$C}(
                $C[data($(codegen(G, Iso{B}, :( Iso{$B}($Y) ), Iso{B}))) for $Y in $(codegen(F, Iso{A}, X, Iso{A}))])
        end
    end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# Check that two vectors or matrices have matching dimensions
function check_dimensions(X::AbstractVecOrMat, Y::AbstractVecOrMat)
    size(X, 1) == size(Y, 1) ||
        throw(DimensionMismatch("The two objects don't have the same " *
                                "number of rows."))
    return nothing
end


text:
abstract type Element <: OrgComponent end

include("lesserelements.jl")
include("greaterelements.jl")
include("affiliatedkeywords.jl")
include("sectioning.jl")

# Taken out of greaterelements.jl for load order reasons
mutable struct InlineTask <: GreaterElement
    heading::Heading
end


text:
struct Identity <: LocalOperator
end

kernelvals(biop::Identity, x) = Nothing
integrand(op::Identity, kernel, x, g, f) = dot(f[1], g[1])
# scalartype(op::Identity) = Union{}
scalartype(op::Identity) = Float64

struct NCross <: LocalOperator
end

kernelvals(op::NCross, mp) = Nothing()
integrand(op::NCross, kernel, x, g, f) = dot(g[1], normal(x) × f[1])
scalartype(op::NCross) = Union{}

function quaddata(op::LocalOperator, g::RTRefSpace, f::RTRefSpace, tels, bels)
    u, w = trgauss(6)
    return [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
end

function quaddata(op::LocalOperator, g::NDRefSpace, f::NDRefSpace, tels, bels)
    u, w = trgauss(6)
    return [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
end

function quaddata(op::LocalOperator, g::subReferenceSpace, f::subReferenceSpace, tels, bels)
    u, w = trgauss(6)
    return [(w[i],SVector(u[1,i],u[2,i])) for i in 1:length(w)]
end

# function quaddata(op::LocalOperator, g::Nd4DRefSpace, f::Nd4DRefSpace, tels, bels)
#     o, x, y, z = CompScienceMeshes.euclidianbasis(3)
#     reftet = simplex(o,x,y,z)
#     qps = quadpoints(reftet, rul)
#     [(w, parametric(p)) for (p,w) in qps]
# end

function quaddata(op::LocalOperator, g::NDLCCRefSpace, f::NDLCCRefSpace, tels, bels)
     o, x, y, z = CompScienceMeshes.euclidianbasis(3)
     reftet = simplex(x,y,z,o)
     qps = quadpoints(reftet, 6)
     [(w, parametric(p)) for (p,w) in qps]
end

function quaddata(op::LocalOperator, g::NDLCDRefSpace, f::NDLCDRefSpace, tels, bels)
     o, x, y, z = CompScienceMeshes.euclidianbasis(3)
     reftet = simplex(x,y,z,o)
     qps = quadpoints(reftet, 6)
     [(w, parametric(p)) for (p,w) in qps]
end

function quaddata(op::LocalOperator, g::NDLCDRefSpace, f::NDLCCRefSpace, tels, bels)
     o, x, y, z = CompScienceMeshes.euclidianbasis(3)
     reftet = simplex(x,y,z,o)
     qps = quadpoints(reftet, 6)
     [(w, parametric(p)) for (p,w) in qps]
end

quaddata(op::LocalOperator, g::LagrangeRefSpace, f::LagrangeRefSpace,
        tels::Vector, bels::Vector) = quaddata(op, g, f, tels, bels, Val{dimension(tels[1])})

function quaddata(op::LocalOperator, g::LagrangeRefSpace, f::LagrangeRefSpace,
        tels::Vector, bels::Vector, dim::Type{Val{1}})

    u, w = legendre(6, 0.0, 1.0)
    [(w[i],u[i]) for i in eachindex(w)]
end

function quaddata(op::LocalOperator, g::LagrangeRefSpace, f::LagrangeRefSpace,
        tels::Vector, bels::Vector, dim::Type{Val{2}})

    u, w = trgauss(6)
    [(w[i], SVector(u[1,i], u[2,i])) for i in 1:length(w)]
end

function quadrule(op::LocalOperator, ψ::RefSpace, ϕ::RefSpace, τ, qd)
    q = qd[1]
    w, p = q[1], neighborhood(τ,q[2])
    a = (w, p, ψ(p), ϕ(p))
    A = Vector{typeof(a)}(undef,length(qd))
    A[1] = a
    for i in 2:length(qd)
        q = qd[i]
        w, p = q[1], neighborhood(τ,q[2])
        A[i] = (w, p, ψ(p), ϕ(p))
    end
    return A
end


text:
using Test
using Multilane
using MCTS
using RobustMCTS
using POMDPs
using POMDPModelTools
using POMDPSimulators
using Random
using Test
using Cairo

#Set up problem configuration
nb_lanes = 4
pp = PhysicalParam(nb_lanes,lane_length=100.) #2.=>col_length=8
_discount = 1.
nb_cars=10

rmodel = NoCrashRewardModel()

dmodel = NoCrashIDMMOBILModel(nb_cars, pp)

mdp = NoCrashMDP(dmodel, rmodel, _discount, true);

rng = MersenneTwister(5)

s = initial_state(mdp::NoCrashMDP, rng)

rsolver = RobustMCTSSolver(
    depth=20,
    c = 100.0,
    n_iterations=500,
    k_state=4.0,
    alpha_state=1/8,
    k_nature=4.0,
    alpha_nature=1/8,
    rollout_solver=SimpleSolver(),
    rollout_nature=StochasticBehaviorNoCrashMDP(mdp))

solver = RobustMLSolver(rsolver)
policy = solve(solver, mdp)

sim = HistoryRecorder(rng=rng, max_steps=100) # initialize a random number generator

simulate(sim, mdp, policy, s)

# check for crashes
for i in 1:length(sim.state_hist)-1
    if is_crash(mdp, sim.state_hist[i], sim.state_hist[i+1])
        println("""
        Crash!
        mdp = $mdp
        s = $(sim.state_hist[i])
        a = $(sim.action_hist[i])
        sp = $(sim.state_hist[i+1])
        """)
    end
    @test !is_crash(mdp, sim.state_hist[i], sim.state_hist[i+1])
end


text:
using DataFrames
using CSV
using Statistics
using Random
using Plots.Measures
using Plots, StatsPlots

const FULL_FILE_PROCESSED = joinpath("data","Melbourne_housing_FULL_processed_2.csv")

"""
    currency_format(num)::String

Takes a number, returns it as a string in currency (\$AUD) form. Taken from [here](https://stackoverflow.com/questions/52213829/in-julia-insert-commas-into-integers-for-printing-like-python-3-6/52213830#52213830).

# Examples
```jldoctest
julia>a = 10000
10000
julia>currency_format(a)
\$10,000
```
"""
currency_format(num)::String = "\$" * replace(num |> string, r"(?<=[0-9])(?=(?:[0-9]{3})+(?![0-9]))" => ",")

###
# Distance to CBD
###

"""
    plot_distance_price(df_full::AbstractDataFrame)

Plots price in relation to distance as a box-and-whisker plot (range=5), the latter grouped in 4km segments.

X = distance from CBD (km), Y = sale price for properties within the distance group.
"""
function plot_distance_price(df_full::DataFrame)
    df_dp = transform!(dropmissing!(select(df_full, :Price, :Distance), [:Price, :Distance]), :Distance => ByRow(x->convert(Int, x ÷ 4 * 4)) => :Distance_grouped)

    y_ticks = 2500000:2500000:12500000

    return @df df_dp boxplot(:Distance_grouped, :Price, xlabel= "Distance from CBD (km)", ylabel = "Sale Price (\$AUD)", yticks=(y_ticks, currency_format.(y_ticks)), outliers=false, range=5, width=1, label="", plot_title="Price vs. Distance From CBD")
end

###
# Price vs. Rooms
###

"""
    plot_rooms_price(df_full::AbstractDataFrame)

Plots sale price against number of rooms as a box-and-whisker plot (range=5).

X = number of rooms, Y = sale price for properties with that number of rooms. Y scale is log10.
"""
function plot_rooms_price(df_full::DataFrame)
    df_rp_init = sort!(transform!(select(df_full, :Price, :Rooms) |> dropmissing!, :Price => ByRow(x->x÷500000*500000+1) => :Price_adj), :Price_adj)

    yticks=10 .^ collect(5:7)

    return @df df_rp_init plot(:Rooms, :Price, seriestype=:box, outliers=false, range=5, label="", xlabel = "Rooms in Property", ylabel = "Price (\$AUD)", xticks=2:2:20, yscale=:log10, yticks=(yticks, currency_format.(yticks)), ylims=(10^4.9, 10^7), plot_title="Rooms vs. Price")
end

###
# Price vs. Rooms and CBD Distance
###

"""
    plot_rooms_distance_price(df_full::AbstractDataFrame)

Plots sale price against number of rooms and distance from CBD as three-dimensional surface plot.

X = distance from CBD (km), Y = num. of rooms, Z = mean sale price for props. w/ that num. of rooms and distance. Z scale is log10.

Note that Z tick labels are incorrectly placed. Their true position should correspond to the grid.
"""
function plot_rooms_distance_price(df_full::DataFrame)
    df_rdp = combine(groupby(select(df_full, :Rooms, :Distance => ByRow(floor) => :Distance, :Price) |> dropmissing!, [:Rooms, :Distance]), :Price => mean => :Price_mean)
    zticks=10 .^ collect(5:7)
    return @df df_rdp plot(:Distance, :Rooms, :Price_mean, seriestype=:surface, zscale=:log10, camera = (50,70), zlims=(5000, 10000000), size=(800, 500), xlims=(0, 50), ylims=(0, 16), yflip=true, xflip=true, xlabel="Distance from CBD (km)                         ", ylabel="          Number of Rooms", zlabel="Average Sale Price (\$AUD)", plot_title="Room and CBD Distance as a Factor in Price", left_margin=5mm, right_margin=5mm, zticks=(zticks, currency_format.(zticks)))
end

"""
    generate_task2_df(filename=FULL_FILE_PROCESSED)

Generates a `DataFrame` from a CSV file for use throughout this task.
"""
generate_task2_df(filename=FULL_FILE_PROCESSED) = DataFrame(CSV.File(filename))


text:
module MyMatching

function matching(prop_prefs,resp_prefs,reverse)

    if reverse == false

        m_prefs = prop_prefs
        f_prefs = resp_prefs
        mm = length(m_prefs)
        nn = length(f_prefs)

        for o in 1:mm
            push!(m_prefs[o],0)
        end

        for k in 1:nn
            push!(f_prefs[k],0)
        end

        matched_m = zeros(Int64,mm)
        matched_f = zeros(Int64,nn)
        accepted_m = zeros(Int64,mm)
        order = ones(Int64,mm)

        while sum(accepted_m) < mm          #男性全員がマッチするまで

            for i in 1:mm

                if accepted_m[i] == 1     #男性iがマッチ済み

                else　　　#男性iがマッチ済みでない

                    if m_prefs[i][order[i]] == 0       #男性が独身を希望

                        accepted_m[i] = 1

                    else　　　　#男性が女性を希望

                        if findfirst(f_prefs[m_prefs[i][order[i]]],i) == 0      #男性の希望する女性の選択肢の中に男性がいない

                            order[i] +=1           #男性の志望順位を一つ下げる

                        else　　　　　　#男性が女性の評価対象

                            if matched_f[m_prefs[i][order[i]]] == 0       #男性の希望する女性のもとに誰もいないとき

                                matched_m[i] = m_prefs[i][order[i]]
                                matched_f[m_prefs[i][order[i]]] = i
                                accepted_m[i] = 1

                            else　　　　　　　　　#女性が別の男性とマッチ済み

                                if findfirst(f_prefs[m_prefs[i][order[i]]],i) < findfirst(f_prefs[m_prefs[i][order[i]]],matched_f[m_prefs[i][order[i]]])

                                    order[matched_f[m_prefs[i][order[i]]]] +=1         #取って代わられた男性の志望順位を一つ下げる
                                    accepted_m[i] = 1
                                    matched_m[matched_f[m_prefs[i][order[i]]]] = 0
                                    accepted_m[matched_f[m_prefs[i][order[i]]]] = 0
                                    matched_m[i] = m_prefs[i][order[i]]
                                    matched_f[m_prefs[i][order[i]]] = i

                                else　　　　　#マッチ済みの男性のほうがいい

                                    order[i] +=1

                                end

                            end

                        end

                    end

                end

            end

        end

        return matched_m, matched_f

    elseif reverse == true

        m_prefs = prop_prefs
        f_prefs = resp_prefs
        mm = length(m_prefs)
        nn = length(f_prefs)

        for o in 1:mm
            push!(m_prefs[o],0)
        end

        for k in 1:nn
            push!(f_prefs[k],0)
        end

        matched_m = zeros(Int64,mm)
        matched_f = zeros(Int64,nn)
        accepted_f = zeros(Int64,nn)
        order = ones(Int64,nn)

        while sum(accepted_f) < nn          #女性全員がマッチするまで

            for i in 1:nn

                if accepted_f[i] == 1     #女性iがマッチ済み

                else　　　#女性iがマッチ済みでない

                    if f_prefs[i][order[i]] == 0       #女性iが独身を希望

                        accepted_f[i] = 1

                    else　　　　#女性が男性を希望

                        if findfirst(m_prefs[f_prefs[i][order[i]]],i) == 0      #男性の希望する女性の選択肢の中に男性がいない

                            order[i] +=1           #男性の志望順位を一つ下げる

                        else　　　　　　#男性が女性の評価対象

                            if matched_m[f_prefs[i][order[i]]] == 0       #男性の希望する女性のもとに誰もいないとき

                                matched_f[i] = f_prefs[i][order[i]]
                                matched_m[f_prefs[i][order[i]]] = i
                                accepted_f[i] = 1

                            else　　　　　　　　　#女性が別の男性とマッチ済み

                                if findfirst(m_prefs[f_prefs[i][order[i]]],i) < findfirst(m_prefs[f_prefs[i][order[i]]],matched_m[f_prefs[i][order[i]]])

                                    order[matched_m[f_prefs[i][order[i]]]] +=1         #取って代わられた男性の志望順位を一つ下げる
                                    accepted_f[i] = 1
                                    matched_f[matched_m[f_prefs[i][order[i]]]] = 0
                                    accepted_f[matched_m[f_prefs[i][order[i]]]] = 0
                                    matched_f[i] = f_prefs[i][order[i]]
                                    matched_m[f_prefs[i][order[i]]] = i

                                else　　　　　#マッチ済みの男性のほうがいい

                                    order[i] +=1

                                end

                            end

                        end

                    end

                end

            end

        end

        return matched_m, matched_f

    end

end

function matching(prop_prefs,resp_prefs,caps,reverse)

    if reverse == false           #学生　→　大学

        mm = length(prop_prefs)
        nn = length(resp_prefs)
        for o in 1:mm
            push!(prop_prefs[o],0)
        end

    #resp_prefsに0を入れる必要はなかった

        indptr = Array{Int}(nn+1)
        indptr[1] = 1
        for k in 1:nn
            indptr[k+1] = indptr[k] + caps[k]
        end
        #indptrを関数内に内蔵！

        prop_matched = zeros(Int64,mm)
        resp_matched = zeros(Int64,sum(caps))
        prop_accepted = zeros(Int64,mm)
        order = ones(Int64,mm)

        while sum(prop_accepted) < mm       #学生全員がマッチするまで

            for i in 1:mm

                if prop_accepted[i] == 1    #学生iがマッチ済み

                else

                    if prop_prefs[i][order[i]] == 0   #学生iの第一志望が浪人

                        prop_accepted[i] = 1

                    else

                        if findfirst(resp_prefs[prop_prefs[i][order[i]]],i) == 0  #学生iが志望する大学の選択肢の中に学生iがいない

                            order[i] +=1    #学生iの志望順位を一つ下げる

                        else

                            if findfirst(resp_matched[indptr[prop_prefs[i][order[i]]]:indptr[prop_prefs[i][order[i]]+1]-1],0) !== 0
                               #学生iが志望する大学の受け入れ枠に空きがあるとき(0が受け入れ枠にある)

                                z = findfirst(resp_matched[indptr[prop_prefs[i][order[i]]]:indptr[prop_prefs[i][order[i]]+1]-1],0)
                                #0が左から何番目にあるか

                                resp_matched[indptr[prop_prefs[i][order[i]]]+z-1] = i     #0...0にiを左詰めで入れる
                                prop_matched[i] = prop_prefs[i][order[i]]
                                prop_accepted[i] = 1

                            else

                                passed = resp_matched[indptr[prop_prefs[i][order[i]]]:indptr[prop_prefs[i][order[i]]+1]-1]
                                #受け入れ保留されている学生たちのarray

                                ranking = [findfirst(resp_prefs[prop_prefs[i][order[i]]],passed[p]) for p in 1:length(passed)]
                                #彼らのランキング

                                worst = resp_prefs[prop_prefs[i][order[i]]][maximum(ranking)]
                                #保留されている中で最もランキングの低い学生の”名前”

                                if findfirst(resp_prefs[prop_prefs[i][order[i]]],i) > maximum(ranking)
                                   #最下位の人よりも低い

                                    order[i] += 1

                                else

                                    ww = findfirst(resp_matched[indptr[prop_prefs[i][order[i]]]:indptr[prop_prefs[i][order[i]]+1]-1],worst)
                                    #最下位の人が何番目にいるのか

                                    prop_accepted[worst] = 0
                                    prop_accepted[i] = 1
                                    prop_matched[worst] = 0
                                    prop_matched[i] = prop_prefs[i][order[i]]
                                    order[worst] += 1
                                    resp_matched[indptr[prop_prefs[i][order[i]]]+ww-1] = i

                                end

                            end

                        end

                    end

                end

            end

        end

        return prop_matched, resp_matched, indptr

    elseif reverse == true            #大学　→　学生

        mm = length(prop_prefs)
        nn = length(resp_prefs)
        for o in 1:nn
            push!(resp_prefs[o],0)
        end

        #今回はprop_prefsに0を入れる必要がない

        indptr = Array{Int}(nn+1)
        indptr[1] = 1
        for k in 1:nn
            indptr[k+1] = indptr[k] + caps[k]
        end
        #indptrを関数内に内蔵！

        prop_matched = zeros(Int64,mm)
        resp_matched = zeros(Int64,sum(caps))
        resp_accepted = zeros(Int64,nn)
        order = ones(Int64,nn)

        while sum(resp_accepted) < nn       #学生全員がマッチするまで

            for i in 1:nn

                if resp_accepted[i] == 1    #学生iがマッチ済み

                else

                    if resp_prefs[i][order[i]] == 0   #大学iの第一志望がいないとき

                        resp_accepted[i] = 1

                    elseif findfirst(prop_prefs[resp_prefs[i][order[i]]],i) == 0 && resp_prefs[i][order[i]] !== 0
                        #大学が志望する学生の選択肢の中に大学がない

                        order[i] += 1

                    else

                        if prop_matched[resp_prefs[i][order[i]]] == 0      #希望の学生が空いているとき

                            prop_matched[resp_prefs[i][order[i]]] = i
                            vacant = findfirst(resp_matched[indptr[i]:indptr[i+1]-1],0)   #大学の空いている番号
                            resp_matched[indptr[i]+vacant-1] = resp_prefs[i][order[i]]

                            if findfirst(resp_matched[indptr[i]:indptr[i+1]-1],0) == 0  #大学の定員が埋まったとき

                                resp_accepted[i] = 1

                            else

                                order[i] += 1

                            end

                        else

                            rival = prop_matched[resp_prefs[i][order[i]]]     #希望の学生がすでにマッチング済みの時の相手

                            if findfirst(prop_prefs[resp_prefs[i][order[i]]],rival) < findfirst(prop_prefs[resp_prefs[i][order[i]]],i)
                                #大学iはライバル大学よりも学生に好まれていないとき

                                order[i] += 1

                            else
                                #大学iがライバル大学よりも人気だったとき
                                losestu = findfirst(resp_matched[indptr[rival]:indptr[rival+1]-1],resp_prefs[i][order[i]])
                                #ライバルの定員のどこに学生が位置していたか
                                resp_matched[indptr[rival]+losestu-1] = 0
                                vacant = findfirst(resp_matched[indptr[i]:indptr[i+1]-1],0)   #大学の空いている番号
                                resp_matched[indptr[i]+vacant-1] = resp_prefs[i][order[i]]
                                resp_accepted[rival] = 0

                                if order[rival] < length(resp_prefs[rival])

                                    order[rival] += 1

                                end

                                prop_matched[resp_prefs[i][order[i]]] = i    # 最後に学生のマッチング相手を置き換える

                                if findfirst(resp_matched[indptr[i]:indptr[i+1]-1],0) == 0  #大学の定員が埋まったとき

                                    resp_accepted[i] = 1

                                else

                                    order[i] += 1

                                end

                            end

                        end

                    end

                end

            end

        end

        return prop_matched, resp_matched, indptr

    end

end

export matching

end


text:
# constants
const FAC = 1e-9 / (H_PLANCK() * LIGHT_SPEED() * AVOGADRO());

#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Oct-22: rename the function to photon
#     2022-Feb-02: fix documentation
#
#######################################################################################################################################################################################################
"""
This function convert energy to photons based the wavelength. Supported methods are

$(METHODLIST)

"""
function photon end

#######################################################################################################################################################################################################
#
# Changes made to this method
# General
#     2021-Oct-22: add a method to convert direct from number to number
#     2022-Feb-02: fix documentation
#
#######################################################################################################################################################################################################
"""

    photon(λ::FT, E::FT) where {FT<:AbstractFloat}

Return the number of moles of photons, given
- `λ` Wave length in `[nm]`, converted to `[m]` by FAC
- `E` Joules of energy

---
# Examples
```julia
phot = photon(500.0, 0.1);
phots = photon.([500.0, 400.0], [0.1, 0.1]);
```
"""
photon(λ::FT, E::FT) where {FT<:AbstractFloat} = E * λ * FT(FAC);

#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Oct-22: define function to convert photon back to energy
#     2022-Feb-02: fix documentation
#
#######################################################################################################################################################################################################
"""
This function convert photons to energy based the wavelength. Supported methods are

$(METHODLIST)

"""
function energy end

#######################################################################################################################################################################################################
#
# Changes made to this function
# General
#     2021-Oct-22: add a method to convert photon back to energy
#     2022-Feb-02: fix documentation
#
#######################################################################################################################################################################################################
"""

    energy(λ::FT, phot::FT) where {FT<:AbstractFloat}

Return the energy, given
- `λ` Wave length in `[nm]`, converted to `[m]` by FAC
- `phot` Number of moles of photon

---
# Examples
```julia
e = energy(500.0, 0.1);
es = energy.([500.0, 400.0], [0.1, 0.1]);
```
"""
energy(λ::FT, phot::FT) where {FT<:AbstractFloat} = phot / (λ * FT(FAC));


text:
"""
    test_frequencies() :: Nothing

Run consistency checks for current frequency implementation for meshes with 30 positive frequencies.
"""
function test_frequencies() :: Nothing

    # init test dummys
    list = get_mesh(rand(), 1.0, 30, 0.4)
    m    = Mesh(31, 31, 31, 31, list, list, list, list, list, list)
    w    = rand()
    v    = rand()
    vp   = rand()

    # test list bounds
    @testset "list bounds" begin
        @test list[1]   ≈ 0.0
        @test list[end] ≈ 1.0
    end

    # test frequency buffers
    @testset "buffers" begin
        @testset "s channel" begin
            b = get_buffer_s(w,     v,  vp, m); @test b.kernel == 4
            b = get_buffer_s(w,   Inf,  vp, m); @test b.kernel == 3
            b = get_buffer_s(w,     v, Inf, m); @test b.kernel == 2
            b = get_buffer_s(w,   Inf, Inf, m); @test b.kernel == 1
            b = get_buffer_s(Inf, Inf, Inf, m); @test b.kernel == 0
        end

        @testset "t channel" begin
            b = get_buffer_t(w,     v,  vp, m); @test b.kernel == 4
            b = get_buffer_t(w,   Inf,  vp, m); @test b.kernel == 3
            b = get_buffer_t(w,     v, Inf, m); @test b.kernel == 2
            b = get_buffer_t(w,   Inf, Inf, m); @test b.kernel == 1
            b = get_buffer_t(Inf, Inf, Inf, m); @test b.kernel == 0
        end

        @testset "u channel" begin
            b = get_buffer_u(w,     v,  vp, m); @test b.kernel == 4
            b = get_buffer_u(w,   Inf,  vp, m); @test b.kernel == 3
            b = get_buffer_u(w,     v, Inf, m); @test b.kernel == 2
            b = get_buffer_u(w,   Inf, Inf, m); @test b.kernel == 1
            b = get_buffer_u(Inf, Inf, Inf, m); @test b.kernel == 0
        end
    end

    return nothing
end


text:
using GeoClustering
using Test

@testset "GeoClustering.jl" begin
    # Write your tests here.
end


text:
using PyPlot
using Statistics

using SubmodularMaximization

num_trials = 100
num_agents = 20

f(x) = mean_area_coverage(x, floor(Integer, 20 * num_agents / 4))

struct CoverageParams
  num_sensors
  nominal_area
end

num_range =  collect(2:6)
area_range = collect(0.5:0.1:2.0)

param_set = [CoverageParams(num_sensors, nominal_area) for num_sensors in num_range,
             nominal_area in area_range]

function coverage_performance_ratio(params::CoverageParams)
  println("Testing: sensors ($(params.num_sensors)), nominal area ($(params.nominal_area))")

  sensor_radius = sqrt(params.nominal_area / (num_agents * pi))
  station_radius = 2 * sensor_radius

  agent_specification = CircleAgentSpecification(sensor_radius, station_radius,
                                           params.num_sensors)

  problems = map(1:num_trials) do unused
    agents = generate_agents(agent_specification, num_agents)
    ExplicitPartitionProblem(f, agents)
  end

  myopic_mean = mean(map(x->solve_myopic(x).value, problems))
  sequential_mean = mean(map(x->solve_sequential(x).value, problems))

  myopic_mean / sequential_mean
end

performance_ratios = map(coverage_performance_ratio, param_set)

_, best_index = findmin(performance_ratios)
best_params = param_set[best_index]

println("Best ratio: $(performance_ratios[best_index])")
println("Params: sensors ($(best_params.num_sensors)), nominal area ($(best_params.nominal_area))")

figure()
PyPlot.plt.hist(performance_ratios[:], 10)

figure()
pcolor(performance_ratios, cmap = PyPlot.cm.Blues)
colorbar()

nothing

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# A Frisch-Newton algorithm as described in http://projecteuclid.org/euclid.ss/1030037960
# qreg_ip_coef is a modified lp_fnm of Daniel Morillo & Roger Koenker
# Translated from Ox to Matlab by Paul Eilers 1999
# Modified by Roger Koenker 2000
# More changes by Paul Eilers 2004
# QuantileRegressions.jl developers 2015
#
# The functions in this file is released under the BSD-3 license, although the original
# versions were not. Permissing was given to Patrick Kofod Mogensen by Roger Koenker and
# Paul Eilers on April 24 2015 via E-mail. Original mail correspondences can be provided
# upon request.

struct IP <: Solver
    cholesky::Bool
end
IP() = IP(false)

function bound!(d, x, dx)
# Fill vector with allowed step lengths
# Replace with -x/dx for negative dx
    d .= 1e20
    @inbounds for i = 1:length(dx)
        if dx[i] < 0.0
            d[i] = -x[i] / dx[i]
        end
    end
    return d
end

function qreg_coef(Y, X::Matrix, p, method::IP)
# input: X is an n x k matrix of exogenous regressors,
#        Y is an n x 1 vector of outcome variables
#        p \in (0,1) is the quantile of interest
# Output: p^{th} regression quantiles.
# Construct the dual problem of quantile regression
c = -Y
T = eltype(X)

m = size(X, 1)

x = fill((1 - p), m)
b = X'x
# Solve a linear program by the interior point method:
# min(c * u), s.t. A * x = b and 0 < x < u
# An initial feasible solution has to be provided as x

# Set some constants
beta = 0.99995
small = sqrt(eps(eltype(c)))
max_it = 500
n, m = size(X)
# Generate inital feasible point
s = 1 .- x
if method.cholesky
    y = -((X'X)\(X'Y))
else
    y = -X\Y
end
dy = copy(y)
r = c - X*y
# z and w are the "residuals" but on either
# side of the "check" we need to flip the sign
#    BLAS.axpy!(small, (r .== 0.0).*1.0, r)
#    z = r .* (r .> small)
#    w = abs.(z .- r)
z = copy(r)
w = copy(r)
for i ∈ eachindex(r)
    ri = r[i]
    z[i] = max(ri, 0)
    w[i] = max(-ri, 0)
    if abs(ri) <= small
        z[i] = z[i] + small
        w[i] = w[i] + small
    end
end

gap = dot(z, x) + dot(s, w)

# set up caches
Xtmp = copy(X)
if method.cholesky
    XtX  = similar(X, m, m)
end
xinv, xi = copy(x), copy(x)
sinv = copy(s)
q = copy(z)
dx, ds, dz, dw = copy(w), copy(w), copy(w), copy(w)
fx, fs, fz, fw = copy(w), copy(w), copy(w), copy(w)
dxdz, dsdw = copy(w), copy(w)
tmp = copy(w)
if method.cholesky
    rtmp = copy(r)
    Xtqr = similar(r, m)
else
    qr_tmp = copy(q)
    rhs = copy(r)
end
# Start iterations
for it = 1:max_it
    #   Compute affine step
    @. q = 1 / (z / x + w / s)
    @. r = z - w

    if method.cholesky
        @. Xtmp = q * X
        mul!(XtX, Xtmp', X)
        F = cholesky!(Symmetric(XtX))
        mul!(Xtqr, Xtmp', r)
        dy .= F\Xtqr
        #ldiv!(dy, F, Xtqr)
    else
        @. qr_tmp = sqrt(q)
        Q = Diagonal(qr_tmp) # Very efficient to do since Q diagonal
        AQtF = qr(mul!(Xtmp, Q, X)) # PE 2004
        @. rhs = qr_tmp*r
        dy .= AQtF\rhs
    end

    mul!(tmp, X, dy)
    @. dx = q*(tmp - r)
    @. ds = -dx
    @. dz = -z * (1 + dx / x)
    @. dw = -w * (1 + ds / s)
    # Compute maximum allowable step lengths
    bound!(fx, x, dx)
    bound!(fs, s, ds)
    bound!(fw, w, dw)
    bound!(fz, z, dz)
    @. tmp = min(fx, fs)
    fp = min(beta*minimum(tmp), 1)
    @. tmp = min(fw, fz)
    fd = min(beta*minimum(tmp), 1)
    # If full step is feasible, take it. Otherwise modify it
    if min(fp, fd) < 1.0
        # Update mu
        g0 = dot(z, x) + dot(w, s)
        gγPγD = T(0)
        for i = 1:length(z)
          @inbounds gγPγD += (z[i] + fd*dz[i])*(x[i] + fp*dx[i]) + (w[i] + fd*dw[i])*(s[i]+fp*ds[i])
        end
        mu = (gγPγD / g0)^3 * g0/(2n)
        # Compute modified step
        @. dxdz = dx * dz
        @. dsdw = ds * dw
        @. xinv = 1 / x
        @. sinv = 1 / s
        @. xi = mu * (xinv - sinv)
        if method.cholesky
            @. rtmp = r + dxdz - dsdw - xi
            mul!(Xtqr, Xtmp', rtmp)
            dy .= F\Xtqr
        else
            @. qr_tmp = qr_tmp * (dxdz .- dsdw .- xi)
            BLAS.axpy!(1.0, qr_tmp, rhs) # no gemv-wrapper gemv(Q, (dxdz - dsdw - xi), rhs,1,1,n)?
            dy .= AQtF\rhs
        end

        mul!(tmp, X, dy)
        @. dx = q * (tmp + xi - r - dxdz + dsdw)
        @. ds = -dx
        @. dz = (mu - z * dx)*xinv - z - dxdz
        @. dw = (mu - w * ds)*sinv - w - dsdw

        # Compute maximum allowable step lengths
        bound!(fx, x, dx)
        bound!(fs, s, ds)
        bound!(fw, w, dw)
        bound!(fz, z, dz)
        @. tmp = min(fx, fs)
        fp = min(beta*minimum(tmp), 1)
        @. tmp = min(fw, fz)
        fd = min(beta*minimum(tmp), 1)
    end

    # Take the steps
    BLAS.axpy!(fp, dx, x)
    BLAS.axpy!(fp, ds, s)
    BLAS.axpy!(fd, dy, y)
    BLAS.axpy!(fd, dw, w)
    BLAS.axpy!(fd, dz, z)

    gap = dot(z, x) + dot(s, w)

    if gap < small || isnan(gap)
        break
    end
end

return -y
end


text:
module Day10

using Advent: save_input
using DataStructures: CircularBuffer

function get_jolts(filename)
    nums = [parse(Int, x) for x in eachline(filename)] |> sort!
    pushfirst!(nums, 0)
    push!(nums, nums[end] + 3)
end

function solve_a(filename)
    nums = get_jolts(filename) |> diff
    count(==(1), nums) * count(==(3), nums)
end

function solve_b(filename)
    nums = get_jolts(filename)
    memo = CircularBuffer{Int}(3)
    append!(memo, [1, 1, 1])
    for i in (lastindex(nums) - 3):-1:1
        push!(memo, sum(memo[4 - j] for j in Base.OneTo(3) if (nums[i + j] - nums[i]) < 4))
    end
    memo[end]
end

solve(filename = save_input(2020, 10)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# Adaptive MvNormal random walk.
# See: https://www.tandfonline.com/doi/abs/10.1198/jcgs.2009.06134
mutable struct MvAdaptiveRWM{M <: AbstractVector{<:Real}, V <: AbstractMatrix{<:Real},
                             B <: Real} <: Metropolis
  sample_mean::M
  sample_cov::V
  beta::B
  iter::Int
  d::Int
end
function MvAdaptiveRWM(v::AbstractVector{<:Real}; beta=0.05, iter=1)
  d = length(v)
  return MvAdaptiveRWM(v, eye(d) * 0.01 / d, beta, iter, d)
end
function update_stats!(rwm::MvAdaptiveRWM, x::AbstractVector{<:Real})
  rwm.iter += 1
  update_mean!(rwm.sample_mean, x, rwm.iter)
  update_cov!(rwm.sample_cov, rwm.sample_mean, x, rwm.iter)
end
function _update(rng::Random.AbstractRNG, rwm::MvAdaptiveRWM, curr::AbstractVector{<:Real},
                 logprob::Function)
  # Update summary stats.
  update_stats!(rwm, curr)

  # Construct proposal covariance.
  proposal_cov = if rwm.iter <= 2 * rwm.d || rwm.beta > rand()
    0.01 * eye(rwm.d) / rwm.d
  else
    5.6644 * rwm.sample_cov / rwm.d
  end
  proposal_cov .= Matrix(LinearAlgebra.Symmetric(proposal_cov))

  return _update(rng, StaticRWM(MvNormal(proposal_cov)), curr, logprob)
end
const mvarwm = MvAdaptiveRWM


text:
module MarkElide
using Cassette
export @mark
export @elide

ismark_begin(x) = false
ismark_end(x) = false
function ismark_begin(ex::Expr)
    ex.head == :meta &&
    length(ex.args) == 3 &&
    ex.args[1] == :mark &&
    ex.args[3] == :begin
end
function ismark_end(ex::Expr)
    ex.head == :meta &&
    length(ex.args) == 3 &&
    ex.args[1] == :mark &&
    ex.args[3] == :end
end

function mark_begin(s::Symbol)
    Expr(:meta, :mark, s, :begin)
end
function mark_end(s::Symbol)
    Expr(:meta, :mark, s, :end)
end

Cassette.@context MarkElideCtx
function Cassette.getpass(::Type{MarkElideCtx})
    function inner(signature, method_body)::CodeInfo
        ret = []
        elide = 0
        for ex in method_body.code
            if ismark_begin(ex)
                elide += 1
            elseif ismark_end(ex)
                @assert elide > 0
                elide -= 1
            elseif elide > 0
                nothing
            else
                push!(ret, ex)
            end
        end
        @assert elide == 0
        method_body = deepcopy(method_body)
        method_body.code = code
        method_body
    end
end

macro mark(symbol, code)
    @assert symbol isa QuoteNode
    @assert symbol.value isa Symbol
    Expr(:block,
         mark_begin(symbol.value),
         Expr(:(=), :ret, esc(code)),
         mark_end(symbol.value),
         :ret)
end

macro elide(symbol, code)
    @assert symbol isa QuoteNode
    @assert symbol.value isa Symbol
    @assert code isa Expr
    @assert code.head == :call
    f = esc(code.args[1])
    args = map(esc, code.args[2:end])
    :(Cassette.overdub(MarkElideCtx,$f)($(args...)))
end # module
end


text:
using BinaryBuilder

name = "Darknet"
version = v"1.0.0"

# Collection of sources required to build HDF5.  Use the CMake download because
# it includes appropriate zlib and szip sources, letting us build them into libhdf5.
sources = [
    "https://github.com/AlexeyAB/darknet/archive/88cccfcad4f9591a429c1e71c88a42e0e81a5e80.zip" =>
    "ee54ccb1a648f0ae4cfa41c6f6e312c5d4f59bb907c985e705c7eae06f81f071"
]

GPU = true
CUDNN = true
CUDNN_HALF = false
OPENCV = false
DEBUG = false
OPENMP = false
LIBSO = true
ZED_CAMERA = false

# ./configure --prefix=$prefix --host=$target
# make -j${nproc}
# make install

# sed -i 's/GPU=0/GPU=1/g' Makefile
# sed -i 's/CUDNN=0/CUDNN=1/g' Makefile
# sed -i 's/CUDNN_HALF=0/CUDNN_HALF=1/g' Makefile
# sed -i 's/OPENCV=0/OPENCV=1/g' Makefile
# sed -i 's/DEBUG=0/DEBUG=1/g' Makefile
# sed -i 's/OPENMP=0/OPENMP=1/g' Makefile
# sed -i 's/LIBSO=0/LIBSO=1/g' Makefile
# sed -i 's/ZED_CAMERA=0/ZED_CAMERA=1/g' Makefile

# cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/darknet-*

./configure --prefix=$prefix --host=$target
make -j${nproc}
make install
"""

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

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

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

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

max_stars_repo_name: ianshmean/DarknetBuilder

max_issues_repo_name: ianshmean/DarknetBuilder

max_forks_repo_name: ianshmean/DarknetBuilder


text: outdirpath = "/Users/frankvandermeulen/Sync/DOCUMENTS/onderzoek/BL_output"