text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m7_opt_rho6_4(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m7_opt_rho6_4!(A_copy)
end

@inline function exp_ps_m7_opt_rho6_4!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 T2k10
    # Computing Bb2 = x*I+x*A
    coeff1=0.42031771622120523
    coeff2=0.309252105381124
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.42031771622120523
    coeff2=0.309252105381124
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.09654809779701808
    coeff2=0.9314722195252554
    coeff3=0.19731022777550086
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.9062397883450928
    coeff2=0.11339252602614255
    coeff3=0.12302744324936216
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.06178357797354583
    coeff2=0.13477081084608053
    coeff3=0.7329752440067754
    coeff4=0.5453304515439815
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.7648819021317048
    coeff2=0.3110254108862011
    coeff3=0.13872927228389606
    coeff4=0.031072070853003377
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0037226968956481066
    coeff2=0.3027942237668864
    coeff3=0.16951280499357504
    coeff4=0.7173684223758208
    coeff5=0.6151562256858113
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.4984122280039027
    coeff2=0.6168934079023138
    coeff3=0.32464232547450916
    coeff4=0.11022586483948126
    coeff5=0.00562561032811143
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-2.47688825760804e-9
    coeff2=0.0004847383569867066
    coeff3=-0.008759615364757254
    coeff4=-0.024971360914667298
    coeff5=0.20141416233600964
    coeff6=0.9838034227912541
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-7.332920860083383e-8
    coeff2=7.255106883862399e-10
    coeff3=1.599720660809858e-11
    coeff4=5.386947932212966e-11
    coeff5=3.696432230206371e-12
    coeff6=1.7604767541951467e-14
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=5.241383560289876e-5
    coeff2=0.014487073789250026
    coeff3=0.01054251273835961
    coeff4=-0.12590635459023522
    coeff5=-0.35076815913299725
    coeff6=0.9784008568588043
    coeff7=4.756683002701541e-7
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.004732216406789163
    coeff2=2.4676441956771072e-5
    coeff3=1.805471807800756e-5
    coeff4=6.965400723989983e-6
    coeff5=1.4212591399083807e-6
    coeff6=1.2164967481538e-7
    coeff7=1.0000000002008125
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9900217594371604
    coeff2=-0.0426544708742059
    coeff3=0.0895038182070465
    coeff4=0.4332397850086141
    coeff5=0.6230790774477177
    coeff6=0.836798518267107
    coeff7=6.715240076479665e-7
    coeff8=0.03410244422083362
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.09171588796170577
    coeff2=0.01956469387438758
    coeff3=0.01590736865892108
    coeff4=0.010821193755227018
    coeff5=0.002975380215829149
    coeff6=0.005577510585906947
    coeff7=-7.0435154500906185e-6
    coeff8=1.0199299719677057
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing T2k10 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=0.7602626115113443
    coeff2=0.2868135868803526
    coeff3=0.36416792067469345
    coeff4=0.1470803112875071
    coeff5=0.301249740245868
    coeff6=0.05977870243645669
    coeff7=3.4088571170380866e-7
    coeff8=0.00148688702927498
    coeff9=1.021358342405418
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    return memslots1 # Returning T2k10
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
function build_jac_func(symjac,indvar_dict,params;params_from_function=true)
  Jex = :()
  for i in 1:size(symjac,1)
    for j in 1:size(symjac,2)
      ex = Meta.parse(string(symjac[i,j]))
      if typeof(ex) <: Expr
        ode_findreplace(ex,copy(ex),indvar_dict,params,params_from_function=params_from_function)
      else
        ex = ode_symbol_findreplace(ex,indvar_dict,params,params_from_function=params_from_function)
      end
      push!(Jex.args,:(internal_var___J[$i,$j] = $ex))
    end
  end
  Jex.head = :block
  push!(Jex.args,nothing)
  Jex
end

function build_tgrad_func(symtgrad,indvar_dict,params;params_from_function=false)
  tgradex = :()
  for i in 1:length(symtgrad)
    ex = Meta.parse(string(symtgrad[i]))
    if typeof(ex) <: Expr
      ode_findreplace(ex,copy(ex),indvar_dict,params,params_from_function=params_from_function)
    else
      ex = ode_symbol_findreplace(ex,indvar_dict,params,params_from_function=params_from_function)
    end
    push!(tgradex.args,:(internal_var___grad[$i] = $ex))
  end
  tgradex.head = :block
  push!(tgradex.args,nothing)
  tgradex
end

function build_p_funcs(paramfuncs,indvar_dict,params)
  pfuncs = Vector{Expr}(undef,length(params))
  params_type = typeof(params)
  for i in 1:length(params)
    pfunc = :()
    params_drop_cur = copy(params)
    deleteat!(params_drop_cur,i)
    for j in 1:length(paramfuncs[1])
      ex = paramfuncs[i][j]
      if typeof(ex) <: Expr
        ode_findreplace(ex,copy(ex),indvar_dict,params_drop_cur)
      else
        ex = ode_symbol_findreplace(ex,indvar_dict,params_drop_cur)
      end
      push!(pfunc.args,:(internal_var___du[$j] = $ex))
    end
    pfunc.head = :block
    push!(pfunc.args,nothing)
    pfuncs[i] = pfunc
  end
  pfuncs
end

function build_component_funcs(symex)
  funcs = Vector{Expr}(undef,0) # Get all of the functions for symbolic computation
  for (i,arg) in enumerate(symex.args)
    if i%2 == 0
      ex = arg.args[2]
      if (typeof(ex) <: Symbol) || (typeof(ex) <: Number)
        push!(funcs,:($ex*1))
      else # It's an expression, just push
        fix_ex = copy(ex)
        flip_mult!(fix_ex)
        push!(funcs,fix_ex)
      end
    end
  end
  funcs
end

ext: jl


text:
boson = @theory begin :c * :cdag => :cdag * :c + 1 end

distr_assoc = distrib(:*, :+) ∪ associativity(:*)

ident = @theory begin
    1 * x => x
end

function boson_expand(ex)
    params = SaturationParams(timeout=4)
    G = EGraph(ex)
    saturate!(G, boson ∪ distr_assoc, params)
    ex = extract!(G, astsize_inv)

    G = EGraph(ex)
    saturate!(G, distr_assoc ∪ ident, params)
    ex = extract!(G, astsize)
    # ex = rewrite(ex, ident)
    # println(ex)
    return ex
end

e1 = :(c * c * cdag * cdag)

# TODO add expected result
# println(boson_expand(e1))
# println(normalize_nocycle(boson_expand, e1))


text:
function IMPES!(sim::Simulation, res::Reservoir, kr::RelPerms, fluid::Fluids, grid::Grid,
                bc::BoundaryConditions, post::FlowResults)

    # Define flow functions
    # krw(s) = water_rel_perm(s, res.kr)
    # kro(s) = oil_rel_perm(s, res.kr)

    # λw(s) = krw(s) / fluid.μw
    # λo(s) = kro(s) / fluid.μo
    # λ(s) = λo(s) + λw(s)

    # fw(s) = λw(s) ./ (λw(s) .+ λo(s))

    # Get data
    n = grid.Nx
    A = res.A
    k = res.k
    ϕ = res.phi
    Δx = grid.dx
    Δt = sim.dt
    s = sim.s
    θ = sim.θ

    # Initialize velocity vectors
    v = zeros(n+1)
    v[1] = bc.value[1] / A

    ∇v = zeros(n)

    # Calculate all mobilities
    krw, kro = kr(s, θ)
    Λw = krw ./ fluid.μw
    Λo = kro ./ fluid.μo
    Λ = Λw + Λo

    # Create T Vector
    T = zeros(n,n)
    for i in 2:n-1
        T[i,i-1] = Λ[i-1]
        T[i,i+1] = Λ[i]
        T[i,i] = - (T[i, i-1] + T[i,i+1])
    end

    # Left - Neumann
    T[1,1] = -Λ[1]
    T[1,2] = Λ[1]
    # Right - Dirichlet
    T[n,n] = - 9*Λ[n] - 3*Λ[n-1]
    T[n,n-1] = 3*Λ[n-1] + Λ[n]

    ## Get Q vector
    Q = zeros(n)
    # Left - Neumann
    Q[1] = - bc.value[1] / A / k * Δx
    # Right - Dirichlet
    Q[n] = -8 * (Λw[n] + Λo[n]) * bc.value[2]

    ## Solve for Pressure
    p = T \ Q

    ## Solve for Saturation
    # Calculate water velocity
    v[2:end-1] = [Λw[l-1] * (p[l-1] - p[l]) for l in 2:n] / Δx * k
    v[end] = 2 * Λw[n]*(p[n]-bc.value[2]) / Δx * k
    # Calculate gradient of water velocity
    ∇v = [v[l-1] - v[l] for l in 2:n+1] / Δx

    s = s .+ ∇v / ϕ * Δt

    # Update Simulation
    sim.p = p
    sim.s = s
    sim.v = v
    sim.i += 1
    sim.t += sim.dt

    # Save Results
    post.p[:,sim.i] = p
    post.s[:,sim.i] = s
    post.θ[:,sim.i] = θ

end

function IMPEC!(ctrans::CompTransport, sim, res, grid, post)
    n = grid.Nx
    Δx = grid.dx
    ϕ = res.phi
    Δt = sim.dt
    v = sim.v
    s = sim.s
    i = sim.i
    Δs = sim.s .- post.s[:,i-1]

    ncomps = ctrans.ncomps

    for j in 1:ncomps
        c = ctrans.c[:,j,i-1]

        cinj = ctrans.cinj[j]
        ∇cuw  = zeros(n)
        ∇cuw[2:n] = [c[l] * v[l+1] - c[l-1] * v[l]  for l in 2:n] / Δx
        ∇cuw[1] = (c[1] * v[2] - cinj * v[1])  / Δx

        c = c  .- Δt * ∇cuw / ϕ ./ s - c .* Δs ./ s

        ctrans.c[:,j,i] = c
    end

end

# function IMPEC!(ctrans::CompTransport, sim::Simulation,  res::Reservoir, fluid::Fluids, grid::Grid, bc::BoundaryConditions, wells::Wells)
#     Nx = grid.Nx
#     Ny = grid.Ny
#     p = sim.p
#     s = sim.s
#     ϕ = res.phi
#     k = res.k
#     Δx = grid.dx
#     Δt = sim.dt
#     μw = fluid.μw
#     Vi = res.h * res.W * Δx
#     qinj = wells.Q[1]
#     kr = res.kr
#     c = ctrans.c
#     cinj = ctrans.cinj
#     sold = ctrans.s_old

#     η = k * Δt / (ϕ * μw * Δx^2) * 6894757 * 1.06e-14 * 3600 * 24
#     A = zeros(Nx * Ny, Nx * Ny)
#     krw = [water_rel_perm(s[l], kr) for l in 1:Nx]

#     # Upwinding
#     for l in 2:Nx-1
#         # krw1 = water_rel_perm(s[l], kr)
#         # krw0 = water_rel_perm(s[l-1], kr)

#         A[l, l-1] =  krw[l-1] * (p[l-1] - p[l]) * η / s[l]
#         A[l, l] = krw[l] * (p[l+1] - p[l])  * η / s[l]
#     end

#     A[Nx, Nx-1] =  krw[Nx-1] * (p[Nx-1] - p[Nx]) * η / s[Nx]
#     A[Nx, Nx] =  - krw[Nx] * (2* p[Nx]) * η  / s[Nx]

#     A[1, 1] =  krw[1] * (p[2] - p[1]) * η  / s[1]

#     Q = zeros(Nx * Ny)
#     Q[1] = qinj * Δt / (s[1] * Vi * ϕ) * cinj

#     Δs = (s .- sold) ./ sold

#     ctrans.c = c + A * c + Q - Δs .* c
#     ctrans.s_old = sim.s
# end

# function IMPEC2!(ctrans::CompTransport, sim::Simulation,  res::Reservoir, fluid::Fluids, grid::Grid, bc::BoundaryConditions, wells::Wells)
#     Nx = grid.Nx
#     Ny = grid.Ny
#     p = sim.p
#     s = sim.s
#     ϕ = res.phi
#     k = res.k
#     Δx = grid.dx
#     Δt = sim.dt
#     μw = fluid.μw
#     Vi = res.h * res.W * Δx
#     qinj = wells.Q[1]
#     kr = res.kr
#     c = ctrans.c
#     cinj = ctrans.cinj
#     sold = ctrans.s_old

#     η = k * Δt / (ϕ * μw * Δx^2) * 6894757 * 1.06e-14 * 3600 * 24
#     A = zeros(Nx * Ny, Nx * Ny)

#     krw = [water_rel_perm(s[l], kr) for l in 1:Nx]

#     u = zeros(Nx)
#     for l in 2:Nx-1
#         u[l] = krw[l] * (p[l] - p[l+1])
#     end
#     u[1] = krw[1] * (p[2] - p[1])
#     u[Nx] = krw[Nx] * (p[Nx] - p[Nx-1])

#     # Upwinding
#     for l in 3:Nx - 1
#         A[l, l] = -3 * u[l] * η / 2 / s[l]
#         A[l, l-1] =  4 * u[l-1] * η / 2 / s[l]
#         A[l, l-2] =  -u[l-2] * η / 2 / s[l]
#     end

#     A[2, 1] = krw[1] * (p[1] - p[2]) * η / s[2]
#     A[2, 2] = krw[2] * (p[3] - p[2])  * η / s[2]

#     A[1, 1] =  krw[1] * (p[2] - p[1]) * η  / s[1]

#     A[Nx, Nx-1] =  krw[Nx-1] * (p[Nx-1] - p[Nx]) * η / s[Nx]
#     A[Nx, Nx] =  - krw[Nx] * (2* p[Nx]) * η  / s[Nx]

#     Q = zeros(Nx * Ny)
#     Q[1] = qinj * Δt / (s[1] * Vi * ϕ) * cinj

#     Δs = (s .- sold) ./ sold

#     ctrans.s_old = sim.s
#     ctrans.c = c + A * c + Q - Δs .* c
# end


text:
using SymbolicDiff
using SparseMatrix
using SparseArrays
using LinearAlgebra
using Test

#include("test_eval.jl")
#include("test_deriv.jl")
#include("test_vec.jl")
include("test_mat.jl")

#include("test_macro_promotion.jl")


text:
using Base.LinAlg: axpy!, scale!
using Knet, MPI
@useifgpu CUDArt
include("Problem.jl")

abstract AbsParticle

type Particle <: AbsParticle
    x::Array{Float64,1}
    v::Array{Float64,1}
    pbest::Array{Float64, 1}
    pbest_fitness::Float64

    Particle(x::Array{Float64, 1}) = new(x, similar(x), similar(x), 1e15)
    Particle(x::Array{Float64, 1}, v::Array{Float64, 1}) = new(x, v, similar(x), 1e12)
end

type NNParticle <: AbsParticle
    x::Knet.Net
    v::Knet.Net
    pbest::Knet.Net
    pbest_fitness::Float64

    NNParticle(x::Knet.Net) = new(x, Knet.Net(similar(x.reg)), Knet.Net(similar(x.reg)), 1e15)
    NNParticle(x::Knet.Net, v::Knet.Net) = new(x, v, deepcopy(x), 1e12)
    NNParticle(x::Knet.Net, v::Knet.Net, pbest::Knet.Net) = new(x, v, pbest, 1e12)
end

fitness(p::Particle, problem::RegularProblem) = problem.f(p.x)
fitness(p::NNParticle, problem::NNProblem) = fitness(p, problem, problem.data)
fitness(p::NNParticle, problem::NNProblem, data::Single) = problem.loss(forw(p.x, problem.data.x), problem.data.y)

function fitness(p::NNParticle, problem::NNProblem, data::MiniBatch)
    sumloss = numloss = 0
    for (x,ygold) in data.trn
        ypred = forw(p.x, x)
        sumloss += problem.loss(ypred, ygold)
        numloss += 1
    end
    return sumloss / numloss
end

function uniform(dim, l, h)
    arr = rand(dim)
    arr = arr * (h - l) + l
end

function uniform(dim, l::Array{Float64, 1}, h::Array{Float64, 1})
    arr = rand(dim)
    for i=1:length(arr); arr[i] = arr[i] * (h[i] - l[i]) + l[i] end
    return arr
end

function init_particles(num, problem::RegularProblem)
    ps = AbsParticle[]
    for p=1:num
        push!(ps, Particle(uniform(problem.dim, problem.l, problem.h), uniform(problem.dim, -1.0 , 1.0)))
    end
    gbest = Particle(uniform(problem.dim, problem.l, problem.h))
    return ps, gbest
end

init_particles(num, problem::NNProblem) = init_particles(num, problem, problem.data)

function init_particles(num, problem::NNProblem, Data::Single)
    ps = AbsParticle[]
    for p=1:num
        net1 = compile(problem.f);
        forw(net1, problem.data.x[:, 1])
        net2 = compile(problem.f);
        forw(net2, problem.data.x[:, 1])
        net3 = compile(problem.f);
        forw(net3, problem.data.x[:, 1])
        push!(ps, NNParticle(net1, net2, net3))
    end

    net = compile(problem.f);
    forw(net, problem.data.x[:,1])

    gbest = NNParticle(net)
    return ps, gbest
end

function init_particles(num, problem::NNProblem, Data::MiniBatch)
    ps = AbsParticle[]
    for p=1:num
        net1 = compile(problem.f);
        forw(net1, problem.data.trn[1][1])
        net2 = compile(problem.f);
        forw(net2, problem.data.trn[1][1])
        net3 = compile(problem.f);
        forw(net3, problem.data.trn[1][1])
        push!(ps, NNParticle(net1, net2, net3))
    end

    net = compile(problem.f);
    forw(net, problem.data.trn[1][1])
    gbest = NNParticle(net)
    return ps, gbest
end

function update!(px, pv, pbest, gbest, w, l::Float64, h::Float64, vmin, vmax, r1, r2)
    scale!(pv, w)
    axpy!(1, r1 * (pbest - px), pv)
    axpy!(1, r2 * (gbest - px), pv)
    for i=1:length(pv); pv[i] = pv[i] < vmin ? vmin : pv[i] > vmax ? vmax : pv[i] end
    axpy!(1, pv, px)
    for i=1:length(px); px[i] = px[i] < l ? l : px[i] > h ? h : px[i] end
end

function update!(px, pv, pbest, gbest, w, l::Array{Float64, 1}, h::Array{Float64, 1}, vmin, vmax, r1, r2)
    scale!(pv, w)
    axpy!(1, r1 * (pbest - px), pv)
    axpy!(1, r2 * (gbest - px), pv)
    for i=1:length(pv); pv[i] = pv[i] < vmin ? vmin : pv[i] > vmax ? vmax : pv[i] end
    axpy!(1, pv, px)
    for i=1:length(px); px[i] = px[i] < l[i] ? l[i] : px[i] > h[i] ? h[i] : px[i] end
end

const libpso = Libdl.find_library(["libpso"], [Pkg.dir("ParticleSwarmOptimization/src")])

@gpu update!(px::CudaArray{Float32}, pv::CudaArray{Float32}, pbest::CudaArray{Float32}, gbest::CudaArray{Float32}, w, l, h, vmin, vmax, r1, r2)=ccall((:update32,libpso),Void,(Cint,Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble),length(px),px, pv, pbest, gbest, w, l, h, vmin, vmax, r1, r2); gpusync();

@gpu update!(px::CudaArray{Float64}, pv::CudaArray{Float64}, pbest::CudaArray{Float64}, gbest::CudaArray{Float64}, w, l, h, vmin, vmax, r1, r2)=ccall((:update64,libpso),Void,(Cint,Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble),length(px),px, pv, pbest, gbest, w, l, h, vmin, vmax, r1, r2); gpusync();

@gpu update!(px::CudaArray{Float32}, pv::CudaArray{Float32}, pbest::CudaArray{Float32}, gbest::CudaArray{Float32}, w, l::CudaArray{Float64}, h::CudaArray{Float64}, vmin, vmax, r1, r2)=ccall((:update_32,libpso),Void,(Cint,Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Ptr{Float32}, Cdouble, Ptr{Float64}, Ptr{Float64}, Cdouble, Cdouble, Cdouble, Cdouble),length(px),px, pv, pbest, gbest, w, l, h, vmin, vmax, r1, r2); gpusync();

@gpu update!(px::CudaArray{Float64}, pv::CudaArray{Float64}, pbest::CudaArray{Float64}, gbest::CudaArray{Float64}, w, l::CudaArray{Float64}, h::CudaArray{Float64}, vmin, vmax, r1, r2)=ccall((:update_64,libpso),Void,(Cint,Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Cdouble, Ptr{Float64}, Ptr{Float64}, Cdouble, Cdouble, Cdouble, Cdouble),length(px),px, pv, pbest, gbest, w, l, h, vmin, vmax, r1, r2); gpusync();

function update!(problem::RegularProblem, p::Particle, gbest::Particle, phi1, phi2, w, vmin, vmax)
    l = problem.l
    h = problem.h
    update!(p.x, p.v, p.pbest, gbest.x, w, l, h, vmin, vmax, rand() * phi1, rand() * phi2)
end

function update!(problem::NNProblem, p::NNParticle, gbest::NNParticle, phi1, phi2, w, vmin, vmax)
    l = problem.l
    h = problem.h
    for i=1:length(p.v.reg)
        if isa(p.v.reg[i].op, Knet.Par)
            update!(p.x.reg[i].out, p.v.reg[i].out, p.pbest.reg[i].out, gbest.x.reg[i].out, w, l, h, vmin, vmax, rand() * phi1, rand() * phi2)
        end
    end
end

function update!(problem::Problem, ps::Array{AbsParticle, 1}, gbest::AbsParticle, phi1, phi2, w, vmin, vmax)
    for p in ps
        update!(problem, p, gbest, phi1, phi2, w, vmin, vmax)
    end
end

Bcast!(gbest::Particle, root::Int, comm::MPI.Comm) = MPI.Bcast!(gbest.x, length(gbest.x), root, comm)

addbuffer(buffer, arr) = (length(buffer) == 0) ? vec(arr) : vcat(buffer, vec(arr))
@gpu addbuffer(buffer, arr::CudaArray) = (length(buffer) == 0) ? vec(to_host(arr)) : vcat(buffer, vec(to_host(arr)))

function Bcast!(gbest::NNParticle, root::Int, comm::MPI.Comm)
    buffer = Any[]
    for i=1:length(gbest.x.reg)
        if isa(gbest.x.reg[i].op, Knet.Par)
            buffer = addbuffer(buffer, gbest.x.reg[i].out)
        end
    end
    #(MPI.Comm_rank(MPI.COMM_WORLD) == 0) && println("Root: $root Buffer1: $(buffer[1:5])")
    MPI.Bcast!(buffer, length(buffer), root, comm)
    #(MPI.Comm_rank(MPI.COMM_WORLD) == 0) && println("Root: $root Buffer2: $(buffer[1:5])")
    start = 1
    for i=1:length(gbest.x.reg)
        if isa(gbest.x.reg[i].op, Knet.Par)
            rs = size(gbest.x.reg[i].out)
            l = prod(rs)
            copy!(gbest.x.reg[i].out, reshape(buffer[start:(start+l-1)], rs))
            start += l
        end
    end
end

copypbest!(p::Particle) = copy!(p.pbest, p.x)

function copypbest!(p::NNParticle)
    for i=1:length(p.x.reg)
        if isa(p.x.reg[i].op, Knet.Par)
            copy!(p.pbest.reg[i].out, p.x.reg[i].out)
        end
    end
end

copygbest!(gbest::Particle, p::Particle) = copy!(gbest.x, p.x)

function copygbest!(gbest::NNParticle, p::NNParticle)
    for i=1:length(p.x.reg)
        if isa(p.x.reg[i].op, Knet.Par)
            copy!(gbest.x.reg[i].out, p.x.reg[i].out)
        end
    end
end


text:
using Burgers
using Test

@testset "Burgers" begin
    include("data.jl")
end


text:
module C2_2
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    l1_dot.*(4.99826877401903e-11*l1 - 3.06236268282961e-11*l2 - 1.96418867892903e-11*l3 + 4.29485424042179e-14) - l2_dot.*(3.06236268282961e-11*l1 - 3.88828793078652e-11*l2 + 8.05356093237194e-12*l3 + 3.15224925981314e-14) - l3_dot.*(1.96418867892903e-11*l1 + 8.05356093237194e-12*l2 - 2.77725822370262e-11*l3 + 1.19248520857872e-14)
end
end


text:
push!(LOAD_PATH, joinpath(pwd(), "src"))
using JuML
using Test

traintest_df = DataFrame(joinpath("data", "airlinetraintest"); preload=true)
train_df = DataFrame(joinpath("data", "airlinetrain"); preload=true)
test_df = DataFrame(joinpath("data", "airlinetest"))

distance = traintest_df["Distance"]
deptime = traintest_df["DepTime"]
label = covariate(traintest_df["dep_delayed_15min"], level -> level == "Y" ? 1.0 : 0.0)
deptime = factor(traintest_df["DepTime"])
distance = factor(traintest_df["Distance"])

factors = [traintest_df.factors; [deptime, distance]]

trainsel = BoolVariate("trainsel", (1:1100000) .<= 1000000)
validsel = BoolVariate("validsel", (1:1100000) .> 1000000)

model1 = xgblogit(label, factors; trainselector = trainsel, validselector = validsel,  η = 1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 1, maxdepth = 1, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = true, slicelength = 0);
_, testauc1 = getauc(model1.pred, label, trainsel, validsel)
@test testauc1 ≈ 0.634631 atol = 0.0000001

model2 = xgblogit(label, factors; trainselector = trainsel, validselector = validsel, η = 0.1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 2, maxdepth = 2, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = true, slicelength = 0);
_, testauc2 = getauc(model2.pred, label, trainsel, validsel)
@test testauc2 ≈ 0.6749474 atol = 0.0000001

model3 = xgblogit(label, factors; trainselector = trainsel, validselector = validsel, η = 1, λ = 1.0, γ = 0.0, minchildweight = 1000.0, nrounds = 1, maxdepth = 6, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = false, slicelength = 0);
_, testauc3 = getauc(model3.pred, label, trainsel, validsel)
@test testauc3 ≈ 0.7002925 atol = 0.0000001

model4 = xgblogit(label, factors; trainselector = trainsel, validselector = validsel, η = 1, λ = 1.0, γ = 500.0, minchildweight = 1.0, nrounds = 1, maxdepth = 4, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = false, slicelength = 0);
_, testauc4 = getauc(model4.pred, label, trainsel, validsel)
@test testauc4 ≈ 0.6888606 atol = 0.0000001

model5 = xgblogit(label, factors; trainselector = trainsel, validselector = validsel, η = 0.1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 10, maxdepth = 10, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = true, slicelength = 0);
_, testauc5 = getauc(model5.pred, label, trainsel, validsel)
@test testauc5 ≈ 0.7255029 atol = 0.0002

@time model6 = xgblogit(label, factors; trainselector = trainsel, validselector = validsel, η = 0.1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 1, maxdepth = 5, ordstumps = false, pruning = false, caching = false, usefloat64 = false, singlethread = true, slicelength = 0);
print(model6.pred[1:20])
trainauc6, testauc6 = getauc(model6.pred, label, trainsel, validsel)

distance = train_df["Distance"]
deptime = train_df["DepTime"]
label = covariate(train_df["dep_delayed_15min"], level -> level == "Y" ? 1.0 : 0.0)
deptime = factor(train_df["DepTime"])
distance = factor(train_df["Distance"])

factors = [train_df.factors; [deptime, distance]]

@time model7 = xgblogit(label, factors; η = 0.1, λ = 1.0, γ = 0.0, minchildweight = 1.0, nrounds = 1, maxdepth = 5, ordstumps = false, pruning = false, caching = true, usefloat64 = false, singlethread = true, slicelength = 0);
auc = getauc(model7.pred, label)
model7.pred[1:5]
p = predict(model7, test_df)
@test model6.pred[1000001:end] ≈ p atol = 1e-15

# XGBoost R script to compare:
# Data:
# wget https://s3.amazonaws.com/benchm-ml--main/train-1m.csv
# wget https://s3.amazonaws.com/benchm-ml--main/train-10m.csv
# wget https://s3.amazonaws.com/benchm-ml--main/test.csv

# suppressMessages({
# library(data.table)
# library(ROCR)
# library(xgboost)
# library(MLmetrics)
# library(Matrix)
# })

# d_train <- fread("C:\\Users\\adamm\\Documents\\julia\\airlinetrain1m.csv", showProgress=FALSE, stringsAsFactors=TRUE)
# d_test <- fread("C:\\Users\\adamm\\Documents\\julia\\airlinetest.csv", showProgress=FALSE, stringsAsFactors=TRUE)
# X_train_test <- sparse.model.matrix(dep_delayed_15min ~ .-1, data = rbind(d_train, d_test))
# n1 <- nrow(d_train)
# n2 <- nrow(d_test)
# X_train <- X_train_test[1:n1,]
# X_test <- X_train_test[(n1+1):(n1+n2),]
# dxgb_train <- xgb.DMatrix(data = X_train, label = ifelse(d_train$dep_delayed_15min=='Y',1,0))
# md <- xgb.train(data = dxgb_train, objective = "binary:logistic", nround = 1, max_depth = 1, eta = 1.0, tree_method='exact')
# # md <- xgb.train(data = dxgb_train, objective = "binary:logistic", nround = 10, max_depth = 10, eta = 0.1, tree_method='hist', grow_policy='lossguide', max_leaves=256, max_bin=5000)
# phat <- predict(md, newdata = X_test)
# testlabel = ifelse(d_test$dep_delayed_15min=='Y',1,0)
# AUC(phat, testlabel)

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
using Test
import Thermodynamics
const TD = Thermodynamics

import UnPack
import BenchmarkTools

import RootSolvers
const RS = RootSolvers

import CLIMAParameters
const CP = CLIMAParameters

struct EarthParameterSet <: CP.AbstractEarthParameterSet end
const param_set = EarthParameterSet()

ArrayType = Array{Float64}
profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, ArrayType)
UnPack.@unpack e_int, T, ρ, p, θ_liq_ice, q_tot = profiles

