text:
using TinyNERF, Images, Flux, CUDA, EllipsisNotation, Statistics, StatsBase, JSON
using Flux: @nograd, Zygote.ignore, params, glorot_uniform
CUDA.allowscalar(false)

transforms = JSON.parsefile("local/transforms_train.json")

W, H = 800, 800

const focal = W / (2*tan(transforms["camera_angle_x"]/2))
ṙ = tan(transforms["camera_angle_x"]/2) * 2 / W * 2/√12

poses = []
images = []

for frame in transforms["frames"]
    push!(poses, hcat(frame["transform_matrix"]...)')
    push!(images, float32.(channelview(load("./local/$(frame["file_path"]).png"))[1:3, :, :]))
end

poses = cat(poses..., dims=3)
images = cat(images..., dims=4) # channels × H × W × images

testpose = poses[:, :, 100]
testimg = gpu(images[:, :, :, 100])

poses = poses[:, :, 1:99]
images = images[:, :, :, 1:99]
colorview(RGB, cpu(testimg))

# function multiscale(imagetensor, n_scales, downsample, sigma=1)
#     out = [gaussian_pyramid(colorview(RGB, img), n_scales, downsample, sigma) for img in eachslice(imagetensor, dims=4)]
#     [cat(channelview.(getindex.(out, i))..., dims=4) for i in 1:n_scales+1]
# end

# multi_images = multiscale(images, 3, 2)

function posenc(x, L_embed)
    out = similar(x, (3 + 3 * 2 * L_embed, size(x)[2:end]...))
    out[1:3, ..] .= x
    out[4:end, ..] .= vcat((f.(2^i .* x) for i in 0:(L_embed - 1) for f in (sin, cos))...)

    return out
end
@nograd posenc

struct MLP{T}
    arch::T
end
Flux.@functor MLP

shifted_softplus(x) = softplus(x-1)

widened_sigmoid(x, ϵ) = (1+2ϵ)*sigmoid(x)-ϵ
widened_sigmoid(x) = widened_sigmoid(x, eltype(x)(0.001))

function MLP(W::Integer=256, size_pos=96, size_dir=3*(1+2*4))
    arch = Chain(
        SkipConnection(
            Chain(
                Dense(size_pos, W, relu),
                Dense(W, W, relu),
                Dense(W, W, relu),
                Dense(W, W, relu),
                Dense(W, W, relu)),
            (mx, x) -> vcat(mx, x)),
        Dense(W + size_pos, W, relu),
        Dense(W, W, relu),
        Dense(W, W, relu),
        Dense(W, W, identity),
        Dense(W-1+size_dir, W÷2, relu),
        Dense(W÷2, 3, widened_sigmoid))
    return MLP(arch)
end

function (m::MLP)(𝐱, normed_directions)
    x = reshape(𝐱, (size(𝐱, 1), :))
    x = m.arch[1:5](x)
    σ = reshape(shifted_softplus.(x[1, :]), (1, size(𝐱, 2), size(𝐱, 3)) )
    x = vcat(x[2:end, :], reshape(normed_directions, (size(normed_directions, 1), :) ))
    x = m.arch[6:7](x)
    rgb = reshape(x, (3, size(𝐱, 2), size(𝐱, 3)) )

    return rgb, σ
end

using BenchmarkTools
mlp = MLP()|>gpu

a = CUDA.rand(96, 1024, 128)
b = CUDA.rand(27, 1024, 128)

CUDA.@time mlp(a, b)

CUDA.@time begin
    gradient(mlp) do mlp
        rgb, s = Flux.Zygote.checkpointed(mlp, a, b)
        rgb2, s2 = Flux.Zygote.checkpointed(mlp, 2 .* a, b)

        sum(rgb .+ rgb2)
    end;
end

a = CUDA.rand(96, 1024, 20)
b = CUDA.rand(27, 1024, 20)

gs1 = gradient(Flux.params(mlp)) do
    sum(sum.(mlp(a, b)))
end

gs2 = gradient(Flux.params(mlp)) do
    sum(sum.(mlp(a, b)))
end

gs1 .+ 0.5 .* gs2

map!

@benchmark CUDA.@sync sum(sum.(mlp(a, b)))

@benchmark CUDA.@sync blocking=false gradient(mlp) do mlp
    sum(sum.(mlp(a, b)))
end

CUDA.@time gradient(mlp) do mlp
    sum(sum.(mlp(a, b)))
end

@btime CUDA.@sync blocking=false gradient(mlp) do mlp
    sum(sum.(mlp(a, b)))
end

"""
    get_rays(pixels, H, W, focal, cam2world)

Get origin and directions of rays through `pixels`.
`pixels` should be a vector of CartesianIndices.
"""
function get_rays(pixels, H, W, focal, cam2world)
    origin = similar(cam2world, 3)
    directions = similar(cam2world, (3, length(pixels)))

    directions[1, :] = (getindex.(pixels, 2) .- W/2 .- 1/2) ./ focal
    directions[2, :] = -(getindex.(pixels, 1) .- H/2 .- 1/2) ./ focal
    directions[3, :] .= -1
    directions = cam2world[1:3, 1:3] * directions

    origin = cam2world[1:3, 4]

    return origin, directions
end

gaussian_2d(x,y, σx, σy) = 1/(2π*σx*σy)*exp(-(((x-W/2)/σx)^2+((y-H/2)/σy)^2) / 2)

# using GLMakie
# begin
#     s = Scene()
#     cam3d!(s)
#     for p in eachslice(poses, dims=3)
#         weights = ProbabilityWeights(map((x)->gaussian_2d(x..., W/5, H/5), Tuple.(CartesianIndices((H, W))))[:])
#         pixels = CartesianIndices((H, W))[sample(1:H*W, weights, 10, replace=false)]
#         origin, directions = get_rays(pixels, H, W, focal, p)
#         # origin, directions = get_rays(CartesianIndices((H, W))[:], H, W, focal, p)
#         # origin, dirs = get_rays(H, W, focal, p)
#         arrows!(s,
#             [Point3(origin) for _ in 1:size(directions, 2)],
#             [Vec3(directions[:, i]) for i in 1:size(directions, 2)],
#             linewidth=0.01,
#             arrowsize=0,
#             lengthscale=0.5)
#     end
#     display(s)
# end

function get_t(directions, near, far, n_samples; randomized=false)
    t = reshape(range(near, stop=far, length=n_samples+1), (1, 1, :))
    if randomized
        t = t .+ rand(1, size(directions, 2), n_samples+1) .* ((far - near)/(n_samples+1))
    end
    return t
end

function get_t(directions::CuArray{T}, near, far, n_samples; randomized=false) where T
    t = CuArray(reshape(range(T(near), stop=far, length=n_samples+1), (1, 1, :)))
    if randomized
        return t .+ (CUDA.rand(1, size(directions, 2), n_samples+1) .* T((far - near)/(n_samples+1)))
    end
    return t
end

function cast(origin, directions, ṙ, t)
    T = eltype(t)

    t₀, t₁ = t[:, :, 1:(end-1)], t[:, :, 2:end]
    midpoint, halfwidth = (t₀ .+ t₁)./2, (t₁ .- t₀)./2

    μₜ = @. midpoint + (2*midpoint * (halfwidth^2)) / (3*(midpoint^2) + halfwidth^2)

    # https://github.com/JuliaGPU/CUDA.jl/issues/1044
    # varₜ = @. (halfwidth^2)/3 -  T(4/15) * ((halfwidth^4) * (12*(midpoint^2) - halfwidth^2) / (3*(midpoint^2) + halfwidth^2)^2)
    # varᵣ = @. ṙ^2 * ((midpoint^2)/4 + T(5/12)*(halfwidth^2) - T(4/15)*(halfwidth^4) / (3*midpoint^2 + halfwidth^2))

    varₜ = (halfwidth.^2)./3 .-  T(4/15) .* ((halfwidth.^4) .* (12 .*(midpoint.^2) .- halfwidth.^2) ./ (3 .*(midpoint.^2) + halfwidth.^2).^2)
    varᵣ = T(ṙ) .^ 2 .* ((midpoint .^ 2)./4 .+ T(5/12).*(halfwidth .^ 2) .- T(4/15).*(halfwidth .^ 4) ./ (3 .* (midpoint .^ 2)  .+ halfwidth .^ 2))

    μ = @. origin + μₜ * directions
    Σ_diag = varₜ .* (directions .^ 2) .+ varᵣ .* ((1 .- (directions .^ 2)) ./ sum(directions .^ 2, dims=1))

    normed_directions = similar(μ)
    normed_directions .= directions ./ sqrt.(sum(directions .^ 2, dims=1))

    return μ, Σ_diag, normed_directions
end

function IPE(μ, Σ_diag, degree=16)
    elt = eltype(Σ_diag)

    μᵧ = reduce(vcat, [2^i .* μ for i in 0:degree-1])
    Σᵧ_diag = reduce(vcat, [exp.((elt(-1/2) * 4^i) .* Σ_diag) for i in 0:degree-1])

    return [sin.(μᵧ) .* Σᵧ_diag; cos.(μᵧ) .* Σᵧ_diag]
end
@nograd IPE

function render(rgb, σ, t, directions)
    δ, midpoint, ϵ = ignore() do
        t₀, t₁ = t[:, :, 1:(end-1)], t[:, :, 2:end]
        midpoint, distance = (t₀ .+ t₁)./2, (t₁ .- t₀) # 1 × n_rays × (1 or n_samples)

        δ = distance .* sqrt.(sum(directions .^ 2, dims=1)) # 1 × n_rays × n_samples

        ϵ = eps(eltype(σ))
        return δ, midpoint, ϵ
    end

    α = 1 .- exp.(-σ .* δ) # 1 × n_rays × n_samples
    weights = α .* cat(gpu(ones(1, size(α, 2))), cumprod((1 .- α .+ ϵ)[:, :, 1:end-1], dims=3), dims=3) # 1 × n_rays × n_samples

    rgb = sum(weights .* rgb, dims=3) # 3 × n_rays × 1
    depth = sum(weights .* midpoint, dims=3) # 1 × n_rays × 1
    acc = sum(weights, dims=3) # 1 × n_rays × 1

    return rgb, depth, acc, weights
end

function train!(mlp, opt, images, poses; n_samples=128, batch_size=1024, n_iters=1000, i_plot=25)
    @assert batch_size <= H*W

    ps = params(mlp)
    local ℒ

    for i in 1:n_iters
        img_index = rand(1:size(images, 4))
        pose = gpu(poses[:, :, img_index])

        H, W = size(images)[[2, 3]]

        weights = ProbabilityWeights(map((x)->gaussian_2d(x..., W/6, H/6), Tuple.(CartesianIndices((H, W))))[:])

        pixels = CartesianIndices((H, W))[sample(1:H*W, weights, batch_size, replace=false)]

        # pixels = CartesianIndices((H, W))[rand(1:H*W, batch_size)]

        target = gpu(images[:, :, :, img_index][:, pixels])

        origin, directions = get_rays(pixels, H, W, focal, pose)
        t = get_t(directions, 2, 5, n_samples, randomized=true)
        μ, Σ_diag, normed_directions = cast(origin, directions, ṙ, t)

        gs = gradient(ps) do
            rgb, σ = mlp(IPE(μ, Σ_diag), posenc(normed_directions, 4))
            rgb, depth, acc, weights = render(rgb, σ, t, directions)
            ℒ = mean((rgb .- target).^2)
            return ℒ
        end

        Flux.update!(opt, ps, gs)

        if i%i_plot == 0
            rgbs = []
            for chunk in 0:ceil(Integer, (H*W)/batch_size)-1
                pixels = CartesianIndices((H, W))[chunk*batch_size+1:min(end, (chunk+1)*batch_size)]
                origin, directions = get_rays(pixels, H, W, focal, gpu(testpose))
                t = get_t(directions, 2, 5, n_samples, randomized=true)
                μ, Σ_diag, normed_directions = cast(origin, directions, ṙ, t)

                rgb, σ = mlp(IPE(μ, Σ_diag), posenc(normed_directions, 4))
                rgb, depth, acc, weights = render(rgb, σ, t, directions)
                push!(rgbs, rgb)
            end
            rgb = reshape(reduce(hcat, rgbs), (3, H, W))
            display(colorview(RGB, cpu(rgb)))
            testing_loss = mean((rgb .- testimg) .^ 2)
            PSNR = -10 * log10(testing_loss)
            println("training loss: $ℒ\ttesting loss: $testing_loss\tPSNR: $PSNR")
        end
    end
end

opt = ADAM(5e-3)
mlp = MLP()|>gpu

using BenchmarkTools

# multi_images[3]

# testimg = gpu(multiscale(reshape(cpu(testimg), (3, H, W, 1)), 3, 2)[1])

# ṙ = tan(transforms["camera_angle_x"]/2) * 2 / 200 * 2/√12

CUDA.@time train!(mlp, opt, images, poses, n_samples=128, batch_size=1024, n_iters=49, i_plot=50)
CUDA.@time train!(mlp, opt, images, poses; n_samples=64, batch_size=1024, n_iters=49, i_plot=50)

function test(mlp, opt, images, poses; n_samples=128, batch_size=1024, n_iters=1000)
    ps = params(mlp)
    for i in 1:n_iters
        img_index = rand(1:size(images, 4))
        pose = gpu(poses[:, :, img_index])

        H, W = size(images)[[2, 3]]

        weights = ProbabilityWeights(map((x)->gaussian_2d(x..., W/4, H/4), Tuple.(CartesianIndices((H, W))))[:])

        pixels = CartesianIndices((H, W))[sample(1:H*W, weights, batch_size, replace=false)]

        target = gpu(images[:, :, :, img_index][:, pixels])

        origin, directions = get_rays(pixels, H, W, focal, pose)
        t = get_t(directions, 2, 5, n_samples, randomized=true)
        μ, Σ_diag, normed_directions = cast(origin, directions, ṙ, t)

        gs = gradient(ps) do
            rgb, σ = mlp(IPE(μ, Σ_diag), posenc(normed_directions, 4))
            rgb, depth, acc, weights = render(rgb, σ, t, directions)
            ℒ = mean((rgb .- target).^2)
            return ℒ
        end

        ps_cpy, re = Flux.destructure(mlp)

        Flux.update!(opt, ps, gs)
        if any(cpu(isnan.(Flux.destructure(mlp)[1])))
            return re(ps_cpy), IPE(μ, Σ_diag), posenc(normed_directions, 4), t, directions, target
        end
        println(i)
    end
end

opt = ADAM(5e-4)
mlp = MLP()|>gpu

mlp, a, b, t, directions, target = test(mlp, opt, images, poses, n_samples=64, batch_size=2048, n_iters=1000)

any(cpu(isnan.(Flux.destructure(mlp)[1])))

target

gs = gradient(params(Flux.destructure(mlp)[1])) do
    rgb, σ = mlp(a, b)
    rgb, depth, acc, weights = render(rgb, σ, t, directions)
    ℒ = mean((rgb .- target).^2)
    return ℒ
end

gs.params[1]

extrema(cpu(abs.(gs.params[1])))

opt.state

any(cpu(isnan.(gs.params[1])))

c = nothing

for p in gs.params
    if any(cpu(isnan.(p)))
        c = p
        break
    end
end

collect(values(opt.state))[1][1]

c

findall(isnan, cpu(c))

extrema.(cpu.(gs.params))

[extrema(abs.(p)) for p in cpu.(gs.params)]

cpu.(gs.params)

params(mlp)

Flux.update!(opt, params(mlp), gs)

#/------------------------------------------------------------------------------------------------------\

img_index = rand(1:size(images, 4))
pose = gpu(poses[:, :, img_index])

weights = ProbabilityWeights(map((x)->gaussian_2d(x..., W/6, H/6), Tuple.(CartesianIndices((H, W))))[:])

pixels = CartesianIndices((H, W))[sample(1:H*W, weights, 1024, replace=false)]

target = gpu(images[:, :, :, img_index][:, pixels])

CUDA.@time origin, directions = get_rays(pixels, H, W, focal, pose)
CUDA.@time t = get_t(directions, 2, 5, 128, randomized=true)
CUDA.@time μ, Σ_diag, normed_directions = cast(origin, directions, ṙ, t)
@benchmark CUDA.@sync begin
    rgb, s = mlp(IPE(μ, Σ_diag), posenc(normed_directions, 4))
    # return sum(rgb)
    rgb, depth, acc, weights = render(rgb, s, t, directions)
    ℒ = mean((rgb .- target).^2)
end

a, b = IPE(μ, Σ_diag), posenc(normed_directions, 4)

CUDA.@time gs = Flux.gradient(mlp) do mlp
    rgb, s = mlp(IPE(μ, Σ_diag), posenc(normed_directions, 4))
    # return sum(rgb)
    rgb, depth, acc, weights = render(rgb, s, t, directions)
    ℒ = mean((rgb .- target).^2)
end
#\------------------------------------------------------------------------------------------------------/

using Printf
function lr_finder(mlp, opt, images, poses; start_lr=1e-6, end_lr=1, factor=2, n_samples=128, batch_size=1024)
    losses = []
    local ℒ
    ps = params(mlp)
    i = 1
    opt.eta = start_lr

    img_index = rand(1:size(images, 4))
    pose = gpu(poses[:, :, img_index])

    weights = ProbabilityWeights(map((x)->gaussian_2d(x..., W/4, H/4), Tuple.(CartesianIndices((H, W))))[:])

    pixels = CartesianIndices((H, W))[sample(1:H*W, weights, batch_size, replace=false)]

    target = gpu(images[:, :, :, img_index][:, pixels])

    origin, directions = get_rays(pixels, H, W, focal, pose)
    t = get_t(directions, 2, 5, n_samples, randomized=true)
    μ, Σ_diag, normed_directions = cast(origin, directions, ṙ, t)

    while opt.eta <= end_lr
        gs = gradient(ps) do
            rgb, σ = mlp(IPE(μ, Σ_diag), posenc(normed_directions, 4))
            rgb, depth, acc, weights = render(rgb, σ, t, directions)
            ℒ = mean((rgb .- target).^2)
            return ℒ
        end
        @printf "iteration %d\tη: %.2e\tloss: %.3f\n" i opt.eta ℒ
        if isnan(ℒ); break; end
        push!(losses, (opt.eta, ℒ))
        Flux.update!(opt, ps, gs)
        i += 1
        opt.eta *= factor
    end
    return losses
end


text:
using Benchmark
################################################################################
# Benchmark to see whether loop is faster than matrix operation
#
# | Row | Function     | Average   | Relative | Replications |
# |-----|--------------|-----------|----------|--------------|
# | 1   | "relu_loop"  | 0.0283602 | 1.0      | 30           |
# | 2   | "relu_whole" | 0.0611887 | 2.15755  | 30           |
#
# | Row | Function         | Average   | Relative | Replications |
# |-----|------------------|-----------|----------|--------------|
# | 1   | "relu_bwd_loop"  | 0.0314793 | 1.0      | 30           |
# | 2   | "relu_bwd_whole" | 0.0696191 | 2.21158  | 30           |
################################################################################

function relu_loop(input)
  input = copy(input)
  for i = 1:length(input)
    input[i] = max(input[i], 0)
  end
  return input
end

function relu_whole(input)
  input = copy(input)
  input = max(input, 0)
  return input
end

function relu_bwd_loop(input, grad)
  grad = copy(grad)
  for i = 1:length(input)
    if input[i] <= 0
      grad[i] = 0
    end
  end
  return grad
end
function relu_bwd_whole(input, grad)
  grad = copy(grad)
  grad .*= (input .> 0)
  return grad
end

input = rand(28, 28, 50, 128)
i1 = relu_loop(input)
i2 = relu_whole(input)
@assert all(abs(i1-i2) .< 1e-10)

grad = rand(size(input))

relu_loop() = relu_loop(input)
relu_whole() = relu_whole(input)
relu_bwd_loop() = relu_bwd_loop(i1, grad)
relu_bwd_whole() = relu_bwd_whole(i1, grad)

df = compare([relu_loop, relu_whole], 30)
println("$df")

df = compare([relu_bwd_loop, relu_bwd_whole], 30)
println("$df")

ext: jl


text:
# utility functions

tofloat(array::AbstractArray{T}) where T = map(Float32, array)
tofloat(array::AbstractArray{Float32}) = array

fft_frequencies(samplerate::Real, nfft::Int) = collect(linspace(0f0, samplerate / 2f0, (nfft >> 1) + 1))

"""returns the number of frames when the signal is partitioned into overlapping frames"""
nframes(length::Int, framesize::Int, hopsize::Int) = div(length - framesize, hopsize) + 1

"""Provides a view of a signal partitioned into overlapping frames, similar to MATLAB's buffer"""
struct FrameView{T<:AbstractVector} <: AbstractVector{Vector}
    original::T
    framesize::Int
    hopsize::Int
    nframes::Int

    function FrameView(original, framesize, hopsize)
        # n = noverlap is a problem - the algorithm will not terminate.
        #@boundscheck((0 < hopsize <= framesize), error("hopsize must be between zero and framesize"))
        new{T}(original, framesize, hopsize, nframes(length(original), framesize, hopsize))
    end
end

function Base.getindex(x::FrameView, i::Int)
    #@boundscheck((i >= 1 && i <= x.k), Base.throw_boundserror(x, i))
    slice(x.original, (i-1) * x.hopsize + (1:x.framesize))
end

Base.size(x::FrameView) = (x.nframes,)

function Base.iterate(x::FrameView, i::Int)
   if i > x.nframes
       nothing
   else
       (x[i], i+1)
   end
end

"""returns the DCT filters"""
function dct(nfilters::Int, ninput::Int)
    basis = Array(Float32, nfilters, ninput)
    samples = (1f0:2f0:2ninput) * π / 2ninput
    for i = 1:nfilters
        basis[i, :] = cos(i * samples)
    end

    basis *= sqrt(2f0/ninput)
    basis
end

"""returns the number of times that the value changes its sign"""
function zero_crossings(array::AbstractVector{T}, size::Int = length(array), offset::Int = 0) where {T<:Real}
    result = 0
    previous = 0
    for i in offset + (1:size)
        number = array[i]
        sgn = number == 0 ? 0 : number > 0 ? 1 : -1
        if sgn != previous
            result += 1
        end
        previous = sgn
    end
    result
end


text:
isapproxstr(s1::AbstractString, s2::AbstractString) =
    isequal(map(s->replace(s, r"\s|\n"=>""), String.((s1, s2)))...)


text:
@testset "operators" begin
    ### longitudinal magnetic field
    @test on_site_term(zeros(2), σᶻ) == spzeros(4, 4)

    @test on_site_term([0,1], σᶻ) == sparse(Diagonal([1,-1,1,-1]))
    @test on_site_term([1,0], σᶻ) == sparse(Diagonal([1,1,-1,-1]))
    @test on_site_term(ones(2), σᶻ) == sparse(Diagonal([2,0,0,-2]))

    @test on_site_term([2,0], σᶻ) == sparse(Diagonal([2,2,-2,-2]))
    @test on_site_term([2.1,0], σᶻ) == sparse(Diagonal([2.1,2.1,-2.1,-2.1]))

    ### spin-spin interaction term
    # sanity check
    @test spin_spin_interaction_term(zeros(2,2), σᶻ) == spzeros(4, 4)

    A = [
        0 1 0
        1 0 0
        0 0 0
        ]
    @test spin_spin_interaction_term(A, σᶻ) == sparse(Diagonal([1,1,-1,-1,-1,-1,1,1]))

    A = [
        0 2 0
        2 0 0
        0 0 0
        ]
    @test spin_spin_interaction_term(A, σᶻ) == sparse(Diagonal([2,2,-2,-2,-2,-2,2,2]))

    # 1d Ising model wiht the periodic boundary condition
    A = [
        0 1 1
        1 0 1
        1 1 0
        ]
    Ising1d = sparse(Diagonal([3,-1,-1,-1,-1,-1,-1,3]))
    @test spin_spin_interaction_term(A, σᶻ) == Ising1d
    @test spin_spin_interaction_term(sparse(A), σᶻ) == Ising1d

    # 2d Ising model
    A = [
     0  1  1  0
     1  0  0  1
     1  0  0  1
     0  1  1  0
    ]
    Ising2d = Diagonal([1,-1,-1,1, 1,-1,-1,1, 1,-1,-1,1, 1,-1,-1,1])
    Ising2d += Diagonal([1,1,1,1, -1,-1,-1,-1, -1,-1,-1,-1, 1,1,1,1])
    Ising2d += Diagonal([1,-1,1,-1, -1,1,-1,1, 1,-1,1,-1, -1,1,-1,1])
    Ising2d += Diagonal([1,1,-1,-1, 1,1,-1,-1, -1,-1,1,1, -1,-1,1,1])
    @test spin_spin_interaction_term(A, σᶻ) == sparse(Ising2d)
    @test spin_spin_interaction_term(A, σᶻ) == hamiltonian(zeros(4,4), zeros(4,4), A, zeros(4), zeros(4), zeros(4))

    # Assertion
    A = [
     0  1  1  0
     1  0  0  1
     1  0  0  1
     1  1  1  0
    ] # non-symmetrix matrix
    @test_throws AssertionError spin_spin_interaction_term(A, σᶻ)
end


text:
module Codim2

using LinearAlgebra
using Statistics: mean

using ..ArrayUtils: cast_container

using ..Continuations: AbstractContinuationProblem, AbstractProblemCache
import ..Continuations: get_prob_cache, get_u0, residual!, residual_jacobian!,
    residual, isindomain

# TimeKind trait:
using ..BifurcationsBase: timekind, Continuous, Discrete
import ..BifurcationsBase: TimeKind

# StateKind trait:
using ..BifurcationsBase: statekind, MutableState, ImmutableState
import ..BifurcationsBase: StateKind

# ContinuationKind trait:
using ..BifurcationsBase: contkind, SaddleNodeCont, HopfCont
using ..BifurcationsBase: ContinuationKind

import ..BifurcationsBase: BifurcationProblem

using Setfield: Lens
using StaticArrays: SVector

# TODO: Maybe rename (to be compatible with codim1lc):
#   problem.jl -> factories.jl
#   diffeq.jl  -> problem.jl

include("utils.jl")
include("augmented_systems.jl")
include("problem.jl")
include("diffeq.jl")
include("solver.jl")
include("analysis.jl")
include("switch.jl")
include("api.jl")
include("display.jl")
include("diagnosis.jl")

end  # module


text:
using Arblib

function sin_naive!(res, x)
    s, t, u = zero(x), zero(x), zero(x)
    tol = one(x)
    Arblib.mul_2exp!(tol, tol, -precision(tol))
    # @show tol
    k = 0
    while true
        Arblib.pow!(t, x, 2k + 1)
        Arblib.fac!(u, UInt(2k + 1))
        Arblib.div!(t, t, u)
        Arblib.abs!(u, t)

        if u ≤ tol
            # @show u
            Arblib.add_error!(s, u)
            break
        end
        if iseven(k)
            Arblib.add!(s, s, t)
        else
            Arblib.sub!(s, s, t)
        end
        k += 1
    end
    Arblib.set!(res, s)
end

global prec = 64
while true
    x = Arb("2016.1"; prec = prec)
    y = zero(x)
    y = sin_naive!(y, x)
    print("Using $(lpad(prec, 5)) bits, sin(x) = ")
    println(Arblib.string_nice(y, 10))
    y < zero(y) && break
    global prec *= 2
end


text:
using SparseIR
using Test

@testset "basis_set.consistency" begin
    beta = 2
    wmax = 5
    eps = 1e-5

    basis_f_org = FiniteTempBasis(fermion, beta, wmax, eps)
    basis_f, basis_b = finite_temp_bases(beta, wmax, eps, basis_f_org.sve_result)
    smpl_tau_f = TauSampling(basis_f)
    smpl_tau_b = TauSampling(basis_b)
    smpl_wn_f = MatsubaraSampling(basis_f)
    smpl_wn_b = MatsubaraSampling(basis_b)

    bs = FiniteTempBasisSet(beta, wmax, eps; sve_result=basis_f_org.sve_result)
    @test all(smpl_tau_f.sampling_points .== smpl_tau_b.sampling_points)
    @test all(bs.smpl_tau_f.matrix.a .== smpl_tau_f.matrix.a)
    @test all(bs.smpl_tau_b.matrix.a .== smpl_tau_b.matrix.a)
    @test all(bs.smpl_wn_f.matrix.a .== smpl_wn_f.matrix.a)
    @test all(bs.smpl_wn_b.matrix.a .== smpl_wn_b.matrix.a)

    bs2 = FiniteTempBasisSet(beta, wmax, eps; sve_result=bs.sve_result)
    @test all(bs2.smpl_tau_f.matrix.a .== smpl_tau_f.matrix.a)
end


text:
using RegressionTables: _escape
using Test
@testset "label transforms" begin
  # ampersand
  @test _escape(:ampersand)("X1 & X2") == "X1 \\& X2"

  # underscores
  @test _escape(:underscore2space)("long_var") == "long var"
  @test _escape(:underscore)("long_var") == "long\\_var"

  # latex
  @test _escape(:latex)("& % \$ # _ { }") == "\\& \\% \\\$ \\# \\_ \\{ \\}"

  # Dictionary
  @test _escape(Dict("bla" => "blaaaa"))("bla bla bla") == "blaaaa blaaaa blaaaa"
end


text:
using CompileBot

# using runtests:
snoop_bot(BotConfig("TestPackage0", tmin =0.0))


text:
# SPDX-License-Identifier: X11
# 2020-11-14

using Random

function geninput(bound::Integer, fn::AbstractString)
  X = collect(-bound:bound)
  Y = collect(-bound:bound)
  n = 2bound + 1

  println("Shuffling...")

  shuffle!(X)
  shuffle!(Y)

  println("Printing to $fn...")
  open(fn, "w") do f
    for i ∈ 1:n
      println(f, join([string(i, base=16),
                       string(X[i]),
                       string(Y[i])],
                      " "))
    end
  end
end

max_stars_repo_name: jthvai/fun-in-z

max_issues_repo_name: jthvai/fun-in-z

max_forks_repo_name: jthvai/fun-in-z


text:
# FASTA Writer
# ============

"""
    FASTAWriter(output::IO; width=70)

Create a data writer of the FASTA file format.

# Arguments
* `output`: data sink
* `width=70`: wrapping width of sequence characters
"""
type FASTAWriter{T<:IO} <: Bio.IO.AbstractWriter
    output::T
    # maximum sequence width (no limit when width ≤ 0)
    width::Int
end

function Bio.IO.stream(writer::FASTAWriter)
    return writer.output
end

function FASTAWriter(output::IO; width::Integer=70)
    return FASTAWriter(output, width)
end

function Base.write(writer::FASTAWriter, seqrec::FASTASeqRecord)
    output = writer.output
    n = 0

    # header
    n += write(output, '>', seqrec.name)
    if !isempty(seqrec.metadata.description)
        n += write(output, ' ', seqrec.metadata.description)
    end
    n += write(output, '\n')

    # sequence
    w = writer.width
    for x in seqrec.seq
        if writer.width > 0 && w == 0
            n += write(output, '\n')
            w = writer.width
        end
        n += write(output, Char(x))
        w -= 1
    end
    n += write(output, '\n')

    return n
end


text:
mutable struct SumFun <: Function# OscLaurent{D<:SPeriodicLine,R}?
    funs::Vector{Fun}
end

evaluate(f::SumFun,x) = sum(map(y -> y(x),f.funs))

SumFun(f::Fun{T}) where T = SumFun([f])
SumFun(f::Fun{T}) where T<: SumSpace = SumFun(components(f))
SumFun(f::Fun{T}) where T<: PiecewiseSpace = SumFun(components(f))

conj(s::SumFun) = SumFun(map(conj,s.funs))
conj(s::Array{T}) where T<:SumFun = map(conj,s)

#sum(s::SumFun) = sum(map(sum,s.funs))
sum(s::SumFun) = s.funs == [] ? 0.0 : sum(map(sum,s.funs))
sum(s::Array{T}) where T<:SumFun = map(sum,s)

function combine!(s::SumFun)
    s.funs = combine!(s.funs)
    s
end
combine!(s::Vector{SumFun}) = map(combine!,s)
combine!(s::Array{SumFun}) = map(combine!,s)
components(f::Fun) = f

function combine!(f::Vector{T}) where T <: Fun
    f = vcat(map(components,f))
    j = 0
    @inbounds while j < length(f)
        j += 1
        if (length(f[j].coefficients) <= 4 && sum(abs.(f[j].coefficients)) < 1e-15 && length(f) > 2) || maximum(abs.(f[j].coefficients)) < 1e-15
            deleteat!(f,j)
            j -= 1
        end
    end
    i = 0
    @inbounds while i < length(f)
        i += 1
        j = i
        while j < length(f)
            j += 1
            if spacescompatible(f[i],f[j])
                f[i] += f[j]
                deleteat!(f,j)
                j -= 1
            end
        end
    end
    f
end

+(f::SumFun,g::SumFun) = SumFun(vcat(g.funs,f.funs))
-(f::SumFun,g::SumFun) = SumFun(vcat(-g.funs,f.funs))

+(f::SumFun,g::Fun) = f + SumFun(g)
+(g::Fun,f::SumFun) = f + SumFun(g)

-(f::SumFun,g::Fun) = f - SumFun(g)
-(g::Fun,f::SumFun) = -f + SumFun(g)

-(f::SumFun) = SumFun(map(y -> -y,f.funs))

*(f::Fun,g::SumFun) = SumFun(map(y -> f*y,g.funs))
*(f::SumFun,g::Fun) = SumFun(map(y -> g*y,f.funs))
*(f::SumFun,g::Number) = SumFun(map(y -> g*y,f.funs))
*(f::Number,g::SumFun) = SumFun(map(y -> f*y,g.funs))
*(f::SumFun,g::Array) = map(x -> f*x,g)

*(f::SumFun,g::SumFun) = combine!(+(map(y -> y*g,f.funs)...))

*(c::Float64,g::SumFun) = SumFun(map(y -> c*y,g.funs))
*(c::Complex{Float64},g::SumFun) = SumFun(map(y -> c*y,g.funs))

(f::SumFun)(x) = evaluate(f,x)
+(f::SumFun) = f

(f::Array)(x) = map(y -> evaluate(y,x),f)
chop!(f::SumFun) = SumFun(map(chop!,f.funs))
chop!(f::Array) = map(chop!,f)

chop!(f::SumFun,tol::Float64) = SumFun(map(x -> chop!(x,tol),f.funs))
chop!(f::Array,tol::Float64) = map(x -> chop!(x,tol),f)
copy(f::SumFun) = SumFun(copy(f.funs))


text:
using Test

@testset "metric" begin
	include("metric.jl")
end

@testset "kernel" begin
	include("kernel.jl")
end

@testset "negloglik" begin
	include("negloglik.jl")
end

@testset "utils" begin
	include("utils.jl")
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
using LinearAlgebra, MessyTimeSeries, Test;

"""
    read_test_input(filepath::String)

Read input data necessary to run the test for the Kalman routines. It does not use external dependencies to read input files.
"""
function read_test_input(filepath::String)

    # Load CSV into Array{SubString{String},1}
    data_str = split(read(open("$filepath.txt"), String), "\n");
    deleteat!(data_str, findall(x->x=="", data_str));

    # Return output
    data = eval(Meta.parse(data_str[1]));
    return data;
end


text:
function readrawim_Mono8(file::Array,width,height;n=length(file),T=Float64)
    n == 0 && error("No files to read.")
    imgs = Array{T}(undef,height,width,n)
    @inbounds for i in eachindex(file)
        imgs[:,:,i] = readrawim_Mono8(file[i],width,height)
    end
    imgs
end
function readrawim_Mono12Packed(file::Array,width,height;n=length(file),T=Float64)
    n == 0 && error("No files to read.")
    imgs = Array{T}(undef,height,width,n)
    raw = Vector{UInt8}(undef,Int(width*height*1.5))
    npack = Int(length(raw)/3)
    img = Vector{UInt8}(undef,4npack)
    @inbounds for i in eachindex(file)
        imgs[:,:,i] = readrawim_Mono12Packed(file[i],width,height;raw,npack,img)
    end
    imgs
end

readrawim_Mono8(file::String,width,height) = Mmap.mmap(file,Matrix{UInt8},(width,height))'
function readrawim_Mono12Packed(file::String,width,height;raw = Vector{UInt8}(undef,Int(width*height*1.5)),
                                npack = Int(length(raw)/3), img = Vector{UInt8}(undef,4npack))
    # Unpack Mono12Packed 3 bytes to 2 UInt16
    # The two 12bits pixels | A | B | are packed in | A11 A10 A9 A8 A7 A6 A5 A4 | B3 B2 B1 B0 A3 A2 A1 A0 | B11 B10 B9 B8 B7 B6 B5 B4 | 3 bytes
    read!(file,raw)
    @turbo unroll=1 thread=true for i in 0:npack-1
        img[1+4i] = raw[2+3i] << 4 # second byte shift 4 bits left | A3 A2 A1 A0 0 0 0 0 |
        img[2+4i] = raw[1+3i]
        img[3+4i] = raw[2+3i]
        img[4+4i] = raw[3+3i]
    end
    img = reinterpret(UInt16,img)
    @turbo unroll=1 thread=true for i in 1:2npack
        img[i] >>>= 4 # shift UInt16s 4 bits right
    end
    reshape(img,width,height)'
end


text:
using ADFPEPS
using ADFPEPS:double_ipeps_energy
using CUDA
using Random

CUDA.allowscalar(false)
Random.seed!(100)
model = hop_pair(1.0,1.0)
folder = "E:/1 - research/4.9 - AutoDiff/data/ADFPEPS/"
ipeps, key = init_ipeps(model; Ni = 2, Nj = 2, atype = Array, folder = folder, D=2, χ=20, tol=1e-10, maxiter=10)
# folder, model, Ni, Nj, atype, D, χ, tol, maxiter = key
# double_ipeps_energy(atype(ipeps), model; Ni=Ni,Nj=Nj,χ=χ,maxiter=10,infolder=folder,outfolder=folder)
res = optimiseipeps(ipeps, key; f_tol = 1e-10, opiter = 100, verbose = true)


text:
using JuMP, Mosek, SCS
import MathProgComplex

function run_SDPpb(cst, optsense)
    # m = Model(solver=MosekSolver(LOG=0))
    m = Model(solver=SCSSolver(verbose=0))

    @variable(m, Xs1[1:3,1:3], SDP)
    @variable(m, Xs2[1:1,1:1], SDP)
    @variable(m, Xs3[1:1,1:1], SDP)

    @constraint(m, Xs1[1,1] == 1)
    @constraint(m, Xs2[1,1] + Xs1[2,2] + Xs1[3,3] == 4)
    @constraint(m, Xs3[1,1] - Xs1[2,1] - Xs1[3,1] == 0)

    # Find upper bound
    @objective(m, optsense, -Xs1[1,3]+cst)
    solve(m)
    return getobjectivevalue(m)
end

@show run_SDPpb(0., :Min)
@show run_SDPpb(0.5, :Min)

@show run_SDPpb(0., :Max)
@show run_SDPpb(0.5, :Max)


text:
# 1D regression.

import Random
import PyPlot

using LinearAlgebra
#import Interpolations
#using Revise

import RKHSRegularization # https://github.com/RoyCCWang/RKHSRegularization

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

PyPlot.close("all")

Random.seed!(25)

fig_num = 1

##### regression
θ = RKHSRegularization.Spline34KernelType(0.04)

# user inputs.
σ² = 1e-3
N = 15

x_range = LinRange(-24.0, 45.0, N)
Nq = 1000
#xq_range = LinRange(-30.0, 50.0, Nq)
Nq = 1000
xq_range = LinRange(x_range[1], x_range[end], Nq)

#DelimitedFiles.writedlm( "query.txt",  xq_range, '\n')

# if 0 and 1 are included, we have posdef error, or rank = N - 2.
X = collect( [x_range[n]] for n = 1:N )

f = xx->sinc(4*xx)*xx^3
y = f.(x_range)

# check posdef.
K = RKHSRegularization.constructkernelmatrix(X, θ)
println("rank(K) = ", rank(K))

println("isposdef = ", isposdef(K))

# fit RKHS.
η = RKHSRegularization.RKHSProblemType( zeros(Float64,length(X)),
                     X,
                     θ,
                     σ²)
RKHSRegularization.fitRKHS!(η, y)

# query.
xq = collect( [xq_range[n]] for n = 1:Nq )

f_xq = f.(xq_range) # generating function.

yq = Vector{Float64}(undef, Nq)
RKHSRegularization.query!(yq,xq,η)

# Visualize regression result.
PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(X, y, ".", label = "observed")
PyPlot.plot(xq, yq, label = "GP mean function")

title_string = "1D GP regression"
PyPlot.title(title_string)
PyPlot.legend()
##### end of regression

max_stars_repo_name: RoyCCWang/RKHSRegularization

max_issues_repo_name: RoyCCWang/RKHSRegularization

max_forks_repo_name: RoyCCWang/RKHSRegularization


text:
@testset "detrend()" begin
    # use detrend on polynomials and see if results are close to 0
    # (add 1. and see if results ≈ 1. as a workaround)
    for n ∈ 2:4
        @test begin
            x = collect(Float64, 1:20)
            y = x.^n
            all(yᵢ + 1.0 ≈ 1.0 for yᵢ ∈ nsi.detrend(x, y, n)) == true
        end
    end
    @test begin
        x = collect(Float64, 1:20)
        y = @. 0.1 + 2x - 0.01x^2
        all(yᵢ + 1.0 ≈ 1.0 for yᵢ ∈ nsi.detrend(x, y, 2)) == true
    end
end

@testset "NoiseSample" begin

    # test construction

    let
        ns1 = NoiseSample{Int}([1, 2, 3], [4, 5, 6])
        ns2 = NoiseSample([1, 2, 3], [4, 5, 6])
        ns3 = NoiseSample(Curve([1, 2, 3], [4, 5, 6]))
        ns4 = NoiseSample([4, 5, 6])
        for ns in (ns1, ns2, ns3, ns4)
            @test length(ns) == 3
            @test eltype(ns) == Int
            @test ns.x == [1, 2, 3]
            @test ns.y == [-1, 0, 1]
        end
    end

    # test basic operation

    noise = NoiseSample([1.0, 2.0], [4.0, 4.0])
    @test length(noise) == 2
    @test noise[end][1] == 2.0
    @test noise[end][2] + 1 ≈ 1.0
    @test [noise; noise] == NoiseSample([1.0, 2.0, 1.0, 2.0], [4.0, 4.0, 4.0, 4.0])

    # x and y need to have same length
    @test_throws ArgumentError NoiseSample([1.0], [4.0, 4.0])
    @test_throws ArgumentError NoiseSample([1.0, 2.0], [4.0])
end;

@testset "GaussianNoiseModel" begin
    @test eltype(GaussianNoiseModel(1)) == Int
    @test eltype(GaussianNoiseModel(1.)) == Float64
    @test_throws MethodError GaussianNoiseModel(1 + 1im)
end

@testset "MvGaussianNoiseModel" begin
    @test eltype(MvGaussianNoiseModel(1, 1, 1)) == Int
    @test eltype(MvGaussianNoiseModel(1., 1., 1.)) == Float64
    @test_throws MethodError MvGaussianNoiseModel(1 + 1im, 1, 1)
end

@testset "gauss_kernel()" begin
    @test nsi.gauss_kernel(0, [2, 1]) == 4.0
    @test nsi.gauss_kernel(2.0, [2, 2]) == 4.0 * exp(-0.5)
    @test nsi.gauss_kernel(4.0, [2, 1]) == 4.0 * exp(-8)
    @test nsi.gauss_kernel(4.0, [0, 1]) == 0.0
end

@testset "get_cov()" begin
    let f = 1.0000001
        @test_throws ArgumentError nsi.get_cov(1.0, -1, 2.0, 1.0)
        @test nsi.get_cov(1.0, 3, 2.0, 1.0) == 4 .* [f         exp(-0.5) exp(-2)  ;
                                                     exp(-0.5) f         exp(-0.5);
                                                     exp(-2)   exp(-0.5) f         ]
    end
end

@testset "add_noise()" begin
    let # use MvGaussianNoiseModel
        c = Curve([1., 2., 3.], [4., 5., 6.])
        nm = MvGaussianNoiseModel(1., 3., 5.)
        noise = add_noise(c, nm)
        @test length(noise) == 3
        @test noise.y |> first |> x -> x.particles |> length == 10_000
        noise = add_noise(c, nm, 1000)
        @test noise.y |> first |> x -> x.particles |> length == 1000
    end
    let # use GaussianNoiseModel
        c = Curve([1., 2., 3.], [4., 5., 6.])
        nm = GaussianNoiseModel(1.)
        noise = add_noise(c, nm)
        @test length(noise) == 3
        @test noise.y |> first |> x -> x.particles |> length == 10_000
        noise = add_noise(c, nm, 1000)
        @test noise.y |> first |> x -> x.particles |> length == 1000
    end
end

@testset "generate_noise() (correlated)" begin
    let
        nm = MvGaussianNoiseModel(1.0, 3.0, 5.0)
        noise = nsi.generate_noise(nm, 100, 2000)
        @test length(noise) == 100
        @test noise |> first |> x -> x.particles |> length == 2000
        ns = NoiseSample(nsi.get_draw.(1, noise))
        noise2 = nsi.generate_noise(ns, 100)
        @test noise2 |> first |> x -> x.particles |> length == 100
    end
end

@testset "generate_noise() (uncorrelated)" begin
    let
        nm = GaussianNoiseModel(1.0)
        noise = nsi.generate_noise(nm, 100, 2000)
        @test length(noise) == 100
        @test noise |> first |> x -> x.particles |> length == 2000
    end
end

@testset "fit_noise()" begin
    # this also tests estimate_autocov()
    let
        seed!(2)
        samples = 100
        len = 2_000
        for (δxᵢ, αᵢ, λᵢ) ∈ [(1.0, 3.0, 5.0),
                             (0.5, 2.0, 0.5),
                             (2.0, 0.5, 3.0)]
            nm = MvGaussianNoiseModel(δxᵢ, αᵢ, λᵢ)
            noise = nsi.generate_noise(nm, len, samples)
            A = Float64[]
            Λ = Float64[]
            for i in 1:samples
                x = δxᵢ:δxᵢ:(len * δxᵢ) |> collect
                ns = NoiseSample(x, nsi.get_draw.(i, noise))
                noise_param = fit_noise(ns)
                push!(A, noise_param.α)
                push!(Λ, noise_param.λ)
            end
            # To approach the true values, we need to have very long
            # noise samples, which gets computationally expensive
            # quickly. As a workaround, shorter noise samples
            # are used here, and the fitted parameters are averaged
            # over several trials.
            # When using a noise sample length of 2000 and averaging
            # over 100 trials, the fitted parameters should reproduce
            # the input parameters with 5% accuray:
            @test mean(A) ≈ αᵢ atol=5e-2
            @test mean(Λ) ≈ λᵢ atol=5e-2
        end
    end
end


text:
#=

  The BUGS Book
  Example 3.4.1 Three coins (page 45)
  """
  Suppose I have 3 coins in my pocket. The coins may be either fair, biased 3:1 in
  favour of heads, or 3:1 in favour of tails, but I don't know how many of
  each type there are among the 3 coins. I randomly select 1 coin att toss it once,
  observing a head. What is the posterior distribution of the probability of a head?
  """

  Note that JAGS has the somewhat weird structure of dcat etc.

  ~/jags/bugs_book_3_4_1.R
  """
  data {
    y <- 1
  }

  model {

    y ~ dbern(theta.true)
    theta.true <- theta[coin]
    coin ~ dcat(p[])
    for(i in 1:3) {
      p[i] <- 1/3
      theta[i] <- 0.25*i
      coin.prob[i] <- equals(coin, i)
    }
  }

  Output:
               Mean     SD  Naive SE Time-series SE
coin         2.3320 0.7458 0.0026369      0.0026369
coin.prob[1] 0.1672 0.3732 0.0013194      0.0013194
coin.prob[2] 0.3335 0.4715 0.0016669      0.0016669
coin.prob[3] 0.4993 0.5000 0.0017678      0.0017590
theta.true   0.5830 0.1865 0.0006592      0.0006592
  """

  Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

          coin    2.3303    0.7518     0.0075    0.0101   4837.8060    1.0000      321.8124
      theta[1]    0.2500    0.0000     0.0000    0.0000         NaN       NaN           NaN
      theta[2]    0.5000    0.0000     0.0000    0.0000         NaN       NaN           NaN
      theta[3]    0.7500    0.0000     0.0000    0.0000         NaN       NaN           NaN
  coin_prob[1]    0.1720    0.3774     0.0038    0.0049   4883.5182    1.0000      324.8532
  coin_prob[2]    0.3257    0.4687     0.0047    0.0064   4998.4834    0.9999      332.5007
  coin_prob[3]    0.5023    0.5000     0.0050    0.0069   4914.2258    1.0000      326.8959
    theta_true    0.5826    0.1880     0.0019    0.0025   4837.8060    1.0000      321.8124
             y    1.0000    0.0000     0.0000    0.0000         NaN       NaN           NaN

  Distributions of variable coin (num:0)
  3.00000 =>    4996  (0.499600)
  2.00000 =>    3355  (0.335500)
  1.00000 =>    1649  (0.164900)

  Cf ~/webppl/bugs_book_3_4_1.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function bugs_book_3_4_1()
    n = 3
    coin ~ DiscreteUniform(1,3) #  Select one coin. Note: JAGS has 1-based indexing.
    p = 1/3
    coin_prob = tzeros(3)
    theta = tzeros(3)
    for i in 1:n
        theta[i] ~ Dirac(0.25*i)
        coin_prob[i] ~ Dirac(coin == i) # Is this is the coin
    end

    theta_true ~ Dirac(theta[coin])
    y ~ flip(theta_true)

    true ~ Dirac(y == 1) #  We observe a head.

end

model = bugs_book_3_4_1()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 100_000)
# chns = sample(model, PG(15), 10_000)
chns = sample(model, SMC(10_000), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, NUTS(), 10_000)
# chns = sample(model, HMC(0.1,5), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :coin)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
function bin(x::Integer,pad::Int,id::Int)
    return (x>>(pad-id))&1
