text:
export modifiedMatsuno

"""
    modifiedMatsuno(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the modified Matsuno model

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points as `mesh = Mesh(param)`

## Optional keyword arguments
- `ν`: shallow/deep water multiplication factor. By default, `ν=1` if `μ≦1` and `ν=1/√μ` otherwise. Set the infinite-layer case if `ν=0` (or `μ=Inf`).
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `verbose`: prints information if `true` (default is `true`).

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!` and in particular
1. a function `modifiedMatsuno.f!` to be called in the time-integration solver;
2. a function `modifiedMatsuno.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `modifiedMatsuno.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v)`, where
    - `η` is the surface deformation;
    - `v` is the derivative of the trace of the velocity potential.

"""
mutable struct modifiedMatsuno <: AbstractModel

    label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	param	:: NamedTuple
	kwargs	:: NamedTuple

    function modifiedMatsuno(param::NamedTuple;
							ν=nothing,ktol=0,dealias=0, verbose=true)

		label = "modified Matsuno"

		μ 	= param.μ
		ϵ 	= param.ϵ
		if ν == nothing
			if μ > 1
				ν = 1/sqrt(μ)
			else
				ν = 1
			end
		end

		mesh = Mesh(param)

		param = ( ϵ = ϵ, μ = μ, xmin = mesh.xmin, xmax = mesh.xmax, N = mesh.N )
		kwargs = (ν=ν,dealias=dealias,ktol=ktol,verbose=verbose)

		x = mesh.x
		k = copy(mesh.k)
		K = mesh.kmax * (1-dealias/(2+dealias))
		Π⅔ 	= abs.(mesh.k) .<= K # Dealiasing low-pass filter
		if dealias == 0
			if verbose @info "no dealiasing" end
			Π⅔ 	= ones(size(mesh.k))
		elseif verbose
			@info string("dealiasing : spectral scheme for power ", dealias + 1," nonlinearity ")
		end
		if μ == Inf || ν==0
			∂ₓF₀ 	= 1im * sign.(mesh.k)
			G₀ 	= abs.(mesh.k)
			μ = 1
			ν = 1
		else
        	∂ₓF₀ 	= 1im* sign.(mesh.k) .* tanh.(sqrt(μ)*abs.(mesh.k))
			G₀ 	= sqrt(μ)*abs.(mesh.k).*tanh.(sqrt(μ)*abs.(mesh.k))
		end
		∂ₓ	=  1im * sqrt(μ)* mesh.k            # Differentiation
		z = zeros(Complex{Float64}, mesh.N)
		η = copy(z) ; v = copy(z) ; F = copy(z) ;
		fftη = copy(z) ; fftv = copy(z) ; Q = copy(z) ; R = copy(z) ;

		# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
		function mapto(data::InitialData)
			U = [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
			U[abs.(U).< ktol ].=0
			return U
		end

		# Return `(η,v)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential.
		# Inverse Fourier transform and takes the real part.
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2]))
		end

		# Evolution equations are ∂t U = f!(U)
		function f!(U)

			fftη .= U[:,1]
		    fftv .= U[:,2]
			η  .= ifft(U[:,1])
			v  .= ifft(U[:,2])
			Q .= -∂ₓF₀.*fftv
			Q += -ϵ *∂ₓ.*fft(η.*ifft(fftv)) .+ ϵ*∂ₓF₀.*fft(η.*ifft(G₀.*fftv))
			F .= exp.(-ϵ*ifft(G₀.*fftη))
			R .= -fft(F.*ifft(∂ₓ.*fftη))*ν
			R += ϵ/2*∂ₓ.*fft(-v.^2)

			U[:,1] .= Π⅔.*Q/sqrt(μ)/ν
		   	U[:,2] .= Π⅔.*R/sqrt(μ)/ν
			U[abs.(U).< ktol ].=0

		end

        new(label, f!, mapto, mapfro, param, kwargs )
    end
end


text:
function correlation_index(
    spike_vactor1::Array,
    spike_vactor2::Array,
    delta_t::Float64,
    Time::Float64,
)
    total_spikes1 = length(spike_vactor1)
    total_spikes2 = length(spike_vactor2)

    interval_spikes = 0
    spike_index2 = 1

    for spike_index1 = 1:total_spikes1
        # for spike_index2=1:total_spikes2
        # TODO Shouldn't the spike_index2 be reset here?
        while (spike_index2 < total_spikes2)
            spike_interval =
                spike_vactor1[spike_index1] - spike_vactor2[spike_index2]
            if spike_itnerval > delta_t
                spike_index2 = spike_index2 + 1

            elseif fabs(spike_interval) <= delta_t
                interval_spikes += 1

                u = spike_index2 + 1

                while (
                    abs(spike_vactor1[spike_index1] - spike_vactor2[u]) <=
                    delta_t
                )
                    interval_spikes += 1
                    u += 1
                end
                break
            else
                break
            end
        end
    end

    signal_duration = Time[1] - Time[0]
    return index =
        (interval_spikes * signal_duration) /
        (total_spikes1 * total_spikes2 * 2 * delta_t)
end

ext: jl


text:
# Generates the file libleveldb_api.jl containing all calls to libleveldb and
# libleveldb_common.jl, containing all custom types
# this file should NOT be included inside LevelDB.jl

using Clang
using LevelDB

const LIBCLANG_INCLUDE =
    joinpath(pathof(Clang), "..", "..", "deps",
             "usr", "include", "clang-c") |>
    normpath
const LIBCLANG_HEADERS =
    [joinpath(LIBCLANG_INCLUDE, header)
     for header in readdir(LIBCLANG_INCLUDE)
     if endswith(header, ".h")]

const HEADERS =
    joinpath(pathof(LevelDB), "..", "..",
             "deps", "src", "leveldb-1.20",
             "include", "leveldb", "c.h") |>
                 normpath

wc = init(
    headers           = [HEADERS],
    output_file       = joinpath(pathof(LevelDB), "src", "libleveldb_api.jl"),
    common_file       = joinpath(pathof(LevelDB), "src", "libleveldb_common.jl"),
    clang_includes    = vcat(LIBCLANG_INCLUDE, CLANG_INCLUDE),
    clang_args        = ["-I", joinpath(LIBCLANG_INCLUDE, "..")],
    header_wrapped    = (root, current)->root == current,
    header_library    = x -> "libleveldb",
    clang_diagnostics = true,
)

run(wc)


text:
# Type that contains the data and functions associated to an ODE model in state-space form
type OdeModel
  # Initial values for the state variables
    States::NamedArray{Float64,1,Array{Float64,1},(Dict{ASCIIString,Int64},)}
  # Fixed parameters
    Parameters::NamedArray{Float64,1,Array{Float64,1},(Dict{ASCIIString,Int64},)}
  # Dynamic inputs
    Inputs::NamedArray{Function,1,Array{Function,1},(Dict{ASCIIString,Int64},)}
  # Functions implementing time derivatives and observed variables
    Model::Function
    Observer::Function
end

# Example: NamedArray([1,2.0], (["A" => 1,"B" => 2],),("States",))

# Getters and Setters

#=# Pretty-printing of lists of names
function pretty!(names, element)
    if(length(names) == 1)
        return "There is no $element called \"$(names[1])\""
    else
        output = "There are no $(element)s called \"$(names[1])\""
        if(length(names) > 2)
            for i in names[2:(end - 1)]
                output *= ", \"$i\""
            end
        end
        output *= " or \"$(names[end])\""
        return output
    end
end

# Check if all names_required are contained in actual_names and, if not, raise an informative error
# pointing to the element of OdeModel where the mismatch occured.
function check_names!(names_required, actual_names, element)
  not_contained = filter(x -> !in(x, actual_names), names_required)
  if(length(not_contained) > 0)
        error(pretty!(not_contained, element))
  end
end
=#


text:
function generate(
  ::Type{Transmission},
  tr::TransmissionRates,
  tnd::Nothing,
  id::Int64)
  external_or_internal = Weights([tr.external[id];
                                  sum(tr.internal[:, id])])
  if sum(external_or_internal) == 0.0
    @error "All transmission rates = 0.0, No transmission generated" external_or_internal
    return NoTransmission()
  elseif sample([true; false], external_or_internal)
    @debug "Exogenous tranmission generated" external_or_internal
    return ExogenousTransmission(id)
  else
    source = sample(1:individuals(tr), Weights(tr.internal[:, id]))
    @debug "Endogenous transmission generated (source id = $source)" external_or_internal
    return EndogenousTransmission(id, source)
  end
end

function generate(
  ::Type{Transmission},
  tr::TransmissionRates,
  tnd::TNDistribution,
  id::Int64)
  internalweights = tr.internal[:, id] .* tnd.internal[:, id]
  external_or_internal = Weights([tr.external[id] * tnd.external[id];
                                 sum(internalweights)])
  if sum(external_or_internal) == 0.0
    @error "All transmission rates = 0.0, No transmission generated" external_or_internal
    return NoTransmission()
  elseif sample([true; false], external_or_internal)
    @debug "Exogenous tranmission generated" external_or_internal
    return ExogenousTransmission(id)
  else
    source = sample(1:individuals(tr), Weights(internalweights))
    @debug "Endogenous transmission generated (source id = $source)" external_or_internal
    return EndogenousTransmission(id, source)
  end
end

function generate(
  ::Type{TransmissionNetwork},
  starting_states::DiseaseStates,
  tr::TransmissionRates,
  tnd::Union{Nothing, TNDistribution},
  ids::Vector{Int64})
  tn = TransmissionNetwork(starting_states)
  for i in ids
    tx = generate(Transmission, tr, tnd, i)
    update!(tn, tx)
  end
  return tn
end

function generate(
  ::Type{TransmissionNetwork},
  tr::TransmissionRates,
  mcmc::MCMC,
  ids::Vector{Int64})
  return generate(TransmissionNetwork, mcmc.starting_states, tr, mcmc.transmission_network_prior, ids)
end

function generate(
  ::Type{TransmissionNetwork},
  lasttn::TransmissionNetwork,
  tr::TransmissionRates,
  tnd::Union{Nothing, TNDistribution},
  ids)
  tn = copy(lasttn)
  for i in ids
    tn.external[i] = 0
    tn.internal[:, i] .= 0
    tx = generate(Transmission, tr, tnd, i)
    update!(tn, tx)
  end
  return tn
end

function generate(
  ::Type{TransmissionNetwork},
  lasttn::TransmissionNetwork,
  tr::TransmissionRates,
  mcmc::MCMC,
  ids)
  return generate(TransmissionNetwork, lasttn, tr, mcmc.transmission_network_prior, ids)
end

function generate(
  ::Type{Events},
  obs::EventObservations{T, M},
  extents::EventExtents{T}) where {
  T <: DiseaseStateSequence,
  M <: ILM}
  events = Events{T}(individuals(obs))
  exposed_state = State_E ∈ T
  removed_state = State_R ∈ T
  for i = 1:individuals(obs)
    if obs.infection[i] == -Inf
      update!(events, Event{T}(-Inf, i, State_I))
      if exposed_state
        update!(events, Event{T}(-Inf, i, State_E))
      end
      if removed_state && obs.removal[i] == -Inf
        update!(events, Event{T}(-Inf, i, State_R))
      end
    elseif !isnan(obs.infection[i])
      i_lb = obs.infection[i] - extents.infection[2]
      i_ub = obs.infection[i] - extents.infection[1]
      i_time = rand(Uniform(i_lb, i_ub))
      update!(events, Event{T}(i_time, i, State_I))
      if exposed_state
        e_lb = i_time - extents.exposure[2]
        e_ub = i_time - extents.exposure[1]
        e_time = rand(Uniform(e_lb, e_ub))
        update!(events, Event{T}(e_time, i, State_E))
      end
      if removed_state && !isnan(obs.removal[i])
        r_lb = maximum([obs.infection[i]; obs.removal[i] - extents.removal[2]])
        r_ub = obs.removal[i] - extents.removal[1]
        r_time = rand(Uniform(r_lb, r_ub))
        update!(events, Event{T}(r_time, i, State_R))
      end
    end
  end
  return events
end

function generate(
  ::Type{Events},
  mcmc::MCMC{T}) where {
  T <: DiseaseStateSequence}
  return generate(Events, mcmc.event_observations, mcmc.event_extents)
end

function generate(
  ::Type{Event},
  last_event::Event{T},
  σ::Float64,
  extents::EventExtents{T},
  obs::EventObservations,
  events::Events{T}) where {
  T <: DiseaseStateSequence}
  lowerbound, upperbound = _bounds(last_event, extents, obs, events)
  time = rand(truncated(Normal(_time(last_event), σ),
                        lowerbound,
                        upperbound))
  return Event(time, last_event)
end

function generate(
  ::Type{Event},
  last_event::Event{T},
  σ::Float64,
  extents::EventExtents{T},
  obs::EventObservations,
  events::Events{T},
  network::TransmissionNetwork) where {
  T <: DiseaseStateSequence}
  lowerbound, upperbound = _bounds(last_event, extents, obs, events, network)
  time = rand(truncated(Normal(_time(last_event), σ),
                        lowerbound,
                        upperbound))
  return Event(time, last_event)