function thermo_state_ρeq()
    ts = TD.PhaseEquil_ρeq.(param_set, ρ, e_int, q_tot)
end

function thermo_state_pθq()
    ts = TD.PhaseEquil_pθq.(param_set, p, θ_liq_ice, q_tot)
end

function thermo_state_pTq()
    ts = TD.PhaseEquil_pTq.(param_set, p, T, q_tot)
end


text:
struct Figure6
    θ1::Vector{Float64}
    θ2::Vector{Float64}
    θ3::Vector{Float64}
    function Figure6(
            θ1 = [0.7, 0.8, 3.3, 5, π / 2, 4, 2.7, 0.55, 0.26],
            θ2 = [0.7, 1.1, 3.3, 5, π / 2, 4, 2.7, 0.55, 0],
            θ3 = [0.7, 1.0, 1, 5, π / 2, 4, 2.7, 0.55, 0.26]
        )
        new(θ1,θ2,θ3)
    end
end

@recipe function f(f::Figure6)
    layout := (1,3)
    for (i,θ) in enumerate((f.θ1,f.θ2,f.θ3))
        @series begin
            ylims --> (0.4,1.9)
            yticks --> (0.6:0.3:1.5)
            clims --> (0,1.25)
            levels --> 50
            seriescolor := :thermal
            subplot := i
            FrequencyDirectionSpectra.FDJONSWAP(θ),true
        end
    end
