text:
"""
    plot(path::RegularizationPath, args...; <keyword arguments>)

Plots a `RegularizationPath` fitted with the `Lasso` package.

The minimum AICc segment is represented by a solid vertical line and the CVmin
and CV1se cross-validation selected segments by dashed vertical lines.

By default it shows nonzero coefficients at the AICc in color and the rest grayed out.

LassoPlot uses Plots.jl, so you can set any supported backend before plotting,
and add any features of the plot afterwards.

# Example:
```julia
    using Lasso, LassoPath, Plots
    path = fit(LassoPath, X, y, dist, link)
    plot(path)
```
# Arguments
- `args...` additional arguments passed to Plots.plot()

# Keywords
- `x=:segment` one of (:segment, :λ, :logλ)
- `varnames=nothing` specify variable names
- `select=:AICc` Selection criteria in (:AIC, :CVmin, :CV1se) for which coefficients
    will be shown in color. The rest are grayed out.
- `showselectors=[:AICc,:CVmin,:CV1se]` shown vertical lines
- `selectedvars=[]` Subset of the variables to present, or empty vector for all
- `kwargs...` additional keyword arguments passed along to fit(GammaLassoPath,...)
"""
function Plots.plot(path::RegularizationPath, args...;
    x=:segment, varnames=nothing, selectedvars=[], select=:AICc, showselectors=[:AICc,:CVmin,:CV1se], nCVfolds=10)
    β=coef(path)
    if hasintercept(path)
        β = β[2:end,:]
    end

    (p,nλ)=size(β)

    if varnames==nothing
        varnames=[Symbol("x$i") for i=1:p]
    end

    indata=DataFrame()
    if x==:λ
        indata[x]=path.λ
    elseif x==:logλ
        indata[x]=log.(path.λ)
    else
        x=:segment
        indata[x]=1:nλ
    end
    outdata = deepcopy(indata)

    # automatic selectors
    # xintercept = Float64[]
    dashed_vlines=Float64[]
    solid_vlines=Float64[]

    if select == :AICc || :AICc in showselectors
        minAICcix=minAICc(path)
        if select == :AICc
            push!(solid_vlines,indata[minAICcix,x])
        else
            push!(dashed_vlines,indata[minAICcix,x])
        end
    end

    if select == :CVmin || :CVmin in showselectors
        gen = Kfold(length(path.m.rr.y),nCVfolds)
        segCVmin = cross_validate_path(path;gen=gen,select=:CVmin)
        if select == :CVmin
            push!(solid_vlines,indata[segCVmin,x])
        else
            push!(dashed_vlines,indata[segCVmin,x])
        end
    end

    if select == :CV1se || :CV1se in showselectors
        gen = Kfold(length(path.m.rr.y),nCVfolds)
        segCV1se = cross_validate_path(path;gen=gen,select=:CV1se)
        if select == :CV1se
            push!(solid_vlines,indata[segCV1se,x])
        else
            push!(dashed_vlines,indata[segCV1se,x])
        end
    end

    if length(selectedvars) == 0
        if select == :all
            selectedvars = 1:p
        elseif select == :AICc
            selectedvars = findall(!iszero, β[:,minAICcix])
        elseif select == :CVmin
            selectedvars = findall(!iszero, β[:,segCVmin])
        elseif select == :CV1se
            selectedvars = findall(!iszero, β[:,segCV1se])
        else
            error("unknown selector $select")
        end
    end

    # colored paths
    for j in selectedvars
        indata[varnames[j]]=Vector(β[j,:])
    end

    # grayed out paths
    for j in setdiff(1:p,selectedvars)
        outdata[varnames[j]]=Vector(β[j,:])
    end

    inmdframe=melt(indata,x)
    outmdframe=melt(outdata,x)
    rename!(inmdframe,:value=>:coefficients)
    rename!(outmdframe,:value=>:coefficients)
    inmdframe = inmdframe[convert(BitArray,map(b->!isnan(b),inmdframe[:coefficients])),:]
    outmdframe = outmdframe[convert(BitArray,map(b->!isnan(b),outmdframe[:coefficients])),:]

    p = plot(xlabel=string(x), ylabel="Coefficient", args...)
    if size(inmdframe,1) > 0
      @df inmdframe plot!(cols(x), :coefficients, group=:variable)
    end
    if size(outmdframe,1) > 0
      @df outmdframe plot!(cols(x), :coefficients, group=:variable, palette=:grays)
    end
    if length(dashed_vlines) > 0
        vline!(dashed_vlines, line = (:dash, 0.5, 2, :black), label="")
    end
    if length(solid_vlines) > 0
        vline!(solid_vlines, line = (:solid, 0.5, 2, :black), label="")
    end

    p
end


text:
"""
    ProductOfSets{T} <: MOI.Utilities.OrderedProductOfSets{T}

The `MOI.Utilities.@product_of_sets` macro requires to know the list of sets
at compile time. In DiffOpt however, the list depends on what the user is going
to use as set as DiffOpt supports any set as long as it implements the
required function of MathOptSetDistances.
For this type, the list of sets can be given a run-time.
"""
mutable struct ProductOfSets{T} <: MOI.Utilities.OrderedProductOfSets{T}
    """
    During the copy, this counts the number of rows corresponding to
    each set. At the end of copy, `final_touch` is called, which
    converts this list into a cumulative ordering.
    """
    num_rows::Vector{Int}

    """
    A dictionary which maps the `set_index` and `offset` of a set to the
    dimension, i.e., `dimension[(set_index,offset)] → dim`.
    """
    dimension::Dict{Tuple{Int,Int},Int}

    """
    A sanity bit to check that we don't call functions out-of-order.
    """
    final_touch::Bool

    set_types::Vector{Type}
    set_types_dict::Dict{Type,Int}

    function ProductOfSets{T}() where {T}
        return new(
            Int[],
            Dict{Tuple{Int,Int},Int}(),
            false,
            Type[],
            Dict{Type,Int}(),
        )
    end
end

function MOI.Utilities.set_index(set::ProductOfSets, S::Type{<:MOI.AbstractSet})
    return get(set.set_types_dict, S, nothing)
end
MOI.Utilities.set_types(set::ProductOfSets) = set.set_types
function set_set_types(set::ProductOfSets, set_types)
    resize!(set.num_rows, length(set_types))
    fill!(set.num_rows, 0)
    resize!(set.set_types, length(set_types))
    copy!(set.set_types, set_types)
    empty!(set.set_types_dict)
    for i in eachindex(set_types)
        set.set_types_dict[set_types[i]] = i
    end
end
function add_set_types(set::ProductOfSets, S::Type)
    if !haskey(set.set_types_dict, S)
        push!(set.num_rows, 0)
        push!(set.set_types, S)
        set.set_types_dict[S] = length(set.set_types)
        return true
    end
    return false
end

ext: jl


text:
f(α::Array; p=3.) = α.^p # shape of individual cost functions

E(t, E0, γ) = E0 * (1. .+ γ).^(t .- t[1])
E(m) = E(t(m), m.economics.E0, m.economics.γ)

discount(t, ρ, tp) = .~past_mask(t, tp) .* (1. .+ ρ) .^ (- (t .- tp))
discount(m::ClimateModel) = discount(t(m), m.economics.ρ, m.domain.present_year)

T_adapt(T, A) = T .* sqrt.(1 .- A)

T_adapt(m::ClimateModel; M=false, R=false, G=false, A=false) = T_adapt(
    T(m, M=M, R=R, G=G),
    m.controls.adapt .* (1. .- .~past_mask(m) * ~A),
)

damage(β, E, Ta; discount=1.) = (β .* E .* Ta.^2) .* discount

damage(m; discounting=false, M=false, R=false, G=false, A=false) = damage(
    m.economics.β,
    E(m),
    T_adapt(m, M=M, R=R, G=G, A=A),
    discount=1. .+ discounting * (discount(m) .- 1.)
)

cost(CM, CR, CG, CA, ϵCG, E, q, M, R, G, A; discount=1., p=3.) = (
    ( ppm_to_GtCO2(q).*CM.*f(M, p=p) +
      E.*(CG.*f(G, p=p) .+ ϵCG*(G.>1.e-3)) +
      CR*f(R, p=p) +
      E.*CA.*f(A, p=p)
    ) .* discount
)
cost(m::ClimateModel; discounting=false, p=3., M=false, R=false, G=false, A=false) = cost(
    m.economics.mitigate_cost,
    m.economics.remove_cost,
    m.economics.geoeng_cost,
    m.economics.adapt_cost,
    m.economics.epsilon_cost,
    E(m),
    m.economics.baseline_emissions,
    m.controls.mitigate .* M,
    m.controls.remove .* R,
    m.controls.geoeng .* G,
    m.controls.adapt .* A,
    discount=1. .+ discounting * (discount(m) .- 1.),
    p=p
)

function cost(m::ClimateModel, controls::String; discounting=false)
    vars = Dict("M"=>false, "R"=>false, "G"=>false, "A"=>false)
    for (key, value) in vars
        if occursin(key, controls)
            vars[key] = true
        end
    end
    return cost(m, discounting=discounting; M=vars["M"], R=vars["R"], G=vars["G"], A=vars["A"])
end

benefit(damage_baseline, damage) = damage_baseline .- damage
benefit(m::ClimateModel; discounting=false, M=false, R=false, G=false, A=false) = benefit(
    damage(m, discounting=discounting),
    damage(m, discounting=discounting, M=M, R=R, G=G, A=A)
)

net_benefit(benefit, cost) = benefit .- cost
net_benefit(m::ClimateModel; discounting=true, M=false, R=false, G=false, A=false) = net_benefit(
    benefit(m, discounting=discounting, M=M, R=R, G=G, A=A),
    cost(m, discounting=discounting, M=M, R=R, G=G, A=A)
)

net_present_cost(cost, dt) = sum(cost*dt)
function net_present_cost(
        m::ClimateModel;
        discounting=true, M=false, R=false, G=false, A=false
    )
    return net_present_cost(
        cost(m, discounting=discounting, M=M, R=R, G=G, A=A),
        m.domain.dt
    )
end

net_present_benefit(net_benefit, dt) = sum(net_benefit*dt)
net_present_benefit(m::ClimateModel; discounting=true, M=false, R=false, G=false, A=false) = net_present_benefit(
    net_benefit(m, discounting=discounting, M=M, R=R, G=G, A=A),
    m.domain.dt
)


text:
@testset "test BasisMatrices.lookup" begin
    table1 = [1.0, 4.0]
    table2 = [1.0, 1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0, 4.0]

    x = [0.5, 1.0, 1.5, 4.0, 5.5]
    x2 = [0.5, 2.0]
    @test BasisMatrices.lookup(table1, x, 0) == [0, 1, 1, 2, 2]
    @test BasisMatrices.lookup(table1, x, 1) == [1, 1, 1, 2, 2]
    @test BasisMatrices.lookup(table1, x, 2) == [0, 1, 1, 1, 1]
    @test BasisMatrices.lookup(table1, x, 3) == [1, 1, 1, 1, 1]

    @test BasisMatrices.lookup(table2, x, 0) == [0, 3, 3, 12, 12]
    @test BasisMatrices.lookup(table2, x, 1) == [3, 3, 3, 12, 12]
    @test BasisMatrices.lookup(table2, x, 2) == [0, 3, 3, 8, 8]
    @test BasisMatrices.lookup(table2, x, 3) == [3, 3, 3, 8, 8]

    @test BasisMatrices.lookup([1.0], x2, 0) == [0, 1]
    @test BasisMatrices.lookup([1.0], x2, 1) == [1, 1]
    @test BasisMatrices.lookup([1.0], x2, 2) == [0, 1]
    @test BasisMatrices.lookup([1.0], x2, 3) == [1, 1]

    # test scalar version of lookup
    x2 = collect(range(-2.0, stop=4.0, length=10))
    @test [BasisMatrices.lookup(x2, -3.0, i) for i=0:3] == [0, 1, 0, 1]
    @test [BasisMatrices.lookup(x2, 5.0, i) for i=0:3] == [10, 10, 9, 9]
    @test [BasisMatrices.lookup(x2, i, 0) for i=x2] == collect(0:length(x2)-1)
    @test [BasisMatrices.lookup(x2, i, 1) for i=x2] == [1; 1:length(x2)-1]
end

@testset "test BasisMatrices._check_order" begin
    # check (::Int, ::Int) method
    @test BasisMatrices._check_order(10, 0) == fill(0, 1, 10)
    @test BasisMatrices._check_order(1, 0) == fill(0, 1, 1)

    # check (::Int, ::Vector) method
    ov = [0, 0]
    @test BasisMatrices._check_order(2, ov) == reshape(ov, 1, 2)
    @test_throws DimensionMismatch BasisMatrices._check_order(3, ov)

    # check (::Int, ::Matrix) method
    om = [0 0]
    @test BasisMatrices._check_order(2, om) == om
    @test BasisMatrices._check_order(1, om) == om'
    @test_throws DimensionMismatch BasisMatrices._check_order(3, ov)
end

