text:
macro temp_grad_min(c1, x1, y1, c2, x2, y2)
    sym_grad11 = Symbol(String(x1)*"_grad")
    sym_grad12 = Symbol(String(y1)*"_grad")
    sym_grad21 = Symbol(String(x2)*"_grad")
    sym_grad22 = Symbol(String(y2)*"_grad")
    esc(quote
        grad_temp = $c1*$x1 < $c1*$y1 ? $c1*$sym_grad11 : $c1*$sym_grad12
        grad_temp += $c2*$x2 < $c2*$y2 ? $c2*$sym_grad21 : $c2*$sym_grad22
        end)
end

macro temp_grad_max(c1, x1, y1, c2, x2, y2)
    sym_grad11 = Symbol(String(x1)*"_grad")
    sym_grad12 = Symbol(String(y1)*"_grad")
    sym_grad21 = Symbol(String(x2)*"_grad")
    sym_grad22 = Symbol(String(y2)*"_grad")
    esc(quote
        grad_temp = $c1*$x1 > $c1*$y1 ? $c1*$sym_grad11 : $c1*$sym_grad12
        grad_temp += $c2*$x2 > $c2*$y2 ? $c2*$sym_grad21 : $c2*$sym_grad22
        end)
end

min9(args...) = findmin(args)
max9(args...) = findmax(args)

@inline function mult_apriori_kernel(x1::MC{N,T}, x2::MC{N,T}, z::Interval{Float64},
                                    u1::Float64, u2::Float64, a1::Float64, a2::Float64,
                                    v1::Float64, v2::Float64, b1::Float64, b2::Float64,
                                    u1grad::SVector{N,Float64}, u2grad::SVector{N,Float64},
                                    v1grad::SVector{N,Float64}, v2grad::SVector{N,Float64}) where {N, T<:Union{NS, MV}}

    x1L = lo(x1);  x1U = hi(x1);  x2L = lo(x2);  x2U = hi(x2)
    x1cv = x1.cv;  x1cc = x1.cc;  x2cv = x2.cv; x2cc = x2.cc
    x1cv_grad = x1.cv_grad;  x1cc_grad = x1.cc_grad;
    x2cv_grad = x2.cv_grad;  x2cc_grad = x2.cc_grad

    w = x1*x2
    w0cv = w.cv
    w0cc = w.cc

    # additional constraints from underestimators (convex)
    w1cv = (x2U - a2)*u1 + (x1U - a1)*u2 + min(a2*x1cv, a2*x1cc) + min(a1*x2cv, a1*x2cc) + a1*a2 - a1*x2U - a2*x1U # GOOD
    w2cv = (x2U - x2L)*u1 + min(x2L*x1cv, x2L*x1cc) + min(a1*x2cv, a1*x2cc) - a1*x2U #BAD
    w3cv = (x1U - x1L)*u2 + min(a2*x1cv, a2*x1cc) + min(x1L*x2cv, x1L*x2cc) - x1U*a2
    w4cv = (a2 - x2L)*u1 + (a1 - x1L)*u2 + min(x2L*x1cv, x2L*x1cc) + min(x1L*x2cv, x1L*x2cc) - a1*a2

    # additional constraints from underestimators (concave)
    w1cc = (x2L - a2)*u1 + (a1 - x1U)*u2 + max(a2*x1cv, a2*x1cc) + max(x1U*x2cv, x1U*x2cc) - a1*x2L
    w2cc = (x2L - x2U)*u1 + max(a1*x2cv, a1*x2cc) + max(x2U*x1cv, x2U*x1cc) - x2L*a1
    w3cc = (x1L - x1U)*u2 + max(a2*x1cv, a2*x1cc) + max(x1U*x2cv, x1U*x2cc) - x1L*a2
    w4cc = (a2 - x2U)*u1 + (x1L - a1)*u2 + max(x2U*x1cv, x2U*x1cc) + max(a1*x2cv, a1*x2cc) - x1L*a2

    # additional constraints from overestimators (convex)
    # 5 fails,
    w5cv = (x2L - b2)*v1 + (x1L - b1)*v2 + min(b2*x1cv, b2*x1cc) + min(x1L*x2cv, x1L*x2cc) + b1*b2 - b1*x2L - x1L*b2
    w6cv = (x2L - x2U)*v1 + min(x2U*x1cv, x2U*x1cc) + min(b1*x2cv, b1*x2cc) - b1*x2L
    w7cv = (x1L - x1U)*v2 + min(b2*x1cv, b2*x1cc) + min(x1U*x2cv, x1U*x2cc) - b2*x1L
    w8cv = (b2 - x2U)*v1 + (b1 - x1U)*v2 + min(x2U*x1cv, x2U*x1cc) + min(x1U*x2cv, x1U*x2cc) - b1*b2

    # additional constraints from overestimators (concave)
    w5cc = (x2U - b2)*v1 + (b1 - x1L)*v2 + max(b2*x1cv, b2*x1cc) + max(x1L*x2cv, x1L*x2cc) - b1*x2U
    w6cc = (x2U - x2L)*v1 + max(b1*x2cv, b1*x2cc) + max(x2L*x1cv, x2L*x1cc) - x1U*b1
    w7cc = (x1U - x1L)*v2 + max(b2*x1cv, b2*x1cc) + max(x1L*x2cv, x1L*x2cc) - x1U*b2
    w8cc = (b2 - x2L)*v1 + (x1U - b1)*v2 + max(x2L*x1cv, x2L*x1cc) + max(b1*x2cv, b1*x2cc) - x1U*b2

    cv, cvind = max9(w0cv, w1cv, w2cv, w3cv, w4cv, w6cv, w7cv, w8cv)
    #cv, cvind = max9(w0cv, w1cv, w2cv, w3cv, w4cv, w5cv, w6cv, w7cv, w8cv)

    #cv = w8cv
    #cvind = 9

    if cvind == 1
        cv_grad = w.cv_grad
    elseif cvind == 2
        @temp_grad_min(a2, x1cv, x1cc, a1, x2cv, x2cc)
        cv_grad = (x2U - a2)*u2grad + (x1U - a1)*u2grad + grad_temp
    elseif cvind == 3
        @temp_grad_min(x2L, x1cv, x1cc, a1, x2cv, x2cc)
        cv_grad = (x2U - x2L)*u2grad + grad_temp
    elseif cvind == 4
        @temp_grad_min(a2, x1cv, x1cc, x1L, x2cv, x2cc)
        cv_grad = (x1U - x1L)*u2grad + grad_temp
    elseif cvind == 5
        @temp_grad_min(x2L, x1cv, x1cc, x1L, x2cv, x2cc)
        cv_grad = (a2 - x2L)*u1grad + (a1 - x1L)*u2grad + grad_temp
    elseif cvind == 6
        @temp_grad_min(b2, x1cv, x1cc, x1L, x2cv, x2cc)
        cv_grad = (x2L - b2)*v1grad + (x1L - b1)*v2grad + grad_temp
    elseif cvind == 7
        @temp_grad_min(x2U, x1cv, x1cc, b1, x2cv, x2cc)
        cv_grad = (x2L - x2U)*v1grad + grad_temp
    elseif cvind == 8
        @temp_grad_min(b2, x1cv, x1cc, b1, x2cv, x2cc)
        cv_grad = (x1L - x1U)*v2grad + grad_temp
    elseif cvind == 9
        @temp_grad_min(x2U, x1cv, x1cc, x1U, x2cv, x2cc)
        cv_grad = (b2 - x2U)*v1grad + (b1 - x1U)*v2grad + grad_temp
    end

    # w5cc
    cc, ccind = min9(w0cc, w1cc, w2cc, w3cc, w4cc, w6cc, w7cc, w8cc)

    if ccind == 1
        cc_grad = w.cc_grad
    elseif ccind == 2
        @temp_grad_max(a2, x1cv, x1cc, x1U, x2cv, x2cc)
        cc_grad = (x2L - a2)*u1grad + (a1 - x1U)*u2grad + grad_temp
    elseif ccind == 3
        @temp_grad_max(a1, x1cv, x1cc, x2U, x1cv, x1cc)
        cc_grad = (x2L - x2U)*u1grad + grad_temp
    elseif ccind == 4
        @temp_grad_max(a2, x1cv, x1cc, x1U, x2cv, x2cc)
        cc_grad = (x1L - x1U)*u2grad + grad_temp
    elseif ccind == 5
        @temp_grad_max(x2U, x1cv, x1cc, a1, x2cv, x2cc)
        cc_grad = (a2 - x2U)*u1grad + (x1L - a1)*u2grad + grad_temp
    elseif ccind == 6
        @temp_grad_max(b2, x1cv, x1cc, x1L, x2cv, x2cc)
        cc_grad = (x2U - b2)*v1grad + (b1 - x1L)*v2grad + grad_temp
    elseif ccind == 7
        @temp_grad_max(b2, x2cv, x2cc, x2L, x1cv, x1cc)
        cc_grad = (x2U - x2L)*v1grad + grad_temp
    elseif ccind == 8
        @temp_grad_max(b2, x1cv, x1cc, x1L, x2cv, x2cc)
        cc_grad = (x1U - x1L)*v2grad + grad_temp
    elseif ccind == 9
        @temp_grad_max(x2L, x1cv, x1cc, b1, x2cv, x2cc)
        cc_grad = (b2 - x2L)*v1grad + (x1U - b1)*v2grad + grad_temp
    end

    cv, cc, cv_grad, cc_grad = cut(z.lo, z.hi, cv, cc, cv_grad, cc_grad)
    MC{N,T}(cv, cc, z, cv_grad, cc_grad, x1.cnst && x2.cnst)
end


text:
function __init__()
  @require Distributions="31c24e10-a181-5473-b8eb-7969acd0382f" include("distributions.jl")
  @require StatsFuns="4c63d2b9-4356-54db-8cca-17b64c39e42c" include("statsfuns.jl")
end

ext: jl


text:
using AdventOfCodeSolutions
using Combinatorics
using Test

function input(puzzle::Puzzle{2020, 9, n}) where n
    io = openInput(puzzle)
    split(read(io, String), "\n", keepempty=false)
end

function parseInput(input)
    parse.(Int, input)
end

function findInvalidXmas(numbers, preambleSize, idx)
    n = numbers[idx]
    preamble = @view numbers[idx - preambleSize:idx - 1]

    if (n in map(sum, combinations(preamble, 2)))
        return findInvalidXmas(numbers, preambleSize, idx + 1)
    end

    return n
end

function solve(::Puzzle{2020, 9, 1}, input, preambleSize)
    numbers = map(parseInput, input)
    findInvalidXmas(numbers, preambleSize, preambleSize + 1)
end

function findWeakness(numbers, target, left, right)
    r = @view numbers[left:right]
    s = sum(r)
    if (s > target)
        return findWeakness(numbers, target, left + 1, right)
    elseif (s < target)
        return findWeakness(numbers, target, left, right + 1)
    else
        return sum(extrema(r))
    end
end

function solve(::Puzzle{2020, 9, 2}, input, target)
    numbers = map(parseInput, input)
    return findWeakness(numbers, target, 1, 2)
end

@testset "Day 9" begin
    TEST = split("""
35
20
15
25
47
40
62
55
65
95
102
117
150
182
127
219
299
277
309
576
""", '\n', keepempty=false)

    @test(solve(Puzzle(9, 1), TEST, 5) == 127)
    @test(solve(Puzzle(9, 1), input(Puzzle(9, 1)), 25) == 10884537)

    @test(solve(Puzzle(9, 2), TEST, 127) == 62)
    @test(solve(Puzzle(9, 2), input(Puzzle(9, 2)), 10884537) == 1261309)
end


text:
#
# Plot the results of a SIRD simulation
#

using Plots