end

##
struct ColorFig6
    fig6::Figure6
end
@recipe function f(fc::ColorFig6)
    legend := true
    layout := 1
    ylims --> (0.4,1.9)
    yticks --> (0.6:0.3:1.5)
    clims --> (0,1.25)
    levels --> 50
    seriescolor := :thermal
    FrequencyDirectionSpectra.FDJONSWAP(fc.fig6.θ1),true
end


text:
"""
    AbstractLattice

Abstract type for general lattices.
For a more concrete definition please refer the following material:

- Lattice (group): https://en.wikipedia.org/wiki/Lattice_(group)
- Lattice Model: https://en.wikipedia.org/wiki/Lattice_model_(physics)
- Lattice Graph: https://en.wikipedia.org/wiki/Lattice_graph
"""
abstract type AbstractLattice end

unitcellsize(l::AbstractLattice) = 1
hasunitcell(::AbstractLattice) = Val{false}()

_ncoords(::Val{true}, l::AbstractLattice) = 1 + ndims(l)
_ncoords(::Val{false}, l::AbstractLattice) = ndims(l)
ncoordinates(l::AbstractLattice) = _ncoords(hasunitcell(l), l)

function size(::Val{true}, l::AbstractLattice, dim::Int)
    if dim > ncoordinates(l)
        return 1  # copying behaviour of `size` for multidimensional arrays
    elseif dim == ncoordinates(l)
        return unitcellsize(l)
    else
        return l.dims[dim]
    end
end
size(::Val{false}, l::AbstractLattice, dim::Int) = (dim > ndims(l)) ? 1 : l.dims[dim]
size(l::AbstractLattice, dim::Int) = size(hasunitcell(l), l, dim)

size(::Val{true}, l::AbstractLattice) = (l.dims..., unitcellsize(l))
size(::Val{false}, l::AbstractLattice) = l.dims
size(l::AbstractLattice) = size(hasunitcell(l), l)

axes(l::AbstractLattice, dim::Int) = Base.OneTo(size(l, dim))
axes(l::AbstractLattice) = map(d -> Base.OneTo(d), size(l))

strides(l::AbstractLattice) = (1, cumprod(size(l)[1:end-1])...)
nsites(l::AbstractLattice) = prod(size(l))

struct WeightedLattice{L <: AbstractLattice, W} <: AbstractLattice
    lattice::L
    weights::W
end

struct CoordinateLattice{L <: AbstractLattice, Coordinate} <: AbstractLattice
    lattice::L
    coordinates::Vector{Coordinate}
end

struct GraphLattice{Graph} <: AbstractLattice
    graph::Graph
end

# struct HasCoordinate end

# function coordinate(::HoneyComb, id)
# end

# function coordinate(lattice, id)
# end

# function coordinate(::HasCoordinate, lattice, id)
# end

# for (i, j) in edges(lattice; bond=k)
# end

# # <=>
# for i in sites(lattice)
#     for j in neighbors(lattice, i, k)
#         if i < j
#             value = get_term(H, i, j)
#         end
#     end
# end

# struct SitesIt{L}
#     lattice::L
# end

# struct BondsIt{L}
#     lattice::L
#     bond::Int
# end

# function sites end
# function neighbors end
# function bonds end

# Base.iterate(it::BondsIt, st) = iterate_lattice(it, it.lattice, st)
# function iterate_lattice(::BondsIt, lattice, st) end


text:
using Test
using CellMLToolkit
using OrdinaryDiffEq
using ModelingToolkit

@testset "CellMLToolkit.jl" begin
    @testset "beeler.jl" begin include("beeler.jl") end
    @testset "noble_1962.jl" begin include("noble_1962.jl") end
end


text:
function FiniteDifferenceAtParameters(P,ϵs)
    y = CF(P)
    ωᵢ = P.ω
    Almᵢ = P.Alm

    P.ω = ωᵢ + ϵs+im*0
    P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), Almᵢ,P.lmax)
    yʼr = (CF(P) - y)/ϵs

    P.ω = ωᵢ + im*ϵs;
    P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), Almᵢ,P.lmax)
    yʼi = (CF(P) - y)/ϵs

    P.ω = ωᵢ
    P.Alm = Almᵢ
    return yʼr,yʼi,y
end

