text:
abstract type AbstractForcing end

struct NoForcing <: AbstractForcing end

statsheader(a::NoForcing) = ""

stats(a::NoForcing,s::AbstractSimulation) = (),()

msg(a::NoForcing) = "\nForcing: No forcing\n"

add_forcing!(u,::NoForcing) = nothing

struct RfForcing{T<:AbstractFloat} <: AbstractForcing
    Tf::T
    α::T
    Kf::T
    maxDk::T
    avgK::Vector{T}
    Zf::Vector{T}
    Ef::Vector{T} # Velocity Field Spectrum
    Em::Vector{T} # Target Spectrum
    R::Vector{T} # Solution to ODE
    factor::Vector{T} # Factor to multiply velocity Field
    forcex::PaddedArray{T,3,2,false} # Final force
    forcey::PaddedArray{T,3,2,false} # Final force
    init::Bool # Tell if the initial condition spectra should be used instead of from data
    hp::Base.RefValue{T}
    vp::Base.RefValue{T}
end

statsheader(a::RfForcing) = "hp,vp,tp"

stats(a::RfForcing,s::AbstractSimulation) = ((a.hp[],a.vp[],a.hp[]+a.vp[]),(a.hp[],a.vp[],a.hp[]+a.vp[]))

msg(a::RfForcing) = "\nForcing:  Rf forcing\nTf: $(getTf(a))\nalphac: $(getalpha(a))\nKf: $(getKf(a))\n"

@inline getTf(f::RfForcing) = f.Tf

@inline getalpha(f::RfForcing) = f.α

@inline getKf(f::RfForcing) = f.Kf

@inline getmaxdk(f::RfForcing) = f.maxDk

@inline getavgk(f::RfForcing) = f.avgK

@inline getZf(f::RfForcing) = f.Zf

function add_forcing!(u,f::RfForcing)
    add_Rf_forcing!(u.c.x,f.forcex)
    add_Rf_forcing!(u.c.y,f.forcey)
end

function add_Rf_forcing!(ui,f)
    @mthreads for j in YRANGE
        @inbounds for i in XRANGE
            ui[i,j,1] += f[i,j,1]
        end
    end

    for k in (2,3,4,5,6,NZ-4,NZ-3,NZ-2,NZ-1,NZ)
        @inbounds ui[1,1,k] += f[1,1,k]
    end
end

function initialize!(f::RfForcing,s)
    if f.init
        calculate_u1u2_spectrum!(f.Em,s.u,1)
    end
    if isfile("R.$(s.iteration[])")
        copyto!(f.R, vec(readdlm("R.$(s.iteration[])",Float64)))
    end
    return nothing
end

function calculate_Em!(Em,kh)
    data = readdlm("targSpectrum.dat",skipstart=1)
    E = interpolate((data[:,1],),data[:,2],Gridded(Linear()))
    for i in eachindex(Em)
        @inbounds Em[i] = E(kh[i])
    end
    return Em
end

#############################################################3
struct AForcing{T<:AbstractFloat} <: AbstractForcing
    Tf::T
    α::T
    Kf::T
    Zf::Matrix{T}
    Em::Matrix{T} # Target Spectrum
    R::Matrix{T} # Solution to ODE
    factor::Matrix{T} # Factor to multiply velocity Field
    forcex::PaddedArray{T,3,2,false} # Final force
    forcey::PaddedArray{T,3,2,false} # Final force
    hp::Base.RefValue{T}
    vp::Base.RefValue{T}
end

statsheader(a::AForcing) = "hp,vp"

stats(a::AForcing,s::AbstractSimulation) = (a.hp[],a.vp[]),(a.hp[],a.vp[])

msg(a::AForcing) = "\nForcing:  Aforcing\nTf: $(a.Tf)\nalphac: $(a.α)\nKf: $(a.Kf)\n"

function add_forcing!(u,f::AForcing)
    add_Rf_forcing!(u.c.x,f.forcex)
    add_Rf_forcing!(u.c.y,f.forcey)
end

function initialize!(f::AForcing,s)
    calculate_abs2!(f.Em,s.u,1)
    if isfile("R.$(s.iteration[])")
        copyto!(f.R, vec(readdlm("R.$(s.iteration[])",Float64)))
    end
    return nothing
end
#################################################################################################3
struct NRfForcing{T<:AbstractFloat} <: AbstractForcing
    Tf::T
    α::T
    Kf::T
    maxDk::T
    avgK::Vector{T}
    Zf::Vector{T}
    Ef::Vector{T} # Velocity Field Spectrum
    Em::Vector{T} # Target Spectrum
    R::Vector{T} # Solution to ODE
    N::Vector{T}
    Nm1::Vector{T}
    factor::Vector{T} # Factor to multiply velocity Field
    forcex::PaddedArray{T,3,2,false} # Final force
    forcey::PaddedArray{T,3,2,false} # Final force
    init::Bool # Tell if the initial condition spectra should be used instead of from data
    hp::Base.RefValue{T}
    vp::Base.RefValue{T}
end

statsheader(a::NRfForcing) = "hp,vp,tp"

stats(a::NRfForcing,s::AbstractSimulation) = (a.hp[],a.vp[],a.hp[]+a.vp[]),(a.hp[],a.vp[],a.hp[]+a.vp[])

msg(a::NRfForcing) = "\nForcing:  NRf forcing\nTf: $(a.Tf)\nalphac: $(a.α)\nKf: $(a.Kf)\n"

function add_forcing!(u,f::NRfForcing)
    add_Rf_forcing!(u.c.x,f.forcex)
    add_Rf_forcing!(u.c.y,f.forcey)
end

function initialize!(f::NRfForcing,s)
    if f.init
        calculate_u1u2_spectrum!(f.Em,s.u,1)
        calculate_nn_spectrum!(f.N,s.u,s.rhs,1)
        copyto!(f.Nm1,f.N)
    end
    if isfile("R.$(s.iteration[])")
        copyto!(f.R, vec(readdlm("R.$(s.iteration[])",Float64)))
    end
    if isfile("N.$(s.iteration[])")
        copyto!(f.N, vec(readdlm("N.$(s.iteration[])",Float64)))
    end
    return nothing
end
####################################################################################
struct CRfForcing{T<:AbstractFloat} <: AbstractForcing
    Tf::T
    α::T
    Kf::T
    maxDk::T
    avgK::Vector{T}
    Zf::Vector{T}
    Ef::Vector{T} # Velocity Field Spectrum
    Em::Vector{T} # Target Spectrum
    R::Vector{T} # Solution to ODE
    factor::Vector{T} # Factor to multiply velocity Field
    forcex::PaddedArray{T,3,2,false} # Final force
    forcey::PaddedArray{T,3,2,false} # Final force
    init::Bool # Tell if the initial condition spectra should be used instead of from data
    hp::Base.RefValue{T}
    vp::Base.RefValue{T}
end

statsheader(a::CRfForcing) = "hp,vp,tp"

stats(a::CRfForcing,s::AbstractSimulation) = (a.hp[],a.vp[],a.hp[]+a.vp[]),(a.hp[],a.vp[],a.hp[]+a.vp[])

msg(a::CRfForcing) = "\nForcing:  CRf forcing\nTf: $(getTf(a))\nalphac: $(getalpha(a))\nKf: $(getKf(a))\n"

@inline getTf(f::CRfForcing) = f.Tf

@inline getalpha(f::CRfForcing) = f.α

@inline getKf(f::CRfForcing) = f.Kf

@inline getmaxdk(f::CRfForcing) = f.maxDk

@inline getavgk(f::CRfForcing) = f.avgK

@inline getZf(f::CRfForcing) = f.Zf

function add_forcing!(u,f::CRfForcing)
    add_Rf_forcing!(u.c.x,f.forcex)
    add_Rf_forcing!(u.c.y,f.forcey)
end

function initialize!(f::CRfForcing,s)
    if f.init
        calculate_u1u2_spectrum!(f.Em,s.u,1)
    end
    if isfile("R.$(s.iteration[])")
        copyto!(f.R, vec(readdlm("R.$(s.iteration[])",Float64)))
    end
    return nothing
end

include("cforcing.jl")

#######################################################################
struct CNRfForcing{T<:AbstractFloat} <: AbstractForcing
    Tf::T
    α::T
    Kf::T
    maxDk::T
    avgK::Vector{T}
    Zf::Vector{T}
    Ef::Vector{T} # Velocity Field Spectrum
    Em::Vector{T} # Target Spectrum
    R::Vector{T} # Solution to ODE
    N::Vector{T}
    Nm1::Vector{T}
    factor::Vector{T} # Factor to multiply velocity Field
    forcex::PaddedArray{T,3,2,false} # Final force
    forcey::PaddedArray{T,3,2,false} # Final force
    init::Bool # Tell if the initial condition spectra should be used instead of from data
    hp::Base.RefValue{T}
    vp::Base.RefValue{T}
end

statsheader(a::CNRfForcing) = "hp,vp,tp"

stats(a::CNRfForcing,s::AbstractSimulation) = (a.hp[],a.vp[],a.hp[]+a.vp[]),(a.hp[],a.vp[],a.hp[]+a.vp[])

msg(a::CNRfForcing) = "\nForcing:  CNRf forcing\nTf: $(a.Tf)\nalphac: $(a.α)\nKf: $(a.Kf)\n"

function add_forcing!(u,f::CNRfForcing)
    add_Rf_forcing!(u.c.x,f.forcex)
    add_Rf_forcing!(u.c.y,f.forcey)
end

function initialize!(f::CNRfForcing,s)
    if f.init
        calculate_u1u2_spectrum!(f.Em,s.u,1)
        calculate_nn_spectrum!(f.N,s.u,s.rhs,1)
        copyto!(f.Nm1,f.N)
    end
    if isfile("R.$(s.iteration[])")
        copyto!(f.R, vec(readdlm("R.$(s.iteration[])",Float64)))
    end
    if isfile("N.$(s.iteration[])")
        copyto!(f.N, vec(readdlm("N.$(s.iteration[])",Float64)))
    end
    return nothing
end

include("cnrforcing.jl")
#######################################################################
struct QfForcing{T<:AbstractFloat} <: AbstractForcing
    alveliusP::T
    alveliusC::T
    alveliusKf::T
    Kf::T
    maxDk::T
    avgK::Vector{T}
    Ef::Vector{T} # Velocity Field Spectrum
    R::Vector{T} # Solution to ODE
    numPtsInShell2D::Vector{Int} # Factor to multiply velocity Field
    forcex::PaddedArray{T,3,2,false} # Final force
    forcey::PaddedArray{T,3,2,false} # Final force
    hp::Base.RefValue{T}
    vp::Base.RefValue{T}
end

statsheader(a::QfForcing) = "hp,vp,tp"

stats(a::QfForcing,s::AbstractSimulation) = (a.hp[],a.vp[],a.hp[]+a.vp[]),(a.hp[],a.vp[],a.hp[]+a.vp[])

msg(a::QfForcing) = "\nForcing:  Qf forcing\nP: $(a.alveliusP)\nc: $(a.alveliusC)\nalveliusKf: $(a.alveliusKf)\nKf: $(a.Kf)\n"

function add_forcing!(u,f::QfForcing)
    add_Rf_forcing!(u.c.x,f.forcex)
    add_Rf_forcing!(u.c.y,f.forcey)
end

function initialize!(f::QfForcing,s)
    return nothing
end

@inline getKf(f::QfForcing) = f.Kf
@inline getmaxdk(f::QfForcing) = f.maxDk
@inline getavgk(f::QfForcing) = f.avgK

function build_Qf_power_spectrum!(R,avgK,alveliusP, alveliusC,alveliusKf,kf,filter_width,lfilter)

    nShells2d = length(avgK)
    integral=0.0
    R .= 0.0
    @inbounds for i=2:nShells2d
        kh = avgK[i]
        if kh <= kf
            term = exp(-((kh - alveliusKf)/alveliusC)^2)
            R[i] = alveliusP*term
            integral += term#*(kh - avgK[i-1])
        end
    end
    R ./= integral

    if lfilter
        Δ² = filter_width^2
        @inbounds for i=2:nShells2d
            kh = avgK[i]
            if kh <= kf
                G = Gaussfilter(Δ²,kh*kh)
                R[i] *= G*G
            end
        end
    end

    return R
end
include("qfforcing.jl")
#################################################################################################