end

function bin(x::Integer,pad::Int,ids::Vector{Int})
    ret = 0
    subid = 1;subpad = length(ids)
    for id in ids
        ret+=(bin(x,pad,id)<<(subpad-subid))
        subid+=1
    end
    return ret
end

function flip(x::Integer,pad::Int,id::Int)
    return x$(1<<(pad-id))
end

function flip(x::Integer,pad::Int,ids::Vector{Int},assign::Int)
    ret = x
    n = length(ids)
    for id = 1:n
        if bin(assign,n,id)==1
            ret = flip(ret,pad,id)
        end
    end
    return ret
end

function set(x::Integer,pad::Int,ids::Vector{Int},assign::Integer)
    ret = x
    n = length(ids)
    for i = 1:n
        if bin(assign,n,i) != bin(x,n,ids[i])
            ret = flip(ret,pad,ids[i])
        end
    end
    return ret
end

function process!{N}(unit::AbstractGateUnit{N},input::AbstractSparseArray)
    len = length(unit.pos)-1
    eigens = unit.gate()
    ret = spzeros(Complex,size(input)...)
    for i in rowvals(input)
        subidx = bin(i-1,len,unit.pos[2:end])
        rett = spzeros(Complex,size(input)...)
        for j = 0:(2^len-1)
            idx = set(i-1,N,unit.pos[2:end],j)
            rett[idx+1] = input[i]*eigens[subidx+1][j+1]
        end
        ret += rett
    end
    input[:] = ret
    return ret
