text:
function get_af0(obj::GNSSOrbitalElements)
    return jcall(obj, "getAf0", jdouble, ())
end

function get_af1(obj::GNSSOrbitalElements)
    return jcall(obj, "getAf1", jdouble, ())
end

function get_af2(obj::GNSSOrbitalElements)
    return jcall(obj, "getAf2", jdouble, ())
end

function get_aodc(obj::BeidouOrbitalElements)
    return jcall(obj, "getAODC", jint, ())
end

function get_aode(obj::BeidouOrbitalElements)
    return jcall(obj, "getAODE", jint, ())
end

function get_cic(obj::GNSSOrbitalElements)
    return jcall(obj, "getCic", jdouble, ())
end

function get_cis(obj::GNSSOrbitalElements)
    return jcall(obj, "getCis", jdouble, ())
end

function get_crc(obj::GNSSOrbitalElements)
    return jcall(obj, "getCrc", jdouble, ())
end

function get_crs(obj::GNSSOrbitalElements)
    return jcall(obj, "getCrs", jdouble, ())
end

function get_cuc(obj::GNSSOrbitalElements)
    return jcall(obj, "getCuc", jdouble, ())
end

function get_cus(obj::GNSSOrbitalElements)
    return jcall(obj, "getCus", jdouble, ())
end

function get_date(obj::TimeStamped)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::GNSSOrbitalElements)
    return jcall(obj, "getE", jdouble, ())
end

function get_i0(obj::GNSSOrbitalElements)
    return jcall(obj, "getI0", jdouble, ())
end

function get_i_dot(obj::GNSSOrbitalElements)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_iod(obj::BeidouOrbitalElements)
    return jcall(obj, "getIOD", jint, ())
end

function get_m0(obj::GNSSOrbitalElements)
    return jcall(obj, "getM0", jdouble, ())
end

function get_mean_motion(obj::GNSSOrbitalElements)
    return jcall(obj, "getMeanMotion", jdouble, ())
end

function get_omega0(obj::GNSSOrbitalElements)
    return jcall(obj, "getOmega0", jdouble, ())
end

function get_omega_dot(obj::GNSSOrbitalElements)
    return jcall(obj, "getOmegaDot", jdouble, ())
end

function get_pa(obj::GNSSOrbitalElements)
    return jcall(obj, "getPa", jdouble, ())
end

function get_prn(obj::GNSSOrbitalElements)
    return jcall(obj, "getPRN", jint, ())
end

function get_sma(obj::GNSSOrbitalElements)
    return jcall(obj, "getSma", jdouble, ())
end

function get_tgd1(obj::BeidouOrbitalElements)
    return jcall(obj, "getTGD1", jdouble, ())
end

function get_tgd2(obj::BeidouOrbitalElements)
    return jcall(obj, "getTGD2", jdouble, ())
end

function get_time(obj::GNSSOrbitalElements)
    return jcall(obj, "getTime", jdouble, ())
end

function get_toc(obj::GNSSOrbitalElements)
    return jcall(obj, "getToc", jdouble, ())
end

function get_week(obj::GNSSOrbitalElements)
    return jcall(obj, "getWeek", jint, ())
end


text:
function get_svn(url::String, dest::String)

  if isdir(dest)
    println(dest * " exists; not downloading.")
  else
    println("dowloading via SVN...")
    status = 1
    if Sys.iswindows()
      status = (try
        p = run(`cmd /c svn export $url $dest`)
        0
      catch err
        @warn(string("error thrown: ", err))
        1
      end)
    else
      status = (try
        p = run(`svn export $url $dest`)
        p.exitcode
      catch err
        @warn(string("error thrown: ", err))
        1
      end)
    end
    if status != 0
      err_string = "download failed. Is a command-line SVN client installed?

      (type \"run(`svn --version`)\"; if an error occurs, SVN isn't installed.)

      Subversion for Ubuntu: sudo apt install subversion
      Subversion for OS X: pkg_add subversion
      SlikSVN Windows client: https://sliksvn.com/download/
      "
      error(err_string)
    end
  end
  return nothing
end

ext: jl


text:
### ============================================================================
### Method related to lasso path associated with type `LassoRoot`
### ============================================================================

"""
    Get the lassopath for the instance made, the function will modify the
    λs field, and the LassoPath field of the instance.

"""
function GetLassoPath(this::LassoRoot, decay_by=0.9)

    @assert decay_by != 1 "lambda decay value cannot be 1. "
    if decay_by < 1
        decay_by = 1/decay_by
    end

    u = this.u
    A = this.Z
    y = this.l
    Results = Vector{Vector}()

    function λMax(A, y)
        ToMax = A'*(y .- u)
        ToMax *= 2
        return maximum(abs.(ToMax))
    end

    λ = λMax(A, y)
    λs = Vector{Float64}()
    Changeλ!(this, λ)
    x = SolveLasso(this)
    dx = Inf
    push!(Results, x)
    push!(λs, λ)
    pb = Pm.ProgressThresh(this.Tol, "inf norm of δx: ")
    while dx >= this.Tol && λ >= this.λMin
        λ /= decay_by
        push!(λs, λ)
        Changeλ!(this, λ)
        x = SolveLasso(this, x0=x)
        push!(Results, value.(x))
        dx = norm(Results[end - 1] - Results[end], 1)/norm(Results[end - 1], 1)
        Pm.update!(pb, dx)
    end

    ResultsMatrix = zeros(length(x), length(Results))
    for II ∈ 1:length(Results)
        ResultsMatrix[:, II] = Results[II]
    end

    # Store it.
    this.LassoPath = ResultsMatrix
    this.λs = λs
    return ResultsMatrix, λs
end

function VisualizeLassoPath(
        this::LassoRoot;
        fname::Union{String, Nothing}=nothing,
        title::Union{String, Nothing}=nothing
    )
    """
        Make a plots for the lasso path and save it in the pwd.

    """
    @assert isdefined(this, :LassoPath) "Lasso Path not defined for the object"*
    "yet, use GetLassoPath() function on it first".
    error("Haven't implemented it yet.")
    λs = this.λs
    Paths = this.LassoPath
    Loggedλ = log2.(λs)
    p1 = Plots.heatmap(
            Paths[end:-1:begin, :],
            title=title===nothing ? "Lasso Path" : title
        )
    p2 = Plots.plot(Loggedλ, Paths', label=nothing)
    Plots.xlabel!(p2, "log_2(λ)")
    p = Plots.plot(p1, p2, layout=(2, 1))
    Plots.plot!(size=(600, 800))
    Plots.plot!(dpi=400)
    Plots.savefig(p, fname===nothing ? "plot.png" : fname)
    return
end

"""
    Caputre the indices for the most important predictors from the
    regression. Returns the indices of the important weights.

    It decreases the λ geometrically and count the number of non-zero
    weigthts at for each λ. Once the number of non-zero weights reached
    the threshold (top_k), it will return the indices of these non-zero
    weights.

    ---
    this::LassoRoot:
        An instance of the Lasso Analyzer.

    top_k::Union{Float64, Int64}:
        The number of non-zero weights you want, if it's a float between 0, 1, then
        it's interpreted as the ratio of non-zero weights you want.

    threshold::Float64=1e-4
        Set a threshold for the weights that you deem to be "important", only
        weights larger than this threshold will be viiwed as "non-zero".

    ---
    returns:
        * indices of the important predictors
        * values of the regularization parameters
        * the actual weigths
"""
function CaptureImportantWeights(
        this::LassoRoot,
        top_k::Union{Float64, Int64} = 0.5,
        threshold::Float64=1e-10
    )
    @assert isdefined(this, :LassoPath) "Lasso Path not defined for this"*
    "Object yet. "
    @assert top_k >= 0 "this parameters, should be a positive number"
    @assert threshold >= 0 "This parameters shouold be a positive number"

    Paths = this.LassoPath
    top_k::Int64 = top_k <= 1 ? ceil(size(Paths, 1)*0.5) : round(top_k)
    for JJ in 1:size(Paths, 2)
        Col = view(Paths, :, JJ)
        BigEnough = sum(abs.(Col) .>= threshold)
        # rank them by abs and returns the indices for top k weights
        if BigEnough >= top_k
            Indices = sortperm(abs.(Col), rev=true)
            Indices = Indices[begin:top_k]
            return Indices, this.λs[JJ], this.LassoPath[Indices, JJ]
        end
    end
    return Vector{Int64}()
end

"""
    Find the weights that always monotinically increases as the value of
    λ decreases.

    Monotone weights are likely to be independent and have less covariances
    with other predictors.
"""
function CaptureMonotoneWeights(this::LassoRoot)
    # TODO: Implement his

end

function WeightCorrelations(this::LassoRoot, index::Int64)
    # TODO: implement this
end


text:
# Average fields along the field directions (i.e., Fw along the w-direction).  See the
# description of matrix/mean.jl/create_m̂().  Technically, apply_m̂!() can be used to average
# fields normal to the direction of averaging, but it is not used that way in apply_mean!().

# Assumes the space dimension and field dimension are the same.  In other words, when the
# space coordinate indices are (i,j,k), then the field has three vector components.
# Therefore, for the input field array F[i,j,k,w], we assume w = 1:3.

# The functions calculate the average and add to the output array, instead of replacing the
# values stored in the output array.  Therefore, if the derivative values themselves are
# desired, pass the output array initialized with zeros.

# Inside @threads, avoid defining β, which requires another let block to avoid unnecessary
# allocations.

export apply_m̂!, apply_mean!

# Wrapper for converting AbstractVector's to SVec's for arithmetic averaging
# This corresponds to the wrappers to create discrete versions of differential operators,
# but the nonzero enries in the matrix are 0.5 rather than 1.0.
apply_mean!(G::AbsArrNumber,  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
            F::AbsArrNumber,  # input field; G[i,j,k,w] is w-component of G at (i,j,k)
            ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
            isfwd::AbsVecBool,  # isfwd[w] = true|false: m̂_w is forward|backward averaging
            isbloch::AbsVecBool=fill(true,length(isfwd)),  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(length(isfwd));  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
            α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α mean(F)
            ) where {OP} =
    (K = length(isfwd); apply_mean!(G, F, Val(OP), SBool{K}(isfwd), SBool{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ); α))

# Wrapper for converting AbstractVector's to SVec's for weighted averaging
apply_mean!(G::AbsArrNumber,  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
            F::AbsArrNumber,  # input field; G[i,j,k,w] is w-component of G at (i,j,k)
            ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
            isfwd::AbsVecBool,  # isfwd[w] = true|false: m̂_w is forward|backward averaging
            ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with; vectors of length N
            ∆l′⁻¹::NTuple{K,AbsVecNumber},  # inverse of line segments to divide by; vectors of length N
            isbloch::AbsVecBool=fill(true,K),  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
            e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
            α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α mean(F)
            ) where {K,OP} =
    apply_mean!(G, F, Val(OP), SBool{K}(isfwd), ∆l, ∆l′⁻¹, SBool{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ); α)

# Concrete implementation for arithmetic averaging
# Initially this was implemented by passing a tuple of ones() os ∆l and ∆l′⁻¹, but because
# arithmetic averaging is used frequently, the optimized version is implemented separately.
function apply_mean!(G::AbsArrNumber{K₊₁},  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
                     F::AbsArrNumber{K₊₁},  # input field; G[i,j,k,w] is w-component of G at (i,j,k)
                     ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                     isfwd::SBool{K},  # isfwd[w] = true|false: m̂_w is forward|backward averaging
                     isbloch::SBool{K},  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
                     e⁻ⁱᵏᴸ::SNumber{K};  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
                     α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α mean(F)
                     ) where {K,K₊₁,OP}  #  K is space dimension; K₊₁ = K + 1
    @assert K₊₁==K+1
    for nw = 1:K  # direction of averaging
        nv = nw  # component of output field generated by w-directional averaging operator m̂_w (row index of matrix block)
        Gv = selectdim(G, K₊₁, nv)  # nv-th component of output field

        nu = nw  # component of input field to feed to w-directional averaging operator m̂_w (column index of matrix block)
        Fu = selectdim(F, K₊₁, nu)  # nu-th component of input field

        apply_m̂!(Gv, Fu, Val(OP), nw, isfwd[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]; α)
    end
end

# Concrete implementation for weighted averaging
# For the implementation details, see the comments in matrix/mean.jl.
function apply_mean!(G::AbsArrNumber{K₊₁},  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
                     F::AbsArrNumber{K₊₁},  # input field; G[i,j,k,w] is w-component of G at (i,j,k)
                     ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                     isfwd::SBool{K},  # isfwd[w] = true|false: m̂_w is forward|backward averaging
                     ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with; vectors of length N
                     ∆l′⁻¹::NTuple{K,AbsVecNumber},  # inverse of line segments to divide by; vectors of length N
                     isbloch::SBool{K},  # isbloch[w] = true|false : w-boundary condition is Bloch|symmetric
                     e⁻ⁱᵏᴸ::SNumber{K};  # e⁻ⁱᵏᴸ[w]: Bloch phase factor in w-direction
                     α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α mean(F)
                     ) where {K,K₊₁,OP}  #  K is space dimension; K₊₁ = K + 1
    @assert K₊₁==K+1
    for nw = 1:K  # direction of averaging
        nv = nw  # component of output field generated by w-directional averaging operator m̂_w (row index of matrix block)
        Gv = selectdim(G, K₊₁, nv)  # nv-th component of output field

        nu = nw  # component of input field to feed to w-directional averaging operator m̂_w (column index of matrix block)
        Fu = selectdim(F, K₊₁, nu)  # nu-th component of input field

        apply_m̂!(Gv, Fu, Val(OP), nw, isfwd[nw], ∆l[nw], ∆l′⁻¹[nw], isbloch[nw], e⁻ⁱᵏᴸ[nw]; α)
    end
end

## Field-averaging operators "m̂" (as used in de Moerloose and de Zutter)
#
# This applies the averaging operator for a single Cartesian component.  For the operator
# for all three Cartesian components, use apply_mean!.

# Below, I set the default value for e⁻ⁱᵏᴸ below, even though I don't usually set default
# values for the arguments in concrete implementation (except for keyword arguments, which
# require default values).  This is to make the arguments of apply_m̂!() to have default
# values when the corresponding arguments of create_m̂() do.

# Arithmetic m̂ for 3D
function apply_m̂!(Gv::AbsArrNumber{3},  # v-component of output field (v = x, y, z)
                  Fu::AbsArrNumber{3},  # u-component of input field (u = x, y, z)
                  ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                  nw::Integer,  # 1|2|3 for averaging along x|y|z
                  isfwd::Bool,  # true|false for forward|backward averaging
                  isbloch::Bool,  # boundary condition in w-direction
                  e⁻ⁱᵏᴸ::Number=1.0;  # Bloch phase factor
                  α::Number=1.0  # scale factor to multiply to result before adding it to Gv: Gv += α m(Fu)
                  ) where {OP}
    @assert size(Gv)==size(Fu)

    Nx, Ny, Nz = size(Fu)
    α2 = 0.5α

    # Make sure not to include branches inside for loops.
    if isfwd
        if nw == 1  # w = x
            if isbloch
                # 1. At locations except for the positive end of the x-direction
                @threads for k = 1:Nz
                    for j = 1:Ny, i = 1:Nx-1
                        @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i+1,j,k] + Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the positive end of the x-direction (where the boundary fields are
                # taken from the negative-end boundary)
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j,k), α2 * (e⁻ⁱᵏᴸ*Fu[1,j,k] + Fu[Nx,j,k]), Val(OP))  # Fu[Nx+1,j,k] = exp(-i kx Lx) * Fu[1,j,k]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # x-direction
                @threads for k = 1:Nz
                    for j = 1:Ny, i = 2:Nx-1
                        @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i+1,j,k] + Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the negative end of the x-direction (where the boundary fields are
                # assumed zero)
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j,k), α2 * Fu[2,j,k], Val(OP))  # Fu[1,j,k] == 0
                end

                # 3. At the positive end of the x-direction (where the boundary fields are
                # assumed zero)
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j,k), α2 * Fu[Nx,j,k], Val(OP))  # Fu[Nx+1,j,k] == 0
                end
            end
        elseif nw == 2  # w = y
            if isbloch
                # 1. At locations except for the positive end of the y-direction
                @threads for k = 1:Nz
                    for j = 1:Ny-1, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j+1,k] + Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the positive end of the y-direction (where the boundary fields are
                # taken from the negative-end boundary)
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny,k), α2 * (e⁻ⁱᵏᴸ*Fu[i,1,k] + Fu[i,Ny,k]), Val(OP))  # Fu[i,Ny+1,k] = exp(-i ky Ly) * Fu[i,1,k]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # y-direction
                @threads for k = 1:Nz
                    for j = 2:Ny-1, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j+1,k] + Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the negative end of the y-direction (where the boundary fields are
                # assumed zero)
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1,k), α2 * Fu[i,2,k], Val(OP))  # Fu[i,1,k] == 0
                end

                # 3. At the positive end of the y-direction (where the boundary fields are
                # assumed zero)
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny,k), α2 * Fu[i,Ny,k], Val(OP))  # Fu[i,Ny+1,k] == 0
                end
            end
        else  # nw == 3; w = z
            if isbloch
                # 1. At locations except for the positive end of the z-direction
                @threads for k = 1:Nz-1
                    for j = 1:Ny, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j,k+1] + Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the positive end of the z-direction (where the boundary fields are
                # taken from the negative-end boundary)
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,Nz), α2 * (e⁻ⁱᵏᴸ*Fu[i,j,1] + Fu[i,j,Nz]), Val(OP))  # Fu[i,j,Nz+1] = exp(-i kz Lz) * Fu[i,j,1]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # z-direction
                @threads for k = 2:Nz-1
                    for j = 1:Ny, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j,k+1] + Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the negative end of the z-direction (where the boundary fields are
                # assumed zero)
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,1), α2 * Fu[i,j,2], Val(OP))  # Fu[i,j,1] == 0
                end

                # 3. At the positive end of the z-direction (where the boundary fields are
                # assumed zero)
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,Nz), α2 * Fu[i,j,Nz], Val(OP))  # Fu[i,j,Nz+1] == 0
                end
            end
        end  # if nw == ...
    else  # backward averaging
        if nw == 1  # w = x
            # 1. At the locations except for the negative end of the x-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for k = 1:Nz
                for j = 1:Ny, i = 2:Nx  # not i = 2:Nx-1
                    @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j,k] + Fu[i-1,j,k]), Val(OP))
                end
            end

            # 2. At the negative end of the x-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j,k), α2 * (Fu[1,j,k] + Fu[Nx,j,k]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[0,j,k] = Fu[Nx,j,k] / exp(-i kx Lx)
                end
            else  # symmetry boundary
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j,k), α * Fu[1,j,k], Val(OP))  # Fu[0,j,k] = Fu[1,j,k]; use α instead of α2 as field is doubled
                end
            end
        elseif nw == 2  # w = y
            # 1. At the locations except for the negative end of the y-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for k = 1:Nz
                for j = 2:Ny, i = 1:Nx  # not j = 2:Ny-1
                    @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j,k] + Fu[i,j-1,k]), Val(OP))
                end
            end

            # 2. At the negative end of the y-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1,k), α2 * (Fu[i,1,k] + Fu[i,Ny,k]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[i,0,k] = Fu[0,Ny,k] / exp(-i ky Ly)
                end
            else  # symmetry boundary
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1,k), α * Fu[i,1,k], Val(OP))  # Fu[i,0,k] = Fu[i,1,k]; use α instead of α2 as field is doubled
                end
            end
        else  # nw == 3; w = z
            # 1. At the locations except for the negative end of the z-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for k = 2:Nz
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,k), α2 * (Fu[i,j,k] + Fu[i,j,k-1]), Val(OP))
                end
            end

            # 2. At the negative end of the z-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,1), α2 * (Fu[i,j,1] + Fu[i,j,Nz]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[i,j,0] = Fu[i,j,Nz] / exp(-i kz Lz)
                end
            else  # symmetry boundary
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,1), α * Fu[i,j,1], Val(OP))  # Fu[i,j,0] = Fu[i,j,Nz]; use α instead of α2 as field is doubled
                end
            end
        end  # if nw == ...
    end  # if isfwd

    return nothing
