text:
import Base: append!, ndims, length, isempty, getindex
struct EGO{S,KRG<:Kriging,AQ<:EGOAquisition,AC<:Acquire,LB,UB,NC}
    krg::KRG
    tuner::KrigingTuner
    aquisition::AQ
    acquire::AC
    lowerbounds::LB
    upperbounds::UB
    constrains::NTuple{NC,Function}
    function EGO(
        kriging::Kriging{N,1},
        tuner::KrigingTuner,
        lowerbounds::AbstractVector,
        upperbounds::AbstractVector;
        acquire::Acquire=Acquire(),
        aquisition::SingleObjAquisition{S}=ExpectedImprovement(:Min),
        constrains::NTuple{NC,Function}=(),
    ) where {N,S,NC}
        length(lowerbounds) != N &&
            throw(ArgumentError("length of `lowerbounds` must be $N"))
        length(upperbounds) != N &&
            throw(ArgumentError("length of `upperbounds` must be $N"))
        lob = copy(lowerbounds)
        upb = copy(upperbounds)
        KRG = typeof(kriging)
        aq = NC == 0 ? aquisition : ConstrainedAquisition(aquisition, constrains)
        AQ = typeof(aq)
        AC = typeof(acquire)
        LB = typeof(lob)
        UB = typeof(upb)
        return new{S,KRG,AQ,AC,LB,UB,NC}(kriging, tuner, aq, acquire, lob, upb, constrains)
    end
end

Base.length(ego::EGO) = length(ego.krg.gps[1].y)
Base.isempty(ego::EGO) = isempty(ego.krg.gps[1].y)
lowerbounds(ego::EGO) = ego.lowerbounds
upperbounds(ego::EGO) = ego.upperbounds
constrains(ego::EGO) = ego.constrains

function Base.append!(ego::EGO, x::AbstractMatrix, y::AbstractVector)
    try
        append!(ego.krg, x, transpose(y))
        update_parameters!(ego.aquisition, ego.krg)
        tunekrg!(ego.krg, ego.tuner)
    catch e
        @info "`append!` failed, error=$e, calling `forcetune!`"
        forcetune!(ego.krg, ego.tuner)
    end
end

Base.append!(ego::EGO, x::AbstractVector, y::Float64) = append!(ego, reshape(x, :, 1), [y])

optimum(ego::EGO{Min}) = ego.krg[argmin(ego.krg.gps[1].y)]

optimum(ego::EGO{Max}) = ego.krg[argmax(ego.krg.gps[1].y)]

function acquire(
    ego::EGO, x_cache::AbstractVector...; retry::Integer=10, verbose::Bool=false
)
    try
        return acquire(
            ego.krg,
            ego.aquisition,
            ego.acquire,
            ego.lowerbounds,
            ego.upperbounds,
            x_cache...,
        )
    catch e
        if retry > 1
            verbose && @info "acquire() failed, error=$e, retry($(retry-1))"
            return acquire(ego, x_cache...; retry=retry - 1, verbose=verbose)
        else
            if isa(e, NoFeasibleInfill)
                rethrow(e)
            else
                throw(NoFeasibleInfill(sprint(print, e)))
            end
        end
    end
end

history(ego::EGO) = get_samples(ego.krg)


text:
function gaussPoints1D(n::Int64)
    # Gauss points
    if n == 2
        gpCoord = [-1/sqrt(3) 1/sqrt(3)]
        gpW = [1.0 1.0]
    elseif n == 3
        gpCoord = [-sqrt(3/5) 0 sqrt(3/5)]
        gpW = [5/9 8/9 5/9]
    elseif n == 4
        gpCoord = [-sqrt(3/7 + 2/7*sqrt(6/5)) -sqrt(3/7 - 2/7*sqrt(6/5)) sqrt(3/7 - 2/7*sqrt(6/5)) sqrt(3/7 + 2/7*sqrt(6/5))]
        gpW = [(18 - sqrt(30))/36 (18 + sqrt(30))/36 (18 + sqrt(30))/36 (18 - sqrt(30))/36]
    end
    return hcat(gpW', gpCoord')
end

function gaussPoints2D(n::Int64)
    # Gauss points
    coords, weights, x, w = zeros(n, n, 2), zeros(n, n), [], []
    if n == 2
        x = [-1 / sqrt(3), 1 / sqrt(3)]
        w = [1.0, 1.0]
    elseif n == 3
        x = [-sqrt(3 / 5), 0, sqrt(3 / 5)]
        w = [5 / 9, 8 / 9, 5 / 9]
    elseif n == 4
        x = [-sqrt(3 / 7 + 2 / 7 * sqrt(6 / 5)), -sqrt(3 / 7 - 2 / 7 * sqrt(6 / 5)), sqrt(3 / 7 - 2 / 7 * sqrt(6 / 5)), sqrt(3 / 7 + 2 / 7 * sqrt(6 / 5))]
        w = [(18 - sqrt(30)) / 36, (18 + sqrt(30)) / 36, (18 + sqrt(30)) / 36, (18 - sqrt(30)) / 36]
    end
    for i in 1:n
        for j = 1:n
            coords[i, j, :] = [x[i], x[j]]
            weights[i, j] = w[i] * w[j]
        end
    end
    gpCoord = reshape(coords, n * n, 2)
    gpW = reshape(weights, n * n, 1)
    return hcat(gpW, gpCoord)
end

ext: jl

max_stars_repo_name: likemaestro/Nonlinear-Homogenization

max_issues_repo_name: likemaestro/Nonlinear-Homogenization

max_forks_repo_name: likemaestro/Nonlinear-Homogenization


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 = "PGPLOT"
version = v"5.2"

sources = [ArchiveSource("ftp://ftp.astro.caltech.edu/pub/pgplot/pgplot5.2.tar.gz", "a5799ff719a510d84d26df4ae7409ae61fe66477e3f1e8820422a9a4727a5be4")]

script = raw"""
cd $WORKSPACE/srcdir
mkdir pgplot_build && cd pgplot_build/
cat ../pgplot/drivers.list | sed 's|! PSDRIV|  PSDRIV|g' | sed 's|! GIDRIV|  GIDRIV|g' > drivers.list
../pgplot/makemake ../pgplot/ linux g77_gcc
sed -i 's|FCOMPL=g77|FCOMPL=gfortran|' makefile
make lib SHARED_LD="${FC} -shared  -o libpgplot.${dlext}"
cp libpgplot.so "${libdir}/libpgplot.so"
install_license ../pgplot/copyright.notice
"""

platforms = [
    Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("i686", "linux"; libc="glibc"),
]
platforms = expand_gfortran_versions(platforms)

products = [LibraryProduct("libpgplot", :libpgplot)]

dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# TODO: Move to the MPIStateArrays
# {{{ MPI Buffer handling
fillsendbuf!(h, d, b::MPIStateArray, e) = fillsendbuf!(h, d, b.Q, e)
transferrecvbuf!(h, d, b::MPIStateArray, e) = transferrecvbuf!(h, d, b.Q, e)

function fillsendbuf!(host_sendbuf, device_sendbuf::Array, buf::Array, sendelems)
  host_sendbuf[:, :, :] .= buf[:, :, sendelems]
end

function transferrecvbuf!(device_recvbuf::Array, host_recvbuf, buf::Array,
                          nrealelem)
  buf[:, :, nrealelem+1:end] .= host_recvbuf[:, :, :]
end
# }}}

# TODO 2 dimensions convert z to the vertical coordinate
# {{{ Volume Gradient for 2-D
function volumegrad!(::Val{2}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                     grad::Array, Q, vgeo, gravity, D,
                     elems) where {N, nmoist, ntrace}
  dim = 2
  DFloat = eltype(Q)

  nvar = _nstate + nmoist + ntrace
  ngrad = _nstategrad + 3nmoist

  Nq = N + 1

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, nvar, nelem)
  grad = reshape(grad, Nq, Nq, ngrad, nelem)
  vgeo = reshape(vgeo, Nq, Nq, size(vgeo,2), nelem)

  s_q = s_ρ = Array{DFloat}(undef, Nq, Nq)
  s_u = Array{DFloat}(undef, Nq, Nq)
  s_v = Array{DFloat}(undef, Nq, Nq)
  s_T = Array{DFloat}(undef, Nq, Nq)

  #Initialise q_m vector of moist variables
  #Requires at least 3 elements for q_t, q_liq, q_ice
  q_m = zeros(DFloat, max(3, nmoist))

  @inbounds for e in elems
    for j = 1:Nq, i = 1:Nq
      U, V = Q[i, j, _U, e], Q[i, j, _V, e]
      ρ, E = Q[i, j, _ρ, e], Q[i, j, _E, e]
      y = vgeo[i,j,_y,e]
      E_int = E - (U^2 + V^2)/(2*ρ) - ρ * gravity * y
      # Get specific humidity quantities from state vector
      # per unit mass conversion required for sat_adjust only
      for m = 1:nmoist
          s = _nstate + m
          q_m[m]  = Q[i, j, s, e] / ρ
      end
      # Saturation temperature to obtain temperature assuming thermodynamic equilibrium
      T = saturation_adjustment(E_int/ρ, ρ, q_m[1])
      # TODO: Possibility of carrying q_liq and q_ice through state vector to include non-equilibrium thermodynamics (?)
      q = PhasePartition_equil(T, ρ, q_m[1])
      s_ρ[i, j] = ρ
      s_u[i, j] = U/ρ
      s_v[i, j] = V/ρ
      s_T[i, j] = T
    end

    for j = 1:Nq, i = 1:Nq
      ρξ = ρη = zero(DFloat)
      uξ = uη = zero(DFloat)
      vξ = vη = zero(DFloat)
      Tξ = Tη = zero(DFloat)

      for n = 1:Nq
        ρξ += D[i, n] * s_ρ[n, j]
        ρη += D[j, n] * s_ρ[i, n]

        uξ += D[i, n] * s_u[n, j]
        uη += D[j, n] * s_u[i, n]

        vξ += D[i, n] * s_v[n, j]
        vη += D[j, n] * s_v[i, n]

        Tξ += D[i, n] * s_T[n, j]
        Tη += D[j, n] * s_T[i, n]
      end

      ξx, ξy = vgeo[i,j,_ξx,e], vgeo[i,j,_ξy,e]
      ηx, ηy = vgeo[i,j,_ηx,e], vgeo[i,j,_ηy,e]

      ρx = ξx*ρξ + ηx*ρη
      ρy = ξy*ρξ + ηy*ρη

      ux = ξx*uξ + ηx*uη
      uy = ξy*uξ + ηy*uη

      vx = ξx*vξ + ηx*vη
      vy = ξy*vξ + ηy*vη

      Tx = ξx*Tξ + ηx*Tη
      Ty = ξy*Tξ + ηy*Tη

      grad[i, j, _ρx, e] = ρx
      grad[i, j, _ρy, e] = ρy
      grad[i, j, _ρz, e] = zero(DFloat)

      grad[i, j, _ux, e] = ux
      grad[i, j, _uy, e] = uy
      grad[i, j, _uz, e] = zero(DFloat)

      grad[i, j, _vx, e] = vx
      grad[i, j, _vy, e] = vy
      grad[i, j, _vz, e] = zero(DFloat)

      grad[i, j, _wx, e] = zero(DFloat)
      grad[i, j, _wy, e] = zero(DFloat)
      grad[i, j, _wz, e] = zero(DFloat)

      grad[i, j, _Tx, e] = Tx
      grad[i, j, _Ty, e] = Ty
      grad[i, j, _Tz, e] = zero(DFloat)
    end

    # loop over moist variables
    for m = 1:nmoist
      s = _nstate + m
      sout = _nstategrad + 3*(m-1)

      for j = 1:Nq, i = 1:Nq
        s_q[i, j] =  Q[i, j, s, e]
      end

      for j = 1:Nq, i = 1:Nq
        qξ = qη = zero(DFloat)

        for n = 1:Nq
          qξ += D[i, n] * s_q[n, j]
          qη += D[j, n] * s_q[i, n]
        end

        ξx, ξy = vgeo[i, j, _ξx, e], vgeo[i, j, _ξy, e]
        ηx, ηy = vgeo[i, j, _ηx, e], vgeo[i, j, _ηy, e]

        qx = ξx*qξ + ηx*qη
        qy = ξy*qξ + ηy*qη

        grad[i, j, sout + 1, e] = qx
        grad[i, j, sout + 2, e] = qy
        grad[i, j, sout + 3, e] = zero(DFloat)
      end
    end
  end
end
# }}}

