text:
using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Atmos
using CLIMA.VariableTemplates
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray, )
else
  const ArrayTypes = (Array, )
end

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

include("mms_solution_generated.jl")

using CLIMA.Atmos
import CLIMA.Atmos: MoistureModel, temperature, pressure, soundspeed, total_specific_enthalpy

"""
    MMSDryModel

Assumes the moisture components is in the dry limit.
"""
struct MMSDryModel <: MoistureModel
end

function total_specific_enthalpy(moist::MoistureModel, orientation::Orientation, state::Vars, aux::Vars)
  zero(eltype(state))
end
function pressure(m::MMSDryModel, orientation::Orientation, state::Vars, aux::Vars)
  T = eltype(state)
  γ = T(7)/T(5)
  ρinv = 1 / state.ρ
  return (γ-1)*(state.ρe - ρinv/2 * sum(abs2, state.ρu))
end

function soundspeed(m::MMSDryModel, orientation::Orientation, state::Vars, aux::Vars)
  T = eltype(state)
  γ = T(7)/T(5)
  ρinv = 1 / state.ρ
  p = pressure(m, orientation, state, aux)
  sqrt(ρinv * γ * p)
end

function mms2_init_state!(state::Vars, aux::Vars, (x1,x2,x3), t)
  state.ρ = ρ_g(t, x1, x2, x3, Val(2))
  state.ρu = SVector(U_g(t, x1, x2, x3, Val(2)),
                     V_g(t, x1, x2, x3, Val(2)),
                     W_g(t, x1, x2, x3, Val(2)))
  state.ρe = E_g(t, x1, x2, x3, Val(2))
end

function mms2_source!(source::Vars, state::Vars, aux::Vars, t::Real)
  x1,x2,x3 = aux.coord
  source.ρ  = Sρ_g(t, x1, x2, x3, Val(2))
  source.ρu = SVector(SU_g(t, x1, x2, x3, Val(2)),
                      SV_g(t, x1, x2, x3, Val(2)),
                      SW_g(t, x1, x2, x3, Val(2)))
  source.ρe = SE_g(t, x1, x2, x3, Val(2))
end

function mms3_init_state!(state::Vars, aux::Vars, (x1,x2,x3), t)
  state.ρ = ρ_g(t, x1, x2, x3, Val(3))
  state.ρu = SVector(U_g(t, x1, x2, x3, Val(3)),
                     V_g(t, x1, x2, x3, Val(3)),
                     W_g(t, x1, x2, x3, Val(3)))
  state.ρe = E_g(t, x1, x2, x3, Val(3))
end

function mms3_source!(source::Vars, state::Vars, aux::Vars, t::Real)
  x1,x2,x3 = aux.coord
  source.ρ  = Sρ_g(t, x1, x2, x3, Val(3))
  source.ρu = SVector(SU_g(t, x1, x2, x3, Val(3)),
                      SV_g(t, x1, x2, x3, Val(3)),
                      SW_g(t, x1, x2, x3, Val(3)))
  source.ρe = SE_g(t, x1, x2, x3, Val(3))
end

# initial condition

function run(mpicomm, ArrayType, dim, topl, warpfun, N, timeend, FT, dt)

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = warpfun,
                                         )

  if dim == 2
    model = AtmosModel(NoOrientation(),
                       NoReferenceState(),
                       ConstantViscosityWithDivergence(FT(μ_exact)),
                       MMSDryModel(),
                       NoRadiation(),
                       mms2_source!,
                       InitStateBC(),
                       mms2_init_state!)
  else
    model = AtmosModel(NoOrientation(),
                       NoReferenceState(),
                       ConstantViscosityWithDivergence(FT(μ_exact)),
                       MMSDryModel(),
                       NoRadiation(),
                       mms3_source!,
                       InitStateBC(),
                       mms3_init_state!)
  end

  dg = DGModel(model,
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

  Q = init_ode_state(dg, FT(0))

  lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, ))
  # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  # Print some end of the simulation information
  engf = norm(Q)
  Qe = init_ode_state(dg, FT(timeend))

  engfe = norm(Qe)
  errf = euclidean_distance(Q, Qe)
  @info @sprintf """Finished
  norm(Q)                 = %.16e
  norm(Q) / norm(Q₀)      = %.16e
  norm(Q) - norm(Q₀)      = %.16e
  norm(Q - Qe)            = %.16e
  norm(Q - Qe) / norm(Qe) = %.16e
  """ engf engf/eng0 engf-eng0 errf errf / engfe
  errf
end

using Test
let
  MPI.Initialized() || MPI.Init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
      device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end

  polynomialorder = 4
  base_num_elem = 4

  expected_result = [1.5606126271800694e-01 5.3315235468477966e-03 2.2572701271274964e-04;
                     2.5754410198970158e-02 1.1781217145186288e-03 6.1752962472904071e-05]
lvls = integration_testing ? size(expected_result, 2) : 1

@testset "$(@__FILE__)" for ArrayType in ArrayTypes
for FT in (Float64,) #Float32)
  result = zeros(FT, lvls)
  for dim = 2:3
    for l = 1:lvls
      if dim == 2
        Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
        brickrange = (range(FT(0); length=Ne[1]+1, stop=1),
                      range(FT(0); length=Ne[2]+1, stop=1))
        topl = BrickTopology(mpicomm, brickrange,
                             periodicity = (false, false))
        dt = 1e-2 / Ne[1]
        warpfun = (x1, x2, _) -> begin
          (x1 + sin(x1*x2), x2 + sin(2*x1*x2), 0)
        end

      elseif dim == 3
        Ne = (2^(l-1) * base_num_elem, 2^(l-1) * base_num_elem)
        brickrange = (range(FT(0); length=Ne[1]+1, stop=1),
                      range(FT(0); length=Ne[2]+1, stop=1),
        range(FT(0); length=Ne[2]+1, stop=1))
        topl = BrickTopology(mpicomm, brickrange,
                             periodicity = (false, false, false))
        dt = 5e-3 / Ne[1]
        warpfun = (x1, x2, x3) -> begin
          (x1 + (x1-1/2)*cos(2*π*x2*x3)/4,
           x2 + exp(sin(2π*(x1*x2+x3)))/20,
          x3 + x1/4 + x2^2/2 + sin(x1*x2*x3))
        end
      end
      timeend = 1
      nsteps = ceil(Int64, timeend / dt)
      dt = timeend / nsteps

      @info (ArrayType, FT, dim)
      result[l] = run(mpicomm, ArrayType, dim, topl, warpfun,
                      polynomialorder, timeend, FT, dt)
      @test result[l] ≈ FT(expected_result[dim-1, l])
    end
    if integration_testing
      @info begin
        msg = ""
        for l = 1:lvls-1
          rate = log2(result[l]) - log2(result[l+1])
          msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
        end
        msg
      end
    end
  end
end
end
end

#nothing

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
####
# Matrix memory layout traits
#
# if MemoryLayout(A) returns BandedColumnMajor, you must override
# pointer and leadingdimension
# in addition to the banded matrix interface
####

struct BandedLayout <: AbstractBandedLayout end
struct BandedColumns{ML} <: AbstractBandedLayout end
struct BandedRows{ML} <: AbstractBandedLayout end

const BandedColumnMajor = BandedColumns{<:AbstractColumnMajor}
const BandedRowMajor = BandedRows{<:AbstractColumnMajor}
BandedColumnMajor() = BandedColumns{DenseColumnMajor}()
BandedRowMajor() = BandedRows{DenseRowMajor}()

transposelayout(M::BandedColumns{ML}) where ML = BandedRows{ML}()
transposelayout(M::BandedRows{ML}) where ML = BandedColumns{ML}()
conjlayout(::Type{<:Complex}, ::M) where M<:AbstractBandedLayout = ConjLayout{M}()

# bandedrowsdata(Ac::Adjoint) = bandeddata(parent(Ac))'
# bandedrowsdata(Ac::Transpose) = transpose(bandeddata(parent(Ac)))
# bandedrowsdata(A) = bandedrowsdata(MemoryLayout(A), A)

# Here we override broadcasting for banded matrices.
# The design is to to exploit the broadcast machinery so that
# banded matrices that conform to the banded matrix interface but are not
# <: AbstractBandedMatrix can get access to fast copyto!, lmul!, rmul!, axpy!, etc.
# using broadcast variants (B .= A, B .= 2.0 .* A, etc.)

struct BandedStyle <: AbstractArrayStyle{2} end
BandedStyle(::Val{2}) = BandedStyle()
BroadcastStyle(::Type{<:AbstractBandedMatrix}) = BandedStyle()
BroadcastStyle(::Type{<:Adjoint{<:Any,<:AbstractBandedMatrix}}) = BandedStyle()
BroadcastStyle(::Type{<:Transpose{<:Any,<:AbstractBandedMatrix}}) = BandedStyle()
BroadcastStyle(::Type{<:SubArray{<:Any,2,<:AbstractBandedMatrix,<:NTuple{2,AbstractUnitRange{Int}}}}) = BandedStyle()
BroadcastStyle(::DefaultArrayStyle{2}, ::BandedStyle) = BandedStyle()
BroadcastStyle(::BandedStyle, ::DefaultArrayStyle{2}) = BandedStyle()

size(bc::Broadcasted{BandedStyle}) = length.(axes(bc))
isbanded(bc::Broadcasted{BandedStyle}) = true

####
# Default to standard Array broadcast
#
# This is because, for example, exp.(B) is not banded
####

copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle}) =
   copyto!(dest, Broadcasted{DefaultArrayStyle{2}}(bc.f, bc.args, bc.axes))
##
# copyto!
##

_copyto!(dest_L::AbstractBandedLayout, src_L::AbstractBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) =
    _banded_broadcast!(dest, identity, src, dest_L, src_L)

function checkbroadcastband(dest, sizesrc, bndssrc)
    size(dest) == sizesrc || throw(DimensionMismatch())
    min(sizesrc[1],bndssrc[1]+1) ≤ min(bandwidth(dest,1)+1,size(dest,1))  &&
        min(sizesrc[2],bndssrc[2]+1) ≤ min(bandwidth(dest,2)+2,size(dest,2)) || throw(BandError(dest,size(dest,2)-1))
end

###########
# matrix broadcast
###########

function checkzerobands(dest, f, A::AbstractMatrix)
    m,n = size(A)
    d_l, d_u = bandwidths(dest)
    l, u = bandwidths(A)

    if (l,u) ≠ (d_l,d_u)
        for j = 1:n
            for k = max(1,j-u) : min(j-d_u-1,m)
                iszero(f(A[k,j])) || throw(BandError(dest,j-k))
            end
            for k = max(1,j+d_l+1) : min(j+l,m)
                iszero(f(A[k,j])) || throw(BandError(dest,j-k))
            end
        end
    end
end