function Newton!(P,ϵs)
    ω₀ = P.ω
    Alm₀ = P.Alm
    function f!(F,x)
        ωᵢ = P.ω
        Almᵢ = P.Alm
        P.ω = x[1]+im*x[2]
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
        val = CF(P)
        P.ω = ωᵢ
        P.Alm = Almᵢ
        F[1] = real(val)
        F[2] = imag(val)
    end

    function j!(J,x)
        ω = x[1]+im*x[2]
        P.ω = ω
        ydr,ydi,y = FiniteDifferenceAtParameters(P,ϵs)
        J[1,1] = real(ydr)
        J[1,2] = real(ydi)
        J[2,1] = imag(ydr)
        J[2,2] = imag(ydi)
    end

    function fj!(F,J,x)
        ωᵢ = P.ω
        Almᵢ = P.Alm
        P.ω = x[1]+im*x[2]
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), Almᵢ,P.lmax)
        val = CF(P)
        F[1] = real(val)
        F[2] = imag(val)

        P.ω = x[1]+im*x[2] + ϵs+im*0
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), Almᵢ,P.lmax)
        yʼr = (CF(P) - val)/ϵs

        P.ω = x[1]+im*x[2] + im*ϵs
        P.Alm,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), Almᵢ,P.lmax)
        yʼi = (CF(P) - val)/ϵs

        P.ω = ωᵢ
        P.Alm = Almᵢ

        J[1,1] = real(yʼr )
        J[1,2] = real(yʼi)
        J[2,1] = imag(yʼr )
        J[2,2] = imag(yʼi)
    end

    x0 = [real(ω₀),imag(ω₀)]
    df = OnceDifferentiable(f!, j!, fj!,x0,similar(x0))
    sol = nlsolve(df, x0, ftol = 1e-14)

    ωnew = sol.zero
    P.ω = ωnew[1] + im*ωnew[2]
    Almnew,_,_ = ComputeAₗₘ(P.s,P.m,P.a*(P.ω), P.Alm,P.lmax)
    P.Alm = Almnew
    P
end


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

name = "Exiv2"
version = v"0.27.4"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/Exiv2/exiv2/archive/refs/tags/v$(version).tar.gz", "9fb2752c92f63c9853e0bef9768f21138eeac046280f40ded5f37d06a34880d9"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""

cd $WORKSPACE/srcdir/exiv2-*

if [[ "${target}" == x86_64-linux-musl ]]; then
    # Delete libexpat to prevent it from being picked up by mistake
    rm /usr/lib/libexpat.so*

elif [[ "${target}" == i686-linux-musl ]]; then

    #otherwise, patch fails with different line endings message
    dos2unix cmake/compilerFlags.cmake

    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/i686-musl-comment-stack-protector-strong.patch

fi

mkdir build
cd build/

cmake .. \
-DCMAKE_INSTALL_PREFIX=$prefix \
-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
-DCMAKE_BUILD_TYPE=Release \
-DEXIV2_BUILD_SAMPLES=OFF \
-DIconv_INCLUDE_DIR=${includedir}

make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency("Expat_jll"; compat="2.2.10")
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using GRcoder
using Test

@testset "GRcoder.jl" begin
    # Write your tests here.
end


text:
# This is pretty naive and probably needs to be thought about a bit more.
# Right now this is sufficient to reproduce uniform fields for isolated systems.

@doc raw"""
Magnetic term ``A⋅(-i∇)``. It is assumed (but not checked) that ``∇⋅A = 0``.
"""
struct Magnetic
    Afunction::Function  # A(x,y,z) returns [Ax,Ay,Az]
end
function (M::Magnetic)(basis)
    TermMagnetic(basis, M.Afunction)
end

struct TermMagnetic <: Term
    basis::PlaneWaveBasis
    Apotential::AbstractArray  # Apotential[α] is an array of size fft_size for α=1:3
end
function TermMagnetic(basis::PlaneWaveBasis{T}, Afunction::Function) where T
    lattice = basis.model.lattice
    Apotential = [zeros(T, basis.fft_size) for α = 1:3]
    N1, N2, N3 = basis.fft_size
    rvecs = collect(r_vectors(basis))
    for i = 1:N1
        for j = 1:N2
            for k = 1:N3
                Apotential[1][i, j, k],
                Apotential[2][i, j, k],
                Apotential[3][i, j, k] = Afunction(lattice * rvecs[i, j, k])
            end
        end
    end
    TermMagnetic(basis, Apotential)
end

function ene_ops(term::TermMagnetic, ψ, occ; kwargs...)
    basis = term.basis
    T = eltype(basis)

    ops = [MagneticFieldOperator(basis, kpoint, term.Apotential)
           for (ik, kpoint) in enumerate(basis.kpoints)]
    ψ === nothing && return (E=T(Inf), ops=ops)

    E = zero(T)
    for (ik, k) in enumerate(basis.kpoints)
        for iband = 1:size(ψ[1], 2)
            ψnk = @views ψ[ik][:, iband]
            # TODO optimize this
            E += basis.kweights[ik] * occ[ik][iband] * real(dot(ψnk, ops[ik] * ψnk))
        end
    end

    (E=E, ops=ops)
end


text:
@testset "numerical gradients" begin
    f = a -> 2a .+ 1 .+ sin.(a)
    @test gradtest(f, 1.0)
    @test gradtest(f, 2.0)
    @test gradtest(f, randn(2, 2))
    b = 4
    f = a -> 2a .+ b .+ 1 .+ sin.(a.*b)
    @test gradtest(f, 1.0)
    @test gradtest(f, 2.0)
    @test gradtest(f, randn(2, 2))

    m1 = Dense(2, 2, relu) |> f64
    m2 = Dense(2, 2, relu) |> f64
    x = randn(2, 10)
    @test gradtest(() -> m1(x), Flux.params(m1))
    @test gradtest(() -> m2(x), Flux.params(m2))
end

@testset "gradf" begin
    b = 4
    f = gradf(a -> 2a .+ b .+ 1 .+ sin.(a.*b), 1)
    # outputs are not constant
    @test f(1) ≠ f(2)
    @test f(2) ≠ f(3)
    # gradients do differ
    @test gradient(f, 1) ≠ gradient(f, 2)
    @test gradient(f, 2) ≠ gradient(f, 3)
    # 2nd derivative
    df = gradf(f, 1)
    @test df(1) ≠ df(2)
    @test df(2) ≠ df(3)
    @test gradient(df, 1) ≠ gradient(df, 2)
    @test gradient(df, 2) ≠ gradient(df, 3)

    B = randn(2, 2)
    A1 = randn(2, 2)
    A2 = randn(2, 2)
    f = gradf(A -> 2A .+ B .+ 1 .+ sin.(A.*B), A1)
    @test f(A1) ≠ f(A2)
    @test gradient(f, A1) ≠ gradient(f, A2)
    df = gradf(f, A1)
    @test df(A1) ≠ df(A2)
    @test gradient(df, A1) ≠ gradient(df, A2)

    m1 = Dense(2, 2, relu) |> f64
    m2 = Dense(2, 2, relu) |> f64
    x = randn(2, 10)
    f1 = gradf(() -> m1(x), Flux.params(m1))
    f2 = gradf(() -> m2(x), Flux.params(m2))
    @test f1() ≠ f2()
    @test gradient(f1, Flux.params(m1)) ≠ gradient(f2, Flux.params(m2))
    @test gradient(f1, Flux.params(m1)) ≠ gradient(f1, Flux.params(m2))
    @test gradient(f2, Flux.params(m2)) ≠ gradient(f2, Flux.params(m1))
    df1 = gradf(f1, Flux.params(m1))
    df2 = gradf(f2, Flux.params(m2))
    @test df1() ≠ df2()
end

@testset "not implemented" begin
    f(x, y) = x + y
    function ChainRulesCore.rrule(::typeof(f), x, y)
        x+y, Δ -> (NoTangent(), Δ * 1, @not_implemented("Not implemented"))
    end
    @test_throws NotImplementedException gradtest(f, 1.0, 1.0)
    @test_throws NotImplementedException gradtest(y -> f(1, y), 1.0)
    @test gradtest(x -> f(x, 1), 1.0)
end


text:
module mznmem
using Glob

type RetString
    length::Csize_t
    content::Ptr{UInt8}
end

function readalldzn()
    ret = ""
    for f in glob("*.dzn")
        tmp = strip(readall(f))
        # Replace EOF with ';' when appending .dzn files together
        tmp[end] == ';' ? ret = string(ret, tmp) : ret = string(ret, tmp, ";\n")
    end
    ret
end

"""Read the minizinc code from a file without any seperate data"""
solvefile(filename::AbstractString, numsolutions = 1, verbose = false) =
    solve(readall(filename), numsolutions, verbose)

"""Read the minizinc code from a file and pass the data"""
function solvefile(filename::AbstractString, data::AbstractString, numsolutions = 1,
        verbose = false)
    tmp = strip(readall(filename))
    # Replace EOF with ';' when appending files together.
    # FIXME: Fails if last char is part of a comment
    tmp[end] == ';' ? tmp = string(tmp, data) : tmp = string(tmp, ";\n", data)
    solve(tmp, numsolutions, verbose)
end

"""Read the minizinc code and data seperately"""
function solve(mzn::AbstractString, data::AbstractString, numsolutions = 1, verbose = false)
    mzn[end] == ';' ? mzn = string(mzn, data) : mzn = string(mzn, ";\n", data)
    solve(mzn, numsolutions, verbose)
end

"""Read the minizinc code and solve it"""
function solve(data::AbstractString,  numsolutions = 1, verbose = false)
    ret = ccall((:solve, "libmznmem.so"), RetString, (Ptr{UInt8}, Cint, Bool, Ptr{UInt8}), data,
        numsolutions, verbose, std_lib_dir)
    str = utf8(pointer_to_array(ret.content, ret.length))
    ccall((:release_retstring, "libmznmem.so"), Void, (RetString, ), ret)
    split(str, '\n')
end