end

# Arithmetic m̂ for 2D
function apply_m̂!(Gv::AbsArrNumber{2},  # v-component of output field (v = x, y)
                  Fu::AbsArrNumber{2},  # u-component of input field (u = x, y)
                  ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                  nw::Integer,  # 1|2 for averaging along x|y
                  isfwd::Bool,  # true|false for forward|backward averaging
                  isbloch::Bool,  # boundary condition in w-direction
                  e⁻ⁱᵏᴸ::Number=1.0;  # Bloch phase factor
                  α::Number=1.0  # scale factor to multiply to result before adding it to Gv: Gv += α m(Fu)
                  ) where {OP}
    @assert size(Gv)==size(Fu)

    Nx, Ny = size(Fu)
    α2 = 0.5 * α

    # Make sure not to include branches inside for loops.
    if isfwd
        if nw == 1  # w = x
            if isbloch
                # 1. At locations except for the positive end of the x-direction
                @threads for j = 1:Ny
                    for i = 1:Nx-1
                        @inbounds set_or_add!(Gv, (i,j), α2 * (Fu[i+1,j] + Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the positive end of the x-direction (where the boundary fields are
                # taken from the negative-end boundary)
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j), α2 * (e⁻ⁱᵏᴸ*Fu[1,j] + Fu[Nx,j]), Val(OP))  # Fu[Nx+1,j] = exp(-i kx Lx) * Fu[1,j]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # x-direction
                @threads for j = 1:Ny
                    for i = 2:Nx-1
                        @inbounds set_or_add!(Gv, (i,j), α2 * (Fu[i+1,j] + Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the negative end of the x-direction (where the boundary fields are
                # assumed zero)
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j), α2 * Fu[2,j], Val(OP))  # Fu[1,j] == 0
                end

                # 3. At the positive end of the x-direction (where the boundary fields are
                # assumed zero)
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j), α2 * Fu[Nx,j], Val(OP))  # Fu[Nx+1,j] == 0
                end
            end
        else  # nw == 2; w = y
            if isbloch
                # 1. At locations except for the positive end of the y-direction
                @threads for j = 1:Ny-1
                    for i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j), α2 * (Fu[i,j+1] + Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the positive end of the y-direction (where the boundary fields are
                # taken from the negative-end boundary)
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny), α2 * (e⁻ⁱᵏᴸ*Fu[i,1] + Fu[i,Ny]), Val(OP))  # Fu[i,Ny+1] = exp(-i ky Ly) * Fu[i,1]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # y-direction
                @threads for j = 2:Ny-1
                    for i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j), α2 * (Fu[i,j+1] + Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the negative end of the y-direction (where the boundary fields are
                # assumed zero)
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1), α2 * Fu[i,2], Val(OP))  # Fu[i,1] == 0
                end

                # 3. At the positive end of the y-direction (where the boundary fields are
                # assumed zero)
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny), α2 * Fu[i,Ny], Val(OP))  # Fu[i,Ny+1] == 0
                end
            end
        end  # if nw == ...
    else  # backward averaging
        if nw == 1  # w = x
            # 1. At the locations except for the negative end of the x-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for j = 1:Ny
                for i = 2:Nx  # not i = 2:Nx-1
                    @inbounds set_or_add!(Gv, (i,j), α2 * (Fu[i,j] + Fu[i-1,j]), Val(OP))
                end
            end

            # 2. At the negative end of the x-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j), α2 * (Fu[1,j] + Fu[Nx,j]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[0,j] = Fu[Nx,j] / exp(-i kx Lx)
                end
            else  # symmetry boundary
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j), α * Fu[1,j], Val(OP))  # Fu[0,j] = Fu[1,j]; use α instead of α2 as field is doubled
                end
            end
        elseif nw == 2  # w = y
            # 1. At the locations except for the negative end of the y-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for j = 2:Ny
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j), α2 * (Fu[i,j] + Fu[i,j-1]), Val(OP))
                end
            end

            # 2. At the negative end of the y-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1), α2 * (Fu[i,1] + Fu[i,Ny]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[i,0] = Fu[0,Ny] / exp(-i ky Ly)
                end
            else  # symmetry boundary
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1), α * Fu[i,1], Val(OP))  # Fu[i,0] = Fu[i,1]; use α instead of α2 as field is doubled
                end
            end
        end  # if nw == ...
    end  # if isfwd

    return nothing