function _banded_broadcast!(dest::AbstractMatrix, f, src::AbstractMatrix{T}, _1, _2) where T
    z = f(zero(T))
    iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, src)))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    s_l, s_u = bandwidths(src)
    (d_l ≥ min(s_l,m-1) && d_u ≥ min(s_u,n-1)) || throw(BandError(dest))

    for j=1:n
        for k = max(1,j-d_u):min(j-s_u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-s_u):min(j+s_l,m)
            inbands_setindex!(dest, f(inbands_getindex(src, k, j)), k, j)
        end
        for k = max(1,j+s_l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

# columns to the left (start) and right (stop) of the non-zero columns _nzcols
@inline _startzcols((m,n), (l,u)) = 1:-l
@inline _nzcols((m,n), (l,u)) = max(1-l,1):min(m+u,n)
@inline _stopzcols((m,n), (l,u)) = m+u+1:n

# columns in the "start up", bulk and "stop
@inline _startcols((m,n), (l,u)) = max(1-l,1):min(u,n)
@inline _bulkcols((m,n), (l,u)) = max(u+1,1):min(m-l,n)
@inline _stopcols((m,n), (l,u)) = max(m-l+1,u+1,1):min(m+u,n)

# the non-zero rows in the data in column j
@inline _startcolrange((m,n), (l,u), j) = u-j+2:min(u+m+1-j,l+u+1)
@inline _bulkcolrange((m,n), (l,u), j) = 1:l+u+1
@inline _stopcolrange((m,n), (l,u), j) = 1:u+m+1-j

# how many rows we shift down
@inline _bulkshift((l,u), j) = j-u-1

@inline _startzcols(A) = _startzcols(size(A), bandwidths(A))
@inline _nzcols(A) = _nzcols(size(A), bandwidths(A))
@inline _stopzcols(A) = _stopzcols(size(A), bandwidths(A))
@inline _startcols(A) = _startcols(size(A), bandwidths(A))
@inline _bulkcols(A) = _bulkcols(size(A), bandwidths(A))
@inline _stopcols(A) = _stopcols(size(A), bandwidths(A))
@inline _startcolrange(A, j) = _startcolrange(size(A), bandwidths(A), j)
@inline _bulkcolrange(A, j) = _bulkcolrange(size(A), bandwidths(A), j)
@inline _stopcolrange(A, j) = _stopcolrange(size(A), bandwidths(A), j)
@inline _colrange(A, j) = _colrange(size(A), bandwidths(A), j)
@inline _colshift(A::AbstractMatrix, j) = _colshift(bandwidths(A), j)
@inline _bulkshift(A::AbstractMatrix, j) = _bulkshift(bandwidths(A), j)

@inline function _colrange((m,n), (l,u), j)
    j ≤ u && return _startcolrange((m,n), (l,u), j)
    j ≤ m-l && return _bulkcolrange((m,n), (l,u), j)
    return _stopcolrange((m,n), (l,u), j)
end
@inline function _colshift((l,u), j)
    j ≤ u && return 0
    _bulkshift((l,u),j)
end

function _intersectcolrange(kr_A,kr_d,sh)
    if sh ≤ 0
        kr_dA = first(kr_d):min(last(kr_d),first(kr_d)+length(kr_A)+sh-1)
        kr_Ad = range(first(kr_A)-sh; length=length(kr_dA))
    else
        kr_dA = first(kr_d)+sh:min(last(kr_d),first(kr_d)+length(kr_A)+sh-1)
        kr_Ad = range(first(kr_A); length=length(kr_dA))
    end
    kr_dA, kr_Ad
end

function _banded_broadcast!(dest::AbstractMatrix, f, A::AbstractMatrix{T}, ::BandedColumns, ::BandedColumns) where T
    z = f(zero(T))
    iszero(z) || checkbroadcastband(dest, size(A), bandwidths(broadcasted(f, A)))

    m,n = size(dest)
    (m == 0 || n == 0) && return dest
    data_d,data_A = bandeddata(dest), bandeddata(A)

    checkzerobands(dest, f, A)
    # copy top left entries
    jr = _bulkcols(A) ∩ _bulkcols(dest)
    for j = last(_startzcols(dest))+1:min(first(jr)-1,n)
        # don't bother optimising
        kr_A = _colrange(A,j)
        kr_d = _colrange(dest,j)
        sh = _colshift(A,j)-_colshift(dest,j)
        kr_dA,kr_Ad = _intersectcolrange(kr_A,kr_d,sh)
        if isempty(kr_Ad)
            data_d[kr_d,j] .= z
        else
            data_d[first(kr_d):first(kr_dA)-1,j] .= z
            data_d[kr_dA,j] .= f.(view(data_A,kr_Ad,j))
            data_d[last(kr_dA)+1:last(kr_d),j] .= z
        end
    end

    # the bulk
    kr_A = _bulkcolrange(A,first(jr))
    kr_d = _bulkcolrange(dest,first(jr))
    sh = _bulkshift(A,first(jr))-_bulkshift(dest,first(jr))
    kr_dA,kr_Ad = _intersectcolrange(kr_A,kr_d,sh)
    if isempty(kr_Ad)
        data_d[kr_d,jr] .= z
    else
        data_d[first(kr_d):first(kr_dA)-1,jr] .= z
        data_d[kr_dA,jr] .= f.(view(data_A,kr_Ad,jr))
        data_d[last(kr_dA)+1:last(kr_d),jr] .= z
    end

    # bottom right
    for j = last(jr)+1:min(last(_stopcols(dest)),n)
        kr_A = _colrange(A,j)
        kr_d = _colrange(dest,j)
        sh = _colshift(A,j)-_colshift(dest,j)
        kr_dA,kr_Ad = _intersectcolrange(kr_A,kr_d,sh)
        if isempty(kr_Ad)
            data_d[kr_d,j] .= z
        else
            data_d[first(kr_d):first(kr_dA)-1,j] .= z
            data_d[kr_dA,j] .= f.(view(data_A,kr_Ad,j))
            data_d[last(kr_dA)+1:last(kr_d),j] .= z
        end
    end

    dest
end

function copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix}})
    (A,) = bc.args
    _banded_broadcast!(dest, bc.f, A, MemoryLayout(typeof(dest)), MemoryLayout(typeof(A)))
end

###########
# matrix-number broadcast
###########

function _banded_broadcast!(dest::AbstractMatrix, f, (src,x)::Tuple{AbstractMatrix{T},Number}, _1, _2) where T
    z = f(zero(T), x)
    iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, src,x)))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    s_l, s_u = bandwidths(src)
    (d_l ≥ min(s_l,m-1) && d_u ≥ min(s_u,n-1)) || throw(BandError(dest))

    for j = rowsupport(dest)
        for k = max(1,j-d_u):min(j-s_u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-s_u):min(j+s_l,m)
            inbands_setindex!(dest, f(inbands_getindex(src, k, j), x), k, j)
        end
        for k = max(1,j+s_l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

function _banded_broadcast!(dest::AbstractMatrix, f, (x,src)::Tuple{Number,AbstractMatrix{T}}, _1, _2) where T
    z = f(x, zero(T))
    iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, x,src)))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    s_l, s_u = bandwidths(src)
    (d_l ≥ min(s_l,m-1) && d_u ≥ min(s_u,n-1)) || throw(BandError(dest))

    for j = rowsupport(dest)
        for k = max(1,j-d_u):min(j-s_u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-s_u):min(j+s_l,m)
            inbands_setindex!(dest, f(x, inbands_getindex(src, k, j)), k, j)
        end
        for k = max(1,j+s_l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

function _banded_broadcast!(dest::AbstractMatrix, f, (src,x)::Tuple{AbstractMatrix{T},Number}, ::BandedColumns, ::BandedColumns) where T
    z = f(zero(T),x)
    iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, src,x)))

    l,u = bandwidths(src)
    λ,μ = bandwidths(dest)
    m,n = size(src)
    data_d,data_s = bandeddata(dest), bandeddata(src)
    d_m, d_n = size(data_d)
    if (l,u) == (λ,μ)
        data_d .= f.(data_s, x)
    elseif μ > u && λ > l
        fill!(view(data_d,1:min(μ-u,d_m),:), z)
        fill!(view(data_d,μ+l+2:μ+λ+1,:), z)
        view(data_d,μ-u+1:μ+l+1,:) .= f.(data_s,x)
    elseif μ > u
        fill!(view(data_d,1:min(μ-u,d_m),:), z)
        for b = λ+1:l
            any(!iszero, view(data_s,u+b+1,1:min(m-b,n))) && throw(BandError(dest,b))
        end
        view(data_d,μ-u+1:μ+λ+1,:) .= f.(view(data_s,1:u+λ+1,:),x)
    elseif λ > l
        for b = μ+1:u
            any(!iszero, view(data_s,u-b+1,b+1:n)) && throw(BandError(dest,b))
        end
        fill!(view(data_d,μ+l+2:μ+λ+1,:), z)
        view(data_d,1:μ+l+1,:) .= f.(view(data_s,u-μ+1:u+l+1,:),x)
    else # μ < u && λ < l
        for b = μ+1:u
            any(!iszero, view(data_s,u-b+1,b+1:n)) && throw(BandError(dest,b))
        end
        for b = λ+1:l
            any(!iszero, view(data_s,u+b+1,1:min(m-b,n))) && throw(BandError(dest,b))
        end
        data_d .= f.(view(data_s,u-μ+1:u+λ+1,:),x)
    end

    dest
end

function _banded_broadcast!(dest::AbstractMatrix, f, (x, src)::Tuple{Number,AbstractMatrix{T}}, ::BandedColumns, ::BandedColumns) where T
    z = f(x, zero(T))
    iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, x,src)))

    l,u = bandwidths(src)
    λ,μ = bandwidths(dest)
    m,n = size(src)
    data_d,data_s = bandeddata(dest), bandeddata(src)
    d_m, d_n = size(data_d)
    if (l,u) == (λ,μ)
        data_d .= f.(x, data_s)
    elseif μ > u && λ > l
        fill!(view(data_d,1:min(μ-u,d_m),:), z)
        fill!(view(data_d,μ+l+2:μ+λ+1,:), z)
        view(data_d,μ-u+1:μ+l+1,:) .= f.(x, data_s)
    elseif μ > u
        fill!(view(data_d,1:min(μ-u,d_m),:), z)
        for b = λ+1:l
            any(!iszero, view(data_s,u+b+1,1:min(m-b,n))) && throw(BandError(dest,b))
        end
        view(data_d,μ-u+1:μ+λ+1,:) .= f.(x, view(data_s,1:u+λ+1,:))
    elseif λ > l
        for b = μ+1:u
            any(!iszero, view(data_s,u-b+1,b+1:n)) && throw(BandError(dest,b))
        end
        fill!(view(data_d,μ+l+2:μ+λ+1,:), z)
        view(data_d,1:μ+l+1,:) .= f.(x, view(data_s,u-μ+1:u+l+1,:))
    else # μ < u && λ < l
        for b = μ+1:u
            any(!iszero, view(data_s,u-b+1,b+1:n)) && throw(BandError(dest,b))
        end
        for b = λ+1:l
            any(!iszero, view(data_s,u+b+1,1:min(m-b,n))) && throw(BandError(dest,b))
        end
        data_d .= f.(x, view(data_s,u-μ+1:u+λ+1,:))
    end

    dest
end

function copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:Number}})
    (A,x) = bc.args
    _banded_broadcast!(dest, bc.f, (A, x), MemoryLayout(typeof(dest)), MemoryLayout(typeof(A)))
end

function copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:Number,<:AbstractMatrix}})
    (x,A) = bc.args
    _banded_broadcast!(dest, bc.f, (x,A), MemoryLayout(typeof(dest)), MemoryLayout(typeof(A)))
end

###############
# matrix-vector broadcast
###############

_banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractVector{T},AbstractMatrix{V}}, _1, _2) where {T,V} =
    _left_colvec_banded_broadcast!(dest, f, (A,B), _1, _2)

function _left_colvec_banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractVecOrMat{T},AbstractMatrix{V}}, _1, _2) where {T,V}
    z = f(zero(T), zero(V))
    bc = broadcasted(f, A, B)
    l, u = bandwidths(bc)
    iszero(z) || checkbroadcastband(dest, size(bc), (l,u))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    A_l, A_u = _broadcast_bandwidths((m-1,n-1),A)
    B_l, B_u = bandwidths(B)
    (d_l ≥ min(l,m-1) && d_u ≥ min(u,n-1)) || throw(BandError(dest))

    for j=1:n
        for k = max(1,j-d_u):min(j-u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-d_u,j-A_u):min(j-B_u-1,j+d_l,m)
            inbands_setindex!(dest, f(A[k], zero(V)), k, j)
        end
        for k = max(1,j-d_u,j-B_u):min(j-A_u-1,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j-min(A_u,B_u)):min(j+min(A_l,B_l),m)
            inbands_setindex!(dest, f(A[k], inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j-d_u,j+B_l+1):min(j+A_l,j+d_l,m)
            inbands_setindex!(dest, f(A[k], zero(V)), k, j)
        end
        for k = max(1,j-d_u,j+A_l+1):min(j+B_l,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j+l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

_banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix{T},AbstractVector{V}}, _1, _2) where {T,V} =
    _right_colvec_banded_broadcast!(dest, f, (A,B), _1, _2)