@testset "test BasisMatrices.ckronx" begin
    # will test by constructing an interpoland, then evaluating at the nodes
    # and verifying that we get back our original function
    basis = Basis(Basis(Spline(), 13, -1.0, 1.0, 3),
                  Basis(Spline(), 18, -5.0, 3.0, 3))
    X, x12 = nodes(basis);

    # make up a funciton and evaluate at the nodes
    f(x1, x2) = cos.(x1) ./ exp.(x2)
    f(X::Matrix) = f(X[:, 1], X[:, 2])
    y = f(X)

    # fit the interpoland in Tensor form (tensor b/c using x12)
    c, bs = funfitxy(basis, x12, y);

    # verify that we are actually interpolating -- all heavy lifting in funeval
    # is done by ckronx so this is effectively testing that we wrote that
    # function properly
    @test maximum(abs, funeval(c, bs, [0 0]) - y) <= 1e-13
end

@testset "test row_kron" begin
    h = ["a" "b"; "c" "d"]
    z = ["1" "2" "3"; "4" "5" "6"]
    want = ["a1" "a2" "a3" "b1" "b2" "b3"; "c4" "c5" "c6" "d4" "d5" "d6"]
    @test row_kron(h, z) == want

    # now test on some bigger matrices
    a = randn(400, 3)
    b = randn(400, 5)
    out = row_kron(a, b)
    @test size(out) == (400, 15)

    rows_good = true
    for row=1:400
        rows_good &= out[row, :] == kron(a[row, :], b[row, :])
    end
    @test rows_good == true

    for i in 1:100
        A = SparseArrays.sprandn(30, 5, 0.3)
        B = SparseArrays.sprandn(30, 3, 0.3)

        want = row_kron(A, B)

        @test maximum(abs, want - row_kron(Array(A), B)) == 0
        @test maximum(abs, want - row_kron(A, Array(B))) == 0
        @test maximum(abs, want - row_kron(Array(A), Array(B))) == 0
    end
end

@testset "RowKron" begin
    # throws when try ot pass a non matrix
    @test_throws MethodError RowKron((rand(2, 2), "foo"))

    rk = RowKron(Matrix(1.0I, 3, 3), Matrix(1.0I, 3, 3), Matrix(1.0I, 3, 100))
    @test size(rk, 1) == 3
    @test size(rk, 2) == 900
    @test size(rk) == (3, 900)
    @test BasisMatrices.sizes(rk, 1) == [3, 3, 3]
    @test BasisMatrices.sizes(rk, 2) == [3, 3, 100]
    for i in 3:10
        @test BasisMatrices.sizes(rk, i) == [1, 1, 1]
    end

    bs = AbstractMatrix[SparseArrays.sprandn(5, rand(5:13), 0.8) for _ in 1:3]
    rk_last_sparse = RowKron(bs...)
    big_last_sparse = reduce(row_kron, bs)

    push!(bs, Matrix(1.0I, 5, 5))
    rk_last_full = RowKron(bs...)
    big_last_full = reduce(row_kron, bs)

    for (rk, _big) in [(rk_last_sparse, big_last_sparse), (rk_last_full, big_last_full)]
        c = rand(size(rk, 2))
        c2 = rand(size(rk, 1))

        # non-mutating
        @test rk * c == _big * c
        @test rk * [c c] == _big * [c c]

        # mutating
        mul!(c2, rk, c)
        @test c2 == _big*c

        # mutating matrix
        c2mat = [c2 c2]
        mul!(c2mat, rk, [c c])
        @test c2mat == _big * [c c]

        ## transpose op
        # reuse c and c2 for this...
        mul!(c, Transpose(rk), c2)
        @test c == _big'*c2

        # vector and matrix versions
        @test *(Transpose(rk), c2) == c
        @test *(Transpose(rk), [c2 c2]) == [c c]

        # mutating matrix
        c1mat = [c c]
        mul!(c1mat, Transpose(rk), [c2 c2])
        @test c1mat == _big'*[c2 c2]

    end

end

@testset "cdprodx" begin
    for nrow in 3:3:100
        for nb in 2:5
            b = [SparseArrays.sprand(nrow, rand(5:13), 0.3) for _ in 1:nb]
            c = rand(prod([size(A, 2) for A in b]))

            full_b = reduce(row_kron, b)
            want = full_b * c
            have = BasisMatrices.cdprodx(b, c)

            @test maximum(abs, want - have) < 1e-12

            # test RowKron object
            rk = RowKron(b...)
            @test maximum(abs, want - rk*c) < 1e-12

            # now test transpose
            c2 = rand(size(rk, 1))
            want2 = full_b'c2
            @test maximum(abs, want2 - *(Transpose(rk), c2)) < 1e-12
        end
    end

    x = rand(10)
    @test BasisMatrices.cdprodx(Matrix(1.0I, 10, 10), x) == x
end

@testset "nodeunif" begin
    X, x12 = BasisMatrices.nodeunif([5, 5], [0, 1], [3, 3])
    @test x12[1] == range(0, stop=3, length=5)
    @test x12[2] == range(1, stop=3, length=5)
    @test size(X) == (25, 2)
    @test X[:, 1] == repeat(x12[1], 5)
    @test X[:, 2] == repeat(x12[2], inner=[5])

    x, x1 = BasisMatrices.nodeunif(5, 0, 3)
    @test x == range(0, stop=3, length=5)
    @test x1 == range(0, stop=3, length=5)
end


text:
int main () {
int i = 0 ;
while (i < 10) {
if (i % 2 == 0) printInt(i) ;
i++ ;
}
return 0 ;
}

max_stars_repo_name: fredefox/javalette

max_issues_repo_name: fredefox/javalette

max_forks_repo_name: fredefox/javalette


text:
struct ProbIntsCache{T}
  σ::T
  order::Int
end
function (p::ProbIntsCache)(integrator)
  integrator.u .= integrator.u .+ p.σ*sqrt(integrator.dt^(2*p.order))*randn(size(integrator.u))
end

"""
Conrad P., Girolami M., Särkkä S., Stuart A., Zygalakis. K, Probability
Measures for Numerical Solutions of Differential Equations, arXiv:1506.04592
"""
function ProbIntsUncertainty(σ,order,save=true)
  affect! = ProbIntsCache(σ,order)
  condtion = (t,u,integrator) -> true
  save_positions = (save,false)
  DiscreteCallback(condtion,affect!,save_positions=save_positions)
end

struct AdaptiveProbIntsCache
  order::Int
end
function (p::AdaptiveProbIntsCache)(integrator)
  integrator.u .= integrator.u .+ integrator.EEst*sqrt(integrator.dt^(2*p.order))*randn(size(integrator.u))
end

function AdaptiveProbIntsUncertainty(order,save=true)
  affect! = AdaptiveProbIntsCache(order)
  condtion = (t,u,integrator) -> true
  save_positions = (save,false)
  DiscreteCallback(condtion,affect!,save_positions=save_positions)
end


text:
function load_mets_map()
    ## ---------------------------------------------------------------
    # This create a map between the metabolites ids of the Human1 model
    # and the ids used in Rath 2017 (https://pure.mpg.de/pubman/item/item_2508673_4)

    mets_map = Dict{AbstractString, AbstractString}()
    # Mammalian cells only contain L-LDH so that in humans the lactate produced is almost exclusively L-lactate.
    # https://acutecaretesting.org/en/articles/l-lactate-and-d-lactate-clinical-significance-of-the-difference
    mets_map["LAC"] = "m02403s" # "L-lactate[s]" # "m01716s" # "D-lactate[s]"
    mets_map["NH4"] = "m02579s"
    mets_map["HCO3"] = "m02046s"
    mets_map["ASN"] = "m01369s"
    mets_map["MET"] = "m02471s"
    mets_map["GLC"] = "m01965s"
    mets_map["GLY"] = "m01986s"
    mets_map["TYR"] = "m03101s"
    mets_map["LEU"] = "m02360s"
    mets_map["PHE"] = "m02724s"
    mets_map["GLU"] = "m01974s"
    mets_map["ASP"] = "m01369s"
    mets_map["GLN"] = "m01975s"
    mets_map["CYS"] = "m01628s"
    mets_map["ILE"] = "m02184s"
    mets_map["PYR"] = "m02819s"
    mets_map["H"] = "m02039s"
    mets_map["LYS"] = "m02426s"
    mets_map["TRP"] = "m03089s"
    mets_map["SER"] = "m02896s"
    mets_map["HIS"] = "m02125s"
    mets_map["VAL"] = "m03135s"
    mets_map["ALA"] = "m01307s"
    mets_map["ARG"] = "m01365s"
    mets_map["PRO"] = "m02770s"
    mets_map["GAL"] = "m01910s"
    mets_map["THR"] = "m02993s"

    # making it both ways
    for (k, v) in mets_map
        mets_map[v] = k
    end

    return mets_map
end


text:
#=
NOTE: a greedy match takes

- head_chars::Vector{Vector{Char}}
- tail_chars::Vector{Char}
- check::Regex

The head_chars indicate the accepted chars in order after the first triggering
char (so the first vector is for the first character after trigger, second
for second etc), as soon as we're beyond the range of head_chars, we check
whether the characters are in `tail_chars`

Finally the check regex allows to discard expressions that match the head
and char but might not overall meet the required regex. This is only
required for things where head/tail are insufficient.
The regex includes the trigger character.
=#

"""
    F_LINE_RETURN

Finder for a line return (`\n`) followed by any number of whitespaces or tabs.
These will subsequently be checked to see if they are followed by something
that constitutes a list item or not.
"""
const F_LINE_RETURN = greedy_match(
    tail_chars=[' ', '\t']
)

"""
    F_DIV_OPEN

Finder for `@@div` checking that `div` matches a simplified rule for allowed
CSS class names. The complete rule being `-?[_a-zA-Z]+[_a-zA-Z0-9-]*` which
we simplify here to `[a-zA-Z]+[_a-zA-Z0-9-]` and allow `,` for separation so
`@@d1,d2` is allowed and corresponds to a setting where we pass two classes
`class="d1 d2"`.
"""
const F_DIV_OPEN = greedy_match(
    head_chars=[['@'], ALPHA_LATIN],
    tail_chars=vcat(ALPHANUM_LATIN, ['-', '_', ','])
)

"""
    F_LX_COMMAND

Finder for latex command. First character is `[a-zA-Z]`.
We do allow numbers (there's no ambiguity because `\\com1` is not allowed to
mean `\\com{1}` unlike in LaTeX).
Underscores are allowed *inside* the command but not at the very start or very
end to avoid confusion respectively with the escaped `_` character or the
emphasis in markdown, `*` are not allowed anywhere (including at the end).
See also the check pattern.
"""
const F_LX_COMMAND = greedy_match(
    head_chars=[vcat(ALPHA_LATIN, ['_'])],
    tail_chars=vcat(ALPHANUM_LATIN, ['_']),
    check=LX_COMMAND_PAT
)

"""
    F_LANG_*

Finder for code blocks. I.e. something like a sequence of 3, 4 or 5 backticks
followed by a valid combination of letter defining a language.
"""
const F_LANG_3 = greedy_match(
    head_chars=[['`'], ['`']],
    tail_chars=vcat(ALPHANUM_LATIN, ['-']),
    check=CODE_LANG3_PAT
)
const F_LANG_4 = greedy_match(
    head_chars=[['`'], ['`'], ['`']],
    tail_chars=vcat(ALPHANUM_LATIN, ['-']),
    check=CODE_LANG4_PAT
)
const F_LANG_5 = greedy_match(
    head_chars=[['`'], ['`'], ['`'], ['`']],
    tail_chars=vcat(ALPHANUM_LATIN, ['-']),
    check=CODE_LANG5_PAT
)

"""
    F_HTML_ENTITY

Finder for html entities.
"""
const F_HTML_ENTITY = greedy_match(
    tail_chars=vcat(ALPHANUM_LATIN, ['#', ';']),
    check=HTML_ENTITY_PAT
)

"""
    F_EMOJI

Finder for emojis (those will have to be validated separately to check Julia
recognises them).
"""
const F_EMOJI = greedy_match(
    tail_chars=vcat(ALPHANUM_LATIN, [':', '+', '-', '_']),
    check=EMOJI_PAT
)

"""
    F_FOOTNOTE

Finder for footnotes.
"""
const F_FOOTNOTE = greedy_match(
    head_chars=[['^'], ALPHANUM_ALL],
    tail_chars=vcat(ALPHANUM_ALL, [']', ':', '_']),
    check=FOOTNOTE_PAT
)


text:
using ObaraSaika
using Test