end

# Arithmetic m̂ for 1D
function apply_m̂!(Gv::AbsArrNumber{1},  # v-component of output field (v = x, y)
                  Fu::AbsArrNumber{1},  # u-component of input field (u = x, y)
                  ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                  nw::Integer,  # 1|2 for averaging along x|y
                  isfwd::Bool,  # true|false for forward|backward averaging
                  isbloch::Bool,  # boundary condition in w-direction
                  e⁻ⁱᵏᴸ::Number=1.0;  # Bloch phase factor
                  α::Number=1.0  # scale factor to multiply to result before adding it to Gv: Gv += α m(Fu)
                  ) where {OP}
    @assert size(Gv)==size(Fu)

    Nx = length(Fu)  # not size(Fu) unlike code for 2D and 3D
    α2 = 0.5 * α

    # Make sure not to include branches inside for loops.
    if isfwd
        if isbloch
            # 1. At locations except for the positive end of the x-direction
            @threads for i = 1:Nx-1
                @inbounds set_or_add!(Gv, (i,), α2 * (Fu[i+1] + Fu[i]), Val(OP))
            end

            # 2. At the positive end of the x-direction (where the boundary fields are taken
            # from the negative-end boundary)
            @inbounds set_or_add!(Gv, (Nx,), α2 * (e⁻ⁱᵏᴸ*Fu[1] + Fu[Nx]), Val(OP))  # Fu[Nx+1] = exp(-i kx Lx) * Fu[1]
        else  # symmetry boundary
            # 1. At the locations except for the positive and negative ends of the x-direction
            @threads for i = 2:Nx-1
                @inbounds set_or_add!(Gv, (i,), α2 * (Fu[i+1] + Fu[i]), Val(OP))
            end

            # 2. At the negative end of the x-direction (where the boundary fields are assumed
            # zero)
            @inbounds set_or_add!(Gv, (1,), α2 * Fu[2], Val(OP))  # Fu[1] == 0

            # 3. At the positive end of the x-direction (where the boundary fields are assumed
            # zero)
            @inbounds set_or_add!(Gv, (Nx,), α2 * Fu[Nx], Val(OP))  # Fu[Nx+1] == 0
        end
    else  # backward averaging
        # 1. At the locations except for the negative end of the x-direction; unlike for the
        # forward difference, for the backward difference this part of the code is common
        # for both the Bloch and symmetry boundary conditions.
        @threads for i = 2:Nx
            @inbounds set_or_add!(Gv, (i,), α2 * (Fu[i] + Fu[i-1]), Val(OP))
        end

        # 2. At the negative end of the x-direction (where the boundary fields are taken
        # from the positive-end boundary for the Bloch boundary condition)
        if isbloch
            @inbounds set_or_add!(Gv, (1,), α2 * (Fu[1] + Fu[Nx]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[0] = Fu[Nx] / exp(-i kx Lx)
        else  # symmetry boundary
            @inbounds set_or_add!(Gv, (1,), α * Fu[1], Val(OP))  # Fu[0] = Fu[1]; use α instead of α2 as field is doubled
        end
    end  # if isfwd

    return nothing
end

# Weighted m̂ for 3D
function apply_m̂!(Gv::AbsArrNumber{3},  # v-component of output field (v = x, y, z)
                  Fu::AbsArrNumber{3},  # u-component of input field (u = x, y, z)
                  ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                  nw::Integer,  # 1|2|3 for averaging along x|y|z
                  isfwd::Bool,  # true|false for forward|backward averaging
                  ∆w::AbsVecNumber,  # line segments to multiply with; vector of length N[nw]
                  ∆w′⁻¹::AbsVecNumber,  # inverse of line segments to divide by; vector of length N[nw]
                  isbloch::Bool,  # boundary condition in w-direction
                  e⁻ⁱᵏᴸ::Number=1.0;  # Bloch phase factor
                  α::Number=1.0  # scale factor to multiply to result before adding it to Gv: Gv += α m(Fu)
                  ) where {OP}
    @assert size(Gv)==size(Fu)
    @assert size(Fu,nw)==length(∆w)
    @assert length(∆w)==length(∆w′⁻¹)

    Nx, Ny, Nz = size(Fu)
    α2 = 0.5 * α

    # Make sure not to include branches inside for loops.
    if isfwd
        if nw == 1  # w = x
            if isbloch
                # 1. At locations except for the positive end of the x-direction
                @threads for k = 1:Nz
                    for j = 1:Ny, i = 1:Nx-1
                        @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[i]) * (∆w[i+1]*Fu[i+1,j,k] + ∆w[i]*Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the positive end of the x-direction (where the boundary fields are
                # taken from the negative-end boundary)
                β = α2 * ∆w′⁻¹[Nx]
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j,k), β * (∆w[1]*e⁻ⁱᵏᴸ*Fu[1,j,k] + ∆w[Nx]*Fu[Nx,j,k]), Val(OP))  # Fu[Nx+1,j,k] = exp(-i kx Lx) * Fu[1,j,k]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # x-direction
                @threads for k = 1:Nz
                    for j = 1:Ny, i = 2:Nx-1
                        @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[i]) * (∆w[i+1]*Fu[i+1,j,k] + ∆w[i]*Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the negative end of the x-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[1]
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j,k), β * ∆w[2]*Fu[2,j,k], Val(OP))  # Fu[1,j,k] == 0
                end

                # 3. At the positive end of the x-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[Nx]
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j,k), β * ∆w[Nx]*Fu[Nx,j,k], Val(OP))  # Fu[Nx+1,j,k] == 0
                end
            end
        elseif nw == 2  # w = y
            if isbloch
                # 1. At locations except for the positive end of the y-direction
                @threads for k = 1:Nz
                    for j = 1:Ny-1, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[j]) * (∆w[j+1]*Fu[i,j+1,k] + ∆w[j]*Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the positive end of the y-direction (where the boundary fields are
                # taken from the negative-end boundary)
                β = α2 * ∆w′⁻¹[Ny]
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny,k), β * (∆w[1]*e⁻ⁱᵏᴸ*Fu[i,1,k] + ∆w[Ny]*Fu[i,Ny,k]), Val(OP))  # Fu[i,Ny+1,k] = exp(-i ky Ly) * Fu[i,1,k]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # y-direction
                @threads for k = 1:Nz
                    for j = 2:Ny-1, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[j]) * (∆w[j+1]*Fu[i,j+1,k] + ∆w[j]*Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the negative end of the y-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[1]
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1,k), β * ∆w[2]*Fu[i,2,k], Val(OP))  # Fu[i,1,k] == 0
                end

                # 3. At the positive end of the y-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[Ny]
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny,k), β * ∆w[Ny]*Fu[i,Ny,k], Val(OP))  # Fu[i,Ny+1,k] == 0
                end
            end
        else  # nw == 3; w = z
            if isbloch
                # 1. At locations except for the positive end of the z-direction
                @threads for k = 1:Nz-1
                    for j = 1:Ny, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[k]) * (∆w[k+1]*Fu[i,j,k+1] + ∆w[k]*Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the positive end of the z-direction (where the boundary fields are
                # taken from the negative-end boundary)
                β = α2 * ∆w′⁻¹[Nz]
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,Nz), β * (∆w[1]*e⁻ⁱᵏᴸ*Fu[i,j,1] + ∆w[Nz]*Fu[i,j,Nz]), Val(OP))  # Fu[i,j,Nz+1] = exp(-i kz Lz) * Fu[i,j,1]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # z-direction
                @threads for k = 2:Nz-1
                    for j = 1:Ny, i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[k]) * (∆w[k+1]*Fu[i,j,k+1] + ∆w[k]*Fu[i,j,k]), Val(OP))
                    end
                end

                # 2. At the negative end of the z-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[1]
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,1), β * ∆w[2]*Fu[i,j,2], Val(OP))  # Fu[i,j,1] == 0
                end

                # 3. At the positive end of the z-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[Nz]
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,Nz), β * ∆w[Nz]*Fu[i,j,Nz], Val(OP))  # Fu[i,j,Nz+1] == 0
                end
            end
        end  # if nw == ...
    else  # backward averaging
        if nw == 1  # w = x
            # 1. At the locations except for the negative end of the x-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for k = 1:Nz
                for j = 1:Ny, i = 2:Nx  # not i = 2:Nx-1
                    @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[i]) * (∆w[i]*Fu[i,j,k] + ∆w[i-1]*Fu[i-1,j,k]), Val(OP))
                end
            end

            # 2. At the negative end of the x-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                β = α2 * ∆w′⁻¹[1]
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j,k), β * (∆w[1]*Fu[1,j,k] + ∆w[Nx]*Fu[Nx,j,k]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[0,j,k] = Fu[Nx,j,k] / exp(-i kx Lx)
                end
            else  # symmetry boundary
                β = α * ∆w′⁻¹[1]
                for k = 1:Nz, j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j,k), β * ∆w[1]*Fu[1,j,k], Val(OP))  # Fu[0,j,k] = Fu[1,j,k]
                end
            end
        elseif nw == 2  # w = y
            # 1. At the locations except for the negative end of the y-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for k = 1:Nz
                for j = 2:Ny, i = 1:Nx  # not j = 2:Ny-1
                    @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[j]) * (∆w[j]*Fu[i,j,k] + ∆w[j-1]*Fu[i,j-1,k]), Val(OP))
                end
            end

            # 2. At the negative end of the y-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                β = α2 * ∆w′⁻¹[1]
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1,k), β * (∆w[1]*Fu[i,1,k] + ∆w[Ny]*Fu[i,Ny,k]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[i,0,k] = Fu[0,Ny,k] / exp(-i ky Ly)
                end
            else  # symmetry boundary
                β = α * ∆w′⁻¹[1]
                for k = 1:Nz, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1,k), β * ∆w[1]*Fu[i,1,k], Val(OP))  # Fu[i,0,k] = Fu[i,1,k]
                end
            end
        else  # nw == 3; w = z
            # 1. At the locations except for the negative end of the z-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for k = 2:Nz
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,k), (α2 * ∆w′⁻¹[k]) * (∆w[k]*Fu[i,j,k] + ∆w[k-1]*Fu[i,j,k-1]), Val(OP))
                end
            end

            # 2. At the negative end of the z-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                β = α2 * ∆w′⁻¹[1]
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,1), β * (∆w[1]*Fu[i,j,1] + ∆w[Nz]*Fu[i,j,Nz]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[i,j,0] = Fu[i,j,Nz] / exp(-i kz Lz)
                end
            else  # symmetry boundary
                β = α * ∆w′⁻¹[1]
                for j = 1:Ny, i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j,1), β * ∆w[1]*Fu[i,j,1], Val(OP))  # Fu[i,j,0] = Fu[i,j,Nz]
                end
            end
        end  # if nw == ...
    end  # if isfwd

    return nothing
end

