text:
using DINA
import CSV
using Test

@testset "DINA.jl" begin

    tbl = get_dina(1980)
    @test tbl isa DINA.StatFiles.StatFile

    @test sort(dina_years()) == [1962; 1964; 1966:2019]

    var = [:fiinc, :fninc, :ownermort, :ownerhome, :rentalmort, :rentalhome]
    df = dina_quantile_panel(var, :fiinc, 10)
    @test df isa DINA.DataFrames.DataFrame
    @test size(df) == (1530, 11)

    filename = joinpath(@__DIR__(), "dina-aggregated.csv")
    df |> CSV.write(filename)

    @show filename
end


text:
@testset "sequence_length(::Vector{FASTX.FASTA.Record})" begin
    fasta_file = joinpath(@__DIR__, "data", "test.fasta")
    records = FASTX.FASTA.Record[]

    open(FASTA.Reader, fasta_file) do reader
        for record in reader
            push!(records, record)
        end
    end

    @test Escape.sequence_length(records) == 14
end

@testset "consensus_sequence(::Vector{FASTX.FASTA.Record})" begin
    fasta_file = joinpath(@__DIR__, "data", "test.fasta")
    records = FASTX.FASTA.Record[]

    open(FASTA.Reader, fasta_file) do reader
        for record in reader
            push!(records, record)
        end
    end

    @test Escape.consensus_sequence(records) == "MSARASVMGSRASV"
end

@testset "epitope_prediction(::String)" begin
    peptide = "ARGDEFFPE"
    alleles = Escape.valid_alleles()
    prediction = Escape.epitope_prediction(peptide, alleles[1:10])

    @test prediction isa DataFrame
end

@testset "epitope_prediction(::AbstractHLAData)" begin
    hla_types = [HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01")),
                 HLAType(parse_allele("A01", "A01", "B01", "B01", "C01", "C01"))]
    fasta_file = joinpath(@__DIR__, "data", "test.fasta")
    data = HLAData("test", fasta_file, hla_types, missing, missing)

    df = @suppress Escape.epitope_prediction(data, rank_threshold = 100)

    @test nrow(df) > 0
end

@testset "epitope_feature_matrix(::AbstractHLAData)" begin
    hla_types = [HLAType(parse_allele("A01", "A02", "B05", "B06", "C03", "C07")),
                 HLAType(parse_allele("A01", "A03", "B03", "B05", "C03", "C07")),
                 HLAType(parse_allele("A01", "A04", "B03", "B05", "C03", "C07")),
                 HLAType(parse_allele("A02", "A04", "B03", "B04", "C03", "C07")),
                 HLAType(parse_allele("A02", "A04", "B03", "B04", "C03", "C07"))]

    fasta_file = joinpath(@__DIR__, "data", "test.fasta")
    hla_data = HLAData("test", fasta_file, hla_types, missing, missing)

    @test Escape.epitope_feature_matrix(hla_data) isa Matrix
end

ext: jl


text:
using StaticArrays: SA

const atomic_symbols = SA["H","He","Li","Be","B","C","N","O","F","Ne",
        "Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca",
        "Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn",
        "Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr",
        "Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd","In","Sn",
        "Sb","Te","I","Xe","Cs","Ba","La","Ce","Pr","Nd",
        "Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb",
        "Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg",
        "Tl","Pb","Bi","Po","At","Rn","Fr","Ra","Ac","Th",
        "Pa","U"] #,"Np","Pu","Am"] # we don't have partition funcs for these last three

const Natoms = UInt8(length(atomic_symbols))
const atomic_numbers = Dict(atomic_symbols .=> UInt8.(1:Natoms))

#in grams
const atomic_masses = SA[
     1.008,4.003,6.941,9.012,10.81,12.01,14.01,16.00,19.00,20.18,
     22.99,24.31,26.98,28.08,30.97,32.06,35.45,39.95,39.10,40.08,
     44.96,47.90,50.94,52.00,54.94,55.85,58.93,58.71,63.55,65.37,
     69.72,72.59,74.92,78.96,79.90,83.80,85.47,87.62,88.91,91.22,
     92.91,95.94,98.91,101.1,102.9,106.4,107.9,112.4,114.8,118.7,
     121.8,127.6,126.9,131.3,132.9,137.3,138.9,140.1,140.9,144.2,
     145.0,150.4,152.0,157.3,158.9,162.5,164.9,167.3,168.9,173.0,
     175.0,178.5,181.0,183.9,186.2,190.2,192.2,195.1,197.0,200.6,
     204.4,207.2,209.0,210.0,210.0,222.0,223.0,226.0,227.0,232.0,
     231.0,238.0].* amu_cgs #,237.0,244.0,243.0]

#solar/meteoritic abundances per Asplund et al. (2009, Ann. Rev. Ast. Ap., 47, 481).
const solar_abundances = SA[
      12.00,10.93, 1.05, 1.38, 2.70, 8.43, 7.83, 8.69, 4.56, 7.93,
      6.24, 7.60, 6.45, 7.51, 5.41, 7.12, 5.50, 6.40, 5.03, 6.34,
      3.15, 4.95, 3.93, 5.64, 5.43, 7.50, 4.99, 6.22, 4.19, 4.56,
      3.04, 3.65, 2.30, 3.34, 2.54, 3.25, 2.52, 2.87, 2.21, 2.58,
      1.46, 1.88,-5.00, 1.75, 0.91, 1.57, 0.94, 1.71, 0.80, 2.04,
      1.01, 2.18, 1.55, 2.24, 1.08, 2.18, 1.10, 1.58, 0.72, 1.42,
     -5.00, 0.96, 0.52, 1.07, 0.30, 1.10, 0.48, 0.92, 0.10, 0.84,
      0.10, 0.85,-0.12, 0.85, 0.26, 1.40, 1.38, 1.62, 0.92, 1.17,
      0.90, 1.75, 0.65,-5.00,-5.00,-5.00,-5.00,-5.00,-5.00, 0.02,
      -5.00,-0.54]#,-5.00,-5.00,-5.00]


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
name = "libvorbis"
version = v"1.3.6"