end

function process!{N}(circuit::Circuit{N},input::AbstractSparseArray)
    ret = deepcopy(input)
    for unit in circuit.gates
        if N==bitnum(unit)
            ret = unit.gate(ret)
        else
            ret = process!(unit,ret)
        end
    end
    return ret
end

export process!


text:
using Test

include("perfect-numbers.jl")

@testset "Perfect numbers" begin

    @testset "Smallest perfect number is classified correctly" begin
        @test isperfect(6)
    end

    @testset "Medium perfect number is classified correctly" begin
        @test isperfect(28)
    end

    @testset "Large perfect number is classified correctly" begin
        @test isperfect(33550336)
    end
end

@testset "Abundant numbers" begin

    @testset "Smallest abundant number is classified correctly" begin
        @test isabundant(12)
    end

    @testset "Medium abundant number is classified correctly" begin
        @test isabundant(30)
    end

    @testset "Large abundant number is classified correctly" begin
        @test isabundant(33550335)
    end
end

@testset "Deficient numbers" begin

    @testset "Smallest prime deficient number is classified correctly" begin
        @test isdeficient(2)
    end

    @testset "Smallest non-prime deficient number is classified correctly" begin
        @test isdeficient(4)
    end

    @testset "Medium deficient number is classified correctly" begin
        @test isdeficient(32)
    end

    @testset "Large deficient number is classified correctly" begin
        @test isdeficient(33550337)
    end

    @testset "Edge case (no factors other than itself) is classified correctly" begin
        @test isdeficient(1)
    end