# {{{ Volume Gradient for 3-D
function volumegrad!(::Val{3}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                     grad::Array, Q, vgeo, gravity, D,
                     elems) where {N, nmoist, ntrace}
  dim = 3
  DFloat = eltype(Q)

  nvar = _nstate + nmoist + ntrace
  ngrad = _nstategrad + 3nmoist

  Nq = N + 1

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, Nq, nvar, nelem)
  grad = reshape(grad, Nq, Nq, Nq, ngrad, nelem)
  vgeo = reshape(vgeo, Nq, Nq, Nq, size(vgeo,2), nelem)

  s_q = s_ρ = Array{DFloat}(undef, Nq, Nq, Nq)
  s_u = Array{DFloat}(undef, Nq, Nq, Nq)
  s_v = Array{DFloat}(undef, Nq, Nq, Nq)
  s_w = Array{DFloat}(undef, Nq, Nq, Nq)
  s_T = Array{DFloat}(undef, Nq, Nq, Nq)

  q_m = zeros(DFloat, max(3, nmoist))

  @inbounds for e in elems
    for k = 1:Nq, j = 1:Nq, i = 1:Nq
      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]
      z = vgeo[i, j, k, _z, e]
      E_int = E - (U^2 + V^2 + W^2)/(2*ρ) - ρ * gravity * z

      for m = 1:nmoist
         s = _nstate + m
         q_m[m] = Q[i, j, k, s, e] / ρ
      end

      T = saturation_adjustment(E_int/ρ, ρ, q_m[1])
      # TODO: Possibility of carrying q_liq and q_ice through state vector to include non-equilibrium thermodynamics
      q = PhasePartition_equil(T, ρ, q_m[1])
      s_ρ[i, j, k] = ρ
      s_u[i, j, k] = U/ρ
      s_v[i, j, k] = V/ρ
      s_w[i, j, k] = W/ρ
      s_T[i, j, k] = air_temperature(E_int/ρ, q)
    end

    for k = 1:Nq, j = 1:Nq, i = 1:Nq
      ρξ = ρη = ρζ = zero(DFloat)
      uξ = uη = uζ = zero(DFloat)
      vξ = vη = vζ = zero(DFloat)
      wξ = wη = wζ = zero(DFloat)
      Tξ = Tη = Tζ = zero(DFloat)

      for n = 1:Nq
        ρξ += D[i, n] * s_ρ[n, j, k]
        ρη += D[j, n] * s_ρ[i, n, k]
        ρζ += D[k, n] * s_ρ[i, j, n]

        uξ += D[i, n] * s_u[n, j, k]
        uη += D[j, n] * s_u[i, n, k]
        uζ += D[k, n] * s_u[i, j, n]

        vξ += D[i, n] * s_v[n, j, k]
        vη += D[j, n] * s_v[i, n, k]
        vζ += D[k, n] * s_v[i, j, n]

        wξ += D[i, n] * s_w[n, j, k]
        wη += D[j, n] * s_w[i, n, k]
        wζ += D[k, n] * s_w[i, j, n]

        Tξ += D[i, n] * s_T[n, j, k]
        Tη += D[j, n] * s_T[i, n, k]
        Tζ += D[k, n] * s_T[i, j, n]
      end

      ξ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]

      ρx = ξx*ρξ + ηx*ρη + ζx*ρζ
      ρy = ξy*ρξ + ηy*ρη + ζy*ρζ
      ρz = ξz*ρξ + ηz*ρη + ζz*ρζ

      ux = ξx*uξ + ηx*uη + ζx*uζ
      uy = ξy*uξ + ηy*uη + ζy*uζ
      uz = ξz*uξ + ηz*uη + ζz*uζ

      vx = ξx*vξ + ηx*vη + ζx*vζ
      vy = ξy*vξ + ηy*vη + ζy*vζ
      vz = ξz*vξ + ηz*vη + ζz*vζ

      wx = ξx*wξ + ηx*wη + ζx*wζ
      wy = ξy*wξ + ηy*wη + ζy*wζ
      wz = ξz*wξ + ηz*wη + ζz*wζ

      Tx = ξx*Tξ + ηx*Tη + ζx*Tζ
      Ty = ξy*Tξ + ηy*Tη + ζy*Tζ
      Tz = ξz*Tξ + ηz*Tη + ζz*Tζ

      grad[i, j, k, _ρx, e] = ρx
      grad[i, j, k, _ρy, e] = ρy
      grad[i, j, k, _ρz, e] = ρz

      grad[i, j, k, _ux, e] = ux
      grad[i, j, k, _uy, e] = uy
      grad[i, j, k, _uz, e] = uz

      grad[i, j, k, _vx, e] = vx
      grad[i, j, k, _vy, e] = vy
      grad[i, j, k, _vz, e] = vz

      grad[i, j, k, _wx, e] = wx
      grad[i, j, k, _wy, e] = wy
      grad[i, j, k, _wz, e] = wz

      grad[i, j, k, _Tx, e] = Tx
      grad[i, j, k, _Ty, e] = Ty
      grad[i, j, k, _Tz, e] = Tz
    end

    # loop over moist variables
    for m = 1:nmoist
      s = _nstate + m
      sout = _nstategrad + 3*(m-1)

      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        s_q[i, j, k] =  Q[i, j, k, s, e]
      end

      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        qξ = qη = qζ = zero(DFloat)

        for n = 1:Nq
          qξ += D[i, n] * s_q[n, j, k]
          qη += D[j, n] * s_q[i, n, k]
          qζ += D[k, n] * s_q[i, j, n]
        end

        ξ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]

        qx = ξx*qξ + ηx*qη + ζx*qζ
        qy = ξy*qξ + ηy*qη + ζy*qζ
        qz = ξz*qξ + ηz*qη + ζz*qζ

        grad[i, j, k, sout + 1, e] = qx
        grad[i, j, k, sout + 2, e] = qy
        grad[i, j, k, sout + 3, e] = qz
      end
    end
  end
end
# }}}

# {{{ Face Gradient (all dimensions)
function facegrad!(::Val{dim}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                   grad::Array, Q, vgeo, sgeo, gravity, elems, vmapM,
                   vmapP, elemtobndy) where {N, dim, nmoist, ntrace}
  DFloat = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

    q_mM = zeros(DFloat, max(3, nmoist))
    q_mP = zeros(DFloat, max(3, nmoist))

  @inbounds for e in elems
    for f = 1:nface
      for n = 1:Nfp
        (nxM, nyM, nzM, sMJ, vMJI) = sgeo[:, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        ρM = Q[vidM, _ρ, eM]
        UM = Q[vidM, _U, eM]
        VM = Q[vidM, _V, eM]
        WM = Q[vidM, _W, eM]
        EM = Q[vidM, _E, eM]
        yorzM = (dim == 2) ? vgeo[vidM, _y, eM] : vgeo[vidM, _z, eM]
        E_intM = EM - (UM^2 + VM^2+ WM^2)/(2*ρM) - ρM * gravity * yorzM

	for m = 1:nmoist
            s = _nstate + m
            q_mM[m] = Q[vidM, s, eM] / ρM
        end

        uM=UM/ρM
        vM=VM/ρM
        wM=WM/ρM

        # Saturation Adjustment
        TM = saturation_adjustment(E_intM/ρM, ρM, q_mM[1])
        qM = PhasePartition_equil(TM, ρM, q_mM[1])
        PM = air_pressure(TM, ρM, qM)
        bc = elemtobndy[f, e]

        if bc == 0
          ρP = Q[vidP, _ρ, eP]
          UP = Q[vidP, _U, eP]
          VP = Q[vidP, _V, eP]
          WP = Q[vidP, _W, eP]
          EP = Q[vidP, _E, eP]
          yorzP = (dim == 2) ? vgeo[vidP, _y, eP] : vgeo[vidP, _z, eP]
          for m = 1:nmoist
              s = _nstate + m
              q_mP[m] = Q[vidP, s, eP] / ρP
           end

          E_intP = EP - (UP^2 + VP^2+ WP^2)/(2*ρP) - ρP * gravity * yorzP

          # Saturation Adjustment
          TP = saturation_adjustment(E_intP/ρP, ρP, q_mP[1])
          qP = PhasePartition_equil(TP, ρP, q_mP[1])
          PP = air_pressure(TP, ρP, qP)
          uP = UP/ρP
          vP = VP/ρP
          wP = WP/ρP

        elseif bc == 1

          UnM = nxM * UM + nyM * VM + nzM * WM
          UP = UM - 2 * UnM * nxM
          VP = VM - 2 * UnM * nyM
          WP = WM - 2 * UnM * nzM
          ρP = ρM
          EP = EM
          PP = PM
          uP = UP/ρP
          vP = VP/ρP
          wP = WP/ρP
          TP = TM

        else
          error("Invalid boundary conditions $bc on face $f of element $e")

        end

        fluxρS = (ρP - ρM)/2
        fluxuS = (uP - uM)/2
        fluxvS = (vP - vM)/2
        fluxwS = (wP - wM)/2
        fluxTS = (TP - TM)/2

        grad[vidM, _ρx, eM] += vMJI*sMJ*nxM*fluxρS
        grad[vidM, _ρy, eM] += vMJI*sMJ*nyM*fluxρS
        grad[vidM, _ρz, eM] += vMJI*sMJ*nzM*fluxρS
        grad[vidM, _ux, eM] += vMJI*sMJ*nxM*fluxuS
        grad[vidM, _uy, eM] += vMJI*sMJ*nyM*fluxuS
        grad[vidM, _uz, eM] += vMJI*sMJ*nzM*fluxuS
        grad[vidM, _vx, eM] += vMJI*sMJ*nxM*fluxvS
        grad[vidM, _vy, eM] += vMJI*sMJ*nyM*fluxvS
        grad[vidM, _vz, eM] += vMJI*sMJ*nzM*fluxvS
        grad[vidM, _wx, eM] += vMJI*sMJ*nxM*fluxwS
        grad[vidM, _wy, eM] += vMJI*sMJ*nyM*fluxwS
        grad[vidM, _wz, eM] += vMJI*sMJ*nzM*fluxwS
        grad[vidM, _Tx, eM] += vMJI*sMJ*nxM*fluxTS
        grad[vidM, _Ty, eM] += vMJI*sMJ*nyM*fluxTS
        grad[vidM, _Tz, eM] += vMJI*sMJ*nzM*fluxTS

        # loop over moist variables
        for m = 1:nmoist
          s = _nstate + m
          ss = _nstategrad + 3*(m-1)
          qM, qP = Q[vidM, s, eM], Q[vidP, s, eP]

          fluxqS = (qP - qM)/2

          grad[vidM, ss+1, eM] += vMJI*sMJ*nxM*fluxqS
          grad[vidM, ss+2, eM] += vMJI*sMJ*nyM*fluxqS
          grad[vidM, ss+3, eM] += vMJI*sMJ*nzM*fluxqS
        end
      end
    end
  end
end
# }}}

# {{{ Volume RHS for 2-D
function volumerhs!(::Val{2}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                    rhs::Array, Q, grad, vgeo, gravity, viscosity, D,
                    elems) where {N, nmoist, ntrace}
  DFloat = eltype(Q)

  nvar = _nstate + nmoist + ntrace
  ngrad = _nstategrad + 3nmoist

  Nq = N + 1

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, nvar, nelem)
  grad = reshape(grad, Nq, Nq, ngrad, nelem)
  rhs = reshape(rhs, Nq, Nq, nvar, nelem)
  vgeo = reshape(vgeo, Nq, Nq, _nvgeo, nelem)

  s_F = Array{DFloat}(undef, Nq, Nq, _nstate)
  s_G = Array{DFloat}(undef, Nq, Nq, _nstate)
  l_u = Array{DFloat}(undef, Nq, Nq)
  l_v = Array{DFloat}(undef, Nq, Nq)

    q_m = zeros(DFloat, max(3, nmoist))

  @inbounds for e in elems
    for j = 1:Nq, i = 1:Nq
      MJ = vgeo[i, j, _MJ, e]
      ξx, ξy = vgeo[i,j,_ξx,e], vgeo[i,j,_ξy,e]
      ηx, ηy = vgeo[i,j,_ηx,e], vgeo[i,j,_ηy,e]
      y = vgeo[i,j,_y,e]

      U, V = Q[i, j, _U, e], Q[i, j, _V, e]
      ρ, E = Q[i, j, _ρ, e], Q[i, j, _E, e]
      E_int = E - (U^2 + V^2)/(2*ρ) - ρ * gravity * y

      for m = 1:nmoist
          s = _nstate + m
          q_m[m] = Q[i, j, s, e] / ρ
      end

      # Returns temperature after saturation adjustment
      # Required for phase-partitioning to find q_liq, q_ice
      T = saturation_adjustment(E_int/ρ, ρ, q_m[1])
      q = PhasePartition_equil(T, ρ, q_m[1])
      P = air_pressure(T, ρ, q)

      ρx, ρy = grad[i,j,_ρx,e], grad[i,j,_ρy,e]
      ux, uy = grad[i,j,_ux,e], grad[i,j,_uy,e]
      vx, vy = grad[i,j,_vx,e], grad[i,j,_vy,e]
      Tx, Ty = grad[i,j,_Tx,e], grad[i,j,_Ty,e]

      ρinv = 1 / ρ

      ldivu = λ_stokes*(ux + vy)
      u = ρinv*U
      v = ρinv*V

      vfluxρ_x = 0*ρx
      vfluxρ_y = 0*ρy

      vfluxU_x = 2*ux + ldivu
      vfluxU_y = uy + vx

      vfluxV_x = vx + uy
      vfluxV_y = 2*vy + ldivu

      vfluxE_x = u*(2*ux + ldivu) + v*(uy + vx) + k_μ*Tx
      vfluxE_y = u*(vx + uy) + v*(2*vy + ldivu) + k_μ*Ty

      fluxρ_x = U            - viscosity*vfluxρ_x
      fluxU_x = u * U + P    - viscosity*vfluxU_x
      fluxV_x = u * V        - viscosity*vfluxV_x
      fluxE_x = u * (E + P)  - viscosity*vfluxE_x

      fluxρ_y = V            - viscosity*vfluxρ_y
      fluxU_y = v * U        - viscosity*vfluxU_y
      fluxV_y = v * V + P    - viscosity*vfluxV_y
      fluxE_y = v * (E + P)  - viscosity*vfluxE_y

      s_F[i, j, _ρ] = MJ * (ξx * fluxρ_x + ξy * fluxρ_y)
      s_F[i, j, _U] = MJ * (ξx * fluxU_x + ξy * fluxU_y)
      s_F[i, j, _V] = MJ * (ξx * fluxV_x + ξy * fluxV_y)
      s_F[i, j, _W] = 0
      s_F[i, j, _E] = MJ * (ξx * fluxE_x + ξy * fluxE_y)

      s_G[i, j, _ρ] = MJ * (ηx * fluxρ_x + ηy * fluxρ_y)
      s_G[i, j, _U] = MJ * (ηx * fluxU_x + ηy * fluxU_y)
      s_G[i, j, _V] = MJ * (ηx * fluxV_x + ηy * fluxV_y)
      s_G[i, j, _W] = 0
      s_G[i, j, _E] = MJ * (ηx * fluxE_x + ηy * fluxE_y)

      # buoyancy term
      rhs[i, j, _V, e] -= ρ * gravity

      # Store velocity
      l_u[i, j], l_v[i, j] = u, v
    end

    # loop of ξ-grid lines
    for s = 1:_nstate, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, _MJI, e]
      for n = 1:Nq
        rhs[i, j, s, e] += MJI * D[n, i] * s_F[n, j, s]
      end
    end
    # loop of η-grid lines
    for s = 1:_nstate, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, _MJI, e]
      for n = 1:Nq
        rhs[i, j, s, e] += MJI * D[n, j] * s_G[i, n, s]
      end
    end

    # loop over moist variables
    # FIXME: Currently just passive advection
    for m = 1:nmoist
      s = _nstate + m
      ss = _nstategrad + 3*(m-1)

      for j = 1:Nq, i = 1:Nq
        MJ = vgeo[i, j, _MJ, e]
        ξx, ξy = vgeo[i,j,_ξx,e], vgeo[i,j,_ξy,e]
        ηx, ηy = vgeo[i,j,_ηx,e], vgeo[i,j,_ηy,e]
        u, v = l_u[i, j], l_v[i, j]

        q = Q[i, j, s, e]
        qx = grad[i, j, ss+1, e]
        qy = grad[i, j, ss+2, e]

        fx = u * q - viscosity*0*qx
        fy = v * q - viscosity*0*qy

        s_F[i, j, 1] = MJ * (ξx * fx + ξy * fy)
        s_G[i, j, 1] = MJ * (ηx * fx + ηy * fy)
      end
      for j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, _MJI, e]
        for n = 1:Nq
          rhs[i, j, s, e] += MJI * D[n, i] * s_F[n, j, 1]
        end
      end
      for j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, _MJI, e]
        for n = 1:Nq
          rhs[i, j, s, e] += MJI * D[n, j] * s_G[i, n, 1]
        end
      end
    end

    # loop over tracer variables
    for t = 1:ntrace
      s = _nstate + nmoist + t

      for j = 1:Nq, i = 1:Nq
        MJ = vgeo[i, j, _MJ, e]
        ξx, ξy = vgeo[i,j,_ξx,e], vgeo[i,j,_ξy,e]
        ηx, ηy = vgeo[i,j,_ηx,e], vgeo[i,j,_ηy,e]
        u, v = l_u[i, j], l_v[i, j]

        fx = u * Q[i, j, s, e]
        fy = v * Q[i, j, s, e]

        s_F[i, j, 1] = MJ * (ξx * fx + ξy * fy)
        s_G[i, j, 1] = MJ * (ηx * fx + ηy * fy)
      end
      for j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, _MJI, e]
        for n = 1:Nq
          rhs[i, j, s, e] += MJI * D[n, i] * s_F[n, j, 1]
        end
      end
      for j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, _MJI, e]
        for n = 1:Nq
          rhs[i, j, s, e] += MJI * D[n, j] * s_G[i, n, 1]
        end
      end
    end
  end