# Collection of sources required to build libvorbis
sources = [
    ArchiveSource("https://downloads.xiph.org/releases/vorbis/libvorbis-$(version).tar.xz",
                  "af00bb5a784e7c9e69f56823de4637c350643deedaf333d0fa86ecdba6fcb415"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libvorbis-*

# Patch in explicit dependency on `stdint.h` to fix libogg
# This will be fixed in a future libogg release:
# https://github.com/xiph/ogg/commit/c8fca6b4a02d695b1ceea39b330d4406001c03ed
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/stdint.patch"

# Force `./configure` to repent from using `-ffast-math`
sed -ie 's/-ffast-math//g' ./configure

./configure --prefix=$prefix --host=$target --build=${MACHTYPE} --disable-static
make -j${nproc}
make install
# Remove large docs directory
rm -r "${prefix}/share/doc"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# The products that we will ensure are always built
products = [
    LibraryProduct("libvorbis", :libvorbis),
    LibraryProduct("libvorbisenc", :libvorbisenc),
    LibraryProduct("libvorbisfile", :libvorbisfile),
]

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

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
function exp_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,exp(xL),xU,exp(xU)),dline_seg(x,xL,exp(xL),xU,exp(xU),exp(x))
end
function exp_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return exp(x),exp(x)
end
function exp(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = exp_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = exp_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = exp_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = exp_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = exp_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = exp_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, exp(x.Intv),x.cnst, x.IntvBox, x.xref)
end

function exp2_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,exp2(xL),xU,exp2(xU)),dline_seg(x,xL,exp2(xL),xU,exp2(xU),exp2(x)*log(2))
end
function exp2_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return exp2(x),exp2(x)*log(2)
end
function exp2(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = exp2_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = exp2_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = exp2_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = exp2_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = exp2_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = exp2_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, exp2(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function exp10_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,exp10(xL),xU,exp10(xU)),dline_seg(x,xL,exp10(xL),xU,exp10(xU),exp10(x)*log(10))
end
function exp10_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return exp10(x),exp10(x)*log(10)
end
function exp10(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = exp10_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = exp10_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = exp10_cc(x.cv,x.Intv.lo,x.Intv.hi)
    gcv1::T,gdcv1::T = exp10_cv(x.cv,x.Intv.lo,x.Intv.hi)
    gcc2::T,gdcc2::T = exp10_cc(x.cc,x.Intv.lo,x.Intv.hi)
    gcv2::T,gdcv2::T = exp10_cv(x.cc,x.Intv.lo,x.Intv.hi)
    cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, exp10(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function sqrt_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return sqrt(x),one(T)/(2*sqrt(x))
end
function sqrt_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,sqrt(xL),xU,sqrt(xU)),dline_seg(x,xL,sqrt(xL),xU,sqrt(xU),one(T)/(2*sqrt(x)))
end
function sqrt(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = sqrt_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = sqrt_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = sqrt_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = sqrt_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = sqrt_cc(x.cc,x.Intv.lo,x.Intv.hi)
    gcv2::T,gdcv2::T = sqrt_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, sqrt(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function log_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return log(x),one(T)/x
end
function log_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,log(xL),xU,log(xU)),dline_seg(x,xL,log(xL),xU,log(xU),one(T)/x)
end
function log(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = log_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = log_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1,gdcc1 = log_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1,gdcv1 = log_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2,gdcc2 = log_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2,gdcv2 = log_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
    cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, log(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function log2_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return log2(x),one(T)/(log(2)*x)
end
function log2_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,log2(xL),xU,log2(xU)),dline_seg(x,xL,log2(xL),xU,log2(xU),one(T)/(log(2)*x))
end
function log2(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = log2_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = log2_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = log2_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = log2_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = log2_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = log2_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, log2(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function log10_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return log10(x),one(T)/(log(10)*x)
end
function log10_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,log10(xL),xU,log10(xU)),dline_seg(x,xL,log10(xL),xU,log10(xU),one(T)/(log(10)*x))
end
function log10(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = log10_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = log10_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = log10_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = log10_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = log10_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = log10_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, log10(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function acosh_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return acosh(x),one(T)/sqrt(x^2 - one(T))
end
function acosh_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return line_seg(x,xL,acosh(xL),xU,acosh(xU)),dline_seg(x,xL,acosh(xL),xU,acosh(xU),one(T)/sqrt(x^2 - one(T)))
end
function acosh(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_max::T = x.Intv.hi
  eps_min::T = x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = acosh_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = acosh_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = acosh_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = acosh_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = acosh_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = acosh_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, acosh(x.Intv),x.cnst,x.IntvBox, x.xref)
end

function abs_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (x>=zero(x))
    return xU*(x/xU)^(MC_param.mu+1),(MC_param.mu+1)*(x/xU)^MC_param.mu
  else
    return -xL*(x/xL)^(MC_param.mu+1), -(MC_param.mu+1)*(x/xL)^MC_param.mu
  end
end
function abs_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  if (x>=zero(x))
    d = (MC_param.mu+1)*(x/xU)^MC_param.mu
  else
    d = -(MC_param.mu+1)*(x/xL)^MC_param.mu
  end
  return line_seg(x,xL,abs(xL),xU,abs(xU)),dline_seg(x,xL,abs(xL),xU,abs(xU),d)
end
function abs_cc_NS(x::T,lo::T,hi::T) where {T<:AbstractFloat}
  return line_seg(x,lo,abs(lo),hi,abs(hi)),dline_seg(x,lo,abs(lo),hi,abs(hi),sign(x))
end
function abs_cv_NS(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  return abs(x),sign(x)
end
function abs(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  #println("abs x in: ", x)
  eps_min::T,blank = mid3(x.Intv.lo,x.Intv.hi,zero(x.Intv.lo))
  eps_max::T = (abs(x.Intv.hi)>=abs(x.Intv.lo)) ? x.Intv.hi : x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  if (MC_param.mu >= 1)
    cc::T,dcc::T = abs_cc(midcc,x.Intv.lo,x.Intv.hi)
    cv::T,dcv::T = abs_cv(midcv,x.Intv.lo,x.Intv.hi)
    gcc1::T,gdcc1::T = abs_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = abs_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = abs_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = abs_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc,dcc = abs_cc_NS(midcc,x.Intv.lo,x.Intv.hi)
    cv,dcv = abs_cv_NS(midcv,x.Intv.lo,x.Intv.hi)
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, abs(x.Intv),x.cnst ,x.IntvBox, x.xref)
end

function cosh_cv(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  cosh(x),-sinh(x)
end
function cosh_cc(x::T,xL::T,xU::T) where {T<:AbstractFloat}
  line_seg(x,xL,cosh(xL),xU,cosh(xU)),dline_seg(x,xL,cosh(xL),xU,cosh(xU),-sinh(x))
end
function cosh(x::SMCg{N,T}) where {N,T<:AbstractFloat}
  eps_min::T,blank = mid3(x.Intv.lo,x.Intv.hi,zero(x.Intv.lo))
  eps_max::T = (abs(x.Intv.hi)>=abs(x.Intv.lo)) ? x.Intv.hi : x.Intv.lo
  midcc::T,cc_id::Int64 = mid3(x.cc,x.cv,eps_max)
  midcv::T,cv_id::Int64 = mid3(x.cc,x.cv,eps_min)
  cc::T,dcc::T = cosh_cc(midcc,x.Intv.lo,x.Intv.hi)
  cv::T,dcv::T = cosh_cv(midcv,x.Intv.lo,x.Intv.hi)
  if (MC_param.mu >= 1)
    gcc1::T,gdcc1::T = cosh_cc(x.cv,x.Intv.lo,x.Intv.hi)
	  gcv1::T,gdcv1::T = cosh_cv(x.cv,x.Intv.lo,x.Intv.hi)
	  gcc2::T,gdcc2::T = cosh_cc(x.cc,x.Intv.lo,x.Intv.hi)
	  gcv2::T,gdcv2::T = cosh_cv(x.cc,x.Intv.lo,x.Intv.hi)
	  cv_grad::SVector{N,T} = max(zero(T),gdcv1)*x.cv_grad_grad + min(zero(T),gdcv2)*x.cc_grad
	  cc_grad::SVector{N,T} = min(zero(T),gdcc1)*x.cv_grad + max(zero(T),gdcc2)*x.cc_grad
  else
    cc_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*dcc
    cv_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
  end
  return SMCg{N,T}(cc, cv, cc_grad, cv_grad, cosh(x.Intv), x.cnst ,x.IntvBox, x.xref)
end


text:
mutable struct SVAEvampmeans <: SVAE
	q
	g
	zdim
	hue
	μzfromhidden
	κzfromhidden
	pseudo_inputs

	"""
	SVAEvampmeans(q, g, hdim, zdim, T) Constructor of the S-VAE where `zdim > 3` and T determines the floating point type (default Float32)
	"""
	SVAEvampmeans(q, g, hdim::Integer, zdim::Integer, num_pseudoinputs::Integer, T = Float32) = new(q, g, zdim, convert(T, huentropy(zdim)), Adapt.adapt(T, Chain(Dense(hdim, zdim), x -> normalizecolumns(x))), Adapt.adapt(T, Dense(hdim, 1, x -> σ.(x) .* 100)), Flux.param(Adapt.adapt(T, normalizecolumns(randn(size(q[1].W, 2), num_pseudoinputs)))))
end

function SVAEvampmeans(inputDim, hiddenDim, latentDim, numLayers, nonlinearity, layerType, num_pseudoinputs, T = Float32)
	encoder = Adapt.adapt(T, FluxExtensions.layerbuilder(inputDim, hiddenDim, hiddenDim, numLayers, nonlinearity, "", layerType))
    decoder = Adapt.adapt(T, FluxExtensions.layerbuilder(latentDim, hiddenDim, inputDim, numLayers + 1, nonlinearity, "linear", layerType))
	return SVAEvampmeans(encoder, decoder, hiddenDim, latentDim, num_pseudoinputs, T)
end

Flux.@treelike(SVAEvampmeans)

function sampleVamp(m::SVAEvampmeans, n::Integer)
	components = rand(1:size(m.pseudo_inputs, 2), n)
	(μz, κz) = zparams(m, m.pseudo_inputs[:, components])
	return samplez(m, μz, κz)
end

function wloss(m::SVAEvampmeans, x, β, d)
	(z, _) = zparams(m, x)
	zp = sampleVamp(m, size(z, 2))
	#prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, zp)
	xgivenz = m.g(z)
	return Flux.mse(x, xgivenz) + β * Ω
end

function decomposed_wloss(m::SVAEvampmeans, x, β, d)
	(z, _) = zparams(m, x)
	zp = sampleVamp(m, size(z, 2))
	#prior = samplez(m, ones(size(μz)) .* normalizecolumns(m.priorμ), ones(size(κz)) .* m.priorκ)
	Ω = d(z, zp)
	xgivenz = m.g(z)
	Flux.mse(x, xgivenz), Ω
end

log_pz(m::SVAEvampmeans, x) = log.(pz(m, x))
log_pz_from_z(m::SVAEvampmeans, z) = log.(pz_from_z(m, z))

function pz_from_z(m::SVAEvampmeans, z)
	if size(z, 2) > 1
		zs = [z[:, i] for i in 1:size(z, 2)]
		return map(x -> pz_singleinstance(m, x), zs)
	else
		return pz_singleinstance(m, z)
	end
end

pz(m::SVAEvampmeans, x) = pz_from_z(m, zparams(m, x)[1].data)

function pz_singleinstance(m::SVAEvampmeans, z)
	@assert size(z, 2) == 1
	pseudoz = zparams(m, m.pseudo_inputs)[1].data
	p = 0
	for i in 1:size(m.pseudo_inputs, 2)
		p += exp.(log_vmf_wo_c(z, pseudoz[:, i], 1))
	end
	p /= size(m.pseudo_inputs, 2)
end

function log_px(m::SVAEvampmeans, x::Matrix, k::Int = 100)
	x = [x[:, i] for i in 1:size(x, 2)]
	return map(a -> log_px(m, a, k), x)
end

function log_px(m::SVAEvampmeans, x::Vector, k::Int = 100)
	μz, κz = zparams(m, x)
	μz = repeat(Flux.Tracker.data(μz), 1, k)
	κz = repeat(Flux.Tracker.data(κz), 1, k)
	z = Flux.Tracker.data(samplez(m, μz, κz))
	xgivenz = Flux.Tracker.data(m.g(z))

	log_pxgivenz = log_normal(repeat(x, 1, k), xgivenz)
	log_pz = log_pz_from_z(m, z)
	log_qzgivenx = log_vmf_wo_c(z, μz[:, 1], κz[1])

	return log(sum(exp.(Flux.Tracker.data(log_pxgivenz .+ log_pz .- log_qzgivenx))))
end

max_stars_repo_name: honzabim/FewShotAnomalyDetection

max_issues_repo_name: honzabim/FewShotAnomalyDetection

max_forks_repo_name: honzabim/FewShotAnomalyDetection


text:
#=
Algorithms:
- Julia version: 1.4.2
- Author: akhil
- Date: 2020-05-28
=#
module Algorithms
include("Maths/Maths.jl") # Maths Directory
include("Algebra/Algebra.jl") # Algebra Directory
include("ArithematicAnalysis/ArithematicAnalysis.jl") #ArithematicAnalysis

export Maths
export Algebra
export ArithematicAnalysis
end


text:
using DataStructures: nlargest

@doc raw"""
    GlobalPool(aggr)

Global pooling layer for graph neural networks.
Takes a graph and feature nodes as inputs
and performs the operation

```math
\mathbf{u}_V = \square_{i \in V} \mathbf{x}_i
```

where ``V`` is the set of nodes of the input graph and
the type of aggregation represented by ``\square`` is selected by the `aggr` argument.
Commonly used aggregations are `mean`, `max`, and `+`.

See also [`reduce_nodes`](@ref).

# Examples

```julia
using Flux, GraphNeuralNetworks, Graphs

pool = GlobalPool(mean)

g = GNNGraph(erdos_renyi(10, 4))
X = rand(32, 10)
pool(g, X) # => 32x1 matrix

g = Flux.batch([GNNGraph(erdos_renyi(10, 4)) for _ in 1:5])
X = rand(32, 50)
pool(g, X) # => 32x5 matrix
```
"""
struct GlobalPool{F} <: GNNLayer
    aggr::F
end

function (l::GlobalPool)(g::GNNGraph, x::AbstractArray)
    return reduce_nodes(l.aggr, g, x)
end

(l::GlobalPool)(g::GNNGraph) = GNNGraph(g, gdata=l(g, node_features(g)))

@doc raw"""
    GlobalAttentionPool(fgate, ffeat=identity)

Global soft attention layer from the [Gated Graph Sequence Neural
Networks](https://arxiv.org/abs/1511.05493) paper

```math
\mathbf{u}_V = \sum_{i\in V} \alpha_i\, f_{feat}(\mathbf{x}_i)
```

where the coefficients ``\alpha_i`` are given by a [`softmax_nodes`](@ref)
operation:

```math
\alpha_i = \frac{e^{f_{gate}(\mathbf{x}_i)}}
                {\sum_{i'\in V} e^{f_{gate}(\mathbf{x}_{i'})}}.
```

# Arguments

- `fgate`: The function ``f_{gate}: \mathbb{R}^{D_{in}} \to \mathbb{R}``.
           It is tipically expressed by a neural network.

- `ffeat`: The function ``f_{feat}: \mathbb{R}^{D_{in}} \to \mathbb{R}^{D_{out}}``.
           It is tipically expressed by a neural network.

# Examples

```julia
chin = 6
chout = 5

fgate = Dense(chin, 1)
ffeat = Dense(chin, chout)
pool = GlobalAttentionPool(fgate, ffeat)

g = Flux.batch([GNNGraph(random_regular_graph(10, 4),
                         ndata=rand(Float32, chin, 10))
                for i=1:3])

u = pool(g, g.ndata.x)

@assert size(u) == (chout, g.num_graphs)
```
"""
struct GlobalAttentionPool{G,F}
    fgate::G
    ffeat::F
end

@functor GlobalAttentionPool

GlobalAttentionPool(fgate) = GlobalAttentionPool(fgate, identity)

function (l::GlobalAttentionPool)(g::GNNGraph, x::AbstractArray)
    α = softmax_nodes(g, l.fgate(x))
    feats = α .* l.ffeat(x)
    u = reduce_nodes(+, g, feats)
    return u
end

(l::GlobalAttentionPool)(g::GNNGraph) = GNNGraph(g, gdata=l(g, node_features(g)))

"""
    TopKPool(adj, k, in_channel)

Top-k pooling layer.

# Arguments

- `adj`: Adjacency matrix  of a graph.
- `k`: Top-k nodes are selected to pool together.
- `in_channel`: The dimension of input channel.
"""
struct TopKPool{T,S}
    A::AbstractMatrix{T}
    k::Int
    p::AbstractVector{S}
    Ã::AbstractMatrix{T}
end

function TopKPool(adj::AbstractMatrix, k::Int, in_channel::Int; init=glorot_uniform)
    TopKPool(adj, k, init(in_channel), similar(adj, k, k))
end

function (t::TopKPool)(X::AbstractArray)
    y = t.p' * X / norm(t.p)
    idx = topk_index(y, t.k)
    t.Ã .= view(t.A, idx, idx)
    X_ = view(X, :, idx) .* σ.(view(y, idx)')
    return X_
end

function topk_index(y::AbstractVector, k::Int)
    v = nlargest(k, y)
    return collect(1:length(y))[y .>= v[end]]
end

topk_index(y::Adjoint, k::Int) = topk_index(y', k)


text:
# import MISPTools
import PyCall

# PyCall.PyObject(::typeof(MISPTools.solve)) =
#     PyCall.pyfunctionret(MISPTools.solve,Any,Vararg{PyCall.PyAny})

# function MISPTools.numargs(f::PyCall.PyObject)
#     inspect = PyCall.pyimport("inspect")
#     haskey(f,:py_func) ? _f = f[:py_func] : _f = f
#     if PyCall.pyversion < v"3.0.0"
#         return length(first(inspect[:getargspec](_f)))
#     else
#         return length(first(inspect[:getfullargspec](_f)))
#     end
# end

max_stars_repo_name: jgoldfar/misppy

max_issues_repo_name: jgoldfar/misppy

max_forks_repo_name: jgoldfar/misppy


text:
#!/usr/bin/env julia

using Pkg

cfgribjl_dir = normpath(joinpath(@__DIR__, "..", ".."))

Pkg.activate(cfgribjl_dir); print("\n")

println(
    """
          Running pre-commit hook
    -------------------------------------
    | Testing code imports to catch any |
    |   syntax or dependency problems   |
    -------------------------------------
    """
)

include(joinpath(cfgribjl_dir, "./src/cfgrib.jl"))


text:
# convenience and indirect construction

function Base.fill!(A::AnyGPUArray{T}, x) where T
    length(A) == 0 && return A
    gpu_call(A, convert(T, x)) do ctx, a, val
        idx = @linearidx(a)
        @inbounds a[idx] = val
        return
    end
    A
end

## identity matrices

function identity_kernel(ctx::AbstractKernelContext, res::AbstractArray{T}, stride, val) where T
    i = linear_index(ctx)
    i > stride && return
    ilin = (stride * (i - 1)) + i
    @inbounds res[ilin] = val
    return
end

function (T::Type{<: AnyGPUArray{U}})(s::UniformScaling, dims::Dims{2}) where {U}
    res = similar(T, dims)
    fill!(res, zero(U))
    gpu_call(identity_kernel, res, size(res, 1), s.λ; total_threads=minimum(dims))
    res
end

(T::Type{<: AnyGPUArray})(s::UniformScaling{U}, dims::Dims{2}) where U = T{U}(s, dims)

(T::Type{<: AnyGPUArray})(s::UniformScaling, m::Integer, n::Integer) = T(s, Dims((m, n)))

function Base.copyto!(A::AbstractGPUMatrix{T}, s::UniformScaling) where T
    fill!(A, zero(T))
    gpu_call(identity_kernel, A, size(A, 1), s.λ; total_threads=minimum(size(A)))
    A
end

function _one(unit::T, x::AbstractGPUMatrix) where {T}
    m,n = size(x)
    m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    I = similar(x, T)
    fill!(I, zero(T))
    gpu_call(identity_kernel, I, m, unit; total_threads=m)
    I
end

Base.one(x::AbstractGPUMatrix{T}) where {T} = _one(one(T), x)
Base.oneunit(x::AbstractGPUMatrix{T}) where {T} = _one(oneunit(T), x)

## collect & convert

function indexstyle(x::T) where T
    style = try
        Base.IndexStyle(x)
    catch
        nothing
    end
    style
end

function collect_kernel(ctx::AbstractKernelContext, A, iter, ::IndexCartesian)
    idx = @cartesianidx(A)
    @inbounds A[idx...] = iter[idx...]
    return
end

function collect_kernel(ctx::AbstractKernelContext, A, iter, ::IndexLinear)
    idx = linear_index(ctx)
    @inbounds A[idx] = iter[idx]
    return
end

eltype_or(::Type{<: AbstractGPUArray}, or) = or
eltype_or(::Type{<: AbstractGPUArray{T}}, or) where T = T
eltype_or(::Type{<: AbstractGPUArray{T, N}}, or) where {T, N} = T

function Base.convert(AT::Type{<: AbstractGPUArray}, iter)
    isize = Base.IteratorSize(iter)
    style = indexstyle(iter)
    ettrait = Base.IteratorEltype(iter)
    if isbits(iter) && isa(isize, Base.HasShape) && style != nothing && isa(ettrait, Base.HasEltype)
        # We can collect on the GPU
        A = similar(AT, eltype_or(AT, eltype(iter)), size(iter))
        gpu_call(collect_kernel, A, iter, style)
        A
    else
        convert(AT, collect(iter))
    end
end

function Base.convert(AT::Type{<: AbstractGPUArray{T, N}}, A::DenseArray{T, N}) where {T, N}
    copyto!(AT(undef, size(A)), A)
end

function Base.convert(AT::Type{<: AbstractGPUArray{T1}}, A::DenseArray{T2, N}) where {T1, T2, N}
    copyto!(similar(AT, size(A)), convert(Array{T1, N}, A))
end

function Base.convert(AT::Type{<: AbstractGPUArray}, A::DenseArray{T2, N}) where {T2, N}
    copyto!(similar(AT{T2}, size(A)), A)
end

function Base.convert(AT::Type{Array{T, N}}, A::AbstractGPUArray{CT, CN}) where {T, N, CT, CN}
    convert(AT, copyto!(Array{CT, CN}(undef, size(A)), A))
end


text:
using JuliaPackageDemo
using MittagLeffler
using Test

const testx = 0.5

function _coolfunction1(x)
    (2*x + x^2) == coolfunction1(x)
end
@test _coolfunction1(testx)

function _coolfunction2(x)
    (50*mittleff(0.5, -x) - 10) == coolfunction2(x)
end
@test _coolfunction2(testx)

include(joinpath(@__DIR__, "..", "docs", "src", "simpleexamples.jl"))

max_stars_repo_name: moustachio-belvedere/JuliaPackageDemo

max_issues_repo_name: moustachio-belvedere/JuliaPackageDemo

max_forks_repo_name: moustachio-belvedere/JuliaPackageDemo


text:
abstract type Transform end
#TODO:code is repeated for error
#TODO: create reusable func
function encode(item)
    if item===nothing
        error("A Transform must be provided")
    end
    item
end

function decode(item)
    if item===nothing
        error("A Transform must be provided to decode")
    end
    item
end


text:
@system Sheath(Organ) begin
    rank ~ ::int(override) # preserve

    #FIXME sheath biomass
    mass => 0 ~ track(u"g")
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@doc raw"""
    ShallowWaterEquations2D(gravity, H0)

Shallow water equations (SWE) in two space dimensions. The equations are given by
```math
\begin{aligned}
  \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h v_1)
    + \frac{\partial}{\partial y}(h v_2) &= 0 \\
    \frac{\partial}{\partial t}(h v_1) + \frac{\partial}{\partial x}\left(h v_1^2 + \frac{g}{2}h^2\right)
    + \frac{\partial}{\partial y}(h v_1 v_2) + g h \frac{\partial b}{\partial x} &= 0 \\
    \frac{\partial}{\partial t}(h v_2) + \frac{\partial}{\partial x}(h v_1 v_2)
    + \frac{\partial}{\partial y}\left(h v_2^2 + \frac{g}{2}h^2\right) + g h \frac{\partial b}{\partial y} &= 0.
\end{aligned}
```
The unknown quantities of the SWE are the water height ``h`` and the velocities ``\mathbf{v} = (v_1, v_2)^T``.
The gravitational constant is denoted by `g` and the (possibly) variable bottom topography function ``b(x,y)``.
Conservative variable water height ``h`` is measured from the bottom topography ``b``, therefore one
also defines the total water height as ``H = h + b``.

The additional quantity ``H_0`` is also available to store a reference value for the total water height that
is useful to set initial conditions or test the "lake-at-rest" well-balancedness.

The bottom topography function ``b(x,y)`` is set inside the initial condition routine
for a particular problem setup. To test the conservative form of the SWE one can set the bottom topography
variable `b` to zero.

In addition to the unknowns, Trixi currently stores the bottom topography values at the approximation points
despite being fixed in time. This is done for convenience of computing the bottom topography gradients
on the fly during the approximation as well as computing auxiliary quantities like the total water height ``H``
or the entropy variables.
This affects the implementation and use of these equations in various ways:
* The flux values corresponding to the bottom topography must be zero.
* The bottom topography values must be included when defining initial conditions, boundary conditions or
  source terms.
* [`AnalysisCallback`](@ref) analyzes this variable.
* Trixi's visualization tools will visualize the bottom topography by default.

References for the SWE are many but a good introduction is available in Chapter 13 of the book:
- Randall J. LeVeque (2002)
  Finite Volume Methods for Hyperbolic Problems
  [DOI: 10.1017/CBO9780511791253](https://doi.org/10.1017/CBO9780511791253)
"""
struct ShallowWaterEquations2D{RealT<:Real} <: AbstractShallowWaterEquations{2, 4}
  gravity::RealT # gravitational constant
  H0::RealT      # constant "lake-at-rest" total water height
end

# Allow for flexibility to set the gravitational constant within an elixir depending on the
# application where `gravity_constant=1.0` or `gravity_constant=9.81` are common values.
# The reference total water height H0 defaults to 0.0 but is used for the "lake-at-rest"
# well-balancedness test cases
function ShallowWaterEquations2D(; gravity_constant, H0=0.0)
  ShallowWaterEquations2D(gravity_constant, H0)
end

have_nonconservative_terms(::ShallowWaterEquations2D) = Val(true)
varnames(::typeof(cons2cons), ::ShallowWaterEquations2D) = ("h", "h_v1", "h_v2", "b")
# Note, we use the total water height, H = h + b, as the first primitive variable for easier
# visualization and setting initial conditions
varnames(::typeof(cons2prim), ::ShallowWaterEquations2D) = ("H", "v1", "v2", "b")

# Set initial conditions at physical location `x` for time `t`
"""
    initial_condition_convergence_test(x, t, equations::ShallowWaterEquations2D)

A smooth initial condition used for convergence tests in combination with
[`source_terms_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).
"""
function initial_condition_convergence_test(x, t, equations::ShallowWaterEquations2D)
  # some constants are chosen such that the function is periodic on the domain [0,sqrt(2)]^2
  c  = 7.0
  omega_x = 2.0 * pi * sqrt(2.0)
  omega_t = 2.0 * pi

  x1, x2 = x

  H = c + cos(omega_x * x1) * sin(omega_x * x2) * cos(omega_t * t)
  v1 = 0.5
  v2 = 1.5
  b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2)
  return prim2cons(SVector(H, v1, v2, b), equations)
end

"""
    source_terms_convergence_test(u, x, t, equations::ShallowWaterEquations2D)

Source terms used for convergence tests in combination with
[`initial_condition_convergence_test`](@ref)
(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains).

This manufactured solution source term is specifically designed for the bottom topography function
`b(x,y) = 2 + 0.5 * sin(sqrt(2)*pi*x) + 0.5 * sin(sqrt(2)*pi*y)`
as defined in [`initial_condition_convergence_test`](@ref).
"""
@inline function source_terms_convergence_test(u, x, t, equations::ShallowWaterEquations2D)
  # Same settings as in `initial_condition_convergence_test`. Some derivative simplify because
  # this manufactured solution velocities are taken to be constants
  c  = 7.0
  omega_x = 2.0 * pi * sqrt(2.0)
  omega_t = 2.0 * pi
  omega_b = sqrt(2.0) * pi
  v1 = 0.5
  v2 = 1.5

  x1, x2 = x

  sinX, cosX = sincos(omega_x * x1)
  sinY, cosY = sincos(omega_x * x2)
  sinT, cosT = sincos(omega_t * t )

  H = c + cosX * sinY * cosT
  H_x = -omega_x * sinX * sinY * cosT
  H_y =  omega_x * cosX * cosY * cosT
  # this time derivative for the water height exploits that the bottom topography is
  # fixed in time such that H_t = (h+b)_t = h_t + 0
  H_t = -omega_t * cosX * sinY * sinT

  # bottom topography and its gradient
  b = 2.0 + 0.5 * sin(sqrt(2.0) * pi * x1) + 0.5 * sin(sqrt(2.0) * pi * x2)
  tmp1 = 0.5 * omega_b
  b_x = tmp1 * cos(omega_b * x1)
  b_y = tmp1 * cos(omega_b * x2)

  du1 = H_t + v1 * (H_x - b_x) + v2 * (H_y - b_y)
  du2 = v1 * du1 + equations.gravity * (H - b) * H_x
  du3 = v2 * du1 + equations.gravity * (H - b) * H_y
  return SVector(du1, du2, du3, 0.0)
end

"""
    initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquations2D)

A weak blast wave discontinuity useful for testing, e.g., total energy conservation.
Note for the shallow water equations to the total energy acts as a mathematical entropy function.
"""
function initial_condition_weak_blast_wave(x, t, equations::ShallowWaterEquations2D)
  # Set up polar coordinates
  inicenter = SVector(0.7, 0.7)
  x_norm = x[1] - inicenter[1]
  y_norm = x[2] - inicenter[2]
  r = sqrt(x_norm^2 + y_norm^2)
  phi = atan(y_norm, x_norm)
  sin_phi, cos_phi = sincos(phi)

  # Calculate primitive variables
  H = r > 0.5 ? 3.25 : 4.0
  v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
  v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
  b = 0.0 # by default assume there is no bottom topography

  return prim2cons(SVector(H, v1, v2, b), equations)
end

"""
    boundary_condition_slip_wall(u_inner, normal_direction, x, t, surface_flux_function,
                                 equations::ShallowWaterEquations2D)

Create a boundary state by reflecting the normal velocity component and keep
the tangential velocity component unchanged. The boundary water height is taken from
the internal value.

For details see Section 9.2.5 of the book:
- Eleuterio F. Toro (2001)
  Shock-Capturing Methods for Free-Surface Shallow Flows
  1st edition
  ISBN 0471987662
"""
@inline function boundary_condition_slip_wall(u_inner, normal_direction::AbstractVector,
                                              x, t,
                                              surface_flux_function,
                                              equations::ShallowWaterEquations2D)
  # normalize the outward pointing direction
  normal = normal_direction / norm(normal_direction)

  # compute the normal velocity
  u_normal = normal[1] * u_inner[2] + normal[2] * u_inner[3]

  # create the "external" boundary solution state
  u_boundary = SVector(u_inner[1],
                       u_inner[2] - 2.0 * u_normal * normal[1],
                       u_inner[3] - 2.0 * u_normal * normal[2],
                       u_inner[4])

  # calculate the boundary flux
  flux = surface_flux_function(u_inner, u_boundary, normal_direction, equations)

  return flux
end

# Calculate 1D flux for a single point
# Note, the bottom topography has no flux
@inline function flux(u, orientation::Integer, equations::ShallowWaterEquations2D)
  h, h_v1, h_v2, _ = u
  v1, v2 = velocity(u, equations)

  p = 0.5 * equations.gravity * h^2
  if orientation == 1
    f1 = h_v1
    f2 = h_v1 * v1 + p
    f3 = h_v1 * v2
  else
    f1 = h_v2
    f2 = h_v2 * v1
    f3 = h_v2 * v2 + p
  end
  return SVector(f1, f2, f3, zero(eltype(u)))
end

# Calculate 1D flux for a single point in the normal direction
# Note, this directional vector is not normalized and the bottom topography has no flux
@inline function flux(u, normal_direction::AbstractVector, equations::ShallowWaterEquations2D)
  h = waterheight(u, equations)
  v1, v2 = velocity(u, equations)

  v_normal = v1 * normal_direction[1] + v2 * normal_direction[2]
  h_v_normal = h * v_normal
  p = 0.5 * equations.gravity * h^2

  f1 = h_v_normal
  f2 = h_v_normal * v1 + p * normal_direction[1]
  f3 = h_v_normal * v2 + p * normal_direction[2]
  return SVector(f1, f2, f3, zero(eltype(u)))
end

"""
    flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                          equations::ShallowWaterEquations2D)
    flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                          normal_direction_ll     ::AbstractVector,
                                          normal_direction_average::AbstractVector,
                                          equations::ShallowWaterEquations2D)

Non-symmetric two-point volume flux discretizing the nonconservative (source) term
that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).

On curvilinear meshes, this nonconservative flux depends on both the
contravariant vector (normal direction) at the current node and the averaged
one. This is different from numerical fluxes used to discretize conservative
terms.

Further details are available in the paper:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
  An entropy stable nodal discontinuous Galerkin method for the two dimensional
  shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
  [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
"""
@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr, orientation::Integer,
                                                       equations::ShallowWaterEquations2D)
  # Pull the necessary left and right state information
  h_ll = waterheight(u_ll, equations)
  b_rr = u_rr[4]

  z = zero(eltype(u_ll))
  # Bottom gradient nonconservative term: (0, g h b_x, g h b_y, 0)
  if orientation == 1
    f = SVector(z, equations.gravity * h_ll * b_rr, z, z)
  else # orientation == 2
    f = SVector(z, z, equations.gravity * h_ll * b_rr, z)
  end
  return f
end

@inline function flux_nonconservative_wintermeyer_etal(u_ll, u_rr,
                                                       normal_direction_ll::AbstractVector,
                                                       normal_direction_average::AbstractVector,
                                                       equations::ShallowWaterEquations2D)
  # Pull the necessary left and right state information
  h_ll = waterheight(u_ll, equations)
  b_rr = u_rr[4]
  # Note this routine only uses the `normal_direction_average` and the average of the
  # bottom topography to get a quadratic split form DG gradient on curved elements
  return SVector(zero(eltype(u_ll)),
                 normal_direction_average[1] * equations.gravity * h_ll * b_rr,
                 normal_direction_average[2] * equations.gravity * h_ll * b_rr,
                 zero(eltype(u_ll)))
end

"""
    flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                        equations::ShallowWaterEquations2D)
    flux_nonconservative_fjordholm_etal(u_ll, u_rr,
                                        normal_direction_ll     ::AbstractVector,
                                        normal_direction_average::AbstractVector,
                                        equations::ShallowWaterEquations2D)

Non-symmetric two-point surface flux discretizing the nonconservative (source) term of
that contains the gradient of the bottom topography [`ShallowWaterEquations2D`](@ref).

On curvilinear meshes, this nonconservative flux depends on both the
contravariant vector (normal direction) at the current node and the averaged
one. This is different from numerical fluxes used to discretize conservative
terms.

This contains additional terms compared to [`flux_nonconservative_wintermeyer_etal`](@ref)
that account for possible discontinuities in the bottom topography function.
Thus, this flux should be used in general at interfaces. For flux differencing volume terms,
[`flux_nonconservative_wintermeyer_etal`](@ref) is analytically equivalent but slightly
cheaper.

Further details for the original finite volume formulation are available in
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
  Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
  [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
and for curvilinear 2D case in the paper:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
  An entropy stable nodal discontinuous Galerkin method for the two dimensional
  shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
  [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
"""
@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr, orientation::Integer,
                                                     equations::ShallowWaterEquations2D)
  # Pull the necessary left and right state information
  h_ll, _, _, b_ll = u_ll
  h_rr, _, _, b_rr = u_rr

  h_average = 0.5 * (h_ll + h_rr)
  b_jump = b_rr - b_ll

  # Includes two parts:
  #   (i)  Diagonal (consistent) term from the volume flux that uses `b_ll` to avoid
  #        cross-averaging across a discontinuous bottom topography
  #   (ii) True surface part that uses `h_average` and `b_jump` to handle discontinuous bathymetry
  z = zero(eltype(u_ll))
  if orientation == 1
    f = SVector(z,
                equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump,
                z, z)
  else # orientation == 2
    f = SVector(z, z,
                equations.gravity * h_ll * b_ll + equations.gravity * h_average * b_jump,
                z)
  end

  return f
end

@inline function flux_nonconservative_fjordholm_etal(u_ll, u_rr,
                                                     normal_direction_ll::AbstractVector,
                                                     normal_direction_average::AbstractVector,
                                                     equations::ShallowWaterEquations2D)
  # Pull the necessary left and right state information
  h_ll, _, _, b_ll = u_ll
  h_rr, _, _, b_rr = u_rr

  # Comes in two parts:
  #   (i)  Diagonal (consistent) term from the volume flux that uses `normal_direction_average`
  #        but we use `b_ll` to avoid cross-averaging across a discontinuous bottom topography
  f2 = normal_direction_average[1] * equations.gravity * h_ll * b_ll
  f3 = normal_direction_average[2] * equations.gravity * h_ll * b_ll

  #   (ii) True surface part that uses `normal_direction_ll`, `h_average` and `b_jump`
  #        to handle discontinuous bathymetry
  h_average = 0.5 * (h_ll + h_rr)
  b_jump = b_rr - b_ll

  f2 += normal_direction_ll[1] * equations.gravity * h_average * b_jump
  f3 += normal_direction_ll[2] * equations.gravity * h_average * b_jump

  # First and last equations do not have a nonconservative flux
  f1 = f4 = zero(eltype(u_ll))

  return SVector(f1, f2, f3, f4)
end

"""
    flux_fjordholm_etal(u_ll, u_rr, orientation_or_normal_direction,
                        equations::ShallowWaterEquations2D)

Total energy conservative (mathematical entropy for shallow water equations). When the bottom topography
is nonzero this should only be used as a surface flux otherwise the scheme will not be well-balanced.
For well-balancedness in the volume flux use [`flux_wintermeyer_etal`](@ref).

Details are available in Eq. (4.1) in the paper:
- Ulrik S. Fjordholm, Siddhartha Mishr and Eitan Tadmor (2011)
  Well-balanced and energy stable schemes for the shallow water equations with discontinuous topography
  [DOI: 10.1016/j.jcp.2011.03.042](https://doi.org/10.1016/j.jcp.2011.03.042)
"""
@inline function flux_fjordholm_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D)
  # Unpack left and right state
  h_ll = waterheight(u_ll, equations)
  v1_ll, v2_ll = velocity(u_ll, equations)
  h_rr = waterheight(u_rr, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)

  # Average each factor of products in flux
  h_avg  = 0.5 * (h_ll   + h_rr  )
  v1_avg = 0.5 * (v1_ll  + v1_rr )
  v2_avg = 0.5 * (v2_ll  + v2_rr )
  p_avg  = 0.25 * equations.gravity * (h_ll^2 + h_rr^2)

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = h_avg * v1_avg
    f2 = f1 * v1_avg + p_avg
    f3 = f1 * v2_avg
  else
    f1 = h_avg * v2_avg
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_avg
  end

  return SVector(f1, f2, f3, zero(eltype(u_ll)))
end

@inline function flux_fjordholm_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::ShallowWaterEquations2D)
  # Unpack left and right state
  h_ll = waterheight(u_ll, equations)
  v1_ll, v2_ll = velocity(u_ll, equations)
  h_rr = waterheight(u_rr, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)

  v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
  v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

  # Average each factor of products in flux
  h_avg  = 0.5 * (h_ll   + h_rr  )
  v1_avg = 0.5 * (v1_ll  + v1_rr )
  v2_avg = 0.5 * (v2_ll  + v2_rr )
  h2_avg = 0.5 * (h_ll^2 + h_rr^2)
  p_avg  = 0.5 * equations.gravity * h2_avg
  v_dot_n_avg = 0.5 * (v_dot_n_ll + v_dot_n_rr)

  # Calculate fluxes depending on normal_direction
  f1 = h_avg * v_dot_n_avg
  f2 = f1 * v1_avg + p_avg * normal_direction[1]
  f3 = f1 * v2_avg + p_avg * normal_direction[2]

  return SVector(f1, f2, f3, zero(eltype(u_ll)))
end

"""
    flux_wintermeyer_etal(u_ll, u_rr, orientation_or_normal_direction,
                          equations::ShallowWaterEquations2D)

Total energy conservative (mathematical entropy for shallow water equations) split form.
When the bottom topography is nonzero this scheme will be well-balanced when used as a `volume_flux`.
The `surface_flux` should still use, e.g., [`flux_fjordholm_etal`](@ref).

Further details are available in Theorem 1 of the paper:
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and David A. Kopriva (2017)
  An entropy stable nodal discontinuous Galerkin method for the two dimensional
  shallow water equations on unstructured curvilinear meshes with discontinuous bathymetry
  [DOI: 10.1016/j.jcp.2017.03.036](https://doi.org/10.1016/j.jcp.2017.03.036)
"""
@inline function flux_wintermeyer_etal(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D)
  # Unpack left and right state
  h_ll, h_v1_ll, h_v2_ll, _ = u_ll
  h_rr, h_v1_rr, h_v2_rr, _ = u_rr

  # Get the velocities on either side
  v1_ll, v2_ll = velocity(u_ll, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)

  # Average each factor of products in flux
  v1_avg = 0.5 * (v1_ll + v1_rr)
  v2_avg = 0.5 * (v2_ll + v2_rr)
  p_avg  = 0.5 * equations.gravity * h_ll * h_rr

  # Calculate fluxes depending on orientation
  if orientation == 1
    f1 = 0.5 * (h_v1_ll + h_v1_rr)
    f2 = f1 * v1_avg + p_avg
    f3 = f1 * v2_avg
  else
    f1 = 0.5 * (h_v2_ll + h_v2_rr)
    f2 = f1 * v1_avg
    f3 = f1 * v2_avg + p_avg
  end

  return SVector(f1, f2, f3, zero(eltype(u_ll)))
end

@inline function flux_wintermeyer_etal(u_ll, u_rr, normal_direction::AbstractVector, equations::ShallowWaterEquations2D)
  # Unpack left and right state
  h_ll, h_v1_ll, h_v2_ll, _ = u_ll
  h_rr, h_v1_rr, h_v2_rr, _ = u_rr

  # Get the velocities on either side
  v1_ll, v2_ll = velocity(u_ll, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)

  # Average each factor of products in flux
  h_v1_avg = 0.5 * (h_v1_ll + h_v1_rr )
  h_v2_avg = 0.5 * (h_v2_ll + h_v2_rr )
  v1_avg   = 0.5 * (v1_ll   + v1_rr   )
  v2_avg   = 0.5 * (v2_ll   + v2_rr   )
  p_avg    = 0.5 * equations.gravity * h_ll * h_rr

  # Calculate fluxes depending on normal_direction
  f1 = h_v1_avg * normal_direction[1] + h_v2_avg * normal_direction[2]
  f2 = f1 * v1_avg + p_avg * normal_direction[1]
  f3 = f1 * v2_avg + p_avg * normal_direction[2]

  return SVector(f1, f2, f3, zero(eltype(u_ll)))
end

# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation as the
# maximum velocity magnitude plus the maximum speed of sound
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, equations::ShallowWaterEquations2D)
  # Get the velocity quantities in the appropriate direction
  if orientation == 1
    v_ll, _ = velocity(u_ll, equations)
    v_rr, _ = velocity(u_rr, equations)
  else
    _, v_ll = velocity(u_ll, equations)
    _, v_rr = velocity(u_rr, equations)
  end

  # Calculate the wave celerity on the left and right
  h_ll = waterheight(u_ll, equations)
  h_rr = waterheight(u_rr, equations)
  c_ll = sqrt(equations.gravity * h_ll)
  c_rr = sqrt(equations.gravity * h_rr)

  λ_max = max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr)
end

@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector, equations::ShallowWaterEquations2D)
  # Extract and compute the velocities in the normal direction
  v1_ll, v2_ll = velocity(u_ll, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)
  v_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
  v_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

  # Compute the wave celerity on the left and right
  h_ll = waterheight(u_ll, equations)
  h_rr = waterheight(u_rr, equations)
  c_ll = sqrt(equations.gravity * h_ll)
  c_rr = sqrt(equations.gravity * h_rr)

  # The normal velocities are already scaled by the norm
  return max(abs(v_ll), abs(v_rr)) + max(c_ll, c_rr) * norm(normal_direction)
end

# Specialized `DissipationLocalLaxFriedrichs` to avoid spurious dissipation in the bottom topography
@inline function (dissipation::DissipationLocalLaxFriedrichs)(u_ll, u_rr, orientation_or_normal_direction,
                                                              equations::ShallowWaterEquations2D)
  λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)
  diss = -0.5 * λ * (u_rr - u_ll)
  return SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll)))