function __init__()
    global std_lib_dir
    try
        std_lib_dir = ENV["MZN_STDLIB_DIR"]
        isfile(joinpath(std_lib_dir, "std/builtins.mzn")) || throw()
    catch
        mypath = "/usr/local"
        if isfile(joinpath(mypath, "share/minizinc/std/builtins.mzn"))
            std_lib_dir = joinpath(mypath,  "share/minizinc")
        elseif isfile(string(mypath, "../share/minizinc/std/builtins.mzn"))
            std_lib_dir = joinpath(mypath, "../share/minizinc")
        else
            throw("Error: unknown minizinc standard library directory. Set the MZN_STDLIB_DIR
                environment variable.")
        end
    end

end

end

max_stars_repo_name: vettujayan/minizinc-julia

max_issues_repo_name: vettujayan/minizinc-julia

max_forks_repo_name: vettujayan/minizinc-julia


text:
using NativeSVG

dr =
    Drawing(width = "5cm", height = "4cm") do
        desc() do
            str("Four separate rectangles")
        end
        rect(x = "0.5cm", y = "0.5cm", width = "2cm", height = "1cm")
        rect(x = "0.5cm", y = "2cm", width = "1cm", height = "1.5cm")
        rect(x = "3cm", y = "0.5cm", width = "1.5cm", height = "2cm")
        rect(x = "3.5cm", y = "3cm", width = "1cm", height = "0.5cm")

    #Show outline of viewport using 'rect' element
        rect(
            x = ".01cm",
            y = ".01cm",
            width = "4.98cm",
            height = "3.98cm",
            fill = "none",
            stroke = "blue",
            stroke_width = ".02cm"
        )
    end
display(dr)

dr =
    Drawing(width = "5cm", height = "5cm") do
        desc() do
            str("Two groups, each of two rectangles")
        end
        g(id = "group1", fill = "red") do
            rect(x = "1cm", y = "1cm", width = "1cm", height = "1cm")
            rect(x = "3cm", y = "1cm", width = "1cm", height = "1cm")
        end
        g(id = "group2", fill = "blue") do
            rect(x = "1cm", y = "3cm", width = "1cm", height = "1cm")
            rect(x = "3cm", y = "3cm", width = "1cm", height = "1cm")
        end

    # Show outline of viewport using 'rect' element
        rect(
            x = ".01cm",
            y = ".01cm",
            width = "4.98cm",
            height = "4.98cm",
            fill = "none",
            stroke = "blue",
            stroke_width = ".02cm"
        )
    end
display(dr)

dr =
    Drawing(
        xmlns!xlink = "http://www.w3.org/1999/xlink",
        width = "200",
        height = "100",
        viewBox = "0 0 200 100"
    ) do
        title() do
            str("Style inheritance and the use element")
        end
        desc() do
            str("""Two circles, one of which is a re-styled clone of the other.
                   This file demonstrates one of the cases where
                   the shadow-DOM style matching rules in SVG 2
                   have a different effect than the SVG 1.1 style cloning rules.
                   The original circle on the left
                   should have blue fill
                   and green stroke.
                   In a conforming SVG 1.1 user agent,
                   the re-used circle on the right
                   should have orange fill and green stroke.
                   In a conforming SVG 2 user agent,
                   the re-used circle should have orange fill and purple stroke.
                   In all cases,
                   the stroke should be partially transparent
                   and 20 units wide,
                   relative to a total circle diameter of 100 units.""")
        end
        style(type = "text/css") do
            str("""circle          { stroke-opacity: 0.7; }
                   .special circle { stroke: green; }
                   use             { stroke: purple;
                                     fill: orange; }""")
        end
        g(class = "special", style = "fill: blue") do
            circle(
                id = "c",
                cy = "50",
                cx = "50",
                r = "40",
                stroke_width = "20"
            )
        end
        use(xlink!href = "#c", x = "100")
    end

max_stars_repo_name: BenLauwens/NativeSVG

max_issues_repo_name: BenLauwens/NativeSVG

max_forks_repo_name: BenLauwens/NativeSVG


text:
import DataFrames
import MLBase
import StatsBase

"""
"""
function singlelabelbinaryytrue(
        labels::AbstractVector,
        positive_class::AbstractString;
        inttype::Type=Int,
        )
    if !(inttype<:Integer)
        error("!(inttype<:Integer)")
    end
    result = inttype.(labels .== positive_class)
    return result
end

"""
"""
function singlelabelbinaryyscore(
        single_labelprobabilities::AbstractDict,
        positive_class::AbstractString;
        float_type::Type{<:AbstractFloat} = Float64,
        )
    result = float_type.(single_labelprobabilities[positive_class])
    return result
end

"""
"""
function singlelabelbinaryclassificationmetrics_tunableparam(
        kwargsassoc::AbstractDict,
        )
    tunableparams = [
        :threshold,
        :sensitivity,
        :specificity,
        :maximize,
        ]
    maximizableparams = [
        :f1score,
        :cohen_kappa,
        ]
    kwargshastunableparam = [
        haskey(kwargsassoc, x) for x in tunableparams
        ]
    if sum(kwargshastunableparam) != 1
        msg = "you must specify one (and only one) of the following: " *
            join(tunableparams, ", ")
        error(msg)
    end
    if length(tunableparams[kwargshastunableparam]) != 1
        error("oh boy you definitely should never see this error message")
    end
    selectedtunableparam =
        tunableparams[kwargshastunableparam][1]
    if selectedtunableparam == :maximize
        selectedparamtomax = kwargsassoc[:maximize]
        if !in(selectedparamtomax, maximizableparams)
            msg = "Cannot max $(selectedparamtomax). Select one " *
                "of the following: " * join(maximizableparams, ", ")
            error(msg)
        end
    else
        selectedparamtomax = :notapplicable
    end
    #
    metricdisplaynames = Dict()
    metricdisplaynames[:AUPRC] = string("AUPRC")
    metricdisplaynames[:AUROCC] = string("AUROCC")
    metricdisplaynames[:AveragePrecision] = string("Average precision")
    if selectedtunableparam == :threshold
        metricdisplaynames[:threshold] = string("[fix] * Threshold")
    else
        metricdisplaynames[:threshold] = string("* Threshold")
    end
    metricdisplaynames[:accuracy] = string("* Accuracy")
    if selectedtunableparam == :maximize && selectedparamtomax ==
            :cohen_kappa
        metricdisplaynames[:cohen_kappa] =
            string("[max] * Cohen's Kappa statistic")
    else
        metricdisplaynames[:cohen_kappa] =
            string("* Cohen's Kappa statistic")
    end
    if selectedtunableparam == :maximize && selectedparamtomax ==
            :f1score
        metricdisplaynames[:f1score] = string("[max] * F1 score")
    else
        metricdisplaynames[:f1score] = string("* F1 Score")
    end
    metricdisplaynames[:precision] =
        string("* Precision (positive predictive value)")
    metricdisplaynames[:negative_predictive_value] =
        string("* Negative predictive value")
    metricdisplaynames[:recall] =
        string("* Recall (sensitivity, true positive rate)")
    if selectedtunableparam == :sensitivity
        metricdisplaynames[:sensitivity] =
            string("[fix] * Sensitivity (recall, true positive rate)")
    else
        metricdisplaynames[:sensitivity] =
            string("* Sensitivity (recall, true positive rate)")
    end
    if selectedtunableparam == :specificity
        metricdisplaynames[:specificity] =
            string("[fix] * Specificity (true negative rate)")
    else
        metricdisplaynames[:specificity] =
            string("* Specificity (true negative rate)")
    end
    metricdisplaynames = fix_type(metricdisplaynames)
    return selectedtunableparam, selectedparamtomax, metricdisplaynames
end

"""
"""
function singlelabelbinaryclassificationmetrics_resultdict(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        kwargs...
        )
    #
    kwargsdict = Dict(kwargs)
    kwargsdict = fix_type(kwargsdict)
    selectedtunableparam, selectedparamtomax, metricdisplaynames =
        singlelabelbinaryclassificationmetrics_tunableparam(kwargsdict)
    #
    predictedprobabilitiesalllabels = predict_proba(estimator, features_df)
    yscore = Float64.(
        singlelabelbinaryyscore(
            predictedprobabilitiesalllabels[single_label_name],
            positive_class,
            )
        )
    ytrue = Int.(
        singlelabelbinaryytrue(
            labels_df[single_label_name],
            positive_class,
            )
        )
    results = Dict()
    results[:ytrue] = ytrue
    results[:yscore] = yscore
    results[:AUROCC] = aurocc(ytrue, yscore)
    results[:AUPRC] = auprc(ytrue, yscore)
    results[:AveragePrecision] = averageprecisionscore(ytrue, yscore)
    if selectedtunableparam == :threshold
        additionalthreshold = kwargsdict[:threshold]
    else
        additionalthreshold = 0.5
    end
    allrocnums, allthresholds = getallrocnums(
        ytrue,
        yscore;
        additionalthreshold = additionalthreshold,
        )
    if selectedtunableparam == :threshold
        selectedthreshold = kwargsdict[:threshold]
        bestindex = argmin(abs.(allthresholds .- selectedthreshold))
    elseif selectedtunableparam == :sensitivity
        selectedsensitivity = kwargsdict[:sensitivity]
        allsensitivity = [sensitivity(x) for x in allrocnums]
        bestindex = argmin(abs.(allsensitivity .- selectedsensitivity))
    elseif selectedtunableparam == :specificity
        selectedspecificity = kwargsdict[:specificity]
        allspecificity = [specificity(x) for x in allrocnums]
        bestindex = argmin(abs.(allspecificity .- selectedspecificity))
    elseif selectedtunableparam == :maximize
        selectedparamtomax = kwargsdict[:maximize]
        if selectedparamtomax == :f1score
            allf1score = [fbetascore(x, 1) for x in allrocnums]
            bestindex = argmin(allf1score)
        elseif selectedparamtomax == :cohen_kappa
            allcohen_kappa = [cohen_kappa(x) for x in allrocnums]
            bestindex = argmin(allcohen_kappa)
        else
            error("this is an error that should never happen")
        end
    else
        error("this is another error that should never happen")
    end

    results[:allrocnums] = allrocnums
    results[:allthresholds] = allthresholds
    results[:bestindex] = bestindex
    bestrocnums = allrocnums[bestindex]
    bestthreshold = allthresholds[bestindex]
    results[:threshold] = bestthreshold
    results[:accuracy] = accuracy(bestrocnums)
    results[:sensitivity] = sensitivity(bestrocnums)
    results[:specificity] = specificity(bestrocnums)
    results[:precision] = precision(bestrocnums)
    results[:negative_predictive_value] =
        negative_predictive_value(bestrocnums)
    results[:recall] = recall(bestrocnums)
    results[:f1score] = f1score(bestrocnums)
    results[:cohen_kappa] = cohen_kappa(bestrocnums)
    results = fix_type(results)
    return results
end

"""
"""
function singlelabelbinaryclassificationmetrics(
        estimator::AbstractFittable,
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        kwargs...
        )
    vectorofestimators = AbstractFittable[estimator]
    result = singlelabelbinaryclassificationmetrics(
        vectorofestimators,
        features_df,
        labels_df,
        single_label_name,
        positive_class;
        kwargs...
        )
    return result
end

"""
"""
function singlelabelbinaryclassificationmetrics(
        vectorofestimators::AbstractVector{AbstractFittable},
        features_df::DataFrames.AbstractDataFrame,
        labels_df::DataFrames.AbstractDataFrame,
        single_label_name::Symbol,
        positive_class::AbstractString;
        kwargs...
        )
    kwargsdict = Dict(kwargs)
    kwargsdict = fix_type(kwargsdict)
    selectedtunableparam, selectedparamtomax, metricdisplaynames =
        singlelabelbinaryclassificationmetrics_tunableparam(kwargsdict)
    metricsforeachestimator = [
        singlelabelbinaryclassificationmetrics_resultdict(
            est,
            features_df,
            labels_df,
            single_label_name,
            positive_class;
            kwargs...
            )
            for est in vectorofestimators
        ]
    result = DataFrames.DataFrame()
    result[:metric] = [
        metricdisplaynames[:AUPRC],
        metricdisplaynames[:AUROCC],
        metricdisplaynames[:AveragePrecision],
        metricdisplaynames[:threshold],
        metricdisplaynames[:accuracy],
        metricdisplaynames[:cohen_kappa],
        metricdisplaynames[:f1score],
        metricdisplaynames[:precision],
        metricdisplaynames[:negative_predictive_value],
        metricdisplaynames[:recall],
        metricdisplaynames[:sensitivity],
        metricdisplaynames[:specificity],
        ]
    for i = 1:length(vectorofestimators)
        result[Symbol(vectorofestimators[i].name)] = [
            metricsforeachestimator[i][:AUPRC],
            metricsforeachestimator[i][:AUROCC],
            metricsforeachestimator[i][:AveragePrecision],
            metricsforeachestimator[i][:threshold],
            metricsforeachestimator[i][:accuracy],
            metricsforeachestimator[i][:cohen_kappa],
            metricsforeachestimator[i][:f1score],
            metricsforeachestimator[i][:precision],
            metricsforeachestimator[i][:negative_predictive_value],
            metricsforeachestimator[i][:recall],
            metricsforeachestimator[i][:sensitivity],
            metricsforeachestimator[i][:specificity],
            ]
    end
    return result
end


text:
@defcomp GDP begin
    # GDP: Gross domestic product $M
    # GRW: GDP growth rate %/year
    region            = Index()

    # Variables
    gdp               = Variable(index=[time, region], unit="\$M")
    cons_consumption  = Variable(index=[time, region], unit="\$million")
    cons_percap_consumption = Variable(index=[time, region], unit="\$/person")
    cons_percap_consumption_0 = Variable(index=[region], unit="\$/person")
    yagg_periodspan = Variable(index=[time], unit="year")

    # Parameters
    y_year_0          = Parameter(unit="year")
    y_year            = Parameter(index=[time], unit="year")
    grw_gdpgrowthrate = Parameter(index=[time, region], unit="%/year") # From p.32 of Hope 2009
    gdp_0             = Parameter(index=[region], unit="\$M") # GDP in y_year_0
    save_savingsrate  = Parameter(unit="%", default=15.00) # pp33 PAGE09 documentation, "savings rate".
    pop0_initpopulation = Parameter(index=[region], unit="million person")
    pop_population    = Parameter(index=[time,region], unit="million person")

    # Saturation, used in impacts
    isat0_initialimpactfxnsaturation = Parameter(unit="unitless", default=20.0) # pp34 PAGE09 documentation
    isatg_impactfxnsaturation = Variable(unit="unitless")

    ###############################################
    # Growth Effects - additional variables and parameters
    ###############################################
    # parameters and variables for growth effects
    gdp_leveleffect   = Variable(index=[time, region], unit="\$M")
    isat_ImpactinclSaturationandAdaptation = Parameter(index=[time,region], unit="\$")
    lgdp_gdploss =  Variable(index=[time, region], unit="\$M")
    ge_growtheffects = Parameter(unit = "none", default =  0.)
    geadrate_growtheffects_adaptationrate = Parameter(unit = "none", default = 0.)
    geadpt_growtheffects_adapted = Variable(index = [time], unit = "none")
    grwnet_realizedgdpgrowth = Variable(index=[time, region], unit = "%/year")
    # bound variables
    use_convergence = Parameter(unit="none", default=1.)
    cbabs_pcconsumptionbound = Parameter(unit="\$/person", default=740.65)
    cbabsn_pcconsumptionbound_neighbourhood = Variable(unit="\$/person")
    cbaux1_pcconsumptionbound_auxiliary1 = Variable(unit="none")
    cbaux2_pcconsumptionbound_auxiliary2 = Variable(unit="none")
    cbreg_regionsatbound = Variable(index=[time, region], unit="regions")
    cons_consumption_noconvergence = Variable(index=[time, region], unit="\$million")
    cons_percap_consumption_noconvergence = Variable(index=[time, region], unit="\$/person")
    # switch to overwrite the growth effects parameter with empirical distribution
    ge_empirical_distribution = Parameter(index=[draw], unit="none")
    ge_use_empiricaldistribution = Parameter(unit="none", default=0.)
    ge_seed_empiricaldistribution = Parameter(unit="none", default=1.)
    # switch to shut off persistence for certain regions
    ge_regionswitch = Parameter(index=[region], unit="none")
    ge_use_regionswitch = Parameter(unit="none", default = 0.)
    ###############################################

    function init(p, v, d)

        v.isatg_impactfxnsaturation = p.isat0_initialimpactfxnsaturation * (1 - p.save_savingsrate / 100)
        for rr in d.region
            v.cons_percap_consumption_0[rr] = (p.gdp_0[rr] / p.pop0_initpopulation[rr]) * (1 - p.save_savingsrate / 100)
        end
    end

    function run_timestep(p, v, d, t)

        # Analysis period ranges - required for abatemnt costs and equity weighting, from Hope (2006)
        if is_first(t)
            ylo_periodstart = p.y_year_0
        else
            ylo_periodstart = (p.y_year[t] + p.y_year[t - 1]) / 2
        end

        if t.t == length(p.y_year)
            yhi_periodend = p.y_year[t]
        else
            yhi_periodend = (p.y_year[t] + p.y_year[t + 1]) / 2
        end

        v.yagg_periodspan[t] = yhi_periodend - ylo_periodstart

        if is_first(t)
            # calculate the lower consumption bound and the consumption level which triggers the convergence system
            v.cbabsn_pcconsumptionbound_neighbourhood = p.cbabs_pcconsumptionbound * 1.5

            # define auxiliary parameters for the convergence system to make formulas easier to read
            v.cbaux1_pcconsumptionbound_auxiliary1 = 2 * (p.cbabs_pcconsumptionbound - v.cbabsn_pcconsumptionbound_neighbourhood)
            v.cbaux2_pcconsumptionbound_auxiliary2 = 4 / v.cbaux1_pcconsumptionbound_auxiliary1

            # if the switch is set to one, overwrite the growth effects parameter with an empirical distribution, using the seed parameter
            if p.ge_use_empiricaldistribution == 1.
                Random.seed!(trunc(Int, p.ge_seed_empiricaldistribution))
                p.ge_growtheffects = p.ge_empirical_distribution[Random.rand(1:10^6)]
            elseif p.ge_use_empiricaldistribution == 2. # set negative rho values to zero if switch set to two
                Random.seed!(trunc(Int, p.ge_seed_empiricaldistribution))
                p.ge_growtheffects = max(p.ge_empirical_distribution[Random.rand(1:10^6)], 0)
            end
        end

        # calculate the growth effects parameter net of the assumed adaptation to damage persistence
        v.geadpt_growtheffects_adapted[t] = p.ge_growtheffects * (1 - p.geadrate_growtheffects_adaptationrate)^(p.y_year[t] - p.y_year_0)

        for r in d.region
            # eq.28 in Hope 2002
            if is_first(t)
                v.gdp[t, r] = p.gdp_0[r] * (1 + (p.grw_gdpgrowthrate[t,r] / 100))^(p.y_year[t] - p.y_year_0)

                v.grwnet_realizedgdpgrowth[t,r] = p.grw_gdpgrowthrate[t,r]
                v.gdp_leveleffect[t,r] = v.gdp[t,r]

                v.cons_consumption[t, r] = v.gdp[t, r] * (1 - p.save_savingsrate / 100)
                v.cons_percap_consumption[t, r] = v.cons_consumption[t, r] / p.pop_population[t, r]
            else
                # if region switch is used, multiply the growth effect by the switch; otherwise, multiply by one
                v.grwnet_realizedgdpgrowth[t,r] = p.grw_gdpgrowthrate[t,r] - ifelse(p.ge_use_regionswitch == 1., p.ge_regionswitch[r], 1.) * v.geadpt_growtheffects_adapted[t] * p.isat_ImpactinclSaturationandAdaptation[t-1,r]
                v.gdp[t, r] = v.gdp[t-1, r] * (1 + (v.grwnet_realizedgdpgrowth[t,r]/100))^(p.y_year[t] - p.y_year[t-1])
                v.gdp_leveleffect[t,r] = v.gdp_leveleffect[t-1, r] *  (1 + (p.grw_gdpgrowthrate[t,r]/100))^(p.y_year[t] - p.y_year[t-1])

                v.cons_consumption[t, r] = v.gdp[t, r] * (1 - p.save_savingsrate / 100)
                v.cons_percap_consumption[t, r] = v.cons_consumption[t, r] / p.pop_population[t, r]

                # let boundary take effect if pc consumption is in the neighbourhood of the boundary
                if p.use_convergence == 1.
                    if v.cons_percap_consumption[t,r] >= v.cbabsn_pcconsumptionbound_neighbourhood
                        v.cbreg_regionsatbound[t,r] = 0.
                        v.cons_consumption_noconvergence[t,r] = v.cons_consumption[t,r]
                        v.cons_percap_consumption_noconvergence[t,r] = v.cons_consumption_noconvergence[t,r] / p.pop_population[t,r]
                    else
                        # calculate the consumption level if there was no convergence system
                        v.cons_consumption_noconvergence[t,r] = v.cons_consumption_noconvergence[t - 1, r] * (1 + (v.grwnet_realizedgdpgrowth[t,r] / 100))^(p.y_year[t] - p.y_year[t - 1])
                        v.cons_percap_consumption_noconvergence[t,r] = v.cons_consumption_noconvergence[t,r] / p.pop_population[t,r]

                        # send the pc cconsumption on a logistic path convergenging against the bound
                        v.cons_percap_consumption[t,r] = v.cbabsn_pcconsumptionbound_neighbourhood - 0.5 * v.cbaux1_pcconsumptionbound_auxiliary1 +
                                            v.cbaux1_pcconsumptionbound_auxiliary1 * exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence[t,r] - v.cbabsn_pcconsumptionbound_neighbourhood)) /
                                                                                (1 + exp(v.cbaux2_pcconsumptionbound_auxiliary2 *
                                                                                (v.cons_percap_consumption_noconvergence[t,r] - v.cbabsn_pcconsumptionbound_neighbourhood)))

                    # recalculate all variables accordingly
                        v.cons_consumption[t, r] = v.cons_percap_consumption[t,r] * p.pop_population[t,r]
                        v.gdp[t, r] = v.cons_consumption[t, r] / (1 - p.save_savingsrate / 100)
                        v.grwnet_realizedgdpgrowth[t,r] = 100 * ((v.gdp[t, r] / v.gdp[t - 1, r])^(1 / (p.y_year[t] - p.y_year[t - 1])) - 1)

                        v.cbreg_regionsatbound[t,r] = 1.

                    end
                else
                    if v.cons_percap_consumption[t,r] < p.cbabs_pcconsumptionbound
                        v.cons_percap_consumption[t,r] = p.cbabs_pcconsumptionbound

                        # recalculate all variables accordingly
                        v.cons_consumption[t, r] = v.cons_percap_consumption[t,r] * p.pop_population[t,r]
                        v.gdp[t, r] = v.cons_consumption[t, r] / (1 - p.save_savingsrate / 100)
                        v.grwnet_realizedgdpgrowth[t,r] = 100 * ((v.gdp[t, r] / v.gdp[t - 1, r])^(1 / (p.y_year[t] - p.y_year[t - 1])) - 1)

                        v.cbreg_regionsatbound[t,r] = 1.
                    end
                end
            end

            v.lgdp_gdploss[t,r] = v.gdp_leveleffect[t,r] - v.gdp[t,r]

        end
    end