end

function generate(
  ::Type{RiskParameters{T}},
  rpriors::RiskPriors{T}) where {
  T <: DiseaseStateSequence}
  sparks = Float64[rand(x) for x in rpriors.sparks]
  susceptibility = Float64[rand(x) for x in rpriors.susceptibility]
  infectivity = Float64[rand(x) for x in rpriors.infectivity]
  transmissibility = Float64[rand(x) for x in rpriors.transmissibility]
  if State_E ∈ T
    latency = Float64[rand(x) for x in rpriors.latency]
  end
  if State_R ∈ T
    removal = Float64[rand(x) for x in rpriors.removal]
  end
  if T == SEIR
    return RiskParameters{T}(sparks,
                             susceptibility,
                             infectivity,
                             transmissibility,
                             latency,
                             removal)
  elseif T == SEI
    return RiskParameters{T}(sparks,
                             susceptibility,
                             infectivity,
                             transmissibility,
                             latency)
  elseif T == SIR
    return RiskParameters{T}(sparks,
                             susceptibility,
                             infectivity,
                             transmissibility,
                             removal)
  elseif T == SI
    return RiskParameters{T}(sparks,
                             susceptibility,
                             infectivity,
                             transmissibility)
  end
end

function generate(::Type{RiskParameters}, mcmc::MCMC{T}) where T <: DiseaseStateSequence
  return generate(RiskParameters{T}, mcmc.risk_priors)
end

function generate(
  ::Type{RiskParameters{T}},
  last_rparams::RiskParameters{T},
  Σ::Array{Float64, 2}) where {
  T <: DiseaseStateSequence}
  rparams_vector = rand(MvNormal(convert(Vector{Float64}, last_rparams), Σ))
  return _like(last_rparams, rparams_vector)
end

function generate(
  ::Type{RiskParameters{T}},
  mc::MarkovChain{T},
  Σ::Array{Float64, 2}) where {
  T <: DiseaseStateSequence}
  return generate(RiskParameters{T}, mc.risk_parameters[end], Σ)
end

function generate(
  ::Type{SM},
  priors::Vector{UnivariateDistribution}) where {
  SM <: NucleicAcidSubstitutionModel}
  return SM([rand(x) for x in priors])
end

function generate(
  ::Type{SM},
  last_sm::SM,
  Σ::Array{Float64, 2}) where {
  SM <: NucleicAcidSubstitutionModel}
  if size(Σ) == (0,0)
    return last_sm
  else
    return SM(rand(MvNormal([getproperty(last_sm, θ) for θ in [fieldnames(SM)...]], Σ)), safe=false)
  end
end


text:
"""
    Point(x)

represents a single point at `x`.
"""
struct Point{T} <: Domain{T}
    x::T
end

similardomain(d::Point, ::Type{T}) where {T} = Point{T}(d.x)

convert(::Type{Number}, d::Point{<:Number}) = d.x
convert(::Type{N}, d::Point{<:Number}) where N<:Number = convert(N, convert(Number, d.x))
Number(d::Point) = convert(Number, d)

convert(::Type{Domain}, c::Number) = Point(c)
convert(::Type{Domain{T}}, c::Number) where T = Point{T}(c)

==(a::Point,b::Point) = a.x == b.x
indomain(x, d::Point) = x == d.x
isempty(::Point) = false

approx_indomain(x, d::Point, tolerance) = norm(x-d.x) <= tolerance

dimension(d::Point{Vector{T}}) where {T} = length(d.x)

canonicaldomain(d::Point{T}) where {T<:StaticTypes} = Point(zero(T))
canonicaldomain(d::Point{T}) where {T<:AbstractVector} =
    Point(zeros(eltype(T),dimension(d)))

mapfrom_canonical(d::Point) = Translation(d.x)

isopenset(d::Point) = false
isclosedset(d::Point) = true

boundary(d::Point) = d
boundingbox(d::Point) = d.x..d.x

interior(d::Point{T}) where {T} = EmptySpace{T}()
closure(d::Point) = d

point_in_domain(d::Point) = d.x

function map_domain(map, p::Point)
    x = applymap(map, p.x)
    Point(x)
end

mapped_domain(invmap, p::Point) = map_domain(inverse(invmap), p)

for op in (:+,:-)
    @eval $op(a::Point, b::Point) = Point($op(a.x,b.x))
end

# Interval minus a point:
setdiffdomain(d::Interval, x::Number) = setdiffdomain(d, Point(x))
setdiffdomain(d::Interval, p::Point) = setdiffdomain(promote_domains((d,p))...)
function setdiffdomain(d::Interval{L,R,T}, p::Point{T}) where {L,R,T}
    a = leftendpoint(d)
    b = rightendpoint(d)
    x = p.x

    a == x && return Interval{:open,R,T}(a,b)
    a < x < b && return UnionDomain(Interval{L,:open,T}(a,p.x), Interval{:open,R,T}(p.x,b))
    b == x && return Interval{L,:open,T}(a,b)
    return d
end

issubset1(d1::Point, d2) = d1.x ∈ d2

setdiffdomain1(p::Point, d2) = issubset(p, d2) ? EmptySpace{eltype(p)}() : p


text:
using PkgDev
using Test

@testset "PkgDev" begin
end


text:
## SCENARIO PARTITION
function GAPM_Vec(Ind_old, duals; α_val = 0.)
    NP = maximum(Ind_old)

    Scenarios = [Ind_old collect(1:N)]
        Ind,NsubParts = ones(N), 1

        dist = collect(abs(duals[i,t]-duals[j,t]) for i in 1:N, j in 1:N, t in T)
    for p in 1:NP
        DualPart = Scenarios[Scenarios[:,1].==p,:]

        sizePart = size(DualPart, 1)
        Rem = Integer.(DualPart[:,2]) # Remaining scenarios to assign
        while length(Rem)≥1

            indsP = collect(i for i in 2:length(Rem) if maximum(dist[Rem[1],Rem[i],:])≤α_val)
            indsP = [Rem[1]; Rem[indsP]]

            Ind[indsP] .= NsubParts

            RemInds = collect(i for i in 1:length(Rem) if !(Rem[i] in indsP))
            Rem = Rem[RemInds]
            NsubParts += 1
        end
        end
    return Integer.(Ind)
end

max_stars_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_issues_repo_name: gonzalez-alvaro/MultidimensionalGAPM

max_forks_repo_name: gonzalez-alvaro/MultidimensionalGAPM


text:
using MPIClusterManagers # to import MPIManager
using Distributed #provides addprocs()

manager=MPIManager(np=2)
addprocs(manager)

#Setup worker environments
@everywhere using Pkg
@everywhere Pkg.activate((@__DIR__)*"/..")
@everywhere using Plasmo
@everywhere using PipsNLP

graph = OptiGraph()

n1 = @optinode(graph)
n2 = @optinode(graph)

@variable(n1, 0 <= x <= 2)
@variable(n1, 0 <= y <= 3)
@variable(n1, z >= 0)
@constraint(n1, x+y+z >= 4)
@objective(n1, Min, y)

@variable(n2,x)
@NLconstraint(n2,ref,exp(x) >= 2)
@variable(n2,z >= 0)
@constraint(n2,z + x >= 4)
@objective(n2,Min,x)

@linkconstraint(graph,n1[:x] == n2[:x])

#distribute the graph to workers
julia_workers = collect(values(manager.mpi2j))
remote_references = PipsNLP.distribute_optigraph(graph,julia_workers,remote_name = :pipsgraph)

# The remote optigraphs can be queried if they are fetched from the other workers
# NOTE: this will copy the optigraph to the main process
# r1 = fetch(remote_references[1])
# r2 = fetch(remote_references[2])

#Solve with pips-nlp
@mpi_do manager begin
    using MPI
    PipsNLP.pipsnlp_solve(pipsgraph)
end

#fill the local solution with the result on worker 2 (worker 2 is the root MPI rank)
PipsNLP.fill_solution!(graph,:pipsgraph,2)


text:
# DSANet.jl
#
# Layers for a DSANet network
# Based on paper Siteng Huang et. al. "DSANet: Dual Self-Attention Network for Multivariate
# Time Series Forecasting" https://kyonhuang.top/files/Huang-DSANet.pdf

# Helper Functions
"""
    Scaled_Dot_Product_Attention(q, k, v, temperature)

Scaled dot product attention function with query `q`, keys `k` and values `v`. Normalisation
is given by `temperature`. Outputs
``\\mathrm{softmax}\\left( \\frac{q \\cdot k^T}{\\mathrm{temperature}} \\right)\\cdot v``.
"""
function Scaled_Dot_Product_Attention(q, k, v, temperature, attn_dropout=0.1)
    #attn1 = NNlib.batched_mul(q, Flux.batched_transpose(k)) / temperature
	attn1 = NNlib.batched_mul(q, permutedims(k,[2,1,3])) / temperature
	attn2 = Flux.softmax(attn1, dims=2)
    attn3 = Dropout(attn_dropout)(attn2)
    return NNlib.batched_mul(attn3,v)
end

function regularized_normalise(x::AbstractArray; dims=1)
  μ′ = Flux.mean(x, dims = dims)
  σ′ = Flux.std(x, dims = dims, mean = μ′, corrected=false)
  ϵ = sqrt(eps(Float32))
  return (x .- μ′) ./ (σ′ .+ ϵ)
end

"""
    Reg_LayerNorm(h::Integer)
A [normalisation layer](https://arxiv.org/pdf/1607.06450.pdf) designed to be
used with recurrent hidden states of size `h`. Normalises the mean and standard
deviation of each input before applying a per-neuron gain/bias. To avoid numeric
overflow, the division by the standard deviation has been regularised by
adding `ϵ = 1E-5`.
"""
struct Reg_LayerNorm{T}
  diag::Flux.Diagonal{T}
end

Reg_LayerNorm(h::Integer) =
  Reg_LayerNorm(Flux.Diagonal(h))

Flux.@functor Reg_LayerNorm

(a::Reg_LayerNorm)(x) = a.diag(regularized_normalise(x))

function Base.show(io::IO, l::Reg_LayerNorm)
  print(io, "Reg_LayerNorm(", length(l.diag.α), ")")
end

# Encoder structure
mutable struct SelfAttn_EncoderCell{A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P}
    # General
    dropout::A
    # Multi head attention
    layer_norm::B  # currently not used
    mha_w_qs::C
    mha_w_ks::D
    mha_w_vs::E
    mha_fc::F
    # Positionwise Feedforward
    pwff_w_1::G
    pwff_w_2::H
	pwff_layer_norm::P
    # Metadata
    inp::I
    n_head::J
    d_model::K
    d_hid::L
    d_k::M
    d_v::N
    σ::O
end

Flux.@functor SelfAttn_EncoderCell
Flux.trainable(m::SelfAttn_EncoderCell) = (m.layer_norm, m.mha_w_qs, m.mha_w_ks, m.mha_w_vs,
            m.mha_fc, m.pwff_w_1, m.pwff_w_2, m.pwff_layer_norm)

"""
    SelfAttn_Encoder(inp, d_model, n_head, d_hid, drop_prob = 0.1f0, σ = Flux.relu)

Encoder part for the self attention networks that comprise the `DSANet`. For parameters see
[DSANet](@ref).
"""
function SelfAttn_Encoder(inp, d_model, n_head, d_hid, drop_prob = 0.1f0, σ = Flux.relu)
    (d_model % n_head == 0) || throw(ArgumentError("d_model = $d_model must be divisible by n_head = $n_head"))

    d_k = d_model ÷ n_head
    d_v = d_model ÷ n_head

    # Multi head attention
    mha_w_qs = Dense(d_model, n_head * d_k)
    mha_w_ks = Dense(d_model, n_head * d_k)
    mha_w_vs = Dense(d_model, n_head * d_v)
    mha_fc = Dense(n_head*d_v, d_model)
    layer_norm = Reg_LayerNorm(d_model)
    # Positionwise Feedforward
    pwff_w_1 = Conv((1,), d_model => d_hid)
    pwff_w_2 = Conv((1,), d_hid => d_model)
	pwff_layer_norm = Reg_LayerNorm(d_model)

    return SelfAttn_EncoderCell(Dropout(drop_prob), layer_norm, mha_w_qs, mha_w_ks, mha_w_vs, mha_fc,
                        pwff_w_1, pwff_w_2, pwff_layer_norm, inp, n_head, d_model, d_hid, d_k, d_v, σ)
end

function Base.show(io::IO, m::SelfAttn_EncoderCell)
	print(io, "SelfAttn_Encoder(", m.inp)
	print(io, ", ", m.d_model)
	print(io, ", ", m.n_head)
	print(io, ", ", m.d_hid)
	m.dropout.p == 0.1 || print(io, ", ", m.dropout.p)
	m.σ == Flux.relu || print(io, ", ", m.σ)
	print(io, ")")
end