end

@testset "Invalid inputs" begin

    @testset "Zero is rejected (not a natural number)" begin
        @test_throws DomainError isdeficient(0)
        @test_throws DomainError isperfect(0)
        @test_throws DomainError isabundant(0)
    end

    @testset "Negative integer is rejected (not a natural number)" begin
        @test_throws DomainError isdeficient(-1)
        @test_throws DomainError isperfect(-1)
        @test_throws DomainError isabundant(-1)
    end
end


text:
using ParetoSmooth
using Test

using AxisKeys
using NamedDims
using Statistics

import RData

@testset "ParetoSmooth.jl" begin

    # Tests must be kept in this order
    include("BasicTests.jl")
    include("TuringTests.jl")
    include("ComparisonTests.jl")

end


text:
using Plots, Test

@testset "Shorthands" begin
    @testset "Set Lims" begin
        p = plot(rand(10))

        xlims!((1, 20))
        @test xlims(p) == (1, 20)

        xlims!(p, (1, 21))
        @test xlims(p) == (1, 21)

        ylims!((-1, 1))
        @test ylims(p) == (-1, 1)

        ylims!(p, (-2, 2))
        @test ylims(p) == (-2, 2)

        zlims!((-1, 1))
        @test zlims(p) == (-1, 1)

        zlims!(p, (-2, 2))
        @test zlims(p) == (-2, 2)

        xlims!(-1, 11)
        @test xlims(p) == (-1, 11)

        xlims!(p, -2, 12)
        @test xlims(p) == (-2, 12)

        ylims!((-10, 10))
        @test ylims(p) == (-10, 10)

        ylims!(p, (-11, 9))
        @test ylims(p) == (-11, 9)

        zlims!((-10, 10))
        @test zlims(p) == (-10, 10)

        zlims!(p, (-9, 8))
        @test zlims(p) == (-9, 8)
    end

    @testset "Set Title / Labels" begin
        p = plot()
        title!(p, "Foo")
        sp = p[1]
        @test sp[:title] == "Foo"
        xlabel!(p, "xlabel")
        @test sp[:xaxis][:guide] == "xlabel"
        ylabel!(p, "ylabel")
        @test sp[:yaxis][:guide] == "ylabel"
    end

    @testset "Misc" begin
        p = plot()
        sp = p[1]

        xflip!(p)
        @test sp[:xaxis][:flip]

        yflip!(p)
        @test sp[:yaxis][:flip]

        xgrid!(p, true)
        @test sp[:xaxis][:grid]

        xgrid!(p, false)
        @test !sp[:xaxis][:grid]

        ygrid!(p, true)
        @test sp[:yaxis][:grid]

        ygrid!(p, false)
        @test !sp[:yaxis][:grid]

        ann = [(7, 3, "(7,3)"), (3, 7, text("hey", 14, :left, :top, :green))]
        annotate!(p, ann)
        annotate!(p, ann...)

        xaxis!(p, true)
        @test sp[:xaxis][:showaxis]

        xaxis!(p, false)
        @test !sp[:xaxis][:showaxis]

        yaxis!(p, true)
        @test sp[:yaxis][:showaxis]

        yaxis!(p, false)
        @test !sp[:yaxis][:showaxis]

        p = plot3d([1, 2], [1, 2], [1, 2])
        plot3d!(p, [3, 4], [3, 4], [3, 4])
        @test Plots.series_list(p[1])[1][:seriestype] == :path3d
    end

    @testset "Set Ticks" begin
        p = plot([0, 2, 3, 4, 5, 6, 7, 8, 9, 10])
        sp = p[1]

        xticks = 2:6
        xticks!(xticks)
        @test sp.attr[:xaxis][:ticks] == xticks

        xticks = 1:5
        xticks!(p, xticks)
        @test sp.attr[:xaxis][:ticks] == xticks

        yticks = 0.2:0.1:0.7
        yticks!(yticks)
        @test sp.attr[:yaxis][:ticks] == yticks

        yticks = 0.1:0.5
        yticks!(p, yticks)
        @test sp.attr[:yaxis][:ticks] == yticks

        xticks = [5, 6, 7.5]
        xlabels = ["a", "b", "c"]
        xticks!(xticks, xlabels)
        @test sp.attr[:xaxis][:ticks] == (xticks, xlabels)

        xticks = [5, 2]
        xlabels = ["b", "a"]
        xticks!(p, xticks, xlabels)
        @test sp.attr[:xaxis][:ticks] == (xticks, xlabels)

        yticks = [0.5, 0.6, 0.75]
        ylabels = ["z", "y", "x"]
        yticks!(yticks, ylabels)
        @test sp.attr[:yaxis][:ticks] == (yticks, ylabels)

        yticks = [0.5, 0.1]
        ylabels = ["z", "y"]
        yticks!(p, yticks, ylabels)
        @test sp.attr[:yaxis][:ticks] == (yticks, ylabels)
    end
