text:
""" `numerical_deriv( x, y)`
Calculate simple estimate of numerical derivative
"""
function numerical_deriv( x::AA1, y::AA2 )   where { T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Real, AA2<:AbstractArray{T2,1}  }
	@assert length(x) == length(y)
	dfluxdlnλ = zeros(size(y))
	dfluxdlnλ[1] = (y[2]-y[1])/(x[2]-x[1])
	dfluxdlnλ[2:end-1] .= (y[3:end].-y[1:end-2])./(x[3:end].-x[1:end-2])
	dfluxdlnλ[end] = (y[end]-y[end-1])/(x[end]-x[end-1])
	return dfluxdlnλ
end

""" Estimate numerical derivative of fluxes given wavelengths. """
function calc_dfluxdlnlambda(flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real }
    @assert size(flux) == size(λ)
    @assert length(flux) >= 3
    dfdlogλ = Array{T1,1}(undef,length(flux))
    calc_dfluxdlnlambda!(dfdlogλ,flux,λ)
    #dfdlogλ[1] = 0.5 * (flux[2]-flux[1])/(λ[2]-λ[1])*(λ[2]+λ[1])
    #dfdlogλ[2:end-1] .= 0.5 .* (flux[3:end].-flux[1:end-2])./(λ[3:end].-λ[1:end-2]).*(λ[3:end].+λ[1:end-2])
    #dfdlogλ[end] = 0.5 * (flux[end]-flux[end-1])/(λ[end]-λ[end-1])*(λ[end]+λ[end-1])
    return dfdlogλ
end

""" Estimate numerical derivative of fluxes given wavelengths. """
function calc_dfluxdlnlambda!(dfdlogλ::AbstractArray{T1,1}, flux::AbstractArray{T2,1}, λ::AbstractArray{T3,1}) where { T1<:Real, T2<:Real, T3<:Real }
    @assert size(flux) == size(λ)
    @assert size(dfdlogλ) == size(flux)
    @assert length(flux) >= 3
    #dfdlogλ = Array{T1,1}(undef,length(flux))
    dfdlogλ[1] = 0.5 * (flux[2]-flux[1])/(λ[2]-λ[1])*(λ[2]+λ[1])
    dfdlogλ[2:end-1] .= 0.5 .* (flux[3:end].-flux[1:end-2])./(λ[3:end].-λ[1:end-2]).*(λ[3:end].+λ[1:end-2])
    dfdlogλ[end] = 0.5 * (flux[end]-flux[end-1])/(λ[end]-λ[end-1])*(λ[end]+λ[end-1])
    return dfdlogλ
end

""" Estimate numerical second derivative of fluxes given wavelengths. """
function calc_d2fluxdlnlambda2(flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real }
    @assert size(flux) == size(λ)
    @assert length(flux) >= 3
    logλ = log.(λ)
    d2fdlogλ2 = Array{T1,1}(undef,length(flux))
    #d2fdlogλ2[2:end-1] .= 0.25*(flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1])./(λ[3:end].+λ[1:end-2].-2.0*λ[2:end-1]).*(λ[3:end].+λ[end-2]).^2
    #d2fdlogλ2[2:end-1] .= 0.5 * (flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1]).* ((λ[3:end].+λ[1:end-2])./(logλ[3:end].-logλ[1:end-2])).^2
    #d2fdlogλ2[end] = d2fdlogλ2[end-1]
    #d2fdlogλ2[1] = d2fdlogλ2[2]
    calc_d2fluxdlnlambda2!(d2fdlogλ2,flux,λ)
    return d2fdlogλ2
end

""" Estimate numerical second derivative of fluxes given wavelengths. """
function calc_d2fluxdlnlambda2!(d2fdlogλ2::AbstractArray{T1,1}, flux::AbstractArray{T1,1}, λ::AbstractArray{T2,1}) where { T1<:Real, T2<:Real, T3<:Real }
    @assert size(flux) == size(λ)
    @assert size(d2fdlogλ2) == size(flux)
    @assert length(flux) >= 3
    #logλ = log.(λ)
    #d2fdlogλ2 = Array{T1,1}(undef,length(flux))
    #d2fdlogλ2[2:end-1] .= 0.25*(flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1])./(λ[3:end].+λ[1:end-2].-2.0*λ[2:end-1]).*(λ[3:end].+λ[end-2]).^2
    d2fdlogλ2[2:end-1] .= 0.5 * (flux[3:end].+flux[1:end-2].-2.0.*flux[2:end-1]).* ((λ[3:end].+λ[1:end-2]).^2 ./ (2.0.*log.(λ[3:end]./λ[1:end-2])))
    d2fdlogλ2[end] = d2fdlogλ2[end-1]
    d2fdlogλ2[1] = d2fdlogλ2[2]
    return d2fdlogλ2
end


text:
"Global Lax Friedrichs Flux"
function glf_num_flux(ul, ur, f, α)
    return 0.5*(f(ul)+f(ur))-α/2*(ur-ul)
end

"Flux for advection equation (simple upwinding)"
advection_num_flux(ul, ur, f, α) = ul

"""
Returns the Rusanov interface flux for the Euler equations
V. V. Rusanov, Calculation of Interaction of Non-Steady Shock Waves with
Obstacles, J. Comput. Math. Phys. USSR, 1, pp. 267-279, 1961.
"""
function rusanov_euler_num_flux(ul, ur, f, α)
  @assert(size(ul,1) == size(ur,1),
  "ul and ur vector have different number of components")
  F = zeros(ul)

  # Fisical variables and sound speeds
  ρL = ul[1];vL=ul[2]/ρL;EL = ul[3]
  pL = (γ - 1)*(EL - 0.5*ρL*vL*vL)
  aL = sqrt(γ*max(pL/ρL,1.0))

  ρR = ur[1];vR=ur[2]/ρR;ER = ur[3]
  pR = (γ - 1)*(ER - 0.5*ρR*vR*vR)
  aR = sqrt(γ*max(pR/ρR,1.0))

  # Find the maximum eigenvalue for each interface
  maxρ = max(abs(vL)+aL, abs(vR)+aR)

  # f₁ = ρ*u
  F[1] = 0.5*(ρL*vL + ρR*vR - maxρ*(ρR - ρL))
  # f₂ = rho*u*u+p
  F[2] = 0.5*(ρL*vL*vL + pL + ρR*vR*vR +
                   pR - maxρ*(ρR*vR - ρL*vL))
  # f₃ = (E+p)*u
  F[3] = 0.5*((EL + pL)*vL + (ER + pR)*vR - maxρ*(ER - EL))
  return F
end

ext: jl

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
module RTESolver

using StaticArrays
using CUDA
using Adapt
using UnPack

using ..Device: array_type, array_device, CUDADevice, CPU
using ..AngularDiscretizations
using ..Vmrs
using ..AtmosphericStates
using ..Fluxes
using ..BCs
using ..Sources
using ..RTE
using ..Optics
using ..LookUpTables

export solve_lw!, solve_sw!

include("RTESolverKernels.jl")

"""
    solve_lw!(
        slv::Solver{FT,I},
        max_threads::I,
        lkp_args...,
    ) where {I<:Int,FT<:AbstractFloat}

Solver for the longwave radiation problem
"""
function solve_lw!(
    slv::Solver{FT,I},
    max_threads::I,
    lkp_args...,
) where {I<:Int,FT<:AbstractFloat}
    @unpack as, op, bcs_lw, src_lw, flux_lw, fluxb_lw = slv
    @unpack nlay, ncol = as

    nargs = length(lkp_args)
    @assert nargs < 3
    if nargs > 0
        @assert lkp_args[1] isa LookUpLW
        if nargs == 2
            @assert lkp_args[2] isa LookUpCld
        end
        n_gpt = lkp_args[1].n_gpt
        major_gpt2bnd = Array{I,1}(lkp_args[1].major_gpt2bnd) #TODO temp fix to avoid scalar indexing
        flux = fluxb_lw
    else
        n_gpt = 1
        flux = flux_lw
    end
    ibnd = 1
    set_flux_to_zero!(flux_lw)

    for igpt = 1:n_gpt
        if nargs > 0 && lkp_args[1] isa LookUpLW
            ibnd = major_gpt2bnd[igpt]
        end
        # computing optical properties
        compute_optical_props!(op, as, src_lw, igpt, lkp_args...)

        if op isa OneScalar
            rte_lw_noscat_solve!(
                flux,
                src_lw,
                bcs_lw,
                op,
                nlay,
                ncol,
                igpt,
                ibnd,
                max_threads,
            ) # no-scattering solver
        else
            rte_lw_2stream_solve!(
                flux,
                src_lw,
                bcs_lw,
                op,
                nlay,
                ncol,
                igpt,
                ibnd,
                max_threads,
            ) # 2-stream solver
        end

        nargs > 0 && add_to_flux!(flux_lw, flux)
    end
    flux_lw.flux_net .= flux_lw.flux_up .- flux_lw.flux_dn
    return nothing
end