function (m::SelfAttn_EncoderCell)(src_seq)
    # Multi head attention
    q = k = v = src_seq
    mha_residual = q
    q1 = reshape(q, (m.d_model, :))
    q2 = reshape(m.mha_w_qs(q1), (m.n_head, m.d_k, :) )
	k1 = reshape(k, (m.d_model, :))
	k2 = reshape(m.mha_w_ks(k1), (m.n_head, m.d_k,:) )
	v1 = reshape(v, (m.d_model, :))
	v2 = reshape(m.mha_w_ks(v1), (m.n_head, m.d_v,:) )

    q3 = permutedims(q2,[1,3,2])
    k3 = permutedims(k2,[1,3,2])
    v3 = permutedims(v2,[1,3,2])

    output1 = Scaled_Dot_Product_Attention(q3, k3, v3, sqrt(Float32(m.d_k)))
	output2 = permutedims(output1, [1,3,2])
    output3 = reshape(output2, (m.d_model, :))
    output4 = m.mha_fc(output3)
    output5 = reshape(output4,(m.d_model, m.inp,:))
    output6 = m.dropout(output5)
	# this line crashes the training by creating NaNs in the Conv layer parameters
	# output7 = m.layer_norm(output6 .+ mha_residual)

    # Positionwise feedforward
    pwff_residual = output6 .+ mha_residual  # output7
    #output1 = Flux.batched_transpose(pwff_residual)
	output1 = permutedims(pwff_residual, [2,1,3])
    output2 = σ.(m.pwff_w_1(output1))
    #output3 = Flux.batched_transpose(m.pwff_w_2(output2))
	output3 = permutedims(m.pwff_w_2(output2), [2,1,3])
    output4 = m.dropout(output3)
    output5 = m.pwff_layer_norm(output4 .+ pwff_residual)
    return output5
end

# Global Self Attention
mutable struct Global_SelfAttn_Cell{A, B, C, D, E, F, G, H, I, J}
    conv::A
    dropout::B
    in_linear::C
    out_linear::D
    encoder_chain::E
    # Metadata
    inp::F
    d_model::G
    n_kernels::H
    window::I
    σ::J
end

Flux.@functor Global_SelfAttn_Cell
Flux.trainable(m::Global_SelfAttn_Cell) = (m.conv, m.in_linear, m.out_linear,
            m.encoder_chain)

"""
    Global_SelfAttn(inp, window, n_kernels, w_kernel, d_model, d_hid, n_layers, n_head)

Global self attention module for `DSANet`. For parameters see [DSANet](@ref).
"""
function Global_SelfAttn(inp, window, n_kernels, w_kernel, d_model, d_hid,
                n_layers, n_head, drop_prob = 0.1f0, σ = Flux.relu)

    conv = Conv((w_kernel, window), 1 => n_kernels)
    dropout = Dropout(drop_prob)
    in_linear = Dense(n_kernels, d_model)
    out_linear = Dense(d_model, n_kernels)
    encoder_chain = Chain([SelfAttn_Encoder(inp, d_model, n_head, d_hid, drop_prob, σ) for i=1:n_layers]...)
    return Global_SelfAttn_Cell(conv, dropout, in_linear, out_linear, encoder_chain,
                inp, d_model, n_kernels, window, σ)
end

function Base.show(io::IO, m::Global_SelfAttn_Cell)
	print(io, "Global_SelfAttn(", m.inp)
	print(io, ", ", m.window)
	print(io, ", ", m.n_kernels)
	print(io, ", ", size(m.conv.weight,1))
	print(io, ", ", m.d_model)
	print(io, ", ", m.encoder_chain[1].d_hid)
	print(io, ", ", length(m.encoder_chain))
	print(io, ", ", m.encoder_chain[1].n_head)
	m.dropout.p == 0.1 || print(io, ", ", m.dropout.p)
	m.σ == Flux.relu || print(io, ", ", m.σ)
	print(io, ")")
end

function (m::Global_SelfAttn_Cell)(x)
	x1 = m.conv(x)
    x2 = m.dropout(m.σ.(x1))
    x3 = dropdims(x2;dims=2)
    #x4 = Flux.batched_transpose(x3)
	x4 = permutedims(x3, [2,1,3])
    x5 = reshape(x4, (m.n_kernels, :))
    x6 = m.in_linear(x5)
    src_seq = reshape(x6, (m.d_model,m.inp,:))
    x8 = m.encoder_chain(src_seq)
    x9 = reshape(x8, (m.d_model, :))
    x10 = m.out_linear(x9)
    x11 = reshape(x10, (m.n_kernels, m.inp,:))
    return x11
end

# Local Self Attention
mutable struct Local_SelfAttn_Cell{A, B, C, D, E, F, G, H, I, J, K}
    conv::A
    dropout::B
    in_linear::C
    out_linear::D
    encoder_chain::E
    # Metadata
    inp::F
    d_model::G
    n_kernels::H
    window::I
    local_length::J
    σ::K
end

Flux.@functor Local_SelfAttn_Cell
Flux.trainable(m::Local_SelfAttn_Cell) = (m.conv, m.in_linear, m.out_linear,
                m.encoder_chain)

"""
    Local_SelfAttn(inp, window, local_length, n_kernels, w_kernel, d_model, d_hid, n_layers, n_head)

Local self attention module for `DSANet`. For parameters see [DSANet](@ref).
"""
function Local_SelfAttn(inp, window, local_length, n_kernels, w_kernel, d_model, d_hid,
                n_layers, n_head, drop_prob = 0.1f0, σ = Flux.relu)

    conv = Conv((w_kernel, local_length), 1 => n_kernels)
    dropout = Dropout(drop_prob)
    in_linear = Dense(n_kernels, d_model)
    out_linear = Dense(d_model, n_kernels)
    encoder_chain = Chain([SelfAttn_Encoder(inp, d_model, n_head, d_hid, drop_prob, σ) for i=1:n_layers]...)
    return Local_SelfAttn_Cell(conv, dropout, in_linear, out_linear, encoder_chain,
                inp, d_model, n_kernels, window, local_length, σ)
end

function Base.show(io::IO, m::Local_SelfAttn_Cell)
	print(io, "Local_SelfAttn(", m.inp)
	print(io, ", ", m.window)
	print(io, ", ", m.local_length)
	print(io, ", ", m.n_kernels)
	print(io, ", ", size(m.conv.weight,1))
	print(io, ", ", m.d_model)
	print(io, ", ", m.encoder_chain[1].d_hid)
	print(io, ", ", length(m.encoder_chain))
	print(io, ", ", m.encoder_chain[1].n_head)
	m.dropout.p == 0.1 || print(io, ", ", m.dropout.p)
	m.σ == Flux.relu || print(io, ", ", m.σ)
	print(io, ")")
end

function (m::Local_SelfAttn_Cell)(x)
	x1 = m.σ.(m.conv(x))
    x2 = MaxPool((1,m.window-m.local_length+1))(x1)
    x3 = m.dropout(x2)
    #x4 = Flux.batched_transpose(dropdims(x3, dims=(2)))
	x4 = permutedims(dropdims(x3, dims=(2)), [2,1,3])
    x5 = reshape(x4, (m.n_kernels, :))
    src_seq = reshape(m.in_linear(x5),(m.d_model, m.inp, :))
    x6 = m.encoder_chain(src_seq)
    x7 = reshape(x6, (m.d_model, :))
    x8 = m.out_linear(x7)
    x9 = reshape(x8, (m.n_kernels, m.inp,:))
    return x9
end

# Full DSANet
mutable struct DSANetCell{A, B, C, D, E, F, G, H, I, J}
    gsa::A
    lsa::B
    AR_linear::C
    output::D
	final_output::E
    dropout::F
    # Metadata
    inp::G
    n_kernels::H
    window::I
    output_dimensions::J
end

Flux.@functor DSANetCell
Flux.trainable(m::DSANetCell) = (m.gsa, m.lsa, m.AR_linear, m.output, m.final_output)

"""
    DSANet(inp, window, local_length, n_kernels, d_model, d_hid, n_layers, n_head, out=1, drop_prob = 0.1f0, σ = Flux.relu)

Creates a `DSANet` network based on the architecture described in
[Siteng Huang et. al.](https://kyonhuang.top/publication/dual-self-attention-network). The
code follows the [PyTorch implementation](https://github.com/bighuang624/DSANet).
`inp` specifies the number of input features. `window` gives the length of the window
for the pooled input data. `local_length` defines the length of the convolution window for
the local self attention mechanism. `n_kernel` defines the number of convolution kernels
for both the local and global self attention mechanism. `d_hid` defines the number of
"hidden" convolution kernels in the self attention encoder structure.
`n_layers` gives the number of self attention encoders used in the network, and `n_head`
defines the number of attention heads. `out` gives the number of output time series,
`drop_prob` is the dropout probability for the `Dropout` layers, and `σ` defines the
network's activation function.

Data is expected as array with dimensions `features x poolsize x 1 x data`, i.e.
for 1000 data points containing 31 features that have been windowed over 6
timesteps, `DSANet` expects an input size of `(31, 6, 1, 1000)`.
"""
function DSANet(inp, window, local_length, n_kernels, d_model, d_hid,
                n_layers, n_head, out=1, drop_prob = 0.1f0, σ = Flux.relu)
    out <= inp || throw(ArgumentError("Number of inputs ($inp) needs to be larger or equal to outputs ($out) "))
	w_kernel = 1
	gsa = Global_SelfAttn(inp, window, n_kernels, w_kernel, d_model, d_hid,
                    n_layers, n_head, drop_prob, σ)
    lsa = Local_SelfAttn(inp, window, local_length, n_kernels, w_kernel, d_model, d_hid,
                    n_layers, n_head, drop_prob, σ)
    AR_linear = Dense(window, 1)
    output = Dense(2*n_kernels, 1, σ)
	final_output = Dense(inp, out)
    dropout = Dropout(drop_prob)
    return DSANetCell(gsa, lsa, AR_linear, output, final_output, dropout, inp, n_kernels,
                    window, out)
end

function Base.show(io::IO, m::DSANetCell)
	print(io, "DSANet(", m.inp)
	print(io, ", ", m.window)
	print(io, ", ", m.lsa.local_length)
	print(io, ", ", m.lsa.n_kernels)
	#print(io, ", ", size(m.lsa.conv.weight,1))  # w_kernel
	print(io, ", ", m.lsa.d_model)
	print(io, ", ", m.lsa.encoder_chain[1].d_hid)
	print(io, ", ", length(m.lsa.encoder_chain))
	print(io, ", ", m.lsa.encoder_chain[1].n_head)
	m.output_dimensions == 1 || print(io, ", ", m.output_dimensions)
	m.dropout.p == 0.1f0 || print(io, ", ", m.dropout.p)
	m.lsa.σ == Flux.relu || print(io, ", ", m.lsa.σ)
	print(io, ")")
end

function (m::DSANetCell)(x)
	# Autoregressive Part
    x1 = dropdims(x, dims=3)
    #x2 = Flux.batched_transpose(x1)
	x2 = permutedims(x1, [2,1,3])
    x3 = reshape(x2, (m.window,:))
    x4 = m.AR_linear(x3)
    ar_output = reshape(x4,(m.inp,:))

    # Self Attention Modules
    gsa_output = m.gsa(x)
    lsa_output = m.lsa(x)

    # Full Net
    sf_output = cat(gsa_output, lsa_output, dims=1)
    sf_output1 = m.dropout(sf_output)
    sf_output2 = reshape(sf_output1, (2*m.n_kernels,:))
    sf_output3 = m.output(sf_output2)
    sf_output4 = reshape(sf_output3, (m.inp,:))

    out = sf_output4 .+ ar_output
    return m.final_output(out)
end

max_stars_repo_name: cbdavis/FluxArchitectures

max_issues_repo_name: cbdavis/FluxArchitectures

max_forks_repo_name: cbdavis/FluxArchitectures


text:
module TestStructuredBroadcast
using Test, LinearAlgebra

@testset "broadcast[!] over combinations of scalars, structured matrices, and dense vectors/matrices" begin
    N = 10
    s = rand()
    fV = rand(N)
    fA = rand(N, N)
    Z = copy(fA)
    D = Diagonal(rand(N))
    B = Bidiagonal(rand(N), rand(N - 1), :U)
    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
    U = UpperTriangular(rand(N,N))
    L = LowerTriangular(rand(N,N))
    structuredarrays = (D, B, T, U, L)
    fstructuredarrays = map(Array, structuredarrays)
    for (X, fX) in zip(structuredarrays, fstructuredarrays)
        @test (Q = broadcast(sin, X); typeof(Q) == typeof(X) && Q == broadcast(sin, fX))
        @test broadcast!(sin, Z, X) == broadcast(sin, fX)
        @test (Q = broadcast(cos, X); Q isa Matrix && Q == broadcast(cos, fX))
        @test broadcast!(cos, Z, X) == broadcast(cos, fX)
        @test (Q = broadcast(*, s, X); typeof(Q) == typeof(X) && Q == broadcast(*, s, fX))
        @test broadcast!(*, Z, s, X) == broadcast(*, s, fX)
        @test (Q = broadcast(+, fV, fA, X); Q isa Matrix && Q == broadcast(+, fV, fA, fX))
        @test broadcast!(+, Z, fV, fA, X) == broadcast(+, fV, fA, fX)
        @test (Q = broadcast(*, s, fV, fA, X); Q isa Matrix && Q == broadcast(*, s, fV, fA, fX))
        @test broadcast!(*, Z, s, fV, fA, X) == broadcast(*, s, fV, fA, fX)
        for (Y, fY) in zip(structuredarrays, fstructuredarrays)
            @test broadcast(+, X, Y) == broadcast(+, fX, fY)
            @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
            @test broadcast(*, X, Y) == broadcast(*, fX, fY)
            @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
        end
    end
    diagonals = (D, B, T)
    fdiagonals = map(Array, diagonals)
    for (X, fX) in zip(diagonals, fdiagonals)
        for (Y, fY) in zip(diagonals, fdiagonals)
            @test broadcast(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY)
            @test broadcast!(+, Z, X, Y) == broadcast(+, fX, fY)
            @test broadcast(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY)
            @test broadcast!(*, Z, X, Y) == broadcast(*, fX, fY)
        end
    end
