text:
module Demo
using PolygonInbounds
using PolygonInbounds: PolygonMesh, all_areas, edges_dict_for_area, vertex

export polydemo
struct Dummy end

Base.getproperty(::Dummy,::Symbol) = dummy(a...;kw...) = nothing
setplot(x) = global Plot = x

Plot = Dummy()
plot(a...;kw...) = Plot.plot(a...;kw...)
plot!(a...;kw...) = Plot.plot!(a...;kw...)
scatter(a...;kw...) = Plot.scatter(a...;kw...)
scatter!(a...;kw...) = Plot.scatter!(a...;kw...)
savefig(a...) = Plot.savefig(a...)
display(a...) = Plot.display(a...)

const PLOTARG = (legend = :topleft, size = (1024, 780))

"""
%POLYDEMO run a series of point-in-polygon demos.
%   POLYDEMO(II) runs the II-th demo, where +1 <= II <= +3.
%   Demo problems illustrate varying functionality for the
%   INPOLY2 routine.
%
%   See also INPOLY2, INPOLYGON

%   Darren Engwirda : 2018 --
%   Email           : darren.engwirda@columbia.edu
%   Last updated    : 28/10/2018
"""
function polydemo(id::Integer=1; args...)
    [demo1, demo2, demo3, demo4, demo5][id](;args...)
end

#-----------------------------------------------------------
function demo1(;tol=1e-2, r=10^5, do_plot=true)
    do_plot && println("""
        INPOLY2 provides fast point-in-polygon queries for ob-\n
        jects in R^2. Like INPOLYGON, it detects points inside\n
        and on the boundary of polygonal geometries.\n
""")

    nodes = [ 2.0 0; 8 4; 4 8; 0 4; 5 0; 7 4; 3 6; 6 1 # inner nodes
             3 3; 5 3; 5 5; 3 5 ]                     # outer nodes
    edges = [ 1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1
             9 10; 10 11; 11 12; 12 9 ]

    poly = PolygonMesh(nodes, edges)
    dx = 10 * sqrt(1/r)
    a = [(x,y) for x in -1:dx:9 for y in -1:dx:9]
    xpos, ypos = [p[1] for p in a], [p[2] for p in a]

    stat = @time inpoly2([xpos ypos], nodes, edges, atol=tol)
    if do_plot
        p = plot(title="demo1 (r $r, atol $tol)"; PLOTARG...)
        p = plotpoints(p, xpos, ypos, stat)
        p = plotpolygon(p, poly)
        savefig(p, "demo1.png")
        display(p)
    end
    stat
end

#-----------------------------------------------------------
function demo2(;r=2500, tol=1e-3, do_plot=true)
    do_plot && println("""
        INPOLY2 supports multiply-connected geometries, consi-\n
        sting of arbitrarily nested sets of outer + inner bou-\n
        ndaries.\n
""")

    xpos, ypos, nodes, edges = testdata("lakes.msh", r)
    poly = PolygonMesh(nodes, edges)
    stat = @time inpoly2([xpos ypos], nodes, edges, rtol=tol)
    if do_plot
        p = plot(title="demo2 - lakes (r $r, rtol $tol)"; PLOTARG...)
        p = plotpoints(p, xpos, ypos, stat)
        p = plotpolygon(p, poly)
        savefig(p, "demo2.png")
        display(p)
    end
    stat
end

#-----------------------------------------------------------
function demo3(;r=2500, tol=1e-3, do_plot=true)
    do_plot && println("""
        INPOLY2 implements a "pre-sorted" variant of the cros-\n
        sing-number test - returning queries in approximately \n
        O((N+M)*LOG(N)) time for configurations consisting of \n
        N test points and M edges. This is often considerably \n
        faster than conventional approaches that typically re-\n
        quire O(N*M) operations.\n
        """);

    xpos, ypos, nodes, edges = testdata("coast.msh", r)
    poly = PolygonMesh(nodes, edges)

    stat = @time inpoly2([xpos ypos], nodes, edges, rtol=tol)
    if do_plot
        p = plot(title="demo3 - coast (r $r, rtol $tol)"; PLOTARG...)
        p = plotpoints(p, xpos, ypos, stat)
        p = plotpolygon(p, poly)
        savefig(p, "demo3.png")
        display(p)
    end
    stat
end

#-----------------------------------------------------------
function demo4(;r=5*10^5, tol=1e-2, do_plot=true)
    do_plot && println("""
        INPOLY2 provides fast point-in-polygon queries for ob-\n
        jects in R^2. Like INPOLYGON, it detects points inside\n
        and on the boundary of polygonal geometries.\n
    """)

    nodes = [0.0 0; 3 4; 7 5; 10 0]
    edges = [1 2; 2 3; 3 4; 4 1]
    poly = PolygonMesh(nodes, edges)

    rpts = testbox(r, [-1.0 -1.0], [ 11.0 6.0])
    xpos, ypos = rpts[:,1], rpts[:,2]

    stat = @time inpoly2(rpts, nodes, edges, atol=tol)
    if do_plot
        p = plot(title="demo4 (r $r, atol $tol)"; PLOTARG...)
        p = plotpoints(p, xpos, ypos, stat)
        p = plotpolygon(p, poly)
        savefig(p, "demo4.png")
        display(p)
    end
    stat
end

#-----------------------------------------------------------
function demo5(;r=5*10^5, tol=1e-2, do_plot=true)
    do_plot && println("""
        INPOLY2 provides fast point-in-polygon queries for ob-\n
        jects in R^2. Like INPOLYGON, it detects points inside\n
        and on the boundary of polygonal geometries.\n
    """)

    nodes = [0.0 0; 3 4; 7 5; 10 0; 11 4]
    edges = [1 2 1; 2 3 1; 3 4 1; 4 1 1; 3 4 2; 4 5 2; 5 3 2]
    poly = PolygonMesh(nodes, edges)

    rpts = testbox(r, [-1.0 -1.0], [ 11.0 6.0])
    xpos, ypos = rpts[:,1], rpts[:,2]

    stat = @time inpoly2(rpts, nodes, edges, atol=tol)
    if do_plot
        p = plot(title="demo5 (r $r, atol $tol)"; PLOTARG...)
        p = plotpoints(p, xpos, ypos, stat)
        p = plotpolygon(p, poly)
        savefig(p, "demo5.png")
        display(p)
    end
    stat
end

#----- Utility functions

function testdata(dataset, r::Integer)
    filepath = @__DIR__
    nodes, edges = loadmsh(joinpath(filepath, "..", "test-data", dataset))

    emid = .5 * nodes[edges[:,1],:] + .5 * nodes[edges[:,2],:]

    ma, mi =  maximum(nodes, dims=1), minimum(nodes, dims=1)
    rpts = testbox(r, mi, ma)

    vert = r <= 100 ? [rpts;] : [nodes; emid; rpts]
    vert[:,1], vert[:,2], nodes, edges
end

function testbox(r::Integer, mi, ma)
    mi = reshape(mi, 1, 2)
    ma = reshape(ma, 1, 2)
    half = (ma + mi) / 2
    scal = ma - mi
    rpts = (rand(r, 2) .- 0.5).* scal .* 1.1 .+ half
end

function plotpolygon(p, poly::PolygonMesh)
    areas = all_areas(poly)
    T = Float64

    for a in areas

        d = edges_dict_for_area(poly, a)
        while !isempty(d)
            start = j = minimum(keys(d))
            n = length(d)
            xpos = Vector{T}(undef, n+1)
            ypos = Vector{T}(undef, n+1)
            k = 0
            while k <= n
                k += 1
                xpos[k] = vertex(poly, j,1)
                ypos[k] = vertex(poly, j,2)
                j == start && k > 1 && break
                j = pop!(d, j)
            end
            resize!(xpos, k)
            resize!(ypos, k)
            p = plot!(p, xpos, ypos, color=:purple, width=1.5)
        end
    end
    p
end

function plotpoints(p, xpos, ypos, stat)
    function sca(ina, color=:black, mshape=:cross)
        scatter!(p, xpos[ina], ypos[ina]; markershape=mshape, markersize=ms, color=color)
    end

    ms = 1.0
    inside = onbound = nothing
    for area in axes(stat,3)
        inside = stat[:,1,area]
        onbound = stat[:,2,area]
        ina = inside .& (!).(onbound)
        inb = inside .& onbound

        sca(ina, area == 1 ? :green : :brown)
        sca(inb, :yellow)
    end
    if size(stat,3) > 1
        inside = BitVector([any(view(stat,i,1,:)) for i in axes(stat,1)])
        onbound = BitVector([any(view(stat,i,2,:)) for i in axes(stat,1)])
    end
    outa = (!).(inside) .& (!).(onbound)
    outb = (!).(inside) .& onbound
    sca(outb, :blue)
    sca(outa, :lightblue)

    p
end

function loadmsh(file)
    nodes = Float64[]
    edges = Int[]
    m = 0
    n = 0
    open(file) do io
        name =""
        while !eof(io)
            line = readline(io)
            startswith(line, "#") && continue
            ls = split(line, "=")
            if length(ls) == 2
                name = ls[1]
                if name == "point"
                    n = parse(Int, ls[2])
                    nodes = Vector{Float64}(undef, 2n)
                    m = 0
                elseif name == "edge2"
                    n = parse(Int, ls[2])
                    edges = Vector{Int}(undef, 2n)
                    m = 0
                end
                continue
            end
            ls = split(line, ";")
            if length(ls) >= 2
                m >= 2n - 1 && continue
                if name == "point"
                    nodes[m+=1] = parse(Float64, ls[1])
                    nodes[m+=1] = parse(Float64, ls[2])
                elseif name == "edge2"
                    edges[m+=1] = parse(Int, ls[1]) + 1
                    edges[m+=1] = parse(Int, ls[2]) + 1
                end
            end
        end
    end
    res(nodes) = permutedims(reshape(nodes, 2, length(nodes)÷2))
    res(nodes), res(edges)
end

end # Demo


text:
using Gadfly
# using Sigma
# REVIEW - Decide whether plotting functions should return plot or layer
# REVIEW - Add 3D preimages. Important, not just because the look cool
#        - but because they will help with intuitions about independence
#        - and symmetry

function plot_preimage(X::RandVar{Bool};prob_args...)
  satsets, mixedsets = Sigma.pre_bfs(X, t, LazyOmega();prob_args...)
  plot_preimage(satsets)
end

function plot_preimage(X::RandVar{Bool}, dims::Vector; prob_args...)
  satsets, mixedsets = Sigma.pre_bfs(X, t, LazyOmega();prob_args...)
  plot_preimage(satsets, dims)
end

function plot_preimage(os::Vector{Omega{Interval}})
  intervals = Vector{Interval}[convert(Vector{Interval},o) for o in os]
  plot_2d_boxes(intervals)
end

function plot_preimage(os::Vector{Omega{Interval}}, dims::Vector)
  intervals = Vector{Interval}[convert(Vector{Interval},o, dims) for o in os]
  @show typeof(intervals)
  plot_2d_boxes(intervals)
end

# Returns mathematica function
function plot_preimage3D(X::RandVar{Bool};prob_args...)
  satsets, mixedsets = Sigma.pre_bfs(X, t, LazyOmega();prob_args...)
  plot_preimage3D(satsets,mixedsets)
end

function plot_preimage3D(satos::Vector{Omega{Interval}}, mixedos::Vector{Omega{Interval}})
  satintervals = Vector{Interval}[convert(Vector{Interval},o) for o in satos]
  mixedintervals = Vector{Interval}[convert(Vector{Interval},o) for o in mixedos]

  satcuboids = ["Cuboid[{$(box[1].l), $(box[2].l), $(box[3].l)}, {$(box[1].u), $(box[2].u), $(box[3].u)}]" for box in satintervals]
  mixedos = ["Cuboid[{$(box[1].l), $(box[2].l), $(box[3].l)}, {$(box[1].u), $(box[2].u), $(box[3].u)}]" for box in mixedintervals]
  "Graphics3D[{Green, Opacity[0.3],$(join(satcuboids, ",")),Blue, Opacity[1.0],$(join(mixedos, ","))},{Axes->True},{PlotRange -> {{0, 1}, {0, 1}, {0, 1}}}, {AxesLabel->{ω1,ω2,ω3}}, {ViewPoint -> Dynamic[vp]}]"
end

# Take as argument canvas range
function plot_2d_boxes(bs::Vector{Vector{Interval}};
                       abs_x_min = 0.0, abs_y_min = 0.0,
                       abs_x_max = 1.0, abs_y_max = 1.0)
  x_min = Float64[]
  y_min = Float64[]
  x_max = Float64[]
  y_max = Float64[]

  for b in bs
    push!(x_min, b[1].l)
    push!(x_max, b[1].u)
    push!(y_min, b[2].l)
    push!(y_max, b[2].u)
  end

  # Update absolumte min/max values for canvas
  abs_x_min = min(abs_x_min, minimum(x_min))
  abs_y_min = min(abs_y_min, minimum(y_min))
  abs_y_max = min(abs_y_max, maximum(y_max))
  abs_x_max = min(abs_x_max, maximum(x_max))

  plot(x_min=x_min, x_max=x_max, y_min=y_min, y_max=y_max,
       color=linspace(0,1,length(bs)),
       Theme(line_width=10.0mm,highlight_width=20.0mm),
       Geom.rectbin,
       Scale.x_continuous(minvalue=abs_x_min, maxvalue=abs_y_max),
       Scale.y_continuous(minvalue=abs_y_min, maxvalue=abs_y_max),
       Coord.cartesian(fixed=true))
end

# FIXME, REVIEW: This doesn't seem to be working
function plot_volume_hist{T}(os::Vector{Omega{T}})
  vols = Float64[Sigma.measures(o) for o in os]
  plot(x=vols, Geom.histogram, Scale.y_log10)