end
# }}}

# {{{ Volume RHS for 3-D
function volumerhs!(::Val{3}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                    rhs::Array, Q, grad, vgeo, gravity, viscosity, D,
                    elems) where {N, nmoist, ntrace}
  DFloat = eltype(Q)

  nvar = _nstate + nmoist + ntrace
  ngrad = _nstategrad + 3nmoist

  Nq = N + 1

  nelem = size(Q)[end]

  Q = reshape(Q, Nq, Nq, Nq, nvar, nelem)
  grad = reshape(grad, Nq, Nq, Nq, ngrad, nelem)
  rhs = reshape(rhs, Nq, Nq, Nq, nvar, nelem)
  vgeo = reshape(vgeo, Nq, Nq, Nq, _nvgeo, nelem)

  s_F = Array{DFloat}(undef, Nq, Nq, Nq, _nstate)
  s_G = Array{DFloat}(undef, Nq, Nq, Nq, _nstate)
  s_H = Array{DFloat}(undef, Nq, Nq, Nq, _nstate)
  l_u = Array{DFloat}(undef, Nq, Nq, Nq)
  l_v = Array{DFloat}(undef, Nq, Nq, Nq)
  l_w = Array{DFloat}(undef, Nq, Nq, Nq)

  q_m = zeros(DFloat, max(3, nmoist))

  @inbounds for e in elems
    for k = 1:Nq, j = 1:Nq, i = 1:Nq
      MJ = vgeo[i, j, k, _MJ, e]
      MJI = vgeo[i, j, k, _MJI, 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]
      E_int = E - (U^2 + V^2 + W^2)/(2*ρ) - ρ * gravity * z

      for m = 1:nmoist
          s = _nstate + m
          q_m[m] = Q[i, j, k, s, e] / ρ
      end

      # Calculate temperpature after saturation adjustment
      # Required for phase partitioning to get q_liq, q_ice
      T = saturation_adjustment(E_int/ρ, ρ, q_m[1])

      q = PhasePartition_equil(T, ρ, q_m[1])
      P = air_pressure(T, ρ, q)

      ρx, ρy, ρz = grad[i,j,k,_ρx,e], grad[i,j,k,_ρy,e], grad[i,j,k,_ρz,e]
      ux, uy, uz = grad[i,j,k,_ux,e], grad[i,j,k,_uy,e], grad[i,j,k,_uz,e]
      vx, vy, vz = grad[i,j,k,_vx,e], grad[i,j,k,_vy,e], grad[i,j,k,_vz,e]
      wx, wy, wz = grad[i,j,k,_wx,e], grad[i,j,k,_wy,e], grad[i,j,k,_wz,e]
      Tx, Ty, Tz = grad[i,j,k,_Tx,e], grad[i,j,k,_Ty,e], grad[i,j,k,_Tz,e]

      ρinv = 1 / ρ

      ldivu = λ_stokes*(ux + vy + wz)
      u = ρinv*U
      v = ρinv*V
      w = ρinv*W

      vfluxρ_x = 0*ρx
      vfluxρ_y = 0*ρy
      vfluxρ_z = 0*ρz

      vfluxU_x = 2*ux + ldivu
      vfluxU_y = uy + vx
      vfluxU_z = uz + wx

      vfluxV_x = vx + uy
      vfluxV_y = 2*vy + ldivu
      vfluxV_z = vz + wy

      vfluxW_x = wx + uz
      vfluxW_y = wy + vz
      vfluxW_z = 2*wz + ldivu

      vfluxE_x = u*(2*ux + ldivu) + v*(uy + vx) + w*(uz + wx) + k_μ*Tx
      vfluxE_y = u*(vx + uy) + v*(2*vy + ldivu) + w*(vz + wy) + k_μ*Ty
      vfluxE_z = u*(wx + uz) + v*(wy + vz) + w*(2*wz + ldivu) + k_μ*Tz

      fluxρ_x = U           - viscosity*vfluxρ_x
      fluxU_x = u * U  + P  - viscosity*vfluxU_x
      fluxV_x = u * V       - viscosity*vfluxV_x
      fluxW_x = u * W       - viscosity*vfluxW_x
      fluxE_x = u * (E + P) - viscosity*vfluxE_x

      fluxρ_y = V           - viscosity*vfluxρ_y
      fluxU_y = v * U       - viscosity*vfluxU_y
      fluxV_y = v * V + P   - viscosity*vfluxV_y
      fluxW_y = v * W       - viscosity*vfluxW_y
      fluxE_y = v * (E + P) - viscosity*vfluxE_y

      fluxρ_z = W           - viscosity*vfluxρ_z
      fluxU_z = w * U       - viscosity*vfluxU_z
      fluxV_z = w * V       - viscosity*vfluxV_z
      fluxW_z = w * W + P   - viscosity*vfluxW_z
      fluxE_z = w * (E + P) - viscosity*vfluxE_z

      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)

      # buoyancy term
      rhs[i, j, k, _W, e] -= ρ * gravity

      # Store velocity
      l_u[i, j, k], l_v[i, j, k], l_w[i, j, k] = u, v, w
    end

    # loop of ξ-grid lines
    for s = 1:_nstate, k = 1:Nq, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, k, _MJI, e]
      for n = 1:Nq
        rhs[i, j, k, s, e] += MJI * D[n, i] * s_F[n, j, k, s]
      end
    end
    # loop of η-grid lines
    for s = 1:_nstate, k = 1:Nq, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, k, _MJI, e]
      for n = 1:Nq
        rhs[i, j, k, s, e] += MJI * D[n, j] * s_G[i, n, k, s]
      end
    end
    # loop of ζ-grid lines
    for s = 1:_nstate, k = 1:Nq, j = 1:Nq, i = 1:Nq
      MJI = vgeo[i, j, k, _MJI, e]
      for n = 1:Nq
        rhs[i, j, k, s, e] += MJI * D[n, k] * s_H[i, j, n, s]
      end
    end

    # loop over moist variables
    # FIXME: Currently just passive advection
    for m = 1:nmoist
      s = _nstate + m
      ss = _nstategrad + 3*(m-1)

      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        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]
        u, v, w = l_u[i, j, k], l_v[i, j, k], l_w[i, j, k]

        q = Q[i, j, k, s, e]
        qx = grad[i, j, k, ss+1, e]
        qy = grad[i, j, k, ss+2, e]
        qz = grad[i, j, k, ss+3, e]

        fx = u*q - viscosity*0*qx
        fy = v*q - viscosity*0*qy
        fz = w*q - viscosity*0*qz

        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
      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nq
          rhs[i, j, k, s, e] += MJI * D[n, i] * s_F[n, j, k, 1]
        end
      end
      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nq
          rhs[i, j, k, s, e] += MJI * D[n, j] * s_G[i, n, k, 1]
        end
      end
      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nq
          rhs[i, j, k, s, e] += MJI * D[n, k] * s_H[i, j, n, 1]
        end
      end
    end

    # loop over tracer variables
    for t = 1:ntrace
      s = _nstate + nmoist + t

      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        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]
        u, v, w = l_u[i, j, k], l_v[i, j, k], l_w[i, j, k]

        fx = u * Q[i, j, k, s, e]
        fy = v * Q[i, j, k, s, e]
        fz = w * Q[i, j, k, s, e]

        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
      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nq
          rhs[i, j, k, s, e] += MJI * D[n, i] * s_F[n, j, k, 1]
        end
      end
      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nq
          rhs[i, j, k, s, e] += MJI * D[n, j] * s_G[i, n, k, 1]
        end
      end
      for k = 1:Nq, j = 1:Nq, i = 1:Nq
        MJI = vgeo[i, j, k, _MJI, e]
        for n = 1:Nq
          rhs[i, j, k, s, e] += MJI * D[n, k] * s_H[i, j, n, 1]
        end
      end
    end
  end
end
# }}}

