text:
abstract type ApplyStyle end
abstract type AbstractArrayApplyStyle <: ApplyStyle end
struct DefaultApplyStyle <: ApplyStyle end
struct DefaultArrayApplyStyle <: AbstractArrayApplyStyle end

@inline ApplyStyle(f, args...) = DefaultApplyStyle()
@inline ApplyStyle(f, ::Type{<:AbstractArray}, args::Type{<:AbstractArray}...) = DefaultArrayApplyStyle()

struct Applied{Style, F, Args<:Tuple}
    f::F
    args::Args
end

@inline Applied{Style}(f::F, args::Args) where {Style,F,Args<:Tuple} = Applied{Style,F,Args}(f, args)
@inline Applied{Style}(A::Applied) where Style = Applied{Style}(A.f, A.args)

call(a) = a.f
call(_, a) = a.f
call(a::AbstractArray) = call(MemoryLayout(typeof(a)), a)
arguments(a) = a.args
arguments(_, a) = a.args
arguments(a::AbstractArray) = arguments(MemoryLayout(typeof(a)), a)

@inline check_applied_axes(A::Applied) = nothing

function instantiate(A::Applied{Style}) where Style
    check_applied_axes(A)
    Applied{Style}(A.f, map(instantiate, A.args))
end

@inline _typesof() = ()
@inline _typesof(a, b...) = tuple(typeof(a), _typesof(b...)...)
@inline _typesof(a, b) = tuple(typeof(a), typeof(b))
@inline _typesof(a, b, c) = tuple(typeof(a), typeof(b), typeof(c))
@inline combine_apply_style(f, a...) = ApplyStyle(f, _typesof(a...)...)
@inline combine_apply_style(f, a, b) = ApplyStyle(f, typeof(a), typeof(b))
@inline combine_apply_style(f, a, b, c) = ApplyStyle(f, typeof(a), typeof(b), typeof(c))

@inline Applied(f, args...) = Applied{typeof(combine_apply_style(f, args...))}(f, args)
@inline applied(f, args...) = Applied(f, args...)
@inline apply(f, args...) = materialize(applied(f, args...))
@inline apply!(f, args...) = materialize!(applied(f, args...))

materialize(A::Applied) = copy(instantiate(A))
materializeargs(A::Applied) = applied(A.f, materialize.(A.args)...)

# the following materialzes the args and calls materialize again, unless it hasn't
# changed in which case it falls back to the default
_default_materialize(A::App, ::App) where App = A.f(A.args...)
_default_materialize(A, _) = materialize(A)
# copy(A::Applied{DefaultApplyStyle}) = A.f(A.args...)
copy(A::Applied) = _default_materialize(materializeargs(A), A)

@inline copyto!(dest, M::Applied) = copyto!(dest, materialize(M))
@inline copyto!(dest::AbstractArray, M::Applied) = copyto!(dest, materialize(M))

broadcastable(M::Applied) = M

similar(M::Applied{<:AbstractArrayApplyStyle}, ::Type{T}, axes) where {T,N} = Array{T}(undef, length.(axes))
similar(M::Applied{<:AbstractArrayApplyStyle}, ::Type{T}) where T = similar(M, T, axes(M))
similar(M::Applied) = similar(M, eltype(M))

axes(A::Applied, j) = axes(A)[j]

struct ApplyArrayBroadcastStyle{N} <: Broadcast.AbstractArrayStyle{N} end
ApplyArrayBroadcastStyle{N}(::Val{N}) where N = ApplyArrayBroadcastStyle{N}()

BroadcastStyle(::Type{<:Applied}) = ApplyBroadcastStyle()
# temporary fix
BroadcastStyle(::Type{<:Applied{<:Any,typeof(*),<:Tuple{<:AbstractMatrix,<:AbstractVector}}}) = ApplyArrayBroadcastStyle{1}()
BroadcastStyle(::Type{<:Applied{<:Any,typeof(*),<:Tuple{<:AbstractMatrix,<:AbstractMatrix}}}) = ApplyArrayBroadcastStyle{2}()
BroadcastStyle(::Type{<:Applied{<:Any,typeof(*),<:Tuple{<:AbstractVector,<:AbstractMatrix}}}) = ApplyArrayBroadcastStyle{2}()

BroadcastStyle(::ApplyArrayBroadcastStyle{Any}, b::DefaultArrayStyle) = b
BroadcastStyle(::ApplyArrayBroadcastStyle{N}, b::DefaultArrayStyle{N}) where N = b
BroadcastStyle(::ApplyArrayBroadcastStyle{M}, b::DefaultArrayStyle{N}) where {M,N} =
    typeof(b)(Val(max(M, N)))

similar(bc::Broadcasted{ApplyArrayBroadcastStyle{N}}, ::Type{ElType}) where {N,ElType} =
    similar(Array{ElType}, axes(bc))

@inline function copyto!(dest::AbstractArray, bc::Broadcasted{ApplyArrayBroadcastStyle{N}}) where N
    if length(bc.args) == 1
        copyto!(dest, first(bc.args))
        if bc.f !== identity
            dest .= bc.f.(dest)
        end
    else
        bc′ = mapbc(bc) do x
            if x isa Applied
                materialize(x)
            else
                x
            end
        end
        materialize!(dest, bc′)
    end
    return dest
end

# Map over all nested Broadcasted and their arguments.  Using `broadcasted`
# instead of `Broadcasted` to re-process arguments via `broadcastable`.
@inline mapbc(f, bc::Broadcasted) =
    f(broadcasted(bc.f, map(a -> mapbc(f, a), bc.args)...))
@inline mapbc(f, x) = f(x)

struct MatrixFunctionStyle{F} <: AbstractArrayApplyStyle end

for f in (:exp, :sin, :cos, :sqrt)
    @eval ApplyStyle(::typeof($f), ::Type{<:AbstractMatrix}) = MatrixFunctionStyle{typeof($f)}()
end

function check_applied_axes(A::Applied{<:MatrixFunctionStyle})
    length(A.args) == 1 || throw(ArgumentError("MatrixFunctions only defined with 1 arg"))
    axes(A.args[1],1) == axes(A.args[1],2) || throw(DimensionMismatch("matrix is not square: dimensions are $axes(A.args[1])"))
