text:
Pkg.add("Cxx")
Pkg.add("Nemo")
# Pkg.add("LLLplus")
Pkg.add("Plots")

########################################
# Int64 bases: plot of time vs length of smallest basis as delta is varied
########################################

using Cxx
using Libdl
import Nemo
import LLLplus
import Statistics
using Plots
gr();

# following assumes fplll is in a sibling dir. It could also be in a system directory.
const path_to_lib = "../fplll/"
addHeaderDir(path_to_lib, kind=C_System)
addHeaderDir("/usr/local/include", kind=C_System)
Libdl.dlopen(path_to_lib * "/lib/libfplll.dylib", Libdl.RTLD_GLOBAL)

cxxinclude("fplll.h")
cxxinclude("fplll/main.h")
cxx"""#include <iostream>"""

cxx"""ZZ_mat<long> b;"""
cxx"""ZZ_mat<long> u;"""
cxx"""ZZ_mat<long> u_inv;"""
cxx"""Options o;"""
cxx"""int status, flags;"""

icxx"""flags = LLL_DEFAULT;""";
icxx"""o.action = ACTION_LLL;""";
icxx"""o.eta = LLL_DEF_ETA;""";     # 0.51; see defs.h
icxx"""o.method = LM_FAST;""";    # see defs.h
icxx"""o.int_type = ZT_LONG;""";      # use double  (Int64?). See defs.h line 197
icxx"""o.float_type = FT_DOUBLE;""";  # use double  (Float64?). See defs.h line 206
# ?does generating or passing out unimodular slow the lll down?
function fplll()
    icxx"""status=lll_reduction(b,u,u_inv,o.delta,o.eta,o.method,o.float_type,o.precision,flags);"""
end

# I typically copy-and-pasted the code from here to the plot multiple times
# when changing things like the bit depth or largest dimension.

δv=[collect(.49:.05:.99); .995]
Nδ = length(δv);
bits = 25;
Nsamp = 8
N = 100
k = Int64(floor(N/2))
times = zeros(2,Nδ,Nsamp);
small = zeros(2,Nδ,Nsamp);
icxx"""b.resize($N,$N);""";
icxx"""u.resize($N,$N);""";
b = zeros(Int64,N,N);
bfp = zeros(Int64,N,N);