end


text:
using DelimitedFiles, Statistics, LinearAlgebra

include("HRP.jl")

z = readdlm(ARGS[1], ',')

w = size(z, 1)
if length(ARGS) > 1
    w = min(w, parse(Int, ARGS[2]))
end

_, whrp = HRP.hrp(z[end-w+1:end, :])

open("weights.csv", "w") do io
    writedlm(io, whrp, ',')
end

println(whrp)

max_stars_repo_name: banachtech/hrp

max_issues_repo_name: banachtech/hrp

max_forks_repo_name: banachtech/hrp


text:
using LinearAlgebra, Logging, Random
using CUDAdrv
using CUDAnative
using CuArrays
using GPUifyLoops
using StaticArrays
using ArgParse

# {{{ constants
# note the order of the fields below is also assumed in the code.
const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate
const stateid = (ρ = _ρ, U = _U, V = _V, W = _W, E = _E)

const _nvgeo = 14
const _ξx, _ηx, _ζx, _ξy, _ηy, _ζy, _ξz, _ηz, _ζz, _MJ, _MJI,
_x, _y, _z = 1:_nvgeo
const vgeoid = (ξx = _ξx, ηx = _ηx, ζx = _ζx,
                ξy = _ξy, ηy = _ηy, ζy = _ζy,
                ξz = _ξz, ηz = _ηz, ζz = _ζz,
                MJ = _MJ, MJI = _MJI,
                x = _x,   y = _y,   z = _z)

Base.@irrational grav 9.81 BigFloat(9.81)
Base.@irrational R_d 287.0024093890231 BigFloat(287.0024093890231)
Base.@irrational cp_d 1004.5084328615809 BigFloat(1004.5084328615809)
Base.@irrational cv_d 717.5060234725578 BigFloat(717.5060234725578)
Base.@irrational gamma_d 1.4 BigFloat(1.4)
Base.@irrational gdm1 0.4 BigFloat(0.4)
# }}}