function forcing_model(d::AbstractDict,nx::Integer,ny::Integer,nz::Integer,ncx::Integer,lx::Real)

    forcingtype = NoForcing()

    if haskey(d,:forcing)
        if d[:forcing] == "rfForcing"
            TF = parse(Float64,d[:TF])
            alphac = parse(Float64,d[:alphac])
            kf = parse(Float64,d[:kf])
            nShells2D, maxdk2D, numPtsInShell2D, kh = compute_shells2D(KX,KY,ncx,ny)
            Ef = zeros(nShells2D)
            Em = zeros(nShells2D)
            R = zeros(nShells2D)
            factor = zeros(length(kh))
            forcex = PaddedArray((nx,ny,nz))
            forcey = PaddedArray((nx,ny,nz))
            Zf = calculate_Zf(d,kf,kh)
            if !isfile("targSpectrum.dat")
                forcingtype = RfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,factor,forcex,forcey,true,Ref(0.0),Ref(0.0))
            else
                calculate_Em!(Em,kh)
                forcingtype = RfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,factor,forcex,forcey,false,Ref(0.0),Ref(0.0))
            end
        elseif d[:forcing] == "aForcing"
            TF = parse(Float64,d[:TF])
            alphac = parse(Float64,d[:alphac])
            kf = parse(Float64,d[:kf])
            Em = zeros(length(KX),ny)
            R = zeros(length(KX),ny)
            factor = zeros(length(KX),ny)
            forcex = PaddedArray((nx,ny,nz))
            forcey = PaddedArray((nx,ny,nz))
            Zf = calculate_Zf(d,kf,KX,KY)
            forcingtype = AForcing{Float64}(TF, alphac, kf, Zf,Em,R,factor,forcex,forcey,Ref(0.0),Ref(0.0))
        elseif d[:forcing] == "nRfForcing"
            TF = parse(Float64,d[:TF])
            alphac = parse(Float64,d[:alphac])
            kf = parse(Float64,d[:kf])
            nShells2D, maxdk2D, numPtsInShell2D, kh = compute_shells2D(KX,KY,ncx,ny)
            Ef = zeros(nShells2D)
            Em = zeros(nShells2D)
            R = zeros(nShells2D)
            N = zeros(nShells2D)
            Nm1 = zeros(nShells2D)
            factor = zeros(length(kh))
            forcex = PaddedArray((nx,ny,nz))
            forcey = PaddedArray((nx,ny,nz))
            Zf = calculate_Zf(d,kf,kh)
            if !isfile("targSpectrum.dat")
                forcingtype = NRfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,N,Nm1,factor,forcex,forcey,true,Ref(0.0),Ref(0.0))
            else
                calculate_Em!(Em,kh)
                forcingtype = NRfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,N,Nm1,factor,forcex,forcey,false,Ref(0.0),Ref(0.0))
            end
        elseif d[:forcing] == "cRfForcing"
            TF = parse(Float64,d[:TF])
            alphac = parse(Float64,d[:alphac])
            kf = parse(Float64,d[:kf])
            nShells2D, maxdk2D, numPtsInShell2D, kh = compute_shells2D(KX,KY,ncx,ny)
            Ef = zeros(nShells2D)
            Em = zeros(nShells2D)
            R = zeros(nShells2D)
            factor = zeros(length(kh))
            forcex = PaddedArray((nx,ny,nz))
            forcey = PaddedArray((nx,ny,nz))
            Zf = calculate_Zf(d,kf,kh)
            if !isfile("targSpectrum.dat")
                forcingtype = CRfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,factor,forcex,forcey,true,Ref(0.0),Ref(0.0))
            else
                calculate_Em!(Em,kh)
                forcingtype = CRfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,factor,forcex,forcey,false,Ref(0.0),Ref(0.0))
            end
        elseif d[:forcing] == "cNRfForcing"
            TF = parse(Float64,d[:TF])
            alphac = parse(Float64,d[:alphac])
            kf = parse(Float64,d[:kf])
            nShells2D, maxdk2D, numPtsInShell2D, kh = compute_shells2D(KX,KY,ncx,ny)
            Ef = zeros(nShells2D)
            Em = zeros(nShells2D)
            R = zeros(nShells2D)
            N = zeros(nShells2D)
            Nm1 = zeros(nShells2D)
            factor = zeros(length(kh))
            forcex = PaddedArray((nx,ny,nz))
            forcey = PaddedArray((nx,ny,nz))
            Zf = calculate_Zf(d,kf,kh)
            if !isfile("targSpectrum.dat")
                forcingtype = CNRfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,N,Nm1,factor,forcex,forcey,true,Ref(0.0),Ref(0.0))
            else
                calculate_Em!(Em,kh)
                forcingtype = CNRfForcing{Float64}(TF, alphac, kf, maxdk2D, kh,Zf,Ef,Em,R,N,Nm1,factor,forcex,forcey,false,Ref(0.0),Ref(0.0))
            end
        elseif d[:forcing] == "QfForcing"
            alveliusP = parse(Float64,d[:alveliusP])
            alveliusC = parse(Float64,d[:alveliusC])
            alveliusKf = parse(Float64,d[:alveliusKf])
            kf = parse(Float64,d[:kf])
            nShells2D, maxdk2D, numPtsInShell2D, kh = compute_shells2D(KX,KY,ncx,ny)
            Ef = zeros(nShells2D)
            R = zeros(nShells2D)
            if haskey(d,:lesModel)
                lfilter = d[:lesModel] != "noLesModel" ? true : false
                filter_width = haskey(d,:filterWidth) ? Float64(eval(Meta.parse(d[:filterWidth]))) : (lx*2π/nx)/Globals.cutoffr
            else
                lfilter = false
                filter_width=0.0
            end
            build_Qf_power_spectrum!(R,kh,alveliusP,alveliusC,alveliusKf,kf,filter_width,lfilter)
            forcex = PaddedArray((nx,ny,nz))
            forcey = PaddedArray((nx,ny,nz))
            forcingtype = QfForcing{Float64}(alveliusP,alveliusC,alveliusKf,kf, maxdk2D, kh,Ef,R,numPtsInShell2D,forcex,forcey,Ref(0.0),Ref(0.0))
        end
    end

    return forcingtype
end

include("forcing.jl")
include("aforcing.jl")
include("nrforcing.jl")


text:
#' ## Run this cell to load the graphics packages
#+ results="hidden"
using Plots
#gr()
#plotly()
pyplot()
#+

#' ## Adding a function parameter
#' In the last notebook, we saw an example of adding **parameters** to
#' functions. These are values that control the behavior of a function.
#' Let's look at that in some more detail.
#' Let's go back to our original version of the σ function:
#+ results="hidden"
sigma(x) = 1/(1+exp(-x))
#+
#' Instead of working with a single function, we can work with a whole class
#' (set) of functions that look similar but differ in the value of a
#' **parameter**. Let's make a new function that uses the previous $\sigma$
#' function, but also has a parameter, $w$. Mathematically, we could write

#' $$f_w(x) = f(x; w) = \sigma(w \, x).$$

#' (Here, $w$ and $x$ are multiplied in the argument of $\sigma$; we could
#' write $w \times x$ or $w * x$, or even $w \cdot x$, but usually the symbols
#' are not written.)

#' Mathematically speaking, we can think of $f_w$ as a different function for
#' each different value of the parameter $w$.

#' In Julia, we write this as follows:
#+ results="hidden"
f(x, w) = sigma(w * x)
#+
#' Note that Julia just treats parameters as additional *arguments* to the
#' function; the function `f` has two arguments, the value of `x` and the
#' value of `w` that we want to use.

#' We can now investigate the effect of $w$ interactively. To do so, we need
#' a way of writing in Julia "the function of one variable $x$ that we obtain
#' when we fix the value of $w$". We write this as an "anonymous function", as
#' we saw in the notebook on functions:
#'
#'     x -> f(x, w)
#'
#' We can read this as "the function that maps $x$ to the value of $f(x, w)$,
#' for a value of $w$ that was previously given".

#' Now we are ready to draw the function. For each plot, we *fix* a value of
#' the parameter $w$ and draw the resulting function as a function of $x$.
#' We can interactively modify the value of $w$, and plot the new function that
#' comes out:
#' ## Try manipulating `w` between -2 and 2
w = -1
plot(x->f(x, w), -5, 5, ylims=(0,1), label="sigmoid")
plot!(x->(x>0), -5,5, label="Square Wave")

#' #### Exercise 1
#'
#' Try writing your own function that takes a parameter. Start by copying and
#' executing
#'
square(x) = x^2
#'
#' Then use `square` to declare a new function `square_and_scale` that takes
#' two inputs, `a` and `x` such that
#'
#' $$\mathrm{square\_and\_scale}(x; a) := a \cdot x^2$$

#' ####
function square_and_scale(x, a)
    return a*x^2
end

#' We can write some tests as follows, using the approximate equality operator,
# `≈`, written as `\approx<TAB>`:
square_and_scale(1, 1) ≈ 1
square_and_scale(1, 2) ≈ 2
square_and_scale(2, 1) ≈ 4
square_and_scale(2, 2) ≈ 8

# #### Exercise 2

#' Once you have declared `square_and_scale`, uncomment the code below and see
#' how the parameter `a` scales the function `square` :

x = -10:10
a = 5.0 # Try manipulating a here

plot(x, square.(x), label="x^2")
plot!(x, square_and_scale.(x, a), ls=:dash, label="ax^2")

#' ## Fitting a function to data
#' As we saw in the previous notebook, what we would like to do is use the
#' fact that we now have a parameter in our function in order to do something
#' useful! Namely, we want to model data with it.

#' So suppose that we are given a single data point $(x_0, y_0) = (2, 0.8)$.
#' We can try to "fit" the function $f_w$ by adjusting the parameter $w$ until
#' the function passes through the data.
#'
#' **Game**: Change `w` until the graph of the function hits the data point.
#' Which value of $w$ does that correspond to?
x0, y0 = 2, 0.8

#' f(x, w) = sigma(w * x)
w = 0.7 # try mainipulatoing `w` between -2 and 2
plot(x->f(x, w), -5, 5, ylims=(0, 1), label="f")
scatter!([x0], [y0], label="data")
#' ## Quantifying how far we are from the goal: the *loss function*
#' We can see visually when the graph of the function passes through the data
#' point. But the goal is to be able to automate this process so that the
#' computer can do it unaided.
#'
#' So we will need a more precise way of deciding and quantifying
#' (i.e. measuring with a number) *how far away we are from the goal*;
#' here, the goal means hitting the data point with the function.

#' #### Exercise 3
#' Can you think of a way of measuring how far away our function is from
#' the data point?

#' ### Defining the loss function
#' We need to measure how far away the curve is from the data point when we
#' choose a particular value of $w$.
#' One way to do this is by finding the vertical distance $d$ from the curve
#' to the data point.
#'
#' Instead of just taking the distance, it is common to take the *square* of
#' the distance, $d^2$.
#'
#' Since we are taking the vertical distance, we need the distance at the given
#' value of $x_0$ where the data point lies. For a given value of the parameter
#' $w$, the height of the point on the curve with that value of $x_0$ is
#' $f(x_0, w)$.
#'
#' So we take
#'
#' $$d := y_0 - f(x_0, w)$$
#'
#' and
#'
#' $$d^2 = [y_0 - f(x_0, w)]^2.$$
#'
#' This is our measure of distance. It will change when $w$ changes -- in other
#' words, it is itself a *function of $w$*; we will denote this function by
#' $L(w)$, and call it the **loss function**:
#'
#' $$L(w) := [y_0 - f(x_0, w)]^2.$$
#'
#' So the goal is to find the value $w^*$ of $w$ where the loss function is *least*; in other words, we need to *minimize* the loss function!
#'
#' (Another name for a loss function is a *cost function*.)

#' #### Exercise 4
#'
#' (a) Define the loss function `L(w)` in Julia.
#+ results="hidden"
x0, y0 = 2, 0.8
L(w) = (y0 - f(x0,w))^2
#+

#' (b) Draw the data point and the function `x -> f(x, w)`.
#' Also draw a vertical line from the data point to the function `x -> f(x, w)`.
w = 0.5
pl = plot(x->f(x,w), -5, 5, ylims=(0,1), label="f_w", lw=3)
scatter!(pl, [x0], [y0], label="data")
plot!([x0, x0], [y0, f(x0, w)], label="loss", legend=(0.2,0.9))
#' (c) Make the plot interactive.
#'
#' (d) Add as the plot title the value of the loss function for the current value
#' of $w$.
title!("L(w) = $(round(L(w); sigdigits = 5))")
#'
#' (e) Use the slider to find the value $w^*$ of $w$ for which the loss function
#' reaches its minimum value. What is $w^*$? What is the value of the loss
#' function there, $L(w^*)$?
w = 0.7
pl = plot(x->f(x,w), -5, 5, ylims=(0,1), label="f_w", lw=3)
scatter!(pl, [x0], [y0], label="data")
plot!([x0, x0], [y0, f(x0, w)], label="loss", legend=(0.1,0.9))
title!("L(w) = $(round(L(w); sigdigits = 5))")

#' ## What does the loss function look like?
#' The loss function $L(w)$ tells us how far away the function $f_w$ is from
#' the data when the parameter value is $w$, represented visually as the
#' vertical line in the previous plot. When the data are fixed, this is a
#' function only of the parameter $w$. What does this function look like as
#' a function of $w$? Let's draw it!

#' #### Exercise 5
#' Draw the function $L(w)$ as a function of $w$.

#' #### Solution
pl = plot(L, -2, 2, xlabel="w", ylabel="L(w)", ylims=(0, 0.7), label="L", legend=(0.8,0.9))
#' This graph quantifies how far we are from the data point for a given value
#' of $w$.
#' What features can we see from the graph?
#'
#' Firstly, we see that $L(w)$ is always bigger than $0$, for any value of $w$.
#' This is because we want $L$ to be some kind of measure of *distance*, and
#' distances cannot be negative.
#'
#' Secondly, we see that there is a special value $w^*$ of $w$ where the
#' function $L$ reaches its minimum value. In this particular case, it actually
#' reaches all the way down to $0$!
#' This means that the original function $f$ (the one we manipulated above)
#' passes exactly through the data point $(x_0, y_0)$.

#' #### Exercise 6
#' Draw a zoomed-in version of the graph to find the place $w^*$ where the
#' function hits $0$ more precisely.
pl = plot(L, 0.5, 1, xlabel="w", ylabel="L(w)", ylims=(0, 0.01), label="L", legend=(0.8,0.9))

#' ### A different way of defining the loss function
#' **Why did we use such a complicated function $L$ with those squares inside?**
#' We could instead just have used the absolute distance, instead of the
#' distance squared, using the *absolute value* function, written mathematically
#' as $| \cdot |$, and in Julia as `abs`.

#' #### Exercise 7
#'
#' Define a new loss function, `L_abs`, using the absolute value, and see what
#' it looks like.

#' #### Solution
L_abs(w) = abs(y0 - f(x0,w))

pl = plot(L_abs, 0.5, 1, xlabel="w", ylabel="L_abs(w)", ylims=(0, 0.01), label="L", legend=(0.8,0.9))