for nx=1:Nδ
    δ = δv[nx]
    ctx = Nemo.lll_ctx(δ, 0.51)

    icxx"""o.delta = $(δ);""";

    for sx = 1:Nsamp

        icxx"""b.gen_qary($k,$bits);""";
        for ix = 1:N, jx = 1:N
            b[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        end

        times[1,nx,sx] = @elapsed B,T = LLLplus.lll(b,δ)
        small[1,nx,sx] = B[:,1]'*B[:,1]
        times[2,nx,sx] = @elapsed fplll()
        for ix = 1:N, jx = 1:N
            bfp[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        end
        small[2,nx,sx] = bfp[:,1]'*bfp[:,1]

    end

tave  = Statistics.mean(times,dims=3)[:,:,1];
save  = Statistics.mean(small,dims=3)[:,:,1];
display([δv'; tave])
display([δv'; save])
end

plot(save',tave',label=["LLLplus.lll Int64",
                       "fplll.lll_w_transform Int64"],
                       linewidth=3,legend=:topleft)
xaxis!("norm of smallest vector", :log10)
yaxis!("time", :log10)
title!("fplll's gen_qary(floor(N/2),$bits) used")

savefig("timeVsmallest_$(bits)bitsInt64.png")

########################################
# Int64 bases: plot of time vs dimenstion
########################################

using Cxx
using Libdl
import Nemo
import LLLplus
import Statistics
using Plots
gr();

# following assumes fplll is in a sibling dir. It could also be in a system directory.
const path_to_lib = "../fplll/"
addHeaderDir(path_to_lib, kind=C_System)
addHeaderDir("/usr/local/include", kind=C_System)
Libdl.dlopen(path_to_lib * "/lib/libfplll.dylib", Libdl.RTLD_GLOBAL)

cxxinclude("fplll.h")
cxxinclude("fplll/main.h")
cxx"""#include <iostream>"""

cxx"""ZZ_mat<long> b;"""
cxx"""ZZ_mat<long> u;"""
cxx"""ZZ_mat<long> u_inv;"""
cxx"""Options o;"""
cxx"""int status, flags;"""

icxx"""flags = LLL_DEFAULT;""";
icxx"""o.action = ACTION_LLL;""";
icxx"""o.delta = LLL_DEF_DELTA;"""; # 0.99
icxx"""o.eta = LLL_DEF_ETA;""";     # 0.51; see defs.h
icxx"""o.method = LM_FAST;""";    # see defs.h
icxx"""o.int_type = ZT_LONG;""";      # use double  (Int64?). See defs.h line 197
icxx"""o.float_type = FT_DOUBLE;""";  # use double  (Float64?). See defs.h line 206
# ?does generating or passing out unimodular slow the lll down?
function fplll()
    icxx"""status=lll_reduction(b,u,u_inv,o.delta,o.eta,o.method,o.float_type,o.precision,flags);"""
end

# I typically copy-and-pasted the code from here to the plot multiple times
# when changing things like the bit depth or largest dimension.

δ=.99;
bits = 35;
ctx = Nemo.lll_ctx(δ, 0.51)
Nsamp = 8
Nvec = Int64.(round.(2 .^(2:7)))
Nnum = length(Nvec);
times = zeros(2,Nnum,Nsamp);
small = zeros(2,Nnum,Nsamp);
for nx=1:Nnum
    N = Nvec[nx]
    k = Int64(floor(N/2))

    icxx"""b.resize($N,$N);""";
    icxx"""u.resize($N,$N);""";

    b = zeros(Int64,N,N);
    bfp = zeros(Int64,N,N);

    for sx = 1:Nsamp
        icxx"""b.gen_qary($k,$bits);""";
        for ix = 1:N, jx = 1:N
            b[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        end

        times[1,nx,sx] = @elapsed B,T = LLLplus.lll(b,δ)
        small[1,nx,sx] = B[:,1]'*B[:,1]
        times[2,nx,sx] = @elapsed fplll()
        for ix = 1:N, jx = 1:N
            bfp[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        end
        small[2,nx,sx] = bfp[:,1]'*bfp[:,1]

    end
end

tave  = Statistics.mean(times,dims=3)[:,:,1];
save  = Statistics.mean(small,dims=3)[:,:,1];
[Nvec'; tave]
[Nvec'; save]

plot(Nvec,tave',label=["LLLplus.lll Int64",
                       "fplll.lll_w_transform Int64"],
                       linewidth=3,legend=:topleft)
xaxis!("log10(N)", :log10)
yaxis!("time", :log10)
title!("fplll's gen_qary(floor(N/2),$bits) used")

savefig("timeVdim_$(bits)bitsInt64.png")

########################################
# GMP bases: plot of time vs dimenstion
########################################

using Cxx
using Libdl
import Nemo
import LLLplus
import Statistics
using Plots
gr();

# following assumes fplll is in a sibling dir. It could also be in a system directory.
const path_to_lib = "../fplll/"
addHeaderDir(path_to_lib, kind=C_System)
addHeaderDir("/usr/local/include", kind=C_System)
Libdl.dlopen(path_to_lib * "/lib/libfplll.dylib", Libdl.RTLD_GLOBAL)

cxxinclude("fplll.h")
cxxinclude("fplll/main.h")
cxx"""#include <iostream>"""

cxx"""ZZ_mat<mpz_t> b;"""
cxx"""ZZ_mat<mpz_t> u;"""
cxx"""ZZ_mat<mpz_t> u_inv;"""
# cxx"""ZZ_mat<long> b;"""
# cxx"""ZZ_mat<long> u;"""
# cxx"""ZZ_mat<long> u_inv;"""
cxx"""Options o;"""
cxx"""int status, flags;"""

icxx"""flags = LLL_DEFAULT;""";
icxx"""o.action = ACTION_LLL;""";
icxx"""o.delta = LLL_DEF_DELTA;"""; # 0.99
icxx"""o.eta = LLL_DEF_ETA;""";     # 0.51; see defs.h
#icxx"""o.method = LM_WRAPPER;""";    # see defs.h
#icxx"""o.method = LM_HEURISTIC;""";    # see defs.h
icxx"""o.method = LM_PROVED;""";    # see defs.h
icxx"""o.int_type = ZT_MPZ;""";      # use GMP
icxx"""o.float_type = FT_MPFR;""";
# icxx"""o.method = LM_FAST;""";    # see defs.h
# icxx"""o.int_type = ZT_LONG;""";      # use double  (Int64?). See defs.h line 197
# icxx"""o.float_type = FT_DOUBLE;""";  # use double  (Float64?). See defs.h line 206
# does generating or passing out unimodular slow the lll down?
function fplll()
    icxx"""status=lll_reduction(b,u,u_inv,o.delta,o.eta,o.method,o.float_type,o.precision,flags);"""
end

δ=.99;
bits = 35;
ctx = Nemo.lll_ctx(δ, 0.51)
Nsamp = 8
Nvec = Int64.(round.(2 .^(2:6)))
Nnum = length(Nvec);
times = zeros(3,Nnum,Nsamp);
small = zeros(3,Nnum,Nsamp);
icxx"""o.precision = 100;""";
setprecision(BigFloat,100)
for nx=1:Nnum
    N = Nvec[nx]
    k = Int64(floor(N/2))

    icxx"""b.resize($N,$N);""";
    icxx"""u.resize($N,$N);""";

    S = Nemo.MatrixSpace(Nemo.ZZ, N,N);

    # b = zeros(Int64,N,N);
    # bfp = zeros(Int64,N,N);
    b = zeros(BigInt,N,N);
    bfp = zeros(BigInt,N,N);

    for sx = 1:Nsamp
        icxx"""b.gen_qary($k,$bits);""";
        for ix = 1:N, jx = 1:N
            b[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        end
        bn = S(Nemo.fmpz.(b')); # fmpz is GMP

        times[1,nx,sx] = @elapsed B,T = LLLplus.lll(b,δ)
        small[1,nx,sx] = B[:,1]'*B[:,1]
        times[2,nx,sx] = @elapsed fplll()
        for ix = 1:N, jx = 1:N
            bfp[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        end
        small[2,nx,sx] = bfp[:,1]'*bfp[:,1]
        times[3,nx,sx] = @elapsed Bn, Tn = Nemo.lll_with_transform(bn,ctx);
        small[3,nx,sx] = Float64((Bn[1,:]*Bn[1,:]')[1,1])

    end
end

tave  = Statistics.mean(times,dims=3)[:,:,1];
save  = Statistics.mean(small,dims=3)[:,:,1];
[Nvec'; tave]
[Nvec'; save]

# plot(Nvec,tave',label=["LLLplus.lll Int64",
#                        "fplll.lll_w_transform Int64","Nemo.lll GMP"],
#                        linewidth=3,legend=:topleft)
plot(Nvec,tave',label=["LLLplus.lll GMP",
                       "fplll.lll_w_transform GMP","Nemo.lll GMP"],
                       linewidth=3,legend=:topleft)
xaxis!("log10(N)", :log10)
yaxis!("time", :log10)
title!("fplll's gen_qary(floor(N/2),$bits) used")

savefig("timeVdim_$(bits)bitsGMP.png")
#savefig("timeVdim_$(bits)bits.png")

#logt = log10.(tave*1e6)./ repeat(log10.(Nvec'),2,1);
#plot(Nvec,logt')
#xaxis!("log10(N)", :log10)

    bfp = zeros(Int64,N,N);
    for ix = 1:N, jx = 1:N
        bfp[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
    end

# ratio = tave[2,:] ./ tave[3,:]
# plot(Nvec,ratio)
# xaxis!("log10(N)", :log10)
# yaxis!("ratio", :log10)

δ=.99;
ctx = Nemo.lll_ctx(δ, 0.51)
Nsamp = 8
N = 20
Bvec = 4:12:120;
#Bvec = 5:50:505;
Bnum = length(Bvec);
icxx"""b.resize($N,$N);""";
icxx"""u.resize($N,$N);""";
S = Nemo.MatrixSpace(Nemo.ZZ, N,N);
times = ones(3,Bnum,Nsamp);

#b = zeros(BigInt,N,N);
#b = zeros(Int64,N,N);

for bx=1:Bnum
    Bt = Bvec[bx]

    for sx = 1:Nsamp
        icxx"""b.gen_uniform($Bt);""";
        # for ix = 1:N, jx = 1:N
        #     b[jx,ix] = icxx"""b[$(ix-1)][$(jx-1)].get_si();"""
        # end
        # bn = S(Nemo.fmpz.(b')); # fmpz is also GMP
        setprecision(BigFloat,Bt+1)

        b = rand(0:big(2)^Bt,N,N)
        bn = S(Nemo.fmpz.(b')); # fmpz is also GMP

        times[1,bx,sx] = @elapsed B,T = LLLplus.lll(b,δ)
        times[2,bx,sx] = @elapsed fplll()
        times[3,bx,sx] = @elapsed Bn, Tn = Nemo.lll_with_transform(bn,ctx);
    end
end
tave  = Statistics.mean(times,dims=3)[:,:,1]

plot(Bvec,tave',label=["LLLplus GMP","fplll GMP","Nemo GMP"],
                       linewidth=3,legend=:topleft)
xaxis!("B (Number of bits)")
yaxis!("time", :log10)
title!("$(N)x$(N) matrices were used")

savefig("timeVnbits_Neq25.png")


text:
using GoogleSheets

# Example based upon: # https://developers.google.com/sheets/api/quickstart/python

client = sheets_client(AUTH_SCOPE_READONLY)

# The ID and range of a sample spreadsheet.
SAMPLE_SPREADSHEET_ID = "1BxiMVs0XRA5nFMdKvBdBZjgmUUqptlbs74OgvE2upms"
SAMPLE_RANGE_NAME = "Class Data!A2:E"

sheet = Spreadsheet(SAMPLE_SPREADSHEET_ID)
range = CellRange(sheet, SAMPLE_RANGE_NAME)
result = get(client, range)
println("RESULT: $(result)")

if isnothing(result.values)
    println("No data found.")
else
    for row in eachrow(result.values)
        println("ROW: $row")
    end

    println("")
    println("Name, Major:")
    for row in eachrow(result.values)
        # Print columns A and E, which correspond to indices 1 and 5.
        println("ROW: $(row[1]), $(row[5])")
    end
end

ext: jl


text:
@testset "Identity inflation" begin
    A = IdentityInflation()

    Nx = 10
    Ne = 20
    X = zeros(Nx, Ne)
    X .= randn(Nx, Ne)

    Xinflation = deepcopy(X)
    A(X)

    @test norm(X - Xinflation)<1e-14
end

@testset "Additive inflation I" begin
    Ny = 10
    Nx= 10
    Ne = 20
    A = AdditiveInflation(Nx)

    @test A.Nx==Nx
    @test mean(A) == zeros(Nx)
    @test cov(A) == Diagonal(ones(Nx))

    X = zeros(Ny + Nx, Ne)
    X .= randn(Ny + Nx, Ne)

    Xadd = deepcopy(X)
    A(Xadd, Ny+1, Ny+Nx)
    @test size(Xadd) == (Ny + Nx, Ne)
    @test Xadd != X
    @test norm(Xadd[1:Ny,:] -  X[1:Ny,:]) < 10*eps()
    @test norm(Xadd[Ny+1:Ny+Nx,:] -  X[Ny+1:Ny+Nx,:]) > 10*eps()

    Σ = zeros(Nx,Nx)

    for i=1:Nx
        for j=1:Nx
            Σ[i,j] = 1.0/(i + j)
        end
    end

    A = AdditiveInflation(Nx, zeros(Nx), Σ)

    @test A.Nx==10
    @test mean(A) == zeros(Nx)
    @test cov(A) == Σ
end

@testset "Additive inflation II" begin
    Nx= 10
    Ne = 20
    A = AdditiveInflation(Nx)

    @test A.Nx==10
    @test mean(A) == zeros(Nx)
    @test cov(A) == Diagonal(ones(Nx))

    X = zeros(Nx, Ne)
    X .= randn(Nx, Ne)

    Xadd = deepcopy(X)
    A(Xadd)
    @test size(Xadd) == (Nx, Ne)
    @test Xadd != X

    Σ = zeros(Nx,Nx)

    for i=1:Nx
        for j=1:Nx
            Σ[i,j] = 1.0/(i + j)
        end
    end

    A = AdditiveInflation(Nx, zeros(Nx), Σ)

    @test A.Nx==10
    @test mean(A) == zeros(Nx)
    @test cov(A) == Σ
end

@testset "Multiplicative inflation I" begin
    Ny = 10
    Nx= 10
    Ne = 20
    A = MultiplicativeInflation(1.2)

    @test A.β==1.2

    X = zeros(Ny + Nx, Ne)
    X .= randn(Ny + Nx, Ne)

    Xmul = deepcopy(X)
    A(Xmul, Ny+1, Ny+Nx)

    for i=1:Ne
        x̄ = deepcopy(mean(X, dims = 2))[Ny+1:Ny+Nx,1]
        @test norm(Xmul[1:Ny,i] - X[1:Ny,i])<1e-10
        @test norm(Xmul[Ny+1:Ny+Nx,i] - (x̄ + A.β*(X[Ny+1:Ny+Nx,i] - x̄)))<1e-10
    end
end

@testset "Multiplicative inflation II" begin
    Nx= 10
    Ne = 20
    A = MultiplicativeInflation(1.2)

    @test A.β==1.2

    X = zeros(Nx, Ne)
    X .= randn(Nx, Ne)

    Xmul = deepcopy(X)
    A(Xmul)

    for i=1:Ne
        x̄ = deepcopy(mean(X, dims = 2))[:,1]
        @test norm(Xmul[:,i] - (x̄ + A.β*(X[:,i] - x̄)))<1e-10
    end
end

@testset "MultiAdditive Inflation" begin

    T = MultiAddInflation(10)

    @test size(T)== 10
    @test T.β == 1.0

    @test mean(T) == zeros(10)
    @test cov(T) ==  Diagonal(ones(10))

    T = MultiAddInflation(10, 1.01, zeros(10), 1.01)

    @test T.β == 1.01
    @test mean(T) == zeros(10)
    @test cov(T) == Diagonal(1.01^2*ones(10))
end


text:
module GmshTools

using MLStyle

const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("libgmsh not installed properly, run Pkg.build(\"GmshTools\"), restart Julia and try again")
end
include(depsjl_path)

gmshmodule = joinpath(dirname(libgmsh), "gmsh.jl")
include(gmshmodule)

function __init__()
    check_deps()
end

export gmsh, @gmsh_do, @gmsh_open

macro gmsh_do(f)
    quote
        try
            gmsh.initialize()
            $(esc(f))
        finally
            gmsh.finalize()
        end
    end
end

macro gmsh_open(name, f)
    quote
        try
            gmsh.initialize()
            gmsh.open($(esc(name)))
            $(esc(f))
        finally
            gmsh.finalize()
        end
    end
end

include("macros.jl")

end


text:
function debug(io,pkg::String)
    info("Reading build script...")
    dir = Pkg2.dir(pkg)
    file = joinpath(dir,"deps/build.jl")
    context = BinDeps.PackageContext(false,dir,pkg,{})
    m = Module(:__anon__)
    body = Expr(:toplevel,:(ARGS=[$context]),:(include($file)))
    eval(m,body)
    println(io,"The package declares $(length(context.deps)) dependencies.")
    for dep in context.deps
        print(io," - Library \"$(dep.name)\"")
        if !applicable(dep)
            println(io," (not applicable to this system)")
        else
            lib = _find_library(dep)
            if !isempty(lib)
                println(io," (satisfied by $lib)")
            else
                println(io)
            end
            if length(dep.providers) > 0
                println(io,"    - Providers:")
                for (dep,opts) in dep.providers
                    show_indented(io,dep,6)
                    println(io)
                end
            end
        end
    end
end
debug(pkg::String) = debug(STDOUT,pkg)


text:
using Base.Iterators: partition
using Flux
using Flux.Optimise: update!
using Flux: binarycrossentropy, mse, logitbinarycrossentropy
using Images
using MLDatasets
using Statistics, StatsBase
using Parameters: @with_kw
using Random
using Printf
using CUDAapi
using Zygote
using Augmentor
if has_cuda()		# Check if CUDA is available
    @info "CUDA is on"
    import CuArrays		# If CUDA is available, import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw struct HyperParams
    batch_size::Int = 128
    latent_dim::Int = 100
    nclasses::Int = 10
    epochs::Int = 30
    verbose_freq::Int = 500
    output_x::Int = 4        # No. of sample images to concatenate along x-axis
    output_y::Int = 3        # No. of sample images to concatenate along y-axis
    lr_dscr::Float64 = 0.0001
    lr_gen::Float64 = 0.0001
end

function load_data(hparams)
    # Load MNIST dataset
    images, labels = MLDatasets.MNIST.traindata(Float32)
    # Normalize to [-1, 1] and convert it to WHCN
    image_tensor = permutedims(reshape(@.(2f0 * images - 1f0), 28, 28, 1, :), (2, 1, 3, 4))
    img_rotated = similar(image_tensor)
    n = size(image_tensor,4)
    for i in 1:n
        img_rotated[:,:,1,i] = augment(image_tensor[:,:,1,i], Rotate90())
    end
    img_tensor = image_tensor[:,:,:,1:Int(n/2)]
    img_rotated = img_rotated[:,:,:,Int(n/2)+1:end]
    img_tensor = reshape(img_tensor, :, size(img_tensor,4))
    img_rotated = reshape(img_rotated, :, size(img_rotated,4))
    data = [(img_tensor[:,r], img_rotated[:,r]) |> gpu for r in partition(1:30000, hparams.batch_size)]
    return data
end

mutable struct gan_model
    m_common
    m1
    m2
end

function generator(args)
    m_common = Chain(Dense(args.latent_dim, 256),x-> leakyrelu.(x,0.2f0), BatchNorm(256), Dense(256, 512), x-> leakyrelu.(x,0.2f0), BatchNorm(512)) |> gpu
    m1 = Chain(Dense(512, 1024),x-> leakyrelu.(x,0.2f0), BatchNorm(1024),Dense(1024, 784, tanh)) |> gpu
    m2 = Chain(Dense(512, 1024),x-> leakyrelu.(x,0.2f0), BatchNorm(1024),Dense(1024, 784, tanh)) |> gpu
    gan_model(m_common,m1, m2)
end

function discriminator(args)
    m1 = Chain(Dense(784, 512),x-> leakyrelu.(x,0.2f0), Dense(512, 256),x-> leakyrelu.(x,0.2f0)) |> gpu
    m2 = Chain(Dense(784, 512),x-> leakyrelu.(x,0.2f0), Dense(512, 256),x-> leakyrelu.(x,0.2f0)) |> gpu
    m_common = Chain(Dense(256,1, sigmoid)) |> gpu
    gan_model(m_common, m1, m2)
end

# Loss functions
function discr_loss(real_output, fake_output)
    real_loss = mean(binarycrossentropy.(real_output, 1f0))
    fake_loss = mean(binarycrossentropy.(fake_output, 0f0))
    return real_loss + fake_loss
end

generator_loss(fake_output) = mean(binarycrossentropy.(fake_output, 1f0))

function train_discr(discr, img1, img2, fake_img1, fake_img2, opt_discr)
    ps_d = Flux.params(discr.m_common, discr.m1, discr.m2)
    loss = 0.0
    loss, back = Zygote.pullback(ps_d) do
                    x1_real = discr.m_common(discr.m1(img1))
                    x2_real = discr.m_common(discr.m2(img2))
                    x1_fake = discr.m_common(discr.m1(fake_img1))
                    x2_fake = discr.m_common(discr.m2(fake_img2))
                    discr_loss(x1_real, x1_fake) + discr_loss(x2_real, x2_fake)
    end
    grads = back(1f0)
    for i in params(discr.m_common)
        grads.grads[i] ./=2
    end
    update!(opt_discr, ps_d, grads)
    return loss
end

Zygote.@nograd train_discr

function train_gan(gen, discr, img1, img2, opt_gen, opt_discr, hparams)
    noise = randn!(similar(img1, (hparams.latent_dim, hparams.batch_size))) |> gpu
    ps = Flux.params(gen.m_common, gen.m1, gen.m2)
    loss = Dict()
    loss["gen"], back = Zygote.pullback(ps) do
                            x = gen.m_common(noise)
                            fake_1 = gen.m1(x)
                            fake_2 = gen.m2(x)
                            loss["discr"] = train_discr(discr, img1, img2, fake_1, fake_2, opt_discr)
                            generator_loss(discr.m_common(discr.m1(fake_1))) + generator_loss(discr.m_common(discr.m2(fake_2)))
    end
    grads = back(1f0)
    for i in params(gen.m_common)
         grads.grads[i] ./=2
    end
    update!(opt_gen, ps, grads)
    return loss
end

function create_output_image(gen, fixed_noise, hparams)
    @eval Flux.istraining() = false
    fake_i = @. cpu(gen.m1(gen.m_common(fixed_noise)))
    fake_i = @. reshape(fake_i, 28,28,1,size(fake_i,2))
    fake_i_2 = @. cpu(gen.m2(gen.m_common(fixed_noise)))
    fake_i_2 = @. reshape(fake_i_2, 28,28,1,size(fake_i_2,2))
    fake_images =@. hcat(fake_i, fake_i_2)
    @eval Flux.istraining() = true
    image_array = dropdims(reduce(vcat, reduce.(hcat, partition(fake_images, hparams.output_y))); dims=(3, 4))
    image_array = @. Gray(image_array + 1f0) / 2f0
    return image_array
end

function train()
    hparams = HyperParams()

    data = load_data(hparams)

    fixed_noise = [randn(hparams.latent_dim, 1) |> gpu for _=1:hparams.output_x*hparams.output_y]

    # Discriminator
    dscr = discriminator(hparams)

    # Generator
    gen =  generator(hparams)

    # Optimizers
    opt_dscr = ADAM(hparams.lr_dscr, (0.5,0.99))
    opt_gen = ADAM(hparams.lr_gen, (0.5,0.99))

    isdir("output")||mkdir("output")

    # Training
    train_steps = 0
    for ep in 1:hparams.epochs
        @info "Epoch $ep"
        for (img1, img2) in data
            loss = train_gan(gen, dscr, img1, img2, opt_gen, opt_dscr, hparams)

            if train_steps % hparams.verbose_freq == 0
                @info("Train step $(train_steps), Discriminator loss = $(loss["discr"]), Generator loss = $(loss["gen"])")
                # Save generated fake image
                output_image = create_output_image(gen, fixed_noise, hparams)
                save(@sprintf("output/cogan_steps_%06d.png", train_steps), output_image)
            end
            train_steps += 1
        end
    end

    output_image = create_output_image(gen, fixed_noise, hparams)
    save(@sprintf("output/cogan_steps_%06d.png", train_steps), output_image)
end

cd(@__DIR__)
train()


text:
"""
    dummy_matrix(df::DataFrame, column)

Takes a DataFrame and column name and expands factors
into dummy columns - helps to convert columns to "colname_colvalue" since values
are used as column names in the dummy DataFrame.
"""
function dummy_matrix(df::DataFrame, column)
    colnames = permutedims(df[!, column])
    dummy_m = unique(df[!, column]) .== colnames
    dummy_df = DataFrame(Matrix(permutedims(dummy_m)), convert(Vector{String}, string.(unique(colnames[:]))))
    return dummy_df
end

"""
    tidy(fit::LinearModelFit)
Inspired by R's broom::tidy. Returns a dataframe where each row is a variable
and its associated estimate.
"""
function tidy(fit::T) where {T <: Fit}
    X_names = fit.X_names
    β = fit.modelfit.β
    se_β = fit.modelfit.se_β
    pval = fit.modelfit.pval
    df = DataFrame(
        term = X_names,
        β = β,
        se_β = se_β,
        t_stat = β ./ se_β,
        p_val = pval
        )
    return df
end

"""
    tidy(model_list::Vector{LinearModelFit}, model_name::Vector{String})
Method dispatch over a list of models and a vector of model identifiers
so we can just tidy over multiple broadcasted models.
"""
function tidy(model_list::Vector{T} where {T <: Fit}, model_names::Vector)
    length(model_list) == length(model_names) || error("Model - Name Dimension Mismatch")
    tidy_model_list = Vector{DataFrame}(undef, length(model_list))
    for i in 1:length(model_list)
        df = tidy(model_list[i])
        df[!, "model"] .= model_names[i]
        tidy_model_list[i] =  df
    end
    tidy_df = reduce(vcat, tidy_model_list)
    return tidy_df
end
"""
    tidy(model_list::Vector{T} where {T <: Fit})
Tidy a list of models and generate model_i names automatically.
"""
function tidy(model_list::Vector{T} where {T <: Fit})
    tidy_model_list = Vector{DataFrame}(undef, length(model_list))
    for i in 1:length(model_list)
        df = tidy(model_list[i])
        df[!, "model"] .= "model_" * string(i)
        tidy_model_list[i] =  df
    end
    tidy_df = reduce(vcat, tidy_model_list)
    return tidy_df
end

function TableCol(header,
                  m::LinearModelFit;
                  drop = [""],
                  keep = [""],
                  stats=(:N=>Int∘nobs, "\$R^2\$"=>r2),
                  meta=())
    # Initialize the column
    col  = RegCol(header)
    keep_used = !(keep == [""])
    # if we use keep, add everything to drop and remove elements in keep
    if keep_used
        drop = deepcopy(m.X_names)
        drop = filter(e -> !(e in keep), drop)
    end
    # Add the coefficients
    for (name, val, se, p) in zip(m.X_names, coef(m), m.modelfit.se_β, m.modelfit.pval)
        if !(name in drop)
            setcoef!(col, name, val, se)
            0.05 <  p <= .1  && star!(col[name], 1)
            0.01 <  p <= .05 && star!(col[name], 2)
                    p <= .01 && star!(col[name], 3)
        end
    end

    # Add in the fit statistics
    setstats!(col, OrderedDict(p.first=>p.second(m) for p in stats))

    # Add in the metadata
    setmeta!(col, OrderedDict(p.first=>p.second(m) for p in meta))

    return col

end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
using Statistics
using BlackBoxOptim
using NLopt
using Distributions

## Optimization
# Generate random index
indfast=rand(1:Param.nsim, Param.nfast)
indfast[1]=1 # Keeping the first element fixed.

# schain=SharedArray{Float64,2}(mchain[indfast,:])
mchain_op=mchain[indfast,:]

## Objective function

function objfun(gamma0...; Par::Par=Param, mchain=mchain_op)
    @unpack nsim, nfast = Par
    gamma = gamma0
    N = size(mchain[1,:])
    val=0
    logunif=log.(rand(nfast))
    geta=mchain[1,:]
    hM=zeros(N)
    for r=1:nfast
        gtry=mchain[r,:]
        for j=1:size(gtry,2)
            val=+gtry[:,j].*gamma[j]-geta[:,j].*gamma[j]
        end
        geta=logunif[r] < val[1] ? gtry : geta
        hM=+geta./nfast
    end
    mhM=mean(hM)
    var=Statistics.var(hM)

    obj=(mhM^2)/var

    return obj
end

function objfun_bbo(gamma0; Par::Par=Param, mchain=mchain_op)
    @unpack nsim, nfast = Par
    gamma = gamma0
    N = size(mchain[1,:])
    val=0
    logunif=log.(rand(nfast))
    geta=mchain[1,:]
    hM=zeros(N)
    for r=1:nfast
        gtry=mchain[r,:]
        for j=1:size(gtry,2)
            val=+gtry[:,j].*gamma[j]-geta[:,j].*gamma[j]
        end
        geta=logunif[r] < val[1] ? gtry : geta
        hM=+geta./nfast
    end
    mhM=mean(hM)
    var=Statistics.var(hM)

    obj=(mhM^2)/var

    return obj
end
## testing obj fun

objfun(collect(1:6))

objfun_bbo(repeat([200],6))
objfun_bbo(collect(1:6))
#
fop=nothing
## Solvin BBO
res = bboptimize(x->objfun_bbo(x); Sea = (-10e300,10e300), NumDimensions = 6)

## Getting results

minbb = best_fitness(res)
γ0 = best_candidate(res)
TS_bb = minbb*6
p_val_bb = ccdf(Chisq(4),TS_bb)
## Optimizing

fop=nothing
fop=Model(NLopt.Optimizer)
set_optimizer_attribute(fop, "algorithm", :LN_BOBYQA)
set_optimizer_attribute(fop, "xtol_rel", 1e-8)
@variable(fop, gamma0[1:6])

register(fop, :objfun, 6, objfun; autodiff=true )
@NLobjective(fop, Min, objfun(gamma0...))

for i=1:6
    set_start_value(gamma0[i],γ0[i])
end

## Solving in NLopt

JuMP.optimize!(fop)

## Getting Results

termination_status(fop)
primal_status(fop)

value.(gamma0)
obj_v=objective_value(fop)
# n=size(q,2)
TS=6*obj_v
p_val = ccdf(Chisq(4),TS)

## Testing function

function RP_Cournot_test(;Par::Par=Param, mchain=mchain_op)
    N = size(mchain,2)
    # First optimization BlackBox
    res = bboptimize(x->objfun_bbo(x); Sea = (-10e300,10e300), NumDimensions = 6)
    ## Getting results
    minbb = best_fitness(res)
    γ0 = best_candidate(res)
    # TS_bb = minbb*N
    # p_val_bb = ccdf(Chisq(N),TS_bb)

    ## Second optimization NLopt
    fop=nothing
    fop=Model(NLopt.Optimizer)
    set_optimizer_attribute(fop, "algorithm", :LN_BOBYQA)
    set_optimizer_attribute(fop, "xtol_rel", 1e-8)
    @variable(fop, gamma0[1:N])

    register(fop, :objfun, N, objfun; autodiff=true )
    @NLobjective(fop, Min, objfun(gamma0...))

    for i=1:N
        set_start_value(gamma0[i],γ0[i])
    end

    ## Solving in NLopt
    JuMP.optimize!(fop)

    ## Getting Results
    termination_status(fop)
    primal_status(fop)

    γ1=value.(gamma0)
    obj_v=objective_value(fop)
    TS=N*obj_v
    p_val = ccdf(Chisq(N),TS)

    ## Second optimization for refinement
    for i=1:N
        set_start_value(gamma0[i],γ1[i])
    end
    ## Solving in NLopt
    JuMP.optimize!(fop)
    obj_v=objective_value(fop)
    TS=N*obj_v
    p_val = ccdf(Chisq(N),TS)

    fop=nothing

    return TS, p_val
end

max_stars_repo_name: hans-mtz/SCR

max_issues_repo_name: hans-mtz/SCR

max_forks_repo_name: hans-mtz/SCR


text:
export trees2graph

"""
	trees2graph(t::Vararg{TreeTools.Tree})

Create a `Graph` object `g` from a collection of trees. All trees should have the same leaf nodes.
1. An array of `LeafNode` from the leaves of trees
2. Add internal nodes of each tree `t` to `g` calling `tree2graph!(g,t,k)`, where `k` is the position of `t` in the list
"""
function trees2graph(t::Vararg{TreeTools.Tree})
	treelist = t
	# checklabels(treelist)

	# labels = collect(keys(first(treelist).lleaves))
	# labels_to_int = Dict(k=>findfirst(==(k), labels) for k in labels)
	K = length(treelist)
	labels = collect(keys(first(treelist).lleaves))
	labels_to_int = Dict(l=>k for (k,l) in enumerate(labels))
	leaves = [
		LeafNode(;
			index=i,
			conf=Int[i],
			anc=Array{SplitNode,1}(undef,K)
		) for i in 1:length(labels)
	]
	lleaves = Dict(label => leaf for (leaf, label) in zip(leaves, labels))
	g = Graph(;
		labels=labels,
		labels_to_int=labels_to_int,
		leaves=leaves,
		lleaves=lleaves,
		K=K,
	)
	for (k,t) in enumerate(treelist)
		tree2graph!(g, t, k) # Adding tree t to color k of graph
	end

	return g
end

"""
"""
function trees2graph(treelist)
	return trees2graph(treelist...)
end

"""
	tree2graph!(g::Graph, t::TreeTools.Tree, k::Int64)

Add internal nodes of `t` to `Graph` `g`. `g` should already have fields `leaves`, `labels` and `labels_to_int` initialized.
Call `tree2splitnodes(g, t.root, SplitNode(), k)` to start the process.
"""
function tree2graph!(g::Graph, t::TreeTools.Tree, k::Int)
	if length(t.lleaves) != length(g.leaves)
		@error("`t.lleaves` and `g.leaves` do not have the same length.")
	end
	tree2splitnodes!(g, t.root, SplitNode(), k)
end

"""
	tree2splitnodes!(g::Graph, r::TreeTools.TreeNode, sr::SplitNode, k::Int)

Add children of `TreeNode` `r` to internal nodes of `g`. The `SplitNode` object `sr` corresponding to `r` should exist, except if `r` is root.
Recursively call `tree2splitnodes!` on children of `r`.
"""
function tree2splitnodes!(g::Graph, r::TreeTools.TreeNode, sr::SplitNode, k::Int)
	# If r is root, create a custom `sr`. Otherwise, use input one.
	if r.isroot
		sr = SplitNode(;
			anc = nothing,
			child = Array{GraphNode,1}(undef, length(r.child)),
			color = k,
			conf = treenode2conf(g,r),
			isroot = true,
		)
		push!(g.internals, sr)
	end
	# Special case of single leaf tree
	if r.isleaf
		push!(sr.child, g.lleaves[r.label])
		g.lleaves[r.label].anc[k] = sr
	end
	#
	for (i,c) in enumerate(r.child)
		if !c.isleaf
			sc = SplitNode(;
				anc = sr,
				child = Array{GraphNode,1}(undef, length(c.child)),
				color = k,
				conf = treenode2conf(g,c),
				isroot = false
			)
			tree2splitnodes!(g, c, sc, k)
			sr.child[i] = sc
			push!(g.internals, sc)
		else
			sr.child[i] = g.lleaves[c.label]
			if isassigned(g.lleaves[c.label].anc, k)
				error("Leaf $(c.label) has ancestor of color $k already assigned")
			end
			g.lleaves[c.label].anc[k] = sr
		end
	end
end

"""
COULD BE MADE FASTER WITH A MESSAGE PASSING ALG
"""
treenode2conf(g::Graph, n) = Int[g.labels_to_int[x.label] for x in POTleaves(n)]

"""
"""
function checklabels(treelist)
	t1 = first(treelist)
	for t in treelist
		if !TreeTools.share_labels(t1,t)
			@error("Trees do not share leaves nodes labels")
		end
	end
	return nothing
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
# This file was generated, do not modify it. # hide
fit!(ŷ, rows=train)
ypreds = ŷ(rows=test)
rmsl(y[test], ypreds)


text:
"""
    InverseNormalisingFlow{T<:Vector{<:Invertible}}

A normalising flow that operates in reverse. It is thus efficient to compute the log
density of an observation, but not (necessarily) efficent to compute the
"""
struct InverseNormalisingFlow{T<:Vector{<:Invertible}} <: Flow
    p0::Any
    transforms::T
    function InverseNormalisingFlow(p0, transforms::T) where T<:Vector{<:Invertible}
        @assert all(dim(p0) .== dim.(transforms))
        return new{T}(p0, transforms)
    end
end
const INF = InverseNormalisingFlow

"""
    rand(rng::AbstractRNG, d::INF)

Sample from a INF `d`. Will only work if each transform in `d` has its
`invert` function defined.
"""
function rand(rng::AbstractRNG, d::INF, N::Int=1)
    y = rand(rng, d.p0, N)
    for f in reverse(d.transforms)
        y = invert(f, y)
    end
    return y
end

"""
    logpdf(d::INF, y::RealOrVecOrMat)

Compute the log pdf of an observation `y` of the INF `d`.
"""
function logpdf(d::INF, y::AbstractVecOrMat{<:Real})
    l = zero(eltype(y))
    for k in 1:length(d.transforms)
        l += sum(logdetJ(d.transforms[k], y))
        y = apply(d.transforms[k], y)
    end
    return l + sum(logpdf(d.p0, y))
end

for init in [:naive_init!, :identity_init!]
    @eval begin
    """
        naive_init!(
            θ::AbstractVector{<:Real},
            ::Type{<:INF},
            rng::AbstactRNG,
            D::Int,
            ctors::AbstractVector,
        )

    Fill `θ` with initialisations appropriate for each component of an INF using only
    information available prior to construction.
    """
    function $init(
        θ::AbstractVector{<:Real},
        ::Type{<:INF},
        rng::AbstractRNG,
        D::Int,
        ctors::AbstractVector,
    )
        @assert all(issubtype.(ctors, Invertible))
        @assert length(θ) == nparams(INF, D, ctors)
        pos = 1
        for p in eachindex(ctors)
            Δ = nparams(ctors[p], D)
            $init(view(θ, pos:pos+Δ-1), ctors[p], rng, D)
            pos += Δ
        end
        return θ
    end
    end
end

"""
    INF(::typeof(naive_init), rng::AbstractRNG, D::Int, ctors::AbstractVector)

Initialise an `InverseNormalisingFlow` with the `naive_init`s for each `Invertible`, and
a `DiagonalStandardNormal` as the base distribution.
"""
INF(::typeof(naive_init), rng::AbstractRNG, D::Int, ctors::AbstractVector) =
    INF(DiagStdNormal(D), naive_init(INF, rng, D, ctors), ctors)

"""
    INF(::typeof(identity_init), rng::AbstractRNG, D::Int, ctors::AbstractVector)

Initialise an `InverseNormalisingFlow` with the `identity_init`s for each `Invertible`, and
a `DiagonalStandardNormal` as the base distribution.
"""
INF(::typeof(identity_init), rng::AbstractRNG, D::Int, ctors::AbstractVector) =
    INF(DiagStdNormal(D), identity_init(INF, rng, D, ctors), ctors)


text:
###{{{ phreg

function phreg(X::Matrix, t::Vector, status::Vector, entry::Vector, id=[]; beta=[], opt...)
    pre = EventHistory.coxprep(t,status,X,entry,id)
    p = size(X,2)
    if size(beta,1)==0 beta=vec(zeros(p,1)) end
    pl(beta,indiv=false) =
        EventHistory.coxPL(beta,pre["X"],pre["XX"],pre["sign"],pre["jumps"],indiv)
    val = EventHistory.NR(pl,beta; opt...)
    beta = vec(val[1]); grad = val[2].gradient; steps = val[3];
    U = pl(beta,true)
    I = -pinv(U[3])
    iid = U[2]*I
    V =  iid'iid
    coefmat = [beta sqrt(diag(V)) sqrt(diag(I))]
    coefmat = [coefmat 2*Distributions.cdf(Distributions.Normal(0,1),-abs(coefmat[:,1]./coefmat[:,2]))]
    ## coln = [:Estimate,:SE,:naiveSE,symbol("pvalue")]
    ## cc = convert(DataFrame,coefmat); names!(cc,coln);
    cc = CoefTable(coefmat,["Estimate","S.E","naive S.E.","P-value"],
                   repeat([""],inner=[size(coefmat,1)]))
    chaz = [pre["jumps"] pre["time"][pre["jumps"]] cumsum(1./(U[4][pre["jumps"]]))]
    EventHistoryModel("Cox",Formula(:.,:.),EventHistory.Surv,
                      vec(beta),cc,
                      iid,I,V,
                      [steps],vec(grad),X,[t status],
                      chaz
                      )
end

function phreg(formula::Formula, data::DataFrame; id=[], opt...)
    ## fm = DataFrames.Formula(formula.args[2],formula.args[3])
    M = DataFrames.ModelFrame(formula,data)
    S = M.df[:,1] ## convert(Vector,M.df[:,1])
    time = EventHistory.Time(S)
    status = EventHistory.Status(S)
    entry = try EventHistory.Entry(S) catch [] end
    X = DataFrames.ModelMatrix(M)
    X = X.m[:,collect(2:size(X.m,2))]
    res = EventHistory.phreg(X, time, status, entry; opt...)
    cnames = setdiff(coefnames(M),["(Intercept)"])
    res.coefmat.rownms=cnames
    res.formula = formula
    res.eventtype = typeof(S[1])
    res
end

###}}} phreg

###{{{ coxprep

function coxprep(exit,status,X=[],entry=[],id=[])
    Truncation = size(entry,1)==size(exit,1)
    n = size(exit,1)
    p = size(X,2)
    truncation = size(entry,1)==n
    XX = Array(Number, n, p*p) ## Calculate XX' at each time-point
    ##  XX = Array(Number, p, p, n)
    if size(X,1)>0
        for i=1:n
            Xi = X[i,:]'
            XX[i,:] = vec(Xi'*Xi)
            ## XX[:,:,i] = Xi'Xi
        end
    end
    sgn=[]
    if truncation
        sgn = Array(Int,2*n); fill!(sgn,1)
        status = [status;status]
        for i=1:n
            sgn[i] = -1
            status[i] = 0
        end
        exit = [entry;exit]
        X = [X;X]
        XX = [XX;XX]
    end
    ord0 = sortperm(status*1,rev=true)
    ord = sortperm(exit[ord0])
    ord = ord0[ord]
    if (truncation)
        sgn = sgn[ord]
    end
    if size(X,1)>0
        X = X[ord,:]
        XX = XX[ord,:]
    end
    exit = exit[ord]
    status = status[ord]
    jumps = find(status)
    Dict("X"=>X, "XX"=>XX, "jumps"=>jumps, "sign"=>sgn,
     "ord"=>ord, "time"=>exit, "id"=>[])
end

###}}} coxprep

###{{{ revcumsum
function revcumsum(A,dim=1)
    D = size(A)
    n = D[dim]
    res = similar(A)
    prev = zeros(Number,1,size(A,2))
    for i=1:n
        idx =
        prev += A[n-i+1,:]
        res[n-i+1,:] = prev
    end
    res
end
###}}} revcumsum

###{{{ coxPL
function coxPL(beta::Vector, X::Matrix, XX::Matrix, sgn::Vector, jumps::Vector, indiv=false)
    n = size(X,1)
    p = size(X,2)
    Xb = X*beta
    eXb = map(exp,Xb)
    if size(sgn,1)==n ## Truncation
        eXb = sgn.*eXb
    end
    S0 = revcumsum(eXb)
    E = similar(X)
    D2 = similar(XX)
    S1 = Array(Number,n,p)
    S2 = similar(XX)
    for j=1:p
        if (indiv) S1[:,j] = revcumsum(X[:,j].*eXb); end
        E[:,j] = revcumsum(X[:,j].*eXb)./S0 ## S1/S0(s)
    end
    for j=1:size(XX,2)
        if (indiv) S2[:,j] = revcumsum(XX[:,j].*eXb); end
        D2[:,j] = revcumsum(XX[:,j].* eXb)./S0 ## int S2/S0(s)
    end
    D2 = D2[jumps,:]
    E = E[jumps,:]
    grad = (X[jumps,:]-E) ## Score
    val = Xb[jumps]-log(S0[jumps]) ## Partial log-likelihood
    hess = -(reshape(sum(D2,1),p,p)-E'E)
    if indiv
        return(val,grad,hess,vec(S0),S1,S2,E)
    end
    val = sum(val); grad = sum(grad,1)
    EventHistory.D2Function(val[1],vec(grad),hess)
end

###}}} coxPL

###{{{ predict
# compute the values of a step function,
# i.e. how many of the jumps are smaller or
# equal to the eval points
function sindex(jump::Vector, eval::Vector)
    N = size(eval,1)
    n = size(jump,1)
    index = Array(Int64,N)
    j = 1::Int64
    for t=1:N
        while (j<n) & (jump[j]<eval[t])
            j += 1
        end
        index[t] = j
    end
    return(index)
end

## Returns [time surv-prob]
function predict(mm::EventHistory.EventHistoryModel; X=[]'::Matrix,time=mm.eventtime[:,1]::Vector,surv=true,order=false)
    ord = sortperm(time)
    idx = EventHistory.sindex(mm.chaz[:,2],time[ord,1])
    L0 = mm.chaz[:,3][idx]
    L0 = L0[sortperm(ord)]
    arraytype = size(X,1)==1
    if size(X,2)==0
        X = mm.X
        arraytype = true
    end
    if arraytype
        H = exp(X*coef(mm))
        res = L0.*H
    else
        res = Array(Number,size(time,1),size(X,1))
        for i=1:size(X,1)
            res[:,i] = L0.*exp(X[i,:]'coef(mm))
        end
    end
    if surv res = exp(-res) end
    res = convert(Array{Number},[time res])
    if order
        return(res[sortperm(time),:])
    end
    return(res)
end

###}}} predict

function coef(mm::EventHistoryModel)
    mm.coef
end

function vcov(mm::EventHistoryModel)
    mm.vcov
end

function coeftable(mm::EventHistoryModel)
    mm.coefmat

end


text:
export COC,WL,WR,SL,SR, stateType, actType, ACTIONS, discount_f, RANGES,THETAS,PSIS,OWNSPEEDS,INTRPSEEDS, interp, turns

# ADVISORY INDICES
COC=0
WL=1
WR=2
SL=3
SR=4

# State Type:
stateType = Tuple{Float64,Float64,Float64,Float64,Float64,Int}
actType = Int
ACTIONS = [COC,WL,WR,SL,SR]

# Default parameters
discount_f = 1.0

### STATE CUTPOINTS ###
RANGES = [0.0,25.0,50.0,75.0,100.0,150.0,200.0,300.0,400.0,500.0,510.0,750.0,1000.0,1500.0,2000.0,3000.0,4000.0,5000.0,7000.0,9000.0,11000.0,13000.0,15000.0,17000.0,19000.0,21000.0,25000.0,30000.0,35000.0,40000.0,48000.0,56000.0]
THETAS = Array(LinRange(-pi,pi,41))
PSIS   = Array(LinRange(-pi,pi,41))
OWNSPEEDS = [200.0] #[100.0,200.0,300.0,400.0]  #ft/s
INTRSPEEDS = [200.0] #[100.0,200.0,300.0,400.0]  #ft/s

interp = LocalGIFunctionApproximator(RectangleGrid(RANGES,THETAS,PSIS,OWNSPEEDS,INTRSPEEDS,ACTIONS)) # Create the local function approximator using the grid

### Dictionaries to define transitions ###
probs = [0.5,0.25,0.25]
turns = Dict(COC=>([0.34,0.33,0.33],[0.0,1.5,-1.5].* pi/180.0 ),
              WL=>(probs,[1.5,2.0,1.25].* pi/180.0 ),
              WR=>(probs,[-1.5,-1.25,-2.0].* pi/180.0 ),
              SL=>(probs,[3.0,4.0,2.0].* pi/180.0 ),
              SR=>(probs,[-3.0,-2.0,-4.0].* pi/180.0 ),
              -1=>([0.34,0.33,0.33],[0.0,1.5,-1.5].* pi/180.0 )) # FOR v5, 0, 1, -1

max_stars_repo_name: cherrywoods/HorizontalCAS

max_issues_repo_name: ahmed-irfan/HorizontalCAS

max_forks_repo_name: ahmed-irfan/HorizontalCAS


text:
using InstantiateFromURL
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.3.0")

max_stars_repo_name: yingxiangli/lecture-source-jl

max_issues_repo_name: yingxiangli/lecture-source-jl

max_forks_repo_name: yingxiangli/lecture-source-jl


text:
#### RefState

using OrdinaryDiffEq

"""
    init_ref_state!(aux::StateVec,
                    grid::Grid,
                    params,
                    output_dir::String)

Initializes the reference state variables:
  - `p_0` pressure
  - `ρ_0` density
  - `α_0` specific volume
assuming a hydrostatic balance.
FIXME: add reference
"""
function init_ref_state!(
    aux::StateVec,
    grid::Grid{FT},
    params,
    output_dir::String,
) where {FT}
    @unpack param_set, SurfaceModel = params
    T_g = SurfaceModel.T
    q_tot_g = SurfaceModel.q_tot
    P_g = SurfaceModel.P

    q_pt_g = PhasePartition(q_tot_g)
    θ_liq_ice_g =
        TD.liquid_ice_pottemp_given_pressure(param_set, T_g, P_g, q_pt_g)
    logp = log(P_g)

    function tendencies(logp, u, z)
        p = exp(logp)
        ts = PhaseEquil_pTq(param_set, p, T_g, q_tot_g)
        R_m = gas_constant_air(ts)
        T = air_temperature(ts)
        return -FT(grav(param_set)) / (T * R_m)
    end

    z_span = (grid.zn_min, grid.zn_max)
    prob = ODEProblem(tendencies, logp, z_span)
    sol = solve(prob, Tsit5(), reltol = 1e-12, abstol = 1e-12)
    p_0 = [exp(sol(grid.zc[k])) for k in over_elems_real(grid)]
    assign_real!(aux, :p_0, grid, p_0)
    apply_Neumann!(aux, :p_0, grid, 0.0, Zmin())
    apply_Neumann!(aux, :p_0, grid, 0.0, Zmax())

    @inbounds for k in over_elems(grid)
        ts = PhaseEquil_pTq(param_set, aux[:p_0, k], T_g, q_tot_g)
        q_pt = PhasePartition(ts)
        T = air_temperature(ts)
        aux[:ρ_0, k] = air_density(param_set, T, aux[:p_0, k], q_pt)
        aux[:α_0, k] = 1 / aux[:ρ_0, k]
    end
    extrap!(aux, :ρ_0, grid)
    extrap!(aux, :α_0, grid)
    extrap!(aux, :p_0, grid)
end


text:
module StateSpaceModels

import Base: show, length, isempty

using Distributions
using LinearAlgebra
using ShiftedArrays
using Statistics
using Polynomials
using MatrixEquations
using Printf
using Optim
using OrderedCollections
using RecipesBase
using SeasonalTrendLoess
using SparseArrays
using StatsBase

abstract type StateSpaceModel end

include("datasets.jl")

include("hyperparameters.jl")
include("systems.jl")
include("kalman_filter_and_smoother.jl")
include("statistical_tests/kpss.jl")
include("statistical_tests/canova_hensen.jl")
include("filters/univariate_kalman_filter.jl")
include("filters/multivariate_kalman_filter.jl")
include("filters/scalar_kalman_filter.jl")
include("filters/regression_kalman_filter.jl")
include("filters/sparse_univariate_kalman_filter.jl")

include("smoothers/kalman_smoother.jl")

include("optimizers.jl")
include("fit.jl")
include("prints.jl")
include("forecast.jl")
include("cross_validation.jl")

include("models/common.jl")
include("models/locallevel.jl")
include("models/locallevelcycle.jl")
include("models/locallevelexplanatory.jl")
include("models/locallineartrend.jl")
include("models/basicstructural.jl")
include("models/basicstructural_explanatory.jl")
include("models/basicstructural_multivariate.jl")
include("models/sarima.jl")
include("models/linear_regression.jl")
include("models/unobserved_components.jl")
include("models/exponential_smoothing.jl")
include("models/naive_models.jl")
include("models/dar.jl")

include("visualization/forecast.jl")
include("visualization/components.jl")
include("visualization/cross_validation.jl")
include("visualization/diagnostics.jl")

# Exported types and structs
export BasicStructural
export ExperimentalSeasonalNaive
export BasicStructuralExplanatory
export DAR
export ExponentialSmoothing
export FilterOutput
export LinearMultivariateTimeInvariant
export LinearMultivariateTimeVariant
export LinearRegression
export LinearUnivariateTimeInvariant
export LinearUnivariateTimeVariant
export LocalLevel
export LocalLevelCycle
export LocalLevelExplanatory
export LocalLinearTrend
export MultivariateBasicStructural
export Naive
export Optimizer
export SARIMA
export ScalarKalmanFilter
export SeasonalNaive
export SmootherOutput
export SparseUnivariateKalmanFilter
export StateSpaceModel
export UnivariateKalmanFilter
export UnobservedComponents

# Exported functions
export auto_arima
export auto_ets
export cross_validation
export constrain_box!
export constrain_identity!
export constrain_variance!
export fit!
export fix_hyperparameters!
export forecast
export forecast_expected_value
export get_constrained_value
export get_filtered_state
export get_filtered_state_variance
export get_hyperparameters
export get_innovations
export get_innovations_variance
export get_names
export get_predictive_state
export get_predictive_state_variance
export get_smoothed_state
export get_smoothed_state_variance
export has_fit_methods
export isfitted
export isunivariate
export kalman_filter
export kalman_smoother
export loglike
export num_states
export number_hyperparameters
export results
export set_initial_hyperparameters!
export simulate
export simulate_scenarios
export unconstrain_box!
export unconstrain_identity!
export unconstrain_variance!

end


text:
module TestVecnorm
using ParallelAccelerator

@acc norm_t(y) = norm(y,1)
@acc norm_t2(y) = norm(y,2)

function test()
    y = [1.,2.,3.]
    z = norm_t(y)
    return z
end

function test2()
    y = [1.,2.,3.]
    z = norm_t2(y)
    return z
end

end

using Base.Test
println("testing vecnorm...")
@test_approx_eq TestVecnorm.test() 6.0
@test_approx_eq TestVecnorm.test2() 3.7416573867739413
println("Done testing vecnorm.")


text:
# ----------------------------------------------------
# --- Mutators.jl
# ----------------------------------------------------
# Functions to update radio configuration and parameters
"""
Update carrier frequency of current radio device, and update radio object with the new obtained sampling frequency.
# --- Syntax
updateCarrierFreq!(radio,carrierFreq)
# --- Input parameters
- radio	  : SDR device
- carrierFreq	: New desired carrier frequency
# --- Output parameters
- carrierFreq : Effective carrier frequency
"""
updateCarrierFreq!(obj::SDROverNetwork,tul...) = SDROverNetworks.updateCarrierFreq!(obj,tul...);
updateCarrierFreq!(obj::UHDBinding,tul...) = UHDBindings.updateCarrierFreq!(obj,tul...);
updateCarrierFreq!(obj::RadioSim,tul...) = RadioSims.updateCarrierFreq!(obj,tul...);
updateCarrierFreq!(obj::RTLSDRBinding,tul...) = RTLSDRBindings.updateCarrierFreq!(obj,tul...);
function updateCarrierFreq!(obj::PlutoSDR,tul...)
    # In pluto we only get a flag so we need to call to the accessor
    AdalmPluto.updateCarrierFreq!(obj,_toInt.(tul)...);
    return getCarrierFreq(obj)
end

"""
Update sampling rate of current radio device, and update radio object with the new obtained sampling frequency.
# --- Syntax
updateSamplingRate!(radio,samplingRate)
# --- Input parameters
- radio	  : SDR device
- samplingRate	: New desired sampling rate
# --- Output parameters
- samplingRate : Effective sampling rate
"""
updateSamplingRate!(obj::SDROverNetwork,tul...) = SDROverNetworks.updateSamplingRate!(obj,tul...);
updateSamplingRate!(obj::UHDBinding,tul...) = UHDBindings.updateSamplingRate!(obj,tul...);
updateSamplingRate!(obj::RadioSim,tul...) = RadioSims.updateSamplingRate!(obj,tul...);
updateSamplingRate!(obj::RTLSDRBinding,tul...) = RTLSDRBindings.updateSamplingRate!(obj,tul...);
function updateSamplingRate!(obj::PlutoSDR,tul...)
    # For Adalm Pluto we should also update the RF filter band
    AdalmPluto.updateSamplingRate!(obj,_toInt.(tul)...);
    AdalmPluto.updateBandwidth!(obj,_toInt.(tul)...);
    return getSamplingRate(obj)
end

"""
Update gain of current radio device, and update radio object with the new obtained gain.
If the input is a [UHDRx] or a [UHDTx] object, it updates only the Rx or Tx gain
# --- Syntax
updateGain!(radio,gain)
# --- Input parameters
- radio	  : SDR device
- gain	: New desired gain
# --- Output parameters
- gain : New gain value
"""
updateGain!(obj::SDROverNetwork,tul...) = SDROverNetworks.updateGain!(obj,tul...);
updateGain!(obj::UHDBinding,tul...) = UHDBindings.updateGain!(obj,tul...);
updateGain!(obj::RadioSim,tul...) = RadioSims.updateGain!(obj,tul...);
updateGain!(obj::RTLSDRBinding,tul...) = RTLSDRBindings.updateGain!(obj,tul...);
function updateGain!(obj::PlutoSDR,tul...)
    # In pluto we only get a flag so we have to access to gain value to return the updated gain value
    AdalmPluto.updateGain!(obj,_toInt.(tul)...);
    return getGain(obj)
end

"""
Define Gain policy for the SDR radio. Only supported on AdalmPluto
"""
updateGainMode!(sdr::AbstractSDR) = "manual"
# No need to redefine  for Pluto backend


text:
using ChainRulesCore: canonicalize, backing, Tangent, AbstractZero
base(dx::Tangent) = backing(canonicalize(dx))
base(dx) = dx
const Zero = Union{Nothing, AbstractZero}  # Union{Zygote, Diffractor}

struct Leaf{R,S}
  rule::R
  state::S
end

function setup(rule, x; seen = Base.IdSet())
  if isnumeric(x)
    x in seen && throw(ArgumentError("Optimisers.jl does not at present handle tied weights, sorry."))
    isbits(x) || push!(seen, x)
    return Leaf(rule, init(rule, x))
  elseif isleaf(x)
    return nothing
  else
    return map(xᵢ -> setup(rule, xᵢ; seen), _trainable(x))
  end
end

subtract!(x, x̄) = iswriteable(x) ? (x .= x .- x̄) : eltype(x).(x .- x̄)

update!(::Nothing, x, ::Zero, ::Zero...) = nothing, x
update!(::Nothing, x, x̄s...) = nothing, x

update!(ℓ::Leaf, x, ::Zero, ::Zero...) = ℓ, x
function update!(ℓ::Leaf, x, x̄s...)
  s′, x̄′ = apply!(ℓ.rule, ℓ.state, x, base.(x̄s)...)
  Leaf(ℓ.rule, s′), subtract!(x, x̄′)
end

update!(tree, x, ::Zero, ::Zero...) = tree, x
function update!(tree, x, x̄s...)
  x̄s′ = map(x̄ -> functor(typeof(x), base(x̄))[1], x̄s)
  x′, re = functor(typeof(x), x)
  xtree = map((stᵢ, xᵢ, x̄sᵢ...) -> update!(stᵢ, xᵢ, x̄sᵢ...), tree, x′, x̄s′...)
  map(first, xtree), re(map(last, xtree))
end

function update(tree, x, x̄s...)
  t′ = fmap(copy, tree; exclude = iswriteable)
  x′ = fmap(copy, x; exclude = iswriteable)
  update!(t′, x′, x̄s...)
end

# default all rules to first order calls
apply!(o, state, x, dx, dxs...) = apply!(o, state, x, dx)

isnumeric(x::AbstractArray{<:Number}) = isleaf(x)  # isleaf to allow for e.g. transposed shared weights
isnumeric(x::AbstractArray{<:Integer}) = false
isnumeric(x) = false

iswriteable(::DenseArray) = true  # more elaborate versions are possible, wait until needed?
iswriteable(_) = false

"""
    trainable(x::Layer) -> NamedTuple

This should be overloaded to make optimisers ignore some fields of
every `Layer`, which would otherwise contain trainable parameters.
(Elements such as functions and sizes are always ignored.)

The default is `Functors.children(x)`, usually a NamedTuple of all fields,
and `trainable(x)` must contain a subset of these.
"""
trainable(x) = functor(x)[1]

_trainable(x) = _trainable(functor(x)[1], trainable(x))
_trainable(ch::NamedTuple, tr::NamedTuple) = merge(map(_ -> nothing, ch), tr)
_trainable(ch::Tuple{Vararg{Any,N}}, tr::Tuple{Vararg{Any,N}}) where N = tr
_trainable(ch::AbstractArray, tr::AbstractArray) = tr
function _trainable(ch::NamedTuple, tr::Tuple)  # for old Flux-style no-names tuple
  @warn "trainable(x) should now return a NamedTuple with the field names, not a Tuple" maxlog=3
  map(c -> c in tr ? c : nothing, ch)
end

"""
    @.. x = x + y

Sometimes in-place broadcasting macro, for use in `apply!` rules.
If `iswriteable(x)` then it is just `@. x = rhs`, but if not, it becomes `x = @. rhs`.
"""
macro var".."(ex)
  Meta.isexpr(ex, :(=)) || throw("the macro @.. only accepts assignment, like @.. x = y + z")
  dst = esc(ex.args[1])
  src = esc(Broadcast.__dot__(ex.args[2]))
  :($dst = if $iswriteable($dst)
      $dst .= $src
    else
      $src
    end)
end

"""
    x = @lazy y + z

Lazy broadcasting macro, for use in `apply!` rules. It broadcasts like `@.`
but does not materialise, returning a `Broadcasted` object for later use.
Beware that mutation of arguments will affect the result,
and that if it is used in two places, work will be done twice.
"""
macro lazy(ex)
  bc = esc(Broadcast.__dot__(ex))
  :($lazy.($bc))
end

function lazy end
Broadcast.broadcasted(::typeof(lazy), x) = Lazy(x)
struct Lazy{T}; bc::T; end
Broadcast.materialize(x::Lazy) = Broadcast.instantiate(x.bc)

onevalue(λ::T, x::AbstractArray{T}) where T = map(_ -> λ, x)
onevalue(λ, x::AbstractArray{T}) where T = onevalue(convert(float(T), λ), x)

function Base.show(io::IO, ℓ::Leaf)  # show method is mostly to hide its long type!
  ioc = IOContext(io, :compact => true)
  print(ioc, "Leaf(", ℓ.rule, ", ")
  show(ioc, ℓ.state)
  print(io, ")")
end


text:
"""
    timeScale(eph)

   Retrieve the timescale associated with ephemeris handler eph
   Returns 1 for TDB and 2 for TCB.

"""
function timeScale(eph::Ephem)
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    return ccall((:calceph_gettimescale , libcalceph), Cint,
                  (Ptr{Cvoid},),eph.data)
end

"""
    PositionRecord

    stores position record metadata.
"""
struct PositionRecord
   " Naif Id of target "
   target::Int
   " Naif Id of center "
   center::Int
   " start of epoch span "
   startEpoch::Float64
   "end of epoch span"
   stopEpoch::Float64
   "frame : 1 for ICRF"
   frame::Int
end

"""
    positionRecords(eph)

   Retrieve position records metadata in ephemeris associated to
   handler eph .

"""
function positionRecords(eph::Ephem)
    res = Array{PositionRecord,1}()
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    NR::Int = ccall((:calceph_getpositionrecordcount , libcalceph), Cint,
    (Ptr{Cvoid},),eph.data)
    (NR == 0) && throw(CALCEPHException("Could not find any position records!"))
    target = Ref{Cint}(0)
    center = Ref{Cint}(0)
    startEpoch = Ref{Cdouble}(0.0)
    stopEpoch = Ref{Cdouble}(0.0)
    frame = Ref{Cint}(0)
    for i=1:NR
       stat = ccall((:calceph_getpositionrecordindex , libcalceph), Cint,
                  (Ptr{Cvoid},Cint,Ref{Cint},Ref{Cint},Ref{Cdouble},Ref{Cdouble},Ref{Cint}),
                  eph.data, i ,target,center,startEpoch,stopEpoch,frame)
       if (stat!=0)
          push!(res,PositionRecord(target[],center[],startEpoch[],stopEpoch[],frame[]))
       end
    end

    return res
end

"""
    OrientationRecord

    stores orientation record metadata.
"""
struct OrientationRecord
   " Naif Id of target "
   target::Int
   " start of epoch span "
   startEpoch::Float64
   "end of epoch span"
   stopEpoch::Float64
   "frame : 1 for ICRF"
   frame::Int
end

"""
    orientationRecords(eph)

   Retrieve orientation records metadata in ephemeris associated to
   handler eph .

"""
function orientationRecords(eph::Ephem)
    res = Array{OrientationRecord,1}()
    @_checkPointer eph.data "Ephemeris is not properly initialized!"
    NR::Int = ccall((:calceph_getorientrecordcount , libcalceph), Cint,
    (Ptr{Cvoid},),eph.data)
    (NR == 0) && throw(CALCEPHException("Could not find any orientation records!"))
    target = Ref{Cint}(0)
    startEpoch = Ref{Cdouble}(0.0)
    stopEpoch = Ref{Cdouble}(0.0)
    frame = Ref{Cint}(0)
    for i=1:NR
       stat = ccall((:calceph_getorientrecordindex , libcalceph), Cint,
                  (Ptr{Cvoid},Cint,Ref{Cint},Ref{Cdouble},Ref{Cdouble},Ref{Cint}),
                  eph.data, i ,target,startEpoch,stopEpoch,frame)
       if (stat!=0)
          push!(res,OrientationRecord(target[],startEpoch[],stopEpoch[],frame[]))
       end
    end

    return res
end


text:
#========================================================================================#
#	Laboratory 4
#
# Tuples, Pairs, Dict, Symbol, files, Dates, Random numbers and downloading.
# Files and Dates are available in Julia Programming Cookbook
#
# Author: Niall Palfreyman, 13/01/2022
#========================================================================================#
[
	Activity(
		"""
		In this laboratory we introduce several helpful Julia tools that make life a
		little easier: Tuples, Pairs, Dicts, Symbols, Files, DateTimes, Random and downloading.

		A Tuple is an immutable container that can contain several different types. We
		construct Tuples using round brackets, for example:

			my_tuple = (5, 2.718, "Niall")

		What is my_tuple[2]?
		""",
		"my_tuple[2]",
		x -> x==Main.my_tuple[2]
	),
	Activity(
		"""
		You have already seen that size() returns a Tuple. Enter the following code:

			ret_size = size(zeros(3,4))

		What is the value of ret_size[1]? Try changing the value of ret_size[2].
		Finally, what is the type of ret_size?
		""",
		"typeof(ret_size)",
		x -> x==typeof(size(zeros(3,4)))
	),
	Activity(
		"""
		Tuples are especially useful when we want to define anonymous functions with
		more than one argument:

			map((x,y)->3x+2y,4,5)

		Construct a single line mapping that calculates sin(x*y) for corresponding
		elements in the two ranges x = 1:5 and y = 5:-1:1
		""",
		"map((x,y)->sin(x*y),1:5,5:-1:1)",
		x -> x==map((x,y)->sin(x*y),1:5,5:-1:1)
	),
	Activity(
		"""
		A Pair is a structure that contains two objects - typically a key and its
		entry in a dictionary. We construct a Pair like this:

			my_pair = "Yellow Submarine" => "Beatles"

		Construct my_pair and then find the value of last(my_pair):
		""",
		"",
		x -> x==last("Yellow Submarine" => "Beatles")
	),
	Activity(
		"""
		A Dict(ionary) is a hashed database of Pairs. Dicts are very lightweight, so
		we can easily use them in everyday code. We construct a Dict by passing a
		sequence of Pairs to the constructor:

			my_dict = Dict( "pi" => 3.142, "e" => 2.718)

		Construct my_dict, and notice that the Pairs are not necessarily in the same
		order that you entered them. Now look up the value of "pi" in my_dict:
		""",
		"my_dict[\"pi\"]",
		x -> x==3.142
	),
	Activity(
		"""
		We can find out whether our Dict contains the key "e" by using the keys() function:

			"e" in keys(my_dict)

		Delete the entry for "e" from my_dict using the delete!() function. What word do
		you now see in red if you enter my_dict["e"]?
		""",
		"delete!(my_dict,\"e\")",
		x -> x=="ERROR"
	),
	Activity(
		"""
		Now add two extra entries to my_dict:

			my_dict["root2"] = 1.414
			my_dict["epsilon0"] = 8.854e-12

		What is the result of calling haskey(my_dict,"root2")?
		""",
		"",
		x -> x==true
	),
	Activity(
		"""
		We have already met Symbols - they define components of the Julia language, and
		we use them to extend the language with new components. We construct Symbols
		using the colon ':'. Enter the following lines:

			a,b = 2,3
			expr = :(a+b)
			dump(expr)

		What is the value of expr.args?
		""",
		"",
		x -> x==Main.expr.args
	),
	Activity(
		"""
		You can see that the arguments of expr are Symbols waiting to be evaluated. Try:

			typeof(expr.args[2])

		What do you get if you enter string(expr)?
		""",
		"",
		x -> occursin("a+b",replace(x," "=>""))
	),
	Activity(
		"""
		Symbols and Strings are very similar to each other. You will see that graphics
		functions often use Symbols to define special switches such as :red or :filled.

		In fact, a Symbol is a String that has been prepared for being evaluated.
		What do you get if you apply the function eval() to expr?
		""",
		"eval(expr)",
		x -> x==5
	),
	Activity(
		"""
		We use the "splat" operator (...) to convert a collection into a set of
		arguments for a function. Define this function and test it with a few numbers:

			my_function(x,y,z) = x * (y+z)

		Now define this vector: v = [2,3,4]. Suppose we want to use the three numbers
		in v as arguments for my_function. First try it like this: my_function(v), and
		see what happens ...

		It didn't work, did it? Now tell me what you get when you enter this:

			my_function(v...)
		""",
		"",
		x -> x==14
	),
	Activity(
		"""
		Files are an important part of everyday data science. First of all, we need to
		be able to find files in the background filesystem. Enter: pwd() at the Julia
		prompt. This stands for Present Working Directory ...

		What you get back is the path of the current folder. Now enter: readdir() ...

		This gives you a list of filenames in the current folder. Choose a filname in
		this list - for example maybe "filename.ext". Ask whether your file is in the
		current folder. What answer do you get back?
		""",
		"\"filename.ext\" in readdir()",
		x -> x==true
	),
	Activity(
		"""
		pwd() gives us the path of the current folder, and readdir() gives us a vector
		of filenames in it. We can put these together using joinpath(). Find out the
		answer returned by the following function call:

			isfile(joinpath(pwd(),"filename.ext"))

		where filename.ext is the file you previously found in the current directory.
		""",
		"Investigate the value returned by the joinpath() call",
		x -> true
	),
	Activity(
		"""
		Now let's create a new file. The following command creates a file named
		"my_data.txt", and sets up a FILESTREAM named "file" for "w"riting to it:

			file = open("my_data.txt","w")

		Now write some information to the file:

			write( file, "This is my file\nIt belongs to me!\n")

		The value returned by write() is the number of bytes you have written to the
		file. We can also us the print() and println() functions, for example:

			println( file, "It really does!")

		Finally, close() file and tell me the result of isfile("my_data.txt"):
		""",
		"close(file)",
		x -> x==true
	),
	Activity(
		"""
		Congratulations! You have created your first file! Now let's trying
		r(eading) from the file:

			file = open("my_data.txt","r")

		Now enter: readline(file) several times to read the lines of the file.
		Once you have read all the lines, the file is in an end-of-file state:

			eof(file)

		What value does readline() return if you continue to read lines after
		the end-of-file?
		""",
		"readline(file)",
		x -> isempty(x)
	),
	Activity(
		"""
		Now close() the file, then reopen it again to start reading at the
		beginning of the file. We can read all lines of the file at once.
		What is the type of the structure returned by the following code?

			file = open("my_data.txt","r")
			readlines( file)
		""",
		"",
		x -> x <: Vector
	),
	Activity(
		"""
		If our file contained binary data, it would not be possible to read
		it in separate lines - in this case we use the read() function.

		Rewind the file to the beginning using: seekstart(file).
		Enter read(file) to see the characters in the file, and tell me the
		first character:
		""",
		"0x54: Scroll the Julia console back up to see the beginning",
		x -> x==0x54
	),
	Activity(
		"""
		Well, that wasn't very pleasant, was it, with all those characters
		screaming across the screen? Let's do it in a more civilised way this time...

		Rewind the file to the beginning using seekstart().
		Now enter:

			data = read(file);

		Did you remember to write ';' at the end of the line? If not, you
		had the "screaming characters" problem. ';' at the end of a line
		stops the return value being written to the console. Now tell me
		the value of the fifth element of data:
		""",
		"data[5]",
		x -> x==0x20
	),
	Activity(
		"""
		This hex code represents a character. Can you convert the code to
		a character?
		""",
		"Char(data[5])",
		x -> x==' '
	),
	Activity(
		"""
		We can even convert the data entirely to a String like this:

			str = String(data)

		However, this conversion uses up the data values. What value is
		now returned by the function call isempty(data)?
		""",
		"",
		x -> x==true
	),
	Activity(
		"""
		Finally, we must always close a filestream after we have finished
		with it: close(file).
		Also, we should clean up afterwards, so now remove the file we have
		created using rm("my_data.txt"), and tell me the return type of rm():
		""",
		"First call rm(), and then ask: typeof(ans)",
		x -> x==Nothing
	),
	Activity(
		"""
		Now we investigate DateTimes in Julia. Support for date and time
		handling is provided by the Dates package, which we must first load:

			using Dates

		We can access the current time using the now() function:

			datim = Dates.now()

		What is the type of datim?
		""",
		"",
		x -> x==Main.DateTime
	),
	Activity(
		"""
		To create a new date, we pass year, month and day to the constructor:

			Date( 1996, 7, 16)
			Date( 2020, 6)

		What Date value is constructed by the call Date(2022)?
		""",
		"",
		x -> x==Main.Date(2022)
	),
	Activity(
		"""
		We can also create times. Use the minute() function to find the number
		of minutes past the hour in this time:

			DateTime(1992,10,13,6,18)
		""",
		"minute(ans)",
		x -> x==Main.minute(Main.DateTime(1992,10,13,6,18))
	),
	Activity(
		"""
		The module Dates also makes available Periods of time. Use the subtypes()
		function to find the subtypes of Period and also the subtypes of these
		subtypes. How many subtypes does the type TimePeriod have?
		""",
		"subtypes(TimePeriod)",
		x -> x==length(Main.subtypes(Main.TimePeriod))
	),
	Activity(
		"""
		However, we don't just want to construct dates and times - we usually
		want to PARSE (i.e., analyse) them. We can construct a Date from a
		String by passing a DateTime format argument:

			Date("19760915","yyyymmdd")

		For DateTimes, this format gets a little more complicated, so you may
		wish to define your own format:

			format = DateFormat("HH:MM, dd.mm.yyyy")

		Use this format to parse the DateTime "06:18, 13.10.1992". What
		character separates the date from the time in the result?
		""",
		"",
		x -> x=='T'
	),
	Activity(
		"""
		DateTimes contains many useful functions that you can look up at
		docs.julialang.org. For example, use the function dayname() to find
		out the day on which you were born ...

		When you've finished experimenting, just enter reply() to move on
		to the next activity.
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Find out your age in days by subtracting your birthday Date from
		today(), then move on with reply():
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Form a list of the past 8 days by collecting this Range into a Vector:

			today()-Week(1):Day(1):today()
		""",
		"",
		x -> x==Main.eval(:(collect((today()-Week(1)):Day(1):today())))
	),
	Activity(
		"""
		Next we'll look at a very important tool of data science: Random numbers.
		First load the functions that we'll be using:

			using Random: rand, randn, seed!

		By itself, the rand() function returns a pseudo-random Float64 number in the
		half-open interval [0.0,1.0). Try this now.

		If you pass rand() a range as its first argument, it returns a random element
		from that range. If you pass it a tuple of integers, it returns an array with
		the size specified by those integers. Use rand() to produce a (2,3) array of
		integer values between -1 and +1:
		""",
		"rand(Range,Size)",
		x -> (typeof(x) <: Matrix{Int}) && size(x) == (2,3) && all(map(y->in(y,-1:1),x))
	),
	Activity(
		"""
		randn() works just the same as rand(), but returns normally distributed values
		with mean 0.0 and standard deviation 1.0. Create a (5,7) array of such normally
		distributed numbers:
		""",
		"",
		x -> size(x) == (5,7) && sum(x)/length(x) < 0.5
	),
	Activity(
		"""
		When we perform simulations with random numbers, we never know in advance how
		our program will behave. On the one hand, this is certainly good, because many
		biological processes are essentially random. On the other hand, it can be very
		frustrating if you observe some special problem or phenomenon, because you may
		never again be able to reproduce that special situation. Because of this, we
		need to be able to make random-number generation REPRODUCIBLE. We do this by
		SEEDING the random-number generator (rng).

		To do this, we use the function seed!() at the beginning of our program to make
		sure all random numbers follow an identical pattern across separate runs of the
		program. Run the following code several times, then tell me what result you get:

			seed!(123); rand(5)
		""",
		"",
		x -> x==(Main.seed!(123); Main.rand(5))
	),
	Activity(
		"""
		OK, and fi-inally at the end of this very long (phew!) laboratory, we look
		briefly at how to download resources from the internet. First, we load the
		download() function from the package Downloads:

			using Downloads: download

		Next we define the url of our resource:

			url = "https://raw.githubusercontent.com/NiallPalfreyman/Ingolstadt.jl/main/src/Ingolstadt.jl"

		Next, we download this page into a local file:

			file = download(url)

		Use readlines() (don't forget the ';'!) to discover the Date on which Niall
		Palfreyman started writing the Ingolstadt project:
		""",
		"data = readlines(file);",
		x -> x == Main.Date("7/12/2021","d/mm/yyyy")
	),
	Activity(
		"""
		OK, that's the end of this laboratory. The resource we have downloaded is my
		source code - feel free to explore it and use it as much as you like. Bye! :-)
		""",
		"",
		x -> true
	),
]


text:
function learnable_normal(x::Float64, y::Float64)
    l = learnable(:l)
    m = learnable(:m)
    q = rand(:q, Normal(l, 3 + m^2))
    return q
end

@testset "Convergence for learning - MAP 1" begin
    sel = selection([(:q, ) => 6.0])
    params = parameters([(:l, ) => 3.0,
                         (:m, ) => 10.0])
    for i in 1:200
        ret, cl, w = generate(sel, params, learnable_normal, 5.0, 3.0)
        param_grads = get_parameter_gradients(params, cl, 1.0)
        params = update_parameters(ADAM(0.05, (0.9, 0.8)), params, param_grads)
    end
    @test params[:l] ≈ 6.0 atol = 1e-2
    @test params[:m] ≈ 0.0 atol = 1e-2

    # Train.
    params = train(sel, params, learnable_normal, 5.0, 3.0; opt = ADAM(0.05, (0.9, 0.8)), iters = 2000)
    @test params[:l] ≈ 6.0 atol = 1e-2
    @test params[:m] ≈ 0.0 atol = 1e-2
end

function model()
    slope = rand(:slope, Normal(-1, exp(0.5)))
    intercept = rand(:intercept, Normal(1, exp(2.0)))
end

function var()
    slope_mu = learnable(:slope_mu)
    slope_log_std = learnable(:slope_log_std)
    intercept_mu = learnable(:intercept_mu)
    intercept_log_std = learnable(:intercept_log_std)
    slope = rand(:slope, Normal(slope_mu, exp(slope_log_std)))
    intercept = rand(:intercept, Normal(intercept_mu, exp(intercept_log_std)))
end

@testset "Convergence for learning - VI 1" begin
    sel = selection()
    initial_params = parameters([(:slope_mu, ) => 0.0,
                                 (:slope_log_std, ) => 0.0,
                                 (:intercept_mu, ) => 0.0,
                                 (:intercept_log_std, ) => 0.0])
    params, _, _ = advi(sel, initial_params, var, (), model, (); iters = 2000, opt = ADAM(0.005, (0.9, 0.999)))
    @test params[:slope_mu] ≈ -1.0 atol = 7e-2
    @test params[:slope_log_std] ≈ 0.5 atol = 7e-2
    @test params[:intercept_mu] ≈ 1.0 atol = 7e-2
    @test params[:intercept_log_std] ≈ 2.0 atol = 7e-2
end


text:
@test first(evaluate5(A♣, A♡, A♢, A♠, K♠)) == 11 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, Q♠)) == 12 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, J♠)) == 13 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, T♠)) == 14 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 9♠)) == 15 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 8♠)) == 16 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 7♠)) == 17 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 6♠)) == 18 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 5♠)) == 19 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 4♠)) == 20 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 3♠)) == 21 # Four Aces
@test first(evaluate5(A♣, A♡, A♢, A♠, 2♠)) == 22 # Four Aces
@test first(evaluate5(K♣, K♡, K♢, K♠, A♠)) == 23 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, Q♠)) == 24 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, J♠)) == 25 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, T♠)) == 26 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 9♠)) == 27 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 8♠)) == 28 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 7♠)) == 29 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 6♠)) == 30 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 5♠)) == 31 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 4♠)) == 32 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 3♠)) == 33 # Four Kings
@test first(evaluate5(K♣, K♡, K♢, K♠, 2♠)) == 34 # Four Kings
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, A♠)) == 35 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, K♠)) == 36 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, J♠)) == 37 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, T♠)) == 38 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 9♠)) == 39 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 8♠)) == 40 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 7♠)) == 41 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 6♠)) == 42 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 5♠)) == 43 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 4♠)) == 44 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 3♠)) == 45 # Four Queens
@test first(evaluate5(Q♣, Q♡, Q♢, Q♠, 2♠)) == 46 # Four Queens
@test first(evaluate5(J♣, J♡, J♢, J♠, A♠)) == 47 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, K♠)) == 48 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, Q♠)) == 49 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, T♠)) == 50 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 9♠)) == 51 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 8♠)) == 52 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 7♠)) == 53 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 6♠)) == 54 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 5♠)) == 55 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 4♠)) == 56 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 3♠)) == 57 # Four Jacks
@test first(evaluate5(J♣, J♡, J♢, J♠, 2♠)) == 58 # Four Jacks
@test first(evaluate5(T♣, T♡, T♢, T♠, A♠)) == 59 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, K♠)) == 60 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, Q♠)) == 61 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, J♠)) == 62 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 9♠)) == 63 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 8♠)) == 64 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 7♠)) == 65 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 6♠)) == 66 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 5♠)) == 67 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 4♠)) == 68 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 3♠)) == 69 # Four Tens
@test first(evaluate5(T♣, T♡, T♢, T♠, 2♠)) == 70 # Four Tens
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, A♠)) == 71 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, K♠)) == 72 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, Q♠)) == 73 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, J♠)) == 74 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, T♠)) == 75 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 8♠)) == 76 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 7♠)) == 77 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 6♠)) == 78 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 5♠)) == 79 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 4♠)) == 80 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 3♠)) == 81 # Four Nines
@test first(evaluate5(9♣, 9♡, 9♢, 9♠, 2♠)) == 82 # Four Nines
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, A♠)) == 83 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, K♠)) == 84 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, Q♠)) == 85 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, J♠)) == 86 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, T♠)) == 87 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 9♠)) == 88 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 7♠)) == 89 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 6♠)) == 90 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 5♠)) == 91 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 4♠)) == 92 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 3♠)) == 93 # Four Eights
@test first(evaluate5(8♣, 8♡, 8♢, 8♠, 2♠)) == 94 # Four Eights
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, A♠)) == 95 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, K♠)) == 96 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, Q♠)) == 97 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, J♠)) == 98 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, T♠)) == 99 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 9♠)) == 100 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 8♠)) == 101 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 6♠)) == 102 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 5♠)) == 103 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 4♠)) == 104 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 3♠)) == 105 # Four Sevens
@test first(evaluate5(7♣, 7♡, 7♢, 7♠, 2♠)) == 106 # Four Sevens
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, A♠)) == 107 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, K♠)) == 108 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, Q♠)) == 109 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, J♠)) == 110 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, T♠)) == 111 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 9♠)) == 112 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 8♠)) == 113 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 7♠)) == 114 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 5♠)) == 115 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 4♠)) == 116 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 3♠)) == 117 # Four Sixes
@test first(evaluate5(6♣, 6♡, 6♢, 6♠, 2♠)) == 118 # Four Sixes
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, A♠)) == 119 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, K♠)) == 120 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, Q♠)) == 121 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, J♠)) == 122 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, T♠)) == 123 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 9♠)) == 124 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 8♠)) == 125 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 7♠)) == 126 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 6♠)) == 127 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 4♠)) == 128 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 3♠)) == 129 # Four Fives
@test first(evaluate5(5♣, 5♡, 5♢, 5♠, 2♠)) == 130 # Four Fives
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, A♠)) == 131 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, K♠)) == 132 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, Q♠)) == 133 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, J♠)) == 134 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, T♠)) == 135 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 9♠)) == 136 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 8♠)) == 137 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 7♠)) == 138 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 6♠)) == 139 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 5♠)) == 140 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 3♠)) == 141 # Four Fours
@test first(evaluate5(4♣, 4♡, 4♢, 4♠, 2♠)) == 142 # Four Fours
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, A♠)) == 143 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, K♠)) == 144 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, Q♠)) == 145 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, J♠)) == 146 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, T♠)) == 147 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 9♠)) == 148 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 8♠)) == 149 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 7♠)) == 150 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 6♠)) == 151 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 5♠)) == 152 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 4♠)) == 153 # Four Treys
@test first(evaluate5(3♣, 3♡, 3♢, 3♠, 2♠)) == 154 # Four Treys
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, A♠)) == 155 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, K♠)) == 156 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, Q♠)) == 157 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, J♠)) == 158 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, T♠)) == 159 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 9♠)) == 160 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 8♠)) == 161 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 7♠)) == 162 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 6♠)) == 163 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 5♠)) == 164 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 4♠)) == 165 # Four Deuces
@test first(evaluate5(2♣, 2♡, 2♢, 2♠, 3♠)) == 166 # Four Deuces


text:
const EXPONENTS = [
    1e0,   1e1,   1e2,   1e3,   1e4,   1e5,   1e6,   1e7,   1e8,    1e9,
    1e10,  1e11,  1e12,  1e13,  1e14,  1e15,  1e16,  1e17,  1e18,  1e19,
    1e20,  1e21,  1e22,  1e23,  1e24,  1e25,  1e26,  1e27,  1e28,  1e29,
    1e30,  1e31,  1e32,  1e33,  1e34,  1e35,  1e36,  1e37,  1e38,  1e39,
    1e40,  1e41,  1e42,  1e43,  1e44,  1e45,  1e46,  1e47,  1e48,  1e49,
    1e50,  1e51,  1e52,  1e53,  1e54,  1e55,  1e56,  1e57,  1e58,  1e59,
    1e60,  1e61,  1e62,  1e63,  1e64,  1e65,  1e66,  1e67,  1e68,  1e69,
    1e70,  1e71,  1e72,  1e73,  1e74,  1e75,  1e76,  1e77,  1e78,  1e79,
    1e80,  1e81,  1e82,  1e83,  1e84,  1e85,  1e86,  1e87,  1e88,  1e89,
    1e90,  1e91,  1e92,  1e93,  1e94,  1e95,  1e96,  1e97,  1e98,  1e99,
    1e100, 1e101, 1e102, 1e103, 1e104, 1e105, 1e106, 1e107, 1e108, 1e109,
    1e110, 1e111, 1e112, 1e113, 1e114, 1e115, 1e116, 1e117, 1e118, 1e119,
    1e120, 1e121, 1e122, 1e123, 1e124, 1e125, 1e126, 1e127, 1e128, 1e129,
    1e130, 1e131, 1e132, 1e133, 1e134, 1e135, 1e136, 1e137, 1e138, 1e139,
    1e140, 1e141, 1e142, 1e143, 1e144, 1e145, 1e146, 1e147, 1e148, 1e149,
    1e150, 1e151, 1e152, 1e153, 1e154, 1e155, 1e156, 1e157, 1e158, 1e159,
    1e160, 1e161, 1e162, 1e163, 1e164, 1e165, 1e166, 1e167, 1e168, 1e169,
    1e170, 1e171, 1e172, 1e173, 1e174, 1e175, 1e176, 1e177, 1e178, 1e179,
    1e180, 1e181, 1e182, 1e183, 1e184, 1e185, 1e186, 1e187, 1e188, 1e189,
    1e190, 1e191, 1e192, 1e193, 1e194, 1e195, 1e196, 1e197, 1e198, 1e199,
    1e200, 1e201, 1e202, 1e203, 1e204, 1e205, 1e206, 1e207, 1e208, 1e209,
    1e210, 1e211, 1e212, 1e213, 1e214, 1e215, 1e216, 1e217, 1e218, 1e219,
    1e220, 1e221, 1e222, 1e223, 1e224, 1e225, 1e226, 1e227, 1e228, 1e229,
    1e230, 1e231, 1e232, 1e233, 1e234, 1e235, 1e236, 1e237, 1e238, 1e239,
    1e240, 1e241, 1e242, 1e243, 1e244, 1e245, 1e246, 1e247, 1e248, 1e249,
    1e250, 1e251, 1e252, 1e253, 1e254, 1e255, 1e256, 1e257, 1e258, 1e259,
    1e260, 1e261, 1e262, 1e263, 1e264, 1e265, 1e266, 1e267, 1e268, 1e269,
    1e270, 1e271, 1e272, 1e273, 1e274, 1e275, 1e276, 1e277, 1e278, 1e279,
    1e280, 1e281, 1e282, 1e283, 1e284, 1e285, 1e286, 1e287, 1e288, 1e289,
    1e290, 1e291, 1e292, 1e293, 1e294, 1e295, 1e296, 1e297, 1e298, 1e299,
    1e300, 1e301, 1e302, 1e303, 1e304, 1e305, 1e306, 1e307, 1e308,
]