end

# Specialized `FluxHLL` to avoid spurious dissipation in the bottom topography
@inline function (numflux::FluxHLL)(u_ll, u_rr, orientation_or_normal_direction,
                  equations::ShallowWaterEquations2D)
  λ_min, λ_max = numflux.min_max_speed(u_ll, u_rr, orientation_or_normal_direction, equations)

  if λ_min >= 0 && λ_max >= 0
    return flux(u_ll, orientation_or_normal_direction, equations)
  elseif λ_max <= 0 && λ_min <= 0
    return flux(u_rr, orientation_or_normal_direction, equations)
  else
    f_ll = flux(u_ll, orientation_or_normal_direction, equations)
    f_rr = flux(u_rr, orientation_or_normal_direction, equations)
    inv_λ_max_minus_λ_min = inv(λ_max - λ_min)
    factor_ll = λ_max * inv_λ_max_minus_λ_min
    factor_rr = λ_min * inv_λ_max_minus_λ_min
    factor_diss = λ_min * λ_max * inv_λ_max_minus_λ_min
    diss = u_rr - u_ll
    return factor_ll * f_ll - factor_rr * f_rr + factor_diss * SVector(diss[1], diss[2], diss[3], zero(eltype(u_ll)))
  end
end

# Calculate minimum and maximum wave speeds for HLL-type fluxes
@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer,
                                     equations::ShallowWaterEquations2D)
  h_ll = waterheight(u_ll, equations)
  v1_ll, v2_ll = velocity(u_ll, equations)
  h_rr = waterheight(u_rr, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)

  if orientation == 1 # x-direction
    λ_min = v1_ll - sqrt(equations.gravity * h_ll)
    λ_max = v1_rr + sqrt(equations.gravity * h_rr)
  else # y-direction
    λ_min = v2_ll - sqrt(equations.gravity * h_ll)
    λ_max = v2_rr + sqrt(equations.gravity * h_rr)
  end

  return λ_min, λ_max