@testset "ObaraSaika.jl" begin
    @test ObaraSaika.get_ijk_list(0) == Array{Int, 1}[[0,0,0]]
    @test ObaraSaika.get_ijk_list(1) == Array{Int, 1}[[1,0,0],[0,1,0],[0,0,1]]
    @test ObaraSaika.get_ijk_list(2) == Array{Int, 1}[[2,0,0],[1,1,0],[1,0,1],[0,2,0],[0,1,1],[0,0,2]]

    @test ObaraSaika.get_shell2(0, 0) == Array{Int, 1}[
        [0, 0, 0, 0, 0, 0]
    ]
    @test ObaraSaika.get_shell2(0, 1) == Array{Int, 1}[
        [0, 0, 0, 1, 0, 0],
        [0, 0, 0, 0, 1, 0],
        [0, 0, 0, 0, 0, 1]
    ]
    @test ObaraSaika.get_shell2(1, 0) == Array{Int, 1}[
        [1, 0, 0, 0, 0, 0],
        [0, 1, 0, 0, 0, 0],
        [0, 0, 1, 0, 0, 0]
    ]
    @test ObaraSaika.get_shell2(0, 2) == Array{Int, 1}[
        [0, 0, 0, 2, 0, 0],
        [0, 0, 0, 1, 1, 0],
        [0, 0, 0, 1, 0, 1],
        [0, 0, 0, 0, 2, 0],
        [0, 0, 0, 0, 1, 1],
        [0, 0, 0, 0, 0, 2]
    ]
end


text:
using Nettle

msg = "Rosetta Code"

h = HashState(MD4)
update!(h, msg)
h = hexdigest!(h)

println("\"", msg, "\" => ", h)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function index_page(app::DashApp; debug = false)
    metas = ["""<meta http-equiv="X-UA-Compatible" content="IE=edge">""", """<meta charset="UTF-8">"""]
    title = app.name
    css = join(map(app.external_stylesheets) do s
        """<link rel="stylesheet" href="$(s)">"""
    end, "\n"
    )
    external_scripts = join(map(app.external_scripts) do s
        """"""
    end, "\n")

    app_entry = """

                Loading...

    """
    config = (
        url_base_pathname = "$(app.url_base_pathname)",
        requests_pathname_prefix = "$(app.url_base_pathname)",
        ui = true,
        props_check = debug,
        show_undo_redo = false
    )

    scripts = ComponentPackages.components_js_include(app.url_base_pathname, debug = debug)

    """<!DOCTYPE html>

            $(join(metas, ""))
            <title>$(title)</title>
            $(css)

            $(app_entry)

            $external_scripts
            $scripts

    """
end

function dependencies_json(app::DashApp)
    id_prop_named(p::IdProp) = (id = p[1], property = p[2])
    result = map(values(app.callbacks)) do dep
        (inputs = id_prop_named.(dep.id.input),
        state = id_prop_named.(dep.id.state),
        output = output_string(dep.id)
        )
    end
    return JSON2.write(result)
end

function process_callback(app::DashApp, body::String)
    params = JSON2.read(body)
    output = Symbol(params[:output])
    if !haskey(app.callbacks, output)
        return []
    end
    convert_values(inputs) = map(inputs) do x
        return get(x, :value, nothing)
    end
    args = []
    if app.callbacks[output].pass_changed_props
        push!(args, params[:changedPropIds])
    end
    if haskey(params, :state)
        append!(args, convert_values(params.state))
    end
    if haskey(params, :inputs)
        append!(args, convert_values(params.inputs))
    end

    res = app.callbacks[output].func(args...)
    if length(app.callbacks[output].id.output) == 1
        if !(res isa NoUpdate)
            return Dict(
                :response => Dict(
                    :props => Dict(
                        Symbol(app.callbacks[output].id.output[1][2]) => Front.to_dash(res)
                    )
                )
            )
        else
            throw(PreventUpdate())
        end
    end
    response = Dict{Symbol, Any}()
    for (ind, out) in enumerate(app.callbacks[output].id.output)
        if !(res[ind] isa NoUpdate)
            push!(response,
            Symbol(out[1]) => Dict(
                Symbol(out[2]) => Front.to_dash(res[ind])
            )
            )
        end
    end
    if length(response) == 0
        throw(PreventUpdate())
    end
    return Dict(:response=>response, :multi=>true)

end

function process_assets(app::DashApp, path)
    assets_path = "$(app.url_base_pathname)assets/"
    filename = joinpath(app.assets_folder, replace(path, assets_path=>""))
    try
        return HTTP.Response(200, [], body = read(filename))
    catch
        return HTTP.Response(404)
    end
end

function make_handler(app::DashApp; debug::Bool = false)
    function (req::HTTP.Request)
        uri = HTTP.URI(req.target)
        ComponentPackages.@register_js_sources(uri.path, app.url_base_pathname)
        if uri.path == "$(app.url_base_pathname)"
            return HTTP.Response(200, index_page(app, debug = debug))
        end
        if uri.path == "$(app.url_base_pathname)_dash-layout"
            return HTTP.Response(200, ["Content-Type" => "application/json"], body = JSON2.write(app.layout))
        end
        if uri.path == "$(app.url_base_pathname)_dash-dependencies"
            return HTTP.Response(200, ["Content-Type" => "application/json"], body = dependencies_json(app))
        end
        if startswith(uri.path, "$(app.url_base_pathname)assets/")
            return process_assets(app, uri.path)
        end
        if uri.path == "$(app.url_base_pathname)_dash-update-component" && req.method == "POST"
            try
                return HTTP.Response(200, ["Content-Type" => "application/json"],
                    body = JSON2.write(
                        process_callback(app, String(req.body))
                    )
                )
            catch e
                if isa(e,PreventUpdate)
                    return HTTP.Response(204)
                else
                    rethrow(e)
                end
            end
        end
        return HTTP.Response(404)
    end
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

debug = false
using Base.Test

using Base.LinAlg: BlasComplex, BlasFloat, BlasReal, QRPivoted

n = 10

# Split n into 2 parts for tests needing two matrices
n1 = div(n, 2)
n2 = 2*n1

srand(1234321)

areal = randn(n,n)/2
aimg  = randn(n,n)/2
a2real = randn(n,n)/2
a2img  = randn(n,n)/2
breal = randn(n,2)/2
bimg  = randn(n,2)/2

for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
    a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
    a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    asym = a'+a                  # symmetric indefinite
    apd  = a'*a                 # symmetric positive-definite
    ε = εa = eps(abs(float(one(eltya))))

    for eltyb in (Float32, Float64, Complex64, Complex128, Int)
        b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex.(breal, bimg) : breal)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

        α = rand(eltyb)
        aα = fill(α,1,1)
        @test qrfact(α)[:Q]*qrfact(α)[:R] ≈ qrfact(aα)[:Q]*qrfact(aα)[:R]
        @test abs(qrfact(α)[:Q][1,1]) ≈ one(eltyb)
        tab = promote_type(eltya,eltyb)

debug && println("\ntype of a: ", eltya, " type of b: ", eltyb, "\n")
debug && println("QR decomposition (without pivoting)")
        for i = 1:2
            let a = i == 1 ? a : view(a, 1:n - 1, 1:n - 1), b = i == 1 ? b : view(b, 1:n - 1), n = i == 1 ? n : n - 1
                qra   = @inferred qrfact(a)
                @inferred qr(a)
                q, r  = qra[:Q], qra[:R]
                @test_throws KeyError qra[:Z]
                @test q'*full(q, thin=false) ≈ eye(n)
                @test q*full(q, thin=false)' ≈ eye(n)
                @test q'*eye(n)' ≈ full(q, thin=false)'
                @test full(q, thin=false)'q ≈ eye(n)
                @test eye(n)'q' ≈ full(q, thin=false)'
                @test q*r ≈ a
                @test a*(qra\b) ≈ b atol=3000ε
                @test full(qra) ≈ a
                @test A_mul_Bc(eye(eltyb,size(q.factors,2)),q)*full(q,thin=false) ≈ eye(n) atol=5000ε
                if eltya != Int
                    @test eye(eltyb,n)*q ≈ convert(AbstractMatrix{tab},q)
                    ac = copy(a)
                    @test qrfact!(a[:, 1:5])\b == qrfact!(view(ac, :, 1:5))\b
                end

debug && println("Thin QR decomposition (without pivoting)")
                qra   = @inferred qrfact(a[:,1:n1], Val{false})
                @inferred qr(a[:,1:n1], Val{false})
                q,r   = qra[:Q], qra[:R]
                @test_throws KeyError qra[:Z]
                @test q'*full(q, thin=false) ≈ eye(n)
                @test q'*full(q) ≈ eye(n,n1)
                @test q*r ≈ a[:,1:n1]
                @test q*b[1:n1] ≈ full(q)*b[1:n1] atol=100ε
                @test q*b ≈ full(q,thin=false)*b atol=100ε
                @test_throws DimensionMismatch q*b[1:n1 + 1]
                @test_throws DimensionMismatch b[1:n1 + 1]*q'
                @test A_mul_Bc(UpperTriangular(eye(eltyb,size(q.factors,2))),q)*full(q,thin=false) ≈ eye(n1,n) atol=5000ε
                if eltya != Int
                    @test eye(eltyb,n)*q ≈ convert(AbstractMatrix{tab},q)
                end