function plot( input :: SIRDInput, traj :: SIRDTraj ;
              size=(400,400) )

  time = traj.time
  U = traj.U
  S = traj.S
  D = traj.D
  I = traj.I
  nsteps = traj.nsteps

  xmax = time[nsteps]

  plot(size=(400,400))
  plot!(time,U,label="Susceptible: $(@sprintf("%2.1f",U[n]*100))%",linewidth=2,color=:blue)
  plot!(time,S,label="Sick: $(@sprintf("%2.1f",D[n]*100))%",linewidth=2,color=:red)
  plot!(time,I,label="Immune: $(@sprintf("%2.1f",I[n]*100))%",linewidth=2,color=:darkgreen)
  plot!(time,D,label="Dead: $(@sprintf("%2.1f",M[n]*100))%",linewidth=2,color=:black)
  plot!(xlabel="time")
  plot!(ylabel="Fraction of the population")
  plot!(xlim=[0,xmax,ylim=[0,1])

  plot!(rectangle((7/20)*xmax, # xlength
                  0.210*traj.n, #ylength
                  0, #xmin
                  0.815*traj.n, #ymin
                  ), opacity=0.9,label="",color=:white)
  plot!(rectangle((9/20)*xmax, # xlength
                  0.110*traj.n, #ylength
                  (11/20)*xmax, #xmin
                  0.910*traj.n, #ymin
                  ), opacity=0.9,label="",color=:white)

  fontsize=8
  annotate!(0.05*xmax,input.n-0.00*input.n,text("Susceptible: $(U[i])",:left,fontsize,:serif,:blue))
  annotate!(0.05*xmax,input.n-0.05*input.n,text("Sick: $(S[i])",:left,fontsize,:serif,:red))
  annotate!(0.05*xmax,input.n-0.10*input.n,text("Dead: $(D[i])",:left,fontsize,:serif,:black))
  annotate!(0.05*xmax,input.n-0.15*input.n,text("Immune: $(I[i])",:left,fontsize,:serif,:darkgreen))
  annotate!(0.98*xmax,input.n-0.00*input.n,text("kc = $kc ; kd = $kd",:right,fontsize,:serif,:black))
  annotate!(0.98*xmax,input.n-0.05*input.n,text("ki = $ki ; kis = $kis",:right,fontsize,:serif,:black))

    next!(p)

  end

  gif(anim, filename, fps = fps)

end


text:
function allocate_matrix end
function allocate_matrix_and_vector end

"""
    allocate_vector(::Type{V},indices) where V

Allocate a vector of type `V` indexable at the indices `indices`
"""
function allocate_vector(::Type{V},indices) where V
  n = length(indices)
  allocate_vector(V,n)
end

function allocate_vector(::Type{V},n::Integer) where V
  T = eltype(V)
  zeros(T,n)
end

"""
    allocate_in_range(::Type{V},matrix) where V

Allocate a vector of type `V` in the range of matrix `matrix`.
"""
function allocate_in_range(::Type{V},matrix) where V
  n = size(matrix,1)
  allocate_vector(V,n)
end

"""
    allocate_in_domain(::Type{V},matrix) where V

Allocate a vector of type `V` in the domain of matrix `matrix`.
"""
function allocate_in_domain(::Type{V},matrix) where V
  n = size(matrix,2)
  allocate_vector(V,n)
end

"""
    fill_entries!(a,v)

Fill the entries of array `a` with the value `v`. Returns `a`.
"""
function fill_entries!(a,v)
  fill!(a,v)
  a
end

"""
    copy_entries!(a,b)

Copy the entries of array `b` into array `a`. Returns `a`.
"""
function copy_entries!(a,b)
  if a !== b
    copyto!(a,b)
  end
  a
end

"""
    add_entries!(a,b,combine=+)

Perform the operation `combine` element-wise in the entries of arrays `a` and `b`
and store the result in array `a`. Returns `a`.
"""
function add_entries!(a,b,combine=+)
  @assert length(a) == length(b)
  @inbounds for i in eachindex(a)
    a[i] = combine(a[i],b[i])
  end
  a
end

"""
    add_entry!(A,v,i,j,combine=+)

Add an entry given its position and the operation to perform.
"""
function add_entry!(A,v,i::Integer,j::Integer,combine=+)
  aij = A[i,j]
  A[i,j] = combine(aij,v)
end

"""
    add_entry!(A,v,i,combine=+)

Add an entry given its position and the operation to perform.
"""
function add_entry!(A,v,i::Integer,combine=+)
  ai = A[i]
  A[i] = combine(ai,v)
end

"""
    scale_entries!(a,v)

Scale the entries of array `a` with the value `v`. Returns `a`.
"""
function scale_entries!(a,b)
  @inbounds for i in eachindex(a)
    a[i] = b*a[i]
  end
  a
end

# Base.mul!

"""
    muladd!(c,a,b)

Matrix multiply a*b and add to result to c. Returns c.
"""
function muladd!(c,a,b)
  _muladd!(c,a,b)
  c
end

@static if VERSION >= v"1.3"
  function _muladd!(c,a,b)
    mul!(c,a,b,1,1)
  end
else
  function _muladd!(c,a,b)
    @assert length(c) == size(a,1)
    @assert length(b) == size(a,2)
    @inbounds for j in 1:size(a,2)
      for i in 1:size(a,1)
        c[i] += a[i,j]*b[j]
      end
    end
  end
end


text:
using Test

using PlutoPDF

testfile = download(
    "https://raw.githubusercontent.com/fonsp/Pluto.jl/main/sample/Tower%20of%20Hanoi.jl",
)

is_CI = get(ENV, "CI", "no") == "no"

outfile = pluto_to_pdf(testfile; open=is_CI)

@test isfile(outfile)
@test endswith(outfile, ".pdf")

output_dir = get(ENV, "TEST_OUTPUT_DIR", nothing)

if @show(output_dir) isa String
    mkpath(output_dir)
    @assert isdir(output_dir)
    cp(outfile, joinpath(output_dir, "hanoi.pdf"))
end

max_stars_repo_name: JuliaPluto/pluto-export

max_issues_repo_name: JuliaPluto/pluto-export

max_forks_repo_name: JuliaPluto/pluto-export


text:
using LoggingProfiler
function juliaset_pixel(z₀, c)
    z = z₀
    for i in 1:255
        abs2(z)> 4.0 && return (i - 1)%UInt8
        z = z*z + c
    end
    return UInt8(255)
end

function juliaset_column!(img, c, n, j)
    x = -2.0 + (j-1)*4.0/(n-1)
    for i in 1:n
        y = -2.0 + (i-1)*4.0/(n-1)
        @inbounds img[i,j] = juliaset_pixel(x+im*y, c)
    end
    nothing
end

function juliaset_single!(img, c, n)
    for j in 1:n
        juliaset_column!(img, c, n, j)
    end
    nothing
end

function juliaset_static!(img, c, n)
    Threads.@threads for j in 1:n
        juliaset_column!(img, c, n, j)
    end
    nothing
end

function juliaset(x, y, n=1000, method = juliaset_single!, extra...)
    c = x + y*im
    img = Array{UInt8,2}(undef,n,n)
    method(img, c, n, extra...)
    return img
end

# @record juliaset(-0.79, 0.15, 10)
# without any filter this makes
# LoggingProfiler.recorded()
# 144

# blacklist!(LoggingProfiler.timable_list, Core)
# blacklist!(LoggingProfiler.timable_list, Base)
# whitelist!(LoggingProfiler.recursable_list, :juliaset_column!)
# whitelist!(LoggingProfiler.recursable_list, :juliaset_single!)
# whitelist!(LoggingProfiler.timable_list, :juliaset_column!)
# whitelist!(LoggingProfiler.timable_list, :juliaset_single!)
# whitelist!(LoggingProfiler.timable_list, :juliaset_pixel!)
LoggingProfiler.initbuffer!(1000)
@record juliaset(-0.79, 0.15, 10, juliaset_static!)
LoggingProfiler.recorded()
LoggingProfiler.adjustbuffer!()
@record juliaset(-0.79, 0.15, 10, juliaset_static!)
# 53828
LoggingProfiler.export2svg("/tmp/profile.svg")

max_stars_repo_name: pevnak/LoggingProfiler

max_issues_repo_name: pevnak/LoggingProfiler

max_forks_repo_name: pevnak/LoggingProfiler


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script to compile data for some descriptive statistics (VALUES)
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, CSV, XLSX, LinearAlgebra, Statistics, StatsBase, StatsPlots, Dates, Pipe

# other scripts
dir_home = "x:/VIVES/1-Personal/Florian/git/2022_01_comext/src/"
include(dir_home * "functions.jl")

# location of data input/output (io)
dir_io = "C:/Users/u0148308/data/comext/"
dir_dropbox = "C:/Users/u0148308/Dropbox/BREXIT/"

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Raw data
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

df_VLAIO = CSV.read(dir_dropbox * "rawdata/" * "df_VLAIO" * ".csv", DataFrame)
transform!(df_VLAIO, ["PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "DECLARANT_ISO", "TRADE_TYPE", "PARTNER_ISO", "FLOW"] .=> ByRow(string), renamecols=false)
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> convert(Union{Missing, Float64}, x)), renamecols=false)
df_VLAIO.PRODUCT_NC = lpad.(string.(df_VLAIO.PRODUCT_NC), 8, '0') # needs to be done again

# remove missing, only for VALUES
#   - drop missing [VALUE, QUANTITY]: lose 1-16027346/21667553 ~26% observations
#   - drop missing VALUE: lose 1-21666543/21667553 ~ 0.1% observations
subset!(df_VLAIO, :VALUE_IN_EUROS .=> ByRow(x -> !ismissing(x)))

# not needed in case of VALUES!!
# MAD adjustment
#   - cutoff = 3: lose further 1-15177297/21666543 ~30% observations (~30% of total)
#   - cutoff = 4: lose further 1-15812254/21666543 ~27% observations (~27% of total)
# outlier_cutoff = 4
# gdf = groupby(df_VLAIO, "PRODUCT_NC")
# df_VLAIO = transform(gdf, :VALUE_IN_EUROS => MAD_method => :MAD)
# subset!(df_VLAIO, :MAD => ByRow(x -> x < outlier_cutoff))
# df_VLAIO = df_VLAIO[:, Not(:MAD)]

# add EU/WORLD again
EU27 = ["Roemenië", "Griekenland", "Oostenrijk", "Polen", "Duitsland", "Spanje", "Hongarije", "Slovakije", "Italië", "Nederland",
       "Frankrijk", "Letland", "Kroatië", "Cyprus", "Malta", "Litouwen", "Slovenië", "Estland", "Portugal", "Finland", "Tsjechië",
       "Luxemburg", "Zweden", "Denemarken", "Bulgarije", "Ierland", "België"]

# slightly modified functions
function append_WORLD(df::DataFrame)

    # sum over PARTNER
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df, cols_grouping)
    df_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_WORLD.TRADE_TYPE .= "total"
    df_WORLD.PARTNER_ISO .= "WORLD"

    df = vcat(df, df_WORLD)

    return df
end

function append_EU(df::DataFrame, EU::Vector{String})

    # subset and sum over EU ctrys
    df_EU = subset(df, :PARTNER_ISO => ByRow(x -> x in EU))
    cols_grouping = ["DECLARANT_ISO", "PRODUCT_NC", "PRODUCT_NC_LAB1", "PRODUCT_NC_LAB2", "TRADE_TYPE", "FLOW", "PERIOD"]
    gdf = groupby(df_EU, cols_grouping)
    df_EU = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

    df_EU.TRADE_TYPE .= "intra"
    df_EU.PARTNER_ISO .= "EU"

    df = vcat(df, df_EU)

    return df
end

df_VLAIO = append_WORLD(df_VLAIO)
df_VLAIO = append_EU(df_VLAIO, EU27)

# double check if no zeros/missing introduced by WORLD/EU
transform!(df_VLAIO, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> ByRow(x -> ifelse(iszero(x) | ismissing(x), missing, x)), renamecols=false)
subset!(df_VLAIO, :VALUE_IN_EUROS => ByRow(x -> !ismissing(x))) # lose 0% observations

cols_subset = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "PRODUCT_NC", "VALUE_IN_EUROS", "QUANTITY_IN_KG"]
df_VLAIO = df_VLAIO[:, cols_subset]

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 1
#   - Export/import value
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# formatting
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)

# need to sort before computing HP and 3MMA
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :VALUE_IN_EUROS => (x -> HP(x, λ)) => :VALUE_IN_EUROS_HP)

# plotting
#   TO DO:
#       - translate title!
line_color = [palette(:tab10)[1:3]; palette(:tab10)[7]; palette(:tab10)[4]]';
tick_years = Date.(2014:2021)
DateTick = Dates.format.(tick_years, "yyyy")

for flow in ["imports", "exports"]

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :VALUE_IN_EUROS_HP/1e9, linecolor=line_color, xticks=false,
            group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="euro (miljard)", title="Flemish "*flow*": values (HP, λ=$λ)")
    plot!(xticks=(tick_years, DateTick))
    vline!([Date(2016,6,23)], label="", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/formatted/" * "fig1_" * flow * "_values" * "_HP" * ".png")

end

df_export = subset(df, :FLOW => ByRow(x -> x == "exports"))[:,Not([:VALUE_IN_EUROS])]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure1" * ".xlsx", df_export, overwrite=true)

df_export = subset(df, :FLOW => ByRow(x -> x == "imports"))[:,Not([:VALUE_IN_EUROS])]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure2" * ".xlsx", df_export, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 3
#   - Export/import share
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "WORLD"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

df_WORLD = subset(df, :PARTNER_ISO => ByRow(x -> x == "WORLD"))
cols_grouping = ["DECLARANT_ISO", "FLOW", "PERIOD"]
gdf = groupby(df_WORLD, cols_grouping)
df_WORLD = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_WORLD)