end

# FIXME, REVIEW: This doesn't seem to be working
function plot_volume_hist(t::Sigma.Tree)
  sat_nodes::Vector{Omega} = map(n->n.data,t.nodes)
  volumes = measures(sat_nodes)
  plot(x = volumes, Geom.histogram, Scale.y_log10)
end

function plot_density(X::RandVar, lower::Real, upper::Real;
                      n_bars = 20, progargs...)
  xs = Array{Float64}(n_bars - 1)
  ys = Array{Float64}(n_bars - 1)
  ymin = Array{Float64}(n_bars - 1)
  ymax = Array{Float64}(n_bars - 1)
  for i in 1:(n_bars - 1)
    j = i + 1
    l,u = linspace(lower,upper,n_bars)[i],linspace(lower,upper,n_bars)[j]
    prob_bounds = prob((X > l) & (X <= u); progargs...)
    ys[i] = Sigma.mid(prob_bounds)
    ymin[i] = prob_bounds.l
    ymax[i] = prob_bounds.u
    xs[i] = (u-l)/2 + l
  end
  plot(x=xs, y=ys, Scale.x_continuous(minvalue=lower, maxvalue=upper),
       Scale.y_continuous(minvalue=0),
       ymin = ymin, ymax = ymax, Geom.line, Geom.errorbar)
end

# REVIEW - DRY
function plot_cond_density(X::RandVar, Y::RandVar, lower::Real, upper::Real;
                           n_bars = 20, progargs...)
  xs = Array{Float64}(n_bars - 1)
  ys = Array{Float64}(n_bars - 1)
  ymin = Array{Float64}(n_bars - 1)
  ymax = Array{Float64}(n_bars - 1)
  for i in 1:(n_bars - 1)
    j = i + 1
    l,u = linspace(lower,upper,n_bars)[i],linspace(lower,upper,n_bars)[j]
    prob_bounds = cond_prob((X > l) & (X <= u),Y; progargs...)
    ys[i] = Sigma.mid(prob_bounds)
    ymin[i] = prob_bounds.l
    ymax[i] = prob_bounds.u
    xs[i] = (u-l)/2 + l
  end
  plot(x=xs, y=ys, Scale.x_continuous(minvalue=lower, maxvalue=upper),
       Scale.y_continuous(minvalue=0),
       ymin = ymin, ymax = ymax, Geom.line, Geom.errorbar)
end

# REVIEW - Test,
function plot_sample_density(X::RandVar, nsamples::Int64)
  samples = [rand(X) for i = 1:nsamples]
  Gadfly.plot(x = samples, Gadfly.Geom.density, Gadfly.Geom.Scale.y_continuous(minvalue=0))
end

function plot_sample_histogram(X::RandVar, nsamples::Int64)
  samples = [rand(X) for i = 1:nsamples]
  Gadfly.plot(x = samples, Gadfly.Geom.histogram, Gadfly.Geom.Scale.y_continuous(minvalue=0))
end

function plot_sample_cond_histogram(X::RandVar, Y::RandVar{Bool}, nsamples::Int64; probargs...)
  condX = conditional(X,Y; probargs...)
  samples = [rand(condX) for i = 1:nsamples]
  Gadfly.plot(x = samples, Gadfly.Geom.histogram, Gadfly.Geom.Scale.y_continuous(minvalue=0))
end

function plot_sample_cond_density(X::RandVar, Y::RandVar{Bool}, nsamples::Int64; probargs...)
  condX = conditional(X,Y; probargs...)
  samples = [rand(condX) for i = 1:nsamples]
  Gadfly.plot(x = samples, Gadfly.Geom.density, Gadfly.Geom.Scale.y_continuous(minvalue=0))
end

function plot_sat_distribution(t::Sigma.Tree)
  sat_nodes = filter(x->x.status==SAT,t.nodes)
  unsat_nodes = filter(x->x.status==UNSAT,t.nodes)
  mixed_nodes = filter(x->x.status==PARTIALSAT,t.nodes)
  println("Relative counts of SAT UNSAT AND MIXED SAT are: ", map(length, Array[sat_nodes, unsat_nodes, mixed_nodes]))
  plot(y = map(length, Array[sat_nodes, unsat_nodes, mixed_nodes]), x = ["SAT", "UNSAT", "MIXED"], Geom.bar)
end

ext: jl


text:
export Tile, render

@doc """
A `Tile` is the basic currency in Canvas.
Most of the functions in the Canvas API take `Tile`s
among other things as arguments, and return a `Tile` as the result.

Tiles are immutable: once created there is no way to mutate them.
""" ->
abstract Tile

render{T <: Tile}(x::T) =
    error("$T cannot be rendered.")

immutable Leaf <: Tile
    element::Elem
end
render(x::Elem) = x
render(x::Leaf) = x.element

convert(::Type{Tile}, x::String) = Leaf(Elem(:span, x))
convert{ns, tag}(::Type{Tile}, x::Elem{ns, tag}) = Leaf(x)

function bestmime(val)
  for mime in ("text/html", "image/svg+xml", "image/png", "text/plain")
    mimewritable(mime, val) && return MIME(symbol(mime))
  end
  error("Cannot render $val.")
end

render_fallback(m::MIME"text/plain", x) = Elem(:div, stringmime(m, x))
render_fallback(m::MIME"text/html", x)  = Elem(:div, innerHTML=stringmime(m, x))
render_fallback(m::MIME"text/svg", x)   = Elem(:div, innerHTML=stringmime(m, x))
render_fallback(m::MIME"image/png", x)  = Elem(:img, src="data:image/png;base64," * stringmime(m, x))

render(x::FloatingPoint) = @sprintf "%0.3f" x
render(x::Symbol) = string(x)
render(x::String) = Elem(:span, x)
render(xs::AbstractArray, tag="div") = Elem(tag, map(render, xs))

immutable AnyWrap <: Tile
    value
end

render{T}(x::T) =
    # Try to convert first
    method_exists(convert, (Type{Tile}, T)) ?
        render(convert(Tile, x)) :
        render(AnyWrap(x))

# Catch-all render
function render(x::AnyWrap)
    render_fallback(bestmime(x.value), x.xvalue)
end

@doc """
`Empty` is handy tile that is well... Empty.
use `empty` constant exported by Canvas in your code.
""" ->
immutable Empty <: Tile
end

@doc """
An Empty element.
""" ->
const empty = Empty()

render(t::Empty) = Elem(:div)

writemime(io::IO, m::MIME"text/html", x::Tile) =
    writemime(io, m, Elem(:div, Canvas.render(x), className="canvasRoot"))

writemime{T <: Tile}(io::IO, m::MIME"text/html", x::Signal{T}) =
    writemime(io, m, lift(Canvas.render, Patchwork.Elem, x))

render{T <: Tile}(s::Signal{T}) =
    render(value(s))

# Note a TileList is NOT a Tile
immutable TileList
    tiles::AbstractArray
end

convert(::Type{TileList}, xs::AbstractArray) =
    TileList(xs)
convert(::Type{TileList}, xs::Tuple) =
    TileList([x for x in xs])
convert(::Type{TileList}, x) =
    TileList([x])

render(t::TileList) =
    map(render, t.tiles)

render(t::TileList, wrap) =
    Elem(wrap, map(render, t.tiles))


text:
export simulate_sibship

"""
    _cross(parent1::Vector{Vector{T}}, parent2::Vector{Vector{T}}) where T <: Signed
Simulate a mating cross between two parents, generating one offspring with the same
ploidy as `parent1`. This variant of `cross` is used internally for `simulate_sibship`.
"""
function _cross(parent1::Vector{Vector{T}}, parent2::Vector{Vector{T}}) where T <: Signed
    ploidy = length(first(parent1))
    ploidy == 1 && error("Haploid crosses are not yet supported. Please file and issue or pull request")
    if ploidy == 2
        p1_contrib = rand.(parent1)
        p2_contrib = rand.(parent2)
        geno_out = sort.(zip(p1_contrib, p2_contrib))
    elseif iseven(ploidy)
        n_allele = ploidy ÷ 2
        p1_contrib = sample.(parent1, n_allele, replace = false)
        p2_contrib = sample.(parent2, n_allele, replace = false)
        geno_out = Tuple.(sort!.(append!.(p1_contrib, p2_contrib)))
    else
        # special method to provide a 50% chance of one parent giving more alleles than the other
        rng = rand()
        contrib_1 = ploidy ÷ 2
        contrib_2 = ploidy - contrib_1
        p1_contrib = rng > 0.5 ? sample.(parent1, contrib_1, replace = false) : sample.(parent1, contrib_2, replace = false)
        p2_contrib = rng > 0.5 ? sample.(parent2, contrib_2, replace = false) : sample.(parent2, contrib_1, replace = false)
        geno_out = Tuple.(sort!.(append!.(p1_contrib, p2_contrib)))
    end
    return geno_out
end

function parentoffspring(data::PopData; n::Int, ploidy::Int)
    loc, alleles = allele_pool(data)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim$i"
        p1,p2 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:2]
        insertcols!(out_df, Symbol(prefix * "_parent") =>  Tuple.(p1))
        insertcols!(out_df, Symbol(prefix * "_offspring") => _cross(p1, p2))
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "parent_offspring")
    out_df.name = PooledArray(out_df.name)
    out_df.population = PooledArray(out_df.population)
    out_df.locus = PooledArray(out_df.locus)

    meta_df = DataFrame(
        :name => unique(out_df.name),
        :population => "parent_offspring",
        :ploidy => ploidy,
        :longitude => Vector{Union{Missing, Float32}}(undef, n*2),
        :latitude => Vector{Union{Missing, Float32}}(undef, n*2),
        )
    return PopData(meta_df, out_df)
end

function fullsib(data::PopData; n::Int, ploidy::Int)
    loc, alleles = allele_pool(data)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim$i"
        p1,p2 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:2]
        [insertcols!(out_df, Symbol(prefix * "_fs_off$j") => _cross(p1, p2)) for j in 1:2]
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "fullsib")
    out_df.name = PooledArray(out_df.name)
    out_df.population = PooledArray(out_df.population)
    out_df.locus = PooledArray(out_df.locus)

    meta_df = DataFrame(
        :name => unique(out_df.name),
        :population => "fullsib",
        :ploidy => ploidy,
        :longitude => Vector{Union{Missing, Float32}}(undef, n*2),
        :latitude => Vector{Union{Missing, Float32}}(undef, n*2),
        )
    return PopData(meta_df, out_df)
end

function halfsib(data::PopData; n::Int, ploidy::Int)
    loc, alleles = allele_pool(data)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim$i"
        p1,p2,p3 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:3]
        insertcols!(out_df, Symbol(prefix * "_hs_off1") => _cross(p1, p2))
        insertcols!(out_df, Symbol(prefix * "_hs_off2") => _cross(p1, p3))
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "halfsib")
    out_df.name = PooledArray(out_df.name)
    out_df.population = PooledArray(out_df.population)
    out_df.locus = PooledArray(out_df.locus)

    meta_df = DataFrame(
        :name => unique(out_df.name),
        :population => "halfsib",
        :ploidy => ploidy,
        :longitude => Vector{Union{Missing, Float32}}(undef, n*2),
        :latitude => Vector{Union{Missing, Float32}}(undef, n*2),
        )
    return PopData(meta_df, out_df)
end

function unrelated(data::PopData; n::Int, ploidy::Int)
    loc, alleles = allele_pool(data)
    out_df = DataFrame(:locus => loc)
    for i in 1:n
        prefix = "sim$i"
        p1,p2 = [simulate_sample(alleles, loc, ploidy = ploidy) for j in 1:2]
        insertcols!(out_df, Symbol(prefix * "_un_off1") => Tuple.(p1))
        insertcols!(out_df, Symbol(prefix * "_un_off2") => Tuple.(p2))
    end
    out_df = rename!(select!(DataFrames.stack(out_df, Not(:locus)), 2, 1, 3), [:name, :locus, :genotype])
    insertcols!(out_df, 2, :population => "unrelated")
    out_df.name = PooledArray(out_df.name)
    out_df.population = PooledArray(out_df.population)
    out_df.locus = PooledArray(out_df.locus)

    meta_df = DataFrame(
        :name => unique(out_df.name),
        :population => "unrelated",
        :ploidy => ploidy,
        :longitude => Vector{Union{Missing, Float32}}(undef, n*2),
        :latitude => Vector{Union{Missing, Float32}}(undef, n*2),
        )
    return PopData(meta_df, out_df)
end