# Weighted m̂ for 2D
function apply_m̂!(Gv::AbsArrNumber{2},  # v-component of output field (v = x, y)
                  Fu::AbsArrNumber{2},  # u-component of input field (u = x, y)
                  ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                  nw::Integer,  # 1|2 for averaging along x|y
                  isfwd::Bool,  # true|false for forward|backward averaging
                  ∆w::AbsVecNumber,  # line segments to multiply with; vector of length N[nw]
                  ∆w′⁻¹::AbsVecNumber,  # inverse of line segments to divide by; vector of length N[nw]
                  isbloch::Bool,  # boundary condition in w-direction
                  e⁻ⁱᵏᴸ::Number=1.0;  # Bloch phase factor
                  α::Number=1.0  # scale factor to multiply to result before adding it to Gv: Gv += α m(Fu)
                  ) where {OP}
    @assert size(Gv)==size(Fu)
    @assert size(Fu,nw)==length(∆w)
    @assert length(∆w)==length(∆w′⁻¹)

    Nx, Ny = size(Fu)
    α2 = 0.5 * α

    # Make sure not to include branches inside for loops.
    if isfwd
        if nw == 1  # w = x
            if isbloch
                # 1. At locations except for the positive end of the x-direction
                @threads for j = 1:Ny
                    for i = 1:Nx-1
                        @inbounds set_or_add!(Gv, (i,j), (α2 * ∆w′⁻¹[i]) * (∆w[i+1]*Fu[i+1,j] + ∆w[i]*Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the positive end of the x-direction (where the boundary fields are
                # taken from the negative-end boundary)
                β = α2 * ∆w′⁻¹[Nx]
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j), β * (∆w[1]*e⁻ⁱᵏᴸ*Fu[1,j] + ∆w[Nx]*Fu[Nx,j]), Val(OP))  # Fu[Nx+1,j] = exp(-i kx Lx) * Fu[1,j]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # x-direction
                @threads for j = 1:Ny
                    for i = 2:Nx-1
                        @inbounds set_or_add!(Gv, (i,j), (α2 * ∆w′⁻¹[i]) * (∆w[i+1]*Fu[i+1,j] + ∆w[i]*Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the negative end of the x-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[1]
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j), β * ∆w[2]*Fu[2,j], Val(OP))  # Fu[1,j] == 0
                end

                # 3. At the positive end of the x-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[Nx]
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (Nx,j), β * ∆w[Nx]*Fu[Nx,j], Val(OP))  # Fu[Nx+1,j] == 0
                end
            end
        else  # nw == 2; w = y
            if isbloch
                # 1. At locations except for the positive end of the y-direction
                @threads for j = 1:Ny-1
                    for i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j), (α2 * ∆w′⁻¹[j]) * (∆w[j+1]*Fu[i,j+1] + ∆w[j]*Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the positive end of the y-direction (where the boundary fields are
                # taken from the negative-end boundary)
                β = α2 * ∆w′⁻¹[Ny]
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny), β * (∆w[1]*e⁻ⁱᵏᴸ*Fu[i,1] + ∆w[Ny]*Fu[i,Ny]), Val(OP))  # Fu[i,Ny+1] = exp(-i ky Ly) * Fu[i,1]
                end
            else  # symmetry boundary
                # 1. At the locations except for the positive and negative ends of the
                # y-direction
                @threads for j = 2:Ny-1
                    for i = 1:Nx
                        @inbounds set_or_add!(Gv, (i,j), (α2 * ∆w′⁻¹[j]) * (∆w[j+1]*Fu[i,j+1] + ∆w[j]*Fu[i,j]), Val(OP))
                    end
                end

                # 2. At the negative end of the y-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[1]
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1), β * ∆w[2]*Fu[i,2], Val(OP))  # Fu[i,1] == 0
                end

                # 3. At the positive end of the y-direction (where the boundary fields are
                # assumed zero)
                β = α2 * ∆w′⁻¹[Ny]
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,Ny), β * ∆w[Ny]*Fu[i,Ny], Val(OP))  # Fu[i,Ny+1] == 0
                end
            end
        end  # if nw == ...
    else  # backward averaging
        if nw == 1  # w = x
            # 1. At the locations except for the negative end of the x-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for j = 1:Ny
                for i = 2:Nx  # not i = 2:Nx-1
                    @inbounds set_or_add!(Gv, (i,j), (α2 * ∆w′⁻¹[i]) * (∆w[i]*Fu[i,j] + ∆w[i-1]*Fu[i-1,j]), Val(OP))
                end
            end

            # 2. At the negative end of the x-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                β = α2 * ∆w′⁻¹[1]
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j), β * (∆w[1]*Fu[1,j] + ∆w[Nx]*Fu[Nx,j]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[0,j] = Fu[Nx,j] / exp(-i kx Lx)
                end
            else  # symmetry boundary
                β = α * ∆w′⁻¹[1]
                for j = 1:Ny
                    @inbounds set_or_add!(Gv, (1,j), β * ∆w[1]*Fu[1,j], Val(OP))  # Fu[0,j] = Fu[1,j]
                end
            end
        elseif nw == 2  # w = y
            # 1. At the locations except for the negative end of the y-direction; unlike for
            # the forward difference, for the backward difference this part of the code is
            # common for both the Bloch and symmetry boundary conditions.
            @threads for j = 2:Ny
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,j), (α2 * ∆w′⁻¹[j]) * (∆w[j]*Fu[i,j] + ∆w[j-1]*Fu[i,j-1]), Val(OP))
                end
            end

            # 2. At the negative end of the y-direction (where the boundary fields are taken
            # from the positive-end boundary for the Bloch boundary condition)
            if isbloch
                β = α2 * ∆w′⁻¹[1]
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1), β * (∆w[1]*Fu[i,1] + ∆w[Ny]*Fu[i,Ny]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[i,0] = Fu[0,Ny] / exp(-i ky Ly)
                end
            else  # symmetry boundary
                β = α * ∆w′⁻¹[1]
                for i = 1:Nx
                    @inbounds set_or_add!(Gv, (i,1), β * ∆w[1]*Fu[i,1], Val(OP))  # Fu[i,0] = Fu[i,1]
                end
            end
        end  # if nw == ...
    end  # if isfwd

    return nothing
end

# Weighted m̂ for 1D
function apply_m̂!(Gv::AbsArrNumber{1},  # v-component of output field (v = x, y)
                  Fu::AbsArrNumber{1},  # u-component of input field (u = x, y)
                  ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                  nw::Integer,  # 1|2 for averaging along x|y
                  isfwd::Bool,  # true|false for forward|backward averaging
                  ∆w::AbsVecNumber,  # line segments to multiply with; vector of length N[nw]
                  ∆w′⁻¹::AbsVecNumber,  # inverse of line segments to divide by; vector of length N[nw]
                  isbloch::Bool,  # boundary condition in w-direction
                  e⁻ⁱᵏᴸ::Number=1.0;  # Bloch phase factor
                  α::Number=1.0  # scale factor to multiply to result before adding it to Gv: Gv += α m(Fu)
                  ) where {OP}
    @assert size(Gv)==size(Fu)
    @assert size(Fu,nw)==length(∆w)
    @assert length(∆w)==length(∆w′⁻¹)

    Nx = length(Fu)  # not size(Fu) unlike code for 2D and 3D
    α2 = 0.5 * α

    # Make sure not to include branches inside for loops.
    if isfwd
        if isbloch
            # 1. At locations except for the positive end of the x-direction
            @threads for i = 1:Nx-1
                @inbounds set_or_add!(Gv, (i,), (α2 * ∆w′⁻¹[i]) * (∆w[i+1]*Fu[i+1] + ∆w[i]*Fu[i]), Val(OP))
            end

            # 2. At the positive end of the x-direction (where the boundary fields are taken
            # from the negative-end boundary)
            β = α2 * ∆w′⁻¹[Nx]
            @inbounds set_or_add!(Gv, (Nx,), β * (∆w[1]*e⁻ⁱᵏᴸ*Fu[1] + ∆w[Nx]*Fu[Nx]), Val(OP))  # Fu[Nx+1] = exp(-i kx Lx) * Fu[1]
        else  # symmetry boundary
            # 1. At the locations except for the positive and negative ends of the x-direction
            @threads for i = 2:Nx-1
                @inbounds set_or_add!(Gv, (i,), (α2 * ∆w′⁻¹[i]) * (∆w[i+1]*Fu[i+1] + ∆w[i]*Fu[i]), Val(OP))
            end

            # 2. At the negative end of the x-direction (where the boundary fields are assumed
            # zero)
            β = α2 * ∆w′⁻¹[1]
            @inbounds set_or_add!(Gv, (1,), β * ∆w[2]*Fu[2], Val(OP))  # Fu[1] == 0

            # 3. At the positive end of the x-direction (where the boundary fields are assumed
            # zero)
            β = α2 * ∆w′⁻¹[Nx]
            @inbounds set_or_add!(Gv, (Nx,), β * ∆w[Nx]*Fu[Nx], Val(OP))  # Fu[Nx+1] == 0
        end
    else  # backward averaging
        # 1. At the locations except for the negative end of the x-direction; unlike for the
        # forward difference, for the backward difference this part of the code is common
        # for both the Bloch and symmetry boundary conditions.
        @threads for i = 2:Nx
            @inbounds set_or_add!(Gv, (i,), (α2 * ∆w′⁻¹[i]) * (∆w[i]*Fu[i] + ∆w[i-1]*Fu[i-1]), Val(OP))
        end

        # 2. At the negative end of the x-direction (where the boundary fields are taken
        # from the positive-end boundary for the Bloch boundary condition)
        if isbloch
            β = α2 * ∆w′⁻¹[1]
            @inbounds set_or_add!(Gv, (1,), β * (∆w[1]*Fu[1] + ∆w[Nx]*Fu[Nx]/e⁻ⁱᵏᴸ), Val(OP))  # Fu[0] = Fu[Nx] / exp(-i kx Lx)
        else  # symmetry boundary
            β = α * ∆w′⁻¹[1]
            @inbounds set_or_add!(Gv, (1,), β * ∆w[1]*Fu[1], Val(OP))  # Fu[0] = Fu[1]
        end
    end  # if isfwd

    return nothing
end


text:
function mkent2elem(elcon)

ne = size(elcon,2);
ent = unique(elcon[:]);
ent = ent[ent.>0];
ndof = length(ent);
entmax = ent[end]; # largest entity

#display([ndof entmax])
#display(size(elcon))
# display(ent)
# display(entmax)

# entity-to-index mapping
ent2ind = zeros(Int,entmax,1);
ent2ind[ent] = 1:ndof;

 # store number of neighboring elements for each entity
rowent2elem = zeros(Int,ndof,1);
for i = 1:ne  # for each element i
    elc = elcon[:,i];     # entities on element i
    k   = unique(elc[:]); # remove duplicate entities on \element i
    k   = k[k.>0];
    ind = ent2ind[k];   # get entity indices
    rowent2elem[ind] = rowent2elem[ind] .+ 1;
end
rowent2elem = [0; cumsum(rowent2elem[:])];

 # store neighboring-element indices for each entity
colent2elem = zeros(Int,rowent2elem[end],1);
inc = ones(Int,ndof,1);
for i = 1:ne
    elc = elcon[:,i];   # entities on element i
    k = unique(elc[:]); # remove duplicate entities on element i
    k   = k[k.>0];
    ind = ent2ind[k];   # get entity indices
    # rowent2elem(ind): pointer to the list of entities k
    colent2elem[rowent2elem[ind].+inc[ind]] .= i;
    inc[ind] = inc[ind] .+ 1;
end

return rowent2elem,colent2elem,ent2ind,ent,ndof;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# tutorial on conventional 2D SLAM example

addprocs(3)

# This tutorial shows how to use some of the commonly used factor types
# This tutorial follows from the ContinuousScalar example from IncrementalInference

@everywhere begin
using RoME, Distributions
using IncrementalInference
using TransformUtils

import IncrementalInference: getSample

# struct Pose2 <: IncrementalInference.InferenceVariable
#   dims::Int
#   Pose2() = new(3)
# end
# struct Point2 <: IncrementalInference.InferenceVariable
#   dims::Int
#   Point2() = new(2)
# end

struct Prior{T} <: IncrementalInference.FunctorSingleton where {T <: Distribution}
  z::T
end
getSample(s::Prior, N::Int=1) = (rand(s.z,N), )

struct Pose2Pose2_NEW{T} <: IncrementalInference.FunctorPairwise where {T <: Distribution}
  z::T
  Pose2Pose2_NEW() where {T <: Distribution} = new{T}()
  Pose2Pose2_NEW(z1::T) where {T <: Distribution} = new{T}(z1)
  Pose2Pose2_NEW(mean::Vector{Float64}, cov::Array{Float64,2}) where {T <: Distribution} = new{Distributions.MvNormal}(MvNormal(mean, cov))
end
getSample(s::Pose2Pose2_NEW, N::Int=1) = (rand(s.z,N), )
function (s::Pose2Pose2_NEW)(res::Array{Float64},
      idx::Int,
      meas::Tuple,
      wxi::Array{Float64,2},
      wxj::Array{Float64,2}  )
  # res[1] = meas[1][idx] - (X2[1,idx] - X1[1,idx])
  wXjhat = SE2(wxi[:,idx])*SE2(meas[1][:,idx]) #*SE2(pp2.Zij[:,1])*SE2(meas[1][:,idx])
  jXjhat = SE2(wxj[:,idx]) \ wXjhat
  se2vee!(res, jXjhat)
  nothing
end
end # everywhere

# Start with an empty graph
fg = initfg(sessionname="SLAM2D_TUTORIAL")

# also add a PriorPose2 to pin the first pose at a fixed location
addNode!(fg, :x0, Pose2, labels=["POSE"])
addFactor!(fg, [:x0], Prior(MvNormal([0.0;0.0;0], diagm([1.0;1.0;0.01].^2))) )

# Drive around in a hexagon
for i in 0:5
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  addNode!(fg, nsym, Pose2, labels=["POSE"])
  addFactor!(fg, [psym;ensureAllInitialized!nsym], Pose2Pose2_NEW(MvNormal([10.0;0;pi/3], diagm([0.1;0.1;0.1].^2))))
end

# Graphs.plot(fg.g)

# Add landmarks with Bearing range measurements
addNode!(fg, :l1, Point2, labels=["LANDMARK"])
p2br = Pose2DPoint2DBearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x0; :l1], p2br)