function _right_colvec_banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix{T},AbstractVecOrMat{V}}, _1, _2) where {T,V}
    z = f(zero(T), zero(V))
    bc = broadcasted(f, A, B)
    l, u = bandwidths(bc)
    iszero(z) || checkbroadcastband(dest, size(bc), (l,u))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    A_l, A_u = bandwidths(A)
    B_l, B_u = _broadcast_bandwidths((m-1,n-1),B)
    (d_l ≥ min(l,m-1) && d_u ≥ min(u,n-1)) || throw(BandError(dest))

    for j=1:n
        for k = max(1,j-d_u):min(j-u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-d_u,j-A_u):min(j-B_u-1,j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), zero(V)), k, j)
        end
        for k = max(1,j-d_u,j-B_u):min(j-A_u-1,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), B[k]), k, j)
        end
        for k = max(1,j-min(A_u,B_u)):min(j+min(A_l,B_l),m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), B[k]), k, j)
        end
        for k = max(1,j-d_u,j+B_l+1):min(j+A_l,j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), zero(V)), k, j)
        end
        for k = max(1,j-d_u,j+A_l+1):min(j+B_l,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), B[k]), k, j)
        end
        for k = max(1,j+l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

function _left_rowvec_banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix{T},AbstractMatrix{V}}, _1, _2) where {T,V}
    @assert size(A,1) == 1
    z = f(zero(T), zero(V))
    bc = broadcasted(f, A, B)
    l, u = bandwidths(bc)
    iszero(z) || checkbroadcastband(dest, size(bc), (l,u))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    A_l, A_u = _broadcast_bandwidths((m-1,n-1),A)
    B_l, B_u = bandwidths(B)
    (d_l ≥ min(l,m-1) && d_u ≥ min(u,n-1)) || throw(BandError(dest))

    for j=1:n
        for k = max(1,j-d_u):min(j-u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-d_u,j-A_u):min(j-B_u-1,j+d_l,m)
            inbands_setindex!(dest, f(A[j], zero(V)), k, j)
        end
        for k = max(1,j-d_u,j-B_u):min(j-A_u-1,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j-min(A_u,B_u)):min(j+min(A_l,B_l),m)
            inbands_setindex!(dest, f(A[j], inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j-d_u,j+B_l+1):min(j+A_l,j+d_l,m)
            inbands_setindex!(dest, f(A[j], zero(V)), k, j)
        end
        for k = max(1,j-d_u,j+A_l+1):min(j+B_l,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j+l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

function _right_rowvec_banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix{T},AbstractMatrix{V}}, _1, _2) where {T,V}
    @assert size(B,1) == 1
    z = f(zero(T), zero(V))
    bc = broadcasted(f, A, B)
    l, u = bandwidths(bc)
    iszero(z) || checkbroadcastband(dest, size(bc), (l,u))
    m,n = size(dest)

    d_l, d_u = bandwidths(dest)
    A_l, A_u = bandwidths(A)
    B_l, B_u = _broadcast_bandwidths((m-1,n-1),B)
    (d_l ≥ min(l,m-1) && d_u ≥ min(u,n-1)) || throw(BandError(dest))

    for j=1:n
        for k = max(1,j-d_u):min(j-u-1,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-d_u,j-A_u):min(j-B_u-1,j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), zero(V)), k, j)
        end
        for k = max(1,j-d_u,j-B_u):min(j-A_u-1,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), B[j]), k, j)
        end
        for k = max(1,j-min(A_u,B_u)):min(j+min(A_l,B_l),m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), B[j]), k, j)
        end
        for k = max(1,j-d_u,j+B_l+1):min(j+A_l,j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), zero(V)), k, j)
        end
        for k = max(1,j-d_u,j+A_l+1):min(j+B_l,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), B[j]), k, j)
        end
        for k = max(1,j+l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractVector,<:AbstractMatrix}}) =
    _banded_broadcast!(dest, bc.f, bc.args, MemoryLayout(dest), MemoryLayout.(bc.args))

copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:AbstractVector}}) =
    _banded_broadcast!(dest, bc.f, bc.args, MemoryLayout(dest), MemoryLayout.(bc.args))

################
# matrix-matrix broadcast
################

function _banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix{T},AbstractMatrix{V}}, _1, _2) where {T,V}
    z = f(zero(T), zero(V))
    bc = broadcasted(f, A, B)
    l, u = bandwidths(bc)
    iszero(z) || checkbroadcastband(dest, size(bc), (l,u))
    m,n = size(dest)
    if size(A) ≠ (m,n)
        size(A,2) == 1 && return _left_colvec_banded_broadcast!(dest, f, (A,B), _1, _2)
        return _left_rowvec_banded_broadcast!(dest, f, (A,B), _1, _2)
    end
    if size(B) ≠ (m,n)
        size(B,2) == 1 && return _right_colvec_banded_broadcast!(dest, f, (A,B), _1, _2)
        return _right_rowvec_banded_broadcast!(dest, f, (A,B), _1, _2)
    end

    d_l, d_u = bandwidths(dest)
    A_l, A_u = bandwidths(A)
    B_l, B_u = bandwidths(B)
    (d_l ≥ min(l,m-1) && d_u ≥ min(u,n-1)) || throw(BandError(dest))

    for j=1:n
        for k = max(1,j-d_u):min(j-u-1,j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
        for k = max(1,j-A_u,j-d_u):min(j-B_u-1,j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), zero(V)), k, j)
        end
        for k = max(1,j-B_u,j-d_u):min(j-A_u-1,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j-min(A_u,B_u),j-d_u):min(j+min(A_l,B_l),j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j+B_l+1,j-d_u):min(j+A_l,j+d_l,m)
            inbands_setindex!(dest, f(inbands_getindex(A, k, j), zero(V)), k, j)
        end
        for k = max(1,j+A_l+1,j-d_u):min(j+B_l,j+d_l,m)
            inbands_setindex!(dest, f(zero(T), inbands_getindex(B, k, j)), k, j)
        end
        for k = max(1,j-d_u,j+l+1):min(j+d_l,m)
            inbands_setindex!(dest, z, k, j)
        end
    end
    dest
end

function checkzerobands(dest, f, (A,B)::Tuple{AbstractMatrix,AbstractMatrix})
    m,n = size(A)
    d_l, d_u = bandwidths(dest)
    A_l, A_u = bandwidths(A)
    B_l, B_u = bandwidths(B)
    l, u = max(A_l,B_l), max(A_u,B_u)

    for j = 1:n
        for k = max(1,j-u) : min(j-d_u-1,m)
            iszero(f(A[k,j], B[k,j])) || throw(BandError(dest,b))
        end
        for k = max(1,j+d_l+1) : min(j+l,m)
            iszero(f(A[k,j], B[k,j])) || throw(BandError(dest,b))
        end
    end
end

function _banded_broadcast!(dest::AbstractMatrix, f, (A,B)::Tuple{AbstractMatrix,AbstractMatrix}, ::BandedColumns, ::Tuple{<:BandedColumns,<:BandedColumns})
    z = f(zero(eltype(A)), zero(eltype(B)))
    bc = broadcasted(f, A, B)
    l, u = bandwidths(bc)
    (-u ≥ size(dest,1) || -l ≥ size(dest,2)) && return fill!(dest, z)
    iszero(z) || checkbroadcastband(dest, size(bc), (l,u))

    if size(A) ≠ size(dest) || size(B) ≠ size(dest)
        # special broadcast
        return _banded_broadcast!(dest, f, (A,B), UnknownLayout(), UnknownLayout())
    end

    A_l,A_u = bandwidths(A)
    B_l,B_u = bandwidths(B)
    -A_l > A_u && return broadcast!(f, dest, zero(eltype(A)), B)
    -B_l > B_u && return broadcast!(f, dest, A, zero(eltype(B)))
    d_l,d_u = bandwidths(dest)
    m,n = size(A)
    data_d,data_A, data_B = bandeddata(dest), bandeddata(A), bandeddata(B)

    if (d_l,d_u) == (A_l,A_u) == (B_l,B_u)
        data_d .= f.(data_A,data_B)
    else
        max_l,max_u = max(A_l,B_l,d_l),max(A_u,B_u,d_u)
        min_l,min_u = min(A_l,B_l,d_l),min(A_u,B_u,d_u)
        checkzerobands(dest, f, (A,B))

        # fill extra bands in dest
        fill!(view(data_d,1:d_u-max(A_u,B_u),:), z)
        fill!(view(data_d,d_u+max(A_l,B_l)+2:size(data_d,1),:), z)
        # between A and B if non-overlapping
        bs = max(-d_u,B_l+1):min(-A_u-1,d_l)
        fill!(view(data_d, bs .+ d_u .+ 1, :), z)
        bs = max(-d_u,A_l+1):min(-B_u-1,d_l)
        fill!(view(data_d, bs .+ d_u .+ 1, :), z)

        # construct where B upper is zero
        # this is from band A_u:B_u+1
        bs = max(-d_u,-A_u):min(-1-B_u,A_l,d_l)
        data_d_u_A = view(data_d, bs .+ d_u .+ 1, :)
        data_A_u_A = view(data_A, bs .+ A_u .+ 1, :)
        data_d_u_A .= f.(data_A_u_A, zero(eltype(B)))

        # construct where A upper is zero
        # this is from band B_u:min(A_u+1,-B_l)
        bs = max(-d_u,-B_u):min(-1-A_u,B_l,d_l)
        data_d_u_B = view(data_d, bs .+ d_u .+ 1, :)
        data_B_u_B = view(data_B, bs .+ B_u .+ 1, :)
        data_d_u_B .= f.(zero(eltype(A)), data_B_u_B)

        # construct where A and B are non-zero
        data_d̃ = view(data_d, d_u-min_u+1 : d_u+min_l+1, :)
        data_Ã = view(data_A, A_u-min_u+1 : A_u+min_l+1, :)
        data_B̃ = view(data_B, B_u-min_u+1 : B_u+min_l+1, :)
        data_d̃ .= f.(data_Ã,data_B̃)

        # construct where A lower is zero
        bs = max(A_l+1,-B_u,-d_u):min(B_l,d_l)
        data_d_l_B = view(data_d, bs .+ d_u .+ 1, :)
        data_B_l_B = view(data_B, bs .+ B_u .+ 1, :)
        data_d_l_B .= f.(zero(eltype(A)), data_B_l_B)

        # construct where B lower is zero
        bs = max(B_l+1,-A_u,-d_u):min(A_l,d_l)
        data_d_l_A = view(data_d, bs .+ d_u .+ 1, :)
        data_A_l_A = view(data_A, bs .+ A_u .+ 1, :)
        data_d_l_A .= f.(data_A_l_A, zero(eltype(B)))
    end

    dest
end

copyto!(dest::AbstractArray, bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:AbstractMatrix}}) =
    _banded_broadcast!(dest, bc.f, bc.args, MemoryLayout(typeof(dest)), map(MemoryLayout,bc.args))

# override copy in case data has special broadcast
_default_banded_broadcast(bc::Broadcasted{Style}, _) where Style = Base.invoke(copy, Tuple{Broadcasted{Style}}, bc)
_default_banded_broadcast(bc::Broadcasted) = _default_banded_broadcast(bc, axes(bc))

_banded_broadcast(f, args::Tuple, _...) = _default_banded_broadcast(broadcasted(f, args...))
_banded_broadcast(f, arg, _...) = _default_banded_broadcast(broadcasted(f, arg))

function _banded_broadcast(f, A::AbstractMatrix{T}, ::BandedColumns) where T
    iszero(f(zero(T))) || return _default_banded_broadcast(broadcasted(f, A))
    _BandedMatrix(f.(bandeddata(A)), axes(A,1), bandwidths(A)...)
end
function _banded_broadcast(f, (src,x)::Tuple{AbstractMatrix{T},Number}, ::BandedColumns) where T
    iszero(f(zero(T),x)) || return _default_banded_broadcast(broadcasted(f, src,x))
    _BandedMatrix(f.(bandeddata(src),x), axes(src,1), bandwidths(src)...)
end
function _banded_broadcast(f, (x, src)::Tuple{Number,AbstractMatrix{T}}, ::BandedColumns) where T
    iszero(f(x, zero(T))) || return _default_banded_broadcast(broadcasted(f, x,src))
    _BandedMatrix(f.(x, bandeddata(src)), axes(src,1), bandwidths(src)...)
end

function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix}})
    (A,) = bc.args
    _banded_broadcast(bc.f, A, MemoryLayout(typeof(A)))
end

function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:Number}})
    (A,x) = bc.args
    _banded_broadcast(bc.f, (A, x), MemoryLayout(typeof(A)))
end

function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:Number,<:AbstractMatrix}})
    (x,A) = bc.args
    _banded_broadcast(bc.f, (x,A), MemoryLayout(typeof(A)))
end

function copy(bc::Broadcasted{BandedStyle, <:Any, <:Any, <:Tuple{<:AbstractMatrix,<:AbstractMatrix}})
    _banded_broadcast(bc.f, bc.args, MemoryLayout.(typeof.(bc.args)))
end

###
# broadcast bandwidths
###

_bcsize(A) = size(A)

_broadcast_bandwidths(bnds) = bnds
_broadcast_bandwidths(bnds, _) = bnds
_broadcast_bandwidths((l,u), a::AbstractVector) = (bandwidth(a,1),u)
function __broadcast_bandwidths((l,u), A)
    sz = _bcsize(A)
    (length(sz) == 1 || sz[2] == 1) && return (bandwidth(A,1),u)
    sz[1] == 1 && return (l, bandwidth(A,2))
    bandwidths(A) # need to special case vector broadcasting
end
_broadcast_bandwidths((l,u), A::AbstractArray) = __broadcast_bandwidths((l,u), A)
_broadcast_bandwidths((l,u), A::Broadcasted) = __broadcast_bandwidths((l,u), A)

_band_eval_args() = ()
_band_eval_args(a, b...) = (a, _band_eval_args(b...)...)
_band_eval_args(::Base.RefValue{Type{T}}, b...) where T = (T, _band_eval_args(b...)...)
_band_eval_args(a::AbstractMatrix{T}, b...) where T = (zero(T), _band_eval_args(b...)...)
_band_eval_args(a::AbstractVector{T}, b...) where T = (one(T), _band_eval_args(b...)...)
_band_eval_args(a::Broadcasted, b...) = (zero(Base.Broadcast.combine_eltypes(a.f, a.args)), _band_eval_args(b...)...)