"""
    simulate_sibship(data::PopData; n::Int, relationship::String, ploidy::Int)
Simulate mating crosses to generate `n` sample pairs (default: `500`) having the specified `relationship`,
returning a `PopData` object. The simulations will first generate parents of a given `ploidy` (inferred or specified)
by drawing alleles from a global allele pool derived from the given `data` (i.e. weighted by their frequencies).

#### Relationship
Simulated parents will be crossed to generate offspring depending on the `relationship`:
- `"fullsib"` : 2 parents generate 2 full-sibling offspring, return 2 offspring
- `"halfsib` : 3 parents generate 2 half-sibling offspring, returns 2 offspring
- `"unrelated"` : returns 2 randomly generated individuals from the global allele pools
- `"parent-offspring"` : 2 parents generate 1 offspring, returns 1 offspring and 1 parent

#### Identifying pairs
The relationship between the newly generated samples can be identified by:
- Sample `name`s will specify their simulation number, relationship, and whether parent or offspring
- Their `population` name will be that of their relationship (e.g. "fullsib")

#### Ploidy
If the samples in your `PopData` are of a single ploidy, then `ploidy = 0` (the default) will infer the ploidy
and generate parents and offspring according to the ploidy of your data. If you have mixed-ploidy data or wish
to generate parents and offspring of a ploidy different than the source `PopData` you can specify the ploidy
with which to simulate parents and offspring. For example, if your `PopData` is diploid, but you wish to generate
triploid or octoploid parents and offspring, you would specify `ploidy = 3` or `ploidy = 8` repectively.
#### Odd ploidy
If trying to create offspring with an odd ploidy (3,5, etc.), each parent has a 50% chance of
contributing (½ × ploidy) + 1 alleles for all loci to the offspring. In other words, if ploidy = 3,
there's a 50% chance parent_1 will give 2 alleles for every locus for that simulated cross.

**Example**
```
julia> cats = nancycats() ;

julia> fullsib_sims = simulate_sibship(cats, n = 50, relationship= "fullsib")
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 100
  Loci: 9
  Populations: 1
  Longitude: absent
  Latitude: absent

julia> fullsib_sims.meta_df100×5 DataFrame
│ Row │ name          │ population │ ploidy │ longitude │ latitude │
│     │ String        │ String     │ Int64  │ Float32?  │ Float32? │
├─────┼───────────────┼────────────┼────────┼───────────┼──────────┤
│ 1   │ sim1_fs_off1  │ fullsib    │ 2      │ missing   │ missing  │
│ 2   │ sim1_fs_off2  │ fullsib    │ 2      │ missing   │ missing  │
│ 3   │ sim2_fs_off1  │ fullsib    │ 2      │ missing   │ missing  │
│ 4   │ sim2_fs_off2  │ fullsib    │ 2      │ missing   │ missing  │
│ 5   │ sim3_fs_off1  │ fullsib    │ 2      │ missing   │ missing  │
⋮
│ 95  │ sim48_fs_off1 │ fullsib    │ 2      │ missing   │ missing  │
│ 96  │ sim48_fs_off2 │ fullsib    │ 2      │ missing   │ missing  │
│ 97  │ sim49_fs_off1 │ fullsib    │ 2      │ missing   │ missing  │
│ 98  │ sim49_fs_off2 │ fullsib    │ 2      │ missing   │ missing  │
│ 99  │ sim50_fs_off1 │ fullsib    │ 2      │ missing   │ missing  │
│ 100 │ sim50_fs_off2 │ fullsib    │ 2      │ missing   │ missing  │
```
"""
function simulate_sibship(data::PopData; n::Int = 500, relationship::String = "nothing", ploidy::Int = 0)
    if relationship == "nothing"
        error("Please use the keyword \'relationship\' and specify one of: \n- \"fullsib\" \n- \"halfsib\" \n- \"unrelated\"\n- \"parent-offspring\"")
    elseif relationship ∉ ["fullsib", "halfsib", "unrelated", "parent-offspring"]
        error("relationship = \"$relationship\" is invalid, please specify one of: \n- \"fullsib\" \n- \"halfsib\" \n- \"unrelated\"\n- \"parent-offspring\"")
    end
    # automatic ploidy finding
    if ploidy == 0
        ploids = unique(data.meta.ploidy)
        if length(ploids) != 1
            error("For mixed ploidy data, please specify a single ploidy with which to generate parents and offspring")
        else
            ploidy += first(ploids)
        end
    end
    if relationship == "fullsib"
        fullsib(data, n = n, ploidy = ploidy)
    elseif relationship == "halfsib"
        halfsib(data, n = n, ploidy = ploidy)
    elseif relationship == "unrelated"
        unrelated(data, n = n, ploidy = ploidy)
    else relationship == "halfsib"
        parentoffspring(data, n = n, ploidy = ploidy)
    end
end


text:
# Utilities we need. unzip is being contributed to ZipFile,
# the more memory-efficient lda is being offered to MultivariateStats.
function unzip(inputfilename, outputpath=pwd())
    r = ZipFile.Reader(inputfilename)
    for f in r.files
        outpath = joinpath(outputpath, f.name)
        if isdirpath(outpath)
            mkpath(outpath)
        else
            open(outpath, "w") do io
                write(io, read(f))
            end
        end
    end
    nothing
end

## Memory-efficient LDA
# These are large vectors, and storing the full covariance matrix is
# problematic. This algorithm works in terms of a projection onto a
# subspace spanning the means, which also serves to regularize. See
#   Yu, Hua, and Jie Yang. "A direct LDA algorithm for
#      high-dimensional data—with application to face recognition."
#      Pattern recognition 34.10 (2001): 2067-2070.
function lda(X::AbstractMatrix{T}, group::Vector{Int}) where {T}
    nd = size(X,1)
    dmeans, dX = lda_prepare(X, group)
    UB, SB, _ = svd(dmeans, thin=true)
    ikeep = SB .>= sqrt(eps(T))*maximum(SB)
    Y = UB[:,ikeep]
    Z = scale(Y, 1 ./ SB[ikeep])
    projW = Z'*dX
    UW, SW, _ = svd(projW, thin=true)
    eigenval = 1 ./ SW.^2
    eigenvec = Z*scale(UW, 1 ./ SW)
    # Normalize
    for j = 1:size(eigenvec,2)
        evn = zero(T)
        for i = 1:nd
            evn += eigenvec[i,j]^2
        end
        evn = sqrt(evn)
        for i = 1:nd
            eigenvec[i,j] /= evn
        end
    end
    # Start with largest eigenvalues first
    return eigenvec[:,end:-1:1], eigenval[end:-1:1]
end

function lda_prepare(X::AbstractMatrix{T}, group::Vector{Int}) where {T}
    nd = size(X,1)
    npoints = size(X,2)
    xbar = zeros(T,nd)
    ngroups = maximum(group)
    if ngroups > npoints
        error("The group indices must be consecutive, starting at 1")
    end
    ngroups > 1 || error("Must have more than one group")
    # Calculate the group means
    means = zeros(T,nd,ngroups)
    npergroup = zeros(Int,ngroups)
    for j = 1:npoints
        k = group[j]
        for i = 1:nd
            means[i,k] += X[i,j]
        end
        npergroup[k] += 1
    end
    keep = npergroup .> 0
    if !all(keep)
        npergroup = npergroup[keep]
        means = means[:,keep]
    end
    gindex = cumsum(keep)
    ngroups = length(npergroup)
    for j = 1:ngroups
        n = npergroup[j]
        for i = 1:nd
            xbar[i] += means[i,j]
            means[i,j] /= n
        end
    end
    for i = 1:nd
        xbar[i] /= npoints
    end
    # Compute differences from the means
    dX = Array(T, nd, npoints)
    for j = 1:npoints
        k = gindex[group[j]]
        for i = 1:nd
            dX[i,j] = X[i,j] - means[i,k]
        end
    end
    dmeans = means .- xbar
    dmeans, dX
end


text:
if basename(pwd()) == "aoc"
    cd("2017/5")
end

function part1(filename::AbstractString)
    xs = parse.(Int, eachline(filename))
    i = 1
    total = 0
    while checkbounds(Bool, xs, i)
        x = xs[i]
        xs[i] += 1
        i += x
        total += 1
    end
    total
end

@assert part1("example.txt") == 5
part1("input.txt")

function part2(filename::AbstractString)
    xs = parse.(Int, eachline(filename))
    i = 1
    total = 0
    while checkbounds(Bool, xs, i)
        x = xs[i]
        if x >= 3
            xs[i] -= 1
        else
            xs[i] += 1
        end
        i += x
        total += 1
    end
    total
end

@assert part2("example.txt") == 10
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
include(string(pwd(),"/src/gensys.jl"))

theta = 0.36
delta = 0.025
beta = 0.99
A = 1.72
h0 = 0.58
gamma = 0.95
sig = 0.0712

B = A*log(1-h0)/h0

r_bar = 1/beta-(1-delta)

G0 = zeros(7,7)
G1 = zeros(7,7)
Pi = zeros(7,2)
Psi = zeros(7)

G0[1,1] = 1
G0[1,2] = -r_bar*beta

G0[3,5] = 1

G0[7,7] = 1

###G1####

G1[1,1] = 1

G1[2,1] = 1
G1[2,3] = 1
G1[2,4] = -1

G1[3,1] = -(r_bar/theta-delta)
G1[3,4] = r_bar/theta
G1[3,5] = (1-delta)

G1[4,3] = 1-theta
G1[4,4] = -1
G1[4,5] = theta
G1[4,7] = 1

G1[5,2] = -1
G1[5,4] = 1
G1[5,5] = -1

G1[6,3] = -1
G1[6,4] = 1
G1[6,6] = -1

G1[7,7] = gamma

Pi[1,1] = 1
Pi[1,2] = -beta*r_bar

Psi[7] = 1

sol1 = gensys(G0,G1,Psi,Pi)

irf1 = irf(sol1,100,0.01)

using Plots

plot(irf1[:,1], w = 2, label = "C")
plot!(irf1[:,2], w = 2, label = "r")
plot!(irf1[:,3], w = 2, label = "h")
plot!(irf1[:,4], w = 2, label = "y")
plot!(irf1[:,5], w = 2, label = "k")
plot!(irf1[:,6], w = 2, label = "w", line = :dash, color = "red")
hline!([0], color = "black", w = 2, label = "0")

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

G0 = zeros(7,7)
G1 = zeros(7,7)
Pi = zeros(7,2)
Psi = zeros(7)

h_bar = - (1-theta)/(B*(1-delta*theta*beta/(1-beta*(1-delta))-1))
k_bar = (theta*beta/(1-beta*(1-delta)))^(1/(1-theta))*h_bar
y_bar = r_bar*k_bar/theta
c_bar = y_bar - delta*k_bar

G0[1,1] = 1
G0[1,2] = -r_bar*beta

G0[2,1] = 1
G0[2,3] = 1
G0[2,4] = -1

G0[3,5] = k_bar

G0[4,3] = -(1-theta)
G0[4,4] = 1
G0[4,5] = -theta
G0[4,7] = -1

G0[5,2] = -1
G0[5,4] = 1
G0[5,5] = -1

G0[6,3] = theta
G0[6,5] = -theta
G0[6,6] = 1
G0[6,7] = -1

G0[7,7] = 1

### G1 ####

G1[1,1] = 1

G1[3,1] = -c_bar
G1[3,4] = y_bar
G1[3,5] = (1-delta)*k_bar

G1[7,7] = gamma

Pi[1,1] = 1
Pi[1,2] = -beta*r_bar

Psi[7] = 1

sol2 = gensys(G0,G1,Psi,Pi)

irf1 = irf(sol1,100,0.01)
irf2 = irf(sol2,100,0.01)

plot(irf1[:,1])
plot!(irf2[:,1])

plot(irf1[:,2])
plot!(irf2[:,2])

plot(irf1[:,3])
plot!(irf2[:,3])

plot(irf1[:,4])
plot!(irf2[:,4])

plot(irf1[:,5])
plot!(irf2[:,5])

plot(irf1[:,6])
plot!(irf2[:,6])

G0 = zeros(6,6)
G1 = zeros(6,6)
Pi = zeros(6,2)
Psi = zeros(6)

G0[1,1] = 1
G0[1,2] = -r_bar*beta

G0[3,5] = 1

G0[6,6] = 1

###G1####

G1[1,1] = 1

G1[2,1] = 1
G1[2,3] = 1
G1[2,4] = -1

G1[3,1] = r_bar/theta-delta
G1[3,4] = r_bar/theta
G1[3,5] = (1-delta)

G1[4,3] = 1-theta
G1[4,4] = -1
G1[4,5] = theta
G1[4,6] = 1

G1[5,2] = -1
G1[5,4] = 1
G1[5,5] = -1

G1[6,6] = gamma

Pi[1,1] = 1
Pi[1,2] = -beta*r_bar

Psi[6] = sig

sol3 = gensys(G0,G1,Psi,Pi)
irf3 = irf(sol3,20,0.5)

plot(irf1[:,1])
plot!(irf2[:,1])
plot!(irf3[:,1])

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

G0 = zeros(7,7)
G1 = zeros(7,7)
Pi = zeros(7,2)
Psi = zeros(7)

G0[1,1] = 1
G0[1,2] = -r_bar*beta

G0[2,1] = -1
G0[2,3] = -1
G0[2,4] = 1

G0[3,5] = 1

G0[4,3] = -(1-theta)
G0[4,4] = 1
G0[4,5] = -theta
G0[4,7] = -1

G0[5,2] = 1
G0[5,4] = -1
G0[5,5] = 1

G0[6,3] = 1
G0[6,4] = -1
G0[6,6] = 1

G0[7,7] = 1

##### G1 #####

G1[1,1] = 1

G1[3,1] = -(r_bar/theta-delta)
G1[3,4] = r_bar/theta
G1[3,5] = 1-delta

G1[7,7] = gamma

Pi[1,1] = 1
Pi[1,2] = -beta*r_bar

Psi[7] = 1

sol4 = gensys(G0,G1,Psi,Pi)

irf4 = irf(sol4,100,0.5)

plot(irf1[:,1])
plot!(irf4[:,1])

plot(irf1[:,2])
plot!(irf4[:,2])

plot(irf1[:,3])
plot!(irf4[:,3])

plot(irf1[:,4])
plot!(irf4[:,4])

plot(irf1[:,5])
plot!(irf4[:,5])

plot(irf1[:,6])
plot!(irf4[:,6])

plot(irf1[:,7])
plot!(irf4[:,7])

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


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

function getCallResponseDurations(sim::Simulation)
	@assert(sim.complete)
	@assert(all(call -> call.responseDuration >= 0, sim.calls))
	return map(call -> call.responseDuration, sim.calls)
end