# {{{ Face RHS (all dimensions)
function facerhs!(::Val{dim}, ::Val{N}, ::Val{nmoist}, ::Val{ntrace},
                  rhs::Array, Q, grad, vgeo, sgeo, gravity, viscosity,
                  elems, vmapM, vmapP, elemtobndy) where {N, dim, nmoist,
                                                          ntrace}
  DFloat = eltype(Q)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

    q_mM = zeros(DFloat, max(3, nmoist))
    q_mP = zeros(DFloat, max(3, nmoist))

  @inbounds for e in elems
    for f = 1:nface
      for n = 1:Nfp
        (nxM, nyM, nzM, sMJ, vMJI) = sgeo[:, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        ρxM = grad[vidM, _ρx, eM]
        ρyM = grad[vidM, _ρy, eM]
        ρzM = grad[vidM, _ρz, eM]
        uxM = grad[vidM, _ux, eM]
        uyM = grad[vidM, _uy, eM]
        uzM = grad[vidM, _uz, eM]
        vxM = grad[vidM, _vx, eM]
        vyM = grad[vidM, _vy, eM]
        vzM = grad[vidM, _vz, eM]
        wxM = grad[vidM, _wx, eM]
        wyM = grad[vidM, _wy, eM]
        wzM = grad[vidM, _wz, eM]
        TxM = grad[vidM, _Tx, eM]
        TyM = grad[vidM, _Ty, eM]
        TzM = grad[vidM, _Tz, eM]

        ρM = Q[vidM, _ρ, eM]
        UM = Q[vidM, _U, eM]
        VM = Q[vidM, _V, eM]
        WM = Q[vidM, _W, eM]
        EM = Q[vidM, _E, eM]
        yorzM = (dim == 2) ? vgeo[vidM, _y, eM] : vgeo[vidM, _z, eM]

        ρMinv = 1 / ρM
        uM, vM, wM = ρMinv * UM, ρMinv * VM, ρMinv * WM
        bc = elemtobndy[f, e]

        for m = 1:nmoist
            s = _nstate + m
            q_mM[m] = Q[vidM, s, eM] / ρM
        end

        E_intM = EM - (UM^2 + VM^2+ WM^2)/(2*ρM) - ρM * gravity * yorzM

        # get adjusted temperature and liquid and ice specific humidities
        TM = saturation_adjustment(E_intM/ρM , ρM, q_mM[1])
        qM = PhasePartition_equil(TM, ρM, q_mM[1])
        PM = air_pressure(TM, ρM, qM)

        if bc == 0

          ρP = Q[vidP, _ρ, eP]
          UP = Q[vidP, _U, eP]
          VP = Q[vidP, _V, eP]
          WP = Q[vidP, _W, eP]
          EP = Q[vidP, _E, eP]

          for m = 1:nmoist
            s = _nstate + m
            q_mP[m] = Q[vidP, s, eP] / ρP
          end

          yorzP = (dim == 2) ? vgeo[vidP, _y, eP] : vgeo[vidP, _z, eP]
          E_intP= EP - (UP^2 + VP^2+ WP^2)/(2*ρP) - ρP * gravity * yorzP

          TP = saturation_adjustment(E_intP/ρP, ρP, q_mP[1])
          qP = PhasePartition_equil(TP, ρP, q_mP[1])
          PP = air_pressure(TP, ρP, qP)

          ρxP = grad[vidP, _ρx, eP]
          ρyP = grad[vidP, _ρy, eP]
          ρzP = grad[vidP, _ρz, eP]
          uxP = grad[vidP, _ux, eP]
          uyP = grad[vidP, _uy, eP]
          uzP = grad[vidP, _uz, eP]
          vxP = grad[vidP, _vx, eP]
          vyP = grad[vidP, _vy, eP]
          vzP = grad[vidP, _vz, eP]
          wxP = grad[vidP, _wx, eP]
          wyP = grad[vidP, _wy, eP]
          wzP = grad[vidP, _wz, eP]
          TxP = grad[vidP, _Tx, eP]
          TyP = grad[vidP, _Ty, eP]
          TzP = grad[vidP, _Tz, eP]
        elseif bc == 1
          UnM = nxM * UM + nyM * VM + nzM * WM
          UP = UM - 2 * UnM * nxM
          VP = VM - 2 * UnM * nyM
          WP = WM - 2 * UnM * nzM
          ρP = ρM
          EP = EM
          PP = PM
          TP = TM

          ρnM = nxM * ρxM + nyM * ρyM + nzM * ρzM
          ρxP = ρxM - 2 * ρnM * nxM
          ρyP = ρyM - 2 * ρnM * nyM
          ρzP = ρzM - 2 * ρnM * nzM

          unM = nxM * uxM + nyM * uyM + nzM * uzM
          uxP = uxM - 2 * unM * nxM
          uyP = uyM - 2 * unM * nyM
          uzP = uzM - 2 * unM * nzM

          vnM = nxM * vxM + nyM * vyM + nzM * vzM
          vxP = vxM - 2 * vnM * nxM
          vyP = vyM - 2 * vnM * nyM
          vzP = vzM - 2 * vnM * nzM

          wnM = nxM * wxM + nyM * wyM + nzM * wzM
          wxP = wxM - 2 * wnM * nxM
          wyP = wyM - 2 * wnM * nyM
          wzP = wzM - 2 * wnM * nzM

          TxP, TyP, TzP = TxM, TyM, TzM
        else
          error("Invalid boundary conditions $bc on face $f of element $e")
        end

        ρPinv = 1 / ρP
        uP, vP, wP = ρPinv * UP, ρPinv * VP, ρPinv * WP

        fluxρM_x = UM
        fluxUM_x = uM * UM + PM
        fluxVM_x = uM * VM
        fluxWM_x = uM * WM
        fluxEM_x = uM * (EM + PM)

        fluxρM_y = VM
        fluxUM_y = vM * UM
        fluxVM_y = vM * VM + PM
        fluxWM_y = vM * WM
        fluxEM_y = vM * (EM + PM)

        fluxρM_z = WM
        fluxUM_z = wM * UM
        fluxVM_z = wM * VM
        fluxWM_z = wM * WM + PM
        fluxEM_z = wM * (EM + PM)

        fluxρP_x = UP
        fluxUP_x = uP * UP + PP
        fluxVP_x = uP * VP
        fluxWP_x = uP * WP
        fluxEP_x = uP * (EP + PP)

        fluxρP_y = VP
        fluxUP_y = vP * UP
        fluxVP_y = vP * VP + PP
        fluxWP_y = vP * WP
        fluxEP_y = vP * (EP + PP)

        fluxρP_z = WP
        fluxUP_z = wP * UP
        fluxVP_z = wP * VP
        fluxWP_z = wP * WP + PP
        fluxEP_z = wP * (EP + PP)

        λM = abs(nxM * uM + nyM * vM + nzM * wM) + soundspeed_air(TM)
        λP = abs(nxM * uP + nyM * vP + nzM * wP) + soundspeed_air(TP)
        λ  = max(λM, λP)

        #Compute Numerical Flux
        fluxρS = (nxM * (fluxρM_x + fluxρP_x) + nyM * (fluxρM_y + fluxρP_y) +
                  nzM * (fluxρM_z + fluxρP_z) - λ * (ρP - ρM)) / 2
        fluxUS = (nxM * (fluxUM_x + fluxUP_x) + nyM * (fluxUM_y + fluxUP_y) +
                  nzM * (fluxUM_z + fluxUP_z) - λ * (UP - UM)) / 2
        fluxVS = (nxM * (fluxVM_x + fluxVP_x) + nyM * (fluxVM_y + fluxVP_y) +
                  nzM * (fluxVM_z + fluxVP_z) - λ * (VP - VM)) / 2
        fluxWS = (nxM * (fluxWM_x + fluxWP_x) + nyM * (fluxWM_y + fluxWP_y) +
                  nzM * (fluxWM_z + fluxWP_z) - λ * (WP - WM)) / 2
        fluxES = (nxM * (fluxEM_x + fluxEP_x) + nyM * (fluxEM_y + fluxEP_y) +
                  nzM * (fluxEM_z + fluxEP_z) - λ * (EP - EM)) / 2

        #Compute Viscous Numerical Flux
        #Left Fluxes
        ldivuM = λ_stokes*(uxM + vyM + wzM)
        vfρM_x = 0*ρxM
        vfρM_y = 0*ρyM
        vfρM_z = 0*ρzM
        vfUM_x = 2*uxM + ldivuM
        vfUM_y = uyM + vxM
        vfUM_z = uzM + wxM
        vfVM_x = vxM + uyM
        vfVM_y = 2*vyM + ldivuM
        vfVM_z = vzM + wyM
        vfWM_x = wxM + uzM
        vfWM_y = wyM + vzM
        vfWM_z = 2*wzM + ldivuM
        vfEM_x = uM*(2*uxM + ldivuM) + vM*(uyM + vxM) + wM*(uzM + wxM) + k_μ*TxM
        vfEM_y = uM*(vxM + uyM) + vM*(2*vyM + ldivuM) + wM*(vzM + wyM) + k_μ*TyM
        vfEM_z = uM*(wxM + uzM) + vM*(wyM + vzM) + wM*(2*wzM + ldivuM) + k_μ*TzM

        #Right Fluxes
        ldivuP = λ_stokes*(uxP + vyP + wzP)
        vfρP_x = 0*ρxP
        vfρP_y = 0*ρyP
        vfρP_z = 0*ρzP
        vfUP_x = 2*uxP + ldivuP
        vfUP_y = uyP + vxP
        vfUP_z = uzP + wxP
        vfVP_x = vxP + uyP
        vfVP_y = 2*vyP + ldivuP
        vfVP_z = vzP + wyP
        vfWP_x = wxP + uzP
        vfWP_y = wyP + vzP
        vfWP_z = 2*wzP + ldivuP
        vfEP_x = uP*(2*uxP + ldivuP) + vP*(uyP + vxP) + wP*(uzP + wxP) + k_μ*TxP
        vfEP_y = uP*(vxP + uyP) + vP*(2*vyP + ldivuP) + wP*(vzP + wyP) + k_μ*TyP
        vfEP_z = uP*(wxP + uzP) + vP*(wyP + vzP) + wP*(2*wzP + ldivuP) + k_μ*TzP

        #Compute Numerical Flux
        vfluxρS = (nxM*(vfρM_x + vfρP_x) + nyM*(vfρM_y + vfρP_y) +
                   nzM*(vfρM_z + vfρP_z))/2
        vfluxUS = (nxM*(vfUM_x + vfUP_x) + nyM*(vfUM_y + vfUP_y) +
                   nzM*(vfUM_z + vfUP_z))/2
        vfluxVS = (nxM*(vfVM_x + vfVP_x) + nyM*(vfVM_y + vfVP_y) +
                   nzM*(vfVM_z + vfVP_z))/2
        vfluxWS = (nxM*(vfWM_x + vfWP_x) + nyM*(vfWM_y + vfWP_y) +
                   nzM*(vfWM_z + vfWP_z))/2
        vfluxES = (nxM*(vfEM_x + vfEP_x) + nyM*(vfEM_y + vfEP_y) +
                   nzM*(vfEM_z + vfEP_z))/2

        #Update RHS
        rhs[vidM, _ρ, eM] -= vMJI * sMJ * (fluxρS - viscosity*vfluxρS)
        rhs[vidM, _U, eM] -= vMJI * sMJ * (fluxUS - viscosity*vfluxUS)
        rhs[vidM, _V, eM] -= vMJI * sMJ * (fluxVS - viscosity*vfluxVS)
        rhs[vidM, _W, eM] -= vMJI * sMJ * (fluxWS - viscosity*vfluxWS)
        rhs[vidM, _E, eM] -= vMJI * sMJ * (fluxES - viscosity*vfluxES)

        # loop over moist variables
        # FIXME: Currently just passive advection
        for m = 1:nmoist
          s = _nstate + m
          ss = _nstategrad + 3*(m-1)

          QmoistM = Q[vidM, s, eM]

          qxM = grad[vidM, ss + 1, eM]
          qyM = grad[vidM, ss + 2, eM]
          qzM = grad[vidM, ss + 3, eM]

          if bc == 0
            QmoistP = Q[vidP, s, eP]
            qxP = grad[vidP, ss + 1, eP]
            qyP = grad[vidP, ss + 2, eP]
            qzP = grad[vidP, ss + 3, eP]
          elseif bc == 1
            QmoistP = QmoistM
            qnM = nxM * qxM + nyM * qyM + nzM * qzM
            qxP = qxM - 2 * qnM * nxM
            qyP = qyM - 2 * qnM * nyM
            qzP = qzM - 2 * qnM * nzM
          else
            error("Invalid boundary conditions $bc on face $f of element $e")
          end

          fluxM_x, fluxP_x = uM * QmoistM, uP * QmoistP
          fluxM_y, fluxP_y = vM * QmoistM, vP * QmoistP
          fluxM_z, fluxP_z = wM * QmoistM, wP * QmoistP

          fluxS = (nxM * (fluxM_x + fluxP_x) + nyM * (fluxM_y + fluxP_y) +
                   nzM * (fluxM_z + fluxP_z) - λ * (QmoistP - QmoistM)) / 2

          # Moist Variable viscous flux
          vfqM_x, vfqP_x = 0*qxM, 0*qxP
          vfqM_y, vfqP_y = 0*qyM, 0*qyP
          vfqM_z, vfqP_z = 0*qzM, 0*qzP
          vfluxqS = (nxM*(vfqM_x + vfqP_x) + nyM*(vfqM_y + vfqP_y) +
                     nzM*(vfqM_z + vfqP_z))/2

          rhs[vidM, s, eM] -= vMJI * sMJ * (fluxS - viscosity*vfluxqS)
        end

        # FIXME: Will need to be updated for other bcs...
        vidP = bc == 0 ? vidP : vidM

        # loop over tracer variables
        for t = 1:ntrace
          s = _nstate + nmoist + t
          QtraceM, QtraceP = Q[vidM, s, eM], Q[vidP, s, eP]

          fluxM_x, fluxP_x = uM * QtraceM, uP * QtraceP
          fluxM_y, fluxP_y = vM * QtraceM, vP * QtraceP
          fluxM_z, fluxP_z = wM * QtraceM, wP * QtraceP

          fluxS = (nxM * (fluxM_x + fluxP_x) + nyM * (fluxM_y + fluxP_y) +
                   nzM * (fluxM_z + fluxP_z) - λ * (QtraceP - QtraceM)) / 2

          rhs[vidM, s, eM] -= vMJI * sMJ * fluxS
        end
      end
    end
  end
end
# }}}

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
#!/usr/bin/env julia
using NeRCA
using Plots
using PlotThemes
using Dates
using Measures
GR.inline("png")

if length(ARGS) < 2
    println("Usage: ./live_royfit.jl DETX TIME_RES")
    exit(1)
end

const calib = NeRCA.read_calibration(ARGS[1])
const LIGIER_PORT = parse(Int, ARGS[2])
const TIME_RES = ARGS[3]

function main()
    println("Starting live ROyFit")

    for message in CHClient(ip"127.0.0.1", LIGIER_PORT, ["IO_EVT"])
        event = NeRCA.read_io(IOBuffer(message.data), NeRCA.DAQEvent)

        hits = calibrate(event.hits, calib)
        triggered_hits = filter(h->h.triggered, hits)
        dus = sort(unique(map(h->h.du, hits)))
        n_dus = length(dus)
        n_doms = length(unique(h->h.dom_id, triggered_hits))

        if n_doms < 4
            continue
        end

        colours = palette(:default)
        plot()
        Q = []
        for (idx, du) in enumerate(dus)
            du_hits = filter(h->h.du == du, hits)
            if length(triggered(du_hits))== 0
                println("No triggered hits")
                continue
            end
            fit = NeRCA.single_du_fit(du_hits)
            push!(Q, fit.Q)
            plot!(du_hits, fit, markercolor=colours[idx], label="DU $(du)", max_z=calib.max_z)
            write_time_residuals(TIME_RES, event, du_hits, fit)
        end
        if sum(Q) < 200 && n_doms > 12 && n_dus > 1
            println("Plotting...")
            fit_params = "ROy live reconstruction (combined single line): Q=$([round(_Q,digits=2) for _Q in Q])"
            event_params = "Det ID $(event.det_id), Run $(event.run_id), FrameIndex $(event.timeslice_id), TriggerCounter $(event.trigger_counter), Overlays $(event.overlays)"
            time_params = "$(unix2datetime(event.timestamp)) UTC"
            trigger_params = "Trigger: $(is_mxshower(event) ? "MX " : "")$(is_3dmuon(event) ? "3DM " : "")$(is_3dshower(event) ? "3DS " : "")"
            time_params = "$(unix2datetime(event.timestamp)) UTC"

            plot!(title="$(fit_params)\n$(event_params), $(trigger_params)\n$(time_params)", titlefontsize=8, margin=5mm)

            savefig("plots/ztplot_roy.png")
        end
    end
end

function write_time_residuals(filename, event, hits, fit)
    if !isfile(filename)
        fobj = open(filename, "w")
        write(fobj, "run,timestamp,du,floor,dom_id,t_res,Q\n")
    else
        fobj = open(filename, "a")
    end
    dγ, ccalc = NeRCA.make_cherenkov_calculator(fit.sdp)
    for hit in hits
        Δt = hit.t - ccalc(hit.pos.z)
        write(fobj, "$(event.run_id),$(event.timestamp),$(hit.du),$(hit.floor),$(hit.dom_id),$(Δt),$(fit.Q)\n")
    end
    close(fobj)
end

main()


text:
"""
Free Variable

Like random variables, free variables are uncertain values.
Unlike random variables the uncertainty is not quantified,
no more value is more or less likely than any other.

"""
abstract type Var <: NonDetVar end

# Primitive Free Variables
struct Unit{T <: Real} <: Var
  id::ID
end

"`unit(T)` is a random variable over [0, 1]"
unit(::Type{T} = Float64; id = uid()) where T = Unit{T}(id)

solve(::Unit{T}) where T = rand(T)

"Find a solution to all the variables in the model that satisfies"
function solve(x::Var) end

"Find a solution that maps the variable values to free"
function optim(x::Var, ℓ) end


text:
# helper functions for vector inputs

@inline isweekend(dt::Vector{Dates.Date}) = isweekend.(dt)

function isbday(hc::HolidayCalendar, dt::Vector{Dates.Date})
    result = Vector{Bool}(undef, length(dt))

    for i in eachindex(dt)
        @inbounds result[i] = isbday(hc, dt[i])
    end

    return result
end

function isbday(hc::Vector{HolidayCalendar}, dt::Vector{Dates.Date})
    l_hc = length(hc)
    l_dt = length(dt)

    @assert l_hc == l_dt "Input vectors must have the same size. $(l_hc) != $(l_dt)"

    result = Vector{Bool}(undef, length(dt))

    for i in 1:l_hc
        @inbounds result[i] = isbday(hc[i], dt[i])
    end

    return result
end

isbday(calendar, dt::Vector{Dates.Date}) = isbday(convert(HolidayCalendar, calendar), dt)
isbday(calendars::A, dt::Vector{Dates.Date}) where {A<:AbstractArray} = isbday(convert(Vector{HolidayCalendar}, calendars), dt)

function tobday(hc::HolidayCalendar, dt::Vector{Dates.Date}; forward::Bool = true)
    result = Vector{Dates.Date}(undef, length(dt))

    for i in eachindex(dt)
        @inbounds result[i] = tobday(hc, dt[i]; forward=forward)
    end

    return result
end

function tobday(hc::Vector{HolidayCalendar}, dt::Vector{Dates.Date}; forward::Bool = true)
    l_hc = length(hc)
    l_dt = length(dt)

    @assert l_hc == l_dt "Input vectors must have the same size. $(l_hc) != $(l_dt)"

    result = Vector{Dates.Date}(undef, l_hc)

    for i in 1:l_hc
        @inbounds result[i] = tobday(hc[i], dt[i]; forward=forward)
    end

    return result
end

tobday(calendar, dt::Vector{Dates.Date}; forward::Bool = true) = tobday(convert(HolidayCalendar, calendar), dt; forward=forward)
tobday(calendars::A, dt::Vector{Dates.Date}; forward::Bool = true) where {A<:AbstractArray} = tobday(convert(Vector{HolidayCalendar}, calendars), dt; forward=forward)

function bdays(hc::HolidayCalendar, base_date::Dates.Date, dt_vec::Vector{Dates.Date})
    len = length(dt_vec)
    result = Vector{Dates.Day}(undef, len)

    for i in 1:len
        @inbounds result[i] = bdays(hc, base_date, dt_vec[i])
    end

    return result
end

function bdayscount(hc::HolidayCalendar, base_date::Dates.Date, dt_vec::Vector{Dates.Date})
    len = length(dt_vec)
    result = Vector{Int}(undef, len)

    for i in 1:len
        @inbounds result[i] = bdayscount(hc, base_date, dt_vec[i])
    end

    return result
end

function bdays(hc::HolidayCalendar, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date})
    l0 = length(dt0)
    l1 = length(dt1)

    @assert l0 == l1 "Input vectors must have the same size. $(l0) != $(l1)"

    result = Vector{Dates.Day}(undef, l0)

    for i in 1:l0
        @inbounds result[i] = bdays(hc, dt0[i], dt1[i])
    end

    return result