end

@inline function min_max_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
                                     equations::ShallowWaterEquations2D)
  h_ll = waterheight(u_ll, equations)
  v1_ll, v2_ll = velocity(u_ll, equations)
  h_rr = waterheight(u_rr, equations)
  v1_rr, v2_rr = velocity(u_rr, equations)

  v_normal_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
  v_normal_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]

  norm_ = norm(normal_direction)
  # The v_normals are already scaled by the norm
  λ_min = v_normal_ll - sqrt(equations.gravity * h_ll) * norm_
  λ_max = v_normal_rr + sqrt(equations.gravity * h_rr) * norm_

  return λ_min, λ_max
end

@inline function max_abs_speeds(u, equations::ShallowWaterEquations2D)
  h = waterheight(u, equations)
  v1, v2 = velocity(u, equations)

  c = equations.gravity * sqrt(h)
  return abs(v1) + c, abs(v2) + c
end

# Helper function to extract the velocity vector from the conservative variables
@inline function velocity(u, equations::ShallowWaterEquations2D)
  h, h_v1, h_v2, _ = u

  v1 = h_v1 / h
  v2 = h_v2 / h
  return SVector(v1, v2)
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::ShallowWaterEquations2D)
  h, _, _, b = u

  H = h + b
  v1, v2 = velocity(u, equations)
  return SVector(H, v1, v2, b)