debug && println("(Automatic) Fat (pivoted) QR decomposition")
                @inferred qrfact(a, Val{true})
                @inferred qr(a, Val{true})

                qrpa  = factorize(a[1:n1,:])
                q,r = qrpa[:Q], qrpa[:R]
                @test_throws KeyError qrpa[:Z]
                p = qrpa[:p]
                @test q'*full(q, thin=false) ≈ eye(n1)
                @test q*full(q, thin=false)' ≈ eye(n1)
                @test (UpperTriangular(eye(eltya,size(q,2)))*q')*full(q, thin=false) ≈ eye(n1)
                @test q*r ≈ (isa(qrpa,QRPivoted) ? a[1:n1,p] : a[1:n1,:])
                @test q*r[:,invperm(p)] ≈ a[1:n1,:]
                @test q*r*qrpa[:P].' ≈ a[1:n1,:]
                @test a[1:n1,:]*(qrpa\b[1:n1]) ≈ b[1:n1] atol=5000ε
                @test full(qrpa) ≈ a[1:5,:]
                @test_throws DimensionMismatch q*b[1:n1+1]
                @test_throws DimensionMismatch b[1:n1+1]*q'
                if eltya != Int
                    @test eye(eltyb,n1)*q ≈ convert(AbstractMatrix{tab},q)
                end

debug && println("(Automatic) Thin (pivoted) QR decomposition")
                qrpa  = factorize(a[:,1:n1])
                q,r = qrpa[:Q], qrpa[:R]
                @test_throws KeyError qrpa[:Z]
                p = qrpa[:p]
                @test q'*full(q, thin=false) ≈ eye(n)
                @test q*full(q, thin=false)' ≈ eye(n)
                @test q*r ≈ a[:,p]
                @test q*r[:,invperm(p)] ≈ a[:,1:n1]
                @test full(qrpa) ≈ a[:,1:5]
                @test_throws DimensionMismatch q*b[1:n1+1]
                @test_throws DimensionMismatch b[1:n1+1]*q'
                @test A_mul_Bc(UpperTriangular(eye(eltyb,size(q.factors,2))),q)*full(q,thin=false) ≈ eye(n1,n) atol=5000ε
                if eltya != Int
                    @test eye(eltyb,n)*q ≈ convert(AbstractMatrix{tab},q)
                end
            end
        end

debug && println("Matmul with QR factorizations")
        if eltya != Int
            qrpa = factorize(a[:,1:n1])
            q, r = qrpa[:Q], qrpa[:R]
            @test A_mul_B!(full(q, thin=false)',q) ≈ eye(n)
            @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q)
            @test A_mul_Bc!(full(q, thin=false),q) ≈ eye(n)
            @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q)
            @test_throws BoundsError size(q,-1)
            @test_throws DimensionMismatch Base.LinAlg.A_mul_B!(q,zeros(eltya,n1+1))
            @test_throws DimensionMismatch Base.LinAlg.Ac_mul_B!(q,zeros(eltya,n1+1))

            qra = qrfact(a[:,1:n1], Val{false})
            q, r = qra[:Q], qra[:R]
            @test A_mul_B!(full(q, thin=false)',q) ≈ eye(n)
            @test_throws DimensionMismatch A_mul_B!(eye(eltya,n+1),q)
            @test A_mul_Bc!(full(q, thin=false),q) ≈ eye(n)
            @test_throws DimensionMismatch A_mul_Bc!(eye(eltya,n+1),q)
            @test_throws BoundsError size(q,-1)
            @test_throws DimensionMismatch q * eye(Int8,n+4)
        end
    end
end

# Because transpose(x) == x
@test_throws ErrorException transpose(qrfact(randn(3,3)))
@test_throws ErrorException ctranspose(qrfact(randn(3,3)))
@test_throws ErrorException transpose(qrfact(randn(3,3), Val{false}))
@test_throws ErrorException ctranspose(qrfact(randn(3,3), Val{false}))
@test_throws ErrorException transpose(qrfact(big.(randn(3,3))))
@test_throws ErrorException ctranspose(qrfact(big.(randn(3,3))))

# Issue 7304
let
    A = [-√.5 -√.5; -√.5 √.5]
    Q = full(qrfact(A)[:Q])
    @test vecnorm(A-Q) < eps()
end

let
    debug && println("qr on AbstractVector")

    vr = [3.0, 4.0]
    for Tr in (Float32, Float64)
        for T in (Tr, Complex{Tr})
            v = convert(Vector{T}, vr)
            nv, nm = qr(v)
            @test norm(nv - [0.6, 0.8], Inf) < eps(Tr)
            @test nm == 5.0
        end
    end
end

@test qr(Int[]) == (Int[],1)
@test Base.LinAlg.qr!(Int[1]) == (Int[1],1)

B = rand(7,2)
@test (1:7)\B ≈ collect(1:7)\B

# Issue 16520
@test_throws DimensionMismatch ones(3,2)\(1:5)


text:
"""
    Position(llh, xzy_earth, xyz_ecef)

Position information:
    - llh (latitude, longitude, height) [rad, rad, m]
    - Earth (x, y, z) [m]
    - ECEF (x, y, z) [m]

# Constructors
    PositionLLH(lat, lon, h, xe=0., ye=0., ze=0.)
    PositionECEF(x, y, z, xe=0., ye=0., ze=0.)
    PositionEarth(x, y, z, lat=0., lon=0., h=0.)
"""
struct Position
    # latitude [rad], longitude [rad], height [m]
    llh::Array{T, 1} where T<:Number
    # x, y, z [m] Earth
    xyz_earth::Array{T, 1} where T<:Number
    # x, y, z [m] ECEF
    xyz_ecef::Array{T, 1} where T<:Number
end

# Constructors
PositionLLH(lat, lon, h, xe=0., ye=0.) = Position([lat, lon, h], [xe, ye, -h], llh2ecef(lat, lon, h))
PositionECEF(x, y, z, xe=0., ye=0.) = Position(ecef2llh(x, y, z), [xe, ye, -ecef2llh(x, y, z)[3]], [x, y, z])
PositionEarth(x, y, z, lat=0., lon=0.) = Position([lat, lon, -z], [x, y, z], llh2ecef(lat, lon, -z))
Position() = Position([0, 0, 0], [0, 0, 0], llh2ecef(0, 0, 0))

get_llh(pos::Position) = pos.llh
get_xyz_earth(pos::Position) = pos.xyz_earth
get_xyz_ecef(pos::Position) = pos.xyz_ecef
get_height(pos::Position) = pos.llh[3]


text:
export msec_to_sample
export sec_to_sample
export sample_to_sec
export sample_to_msec
export sample_to_hr

msec_to_sample(msec:: Real, fs::Real, to_int::Bool = true) =
    to_int ? convert(Integer, msec*fs/1000) : msec*fs/1000
sec_to_sample(sec:: Real, fs::Real, to_int::Bool = true) =
    to_int ? convert(Integer, sec*fs) : sec*fs
sample_to_sec(sample:: Real, fs::Real) = sample/fs
sample_to_msec(sample:: Real, fs::Real) = sample/fs * 1000
hr_to_sample(hr::Real, fs::Real, to_int::Bool = true) = sec_to_sample(60/hr, fs, to_int)
sample_to_hr(sample::Real, fs::Real) = 60/sample_to_sec(sample, fs)

msec_to_sample(500.0, 500)


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

struct Variable
    name::String
    defaultValue::Float64
    range::Tuple{Float64, Float64}
    units::String # For now it's a string
    description::String

    function Variable(name, defaultValue, range, units, description)
        if match(r"^[a-zA-Z_][a-zA-Z0-9_]*$", name) === nothing
            error("Name must be a valid identifier: $name");
        end
        if range[1] > range[2]
            error("Invalid range in $name: $(range[1]) > $(range[2])");
        end
        new(name, defaultValue, range, units, description)
    end
end

Base.show(io::IO, v::Variable) = print(io, "Var $(v.name) in $(v.range) default: $(v.defaultValue) units: $(v.units) $(v.description)");

struct Variables
    values::Vector{Any}
    variables::Array{Variable,1}
    nameToIndex::Dict{String,Int32}
end

function Variables(variables::Array{Variable,1})
    values = Vector{Any}([variable.defaultValue for variable = variables])
    nameToIndex = Dict{String,Int32}(variables[i].name => i for i = 1:length(variables));
    if length(nameToIndex) != length(variables)
        error("Duplicate variable names");
    end
    return Variables(values, copy(variables), nameToIndex);
end

function Variables(variables::Variables, values::Vector{Any})
    if length(variables.values) != length(values)
        error("Cannot create $(length(variables.values)) variables with a vector of length $(length(values))");
    end
    return Variables(values, variables.variables, variables.nameToIndex);
end

function Variables(variables::Variables, values::Vector)
    return Variables(variables, Vector{Any}(values));
end

function variablesToMatlab(variables::Variables, filename::String)
    v = Dict{String,Any}(variables.variables[i].name => variables.values[i] for i = 1:length(variables.values));
    file = matopen(filename, "w");
    write(file, "v", v);
    close(file);
end

function Base.getproperty(variables::Variables, symbol::Symbol)
    if symbol === :values || symbol === :nameToIndex || symbol === :variables
        return getfield(variables, symbol)
    end
    variables.values[variables.nameToIndex[String(symbol)]]
end

function Base.getindex(variables::Variables, name::String)
    variables.values[variables.nameToIndex[name]]
end

function Base.setproperty!(variables::Variables, symbol::Symbol, value)
    if symbol === :values
        variables.values[1:end] = value;
        return
    end
    variables.values[variables.nameToIndex[String(symbol)]] = value;
end

function Base.setindex!(variables::Variables, value, name::String)
    variables.values[variables.nameToIndex[name]] = value
end

function Base.show(io::IO, ::MIME"text/plain", vars::Variables)
    println("Variables:")
    maxLen = maximum(v -> length(v.name), vars.variables);
    for i = 1:length(vars.variables)
        v = vars.variables[i];
        value = vars.values[i];
        println(io, "$(lpad(v.name, maxLen + 2)) = $value $(v.units) in $(v.range) $(v.description)");
    end
end

# Combining variables
function variablesUnion(v1::Variables, v2::Variables)::Variables
    newSize = length(v1.variables);
    mergedNameToIndex = copy(v1.nameToIndex);
    for var in v2.variables
        if !haskey(mergedNameToIndex, var.name)
            mergedNameToIndex[var.name] = newSize += 1;
        end
    end
    newVariables = Array{Variable}(undef, newSize);
    newValues = Vector{Any}(undef, newSize);
    newVariables[1:length(v1.variables)] = v1.variables;
    newValues[1:length(v1.variables)] = v1.values;
    for (name, index) in v2.nameToIndex
        newIndex = mergedNameToIndex[name];
        newValues[newIndex] = v2.values[index];
        newVariables[newIndex] = v2.variables[index];
    end
    return Variables(newValues, newVariables, mergedNameToIndex);
end

function variablesSubtract(v::Variables, toRemove::AbstractSet{String})::Variables
    newValues = Vector{Any}();
    newVariables = Array{Variable,1}();
    newNameToIndex = Dict{String,Int32}();
    for var in v.variables
        if in(var.name, toRemove)
            continue;
        end
        index = v.nameToIndex[var.name];
        push!(newVariables, var);
        push!(newValues, v.values[index]);
        newNameToIndex[var.name] = length(newValues);
    end
    return Variables(newValues, newVariables, newNameToIndex);
end

function renameVariables(variables::Variables, renameFunction::Function)::Variables
    newVariables = [Variable(renameFunction(v.name), v.defaultValue, v.range, v.units, v.description) for v in variables.variables];
    result = Variables(newVariables);
    result.values = variables.values;
    return result;
end

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
module PointArrays

export PointArray

const PointArray = Array{UInt32, 2}

const ZERO_UINT32 = zero(UInt32)
const ONE_UINT32  = one(UInt32)
const DEFAULT_OFFSET = (ZERO_UINT32, ZERO_UINT32, ZERO_UINT32)

const MAX_BOUNDARY_DISTANCE = 100000

"""
parameter:
    bin_im: binary array. object voxels are false, non-object voxels are true!
"""
function from_binary_image(bin_im::Union{Array{Bool,3}, BitArray{3}})
    N = length(bin_im) - count(bin_im)
    ret = zeros(UInt32, (N,3))
    i = 0
    for z in ONE_UINT32:UInt32(size(bin_im, 3))
        for y in ONE_UINT32:UInt32(size(bin_im, 2))
            for x in ONE_UINT32:UInt32(size(bin_im, 1))
                if !bin_im[x,y,z]
                    i += 1
                    ret[i,:] = [x,y,z]
                end
            end
        end
    end
    ret
end

"""
find points inside an object from a segmentation array.
"""
function from_seg(seg::Array{T,3}; obj_id::T=convert(T,1),
                     offset::NTuple{3,UInt32}=DEFAULT_OFFSET) where T
    # compute the number of object voxels
    nov = 0
    for i in eachindex(seg)
        if seg[i] == obj_id
            nov += 1
        end
    end
    # initialize the points array with offset
    points = Array{UInt32}(undef, nov, 3)
    points[:,1] .= offset[1]
    points[:,2] .= offset[2]
    points[:,3] .= offset[3]

    # position of current row
    pos = 0
    for z in ONE_UINT32 : convert(UInt32, size(seg, 3))
        for y in ONE_UINT32 : convert(UInt32, size(seg, 2))
            for x in ONE_UINT32 : convert(UInt32, size(seg, 1))
                if seg[x,y,z] == obj_id
                    pos += 1
                    points[pos,:] .+= [x,y,z]
                end
            end
        end
    end
    return points
end

"""
find out the boundary voxels and represent them as indexes in the point array
"""
function get_boundary_point_indexes(self::Array{T,2}, seg::Array{TSeg,3},
                                    obj_id::TSeg = TSeg(1)) where {T, TSeg}
    # compute the indexes of boundary voxels
    ret = Vector{T}()
    for i in 1:size(self,1)
        point = self[i,:]
        x,y,z = (point...,)
        if  z==1 || z==size(seg, 3) ||
            y==1 || y==size(seg, 2) ||
            x==1 || x==size(seg, 1)
            push!(ret, i)
        else
            if  seg[x-1,y,z]!=obj_id || seg[x+1,y,z]!=obj_id ||
                seg[x,y-1,z]!=obj_id || seg[x,y+1,z]!=obj_id ||
                seg[x,y,z-1]!=obj_id || seg[x,y,z+1]!=obj_id
                push!(ret, i)
            end
        end
    end
    return ret
end

"""
add offset to points
"""
function add_offset!(self::Array{T,2}, offset::NTuple{3,T}) where T
    self[:, 1] .+= offset[1]
    self[:, 2] .+= offset[2]
    self[:, 3] .+= offset[3]
end

"""
    merge(self::Array{T,3}, other::Array{T,3})
"""
function merge(self::Array{T,2}, other::Array{T,2}) where T
    vcat(self, other)
end

end # module


text:
@testset "fastio" begin
  data = [1/3*ones(500); pi*ones(500); nextfloat(0.0)*ones(500)]
  buf = IOBuffer()
  for i in data
    print(buf, i)
    print(buf, "  ")
  end
  seek(buf, 0) # reset buffer
  a = DGFun.tseq2(buf, zeros(0))
  @test a == data

  data = convert(Array{UInt8}, "1.0"*("0"^2048)*"      2.0")
  buf = IOBuffer(data)
  @test_throws ArgumentError DGFun.readarray(buf)
  seek(buf,0)
  @test_throws ArgumentError DGFun.readarray(buf; maxbuf=2048)
  seek(buf,0)
  a = DGFun.readarray(buf; maxbuf=2052)
  @test a == [1.0 2.0]

  data = [1/3*ones(500); pi*ones(500); nextfloat(0.0)*ones(500)]
  buf = IOBuffer()
  for i in data
    print(buf, i)
    print(buf, " ")
  end
  seek(buf, 0) # reset buffer
  a = DGFun.tseq(buf, zeros(0))
  @test a == data

  data = [1/3*ones(500); pi*ones(500); nextfloat(0.0)*ones(500)]
  buf = IOBuffer()
  for i in data
    print(buf, i)
    print(buf, "  ")
  end
  seek(buf, 0) # reset buffer
  a = DGFun.tseq(buf, zeros(0))
  @test a == data

  data = randn(3,2)
  buf = IOBuffer()
  writedlm(buf, data, "\t ")
  seek(buf, 0) # reset
  a = DGFun.readarray(buf)
  @test a == data

  buf = IOBuffer("1.0")
  a = DGFun.readarray(buf)
  @test a == ones(1,1)

  data = convert(Array{UInt8}, "1.0"*("0"^2048))
  buf = IOBuffer(data)
  @test_throws ArgumentError a = DGFun.readarray(buf)

  data = convert(Array{UInt8}, "1.0 1.0\n1.0"*("0"^2048)*"      2.0")
  buf = IOBuffer(data)
  a = DGFun.readarray(buf; maxbuf = 2052)
  @test a == [1.0 1.0; 1.0 2.0]

end


text:
module Sdistance
include("./distance_function.jl")
include("./floodfill.jl")
include("./utils/utils.jl")
include("./environments.jl")
import .DistanceFunction: create_signed_distance_function_multiprocess, create_signed_distance_function, distanceToCurve, create_distance_function, create_distance_function_multiprocess
import .Floodfill: signining_field
import .Utils: is_jordan_curve, interpolation
import CSV, DataFrames, Plots, DelimitedFiles, Luxor, BenchmarkTools, TimerOutputs
using CSV, DataFrames, Plots, DelimitedFiles, Luxor, BenchmarkTools, TimerOutputs
const tmr = TimerOutput()

function P(_phi, _x, _y, N, L, _gamma)
    _phi = create_distance_function_multiprocess(_x, _y, _gamma)
    signining_field(_phi, N + 1, L)
    return _phi
end
precompile(P, (Array, Array, Array, Int, Float64, Array))

"""
    benchmark用のためのmethod
    - floodfill使用時のmultiprocessとsingleprocessの比較
"""
function benchmark_floodfill(N::Int = 1000, _csv_datafile::String = "../test/mock_csv_data/interface.csv", multiprocess::Bool = true)
    exetimes = 3
    runtime = 0
    if multiprocess
        for i = 1:exetimes
            _phi, time = @timed signedDistance2D(_csv_datafile, N, "multi")
            runtime += time
        end
    else
        for i = 1:exetimes
            _phi, time = @timed signedDistance2D_singleprocess(_csv_datafile, N, "multi")
            runtime += time
        end
    end

    return (runtime / exetimes)
end

"""
benchmark用のためのmethod
- single curveの時のisinsideメソッド使用時ののmultiprocessとsingleprocessの比較
"""
function benchmark_singlecurves_isinside(N::Int = 1000, _csv_datafile::String = "../test/mock_csv_data/interface.csv", multiprocess::Bool = true)
    exetimes = 3
    runtime = 0
    if multiprocess
        for i = 1:exetimes
            _phi, time = @timed signedDistance2D(_csv_datafile, N, "single")
            runtime += time
        end
    else
        for i = 1:exetimes
            _phi, time = @timed signedDistance2D_singleprocess(_csv_datafile, N, "single")
            runtime += time
        end
    end

    return (runtime / exetimes)
end

"""
    csv_datafile::Union{String, DataFrame}
    N::Int
    curves::Union{String, Nothing}
"""
function signedDistance2D(csv_datafile::Union{String,DataFrame}, N::Int = 100, curves::Union{String,Nothing} = nothing)
    #===  case: multiple circles ===#
    if curves == "multi"
        # create the computational domain
        _phi = zeros(Float64, N + 1, N + 1)
        # loading the csv file(the circle data)
        _gamma = Array{Any}(undef, 0, 2)
        # ganma曲線 のデータの読み込み
        _gamma = readdlm(csv_datafile, ',', Float64)
        if length(_gamma) < 200
            _gamma = interpolation(_gamma, Int(floor(log(length(_gamma) / 2, N^1.5))) + 2, true)
        elseif Int(floor(log(length(_gamma) / 2, N^1.5))) > 0
            _gamma = interpolation(_gamma, Int(floor(log(length(_gamma) / 2, N^1.5))), true)
        end

        _x = [i for i = -L:2*L/N:L] # len:N+1
        _y = [i for i = -L:2*L/N:L] # len:N+1
        println("Data format: Multi curves\nThe CSV data size: ", size(_gamma))
        _phi = create_distance_function_multiprocess(_x, _y, _gamma)
        signining_field(_phi, N + 1, L)
        return _phi

        #=== case: simple circle ===#
    else
        # create the computational domain
        _phi = zeros(Float64, N + 1, N + 1)
        # loading the csv file(the circle data)
        _gamma = readdlm(csv_datafile, ',', Float64)
        _x = [i for i = -L:2*L/N:L] # len:N+1
        _y = [i for i = -L:2*L/N:L] # len:N+1

        is_jordan_curve(_gamma) # TODO: 丁寧なError messageを付与
        _gamma = interpolation(_gamma, Int(floor(log(length(_gamma) / 2, 2 * N))) + 1, false)
        println("Data format: The jordan curve\nThe CSV data size: ", size(_gamma))
        _phi = create_signed_distance_function_multiprocess(_x, _y, _gamma) # parallel processing
        return _phi
    end
end

function signedDistance2D_singleprocess(csv_datafile::Union{String,DataFrame}, N::Int = 100, curves::Union{String,Nothing} = nothing)
    #===  case: double circle ===#
    if curves == "multi"
        # create the computational domain
        _phi = zeros(Float64, N + 1, N + 1)
        _gamma = Array{Any}(undef, 0, 2)

        # ganma曲線 のデータの読み込み
        _gamma = readdlm(csv_datafile, ',', Float64)
        # _gamma = interpolation(_gamma, 3 + floor(Int, N / 250), true)
        if length(_gamma) < 200
            _gamma = interpolation(_gamma, Int(floor(log(length(_gamma) / 2, N^1.5))) + 2, true)
        elseif Int(floor(log(length(_gamma) / 2, N^1.5))) > 0
            _gamma = interpolation(_gamma, Int(floor(log(length(_gamma) / 2, N^1.5))), true)
        end

        println("Data format: Multi curves\nThe CSV data size: ", size(_gamma))

        _x = [i for i = -L:2*L/N:L] # len:N+1
        _y = [i for i = -L:2*L/N:L] # len:N+1

        _phi = create_distance_function(_x, _y, _gamma)
        signining_field(_phi, N + 1, L, false)
        return _phi

        #=== case: simple circle ===#
    else
        # create the computational domain
        _phi = zeros(Float64, N + 1, N + 1)
        # ganma曲線 のデータの読み込み
        _gamma = readdlm(csv_datafile, ',', Float64)
        _x = [i for i = -L:2*L/N:L] # len:N+1
        _y = [i for i = -L:2*L/N:L] # len:N+1

        is_jordan_curve(_gamma) # TODO: 丁寧なError messageを付与
        _gamma = interpolation(_gamma, Int(floor(log(length(_gamma) / 2, 2 * N))) + 1, false)

        println("Data format: The jordan curve\nThe CSV data size: ", size(_gamma))
        _phi = create_signed_distance_function(_x, _y, _gamma) # single processing
        return _phi
    end
end
precompile(signedDistance2D, (Union{String,DataFrame}, Int, Union{String,Nothing}))

export signedDistance2D, benchmark_floodfill, benchmark_singlecurves_isinside
end

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
"""
    StackFrameCategory(modcat=FlameGraphs.default_modcat,
                       loccat=FlameGraphs.default_loccat,
                       colorbg=colorant"white",
                       colorfont=colorant"black")

Colorize stackframes based on their category.

`modcat(mod::Module)` should return a color based on the stackframe's module, or `nothing`
if it cannot categorize the stack frame based on the module.

`loccat(sf::StackFrame)` must return a color. It can use any of the fields of the stackframe,
but `func`, `file`, `line`, and `from_c` might be common choices.

`colorbg` is the background color, and `colorfont` stores the choice of font color.

# Examples

```julia
using Plots, Profile, FlameGraphs
@profile plot(rand(5))    # "time to first plot"
g = flamegraph(C=true)
img = flamepixels(StackFrameCategory(), g)
```

Or you can tweak the coloration yourself:
```julia
function modcat(mod)
    mod == Plots && return colorant"purple"
    return nothing
end
img = flamepixels(StackFrameCategory(modcat), g)
```
"""
struct StackFrameCategory
    modcat   # categorization based on module (color if categorized, otherwise `nothing`)
    loccat   # categorization based on location
    colorbg::RGB{N0f8}
    colorfont::RGB{N0f8}
end

"""
    default_modcat(mod::Module)

Returns dark gray for `Core.Compiler`, light gray for `Core`, light blue for `Base`,
and otherwise returns `nothing`.
"""
function default_modcat(mod::Module)
    (mod === Core.Compiler || pm(mod) === Core.Compiler) && return colorant"gray60"
    (mod === Core || pm(mod) === Core) && return colorant"gray30"
    (mod === Base || pm(mod) === Base) && return colorant"lightblue"
    return nothing
end

"""
    default_loccat(sf::StackFrame)

Returns yellow for LLVM, orange for any other `ccall`, dark gray for anything in `./compiler`,
light blue for Base code, and red for anything else.
"""
function default_loccat(sf::StackFrame)
    file = String(sf.file)
    occursin("LLVM", file) && return colorant"yellow"
    sf.from_c && return colorant"orange"
    occursin("compiler", file) && return colorant"gray60"
    startswith(file, "./") && return colorant"lightblue"
    return colorant"red"   # uncategorized
end

StackFrameCategory(modcat = default_modcat,
                   loccat = default_loccat) =
    StackFrameCategory(modcat, loccat, colorant"white", colorant"black")

function (sfc::StackFrameCategory)(s::Symbol)
    s === :bg && return sfc.colorbg
    s === :font && return sfc.colorfont
    throw(ArgumentError("unrecognized color id ", s))
end

function (sfc::StackFrameCategory)(nextidx, j, nodedata)
    sf = nodedata.sf
    if isdefined(sf, :linfo)
        mi = sf.linfo
        if isa(mi, Core.MethodInstance)
            def = mi.def
            if isa(def, Module)
                mod = def
            else
                mod = def.module
            end
            col = sfc.modcat(mod)
            col !== nothing && return col
        end
    end
    return sfc.loccat(sf)
end

function pm(mod)
    pmod = parentmodule(mod)
    while mod !== pmod && mod !== Base && mod !== Core && mod !== Core.Compiler
        mod = pmod
        pmod = parentmodule(mod)
    end
    return mod
end


text:
using DiffEqBayesTuring, BenchmarkTools, Random
using OrdinaryDiffEq, RecursiveArrayTools, Distributions
using ParameterizedFunctions, Turing
using Plots
gr(fmt=:png)

Random.seed!(123)

### Defining the problem.

# The [FitzHugh-Nagumo model](https://en.wikipedia.org/wiki/FitzHugh%E2%80%93Nagumo_model) is a simplified version of [Hodgkin-Huxley model](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model) and is used to describe an excitable system (e.g. neuron).

fitz = @ode_def FitzhughNagumo begin
  dv = v - v^3/3 -w + l
  dw = τinv*(v +  a - b*w)
end a b τinv l

prob_ode_fitzhughnagumo = ODEProblem(fitz,[1.0,1.0],(0.0,10.0),[0.7,0.8,1/12.5,0.5])
sol = solve(prob_ode_fitzhughnagumo, Tsit5())

# Data is genereated by adding noise to the solution obtained above.

t = collect(range(1,stop=10,length=10))
sig = 0.20
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))