end

function bdayscount(hc::HolidayCalendar, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date})
    l0 = length(dt0)
    l1 = length(dt1)

    @assert l0 == l1 "Input vectors must have the same size. $(l0) != $(l1)"

    result = Vector{Int}(undef, l0)

    for i in 1:l0
        @inbounds result[i] = bdayscount(hc, dt0[i], dt1[i])
    end

    return result
end

function bdays(hc::Vector{HolidayCalendar}, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date})
    l_hc = length(hc)
    l0 = length(dt0)
    l1 = length(dt1)

    @assert l_hc == l0 && l0 == l1 "Input vectors must have the same size. $(l_hc), $(l0), $(l1)"

    result = Vector{Dates.Day}(undef, l0)

    for i in 1:l0
        @inbounds result[i] = bdays(hc[i], dt0[i], dt1[i])
    end

    return result
end

function bdayscount(hc::Vector{HolidayCalendar}, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date})
    l_hc = length(hc)
    l0 = length(dt0)
    l1 = length(dt1)

    @assert l_hc == l0 && l0 == l1 "Input vectors must have the same size. $(l_hc), $(l0), $(l1)"

    result = Vector{Int}(undef, l0)

    for i in 1:l0
        @inbounds result[i] = bdayscount(hc[i], dt0[i], dt1[i])
    end

    return result
end

bdays(calendars::A, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date}) where {A<:AbstractArray} = bdays(convert(Vector{HolidayCalendar}, calendars), dt0, dt1)
bdayscount(calendars::A, dt0::Vector{Dates.Date}, dt1::Vector{Dates.Date}) where {A<:AbstractArray} = bdayscount(convert(Vector{HolidayCalendar}, calendars), dt0, dt1)

function advancebdays(hc::HolidayCalendar, dt::Dates.Date, bdays_count_vec::Vector{T}) where {T<:Union{Int, Dates.Day}}
    l = length(bdays_count_vec)
    result = Vector{Dates.Date}(undef, l)
    for i in 1:l
        @inbounds result[i] = advancebdays(hc, dt, bdays_count_vec[i])
    end
    return result
end

function advancebdays(calendar, dt::Dates.Date, bdays_count_vec::Vector{T}) where {T<:Union{Int, Dates.Day}}
    advancebdays(convert(HolidayCalendar, calendar), dt, bdays_count_vec)
end

advancebdays(hc, dt::Dates.Date, bdays_range::AbstractRange) = advancebdays(hc, dt, collect(bdays_range))


text:
module URLParse

warn("""
The URLParse.jl package has been deprecated in favour of JuliaWeb/URIParser.jl
https://github.com/JuliaWeb/URIParser.jl
As of Julia 0.4 this package will no longer be installable
through `Pkg.add`. Please convert your code accordingly.""")

import Base.copy

export URLComponents,
        urlsplit, urlunsplit, urlparse, urlunparse, urldefrag,
        username, password, hostname, port,
        copy,
        escape, escape_form, escape_with, unescape, unescape_form

#include("strutils.jl")
include("urlescape.jl")

# A classification of schemes ('' means apply by default)
const uses_relative = ["hdfs", "ftp", "http", "gopher", "nntp", "imap", "wais", "file", "https", "shttp", "mms", "prospero", "rtsp", "rtspu", "", "sftp", "svn", "svn+ssh"]
const uses_netloc = ["hdfs", "ftp", "http", "gopher", "nntp", "telnet", "imap", "wais", "file", "mms", "https", "shttp", "snews", "prospero", "rtsp", "rtspu", "rsync", "", "svn", "svn+ssh", "sftp" ,"nfs", "git", "git+ssh"]
const uses_params = ["ftp", "hdl", "prospero", "http", "imap", "https", "shttp", "rtsp", "rtspu", "sip", "sips", "mms", "", "sftp", "tel"]

# These are not actually used anymore, but should stay for backwards
# compatibility.  (They are undocumented, but have a public-looking name.)
const non_hierarchical = ["gopher", "hdl", "mailto", "news", "telnet", "wais", "imap", "snews", "sip", "sips"]
const uses_query = ["http", "wais", "imap", "https", "shttp", "mms", "gopher", "rtsp", "rtspu", "sip", "sips", ""]
const uses_fragment = ["hdfs", "ftp", "hdl", "http", "gopher", "news", "nntp", "wais", "https", "shttp", "snews", "file", "prospero", ""]
const scheme_chars = "abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789+-."

type URLComponents
    scheme::String
    netloc::String
    url::String
    params::String
    query::String
    fragment::String

    _username::Union(String, Bool, Nothing)
    _password::Union(String, Bool, Nothing)
    _hostname::Union(String, Bool, Nothing)
    _port::Union(Int, Bool, Nothing)

    function URLComponents(scheme::String, netloc::String, url::String, params::String, query::String, fragment::String)
        new(scheme, netloc, url, params, query, fragment, nothing, nothing, nothing, nothing)
    end
end

const MAX_CACHE_SIZE = 20
const _parse_cache = Dict{(String, String, Bool),URLComponents}()

clear_cache() = empty!(_parse_cache)
copy(up::URLComponents) = URLComponents(up.scheme, up.netloc, up.url, up.params, up.query, up.fragment)

function _parse_user_name_password(up::URLComponents)
    netloc = up.netloc
    if('@' in netloc)
        userinfo = rsplit(netloc, '@', 2)[1]
        if(':' in userinfo)
            spl = split(userinfo, ':', 2)
            up._username = spl[1]
            up._password = (length(spl) > 1) ? spl[2] : false
        else
            up._username = userinfo
            up._password = false
        end
    else
        up._username = up._password = false
    end
end

function _parse_hostname_port(up::URLComponents)
    netloc = rsplit(up.netloc, '@', 2)[end]
    lb = search(netloc, '[')
    rb = search(netloc, ']')
    cpos = rsearch(netloc, ':')
    if((lb > 0) && (rb > 0))
        up._hostname = lowercase(netloc[lb+1:rb-1])
    elseif(cpos > 0)
        up._hostname = lowercase(netloc[1:cpos-1])
    elseif(length(netloc) > 0)
        up._hostname = lowercase(netloc)
    else
        up._hostname = up._port = false
    end

    if((false != up._port) && (cpos > 0) && ((rb == 0) || (cpos > rb)))
        up._port = int(netloc[cpos+1:end])
    else
        up._port = false
    end
end

# returns the username string or nothing if it was not specified in the url
function username(up::URLComponents)
    (nothing != up._username) && (return (false == up._username) ? nothing : up._username)
    _parse_user_name_password(up)
    return (false == up._username) ? nothing : up._username
end

# returns the password string or nothing if it was not specified in the url
function password(up::URLComponents)
    (nothing != up._username) && (return (false == up._password) ? nothing : up._password)
    _parse_user_name_password(up)
    return (false == up._password) ? nothing : up._password
end

# returns the hostname string or nothing if it was not specified in the url
function hostname(up::URLComponents)
    (nothing != up._hostname) && (return (false == up._hostname) ? nothing: up._hostname)
    _parse_hostname_port(up)
    return (false == up._hostname) ? nothing : up._hostname
end

# returns the port (int) or nothing if it was not specified in the url
function port(up::URLComponents)
    (nothing != up._port) && (return (false == up._port) ? nothing: up._port)
    _parse_hostname_port(up)
    return (false == up._port) ? nothing : up._port
end

function _splitparams(url)
    slash_pos = rsearch(url, '/')
    sc_pos = 0
    if(slash_pos > 0)
        sc_pos = search(url, ';', slash_pos)
        (sc_pos < 1) && (return (url, ""))
    else
        sc_pos = search(url, ';')
    end
    (sc_pos < 1) ? (url, url) : (url[1:sc_pos-1], url[sc_pos+1:end])
end

function _splitnetloc(url::String, start::Int=1)
    delim = length(url)+1                     # position of end of domain part of url, default is end
    for c in "/?#"                            # look for delimiters; the order is NOT important
        wdelim = search(url, c, start)        # find first of this delim
        if (wdelim > 0)                        # if found
            delim = min(delim, wdelim)        # use earliest delim position
        end
    end
    return _chk_invalid_ipv6(url[start:delim-1]), url[delim:end]   # return (domain, rest)
end

function _chk_invalid_ipv6(netloc::String)
    lb = search(netloc, '[')
    rb = search(netloc, ']')
    if( ((lb > 0) && (rb == 0)) || ((lb == 0) && (rb > 0)))
        error("Invalid IPV6 URL")
    end
    netloc
end

##
# Parse a URL into 6 components:
# <scheme>://<netloc>/<path>;<params>?<query>#<fragment>
function urlparse(url::String, scheme::String="", allow_fragments::Bool=true)
    up = urlsplit(url, scheme, allow_fragments)
    if (up.scheme in uses_params) && (';' in url)
        up.url, up.params = _splitparams(up.url)
    end
    up
end

##
# Parse a URL into 5 components:
# <scheme>://<netloc>/<path>?<query>#<fragment>
function urlsplit(url::String, scheme::String="", allow_fragments::Bool=true)
    key = (url, scheme, allow_fragments)
    if(haskey(_parse_cache, key))
        return _parse_cache[key]
    end
    (length(_parse_cache) >= MAX_CACHE_SIZE) && empty!(_parse_cache)

    netloc = ""; query = ""; fragment = ""
    i = search(url, ':')
    if(i > 1)
        if(url[1:(i-1)] == "http")  # optimize the common case
            scheme = lowercase(url[1:(i-1)])
            url = url[i+1:end]
            if(beginswith(url, "//"))
                netloc, url = _splitnetloc(url, 3)
            end
            if(allow_fragments && ('#' in url))
                url, fragment = tuple(split(url, '#', 2)...)
            end
            if('?' in url)
                url, query = tuple(split(url, '?', 2)...)
            end
            return _parse_cache[key] = URLComponents(scheme, netloc, url, "", query, fragment)
        end
        urlinvalidchars = false
        for ci in 1:(i-1)
            !(url[ci] in scheme_chars) && (urlinvalidchars = true) && break
        end
        if(!urlinvalidchars)
            try
                _testportnum = int(url[i+1:end])
            catch
                scheme = lowercase(url[1:i-1])
                url = url[i+1:end]
            end
        end
    end
    if(beginswith(url, "//"))
        netloc, url = _splitnetloc(url, 3)
    end

    scheme_uses_fragment = (scheme in uses_fragment)
    scheme_uses_query = (scheme in uses_query)
    if (allow_fragments && scheme_uses_fragment && ('#' in url))
        url, fragment = tuple(split(url, '#', 2)...)
    end
    if (scheme_uses_query && ('?' in url))
        url, query = tuple(split(url, '?', 2)...)
    end
    return _parse_cache[key] = URLComponents(scheme, netloc, url, "", query, fragment)
end

urlunparse(up::URLComponents) = urlunsplit(up, true)
function urlunsplit(up::URLComponents, add_params::Bool=false)
    has_netloc = (length(up.netloc) > 0)
    has_scheme = (length(up.scheme) > 0)
    scheme_uses_netloc = has_scheme && (up.scheme in uses_netloc)
    url = (add_params && (length(up.params) > 0)) ? (up.url * ";" * up.params) : up.url
    if (has_netloc || (has_scheme && scheme_uses_netloc && !beginswith(url, "//")))
        if((length(url) > 0) && (url[1] != '/'))
            url = "/" * url
        end
        url = "//" * up.netloc * url
    end
    has_scheme && (url = up.scheme * ":" * url)
    (length(up.query) > 0) && (url = url * "?" * up.query)
    (length(up.fragment) > 0) && (url = url * "#" * up.fragment)
    url
end

##
# Removes any existing fragment from URL.
# Returns a tuple of the defragmented URL and the fragment.  If
#    the URL contained no fragments, the second element is the
#    empty string.
function urldefrag(url::String)
    if('#' in url)
        up = copy(urlparse(url))
        frag = up.fragment
        up.fragment = ""
        defrag = urlunparse(up)
        return defrag, frag
    end
    return url, ""
end

end


text:
# txt = readlines("test_day1.txt")

# txt = map(x->parse(Int64,x),txt)

# count = 0
# for i in 2:size(txt, 1)
#     if txt[i] > txt[i-1]
#         count += 1
#     end
# end

# print(count)

txt = map(x->parse(Int64,x),readlines("day1.txt"))
sum([1 for i in 2:size(txt, 1) if txt[i] > txt[i-1]])


text:
# julia1_2.1.jl Syntax and some simple I/O
#
# Shell Game example (Coconut shells and 1 pea)
# shells are variable names
# under each shell, type with value
# This example shows that variable names are just labels
# they can be made concrete and that will be shown later
# type and data go together as the following will illustrate
a=3
println("a is ",typeof(a))
b=[1 2; 3 4]
println("b is ",typeof(b))
c="pea"
println("c is ",typeof(c))
println("Before: $a, $b, $c")
# we can shuffle these around
a,b,c = b,c,a  #shell game
a,b,c = b,a,c  #where's the pea?
println("After: $a, $b, $c")

