text:
using LazyFym
using Transducers

# using Lazy
using InfiniteArrays
using StructArrays
using Random
using Plots, LaTeXStrings

function ẋ(env::LazyFym.InputAffineQuadraticCostEnv, x, t)
    u = command(env, x)
    ẋ = LazyFym.ẋ(env, x, t, u)
end
_env = LazyFym.InputAffineQuadraticCostEnv()
command(env, x) = LazyFym.u_optimal(_env, x)  # you can customise it

function initial_condition(env::LazyFym.InputAffineQuadraticCostEnv)
    2*(rand(2) .- 0.5)
end

function postprocess(env::LazyFym.InputAffineQuadraticCostEnv)
    function _postprocess(_datum)
        t = _datum.t
        x = _datum.x
        x1 = _datum.x[1]
        x2 = _datum.x[2]
        u = command(env, x)
        datum = (; t = t, x = x, x1=x1, x2=x2, u = u)
    end
end

function single()
    Random.seed!(1)
    env = LazyFym.InputAffineQuadraticCostEnv()
    t0 = 0.0
    t1 = 5.00
    Δt = 0.01
    ts = t0:Δt:∞
    x0 = initial_condition(env)
    sim(x0) = t -> Sim(env, x0, ts, ẋ) |> TakeWhile(datum -> datum.t <= t) |> Map(postprocess(env)) |> collect |> StructArray
    traj_x0 = sim(x0)
    data = traj_x0(t1)
    # data = @lazy traj_x0(t1);  # for lazy evaluation, see Lazy.jl
    l = @layout [a; b; c]
    p_x1 = plot(data.t, data.x1,
                ylabel=L"x_{1}", label=nothing, ylim=(-1.5, 1.5))
    p_x2 = plot(data.t, data.x2,
                ylabel=L"x_{2}", label=nothing, ylim=(-1.5, 1.5))
    p_u = plot(data.t, data.u,
               xlabel=L"t", ylabel=L"u", label=nothing, ylim=(-3, 3))
    p = plot(p_x1, p_x2, p_u, layout = l)
    savefig(p, "figures/single.png")
end

function parallel()
    Random.seed!(1)
    env = LazyFym.InputAffineQuadraticCostEnv()
    t0 = 0.0
    t1 = 5.00
    Δt = 0.01
    ts = t0:Δt:∞
    num = 10
    x0s = 1:num |> Map(i -> initial_condition(env))
    traj(x0) = Sim(env, x0, ts, ẋ) |> TakeWhile(datum -> datum.t <= t1) |> Map(postprocess(env)) |> collect
    data_parallel = x0s |> Map(x0 -> traj(x0)) |> Map(StructArray) |> tcollect
    # data_parallel_whole = data_parallel |> TCat(Threads.nthreads()) |> collect |> StructArray   # merge data
    l = @layout [a; b; c]
    p_x1 = plot()
    _ = data_parallel |> Map(data -> plot!(p_x1, data.t, data.x1,
                                           ylabel=L"x_{1}", label=nothing,
                                           ylim=(-1.5, 1.5))) |> collect
    p_x2 = plot()
    _ = data_parallel |> Map(data -> plot!(p_x2, data.t, data.x2,
                                           ylabel=L"x_{2}", label=nothing,
                                           ylim=(-1.5, 1.5))) |> collect
    p_u = plot()
    _ = data_parallel |> Map(data -> plot!(p_u, data.t, data.u,
                                           xlabel=L"t", ylabel=L"u", label=nothing,
                                           ylim=(-3, 3))) |> collect
    p = plot(p_x1, p_x2, p_u, layout = l)
    savefig(p, "figures/parallel.png")
end

# single()
# parallel()


text:
using Documenter
using NKFtool

makedocs(
    sitename = "NKFtool.jl",
    authors = "Hiroharu Sugawara <hsugawa@gmail.com>",
    format = Documenter.HTML(lang="ja"),
    modules = [NKFtool],
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide"=>"man/guide.md",
            "手引き"=>"man/guideja.md",],
        "Library" => Any["Public"=>"lib/public.md",],
    ],
    repo="https://github.com/hsugawa8651/NKFtool.jl/blob/{commit}{path}#L{line}",
    #assets=String[],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(;
    repo="github.com/hsugawa8651/NKFtool.jl",
)

ext: jl


text:
module CombUtils

export swap!
function swap!{T<:Any, N<:Integer}(a::Array{T, 1}, i::N, j::N)
    tmp = a[i]
    a[i] = a[j]
    a[j] = tmp
end

export next_permutation!
# Precondition: Input is array consists of unique integers, initially sorted in
# strictly ascending order
function next_permutation!{N<:Integer}(ar::Array{N, 1})
    # 1. Find the highest index, i1 such that ar[i1] is the first of a pair of
    # elements in ascending order. If there isn't one, the sequence is the
    # highest permutation, so reverse the whole thing to begin again.
    alen = length(ar)
    i1 = alen - 1
    while i1 > 0
        if ar[i1] < ar[i1 + 1]
            break
        end
        i1 -= 1
    end
    if i1 > 0
        # 2. Find the highest index i2, such that i2 > i1 and ar[i2] > ar[i1].
        i2 = alen
        while i2 > i1
            if ar[i2] > ar[i1]
                break
            end
            i2 -= 1
        end
        # 3. Swap ar[i1] and ar[i2].
        swap!(ar, i1, i2)
        # 4. The elements from ar[i1 + 1] to the end are now in descending order
        # (a later permutation), so reverse them.
        reverse!(ar, i1 + 1, alen)
    else
        reverse!(ar, 1, alen)
    end
end

end


text:
#!/usr/bin/env julia

# this test isn't run by Pkg.test

# it can be run to test Juno/Jupyter integration...

using Luxor

using Test

using Random
Random.seed!(42)

function draw_hypotrochoids()
    background("grey10")
    sethue("antiquewhite")
    setopacity(0.8)
    tiles = Tiler(1200, 1200, 6, 6)
    for (pos, n) in tiles
        gsave()
        setline(.5)
        translate(pos)

        R = 45; r = n; d = 20
        sethue("antiquewhite")
        poly(hypotrochoid(R, n, d, :stroke, stepby=pi/300, vertices=true), :stroke, close=true)
        text(string("R = $(R), r = $(r), d = $(d)"), halign=:center, O.x, O.y+tiles.tileheight/2)

        # test perimeterlength
        ptslist = hypotrochoid(R, n, d, :stroke, stepby=pi/300, vertices=true)
        polyperimeter(ptslist)

        R = 2n; r = 12; d = 5
        sethue("orange")
        (R == r) && (r += 1)
        p = hypotrochoid(R, r, d, :none, stepby=pi/305, vertices=true)
        polydistances(p) # just testing
        poly(p, :stroke, close=true)
        text(string("R = $(R), r = $(r), d = $(d)"), halign=:center, O.x, O.y+tiles.tileheight/2 + 10)

        # testing portions of polys
        portion = rand()
        setline(.25)
        sethue("magenta")
        section = polyportion(p, portion)
        poly(section, :stroke, close=false)
        sethue("green")
        circle(section[end], 2, :stroke)
        poly(polyremainder(p, portion), :stroke, close=false)

        text(string("portion ", round(portion, digits=2)), halign=:center, O.x, O.y+tiles.tileheight/2 + 30)

        sethue("cyan")
        setline(0.3)
        poly(offsetpoly(p, 15), :stroke, close=true)
        text("offset 15", halign=:center, O.x, O.y+tiles.tileheight/2 + 20)

        grestore()
    end
end

function test_hypotrochoid_at_png()
    @png begin
        draw_hypotrochoids()
    end 1200 1200
end
function test_hypotrochoid_at_svg()
    @svg begin
        draw_hypotrochoids()
    end 1200 1200
end
function test_hypotrochoid_at_pdf()
    @pdf begin
        draw_hypotrochoids()
    end 1200 1200
end

Random.seed!(42)

test_hypotrochoid_at_png() ; println("output PNG")
test_hypotrochoid_at_svg() ; println("output SVG")
test_hypotrochoid_at_pdf() ; println("output PDF")

#testing @manipulate

using Interact

function makecircle(r)
    d = Drawing(300, 300, :svg)
    sethue("black")
    setline(0.1)
    origin()
    hypotrochoid(150, r, 35, :stroke)
    finish()
    return d
end

@manipulate for r in 5:150
    makecircle(r)
end


text:
using SIMD

function chipsort_merge_medium(input::AbstractArray{T,1}, ::Val{V}, ::Val{J}, ::Val{K}) where {T,V,J,K}
    output = valloc(T, div(32,sizeof(T)), V*J*K)

    if J>1
        output .= input
        sort_blocks!(output, Val(J), Val(V))
    else
        for cc in 1:K
            block = ntuple(v->input[v+(cc-1)*V], Val(V))
            srt = Vec(sort_net(block...)) ::Vec{V*J,T}
            vstorea(srt, output, 1+(cc-1)*V*J)
        end
    end

    new_input = reshape(output, V, J, K)

    blocks_a = [(@view new_input[:,1:end,c*2-1]) for c in 1:K>>1]
    blocks_b = [(@view new_input[:,1:end,c*2]) for c in 1:K>>1]

    do_merge_pass(
        new_input,
        reshape(valloc(T, div(32,sizeof(T)), V*J*K), V, J<<1, K>>1),
        blocks_a, blocks_b,
        Val(V), Val(J<<1), Val(K>>1)
    )
end

@inline function do_merge_pass(input::AbstractArray{T,3}, output::AbstractArray{T,3}, blocks_a, blocks_b,::Val{V}, ::Val{J}, ::Val{K}) where {T,V,J,K}

    for c in 1:K
        output[:,1,c] .= blocks_a[c][:,1]
    end
    next_inputs = [pointer(c) for c in blocks_b]

    bitonic_merge_interleaved(
        output,
        next_inputs,
        Val(V), 1, Val(K)
    )

    for c in 1:K
        blocks_a[c] = (@view (blocks_a[c][:,2:end]))
        blocks_b[c] = (@view (blocks_b[c][:,2:end]))
    end

    for iter in 2:(J-1)
        for c in 1:K
            if length(blocks_a[c]) > 0 && (length(blocks_b[c]) == 0 || blocks_a[c][1,1] < blocks_b[c][1,1])
                next_inputs[c] = pointer(blocks_a[c], 1)
                blocks_a[c] = (@view (blocks_a[c][:,2:end]))
            else
                next_inputs[c] = pointer(blocks_b[c], 1)
                blocks_b[c] = (@view (blocks_b[c][:,2:end]))
            end
        end

        bitonic_merge_interleaved(
            output,
            next_inputs,
            Val(V), iter, Val(K)
        )
    end

    if K == 1
        reshape(output, :)
    else

        for c in 1:K>>1
            blocks_a[c] = @view output[:,1:end,c*2-1]
            blocks_b[c] = @view output[:,1:end,c*2]
        end

        do_merge_pass(
            output,
            reshape(input, V,J<<1,K>>1),
            blocks_a, blocks_b,
            Val(V), Val(J<<1), Val(K>>1)
        )
    end
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module Nervoso

import Base: ==, show, length, size, endof, getindex, setindex!,
             start, next, done, tanh

include("utils.jl")
include("feedforward/layer.jl")
include("activation.jl")
include("cost.jl")
include("feedforward/network.jl")
include("feedforward/performance.jl")

function __init__()
    # Put every derivative into a Dict in import-time
    derivatives[tanh]      = tanhprime
    derivatives[softmax]   = softmaxprime
    derivatives[quaderror] = quaderrorprime
    derivatives[ceerror]   = ceerrorprime
end

end # module Nervoso


text:
# Silhouette

# this function returns r of size (k, n), such that
# r[i, j] is the sum of distances of all points from cluster i to sample j
#
function sil_aggregate_dists{T<:Real}(k::Int, a::AbstractVector{Int}, dists::DenseMatrix{T})
    n = length(a)
    r = zeros(k, n)
    for j = 1:n
        (1 <= a[j] <= k) || error("a[j] should have 1 <= a[j] <= k.")
    end
    @inbounds for j = 1:n
        for i = 1:j-1
            r[a[i],j] += dists[i,j]
        end
        for i = j+1:n
            r[a[i],j] += dists[i,j]
        end
    end
    return r
end

function silhouettes{T<:Real}(assignments::Vector{Int},
                              counts::AbstractVector{Int},
                              dists::DenseMatrix{T})

    n = length(assignments)
    k = length(counts)
    size(dists) == (n, n) || throw(DimensionMismatch("Inconsistent array dimensions."))

    # compute average distance from each cluster to each point --> r
    r = sil_aggregate_dists(k, assignments, dists)
    # from sum to average
    @inbounds for j = 1:n
        for i = 1:k
            c = counts[i]
            if i == assignments[j]
                c -= 1
            end
            if c == 0
                r[i,j] = 0.0
            else
                r[i,j] /= c
            end
        end
    end

    # compute a and b
    # a: average distance w.r.t. the assigned cluster
    # b: the minimum average distance w.r.t. other cluster
    a = Array(Float64, n)
    b = Array(Float64, n)

    for j = 1:n
        l = assignments[j]
        a[j] = r[l, j]

        v = Inf
        p = -1
        for i = 1:k
            @inbounds rij = r[i,j]
            if (i != l) && (rij < v)
                v = rij
                p = i
            end
        end
        b[j] = v
    end

    # compute silhouette score
    sil = a   # reuse the memory of a for sil
    for j = 1:n
        @inbounds sil[j] = (b[j] - a[j]) / max(a[j], b[j])
    end
    return sil
end

silhouettes(R::ClusteringResult, dists::DenseMatrix) =
    silhouettes(assignments(R), counts(R), dists)


text:
import BinaryProvider: satisfied
export BuildStep, Dependency, satisfied, build

"""
    Dependency

A `Dependency` represents a set of `Product`s that must be satisfied before a
package can be run.  These `Product`s can be libraries, basic files,
executables, etc...

To build a `Dependency`, construct it and use `build()`.  To check to see if it
is already satisfied, use `satisfied()`.
"""
struct Dependency
    # The "name" of this dependency (e.g. "cairo")
    name::String

    # The resultant objects that must be present by this recipe's Products
    # have different rules that are automatically applied when verifying that
    # this dependency is satisfied
    results::Vector{Product}

    # The bash script to build this dependency
    script::String

    # The platform this dependency is built for
    platform::Platform

    # The parent "prefix" this dependency must be installed within
    prefix::Prefix

    """
        Dependency(name::AbstractString,
                   results::Vector{Product},
                   cmds::Vector{Cmd},
                   platform::Platform,
                   prefix::Prefix = global_prefix)

    Defines a new dependency that must be built by its name, the binary objects
    that will be built, the commands that must be run to build them, and the
    prefix within which all of this will be installed.

    You may also provide a `Vector` of tuples containing the URLs and hashes
    """
    function Dependency(name::AbstractString,
                    results::Vector{P},
                    script::AbstractString,
                    platform::Platform,
                    prefix::Prefix = BinaryProvider.global_prefix) where {P <: Product}
        return new(name, results, script, platform, prefix)
    end
end

"""
    satisfied(dep::Dependency; verbose::Bool = false, isolate::Bool = true)

Return true if all results are satisfied for this dependency.
"""
function satisfied(dep::Dependency; verbose::Bool = false, isolate::Bool = true)
    s = result -> satisfied(result; verbose=verbose, platform=dep.platform, isolate=isolate)
    return all(s(result) for result in dep.results)
end

"""
    build(runner, dep::Dependency; verbose::Bool = false, force::Bool = false,
          autofix::Bool = false, ignore_audit_errors::Bool = true,
          debug::Bool = false)

Build the dependency for given `platform` (defaulting to the host platform)
unless it is already satisfied.  If `force` is set to `true`, then the
dependency is always built.  Runs an audit of the built files, printing out
warnings if hints of unrelocatability are found.  These warnings are, by
default, ignored, unless `ignore_audit_errors` is set to `false`.  Some
warnings can be automatically fixed, and this will be attempted if `autofix` is
set to `true`.
"""
function build(runner, dep::Dependency; verbose::Bool = false, force::Bool = false,
               autofix::Bool = false, ignore_audit_errors::Bool = true,
               ignore_manifests::Vector = [], debug::Bool = false)
    # First, look to see whether this dependency is satisfied or not
    should_build = !satisfied(dep; isolate=true)

    # If it is not satisfied, (or we're forcing the issue) build it
    if force || should_build
        # Verbose mode tells us what's going on
        if verbose
            if !should_build
                Compat.@info("Force-building $(dep.name) despite its satisfaction")
            else
                Compat.@info("Building $(dep.name) as it is unsatisfied")
            end
        end

        # This `bash trap` snippet causes each command to get written out to bash
        # history so that, if we fail and fall into a debug shell, we can just
        # up-arrow to get to the last run command.
        trapped_script = """
        save_history() {
            if [[ "$(verbose)" == "true" ]]; then
                echo " ---> \$BASH_COMMAND" >&2
            fi
             history -s "\$BASH_COMMAND"
            history -a
        }

        save_env() {
            set +x
            set > /meta/.env
            # Ignore read-only variables
            for l in BASHOPTS BASH_VERSINFO UID EUID PPID SHELLOPTS; do
                grep -v "^\$l=" /meta/.env > /meta/.env2
                mv /meta/.env2 /meta/.env
            done
            echo "cd \$(pwd)" >> /meta/.env
        }

        # If /meta is mounted, then we want to save history and environment
        if [[ -d /meta ]]; then
            trap save_history DEBUG
            trap save_env INT TERM EXIT ERR
        fi

        # Stop if we hit any errors.
        set -e

        $(dep.script)
        """

        logpath = joinpath(logdir(dep.prefix), "$(dep.name).log")
        did_succeed = run(runner, `/bin/bash -c $(trapped_script)`, logpath; verbose=verbose)
        if !did_succeed
            if debug
                Compat.@warn("Build failed, launching debug shell")
                run_interactive(runner, `/bin/bash --init-file /meta/.env`)
            end
            msg = "Build for $(dep.name) on $(triplet(dep.platform)) did not complete successfully\n"
            error(msg)
        end

        # Run an audit of the prefix to ensure it is properly relocatable
        audit_result = audit(dep.prefix; platform=dep.platform,
                             verbose=verbose, autofix=autofix,
                             ignore_manifests=ignore_manifests)
        if !audit_result && !ignore_audit_errors
            msg = replace("""
            Audit failed for $(dep.prefix.path).
            Address the errors above to ensure relocatability.
            To override this check, set `ignore_audit_errors = true`.
            """, '\n' => ' ')
            error(strip(msg))
        end

        # Finally, check to see if we are now satisfied
        if !satisfied(dep; verbose=verbose, isolate=true)
            if verbose
                Compat.@warn("Built $(dep.name) but still unsatisfied!")
            end
        end
    elseif !should_build && verbose
        Compat.@info("Not building as $(dep.name) is already satisfied")
    end
    return true
end


text:
function restrict(A::Poisson2{T}) where {T}
    interiorc = restrict(A.interior)
    sz = size(A.interior)
    szc = size(interiorc)
    dxc = ((sz .- 1) ./ (szc .- 1)) .* A.dx
    return Poisson2(T, interiorc, dxc)
end

function restrict(interior::AbstractArray{Bool, 3})
    sz = size(interior)
    szc = restrict_size(sz)
    return restrict!(tzero(interior, szc), interior)
end

function restrict!(mc::AbstractArray{Bool, 3}, m::AbstractArray{Bool, 3})
    nxc, nyc, nzc = size(mc)

    @inbounds @batch minbatch=8 for kc in 2:nzc-1
        k = (kc << 1) - 1
        for jc in 2:nyc-1
            j = (jc << 1) - 1
            for ic in 2:nxc-1
                i = (ic << 1) - 1

                mc[ic,jc,kc] =
                    m[i,j,k] &&
                    m[i+1,j,k] &&
                    m[i,j+1,k] &&
                    m[i+1,j+1,k] &&
                    m[i,j,k+1]   &&
                    m[i+1,j,k+1] &&
                    m[i,j+1,k+1] &&
                    m[i+1,j+1,k+1]
            end
        end
    end

    return mc
end

function restrict!(
    xc::AbstractArray{<:AbstractFloat, 3},
    Ac::Poisson2{<:AbstractFloat, 3},
    x::AbstractArray{T, 3},
) where {T<:AbstractFloat}
    a0 = convert(T, 8//64)
    a1 = convert(T, 4//64)
    a2 = convert(T, 2//64)
    a3 = convert(T, 1//64)

    @inbounds @batch minbatch=8 for (Ic, Jc, Kc) in Ac.R27
        for kc in Kc
            k = (kc << 1) - 1
            for jc in Jc
                j = (jc << 1) - 1
                for ic in Ic
                    i = (ic << 1) - 1

                    x000 = x[i-1,j-1,k-1]
                    x100 = x[i,j-1,k-1]
                    x200 = x[i+1,j-1,k-1]

                    x010 = x[i-1,j,k-1]
                    x110 = x[i,j,k-1]
                    x210 = x[i+1,j,k-1]

                    x020 = x[i-1,j+1,k-1]
                    x120 = x[i,j+1,k-1]
                    x220 = x[i+1,j+1,k-1]

                    x001 = x[i-1,j-1,k]
                    x101 = x[i,j-1,k]
                    x201 = x[i+1,j-1,k]

                    x011 = x[i-1,j,k]
                    x111 = x[i,j,k]
                    x211 = x[i+1,j,k]

                    x021 = x[i-1,j+1,k]
                    x121 = x[i,j+1,k]
                    x221 = x[i+1,j+1,k]

                    x002 = x[i-1,j-1,k+1]
                    x102 = x[i,j-1,k+1]
                    x202 = x[i+1,j-1,k+1]

                    x012 = x[i-1,j,k+1]
                    x112 = x[i,j,k+1]
                    x212 = x[i+1,j,k+1]

                    x022 = x[i-1,j+1,k+1]
                    x122 = x[i,j+1,k+1]
                    x222 = x[i+1,j+1,k+1]

                    if Ac.interior[ic,jc,kc]
                        x3 = x000 + x200 + x020 + x220 +
                             x002 + x202 + x022 + x222

                        x2 = x100 + x010 + x210 + x120 +
                             x001 + x201 + x021 + x221 +
                             x102 + x012 + x212 + x122

                        x1 = x110 + x101 + x011 + x211 +
                             x121 + x112

                        c = a0 * x111
                        c = muladd(a3, x3, c)
                        c = muladd(a2, x2, c)
                        c = muladd(a1, x1, c)

                        xc[ic,jc,kc] = c
                    end
                end
            end
        end
    end

    return xc
end

function prolong!(
    x::AbstractArray{<:AbstractFloat, 3},
    A::Poisson2{<:AbstractFloat, 3},
    xc::AbstractArray{T, 3},
) where {T<:AbstractFloat}
    a1 = convert(T, 4//8)
    a2 = convert(T, 2//8)
    a3 = convert(T, 1//8)

    nxc, nyc, nzc = size(xc)
    m = A.interior

    @inbounds @batch for kc in 1:nzc-1
        k = (kc << 1) - 1
        for jc in 1:nyc-1
            j = (jc << 1) - 1
            for ic in 1:nxc-1
                i = (ic << 1) - 1

                x000 = xc[ic,jc,kc]
                x100 = xc[ic+1,jc,kc]
                x010 = xc[ic,jc+1,kc]
                x110 = xc[ic+1,jc+1,kc]
                x001 = xc[ic,jc,kc+1]
                x101 = xc[ic+1,jc,kc+1]
                x011 = xc[ic,jc+1,kc+1]
                x111 = xc[ic+1,jc+1,kc+1]

                x1 = x000 + x100
                x2 = x000 + x010
                x3 = x000 + x001

                x4 = x010 + x110
                x5 = x001 + x101
                x6 = x001 + x011

                if m[i,j,k]
                    x[i,j,k] = x000
                end

                if m[i+1,j,k]
                    x[i+1,j,k] = a1*x1
                end

                if m[i,j+1,k]
                    x[i,j+1,k] = a1*x2
                end

                if m[i+1,j+1,k]
                    x[i+1,j+1,k] = a2*(x1 + x4)
                end

                if m[i,j,k+1]
                    x[i,j,k+1] = a1*x3
                end

                if m[i+1,j,k+1]
                    x[i+1,j,k+1] = a2*(x1 + x5)
                end

                if m[i,j+1,k+1]
                    x[i,j+1,k+1] = a2*(x2 + x6)
                end

                if m[i+1,j+1,k+1]
                    x[i+1,j+1,k+1] = a3*(x1 + x4 + x5 + x011 + x111)
                end
            end
        end
    end

    return x
end

function correct_prolong!(
    x::AbstractArray{<:AbstractFloat, 3},
    A::Poisson2{<:AbstractFloat, 3},
    xc::AbstractArray{T, 3},
) where {T}
    a1 = convert(T, 4//8)
    a2 = convert(T, 2//8)
    a3 = convert(T, 1//8)

    nxc, nyc, nzc = size(xc)
    m = A.interior

    @inbounds @batch for kc in 1:nzc-1
        k = (kc << 1) - 1
        for jc in 1:nyc-1
            j = (jc << 1) - 1
            for ic in 1:nxc-1
                i = (ic << 1) - 1

                x000 = xc[ic,jc,kc]
                x100 = xc[ic+1,jc,kc]
                x010 = xc[ic,jc+1,kc]
                x110 = xc[ic+1,jc+1,kc]
                x001 = xc[ic,jc,kc+1]
                x101 = xc[ic+1,jc,kc+1]
                x011 = xc[ic,jc+1,kc+1]
                x111 = xc[ic+1,jc+1,kc+1]

                x1 = x000 + x100
                x2 = x000 + x010

                x4 = x010 + x110
                x5 = x001 + x101
                x6 = x001 + x011

                if m[i,j,k]
                    x[i,j,k] += x000
                end

                if m[i+1,j,k]
                    x[i+1,j,k] = muladd(a1, x1, x[i+1,j,k])
                end

                if m[i,j+1,k]
                    x[i,j+1,k] = muladd(a1, x2, x[i,j+1,k])
                end

                if m[i+1,j+1,k]
                    xx = x1 + x4
                    x[i+1,j+1,k] = muladd(a2, xx, x[i+1,j+1,k])
                end

                if m[i,j,k+1]
                    xx = x000 + x001
                    x[i,j,k+1] = muladd(a1, xx, x[i,j,k+1])
                end

                if m[i+1,j,k+1]
                    xx = x1 + x5
                    x[i+1,j,k+1] = muladd(a2, xx, x[i+1,j,k+1])
                end

                if m[i,j+1,k+1]
                    xx = x2 + x6
                    x[i,j+1,k+1] = muladd(a2, xx, x[i,j+1,k+1])
                end

                if m[i+1,j+1,k+1]
                    xx = x1 + x4 + x5 + x011 + x111
                    x[i+1,j+1,k+1] = muladd(a3, xx, x[i+1,j+1,k+1])
                end
            end
        end
    end

    return x
end

restrict_size(n::Integer) = (n + 1) >> 1
restrict_size(ns::Tuple) = restrict_size.(ns)
restrict_size(ns...) = restrict_size(ns)


text:
import Iterators: product

"Helper type holding information about a fixture definition"
type Fixture
 s::Symbol
 f::Function
 fargs::Array{Symbol, 1}  # all dependency-fixtures' symbols
 fixtures_dict  # maps dependency-fixtures' symbols to their resp. Fixtures
 kwargs
end

"Fixture's params as an array of (Symbol, param) pairs"
function parametrization_as_pairs(fixture::Fixture)
  [(fixture.s, param) for param in fixture.kwargs[:params]]
end

"Gets a matrix of parametrizations (cartesian product) for an iterable of fixtures"
function get_param_matrix(fixtures)
  consumed = Set{Symbol}()
  parametrized = Array{Fixture, 1}()
  sift_for_parametrized_fixtures!(fixtures, parametrized, consumed)
  product([parametrization_as_pairs(f) for f in parametrized]...)
end

# FIXME: get rid of recursions?
"Searches recursively an iterable of fixtures in search of parametrized ones"
function sift_for_parametrized_fixtures!(fixtures, parametrized::Array{Fixture, 1},
                                         consumed::Set{Symbol})
  for f in fixtures
    if :params in keys(f.kwargs) && !(f.s in consumed)
      push!(parametrized, f)
      push!(consumed, f.s)
    end
    sift_for_parametrized_fixtures!(values(f.fixtures_dict), parametrized, consumed)
  end
end


text:
import Knet.Ops21: elu, gelu, hardsigmoid, hardswish, relu, selu, sigm, swish, tanh_
import Knet.Ops21: eluback, geluback, hardsigmoidback, hardswishback, reluback, seluback, sigmback, swishback, tanh_back
import Base.Broadcast: broadcasted
import Knet
using Knet.KnetArrays: KnetArray, Bcasted
using CUDA: CUDA, CuArray, CuPtr
using Knet.LibKnet8: @knet8

for (A,P) in ((KnetArray,Ptr), (CuArray,CuPtr)), T in (Float32,Float64); S = sizeof(T) * 8
    for f in ("elu", "gelu", "hardsigmoid", "hardswish", "selu", "sigm", "swish", "tanh_")
        J, Jback = Symbol(f), Symbol("$(f)back")
        M, Mback = which(@__MODULE__,J), which(@__MODULE__,Jback)
        F, Fback = "$(f)_$S", "$(f)back_$(S)_111"
        @eval begin

            function ($M).$J(x::$A{$T})
                y = similar(x)
                @knet8($F,(Cint,$P{$T},$P{$T}),length(y),x,y)
                return y
            end

            function ($Mback).$Jback(x::X, y::X, dy::X) where {X<:$A{$T}}
                dx = similar(x)
                @knet8($Fback,(Cint,$P{$T},$P{$T},$P{$T},$P{$T}),length(dx),x,y,dy,dx)
                return dx
            end
        end
    end

    for f in ("relu",)
        J, Jback = Symbol(f), Symbol("$(f)back")
        M, Mback = which(@__MODULE__,J), which(@__MODULE__,Jback)
        F, Fback = "$(f)_$(S)_1", "$(f)back_$(S)_1"
        @eval begin

            function ($M).$J(x::$A{$T}; max_value=Inf, negative_slope=0, threshold=0)
                y = similar(x)
                @knet8($F,(Cint,$T,$T,$T,$P{$T},$P{$T}),length(y),max_value,negative_slope,threshold,x,y)
                return y
            end

            function ($Mback).$Jback(x::X, y::X, dy::X; max_value=Inf, negative_slope=0, threshold=0) where {X<:$A{$T}}
                dx = similar(x)
                @knet8($Fback,(Cint,$T,$T,$T,$P{$T},$P{$T},$P{$T},$P{$T}),length(dx),max_value,negative_slope,threshold,x,y,dy,dx)
                return dx
            end
        end
    end
end


text:
using Makie

function main()
    s = Scene()
    x = -10:0.01:-0.1
    lines!(s, x, 1 ./ x)
    x = 0.1:0.01:10
    lines!(s, x, 1 ./ x)
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Toxcore

type JuliaBot
    quit::Bool
    svg_file_bytes::Vector{Uint8}
end

function JuliaBot()
    JuliaBot(false, [])
end

bot = JuliaBot()

function OnToxFileChunkRequest(tox::Ptr{Tox}, friend_number::Uint32, file_number::Uint32, position::Uint64, chunk_length::Csize_t, user_data::Ptr{Void})
    global bot

    chunk_length = max(chunk_length, length(bot.svg_file_bytes)-position)
    data = bot.svg_file_bytes[position+1:position+chunk_length]

    if tox_file_send_chunk(tox, friend_number, file_number, position, data)
        println("File chunk send.")
    else
        println("Failed to send file chunk.")
    end

    nothing
end

function OnToxFileRecv(tox::Ptr{Tox}, friend_number::Uint32, file_number::Uint32, kind::TOX_FILE_KIND, file_size::Uint64, filename::Ptr{Uint8}, filename_length::Csize_t, user_data::Ptr{Void})
    println("File Recv Callback")

    nothing
end

function OnToxFileRecvChunk(tox::Ptr{Tox}, friend_number::Uint32, file_number::Uint32, position::Uint64, data::Ptr{Uint8}, data_length::Csize_t, user_data::Ptr{Void})
    println("File Recv Chunk Callback")

    nothing
end

function OnToxFileRecvControl(tox::Ptr{Tox}, friend_number::Uint32, file_number::Uint32, control::TOX_FILE_CONTROL, user_data::Ptr{Void})
    if control == TOX_FILE_CONTROL_RESUME
        println("Resumming file transfer.")
    elseif control == TOX_FILE_CONTROL_PAUSE
        println("File transfer paused.")
    else # TOX_FILE_CONTROL_CANCEL
        println("File transfer cancelled.")
    end

    nothing
end

function OnToxFriendRequest(tox::Ptr{Tox}, public_key::Ptr{Uint8}, message::Ptr{Uint8}, length::Csize_t, user_data::Ptr{Void})
    tox_friend_add_norequest(tox, ToxPublicKey(public_key))
    println("Accepted a friend request.")
    nothing
end

function OnToxFriendMessage(tox::Ptr{Tox}, friend_number::Uint32, typ::TOX_MESSAGE_TYPE, message::Ptr{Uint8}, message_length::Csize_t, user_data::Ptr{Void})
    global bot

    msg = utf8(pointer_to_array(message, message_length))
    println(msg)

    if (msg == "/h")
        tox_friend_send_message(tox, friend_number, utf8("Hi, I'm julia!"))
        tox_friend_send_message(tox, friend_number, utf8("I know these commands:"))
        tox_friend_send_message(tox, friend_number, utf8(" '/h' - some help"))
        tox_friend_send_message(tox, friend_number, utf8(" '/f' - I send you a file"))
        tox_friend_send_message(tox, friend_number, utf8(" '/q' - I quit"))
    elseif (msg == "/q")
        tox_friend_send_message(tox, friend_number, utf8("Ok, I’ll quit!"))
        bot.quit = true
    elseif msg == "/f"
        #tox_friend_send_message(tox, friend_number, utf8("Ok dude, I'll send you a nice file!"))

        #if tox_file_send(tox, friend_number, TOX_FILE_KIND_DATA, length(bot.svg_file_bytes), utf8("julia.svg")) == typemax(Uint32)
        #    println("Attempt to send file failed.")
        #end
    else
        tox_friend_send_message(tox, friend_number, utf8("Thanks for your message, dude!"))
    end

    nothing
end

function main()
    global bot

    info("This is the Julia Tox bot")

    # Load the file that the bot always sends
    julia_svg_file = open(Pkg.dir("Toxcore", "test/julia.svg"), "r")
    bot.svg_file_bytes = readbytes(julia_svg_file)
    close(julia_svg_file)

    # Try to load the tox settings
    my_tox = 0

    try
        savefile = open(Pkg.dir("Toxcore", "test", "bot_savedata.binary"), "r")
        savedata = readbytes(savefile)
        close(savefile)

        default_options = tox_options_default()

        options = Tox_Options(default_options.ipv6_enabled,
                            default_options.udp_enabled,
                            default_options.proxy_type,
                            default_options.proxy_host,
                            default_options.proxy_port,
                            default_options.start_port,
                            default_options.end_port,
                            default_options.tcp_port,
                            TOX_SAVEDATA_TYPE_TOX_SAVE,
                            pointer(savedata),
                            length(savedata) )

        my_tox = tox_new(options)

        info("Previous bot instance found. Reusing it!")
    catch e
        println(e)
        #Create a default Tox
        my_tox = tox_new()

        info("Created new bot instance")
    end

    # register the callbacks
    tox_callback_file_chunk_request(my_tox, OnToxFileChunkRequest, C_NULL)
    tox_callback_friend_request(my_tox, OnToxFriendRequest, C_NULL)
    tox_callback_friend_message(my_tox, OnToxFriendMessage, C_NULL)
    tox_callback_file_recv(my_tox, OnToxFileRecv, C_NULL)
    tox_callback_file_recv_chunk(my_tox, OnToxFileRecvChunk, C_NULL)
    tox_callback_file_recv_control(my_tox, OnToxFileRecvControl, C_NULL)

    # print own address
    info("Here is my address")
    println(convert(ASCIIString, tox_self_get_address(my_tox)))

    # define user details
	tox_self_set_name(my_tox, "Julia")
    tox_self_set_status_message(my_tox, "I am Julia, a high-level, high-performance dynamic programming language for technical computing.")

    # bootstrap from the node defined above
    if !tox_bootstrap(my_tox)
        println("Failed to bootstrap.")
        exit()
    end

    # get the friend list
    friendlist = tox_self_get_friend_list(my_tox)
    info("I have $(length(friendlist)) friend(s)")

    for friend in friendlist
        println(tox_friend_get_name(my_tox, friend))
    end

    while true
        tox_iterate(my_tox)
        sleep(tox_iteration_interval(my_tox)/1000)

        if bot.quit
            break
        end
    end

    savedata = tox_get_savedata(my_tox)
    tox_kill(my_tox)

    savefile = open(Pkg.dir("Toxcore", "test", "bot_savedata.binary"), "w")
    write(savefile, savedata)
    close(savefile)

    info("Bot instance saved")
end

main()


text:
export kde

import ..KDEModel

using mlpack._Internal.io

import mlpack_jll
const kdeLibrary = mlpack_jll.libmlpack_julia_kde

# Call the C binding of the mlpack kde binding.
function kde_mlpackMain()
  success = ccall((:kde, kdeLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module kde_internal
  import ..kdeLibrary

import ...KDEModel

# Get the value of a model pointer parameter of type KDEModel.
function IOGetParamKDEModel(paramName::String)::KDEModel
  KDEModel(ccall((:IO_GetParamKDEModelPtr, kdeLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type KDEModel.
function IOSetParamKDEModel(paramName::String, model::KDEModel)
  ccall((:IO_SetParamKDEModelPtr, kdeLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeKDEModel(stream::IO, model::KDEModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeKDEModelPtr, kdeLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeKDEModel(stream::IO)::KDEModel
  buffer = read(stream)
  KDEModel(ccall((:DeserializeKDEModelPtr, kdeLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    kde(; [abs_error, algorithm, bandwidth, initial_sample_size, input_model, kernel, mc_break_coef, mc_entry_coef, mc_probability, monte_carlo, query, reference, rel_error, tree, verbose])

This program performs a Kernel Density Estimation. KDE is a non-parametric way
of estimating probability density function. For each query point the program
will estimate its probability density by applying a kernel function to each
reference point. The computational complexity of this is O(N^2) where there are
N query points and N reference points, but this implementation will typically
see better performance as it uses an approximate dual or single tree algorithm
for acceleration.

Dual or single tree optimization avoids many barely relevant calculations (as
kernel function values decrease with distance), so it is an approximate
computation. You can specify the maximum relative error tolerance for each query
value with `rel_error` as well as the maximum absolute error tolerance with the
parameter `abs_error`. This program runs using an Euclidean metric. Kernel
function can be selected using the `kernel` option. You can also choose what
which type of tree to use for the dual-tree algorithm with `tree`. It is also
possible to select whether to use dual-tree algorithm or single-tree algorithm
using the `algorithm` option.

Monte Carlo estimations can be used to accelerate the KDE estimate when the
Gaussian Kernel is used. This provides a probabilistic guarantee on the the
error of the resulting KDE instead of an absolute guarantee.To enable Monte
Carlo estimations, the `monte_carlo` flag can be used, and success probability
can be set with the `mc_probability` option. It is possible to set the initial
sample size for the Monte Carlo estimation using `initial_sample_size`. This
implementation will only consider a node, as a candidate for the Monte Carlo
estimation, if its number of descendant nodes is bigger than the initial sample
size. This can be controlled using a coefficient that will multiply the initial
sample size and can be set using `mc_entry_coef`. To avoid using the same amount
of computations an exact approach would take, this program recurses the tree
whenever a fraction of the amount of the node's descendant points have already
been computed. This fraction is set using `mc_break_coef`.

For example, the following will run KDE using the data in `ref_data` for
training and the data in `qu_data` as query data. It will apply an Epanechnikov
kernel with a 0.2 bandwidth to each reference point and use a KD-Tree for the
dual-tree optimization. The returned predictions will be within 5% of the real
KDE value for each query point.

```julia
julia> using CSV
julia> ref_data = CSV.read("ref_data.csv")
julia> qu_data = CSV.read("qu_data.csv")
julia> _, out_data = kde(bandwidth=0.2, kernel="epanechnikov",
            query=qu_data, reference=ref_data, rel_error=0.05, tree="kd-tree")
```

the predicted density estimations will be stored in `out_data`.
If no `query` is provided, then KDE will be computed on the `reference` dataset.
It is possible to select either a reference dataset or an input model but not
both at the same time. If an input model is selected and parameter values are
not set (e.g. `bandwidth`) then default parameter values will be used.

In addition to the last program call, it is also possible to activate Monte
Carlo estimations if a Gaussian kernel is used. This can provide faster results,
but the KDE will only have a probabilistic guarantee of meeting the desired
error bound (instead of an absolute guarantee). The following example will run
KDE using a Monte Carlo estimation when possible. The results will be within a
5% of the real KDE value with a 95% probability. Initial sample size for the
Monte Carlo estimation will be 200 points and a node will be a candidate for the
estimation only when it contains 700 (i.e. 3.5*200) points. If a node contains
700 points and 420 (i.e. 0.6*700) have already been sampled, then the algorithm
will recurse instead of keep sampling.

```julia
julia> using CSV
julia> ref_data = CSV.read("ref_data.csv")
julia> qu_data = CSV.read("qu_data.csv")
julia> _, out_data = kde(bandwidth=0.2, initial_sample_size=200,
            kernel="gaussian", mc_break_coef=0.6, mc_entry_coef=3.5,
            mc_probability=0.95, monte_carlo=, query=qu_data,
            reference=ref_data, rel_error=0.05, tree="kd-tree")
```

# Arguments

 - `abs_error::Float64`: Relative error tolerance for the prediction.
      Default value `0`.

 - `algorithm::String`: Algorithm to use for the prediction.('dual-tree',
      'single-tree').  Default value `dual-tree`.

 - `bandwidth::Float64`: Bandwidth of the kernel.  Default value `1`.

 - `initial_sample_size::Int`: Initial sample size for Monte Carlo
      estimations.  Default value `100`.

 - `input_model::KDEModel`: Contains pre-trained KDE model.
 - `kernel::String`: Kernel to use for the prediction.('gaussian',
      'epanechnikov', 'laplacian', 'spherical', 'triangular').  Default value
      `gaussian`.

 - `mc_break_coef::Float64`: Controls what fraction of the amount of
      node's descendants is the limit for the sample size before it recurses.
      Default value `0.4`.

 - `mc_entry_coef::Float64`: Controls how much larger does the amount of
      node descendants has to be compared to the initial sample size in order to
      be a candidate for Monte Carlo estimations.  Default value `3`.

 - `mc_probability::Float64`: Probability of the estimation being bounded
      by relative error when using Monte Carlo estimations.  Default value
      `0.95`.

 - `monte_carlo::Bool`: Whether to use Monte Carlo estimations when
      possible.  Default value `false`.

 - `query::Array{Float64, 2}`: Query dataset to KDE on.
 - `reference::Array{Float64, 2}`: Input reference dataset use for KDE.
 - `rel_error::Float64`: Relative error tolerance for the prediction.
      Default value `0.05`.

 - `tree::String`: Tree to use for the prediction.('kd-tree', 'ball-tree',
      'cover-tree', 'octree', 'r-tree').  Default value `kd-tree`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::KDEModel`: If specified, the KDE model will be saved
      here.
 - `predictions::Array{Float64, 1}`: Vector to store density predictions.

"""
function kde(;
             abs_error::Union{Float64, Missing} = missing,
             algorithm::Union{String, Missing} = missing,
             bandwidth::Union{Float64, Missing} = missing,
             initial_sample_size::Union{Int, Missing} = missing,
             input_model::Union{KDEModel, Missing} = missing,
             kernel::Union{String, Missing} = missing,
             mc_break_coef::Union{Float64, Missing} = missing,
             mc_entry_coef::Union{Float64, Missing} = missing,
             mc_probability::Union{Float64, Missing} = missing,
             monte_carlo::Union{Bool, Missing} = missing,
             query = missing,
             reference = missing,
             rel_error::Union{Float64, Missing} = missing,
             tree::Union{String, Missing} = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, kdeLibrary), Nothing, ());

  IORestoreSettings("Kernel Density Estimation")

  # Process each input argument before calling mlpackMain().
  if !ismissing(abs_error)
    IOSetParam("abs_error", convert(Float64, abs_error))
  end
  if !ismissing(algorithm)
    IOSetParam("algorithm", convert(String, algorithm))
  end
  if !ismissing(bandwidth)
    IOSetParam("bandwidth", convert(Float64, bandwidth))
  end
  if !ismissing(initial_sample_size)
    IOSetParam("initial_sample_size", convert(Int, initial_sample_size))
  end
  if !ismissing(input_model)
    kde_internal.IOSetParamKDEModel("input_model", convert(KDEModel, input_model))
  end
  if !ismissing(kernel)
    IOSetParam("kernel", convert(String, kernel))
  end
  if !ismissing(mc_break_coef)
    IOSetParam("mc_break_coef", convert(Float64, mc_break_coef))
  end
  if !ismissing(mc_entry_coef)
    IOSetParam("mc_entry_coef", convert(Float64, mc_entry_coef))
  end
  if !ismissing(mc_probability)
    IOSetParam("mc_probability", convert(Float64, mc_probability))
  end
  if !ismissing(monte_carlo)
    IOSetParam("monte_carlo", convert(Bool, monte_carlo))
  end
  if !ismissing(query)
    IOSetParamMat("query", query, points_are_rows)
  end
  if !ismissing(reference)
    IOSetParamMat("reference", reference, points_are_rows)
  end
  if !ismissing(rel_error)
    IOSetParam("rel_error", convert(Float64, rel_error))
  end
  if !ismissing(tree)
    IOSetParam("tree", convert(String, tree))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("predictions")
  # Call the program.
  kde_mlpackMain()

  return kde_internal.IOGetParamKDEModel("output_model"),
         IOGetParamCol("predictions")
end


text:
"""
    DataEntity

Abstract type for an SMLM data entity that has an index, channel, and x/y/z coordinates.
"""
abstract type DataEntity
end

extractcoordinates(entity::DataEntity) = [entity.x; entity.y; entity.z]
extractcoordinates(entities::Vector{<:DataEntity}) =
    mapreduce(extractcoordinates, hcat, entities) |> (length(entities) > 1 ? identity : make2d)
make2d(array1d::AbstractVector) = reshape(array1d, length(array1d), 1)


text:
"""
Description
===========

An instance of `Inversion` holds the current state of an inversion. This
includes the population, the current generation, the boundaries of the inversion
and the inversion parameters.

Fields
======

- `data::Matrix`: The data being inverted. Columns should consist of:
    1. Period
    2. Apparent resistivity
    3. Phase
    4. Apparent resistivity uncertainty
    5. Phase uncertainty
- `popsize::Integer`: The target size of the population.
- `depth_lims::Array{LayerBC}`: Array of boundaries for each layer's depth.
- `res_lims::Array{LayerBC}`: Array of boundaries for each layer's resistivity.
- `nelitists::Integer`: Number of elitist clones to make per generation.
- `pmutation::Real`: Probability of mutation. 0 <= Pm <= 1.
- `pselection::Real`: Probability of selection. 0 <= Ps <= 1.
- `tourn_size::Integer`: Number of competitors in tournament selection
- `η::Integer`: Distribution index for SBX. Must be non-negative.
- `pop::Array{Model}`: Population for the genetic algorithm. Sorted by fitness.
- `gen::Integer`: Current generation. Starts at 1.

Usage
=====

- `I = Inversion(data, popsize, depth_lims, res_lims)`
- `I = Inversion(data, popsize, depth_lims, res_lims, nelitists, pmutation, pselection, tourn_size, η)`

Notes
=====

- If `nelitists` is unspecified, it defaults to 5% of the population.
- If `pmutation` is unspecified, it defaults to 0.01.
- If `pselection` is unspecified, it defaults to 0.65.
- If `tourn_size` is unspecified, it defaults to 2.
- If `η` is unspecified, it defaults to 2.
"""
type Inversion
    data::Matrix
    popsize::Integer
    depth_lims::Array{LayerBC}
    res_lims::Array{LayerBC}

    nelitists::Integer
    pmutation::Real
    pselection::Real
    tourn_size::Integer

    η::Integer

    pop::Array{Model}
    gen::Int

    function Inversion(data::Matrix, popsize::Integer, depth_lims::Array{LayerBC},
                       res_lims::Array{LayerBC})
        Inversion(data, popsize, depth_lims, res_lims, ceil(Integer, 0.05 * popsize), 0.01, 0.65, 2, 2)
    end

    function Inversion(data::Matrix, popsize::Integer, depth_lims::Array{LayerBC},
                       res_lims::Array{LayerBC}, nelitists::Integer, pmutation::Real, pselection::Real,
                       tourn_size::Int, η::Integer)
        if length(depth_lims) != length(res_lims)
            error("Length of depth_lims != length of res_lims")
        end

        nLayers = length(depth_lims)
        pop = vec(Array{Model}(popsize,1))
        for i in 1:popsize
            pop[i] = rand_model(nLayers, depth_lims, res_lims)
            calcfitness!(pop[i], data)
        end
        sortpop!(pop)
        new(data, popsize, depth_lims, res_lims, nelitists, pmutation, pselection, tourn_size, η, pop, 1)
    end
end

"""
Description
===========

Evolve an instance of `Inversion` through a specified number of
generations. This involves carrying out all the steps of the genetic
algorithm. Evolution will be terminated if the fitness falls within the stopping
conditions.

Arguments
=========

- `I::Inversion`: Inversion instance to evolve.
- `ngen::Integer=100`: Number of generations to advance `I` by.
- `verbose::Bool=false`: Print the best RMS misfit every 100 generations.

Side Effects
============

- Mutates the `pop` and `gen` fields of `I`.

Returns
=======

- `nothing`
"""
function evolve!(I::Inversion, ngen=100; verbose=false)
    if verbose
        @printf("Gen: %d\tBest RMS:%.2f\tSize: %d\n", I.gen, I.pop[1].fitness, length(I.pop))
    end
    for i in 1:ngen
        make_new_generation!(I)

        for idx in 1:I.popsize
            mutate!(I.pop[idx], I.pmutation)
            calcfitness!(I.pop[idx], I.data)
        end

        sortpop!(I.pop)
        I.gen += 1

        if 0.5 <= I.pop[1].fitness <= 1.5
            break
        end

        if verbose && I.gen % 100 == 0
            @printf("Gen: %d\tBest RMS:%.2f\tSize: %d\n", I.gen, I.pop[1].fitness, length(I.pop))
        end
    end
    return
end

"""
Description
===========

Create the next generation for an inversion by making elitist clones, selecting
parents and mating the selections.

Arguments
=========

- `I::Inversion`: Inversion instance to create the next generation for.

Side Effects
============

- Mutates the `pop` field of `I`.

Returns
=======

- `nothing`
"""
function make_new_generation!(I::Inversion)
    nextgen = I.pop[1:I.nelitists]

    # Create new children to fill the next generation
    while length(nextgen) < I.popsize
        p1 = select_with_tournament(I.pop, I.tourn_size, I.pselection)
        p2 = select_with_tournament(I.pop, I.tourn_size, I.pselection)

        childA, childB = crossover(p1, p2, I.η)
        calcfitness!(childA, I.data)
        calcfitness!(childB, I.data)
        push!(nextgen, childA, childB)
    end

    # At most, the next generation can have 1 child too many. Kill
    # this child if this is the case.
    if length(nextgen) > I.popsize
        I.pop = nextgen[1:end-1]
    else
        I.pop = nextgen
    end

    return
end


text:
mark(x) = ifelse(iseven(x), "#", ".")
function print_board(H, W)
    for i in 1:H
        for j in 1:W
            print(mark(i+j))
        end
        println()
    end
end

while true
    H, W = parse.(Int, split(readline()))
    (H == 0 && W == 0) && break
    print_board(H, W)
    println()
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Unary operators ##

conj!(A::AbstractQuasiArray{<:Number}) = (@inbounds broadcast!(conj, A, A); A)

for f in (:-, :conj, :real, :imag)
    @eval ($f)(A::AbstractQuasiArray) = broadcast_preserving_zero_d($f, A)
end

## Binary arithmetic operators ##

for f in (:+, :-)
    @eval function ($f)(A::AbstractQuasiArray, B::AbstractQuasiArray)
        promote_shape(A, B) # check size compatibility
        broadcast_preserving_zero_d($f, A, B)
    end
end

function +(A::QuasiArray, Bs::QuasiArray...)
    for B in Bs
        promote_shape(A, B) # check size compatibility
    end
    broadcast_preserving_zero_d(+, A, Bs...)
end

for f in (:/, :\, :*)
    if f != :/
        @eval ($f)(A::Number, B::AbstractQuasiArray) = broadcast_preserving_zero_d($f, A, B)
    end
    if f != :\
        @eval ($f)(A::AbstractQuasiArray, B::Number) = broadcast_preserving_zero_d($f, A, B)
    end
end


text:
module VariantCallFormat

import Automa
import Automa.RegExp: @re_str
import BioCore
import BioCore: isfilled, metainfotag, metainfoval, header
import BioCore.Exceptions: missingerror
import BufferedStreams

export
    VCF,
    BCF,
    header,
    metainfotag,
    metainfoval,
    isfilled,
    MissingFieldException

include("record.jl")
include("metainfo.jl")
include("header.jl")
include("reader.jl")
include("writer.jl")

# This module is here for backwards compatibility and to serve as a way to conveniently access unexported names
module VCF
	import VariantCallFormat

	const alt = VariantCallFormat.alt
	const chrom = VariantCallFormat.chrom
	const filter = VariantCallFormat.filter
	const format = VariantCallFormat.format
	const genotype = VariantCallFormat.genotype
	const hasalt = VariantCallFormat.hasalt
	const haschrom = VariantCallFormat.haschrom
	const hasfilter = VariantCallFormat.hasfilter
	const hasformat = VariantCallFormat.hasformat
	const hasid = VariantCallFormat.hasid
	const hasinfo = VariantCallFormat.hasinfo
	const haspos = VariantCallFormat.haspos
	const hasqual = VariantCallFormat.hasqual
	const hasref = VariantCallFormat.hasref
	const Header = VariantCallFormat.Header
	const id = VariantCallFormat.id
	const info = VariantCallFormat.info
	const infokeys = VariantCallFormat.infokeys
	const initialize! = VariantCallFormat.initialize!
	const isfilled = VariantCallFormat.isfilled
	const MetaInfo = VariantCallFormat.MetaInfo
	const pos = VariantCallFormat.pos
	const qual = VariantCallFormat.qual
	const Reader = VariantCallFormat.Reader
	const Record = VariantCallFormat.Record
	const ref = VariantCallFormat.ref
	const Writer = VariantCallFormat.Writer
end

include("bcf/bcf.jl")

end


text:
const cache = Dict{AbstractString, Function}()

"""
low-level API; template string + argnames -> (args -> rendered string)
"""
function compile(x::AbstractString, argnames::Symbol...) # TODO: escape `"""`
    raw = split(x, "\$\$")
    buf = IOBuffer(length(x.data) + length(raw)*32 + 128)

    iscode = false # indicating if i is code

    arglist = join(("$x=throw(error(\"missing required argument $x\"))" for x in argnames), ',')

    write(buf, "(_out::IO; ", arglist, ") -> begin\n")

    for i in raw
        if iscode
            write(buf, '\n', i, '\n')
        else
            write(buf, "write(_out, \"\"\"", i, "\"\"\")")
        end

        iscode = !iscode
    end

    write(buf, "\n_out\nend")

    buf |> takebuf_string |> parse |> eval # eval in the Restful module
end

"""
render with cache; automatically determin argnames by the arguments of first call
filepath and args -> rendered string
"""
function render(file::AbstractString; kwargs...) # maybe we can use @generated to perform cache
    if haskey(cache, file)
        buf = IOBuffer()
        cache[file](buf; kwargs...) |> takebuf_string
    else
        cache[file] = compile(open(readstring, file), map(x->x[1], kwargs)...)
        render(file; kwargs...)
    end
end


text:
module NetCDFIO

using NCDatasets
const NC = NCDatasets

using Statistics
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions

using ..ReferenceModels
using ..ReferenceStats
using ..Diagnostics
using ..HelperFuncs

export NetCDFIO_Diags
export open_files, close_files, write_iteration
export init_iteration_io, init_particle_diags, init_metrics, init_ensemble_diags
export init_val_diagnostics
export io_prior, io_reference, io_diagnostics, io_val_diagnostics

mutable struct NetCDFIO_Diags
    root_grp::NC.NCDataset{Nothing}
    ensemble_grp::NC.NCDataset{NC.NCDataset{Nothing}}
    particle_grp::NC.NCDataset{NC.NCDataset{Nothing}}
    metric_grp::NC.NCDataset{NC.NCDataset{Nothing}}
    outdir_path::String
    filepath::String
    vars::Dict{String, Any} # Hack to avoid https://github.com/Alexander-Barth/NCDatasets.jl/issues/135
    function NetCDFIO_Diags(
        config::Dict{Any, Any},
        outdir_path::String,
        ref_stats::ReferenceStatistics,
        N_ens::IT,
        priors::ParameterDistribution,
        val_ref_stats::Union{ReferenceStatistics, Nothing} = nothing,
    ) where {IT <: Integer}

        # Initialize properties with valid type:
        tmp = tempname()
        root_grp = NC.Dataset(tmp, "c")
        NC.defGroup(root_grp, "ensemble_diags")
        NC.defGroup(root_grp, "particle_diags")
        NC.defGroup(root_grp, "metrics")
        ensemble_grp = root_grp.group["ensemble_diags"]
        particle_grp = root_grp.group["particle_diags"]
        metric_grp = root_grp.group["metrics"]
        close(root_grp)

        filepath = joinpath(outdir_path, "Diagnostics.nc")

        # Remove the NC file if it already exists.
        isfile(filepath) && rm(filepath; force = true)

        NC.Dataset(filepath, "c") do root_grp

            # Parameter dimension
            p = ndims(priors)
            # Output dimension for global problem: full and low-dim encoding
            d_full = full_length(ref_stats)
            d = pca_length(ref_stats)

            # Number of configurations
            C = length(ref_stats.pca_vec)

            # Max number of fields per configuration
            num_vars = [length(norm_scale) for norm_scale in ref_stats.norm_vec]
            f = maximum(num_vars)
            # Max number of variable degrees of freedom per configuration
            var_dof_max = maximum(Int.([size(P_pca, 1) for P_pca in ref_stats.pca_vec] ./ num_vars))

            # Number of configuration per batch
            batch_size = get_entry(config["reference"], "batch_size", length(ref_stats.pca_vec))
            batch_size = isnothing(batch_size) ? length(ref_stats.pca_vec) : batch_size

            particle = Array(1:N_ens)
            param = priors.name
            out_full = Array(1:d_full)
            out = Array(1:d)
            configuration = Array(1:C)
            field = Array(1:f)
            batch_index = Array(1:batch_size)
            dof = Array(1:var_dof_max)

            # Ensemble diagnostics (over all particles)
            ensemble_grp = NC.defGroup(root_grp, "ensemble_diags")
            NC.defDim(ensemble_grp, "param", p)
            NC.defVar(ensemble_grp, "param", param, ("param",))
            NC.defDim(ensemble_grp, "iteration", Inf)
            NC.defVar(ensemble_grp, "iteration", Int16, ("iteration",))

            # Parameter prior diagnostics
            prior_grp = NC.defGroup(root_grp, "prior")
            NC.defDim(prior_grp, "param", p)
            NC.defVar(prior_grp, "param", param, ("param",))

            # Reference model and stats diagnostics
            reference_grp = NC.defGroup(root_grp, "reference")
            NC.defDim(reference_grp, "out_full", d_full)
            NC.defVar(reference_grp, "out_full", out_full, ("out_full",))
            NC.defDim(reference_grp, "out", d)
            NC.defVar(reference_grp, "out", out, ("out",))
            NC.defDim(reference_grp, "config", C)
            NC.defVar(reference_grp, "config", configuration, ("config",))
            NC.defDim(reference_grp, "config_field", f)
            NC.defVar(reference_grp, "config_field", field, ("config_field",))
            NC.defDim(reference_grp, "dof", var_dof_max)
            NC.defVar(reference_grp, "dof", dof, ("dof",))

            augmented = get_entry(config["process"], "augmented", false)
            if augmented
                d = d + p
            end

            # Particle diagnostics
            particle_grp = NC.defGroup(root_grp, "particle_diags")
            NC.defDim(particle_grp, "particle", N_ens)
            NC.defVar(particle_grp, "particle", particle, ("particle",))
            NC.defDim(particle_grp, "out_full", d_full)
            NC.defVar(particle_grp, "out_full", out_full, ("out_full",))
            NC.defDim(particle_grp, "out_aug", d)
            NC.defVar(particle_grp, "out_aug", Array(1:d), ("out_aug",))
            NC.defDim(particle_grp, "param", p)
            NC.defVar(particle_grp, "param", param, ("param",))
            NC.defDim(particle_grp, "config", C)
            NC.defVar(particle_grp, "config", configuration, ("config",))
            NC.defDim(particle_grp, "batch_index", batch_size)
            NC.defVar(particle_grp, "batch_index", batch_index, ("batch_index",))
            NC.defDim(particle_grp, "iteration", Inf)
            NC.defVar(particle_grp, "iteration", Int16, ("iteration",))

            if !isnothing(val_ref_stats)
                d_full_val = full_length(val_ref_stats)
                d_val = pca_length(val_ref_stats)
                C_val = length(val_ref_stats.pca_vec)

                # Max number of fields per configuration
                num_vars_val = [length(norm_vec) for norm_vec in val_ref_stats.norm_vec]
                f_val = maximum(num_vars_val)
                # Max number of variable degrees of freedom per configuration
                var_dof_max_val = maximum(Int.([size(P_pca, 1) for P_pca in val_ref_stats.pca_vec] ./ num_vars_val))

                batch_size_val = get_entry(config["validation"], "batch_size", length(val_ref_stats.pca_vec))
                batch_size_val = isnothing(batch_size_val) ? length(val_ref_stats.pca_vec) : batch_size_val

                out_full_val = Array(1:d_full_val)
                out_val = Array(1:d_val)
                configuration_val = Array(1:C_val)
                field_val = Array(1:f_val)
                batch_index_val = Array(1:batch_size_val)
                dof_val = Array(1:var_dof_max_val)

                NC.defDim(reference_grp, "out_full_val", d_full_val)
                NC.defVar(reference_grp, "out_full_val", out_full_val, ("out_full_val",))
                NC.defDim(reference_grp, "out_val", d_val)
                NC.defVar(reference_grp, "out_val", out_val, ("out_val",))
                NC.defDim(reference_grp, "config_val", C_val)
                NC.defVar(reference_grp, "config_val", configuration_val, ("config_val",))
                NC.defDim(reference_grp, "config_field_val", f_val)
                NC.defVar(reference_grp, "config_field_val", field_val, ("config_field_val",))
                NC.defDim(reference_grp, "dof_val", var_dof_max_val)
                NC.defVar(reference_grp, "dof_val", dof_val, ("dof_val",))

                if augmented
                    d_val = d_val + p
                end

                NC.defDim(particle_grp, "out_full_val", d_full_val)
                NC.defVar(particle_grp, "out_full_val", out_full_val, ("out_full_val",))
                NC.defDim(particle_grp, "out_aug_val", d_val)
                NC.defVar(particle_grp, "out_aug_val", Array(1:d_val), ("out_aug_val",))
                NC.defDim(particle_grp, "batch_index_val", batch_size_val)
                NC.defVar(particle_grp, "batch_index_val", batch_index_val, ("batch_index_val",))
                NC.defDim(particle_grp, "config_val", C_val)
                NC.defVar(particle_grp, "config_val", configuration_val, ("config_val",))

            end

            # Calibration metrics
            metric_grp = NC.defGroup(root_grp, "metrics")
            NC.defDim(metric_grp, "iteration", Inf)
            NC.defVar(metric_grp, "iteration", Int16, ("iteration",))
        end
        vars = Dict{String, Any}()
        return new(root_grp, ensemble_grp, particle_grp, metric_grp, outdir_path, filepath, vars)
    end

    function NetCDFIO_Diags(filepath::String)
        vars = Dict{String, Any}()
        diags = nothing
        NC.Dataset(filepath, "r") do root_grp
            ensemble_grp = root_grp.group["ensemble_diags"]
            particle_grp = root_grp.group["particle_diags"]
            metric_grp = root_grp.group["metrics"]
            diags = new(root_grp, ensemble_grp, particle_grp, metric_grp, dirname(filepath), filepath, vars)
        end
        return diags
    end
end

function open_files(diags::NetCDFIO_Diags)
    diags.root_grp = NC.Dataset(diags.filepath, "a")
    diags.ensemble_grp = diags.root_grp.group["ensemble_diags"]
    diags.particle_grp = diags.root_grp.group["particle_diags"]
    diags.metric_grp = diags.root_grp.group["metrics"]
    vars = diags.vars

    # Hack to avoid https://github.com/Alexander-Barth/NCDatasets.jl/issues/135
    vars["ensemble_diags"] = Dict{String, Any}()
    for k in keys(diags.ensemble_grp)
        vars["ensemble_diags"][k] = diags.ensemble_grp[k]
    end
    vars["particle_diags"] = Dict{String, Any}()
    for k in keys(diags.particle_grp)
        vars["particle_diags"][k] = diags.particle_grp[k]
    end
    vars["metrics"] = Dict{String, Any}()
    for k in keys(diags.metric_grp)
        vars["metrics"][k] = diags.metric_grp[k]
    end
end

function close_files(diags::NetCDFIO_Diags)
    close(diags.root_grp)
end

"""Adds a given field to an existing NetCDF Dataset."""
function add_field(diags::NetCDFIO_Diags, var_name::String; dims, group, type)
    NC.Dataset(diags.filepath, "a") do root_grp
        grp = root_grp.group[group]
        new_var = NC.defVar(grp, var_name, type, dims)
    end
end

"""Writes current field `data` to an existing variable in a NetCDF Dataset."""
function write_current(diags::NetCDFIO_Diags, var_name::String, data; group)
    var = diags.vars[group][var_name]
    last_dim = length(size(var))
    last_dim_end = size(var, last_dim)
    try
        selectdim(var, last_dim, last_dim_end) .= data
    catch e
        @error string(
            "Failed to write array of dimension $(size(data)) as $var_name",
            " to NetCDF file. Expected array of dimension $(size(var)[1:end-1]).",
        )
        e
    end
end

"""Writes all current fields in a given io_dict to an existing NetCDF Dataset."""
function write_current_dict(diags::NetCDFIO_Diags, io_dict)
    for var in keys(io_dict)
        write_current(diags, var, io_dict[var].field; group = io_dict[var].group)
    end
end

function init_io_dict(diags::NetCDFIO_Diags, io_dict)
    for var in keys(io_dict)
        add_field(diags, var; dims = io_dict[var].dims, group = io_dict[var].group, type = io_dict[var].type)
    end
end

function io_prior(diags::NetCDFIO_Diags, prior::ParameterDistribution)
    io_dict = io_dictionary_prior(prior)
    for var in keys(io_dict)
        add_field(diags, var; dims = io_dict[var].dims, group = io_dict[var].group, type = io_dict[var].type)
        NC.Dataset(diags.filepath, "a") do root_grp
            prior_grp = root_grp.group["prior"]
            var_val = prior_grp[var]
            var_val .= io_dict[var].field
        end
    end
end

function write_ref(diags::NetCDFIO_Diags, var_name::String, data)
    NC.Dataset(diags.filepath, "a") do root_grp
        reference_grp = root_grp.group["reference"]
        var = reference_grp[var_name]
        var .= data
    end
end

function io_reference(
    diags::NetCDFIO_Diags,
    ref_stats::ReferenceStatistics,
    ref_models::Vector{ReferenceModel},
    write_full_stats::Bool = true,
)
    io_dict = io_dictionary_reference(ref_stats, ref_models, write_full_stats)
    for var in keys(io_dict)
        add_field(diags, var; dims = io_dict[var].dims, group = io_dict[var].group, type = io_dict[var].type)
        write_ref(diags, var, io_dict[var].field)
    end
end

function io_val_reference(
    diags::NetCDFIO_Diags,
    ref_stats::ReferenceStatistics,
    ref_models::Vector{ReferenceModel},
    write_full_stats::Bool = true,
)
    io_dict = io_dictionary_val_reference(ref_stats, ref_models, write_full_stats)
    for var in keys(io_dict)
        add_field(diags, var; dims = io_dict[var].dims, group = io_dict[var].group, type = io_dict[var].type)
        write_ref(diags, var, io_dict[var].field)
    end
end

function init_ensemble_diags(diags::NetCDFIO_Diags, ekp::EnsembleKalmanProcess, priors::ParameterDistribution)
    io_dict = io_dictionary_ensemble(ekp, priors)
    init_io_dict(diags, io_dict)
    # Write initial ensemble_state to file.
    open_files(diags)
    write_current_dict(diags, io_dict)
    close_files(diags)
end

function init_particle_diags(diags::NetCDFIO_Diags, ekp::EnsembleKalmanProcess, priors::ParameterDistribution)
    io_dict = io_dictionary_particle_eval()
    init_io_dict(diags, io_dict)
    io_dict = io_dictionary_particle_state(ekp, priors)
    init_io_dict(diags, io_dict)
    # Write initial particle_state to file.
    open_files(diags)
    write_current_dict(diags, io_dict)
    close_files(diags)
end

function init_metrics(diags::NetCDFIO_Diags)
    io_dict = io_dictionary_metrics()
    init_io_dict(diags, io_dict)
end

function init_val_metrics(diags::NetCDFIO_Diags)
    io_dict = io_dictionary_val_metrics()
    init_io_dict(diags, io_dict)
end

function init_val_particle_diags(diags::NetCDFIO_Diags)
    io_dict = io_dictionary_val_particle_eval()
    init_io_dict(diags, io_dict)
end

function init_val_diagnostics(
    diags::NetCDFIO_Diags,
    val_ref_stats::ReferenceStatistics,
    val_ref_models::Vector{ReferenceModel},
    write_full_stats::Bool = true,
)
    io_val_reference(diags, val_ref_stats, val_ref_models, write_full_stats)
    init_val_metrics(diags)
    init_val_particle_diags(diags)
end

function io_metrics(diags::NetCDFIO_Diags, ekp::EnsembleKalmanProcess, mse_full::Vector{FT}) where {FT <: Real}
    io_dict = io_dictionary_metrics(ekp, mse_full)
    write_current_dict(diags, io_dict)
end

function io_val_metrics(
    diags::NetCDFIO_Diags,
    ekp::EnsembleKalmanProcess,
    val_ref_stats::ReferenceStatistics,
    g_val::Matrix{FT},
    val_mse_full::Vector{FT},
) where {FT <: Real}
    io_dict = io_dictionary_val_metrics(ekp, val_ref_stats, g_val, val_mse_full)
    write_current_dict(diags, io_dict)
end

function io_ensemble_diags(diags::NetCDFIO_Diags, ekp::EnsembleKalmanProcess, priors::ParameterDistribution)
    io_dict = io_dictionary_ensemble(ekp, priors)
    write_current_dict(diags, io_dict)
end

function io_val_particle_diags(
    diags::NetCDFIO_Diags,
    mse_full::Vector{FT},
    g::Matrix{FT},
    g_full::Matrix{FT},
    batch_indices::Union{Vector{Int}, Nothing},
) where {FT <: Real}
    # Dimension of the outputs
    d_aug_val = length(diags.particle_grp["out_aug_val"])
    d_full_val = length(diags.particle_grp["out_full_val"])
    # If not minibatching - training set size
    batch_indices = isnothing(batch_indices) ? Array(diags.particle_grp["config_val"]) : batch_indices

    # Write eval diagnostics to file
    io_dict = io_dictionary_val_particle_eval(g, g_full, mse_full, d_aug_val, d_full_val, batch_indices)
    write_current_dict(diags, io_dict)
end

function io_particle_diags_eval(
    diags::NetCDFIO_Diags,
    ekp::EnsembleKalmanProcess,
    mse_full::Vector{FT},
    g_full::Matrix{FT},
    batch_indices::Union{Vector{Int}, Nothing},
) where {FT <: Real}
    # Dimension of the outputs
    d = length(diags.particle_grp["out_aug"])
    d_full = length(diags.particle_grp["out_full"])
    # If not minibatching - training set size
    batch_indices = isnothing(batch_indices) ? Array(diags.particle_grp["config"]) : batch_indices

    # Write eval diagnostics to file
    io_dict = io_dictionary_particle_eval(ekp, g_full, mse_full, d, d_full, batch_indices)
    write_current_dict(diags, io_dict)
end

function io_particle_diags_state(
    diags::NetCDFIO_Diags,
    ekp::EnsembleKalmanProcess,
    priors::ParameterDistribution,
) where {FT <: Real}
    io_dict = io_dictionary_particle_state(ekp, priors)
    write_current_dict(diags, io_dict)
end

function io_diagnostics(
    diags::NetCDFIO_Diags,
    ekp::EnsembleKalmanProcess,
    priors::ParameterDistribution,
    mse_full::Vector{FT},
    g_full::Matrix{FT},
    batch_indices::Union{Vector{Int}, Nothing},
) where {FT <: Real}
    open_files(diags)
    # Eval diagnostics
    io_metrics(diags, ekp, mse_full)
    io_particle_diags_eval(diags, ekp, mse_full, g_full, batch_indices)
    write_iteration(diags)
    # State diagnostics
    io_particle_diags_state(diags, ekp, priors)
    io_ensemble_diags(diags, ekp, priors)
    close_files(diags)
end

function io_val_diagnostics(
    diags::NetCDFIO_Diags,
    ekp::EnsembleKalmanProcess,
    mse_full_val::Vector{FT},
    g_val::Matrix{FT},
    g_full_val::Matrix{FT},
    val_ref_stats::ReferenceStatistics,
    val_batch_indices::Union{Vector{Int}, Nothing},
) where {FT <: Real}
    open_files(diags)
    io_val_metrics(diags, ekp, val_ref_stats, g_val, mse_full_val)
    io_val_particle_diags(diags, mse_full_val, g_val, g_full_val, val_batch_indices)
    close_files(diags)
end

function init_iteration_io(diags::NetCDFIO_Diags)
    open_files(diags)
    ensemble_t = diags.ensemble_grp["iteration"]
    @inbounds ensemble_t[1] = 0
    particle_t = diags.particle_grp["iteration"]
    @inbounds particle_t[1] = 0
    metric_t = diags.metric_grp["iteration"]
    @inbounds metric_t[1] = 0
    close_files(diags)
end

function write_iteration(diags::NetCDFIO_Diags)
    ensemble_t = diags.ensemble_grp["iteration"]
    @inbounds ensemble_t[end + 1] = ensemble_t[end] + 1

    particle_t = diags.particle_grp["iteration"]
    @inbounds particle_t[end + 1] = particle_t[end] + 1

    metric_t = diags.metric_grp["iteration"]
    @inbounds metric_t[end + 1] = metric_t[end] + 1
end

end # module


text:
module KernelFunctions

export kernelmatrix, kernelmatrix!, kernelmatrix_diag, kernelmatrix_diag!
export duplicate, set! # Helpers

export Kernel, MOKernel
export ConstantKernel, WhiteKernel, EyeKernel, ZeroKernel, WienerKernel
export CosineKernel
export SqExponentialKernel, RBFKernel, GaussianKernel, SEKernel
export LaplacianKernel, ExponentialKernel, GammaExponentialKernel
export ExponentiatedKernel
export FBMKernel
export MaternKernel, Matern12Kernel, Matern32Kernel, Matern52Kernel
export LinearKernel, PolynomialKernel
export RationalKernel, RationalQuadraticKernel, GammaRationalKernel
export PiecewisePolynomialKernel
export PeriodicKernel, NeuralNetworkKernel
export KernelSum, KernelProduct, KernelTensorProduct
export TransformedKernel, ScaledKernel, NormalizedKernel

export Transform,
    SelectTransform,
    ChainTransform,
    ScaleTransform,
    LinearTransform,
    ARDTransform,
    IdentityTransform,
    FunctionTransform,
    PeriodicTransform
export with_lengthscale

export NystromFact, nystrom

export gaborkernel
export spectral_mixture_kernel, spectral_mixture_product_kernel

export ColVecs, RowVecs

export MOInput, prepare_isotopic_multi_output_data
export IndependentMOKernel,
    LatentFactorMOKernel, IntrinsicCoregionMOKernel, LinearMixingModelKernel

# Reexports
export tensor, ⊗, compose

using Compat
using ChainRulesCore: ChainRulesCore, Tangent, ZeroTangent, NoTangent
using ChainRulesCore: @thunk, InplaceableThunk
using CompositionsBase
using Distances
using FillArrays
using Functors
using LinearAlgebra
using Requires
using SpecialFunctions: loggamma, besselk, polygamma
using StatsFuns: logtwo, twoπ, softplus
using StatsBase
using TensorCore
using ZygoteRules: ZygoteRules

abstract type Kernel end
abstract type SimpleKernel <: Kernel end

include("utils.jl")
include(joinpath("distances", "pairwise.jl"))
include(joinpath("distances", "dotproduct.jl"))
include(joinpath("distances", "delta.jl"))
include(joinpath("distances", "sinus.jl"))

include(joinpath("transform", "transform.jl"))
include(joinpath("transform", "scaletransform.jl"))
include(joinpath("transform", "ardtransform.jl"))
include(joinpath("transform", "lineartransform.jl"))
include(joinpath("transform", "functiontransform.jl"))
include(joinpath("transform", "selecttransform.jl"))
include(joinpath("transform", "chaintransform.jl"))
include(joinpath("transform", "periodic_transform.jl"))
include(joinpath("kernels", "transformedkernel.jl"))
include(joinpath("transform", "with_lengthscale.jl"))

include(joinpath("basekernels", "constant.jl"))
include(joinpath("basekernels", "cosine.jl"))
include(joinpath("basekernels", "exponential.jl"))
include(joinpath("basekernels", "exponentiated.jl"))
include(joinpath("basekernels", "fbm.jl"))
include(joinpath("basekernels", "gabor.jl"))
include(joinpath("basekernels", "matern.jl"))
include(joinpath("basekernels", "nn.jl"))
include(joinpath("basekernels", "periodic.jl"))
include(joinpath("basekernels", "piecewisepolynomial.jl"))
include(joinpath("basekernels", "polynomial.jl"))
include(joinpath("basekernels", "rational.jl"))
include(joinpath("basekernels", "sm.jl"))
include(joinpath("basekernels", "wiener.jl"))

include(joinpath("kernels", "scaledkernel.jl"))
include(joinpath("kernels", "normalizedkernel.jl"))
include(joinpath("matrix", "kernelmatrix.jl"))
include(joinpath("kernels", "kernelsum.jl"))
include(joinpath("kernels", "kernelproduct.jl"))
include(joinpath("kernels", "kerneltensorproduct.jl"))
include(joinpath("kernels", "overloads.jl"))
include(joinpath("kernels", "neuralkernelnetwork.jl"))
include(joinpath("approximations", "nystrom.jl"))
include("generic.jl")

include(joinpath("mokernels", "mokernel.jl"))
include(joinpath("mokernels", "moinput.jl"))
include(joinpath("mokernels", "independent.jl"))
include(joinpath("mokernels", "slfm.jl"))
include(joinpath("mokernels", "intrinsiccoregion.jl"))
include(joinpath("mokernels", "lmm.jl"))

include("chainrules.jl")
include("zygoterules.jl")

include("test_utils.jl")

function __init__()
    @require Kronecker = "2c470bb0-bcc8-11e8-3dad-c9649493f05e" begin
        include(joinpath("matrix", "kernelkroneckermat.jl"))
    end
    @require PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150" begin
        include(joinpath("matrix", "kernelpdmat.jl"))
    end
end

end


text:
using Documenter
using Bootstrap

makedocs()

deploydocs(
    repo = "github.com/juliangehring/Bootstrap.jl.git",
    julia  = "0.4",
    deps = Deps.pip("mkdocs", "python-markdown-math", "pygments", "mkdocs-material")
)


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

Makie.plottype(::Partition) = Viz{<:Tuple{Partition}}

function Makie.plot!(plot::Viz{<:Tuple{Partition}})
  # retrieve partition object
  partition = plot[:object][]

  colors = distinguishable_colors(length(partition))

  for (i, subset) in Iterators.enumerate(partition)
    # fallback to collection recipe
    viz!(plot, subset,
      size         = plot[:size],
      color        = colors[i],
      facetcolor   = plot[:boundarycolor],
      showfacets   = plot[:showboundary],
      decimation   = plot[:decimation],
    )
  end
end


text:
using RecyclePermanent
using RNGPool
using JLD2
using StaticArrays
using LinearAlgebra
using Random
import MonteCarloMarkovKernels: simulateChainProgress, makeAMKernel

include("../models/bpModel.jl")

@load "data/bpData.jld" ys θ
θ0 = θ

d = 3

@inline function toBPθ(v::SVector{d, Float64})
  return BPθ(v[1], v[2], v[3])
end

t0 = SVector{d, Float64}(θ0.λ, θ0.a, θ0.b)

sigmaProp = 1.0*SMatrix{d, d, Float64}(Matrix{Float64}(I, d, d))

@inline function lglogprior(θ::BPθ)
  if θ.λ < 0 return -Inf end
  if θ.a < 0 return -Inf end
  if θ.b < 0 return -Inf end
  v::Float64 = 0.0
  v -= 0.001*θ.λ
  v -= 0.1*θ.a
  v -= 0.1*θ.b
  return v
end

N = 40 * length(ys)
chainLength = 1000000


text:
struct ParameterServer
    params
end

function (ps::ParameterServer)(gs)
    for (p, g) in zip(ps.params, gs)
        p .-= g
    end
end

(ps::ParameterServer)() = deepcopy(p.params)

# Example usage
# ```julia
# ps = @actor ParameterServer([zeros(Float32, 3, 4), zeros(Float32, 3)])
# for c in clients
#     params = ps()[]
#     gs = calc_gradients(params)
#     ps[gs]
# end
# ```


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

module Math

using LinearAlgebra

export unitary_error, matrix_norm, eye, observe_unitary_error,
       min_observe_unitary_error

"Calculate the unitary matrix error"
unitary_error(exp_unit, mat_unit) = _unitary_error(exp_unit, mat_unit)
function unitary_error(exp_unit, mat_unit, exp_scale, mat_scale)
    return _unitary_error(exp_unit, mat_unit, x -> x * exp_scale, x -> x * mat_scale)
end
function _unitary_error(exp_unit, mat_unit, fa=identity, fb=identity)
    return sum(((a, b),) -> abs2(fa(a) - fb(b)), zip(exp_unit, mat_unit))
end

"Calculate the observe unitary matrix error, the error which we can obserwe
from quantum state"
function observe_unitary_error(exp_unit, mat_unit, index=1)
    # Remove phase
    exp_scale = cis(-angle(exp_unit[index]))
    mat_scale = cis(-angle(mat_unit[index]))

    return unitary_error(exp_unit, mat_unit, exp_scale, mat_scale)
end

"Calculate the minimum observe unitary matrix error."
function min_observe_unitary_error(exp_unit, mat_unit)
    l = length(exp_unit)
    return minimum(i -> observe_unitary_error(exp_unit, mat_unit, i), 1:l)
end

"Calculate the matrix norm"
function matrix_norm(mat)
    return abs(tr(adjoint(mat) * mat))
end

"Return the identity matrix."
eye(n) = Matrix(I, n, n)

end  # module Math


text:
facts("Evaluator") do

  # Set up a small example problem
  f(x) = sum(x.^2)
  P = BlackBoxOptim.fixeddim_problem(f; dims = 2)

  context("Basic evaluation with a single-objective function to be minimized") do

    e = ProblemEvaluator(P)

    @fact BlackBoxOptim.num_evals(e) => 0

    f = BlackBoxOptim.evaluate(e, [0.0, 1.0])
    @fact BlackBoxOptim.num_evals(e) => 1
    @fact BlackBoxOptim.last_fitness(e) => f

    f2 = BlackBoxOptim.evaluate(e, [2.0, 1.0])
    @fact BlackBoxOptim.num_evals(e) => 2
    @fact BlackBoxOptim.last_fitness(e) => f2

    @fact BlackBoxOptim.numdims(e) => 2

    @fact BlackBoxOptim.is_better(e, 0.0, 1.0) => true
    @fact BlackBoxOptim.is_better(e, 0.0, -1.0) => false

    a, b = [0, 0.5], [1, 1]

    @fact BlackBoxOptim.is_better(e, a, 1.0) => true
    @fact BlackBoxOptim.last_fitness(e) => 0.25

    @fact BlackBoxOptim.is_better(e, b, 1.0) => false
    @fact BlackBoxOptim.last_fitness(e) => 2.0

    @fact BlackBoxOptim.best_of(e, a, b) => (a, 0.25)
    @fact BlackBoxOptim.best_of(e, b, a) => (a, 0.25)

  end

end


text:
# 0.7 deprecations

import Base.colon
import Base.range

# deprecate remaining vectorized methods from Dates
@deprecate(
    DateTime(Y::AbstractArray{<:AbstractString}, f::AbstractString; locale::Locale=ENGLISH),
    DateTime.(Y, f; locale=locale) )
@deprecate(
    DateTime(Y::AbstractArray{<:AbstractString}, df::DateFormat=ISODateTimeFormat),
    DateTime.(Y, df) )
@deprecate(
    Date(Y::AbstractArray{<:AbstractString}, f::AbstractString; locale::Locale=ENGLISH),
    Date.(Y, f; locale=locale) )
@deprecate(
    Date(Y::AbstractArray{<:AbstractString}, df::DateFormat=ISODateFormat),
    Date.(Y, df) )
@deprecate(
    format(Y::AbstractArray{<:TimeType}, f::AbstractString; locale::Locale=ENGLISH),
    format.(Y, f; locale=locale) )
@deprecate(
    format(Y::AbstractArray{T}, df::DateFormat=default_format(T)) where {T<:TimeType},
    format.(Y, df) )

@deprecate +(a::GeneralPeriod, b::StridedArray{<:GeneralPeriod}) broadcast(+, a, b)
@deprecate +(a::StridedArray{<:GeneralPeriod}, b::GeneralPeriod) broadcast(+, a, b)
@deprecate -(a::GeneralPeriod, b::StridedArray{<:GeneralPeriod}) broadcast(-, a, b)
@deprecate -(a::StridedArray{<:GeneralPeriod}, b::GeneralPeriod) broadcast(-, a, b)

# #24258
# Physical units define an equivalence class: there is no such thing as a step of "1" (is
# it one day or one second or one nanosecond?). So require the user to specify the step
# (in physical units).
@deprecate colon(start::T, stop::T) where {T<:DateTime}   start:Day(1):stop
@deprecate colon(start::T, stop::T) where {T<:Date}       start:Day(1):stop
@deprecate colon(start::T, stop::T) where {T<:Time}       start:Second(1):stop

@deprecate range(start::DateTime, len::Integer)  range(start, Day(1), len)
@deprecate range(start::Date, len::Integer)      range(start, Day(1), len)

max_stars_repo_name: ihnorton/julia

max_issues_repo_name: ihnorton/julia

max_forks_repo_name: ihnorton/julia


text:
"""
`profile` is truncated to eliminate end-effects. This is equivalent to zero-padding described in [DIN EN ISO 16610-28:2017]().
"""
function zeropad(profile::SurfaceProfile, component::Symbol)
    c = profile.component
    bc = profile.λc * profile.resolution
    cut = round(Int, bc/2)
    c = c[cut:end-cut]
end


text:
using Dates, DataFrames, TimeSeries, GetColumn
using Test

colA = [4,3,2,1]
colB = [7,8,9,10]

colsAB = hcat(colA, colB)

colsABC = reshape(colsAB, (2,2,2))

df = DataFrame(ColA = colA, ColB = colB)

dates = collect((Date(Dates.now())-Day(3)):Day(1):Date(Dates.now()))
ta = TimeArray(dates, colsAB, ["ColA","ColB"])

@test getcolumn(colA, 1)    == colA
@test getcolumn(colB, 1)    == colB
@test getcolumn(colsAB, 1)  == colA
@test getcolumn(colsAB, 2)  == colB
@test getcolumn(colsABC, (1,2,0)) == colsABC[1,2,:]
@test getcolumn(colsABC, (2,0,1)) == colsABC[2,:,1]
@test getcolumn(df, "ColA") == colA
@test getcolumn(df, "ColB") == colB
@test getcolumn(df, :ColA)  == colA
@test getcolumn(df, :ColB)  == colB
@test getcolumn(ta, "ColA") == colA
@test getcolumn(ta, "ColB") == colB
@test getcolumn(ta, :ColA)  == colA
@test getcolumn(ta, :ColB)  == colB


text:
module Reader

using GZip
using BinDeps

export
    readimages,
    readlabels

# Constants

const IMAGEOFFSET = 16
const LABELOFFSET = 8

# Includes

include("readheader.jl")
include("readimages.jl")
include("readlabels.jl")

end


text:
module amis

using Formatting: format
using ..gmm: GMM, gmm_fit, importance_sample, update, rmcomponents, add_noise_comp
using AxUtil.MCDiagnostic: is_eff_ss
using AxUtil: dropdim1, dropdim2
using Distributions: logpdf, ncomponents, partype
import StatsBase: sample, weights
using StatsFuns: logaddexp
using NNlib: softmax
using Logging: LogLevel, @logmsg

const DETAILED = LogLevel(-500)

@inline is_tilt(d::GMM, tilt_amt::AbstractFloat) = isapprox(tilt_amt, 1.0) ? d : update(d; sigmas = d.sigmas * tilt_amt)
@inline range_j(j::Signed, chunksize::Signed) = UnitRange((j-1)*chunksize + 1, j*chunksize)

# Changed from Adaptive *Mixture* to Adaptive *Multiple* Importance Sampling (Cornuet et al. '09)
function adamult(dGMM::GMM, log_f::Function; max_iter=10, nsmps=1000,
                    IS_tilt=1.0, mle_max_iter=30, final_smps=1000)

    T = partype(dGMM)
    d = size(dGMM)
    n_total = max_iter*nsmps
    iter_fin = max_iter       # for logging: assume complete all iters (unless terminate early).
    X, U, L = Matrix{T}(undef, d, n_total), Vector{T}(undef, n_total), Vector{T}(undef, n_total)

    Qs = Vector{GMM}(undef, max_iter)
    cQ = is_tilt(dGMM, IS_tilt)  # exponential tilt if reqd

    for j = 1:max_iter
        Qs[j] = cQ

        cX = rand(cQ, nsmps)
        U[range_j(j, nsmps)] = log_f(cX)

        # Calculate log density of CURRENT SAMPLE under MIXTURE OF ALL IS DENSITIES
        cL = logpdf(Qs[1], cX)
        for j_prev = 2:j
            cL = logaddexp.(cL, logpdf(Qs[j_prev], cX))   # log(exp(Lprev) + exp(Lcur))
        end

        # Increment log density of ALL PREV SAMPLES under CURRENT MIXTURE DENSITY
        for j_prev = 1:j-1
            L_jprev = logpdf(cQ, X[:, range_j(j_prev, nsmps)])
            L[range_j(j_prev, nsmps)] = logaddexp.(L[range_j(j_prev, nsmps)], L_jprev)
        end

        # push samples / log proposal density to "stack"
        L[range_j(j, nsmps)] = cL
        X[:, range_j(j, nsmps)] = cX

        # current weights of ALL MC samples so far...
        W = softmax(U[1:j*nsmps] - L[1:j*nsmps])
        Wcur = softmax(U[range_j(j, nsmps)] - L[range_j(j, nsmps)])
        ess = is_eff_ss(Wcur) # .. but only interested in ESS of new particles
        # bson(format("{:s}_{:d}.bson", save_nm, j), dGMM=cQ, S=X[:, 1:j*nsmps], W=W)

        if ess / nsmps > 0.8
            iter_fin = j
            break
        end

        @logmsg DETAILED format("          | +++ AMIS ({:3d}/{:3d}) NOT TERMINATED, " *
                      "ncls={:4d}, ESS={:.2f}", j, max_iter, ncomponents(cQ), ess)

        # => When fitting the new density, start from the *original* density
        # so we don't lose any components due to _surprises_ in some iterations.
        gmm_iter = max(mle_max_iter - 2*j,5)
        Qbase = add_noise_comp(cQ; n = ncomponents(Qs[1]) - ncomponents(cQ), pi_eat=0.1)
        cQ = gmm_fit(X[:, 1:j*nsmps], W, Qbase; max_iter=gmm_iter, tol=1e-3, rm_inactive=true)
        cQ = is_tilt(cQ, IS_tilt)  # exponential tilt if reqd
    end

    S, logW = importance_sample(cQ, final_smps, log_f, shuffle=false)
    W = softmax(logW)
    ess = is_eff_ss(W)
    @logmsg DETAILED format("          | +++ AMIS ({:3d}/{:3d}) TERMINATED, " *
                  "ncls={:4d}, ESS={:.2f}", iter_fin, max_iter, ncomponents(cQ), ess)

    return cQ, S, W, ess, iter_fin
end

function adamult(S::Matrix{T}, W::Vector{T}, dGMM::GMM, log_f::Function; max_iter=10,
    nsmps=1000, IS_tilt=1.0, mle_max_iter=30, final_smps=1000) where T <: AbstractFloat

    cQ = gmm_fit(S, W, dGMM; max_iter=mle_max_iter, tol=1e-3, rm_inactive=true)
    adamult(cQ, log_f; max_iter=max_iter, nsmps=nsmps,
                        IS_tilt=IS_tilt, mle_max_iter=mle_max_iter, final_smps=final_smps)
end

end


text:
using MultivariateStats
using Base.Test

if !isdefined(Base, :normalize)
    normalize(x) = x/norm(x)
end

#Test equivalence of eigenvectors/singular vectors taking into account possible
#phase (sign) differences
#This may end up in Base.Test; see JuliaLang/julia#10651
function test_approx_eq_vecs{S<:Real,T<:Real}(a::StridedVecOrMat{S},
    b::StridedVecOrMat{T}, error=nothing)
    n = size(a, 2)
    @test n==size(b, 2) && size(a, 1)==size(b, 1)
    error==nothing && (error=n^3*(eps(S)+eps(T)))
    for i=1:n
        ev1, ev2 = a[:, i], b[:, i]
        deviation = min(abs(norm(ev1-ev2)), abs(norm(ev1+ev2)))
        if !isnan(deviation)
            @test_approx_eq_eps deviation 0.0 error
        end
    end
end

srand(34568)

## prepare data

d = 5
ns = [10, 15, 20]
nc = length(ns)
n = sum(ns)
Xs = Matrix{Float64}[]
ys = Vector{Int}[]
Ss = Matrix{Float64}[]
cmeans = zeros(d, nc)

for k = 1:nc
	R = qr(randn(d, d))[1]
	nk = ns[k]

	Xk = R * Diagonal(2 * rand(d) + 0.5) * randn(d, nk) .+ randn(d)
	yk = fill(k, nk)
	uk = vec(mean(Xk, 2))
	Zk = Xk .- uk
	Sk = Zk * Zk'

	push!(Xs, Xk)
	push!(ys, yk)
	push!(Ss, Sk)
	cmeans[:,k] = uk
end

X = hcat(Xs...)
y = vcat(ys...)
mv = vec(mean(X, 2))
Sw = zeros(d, d)
for k = 1:nc
	Sw += Ss[k]
end

Sb = zeros(d, d)
for k = 1:nc
	dv = cmeans[:,k] - mv
	Sb += ns[k] * (dv * dv')
end

Z = X .- mv
Sa = Z * Z'
@test_approx_eq Sw + Sb Sa

@assert size(X) == (d, n)
@assert size(y) == (n,)

## Stats

S = multiclass_lda_stats(nc, X, y)

@test S.dim == d
@test S.nclasses == nc

@test classweights(S) == ns
@test_approx_eq classmeans(S) cmeans
@test_approx_eq mean(S) mv

@test_approx_eq withclass_scatter(S) Sw
@test_approx_eq betweenclass_scatter(S) Sb

## Solve

emax = maximum(eigvals(Sw))
Sw_c = copy(Sw)
Sb_c = copy(Sb)

lambda = 1.0e-3
Sw_r = Sw + (lambda * emax) * eye(d)

P1 = mclda_solve(Sb, Sw, :gevd, nc-1, lambda)
@test_approx_eq P1' * Sw_r * P1 eye(nc-1)
U = Sb * P1
V = Sw_r * P1
# test whether U is proportional to V,
# which indicates P is the generalized eigenvectors
@test_approx_eq U V*Diagonal(vec(mean(U./V, 1)))

P2 = mclda_solve(Sb, Sw, :whiten, nc-1, lambda)
@test_approx_eq P2' * Sw_r * P2 eye(nc-1)

test_approx_eq_vecs(P1, P2)

## LDA

M = fit(MulticlassLDA, nc, X, y; method=:gevd, regcoef=lambda)
@test indim(M) == d
@test outdim(M) == nc - 1
@test_approx_eq projection(M) P1
@test_approx_eq M.pmeans M.proj'cmeans
@test_approx_eq transform(M, X) M.proj'X

M = fit(MulticlassLDA, nc, X, y; method=:whiten, regcoef=lambda)
@test indim(M) == d
@test outdim(M) == nc - 1
# @test_approx_eq projection(M) P2  # signs may change
@test_approx_eq M.pmeans M.proj'cmeans
@test_approx_eq transform(M, X) M.proj'X

## High-dimensional LDA (subspace LDA)

# Low-dimensional case (no singularities)

dX = randn(5,1500);
# 3 groups of 500 with zero mean
for i = 0:500:1000
    dX[:,(1:500)+i] = dX[:,(1:500)+i] .- mean(dX[:,(1:500)+i],2)
end
centers = [zeros(5) [10.0;zeros(4)] [0.0;10.0;zeros(3)]]
X = [dX[:,1:500].+centers[:,1] dX[:,501:1000].+centers[:,2] dX[:,1001:1500].+centers[:,3]]
label = [fill(1,500); fill(2,500); fill(3,500)]
M = fit(SubspaceLDA, X, label)
@test indim(M) == 5
@test outdim(M) == 2
@test_approx_eq mean(M) vec(mean(centers,2))
@test_approx_eq classmeans(M) centers
@test classweights(M) == [500,500,500]
x = rand(5)
@test_approx_eq transform(M, x) projection(M)'*x
dcenters = centers .- mean(centers,2)
Hb = dcenters*sqrt(500)
Sb = Hb*Hb'
Sw = dX*dX'
# When there are no singularities, we need to solve the "regular" LDA eigenvalue equation
proj = projection(M)
@test_approx_eq Sb*proj Sw*proj*Diagonal(M.λ)
@test_approx_eq proj'*Sw*proj eye(2,2)

# also check that this is consistent with the conventional algorithm
Mld = fit(MulticlassLDA, 3, X, label)
for i = 1:2
    @test_approx_eq abs(dot(normalize(proj[:,i]), normalize(Mld.proj[:,i]))) 1
end

# High-dimensional case (undersampled => singularities)
X = randn(10^6, 9)
label = rand(1:3, 9)
M = fit(SubspaceLDA, X, label)
centers = M.cmeans
for i = 1:3
    flag = label.==i
    @test_approx_eq centers[:,i] mean(X[:,flag], 2)
    @test M.cweights[i] == sum(flag)
end
dcenters = centers .- mean(X, 2)
Hb = dcenters*Diagonal(sqrt(M.cweights))
Hw = X - centers[:,label]
@test_approx_eq (M.projw'*Hb)*(Hb'*M.projw)*M.projLDA (M.projw'*Hw)*(Hw'*M.projw)*M.projLDA*Diagonal(M.λ)
Gw = projection(M)'*Hw
@test_approx_eq Gw*Gw' eye(2,2)

# Test that nothing breaks if one class has no members
label[1:4] = 1
label[5:9] = 3  # no data points have label == 2
M = fit(SubspaceLDA, X, label)
centers = M.cmeans
dcenters = centers .- mean(X, 2)
Hb = dcenters*Diagonal(sqrt(M.cweights))
Hw = X - centers[:,label]
@test_approx_eq (M.projw'*Hb)*(Hb'*M.projw)*M.projLDA (M.projw'*Hw)*(Hw'*M.projw)*M.projLDA*Diagonal(M.λ)


text:
export Contingency

"""
Description:

    Contingency factor node

    The contingency distribution is a multivariate generalization of
    the categorical distribution. As a bivariate distribution, the
    contingency distribution defines the joint probability
    over two unit vectors. The parameter p encodes a contingency matrix
    that specifies the probability of co-occurrence.

    out1 ∈ {0, 1}^d1 where Σ_j out1_j = 1
    out2 ∈ {0, 1}^d2 where Σ_k out2_k = 1
    p ∈ [0, 1]^{d1 × d2}, where Σ_jk p_jk = 1

    f(out1, out2, p) = Con(out1, out2 | p)
                     = Π_jk p_jk^{out1_j * out2_k}

    A Contingency distribution over more than two variables requires
    higher-order tensors as parameters; these are not implemented in ForneyLab.

Interfaces:

    1. out1
    2. out2
    3. p

Construction:

    Contingency(id=:some_id)
"""
mutable struct Contingency <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function Contingency(out1, out2, p; id=generateId(Contingency))
        @ensureVariables(out1, out2, p)
        self = new(id, Array{Interface}(undef, 3), Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)
        self.i[:out1] = self.interfaces[1] = associate!(Interface(self), out1)
        self.i[:out2] = self.interfaces[2] = associate!(Interface(self), out2)
        self.i[:p] = self.interfaces[3] = associate!(Interface(self), p)

        return self
    end
end

slug(::Type{Contingency}) = "Con"

format(dist::ProbabilityDistribution{Multivariate, Contingency}) = "$(slug(Contingency))(p=$(format(dist.params[:p])))"

ProbabilityDistribution(::Type{Multivariate}, ::Type{Contingency}; p=1/9*ones(3,3)) = ProbabilityDistribution{Multivariate, Contingency}(Dict(:p=>p))
ProbabilityDistribution(::Type{Contingency}; p=1/9*ones(3,3)) = ProbabilityDistribution{Multivariate, Contingency}(Dict(:p=>p))

dims(dist::ProbabilityDistribution{Multivariate, Contingency}) = (length(size(dist.params[:p])),)

vague(::Type{Contingency}, n_factors::Tuple{Int64, Int64}=(3,3)) = ProbabilityDistribution(Multivariate, Contingency, p=(1/prod(n_factors))*ones(n_factors))

isProper(dist::ProbabilityDistribution{Multivariate, Contingency}) = (abs(sum(dist.params[:p])-1.) < 1e-6)

# Entropy functional
function differentialEntropy(dist::ProbabilityDistribution{Multivariate, Contingency})
    B = dist.params[:p]
    if isa(B, Vector{<:Matrix}) # Contingency tensor
        H = 0.0
        for k=1:length(B)
            H += -sum(B[k].*log.(clamp.(B[k], tiny, Inf)))
        end
    else # Contingency matrix
        H = -sum(B.*log.(clamp.(B, tiny, Inf)))
    end

    return H
end