# Add landmarks with Bearing range measurements
p2br2 = Pose2DPoint2DBearingRange(Normal(0,0.1),Normal(20.0,1.0))
addFactor!(fg, [:x6; :l1], p2br2)

isInitialized(fg, :l1)

# Replace this call with an initialization service that traverses the entire graph
ensureAllInitialized!(fg)

# get all initialized nodes

using RoMEPlotting, Gadfly

pl = plotKDE(fg, [:x0; :x1; :x2; :x3; :x4; :x5; :x6]);

Gadfly.draw(PDF("tmpX0123456.pdf", 15cm, 15cm), pl)

@async run(`evince tmpX0123456.pdf`)

# pl = drawPoses(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(PDF("tmpPosesFg.pdf", 15cm, 15cm), pl)
@async run(`evince tmpPosesFg.pdf`)

# Perform inference over the entire graph

tree = wipeBuildNewTree!(fg)
# @async Graphs.plot(tree.bt)
@time inferOverTree!(fg, tree)

#


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using JuMP
using MIPLearn

Base.@kwdef struct KnapsackData
    weights = [1.0, 2.0, 3.0]
    prices = [5.0, 6.0, 7.0]
    capacity = 3.0
end

function build_knapsack_model(data = KnapsackData())
    model = Model()
    n = length(data.weights)
    @variable(model, x[1:n], Bin)
    @objective(model, Max, sum(x[i] * data.prices[i] for i = 1:n))
    @constraint(model, c1, sum(x[i] * data.weights[i] for i = 1:n) <= data.capacity)

    # # Add ML information to the model
    # @feature(model, [5.0])
    # @feature(c1, [1.0, 2.0, 3.0])
    # @category(c1, "c1")
    # for i = 1:n
    #     @feature(x[i], [weights[i]; prices[i]])
    #     @category(x[i], "type-$i")
    # end

    # # Should store ML information
    # @test model.ext[:miplearn]["variable_features"]["x[1]"] == [1.0, 5.0]
    # @test model.ext[:miplearn]["variable_features"]["x[2]"] == [2.0, 6.0]
    # @test model.ext[:miplearn]["variable_features"]["x[3]"] == [3.0, 7.0]
    # @test model.ext[:miplearn]["variable_categories"]["x[1]"] == "type-1"
    # @test model.ext[:miplearn]["variable_categories"]["x[2]"] == "type-2"
    # @test model.ext[:miplearn]["variable_categories"]["x[3]"] == "type-3"
    # @test model.ext[:miplearn]["constraint_features"]["c1"] == [1.0, 2.0, 3.0]
    # @test model.ext[:miplearn]["constraint_categories"]["c1"] == "c1"
    # @test model.ext[:miplearn]["instance_features"] == [5.0]

    return model
end

function build_knapsack_file_instance()
    data = KnapsackData()
    filename = tempname()
    MIPLearn.save_data(filename, data)
    return FileInstance(filename, build_knapsack_model)
end


text:
"""
Noop index type for storing vectors. Returns empty vectors of indexes, scores.
Useful when search is done only in the db.
"""
struct NoopIndex{E} <: Garamond.AbstractIndex
    length::Int
end

NoopIndex(data, args...; kwargs...) = NoopIndex{eltype(data)}(size(data,2))

# Nearest neighbor search method
function knn_search(index::NoopIndex{E},
                    point::AbstractVector,
                    k::Int,
                    keep::Vector{Int}=collect(1:length(index));
                    kwargs...
                   ) where {E<:AbstractFloat}
    Vector{Int}(), Vector{E}()
end

# Length method
length(index::NoopIndex) = index.length

# push!, pushfirst!, pop!, popfirst!, delete_from_index!
Base.push!(index::NoopIndex, point) = nothing

Base.pushfirst!(index::NoopIndex, point) = nothing

Base.pop!(index::NoopIndex) = nothing

Base.popfirst!(index::NoopIndex) = nothing

delete_from_index!(index::NoopIndex, pos) = nothing


text:
# Returns a partitioning of `range` into sub-ranges of length `partition_size`.
# If `range` is not divisible by `partition_size`, the last sub-range will be
# of length `range` mod `partition_size`.
function partition_range(range::UnitRange, partition_size::Int)
    n_partitions = ceil(Int, last(range) / partition_size)
    partitions = Vector{UnitRange{Int}}(undef, n_partitions)

    for p in 1:n_partitions
        l::Int = first(range) + (p-1) * partition_size
        r::Int = min(l - 1 + partition_size, last(range))
        partitions[p] = l:r
    end

    return partitions
end


text:
function isPrime(n::Int)
	for i in 2:ceil(sqrt(n))
		if n%i==0 && n!=i return false end
	end
	return true
end

function nextPrime!(l)
	x = l[end]+2
	while !isPrime(x)
		x += 2
	end
	push!(l, x)
end

l = [2,3,5,7,11,13]
while length(l) < 10001
	nextPrime!(l)
end
println(length(l))
println(l[end])

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
#! /usr/bin/julia

# Rosetta Code, Penney's game

const SLEN = 3

autobet() = randbool(SLEN)
function autobet(ob::BitArray{1})
    rlen = length(ob)
    2 < rlen || return ~ob
    3 < rlen || return [~ob[2], ob[1:2]]
    opt = falses(rlen)
    opt[1] = true
    opt[end-1:end] = true
    ob != opt || return ~opt
    return opt
end
autobet(ob::Array{Bool,1}) = autobet(convert(BitArray{1}, ob))

function pgencode{T<:String}(a::T)
    b = uppercase(a)
    0 < length(b) || return trues(0)
    !ismatch(r"[^HT]+", b) || error(@sprintf "%s is not a HT sequence" a)
    b = split(b, "")
    b .== "H"
end
pgdecode(a::BitArray{1}) = join([i ? "H" : "T" for i in a], "")

function humanbet()
    b = ""
    while length(b) != SLEN  || ismatch(r"[^HT]+", b)
        print("Your bet? ")
        b = uppercase(chomp(readline()))
    end
    return b
end

println("Playing Penney's Game Against the computer.")

if randbool()
    mach = autobet()
    println(@sprintf "The computer bet first, chosing %s." pgdecode(mach))
    println("Now you can bet.")
    human = pgencode(humanbet())
else
    println("You bet first.")
    human = pgencode(humanbet())
    mach = autobet(human)
end
print(@sprintf "You bet %s " pgdecode(human))
println(@sprintf "and the computer bet %s." pgdecode(mach))

pg = randbool(SLEN)
pgtail = copy(pg)
while pgtail != mach && pgtail != human
    push!(pg, randbool())
    pgtail = [pgtail[2:end], pg[end]]
end

println(@sprintf("This game lasted %d turns yielding\n    %s",
                 length(pg), pgdecode(pg)))

if human == mach
    println("so you and the computer tied.")
elseif pgtail == mach
    println("so the computer won.")