# How deep is a well? Use a stopwatch
# dropping a rock and listening for the splash
# depth=.5g*td^2 td=time to reach depth (s=.5a*t^2)
# depth=vs*tu    tu=time sound take from bottom to top
g=9.8065    # acceleration due to gravity
vs=343.0    # m/s velocity of sound
print("drop to Splash time: ") # seconds from drop to ear
st=readline();       #Arg STDIN assumed if left out
t=parse(Float32,st)  #string to Float32 conversion for t
# solve equations by
# EQ1: .5g*td^2 - vs*tu = 0
# EQ2: td+tu=t => tu=t-td subsitute in EQ1 and get
# quadratic .5g*td^2 + vs*td -vs*t = 0  in t1
td = (-vs + sqrt(vs^2 + 2.0g*vs*t))/g
tu = t - td
depthu = vs*tu
println("depth $depthu meters\n sound time $tu seconds")
# confirm
depthd = .5g*td^2
println("confirming depth $depthd meters\n fall time $td seconds")
# terminal velocity = g*td
termv = g*td
println("terminal velocity $termv meters/second")

# some simple files IO to read and write CSV text files.
fi1=open("str1.csv")
fo1=open("str2.csv","w")
cnt=1
while(!eof(fi1))
  global cnt
  row=readline(fi1)
  row="$cnt,"*row*"\n"
  write(fo1,row)
  print(row)
  cnt+=1
end
close(fi1)
close(fo1)
println("Done")


text: # import MNIST


text:
include("time.jl")
include("sleep.jl")

"""
`AbstractClock` is an abstract type for clocks

Clocks return instant (`DateTime`) when asked using `now(clock)`

Clock can be real (ie using system time) but can also be fake
(for simulation purpose).

Time can be set on a simulutated clock using `set(clock, new_datetime)`
"""
abstract type AbstractClock end

"""
An exception that a clock can throw.

Generally a clock run an exception when user is trying to set time
on a system clock (which is not allowed).
"""
struct ClockException <: Exception
    s
end

include("rtclock.jl")
include("simclock.jl")


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
import Base.show
using ViewerGL,SparseArrays
GL = ViewerGL

function show(filename)
	V, EV = Lar.svg2lar(filename)
	GL.VIEW([ GL.GLLines(V,EV) ])
	return V, EV
end

# show("./test/svg/new.svg")
# show("./test/svg/curved.svg")
# show("./test/svg/twopaths.svg")
# show("./test/svg/paths.svg")
# show("./test/svg/boundarytest2.svg")
# show("./test/svg/tile.svg")
# show("./test/svg/interior.svg")
V,EV = show("./test/svg/holes.svg")
# V,EV = show("./test/svg/Lar.svg")

# subdivision of input edges
W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV, copFE = Lar.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# visualization of component graphs
EW = Lar.cop2lar(copEV)
W = convert(Lar.Points, V')
V,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells)) # 2-connected components
VV = [[k] for k=1:size(V,2)]
meshes = [ GL.numbering(.1)((V, (VV,EVs[k])), GL.COLORS[k],1) for k=1:length(EVs) ]
GL.VIEW( cat(meshes) );

# final solid visualizations
FE = [SparseArrays.findnz(copFE[k,:])[1] for k=1:size(copFE,1)]
FV = [collect(Set(cat([EW[e] for e in FE[f]]))) for f=1:length(FE)]
V = convert(Lar.Points, V')
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
V = convert(Lar.Points, V')
FVs = convert(Array{Lar.Cells}, triangulated_faces)
GL.VIEW(GL.GLExplode(V,FVs,1,1,1,99,1));  # TODO:  fix inner holes

V,FVs,EVs = Lar.arrange2D(V,EV)
GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,99,1));
GL.VIEW(GL.GLExplode(V,FVs,1,1,1,99,1));  # TODO:  fix inner holes

EVs = Lar.FV2EVs(copEV, copFE) # polygonal face boundaries
EVs = convert(Array{Array{Array{Int64,1},1},1}, EVs)
GL.VIEW(GL.GLExplode(V,EVs,1,1,1,3,1));  # TODO:  fix inner holes


text:
module GMAT

using Pkg.Artifacts

using GMAT_jll

export run_script

include("startup_file.jl")

"""
	run_script(script;
		data=joinpath(artifact"data", "gmat-data-2020a"),
		output=joinpath(pwd(), "gmat-output"),
	)

Run a GMAT `script` while optionally specifying the path to the `data` package and an
`output` directory. The function will return a list of output files.

### References

- [GMAT Documentation](http://gmat.sourceforge.net/docs/R2020a/html/index.html)
"""
function run_script(script;
	data=joinpath(artifact"data", "gmat-data-2020a"),
	output=joinpath(pwd(), "gmat-output"),
)
	# Create output dir if it does not exist
	isdir(output) || mkpath(output)
	# Write the GMAT startup file to a temp file
	startup = tempname()
	open(startup, "w") do f
		write(f, startup_file(data, output))
	end
	gmatconsole() do exe
		run(`$exe -s $startup -r $script`)
	end
	return readdir(output)
end

end # module


text:
using Franklin
const F = Franklin

html(s) = "\n~~~$s~~~\n"

function hfun_adddescription()
    d = locvar(:description)
    isnothing(d) ? "" : F.fd2html(d, internal=true)
end

function hfun_frontmatter()
    fm = locvar(:front_matter)
    if isnothing(fm)
        ""
    else
        """
        <d-front-matter>

        </d-front-matter>
        """
    end
end

function hfun_byline()
    fm = locvar(:front_matter)
    if isnothing(fm)
        ""
    else
        "<d-byline></d-byline>"
    end
end

function hfun_dtoc()
    is_enable_toc = locvar(:is_enable_toc)
    minlevel = locvar("mintoclevel")
    maxlevel = locvar("maxtoclevel")
    toc = F.hfun_toc([string(minlevel), string(maxlevel)])
    if isnothing(is_enable_toc) || !is_enable_toc
        ""
    else
        """

        <d-article class="franklin-content">

Table of content

        $toc
        </d-article>

        """
    end
end

function hfun_appendix()
    ap = locvar(:appendix)
    ap = isnothing(ap) ? "" : F.md2html(ap)
    bib = locvar(:bibliography)
    if isnothing(bib)
        bib = ""
    else
        bib_in_cur_folder = joinpath(splitdir(locvar("fd_rpath"))[1], bib)
        if isfile(bib_in_cur_folder)
            bib_resolved = F.parse_rpath("/" * bib_in_cur_folder)
        else
            bib_resolved = F.parse_rpath(bib; canonical=false, code=true)
        end
        bib = "<d-bibliography src=\"$bib_resolved\"></d-bibliography>"
    end
    """
    <d-appendix>
        $ap
        $bib
    </d-appendix>
    """
end

function lx_dcite(lxc,_)
    content = F.content(lxc.braces[1])
    "<d-cite key=\"$content\"></d-cite>" |> html
end

function dfigure(layout, src, caption)
    """

$caption

    """ |> html
end

# https://github.com/distillpub/template/blob/b854bd0124911e1be4744e472b67832e3251b96c/src/styles/styles-layout.css#L137-L166
"""
Possible layouts:
"""
function lx_dfig(lxc,lxd)
    content = F.content(lxc.braces[1])
    info = split(content, ';')
    layout = info[1]
    src = info[2]
    caption = F.reprocess(join(info[3:end], ';'), lxd)

    # (case 1) it is a figure from network
    startswith(src, "http") && return dfigure(layout, src, caption)

    # (case 2) it's in current folder
    src_in_cur_folder = joinpath(splitdir(locvar("fd_rpath"))[1], src)
    if isfile(src_in_cur_folder)
        src_resolved = F.parse_rpath("/" * src_in_cur_folder)
        @info "resolving fig path" src src_resolved
        return dfigure(layout, src_resolved, caption)
    end

    # (case 3) assume it is generated by code
    src = F.parse_rpath(src; canonical=false, code=true)

    # !!! directly take from `lx_fig` in Franklin.jl
    fdir, fext = splitext(src)

    # there are several cases
    # A. a path with no extension --> guess extension
    # B. a path with extension --> use that
    # then in both cases there can be a relative path set but the user may mean
    # that it's in the subfolder /output/ (if generated by code) so should look
    # both in the relpath and if not found and if /output/ not already last dir
    candext = ifelse(isempty(fext),
                     (".png", ".jpeg", ".jpg", ".svg", ".gif"), (fext,))
    for ext ∈ candext
        candpath = fdir * ext
        if F.FD_ENV[:STRUCTURE] < v"0.2"
            syspath  = joinpath(F.PATHS[:folder], split(candpath, '/')...)
        else
            syspath  = joinpath(F.PATHS[:site], split(candpath, '/')...)
        end
        isfile(syspath) && return dfigure(layout, candpath, caption)
    end
    # now try in the output dir just in case (provided we weren't already
    # looking there)
    p1, p2 = splitdir(fdir)
    @debug "TEST" p1 p2
    if splitdir(p1)[2] != "output"
        for ext ∈ candext
            candpath = joinpath(splitdir(p1)[1], "output", p2 * ext)
            if F.FD_ENV[:STRUCTURE] < v"0.2"
                syspath  = joinpath(F.PATHS[:folder], split(candpath, '/')...)
            else
                syspath  = joinpath(F.PATHS[:site], split(candpath, '/')...)
            end
            isfile(syspath) && return dfigure(layout, candpath, caption)
        end
    end
end

function lx_aside(lxc,lxd)
    content = F.reprocess(F.content(lxc.braces[1]), lxd)
    "$content" |> html
end

function lx_footnote(lxc,lxd)
    content = F.reprocess(F.content(lxc.braces[1]), lxd)
    # workaround
    if startswith(content, "
")
        content = content[4:end-5]
    end
    "<d-footnote>$content</d-footnote>" |> html
end

function lx_appendix(lxc,lxd)
    content = F.reprocess(F.content(lxc.braces[1]), lxd)
    "<d-appendix>$content</d-appendix>" |> html
end


text:
function get_block(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    radius=0.0,
    side=0.5,
    contact=true,
    contact_type=:nonlinear,
    color=RGBA(0.0, 0.0, 0.0, 1.0),
    mode=:box,
    T=Float64)

    # Parameters
    origin = Origin{T}()
    pbody = Box(side, side, side, 1.0,
        color=color)
    joint0to1 = JointConstraint(Floating(origin, pbody))
    bodies = [pbody]
    joints = [joint0to1]

    if contact
        # Corner vectors
        if mode == :particle
            corners = [[0.0, 0.0, 0.0]]
        elseif mode == :box
            corners = [
                [[ side / 2.0;  side / 2.0; -side / 2.0]]
                [[ side / 2.0; -side / 2.0; -side / 2.0]]
                [[-side / 2.0;  side / 2.0; -side / 2.0]]
                [[-side / 2.0; -side / 2.0; -side / 2.0]]
                [[ side / 2.0;  side / 2.0;  side / 2.0]]
                [[ side / 2.0; -side / 2.0;  side / 2.0]]
                [[-side / 2.0;  side / 2.0;  side / 2.0]]
                [[-side / 2.0; -side / 2.0;  side / 2.0]]
            ]
        else
            @error "incorrect mode specified, try :particle or :box"
        end
        n = length(corners)
        normal = [[0.0, 0.0, 1.0] for i = 1:n]
        contact_radius = [radius for i = 1:n]
        friction_coefficient = friction_coefficient * ones(n)

        contacts = contact_constraint(pbody, normal,
            friction_coefficient=friction_coefficient,
            contact_origins=corners,
            contact_radius=contact_radius,
            contact_type=contact_type)

        mech = Mechanism(origin, bodies, joints, contacts,
            gravity=gravity,
            timestep=timestep)
    else
        mech = Mechanism(origin, bodies, joints,
            gravity=gravity,
            timestep=timestep)
    end
    return mech
end

function initialize_block!(mechanism::Mechanism{T};
        x=[0.0, 0.0, 1.0],
        q=Quaternion(1.0, 0.0, 0.0, 0.0),
        v=[1.0, 0.3, 0.2],
        ω=[2.5, -1.0, 2.0]) where T

    body = mechanism.bodies[1]

    halfside = body.shape.xyz[1] / 2.0

    if length(mechanism.contacts) > 0
        model = mechanism.contacts[1].model
        offset = model.collision.contact_radius
    else
        offset = 0.0
    end

    z = halfside + offset

    set_maximal_configurations!(body,
        x=x + [0.0, 0.0, z],
        q=q)
    set_maximal_velocities!(body,
        v=v,
        ω=ω)
end


text:
@testset "sigma clipping" begin
    x = [1, 2, 3]
    @test sigma_clip(x, 1, 1) ≈ [1.0, 2.0, 3.0] rtol = 1e-4
    @test sigma_clip(x, 1) == sigma_clip(x, 1, 1)

    y = [1, 2, 3, 4, 5, 6]
    @test sigma_clip(y, 1, 1) ≈ [1.62917, 2.0, 3.0, 4.0, 5.0, 5.37083] rtol = 1e-4

    # using different center
    @test sigma_clip(y, 1, center = 4) ≈ [2.1291713066130296, 2.1291713066130296, 3, 4, 5, 5.87082869338697]
    # using different std
    @test sigma_clip(y, 1, std = 1) ≈ [2.5, 2.5, 3, 4, 4.5, 4.5]

    @test sum(sigma_clip(y, 1, fill = NaN) .=== NaN) == 2
    @test sum(sigma_clip(y, 1, fill = NaN) .=== NaN) == 2
end

@testset "estimate_background interface" begin
    data = ones(100, 100)

    @test all(estimate_background(data, 20) .== estimate_background(data, (20, 20)))
    @test all(estimate_background(data, 20, filter_size = 3) .== estimate_background(data, (20, 20), filter_size = (3, 3)))
    @test size(estimate_background(data, 19, edge_method = :pad)[1]) == (114, 114)
    @test size(estimate_background(data, 19, edge_method = :crop)[1]) == (95, 95)
    X = rand(100, 100)
    @test estimate_background(data, location = median, rms = mean) == estimate_background(data, location = mean, rms = median)[[2, 1]]
    @test_throws MethodError estimate_background(data, (4, 4), edge_method = :yeet)
    @test_throws MethodError estimate_background(data, (4, 4, 4))
    @test_throws ErrorException estimate_background(data, 4, filter_size = 2)

    nan_bkg, nan_rms = estimate_background(fill(NaN, 100, 100), 20)
    @test nan_bkg ≈ zeros(100, 100)
    @test nan_rms ≈ zeros(100, 100)
end

@testset "flat background - $B, $S" for B in [median, mean, MMMBackground(), BiweightLocationBackground(), SourceExtractorBackground()], S in [StdRMS(), MADStdRMS(), BiweightScaleRMS()]
    data = ones(100, 100)

    bk, rms = estimate_background(data, location = B, rms = S)
    @test bk ≈ 1
    @test rms ≈ 0

    bk, rms = estimate_background(data, (25, 25), location = B, rms = S)
    @test median(bk) ≈ 1
    @test median(rms) ≈ 0
end

# regression test for failing to recognize NaN32 (or NaN16)
data = ones(Float32, 100, 100)
bk, rms = estimate_background(data, 3, edge_method = :pad)
@test !any(isnan.(bk))
@test !any(isnan.(rms))

@testset "interpolators" begin

    @testset "zoom interface" begin
        @test ZoomInterpolator(3) == ZoomInterpolator(3, 3) == ZoomInterpolator((3, 3))
    end

    @testset "trivial ones" begin
        z = ZoomInterpolator(4, 3)
        @test z(ones(3, 3)) ≈ ones(12, 9)
    end

end


text:
#Wigner 3j-symbol calculator
#17/12/2015

function Wigner3j(aa,bb,cc,alpha,beta,gamma)

# Find triangle coefficient using A. B. Deb (2007) script
triang = triangle_coeff(aa, bb, cc)

#Second term
sqrtterm = sqrt(factorial(aa+alpha)*factorial(aa-alpha)*
    factorial(bb+beta)*factorial(bb-beta)*factorial(cc+gamma)
    *factorial(cc-gamma))

#Finding the range of summation in Racah formula
tmin = zeros(3)
tmin[1] = bb - cc - alpha
tmin[2] = aa + beta - cc
tmin[3] = 0

tmax = zeros(3)
tmax[1] = aa + bb - cc
tmax[2] = aa - alpha
tmax[3] = bb + beta

rangei = maximum(tmin)
rangef = minimum(tmax)

#Sum over the function of t
sumscalar = 0
for tt = rangei:rangef

    sumscalar = sumscalar +
        ((-1)^tt)/tfunction3j(tt,aa,bb,cc,alpha,beta)

end

Wigner = ((-1)^float(aa-bb-gamma))*sqrt(triang)*sqrtterm*sumscalar

return Wigner
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
module DecidableDependentTypes

"""
TODO
----

Simply-typed λ calculus
-----------------------

τ ::= τ → τ
e ::= x | λx:τ.e | e e
Γ ::= x:τ, ... (types, terms)

         x:τ ϵ Γ
(T-VAR) ---------
        Γ ⊢ x : τ

          Γ, x:τ₁ ⊢ e : τ₂
(T-ABS) ---------------------
        Γ ⊢ λx:τ₁.e : τ₁ → τ₂

        Γ ⊢ e₁ : τ₁ → τ₂    Γ ⊢ e₂ : τ₁
(T-APP) -------------------------------
                 Γ ⊢ e₁e₂ : τ₂

(erase) er(x)=x, er(λx:τ.e)=λx.er(e), er(e e')=er(e)er(e')
"""

export stx_ob

#Types (just type synonyms for now)
const Ex = Expr
const Ty = Symbol

"""Ugly hack for struct comparison.
Tests syntactic (and therefore type) equality.
TODO: find better solution
Structs are ≂ if their fields are =="""
@generated function ≂(x, y)
    if !isempty(fieldnames(x)) && x == y
        mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), fieldnames(x))
    else
        :(x == y)
    end
