text:
immutable Exponential{T} <: Model
    a::T
    δ::T
end

Base.show(io::IO, ::Type{Exponential}) = print(io, "Classical Exponential")

function Base.show(io::IO, m::Exponential)
    @printf(
        io,
        "Classical Exponential (a = %.3f, δ = %.3f)",
        m.a,
        m.δ,
    )
end

function initial_parameters(::Type{Exponential})
    return Float64[log(1.0), logit(0.5)]
end

canonical_parameters(m::Exponential) = (m.a, m.δ)

function Exponential(parameters::Vector)
    a = exp(parameters[1])
    δ = invlogit(parameters[2])
    return Exponential(a, δ)
end

@inline function predict(
    m::Exponential,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, δ = m.a, m.δ
    u1 = x1 * δ^t1
    u2 = x2 * δ^t2
    z = a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::Exponential,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    a, δ = m.a, m.δ
    u2 = x2 * δ^t2
    u1 = x1 * δ^t1
    gr = (
        (u2 - u1) * a,
        a * ((x2 * t2 * δ^t2 - x1 * t1 * δ^t1) / δ) * δ * (1 - δ),
    )
    return gr
end


text:
module Crypto

include("crypto_common_h.jl")

@ctypedef EVP_MD Void
@ctypedef EVP_MD_CTX Void
@ctypedef ENGINE Void

typealias size_t Csize_t

include("crypto_evp_h.jl")
include("crypto_hmac_h.jl")
include("crypto_md5_h.jl")
include("crypto_sha_h.jl")
include("crypto_defines_h.jl")

include("crypto_exports_h.jl")

hmacsha256_digest(s::String, k::Union(String, Vector{Uint8})) =  hmacsha_digest(s, k, EVP_sha256(), 32)
export hmacsha256_digest

hmacsha1_digest(s::String, k::Union(String, Vector{Uint8})) = hmacsha_digest(s, k, EVP_sha1(), 20)
export hmacsha1_digest

function hmacsha_digest(s::String, k::Union(String, Vector{Uint8}), evp, dgst_len)
    if evp == C_NULL error("EVP_sha1() failed!") end

    sig = zeros(Uint8, dgst_len)
    sig_len = zeros(Uint32, 1)

    if isa(k, String)
        k = convert(Array{Uint8}, k)
    end

    if HMAC(evp, k, length(k), s, length(s), sig, sig_len) == C_NULL error("HMAC() failed!") end
    if (sig_len[1] != dgst_len) error("Wrong length of signature!") end

    return sig
end

function md5(s::String)
    md = zeros(Uint8, 16)
    assert (MD5(s, length(s), md) != C_NULL)
    return md
end

function md5(s::IO)
    evp_md_ctx = EVP_MD_CTX_create()
    assert (evp_md_ctx != C_NULL)

    md = zeros(Uint8, 16)
    try
        evp_md = EVP_md5()
        assert (evp_md != C_NULL)

        rc = EVP_DigestInit_ex(evp_md_ctx, evp_md, C_NULL)
        assert(rc == 1)

        while (!eof(s))
            b = read(s, Uint8, min(nb_available(s), 65536))    # Read in 64 K chunks....

            rc = EVP_DigestUpdate(evp_md_ctx, b, length(b));
            assert(rc == 1)
        end

        rc = EVP_DigestFinal_ex(evp_md_ctx, md, C_NULL)
        assert(rc == 1)

    finally
        EVP_MD_CTX_destroy(evp_md_ctx)
    end

    return md
end

end # Module end

ext: jl


text:
using OrdinaryDiffEq
using Debugger
f(u,p,t) = 1.01*u
u0=1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
@bp
sol = solve(prob,Midpoint(),reltol=1e-8,abstol=1e-8)


text:
#   simple QWTwPlot package example
#   for the introduction please see readme https://github.com/ig-or/QWTWPlot.jl
#

using QWTWPlot
using Random # just for data generation

qstart() # sorry, have to call this explicity here, not in __init__
#qstart(debug=true, qwtw_test = true, libraryName = "qwtw")

# draw thin blue 'sinus':
tMax = 10.0 # let it be `maximum time`
t = Array(range(0.0,stop=tMax, length=10000)); # our `time info`
n = length(t);
y = sin.(t);     # we will draw this signal
qfigure(1); # "1" the number of the first plot window;
# actually we can skip this number, if just creating new plot window (like qfigure() )

 # parameters: 'x' and 'y' data vectors,
 # then 'name of this line',
 #then 'style description', then 'line width'
qplot(t, y, "blue line", "-b", 1)

# add green thick line (another sinus) on the same plot:
y = sin.(t .* 4.0) + cos.(t * 0.8);
qplot(t, y, "thick green line", "-g", 4) #'x' and 'y' data vectors, then 'name of this line', then 'style description', then 'line width'
qtitle("first plot window") # add a title for the first plot
qxlabel("time (in seconds)") # put a label on X axis
qylabel("happiness") # put a label on Y axis

#= By default, "ZOOM" mode is active. This means that
	now you can select some part of the plot with left mouse button
	(press left mouse button, select what you need on a plot,
	release left mouse button).
	Right mouse button will return you back to the previous `ZOOM` state.
	And try how mouse wheel is working also (should work for ZOOM in different modes).

	Button with a "hand" - it's a "pan mode".  It will shift all the picture,
	but it will not change scale.

	"disk" button - allow you to save image to a (png) file

	"[]" buttom means "make square axis" i.e. make equal scale for X and Y axis.
	this `[]` button is useful when you draw something like a `top view plot`

	Real magic is in left and right buttons; see below about it.
    and also look at readme https://github.com/ig-or/QWTWPlot.jl
=#

# create another plot, with high frequency signal and noise
noise = rand(n);
y = sin.(t * 100.0) + noise;
qfigure(2)   # make a plot with another plot ID
qplot(t, y, "sinus + noise", "-m", 2)
qtitle("frequency test")

#= what is the frequency of our signal from last window?
 pless "f" button, another small window will appear;
 then - on new small window - select "sinus + noise",
 select "4" from combo box and close this small window (with Alt-F4 (for example))

 after this, "frequency" plot will be created,
 actually, it is called "Power Spectral Dencity plot"
=#

# add another line to the first plot:
t1 = Array(range(0.,stop = 10., length=5))
y1 = cos.(t1 * 0.5)
qfigure(1); # switch back to the first plot

# if we do not need the lines, only symbols:
qplot1(t1, y1, "points and line #1", " eb",  8)

#= the parameter #4 in qplot or qplot1 functions is a "line style"
	it is a string consisting from 1 or 2 or 3 characters
	the last symbol is a color description;
	first is always a line style:
	middle is symbol type

	for additional info about this line style,
	see info here:
	https://github.com/ig-or/QWTWPlot.jl/blob/master/docs/line-styles.md
=#

 # parameters: 'x' and 'y' data vectors, then 'name of this line',
 # then 'style description', then 'line width', and "symbol size"
qplot(t1, y1, "points and line #2", "-tm", 2, 10)

#= now, try to use an arrow burtton: it will draw a marker on all the plots.
	Press "Arrow" button, then left mouse button somewhere inside the plot.
	Notice that markers will appear on all the plots, and all the lines.

	Next, select some region on one of the plots using ZOOM tool;  and press right "clip" button.
	All the plots will (try to) zoom to the same region.
=#

# create a circle on another plot window:
x = sin.(t .* (2.0*pi/tMax)); y = cos.(t .* (2.0*pi/tMax));
qfigure(3);
qplot2(x, y, t, "circle #1", "-r", 2)
qtitle("circle")

#=  if you do not need _all_ the plots to adjust its range, when you press the 'clip' button,
	then you can use qclipgrp(gr) function. All the plot belong to some 'clip group'.
	after call to 'qclipgrp', new plots will be created with another 'clip group'.
	So, 'clip' button will work only for plots in the same group.

=#
#qclipgrp(1) # next plots will belong to another 'clip group #1'

# now press "[]" buttol in order to make a circle look like circle

# try to use left button ("ARROW") on a circle plot - it also works

# draw one more circle on the same plot:
t1 = Array(range(0.,  stop=2. * pi, length=8));
x1 = 0.5*sin.(t1);
y1 = 0.5*cos.(t1);
qplot2(x1, y1, t1, "circle #2", " ec", 1, 20)

# also, you can remove lines from plots, using values returned from qplot, qplot1 and qplot2 functions.
#  this not working for maps/3D lines (does anybody needs it?)
# lets draw one more circle on the same plot:
id3 = qplot2(x1 .+ 1.0, y1 .+ 0.5, t1, "circle #3", " er", 1, 15)
sleep(1.5)
qremove(id3) # this supposed to remove the last circle
sleep(1.1)

# and it's possible to change existing lines
x8 = collect(0.0:0.1:10)
y8 =  sin.(x8)
qfigure()
id123 = qplot(x8, y8, "changing line", "-b", 2)
qtitle("changing line demo")
for i=1:50
	global y8
	sleep(0.05)
	y8 = sin.(x8 .* (i / 8.))
	qchange(id123, x8, y8)
end

## small 3D example: - disabled (not supported now because of the license limitations)
##qf3d(124)
##qplot3d(x, y, x .* y, "3D", "-g", 1, t)
##qplot3d(x1, y1, x1 .* y1, "3D", "-r", 8, t1)
#
#= show "main window" in order to control all other windows:
it will allow you to switch between all other plot windows

When you have more than 15 plots, it's very convinient to switch
between plots using this helpful window
=#
qsmw()

# another 2000000 points plot:
N = 1000000.0
R1 = 10.0
R2 = 4.0
R3 = 1.2
t = collect(1.0:N) .* (2.0 * pi / N)
x = similar(t)
y = similar(t)
x1 = similar(t)
y1 = similar(t)
@. x = sin(t) * R1 + sin(t * 10.) * R2 + sin(t * 24.) * R3 + rand() * R3*0.25
@. y = cos(t) * R1 + cos(t * 9.) * R2 + cos(t * 21.) * R3 + rand() * R3*0.25
@. x1 = sin(t) * R1 + sin(t * 11.) * R2 + sin(t * 31.) * R3 + randn() * R3*0.015
@. y1 = cos(t) * R1 + cos(t * 8.) * R2 + cos(t * 23.) * R3 + randn() * R3*0.015
qfigure()
qplot2(x, y, t, "$N test 1", " tb", 1, 5)
qplot2(x1, y1, t, "$N test 2", " eg", 1, 3)
qxlabel("?"); qylabel("?"); qtitle("$N points test");

#= qfigure can return a Int32 value, "plot ID", which can be used
 for example
	f1 = qfigure(); qplot(....)
	f2 = qfigure(); qplot(....)
	f3 = qfigure(); qplot(....)
	now we can add something to the plot #2:
	qfigure(f2); qplot(....)
=#

#= log scale supported
=#
x = collect(1.0:0.1:25.)
y = sin.(x)
# this is (default) linear scale:
f1 = qfigure(); qplot1(x, y, "test normal scale", "-eb", 4)
qtitle("normal linear scale")

# x axis log scale:
f2 = qfigure(xAxisType = :aLog); qplot1(x, y, "test x axis log scale", "-eb", 4)
qtitle("x axis log scale")

# both axis log scale:
f3 = qfigure(xAxisType = :aLog, yAxisType = :aLog)
qplot1(x, y .^2 .+ 0.001, "log scale", "-em", 4)
qtitle("x and y axis log scale")

# now lets try to draw a map;
mwn = 4 # for this small example, we have only 4 points
north = [55.688713, 55.698713, 55.678713, 55.60]; # coords in degrees
east = [37.901073, 37.911073, 37.905073, 37.9]; # coords in degrees
t4 = Array(range(0.,stop=tMax, length=4)); # create corresponding time info (for our cool magic markers)

qfmap(5)
qplot2(east, north, t4, "trajectory #1", "-rb",  3, 16);

qtitle("top view test");

#another map:
north1 = [65.688713, 65.698713, 65.678713, 65.60]; # coords in degrees
east1 = [27.901073, 28.111073, 28.005073, 27.9]; # coords in degrees
qfmap(6)
qplot2(east1, north1, t4, "trajectory #2", "-rb",  3);
qplot2(east1, north1, t4, "points", " er",  16);
qtitle("top view test #2");

# ========    draw a few 3D lines ====================
N = 100  # number of points
R = 10.0   # radius
t = [i / N * 2. * pi * 3.5 for i=1:N] # parameter
x = sin.(t) .* R
y = cos.(t) .* R
z = t
z1 = x + y

qmgl() # create a window for 3D
qmgline(x, y, z, "", name = "first 3D line") # draw line #1 with some default style
qmgline(x, y, z1, "-or", name = "second 3D line")  # draw line #2

# ==========  draw a surface =========================

f = [sin((i - 10.0) / 5.0) * sin((j - 10) / 7.0) * 10.0 for i=1:20, j = 1:20]
f1 = [sin((i - 10.0) / 2.0) * sin((j - 10) / 3.0) * 10.0 for i=1:20, j = 1:20]

qmglmesh(f, -10.0, 10., -10., 10.0, type = 1)
qmglmesh(f1, -10.0, 10., -10., 10.0)

# if everything is working as you need, you can
# close all the plots with following command:
#qclear()

# if you need bugs to be fixed, or add some features, do not hesitate to fix them %)
# and create a pull request
# and/or create an issue on github - maybe it'll be
# much faster for me to fix everything or make some  changes or add functionality
#


text:
# Some initial support for Healpix fields.
# The main functionality of broadcasting, indexing, and projection for
# a few field types is implemented, but not much beyond that.

@init global hp = lazy_pyimport("healpy")

struct ProjHealpix <: FieldMetadata
    Nside :: Int
end
typealias_def(::Type{<:ProjHealpix}) = "ProjHealpix"

# spin-0
const HealpixMap{        M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{Map,        M, T, A}
const HealpixFourier{    M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{Fourier,    M, T, A}
# spin-2
const HealpixQUMap{      M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{QUMap,      M, T, A}
const HealpixQUFourier{  M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{QUFourier,  M, T, A}
const HealpixEBMap{      M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{EBMap,      M, T, A}
const HealpixEBFourier{  M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{EBFourier,  M, T, A}
# spin-(0,2)
const HealpixIQUMap{     M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{IQUMap,     M, T, A}
const HealpixIQUFourier{ M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{IQUFourier, M, T, A}
const HealpixIEBMap{     M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{IEBMap,     M, T, A}
const HealpixIEBFourier{ M<:ProjHealpix, T, A<:AbstractArray{T} } = BaseField{IEBFourier, M, T, A}

const HealpixField{B, M<:ProjHealpix, T, A<:AbstractArray{T}} = BaseField{B, M, T, A}

## constructing from arrays
# spin-0
function HealpixMap(I::A) where {T, A<:AbstractArray{T}}
    HealpixMap(I, ProjHealpix(hp.npix2nside(length(I))))
end
# spin-2
function HealpixField{B}(X::A, Y::A) where {T, A<:AbstractArray{T}, B<:Basis2Prod{<:Union{𝐐𝐔,𝐄𝐁},Map}}
    HealpixField{B}(cat(X, Y, dims=Val(2)), ProjHealpix(hp.npix2nside(length(X))))
end
# spin-(0,2)
function HealpixField{B}(I::A, X::A, Y::A) where {T, A<:AbstractArray{T}, B<:Basis3Prod{𝐈,<:Union{𝐐𝐔,𝐄𝐁},Map}}
    HealpixField{B}(cat(I, X, Y, dims=Val(2)), ProjHealpix(hp.npix2nside(length(I))))
end

### pretty printing
typealias_def(::Type{F}) where {B,M<:ProjHealpix,T,A,F<:HealpixField{B,M,T,A}} = "Healpix$(typealias(B)){$(typealias(A)),$(typealias(M))}"
function Base.summary(io::IO, f::HealpixField)
    @unpack Nside = f
    print(io, "$(length(f))-element Nside=$Nside ")
    Base.showarg(io, f, true)
end

# useful for enumerating some cases below
_healpix_sub_components = [
    (:HealpixMap,        ("Ix"=>:, "I"=>:)),
    (:HealpixFourier,    ("Il"=>:, "I"=>:)),
    (:HealpixQUMap,      ("Qx"=>1, "Ux"=>2, "Q" =>1, "U"=>2, "P"=>:)),
    (:HealpixQUFourier,  ("Ql"=>1, "Ul"=>2, "Q" =>1, "U"=>2, "P"=>:)),
    (:HealpixEBMap,      ("Ex"=>1, "Bx"=>2, "E" =>1, "B"=>2, "P"=>:)),
    (:HealpixEBFourier,  ("El"=>1, "Bl"=>2, "E" =>1, "B"=>2, "P"=>:)),
    (:HealpixIQUMap,     ("Ix"=>1, "Qx"=>2, "Ux"=>3, "I"=>1, "Q"=>2, "U"=>3, "P"=>2:3, "IP"=>:)),
    (:HealpixIQUFourier, ("Il"=>1, "Ql"=>2, "Ul"=>3, "I"=>1, "Q"=>2, "U"=>3, "P"=>2:3, "IP"=>:)),
    (:HealpixIEBMap,     ("Ix"=>1, "Ex"=>2, "Bx"=>3, "I"=>1, "E"=>2, "B"=>3, "P"=>2:3, "IP"=>:)),
    (:HealpixIEBFourier, ("Il"=>1, "El"=>2, "Bl"=>3, "I"=>1, "E"=>2, "B"=>3, "P"=>2:3, "IP"=>:)),
]
# sub-components
for (F, props) in _healpix_sub_components
    for (k,I) in props
        body = if k[end] in "xl"
            I==(:) ? :(getfield(f,:arr)) : :(view(getfield(f,:arr), :, $I))
        else
            I==(:) ? :f : :($(HealpixField{k=="P" ? Basis2Prod{basis(@eval($F)).parameters[end-1:end]...} : basis(@eval($F)).parameters[end]})(view(getfield(f,:arr), :, $I), f.metadata))
        end
        @eval getproperty(f::$F, ::Val{$(QuoteNode(Symbol(k)))}) = $body
    end
end
getproperty(f::HealpixField, ::Val{:proj}) = getfield(f,:metadata)

### broadcasting

function promote_metadata_strict(metadata₁::ProjHealpix, metadata₂::ProjHealpix)
    if metadata₁ === metadata₂
        metadata₁
    else
        error("Can't broadcast Healpix maps with two different Nsides ($metadata₁.Nside, $metadata₂.Nside).")
    end
end

### Projection

## Coordinate transformations

# choices of negative signs are such that the default projection is
# centered on the center of a healpy.mollview plot and oriented the
# same way.

function xy_to_θϕ(::ProjLambert, x, y)
    r = sqrt(x^2+y^2)
    θ = 2*acos(r/2)
    ϕ = atan(-x,-y)
    θ, ϕ
end

function θϕ_to_xy(::ProjLambert, θ, ϕ)
    r = 2cos(θ/2)
    x = -r*sin(ϕ)
    y = -r*cos(ϕ)
    x, y
end

# adding more flat projections in the future should be as easy as
# defining pairs of functions like the above two for other cases

function xy_to_θϕ(flat_proj::FlatProj; rotator)
    @unpack Δx, Ny, Nx = flat_proj
    xs = @. Δx * ((-Nx÷2:Nx÷2-1) + 0.5)
    ys = @. Δx * ((-Ny÷2:Ny÷2-1) + 0.5)
    θϕs = xy_to_θϕ.(Ref(flat_proj), xs, ys')
    (eachrow(rotator.get_inverse()(first.(θϕs)[:], last.(θϕs)[:]))...,)
end

function θϕ_to_xy((hpx_proj, flat_proj)::Pair{<:ProjHealpix, <:FlatProj}; rotator)
    @unpack Δx, Ny, Nx = flat_proj
    @unpack Nside = hpx_proj
    (θs, ϕs) = hp.pix2ang(Nside, 0:(12*Nside^2-1))
    (θs′, ϕs′) = eachrow(rotator(θs, ϕs))
    xys = θϕ_to_xy.(Ref(flat_proj), θs′, ϕs′)
    xs = @. first(xys) / Δx + Nx÷2 + 0.5
    ys = @.  last(xys) / Δx + Ny÷2 + 0.5
    (xs, ys)
end

## Healpix => Flat

"""
    project(healpix_map::HealpixField => flat_proj::FlatProj; [rotator])

Project `healpix_map` to a flat projection specified by `flat_proj`.
E.g.:

    healpix_map = HealpixMap(rand(12*2048^2))
    flat_proj = ProjLambert(Ny=128, Nx=128, θpix=3, T=Float32)
    f = project(healpix_map => flat_proj; rotator=pyimport("healpy").Rotator((0,28,23)))

The use of `=>` is to help remember in which order the arguments are
specified. Optional keyword argument `rotator` is a `healpy.Rotator`
object specifying a rotation which rotates the north pole to the
center of the desired field. If `healpix_map` is a `HealpixQU`, Q/U
polarization angles are rotated to be aligned with the local
coordinates (sometimes called "polarization flattening").
"""
function project((hpx_map, flat_proj)::Pair{<:HealpixMap,<:ProjLambert{T}}; rotator=hp.Rotator((0,90,0))) where {T}
    _project(xy_to_θϕ(flat_proj; rotator), hpx_map => flat_proj)
end

function project((hpx_map, flat_proj)::Pair{<:HealpixQUMap,<:ProjLambert{T}}; rotator=hp.Rotator((0,90,0))) where {T}
    @unpack Ny, Nx = flat_proj
    (θs, ϕs) = xy_to_θϕ(flat_proj; rotator)
    Q = _project((θs, ϕs), hpx_map.Q => flat_proj)
    U = _project((θs, ϕs), hpx_map.U => flat_proj)
    ψpol = permutedims(reshape((hp.Rotator((0,-90,0)) * rotator).angle_ref(θs, ϕs), Nx, Ny))
    QU_flat = @. (Q.arr + im * U.arr) * exp(im * 2 * ψpol)
    FlatQUMap(real.(QU_flat), imag.(QU_flat), flat_proj)
end

function _project((θs,ϕs)::Tuple, (hpx_map, flat_proj)::Pair{<:HealpixMap,<:ProjLambert{T}}) where {T}
    @unpack Ny, Nx = flat_proj
    FlatMap(T.(permutedims(reshape(hp.get_interp_val(collect(hpx_map), θs, ϕs), Nx, Ny))), flat_proj)
end

## Flat => Healpix

"""
    project(flat_map::FlatField => healpix_proj::ProjHealpix; [rotator])

Reproject a `flat_map` back onto the sphere in a Healpix projection
specified by `healpix_proj`. E.g.

    flat_map = FlatMap(rand(128,128))
    f = project(flat_map => ProjHealpix(2048); rotator=pyimport("healpy").Rotator((0,28,23)))

The use of `=>` is to help remember in which order the arguments are
specified. Optional keyword argument `rotator` is a `healpy.Rotator`
object specifying a rotation which rotates the north pole to the
center of the desired field.
"""
function project((flat_map, hpx_proj)::Pair{<:FlatS0, <:ProjHealpix}; rotator=hp.Rotator((0,90,0)))
    (xs, ys) = θϕ_to_xy(hpx_proj => flat_map.proj; rotator)
    HealpixMap(@ondemand(Images.bilinear_interpolation).(Ref(cpu(flat_map[:Ix])), ys, xs), hpx_proj)
end


text:
# Gaussian
using Distributions
using LinearAlgebra: norm_sqr

import Base: rand
import Distributions: pdf, logpdf, sqmahal, cdf, quantile
import Base: size

"""
    PSD{T}

Simple wrapper for the lower triangular Cholesky root of a positive (semi-)definite element `σ`.
"""
mutable struct PSD{T}
    σ::T
    PSD(σ::T) where {T} = istril(σ) ? new{T}(σ) : throw(ArgumentError("Argument not lower triangular"))
end
cholupper(P::PSD) = P.σ'

sumlogdiag(Σ::Float64, d=1) = log(Σ)
sumlogdiag(Σ,d) = sum(log.(diag(Σ)))
sumlogdiag(J::UniformScaling, d)= log(J.λ)*d

_logdet(Σ::PSD, d) = 2*sumlogdiag(Σ.σ, d)

_logdet(Σ, d) = logdet(Σ)
_logdet(J::UniformScaling, d) = log(J.λ) * d

_symmetric(Σ) = Symmetric(Σ)
_symmetric(J::UniformScaling) = J

"""
    Gaussian(μ, Σ) -> P

Gaussian distribution with mean `μ`` and covariance `Σ`. Defines `rand(P)` and `(log-)pdf(P, x)`.
Designed to work with `Number`s, `UniformScaling`s, `StaticArrays` and `PSD`-matrices.

Implementation details: On `Σ` the functions `logdet`, `whiten` and `unwhiten`
(or `cholupper` as fallback for the latter two) are called.
"""
struct Gaussian{T,S}
    μ::T
    Σ::S
    Gaussian(μ::T, Σ::S) where {T,S} = new{T,S}(μ, Σ)
end
Gaussian() = Gaussian(0.0, 1.0)

dim(P::Gaussian) = length(P.μ)
whiten(Σ::PSD, z) = Σ.σ\z
whiten(Σ, z) = cholupper(Σ)'\z
whiten(Σ::UniformScaling, z) = z/sqrt(Σ.λ)
sqmahal(P::Gaussian, x) = norm_sqr(whiten(P.Σ, x - P.μ))

rand(P::Gaussian) = P.μ + cholupper(P.Σ)'*randn(typeof(P.μ))
rand(P::Gaussian{Vector}) = P.μ + cholupper(P.Σ)'*randn(eltype(P.μ), length(P.μ))

logpdf(P::Gaussian, x) = -(sqmahal(P,x) + _logdet(P.Σ, dim(P)) + dim(P)*log(2pi))/2
pdf(P::Gaussian, x) = exp(logpdf(P::Gaussian, x))
cdf(P::Gaussian{Float64,Float64}, x) = Distributions.normcdf(P.μ, sqrt(P.\Sigma), x)

"""
    logpdfnormal(x, Σ)

logpdf of centered Gaussian with covariance Σ
"""
function logpdfnormal(x, Σ)

    S = cholupper(_symmetric(Σ))'

    d = length(x)
     -((norm(S\x))^2 + 2sumlogdiag(S,d) + d*log(2pi))/2
end
function logpdfnormal(x::Float64, Σ)
     -(x^2/Σ + log(Σ) + log(2pi))/2
end

function conditional(P::Gaussian, A, B, xB)
    assert(length(b) == length(xB))
    Z = P.Σ[A,B]*inv(P.Σ[B,B])
    Gaussian(P.μ[A] + Z*(xB - P.μ[B]), P.Σ[A,A] - Z*P.Σ[B,A])
end


text:
#=
    Scenario reduction functions
=#
abstract type AbstractScenarios end

mutable struct Scenarios{T, O, I} <: AbstractScenarios
    demands::Vector{NamedTuple{(:t, :power),Tuple{T,O}}}
    generations::Vector{NamedTuple{(:t, :power, :cost), Tuple{T, O, I}}}
    storages::Vector{NamedTuple{(:cost,), Tuple{I}}}
    converters::Vector{NamedTuple{(:cost,), Tuple{I}}}
    grids::Vector{NamedTuple{(:cost_in, :cost_out), Tuple{O, O}}}
end

# Constructor
function Scenarios(mg::Microgrid, d::Dict{})
    # Utils to simplify the writting
    h, y, s = 1:mg.parameters.nh, 1:mg.parameters.ny, 1:mg.parameters.ns
    T, O, I = Array{DateTime,3}, Array{Float64, 3}, Array{Float64, 2}
    # Initialize
    demands = Vector{NamedTuple{(:t, :power),Tuple{T,O}}}(undef, length(mg.demands))
    generations = Vector{NamedTuple{(:t, :power, :cost), Tuple{T, O, I}}}(undef, length(mg.generations))
    storages = Vector{NamedTuple{(:cost,), Tuple{I}}}(undef, length(mg.storages))
    converters = Vector{NamedTuple{(:cost,), Tuple{I}}}(undef, length(mg.converters))
    grids = Vector{NamedTuple{(:cost_in, :cost_out), Tuple{O, O}}}(undef, length(mg.grids))
    # Demands
    for (k, a) in enumerate(mg.demands)
        if a.carrier isa Electricity
            demands[k] = (t = d["ld_E"].t[h, y, s], power = d["ld_E"].power[h, y, s])
        elseif a.carrier isa Heat
            demands[k] = (t = d["ld_H"].t[h, y, s], power = d["ld_H"].power[h, y, s])
        end
    end
    # Generation
    for (k, a) in enumerate(mg.generations)
        if a isa Solar
            generations[k] = (t = d["pv"].t[h, y, s], power = d["pv"].power[h, y, s], cost = d["pv"].cost[y, s])
        end
    end
    # Storages
    for (k, a) in enumerate(mg.storages)
        if a isa Liion
            storages[k] = (cost = d["liion"].cost[y, s],)
        elseif a isa ThermalStorage
            storages[k] = (cost = d["tes"].cost[y, s],)
        elseif a isa H2Tank
            storages[k] = (cost = d["h2tank"].cost[y, s],)
        end
    end
    # Converters
    for (k, a) in enumerate(mg.converters)
        if a isa Electrolyzer
            converters[k] = (cost = d["elyz"].cost[y, s],)
        elseif a isa FuelCell
            converters[k] = (cost = d["fc"].cost[y, s],)
        elseif a isa Heater
            converters[k] = (cost = d["heater"].cost[y, s],)
        end
    end
    # Grids
    for (k, a) in enumerate(mg.grids)
        if a.carrier isa Electricity
            grids[k] = (cost_in = d["grid"].cost_in[h, y, s], cost_out = d["grid"].cost_out[h, y, s])
        end
    end

    return Scenarios(demands, generations, storages, converters, grids)
end


text:
using Ipopt
using PowerModelsPrivacyPreserving
using PowerModelsDistribution
using Distributions
using Random

# const PMs = PowerModels
const PMD = PowerModelsDistribution
const PMPP = PowerModelsPrivacyPreserving

ipopt = Ipopt.Optimizer
file =  "test/data/opendss/case3_unbalanced.dss"
data_unpert = PMD.parse_file(file)

result_unpert_cost_mc = PMPP.run_mc_opf_bf_dvorkin(data_unpert, LinDist3FlowPowerModel, ipopt)

##


text:
using DataFrames, Resampling

## File that runs the entire cross-val analysis and dumps it for the traffic example
type Arc
	initNode::Int
	termNode::Int
    capacity::Float64
    freeflowtime::Float64
    trueflow::Float64
    obsflow::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64,freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0., 0.)

include("fitTraffic.jl")

############
#Read in the demand file
file = open("SiouxFalls_trips.txt")
demands = Dict{(Int64,Int64), Float64}()
s = 0
for line in eachline(file)
    if contains(line, "Origin")
        s = int(split(line)[2])
    else
        pairs = split(line, ";")
        for pair in pairs
            if !contains(pair, "\n")
                pair_vals = split(pair, ":")
                t, demand = int(pair_vals[1]), float(pair_vals[2])
                demands[(s,t)] = demand
            end
        end
    end
end
close(file)

############
#read in the arc files
arcs = Dict{(Int, Int), Arc}()
file = open("SiouxFalls_net.txt")
inHeader=true
for line in eachline(file)
    if inHeader
        inHeader = !contains(line, "Init node")
        continue
    end
    vals = split(line, )
    arcs[(int(vals[1]), int(vals[2]))] = Arc(int(vals[1]), int(vals[2]),
                                              float(vals[3]), float(vals[5]))
end
close(file)

###########
#read in the true costs
file = open("SiouxFallsFlow.txt")
ix = 0; test_cost = 0.0
for line in eachline(file)
    ix +=1
    if ix ==1
        continue
    end
    vals = split(line)
    arcs[(int(vals[1]), int(vals[2]))].trueflow = float(vals[3])
    test_cost += float(vals[4])
end
tot_cost = sum(map(a -> bpacost(a.trueflow, a), values(arcs) ))
assert( abs( tot_cost - test_cost) <= 1e-8 ) #note, these costs are not scaled by x!!!

##########
#Generate the simulated data
##########
numData = 20 ; sigma = .1
flow_data = Array(Float64, length(arcs), numData)
demand_data = Dict{(Int, Int), Array{Float64, 1}}()

numNodes = maximum(map(pair->pair[1], keys(demands)))
g = simple_inclist(numNodes, is_directed=true)
vArcs = Arc[]
for arc in values(arcs)
    arc.obsflow = arc.trueflow
    add_edge!(g, arc.initNode, arc.termNode)
    push!(vArcs, arc)
end

srand(8675309)
for iRun = 1:numData
    #perturb the demand_data
    for odpair in keys(demands)
        if ! haskey(demand_data, odpair)
            demand_data[odpair] = [demands[odpair] * (1 + sigma * rand()), ]
        else
            push!(demand_data[odpair], demands[odpair] * (1  + sigma * rand()))
        end
    end
    #solve using FW and record
    conv_tol, flow_data[:, iRun] = frank_wolfe(g, vArcs, demand_data, iRun)
    println(conv_tol)
end

#Randomize the flow data as well
for i = 1:size(flow_data, 1)
    for j = 1:size(flow_data, 2)
        flow_data[i, j] *= (1 + sigma * rand() )
    end
end

##########
#Fitting Funcs
##########
#build a little train function that just takes indices
function train(indices, lam::Float64, deg::Int, c::Float64,
                demand_data, flow_data, arcs; fcoeffs=nothing)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))
    m, coeffs, reg_term = setUpFitting(deg, c, [k for k=keys(demand_data)], 1:numNodes)

    for a in values(arcs)
        a.obsflow = a.trueflow
    end
    addIncreasingCnsts(m, coeffs, arcs, TOL=1e-8)  #uses the original obs flows

    avgCost = mean( [bpacost(a.obsflow, a.capacity, 1.0) for a in values(arcs)] )
    normalize(m, coeffs, [a.obsflow / a.capacity for a in values(arcs)], avgCost)

    # normalize(m, coeffs, 1., bpacost(1., 1., 1.))

    resids = Variable[]
    for i = indices
        #copy the flow data over to the arcs, demand data to demands (slow)
        for (ix, a) in enumerate(vArcs)
            a.obsflow = flow_data[ix, i]
        end
        for odpair in keys(demands)
            demands[odpair] = demand_data[odpair][i]
        end

        #Dual Feasibility
        ys = addNetworkCnsts(m, coeffs, demands, arcs, numNodes)

        #add the residual for this data point
        push!(resids, addResid(m, coeffs, ys, demands, arcs, 1e6))
    end

    if fcoeffs != nothing
        fixCoeffs(m, fcoeffs, coeffs)
    end
    @setObjective(m, Min, sum{resids[i], i = 1:length(resids)} + lam*reg_term)
    solve(m)
    println(getObjectiveValue(m) - lam * getValue(reg_term) )
    return [getValue(coeffs[i]) for i =1:length(coeffs)]
end

# a function for testing.
function test(fcoeffs, indices, demand_data, flow_data, arcs, g, vArcs)
    numNodes = maximum(map(pair->pair[1], keys(arcs)))

    l1_eps = 0.
    for i = indices
        costs = [polyEval(fcoeffs, flow_data[ix, i]/a.capacity)*a.freeflowtime
                    for (ix, a) in enumerate(vArcs)]
        primal = dot(costs, flow_data[:, i])

        #solve a bunch of shortest path problems to get dual cost
        dual = 0.
        for odpair = keys(demand_data)
            #solve the shortest path problems, and update the total flow
            r = dijkstra_shortest_paths(g, costs, odpair[1] )
            dual += r.dists[odpair[2]] * demand_data[odpair][i]
        end
        l1_eps += abs(primal - dual)
    end
    return l1_eps / length(indices)
end

df = DataFrame()
df["A"] = 1:numData

lamb_grid = [10. .^(-7:1)]
c_grid = linspace(1, 5, 6)  # This choice of c_grid yields no distinguishable difference. Try: c_grid = 2. .^(1:5)
deg_grid = [2:6]
N = length(lamb_grid) * length(c_grid) * length(deg_grid)

res = Array(Float64, N, 5)
ix = 1
for l in lamb_grid
    for c in c_grid
        for d in deg_grid
            train(df) = train(df["A"], l, d, c, demand_data, flow_data, arcs)
            test(df, fit) = test(fit, df["A"], demand_data, flow_data, arcs, g, vArcs)
            rtrain, rtest = kfold_crossvalidate(df, train, test, 5)

            res[ix, 1] = l
            res[ix, 2] = c
            res[ix, 3] = float(d)
            res[ix, 4] = mean(rtest) / 1e6
            res[ix, 5] = std(rtest) / 1e6
            ix +=1

            show(res)
        end
    end
end

writetable("trafficCVal.csv", DataFrame(res))

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs


text:
# HTS
# ===
#
# High-throughtput sequencing.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

include("sam/sam.jl")
include("bam/bam.jl")

function Base.show(io::IO, rec::Union{SAMRecord,BAMRecord})
    name = seqname(rec)
    if isempty(name)
        name = "(empty name)"
    end
    if get(io, :compact, false)
        if ismapped(rec)
            range = string(refname(rec), ':', leftposition(rec), '-', rightposition(rec))
            print(io, name, "\t$(range)\t", cigar(rec))
        else
            print(io, name, "\tunmapped")
        end
    else
        println(io, summary(rec), ':')
        println(io, "  sequence name: ", name)
        println(io, "  reference name: ", refname(rec))
        println(io, "  leftmost position: ", leftposition(rec))
        println(io, "  next reference name: ", nextrefname(rec))
        println(io, "  next leftmost position: ", nextleftposition(rec))
        println(io, "  mapping quality: ", mappingquality(rec))
        println(io, "  flag: ", flag(rec))
        println(io, "  template length: ", templatelength(rec))
        println(io, "  CIGAR string: ", cigar(rec))
        println(io, "  sequence: ", String(sequence(rec)))
        println(io, "  base qualities: ", qualities(rec))
          print(io, "  optional fields: ", optional_fields(rec))
    end
end


text:
module TestArrayCommon

using Base.Test
using CategoricalArrays
using NullableArrays
using CategoricalArrays: DefaultRefType, index
using Compat

# == currently throws an error for Nullables
(==) = isequal

# Test that mergelevels handles mutually compatible orderings
@test CategoricalArrays.mergelevels(true, [6, 2, 4, 8], [2, 3, 5, 4], [2, 4, 8]) ==
    ([6, 2, 3, 5, 4, 8], true)
@test CategoricalArrays.mergelevels(true, ["A", "B", "D"], ["C", "D", "E", "F"], ["B", "C", "E"]) ==
    (["A", "B", "C", "D", "E", "F"], true)
@test CategoricalArrays.mergelevels(true, ["A", "B", "D"], ["C", "D", "E", "F"], ["A", "B", "C", "F"]) ==
    (["A", "B", "C", "D", "E", "F"], true)
@test CategoricalArrays.mergelevels(true, ["B", "C", "D"], ["A", "B"]) ==
    (["A", "B", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, ["B", "C", "D"], []) ==
    (["B", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, [], ["A", "C"]) ==
    (["A", "C"], true)
@test CategoricalArrays.mergelevels(true, [], []) ==
    ([], true)
@test CategoricalArrays.mergelevels(true, ["A", "B", "C"], [], ["A", "C", "D"]) ==
    (["A", "B", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, [], ["A", "B", "C"], ["A", "C", "D"]) ==
    (["A", "B", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, ["A", "C", "D"], ["A", "B", "C"], []) ==
    (["A", "B", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, [], ["A", "C", "D"], []) ==
    (["A", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, ["A", "C", "D"], [], []) ==
    (["A", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, [], [], ["A", "C", "D"]) ==
    (["A", "C", "D"], true)
@test CategoricalArrays.mergelevels(true, [], [], []) ==
    ([], true)
@test CategoricalArrays.mergelevels(true, ["A", "C", "D"], ["A", "B", "C"], ["A", "D", "E"], ["C", "D"]) ==
    (["A", "B", "C", "D", "E"], true)

# Test that mergelevels handles mutually incompatible orderings
@test CategoricalArrays.mergelevels(true, [6, 3, 4, 7], [2, 3, 6, 5, 4], [2, 4, 8]) ==
    ([6, 2, 3, 5, 4, 7, 8], false)
@test CategoricalArrays.mergelevels(true, ["A", "C", "D"], ["D", "C"], []) ==
    (["A", "C", "D"], false)
@test CategoricalArrays.mergelevels(true, ["A", "D", "C"], ["A", "B", "C"], ["A", "D", "E"], ["C", "D"]) ==
    (["A", "D", "B", "C", "E"], false)

# Test that mergelevels handles incomplete orderings
@test CategoricalArrays.mergelevels(true, ["B", "C", "D"], ["A", "C"]) ==
    (["B", "A", "C", "D"], false)
@test CategoricalArrays.mergelevels(true, ["A", "B", "D"], ["C", "D", "E", "F"]) ==
    (["A", "B", "C", "D", "E", "F"], false)
@test CategoricalArrays.mergelevels(true, ["A", "B", "E", "G"], ["C", "D", "E", "F"]) ==
    (["A", "B", "C", "D", "E", "G", "F"], false)
@test CategoricalArrays.mergelevels(true, ["A", "B", "D"], ["C", "D", "E", "F"], ["A", "B", "D", "F"]) ==
    (["A", "B", "C", "D", "E", "F"], false)

# Test with ordered=false (simpler, almost a subset of code paths tested above)
@test CategoricalArrays.mergelevels(false, ["A", "B", "E", "G"], ["C", "D", "E", "F"]) ==
    (["A", "B", "C", "D", "E", "G", "F"], false)

for (CA, A) in ((CategoricalArray, Array), (NullableCategoricalArray, NullableArray))
    # Test vcat()

    # Test that vcat of compress arrays uses a reftype that doesn't overflow
    a1 = 3:200
    a2 = 300:-1:100
    ca1 = CA(a1)
    ca2 = CA(a2)
    cca1 = compress(ca1)
    cca2 = compress(ca2)
    r = vcat(cca1, cca2)
    @test r == A(vcat(a1, a2))
    @test isa(cca1, CA{Int, 1, UInt8})
    @test isa(cca2, CA{Int, 1, UInt8})
    @test isa(r, CA{Int, 1, CategoricalArrays.DefaultRefType})
    @test isa(vcat(cca1, ca2), CA{Int, 1, CategoricalArrays.DefaultRefType})
    @test isordered(r) == false
    @test levels(r) == collect(3:300)

    # Test vcat of multidimensional arrays
    a1 = Array{Int}(2, 3, 4, 5)
    a2 = Array{Int}(3, 3, 4, 5)
    a1[1:end] = (length(a1):-1:1) + 2
    a2[1:end] = (1:length(a2)) + 10
    ca1 = CA(a1)
    ca2 = CA(a2)
    cca1 = compress(ca1)
    cca2 = compress(ca2)
    r = vcat(cca1, cca2)
    @test r == A(vcat(a1, a2))
    @test isa(r, CA{Int, 4, CategoricalArrays.DefaultRefType})
    @test isordered(r) == false
    @test levels(r) == collect(3:length(a2)+10)

    # Test concatenation of mutually compatible levels
    a1 = ["Young", "Middle"]
    a2 = ["Middle", "Old"]
    ca1 = CA(a1, ordered=true)
    ca2 = CA(a2, ordered=true)
    @test levels!(ca1, ["Young", "Middle"]) === ca1
    @test levels!(ca2, ["Middle", "Old"]) === ca2
    r = vcat(ca1, ca2)
    @test r == A(vcat(a1, a2))
    @test levels(r) == ["Young", "Middle", "Old"]
    @test isordered(r) == true

    # Test concatenation of conflicting ordering. This drops the ordering
    a1 = ["Old", "Young", "Young"]
    a2 = ["Old", "Young", "Middle", "Young"]
    ca1 = CA(a1, ordered=true)
    ca2 = CA(a2, ordered=true)
    levels!(ca1, ["Young", "Middle", "Old"])
    # ca2 has another order
    r = vcat(ca1, ca2)
    @test r == A(vcat(a1, a2))
    @test levels(r) == ["Young", "Middle", "Old"]
    @test isordered(r) == false

    # Test similar()

    x = CA(["Old", "Young", "Middle", "Young"])
    y = similar(x)
    @test typeof(x) === typeof(y)
    @test size(y) == size(x)

    x = CA(["Old", "Young", "Middle", "Young"])
    y = similar(x, 3)
    @test typeof(x) === typeof(y)
    @test size(y) == (3,)

    x = CA{String, 1, UInt8}(["Old", "Young", "Middle", "Young"])
    y = similar(x, Int)
    @test isa(y, CA{Int, 1, UInt8})
    @test size(y) == size(x)

    x = CA(["Old", "Young", "Middle", "Young"])
    y = similar(x, Int, 3, 2)
    @test isa(y, CA{Int, 2, CategoricalArrays.DefaultRefType})
    @test size(y) == (3, 2)

    # Test copy!()

    x = CA(["Old", "Young", "Middle", "Young"])
    levels!(x, ["Young", "Middle", "Old"])
    ordered!(x, true)
    y = CA(["X", "Z", "Y", "X"])
    @test copy!(x, y) === x
    @test x == y
    @test levels(x) == ["Young", "Middle", "Old", "X", "Y", "Z"]
    @test !isordered(x)

    x = CA(["Old", "Young", "Middle", "Young"])
    levels!(x, ["Young", "Middle", "Old"])
    ordered!(x, true)
    y = CA(["X", "Z", "Y", "X"])
    a = A(["Z", "Y", "X", "Young"])
    # Test with null values
    if CA === NullableCategoricalArray
        x[3] = Nullable()
        y[3] = a[2] = Nullable()
    end
    @test copy!(x, 1, y, 2) === x
    @test x == a
    @test levels(x) == ["Young", "Middle", "Old", "X", "Y", "Z"]
    @test !isordered(x)

    # Test that no corruption happens in case of bounds error
    @test_throws BoundsError copy!(x, 10, y, 2)
    @test x == a
    @test_throws BoundsError copy!(x, 1, y, 10)
    @test x == a
    @test_throws BoundsError copy!(x, 10, y, 20)
    @test x == a
    @test_throws BoundsError copy!(x, 10, y, 2)
    @test x == a
    @test_throws BoundsError copy!(x, 1, y, 2, 10)
    @test x == a
    @test_throws BoundsError copy!(x, 4, y, 1, 2)
    @test x == a
    @test_throws BoundsError copy!(x, 1, y, 4, 2)
    @test x == a

    # Test resize!()
    x = CA(["Old", "Young", "Middle", "Young"])
    @test resize!(x, 3) === x
    @test x == A(["Old", "Young", "Middle"])
    @test resize!(x, 4) === x
    if CA === NullableCategoricalArray
        @test x == A(["Old", "Young", "Middle", Nullable()])
    else
        @test x[1:3] == A(["Old", "Young", "Middle"])
        @test !isassigned(x, 4)
    end

    for (sstart, dstart, n) in ((1, 1, 4), (1, 2, 3))
        # Conflicting orders: check that the destination wins and that result is not ordered
        x = CA(["Old", "Young", "Middle", "Young"])
        levels!(x, ["Young", "Middle", "Old"])
        ordered!(x, true)
        y = CA(["Middle", "Middle", "Old", "Young"])
        levels!(y, ["Old", "Middle", "Young"])
        ordered!(x, true)
        @test copy!(x, y) === x
        @test levels(x) == ["Young", "Middle", "Old"]
        @test !isordered(x)

        # Destination ordered, but not origin: check that destination wins
        x = CA(["Old", "Young", "Middle", "Young"])
        levels!(x, ["Young", "Middle", "Old"])
        ordered!(x, true)
        y = CA(["Middle", "Middle", "Old", "Young"])
        levels!(y, ["Young", "Middle", "Old"])
        @test copy!(x, y) === x
        @test levels(x) == ["Young", "Middle", "Old"]
        @test isordered(x)

        # Origin ordered, but not destination: check that destination wins
        x = CA(["Old", "Young", "Middle", "Young"])
        levels!(x, ["Young", "Middle", "Old"])
        y = CA(["Middle", "Middle", "Old", "Young"])
        ordered!(y, true)
        levels!(y, ["Young", "Middle", "Old"])
        @test copy!(x, y) === x
        @test levels(x) == ["Young", "Middle", "Old"]
        @test !isordered(x)

        # Origin ordered, destination ordered with no levels: check that result is ordered
        y = CA(["Middle", "Middle", "Old", "Young"])
        ordered!(y, true)
        levels!(y, ["Young", "Middle", "Old"])
        x = similar(x)
        ordered!(x, true)
        @test copy!(x, y) === x
        @test levels(x) == ["Young", "Middle", "Old"]
        @test isordered(x)

        # Destination ordered, but not origin, and new levels: check that result is unordered
        x = CA(["Old", "Young", "Middle", "Young"])
        levels!(x, ["Young", "Middle", "Old"])
        ordered!(x, true)
        y = CA(["Middle", "Middle", "Old", "Young"])
        levels!(y, ["X", "Young", "Middle", "Old"])
        @test copy!(x, y) === x
        @test levels(x) == ["X", "Young", "Middle", "Old"]
        @test !isordered(x)
    end

    # Test that overflow of reftype is detected and doesn't corrupt data and levels

    res = @test_throws LevelsException{Int, UInt8} CA{Int, 1, UInt8}(256:-1:1)
    @test res.value.levels == [1]
    @test sprint(showerror, res.value) == "cannot store level(s) 1 since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."

    x = CA{Int, 1, UInt8}(254:-1:1)
    x[1] = 1000
    res = @test_throws LevelsException{Int, UInt8} x[1] = 1001
    @test res.value.levels == [1001]
    @test sprint(showerror, res.value) == "cannot store level(s) 1001 since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."
    @test x == A(vcat(1000, 253:-1:1))
    @test levels(x) == vcat(1:254, 1000)

    x = CA{Int, 1, UInt8}(1:254)
    res = @test_throws LevelsException{Int, UInt8} x[1:2] = 1000:1001
    @test res.value.levels == [1001]
    @test sprint(showerror, res.value) == "cannot store level(s) 1001 since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."
    @test x == A(vcat(1000, 2:254))
    @test levels(x) == vcat(1:254, 1000)

    x = CA{Int, 1, UInt8}([1, 3, 256])
    res = @test_throws LevelsException{Int, UInt8} levels!(x, collect(1:256))
    @test res.value.levels == [255]
    @test sprint(showerror, res.value) == "cannot store level(s) 255 since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."

    x = CA(30:2:131115)
    res = @test_throws LevelsException{Int, UInt16} CategoricalVector{Int, UInt16}(x)
    @test res.value.levels == collect(131100:2:131114)
    @test sprint(showerror, res.value) == "cannot store level(s) 131100, 131102, 131104, 131106, 131108, 131110, 131112 and 131114 since reference type UInt16 can only hold 65535 levels. Use the decompress function to make room for more levels."

    x = CA{String, 1, UInt8}(string.(Char.(65:318)))
    res = @test_throws LevelsException{String, UInt8} levels!(x, vcat(levels(x), "az", "bz", "cz"))
    @test res.value.levels == ["bz", "cz"]
    @test sprint(showerror, res.value) == "cannot store level(s) \"bz\" and \"cz\" since reference type UInt8 can only hold 255 levels. Use the decompress function to make room for more levels."
    @test x == A(string.(Char.(65:318)))
    lev = copy(levels(x))
    levels!(x, vcat(lev, "az"))
    @test levels(x) == vcat(lev, "az")

    x = compress(CA([1, 3, 736251]))
    ux = decompress(x)
    @test x == ux
    @test typeof(x) == CA{Int, 1, UInt8}
    @test typeof(ux) == CA{Int, 1, CategoricalArrays.DefaultRefType}

    # Test reshape()
    x = CA(["Old", "Young", "Middle", "Young"])
    levels!(x, ["Young", "Middle", "Old"])
    ordered!(x, true)

    y = reshape(x, 1, 4)
    @test isa(y, CA{String, 2, CategoricalArrays.DefaultRefType})
    @test y == A(["Old" "Young" "Middle" "Young"])
    @test levels(x) == levels(y)
    @test isordered(x)

    y = reshape(x, 2, 2)
    @test isa(y, CA{String, 2, CategoricalArrays.DefaultRefType})
    @test y == A(["Old"  "Middle"; "Young" "Young"])
    @test levels(x) == levels(y)
    @test isordered(x)

    # Test with null values
    if CA === NullableCategoricalArray
        x[3] = Nullable()
        y = reshape(x, 1, 4)
        @test isa(y, CA{String, 2, CategoricalArrays.DefaultRefType})
        @test y == A(["Old" "Young" Nullable() "Young"])
        @test levels(x) == levels(y)
        @test isordered(x)
    end

    # Test sort() on both unordered and ordered arrays
    x = CA(["Old", "Young", "Middle", "Young"])
    levels!(x, ["Young", "Middle", "Old"])
    @test sort(x) == A(["Young", "Young", "Middle", "Old"])
    ordered!(x, true)
    @test sort(x) == A(["Young", "Young", "Middle", "Old"])
    @test sort!(x) === x
    @test x == A(["Young", "Young", "Middle", "Old"])

    # Test constructors and convert() between categorical arrays:
    # check that they preserve levels and ordering by default even across types,
    # do not modify original array when changing ordering, and make copies
    # only when necessary (though for construcors the pool is always copied because of ordered)
    for ordered_orig in (true, false),
        ordered in (true, false),
        R in (DefaultRefType, UInt8, UInt, Int8, Int),
        (A2, CA2, CVM2, N) in ((Array, CategoricalArray, CategoricalVector, 1),
                               (Array, CategoricalArray, CategoricalMatrix, 2),
                               (NullableArray, NullableCategoricalArray, NullableCategoricalVector, 1),
                               (NullableArray, NullableCategoricalArray, NullableCategoricalMatrix, 2))
        if CVM2 <: AbstractVector
            x = CA{String, 1, R}(["A", "B"], ordered=ordered_orig)
        else
            x = CA{String, 2, R}(["A" "B"], ordered=ordered_orig)
        end
        # Do not use lexicographic order in order to catch bugs about order preservation
        levels!(x, ["B", "A"])

        for y in (CA2(x),
                  CA2{String}(x),
                  CA2{String, N}(x),
                  CA2{String, N, R}(x),
                  CA2{String, N, DefaultRefType}(x),
                  CA2{String, N, UInt8}(x),
                  CA2(x),
                  CA2{String}(x),
                  CVM2{String, R}(x),
                  CVM2{String, DefaultRefType}(x),
                  CVM2{String, UInt8}(x))
            @test isa(y, CVM2{String})
            @test isordered(y) === isordered(x)
            @test isordered(x) === ordered_orig
            @test y.refs == x.refs
            @test index(y.pool) == index(x.pool)
            @test levels(y) == levels(x)
            @test (y.refs === x.refs) == (eltype(x.refs) === eltype(y.refs))
            @test y.pool !== x.pool
        end
        for y in (categorical(x),
                  categorical(x, false),
                  categorical(x, true))
            @test isa(y, CA{String})
            @test isordered(y) === isordered(x)
            @test isordered(x) === ordered_orig
            @test y.refs == x.refs
            @test index(y.pool) == index(x.pool)
            @test levels(y) == levels(x)
            @test (y.refs === x.refs) == (eltype(x.refs) === eltype(y.refs))
            @test y.pool !== x.pool
        end
        for y in (CA2(x, ordered=ordered),
                  CA2{String}(x, ordered=ordered),
                  CA2{String, N}(x, ordered=ordered),
                  CA2{String, N, R}(x, ordered=ordered),
                  CA2{String, N, DefaultRefType}(x, ordered=ordered),
                  CA2{String, N, UInt8}(x, ordered=ordered),
                  CVM2(x, ordered=ordered),
                  CVM2{String}(x, ordered=ordered),
                  CVM2{String, R}(x, ordered=ordered),
                  CVM2{String, DefaultRefType}(x, ordered=ordered),
                  CVM2{String, UInt8}(x, ordered=ordered))
            @test isa(y, CVM2{String})
            @test isordered(y) === ordered
            @test isordered(x) === ordered_orig
            @test y.refs == x.refs
            @test index(y.pool) == index(x.pool)
            @test levels(y) == levels(x)
            @test (y.refs === x.refs) == (eltype(x.refs) === eltype(y.refs))
            @test y.pool !== x.pool
        end
        for y in (categorical(x, ordered=ordered),
                  categorical(x, false, ordered=ordered),
                  categorical(x, true, ordered=ordered))
            @test isa(y, CA{String})
            @test isordered(y) === ordered
            @test isordered(x) === ordered_orig
            @test y.refs == x.refs
            @test index(y.pool) == index(x.pool)
            @test levels(y) == levels(x)
            @test (y.refs === x.refs) == (eltype(x.refs) === eltype(y.refs))
            @test y.pool !== x.pool
        end
        for y in (convert(CA2, x),
                  convert(CA2{String}, x),
                  convert(CA2{String, N}, x),
                  convert(CA2{String, N, R}, x),
                  convert(CA2{String, N, DefaultRefType}, x),
                  convert(CA2{String, N, UInt8}, x))
            @test isa(y, CVM2{String})
            @test isordered(y) === isordered(x)
            @test isordered(x) === ordered_orig
            @test y.refs == x.refs
            @test index(y.pool) == index(x.pool)
            @test levels(y) == levels(x)
            @test (y.refs === x.refs) == (eltype(x.refs) === eltype(y.refs))
            @test (y.pool === x.pool) == (eltype(x.refs) === eltype(y.refs))
        end
    end
end

# Check that converting from NullableCategoricalArray to CategoricalArray fails with nulls
x = NullableCategoricalArray(1)
@test_throws NullException CategoricalArray(x)
@test_throws NullException convert(CategoricalArray, x)

# Test in()
ca1 = CategoricalArray([1, 2, 3])
ca2 = CategoricalArray([4, 3, 2])

@test (ca1[1] in ca1) === true
@test (ca2[2] in ca1) === true
@test (ca2[1] in ca1) === false

@test (1 in ca1) === true
@test (5 in ca1) === false

nca1 = NullableCategoricalArray([1, 2, 3])
nca2 = NullableCategoricalArray([4, 3, 2])

@test (ca1[1] in nca1) === false
@test (1 in nca1) === false

@test (nca1[1] in nca1) === true
@test (nca2[2] in nca1) === true
@test (nca2[1] in nca1) === false

@test (Nullable(1) in nca1) === true
@test (Nullable(5) in nca1) === false

end


text:
#
# Storage Driver implementation for trees which are rooted in the file system
# (in git terminology, there exists a "working copy")
#
abstract type AbstractFileSystemRoot end

# These functions sys_abspath and sys_joinpath generate/joins OS-specific
# _local filesystem paths_ out of logical paths. They should be defined only
# for trees which are rooted in the actual filesystem.
function sys_abspath(root::AbstractFileSystemRoot, path::RelPath)
    rootpath = sys_abspath(root)
    return isempty(path.components) ? rootpath : joinpath(rootpath, sys_joinpath(path))
end

sys_joinpath(path::RelPath) = isempty(path.components) ? "" : joinpath(path.components...)
sys_abspath(path::AbsPath) = sys_abspath(path.root, path.path)
sys_abspath(tree::BlobTree) = sys_abspath(tree.root, tree.path)
sys_abspath(file::Blob) = sys_abspath(file.root, file.path)

#--------------------------------------------------
# Storage data interface for trees
#
# TODO: Formalize this interface!

## 1. Query

# TODO: would it be better to express the following dispatch in terms of
# AbsPath{<:AbstractFileSystemRoot} rather than usin double dispatch?

Base.isdir(root::AbstractFileSystemRoot, path::RelPath) = isdir(sys_abspath(root, path))
Base.isfile(root::AbstractFileSystemRoot, path::RelPath) = isfile(sys_abspath(root, path))
Base.ispath(root::AbstractFileSystemRoot, path::RelPath) = ispath(sys_abspath(root, path))

Base.summary(io::IO, root::AbstractFileSystemRoot) = print(io, sys_abspath(root))

Base.readdir(root::AbstractFileSystemRoot, path::RelPath) = readdir(sys_abspath(root, path))

## 2. Mutation
#
# TODO: Likely requires rework!

function Base.mkdir(root::AbstractFileSystemRoot, path::RelPath; kws...)
    if !iswriteable(root)
        error("Cannot make directory in read-only tree root at $(sys_abspath(p.root))")
    end
    mkdir(sys_abspath(root, path), args...)
    return BlobTree(root, path)
end

function Base.rm(root::AbstractFileSystemRoot, path::RelPath; kws...)
    rm(sys_abspath(root,path); kws...)
end

#--------------------------------------------------
# Storage data interface for Blob

# TODO: Make this the generic implementation for AbstractDataStorage
function Base.open(f::Function, as_type::Type{IO},
                   root::AbstractFileSystemRoot, path; kws...)
    @context f(@! open(as_type, root, path; kws...))
end

@! function Base.open(::Type{IO}, root::AbstractFileSystemRoot, path;
                      write=false, read=!write, kws...)
    if !iswriteable(root) && write
        error("Error writing file at read-only path $path")
    end
    @! open(sys_abspath(root, path); read=read, write=write, kws...)
end

Base.read(root::AbstractFileSystemRoot, path::RelPath, ::Type{T}) where {T} =
    read(sys_abspath(root, path), T)
Base.read(root::AbstractFileSystemRoot, path::RelPath) =
    read(sys_abspath(root, path))

#--------------------------------------------------
"""

## Metadata spec

For Blob:
```
    [datasets.storage]
    driver="FileSystem"
    type="Blob"
    path=\$(path_to_file)
```

For BlobTree:
```
    [datasets.storage]
    driver="FileSystem"
    type="BlobTree"
    path=\$(path_to_directory)
```
"""
struct FileSystemRoot <: AbstractFileSystemRoot
    path::String
    read::Bool
    write::Bool
end

function FileSystemRoot(path::AbstractString; write=false, read=true)
    path = abspath(path)
    FileSystemRoot(path, read, write)
end

iswriteable(root::FileSystemRoot) = root.write

sys_abspath(root::FileSystemRoot) = root.path

function Base.abspath(relpath::RelPath)
    Base.depwarn("""
        `abspath(::RelPath)` defaults to using `pwd()` as the root of the path
        but this leads to fragile code so will be removed in the future""",
        :abspath)
    AbsPath(FileSystemRoot(pwd(); write=true, read=true), relpath)
end

#-------------------------------------------------------------------------------
# Infrastructure for a somewhat more functional interface for creating file
# trees than the fully mutable version we usually use.

mutable struct TempFilesystemRoot <: AbstractFileSystemRoot
    path::Union{Nothing,String}
    function TempFilesystemRoot(path)
        root = new(path)
        finalizer(root) do r
            if !isnothing(r.path)
                rm(r.path, recursive=true, force=true)
            end
        end
        return root
    end
end

function Base.readdir(root::TempFilesystemRoot, path::RelPath)
    return isnothing(root.path) ? [] : readdir(sys_abspath(root, path))
end

iswriteable(root::TempFilesystemRoot) = true
sys_abspath(root::TempFilesystemRoot) = root.path

"""
    newdir()

Create a new temporary `BlobTree` which can have files assigned into it and may
be assigned to a permanent location in a persistent `BlobTree`. If not assigned
to a permanent location, the temporary tree is cleaned up during garbage
collection.
"""
function newdir(ctx::AbstractFileSystemRoot=FileSystemRoot(tempdir(), write=true))
    # cleanup=false: we manage our own cleanup via the finalizer
    path = mktempdir(sys_abspath(ctx), cleanup=false)
    return BlobTree(TempFilesystemRoot(path))
end
newdir(ctx::BlobTree) = newdir(ctx.root)

function newfile(ctx::AbstractFileSystemRoot=FileSystemRoot(tempdir(), write=true))
    path, io = mktemp(sys_abspath(ctx), cleanup=false)
    close(io)
    return Blob(TempFilesystemRoot(path))
end
newfile(ctx::BlobTree) = newfile(ctx.root)

"""
    newfile(func)
    newfile(func, ctx)

Create a new temporary `Blob` object which may be later assigned to a permanent
location in a `BlobTree`. If not assigned to a permanent location, the
temporary file is cleaned up during garbage collection.

# Example

```
tree[path"some/demo/path.txt"] = newfile() do io
    println(io, "Hi there!")
end
```
"""
function newfile(f::Function, ctx=FileSystemRoot(tempdir(), write=true))
    path, io = mktemp(sys_abspath(ctx), cleanup=false)
    try
        f(io)
    catch
        rm(path)
        rethrow()
    finally
        close(io)
    end
    return Blob(TempFilesystemRoot(path))
end

# Move srcpath to destpath, making all attempts to preserve the original
# content of `destpath` if anything goes wrong. We assume that `srcpath` is
# temporary content which doesn't need to be protected.
function mv_force_with_dest_rollback(srcpath, destpath, tempdir_parent)
    holding_area = nothing
    held_path = nothing
    if ispath(destpath)
        # If the destination path exists, improve the atomic nature of the
        # update by first moving existing data to a temporary directory.
        holding_area = mktempdir(tempdir_parent, prefix="jl_to_remove_", cleanup=false)
        name = basename(destpath)
        held_path = joinpath(holding_area,name)
        mv(destpath, held_path)
    end
    try
        mv(srcpath, destpath)
    catch
        try
            if !isnothing(holding_area)
                # Attempt to put things back as they were!
                mv(held_path, destpath)
            end
        catch
            # At this point we've tried our best to preserve the user's data
            # but something has gone wrong, likely at the OS level. The user
            # will have to clean up manually if possible.
            error("""
                  Something when wrong while moving data to path $destpath.

                  We tried restoring the original data to $destpath, but were
                  met with another error. The original data is preserved in
                  $held_path

                  See the catch stack for the root cause.
                  """)
        end
        rethrow()
    end
    if !isnothing(holding_area)
        # If we get to here, it's safe to remove the holding area
        rm(holding_area, recursive=true)
    end
end

function Base.setindex!(tree::BlobTree{<:AbstractFileSystemRoot},
                        tmpdata::Union{Blob{TempFilesystemRoot},BlobTree{TempFilesystemRoot}},
                        name::AbstractString)
    if !iswriteable(tree.root)
        error("Attempt to move to a read-only tree $tree")
    end
    if isnothing(tmpdata.root.path)
        type = isdir(tmpdata) ? "directory" : "file"
        error("Attempted to root a temporary $type which has already been moved to $(tree.path)/$name ")
    end
    if !isempty(tree.path)
        # Eh, the number of ways the user can misuse this isn't really funny :-/
        error("Temporary trees must be moved in full. The tree had non-empty path $(tree.path)")
    end
    destpath = sys_abspath(joinpath(tree, name))
    srcpath = sys_abspath(tmpdata)
    tempdir_parent = sys_abspath(tree)
    mv_force_with_dest_rollback(srcpath, destpath, tempdir_parent)
    # Transfer ownership of the data to `tree`. This is ugly to be sure, as it
    # leaves `tmpdata` empty! However, we'll have to live with this wart unless
    # we want to be duplicating large amounts of data on disk.
    tmpdata.root.path = nothing
    return tree
end

# It's interesting to read about the linux VFS interface in regards to how the
# OS actually represents these things. For example
# https://stackoverflow.com/questions/36144807/why-does-linux-use-getdents-on-directories-instead-of-read

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

# Filesystem storage driver
function connect_filesystem(f, config, _)
    path = config["path"]
    type = config["type"]
    if type == "Blob"
        isfile(path) || throw(ArgumentError("$(repr(path)) should be a file"))
        storage = Blob(FileSystemRoot(path))
    elseif type == "BlobTree"
        isdir(path)  || throw(ArgumentError("$(repr(path)) should be a directory"))
        storage = BlobTree(FileSystemRoot(path))
    else
        throw(ArgumentError("DataSet type $type not supported on the filesystem"))
    end
    f(storage)
end
add_storage_driver("FileSystem"=>connect_filesystem)


text:
using Pkg

pkg"""
add Documenter Literate Plots Yao
"""

using Documenter, Literate, QuAlgorithmZoo

const Examples = ["Grover", "VQE", "Shor"]

const PATH = (
    tutorial = joinpath(@__DIR__, "src/tutorial"),
    examples = joinpath(@__DIR__, "..", "examples")
)

function process_literate_scripts()
    TUTORIALS = []
    for token in Examples
        file = "$token.jl"
        filepath = joinpath(PATH.examples, token, file)
        Literate.markdown(filepath, PATH.tutorial)

        filename, _ = splitext(file)
        mdfile = join([filename, ".md"])
        # TODO: use PATH.tutorial rather then manual path
        push!(TUTORIALS, relpath(joinpath("tutorial", mdfile)))
    end
    TUTORIALS
end

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

function generate(islocal::Bool="local" in ARGS)
    makedocs(
        modules = [QuAlgorithmZoo],
        clean = false,
        format = :html,
        sitename = "Quantum Algorithm Zoo",
        linkcheck = !("skiplinks" in ARGS),
        analytics = "UA-89508993-1",
        pages = [
            "Home" => "index.md",
            "Algorithms" => process_literate_scripts(),
            "Manual" => Any[
                "man/zoo.md",
            ],
        ],
        html_prettyurls = !islocal,
        html_canonical = "https://quantumbfs.github.io/QuAlgorithmZoo.jl/latest/",
    )

    deploydocs(
        repo = "github.com/QuantumBFS/QuAlgorithmZoo.jl.git",
        target = "build",
        julia = "1.0",
        deps = nothing,
        make = nothing,
    )
end

generate(false)


text:
import Base:*

"""

# stan

Execute a Stan model.

### Method
```julia
rc, sim, cnames = stan(
  model::Stanmodel,
  data=Nothing,
  ProjDir=pwd();
  init=Nothing,
  summary=true,
  diagnostics=false,
  CmdStanDir=CMDSTAN_HOME,
  file_run_log=true
)
```
### Required arguments
```julia
* `model::Stanmodel`              : See ?Method
```

### Optional positional arguments

```julia
* `data=Nothing`                  : Observed input data dictionary
* `ProjDir=pwd()`                 : Project working directory
```

### Keyword arguments
```julia
* `init=Nothing`                  : Initial parameter value dictionary
* `summary=true`                  : Use CmdStan's stansummary to display results
* `diagnostics=false`             : Generate diagnostics file
* `CmdStanDir=CMDSTAN_HOME`       : Location of cmdstan directory
* `file_run_log=true`             : Create run log file (if false, write log to stdout)
* `file_make_log=true`            : Create make log file (if false, write log to stdout)
```

If the type the `data` or `init` arguments are an AbstartcString this is interpreted as
a path name to an existing .R file. This file is copied to the coresponding .R data
and/or init files for for each chain.

### Return values
```julia
* `rc::Int`                       : Return code from stan(), rc == 0 if all is well
* `sim`                           : Chain results
* `cnames`                        : Vector of variable names
```

### Examples

```julia
# no data, use default ProjDir (pwd)
stan(mymodel)
# default ProjDir (pwd)
stan(mymodel, mydata)
# specify ProjDir
stan(mymodel, mydata, "~/myproject/")
# keyword arguments
stan(mymodel, mydata, "~/myproject/", diagnostics=true, summary=false)
# use default ProjDir (pwd), with keyword arguments
stan(mymodel, mydata, diagnostics=true, summary=false)
```

### Related help
```julia
?Stanmodel                         : Create a StanModel
```
"""
function stan(
  model::Stanmodel,
  data=Nothing,
  ProjDir=pwd();
  init=Nothing,
  summary=true,
  diagnostics=false,
  CmdStanDir=CMDSTAN_HOME,
  file_run_log=true,
  file_make_log=true)

  old = pwd()
  local rc = 0
  local res = Dict[]

  if length(model.model) == 0
    println("\nNo proper model specified in \"$(model.name).stan\".")
    println("This file is typically created from a String passed to Stanmodel().\n")
    return((-1, res))
  end

  @assert isdir(ProjDir) "Incorrect ProjDir specified: $(ProjDir)"
  @assert isdir(model.tmpdir) "$(model.tmpdir) not created"

  cd(model.tmpdir)
  isfile("$(model.name)_build.log") && rm("$(model.name)_build.log")
  isfile("$(model.name)_make.log") && rm("$(model.name)_make.log")
  isfile("$(model.name)_run.log") && rm("$(model.name)_run.log")

  cd(CmdStanDir)
  local tmpmodelname::String
  tmpmodelname = joinpath(model.tmpdir, model.name)
  if @static Sys.iswindows() ? true : false
    tmpmodelname = replace(tmpmodelname*".exe", "\\" => "/")
  end
  try
    if file_make_log
      run(pipeline(`make $(tmpmodelname)`,
        stdout="$(tmpmodelname)_make.log",
        stderr="$(tmpmodelname)_build.log"))
    else
      run(pipeline(`make $(tmpmodelname)`,
        stderr="$(tmpmodelname)_build.log"))
    end
  catch
    println("\nAn error occurred while compiling the Stan program.\n")
    print("Please check your Stan program in variable '$(model.name)' ")
    print("and the contents of $(tmpmodelname)_build.log.\n")
    println("Note that Stan does not handle blanks in path names.")
    error("Return code = -3");
  end

  cd(model.tmpdir)

  if data != Nothing && (typeof(data) <: AbstractString || check_dct_type(data))
    if typeof(data) <: AbstractString
      for i in 1:model.nchains
        cp(data, "$(model.name)_$(i).data.R", force=true)
      end
    else
      if typeof(data) <: Array && length(data) == model.nchains
        for i in 1:model.nchains
          if length(keys(data[i])) > 0
            update_R_file("$(model.name)_$(i).data.R", data[i])
          end
        end
      else
        if typeof(data) <: Array
          for i in 1:model.nchains
            if length(keys(data[1])) > 0
              update_R_file("$(model.name)_$(i).data.R", data[1])
            end
          end
        else
          for i in 1:model.nchains
            if length(keys(data)) > 0
              update_R_file("$(model.name)_$(i).data.R", data)
            end
          end
        end
      end
    end
  end

  if init != Nothing && (typeof(init) <: AbstractString || check_dct_type(init))
    if typeof(init) <: AbstractString
      for i in 1:model.nchains
        cp(init, "$(model.name)_$(i).init.R", force=true)
      end
    else
      if typeof(init) <: Array && length(init) == model.nchains
        for i in 1:model.nchains
          if length(keys(init[i])) > 0
            update_R_file("$(model.name)_$(i).init.R", init[i])
          end
        end
      else
        if typeof(init) <: Array
          for i in 1:model.nchains
            if length(keys(init[1])) > 0
              update_R_file("$(model.name)_$(i).init.R", init[1])
            end
          end
        else
          for i in 1:model.nchains
            if length(keys(init)) > 0
              update_R_file("$(model.name)_$(i).init.R", init)
            end
          end
        end
      end
    end
  end

  for i in 1:model.nchains
    model.id = i
    model.data_file ="$(model.name)_$(i).data.R"
    if !(model.init == Nothing)
        model.init_file = "$(model.name)_$(i).init.R"
    end
    if isa(model.method, Sample)
      model.output.file = model.name*"_samples_$(i).csv"
      isfile("$(model.name)_samples_$(i).csv") && rm("$(model.name)_samples_$(i).csv")
      if diagnostics
        model.output.diagnostic_file = model.name*"_diagnostics_$(i).csv"
        isfile("$(model.name)_diagnostics_$(i).csv") && rm("$(model.name)_diagnostics_$(i).csv")
      end
    elseif isa(model.method, Optimize)
      isfile("$(model.name)_optimize_$(i).csv") && rm("$(model.name)_optimize_$(i).csv")
      model.output.file = model.name*"_optimize_$(i).csv"
    elseif isa(model.method, Variational)
      isfile("$(model.name)_variational_$(i).csv") && rm("$(model.name)_variational_$(i).csv")
      model.output.file = model.name*"_variational_$(i).csv"
    elseif isa(model.method, Diagnose)
      isfile("$(model.name)_diagnose_$(i).csv") && rm("$(model.name)_diagnose_$(i).csv")
      model.output.file = model.name*"_diagnose_$(i).csv"
    end
    model.command[i] = cmdline(model)
  end

  try
    if file_run_log
      run(pipeline(par(model.command), stdout="$(model.name)_run.log"))
    else
      run(par(model.command))
    end
  catch e
    println("\nAn error occurred while running the previously compiled Stan program.\n")
    print("Please check the contents of file $(tmpmodelname)_run.log and the")
    println("'command' field in the Stanmodel, e.g. stanmodel.command.\n")
    cd(old)
    error("Return code = -5")
  end

  local samplefiles = String[]
  local ftype
  local cnames = String[]

  if typeof(model.method) in [Sample, Variational]
    if isa(model.method, Sample)
      ftype = diagnostics ? "diagnostics" : "samples"
    else
      ftype = lowercase(string(typeof(model.method)))
    end

    for i in 1:model.nchains
      push!(samplefiles, "$(model.name)_$(ftype)_$(i).csv")
    end

    if summary
      stan_summary(model, par(samplefiles), CmdStanDir=CmdStanDir)
    end

    (res, cnames) = read_samples(model, diagnostics)

  elseif isa(model.method, Optimize)
    res, cnames = read_optimize(model)

  elseif isa(model.method, Diagnose)
    res, cnames = read_diagnose(model)

  else
    println("\nAn unknown method is specified in the call to stan().")
    cd(old)
    error("Return code = -10")
  end

  if model.output_format != :array
    start_sample = 1
    if isa(model.method, Sample) && !model.method.save_warmup
      start_sample = model.method.num_warmup+1
    end
    res = convert_a3d(res, cnames, Val(model.output_format);
      start=start_sample)
  end

  cd(old)

  (rc, res, cnames)

end


text:
# If you want a scatter plot of the spectra, define `plotresults = true`
# prior to running this script.

using CompScienceMeshes, BEAST
using LinearAlgebra

# Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
Γ = meshrectangle(1.0, 1.0, 0.2, 3)
println("Mesh with $(numvertices(Γ)) vertices and $(numcells(Γ)) cells.")
X = raviartthomas(Γ)
Y = BEAST.buffachristiansen3(Γ)

κ = ω = 1.0; γ = κ*im
T = MWSingleLayer3D(γ)
N = NCross()

Txx = assemble(T,X,X); println("primal discretisation assembled.")
Tyy = assemble(T,Y,Y); println("dual discretisation assembled.")
Nxy = assemble(N,X,Y); println("duality form assembled.")

iNxy = inv(Nxy); println("duality form inverted.")
A = iNxy' * Tyy * iNxy * Txx

@show cond(Txx)
@show cond(Tyy)
@show cond(A)

wx, wy, wp = eigvals(Txx), eigvals(Tyy), eigvals(A); println("eigvals found.")

(@isdefined plotresults) || (plotresults = false)
if plotresults
    @eval begin
        using Plots
        plot()
        scatter!(wx,c=:blue,label="primal")
        scatter!(wy,c=:red,label="dual")
        scatter!(wp,c=:green,label="CMP")
        plot!(xlim=(-0.5,0.5),ylim=(-10.0,10.0))
    end
end


text:
# File Names
DataDir = joinpath(dirname(@__FILE__),"..","data")
CompressFile = joinpath(DataDir,"bindata.tar.gz")
CifarOrigDir = joinpath(DataDir,"cifar-10-batches-bin")
CifarFinalDir = joinpath(DataDir,"bin")

# Download Compressed Files
download("http://cbcl.mit.edu/projects/cbcl/software-datasets/faces.tar.gz",CompressFile)

# # Decompress Dataset
cd(DataDir)
run(`tar -zxvf $CompressFile`)
run(`tar -zxvf face.test.tar.gz`)
run(`tar -zxvf face.train.tar.gz`)


text:
abstract type RandBase{T} <: UnaryNodeOp{T} end

mutable struct RandState{RNG<:Random.AbstractRNG} <: NodeEvaluationState
    rng::RNG
end

always_ticks(::RandBase) = true
time_agnostic(::RandBase) = true
value_agnostic(::RandBase) = true
# Note that we should never mutate the random state on the node op itself.
create_operator_evaluation_state(::Tuple{Node}, op::RandBase) = RandState(copy(op.rng))

struct Rand{T<:Real} <: RandBase{T}
    rng::Random.AbstractRNG
end

Base.hash(x::Rand, h::UInt64) = hash(x.rng, hash(:Rand, h))
Base.:(==)(x::Rand{T}, y::Rand{T}) where {T} = x.rng == y.rng

operator!(::Rand{T}, state::RandState) where {T} = rand(state.rng, T)

struct RandArray{T,N} <: RandBase{Array{T,N}}
    rng::Random.AbstractRNG
    dims::Dims{N}
end

Base.hash(x::RandArray, h::UInt64) = hash(x.rng, hash(x.dims, hash(:RandArray, h)))
Base.:(==)(x::RandArray{T}, y::RandArray{T}) where {T} = x.rng == y.rng && x.dims == y.dims

operator!(op::RandArray{T}, state::RandState) where {T} = rand(state.rng, T, op.dims)

"""
    rand([rng=MersenneTwister(),] alignment::Node[, S, dims...])

Generate random numbers aligned to `alignment`, with the given `rng` if provided.

Semantics are otherwise very similar to the usual `Base.rand`:
* If specified, `S` will be the element type, and will default to `Float64` otherwise.
* If specified, `dims` should be a tuple or vararg of integers representing the dimensions
    of an array.

**NB** The values of `alignment` will be ignored.

!!! tip
    If provided, `rng` will be copied before it is used. This is to ensure reproducability
    when evaluating a node multiple times.
"""
Base.rand(alignment::Node) = rand(MersenneTwister(), alignment)
Base.rand(rng::Random.AbstractRNG, alignment::Node) = rand(rng, alignment, Float64)
function Base.rand(rng::Random.AbstractRNG, alignment::Node, ::Type{X}) where {X}
    return obtain_node((alignment,), Rand{X}(copy(rng)))
end

function Base.rand(alignment::Node, dims::Integer...)
    return rand(MersenneTwister(), alignment, Dims(dims))
end
Base.rand(alignment::Node, dims::Dims) = rand(MersenneTwister(), alignment, dims)
function Base.rand(rng::Random.AbstractRNG, alignment::Node, dims::Dims)
    return rand(rng, alignment, Float64, dims)
end
function Base.rand(
    rng::Random.AbstractRNG, alignment::Node, ::Type{X}, dims::Integer...
) where {X}
    return rand(rng, alignment, X, Dims(dims))
end
function Base.rand(
    rng::Random.AbstractRNG, alignment::Node, ::Type{X}, dims::Dims
) where {X}
    return obtain_node((alignment,), RandArray{X,length(dims)}(copy(rng), dims))
end


text:
@testset "780.reaching-points.jl" begin
    @test reaching_points(1,4,999999997,4) == true
    @test reaching_points(9,5,12,8) == false
    @test reaching_points(1,1,3,5) == true
    @test reaching_points(3,3,12,9) == true
    @test reaching_points(9, 10, 9, 19) == true
end


text:
benchmarksdir() = joinpath(@__DIR__, "..", "benchmark")

function run_benchmarks(
    outdir = default_benchmarkstore();
    nthreadsrange = 2 .^ (1:floor(Int, log2(cld(Sys.CPU_THREADS, 2)))),
)
    julia = `$(Base.julia_cmd()) --startup-file=no`
    for name in ["task_overheads", "notify_overheads"]
        for nthreads in nthreadsrange
            @info "Executing benchmark `$name` with $nthreads threads"
            outfile = joinpath(outdir, string(nthreads), name * ".json")
            script = joinpath(benchmarksdir(), name * ".jl")
            @time run(`$julia --threads=$nthreads $script $outfile`)
        end
    end
end


text:
using Stheno: FiniteMean, ConstantMean, AM, AV, pairwise

@testset "finite" begin

    # Tests for FiniteMean.
    let
        rng, N, D = MersenneTwister(123456), 5, 2
        μ, X, x = ConstantMean(1.5), ColsAreObs(randn(rng, D, N)), randn(rng, N)
        μX, μx = FiniteMean(μ, X), FiniteMean(μ, x)

        # Recusion depth 1.
        for (μ′, X′) in zip([μX, μx], [X, x])
            @test length(μ′) == N
            @test eachindex(μ′) == 1:N
            @test AbstractVector(μ′) == map(μ′, eachindex(μ′))
            @test AbstractVector(μ′) == map(μ, X′)
            mean_function_tests(μ′, eachindex(μ′))
            mean_function_tests(μ′, 1:N-1)
            mean_function_tests(μ′, BlockData([1:N-1, eachindex(μ′)]))

            # Recursion depth 2.
            r = 1:N-1
            μ2 = FiniteMean(μ′, r)
            @test length(μ2) == N - 1
            @test eachindex(μ2) == 1:N-1
            @test AbstractVector(μ′)[r] == AbstractVector(μ2)
            @test AbstractVector(μ′)[r] == map(μ2, eachindex(μ2))
            mean_function_tests(μ2, eachindex(μ2))
            mean_function_tests(μ2, 1:N-2)
            mean_function_tests(μ2, 1:N-2)
        end

        @test map(μX, 1:N-1) == map(μ, X[1:N-1])
        @test map(μX, 2:N) == map(μ, X[2:N])
        @test map(μX, 2:N-2) == map(μ, X[2:N-2])

        @test map(μx, 1:N-1) == map(μ, x[1:N-1])
        @test map(μx, 2:N) == map(μ, x[2:N])
        @test map(μx, 2:N-2) == map(μ, x[2:N-2])

        show(IOBuffer(), μX)
    end

    # Tests for FiniteKernel.
    let
        rng, N, D = MersenneTwister(123456), 5, 2
        k, X, x = EQ(), ColsAreObs(randn(rng, D, N)), randn(rng, N)
        kX, kx = FiniteKernel(k, X), FiniteKernel(k, x)

        for (k′, X′) in zip((kX, kx), (X, x))

            # Check for correctness relative to base kernel.
            @test size(k′) == (N, N)
            @test size(k′, 1) == N
            @test size(k′, 2) == N
            @test isstationary(k′) == false
            @test eachindex(k′) == 1:N
            @test AM(k′) == pairwise(k′, eachindex(k′))
            @test Matrix(AM(k′)) == pairwise(k′, eachindex(k′), eachindex(k′))
            @test AM(k′) ≈ pairwise(k, X′)
            @test AM(k′) ≈ pairwise(k, X′)

            @test pairwise(k′, 1:N-1, 2:N) ≈ pairwise(k, X′[1:N-1], X′[2:N])
            @test pairwise(k′, eachindex(k′), 1:N-1) ≈ pairwise(k, X′, X′[1:N-1])
            @test pairwise(k′, 2:N, eachindex(k′)) ≈ pairwise(k, X′[2:N], X′)

            kernel_tests(k′, 1:N-1, 2:N, eachindex(k′))
            d0 = BlockData([1:N-1, 2:N])
            d1 = BlockData([2:N, 1:N-1])
            d2 = BlockData([eachindex(k′)])
            kernel_tests(k′, d0, d1, d2)

            # Recursion depth 2.
            r = 1:N-1
            k2 = FiniteKernel(k′, r)
            @test size(k2) == (N-1, N-1)
            @test eachindex(k2) == r
            @test AM(k2) == pairwise(k′, r)

            # Known issue: offsetting doesn't really work.
            binary_map_tests(k2, r, r)
            pairwise_tests(k′, r, r)
        end

        show(IOBuffer(), kX)
    end

    # Tests for LhsFiniteCrossKernel.
    let
        rng, N, N′, D = MersenneTwister(123456), 4, 5, 2
        k, X, X′ = EQ(), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        x, x′ = randn(rng, N), randn(rng, N′)
        k′, kx = LhsFiniteCrossKernel(k, X), LhsFiniteCrossKernel(k, x)

        @test size(k′) == (N, Inf)
        @test size(k′, 1) == N
        @test size(k′, 2) == Inf
        @test isstationary(k′) == false
        @test eachindex(k′, 1) == 1:N
        @test pairwise(k′, eachindex(k′, 1), X′) ≈ pairwise(k, X, X′)
        @test pairwise(k′, 1:N-1, X′) ≈ pairwise(k, X[1:N-1], X′)

        cross_kernel_tests(k′, eachindex(k′, 1), X, X′)
        cross_kernel_tests(kx, eachindex(kx, 1), x, x′)
        cross_kernel_tests(k′, eachindex(k′, 1), BlockData([X]), BlockData([X′]))
        cross_kernel_tests(kx, eachindex(kx, 1), BlockData([x]), BlockData([x′]))

        show(IOBuffer(), k′)
    end

    # Tests for RhsFiniteCrossKernel.
    let
        rng, N, N′, D = MersenneTwister(123456), 4, 5, 2
        X, Y, X′ = ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        x, y, x′ = randn(rng, N), randn(rng, N), randn(rng, N′)
        k = EQ()
        k′, kx = RhsFiniteCrossKernel(k, X′), RhsFiniteCrossKernel(k, x′)

        @test size(k′) == (Inf, N′)
        @test size(k′, 1) == Inf
        @test size(k′, 2) == N′
        @test isstationary(k′) == false
        @test eachindex(k′, 2) == 1:N′
        @test pairwise(k′, X, eachindex(k′, 2)) ≈ pairwise(k, X, X′)
        @test pairwise(k′, X, 1:N′-1) ≈ pairwise(k, X, X′[1:N′-1])

        d1, d2 = 1:N, eachindex(k′, 2)
        cross_kernel_tests(k′, X, d1, d2)
        cross_kernel_tests(kx, x, d1, d2)
        cross_kernel_tests(k′, BlockData([X]), BlockData([d1]), BlockData([d2]))
        cross_kernel_tests(kx, BlockData([x]), BlockData([d1]), BlockData([d2]))

        show(IOBuffer(), k′)
    end

    # Tests for FiniteCrossKernel.
    let
        rng, N, N′, D = MersenneTwister(123456), 5, 7, 2
        x, x′ = randn(rng, N), randn(rng, N′)
        k, X, X′ = EQ(), ColsAreObs(randn(rng, D, N)), ColsAreObs(randn(rng, D, N′))
        k′, kx = FiniteCrossKernel(k, X, X′), FiniteCrossKernel(k, x, x′)

        @test size(k′) == (N, N′)
        @test size(k′, 1) == N
        @test size(k′, 2) == N′
        @test isstationary(k′) == false
        @test eachindex(k′, 1) == 1:N
        @test eachindex(k′, 2) == 1:N′
        @test AM(k′) == pairwise(k′, eachindex(k′, 1), eachindex(k′, 2))
        @test AM(k′) ≈ pairwise(k, X, X′)
        @test pairwise(k′, 1:N-1, 2:N) ≈ pairwise(k, X[1:N-1], X′[2:N])
        @test pairwise(k′, eachindex(k′, 1), 1:N-1) ≈ pairwise(k, X, X′[1:N-1])
        @test pairwise(k′, 2:N, eachindex(k′, 2)) ≈ pairwise(k, X[2:N], X′)

        cross_kernel_tests(k′, 1:N, 1:N, 1:N′)
        cross_kernel_tests(kx, 1:N, 1:N, 1:N′)
        cross_kernel_tests(k′, BlockData([1:N]), BlockData([1:N]), BlockData([1:N′]))
        cross_kernel_tests(kx, BlockData([1:N]), BlockData([1:N]), BlockData([1:N′]))

        # Recursion depth == 2.
        r, c = 1:N-1, 1:N′-2
        k2 = FiniteCrossKernel(k′, r, c)
        @test size(k2) == (length(r), length(c))
        @test AM(k2) == pairwise(k2, r, c)
        @test AM(k2) == pairwise(k′, r, c)
        @test pairwise(k2, 2:N-1, c) == pairwise(k′, 2:N-1, c)

        cross_kernel_tests(k2, r, r, c)
        cross_kernel_tests(k2, BlockData([r]), BlockData([r]), BlockData([c]))

        show(IOBuffer(), k′)
    end
end


text:
# using Pkg

# pkg"add SparseArrays, DifferentialEquations, CSV, DataFrames"
# pkg"add FileIO, MeshIO, Makie"

using CSV, DataFrames, SparseArrays

using GeometryBasics

using CUDA
using FileIO, MeshIO
using Makie # doesn't work if no GPU :(

# citations
# "Human Femur" (https://skfb.ly/6ursH) by Eric Bauer is licensed under Creative Commons Attribution (http://creativecommons.org/licenses/by/4.0/).

# https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.2002811
# Citation: Murphy AC, Muldoon SF, Baker D, Lastowka A, Bennett B, Yang M, et al. (2018) Structure, function, and control of the human musculoskeletal network. PLoS Biol 16(1): e2002811. https://doi.org/10.1371/journal.pbio.2002811

# not adj mat b/c edges can connect more than two nodes
function get_inc_mats()
    df, df2 = CSV.read.(["data/s9.csv", "data/s10.csv"], DataFrame)
    v1, v2 = df[1:end, 3:end], df2[1:end, 3:end]
    m1, m2 = Matrix.([v1, v2])
    replace!.([m1, m2], missing=>0, "1"=>1, " "=>0)
    bm1 = Bool.(m1)
    bm2 = Bool.(m2)
    sparse.([bm1, bm2])
end

m, m2 = get_inc_mats()

x = load("data/Femur.stl") # 173 MB

typeof(x)

plot(x)


text:
# # Eady turbulence example
#
# In this example, we initialize a random velocity field and observe its viscous,
# turbulent decay in a two-dimensional domain. This example demonstrates:
#
#   * How to use a tuple of turbulence closures
#   * How to use biharmonic diffusivity
#   * How to implement a background flow (a background geostrophic shear)

using Random, Printf
using Oceananigans, Oceananigans.Diagnostics, Oceananigans.OutputWriters,
      Oceananigans.AbstractOperations, Oceananigans.Utils

# These imports from Oceananigans's `TurbuleneClosures` module are needed for imposing
# a background flow as a user-defined forcing.

using Oceananigans.Operators: ∂xᶠᵃᵃ, ∂xᶜᵃᵃ, ℑxᶠᵃᵃ, ℑyᵃᶜᵃ, ℑxᶜᵃᵃ, ℑxzᶠᵃᶜ

using Oceananigans.Fields: Face, Cell

# # Parameters
#
# Here we define numerical and physical parameters for an Eady problem.
# We use a very fast baroclinic growth rate for the sake of this coarse 3D
# example.

Nh = 64           # horizontal resolution
Nz = 32           # vertical resolution
Lh = 2e6          # [meters] horizontal domain extent
Lz = 1e3          # [meters] vertical domain extent
Rᵈ = Lh / 10      # [m] Deformation radius
σᵇ = 0.1day       # [s] Growth rate for baroclinic instability
τᵏ = 1.0day       # [s] biharmonic / viscous damping timescale
 μ = 1/30day      # [s⁻¹] linear drag decay scale
 f = 1e-4         # [s⁻¹] Coriolis parameter

 Δh = Lh / Nh     # [meters] horizontal grid spacing for diffusivity calculations
 Δz = Lz / Nz     # [meters] vertical grid spacing for diffusivity calculations
κ₄h = Δh^4 / τᵏ   # [m⁴ s⁻¹] Biharmonic horizontal diffusivity
 κᵥ = Δz^2 / 20τᵏ # [m² s⁻¹] Laplacian vertical diffusivity

@show N² = (Rᵈ * f / Lz)^2      # [s⁻²] Initial buoyancy gradient
@show  α = sqrt(N²) / (f * σᵇ)  # [s⁻¹] background shear

end_time = 3day # Simulation end time

# # Boundary conditions
#
# These boundary conditions prescribe a linear drag at the bottom as a flux
# condition. We also fix the surface and bottom buoyancy to enforce a buoyancy
# gradient `N²`.

grid = RegularCartesianGrid(size=(Nh, Nh, Nz), halo=(2, 2, 2),
                            x=(-Lh/2, Lh/2), y=(-Lh/2, Lh/2), z=(-Lz, 0))

bc_parameters = (μ=μ, H=Lz)

@inline τ₁₃_linear_drag(i, j, grid, time, iter, U, C, p) = @inbounds p.μ * p.H * U.u[i, j, 1]
@inline τ₂₃_linear_drag(i, j, grid, time, iter, U, C, p) = @inbounds p.μ * p.H * U.v[i, j, 1]

ubcs = UVelocityBoundaryConditions(grid, bottom = BoundaryCondition(Flux, τ₁₃_linear_drag))
vbcs = VVelocityBoundaryConditions(grid, bottom = BoundaryCondition(Flux, τ₂₃_linear_drag))
bbcs = TracerBoundaryConditions(grid,    top = BoundaryCondition(Value, 0),
                                      bottom = BoundaryCondition(Value, -N² * Lz))

# # Forcing functions
#
# The Eady problem is non-linearized around a geostrophic flow
#
# $ ψ = - α y (z + H) $,
#
# where α is the geostrophic shear and horizontal buoyancy gradient.
# The background buoyancy associated with this geostrophic flow is
#
# $ b = f ∂_z ψ = - α f y $
#
# To impose this background flow with user-defined forcing functions we
# require parameters that specify the veritcal shear, the Coriolis parameter
# and the domain depth.

forcing_parameters = (α=α, f=f, H=Lz)

# The $x$-momentum forcing
#
# $ F_u = - α w - α (z + H) ∂ₓu $
#
# is applied at location `(f, c, c)`.

Fu_eady(i, j, k, grid, time, U, C, p) = @inbounds (- p.α * ℑxzᶠᵃᶜ(i, j, k, grid, U.w)
                                                   - p.α * (grid.zC[k] + p.H) * ∂xᶠᵃᵃ(i, j, k, grid, ℑxᶜᵃᵃ, U.u))

# The $y$-momentum forcing
#
# $ F_v = - α (z + H) ∂ₓv
#
# is applied at location `(c, f, c)`.

Fv_eady(i, j, k, grid, time, U, C, p) = @inbounds -p.α * (grid.zC[k] + p.H) * ∂xᶜᵃᵃ(i, j, k, grid, ℑxᶠᵃᵃ, U.v)

# The $z$-momentum forcing
#
# $ F_w = - α (z + H) ∂ₓw
#
# is applied at location `(c, c, f)`.

Fw_eady(i, j, k, grid, time, U, C, p) = @inbounds -p.α * (grid.zF[k] + p.H) * ∂xᶜᵃᵃ(i, j, k, grid, ℑxᶠᵃᵃ, U.w)

# The buoyancy forcing
#
# $ F_b = - α (z + H) ∂ₓb + α f v $
#
# is applied at location `(c, c, c)`.

Fb_eady(i, j, k, grid, time, U, C, p) = @inbounds (- p.α * (grid.zC[k] + p.H) * ∂xᶜᵃᵃ(i, j, k, grid, ℑxᶠᵃᵃ, C.b)
                                                   + p.f * p.α * ℑyᵃᶜᵃ(i, j, k, grid, U.v))

# # Turbulence closures
#
# We use a horizontal biharmonic diffusivity and a Laplacian vertical diffusivity
# to dissipate energy in the Eady problem.
# To use both of these closures at the same time, we set the keyword argument
# `closure` a tuple of two closures. Note that the "2D Leith" parameterization may
# also be a sensible choice to pair with a Laplacian vertical diffusivity for this problem.

closure = (ConstantAnisotropicDiffusivity(νh=0, κh=0, νv=κᵥ, κv=κᵥ),
           AnisotropicBiharmonicDiffusivity(νh=κ₄h, κh=κ₄h))
           #TwoDimensionalLeith())

# Form a prefix from chosen resolution, boundary condition, and closure name

output_filename_prefix = string("eady_turb_Nh", Nh, "_Nz", Nz)

# # Model instantiation
#
# Our use of biharmonic diffusivity means we must instantiate the model grid
# with halos of size `(2, 2, 2)` in `(x, y, z)`.

model = Model( grid = grid,
       architecture = CPU(),
           coriolis = FPlane(f=f),
           buoyancy = BuoyancyTracer(), tracers = :b,
            forcing = ModelForcing(u=Fu_eady, v=Fv_eady, w=Fw_eady, b=Fb_eady),
            closure = closure,
boundary_conditions = SolutionBoundaryConditions(grid, u=ubcs, v=vbcs, b=bbcs),
# "parameters" is a NamedTuple of user-defined parameters that can be used in boundary condition and forcing functions.
         parameters = merge(bc_parameters, forcing_parameters))

# # Initial conditions
#
# For initial conditions we impose a linear stratifificaiton with some
# random noise.

## A noise function, damped at the boundaries
Ξ(z) = rand() * z/Lz * (z/Lz + 1)

## Buoyancy: linear stratification plus noise
b₀(x, y, z) = N² * z + 1e-2 * Ξ(z) * (N² * Lz + α * f * Lh)
u₀(x, y, z) = 1e-2 * α * Lz

set!(model, u=u₀, v=u₀, b=b₀)

# # Diagnostics and output
#
# Diagnostics that return the maximum absolute value of `u, v, w` by calling
# `umax(model), vmax(model), wmax(model)`:

umax = FieldMaximum(abs, model.velocities.u)
vmax = FieldMaximum(abs, model.velocities.v)
wmax = FieldMaximum(abs, model.velocities.w)

# Set up output. Here we output the velocity and buoyancy fields at intervals of one day.
fields_to_output = merge(model.velocities, (b=model.tracers.b,))
output_writer = JLD2OutputWriter(model, FieldOutputs(fields_to_output);
                                 interval=day, prefix=output_filename_prefix,
                                 force=true, max_filesize=10GiB)

# # The `TimeStepWizard`
#
# The TimeStepWizard manages the time-step adaptively, keeping the CFL close to a
# desired value.

wizard = TimeStepWizard(cfl=0.05, Δt=20.0, max_change=1.1, max_Δt=min(1/10f, σᵇ/10))

# # Vertical vorticity and divergence for plotting purposes
#
# We also create objects for computing the vertical vorticity and divergence
# for plotting purposes.

u, v, w = model.velocities
ζ = Field(Face, Face, Cell, model.architecture, model.grid)
δ = Field(Cell, Cell, Cell, model.architecture, model.grid)

vertical_vorticity = Computation(∂x(v) - ∂y(u), ζ)
        divergence = Computation(-∂z(w), δ)

ζmax = FieldMaximum(abs, ζ)
δmax = FieldMaximum(abs, δ)

# # Plot setup
#
# Set `makeplot = true` to make plots as the model runs.

makeplot = false

using PyPlot, PyCall

GridSpec = pyimport("matplotlib.gridspec").GridSpec

fig = figure(figsize=(12, 8))

gs = GridSpec(2, 2, height_ratios=[2, 1])

axs = ntuple(4) do i
    fig.add_subplot(get(gs, i-1))
end

function makeplot!(axs, model)
    nx, ny, nz = size(model.grid)

    xC_xy = repeat(reshape(model.grid.xC, nx, 1), 1, ny)
    xF_xy = repeat(reshape(model.grid.xF[1:end-1], nx, 1), 1, ny)

    yC_xy = repeat(reshape(model.grid.yC, 1, ny), nx, 1)
    yF_xy = repeat(reshape(model.grid.yF[1:end-1], 1, ny), nx, 1)

    xC_xz = repeat(reshape(model.grid.xC, nx, 1), 1, nz)
    xF_xz = repeat(reshape(model.grid.xF[1:end-1], nx, 1), 1, nz)

    zC_xz = repeat(reshape(model.grid.zC, 1, nz), nx, 1)
    zF_xz = repeat(reshape(model.grid.zF[1:end-1], 1, nz), nx, 1)

    compute!(vertical_vorticity)
    compute!(divergence)

    @printf("\nmax ζ/f: %.2e, max δ/f: %.2e\n\n", ζmax()/f, δmax()/f)

    sca(axs[1]); cla()
    pcolormesh(xF_xy/1e3, yF_xy/1e3, Array(interior(ζ)[:, :, Nz]))
    title("Surface vertical vorticity (1/s)"); xlabel("\$ x \$ (km)"); ylabel("\$ y \$ (km)")
    axs[1].set_aspect(1)

    sca(axs[2]); cla()
    pcolormesh(xC_xy/1e3, yC_xy/1e3, Array(interior(δ)[:, :, Nz]))
    title("Surface divergence (1/s)"); xlabel("\$ x \$ (km)"); ylabel("\$ y \$ (km)")
    axs[2].set_aspect(1)

    sca(axs[3]); cla()
    pcolormesh(xF_xz/1e3, zC_xz, Array(interior(ζ)[:, Int(Nh/2), :]))
    title("Vertical vorticity (1/s)"); xlabel("\$ x \$ (km)"); ylabel("\$ z \$ (m)")

    sca(axs[4]); cla()
    pcolormesh(xC_xz/1e3, zF_xz, Array(interior(w)[:, Int(Nh/2), :]))
    title("Vertical velocity (m/s)"); xlabel("\$ x \$ (km)"); ylabel("\$ z \$ (m)")

    tight_layout()
    pause(0.1)

    return nothing
end

# # Plot setup
#
# This time-stepping loop runs until end_time is reached. It prints a progress statement
# every 100 iterations.

while model.clock.time < end_time

    ## Update the time step associated with `wizard`.
    update_Δt!(wizard, model)

    ## Time step the model forward
    walltime = Base.@elapsed time_step!(model, 10, wizard.Δt)

    ## Print a progress message
    @printf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
            model.clock.iteration, prettytime(model.clock.time), prettytime(wizard.Δt),
            umax(), vmax(), wmax(), prettytime(walltime))

    if model.clock.iteration % 100 == 0 && makeplot
        makeplot!(axs, model)
    end
end

# Make a plot at the end

makeplot!(axs, model)
gcf()


text:
include("./test_file_specs.jl")

@testset "Data Types" begin
end