end


text:
# A Lexer / Tokenizer for tokens found in plists

# Token types
const NUMBER = 1
const STRING = 2
const UNKNOWN = 3
const EOF = 4
const LPAREN = Int('(')
const RPAREN = Int(')')
const LBRACE = Int('{')
const RBRACE = Int('}')
const COMMA = Int(',')
const EQUAL = Int('=')
const SEMICOLON = Int(';')

# mapping from token type to a string representation
const token_name = Dict(NUMBER => "Number", STRING => "String", UNKNOWN => "Unknown", EOF => "EOF")

"""
A string of code is turned into an array of Tokens by the lexer. Each
symbol or word in code is represented as a Token.
"""
struct Token
    typ::Int
    lexeme::String
end

function Token(typ::Int)
    lexeme = ""
    if !haskey(token_name, typ)
       lexeme = string(Char(typ))
    end
    Token(typ, lexeme)
end

function show(io::IO, token::Token)  # avoids recursion into prev and next
    token_str = if haskey(token_name, token.typ)
        token_name[token.typ]
    else
        string(Char(token.typ))
    end

    lex_str = if token.lexeme != "" && token.lexeme != string(Char(token.typ))
        "($(token.lexeme))"
    else
        ""
    end
    print(io, "$token_str $lex_str")
end

==(t1::Token, t2::Token) = t1.typ == t2.typ && t1.lexeme == t2.lexeme

"Keeps track of string of code we want to turn into array of tokens"
mutable struct Lexer
   buffer::String
end

# define iterator interface to tokens in lexer. Requires
# implementation of start, done and next
start(lex::Lexer) = 1
done(lex::Lexer, pos::Int) = length(lex.buffer) < pos
iteratorsize(lex::Lexer) =  Base.SizeUnknown()

function next(lex::Lexer, pos::Int)
    if done(lex, pos)
        return (Token(EOF, ""), pos)
    end

    buf = lex.buffer
    ch = buf[pos]

    while ch in " \t\n"
        pos += 1
        if done(lex, pos)
            return (Token(EOF, ""), pos)
        end
        ch = buf[pos]
    end

    if in(ch, "{}(),=;")
        (Token(Int(ch), string(ch)), pos + 1)
    elseif isdigit(ch)
        range = search(buf, r"\d+", pos)
        (Token(NUMBER, buf[range]), last(range) + 1)
    elseif ch == '"'
        npos = search(buf, '"', pos + 1) # TODO: Does not handle escaped quotes
        (Token(STRING, buf[pos + 1:npos-1]), npos + 1)
    elseif isalpha(ch)
        range = search(buf, r"\w+", pos)
        (Token(STRING, buf[range]), last(range) + 1)
    else
        error("Unknown character '$ch' at position: $(pos)")
    end
end


text:
## API for menus
type Action <: Widget
    o
    toolkit
    label
    tooltip
    shortcut
    icon
    command
end

function action(parent, label::String, shortcut::Union(Nothing,String), icon::Union(Nothing, Icon), command::Function)
    o = action(parent.toolkit, parent)
    a = Action(o, parent.toolkit, label, nothing, shortcut, icon, command)
    a[:label] = label
    if !isa(shortcut, Nothing)  a[:shortcut] = shortcut end
    if !isa(icon, Nothing)  a[:icon] = icon end
    a[:command] = command
    a
end

action(parent, label::String, icon::Icon, command::Function) = action(parent, label, nothing, icon, command)
action(parent, label::String, command::Function) = action(parent, label, nothing, nothing, command)

## properties
getWidget(action::Action) = action.o
getEnabled(action::Action) = getEnabled(action.toolkit, action)
setEnabled(action::Action, value::Bool) = setEnabled(action.toolkit, action, value)

getLabel(action::Action)  = action.label
function setLabel(action::Action, label::String)
    action.label = label
    setLabel(action.toolkit, action, label)
end

getIcon(action::Action) = a.icon
function setIcon(action::Action, icon::Icon)
    action.icon = icon
    setIcon(action.toolkit, action, icon)
end

getShortcut(action::Action) = a.shortcut
function setShortcut(action::Action, shortcut::String)
    action.shortcut = shortcut
    setShortcut(action.toolkit, action, shortcut)
end

getTooltip(action::Action) = action.tooltip
function setTooltip(action::Action, tooltip::String)
    action.tooltip = tooltip
    setTooltip(action.toolkit, action, tooltip)
end

function setCommand(action::Action, command::Function)
    action.command = command
    setCommand(action.toolkit, action, command)
end

## MenuBar for attaching to a window
type MenuBar <: Container
    o
    toolkit
end

## Menu is flexible (attached to menubar, menu or widget (for popup))
type Menu <: Container
    o
    toolkit
end

## This is also menu
function menubar(parent::Window)
    mb = menubar(parent.toolkit, parent)
    MenuBar(mb, parent.toolkit)
end

menu(parent::Window) = menubar(parent)

## return Menu Instance
## not primary interface, use menu() instead.
## This is here, as Qt has it.
function addMenu(menubar::MenuBar, label::String)
    menu(menubar, label)
end

## menu is used with menubar and for popups
## parent is
## - MenuBar instance -- toplevel submenu
## - Menu instance -- will be submenu
## - widget -- menu will be popup menu, actions may be passed (i) -> if storeview, ...
function menu(parent::MenuBar, label::String)
    widget = menu(parent.toolkit, parent, label)
    Menu(widget, parent.toolkit)
end

## Add a sub menu to a parent menu
function menu(parent::Menu, label::String)
    widget = menu(parent.toolkit, parent, label)
    Menu(widget, parent.toolkit)
end

## Used for context menus
## Context, if any, is set via `widget[:context]`
##
function menu(parent::Widget)
    widget = menu(parent.toolkit, parent)
    Menu(widget, parent.toolkit)
end

## addAction for adding to a menu
function addAction(menu::Menu, action::Action)
    addAction(menu.toolkit, menu, action)
end

function addAction(menu::Menu, sep::Separator)
    addAction(menu.toolkit, menu, sep)
end

function addAction(menu::Menu, rg::RadioGroup)
    addAction(menu.toolkit, menu, rg)
end

function addAction(menu::Menu, cb::CheckBox)
    addAction(menu.toolkit, menu, cb)
end

## if a Menu is like a queue, the addAction is like push!
push!(menu::Menu, obj) = addAction(menu, obj)


text:
######################################################################
# scoring.jl
# Computes ergodic, total, and control scores
######################################################################

# These are required because multi-agent is different
function ergodic_score(em::ErgodicManager, traj::VVF, d::Dynamics)
	return ergodic_score(em, traj)
end
function ergodic_score(em::ErgodicManager, traj::VVF, gd::GroupDynamics)
	xds = vvf2vvvf(traj, gd)
	ck = decompose(em, xds)
	return ergodic_score(em, ck)
end
function ergodic_score(em::ErgodicManager, traj::VVF, vtm::VTM)
	xds = vvf2vvvf(traj, vtm)
	ck = decompose(em, xds)
	return ergodic_score(em, ck)
end

# TODO: uncomment and test for time, comparing to new method in
#  ergodic_manager.jl
# It didn't matter for other types, but maybe it will here because of
#  complex numbers and stuff?
function ergodic_score(em::ErgodicManagerSE2, ck::Array{Complex{Float64},3})
	val = 0.0

    # TODO: Must be scrubbed since em.M became range and not int
	#for m = 0:em.M
	#	for n = 0:em.N
	#		for p = 0:em.P
    for (mi,m) in enumerate(em.M)
        for (ni,n) in enumerate(em.N)
            for (pl,p) in enumerate(em.P)
				d = em.phik[mi,ni,pl] - ck[mi,ni,pl]
				dr = real(d)
				di = imag(d)
				val += em.Lambda[mi,ni,pl] * (dr*dr + di*di)
			end
		end
	end
	return val
end

"""
`control_score(ud::VVF, R, h)`

Assumes only non-zero elements of `R` are corners (not anymore I think)

`control_score(tm::TrajectoryManager,u::VVF) = control_score(u, tm.R, tm.h)`
"""
function control_score(ud::VVF, R::Matrix{Float64}, h::Float64)
	cs = 0.0
	num_u = length(ud[1])
	for ui in ud

		for j = 1:num_u
		   cs += R[j,j] * ui[j] * ui[j]
		end

		# TODO: this is how I do it. Benchmark it
		#cs += dot(ui, R*ui)

		# old way I suppose. delete this dumb stuff
		#cs += R[1,1] * ui[1] * ui[1]
		#cs += R[2,2] * ui[2] * ui[2]
	end
	return 0.5 * h * cs
end

function control_score(tm::TrajectoryManager, ud::VVF)
	return control_score(ud, tm.R, tm.h)
end

# TODO: do I even use these functions?
control_score(ud::VVF) = control_score(ud, eye(2), 1.0)
function control_score(ud::VVF, N::Int)
	cs = 0.0
	for n = 1:N
		num_u = length(ud[n])
		for j = 1:num_u
			cs += ud[n][j] * ud[n][j]
		end
		#cs += ud[n][1] * ud[n][1]
		#cs += ud[n][2] * ud[n][2]
	end
	return cs
end

function control_effort(ud::VVF, N::Int=length(ud))
	cs = 0.0
	for n = 1:N
		udx = ud[n][1]
		udy = ud[n][1]
		cs += sqrt(udx*udx + udy*udy)
	end
	return cs
end

function barrier_score(em::ErgodicManager, xd::VVF, c::Float64)
	if c == 0.0; return 0.0; end

	bs = 0.0
	xmax = x_max(em)
	ymax = y_max(em)
	xmin = x_min(em)
	ymin = y_min(em)

	for xi in xd
		if (xi[1] > xmax)
			dx = xi[1] - xmax
			bs += c * dx * dx
		elseif (xi[1] < xmin)
			dx = xi[1] - xmin
			bs += c * dx * dx
		end
		if (xi[2] > ymax)
			dy = xi[2] - ymax
			bs += c * dy * dy
		elseif (xi[2] < ymin)
			dy = xi[2] - ymin
			bs += c * dy * dy
		end
	end
	return bs
end

"""
`total_score(em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF)`

Computes the total score `q*ergodic_score + sum_n h/2 un'Rn un`
"""
function total_score(em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF)
	es = tm.q * ergodic_score(em, xd, tm.dynamics)
	cs = control_score(ud, tm.R, tm.h)
	bs = barrier_score(em, xd, tm.barrier_cost)
	return es + cs + bs
end

"""
`all_scores(em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF)`
"""
function all_scores(em::ErgodicManager, tm::TrajectoryManager, xd::VVF, ud::VVF)
	es = tm.q * ergodic_score(em, xd, tm.dynamics)
	cs = control_score(ud, tm.R, tm.h)
	ts = es + cs
	return es, cs, ts
end


text:
"""
    Graph
A generic directed graph capable of storing parallel arcs and its attributes.
"""
mutable struct Graph
    nn::Integer
    na::Integer
    fadj::Dictionary
    badj::Dictionary
end

# Accessors
nn(g::Graph) = g.nn
na(g::Graph) = g.na
fadj(g::Graph) = g.fadj
badj(g::Graph) = g.badj

# convenience
function ==(g::Graph, h::Graph)
    return g.fadj == h.fadj && g.badj == h.badj && g.nn == h.nn && g.na == h.na
end

function empty!(g::Graph)
    g.nn = 0; g.na = 0
    empty!(g.fadj); empty!(g.badj)
    return g
end

"""
    Graph(n)
Constructs a graph with `n` number of nodes.
"""
function Graph(n::Integer=0)
    fadj = dictionary(Pair{Any,Dictionary}(i, Dictionary()) for i in 1:n)
    badj = dictionary(Pair{Any,Dictionary}(i, Dictionary()) for i in 1:n)
    return Graph(n, 0, fadj, badj)
end

"""
    Graph(abunch)
Constructs a graph from a bunch of arcs `abunch` defined by tuples.
"""
function Graph(abunch::Vector)
    G = Graph(); add_arcs_from!(G, abunch)
    return Graph(G.nn, G.na, G.fadj, G.badj)
end

"""
    add_arc!(g,u,v,key)
Add an arc (u,v,key) to graph g. The default for key is nothing.
Node existence will be first checked then the keys.
"""
function add_arc!(g::Graph, u, v, key=nothing)
    if !(u in keys(g.fadj))
        insert!(g.fadj, u, Dictionary())
        insert!(g.badj, u, Dictionary())
        g.nn += 1 # add node
    end
    if !(v in keys(g.fadj))
        insert!(g.fadj, v, Dictionary())
        insert!(g.badj, v, Dictionary())
        g.nn += 1 # add node
    end
    if isnothing(key)
        key = new_arc_key(g, u, v)
    end
    if v in keys(g.fadj[u])
        innerdict = g.fadj[u][v]
        attribdict = get(innerdict, key, Dictionary())
        if !(key in keys(innerdict))
            g.na += 1 # add arc
        end
        set!(innerdict, key, attribdict)
    else
        attribdict = Dictionary()
        innerdict = Dictionary()
        set!(innerdict, key, attribdict)
        set!(g.fadj[u], v, innerdict)
        set!(g.badj[v], u, innerdict)
        g.na += 1 # add arc
    end
    return key
end