end

# The syntax object
struct stx_ob
    e::Ex
    τ::Ty
end

#BOOL lang
# (bool->lit arg:boolean)
macro bool_to_lit(a)
    eval(a) ? (x -> x) : (x -> !x)
end
#DecidableDependentTypes.@bool_to_lit(:false)
#DecidableDependentTypes.@bool_to_lit(false)

#(λ (x y) (and #t ((bool->lit #f) x) ((bool->lit #f) y)))
macro make_clause(xs, bools, res)
    lits = map((x -> @eval @bool_to_lit $x), bools.args)
    clause_elements = map.(lits, xs.args)
    and_res = all(push!(clause_elements, res))
    return and_res
end
#DecidableDependentTypes.@make_clause((false, false), (false, false), true)
#DecidableDependentTypes.@make_clause((true, false), (true, false), false)
#DecidableDependentTypes.@make_clause((false, true), (false, true), true)
#DecidableDependentTypes.@make_clause((true, true), (true, true), true)

"""
xs : Tuple{Vararg{Bool}}
tbl : Tuple{Tuple{Vararg{Bool}}, Tuple{Vararg{Bool}}}
TODO: learn how to specify those types in the macro args.
"""
macro truth_table(xs, vars, r)
    clauses = []
    for i=1:length(vars.args)
        push!(clauses, @eval @make_clause($xs, $(vars.args[i]), $(r.args[i])))
    end
    return any(clauses)
end
#DecidableDependentTypes.@truth_table((false, true),((true, false), (false, true), (false, false)), (false, true, false))

function subset(x, y)
    @eval @truth_table(
        ($x, $y),
        ((false, false),
         (true, false),
         (false, true),
         (true, true)),
        (true, false, true, true))
end
#DecidableDependentTypes.subset(true, false)
#DecidableDependentTypes.subset(true, true)
#DecidableDependentTypes.subset(false, false)
#DecidableDependentTypes.subset(false, true)

end # module


text:
using Pkg

Pkg.add("PackageCompiler")
Pkg.add("Libdl")

using PackageCompiler, Libdl

f = open("/tmp/CompileHecke.jl", "w")
println(f, "using Hecke")
println(f, "using Pkg")
println(f, "Hecke.system(\"precompile.jl\")")
close(f)

PackageCompiler.create_sysimage([:Hecke], sysimage_path="/tmp/Hecke.$(Libdl.dlext)", precompile_execution_file="/tmp/CompileHecke.jl")

println("(re)start julia as")
println("\tjulia -J /tmp/Hecke.$(Libdl.dlext)")


text:
"""
Plan to cross-correlate every possible pair of columns.
"""
mutable struct P_xcorr{T<:Real}
	gsize::Vector{Int}
	cgsize::Vector{Vector{Int}}
	cg::Vector{Matrix{T}}
	g::Array{T,2}
	paconv::P_conv{T,1,1,1}
	cg_indices::Vector{Int64}
	norm_flag::Bool
end

function P_xcorr(nt::Int64, nr::Int64, T=Float64;
		 cg_indices::Vector{Int64}=collect(1:nr),
		 cglags=[nt-1, nt-1],
		 norm_flag=true)
	nts=sum(cglags)+1
	s=zeros(nts);
	g=zeros(nt);
	d=zeros(nt);

	# allocate g
	gsize=[nt, nr]
	gmat=zeros(gsize...)

	# allocate cg
	cgsize=[[nts, nr-cg_indices[i]+1] for i in eachindex(cg_indices)]
	cg=[zeros(cgsize[i]...) for i in eachindex(cg_indices)]

	# create a plan for conv
	paconv=P_conv(gsize=[nt], dsize=[nt], ssize=[nts], g=g, s=s, d=d, slags=cglags)

	# create a plan for xcorr
	pa=P_xcorr(gsize, cgsize, cg, gmat, paconv, cg_indices, norm_flag)

	# return plan
	return pa
end

function xcorr(g::AbstractArray, pa=P_xcorr(size(g,1), size(g,2), eltype(g)))

	mod!(pa, g=g)

	return pa.cg
end

"""
Use first colomn of g and cross-correlate with rest of columns of it.
And store results in cg.
After xcorr, the coeffcients are normalized with norm(g[:,1])
By default, cg has almost same positive and negative lags.
"""
function mod!(pa::P_xcorr; cg=pa.cg, g=pa.g)

	nr=size(g,2)
	cg_indices=pa.cg_indices
	lags=pa.paconv.slags
	norm_flag=pa.norm_flag

	any([(size(cg[ir]) ≠ (sum(lags)+1,nr-ir+1)) for ir in pa.cg_indices]) && error("size cg")

	irrr=0
	α=1.0
	nt=size(g,1)
	nts=sum(lags)+1
	ir1=0
	for ir in cg_indices
		ir1+=1
		cgx=cg[ir1]
		for i in 1:nt
			pa.paconv.d[i]=g[i,ir]
		end
		if(ir1==1 && norm_flag) # save scale reference for first column only
			for i in 1:nt
				α+=g[i,ir]*g[i,ir]
			end
			α = (iszero(α)) ? 1.0 : inv(α) # take inverse if not zero
		end
		for (iir2,ir2) in enumerate(ir:nr)
			for i in eachindex(pa.paconv.g)
				pa.paconv.g[i]=g[i,ir2]
			end
			mod!(pa.paconv, S())
		        for i in 1:nts
		        	cgx[i,iir2]=pa.paconv.s[i]
		        end
		        if(norm_flag)
		        	for i in 1:nts
		        		cgx[i,iir2]*=α
		        	end
			end
		end
	end
	return pa
end

"""
given dJdcg and g, computes dJdg
"""
function mod_grad!(dg::AbstractArray{T}, pa::P_xcorr{T}; dcg=pa.cg, g=pa.g) where {T}
	nr=size(g,2)
	nt=size(g,1)
	cg_indices=pa.cg_indices
	lags=pa.paconv.slags
	nts=sum(lags)+1
	for i in eachindex(dg)
		dg[i]=0.0
	end

	ir1=0
	for ir in cg_indices
		ir1+=1
		dcgx=dcg[ir1]

		for (iir2,ir2) in enumerate(ir:nr)

			for i in 1:nt
				pa.paconv.g[i]=g[i,ir2]
			end
			for i in 1:nts
				pa.paconv.s[i]=dcgx[i,iir2]
			end

			mod!(pa.paconv, D())

			for i in 1:nt
				dg[i,ir]+=pa.paconv.d[i]
			end
		end

		for i in 1:nt
			pa.paconv.d[i]=g[i,ir]
		end

		for (iir2,ir2) in enumerate(ir:nr)
			for i in 1:nts
				pa.paconv.s[i]=dcgx[i,iir2]
			end

			mod!(pa.paconv, G())
			for i in 1:nt
				dg[i,ir2] +=pa.paconv.g[i]
			end
		end
	end
	return dg

end

"""
Convert Array{Array{Float,2},1} to
Array{Float,2} and vice versa
"""
function cgmat!(cgmat, cg, flag; cg_indices::Vector{Int}=collect(1:length(cg)),
					cgmat_indices::Vector{Int}=collect(1:length(cg)))
	nr=length(cg);
	nt=size(cgmat,1)

	temp=zeros(nt)
	if(flag==1)
		for (ir,ircg) in enumerate(cg_indices)
			nr1=nr-ir+1
			cgx=cg[ir]
			for irr in 1:nr1
				irm1=findall(in(ircg), cgmat_indices)[1]
				irm2=findall(in(cg_indices[ir+irr-1]), cgmat_indices)[1]
				if(irm1<=irm2)
					igmat=sum([nr-i+1 for i in 1:irm1-1])+1+(irm2-irm1)
					for i in eachindex(temp)
						temp[i]=cgmat[i, igmat]
					end
				else
					igmat=sum([nr-i+1 for i in 1:irm2-1])+1+(irm1-irm2)
					for i in eachindex(temp)
						temp[i]=cgmat[nt-i+1, igmat]
					end
				end
				for i in eachindex(temp)
					cgx[i, irr]=temp[i]
				end
			end
		end
	end
	if(flag==-1)
		ntot=0
		for (ir,ircgmat) in enumerate(cgmat_indices)
			nr1=nr-ir+1
			for irr in 1:nr1
				irm1=findall(in(ircgmat), cg_indices)[1]
				irm2=findall(in(cgmat_indices[ir+irr-1]), cg_indices)[1]
				if(irm1<=irm2)
					cgx=cg[irm1]
					for i in eachindex(temp)
						temp[i]=cgx[i,irm2-irm1+1]
					end
				else
					cgx=cg[irm2]
					for i in eachindex(temp)
						temp[i]=cgx[nt-i+1,irm1-irm2+1]
					end
				end
				for i in eachindex(temp)
					cgmat[i, ntot+irr]=temp[i]
				end
			end
			ntot+=nr1
		end
	end
end

function cgmat(cgmat::AbstractMatrix, nr::Int; cg_indices=collect(1:nr))
	nts=size(cgmat,1)
	cg=[zeros(nts, nr-i+1) for i in 1:nr]
	cgmat!(cgmat, cg, 1, cg_indices=cg_indices)
	return cg
end

function cgmat(cg::Vector{Matrix}, nr::Int;cgmat_indices=collect(1:nr))
	nts=size(cg[1],1)
	cgmat=zeros(nts,binomial(nr, 2)+nr)
	cgmat!(cgmat, cg, -1, cgmat_indices=cgmat_indices)
	return cgmat
end


text:
module Linting

using FileIO

function lintfile(UserEnteredFile::String)
    f = open(UserEnteredFile)
    linesFromFIle = readlines(f)

    counter = 1

    for line in linesFromFIle
        if occursin("function", line)
            println("Found a function!")
        end
        println("$counter $line")
        counter += 1
    end
end
lintfile("/Users/logankilpatrick/.julia/dev/Linting.jl/src/Linting.jl")

end # module


text:
"""
    PlaneWaveBasis(b1::Vector{Float64}, b2::Vector{Float64},
            ps::Vector{Int}, qs::Vector{Int})

A 2D basis of plane waves, `ks = ps*b1 + qs*b2`.
"""
struct PlaneWaveBasis
    b1::Vector{Float64}
    b2::Vector{Float64}
    ps::Vector{Int}
    qs::Vector{Int}
    kxs::Vector{Float64}
    kys::Vector{Float64}
    function PlaneWaveBasis(b1::AbstractVector{Float64}, b2::AbstractVector{Float64},
                                        ps::AbstractVector{Int}, qs::AbstractVector{Int})
        kxs = [p*b1[1]+q*b2[1] for (p,q) in zip(ps,qs)]
        kys = [p*b1[2]+q*b2[2] for (p,q) in zip(ps,qs)]
        return new(b1, b2, ps, qs, kxs, kys)
    end
end