#' Now we see why it was previously generally preferred to use squares: using
#' the absolute value gives a cost function that is *not smooth*. This makes it
#' difficult to use methods from calculus to find the minimum. Nonetheless,
#' using non-smooth functions is very common in machine learning nowadays.

ext: jl


text:
"""
    vel_copy!(vel::NamedTuple, u, v, w, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
Copy external velocities into `PlanktonModel`
"""
function vel_copy!(vel::NamedTuple, u, v, w, g::AbstractGrid{TX, TY, TZ}) where {TX, TY, TZ}
    copy_interior_u!(vel.u.data, u, g, TX())
    copy_interior_v!(vel.v.data, v, g, TY())
    copy_interior_w!(vel.w.data, w, g, TZ())

    fill_halo_vel!(vel, g)
end

@inline function copy_interior!(c, t, g::AbstractGrid)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_u!(c, t, g::AbstractGrid, ::Periodic)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_v!(c, t, g::AbstractGrid, ::Periodic)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_w!(c, t, g::AbstractGrid, ::Periodic)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end

@inline function copy_interior_u!(c, t, g::AbstractGrid, ::Bounded)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx+1, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_v!(c, t, g::AbstractGrid, ::Bounded)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny+1, g.Hz+1:g.Hz+g.Nz), t)
end
@inline function copy_interior_w!(c, t, g::AbstractGrid, ::Bounded)
    copyto!(view(c, g.Hx+1:g.Hx+g.Nx, g.Hy+1:g.Hy+g.Ny, g.Hz+1:g.Hz+g.Nz+1), t)
end


text: error("Expected exception while running include")

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using GZip

const ACTIVITIES = "activities"
const FINDER_JAR = joinpath(@__DIR__, "..", "artifacts", "ActivityFinder.jar") #TODO: make proper artifact

function import_strava(path::String, dest_dir::String)
    mkpath(dest_dir)
    if endswith(path, ".zip")
        import_strava_zip(path, dest_dir)
    else
        error("not yet implemented")
    end
end

function import_strava_zip(path::String, dest_dir::String)

    bn = basename(path)
    bn, _ = splitext(bn)

    tmp_dir = mktempdir()
    unzip_dir = joinpath(tmp_dir, bn)
    unzip(path, unzip_dir)
    activities_dir = joinpath(unzip_dir, ACTIVITIES)
    for file in readdir(activities_dir, join=true)
        process_file(file, dest_dir)
    end
end

function import_garmin(path::String, dest_dir::String)
    mkpath(dest_dir)
    if endswith(path, ".zip")
        import_garmin_zip(path, dest_dir)
    else
        error("not yet implemented")
    end
end

function import_garmin_zip(path::String, dest_dir::String)
    bn, _ = splitext(basename(path))
    tmp_dir_fit = mktempdir()
    tmp_dir_zip = joinpath(mktempdir(), bn)

    unzip(path, tmp_dir_zip)
    uploaded_dir = joinpath(tmp_dir_zip, "DI_CONNECT",
            "DI-Connect-Fitness-Uploaded-Files")

    for file in readdir(uploaded_dir, join = true)
        _, ext = splitext(file)
        if ext == ".zip"
            unzip(file, tmp_dir_fit)
        end
    end

    activity_files = find_activity_files(tmp_dir_fit)
    for file in activity_files
        process_file(file, dest_dir)
    end
end

function process_file(path::String, dir::String)
    bn = basename(path)
    this_dir = dirname(path)
    nm, ext = splitext(bn)
    tcx_path = joinpath(dir, "$nm.tcx")
    if isfile(tcx_path)
        println("$tcx_path already exists, skipping.")
        return nothing
    end
    try
        if ext == ".tcx"
            mv(path, joinpath(dir, bn))
        elseif ext == ".fit"
            tcx_output = joinpath(this_dir, "$nm.tcx")
            fit_to_tcx(path, tcx_output)
            mv(tcx_output, tcx_path)
        elseif ext == ".gz"
            new_path = joinpath(this_dir, nm)

            GZip.open(path) do f
                contents = read(f)
                write(new_path, contents)
            end

            process_file(new_path, dir)
        else
            @warn "cannot process file $(basename(path)). Unknown file extension"
        end
    catch
        failed_dir = joinpath(dir, "failed")
        mkpath(failed_dir)
        failed_path = joinpath(failed_dir, bn)
        @warn "could not process file $bn. moving to $failed_path"
        mv(path, failed_path, force = true)
    end

    return nothing
end

function find_activity_files(dir::String)

    activity_files = joinpath(dir, "files.txt")

    cmd = Cmd(["java", "-jar", FINDER_JAR, dir])
    open(activity_files, "w") do file
        redirect_stdout(file) do
            run(cmd)
        end
    end

    return(readlines(activity_files))
end

function unzip(src::String, dest::String)
    if Sys.iswindows()
        cmd = Cmd(["PowerShell", "-Command", "Expand-Archive", "-LiteralPath", "'$src'", "-DestinationPath", "'$dest'"])
        run(cmd)
    else
        error("Not currently implemented for non-Windows machines.")
    end
end

function unzip(path::String)
    unzip(path, splitext(path)[1])
end

function fit_to_tcx(fit::String, tcx::String)
    cmd_strings = ["fittotcx", fit, tcx]
    if Sys.iswindows()
        cmd_strings = [["cmd", "/C"]; cmd_strings]
    end
    run(Cmd(cmd_strings))
end


text:
using Stheno: GPC

@timedtestset "indexing" begin
    @timedtestset "no noise" begin
        f, x = GP(EQ(), GPC()), randn(17)
        fx = f(x)
        @test mean(fx) == mean_vector(f, x)
        @test cov(fx) == cov(f, x)
    end
    @timedtestset "(Symmetric) Matrix-valued noise" begin
        rng, N = MersenneTwister(123456), 13
        f, x, A = GP(EQ(), GPC()), randn(rng, N), randn(rng, N, N)
        C = Symmetric(A * A' + I)
        fx = f(x, C)
        @test mean(fx) == mean_vector(f, x)
        @test cov(fx) == cov(f, x) + C
    end
    @timedtestset "Vector-valued noise" begin
        rng, N = MersenneTwister(123456), 11
        f, x, a = GP(EQ(), GPC()), randn(rng, N), exp.(randn(rng, N))
        fx = f(x, a)
        @test mean(fx) == mean_vector(f, x)
        @test cov(fx) == cov(f, x) + Diagonal(a)
    end
    @timedtestset "Real-valued noise" begin
        rng, N = MersenneTwister(123456), 13
        f, x, σ² = GP(EQ(), GPC()), randn(rng, N), exp(randn(rng))
        fx = f(x, σ²)
        @test mean(fx) == mean_vector(f, x)
        @test cov(fx) == cov(f, x) + Diagonal(fill(σ², N))
    end
end


text:
using FractionalTransforms
using Test

@testset "FractionalTransforms.jl" begin
    include("frft.jl")
    include("frst.jl")
    include("frct.jl")
end


text:
module Utils
using ..Types, ..scenario
using LinearAlgebra

export rotation, machNumber, temperature, density, gravity, T_MatrixQ, euler2quartenions, T_matrixE

function machNumber(v::Float64)
    ma = v / 340.0
    return ma
end
function machNumber(v::Float64, y::Float64)
    t = temperature(y)
    a = speedOfSound(t)
    ma = v / a
    return ma
end

function temperature(y::Float64;y0 = 0.0, t0 = 288.15, β = -6.5)
    t = (t0+β * (y-y0)*1e-3)#-273.15 #deg Celsius
    return t
end

function speedOfSound(t::Float64; γ = 1.4, R = 287.05287)
    a = sqrt(γ*R*t)
    return a
end
function pressure(y::Float64;p0 = 101325.0, y0 = 0.0, t0 = 288.15, β = -6.5, R = 287.05287)
    p = p0*(1+β/t0*(y-y0)*1e-3)^(-gravity(y)/(β*R))
    return p
end
function density(y::Float64; R = 287.05287)
    t = temperature(y)
    p = pressure(y)
    ρ = p/(R*t)
    return ρ
end

function gravity(y::Float64; g0 = 9.80665, r = 6356766.0)
    g = g0*(r/(r+y))^2
    return g
end

function Δrange(v0::Float64, θ::Float64, target::Target2D, trajectory::Function,proj::Projectile2D)
p= deepcopy(proj)
#println(proj.position)
#println(p.position)
vx = v0*cosd(θ)
vy = v0*sind(θ)
#p = createProjectile(0.03, 0.04,[vx,vy], [0.0,0.0])
p.velocity_m = [vx,vy]
trajectory(p,target)
#println(p.position)
Δrange = target.position[2]-p.position[2]

return Δrange
end

function bissection(target::Target2D, trajectory::Function, p::Projectile2D; θmin = 0.0, θmax = 45.0, ϵ = 0.001 )

#t2D = scenario.createTarget(0.2,[200.0,50.0])

#global θmin = 0.0
#global θmax = 45.0
#v0 = 900.0
#ϵ = 0.001
#global fmin = Δrange(inputs.v0, θmin,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
#global fmax = Δrange(inputs.v0, θmax,t2D, ExternalBallistics.PmmZeroDrag.timeOfFlight!)
#global error = min(abs(fmin),abs(fmax))
fmin = Δrange(inputs.v0, θmin,target, trajectory,p)
#println(θmin," ", θmax)
fmax = Δrange(inputs.v0, θmax,target, trajectory,p)
#println(fmin," ", fmax)
error = min(abs(fmin),abs(fmax))
θ = 0.0

if sign(fmin)==sign(fmax)
  println("Bisection method fails.")
else

  while error>ϵ

    #global θmin
    #global θmax
    #global fmin
    #global fmax
    #global error
    #global θ = 1/2*(θmin+θmax)
    θ = 1/2*(θmin+θmax)
    #println(error)

    f = Δrange(inputs.v0, θ,target, trajectory,p)
    error = abs(f)
    if f*fmin < 0
            θmax = θ
        else f*fmax < 0
            θmin = θ
      end
    end

  end
  return θ

end

function rotation(v::Array{Float64,1}, T::Array{Float64,2})
    vrot = T*v
    return vrot
end

function arctangent(A::Float64, B::Float64)
    if B>0
        return atan(A/B)
    elseif B==0 && A>0
        return pi/2
    elseif B==0 && A<0
        return -pi/2
    elseif B<0
        return pi+atan(A/B)
    end
end

function quaternions(λ::String, T)
    if λ == "λ0"
        λ0 = 1/2*(1+tr(T))^(1/2)
        λ1 = 1/λ0*1/4*(T[3,2]-T[2,3])
        λ2 = 1/λ0*1/4*(T[1,3]-T[3,1])
        λ3 = 1/λ0*1/4*(T[2,1]-T[1,2])
    elseif λ =="λ1"
        λ1 = 1/2*(1+2*T[1,1]-tr(T))^(1/2)
        λ0 = 1/λ1*1/4*(T[3,2]-T[2,3])
        λ2 = 1/λ1*1/4*(T[1,2]+T[2,1])
        λ3 = 1/λ1*1/4*(T[1,3]+T[3,1])
    elseif λ == "λ2"
        λ2 = 1/2*(1+2*T[2,2]-tr(T))^(1/2)
        λ0 = 1/λ2*1/4*(T[1,3]-T[3,1])
        λ1 = 1/λ2*1/4*(T[1,2]+T[2,1])
        λ3 = 1/λ2*1/4*(T[2,3]+T[3,2])
    elseif λ == "λ3"
        λ3 = 1/2*(1+2*T[3,3]-tr(T))^(1/2)
        λ0 = 1/λ3*1/4*(T[2,1]-T[1,2])
        λ1 = 1/λ3*1/4*(T[1,3]+T[3,1])
        λ2 = 1/λ3*1/4*(T[2,3]-T[3,2])
    end
    if sign(λ0)==-1
        λ0 = -1* λ0
        λ1 = -1* λ1
        λ2 = -1* λ2
        λ3 = -1* λ3
    end

    return λ0, λ1, λ2, λ3

end

function quaternions2euler(λ0::Float64, λ1::Float64, λ2::Float64, λ3::Float64)
θ = asin(2*(λ0*λ2-λ1*λ3))
ψ_num = 2*(λ1*λ2+λ0*λ3)
ψ_denom = (λ0^2+λ1^2-λ2^2-λ3^2)
ψ = arctangent(ψ_num,ψ_denom)
Φ_num = 2*(λ2*λ3+λ0*λ1)
Φ_denom = (λ0^2-λ1^2-λ2^2+λ3^2)
Φ = arctangent(Φ_num,Φ_denom)
return θ, ψ, Φ
end
function T_matrixE(θ::Float64, ψ::Float64, Φ::Float64)
    T = Array{Float64}(undef, 3,3)
    T[1,1] = cos(θ)*cos(ψ)
    T[1,2] = -cos(Φ)*sin(ψ) + sin(Φ)*sin(θ)*cos(ψ)
    T[1,3] = sin(Φ)*sin(ψ) + cos(Φ)*sin(θ)*cos(ψ)
    T[2,1] = cos(θ)*sin(ψ)
    T[2,2] = cos(Φ)*cos(ψ) + sin(Φ)*sin(θ)*sin(ψ)
    T[2,3] = -sin(Φ)*cos(ψ) + cos(Φ)*sin(θ)*sin(ψ)
    T[3,1] = -sin(θ)
    T[3,2] = sin(Φ)*cos(θ)
    T[3,3] = cos(Φ)*cos(θ)
    return T
end

function dominant(T)
    T00 = tr(T)
    trace = Dict("λ0"=>T00,"λ1"=>T[1,1],"λ2"=>T[2,2],"λ3"=>T[3,3])
    max_trace = maximum(values(trace))
    quater = "none"
    for (k,v) in trace
        if v == max_trace
            quater = k
        end
    end
    return quater
end

function euler2quartenions(θ::Float64, ψ::Float64, Φ::Float64)

    T = T_matrixE(θ, ψ, Φ)
    dominant_quater = dominant(T)

    λ0, λ1, λ2, λ3 = quaternions(dominant_quater, T)

    return λ0, λ1, λ2, λ3

end

function T_MatrixQ(λ0::Float64, λ1::Float64, λ2::Float64, λ3::Float64)
    T = Array{Float64}(undef, 3,3)
    T[1,1] = (λ0^2+λ1^2-λ2^2-λ3^2)
    T[1,2] = 2*(λ1*λ2-λ0*λ3)
    T[1,3] = 2*(λ1*λ3+λ0*λ2)
    T[2,1] = 2*(λ1*λ2+λ0*λ3)
    T[2,2] = (λ0^2-λ1^2+λ2^2-λ3^2)
    T[2,3] = 2*(λ2*λ3-λ0*λ1)
    T[3,1] = 2*(λ1*λ3-λ0*λ2)
    T[3,2] = 2*(λ2*λ3+λ0*λ1)
    T[3,3] = (λ0^2-λ1^2-λ2^2+λ3^2)

    return T

end

end


text:
module LightChem

include("molecule.jl")
# include("eval.jl")

end


text:
nheads = @parallel (+) for i=1:100000000
  Int(rand(Bool))
end

max_stars_repo_path: Chapter 1/2.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 1/2.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 1/2.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science


text:
module DFA

using ArrayViews
export dfa, integrate, polyfit

function integrate{T<:Real}(x::AbstractArray{T})
    N = length(x)
    y = zeros(N)
    sumx = 0.0
    for i=1:N
        sumx += x[i]
        y[i] = sumx
    end
    y
end

function integrate1{T<:Real}(x::AbstractArray{T})
    N = length(x)
    y = zeros(N)
    sumx = 0.0
    meanx = 0.0
    for i=1:N
        sumx += x[i]
        meanx += sumx/i
        y[i] = sumx - meanx
    end
    y
end

function polyfit{S<:Real, T<:Real}(x::AbstractArray{S},
                                   y::AbstractArray{T},
                                   order::Int = 1)
  A = [ float(x[i])^p for i = 1:length(x), p = 0:order ]
  A \ y
end

function dfa{T<:Real}(x::AbstractArray{T},
                      boxsize::Int;
                      order::Int = 1,
                      overlap::Real = 0.0,
                      integrated::Bool = false)
    if boxsize > 1 && order >= 0 && 0.0 <= overlap < 1.0
        if !integrated
            x = integrate(x)
        end
        overlapnum = ifloor(overlap*boxsize)
        N = length(x)
        tr = zeros(N)
        i = 0
        fluc = 0.0
        while i*(boxsize-overlapnum)+boxsize <= N
            boxstart = i*(boxsize - overlapnum) + 1
            boxend = i*(boxsize - overlapnum) + boxsize
            p = polyfit(boxstart:boxend, view(x, boxstart:boxend), order)
            for j=boxstart:boxend
                Δ = x[j] - sum([p[k+1]*j^k for k=0:order])
                fluc += Δ*Δ
            end
            i += 1
        end
        sqrt(fluc/i/boxsize)
    else
        error("boxsize at least 2, fit order at least 0, overlap must between [0, 1)")
    end
end

function dfa{T<:Real, S<:Real}(x::AbstractArray{T},
                               boxes::AbstractArray{S};
                               order::Int = 1,
                               overlap::Real = 0.0,
                               integrated::Bool = false)
    if !integrated
        x = integrate(x)
        integrated = true
    end
    N = length(boxes)
    fluc = zeros(N)
    for i=1:N
        fluc[i] = dfa(x, boxes[i], order=order, overlap=overlap, integrated=integrated)
    end
    fluc
end

function dfa{T<:Real}(x::AbstractArray{T};
                      order::Int = 1,
                      overlap::Real = 0.0,
                      boxmax::Int = div(length(x), 2),
                      boxmin::Int = 2*(order + 1),
                      boxratio::Real = 2,
                      integrated::Bool = false)
    if boxmax > boxmin
        boxes = unique(int(boxratio.^[log(boxratio, boxmin):log(boxratio, boxmax)]))
        if length(boxes) > 1
            if !integrated
                x = integrate(x)
                integrated = true
            end
            boxes, dfa(x, boxes, order=order, overlap=overlap, integrated=integrated)
        else
            error("box number at least 2")
        end
    else
        error("boxmax must be greater than boxmin, or length(x)/2 must be greater than 2(order+1)")
    end
end

end # module


text:
#=
This is an implementation of a priority_queue data structure, customized for modern cpu architectures

Special Thanks to Gandalf and other people in the Hoj discord server for helping me out.

The purpose of this implementation:
    priority_queue, despite being one crucial data structure in many applications
    (pathfinding, best-first-search, etc) has not been optimized for modern cpu architectures.
    Commonly used heaps (pairing heap, binary heap) has poor data locality.
    Moreover, this makes heavy use of inlining, making it fast.

This is a minimum viable version.

#Decrease_key is possible, but would require altering the elements/heaps to be aware of their locations.
There are several ways of doing this.
1) Making a handle class which contains the information and the data, then specialize the dispatch in case the heap stores that handle.
2) Making a pointer back to the parent heap.
3) Storing the information of the position somewhere.
4) Make a lookup table storing position for key positions.
2) is likely the most efficient option, but not gonna do it here.

There are also many ways to "defrag", such as by "transposing", merging, etc... but not doing it here either.
Even merging likely not implemented.
=#

struct summary_key{sub_heap_type,Dtype}
    key::sub_heap_type
    data::Dtype
    summary_key(a::sub_heap_type, b::Dtype) where{sub_heap_type, Dtype} = new{sub_heap_type,Dtype}(a,b)
end

#=
syms = [:<, :(==) , :>, :<=, :>=]
for i in syms
    @eval Base.$i(a::summary_key, b::summary_key) = $i(a.data, b.data)
end
No longer needed, the cmp function will now be generated during type compute.
=#
mutable struct Base_layer{Dtype, Base_size, cmp}
    data::Vector{Dtype}
    size::UInt8
    function Base_layer{Dtype, Base_size, cmp}() where {Dtype,Base_size, cmp}
        x = new{Dtype,Base_size, cmp}()
        x.data = Vector{Dtype}(undef,Base_size)
        x.size = 0
        x
    end
end

mutable struct Heap_layer{sub_layer, summary_layer}
    #TODO, fix summary to make it point to appropriate vector.
    summary::summary_layer
    #data::Vector{sub_layer} Fixing the issue.
    #size::UInt
    function Heap_layer{sub_layer, summary_layer}() where {sub_layer, summary_layer}
        x = new{sub_layer, summary_layer}()
        x.summary = summary_layer()
        #x.data = Vector{sub_layer}(undef, max_size(x.summary))
        #x.size = 0
        x
    end
end

function Heap_Heap_type(Dtype::Type, x::Integer, cmp::Function, Base_size::Integer = UInt8(32))
    #TODO: Fix summary key thing.
    if x==0
        return Base_layer{Dtype, Base_size, cmp}
    else
        @inline cmp_key(a,b) = cmp(a.data, b.data)
        sub_heap_type = Heap_Heap_type(Dtype, x-1, cmp, Base_size)
        return Heap_layer{sub_heap_type, Heap_Heap_type(summary_key{sub_heap_type,Dtype}, x-1, cmp_key, Base_size)}
    end
end

@inline function max_size(::Base_layer{Dtype, Base_size, cmp})::UInt where {Dtype, Base_size, cmp}
    return Base_size
end

@inline function max_size(A::Heap_layer{sub_layer, summary_layer})::UInt where {sub_layer, summary_layer}
    return max_size(A.summary)^2
end

@inline function need_splitting(X::Base_layer{Dtype, Base_size, cmp}) where {Dtype, Base_size, cmp}
    return X.size >= Base_size #Maybe replace this with == for improved performance? Considering that X.size only increases by one each time.
end
@inline function push!(X::Base_layer{Dtype, Base_size, cmp}, content::Dtype) where {Dtype, Base_size, cmp}
    #This does not need splitting.
    #Binary search the array to insert.
    #println("Working in progress.")
    #Wait, can I also unroll this too? #Later.

    if X.size == 0
        X.size = 1
        X.data[1] = content
        return
    elseif !cmp(content,X.data[X.size])
        X.size += 1
        X.data[X.size] = content
        return

    end

    mindex::UInt8 = 1
    maxdex::UInt8 = X.size
    while mindex < maxdex
        middex = (mindex+maxdex)>>1
        if cmp(content, X.data[middex])
            maxdex = middex
        else
            mindex = middex+1
        end
    end
    #May replace "end" with "X.size" (careful consideration needed)
    X.size += 1
    X.data[mindex+1:end] =  X.data[mindex:end-1]
    X.data[mindex] = content
    return
end

@inline function length(X::Base_layer)
    return X.size
end

@inline function empty(X::Base_layer)
    return length(X) == 0
end

@inline function pop!(X::Base_layer)
    data = X.data[1]
    #May replace "end" with "X.size"
    X.data[1:end-1] = X.data[2:end] #May need optimiza tion.I don't know how Julia arrays are implemented.
    X.size -= 1
    #X.data[end] = undef #Free the memory.
    return data
    #return the first element.
end

@inline function peek(X::Base_layer)
    return X.data[1]
end

@inline function replace!(X::Base_layer{Dtype, Base_size, cmp}, content::Dtype) where {Dtype, Base_size, cmp}
    #pop one element, then push one element in, like python's heapq.
    Out::Dtype = X.data[1]
    if !cmp(content,X.data[X.size])
        #In last element, needs to be taken care of separately.
        X.data[1:X.size-1] = X.data[2:X.size]
        X.data[X.size] = content
        return Out
    end

    mindex::UInt8 = 2 #The 1st element will be popped out regardless. The element this can come before is the 2nd element onward.
    #No need for special case of size 1: this will be tackled.
    maxdex::UInt8 = X.size
    #The 1st element shouldn't matter because it gets popped anyway?
    while mindex < maxdex
        middex = (mindex+maxdex)>>1
        if cmp(content, X.data[middex])
            maxdex = middex
        else
            mindex = middex+1
        end
    end
    X.data[1:mindex-2] = X.data[2:mindex-1]
    X.data[mindex-1] = content
    #What to do next?
    return Out
end

@inline function split!(X::Base_layer{Dtype, Base_size, cmp}) where {Dtype, Base_size, cmp}
    #create one new base layers.
    #Assumes full X.
    new_one = Base_layer{Dtype, Base_size, cmp}()
    cutoff = div(Base_size, 2)
    new_one.size = Base_size-cutoff
    new_one.data[1:new_one.size] = X.data[cutoff+1:end]
    X.size = cutoff
    return new_one
end

@inline function need_splitting(X::Heap_layer)
    return need_splitting(X.summary)
end

@inline function push!(X::Heap_layer{sub_layer, summary_layer}, content) where {sub_layer, summary_layer}
    #X.size += 1
    if empty(X.summary)
        new_one = sub_layer()
        push!(new_one, content)
        push!(X.summary, summary_key(new_one, content))
        return
        #println("Working in progress.")
    end

    summary_peek = peek(X.summary)
    sub_heap = summary_peek.key
    if need_splitting(sub_heap)
        #Do_something
        new_one = split!(sub_heap)
        push!(sub_heap, content)
        #Get new summary key.

        push!(X.summary, summary_key(new_one, peek(new_one)))
        #println("Working in progress")
    else
        push!(sub_heap, content)
    end
    return
end

@inline function length(X::Heap_layer)
    #Might not even store X, in which case, it is 0.
    #return X.size
end

@inline function empty(X::Heap_layer)
    return empty(X.summary)
end

@inline function pop!(X::Heap_layer)
    #Pop the 1st one, then the summary heap.
    #X.size -= 1
    summary_peek = peek(X.summary)
    sub_heap = summary_peek.key
    out = pop!(sub_heap)
    if empty(sub_heap)
        pop!(X.summary)
        #println("Working in progress")
        #Todo: make sure the data stays in place.
    else
        #println(X.summary)
        replace!(X.summary, summary_key(sub_heap, peek(sub_heap)))
    end
    return out

end

@inline function peek(X::Heap_layer)
    return peek(X.summary).data
end

@inline function replace!(X::Heap_layer, content)
    summary_peek = peek(X.summary)
    sub_heap = summary_peek.key
    out = replace!(sub_heap, content)
    replace!(X.summary, summary_key(sub_heap, peek(sub_heap)))
    return out
    # Replace the pointed heap then replace the summary.
end

@inline function split!(X::Heap_layer{sub_layer, summary_layer}) where {sub_layer, summary_layer}
    new_summary = split!(X.summary)
    new_heap = Heap_layer{sub_layer, summary_layer}()
    new_heap.summary = new_summary
    return new_heap
end

#Trust the compiler.
@inline function debug_print(X)
    print(X, " ")
end
@inline function debug_print(X::summary_key)
    #println("Summary key with key:")
    #print("key = ")
    #debug_print(X.key)
    print("data = ")
    debug_print(X.data)
end
@inline function debug_print(X::Base_layer)
    #println("Base layer:")
    #println(typeof(X))
    println("[")
    for i in 1:X.size
        debug_print(X.data[i])
        print(",")
    end
    println("]")
end

@inline function debug_print(X::Heap_layer)
    #println("Heap layer with summary:")
    print("summary = ")
    debug_print(X.summary)
end

function New_heap(Dtype::Type; layer::Integer=2, cmp::Function = <, base_size::Integer = 64)
    type = Heap_Heap_type(Dtype,layer, cmp, base_size)
    return type()
end


text:
using DocSeeker
import Markdown

handle("searchdocs") do data
  @destruct [mod || Main, nameOnly || false, exportedOnly || false, allPackages || false, query] = data
  items = @errs DocSeeker.searchdocs(query, mod = mod, exportedonly = exportedOnly,
                                     loaded = !allPackages, name_only = nameOnly)
  items isa EvalError ?
    Dict(:error => true, :errmsg => sprint(showerror, items.err)) :
    Dict(:error => false, :items => [renderitem(i[2]) for i in items], :scores => [i[1] for i in items])
end

function renderitem(x)
  r = Dict(f => getfield(x, f) for f in fieldnames(DocSeeker.DocObj))
  r[:html] = view(renderMD(x.html))
  r
end

handle("moduleinfo") do data
  @destruct [mod] = data
  d, items = getmoduleinfo(mod)
  items = [renderitem(i) for i in items]
  Dict(:doc => view(d), :items => items)
end

handle("regenerateCache") do
  DocSeeker.createdocsdb()
end

ispackage(mod) = Base.find_package(mod) ≠ nothing

function modulesymbols(mod)
  syms = filter(x -> x.mod == mod, DocSeeker.alldocs())
  sort(syms, by = x -> x.name)[1:min(100, length(syms))]
end

function packageinfo(mod)
  path = DocSeeker.readmepath(mod)
  readme = ispath(path) ? String(read(path)) : ""

  return  Hiccup.div(
            renderMD(Markdown.parse(readme)),
            Hiccup.Node(:hr),
            Hiccup.h2("defined symbols:")
          ), modulesymbols(mod)
end

function moduleinfo(mod)
  header = if mod == "Core"
    renderMD("## Julia `Core`")
  elseif first(split(mod, '.')) == "Base"
    renderMD("## Julia Standard Library: `$mod`")
  else
    renderMD("## Module `$mod`")
  end

  header, modulesymbols(mod)
end

getmoduleinfo(mod) = ispackage(mod) ? packageinfo(mod) : moduleinfo(mod)


text:
"""
    getknn(Xtrain, X; k = 1, metric = "eucl")
Return the k nearest neighbors in Xtrain of each row of `X`.
* `Xtrain` : Training X-data.
* `X` : Query X-dta.
* `metric` : Type of distance used for the query.
    Possible values are "eucl" or "mahal".

The distances (not the squared distances) are also returned.
"""
function getknn(Xtrain, X; k = 1, metric = "eucl")
    Xtrain = ensure_mat(Xtrain)
    Xt = ensure_mat(X')
    n, p = size(Xtrain)
    k > n ? k = n : nothing
    if metric == "eucl"
        ztree = BruteTree(Xtrain', Euclidean())
        ind, d = knn(ztree, Xt, k, true)
    elseif metric == "mahal"
        S = Statistics.cov(Xtrain, corrected = false)
        if p == 1
            Uinv = inv(sqrt(S))
        else
            #S = S + Diagonal(1e-10 * ones(p))
            Uinv = LinearAlgebra.inv!(cholesky!(Hermitian(S)).U)
        end
        zXtrain = Xtrain * Uinv
        zX = X * Uinv
        ztree = BruteTree(zXtrain', Euclidean())
        # ztree = BruteTree(Xtraint, Mahalanobis(Sinv))
        # is very slow
        ind, d = knn(ztree, zX', k, true)    # ind and d = lists
    end
    #ind = reduce(hcat, ind)'
    (ind = ind, d = d)
end


text:
using ITensors
using ITensorVisualization
using GeometryBasics
using Graphs
using NetworkLayout

N = 10
g = grid((N,))
tn = itensornetwork(g; linkspaces=10, sitespaces=2)
@visualize fig tn siteinds_direction = Point(1, -0.5) layout = SquareGrid(; cols=1) width =
  20 height = 50
fig


text:
const INSERTION_VALUE = 60

"""
    selectionsort!(a::AbstractVector{T})

A selection sort algorithm. Time complexity O(n^2).
"""
function selectionsort!(a::AbstractVector{T}) where T <: Number
    for i in 1:length(a)
        min = a[i]
        k = i
        for j in i+1:length(a)
            if a[j] < min
                min = a[j]
                k = j
            end
        end
        swap!(a, i, k)
    end

    return a
end

"""
    insertionsort!(a::AbstractVector{T}, lo::Int = 2, hi::Int = length(a))

Insertion sort with linear search. Time complexity O(n^2), best case O(n).
Fast for small arrays.
"""
function insertionsort!(a::AbstractVector{T}, lo::Int = 2, hi::Int = length(a)) where T
    i = lo
    while i <= hi
        j = i
        while j > lo - 1 && a[j - 1] > a[j]
            swap!(a, j, j - 1)
            j -= 1
        end
        i += 1
    end

    return a
end

"""
    mergesort(a::AbstractVector{T})

Mergesort algorithm with insertion sort at lowest level. Time complexity O(n log n).
"""
function mergesort(a::AbstractVector{T}) where T
    return mergesort(a, 1, length(a))
end

"""
    mergesort(a::AbstractVector{T}, lo::Int, hi::Int)

Recursive mergesort function. Something should be done about all the allocations.
"""
function mergesort(a::AbstractVector{T}, lo::Int, hi::Int) where T
    if hi - lo <= INSERTION_VALUE
        return insertionsort!(a, lo + 1, hi)[lo:hi]
    else
        mid = lo + (hi - lo) ÷ 2
        a1 = mergesort(a, lo, mid)
        a2 = mergesort(a, mid + 1, hi)

        a_ = Array{T}(undef, length(a1) + length(a2))
        merge!(a_, a1, a2)
        return a_
    end
end

"""
    merge(a::AbstractVector{T}, a1::AbstractVector{T}, a2::AbstractVector{T})

Merge two arrays.
"""
function merge!(a::AbstractVector{T}, a1::AbstractVector{T}, a2::AbstractVector{T}) where T
    i = 1; j = 1; k = 1
    while j <= length(a1) && k <= length(a2)
        if a1[j] < a2[k]
            a[i] = a1[j]
            j += 1
        else
            a[i] = a2[k]
            k += 1
        end

        i += 1
    end
    a[i:end] = [a1[j:end]; a2[k:end]]

    return nothing
end

"""
    swap(a::AbstractVector, i::Int, j::Int)

Swap to element i and j in array a.
"""
function swap!(a::AbstractVector, i::Int, j::Int)
    tmp = a[j]
    a[j] = a[i]
    a[i] = tmp
    return nothing
end

"""
    quicksort!(a::AbstractVector)

Basic quicksort algorithm. Time complexity O(nlogn), worst case with bad pivor O(n^2).
"""
quicksort!(a::AbstractVector) = quicksort!(a, 1, length(a))

"""
    quicksort!(a::AbstractVector, lo::Int, hi::Int)

Basic quicksort algorithm. Time complexity O(nlogn), worst case with bad pivor O(n^2).
"""
function quicksort!(a::AbstractVector, lo::Int, hi::Int)
    if lo >= hi return a end

    if hi - lo <= INSERTION_VALUE
        return insertionsort!(a, lo + 1, hi)
    else
        pivot = partition(a, lo, hi)
        quicksort!(a, lo, pivot - 1)
        quicksort!(a, pivot + 1, hi)
    end
end

"""
    partition(a, lo::Int, hi::Int)

Partition function for quicksort.
"""
function partition(a::AbstractVector, lo::Int, hi::Int)
    pivot = pivot_mo3(a, lo, hi) # median of three
    lo0 = lo
    lo += 1

    while true
        while lo <= hi && a[lo] < pivot
            lo += 1
        end
        while hi >= lo && a[hi] > pivot
            hi -= 1
        end
        if lo > hi
            break
        end

        swap!(a, lo, hi)
        lo += 1
        hi -= 1
    end
    swap!(a, lo0, hi)

    return hi
end

"""
    pivot_mo3(a::AbstractVector, lo::Int, hi::Int)

Use median of three elements as pivot element, i.e. middle element of three.
"""
function pivot_mo3(a::AbstractVector, lo::Int, hi::Int)
    mid = lo + (hi - lo) ÷ 2
    if a[mid] < a[lo]
        swap!(a, lo, mid)
    end
    if a[hi] < a[lo]
        swap!(a, lo, hi)
    end
    if a[mid] < a[hi]
        swap!(a, mid, hi)
    end
    swap!(a, lo, hi)
    pivot = a[lo]
end


text:
module cache

import SynapseClient: @cachefunction, @standalonefunction, @synapsefunction, @createtype, AbstractSynapse, synapsecall, synapseclient, pycache

using PyCall

export
	add,
	contains,
	get,
	get_cache_dir,
	purge,
	remove

@createtype(Cache,AbstractSynapse,pycache.Cache,PyObject,false)

@synapsefunction Cache add
@synapsefunction Cache contains
@synapsefunction Cache get
@synapsefunction Cache get_cache_dir
@synapsefunction Cache purge
@synapsefunction Cache remove

end


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/julia/1.6.0/bin/julia
#SBATCH --mail-user=mmcleod2@ualberta.ca
#SBATCH --mail-type=ALL
#SBATCH -o twod_hail.out # Standard output
#SBATCH -e twod_hail.err # Standard error
#SBATCH --time=02:30:00
#SBATCH --nodes=2
#SBATCH --ntasks-per-node=48
#SBATCH --mem=0
#SBATCH --account=def-amw8

using Pkg
Pkg.activate(".")

include(joinpath(ENV["SLURM_SUBMIT_DIR"], "parallel/parallel_config.jl"))
reproduce_config_experiment("configs/twodgridworld_hailmary/twodgridworld_gpi.toml", save_path="~/scratch/curiosity")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
    AbstractIndices{I} <: AbstractDictionary{I, I}

Abstract type for the unique keys of an `AbstractDictionary`. It is itself an `AbstractDictionary` for
which `getindex` is idempotent, such that `indices[i] = i`. (This is a generalization of
`Base.Slice`).

At minimum, an `AbstractIndices{I}` must implement:

 * The `iterate` protocol, returning unique values of type `I`.
 * `in`, such that `in(i, indices)` implies there is an element of `indices` which `isequal` to `i`.
 * Either `length`, or override `IteratorSize` to `SizeUnknown`.

While an `AbstractIndices` object is a dictionary, the value corresponding to each index is
fixed, so `issettable(::AbstractIndices) = false` and `setindex!` is never defined.

If arbitrary indices can be added or removed from the set, implement:

* `isinsertable(::AbstractIndices)` (returning `true`)
* `insert!(indices::AbstractIndices{I}, ::I}` (returning `indices`)
* `delete!(indices::AbstractIndices{I}, ::I}` (returning `indices`)
"""
abstract type AbstractIndices{I} <: AbstractDictionary{I, I}; end

@inline function Base.getindex(indices::AbstractIndices{I}, i::I) where {I}
    @boundscheck checkindex(indices, i)
    return i
end

function Base.setindex!(i::AbstractIndices{I}, ::I, ::I) where {I}
    error("Indices are not settable: $(typeof(i))")
end

function Base.isassigned(indices::AbstractIndices{I}, i::I) where {I}
    return i in indices
end

Base.keys(i::AbstractIndices) = i

@inline function Base.iterate(inds::AbstractIndices, s...)
    if istokenizable(inds)
        tmp = iteratetoken(keys(inds), s...)
        tmp === nothing && return nothing
        (t, s2) = tmp
        return (gettokenvalue(inds, t), s2)
    else
        error("All AbstractIndices must define `iterate`: $(typeof(inds))")
    end
end

function Base.in(i::I, indices::AbstractIndices{I}) where I
    if !istokenizable(indices)
        # A fallback definition based on iteration would be rediculously slow. There
        # shouldn't be many `AbstractIndex` types that rely on iteration for this.
        error("All AbstractIndices must define `in`: $(typeof(indices))")
    end

    return gettoken(indices, i)[1]
end

function Base.in(i, indices::AbstractIndices{I}) where I
    return convert(I, i) in indices
end

# Match the default setting from Base - the majority of containers will know their size
Base.IteratorSize(::AbstractIndices) = Base.HasLength()

function Base.length(indices::AbstractIndices)
    if Base.IteratorSize(indices) isa Base.SizeUnknown
        out = 0
        for _ in tokens(indices)
            out += 1
        end
        return out
    end

    error("All AbstractIndices must define `length` or else have `IteratorSize` of `SizeUnknown`: $(typeof(indices))")
end

Base.unique(i::AbstractIndices) = i

"""
    copy(inds::AbstractIndices)
    copy(inds::AbstractIndices, I::Type)

Construct a shallow copy of `inds`, possibly specifying a new element type `I`. The output
container is not guaranteed to be the same type as the input.
"""
Base.copy(inds::AbstractIndices) = copy(inds, eltype(inds))

function Base.copy(inds::AbstractIndices, ::Type{I}) where I
    out = empty(inds, I)
    for i in inds
        insert!(out, i)
    end
    return out
end

Base.empty(::AbstractIndices, ::Type{I}) where {I} = Indices{I}()

"""
    distinct(itr)

Collect the distinct elements of iterator `itr` into a new collection. Similar to
`Base.unique`, except returning a set (`Indices`) instead of an array.

# Example

```julia
julia> distinct([1,2,3,3])
3-element Indices{Int64}
 1
 2
 3
```
"""
distinct(itr) = _distinct(Indices, itr)
distinct(inds::AbstractIndices) = inds

function _distinct(::Type{T}, itr) where T
    out = T()
    union!(out, itr)
    return out
end

struct IndexError <: Exception
	msg::String
end

function Base.checkindex(indices::AbstractIndices{I}, i::I) where {I}
	if i ∉ indices
		throw(IndexError("Index $i not found in indices $indices"))
	end
end
Base.checkindex(indices::AbstractIndices{I}, i) where {I} = checkindex(indices, convert(I, i))

function checkindices(indices::AbstractIndices, inds)
    if !(inds ⊆ indices)
        throw(IndexError("Indices $inds are not a subset of $indices"))
    end
end

# Indices are isequal if they iterate in the same order
function Base.isequal(i1::AbstractIndices, i2::AbstractIndices)
    if sharetokens(i1, i2)
        return true
    end

    if length(i1) != length(i2)
        return false
    end

    for (j1, j2) in zip(i1, i2)
        if !isequal(j1, j2)
            return false
        end
    end

    return true
end

# The indices must be isequal and the values ==, same ordering
function Base.:(==)(i1::AbstractIndices, i2::AbstractIndices)
    out = true

    if sharetokens(i1, i2)
        # TODO - can we get rid of this loop for reflexive == element types?
        @inbounds for t in tokens(i1)
            # make sure it works for `missing`
            out &= gettokenvalue(i1, t) == gettokenvalue(i2, t)
            if out === false
                return false
            end
        end
        return out
    end

    if length(i1) != length(i2)
        return false
    end

    for (j1, j2) in zip(i1, i2)
        if !isequal(j1, j2)
            return false
        end
        out &= j1 == j2 # make sure it works for `missing`
        if out === false
            return false
        end
    end

    return out
end

# Lexical ordering based on iteration
function Base.isless(inds1::AbstractIndices, inds2::AbstractIndices)
    if sharetokens(inds1, inds2)
        return false # they are isequal
    end

    # We want to iterate... until one is longer than the other (no `zip`)
    tmp1 = iterate(inds1)
    tmp2 = iterate(inds2)
    while tmp1 !== nothing
        if tmp2 === nothing
            return false # shorter collections are isless in lexical ordering
        end
        (i1, s1) = tmp1
        (i2, s2) = tmp2
        !isequal(i1, i2)
        c = cmp(i1, i2)

        if c == -1
            return true
        elseif c == 1
            return false
        end

        tmp1 = iterate(inds1, s1)
        tmp2 = iterate(inds2, s2)
    end
    return tmp2 !== nothing # shorter collections are isless in lexical ordering, otherwise isequal
end

function Base.cmp(inds1::AbstractIndices, inds2::AbstractIndices)
    if sharetokens(inds1, inds2)
        return 0 # they are isequal
    end

    # We want to iterate... until one is longer than the other (no `zip`)
    tmp1 = iterate(inds1)
    tmp2 = iterate(inds2)
    while tmp1 !== nothing
        if tmp2 === nothing
            return 1 # shorter collections are isless in lexical ordering
        end
        (i1, s1) = tmp1
        (i2, s2) = tmp2
        !isequal(i1, i2)
        c = cmp(i1, i2)

        if c == -1
            return -1
        elseif c == 1
            return 1
        end

        tmp1 = iterate(inds1, s1)
        tmp2 = iterate(inds2, s2)
    end
    if tmp2 === nothing
        return 0
    else
        return -1 # shorter collections are isless in lexical ordering
    end
end

## Hashing - matches the dictionary case (matching keys and matching values)
function Base.hash(inds::AbstractIndices, h::UInt)
    h1 = h
    for i in inds
        h1 = hash(i, h1)
    end

    return hash(hash(UInt === UInt64 ? 0x8955a87bc313a509 : 0xa9cff5d1, h1), h1)
end

function Base.union(i::AbstractIndices, itr)
    if isinsertable(i)
        out = copy(i)
        union!(out, itr)
    else
        out = empty(i)
        union!(out, i)
        union!(out, itr)
    end
    return out
end

function Base.intersect(i::AbstractIndices, itr)
    if isinsertable(i)
        out = copy(i)
        intersect!(out, itr)
    else
        out = empty(i)
        union!(out, i)
        intersect!(out, itr)
    end
    return out
end

function Base.setdiff(i::AbstractIndices, itr)
    if isinsertable(i)
        out = copy(i)
        setdiff!(out, itr)
    else
        out = empty(i)
        union!(out, i)
        setdiff!(out, itr)
    end
    return out
end

function Base.symdiff(i::AbstractIndices, itr)
    if isinsertable(i)
        out = copy(i)
        symdiff!(out, itr)
    else
        out = empty(i)
        union!(out, i)
        symdiff!(out, itr)
    end
    return out
end

# issetequal and issubset(equal) should work already

function randtoken(i::AbstractIndices)
    return randtoken(Random.GLOBAL_RNG, i)
end

function randtoken(::Random.AbstractRNG, i::AbstractIndices)
    error("randtoken is not implemented for $(typeof(i))")
end

@static if VERSION < v"1.5-"
    """
        disjoint(set1, set2)

    Return `true` if `set1` and `set2` are disjoint or `false`. Two sets are disjoint if no
    elements of `set1` is in `set2`, and vice-versa. Somewhat equivalent to, but faster than,
    `isempty(intersect(set1, set2))`.
    """
    function disjoint(set1, set2)
        for i in set1
            if i in set2
                return false
            end
        end
        return true
    end
else
    @deprecate disjoint(set1, set2) isdisjoint(set1, set2)
end


text:
# TODO: parse expressions with MacroTools

"""
    @cfunc ex

Define a "calculation" function with an "adapter" methods used to parse args from model.
"calculate" functions take arguments from system state
and calculate "rate"s determining the probability weight of the reaction be selected.

# Example

For function definition:
```julia
@cfunc @inline growth_c(r, x::Vector) = r * x # function to calculate "growth rate"
```
This macro creates two methods, an "adapter" method
```julia
growth_c(args::NamedTuple) = growth_c(args.r, args.x)
```
and the origin method
```julia
@inline growth_c(r, x::Vector) = r * x
```

!!! info
    For function definition with other macros, put those macros after this macro.

!!! warning
    The argument name `t` is reserved for time.
    Don't use those variable names for other usage
"""
macro cfunc(ex)
    fname, fargs = parse_func(ex, (:t,))
    adapter = :(Base.@propagate_inbounds $fname(t, args::NamedTuple) = $fname($(fargs...)))
    return esc(Expr(:block, adapter, ex))
end

"""
    @ufunc ex

Define a "update" function with an "adapter" methods used to parse args from model.
"update" functions take arguments from system state and update system state.

# Example

For function definition:
```julia
@ufunc Base.@propagate_inbounds growth_u!(ind::CartesianIndex{1}, x) = x[ind] += 1
```
This macro creates two methods, an "adapter" method
```julia
growth_u!(ind, args::NamedTuple) = growth_u(ind, args.x)
```
and the origin method
```julia
Base.@propagate_inbounds growth_u!(ind::CartesianIndex{1}, x) = x[ind] += 1
```

!!! info
    For function definition with other macros, put those macros after this macro.

!!! warning
    The argument name `t` is reserved for time,
    arguments name `rng` is reserved for random number generator,
    and the argument name `ind` is preserved for index of "reaction".
    Don't use those variable names for other usage
"""
macro ufunc(ex)
    fname, fargs = parse_func(ex, (:t, :ind, :rng))
    adapter = :(Base.@propagate_inbounds $fname(rng, t, ind, args::NamedTuple) =
        $fname($(fargs...)))
    return esc(Expr(:block, adapter, ex))
end

"""
    @reaction name ex

Define a [`Reaction`](@ref) with the given `name` and `ex`.

```julia
@reaction growth begin
    r * x # "calculate" expression
    begin
        i = ind[1]
        x[i] += 1 # "update" expression
    end
end
```
will create a `Reaction` with a "calculation" function:
```julia
@cfunc Base.@propagate_inbounds growth_c(r, x) = r * x
```
and an "update" function:
```julia
@ufunc Base.@propagate_inbounds growth_u!(ind, x) = begin
    i = ind[1]
    x[i] += 1
end
```
where arguments of functions were collected from given expression automatically.

!!! note
    If there are global variables, this macro may can not collect arguments correctly.
    Especially, for functions which accept function and types as arguments,
    those functions and types may also be collect as an arguments.
    Thus, these variables must be marked as global variables by `global` before use them,
    even functions and types.
    Besides, type annotation expressions like `x::Vector{Int}`,
    types `Vector` and `Int` will not be collected.
    Avoid to defined your reaction with a type arguments for type annotation.

!!! warning
    The expression follow the same name preserve rule as `@cfunc` and `@ufunc`,
    don't use those variable names for other usage.

This macro cannot parse macro expression.
Thus in some cases, arguments of expression with macro may not be collected correctly.
To avoid this, define reaction with anonymous functions may helpfully:
```julia
@reaction reaction begin
    (A, B) -> @einsum C[i, k] := A[i, j] * B[j, k] # where C, i, j, k should not be collected
    begin
        # do something
    end
end
```
in this case, only arguments of anonymous function (`A` and `B`) will be collected.

Besides, there are another macro [`@quickloop`](@ref), which is similar to `@einsum`,
but it don't follow the einstein convenience which just generate nest loops out of the expression.
Because the `@quickloop` generate a function and arguments are provided,
indexing variables and return variable will not be collected.
"""
macro reaction(name::Symbol, block::Expr)
    block = rmlines(block)
    if block.head == :block && length(block.args) == 2
        cargs, cbody = parse_args_body(block.args[1])
        uargs, ubody = parse_args_body(block.args[2])
        return esc(gen_reaction(name, cargs, cbody, uargs, ubody))
    else
        throw(ArgumentError("can't parse expression"))
    end
end

function gen_reaction(name, cargs, cbody, uargs, ubody)
    cname = Symbol("$(name)_c")
    uname = Symbol("$(name)_u!")
    return Expr(
        :block,
        :(Base.@propagate_inbounds $cname(t, args::NamedTuple) =
            $cname($((_warparg(arg, (:t,)) for arg in cargs)...))),
        :($cname($(cargs...)) = $cbody),
        :(Base.@propagate_inbounds $uname(rng, t, ind, args::NamedTuple) =
            $uname($((_warparg(arg, (:t, :ind, :rng)) for arg in uargs)...))),
        :($uname($(uargs...)) = $ubody),
        Expr(:(=), name, :(Reaction($cname, $uname))),
    )
end

# extract name and arguments of given function expression
function parse_func(ex::Expr, preserve)
    head = ex.head
    if head in (:function, :(=), :where) # unpack function expression
        return parse_func(ex.args[1], preserve)
    elseif head == :macrocall # unpack macro expression
        return parse_func(ex.args[end], preserve)
    elseif head == :call # the args of :call expression is the function name and the arguments
        fname = ex.args[1]
        fargs = map(arg -> parse_arg(arg, preserve), ex.args[2:end])
        return fname, fargs
    else
        error("Unknown expr")
    end
end

# unpack argument names
parse_arg(arg::Symbol, preserve) = _warparg(arg, preserve)
parse_arg(arg::Expr, preserve) = _warparg(arg.args[1]::Symbol, preserve)
function _warparg(sym::Symbol, preserve)
    if sym in preserve
        return sym
    else
        return Expr(:., :args, QuoteNode(sym))
    end
end

# collect arguments from given expression
function parse_args_body(ex)
    if ex.head == :(->)
        return unbox_args!(Set{Symbol}(), ex.args[1]), unblock(ex.args[2])
    elseif ex.head == :macrocall && ex.args[1] == Symbol("@quickloop")
        args = filter(!isline, ex.args)
        return _quickloop(args[2:end]...)
    else
        return collectargs!(Set{Symbol}(), Set{Symbol}(), ex), ex
    end
end

function unbox_args!(args, ex::Expr)
    if ex.head == :(::)
        push!(args, ex.args[1])
    elseif ex.head == :tuple
        for arg in ex.args
            unbox_args!(args, arg)
        end
    else
        error("can't parse expression")
    end
    return args
end
unbox_args!(args, s::Symbol) = (push!(args, s), args)

# `args` is a list of symbols will be treat as arguments name
# `vars` is a list of symbols which are names of internal variable or global variables
collectargs!(args, _, _) = args # ignore when third argument is neither a symbol nor a expr
function collectargs!(args, vars, sym::Symbol) # symbols neither in syms or iv will be collected
    if !in(sym, args) && !in(sym, vars) && sym != :(:)
        push!(args, sym)
    end
    return args
end
function collectargs!(args, vars, _ex::Expr)
    ex = shortdef(_ex) # force convert function definition to short form
    head = ex.head
    if head == :(=) # assignment or function definition
        arg1 = ex.args[1]
        if (arg1 isa Symbol || (arg1 isa Expr && arg1.head == :tuple)) # assignment to one or multi new symbol
            collectargs!(args, vars, ex.args[2]) # must collect exprs right of equal firstly #2
            collectargs!(vars, args, ex.args[1]) # collect internal variables
            return args
        elseif arg1.head == :call # local function definition as short form
            push!(vars, arg1.args[1]) # local function name will be internal variable
            func_args = Set{Symbol}() # local function arguments
            collectargs!(func_args, Set{Symbol}(), arg1) # collect local function arguments
            collectargs!(args, union(vars, func_args), ex.args[2]) # collect variables used in local function
            return args
        else
            collectargs!(args, vars, ex.args[1]) # collect left side of equal
            collectargs!(args, vars, ex.args[2]) # collect right side of equal
            return args
        end
    elseif head == :-> # anonymous function
        func_args = Set{Symbol}() # anonymous function arguments
        collectargs!(func_args, Set{Symbol}(), ex.args[1]) # collect function arguments
        collectargs!(args, union(vars, func_args), ex.args[2]) # collect variables used in anonymous function
        return args
    elseif head == :. # for x.y only collect x, y will be ignore
        return collectargs!(args, vars, ex.args[1])
    elseif head == :global # global variables
        return union!(vars, ex.args)
    elseif head == :(::) # type annotation
        return collectargs!(args, vars, ex.args[1]) # collect left side of :: but ignore right side
    elseif head == :ref # ignore begin and end in reference
        collectargs!(args, vars, ex.args[1])
        foreach(
            arg -> collectargs!(args, TempUnion(vars, (:begin, :end)), arg),
            ex.args[2:end],
        )
    elseif head in (:call, :macrocall) # function and macro will not be treat as a args
        foreach(arg -> collectargs!(args, vars, arg), ex.args[2:end])
    else # in other cases, collect all arguments
        foreach(arg -> collectargs!(args, vars, arg), ex.args)
    end
    return args
end

# union two collections
# but push! to the first collection
struct TempUnion{S,T}
    s::S
    t::T
end
Base.in(item, s::TempUnion) = (item in s.s) || (item in s.t)
Base.push!(s::TempUnion, item) = push!(s.s, item)

"""
@quickloop ex opts...

Generate a function to calculate given expression with nested loops.

# Example
```julia
@quickloop x[i] * m[i, j] * y[j]
```
will generate a function
```julia
(m, y, x)->begin
    ret = Array{...}(undef, size(m, 1), size(m, 2)) # expression for type inference is omitted
    @fastmath @inbounds for j = axes(m, 2), i = axes(m, 1)
        ret[i, j] = x[i] * m[i, j] * y[j]
    end
    return ret
end
```
where the size of return array is determined automatically.
The size can be specified by
```
@quickloop z[j, i] := x[i] * m[i, j] * y[j]
```
which will generate a function
```julia
(m, y, x)->begin
    z = Array{...}(undef, size(m, 2), size(m, 1)) # expression for type inference is omitted
    @fastmath @inbounds for i = axes(m, 1), j = axes(m, 2)
        z[j, i] = x[i] * m[i, j] * y[j]
    end
    return z
end
```
where the return value is named by `z` and the size of return array is
`(size(m, 2), size(m, 1))` instead of `(size(m, 1), size(m, 2))`.

# Options
Options should be specified like `@quickloop ex opt=true/false` or `@quickloop ex opt`,
if the value of option is omitted, it will be treated as `true`.

- `avx` or `turbo`: use `LoopVectorization` to accelerate the loop, default is `false`,
    unless `LoopVectorization` have be loaded.
- `inbounds`: use `@inbounds` to accelerate the loop, default is `true`, this option is ignored
    if `avx` or `turbo` is enabled.
- `fastmath`: use `@fastmath` to accelerate the loop, default is `true`, this option is ignored
    if `avx` or `turbo` is enabled.
- `offset`: if `true`, input arrays will be treat as offset arrays, default is `false`.

!!! note
    The order of function arguments generated by this macro is undetermined.
    Thus, it's recommender use this macro inside of `@reaction`,
    which detect the order of function arguments automatically.
"""
macro quickloop(args...)
    args, fn_body = _quickloop(args...)
    fn_args = Expr(:tuple, args...)
    return esc(Expr(:(->), fn_args, fn_body))
end

function _quickloop(ex, opts...)
    lhs, rhs = _split_eq(ex)
    inds, inds_ax = collectaxes(lhs, rhs)
    fn_args = collectargs!(Set{Symbol}(), Set{Symbol}(inds), rhs)
    fn_body = gen_loop(inds, inds_ax, fn_args, lhs, rhs, opts...)
    return fn_args, fn_body
end

function _split_eq(ex)
    if ex.head == :(=) || ex.head == :(:=)
        return ex.args[1], ex.args[2]
    else
        return Expr(:ref, gensym(:ret)), ex
    end
end

function collectaxes(lhs, rhs)
    # Dict of index symbol and dict of array and dim of the index
    inds_ax = Dict{Symbol,Dict{Symbol,Vector{Int}}}()
    collectaxes!(rhs, inds_ax)
    if length(lhs.args) == 1
        inds = collect(keys(inds_ax))
        max_pos(ind) = maximum(Iterators.flatten(values(inds_ax[ind])))
        sort!(inds; by=max_pos)
        append!(lhs.args, inds)
    else
        inds = lhs.args[2:end]
        for ind in keys(inds_ax)
            if !in(ind, inds)
                throw(ArgumentError(
                    "index symbol `$ind` don't occur in left hand side `$lhs`"
                ))
            end
        end
    end
    return inds, inds_ax
end

function collectaxes!(ex::Expr, inds_ax)
    if ex.head == :ref
        array_name = ex.args[1]
        for i in 2:length(ex.args)
            ind = ex.args[i]
            if !haskey(inds_ax, ind)
                inds_ax[ind] = Dict{Symbol,Vector{Int}}()
            end
            if !haskey(inds_ax[ind], array_name)
                inds_ax[ind][array_name] = Vector{Int}()
            end
            push!(inds_ax[ex.args[i]][array_name], i - 1)
        end
    else
        for arg in ex.args
            collectaxes!(arg, inds_ax)
        end
    end
    return inds_ax
end
collectaxes!(::Any, axs) = axs

function gen_loop(inds, inds_ax, fn_args, lhs, rhs, opts...)
    if isempty(inds) # scalar expression
        return Expr(:return, rhs)
    end
    mode = copy(DEFAULT_LOOP_MODE)
    for ex in opts
        opt, flag = _parse_opt(ex)
        if opt == :turbo
            mode[1] = flag
        elseif opt == :avx
            mode[1] = flag
        elseif opt == :inbounds
            mode[2] = flag
        elseif opt == :fastmath
            mode[3] = flag
        # mode[4] is used to determine of check axes of given array is matched,
        # and it's removed because of the introduction of `ArrayInterface.indices`
        # this mode[4] maybe reused in the future
        elseif opt == :offset
            mode[5] = flag
        else
            error("unknown loop mode: $opt")
        end
    end
    ret_name = lhs.args[1]
    ret_size = Expr[]
    ret_type = _infer_ret_type(inds, inds_ax, fn_args, rhs)
    loop_head = Expr(:block)
    for (i, ind) in enumerate(inds)
        array_name, dims = first(inds_ax[ind])
        dim = first(dims)
        push!(ret_size, :(size($array_name, $dim)))
        array_tuple = Expr(:tuple, ret_name)
        dim_tuple = Expr(:tuple, i)
        for (a, ds) in inds_ax[ind]
            for d in ds
                push!(array_tuple.args, a)
                push!(dim_tuple.args, d)
            end
        end
        ax = Expr(:call, :indices, array_tuple, dim_tuple)
        pushfirst!(loop_head.args, Expr(:(=), ind, ax))
    end
    ret_init = :($ret_name = Array{$ret_type}(undef, $(ret_size...)))
    rhs_offset = if mode[5]
        MacroTools.postwalk(rhs) do ex
            if ex isa Expr && ex.head == :ref
                ex′ = deepcopy(ex)
                for i in 2:length(ex.args)
                    ex′.args[i] = :(axes($(ex.args[1]), $(i-1))[$(ex.args[i])])
                end
                return ex′
            else
                return ex
            end
        end
    else
        rhs
    end
    loop_body = Expr(:block, Expr(:(=), lhs, rhs_offset))
    loop = Expr(:for, loop_head, loop_body)
    if mode[1]
        loop = :(@turbo $loop)
    else
        if mode[2]
            loop = :(@inbounds $loop)
        end
        if mode[3]
            loop = :(@fastmath $loop)
        end
    end
    body = Expr(:block, ret_init, loop, Expr(:return, ret_name))
    return MacroTools.flatten(body)
end

# turbo, inbounds, fastmath, [removed], offset, ...
const DEFAULT_LOOP_MODE = Bool[0, 1, 1, 0, 0, 0, 0, 0]

_parse_opt(opt::Symbol) = opt, true
function _parse_opt(opt::Expr)
    if opt.head == :(=)
        return opt.args[1]::Symbol, opt.args[2]::Bool
    else
        throw(ArgumentError("options should be `opt` or `opt = value::Bool`"))
    end
end

function _infer_ret_type(inds, inds_ax, fn_args, ex)
    ex′ = Expr(:block, ex)
    for ind in inds
        array_name, dims = first(inds_ax[ind])
        dim = first(dims)
        assignment = :($ind = firstindex($array_name, $dim))
        pushfirst!(ex′.args, assignment)
    end
    args = collect(fn_args)
    promote_op_fn = Expr(:(->), Expr(:tuple, args...), ex′)
    args_type = map(arg -> :(typeof($arg)), args)
    ret_type = :(Base.promote_op($promote_op_fn, $(args_type...)))
    return ret_type # expression for inferring the return type
end

"""
    @reaction_eq name param eq [opts...]

Generate a reaction with the given name, parameters and equation.

# Example
```julia
@reaction growth r X[i] --> 2X[i]
```
is equivalent to
```julia
@reaction growth begin
    r * X
    begin
        i = ind[1]
        X[i] += 1
    end
end
```
"""
macro reaction_eq(name::Symbol, param, eq::Expr, opts...)
    lhs_counts, dif_count = parse_eq(eq)
    inds, cargs, cbody = gen_rate(name, param, lhs_counts, opts...)
    update_ex = gen_update(inds, dif_count)
    uargs, ubody = parse_args_body(update_ex)
    return esc(gen_reaction(name, cargs, cbody, uargs, ubody))
end

function parse_eq(eq::Expr)
    if Meta.isexpr(eq, :(-->), 2)
        return _parse_eq(eq.args[1], eq.args[2])
    elseif Meta.isexpr(eq, :call, 3) && eq.args[1] === :(→)
        return _parse_eq(eq.args[2], eq.args[3])
    else
        throw(ArgumentError("unknown expression: $eq, expected `-->` or `→`"))
    end
end
function _parse_eq(lhs, rhs)
    lhs_counts = countsyms(lhs)
    rhs_counts = countsyms(rhs)
    dif_counts = Dict{Any,Int}()
    syms = union(keys(lhs_counts), keys(rhs_counts))
    for sym in syms
        if !haskey(dif_counts, sym)
            dif_counts[sym] = 0
        end
        rhs_count = get(rhs_counts, sym, 0)
        lhs_count = get(lhs_counts, sym, 0)
        dif_counts[sym] = rhs_count - lhs_count
    end
    return lhs_counts, dif_counts
end

function gen_rate(name::Symbol, param, counts, opts...)
    ex = Expr(:call, :*, param)
    for (sym, n) in counts
        if n != 1
            push!(ex.args, :($sym^$n))
        else
            push!(ex.args, sym)
        end
    end
    lhs = Expr(:ref, Symbol(name, :_rate))
    inds, inds_ax = collectaxes(lhs, ex)
    fn_args = collectargs!(Set{Symbol}(), Set{Symbol}(inds), ex)
    fn_body = gen_loop(inds, inds_ax, fn_args, lhs, ex, opts...)
    return inds, fn_args, fn_body
end

function gen_update(inds, counts)
    update = quote
        $(Expr(:tuple, inds...)) = ind.I
    end
    for (sym, n) in counts
        if n != 0
            push!(update.args, :($(_maybe_ref(sym)) += $n))
        end
    end
    return update
end

_maybe_ref(sym::Symbol) = Expr(:ref, sym)
_maybe_ref(sym::Expr) = sym

function countsyms(ex)
    count_dict = Dict{Any, Int}()
    return countsyms!(count_dict, ex, 1)
end

function countsyms!(count_dict, ex::Expr, n::Int)
    if ex.head === :ref
        if haskey(count_dict, ex)
            count_dict[ex] += n
        else
            count_dict[ex] = n
        end
    elseif ex.head === :call && ex.args[1] == :+
        for arg in ex.args[2:end]
            countsyms!(count_dict, arg, n)
        end
    elseif ex.head === :call && ex.args[1] == :* && length(ex.args) == 3 && ex.args[2] isa Int
        countsyms!(count_dict, ex.args[3], ex.args[2] * n)
    else
        throw(ArgumentError("unexpected expression: $ex"))
    end
    return count_dict
end
function countsyms!(count_dict, sym::Symbol, n::Int)
    if haskey(count_dict, sym)
        count_dict[sym] += n
    elseif sym != :∅
        count_dict[sym] = n
    end
    return count_dict
end
countsyms!(count_dict, ::Any, ::Int) = count_dict


text:
layers(::Type{WorldClim{Climate}}) = (:tmin, :tmax, :tavg, :prec, :srad, :wind, :vapr)

"""
    getraster(T::Type{WorldClim{Climate}}, [layer::Union{Tuple,Symbol}]; month, res::String="10m") => Vector{String}

Download [`WorldClim`](@ref) [`Climate`](@ref) data.

# Arguments
- `layer` `Symbol` or `Tuple` of `Symbol` from `$(layers(WorldClim{Climate}))`.
    Without a `layer` argument, all layers will be downloaded, and a `NamedTuple` of paths returned.

# Keywords
- `month`: `Integer` or `AbstractArray` of `Integer`. Chosen from `1:12`.
- `res`: `String` chosen from $(resolutions(WorldClim{Climate})), "10m" by default.

Returns the filepath/s of the downloaded or pre-existing files.
"""
function getraster(T::Type{WorldClim{Climate}}, layers::Union{Tuple,Symbol};
    month, res::String=defres(T)
)
    _getraster(T, layers, month, res)
end

function _getraster(T::Type{WorldClim{Climate}}, layers, month::AbstractArray, res::String)
    _getraster.(T, Ref(layers), month, Ref(res))
end
function _getraster(T::Type{WorldClim{Climate}}, layers::Tuple, month::Integer, res::String)
    _map_layers(T, layers, month, res)
end
function _getraster(T::Type{WorldClim{Climate}}, layer::Symbol, month::Integer, res::String)
    _check_layer(T, layer)
    _check_res(T, res)
    raster_path = rasterpath(T, layer; res, month)
    if !isfile(raster_path)
        zip_path = zippath(T, layer; res, month)
        _maybe_download(zipurl(T, layer; res, month), zip_path)
        zf = ZipFile.Reader(zip_path)
        mkpath(dirname(raster_path))
        raster_name = rastername(T, layer; res, month)
        write(raster_path, read(_zipfile_to_read(raster_name, zf)))
        close(zf)
    end
    return raster_path
end

# Climate layers don't get their own folder
rasterpath(T::Type{<:WorldClim{Climate}}, layer; res, month) =
    joinpath(_rasterpath(T, layer), rastername(T, layer; res, month))
_rasterpath(T::Type{<:WorldClim{Climate}}, layer) = joinpath(rasterpath(T), string(layer))
rastername(T::Type{<:WorldClim{Climate}}, layer; res, month) =
    "wc2.1_$(res)_$(layer)_$(_pad2(month)).tif"
zipname(T::Type{<:WorldClim{Climate}}, layer; res, month=1) =
    "wc2.1_$(res)_$(layer).zip"
zipurl(T::Type{<:WorldClim{Climate}}, layer; res, month=1) =
    joinpath(WORLDCLIM_URI, "base", zipname(T, layer; res, month))
zippath(T::Type{<:WorldClim{Climate}}, layer; res, month=1) =
    joinpath(rasterpath(T), "zips", zipname(T, layer; res, month))

_pad2(month) = lpad(month, 2, '0')

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
# Structs
"""
    Limits

GitHub API limits.

It includes how many remaining queries are available for the current time period and when it resets.

# Fields
- `limit::UInt16`
- `remaining::UInt16`
- `reset::DateTime`
"""
mutable struct Limits
    limit::UInt16
    remaining::UInt16
    reset::DateTime
end
"""
    GitHubPersonalAccessToken(login::AbstractString,
                              token::AbstractString,
                              )::GitHubPersonalAccessToken

A GitHub Personal Access Token

# Fields

- `login::String`
- `token::String`
- `client::Client`
- `limits::Limits`
"""
struct GitHubPersonalAccessToken
    login::String
    token::String
    client::Client
    limits::NamedTuple{(:core, :search, :graphql),NTuple{3,Limits}}
    function GitHubPersonalAccessToken(login::AbstractString, token::AbstractString)
        client = GraphQLClient(
            GITHUB_GRAPHQL_ENDPOINT,
            auth = "bearer $token",
            headers = Dict("User-Agent" => login),
            )
        # Dummy values
        limit = Limits(0, 0, DateTime(now(utc_tz), UTC))
        limits = (core = limit, search = limit, graphql = limit)
        output = new(login, token, client, limits)
        # Update dummy values for actual ones
        update!(output)
    end
end
summary(io::IO, obj::GitHubPersonalAccessToken) = println(io, "GitHub Personal Access Token")
function show(io::IO, obj::GitHubPersonalAccessToken)
    summary(io, obj)
    println(io, "  login: $(obj.login)")
    println(io, "  core remaining: $(obj.limits.core.remaining)")
    println(io, "  core reset: $(obj.limits.core.reset)")
    println(io, "  graphql remaining: $(obj.limits.graphql.remaining)")
    println(io, "  graphql reset: $(obj.limits.graphql.reset)")
end
function update!(obj::GitHubPersonalAccessToken)
    response = request(
        "GET",
        "$GITHUB_REST_ENDPOINT/rate_limit",
        [
            "Accept" => "application/vnd.github.v3+json",
            "User-Agent" => obj.login,
            "Authorization" => "token $(obj.token)",
        ],
    )
    json = JSON3.read(response.body).resources
    obj.limits.core.remaining = json.core.remaining
    obj.limits.core.reset = DateTime(ZonedDateTime(unix2datetime(json.core.reset), utc_tz), UTC)
    obj.limits.search.remaining = json.search.remaining
    obj.limits.search.reset = DateTime(ZonedDateTime(unix2datetime(json.search.reset), utc_tz), UTC)
    obj.limits.graphql.remaining = json.graphql.remaining
    obj.limits.graphql.reset = DateTime(ZonedDateTime(unix2datetime(json.graphql.reset), utc_tz), UTC)
    obj
end
"""
    graphql(obj::GitHubPersonalAccessToken,
            operationName::AbstractString,
            vars::Dict{String};
            max_retries::Integer = 3)

Return JSON of the GraphQL query.
"""
function graphql(
    # obj::GitHubPersonalAccessToken = PARALLELENABLER.pat,
    query::AbstractString = query,
    operationName::AbstractString = string(match(r"(?<=query )\w+(?=[\(|\{])", query).match);
    vars::Dict{String} = Dict{String,Any}(),
    max_retries::Integer = 3,
    )
    obj = PARALLELENABLER.pat
    # operationName = match(r"(?<=query )\w+", query).match
    update!(obj)
    if iszero(obj.limits.graphql.remaining)
        w = obj.limits.graphql.reset - DateTime(now(utc_tz), UTC)
        sleep(max(w, zero(w)))
        obj.limits.graphql.remaining = obj.limits.graphql.limit
    end
    result = try
        obj.client.Query(query, operationName = operationName, vars = vars)
    catch err
        err
    end
    retries = 0
    # Either it exceeded the rate limit or it timeout (for those we retry threee times)
    isok = isa(result, Result) &&
        result.Info.status == 200 &&
        result.Data ≠ """{"errors":[{"type":"RATE_LIMITED","message":"API rate limit exceeded"}]}"""
    while !isok && retries < max_retries
        sleep(61)
        result = try
            obj.client.Query(query, operationName = operationName, vars = vars)
        catch err
            err
        end
        isok = isa(result, Result) &&
            result.Info.status == 200 &&
            result.Data ≠ """{"errors":[{"type":"RATE_LIMITED","message":"API rate limit exceeded"}]}"""
        retries += 1
    end
    # @assert isok (query = query, vars = vars, result = result)
    update!(obj)
    result
end
"""
    setup(;host::AbstractString = get(ENV, "PGHOST", "localhost"),
           port::AbstractString = get(ENV, "PGPORT", "5432"),
           dbname::AbstractString = get(ENV, "PGDATABASE", "postgres"),
           user::AbstractString = get(ENV, "PGUSER", "postgres"),
           password::AbstractString = get(ENV, "PGPASSWORD", "postgres"),
           schema::AbstractString = "gh_2007_\$(year(floor(now(utc_tz), Year) - Day(1)))",
           pats::Union{Nothing, Vector{GitHubPersonalAccessToken}} = nothing)

Sets up your PostgreSQL database for the project.

# Example

```julia-repl
julia> setup(pats = [GitHubPersonalAccessToken("MyGH_Login", ENV["GH_PAT"])])

```

```julia-repl
julia> setup()

```
"""
function setup(;host::AbstractString = get(ENV, "PGHOST", "localhost"),
                port::AbstractString = get(ENV, "PGPORT", "5432"),
                dbname::AbstractString = get(ENV, "PGDATABASE", "postgres"),
                user::AbstractString = get(ENV, "PGUSER", "postgres"),
                password::AbstractString = get(ENV, "PGPASSWORD", "postgres"),
                schema::AbstractString = "gh_2007_$(year(floor(now(utc_tz), Year) - Day(1)))",
                pats::Union{Nothing, Vector{GitHubPersonalAccessToken}} = nothing)
    GHOST.PARALLELENABLER.conn = Connection("host = $host port = $port dbname = $dbname user = $user password = $password")
    GHOST.PARALLELENABLER.schema = schema
    schema = "gh_2007_$(year(floor(now(utc_tz), Year) - Day(1)))"
    execute(GHOST.PARALLELENABLER.conn,
            replace(join(["CREATE EXTENSION IF NOT EXISTS btree_gist; CREATE SCHEMA IF NOT EXISTS schema;",
                          String(read(joinpath(@__DIR__, "assets", "sql", "pats.sql"))),
                          String(read(joinpath(@__DIR__, "assets", "sql", "licenses.sql"))),
                          String(read(joinpath(@__DIR__, "assets", "sql", "queries.sql"))),
                          String(read(joinpath(@__DIR__, "assets", "sql", "repos.sql"))),
                          String(read(joinpath(@__DIR__, "assets", "sql", "commits.sql"))),
                         ],
                        ' '),
                    "schema" => schema))
    if isnothing(pats)
        try
            pat = DataFrame(execute(GHOST.PARALLELENABLER.conn, "SELECT login, token FROM $schema.pats ORDER BY login LIMIT 1;"))
            GHOST.PARALLELENABLER.pat = only(GitHubPersonalAccessToken.(pat.login, pat.token))
        catch err
            throw(ArgumentError("No PAT was provided nor available in the database. You can provide PAT directly through `setup` using the keyword argument `pats`."))
        end
    else
        pats = DataFrame((login = pat.login, token = pat.token) for pat in pats)
        load!(pats, GHOST.PARALLELENABLER.conn, "INSERT INTO $(schema).pats VALUES(\$1, \$2) ON CONFLICT DO NOTHING;")
        GHOST.PARALLELENABLER.pat = GitHubPersonalAccessToken(values(first(pats))...)
    end
    nothing
end

"""
    setup_parallel(limit::Integer = 0; password::AbstractString = get(ENV, "PGPASSWORD", "postgres"))::Nothing

Setup workers.
"""
function setup_parallel(limit::Integer = 0; password::AbstractString = get(ENV, "PGPASSWORD", "postgres"))
    @unpack conn, schema = PARALLELENABLER
    io = IOBuffer()
    show(io, conn)
    s = String(take!(io))
    lns = strip.(split(s, '\n'))
    host = lns[findfirst(ln -> startswith(ln, "host = "), lns)]
    port = lns[findfirst(ln -> startswith(ln, "port = "), lns)]
    dbname = lns[findfirst(ln -> startswith(ln, "dbname = "), lns)]
    user = lns[findfirst(ln -> startswith(ln, "user = "), lns)]
    if limit > 0
        pat = DataFrame(execute(conn, "SELECT login, token FROM $schema.pats ORDER BY login LIMIT $limit;"))
    else
        pat = DataFrame(execute(conn, "SELECT login, token FROM $schema.pats ORDER BY login;"))
    end
    pats = GitHubPersonalAccessToken.(pat.login, pat.token)
    npats = length(pats)
    addprocs(npats, exeflags = `--proj`)
    remotecall_eval(Main, workers(), :(using GHOST))
    @everywhere workers() host = $host
    @everywhere workers() port = $port
    @everywhere workers() dbname = $dbname
    @everywhere workers() user = $user
    @everywhere workers() password = $password
    @everywhere workers() GHOST.PARALLELENABLER.conn = Connection("$host $port $dbname $user password = $password")
    @everywhere workers() GHOST.PARALLELENABLER.schema = $schema
    GHOST.READY.x = Vector{Future}(undef, npats)
    for proc ∈ workers()
        GHOST.READY.x[proc - 1] = GHOST.@spawnat proc nothing
        pat = pats[proc - 1]
        expr = :(GHOST.PARALLELENABLER.pat = $pat)
        remotecall_eval(Main, proc, expr)
    end
end


text:
using PyCall:PyObject, PyReverseDims, pyimport_conda

struct PyModel <: DetectorModel
    pyobject::PyObject
end

# lazily import a python model
pyod_import(name::String) = () -> begin
    raw_name = name[1:end-length("Detector")] # name without Detector ending
    module_name = lowercase(raw_name)
    # names are not consistently uppercased
    model_name = raw_name == "IForest" ? raw_name : uppercase(raw_name)
    getproperty(pyimport_conda("pyod.models.$(module_name)", "pyod", "conda-forge"), model_name)
end
pyod_import(name::Symbol) = pyod_import(String(name))

make_docs_link(name::String) =
    "<https://pyod.readthedocs.io/en/latest/pyod.models.html#module-pyod.models.$(lowercase(name))>"

# Implements the `fit` method for an underlying PyOD model.
function pyod_fit(modelname, params)
    pymodelname = String(modelname)
    quote
        function OD.fit(model::$modelname, X::Data; verbosity)::Fit
            Xt = PyReverseDims(X) # from column-major to row-major
            # load the underlying python model with key-word arguments
            detector = pyod_import($pymodelname)()($((Expr(:kw, p, :(model.$p)) for p in params)...))
            detector.fit(Xt)
            # the underlying python model is out model
            return PyModel(detector), detector.decision_scores_
        end
    end
end

# Implements the `score` method for an underlying PyOD model.
function pyod_score(modelname)
    quote
        function OD.transform(_::$modelname, model::DetectorModel, X::Data)::Scores
            Xt = PyReverseDims(X) # change from column-major to row-major
            scores_test = model.pyobject.decision_function(Xt)
            return scores_test
        end
    end
end

# Extracts the relevant information from the expr and build the expression corresponding to the model constructor.
function py_constructor(ex)
    # similar to @detector
    ex, modelname, params, defaults, constraints = OD._process_model_def(@__MODULE__, ex)

    # keyword constructor
    const_ex = OD._model_constructor(modelname, params, defaults)

    # associate the constructor with the definition of the struct
    push!(ex.args[3].args, const_ex)

    # cleaner
    clean_ex = OD._model_cleaner(modelname, defaults, constraints)

    # return
    return modelname, params, clean_ex, ex
end

macro pymodel(ex)
    modelname, params, clean_ex, ex = py_constructor(ex)
    expr = quote
        Base.@__doc__($ex)
        $clean_ex
        $(pyod_fit(modelname, params))
        $(pyod_score(modelname))
    end
    esc(expr)
end