# zero dominates. Take the minimum bandwidth. Always treat like matrix.
__bnds(m, n) = (m-1, n-1)
__bnds(m) = (m-1, 0)
_bnds(bc) = __bnds(size(bc)...)

bandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(*)}) =
    min.(_broadcast_bandwidths.(Ref(_bnds(bc)), bc.args)...)

bandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(/)}) = _broadcast_bandwidths(_bnds(bc), first(bc.args))
bandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(\)}) = _broadcast_bandwidths(_bnds(bc), last(bc.args))

# zero is preserved. Take the maximum bandwidth
_isweakzero(f, args...) =  iszero(f(_band_eval_args(args...)...))

function bandwidths(bc::Broadcasted)
    sz = _bcsize(bc)
    a,b = length(sz) == 1 ? (sz[1],1) : sz
    bnds = (a-1,b-1)
    _isweakzero(bc.f, bc.args...) && return min.(bnds, max.(_broadcast_bandwidths.(Ref(bnds), bc.args)...))
    bnds
end

banded_similar(T, sz, bnds) = BandedMatrix{T}(undef, sz, bnds)
similar(bc::Broadcasted{BandedStyle}, ::Type{T}) where T = banded_similar(T, size(bc), bandwidths(bc))

##
# lmul!/rmul!
##

function _banded_lmul!(α, A::AbstractMatrix, _)
    for j = rowsupport(A), k = colsupport(A,j)
        inbands_setindex!(A, α*inbands_getindex(A,k,j), k,j)
    end
    A
end

function _banded_rmul!(A::AbstractMatrix, α, _)
    for j = rowsupport(A), k = colsupport(A,j)
        inbands_setindex!(A, inbands_getindex(A,k,j)*α, k,j)
    end
    A
end

function _banded_lmul!(α::Number, A::AbstractMatrix, ::BandedColumns)
    lmul!(α, bandeddata(A))
    A
end

function _banded_rmul!(A::AbstractMatrix, α::Number, ::BandedColumns)
    rmul!(bandeddata(A), α)
    A
end

banded_lmul!(α::Number, A::AbstractMatrix) = _banded_lmul!(α, A, MemoryLayout(typeof(A)))
banded_rmul!(A::AbstractMatrix, α::Number) = _banded_rmul!(A, α, MemoryLayout(typeof(A)))

materialize!(L::Lmul{ScalarLayout,<:AbstractBandedLayout}) = banded_lmul!(L.A, L.B)
materialize!(L::Rmul{<:AbstractBandedLayout,ScalarLayout}) = banded_rmul!(L.A, L.B)

##
# axpy!
##

# these are the routines of the banded interface of other AbstractMatrices
banded_axpy!(a::Number, X::AbstractMatrix, Y::AbstractMatrix) = _banded_axpy!(a, X, Y, MemoryLayout(typeof(X)), MemoryLayout(typeof(Y)))
_banded_axpy!(a::Number, X::AbstractMatrix, Y::AbstractMatrix, ::BandedColumns, ::BandedColumns) =
    banded_generic_axpy!(a, X, Y)
_banded_axpy!(a::Number, X::AbstractMatrix, Y::AbstractMatrix, notbandedX, notbandedY) =
    banded_dense_axpy!(a, X, Y)

# additions and subtractions
@propagate_inbounds function banded_generic_axpy!(a::Number, X::AbstractMatrix, Y::AbstractMatrix)
    n,m = size(X)
    if (n,m) ≠ size(Y)
        throw(BoundsError())
    end
    Xl, Xu = bandwidths(X)
    Yl, Yu = bandwidths(Y)

    @boundscheck if Xl > Yl
        # test that all entries are zero in extra bands
        for j=1:size(X,2),k=max(1,j+Yl+1):min(j+Xl,n)
            if inbands_getindex(X, k, j) ≠ 0
                throw(BandError(X, (k,j)))
            end
        end
    end
    @boundscheck if Xu > Yu
        # test that all entries are zero in extra bands
        for j=1:size(X,2),k=max(1,j-Xu):min(j-Yu-1,n)
            if inbands_getindex(X, k, j) ≠ 0
                throw(BandError(X, (k,j)))
            end
        end
    end

    l = min(Xl,Yl)
    u = min(Xu,Yu)

    @inbounds for j=1:m,k=max(1,j-u):min(n,j+l)
        inbands_setindex!(Y, a*inbands_getindex(X,k,j) + inbands_getindex(Y,k,j) ,k, j)
    end
    Y
end

function banded_dense_axpy!(a::Number, X::AbstractMatrix, Y::AbstractMatrix)
    if size(X) != size(Y)
        throw(DimensionMismatch("+"))
    end
    @inbounds for j=1:size(X,2),k=colrange(X,j)
        Y[k,j] += a*inbands_getindex(X,k,j)
    end
    Y
end