function getAvgCallResponseDuration(sim::Simulation; useMinutes::Bool = false)
	@assert(sim.complete)
	@assert(all(call -> call.responseDuration >= 0, sim.calls))
	return mean(call -> call.responseDuration, sim.calls) * (useMinutes ? 60*24 : 1) # if useMinutes, convert days to minutes
end

function getCallsReachedInTime(sim::Simulation;
	targetResponseDurations::Vector{Float} = sim.targetResponseDurations)
	@assert(sim.complete)
	@assert(all(call -> call.responseDuration >= 0, sim.calls))
	return map(call -> call.responseDuration <= targetResponseDurations[Int(call.priority)], sim.calls)
end

function countCallsReachedInTime(sim::Simulation;
	targetResponseDurations::Vector{Float} = sim.targetResponseDurations)
	@assert(sim.complete)
	@assert(all(call -> call.responseDuration >= 0, sim.calls))
	return count(call -> call.responseDuration <= targetResponseDurations[Int(call.priority)], sim.calls)
end

printDays(t::Float) = string(round(t, digits = 2), " days")
printHours(t::Float) = string(round(t*24, digits = 2), " hours") # convert days to hours, round
printMinutes(t::Float) = string(round(t*24*60, digits = 2), " minutes") # convert days to minutes, round
printSeconds(t::Float) = string(round(t*24*60*60, digits = 2), " seconds") # convert days to seconds, round
printPercent(p::Float) = string(round(p*100, digits = 2), "%")

printKilometres(d::Float) = string(round(d, digits = 2), " km")

function printSimStats(sim::Simulation)
	println("=== Simulation statistics ===")
	println()
	println("Simulation duration = ", printDays(sim.endTime - sim.startTime))
	println()
	printAmbsStats(sim)
	println()
	printCallsStats(sim)
	println()
	printHospitalsStats(sim)
end

function printAmbsStats(sim::Simulation)
	@assert(sim.complete)

	ambsAvgDailyTravelTimes = [amb.statusDurations[ambTravelling] for amb in sim.ambulances] ./ (sim.time - sim.startTime)
	ambsAvgDailyTravelDists = [amb.statusDistances[ambTravelling] for amb in sim.ambulances] ./ (sim.time - sim.startTime)

	# Gadfly.plot(ecdf(ambsAvgDailyTravelTimes*24*60),
		# x="Average daily travel time (minutes)",
		# y="Cdf")

	println("Ambulance statistics:")

	println("Number of ambulances = ", sim.numAmbs)

	println("Average daily travel time: ")
	println(" mean = ", printHours(mean(ambsAvgDailyTravelTimes)))
	println(" std = ", printHours(std(ambsAvgDailyTravelTimes, corrected = false)))
	println(" min = ", printHours(minimum(ambsAvgDailyTravelTimes)))
	println(" max = ", printHours(maximum(ambsAvgDailyTravelTimes)))

	println("Average daily travel distance: ")
	println(" mean = ", printKilometres(mean(ambsAvgDailyTravelDists)))
	println(" std = ", printKilometres(std(ambsAvgDailyTravelDists, corrected = false)))
	println(" min = ", printKilometres(minimum(ambsAvgDailyTravelDists)))
	println(" max = ", printKilometres(maximum(ambsAvgDailyTravelDists)))
end

function printCallsStats(sim::Simulation)
	@assert(sim.complete)

	responseDurations = getCallResponseDurations(sim)
	callsReachedInTime = getCallsReachedInTime(sim)

	println("Call statistics:")

	# # print, then remove, unanswered calls
	# callAnswered = [call.responseDuration != nullTime for call in sim.calls] # may be false if sim not finished, or call cancelled
	# println("Unanswered calls: ", sim.numCalls - sum(callAnswered))
	# responseDurations = responseDurations[callAnswered] # no longer consider unanswered calls

	# Gadfly.plot(ecdf(responseDurations*24*60), x="Call response duration (minutes)", y="Cdf")

	println("Number of calls = ", sim.numCalls)

	println("Response duration: ")
	println(" mean = ", printMinutes(mean(responseDurations)))
	println(" std = ", printMinutes(std(responseDurations, corrected = false)))
	println(" min = ", printMinutes(minimum(responseDurations)))
	println(" max = ", printMinutes(maximum(responseDurations)))

	println("Reached in target response duration: ")
	println(" mean = ", printPercent(mean(callsReachedInTime)))
	println(" sem = ", printPercent(sem(callsReachedInTime)))
end

function printHospitalsStats(sim::Simulation)
	@assert(sim.complete)

	println("Hospital statistics:")
	println("Number of hospitals = ", sim.numHospitals)
	println("Calls:")
	hospitalsNumCalls = map(h -> h.numCalls, sim.hospitals)
	println(" mean = ", round(mean(hospitalsNumCalls), digits = 2))
	println(" std = ", round(std(hospitalsNumCalls, corrected = false), digits = 2))
end

"""
	function calcBatchMeans(values::Vector{Float}, batchSize::Int;
		batchGapSize::Int = 0, rmPartialBatch::Bool = false, returnBatchSizes::Bool = false)
Returns the batch means of `values` batched by size `batchSize`.

# Keyword arguments
- `batchGapSize` is the gap (number of values to ignore) between batches
- `rmPartialBatch` can be set to `true` to remove the last batch if it has size < `batchSize`
- `returnBatchSizes` can be set to `true` to also return the number of values in each batch
"""
function calcBatchMeans(values::Vector{Float}, batchSize::Int;
	batchGapSize::Int = 0, rmPartialBatch::Bool = false, returnBatchSizes::Bool = false)

	@assert(batchSize >= 1)
	@assert(batchGapSize >= 0)

	# calculate number of batches
	numValues = length(values)
	interBatchSize = batchSize + batchGapSize # number of values between start of subsequent batches
	numBatches = max(0, ceil(Int, numValues / interBatchSize))
	if rmPartialBatch
		numBatches = max(0, div(numValues + batchGapSize, interBatchSize))
	end

	# calculate batch means by time
	batchTotals = zeros(Float,numBatches) # total of values in each batch
	batchSizes = zeros(Int,numBatches) # total number of values in each batch
	for i = 1:numValues
		batchIndex = div(i-1, interBatchSize) + 1 # = ceil(Int, i / interBatchSize)
		isInGap = (i > batchIndex * interBatchSize - batchGapSize)
		if batchIndex >= 1 && batchIndex <= numBatches && !isInGap
			batchTotals[batchIndex] += values[i]
			batchSizes[batchIndex] += 1
		end
	end
	batchMeans = batchTotals ./ batchSizes

	# check batch sizes
	if rmPartialBatch
		@assert(all(n -> n == batchSize, batchSizes))
	else
		@assert(all(n -> n == batchSize, batchSizes[1:end-1]))
		@assert(numBatches == 0 || batchSizes[end] == rem(numValues, interBatchSize))
	end

	return returnBatchSizes ? (batchMeans, batchSizes) : batchMeans
end

"""
	function calcBatchMeans(values::Vector{Float}, times::Vector{Float}, batchTime::Float;
		startTime::Float = minimum(times), endTime::Float = maximum(times)*(1+eps(Float)),
		batchGapTime::Float = 0.0, rmPartialBatch::Bool = false, returnBatchSizes::Bool = false)
Returns the batch means of `values` batched by time, where `values[i]` corresponds with `times[i]`.
Each batch will batch values in a time interval `[t, t + batchTime)` for some `t`.
Empty batches have mean `NaN`.

# Keyword arguments
- `startTime` is the time at which batching starts; values with times before this are omitted
- `endTime` is the time at which batching ends; values with times at or after this are omitted
- `batchGapTime` is the time (duration) gap between batches
- `rmPartialBatch` can be set to `true` to remove the last batch if it does not end before `endTime - batchGapTime`.
- `returnBatchSizes` can be set to `true` to also return the number of values in each batch
"""
function calcBatchMeans(values::Vector{Float}, times::Vector{Float}, batchTime::Float;
	startTime::Float = minimum(times), endTime::Float = maximum(times)*(1+eps(Float)),
	batchGapTime::Float = 0.0, rmPartialBatch::Bool = false, returnBatchSizes::Bool = false)

	@assert(length(values) == length(times))
	@assert(batchTime > 0)
	@assert(startTime <= endTime)
	@assert(batchGapTime >= 0)

	# calculate number of batches
	interBatchTime = batchTime + batchGapTime # time difference between start of subsequent batches
	numBatches = max(0, ceil(Int, (endTime - startTime) / interBatchTime))
	if rmPartialBatch
		numBatches = max(0, floor(Int, (endTime - startTime + batchGapTime) / interBatchTime))
	end

	# calculate batch means by time
	batchTotals = zeros(Float,numBatches) # total of values in each batch
	batchSizes = zeros(Int,numBatches) # total number of values in each batch
	for (i,t) in enumerate(times)
		batchIndex = floor(Int, (t - startTime) / interBatchTime) + 1
		isInGap = (t - startTime >= batchIndex * interBatchTime - batchGapTime)
		if batchIndex >= 1 && batchIndex <= numBatches && !isInGap
			batchTotals[batchIndex] += values[i]
			batchSizes[batchIndex] += 1
		end
	end
	batchMeans = batchTotals ./ batchSizes

	return returnBatchSizes ? (batchMeans, batchSizes) : batchMeans
end

# calculate statistics on average response duration based on batch means
function calcBatchMeanResponseDurations(sim::Simulation;
	batchTime::Float = nullTime, warmUpTime::Float = nullTime, coolDownTime::Float = nullTime,
	rmPartialBatch::Bool = true, returnBatchSizes::Bool = true)

	@assert(sim.complete == true)

	# collate data for use by calcBatchMeans function
	n = sim.numCalls
	times = [sim.calls[i].arrivalTime for i = 1:n]
	values = [sim.calls[i].responseDuration for i = 1:n]

	return calcBatchMeans(values, times, batchTime;
		startTime = sim.startTime + warmUpTime, endTime = sim.endTime - coolDownTime,
		rmPartialBatch = rmPartialBatch, returnBatchSizes = returnBatchSizes)
end