pow10(exp) = (@inbounds v = EXPONENTS[exp+1]; return v)

maxexponent(::Type{Int16}) = 4
maxexponent(::Type{Int32}) = 38
maxexponent(::Type{Int64}) = 308

minexponent(::Type{Int16}) = -5
minexponent(::Type{Int32}) = -38
minexponent(::Type{Int64}) = -308

inttype(::Type{Float16}) = Int16
inttype(::Type{Float32}) = Int32
inttype(::Type{Float64}) = Int64

const BIGN = UInt8('N')
const LITTLEN = UInt8('n')
const BIGA = UInt8('A')
const LITTLEA = UInt8('a')
const BIGI = UInt8('I')
const LITTLEI = UInt8('i')
const BIGF = UInt8('F')
const LITTLEF = UInt8('f')
const BIGT = UInt8('T')
const LITTLET = UInt8('t')
const BIGY = UInt8('Y')
const LITTLEY = UInt8('y')
const BIGE = UInt8('E')
const LITTLEE = UInt8('e')

ParsingException(::Type{<:AbstractFloat}, exp::Signed, row, col) = ParsingException("error parsing a `$T` value on column $col, row $row; exponent out of range: $exp")

function scale(exp, v::T, frac, row, col) where T
    if exp >= 0
        max_exp = maxexponent(T)
        exp > max_exp && throw(ParsingException(T, exp, row, col))
        if exp > 15
            return Float64(Base.TwicePrecision{Float64}(v) * Base.TwicePrecision{Float64}(pow10(exp)))
        else
            return v * pow10(exp)
        end
    else
        min_exp = minexponent(T)
        if exp < min_exp
            -exp + min_exp > -min_exp && throw(ParsingException(T, exp, row, col))
            return Float64(Base.TwicePrecision{Float64}(v) / Base.TwicePrecision{Float64}(pow10(-exp + min_exp)))
        else
            if exp > 15
                return Float64(Base.TwicePrecision{Float64}(v) / Base.TwicePrecision{Float64}(pow10(-exp)))
            else
                return v / pow10(-exp)
            end
        end
    end