else
    println("so you won.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
coupler_atmos_boundary_flux(_, ..) = nothing

# momentum
struct CouplerVectorFlux{C, F} <: AbstractBoundary
    coefficients::C
    flux::F
end
function get_boundary_flux(model, bc::CouplerVectorFlux, var::Fields.Field, Y, Ya) # TODO: adapt for Field BC (CC#325)
    flux = Geometry.Covariant3Vector(FT(1)) ⊗ Geometry.Covariant12Vector(FT(1), FT(1)) * parent(bc.flux)[1]
end
"""
get_boundary_flux(
    model,
    bc::CouplerVectorFlux,
    ρc::Fields.Field,
    Y,
    Ya,
)

This is an extenison of ClimaAtmos.get_boundary_flux for momentum
"""
function coupler_atmos_boundary_flux(bc::CouplerVectorFlux, atmos_sim, slab_sim, F_a_space, F_s_space)

    Y_atm = atmos_sim.integrator.u
    FT = eltype(Y_atm)

    uv = Y_atm.base.uh

    uh1_cov = Fields.level(uv.components.data.:1, 1) # TODO: change to Field - pending ClimaCore PR #325
    uh2_cov = Fields.level(uv.components.data.:2, 1) # TODO: change to Field - pending ClimaCore PR #325

    # physical scale (wind * coeff)
    uv_1 = Fields.level(Geometry.UVVector.(uv), 1)
    u_wind = norm(uv_1) # TODO: this will need to be spatially variable

    # unit vector in W direction on central space
    normal_unit_vector = Fields.level(similar(Geometry.WVector.(Y_atm.thermodynamics.ρe_tot)), 1) # TODO: for topography, needs to be perp to surface
    parent(normal_unit_vector) .= FT(1)
    local_geometry = Fields.local_geometry_field(axes(normal_unit_vector))

    # contravariant scale
    scale_con =
        Geometry.Contravariant3Vector.(
            Geometry.WVector.(normal_unit_vector .* bc.coefficients .* u_wind),
            local_geometry,
        )

    parent(bc.flux) .=
        parent(Geometry.Contravariant3Vector.(scale_con) .⊗ Geometry.Covariant12Vector.(uh1_cov, uh2_cov))
end

# energy
struct CouplerEnergyFlux{C, F} <: AbstractBoundary
    coefficients::C
    flux::F
end

"""
get_boundary_flux(
    model,
    bc::CouplerEnergyFlux,
    ρc::Fields.Field,
    Y,
    Ya,
)

This is an extenison of ClimaAtmos.get_boundary_flux for enthalpy
"""
function get_boundary_flux(model, bc::CouplerEnergyFlux, ρc::ClimaCore.Fields.Field, Y, Ya)
    bulk_flux = bc.flux
    ClimaCore.Geometry.WVector.(parent(bulk_flux)[1]) #flux = Geometry.Contravariant3Vector.(bulk_flux) # need CC extension - waiting for PR (https://github.com/CliMA/ClimaCore.jl/pull/325)
end

"""
coupler_boundary_flux(
    model,
    bc::CouplerEnergyFlux,
    ρc::Fields.Field,
    Y,
    Ya,
)

Vertical fluxes for arbitrary variables with the bulk aerodynamic turbulent formula,
given variable exchange coefficients. (e.g. for energy, or tracers)
Currently this is a simplified version of the bulk formula for testing.
"""
function coupler_atmos_boundary_flux(bc::CouplerEnergyFlux, atmos_sim, slab_sim, F_a_space, F_s_space)

    Y_atm = atmos_sim.integrator.u
    FT = eltype(Y_atm)

    T_1 = ClimaCore.Fields.level(calculate_temperature(atmos_sim, FT), 1)

    ρ_1 = ClimaCore.Fields.level(Y_atm.base.ρ, 1)

    # remap slab > atmos
    T_sfc_atmos_grid = ClimaCore.Fields.zeros(F_a_space)
    T_sfc_slab_grid = copy(slab_sim.integrator.u.T_sfc)
    # ClimaCoreTempestRemap.remap!(T_sfc_atmos_grid, T_sfc_slab_grid, R_slab2atm)
    dummmy_remap!(T_sfc_slab_grid, T_sfc_atmos_grid)

    # save in atmos BCs
    c_pd = CLIMAParameters.Planet.cp_d(atmos_sim.model.parameters)
    # parent(bc.flux) .= bc.coefficients .* c_pd .* parent(ρ_1) .* (parent(T_1) .- parent(T_sfc_atmos_grid)) # parent(bc.flux) .= bc.coefficients .* parent(ρ_1 .* u_wind) .* (parent(c_1) .- parent(c_sfc)) # TODO: make neater - same space but different instance error otherwise
    parent(bc.flux) .= FT(10)
    # TODO: make neater - same space but different instance error otherwise
    nothing
end

# code below should be moved to CA / thermodyn
using LinearAlgebra: norm_sqr
using Thermodynamics

function calculate_temperature(atmos_sim, FT)
    model = atmos_sim.model
    Y = atmos_sim.integrator.u
    p = calculate_pressure(Y, Y, model.thermodynamics, model.moisture, model.moisture, model.parameters, FT)
    R = FT(287)
    T = p ./ R ./ Y.base.ρ
end

@inline function calculate_gravitational_potential(Y, Ya, params, FT)
    g::FT = CLIMAParameters.Planet.grav(params)
    ρ = Y.base.ρ
    z = Fields.coordinate_field(axes(ρ)).z

    return @. g * z
end

@inline function calculate_pressure(Y, Ya, i, ii, iii, params, FT)
    ρ = Y.base.ρ
    uh = Y.base.uh
    w = Y.base.w
    ρe_tot = Y.thermodynamics.ρe_tot

    interp_f2c = Operators.InterpolateF2C()

    z = Fields.coordinate_field(axes(ρ)).z
    K = calculate_kinetic_energy(Y, Y, params, FT)
    Φ = calculate_gravitational_potential(Y, Y, params, FT)

    e_int = @. ρe_tot / ρ - Φ - K
    p = Thermodynamics.air_pressure.(Thermodynamics.PhaseDry.(params, e_int, ρ))

    return p
end

@inline function calculate_kinetic_energy(Y, Ya, params, FT)
    cuₕ = Y.base.uh # Covariant12Vector on centers
    fw = Y.base.w # Covariant3Vector on faces
    If2c = Operators.InterpolateF2C()

    cuvw = Geometry.Covariant123Vector.(cuₕ) .+ Geometry.Covariant123Vector.(If2c.(fw))
    cK = @. norm_sqr(cuvw) / 2
    return cK
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
module RandomFourierFeatures

using KernelFunctions
using Distributions
using LinearAlgebra
using FillArrays

export sample_basis, gp_rff_approx

include("random_fourier_features.jl")

end


text:
using Documenter
using AlphaZero

const PRETTY_URLS = get(ENV, "CI", nothing) == "true"

makedocs(
  sitename = "AlphaZero",
  authors="Jonathan Laurent",
  format = Documenter.HTML(prettyurls=PRETTY_URLS),
  modules = [AlphaZero],
  pages = [
    "Home" => "index.md",
    "Tutorial" => [
      "tutorial/alphazero_intro.md",
      "tutorial/connect_four.md"
      #"tutorial/codebase_overview.md"
    ],
    "Reference" => [
      "reference/params.md",
      "reference/game_interface.md",
      "reference/mcts.md",
      "reference/network.md",
      "reference/networks_library.md",
      "reference/player.md",
      "reference/memory.md",
      "reference/environment.md",
      "reference/benchmark.md",
      "reference/session.md",
      "reference/reports.md",
    ],
    "Contributing" => [
      "contributing/guide.md",
      "contributing/add_game.md"
    ]
  ],
  repo="https://github.com/jonathan-laurent/AlphaZero.jl/blob/{commit}{path}#L{line}"
)

deploydocs(;
  repo="github.com/jonathan-laurent/AlphaZero.jl",
)


text:
using Test
using Dash
using Dash: ApplicationResources, main_registry, HandlerState, make_handler
using HTTP

@testset "external_stylesheets" begin
    app = dash()
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)

    @test isnothing(findfirst("link rel=\"stylesheet\"", index_page))

    app = dash(external_stylesheets = ["https://test.css"])
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test !isnothing(findfirst("<link rel=\"stylesheet\" href=\"https://test.css\">", index_page))

    app = dash(external_stylesheets = [
        "https://test.css",
         Dict("href" => "https://test2.css", "rel" => "stylesheet")
         ]
         )
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)

    @test !isnothing(findfirst("<link rel=\"stylesheet\" href=\"https://test.css\">", index_page))
    @test !isnothing(findfirst("href=\"https://test2.css\"", index_page))

end

@testset "external_scripts" begin

    app = dash(external_scripts = ["https://test.js"])
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test !isnothing(findfirst("""<script src="https://test.js"></script>""", index_page))

    app = dash(external_scripts = [
        "https://test.js",
        Dict("src" => "https://test2.js", "crossorigin" => "anonymous")
        ])
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)

    @test !isnothing(findfirst("""<script src="https://test.js"></script>""", index_page))

    @test !isnothing(findfirst("""<script src="https://test2.js" crossorigin="anonymous"></script>""", index_page))

end

@testset "url paths" begin
    app = dash(requests_pathname_prefix = "/reg/prefix/", routes_pathname_prefix = "/prefix/")
    @test app.config.requests_pathname_prefix == "/reg/prefix/"
    @test app.config.routes_pathname_prefix == "/prefix/"

    app = dash(routes_pathname_prefix = "/prefix/")
    @test app.config.routes_pathname_prefix == "/prefix/"
    @test app.config.requests_pathname_prefix == "/prefix/"

    app = dash()
    @test app.config.routes_pathname_prefix == "/"
    @test app.config.requests_pathname_prefix == "/"

    ENV["DASH_APP_NAME"] = "test-app"
    app = dash(routes_pathname_prefix = "/prefix/")
    @test app.config.routes_pathname_prefix == "/prefix/"
    @test app.config.requests_pathname_prefix == "/test-app/prefix/"

    app = dash()
    @test app.config.routes_pathname_prefix == "/"
    @test app.config.requests_pathname_prefix == "/test-app/"
    delete!(ENV, "DASH_APP_NAME")

end

@testset "assets paths" begin
    app = dash()
    app.layout = html_div()
    handler = make_handler(app)
    request = HTTP.Request("GET", "/assets/test.png")
    res = HTTP.handle(handler, request)
    @test res.status == 200
    request = HTTP.Request("GET", "/assets/test3.png")
    res = HTTP.handle(handler, request)
    @test res.status == 404
    request = HTTP.Request("GET", "/images/test.png")
    res = HTTP.handle(handler, request)
    @test startswith(HTTP.header(res, "Content-Type"), "text/html")

    app = dash(url_base_pathname = "/test/")
    app.layout = html_div()
    handler = make_handler(app)
    request = HTTP.Request("GET", "/assets/test.png")
    res = HTTP.handle(handler, request)
    @test res.status == 404
    request = HTTP.Request("GET", "/test/assets/test.png")
    res = HTTP.handle(handler, request)
    @test res.status == 200
    request = HTTP.Request("GET", "/images/test.png")
    res = HTTP.handle(handler, request)
    @test res.status == 404

    app = dash(assets_url_path = "ass")
    app.layout = html_div()
    handler = make_handler(app)
    request = HTTP.Request("GET", "/ass/test.png")
    res = HTTP.handle(handler, request)
    @test res.status == 200
    request = HTTP.Request("GET", "/ass/test3.png")
    res = HTTP.handle(handler, request)
    @test res.status == 404
    request = HTTP.Request("GET", "/assets/test3.png")
    res = HTTP.handle(handler, request)
    @test startswith(HTTP.header(res, "Content-Type"), "text/html")
    request = HTTP.Request("GET", "/images/test.png")
    res = HTTP.handle(handler, request)
    @test startswith(HTTP.header(res, "Content-Type"), "text/html")

    app = dash(assets_folder = "images")
    app.layout = html_div()
    handler = make_handler(app)
    request = HTTP.Request("GET", "/assets/test.png")
    res = HTTP.handle(handler, request)
    @test res.status == 404
    request = HTTP.Request("GET", "/assets/test_images.png")
    res = HTTP.handle(handler, request)
    @test res.status == 200
    request = HTTP.Request("GET", "/images/test.png")
    res = HTTP.handle(handler, request)
    @test startswith(HTTP.header(res, "Content-Type"), "text/html")
end

@testset "suppress_callback_exceptions" begin
    app = dash()
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test !isnothing(findfirst("\"suppress_callback_exceptions\":false", index_page))
    @test isnothing(findfirst("\"suppress_callback_exceptions\":true", index_page))

    app = dash(suppress_callback_exceptions = true)
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test isnothing(findfirst("\"suppress_callback_exceptions\":false", index_page))
    @test !isnothing(findfirst("\"suppress_callback_exceptions\":true", index_page))
end

@testset "meta_tags" begin
    app = dash()
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)

    @test !isnothing(
        findfirst(
            "<meta http-equiv=\"X-UA-Compatible\" content=\"IE=edge\">",
            index_page)
        )

    @test !isnothing(
        findfirst(
            "<meta charset=\"UTF-8\">",
            index_page)
        )

    app = dash(meta_tags = [Dict("type" => "tst", "rel" => "r")])
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)

    @test !isnothing(
            findfirst(
                "<meta http-equiv=\"X-UA-Compatible\" content=\"IE=edge\">",
                index_page)
            )

    @test !isnothing(
            findfirst(
                "<meta charset=\"UTF-8\">",
                index_page)
            )

    @test !isnothing(
        findfirst(
            Dash.format_tag("meta", Dict("type" => "tst", "rel" => "r"), opened = true),
            index_page)
        )

 app = dash(meta_tags = [Dict("charset" => "Win1251"), Dict("type" => "tst", "rel" => "r")])
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)

    @test isnothing(
            findfirst(
                "<meta charset=\"UTF-8\">",
                index_page)
            )
    @test !isnothing(
        findfirst(
            "<meta charset=\"Win1251\">",
            index_page)
        )

    @test !isnothing(
        findfirst(
            Dash.format_tag("meta", Dict("type" => "tst", "rel" => "r"), opened = true),
            index_page)
        )

    app = dash(meta_tags = [Dict("http-equiv" => "X-UA-Compatible", "content" => "IE"), Dict("type" => "tst", "rel" => "r")])
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test isnothing(
        findfirst(
            "<meta http-equiv=\"X-UA-Compatible\" content=\"IE=edge\">",
            index_page)
        )
    @test !isnothing(
            findfirst(
                Dash.format_tag("meta", Dict("http-equiv" => "X-UA-Compatible", "content" => "IE"), opened = true),
                index_page)
            )

end

@testset "index_string" begin
    index_string = "test test test, {%metas%},{%title%},{%favicon%},{%css%},{%app_entry%},{%config%},{%scripts%},{%renderer%}"
    app = dash(index_string = index_string)
    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test startswith(index_page, "test test test,")