"""
    new_arc_key(g,u,v)
New key generator for arc (u,v). The key will be the lowest integer not used.
The function is internally used for adding arcs.
"""
function new_arc_key(g::Graph, u, v)
    keydict = try
        g.fadj[u][v] # no arc yet
    catch KeyError
        return 1 # lowest key is 1
    end
    key = length(keydict)
    while key in keys(keydict)
        key += 1
    end # find lowest integer not used
    return key
end

"""
    rem_arc!(g, u, v, key)
remove an arc u,v,key from the graph. u and v are nodes.
key is used to distinguish multiple arcs between pair of nodes.
"""
function rem_arc!(g::Graph, u, v, key=nothing)
    d = try
        g.fadj[u][v]
    catch KeyError
        return false # no key uv on graph
    end
    if isnothing(key)
        unset!(d, Dictionaries.lastindex(d)) # pop!
    else
        if haskey(d, key)
            delete!(d, key) # remove the item key
        else # doesn't have key
            return false
        end
    end
    if length(d) == 0
        delete!(g.fadj[u], v)
        delete!(g.badj[v], u)
    end
    g.na -= 1
    return true
end

"""
    add_arcs_from!(g, abunch)
add arcs from an iterable input. abunch is the container of the arc.
each arc can be one of 2-tuples (u,v) or 3-tuples (u,v,k).
adding an already existing arc will update its list of keys.
all non-compliant arc will be silently not added.
"""
function add_arcs_from!(g::Graph, abunch::Vector)
    keylist = []
    for a in abunch
        la = length(a)
        if la == 2 # 2-tuple
            u, v = a; key = nothing
        elseif la == 3 # 3-tuple
            u, v, key = a
        end
        res = add_arc!(g, u, v, key)
        push!(keylist, res)
    end
    return keylist
end

"""
    rem_arcs_from!(g, abunch)
remove several arcs in abunch from g. the arcs in abunch can be one of:
    - 2-tuples (u,v)
    - 3-tuples (u,v,k)
nonexisting arcs in abunch will silently not be removed.
"""
function rem_arcs_from!(g::Graph, abunch::Vector)
    removed = 0
    for a in abunch
        if length(a) >= 2
            remove = try
                rem_arc!(g, a...) # if a doesn't comply,
            catch end # just pass.
            if remove
                removed += 1
            end
        end
    end
    return removed
end

"""
    add_node!(g, u)
add a node u to graph g. silently doesn't add anything if object is already a node.
"""
function add_node!(g::Graph, u)
    if !haskey(g.fadj, u) # operation only for nonexisting
        insert!(g.fadj, u, Dictionary())
        insert!(g.badj, u, Dictionary())
        g.nn += 1
        return u # something done
    else
        return nothing
    end
end

"""
    rem_node!(g, n)
Removes the node n and all adjacent arcs.
Attempting to remove a non-existent node will raise an exception.
"""
function rem_node!(g::Graph, n)
    succs = try
        g.fadj[n]
    catch KeyError
        return false # node n not in graph
    end
    for u in keys(g.fadj[n])
        delete!(g.badj[u], n) # delete arc n-u on pred
    end
    g.na -= length(g.fadj[n]) # delete all arc from n on succ
    delete!(g.fadj, n)
    for u in keys(g.badj[n])
        g.na -= length(g.fadj[u][n])
        delete!(g.fadj[u], n) # delete arc u-n on succ
    end
    delete!(g.badj, n)
    g.nn -= 1
    return true # successful operation returns true
end

"""
    add_nodes_from!(g, nbunch)
add nodes from nbunch to g. nbunch is a container of nodes to add.
"""
function add_nodes_from!(g::Graph, nbunch::Vector)
    added = []
    for n in nbunch
        add = add_node!(g, n)
        if !isnothing(add)
            push!(added, add)
        end
    end
    return added
end

"""
    rem_nodes_from!(g, nbunch)
remove several nodes at once. nbunch are the containers of nodes to remove.
If a node in the container is not in the graph it is silently ignored.
"""
function rem_nodes_from!(g::Graph, nbunch::Vector)
    removed = 0
    for n in nbunch
        remove = try
            rem_node!(g, n) # if any key error,
        catch end # just pass.
        if remove
            removed += 1
        end
    end
    return removed
end


text:
"""
	$(TYPEDEF)
Struct to store tag id with translations from body to landmark
"""
struct TagkTl
	tagID::Symbol
	ktl::Translation
	kQl::SVector{4}#Quat #Quat does not serialize as wanted out of the box so use SVector instead
end

# ============================================================
# --------------------TagOnPose--------------------
# ============================================================
# import Arena.Amphitheatre: visualize!

struct TagOnPose <: AbstractAmphitheatre
	robotId::String
	sessionId::String

	tags::Dict{Symbol,Vector{TagkTl}}

	isdrawn::Dict{Symbol, Bool}
	#tag on pose drawing propeties
	tagProp::plDrawProp
end

TagOnPose(robotId::String, sessionId::String, tagsOnPoses::Dict{Symbol,Vector{TagkTl}};
								tagProp::plDrawProp = plDrawProp(0.2, 0.1, RGBA(1,0,1,0.5))) =
			TagOnPose(robotId, sessionId, tagsOnPoses, Dict{Symbol, Bool}(), tagProp)

"""
    $(SIGNATURES)
Basic visualizer object visualize! function.
"""
function visualize!(vis::Visualizer, top::TagOnPose)::Nothing
	#TODO maybe improve this function to lower memmory allocations

	robotId = top.robotId
	sessionId = top.sessionId

	tagsOnPoses = top.tags
	isdrawn = top.isdrawn
	tagProp = top.tagProp

	for poseKey in keys(tagsOnPoses)

		if !get(isdrawn, poseKey, false)
			haskey(isdrawn, poseKey) ? isdrawn[poseKey] = true : push!(isdrawn, poseKey=>true)

			for tag = tagsOnPoses[poseKey]

				tform = tag.ktl ∘ LinearMap(Quat(tag.kQl...))

				visPose!(vis[robotId][sessionId][:poses][poseKey][tag.tagID], tform,
							scale = tagProp.scale,
							sphereScale = tagProp.sphereScale,
							color = tagProp.color)
			end
		end
	end

    return nothing
end

# ============================================================
# --------------------GraffTagOnPose--------------------
# ============================================================
# import Arena.Amphitheatre: visualize!

struct GraffTagOnPose <: AbstractAmphitheatre
	robotId::String
	sessionId::String
	config::GraffConfig

	tags::Dict{Symbol,Vector{TagkTl}}

	isdrawn::Dict{Symbol, Bool}
	#tag on pose drawing propeties
	tagProp::plDrawProp
end

GraffTagOnPose(config::GraffConfig; tagProp::plDrawProp = plDrawProp(0.2, 0.1, RGBA(1,0,1,0.5))) =
			GraffTagOnPose(config.robotId, config.sessionId, config, Dict{Symbol,Vector{TagkTl}}(), Dict{Symbol, Bool}(), tagProp)

"""
    $(SIGNATURES)
Basic visualizer object visualize! function.
"""
function visualize!(vis::Visualizer, top::GraffTagOnPose)::Nothing

	# only read and update from graff once
	isdrawn = top.isdrawn

	if length(isdrawn) == 0

		robotId = top.robotId
		sessionId = top.sessionId

		tagsOnPoses = top.tags

		tagProp = top.tagProp
		config = top.config

		sessionDataEntries = getSessionDataEntries(config.robotId, config.sessionId)

		for nodeKey = keys(sessionDataEntries)

			dataEntries = sessionDataEntries[nodeKey]

			!("TagkTl" in [de.id for de in dataEntries]) && continue

			d = GraffSDK.getData(robotId, sessionId, dataEntries[1].neoNodeId, "TagkTl")
			jtags = JSON2.read(d.data)
			poseKey = Symbol(nodeKey)

			if !get(isdrawn, poseKey, false)
				haskey(isdrawn, poseKey) ? isdrawn[poseKey] = true : push!(isdrawn, poseKey=>true)

				for tag = jtags

					tform = Translation(tag.ktl.translation) ∘ LinearMap(Quat(tag.kQl...))

					visPose!(vis[robotId][sessionId][:poses][poseKey][Symbol(tag.tagID)], tform,
								scale = tagProp.scale,
								sphereScale = tagProp.sphereScale,
								color = tagProp.color)
				end
			end
		end
	end
    return nothing
end


text:
module Utils

macro inject(x::Symbol)
    t = eval(x)
    @assert isa(t, DataType)
    returnValue = Expr[]
    for i =1:length(t.types)
        push!(returnValue, Expr(:(::), t.names[i], t.types[i].name.name))
    end
    esc((Expr(:block, returnValue...)))
end

foreach(func::Function, collection) = for elem in collection; func(elem); end

function mapvalues(func::Union(Function, Base.Func), collection::Dict)
   [key => func(value) for (key, value) in collection]
end
function mapkeys(func::Union(Function, Base.Func), collection::Dict)
   [func(key) => value for (key, value) in collection]
end
# Simple file wrapper, which encodes the type of the file in its parameter
# Usefull for file IO
immutable File{Ending}
  abspath::UTF8String
end
File(folders...) = File(joinpath(folders...))
function File(file)
  @assert !isdir(file) "file string refers to a path, not a file. Path: $file"
  file  = abspath(file)
  path  = dirname(file)
  name  = file[length(path):end]
  ending  = rsearch(name, ".")
  ending  = isempty(ending) ? "" : name[first(ending)+1:end]
  File{symbol(ending)}(file)
end

Base.open(x::File)    = open(abspath(x))
Base.abspath(x::File) = x.abspath

function print_with_lines(text::AbstractString)
    for (i,line) in enumerate(split(shadercode, "\n"))
        @printf("%-4d: %s\n", i, line)
    end
end

end # module


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>100.0, :alg=>TMJets(abstol=1e-15)),
    "X0"   => BallInf([0.1, 0.0, 0.0, 0.5], 0.001),
    "info" => Dict(
        "name" => "henon_heiles",
        "dim" => 4,
        "linear" => false,
        "hybrid" => false
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using NumericalMethodsforEngineers

a = [1. 2. -2.; 2. 5. -4.; -2. -4. 5.]

println("\n|a[1, 1]| = $(det(a[1, 1]))\n")
println("|a[1:2, 1:2]| = $(det(a[1:2, 1:2])))\n")
println("|a| = $(det(a))")
println("\nchol(a):")
chol(a) |> display

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]

println("\n|a[1, 1]| = $(det(a[1, 1]))\n")
println("|a[1:2, 1:2]| = $(det(a[1:2, 1:2])))\n")
println("|a| = $(det(a))")

println("\nchol(a) - upper:")
upper = chol(a)
upper |> display

println("\nchol(a) - lower:")
lower = chol(a)'
lower |> display

println("\nlower * upper:")
lower * upper |> display

@test a == lower * upper


text:
load_result_comment(res, ver::Union{Void, Int, String} = nothing) =
    return Core.get_desc_comment(Core.load_res_desc(get_loadable_repo(res, ver)));

load_result_log(res, ver::Union{Void, Int, String} = nothing) =
    return Core.get_desc_log(Core.load_res_desc(get_loadable_repo(res, ver)));

load_result_data(res, ver::Union{Void, Int, String} = nothing) =
    return Core.load_res_data(get_loadable_repo(res, ver));

function get_loadable_repo(res, ver::Union{Void, Int, String} = nothing)
    if Core.is_valid_res_repo_name(res) # is_valid_res_repo_name (true)
        if ver != nothing # ver (true)
            if Core.is_valid_time_tag(ver) # is_valid_time_tag (true)
                repo = parse_repo(res, ver);
                if Core.repo_exist(repo)
                    return repo;
                else
                    error("Repo $(repo.res_dir) not found!!!")
                end
            else # is_valid_time_tag (false)
                if has_versions(res)
                    if isa(ver, Int)
                        return get_version(res, ver)
                    else
                        error("'$(ver)' must be an Int for be a ver_num")
                    end
                else
                    error("There are any version repo called '$(res)'!!!")
                end
            end # end is_valid_time_tag
        else # ver (false)
            error("You most specify a 'ver_num' or a valid 'time_tag'!!!")
        end # end ver
    else # is_valid_res_repo_name (false)
        if ver != nothing Core.validate_res_name(res) end;
        res = extract_basename(res);
        if Core.is_valid_res_dir(res) # is_valid_res_dir (true)
            repo = parse_repo(res);
            if Core.repo_exist(repo)
                return repo;
            else
                error("Repo $(repo.res_dir) not found!!!")
            end
        else # is_valid_res_dir (false)
            Core.validate_res_dir(res);
        end # end is_valid_res_dir
    end
end


text:
#   Unit examples of the FourierAnalysis Package for julia language
#   v 0.0.1 - last update 24th of September 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This example shows how to compute cross-spectra
#   and how to extract inforation from them.

#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #
#                                                                             #
#   ~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~~¤~  #

using FourierAnalysis, FFTW, LinearAlgebra, Statistics, Plots, Plots.Measures

# add module for reading the two EEG text files to be used ater
push!(LOAD_PATH, @__DIR__)
using IOtxt

# Get EEG file names with complete path (they have extension .txt)
S=getFilesInDir(@__DIR__; ext=(".txt",))

# read the two EEG data files and put them in a Matrix object
X1=readEEG(S[1])
X2=readEEG(S[2])

# Cross-Spectra of EEG data
##########################################

t, sr, slide, tapering = 512, 128, 64, harris4

# gather some attributes to obtain nice spectra plots
spectraArgs=(left_margin   = 2mm,
             bottom_margin = 2mm,
             xtickfont     = font(10, "Times"),
             ytickfont     = font(10, "Times"))

# spectra
S=spectra(X1, sr, t; tapering=tapering, func=√)
plot(S; ytitle="Amplitude (\\muV)", spectraArgs...)
# smoothed spectra
S2=spectra(X1, sr, t; tapering=tapering, smoothing=blackmanSmoother, func=√)
plot(S2; ytitle="Amplitude (\\muV)", spectraArgs...)

#cross-spectra
𝙎=crossSpectra(X1, sr, t; tapering=tapering, tril=true)
# smooth a-posteriori the cross-spectra
𝙎2=smooth(blackmanSmoother, 𝙎)

# mean cross-spectra in 8Hz-12Hz range
alpha=mean(𝙎, (8, 12))

# extract all cross-spectra in 8Hz-12Hz range
E=extract(𝙎, (8, 12))