# {{{ volumerhs_v2!
function volumerhs_v2!(::Val{3}, ::Val{N}, ::Val{nmoist},
                       ::Val{ntrace}, rhs, Q, vgeo, gravity, D,
                       nelem) where {N, nmoist, ntrace}
    DFloat = eltype(Q)

    nvar = _nstate + nmoist + ntrace

    Nq = N + 1

    s_D = @shmem DFloat (Nq, Nq)
    s_F = @shmem DFloat (Nq, Nq, Nq, _nstate)
    s_G = @shmem DFloat (Nq, Nq, Nq, _nstate)
    s_H = @shmem DFloat (Nq, Nq, Nq, _nstate)
    l_ρinv = @shmem DFloat (Nq, Nq, Nq)

    @inbounds @loop for e in (1:nelem; blockIdx().x)
        @loop for k in (1:Nq; threadIdx().z)
            @loop for j in (1:Nq; threadIdx().y)
                @loop for i in (1:Nq; threadIdx().x)

                    if k == 1
                        s_D[i, j] = D[i, j]
                    end

                    # Load values will need into registers
                    MJ = vgeo[i, j, k, _MJ, e]
                    ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
                    ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
                    ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]
                    z = vgeo[i,j,k,_z,e]

                    U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
                    ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

                    P = gdm1*(E - (U^2 + V^2 + W^2)/(2*ρ) - ρ*gravity*z)

                    #          l_ρinv[i, j, k] = ρinv = 1 / ρ
                    l_ρinv[i, j, k] = ρinv = 1 / ρ
                    fluxρ_x = U
                    fluxU_x = ρinv * U * U + P
                    fluxV_x = ρinv * U * V
                    fluxW_x = ρinv * U * W
                    fluxE_x = ρinv * U * (E + P)

                    fluxρ_y = V
                    fluxU_y = ρinv * V * U
                    fluxV_y = ρinv * V * V + P
                    fluxW_y = ρinv * V * W
                    fluxE_y = ρinv * V * (E + P)

                    fluxρ_z = W
                    fluxU_z = ρinv * W * U
                    fluxV_z = ρinv * W * V
                    fluxW_z = ρinv * W * W + P
                    fluxE_z = ρinv * W * (E + P)

                    s_F[i, j, k, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
                    s_F[i, j, k, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
                    s_F[i, j, k, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
                    s_F[i, j, k, _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
                    s_F[i, j, k, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

                    s_G[i, j, k, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
                    s_G[i, j, k, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
                    s_G[i, j, k, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
                    s_G[i, j, k, _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
                    s_G[i, j, k, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

                    s_H[i, j, k, _ρ] = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
                    s_H[i, j, k, _U] = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
                    s_H[i, j, k, _V] = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
                    s_H[i, j, k, _W] = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
                    s_H[i, j, k, _E] = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)
                end
            end
        end

        @synchronize

        @loop for k in (1:Nq; threadIdx().z)
            @loop for j in (1:Nq; threadIdx().y)
                @loop for i in (1:Nq; threadIdx().x)
                    # TODO: Prefetch MJI and rhs

                    rhsU = rhsV = rhsW = rhsρ = rhsE = zero(DFloat)
                    MJI = vgeo[i, j, k, _MJI, e]

                    # buoyancy term
                    ρ = Q[i, j, k, _ρ, e]
                    rhsW -= ρ * gravity

                    # loop of ξ-grid lines
                    for n = 1:Nq
                        MJI_Dni = MJI * s_D[n, i]
                        MJI_Dnj = MJI * s_D[n, j]
                        MJI_Dnk = MJI * s_D[n, k]

                        rhsρ += MJI_Dni * s_F[n, j, k, _ρ]
                        rhsρ += MJI_Dnj * s_G[i, n, k, _ρ]
                        rhsρ += MJI_Dnk * s_H[i, j, n, _ρ]

                        rhsU += MJI_Dni * s_F[n, j, k, _U]
                        rhsU += MJI_Dnj * s_G[i, n, k, _U]
                        rhsU += MJI_Dnk * s_H[i, j, n, _U]

                        rhsV += MJI_Dni * s_F[n, j, k, _V]
                        rhsV += MJI_Dnj * s_G[i, n, k, _V]
                        rhsV += MJI_Dnk * s_H[i, j, n, _V]

                        rhsW += MJI_Dni * s_F[n, j, k, _W]
                        rhsW += MJI_Dnj * s_G[i, n, k, _W]
                        rhsW += MJI_Dnk * s_H[i, j, n, _W]

                        rhsE += MJI_Dni * s_F[n, j, k, _E]
                        rhsE += MJI_Dnj * s_G[i, n, k, _E]
                        rhsE += MJI_Dnk * s_H[i, j, n, _E]
                    end

                    rhs[i, j, k, _U, e] += rhsU
                    rhs[i, j, k, _V, e] += rhsV
                    rhs[i, j, k, _W, e] += rhsW
                    rhs[i, j, k, _ρ, e] += rhsρ
                    rhs[i, j, k, _E, e] += rhsE
                end
            end
        end

        # loop over moist variables
        # FIXME: Currently just passive advection
        # TODO: This should probably be unrolled by some factor
        rhsmoist = zero(eltype(rhs))
        for m = 1:nmoist
            s = _nstate + m

            @synchronize

            @loop for k in (1:Nq; threadIdx().z)
                @loop for j in (1:Nq; threadIdx().y)
                    @loop for i in (1:Nq; threadIdx().x)
                        MJ = vgeo[i, j, k, _MJ, e]
                        ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
                        ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
                        ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

                        Qmoist = Q[i, j, k, s, e]
                        U = Q[i, j, k, _U, e]
                        V = Q[i, j, k, _V, e]
                        W = Q[i, j, k, _W, e]

                        fx = U * l_ρinv[i, j, k] * Qmoist
                        fy = V * l_ρinv[i, j, k] * Qmoist
                        fz = W * l_ρinv[i, j, k] * Qmoist

                        s_F[i, j, k, 1] = MJ * (ξx * fx + ξy * fy + ξz * fz)
                        s_G[i, j, k, 1] = MJ * (ηx * fx + ηy * fy + ηz * fz)
                        s_H[i, j, k, 1] = MJ * (ζx * fx + ζy * fy + ζz * fz)
                    end
                end
            end

            @synchronize

            @loop for k in (1:Nq; threadIdx().z)
                @loop for j in (1:Nq; threadIdx().y)
                    @loop for i in (1:Nq; threadIdx().x)
                        # TODO: Prefetch MJI and rhs
                        MJI = vgeo[i, j, k, _MJI, e]

                        rhsmoist = zero(DFloat)
                        for n = 1:Nq
                            MJI_Dni = MJI * s_D[n, i]
                            MJI_Dnj = MJI * s_D[n, j]
                            MJI_Dnk = MJI * s_D[n, k]

                            rhsmoist += MJI_Dni * s_F[n, j, k, 1]
                            rhsmoist += MJI_Dnj * s_G[i, n, k, 1]
                            rhsmoist += MJI_Dnk * s_H[i, j, n, 1]
                        end
                        rhs[i, j, k, s, e] += rhsmoist
                    end
                end
            end
        end

        # Loop over trace variables
        # TODO: This should probably be unrolled by some factor
        rhstrace = zero(eltype(rhs))
        for m = 1:ntrace
            s = _nstate + nmoist + m

            @synchronize

            @loop for k in (1:Nq; threadIdx().z)
                @loop for j in (1:Nq; threadIdx().y)
                    @loop for i in (1:Nq; threadIdx().x)
                        MJ = vgeo[i, j, k, _MJ, e]
                        ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
                        ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
                        ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]

                        Qtrace = Q[i, j, k, s, e]
                        U = Q[i, j, k, _U, e]
                        V = Q[i, j, k, _V, e]
                        W = Q[i, j, k, _W, e]

                        fx = U * l_ρinv[i, j, k] * Qtrace
                        fy = V * l_ρinv[i, j, k] * Qtrace
                        fz = W * l_ρinv[i, j, k] * Qtrace

                        s_F[i, j, k, 1] = MJ * (ξx * fx + ξy * fy + ξz * fz)
                        s_G[i, j, k, 1] = MJ * (ηx * fx + ηy * fy + ηz * fz)
                        s_H[i, j, k, 1] = MJ * (ζx * fx + ζy * fy + ζz * fz)
                    end
                end
            end

            @synchronize

            @loop for k in (1:Nq; threadIdx().z)
                @loop for j in (1:Nq; threadIdx().y)
                    @loop for i in (1:Nq; threadIdx().x)
                        # TODO: Prefetch MJI and rhs
                        MJI = vgeo[i, j, k, _MJI, e]

                        rhstrace = zero(DFloat)
                        for n = 1:Nq
                            MJI_Dni = MJI * s_D[n, i]
                            MJI_Dnj = MJI * s_D[n, j]
                            MJI_Dnk = MJI * s_D[n, k]

                            rhstrace += MJI_Dni * s_F[n, j, k, 1]
                            rhstrace += MJI_Dnj * s_G[i, n, k, 1]
                            rhstrace += MJI_Dnk * s_H[i, j, n, 1]
                        end
                        rhs[i, j, k, s, e] += rhstrace
                    end
                end
            end
        end
    end
    nothing
end
# }}}

# {{{ volumerhs_v3!
function volumerhs_v3!(::Val{3},
                       ::Val{N},
                       ::Val{nmoist},
                       ::Val{ntrace},
                       rhs,
                       Q,
                       vgeo,
                       gravity,
                       D,
                       nelem) where {N, nmoist, ntrace}
    nvar = _nstate + nmoist + ntrace

    Nq = N + 1

    s_D = @shmem eltype(D) (Nq, Nq)
    s_F = @shmem eltype(Q) (Nq, Nq, _nstate)
    s_G = @shmem eltype(Q) (Nq, Nq, _nstate)
    s_H = @shmem eltype(Q) (Nq, Nq, _nstate)

    r_rhsρ = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsU = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsV = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsW = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsE = @scratch eltype(rhs) (Nq, Nq, Nq) 2

    @inbounds @loop for e in (1:nelem; blockIdx().x)
        @loop for j in (1:Nq; threadIdx().y)
            @loop for i in (1:Nq; threadIdx().x)
                 for k in 1:Nq
                     r_rhsρ[k, i, j] = zero(eltype(rhs))
                     r_rhsU[k, i, j] = zero(eltype(rhs))
                     r_rhsV[k, i, j] = zero(eltype(rhs))
                     r_rhsW[k, i, j] = zero(eltype(rhs))
                     r_rhsE[k, i, j] = zero(eltype(rhs))
                end

                # fetch D into shared
                s_D[i, j] = D[i, j]
            end
        end

        @unroll for k in 1:Nq
            @synchronize
            @loop for j in (1:Nq; threadIdx().y)
                @loop for i in (1:Nq; threadIdx().x)

                    # Load values will need into registers
                    MJ = vgeo[i, j, k, _MJ, e]
                    ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
                    ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
                    ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]
                    z = vgeo[i,j,k,_z,e]

                    U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
                    ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

                    P = gdm1*(E - (U^2 + V^2 + W^2)/(2*ρ) - ρ*gravity*z)

                    ρinv = 1 / ρ

                    fluxρ_x = U
                    fluxU_x = ρinv * U * U + P
                    fluxV_x = ρinv * U * V
                    fluxW_x = ρinv * U * W
                    fluxE_x = ρinv * U * (E + P)

                    fluxρ_y = V
                    fluxU_y = ρinv * V * U
                    fluxV_y = ρinv * V * V + P
                    fluxW_y = ρinv * V * W
                    fluxE_y = ρinv * V * (E + P)

                    fluxρ_z = W
                    fluxU_z = ρinv * W * U
                    fluxV_z = ρinv * W * V
                    fluxW_z = ρinv * W * W + P
                    fluxE_z = ρinv * W * (E + P)

                    s_F[i, j,  _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
                    s_F[i, j,  _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
                    s_F[i, j,  _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
                    s_F[i, j,  _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
                    s_F[i, j,  _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

                    s_G[i, j,  _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
                    s_G[i, j,  _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
                    s_G[i, j,  _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
                    s_G[i, j,  _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
                    s_G[i, j,  _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

                    r_Hρ = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
                    r_HU = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
                    r_HV = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
                    r_HW = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
                    r_HE = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)

                    # one shared access per 10 flops
                    for n = 1:Nq
                        Dkn = s_D[k, n]

                        r_rhsρ[n,i,j] += Dkn * r_Hρ
                        r_rhsU[n,i,j] += Dkn * r_HU
                        r_rhsV[n,i,j] += Dkn * r_HV
                        r_rhsW[n,i,j] += Dkn * r_HW
                        r_rhsE[n,i,j] += Dkn * r_HE
                    end

                    r_rhsW[k,i,j] -= MJ * ρ * gravity
                end
            end

            @synchronize

            @loop for j in (1:Nq; threadIdx().y)
                @loop for i in (1:Nq; threadIdx().x)

                    # loop of ξ-grid lines
                    for n = 1:Nq
                        Dni = s_D[n, i]
                        Dnj = s_D[n, j]

                        r_rhsρ[k,i,j] += Dni * s_F[n, j, _ρ]
                        r_rhsρ[k,i,j] += Dnj * s_G[i, n, _ρ]

                        r_rhsU[k,i,j] += Dni * s_F[n, j, _U]
                        r_rhsU[k,i,j] += Dnj * s_G[i, n, _U]

                        r_rhsV[k,i,j] += Dni * s_F[n, j, _V]
                        r_rhsV[k,i,j] += Dnj * s_G[i, n, _V]

                        r_rhsW[k,i,j] += Dni * s_F[n, j, _W]
                        r_rhsW[k,i,j] += Dnj * s_G[i, n, _W]

                        r_rhsE[k,i,j] += Dni * s_F[n, j, _E]
                        r_rhsE[k,i,j] += Dnj * s_G[i, n, _E]
                    end
                end
            end
        end # k

        @loop for j in (1:Nq; threadIdx().y)
            @loop for i in (1:Nq; threadIdx().x)
                for k in 1:Nq
                    MJI = vgeo[i, j, k, _MJI, e]

                    rhs[i, j, k, _U, e] += MJI*r_rhsU[k, i, j]
                    rhs[i, j, k, _V, e] += MJI*r_rhsV[k, i, j]
                    rhs[i, j, k, _W, e] += MJI*r_rhsW[k, i, j]
                    rhs[i, j, k, _ρ, e] += MJI*r_rhsρ[k, i, j]
                    rhs[i, j, k, _E, e] += MJI*r_rhsE[k, i, j]
                end
            end
        end
    end
    nothing
end
# }}}

# {{{ volumerhs_v3_1
function volumerhs_v3_1!(::Val{3},
                         ::Val{N},
                         ::Val{nmoist},
                         ::Val{ntrace},
                         rhs,
                         Q,
                         vgeo,
                         gravity,
                         D,
                         nelem) where {N, nmoist, ntrace}

    nvar = _nstate + nmoist + ntrace

    Nq = N + 1

    s_D = @shmem eltype(D) (Nq, Nq)
    s_F = @shmem eltype(Q) (Nq, Nq, _nstate)
    s_G = @shmem eltype(Q) (Nq, Nq, _nstate)
    s_H = @shmem eltype(Q) (Nq, Nq, _nstate)

    r_rhsρ = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsU = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsV = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsW = @scratch eltype(rhs) (Nq, Nq, Nq) 2
    r_rhsE = @scratch eltype(rhs) (Nq, Nq, Nq) 2

    @inbounds for e = blockIdx().x
        for j = threadIdx().y
            for i = threadIdx().x
                for k in 1:Nq
                    r_rhsρ[k, i, j] = zero(eltype(rhs))
                    r_rhsU[k, i, j] = zero(eltype(rhs))
                    r_rhsV[k, i, j] = zero(eltype(rhs))
                    r_rhsW[k, i, j] = zero(eltype(rhs))
                    r_rhsE[k, i, j] = zero(eltype(rhs))
                end

                # fetch D into shared
                s_D[i, j] = D[i, j]
            end
        end

        @unroll for k in 1:Nq
            @synchronize
            for j = threadIdx().y
                for i = threadIdx().x

                    # Load values will need into registers
                    MJ = vgeo[i, j, k, _MJ, e]
                    ξx, ξy, ξz = vgeo[i,j,k,_ξx,e], vgeo[i,j,k,_ξy,e], vgeo[i,j,k,_ξz,e]
                    ηx, ηy, ηz = vgeo[i,j,k,_ηx,e], vgeo[i,j,k,_ηy,e], vgeo[i,j,k,_ηz,e]
                    ζx, ζy, ζz = vgeo[i,j,k,_ζx,e], vgeo[i,j,k,_ζy,e], vgeo[i,j,k,_ζz,e]
                    z = vgeo[i,j,k,_z,e]

                    U, V, W = Q[i, j, k, _U, e], Q[i, j, k, _V, e], Q[i, j, k, _W, e]
                    ρ, E = Q[i, j, k, _ρ, e], Q[i, j, k, _E, e]

                    P = gdm1*(E - (U^2 + V^2 + W^2)/(2*ρ) - ρ*gravity*z)

                    ρinv = 1 / ρ

                    fluxρ_x = U
                    fluxU_x = ρinv * U * U + P
                    fluxV_x = ρinv * U * V
                    fluxW_x = ρinv * U * W
                    fluxE_x = ρinv * U * (E + P)

                    fluxρ_y = V
                    fluxU_y = ρinv * V * U
                    fluxV_y = ρinv * V * V + P
                    fluxW_y = ρinv * V * W
                    fluxE_y = ρinv * V * (E + P)

                    fluxρ_z = W
                    fluxU_z = ρinv * W * U
                    fluxV_z = ρinv * W * V
                    fluxW_z = ρinv * W * W + P
                    fluxE_z = ρinv * W * (E + P)

                    s_F[i, j,  _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y + ξz * fluxρ_z)
                    s_F[i, j,  _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y + ξz * fluxU_z)
                    s_F[i, j,  _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y + ξz * fluxV_z)
                    s_F[i, j,  _W] = MJ * (ξx * fluxW_x + ξy * fluxW_y + ξz * fluxW_z)
                    s_F[i, j,  _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y + ξz * fluxE_z)

                    s_G[i, j,  _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y + ηz * fluxρ_z)
                    s_G[i, j,  _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y + ηz * fluxU_z)
                    s_G[i, j,  _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y + ηz * fluxV_z)
                    s_G[i, j,  _W] = MJ * (ηx * fluxW_x + ηy * fluxW_y + ηz * fluxW_z)
                    s_G[i, j,  _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y + ηz * fluxE_z)

                    r_Hρ = MJ * (ζx * fluxρ_x + ζy * fluxρ_y + ζz * fluxρ_z)
                    r_HU = MJ * (ζx * fluxU_x + ζy * fluxU_y + ζz * fluxU_z)
                    r_HV = MJ * (ζx * fluxV_x + ζy * fluxV_y + ζz * fluxV_z)
                    r_HW = MJ * (ζx * fluxW_x + ζy * fluxW_y + ζz * fluxW_z)
                    r_HE = MJ * (ζx * fluxE_x + ζy * fluxE_y + ζz * fluxE_z)

                    # one shared access per 10 flops
                    for n = 1:Nq
                        Dkn = s_D[k, n]

                        r_rhsρ[n,i,j] += Dkn * r_Hρ
                        r_rhsU[n,i,j] += Dkn * r_HU
                        r_rhsV[n,i,j] += Dkn * r_HV
                        r_rhsW[n,i,j] += Dkn * r_HW
                        r_rhsE[n,i,j] += Dkn * r_HE
                    end

                    r_rhsW[k,i,j] -= MJ * ρ * gravity
                end
            end

            @synchronize

            for j = threadIdx().y
                for i = threadIdx().x

                    # loop of ξ-grid lines
                    for n = 1:Nq
                        Dni = s_D[n, i]
                        Dnj = s_D[n, j]

                        r_rhsρ[k,i,j] += Dni * s_F[n, j, _ρ]
                        r_rhsρ[k,i,j] += Dnj * s_G[i, n, _ρ]

                        r_rhsU[k,i,j] += Dni * s_F[n, j, _U]
                        r_rhsU[k,i,j] += Dnj * s_G[i, n, _U]

                        r_rhsV[k,i,j] += Dni * s_F[n, j, _V]
                        r_rhsV[k,i,j] += Dnj * s_G[i, n, _V]

                        r_rhsW[k,i,j] += Dni * s_F[n, j, _W]
                        r_rhsW[k,i,j] += Dnj * s_G[i, n, _W]

                        r_rhsE[k,i,j] += Dni * s_F[n, j, _E]
                        r_rhsE[k,i,j] += Dnj * s_G[i, n, _E]
                    end
                end
            end
        end # k

        for j = threadIdx().y
            for i = threadIdx().x
                for k in 1:Nq
                    MJI = vgeo[i, j, k, _MJI, e]

                    rhs[i, j, k, _U, e] += MJI*r_rhsU[k, i, j]
                    rhs[i, j, k, _V, e] += MJI*r_rhsV[k, i, j]
                    rhs[i, j, k, _W, e] += MJI*r_rhsW[k, i, j]
                    rhs[i, j, k, _ρ, e] += MJI*r_rhsρ[k, i, j]
                    rhs[i, j, k, _E, e] += MJI*r_rhsE[k, i, j]
                end
            end
        end
    end
    nothing
end
# }}}

# {{{ main
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--Float64"
            help = "Run with Float64 data (otherwise run with Float32)"
            action = :store_true
        "--ntrials"
            help = "number trials for each kernel"
            arg_type = Int
            default = 100
        "N"
            help = "polynomial order"
            required = true
            arg_type = Int
        "nelem"
            help = "number of elements"
            required = true
            arg_type = Int
    end

    return parse_args(s)
end

function main()
    CUDAnative.timings!()

    parsed_args = parse_commandline()

    N = parsed_args["N"]
    nelem = parsed_args["nelem"]
    ntrials = parsed_args["ntrials"]
    DFloat = parsed_args["Float64"] ? Float64 : Float32

    rnd = MersenneTwister(0)
    nmoist = 0
    ntrace = 0
    nvar = _nstate + nmoist + ntrace

    Nq = N + 1
    Q = 1 .+ CuArray(rand(rnd, DFloat, Nq, Nq, Nq, nvar, nelem))
    Q[:, :, :, _E, :] .+= 20
    vgeo = CuArray(rand(rnd, DFloat, Nq, Nq, Nq, _nvgeo, nelem))

    # Make sure the entries of the mass matrix satisfy the inverse relation
    vgeo[:, :, :, _MJ, :] .+= 3
    vgeo[:, :, :, _MJI, :] .= 1 ./ vgeo[:, :, :, _MJ, :]

    D = CuArray(rand(rnd, DFloat, Nq, Nq))

    rhs = CuArray(zeros(DFloat, Nq, Nq, Nq, nvar, nelem))

    #
    # volumerhs_v2!
    #
    fill!(rhs, 0)
    @launch(CUDA(), threads=(N+1, N+1, N+1), blocks=nelem,
          volumerhs_v2!(Val(3), Val(N), Val(nmoist), Val(ntrace),
                        rhs, Q, vgeo, DFloat(grav), D, nelem))
    norm_v2 = norm(rhs)
    CUDAdrv.@profile for _ = 1:ntrials
         @launch(CUDA(), threads=(N+1, N+1, N+1), blocks=nelem,
               volumerhs_v2!(Val(3), Val(N), Val(nmoist), Val(ntrace), rhs,
                             Q, vgeo, DFloat(grav), D, nelem))
    end

    #
    # volumerhs_v3!
    #
    fill!(rhs, 0)
    @launch(CUDA(), threads=(N+1, N+1), blocks=nelem,
            volumerhs_v3!(Val(3),
                       Val(N),
                       Val(nmoist),
                       Val(ntrace),
                       rhs,
                       Q,
                       vgeo,
                       DFloat(grav),
                       D,
                       nelem))
    norm_v3 = norm(rhs)
    CUDAdrv.@profile for _ = 1:ntrials
        @launch(CUDA(), threads=(N+1, N+1), blocks=nelem, maxregs=255,
              volumerhs_v3!(Val(3), Val(N), Val(nmoist), Val(ntrace), rhs,
                            Q, vgeo, DFloat(grav), D, nelem))
    end

    #
    # volumerhs_v3_1!
    #
    fill!(rhs, 0)
    @launch(CUDA(), threads=(N+1, N+1), blocks=nelem, maxregs=255,
          volumerhs_v3_1!(Val(3), Val(N), Val(nmoist), Val(ntrace),
                          rhs, Q, vgeo, DFloat(grav), D, nelem))
    norm_v3_1 = norm(rhs)
    CUDAdrv.@profile for _ = 1:ntrials
        @launch(CUDA(), threads=(N+1, N+1), blocks=nelem, maxregs=255,
              volumerhs_v3_1!(Val(3), Val(N), Val(nmoist), Val(ntrace),
                              rhs, Q, vgeo, DFloat(grav), D, nelem))
    end

    @assert norm_v2 ≈ norm_v3
    @assert norm_v2 ≈ norm_v3_1

    CUDAnative.timings()

    nothing
end
# }}}

main()


text:
@testset "ResNet" begin
  for MODEL in [
    ResNet18,
    ResNet34,
    ResNet50,
    ResNet101,
    ResNet152,
    ]
    model = MODEL()
    @test size(model(rand(256, 256, 3, 50))) == (1000, 50)
  end
end