end

@testset "show_undo_redo" begin

    app = dash()

    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test !isnothing(findfirst("\"show_undo_redo\":false", index_page))

    app = dash(show_undo_redo = true)

    resources = ApplicationResources(app, main_registry())
    index_page = Dash.index_page(app, resources)
    @test !isnothing(findfirst("\"show_undo_redo\":true", index_page))
end


text:
#   The extended Woods problem.

#   This problem is a sum of n/4 sets of 6 terms, each of which is
#   assigned its own group.  For a given set i, the groups are
#   A(i), B(i), C(i), D(i), E(i) and F(i). Groups A(i) and C(i) contain 1
#   nonlinear element each, denoted Y(i) and Z(i).

#   The problem dimension is defined from the number of these sets.
#   The number of problem variables is then 4 times larger.

#   This version uses a slightly unorthodox expression of Woods
#   function as a sum of squares (see Buckley)

#   Source:  problem 14 in
#   J.J. More', B.S. Garbow and K.E. Hillstrom,
#   "Testing Unconstrained Optimization Software",
#   ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981.

#   See also Toint#27, Buckley#17 (p. 101), Conn, Gould, Toint#7

#   classification SUR2-AN-V-0

# See also
#
#   problem 58 in
#   L. Luksan, C. Matonoha and J. Vlcek
#   Modified CUTE problems for sparse unconstrained optimization,
#   Technical Report 1081,
#   Institute of Computer Science,
#   Academy of Science of the Czech Republic
#
#   http://www.cs.cas.cz/matonoha/download/V1081.pdf
#
#   classification SUR2-AN-V-0
#
# J.-P. Dussault, Clermont-Ferrand, 05/2016.

export woods

"The extended Woods problem `n` "
function woods(n :: Int=100)

    (n % 4 == 0) || warn("woods: number of variables adjusted to be a multiple of 4")
    n = 4 * max(1, div(n, 4))

    nlp = Model()

    x0 = -3 * ones(n)
    x0[2*(collect(1:div(n,2)))] .= -1.0

    @variable(nlp, x[i=1:n], start=-2)

    @NLobjective(
      nlp,
      Min,
      1.0 + sum(
        100 * (x[4*i-2] - x[4*i-3]^2)^2 + (1 - x[4*i-3])^2 +
        90 * (x[4*i] - x[4*i-1]^2)^2 + (1 - x[4*i-1])^2 +
        10 * (x[4*i-2] + x[4*i] - 2)^2 + 0.1 * (x[4*i-2] - x[4*i])^2 for i=1:div(n,4))
    )

    return nlp
end


text:
n = rand(1:10)
@test matrixdepot("binomial", n) == matrixdepot("binomial", Float64, n)

A = matrixdepot("binomial", n)

@test_approx_eq A*A 2^(n-1)*eye(n)
println("'binomial' passed test...")


text:
export SCPIInstrument, SCPIInstrumentError, command, query

struct SCPIInstrumentError <: Exception
    msg::String
end

Base.showerror(io::IO, e::SCPIInstrumentError) = print(io, e.msg)

abstract type SCPIInstrument end

@mustimplement command(inst::SCPIInstrument, command_::String)
@mustimplement query(inst::SCPIInstrument, query_::String)

include("SerialSCPIInstrument.jl")
include("TCPSCPIInstrument.jl")


text:
using SegwayRobot

dt = 1/25
t = 0:dt:20

# Simulate
@time for j = 1:200
    xt = randn(7)*0.01;
    for i = 1:length(t)
        segway_rk4!(xt, (2*xt[1]+50*xt[4]+30*xt[5])*ones(2), dt)
    end
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
confoundersfile = joinpath("data", "confounders.csv")
queryfile = joinpath("config", "query.toml")
phenotypefile = joinpath("data", "phenotypes.csv")

function build_query_file(;path=queryfile, threshold=0.9, het="AG")
    queries = Dict(
        "threshold" => threshold,
        "SNPS" => Dict(
            "RSID_10"  => BGEN.datadir("example.8bits.bgen"),
            "RSID_100" => BGEN.datadir("example.8bits.bgen")),
        "QUERY_1" => Dict(
            "RSID_10"  => het*" -> GG",
            "RSID_100" => het*" -> GG"),
        "QUERY_2" => Dict(
            "RSID_10"  => het*" -> GG",
            "RSID_100" => "AA -> GG"
        )
    )
    open(path, "w") do io
        TOML.print(io, queries)
    end
end

function build_ate_query_file(;path=queryfile, threshold=0.9)
    queries = Dict(
        "threshold" => threshold,
        "SNPS" => Dict("RSID_10"  => BGEN.datadir("example.8bits.bgen")),
        "QUERY_1" => Dict("RSID_10"  => "AG -> GG"),
        "QUERY_2" => Dict("RSID_10"  => "AA -> GG")
    )
    open(path, "w") do io
        TOML.print(io, queries)
    end
end


text:
using MAT
using DataStructures
using StatsBase
using SparseArrays

include("../src/pdensity_algorithms.jl")

snaps = ["amazon-communities","ca-AstroPhA","condmat2005A","loc-Brightkite","EmailEnronA","roadNet-CA","roadNet-TX","web-Stanford","com-Youtube","web-Google","web-BerkStan"]

Ps = [collect(0.5:0.5:2); 1.05]
for i = 1:length(snaps)
    graph = snaps[i]
    F = matread("../data/snap/standard_$graph.mat")
    A = F["A"]
    n = size(A,1)
    @assert(SparseArrays.issymmetric(A))
    Ssets = zeros(n,length(Ps))
    Times = zeros(length(Ps))
    Objectives = zeros(length(Ps))
    Sizes = zeros(length(Ps))
    EdgeDens = zeros(length(Ps))
    MeanDegree = zeros(length(Ps))
    MaxDegree = zeros(length(Ps))
    MinDegree = zeros(length(Ps))
    Mean2Degree = zeros(length(Ps))
    S, density, degeneracy, maxcore, Objs, ranking, degrees = SimplePeel(A)
    for i = 1:length(Ps)
        p = Ps[i]
        tic = time()
        S, objS, Objs, ranking = GenPeel(A,p)
        timers = time()-tic
        Times[i] = timers
        Ssets[S,i] .= 1
        Objectives[i] = maximum(Objs)
        Sizes[i] = length(S)

        As = A[S,S]
        ds = vec(sum(As,dims = 2))
        edges = sum(As)/2
        posedges = binomial(length(S),2)
        EdgeDens[i] = edges/posedges
        MeanDegree[i] = StatsBase.mean(ds)
        MaxDegree[i] = maximum(ds)
        Mean2Degree[i] = StatsBase.mean(ds.^2)
        MinDegree[i] = minimum(ds)
        println("$graph $n \t $p \t $(length(S)) \t $(StatsBase.mean(ds)) \t $(EdgeDens[i])")
    end

    matwrite("Output/$(graph)_output.mat",Dict(
    "Ps"=>Ps,"Times"=>Times,"Objectives"=>Objectives,"MinDegree"=>MinDegree,
    "Mean2Degree"=>Mean2Degree,
    "Sizes"=>Sizes,"MeanDegree"=>MeanDegree,"MaxDegree"=>MaxDegree,
    "EdgeDens"=>EdgeDens,"n"=>n,"degeneracy"=>degeneracy,"maxcore"=>maxcore))
    matwrite("Output/$(graph)_sets.mat", Dict("Ssets"=>Ssets))
end

## Get stats for the maxcore as well

for i = 1:length(snaps)
    graph = snaps[i]
    F = matread("../data/snap/standard_$graph.mat")
    A = F["A"]
    n = size(A,1)
    M = matread("Output/$(graph)_output.mat")
    S =  M["maxcore"]

    As = A[S,S]
    ds = vec(sum(As,dims = 2))
    edges = sum(As)/2
    posedges = binomial(length(S),2)
    EdgeDen = edges/posedges
    MeanDeg = StatsBase.mean(ds)
    MaxDeg = maximum(ds)
    Mean2 = StatsBase.mean(ds.^2)
    MinDeg = minimum(ds)
    stats = [length(S) EdgeDen MinDeg MeanDeg Mean2 MaxDeg]

    matwrite("Output/$(graph)_maxcore_stats.mat",Dict("stats"=>stats))

end

max_stars_repo_name: nveldt/GenMeanDSG

max_issues_repo_name: nveldt/GenMeanDSG

max_forks_repo_name: nveldt/GenMeanDSG


text:
using DataFrames

"""
    LabeledValues

A `LabeledValues` object represents an efficient way to deal with
arrays of `Label`ed values (`Float64`). Each value is indexed
via a `Label` so storing them in `Vector`s with be O(N) whereas
storing the `Label` and value in a `Dict` is O(1).  However, storing
them in a `Dict` loses the ordering which we would also like
to retain. `LabeledValues` are O(1) while retaining the order of
the labels.
"""
struct LabeledValues
    values::Vector{Float64}
    index::Dict{<:Label,Int}

    function LabeledValues(labels::AbstractVector{<:Label}, values::AbstractVector{Float64})
        @assert length(unique(labels)) == length(labels) "The labels [$labels] are not all unique."
        index = Dict{Label,Int}(labels[i] => i for i in eachindex(labels))
        return new(copy(values), index)
    end
    function LabeledValues(
        labels::NTuple{N,Label},
        values::NTuple{N,Float64},
    ) where {N<:Any}
        @assert length(unique(labels)) == length(labels) "The labels are not all unique."
        index = Dict{Label,Int}(labels[i] => i for i in eachindex(labels))
        return new([values...], index)
    end
    function LabeledValues(prs::Pair{<:Label,Float64}...)
        labels = [pr.first for pr in prs]
        values = [pr.second for pr in prs]
        return LabeledValues(labels, values)
    end
end

Base.eachindex(lv::LabeledValues) = eachindex(lv.values)

Base.show(io::IO, lv::LabeledValues) = print(
    io,
    "LabeledValues[\n" *
    join(("\t" * repr(lbl) * " => " * repr(lv[lbl]) for lbl in labels(lv)), "\n") *
    "\n]",
)

Base.copy(lv::LabeledValues) = LabeledValues(copy(lv.values), copy(lv.index))

Base.merge(lvs::LabeledValues...) =
    LabeledValues(vcat(map(labels, lvs)...), vcat(map(values, lvs)...))

"""
    value(lv::LabeledValues, lbl::Label)::Float64

Returns the value associated with the specified `Label`.
"""
value(lv::LabeledValues, lbl::Label)::Float64 = lv.values[indexin(lv, lbl)]

Base.length(lv::LabeledValues) = length(lv.values)
Base.getindex(lv::LabeledValues, lbl::Label)::Float64 = lv.values[indexin(lv, lbl)]
function Base.get(lv::LabeledValues, lbl::Label, def)
    idx = get(lv.index, lbl, nothing)
    return isnothing(idx) ? def : lv.values[idx]
end

function Base.setindex!(lv::LabeledValues, val::Real, lbl::Label)
    @assert !Base.haskey(lv.index, lbl) "$lbl already exists in LabeledValues - $lv"
    push!(lv.values, val)
    lv.index[lbl] = length(lv.values)
end

"""
    indexin(lv::LabeledValues, lbl::Label)::Int

Returns the index associated with the specified `Label`.
"""
Base.indexin(lv::LabeledValues, lbl::Label)::Int = lv.index[lbl]

Base.haskey(lv::LabeledValues, lbl::Label) = haskey(lv.index, lbl)
Base.keys(lv::LabeledValues) = keys(lv.index)

"""
    labels(lv::LabeledValues)::Vector{Label}

A `Vector` of the `Label`s in order.
"""
function labels(lv::LabeledValues)::Vector{Label}
    res = Array{Label}(undef, length(lv.index))
    for (lbl, idx) in lv.index
        res[idx] = lbl
    end
    return res
end