end

# Convert conservative variables to entropy
# Note, only the first three are the entropy variables, the fourth entry still
# just carries the bottom topography values for convenience
@inline function cons2entropy(u, equations::ShallowWaterEquations2D)
  h, h_v1, h_v2, b = u

  v1, v2 = velocity(u, equations)
  v_square = v1^2 + v2^2

  w1 = equations.gravity * (h + b) - 0.5 * v_square
  w2 = v1
  w3 = v2
  return SVector(w1, w2, w3, b)
end

# Convert entropy variables to conservative
@inline function entropy2cons(w, equations::ShallowWaterEquations2D)
  w1, w2, w3, b = w

  h = (w1 + 0.5 * (w2^2 + w3^2)) / equations.gravity - b
  h_v1 = h * w2
  h_v2 = h * w3
  return SVector(h, h_v1, h_v2, b)
end

# Convert primitive to conservative variables
@inline function prim2cons(prim, equations::ShallowWaterEquations2D)
  H, v1, v2, b = prim

  h = H - b
  h_v1 = h * v1
  h_v2 = h * v2
  return SVector(h, h_v1, h_v2, b)
end

@inline function waterheight(u, equations::ShallowWaterEquations2D)
  return u[1]
end

@inline function pressure(u, equations::ShallowWaterEquations2D)
  h = waterheight(u, equations)
  p = 0.5 * equations.gravity * h^2
  return p
end

@inline function waterheight_pressure(u, equations::ShallowWaterEquations2D)
  return waterheight(u, equations) * pressure(u, equations)
end

# Entropy function for the shallow water equations is the total energy
@inline entropy(cons, equations::ShallowWaterEquations2D) = energy_total(cons, equations)

# Calculate total energy for a conservative state `cons`
@inline function energy_total(cons, equations::ShallowWaterEquations2D)
  h, h_v1, h_v2, b = cons

  e = (h_v1^2 + h_v2^2) / (2 * h) + 0.5 * equations.gravity * h^2 + equations.gravity * h * b
  return e
end

# Calculate kinetic energy for a conservative state `cons`
@inline function energy_kinetic(u, equations::ShallowWaterEquations2D)
  h, h_v1, h_v2, _ = u
  return (h_v1^2 + h_v2^2) / (2 * h)
end

# Calculate potential energy for a conservative state `cons`
@inline function energy_internal(cons, equations::ShallowWaterEquations2D)
  return energy_total(cons, equations) - energy_kinetic(cons, equations)
end

# Calculate the error for the "lake-at-rest" test case where H = h+b should
# be a constant value over time
@inline function lake_at_rest_error(u, equations::ShallowWaterEquations2D)
  h, _, _, b = u
  return abs(equations.H0 - (h + b))
end

end # @muladd


text:
#
# Fast API module. Contains all available data types lists splited by
# Void type. For example: type T i::Int end will be stored like:
# DataType[Int, Void]. Void type is a separator between types. So
# if you have no parameters you may set dataIndex to Void type.
#
# @author DeadbraiN
#
module FastApi
  export API_UINT64
  export API_ASCIISTRING
  export API_UINT8
  export API_UINT16
  export API_DOT_COLOR
  export API_ORG_COLOR
  export API_FLOAT64
  #
  # Indexes in _api field, which are used in request() method in
  # fast mode.
  #
  const API_UINT64        = UInt8(1)
  const API_ASCIISTRING   = UInt8(3)
  const API_UINT8         = UInt8(5)
  const API_UINT16        = UInt8(7)
  const API_DOT_COLOR     = UInt8(9)
  const API_ORG_COLOR     = UInt8(13)
  const API_FLOAT64       = UInt8(19)
  #
  # Just returns all available types API
  # @return {Array{DataType, 1}}
  #
  function getTypes()
    _api::Array{DataType, 1}
  end
  #
  # Available data types with indexes in comments
  #
  const _api = DataType[
    UInt64, Void,                              # 1
    String, Void,                              # 3
    UInt8, Void,                               # 5
    UInt16, Void,                              # 7
    UInt16, UInt16, UInt16, Void,              # 9  x, y, color
    UInt16, UInt16, UInt16, UInt, UInt8, Void, # 13 x, y, dir|color, orgId, infoBits
    Float64, Void                              # 19
  ]
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
cd(Pkg.dir("Latte"))
Pkg.add("Coverage")
using Coverage