end

for op in (:axes, :size)
    @eval begin
        $op(A::Applied{<:MatrixFunctionStyle}) = $op(first(A.args))
        $op(A::Applied{<:MatrixFunctionStyle}, j) = $op(first(A.args), j)
    end
end

ndims(A::Applied{<:MatrixFunctionStyle}) = ndims(first(A.args))

eltype(A::Applied{<:MatrixFunctionStyle}) = float(eltype(first(A.args)))

getindex(A::Applied, kj...) = materialize(A)[kj...]

"""
    LazyArray(x::Applied) :: ApplyArray
    LazyArray(x::Broadcasted) :: BroadcastArray

Wrap a lazy object that wraps a computation producing an array to an
array.
"""
abstract type LazyArray{T,N} <: AbstractArray{T,N} end

const LazyMatrix{T} = LazyArray{T,2}
const LazyVector{T} = LazyArray{T,1}

struct ApplyArray{T, N, F, Args<:Tuple} <: LazyArray{T,N}
    f::F
    args::Args
end

const ApplyVector{T, F, Args<:Tuple} = ApplyArray{T, 1, F, Args}
const ApplyMatrix{T, F, Args<:Tuple} = ApplyArray{T, 2, F, Args}

LazyArray(A::Applied) = ApplyArray(A)

ApplyArray{T,N,F,Args}(M::Applied) where {T,N,F,Args} = ApplyArray{T,N,F,Args}(M.f, M.args)
ApplyArray{T,N}(M::Applied{Style,F,Args}) where {T,N,Style,F,Args} = ApplyArray{T,N,F,Args}(instantiate(M))
ApplyArray{T}(M::Applied) where {T} = ApplyArray{T,ndims(M)}(M)
ApplyArray(M::Applied) = ApplyArray{eltype(M)}(M)
ApplyVector(M::Applied) = ApplyVector{eltype(M)}(M)
ApplyMatrix(M::Applied) = ApplyMatrix{eltype(M)}(M)

ApplyArray(f, factors...) = ApplyArray(applied(f, factors...))
ApplyArray{T}(f, factors...) where T = ApplyArray{T}(applied(f, factors...))
ApplyArray{T,N}(f, factors...) where {T,N} = ApplyArray{T,N}(applied(f, factors...))

ApplyVector(f, factors...) = ApplyVector(applied(f, factors...))
ApplyMatrix(f, factors...) = ApplyMatrix(applied(f, factors...))

convert(::Type{AbstractArray{T}}, A::ApplyArray{T}) where T = A
convert(::Type{AbstractArray{T}}, A::ApplyArray{<:Any,N}) where {T,N} = ApplyArray{T,N}(A.f, A.args...)
convert(::Type{AbstractArray{T,N}}, A::ApplyArray{T,N}) where {T,N} = A
convert(::Type{AbstractArray{T,N}}, A::ApplyArray{<:Any,N}) where {T,N} = ApplyArray{T,N}(A.f, A.args...)

AbstractArray{T}(A::ApplyArray{T}) where T = copy(A)
AbstractArray{T}(A::ApplyArray{<:Any,N}) where {T,N} = ApplyArray{T,N}(A.f, map(copy,A.args)...)
AbstractArray{T,N}(A::ApplyArray{T,N}) where {T,N} = copy(A)
AbstractArray{T,N}(A::ApplyArray{<:Any,N}) where {T,N} = ApplyArray{T,N}(A.f, map(copy,A.args)...)

@inline axes(A::ApplyArray) = axes(Applied(A))
@inline size(A::ApplyArray) = map(length, axes(A))
@inline copy(A::ApplyArray) = ApplyArray(A.f, map(copy,A.args)...)

struct LazyArrayApplyStyle <: AbstractArrayApplyStyle end
copy(A::Applied{LazyArrayApplyStyle}) = ApplyArray(A)

@propagate_inbounds getindex(A::ApplyArray{T,N}, kj::Vararg{Integer,N}) where {T,N} = convert(T, Applied(A)[kj...])::T
@propagate_inbounds getindex(A::Applied{LazyArrayApplyStyle}, kj...) = materialize(Applied{DefaultArrayApplyStyle}(A))[kj...]

for F in (:exp, :log, :sqrt, :cos, :sin, :tan, :csc, :sec, :cot,
            :cosh, :sinh, :tanh, :csch, :sech, :coth,
            :acosh, :asinh, :atanh, :acsch, :asech, :acoth,
            :acos, :asin, :atan, :acsc, :asec, :acot)
    @eval begin
        ndims(M::Applied{LazyArrayApplyStyle,typeof($F)}) = ndims(first(M.args))
        axes(M::Applied{LazyArrayApplyStyle,typeof($F)}) = axes(first(M.args))
        size(M::Applied{LazyArrayApplyStyle,typeof($F)}) = size(first(M.args))
        eltype(M::Applied{LazyArrayApplyStyle,typeof($F)}) = eltype(first(M.args))
    end
end

abstract type AbstractLazyLayout <: MemoryLayout end
struct LazyLayout <: AbstractLazyLayout end

MemoryLayout(::Type{<:LazyArray}) = LazyArrayLayout()

transposelayout(L::LazyLayout) = L
conjlayout(L::LazyLayout) = L
sublayout(L::LazyLayout, _) = L
reshapedlayout(::LazyLayout, _) = LazyLayout()

combine_mul_styles(::LazyLayout) = LazyArrayApplyStyle()
result_mul_style(::LazyArrayApplyStyle, ::LazyArrayApplyStyle) = LazyArrayApplyStyle()
result_mul_style(::LazyArrayApplyStyle, _) = LazyArrayApplyStyle()
result_mul_style(_, ::LazyArrayApplyStyle) = LazyArrayApplyStyle()

struct  ApplyLayout{F} <: MemoryLayout end

applylayout(::Type{F}, args...) where F = ApplyLayout{F}()