### Plot of the data and the solution.

scatter(t, data[1,:])
scatter!(t, data[2,:])
plot!(sol)

### Priors for the parameters which will be passed for the Bayesian Inference

priors = [truncated(Normal(1.0,0.5),0,1.5),truncated(Normal(1.0,0.5),0,1.5),truncated(Normal(0.0,0.5),0.0,0.5),truncated(Normal(0.5,0.5),0,1)]

### Benchmarks

#### Turing.jl backend

brt = turing_inference(prob_ode_fitzhughnagumo,
  Tsit5(), t,data,priors;num_samples = 2500, num_chains=4)

@btime bayesian_result_turing = turing_inference(prob_ode_fitzhughnagumo,
  Tsit5(), t,data,priors;num_samples = 2500, num_chains=4)

t = @elapsed turing_inference(prob_ode_fitzhughnagumo,
  Tsit5(), t,data,priors;num_samples = 2500, num_chains=4)


text:
# Autogenerated wrapper script for USRPHardwareDriver_jll for x86_64-linux-gnu-cxx11
export libuhd, uhd_adc_self_cal, uhd_cal_rx_iq_balance, uhd_cal_tx_dc_offset, uhd_cal_tx_iq_balance, uhd_config_info, uhd_find_devices, uhd_image_loader, uhd_usrp_probe