# mean smoothed cross-spectra in 8Hz-12Hz range
alpha=mean(𝙎2, (8, 12))

# extract smoothed cross-spectra in 8Hz-12Hz range
E=extract(𝙎2, (8, 12))

# get amplitude spectra from cross-spectra and compare with S (the long way)
S__=sqrt.(Real.([𝙎.y[i][j, j] for i=1:length(𝙎.y), j=1:size(X1, 2)]))
# the short way using FourierAnalysis
S_=Spectra(𝙎; func=√)

# check they are the same
norm(S.y-S_.y)

# cross-spectra of several data matrix at once
t, sr, slide, tapering = 1024, 128, 512, harris4
𝐗=[X1, X2]

# Compute the cross-spectra altogether
𝓢=crossSpectra(𝐗, sr, t; tapering=slepians(sr, t))

# mean cross-spectrum in 8Hz-12Hz range for each cross-spectra in 𝓢
alpha=mean(𝓢, (8, 12))
# mean of the above cross-spectra in one pass
alpha=mean(mean(𝓢, (8, 12)))
# the above computation without checking for homogeneity of elements in 𝓢 (faster)
alpha=mean(mean(𝓢, (8, 12), check=false))

# extract cross-spectra in 8Hz-12Hz range for all cross-spectra in 𝓢
E=extract(𝓢, (8, 12))

# average cross-spectra in 8Hz-12Hz range for 𝓢[1] only
Y=mean(𝓢[1], (8, 12))

# This (8 is the frequency in Hz)
Y=mean(𝓢[1], 8)
# is equivalent to this (64 is the bin where discrete Fourier Frequency 8Hz is)
𝓢[1].y[64]
# and equivalent to this
Y=extract(𝓢[1], 8)
# and equivalent to this
Y=extract(𝓢[1], (8, 8))

# Compute the lower-trianguler part of the cross-spectra
𝓢=crossSpectra(𝐗, sr, t; tapering=slepians(sr, t), tril=true)

# average cross-spectra in 8Hz-12Hz range for all 𝕃Vector in 𝓢
Y=mean(𝓢, (8, 12))

# extract cross-spectra in 8Hz-12Hz range for all 𝕃Vector in 𝓢
E=extract(𝓢, (8, 12))

# cross-spectra 𝓢[1] averaged in 1Hz band-pass regions
Z=bands(𝓢[1], 1)

# all cross-spectra in 𝓢 averaged in 1Hz band-pass regions
Z=bands(𝓢, 1)

# get amplitude spectra from cross-spectra (the long way)
S_=sqrt.(Real.([𝓢[1].y[i][j, j] for i=1:length(𝓢[1].y), j=1:size(X1, 2)]))*1.5687746
plot(S_)

# the short way using FourierAnalysis
S_sl=Spectra(𝓢[1]; func=x->sqrt(x)*1.5687746)
norm(S_sl.y-S_) # must be zero

# Compute non-linear cross-spectra
𝓢=crossSpectra(𝐗, sr, t; tapering=slepians(sr, t), tril=true, nonlinear=true)


text:
using BlockArrays
using LinearAlgebra
using JuMP
using Ipopt
using Plots
using PyPlot
using Plots.PlotMeasures
using Formatting
using Printf

pyplot()
PyPlot.matplotlib.rc("text", usetex=true)
PyPlot.matplotlib.rc("font", family="serif")
Plots.scalefontsizes(1.5)

φ(x, n, σ) = exp(-x^2 / (2*σ^2)) * x^n

function dφ(x, n, σ)
    if n == 0
        exp(-x^2 / (2*σ^2)) * (-x/σ^2)
    else
        exp(-x^2 / (2*σ^2)) * (n * x^(n-1) - x^(n+1)/σ^2)
    end
end

K(x, y, σ, N) = sum(φ(x, n, σ) * φ(y, n, σ) for n=0:N)
∂K(x, y, σ, N) = sum(φ(x, n, σ) * dφ(y, n, σ) for n=0:N)
∂∂K(x, y, σ, N) = sum(dφ(x, n, σ) * dφ(y, n, σ) for n=0:N)

function D(x, σ, N)
    s = length(x)
    D0 = [K(x[i], x[j], σ, N) for i=1:s, j=1:s]
    D1 = [∂K(x[i], x[j], σ, N) for i=1:s, j=1:s]
    D2 = [∂∂K(x[i], x[j], σ, N) for i=1:s, j=1:s]
    [D0 D1; D1 D2]
end

function compute_ηV(x0, σ, N)
    s = length(x0)

    model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
    @variable(model, p[1:N+1])

    for i=1:s
        @constraint(model, sum(p[n+1] * φ(x0[i], n, σ) for n=0:N) == 1)
        @constraint(model, sum(p[n+1] * dφ(x0[i], n, σ) for n=0:N) == 0)
    end

    @objective(model, Min, sum(p[n]^2 for n=1:N+1))
    optimize!(model)

    pV = [value(p[n]) for n=1:N+1]
    ηV(x) = sum(φ(x, n, σ) * pV[n+1] for n=0:N)
end

fetitle = FormatExpr(L"$\sigma={}$")

function plot_ηV(x0, ηV, σ)
    s = length(x0)
    numpoints = 500

    x_tab = range(-1, stop=1, length=numpoints)
    y_tab = [ηV(x) for x in x_tab]

    Plots.plot(x_tab, y_tab, gridcolor="lightgrey", title=format(fetitle, σ))
    hline!([1], linestyle=:dash, linecolor=:black)

    if min(y_tab...) < -0.05
        hline!([-1], linestyle=:dash, linecolor=:black)
    end

    plot!(x0, ones(s), line=:stem, linestyle=:dash, marker=:circle, markersize=5,
          color="red")
    plot!(xlabel=L"$x$", legend=false, yformatter = yi -> @sprintf("%0.2f", yi))
    plot!(bottom_margin=5mm, right_margin=10mm, left_margin=5mm, top_margin=5mm)
end

N = 9

l = @layout[a; [b c d]; [e f g]]