folder = "src"
println("""Coverage.process_folder: Searching $folder for .jl files...""")
source_files = FileCoverage[]
files = readdir(folder)
for file in files
    fullfile = joinpath(folder,file)
    if isfile(fullfile)
        # Is it a Julia file?
        if splitext(fullfile)[2] == ".jl"
            push!(source_files, process_file(fullfile,folder))
        else
            println("Coverage.process_folder: Skipping $file, not a .jl file")
        end
    else isdir(fullfile)
        if file == "stdlib"
            # Skip coverage in stdlib because layer functions are never actually executed
            continue
        end
        # If it is a folder, recursively traverse
        append!(source_files, process_folder(fullfile))
    end
end

Coveralls.submit(source_files)


text:
using TSML
using TSMLextra
using Distributed
using Base.Threads
nthreads()

nprocs()==1 && addprocs()
nworkers()

@everywhere function predict(learner,data,train_ind,test_ind)
    features = convert(Matrix,data[:, 1:(end-1)])
    labels = convert(Array,data[:, end])
    # Create pipeline
    pipeline = Pipeline(
       Dict(
         :transformers => [
           OneHotEncoder(), # Encodes nominal features into numeric
           Imputer(), # Imputes NA values
           StandardScaler(),
           learner # Predicts labels on instances
         ]
       )
    )
    # Train
    fit!(pipeline, features[train_ind, :], labels[train_ind]);
    # Predict
    predictions = transform!(pipeline, features[test_ind, :]);
    # Assess predictions
    result = score(:accuracy, labels[test_ind], predictions)
    return result,pipeline
end

function parallelmodel(learners::Dict,data::DataFrame;trials=5)
    models=collect(keys(learners))
    ctable=@distributed (vcat) for i=1:trials
        # Split into training and test sets
        Random.seed!(3i)
        (train_ind, test_ind) = holdout(size(data, 1), 0.20)
        acc=@distributed (vcat) for model in models
            res,_=predict(learners[model],data,train_ind,test_ind)
            println("trial ",i,", ",model," => ",round(res))
            [model res i]
        end
        acc
    end
    df = ctable |> DataFrame
    rename!(df,:x1=>:model,:x2=>:acc,:x3=>:trial)
    gp=by(df,:model) do x
       DataFrame(mean=mean(x.acc),std=std(x.acc),n=length(x.acc))
    end
    sort!(gp,:mean,rev=true)
    return gp
end

@everywhere using TSML
@everywhere using TSMLextra

# Caret ML
@everywhere caret_svmlinear = CaretLearner(Dict(:learner=>"svmLinear"))
@everywhere caret_treebag = CaretLearner(Dict(:learner=>"treebag"))
@everywhere caret_rpart = CaretLearner(Dict(:learner=>"rpart"))
@everywhere caret_rf = CaretLearner(Dict(:learner=>"rf"))

# ScikitLearn ML
@everywhere sk_ridge = SKLearner(Dict(:learner=>"RidgeClassifier"))
@everywhere sk_sgd = SKLearner(Dict(:learner=>"SGDClassifier"))
@everywhere sk_knn = SKLearner(Dict(:learner=>"KNeighborsClassifier"))
@everywhere sk_gb = SKLearner(Dict(:learner=>"GradientBoostingClassifier",:impl_args=>Dict(:n_estimators=>10)))
@everywhere sk_extratree = SKLearner(Dict(:learner=>"ExtraTreesClassifier",:impl_args=>Dict(:n_estimators=>10)))
@everywhere sk_rf = SKLearner(Dict(:learner=>"RandomForestClassifier",:impl_args=>Dict(:n_estimators=>10)))

# Julia ML
@everywhere jrf = RandomForest(Dict(:impl_args=>Dict(:num_trees=>300)))
@everywhere jpt = PrunedTree()
@everywhere jada = Adaboost()

# Julia Ensembles
@everywhere jvote_ens=VoteEnsemble(Dict(:learners=>[jrf,jpt,sk_gb,sk_extratree,sk_rf]))
@everywhere jstack_ada=StackEnsemble(Dict(:stacker=>Adaboost(),:learners=>[jrf,jpt,sk_gb,sk_extratree,sk_rf]))
@everywhere jstack_rf=StackEnsemble(Dict(:learners=>[jrf,jpt,sk_gb,sk_extratree,sk_rf]))
@everywhere jbest_ens=BestLearner(Dict(:learners=>[jrf,sk_gb,sk_rf]))
@everywhere jsuper_ens=VoteEnsemble(Dict(:learners=>[jvote_ens,jstack_ada,jstack_rf,jbest_ens,sk_rf,sk_gb]))

learners=Dict(
      :jvote_ens=>jvote_ens,:jstack_rf=>jstack_rf,:jbest_ens=>jbest_ens, :jstack_ada=>jstack_ada,
      :jrf => jrf,:jada=>jada,:jsuper_ens=>jsuper_ens,:crt_rpart=>caret_rpart,
      :crt_svmlinear=>caret_svmlinear,:crt_treebag=>caret_treebag,:crt_rf=>caret_rf,
      :skl_knn=>sk_knn,:skl_gb=>sk_gb,:skl_extratree=>sk_extratree,
      :sk_rf => sk_rf
);

# use iris dataset
using RCall
iris = R"iris"|> rcopy
first(iris,5)

df = parallelmodel(learners,iris;trials=3)

(train_ind, test_ind) = holdout(size(iris, 1), 0.20)
res,workflow=predict(learners[:jstack_rf],iris,train_ind,test_ind)
res

showtree(workflow)

function threadedmodel(learners::Dict,data::DataFrame;trials=5)
    models=collect(keys(learners))
    global ctable = DataFrame()
    mtx = ReentrantLock()
    @threads for i=1:trials
        # Split into training and test sets
        Random.seed!(3i)
        (train_ind, test_ind) = holdout(size(data, 1), 0.20)
        @threads for model in models
            res=predict(learners[model],data,train_ind,test_ind)
            println("thread ",threadid(),", ",model," => ",round(res))
            lock(mtx)
            global ctable = vcat(ctable,DataFrame(model=model,acc=res))
            unlock(mtx)
        end
    end
    df = ctable |> DataFrame
    gp=by(df,:model) do x
       DataFrame(mean=mean(x.acc),std=std(x.acc),n=nrow(x))
    end
    sort!(gp,:mean,rev=true)
    return gp
end

df = threadedmodel(learners,tsdata;trials=5)

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
module CUTENSOR

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cudaDataType
using ..CUDA: libcutensor,  @retry_reclaim

using CEnum

const cudaDataType_t = cudaDataType

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

# thread cache for task-local library handles
const thread_handles = Vector{Union{Nothing,Base.RefValue{cutensorHandle_t}}}()

function handle()
    tid = Threads.threadid()
    if @inbounds thread_handles[tid] === nothing
        ctx = context()
        thread_handles[tid] = get!(task_local_storage(), (:CUTENSOR, ctx)) do
            handle = Ref{cutensorHandle_t}()
            cutensorInit(handle)
            handle
        end
    end
    something(@inbounds thread_handles[tid])
end

function __init__()
    resize!(thread_handles, Threads.nthreads())
    fill!(thread_handles, nothing)

    CUDA.atdeviceswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end

    CUDA.attaskswitch() do
        tid = Threads.threadid()
        thread_handles[tid] = nothing
    end
end

end


text:
mutable struct RunningStats
  num_infected::Int
  num_detected::Int
  num_dead::Int
  #daily_detections::Vector{UInt32}

  RunningStats() = new(0, 0, 0)
end

numinfected(stats::RunningStats) = stats.num_infected
numdetected(stats::RunningStats) = stats.num_detected
numdead(stats::RunningStats) = stats.num_dead
#dailydetections(stats::RunningStats) = stats.daily_detections

function update!(stats::RunningStats, event::Event)
  ek = kind(event)
  if isdetection(ek)
    stats.num_detected += 1

    #day = floor(UInt, time(event)) + 1
    #while day > length(stats.daily_detections)
    #  push!(stats.daily_detections, 0)
    #end
    #stats.daily_detections[day] += 1

  elseif istransmission(ek)
    stats.num_infected += 1
  elseif isdeath(ek)
    stats.num_dead += 1
  end
  nothing
end

function reset!(stats::RunningStats)
  stats.num_infected = 0
  stats.num_detected = 0
  stats.num_infected = 0
  #empty!(stats.daily_detections)
end

function sizehint!(stats::RunningStats, days::Integer)
  sizehint!(stats.daily_detections, days)
end


text:
struct Rect
    top::Int
    left::Int
    height::Int
    width::Int
    resolution::Tuple{Int, Int}
end
Rect((top, left), (height, width); resolution = (0, 0)) = Rect(
    top, left, height, width, resolution
)

struct SubImage{T, TSA <: SubArray{T}} <: AbstractMatrix{T}
    image::TSA
    rect::Rect
end
@forward SubImage.image Base.setindex!, Base.getindex, Base.length, Base.size
image(s::SubImage) = s.image
Base.copy(s::SubImage) = copy(s.image)

"""
    subimage(image, (top, left), (height, width); resolution = size(image))
    subimage(image, rect)

Return a sub-image of the image. The rect is in the given `resolution`.
"""
function subimage(img::AbstractMatrix, rect::Rect)
    resolution = rect.resolution
    top, left, height, width = if resolution ≠ (0, 0)
        ih, iw = size(parent(img))
        dh, dw = resolution
        rect.top * ih ÷ dh,
        rect.left * iw ÷ dw,
        rect.height * ih ÷ dh,
        rect.width * iw ÷ dw
    else
        rect.top, rect.left, rect.height, rect.width
    end
    h, w = size(img)
    top = max(top, 1)
    bottom = min(top + height - 1, h)
    left = max(left, 1)
    right = min(left + width - 1, w)
    s = @view img[top:bottom, left:right]
    SubImage(s, rect)
end

subimage(
    img::AbstractMatrix,
    (top, left),
    (height, width);
    resolution = (0, 0)
) = subimage(
    img,
    Rect((top, left), (height, width), resolution = resolution)
)

subimage(img::SubImage, rect::Rect) = let pr = img.rect
    @assert rect.resolution == (0, 0) || pr.resolution == rect.resolution
    subimage(
        parent(img.image),
        (pr.top + rect.top - 1, pr.left + rect.left - 1),
        (rect.height, rect.width),
        resolution = pr.resolution
    )
end

Base.getindex(img::AbstractMatrix, rect::Rect) = subimage(img, rect)
Base.getindex(img::SubImage, rect::Rect) = subimage(img, rect)
Base.getindex(img::VsFrame, rect::Rect) = subimage(img, rect)
Base.size(rect::Rect) = (rect.height, rect.width)
topleft(rect::Rect) = (rect.top, rect.left)

is_gray(img::AbstractMatrix{<:Gray}) = true
is_gray(img::AbstractMatrix) = false


text:
import DataFrames
# using InvertedIndices
import Luxor

import Base
import CSV
import Printf

import Combinatorics
import Statistics
import StatsBase
import Random
import Distributions
import KernelDensity
import SparseArrays

import LinearAlgebra, LinearAlgebra.I

include(joinpath(WATERFALL_DIR,"src","Figure","Color","Coloring.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Color","Blending.jl"))

abstract type Shape end
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Box.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Line.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Point.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Poly.jl"))

include(joinpath(WATERFALL_DIR,"src","Figure","Text","Label.jl"))