"""
    PlaneWaveBasis(geometry::Geometry, cutoff::Int)

Approximate a basis of plane waves truncated in a circle.

The circle has a diameter of `cutoff` Brillouin zones. Increasing the `cutoff`
will increase the number of plane waves leading to a more accurate solution.
It is assumed that `norm(b1) == norm(b2)`.
"""
function PlaneWaveBasis(geometry::Geometry, cutoff::Int)
    @assert isodd(cutoff)
    b1, b2 = as_to_bs(geometry.a1, geometry.a2)
    @assert norm(b1) ≈ norm(b2) # for now
    ps, qs = Int[], Int[]
    cutoff_radius = norm(b1) * cutoff / 2
    for p in -cutoff:cutoff, q in -cutoff:cutoff
        k = p*b1 + q*b2
        if norm(k) <= cutoff_radius
            push!(ps, p)
            push!(qs, q)
        end
    end
    return PlaneWaveBasis(b1, b2, ps, qs)
end

"""
    Solver(geometry::Geometry, cutoff_b1::Int, cutoff_b2::Int)

Approximate the geometry using a basis of plane waves truncated in a rhombus.

The rhombus has lengths `cutoff_b1` and `cutoff_b2` in the `b1` and `b2`
directions, respectively.
"""
function PlaneWaveBasis(geometry::Geometry, cutoff_b1::Int, cutoff_b2::Int)
    @assert isodd(cutoff_b1)
    @assert isodd(cutoff_b2)
    P = div(cutoff_b1, 2) # rounds down
    Q = div(cutoff_b2, 2) # rounds down
    ps = [p for p in -P:P for q in -Q:Q]
    qs = [q for p in -P:P for q in -Q:Q]
    b1, b2 = as_to_bs(geometry.a1, geometry.a2)
    return PlaneWaveBasis(b1, b2, ps, qs)
end

"""
    BrillouinZoneCoordinate(p::Float64, q::Float64, label::String="")

A labelled coordinate in the Brillouin zone.

The arguments `p` and `q` are the coefficients of reciprocal lattice vectors
`b1` and `b2`. The k-space coordinate, `k = p * b1 + q * b2`, is generated by
`get_k(coord::BrillouinZoneCoordinate, basis::PlaneWaveBasis)`. For example,
`BrillouinZoneCoordinate(0.5,0)` is on the edge of the first Brillouin zone.
"""
struct BrillouinZoneCoordinate
    p::Float64
    q::Float64
    label::String
    function BrillouinZoneCoordinate(p::Real, q::Real, label::String="")
        return new(p, q, label)
    end
end

"""
    get_k(coord::BrillouinZoneCoordinate, basis::PlaneWaveBasis)

Return the k-space coordinate of the `BrillouinZoneCoordinate` in a particular
`PlaneWaveBasis`, ie `k = p*b1 + q*b2`.
"""
function get_k(coord::BrillouinZoneCoordinate, basis::PlaneWaveBasis)
    return coord.p*basis.b1 + coord.q*basis.b2
end


text:
#unitary test GenericStatemod
x0 = ones(6)
state0 = GenericState(x0)

@test state0.current_time == nothing #Default value of start_time is void
@test state0.current_score == nothing
x1 = [1.0]
update!(state0, x = x1) #Check the update of state0
@test state0.x == x1
@test state0.current_time == nothing #start_time must be unchanged
@test state0.current_score == nothing

update!(state0, current_time = 1.0)
@test state0.x == x1 #must be unchanged
@test state0.current_time == 1.0
@test state0.current_score == nothing

reinit!(state0, x0)
@test state0.x == x0
@test state0.current_time == nothing
@test state0.current_score == nothing

update!(state0, x = x1)
reinit!(state0, current_time = 0.5)
@test state0.x == x1
@test state0.current_time == 0.5
@test state0.current_score == nothing


text:
################################################################################
#                               Linear                                         #
################################################################################

mutable struct Linear
    W::AbstractMatrix{T} where T
    b::AbstractVector{T} where T
end

function Linear(in_features::Int, out_features::Int)
    k_sqrt = sqrt(1 / in_features)
    d = Uniform(-k_sqrt, k_sqrt)
    return Linear(rand(d, out_features, in_features), rand(d, out_features))
end
Linear(in_out::Pair{Int, Int}) = Linear(in_out[1], in_out[2])

function Base.show(io::IO, l::Linear)
    o, i = size(l.W)
    print(io, "Linear($i=>$o)")
end

(l::Linear)(x::Union{AbstractVector{T}, AbstractMatrix{T}}) where T = l.W * x .+ l.b

################################################################################
#                               Sequential                                     #
################################################################################

mutable struct Sequential
    seq::Tuple
end

Sequential(args...) = Sequential(args)

function Base.show(io::IO, s::Sequential)
    println(io, "Sequential(")
    for m in s.seq
        println(io, "  $m,")
    end
    print(io, ")")
end

function (s::Sequential)(x)
    for m in s.seq
        x = m(x)
    end
    return x
end

################################################################################
#                           Convolutions                                       #
################################################################################

# Although internall 1D, 2D and 3D convolutions use the same convNd method,
# we keep different wrappers to make API easier to understand

################################# Conv1D #######################################

"""
1D convolutional layer.

Input and output dimensions: (sequence_length, channels, batch_size)
"""
mutable struct Conv1d
    W::AbstractArray{T,3} where T
    b::Union{AbstractVector{T} where T, Nothing}
    stride
    padding
    dilation
end

function Conv1d(in_channels::Int, out_channels::Int, kernel_size::Int;
                stride=1, padding=0, dilation=1, bias=true)
    kernel_tuple = (kernel_size,)
    # init weights same as in https://pytorch.org/docs/stable/nn.html#conv2d
    k_sqrt = sqrt(1 / (in_channels * prod(kernel_tuple)))
    d = Uniform(-k_sqrt, k_sqrt)
    W = rand(d, kernel_tuple..., in_channels, out_channels)
    b = bias ? rand(d, out_channels) : nothing
    return Conv1d(W, b, stride, padding, dilation)
end
Conv1d(in_out::Pair{Int, Int}, k) = Conv1d(in_out[1], in_out[2], k)

function Base.show(io::IO, c::Conv1d)
    k1, i, o = size(c.W)
    print(io, "Conv1d($i=>$o, $(k1))")
end

function (c::Conv1d)(x::AbstractArray{T,3}) where T
    y = convNd(x, c.W; stride=c.stride, padding=c.padding, dilation=c.dilation)
    if c.b != nothing
        y = y .+ reshape(c.b, (1, length(c.b), 1))  # TODO: check for 1D and 3D!
    end
    return y
end

################################# Conv2D #######################################

mutable struct Conv2d
    W::AbstractArray{T,4} where T
    b::Union{AbstractVector{T} where T, Nothing}
    stride
    padding
    dilation
end

function Conv2d(in_channels::Int, out_channels::Int, kernel_size::Union{Int, NTuple{2, Int}};
                stride=1, padding=0, dilation=1, bias=true)
    kernel_tuple = kernel_size isa Tuple ? kernel_size : (kernel_size, kernel_size)
    # init weights same as in https://pytorch.org/docs/stable/nn.html#conv2d
    k_sqrt = sqrt(1 / (in_channels * prod(kernel_tuple)))
    d = Uniform(-k_sqrt, k_sqrt)
    W = rand(d, kernel_tuple..., in_channels, out_channels)
    b = bias ? rand(d, out_channels) : nothing
    return Conv2d(W, b, stride, padding, dilation)
end
Conv2d(in_out::Pair{Int, Int}, k) = Conv2d(in_out[1], in_out[2], k)

function Base.show(io::IO, c::Conv2d)
    k1, k2, i, o = size(c.W)
    print(io, "Conv2d($i=>$o, $(k1)x$(k2))")
end

function (c::Conv2d)(x::AbstractArray{T,4}) where T
    y = convNd(x, c.W; stride=c.stride, padding=c.padding, dilation=c.dilation)
    if c.b != nothing
        y = y .+ reshape(c.b, (1, 1, length(c.b), 1))
    end
    return y
end

################################# Conv3D #######################################

mutable struct Conv3d
    W::AbstractArray{T,5} where T
    b::Union{AbstractVector{T} where T, Nothing}
    stride
    padding
    dilation
end

function Conv3d(in_channels::Int, out_channels::Int, kernel_size::Union{Int, NTuple{2, Int}};
                stride=1, padding=0, dilation=1, bias=true)
    kernel_tuple = kernel_size isa Tuple ? kernel_size : (kernel_size, kernel_size, kernel_size)
    # init weights same as in https://pytorch.org/docs/stable/nn.html#conv2d
    k_sqrt = sqrt(1 / (in_channels * prod(kernel_tuple)))
    d = Uniform(-k_sqrt, k_sqrt)
    W = rand(d, kernel_tuple..., in_channels, out_channels)
    b = bias ? rand(d, out_channels) : nothing
    return Conv3d(W, b, stride, padding, dilation)
end
Conv3d(in_out::Pair{Int, Int}, k) = Conv3d(in_out[1], in_out[2], k)

function Base.show(io::IO, c::Conv3d)
    k1, k2, k3, i, o = size(c.W)
    print(io, "Conv3d($i=>$o, $(k1)x$(k2)x$(k3))")
end

function (c::Conv3d)(x::AbstractArray{T,5}) where T
    y = convNd(x, c.W; stride=c.stride, padding=c.padding, dilation=c.dilation)
    if c.b != nothing
        y = y .+ reshape(c.b, (1, 1, length(c.b), 1))
    end
    return y
end

################################################################################
#                           Loss Functions                                     #
################################################################################

mutable struct NLLLoss
end

Base.show(io::IO, loss::NLLLoss) = print(io, "NLLLoss()")

(loss::NLLLoss)(logp::AbstractMatrix, c::Union{AbstractMatrix{<:Real}, AbstractVector{<:Real}}) =
    nllloss(logp, c)

mutable struct CrossEntropyLoss
end

Base.show(io::IO, loss::CrossEntropyLoss) = print(io, "CrossEntropyLoss()")

(loss::CrossEntropyLoss)(x::AbstractMatrix, c::Union{AbstractMatrix{<:Real}, AbstractVector{<:Real}}) =
    crossentropyloss(x, c)

mutable struct MSELoss
end

Base.show(io::IO, loss::MSELoss) = print(io, "MSELoss()")

(loss::MSELoss)(x::AbstractArray, x_target::AbstractArray) = mseloss(x, x_target)

################################################################################
#                                 Common                                       #
################################################################################

function trainmode!(m, train::Bool=true)
    # by default, recursively call trainmode!() on all fields
    T = typeof(m)
    for fld in fieldnames(T)
        trainmode!(getfield(m, fld), train)
    end
end
testmode!(m) = trainmode!(m, false)


text:
"""
    TraceRecord <: Record

Stores extra metrics (e.g., execution time, memory allocated) about a function
call for trace records.
"""
type TraceRecord <: Record
    date::Attribute
    level::Attribute
    levelnum::Attribute
    msg::Attribute
    name::Attribute
    pid::Attribute
    lookup::Attribute
    stacktrace::Attribute
    time::Attribute
    alloc::Attribute
    gctime::Attribute
    nalloc::Attribute

    function TraceRecord(
        name::AbstractString,
        level::AbstractString,
        levelnum::Int,
        msg::AbstractString,
        metrics::Dict=Dict{Symbol, Real}()
    )
        time = now()
        trace = Attribute(StackTrace, Memento.get_trace)

        new(
            Attribute(DateTime, () -> round(time, Base.Dates.Second)),
            Attribute(level),
            Attribute(levelnum),
            Attribute(AbstractString, Memento.get_msg(msg)),
            Attribute(name),
            Attribute(myid()),
            Attribute(StackFrame, Memento.get_lookup(trace)),
            trace,
            Attribute(get(metrics, :time, "NaN seconds")),
            Attribute(get(metrics, :alloc, "NaN bytes")),
            Attribute(get(metrics, :gctime, "NaN % gc time")),
            Attribute(get(metrics, :nalloc, "NaN")),
        )
    end
end

"""
    trace(f, logger, args...; kwargs...)

Wrapper function which measures the executions time and memory allocation
of function `f(args...; kwargs...)` and logs the resulting `TraceRecord`.
"""
function trace(f::Function, logger::Logger, args...; kwargs...)
    stats = Base.gc_num()
    elapsedtime = Base.time_ns()

    val = f(args...; kwargs...)

    elapsedtime = Base.time_ns() - elapsedtime
    diff = Base.GC_Diff(Base.gc_num(), stats)

    args = Dict(
        :time => "$(round(elapsedtime / 1e9, 6)) seconds",
        :alloc => datasize(diff.allocd, style=:bin, format="%.3f"),
        :gctime => "$(round(100 * diff.total_time / elapsedtime, 4)) % gc time",
        :nalloc => _humanize_alloc(Base.gc_alloc_count(diff))
    )

    msg = string(
        "`$(string(f))` <",
        "Time: $(args[:time]) ($(args[:gctime])), ",
        "Allocated: $(args[:alloc]) ($(args[:nalloc]))>"
    )

    rec = TraceRecord(logger.name, "trace", logger.levels["trace"], msg, args)
    @sync log(logger, rec)

    return val
end

trace(logger::Logger, f::Function, args...; kwargs...) = trace(f, logger, args...; kwargs...)

# Simple utility method for humanizing the number of allocations.
function _humanize_alloc(n)
    alloc, ma = Base.prettyprint_getunits(n, length(Base._cnt_units), Int64(1000))
    return "$(alloc)$(Base._cnt_units[ma])"
end


text:
# https://adventofcode.com/2021/day/15
using AdventOfCode
using Graphs
using SimpleWeightedGraphs

input = split.(readlines("2021/data/day_15.txt"),"") |> Base.Fix1(reduce,hcat) .|> parse
function makegraph(mat)
    g=SimpleWeightedGraph(Graphs.SimpleGraphs.grid(size(mat)), 1)
    for edge in edges(g)
        w=mat[src(edge)]+mat[dst(edge)]
        add_edge!(g,src(edge),dst(edge),w)
    end
    g
end

function part_1(input)
    g=makegraph(input)
    ds=dijkstra_shortest_paths(g,1)
    (ds.dists[end]-input[begin]+input[end])÷2
end

function part_2(input)
    fullmap=(repeat(input,5,5).+
        repeat((0:4).+(0:4)',inner=size(input)).-1) .%9 .+1
    g=makegraph(fullmap)
    ds=dijkstra_shortest_paths(g,1)
    (ds.dists[end]-input[begin]+input[end])÷2
end