using boost_jll
using libusb_jll
JLLWrappers.@generate_wrapper_header("USRPHardwareDriver")
JLLWrappers.@declare_library_product(libuhd, "libuhd.so.4.1.0")
JLLWrappers.@declare_executable_product(uhd_adc_self_cal)
JLLWrappers.@declare_executable_product(uhd_cal_rx_iq_balance)
JLLWrappers.@declare_executable_product(uhd_cal_tx_dc_offset)
JLLWrappers.@declare_executable_product(uhd_cal_tx_iq_balance)
JLLWrappers.@declare_executable_product(uhd_config_info)
JLLWrappers.@declare_executable_product(uhd_find_devices)
JLLWrappers.@declare_executable_product(uhd_image_loader)
JLLWrappers.@declare_executable_product(uhd_usrp_probe)
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, libusb_jll)
    JLLWrappers.@init_library_product(
        libuhd,
        "lib/libuhd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        uhd_adc_self_cal,
        "bin/uhd_adc_self_cal",
    )

    JLLWrappers.@init_executable_product(
        uhd_cal_rx_iq_balance,
        "bin/uhd_cal_rx_iq_balance",
    )

    JLLWrappers.@init_executable_product(
        uhd_cal_tx_dc_offset,
        "bin/uhd_cal_tx_dc_offset",
    )

    JLLWrappers.@init_executable_product(
        uhd_cal_tx_iq_balance,
        "bin/uhd_cal_tx_iq_balance",
    )

    JLLWrappers.@init_executable_product(
        uhd_config_info,
        "bin/uhd_config_info",
    )

    JLLWrappers.@init_executable_product(
        uhd_find_devices,
        "bin/uhd_find_devices",
    )

    JLLWrappers.@init_executable_product(
        uhd_image_loader,
        "bin/uhd_image_loader",
    )

    JLLWrappers.@init_executable_product(
        uhd_usrp_probe,
        "bin/uhd_usrp_probe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using InteractionNetworkModels, Plots, Distributions
using StatsPlots, Plots.Measures, StatsBase, LaTeXStrings

E = [[1,2,1,2],
    [1,2,1],
    [3,4,3],
    [3,4],
    [1,2]
]

d = MatchingDist(FastLCS(100))

K_inner = DimensionRange(1,50)
K_outer = DimensionRange(1,50)
V = 1:30
model = SIM(
    E, 4.4,
    d,
    V, K_inner, K_outer)

mcmc_sampler = SimMcmcInsertDelete(
    ν_ed=7, ν_td=1, β=0.7,
    len_dist=TrGeometric(0.9, 1, K_inner.u),
    burn_in=2000, lag=50
)

@time test = mcmc_sampler(
    model, desired_samples=500)

plot(test)
summaryplot(test)

@time mcmc_out = mcmc_sampler(
    model,
    desired_samples=50,
    lag=500,
    burn_in=10000
    )

plot(mcmc_out)
summaryplot(mcmc_out)

data = mcmc_out.sample
E_prior = SIM(E, 0.1, model.dist, model.V, model.K_inner, model.K_outer)
γ_prior = Uniform(0.5,7.0)

target = SimPosterior(data, E_prior, γ_prior)

# Construct posterior sampler
posterior_sampler = SimIexInsertDelete(
    mcmc_sampler,
    len_dist=TrGeometric(0.9,1,100),
    ν_ed=1, ν_td=1,
    β=0.7, ε=0.2
)

# Dispersion conditional

E_restr = E[1:2]
posterior_out = posterior_sampler(
    target,
    desired_samples=200, lag=1, burn_in=0,
    γ_init=4.6, S_init=E[1:1]
)
posterior_out.S_sample
plot(posterior_out)

# Nested means

x_num_inters = [zeros(200) for i in eachindex(E)]

for i in eachindex(E)
    draw_sample_gamma!(
        x[i],
        posterior_sampler,
        target,
        E[1:i],
        burn_in=250, lag=5,
        γ_init = 4.4
    )
end

density(
    x,
    label=[L"\mathcal{E}_{1:%$(i)}" for j=1:1, i=eachindex(E)],
    legend_font_pointsize=10,
    legend_font_valign=:bottom,
    xlabel="γ",
    ylabel="(Conditional) Posterior Density"
)

# Growing dimensions

modes = InteractionSequenceSample{Int}()
for i in eachindex(x)
    tmp_count = 0
    path_ind = findlast(cumsum(length.(E)) .< i)
    if path_ind == nothing
        E_restr = [E[1][1:i]]
    else
        entry_ind = i - cumsum(length.(E))[path_ind]
        E_restr = E[1:path_ind]
        push!(E_restr, E[path_ind+1][1:entry_ind])
    end
    push!(modes, E_restr)
end
x = [zeros(200) for i in eachindex(modes)]
for i in eachindex(modes)
    draw_sample_gamma!(
        x[i],
        posterior_sampler,
        target,
        modes[i],
        burn_in=250, lag=5,
        γ_init = 4.4
    )
end

density(
    x,
    label=["Total dim. = $(i)" for j=1:1, i=eachindex(x)],
    legend=:outerright,
    xlabel="γ",
    ylabel="(Conditional) Posterior Density"
)

# Using dependent sampler
posterior_sampler_dep = SimIexInsertDeleteDependent(
    mcmc_sampler,
    len_dist=TrGeometric(0.9,1,100),
    ν_ed=1, ν_td=1,
    β=0.7, ε_ed=0.07, ε_td=0.2
)

posterior_out = posterior_sampler_dep(
    target,
    desired_samples=400, lag=1, burn_in=0,
    γ_init=4.6, S_init=E[1:1]
)
plot(posterior_out, E)
posterior_out.S_sample


text:
#
# Copyright: Rob J Goedman (2015)
#
# Defines the composite Parameters type. The params variable will be created
# after checking the input Dictionary for correctness and completeness.

struct Parameters
  #
  # CONSTANTS
  #
  # Conversion factor from degrees to radians
  #
	dtor::Float64
  #
  # Gravitational acceleration
  #
	g::Float64
  #
  # PARAMETERS
  #
  # Used in construction of element mesh.
  # ratio = segment length / segment outer diameter
  #
	ratio::Float64
  #
  # INPUTS
  #
  wob::Float64
  maxiters::Int64
  remove_bit_constraints::Bool
  use_geometric_correction::Bool
  #
  # Derived from inputs
  #
  has_bit::Bool
  has_rig::Bool

  noofelements::Int64
  noofnodes::Int64
  noofiters::Vector{Int64}

  elements::Array{BHAtp.Element, 1}
  nodes::Array{BHAtp.Node, 1}

  f::Array{Float64, 3}
  bf::Array{Float64, 3}
  fn::Vector{Float64}
  fs::Vector{Float64}
  friction::Vector{Float64}
  del::Array{Float64, 3}
  fvec::Array{Float64, 3}

  sainit::Array{Float64, 3}
  sb1init::Array{Float64, 3}
  sb2init::Array{Float64, 3}
  scinit::Array{Float64, 3}
  siinit::Array{Float64, 3}
  sa::Array{Float64, 3}
  sb1::Array{Float64, 3}
  sb2::Array{Float64, 3}
  sc::Array{Float64, 3}
  si::Array{Float64, 3}
  gm::Array{Float64, 3}

  fdel::Matrix{Float64}
  rdel::Matrix{Float64}
  exc::Matrix{Float64}
  sasc::Matrix{Float64}
  delta::Matrix{Float64}

  dlerrormax::Float64

end

function Parameters(data::Dict{Symbol, Any})

	ratio =
	if :ratio in keys(data)
		ratio = data[:ratio]::Float64
	else
		ratio = 0.6
	end

  noofelements = data[:noofelements]
  noofnodes = data[:noofelements] + 1

  elements = [BHAtp.Element() for i in 1:noofelements]
  nodes = [BHAtp.Node() for i in 1:noofnodes]

  f = zeros(Float64, 6, 3, noofnodes)
  bf = zeros(Float64, 2, 3, noofnodes)
  fn = zeros(Float64, noofnodes)
  fs = zeros(Float64, noofnodes)
  friction = zeros(Float64, noofnodes)
  del = zeros(Float64, 6, 3, noofnodes)
  fvec = zeros(Float64, 6, 3, noofnodes)

  sainit = zeros(Float64, 6, 6, noofnodes)
  sb1init = zeros(Float64, 6, 6, noofnodes)
  sb2init = zeros(Float64, 6, 6, noofnodes)
  scinit = zeros(Float64, 6, 6, noofnodes)
  siinit = zeros(Float64, 6, 6, noofnodes)
  sa = zeros(Float64, 6, 6, noofnodes)
  sb1 = zeros(Float64, 6, 6, noofnodes)
  sb2 = zeros(Float64, 6, 6, noofnodes)
  sc = zeros(Float64, 6, 6, noofnodes)
  si = zeros(Float64, 6, 6, noofnodes)
  gm = zeros(Float64, 12, 12, noofnodes)

  fdel = zeros(Float64, 6, 3)
  rdel = zeros(Float64, 6, 3)
  exc = zeros(Float64, 2, 3)
  sasc = zeros(Float64, 6, 12)
  delta = zeros(Float64, 12, 3)

  dtor = pi / 180.0

  Parameters(
  	# dtor: 1 degree = 2 * pi / 360 radians = 0.0174532925 radians
  	dtor,

  	# g: gravitational constant [m/sec^2]
  	9.80665,

  	# ratio used to subdivide segments into (finite) elements
  	ratio,

    # wob, already converted from tonnes
    data[:wob],

    # Max number of iterations
    data[:maxiters],

    # Set remove_bit_constraints in params (or better still, use tp!)
    :remove_bit_constraints in keys(data) ? data[:removebitcontsraints] : false,

    # Flag that geometriccorrection is to be applied
    :geometriccorrection in keys(data) ? data[:geometriccorrection] : false,

    # Set has_bit field in params, CheckInput.jl checked that if a bit is present,
    # it must be the first segment defined in the BHAJ.
    typeof(data[:segments][1]) == Bit ? true : false,

    # Set has_rig to indicate top element needs to be fix in central position
    typeof(data[:segments][data[:noofsegments]]) == Rig ? true : false,

    noofelements,
    noofnodes,
    [0],
    elements, nodes,
    f, bf, fn, fs, friction, del, fvec,
    sainit, sb1init, sb2init, scinit, siinit,
    sa, sb1, sb2, sc, si, gm,
    fdel, rdel, exc, sasc, delta,

    data[:dlerrormax]
  )
end

function parameters_show(io::IO, p::Parameters, compact::Bool)
  if compact
    println(io, [getfield(p, f) for f in fieldnames(p)])
  else
    println("Parameters(")
    println("  # Conversion factors:")
    println("  dtor                     = ", p.dtor)
    println("  g                        = ", p.g)
    println()
    println("  # Inputs:")
    println("  ratio                    = ", p.ratio)
    println("  wob                      = ", p.wob)
    println("  dlerrormax               = ", p.dlerrormax)
    println("  drerrormax               = ", p.drerrormax)
    println("  maxiters                 = ", p.maxiters)
    println("  noofiters                = ", p.noofiters[1])
    println()
    println("  # Mesh creation:")
    println("  noofelements             = ", p.noofelements)
    println("  noofnodes                = ", p.noofnodes)
    println()
    println("  # Execution flags:")
    println("  has_bit                  = ", p.has_bit)
    println("  remove_bit_constraints   = ", p.remove_bit_constraints)
    println("  has_rig                  = ", p.has_rig)
    println("  use_geometric_correction = ", p.use_geometric_correction)
    println(")")
  end
end

show(io::IO, p::Parameters) = parameters_show(io, p, false)
showcompact(io::IO, p::Parameters) = parameters_show(io, p, true)


text:
function checkdownload_ionosphere()
    if !isfile(joinpath(instances_dir, "ionosphere.data"))
                println("Downloading ionosphere dataset to: ", pwd(), "...")

        run(`wget https://archive.ics.uci.edu/ml/machine-learning-databases/ionosphere/ionosphere.data`)

        mv("ionosphere.data", joinpath(instances_dir, "ionosphere.data"))
        println("Moving ionosphere dataset to: ", joinpath(instances_dir, "ionosphere.data"), "...")

        println("All done.")
    end
    return
end


text:
module Batched

include("batched_arrays.jl")
include("broadcast.jl")

end # module


text:
export  DropTable,
        addtables!,
        droptable

"""
```
DropTable( names::Vector{Var} )
DropTable( names::Vector{T} ) where T <: AbstractString
DropTable( names::Union{Var, AbstractString}... )
```
This mutable struct represents a `DROP TABLE` SQL statement where `names` are the names of the tables to drop.
"""
mutable struct DropTable <: SQLStatement
    names::Vector{Var}
end

DropTable( names::Vector{T} ) where T <: AbstractString = Var.(names) |> DropTable
DropTable( names::SVar... ) = Var.( names |> collect ) |> DropTable

function Base.show( io::IO, dt::DropTable )
    isempty(dt.names) && return
    print( io, "DROP TABLE ", join( unique(dt.names), ", " ) )
end

"""
```
addtables!( names::Vector{Var} )
addtables!( names::Vector{T} ) where T <: AbstractString
addtables!( names::Union{Var, AbstractString}... )
```
This function permits the user to add a command to add the tables in `names` to the list of tables to drop of a `DropTable` struct.

The output of this function is a function that takes a `DropTable` as argument and returns the updated argument, such that Julia's function chaining syntax can be used.

Example:
```
sq |> addtables!( "table1", "table2" )
```
adds two tables to drop to the `DropTable` struct: `table1` and `table2`.
"""
addtables!( names::Vector{Var} ) = dt::DropTable -> begin
    append!( dt.names, names )
    dt
end
addtables!( names::Vector{T} ) where T <: AbstractString = Var.(names) |> addtables!
addtables!( names::SVar... ) = Var.( names |> collect ) |> addtables!

function sanitise!( dt::DropTable, conn::Conn, muteerror::Bool )
    unique!(dt.names)
    muteerror || return

    dbtables = Var.( conn |> gettables )
    filter!( name -> name ∈ dbtables, dt.names )
end

"""
```
droptable( conn::Conn; muteerror::Bool=true )
```
This function permits the user to execute the SQL equivalent of a `DropTable` struct on the database that `conn` connects to. If the flag `muteerror` is `true`, the `DropTable` struct gets sanitised so only existing tables

The output of this function is a function that takes a `DropTable` as argument and returns the result of the executed SQL statement.

Example:
```
sq |> droptable(conn)
```
"""
droptable( conn::Conn; muteerror::Bool=true ) = dt::DropTable -> begin
    sanitise!( dt, conn, muteerror )
    LibPQ.execute( conn, string( dt, ";" ) )
end


text:
#==============================================================================
    Open an HTML string in browser as a temp file

    author:     mahiki@users.noreply.github.com
=============================================================================#

"""
    browserview_html(html_string::String)

Open an html string as a temp file in local system, usually defaul app is a browser.
Useful for interactive work processing html data.
"""
function browserview_html(html_string::String)
    f = tempname()
    io = open("$f.html", "w")
    print(io, html_string)
    close(io)
    run(`open $f.html`)
end

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/implicit_routines/contact.jl
# Definitions for contract functions used in implicit relaxations.
#############################################################################

"""
$(SIGNATURES)

Performs a single step of the parametric method associated with `t` assumes that
the inputs have been preconditioned.
"""
function contract! end

"""
$(TYPEDSIGNATURES)

Applies the Gauss-Siedel variant of the Newton type contractor.
"""
function contract!(t::NewtonGS, d::MCCallback{FH,FJ,C,PRE,N,T}) where {FH <: Function,
                                                                       FJ <: Function,
                                                                       C, PRE, N,
                                                                       T<:RelaxTag}
    S = zero(MC{N,T})
    @. d.x0_mc = d.x_mc
    for i = 1:d.nx
        S = zero(MC{N,T})
        for j = 1:d.nx
            if (i !== j)
                @inbounds S += d.J[i,j]*(d.x_mc[j] - d.z_mc[j])
            end
        end
        @inbounds d.x_mc[i] = d.z_mc[i] - (d.H[i] + S)*McCormick.inv1(d.J[i,i], 1.0/d.J[i,i].Intv)
        @inbounds d.x_mc[i] = final_cut(d.x_mc[i], d.x0_mc[i])
    end
    return
end

"""
$(TYPEDSIGNATURES)

Applies the componentwise variant of the Krawczyk type contractor.
"""
function contract!(t::KrawczykCW, d::MCCallback{FH,FJ,C,PRE,N,T}) where {FH <: Function,
                                                                         FJ <: Function,
                                                                         C, PRE, N,
                                                                         T<:RelaxTag}
    S::MC{N,T} = zero(MC{N,T})
    @. d.x0_mc = d.x_mc
    for i=1:d.nx
        S = zero(MC{N,T})
        for j=1:d.nx
            if (i !== j)
                @inbounds S -= (d.J[i,j])*(d.x_mc[j] - d.z_mc[j])
            else
                @inbounds S += (one(MC{N,T}) - d.J[i,j])*(d.x_mc[j] - d.z_mc[j])
            end
        end
        @inbounds d.x_mc[i] = d.z_mc[i] - d.H[i] + S
        @inbounds d.x_mc[i] = final_cut(d.x_mc[i], d.x0_mc[i])
    end
    return
end


text:
@testset "Continous Weight" begin

    T = Float64

    ## known set of modified moments
    P = ShiftedLegendre{T}
    f(x) = -log(x)
    q0 = ContinuousWeight(P,f,[1], :x)
    WW = typeof(q0)
    function SpecialPolynomials.modified_moment(w::W, j) where {W <: WW}
        iszero(j) && return 1.0
        (-1)^j * gamma(j+1)^2/(j*(j+1)*gamma(2j+1))
    end

    q0 = ContinuousWeight(P, f, [1], :x)
    q1 = ContinuousWeight(P, f, [0,1], :x)
    q2 = ContinuousWeight(P, f, [0,0,1], :x)
    q3 = ContinuousWeight(P, f, [0,0,0,1], :x)
    q4 = ContinuousWeight(P, f, [0,0,0,0,1], :x)
    q5 = ContinuousWeight(P, f, [0,0,0,0,0,1], :x)
    qs = convert.(Polynomial, (q0,q1,q2,q3,q4,q5))
    a, b = first(domain(q0)), last(domain(q0))
    for i in eachindex(qs)
        for j in eachindex(qs)
            j <= i  && continue
            @test abs(SpecialPolynomials._quadgk(x -> qs[i](x) * qs[j](x) * f(x),  first(domain(q0)), last(domain(q0)))) <= 100*sqrt(eps(T))
        end
    end

    ## non-known set
    T = Float64
    P = Chebyshev{T}
    f(x) = SpecialPolynomials.weight_function(ChebyshevU{T})(x)
    q0 = ContinuousWeight(P, f, [1], :x)
    q1 = ContinuousWeight(P, f, [0,1], :x)
    q2 = ContinuousWeight(P, f, [0,0,1], :x)
    q3 = ContinuousWeight(P, f, [0,0,0,1], :x)
    q4 = ContinuousWeight(P, f, [0,0,0,0,1], :x)
    q5 = ContinuousWeight(P, f, [0,0,0,0,0,1], :x)
    qs = convert.(Polynomial, (q0,q1,q2,q3,q4,q5))
    a, b = first(domain(q0)), last(domain(q0))
    for qi in qs
        for qj in qs
            j <= i && continue
            @test abs(SpecialPolynomials._quadgk(x -> qi(x) * qj(x) * f(x), a, b)) <= 100*sqrt(eps(T))
        end
    end

end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

# get_receive = load_op_and_grad("build/libGetReceive", "get_receive")
get_receive = load_op_and_grad("build/libGetReceive", "get_receive")

################## End Load Operator ##################
vx = constant(ones(2, 9))
vy = constant(2ones(2, 9))
sigmaxx = constant(3ones(2, 9))
sigmayy = constant(4ones(2, 9))
sigmaxy = constant(5ones(2, 9))

rcvi = constant([1])
rcvj = constant([1])
rcvtype = constant([1])
nx = constant(1)
ny = constant(1)

# TODO: specify your input parameters
u = get_receive(vx,vy,sigmaxx,sigmayy,sigmaxy,rcvi,rcvj,rcvtype,nx,ny)
sess = tf.Session()
init(sess)
run(sess, u)
# error()

# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(m)
    return sum(get_receive(vx, m, sigmaxx,sigmayy,sigmaxy,rcvi,rcvj,rcvtype,nx,ny)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(ones(2,9))
v_ = rand(2,9)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
export Shape, Tria, Rect, Hexa
export Cell, Grid
export TriaGrid, PeriodicTriaGrid, RectGrid, HexaGrid, PeriodicRectGrid, PeriodicHexaGrid
export connect!, disconnect!, is_connected, get_deadends, remove_deadends!
export north, south, west, east
export recursive_backtracker!, wall_locations
export draw_grid, draw_solution, draw_distance
export dijkstra, astar


text:
export create_paramop, apply_paramop!

create_paramop(paramKd::AbsArrComplexF,  # material parameter array
               isfwd_in::AbsVecBool,  # isfwd_in[w] = true if w-component of input field is forward-averaged to be interpolated at voxel corners
               ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with (defined at field voxel edges); vectors of length N
               ∆l′⁻¹::NTuple{K,AbsVecNumber},  # inverse of line segments to divide by (defined at field voxel corners); vectors of length N
               isbloch::AbsVecBool,  # for K = 3, boundary conditions in x, y, z
               e⁻ⁱᵏᴸ::AbsVecNumber=ones(K);  # for K = 3, Bloch phase factor in x, y, z
               order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
               ) where {K} =
    create_paramop(paramKd, SVec{K}(isfwd_in), ∆l, ∆l′⁻¹, SVec{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ); order_cmpfirst)

# Version for scalar material parameters
function create_paramop(paramKd::AbsArrComplexF{K₊₂};  # scalar material parameter array
                        order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
                        ) where {K₊₂}
    K = K₊₂ - 2
    N = SVec(ntuple(k->size(paramKd,k)-1, Val(K)))  # for -1, see comments about size of paramKd in assignment.jl
    kdiag = 0
    paramop = create_paramop_supdiag(paramKd, kdiag, N; order_cmpfirst)

    return paramop
end

# Below, N can be retrieved from the size of paramKd or ∆l, but I pass it to make the size
# more explicit.  This also makes the argument list similar to create_mean()'s, which
# includes N because it allows omitting ∆l and ∆l′ and using ∆l = ∆l′ = ones.((N...,)) as
# the default arguments.
function create_paramop(paramKd::AbsArrComplexF{K₊₂},  # material parameter array
                        isfwd_in::SBool{K},  # isfwd_in[w] = true if w-component of input field is forward-averaged to be interpolated at voxel corners
                        ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with (defined at field voxel edges); vectors of length N
                        ∆l′⁻¹::NTuple{K,AbsVecNumber},  # inverse of line segments to divide by (defined at field voxel corners); vectors of length N
                        isbloch::SBool{K},  # for K = 3, boundary conditions in x, y, z
                        e⁻ⁱᵏᴸ::SNumber{K};  # for K = 3, Bloch phase factor in x, y, z
                        order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
                        ) where {K,K₊₂}
    @assert K₊₂==K+2

    Kf = size(paramKd, K+1)  # field dimension
    @assert size(paramKd,K+2)==Kf

    N = SVec(length.(∆l))
    length.(∆l′⁻¹)==N.data || error("Corresponding entries of ∆l and ∆l′⁻¹ should have same length.")
    M = prod(N)  # number voxels

    # Make sure Kf ≥ 2; the material parameters should be tensorial.
    Kf==1 && @error "paramKd should represent array of tensorial material parameters"

    # Following Oskooi et al.'s 2009 Optics Letters paper, off-diagonal entries of
    # material parameter tensors (e.g., ε) are evaluated at the corners of voxels whose
    # edges are field lines (e.g., E).  Therefore, if w-normal voxel faces are primal
    # grid planes, the w-component of the input fields need to be averaged in the
    # backward direction to be interpolated at voxel corners.
    #
    # Once these interpolated input fields (e.g., Ex) are multiplied with off-diagonal
    # entries (e.g., εzx) of the material parameter tensor, we get the output fields
    # (e.g., Dzx of Dz = Dzx + Dzy + Dzz = εzx Ex + εzy Ey + εzz Ez) in the direction
    # normal (say v-direction) to the input fields.  These output fields are still at
    # the voxel corners.  Now, if the v-normal voxel faces are primal grid planes, the
    # v-component of the output fields need to be averaged in the forward direction to
    # be interpolated at voxel edges.
    #
    # In summary,
    #
    # - to obtain the fields to feed to the off-diagonal entries of material parameter
    # tensors, the w-component of the input fields need to be forward(backward)-averaged
    # along the w-direction if the w-normal voxel faces are dual (primal) grid planes,
    #
    # - to distribute the resulting output fields in the v(≠w)-direction back to voxel
    # edges, the output fields need to be forward(backward)-averaged along the
    # v-direction if the v-normal voxel faces are primal (dual) grid planes.
    #
    # This means if the w-boundaries are E-field boundaries,
    # - the input E(H)w-field needs to be backward(forward)-averaged in the w-direction,
    # - the output E(H)w-field needs to be forward(backward)-averaged in the w-direction.
    isfwd_out = .!isfwd_in  # SBool{K}; true if output fields need to be backward-averaged

    # For the output averaging, ∆l and ∆l′ are not supplied to create_mean in order to
    # create a simple arithmetic averaging operator.  This is because the area factor
    # matrix multiplied for symmetry to the left of the material parameter matrix
    # multiplies the same area factor to the two fields being averaged.  (See my notes
    # on Jul/16/2018 in MaxwellFDFD in Agenda.)
    Mout = create_mean(isfwd_out, N, isbloch, e⁻ⁱᵏᴸ; order_cmpfirst)

    # For the input averaging, ∆l and ∆l′ are supplied to create_mean in order to create
    # a line integral averaging operator.  This is because the inverse of the length
    # factor matrix multiplied for symmetry to the right of the material parameter
    # matrix divides the two fields being averaged by different (= nonuniform) line
    # segments.  The ∆l factors multiplied inside create_minfo cancel the effect of this
    # multiplication with the nonuniform line segments.  (See my notes on Jul/16/2018 in
    # MaxwellFDFD in Agenda.)
    Min = create_mean(isfwd_in, ∆l, ∆l′⁻¹, isbloch, e⁻ⁱᵏᴸ; order_cmpfirst)

    # Off-diagonal components of material parameter tensor
    # Min and Mout are block-diagonal matrices, but create_paramop_supdiag() generates
    # block-off-diagonal matrices for 1 ≤ kdiag ≤ Kf-1.  See my bullet point entitled
    # [Update (May/13/2018)] in RN - Subpixel Smoothing.
    kdiag = 1
    paramop = create_paramop_supdiag(paramKd, kdiag, N; order_cmpfirst)
    for kdiag = 2:Kf-1  # index of diagonal of material parameter tensor
        paramop += create_paramop_supdiag(paramKd, kdiag, N; order_cmpfirst)
    end
    paramop = Mout * paramop * Min

    # Diagonal components of material parameter tensor
    kdiag = 0
    paramop += create_paramop_supdiag(paramKd, kdiag, N; order_cmpfirst)  # diagonal components of material parameter tensor

    return paramop
end

# Create the kth superdiagonal components of the parameter operator (0 ≤ k ≤ Kf-1)
#
# About the meaning of kdiag
#
# Suppose material parameter tensor is Kf-by-Kf.  Different values of kdiag take different
# entries of the material parameter tensor and used to create a matrix.  For Kf = 3, here
# are the entries (indicated by X) taken for the three different values of kdiag.
#
# kdiag = 0:
# ⎡X O O⎤
# ⎢O X O⎥
# ⎣O O X⎦
#
# kdiag = 1:
# ⎡O X O⎤
# ⎢O O X⎥
# ⎣X O O⎦
#
# kdiag = 2 (= Kf-1):
# ⎡O O X⎤
# ⎢X O O⎥
# ⎣O X O⎦
#
# So, kdiag is the index counted from the main diagonal (kdiag = 0) along the superdiagonal
# direction.  See RN - Subpixel Smoothing > [Update (May/13/2018)].  Note that if the X is
# the vw-block, then it is constructed with the vw-component of the material parameter (e.g.,
# for ε, if X is the xy-block, then it is constructed with εxy).  Therefore, whether the
# direction of the superdiagonal shift is along the row- or column-direction is not relevant.
function create_paramop_supdiag(paramKd::AbsArrComplexF{K₊₂},  # size of last two dimensions is Kf-by-Kf
                                kdiag::Integer,  # 0 ≤ kdiag ≤ Kf-1; index of diagonal of material parameter tensor to set (kdiag = 0: main diagonal)
                                N::SInt{K};  # size of grid
                                order_cmpfirst::Bool=true  # true to use Cartesian-component-major ordering for more tightly banded matrix
                                ) where {K,K₊₂}
    @assert K₊₂==K+2

    Kf = size(paramKd, K+1)  # field dimension
    @assert size(paramKd,K+2)==Kf

    # Note that paramKd's i, j, k indices run from 1 to N+1 rather than to N, so we should
    # not iterate those indices from 1 to end (= N+1).
    M = prod(N)  # number of voxels
    KfM = Kf * M
    I = VecInt(undef, KfM)
    J = VecInt(undef, KfM)
    V = VecComplexF(undef, KfM)
    n = 0

    CI = CartesianIndices(N.data)
    LI = LinearIndices(N.data)
    for nv = 1:Kf  # row index of material parameter tensor
        istr, ioff = order_cmpfirst ? (Kf, nv-Kf) : (1, M*(nv-1))  # (row stride, row offset)
        nw = mod1(nv+kdiag, Kf)  # column index of material parameter tensor
        jstr, joff = order_cmpfirst ? (Kf, nw-Kf) : (1, M*(nw-1))  # (column stride, column offset)
        for ci = CI
            n += 1
            @inbounds ind = LI[ci]  # linear index of Yee's cell

            I[n] = istr * ind + ioff
            J[n] = jstr * ind + joff
            @inbounds V[n] = paramKd[ci,nv,nw]
        end
    end

    return sparse(I, J, V, KfM, KfM)
end

apply_paramop!(G::AbsArrNumber,  # output field; in 3D, G[i,j,k,w] is w-component of G at (i,j,k)
               F::AbsArrNumber,  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
               G′::AbsArrNumber,  # first temporary storage for intermediate output; in 3D, G′[i,j,k,w] is w-component of G′ at (i,j,k)
               G′′::AbsArrNumber,  # second temporary storage for intermediate output; in 3D, 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
               paramKd::AbsArrComplexF,  # material parameter array
               isfwd_in::AbsVecBool,  # isfwd_in[w] = true if w-component of input field is forward-averaged to be interpolated at voxel corners
               ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with (defined at field voxel edges); vectors of length N
               ∆l′⁻¹::NTuple{K,AbsVecNumber},  # line segments to divide by (defined at field voxel corners); vectors of length N
               isbloch::AbsVecBool,  # for K = 3, boundary conditions in x, y, z
               e⁻ⁱᵏᴸ::AbsVecNumber=ones{K};  # for K = 3, Bloch phase factor in x, y, z
               α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
               ) where {K,OP} =
apply_paramop!(G, F, G′, G′′, Val(OP), paramKd, SVec{K}(isfwd_in), ∆l, ∆l′⁻¹, SVec{K}(isbloch), SVec{K}(e⁻ⁱᵏᴸ); α)

# Version for scalar material parameters
function apply_paramop!(G::AbsArrNumber{K₊₁},  # output field; in 3D, G[i,j,k,w] is w-component of G at (i,j,k)
                        F::AbsArrNumber{K₊₁},  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
                        ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                        paramKd::AbsArrComplexF{K₊₂};  # material parameter array
                        α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
                        ) where {K₊₁,K₊₂,OP}
    kdiag = 0
    apply_paramop_supdiag!(G, F, Val(OP), paramKd, kdiag; α)  # diagonal components of material parameter tensor

    return nothing
end

function apply_paramop!(G::AbsArrNumber{K₊₁},  # output field; in 3D, G[i,j,k,w] is w-component of G at (i,j,k)
                        F::AbsArrNumber{K₊₁},  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
                        G′::AbsArrNumber{K₊₁},  # first temporary storage for intermediate output; in 3D, G′[i,j,k,w] is w-component of G′ at (i,j,k)
                        G′′::AbsArrNumber{K₊₁},  # second temporary storage for intermediate output; in 3D, 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
                        paramKd::AbsArrComplexF{K₊₂},  # material parameter array
                        isfwd_in::SBool{K},  # isfwd_in[w] = true if w-component of input field is forward-averaged to be interpolated at voxel corners
                        ∆l::NTuple{K,AbsVecNumber},  # line segments to multiply with (defined at field voxel edges); vectors of length N
                        ∆l′⁻¹::NTuple{K,AbsVecNumber},  # line segments to divide by (defined at field voxel corners); vectors of length N
                        isbloch::SBool{K},  # for K = 3, boundary conditions in x, y, z
                        e⁻ⁱᵏᴸ::SNumber{K};  # for K = 3, Bloch phase factor in x, y, z
                        α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
                        ) where {K,K₊₁,K₊₂,OP}
    @assert K₊₁==K+1 && K₊₂==K+2

    Kf = size(paramKd, K+1)  # field dimension
    @assert size(paramKd,K+2)==Kf

    # Make sure Kf ≥ 2; the material parameters should be tensorial.
    Kf==1 && @error "paramKd should represent array of tensorial material parameters"

    isfwd_out = .!isfwd_in  # SBool{K}; true if output fields need to be backward-averaged

    apply_mean!(G′, F, Val(:(=)), isfwd_in, ∆l, ∆l′⁻¹, isbloch, e⁻ⁱᵏᴸ)

    # Off-diagonal components of material parameter tensor
    kdiag = 1
    apply_paramop_supdiag!(G′′, G′, Val(:(=)), paramKd, kdiag)
    for kdiag = 2:Kf-1  # index of diagonal of material parameter tensor
        apply_paramop_supdiag!(G′′, G′, Val(:(+=)), paramKd, kdiag)
    end
    apply_mean!(G, G′′, Val(OP), isfwd_out, isbloch, e⁻ⁱᵏᴸ; α)

    # Diagonal components of material parameter tensor
    kdiag = 0
    apply_paramop_supdiag!(G, F, Val(:(+=)), paramKd, kdiag; α)

    return nothing
end

# 3D
function apply_paramop_supdiag!(G::AbsArrNumber{4},  # output field; G[i,j,k,w] is w-component of G at (i,j,k)
                                F::AbsArrNumber{4},  # input field; F[i,j,k,w] is w-component of F at (i,j,k)
                                ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                                paramKd::AbsArrComplexF{5},  # size of last two dimensions is Kf-by-Kf
                                kdiag::Integer;  # 0 ≤ kdiag ≤ Kf-1; index of diagonal of material parameter tensor to set (kdiag = 0: main diagonal)
                                α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
                                ) where {OP}
    Nx, Ny, Nz, Kf = size(G)  # Kf: field dimension
    @assert size(paramKd)[end] == size(paramKd)[end-1] == Kf

    for nv = 1:Kf  # row index of material parameter tensor
        Gv = @view G[:,:,:,nv]
        nw = mod1(nv+kdiag, Kf)  # column index of material parameter tensor
        @threads for k = 1:Nz
            for j = 1:Ny, i = 1:Nx
                @inbounds set_or_add!(Gv, (i,j,k), paramKd[i,j,k,nv,nw] * F[i,j,k,nw], Val(OP))
                # @inbounds set_or_add!(G, (i,j,k,nv), paramKd[i,j,k,nv,nw] * F[i,j,k,nw], Val(OP))
            end
        end
    end

    return nothing
end

# 2D
function apply_paramop_supdiag!(G::AbsArrNumber{3},  # output field; G[i,j,w] is w-component of G at (i,j)
                                F::AbsArrNumber{3},  # input field; F[i,j,w] is w-component of F at (i,j)
                                ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                                paramKd::AbsArrComplexF{4},  # size of last two dimensions is Kf-by-Kf
                                kdiag::Integer;  # 0 ≤ kdiag ≤ Kf-1; index of diagonal of material parameter tensor to set (kdiag = 0: main diagonal)
                                α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
                                ) where {OP}
    Nx, Ny, Kf = size(G)  # Kf: field dimension
    @assert size(paramKd)[end] == size(paramKd)[end-1] == Kf

    for nv = 1:Kf  # row index of material parameter tensor
        Gv = @view G[:,:,nv]
        nw = mod1(nv+kdiag, Kf)  # column index of material parameter tensor
        @threads for j = 1:Ny
            for i = 1:Nx
                @inbounds set_or_add!(Gv, (i,j), paramKd[i,j,nv,nw] * F[i,j,nw], Val(OP))
                # @inbounds set_or_add!(G, (i,j,k,nv), paramKd[i,j,k,nv,nw] * F[i,j,k,nw], Val(OP))
            end
        end
    end

    return nothing
end

# 1D
function apply_paramop_supdiag!(G::AbsArrNumber{2},  # output field; G[i,w] is w-component of G at (i,)
                                F::AbsArrNumber{2},  # input field; F[i,w] is w-component of F at (i,)
                                ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
                                paramKd::AbsArrComplexF{3},  # size of last two dimensions is Kf-by-Kf
                                kdiag::Integer;  # 0 ≤ kdiag ≤ Kf-1; index of diagonal of material parameter tensor to set (kdiag = 0: main diagonal)
                                α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
                                ) where {OP}
    Nx, Kf = size(G)  # Kf: field dimension
    @assert size(paramKd)[end] == size(paramKd)[end-1] == Kf

    for nv = 1:Kf  # row index of material parameter tensor
        Gv = @view G[:,:,nv]
        nw = mod1(nv+kdiag, Kf)  # column index of material parameter tensor
        @threads for j = 1:Ny
            for i = 1:Nx
                @inbounds set_or_add!(Gv, (i,j), paramKd[i,j,nv,nw] * F[i,j,nw], Val(OP))
                # @inbounds set_or_add!(G, (i,j,k,nv), paramKd[i,j,k,nv,nw] * F[i,j,k,nw], Val(OP))
            end
        end
    end

    return nothing
end

# function apply_paramop_supdiag!(G::AbsArrNumber{K₊₁},  # output field; in 3D, G[i,j,k,w] is w-component of G at (i,j,k)
#                                 F::AbsArrNumber{K₊₁},  # input field; in 3D, F[i,j,k,w] is w-component of F at (i,j,k)
#                                 ::Val{OP},  # Val(:(=)) or Val(:(+=)): set (=) or add (+=) operator to use
#                                 paramKd::AbsArrComplexF{K₊₂},  # size of last two dimensions is Kf-by-Kf
#                                 kdiag::Integer;  # 0 ≤ kdiag ≤ Kf-1; index of diagonal of material parameter tensor to set (kdiag = 0: main diagonal)
#                                 α::Number=1.0  # scale factor to multiply to result before adding it to G: G += α param F
#                                 ) where {K₊₁,K₊₂,OP}
#     K = K₊₁ - 1
#     @assert K₊₂==K+2
#
#     Kf = size(paramKd, K+1)  # field dimension
#     @assert size(paramKd,K+2)==Kf
#
#     CI = CartesianIndices(size(G)[1:K])
#     for nv = 1:Kf  # row index of material parameter tensor
#         nw = mod1(nv+kdiag, Kf)  # column index of material parameter tensor
#         for ci = CI
#             @inbounds G[ci,nv] = paramKd[ci,nv,nw] * F[ci,nw]
#         end
#     end
#
#     return nothing
# end