"""
    rte_lw_noscat_solve!(
        flux::FluxLW{FT},
        src_lw::SourceLWNoScat{FT},
        bcs_lw::LwBCs{FT},
        op::OneScalar{FT},
        nlay,
        ncol,
        igpt,
        ibnd,
        max_threads,
    ) where {FT<:AbstractFloat}

No scattering solver for the longwave problem.
LW fluxes, no scattering, mu (cosine of integration angle) specified by column
Does radiation calculation at user-supplied angles; converts radiances to flux
using user-supplied weights. Currently, only n_gauss_angles = 1 supported.
"""
function rte_lw_noscat_solve!(
    flux::FluxLW{FT},
    src_lw::SourceLWNoScat{FT},
    bcs_lw::LwBCs{FT},
    op::OneScalar{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
    max_threads,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    device = array_device(op.τ)
    if device === CUDADevice() # launch CUDA kernel
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (flux, src_lw, bcs_lw, op, nlay, ncol, igpt, ibnd)
        @cuda threads = (tx) blocks = (bx) rte_lw_noscat_solve_CUDA!(args...)
    else # use Julia native multithreading
        @inbounds Threads.@threads for gcol = 1:ncol
            rte_lw_noscat_source!(src_lw, op, gcol, nlay)
            rte_lw_noscat_transport!(
                src_lw,
                bcs_lw,
                op,
                gcol,
                flux,
                igpt,
                ibnd,
                nlay,
                nlev,
            )
        end
    end
    #------------------------------------------------
    return nothing
end
# ---------------------------------------------------------------
function rte_lw_noscat_solve_CUDA!(
    flux::FluxLW{FT},
    src_lw::SourceLWNoScat{FT},
    bcs_lw::LwBCs{FT},
    op::OneScalar{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    if gcol ≤ ncol
        rte_lw_noscat_source!(src_lw, op, gcol, nlay)
        rte_lw_noscat_transport!(
            src_lw,
            bcs_lw,
            op,
            gcol,
            flux,
            igpt,
            ibnd,
            nlay,
            nlev,
        )
    end
    return nothing
end

"""
    rte_lw_2stream_solve!(
        flux::FluxLW{FT},
        src_lw::SourceLW2Str{FT},
        bcs_lw::LwBCs{FT},
        op::TwoStream{FT},
        nlay,
        ncol,
        igpt,
        ibnd,
        max_threads,
    ) where {FT<:AbstractFloat}

Two stream solver for the longwave problem.

RRTMGP provides source functions at each level using the spectral mapping
of each adjacent layer. Combine these for two-stream calculations
lw combine sources (rte_lw_2stream_combine_sources!)
combine RRTMGP-specific sources at levels
compute layer reflectance, transmittance
compute total source function at levels using linear-in-tau (rte_lw_2stream_source!)
transport (adding_lw!)
"""
function rte_lw_2stream_solve!(
    flux::FluxLW{FT},
    src_lw::SourceLW2Str{FT},
    bcs_lw::LwBCs{FT},
    op::TwoStream{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
    max_threads,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    device = array_device(op.τ)
    if device === CUDADevice() # launch CUDA kernel
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (flux, src_lw, bcs_lw, op, nlay, ncol, igpt, ibnd)
        @cuda threads = (tx) blocks = (bx) rte_lw_2stream_solve_CUDA!(args...)
    else # use Julia native multithreading
        Threads.@threads for gcol = 1:ncol
            rte_lw_2stream_combine_sources!(src_lw, gcol, nlev, ncol)
            rte_lw_2stream_source!(op, src_lw, gcol, nlay, ncol)
            adding_lw!(flux, src_lw, bcs_lw, gcol, igpt, ibnd, nlev, ncol)
        end
    end
    #------------------------------------------------
    return nothing
end
# -------------------------------------------------------------------------------------------------
function rte_lw_2stream_solve_CUDA!(
    flux::FluxLW{FT},
    src_lw::SourceLW2Str{FT},
    bcs_lw::LwBCs{FT},
    op::TwoStream{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    if gcol ≤ ncol
        rte_lw_2stream_combine_sources!(src_lw, gcol, nlev, ncol)
        rte_lw_2stream_source!(op, src_lw, gcol, nlay, ncol)
        adding_lw!(flux, src_lw, bcs_lw, gcol, igpt, ibnd, nlev, ncol)
    end
end

"""
    solve_sw!(
        slv::Solver{FT,I},
        max_threads::I,
        lkp_args...,
    ) where {I<:Int,FT<:AbstractFloat}

Solver for the shortwave radiation problem
"""
function solve_sw!(
    slv::Solver{FT,I},
    max_threads::I,
    lkp_args...,
) where {I<:Int,FT<:AbstractFloat}
    @unpack as, op, bcs_sw, src_sw, flux_sw, fluxb_sw = slv
    @unpack nlay, ncol = as

    nargs = length(lkp_args)
    @assert nargs < 3
    if nargs > 0
        @assert lkp_args[1] isa LookUpSW
        if nargs == 2
            @assert lkp_args[2] isa LookUpCld
        end
        n_gpt = lkp_args[1].n_gpt
        major_gpt2bnd = Array{I,1}(lkp_args[1].major_gpt2bnd)#TODO temp fix to avoid scalar indexing
        solar_src_scaled = Array{FT,1}(lkp_args[1].solar_src_scaled)
        flux = fluxb_sw
    else
        n_gpt = 1
        flux = flux_sw
    end
    ibnd = 1
    solar_frac = FT(1)
    set_flux_to_zero!(slv.flux_sw)
    for igpt = 1:n_gpt
        if nargs > 0 && lkp_args[1] isa LookUpSW
            ibnd = major_gpt2bnd[igpt]
            solar_frac = solar_src_scaled[igpt]
        end
        # computing optical properties
        compute_optical_props!(op, as, igpt, lkp_args...)

        # solving radiative transfer equation
        if slv.op isa OneScalar
            rte_sw_noscat_solve!(
                flux,
                op,
                bcs_sw,
                nlay,
                ncol,
                igpt,
                ibnd,
                solar_frac,
                max_threads,
            ) # no-scattering solver
        else
            rte_sw_2stream_solve!(
                flux,
                op,
                bcs_sw,
                src_sw,
                nlay,
                ncol,
                igpt,
                ibnd,
                solar_frac,
                max_threads,
            ) # 2-stream solver
            flux.flux_dn .+= flux.flux_dn_dir
        end
        nargs > 0 && add_to_flux!(flux_sw, flux)
    end
    if slv.op isa TwoStream
        flux_sw.flux_net .= flux_sw.flux_up .- flux_sw.flux_dn
    end

    return nothing
end

"""
    rte_sw_noscat_solve!(
        flux::FluxSW{FT},
        op::OneScalar{FT},
        bcs_sw::SwBCs{FT},
        nlay,
        ncol,
        igpt,
        ibnd,
        solar_frac::FT,
        max_threads,
    ) where {FT<:AbstractFloat}

No-scattering solver for the shortwave problem.
(Extinction-only i.e. solar direct beam)
"""
function rte_sw_noscat_solve!(
    flux::FluxSW{FT},
    op::OneScalar{FT},
    bcs_sw::SwBCs{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
    solar_frac::FT,
    max_threads,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    device = array_device(op.τ)
    if device === CUDADevice() # launching CUDA kernel
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (flux, op, bcs_sw, nlay, ncol, igpt, ibnd, solar_frac)
        @cuda threads = (tx) blocks = (bx) rte_sw_noscat_solve_CUDA!(args...)
    else # launch juila native multithreading
        @inbounds Threads.@threads for gcol = 1:ncol
            rte_sw_noscat_solve_kernel!(
                flux,
                op,
                bcs_sw,
                solar_frac,
                gcol,
                nlev,
            )
        end
    end
    #-----------------------------------------------------------------------------------
    return nothing
end
#--------------------------------------------------------------------------------------------------
function rte_sw_noscat_solve_CUDA!(
    flux::FluxSW{FT},
    op::OneScalar{FT},
    bcs_sw::SwBCs{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
    solar_frac::FT,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    if gcol ≤ ncol
        rte_sw_noscat_solve_kernel!(flux, op, bcs_sw, solar_frac, gcol, nlev)
    end
    return nothing
end
#--------------------------------------------------------------------------------------------------
"""
    rte_sw_2stream_solve!(
        flux::FluxSW{FT},
        op::TwoStream{FT},
        bcs_sw::SwBCs{FT},
        src_sw::SourceSW2Str{FT},
        nlay,
        ncol,
        igpt,
        ibnd,
        solar_frac::FT,
        max_threads,
    ) where {FT<:AbstractFloat}

Two stream solver for the shortwave problem.
"""
function rte_sw_2stream_solve!(
    flux::FluxSW{FT},
    op::TwoStream{FT},
    bcs_sw::SwBCs{FT},
    src_sw::SourceSW2Str{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
    solar_frac::FT,
    max_threads,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    device = array_device(op.τ)
    if device === CUDADevice()
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (flux, op, bcs_sw, src_sw, nlay, ncol, igpt, ibnd, solar_frac)
        @cuda threads = (tx) blocks = (bx) rte_sw_2stream_solve_CUDA!(args...)
    else # launch julia native multithreading
        @inbounds Threads.@threads for gcol = 1:ncol
            sw_two_stream!(op, src_sw, bcs_sw, gcol, nlay) # Cell properties: transmittance and reflectance for direct and diffuse radiation
            sw_source_2str!(src_sw, bcs_sw, gcol, flux, solar_frac, ibnd, nlay) # Direct-beam and source for diffuse radiation
            adding_sw!(src_sw, bcs_sw, gcol, flux, igpt, ibnd, nlev) # Transport
        end
    end
    #--------------------------------------
    return nothing
end
#--------------------------------------------------------------------------------------------------
function rte_sw_2stream_solve_CUDA!(
    flux::FluxSW{FT},
    op::TwoStream{FT},
    bcs_sw::SwBCs{FT},
    src_sw::SourceSW2Str{FT},
    nlay,
    ncol,
    igpt,
    ibnd,
    solar_frac::FT,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    if gcol ≤ ncol
        sw_two_stream!(op, src_sw, bcs_sw, gcol, nlay) # Cell properties: transmittance and reflectance for direct and diffuse radiation
        sw_source_2str!(src_sw, bcs_sw, gcol, flux, solar_frac, ibnd, nlay) # Direct-beam and source for diffuse radiation
        adding_sw!(src_sw, bcs_sw, gcol, flux, igpt, ibnd, nlev) # Transport
    end
end
#--------------------------------------------------------------------------------------------------
end


text:
"""
    macro _make_logging_group(sym::Symbol, doc::String)

This creates a group of functions that allow masking out topic-related logging
actions. A call that goes as follows:

    @_make_logging_tag XYZ

creates the following tools:

- global variable `_XYZ_log_enabled` defaulted to false
- function `log_XYZ` that can be called to turn the logging on/off
- a masking macro `@_XYZ_log` that can be prepended to commands that should
  only happen if the logging of tag XYZ is enabled.

The masking macro is then used as follows:

    @_XYZ_log @info "This is the extra verbose information you wanted!" a b c

The user can direct logging with these:

    log_XYZ()
    log_XYZ(false)

`doc` should be a name of the stuff that is being printed if the corresponding
log_XYZ() is enabled -- it is used to create a friendly documentation for the
logging switch. In this case it could say `"X, Y and Z-related messages"`.
"""
macro _make_logging_tag(sym::Symbol, doc::String)
    enable_flag = Symbol(:_, sym, :_log_enabled)
    enable_fun = Symbol(:log_, sym)
    log_macro = Symbol(:_, sym, :_log)
    # esc() is necessary here because the internal macro processing would
    # otherwise bind the variables incorrectly.
    esc(:(
        begin
            $enable_flag = false

            """
                $(string($enable_fun))(enable::Bool=true)

            Enable (default) or disable (by passing `false`) output of $($doc).
            """
            $enable_fun(x::Bool = true) = begin
                global $enable_flag = x
            end

            macro $log_macro(x)
                $enable_flag ? x : :nothing
            end
        end
    ))
end

@_make_logging_tag models "model-related messages"
@_make_logging_tag io "messages and warnings from model input/output"
@_make_logging_tag perf "performance-related tracing information"


text:
isdefined(Base, :__precompile__) && __precompile__()

module BlsData

using Requests
using DataFrames
import Requests: post
import JSON

export Bls, api_url, set_api_url!, api_key, api_version, requests_made, requests_remaining
export BlsSeries, id, series, catalog
export get_data

const DEFAULT_API_URL            = "http://api.bls.gov/publicAPI/v2/timeseries/data/"
const API_KEY_LENGTH             = 32
const BLS_RESPONSE_SUCCESS       = "REQUEST_SUCCEEDED"
const BLS_RESPONSE_CATALOG_FAIL1 = "unable to get catalog data"
const BLS_RESPONSE_CATALOG_FAIL2 = "catalog has been disabled"
const LIMIT_DAILY_QUERY          = [25, 500]
const LIMIT_SERIES_PER_QUERY     = [25, 50]
const LIMIT_YEARS_PER_QUERY      = [10, 20]

"""
A connection to the BLS API.

Constructors
------------
* `Bls()`
* `Bls(key::AbstractString)`

Arguments
---------
* `key`: Registration key provided by the BLS.

Methods
-------
* `api_url(b::Bls)`: Get the base URL used to connect to the server
* `set_api_url(b::Bls, url::AbstractString)`: Set the base URL used to connect to the server
* `api_key(b::Bls)`: Get the API key
* `api_version(b::Bls)`: Get the API version (v1 or v2) used
* `requests_made(b::Bls)`: Get the number of requests made today
* `requests_remaining(b::Bls)`: Get the number of requests remaining today

Notes
-----
* A default API key can be specified in a ~/.blsdatarc file.
* A valid registration key increases the allowable number of requests per day as well making
  catalog metadata available.

"""
type Bls
    url::AbstractString
    key::AbstractString
    n_requests::Int16
    t_created::DateTime
end
function Bls(key="")
    if isempty(key)
        try
            open(joinpath(homedir(),".blsdatarc"), "r") do f
                key = readall(f)
            end
            key = rstrip(key)
            @printf "API key loaded.\n"
            # Key validation
            if length(key) > API_KEY_LENGTH
                key = key[1:API_KEY_LENGTH]
                warn("Key too long. First ", API_KEY_LENGTH, " chars used.")
            end
            if !isxdigit(key)
                error("Invalid BLS API key: ", key)
            end
        catch err
        end

    end

    url = DEFAULT_API_URL
    n_requests = 0
    t_created = now()
    Bls(url, key, n_requests, t_created)
end
api_url(b::Bls) = b.url
set_api_url!(b::Bls, url::AbstractString) = setfield!(b, :url, url)
api_key(b::Bls) = b.key
api_version(b::Bls) = 1 + !isempty(api_key(b))
requests_made(b::Bls) = b.n_requests
requests_remaining(b::Bls) = LIMIT_DAILY_QUERY[api_version(b)] - requests_made(b)
function increment_requests(b::Bls)
    # Reset request if we are in a new day!
    if Dates.day(now()) ≠ Dates.day(b.t_created)
        b.t_created = now()
        b.n_requests = 0
    end

    b.n_requests += 1
end

function Base.show(io::IO, b::Bls)
    @printf io "BLS API v%d Connection\n"   api_version(b)
    @printf io "\turl: %s\n"                api_url(b)
    @printf io "\tkey: %s\n"                api_key(b)
    @printf io "\trequests made (this connection): %d\n"      requests_made(b)
    @printf io "\trequests remaining (this connection): %d\n" requests_remaining(b)
end

"""
A time series with metadata returned from a `get_data` call.

Prefer to access fields with
```
id(s::BlsSeries)
series(s::BlsSeries)
catalog(s::BlsSeries)
```
"""
type BlsSeries
    id::AbstractString
    df::DataFrame
    catalog::AbstractString
end
id(s::BlsSeries)      = s.id
series(s::BlsSeries)  = s.df
catalog(s::BlsSeries) = s.catalog

function Base.show(io::IO, s::BlsSeries)
    @printf io "BlsSeries\n"
    @printf io "\tid: %s\n" id(s)
    @printf io "\tseries: %dx%d DataFrame with columns %s\n" size(series(s))...  names(series(s))
    @printf io "\tcatalog: %s\n" catalog(s)
end

EMPTY_RESPONSE() = BlsSeries("",DataFrame(),"")
function Base.isempty(s::BlsSeries)
    return id(s) == "" &&
        series(s) == DataFrame() &&
        catalog(s) == ""
end

include("get_data.jl")
end # module


text:
using CmdStan, Turing, ContinuousBenchmarks
using Mamba: describe

const dyes ="
data {
  int BATCHES;
  int SAMPLES;
  real y[BATCHES, SAMPLES];
  // vector[SAMPLES] y[BATCHES];
}
parameters {
  real<lower=0> tau_between;
  real<lower=0> tau_within;
  real theta;
  real mu[BATCHES];
}
transformed parameters {
  real sigma_between;
  real sigma_within;
  sigma_between <- 1/sqrt(tau_between);
  sigma_within <- 1/sqrt(tau_within);
}
model {
  theta ~ normal(0.0, 1E5);
  tau_between ~ gamma(.001, .001);
  tau_within ~ gamma(.001, .001);
  mu ~ normal(theta, sigma_between);
  for (n in 1:BATCHES)
    y[n] ~ normal(mu[n], sigma_within);
}
generated quantities {
  real sigmasq_between;
  real sigmasq_within;

  sigmasq_between <- 1 / tau_between;
  sigmasq_within <- 1 / tau_within;
}
"

const dyesdata = [
  Dict("BATCHES" => 6,
    "SAMPLES" => 5,
    "y" => reshape([
      [1545, 1540, 1595, 1445, 1595];
      [1520, 1440, 1555, 1550, 1440];
      [1630, 1455, 1440, 1490, 1605];
      [1595, 1515, 1450, 1520, 1560];
      [1510, 1465, 1635, 1480, 1580];
      [1495, 1560, 1545, 1625, 1445]
    ], 6, 5)
  )
]

#global stanmodel, rc, sim

#stanmodel = Stanmodel(Sample(algorithm=CmdStan.Hmc(CmdStan.Static(0.38 * 11), CmdStan.diag_e(), 0.38, 0.0),
#  save_warmup=true, adapt=CmdStan.Adapt(engaged=false)),
#  num_samples=2000, num_warmup=0, thin=1,
#  name="dyes", model=dyes, nchains=1);
#rc, sim = stan(stanmodel, dyesdata, CmdStanDir=ContinuousBenchmarks.CMDSTAN_HOME, summary=false)

#stanmodel = Stanmodel(name="dyes", model=dyes, useMamba=false)
#rc, sim = stan(stanmodel, dyesdata, CmdStanDir=ContinuousBenchmarks.CMDSTAN_HOME)

@model dyes_turing(BATCHES, SAMPLES, y) = begin
  theta ~ Normal(0.0, 1E5)
  tau_between ~ Gamma(.001, .001)
  tau_within ~ Gamma(.001, .001)
  mu = tzeros(BATCHES)
  mu ~ [Normal(theta, 1/sqrt(tau_between))]
  for n = 1:BATCHES
    y[n,:] ~ MvNormal(mu[n] .* ones(SAMPLES), 1/sqrt(tau_within) .* ones(SAMPLES))
  end
end

_batches = dyesdata[1]["BATCHES"]
_samples = dyesdata[1]["SAMPLES"]
_y = dyesdata[1]["y"]
chn = sample(dyes_turing(_batches, _samples, _y), NUTS(2000, 0.65))

describe(chn)


text:
## Could tidy up this HTML to make it look nicer
html_templates = Dict()

# thumbs up/down don't show in my editor
grading_partial = """
  if(correct) {
    msgBox.innerHTML = "
 👍  {{#:CORRECT}}{{{:CORRECT}}}{{/:CORRECT}}{{^:CORRECT}}Correct{{/:CORRECT}}
";
    var explanation = document.getElementById("explanation_{{:ID}}")
    if (explanation != null) {
       explanation.style.display = "none";
    }
  } else {
    msgBox.innerHTML = "
👎  {{#:INCORRECT}}{{{:INCORRECT}}}{{/:INCORRECT}}{{^:INCORRECT}}Incorrect{{/:INCORRECT}}
";
    var explanation = document.getElementById("explanation_{{:ID}}")
    if (explanation != null) {
       explanation.style.display = "block";
    }
  }
"""

## Basic question
## has label and hint option.
## Hint is put with label when present; otherwise, it appears at bottom of form.
## this is overridden with input widget in how show method is called
html_templates["question_tpl"] = mt"""
<form class="mx-2 my-3" name='WeaveQuestion' data-id='{{:ID}}' data-controltype='{{:TYPE}}'>

    {{#:LABEL}}
        {{{:LABEL}}}{{#:HINT}} 🎁{{/:HINT}}

    {{/:LABEL}}

    {{{:FORM}}}
    {{^:LABEL}}{{#:HINT}} 🎁{{/:HINT}}{{/:LABEL}}

      {{#:EXPLANATION}}

{{{:EXPLANATION}}}

      {{/:EXPLANATION}}

</form>

"""

html_templates["input_grading_script"] = jmt"""
document.getElementById("{{:ID}}").addEventListener("change", function() {
  var correct = {{{:CORRECT_ANSWER}}};
  var msgBox = document.getElementById('{{:ID}}_message');
  $(grading_partial)
});
"""

##
html_templates["inputq_form"] = mt"""

    <input id="{{:ID}}" type="{{:TYPE}}" class="form-control" placeholder="{{:PLACEHOLDER}}">

  {{#:UNITS}}{{{:UNITS}}}{{/:UNITS}}{{#:HINT}} 🎁{{/:HINT}}

"""

## Multiple choice (one of many)
## XXX add {{INLINE}}
html_templates["Radioq"] = mt"""
{{#:ITEMS}}

    <input class="form-check-input" type="radio" name="radio_{{:ID}}"
              id="radio_{{:ID}}_{{:VALUE}}" value="{{:VALUE}}">

        {{{:LABEL}}}

    </input>

{{/:ITEMS}}
"""

html_templates["radio_grading_script"] = """
document.querySelectorAll('input[name="radio_{{:ID}}"]').forEach(function(rb) {
rb.addEventListener("change", function() {
    var correct = rb.value == {{:CORRECT_ANSWER}};
    var msgBox = document.getElementById('{{:ID}}_message');
    $(grading_partial)
})});
"""
## ----

html_templates["Buttonq"] = mt"""

  {{#:BUTTONS}}
  <button class="toggle-btn" aria-pressed="false" id="button_{{:ID}}_{{:i}}" value="{{:ANSWER}}" style="width:100%;text-align:left; padding-left:10px; {{#:BLUE}}background:{{{:BLUE}}}{{/:BLUE}}" onclick="return false;">
    {{{:TEXT}}
  </button>
  {{/:BUTTONS}}

"""

html_templates["Multiq"] = mt"""
{{#:ITEMS}}

    <input class="form-check-input" type="checkbox" name="check_{{:ID}}"
              id="check_{{:ID}}_{{:VALUE}}" value="{{:VALUE}}">

        {{{:LABEL}}}

    </input>

{{/:ITEMS}}
"""

html_templates["multi_grading_script"] = """
document.querySelectorAll('input[name="check_{{:ID}}"]').forEach(function(rb) {
rb.addEventListener("change", function() {
    var choice_buttons = document.getElementsByName("check_{{:ID}}");
    var selected = [];
    for (var i=0; i < choice_buttons.length; i++) {
        if (choice_buttons[i].checked) {
           selected.push(i+1)
        }
    }
    var a = selected;
    var b = {{{:CORRECT_ANSWER}}};
    // https://stackoverflow.com/questions/7837456/how-to-compare-arrays-in-javascript
    var  correct =  (a.length === b.length && a.find((v,i) => v !== b[i]) === undefined)
    var msgBox = document.getElementById('{{:ID}}_message');
    $(grading_partial)
})});
"""

html_templates["MultiButtonq"] = mt"""

  {{#:BUTTONS}}
  <button class="toggle-btn" aria-pressed="false" id="button_{{:ID}}_{{:i}}" name="{{:i}}" value="unclicked" style="width:100%;text-align:left; padding-left:10px; {{#:BLUE}}background:{{{:BLUE}}}{{/:BLUE}}" onclick="return false;">
      {{{:TEXT}}
  </button>
  {{/:BUTTONS}}

  <button id="button_{{:ID}}-done"
     style="display:block; margin:auto;text-align:center;{{#:BLUE}}background:{{{:BLUE}}}{{/:BLUE}}" onclick="return false;">
      DONE
  </button>

"""

html_templates["multi_button_grading_script"] = """

// toggle button select
document.querySelectorAll('[id^="button_{{:ID}}_"]').forEach(function(btn) {
    btn.addEventListener("click", function(btn) {
	var unclicked = (this.value == "unclicked")
	if (unclicked) {
            this.style.background = "{{{:SELECTED_COLOR}}}";
	    this.value = "clicked";
            this.setAttribute("aria-pressed", "true");
	} else {
	    this.style.background = null;
	    this.value="unclicked";
            this.setAttribute("aria-pressed", "false");
	}
    });
})

// grade question
document.querySelector('[id^="button_{{:ID}}-done"]').addEventListener("click", function() {
                this.disabled = true;
    var multi_choice_buttons = document.querySelectorAll('[id^="button_{{:ID}}_"]');
    var selected = [];
    var correct = {{{:CORRECT_ANSWER}}};
    for (var i=0; i < multi_choice_buttons.length; i++) {
	var btn = multi_choice_buttons[i];
        btn.disabled = true;
	var btn_txt = btn.innerHTML;
	if (multi_choice_buttons[i].value == "clicked") {
	    selected.push(i+1)
	    if (correct.indexOf(i+1) < 0) {
		btn.innerHTML = "{{{:INCORRECT_flag}}}" + btn_txt
	    } else {
		btn.innerHTML =  "{{{:CORRECT_flag}}}" + btn_txt
                btn.style.fontSize = "1.1rem";
                btn.style.borderRadius = "12px";
            }
	} else {
	    if (correct.indexOf(i+1) < 0) {
		btn.innerHTML =  "{{{:INCORRECT_flag}}}" + btn_txt
	    } else {
		btn.innerHTML =  "{{{:CORRECT_flag}}}" + btn_txt
                btn.style.fontSize = "1.1rem";
                btn.style.borderRadius = "12px";
	    }
	}
      }
      var a = selected;
      var b = {{{:CORRECT_ANSWER}}}
          // https://stackoverflow.com/questions/7837456/how-to-compare-arrays-in-javascript
      var  correct =  (a.length === b.length && a.find((v,i) => v !== b[i]) === undefined)
      var msgBox = document.getElementById('{{:ID}}_message');
      $(grading_partial)
})
"""

## ----

html_templates["Matchq"] = mt"""

{{#:ITEMS}}

   {{{:QUESTION}}}

  <select name = "select_{{:ID}}" id="select_{{:ID}}">
     <option value=0 selected="selected">{{{:BLANK}}}</option>
     {{#:ANSWER_CHOICES}}
     <option value="{{:INDEX}}">{{{:LABEL}}}</option>
    {{/:ANSWER_CHOICES}}
  </select>

{{/:ITEMS}}

"""

html_templates["matchq_grading_script"] = """
      function callback(element, iterator) {
	  element.addEventListener("change", function() {
	      var a = [];
	      var selectors = document.querySelectorAll('[id ^= "select_{{:ID}}"]');
	      Array.prototype.forEach.call(selectors, function (element, iterator) {
		 a.push(element.value);
	      })
              var b = {{{:CORRECT_ANSWER}}};
              // https://stackoverflow.com/questions/7837456/how-to-compare-arrays-in-javascript
              var correct =  (a.length === b.length && a.find((v,i) => v !== b[i]) === undefined)
              var msgBox = document.getElementById('{{:ID}}_message');
              $(grading_partial)
	  })
      }
      Array.prototype.forEach.call(document.querySelectorAll('[id ^= "select_{{:ID}}"]'), callback);
"""

## ----

html_templates["fill_in_blank_select"] = """
  <select name = "select_{{:ID}}" id="select_{{:ID}}">
     <option value=0 selected="selected">{{{:BLANK}}}</option>
     {{#:ANSWER_CHOICES}}
     <option value="{{:INDEX}}">{{{:LABEL}}}</option>
    {{/:ANSWER_CHOICES}}
  </select>
"""

html_templates["fill_in_blank_input"] = """
<input id="{{:ID}}" type="{{:TYPE}}" class="form-control" placeholder="{{:PLACEHOLDER}}">
"""

## -------
html_templates["hotspot"] = """

"""

html_templates["hotspot_grading_script"] = """
  document.getElementById("hotspot_{{{:ID}}}").addEventListener("click", function(e) {
      var u = e.offsetX;
      var v = e.offsetY;
      var w = this.offsetWidth;
      var h = this.offsetHeight

      var x = u/w;
      var y = (h-v)/h

      var correct = {{{:CORRECT_ANSWER}}}
      var msgBox = document.getElementById('{{:ID}}_message');
      $(grading_partial)

  })
"""

## -----
## -------
html_templates["plotlylight_grading_script"] = """
document.getElementById("{{{:ID}}}").on("plotly_click", function(e) {
      var x = e.points[0].x
      var y = e.points[0].y

      var correct = {{{:CORRECT_ANSWER}}}
      var msgBox = document.getElementById('{{:ID}}_message');
      $(grading_partial)

  })
"""


text:
module TestRematch

using Compat.Test

@testset "Rematch.jl" begin
    include("rematch.jl")
end

end


text:
@testset "removing outlier observations" begin
    for i in 1:1000
        len = rand(10:80)
        x = Array{Union{Float64, Missing}}(rand(1:100.0,len))
        x[rand(1:len)] = missing
        @test length(outlier_ts(x)) == len
    end
end

@testset "generating outlier mask" begin
    for i in 1:10
        x = rand(10:30)
        y = rand(10:30)
        z = rand(10:30)
        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x, y, z))
        for j in 1:rand(1:20)
            rad[rand(1:x), rand(1:y), rand(1:z)] = missing
        end
        mask = rand(0:1, x, y)
        @test size(outlier_mask(rad, mask)) == (x, y)
    end
end


text:
module CellListsBenchmarks

export benchmark_cell_list_serial
export benchmark_cell_list_threads
export benchmark_brute_force
export benchmark_cell_lists
export benchmark_near_neighbors_serial
export benchmark_near_neighbors_threads
export benchmark_functions
export run_benchmark
export Benchmark

using Base.Sys
using Base.Threads
using Dates
using Random
using BenchmarkTools
using CellLists

function brute_force(p::Array{T, 2}, r::T) where T <: AbstractFloat
    ps = Vector{Tuple{Int, Int}}()
    n, d = size(p)
    for i in 1:(n-1)
        for j in (i+1):n
            if @inbounds distance_condition(p[i, :], p[j, :], r)
                push!(ps, (i, j))
            end
        end
    end
    return ps
end

function cell_lists(p::Array{T, 2}, r::T) where T <: AbstractFloat
    c = CellList(p, r)
    near_neighbors(c, p, r)
end

# --- Benchmark functions ---

function benchmark_cell_list_serial(p, r, seconds)
    @benchmark CellList($p, $r) seconds=seconds
end

function benchmark_cell_list_threads(p, r, seconds)
    @benchmark CellList($p, $r, $Val(:threads)) seconds=seconds
end

function benchmark_brute_force(p, r, seconds)
    @benchmark brute_force($p, $r) seconds=seconds
end

function benchmark_cell_lists(p, r, seconds)
    @benchmark cell_lists($p, $r) seconds=seconds
end

function benchmark_near_neighbors_serial(p, r, seconds)
    c = CellList(p, r)
    @benchmark near_neighbors($c, $p, $r) seconds=seconds
end

function benchmark_near_neighbors_threads(p, r, seconds)
    c = CellList(p, r)
    @benchmark near_neighbors($c, $p, $r, $Val(:threads)) seconds=seconds
end

const benchmark_functions = Dict(
    "cell_list_serial" => benchmark_cell_list_serial,
    "cell_list_threads" => benchmark_cell_list_threads,
    "brute_force" => benchmark_brute_force,
    "cell_lists" => benchmark_cell_lists,
    "near_neighbors_serial" => benchmark_near_neighbors_serial,
    "near_neighbors_threads" => benchmark_near_neighbors_threads,
)

# --- Run benchmark ---

function run_benchmark(benchmark::Function, n::Int, d::Int, r::T, seed::Int, iterations::Int, seconds::Float64) where T <: AbstractFloat
    @info "Function: benchmark_algorithm"
    @info "Arguments" benchmark n d r seed iterations seconds nthreads()
    rng = MersenneTwister(seed)
    trials = BenchmarkTools.Trial[]
    time = Time(now())
    @info "Started: $(time)"
    for i in 1:iterations
        @info "Iteration: $(i) / $(iterations)"
        p = rand(rng, T, n, d)
        t = benchmark(p, r, seconds)
        time2 = Time(now())
        @info "Time: $(Time(time2 - time))"
        time = time2
        push!(trials, t)
    end
    @info "Finished: $(Time(now()))"
    return trials
end

struct Benchmark{T<:AbstractFloat}
    benchmark::String
    n::Int
    d::Int
    r::T
    seed::Int
    trials::Vector{BenchmarkTools.Trial}
    nthreads::Int
    timestamp::DateTime
    version::VersionNumber
    cpu_info::Vector{Sys.CPUinfo}
end

function Benchmark(benchmark::String, n::Int, d::Int, r::T, seed::Int, iterations::Int, seconds::Float64) where T <: AbstractFloat
    trials = run_benchmark(benchmark_functions[benchmark], n, d, r, seed, iterations, seconds)
    Benchmark(benchmark, n, d, r, seed, trials, nthreads(), now(), VERSION, cpu_info())
end

end # module


text:
# To do a better job, one ought to check soundness relative to the
# BibTeX specification (http://www.bibtex.org/Format/): for now, this
# is just tested against a private .bib file, where it works well.
# Currently, one cannot e.g. have @string or @comment entries, though
# it could easily be included.

module SimpleBibTeX

export Citation, Bibliography, parsebibtex

import Base: show, getindex, haskey, get, length, iterate,    # overloads
             summary, ==

const MIN_SEPLEN = 10 # lead length for showing of `Citation`s

# ---- Types ----
struct Citation <: AbstractDict{String, String}
    kind::String
    key::String
    data::Dict{String, String}
end
Citation(kind, key) = Citation(titlecase(kind), key, Dict{String, String}()) # for 'kind', uppercase first letter, lowercase rest
# AbstractDict interface
getindex(C::Citation, key)       = C.data[key]
haskey(C::Citation, key::String) = haskey(C.data, key)
get(C::Citation, key, default)   = get(C.data, key, default)
length(C::Citation)              = length(C.data)
iterate(C::Citation)             = iterate(C.data)
iterate(C::Citation, i::Int)     = iterate(C.data, i)
# equivalence
(==)(C::Citation, C′::Citation)  = C.kind == C′.kind && C.key == C′.key && C.data == C′.data

struct Bibliography <: AbstractDict{String,Citation}
    citations::Dict{String, Citation}
end
Bibliography() = Bibliography(Dict{String, Citation}())
# AbstractDict interface
getindex(B::Bibliography, key)       = B.citations[key]
haskey(B::Bibliography, key::String) = haskey(B.citations, key)
get(B::Bibliography, key, default)   = get(B.citations, key, default)
length(B::Bibliography)              = length(B.citations)
iterate(B::Bibliography)             = iterate(B.citations)
iterate(B::Bibliography, i::Int)     = iterate(B.citations, i)

# ---- Show methods ----
delim(i, L) = i == 1 ? "┌" : i == L ? "└" : "│" # make connections pretty...

function show(io::IO, ::MIME"text/plain", C::Citation)
    println(io, "  ", C.key, " (", C.kind, ")")
    maxlen = !isempty(C.data) ? maximum(length.(keys(C))) : NaN
    seplen = max(MIN_SEPLEN, maxlen)
    entrystart = seplen + 4
    for (i,field) in enumerate(C.data)
        print(io, "  ", delim(i, length(C.data)), " $(field[1]) ",
                  "─"^(seplen-length(field[1])+1), " ")
        print(io, "$(field[2])\n")
    end
    return nothing
end

function summary(io::IO, B::Bibliography)
    n = length(B)
    Base.showarg(io, B, true)
    print(io, " with ", n, (n==1 ? " Citation" : " Citations"))
end
function show(io::IO, ::MIME"text/plain", B::Bibliography)
    summary(io, B)
    println(io, ':')
    i = 1
    for C in values(B)
        if i <= 5 # no need to print any more than at most 5 citations, after that, it's just too much
            println(io); show(io, MIME"text/plain"(), C);
            i+=1
        else
            print(io, "\n  … (+", length(B.citations)-5, " additional entries) …")
            break
        end
    end
    return nothing
end

# ---- Auxiliary, generic functions ----
"""
    _readuntil(s::IO, delims::Vector{T}) where T<:AbstractChar

Same as `readuntil(s::IO, delim::AbstractChar; ...)` but allows `delims` to be a vector or
tuple of characters: the first such match defines the termination condition.

Unlike its single `delim::AbstractChar` sibling, a two-tuple of values is returned:

    1. The read string, excluding the matching delimiter (`::String`).
    2. The matching delimiter (`::Char`). (`nothing` if no match)
"""
function _readuntil(s::IO, delims::Union{Vector{T}, Tuple{T, Vararg{T}}}) where T<:AbstractChar
    out = IOBuffer()
    while !eof(s)
        c = read(s, Char)
        if any(==(c), delims)
            return String(take!(out)), c
        end
        write(out, c)
    end
    return nothing
end

# ---- Functionality to parse a BibTeX file ----
"""
    readfield(io::IO)

Reads a field/preamble entry from a BibTeX file, assuming that `io::IO` (e.g. `IOStream`
or `IOBuffer`) has _just_ read `fieldkey =` or `@preamble`, (where `fieldkey` is e.g.
`author`, `title`, etc). In that case, the field/preamble contents are enclosed in
a set of {...} brances. We do a small dance to allow braces within those braces.
The output is the contents within those braces. The `IOStream`/`IOBuffer` is moved to the
end of the matching braces.
"""
function readfield(io::IO, ispreamble::Bool=false)
    ispreamble || readuntil(io, '{') # move "read-cursor" up to start of contents
    bracesum = 1 # brace counter: `{` adds 1, `}` subtracts 1; "done" when `bracesum = 0`
    fieldentry_io = IOBuffer()
    # start reading contents (loop, because we need to allow for braces within contents)
    while !eof(io)
        r = _readuntil(io, ('{', '}'))
        r !== nothing ? (step=first(r); match_delim=last(r)) : throw("Malformed field content")

        if length(step) == 0 || last(step) ≠ '\\'
            bracesum += match_delim == '}' ? -1 : 1
        end

        if iszero(bracesum) # contents are terminated, cf. brace count; return field string
            write(fieldentry_io, step)
            return String(take!(fieldentry_io))
        else
            write(fieldentry_io, step, match_delim) # "add" what we have read so far into field contents
        end
    end

    throw("Malformed field; unbalanced parenthesis")
end

"""
    parsebibtex(io::IO)

Builds a `Bibliography`, containing `Citation` fields, corresponding
to the BibTeX entries in the `io::IO` object (e.g. `IOStream` or `IOBuffer`).
Also outputs a preamble as a `String` (`""` if empty or absent).
"""
function parsebibtex(io::IO)
    B = Bibliography()
    preamble = ""
    while !eof(io)
        if last(readuntil(io, '@', keep = true)) == '@' # @ marks new citation
            kind = readuntil(io, '{') # type of citation (e.g., Article, Book, etc.)
            if lowercase(kind) == "preamble" # BibTeX entry of the `preamble` kind
                preamble = readfield(io, true)

            else
                # assume this is an "ordinary" BibTeX entry (e.g., kind is not `preamble`,
                # `comment`, or `string`)
                key  = readuntil(io, ',') # key of citation  (e.g., John:2012)
                C = Citation(kind, key)
                # find fields in citation
                while !eof(io)
                    # some citations may have no fields; check for `}` in that case;
                    # if it has any fields, their key terminates with `=`
                    r = _readuntil(io, ('=', '}'))
                    r !== nothing ? (fieldkey=first(r); match_delim=last(r)) : throw("Malformed or unbalanced field key")
                    fieldkey = lowercase(fieldkey)

                    match_delim == '}' && break # no fields in citation; go to next one
                    #otherwise: a field was found; match_delim is '='

                    # remove comma from last line & do some simple clean-up, in case user
                    # entered quasi-garbled format. `fieldkey` is e.g., author, title, etc.
                    fieldkey = replace(fieldkey, r"\n|\r| |\t|,"=>"")

                    # read the field contents (delimited by {...} braces) and clean up in
                    # case user gave quasi-garbled format for field contents
                    fieldentry = replace(readfield(io), r"\n|\r|\t|[ ]{2,}"=>' ')
                    C.data[fieldkey] = fieldentry # add field to the Citation's dict entry
                end
                B.citations[key] = C # add the Citation to the Bibliography
            end
        end
    end

    return B, preamble
end

"""
    parsebibtex(filename_or_string::String; isfilename::Bool=true)

Constructs a `Bibliography` struct, containing `Citation` fields, corresponding
to the BibTeX entries in the `filename_or_string` string.
If `isfilename=true` (default), the input string is interpreted as the
location of a file. If `false`, the string is interpreted and parsed as BibTeX.
Also outputs a preamble as a `String` (`""` if empty or absent).
"""
function parsebibtex(filename_or_string::String; isfilename::Bool=true)
    # treat string-input as the name of a file, relative to caller
    if isfilename == true
        open(filename_or_string) do io
            return parsebibtex(io)
        end

    # treat input as a valid BibTeX entry
    else
        return parsebibtex(IOBuffer(filename_or_string))
    end
end

end # module


text:
include("quadratic.jl")

println(quadratic(x -> 2x^2 + 30x + 9))
# results:
# 1st approximation:
# -0.3062524369827744 + 0.0im,-14.693907564071102 - 0.0im)
# with eps:
# -0.30625271503092044 + 0.0im,-14.69374728496908 - 0


text:
module MuZero

export Mu, MuParams, MuNetworkHP,
  PredictionSimpleHP, DynamicsSimpleHP, RepresentationSimpleHP,
  PredictionResnetHP, DynamicsResnetHP, RepresentationResnetHP,
  MuEnv, mutrain!, self_play_step!, learning_step!, MuNetwork,
  run_duel

# 3x Faster than julia 1.6.2 default implementation #! may break something
Base.hash(x::Array{Float32,3}) = hash(vec(x))
Base.hash(x::Tuple{Array{Float32,3},Int}) = hash(x[1]) >> x[2]

using AlphaZero
using ProgressLogging
using TensorBoardLogger
using Logging
using ParameterSchedulers: Scheduler, Cos
import Flux
import FileIO
import CUDA

import LinearAlgebra
LinearAlgebra.BLAS.set_num_threads(1)

include("mu_game_wrapper.jl")
include("network.jl")
# include("../alphazerolike.jl")
include("trace.jl")
include("play.jl")
include("simulations.jl")
include("training.jl")
include("learning.jl")
include("benchmark.jl")
include("probe-games.jl")

include("ttt_example.jl")
import .TTTexample

end # module


text:
"""
	correctinterval!(datain,corrpar)
Correct time intervals either setting them to NaNs or interpolate value.
In addition, steps in input time series can by corrected.

**Input**
* datain: DataFrame where at least one column contains DateTime
* corrpar: DataFrame with correction parameters or name of a file containing these information
Should contain following parameters (keys):
column = number of a data column or name/key (if number & :datetime in first column: ++column)
id = type of correction: 1=step, 2=set interval to NaN, 3=set interval to interpolated values
5=set interval to given values (range between y1 and y2) (id=4 not yet implemented).
x1,x2 = starting,end time of the interval
y1,y2 = step correction: value at before at after step (id=1) | values used to replace interval (see id=5). If id=0 and y1==0.0 && isnan(y2), than difference is set to zero
* includetime: switch to either to apply correction to time including x1/x2 or just > & < (true: => & <=)

**Output**
* corrected DataFrames. Call the function without ! to get a corrected copy of the input dataframe

**Example**
```
datain = DataFrame(
		   datetime=collect(DateTime(2010,1,1,3):Dates.Minute(30):DateTime(2010,1,2,12)),
		   grav = zeros(Float64,67) .+ 9.8,
		   pres = zeros(Float64,67) .+ 1000.,
		   temp = zeros(Float64,67) .+ 27.);
corrfile = "f:/mikolaj/code/libraries/matlab_octave_library/test/input/correctTimeInterval_inputFile.txt"
dataout = correctinterval(datain,corrfile);
# Set parameters directly and overwrite the input dataframe
corrpar = DataFrame(column=[1,2,3], id = [3,2,1],
					x1 = [DateTime(2010,01,01,04,30,00),
						  DateTime(2010,01,01,08,00,00),
						  DateTime(2010,01,02,04,00,00)],
				  	x2 = [DateTime(2010,01,01,07,30,00),
						  DateTime(2010,01,01,09,30,09),
						  DateTime(2010,01,02,06,30,00)],
					y1 = [NaN,NaN,10.],y2 = [NaN,NaN,0.0]);
correctinterval!(datain,corrpar);
```
"""
function correctinterval!(datain::DataFrame,par::DataFrame;includetime::Bool=true)
	corrpar = deepcopy(par);
	if names(datain)[1] == :datetime && eltype(corrpar[:column]) != Symbol
		corrpar[:column] = corrpar[:column] .+ 1;
	end
	for (i,v) in enumerate(corrpar[:id])
		if eltype(corrpar[:column])!=Symbol ? corrpar[:column][i]<=size(datain,2) : haskey(datain,corrpar[:column][i])
			if v == 1
				correctinterval_step!(datain,corrpar,i,includetime);
			elseif v == 2
				correctinterval_nan!(datain,corrpar,i,includetime);
			elseif v == 3
				correctinterval_interp!(datain,corrpar,i,includetime);
			elseif v == 5
				correctinterval_replace!(datain,corrpar,i,includetime);
			end
		end
	end
	return datain
end
function correctinterval(datain::DataFrame,corrpar::DataFrame;includetime::Bool=true)
	correctinterval!(deepcopy(datain),corrpar,includetime=includetime)
end
function correctinterval(datain::DataFrame,corrfile::String;includetime::Bool=true)
	correctinterval!(deepcopy(datain),FileTools.readcorrpar(corrfile),includetime=includetime)
end
function correctinterval!(datain::DataFrame,corrfile::String;includetime::Bool=true)
	correctinterval!(datain,FileTools.readcorrpar(corrfile),includetime=includetime)
end

"""
Auxiliary function to correct step in time series
"""
function correctinterval_step!(datain,corrpar,i,includetime)
	# find points recorded after the step occur.
	r = includetime ? findall(x->x .>= corrpar[:x2][i], datain[:datetime]) :
					  findall(x->x .> corrpar[:x2][i], datain[:datetime]);
	# check if both values are given (if second not, then set diff. to 0)
	if !isnan(corrpar[:y1][i]) && isnan(corrpar[:y2][i])
		r1 = includetime ? findall(x->x .>= corrpar[:x1][i], datain[:datetime]) :
						  findall(x->x .> corrpar[:x1][i], datain[:datetime]);
		applyDiff = datain[corrpar[:column][i]][r[1]] -
						datain[corrpar[:column][i]][r1[1]] -
							corrpar[:y1][i];
	else
		applyDiff = corrpar[:y2][i] - corrpar[:y1][i];
	end
	for j in r # remove the step by SUBTRACTING the given difference.
        datain[corrpar[:column][i]][j] = datain[corrpar[:column][i]][j] -
						applyDiff;
    end
end

"""
Auxiliary function to set interval to NaNs
"""
function correctinterval_nan!(datain,corrpar,i,includetime)
	# find points recorded in-between given time epochs.
	r = findinterval(datain[:datetime],corrpar[:x1][i],corrpar[:x2][i],
						includetime=includetime)
	for j in r
		datain[corrpar[:column][i]][j] = NaN
	end
end

"""
Auxiliary function to set interval to linearly interpolated values
"""
function correctinterval_interp!(datain,corrpar,i,includetime)
	# find points recorded withing interval.
    r = includetime ? map(x->x .>= corrpar[:x1][i] && x .<= corrpar[:x2][i],datain[:datetime]) :
					  map(x->x .> corrpar[:x1][i] && x .< corrpar[:x2][i],datain[:datetime]);
    if any(r)
		# get values except affected interval and use it in interpolation
		ytemp = datain[corrpar[:column][i]][.!r];
		xtemp = Dates.value.(datain[:datetime][.!r]);
		datain[corrpar[:column][i]][r] = interp1(xtemp,ytemp,Dates.value.(datain[:datetime][r]));
    end
end

"""
Auxiliary function to replace interval using given value
"""
function correctinterval_replace!(datain,corrpar,i,includetime)
	# find points recorded in-between given time epochs.
	r = findinterval(datain[:datetime],corrpar[:x1][i],corrpar[:x2][i],
						includetime=includetime)
	rep_val = range(corrpar[:y1],stop=corrpar[:y2],length=length(r));
	for (j,v) in enumerate(r)
		datain[corrpar[:column][i]][v] = rep_val[j][1];
	end
end

"""
Auxiliary function to find interval between two points
"""
function findinterval(timein,x1,x2;includetime::Bool=true)
	r = includetime ? findall(x->x .>= x1 && x .<= x2,timein) :
					  findall(x->x .> x1 && x .< x2,timein);
end

"""
	prepcorrpar(datain,timein;min_gap,defcol,defid)

Function to prepare correction parameter used in 'correctinterval'.
Will find all blocks of NaNs and create formated DataFrame (that can be then
manually updated)

**Input**
* datain: input vector to be examined. Regular sampling is assumed! (use time2regular)
* timein: input DateTime vector corresponding to datain
* min_gap: minimum gap/time span of NaNs (in index, not time!). See example
* defcol: default column number or symbol (see correctinterval function)
* defid: default correction ID (see correctinterval function)

**Output**
* correction parameter DataFrame. Empty if no NaNs found.

**Example**
```
dfin = DataFrame(Temp = collect(0.:1.:12.),
	 datetime= collect(DateTime(2000,1,1):Dates.Hour(1):DateTime(2000,1,1,12)))
dfin[:Temp][[6,10,11]] .= NaN;
min_gap = 2; # corresponds to 2 hours (see datetime column)
corrpar = prepcorrpar(dfin[:Temp],dfin[:datetime],min_gap=min_gap,
					defcol=:Temp,defid=2);
```
"""
function prepcorrpar(datain,timein::Vector{DateTime};
						min_gap::Int=1,defcol=1,defid::Int=1)
	corrpar = DataFrame();
	nstart,nstop = ResampleAndFit.findnanblocks(datain);
	outlength = length(nstart);
	if !isempty(nstart)
		corrpar = DataFrame(column = repeat([defcol],outlength),
				id = repeat([defid],outlength),
				x1 = timein[nstart],
				x2 = timein[nstop],
				y1 = zeros(outlength).+NaN,
				y2 = zeros(outlength).+NaN;
				comment = repeat(["automatically_generated_using_preparecorrpar"],outlength));
		# remove rows where gap is too short
		deleterows!(corrpar,prepcorrpar_remrow(nstart,nstop,min_gap))
		corrpar = size(corrpar,1) == 0 ? DataFrame() : corrpar; # return empty if all columns deleted
	end
	return corrpar
end

"""
Auxiliary function to find indices where the gap is < as given value
"""
function prepcorrpar_remrow(nstart,nstop,min_gap)
	remrow = Vector{Int}();
	for i in 1:length(nstart)
		(nstop[i] - nstart[i] + 1) < min_gap ? push!(remrow,i) : nothing
	end
	return remrow
end


text:
using Base.Threads, LinearAlgebra
import Base: ndims, getindex, materialize!, broadcasted
export ScalarField

export ndims

abstract type AbstractField end
struct ScalarField{S,T} <: AbstractField
    data::S
    grid::T
end

function (ϕ::ScalarField)(x::Tuple)
    return getvalue(ϕ.data, x, ϕ.grid)
end

function (ϕ::ScalarField)(x::Number, y::Number, z::Number)
    return getvalue(ϕ.data, (x,y,z), ϕ.grid)
end

function (ϕ::ScalarField)(x::Number, y::Number)
    return getvalue(ϕ.data, (x,y), ϕ.grid)
end

getindex(ϕ::ScalarField, i::Int) = ϕ.data[i]

materialize!(ϕ::ScalarField, f::Base.Broadcast.Broadcasted) = materialize!(ϕ.data, f)
broadcasted(identity, ϕ::ScalarField) = broadcasted(Base.identity, ϕ.data)

function (ϕ::ScalarField)(xlist::StepRangeLen,ylist::StepRangeLen,zlist::StepRangeLen; threads = false)
    newfield = zeros(length(xlist), length(ylist), length(zlist))
    if threads
        @threads for k in eachindex(zlist)
            for j in eachindex(ylist)
                for i in eachindex(xlist)
                        newfield[i,j,k] = getvalue(ϕ.data, (xlist[i],ylist[j],zlist[k]), ϕ.grid)
                end
            end
        end
    else
        for k in eachindex(zlist)
            for j in eachindex(ylist)
                for i in eachindex(xlist)
                        newfield[i,j,k] = getvalue(ϕ.data, (xlist[i],ylist[j],zlist[k]), ϕ.grid)
                end
            end
        end
    end
    return newfield
end

function (ϕ::ScalarField)(xlist::StepRangeLen,ylist::StepRangeLen; threads = false)
    newfield = zeros(length(xlist), length(ylist))
    if threads
    @threads for j in eachindex(ylist)
            for i in eachindex(xlist)
                newfield[i,j] = getvalue(ϕ.data, (xlist[i],ylist[j]), ϕ.grid)
            end
        end
    else
        for j in eachindex(ylist)
            for i in eachindex(xlist)
                    newfield[i,j] = getvalue(ϕ.data, (xlist[i],ylist[j]), ϕ.grid)
            end
        end
    end
    return newfield
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
"""
    assemble_matrices(model)

Assemble finite element matrices.
"""
function assemble_matrices(model)
    (; mesh, D, T₂, ρ) = model

    # Deduce sizes
    ncompartment, nboundary = size(mesh.facets)
    npoint_cmpts = size.(mesh.points, 2)

    # Assemble finite element matrices
    volumes = zeros(ncompartment)
    for icmpt = 1:ncompartment
        points = mesh.points[icmpt]
        elements = mesh.elements[icmpt]
        fevolumes, = get_mesh_volumes(points, elements)
        volumes[icmpt] = sum(fevolumes)
    end

    # Assemble finite element matrices compartment-wise
    M_cmpts = []
    S_cmpts = []
    Mx_cmpts = [[] for _ = 1:3]
    G = []
    for icmpt = 1:ncompartment
        # Finite elements
        points = mesh.points[icmpt]
        elements = mesh.elements[icmpt]
        fevolumes, = get_mesh_volumes(points, elements)

        # Assemble mass, stiffness and flux matrices
        push!(M_cmpts, assemble_mass_matrix(elements', fevolumes))
        push!(S_cmpts, assemble_stiffness_matrix(elements', points', D[icmpt]))

        # Assemble first order product moment matrices
        for dim = 1:3
            push!(Mx_cmpts[dim], assemble_mass_matrix(elements', fevolumes, points[dim, :]))
        end

        # Compute surface integrals
        push!(G, zeros(npoint_cmpts[icmpt], 3))
        for iboundary = 1:nboundary
            facets = mesh.facets[icmpt, iboundary]
            if !isempty(facets)
                # Get facet normals
                _, _, normals = get_mesh_surfacenormals(points, elements, facets)

                # Surface normal weighted flux matrix (in each canonical direction)
                for dim = 1:3
                    Q = assemble_flux_matrix(facets', points', normals[dim, :])
                    G[icmpt][:, dim] += sum(Q, dims = 2)
                end
            end
        end

    end

    # Assemble global finite element matrices
    M = blockdiag(M_cmpts...)
    S = blockdiag(S_cmpts...)
    R = blockdiag((M_cmpts ./ T₂)...)
    Mx = [blockdiag(Mx_cmpts[dim]...) for dim = 1:3]
    Q_blocks = assemble_flux_matrices(mesh.points, mesh.facets)
    Q = couple_flux_matrix(model, Q_blocks, false)

    (; M, S, R, Mx, Q, M_cmpts, S_cmpts, Mx_cmpts, G, volumes)
end


text:
d1 = Dict(1 => 4.2, 2 => 5.3)
# Dict{Int64,Float64} with 2 entries:
#   2 => 5.3
#   1 => 4.2
d1 = Dict{Int64,Float64}(1 => 4.2, 2 => 5.3)
d1 = [1 => 4.2, 2 => 5.3]
# 2-element Array{Pair{Int64,Float64},1}:
#  1 => 4.2
#  2 => 5.3

d2 = Dict{Any,Any}("a"=>1, (2,3)=>true)

d3 = Dict(:A => 100, :B => 200)
# Dict{Symbol,Int64} with 2 entries:
#   :A => 100
#   :B => 200
d3 = Dict{Symbol,Int64}(:A => 100, :B => 200)
d3[:B] #> 200
# d3[:Z] #> ERROR: KeyError: key :Z not found
get(d3, :Z, 999) #> 999

d3[:A] = 150 #> d3 is now [:A => 150, :B => 200]
d3[:C] = 300 #> d3 is now [:A => 150, :B => 200, :C => 300]
length(d3) #> 3
# d3["CVO"] = 500 #> ERROR: KeyError: key "CVO" not found
# d3[:CVO] = "Julia" #> ERROR: KeyError: key "CVO" not found

dmus = Dict{Symbol,String}(:first_name => "Louis", :surname => "Armstrong", :occupation => "musician",
 		 :date_of_birth => "4/8/1901")
# Dict{Symbol,String} with 4 entries:
#   :date_of_birth => "4/8/1901"
#   :occupation    => "musician"
#   :surname       => "Armstrong"
#   :first_name    => "Louis"

haskey(d3, :Z) #> false
haskey(d3, :B) #> true

d4 = Dict() #> Dict{Any,Any} with 0 entries
d4["lang"] = "Julia"

d5 = Dict{Float64, Int64}() #> Dict{Float64,Int64} with 0 entries
# d5["c"] = 6 #> ERROR: MethodError: `convert` has no method matching convert(::Type{Float64}, ::ASCIIString)

d3 = Dict(:A => 100, :B => 200)
delete!(d3, :B) #> Dict{Symbol,Int64} with 1 entry:  :A => 100

d3 = Dict(:A => 100, :B => 200)
ki = keys(d3) #> KeyIterator
for k in ki
   println(k)
end #> A   B
:A in ki #> true
:Z in ki #> false

collect(keys(d3)) #> 2-element Array{Symbol,1}:
#  :A
#  :B
vi = values(d3) #> ValueIterator
for v in values(d3)
    println(v)
end #> 100  200

keys1 = ["J.S. Bach", "Woody Allen", "Barack Obama"]
values1 = [ 1685, 1935, 1961]
d5 = Dict(zip(keys1, values1))
#>
# Dict{String,Int64} with 3 entries:
#   "J.S. Bach"    => 1685
#   "Woody Allen"  => 1935
#   "Barack Obama" => 1961

for (k, v) in d5
  println("$k was born in $v")
end
for p in d5
  println("$(p[1]) was born in $(p[2])")
end
#J.S. Bach was born in 1685
#Barack Obama was born in 1961
#Woody Allen was born in 1935

capitals = Dict{String, String}("France"=> "Paris", "China"=>"Beijing")
# Dict{String,String} with 2 entries:
#   "China"  => "Beijing"
#   "France" => "Paris"

# neat tricks:
dict = Dict("a" => 1, "b" => 2, "c" => 3)
arrkey = [key for (key, value) in dict] #> 3-element Array{String,1}: "a" "b" "c"
# same as collect(keys(dict))
arrval = [value for (key, value) in dict] #> 3-element Array{Int64,1}: 1 2 3
# same as collect(values(dict))


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using BenchmarkTools
using Random
function mcpi_f32(n)
    rng = MersenneTwister()
    4count(_ -> rand(rng, Float32)^2 + rand(rng, Float32)^2 ≤ 1, 1:n)/n
end
@show mcpi_f32(10^8)
@btime mcpi_f32(10^8);

# %%
using BenchmarkTools
using CUDA
mcpi_f32_cu(n) = 4count(x -> x^2 + rand(Float32)^2 ≤ 1, CUDA.rand(n))/n
@show mcpi_f32_cu(10^8)
@btime mcpi_f32_cu(10^8);

# %%
using BenchmarkTools
using Random
function mcpi_f64(n)
    rng = MersenneTwister()
    4count(_ -> rand(rng)^2 + rand(rng)^2 ≤ 1, 1:n)/n
end
@show mcpi_f64(10^8)
@btime mcpi_f64(10^8);

# %%
using BenchmarkTools
using CUDA
mcpi_f64_cu(n) = 4count(x -> x^2 + rand(Float64)^2 ≤ 1, CUDA.rand(Float64, n))/n
@show mcpi_f64_cu(10^8)
@btime mcpi_f64_cu(10^8);

# %%
CUDA.versioninfo()

# %%

max_stars_repo_path: 0017/CUDA mcpi.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0017/CUDA mcpi.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0017/CUDA mcpi.jl

max_forks_repo_name: genkuroki/public


text:
module DynamicPricingExamples

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# This file implements the Semi-NCA (SNCA) dominator tree construction
# described in Georgiadis' PhD thesis [LG05], which itself is a simplification
# of the Simple Lenguare-Tarjan (SLT) algorithm [LG79]. This algorithm matches
# the algorithm choice in LLVM and seems to be a sweet spot in implementation
# simplicity and efficiency.
#
# This file also implements an extension of SNCA that supports updating the
# dominator tree with insertion and deletion of edges in the control flow
# graph, described in [GI16] as Dynamic SNCA. DSNCA was chosen over DBS, a
# different algorithm which achieves the best overall performance in [GI16],
# because it is simpler to understand and implement, performs well with edge
# deletions, and is of similar performance overall.
#
# SNCA works by first computing semidominators, then computing immediate
# dominators from them. The semidominator of a node is the node with minimum
# preorder number such that there is a semidominator path from it to the node.
# A semidominator path is a path in which the preorder numbers of all nodes not
# at the endpoints are greater than the preorder number of the last node.
# Intuitively, the semidominator approximates the immediate dominator of a node
# by taking the path (in the CFG) that gets as close to the root as possible
# while avoiding ancestors of the node in the DFS tree.
#
# In computing the semidominators, SNCA performs "path compression" whenever a
# node has a nontrivial semidominator (i.e. a semidominator that is not just
# its parent in the DFS tree). Path compression propagates the "label" of a
# node, which represents a possible semidominator with associated semidominator
# path passing through that node.
#
# For example, path compression will be performed for the following CFG, where
# the edge not in the DFS tree is marked with asterisks. Note that nodes are
# labeled with their preorder numbers, and all edges point downward.
#
#     1
#     |\
#     | \
#     |  4
#     |  |
#     2  5
#     |  |
#     |  6
#     | *
#     |*
#     3
#
# There is a nontrivial semidominator path from 1 to 3, passing through 4, 5,
# and 6. Stepping through the whole algorithm on paper with an example like
# this is very helpful for understanding how it works.
#
# DSNCA runs the whole algorithm from scratch if the DFS tree is invalidated by
# the insertion or deletion, but otherwise recomputes a subset of the
# semidominators (all immediate dominators then need to be recomputed).
#
# [LG05]  Linear-Time Algorithms for Dominators and Related Problems
#         Loukas Georgiadis, Princeton University, November 2005, pp. 21-23:
#         ftp://ftp.cs.princeton.edu/reports/2005/737.pdf
#
# [LT79]  A fast algorithm for finding dominators in a flowgraph
#         Thomas Lengauer, Robert Endre Tarjan, July 1979, ACM TOPLAS 1-1
#         http://www.dtic.mil/dtic/tr/fulltext/u2/a054144.pdf
#
# [GI16]  An Experimental Study of Dynamic Dominators
#         Loukas Georgiadis, Giuseppe F. Italiano, Luigi Laura, Federico
#         Santaroni, April 2016
#         https://arxiv.org/abs/1604.02711

# We could make these real structs, but probably not worth the extra
# overhead. Still, give them names for documentary purposes.
const BBNumber = Int
const PreNumber = Int
const PostNumber = Int

struct DFSTree
    # These map between BB number and pre- or postorder numbers
    to_pre::Vector{PreNumber}
    from_pre::Vector{BBNumber}
    to_post::Vector{PostNumber}
    from_post::Vector{BBNumber}

    # Records parent relationships in the DFS tree
    # (preorder number -> preorder number)
    # Storing it this way saves a few lookups in the snca_compress! algorithm
    to_parent_pre::Vector{PreNumber}
end

function DFSTree(n_blocks::Int)
    return DFSTree(zeros(PreNumber, n_blocks),
                   Vector{BBNumber}(undef, n_blocks),
                   zeros(PostNumber, n_blocks),
                   Vector{BBNumber}(undef, n_blocks),
                   zeros(PreNumber, n_blocks))
end

copy(D::DFSTree) = DFSTree(copy(D.to_pre),
                           copy(D.from_pre),
                           copy(D.to_post),
                           copy(D.from_post),
                           copy(D.to_parent_pre))

function copy!(dst::DFSTree, src::DFSTree)
    copy!(dst.to_pre, src.to_pre)
    copy!(dst.from_pre, src.from_pre)
    copy!(dst.to_post, src.to_post)
    copy!(dst.from_post, src.from_post)
    copy!(dst.to_parent_pre, src.to_parent_pre)
    return dst
end

length(D::DFSTree) = length(D.from_pre)

function DFS!(D::DFSTree, blocks::Vector{BasicBlock})
    copy!(D, DFSTree(length(blocks)))
    to_visit = Tuple{BBNumber, PreNumber, Bool}[(1, 0, false)]
    pre_num = 1
    post_num = 1
    while !isempty(to_visit)
        # Because we want the postorder number as well as the preorder number,
        # we don't pop the current node from the stack until we're moving up
        # the tree
        (current_node_bb, parent_pre, pushed_children) = to_visit[end]

        if pushed_children
            # Going up the DFS tree, so all we need to do is record the
            # postorder number, then move on
            D.to_post[current_node_bb] = post_num
            D.from_post[post_num] = current_node_bb
            post_num += 1
            pop!(to_visit)

        elseif D.to_pre[current_node_bb] != 0
            # Node has already been visited, move on
            pop!(to_visit)
            continue
        else
            # Going down the DFS tree

            # Record preorder number
            D.to_pre[current_node_bb] = pre_num
            D.from_pre[pre_num] = current_node_bb
            D.to_parent_pre[pre_num] = parent_pre

            # Record that children (will) have been pushed
            to_visit[end] = (current_node_bb, parent_pre, true)

            # Push children to the stack
            for succ_bb in blocks[current_node_bb].succs
                push!(to_visit, (succ_bb, pre_num, false))
            end

            pre_num += 1
        end
    end

    # If all blocks are reachable, this is a no-op, otherwise, we shrink these
    # arrays.
    resize!(D.from_pre, pre_num - 1)
    resize!(D.from_post, post_num - 1) # should be same size as pre_num - 1
    resize!(D.to_parent_pre, pre_num - 1)

    return D
end

DFS(blocks::Vector{BasicBlock}) = DFS!(DFSTree(0), blocks)

"""
Keeps the per-BB state of the Semi NCA algorithm. In the original formulation,
there are three separate length `n` arrays, `label`, `semi` and `ancestor`.
Instead, for efficiency, we use one array in a array-of-structs style setup.
"""
struct SNCAData
    semi::PreNumber
    label::PreNumber
end

"Represents a Basic Block, in the DomTree"
struct DomTreeNode
    # How deep we are in the DomTree
    level::Int
    # The BB indices in the CFG for all Basic Blocks we immediately dominate
    children::Vector{BBNumber}
end

DomTreeNode() = DomTreeNode(1, Vector{BBNumber}())

"Data structure that encodes which basic block dominates which."
struct DomTree
    # These can be reused when updating domtree dynamically
    dfs_tree::DFSTree
    snca_state::Vector{SNCAData}

    # Which basic block immediately dominates each basic block, using BB indices
    idoms_bb::Vector{BBNumber}

    # The nodes in the tree (ordered by BB indices)
    nodes::Vector{DomTreeNode}
end

function DomTree()
    return DomTree(DFSTree(0), SNCAData[], BBNumber[], DomTreeNode[])
end

function construct_domtree(blocks::Vector{BasicBlock})
    return update_domtree!(blocks, DomTree(), true, 0)
end

function update_domtree!(blocks::Vector{BasicBlock}, domtree::DomTree,
                         recompute_dfs::Bool, max_pre::PreNumber)
    if recompute_dfs
        DFS!(domtree.dfs_tree, blocks)
    end

    if max_pre == 0
        max_pre = length(domtree.dfs_tree)
    end

    SNCA!(domtree, blocks, max_pre)
    compute_domtree_nodes!(domtree)
    return domtree
end

function compute_domtree_nodes!(domtree::DomTree)
    # Compute children
    copy!(domtree.nodes,
          DomTreeNode[DomTreeNode() for _ in 1:length(domtree.idoms_bb)])
    for (idx, idom) in Iterators.enumerate(domtree.idoms_bb)
        (idx == 1 || idom == 0) && continue
        push!(domtree.nodes[idom].children, idx)
    end
    # Recursively set level
    update_level!(domtree.nodes, 1, 1)
    return domtree.nodes
end

function update_level!(nodes::Vector{DomTreeNode}, node::BBNumber, level::Int)
    worklist = Tuple{BBNumber, Int}[(node, level)]
    while !isempty(worklist)
        (node, level) = pop!(worklist)
        nodes[node] = DomTreeNode(level, nodes[node].children)
        foreach(nodes[node].children) do child
            push!(worklist, (child, level+1))
        end
    end
end

"""
The main Semi-NCA algrithm. Matches Figure 2.8 in [LG05]. Note that the
pseudocode in [LG05] is not entirely accurate. The best way to understand
what's happening is to read [LT79], then the description of SLT in [LG05]
(warning: inconsistent notation), then the description of Semi-NCA.
"""
function SNCA!(domtree::DomTree, blocks::Vector{BasicBlock}, max_pre::PreNumber)
    D = domtree.dfs_tree
    state = domtree.snca_state
    # There may be more blocks than are reachable in the DFS / dominator tree
    n_blocks = length(blocks)
    n_nodes = length(D)

    # `label` is initialized to the identity mapping (though the paper doesn't
    # make that clear). The rationale for this is Lemma 2.4 in [LG05] (i.e.
    # Theorem 4 in [LT79]). Note however, that we don't ever look at `semi`
    # until it is fully initialized, so we could leave it uninitialized here if
    # we wanted to.
    resize!(state, n_nodes)
    for w in 1:max_pre
        # Only reset semidominators for nodes we want to recompute
        state[w] = SNCAData(typemax(PreNumber), w)
    end

    # If we are only recomputing some of the semidominators, the remaining
    # labels should be reset, because they may have become inapplicable to the
    # node/semidominator we are currently processing/recomputing. They can
    # become inapplicable because of path compressions that were triggered by
    # nodes that should only be processed after the current one (but were
    # processed the last time `SNCA!` was run).
    #
    # So, for every node that is not being reprocessed, we reset its label to
    # its semidominator, which is the value that its label assumes once its
    # semidominator is computed. If this was too conservative, i.e. if the
    # label would have been updated before we process the current node in a
    # situation where all semidominators were recomputed, then path compression
    # will produce the correct label.
    for w in max_pre+1:n_nodes
        semi = state[w].semi
        state[w] = SNCAData(semi, semi)
    end

    # Calculate semidominators, but only for blocks with preorder number up to
    # max_pre
    ancestors = copy(D.to_parent_pre)
    for w::PreNumber in reverse(2:max_pre)
        # LLVM initializes this to the parent, the paper initializes this to
        # `w`, but it doesn't really matter (the parent is a predecessor, so at
        # worst we'll discover it below). Save a memory reference here.
        semi_w = typemax(PreNumber)
        last_linked = PreNumber(w + 1)
        for v ∈ blocks[D.from_pre[w]].preds
            # For the purpose of the domtree, ignore virtual predecessors into
            # catch blocks.
            v == 0 && continue

            v_pre = D.to_pre[v]

            # Ignore unreachable predecessors
            v_pre == 0 && continue

            # N.B.: This conditional is missing from the pseudocode in figure
            # 2.8 of [LG05]. It corresponds to the `ancestor[v] != 0` check in
            # the `eval` implementation in figure 2.6
            if v_pre >= last_linked
                # `v` has already been processed, so perform path compression

                # For performance, if the number of ancestors is small avoid
                # the extra allocation of the worklist.
                if length(ancestors) <= 32
                    snca_compress!(state, ancestors, v_pre, last_linked)
                else
                    snca_compress_worklist!(state, ancestors, v_pre, last_linked)
                end
            end

            # The (preorder number of the) semidominator of a block is the
            # minimum over the labels of its predecessors
            semi_w = min(semi_w, state[v_pre].label)
        end
        state[w] = SNCAData(semi_w, semi_w)
    end

    # Compute immediate dominators, which for a node must be the nearest common
    # ancestor in the (immediate) dominator tree between its semidominator and
    # its parent (see Lemma 2.6 in [LG05]).
    idoms_pre = copy(D.to_parent_pre)
    for v in 2:n_nodes
        idom = idoms_pre[v]
        vsemi = state[v].semi
        while idom > vsemi
            idom = idoms_pre[idom]
        end
        idoms_pre[v] = idom
    end

    # Express idoms in BB indexing
    resize!(domtree.idoms_bb, n_blocks)
    for i::BBNumber in 1:n_blocks
        if i == 1 || D.to_pre[i] == 0
            domtree.idoms_bb[i] = 0
        else
            domtree.idoms_bb[i] = D.from_pre[idoms_pre[D.to_pre[i]]]
        end
    end
end

"""
Matches the snca_compress algorithm in Figure 2.8 of [LG05], with the
modification suggested in the paper to use `last_linked` to determine whether
an ancestor has been processed rather than storing `0` in the ancestor array.
"""
function snca_compress!(state::Vector{SNCAData}, ancestors::Vector{PreNumber},
                        v::PreNumber, last_linked::PreNumber)
    u = ancestors[v]
    @assert u < v
    if u >= last_linked
        snca_compress!(state, ancestors, u, last_linked)
        if state[u].label < state[v].label
            state[v] = SNCAData(state[v].semi, state[u].label)
        end
        ancestors[v] = ancestors[u]
    end
    nothing
end

function snca_compress_worklist!(
        state::Vector{SNCAData}, ancestors::Vector{PreNumber},
        v::PreNumber, last_linked::PreNumber)
    # TODO: There is a smarter way to do this
    u = ancestors[v]
    worklist = Tuple{PreNumber, PreNumber}[(u,v)]
    @assert u < v
    while !isempty(worklist)
        u, v = last(worklist)
        if u >= last_linked
            if ancestors[u] >= last_linked
                push!(worklist, (ancestors[u], u))
                continue
            end
            if state[u].label < state[v].label
                state[v] = SNCAData(state[v].semi, state[u].label)
            end
            ancestors[v] = ancestors[u]
        end
        pop!(worklist)
    end
end

"Given updated blocks, update the given dominator tree with an inserted edge."
function domtree_insert_edge!(domtree::DomTree, blocks::Vector{BasicBlock},
                              from::BBNumber, to::BBNumber)
    # `from` is unreachable, so `from` and `to` aren't in domtree
    if bb_unreachable(domtree, from)
        return domtree
    end

    # Implements Section 3.1 of [GI16]
    dt        = domtree.dfs_tree
    from_pre  = dt.to_pre[from]
    to_pre    = dt.to_pre[to]
    from_post = dt.to_post[from]
    to_post   = dt.to_post[to]
    if to_pre == 0 || (from_pre < to_pre && from_post < to_post)
        # The DFS tree is invalidated by the edge insertion, so run from
        # scratch
        update_domtree!(blocks, domtree, true, 0)
    else
        # DFS tree is still valid, so update only affected nodes
        update_domtree!(blocks, domtree, false, to_pre)
    end

    return domtree
end

"Given updated blocks, update the given dominator tree with a deleted edge."
function domtree_delete_edge!(domtree::DomTree, blocks::Vector{BasicBlock},
                              from::BBNumber, to::BBNumber)
    # `from` is unreachable, so `from` and `to` aren't in domtree
    if bb_unreachable(domtree, from)
        return domtree
    end

    # Implements Section 3.1 of [GI16]
    if is_parent(domtree.dfs_tree, from, to)
        # The `from` block is the parent of the `to` block in the DFS tree, so
        # deleting the edge invalidates the DFS tree, so start from scratch
        update_domtree!(blocks, domtree, true, 0)
    elseif on_semidominator_path(domtree, from, to)
        # Recompute semidominators for blocks with preorder number up to that
        # of `to` block. Semidominators for blocks with preorder number greater
        # than that of `to` aren't affected because no semidominator path to
        # the block can pass through the `to` block (the preorder number of
        # `to` would be lower than those of these blocks, and `to` is not their
        # parent in the DFS tree).
        to_pre = domtree.dfs_tree.to_pre[to]
        update_domtree!(blocks, domtree, false, to_pre)
    end
    # Otherwise, dominator tree is not affected

    return domtree
end

"Check if x is the parent of y in the given DFS tree."
function is_parent(dfs_tree::DFSTree, x::BBNumber, y::BBNumber)
    x_pre = dfs_tree.to_pre[x]
    y_pre = dfs_tree.to_pre[y]
    return x_pre == dfs_tree.to_parent_pre[y_pre]
end

"""
Check if x is on some semidominator path from the semidominator of y to y,
assuming there is an edge from x to y.
"""
function on_semidominator_path(domtree::DomTree, x::BBNumber, y::BBNumber)
    x_pre = domtree.dfs_tree.to_pre[x]
    y_pre = domtree.dfs_tree.to_pre[y]

    semi_y = domtree.snca_state[y_pre].semi
    current_block = x_pre

    # Follow the semidominators of `x` up the DFS tree to see if we ever reach
    # the semidominator of `y`. If so, `x` is on a semidominator path between
    # `y` and its semidominator. We can stop if the preorder number of the
    # semidominators becomes less than that of the semidominator of `y`,
    # because it can only decrease further.
    while current_block >= semi_y
        if semi_y == current_block
            return true
        end
        current_block = domtree.snca_state[current_block].semi
    end
    return false
end

"""
Rename basic block numbers in a dominator tree, removing the block if it is
renamed to -1.
"""
function rename_nodes!(domtree::DomTree, rename_bb::Vector{BBNumber})
    # Rename DFS tree
    rename_nodes!(domtree.dfs_tree, rename_bb)

    # `snca_state` is indexed by preorder number, so should be unchanged

    # Rename `idoms_bb` and `nodes`
    old_idoms_bb = copy(domtree.idoms_bb)
    old_nodes = copy(domtree.nodes)
    for (old_bb, new_bb) in enumerate(rename_bb)
        if new_bb != -1
            domtree.idoms_bb[new_bb] = (new_bb == 1) ?
                0 : rename_bb[old_idoms_bb[old_bb]]
            domtree.nodes[new_bb] = old_nodes[old_bb]
            map!(i -> rename_bb[i],
                 domtree.nodes[new_bb].children,
                 domtree.nodes[new_bb].children)
        end
    end

    # length of `to_pre` after renaming DFS tree is new number of basic blocks
    resize!(domtree.idoms_bb, length(domtree.dfs_tree.to_pre))
    resize!(domtree.nodes, length(domtree.dfs_tree.to_pre))
    return domtree
end

"""
Rename basic block numbers in a DFS tree, removing the block if it is renamed
to -1.
"""
function rename_nodes!(D::DFSTree, rename_bb::Vector{BBNumber})
    n_blocks = length(D.to_pre)
    n_reachable_blocks = length(D.from_pre)

    old_to_pre = copy(D.to_pre)
    old_from_pre = copy(D.from_pre)
    old_to_post = copy(D.to_post)
    old_from_post = copy(D.from_post)
    max_new_bb = 0
    for (old_bb, new_bb) in enumerate(rename_bb)
        if new_bb != -1
            D.to_pre[new_bb] = old_to_pre[old_bb]
            D.from_pre[old_to_pre[old_bb]] = new_bb
            D.to_post[new_bb] = old_to_post[old_bb]
            D.from_post[old_to_post[old_bb]] = new_bb

            # Keep track of highest BB number to resize arrays with
            if new_bb > max_new_bb
                max_new_bb = new_bb
            end
        end
    end
    resize!(D.to_pre, max_new_bb)
    resize!(D.to_post, max_new_bb)
    # `to_parent_pre` should be unchanged
    return D
end

"""
Checks if bb1 dominates bb2.
bb1 and bb2 are indexes into the CFG blocks.
bb1 dominates bb2 if the only way to enter bb2 is via bb1.
(Other blocks may be in between, e.g bb1->bbX->bb2).
"""
function dominates(domtree::DomTree, bb1::BBNumber, bb2::BBNumber)
    bb1 == bb2 && return true
    target_level = domtree.nodes[bb1].level
    source_level = domtree.nodes[bb2].level
    source_level < target_level && return false
    for _ in (source_level - 1):-1:target_level
        bb2 = domtree.idoms_bb[bb2]
    end
    return bb1 == bb2
end

bb_unreachable(domtree::DomTree, bb::BBNumber) = bb != 1 && domtree.dfs_tree.to_pre[bb] == 0

"Iterable data structure that walks though all dominated blocks"
struct DominatedBlocks
    domtree::DomTree
    worklist::Vector{BBNumber}
end

"Returns an iterator that walks through all blocks dominated by the basic block at index `root`"
function dominated(domtree::DomTree, root::BBNumber)
    doms = DominatedBlocks(domtree, Vector{BBNumber}())
    push!(doms.worklist, root)
    doms
end

function iterate(doms::DominatedBlocks, state::Nothing=nothing)
    isempty(doms.worklist) && return nothing
    bb = pop!(doms.worklist)
    for dominated in doms.domtree.nodes[bb].children
        push!(doms.worklist, dominated)
    end
    return (bb, nothing)
end

function naive_idoms(blocks::Vector{BasicBlock})
    nblocks = length(blocks)
    # The extra +1 helps us detect unreachable blocks below
    dom_all = BitSet(1:nblocks+1)
    dominators = BitSet[n == 1 ? BitSet(1) : copy(dom_all) for n = 1:nblocks]
    changed = true
    while changed
        changed = false
        for n = 2:nblocks
            if isempty(blocks[n].preds)
                continue
            end
            firstp, rest = Iterators.peel(Iterators.filter(p->p != 0, blocks[n].preds))
            new_doms = copy(dominators[firstp])
            for p in rest
                intersect!(new_doms, dominators[p])
            end
            push!(new_doms, n)
            changed = changed || (new_doms != dominators[n])
            dominators[n] = new_doms
        end
    end
    # Compute idoms
    idoms = fill(0, nblocks)
    for i = 2:nblocks
        if dominators[i] == dom_all
            idoms[i] = 0
            continue
        end
        doms = collect(dominators[i])
        for dom in doms
            i == dom && continue
            hasany = false
            for p in doms
                if p !== i && p !== dom && dom in dominators[p]
                    hasany = true; break
                end
            end
            hasany && continue
            idoms[i] = dom
        end
    end
    idoms
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
export Agent, DynaAgent

include("agent.jl")
include("dyna_agent.jl")


text:
using Dates

add_gigasecond(date::DateTime) = date + Dates.Second(10^9)


text:
using Gadfly

set_default_plot_size(6inch, 8inch)

c = ["red", "blue"]
dodge1 = plot(x=["a","a"], y=[3,2], color=c, Geom.bar(position=:dodge))
dodge2 = plot(x=["a","a","b","b"], y=[3,2,1,2], color=[c;c], Geom.bar(position=:dodge))
stack1 = plot(x=["a","a"], y=[3,2], color=c, Geom.bar(position=:stack))
stack2 = plot(x=["a","a","b","b"], y=[3,2,1,2], color=[c;c], Geom.bar(position=:stack))
identity1 = plot(x=["a","a"], y=[3,2], color=c, alpha=[0.5], Geom.bar(position=:identity))
identity2 = plot(x=["a","a","b","b"], y=[3,2,1,2], color=[c;c], alpha=[0.5],
    Geom.bar(position=:identity))

gridstack([dodge1 dodge2; stack1 stack2; identity1 identity2])


text:
using Documenter, JuliaPackageTemplate

#
include("src/makeplots.jl")

# This function builds the documentation
makedocs(
    modules   = [JuliaPackageTemplate],
    doctest   = false,
    clean     = true,
    linkcheck = false,
    format    = Documenter.HTML(),
    sitename  = "JuliaPackageTemplate.jl",
    authors   = "Duncan Eddy",
    pages     = Any[
        "Home" => "index.md",
        "Modules" => Any[
            "modules/submodule.md", # Use default module name in sidebar
            "Rubber Ducks" => "modules/rubber_ducks.md", # Rename a module in sidebar
        ],
        "Examples" => Any[
            "Plotting Example" => "examples/plotting_example.md"
        ],
        "Library Index" => "library_index.md",
    ]
)

# Generate plots
# Note: Must be called after makedocs so the build folder are created
makeplots()

deploydocs(
    repo = "github.com/sisl/JuliaPackageTemplate.jl",
    devbranch = "master",
    devurl = "latest",
    deps = makeplots,
)