MemoryLayout(::Type{Applied{Style,F,Args}}) where {Style,F,Args} =
    applylayout(F, tuple_type_memorylayouts(Args)...)
MemoryLayout(::Type{ApplyArray{T,N,F,Args}}) where {T,N,F,Args} =
    applylayout(F, tuple_type_memorylayouts(Args)...)

@inline Applied(A::AbstractArray) = Applied(call(A), arguments(A)...)
@inline ApplyArray(A::AbstractArray) = ApplyArray(call(A), arguments(A)...)

function show(io::IO, A::Applied)
    print(io, "Applied(", A.f)
    for a in A.args
        print(io, ',', a)
    end
    print(io, ')')
end

applybroadcaststyle(::Type{<:AbstractArray{<:Any,N}}, _2) where N = DefaultArrayStyle{N}()
applybroadcaststyle(::Type{<:AbstractArray{<:Any,N}}, ::LazyLayout) where N = LazyArrayStyle{N}()
BroadcastStyle(M::Type{<:ApplyArray}) = applybroadcaststyle(M, MemoryLayout(M))

replace_in_print_matrix(A::LazyMatrix, i::Integer, j::Integer, s::AbstractString) =
    i in colsupport(A,j) ? s : replace_with_centered_mark(s)

###
# Number special cases
###

for op in (:+, :-, :*, :\)
    @eval applied(::typeof($op), x::Number, y::Number) = $op(x,y)
end

for op in (:one, :zero)
    @eval applied(::typeof($op), ::Type{T}) where T = $op(T)
end

###
# Lazy getindex
# this uses a lazy-materialize idiom to construct a matrix based
# on the memory layout
###