abstract type Sampling end
abstract type Geometry <: Sampling end
include(joinpath(WATERFALL_DIR,"src","Cascade","Data.jl"))
include(joinpath(WATERFALL_DIR,"src","Cascade","Cascade.jl"))
include(joinpath(WATERFALL_DIR,"src","Cascade","Geometry","Horizontal.jl"))
include(joinpath(WATERFALL_DIR,"src","Cascade","Geometry","Parallel.jl"))
include(joinpath(WATERFALL_DIR,"src","Cascade","Geometry","Vertical.jl"))
include(joinpath(WATERFALL_DIR,"src","Cascade","Geometry","Violin.jl"))

include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Arrow.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Box.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Line.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Point.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Shape","Poly.jl"))

abstract type Axis end
include(joinpath(WATERFALL_DIR,"src","Figure","XAxis.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","YAxis.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Annotation.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Handle.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Legend.jl"))
include(joinpath(WATERFALL_DIR,"src","Figure","Plot.jl"))

include(joinpath(WATERFALL_DIR,"src","definitions.jl"))
include(joinpath(WATERFALL_DIR,"src","options.jl"))

include(joinpath(WATERFALL_DIR,"src","utils","attributes.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","io.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","calc.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","draw.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","random.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","scale.jl"))
include(joinpath(WATERFALL_DIR,"src","utils","utils.jl"))

include(joinpath(WATERFALL_DIR,"src","utils","define.jl"))


text:
export Synapse, Segment, Neuron, NeuronOrSegment,NeuronOrSegmentOrOutput,AbstractNetwork,AbstractSynapse, Network, Input, Logger, BernoulliSynapseWeight
using DataFrames: DataFrame

################################################################################
## model component definitions                                                ##
################################################################################

abstract type AbstractComponent{ID,T,WT,IT} end
abstract type NeuronOrSegmentOrOutput{ID,T,WT,IT} <: AbstractComponent{ID,T,WT,IT} end
abstract type NeuronOrSegment{ID,T,WT,IT} <: NeuronOrSegmentOrOutput{ID,T,WT,IT} end
abstract type AbstractNetwork{ID,T,WT,IT} <:AbstractComponent{ID,T,WT,IT} end
abstract type AbstractSynapse{ID,T,WT,IT} <:AbstractComponent{ID,T,WT,IT} end

@enum VoltageLevel voltage_low=0 voltage_elevated=1 voltage_high=2

mutable struct Segment{ID,T,WT,IT} <: NeuronOrSegment{ID,T,WT,IT}
    id::ID
    θ_syn::IT
    θ_seg::Int
    plateau_duration::T

    state_syn::IT
    state_seg::Int
    state::VoltageLevel
    active::Bool
    last_activated::T

    next_downstream::NeuronOrSegment{ID,T,WT,IT}
    next_upstream::Vector{Segment{ID,T,WT,IT}}
    net::AbstractNetwork{ID,T,WT,IT}

    function Segment(id::ID, root::NeuronOrSegment{ID,T,WT,IT}; θ_syn=1, θ_seg=1, plateau_duration=root.net.default_plateau_duration) where {ID,T,WT,IT}
        this = new{ID,T,WT,IT}(id, θ_syn, θ_seg, plateau_duration, zero(IT), 0, voltage_low, false, zero(T), root, Segment{ID,T,WT,IT}[], root.net)
        push!(root.next_upstream, this)
        return this
    end
end

mutable struct Synapse{ID,T,WT,IT} <: AbstractSynapse{ID,T,WT,IT}
    id::ID
    target::NeuronOrSegmentOrOutput{ID,T,WT,IT}
    delay::T
    spike_duration::T
    weight::WT
    magnitude::IT
    state::Bool

    net::AbstractNetwork{ID,T,WT,IT}

    function Synapse(
                id::ID, source, target::NeuronOrSegmentOrOutput{ID,T,WT,IT};
                delay::T = source.net.default_delay, spike_duration::T = source.net.default_spike_duration, weight::WT=one(WT)
            ) where {ID,T,WT,IT}
        this = new{ID,T,WT,IT}(id, target, delay, spike_duration, weight, zero(IT), false, source.net)
        push!(source.synapses, this)
        return this
    end
end
sample_spike_magnitude(weight)=weight

"""
Fallback: let YAML.jl handle it.
"""
serialize(w) = w

mutable struct Neuron{ID,T,WT,IT} <: NeuronOrSegment{ID,T,WT,IT}
    id::ID
    θ_syn::IT
    θ_seg::Int
    refractory_duration::T

    state_syn::IT
    state_seg::Int
    state::VoltageLevel
    active::Bool

    next_upstream::Vector{Segment{ID,T,WT,IT}}
    synapses::Vector{Synapse{ID,T,WT,IT}}
    net::AbstractNetwork{ID,T,WT,IT}

    function Neuron(id::ID, net::AbstractComponent{ID,T,WT,IT}; θ_syn=one(IT), θ_seg=1, refractory_duration::T=net.default_refractory_duration) where {ID,T,WT,IT}
        this = new{ID,T,WT,IT}(id, θ_syn, θ_seg, refractory_duration, zero(IT), 0, voltage_low, false, Segment{ID,T,WT,IT}[], Synapse{ID,T,WT,IT}[], net)
        push!(net.neurons, this)
        return this
    end
end

mutable struct Input{ID,T,WT,IT} <: AbstractComponent{ID,T,WT,IT}
    id::ID
    synapses::Vector{Synapse{ID,T,WT,IT}}
    net::AbstractNetwork{ID,T,WT,IT}

    function Input(id::ID, net::AbstractComponent{ID,T,WT,IT}) where {ID,T,WT,IT}
        this = new{ID,T,WT,IT}(id::ID,Synapse{ID,T,WT,IT}[],net)
        push!(net.inputs, this)
        return this
    end
end

mutable struct Output{ID,T,WT,IT} <: NeuronOrSegmentOrOutput{ID,T,WT,IT}
    id::ID
    net::AbstractNetwork{ID,T,WT,IT}
    state_syn::IT

    function Output(id::ID, net::AbstractComponent{ID,T,WT,IT}) where {ID,T,WT,IT}
        this = new{ID,T,WT,IT}(id::ID, net, zero(IT))
        push!(net.outputs, this)
        return this
    end
end

mutable struct Network{ID,T,WT,IT} <: AbstractNetwork{ID,T,WT,IT}
    inputs::Vector{Input{ID,T,WT,IT}}
    outputs::Vector{Output{ID,T,WT,IT}}
    neurons::Vector{Neuron{ID,T,WT,IT}}

    default_refractory_duration::T
    default_delay::T
    default_spike_duration::T
    default_plateau_duration::T

    Network(;
        id_type::Type=Symbol, time_type::Type=Float64, weight_type::Type=Int, synaptic_input_type::Type=Int,
        default_refractory_duration = 1.0, default_delay = 1.0, default_spike_duration = 5.0, default_plateau_duration = 100.0
    ) = new{id_type, time_type, weight_type, synaptic_input_type}(
            Input{id_type, time_type, weight_type, synaptic_input_type}[],
            Output{id_type, time_type, weight_type, synaptic_input_type}[],
            Neuron{id_type, time_type, weight_type, synaptic_input_type}[],
            default_refractory_duration, default_delay, default_spike_duration, default_plateau_duration
        )
end

################################################################################
## default simulation handling                                                ##
################################################################################

struct Event{T, V}
    type::Symbol
    created::T
    t::T
    target::V
end
Base.isless(ev1::Event, ev2::Event) = isless(ev1.t, ev2.t)

function handle!(event::Event, queue!, logger!)
    if event.type == :input_spikes
        # definitely turn input on
        # -> this will trigger all synapses to fire
        on!(event.target, event.t, queue!, logger!)
    elseif event.type == :epsp_starts
        # definitively turn synapse on -> this may trigger a cascade of segements / soma to emit plateaus / a spike
        on!(event.target, event.t, queue!, logger!)
    elseif event.type == :epsp_ends
        # definitively turn synapse off
        off!(event.target, event.t, queue!, logger!)
    elseif event.type == :plateau_ends
        # if the event is not relevant any more (the current plateau is not the one that triggered this event), ignore
        # if the event is still relevant, definitely turn the segment off
        # -> this may also turn off upstream segments
        if event.created >= event.target.last_activated
            maybe_off!(event.target, event.t, queue!, logger!)
        end
    elseif event.type == :refractory_period_ends
        # definitively enable neuron to fire again
        # -> this may also trigger another spike rightaway
        # -> a spike will trigger all synapses to fire
        maybe_off!(event.target, event.t, queue!, logger!)
    end
    nothing
end

################################################################################
## model behavior                                                             ##
################################################################################

function reset!(net::Network{ID,T,WT,IT}) where {ID,T,WT,IT}
    foreach(reset!, net.inputs)
    foreach(reset!, net.outputs)
    foreach(reset!, net.neurons)
end

function reset!(x::Neuron{ID,T,WT,IT}) where {ID,T,WT,IT}
    x.state_syn = zero(IT)
    x.state_seg = 0
    x.active = false
    x.state = voltage_low

    foreach(reset!, x.synapses)
    foreach(reset!, x.next_upstream)
    update_state!(x)
end

function reset!(x::Segment{ID,T,WT,IT}) where {ID,T,WT,IT}
    x.state_syn = zero(IT)
    x.state_seg = 0
    x.active = false
    x.state = voltage_low
    x.last_activated = zero(T)

    foreach(reset!, x.next_upstream)
    update_state!(x)
end

function reset!(x::Input{ID,T,WT,IT}) where {ID,T,WT,IT}
    foreach(reset!, x.synapses)
end

function reset!(x::Output{ID,T,WT,IT}) where {ID,T,WT,IT}
    x.state_syn = zero(IT)
end

function reset!(x::Synapse{ID,T,WT,IT}) where {ID,T,WT,IT}
    x.state = false
    x.magnitude = zero(IT)
end

"""Recursively update all upstream branches."""
function backprop!(obj::Segment, now, logger!)
    # update this element
    recurse = update_state!(obj)
    logger!(now, :backprop, obj.id, obj.state)

    # if necessary, also update upstream elements
    if recurse
        for s_up in obj.next_upstream
            backprop!(s_up, now, logger!)
        end
    end
    return nothing
end

function update_state!(obj::Segment)
    old_state = obj.state
    obj.state = if obj.active || (isa(obj.next_downstream, Segment) && obj.next_downstream.state == voltage_high)
        # active or passively held high by downstream segment
        voltage_high
    elseif obj.state_seg >= obj.θ_seg
        # elevated voltage due to enough upstream input
        voltage_elevated
    else
        voltage_low
    end

    return old_state != obj.state
end

"""Check if this segment was just turned on; if so, backpropagate!"""
function maybe_on!(obj::Segment, now, queue!, logger!)
    update_state!(obj)
    return if obj.active && obj.state_syn >= obj.θ_syn #&& (isempty(obj.next_upstream) || obj.state == voltage_elevated)
        # if the segment is already on, but there is sufficient EPSP, re-trigger!
        @debug "$(now): Re-triggered segment $(obj.id) and extended its plateau!"
        logger!(now, :plateau_extended, obj.id, obj.state)

        # Re-schedule plateau-end time
        obj.last_activated = now
        queue!(Event(:plateau_ends, now, now+obj.plateau_duration, obj))
        true

    elseif ~obj.active && obj.state_syn >= obj.θ_syn && (isempty(obj.next_upstream) || obj.state == voltage_elevated)
        # if the segment is currently off, but the plateau conditions are satisfied, turn on
        @debug "$(now): Triggered segment $(obj.id) and it turned on!"

        # The segment is now active
        obj.active = true

        # propagate to the upstream segments & update
        backprop!(obj, now, logger!)
        logger!(now, :plateau_starts, obj.id, obj.state)

        # turning on this segment might have also turned on the next_downstream segment
        obj.next_downstream.state_seg += 1
        updated = update_state!(obj.next_downstream)
        if ~maybe_on!(obj.next_downstream, now, queue!, logger!) && updated
            logger!(now, :upstate_start, obj.next_downstream.id, obj.next_downstream.state)
        end

        # Each segment must turn off
        obj.last_activated = now
        queue!(Event(:plateau_ends, now, now+obj.plateau_duration, obj))
        true
    else
        @debug "$(now): Triggered segment $(obj.id), but it didn't turn on!"
        false
    end