# calculate sem for 2d array, calculating along given dim
function semDim(x::Array{T,2}, dim::Int = 1) where T <: Real
	@assert(dim == 1 || dim == 2)
	dim == 1 && return collect([sem(x[:,j]) for j = 1:size(x,2)]') # sem of each col
	dim == 2 && return [sem(x[i,:]) for i = 1:size(x,1)] # sem of each row
end

# For each x[i], plot mean of y[i,:] with two sided confidence interval.
# Assumes that values in y[i,:] are from a normal distribution with unknown standard deviation.
function meanErrorPlot(x::Vector{Float}, y::Array{Float,2}, conf::Float=0.95)
	@assert(length(x) == size(y,1))
	@assert(0 < conf < 1)
	t = StatsFuns.tdistinvcdf(size(y,2)-1, 1-(1-conf)/2) # t-value, for two sided confidence interval
	# StatsFuns.tdistinvcdf(dof, p) # for n samples, dof = n - 1
	@warn("meanErrorPlot still requires use of Plots package; need to change this to RecipesBase")
	return Plots.scatter(x, mean(y, dims = 2), yerror=repeat(t*semDim(y,2), 1, 2),
		xaxis=("x"), yaxis=("y"), m=(:hline), lab="")
end
function meanErrorPlot(x, y, conf::Float=0.95)
	meanErrorPlot(convert(Vector{Float},x), convert(Array{Float,2},y), conf)
end
# plot y values against indices (y[i,:] plotted at x = i)
function meanErrorPlot(y, conf::Float=0.95)
	meanErrorPlot(1:size(y,1), y, conf)
end

# For a vector of values, fit an AR(0) model, and return a p-value for the Durbin-Watson test.
# From wikipedia: "the Durbin–Watson statistic is a test statistic used to detect the presence
# of autocorrelation at lag 1 in the residuals (prediction errors) from a regression analysis".
function calcAR0DurbinWatsonTestPValue(x::Vector{T}) where T <: Real
	xFit = repeat([mean(x)], length(x))
	residuals = x - xFit
	dwTest = HypothesisTests.DurbinWatsonTest(xFit, residuals)
	return HypothesisTests.pvalue(dwTest)
end

# Return half-width of two-sided confidence interval of estimate of mean.
# Samples should come from a normal distribution, standard deviation is assumed to be unknown.
function tDistrHalfWidth(x::Vector{T}; conf = 0.95) where T <: Real
	t = StatsFuns.tdistinvcdf(length(x)-1, 1-(1-conf)/2) # t-value, for two sided confidence interval
	return t * sem(x)
end

# confidence interval
function confInterval(mhw::MeanAndHalfWidth)
	return (mhw.mean - mhw.halfWidth, mhw.mean + mhw.halfWidth)
end

# Return a list of the main stats periods.
# Assumes that the first period (after any warm-up) is the main period.
# Can only return multiple periods if `sim.stats.periodDurationsIter` is cyclical.
# Useful for batches in a single simulation replication.
function getPeriodStatsList(sim::Simulation)::Vector{SimPeriodStats}
	@assert(sim.complete)
	stats = sim.stats # shorthand
	i = stats.warmUpDuration == 0 ? 1 : 2 # ignore first period if it is for warm-up
	itr = stats.periodDurationsIter.itr # can be empty if stats control not set
	k = isa(itr, Iterators.Cycle) ? length(itr.xs) : length(itr) # gap between period indices
	periods = stats.periods[i:k:end]
	if isempty(periods) return periods end
	periodDuration = periods[1].duration
	if !isapprox(periods[end].duration, periodDuration) pop!(periods) end
	@assert(all(p -> isapprox(p.duration, periodDuration), periods))
	return periods
end

# given the sim replications, return a list with the main stats period for each rep
function getRepsPeriodStatsList(reps::Vector{Simulation}; periodIndex::Int = nullIndex)::Vector{SimPeriodStats}
	@assert(all(rep -> rep.complete, reps))
	if periodIndex == nullIndex
		warmUpDuration = reps[1].stats.warmUpDuration
		@assert(all(rep -> rep.stats.warmUpDuration == warmUpDuration, reps))
		periodIndex = warmUpDuration == 0 ? 1 : 2
	end
	periods = [rep.stats.periods[periodIndex] for rep in reps]
	@assert(all(p -> isapprox(p.duration, periods[1].duration), periods))
	return periods
end

# given a sim replication, return the main stats period
function getRepPeriodStats(rep::Simulation; periodIndex::Int = nullIndex)::SimPeriodStats
	return getRepsPeriodStatsList([rep], periodIndex = periodIndex)[1]
end

# Return a dictionary of statistics from a list of period statistics.
# Periods should be the same duration.
# Confidence intervals assume that samples obtained from periods are IID, and are from a population with a normal distribution and unknown standard deviation.
# The function flatten(dict) may be useful if writing this dict to file.
function statsDictFromPeriodStatsList(periods::Vector{SimPeriodStats}; conf = 0.95)

	if isempty(periods)
		period = SimPeriodStats()
		period.duration = 0.0
		ambStatuses = (instances(AmbStatus)..., instances(AmbStatusSet)...)
		period.ambulance.statusDurations = Dict([s => 0.0 for s in ambStatuses])
		period.ambulance.statusDistances = Dict([s => 0.0 for s in ambStatuses])
		period.callPriorities = Dict([p => CallStats() for p in priorities])
		periods = [period]
	end

	duration = periods[1].duration
	ambDays = length(periods[1].ambulances) * duration
	@assert(all(p -> isapprox(p.duration, duration), periods))

	d = statsDict = Dict{String,Any}()

	meanAndHalfWidth(x; conf = conf) = MeanAndHalfWidth(mean(x), tDistrHalfWidth(x; conf = conf))

	###########
	# ambulance

	d["ambs"] = Dict{String,Any}()

	function getAmbsStat(statName::Symbol)
		x = [getfield(p.ambulance, statName) for p in periods] / ambDays
		return meanAndHalfWidth(x)
	end
	function getAmbsDurationStat(s::Union{AmbStatus,AmbStatusSet})
		x = [p.ambulance.statusDurations[s] for p in periods] / ambDays
		return meanAndHalfWidth(x)
	end
	function getAmbsDistanceStat(s::Union{AmbStatus,AmbStatusSet})
		x = [p.ambulance.statusDistances[s] for p in periods] / ambDays
		return meanAndHalfWidth(x)
	end

	# counts
	d["ambs"]["avgDailyNumCallsTreated"] = getAmbsStat(:numCallsTreated)
	d["ambs"]["avgDailyNumCallsTransported"] = getAmbsStat(:numCallsTransported)
	d["ambs"]["avgDailyNumDispatches"] = getAmbsStat(:numDispatches)
	d["ambs"]["avgDailyNumDispatchesFromStation"] = getAmbsStat(:numDispatchesFromStation)
	d["ambs"]["avgDailyNumDispatchesOnRoad"] = getAmbsStat(:numDispatchesOnRoad)
	d["ambs"]["avgDailyNumDispatchesOnFree"] = getAmbsStat(:numDispatchesOnFree)
	d["ambs"]["avgDailyNumRedispatches"] = getAmbsStat(:numRedispatches)
	d["ambs"]["avgDailyNumMoveUps"] = getAmbsStat(:numMoveUps)
	d["ambs"]["avgDailyNumMoveUpsFromStation"] = getAmbsStat(:numMoveUpsFromStation)
	d["ambs"]["avgDailyNumMoveUpsOnRoad"] = getAmbsStat(:numMoveUpsOnRoad)
	d["ambs"]["avgDailyNumMoveUpsOnFree"] = getAmbsStat(:numMoveUpsOnFree)
	d["ambs"]["avgDailyNumMoveUpsReturnToPrevStation"] = getAmbsStat(:numMoveUpsReturnToPrevStation)

	# durations - status
	d["ambs"]["avgDailySleepingDurationHours"] = getAmbsDurationStat(ambSleeping) * 24
	d["ambs"]["avgDailyIdleAtStationDurationHours"] = getAmbsDurationStat(ambIdleAtStation) * 24
	d["ambs"]["avgDailyMobilisingDurationHours"] = getAmbsDurationStat(ambMobilising) * 24
	d["ambs"]["avgDailyGoingToCallDurationHours"] = getAmbsDurationStat(ambGoingToCall) * 24
	d["ambs"]["avgDailyAtCallDurationHours"] = getAmbsDurationStat(ambAtCall) * 24
	d["ambs"]["avgDailyGoingToHospitalDurationHours"] = getAmbsDurationStat(ambGoingToHospital) * 24
	d["ambs"]["avgDailyAtHospitalDurationHours"] = getAmbsDurationStat(ambAtHospital) * 24
	d["ambs"]["avgDailyReturningToStationDurationHours"] = getAmbsDurationStat(ambReturningToStation) * 24
	d["ambs"]["avgDailyMovingUpToStationDurationHours"] = getAmbsDurationStat(ambMovingUpToStation) * 24

	# durations - status sets
	d["ambs"]["avgDailyWorkingDurationHours"] = getAmbsDurationStat(ambWorking) * 24
	d["ambs"]["avgDailyBusyDurationHours"] = getAmbsDurationStat(ambBusy) * 24
	d["ambs"]["avgDailyFreeDurationHours"] = getAmbsDurationStat(ambFree) * 24
	d["ambs"]["avgDailyTravelDurationHours"] = getAmbsDurationStat(ambTravelling) * 24
	d["ambs"]["avgDailyGoingToStationDurationHours"] = getAmbsDurationStat(ambGoingToStation) * 24

	# distances
	d["ambs"]["avgDailyTravelDistanceKms"] = getAmbsDistanceStat(ambTravelling)
	d["ambs"]["avgDailyGoingToCallDistanceKms"] = getAmbsDistanceStat(ambGoingToCall)
	d["ambs"]["avgDailyGoingToHospitalDistanceKms"] = getAmbsDistanceStat(ambGoingToHospital)
	d["ambs"]["avgDailyReturningToStationDistanceKms"] = getAmbsDistanceStat(ambReturningToStation)
	d["ambs"]["avgDailyMovingUpToStationDistanceKms"] = getAmbsDistanceStat(ambMovingUpToStation)

	###########
	# call

	d["calls"] = Dict{String,Any}()

	for priority in instances(Priority) # will use nullPriority to indicate all priorities
		name = priority == nullPriority ? "all" : string(priority)
		d1 = d["calls"][name] = Dict{String,Any}()

		callStatsList = priority == nullPriority ? [p.call for p in periods] : [p.callPriorities[priority] for p in periods]
		function getCallsStat(dividendName::Symbol, divisorName::Symbol = :numCalls)
			x = [getfield(callStats, dividendName) for callStats in callStatsList]
			y = [getfield(callStats, divisorName) for callStats in callStatsList]
			i = findall(x -> x != 0, y)
			# if isempty(i) return MeanAndHalfWidth(0,0) end # do not do this; lack of samples doesn't imply mean of zero
			return meanAndHalfWidth(x[i]./y[i])
		end

		# fractions (between 0 and 1); don't need to write "avgFrac", as fraction is already an average
		d1["fracQueued"] = getCallsStat(:numQueued)
		d1["fracBumped"] = getCallsStat(:numBumped)
		d1["fracTransported"] = getCallsStat(:numTransports)
		d1["fracResponsesInTime"] = getCallsStat(:numResponsesInTime)

		# durations
		d1["avgDispatchDelayMinutes"] = getCallsStat(:totalDispatchDelay) * (24*60)
		d1["avgOnSceneDurationMinutes"] = getCallsStat(:totalOnSceneDuration) * (24*60)
		d1["avgHandoverDurationMinutes"] = getCallsStat(:totalHandoverDuration, :numTransports) * (24*60)
		d1["avgQueuedDurationMinutes"] = getCallsStat(:totalQueuedDuration, :numQueued) * (24*60)
		d1["avgBumpedDurationMinutes"] = getCallsStat(:totalBumpedDuration, :numBumped) * (24*60)
		d1["avgWaitingForAmbDurationMinutes"] = getCallsStat(:totalWaitingForAmbDuration) * (24*60)
		d1["avgResponseDurationMinutes"] = getCallsStat(:totalResponseDuration) * (24*60)
		d1["avgAmbGoingToCallDurationMinutes"] = getCallsStat(:totalAmbGoingToCallDuration) * (24*60)
		d1["avgTransportDurationMinutes"] = getCallsStat(:totalTransportDuration, :numTransports) * (24*60)
		d1["avgServiceDurationMinutes"] = getCallsStat(:totalServiceDuration) * (24*60)

		# misc
		d1["avgNumBumps"] = getCallsStat(:numBumps)
		d1["avgDailyNumCalls"] = meanAndHalfWidth([callStats.numCalls / duration for callStats in callStatsList])
	end

	###########

	return statsDict
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    TopologicalRelation

A topological relation between faces of a [`Mesh`](@ref) implemented
for a given [`Topology`](@ref).

An object implementing this trait is a functor that can be evaluated
at an integer index representing the face.

## Examples

```julia
# create boundary relation mapping
# 2-faces to 0-faces (i.e. vertices)
∂ = Boundary{2,0}(topology)

# list of vertices for first face
∂(1)
```

## References

* Floriani, L. & Hui, A. 2007. [Shape representations
  based on simplicial and cell complexes]
  (https://diglib.eg.org/handle/10.2312/egst.20071055.063-087)
"""
abstract type TopologicalRelation end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("toporelations/boundary.jl")
include("toporelations/coboundary.jl")
include("toporelations/adjacency.jl")


text:
#Demo 7: Eye diagrams
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))

#==Input data
===============================================================================#
tbit = 1e-9 #Bit period
osr = 20 #samples per bit
nbit_Π = 5 #Π-pulse length, in number of bits
nsamples = 127

#==Computations
===============================================================================#
seq = 1.0*prbs(reglen=7, seed=1, nsamples=nsamples)
tΠ = DataF1(0:(tbit/osr):(nbit_Π*tbit))

#Generate parameter sweeps:
sweeplist = PSweep[
	PSweep("EXTRALVL", [1]) #Add extra level to test more code
	PSweep("tau", tbit.*[1/5, 1/2.5, 1/1.5])
]

#Generate data:
Π = fill(DataHR, sweeplist) do EL, tau
	return pulse(tΠ, Pole(1/tau,:rad), tpw=tbit)
end
pat = (pattern(seq, Π, tbit=tbit)-0.5)*2 #Center data pattern
patRS = DataRS(pat)

#==Generate plot
===============================================================================#
#TODO: supply eyeparam() as a FoldedAxis or Plot constructor???
eyeparam(tbit; teye=1.5*tbit, tstart=0) = cons(:a,
	xfolded = set(tbit, xstart=tstart, xmax=teye),
	xaxis = set(min=0, max=teye), #Force limits on exact data range.
)
eyeaxis = eyeparam(tbit, teye=1.5*tbit, tstart=-.15*tbit)

p1 = push!(cons(:plot, vvst, eyeaxis, title="Eye (DataHR)", legend=false),
	cons(:wfrm, pat, label="eye"),
)
p2 = push!(cons(:plot, vvst, title="Pattern", legend=false),
	cons(:wfrm, pat, label="pat"),
)
p3 = push!(cons(:plot, vvst, eyeaxis, title="Eye (DataRS)", legend=false),
	cons(:wfrm, patRS, label="eye"),
)

pcoll = push!(cons(:plotcoll, title="Eye Diagram Tests"), p1, p2, p3)
	pcoll.ncolumns = 2

#==Display results in pcoll
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll


text:
using Pkg.Artifacts

using CLIMA.ArtifactWrappers

# Get dycoms dataset folder:
dycoms_dataset = ArtifactWrapper(
    joinpath(@__DIR__, "Artifacts.toml"),
    "dycoms",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/bxau6i46y6ikxn2sy9krgz0sw5vuptfo.nc",
        filename = "test_data_PhaseEquil.nc",
    ),],
)
dycoms_dataset_path = get_data_folder(dycoms_dataset)

@testset "Data tests" begin
    FT = Float64
    e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice = MT.tested_convergence_range(50, FT)
    data = joinpath(dycoms_dataset_path, "test_data_PhaseEquil.nc")
    ds_PhaseEquil = Dataset(data, "r")
    e_int = Array{FT}(ds_PhaseEquil["e_int"][:])
    ρ = Array{FT}(ds_PhaseEquil["ρ"][:])
    q_tot = Array{FT}(ds_PhaseEquil["q_tot"][:])

    # ts = PhaseEquil.(e_int, ρ, q_tot) # Fails, but we should get these to pass!
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using Distributions
using Counters, Distributed

"""
`binom_rv(n,p)` generates a single random value according to
the binomial distribution `B(n,p)`.
"""
binom_rv(n::Int, p::Real = 0.5) = rand(Binomial(n, p))

"""
`binomial_counts(n,p,reps)` generates `reps` binomial random values `B(n,p)`
returning a `Counter` that reports how many times each value was observed.
"""
function binomial_counts(n::Int, p::Real, reps::Int)
    c = Counter{Int}()
    for k = 1:reps
        x = binom_rv(n, p)
        c[x] += 1
    end
    return c
end

"""
`parallel_binomial_counts(n,p,reps,rounds)` generates `reps*rounds`
binomial random values `B(n,p)` returning a `Counter` that reports how
many times each value was observed. In this regard, this is the same
as `binomial_counts`. However, this is done by multiple processors by
making `rounds` calls to `binomial_counts(n,p,reps)` and combining
the results.
"""
function parallel_binomial_counts(n::Int, p::Real, reps::Int, rounds::Int)
    counts = @distributed (+) for k = 1:rounds
        binomial_counts(n, p, reps)
    end
    return counts
end


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4
coefs = [0.75;0.5]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1,1, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def, m_step=false)

dsize = parse(Int, ARGS[1])
t_steps = parse(Int, ARGS[2])
t_group = parse(Int, ARGS[3])
ngpus = parse(Int, ARGS[4])
nx = dsize
ny = dsize
nz = dsize
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz);

g1 = ApplyStencil(st_inst, data, t_steps);
g1 = ApplyStencil(st_inst, data, t_steps);
g1 = ApplyStencil(st_inst, data, t_steps);
CUDA.device_reset!.(collect(devices()))
g2 = ApplyMultiGPU(ngpus , st_inst, t_steps, data, t_group=t_group);
CUDA.device_reset!.(collect(devices()))
g2 = ApplyMultiGPU(ngpus , st_inst, t_steps, data, t_group=t_group);
CUDA.device_reset!.(collect(devices()))
g2 = ApplyMultiGPU(ngpus , st_inst, t_steps, data, t_group=t_group);
CUDA.device_reset!.(collect(devices()))
g2 = ApplyMultiGPU(ngpus , st_inst, t_steps, data, t_group=t_group);
CUDA.device_reset!.(collect(devices()))
#
# println("AFTER WARMUP")
#
# CUDA.cuProfilerStart()
# NVTX.@range "single" begin
# g1 = ApplyStencil(st_inst, data, t_steps)end
# NVTX.@range "multi" begin
# g2 = ApplyMultiGPU(ngpus , st_inst, t_steps, data, t_group=t_group) end


text:
function ACQR_simula(nt::Int, seed::Int)
	#
	# simulated model to test the functions
	#
	# nt: number of time instants
	# seed: seed for the random numbers

	# random seed
	# ------------------------------
	Random.seed!(seed)

	# simulated system
	# ----------------------------
	m = simula01()

	# vibration parameters
	# ----------------------------
	# time interval (s)
	dt = m["dt"]

	# state-space matrices
	# --------------------------------------
	A = m["A"]
	B = m["B"]
	C = m["C"]
	D = m["D"]

	# input definition
	# ---------------------------------------
	u = 100*randn(2,nt) # if 100 is not used, Q and R are too small

	# state equation, no noise
	# ---------------------------------
	x = zeros(4,nt)
	w = B*u
	for t=1:nt-1
		x[:,t+1] = A*x[:,t] + w[:,t]
	end
	Q = cov(w')

	# observation equation
	# -------------------------------------
	v = D*u
	y = C*x + v

	# covariance matrices
	Q = cov(w')
	R = cov(v')
	S = cov(w',v')

	# starting values
	m1 = C\y[:,1:1]
	P1 = 0*A

	m = Dict("dt"=>dt,
	"y"=>y,"u"=>u,
	"A"=>A,"B"=>B,"C"=>C,"D"=>D,
	"Q"=>Q,"R"=>R,"S"=>S,
	"m1"=>m1,"P1"=>P1)

	return m
end


text:
# import packages
using PyPlot
using AnsiColor

## add path to project files
dirname = "/Users/dajohnson/Desktop/Julia/Projects/UTest";
push!(LOAD_PATH, dirname);

## include clear
include("clear.jl");

## Clear and Close Figures
clear();
cla();                               # clear current axis
clf();                               # clear current plot
close();                             # close figure window

## include external functions
include("joinString.jl");
#include("gradientDescentMulti.jl");
#include("computeCostMulti.jl");

# prints out array values
showprintln(x) = (show(x); println("\n"));
showprint(x) = (show(x); print());
showinfoln(u,x) = (begin @sprintf("%s ", u); show(x); println("\n") end);

## functions
function Add(num1::Int64, num2::Int64)
	return num1 + num2;
end

# Adds 2 or more numbers (integers)
function addNum(args...)
	begin
		total = 0
		num = 0
		numArgs = size(args,1)
		#@show(numArgs)
		for i = 1:numArgs
			num = args[i]
			total += num
		end
		return total
	end
end

# concatenating two strings together
function stringTogether(args...)
	begin
		newStr = ""
		str = []
		numArgs = size(args,1)
        @show(numArgs)
		for i = 1:numArgs
			push!(str, args[i])
			newStr += str * " "
		end
		return newStr
	end
end

# test suite Data Definitions and Functions
results = Dict();
results["total"] = 0;
#results["good"] = 0;
results["bad"] = 0;

## report
# https://libraries.io/julia/AnsiColor
function report()
   begin
        total = results["total"];
        #good = results["good"];
        bad =  results["bad"];
        passed = total - bad;
        if total >= 1 && passed >= 1
            println(colorize(:light_white,"Of $total tests, $bad failed, $passed passed.", background="light_black", mode="blink"))
        else
            println(colorize(:red,"Of $total test, $bad failed, $passed passed.", background="light_black", mode="blink"))
        end
   end
end

# an anonymous block of code
#=
function codeblock(code)
    return code
end
=#

function checkExpect(f::Function, params::Any, expected::Any, name::AbstractString ="")
    begin
        # run unit test
        result = f(params...);

        if strwidth(chomp(name)) < 0 || chomp(name) != ""
			println("for unit test: " * "$name")
		else
			println("for unit test: " * "test description unknown")
		end

        if result != expected
			println("test: " * colorize(:red, "failed", background="black", mode="bold"))
			println("Expected $expected" * ", but was $result")
            println("info: " * colorize(:light_yellow, "Algorithm design is not well formed!", background="black", mode="bold"))
            println("")
            results["bad"] += 1;
            results["total"] += 1;
		else
			println("test: " * colorize(:green,"passed", background="black", mode="bold"))
			println("Function: $f consumed the following argument(s) $params, and produced $result, which matches $expected the expected result.")
            println("info: " * colorize(:light_blue, "Algorithm is well formed!", background="black", mode="default"))
            println("")
            #results["good"] += 1;
            results["total"] += 1;
		end
        return result
    end
end

# Data Type Definitions

## examples:
# Add 1, 2 should produce 3
# Add 23, 2 should produce 25
param1 = [1,2];
param2 = [23,2];
param4 = [1,2,3];

# stringTogether Darel, Johnson should produce Darel Johnson
param3 = ["Darel","Johnson"];

# use a partial codeblock for the Add function template
# parameters
num1 = 23
num2 = 2
# template
func_for_add = begin
                    ans = 0;
                    ans = num1 + num2;
                end

# use a partial codeblock for the stringTogether function template
# parameters
firstName = "Darel"
lastName = "Johnson"
# template
func_for_stringtogether = begin
                            result = ""
                            if strwidth(firstName) > 1 && strwidth(lastName) > 1
                                result = strip(firstName) * " " * strip(lastName)
                            end
                          end

checkExpect(Add, param1,3);
checkExpect(Add, param1, 1 + 3, "Add two numbers");
checkExpect(Add, param2, 23 + 2, "Add two numbers");
checkExpect(joinString, param3, "Darel Johnson", "Concatenating two strings");
#checkExpect(stringTogether, param3, "Darel Johnson", "Concatenating two strings");
checkExpect(addNum, param4, 6, "Add three numbers");
checkExpect(Add, param2, func_for_add, "Add two numbers with algorithm in func_def");
checkExpect(Add, param2, func_for_add, "testing codeblock")
checkExpect(joinString, param3, func_for_stringtogether, "Concatenating two strings with algorithm in func_def");

# report test results
report();


text:
abstract AbstractPoint
abstract AbstractPoint2D <: AbstractPoint
import Base.convert

immutable Point2D{T <: Real} <: AbstractPoint2D
    x::T
    y::T
end

Base.isless(p1::Point2D, p2::Point2D) = isless(p1.x, p2.x) || (p1.x == p2.x && isless(p1.y, p2.y))

getpoint(points::Array{Point2D,1},i::Int) = points[i]
getx(p::Point2D) = p.x
gety(p::Point2D) = p.y

getpoint{T <: Real}(points::Array{T, 2}, i::Int) = points[i,:]
getx{T <: Real}(p::Array{T,1}) = p[1]
gety{T <: Real}(p::Array{T,1}) = p[2]

# convert{T <: Real}(Point2D, p::Array{T,1}) = Point2D(getx(p), gety(p))

max_stars_repo_name: helgith/GeometricStructures

max_issues_repo_name: helgith/GeometricStructures

max_forks_repo_name: helgith/GeometricStructures


text:
using Documenter, MPSKit

makedocs(modules=[MPSKit],
            sitename="MPSKit.jl",
            format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true",
                                        mathengine = MathJax()),
            pages = [
                "Home" => "index.md",
                "Manual" => ["man/intro.md", "man/states.md", "man/operators.md",
                                "man/algorithms.md","man/environments.md"],
                "Tutorials" => ["tut/xxz_groundstate.md","tut/timeev.md","tut/haldane.md"],
                "Library" => ["lib/lib.md"],
            ])

deploydocs(
    repo = "github.com/maartenvd/MPSKit.jl.git",
)


text:
# We create another version of the Runtime when are evaluating a derived/input function,
# which is able to keep track of the dependencies invoked from that derived function, so
# they can be recorded together with the function's results.
#
# _TracingRuntime is an immutable struct, because it is branched for each derived function
# call with a fresh dependency trace, so that derived functions can be safely called in
# parallel on multiple threads. This makes the Runtime thread-safe, as long as it's used
# correctly.
struct _TracingRuntime{CT,ST<:AbstractSalsaStorage} <: Runtime{CT,ST}
    # NOTE: We store a Ptr{} to the parent Runtime rather than the Storage instance itself,
    # in order to keep this struct as `isbits`. This is super critical because we create a
    # new instance of this struct on _every function call_ in Salsa!  We also use this
    # pointer to allow derived functions to access their Runtime's context.
    #
    # This is guaranteed to be safe, since the parent Runtime will always exist until the
    # function call that created this tracing runtime has completed.
    tl_runtime::Ptr{_TopLevelRuntime{CT,ST}}

    # A trace structure to store the dependencies of derived functions as they are
    # encountered. This field is mutable, but we create a new, empty trace every time we
    # branch the Runtime. It's locked internally to allow spawned derived functions on
    # separate threads to record dependencies to the parent thread.
    #
    # NOTE: We store an index into `tls_trace_pool()`, rather than a `TraceOfDependencyKeys`
    # object itself here, in order to keep this struct as `isbits`. This is super critical
    # because we create a new instance of this struct on _every function call_ in Salsa!
    immediate_dependencies_id::TraceId

    function @__MODULE__().new_trace_runtime!(
        old_rt::_TopLevelRuntime{CT,ST},
        key::DependencyKey,
    )::_TracingRuntime{CT,ST} where {CT,ST<:AbstractSalsaStorage}
        new{CT,ST}(
            reinterpret(Ptr{ST}, pointer_from_objref(old_rt)),
            # Start a new, empty trace (with the provided call stack if in debug mode)
            if Salsa.Debug.debug_enabled()
                get_trace_with_call_stack(SalsaStackFrame(key, nothing))
            else
                get_trace_with_call_stack(nothing)
            end,
        )
    end

    function @__MODULE__().new_trace_runtime!(
        old_rt::_TracingRuntime{CT,ST},
        key::DependencyKey,
    )::_TracingRuntime{CT,ST} where {CT,ST<:AbstractSalsaStorage}
        # Push the new computation onto the current Runtime (if it's not there already)
        push_key!(old_rt, key)
        # Create a new linked list node (pointing to the old stack trace if debug mode).
        new_trace = if Salsa.Debug.debug_enabled()
            get_trace_with_call_stack(SalsaStackFrame(key, trace(old_rt).call_stack))
        else
            get_trace_with_call_stack(nothing)
        end
        new{CT,ST}(old_rt.tl_runtime, new_trace)
    end
end

function push_key!(rt::_TracingRuntime, depkey)
    tr = trace(rt)

    # Test for cycles if in debug mode
    @debug_mode if stack_has_key(tr.call_stack, depkey)
        throw(DependencyCycleException(depkey))
    end

    push_key!(tr, depkey)
    return nothing
end

########## Implementation of Runtime API

context(rt::_TracingRuntime) = unsafe_load(rt.tl_runtime).context

storage(rt::_TracingRuntime) = unsafe_load(rt.tl_runtime).storage

trace(rt::_TracingRuntime) = get_trace(rt.immediate_dependencies_id)

collect_call_stack(rt::_TracingRuntime) = _collect_call_stack_frames(trace(rt).call_stack)
_collect_call_stack_frames(::Nothing) = DependencyKey[]
_collect_call_stack_frames(frame::SalsaStackFrame) = collect(frame)

collect_trace(rt::_TracingRuntime) = collect_trace(trace(rt))

function destruct_trace!(rt::_TracingRuntime)
    release_trace_id(rt.immediate_dependencies_id)
    return nothing
end

# Safety overload.
function unmemoized_input_lookup(rt::_TracingRuntime, key)
    error("Attempted an unmemoized lookup of $key inside of a derived function.")
end

function previous_output(rt::_TracingRuntime)
    dependency_key = trace(rt).call_stack.dp
    return _unwrap_salsa_value(rt, _previous_output_internal(rt, dependency_key))
end


text:
using KhepriAutoCAD
using Test

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


text:
module Aerodynamics
using DataFrames
using CSV
using LinearAlgebra
using StaticArrays
using Interpolations
using SymEngine
using ..RocketlandDefns
import StaticArrays

function load_aerodata(liftdrag::String, finforce::Union{String, Nothing}=nothing)
	lddf = CSV.read(liftdrag) #numbers in newtons
	#normalize
	lddf[:lift_p] = lddf[:lift]
	lddf[:drag_p] = lddf[:drag]
	#mach, aoa
	mach=0.0:0.025:1.5
	aoa=cosd(180):1/90:cosd(0)
	lift_itrp = extrapolate(scale(interpolate(reshape(lddf[:lift_p],length(aoa), length(mach)), BSpline(Cubic(Line(OnGrid())))), aoa, mach), Flat())
	drag_itrp = extrapolate(scale(interpolate(reshape(lddf[:drag_p],length(aoa), length(mach)), BSpline(Cubic(Line(OnGrid())))), aoa, mach), Flat())
	trq_itrp = extrapolate(scale(interpolate(reshape(lddf[:torque],length(aoa), length(mach)), BSpline(Cubic(Line(OnGrid())))), aoa, mach), Flat())

	if !isnothing(finforce)
		ffdf = CSV.read(finforce)

	end
	AtmosphericData(drag_itrp, lift_itrp, trq_itrp, 1.0, 1.0)
end

function rescale_aerodata(data::AtmosphericData, Ul::Float64, Ut::Float64, Um::Float64)
	return AtmosphericData(data.drag_itrp, data.lift_itrp, data.trq_itrp, 1/(Ul * Um/Ut^2), 1/Ul)
end

function rescale_aerodata(data::ExoatmosphericData, Ul::Float64, Ut::Float64, Um::Float64)
	return data
end

function aero_force(data::AtmosphericData, bv::AbstractArray{T}, vel::AbstractArray{T}, spds::Float64)::Tuple{Vector{T}, Vector{T}} where T<:Number
	dp = dot(bv,vel)/norm(vel)
	cos_aoa = clamp(dp/norm(bv), -1.0,1.0)
	mach = norm(vel)/spds
	if abs(dp) >= 0.95
		drag = data.drag_itrp(cos_aoa,mach)::T*data.force_scalar
		dragf = drag*vel/norm(vel)
		return dragf, zeros(T, 3)
	else
		drag = data.drag_itrp(cos_aoa,mach)::T*data.force_scalar
		lift = data.lift_itrp(cos_aoa,mach)::T*data.force_scalar
		trq = data.trq_itrp(cos_aoa,mach)::T*data.length_scalar*data.force_scalar
		trqd = cross(vel, bv)
		liftd = cross(-trqd,vel)
		liftd = liftd/norm(liftd)
		trqd = trqd/norm(trqd)
		dragf = drag*vel/norm(vel)
		liftf = lift*liftd
		return dragf + liftf, trqd*trq
	end
end

function aero_force(data::AtmosphericData, bv::AbstractArray{T}, vel::AbstractArray{T}, spds::Float64)::Tuple{Vector{T}, Vector{T}} where T<:Basic
	dp = sum(bv .* vel)
	cos_aoa = dp #/sqrt(sum(bv .^ 2))
	mach = sqrt(sum(vel .^ 2))/spds

	drag = SymFunction("drag")(cos_aoa,mach,symbols("pinfo"))*data.force_scalar
	lift = SymFunction("lift")(cos_aoa,mach,symbols("pinfo"))*data.force_scalar
	trq = SymFunction("trq")(cos_aoa,mach,symbols("pinfo"))*data.length_scalar*data.force_scalar
	trqd = cross(vel, bv)
	liftd = cross(-trqd,vel)
	liftd_norm = sqrt(sum(liftd .^ 2))#sqrt(sum(vel .^ 2) * sum(bv .^ 2) * (sum(vel .^ 2) - cos_aoa^2))
	drag_norm = sqrt(sum(vel .^ 2))
	liftd = (liftd)/liftd_norm
	liftf = lift*(SymFunction("ifnz").(liftd_norm, liftd))
	dragf = SymFunction("ifnz").(drag_norm, drag*vel/drag_norm)
	total_force = dragf + liftf
	return SymFunction("ifnz").(total_force, total_force), trq*(SymFunction("ifnz").(trqd/sqrt(sum(trqd .^ 2)), trqd))
end

function aero_force(data::ExoatmosphericData, bv::AbstractArray{T} where T, vel::AbstractArray{T} where T, spds::Float64)
	return [0.0,0.0,0.0]
end

function direct_drag(data::AtmosphericData, cos_aoa, mach)
	return data.drag_itrp(cos_aoa, mach)
end

function direct_drag(::Type{Val{:jac}}, data::AtmosphericData, cos_aoa, mach)
	return Interpolations.gradient(data.drag_itrp, cos_aoa, mach)::StaticArrays.SArray{Tuple{2},Float64,1,2}
end

function direct_lift(data::AtmosphericData, cos_aoa, mach)
	return data.lift_itrp(cos_aoa, mach)
end

function direct_lift(::Type{Val{:jac}}, data::AtmosphericData, cos_aoa, mach)
	return Interpolations.gradient(data.lift_itrp, cos_aoa, mach)::StaticArrays.SArray{Tuple{2},Float64,1,2}
end

function direct_trq(data::AtmosphericData, cos_aoa, mach)
	return data.trq_itrp(cos_aoa, mach)
end

function direct_trq(::Type{Val{:jac}}, data::AtmosphericData, cos_aoa, mach)
	return Interpolations.gradient(data.trq_itrp, cos_aoa, mach)::StaticArrays.SArray{Tuple{2},Float64,1,2}
end

function direct_drag(data::ExoatmosphericData, cos_aoa, mach)
	return 0.0
end

function direct_lift(data::ExoatmosphericData, cos_aoa, mach)
	return 0.0
end

function direct_trq(data::ExoatmosphericData, cos_aoa, mach)
	return 0.0
end

function direct_drag(::Type{Val{:jac}}, data::ExoatmosphericData, cos_aoa, mach)
	return StaticArrays.SArray{Tuple{2}}(0.0,0.0)
end

function direct_lift(::Type{Val{:jac}}, data::ExoatmosphericData, cos_aoa, mach)
	return StaticArrays.SArray{Tuple{2}}(0.0,0.0)
end

function direct_trq(::Type{Val{:jac}}, data::ExoatmosphericData, cos_aoa, mach)
	return StaticArrays.SArray{Tuple{2}}(0.0,0.0)
end

export direct_trq, direct_lift, direct_drag
end

max_stars_repo_name: BenChung/SuccessiveConvexification

max_issues_repo_name: BenChung/SuccessiveConvexification

max_forks_repo_name: BenChung/SuccessiveConvexification


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

using Markdown
using InteractiveUtils

# ╔═╡ 80b90b23-c999-4836-ae83-d1b7017daa46
begin
	using Pkg
	Pkg.activate(joinpath(Base.current_project(), "..", ".."))  # Tells Pluto.jl we are going to use local enviroment
	using InteractionNetworkModels, Distributions, Plots
end

# ╔═╡ 0439d40e-768e-44d0-85bc-417b063f4315
md"
# Overview

Here we illustrate how one can sample from the Spherical Interaction Sequence (SIS) models. This amounts to the following steps
1. Define an `model::SIS` model object, specifying the required parameters
2. Construct an `mcmc::SisMcmcSampler` object
3. Call `mcmc(model)` to generate samples
4. Explore outputs manually or via pre-defined visualisations

"

# ╔═╡ 8dda732d-6cd5-4ac9-a749-277b1ba49cea
md"
...constructing model...
"

# ╔═╡ 09b424d8-b54d-4d12-9367-c0413c46475e
model = SIS(
	[[1,2],[1,2],[2,3,4],[2,3,4]],
	5.0,
	FastEditDistance(FastLCS(100),100),
	1:10,
	30, 30
);

# ╔═╡ 48ad33b4-b44c-482e-aa98-baa5a3e890ce
md"
...constructing MCMC sampler...
"

# ╔═╡ 5f990069-ba98-4930-9806-3ed8c25e7253
path_dist = PathPseudoUniform(model.V,TrGeometric(0.8, 1, 10));

# ╔═╡ 699e9669-b105-4d5a-962b-7bf63acb8b49
mcmc = SisMcmcInsertDelete(path_dist);

# ╔═╡ 1504fd48-9e5e-4536-a0ed-7ebc26126b92
md"
...now to get a sample we can use `mcmc` like a function, that is, we can call `mcmc(model)` to obtain an MCMC sample from `model`...
"

# ╔═╡ 508ba9f3-fc56-4012-a447-8bb6c675f271
x = mcmc(model, desired_samples=5000, lag=5, burn_in=1000)

# ╔═╡ ef5fa55b-198c-41a6-acb4-e0bb964c822a
md"...we can now visualise a trace plot by simply calling `plot(x)` where `x::SisMcmcOutput` is the output from our previous call of the sampler...
"

# ╔═╡ ae0503ce-111e-435a-9f6a-a1010d707880
plot(x)

# ╔═╡ Cell order:
# ╠═80b90b23-c999-4836-ae83-d1b7017daa46
# ╟─0439d40e-768e-44d0-85bc-417b063f4315
# ╟─8dda732d-6cd5-4ac9-a749-277b1ba49cea
# ╠═09b424d8-b54d-4d12-9367-c0413c46475e
# ╟─48ad33b4-b44c-482e-aa98-baa5a3e890ce
# ╠═5f990069-ba98-4930-9806-3ed8c25e7253
# ╠═699e9669-b105-4d5a-962b-7bf63acb8b49
# ╟─1504fd48-9e5e-4536-a0ed-7ebc26126b92
# ╠═508ba9f3-fc56-4012-a447-8bb6c675f271
# ╟─ef5fa55b-198c-41a6-acb4-e0bb964c822a
# ╠═ae0503ce-111e-435a-9f6a-a1010d707880


text:
"""
    AbstractQuantumSystem

Any implementation of the `AbstractQuantumSystem` interface, used for
self-consistent field calculations, must implement the functions
[`coefficients`](@ref), [`orbitals`](@ref), [`diff`](@ref),
[`normalize!`](@ref).
"""
abstract type AbstractQuantumSystem end

"""
    coefficients(quantum_system)

Retrieves the mixing coefficients of `quantum_system`. Must return a
`view` that the [`scf!`](@ref) routine can modify. To be overloaded by
the implementation of [`AbstractQuantumSystem`](@ref).
"""
coefficients(::QS) where {QS<:AbstractQuantumSystem} =
    throw(ArgumentError("`coefficients` not implemented for $(QS)"))

"""
    orbitals(quantum_system)

Retrieves the orbitals of `quantum_system`. Must return a
`view` that the [`scf!`](@ref) routine can modify. To be overloaded by
the implementation of [`AbstractQuantumSystem`](@ref).
"""
orbitals(::QS) where {QS<:AbstractQuantumSystem} =
    throw(ArgumentError("`orbitals` not implemented for $(QS)"))

"""
    diff(quantum_system[; kwargs...])

Varies the the `quantum_system` with respect to all orbitals. Used to
derive the multi-configurational Hartree–Fock equations. To be overloaded by
the implementation of [`AbstractQuantumSystem`](@ref).
"""
Base.diff(::QS; kwargs...) where {QS<:AbstractQuantumSystem} =
    throw(ArgumentError("`diff` not implemented for $(QS)"))

"""
    normalize!(quantum_system, v)

Normalize the orbital `v` of `quantum_system`. To be overloaded by the
implementation of [`AbstractQuantumSystem`](@ref).
"""
LinearAlgebra.normalize!(::QS, ::AbstractVector) where {QS<:AbstractQuantumSystem} =
    throw(ArgumentError("`normalize!` not implemented for `$QS`"))

export AbstractQuantumSystem


text:
function missing_cards(col::Collection, subtarget::SingleSetDeckArray)
    for rarity in (4,3,2,1)
        if missing_cards(col[rarity][subtarget.set], subtarget[rarity])
            return true
        end
    end
    return false
end

function missing_cards(col::Collection, subtarget::SingleSetDeckArray, rarity::Int)
    if missing_cards(col[rarity][subtarget.set], subtarget[rarity])
        return true
    end
    return false
end

function missing_cards(A::Vector{Int}, subset::Vector{Tuple{Int,Int}})
    for i in subset
        if A[i[1]] < i[2]
            return true
        end
    end
    return false
end

function count_missing(A::Vector{Int}, subset::Vector{Tuple{Int,Int}})
    missing = 0
    for i in subset
        missing += (i[2] - A[i[1]])
    end
    missing
end

function remaining_dust(col::Collection, x::SingleSetDeckArray)
    dust = 0
    for r in 1:4
        dust += count_missing(col[r][x.set], x[r]) * craft_regular[r]
    end
    dust
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
export UnitVector, UnitSimplex, CorrCholeskyFactor

####
#### building blocks
####

"""
    (y, r, ℓ) = $SIGNATURES

Given ``x ∈ ℝ`` and ``0 ≤ r ≤ 1``, return `(y, r′)` such that

1. ``y² + r′² = r²``,

2. ``y: |y| ≤ r`` is mapped with a bijection from `x`.

`ℓ` is the log Jacobian (whether it is evaluated depends on `flag`).
"""
@inline function l2_remainder_transform(flag::LogJacFlag, x, r)
    z = 2*logistic(x) - 1
    (z * √r, r*(1 - abs2(z)),
     flag isa NoLogJac ? flag : log(2) + logistic_logjac(x) + 0.5*log(r))
end

"""
    (x, r′) = $SIGNATURES

Inverse of [`l2_remainder_transform`](@ref) in `x` and `y`.
"""
@inline l2_remainder_inverse(y, r) = logit((y/√r+1)/2), r-abs2(y)

####
#### UnitVector
####

"""
    UnitVector(n)

Transform `n-1` real numbers to a unit vector of length `n`, under the
Euclidean norm.
"""
@calltrans struct UnitVector <: VectorTransform
    n::Int
    function UnitVector(n::Int)
        @argcheck n ≥ 1 "Dimension should be positive."
        new(n)
    end
end

dimension(t::UnitVector) = t.n - 1

function transform_with(flag::LogJacFlag, t::UnitVector, x::AbstractVector, index)
    @unpack n = t
    T = extended_eltype(x)
    r = one(T)
    y = Vector{T}(undef, n)
    ℓ = logjac_zero(flag, T)
    @inbounds for i in 1:(n - 1)
        xi = x[index]
        index += 1
        y[i], r, ℓi = l2_remainder_transform(flag, xi, r)
        ℓ += ℓi
    end
    y[end] = √r
    y, ℓ, index
end

inverse_eltype(t::UnitVector, y::AbstractVector) = extended_eltype(y)

function inverse_at!(x::AbstractVector, index, t::UnitVector, y::AbstractVector)
    @unpack n = t
    @argcheck length(y) == n
    r = one(eltype(y))
    @inbounds for yi in axes(y, 1)[1:(end-1)]
        x[index], r = l2_remainder_inverse(y[yi], r)
        index += 1
    end
    index
end

####
#### UnitSimplex
####

"""
    UnitSimplex(n)

Transform `n-1` real numbers to a vector of length `n` whose elements are non-negative and sum to one.
"""
@calltrans struct UnitSimplex <: VectorTransform
    n::Int
    function UnitSimplex(n::Int)
        @argcheck n ≥ 1 "Dimension should be positive."
        new(n)
    end
end

dimension(t::UnitSimplex) = t.n - 1

function transform_with(flag::LogJacFlag, t::UnitSimplex, x::AbstractVector, index)
    @unpack n = t
    T = extended_eltype(x)

    ℓ = logjac_zero(flag, T)
    stick = one(T)
    y = Vector{T}(undef, n)
    @inbounds for i in 1:n-1
        xi = x[index]
        index += 1
        z = logistic(xi - log(n-i))
        y[i] = z * stick

        if !(flag isa NoLogJac)
            ℓ += log(stick) - logit_logjac(z)
        end

        stick *= 1 - z
    end

    y[end] = stick

    y, ℓ, index
end

inverse_eltype(t::UnitSimplex, y::AbstractVector) = extended_eltype(y)

function inverse_at!(x::AbstractVector, index, t::UnitSimplex, y::AbstractVector)
    @unpack n = t
    @argcheck length(y) == n

    stick = one(eltype(y))
    @inbounds for i in axes(y, 1)[1:end-1]
        z = y[i]/stick
        x[index] = logit(z) + log(n-i)
        stick -= y[i]
        index += 1
    end
    index
end

####
#### correlation cholesky factor
####

"""
    CorrCholeskyFactor(n)

Cholesky factor of a correlation matrix of size `n`.

Transforms ``n×(n-1)/2`` real numbers to an ``n×n`` upper-triangular matrix `U`, such that
`U'*U` is a correlation matrix (positive definite, with unit diagonal).

# Notes

If

- `z` is a vector of `n` IID standard normal variates,

- `σ` is an `n`-element vector of standard deviations,

- `U` is obtained from `CorrCholeskyFactor(n)`,

then `Diagonal(σ) * U' * z` will be a multivariate normal with the given variances and
correlation matrix `U' * U`.
"""
@calltrans struct CorrCholeskyFactor <: VectorTransform
    n::Int
    function CorrCholeskyFactor(n)
        @argcheck n ≥ 1 "Dimension should be positive."
        new(n)
    end
end

dimension(t::CorrCholeskyFactor) = unit_triangular_dimension(t.n)

function transform_with(flag::LogJacFlag, t::CorrCholeskyFactor, x::AbstractVector, index)
    @unpack n = t
    T = extended_eltype(x)
    ℓ = logjac_zero(flag, T)
    U = Matrix{T}(undef, n, n)
    @inbounds for col in 1:n
        r = one(T)
        for row in 1:(col-1)
            xi = x[index]
            U[row, col], r, ℓi = l2_remainder_transform(flag, xi, r)
            ℓ += ℓi
            index += 1
        end
        U[col, col] = √r
    end
    UpperTriangular(U), ℓ, index
end

inverse_eltype(t::CorrCholeskyFactor, U::UpperTriangular) = extended_eltype(U)

function inverse_at!(x::AbstractVector, index, t::CorrCholeskyFactor, U::UpperTriangular)
    @unpack n = t
    @argcheck size(U, 1) == n
    @inbounds for col in 1:n
        r = one(eltype(U))
        for row in 1:(col-1)
            x[index], r = l2_remainder_inverse(U[row, col], r)
            index += 1
        end
    end
    index
end


text:
#!/usr/bin/julia

# Copyright (C) 2021 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Eleven, Dumbo Octopus."

include("Utils.jl")

using .Utils

"Perform a step and update the `octopuses` energy grid."
function update_octopuses!(octopuses :: Matrix{Int})
    new_energy = ones(Int, size(octopuses))
    flashed = Set{CartesianIndex{2}}([])

    has_flashed = true
    while has_flashed
        has_flashed = false

        for ind in CartesianIndices(octopuses)
            if octopuses[ind] + new_energy[ind] >= 10 && ind ∉ flashed
                push!(flashed, ind)
                has_flashed = true

                neighbors = find_all_neighbors(octopuses, ind.I)
                [new_energy[neighbor...] += 1 for neighbor in neighbors]
            end
        end
    end

    octopuses .+= new_energy

    [octopuses[i] = 0 for i in flashed]

    return length(flashed)
end

"Calculate the number of all `octopuses`’ flashes up to `steps`."
function solve1(octopuses :: Matrix{Int}, steps :: Int)
    octs = copy(octopuses)
    flashes = 0

    for _ in 1:steps
        flashes += update_octopuses!(octs)
    end

    return flashes
end

"Return the first step that all `octopuses` flash together."
function solve2(octopuses :: Matrix{Int})
    octs = copy(octopuses)
    step = 1

    while update_octopuses!(octs) != length(octs)
        step += 1
    end

    return step
end

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readlines(first(ARGS))

    octopuses = hcat([parse.(Int, e) for e in collect.(data)]...)

    println(solve1(octopuses, 100))
    println(solve2(octopuses))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
begin
  groups = split(strip(read("input.txt", String)), "\n\n")

  groupanswers = group -> length(filter(isletter, unique(group)))

  println(mapreduce(groupanswers, +, groups))

  function sharedanswers(group)
    g = string.(split(group, '\n'))
    if length(g) == 1
      return groupanswers(g[1])
    else
      return count(letter -> all(contains(letter), g[2:end]), g[1])
    end
  end

  println(mapreduce(sharedanswers, +, groups))
end


text:
#GracePlot demonstration 3: Generate proper plot & save/export
#(Uses template to avoid specifying too many parameters)
#-------------------------------------------------------------------------------

using GracePlot

#==Input data
===============================================================================#
x=collect(-1:0.02:1).*pi
y1 = sin(x)
y2 = cos(x)

#=="Defaults"
===============================================================================#
smallplot = GracePlot.template("smallplot_mono")

#==Generate plot
===============================================================================#
plot = GracePlot.new(guimode=true, template=smallplot)
g = add(plot, title = "Grace SVG Plot")
	set(g, subtitle = "(\\f{Times-Italic}y\\s1\\N=sin(x), y\\s2\\N=cos(x)\\f{})")
	set(g, xlabel = "Angle (rad)", ylabel = "Amplitude")
	#Add datasets:
		ds = add(g, x, y1)
		ds = add(g, x, y2)
	autofit(g)

#Finalize:
redraw(plot)

#Save plot in multiple formats:
GracePlot._write("sinewaveplot.agr", plot)
GracePlot.write_eps("sinewaveplot.eps", plot)
GracePlot.write_svg("sinewaveplot.svg", plot)
GracePlot.write_png("sinewaveplot.png", plot, dpi=300)
#Last line


text:
# squared Euclidean distance from a set

immutable SqrDistL2{R <: Real} <: ProximableFunction
  ind::IndicatorConvex
  lambda::R
  function SqrDistL2(ind::IndicatorConvex, lambda::R)
    if lambda < 0
      error("parameter λ must be nonnegative")
    else
      new(ind, lambda)
    end
  end
end

SqrDistL2{R <: Real}(ind::IndicatorConvex, lambda::R=1.0) = SqrDistL2{R}(ind, lambda)

function (f::SqrDistL2){T <: RealOrComplex}(x::AbstractArray{T})
  p, = prox(f.ind, x)
  return (f.lambda/2)*vecnorm(x-p)^2
end

function prox!{T <: RealOrComplex}(f::SqrDistL2, x::AbstractArray{T}, y::AbstractArray{T}, gamma::Real=1.0)
  p, = prox(f.ind, x)
  sqrd = (f.lambda/2)*vecnorm(x-p)^2
  c1 = 1/(1+f.lambda*gamma)
  c2 = f.lambda*gamma*c1
  for k in eachindex(p)
    y[k] = c1*x[k] + c2*p[k]
  end
  return sqrd*c1^2
end

fun_name(f::SqrDistL2) = "squared Euclidean distance from a convex set"
fun_dom(f::SqrDistL2) = fun_dom(f.ind)
fun_expr(f::SqrDistL2) = "x ↦ (λ/2) inf { ||x-y||^2 : y ∈ S} "
fun_params(f::SqrDistL2) = string("λ = $(f.lambda), S = ", typeof(f.ind))

function prox_naive{T <: RealOrComplex}(f::SqrDistL2, x::AbstractArray{T}, gamma::Real=1.0)
  p, = prox(f.ind, x)
  sqrd = (f.lambda/2)*vecnorm(x-p)^2
  gamlam = f.lambda*gamma
  return 1/(1+gamlam)*x + gamlam/(1+gamlam)*p, sqrd/(1+gamlam)^2
end


text:
function solution_check_unit_tests()

end


text:
mutable struct NeuronResponse
    name::String
    get_spikes::Function
end

function get_spikes_from_distribution(d)
    N -> (filter(x->first(τ)<x<last(τ), sort(rand(d,2N)[1:N])))
end

# smooth ramp
function ramp(τ)
    mt = maximum(τ)
    f(t) = begin
        if t < 0
            return 0.1
        else
            return t
        end
    end
    [f(t/mt) for t in τ]
end

# oscillate after
function oscillate(τ)
    mt = maximum(τ)
    f(t) = begin
        if t < 0
            return 0.5
        else
            return sin(π*t)+0.5
        end
    end
    [f(t) for t in τ]
end

# short response and refractory period
function refractory(τ)
    mt = maximum(τ)
    f(t) = begin
        if t < 0
            return 0.5
        elseif t < 0.05
            return 0.95
        elseif t < 0.1
            return 0.05
        else
            return 0.5
        end
    end
    [f(t) for t in τ]
end


text:
promote_rule(::Type{DoubleFloat{Float64}}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float64}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float32}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{Float32}
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{Float16}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{T}) where {T<:Integer} = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{T}) where {T<:Integer} = DoubleFloat{Float32}
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{T}) where {T<:Integer} = DoubleFloat{Float16}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{BigInt}) = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{BigInt}) = DoubleFloat{Float32}
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{BigInt}) = DoubleFloat{Float16}
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{Float64}}) = DoubleFloat{Float64}
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{Float32}
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float16}

promote_rule(::Type{DoubleFloat{Float64}}, ::Type{BigFloat}) = DoubleFloat{Float64}
promote_rule(::Type{DoubleFloat{Float32}}, ::Type{BigFloat}) = DoubleFloat{Float32}
promote_rule(::Type{DoubleFloat{Float16}}, ::Type{BigFloat}) = DoubleFloat{Float16}
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{Float64}}) = DoubleFloat{Float64}
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{Float32}
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{Float16}

promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{DoubleFloat{DoubleFloat{Float32}}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{DoubleFloat{DoubleFloat{Float16}}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{DoubleFloat{DoubleFloat{Float16}}}) = DoubleFloat{DoubleFloat{Float32}}

promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{DoubleFloat{Float64}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{DoubleFloat{Float64}}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{DoubleFloat{Float64}}) = DoubleFloat{DoubleFloat{Float16}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{DoubleFloat{Float32}}) = DoubleFloat{DoubleFloat{Float16}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{DoubleFloat{Float16}}) = DoubleFloat{DoubleFloat{Float16}}

promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{T}) where {T<:IEEEFloat} = DoubleFloat{DoubleFloat{Float16}}

promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{T}) where {T<:Integer} = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{T}) where {T<:Integer} = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{T}) where {T<:Integer} = DoubleFloat{DoubleFloat{Float16}}

promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{BigInt}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{BigInt}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{BigInt}) = DoubleFloat{DoubleFloat{Float16}}

promote_rule(::Type{BigInt}, ::Type{DoubleFloat{DoubleFloat{Float64}}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{DoubleFloat{Float32}}}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{BigInt}, ::Type{DoubleFloat{DoubleFloat{Float16}}}) = DoubleFloat{DoubleFloat{Float16}}

promote_rule(::Type{DoubleFloat{DoubleFloat{Float64}}}, ::Type{BigFloat}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float32}}}, ::Type{BigFloat}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{DoubleFloat{DoubleFloat{Float16}}}, ::Type{BigFloat}) = DoubleFloat{DoubleFloat{Float16}}

promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{DoubleFloat{Float64}}}) = DoubleFloat{DoubleFloat{Float64}}
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{DoubleFloat{Float32}}}) = DoubleFloat{DoubleFloat{Float32}}
promote_rule(::Type{BigFloat}, ::Type{DoubleFloat{DoubleFloat{Float16}}}) = DoubleFloat{DoubleFloat{Float16}}