p0 = Plots.plot((1:1)', label=L"$\eta_V$")
plot!(p0, (1:1)', label=L"$x_0$", color="red", linestyle=:dash,
      legend=:bottom, framestyle=:none)

x0 = [-0.8, 0.1, 0.3, 0.9]

σ = 0.285
ηV = compute_ηV(x0, σ, N)
p1 = plot_ηV(x0, ηV, σ)

σ = 0.315
ηV = compute_ηV(x0, σ, N)
p2 = plot_ηV(x0, ηV, σ)

σ = 0.375
ηV = compute_ηV(x0, σ, N)
p3 = plot_ηV(x0, ηV, σ)

x0 = [-0.1, 0.0, 0.05, 0.1]

σ = 0.05
ηV = compute_ηV(x0, σ, N)
p4 = plot_ηV(x0, ηV, σ)

σ = 0.3
ηV = compute_ηV(x0, σ, N)
p5 = plot_ηV(x0, ηV, σ)

σ = 0.5
ηV = compute_ηV(x0, σ, N)
p6 = plot_ηV(x0, ηV, σ)

Plots.plot(p0, p1, p2, p3, p4, p5, p6,
           layout=l, link=:x, size=(1000, 750),
           left_margin=5mm, right_margin=10mm)


text:
using JuliaAcademyData
using Test

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


text:
struct PriorTransformAndLogLikelihood{T,L}
    prior_transform::T
    loglikelihood::L
end

function (f::PriorTransformAndLogLikelihood)(u)
    v = f.prior_transform(u)
    return (v, f.loglikelihood(v))
end

prior_transform(f::PriorTransformAndLogLikelihood, u) = f.prior_transform(u)
function loglikelihood_from_uniform(f::PriorTransformAndLogLikelihood, u)
    return last(prior_transform_and_loglikelihood(f, u))
end
prior_transform_and_loglikelihood(f::PriorTransformAndLogLikelihood, u) = f(u)

"""
    NestedModel(loglike, prior_transform)
    NestedModel(loglike, priors::AbstractVector{<:Distribution})

`loglike` must be callable with a signature `loglike(::AbstractVector)` where the length of the vector must match the number of parameters in your model.

`prior_transform` must be a callable with a signature `prior_transform(::AbstractVector)` that returns the transformation from the unit-cube to parameter space. This is effectively the quantile or ppf of a statistical distribution. For convenience, if a vector of `Distribution` is provided (as a set of priors), a transformation function will automatically be constructed using `Distributions.quantile`.

**Note:**
`loglike` is the only function used for likelihood calculations. This means if you want your priors to be used for the likelihood calculations they must be manually included in the `loglike` function.
"""
struct NestedModel{F} <: AbstractModel
    prior_transform_and_loglikelihood::F
end

function NestedModel(loglike, prior_transform)
    return NestedModel(PriorTransformAndLogLikelihood(prior_transform, loglike))
end

function NestedModel(loglike, priors::AbstractVector{<:UnivariateDistribution})
    prior_transform(X) = quantile.(priors, X)
    return NestedModel(loglike, prior_transform)
end

function prior_transform(model, args...)
    return first(prior_transform_and_loglikelihood(model, args...))
end

function prior_transform(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...)
    return prior_transform(model.prior_transform_and_loglikelihood, args...)
end

function loglikelihood_from_uniform(model, args...)
    return last(prior_transform_and_loglikelihood(model, args...))
end

function loglikelihood_from_uniform(model::NestedModel{<:PriorTransformAndLogLikelihood}, args...)
    return loglikelihood_from_uniform(model.prior_transform_and_loglikelihood, args...)
end

function prior_transform_and_loglikelihood(model::NestedModel, args...)
    return model.prior_transform_and_loglikelihood(args...)
end


text:
struct FileDateTimeVariate <: AbstractDateTimeVariate
    name::String
    length::Int64
    datapath::String
end

Base.length(dtvar::FileDateTimeVariate) = dtvar.length

function slice(dtvar::FileDateTimeVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    sizeT = sizeof(Int64)
    iostream = open(dtvar.datapath)
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{Int64}(undef, slicelength)
    map(Seq(Vector{Int64}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{Int64,1,Array{Int64,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
end

function Base.map(dtvar::FileDateTimeVariate, dataframe::AbstractDataFrame)
    dataframe[getname(dtvar)]
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
function get_option_chain(ticker::AbstractString)
    kwargs = (; symbol = ticker, toDate="2022-01-01")

    uri = construct_api("marketdata/chains", kwargs)
    @show uri
    head = ["Authorization" => "Bearer "*AUTH_KEY.ACCESS_TOKEN]
    # Authorization is a header
    return @pipe HTTP.get(uri, head).body |>
                 JSON3.read  # |> _[Symbol(ticker)]

end

# GET OPTION CHAIN FOR TESLA TO END OF THE YEAR
#GET /v1/marketdata/chains?apikey=LNTJ7SD2D3UP4POVFACAOT608AROD38K&symbol=TSLA&toDate=2022-01-01 HTTP/1.1
# Accept:
# */*
# Accept-Encoding:
# gzip
# Accept-Language:
# en-US
# Authorization:
# Bearer ByvEO4hay1hC..............
# Host:
# api.tdameritrade.com
# NS-Proxy-Client-IP:
# 67.168.170.195
# sec-ch-ua:
# Chromium";v="88
# sec-ch-ua-mobile:
# ?0
# Sec-Fetch-Dest:
# empty
# Sec-Fetch-Mode:
# cors
# Sec-Fetch-Site:
# same-site
# User-Agent:
# Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML
# X-Forwarded-For:
# 10.152.4.23
# X-Forwarded-Port:
# 59002
# X-Forwarded-Proto:
# http

#RESPONSE
# HTTP/1.1 200 OK
# Access-Control-Allow-Headers:
# origin
# Access-Control-Allow-Methods:
# GET
# Access-Control-Allow-Origin:
# https://developer.tdameritrade.com
# Access-Control-Max-Age:
# 3628800
# Cache-Control:
# no-cache
# Connection:
# keep-alive
# Content-Encoding:
# gzip
# Content-Security-Policy:
# frame-ancestors 'self'
# Content-Type:
# application/json;charset=UTF-8
# Date:
# Thu, 25 Feb 2021 03:32:22 GMT
# Strict-Transport-Security:
# max-age=31536000; includeSubDomains
# Transfer-Encoding:
# chunked
# Vary:
# Accept-Encoding
# X-Application-Context:
# application:8080
# X-Content-Type-Options:
# nosniff
# X-Frame-Options:
# SAMEORIGIN
# X-Xss-Protection:
# 1; mode=block
# {
#   "symbol": "TSLA",
#   "status": "SUCCESS",
#   "underlying": null,
#   "strategy": "SINGLE",
#   "interval": 0,
#   "isDelayed": false,
#   "isIndex": false,
#   "interestRate": 0.1,
#   "underlyingPrice": 738.475,
#   "volatility": 29,
#   "daysToExpiration": 0,
#   "numberOfContracts": 4660,
#   "callExpDateMap": {
#     "2021-02-26:2": {
#       "200.0": [
#         {
#           "putCall": "CALL",
#           "symbol": "TSLA_022621C200",
#           "description": "TSLA Feb 26 2021 200 Call (Weekly)",
#           "exchangeName": "OPR",
#           "bid": 541.15,


text:
# src: http://ajkl.github.io/2015/08/10/Titanic-with-Julia/

# DataFrame introduction
# https://en.wikibooks.org/wiki/Introducing_Julia/DataFrames

# Julia documentation index.html
# https://docs.julialang.org/en/v1/manual/documentation/index.html

# Plots.jl
# https://github.com/JuliaPlots/Plots.jl
    # backends: http://docs.juliaplots.org/latest/backends/
        # PyPlot https://github.com/JuliaPy/PyPlot.jl
            # examples: https://docs.juliaplots.org/latest/examples/pyplot/

# Load libraries
using DataArrays
using DataFrames
using CSV
using FreqTables
using StatsBase

using Gadfly

# Snippet to use show() and println() in one line
showln(x) = (show(x); println())
showln("Libraries loaded")

# Load dataset

train_set = CSV.read("titanic/train.csv")
test_set = CSV.read("titanic/test.csv")

showln("Train and Test set loaded")

# Describe loaded train set
println("train_set head, 3 elements:")
showln(first(train_set, 3))
println("train_set tail, 3 elements:")
showln(last(train_set, 3))
println("train_set column names:")
showln(names(train_set))
println("train_set column types:")
showln(eltypes(train_set))
println("train_set shape:")
showln(size(train_set))
println("Describe train_set:")
show(describe(train_set), allcols=true)

# Describe loaded test set
println("test_set shape:")
showln(size(test_set))
println("Describe test_set:")
show(describe(test_set), allcols=true)

# Exploratory data analysis

# Explore unique values in every column

# Way 0: custom
function show_unique_values(dataset::DataFrame, max_unique_vals::Int64=10)
    for column_symbol in names(dataset)
        unique_col_values = unique(dataset[column_symbol])
        if size(unique_col_values, 1) < max_unique_vals
            println(column_symbol, ", unique values are ", unique_col_values)
        else
            println(column_symbol, " has too many unique values to display")
        end
    end
end

println("Showing unique values for train_set:")
show_unique_values(train_set)

println()

println("Showing unique values for test_set:")
show_unique_values(test_set)

# Way 1
survive_zero = count(i -> (i == 0), train_set[:Survived])
survive_one = count(i -> (i == 1), train_set[:Survived])
display(survive_zero)
display(survive_one)

# Way 2: "using FreqTables"
survive_freq_table = freqtable(train_set[:Survived])
display(survive_freq_table)

# Way 3: "using StatsBase"
display(countmap(train_set[:Survived]))

# Way 4: "using StatsBase"
display(counts(train_set[:Survived]))

# View proportions of people survived / not survived
display(proportions(train_set[:Survived]))
display(proportionmap(train_set[:Survived]))

# "counts" does not work for categorical variables.
# Use "countmap()" instead
display(countmap(train_set[:Sex]))

# Create dimension: indicate if person was a child or not
train_set[:Child] = 1
train_set[isna(train_set[:Age]), :Child] = 1

train_set[train_set[:Age] .< 18, :Child] = 1
train_set[train_set[:Age] .> 18, :Child] = 0

display(head(train_set))

# NOTE: this notebook isn't finished

max_stars_repo_name: honchardev/Fun

max_issues_repo_name: honchardev/Fun

max_forks_repo_name: honchardev/Fun


text:
function delta(x::Int64)
   if x == 0
      return 1
   else
      return 0
   end
end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
using IVPTestSuite  # imports ODE, DASSL and Sundials
const S = Solvers # for convenience

# choose test-case (all of them are in IVPTestSuite.tc_all)
tc_name = [:plei
           :vdpol
           :rober
           :threebody
           :bruss1d
           :chemakzo
           :hires][4]
tc = IVPTestSuite.tc_all[tc_name]

# Pick a solver.  S.allsolvers is the list of all solvers.  Generally
# each family of solvers also has a list.
solver = S.allsolvers[[ODE.ode78,
                       Sundials.CVode,
                       Sundials.IDASolve,
                       DASSL.dasslSolve][3]
                      ]

# make a TestRun which combines a TestCase with a Solver + some extras:
##Example protocol For adaptive solvers
abstol = 1e-12
reltol = abstol
dt0 = NaN # size of first step
tr = TestRunAdapt(tc, solver, Dict{Symbol,Any}(), abstol, reltol, dt0)

# For fixed step solvers use TestRunFixedStep instead. See ex3.jl for such a case.

# run it and display results:
re = IVPTestSuite.run_ode_test_throwerror(tr) # This returns a TestResults instance
re = IVPTestSuite.run_ode_test_throwerror(tr) # Run it twice to get accurate timings
show(re)


text:
include("../OnlineAdLearning/src/OnlineAdLearning.jl")

using JuMP, Gurobi

const TIME_LIMIT = 3 * 60.0

const solver = JuMP.with_optimizer(Gurobi.Optimizer, TimeLimit=TIME_LIMIT)
const root_solver = JuMP.with_optimizer(Gurobi.Optimizer, TimeLimit=TIME_LIMIT, NodeLimit=1, Crossover=0)
const relaxation_solver = JuMP.with_optimizer(Gurobi.Optimizer, TimeLimit=TIME_LIMIT, Method=2, Crossover=0)
const dc_solver = JuMP.with_optimizer(Gurobi.Optimizer, TimeLimit=TIME_LIMIT, OutputFlag=0)

const ALGORITHMS = Dict(
    "cp" => OnlineAdLearning.ConstantPolicy(),
    "clairvoyant" => OnlineAdLearning.Clairvoyant(),
    "dc" => OnlineAdLearning.DifferenceOfConvex(dc_solver),
    "mip" => OnlineAdLearning.MIP(solver),
    "mip_root" => OnlineAdLearning.MIP(root_solver),
    "lp" => OnlineAdLearning.LPRelaxation(relaxation_solver),
    "ga" => OnlineAdLearning.GradientAscent()
)

const result_path = "result.csv"

const num_features = 50

const num_train_points = 1000
const num_validation_points = 5000
const num_test_points = 5000
const num_samples = num_train_points + num_validation_points + num_test_points

const num_instances = 3

const HYPERPARAMS = Dict(
    "baseline" => (ρ = 0.9, σ = 0.1, α = 0.1),
    "high-noise" => (ρ = 0.9, σ = 0.5, α = 0.1),
    "low-correlation" => (ρ = 0.5, σ = 0.1, α = 0.1),
    "low-margin" => (ρ = 0.9, σ = 0.1, α = 0.02)
)

const ORDERED_KEYS = ["baseline", "high-noise", "low-correlation", "low-margin"]

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
# add gloabl shape transformation parameters and orthonormalize all vectors
function global_shape_transform(s0, pc)
    npc = size(pc, 2)
    np = @compat Int(length(s0) / 2)
    # columns 1:4 - global transform params
    # columns 5:end - shape principal components
    s_star_pc = zeros(2*np, npc+4)
    s_star_pc[:, 1] = s0
    s_star_pc[1:np, 2] = -s0[np+1:end]
    s_star_pc[np+1:end, 2] = s0[1:np]
    s_star_pc[1:np, 3] = ones(np)
    s_star_pc[np+1:end, 3] = zeros(np)
    s_star_pc[1:np, 4] = zeros(np)
    s_star_pc[np+1:2*np, 4] = ones(np)
    s_star_pc[:, 5:end] = pc
    # orthonormalizing all
    s_star_pc = gs_orthonorm(s_star_pc)
    # splitting back into global transformation params and star
    s_star, S = s_star_pc[:, 1:4], s_star_pc[:, 5:end]
    return s_star, S
end

function build_shape_model(m::AAModel, shapes::Vector{Shape})
    mean_shape, shapes_aligned = align_shapes(shapes)
    mini, minj = minimum(mean_shape[:, 1]), minimum(mean_shape[:, 2])
    maxi, maxj = maximum(mean_shape[:, 1]), maximum(mean_shape[:, 2])
    # move mean shape to upper-left corner
    mean_shape[:, 1] = mean_shape[:, 1] .- (mini - 2)
    mean_shape[:, 2] = mean_shape[:, 2] .- (minj - 2)
    frame = ModelFrame(round(Int, mini), round(Int, minj),
                       round(Int, maxi), round(Int, maxj))
    shape_mat = datamatrix(Shape[shape .- mean_shape for shape in shapes_aligned])
    shape_pca = MultivariateStats.fit(PCA, shape_mat)
    pc = projection(shape_pca)
    # base shape, global transform shape and transformation matrix
    s0 = flatten(mean_shape)
    s_star, S = global_shape_transform(s0, pc)
    return frame, s0, s_star, S
end

function import_shape_model(model_dir)
    s0 = squeeze(matread(joinpath(model_dir, "s0.mat"))["s0"], 2)
    s_star = matread(joinpath(model_dir, "s_star.mat"))["s_star"]
    S = matread(joinpath(model_dir, "S.mat"))["S"]
    mean_shape = reshape(s0, @compat(Int(length(s0) / 2)), 2)
    mini, minj = minimum(mean_shape[:, 1]), minimum(mean_shape[:, 2])
    maxi, maxj = maximum(mean_shape[:, 1]), maximum(mean_shape[:, 2])
    mean_shape[:, 1] = mean_shape[:, 1] .- (mini - 2)
    mean_shape[:, 2] = mean_shape[:, 2] .- (minj - 2)
    frame = ModelFrame(round(Int, mini), round(Int, minj),
                       round(Int, maxi), round(Int, maxj))
    return frame, s0, s_star, S
end

function import_triangulation(model_dir)
    return matread(joinpath(model_dir, "trigs.mat"))["trigs"]
end

function build_app_model{N}(m::AAModel, imgs::Vector{Array{Float64, N}},
                         shapes::Vector{Shape})
    app_mat = zeros(m.frame.h * m.frame.w * m.nc, length(imgs))
    for i=1:length(imgs)
        warped = warp_to_mean_shape(m, imgs[i], shapes[i])
        app_mat[:, i] = flatten(warped)
    end
    A0 = squeeze(mean(app_mat, 2), 2)
    A = projection(MultivariateStats.fit(PCA, app_mat .- A0))
    mean_app = reshape(A0, m.frame.h, m.frame.w, m.nc)
    dA0 = Array(Grad2D, m.nc)
    for i=1:m.nc
        dA0[i] = gradient2d(mean_app[:, :, i], m.wparams.warp_map)
    end
    return A0, A, dA0
end

function jacobians(m::AAModel)
    # jacobians have form (i, j, axis, param_index)
    dW_dp = zeros(m.frame.h, m.frame.w, 2, size(m.S, 2))
    dN_dq = zeros(m.frame.h, m.frame.w, 2, 4)
    for j=1:m.frame.w
        for i=1:m.frame.h
            if m.wparams.warp_map[i, j] != 0
                t = m.wparams.trigs[m.wparams.warp_map[i, j], :]
                # for each vertex
                for k=1:3
                    dik_dp = m.S[t[k], :]
                    djk_dp = m.S[t[k]+m.np, :]

                    dik_dq = m.s_star[t[k], :]
                    djk_dq = m.s_star[t[k]+m.np, :]

                    t2 = copy(t)
                    t2[1] = t[k]
                    t2[k] = t[1]

                    # vertices of the triangle in the mean shape
                    i1 = m.s0[t2[1]]
                    j1 = m.s0[m.np + t2[1]]
                    i2 = m.s0[t2[2]]
                    j2 = m.s0[m.np + t2[2]]
                    i3 = m.s0[t2[3]]
                    j3 = m.s0[m.np + t2[3]]

                    # compute the two barycentric coordinates
                    den = (i2 - i1) * (j3 - j1) - (j2 - j1) * (i3 - i1)
                    alpha = ((i - i1) * (j3 - j1) - (j - j1) * (i3 - i1)) / den
                    beta = ((j - j1) * (i2 - i1) - (i - i1) * (j2 - j1)) / den

                    dW_dij = 1 - alpha - beta

                    dW_dp[i,j,:,:] = (squeeze(dW_dp[i,j,:,:], (1, 2)) +
                                      dW_dij * [dik_dp; djk_dp])
                    dN_dq[i,j,:,:] = (squeeze(dN_dq[i,j,:,:], (1, 2)) +
                                      dW_dij * [dik_dq; djk_dq])
                end
            end
        end
    end
    return dW_dp, dN_dq
end

function sd_images(m::AAModel)
    app_modes = reshape(m.A, m.frame.h, m.frame.w, m.nc, size(m.A, 2))
    SD = zeros(m.frame.h, m.frame.w, m.nc, 4 + size(m.dW_dp, 4))
    # SD images for 4 global transformation parameters
    for i=1:4
        prj_diff = zeros(m.nc, size(m.A, 2))
        for j=1:size(m.A, 2)
            for c=1:m.nc
                prj_diff[c, j] = sum(app_modes[:,:,c,j] .*
                                     (m.dA0[c].di .* m.dN_dq[:,:,1,i] +
                                      m.dA0[c].dj .* m.dN_dq[:,:,2,i]))
            end
        end
        for c=1:m.nc
            SD[:,:,c,i] = (m.dA0[c].di .* m.dN_dq[:,:,1,i] +
                           m.dA0[c].dj .* m.dN_dq[:,:,2,i])
        end
        for j=1:size(m.A, 2)
            for c=1:m.nc
                SD[:,:,c,i] = SD[:,:,c,i] - prj_diff[c,j] * app_modes[:,:,c,j]
            end
        end
    end
    # SD images for shape parameters
    for i=1:size(m.dW_dp, 4)
        prj_diff = zeros(m.nc, size(m.A, 2))
        for j=1:size(m.A, 2)
            for c=1:m.nc
                prj_diff[c,j] = sum(app_modes[:,:,c,j] .*
                                    (m.dA0[c].di .* m.dW_dp[:,:,1,i] +
                                     m.dA0[c].dj .* m.dW_dp[:,:,2,i]))
            end
        end
        for c=1:m.nc
            SD[:,:,c,i+4] = (m.dA0[c].di .* m.dW_dp[:,:,1,i] +
                             m.dA0[c].dj .* m.dW_dp[:,:,2,i])
        end
        for j=1:size(m.A, 2)
            for c=1:m.nc
                SD[:,:,c,i+4] = SD[:,:,c,i+4] - prj_diff[c,j] * app_modes[:,:,c,j]
            end
        end
    end
    SDf = zeros(size(SD, 4), size(m.A, 1))
    for i=1:size(SD, 4)
        SDf[i,:] = flatten(SD[:,:,:,i])
    end
    return SDf
end

function train{N}(m::AAModel,
                  imgs::Vector{Array{Float64, N}},
                  shapes::Vector{Shape})
    @assert(length(imgs) >=  1, "At least one image is required")
    @assert(length(imgs) == length(shapes),
            "Different number of images and landmark sets")
    @assert(0 <= minimum(imgs[1]) && maximum(imgs[1]) <= 1,
            "Images should be in Float64 format with values in [0..1]")
    m.nc = size(imgs[1], 3)
    m.np = size(shapes[1], 1)
    m.frame, m.s0, m.s_star, m.S = build_shape_model(m, shapes)
    mean_shape = reshape(m.s0, m.np, 2)
    m.wparams = pa_warp_params(mean_shape, delaunayindexes(mean_shape),
                            size(m.frame))
    m.A0, m.A, m.dA0 = build_app_model(m, imgs, shapes)
    m.dW_dp, m.dN_dq = jacobians(m)
    m.SD = sd_images(m)
    m.H = m.SD * m.SD'
    m.invH = inv(m.H)
    m.R = m.invH * m.SD
    return m
end