end

"""Check if this segment should switch off, since its plateau has timed out"""
function maybe_off!(obj::Segment, now, queue!, logger!)
    # if the segment is currently on and has no reason to stay on, turn off

    if obj.active
        @debug "$(now): Triggered segment $(obj.id) and it turned off!"

        # The segment is now inactive
        obj.active = false

        # propagate to the upstream segments & update
        backprop!(obj, now, logger!)
        logger!(now, :plateau_ends, obj.id, obj.state)

        # update next_downstream segment
        obj.next_downstream.state_seg -= 1
        if update_state!(obj.next_downstream)
            logger!(now, :upstate_ends, obj.next_downstream.id, obj.next_downstream.state)
        end
    else
        @debug "$(now): Triggered segment $(obj.id), but it didn't turn off!"
    end
    return nothing
end

"""Check if this Input was triggered to spike"""
function on!(obj::Input, now, queue!, logger!)
    @debug "$(now): Triggered input $(obj.id) to fire!"
    for s in obj.synapses
        queue!(Event(:epsp_starts, now, now+s.delay, s))
    end
    return nothing
end

function update_state!(obj::Neuron)
    old_state = obj.state
    obj.state = if obj.active
        voltage_high
    elseif obj.state_seg >= obj.θ_seg
        # elevated voltage due to enough upstream input
        voltage_elevated
    else
        voltage_low
    end

    return old_state != obj.state
end

"""Check if this neuron was triggered to spike"""
function maybe_on!(obj::Neuron, now, queue!, logger!)
    update_state!(obj)

    return if ~obj.active && obj.state_syn >= obj.θ_syn && (isempty(obj.next_upstream) || obj.state == voltage_elevated)
        @debug "$(now): Triggered neuron $(obj.id) and it fired a spike!"

        obj.active = true
        obj.state = voltage_high
        logger!(now, :spikes, obj.id, obj.state)

        # trigger all synapses
        for s in obj.synapses
            queue!(Event(:epsp_starts, now, now+s.delay, s))
        end

        # don't forget to recover from refractory period
        queue!(Event(:refractory_period_ends, now, now+obj.refractory_duration, obj))
        true
    else
        @debug "$(now): Triggered neuron $(obj.id), but didn't fire!"
        false
    end
end

"""Check if this neuron was re-triggered to spike after refractoriness"""
function maybe_off!(obj::Neuron, now, queue!, logger!)
    return if obj.active
        @debug "$(now): Neuron $(obj.id) came out of refractoriness!"
        obj.active = false
        update_state!(obj)
        logger!(now, :refractory_period_ends, obj.id, obj.state)

        # check if the neuron should spike again, right away
        maybe_on!(obj::Neuron, now, queue!, logger!)
    else
        @debug "$(now): Triggered Neuron $(obj.id), but was already out of refractoriness!"
        false
    end
end

"""Check if an incoming spike should trigger an EPSP"""
function on!(obj::Synapse, now, queue!, logger!)
    if ~obj.state
        # set the synapse's state for this spike
        obj.magnitude = sample_spike_magnitude(obj.weight)
        if obj.magnitude ≈ zero(obj.magnitude)
            @debug "$(now): Triggered synapse $(obj.id), but didn't start an EPSP!"
        else
            @debug "$(now): Triggered synapse $(obj.id) and started an EPSP!"
            obj.state = true
            # inform the target segment about this new EPSP
            obj.target.state_syn += obj.magnitude
            logger!(now, :epsp_starts, obj.id, obj.state)

            # don't forget to turn off EPSP
            queue!(Event(:epsp_ends, now, now+obj.spike_duration, obj))

            # if the spike was excitatory this EPSP might have triggered a plateau in the target
            # else if it was inhibitory this EPSP may have destroyed any ongoing plateau
            if obj.magnitude > zero(obj.magnitude)
                maybe_on!(obj.target, now,  queue!, logger!)
            elseif obj.magnitude < zero(obj.magnitude)
                maybe_off!(obj.target, now,  queue!, logger!)
            end
        end
    else
        @debug "$(now): Triggered synapse $(obj.id), but synapse was already active!"
    end
    return nothing
end

"""Turn off the EPSP"""
function off!(obj::Synapse, now, queue!, logger!)
    # only update if the synapse isn't already off (shouldn't happen!)
    if obj.state
        @debug "$(now): Triggered synapse $(obj.id) and turned EPSP off!"
        # inform the target, that the EPSP is over
        obj.target.state_syn -= obj.magnitude

        # if the spike was inhibitory, this EPSP's end may trigger a rebound spike
        if obj.magnitude < zero(obj.magnitude)
            maybe_on!(obj.target, now,  queue!, logger!)
        end

        # reset the synapse's state
        obj.magnitude = zero(obj.magnitude)
        obj.state = false
        logger!(now, :epsp_ends, obj.id, obj.state)

    else
        @debug "$(now): Triggered synapse $(obj.id), but didn't turn EPSP off!"
    end
    return nothing
end

"""Trigger output"""
function maybe_on!(obj::Output, now, queue!, logger!)
    @debug "$(now): Triggered output $(obj.id)!"
    logger!(now, :spikes, obj.id, true)
    return nothing
end

################################################################################
## Logging                                                                    ##
################################################################################

struct Logger
    data::DataFrame
    filter::Function

    function Logger(net::Network{ID,T,WT,IT}; filter=(t,tp,id,x)->true, DT=Union{Bool,VoltageLevel}) where {ID,T,WT,IT}
        data = DataFrame(:t=>T[], :event=>Symbol[], :object=>ID[], :state=>DT[])
        new(data, filter)
    end
end

"""Log state of 'Output' object at each event"""
function (l::Logger)(args...)
    if getfield(l, :filter)(args...)
        push!(getfield(l,:data), args)
    end
end

# Base.getproperty(l::Logger, sym::Symbol) = getproperty(getfield(l,:data),sym)
# Base.propertynames(l::Logger) = propertynames(getfield(l,:data))

################################################################################
## pretty printing                                                            ##
################################################################################

Base.show(io::IO, x::Logger)   = show(io, getfield(x,:data))
Base.show(io::IO, x::Network)  = print(io, "Network with $(length(x.inputs)) inputs and $(length(x.neurons)) neurons")
Base.show(io::IO, x::Neuron)   = print(io, "Neuron '$(x.id)' with $(length(x.next_upstream)) child-segments and $(length(x.synapses)) outgoing synapses")
Base.show(io::IO, x::Input)    = print(io, "Input '$(x.id)' with $(length(x.synapses)) outgoing synapses")
Base.show(io::IO, x::Output)   = print(io, "Output '$(x.id)'")
Base.show(io::IO, x::Segment)  = print(io, "Segment '$(x.id)' with $(length(x.next_upstream)) child-segments")
Base.show(io::IO, x::Synapse)  = print(io, "Synapse '$(x.id)' (connects to $(x.target.id))")

################################################################################
## Probabilistic synapse implementation                                       ##
################################################################################

struct BernoulliSynapseWeight{T}
    magnitude::T
    probability::Float64
end
Base.:>=(w::BernoulliSynapseWeight, other::BernoulliSynapseWeight) = w.magnitude >= other.magnitude
Base.:>=(w::BernoulliSynapseWeight, x) = w.magnitude >= x
Base.one(::Type{BernoulliSynapseWeight{T}}) where T = BernoulliSynapseWeight(one(T), 1.0)
# Add parser for just a number
BernoulliSynapseWeight{T}(x::T) where {T} = BernoulliSynapseWeight(x, 1.0)
# Add parser for a tuple or vector (magnitude,probability)
BernoulliSynapseWeight{T}(x::Union{Tuple,Vector}) where {T} = BernoulliSynapseWeight(x...)
# generate a stochastic EPSP magnitude
sample_spike_magnitude(w::BernoulliSynapseWeight{T}) where {T} = (w.probability ≥ 1.0 || rand() ≤ w.probability) ? w.magnitude : zero(T)
# add serialization
serialize(w::BernoulliSynapseWeight{T}) where {T} = [w.magnitude, w.probability]


text:
using GymSpaces: Discrete, Box

mutable struct MountainCarEnv <: AbstractEnv
    min_position::Float32
    max_position::Float32
    max_speed::Float32
    goal_position::Float32

    force::Float32
    gravity::Float32

    low::Array{Float32, 1}
    high::Array{Float32, 1}

    action_space::Discrete
    observation_space::Box
    state
end

include("vis/mountaincar.jl")

function MountainCarEnv()
    min_position = -1.2f0
    max_position = 0.6f0
    max_speed = 0.07f0
    goal_position = 0.5f0

    force = 0.001f0
    gravity = 0.0025f0

    low = [min_position, -max_speed]
    high = [max_position, max_speed]

    action_space = Discrete(5)
    observation_space = Box(low, high, Float32)
    state

    MountainCarEnv(min_position, max_position, max_speed, goal_position,
                   force, gravity,
                   low, high,
                   action_space, observation_space, nothing)
end

function step!(env::MountainCarEnv, action)
    @assert action ∈ env.action_space "$(action) is unavailable for this environment."

    position, velocity = env.state[1:1], env.state[2:2]
    v = velocity .+ (action.-2)*env.force .+ cos.(3position)*(-env.gravity)
    velocity_ = clamp.(v, -env.max_speed, env.max_speed)
    x = position .+ velocity_
    position_  = clamp.(x, env.min_position, env.max_position)
    if all(position_ .== env.min_position) && all(velocity_ .< 0)
        velocity_ = 0f0velocity_
    end

    done = all(position_ .≥ env.goal_position)
    reward = -1f0

    env.state = vcat(position_, velocity_)
    return env.state, reward, done, Dict()
end

function reset!(env::MountainCarEnv)
    env.state = param([2f-1rand(Float32) - 6f-1, 0f0])
end

Base.show(io::IO, env::MountainCarEnv) = print(io, "MountainCarEnv")


text:
struct RemoteFunction <: Function
    remote_identifier::Int64
    local_function::Function

    function RemoteFunction(remote_identifier::Int64, func::Function)
        new(remote_identifier, func)
    end
end

(r::RemoteFunction)(args...) = r.local_function(args...)

macro remote(pool, ex)
    quote
        register!($(esc(pool)), $(esc(ex)))
    end
end

Serialization.serialize(s::AbstractSerializer, rf::RemoteFunction) = serialize(s, rf.remote_identifier)


text:
"""
    $TYPEDEF

$TYPEDFIELDS
"""
@with_kw mutable struct History
    # Number of completed training epochs
    epochs::Int = 0
    # total number of completed steps
    steps::Int = 0
    # number of completed steps in current epoch
    stepsepoch::Int = 0
    # number of seen samples during training, usually `nsteps * batchsize`
    samples::Int = 0
end

"""
    Recorder(epoch, step, epochstats, stepstats)
"""
struct Recorder <: Callback end

stateaccess(::Recorder) = (
    cbstate = (history = Write(),),
    batch = Read(),
    )

function on(::Init, ::Phase, recorder::Recorder, learner)
    learner.cbstate.history = History()
end

function on(::EpochBegin, ::AbstractTrainingPhase, recorder::Recorder, learner)
    learner.cbstate.history.stepsepoch = 0
end

function on(::BatchEnd, phase::AbstractTrainingPhase, recorder::Recorder, learner)
    history = learner.cbstate.history
    history.steps += 1
    history.stepsepoch += 1
    history.samples += size(learner.batch.xs)[end]
end

function on(::EpochEnd, phase::Phase, recorder::Recorder, learner)
    history = learner.cbstate.history
    if phase isa AbstractTrainingPhase
        history.epochs += 1
    end
end