@inline getindex(A::LazyMatrix, kr::Colon, jr::Colon) = lazy_getindex(A, kr, jr)
@inline getindex(A::LazyMatrix, kr::Colon, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)
@inline getindex(A::LazyMatrix, kr::AbstractUnitRange, jr::Colon) = lazy_getindex(A, kr, jr)
@inline getindex(A::LazyMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = lazy_getindex(A, kr, jr)

@inline copyto!(dest::AbstractArray{T,N}, src::ApplyArray{T,N}) where {T,N} = copyto!(dest, Applied(src))
@inline copyto!(dest::AbstractArray, src::ApplyArray) = copyto!(dest, Applied(src))

# avoid infinite-loop
_base_copyto!(dest::AbstractArray{T,N}, src::AbstractArray{T,N}) where {T,N} = Base.invoke(copyto!, NTuple{2,AbstractArray{T,N}}, dest, src)
_base_copyto!(dest::AbstractArray, src::AbstractArray) = Base.invoke(copyto!, NTuple{2,AbstractArray}, dest, src)
@inline copyto!(dest::AbstractArray, M::Applied{LazyArrayApplyStyle}) = _base_copyto!(dest, materialize(M))

##
# triu/tril
##
for tri in (:tril, :triu)
    for op in (:axes, :size)
        @eval begin
            $op(A::Applied{<:Any,typeof($tri)}) = $op(first(A.args))
            $op(A::Applied{<:Any,typeof($tri)}, j) = $op(first(A.args), j)
        end
    end
    @eval begin
        ndims(::Applied{<:Any,typeof($tri)}) = 2
        eltype(A::Applied{<:Any,typeof($tri)}) = eltype(first(A.args))
    end
end

getindex(A::ApplyMatrix{T,typeof(triu),<:Tuple{<:AbstractMatrix}}, k::Integer, j::Integer) where T =
    j ≥ k ? A.args[1][k,j] : zero(T)

getindex(A::ApplyMatrix{T,typeof(triu),<:Tuple{<:AbstractMatrix,<:Integer}}, k::Integer, j::Integer) where T =
    j ≥ k+A.args[2] ? A.args[1][k,j] : zero(T)

replace_in_print_matrix(A::ApplyMatrix{<:Any,typeof(triu),<:Tuple{<:AbstractMatrix}}, i::Integer, j::Integer, s::AbstractString) =
    j ≥ i ? replace_in_print_matrix(A.args[1], i, j, s) : replace_with_centered_mark(s)
replace_in_print_matrix(A::ApplyMatrix{<:Any,typeof(triu),<:Tuple{<:AbstractMatrix,<:Integer}}, i::Integer, j::Integer, s::AbstractString) =
    j ≥ i+A.args[2] ? replace_in_print_matrix(A.args[1], i, j, s) : replace_with_centered_mark(s)


text:
using PyCall

# Functionality for computing band structures

function high_symmetry_kpath(model; kline_density=20)
    if model.n_dim == 1  # Return fast for 1D model
        # kline_density = Number of kpoint per Angström
        # Length of the kpath is lattice[1, 1] in 1D
        n_points = ceil(Int, kline_density / austrip(1u"Å") * model.lattice[1, 1])
        return (
            kcoords = [[coord, 0, 0] for coord in range(-1//2, 1//2, length=1+n_points)],
            klabels=Dict("Γ" => zeros(3), "-1/2" => [-0.5, 0.0, 0.0], "1/2" => [0.5, 0, 0]),
            kpath=[[raw"-1/2", raw"1/2"]],
        )
    end

    # TODO This is the last function that hard-depends on pymatgen. The way to solve this
    # is to use the julia version implemented in
    # https://github.com/louisponet/DFControl.jl/blob/master/src/structure.jl
    # but for this the best way to go would be to refactor into a small "CrystalStructure"
    # julia module which deals with these sort of low-level details everyone can agree on.
    pystructure = pymatgen_structure(model.lattice, model.atoms)
    symm_kpath = pyimport("pymatgen.symmetry.bandstructure").HighSymmKpath(pystructure)
    kcoords, labels = symm_kpath.get_kpoints(kline_density, coords_are_cartesian=false)

    labels_dict = Dict{String, Vector{eltype(kcoords[1])}}()
    for (ik, k) in enumerate(kcoords)
        if length(labels[ik]) > 0
            labels_dict[detexify_kpoint(labels[ik])] = k
        end
    end

    (kcoords=kcoords, klabels=labels_dict, kpath=symm_kpath.kpath["path"])
end

@timing function compute_bands(basis, ρ, kcoords, n_bands;
                               eigensolver=lobpcg_hyper,
                               tol=1e-3,
                               show_progress=true,
                               kwargs...)
    # Create basis with new kpoints, where we cheat by using any symmetry operations.
    ksymops = [[identity_symop()] for _ in 1:length(kcoords)]
    # For some reason rationalize(2//3) isn't supported (julia 1.6)
    myrationalize(x::T) where {T <: AbstractFloat} = rationalize(x, tol=10eps(T))
    myrationalize(x) = x
    bs_basis = PlaneWaveBasis(basis, [myrationalize.(k) for k in kcoords], ksymops)
    ham = Hamiltonian(bs_basis; ρ=ρ)

    band_data = diagonalize_all_kblocks(eigensolver, ham, n_bands + 3;
                                        n_conv_check=n_bands,
                                        tol=tol, show_progress=show_progress, kwargs...)
    if !band_data.converged
        @warn "Eigensolver not converged" iterations=band_data.iterations
    end
    merge((basis=bs_basis, ), select_eigenpairs_all_kblocks(band_data, 1:n_bands))
end

function split_into_branches(kcoords, data::Dict, klabels::Dict)
    # kcoords in cartesian coordinates, klabels uses cartesian coordinates
    function getlabel(kcoord; tol=1e-4)
        findfirst(c -> norm(c - kcoord) < tol, klabels)
    end

    branches = Any[(kindices = [0], kdistances=[0.0], ), ]
    for (ik, kcoord) in enumerate(kcoords)
        previous_kcoord = ik == 1 ? kcoords[1] : kcoords[ik - 1]
        if !isnothing(getlabel(kcoord)) && !isnothing(getlabel(previous_kcoord))
            # New branch encountered
            previous_distance = branches[end].kdistances[end]
            push!(branches, (kindices=[ik], kdistances=[previous_distance]))
        else
            # Keep adding to current branch
            distance = branches[end].kdistances[end] + norm(kcoord - kcoords[ik - 1])
            push!(branches[end].kdistances, distance)
            push!(branches[end].kindices, ik)
        end
    end

    map(branches[2:end]) do branch
        branch_data = Dict(key => data[key][branch.kindices, :, :] for key in keys(data))
        ret = (klabels=(getlabel(kcoords[branch.kindices[1]]),
                        getlabel(kcoords[branch.kindices[end]])),
               kdistances=branch.kdistances,
               kindices=branch.kindices)
        merge(ret, (; branch_data...))
    end
end

function prepare_band_data(band_data; datakeys=[:λ, :λerror],
                           klabels=Dict{String, Vector{Float64}}())
    basis = band_data.basis
    n_spin   = basis.model.n_spin_components
    n_kcoord = length(basis.kpoints) ÷ n_spin
    n_bands  = nothing

    # Convert coordinates to Cartesian
    kcoords_cart = [basis.kpoints[ik].coordinate_cart for ik in krange_spin(basis, 1)]
    klabels_cart = Dict(lal => basis.model.recip_lattice * vec for (lal, vec) in klabels)

    # Split data into branches
    data = Dict{Symbol, Any}()
    for key in datakeys
        hasproperty(band_data, key) || continue
        n_bands = length(band_data[key][1])
        data_per_kσ = similar(band_data[key][1], n_kcoord, n_bands, n_spin)
        for σ in 1:n_spin, (ito, ik) in enumerate(krange_spin(basis, σ))
            data_per_kσ[ito, :, σ] = band_data[key][ik]
        end
        data[key] = data_per_kσ
    end
    @assert !isnothing(n_bands)
    branches = split_into_branches(kcoords_cart, data, klabels_cart)

    tick_labels    = String[branches[1].klabels[1]]
    tick_distances = Float64[branches[1].kdistances[1]]
    for (i, br) in enumerate(branches)
        # Ignore branches with a single k-point
        branches[i].klabels[1] == branches[i].klabels[2] && continue

        label    = branches[i].klabels[2]
        distance = branches[i].kdistances[end]
        if i != length(branches) && branches[i+1].klabels[1] != label
            # Next branch is not continuous from the current
            label = label * " | " * branches[i+1].klabels[1]
        end
        push!(tick_labels, label)
        push!(tick_distances, distance)
    end

    (branches=branches, ticks=(distances=tick_distances, labels=tick_labels),
     n_bands=n_bands, n_kcoord=n_kcoord, n_spin=n_spin, basis=basis)
end

"""
    is_metal(band_data, εF, tol)

Determine whether the provided bands indicate the material is a metal,
i.e. where bands are cut by the Fermi level.
"""
function is_metal(band_data, εF, tol=1e-4)
    n_bands = length(band_data.λ[1])
    n_kpoints = length(band_data.λ)
    for ib in 1:n_bands
        some_larger = any(band_data.λ[ik][ib] - εF < -tol for ik in 1:n_kpoints)
        some_smaller = any(εF - band_data.λ[ik][ib] < -tol for ik in 1:n_kpoints)
        some_larger && some_smaller && return true
    end
    false
end

function detexify_kpoint(string)
    # For some reason Julia doesn't support this naively: https://github.com/JuliaLang/julia/issues/29849
    replacements = ("\\Gamma" => "Γ",
                    "\\Delta" => "Δ",
                    "\\Sigma" => "Σ",
                    "_1"      => "₁")
    for r in replacements
        string = replace(string, r)
    end
    string
end

"""
Compute and plot the band structure. `n_bands` selects the number of bands to compute.
If this value is absent and an `scfres` is used to start the calculation a default of
`n_bands_scf + 5sqrt(n_bands_scf)` is used. Unlike the rest of DFTK bands energies
are plotted in eV unless a different `unit` (any Unitful unit) is selected.
"""
function plot_bandstructure(basis, ρ, n_bands;
                            εF=nothing, kline_density=20, unit=u"eV", kwargs...)
    mpi_nprocs() > 1 && error("Band structures with MPI not supported yet")
    if !isdefined(DFTK, :PLOTS_LOADED)
        error("Plots not loaded. Run 'using Plots' before calling plot_bandstructure.")
    end

    # Band structure calculation along high-symmetry path
    kcoords, klabels, kpath = high_symmetry_kpath(basis.model; kline_density=kline_density)
    println("Computing bands along kpath:")
    println("       ", join(join.(detexify_kpoint.(kpath), " -> "), "  and  "))
    band_data = compute_bands(basis, ρ, kcoords, n_bands; kwargs...)

    plotargs = ()
    if kline_density ≤ 10
        plotargs = (markersize=2, markershape=:circle)
    end

    plot_band_data(band_data; εF=εF, klabels=klabels, unit=unit, plotargs...)
end
function plot_bandstructure(scfres; n_bands=nothing, kwargs...)
    # Convenience wrapper for scfres named tuples
    n_bands_scf = length(scfres.occupation[1])
    isnothing(n_bands) && (n_bands = ceil(Int, n_bands_scf + 5sqrt(n_bands_scf)))
    plot_bandstructure(scfres.basis, scfres.ρ, n_bands; εF=scfres.εF, kwargs...)
end

ext: jl


text:
mutable struct Flag{Bool} <: State{Bool}
    value::Bool
end

Flag(; _value, _...) = begin
    Flag{Bool}(_value)
end

constructortags(::Val{:Flag}) = ()

genvartype(v::VarInfo, ::Val{:Flag}; _...) = @q Flag{Bool}

gendefault(v::VarInfo, ::Val{:Flag}) = istag(v, :parameter) ? genparameter(v) : false

genupdate(v::VarInfo, ::Val{:Flag}, ::MainStep; kw...) = begin
    @gensym s f q
    if istag(v, :parameter)
        nothing
    elseif istag(v, :once)
        @q let $s = $(symstate(v))
            if !$C.value($s)
                let $f = $(genbody(v))
                    $C.store!($s, $f)
                end
            end
        end
    else
        @q let $s = $(symstate(v)),
               $f = $(genbody(v))
            $C.store!($s, $f)
        end
    end
end


text:
import PkgMirrors
import Pkg
import Test: @test

PkgMirrors.clear()

@test PkgMirrors.availables()[1] == "ZJU"

@test PkgMirrors.current() == nothing

PkgMirrors.setmirror("ZJU")

@test PkgMirrors.current() == PkgMirrors.PkgMirror("ZJU", "https://mirrors.zju.edu.cn/julia")

PkgMirrors.setmirror("USTC")

Pkg.update()

Pkg.add("RandomNumbers")

PkgMirrors.deactivate()

PkgMirrors.clear()


text:
#InspectDR: Defines stylesheets & presets
#-------------------------------------------------------------------------------

#==Important notes
***IEEE Plots***
  - Write out/export `Multiplot` object instead of `Plot` in order to control
    full plot dimensions instead of just dimensions of data area.
==#

#==Constants: Initial default values
===============================================================================#
#Default values for data area dimensions (used to save single plot):
const DEFAULT_DATA_WIDTH = 500.0
const DEFAULT_DATA_HEIGHT = DEFAULT_DATA_WIDTH / MathConstants.φ #Use golden ratio

#Default values for plot dimensions (used to save multi-plot):
const DEFAULT_PLOT_WIDTH = 600.0
const DEFAULT_PLOT_HEIGHT = DEFAULT_PLOT_WIDTH / MathConstants.φ #Use golden ratio

#==Types
===============================================================================#
mutable struct StyleID{T}; end
StyleID(T::Symbol) = StyleID{T}()

mutable struct PlotSylesheet
	plotlayout::PlotLayout
	mplotlayout::MultiplotLayout
end

#==Useful functions
===============================================================================#
#Try to size allocations & offsets appropriately for font selections:
function autofit2font!(lyt::PlotLayout)
	#Allocations for different elements:
	h_title = lyt.font_title._size
	h_axislabel = lyt.font_axislabel._size
	h_ticklabel = lyt.font_ticklabel._size
	w_axislabel = lyt.font_axislabel._size
	w_ticklabel = lyt.font_ticklabel._size

	lyt.valloc_top = h_title + h_ticklabel/2
	lyt.valloc_mid = h_ticklabel * 1.5
	lyt.valloc_bottom = h_axislabel + 1.75*h_ticklabel

	lyt.halloc_left = w_axislabel + 4.25*w_ticklabel
	lyt.halloc_right = w_ticklabel * 3.75 #Room for x10^EXP

	lyt.voffset_title = h_title*0.5
	lyt.voffset_xticklabel = h_ticklabel*0.5
	lyt.hoffset_yticklabel = w_ticklabel*0.25
	lyt.voffset_xaxislabel = h_axislabel*0.5
	lyt.hoffset_yaxislabel = w_axislabel*0.5

	lyt.halloc_legend = lyt.font_legend._size * (100/12)
	return lyt
end

#Try to size allocations & offsets appropriately for font selections:
function autofit2font!(lyt::MultiplotLayout)
	lyt.valloc_title = lyt.font_title._size * 1.5
	return lyt
end

function overwritefont!(lyt::PlotStyle; fontname=nothing, fontscale=1.0, autofit=true)
	proplist = [:font_title, :font_axislabel, :font_ticklabel, :font_annotation, :font_legend, :font_time]

	if fontname!=nothing
		for prop in proplist
			lyt[prop].name = fontname
		end
	end

	for prop in proplist
		ref = lyt[prop]
		ref._size *= fontscale
		lyt[prop] = ref #Updates layout.overwrite
	end

	if autofit
		autofit2font!(lyt.values)
	end
end

#==Preset Stylesheets
===============================================================================#

#Default ":screen" stylesheet:
#-------------------------------------------------------------------------------

#Default ":screen" PlotLayout stylesheet:
function getstyle(::Type{PlotLayout}, ::StyleID{:screen}, fontname::String, fontscale::Float64,
		notation_x::Symbol, notation_y::Symbol, enable_legend::Bool)
	lyt = PlotLayout()

	lyt.enable_legend = enable_legend
	lyt.enable_timestamp = false

	lyt.font_title = Font(DEFAULT_FONTNAME, 14*fontscale, bold=true)
	lyt.font_axislabel = Font(DEFAULT_FONTNAME, 14*fontscale)
	lyt.font_ticklabel = Font(DEFAULT_FONTNAME, 12*fontscale)
	lyt.font_annotation = Font(DEFAULT_FONTNAME, 12*fontscale)
	lyt.font_time = Font(DEFAULT_FONTNAME, 8*fontscale)
	lyt.font_legend = Font(DEFAULT_FONTNAME, 12*fontscale)

	lyt.valloc_data = DEFAULT_DATA_HEIGHT
	lyt.halloc_data = DEFAULT_DATA_WIDTH

	lyt.halloc_legendlineseg = 20
	lyt.hoffset_legendtext = 0.5
	lyt.valloc_legenditemsp = 0.25

	lyt.voffset_xticklabel = 7
	lyt.hoffset_yticklabel = 3

	function getticklabelstyle(notation::Symbol)
		if :ENG == notation
			TickLabelStyle()
		elseif :SI == notation
			TickLabelStyle(NumericIO.UEXPONENT_SI)
		else
			@warn("Unsupported tick label style: :$notation")
		end
	end
	lyt.format_xtick = getticklabelstyle(notation_x)
	lyt.format_ytick = getticklabelstyle(notation_y)

	lyt.frame_data = AreaAttributes(
		line=InspectDR.line(style=:solid, color=COLOR_BLACK, width=2)
	)

	return autofit2font!(lyt) #Compute offsets
end

#Default ":screen" MultiplotLayout stylesheet:
function getstyle(::Type{MultiplotLayout}, ::StyleID{:screen}, fontname::String, fontscale::Float64,
		notation_x::Symbol, notation_y::Symbol, enable_legend::Bool)
	lyt = MultiplotLayout()

	lyt.ncolumns = 1
	lyt.font_title = Font(DEFAULT_FONTNAME, 20*fontscale, bold=true)

	lyt.valloc_plot = DEFAULT_PLOT_HEIGHT
	lyt.halloc_plot = DEFAULT_PLOT_WIDTH

	return autofit2font!(lyt) #Compute offsets
end

#":screen" stylesheet: High-level interface:
getstyle(::Type{T}, ID::StyleID{:screen}; fontname::String=DEFAULT_FONTNAME, fontscale::Real=1.0,
		notation_x::Symbol=:ENG, notation_y::Symbol=:ENG, enable_legend::Bool=true) where T =
	getstyle(T, ID, fontname, Float64(fontscale), notation_x, notation_y, enable_legend)

#":IEEE" stylesheet:
#-------------------------------------------------------------------------------
#":IEEE" MultiplotLayout stylesheet:
function getstyle(::Type{PlotLayout}, ::StyleID{:IEEE}, ppi::Float64, fontscale::Float64, enable_legend::Bool)
	pt2px = ppi/DTPPOINTS_PER_INCH #WANTCONST
	lyt = PlotLayout()

	#IEEE Plot: Must ensure readable axis & tick labels

	fntaxis = Font(DEFAULT_FONTNAME, fontscale*7*pt2px)
	fnttick = fntaxis
	fntannot = Font(DEFAULT_FONTNAME, fontscale*5*pt2px)
	lyt.font_title = fntaxis
	lyt.font_axislabel = fntaxis
	lyt.font_ticklabel = fnttick
	lyt.font_annotation = fntannot
	lyt.font_time = fntannot
	lyt.font_legend = fntannot

	#Allocations for different elements:
	h_title = lyt.font_title._size
	h_axislabel = lyt.font_axislabel._size
	h_ticklabel = lyt.font_ticklabel._size
	w_axislabel = lyt.font_axislabel._size
	w_ticklabel = lyt.font_ticklabel._size

	#Optimize plot for maximum data area:
	lyt.valloc_top = h_ticklabel/2
	lyt.valloc_mid = h_ticklabel
	lyt.valloc_bottom = h_axislabel + 1.75*h_ticklabel
	lyt.halloc_left = w_axislabel + 2.5*h_ticklabel
	lyt.halloc_right = w_ticklabel * 3.75 #Room for x10^EXP

	lyt.voffset_title = 0
	lyt.voffset_xticklabel = h_ticklabel*0.5
	lyt.hoffset_yticklabel = w_ticklabel*0.25
	lyt.voffset_xaxislabel = h_axislabel*0.5
	lyt.hoffset_yaxislabel = w_axislabel*0.5

	lyt.enable_legend = enable_legend
	lyt.halloc_legend = 50*pt2px
	lyt.halloc_legendlineseg = 10*pt2px
	lyt.valloc_legenditemsp = 0.25
	lyt.hoffset_legendtext = 0.5

	lyt.frame_data.line = line(style=:solid, color=COLOR_BLACK, width=1*pt2px)
	return lyt
end

#":IEEE" MultiplotLayout stylesheet:
#(Write/export `Multiplot` instead of `Plot` to control full plot dimensions - not just data area)
function getstyle(::Type{MultiplotLayout}, ::StyleID{:IEEE}, ppi::Float64, fontscale::Float64, enable_legend::Bool)
	pt2px = ppi/DTPPOINTS_PER_INCH #WANTCONST
	wplot = 3.5*ppi #Inches => Pixels #WANTCONST

	lyt = MultiplotLayout()
	lyt.ncolumns = 1
	lyt.font_title = Font(DEFAULT_FONTNAME, 20*fontscale, bold=true)

	lyt.valloc_title = 0
	lyt.halloc_plot = wplot
	lyt.valloc_plot = wplot/MathConstants.φ #Golden ratio
	return lyt
end

#":IEEE" stylesheet: High-level interface:
getstyle(::Type{T}, ID::StyleID{:IEEE}; ppi::Real=300, fontscale::Real=1.0, enable_legend::Bool=true) where T =
	getstyle(T, ID, Float64(ppi), Float64(fontscale), enable_legend)

#==Preset Stylesheets: high-level interface
===============================================================================#
getstyle(::Type{T}, styleid::Symbol; kwargs...) where T =
	getstyle(T, StyleID(styleid); kwargs...)

function setstyle!(p::Plot, styleid::Symbol; refresh::Bool=true, kwargs...)
	s = getstyle(PlotLayout, StyleID(styleid); kwargs...)
	setstyle!(p.layout, s, refresh=refresh)
end
function setstyle!(p::Multiplot, styleid::Symbol; refresh::Bool=true, kwargs...)
	s = getstyle(MultiplotLayout, StyleID(styleid); kwargs...)
	setstyle!(p.layout, s, refresh=refresh)
end

#Last line


text:
# Autogenerated wrapper script for mpv_jll for armv7l-linux-musleabihf
export mpv

using Libiconv_jll
using SDL2_jll
using Zlib_jll
using FFMPEG_jll
using Lua_jll
using JpegTurbo_jll
using Xorg_libXrandr_jll
using Xorg_libXinerama_jll
using Libglvnd_jll
using Xorg_libX11_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `mpv`
const mpv_splitpath = ["bin", "mpv"]

# This will be filled out by __init__() for all products, as it must be done at runtime
mpv_path = ""

# mpv-specific global declaration
function mpv(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(mpv_path)
    end
end

# Inform that the wrapper is available for this platform
wrapper_available = true

"""
Open all libraries
"""
function __init__()
    # This either calls `@artifact_str()`, or returns a constant string if we're overridden.
    global artifact_dir = find_artifact_dir()

    global PATH_list, LIBPATH_list
    # Initialize PATH and LIBPATH environment variable listings
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (Libiconv_jll.PATH_list, SDL2_jll.PATH_list, Zlib_jll.PATH_list, FFMPEG_jll.PATH_list, Lua_jll.PATH_list, JpegTurbo_jll.PATH_list, Xorg_libXrandr_jll.PATH_list, Xorg_libXinerama_jll.PATH_list, Libglvnd_jll.PATH_list, Xorg_libX11_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (Libiconv_jll.LIBPATH_list, SDL2_jll.LIBPATH_list, Zlib_jll.LIBPATH_list, FFMPEG_jll.LIBPATH_list, Lua_jll.LIBPATH_list, JpegTurbo_jll.LIBPATH_list, Xorg_libXrandr_jll.LIBPATH_list, Xorg_libXinerama_jll.LIBPATH_list, Libglvnd_jll.LIBPATH_list, Xorg_libX11_jll.LIBPATH_list,))

    global mpv_path = normpath(joinpath(artifact_dir, mpv_splitpath...))

    push!(PATH_list, dirname(mpv_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()


text:
function _update_left_with_right!(x, y, ranges, allowmissing, mode::F; threads = true, op = nothing) where F
    if op === nothing
        @_threadsfor threads for i in 1:length(x)
            if length(ranges[i]) > 0
                if mode(x[i])
                    if !allowmissing && !ismissing(y[ranges[i].stop])
                        x[i] = y[ranges[i].stop]
                    elseif allowmissing
                        x[i] = y[ranges[i].stop]
                    end
                end
            end
        end
    else
        @_threadsfor threads for i in 1:length(x)
            if length(ranges[i]) > 0
                if mode(x[i])
                    if !allowmissing && !ismissing(y[ranges[i].stop])
                        x[i] = op(x[i], y[ranges[i].stop])
                    elseif allowmissing
                        x[i] = op(x[i], y[ranges[i].stop])
                    end
                end
            end
        end
    end
end

function _update!(dsl::Dataset, dsr::AbstractDataset, ::Val{T}; onleft, onright, check = true, allowmissing = true, mode = :all, mapformats = [true, true], stable = false, alg = HeapSort, accelerate = false, usehash = true, method = :sort, threads = true, op = nothing) where T
    isempty(dsl) && return dsl
    if method == :hash
        ranges, a, idx, minval, reps, sz, right_cols = _find_ranges_for_join_using_hash(dsl, dsr, onleft, onright, mapformats, true, Val(T); threads = threads)
    elseif method == :sort
        oncols_left = onleft
        oncols_right = onright
        right_cols = setdiff(1:length(index(dsr)), oncols_right)

        ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
        if usehash && length(oncols_left) == 1 && nrow(dsr)>1
            success, result = _update!_dict(dsl, dsr, ranges, oncols_left, oncols_right, right_cols, Val(T); mapformats = mapformats, allowmissing = allowmissing, mode = mode, threads = threads, op = op)
            if success
                return result
            end
        end
        idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate, threads = threads)

        for j in 1:length(oncols_left)
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j, threads = threads)
        end
    end

    if mode == :all
        f_mode = x->true
    elseif mode == :missing || mode == :missings
        f_mode = x->ismissing(x)
    else
        f_mode = x->mode(x)::Bool
    end
    for j in 1:length(right_cols)
        if haskey(index(dsl).lookup, _names(dsr)[right_cols[j]])
            left_cols_idx = index(dsl)[_names(dsr)[right_cols[j]]]
            TL = nonmissingtype(eltype(_columns(dsl)[left_cols_idx]))
            TR = nonmissingtype(eltype(_columns(dsr)[right_cols[j]]))
            if promote_type(TR, TL) <: TL
                _update_left_with_right!(_columns(dsl)[left_cols_idx], view(_columns(dsr)[right_cols[j]], idx), ranges, allowmissing, f_mode, threads = threads, op = op)
            end
        end
    end
    _modified(_attributes(dsl))
    dsl
end


text:
using KitBase, Plots, OffsetArrays
using ProgressMeter: @showprogress

begin
    set = Setup(
        "gas",
        "cylinder",
        "2d0f",
        "hll",
        "nothing",
        1, # species
        3, # order of accuracy
        "vanleer", # limiter
        "extra",
        0.1, # cfl
        1.0, # time
    )
    ps = PSpace2D(0.0, 2.0, 100, 0.0, 1.0, 50, 1, 1)
    vs = VSpace2D(-2.0, 2.0, 5, -2.0, 2.0, 5)
    gas = Gas(
        1e-6,
        1.1742169734124273, # Mach
        1.0,
        3.0, # K
        7/5,
        0.81,
        1.0,
        0.5,
    )
    ib = nothing
    ks = SolverSet(set, ps, vs, gas, ib)
end

ctr = OffsetArray{ControlVolume2D}(undef, axes(ks.ps.x, 1), axes(ks.ps.y, 2))
a1face = Array{Interface2D}(undef, ks.ps.nx + 1, ks.ps.ny)
a2face = Array{Interface2D}(undef, ks.ps.nx, ks.ps.ny + 1)

for j in axes(ctr, 2), i in axes(ctr, 1)
    t1 = ib_rh(ks.gas.Ma, ks.gas.γ, rand(3))[2]
    t2 = ib_rh(ks.gas.Ma, ks.gas.γ, rand(3))[6]

    if ps.x[i, j] <= ps.x1*0.25
        prim = [t2[1], t1[2] - t2[2], 0.0, t2[3]]
    else
        prim = [t1[1], 0.0, 0.0, t1[3]]

        s = prim[1]^(1-ks.gas.γ) / (2 * prim[end])

        κ = 0.25 # vortex strength
        μ = 0.204
        rc = 0.05

        x0 = 0.8
        y0 = 0.5

        r = sqrt((ps.x[i, j] - x0)^2 + (ps.y[i, j] - y0)^2)

        η = r / rc

        δu = κ * η * exp(μ * (1-η^2)) * (ps.y[i, j] - y0) / r
        δv = -κ * η * exp(μ * (1-η^2)) * (ps.x[i, j] - x0) / r
        δT = -(ks.gas.γ-1)*κ^2/(8*μ*ks.gas.γ)*exp(2*μ*(1-η^2))

        T0 = 1 / prim[end]

        ρ = prim[1]^(ks.gas.γ-1) * (T0+δT) / T0^(1/(ks.gas.γ-1))
        #@show prim[1]

        #@show (1 / (s * 2 * (prim[end] + δλ)))

        #ρ = (1 / (s * 2 * (prim[end] + δλ)))^(1/(ks.gas.γ-1))
        prim1 = [ρ, prim[2]+δu, prim[3]+δv, 1/(1/prim[4]+δT)]

        if r <= rc * 8
            prim .= prim1
        end
    end
    w = prim_conserve(prim, ks.gas.γ)

    ctr[i, j] = ControlVolume2D(
        ks.ps.x[i, j],
        ks.ps.y[i, j],
        ks.ps.dx[i, j],
        ks.ps.dy[i, j],
        w,
        prim,
    )
end

for j = 1:ks.ps.ny
    for i = 1:ks.ps.nx+1
        a1face[i, j] = Interface2D(ks.ps.dy[i, j], 1.0, 0.0, zeros(4))
    end
end
for i = 1:ks.ps.nx
    for j = 1:ks.ps.ny+1
        a2face[i, j] =
            Interface2D(ks.ps.dx[i, j], 0.0, 1.0, zeros(4))
    end
end

iter = 0
t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(floor(ks.set.maxTime / dt)) + 1
res = zero(4)

@showprogress for iter = 1:nt
    reconstruct!(ks, ctr)
    evolve!(ks, ctr, a1face, a2face, dt)
    #update!(ks, ctr, a1face, a2face, dt, res)

    @inbounds for j = 1:ks.ps.ny
        for i = 1:ks.ps.nx
            KitBase.step!(
                ctr[i, j].w,
                ctr[i, j].prim,
                a1face[i, j].fw,
                a1face[i+1, j].fw,
                a2face[i, j].fw,
                a2face[i, j+1].fw,
                ks.gas.γ,
                ctr[i, j].dx * ctr[i, j].dy,
                zeros(4),
                zeros(4),
            )
        end
    end
    for i = 1:ks.ps.nx
        ctr[i, 0].w .= ctr[i, 1].w
        ctr[i, 0].prim .= ctr[i, 1].prim
        ctr[i, ks.ps.ny+1].w .= ctr[i, ks.ps.ny].w
        ctr[i, ks.ps.ny+1].prim .= ctr[i, ks.ps.ny].prim
    end
    for j = 1:ks.ps.ny
        ctr[ks.ps.nx+1, j].w .= ctr[ks.ps.nx, j].w
        ctr[ks.ps.nx+1, j].prim .= ctr[ks.ps.nx, j].prim
    end
end

begin
    sol = zeros(ps.nx, ps.ny, 4)
    for i = 1:ps.nx, j = 1:ps.ny
        sol[i, j, :] .= conserve_prim(ctr[i, j].w, ks.gas.γ)
        sol[i, j, 4] = 1 / sol[i, j, 4]
    end
    x = ps.x[1:ps.nx, 1]
    plot(x, 0.5 .* (sol[:, end÷2, 1] + sol[:, end÷2+1, 1]))
end

plot_contour(ks, ctr)

using JLD2
cd(@__DIR__)
@save "fvm9.jld2" x sol ctr

max_stars_repo_name: CSMMLab/FRSG

max_issues_repo_name: CSMMLab/FRSG

max_forks_repo_name: CSMMLab/FRSG


text:
using MTKHelpers
using CairoMakie # to load cairomakie functions
using Documenter

DocMeta.setdocmeta!(MTKHelpers, :DocTestSetup, :(using MTKHelpers); recursive=true)

makedocs(;
    #modules=[MTKHelpers], # do not check for not included docstrings
    authors="Thomas Wutzler <twutz@bgc-jena.mpg.de> and contributors",
    repo="https://github.com/bgctw/MTKHelpers.jl/blob/{commit}{path}#{line}",
    sitename="MTKHelpers.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/MTKHelpers.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Update parameters" => "problemparsetter.md",
        "Solution handling" => "solution.md",
        "Embedding a system" => "embed_system.md",
        "Smooth steps" => "smoothstep.md",
        "CairoMakie Helpers" => "cairomakie.md",
        "index" => "zindex.md",
    ],
)

deploydocs(;
    repo="github.com/bgctw/MTKHelpers.jl",
    devbranch="main",
)

max_stars_repo_name: bgctw/MTKHelpers

max_issues_repo_name: bgctw/MTKHelpers

max_forks_repo_name: bgctw/MTKHelpers