end

@testset "broadcast! where the destination is a structured matrix" begin
    N = 5
    A = rand(N, N)
    sA = A + copy(A')
    D = Diagonal(rand(N))
    B = Bidiagonal(rand(N), rand(N - 1), :U)
    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
    @test broadcast!(sin, copy(D), D) == Diagonal(sin.(D))
    @test broadcast!(sin, copy(B), B) == Bidiagonal(sin.(B), :U)
    @test broadcast!(sin, copy(T), T) == Tridiagonal(sin.(T))
    @test broadcast!(*, copy(D), D, A) == Diagonal(broadcast(*, D, A))
    @test broadcast!(*, copy(B), B, A) == Bidiagonal(broadcast(*, B, A), :U)
    @test broadcast!(*, copy(T), T, A) == Tridiagonal(broadcast(*, T, A))
end

@testset "map[!] over combinations of structured matrices" begin
    N = 10
    fA = rand(N, N)
    Z = copy(fA)
    D = Diagonal(rand(N))
    B = Bidiagonal(rand(N), rand(N - 1), :U)
    T = Tridiagonal(rand(N - 1), rand(N), rand(N - 1))
    U = UpperTriangular(rand(N,N))
    L = LowerTriangular(rand(N,N))
    structuredarrays = (D, B, T, U, L)
    fstructuredarrays = map(Array, structuredarrays)
    for (X, fX) in zip(structuredarrays, fstructuredarrays)
        @test (Q = map(sin, X); typeof(Q) == typeof(X) && Q == map(sin, fX))
        @test map!(sin, Z, X) == map(sin, fX)
        @test (Q = map(cos, X); Q isa Matrix && Q == map(cos, fX))
        @test map!(cos, Z, X) == map(cos, fX)
        @test (Q = map(+, fA, X); Q isa Matrix && Q == map(+, fA, fX))
        @test map!(+, Z, fA, X) == map(+, fA, fX)
        for (Y, fY) in zip(structuredarrays, fstructuredarrays)
            @test map(+, X, Y) == map(+, fX, fY)
            @test map!(+, Z, X, Y) == map(+, fX, fY)
            @test map(*, X, Y) == map(*, fX, fY)
            @test map!(*, Z, X, Y) == map(*, fX, fY)
            @test map(+, X, fA, Y) == map(+, fX, fA, fY)
            @test map!(+, Z, X, fA, Y) == map(+, fX, fA, fY)
        end
    end
    diagonals = (D, B, T)
    fdiagonals = map(Array, diagonals)
    for (X, fX) in zip(diagonals, fdiagonals)
        for (Y, fY) in zip(diagonals, fdiagonals)
            @test map(+, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(+, fX, fY)
            @test map!(+, Z, X, Y) == broadcast(+, fX, fY)
            @test map(*, X, Y)::Union{Diagonal,Bidiagonal,Tridiagonal} == broadcast(*, fX, fY)
            @test map!(*, Z, X, Y) == broadcast(*, fX, fY)
        end
    end
end

end

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
@gen function involutive_mcmc_test_model()
    x ~ bernoulli(0.5)
    if x
        y ~ normal(0, 1)
    else
        z ~ normal(1, 1)
    end
end

@gen function involutive_mcmc_test_proposal(trace)
end

@testset "Involutive MCMC" begin

    @transform dsl_involution (model_in, aux_in) to (model_out, aux_out) begin
        x = @read(model_in[:x], :discrete)
        @write(model_out[:x], !x, :discrete)
        if x
            @copy(model_in[:y], model_out[:z])
        else
            @copy(model_in[:z], model_out[:y])
        end
    end

    function julia_involution(trace, forward_choices, forward_retval, q_args)
        argdiffs = map((_) -> NoChange(), get_args(trace))
        constraints = choicemap((:x, !trace[:x]))
        if trace[:x]
            constraints[:z] = trace[:y]
        else
            constraints[:y] = trace[:z]
        end
        (new_trace, weight, _, _) = update(trace, get_args(trace), argdiffs, constraints)
        # NOTE: Jacobian is 1.0
        return (new_trace, choicemap(), weight)
    end

    trace, _ = generate(involutive_mcmc_test_model, (), choicemap((:x, true), (:y, 0.5)))
    for involution in [dsl_involution, julia_involution]

        new_trace, accepted = involutive_mcmc(trace, involutive_mcmc_test_proposal, (), involution; check=true)
        # it should be accepted since the new and old densities are the same (0.5 is halfway between means)
        @test accepted
        @test !new_trace[:x]
        @test new_trace[:z] == 0.5
    end

end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using GLWindow, GLFW
export select, update, delete, addchar, newline, TextField1D, TextField2D, Updated

function delete(s::UTF8String, r::Range{Int})
	i = chr2ind(s, first(r))
	j = chr2ind(s, last(r))
	s[1:prevind(s, i)] * s[nextind(s, j):end]
end

function delete(s::UTF8String, i::Int)
	if i == 0
		return s
	elseif i == length(s)
		return s[1:end-1]
	end
	I = chr2ind(s, i)
	return s[1:prevind(s, I)] * s[nextind(s, I):end]
end

function delete(event::KeyPressed, t::TextField)
	i = first(t.selection)
	if length(t.selection) > 0
		t.text = delete(t.text, t.selection)
	else
		t.text = delete(t.text, i)
	end
	t.selection = max(i-1, 0) : i-2
	update(t)
	publishEvent(Updated(t))
end

function newline(event::KeyPressed, t::TextField)
	addchar(event, t)
	publishEvent(Updated(t))
end

addchar(s::UTF8String, char::Char, i::Int) = addchar(s, utf8(string(char)), i)

function addchar(s::UTF8String, char::UTF8String, i::Int)
	if i == 0
		return char * s
	elseif i == length(s)
		return s * char
	end
	I 		= chr2ind(s, i)
	startI 	= nextind(s, I)
	return s[1:I] *char* s[startI:end]
end

function addchar(event::UnicodeInput, t::TextField)

	char = utf8(string(event.char))
	i = first(t.selection)
	if length(t.selection) > 0
		t.text = delete(t.text, t.selection)
		i -= 1
	end
	t.text = addchar(t.text, char, i)
	i = min(i, length(t.text)-1)
	t.selection = i+1:i
	update(t)
	publishEvent(Updated(t))
end

function update(t::TextField)
	t.newLineIndexes 	= build_line_indexes(t.text)
	defaultStyle 		= Dict{ASCIIString, Any}(["textColor" => Float32[1,1,1,1], "backgroundColor" => Float32[0,0,0,0]])
	t.styles 			= [StyledTextSegment(1:length(t.text), defaultStyle)]
end

function select(event::MouseClicked, t::TextField, f::FontProperties)
	changed, newSelection = select(event.x, event.y, t, f)
	if changed
		t.selection = newSelection
	end
end

function select(event::MouseDragged, t::TextField, f::FontProperties)
	changed, startSelection = select(event.start.x, event.start.y, t, f)
	changed, endSelection 	= select(event.x, event.y, t, f)
	if changed
		t.selection = min(first(startSelection),first(endSelection)):max(first(startSelection),first(endSelection))
	end
end
function select(x::Real, y::Real, t::TextField, f::FontProperties)
	# is x, y inside the textfield?
	if y <= t.y + f.lineHeight && y >= t.y - f.lineHeight * length(t.newLineIndexes)
		#convert into Glyph coordinates
		lineIndex		= max(min(div(t.y + f.lineHeight - y, f.lineHeight) + 1, length(t.newLineIndexes)), 1)
		line 			= t.newLineIndexes[lineIndex]
		xPosCursor::Int	= max(min(div(x - t.x, f.advance) + first(line), last(line)), 0)
		return (true, xPosCursor:xPosCursor-1)
	end
	return (false, 1:0)
end

function select(direction::ASCIIString, t::TextField, f::FontProperties)
	if !isempty(t.text)
		cursor = first(t.selection)
		currentLineIndex, currentLine = findline(t.newLineIndexes, cursor)
		currentLine = first(currentLine)-1:last(currentLine)
		newLine 	= currentLineIndex
		newIndex    = cursor - first(currentLine)
		if direction == "left"
			if newIndex > 0
				newIndex -= 1
			else
				if newLine > 1
					newLine -= 1
					newIndex = length(t.newLineIndexes[newLine])
				end
			end
		elseif direction == "right"

			if newIndex < length(currentLine) - 1
				newIndex += 1
			else
				if newLine < length(t.newLineIndexes)
					newLine += 1
					newIndex = 0
				end
			end
		elseif direction == "up"
			if newLine > 1
				newLine -= 1
				newIndex = min(newIndex, length(t.newLineIndexes[newLine]))
			end
		elseif direction == "down"
			if newLine < length(t.newLineIndexes)
				newLine += 1
				newIndex = min(newIndex, length(t.newLineIndexes[newLine]))
			end
		end
		currentLine = t.newLineIndexes[newLine]
		newIndex 	= first(currentLine) - 1 + newIndex
		return newIndex:newIndex-1
	end
	t.selection
end

scroll(event, t) = t.y += event.yOffset * 30

function render(t::TextField, f::GLFont)
	render(f.gl)
	render(t, f)
end

changefocus(event, t, rect) = t.hasFocus = inside(rect, event.x, event.y)

function TextField2D(id::Symbol, text::String, x, y, area::Shape)
	font 	= getfont()
	t 		= TextField(id, text, x, y, area)

	registerEventAction(MouseClicked{Window}, leftclickdown, changefocus, (t, area))

	registerEventAction(KeyPressed{Window}, (x, tx) -> x.key == GLFW.KEY_BACKSPACE && tx.hasFocus  , (t, ), delete, (t,))

	registerEventAction(UnicodeInput{Window}, (_, tx) -> tx.hasFocus, (t,), addchar, (t,))

	registerEventAction(KeyPressed{Window}, (x,tx) -> tx.hasFocus && x.key == GLFW.KEY_ENTER, (t,), newline, (t,))

	registerEventAction(MouseClicked{Window}, leftclickup, select, (t, font.properties))

	registerEventAction(Scrolled{Window}, (_, tx) -> tx.hasFocus, (t,), scroll, (t,))

	#registerEventAction(KeyPressed, (x,tx) -> tx.hasFocus && (x.key == KEY_UP || x.key == KEY_DOWN || x.key == KEY_RIGHT || x.key == KEY_LEFT), (t,),
	#	select, (t, font.properties))

	registerEventAction(MouseDragged{Window}, x -> true, select, (t, font.properties))

	glDisplay(id, (t, font))

	t
end
function TextField1D(id::Symbol, text::String, x, y, area::Shape)
	font 	= getfont()
	text 	= replace(text, r"\n|\r", "")
	t 		= TextField(id, text, x, y, area)

	registerEventAction(MouseClicked{Window}, leftclickdown, changefocus, (t, area))

	registerEventAction(KeyPressed{Window}, (x, tx) -> x.key == KEY_BACKSPACE && tx.hasFocus, (t, ), delete, (t,))

	registerEventAction(UnicodeInput{Window}, (_, tx) -> tx.hasFocus, (t,), addchar, (t,))

	registerEventAction(MouseClicked{Window}, leftclickup, select, (t, font.properties))

	registerEventAction(MouseDragged{Window}, x -> true, select, (t, font.properties))

	glDisplay(id, (t, font))
	t
end

immutable Updated{T} <: Event
	source::T
end


text:
#### Surface model kernels

using Statistics

"""
    subdomain_surface_values(
        surf::SurfaceModel,
        turbconv::EDMF{FT},
        atmos::AtmosModel{FT},
        state::Vars,
        aux::Vars,
        zLL::FT,
    ) where {FT}

Returns the surface values of updraft area fraction, updraft
liquid water potential temperature (`θ_liq`), updraft total
water specific humidity (`q_tot`), environmental variances of
`θ_liq` and `q_tot`, environmental covariance of `θ_liq` with
`q_tot`, and environmental TKE, given:
 - `surf`, a `SurfaceModel`
 - `turbconv`, an `EDMF` model
 - `atmos`, an `AtmosModel`
 - `state`, state variables
 - `aux`, auxiliary variables
 - `zLL`, height of the lowest nodal level
"""
function subdomain_surface_values(
    surf::SurfaceModel,
    turbconv::EDMF{FT},
    atmos::AtmosModel{FT},
    state::Vars,
    aux::Vars,
    zLL::FT,
) where {FT}

    turbconv = atmos.turbconv
    N_up = n_updrafts(turbconv)
    gm = state
    # TODO: change to new_thermo_state
    ts = recover_thermo_state(atmos, state, aux)
    q = PhasePartition(ts)
    _cp_m = cp_m(ts)
    lv = latent_heat_vapor(ts)
    Π = exner(ts)
    ρ_inv = 1 / gm.ρ
    surface_scalar_coeff = turbconv.surface.scalar_coeff

    θ_liq_surface_flux = surf.shf / Π / _cp_m
    q_tot_surface_flux = surf.lhf / lv
    # these value should be given from the SurfaceFluxes.jl once it is merged
    oblength = turbconv.surface.obukhov_length
    ustar = turbconv.surface.ustar

    unstable = oblength < 0
    fact = unstable ? (1 - surf.ψϕ_stab * zLL / oblength)^(-FT(2 // 3)) : 1
    tke_fact = unstable ? cbrt(zLL / oblength * zLL / oblength) : 0
    ustar² = ustar^2
    θ_liq_cv = 4 * (θ_liq_surface_flux * θ_liq_surface_flux) / (ustar²) * fact
    q_tot_cv = 4 * (q_tot_surface_flux * q_tot_surface_flux) / (ustar²) * fact
    θ_liq_q_tot_cv =
        4 * (θ_liq_surface_flux * q_tot_surface_flux) / (ustar²) * fact
    tke = ustar² * (surf.κ_star² + tke_fact)

    a_up_surf = ntuple(i -> FT(surf.a / N_up), N_up)
    e_int = internal_energy(atmos, state, aux)
    ts_new = new_thermo_state(atmos, state, aux)
    θ_liq = liquid_ice_pottemp(ts_new)

    upd_θ_liq_surf = ntuple(N_up) do i
        θ_liq + surface_scalar_coeff[i] * sqrt(max(θ_liq_cv, 0))
    end

    ρq_tot = atmos.moisture isa DryModel ? FT(0) : gm.moisture.ρq_tot
    q_tot_up_surf = ntuple(N_up) do i
        ρq_tot * ρ_inv + surface_scalar_coeff[i] * sqrt(max(q_tot_cv, 0))
    end

    return (
        a_up_surf = a_up_surf,
        upd_θ_liq_surf = upd_θ_liq_surf,
        q_tot_up_surf = q_tot_up_surf,
        θ_liq_cv = θ_liq_cv,
        q_tot_cv = q_tot_cv,
        θ_liq_q_tot_cv = θ_liq_q_tot_cv,
        tke = tke,
    )
end;

"""
    percentile_bounds_mean_norm(
        low_percentile::FT,
        high_percentile::FT,
        n_samples::Int,
    ) where {FT <: AbstractFloat}

Returns the mean of all instances of a standard Gaussian random
variable that have a CDF higher than low_percentile and lower
than high_percentile, given a total of n_samples of the standard
Gaussian, given:
 - `low_percentile`, lower limit of the CDF
 - `high_percentile`, higher limit of the CDF
 - `n_samples`, the total number of samples drawn from the Gaussian
"""
function percentile_bounds_mean_norm(
    low_percentile::FT,
    high_percentile::FT,
    n_samples::Int,
) where {FT <: AbstractFloat}
    x = rand(Normal(), n_samples)
    xp_low = quantile(Normal(), low_percentile)
    xp_high = quantile(Normal(), high_percentile)
    filter!(y -> xp_low < y < xp_high, x)
    return Statistics.mean(x)
end

max_forks_repo_name: jm-c/CLIMA


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ISCSIPersistentVolumeSource represents an ISCSI disk. ISCSI volumes can only be mounted as read/write once. ISCSI volumes support ownership management and SELinux relabeling.

    IoK8sApiCoreV1ISCSIPersistentVolumeSource(;
        chapAuthDiscovery=nothing,
        chapAuthSession=nothing,
        fsType=nothing,
        initiatorName=nothing,
        iqn=nothing,
        iscsiInterface=nothing,
        lun=nothing,
        portals=nothing,
        readOnly=nothing,
        secretRef=nothing,
        targetPortal=nothing,
    )

    - chapAuthDiscovery::Bool : whether support iSCSI Discovery CHAP authentication
    - chapAuthSession::Bool : whether support iSCSI Session CHAP authentication
    - fsType::String : Filesystem type of the volume that you want to mount. Tip: Ensure that the filesystem type is supported by the host operating system. Examples: \"ext4\", \"xfs\", \"ntfs\". Implicitly inferred to be \"ext4\" if unspecified. More info: https://kubernetes.io/docs/concepts/storage/volumes#iscsi
    - initiatorName::String : Custom iSCSI Initiator Name. If initiatorName is specified with iscsiInterface simultaneously, new iSCSI interface <target portal>:<volume name> will be created for the connection.
    - iqn::String : Target iSCSI Qualified Name.
    - iscsiInterface::String : iSCSI Interface Name that uses an iSCSI transport. Defaults to 'default' (tcp).
    - lun::Int32 : iSCSI Target Lun number.
    - portals::Vector{String} : iSCSI Target Portal List. The Portal is either an IP or ip_addr:port if the port is other than default (typically TCP ports 860 and 3260).
    - readOnly::Bool : ReadOnly here will force the ReadOnly setting in VolumeMounts. Defaults to false.
    - secretRef::IoK8sApiCoreV1SecretReference : CHAP Secret for iSCSI target and initiator authentication
    - targetPortal::String : iSCSI Target Portal. The Portal is either an IP or ip_addr:port if the port is other than default (typically TCP ports 860 and 3260).
"""
mutable struct IoK8sApiCoreV1ISCSIPersistentVolumeSource <: SwaggerModel
    chapAuthDiscovery::Any # spec type: Union{ Nothing, Bool } # spec name: chapAuthDiscovery
    chapAuthSession::Any # spec type: Union{ Nothing, Bool } # spec name: chapAuthSession
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    initiatorName::Any # spec type: Union{ Nothing, String } # spec name: initiatorName
    iqn::Any # spec type: Union{ Nothing, String } # spec name: iqn
    iscsiInterface::Any # spec type: Union{ Nothing, String } # spec name: iscsiInterface
    lun::Any # spec type: Union{ Nothing, Int32 } # spec name: lun
    portals::Any # spec type: Union{ Nothing, Vector{String} } # spec name: portals
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    secretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: secretRef
    targetPortal::Any # spec type: Union{ Nothing, String } # spec name: targetPortal

    function IoK8sApiCoreV1ISCSIPersistentVolumeSource(;chapAuthDiscovery=nothing, chapAuthSession=nothing, fsType=nothing, initiatorName=nothing, iqn=nothing, iscsiInterface=nothing, lun=nothing, portals=nothing, readOnly=nothing, secretRef=nothing, targetPortal=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("chapAuthDiscovery"), chapAuthDiscovery)
        setfield!(o, Symbol("chapAuthDiscovery"), chapAuthDiscovery)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("chapAuthSession"), chapAuthSession)
        setfield!(o, Symbol("chapAuthSession"), chapAuthSession)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("initiatorName"), initiatorName)
        setfield!(o, Symbol("initiatorName"), initiatorName)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("iqn"), iqn)
        setfield!(o, Symbol("iqn"), iqn)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("iscsiInterface"), iscsiInterface)
        setfield!(o, Symbol("iscsiInterface"), iscsiInterface)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("lun"), lun)
        setfield!(o, Symbol("lun"), lun)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("portals"), portals)
        setfield!(o, Symbol("portals"), portals)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("secretRef"), secretRef)
        setfield!(o, Symbol("secretRef"), secretRef)
        validate_property(IoK8sApiCoreV1ISCSIPersistentVolumeSource, Symbol("targetPortal"), targetPortal)
        setfield!(o, Symbol("targetPortal"), targetPortal)
        o
    end
end # type IoK8sApiCoreV1ISCSIPersistentVolumeSource

const _property_map_IoK8sApiCoreV1ISCSIPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("chapAuthDiscovery")=>Symbol("chapAuthDiscovery"), Symbol("chapAuthSession")=>Symbol("chapAuthSession"), Symbol("fsType")=>Symbol("fsType"), Symbol("initiatorName")=>Symbol("initiatorName"), Symbol("iqn")=>Symbol("iqn"), Symbol("iscsiInterface")=>Symbol("iscsiInterface"), Symbol("lun")=>Symbol("lun"), Symbol("portals")=>Symbol("portals"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("secretRef")=>Symbol("secretRef"), Symbol("targetPortal")=>Symbol("targetPortal"))
const _property_types_IoK8sApiCoreV1ISCSIPersistentVolumeSource = Dict{Symbol,String}(Symbol("chapAuthDiscovery")=>"Bool", Symbol("chapAuthSession")=>"Bool", Symbol("fsType")=>"String", Symbol("initiatorName")=>"String", Symbol("iqn")=>"String", Symbol("iscsiInterface")=>"String", Symbol("lun")=>"Int32", Symbol("portals")=>"Vector{String}", Symbol("readOnly")=>"Bool", Symbol("secretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("targetPortal")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1ISCSIPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ISCSIPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ISCSIPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1ISCSIPersistentVolumeSource)
    (getproperty(o, Symbol("iqn")) === nothing) && (return false)
    (getproperty(o, Symbol("lun")) === nothing) && (return false)
    (getproperty(o, Symbol("targetPortal")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ISCSIPersistentVolumeSource }, name::Symbol, val)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e77948a6-1ecc-11eb-3285-7f06dcdb7cf5
using Pkg, DrWatson

# ╔═╡ e7799b9e-1ecc-11eb-311f-035de72c3603
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using StatsBase, Combinatorics, Plots
end

# ╔═╡ 6c38d2ce-1ecc-11eb-281d-0d13a832adea
md"## Listing 2.3"

# ╔═╡ fe1ee6a6-1ecc-11eb-1930-0301940bc5f7
	N = 10^3

# ╔═╡ e77a2a28-1ecc-11eb-0acb-5167936362f5
begin
	matchExists1(n) = 1 - prod([k/365 for k in 365:-1:365-n+1])
	matchExists2(n) = 1- factorial(365,365-big(n))/365^big(n)

	function bdEvent(n)
		birthdays = rand(1:365,n)
		dayCounts = counts(birthdays, 1:365)
		return maximum(dayCounts) > 1
	end

	probEst(n) = sum([bdEvent(n) for _ in 1:N])/N
end

# ╔═╡ e7885e54-1ecc-11eb-03d1-2580aad82a3d
begin
	xGrid = 1:50
	analyticSolution1 = [matchExists1(n) for n in xGrid]
	analyticSolution2 = [matchExists2(n) for n in xGrid]
	("Maximum error: $(round(maximum(abs.(analyticSolution1 - analyticSolution2)), digits=25))")
end

# ╔═╡ e788edce-1ecc-11eb-3682-b1bfc98646a0
begin
	mcEstimates = [probEst(n) for n in xGrid]

	plot(xGrid, analyticSolution1, c=:blue, label="Analytic solution")
	scatter!(xGrid, mcEstimates, c=:red, ms=6, msw=0, shape=:xcross,
		label="MC estimate", xlims=(0,50), ylims=(0, 1),
		xlabel="Number of people in room",
		ylabel="Probability of birthday match",
		legend=:topleft)
end

# ╔═╡ e7965b58-1ecc-11eb-2308-093c9f5a167e
md"## End of listing 2.3"

# ╔═╡ Cell order:
# ╟─6c38d2ce-1ecc-11eb-281d-0d13a832adea
# ╠═e77948a6-1ecc-11eb-3285-7f06dcdb7cf5
# ╠═e7799b9e-1ecc-11eb-311f-035de72c3603
# ╠═fe1ee6a6-1ecc-11eb-1930-0301940bc5f7
# ╠═e77a2a28-1ecc-11eb-0acb-5167936362f5
# ╠═e7885e54-1ecc-11eb-03d1-2580aad82a3d
# ╠═e788edce-1ecc-11eb-3682-b1bfc98646a0
# ╟─e7965b58-1ecc-11eb-2308-093c9f5a167e


text:
function pancam(filter::Symbol)
   filter == :L2 ? TSpec(PancamL2[:,1],PancamL2[:,2]) :
   filter == :L3 ? TSpec(PancamL3[:,1],PancamL3[:,2]) :
   filter == :L4 ? TSpec(PancamL4[:,1],PancamL4[:,2]) :
   filter == :L5 ? TSpec(PancamL5[:,1],PancamL5[:,2]) :
   filter == :L6 ? TSpec(PancamL6[:,1],PancamL6[:,2]) :
   filter == :L7 ? TSpec(PancamL7[:,1],PancamL7[:,2]) : nothing
end

const PancamL2 =
   [729.893842887474	   0.0;
    735.83864118896 	   0.0402542372881356;
    739.235668789809	   0.120762711864407;
    741.783439490446	   0.347457627118644;
    743.481953290871	   0.591101694915254;
    744.755838641189	   0.75635593220339;
    746.454352441614	   0.925847457627119;
    747.303609341826	   0.995762711864407;
    748.152866242038	   1.0;
    752.3991507431  	   0.976694915254237;
    755.796178343949	   0.976694915254237;
    756.645435244161	   0.955508474576271;
    757.494692144374	   0.938559322033898;
    757.91932059448	   0.915254237288136;
    758.768577494692	   0.864406779661017;
    759.617834394904	   0.805084745762712;
    760.891719745223	   0.722457627118644;
    761.740976645435	   0.574152542372881;
    763.014861995754	   0.461864406779661;
    764.713375796178	   0.288135593220339;
    767.261146496815	   0.154661016949153;
    768.95966029724 	   0.0911016949152542;
    771.082802547771	   0.0444915254237288;
    774.055201698514	   0.0233050847457627;
    776.602972399151	   0.0105932203389831;
    780.0           	   0.0]

const PancamL3 =
   [653.460721868365	   0.0;
    656.433121019108	   0.0211864406779661;
    658.980891719745	   0.0741525423728814;
    661.104033970276	   0.199152542372881;
    662.802547770701	   0.415254237288136;
    664.076433121019	   0.586864406779661;
    666.19957537155	   0.720338983050847;
    669.171974522293	   0.775423728813559;
    671.295116772824	   0.832627118644068;
    672.569002123142	   0.889830508474576;
    673.842887473461	   0.959745762711864;
    674.692144373673	   0.991525423728814;
    675.541401273885	   1.0;
    676.815286624204	   0.99364406779661;
    677.23991507431	   0.955508474576271;
    677.664543524416	   0.925847457627119;
    678.513800424629	   0.796610169491525;
    680.212314225053	   0.544491525423729;
    682.335456475584	   0.307203389830509;
    684.033970276009	   0.163135593220339;
    685.307855626327	   0.116525423728814;
    687.855626326964	   0.0550847457627119;
    689.129511677282	   0.0338983050847458;
    691.677282377919	   0.0148305084745763;
    695.074309978769	   0.0]

const PancamL4 =
   [580.000000000000    0.0;
    584.670912951168    0.0105932203389831;
    588.067940552017    0.0444915254237288;
    590.615711252654    0.152542372881356;
    593.588110403397    0.563559322033898;
    595.711252653928    0.930084745762712;
    596.985138004246    0.98728813559322;
    597.409766454352    0.997881355932203;
    598.259023354565    1.0;
    600.382165605096    0.991525423728814;
    602.929936305732    0.974576271186441;
    603.779193205945    0.974576271186441;
    604.628450106157    0.98728813559322;
    605.477707006369    0.98728813559322;
    606.751592356688    0.963983050847458;
    607.600849256900    0.930084745762712;
    608.874734607219    0.802966101694915;
    610.148619957537    0.629237288135593;
    611.847133757962    0.425847457627119;
    612.696390658174    0.298728813559322;
    613.970276008493    0.199152542372881;
    616.093418259023    0.103813559322034;
    618.216560509554    0.0593220338983051;
    619.915074309979    0.0338983050847458;
    622.462845010616    0.0169491525423729;
    625.010615711253    0.00635593220338983;
    630.106157112527    0.0]

const PancamL5 =
   [518.004246284501	   0.0;
    520.127388535032	   0.010593220339;
    521.40127388535	   0.025423728814;
    522.250530785563	   0.055084745763;
    524.7983014862	   0.415254237288;
    525.647558386412	   0.616525423729;
    526.496815286624	   0.690677966102;
    526.92144373673	   0.718220338983;
    529.044585987261	   0.737288135593;
    530.743099787686	   0.783898305085;
    531.592356687898	   0.824152542373;
    532.44161358811	   0.862288135593;
    533.290870488323	   0.904661016949;
    534.989384288747	   0.940677966102;
    537.112526539278	   0.959745762712;
    539.660297239915	   1.0;
    540.509554140127	   1.0;
    541.35881104034	   0.98093220339;
    541.783439490446	   0.96186440678;
    542.208067940552	   0.940677966102;
    542.632696390658	   0.887711864407;
    544.331210191083	   0.580508474576;
    545.605095541401	   0.341101694915;
    546.87898089172	   0.16313559322;
    548.152866242038	   0.091101694915;
    549.426751592357	   0.03813559322;
    550.700636942675	   0.01906779661;
    552.823779193206	   0.006355932203;
    560.042462845011	   0.0]

const PancamL6 =
   [450.063694267516	   0.0;
    455.15923566879	   0.004237288136;
    458.556263269639	   0.014830508475;
    459.405520169851	   0.01906779661;
    460.254777070064	   0.033898305085;
    461.528662420382	   0.080508474576;
    464.501061571125	   0.322033898305;
    465.350318471338	   0.411016949153;
    466.19957537155	   0.506355932203;
    467.473460721868	   0.563559322034;
    468.747346072187	   0.584745762712;
    470.021231422505	   0.622881355932;
    471.295116772824	   0.66313559322;
    472.144373673036	   0.701271186441;
    474.692144373673	   0.866525423729;
    476.390658174098	   0.911016949153;
    477.664543524416	   0.930084745763;
    481.486199575372	   0.930084745763;
    483.184713375796	   0.940677966102;
    484.883227176221	   0.963983050847;
    486.581740976645	   0.987288135593;
    488.28025477707	   1.0;
    489.978768577495	   0.991525423729;
    490.828025477707	   0.972457627119;
    491.677282377919	   0.947033898305;
    492.101910828025	   0.915254237288;
    492.951167728238	   0.877118644068;
    493.375796178344	   0.834745762712;
    493.80042462845	   0.786016949153;
    494.225053078556	   0.716101694915;
    495.498938428875	   0.516949152542;
    496.348195329087	   0.463983050847;
    498.046709129512	   0.28813559322;
    498.895966029724	   0.228813559322;
    500.594479830149	   0.177966101695;
    501.443736730361	   0.14406779661;
    503.142250530786	   0.08686440678;
    505.690021231423	   0.025423728814;
    509.936305732484	   0.008474576271;
    515.031847133758	   0.002118644068;
    524.7983014862	   0.0]

const PancamL7 =
   [380.0     	         0.0
    387.643312101911	   0.023305084746
    392.314225053079	   0.046610169492
    397.409766454352	   0.088983050847
    402.08067940552	   0.150423728814
    410.148619957537	   0.300847457627
    419.915074309979	   0.459745762712
    421.188959660297	   0.483050847458
    422.887473460722	   0.5
    428.407643312102	   0.578389830508
    433.07855626327	   0.677966101695
    436.475583864119	   0.777542372881
    439.023354564756	   0.817796610169
    440.72186836518	   0.83686440678
    447.091295116773	   0.932203389831
    450.488322717622	   1.0
    452.186836518047	   0.987288135593
    453.036093418259	   0.966101694915
    453.885350318471	   0.906779661017
    455.583864118896	   0.819915254237
    457.282377919321	   0.580508474576
    458.980891719745	   0.411016949153
    461.104033970276	   0.218220338983
    462.802547770701	   0.156779661017
    464.925690021231	   0.080508474576
    466.19957537155	   0.048728813559
    468.322717622081	   0.027542372881
    471.295116772824	   0.012711864407
    473.842887473461	   0.004237288136
    479.787685774947	   0.0]


text:
function init_3(t, fft_u0, fft_v0, A1, A2, matr, conjmatr,
           sigma, llambda, Ktaubis, epsilon, Ntaumm)

    h1u, h1v = ftau(t, fft_u0, fft_v0, A1, A2, matr, conjmatr, sigma, llambda)

    h1u_fft = fft(h1u, dims=1)
    h1v_fft = fft(h1v, dims=1)

    dtu0u = h1u_fft[0, :] / Ntaumm
    dtu0v = h1v_fft[0, :] / Ntaumm

    h1u_fft[0, :] = 0 * h1u_fft[0, :]
    h1v_fft[0, :] = 0 * h1v_fft[0, :]

    h1u = epsilon * ifft(h1u_fft / (1im * Ktaubis), dims=1)
    h1v = epsilon * ifft(h1v_fft / (1im * Ktaubis), dims=1)

    w1u = fft_u0 - h1u[0, :]
    w1v = fft_v0 - h1v[0, :]

    h1u, h1v = ftau(t, w1u, w1v, A1, A2, matr, conjmatr, sigma, llambda)

    h1u_fft = fft(h1u, dims=1)
    h1v_fft = fft(h1v, dims=1)

    dtu0u = h1u_fft[0, :] / Ntaumm
    dtu0v = h1v_fft[0, :] / Ntaumm

    h1u_fft[0, :] = 0 * h1u_fft[0, :]
    h1v_fft[0, :] = 0 * h1v_fft[0, :]

    h1u = epsilon * ifft(h1u_fft / (1im * Ktaubis), dims=1)
    h1v = epsilon * ifft(h1v_fft / (1im * Ktaubis), dims=1)

    champu1, champv1 = duftau(t, w1u, w1v, dtu0u, dtu0v,
                              A1, A2, matr, conjmatr, sigma, llambda, Ntaumm)

    champu2, champv2 = dtftau(t, w1u, w1v, A1, A2, matr, conjmatr, sigma, llambda, Ntaumm)

    champu = champu1 + champu2
    champv = champv1 + champv2
    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)
    champu_fft[0, :] = 0 * champu_fft[0, :]
    champv_fft[0, :] = 0 * champv_fft[0, :]
    dth1u = epsilon * ifft(champu_fft / (1im * Ktaubis), dims=1)
    dth1v = epsilon * ifft(champv_fft / (1im * Ktaubis), dims=1)

    h2u, h2v = ftau(t, w1u + h1u, w1v + h1v, A1, A2, matr, conjmatr, sigma,
                    llambda)
    h2u = h2u - dth1u
    h2v = h2v - dth1v
    h2u_fft = fft(h2u, dims=1)
    h2v_fft = fft(h2v, dims=1)
    h2u_fft[0, :] = 0 * h2u_fft[0, :]
    h2v_fft[0, :] = 0 * h2v_fft[0, :]
    h2u = epsilon * ifft(h2u_fft / (1im * Ktaubis), dims=1)
    h2v = epsilon * ifft(h2v_fft / (1im * Ktaubis), dims=1)

    fft_ubar = fft_u0 - h2u[0, :]
    fft_vbar = fft_v0 - h2v[0, :]

    C2u, C2v = C2(t, fft_ubar, fft_vbar, A1, A2, matr, conjmatr,
                  sigma, llambda, Ktaubis, epsilon, Ntaumm)

    fft_ug = fft_u0 - C2u[0, :]
    fft_vg = fft_v0 - C2v[0, :]

    return fft_ubar, fft_vbar, fft_ug, fft_vg

end


text:
"embedding function for the Gupta potential"
type GuptaEmbed <: SimpleFunction
    xi
end
evaluate(p::GuptaEmbed, r) = p.xi * sqrt(r)
evaluate_d(p::GuptaEmbed, r) = 0.5*p.xi ./ sqrt(r)

"""`GuptaPotential`:
    E_i = A ∑_{j ≠ i} v(r_ij) - ξ ∑_i √ ρ_i
        v(r_ij) = exp[ -p (r_ij/r0 - 1) ]
        ρ_i = ∑_{j ≠ i} exp[ -2q (r_ij / r0 - 1) ]
"""
GuptaPotential(A, xi, p, q, r0, TC::Type, TCargs...)  =
    EAMPotential( TC( SimpleExponential(A, p, r0), TCargs... ),      # V
                  TC( SimpleExponential(1.0, 2*q, r0), TCargs...),   # rho
                  GuptaEmbed( xi ) )                                 # embed

gupta_parameters = Dict(
   #          a       D       α      R₀      <<< Gupta notation
   "Cu" => (3.6147, 0.3446, 1.3921, 2.838),
   "Ag" => (4.0862, 0.3294, 1.3939, 3.096),
   "Au" => (4.0785, 0.4826, 1.6166, 3.004),
   "Ni" => (3.5238, 0.4279, 1.3917, 2.793)
)


text:
export load_experiment, load_shot_data, resume_from_checkpoint

"""
   load_experiment(sim_name::String)

   A function that uses the input variables to the experiment, passed by dictionaty
   `parsed_args`, generated via `ArgParse.jl` package, to determine which experiment data to
   load from `datadir(sim_name)`, where `sim_name` is the simulation name.

   In case there are multiple experiments with the same variable valus, the last experiment
   will be returned

    Args:
        parsed_args: A dictionary of input variables.
        return_path: A boolian indicating whether or not return the absolute path to the
            experiment.

   """

function load_experiment(
    parsed_args::Dict,
    keys_to_load::Vector{String};
    return_path = false,
    train = false,
)
    loaded_keys = Dict()
    sim_name = parsed_args["sim_name"]
    !isdir(datadir(sim_name)) && mkpath(datadir(sim_name))
    experiments = collect_results(
        datadir(sim_name),
        black_list = ["Params", "fval", "fval_eval", "opt"],
    )

    # Drop experiments that miss one or some of the input arguments
    if size(experiments, 1) > 0
        # dropmissing!(experiments)
        for (key, value) in parsed_args
            dropmissing!(experiments, key)
        end

        # Keep experiment with same variables
        for (key, value) in parsed_args
            if !train || key != "epoch"
                experiments = experiments[experiments[!, string(key)].==value, :]
            end
        end
    end

    if size(experiments, 1) < 1 && train
        @info "No saved experiments found"
        for key in keys_to_load
            loaded_keys[key] = nothing
        end
        return loaded_keys
    elseif size(experiments, 1) < 1 && !train
        @error "No saved experiments found with such input values"
    end

    @info (string(size(experiments, 1)) * " experiment(s) found — loading the most recent")
    if hasproperty(experiments, :epoch)
        experiment = experiments[partialsortperm(experiments.epoch, 1, rev = true), :]
        experiment_path = experiment[:path]
    else
        experiment_path = experiments[1, :path]
    end
    for key in keys_to_load
        loaded_keys[key] = wload(experiment_path)[key]
    end

    if return_path
        return loaded_keys, experiment_path
    else
        return loaded_keys
    end
end

function load_shot_data(parsed_args::Dict)

    # Upload data/sim_name
    sim_name = parsed_args["sim_name"]
    !isdir(datadir(sim_name)) && mkpath(datadir(sim_name))
    experiments = collect_results(datadir(sim_name), black_list = ["dm", "lin_data", "rtm"])

    # Drop experiments that miss one or some of the input arguments
    if size(experiments, 1) > 0
        dropmissing!(experiments)

        # Keep experiment with same variables
        for (key, value) in parsed_args
            experiments = experiments[experiments[!, string(key)].==value, :]
        end
    end

    if size(experiments, 1) < 1
        @info "No saved experiments found"
        return nothing, nothing, nothing
    end

    @info (string(size(experiments, 1)) * " experiment(s) found — loading the most recent")
    experiment_path = experiments[1, :path]
    dm = wload(experiment_path)["dm"]
    lin_data = wload(experiment_path)["lin_data"]
    rtm = wload(experiment_path)["rtm"]
    return dm, lin_data, rtm
end

function resume_from_checkpoint(args, keys_to_load::Vector{String})
    if args["resume"]
        delete!(args, "resume")
        loaded_keys = load_experiment(args, keys_to_load; train = true)
        if loaded_keys["epoch"] == nothing
            loaded_keys["epoch"] = 1
        else
            loaded_keys["epoch"] = loaded_keys["epoch"] + 1
            @info "Resuming training from epoch " * string(loaded_keys["epoch"])
        end
    else
        loaded_keys = Dict()
        for key in keys_to_load
            loaded_keys[key] = nothing
        end
        loaded_keys["epoch"] = 1
    end
    return loaded_keys
end


text:
using TranscodingStreams: TranscodingStream, Error, Memory,
    test_roundtrip_fileio, test_roundtrip_transcode
using TranscodingStreams

text = b"""
Lorem ipsum dolor sit amet, consectetur adipiscing elit. Aenean sollicitudin
mauris non nisi consectetur, a dapibus urna pretium. Vestibulum non posuere
erat. Donec luctus a turpis eget aliquet. Cras tristique iaculis ex, eu
malesuada sem interdum sed. Vestibulum ante ipsum primis in faucibus orci luctus
et ultrices posuere cubilia Curae; Etiam volutpat, risus nec gravida ultricies,
erat ex bibendum ipsum, sed varius ipsum ipsum vitae dui.
"""

@testset "Transcoding" begin
    @test LZ4FrameCompressorStream <: TranscodingStream
    @test LZ4FrameDecompressorStream <: TranscodingStream

    compressed = transcode(LZ4FrameCompressor, text)
    @test sizeof(compressed) < sizeof(text)

    decompressed = transcode(LZ4FrameDecompressor, compressed)
    @test sizeof(decompressed) > sizeof(compressed)
    @test decompressed == Vector{UInt8}(text)

    test_roundtrip_fileio(LZ4FrameCompressor, LZ4FrameDecompressor)
    test_roundtrip_transcode(LZ4FrameCompressor, LZ4FrameDecompressor)

    file = IOBuffer(text)
    stream = LZ4FrameDecompressorStream(LZ4FrameCompressorStream(file))
    flush(stream)

    @test hash(read(stream)) == hash(text)
    close(stream)
    close(file)

    file = IOBuffer(text)
    stream = LZ4FrameDecompressorStream(LZ4FrameCompressorStream(file; blocksizeid = max64KB))
    flush(stream)

    @test hash(read(stream)) == hash(text)
    close(stream)
    close(file)

    file = IOBuffer("")
    stream = LZ4FrameDecompressorStream(LZ4FrameCompressorStream(file))
    flush(stream)

    @test hash(read(stream)) == hash(b"")
    close(stream)
    close(file)
end

@testset "Errors" begin
    input = Memory(Vector{UInt8}(text))
    output = Memory(Vector{UInt8}(undef, 1280))
    not_initialized = LZ4FrameCompressor()
    @test TranscodingStreams.startproc(not_initialized, :read, Error()) == :error
    @test TranscodingStreams.process(not_initialized, input, output, Error()) == (0, 0, :error)

    compressed = transcode(LZ4FrameCompressor, Vector{UInt8}(text))
    corrupted = copy(compressed)
    corrupted[1] = 0x00
    file = IOBuffer(corrupted)
    stream = LZ4FrameDecompressorStream(file)
    @test_throws CodecLz4.LZ4Exception read(stream)
    @test_throws ArgumentError read(stream)

    output = Memory(Vector{UInt8}(undef, 1))
    compressor = LZ4FrameCompressor()
    @test_nowarn TranscodingStreams.initialize(compressor)
    @test TranscodingStreams.startproc(compressor, :read, Error()) == :ok
    err = Error()
    @test TranscodingStreams.process(compressor, input, output, err) == (0, 0, :error)
    @test err[].msg == "Output buffer too small for header."
    @test_nowarn TranscodingStreams.finalize(compressor)

    codec = LZ4FrameDecompressor()
    @test_throws CodecLz4.LZ4Exception transcode(codec, "not properly formatted")
    @test_nowarn TranscodingStreams.finalize(codec)
end

@testset "keywords" begin
    compressor = LZ4FrameCompressor(
        blocksizeid = max64KB,
        blockmode = block_independent,
        contentchecksum = true,
        blockchecksum = true,
        frametype = skippable_frame,
        contentsize = 100,
        compressionlevel=5,
        autoflush = true
        )

    prefs = compressor.prefs[]
    @test prefs.compressionLevel == Cint(5)
    @test prefs.autoFlush == Cuint(1)
    @test prefs.reserved == (Cuint(0), Cuint(0), Cuint(0), Cuint(0))

    frame = prefs.frameInfo
    @test frame.blockSizeID == Cuint(4)
    @test frame.blockMode == Cuint(1)
    @test frame.contentChecksumFlag == Cuint(1)
    @test frame.frameType == Cuint(1)
    @test frame.contentSize == Culonglong(100)
    @test frame.blockChecksumFlag == Cuint(1)
end

@testset "dst size fix" begin
    teststring = rand(UInt8, 800000)
    io = IOBuffer(teststring)
    stream = LZ4FrameCompressorStream(io)
    result = read(stream)
    @test_nowarn close(stream)
    close(io)

    io = IOBuffer(result)
    stream = LZ4FrameDecompressorStream(io)
    result = read(stream)
    @test result == teststring
    @test_nowarn close(stream)
    close(io)
end


text:
export random_forest

import ..RandomForestModel

using mlpack._Internal.io

import mlpack_jll
const random_forestLibrary = mlpack_jll.libmlpack_julia_random_forest

# Call the C binding of the mlpack random_forest binding.
function random_forest_mlpackMain()
  success = ccall((:random_forest, random_forestLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module random_forest_internal
  import ..random_forestLibrary

import ...RandomForestModel

# Get the value of a model pointer parameter of type RandomForestModel.
function IOGetParamRandomForestModel(paramName::String)::RandomForestModel
  RandomForestModel(ccall((:IO_GetParamRandomForestModelPtr, random_forestLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type RandomForestModel.
function IOSetParamRandomForestModel(paramName::String, model::RandomForestModel)
  ccall((:IO_SetParamRandomForestModelPtr, random_forestLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeRandomForestModel(stream::IO, model::RandomForestModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeRandomForestModelPtr, random_forestLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeRandomForestModel(stream::IO)::RandomForestModel
  buffer = read(stream)
  RandomForestModel(ccall((:DeserializeRandomForestModelPtr, random_forestLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    random_forest(; [input_model, labels, maximum_depth, minimum_gain_split, minimum_leaf_size, num_trees, print_training_accuracy, seed, subspace_dim, test, test_labels, training, verbose])

This program is an implementation of the standard random forest classification
algorithm by Leo Breiman.  A random forest can be trained and saved for later
use, or a random forest may be loaded and predictions or class probabilities for
points may be generated.

The training set and associated labels are specified with the `training` and
`labels` parameters, respectively.  The labels should be in the range [0,
num_classes - 1]. Optionally, if `labels` is not specified, the labels are
assumed to be the last dimension of the training dataset.

When a model is trained, the `output_model` output parameter may be used to save
the trained model.  A model may be loaded for predictions with the
`input_model`parameter. The `input_model` parameter may not be specified when
the `training` parameter is specified.  The `minimum_leaf_size` parameter
specifies the minimum number of training points that must fall into each leaf
for it to be split.  The `num_trees` controls the number of trees in the random
forest.  The `minimum_gain_split` parameter controls the minimum required gain
for a decision tree node to split.  Larger values will force higher-confidence
splits.  The `maximum_depth` parameter specifies the maximum depth of the tree.
The `subspace_dim` parameter is used to control the number of random dimensions
chosen for an individual node's split.  If `print_training_accuracy` is
specified, the calculated accuracy on the training set will be printed.

Test data may be specified with the `test` parameter, and if performance
measures are desired for that test set, labels for the test points may be
specified with the `test_labels` parameter.  Predictions for each test point may
be saved via the `predictions`output parameter.  Class probabilities for each
prediction may be saved with the `probabilities` output parameter.

For example, to train a random forest with a minimum leaf size of 20 using 10
trees on the dataset contained in `data`with labels `labels`, saving the output
random forest to `rf_model` and printing the training error, one could call

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> labels = CSV.read("labels.csv"; type=Int)
julia> rf_model, _, _ = random_forest(labels=labels,
            minimum_leaf_size=20, num_trees=10, print_training_accuracy=1,
            training=data)
```

Then, to use that model to classify points in `test_set` and print the test
error given the labels `test_labels` using that model, while saving the
predictions for each point to `predictions`, one could call

```julia
julia> using CSV
julia> test_set = CSV.read("test_set.csv")
julia> test_labels = CSV.read("test_labels.csv"; type=Int)
julia> _, predictions, _ = random_forest(input_model=rf_model,
            test=test_set, test_labels=test_labels)
```

# Arguments

 - `input_model::RandomForestModel`: Pre-trained random forest to use for
      classification.
 - `labels::Array{Int, 1}`: Labels for training dataset.
 - `maximum_depth::Int`: Maximum depth of the tree (0 means no limit).
      Default value `0`.

 - `minimum_gain_split::Float64`: Minimum gain needed to make a split when
      building a tree.  Default value `0`.

 - `minimum_leaf_size::Int`: Minimum number of points in each leaf node.
      Default value `1`.

 - `num_trees::Int`: Number of trees in the random forest.  Default value
      `10`.

 - `print_training_accuracy::Bool`: If set, then the accuracy of the model
      on the training set will be predicted (verbose must also be specified).
      Default value `false`.

 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `subspace_dim::Int`: Dimensionality of random subspace to use for each
      split.  '0' will autoselect the square root of data dimensionality.
      Default value `0`.

 - `test::Array{Float64, 2}`: Test dataset to produce predictions for.
 - `test_labels::Array{Int, 1}`: Test dataset labels, if accuracy
      calculation is desired.
 - `training::Array{Float64, 2}`: Training dataset.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output_model::RandomForestModel`: Model to save trained random forest
      to.
 - `predictions::Array{Int, 1}`: Predicted classes for each point in the
      test set.
 - `probabilities::Array{Float64, 2}`: Predicted class probabilities for
      each point in the test set.

"""
function random_forest(;
                       input_model::Union{RandomForestModel, Missing} = missing,
                       labels = missing,
                       maximum_depth::Union{Int, Missing} = missing,
                       minimum_gain_split::Union{Float64, Missing} = missing,
                       minimum_leaf_size::Union{Int, Missing} = missing,
                       num_trees::Union{Int, Missing} = missing,
                       print_training_accuracy::Union{Bool, Missing} = missing,
                       seed::Union{Int, Missing} = missing,
                       subspace_dim::Union{Int, Missing} = missing,
                       test = missing,
                       test_labels = missing,
                       training = missing,
                       verbose::Union{Bool, Missing} = missing,
                       points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, random_forestLibrary), Nothing, ());

  IORestoreSettings("Random forests")

  # Process each input argument before calling mlpackMain().
  if !ismissing(input_model)
    random_forest_internal.IOSetParamRandomForestModel("input_model", convert(RandomForestModel, input_model))
  end
  if !ismissing(labels)
    IOSetParamURow("labels", labels)
  end
  if !ismissing(maximum_depth)
    IOSetParam("maximum_depth", convert(Int, maximum_depth))
  end
  if !ismissing(minimum_gain_split)
    IOSetParam("minimum_gain_split", convert(Float64, minimum_gain_split))
  end
  if !ismissing(minimum_leaf_size)
    IOSetParam("minimum_leaf_size", convert(Int, minimum_leaf_size))
  end
  if !ismissing(num_trees)
    IOSetParam("num_trees", convert(Int, num_trees))
  end
  if !ismissing(print_training_accuracy)
    IOSetParam("print_training_accuracy", convert(Bool, print_training_accuracy))
  end
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(subspace_dim)
    IOSetParam("subspace_dim", convert(Int, subspace_dim))
  end
  if !ismissing(test)
    IOSetParamMat("test", test, points_are_rows)
  end
  if !ismissing(test_labels)
    IOSetParamURow("test_labels", test_labels)
  end
  if !ismissing(training)
    IOSetParamMat("training", training, points_are_rows)
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output_model")
  IOSetPassed("predictions")
  IOSetPassed("probabilities")
  # Call the program.
  random_forest_mlpackMain()

  return random_forest_internal.IOGetParamRandomForestModel("output_model"),
         IOGetParamURow("predictions"),
         IOGetParamMat("probabilities", points_are_rows)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Libbpf"
version = v"0.4"

# Collection of sources required to build libbpf
sources = [
    ArchiveSource("https://github.com/libbpf/libbpf/archive/v$version.tar.gz",
                  "21cbee4df093e7fd29e76ed429650d3f3abe3a893f35e346ab9bc3484f6e68c0"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libbpf-*/src/
for patchfile in $WORKSPACE/srcdir/patches/*; do
    atomic_patch -p1 $patchfile
done
make -j${nproc}
make install PREFIX=$prefix
cp -R ../include/uapi "$prefix/include/uapi"
"""

# Only build for Linux
platforms = filter!(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libbpf", :libbpf),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Elfutils_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
using DatagenCopulaBased
using Test
using Distributions
using HypothesisTests
using StatsBase
using Random
using Distributed
using LinearAlgebra
using SharedArrays
using Combinatorics
using SpecialFunctions
using StableRNGs

import DatagenCopulaBased: rand2cop, fncopulagen
import DatagenCopulaBased: logseriescdf, levyel, tiltedlevygen
import DatagenCopulaBased: Ginv, InvlaJ, sampleInvlaJ, elInvlaF, nestedfrankgen
import DatagenCopulaBased: testθ, useρ, useτ, testbivθ, usebivρ, getθ4arch
import DatagenCopulaBased: τ2λ, moρ2τ, norm2unifind
import DatagenCopulaBased: dilog, Debye, frankτ2θ, τ2θ, AMHτ2θ, Ccl, Cg
import DatagenCopulaBased: gumbelθ2ρ, claytonθ2ρ, gumbelρ2θ, claytonρ2θ, frankρ2θ, ρ2θ, AMHρ2θ
import DatagenCopulaBased: mocopula, arch_gen
import DatagenCopulaBased: nestedcopulag
import DatagenCopulaBased: meanΣ, frechet, mean_outer, parameters, are_parameters_good, Σ_theor
import DatagenCopulaBased: getcors_advanced
import DatagenCopulaBased: random_unit_normal_vector
import DatagenCopulaBased: frechet_el!, frechet_el2!, mocopula_el

# axiliary tests
include("tailtest.jl")
include("univdatagentests.jl")
include("marg_cor_tests.jl")

# tests particular copulas generators
include("eliptic_fr_mo_test.jl")
include("archcopulatests.jl")
include("nestedarchcoptest.jl")
include("chaincopulastests.jl")

# test transforming of data by introducing higher order corss-correlations
include("higher_order_cors_tests.jl")


text:
# This file is a part of SimilaritySearch.jl

export Performance, probe, recallscore, timedsearchbatch, macrorecall

"""
    recallscore(gold::Set, res)

Compute recall and precision scores from the result sets.
"""
function recallscore(gold, res)::Float64
    length(intersect(_convert_as_set(gold), _convert_as_set(res))) / length(gold)
end

_convert_as_set(a::Set) = a
_convert_as_set(a::AbstractVector) = Set(a)
_convert_as_set(a::KnnResult) = Set(a.id)

"""
    macrorecall(goldI::AbstractMatrix, resI::AbstractMatrix, k=size(goldI, 1))::Float64

Computes the macro recall score using goldI as gold standard and resI as predictions;
it expects that matrices of integers (identifiers). If `k` is given, then the results are cut to first `k`.
"""
function macrorecall(goldI::AbstractMatrix, resI::AbstractMatrix, k=size(goldI, 1))::Float64
    @assert size(goldI) == size(resI)
    n = size(goldI, 2)
    s = 0.0
    for i in 1:n
        g = view(goldI, 1:k, i)
        r = view(resI, 1:k, i)
        s += recallscore(g, r)
    end

    s / n
end

"""
    macrorecall(goldlist::AbstractVector, reslist::AbstractVector)::Float64

Computes the macro recall score using sets of results (KnnResult objects or vectors of itegers).
"""
function macrorecall(goldlist::AbstractVector, reslist::AbstractVector)::Float64
    @assert size(goldlist) == size(reslist)
    s = 0.0
    n = length(goldlist)
    for i in 1:n
        g = goldlist[i]
        r = reslist[i]
        s += recallscore(g, r)
    end

    s / n
end

"""
    timedsearchbatch(index, Q, ksearch::Integer; parallel=false)

Computes the K nearest neigbors of each object in `Q` and returns two matrices, and the average search time in seconds.
"""
function timedsearchbatch(index, Q, ksearch::Integer; parallel=false, pools=getpools(index))
    m = length(Q)
    I = zeros(Int32, ksearch, m)
    D = Matrix{Float32}(undef, ksearch, m)
    t = @elapsed (I, D = searchbatch(index, Q, I, D; parallel, pools))
    I, D, t/length(Q)
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
# # Parallel evaluation and simple statistics
#
# In this example we will be looking at a way to evaluate several times a given problem
# by exploiting a distributed computing strategy known as **embarrasingly parallel**
# execution.
#
# The idea is very simple: compute the same problem several times, in parallel, and
# collect the results at the end.
# This should work well with metaheuristics due to the following properties:
#
# - Metaheuristics are _self-contained,_ meaning that they possess their own RNG and their own search space.
# - Due to the fact that solutions are searched with a random component, independent runs should be able to give independent results.
#
# With this in mind, we can use the fantastic parallel capabilities of `Julia` and
# show how easy it is to obtain some simple, but meaningful statistics for a given
# optimization problem.
#
# With this we want to also show that this is the **right way** to use this package.
# We should _always_ run several independent solvers to obtain good statistics.
# We must **never** settle with just one run from the solver.
#
# ## Setting up
#
# In order to use the multithreading capabilities of `Julia`, we need to launch the
# REPL using all the threads we can afford.
#
# For instance, the computer where this example is currently being executed has
# 8 physical cores, so I should launch my `Julia` REPL
# like so
#
# ```shell
# JULIA_NUM_THREADS=8 julia
# ```
#
# This will launch the REPL with 8 cores ready to be used. We can check this by entering
# the following in the REPL,
#
# ```julia
# julia> using Base.Threads
# julia> nthreads()
# ```
#
# making sure that this matches the same number introduced with the variable
# `JULIA_NUM_THREADS`.
#
# You should replace the number of threads, in this case 8, with the number of threads
# that you wish.
#
# ## Set up
#
# First, let us import all the necessary modules we will need for this example.

using BioMetaheuristics
using BioMetaheuristics.TestFunctions
using Base.Threads
using Statistics
using BenchmarkTools

# Now, let us create some setup functions.
# First, we will be using the `optimize` interface to solve the [`Levy`](@ref)
# optimization problem. This function just
# wraps the `optimize` function and returns the solution vector from the problem.
#
# We will be using Particle Swarm Optimization (PSO), with 20,000 iterations.
# It is very important to note that no specific RNG is being used.
# This tells the `optimize` function to always use a new, truly randomly seeded RNG.

function solve_problem(dim)
    val = optimize(Levy(), zeros(dim), [-35.0, 35.0], PSO(); iters=20_000)

    return val.x
end;

# ### Some benchmarks
#
# Before we continue, I would like to show a benchmark of the PSO implementation.
# We will use the fantastic [`BenchmarkTools.jl`](https://github.com/JuliaCI/BenchmarkTools.jl)
# to check the time it takes to finish a single optimization problem.
#
# Here, we will solve a 50-dimensional Levy optimization problem. We are not really
# interested if the solution is accurate, just the timing. We will deal with accuracy
# in the next subsection.

b = @benchmark solve_problem(50);
b

# Not bad at all. If it takes between 3-4 seconds to finish I think this is fairly
# good.
# Furthermore, it means that we can take advantage of how little it takes to run
# in order to really push the number of solvers to execute.
#
# ## Parallel evaluation
#
# We now turn to the main event. We will launch several solvers using the previously
# created wrapper function in parallel. First off, we need another function to do so.
#
# This function will use the `@threads` macro that enables a for-loop to run in parallel.
# The function has a matrix called `solutions` where we will be storing the solution
# vectors from each independent run.
#
# We will call the solver in parallel and each time we will save the solutions vector.
# We will then compute the mean, median and standard deviation from all the independent
# runs and report that as the _true_ value, and compare it to the _ground truth._

function parallel_compute(;dim=50, total_iterations=50)
    solutions = Matrix{Float64}(undef, dim, total_iterations)

    @threads for i in 1:total_iterations
        solutions[:, i] = solve_problem(dim)
    end

    mean_result = mean(solutions; dims=2)
    median_result = median(solutions; dims=2)
    std_result = std(solutions; dims=2)

    return (mean_result, median_result, std_result)
end;

# We are ready to go. We will do 200 iterations in order to have good statistics.
#
# !!! warning
#     Depending on your hardware this might take a long time. If that is the case
#     stop the process and reduce the number of iterations.
#
# !!! note
#     Every time you do this exercise, new results will appear. This is normal for
#     metaheuristics, and it is desired. We don't want reproducibility here, we aim
#     for _statistical significance._

μ, med, σ = parallel_compute(;total_iterations=200);

# We can now check the results

# Here is the mean value

μ'

# Here is the median obtained

med'

# And finally, the standard deviation of the solution

σ'

# Let us print the evaluation of the function on these results.
# First, the mean value.

println(evaluate(Levy(), μ))

# Then, the median obtained.

println(evaluate(Levy(), med))

# ## Conclusions
#
# Recall that for the Levy function the solution is
#
# ```math
# \mathbf{x^*} = (1, \dots, 1)
# ```
#
# and it should evaluate to
#
# ```math
# f(\mathbf{x^*}) = 0
# ```
#
# It it safe to say that we have obtained those results, but we are _statistically sure_
# that these results are the _true_ results, because we have several independent runs
# for the answer.
#
# ### On the use of the median
#
# Notice that most of the time, the median result actually performs better than the mean
# result.
# This is somewhat expected, as the median is not really affected by poorly performant runs.
# Sometimes, independent runs will not give good results and when taken into account,
# the mean value will suffer greatly from this.
# The median, however, as a measure of the "value in the middle", will not be so affected
# and a better estimation is expected.
#
# We should always report both results, and conclude that whenever the final function
# evaluation is approximately the same with both the mean and the median, the
# solution is statistically acceptable for that particular problem.
#


text:
using ColorTypes, FixedPointNumbers

# Color space transformation
function rgb_to_ycbcr(pic::Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}},2})::Matrix{YCbCr{Float32}}
    return YCbCr.(pic)
end

function ycbcr_to_rgb(pic::Array{ColorTypes.YCbCr{Float32},2})::Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8,8}},2}
    return RGB.(pic)
end

function split_channels(pic::Array{ColorTypes.YCbCr{Float32},2})::Tuple{Matrix{Float32},Matrix{Float32},Matrix{Float32}}
    ys = []
    cbs = []
    crs = []
    ys = map(x -> x.y, pic)
    cbs = map(x -> x.cb, pic)
    crs = map(x -> x.cr, pic)
    return (ys, cbs, crs)
end

function join_channels(x::Tuple{Matrix{Float32},Matrix{Float32},Matrix{Float32}})::Array{ColorTypes.YCbCr{Float32},2}
    retimg = Matrix{Tuple{Float32,Float32,Float32}}(undef, size(x[1]))
    for i=1:size(x[1])[1], j=1:size(x[1])[2]
        retimg[i,j] = (x[1][i,j], x[2][i,j], x[3][i,j])
    end
    return map(x -> YCbCr(x...), retimg)
end