"""
    values(lv::LabeledValues)::Vector{Float64}

A copy `Vector` of the values in the same order as `labels(lv)`.
"""
Base.values(lv::LabeledValues)::Vector{Float64} = copy(lv.values)

"""
`asa(::Type{DataFrame}, lv::LabeledValues)` extracts a `LabeledValues` object into a `DataFrame` in Label and Value columns.
"""
asa(::Type{DataFrame}, lv::LabeledValues) = DataFrame( Label = labels(lv), Value = values(lv))


text:
export BlockSparseVector
struct BlockSparseVector{B,T,I} <: AbstractSparseVector{T,I}
    elts::SparseVector{SVector{B,T},I}

    function BlockSparseVector{B,T,I}(bx::SparseVector{SVector{B,T},I}) where {B,
                                                                               T,
                                                                               I<:Integer}
        B::Integer
        return new{B,T,I}(bx)
    end
    function BlockSparseVector{B,T}(bx::SparseVector{SVector{B,T},I}) where {B,
                                                                             T,
                                                                             I<:Integer}
        return BlockSparseVector{B,T,I}(bx)
    end
    function BlockSparseVector{B}(bx::SparseVector{SVector{B,T},I}) where {B,T,
                                                                           I<:Integer}
        return BlockSparseVector{B,T,I}(bx)
    end
    function BlockSparseVector(bx::SparseVector{SVector{B,T},I}) where {B,T,
                                                                        I<:Integer}
        return BlockSparseVector{B,T,I}(bx)
    end
end

Base.show(io::IOContext, x::BlockSparseVector) = show(io, SparseVector(x))

###

export blocksparse
function blocksparse(::Val{B}, Is::AbstractVector{I}, Vs::AbstractVector{T},
                     n::Integer) where {B,T,I<:Integer}
    B::Integer
    length(Is) == length(Vs) ||
        throw(ArgumentError("index and value vectors must have the same length"))
    n ≥ 0 || throw(ArgumentError("Number of elements must be non-negative"))
    @assert n % B == 0
    nb = n ÷ B
    Ibs = I[]
    Vbs = SVector{B,T}[]
    curi0 = 0
    v0 = zero(SVector{B,T})
    for (i, v) in zip(Is, Vs)
        @assert 1 ≤ i ≤ n
        i0, i1 = fldmod1(i, B)
        if curi0 ≠ i0
            if curi0 ≠ 0
                push!(Ibs, curi0)
                push!(Vbs, v0)
            end
            curi0 = i0
            v0 = zero(SVector{B,T})
        end
        v0 = setindex(v0, v, i1)
    end
    if curi0 ≠ 0
        push!(Ibs, curi0)
        push!(Vbs, v0)
    end
    return BlockSparseVector(sparsevec(Ibs, Vbs, nb))
end

function blocksparse(::Val{B}, x::AbstractVector{T}) where {B,T}
    n = length(x)
    @assert n % B == 0
    n0 = n ÷ B
    y = Array{SVector{B,T}}(undef, n0)
    for i0 in 1:n0
        y1 = zero(SVector{B,T})
        for i1 in 1:B
            y1 = setindex(y1, x[B * (i0 - 1) + i1], i1)
        end
        y[i0] = y1
    end
    return BlockSparseVector(sparsevec(y))
end

export blocksize
blocksize(x::BlockSparseVector{B}) where {B} = B

for (fun, spfun) in [(:bsprand, :sprand), (:bsprandn, :sprandn)]
    @eval begin
        export $fun
        function $fun(r::AbstractRNG, ::Val{B}, ::Type{T}, n::Integer,
                      p::AbstractFloat) where {B,T}
            @assert n % B == 0
            n0 = n ÷ B
            return BlockSparseVector($spfun(r, SVector{B,T}, n0, p))
        end
        function $fun(r::AbstractRNG, ::Val{B}, n::Integer,
                      p::AbstractFloat) where {B}
            return $fun(r, Val(B), Float64, n, p)
        end
        function $fun(::Val{B}, ::Type{T}, n::Integer,
                      p::AbstractFloat) where {B,T}
            return $fun(Random.default_rng(), Val(B), T, n, p)
        end
        function $fun(::Val{B}, n::Integer, p::AbstractFloat) where {B}
            return $fun(Random.default_rng(), Val(B), n, p)
        end
    end
end

export bspzeros
function bspzeros(::Val{B}, ::Type{T}, ::Type{I},
                  n::Integer) where {B,T,I<:Integer}
    B::Integer
    @assert 0 ≤ n && n % B == 0
    return BlockSparseVector(spzeros(SVector{B,T}, fld1(n, B)))
end
function bspzeros(::Val{B}, ::Type{T}, n::Integer) where {B,T}
    return bspzeros(Val(B), T, Int, n)
end
bspzeros(::Val{B}, n::Integer) where {B} = bspzeros(Val(B), Float64, n)

function (Vector{T})(x::BlockSparseVector{B}) where {T,B}
    r = zeros(T, size(x))
    for j in nzrange(x.elts, 1)
        i0 = rowvals(x.elts)[j]
        v0 = nonzeros(x.elts)[j]
        for i1 in 1:B
            r[B * (i0 - 1) + i1] = v0[i1]
        end
    end
    return r
end
(Array{T})(x::BlockSparseVector) where {T} = Vector{T}(x)
Vector(x::BlockSparseVector{B,T}) where {B,T} = Vector{T}(x)
Array(x::BlockSparseVector) = Vector(x)

### Comparison

Base.:(==)(x::BlockSparseVector, y::BlockSparseVector) = x.elts == y.elts

### Collection

function Base.copy!(dst::BlockSparseVector{B},
                    src::BlockSparseVector{B}) where {B}
    return BlockSparseVector(copy!(dst.elts, src.elts))
end
Base.copy(x::BlockSparseVector) = BlockSparseVector(copy(x.elts))

Base.eltype(x::BlockSparseVector{B,T}) where {B,T} = T
Base.eltype(::Type{<:BlockSparseVector{B,T}}) where {B,T} = T

function Base.getindex(x::BlockSparseVector{B}, i::Integer) where {B}
    i0, i1 = fldmod1(i, B)
    return x.elts[i0][i1]
end

function Base.similar(x::BlockSparseVector; element_type)
    return BlockSparseVector(similar(x.elts; element_type))
end

Base.size(x::BlockSparseVector) = blocksize(x) .* size(x.elts)

### Sparse

function blocksparse(::Val{B}, x::SparseVector) where {B}
    return blocksparse(Val(B), rowvals(x), nonzeros(x), length(x))
end

function (SparseVector{T,I})(x::BlockSparseVector{B}) where {T,B,I}
    Is = I[]
    Vs = T[]
    for k0 in nzrange(x.elts, 1)
        i0 = rowvals(x.elts)[k0]
        v0 = nonzeros(x.elts)[k0]
        for i1 in 1:B
            if !iszero(v0[i1])
                push!(Is, B * (i0 - 1) + i1)
                push!(Vs, v0[i1])
            end
        end
    end
    return sparsevec(Is, Vs, length(x))
end
function (SparseVector{T})(x::BlockSparseVector{B,U,I}) where {T,B,U,I}
    return SparseVector{T,I}(x)
end
SparseVector(x::BlockSparseVector{B,T}) where {B,T} = SparseVector{T}(x)

function SparseArrays.droptol!(x::BlockSparseVector, tol)
    return fkeep!(x, (i, a) -> abs(a) > tol)
end

SparseArrays.dropzeros!(x::BlockSparseVector) = fkeep!(x, (i, a) -> !iszero(a))
SparseArrays.dropzeros(x::BlockSparseVector) = dropzeros!(copy(x))

function SparseArrays.fkeep!(x::BlockSparseVector{B}, f) where {B}
    fkeep!(x.elts, (i0, y) -> any(f(b * (i0 - 1) + i1, y[i1]) for i1 in 1:B))
    return x
end

SparseArrays.indtype(::BlockSparseVector{B,T}) where {B,T,I} = I
SparseArrays.indtype(::Type{<:BlockSparseVector{B,T,I}}) where {B,T,I} = I

SparseArrays.nnz(x::BlockSparseVector) = blocksize(x) * nnz(x.elts)

### Vector space

function Base.zero(x::BlockSparseVector{B,T,I}) where {B,T,I}
    return BlockSparseVector{B,T,I}(spzeros(SVector{B,T}, size(x.elts)))
end
Base.:+(x::BlockSparseVector) = x # BlockSparseVector(+x.elts)
Base.:-(x::BlockSparseVector) = BlockSparseVector(-x.elts)
function Base.:+(x::BlockSparseVector{B}, y::BlockSparseVector{B}) where {B}
    return BlockSparseVector(x.elts + y.elts)
end
function Base.:-(x::BlockSparseVector{B}, y::BlockSparseVector{B}) where {B}
    return BlockSparseVector(x.elts - y.elts)
end
Base.:*(a::Number, x::BlockSparseVector) = BlockSparseVector(a * x.elts)
Base.:\(a::Number, x::BlockSparseVector) = BlockSparseVector(a \ x.elts)
Base.:*(x::BlockSparseVector, a::Number) = BlockSparseVector(x.elts * a)
Base.:/(x::BlockSparseVector, a::Number) = BlockSparseVector(x.elts / a)

### Linear algebra
function LinearAlgebra.norm(A::BlockSparseVector{B,T}, p::Real=2) where {B,T}
    if p == 0
        return zero(norm(zero(T), 0)^0 + norm(one(T), 0)^0) + length(A)
    elseif p == 1
        r = zero(norm(zero(T), 1) + norm(one(T), 1))
        for v0 in nonzeros(A.elts)
            for v in v0
                r += norm(v, 1)
            end
        end
        return r
    elseif p == 2
        r = zero(norm(zero(T), 2)^2 + norm(one(T), 2)^2)
        for v0 in nonzeros(A.elts)
            for v in v0
                r += norm(v, 2)^2
            end
        end
        return sqrt(r)
    elseif p == Inf
        r = zero(max(norm(zero(T), Inf), norm(one(T), Inf)))
        for v0 in nonzeros(A.elts)
            for v in v0
                r = max(r, norm(v, Inf))
            end
        end
        return r
    else
        r = zero(zero(T) + norm(one(T), p)^p)
        for v0 in nonzeros(A.elts)
            for v in v0
                r += norm(v, p)^p
            end
        end
        return r^inv(p)
    end
end


text:
"""
general composite kernel
"""
const ProductCompositeKernel = NeuralKernelNetwork{Tuple{Primitive, typeof(allProduct)}}
ProductCompositeKernel(K::Vararg{AbstractKernel, N}) where {N} = ProductCompositeKernel(Primitive(K...), allProduct)

function Base.show(io::IO, kernel::ProductCompositeKernel)
	print(io, "Product(")
	join(io, kernel.layers[1].kernels, ", ")
	print(io, ")")
end

const AddCompositeKernel = NeuralKernelNetwork{Tuple{Primitive, typeof(allSum)}}
AddCompositeKernel(K::Vararg{AbstractKernel, N}) where {N} = AddCompositeKernel(Primitive(K...), allSum)

function Base.show(io::IO, kernel::AddCompositeKernel)
	print(io, "Add(")
	join(io, kernel.layers[1].kernels, ", ")
	print(io, ")")
end

"""
special composite kernels
"""
SE_mul_PeriodKernel(k1::IsoGaussKernel, k2::IsoPeriodKernel) = ProductCompositeKernel(k1, k2)
Lin_mul_LinKernel(k1::IsoLinearKernel, k2::IsoLinearKernel) = ProductCompositeKernel(k1, k2)
Period_mul_LinKernel(k1::IsoPeriodKernel, k2::IsoLinearKernel) = ProductCompositeKernel(k1, k2)
SE_mul_LinKernel(k1::IsoGaussKernel, k2::IsoLinearKernel) = ProductCompositeKernel(k1, k2)
SE_add_PeriodKernel(k1::IsoGaussKernel, k2::IsoPeriodKernel) = AddCompositeKernel(k1, k2)

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux