text:
# Load galaxy bivariate normals with no PSF.  Used in DeterministicVIImagePSF.
function load_gal_bvn_mixtures{NumType <: Number}(
                    S::Int64,
                    patches::Matrix{SkyPatch},
                    source_params::Vector{Vector{NumType}},
                    active_sources::Vector{Int},
                    n::Int;
                    calculate_gradient::Bool=true,
                    calculate_hessian::Bool=true)
    # To maintain consistency with the rest of the code, use a 4d
    # array.  The first dimension was previously the PSF component.
    gal_mcs = Array(GalaxyCacheComponent{NumType}, 1, 8, 2, S)

    for s in 1:S
        sp  = source_params[s]
        world_loc = sp[lidx.u]
        p = patches[s, n]
        m_pos = linear_world_to_pix(
            p.wcs_jacobian, p.center, p.pixel_center, world_loc)

        for i = 1:2 # i indexes dev vs exp galaxy types.
            e_dev_dir = (i == 1) ? 1. : -1.
            e_dev_i = (i == 1) ? sp[lidx.e_dev] : 1. - sp[lidx.e_dev]

            # Galaxies of type 1 have 8 components, and type 2 have 6 components.
            for j in 1:[8,6][i]
                gal_mcs[1, j, i, s] = GalaxyCacheComponent(
                    e_dev_dir, e_dev_i, galaxy_prototypes[i][j], m_pos,
                    sp[lidx.e_axis], sp[lidx.e_angle], sp[lidx.e_scale],
                    calculate_gradient && (s in active_sources),
                    calculate_hessian)
            end
        end
    end

    gal_mcs
end

# Get a GalaxyCacheComponent with no PSF.
import ..Model.GalaxyCacheComponent
function GalaxyCacheComponent{NumType <: Number}(
    e_dev_dir::Float64, e_dev_i::NumType,
    gc::GalaxyComponent, u::Vector{NumType},
    e_axis::NumType, e_angle::NumType, e_scale::NumType,
    calculate_gradient::Bool, calculate_hessian::Bool)

    # Declare in advance to save memory allocation.
    const empty_sig_sf =
        GalaxySigmaDerivs(Array(NumType, 0, 0), Array(NumType, 0, 0, 0))

    XiXi = get_bvn_cov(e_axis, e_angle, e_scale)
    var_s = gc.nuBar * XiXi

    # d siginv / dsigma is only necessary for the Hessian.
    bmc = BvnComponent{NumType}(
        SVector{2, NumType}(u), var_s, gc.etaBar,
        calculate_gradient && calculate_hessian)

    if calculate_gradient
        sig_sf = GalaxySigmaDerivs(
            e_angle, e_axis, e_scale, XiXi, calculate_hessian)
        sig_sf.j .*= gc.nuBar
        if calculate_hessian
            # The tensor is only needed for the Hessian.
            sig_sf.t .*= gc.nuBar
        end
    else
        sig_sf = empty_sig_sf
    end

    GalaxyCacheComponent(e_dev_dir, e_dev_i, bmc, sig_sf)
end

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

const dir = "/home/rgiordan/Documents/git_repos/CelesteDev.jl/"

"""
Convolve a populated set of SensitiveFloat matrices in fsms with the PSF
and store them in the matching fs*m_conv SensitiveFloat matrices.
"""
function convolve_fs1m_image!(fsms::FSMSensitiveFloatMatrices, s::Int)
    for h in 1:size(fsms.fs1m_image, 1), w in 1:size(fsms.fs1m_image, 2)
        fsms.fs1m_image_padded[h, w] = fsms.fs1m_image[h, w]
    end

    convolve_sensitive_float_matrix!(
        fsms.fs1m_image_padded, fsms.psf_fft, fsms.fs1m_conv_padded)

    for h in 1:size(fsms.fs1m_image, 1), w in 1:size(fsms.fs1m_image, 2)
        fsms.fs1m_conv[h, w] =
            fsms.fs1m_conv_padded[fsms.pad_pix_h + h, fsms.pad_pix_w + w]
    end

    # Set return type
    return true
end

function clear_fs1m!(fsms::FSMSensitiveFloatMatrices)
    for sf in fsms.fs1m_image clear!(sf) end
    for sf in fsms.fs1m_image_padded clear!(sf) end
    for sf in fsms.fs1m_conv clear!(sf) end
    for sf in fsms.fs1m_conv_padded clear!(sf) end
end

function clear_brightness!(fsms::FSMSensitiveFloatMatrices)
    for sf in fsms.E_G clear!(sf) end
    for sf in fsms.var_G clear!(sf) end
end

"""
Populate the fs1m shape matrices and convolve with the PSF for a given
source and band.  Assumes that fsms.psf_fft has already been set.

The result is the sources shapes, convolved with the PSF, stored in fsms.fs1m_conv.

TODO: pass in derivative flags
"""
function populate_gal_fsm_image!(
            ea::ElboArgs{Float64},
            s::Int,
            n::Int,
            gal_mcs::Array{GalaxyCacheComponent{Float64}, 4},
            fsms::FSMSensitiveFloatMatrices)
    clear_fs1m!(fsms)
    is_active_source = s in ea.active_sources
    p = ea.patches[s, n]
    H_patch, W_patch = size(p.active_pixel_bitmap)
    for w_patch in 1:W_patch, h_patch in 1:H_patch
        if !p.active_pixel_bitmap[h_patch, w_patch]
            continue
        end
        h_image = h_patch + p.bitmap_offset[1]
        w_image = w_patch + p.bitmap_offset[2]

        h_fsm = h_image - fsms.h_lower + 1
        w_fsm = w_image - fsms.w_lower + 1

        x = SVector{2, Float64}([h_image, w_image])
        populate_gal_fsm!(ea.elbo_vars.bvn_derivs,
                          fsms.fs1m_image[h_fsm, w_fsm],
                          s, x, is_active_source, Inf,
                          p.wcs_jacobian,
                          gal_mcs)
    end
    convolve_fs1m_image!(fsms, s)
end

"""
Populate the fs1m shape matrices and convolve with the PSF for a given
source and band.  Assumes that fsms.psf_fft has already been set.

The result is the sources shapes, convolved with the PSF, stored in fsms.fs1m_conv.

TODO: pass in derivative flags
"""
function populate_star_fsm_image!(
            ea::ElboArgs{Float64},
            s::Int,
            n::Int,
            psf_image::Matrix{Float64},
            fs0m_conv::fs0mMatrix,
            h_lower::Int, w_lower::Int,
            kernel_fun::Any, # TODO: Is there a type that can be used here?
            kernel_width::Int)
    for sf in fs0m_conv clear!(sf) end
    # The pixel location of the star.
    star_loc_pix =
        linear_world_to_pix(ea.patches[s, n].wcs_jacobian,
                            ea.patches[s, n].center,
                            ea.patches[s, n].pixel_center,
                            ea.vp[s][lidx.u]) -
        Float64[ h_lower - 1, w_lower - 1]
    interpolate!(kernel_fun, kernel_width,
                 fs0m_conv, psf_image, star_loc_pix,
                 ea.patches[s, n].wcs_jacobian,
                 ea.elbo_vars.elbo.has_gradient,
                 ea.elbo_vars.elbo.has_hessian)
end

# TODO: rename this to source image brighntess or something.
"""
Updates fsms.E_G and fsms.var_G in place with the contributions from this
source in this band.
"""
function accumulate_source_image_brightness!(
    ea::ElboArgs{Float64},
    s::Int,
    n::Int,
    fsms::FSMSensitiveFloatMatrices,
    sb::SourceBrightness{Float64})

    is_active_source = s in ea.active_sources
    calculate_hessian =
        ea.elbo_vars.elbo.has_hessian && ea.elbo_vars.elbo.has_gradient &&
        is_active_source

    image_fft = [ sf.v[] for sf in fsms.fs1m_conv ]
    p = ea.patches[s, n]
    H_patch, W_patch = size(p.active_pixel_bitmap)
    for w_patch in 1:W_patch, h_patch in 1:H_patch
        h_fsm = h_patch + p.bitmap_offset[1] - fsms.h_lower + 1
        w_fsm = w_patch + p.bitmap_offset[2] - fsms.w_lower + 1
        if !p.active_pixel_bitmap[h_patch, w_patch]
            continue
        end
        accumulate_source_pixel_brightness!(
                            ea.elbo_vars,
                            ea,
                            fsms.E_G[h_fsm, w_fsm],
                            fsms.var_G[h_fsm, w_fsm],
                            fsms.fs0m_conv[h_fsm, w_fsm],
                            fsms.fs1m_conv[h_fsm, w_fsm],
                            sb, ea.images[n].b, s, is_active_source)
    end
end

"""
Uses the values in fsms to add the contribution from this band to the ELBO.
"""
function accumulate_band_in_elbo!(
    ea::ElboArgs{Float64},
    fsm_mat::Matrix{FSMSensitiveFloatMatrices},
    sbs::Vector{SourceBrightness{Float64}},
    gal_mcs::Array{GalaxyCacheComponent{Float64}, 4},
    n::Int)

    for s in 1:ea.S
        # some sources don't appear in some images
        sum(ea.patches[s, n].active_pixel_bitmap) > 0 || continue

        fsms = fsm_mat[s, n]
        clear_brightness!(fsms)
        populate_star_fsm_image!(
            ea, s, n, fsms.psf, fsms.fs0m_conv,
            fsms.h_lower, fsms.w_lower,
            fsms.kernel_fun, fsms.kernel_width)
        populate_gal_fsm_image!(ea, s, n, gal_mcs, fsms)
        accumulate_source_image_brightness!(ea, s, n, fsms, sbs[s])
    end

    # Iterate only over active sources, since we have already added the
    # contributions from non-active sources to E_G and var_G.
    for s in ea.active_sources
        p = ea.patches[s, n]
        fsms = fsm_mat[s, n]
        H_patch, W_patch = size(p.active_pixel_bitmap)
        for w_patch in 1:W_patch, h_patch in 1:H_patch
            if !p.active_pixel_bitmap[h_patch, w_patch]
                continue
            end
            h_image = h_patch + p.bitmap_offset[1]
            w_image = w_patch + p.bitmap_offset[2]

            image = ea.images[n]
            this_pixel = image.pixels[h_image, w_image]

            if Base.isnan(this_pixel)
                continue
            end

            # These are indices within the fs?m image.
            h_fsm = h_image - fsms.h_lower + 1
            w_fsm = w_image - fsms.w_lower + 1

            E_G = fsms.E_G[h_fsm, w_fsm]
            var_G = fsms.var_G[h_fsm, w_fsm]

            # There are no derivatives with respect to epsilon, so can
            # afely add to the value.
            E_G.v[] += image.epsilon_mat[h_image, w_image]

            if E_G.v[] < 0
                # warn("Image ", n, " sources ", s, " pixel ", (h_image, w_image),
                #      " has negative brightness ", E_G.v[])
                continue
            end

            # Note that with a kernel_width > 1 negative values are
            # possible, and this will result in an error in
            # add_elbo_log_term.

            # Add the terms to the elbo given the brightness.
            iota = image.iota_vec[h_image]
            add_elbo_log_term!(
                ea.elbo_vars, E_G, var_G, ea.elbo_vars.elbo, this_pixel, iota)
            add_scaled_sfs!(ea.elbo_vars.elbo, E_G, -iota)

            # Subtract the log factorial term. This is not a function of the
            # parameters so the derivatives don't need to be updated. Note
            # that even though this does not affect the ELBO's maximum,
            # it affects the optimization convergence criterion, so I will
            # leave it in for now.
            ea.elbo_vars.elbo.v[] -= lfact(this_pixel)
        end
    end
end

function elbo_likelihood_with_fft!(
    ea::ElboArgs,
    fsm_mat::Matrix{FSMSensitiveFloatMatrices})

    sbs = load_source_brightnesses(ea,
        calculate_gradient=ea.elbo_vars.elbo.has_gradient,
        calculate_hessian=ea.elbo_vars.elbo.has_hessian)

    clear!(ea.elbo_vars.elbo)
    for n in 1:ea.N
        gal_mcs = load_gal_bvn_mixtures(
                ea.S, ea.patches, ea.vp, ea.active_sources, n,
                calculate_gradient=ea.elbo_vars.elbo.has_gradient,
                calculate_hessian=ea.elbo_vars.elbo.has_hessian)
        accumulate_band_in_elbo!(ea, fsm_mat, sbs, gal_mcs, n)
    end
end

function initialize_fft_elbo_parameters(
    images::Vector{Image},
    vp::VariationalParams{Float64},
    patches::Matrix{SkyPatch},
    active_sources::Vector{Int};
    use_raw_psf=true,
    use_trimmed_psf=true)

    ea = ElboArgs(images, vp, patches, active_sources, psf_K=1)
    load_active_pixels!(images, ea.patches; exclude_nan=false)
    if use_raw_psf
        psf_image_mat = Array{Matrix{Float64}}(ea.S, ea.N)
        for n in 1:ea.N, s in 1:ea.S
            img = images[n]
            world_loc = ea.vp[s][lidx.u]
            pixel_loc = WCS.world_to_pix(img.wcs, world_loc)
            psf_image_mat[s, n] =
                eval_psf(img.raw_psf_comp, pixel_loc[1], pixel_loc[2])
        end
    else
        psf_image_mat = Matrix{Float64}[
            get_psf_at_point(ea.patches[s, b].psf) for s in 1:ea.S, b in 1:ea.N]
    end

    if use_trimmed_psf
        for n in 1:ea.N, s in 1:ea.S
            psf_image_mat[s, n] = trim_psf(psf_image_mat[s, n])
        end
    end

    fsm_mat = FSMSensitiveFloatMatrices[
        FSMSensitiveFloatMatrices() for s in 1:ea.S, n in 1:ea.N]
    initialize_fsm_sf_matrices!(fsm_mat, ea, psf_image_mat)
    ea, fsm_mat
end

@doc """
Return a function callback for an FFT elbo.
"""
function get_fft_elbo_function{T}(
    ea::ElboArgs{T}, fsm_mat::Matrix{FSMSensitiveFloatMatrices})
    function elbo_fft_opt(ea::ElboArgs)
        @assert ea.psf_K == 1
        elbo = ea.elbo_vars.elbo
        kl_source = SensitiveFloat{T}(length(CanonicalParams), 1,
                                      elbo.has_gradient, elbo.has_hessian)
        elbo_likelihood_with_fft!(ea, fsm_mat)
        subtract_kl_all_sources!(ea, elbo, kl_source,
                                 KL_HELPER_POOL[Base.Threads.threadid()])
        return deepcopy(elbo)
    end