end

function parsefield(io::IO, ::Type{T}, opt::CSV.Options, row, col, state, ifmissing::Function) where {T <: Union{Float16, Float32, Float64}}
    mark(io)
    @checkmissingstart()
    minussign = plussign = false
    if b == MINUS # check for leading '-' or '+'
        minussign = true
        c = peekbyte(io)
        if (NEG_ONE < c < TEN) || c == opt.decimal
            b = readbyte(io)
        end
    elseif b == PLUS
        plussign = true
        c = peekbyte(io)
        if (NEG_ONE < c < TEN) || c == opt.decimal
            b = readbyte(io)
        end
    end
    # float digit parsing
    iT = inttype(T)
    v = zero(iT)
    parseddigits = false
    while NEG_ONE < b < TEN
        parseddigits = true
        # process digits
        v *= iT(10)
        v += iT(b - ZERO)
        eof(io) && (state[] = EOF; result = T(v); @goto done)
        b = readbyte(io)
    end
    # if we didn't get any digits and character isn't leading dot, check for NaN/Inf
    if !parseddigits && b != opt.decimal
        if minussign || plussign # skip sign character, if any
            eof(io) && @goto checkmissingend
            b = readbyte(io)
        end
        if b == LITTLEN || b == BIGN
            eof(io) && @goto checkmissingend
            b = readbyte(io)
            (!(b == LITTLEA || b == BIGA) || eof(io)) && (reset(io); b = readbyte(io); @goto checkmissingend)
            b = readbyte(io)
            !(b == LITTLEN || b == BIGN) && (reset(io); b = readbyte(io); @goto checkmissingend)
            result = T(NaN)
            eof(io) && (state[] = EOF; @goto done)
            b = readbyte(io)
            @goto checkdone
        elseif b == LITTLEI || b == BIGI
            eof(io) && @goto checkmissingend
            b = readbyte(io)
            (!(b == LITTLEN || b == BIGN) || eof(io)) && (reset(io); b = readbyte(io); @goto checkmissingend)
            b = readbyte(io)
            !(b == LITTLEF || b == BIGF) && (reset(io); b = readbyte(io); @goto checkmissingend)
            result = T(Inf)
            eof(io) && (state[] = EOF; @goto done)
            b = readbyte(io)
            if b == LITTLEI || b == BIGI
                # read the rest of INFINITY
                eof(io) && (state[] = EOF; @goto done)
                b = readbyte(io)
                b == LITTLEN || b == BIGN || @goto checkdone
                eof(io) && (state[] = EOF; @goto done)
                b = readbyte(io)
                b == LITTLEI || b == BIGI || @goto checkdone
                eof(io) && (state[] = EOF; @goto done)
                b = readbyte(io)
                b == LITTLET || b == BIGT || @goto checkdone
                eof(io) && (state[] = EOF; @goto done)
                b = readbyte(io)
                b == LITTLEY || b == BIGY || @goto checkdone
                eof(io) && (state[] = EOF; @goto done)
                b = readbyte(io)
            end
            @goto checkdone
        else
            @goto checkmissingend
        end
    end
    # parse fractional part
    frac = 0
    result = T(v)
    if b == opt.decimal
        eof(io) && (state[] = EOF; parseddigits ? @goto(done) : @goto(error))
        b = readbyte(io)
    elseif b == LITTLEE || b == BIGE
        @goto parseexp
    else
        @goto checkdone
    end

    while NEG_ONE < b < TEN
        frac += 1
        # process digits
        v *= iT(10)
        v += iT(b - ZERO)
        eof(io) && (state[] = EOF; result = scale(-frac, v, 0, row, col); @goto done)
        b = readbyte(io)
    end
    # parse potential exp
    if b == LITTLEE || b == BIGE
        @label parseexp
        eof(io) && (state[] = EOF; result = scale(-frac, v, 0, row, col); @goto done)
        b = readbyte(io)
        exp = zero(iT)
        negativeexp = false
        if b == MINUS
            negativeexp = true
            b = readbyte(io)
        elseif b == PLUS
            b = readbyte(io)
        end
        parseddigits = false
        while NEG_ONE < b < TEN
            parseddigits = true
            # process digits
            exp *= iT(10)
            exp += iT(b - ZERO)
            eof(io) && (state[] = EOF; result = scale(ifelse(negativeexp, -exp, exp) - frac, v, frac, row, col); @goto done)
            b = readbyte(io)
        end
        result = parseddigits ? scale(ifelse(negativeexp, -exp, exp) - frac, v, frac, row, col) : scale(-frac, v, 0, row, col)
    else
        result = scale(-frac, v, 0, row, col)
    end

    @label checkdone
    @checkdone(done)
    @goto checkmissingend

    @label checkmissingend
    @checkmissingend()
    @goto error

    @label done
    return T(ifelse(minussign, -result, result))

    @label missing
    return ifmissing(row, col)

    @label error
    throw(ParsingException(T, b, row, col))
end

# (-73.99378204345703, -73.99378204345702)
#  (-73.95227813720703, -73.95227813720702)
#  (-73.98616027832031, -73.98616027832033)
#  (-74.00163269042969, -74.00163269042967)
#  (-73.96940612792969, -73.9694061279297)
#  (-73.96797943115234, -73.96797943115236)
#  (-73.95426940917969, -73.95426940917967)
#  (-73.97286224365234, -73.97286224365236)
#  (-73.99149322509766, -73.99149322509767)
#  (-73.97639465332031, -73.97639465332033)
#  (-73.97297668457031, -73.97297668457033)
#  (-73.98991394042969, -73.98991394042967)
#  (-73.98424530029297, -73.98424530029298)
#  (-73.97872161865234, -73.97872161865236)
#  (-73.99348449707031, -73.99348449707033)
#  (-73.96598815917969, -73.96598815917967)
#  (-73.98743438720703, -73.98743438720702)