df = leftjoin(subset(df, :PARTNER_ISO => ByRow(x -> x != "WORLD")), df_WORLD, on=[:DECLARANT_ISO, :FLOW, :PERIOD])
transform!(df, [:VALUE_IN_EUROS, :VALUE_WORLD] => ByRow((x,s) -> x/s) => :SHARE_VALUE)

# Date
transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :SHARE_VALUE => (x -> HP(x, λ)) => :SHARE_VALUE_HP)

# plotting
#   TO DO:
#       - translate title!
line_color = [:red; :green]';
tick_years = Date.(2014:2021)
DateTick = Dates.format.(tick_years, "yyyy")

# VLA vis-a-vis GB
p = @df subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk")) plot(:DATE, :SHARE_VALUE_HP*100, linecolor=[palette(:tab10)[3] palette(:tab10)[4]], xticks=false,
        group=:FLOW, lw=2, legend=:topright, ylabel="percentages", title="Trade share (values, λ=$λ): \n Flanders vis-a-vis Great Britain")
    plot!(xticks=(tick_years, DateTick))
    vline!([Date(2016,6,23)], label="", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="", color=:black, lw=1.5, ls=:dot) # trans end
savefig(p, dir_dropbox * "results/images/VLAIO/formatted/" * "fig3_" * "VLA_GB" * "_value_share" * "_HP" * ".png") # export image

df_export = subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"))[:,["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD", "DATE", "SHARE_VALUE_HP"]]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure3" * ".xlsx", df_export, overwrite=true)

# --------------
# Figure 8
#   - Export/import share of top 10 industries (2 digit NACE)
# --------------

df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))
transform!(df, :PERIOD => ByRow(x -> Date(string(x)[1:4], DateFormat("yyyy"))) => :YEAR)

# CPA classification
df_CN_PROD = DataFrame(XLSX.readtable(dir_dropbox* "results/correspondence/" * "table_CN_CPA21" * ".xlsx", "Sheet1")...)
transform!(df_CN_PROD, names(df_CN_PROD) .=> ByRow(string), renamecols=false)

df_join = leftjoin(df, df_CN_PROD, on=:PRODUCT_NC)
subset!(df_join, :CPA21_lab => ByRow(x -> !ismissing(x))) # lose 1-7019070/7260237 ~ 3%

# aggregate over products and country
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "CPA21_2digits" , "CPA21_lab", "CPA21_lab_custom"]
gdf = groupby(df_join, cols_grouping)
df = combine(gdf, :VALUE_IN_EUROS => sum, renamecols=false)

# aggregate over products
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "YEAR", "FLOW"]
gdf = groupby(df_join, cols_grouping)
df_total_partner = combine(gdf, :VALUE_IN_EUROS => sum => :VALUE_TOTAL)
transform!(df_total_partner, :YEAR => ByRow(x -> string(x)[1:4]), renamecols=false)

# subset products top 10 and produce share
transform!(df, :YEAR => ByRow(x -> string(x)[1:4]), renamecols=false)
top10_exports = first(sort(subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"), :YEAR => ByRow(x -> x == "2021"),
     :FLOW => ByRow(x -> x == "exports")), order(:VALUE_IN_EUROS, rev=true)).CPA21_2digits, 10)
top10_imports = first(sort(subset(df, :PARTNER_ISO => ByRow(x -> x == "Verenigd Koninkrijk"), :YEAR => ByRow(x -> x == "2021"),
     :FLOW => ByRow(x -> x == "imports")), order(:VALUE_IN_EUROS, rev=true)).CPA21_2digits, 10)

# join
df = leftjoin(df, df_total_partner, on=[:DECLARANT_ISO, :PARTNER_ISO, :FLOW, :YEAR])
transform!(df, [:VALUE_IN_EUROS, :VALUE_TOTAL] => ByRow((x, y) -> x/y) => :SHARE_VALUE)
transform!(df, [:CPA21_2digits, :CPA21_lab_custom] => ByRow((x,y) -> x * " " * y) => :legend)

sort!(df)

# plotting
#   TO DO:
#       - translate title!
tick_years = Date.(2014:2021)
DateTick = Dates.format.(tick_years, "yyyy")
line_color = palette(:tab10)[1:10]';

flow = "exports"
p = @df subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in top10_exports)) groupedbar(:YEAR, :SHARE_VALUE*100, color=line_color,
    group=:legend, bar_position=:stack, legend=:outertopleft, ylabel="percentages", title="Flemish "* flow*" with GB: \n Top 10 industries (2021)", xrotation=35)
savefig(p, dir_dropbox * "results/images/VLAIO/formatted/" * "fig8_" * flow * "_values_annual" * "_top10_2021" * ".png") # export image dropbox

df_export = subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in top10_exports))[:,["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "legend", "SHARE_VALUE"]]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure6" * ".xlsx", df_export, overwrite=true)

flow = "imports"
p = @df subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in top10_imports)) groupedbar(:YEAR, :SHARE_VALUE*100, color=line_color,
    group=:legend, bar_position=:stack, legend=:outertopleft, ylabel="percentages", title="Flemish "* flow*" with GB: \n Top 10 industries (2021)", xrotation=35)
savefig(p, dir_dropbox * "results/images/VLAIO/formatted/" * "fig8_" * flow * "_values_annual" * "_top10_2021" * ".png") # export image dropbox

df_export = subset(df, :FLOW => ByRow(x -> x == flow), :CPA21_2digits => ByRow(x -> x in top10_imports))[:,["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "YEAR", "legend", "SHARE_VALUE"]]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure7" * ".xlsx", df_export, overwrite=true)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Figure 5
#   - Export/import number of products
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# prepare data
df = copy(df_VLAIO)
partners = ["Verenigd Koninkrijk", "Duitsland", "Nederland", "Frankrijk", "Italië"]
subset!(df, :DECLARANT_ISO => ByRow(x -> x == "Vlaanderen"), :PARTNER_ISO => ByRow(x -> x in partners))

cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW", "PERIOD"]
gdf = groupby(df, cols_grouping)
df = combine(gdf, nrow => :COUNT)

transform!(df, :PERIOD => ByRow(x -> Date(string(x), DateFormat("yyyymm"))) => :DATE)
sort!(df)

# HP filter
λ = 20
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :COUNT => (x -> HP(x, λ)) => :COUNT_HP)

# plotting
#   TO DO:
#       - translate title!
line_color = [palette(:tab10)[1:3]; palette(:tab10)[7]; palette(:tab10)[4]]';
tick_years = Date.(2014:2021)
DateTick = Dates.format.(tick_years, "yyyy")