##
# special cases to hack tests where 2.0 .* A .+ B can be constructed in different ways depending on REPL or in a function.
# Should be made more general.
##
function copyto!(dest::AbstractArray{T}, bc::Broadcasted{BandedStyle, <:Any, typeof(+),
                                                        <:Tuple{<:Broadcasted{BandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                                                        <:AbstractMatrix}}) where T
    αA,B = bc.args
    α,A = αA.args
    dest ≡ B || (dest .= B)
    banded_axpy!(α, A, dest)
end

function similar(bc::Broadcasted{BandedStyle, <:Any, typeof(+),
                        <:Tuple{<:Broadcasted{BandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                        <:AbstractMatrix}}, ::Type{T}) where T
    αA,B = bc.args
    α,A = αA.args
    n,m = size(A)
    (n,m) == size(B) || throw(DimensionMismatch())
    Al,Au = bandwidths(A)
    Bl,Bu = bandwidths(B)
    similar(A, T, n, m, max(Al,Bl), max(Au,Bu))
end

axpy!(α, A::AbstractBandedMatrix, dest::AbstractMatrix) = banded_axpy!(α, A, dest)

###
# Fill
###

for op in (:*, :/)
    @eval begin
        broadcasted(::BandedStyle, ::typeof($op), a::Zeros, b::AbstractArray) = _broadcasted_zeros(a, b)
        broadcasted(::BandedStyle, ::typeof($op), a::Ones{T}, b::AbstractArray{V}) where {T,V} = copy_oftype(b, promote_op(*, T, V))
    end
end

for op in (:*, :\)
    @eval begin
        broadcasted(::BandedStyle, ::typeof($op), a::AbstractArray, b::Zeros) = _broadcasted_zeros(a, b)
        broadcasted(::BandedStyle, ::typeof($op), a::AbstractArray{T}, b::Ones{V}) where {T,V} = copy_oftype(a, promote_op(*, T, V))
    end
end

ext: jl


text:
# includet("D:\\projects\\vsCode\\MedPipe\\MedPipe\\src\\structs\\distinctColorsSaved.jl")
# includet("D:\\projects\\vsCode\\MedPipe\\MedPipe\\src\\visualizationUtils\\visualizationFromHdf5.jl")
# includet("D:\\projects\\vsCode\\MedPipe\\MedPipe\\src\\HDF5Utils\\HDF5saveUtils.jl")
# includet("D:\\projects\\vsCode\\MedPipe\\MedPipe\\src\\fromMonai\\LoadFromMonai.jl")


text:
module TestTextAnalysis
    using Test
    using SnowballStemmer

    my_tests = [
        "stemmer.jl",
    ]

    println("Running tests:")

    for my_test in my_tests
        println(" * $(my_test)")
        include(my_test)
    end
end


text:
# example data (udgrade to nside 256 of the Planck 2018 maps)

function planck256_mapdir()
    rootpath = artifact"planck_pr3_lowres_maps"
    return rootpath
end

function planck256_beamdir()
    rootpath = artifact"planck_pr3_beams"
    return rootpath
end

function planck_FITS_name_to_col(name_)
    name = String(name_)
    (name == "I_STOKES") && return 1
    (name == "Q_STOKES") && return 2
    (name == "U_STOKES") && return 3
    (name == "HITS")     && return 4
    (name == "II_COV")   && return 5
    (name == "IQ_COV")   && return 6
    (name == "IU_COV")   && return 7
    (name == "QQ_COV")   && return 8
    (name == "QU_COV")   && return 9
    (name == "UU_COV")   && return 10
end

"""
    planck256_map(freq, split, col, type::Type=Float64) -> HealpixMap{T, RingOrder}

Returns a Planck 2018 half-mission frequency map downgraded to nside 256 in KCMB units.
FITS file column numbers are

1. I_STOKES
2. Q_STOKES
3. U_STOKES
4. HITS
5. II_COV
6. IQ_COV
7. IU_COV
8. QQ_COV
9. QU_COV
10. UU_COV

# Arguments:
- `freq::String`: Planck frequency ∈ {"100", "143", "217"}
- `split::String`: half mission split ∈ {"hm1", "hm2"}
- `col::Int`: FITS file column. Either a number, String, or Symbol above.

# Returns:
- `Map{T, RingOrder}`: the map
"""
function planck256_map(freq::String, split, col, T::Type=Float64)
    if typeof(col) != Int
        col = planck_FITS_name_to_col(col)
    end
    @assert split[1:2] == "hm"
    rootpath = planck256_mapdir()
    fname = "nside256_HFI_SkyMap_$(freq)_2048_R3.01_halfmission-$(split[3]).fits"
    return nest2ring(readMapFromFITS(joinpath(rootpath, "plancklowres", fname), col, T))
end

function planck256_polmap(freq, split, T::Type=Float64)
    return PolarizedHealpixMap(
        planck256_map(freq, split, 1, T),
        planck256_map(freq, split, 2, T),
        planck256_map(freq, split, 3, T))
end

"""
    planck256_mask(freq, split, maptype, T::Type=Float64) -> HealpixMap{T}

# Arguments:
- `freq::String`: Planck frequency ∈ {"100", "143", "217"}
- `split::String`: half mission split ∈ {"hm1", "hm2"}
- `maptype`: pass T or P, as a String or Symbol

# Returns:
- `HealpixMap{T, RingOrder}`: the mask
"""
function planck256_mask(freq, split, maptype, T::Type=Float64)
    maptype_str = String(maptype)
    if maptype_str == "T"
        map_str = "temperature"
    elseif maptype_str == "P"
        map_str = "polarization"
    else
        throw(ArgumentError("maptype must be either \"T\" or \"P\""))
    end

    rootpath = planck256_mapdir()
    fname = "nside256_COM_Mask_Likelihood-$(map_str)-$(freq)-$(split)_2048_R3.00.fits"
    return readMapFromFITS(joinpath(rootpath, "plancklowres", fname), 1, T)
end

function planck_beam_Wl(T::Type, freq1, split1, freq2, split2, spec1_, spec2_;
                          lmax=4000)
    spec1 = String(spec1_)
    spec2 = String(spec2_)

    rootpath = planck256_beamdir()

    if parse(Int, freq1) > parse(Int, freq2)
        freq1, freq2 = freq2, freq1
        split1, split2 = split2, split1
    end
    if (parse(Int, freq1) == parse(Int, freq2)) && ((split1 == "hm2") && (split1 == "hm1"))
        split1, split2 = split2, split1
    end

    fname = "Wl_R3.01_plikmask_$(freq1)$(split1)x$(freq2)$(split2).fits"
    f = FITS(joinpath(rootpath, "planckbeams_plikmask", fname))
    bl = convert(Vector{T}, read(f[spec1], "$(spec1)_2_$(spec2)")[:,1])
    if lmax < 4000
        bl = bl[1:lmax+1]
    else
        bl = vcat(bl, last(bl) * ones(T, lmax - 4000))
    end
    return SpectralVector(bl)
end
planck_beam_Wl(T::Type, freq1, split1, freq2, split2, spec1; kwargs...) =
    lanck_beam_bl(T, freq1, split1, freq2, split2, spec1, spec1; kwargs...)
planck_beam_Wl(freq1::String, split1, freq2, split2, spec1, spec2; kwargs...) =
    planck_beam_Wl(Float64, freq1, split1, freq2, split2, spec1, spec2; kwargs...)

function planck_theory_Dl()
    rootpath = artifact"planck_pr3_theory"
    fname = "COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt"
    table = read_commented_header(joinpath(rootpath, fname))

    syms = (:TT, :TE, :EE, :BB, :PP)
    return NamedTuple{syms}(map(
        XY -> SpectralVector([0.0, 0.0, table[!, XY]...]), syms))
end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
module DifferentiableNAS

using Flux
using Juno
using Base.Iterators
using StatsBase: mean
using Zygote
using LinearAlgebra

include("DARTSModel.jl")

include("DARTSTraining.jl")
include("ActivationTraining.jl")
include("MaskedTraining.jl")
include("ScalingTraining.jl")
include("ADMMTraining.jl")

include("utils.jl")

end # module

max_stars_repo_name: SuReLI/DifferentiableNAS

max_issues_repo_name: SuReLI/DifferentiableNAS

max_forks_repo_name: SuReLI/DifferentiableNAS


text:
function solve_tsp(tsp::TSP, parameters=AlgorithmParameters(); verbose=true, use_dist_mtx=false)
    n = tsp.dimension
    x = tsp.nodes[:, 1]
    y = tsp.nodes[:, 2]
    serv_time = zeros(size(x))
    dem = ones(size(x))
    dem[1] = 0.0

    vehicleCapacity = n
    n_vehicles = 1

    isRoundingInteger = true
    isDurationConstraint = false
    duration_limit = C_DBL_MAX

    if use_dist_mtx
        c_dist_mtx = Matrix(tsp.weights')
        for i in 1:n
            c_dist_mtx[i, i] = 0.0
        end
        # need to input dist_mtx' instead of dist_mtx
        # Julia: column-first indexing
        # C: row-first indexing
        return c_api_solve_cvrp_dist_mtx(
            n, x, y, c_dist_mtx, serv_time, dem, vehicleCapacity, duration_limit, isDurationConstraint,
            n_vehicles, parameters, verbose
        )
    else
        return c_api_solve_cvrp(
            n, x, y, serv_time, dem,
            vehicleCapacity, duration_limit, isRoundingInteger, isDurationConstraint, n_vehicles, parameters, verbose
        )
    end
end

function solve_tsp(dist_mtx::Matrix, parameters=AlgorithmParameters(); verbose=true, x_coordinates=Float64[], y_coordinates=Float64[])
    n = size(dist_mtx, 1)

    serv_time = zeros(n)
    dem = ones(n)
    dem[1] = 0.0

    vehicleCapacity = n
    n_vehicles = 1

    c_dist_mtx = Matrix(dist_mtx')
    for i in 1:n
        c_dist_mtx[i, i] = 0.0
    end

    duration_limit = C_DBL_MAX
    isDurationConstraint = false

    # need to input dist_mtx' instead of dist_mtx
    # Julia: column-first indexing
    # C: row-first indexing
    return c_api_solve_cvrp_dist_mtx(
        n, x_coordinates, y_coordinates, c_dist_mtx, serv_time, dem,
        vehicleCapacity, duration_limit, isDurationConstraint,
        n_vehicles, parameters, verbose
    )
end
function solve_tsp(x::Vector, y:: Vector, parameters=AlgorithmParameters(); verbose=true)
    n = length(x)
    serv_time = zeros(n)
    dem = ones(n)
    dem[1] = 0.0
    vehicleCapacity = n
    n_vehicles = 1

    isRoundingInteger = true
    duration_limit = C_DBL_MAX
    isDurationConstraint = false

    return c_api_solve_cvrp(
        n, x, y, serv_time, dem,
        vehicleCapacity, duration_limit, isRoundingInteger, isDurationConstraint,
        n_vehicles, parameters, verbose
    )
end


text:
using DFTK
using Unitful
using UnitfulAtomic
using AtomsBase
using Test
import DFTK: compute_pspmap

@testset "DFTK -> AbstractSystem -> DFTK" begin
    Si = ElementCoulomb(:Si)
    C  = ElementPsp(:C, psp=load_psp("hgh/pbe/c-q4.hgh"))
    H  = ElementPsp(:H, psp=load_psp("hgh/lda/h-q1.hgh"))

    lattice   = randn(3, 3)
    atoms     = [Si, C, H, C]
    positions = [rand(3) for _ in 1:4]
    magnetic_moments = DFTK.normalize_magnetic_moment.(rand(4))

    system = DFTK.construct_system(lattice, atoms, positions, magnetic_moments)
    @test atomic_symbol(system) == [:Si, :C, :H, :C]
    @test bounding_box(system)  == collect(eachcol(lattice)) * u"bohr"
    @test position(system)      == [lattice * p * u"bohr" for p in positions]

    @test system[1].potential == Si
    @test system[2].potential == C
    @test system[3].potential == H
    @test system[4].potential == C
    @test !hasproperty(system[1], :pseudopotential)
    @test system[2].pseudopotential == "hgh/pbe/c-q4.hgh"
    @test system[3].pseudopotential == "hgh/lda/h-q1.hgh"
    @test system[4].pseudopotential == "hgh/pbe/c-q4.hgh"

    parsed = DFTK.parse_system(system)
    @test parsed.lattice   == lattice
    @test parsed.atoms     == atoms
    @test parsed.positions ≈ positions atol=1e-14
    @test parsed.magnetic_moments == magnetic_moments

    let system = attach_psp(system; family="hgh", functional="lda")
        for i in 1:4
            @test system[i].potential       == atoms[i]
            @test system[i].magnetic_moment == magnetic_moments[i]
        end
    end

    for constructor in (Model, model_atomic, model_LDA, model_PBE, model_SCAN)
        model = constructor(system)
        @test model.spin_polarization == :collinear
        newsys = periodic_system(model, magnetic_moments)

        @test atomic_symbol(system)       == atomic_symbol(newsys)
        @test bounding_box(system)        == bounding_box(newsys)
        @test boundary_conditions(system) == boundary_conditions(newsys)
        @test maximum(maximum, position(system) - position(newsys)) < 1e-12u"bohr"
        for (p, newp) in zip(system, newsys)
            @test p.magnetic_moment == newp.magnetic_moment
            @test p.potential       == newp.potential
        end
    end
end

@testset "AbstractSystem -> DFTK" begin
    lattice     = [12u"bohr" * rand(3) for _ in 1:3]
    atoms       = [:C => rand(3), :Si => rand(3), :H => rand(3), :C => rand(3)]
    pos_units   = last.(atoms)
    pos_lattice = austrip.(reduce(hcat, lattice))
    system      = periodic_system(atoms, lattice; fractional=true)

    let model = Model(system)
        @test model.lattice == pos_lattice
        @test model.positions ≈ pos_units atol=1e-14
        @test model.spin_polarization == :none

        @test length(model.atoms) == 4
        @test model.atoms[1] == ElementCoulomb(:C)
        @test model.atoms[2] == ElementCoulomb(:Si)
        @test model.atoms[3] == ElementCoulomb(:H)
        @test model.atoms[4] == ElementCoulomb(:C)
    end

    let system = attach_psp(system; family="hgh", functional="pbe")
        @test length(system) == 4
        @test system[1].pseudopotential == "hgh/pbe/c-q4.hgh"
        @test system[2].pseudopotential == "hgh/pbe/si-q4.hgh"
        @test system[3].pseudopotential == "hgh/pbe/h-q1.hgh"
        @test system[4].pseudopotential == "hgh/pbe/c-q4.hgh"

        parsed = DFTK.parse_system(system)
        @test parsed.lattice == pos_lattice
        @test parsed.positions ≈ pos_units atol=1e-14
        @test isempty(parsed.magnetic_moments)

        @test length(parsed.atoms) == 4
        @test parsed.atoms[1].psp.identifier == "hgh/pbe/c-q4.hgh"
        @test parsed.atoms[2].psp.identifier == "hgh/pbe/si-q4.hgh"
        @test parsed.atoms[3].psp.identifier == "hgh/pbe/h-q1.hgh"
        @test parsed.atoms[4].psp.identifier == "hgh/pbe/c-q4.hgh"
    end

    let system = attach_psp(system; family="hgh", functional="lda")
        @test length(system) == 4
        @test system[1].pseudopotential == "hgh/lda/c-q4.hgh"
        @test system[2].pseudopotential == "hgh/lda/si-q4.hgh"
        @test system[3].pseudopotential == "hgh/lda/h-q1.hgh"
        @test system[4].pseudopotential == "hgh/lda/c-q4.hgh"

        model = Model(system)
        @test model.lattice == pos_lattice
        @test model.positions ≈ pos_units atol=1e-14
        @test model.spin_polarization == :none

        @test length(model.atoms) == 4
        @test model.atoms[1].psp.identifier == "hgh/lda/c-q4.hgh"
        @test model.atoms[2].psp.identifier == "hgh/lda/si-q4.hgh"
        @test model.atoms[3].psp.identifier == "hgh/lda/h-q1.hgh"
        @test model.atoms[4].psp.identifier == "hgh/lda/c-q4.hgh"
    end
end

@testset "Check attach_psp routine selectively" begin
    symbols = [:Cu, :Au, :Ni]

    let pspmap = compute_pspmap(symbols; functional="lda", family="hgh", core=:semicore)
        @test pspmap[:Cu] == "hgh/lda/cu-q19.hgh"
        @test pspmap[:Au] == "hgh/lda/au-q19.hgh"
        @test pspmap[:Ni] == "hgh/lda/ni-q18.hgh"
    end

    let pspmap = compute_pspmap(symbols; functional="lda", family="hgh", core=:fullcore)
        @test pspmap[:Cu] == "hgh/lda/cu-q11.hgh"
        @test pspmap[:Au] == "hgh/lda/au-q11.hgh"
        @test pspmap[:Ni] == "hgh/lda/ni-q10.hgh"
    end

    let pspmap = compute_pspmap([:Cu, :Au]; functional="pbe", family="hgh", core=:semicore)
        @test pspmap[:Cu] == "hgh/pbe/cu-q19.hgh"
        @test pspmap[:Au] == "hgh/pbe/au-q19.hgh"
    end

    let pspmap = compute_pspmap(symbols; functional="pbe", family="hgh", core=:fullcore)
        @test pspmap[:Cu] == "hgh/pbe/cu-q11.hgh"
        @test pspmap[:Au] == "hgh/pbe/au-q11.hgh"
        @test pspmap[:Ni] == "hgh/pbe/ni-q18.hgh"
    end
end


text:
function transformPolyline(value, index)
end

decodePolyline(str::AbstractString; kwargs...) = decodePolyline(Vector{UInt8}(str); kwargs...)
function decodePolyline(polyline::Vector{UInt8}; precision=5)
	nums = Int32[0]
	count = 0
	for chunk = polyline .- 0x3f

		# accumulate bit chunks starting at LSB
		nums[end] |= eltype(nums)(chunk & 0x1f) << (precision*count)

		# the bit string continues if this bit is set, otherwise move to next
		if chunk & 0x20 == 0x20
			count += 1
		else
			nums[end] = (isodd(nums[end]) ? ~nums[end] : nums[end]) >> 1;
			count = 0
			push!(nums, 0)
		end
	end

	# rescale and accumulate the flat list of numbers
	return cumsum(reshape(nums[1:end-1]./10^precision, 2, :); dims=2)'
end


text:
using Test
using JSON3
using Pinecone

APIKEY = ENV["PINECONE_API_KEY"]
CLOUDENV="us-west1-gcp"

context = Pinecone.init(APIKEY, CLOUDENV)

INDEX = "filter-example"
NAMESPACE = "mynamespace"
index = Pinecone.Index(INDEX);

testdict = Dict{String, Any}("genre"=>"documentary", "year"=>2019);
v1 = [0.3,0.11,0.3,0.3,0.3,0.3,0.3,0.3,0.4,0.3]
v2 = [0.9, 0.8, 0.7, 0.6, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]

@testset "PineconeIndex Type Tests" begin
    pineconeindex = PineconeIndex(TESTINDEX)
    @test typeof(pineconeindex) == PineconeIndex
    @test pineconeindex.indexname == TESTINDEX
    pineconeindex = Pinecone.Index(TESTINDEX)
    @test typeof(pineconeindex) == PineconeIndex
    @test pineconeindex.indexname == TESTINDEX
    println("PineconeIndex is ", pineconeindex)
 end;

 @testset "PineconeContext Type Tests" begin
    KEY = "asdf-asdfadfs-adsfafds"
    ENV = "us.east.gcp"
    PROJECT = "testproject"
    pineconecontext = PineconeContext(KEY, ENV, PROJECT)
    #testing print out for codecov
    println("pineconecontext is ", pineconecontext)
    @test typeof(pineconecontext) == PineconeContext
    @test pineconecontext.apikey == KEY
    @test pineconecontext.cloudenvironment == ENV
    @test pineconecontext.projectname == PROJECT
 end

 @testset "PineconeVector Type Tests" begin
    VID = "vectorid"
    VALUES = [0.1, 0.2, 0.3, 0.4]
    META = Dict{String, Any}("foo"=>"bar", "key"=>"value", "numeric"=>1234)
    pineconevector = PineconeVector(VID, VALUES, META)
    @test typeof(pineconevector) == PineconeVector
    @test pineconevector.id == VID
    @test length(pineconevector.values) == 4
    @test length(pineconevector.metadata) == 3
    @test typeof(pineconevector.metadata) == Dict{String,Any}
    #testing print out for codecov
    println("pineconevector is ", pineconevector)
 end

 @testset "Test Init" begin
    #make a bogus init() call, should return nothing
    BADAPIKEY = "zzzz123456"
    context = Pinecone.init(BADAPIKEY, "badcloud")

    @test context == nothing

    #now check for proper setup
    goodcontext = Pinecone.init(GOODAPIKEY, CLOUDENV)
    @test goodcontext != nothing
    @test typeof(goodcontext) == PineconeContext
 end


text:
const ComplexExpression = Union{Expression2{typeof(+),E1,E2} where {E1,E2},Expression2{typeof(-),E1,E2} where {E1,E2}}

function stringsum(e::ExpressionSum{E,I}) where {E,I}
    res = string(e.inner)
    @simd for i in eachindex(e.es)
        res *= ("+" * string(e.es[i]))
    end
    res
end

function stringsum(e::ExpressionSum{E,Nothing}) where E
    res = string(e.es[1])
    @simd for i in 2:length(e.es)
        res *= ("+" * string(e.es[i]))
    end
    res
end

const printlist = [
    (:(Source{Parameter}),e->"p"),
    (:(Source{Variable}),e->"x"),
    (:(Variable),e-> "x[$(index(e))]"),
    (:(Parameter),e-> "p[$(index(e))]"),
    (:(Constant),e-> "$(refval(e))"),
    (:(Expression2{typeof(+),E1,E2} where {E1,E2}),e-> "$(string(e.e1)) + $(string(e.e2))"),
    (:(Expression2{typeof(-),E1,E2} where {E1,E2}),e-> "$(string(e.e1)) - $(string(e.e2))"),
    (:(Expression2{typeof(*),E1,E2} where {E1,E2}),e-> "$(string(e.e1))*$(string(e.e2))"),
    (:(Expression2{typeof(/),E1,E2} where {E1,E2}),e-> "$(string(e.e1))/$(string(e.e2))"),
    (:(Expression2{typeof(^),E1,E2} where {E1,E2}),e-> "$(string(e.e1))^$(string(e.e2))"),
    (:(Expression2{typeof(+),E1,E2} where {E1<:ComplexExpression,E2}),e-> "$(string(e.e1)) + $(string(e.e2))"),
    (:(Expression2{typeof(-),E1,E2} where {E1<:ComplexExpression,E2}),e-> "$(string(e.e1)) - $(string(e.e2))"),
    (:(Expression2{typeof(*),E1,E2} where {E1<:ComplexExpression,E2}),e-> "($(string(e.e1)))*$(string(e.e2))"),
    (:(Expression2{typeof(/),E1,E2} where {E1<:ComplexExpression,E2}),e-> "($(string(e.e1)))/$(string(e.e2))"),
    (:(Expression2{typeof(^),E1,E2} where {E1<:ComplexExpression,E2}),e-> "($(string(e.e1)))^$(string(e.e2))"),
    (:(Expression2{typeof(+),E1,E2} where {E1,E2<:ComplexExpression}),e-> "$(string(e.e1)) + ($(string(e.e2)))"),
    (:(Expression2{typeof(-),E1,E2} where {E1,E2<:ComplexExpression}),e-> "$(string(e.e1)) - ($(string(e.e2)))"),
    (:(Expression2{typeof(*),E1,E2} where {E1,E2<:ComplexExpression}),e-> "$(string(e.e1))*($(string(e.e2)))"),
    (:(Expression2{typeof(/),E1,E2} where {E1,E2<:ComplexExpression}),e-> "$(string(e.e1))/($(string(e.e2)))"),
    (:(Expression2{typeof(^),E1,E2} where {E1,E2<:ComplexExpression}),e-> "$(string(e.e1))^($(string(e.e2)))"),
    (:(Expression2{typeof(+),E1,E2} where {E1<:ComplexExpression,E2<:ComplexExpression}),e-> "$(string(e.e1)) + $(string(e.e2))"),
    (:(Expression2{typeof(-),E1,E2} where {E1<:ComplexExpression,E2<:ComplexExpression}),e-> "$(string(e.e1)) - ($(string(e.e2)))"),
    (:(Expression2{typeof(*),E1,E2} where {E1<:ComplexExpression,E2<:ComplexExpression}),e-> "($(string(e.e1)))*($(string(e.e2)))"),
    (:(Expression2{typeof(/),E1,E2} where {E1<:ComplexExpression,E2<:ComplexExpression}),e-> "($(string(e.e1)))/($(string(e.e2)))"),
    (:(Expression2{typeof(^),E1,E2} where {E1<:ComplexExpression,E2<:ComplexExpression}),e-> "($(string(e.e1)))^($(string(e.e2)))"),
    (:(ExpressionSum{E,I} where {E,I}),stringsum),
]

for (f,df,ddf) in f_nargs_1
    @eval begin
        push!(printlist,(:(Expression1{$(typeof($f)),E} where E),e->"$(string($f))($(string(e.e1)))"))
    end
end

for (f,df1,df2,ddf11,ddf12,ddf22) in f_nargs_2
    f in [:+,:-,:*,:^,:/] && continue
    @eval begin
        push!(printlist,(:(Expression2{$(typeof($f)),E1,E2} where {E1,E2}),e->"$(string($f))($(string(e.e1)),$(string(e.e2)))"))
    end
end

for (T,f) in printlist
    @eval begin
        print(io::IO,e::$T) = print(io, string(e))
        show(io::IO,::MIME"text/plain",e::$T) = print(io,e)
        string(e::$T) = $f(e)
    end
end

# (:(Model),e-> "NLP model with $(m.n) variables and $(m.m) constraints"


text:
import CoverageTools

const LCOV = CoverageTools.LCOV
const readfile = CoverageTools.LCOV.readfile
const writefile = CoverageTools.LCOV.writefile


text:
module Parallel

using LightGraphs
using LightGraphs: sample, AbstractPathState, JohnsonState, BellmanFordState, FloydWarshallState
using Distributed: @distributed
using Base.Threads: @threads, nthreads, Atomic, atomic_add!, atomic_cas!
using SharedArrays: SharedMatrix, SharedVector, sdata
using Arpack: eigs
using Random:shuffle
import SparseArrays: sparse
import Base: push!, popfirst!, isempty, getindex

include("shortestpaths/bellman-ford.jl")
include("shortestpaths/dijkstra.jl")
include("shortestpaths/floyd-warshall.jl")
include("shortestpaths/johnson.jl")
include("centrality/betweenness.jl")
include("centrality/closeness.jl")
include("centrality/pagerank.jl")
include("centrality/radiality.jl")
include("centrality/stress.jl")
include("distance.jl")
include("traversals/bfs.jl")
include("traversals/greedy_color.jl")
end


text:
using Plots
using ColorSchemes
using Combinatorics: combinations
"""
This file calculates exact burden and fake burden.
"""

"""
    grey_next_flip(n)

A Grey code is a map from a binary number to a different binary number.
The magic is that if you call this `n=2^i` times, you will see every binary
number between `0` and `n-1`, but each one will be only one digit off from
the last. This function, in particular, tells you which bit will flip next.
"""
function grey_next_flip(n)
    for i in 1:1:64
        if n & 1 == 1
            return i
        else
            n >>= 1
        end
    end
    65
end

"""
Calculates (b_j /sum(b_i))(1 - prod(1-b_i)) (prod(p_i) prod(1 - p_j)).
"""
function exact_burden_term(weights, prevalences, which)
    w = weights[which]
    yes = prevalences[which]
    no = prevalences[.!which]
    results = zeros(Float64, length(which))
    results[which] += w * ((1 - prod(1 .- w)) * prod(yes) * prod(1 .- no) / sum(w))
    results
end

"""
Calculates (b_j /sum(b_i))(1 - prod(1-b_i)) (prod(p_i) prod(1 - p_j))
and adds it to the given running sum.
"""
function exact_burden_term!(weights, prevalences, which, running_sum)
    w = weights[which]
    yes = prevalences[which]
    no = prevalences[.!which]
    running_sum[which] += w * ((1 - prod(1 .- w)) * prod(yes) * prod(1 .- no) / sum(w))
end

function exact_burden(weights, prevalences)
    n = length(prevalences)
    which = zeros(Bool, n)
    b = zeros(Float64, n)
    for i in 1:(1<<n - 1)
        to_flip = grey_next_flip(i)
        which[to_flip] = !which[to_flip]
        exact_burden_term!(weights, prevalences, which, b)
    end
    b
end

"""
The burden if we allow up to `n` comorbidities. `n=1` means you can have
one disease at a time.
"""
function order_n_burden(weights, prevalences, n)
    cnt = length(weights)
    burden = similar(prevalences)
    fill!(burden, 0)
    which = zeros(Bool, cnt)
    for depth in 1:n
        for como_idx in combinations(1:cnt, depth)
            fill!(which, false)
            which[como_idx] .= true
            exact_burden_term!(weights, prevalences, which, burden)
        end
    end
    burden
end

"""
A term-by-term gradient calculation with respect to index i
"""
function partial_term!(weights, prevalences, which, running_sum, i)
    w = weights[which]
    p = copy(prevalences)
    p[i] = which[i] ? 1 : 2  # 1 - no = -1, which is the derivative.
    yes = p[which]
    no = p[.!which]
    running_sum[which] += w * ((1 - prod(1 .- w)) * prod(yes) * prod(1 .- no) / sum(w))
end

function order_n_partial(weights, prevalences, n)
    cnt = length(weights)
    partial = zeros(Float64, cnt, cnt)
    which = zeros(Bool, cnt)
    for partial_wrt_j in 1:cnt
        burden = zeros(Float64, cnt)
        fill!(burden, 0)
        # Loop over how many comorbidities a subpopulation has.
        for depth in 1:n
            # Loop over different choices of comorbidities for each subpopulation.
            for como_idx in combinations(1:cnt, depth - 1)
                fill!(which, false)
                which[como_idx] .= true
                partial_term!(weights, prevalences, which, burden, partial_wrt_j)
            end
        end
        partial[:, partial_wrt_j] .= burden
    end
    partial
end

"""
Make random prevalences and burdens to calculate.
"""
function exact_burden_random(cnt)
    weights = rand(cnt)
    prevalences = rand(cnt)
    burden = exact_burden(weights, prevalences)
    weights, prevalences, burden
end

"""
This is the super simple, exact total burden for a population. Given how simple it is,
it's surprising that we don't have a closed form for individual contributions to burden.
"""
function total_burden(weights, prevalences)
    1 - prod(1 .- weights .* prevalences)
end

"""
Partial derivative of total comorbidity with respect to prevalence.
Writes into the `partials` buffer.
"""
function partial_wrt_prevalence!(partials, weights, prevalences)
    for i in 1:length(partials)
        save_p = prevalences[i]
        prevalences[i] = zero(eltype(prevalences))
        partials[i] = weights[i] * prod(1 .- weights .* prevalences)
        prevalences[i] = save_p
    end
end

"""
Maybe we can divide the total burden by some fraction in order to get
individual cause's burden. Nah.
"""
function fake_burden(weights, prevalences)
    total = total_burden(weights, prevalences)
    w = weights .* prevalences
    w * total / sum(w)
end

function fake3_burden(weights, prevalences)
    total = total_burden(weights, prevalences)
    w = weights .* prevalences
    rho = w .* (1 - total)
    rho * total / sum(rho)
end

"""
Another try at weighting the fake burden. Try calculating burden for n-1
and using those values as weights.
"""
function fake2_burden(weights, prevalences)
    total = total_burden(weights, prevalences)
    reduced = similar(weights)
    reduced_p = similar(prevalences)
    for ridx in 1:length(prevalences)
        reduced_p .= prevalences
        reduced_p[ridx] = 0.99 * prevalences[ridx]
        reduced[ridx] = 1 .- total_burden(weights, reduced_p)
    end
    reduced_product = weights .* prevalences .* reduced
    reduced_product * total / sum(reduced_product)
end

relerr(observed, expected) = (observed .- expected) ./ observed

function total_burden_random(cnt)
    w = rand(cnt)
    p = rand(cnt)
    exact = exact_burden(w, p)
    [total_burden(w, p), sum(exact)]
end

"""
Error in the fake burden is above 10%, even for 20 causes.
"""
function fake_burden_random(cnt)
    w = rand(cnt)
    p = rand(cnt)
    exact = exact_burden(w, p)
    fake = fake_burden(w, p)
    [sum(fake), sum(exact), maximum(abs.((fake .- exact) ./ exact))]
end

function plot_compare(weight, prevalence, burden)
    total = sum(burden)
    estimated = (weight .* prevalence) * total / sum(weight .* prevalence)
    rel = relerr(estimated, burden)
    normalized = (rel .- minimum(rel)) / (maximum(rel) - minimum(rel))
    colors = get(ColorSchemes.buda, normalized)
    scatter(weight, prevalence, markercolor = colors)
end

"""
Error in the fake burden is above 10%, even for `cnt` causes.
"""
function fake_burden_basic_plot(fake_func, cnt)
    plotter = scatter
    plot_cnt = 5
    for i in 1:plot_cnt
        w = rand(cnt)
        p = rand(cnt)
        exact = exact_burden(w, p)
        fake = fake_func(w, p)
        rplot = plotter(
                exact,
                fake,
                markercolor = get(ColorSchemes.tab20b, (i-1)/(plot_cnt - 1)),
                legend = false
                )
        plotter = scatter!
        if i == plot_cnt
            display(rplot)
        end
    end
end

function fake_gen(weight, prevalence)
    total = total_burden(weight, prevalence)
    estimated = (weight .* prevalence) * total / sum(weight .* prevalence)
end

function generate_prevalences(burden_gen, cause_cnt, trial_cnt)
    weight = rand(cause_cnt)
    relerr_and_prevalence = zeros(Float64, cause_cnt, 2, trial_cnt)
    for t in 1:trial_cnt
        prevalence = rand(cause_cnt)
        burden = exact_burden(weight, prevalence)
        estimated = burden_gen(weight, prevalence)
        rel = relerr(estimated, burden)
        relerr_and_prevalence[:, 1, t] = rel
        relerr_and_prevalence[:, 2, t] = prevalence
    end
    weight, relerr_and_prevalence
end

"""
This plot asks how the error of a single cause depends on its
prevalence. It plots relative error by prevalence. It shows that
larger prevalence is estimated below what it should be. There is scatter in
this plot. I'll bet that scatter is a function of the total prevalence for
all causes for each trial. Either as a sum or as it appears in the total burden
calculation.
"""
function plot_by_cause(weight, relerr_and_prevalence)
    cause_cnt = length(weight)
    trial_cnt = size(relerr_and_prevalence, 3)
    callit = scatter
    for cidx in 1:5:cause_cnt
        color = get(ColorSchemes.tab20b, (cidx - 1) / (cause_cnt - 1))
        rel = reshape(relerr_and_prevalence[cidx, 1, :], trial_cnt)
        prev = reshape(relerr_and_prevalence[cidx, 2, :], trial_cnt)
        display(callit(prev, rel, markercolor = color))
        callit = scatter!
    end
end

"""
This plots prevalence by weight. It colors by relative error. We see that
larger weights tend to skew towards positive relative error.
"""
function plot_trials(weight, relerr_and_prevalence)
    cause_cnt = length(weight)
    trial_cnt = size(relerr_and_prevalence, 3)
    w = repeat(weight, trial_cnt)
    rel = reshape(relerr_and_prevalence[:, 1, :], cause_cnt * trial_cnt)
    normrel = (rel .- minimum(rel)) / (maximum(rel) - minimum(rel))
    p = reshape(relerr_and_prevalence[:, 2, :], cause_cnt * trial_cnt)
    colors = get(ColorSchemes.viridis, normrel)

    @show length(w)
    @show length(rel)
    @show length(p)
    display(scatter(w, p, markercolor = colors))
end


text:
# Run tests for the NEP transformations

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@testset "transformations" begin

    @bench @testset "Shift and scale" begin
        orgnep=nep_gallery("dep0");

        # Test the shift by solving the orgnep and
        # check with residual of transformed nep.
        σ=-0.3+0.1im;
        nep1=shift_and_scale(orgnep,shift=σ);
        orgλ,orgv=augnewton(orgnep,v=ones(size(orgnep,1)))
        @test norm(compute_Mlincomb(nep1,orgλ-σ,orgv))<100*eps()
        λ1,v1=quasinewton(nep1,v=ones(size(orgnep,1)))
        @test abs(λ1+σ-orgλ)<eps()*100 # check that we get the same eigvals

        # Test shift and scaling
        σ=-0.4+0.01im; α=0.5
        nep2=shift_and_scale(orgnep,shift=σ,scale=α);
        λ2,v2=quasinewton(nep2,v=ones(size(nep2,1)))
        @test abs((α*λ2+σ)-orgλ)<eps()*1000

        # Check that PEP transformations correctly transform coefficients
        pep0=nep_gallery("pep0",10);
        σ=1+0.3im;
        α=3;
        pep1=shift_and_scale(pep0,shift=σ,scale=α)
        λ,v= quasinewton(pep0,λ=-0.3-0.7im,v=ones(size(pep0,1)));
        norm(compute_Mlincomb(pep0, λ, v))
        λv,V=polyeig(pep1);
        @test minimum(abs.(λv .- (λ-σ)/α)) < eps()*1000

        # Check that real PEP with real transformation is still real
        σ=3;
        α=1
        pep2=shift_and_scale(pep0,shift=σ,scale=α)
        @test !isa(Complex,eltype(pep2.A[1])) # Preserve realness

        nep3=nep_gallery("qdep0")
        λ,v= quasinewton(nep3,λ=1+1im,v=ones(size(nep3,1)));
        σ=-3+0.3im
        α=0.9;
        nep3_transf=shift_and_scale(nep3,shift=σ,scale=α);
        @test norm(compute_Mlincomb(nep3_transf,(λ-σ)/α,v))<sqrt(eps())*1000;
        λ,V=iar(nep3_transf,σ=0,neigs=2,maxit=60,v=ones(size(nep3_transf,1)))
        @test norm(compute_Mlincomb(nep3,α*λ[1]+σ,V[:,1]))<sqrt(eps())
        @test norm(compute_Mlincomb(nep3,α*λ[2]+σ,V[:,2]))<sqrt(eps())

        λ,V=iar(nep3_transf,σ=0,neigs=2,maxit=60,v=ones(size(nep3_transf,1)))
        @test norm(compute_Mlincomb(nep3,α*λ[1]+σ,V[:,1]))<sqrt(eps())
        @test norm(compute_Mlincomb(nep3,α*λ[2]+σ,V[:,2]))<sqrt(eps())
        #

    end
    @bench @testset "Möbius transformation" begin
        # Checks mobius_transformation

        pep0=nep_gallery("pep0")
        a=3+0.1im
        b=0.1im;
        c=1;
        d=1-0.3im
        pep0_transf=mobius_transform(pep0,a=a,b=b,c=c,d=d)
        λ,v= quasinewton(pep0_transf,λ=0.3im,v=ones(size(pep0,1)));
        λorg=(a*λ+b)/(c*λ+d)
        @test norm(compute_Mlincomb(pep0,λorg,v))<sqrt(eps());

        nep4=nep_gallery("qdep0")
        a=3-0.01im
        b=0.1im;
        c=1;
        d=1-0.3im
        nep4_transf=mobius_transform(nep4,a=a,b=b,c=c,d=d)
        @test isa(nep4_transf,SPMF_NEP)  # Möbius transformed SPMF-NEP is still a SPMF
        λ,v= quasinewton(nep4_transf,λ=1-0.3im,v=ones(size(nep4,1)));
        λorg=(a*λ+b)/(c*λ+d)
        @test norm(compute_Mlincomb(nep4,λorg,v))<sqrt(eps())*10;

        n=size(nep4,1);
        V=(1:n)*(1:5)'/n; # Arbitrary matrices
        S=(1:5)*(1:5)'/5-3I;
        s = Matrix(3.0*I, 1, 1)
        W1 = compute_MM(nep4, (c*S + d*I) \ (a*S + b*I), V)
        W2=compute_MM(nep4_transf,S,V);
        @test opnorm(W1-W2)<sqrt(eps())

        # Check that mobius_transformation becomes shift_and_scale
        c=0; d=1;
        nep4_transf1=mobius_transform(nep4,a=a,b=b,c=c,d=d)
        nep4_transf2=shift_and_scale(nep4,scale=a,shift=b);
        λ=3;
        w=compute_Mlincomb(nep4,λ*a+b,V[:,1])
        w1=compute_Mlincomb(nep4_transf1,λ,V[:,1])
        w2=compute_Mlincomb(nep4_transf2,λ,V[:,1])
        @test norm(w-w1)<eps()*100
        @test norm(w-w2)<eps()*100
        @test norm(w1-w2)<eps()*100
        W1=compute_Mlincomb(nep4_transf1,λ,V,[0,1,0,1,0])
        W2=compute_Mlincomb(nep4_transf2,λ,V,[0,1,0,1,0]);
        @test norm(W1-W2)<sqrt(eps())

    end

    @bench @testset "taylor exp" begin
        nep=nep_gallery("dep0");
        d=5;
        pep=taylor_expansion_pep(nep,d);
        s=0.1;
        @test norm(compute_Mder(pep,s)-compute_Mder(nep,s)) < abs(s)^d
    end

end


text:
# # Lorenz 63
#
#
# Data assimilation are numerical methods used in geosciences to
# mix the information of observations (noted as ``y``) and a dynamical
# model (noted as `f`) in order to estimate the true/hidden state of
# the system (noted as `x`) at every time step `k`. Usually, they are
# related following a nonlinear state-space model:
#
# ```math
# \begin{aligned}
# x_k & = f(x_{k-1}) + \eta_k \\
# y_k & = H x_k + \epsilon_k
# \end{aligned}
# ```
#
# with ``\eta`` and ``\epsilon`` some independant white
# Gaussian noises respectively respresenting the model forecast error
# and the error of observation.

# In classical data assimilation,
# we require multiple runs of an explicit dynamical model ``f`` with
# possible severe limitations including the computational cost, the
# lack of consistency of the model with respect to the observed data
# as well as modeling uncertainties. Here, an alternative strategy
# is explored by developing a fully data-driven assimilation. No
# explicit knowledge of the dynamical model is required. Only a
# representative catalog of trajectories of the system is assumed to
# be available. Based on this catalog, the Analog Data Assimilation
# (AnDA) is introduced by combining machine learning with the analog
# method (or nearest neighbor search) and stochastic assimilation
# techniques including Ensemble Kalman Filter and Smoother (AnEnKF,
# AnEnKS) and Particle Filter (AnPF). We test the accuracy of the technic
# on different chaotic dynamical models, the Lorenz-63 and Lorenz-96
# systems.  # # This Julia program is derived from the Python library
# is attached to the following publication: # Lguensat, R., Tandeo,
# P., Ailliot, P., Pulido, M., & Fablet, R. (2017). The Analog Data
# Assimilation. *Monthly Weather Review*, 145(10), 4093-4107.  # If
# you use this library, please do not forget to cite this work.

using Plots, DifferentialEquations, AnalogDataAssimilation

# To begin, as dynamical model ``f``, we use the Lorenz-63 chaotic
# system. First, we generate simulated trajectories from this dynamical
# model and store them into the catalog. Then, we use this catalog
# to emulate the dynamical model and we apply the analog data
# assimilation. Finally, we compare the results of this data-driven
# approach to the classical data assimilation (using the true Lorenz-63
# equations as dynamical model.

# ## Generate simulated data

σ = 10.0
ρ = 28.0
β = 8.0 / 3

dt_integration = 0.01
dt_states = 1
dt_obs = 8
parameters = [σ, ρ, β]
var_obs = [1]
nb_loop_train = 100
nb_loop_test = 10
sigma2_catalog = 0.0
sigma2_obs = 2.0

ssm = StateSpaceModel(
    lorenz63,
    dt_integration,
    dt_states,
    dt_obs,
    parameters,
    var_obs,
    nb_loop_train,
    nb_loop_test,
    sigma2_catalog,
    sigma2_obs,
)

# - compute ``u_0`` to be in the attractor space

u0 = [8.0; 0.0; 30.0]
tspan = (0.0, 5.0)
prob = ODEProblem(ssm.model, u0, tspan, parameters)
u0 = last(solve(prob, save_everystep = false))

xt, yo, catalog = generate_data(ssm, u0);

plot(xt.t, vcat(xt.u'...)[:, 1])
scatter!(yo.t, vcat(yo.u'...)[:, 1]; markersize = 2)

# ## Classical data assimilation

regression = :local_linear
sampling = :gaussian
k, np = 100, 500

DA = DataAssimilation(ssm, xt)
x̂_classical = forecast(DA, yo, AnEnKS(np), progress = false)
@time RMSE(xt, x̂_classical)

# ## Analog data assimilation

f = AnalogForecasting(k, xt, catalog; regression = regression, sampling = sampling)
DA = DataAssimilation(f, xt, ssm.sigma2_obs)
@time x̂_analog = forecast(DA, yo, AnEnKS(np), progress = false)
println(RMSE(xt, x̂_analog))

# ## Comparison between classical and analog data assimilation

plot(xt.t, xt[1], label = "true state")
plot!(xt.t, x̂_classical[1], label = "classical")
plot!(xt.t, x̂_analog[1], label = "analog")
scatter!(yo.t, yo[1]; markersize = 2, label = "observations")

# The results show that performances of the data-driven analog data
# assimilation are closed to those of the model-driven data assimilation.
# The error can be reduced by augmenting the size of the catalog
# "nb_loop_train".

# ## Remark
#
# Note that for all the previous experiments, we use the robust
# Ensemble Kalman Smoother (EnKS) with the increment or local linear
# regressions and the Gaussian sampling. If you want to have realistic
# state estimations, we preconize the use of the Particle Filter
# `PF` with the locally constant regression (regression
# = :locally_constant) and the multinomial sampler (sampling =
# :multinomial) with a large number of particles `np`. For more
# details about the different options, see the attached publication:
# Lguensat, R., Tandeo, P., Ailliot, P., Pulido, M., & Fablet, R.
# (2017). The Analog Data Assimilation. *Monthly Weather Review*,
# 145(10), 4093-4107.


text:
export Overwrite

"""
A permission overwrite.
More details [here](https://discordapp.com/developers/docs/resources/channel#overwrite-object).
"""
struct Overwrite <: DiscordObject
    id::Snowflake
    type::Int
    allow::String
    deny::String
end
@boilerplate Overwrite :constructors :docs :lower :merge :mock


text:
# ## Data

using Distributions
using Random

Random.seed!(1)
y = rand(Normal(0,1), 100);

# ## Model

using Turing

@model m8_4(y) = begin
    ## Can't really set a Uniform[-Inf,Inf] on σ
    α₁ ~ Uniform(-3000, 1000)
    α₂ ~ Uniform(-1000, 3000)
    σ ~ truncated(Cauchy(0,1), 0, Inf)

    y ~ Normal(α₁ + α₂, σ)
end

chns = sample(m8_4(y), NUTS(), 2000)

# ## Output

# \defaultoutput{}

# ## Original output

"""
         mean      sd     5.5%   94.5% n_eff Rhat
 a1    -861.15 558.17 -1841.89  -31.04     7 1.43
 a2     861.26 558.17    31.31 1842.00     7 1.43
 sigma    0.97   0.07     0.89    1.09     9 1.17
""";


text:
data = exp.(randn(200)/5)
hist(data; nbins=20)


text:
# Autogenerated wrapper script for CoinUtils_jll for x86_64-linux-musl-cxx03
export libCoinUtils

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("CoinUtils")
JLLWrappers.@declare_library_product(libCoinUtils, "libCoinUtils.so.3")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libCoinUtils,
        "lib/libCoinUtils.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CollegeStratBase
using Test

@testset "All" begin
    include("directories_test.jl");
    include("unit_conversion_test.jl");
    include("notation_test.jl");
    include("display_test.jl");
    include("helpers_test.jl");
    include("school_groups_test.jl");
end

# ---------------


text:
# SPDX-License-Identifier: X11
# 2020-07-30
# Judge Status Summary (200pt)

N = parse(Int,readline())
S = readlines()

O = zeros(Int,4)

for i in 1:N
  if S[i] == "AC"
    O[1] = O[1] + 1
  elseif S[i] == "WA"
    O[2] = O[2] + 1
  elseif S[i] == "TLE"
    O[3] = O[3] + 1
  elseif S[i] == "RE"
    O[4] = O[4] + 1
  end
end

println(string("AC x ", O[1]))
println(string("WA x ", O[2]))
println(string("TLE x ", O[3]))
println(string("RE x ", O[4]))

max_stars_repo_name: jthvai/competitive-programming

max_issues_repo_name: jthvai/competitive-programming

max_forks_repo_name: jthvai/competitive-programming


text:
module GellMannMatrices

using LinearAlgebra, SparseArrays
export gellmann

# ---------------------------------------------------------------------------------------- #
# direct sum (i.e., block-diagonalization of matrices As)
if VERSION < v"1.8"
    # work around for https://github.com/JuliaLang/julia/pull/44211 for Julia ≤ v1.7
    @inline ⊕(As...) = cat(As..., dims=(1,2))::typeof((.+)(As...))
else
    @inline ⊕(As...) = cat(As..., dims=Val((1,2)))
end

# https://en.wikipedia.org/wiki/Matrix_unit
unit(::Type{<:SparseMatrixCSC}, i, j, n, m) = sparse([i], [j], [1], n, m)
unit(T::Type{<:AbstractMatrix}, i, j, n, m) = (A = zeros(eltype(T), n, m); A[i,j] = one(eltype(T)); A)

# ---------------------------------------------------------------------------------------- #
# implementation following https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices

# off-diagonal matrices (f_ij on Wiki)
function offdiag_gellmann(T::Type{<:AbstractMatrix}, i, j, d)
    if i < j
        return unit(T, i, j, d, d) + unit(T, j, i, d, d)
    elseif i > j
        return -1im * (unit(T, j, i, d, d) - unit(T, i, j, d, d))
    else
        throw(DomainError((j, i), "i and j must differ"))
    end
end

# "diagonal"-like matrices (h_i on Wiki)
function diag_gellmann(T::Type{<:AbstractMatrix}, i, d; norm_identity::Bool=false)
    if i == 1
        return norm_identity ? T(I(d) / sqrt(d/2)) : T(I(d))
    elseif i < d
        return diag_gellmann(T, i, d-1) ⊕ zero(eltype(T))
    elseif i == d
        return sqrt(2/(d*(d-1))) * (T(I(d-1)) ⊕ convert(eltype(T), (1-d)))
    else
        throw(DomainError((i,d), "invalid combination of (i,d); must have i ≤ d and d ≥ 1"))
    end
end

# ---------------------------------------------------------------------------------------- #

"""
    gellmann([T::Type{<:AbstractMatrix}=Matrix{ComplexF64},] d::Integer; skip_identity=true)

Return the generalized Gell-Mann matrices in `d` dimensions as a vector of matrices of type
`T`.

`T` can be any mutable `AbstractMatrix`, including `SparseMatrixCSC`, with complex floating
point element type.

## Keyword arguments
- `skip_identity` (default, `true`): if `false`, the identity matrix of dimension `d` is
  included as the last element of the returned vector.
- `norm_identity` (default, `false`): if `true`, the elements of the identity are normalized
   by `sqrt(D/2)` to ensure the the orthogonality relation `tr(A'*A) = 2`.

## Examples
Compute the Pauli matrices ``σ_i`` (`d=2`) and Gell-Mann matrices ``Λ_i`` (`d=3`):
```jl
julia> σᵢ = gellmann(2)
julia> Λᵢ = gellmann(3)
```
"""
gellmann(d::Integer; kwargs...) = gellmann(Matrix{ComplexF64}, d; kwargs...)

function gellmann(T::Type{<:AbstractMatrix}, d::Integer;
                  skip_identity=true, norm_identity=false)
    d > 0 || throw(DomainError(d, "dimension must be a positive integer"))
    eltype(T) == complex(eltype(T)) || throw(DomainError(T, "the element type of T must be complex"))
    # we collect the matrices in the slightly odd-looking manner below in order to be
    # consistent with the conventional number of Pauli and Gell-Mann matrices.
    ms = Vector{T}(undef, d*d - skip_identity)
    k = 0
    for i in 2:d
        for  j in 1:i-1
            ms[k+=1] = offdiag_gellmann(T, j, i, d)
            ms[k+=1] = offdiag_gellmann(T, i, j, d)
        end
        ms[k+=1] = diag_gellmann(T, i, d)
    end
    skip_identity || (ms[k+=1] = diag_gellmann(T, 1, d; norm_identity))
    return ms
end

end


text:
module PowerSimulations

#################################################################################
# Exports

#Network Relevant Exports
export StandardPTDFForm
export CopperPlatePowerModel

#operation_models
export UnitCommitment
export EconomicDispatch
export SCEconomicDispatch
export OptimalPowerFlow

#functions
export solve_op_model!

#################################################################################
# Imports
import JuMP
import ParameterJuMP
#using TimeSeries
import PowerSystems
import PowerModels
import MathOptInterface
import DataFrames
import LinearAlgebra
#import LinearAlgebra.BLAS #needed for the simulation stage
import AxisArrays
import Dates

#################################################################################
# Type Alias From other Packages
const PM = PowerModels
const PSY = PowerSystems
const PSI = PowerSimulations
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const PJ = ParameterJuMP

#Type Alias for JuMP and PJ containers
const JuMPExpressionMatrix = Matrix{<:JuMP.AbstractJuMPScalar}
const PGAE{V} = PJ.ParametrizedGenericAffExpr{Float64,V} where V <: JuMP.AbstractVariableRef
const GAE{V} = JuMP.GenericAffExpr{Float64,V} where V <: JuMP.AbstractVariableRef
const JuMPAffineExpressionArray = Matrix{GAE{V}} where V <: JuMP.AbstractVariableRef

#Type Alias for Unions
const FixResource = Union{PSY.RenewableFix, PSY.HydroFix}

#Type Alias for long type signatures
const MinMax = NamedTuple{(:min, :max),Tuple{Float64,Float64}}
const NamedMinMax = Tuple{String, MinMax}
const UpDown = NamedTuple{(:up, :down),Tuple{Float64,Float64}}

#################################################################################
# Includes

#Abstract Models
include("network/network_models/networks.jl")
include("services/service_models/services.jl")

#Core Models and constructors
include("core/core_structs/device_model.jl")
include("core/core_structs/canonical_model.jl")
include("core/core_structs/service_model.jl")
include("core/core_structs/operation_model.jl")
include("core/device_constructor.jl")
include("core/canonical_constructor.jl")
include("core/operations_constructor.jl")
include("core/core_structs/results_model.jl")
#include("core/core_structs/simulation_model.jl")

#Device Modeling components
include("devices/device_models/common.jl")
include("devices/device_models/renewable_generation.jl")
include("devices/device_models/thermal_generation.jl")
include("devices/device_models/electric_loads.jl")
include("devices/device_models/branches.jl")
include("devices/device_models/storage.jl")
include("devices/device_models/hydro_generation.jl")

#Network models
include("network/network_models/copperplate_model.jl")
include("network/network_models/powermodels_interface.jl")
include("network/network_models/ptdf_model.jl")

#Device constructors
include("devices/device_constructors/device_constructors.jl")

#Network constructors
include("network/network_constructor.jl")

#Services Models
#include("service_models/reserves.jl")

#Services constructors
include("services/services_constructor.jl")

#Operational Model Constructors
include("operation_models/operation_models.jl")

#Utils
include("routines/solve_routines.jl")
include("routines/get_ini_cond.jl")
#include("routines/simulation_routines.jl")
#include("routines/device_retreval.jl")

#################################################################################
##### JuMP methods overloading
JuMP.Model(optimizer::Nothing; kwargs...) = JuMP.Model(kwargs...)

#################################################################################

end


text:
using TeaFiles

function _test_api(
    data::Vector{T},
    lower::Union{Nothing, TeaFiles.TimeLike},
    upper::Union{Nothing, TeaFiles.TimeLike}
) where T
    # Temporary file to use for experiments. We have to use a file to test memory mapping,
    # since this can't be used with a buffer.
    path = tempname()

    # TODO Also test writing items with different field offsets.
    # Write the data to a file by hand.
    metadata = create_metadata_julia_time(T)
    open(path; write=true) do io
        # Create the file and write the header and data to it.
        write(io, metadata)
        write(io, data)
    end

    # Write to a file with the table API
    path_table_api = tempname()
    table = _structs_to_namedtuples(data)
    TeaFiles.write(path_table_api, table)

    # Write to a buffer.
    buffer = IOBuffer()
    write(buffer, metadata)
    write(buffer, data)
    seekstart(buffer)

    result_file = TeaFiles.read(path; lower=lower, upper=upper)
    result_file_table_api = TeaFiles.read(path_table_api; lower=lower, upper=upper)
    result_buffer = TeaFiles.read(buffer; lower=lower, upper=upper)

    # Filter data by time field.
    expected_read_data = data
    i_time_field = nothing
    for i in 1:fieldcount(eltype(data))
        if fieldtype(eltype(data), i) == DateTime
            i_time_field = i
            break
        end
    end
    if !isnothing(lower)
        expected_read_data = [
            x for x in expected_read_data if getfield(x, i_time_field) >= lower
        ]
    end
    if !isnothing(upper)
        expected_read_data = [
            x for x in expected_read_data if getfield(x, i_time_field) < upper
        ]
    end

    namedtuple_data = _structs_to_namedtuples(expected_read_data)
    @test result_file == namedtuple_data
    @test result_file_table_api == namedtuple_data
    @test result_buffer == namedtuple_data

    # Necessary to remove any remnants of the memory mapping on windows, at least. Otherwise
    # we can find that the file remains open.
    GC.gc()
    GC.gc()
end

@testset "simple" begin
    some_data = [Tick(1000 + i, i, 2 * i) for i in 1:100]
    _test_api(some_data, nothing, nothing)
end

@testset "simple datetime" begin
    some_data = [
        DateTimeTick(DateTime(2000, 1, 1) + Day(i), 2 * i, 3 * i)
        for i in 1:100
    ]
    _test_api(some_data, nothing, nothing)
    _test_api(some_data, DateTime(1990, 1, 1), nothing)
    _test_api(some_data, DateTime(2000, 1, 1), nothing)
    _test_api(some_data, DateTime(2000, 3, 1), nothing)
    _test_api(some_data, DateTime(2000, 3, 1), DateTime(2100, 1, 1))
    _test_api(some_data, DateTime(2000, 3, 1), DateTime(2000, 4, 1))
end