end


text:
# # Multiple dispatch

#-

# **Multiple dispatch** is a key feature of Julia, that we will explore in this notebook.
#
# It helps make software fast. It also makes software extensible, programmable, and downright fun to play with.
#
# It may just herald a breakthrough for parallel computation.

#-

# 1. Roman numerals
# 2. Functions
# 3. Parallel computing

#-

# ## 1. Roman numerals (for fun)

#-

# Let's define a **new struct** that represents a Roman numeral. For coding simplicity, we'll just deal with numbers between 0 and 9.
#
# **Exercise**: Extend this to larger numbers. (Recall that Roman numbers are a base-10 system!)

struct Roman
    n::Int
end

Base.show(io::IO, r::Roman) = print(io, 'ⅰ' + (r.n - 1) % 10 )  # nice display; 'ⅰ' is a Unicode Roman numeral

# We can create an object of this type as follows:

Roman(4)

#-

typeof.([5 5.0 Roman(5) "Five" '5'  5//1])

# We would like to display it nicely, in Roman numerals:

x = [7 1 2 5 8 9]
Roman.(x)   # equivalent to map(Roman, x)  or  [Roman(w) for w in x]

# It'd be nice to be able to add Roman numerals together like normal numbers:

Roman(4) + Roman(5)

# But Julia doesn't know how to do that. Let's teach it by `import`ing the `+` function, which then allows us to _extend_ its definition:

import Base: +, *

+(a::Roman, b::Roman) = Roman(a.n + b.n)

#-

Roman(4) + Roman(5)

# This **adds a new method** to the function `+`:

methods(+)

#-

import Base.*
*(i::Roman, j::Roman) = Roman(i.n * j.n)                     # Multiply like a Roman

#-

Roman(3) * Roman(2)

#-

Roman.(1:3) .* [Roman(1) Roman(2) Roman(3)]

# But

Roman(3) * 2

#-

## Complicated mytimes to decide what to do based on type
## not suggested, better way coming soon
function mytimes(i,j)
  if isa(i,Roman) & isa(j,Number)
        return  fill(1, i.n, j)   # i by j matrix with ones
    elseif    isa(i,Number) & isa(j,Roman)
        return "😄"^ (i*j.n)   #  i * j happy faces
    else
        return("I Don't know")
    end
end

#-

mytimes(4,Roman(3)) # Twelve happys

#-

mytimes(Roman(4),3) # 4x3 matrix with ones

# The simplest thing to do is to explicitly define multiplication of a `Roman` by a number. We can do it as we see fit:

*(i::Number, j::Roman) = "😄"^ (i*j.n)        #  i * j happy faces

*(i::Roman, j::Number) =   fill(1, i.n, j)       # i by j matrix

#-

3 * Roman(3) # Nine happys

#-

Roman(3) * 5  # Three by Five matrix of ones

#-

t(x::Roman,y::Roman) = x.n * y.n

#-

t(Roman(5),Roman(4))

#-

##  Notice how tight the assembler is!
@code_native t(Roman(2),Roman(4))

# ## Functions

import Base: *, +, ^

#-

*(α::Number,   g::Function) = x -> α * g(x)   # Scalar times function

*(f::Function, λ::Number)   = x -> f(λ * x)   # Scale the argument

*(f::Function, g::Function) = x -> f(g(x))    # Function composition  -- abuse of notation!  use \circ in Julia 0.6

^(f::Function, n::Integer) = n == 1 ? f : f*f^(n-1) # A naive exponentiation algorithm by recursive multiplication

#-

+(f::Function, g::Function) = x -> f(x) + g(x)

# For example, the exponential function is defined as
#
# $$\exp(x) = \sum_{n=0}^\infty \frac{1}{n!} x^n.$$
#
# We can think of this just in terms of functions:
#
# $$\exp = \sum_{n=0}^\infty \frac{1}{n!} \mathrm{pow}_n,$$
#
# where $\mathrm{pow}_n(x) = x^n$.
#
# (starts to blur the symbolic with the numerical!)

pow(n) = x -> x^n

myexp = sum(1/factorial(big(n)) * pow(n) for n in 0:100)   # taylor series not efficient!

#-

[myexp(1); exp(1); exp(big(1))]

#-

f = x -> x^2
f(10)

#-

g = 3f
g(10)

#-

(f^2)(10)  # since we defined multiplication of functions as composition

#-

using Plots;
gr()

#-

x = pi*(0:0.001:4)

plot(x, sin.(x),    c="black", label="Fun")
plot!(x, (12*sin).(x),    c="green", label="Num * Fun")
plot!(x, (sin*12).(x),    c="red", alpha=0.9, label="Fun * Num")
plot!(x, (5*sin*exp).(x), c="blue", alpha=0.2, label="Num * Fun * Fun")

#-

plot([12*sin, sin*12, 5*sin*exp], 0:.01:4π, α=[1 .9 .2], c=[:green :red :blue])

# <img src="https://lh4.googleusercontent.com/--z5eKJbB7sg/UffjL1iAd4I/AAAAAAAABOc/S_wDVyDOBfQ/gauss.jpg">

#-

# ###  "Sin^2 phi is odious to me, even though Laplace made use of it; should  it be feared that sin^2 phi might become ambiguous, which would perhaps  never occur, or at most very rarely when speaking of sin(phi^2), well  then, let us write (sin phi)^2, but not sin^2 phi, which by analogy  should signify sin(sin phi)." -- Gauss

x=(0:.01:2) * pi;

plot(x, (sin^2).(x), c="blue")     # Squaring just works, y=sin(sin(x)), Gauss would be pleased!
plot!(x, sin.(x).^2,  c="red")

# # Exercise

h(a, b::Any) = "fallback"
h(a::Number, b::Number) = "a and b are both numbers"
h(a::Number, b) = "a is a number"
h(a, b::Number) = "b is a number"
h(a::Integer, b::Integer) = "a and b are both integers"

#-

## Try playing with h

ext: jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# # Gaussian Process Regression (for large data)
#
# ### Loading necessary packages
using AugmentedGaussianProcesses
using Distributions
using Plots

# We create a toy dataset with X ∈ [-20, 20] and y = 5 * sinc(X)
N = 1000
X = reshape((sort(rand(N)) .- 0.5) * 40.0, N, 1)
σ = 0.01

function latent(x)
    return 5.0 * sinc.(x)
end
Y = vec(latent(X) + σ * randn(N));

# Visualization of the data :
scatter(X, Y; lab="")

# ## Gaussian noise

# In this first example we are going to look at the effect of using
# inducing points compared to the true Gaussian Process
# For simplicity we will keep all inducing points and kernel parameters fixed

Ms = [4, 8, 16, 32, 64];
# Create an empty array of GPs
models = Vector{AbstractGPModel}(undef, length(Ms) + 1);
# Chose a kernel
kernel = SqExponentialKernel();#  + PeriodicKernel()
# And Run sparse classification with an increasing number of inducing points
for (index, num_inducing) in enumerate(Ms)
    @info "Training with $(num_inducing) points"
    m = SVGP(
        kernel, # Kernel
        GaussianLikelihood(σ), # Likelihood used
        AnalyticVI(), # Inference usede to solve the problem
        inducingpoints(KmeansAlg(num_inducing), X); # Inducing points initialized with kmeans
        optimiser=false, # Keep kernel parameters fixed
        Zoptimiser=false, # Keep inducing points locations fixed
    )
    @time train!(m, X, Y, 100) # Train the model for 100 iterations
    models[index] = m # Save the model in the array
end

# Train the model without any inducing points (no approximation)
@info "Training with full model"
mfull = GP(
    X,
    Y,
    kernel;
    noise=σ,
    opt_noise=false, # Keep the noise value fixed
    optimiser=false, # Keep kernel parameters fixed
)
@time train!(mfull, 5);
models[end] = mfull;

# Create a grid and compute prediction on it
function compute_grid(model, n_grid=50)
    mins = -20
    maxs = 20
    x_grid = range(mins, maxs; length=n_grid) # Create a grid
    y_grid, sig_y_grid = proba_y(model, reshape(x_grid, :, 1)) # Predict the mean and variance on the grid
    return y_grid, sig_y_grid, x_grid
end;

# Plot the data as a scatter plot
function plotdata(X, Y)
    return Plots.scatter(X, Y; alpha=0.33, msw=0.0, lab="", size=(300, 500))
end

function plot_model(model, X, Y, title=nothing)
    n_grid = 100
    y_grid, sig_y_grid, x_grid = compute_grid(model, n_grid)
    title = if isnothing(title)
        (model isa SVGP ? "M = $(dim(model[1]))" : "full")
    else
        title
    end

    p = plotdata(X, Y)
    Plots.plot!(
        p,
        x_grid,
        y_grid;
        ribbon=2 * sqrt.(sig_y_grid), # Plot 2 std deviations
        title=title,
        color="red",
        lab="",
        linewidth=3.0,
    )
    if model isa SVGP # Plot the inducing points as well
        Plots.plot!(
            p,
            vec(model.f[1].Z),
            zeros(dim(model.f[1]));
            msize=2.0,
            color="black",
            t=:scatter,
            lab="",
        )
    end
    return p
end;

Plots.plot(
    plot_model.(models, Ref(X), Ref(Y))...; layout=(1, length(models)), size=(1000, 200)
) # Plot all models and combine the plots

# ## Non-Gaussian Likelihoods
# We now look at using another noise than Gaussian noise.
# In AGP.jl you can use the Student-T likelihood,
# the Laplace likelihood and the Heteroscedastic likelihood

# We will use the same toy dataset for our experiment

# Create an array of model with different likelihoods:
likelihoods = [
    StudentTLikelihood(3.0), LaplaceLikelihood(3.0), HeteroscedasticLikelihood(1.0)
]
ngmodels = Vector{AbstractGPModel}(undef, length(likelihoods) + 1)
for (i, l) in enumerate(likelihoods)
    @info "Training with the $(l)" # We need to use VGP
    m = VGP(
        X,
        Y, # First arguments are the input and output
        kernel, # Kernel
        l, # Likelihood used
        AnalyticVI(); # Inference usede to solve the problem
        optimiser=false, # Keep kernel parameters fixed
    )
    @time train!(m, 10) # Train the model for 100 iterations
    ngmodels[i] = m # Save the model in the array
end

ngmodels[end] = models[end] # Add the Gaussian model
# We can now repeat the prediction from before :
Plots.plot(
    plot_model.(
        ngmodels, Ref(X), Ref(Y), ["Student-T", "Laplace", "Heteroscedastic", "Gaussian"]
    )...;
    layout=(2, 2),
    size=(1000, 200),
) # Plot all models and combine the plots


text:
# This restarts workers, need this until they fix memory leaks.
# Unfortunately it doesn't work from within the module.
function resetworkers(ncpu)
    rmworkers()
    addprocs(ncpu)
    require("CUDArt")
    @everywhere CUDArt.device((myid()-1) % CUDArt.devcount())
    require("CUBLAS")
    require("Knet")
    require("KUparser")
end # resetworkers

# Really remove all workers
# There is also Base.terminate_all_workers() we can use here...
function rmworkers()
    nprocs() > 1 && rmprocs(workers())
    while nworkers() > 1
        sleep(1)
    end
end

function restartmachines()
    machines = String[]
    for i in keys(Base.map_pid_wrkr)
        i == 1 && continue
        w = Base.map_pid_wrkr[i]
        push!(machines, w.host)
    end
    Base.terminate_all_workers()
    addprocs(machines)
end

# This does not distribute gpus equally:
# @date @everywhere using CUDArt
# @date @everywhere CUDArt.device((myid()-1) % CUDArt.devcount())

function restartcuda()
    @everywhere require("CUDArt")
    d = Dict()
    for i in keys(Base.map_pid_wrkr)
        i == 1 && continue
        a = Base.map_pid_wrkr[i].bind_addr
        n = get!(d, a, 0)
        # @show (a,n)
        @fetchfrom i CUDArt.device(n % CUDArt.devcount())
        d[a] = n+1
    end
end

function restartmachines(host::String)
    wid = 0
    for i in keys(Base.map_pid_wrkr)
        i == 1 && continue
        w = Base.map_pid_wrkr[i]
        w.host == host && (wid = i; break)
    end
    if wid == 0
        warn("Cannot find worker on $host")
        return nothing
    end
    ret = rmprocs(wid; waitfor=0.5)
    if ret != :ok
        warn("Forcibly interrupting busy worker $i")
        # Might be computation bound, interrupt them and try again
        interrupt(wid)
        ret = rmprocs(wid; waitfor=0.5)
        if ret != :ok
            warn("Unable to terminate worker $i")
            return nothing
        end
    end
    @assert ret == :ok
    return addprocs([host])[1]
end

function restartmachines(machines::AbstractVector)
    map(restartmachines, machines)
end


text:
using Plots
using LaTeXStrings

# backend(:pgfplotsx)
backend(:pyplot)
# backend(:gr)

#reset past font scaling
Plots.scalefontsizes()
#scale fonts
Plots.scalefontsizes(5.00)
#Extracting stored varibles

#Add colors
lightblue = colorant"rgb(63%,74%,78%)"
orange = colorant"rgb(85%,55%,13%)"
silver = colorant"rgb(69%,67%,66%)"
rust = colorant"rgb(72%,26%,6%)"

f(x,y) = 2 * x * y - 4 * sum((1 - y).^2) + 100 * Float64(x < 0)
flocal(x,y) = f(0.0, 1.0) + 2 * (x - 0.0) + 2 * (x - 0.0) * (y - 1.0)

x1 = -0.5 : 0.05: 1.5
y1 = - 0.5 : 0.05 : 1.5

x2 = -0.5 : 0.10: 1.5
y2 = -0.5 : 0.10 : 1.5

srf = wireframe(x1, y1, f, zlims=(-3,5), camera=(15,35), linecolor=silver, linealpha=0.2, xlabel=L"x", ylabel=L"y", label=L"\mathrm{true} \ f)", thickness_scaling=1.5, linewidth=0.75)

wireframe!(srf, x2, y2, flocal, linecolor=rust, label=L"\mathrm{approximate} \ f)", legend=:topright)

savefig(srf, "./out/figures/empty_threats.pdf")
run(`pdfcrop --margins '5 5 5 5' ./out/figures/empty_threats.pdf ./out/figures/empty_threats.pdf`)

max_stars_repo_name: f-t-s/CMD

max_issues_repo_name: f-t-s/CMD

max_forks_repo_name: f-t-s/CMD


text:
using Base.Iterators: product
using Distributed
try
    using ProgressMeter
catch
    # pass
end

@everywhere using AstroBase

μ = 3.986004418e14

procs = nprocs()
n = 4
abort = false

if !isempty(ARGS)
    try
        global n
        n = something(tryparse.(Int, ARGS)...)
    catch
        # pass
    end

    abort = "abort" in ARGS
else
    if procs > 1
        n = (procs - 1) * 2
    end
end

@everywhere acceptor(ele, μ) = false

# Reject element sets with an out of range anomaly.
@everywhere function hyperbolic_rejector(ele, μ)
    abs(ele[end]) >= acos(-1/ele[2]) && return true

    false
end

@everywhere function keplerian_roundtrip(ele₀, μ, rejector; abort=false)
    rejected = rejector(ele₀, μ)
    rejected && return false, true, ele₀, ele₀

    ele₁ = keplerian(cartesian(ele₀..., μ)..., μ)
    failed = !(isapprox(ele₁[1], ele₀[1], rtol=1.0) &&
               isapprox(ele₁[2], ele₀[2], atol=sqrt(eps())) &&
               isapprox(ele₁[3], ele₀[3], atol=sqrt(eps())) &&
               isapprox(ele₁[4], ele₀[4], atol=sqrt(eps())) &&
               isapprox(ele₁[5], ele₀[5], atol=sqrt(eps())) &&
               isapprox(abs(ele₁[6]), abs(ele₀[6]), atol=sqrt(eps())))
    if failed && abort
        error("Results do not match!\n",
              "a₁: ", rpad(ele₁[1], 24), "a₀:", ele₀[1], "\n",
              "e₁: ", rpad(ele₁[2], 24), "e₀:", ele₀[2], "\n",
              "i₁: ", rpad(ele₁[3], 24), "i₀:", ele₀[3], "\n",
              "Ω₁: ", rpad(ele₁[4], 24), "Ω₀:", ele₀[4], "\n",
              "ω₁: ", rpad(ele₁[5], 24), "ω₀:", ele₀[5], "\n",
              "ν₁: ", rpad(ele₁[6], 24), "ν₀:", ele₀[6])
    end
    failed, false, ele₁, ele₀
end

function run_tests(func, rngs...)
    iter = product(rngs...)
    total = length(iter)
    println("Testing $total combinations with $procs processes.")

    if isdefined(Main, :progress_map) && total > 1e6
        results = progress_map(func, iter, mapfun=pmap, progress=Progress(total, dt=2.0))
    else
        results = pmap(func, iter)
    end
    failed = reduce(+, map(first, results))
    rejected = reduce(+, map(x->x[2], results))
    passed = total - failed - rejected
    percent_failed = round(Int, failed / total * 100.0)
    percent_rejected = round(Int, rejected / total * 100.0)
    percent_passed = round(Int, passed / total * 100.0)
    println("Passed $passed ($percent_passed%), ",
            "failed $failed ($percent_failed%), ",
            "rejected $rejected ($percent_rejected%).")
end

println("Elliptic Orbits")
println("===============")

a_rng = range(7e6, 9e8, length=n)
e_rng = range(0.001, 0.999, length=n)
i_rng = range(0.001, 0.999π, length=n)
Ω_rng = range(0.001, 1.999π, length=n)
ω_rng = range(0.001, 1.999π, length=n)
ν_rng = range(-π, π, length=n)

t = @elapsed run_tests(x -> keplerian_roundtrip(x, μ, acceptor; abort=abort),
                       a_rng, e_rng, i_rng, Ω_rng, ω_rng, ν_rng)
println("In $t seconds.\n")

println("Hyperbolic Orbits")
println("=================")

a_rng = range(-2e6, -9e8, length=n)
e_rng = range(1.001, 5, length=n)
i_rng = range(0.001, 0.999π, length=n)
Ω_rng = range(0.001, 1.999π, length=n)
ω_rng = range(0.001, 1.999π, length=n)
ν_rng = range(-0.999π, 0.999π, length=n)

t = @elapsed run_tests(x -> keplerian_roundtrip(x, μ, hyperbolic_rejector; abort=abort),
                       a_rng, e_rng, i_rng, Ω_rng, ω_rng, ν_rng)
println("In $t seconds.\n")

println("Circular Orbits")
println("===============")

a_rng = range(7e6, 9e8, length=n)
i_rng = range(0.001, 0.999π, length=n)
Ω_rng = range(0.001, 1.999π, length=n)
ν_rng = range(-π, π, length=n)

t = @elapsed run_tests(x -> keplerian_roundtrip((x[1], 0.0, x[2], x[3], 0.0, x[4]), μ, acceptor; abort=abort),
                       a_rng, i_rng, Ω_rng, ν_rng)
println("In $t seconds.\n")

println("Equatorial Orbits")
println("=================")

a_rng = range(7e6, 9e8, length=n)
e_rng = range(0.001, 0.999, length=n)
ω_rng = range(0.001, 1.999π, length=n)
ν_rng = range(-π, π, length=n)

t = @elapsed run_tests(x -> keplerian_roundtrip((x[1], x[2], 0.0, 0.0, x[3], x[4]), μ, acceptor; abort=abort),
                       a_rng, e_rng, ω_rng, ν_rng)
println("In $t seconds.\n")

println("Circular-Equatorial Orbits")
println("==========================")

a_rng = range(7e6, 9e8, length=n)
ν_rng = range(-π, π, length=n)

t = @elapsed run_tests(x -> keplerian_roundtrip((x[1], 0.0, 0.0, 0.0, 0.0, x[2]), μ, acceptor; abort=abort),
                       a_rng, ν_rng)
println("In $t seconds.\n")


text:
################################################################################
### Tests
################################################################################

using Test

@testset "MyExpr" begin
    @test getfield(MyExpr(:(x + y)), :e) == Expr(:+, :x, :y)

    a = MyExpr(:(x * y))
    b = MyExpr(Expr(:*, :x, :y))
    @test a == b
    @test MyExpr(:*, :x, :y) == b

    e = quote
        a = 1 + c
        b = 2a
    end
    simplified = Expr(:block, Expr(:(=), :a, Expr(:+, 1, :c)), Expr(:(=), :b, Expr(:*, 2, :a)))
    @test getfield(MyExpr(e), :e) == simplified
end

@testset "Operators" begin
    op1 = create(:i, :l, name = :d)
    @test op1.name == :d
    @test op1.daggered == true
    @test op1.index == :i
    @test op1.time == :l

    op2 = annihilate(:j, :k, name = :a)
    @test op2.name == :a
    @test op2.daggered == false
    @test op2.index == :j
    @test op2.time == :k

    @test c(:i, :l, name=:d)' == op1
    @test c(:j, :k, name=:a) == op2

    op = c(:(i + N))
    @test op.index isa MyExpr
    @test getfield(op.index, :e) == Expr(:+, :i, :N)
end

@testset "ArrayElement" begin
    x = ArrayElement(:A, :i, :j)
    @test x.name == :A
    @test x.indices == [:i, :j]

    x = ArrayElement(:vector, :(2i + 1))
    @test x.indices[1] isa MyExpr
    @test getfield(x.indices[1], :e) == Expr(:+, Expr(:*, 2, :i), 1)
end

@testset "Basic math" begin
    a = c(:i)
    b = ArrayElement(:A, :j)
    @test a + b == MyExpr(:+, a, b)
    @test a * b == MyExpr(:*, a, b)
    @test -a == MyExpr(:*, -1, a)

    d = MyExpr(:(x + y))
    @test a + b + d == MyExpr(:+, a, b, d)
    @test a * b * d == MyExpr(:*, a, b, d)

    @test (a + b) * d == MyExpr(:*, MyExpr(:+, a, b), d)
    @test a + b * d == MyExpr(:+, a, MyExpr(:*, b, d))
end

@testset "Expand" begin
    a = c(:i)
    b = ArrayElement(:A, :j)
    e = (a + b) * 2
    @test e == MyExpr(:*, MyExpr(:+, a, b), 2)
    e2 = expand(e)
    @test e2 == MyExpr(:+, MyExpr(:*, a, 2), MyExpr(:*, b, 2))
end

@testset "Wicks" begin
    op1 = c(:(i))'
    op2 = c(:(i))
    op3 = c(:(i + N))'
    op4 = c(:(i + N))
    e = wicks(op1 * op2 * op3 * op4)
    result = MyExpr(
        :+,
        MyExpr(:*, ExpectationValue(op1 * op2), ExpectationValue(op3 * op4)),
        MyExpr(:*, -1, ExpectationValue(op1 * op3), ExpectationValue(op2 * op4)),
        MyExpr(:*, ExpectationValue(op1 * op4), ExpectationValue(op2 * op3))
    )
    @test e == result
end

@testset "to Greens elements" begin
    e = wicks(c(:(i))' * c(:(i)) * c(:(j))' * c(:(j)))
    g = to_greens(e)
    result =
        (1 - ArrayElement(:G00, :i, :i)) * (1 - ArrayElement(:G00, :j, :j)) +
        (ArrayElement(:I, :i, :j) - ArrayElement(:G00, :j, :i)) * ArrayElement(:G00, :i, :j)
    @test g == result
end


text:
export TabularLearner

"""
    TabularLearner{S, T}

Use a `Dict{S,Vector{T}}` to store action probabilities. `S` is the type of
state. `T` is the element type of probabilities.
"""
struct TabularLearner{S,T} <: AbstractLearner
    table::Dict{S,Vector{T}}
end

TabularLearner() = TabularLearner{Int,Float32}()
TabularLearner{S}() where {S} = TabularLearner{S,Float32}()
TabularLearner{S,T}() where {S,T} = TabularLearner(Dict{S,Vector{T}}())

(p::TabularLearner)(env::AbstractEnv) = p(ChanceStyle(env), env)

function (p::TabularLearner)(::ExplicitStochastic, env::AbstractEnv)
    if current_player(env) == chance_player(env)
        prob(env)
    else
        p(DETERMINISTIC, env)  # treat it just like a normal one
    end
end

function (t::TabularLearner)(::RLBase.AbstractChanceStyle, env::AbstractEnv)
    t(ActionStyle(env), env)
end

function (t::TabularLearner)(::FullActionSet, env::AbstractEnv)
    get!(t.table, state(env)) do
        m = legal_action_space_mask(env)
        m ./ sum(m)
    end
end

function (t::TabularLearner)(::MinimalActionSet, env::AbstractEnv)
    get!(t.table, state(env)) do
        n = length(action_space(env))
        fill(1 / n, n)
    end
end

"""
    update!(p::TabularLearner, state => prob)

!!! warn
    For environments of `FULL_ACTION_SET`, `prob` represents the probability
    distribution of `legal_action_space(env)`. For environments of
    `MINIMAL_ACTION_SET`, `prob` should represent the probability distribution
    of `action_space(env)`.
"""
RLBase.update!(p::TabularLearner, experience::Pair) =
    p.table[first(experience)] = last(experience)


text:
module WindowsTimeZoneIDs

using ...TimeZones: DEPS_DIR
using EzXML

if VERSION >= v"1.3"
    using ...TimeZones: @artifact_str
end

const UNICODE_CLDR_VERSION = "release-37"

# A mapping of Windows timezone names to Olson timezone names.
# Details on the contents of this file can be found at:
# http://cldr.unicode.org/development/development-process/design-proposals/extended-windows-olson-zid-mapping
const WINDOWS_ZONE_URL = "https://raw.githubusercontent.com/unicode-org/cldr/$UNICODE_CLDR_VERSION/common/supplemental/windowsZones.xml"
const WINDOWS_ZONE_FILE = joinpath("cldr-$UNICODE_CLDR_VERSION", "common", "supplemental", "windowsZones.xml")

const WINDOWS_XML_DIR = joinpath(DEPS_DIR, "local")
const WINDOWS_XML_FILE = joinpath(WINDOWS_XML_DIR, "windowsZones.xml")

isdir(WINDOWS_XML_DIR) || mkdir(WINDOWS_XML_DIR)

function compile(xml_file::AbstractString)
    # Get the timezone conversions from the file
    doc = readxml(xml_file)

    # Territory "001" is the global default
    map_zones = findall("//mapZone[@territory='001']", doc)

    # TODO: Eliminate the Etc/* POSIX names here? See Windows section of `localzone`

    # Dictionary to store the windows to time zone conversions
    translation = Dict{String,String}()
    for map_zone in map_zones
        win_name = map_zone["other"]
        posix_name = map_zone["type"]
        translation[win_name] = posix_name
    end

    return translation
end

const WINDOWS_TRANSLATION = if isfile(WINDOWS_XML_FILE)
    compile(WINDOWS_XML_FILE)
else
    Dict{AbstractString, AbstractString}()
end

function build(xml_file::AbstractString=WINDOWS_XML_FILE; force::Bool=false)
    fallback_xml_file = joinpath(WINDOWS_XML_DIR, "windowsZones2017a.xml")

    if !isfile(xml_file)
        if isfile(fallback_xml_file) && !force
            cp(fallback_xml_file, xml_file)
        else
            @info "Downloading Windows to POSIX timezone ID XML version: $UNICODE_CLDR_VERSION"
            @static if VERSION >= v"1.3"
                artifact_dir = @artifact_str "unicode-cldr-$UNICODE_CLDR_VERSION"
                xml_file = joinpath(artifact_dir, WINDOWS_ZONE_FILE)
            else
                download(WINDOWS_ZONE_URL, xml_file)
            end
        end
    end

    @info "Compiling Windows time zone name translation"
    translation = compile(xml_file)

    # Copy contents into translation constant
    empty!(WINDOWS_TRANSLATION)
    for (k, v) in translation
        WINDOWS_TRANSLATION[k] = v
    end
end

end


text:
module Memoization

using MacroTools: splitdef, combinedef, splitarg, combinearg, isexpr
export @memoize

# Stores a mapping:
#
#    (func, cache_type) => cache
#
# for each memoized function `func`. For top-level functions, `func` is the
# function's type. For closures or callables, `func` is an individual instance
# of the closure or callable object.
const caches = IdDict()

# Non-function "callables" or any function that has closed-over variables
# attached to it (indicated by having type parameters) are memoized
# per-instance.
memoize_per_instance(::Type{F}) where {F<:Function} = !isempty(F.parameters)
memoize_per_instance(::Type) = true

# Look up (and possibly create) the cache for a given memoized function. Using a
# generated function allows us to move this lookup to compile time for top-level
# functions and improve performance.
@generated function get_cache(func::F, cache_type=IdDict) where {F}
    if memoize_per_instance(F) || !(cache_type <: Type)
        quote
            _get!(cache_type, $caches, (func,cache_type))
        end
    else
        CT, = cache_type.parameters
        _get!(CT, caches, (F,CT))
    end
end

"""
    empty_cache!(arg)

Empties the memoization cache for functions, closures, or callables.

For functions or closures, `arg` should be the name of the function or closure.
For callables, `arg` should be a type and the cache for all callable objects
matching that type will be cleared.
"""
empty_cache!(func) = map(empty!, values(find_caches(func)))
find_caches(func::F) where {F<:Function} =
    filter((((func′,_),_),)->(memoize_per_instance(F) ? (func′ == func) : (func′ == F)), caches)
find_caches(F::Union{DataType,Union,UnionAll}) =
    filter((((func′,_),_),)->(func′ isa F), caches)
empty_all_caches!() = map(empty!, values(caches))

"""
    @memoize f(x) = ...
    @memoize Dict f(x) = ... # with custom cache type
    @memoize LRU(maxsize=4) f(x) = ... # with custom cache type
    @memoize (::Foo)(x) = ... # memoizing a callable
    f(x) = @memoize g(y) = ... # memoizing a closure

Memoize a function, closure, or callable, with respect to all of its arguments
and keyword arguments.

Memoized closures or callables are memoized on a per-instance basis, so closures
are free to use the closed over variables and callables are free to use the
fields of the callable object.

By default, an IdDict is used as a cache. Any dict-like object can be used by
passing a type or an expression to construct the object as the first argment to
the macro before the function definition. For example, if you want to memoize
based on the contents of vectors, you could use a `Dict`.
"""
macro memoize(ex1, ex2=nothing)
    cachetype, funcdef = ex2 == nothing ? ((), ex1) : ((ex1,), ex2)
    sdef = splitdef(funcdef)
    # if cachetype is a call, wrap it in a () -> ...
    if !isempty(cachetype) && isexpr(cachetype[1],:call)
        cachetype = (:(() -> $(cachetype[1])),)
    end
    # give unnamed args a placeholder name:
    sdef[:args] = map(sdef[:args]) do arg
        sarg = splitarg(arg)
        combinearg((sarg[1] == nothing ? gensym() : sarg[1]), sarg[2:end]...)
    end
    arg_signature   = [(issplat ? :($arg...) : arg)          for (arg,_,issplat) in map(splitarg,sdef[:args])]
    kwarg_signature = [(issplat ? :($arg...) : :($arg=$arg)) for (arg,_,issplat) in map(splitarg,sdef[:kwargs])]
    T, getter = (gensym.(("T","getter")))

    # if memoizing just `f(x) = ...` we want to call both `get_cache` and
    # `empty_cache` on `f`, but if memoizing a callable like
    # `(x::Foo{T})(args...) where {T} = ...`, we want to call get_cache on `x`
    # but empty_cache on `Foo{T} where {T}`
    if isexpr(sdef[:name], :(::))
        length(sdef[:name].args)==1 && pushfirst!(sdef[:name].args, gensym())
        cacheid_get   = sdef[:name].args[1]
        cacheid_empty = :($(sdef[:name].args[2]) where {$(sdef[:whereparams]...)})
    else
        cacheid_get = cacheid_empty = sdef[:name]
    end

    # the body of the function definition is replaced with this:
    sdef[:body] = quote
        ($getter)() = $(sdef[:body])
        $T = $(Core.Compiler.return_type)($getter, $Tuple{})
        $_get!($getter, $get_cache($cacheid_get, $(cachetype...)), (($(arg_signature...),),(;$(kwarg_signature...),))) :: $T
    end

    canary = gensym("canary")
    quote
        func = Core.@__doc__ $(esc(combinedef(sdef)))
        # only empty cache if this function definition is at the toplevel
        # see also: https://discourse.julialang.org/t/is-there-a-way-to-determine-whether-code-is-toplevel
        $(esc(canary)) = true
        if isdefined($__module__, $(QuoteNode(canary)))
            $empty_cache!($(esc(cacheid_empty)))
        end
        func
    end
end

_get!(args...) = get!(args...)
if VERSION < v"1.1.0-DEV.752"
    # this was only added in https://github.com/JuliaLang/julia/commit/7ba6c824467d2df51db6e091bbfc9e821e5a6dc2
    function _get!(default::Base.Callable, d::IdDict{K,V}, @nospecialize(key)) where {K, V}
        val = get(d, key, Base.secret_table_token)
        if val === Base.secret_table_token
            val = default()
            setindex!(d, val, key)
        end
        return val
    end
end

end


text:
"""
    FiniteGP{Tf<:AbstractGP, Tx<:AbstractVector, TΣy}

The finite-dimensional projection of the AbstractGP `f` at `x`. Assumed to be observed under
Gaussian noise with zero mean and covariance matrix `Σ`
"""
struct FiniteGP{Tf<:AbstractGP, Tx<:AbstractVector, TΣ} <: ContinuousMultivariateDistribution
    f::Tf
    x::Tx
    Σy::TΣ
end

function FiniteGP(f::AbstractGP, x::AbstractVector, σ²::AbstractVector{<:Real})
    return FiniteGP(f, x, Diagonal(σ²))
end

FiniteGP(f::AbstractGP, x::AbstractVector, σ²::Real) = FiniteGP(f, x, Fill(σ², length(x)))

FiniteGP(f::AbstractGP, x::AbstractVector) = FiniteGP(f, x, 1e-18)

function FiniteGP(
    f::AbstractGP,
    X::AbstractMatrix;
    obsdim::Int = KernelFunctions.defaultobs,
)
    return FiniteGP(f, KernelFunctions.vec_of_vecs(X; obsdim=obsdim))
end

Base.length(f::FiniteGP) = length(f.x)

(f::AbstractGP)(x...) = FiniteGP(f, x...)

"""
    mean(fx::FiniteGP)

Compute the mean vector of `fx`.

```jldoctest
julia> f = GP(Matern52Kernel());

julia> x = randn(11);

julia> mean(f(x)) == zeros(11)
true
```
"""
Statistics.mean(fx::FiniteGP) = mean(fx.f, fx.x)

"""
    cov(f::FiniteGP)

Compute the covariance matrix of `fx`.

## Noise-free observations

```jldoctest cov_finitegp
julia> f = GP(Matern52Kernel());

julia> x = randn(11);

julia> cov(f(x)) == kernelmatrix(Matern52Kernel(), x)
true
```

## Isotropic observation noise

```jldoctest cov_finitegp
julia> cov(f(x, 0.1)) == kernelmatrix(Matern52Kernel(), x) + 0.1 * I
true
```

## Independent anisotropic observation noise

```jldoctest cov_finitegp
julia> s = rand(11);

julia> cov(f(x, s)) == kernelmatrix(Matern52Kernel(), x) + Diagonal(s)
true
```

## Correlated observation noise

```jldoctest cov_finitegp
julia> A = randn(11, 11); S = A'A;

julia> cov(f(x, S)) == kernelmatrix(Matern52Kernel(), x) + S
true
```
"""
Statistics.cov(f::FiniteGP) = cov(f.f, f.x) + f.Σy

"""
    mean_and_cov(f::FiniteGP)

Equivalent to `(mean(f), cov(f))`, but sometimes more efficient to compute them jointly than
separately.

```jldoctest
julia> fx = GP(SqExponentialKernel())(range(-3.0, 3.0; length=10), 0.1);

julia> mean_and_cov(fx) == (mean(fx), cov(fx))
true
```
"""
function mean_and_cov(f::FiniteGP)
    m, C = mean_and_cov(f.f, f.x)
    return m, C + f.Σy
end

"""
    cov(fx::FiniteGP, gx::FiniteGP)

Compute the cross-covariance matrix between `fx` and `gx`.

```jldoctest
julia> f = GP(Matern32Kernel());

julia> x1 = randn(11);

julia> x2 = randn(13);

julia> cov(f(x1), f(x2)) == kernelmatrix(Matern32Kernel(), x1, x2)
true
```
"""
function Statistics.cov(fx::FiniteGP, gx::FiniteGP)
    @assert fx.f == gx.f
    return cov(fx.f, fx.x, gx.x)
end

"""
    marginals(f::FiniteGP)

Compute a vector of Normal distributions representing the marginals of `f` efficiently.
In particular, the off-diagonal elements of `cov(f(x))` are never computed.

```jldoctest
julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> fs = marginals(f(x));

julia> mean.(fs) == mean(f(x))
true

julia> std.(fs) == sqrt.(diag(cov(f(x))))
true
```
"""
marginals(f::FiniteGP) = Normal.(mean(f), sqrt.(cov_diag(f.f, f.x) .+ diag(f.Σy)))

"""
    rand(rng::AbstractRNG, f::FiniteGP, N::Int=1)

Obtain `N` independent samples from the marginals `f` using `rng`. Single-sample methods
produce a `length(f)` vector. Multi-sample methods produce a `length(f)` x `N` `Matrix`.

```jldoctest
julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> rand(f(x)) isa Vector{Float64}
true

julia> rand(MersenneTwister(123456), f(x)) isa Vector{Float64}
true

julia> rand(f(x), 3) isa Matrix{Float64}
true

julia> rand(MersenneTwister(123456), f(x), 3) isa Matrix{Float64}
true
```
"""
function Random.rand(rng::AbstractRNG, f::FiniteGP, N::Int)
    m, C_mat = mean_and_cov(f)
    C = cholesky(Symmetric(C_mat))
    return m .+ C.U' * randn(rng, promote_type(eltype(m), eltype(C)), length(m), N)
end
Random.rand(f::FiniteGP, N::Int) = rand(Random.GLOBAL_RNG, f, N)
Random.rand(rng::AbstractRNG, f::FiniteGP) = vec(rand(rng, f, 1))
Random.rand(f::FiniteGP) = vec(rand(f, 1))

"""
    logpdf(f::FiniteGP, y::AbstractVecOrMat{<:Real})

The logpdf of `y` under `f` if is `y isa AbstractVector`. logpdf of each column of `y` if
`y isa Matrix`.

```jldoctest
julia> f = GP(Matern32Kernel());

julia> x = randn(11);

julia> y = rand(f(x));

julia> logpdf(f(x), y) isa Real
true

julia> Y = rand(f(x), 3);

julia> logpdf(f(x), Y) isa AbstractVector{<:Real}
true
```
"""
function Distributions.logpdf(f::FiniteGP, y::AbstractVector{<:Real})
    return first(logpdf(f, reshape(y, :, 1)))
end

function Distributions.logpdf(f::FiniteGP, Y::AbstractMatrix{<:Real})
    m, C_mat = mean_and_cov(f)
    C = cholesky(Symmetric(C_mat))
    T = promote_type(eltype(m), eltype(C), eltype(Y))
    return -((size(Y, 1) * T(log(2π)) + logdet(C)) .+ diag_Xt_invA_X(C, Y .- m)) ./ 2
end

"""
   elbo(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

The Titsias Evidence Lower BOund (ELBO) [1]. `y` are observations of `f`, and `u`
are pseudo-points, where `u = f(z)` for some `z`.

```jldoctest
julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

julia> z = range(-5.0, 5.0; length=13);

julia> y = rand(f(x, 0.1));

julia> elbo(f(x, 0.1), y, f(z)) < logpdf(f(x, 0.1), y)
true
```

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian
processes". In: Proceedings of the Twelfth International Conference on Artificial
Intelligence and Statistics. 2009.
"""
function elbo(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
    _dtc, chol_Σy, A = _compute_intermediates(f, y, u)
    return _dtc - (tr_Cf_invΣy(f, f.Σy, chol_Σy) - sum(abs2, A)) / 2
end

"""
    dtc(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

The Deterministic Training Conditional (DTC) [1]. `y` are observations of `f`, and `u`
are pseudo-points.

```jldoctest
julia> f = GP(Matern52Kernel());

julia> x = randn(1000);

julia> z = range(-5.0, 5.0; length=256);

julia> y = rand(f(x, 0.1));

julia> isapprox(dtc(f(x, 0.1), y, f(z)), logpdf(f(x, 0.1), y); atol=1e-3, rtol=1e-3)
true
```

[1] - M. Seeger, C. K. I. Williams and N. D. Lawrence. "Fast Forward Selection to Speed Up
Sparse Gaussian Process Regression". In: Proceedings of the Ninth International Workshop on
Artificial Intelligence and Statistics. 2003
"""
function dtc(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
    return first(_compute_intermediates(f, y, u))
end

# Factor out computations common to the `elbo` and `dtc`.
function _compute_intermediates(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
    consistency_check(f, y, u)
    chol_Σy = cholesky(f.Σy)

    A = cholesky(Symmetric(cov(u))).U' \ (chol_Σy.U' \ cov(f, u))'
    Λ_ε = cholesky(Symmetric(A * A' + I))
    δ = chol_Σy.U' \ (y - mean(f))

    tmp = logdet(chol_Σy) + logdet(Λ_ε) + sum(abs2, δ) - sum(abs2, Λ_ε.U' \ (A * δ))
    _dtc = -(length(y) * typeof(tmp)(log(2π)) + tmp) / 2
    return _dtc, chol_Σy, A
end

function consistency_check(f, y, u)
    @assert length(f) == size(y, 1)
end

function tr_Cf_invΣy(f::FiniteGP, Σy::Diagonal, chol_Σy::Cholesky)
    return sum(cov_diag(f.f, f.x) ./ diag(Σy))
end


text:
module Shaders
  using CG, Raytracing, SceneObjects, Lighting;
  export hitShader, lambertShader, blinnPhongShader;

  function hitShader(ray::Ray, scene::Scene, lights::SceneLights)
    intersection, t, = Raytracing.intersect(ray, scene);
    if intersection
      return 1.0f0;
    else
      return 0.0f0;
    end
  end

  function lambertShader(ray::Ray, scene::Scene, lights::SceneLights)
    intersection, t, o = Raytracing.intersect(ray, scene);

    if (!intersection)
      return 0.0f0;
    end

    p = ray.origin + t * ray.direction;
    n = Raytracing.surfaceNormal(ray, t, o);

    L = 0f0;
    for l in lights.lights
      L = L + 1.0f0 * max(0, n * l.position);
    end
    return L;
  end

  function blinnPhongShader(ray::Ray, scene::Scene, lights::SceneLights)
    hit, t, object = Raytracing.intersect(ray, scene)
    if hit
      # compute surface normal
      normal = Raytracing.surfaceNormal(ray,t,object)
      # intersection point of camera ray
      p = ray.origin+t*ray.direction
      # initial shade of object
      shade = 1.0f0
      for light in lights.lights
        l = light.position;
        # direction intersection point to light source
        d = unitize(l-p)
        # dot product of normal vector and direction to light
        s = normal * d;
        # phong
        v = unitize(ray.origin - p);
        h = unitize(l + v);
        phong = 1000;
        specular = 1.0f0 * max(0.0f0, (normal * h)^phong)
        # add light if light is in front of surface
        shade += specular + max(0.0f0,s);
      end
      return shade
    else
      return 0.0f0
    end
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
# Command: Close Window
"""
	window_close!(session::Session)::Nothing

The window handle associated with the current top-level browsing context.
"""
function window_close!(session::Session)::Nothing
    @unpack addr, id = session
    response =
        HTTP.delete("$addr/session/$id/window", [("Content-Type" => "application/json")])
    @assert response.status == 200
    nothing
end


text:
#!/usr/bin/julia
# Renames the FASTQ files in a directory by prepending the sample number based
# on the embedded sample number

# Activate the proper packages
using Tk

# Prompt for the folder containing the FASTQs
fastq_folder = ChooseDirectory()

# Get all of the files
fastqs = readdir(fastq_folder)

# Iterate through each file
for fastq in fastqs
    # Get the full path
    fastq_path = joinpath(fastq_folder, fastq)

    # Delete if this is a JSON file
    if last(fastq,5) == ".json"
        rm(fastq_path)
        println(string("Deleting ", fastq_path))
        continue
    end

    # Find the sample number in the filename
    startS = findfirst("_S", fastq)
    startL = findfirst("_L", fastq)

    # Extract the sample number
    sampleId = fastq[(startS[1]+1):(startL[1]-1)]

    # Convert the number to three digits, no matter what
    sampleNum = parse(Int64, sampleId[2:end])
    sampleId = string("S", lpad(sampleNum, 3, "0"))

    # Construct the new filename
    fastq_newname = string(sampleId, "_", fastq)

    # Rename the file
    mv(fastq_path, joinpath(fastq_folder, fastq_newname))
    println(string("Renaming ", fastq_path, " to ", joinpath(fastq_folder, fastq_newname)))

end

max_forks_repo_name: MillironX/wphl-bioinformatics


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

function Base.show(io::IO, fst::WFST)
  sz = size(fst)
  println("WFST #states: $(sz[1]), #arcs: $(sz[2]), #isym: $(length(fst.isym)), #osym: $(length(fst.osym))")
  if sz[2] < 100 && !isempty(fst)
    int2isym = get_iisym(fst)
    int2osym = get_iosym(fst)
    for (i,s) in enumerate(fst)
      if isinitial(fst,i) && isfinal(fst,i)
        wi = get_initial_weight(fst,i)
        wf = get_final_weight(fst,i)
        println("|(($(i)/$(wi*wf)))|")
      elseif isfinal(fst,i)
        w = get_final_weight(fst,i)
        println("(($(i)/$(w)))")
      elseif isinitial(fst,i)
        w = get_initial_weight(fst,i)
        println("|$(i)/$(w)|")
      else
        println("($(i))")
      end
      show_state(io, s, int2isym, int2osym)
    end
  end
end

function show_state(io::IO, s, int2isym, int2osym)
  for arc in s
    show(io, arc, int2isym, int2osym)
    println()
  end
end

function Base.show(io::IO, p::Path)
  if length(p.ilabels) + length(p.olabels) < 50
    isym = get_iisym(p)
    osym = get_iosym(p)
    print(io, "$([isym[i] for i in p.ilabels]) → $([osym[i] for i in p.olabels]), w:$(p.weight)" )
  else
    print(io, " #$(length(p.ilabels)) → #$(length(p.olabels)), w:$(p.weight)" )
  end
end


text:
# Automatically generated using Clang.jl

const CUBLAS_VER_MAJOR = 11
const CUBLAS_VER_MINOR = 3
const CUBLAS_VER_PATCH = 1
const CUBLAS_VER_BUILD = 68
const CUBLAS_VERSION = CUBLAS_VER_MAJOR * 1000 + CUBLAS_VER_MINOR * 100 + CUBLAS_VER_PATCH

@cenum cublasStatus_t::UInt32 begin
    CUBLAS_STATUS_SUCCESS = 0
    CUBLAS_STATUS_NOT_INITIALIZED = 1
    CUBLAS_STATUS_ALLOC_FAILED = 3
    CUBLAS_STATUS_INVALID_VALUE = 7
    CUBLAS_STATUS_ARCH_MISMATCH = 8
    CUBLAS_STATUS_MAPPING_ERROR = 11
    CUBLAS_STATUS_EXECUTION_FAILED = 13
    CUBLAS_STATUS_INTERNAL_ERROR = 14
    CUBLAS_STATUS_NOT_SUPPORTED = 15
    CUBLAS_STATUS_LICENSE_ERROR = 16
end

@cenum cublasFillMode_t::UInt32 begin
    CUBLAS_FILL_MODE_LOWER = 0
    CUBLAS_FILL_MODE_UPPER = 1
    CUBLAS_FILL_MODE_FULL = 2
end

@cenum cublasDiagType_t::UInt32 begin
    CUBLAS_DIAG_NON_UNIT = 0
    CUBLAS_DIAG_UNIT = 1
end

@cenum cublasSideMode_t::UInt32 begin
    CUBLAS_SIDE_LEFT = 0
    CUBLAS_SIDE_RIGHT = 1
end

@cenum cublasOperation_t::UInt32 begin
    CUBLAS_OP_N = 0
    CUBLAS_OP_T = 1
    CUBLAS_OP_C = 2
    CUBLAS_OP_HERMITAN = 2
    CUBLAS_OP_CONJG = 3
end

@cenum cublasPointerMode_t::UInt32 begin
    CUBLAS_POINTER_MODE_HOST = 0
    CUBLAS_POINTER_MODE_DEVICE = 1
end

@cenum cublasAtomicsMode_t::UInt32 begin
    CUBLAS_ATOMICS_NOT_ALLOWED = 0
    CUBLAS_ATOMICS_ALLOWED = 1
end

@cenum cublasGemmAlgo_t::Int32 begin
    CUBLAS_GEMM_DFALT = -1
    CUBLAS_GEMM_DEFAULT = -1
    CUBLAS_GEMM_ALGO0 = 0
    CUBLAS_GEMM_ALGO1 = 1
    CUBLAS_GEMM_ALGO2 = 2
    CUBLAS_GEMM_ALGO3 = 3
    CUBLAS_GEMM_ALGO4 = 4
    CUBLAS_GEMM_ALGO5 = 5
    CUBLAS_GEMM_ALGO6 = 6
    CUBLAS_GEMM_ALGO7 = 7
    CUBLAS_GEMM_ALGO8 = 8
    CUBLAS_GEMM_ALGO9 = 9
    CUBLAS_GEMM_ALGO10 = 10
    CUBLAS_GEMM_ALGO11 = 11
    CUBLAS_GEMM_ALGO12 = 12
    CUBLAS_GEMM_ALGO13 = 13
    CUBLAS_GEMM_ALGO14 = 14
    CUBLAS_GEMM_ALGO15 = 15
    CUBLAS_GEMM_ALGO16 = 16
    CUBLAS_GEMM_ALGO17 = 17
    CUBLAS_GEMM_ALGO18 = 18
    CUBLAS_GEMM_ALGO19 = 19
    CUBLAS_GEMM_ALGO20 = 20
    CUBLAS_GEMM_ALGO21 = 21
    CUBLAS_GEMM_ALGO22 = 22
    CUBLAS_GEMM_ALGO23 = 23
    CUBLAS_GEMM_DEFAULT_TENSOR_OP = 99
    CUBLAS_GEMM_DFALT_TENSOR_OP = 99
    CUBLAS_GEMM_ALGO0_TENSOR_OP = 100
    CUBLAS_GEMM_ALGO1_TENSOR_OP = 101
    CUBLAS_GEMM_ALGO2_TENSOR_OP = 102
    CUBLAS_GEMM_ALGO3_TENSOR_OP = 103
    CUBLAS_GEMM_ALGO4_TENSOR_OP = 104
    CUBLAS_GEMM_ALGO5_TENSOR_OP = 105
    CUBLAS_GEMM_ALGO6_TENSOR_OP = 106
    CUBLAS_GEMM_ALGO7_TENSOR_OP = 107
    CUBLAS_GEMM_ALGO8_TENSOR_OP = 108
    CUBLAS_GEMM_ALGO9_TENSOR_OP = 109
    CUBLAS_GEMM_ALGO10_TENSOR_OP = 110
    CUBLAS_GEMM_ALGO11_TENSOR_OP = 111
    CUBLAS_GEMM_ALGO12_TENSOR_OP = 112
    CUBLAS_GEMM_ALGO13_TENSOR_OP = 113
    CUBLAS_GEMM_ALGO14_TENSOR_OP = 114
    CUBLAS_GEMM_ALGO15_TENSOR_OP = 115
end

@cenum cublasMath_t::UInt32 begin
    CUBLAS_DEFAULT_MATH = 0
    CUBLAS_TENSOR_OP_MATH = 1
    CUBLAS_PEDANTIC_MATH = 2
    CUBLAS_TF32_TENSOR_OP_MATH = 3
    CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION = 16
end

const cublasDataType_t = cudaDataType

@cenum cublasComputeType_t::UInt32 begin
    CUBLAS_COMPUTE_16F = 64
    CUBLAS_COMPUTE_16F_PEDANTIC = 65
    CUBLAS_COMPUTE_32F = 68
    CUBLAS_COMPUTE_32F_PEDANTIC = 69
    CUBLAS_COMPUTE_32F_FAST_16F = 74
    CUBLAS_COMPUTE_32F_FAST_16BF = 75
    CUBLAS_COMPUTE_32F_FAST_TF32 = 77
    CUBLAS_COMPUTE_64F = 70
    CUBLAS_COMPUTE_64F_PEDANTIC = 71
    CUBLAS_COMPUTE_32I = 72
    CUBLAS_COMPUTE_32I_PEDANTIC = 73
end

const cublasContext = Cvoid
const cublasHandle_t = Ptr{cublasContext}
const cublasLogCallback = Ptr{Cvoid}
const cublasXtContext = Cvoid
const cublasXtHandle_t = Ptr{cublasXtContext}

@cenum cublasXtPinnedMemMode_t::UInt32 begin
    CUBLASXT_PINNING_DISABLED = 0
    CUBLASXT_PINNING_ENABLED = 1
end

@cenum cublasXtOpType_t::UInt32 begin
    CUBLASXT_FLOAT = 0
    CUBLASXT_DOUBLE = 1
    CUBLASXT_COMPLEX = 2
    CUBLASXT_DOUBLECOMPLEX = 3
end

@cenum cublasXtBlasOp_t::UInt32 begin
    CUBLASXT_GEMM = 0
    CUBLASXT_SYRK = 1
    CUBLASXT_HERK = 2
    CUBLASXT_SYMM = 3
    CUBLASXT_HEMM = 4
    CUBLASXT_TRSM = 5
    CUBLASXT_SYR2K = 6
    CUBLASXT_HER2K = 7
    CUBLASXT_SPMM = 8
    CUBLASXT_SYRKX = 9
    CUBLASXT_HERKX = 10
    CUBLASXT_TRMM = 11
    CUBLASXT_ROUTINE_MAX = 12
end


text:
RANDOM_SEED = 809123
rng = MersenneTwister(RANDOM_SEED)

@testset "Particle" begin
    @test_throws AssertionError Particle(zeros(3), zeros(2), zeros(1), 0.0, 1.0)

    # * Test type promotion
    @test let
        p = Particle(15, 20.0, 10, rng)
        eltype(p.x) ≡ eltype(p.v)
    end
end

@testset "Population" begin
    @test let
        pops = typeof(Population(15, 20, -10.0, 10.0, rng))
        actual_type = typeof(Vector{Particle}(undef, 15))
        pops ≡ actual_type
    end

    # * Test random seeds
    @test let
        pops = typeof(Population(15, 20, -10.0, 10.0, rng))
        actual_type = typeof(Vector{Particle}(undef, 15))
        pops ≡ actual_type
    end

    # * Test multiple ranges
    @test let
        range_a = SVector(-10.0, 10.0)
        range_b = SVector(-2.5, 2.0)
        pops = Population(20, [range_a, range_b], rng)
        type_pops = typeof(pops)
        actual_type = typeof(Vector{Particle}(undef, 2))
        type_pops ≡ actual_type
    end
    @test_throws AssertionError Population(15, -1, 1.0, 1.0, rng)
    @test_throws AssertionError Population(0, 5, 1.0, 1.0, rng)
    @test_throws AssertionError Population(-1, 5, 1.0, 1.0, rng)

    @test let
        range_a = SVector(-10.0, 10.0)
        range_b = SVector(-2.5, 2.0)
        pops = Population(20, [range_a, range_b], rng)
        type_pops = typeof(pops)
        actual_type = typeof(Vector{Particle}(undef, 2))
        type_pops ≡ actual_type
    end
end


text:
import GaussianProcesses: GPE, AbstractMatrix, predict_f

function cliff_face(gpT::GPE, gpC::GPE, sentinels::AbstractMatrix)
    pred_T = predict_f(gpT, sentinels; full_cov=true)
    pred_C = predict_f(gpC, sentinels; full_cov=true)
    μposterior = pred_T[1].-pred_C[1]
    Σposterior = pred_T[2]+pred_C[2]
    return μposterior, Σposterior
end

function sim_cliff(gpT::GPE, gpC::GPE, gpNull::GPE, treat::BitVector, X∂::AbstractMatrix)
    μ = mean(gpNull.m, gpNull.x)
    null = MultivariateNormal(μ, gpNull.cK)
    Ysim = rand(null)

    gpT.y = Ysim[treat]
    gpC.y = Ysim[.!treat]

    update_alpha!(gpT)
    update_alpha!(gpC)

    return cliff_face(gpT, gpC, X∂)
end
function nsim_cliff(gpT::GPE, gpC::GPE, X∂::AbstractMatrix, nsim::Int)
    gpT_mod = modifiable(gpT)
    gpC_mod = modifiable(gpC)
    yNull = [gpT.y; gpC.y]
    mNull = gpT.m
    gpNull = GPE([gpT.x gpC.x], yNull, mNull, gpT.kernel, gpT.logNoise)
    treat = BitVector(gpNull.nobs)
    treat[:] = false
    treat[1:gpT.nobs] = true
    mT = gpT_mod.m
    mC = gpC_mod.m
    cliff_sims = [sim_cliff(gpT_mod, gpC_mod, gpNull, treat, X∂)
                  for _ in 1:nsim];
    return cliff_sims
end


text:
using LightXML

" Name mappings from JuliaKara.jl to classic Kara"
const XML_NAMES = Dict(
    :world=> "XmlWorld",
    :tree => "XmlWallPoints",
    :mushroom => "XmlObstaclePoints",
    :leaf => "XmlPaintedfieldPoints",
    :kara => "XmlKaraList"
)

"Directions mappings from classic JuliaKara to Kara.jl"
const XML_DIR = Dict(
    0 => 1,
    1 => 4,
    2 => 3,
    3 => 2,
)

function xml_map_direction(xml_index::Int)
    JuliaKara_noGUI.ActorsWorld.DIRECTIONS[XML_DIR[xml_index]]
end

function kara_map_direction(direction::Symbol)
    for (k,d) in enumerate(JuliaKara_noGUI.ActorsWorld.DIRECTIONS)
        if d == direction
            for p in XML_DIR
                if p.second == k
                    return p.first
                end
            end
        end
    end
end

function xml_map_y(wo::World,y::Int)
    # Original Kara starts counting at 0 in the left upper corner
    wo.size.height - y
end

function xml_map_x(wo::World,x::Int)
    x + 1
end

function kara_map_y(wo::World,y::Int)
    wo.size.height - y
end

function kara_map_x(wo::World,x::Int)
    x - 1
end

function xml_parse_tree!(wo::World,element::XMLElement)
    place_tree(
        wo,
        xml_map_x(wo,parse(Int,attribute(element,"x"))),
        xml_map_y(wo,parse(Int,attribute(element,"y")))
    )
end

function xml_parse_mushroom!(wo::World,element::XMLElement)
    place_mushroom(
        wo,
        xml_map_x(wo,parse(Int,attribute(element,"x"))),
        xml_map_y(wo,parse(Int,attribute(element,"y")))
    )
end

function xml_parse_leaf!(wo::World,element::XMLElement)
    place_leaf(
        wo,
        xml_map_x(wo,parse(Int,attribute(element,"x"))),
        xml_map_y(wo,parse(Int,attribute(element,"y")))
    )
end

function xml_parse_kara!(wo::World,element::XMLElement)
    place_kara(
        wo,
        xml_map_x(wo,parse(Int,attribute(element,"x"))),
        xml_map_y(wo,parse(Int,attribute(element,"y"))),
        parse(Int,attribute(element,"direction")) |> xml_map_direction
    )
end

function xml_generate_world(element::XMLElement)
    World(
        parse(Int,attribute(element,"sizex")),
        parse(Int,attribute(element,"sizey"))
    )
end

function xml_parse_actor!(wo::World,elements::Vector{XMLElement},
                          element_name::String,parser::Function)
    for el in elements
        for p in get_elements_by_tagname(el,element_name)
            parser(wo,p)
        end
    end
end

function xml_load_world(path::AbstractString)
    xworld = parse_file(path)
    xworld_def = root(xworld)

    world = xml_generate_world(xworld_def)

    xtree = get_elements_by_tagname(xworld_def,XML_NAMES[:tree])
    xmushroom = get_elements_by_tagname(xworld_def,XML_NAMES[:mushroom])
    xleaf = get_elements_by_tagname(xworld_def,XML_NAMES[:leaf])
    xkara = get_elements_by_tagname(xworld_def,XML_NAMES[:kara])

    xml_parse_actor!(
        world,
        xtree,
        "XmlPoint",
        xml_parse_tree!
    )

    xml_parse_actor!(
        world,
        xmushroom,
        "XmlPoint",
        xml_parse_mushroom!
    )

    xml_parse_actor!(
        world,
        xleaf,
        "XmlPoint",
        xml_parse_leaf!
    )

    xml_parse_actor!(
        world,
        xkara,
        "XmlKara",
        xml_parse_kara!
    )

    return world
end

function xml_create_world(wo::World,xdoc::XMLDocument)
    xworld_def = create_root(xdoc,XML_NAMES[:world])
    set_attribute(xworld_def,"sizex",wo.size.width)
    set_attribute(xworld_def,"sizey",wo.size.height)
    set_attribute(xworld_def,"version","KaraX 1.0 javakara")
    return xworld_def
end

function xml_save_world(wo::World,path::AbstractString)
    xdoc = XMLDocument()
    # create root element
    xworld_def  = xml_create_world(wo,xdoc)
    # create child groups
    xtree = new_child(xworld_def,XML_NAMES[:tree])
    xmushroom = new_child(xworld_def,XML_NAMES[:mushroom])
    xleaf = new_child(xworld_def,XML_NAMES[:leaf])
    xkara = new_child(xworld_def,XML_NAMES[:kara])
    # iterate over actors
    for ac in wo.actors
        if ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]
            t = new_child(xkara,"XmlKara")
            set_attributes(t,Dict(
                "direction" => kara_map_direction(ac.orientation.value),
                "name" => "Kara",
                "x" => kara_map_x(wo,ac.location.x) |> string,
                "y" => kara_map_y(wo,ac.location.y) |> string
            ))
        elseif ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:mushroom]
            t = new_child(xmushroom,"XmlPoint")
            set_attributes(t,Dict(
                "x" => kara_map_x(wo,ac.location.x) |> string,
                "y" => kara_map_y(wo,ac.location.y) |> string
            ))
        elseif ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:tree]
            t = new_child(xtree,"XmlPoint")
            set_attributes(t,Dict(
                "x" => kara_map_x(wo,ac.location.x) |> string,
                "y" => kara_map_y(wo,ac.location.y) |> string
            ))
        elseif ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:leaf]
            t = new_child(xleaf,"XmlPoint")
            set_attributes(t,Dict(
                "type" => "0",
                "x" =>kara_map_x(wo,ac.location.x) |> string,
                "y" =>kara_map_y(wo,ac.location.y) |> string,
            ))
        else
            error("Missing actor definition, cant create XML Element")
        end
    end
    save_file(xdoc,path)
end


text:
using Zygote, NNlib, Test, Random, LinearAlgebra
using Zygote: gradient
using NNlib: conv
import Random

function ngradient(f, xs::AbstractArray...)
  grads = zero.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

gradcheck(f, xs...) =
  all(isapprox.(ngradient(f, xs...),
                gradient(f, xs...), rtol = 1e-5, atol = 1e-5))

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)

Random.seed!(0)

@test gradient(//, 2, 3) === (1//3, -2//9)

@test gradtest((x, W, b) -> σ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> σ.(W*x .+ b), (5,3), (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), (5,3), (2,5), 2)

@test gradtest((w, x) -> w'*x, randn(10, 2), randn(10))
@test gradtest((w, x) -> transpose(w)*x, randn(5,5), randn(5,5))

@test gradtest(x -> sum(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> prod(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> prod(x), (3,4,5))

@test gradtest(x -> softmax(x).*(1:3), 3)
@test gradtest(x -> softmax(x).*(1:3), (3,5))
@test gradtest(x -> logsoftmax(x).*(1:3), 3)
@test gradtest(x -> logsoftmax(x).*(1:3), (3,5))

@test gradtest(conv, rand(10, 3, 2), randn(Float64,2, 3, 2))
@test gradtest(conv, rand(10, 10, 3, 2), randn(Float64,2, 2, 3, 2))
@test gradtest(conv, rand(10, 10, 10, 3, 2), randn(Float64,2, 2, 2, 3, 2))

@test gradtest(x -> maxpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> maxpool(x, (2,2,2)), rand(10, 10, 10, 3, 2))

@test gradtest(x -> meanpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> meanpool(x, (2,2,2)), rand(5, 5, 5, 3, 2))

@test gradtest(x -> permutedims(x, [3,1,2]), rand(4,5,6))

@test gradtest(x -> repeat(x; inner=2), rand(5))
@test gradtest(x -> repeat(x; inner=2, outer=3), rand(5))
@test gradtest(x -> repeat(x; inner=(2,2,1), outer=(1,1,3)), rand(5,4,3))

@testset "dot" begin
  rng = MersenneTwister(123456)
  @test gradtest((x, y)->dot(x[1], y[1]), [randn(rng)], [randn(rng)])
  @test gradtest(dot, randn(rng, 10), randn(rng, 10))
  @test gradtest(dot, randn(rng, 10, 3), randn(rng, 10, 3))
end

@test gradtest(kron, rand(5), rand(3))
@test gradtest(kron, rand(5), rand(3), rand(8))
@test gradtest(kron, rand(5,1), rand(3,1))
@test gradtest(kron, rand(5,1), rand(3,1), rand(8,1))
@test gradtest(kron, rand(5,2), rand(3,2), rand(8,2))

@testset "maximum" begin
  @test gradtest(maximum, rand(2, 3))

  @test gradtest(x -> maximum(x, dims=1), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=2), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> maximum(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "minimum" begin
  @test gradtest(minimum, rand(2, 3))

  @test gradtest(x -> minimum(x, dims=1), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=2), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> minimum(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "backsolve" begin
  rng, P, Q = MersenneTwister(123456), 10, 9
  X, Y, y = randn(rng, P, P), randn(rng, P, Q), randn(rng, P)

  # \
  @test gradtest(X -> X \ Y, X)
  @test gradtest(Y -> X \ Y, Y)
  @test gradtest(X -> X \ y, X)
  @test gradtest(y -> X \ y, y)

  # /
  @test gradtest(X -> Y' / X, X)
  @test gradtest(Y -> Y' / X, Y)
  @test gradtest(X -> y' / X, X)
  @test gradtest(y -> y' / X, y)

  @testset "Cholesky" begin

    # Check that the forwards pass computes the correct thing.
    @test Zygote.forward(X->cholesky(X * X' + I) \ Y, X)[1] == cholesky(X * X' + I) \ Y
    @test gradtest(X->cholesky(X * X' + I) \ Y, X)
    @test gradtest(Y->cholesky(X * X' + I) \ Y, Y)
    @test gradtest(X->cholesky(X * X' + I) \ y, X)
    @test gradtest(y->cholesky(X * X' + I) \ y, y)
  end
end

@testset "Symmetric" begin
  rng, P = MersenneTwister(123456), 7
  A = randn(rng, P, P)
  @test gradtest(Symmetric, A)
end

@testset "diag" begin
  rng, P = MersenneTwister(123456), 10
  A = randn(rng, P, P)
  @test gradtest(diag, A)
end

@testset "dense + UniformScaling" begin
  rng = MersenneTwister(123456)
  A, λ = randn(rng, 10, 10), randn(rng)
  @test gradtest(A->A + 5I, A)
  @test gradtest(λ->A + λ[1] * I, [λ])
end

@testset "cholesky" begin
  rng, N = MersenneTwister(123456), 5
  A = randn(rng, N, N)
  @test gradtest(A->logdet(cholesky(A' * A + 1e-6I)), A)
  @testset "cholesky - scalar" begin
    y, back = Zygote.forward(cholesky, 5.0 * ones(1, 1))
    y′, back′ = Zygote.forward(cholesky, 5.0)
    C̄ = randn(rng, 1, 1)
    @test back′((factors=C̄,))[1] isa Real
    @test back′((factors=C̄,))[1] ≈ back((factors=C̄,))[1][1, 1]
  end
end

using Distances

Zygote.refresh()

@testset "distances" begin
  rng, P, Q, D = MersenneTwister(123456), 10, 9, 8

   # Check sqeuclidean.
  let
    x, y = randn(rng, D), randn(rng, D)
    @test gradtest(x->sqeuclidean(x, y), x)
    @test gradtest(y->sqeuclidean(x, y), y)
  end

   # Check binary colwise.
  let
    X, Y = randn(rng, D, P), randn(rng, D, P)
    @test gradtest(X->colwise(SqEuclidean(), X, Y), X)
    @test gradtest(Y->colwise(SqEuclidean(), X, Y), Y)
  end

   # Check binary pairwise.
  let
    X, Y = randn(rng, D, P), randn(rng, D, Q)
    @test gradtest(X->pairwise(SqEuclidean(), X, Y), X)
    @test gradtest(Y->pairwise(SqEuclidean(), X, Y), Y)
  end

   # Check unary pairwise.
  @test gradtest(X->pairwise(SqEuclidean(), X), randn(rng, D, P))
end

function cat_test(f, A::Union{AbstractVector, AbstractMatrix}...)
  @test gradtest(f, A...)
  Z, back = Zygote.forward(f, A...)
  Ā = back(randn(size(Z)))
  @test all(map((a, ā)->ā isa typeof(a), A, Ā))
end

@testset "vcat" begin

  # Vector-only.
  cat_test(vcat, randn(1))
  cat_test(vcat, randn(10))

  cat_test(vcat, randn(1), randn(1))
  cat_test(vcat, randn(5), randn(1))
  cat_test(vcat, randn(1), randn(6))
  cat_test(vcat, randn(5), randn(6))

  # Matrix-only.
  for c in [1, 2], r1 in [1, 2], r2 in [1, 2]
    cat_test(vcat, randn(r1, c))
    cat_test(vcat, randn(r1, c), randn(r2, c))
  end

  # Matrix-Vector / Vector-Matrix / Vector-Matrix-Vector.
  for r in [1, 2]
    cat_test(vcat, randn(r, 1), randn(r))
    cat_test(vcat, randn(r), randn(r, 1))
    cat_test(vcat, randn(r), randn(r, 1), randn(r))
  end
end

@testset "hcat" begin

  # Vector-only.
  for r in [1, 2]
    cat_test(hcat, randn(r))
    cat_test(hcat, randn(r), randn(r))
  end

  # Matrix-only.
  for r in [1, 2], c1 in [1, 2], c2 in [1, 2]
    cat_test(hcat, randn(r, c1), randn(r, c2))
  end

  # Matrix-Vector / Vector-Matrix / Vector-Matrix-Vector.
  for r in [1, 2], c in [1, 2]
    cat_test(hcat, randn(r, c), randn(r))
    cat_test(hcat, randn(r), randn(r, c))
    cat_test(hcat, randn(r), randn(r, c), randn(r))
  end
end

@testset "one(s) and zero(s)" begin
  @test Zygote.gradient(x->sum(ones(size(x))), randn(5))[1] isa Nothing
  @test Zygote.gradient(x->sum(one(x)), randn(3, 3))[1] isa Nothing
  @test Zygote.gradient(x->sum(zeros(size(x))), randn(7))[1] isa Nothing
  @test Zygote.gradient(x->sum(zero(x)), randn(3))[1] isa Nothing
end

import StatsFuns

Zygote.refresh()

@testset "xlogx" begin
  @test gradcheck(x->2.5 * StatsFuns.xlogx(x[1]), [1.0])
  @test gradcheck(x->2.5 * StatsFuns.xlogx(x[1]), [2.45])
end

@testset "logistic" begin
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [-5.0])
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [-1.0])
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [-eps()])
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [0.0])
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [eps()])
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [1.0])
  @test gradcheck(x->3.0 * StatsFuns.logistic(x[1]), [5.0])
end

@testset "logit" begin
  @test gradcheck(x->5.0 * StatsFuns.logit(x[1]), [0.1])
  @test gradcheck(x->5.0 * StatsFuns.logit(x[1]), [0.3])
  @test gradcheck(x->5.0 * StatsFuns.logit(x[1]), [0.5])
  @test gradcheck(x->5.0 * StatsFuns.logit(x[1]), [0.7])
  @test gradcheck(x->5.0 * StatsFuns.logit(x[1]), [0.9])
end

function test_log1pexp(T, xs)
  y = T(4.3)
  for x in xs
    @test gradcheck(x->y * StatsFuns.log1pexp(x[1]), [x])
  end
end

@testset "log1pexp" begin
  @testset "Float64" begin
    @testset "x ∈ (-∞, 18.0)" begin
      test_log1pexp(Float64, [-1000.0, -50.0, -25.0, -10.0, 0.0, 10.0, 18.0 - eps()])
    end
    @testset "x ∈ [18.0, 33.3)" begin
      test_log1pexp(Float64, [18.0, 18.0 + eps(), 33.3 - eps()])
    end
    @testset "x ∈ [33.3, ∞)" begin
      test_log1pexp(Float64, [33.3, 33.3 + eps(), 100.0])
    end
  end
end

@testset "log1psq" begin
  rng = MersenneTwister(123456)
  @testset "Float64" begin
    for x in [-10.0, -5.0, -1.0, -eps(), 0.0, eps(), 1.0, 5.0, 10.0]
      @test gradcheck(x->5.1 * StatsFuns.log1psq(x[1]), [x])
    end
  end
end

@testset "logsumexp" begin
  rng = MersenneTwister(123456)
  @testset "Float64" begin
    @test gradtest(StatsFuns.logsumexp, randn(rng, 1))
    @test gradtest(StatsFuns.logsumexp, randn(rng, 1, 1))
    @test gradtest(StatsFuns.logsumexp, randn(rng, 3))
    @test gradtest(StatsFuns.logsumexp, randn(rng, 3, 4, 5))
  end
end

@testset "* sizing" begin
  @test size(Zygote.gradient((x, y)->sum(x * y), randn(1, 1), randn(1, 10))[1]) == (1, 1)
  @test size(Zygote.gradient((x, y)->sum(x * y), randn(1, 1), randn(1, 10))[2]) == (1, 10)
end


text:
# DIGEST-MD5 mechanism
function md5(v...)
    iob = IOBuffer()
    write(iob, v...)
    digest(MD_MD5, take!(iob))
end

function digmd5_decode_challenge(challenge::AbstractString)
    comps = Dict{AbstractString,AbstractString}()
    for comp in split(challenge, ',')
        (n,v) = split(comp, '='; limit=2)
        lv = length(v)
        if (lv > 1) && startswith(v, '"') && endswith(v, '"')
            v = SubString(v, 2, lv-1)
        end
        comps[n] = v
    end
    @debug("decoded challenge components", comps)
    comps
end

function digmd5_challenge_isvalid(challenge::Dict{AbstractString,AbstractString})
    !isempty(get(challenge, "realm", "")) &&
    !isempty(get(challenge, "nonce", "")) &&
    ("auth" == get(challenge, "qop", "")) &&
    ("utf-8" == get(challenge, "charset", "")) &&
    ("md5-sess" == get(challenge, "algorithm", ""))
end

const SASL_DIGESTMD5_CHARSET = "utf-8"
const SASL_DIGESTMD5_NC = "00000001"
const SASL_DIGESTMD5_MAXBUF = "65536"
const SASL_DIGESTMD5_METHOD = "AUTHENTICATE"

function digmd5_respond(username::AbstractString, password::AbstractString, digesturi::AbstractString, challenge::Dict{AbstractString,AbstractString})
    cnonce = randstring(40)
    realm = challenge["realm"]
    nonce = challenge["nonce"]
    qop = challenge["qop"]

    # if qop directive is "auth":
    # HA1=MD5(MD5(username:realm:password):nonce:cnonce)
    # HA2=MD5(method:digestURI)
    # response=MD5(HA1:nonce:nonceCount:clientNonce:qop:HA2)
    HA1_1 = md5(username, ':', realm, ':', password)
    HA1 = md5(HA1_1, ':', nonce, ':', cnonce)
    HA2 = md5(SASL_DIGESTMD5_METHOD, ':', digesturi)
    response = md5(bytes2hex(HA1), ':', nonce, ':', SASL_DIGESTMD5_NC, ':', cnonce, ':', qop, ':', bytes2hex(HA2))

    iob = IOBuffer()
    write(iob, "charset=",      SASL_DIGESTMD5_CHARSET, ',',
               "username=",     '"', username, '"',     ',',
               "realm=",        '"', realm, '"',        ',',
               "nonce=",        '"', nonce, '"',        ',',
               "nc=",           SASL_DIGESTMD5_NC,      ',',
               "cnonce=",       '"', cnonce, '"',       ',',
               "digest-uri=",   '"', digesturi, '"',    ',',
               "maxbuf=",       SASL_DIGESTMD5_MAXBUF,  ',',
               "response=",     bytes2hex(response),    ',',
               "qop=",          qop)
    String(take!(iob))
end

function testdigmd5(username::AbstractString, password::AbstractString, protocol::AbstractString, serverid::AbstractString, challenge::AbstractString)
    parts = digmd5_decode_challenge(challenge)
    digmd5_challenge_isvalid(parts) || error("invalid challenge params")
    digesturi = protocol * "/" * serverid
    digmd5_respond(username, password, digesturi, parts)
end

function digmd5_respond(token::TokenProto, protocol::AbstractString, serverid::AbstractString, challenge::Vector{UInt8})
    isempty(challenge) && error("empty challenge")
    parts = digmd5_decode_challenge(String(challenge))
    digmd5_challenge_isvalid(parts) || error("invalid challenge params")
    digesturi = protocol * "/" * serverid
    digmd5_respond(base64encode(token.identifier), base64encode(token.password), digesturi, parts)
end

buffer_sasl_reqhdr(channel::HadoopRpcChannel) = (channel.sent_call_id = channel.call_id = HRPC_CALL_ID_SASL; buffer_reqhdr(channel, HRPC_CALL_ID_SASL))

function buffer_sasl_message(channel::HadoopRpcChannel, state::Int32, auths::Vector{RpcSaslProto_SaslAuth}=RpcSaslProto_SaslAuth[], token::Vector{UInt8}=UInt8[])
    @debug("buffer SASL message", state, nauths=length(auths), token=!isempty(token))
    saslmsg = RpcSaslProto(version=0, state=state)
    isempty(auths) || setproperty!(saslmsg, :auths, auths)
    isempty(token) || setproperty!(saslmsg, :token, token)
    buffer_size_delimited(channel.iob, saslmsg)
end

function recv_sasl_message(channel::HadoopRpcChannel)
    @debug("recv SASL message")
    resp = RpcSaslProto()
    recv_rpc_message(channel, resp)
    @debug("received SASL message", resp)
    resp
end

token_alias(host::AbstractString, port::Integer) = "$(host):$(port)"
token_alias(channel::HadoopRpcChannel) = token_alias(channel.host, channel.port)
token_alias(nodeid::NodeIdProto) = token_alias(nodeid.host, nodeid.port)

function conditional_sasl_auth(channel::HadoopRpcChannel)
    # select appropriate token from ugi
    tok_alias = token_alias(channel)
    tokens = find_tokens(channel.ugi, alias=tok_alias)
    if isempty(tokens)
        @debug("no token available to authenticate, skipping.", to=tok_alias)
        return false
    end
    @debug("found tokens available to authenticate. authenticating...", ntokens=length(tokens), to=tok_alias)
    token = tokens[1]
    sasl_auth(channel, token)
    true
end

function sasl_auth(channel::HadoopRpcChannel, token::TokenProto)
    # send sasl handshake request
    begin_send(channel)
    buffer_handshake(channel, HRPC_AUTH_PROTOCOL_SASL)
    send_buffered(channel, false)

    # send sasl negotiate request
    begin_send(channel)
    buffer_sasl_reqhdr(channel)
    buffer_sasl_message(channel, RpcSaslProto_SaslState.NEGOTIATE)
    send_buffered(channel, true)
    resp = recv_sasl_message(channel)

    # check if TOKEN/DIGEST-MD5 is one of the supported methods
    nauths = length(resp.auths)
    @debug("server supports auth", nauths)
    idx_auth = 0
    for idx in 1:length(resp.auths)
        auth = resp.auths[idx]
        @debug("supported auth", idx, method=auth.method, mechanism=auth.mechanism)
        ((auth.method == "TOKEN") && (auth.mechanism == "DIGEST-MD5")) || continue
        idx_auth = idx
    end
    if idx_auth == 0
        @debug("no supported auth method found")
        throw(HadoopRpcException("no supported authentication method found"))
    end
    auth = resp.auths[idx_auth]

    # create response to the challenge
    protocol = auth.protocol
    serverId = auth.serverId
    challenge = auth.challenge
    @debug("auth", protocol, serverId, challenge)
    response = digmd5_respond(token, auth.protocol, auth.serverId, auth.challenge)
    @debug("response", response)

    # send response as a sasl initiate request
    respauth = RpcSaslProto_SaslAuth(method = auth.method,
                mechanism = auth.mechanism,
                protocol = auth.protocol,
                serverId = auth.serverId)

    begin_send(channel)
    buffer_sasl_reqhdr(channel)
    buffer_sasl_message(channel, RpcSaslProto_SaslState.INITIATE, [respauth], convert(Vector{UInt8}, codeunits(response)))
    send_buffered(channel, true)
    resp = recv_sasl_message(channel)

    # expect a success response
    if resp.state != RpcSaslProto_SaslState.SUCCESS
        @debug("error completing SASL auth", state=resp.state, expected=RpcSaslProto_SaslState.SUCCESS)
        throw(HadoopRpcException("error completing SASL auth. state: $(resp.state). expected state: $(RpcSaslProto_SaslState.SUCCESS)"))
    end
    nothing
end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

"""
    struct FloydWarshallState{T, U}

An [`AbstractPathState`](@ref) designed for Floyd-Warshall shortest-paths calculations.
"""
struct FloydWarshallState{T,U<:Integer} <: AbstractPathState
    dists::Matrix{T}
    parents::Matrix{U}
end

@doc_str """
floyd_warshall_shortest_paths(g, distmx=weights(g))
Use the [Floyd-Warshall algorithm](http://en.wikipedia.org/wiki/Floyd–Warshall_algorithm)
to compute the shortest paths between all pairs of vertices in graph `g` using an
optional distance matrix `distmx`. Return a [`FloydWarshallState`](@ref) with relevant
traversal information.

### Performance
Space complexity is on the order of ``\\mathcal{O}(|V|^2)``.
"""
function floyd_warshall_shortest_paths(
    g::AbstractGraph,
    distmx::AbstractMatrix{T} = weights(g)
) where T
    U = eltype(g)
    n_v = nv(g)
    dists = fill(typemax(T), (Int(n_v), Int(n_v)))
    parents = zeros(U, (Int(n_v), Int(n_v)))

    # fws = FloydWarshallState(Matrix{T}(), Matrix{Int}())
    for v in 1:n_v
        dists[v, v] = zero(T)
    end
    undirected = !is_directed(g)
    for e in edges(g)
        u = src(e)
        v = dst(e)

        d = distmx[u, v]

        dists[u, v] = min(d, dists[u, v])
        parents[u, v] = u
        if undirected
            dists[v, u] = min(d, dists[v, u])
            parents[v, u] = v
        end
    end
    for w in vertices(g), u in vertices(g), v in vertices(g)
        if dists[u, w] == typemax(T) || dists[w, v] == typemax(T)
            ans = typemax(T)
        else
            ans = dists[u, w] + dists[w, v]
        end
        if dists[u, v] > ans
            dists[u, v] = dists[u, w] + dists[w, v]
            parents[u, v] = parents[w, v]
        end
    end
    fws = FloydWarshallState(dists, parents)
    return fws
end

function enumerate_paths(s::FloydWarshallState{T,U}, v::Integer) where T where U<:Integer
    pathinfo = s.parents[v, :]
    paths = Vector{Vector{U}}()
    for i in 1:length(pathinfo)
        if (i == v) || (s.dists[v, i] == typemax(T))
            push!(paths, Vector{U}())
        else
            path = Vector{U}()
            currpathindex = i
            while currpathindex != 0
                push!(path, currpathindex)
                currpathindex = pathinfo[currpathindex]
            end
            push!(paths, reverse(path))
        end
    end
    return paths
end

enumerate_paths(s::FloydWarshallState) = [enumerate_paths(s, v) for v in 1:size(s.parents, 1)]
enumerate_paths(st::FloydWarshallState, s::Integer, d::Integer) = enumerate_paths(st, s)[d]


text:
using Documenter, TestPackage

makedocs(;
    modules=[TestPackage],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/csimal/TestPackage.jl/blob/{commit}{path}#L{line}",
    sitename="TestPackage.jl",
    authors="Cedric Simal",
    assets=String[],
)

deploydocs(;
    repo="github.com/csimal/TestPackage.jl",
)


text:
const TEST_TREES_NUMBERED = [
                             SINGLE_NODE_1, SINGLE_NODE_3, SINGLE_NODE_4,
                             LINEAR_TREE_1, LINEAR_TREE_2,
                             COMPLETE_BINARY_TREE_1, COMPLETE_BINARY_TREE_2
                            ]

const PREORDERS = [
                   [1],
                   [1],
                   [1],
                   [1,2,3,4],
                   [1,2,3,4],
                   [1,2,4,5,3,6,7],
                   [1,2,4,5,3,6,7]
                  ]

const POSTORDERS = [
                    [1],
                    [1],
                    [1],
                    [4,3,2,1],
                    [4,3,2,1],
                    [4,5,2,6,7,3,1],
                    [4,5,2,6,7,3,1]
                   ]

const LEVELORDERS = [
                     [1],
                     [1],
                     [1],
                     collect(1:4),
                     collect(1:4),
                     collect(1:7),
                     collect(1:7)
                    ]

@testset "preorder node iterator on single tree" for (t,r) in zip(TEST_TREES_NUMBERED, PREORDERS)
    @test r == map(t -> t.n, NodeIterator(t; order=PreOrder())) |> collect
end

@testset "postorder node iterator on single tree" for (t,r) in zip(TEST_TREES_NUMBERED, POSTORDERS)
    @test r == map(t -> t.n, NodeIterator(t; order=PostOrder())) |> collect
end

@testset "levelorder node iterator on single tree" for (t,r) in zip(TEST_TREES_NUMBERED, LEVELORDERS)
    @test r == map(t -> t.n, NodeIterator(t; order=LevelOrder())) |> collect
end

@testset "multiple tree handling" for t in TEST_TREES_NUMBERED
    for o in ORDERS
        @test collect(map(only, PredicateIterator(_node_predicate, (t,), false; order=o))) ==
                    collect(NodeIterator(t; order=o))
        @test collect(PredicateIterator(t -> t.n > 3, t; order=o)) ==
            collect(map(only, PredicateIterator(t -> only(t).n > 3, (t,), false; order=o))) ==
            vcat(map(unique, PredicateIterator(t -> t[1].n > 3, t, t; order=o))...)
    end
end

@testset "type iterator same types" for t in TEST_TREES_NUMBERED
    for o in ORDERS, type in TYPES
        @test collect(TypeIterator(type, t; order=o)) ==
        vcat(map(unique, TypeIterator(type, t, t; order=o))...) ==
        vcat(map(unique, TypeIterator((type, type), t, t; order=o))...)
    end
end

@testset "type iterator different types" for t in TEST_TREES_NUMBERED
    for o in ORDERS, type in TYPES
        @test collect(TypeIterator((type, type), t, t; order=o)) ==
        zip(collect(TypeIterator(type, t; order=o)),
            collect(TypeIterator(type, t; order=o))) |> collect
    end
end

@testset "empty iterator" for o in ORDERS
    @test collect(NodeIterator(; order=o, complete=true)) ==
    collect(NodeIterator(; order=o, complete=false)) == Union{}[]
    @test collect(LeafIterator(; order=o, complete=true)) ==
    collect(LeafIterator(; order=o, complete=false)) == Union{}[]
end

@testset "complete traversals" for o in ORDERS
    FULL_ORDER = Dict(PreOrder() => PREORDERS[end],
                        PostOrder() => POSTORDERS[end],
                        LevelOrder() => LEVELORDERS[end])
    for ts in powerset([T1, T2, T3, T4, T5])
        !isempty(ts) || continue
        res = NodeIterator(ts...; order=o, complete=true)
        if length(ts) == 1
            res = map(tuple, res)
        end
        res = collect(map(ts -> tuple([isnothing(t) ? t : t.n for t in ts]...), res))
        single_res = [collect(map(t -> t.n, NodeIterator(t; order=o))) for t in ts]
        idcs = intersect(FULL_ORDER[o], union(single_res...))
        @test all(irs -> all(r -> isnothing(r) || r == irs[1], irs[2]), zip(idcs, res))
    end
end

@testset "incomplete traversals" for o in ORDERS
    for ts in powerset([T1, T2, T3, T4, T5])
        !isempty(ts) || continue
        res = NodeIterator(ts...; order=o, complete=false)
        if length(ts) == 1
            res = map(tuple, res)
        end
        res = collect(map(ts -> tuple([isnothing(t) ? t : t.n for t in ts]...), res))
        @test all(r -> length(unique(r)) == 1, res)
        res = vcat([unique(r) for r in res]...)
        single_res = [collect(map(t -> t.n, NodeIterator(t; order=o))) for t in ts]
        @test intersect(single_res...) == res
    end
end


text:
#Asymptotic series for the branch 0 of the Lambert W function from Darko Veberic "Having Fun with Lambert W ( x ) Function" (2009)
function lambertWAsymptotic(x::Number)
    a = log(x)
    b = log(a)
    ia = 1 / a
    return a - b +
           b / a * (
        1 +
        ia * (
            0.5 * (-2 + b) +
            ia * (
                (6 + b * (-9 + b * 2)) / 6 +
                ia * (
                    (-12 + b * (36 + b * (-22 + b * 3))) / 12 +
                    ia / 60 * (60 + b * (-300 + b * (350 + b * (-125 + b * 12))))
                )
            )
        )
    )
end

# branch 0, valid for [0.3,7]
@inline lambertWRationalApproximation3(x) = x * (
        1 +
        x * (
            2.4450530707265568 +
            x *
            (1.3436642259582265 + x * (0.14844005539759195 + x * 0.0008047501729129999))
        )
    ) / (
        1 +
        x * (
            3.4447089864860025 +
            x * (3.2924898573719523 + x * (0.9164600188031222 + x * 0.05306864044833221))
        )
    )

# branch 0, valid for [-0.31,0.3]
@inline lambertWRationalApproximation1(x) =
    x * (
        1 +
        x * (
            5.931375839364438 +
            x * (11.392205505329132 + x * (7.338883399111118 + x * 0.6534490169919599))
        )
    ) / (
        1 +
        x * (
            6.931373689597704 +
            x * (16.82349461388016 + x * (16.43072324143226 + x * 5.115235195211697))
        )
    )

#Branch 0 of the Lambert W function
#algorithm follows Darko Veberic "Having Fun with Lambert W ( x ) Function" (2009)
#Fritsch's iteration on top of asymptotic series - machine accuracy for x>=8.70 (practical range for tanhsinh quadrature)
#Fritsch on top Rational approximation for x < 8.70 to reach machine accuracy in this region.
function lambertW(x::Float64)::Float64
    local wn
    if x < 0.14546954290661823
        wn = lambertWRationalApproximation1(x)
    elseif x < 8.706658967856612
        wn = lambertWRationalApproximation3(x)
    else
        wn = lambertWAsymptotic(x)
    end
    zn = log(x / wn) - wn
    qn = 2 * (1 + wn) * (1 + wn + 2 * zn / 3)
    en = (zn / (1 + wn)) * ((qn - zn) / (qn - 2 * zn))
    return wn * (1 + en)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAutoscalingV2beta1PodsMetricStatus <: SwaggerModel
    currentAverageValue::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApiResourceQuantity } # spec name: currentAverageValue
    metricName::Any # spec type: Union{ Nothing, String } # spec name: metricName
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector

    function IoK8sApiAutoscalingV2beta1PodsMetricStatus(;currentAverageValue=nothing, metricName=nothing, selector=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1PodsMetricStatus, Symbol("currentAverageValue"), currentAverageValue)
        setfield!(o, Symbol("currentAverageValue"), currentAverageValue)
        validate_property(IoK8sApiAutoscalingV2beta1PodsMetricStatus, Symbol("metricName"), metricName)
        setfield!(o, Symbol("metricName"), metricName)
        validate_property(IoK8sApiAutoscalingV2beta1PodsMetricStatus, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        o
    end
end # type IoK8sApiAutoscalingV2beta1PodsMetricStatus

const _property_map_IoK8sApiAutoscalingV2beta1PodsMetricStatus = Dict{Symbol,Symbol}(Symbol("currentAverageValue")=>Symbol("currentAverageValue"), Symbol("metricName")=>Symbol("metricName"), Symbol("selector")=>Symbol("selector"))
const _property_types_IoK8sApiAutoscalingV2beta1PodsMetricStatus = Dict{Symbol,String}(Symbol("currentAverageValue")=>"IoK8sApimachineryPkgApiResourceQuantity", Symbol("metricName")=>"String", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1PodsMetricStatus }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1PodsMetricStatus))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1PodsMetricStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1PodsMetricStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1PodsMetricStatus }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1PodsMetricStatus[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1PodsMetricStatus)
    (getproperty(o, Symbol("currentAverageValue")) === nothing) && (return false)
    (getproperty(o, Symbol("metricName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1PodsMetricStatus }, name::Symbol, val)
end


text:
"""
```julia
randUnitary(n::Int)
```
- Generates a n by n random Unitary matrix
- Equivalent to  `rand(Haar(2,n))`, see [`Haar`](@ref) and [`CUE`](@ref)
- For orthogonal matrices, use `randOrthogonal` or `rand(Haar(1,n))` instead
- The algorithm is disccused in [How to Generate Random Matrices from the Classical Compact Groups](http://www.ams.org/notices/200705/fea-mezzadri-web.pdf)

# Examples

Generate a 3 by 3 random Unitary matrix
```julia
randUnitary(3)

3×3 Matrix{ComplexF64}:
 -0.149398+0.0572715im  -0.0935861+0.629201im  -0.257255-0.709625im
  0.337035-0.342606im     -0.36366+0.599236im  -0.100838+0.517231im
  -0.17097+0.845103im   -0.0767105+0.313259im   0.247081+0.3025im
```
"""
function randUnitary(n::Int)
    A = randn(ComplexF64,n,n)
    Q,R = qr(A)
    Q,R = convert(Matrix,Q),convert(Matrix,R)
    Q = Q*Diagonal(sign.(diag(R)))
    return Q
end
"""
```julia
randOrthogonal(n::Int)
```
- Generates a `n` by `n` random Orthogonal matrix
- Equivalent to `rand(Haar(1,n))`, see [`Haar`](@ref) and [`COE`](@ref)
- For unitary matrices, use `randUnitary` or `rand(Haar(2,n))` instead
- The algorithm is disccused in [How to Generate Random Matrices from the Classical Compact Groups](http://www.ams.org/notices/200705/fea-mezzadri-web.pdf)

# Examples

Generates a 3 by 3 random Orthogonal matrix
```julia
randOrthogonal(3)

3×3 Matrix{Float64}:
 -0.875553  0.112448   0.469853
 -0.147915  0.863441  -0.482277
  0.459921  0.491757   0.739356
```
"""
function randOrthogonal(n::Int)
    A = randn(n,n)
    Q,R = qr(A)
    Q = Q*Diagonal(sign.(diag(R)))
    return Q
end
"""
```julia
Haar(beta,n)
```
- Uniform distribution on O(n) (`beta = 1`), U(n) (`beta = 2`)
- `beta`: 1 for Orthogonal, 2 for Unitary
- `n`: dimension

```julia
# Examples

# Generate a 100 by 100 random Unitary Matrix uniformly from U(n)
rand(Haar(2,100))

# Generate a 100 by 100 random Orthogonal Matrix uniformly from O(n)
rand(Haar(1,100))
```
"""
struct Haar <: ContinuousMatrixDistribution
    beta::Int
    n::Int
    Haar(beta,n) = beta in (1,2) ? new(beta,n) : error("Only take beta = 1 for orthogonal, beta = 2 for unitary")
end
Base.size(d::Haar) = d.n
rand(d::Haar) = d.beta == 1 ? randOrthogonal(d.n) : randUnitary(d.n)

######################################################
"""
```julia
COE
```
- `n` : Dimension
- Equivalent to [`Haar`](@ref)(1,n)
"""
const COE = n -> Haar(1,n)

"""
```julia
CUE
```
- `n` : Dimension
- Equivalent to [`Haar`](@ref)(2,n)
"""
const CUE = n -> Haar(2,n)

###################################################
"""
```julia
randPermutation(n; fix)
```
- `n` : dimension
- `fix` : a keyword argument, default is set to `fix = 0`. If `fix = x`, `randPermutation(n,x)` will have atleast `x` fixed points

# Examples

Generates a random 5 by 5 permutation matrix
```julia
randPermutation(5)

5×5 SparseArrays.SparseMatrixCSC{Int8, Int64} with 5 stored entries:
 ⋅  ⋅  ⋅  ⋅  1
 ⋅  1  ⋅  ⋅  ⋅
 1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  1  ⋅
```

Generates  Generates a random 10 by 10 permutation matrix with atleast 7 fix points
```julia
randPermutation(10, fix = 7)

10×10 SparseArrays.SparseMatrixCSC{Int8, Int64} with 10 stored entries:
 1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋮              ⋮
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1  ⋅  ⋅
 ⋅  ⋅  1  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  ⋅  1

```
"""
function randPermutation(n::Int; fix = 0::Int)
    O = sample(1:n,n-fix,replace=false) # randomly choosing flexible indices
    E = shuffle(O)
    M = spzeros(Int8,n,n)
    for i in 1:n-fix
      M[O[i],E[i]] = 1
    end
    for i in setdiff(1:n,O)
      M[i,i] = 1
    end
    return M
end

max_stars_repo_path: src/MatrixModels/Unitary Matrices.jl

max_issues_repo_path: src/MatrixModels/Unitary Matrices.jl

max_forks_repo_path: src/MatrixModels/Unitary Matrices.jl