for flow in ["imports", "exports"]

    # count
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :COUNT, linecolor=line_color, xticks=false,
            group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="number of products", title="Flemish "* flow*": number of products")
    plot!(xticks=(tick_years, DateTick))
    vline!([Date(2016,6,23)], label="", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/formatted/" * "fig5_" * flow * "_product_count" * ".png")

    # HP
    p = @df subset(df, :FLOW => ByRow(x -> x == flow)) plot(:DATE, :COUNT_HP, linecolor=line_color, xticks=false,
        group=:PARTNER_ISO, lw=2, legend=:topleft, ylabel="number of products", title="Flemish "* flow*": number of products (HP, λ=$λ)")
    plot!(xticks=(tick_years, DateTick))
    vline!([Date(2016,6,23)], label="", color=:black, lw=1, ls=:solid) # refer
    vline!([Date(2020,01,31)], label="", color=:black, lw=1, ls=:dash) # exit
    vline!([Date(2020,12,31)], label="", color=:black, lw=1.5, ls=:dot) # trans end
    savefig(p, dir_dropbox * "results/images/VLAIO/formatted/" * "fig5_" * flow * "_product_count_HP" * ".png")

end

df_export = subset(df, :FLOW => ByRow(x -> x == "exports"))[:,Not([:COUNT_HP])]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure4" * ".xlsx", df_export, overwrite=true)

df_export = subset(df, :FLOW => ByRow(x -> x == "imports"))[:,Not([:COUNT_HP])]
XLSX.writetable(dir_dropbox * "results/images/VLAIO/formatted/" * "data_figure5" * ".xlsx", df_export, overwrite=true)


text:
using Random, StatsBase
import StatsBase: sample

function sample(rng::AbstractRNG, df::DataFrame, n;
    replace=true, ordered=false)
  indxs = sample(rng,
    1:size(df,1),
    n,
    replace=replace, ordered=ordered)
  df[indxs, :]
end

function sample(df::DataFrame, n;
    replace=true, ordered=false)
  indxs = sample(Random.GLOBAL_RNG,
    1:size(df,1),
    n,
    replace=replace, ordered=ordered)
  df[indxs, :]
end


text:
abstract type FileIterator end

"""
    interleavestrings(x::Vector{<:AbstractString}, y::Vector{<:AbstractString})

Interleaves two vectors of strings. `x` comes before `y`.
Shout out: Thanks, Don MacMillen from Slack!

"""
function interleavestrings(x::Vector{<:AbstractString}, y::Vector{<:AbstractString})
    lenx, leny = length(x), length(y)
    largest = max(lenx, leny)
    result = fill("", 2*largest)
    result[1:2:2*lenx] = x
    result[2:2:2*leny] = y
    return join(result)
end

struct Expand <: FileIterator
    statictxt::Vector{String}
    mapping::Any
    productiter::Base.Iterators.ProductIterator
    n::Int
end

"""
    Expand( str::String, replace_map::Dict{String,Vector{String}} )

Expand is             an iterator that replaces keywords in a string, `str`, with all
permutations in the `replace_map`. To define a keyword it must be enclosed in
curly brackets `{...}`.

ToDo: Get dictionary type correct...
"""
function Expand( str::String, replace_map::Dict)#{String, Vector{Any}}) #where {T} #<:Union{String, Real}}
    locations, items, cuts = [], [], []
    #TODO: Char-wise/Trie search instead could be more performant?
    for ( curkey, item ) in replace_map
        firstunitofallmatches = first.( findall( "{" * curkey * "}", str ) ) .=> curkey
        @assert(length(firstunitofallmatches) <= 1, "Cannot use keyword $curkey twice in expand statement.")
        if length(firstunitofallmatches) == 1
            locations = vcat( locations, firstunitofallmatches )
            push!( items, all(isa.(item, Real)) ? string.(item) : item )
        end
    end
    #sort the order of the items found in the string by which come first!
    idx = sortperm(locations, by = x -> first( x ) )
    locations = locations[idx]
    #Handle edgecase where string starts with keyword to be replaced
    ( lastloc, tag ) = first( locations )
    firstcut = ( lastloc[ 1 ] > 1 ) ? str[ 1 : ( lastloc - 1 ) ] : ""
    push!( cuts, firstcut )
    for item in 1 : ( length( locations ) - 1 )
        ( loc, tag ) = locations[ item ]
        ( nextloc, tag2 ) = locations[ item + 1 ]
        push!( cuts,  str[ ( loc + length( tag ) + 2 ) : ( nextloc - 1 ) ] )
    end
    #handle end of string, edge case is implicit
    ( lastloc, tag ) = locations[ end ]
    push!( cuts,  str[ ( lastloc + length( tag ) + 2 ) : end] )
    total_len = prod( length.( items[ idx ] ) )
    return Expand( cuts, locations, Iterators.product( items[idx]... ), total_len )
end

Base.length(ex::Expand) = ex.n

function Base.iterate( iter::Expand, state = ( nothing ) )
    if isnothing(state)
        subiter, state = collect( iterate( iter.productiter ) )
    else
        nextiter = iterate( iter.productiter, state )
        isnothing( nextiter ) && return nothing
        subiter, state = collect( nextiter )
    end
    return interleavestrings( iter.statictxt, collect( subiter ) ), state
end

mutable struct VariableGlob
    parsedsegments::Vector{String}
    parsedvariables::Vector{String}
    parsed_values::Vector{Any}#TODO: Type specific...
end

"""
    VariableGlob( path::String )

Given an input `path`, parse for variable names enclosed in curly brackets.
Returns a tuple of nonvariable elements of the path, and variable elements.

"""
function VariableGlob( path::String, expression::String )
    strlen          = length( expression )
    isenclosed      = false
    segments, variablenames = [], []
    cursor,   variable_cursor = 1, 1
    #Parse the text around curly brackets, and the text inside curly brackets into
    # 2 different vectors.
    for (i, chr) in enumerate( collect( expression ) )
        if chr == '{'
            if !isenclosed
                if i > 1
                    push!(segments, expression[ cursor : (i - 1 ) ] )
                    variable_cursor = i
                else
                    push!(segments, "" )
                end
                isenclosed = true
            else
                @error( "Enclosing symbol `{` found without matching `}`. \n Only single nested wild cards supported." )
            end
        elseif chr == '}'
            if isenclosed
                isenclosed = false
                cursor = i + 1
                push!( variablenames, expression[ ( variable_cursor + 1 ) : ( i - 1 ) ] )
            else
                @error( "Enclosing symbol `}` found without matching `{`." )
            end
        elseif i == strlen
            if !isenclosed
                push!( segments, expression[ cursor : strlen ] )
            else
                @error( "Enclosing symbol `{` found without matching `}`." )
            end
        end
    end
    #Safety checks
    @assert( length(unique(variablenames)) == length(variablenames), "Replicate variable names found in string. Cannot proceed." )
    @assert( all( length.(variablenames) .>= 1 ), "Variable name cannot be less then 1 character long!" )
    @assert( ( length(variablenames) < length(segments) ), "Filename cannot end with a wildcard." )
    #Clean up our understanding of the root directory...
    filematches = glob( join(segments, "*"), path  )
    #Now that we have a list of matching file names let's get their metadata
    segmentlengths = length.( segments )
    variables_count = length(variablenames)
    template = [ ]
    templatecount = 0
    for filematch in filematches
        idx = length(path)
        matchlen = length(filematch)
        tmpdict = Dict()
        for ( s, seglen ) in enumerate( segmentlengths )
            if s <= variables_count
                idx += seglen
                findnextseg = matchlen
                if s < length(segments)
                    findnextsegmatch = match( Regex( "(" * segments[s + 1] * ")s?" ), filematch)
                    if !isnothing(findnextsegmatch)
                        findnextseg = findnextsegmatch.offset - 1
                    end
                end
                tmpdict[ variablenames[ s ] ] = filematch[ (idx+1):(findnextseg) ]
                idx = findnextseg
            end
        end
        push!(template, tmpdict)
    end
    return VariableGlob( segments, variablenames, template )
end

max_stars_repo_name: caseykneale/MrPhelps

max_issues_repo_name: caseykneale/MrPhelps

max_forks_repo_name: caseykneale/MrPhelps


text:
struct PipelineColorBlendAttachmentState
    mHandle::vk.VkPipelineColorBlendAttachmentState

    function PipelineColorBlendAttachmentState(;
        blendEnable::vk.VkBool32 = VkExt.VK_FALSE,
        srcColorBlendFactor::vk.VkBlendFactor = vk.VK_BLEND_FACTOR_ZERO,
        dstColorBlendFactor::vk.VkBlendFactor = vk.VK_BLEND_FACTOR_ZERO,
        colorBlendOp::vk.VkBlendOp = vk.VK_BLEND_OP_ADD,
        srcAlphaBlendFactor::vk.VkBlendFactor = vk.VK_BLEND_FACTOR_ZERO,
        dstAlphaBlendFactor::vk.VkBlendFactor = vk.VK_BLEND_FACTOR_ZERO,
        alphaBlendOp::vk.VkBlendOp = vk.VK_BLEND_OP_ADD,
        colorWriteMask::vk.VkColorComponentFlags = 0
    )
        new(
            vk.VkPipelineColorBlendAttachmentState(
                blendEnable, #::VkBool32
                srcColorBlendFactor, #::VkBlendFactor
                dstColorBlendFactor, #::VkBlendFactor
                colorBlendOp, #::VkBlendOp
                srcAlphaBlendFactor, #::VkBlendFactor
                dstAlphaBlendFactor, #::VkBlendFactor
                alphaBlendOp, #::VkBlendOp
                colorWriteMask #::VkColorComponentFlags
            )
        )
    end
end
@class PipelineColorBlendAttachmentState

struct PipelineColorBlendStateCreateInfo
    mHandleRef::Ref{vk.VkPipelineColorBlendStateCreateInfo}
    mPreserve::Vector{Any}

    function PipelineColorBlendStateCreateInfo(;
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineColorBlendStateCreateFlags = vk.VkFlags(0),
        logicOpEnable::vk.VkBool32 = VkExt.VK_FALSE,
        logicOp::vk.VkLogicOp = vk.VK_LOGIC_OP_COPY,
        attachments::Vector{vk.VkPipelineColorBlendAttachmentState} = Vector{vk.VkPipelineColorBlendAttachmentState}(),
        blendConstants::NTuple{4, Cfloat} = (Cfloat(0), Cfloat(0), Cfloat(0), Cfloat(0))
    )

        this = new(Ref(vk.VkPipelineColorBlendStateCreateInfo(
            vk.VK_STRUCTURE_TYPE_PIPELINE_COLOR_BLEND_STATE_CREATE_INFO, #sType::VkStructureType
            pNext, #::Ptr{Cvoid}
            flags, #::VkPipelineColorBlendStateCreateFlags
            logicOpEnable, #::VkBool32
            logicOp, #::VkLogicOp
            length(attachments), #::UInt32
            pointer(attachments), #::Ptr{VkPipelineColorBlendAttachmentState}
            blendConstants #::NTuple{4, Cfloat}
        )), [attachments])
    end
end

function handle(this::PipelineColorBlendAttachmentState)::vk.VkPipelineColorBlendAttachmentState
    return this.mHandle
end

function handleRef(this::PipelineColorBlendStateCreateInfo)::Ref{vk.VkPipelineColorBlendStateCreateInfo}
    this.mHandleRef
end

function addNoBlend(outStates::Vector{vk.VkPipelineColorBlendAttachmentState})
    state = PipelineColorBlendAttachmentState(
                colorWriteMask = (vk.VK_COLOR_COMPONENT_R_BIT | vk.VK_COLOR_COMPONENT_G_BIT |
                vk.VK_COLOR_COMPONENT_B_BIT | vk.VK_COLOR_COMPONENT_A_BIT)
            ).handle()
    push!(outStates, state)
end

function addTransparencyBlend(outStates::Vector{vk.VkPipelineColorBlendAttachmentState})
    state = PipelineColorBlendAttachmentState(
                blendEnable = VkExt.VK_TRUE,
                colorWriteMask = (vk.VK_COLOR_COMPONENT_R_BIT | vk.VK_COLOR_COMPONENT_G_BIT |
                                  vk.VK_COLOR_COMPONENT_B_BIT | vk.VK_COLOR_COMPONENT_A_BIT),
                alphaBlendOp = vk.VK_BLEND_OP_ADD,
                srcAlphaBlendFactor = vk.VK_BLEND_FACTOR_ONE,
                dstAlphaBlendFactor = vk.VK_BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,
                colorBlendOp = vk.VK_BLEND_OP_ADD,
                srcColorBlendFactor = vk.VK_BLEND_FACTOR_SRC_ALPHA,
                dstColorBlendFactor = vk.VK_BLEND_FACTOR_ONE_MINUS_SRC_ALPHA
            ).handle()
    push!(outStates, state)
end


text:
import Pkg
const PT = Pkg.Types

Pkg.activate(pwd())             # current directory as the project
ctx = PT.Context()
pkg = ctx.env.pkg
if pkg ≡ nothing
    @error "Not in a package, I won't generate REQUIRE."
    exit(1)
else
    @info "found package" pkg = pkg
end

deps = PT.get_deps(ctx)
non_std_deps = sort(collect(setdiff(keys(deps), values(ctx.stdlibs))))

open("REQUIRE", "w") do io
    for d in non_std_deps
        println(io, d)
        @info "listing $d"
    end
end


text:
using GeometricAlgebra
using Test
using SafeTestsets

@testset "GeometricAlgebra.jl" begin
    @safetestset "Implementation" begin include("implementations.jl") end
    @safetestset "Identities in 𝓖₄" begin include("algebras/r4.jl") end
    @safetestset "Identities in 𝓖₃,₁" begin include("algebras/spacetime.jl") end
    @safetestset "3D Conformal Geometric Algebra" begin include("algebras/3d_cga.jl") end
    @testset "Aqua" begin include("aqua.jl") end
end


text:
import Base: MethodList

stripparams(t) = replace(t, r"\{([A-Za-z, ]*?)\}", "")

function methodarray(mt::MethodList)
  defs = collect(mt)
  file(m) = m.file |> string |> basename
  line(m) = m.line
  sort!(defs, lt = (a, b) -> file(a) == file(b) ?
                               line(a) < line(b) :
                               file(a) < file(b))
  return defs
end

methodarray(x) = methodarray(methods(x))

interpose(xs, y) = map(i -> iseven(i) ? xs[i÷2] : y, 2:2length(xs))

function view(m::Method)
  tv, decls, file, line = Base.arg_decl_parts(m)
  params = [span(c(x, isempty(T) ? "" : "::", strong(stripparams(T)))) for (x, T) in decls[2:end]]
  params = interpose(params, ", ")
  span(c(string(m.name),
         "(", params..., ")")),
  file == :null ? "not found" : Atom.baselink(string(file), line)
end

@render i::Inline m::Method begin
  sig, link = view(m)
  r(x) = render(i, x)
  span(c(r(sig), " at ", r(link)))
end

# TODO: factor out table view
@render i::Inline m::MethodList begin
  ms = methodarray(m)
  isempty(ms) && return "$(m.mt.name) has no methods."
  r(x) = render(i, x)
  length(ms) == 1 && return r(ms[1])
  Tree(span(c(span(".support.function", string(m.mt.name)),
              " has $(length(ms)) methods:")),
       [table(".methods", [tr(td(c(r(a))), td(c(r(b)))) for (a, b) in map(view, ms)])])
end


text:
# Tests the Arnoldi relation AV = VH when expanding the search subspace

using Test, LinearAlgebra, SparseArrays
using ArnoldiMethod: reinitialize!, Arnoldi, iterate_arnoldi!

@testset "Initialization" begin
    arnoldi = Arnoldi{Float64}(5, 3)
    reinitialize!(arnoldi)
    @test norm(arnoldi.V[:, 1]) ≈ 1
end

@testset "Arnoldi Factorization" begin
    n = 10
    max = 6
    for T in (Float64, BigFloat)
        A = sprand(T, n, n, .1) + I

        arnoldi = Arnoldi{T}(n, max)
        reinitialize!(arnoldi)
        V, H = arnoldi.V, arnoldi.H

        # Do a few iterations
        iterate_arnoldi!(A, arnoldi, 1:3)
        @test A * V[:,1:3] ≈ V[:,1:4] * H[1:4,1:3]
        @test norm(V[:,1:4]' * V[:,1:4] - I) < sqrt(eps(T))/100

        # Do the rest of the iterations.
        iterate_arnoldi!(A, arnoldi, 4:max)
        @test A * V[:,1:max] ≈ V * H
        @test norm(V' * V - I) < sqrt(eps(T))/100
    end
end

@testset "Invariant subspace" begin
    # Generate a block-diagonal matrix A
    for T in (Float64, BigFloat)
        A = [rand(T,4,4)  zeros(T,4,4);
            zeros(T,4,4) rand(T,4,4)]

        # and an initial vector [1; 0; ... 0]
        vh = Arnoldi{T}(8, 5)
        V, H = vh.V, vh.H
        V[:,1] .= zero(T)
        V[1,1] = one(T)

        # Then {v, Av, A²v, A³v}
        # is an invariant subspace
        iterate_arnoldi!(A, vh, 1:5)

        @test norm(V' * V - I) < sqrt(eps(T))/100
        @test iszero(H[5, 4])
    end
end


text:
abstract type AbstractEchoStateNetwork <: AbstractReservoirComputer end
struct ESN{I,S,V,N,T,O,M,B,ST,W,IS} <: AbstractEchoStateNetwork
    res_size::I
    train_data::S
    variation::V
    nla_type::N
    input_matrix::T
    reservoir_driver::O
    reservoir_matrix::M
    bias_vector::B
    states_type::ST
    washout::W
    states::IS
end

"""
    Default()

Sets the type of the ESN as the standard model. No parameters are needed.
"""
struct Default <: AbstractVariation end
struct Hybrid{T,K,O,I,S,D} <: AbstractVariation
    prior_model::T
    u0::K
    tspan::O
    dt::I
    datasize::S
    model_data::D
end

"""
    Hybrid(prior_model, u0, tspan, datasize)

Given the model parameters returns an ```Hybrid``` variation of the ESN. This entails
a different training and prediction. Construction based on [1].

[1] Jaideep Pathak et al. "Hybrid Forecasting of Chaotic Processes: Using Machine
Learning in Conjunction with a Knowledge-Based Model" (2018)
"""
function Hybrid(prior_model, u0, tspan, datasize)
    trange = collect(range(tspan[1], tspan[2], length = datasize))
    dt = trange[2]-trange[1]
    tsteps = push!(trange, dt + trange[end])
    tspan_new = (tspan[1], dt+tspan[2])
    model_data = prior_model(u0, tspan_new, tsteps)
    return Hybrid(prior_model, u0, tspan, dt, datasize, model_data)
end

"""
    ESN(train_data;
        variation = Default(),
        input_layer = DenseLayer(),
        reservoir = RandSparseReservoir(),
        bias = NullLayer(),
        reservoir_driver = RNN(),
        nla_type = NLADefault(),
        states_type = StandardStates())
    (esn::ESN)(prediction::AbstractPrediction,
        output_layer::AbstractOutputLayer;
        initial_conditions=output_layer.last_value,
        last_state=esn.states[:, end])

Constructor for the Echo State Network model. It requires the reservoir size as the input
and the data for the training. It returns a struct ready to be trained with the states
already harvested.

After the training this struct can be used for the prediction following the second
function call. This will take as input a prediction type and the output layer from the
training. The ```initial_conditions``` and ```last_state``` parameters can be left as
they are, unless there is a specific reason to change them. All the components are
detailed in the API documentation. More examples are given in the general documentation.
"""
function ESN(
    train_data;
    variation = Default(),
    input_layer = DenseLayer(),
    reservoir = RandSparseReservoir(100),
    bias = NullLayer(),
    reservoir_driver = RNN(),
    nla_type = NLADefault(),
    states_type = StandardStates(),
    washout = 0,
    matrix_type = typeof(train_data)
)
    if variation isa Hybrid
        train_data = vcat(train_data, variation.model_data[:, 1:end-1])
    end

    if states_type isa AbstractPaddedStates
        in_size = size(train_data, 1) + 1
        train_data = vcat(
            Adapt.adapt(matrix_type, ones(1, size(train_data, 2))), train_data
        )
    else
        in_size = size(train_data, 1)
    end

    input_matrix, reservoir_matrix, bias_vector, res_size = obtain_layers(
        in_size, input_layer, reservoir, bias; matrix_type=matrix_type
    )

    inner_res_driver = reservoir_driver_params(reservoir_driver, res_size, in_size)
    states = create_states(
        inner_res_driver, train_data, washout, reservoir_matrix, input_matrix, bias_vector
    )
    train_data = train_data[:, washout+1:end]

    ESN(sum(res_size), train_data, variation, nla_type, input_matrix,
        inner_res_driver, reservoir_matrix, bias_vector, states_type, washout,
        states)
end

#shallow esn construction
function obtain_layers(
    in_size,
    input_layer,
    reservoir,
    bias;
    matrix_type = Matrix{Float64}
)
    input_res_size = get_ressize(reservoir)
    input_matrix = create_layer(
        input_layer, input_res_size, in_size, matrix_type=matrix_type
    )
    res_size = size(input_matrix, 1) #WeightedInput actually changes the res size
    reservoir_matrix = create_reservoir(reservoir, res_size, matrix_type=matrix_type)
    @assert size(reservoir_matrix, 1) == res_size
    bias_vector = create_layer(bias, res_size, 1, matrix_type=matrix_type)
    return input_matrix, reservoir_matrix, bias_vector, res_size
end

#deep esn construction
#there is a bug going on with WeightedLayer in this construction.
#it works for eny other though
function obtain_layers(
    in_size,
    input_layer,
    reservoir::Vector,
    bias;
    matrix_type=Matrix{Float64}
)
    esn_depth = length(reservoir)
    input_res_sizes = [get_ressize(reservoir[i]) for i=1:esn_depth]
    in_sizes = zeros(Int, esn_depth)
    in_sizes[2:end] = input_res_sizes[1:end-1]
    in_sizes[1] = in_size

    if input_layer isa Array
        input_matrix = [create_layer(
            input_layer[j], input_res_sizes[j], in_sizes[j], matrix_type=matrix_type
        ) for j in 1:esn_depth]
    else
        _input_layer = fill(input_layer, esn_depth)
        input_matrix = [create_layer(
            _input_layer[k], input_res_sizes[k], in_sizes[k], matrix_type=matrix_type
        ) for k in 1:esn_depth]
    end

    res_sizes = [get_ressize(input_matrix[j]) for j in 1:esn_depth]
    reservoir_matrix = [create_reservoir(
        reservoir[k], res_sizes[k], matrix_type=matrix_type
    ) for k in 1:esn_depth]

    if bias isa Array
        bias_vector = [create_layer(
            bias[j], res_sizes[j], 1, matrix_type=matrix_type
        ) for j in 1:esn_depth]
    else
        _bias = fill(bias, esn_depth)
        bias_vector = [create_layer(
            _bias[k], res_sizes[k], 1, matrix_type=matrix_type
        ) for k in 1:esn_depth]
    end

    return input_matrix, reservoir_matrix, bias_vector, res_sizes
end

function (esn::ESN)(
    prediction::AbstractPrediction,
    output_layer::AbstractOutputLayer;
    initial_conditions=output_layer.last_value,
    last_state=esn.states[:, [end]]
)
    variation = esn.variation
    pred_len = prediction.prediction_len

    if variation isa Hybrid
        model = variation.prior_model
        predict_tsteps = [variation.tspan[2]+variation.dt]
        [append!(predict_tsteps, predict_tsteps[end]+variation.dt) for i in 1:pred_len]
        tspan_new = (variation.tspan[2]+variation.dt, predict_tsteps[end])
        u0 = variation.model_data[:, end]
        model_pred_data = model(u0, tspan_new, predict_tsteps)[:, 2:end]
        return obtain_esn_prediction(
            esn, prediction, last_state, output_layer, model_pred_data; initial_conditions=initial_conditions
        )
    else
        return obtain_esn_prediction(
            esn, prediction, last_state, output_layer; initial_conditions=initial_conditions
        )
    end
end

#training dispatch on esn
"""
    train(esn::AbstractEchoStateNetwork, target_data, training_method=StandardRidge(0.0))

Training of the built ESN over the ```target_data```. The default training method is
RidgeRegression. The output is an ```OutputLayer``` object to be fed at the esn call
for the prediction.
"""
function train(
    esn::AbstractEchoStateNetwork,
    target_data,
    training_method=StandardRidge(0.0)
)
    variation = esn.variation

    if esn.variation isa Hybrid
        states = vcat(esn.states, esn.variation.model_data[:, 2:end])
    else
        states=esn.states
    end
    states_new = esn.states_type(esn.nla_type, states, esn.train_data[:,1:end])

    return _train(states_new, target_data, training_method)
end

function pad_esnstate(variation::Hybrid, states_type, x_pad, x, model_prediction_data)
    x_tmp = vcat(x, model_prediction_data)
    x_pad = pad_state!(states_type, x_pad, x_tmp)
end

function pad_esnstate!(variation, states_type, x_pad, x, args...)
    x_pad = pad_state!(states_type, x_pad, x)
end

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
#
# Implementation of AtomsBase interface in an array-of-structs style
#
export FlexibleSystem

struct FlexibleSystem{D, S, L<:Unitful.Length} <: AbstractSystem{D}
    particles::AbstractVector{S}
    box::SVector{D, SVector{D, L}}
    boundary_conditions::SVector{D, BoundaryCondition}
end

function FlexibleSystem(
    particles::AbstractVector{S},
    box::AbstractVector{<:AbstractVector{L}},
    boundary_conditions::AbstractVector{BC}
) where {BC<:BoundaryCondition, L<:Unitful.Length, S}
    D = length(box)
    if !all(length.(box) .== D)
        throw(ArgumentError("Box must have D vectors of length D"))
    end
    FlexibleSystem{D, S, L}(particles, box, boundary_conditions)
end

# Update constructor
function FlexibleSystem(system::AbstractSystem;
                        particles=nothing, atoms=nothing,
                        box=bounding_box(system),
                        boundary_conditions=boundary_conditions(system))
    particles = something(particles, atoms, collect(system))
    FlexibleSystem(particles, box, boundary_conditions)
end

function Base.show(io::IO, system::FlexibleSystem)
    print(io, "FlexibleSystem")
    show_system(io, system)
end

bounding_box(sys::FlexibleSystem)        = sys.box
boundary_conditions(sys::FlexibleSystem) = sys.boundary_conditions
species_type(sys::FlexibleSystem{D, S, L}) where {D, S, L} = S

Base.size(sys::FlexibleSystem)   = size(sys.particles)
Base.length(sys::FlexibleSystem) = length(sys.particles)
Base.getindex(sys::FlexibleSystem, i::Integer) = getindex(sys.particles, i)


text:
__precompile__(true)

module UnalignedVectors

export UnalignedVector, MaybeUnalignedVector, unaligned_reinterpret

"""
    v = UnalignedVector{T}(a::Vector{UInt8})

Create a vector with element type `T` from a memory buffer of bytes
(`UInt8`). In contrast with `reinterpret`, this allows array creation
even if `a` does not have proper pointer alignment for `T`.
"""
struct UnalignedVector{T} <: AbstractArray{T,1}
    a::Vector{UInt8}
    len::Int

    function UnalignedVector{T}(a::Vector{UInt8}) where {T}
        len = length(a) ÷ sizeof(T)
        len*sizeof(T) == length(a) || throw(DimensionMismatch("length(a) should be a multiple of $(sizeof(T)), got $(length(a))"))
        new{T}(a, len)
    end
end

MaybeUnalignedVector{T} = Union{UnalignedVector{T},Vector{T}}

Base.IndexStyle(::Type{UnalignedVector{T}}) where {T} = IndexLinear()
@inline Base.length(a::UnalignedVector) = a.len
@inline Base.size(a::UnalignedVector) = (length(a),)

@inline function Base.getindex(a::UnalignedVector{T}, i::Int) where {T}
    @boundscheck checkbounds(a, i)
    unsafe_load(Ptr{T}(pointer(a.a)), i)
end
@inline function Base.setindex!(a::UnalignedVector{T}, val, i::Int) where {T}
    @boundscheck checkbounds(a, i)
    unsafe_store!(Ptr{T}(pointer(a.a)), val, i)
end

"""
    v = unaligned_reinterpret(T, a::Vector{UInt8})

Reinterprets `a` as an `UnalignedVector{T}`, unless `T == UInt8` in
which case `a` is returned.
"""
unaligned_reinterpret(::Type{T}, a::Vector{UInt8}) where {T} = UnalignedVector{T}(a)
unaligned_reinterpret(::Type{UInt8}, a::Vector{UInt8}) = a

function Base.reinterpret(::Type{T}, a::UnalignedVector, dims::Dims) where {T}
    reshape(unaligned_reinterpret(T, a.a), dims)
end
function Base.reinterpret(::Type{T}, a::Base.ReshapedArray{S,N,UnalignedVector{S}}, dims::Dims) where {T,S,N}
    reshape(unaligned_reinterpret(T, parent(a).a), dims)
end

end


text:
using DataDeps

register(DataDep(
  "congressional-voting-records",
  "http://archive.ics.uci.edu/ml/datasets/Congressional+Voting+Records",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/voting-records/house-votes-84.data",
  "c87c14110a5ba91d4a1e313ec7392824458152bf071fa5f5452340488337936e",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:17,
        categoric_indices=2:17,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using FileTools
using Base.Test

@testset "FileTools" begin
  @testset "JSONFile" begin
    include("jsonfile.jl")
  end

  @testset "TXTFile" begin
    include("txtfile.jl")
  end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e910b348-9e5a-40d9-be9a-e78aacc47743
using Pkg; Pkg.activate("."); Pkg.instantiate()

# ╔═╡ feb48b54-7423-11ec-2426-23da81f3c079
using SatelliteToolbox, Dates

# ╔═╡ c3e74733-4f49-4849-8d64-5def3f5392e9
using LinearAlgebra,Roots

# ╔═╡ 1bd4c734-e245-4027-ac70-2d8043922cbe
using Plots, Plots.PlotMeasures

# ╔═╡ 274cb358-741c-4b43-bd6a-f96f3453c7fa
md"Download some Earth orientation parameters:"

# ╔═╡ 5c813a37-6916-4214-bb0a-380ca33a767f
eop_IAU1980 = get_iers_eop();

# ╔═╡ 5568e66d-8d62-4b78-825a-ec9e382752ea
md"### Some Helpful Definitions"

# ╔═╡ a8e6acad-212e-4cbe-b2b9-03b7292636f0
"""
Speed of light in a vacuum
"""
c = 299792458 # m/s

# ╔═╡ cedb1362-47a2-44ed-abd8-34b37c9d009e
"""
The frequency of the LoRa signal used in the experiment was between 430 and 440 MHz
"""
fLoRa = 435e6 # MHz

# ╔═╡ e271d16d-d600-41c5-9d58-3b56d1683c04
"""
Approximate latitude, longitue, and altitude of the Dwingeloo Radio Telescope used to transmit the signal
"""
dwingelooLLA = (52.8121, 6.3971, 0) # 0 meters elevation

# ╔═╡ 9b2af5c2-42ea-4b4f-a2e7-9b67a98da790
"""
    earthLLA_r(latdeg, lngdeg, altkm, jd)

J2000 position on Julian day `jd` of a point on the Earth's surface specified by a latitude, longitude, and altitude above the geoid in kilometers.
"""
function earthLLA_r(latdeg, lngdeg, altkm, jd)

	rITRFtoJ2K = r_ecef_to_eci(ITRF(), J2000(), jd, eop_IAU1980)
	itrfPos = geodetic_to_ecef(deg2rad(latdeg), deg2rad(lngdeg), altkm*1000)

	posJ2K = rITRFtoJ2K * itrfPos
end

# ╔═╡ b4fc572f-a9e8-4dac-b7b8-023a4229429d
"""
    earthLLA_rv(latdeg, lngdeg, altkm, jd)

Position and velocity of a point on the Earth's surface on the Julian day `jd` in the J2000 frame. Velocity is derived from a forward finite difference of position.
"""
function earthLLA_rv(latdeg, lngdeg, altkm, jd)
	dt = 1 # seconds
	pos1 = earthLLA_r(latdeg, lngdeg, altkm, jd)
	pos2 = earthLLA_r(latdeg, lngdeg, altkm, jd+dt/86400)

	pos1, (pos2-pos1)/dt
end

# ╔═╡ 11e63381-a19c-4202-96bb-3588c9d8d3a5
"""
    moon_rv(jd)

Position and velocity of the Moon at the Julian date `jd` in the J2000 frame (ignoring conversion between UTC time and Barycentric Dynamical Time). Velocity
is derived from a forward finite difference of position.
"""
function moon_rv(jd)
	R = r_eci_to_eci(MOD(), J2000(), jd)

	dt = 1 # seconds

	jd2 = jd + dt/86400

	# convert MOD positions to J2000
	# TODO: convert JD_UTC to JD_TBD
	pos =  R*moon_position_i(jd)
	pos2 = R*moon_position_i(jd2)

	pos, (pos2-pos)/dt
end

# ╔═╡ 57fc4f13-2218-43fe-8193-350f6b6db12e
"""

    moonLLA_r(latdeg, lngdeg, altkm, jd)

J2000 position on Julian day `jd` of a point on the Moon's surface specified by a latitude, longitude, and an altitude above the Moon's average radius in kilometers.

I've defined the Moon frame has +x pointing from Moon center to Earth, +z pointing in the Moon's orbital angular momentum direction, and +y completing the right-handed set.
"""
function moonLLA_r(latdeg, lngdeg, altkm, jd)
	rMoon = 1738.1e3 # meters

	# assumes moon rv is in J2000
	moonpos, moonvel = moon_rv(jd)

	xhat = -normalize(moonpos)
	zhat = normalize(cross(moonpos, moonvel))
	yhat = normalize(cross(zhat, xhat))

	rMoonToJ2000 = [xhat yhat zhat]

	# convert polar angles to a Cartesian position in the Moon frame
	posMoonFrame = (rMoon + altkm*1000) * [cosd(latdeg)*cosd(lngdeg),
										   cosd(latdeg)*sind(lngdeg),
										   sind(latdeg)]
	# conver from my Moon frame to J2000
	posJ2000 = rMoonToJ2000 * posMoonFrame + moonpos
end

# ╔═╡ 73702699-2464-4b22-abc9-422a2285bd4d
"""
    moonLLA_rv(latdeg, lngdeg, altkm, jd)

Position and velocity of a point on the Moons's surface on the Julian day `jd` in the J2000 frame ignoring Barycentric Dynamical Time conversion. Velocity is derived from a forward finite difference of position.
"""
function moonLLA_rv(latdeg, lngdeg, altkm, jd)
	dt = 1 # second
	pos1 = moonLLA_r(latdeg, lngdeg, altkm, jd)
	pos2 = moonLLA_r(latdeg, lngdeg, altkm, jd+dt/86400)

	pos1, (pos2-pos1)/dt
end

# ╔═╡ 461c3f7b-c2a0-4ae2-849b-026978f4dfe3
md"""
## The Main Event: `moonbounce(...)`

With all the position and velocity functions defined, `moonbounce` is able to actually calculate the times of reflection on the Moon and reception back on Earth. The arguments to `moonbounce` are:
- an Earth lat/lng/alt location
- a Moon lat/lng/alt location (relative to my Moon frame defined above)
- a transmission time (as a Julian date)

The functions above make it easy to calculate the positions of (and thus distances and light delays between) points on the Earth and Moon, but we have to be careful because both the Earth and Moon are in motion.

If an Earth point and Moon point are 1.2 light seconds away at a transmission time $t$, it's not correct to say the signal reaches the Moon at $t + (1.2 \text{s})$. Instead, we need to account for the amount that the Moon has moved during the propagation time of the signal! If the Moon is moving slightly away from the transmitter, the arrival time will be something more like $t + (1.20001 \text{s})$. Or if it's moving toward the Earth a bit, the arrival time could be something like $t + (1.19999 \text{s})$. The functions `earthToMoonDelay` and `moonToEarthDelay` below are passed to a quick zero-finding solver to account for this effect.
"""

# ╔═╡ 0e00242a-1023-480f-a187-c5e601712d59
function moonbounce(earthlat, earthlng, earthaltkm,  moonlat, moonlng, moonaltkm,  transmitjd; assumedLightDelaySeconds = 1.0)

	# Cartesian position and velocity of the transmit site at transmission time
	transmitPos, transmitVel = earthLLA_rv(earthlat, earthlng, earthaltkm, transmitjd)

	# calculates the difference between a guess `delaySeconds` and the actual light
	# delay between the transmission site at transmission time and the reflection
	# site `delaySeconds` into the future. Returns zero when `delaySeconds` matches
	# the signal delay between transmission at the transmission site and a
	# reflection at the reflection site `delaySeconds` later
	function earthToMoonDelay(delaySeconds)
		moonPos, moonVel = moonLLA_rv(moonlat, moonlng, moonaltkm, transmitjd + delaySeconds/86400)

		dist = norm(moonPos - transmitPos)

		dist/c - delaySeconds
	end

	# Earth to Moon reflector light delay in seconds
	e2mSeconds = find_zero(earthToMoonDelay, assumedLightDelaySeconds)

	# time of signal reflection by a reflector on the Moon
	# and position and velocity of the reflector at that time
	reflectionTime = transmitjd + e2mSeconds/86400
	reflectPos, reflectVel = moonLLA_rv(moonlat, moonlng, moonaltkm, reflectionTime)

	# like `earthToMoonDelay` above but for the Moon->Earth direction
	function moonToEarthDelay(delaySeconds)
		earthPos, earthVel = earthLLA_rv(earthlat, earthlng, earthaltkm, reflectionTime + delaySeconds/86400)

		dist = norm(earthPos - reflectPos)

		dist/c - delaySeconds
	end

	# Moon to Earth receiver light delay in seconds
	m2eSeconds = find_zero(moonToEarthDelay, assumedLightDelaySeconds)

	# time of reception back on Earth,
	# and the position and velocity of the receiver at that time
	receptionTime = reflectionTime + m2eSeconds/86400
	receptionPos, receptionVel = earthLLA_rv(earthlat, earthlng, earthaltkm, receptionTime)

	(transmission=(time=transmitjd, pos=transmitPos, vel=transmitVel),
	 reflection=(time=reflectionTime, pos=reflectPos,  vel=reflectVel),
	 reception= (time=receptionTime,  pos=receptionPos,vel=receptionVel))
end

# ╔═╡ f3651d01-81bf-4e3e-ab60-cd70541ce24d
"""
    function dopplerFactorBetween(a, b)
Calculates the Doppler effect frequency shift factor of a light wave transmitted from one point and received by another, both with instantaneous velocities.

Based on the equation from the Wikipedia page [Doppler effect](https://en.wikipedia.org/wiki/Doppler_effect#General).
"""
function dopplerFactorBetween(a, b)
	# vector between the two points
	rAtoB = b.pos-a.pos
	rHat = normalize(rAtoB)

	# source velocity relative to the receiver
	vS = dot(a.vel, rHat)

	# receiver velocity, relative to source
	vR = dot(b.vel, -rHat)  # -rHat because we want the direction from b to a

	factor = (c+vR)/(c-vS)
end

# ╔═╡ 5b4fd441-452e-4bf2-b606-b4724909cdd7
md"
## Generating Points on the Moon

I'd like to evenly cover the Moon with points that will act as reflector locations. A great way to evenly distribute a number of points around a sphere is to use, of all things, the Fibonacci spiral. I'm following the procedure described by [extremelearning.com.au](http://extremelearning.com.au/how-to-evenly-distribute-points-on-a-sphere-more-effectively-than-the-canonical-fibonacci-lattice/)

"

# ╔═╡ e019b7ef-3671-4c48-97a0-aed0c0857efc
n = 500

# ╔═╡ c6af0b6f-955f-4ab3-a105-4996aa1a5383
# generate some reflector latitudes and longitudes
refLats, refLngs = let
	ϕ = (1 + 5^0.5)/2
	i = 0:n-1

	lng = rad2deg.(2pi/ϕ * i)
	# colatitude - zero is +z (the north pole), pi is -z (south pole)
	colat = acosd.(1 .- 2/n .* (i.+0.5))

	colat, lng
end

# ╔═╡ 303254ab-7ab0-476e-892c-a2866a48db29
md"""
## Simulate a bounce

With the `moonbounce` function in place, all that needs to be done is turn each of the transmit-reflect-receive sets it produces into a delay-Doppler pair. I'm also returning a "strength" result for each bounce, which hopefully is an approximate indication of how strong a reflection each bounce is. Bouncing off a part of the Moon facing directly at the transmitter should give a strong reflection, while a glancing bounce should be harder to see in the final plot.

Also included here is the handy `angleBetweend` function that returns the angle between two 3D vectors in degrees.
"""

# ╔═╡ 80c79503-c342-4f31-9356-bae9a29aaf4f
"""
    angleBetweend(a,b)
Angle between two (unnormalized) vectors in degrees

Examples:
```jldoctest
julia> using LinearAlgebra

julia> angleBetweend([1,0,0], [0,1,0])
90.0

julia> angleBetweend([1,0,0], [1,1,0])
45.0
```
"""
function angleBetweend(a,b)
	na = normalize(a)
	nb = normalize(b)
	atand(norm(cross(na, nb)), dot(na,nb))
end;

# ╔═╡ e165d46b-af7f-47aa-8fcb-5999606bdb0f
function dopplerDelayPairsAt(jd)

	# generate a Doppler-delay pair a signal transmitted at time t bouncing off each
	# of the reflectors on the Moon (defined above)
	dopplerDelayPairs = map(refLats, refLngs) do moonLat, moonLng
		bounce = moonbounce(dwingelooLLA..., moonLat,moonLng,0, jd)

		moonPosAtReflection = moon_rv(bounce.reflection.time)[1]
		reflectionEarthElevation = 90-angleBetweend(bounce.reflection.pos-moonPosAtReflection, bounce.transmission.pos-bounce.reflection.pos)

		transmitterMoonElevation = 90-angleBetweend(bounce.reflection.pos-bounce.transmission.pos, bounce.transmission.pos)

		if(reflectionEarthElevation > 5 && # skip "reflections" that have a low (or
										   # negative!) elevation on the Moon
			transmitterMoonElevation>10)   # and also skip ones that have too low
										   # an elevation on the Earth

			shift = fLoRa *
					(dopplerFactorBetween(bounce.transmission, bounce.reflection) *
					 dopplerFactorBetween(bounce.reflection, bounce.reception)
					 - 1)

			# convert delay in Julian days to delay in seconds
			delay = (bounce.reception.time - bounce.transmission.time)*86400

			shift, delay, sind(reflectionEarthElevation)
		else
			nothing, nothing, nothing
		end
	end

	# filter out all the passes that had an elevation that was too low
	dopplers = filter(x->!isnothing(x), map(x->x[1], dopplerDelayPairs));
	delays = filter(x->!isnothing(x), map(x->x[2], dopplerDelayPairs));
	strengths = filter(x->!isnothing(x), map(x->x[3], dopplerDelayPairs));

	dopplers, delays, strengths
end

# ╔═╡ 21e56be6-dafc-46e6-9831-b258d2d7aedb
# generate a delay-Doppler plot for a transmission at the given Julian date `jd`
function dopplerDelayPlotAt(jd; plotargs...)
	dopplers, delays, strengths = dopplerDelayPairsAt(jd)

	title = Dates.format(DateTime(jd_to_date(jd)[1:end-1]...), "u. dd HH:MM")

	scatter(dopplers, delays,
		alpha=strengths,     # first cut at Lambertian-type reflection strength TODO
		label=false,
		xlabel="Frequency shift (Hz)",
		ylabel="Delay since start of transmission (s)",
		markerstrokewidth=0,
		guidefont = 8, tickfont = 8,
		title = title;
		plotargs...
	)
end

# ╔═╡ ece89996-7e66-4923-81d2-549bd332752e
let
	transmitTimes = (DateTime(2021, 10, 5, 3, 24, 55):Dates.Hour(1):DateTime(2021, 10, 5, 18, 24, 55))

	bouncePlots = dopplerDelayPlotAt.(date_to_jd.(transmitTimes), aspect_ratio=400, ylabel="Delay (s)", xrot=40)

	p = plot(bouncePlots..., layout=(4,4),size=(1400,1000), bottom_margin=20px, left_margin=20px)

	titleplot = plot(title="Moon Bounce Doppler Delay October 5, 2021", grid=false, showaxis=false, xaxis=nothing, yaxis=nothing, bottom_margin=-20Plots.px)

	plot(titleplot, p, layout=@layout([A{0.01h}; B]))
end

# ╔═╡ 40f45a40-01d2-4569-8825-63a81779cd99
md"
## A Plot that Matches the Original

By guessing and checking, I found that a transimission time of 12:24:55 UTC on October 5, 2021 produces a plot that is a pretty good match to the one shared by Lacuna [here](https://lacuna.space/lora-moon-bounce/)
"

# ╔═╡ 59550d8b-13de-4d99-9164-b5b0a42fea7e
t = date_to_jd(DateTime(2021, 10, 5, 12, 24, 55))

# ╔═╡ 833c3151-9e37-4f72-9ffa-bf1958d9adc7
bounce = moonbounce(dwingelooLLA..., 0,0,0, t)

# ╔═╡ 7059f0f2-c620-4353-a2c7-2070d90320fe
dopplerFactorBetween(bounce.transmission, bounce.reflection)

# ╔═╡ 494659ca-1258-4f04-a8b3-c7889610ab43
dopplerFactorBetween(bounce.reflection, bounce.reception)

# ╔═╡ 356a3352-4a77-41e2-9206-749ca8d40983
# take two Doppler factors and convert to a frequncy shift in Hertz
shift = fLoRa * (dopplerFactorBetween(bounce.transmission, bounce.reflection) *dopplerFactorBetween(bounce.reflection, bounce.reception) - 1)

# ╔═╡ ce1238ea-0897-4d2b-902a-14da86387833
delay = bounce.reception.time - bounce.transmission.time

# ╔═╡ fd28654f-18aa-4e27-a83d-91f2e6cd81c4
dopplerDelayPlotAt(t)

# ╔═╡ Cell order:
# ╠═e910b348-9e5a-40d9-be9a-e78aacc47743
# ╠═feb48b54-7423-11ec-2426-23da81f3c079
# ╠═c3e74733-4f49-4849-8d64-5def3f5392e9
# ╠═1bd4c734-e245-4027-ac70-2d8043922cbe
# ╟─274cb358-741c-4b43-bd6a-f96f3453c7fa
# ╠═5c813a37-6916-4214-bb0a-380ca33a767f
# ╟─5568e66d-8d62-4b78-825a-ec9e382752ea
# ╠═a8e6acad-212e-4cbe-b2b9-03b7292636f0
# ╠═cedb1362-47a2-44ed-abd8-34b37c9d009e
# ╠═e271d16d-d600-41c5-9d58-3b56d1683c04
# ╠═9b2af5c2-42ea-4b4f-a2e7-9b67a98da790
# ╠═b4fc572f-a9e8-4dac-b7b8-023a4229429d
# ╠═11e63381-a19c-4202-96bb-3588c9d8d3a5
# ╠═57fc4f13-2218-43fe-8193-350f6b6db12e
# ╠═73702699-2464-4b22-abc9-422a2285bd4d
# ╟─461c3f7b-c2a0-4ae2-849b-026978f4dfe3
# ╠═0e00242a-1023-480f-a187-c5e601712d59
# ╠═f3651d01-81bf-4e3e-ab60-cd70541ce24d
# ╠═833c3151-9e37-4f72-9ffa-bf1958d9adc7
# ╠═7059f0f2-c620-4353-a2c7-2070d90320fe
# ╠═494659ca-1258-4f04-a8b3-c7889610ab43
# ╠═356a3352-4a77-41e2-9206-749ca8d40983
# ╠═ce1238ea-0897-4d2b-902a-14da86387833
# ╟─5b4fd441-452e-4bf2-b606-b4724909cdd7
# ╠═e019b7ef-3671-4c48-97a0-aed0c0857efc
# ╠═c6af0b6f-955f-4ab3-a105-4996aa1a5383
# ╟─303254ab-7ab0-476e-892c-a2866a48db29
# ╠═80c79503-c342-4f31-9356-bae9a29aaf4f
# ╠═e165d46b-af7f-47aa-8fcb-5999606bdb0f
# ╠═21e56be6-dafc-46e6-9831-b258d2d7aedb
# ╠═ece89996-7e66-4923-81d2-549bd332752e
# ╟─40f45a40-01d2-4569-8825-63a81779cd99
# ╠═59550d8b-13de-4d99-9164-b5b0a42fea7e
# ╠═fd28654f-18aa-4e27-a83d-91f2e6cd81c4

max_stars_repo_name: ThatcherC/moonbounce

max_issues_repo_name: ThatcherC/moonbounce

max_forks_repo_name: ThatcherC/moonbounce


text:
export deck, print_deck

"""
    deck()

Returns a 52-long list of all cards in random order.
Use `deck(false)` to get a list in new-box order (from Ace of Clubs
to King of Spades).
"""
function deck(shuffled::Bool = true)::Vector{Card}
    idx = collect(1:52)
    if shuffled
        idx = randperm(52)
    end

    result = Card.(idx)
end

"""
    print_deck(d)

Here `d` is a 52-long list of `Card`s
(that is, a full deck). Print `d` in four lines
with 13 cards per line.
"""
function print_deck(d::Vector{Card})
    @assert length(d) == 52 "You're not playing with a full deck"
    for j = 1:4
        a = 13 * (j - 1) + 1
        for k = a:a+12
            print(d[k])
            if k < a + 12
                print(" ")
            else
                println()
            end
        end
    end
end


text:
include("orbit/orbit.jl")
#include("external_model/external_model.jl")
#include("dynamic_model/dynamic_model.jl")
#include("dynamics/dynamics.jl")
#include("satellite/satellite.jl")
#include("plot/plot_plots.jl")
#include("plot/plot_makie.jl")
include("coordinates.jl")
include("satellite/target_decision.jl")

using SatelliteToolbox
using Dates
using LinearAlgebra
using Plots
gr()

using HDF5

function main()
	#=
	設定パラメータ
	=#
	DataNum =1500 #シミュレータ反復回数
	dt = 30 ##シミュレータの計算間隔 [s]
	start_time = DateTime(2020, 04, 09, 15, 18, 00)	#シミュレート開始時刻
	TLEFileName = "./orbit/ISS_TLE.txt"

	# 進行方向とSCOFx軸とのずれ（内積）
	dotvs = zeros(DataNum)

	# 撮影地点
	target_geod = [10. 10. 0.]
	limit_time = 1500

	# カメラパラメータ
	cam_viewangle = 40
	sat_axisval = 40

	#=
	# 軌道計算については先に行い、全  時間分を配列に保存する
	JD_log, x_ecef_log, v_ecef_log, x_geod_log = orbit_cal(DataNum,dt,start_time,TLEFileName)

	h5open("orbit.h5", "w") do file
		write(file, "JD_log", JD_log)  # alternatively, say "@write file A"
		write(file, "x_ecef_log", x_ecef_log)
		write(file, "v_ecef_log", v_ecef_log)
		write(file, "x_geod_log", x_geod_log)
	end
	=#
	JD_log = h5read("orbit.h5", "JD_log")
	x_ecef_log = h5read("orbit.h5", "x_ecef_log")
	v_ecef_log = h5read("orbit.h5", "v_ecef_log")
	x_geod_log = h5read("orbit.h5", "x_geod_log")

	println(size(x_ecef_log))
	target_geod = [deg2rad(target_geod[1]), deg2rad(target_geod[2]), target_geod[3]]
	target_ecef = GeodetictoECEF(target_geod[1], target_geod[2], target_geod[3])
	print("target_geod : ")
	println(target_geod)
	print("target_ecef : ")
	println(target_ecef)

	currenttime = 1
	shootingtime, shooting_vec = shootingtime_decision2(x_ecef_log[currenttime:currenttime+limit_time-1, :], target_ecef, currenttime, limit_time, cam_viewangle, sat_axisval)
	println(shootingtime)
	println(shooting_vec)

	M = ecef_to_DCM(x_ecef_log[shootingtime,:],v_ecef_log[shootingtime,:])
	shooting_vec_seof = M*shooting_vec

	current_attqua = SatelliteToolbox.Quaternion(1., 0., 0., 0.)
	targetqua, rotqua = targetqua_dicision(shooting_vec_seof, current_attqua, sat_axisval)
	println(targetqua)
	println(rotqua)

	# 描画のための処理
	target_geod_plot = [target_geod[1] target_geod[2]; target_geod[1] target_geod[2]]
	shoot_geod_plot = [x_geod_log[shootingtime,1] x_geod_log[shootingtime,2]; x_geod_log[shootingtime,1] x_geod_log[shootingtime,2]]

	shootcentor = x_ecef_log[shootingtime, :] + shooting_vec * x_geod_log[shootingtime, 3] / shooting_vec_seof[3]
	shootcentor_geod = ECEFtoGeodetic(shootcentor)
	cθ = dot(shooting_vec, -x_ecef_log[shootingtime, :])/(norm(shooting_vec)*norm(x_ecef_log[shootingtime, :]))
	tθ = sqrt(1/(cθ^2)-1)
	ta = tand(cam_viewangle/2)
	tθa = (tθ+ta)/(1-tθ*ta)
	shootrange = tθa * x_geod_log[shootingtime, 3] - tθ * x_geod_log[shootingtime, 3]
	r = atan(shootrange / norm(shootcentor))
	println(r)
	shootedge = zeros(37,2)
	for i = 1:37
		shootedge[i, :] = [shootcentor_geod[2] + r*cosd(i*10), shootcentor_geod[1] + r*sind(i*10)]
	end

	#=
	cθ = dot(shooting_vec, -x_ecef_log[shootingtime, :])/(norm(shooting_vec)*norm(x_ecef_log[shootingtime, :]))
	tθ = sqrt(1/(cθ^2)-1)
	ta = tand(cam_viewangle/2)
	tθa = (tθ+ta)/(1-tθ*ta)
	shootcentor_xdis = shooting_vec[1]/shooting_vec[3]*x_geod_log[shootingtime, 3]
	shootcentor_ydis = shooting_vec[2]/shooting_vec[3]*x_geod_log[shootingtime, 3]
	shootcentor = ECEFtoGeodetic(x_ecef_log[shootingtime, :] + [shootcentor_xdis, shootcentor_ydis, 0.])
	shootrange = tθa*x_geod_log[shootingtime, 3]
	shootendpoint = ECEFtoGeodetic(x_ecef_log[shootingtime, :] - [shootrange, 0., x_ecef_log[shootingtime, 3]])
	Δx = shootendpoint[1] - shootcentor[1]
	Δy = shootendpoint[2] - shootcentor[2]
	r = sqrt(Δx^2 + Δy^2)
	shootedge = zeros(37,2)
	for i = 1:37
		shootedge[i, :] = [shootcentor[2] + r*cosd(i*10), shootcentor[1] + r*sind(i*10)]
	end
	=#

	# 描画
	plot(x_geod_log[:,2], x_geod_log[:, 1], legend=:topleft)
	plot!(target_geod_plot[:,2], target_geod_plot[:,1], seriestype=:scatter, legend=:topleft)
	plot!(shoot_geod_plot[:, 2], shoot_geod_plot[:,1], seriestype=:scatter, legend=:topleft)
	plot!(shootedge[:, 1], shootedge[:, 2], legend=:topleft)
	xlabel!("x")
	ylabel!("y")
	#zlabel!("z")
	savefig("./figs/TargetDecisionTest2.png")

	plot(x_geod_log[shootingtime - 10:shootingtime + 10, 2],x_geod_log[shootingtime-10:shootingtime+10, 1], legend=:topleft)
	plot!(target_geod_plot[:, 2], target_geod_plot[:,1], seriestype=:scatter, legend=:topleft)
	plot!(shoot_geod_plot[:, 2], shoot_geod_plot[:,1], seriestype=:scatter, legend=:topleft)
	plot!(shootedge[:, 1], shootedge[:, 2], legend=:topleft)
	savefig("./figs/TargetDecisionTest3.png")

end

main()

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
function extrapolate(data::Array{T,2}, basis::Array{String,1}, methods::Array{String,1}) where T <: AbstractFloat
    X = []
    for b in basis
        if occursin("dz", lowercase(b))
            push!(X,2)
        elseif occursin("tz", lowercase(b))
            push!(X,3)
        elseif occursin("qz", lowercase(b))
            push!(X,4)
        else
            m = match(r"(\d)z", lowercase(b))
            if m != nothing
                push!(X,parse(Int64,m.captures[1]))
            else
                error("Cardinality of $b could not be determined")
            end
        end
    end
    push!(X, Inf)

    # Create new data array
    nb = length(basis)
    nm = length(methods)
    extrapolated = zeros(nb+1, nm)
    extrapolated[1:nb, :] .= data
    extrapolated_indexes = Tuple{Int64,Int64}[]

    # Extrapolate SCF energy
    extrapolated[nb+1, 1], _, _ = feller_formula(data[nb-2:end,1], X[nb-2:end])
    push!(extrapolated_indexes, (nb+1,1))

    # Loop through correlated energies

    for m = 2:length(methods)

        corr = filter(x->x!=0, data[:, m])
        if length(corr) > 2
            # Subtract HF energy
            corr = corr .- data[1:length(corr), 1]
            A,B = helgaker_formula(corr, X[1:length(corr)])

            for b = eachindex(X)
                if extrapolated[b,m] == 0.0
                    extrapolated[b,m] = A + B*X[b]^-3 + extrapolated[b,1]
                    push!(extrapolated_indexes,(b,m))
                end
            end
        else
            ΔE = extrapolated[1,m] - extrapolated[1,m-1]
            for b = 2:length(X)
                extrapolated[b,m] = extrapolated[b,m-1] + ΔE
                push!(extrapolated_indexes,(b,m))
            end
        end
    end

    return extrapolated, extrapolated_indexes
end

function feller_formula(E::Array{T,1}, X) where T <: AbstractFloat

    # SCF extrapolation: Escf = A + B*exp(-cX)
    # From: Feller 1993 (https://doi.org/10.1063/1.464749)

    # Get E1, E2, and E3
    l = length(E)
    E1 = E[l-2]
    E2 = E[l-1]
    E3 = E[l]

    X1 = X[l-2]
    X2 = X[l-1]
    X3 = X[l]

    C = log((E2 - E1)/(E3 - E2))
    B = (E3 - E2)/(exp(-C*X3) - exp(-C*X2))
    A = E3 - B*exp(-C*X3)

    return A,B,C
end

function helgaker_formula(E::Array{T,1}, X) where T <: AbstractFloat

    # Correlation energy extrapolation: Ecorr = A + B*X^-3
    # From Helgaker 1997 (https://doi.org/10.1063/1.473863)

    # Get E1, E2, and E3
    l = length(E)
    E1 = E[l-1]
    E2 = E[l]

    X1 = X[l-1]
    X2 = X[l]

    B = (E2 - E1)/(X2^-3 - X1^-3)
    A = E2 - B*X2^-3
    return A,B
end


text:
# common facilities

# tools to check size

function nmf_checksize(X, W::AbstractMatrix, H::AbstractMatrix)

    p = size(X, 1)
    n = size(X, 2)
    k = size(W, 2)

    if !(size(W,1) == p && size(H) == (k, n))
        throw(DimensionMismatch("Dimensions of X, W, and H are inconsistent."))
    end

    return (p, n, k)
end

# the result type

struct Result{T}
    W::Matrix{T}
    H::Matrix{T}
    niters::Int
    converged::Bool
    objvalue::T
    nmfinit::AbstractString
    nmfalg::AbstractString

    function Result{T}(W::Matrix{T}, H::Matrix{T}, niters::Int, converged::Bool, objv, nmfinit::AbstractString, nmfalg::AbstractString) where T
        if size(W, 2) != size(H, 1)
            throw(DimensionMismatch("Inner dimensions of W and H mismatch."))
        end
        new{T}(W, H, niters, converged, objv, nmfinit, nmfalg)
    end
end

# common algorithmic skeleton for iterative updating methods

abstract type NMFUpdater{T} end

function nmf_skeleton!(updater::NMFUpdater{T},
                       X, W::Matrix{T}, H::Matrix{T},
                       maxiter::Int, verbose::Bool, tol, nmfinit::AbstractString, nmfalg::AbstractString) where T
    objv = convert(T, NaN)

    # init
    state = prepare_state(updater, X, W, H)
    preW = Matrix{T}(undef, size(W))
    preH = Matrix{T}(undef, size(H))
    if verbose
        objv = evaluate_objv(updater, state, X, W, H)
        @printf("%-5s     %-13s    %-13s    %-13s\n", "Iter", "objv", "objv.change", "(W & H).change")
        @printf("%5d    %13.6e\n", 0, objv)
    end

    # main loop
    converged = false
    t = 0
    while !converged && t < maxiter
        t += 1
        copyto!(preW, W)
        copyto!(preH, H)

        # update H
        update_wh!(updater, state, X, W, H)

        # determine convergence
        dev = max(maxad(preW, W), maxad(preH, H))
        if dev < tol
            converged = true
        end

        # display info
        if verbose
            preobjv = objv
            objv = evaluate_objv(updater, state, X, W, H)
            @printf("%5d    %13.6e    %13.6e    %13.6e\n",
                t, objv, objv - preobjv, dev)
        end
    end

    if !verbose
        objv = evaluate_objv(updater, state, X, W, H)
    end
    return Result{T}(W, H, t, converged, objv, nmfinit, nmfalg)
end


text:
using DataFrames, CSV, StatsModels, LinearAlgebra, ForwardDiff, ForwardDiff, Optim, Distributions
using NLopt
using Metida
using SnoopCompile
using LineSearches
using BenchmarkTools
path    = dirname(@__FILE__)
cd(path)
df0         = CSV.File(path*"/csv/df0.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame
df1         = CSV.File(path*"/csv/df1.csv"; types = [String, String, String, String, Float64, Float64]) |> DataFrame
ftdf         = CSV.File(path*"/csv/1fptime.csv"; types = [String, String, Float64, Float64]) |> DataFrame
ftdf2        = CSV.File(path*"/csv/1freparma.csv"; types = [String, String, Float64, Float64]) |> DataFrame
ftdf3        = CSV.File(path*"/csv/2f2rand.csv"; types =
[String,  Float64, Float64, String, String, String, String, String]) |> DataFrame

################################################################################
# Metida
################################################################################

lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
)
@benchmark Metida.fit!($lmm, hes = false) seconds = 15

#=
BenchmarkTools.Trial: 527 samples with 1 evaluation.
 Range (min … max):  14.958 ms … 177.005 ms  ┊ GC (min … max):  0.00% … 89.00%
 Time  (median):     22.181 ms               ┊ GC (median):     0.00%
 Time  (mean ± σ):   28.448 ms ±  25.442 ms  ┊ GC (mean ± σ):  18.12% ± 17.00%

  ▂▄▇█▄▁
  ██████▇██▅▆▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▄▁▁▄█▄▁▆▅▅▁▄▄▁▁▄ ▇
  15 ms         Histogram: log(frequency) by time       150 ms <

 Memory estimate: 55.01 MiB, allocs estimate: 209813.
=#

################################################################################
# MetidaNLopt
################################################################################

lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
)
@benchmark Metida.fit!($lmm, hes = false, solver = :nlopt) seconds = 15

#=
BenchmarkTools.Trial: 333 samples with 1 evaluation.
 Range (min … max):  19.021 ms … 198.344 ms  ┊ GC (min … max): 0.00% …  0.00%
 Time  (median):     42.596 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   45.055 ms ±  20.382 ms  ┊ GC (mean ± σ):  7.54% ± 12.83%

          ▁▃▄▆█▆▂
  ▄▄▅▄▄▄▅████████▆▄▂▂▁▁▁▃▂▁▂▁▁▁▃▂▃▂▃▃▁▂▂▁▂▁▁▁▁▂▁▁▂▁▁▁▁▁▁▁▁▁▂▂▂ ▃
  19 ms           Histogram: frequency by time          141 ms <

 Memory estimate: 17.83 MiB, allocs estimate: 125293.
=#
lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
random = Metida.VarEffect(Metida.@covstr(1 + time|factor), Metida.ARH),
)
@benchmark  Metida.fit!(lmm, solver = :nlopt, hes = false)
#=
#No @natch
julia> @benchmark  Metida.fit!(lmm, solver = :nlopt, hes = false)
BenchmarkTools.Trial: 58 samples with 1 evaluation.
 Range (min … max):  80.086 ms … 97.076 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     86.918 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   87.458 ms ±  3.918 ms  ┊ GC (mean ± σ):  2.02% ± 2.75%

     ▄               ▄  ▄▁▁▄██ ▁▁        ▄ ▁                ▁
  ▆▁▆█▁▁▁▁▆▆▁▆▆▆▁▆▁▁▁█▆▁██████▆██▆▆▆▁▁▆▁▁█▁█▁▁▁▁▆▆▆▆▆▆▆▁▁▁▁▆█ ▁
  80.1 ms         Histogram: frequency by time          95 ms <

 Memory estimate: 33.28 MiB, allocs estimate: 2137.
=#
#=
#@batcj
BenchmarkTools.Trial: 101 samples with 1 evaluation.
 Range (min … max):  43.896 ms … 66.160 ms  ┊ GC (min … max): 0.00% … 19.07%
 Time  (median):     47.827 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   49.519 ms ±  4.247 ms  ┊ GC (mean ± σ):  4.05% ±  7.55%

            █
  ▃▃▃▁▃▃▃▄████▇▆▃▃▄▁▃▃▁▃▃▁▃▁▁▁▃▁▄▁▄▃▃▄▃▃▁▁▃▁▃▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▃ ▃
  43.9 ms         Histogram: frequency by time        64.6 ms <

 Memory estimate: 33.42 MiB, allocs estimate: 4009.
=#

################################################################################
# MetidaCu
################################################################################

lmm = Metida.LMM(@formula(response ~1 + factor*time), ftdf;
random = Metida.VarEffect(Metida.@covstr(1 + time|subject&factor), Metida.CSH),
)
@benchmark Metida.fit!($lmm, hes = false, solver = :cuda) seconds = 15

#=

=#

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

################################################################################
# Metida
################################################################################

lmm = Metida.LMM(@formula(tumorsize ~ 1 + CancerStage), hdp;
random = Metida.VarEffect(Metida.@covstr(1|HID), Metida.DIAG),
)
@benchmark  Metida.fit!(lmm, hes = false)

#=
BenchmarkTools.Trial: 1 sample with 1 evaluation.
 Single result which took 10.214 s (1.64% GC) to evaluate,
 with a memory estimate of 3.64 GiB, over 103755 allocations.
=#

################################################################################
# MetidaNLopt
################################################################################

lmm = Metida.LMM(@formula(tumorsize ~ 1 + CancerStage), hdp;
random = Metida.VarEffect(Metida.@covstr(1|HID), Metida.DIAG),
)
@benchmark  Metida.fit!(lmm, solver = :nlopt, hes = false)

#=
BenchmarkTools.Trial: 8 samples with 1 evaluation.
 Range (min … max):  517.141 ms … 822.903 ms  ┊ GC (min … max):  0.00% … 13.60%
 Time  (median):     644.920 ms               ┊ GC (median):    11.81%
 Time  (mean ± σ):   672.774 ms ± 115.055 ms  ┊ GC (mean ± σ):   8.02% ±  6.12%

  ▁            █          ▁ ▁                           ▁▁    ▁
  █▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁█▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁██▁▁▁▁█ ▁
  517 ms           Histogram: frequency by time          823 ms <

 Memory estimate: 919.67 MiB, allocs estimate: 39482.
=#

################################################################################
# MetidaCu
################################################################################

lmm = Metida.LMM(@formula(tumorsize ~ 1 + CancerStage), hdp;
random = Metida.VarEffect(Metida.@covstr(1|HID), Metida.DIAG),
)
@benchmark  Metida.fit!($lmm, solver = :cuda, hes = false)

#=
BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max):  3.331 s …   3.453 s  ┊ GC (min … max): 1.11% … 1.22%
 Time  (median):     3.392 s              ┊ GC (median):    1.16%
 Time  (mean ± σ):   3.392 s ± 85.816 ms  ┊ GC (mean ± σ):  1.16% ± 0.08%

  █                                                       █
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  3.33 s         Histogram: frequency by time        3.45 s <

 Memory estimate: 1.10 GiB, allocs estimate: 7312401
=#