text:
using ADSeismic
using ADCME
using PyPlot
using JLD2
using PyCall
using DelimitedFiles
using Statistics
mpi_init()

M = N = Int(round(sqrt(mpi_size())))
n_total = 1000
@info M
n = n_total÷M
param = MPIAcousticPropagatorParams(NX = n*M, NY = n*N, n = n,
     Rcoef=0.2, DELTAX=10, DELTAY=10, DELTAT=0.05, NSTEP = 100)
compute_PML_Params!(param)

srci = [n*param.M÷5]
srcj = [n*param.N÷2]
srcv = hcat([Ricker(param, 100.0, 500.0)]...)
src = MPIAcousticSource(param, srci,srcj,srcv)

rcvj = Array(20:n*param.M-19)
rcvi = ones(Int64, length(rcvj))*((n*param.N)÷5)
receiver =  MPIAcousticReceiver(param, rcvi, rcvj)

xs = Float64[]
ys = Float64[]
mask = ones(n, n)
for i = 1:n
    for j = 1:n
        x = i*param.DELTAX + (param.II-1) * n * param.DELTAX
        y = j*param.DELTAY + (param.JJ-1) * n * param.DELTAY
        idx = i + (param.II-1) * n
        idy = j + (param.JJ-1) * n
        if idx<=param.NPOINTS_PML || idx>=param.N*n-param.NPOINTS_PML ||
            idy<=param.NPOINTS_PML || idy>=param.M*n-param.NPOINTS_PML
            mask[i,j] = 0.0
        end
        push!(xs, x)
        push!(ys, y)
    end
end

init_guess = fc_init([2, 20, 20, 20, 1])
θ = placeholder(init_guess)
θ_ = mpi_bcast(θ)
c = fc([xs ys], [20, 20, 20, 1], θ_) * 1000.0
c = reshape(c, (n, n))
c = c .* mask + 1000.0 * (1 .- mask)

propagator = MPIAcousticPropagatorSolver(param, src, c)
MPISimulatedObservation!(propagator, receiver)
@load "data/dat$(mpi_rank())-$(mpi_size()).jld2" rcvv

local_loss = sum(propagator.u[end,:]) * 1e-20
if !isnothing(rcvv)
    @info mpi_rank(), "value"
    global local_loss += sum((receiver.rcvv - rcvv )^2)
end

loss = mpi_sum(local_loss)
G = gradients(loss, θ)

sess = Session(); init(sess)

@info "Runing loss..."
l = run(sess, loss)

init(sess)
@info "Runing grads..."
g = run(sess, G, θ=>init_guess)

timer_fwd = []
timer_bwd = []
for i = 1:10
    @info "forward, $i"
    init(sess)
    d = @timed run(sess, loss)
    push!(timer_fwd, d[2])
end

for i = 1:10
    @info "backward, $i"
    init(sess)
    d = @timed run(sess, G)
    push!(timer_bwd, d[2])
end

timer_fwd = mean(timer_fwd)
timer_bwd = mean(timer_bwd)

if mpi_rank()==0
    open("timing.txt", "a") do io
        writedlm(io, [mpi_size() timer_fwd timer_bwd])
    end
end

mpi_finalize()


text:
#TODO: write more tests

@testset "Testing SimpleBoreCollisionModule" begin

# check if delta is correctly calculated
function checkDeltaWithRandomPoints(cm, nPoints=1000)
  for i in 1:nPoints
    pos = [0u"mm", (rand()-0.5)*2*cm.params.scannerDiameter, (rand()-0.5)*2*cm.params.scannerDiameter]

    correct = checkDelta(cm, pos)
    if !correct
      return false
    end

  end
  return true
end

function checkDelta(cm, pos)
  state, delta = checkCoords(cm, pos, returnVerbose=true)
  if !all(state .== :VALID)
    state1, _ = checkCoords(cm, pos.-delta, returnVerbose=true)
    state2, _ = checkCoords(cm, pos.-(delta-sign.(delta).*[0,0.1,0.1]u"mm"), returnVerbose=true)
    if !( all(state1 .== :VALID) & all(state2[2:3] .== :INVALID) )
      @error pos, delta, state1, state2
      MPIMeasurements.plotSafetyErrors(cm, [1,2,3], transpose(hcat(pos, pos.-delta, pos.-(delta-sign.(delta).*[0,0.1,0.1]u"mm"))))
      return false
    end
  end

  return true
end

@testset "SimpleBoreCollisionModule: Circle" begin
  par = SimpleBoreCollisionModuleParams(scannerDiameter=118u"mm", objGeometry = Circle(name="Delta sample", diameter=10u"mm"), minMaxBoreAxis = [-300u"mm",Inf*u"mm"])
  cm = SimpleBoreCollisionModule(deviceID="testCollisionModule", params=par, dependencies=Dict{String, Union{Device, Missing}}())

  #@test checkDeltaWithRandomPoints(cm, 100000)
  @test checkCoords(cm, [0,0,0]u"mm")
  @test !checkCoords(cm, [0,1,0]u"m")
  @test !checkCoords(cm, [-300.1,0,0]u"mm")
  @test checkCoords(cm, [-300,0,0]u"mm")
  @test checkCoords(cm, [0u"mm", par.scannerDiameter/2 - par.clearance.distance - par.objGeometry.diameter/2, 0u"mm"])
  state, dist = checkCoords(cm, [0u"mm", par.scannerDiameter/2 - par.clearance.distance - par.objGeometry.diameter/2 + 0.001u"mm", 0u"mm"], returnVerbose=true)
  @test !all(state .== :VALID)
  @test isapprox(dist[2], 0.001u"mm", atol=1u"nm")

end

@testset "SimpleBoreCollisionModule: Rectangle" begin
  par = SimpleBoreCollisionModuleParams(scannerDiameter=118u"mm", objGeometry = Rectangle(name="Test rectangle", width=10u"mm", height=10u"mm"), minMaxBoreAxis = [-300u"mm",Inf*u"mm"])
  cm = SimpleBoreCollisionModule(deviceID="testCollisionModule", params=par, dependencies=Dict{String, Union{Device, Missing}}())

  #@test checkDeltaWithRandomPoints(cm, 100000)

end

@testset "SimpleBoreCollisionModule: Triangle" begin
  par = SimpleBoreCollisionModuleParams(scannerDiameter=118u"mm", objGeometry = Triangle(name="Test triangle", width=10u"mm", height=10u"mm"), minMaxBoreAxis = [-300u"mm",Inf*u"mm"])

  cm = SimpleBoreCollisionModule(deviceID="testCollisionModule", params=par, dependencies=Dict{String, Union{Device, Missing}}())

 #@test checkDeltaWithRandomPoints(cm, 1000) # this test fails

end

end

ext: jl


text:
function roundtrip_nonlinear_expression(
    expr, variable_to_string, string_to_variable
)
    node_list = MOF.OrderedObject[]
    object = MOF.convert_expr_to_mof(expr, node_list, variable_to_string)
    @test MOF.convert_mof_to_expr(object, node_list, string_to_variable) == expr
end

# hs071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
struct ExprEvaluator <: MOI.AbstractNLPEvaluator
    objective::Expr
    constraints::Vector{Expr}
end
MOI.features_available(::ExprEvaluator) = [:ExprGraph]
MOI.initialize(::ExprEvaluator, features) = nothing
MOI.objective_expr(evaluator::ExprEvaluator) = evaluator.objective
MOI.constraint_expr(evaluator::ExprEvaluator, i::Int) = evaluator.constraints[i]

function HS071(x::Vector{MOI.VariableIndex})
    x1, x2, x3, x4 = x
    return MOI.NLPBlockData(
        MOI.NLPBoundsPair.([25, 40], [Inf, 40]),
        ExprEvaluator(
            :(x[$x1] * x[$x4] * (x[$x1] + x[$x2] + x[$x3]) + x[$x3]),
            [
                :(x[$x1] * x[$x2] * x[$x3] * x[$x4] >= 25),
                :(x[$x1]^2 + x[$x2]^2 + x[$x3]^2 + x[$x4]^2 == 40)
            ]
        ),
        true
    )
end

@testset "Nonlinear functions" begin
    @testset "HS071 via MOI" begin
        model = MOF.Model()
        x = MOI.add_variables(model, 4)
        for (index, variable) in enumerate(x)
            MOI.set(model, MOI.VariableName(), variable, "var_$(index)")
        end
        MOI.add_constraints(model, MOI.SingleVariable.(x),
                            Ref(MOI.Interval(1.0, 5.0)))
        MOI.set(model, MOI.NLPBlock(), HS071(x))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        MOI.write_to_file(model, TEST_MOF_FILE)
        @test replace(read(TEST_MOF_FILE, String), '\r' => "") ==
            replace(read(joinpath(@__DIR__, "nlp.mof.json"), String), '\r' => "")
        MOF.validate(TEST_MOF_FILE)
    end
    @testset "Error handling" begin
        node_list = MOF.Object[]
        string_to_variable = Dict{String, MOI.VariableIndex}()
        variable_to_string = Dict{MOI.VariableIndex, String}()
        # Test unsupported function for Expr -> MOF.
        @test_throws Exception MOF.convert_expr_to_mof(
            :(not_supported_function(x)), node_list, variable_to_string)
        # Test unsupported function for MOF -> Expr.
        @test_throws Exception MOF.convert_mof_to_expr(
            MOF.OrderedObject("type"=>"not_supported_function", "value"=>1),
            node_list, string_to_variable)
        # Test n-ary function with no arguments.
        @test_throws Exception MOF.convert_expr_to_mof(
            :(min()), node_list, variable_to_string)
        # Test unary function with two arguments.
        @test_throws Exception MOF.convert_expr_to_mof(
            :(sin(x, y)), node_list, variable_to_string)
        # Test binary function with one arguments.
        @test_throws Exception MOF.convert_expr_to_mof(
            :(^(x)), node_list, variable_to_string)
        # An expression with something other than :call as the head.
        @test_throws Exception MOF.convert_expr_to_mof(
            :(a <= b <= c), node_list, variable_to_string)
        # Hit the default fallback with an un-interpolated complex number.
        @test_throws Exception MOF.convert_expr_to_mof(
            :(1 + 2im), node_list, variable_to_string)
        # Invalid number of variables.
        @test_throws Exception MOF.substitute_variables(
            :(x[1] * x[2]), [MOI.VariableIndex(1)])
        # Function-in-Set
        @test_throws Exception MOF.extract_function_and_set(
            :(foo in set))
        # Not a constraint.
        @test_throws Exception MOF.extract_function_and_set(:(x^2))
        # Two-sided constraints
        @test MOF.extract_function_and_set(:(1 <= x <= 2)) ==
            MOF.extract_function_and_set(:(2 >= x >= 1)) ==
            (:x, MOI.Interval(1, 2))
        # Less-than constraint.
        @test MOF.extract_function_and_set(:(x <= 2)) ==
            (:x, MOI.LessThan(2))
    end
    @testset "Roundtrip nonlinear expressions" begin
        x = MOI.VariableIndex(123)
        y = MOI.VariableIndex(456)
        z = MOI.VariableIndex(789)
        string_to_var = Dict{String, MOI.VariableIndex}("x"=>x, "y"=>y, "z"=>z)
        var_to_string = Dict{MOI.VariableIndex, String}(x=>"x", y=>"y", z=>"z")
        for expr in [2, 2.34, 2 + 3im, x, :(1 + $x), :($x - 1),
                     :($x + $y), :($x + $y - $z), :(2 * $x), :($x * $y),
                     :($x / 2), :(2 / $x), :($x / $y), :($x / $y / $z), :(2^$x),
                     :($x^2), :($x^$y), :($x^(2 * $y + 1)), :(sin($x)),
                     :(sin($x + $y)), :(2 * $x + sin($x)^2 + $y),
                     :(sin($(3im))^2 + cos($(3im))^2), :($(1 + 2im) * $x),
                     :(ceil($x)), :(floor($x)), :($x < $y), :($x <= $y),
                     :($x > $y), :($x >= $y), :($x == $y), :($x != $y),
                     # :($x && $y), :($x || $y),
                     :(ifelse($x > 0, 1, $y))]
            roundtrip_nonlinear_expression(expr, var_to_string, string_to_var)
        end
    end
    @testset "Reading and Writing" begin
        # Write to file.
        model = MOF.Model()
        (x, y) = MOI.add_variables(model, 2)
        MOI.set(model, MOI.VariableName(), x, "var_x")
        MOI.set(model, MOI.VariableName(), y, "y")
        con = MOI.add_constraint(model,
                 MOF.Nonlinear(:(2 * $x + sin($x)^2 - $y)),
                 MOI.EqualTo(1.0))
        MOI.set(model, MOI.ConstraintName(), con, "con")
        MOI.write_to_file(model, TEST_MOF_FILE)
        # Read the model back in.
        model2 = MOF.Model()
        MOI.read_from_file(model2, TEST_MOF_FILE)
        con2 = MOI.get(model2, MOI.ConstraintIndex, "con")
        foo2 = MOI.get(model2, MOI.ConstraintFunction(), con2)
        # Test that we recover the constraint.
        @test foo2.expr == :(2 * $x + sin($x)^2 - $y)
        @test MOI.get(model, MOI.ConstraintSet(), con) ==
                MOI.get(model2, MOI.ConstraintSet(), con2)
        MOF.validate(TEST_MOF_FILE)
    end
end


text:
# Record view
# -----------

export HTSRecord

"""
A light-weight wrapper of records.
"""
mutable struct HTSRecord
    ptr::Ptr{htslib.bam1_t}
    function HTSRecord()
        ptr = htslib.bam_init1()
        @assert ptr != C_NULL
        record = new(ptr)
        Base.finalizer(record) do record
            ptr = getfield(record, :ptr)
            ptr == C_NULL && error("pointer is invalid")
            htslib.bam_destroy1(ptr)
        end
        record
    end
end

function HTSRecord(x::HTSRecord)
    copy(x)
end

function Base.copy!(dst::HTSRecord, src::HTSRecord)
    GC.@preserve src begin
        p = htslib.bam_copy1(pointer(dst), pointer(src))
        p == C_NULL && error("error copying bam record $src")
    end
    dst
end

function Base.copy(record::HTSRecord)
    ans = HTSRecord()
    copy!(ans, record)
    ans
end

function Base.pointer(x::HTSRecord)
    x.ptr
end

function Base.show(io::IO, record::HTSRecord)
    GC.@preserve record begin
        ptr = record.ptr
        ptr == C_NULL && error("pointer is invalid")
        bam = unsafe_load(ptr)
        @printf(io, "%s(<%p>) %d:%d-%d ", summary(record), ptr, refid(record), leftposition(record), rightposition(record))
    end
    print(io, "\n")
    print(io, "flag     : $(summary(record.flag))\n")
    print(io, "queryname: $(queryname(record))\n")
    print(io, "mapqual  : $(mappingquality(record))\n")
    print(io, "sequence : ")
    for s in sequence(record)
        print(io, s)
    end
    print(io, "\n")
    print(io, "quality  : ")
    for q in quality(record)
        print(io, quality_char(q))
    end
    print(io, "\n")
    print(io, "matepos  : $(next_refid(record)) : $(next_refpos(record))\n")
    record
end

function quality_char(qual::UInt8)
    # taken from jakobnissen
	if qual < 10
		return  ' '
	end
	if qual < 15
		return '▁'
	end
	if qual < 20
		return '▂'
	end
	if qual < 25
		return '▃'
	end
	if qual < 30
		return '▄'
	end
	if qual < 35
		return '▆'
	end
	if qual < 40
		return '▇'
	end
	if qual < 255
		return '█'
	end
	return '?'
end

@inline function Base.getproperty(view::HTSRecord, name::Symbol)
    GC.@preserve view begin
        bam = unsafe_load(getfield(view, :ptr))
    end
    if name == :pos
        return bam.core.pos
    elseif name == :tid
        return bam.core.tid
    elseif name == :bin
        return bam.core.bin
    elseif name == :qual
        return bam.core.qual
    elseif name == :l_extranul
        return bam.core.l_extranul
    elseif name == :flag
        return HTSRecordFlag(bam.core.flag)
    elseif name == :l_qname
        return bam.core.l_qname
    elseif name == :n_cigar
        return bam.core.n_cigar
    elseif name == :l_qseq
        return bam.core.l_qseq
    elseif name == :mtid
        return bam.core.mtid
    elseif name == :mpos
        return bam.core.mpos
    elseif name == :isize
        return bam.core.isize
    else
        return getfield(view, name)
    end
end

function Base.propertynames(record::HTSRecord)
    return (
        "pos",
        "tid",
        "bin",
        "qual",
        "l_extranul",
        "flag",
        "l_qname",
        "n_cigar",
        "l_qseq",
        "mtid",
        "mpos",
        "isize",
    )
end

##### Accessor functions

const seq_nt16_str = ('=','A','C','M','G','R','S','V','T','W','Y','H','K','D','B','N')

function BioGenerics.leftposition(record::HTSRecord)::Int64
    record.pos + 1 # start pos is zero-based
end

function BioGenerics.rightposition(record::HTSRecord)::Int64
    GC.@preserve record begin
        htslib.bam_endpos(pointer(record))
    end
end

function seqlength(record::HTSRecord)::Int32
    record.l_qseq
end

function BioGenerics.sequence(record::HTSRecord, i::Integer)::Char
    GC.@preserve record begin
        seq_ptr = htslib.bam_get_seq(pointer(record))
        idx = htslib.bam_seqi(seq_ptr, i-UInt8(1)) + 1
        seq_nt16_str[idx]
    end
end

function sequence!(record::HTSRecord, array::AbstractVector{Char})
    seql = seqlength(record)
    empty!(array)
    sizehint!(array, seql)
    GC.@preserve record begin
        seq_ptr = htslib.bam_get_seq(pointer(record))
        for i in 1:seql
            push!(array, seq_nt16_str[htslib.bam_seqi(seq_ptr, i-UInt8(1)) + 1])
        end
    end
    array
end

function BioGenerics.sequence(record::HTSRecord)
    ans = Vector{Char}(undef, seqlength(record))
    sequence!(record, ans)
end

function mappingquality(record::HTSRecord)::UInt8
    record.qual
end

function refid(record::HTSRecord)::Int32
    record.tid + Int32(1)
end

# TODO
function refname end

function quality!(record::HTSRecord, array::AbstractVector)
    seql = seqlength(record)
    empty!(array)
    sizehint!(array, seql)
    GC.@preserve record begin
        qual_ptr = htslib.bam_get_qual(pointer(record))
        for i in 1:seql
            push!(array, unsafe_load(qual_ptr, i))
        end
    end
    array
end

function quality(record::HTSRecord)
    #### TODO: quality score 255 represents missing. Do we want to return missing?
    ans = Vector{UInt8}(undef, seqlength(record))
    quality!(record, ans)
end

function setquality!(record::HTSRecord, i::Integer, value::Integer)
    #### TODO: quality score 255 represents missing.?
    value = convert(UInt8, value)
    GC.@preserve record begin
        ## bounds check
        if i <= 0 || i > seqlength(record)
            throw(BoundsError(record, i))
        end
        qual_ptr = htslib.bam_get_qual(pointer(record))
        unsafe_store!(qual_ptr, value, i)
    end
    record
end

function setquality!(record::HTSRecord, value::AbstractVector{<:Integer})
    seqlength(record) != length(value) && throw(BoundsError(value))
    for (i, v) in enumerate(value)
        setquality!(record, i, v)
    end
    record
end

function setsequence!(record::HTSRecord, i::Integer, value::Char)
    GC.@preserve record begin
        ## bounds check
        if i <= 0 || i > seqlength(record)
            throw(BoundsError(record, i))
        end
        idx = findfirst(x -> x == value, seq_nt16_str)
        isnothing(idx) && throw(DomainError(value))
        base_code = UInt8(idx - 1)
        htslib.bam_set_seqi(htslib.bam_get_seq(pointer(record)), i-1, base_code)
    end
    record
end

function queryname(record::HTSRecord)::String
    # Note: l_qname is length of the name _plus_ number of NULs (1-4) for memory alignment
    GC.@preserve record begin
        l_qname = record.l_qname
        l_qname <= 0 && return "*" # If set to 0, the placeholder query name "*" will be used.
        cstr = htslib.bam_get_qname(pointer(record))
        unsafe_string(pointer(cstr))
    end
end

function next_refid(record::HTSRecord)::Union{Int32,Missing}
    ans = record.mtid + Int32(1)
    ans == 0 && return missing
    ans
end

function next_refpos(record::HTSRecord)::Union{Int32,Missing}
    ans = record.mpos + Int32(1)
    ans == 0 && return missing
    ans
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 4b75cde6-db35-11eb-1b95-bba4e05edb2f
begin
    using VoronoiFVM
    using ExtendableGrids
    using GridVisualize
	using PlutoVista
	GridVisualize.default_plotter!(nothing)
    using PlutoUI
	using LinearAlgebra
    if isdefined(Main,:PlutoRunner)
        GridVisualize.default_plotter!(PlutoVista)
    end
end

# ╔═╡ 3fa189e4-9e1c-470c-bf26-15b631945d2d
md"""
# Interfaces conditions and APIs

## Two subdomains, species number based API
For a simple stationary diffusion equation with an interior interface, we discuss possible interface conditions between two subdomains.

Let ``\Omega=\Omega_1\cup\Omega_2`` where ``\Omega_1=(-1,0)`` and ``\Omega_2=(0,1)``.
Let ``\Gamma_1={-1}``,``\Gamma_2={1}``  and ``\Gamma_3={0}``.

Regard the following problem:

``\begin{aligned}
     -\Delta u_1 &= 0 & \text{in}\quad \Omega_1\\
     -\Delta u_2 &= 0 & \text{in}\quad \Omega_2\\
\end{aligned}
``

with exterior boundary conditions

``u_1|_{\Gamma_1} = g_1`` and ``u_2|_{\Gamma_2} = g_2``

For the interior boundary (interface) conditons we set

``\nabla u_1|_{\Gamma_3}+f_1(u_1,u_2)=0``

``-\nabla u_2|_{\Gamma_3}+f_2(u_1,u_2)=0``

where ``f_1``, ``f_2`` are discussed later.
"""

# ╔═╡ 327af4a8-74cc-4834-ab19-d3a1d6873982
TableOfContents(title="")

# ╔═╡ d5a0ee0d-959d-476b-b3c5-79b741059992
md"""
### Set up
"""

# ╔═╡ f03ff283-c989-4b1a-b73e-2e616054e3db
md"""
Create a grid with two subdomins and an interface in the center.
"""

# ╔═╡ 670c78c1-d0be-4362-975b-2c944620681f
nref=2

# ╔═╡ 79193d53-9dfa-47b7-aed2-c7eb43769b5f
begin
	hmax=0.2/2.0^nref
    hmin=0.05/2.0^nref
    X1=geomspace(-1.0,0.0, hmax,hmin)
    X2=geomspace(0.0,1.0, hmin,hmax)
    X=glue(X1,X2)
    grid=VoronoiFVM.Grid(X)

    bfacemask!(grid, [0.0],[0.0],3)
    ## Material 1 left of 0
    cellmask!(grid, [-1.0],[0.0],1)
    ## Material 2 right of 0
    cellmask!(grid, [0.0],[1.0],2)
end;

# ╔═╡ 4cb07222-587b-4a74-a444-43f5433d5b03
gridplot(grid,legend=:rt,resolution=(600,200))

# ╔═╡ 02ec3c0b-6e68-462b-84df-931370cbdcac
md"""
For later use (plotting) extract the two subgrids from the grid
"""

# ╔═╡ 592429a1-108c-4e9e-8961-497c2c31f319
subgrid1=subgrid(grid,[1]);

# ╔═╡ 1fa5d9b3-0558-4558-a5f1-f34f70c8d9a0
subgrid2=subgrid(grid,[2]);

# ╔═╡ 5b539ad2-46d6-43b0-8b3e-f8a7e1ae0a6d
md"""
Define the diffusion flux
"""

# ╔═╡ 6aabfbe1-de7d-49ba-8144-6d364b21b34f
    function flux!(f,u,edge)
        if edge.region==1
            f[1]=u[1,1]-u[1,2]
        end
        if edge.region==2
            f[2]=u[2,1]-u[2,2]
        end
    end

# ╔═╡ dbd27f10-9fd9-450d-9f78-89ea738d605b
md"""
Specify  the outer boundary values.
"""

# ╔═╡ e5b432b2-875e-49a5-8e78-68f7f47f06c3
g_1=1.0

# ╔═╡ 8ae32292-df3e-4190-83a5-ec5ba529299e
g_2=0.1

# ╔═╡ 677853a7-0c43-4d3c-bc74-799535f95aeb
md"""
Create the system. We pass the interface condition function as a parameter.
"""

# ╔═╡ 139058a9-44a0-43d5-a377-4fc72927fa28
function make_system(breaction)
	physics=VoronoiFVM.Physics(
        flux=flux!,
        breaction=breaction
    )

    ## Create system
    sys=VoronoiFVM.System(grid,physics,unknown_storage=:sparse)

    ##  put potential into both regions
    enable_species!(sys,1,[1])
    enable_species!(sys,2,[2])

    ## Set boundary conditions
    for ispec=1:2
        boundary_dirichlet!(sys,ispec,1,g_1)
        boundary_dirichlet!(sys,ispec,2,g_2)
    end
	sys
end

# ╔═╡ de251b03-dd3a-4a44-9440-b7e654c32dac
md"""
Stationary solution with zero initial value
"""

# ╔═╡ 5b3a49d6-30ce-4d6e-94b4-a76645d7d8ce
function mysolve(sys)
	U=solve(unknowns(sys,inival=0),sys)
	U1=view(U[1,:],subgrid1)
    U2=view(U[2,:],subgrid2)
	U1,U2
end

# ╔═╡ 0735c061-68e1-429f-80f1-d8410989a91d
md"""
Plot the results
"""

# ╔═╡ 467dc381-3b3d-4de7-a7f9-bfc51300832b
function plot(U1,U2;title="")
	vis=GridVisualizer(resolution=(600,300))
    scalarplot!(vis,subgrid1,U1,clear=false,show=false, color=:green, label="u1")
    scalarplot!(vis,subgrid2,U2,clear=false,show=true, color=:blue, label="u2",legend=:rt,title=title,flimits=(-0.5,1.5))
end

# ╔═╡ fa4fcc0d-1d3a-45a2-8857-50536bbe39cc
md"""
### No interface reaction

This means we set ``f_1(u_1,u_2)=0`` and ``f_2(u_1,u_2)=0``.
"""

# ╔═╡ 8f210696-fcf4-47bc-a5a2-c561ad7efcbd
function noreaction(f,u,node) end

# ╔═╡ 57e8515e-3be1-4478-af98-430501438ee7
system1=make_system(noreaction);

# ╔═╡ 56136cd1-0c01-449d-9297-68924ac99ee7
plot(mysolve(system1)...)

# ╔═╡ aad305a9-aac6-4aff-9f8e-08d6a2f756c8
md"""
### Mass action law reaction ``u_1 \leftrightharpoons u_2``

This is a rather general ansatz where we assume a backward-forward reaction at the interface with reaction constants ``k_1`` and ``k_2``, respectively.

According to the mass action law, this tranlates to a reaction rate

``r(u_1,u_2)=k_1u_1 - k_2u_2``

and correspondingly

``f_1(u_1,u_2)=r``

``f_2(u_1,u_2)=-r``

Note, that ``f_i`` is monotonically increasing in ``u_i`` and monotonically decreasing in the respective other argument, leading to an M-Property of the overall discretization matrix.

Note that the "no reaction" case is just a special case where ``k_1,k_2=0``.
"""

# ╔═╡ d027ff24-3ad1-4528-b5cf-10814caf30db
k1=0; k2=10;

# ╔═╡ 1328b4bf-2d64-4b02-a910-1995da8be28b
function mal_reaction(f,u,node)
        if node.region==3
            react=k1*u[1]-k2*u[2]
            f[1]= react
            f[2]= -react

        end
end

# ╔═╡ 610a0761-1c23-415d-a187-f7d93a1b7637
system2=make_system(mal_reaction)

# ╔═╡ 87edce1f-df6d-4cd8-bce5-24fb666cd6b5
begin
	k1,k2
	U1,U2=mysolve(system2)
    plot(U1,U2;title="k1=$(k1), k2=$(k2)")
end

# ╔═╡ b82fc6b2-eee1-4a91-a115-61b86621f686
md"""
### Penalty enforcing continuity

Setting ``k_1,k_2`` to a large number leads to another special case of the above reaction - similar to the penalty method to implement the Dirichlet boundary conditions, this lets the reaction equation dominate, which in this case forces
``u_1-u_2=0`` at the interface, and thus continuity.
"""

# ╔═╡ 9eaea813-2628-47d0-9d36-54c367689142
function penalty_reaction(f,u,node)
        if node.region==3
            react=1.0e10*(u[1]-u[2])
            f[1]= react
            f[2]= -react
        end
end

# ╔═╡ 817738c0-f1a3-4779-9075-7ea051a81e73
system3=make_system(penalty_reaction);

# ╔═╡ 80019ef1-bf41-4a55-9262-613a2d20be1f
plot(mysolve(system3)...)

# ╔═╡ 31e00855-8906-4be5-8e69-e2d8d9539e04
md"""
### Interface recombination

Here, we implement an annihilation reaction ``u_1 + u_2 \to \emptyset``
Acoording to the mass action law, this is implemented via

``r(u_1,u_2)=k_r u_1 u_2``

``f_1(u_1,u_2)=r``

``f_2(u_1,u_2)=r``

"""

# ╔═╡ f2490f99-04ca-4f42-af2a-53adae51ca68
k_r=1000

# ╔═╡ 39a0db1b-3a4e-4108-b43f-d4e578c92608
function recombination(f,u,node)
        if node.region==3
            react=k_r*(u[1]*u[2])
            f[1]= react
            f[2]= react
        end
end;

# ╔═╡ 644149fb-2264-42bd-92c9-193ab07c08f6
system4=make_system(recombination);

# ╔═╡ b479402f-ef00-4425-8b0f-45f2dae74d80
plot(mysolve(system4)...)

# ╔═╡ ed068b51-92af-48d5-9230-debc178ec827
md"""
### Thin  conductive interface layer

Let us assume that the interface is of thickness $d$ which is however small with respect to ``\Omega`` that we want to derive an interface condition from the assumption of an exact continuous solution within the interface.

So let ``\Omega_I=(x_l,x_r)`` be  the interface region where
we have ``-\Delta u_I=0`` with values ``u_l``, ``u_r`` at the boundaries.

Then we have for the flux  in the interface region, ``q_I=\nabla u = \frac1{d}(u_r - u_l)``

Continuity of fluxes then gives ``f_1=q_I`` and ``f_2=-q_I``.

Continuity of ``u`` gives ``u_{1,I}=u_l, u_{2,I}=u_r``
This gives

``r=q_I=\frac{1}{d}(u_1-u_{2})``

``f_1(u_1,v_1)=r``

``f_2(u_1,v_1)=-r``

and therefore another special case of the mass action law condition.
"""

# ╔═╡ a2d919a5-a395-40fb-8f93-db742f8a77c2
d=1

# ╔═╡ 58d8831b-ad66-4f77-a33a-933c15c46a52
function thinlayer(f,u,node)
        if node.region==3
            react=(u[1]-u[2])/d
            f[1]= react
            f[2]= -react
        end
end

# ╔═╡ 8c0b4ab5-09da-4d8f-b001-5e15f823423c
system5=make_system(thinlayer);

# ╔═╡ d3d99b9c-ad18-4a04-bb3e-f17dd542f9f3
plot(mysolve(system5)...)

# ╔═╡ eb9abf2e-372c-4f79-afbe-772b90eff9ad
md"""
## Quantity based API

From this discussion it seems that discontinuous interface conditions can be formulated in a rather general way via linear or nonlinear robin boundary conditions for each of the adjacent discontinuous species. Technically, it is necessary to be able to access the adjacent bulk data.
"""

# ╔═╡ 9f3ae7b5-51b3-48bc-b4db-b7236ba30682
md"""
Here, we propose an API layer on top  of the species handling of VoronoiFVM.
For a start, we call these "Meta species" "quantities". There may be a different name, however.
"""

# ╔═╡ 2da8a5b1-b168-41d9-baa8-d65a4ef5c4c0
md"""
We define a grid with N subregions
"""

# ╔═╡ d44407de-8c9c-42fa-b1a2-ae02b826eccc
N=6

# ╔═╡ ae268316-c058-4db8-9b71-57b0d9425274
begin
	XX=collect(0:0.1:1)
	local xcoord=XX
	for i=1:N-1
		xcoord=glue(xcoord,XX.+i)
	end
	grid2=simplexgrid(xcoord)
	for i=1:N
		cellmask!(grid2,[i-1],[i],i)
	end
	for i=1:N-1
		bfacemask!(grid2,[i],[i],i+2)
	end
end

# ╔═╡ b53b9d28-4c25-4fb8-a3e4-599b0e385121
gridplot(grid2,legend=:lt,resolution=(600,200))

# ╔═╡ e7ce7fd4-cfa4-4cc6-84a2-7e20ed2f4e5c
md"""
To work with quantities, we first introduce a new constructor call without the "physics" parameter:
"""

# ╔═╡ 29f36902-e355-4b02-b7b0-c4db12c47d33
system6=VoronoiFVM.System(grid2)

# ╔═╡ 673e9320-ea30-4831-ad85-ba7936293ee2
md"""
First, we introduce a continuous quantity which we name "cspec". Note that the "species number" can be assigned automatically if not given explicitely.
"""

# ╔═╡ f35f419a-94dd-4051-a533-4b1ec9a4c9ec
cspec=ContinuousQuantity(system6,1:N,ispec=1)

# ╔═╡ 9661e4fc-55e1-4c2c-a3ad-515cdac3b514
md"""
A discontinuous quantity can be introduced as well. by default, each reagion gets a new species number. This can be overwritten by the user.
"""

# ╔═╡ 90298676-fda7-4168-8a40-7ff53e7c761b
dspec=DiscontinuousQuantity(system6,1:N; regionspec=[2+i%2 for i=1:N])

# ╔═╡ 7a819522-55e4-4547-a3e0-e047e74cfb6b
system6

# ╔═╡ cebabf33-e769-47bd-b6f1-ddf525fea895
md"""
For both quantities, we define simple diffusion fluxes:
"""

# ╔═╡ 719f206a-5b9f-4d78-8778-1d89edb2bc4d
	function flux2(f,u,edge)
		f[dspec]=u[dspec,1]-u[dspec,2]
		f[cspec]=u[cspec,1]-u[cspec,2]
	end

# ╔═╡ 1d7f442f-c057-4379-8a40-c6ce3646ad5c
md"""
Define a thin layer inteface condition for `dspec` and an interface source for `cspec`.
"""

# ╔═╡ da41b22e-114d-4eee-81d0-73e6f3b45242
md"""
Add physics to the system, set dirichlet bc at both ends, and extract subgrids
for plotting (until there will be a plotting API for this...)
"""

# ╔═╡ 0c3ca093-e011-409c-be58-e617e02f5d4c
function plot2(U,subgrids,system6)
	dvws=VoronoiFVM.views(U,dspec,subgrids,system6)
	cvws=VoronoiFVM.views(U,cspec,subgrids,system6)
	vis=GridVisualizer(resolution=(600,300))
	for i=1:length(dvws)
		scalarplot!(vis,subgrids[i],dvws[i],flimits=(-0.5,1.5),clear=false,color=:red)
		scalarplot!(vis,subgrids[i],cvws[i],flimits=(-0.5,1.5),clear=false,color=:green)
 	end
	reveal(vis)
end

# ╔═╡ b8cd6ad1-d323-4888-bbd1-5deba5a5870d
d1=0.1

# ╔═╡ 441a39a0-a7de-47db-8539-12dee30b8312
q1=0.2

# ╔═╡ d6e1c6c7-060d-4c2f-8054-d8f33f54bd55
function breaction2(f,u,node)
	       if node.region>2
            react=(u[dspec,1]-u[dspec,2])/d1
            f[dspec,1]= react
            f[dspec,2]= -react

		    f[cspec]=-q1
        end
end

# ╔═╡ 59c83a22-a4cc-4b51-a1cc-5eb39588eacd
begin
	physics!(system6,VoronoiFVM.Physics(
        flux=flux2,
        breaction=breaction2
    ))

    ## Set boundary conditions
    boundary_dirichlet!(system6,dspec,1,g_1)
    boundary_dirichlet!(system6,dspec,2,g_2)
    boundary_dirichlet!(system6,cspec,1,0)
    boundary_dirichlet!(system6,cspec,2,0)
	subgrids=VoronoiFVM.subgrids(dspec,system6)
end;

# ╔═╡ de119a22-b695-4b4f-8e04-b7d68ec1e91b
subgrids;sol6=solve(unknowns(system6,inival=0),system6)

# ╔═╡ 5d8aca85-f12d-4a6e-84e4-781d45b65742
plot2(sol6,subgrids,system6)

# ╔═╡ 6203d95b-13d2-48a6-a69f-39e33e2edcdb
md"""

## Open problems with current implementation
- Testfunctions: this will be a another problem to be handeled
- Alternative is "glueing together" systems
- Forwarddiff always sees species all species, leading to some overhead
"""

# ╔═╡ 4d8e81c1-dbec-4379-9ab7-a585a369582d
if d1==0.1
@assert norm(system6,sol6,2)==7.0215437706445245
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PlutoVista = "646e1f28-b900-46d7-9d87-d554eb38a413"
VoronoiFVM = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"

[compat]
ExtendableGrids = "~0.8.3"
GridVisualize = "~0.3.4"
PlutoUI = "~0.7.16"
PlutoVista = "~0.8.1"
VoronoiFVM = "~0.13.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc1"
manifest_format = "2.0"

[[deps.AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "b8d49c34c3da35f220e7295659cd0bab8e739fed"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.33"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Cassette]]
git-tree-sha1 = "6ce3cd755d4130d43bab24ea5181e77b89b51839"
uuid = "7057c7e9-c182-5462-911a-8362d720325c"
version = "0.3.9"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "74e8234fb738c45e2af55fdbcd9bfbe00c2446fa"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.8.0"

[[deps.ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[deps.CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[deps.DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "7220bc21c33e990c14f4a9a319b1d242ebc5b269"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.1"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[deps.ElasticArrays]]
deps = ["Adapt"]
git-tree-sha1 = "a0fcc1bb3c9ceaf07e1d0529c9806ce94be6adf9"
uuid = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
version = "1.2.9"

[[deps.ExtendableGrids]]
deps = ["AbstractTrees", "Dates", "DocStringExtensions", "ElasticArrays", "InteractiveUtils", "LinearAlgebra", "Printf", "Random", "SparseArrays", "Test"]
git-tree-sha1 = "85dbe70afc7153ad510577f158214b45be22593b"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
version = "0.8.3"

[[deps.ExtendableSparse]]
deps = ["DocStringExtensions", "LinearAlgebra", "Printf", "Requires", "SparseArrays", "SuiteSparse", "Test"]
git-tree-sha1 = "0341e41e45e6c5e46be89c33543082a0a867707c"
uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
version = "0.6.5"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "3c041d2ac0a52a12a27af2782b34900d9c3ee68c"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.7"

[[deps.FiniteDiff]]
deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "8b3c09b56acaf3c0e581c66638b85c8650ee9dca"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.8.1"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "c4203b60d37059462af370c4f3108fb5d155ff13"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.20"

[[deps.GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[deps.GridVisualize]]
deps = ["ColorSchemes", "Colors", "DocStringExtensions", "ElasticArrays", "ExtendableGrids", "GeometryBasics", "LinearAlgebra", "Observables", "OrderedCollections", "PkgVersion", "Printf", "Requires", "StaticArrays"]
git-tree-sha1 = "0cf1ee93e7b6e88eb0b7e5122e72c7d88102bbc5"
uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
version = "0.3.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[deps.Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[deps.IterativeSolvers]]
deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
git-tree-sha1 = "1a8c6237e78b714e901e406c096fc8a65528af7d"
uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
version = "0.9.1"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLD2]]
deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "46b7834ec8165c541b0b5d1c8ba63ec940723ffb"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.15"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Observables]]
git-tree-sha1 = "fe29afdef3d0c4a8286128d4e45cc50621b1e43d"
uuid = "510215fc-4207-5dde-b226-833fc4488ee2"
version = "0.4.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "98f59ff3639b3d9485a03a72f3ab35bab9465720"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.6"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PkgVersion]]
deps = ["Pkg"]
git-tree-sha1 = "a7a7e1a88853564e551e4eba8650f8c38df79b37"
uuid = "eebad327-c553-4316-9ea0-9fa01ccd7688"
version = "0.1.1"

[[deps.PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[deps.PlutoVista]]
deps = ["ColorSchemes", "Colors", "GridVisualize", "UUIDs"]
git-tree-sha1 = "05c4f4f1a94ff2cfb6d442cced5965c8ab422e92"
uuid = "646e1f28-b900-46d7-9d87-d554eb38a413"
version = "0.8.1"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

[[deps.RecursiveArrayTools]]
deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "FillArrays", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "ff7495c78a192ff7d59531d9f14db300c847a4bc"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.19.1"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SparseDiffTools]]
deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"]
git-tree-sha1 = "36a4d27a02af48a1eafd2baff58b32deeeb68926"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "1.16.5"

[[deps.SparsityDetection]]
deps = ["Cassette", "DocStringExtensions", "LinearAlgebra", "SparseArrays", "SpecialFunctions"]
git-tree-sha1 = "9e182a311d169cb9fe0c6501aa252983215fe692"
uuid = "684fba80-ace3-11e9-3d08-3bc7ed6f96df"
version = "0.3.4"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "793793f1df98e3d7d554b65a107e9c9a6399a6ed"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.7.0"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.3.3"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.VertexSafeGraphs]]
deps = ["LightGraphs"]
git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e"
uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
version = "0.1.2"

[[deps.VoronoiFVM]]
deps = ["DiffResults", "DocStringExtensions", "ExtendableGrids", "ExtendableSparse", "ForwardDiff", "GridVisualize", "IterativeSolvers", "JLD2", "LinearAlgebra", "Printf", "RecursiveArrayTools", "SparseArrays", "SparseDiffTools", "SparsityDetection", "StaticArrays", "SuiteSparse", "Test"]
git-tree-sha1 = "db381645bec1d798c560fe415593663d4a2194ac"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
version = "0.13.1"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─3fa189e4-9e1c-470c-bf26-15b631945d2d
# ╟─327af4a8-74cc-4834-ab19-d3a1d6873982
# ╟─d5a0ee0d-959d-476b-b3c5-79b741059992
# ╟─f03ff283-c989-4b1a-b73e-2e616054e3db
# ╠═670c78c1-d0be-4362-975b-2c944620681f
# ╠═79193d53-9dfa-47b7-aed2-c7eb43769b5f
# ╠═4cb07222-587b-4a74-a444-43f5433d5b03
# ╟─02ec3c0b-6e68-462b-84df-931370cbdcac
# ╠═592429a1-108c-4e9e-8961-497c2c31f319
# ╠═1fa5d9b3-0558-4558-a5f1-f34f70c8d9a0
# ╟─5b539ad2-46d6-43b0-8b3e-f8a7e1ae0a6d
# ╠═6aabfbe1-de7d-49ba-8144-6d364b21b34f
# ╟─dbd27f10-9fd9-450d-9f78-89ea738d605b
# ╠═e5b432b2-875e-49a5-8e78-68f7f47f06c3
# ╠═8ae32292-df3e-4190-83a5-ec5ba529299e
# ╟─677853a7-0c43-4d3c-bc74-799535f95aeb
# ╠═139058a9-44a0-43d5-a377-4fc72927fa28
# ╟─de251b03-dd3a-4a44-9440-b7e654c32dac
# ╠═5b3a49d6-30ce-4d6e-94b4-a76645d7d8ce
# ╟─0735c061-68e1-429f-80f1-d8410989a91d
# ╠═467dc381-3b3d-4de7-a7f9-bfc51300832b
# ╟─fa4fcc0d-1d3a-45a2-8857-50536bbe39cc
# ╠═8f210696-fcf4-47bc-a5a2-c561ad7efcbd
# ╠═57e8515e-3be1-4478-af98-430501438ee7
# ╠═56136cd1-0c01-449d-9297-68924ac99ee7
# ╟─aad305a9-aac6-4aff-9f8e-08d6a2f756c8
# ╠═1328b4bf-2d64-4b02-a910-1995da8be28b
# ╠═610a0761-1c23-415d-a187-f7d93a1b7637
# ╠═d027ff24-3ad1-4528-b5cf-10814caf30db
# ╟─87edce1f-df6d-4cd8-bce5-24fb666cd6b5
# ╟─b82fc6b2-eee1-4a91-a115-61b86621f686
# ╠═9eaea813-2628-47d0-9d36-54c367689142
# ╠═817738c0-f1a3-4779-9075-7ea051a81e73
# ╠═80019ef1-bf41-4a55-9262-613a2d20be1f
# ╟─31e00855-8906-4be5-8e69-e2d8d9539e04
# ╠═39a0db1b-3a4e-4108-b43f-d4e578c92608
# ╠═644149fb-2264-42bd-92c9-193ab07c08f6
# ╠═f2490f99-04ca-4f42-af2a-53adae51ca68
# ╠═b479402f-ef00-4425-8b0f-45f2dae74d80
# ╟─ed068b51-92af-48d5-9230-debc178ec827
# ╠═58d8831b-ad66-4f77-a33a-933c15c46a52
# ╠═8c0b4ab5-09da-4d8f-b001-5e15f823423c
# ╠═d3d99b9c-ad18-4a04-bb3e-f17dd542f9f3
# ╠═a2d919a5-a395-40fb-8f93-db742f8a77c2
# ╟─eb9abf2e-372c-4f79-afbe-772b90eff9ad
# ╟─9f3ae7b5-51b3-48bc-b4db-b7236ba30682
# ╟─2da8a5b1-b168-41d9-baa8-d65a4ef5c4c0
# ╠═d44407de-8c9c-42fa-b1a2-ae02b826eccc
# ╠═ae268316-c058-4db8-9b71-57b0d9425274
# ╠═b53b9d28-4c25-4fb8-a3e4-599b0e385121
# ╟─e7ce7fd4-cfa4-4cc6-84a2-7e20ed2f4e5c
# ╠═29f36902-e355-4b02-b7b0-c4db12c47d33
# ╟─673e9320-ea30-4831-ad85-ba7936293ee2
# ╠═f35f419a-94dd-4051-a533-4b1ec9a4c9ec
# ╟─9661e4fc-55e1-4c2c-a3ad-515cdac3b514
# ╠═90298676-fda7-4168-8a40-7ff53e7c761b
# ╠═7a819522-55e4-4547-a3e0-e047e74cfb6b
# ╟─cebabf33-e769-47bd-b6f1-ddf525fea895
# ╠═719f206a-5b9f-4d78-8778-1d89edb2bc4d
# ╟─1d7f442f-c057-4379-8a40-c6ce3646ad5c
# ╠═d6e1c6c7-060d-4c2f-8054-d8f33f54bd55
# ╟─da41b22e-114d-4eee-81d0-73e6f3b45242
# ╠═59c83a22-a4cc-4b51-a1cc-5eb39588eacd
# ╠═0c3ca093-e011-409c-be58-e617e02f5d4c
# ╠═de119a22-b695-4b4f-8e04-b7d68ec1e91b
# ╠═5d8aca85-f12d-4a6e-84e4-781d45b65742
# ╠═b8cd6ad1-d323-4888-bbd1-5deba5a5870d
# ╠═441a39a0-a7de-47db-8539-12dee30b8312
# ╟─6203d95b-13d2-48a6-a69f-39e33e2edcdb
# ╠═4d8e81c1-dbec-4379-9ab7-a585a369582d
# ╠═4b75cde6-db35-11eb-1b95-bba4e05edb2f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# *******************************************
#  Coordinate Tools for specific dimensions
# *******************************************
# +++++ Coordinate in X-space ++++++++++
const XKGENS = (
    (:x_Xgen, :k_Kgen, 1, 1),
    (:xy_Xgen, :kl_Kgen, 2, 1),
    (:xy_Ygen, :kl_Lgen, 2, 2),
    (:xyz_Xgen, :klm_Kgen, 3, 1),
    (:xyz_Ygen, :klm_Lgen, 3, 2),
    (:xyz_Zgen, :klm_Mgen, 3, 3)
)

for (fnx, fnk, dim, axis) = XKGENS
    @eval begin
        # X generators
        function $fnx(config::ConfigFFT{T,$dim}) where T
            XFunc(copy(config.Xcoords[$axis]), config)
        end

        # K generators
        function $fnk(config::ConfigFFT{T,$dim}) where T
            KFunc(copy(config.Kcoords[$axis]), config)
        end
    end
end

# *******************************************
#  Low/High-pass Filters
# *******************************************
const PASSFILTERS = (
    (:k_highpass_k, 1, :K_highpass_K),
    (:kl_highpass_k, 2, :K_highpass_K),
    (:klm_highpass_klm, 3, :K_highpass_K),
    (:k_lowpass_k, 1, :K_lowpass_K),
    (:kl_lowpass_kl, 2, :K_lowpass_K),
    (:klm_lowpass_klm, 3, :K_lowpass_K)
)

for (fn, dim, fnK) = PASSFILTERS
    @eval begin
        KF = KFunc{T,$dim} where T
        NT = NTuple{$dim,Int}
        $fn(f::KF, nwaves::NT) = $fnK(f, nwaves)
    end
end

# *******************************************
#  Fourier Transformation
# *******************************************
const FFTS = (
    (:k_x, :x_k, 1),
    (:kl_xy, :xy_kl, 2),
    (:klm_xyz, :xyz_klm, 3)
)
for (kx, xk, dim) = FFTS
    @eval $kx(f::XFunc{T,$dim} where T) = K_X(f)
    @eval $xk(f::KFunc{T,$dim} where T) = X_K(f)
end

# *******************************************
#  Diffrential Operators
# *******************************************
# +++++ basics +++++
const DIFFS = (
    (:k_∂x_k, :x_∂x_x, 1, 1),
    (:kl_∂x_kl, :xy_∂x_xy, 2, 1),
    (:kl_∂y_kl, :xy_∂y_xy, 2, 2),
    (:klm_∂x_klm, :xyz_∂x_xyz, 3, 1),
    (:klm_∂y_klm, :xyz_∂y_xyz, 3, 2),
    (:klm_∂z_klm, :xyz_∂z_xyz, 3, 3)
)

for (opk, opx, dim, axis) = DIFFS
    @eval begin
        KF = KFunc{T,$dim} where T
        $opk(f::KF) = K_∂Xaxis_K(f, $axis)
        XF = XFunc{T,$dim} where T
        $opx(f::XF) = X_∂Xaxis_X(f, $axis)
    end
end

# +++++ laplacian +++++
const LAPLA = (
    (:k_laplacian_k, :x_laplacian_x,
        :k_Δ_k, :x_Δ_x, 1),
    (:kl_laplacian_kl, :xy_laplacian_xy,
        :kl_Δ_kl, :xy_Δ_xy, 2),
    (:klm_laplacian_klm, :xyz_laplacian_xyz,
        :klm_Δ_klm, :xyz_Δ_xyz, 3),
)

for (klap, xlap, kΔ, xΔ, dim) = LAPLA
    @eval begin
        KF = KFunc{T,$dim} where T
        $klap(f::KF) = K_laplacian_K(f)
        XF = XFunc{T,$dim} where T
        $xlap(f::XF) = X_laplacian_X(f)
        $kΔ = $klap
        $xΔ = $xlap
    end
end

# +++++ laplacian +++++
const LAPINV = (
    (:k_laplainv_k, :x_laplainv_x,
        :k_Δ⁻¹_k, :x_Δ⁻¹_x, 1),
    (:kl_laplainv_kl, :xy_laplainv_xy,
        :kl_Δ⁻¹_kl, :xy_Δ⁻¹_xy, 2),
    (:klm_laplainv_klm, :xyz_laplainv_xyz,
        :klm_Δ⁻¹_klm, :xyz_Δ⁻¹_xyz, 3),
)

for (klinv, xlinv, kΔ⁻¹, xΔ⁻¹, dim) = LAPINV
    @eval begin
        KF = KFunc{T,$dim} where T
        $klinv(f::KF) = K_laplainv_K(f)
        XF = XFunc{T,$dim} where T
        $xlinv(f::XF) = X_laplainv_X(f)
        $kΔ⁻¹ = $klinv
        $xΔ⁻¹ = $xlinv
    end
end

# +++++ vector analysis +++++
# 2-dimensional
function kl2_grad_kl(kl_func::KFunc{T,2})::Vector{KFunc{T,2}} where T
    return [
        kl_∂x_kl(kl_func),
        kl_∂y_kl(kl_func)
    ]
end

function kl_rot_kl2(
            kl2_func::Vector{KFunc{T,2}}
        )::KFunc{T,2} where T

    if length(kl2_func) != 2
        return println("ERROR")
    end
    return (
        kl_∂x_kl(kl2_func[2])
        - kl_∂y_kl(kl2_func[1])
    )
end

function kl_div_kl2(
            kl2_func::Vector{KFunc{T,2}}
        )::KFunc{T,2} where T

    if length(kl2_func) != 2
        return println("ERROR")
    end
    return (
        kl_∂x_kl(kl2_func[1])
        + kl_∂y_kl(kl2_func[2])
    )
end

# 3-dimensional
function klm3_grad_klm(
            klm_func::KFunc{T,3}
        )::Vector{KFunc{T,3}} where T

    return [
        klm_∂x_klm(klm_func),
        klm_∂y_klm(klm_func),
        klm_∂z_klm(klm_func)
    ]

end

function klm3_rot_klm3(
            klm3_func::Vector{KFunc{T,3}}
        )::Vector{KFunc{T,3}} where T

    if length(klm3_func) != 3
        return println("ERROR")
    end
    return [
        klm_∂y_klm(klm3_func[3]) - klm_∂z_klm(klm3_func[2]),
        klm_∂z_klm(klm3_func[1]) - klm_∂x_klm(klm3_func[3]),
        klm_∂x_klm(klm3_func[2]) - klm_∂y_klm(klm3_func[1])
    ]

end

function klm_div_klm3(
            klm3_func::Vector{KFunc{T,3}}
        )::KFunc{T,3} where T

    if length(klm3_func) != 3
        return println("ERROR")
    end
    return (
        klm_∂x_klm(klm3_func[1])
        + klm_∂y_klm(klm3_func[2])
        + klm_∂z_klm(klm3_func[3])
    )

end

# *******************************************
#  Analysis Tools
# *******************************************
const TOOLS = (
    (:integ_x, :norm_x, :l2inpr_x_x, 1),
    (:integ_xy, :norm_xy, :l2inpr_xy_xy, 2),
    (:integ_xyz, :norm_xyz, :l2inpr_xyz_xyz, 3)
)

for (integ, norm, inpr, dim) = TOOLS
    @eval begin
        XF = XFunc{T,$dim} where T
        $integ(f::XF) = integ_X(f)
        $norm(f::XF, p::Real=2) = norm_X(f, p)
        $inpr(f::XF, g::XF) = l2inpr_X_X(f, g)
    end
end

max_stars_repo_name: mat-der-D/SpectralMethodFFTW

max_issues_repo_name: mat-der-D/SpectralMethodFFTW

max_forks_repo_name: mat-der-D/SpectralMethodFFTW


text:
load("debug.jl")

module TestAnalysis
using Base, Debug

code = quote
    for i=1:3
        println(i)
    end
end

#translated = translate(code)
#println(translated)

println(argpositions(:(for i=5;end)))
println(argpositions(:{x for x=1:5,y=1:3}))

end


text:
using Test
using Merlin
using Merlin.CUDA

tests = ["functions"]

for t in tests
    path = joinpath(dirname(@__FILE__), "$t.jl")
    println("$path ...")
    include(path)
end


text:
using Base.Test

all_tests = []
for (root, dirs, files) in walkdir(".")
    for filename in files
        !endswith(filename, ".jl") && continue
        "runtests.jl" == filename && continue
        filepath = joinpath(root, filename)[3:end]
        !isempty(ARGS) && !any(x->startswith(filepath, x), ARGS) && continue
        push!(all_tests, filepath)
    end
end

for (idx, filepath) in enumerate(all_tests)
    numbering = string(idx, /, length(all_tests))
    ts = @testset "$numbering $filepath" begin
        include(filepath)
    end
end


text:
# Implementation of main console executable and related I/O.

using YAML
using JLD

import .ArgumentParse
import .Model: Image
import .ParallelRun: infer_box
using .SDSSIO

# read universal and dataset-specific configurations
function read_config(filename)
    data = open(YAML.load, filename)

    # universal settings
    cfg_keys = ["min_radius_pix",
                "num_ais_temperatures",
                "num_ais_samples",
                "num_joint_vi_iters"]
    # dictionary of cfg_keys actually present in configuration file.
    cfg_dict = Dict(k=>data[k] for k in cfg_keys if haskey(data, k))
    cfg = Config(;cfg_dict...)

    # dataset-specific settings
    datasets = Dict()
    for (name, dataset_config) in data["datasets"]
        if name == "sdss"
            basedir = dataset_config["basedir"]
            iostrategy = Symbol(get(dataset_config, "iostrategy", "plain"))
            dirlayout = Symbol(get(dataset_config, "dirlayout", "celeste"))
            slurp = get(dataset_config, "slurp", false)
            compressed = get(dataset_config, "compressed", false)
            datasets[name] = SDSSDataSet(basedir;
                                         dirlayout = dirlayout,
                                         compressed = compressed,
                                         iostrategy = iostrategy,
                                         slurp = slurp)
        else
            error("unrecognized dataset name: $(name)")
        end
    end

    return cfg, datasets
end

"""
Save provided results to a JLD file.
"""
function save_results(outdir, ramin, ramax, decmin, decmax, results)
    fname = @sprintf("%s/celeste-%.4f-%.4f-%.4f-%.4f.jld",
                     outdir, ramin, ramax, decmin, decmax)
    JLD.save(fname, "results", results)
end
save_results(outdir, box::BoundingBox, results) =
    save_results(outdir, box.ramin, box.ramax, box.decmin, box.decmax, results)

# command-line script
function main(console_args::Vector{String}=ARGS)

    parser = ArgumentParse.ArgumentParser(program_name="infer-box")
    ArgumentParse.add_argument(parser, "datasets",
                               help="comma-separated dataset names. Each dataset must appear in config file.")
    ArgumentParse.add_argument(parser, "ramin", arg_type=Float64)
    ArgumentParse.add_argument(parser, "ramax", arg_type=Float64)
    ArgumentParse.add_argument(parser, "decmin", arg_type=Float64)
    ArgumentParse.add_argument(parser, "decmax", arg_type=Float64)
    ArgumentParse.add_argument(parser, "--config",
                               help="YAML configuration file.",
                               default="celeste.yml")
    ArgumentParse.add_argument(parser, "--method",
                               help="inference method: {joint_vi, single_vi, mcmc}",
                               default="joint_vi")

    args = ArgumentParse.parse_args(parser, console_args)

    # read config
    config, known_datasets = read_config(args["config"])

    # initialize data sets
    names = split(args["datasets"], ",")
    datasets = [known_datasets[name] for name in names]

    box = BoundingBox(args["ramin"], args["ramax"],
                      args["decmin"], args["decmax"])

    # load images from all datasets
    images = Image[]
    for dataset in datasets
        append!(images, load_images(dataset, box))
    end

    results = infer_box(images, box; config = config,
                        method = Symbol(args["method"]))
    save_results(".", box, results)
    return 0
end


text:
module ResourcePools

export AbstractPooledResource,
    PooledResource,
    ref_count,
    resource,
    release!,
    retain!

import Base: length, take!
export AbstractResourcePool,
    ResourcePool,
    length,
    take!

export AbstractPooledArray,
    AbstractPooledDenseArray,
    PooledAbstractArray,
    PooledArray

include("pooled_resource.jl")
include("resource_pool.jl")
include("pooled_abstract_array.jl")
include("pooled_dense_array.jl")

end # module


text:
module TestBasic
using ZulipSnippetBot
using ZulipSnippetBot: ZulipRequest, process, OPTS, Message
using JSON3
using Test
using SQLite
using Logging

global_logger(ConsoleLogger(stderr, Logging.Warn))

f = readlines(joinpath(@__DIR__, "data", "requests.json"))
const PRIVATE_MSG = f[1]
const MENTION_MSG = f[2]
const WRONG_TOKEN = f[3]
const NO_TOKEN = f[4]
const NO_SENDER = f[5]
const NO_DATA = f[6]
const PRIVATE_HELP = f[7]
const MENTION_HELP = f[8]

content(x) = JSON3.read(x).content
extract_codeid(x) = match(r"Snippet codeid: `([^`]+)`", content(x))[1]

const TOKEN = "3xP6YwzEQr9dr2TCrzn28Yrvx4FtFgC0"
setupbot!(; token = TOKEN)

@testset "test input" begin
    msg = JSON3.read(PRIVATE_MSG, ZulipRequest)
    @test msg.data == "Hello"
    @test msg.token == "3xP6YwzEQr9dr2TCrzn28Yrvx4FtFgC0"
    @test msg.message.sender_id == 1234

    msg = JSON3.read(MENTION_MSG, ZulipRequest)
    @test msg.data == "@**SnippetBot** Bang"
    @test msg.token == "3xP6YwzEQr9dr2TCrzn28Yrvx4FtFgC0"
    @test msg.message.sender_id == 1234
end

@testset "Simple errors" begin
    msg = JSON3.read(WRONG_TOKEN, ZulipRequest)
    @test content(process(msg, nothing, OPTS[])) == "Incorrect token, verify SnippetBot server configuration"

    msg = JSON3.read(NO_TOKEN, ZulipRequest)
    @test content(process(msg, nothing, OPTS[])) == "Wrong message"

    msg = JSON3.read(NO_SENDER, ZulipRequest)
    @test content(process(msg, nothing, OPTS[])) == "Wrong message"

    msg = JSON3.read(NO_DATA, ZulipRequest)
    @test content(process(msg, nothing, OPTS[])) == "Wrong message"
end

@testset "help" begin
    msg = JSON3.read(PRIVATE_HELP, ZulipRequest)
    @test startswith(content(process(msg, nothing, OPTS[])), "Currently following")

    msg = JSON3.read(MENTION_HELP, ZulipRequest)
    @test startswith(content(process(msg, nothing, OPTS[])), "Currently following")
end

@testset "save/show" begin
    dbpath, _ = mktemp()
    db = SQLite.DB(dbpath)
    ZulipSnippetBot.up(db)

    msg = ZulipRequest("save #qwe abc", TOKEN, Message(1))
    res = content(process(msg, db, OPTS[]))
    @test startswith(res, "Snippet codeid: ")

    cid = match(r"Snippet codeid: `([^`]+)`", res)[1]
    msg = ZulipRequest("show $(cid)", TOKEN, Message(1))
    @test content(process(msg, db, OPTS[])) == "abc"

    msg = ZulipRequest("show xxx", TOKEN, Message(1))
    @test content(process(msg, db, OPTS[])) == "Codeid xxx is not found"
end

@testset "save/list" begin
    dbpath, _ = mktemp()
    db = SQLite.DB(dbpath)
    ZulipSnippetBot.up(db)

    msg = ZulipRequest("save #qwe abc", TOKEN, Message(1))
    cid1 = extract_codeid(process(msg, db, OPTS[]))

    msg = ZulipRequest("save #qwe #zxc abc2", TOKEN, Message(1))
    cid2 = extract_codeid(process(msg, db, OPTS[]))

    msg = ZulipRequest("save #foo abc3", TOKEN, Message(1))
    cid3 = extract_codeid(process(msg, db, OPTS[]))

    msg = ZulipRequest("save #qwe foo", TOKEN, Message(2))
    cid4 = extract_codeid(process(msg, db, OPTS[]))

    msg = ZulipRequest("list #qwe", TOKEN, Message(1))
    # These are dangerous, because order is not determined... Maybe I can use mocks here...
    @test content(process(msg, db, OPTS[])) == "**codeid**: $cid1, **tags**: #qwe\nabc\n\n**codeid**: $cid2, **tags**: #qwe #zxc\nabc2"

    msg = ZulipRequest("list #zxc", TOKEN, Message(1))
    @test content(process(msg, db, OPTS[])) == "**codeid**: $cid2, **tags**: #qwe #zxc\nabc2"

    msg = ZulipRequest("list", TOKEN, Message(1))
    @test content(process(msg, db, OPTS[])) == "**codeid**: $cid1, **tags**: #qwe\nabc\n\n**codeid**: $cid2, **tags**: #qwe #zxc\nabc2\n\n**codeid**: $cid3, **tags**: #foo\nabc3"

    msg = ZulipRequest("list #qwe", TOKEN, Message(2))
    @test content(process(msg, db, OPTS[])) == "**codeid**: $cid4, **tags**: #qwe\nfoo"

    msg = ZulipRequest("list #foo", TOKEN, Message(2))
    @test content(process(msg, db, OPTS[])) == "No snippets found"

    msg = ZulipRequest("list", TOKEN, Message(3))
    @test content(process(msg, db, OPTS[])) == "No snippets found"
end

end # module

max_stars_repo_name: Arkoniak/ZulipSnippetBot

max_issues_repo_name: Arkoniak/ZulipSnippetBot

max_forks_repo_name: Arkoniak/ZulipSnippetBot


text:
let sys = DenseSystem(test_params...),
    beta = 12.34,
    dbeta = 1e-4,
    P = 256,
    num_samples = 100,
    sampling_sys = DiagonalSystem(test_params_sampling...)

    # Only testing that it runs to completion.
    @test_nowarn SamplingFiniteDifference(sys, beta, dbeta, P, num_samples; progress_output=devnull)
    @test_nowarn SamplingFiniteDifference(sys, beta, dbeta, P, num_samples; sampling_sys=sampling_sys, progress_output=devnull)
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    NeighborhoodSearcher(object, neighborhood)

A method for searching neighbors in spatial `object` inside `neighborhood`.
"""
struct NeighborhoodSearcher{O,N,K} <: AbstractNeighborSearcher
  # input fields
  object::O
  neigh::N

  # state fields
  kdtree::K
end

function NeighborhoodSearcher(object::O, neigh::N) where {O,N}
  kdtree = neigh isa BallNeighborhood ? KDTree(coordinates(object), metric(neigh)) : nothing
  NeighborhoodSearcher{O,N,typeof(kdtree)}(object, neigh, kdtree)
end

# search method for any neighborhood
function search(xₒ::AbstractVector, searcher::NeighborhoodSearcher{O,N,T};
                mask=nothing) where {O,N,T}
  object = searcher.object
  neigh  = searcher.neigh
  locs   = mask ≠ nothing ? view(1:npoints(object), mask) : 1:npoints(object)

  x = MVector{ndims(object),coordtype(object)}(undef)

  neighbors = Vector{Int}()
  @inbounds for loc in locs
    coordinates!(x, object, loc)
    if isneighbor(neigh, xₒ, x)
      push!(neighbors, loc)
    end
  end

  neighbors
end

# search method for ball neighborhood
function search(xₒ::AbstractVector, searcher::NeighborhoodSearcher{O,N,T};
                mask=nothing) where {O,N<:BallNeighborhood,T}
  locs = inrange(searcher.kdtree, xₒ, radius(searcher.neigh))
  if mask ≠ nothing
    neighbors = Vector{Int}()
    @inbounds for loc in locs
      if mask[loc]
        push!(neighbors, loc)
      end
    end
    neighbors
  else
    locs
  end
end


text:
using Random, Plots
using DecisionProgramming

scale(x::Float64, low::Float64, high::Float64) = x * (high - low) + low

function random_probability_distribution(rng::AbstractRNG, n_neg::Int, n_pos::Int, low::Real, high::Real)
    low < high || error("")
    u_neg = rand(rng, n_neg)
    u_neg = scale.(u_neg, low, 0.0)
    u_pos = rand(rng, n_pos)
    u_pos = scale.(u_pos, 0.0, high)
    u = [u_neg; u_pos]
    p = rand(rng, n_neg + n_pos)
    p = p / sum(p)
    i = sortperm(u)
    return u[i], p[i]
end

hair(plt, x, y) = plot!(plt, [x, x], [0, y], linewidth=3, color=:grey, alpha=0.5, label=false)

function plot_distribution(x, y, x_α; plt=plot())
    for (x2, y2) in zip(x, y)
        hair(plt, x2, y2)
    end
    tail = x .≤ x_α
    head = x .> x_α
    plot!(plt, ylims = (0, 1.15 * maximum(y)))
    plot!(plt, x[tail], y[tail], linewidth=0, markershape=:circle, markercolor=:darkred,
          label="Tail")
    plot!(plt, x[head], y[head], linewidth=0, markershape=:circle, markercolor=:lightblue,
          label="Head")
    return plt
end

function plot_risk_measures(u, p, α)
    mean = sum(u.*p)
    VaR = value_at_risk(u, p, α)
    CVaR = conditional_value_at_risk(u, p, α)

    plt1 = plot(
        title="Probability distribution",
        legend=false,
        xlabel="x",
        ylabel="f(x)"
    )
    plot_distribution(u, p, VaR; plt=plt1)
    plot!(plt1, [mean], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:green,
        label="E")
    plot!(plt1, [VaR], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:darkblue,
        label="VaR")
    plot!(plt1, [CVaR], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:darkorange,
        label="CVaR")

    plt2 = plot(
        title="Cumulative distribution",
        xlabel="x",
        ylabel="F(x)"
    )
    plot_distribution(u, cumsum(p), VaR; plt=plt2)
    plot!(plt2, u, [α for _ in u],
          legend=false, linewidth=2, label="α", linecolor=:darkorange)
    plot!(plt2, [VaR], [0],
        linewidth=0, markershape=:diamond,
        markersize=6, markercolor=:darkblue,
        label="VaR")

    plt = plot(plt1, plt2, layout=(2, 1), legend=:outerright, size=(720, 400))
end

# rng = MersenneTwister(1)
# u, p = random_probability_distribution(rng, 4, 8, -2.0, 3.0)
# α = 0.25
# plt = plot_risk_measures(u, p, α)
# savefig(plt, "../../docs/src/decision-programming/figures/risk_measures.svg")
# savefig(plt, "risk_measures.svg")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "P4est"
version = v"2.3.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://p4est.github.io/release/p4est-2.3.1.tar.gz", "be66893b039fb3f27aca3d5d00acff42c67bfad5aa09cea9253cdd628b2bdc9a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd p4est-*
if [[ "${target}" == *-freebsd* ]]; then
  export LIBS="-lm"
elif [[ "${target}" == x86_64-linux-musl ]]; then
  # We can't run Fortran programs for the native platform, so a check that the
  # Fortran compiler works would fail.  Small hack: swear that we're
  # cross-compiling.  See:
  # https://github.com/JuliaPackaging/BinaryBuilderBase.jl/issues/50.
  sed -i 's/cross_compiling=no/cross_compiling=yes/' configure
  sed -i 's/cross_compiling=no/cross_compiling=yes/' sc/configure
fi

# Set default preprocessor and linker flags
export CPPFLAGS="-I${includedir}"
export LDFLAGS="-L${libdir}"

# Special Windows treatment
FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
  # Set linker flags only at build time (see https://docs.binarybuilder.org/v0.3/troubleshooting/#Windows)
  FLAGS+=(LDFLAGS="$LDFLAGS -no-undefined")

  # Add manual definitions to fix missing `htonl` according to `INSTALL_WINDOWS` file
  # (see https://github.com/cburstedde/p4est/blob/master/INSTALL_WINDOWS)
  sed -i "1s/^/#define htonl(_val) ( ((uint16_t)(_val) \& 0xff00) >> 8 | ((uint16_t)(_val) \& 0xff) << 8 )\n/" src/p4est_algorithms.c src/p8est_algorithms.c src/p6est.c src/p4est_ghost.c
fi

# Configure, build, install
# Note: BLAS is disabled since it is only needed for SC if it is used outside of p4est
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-static --without-blas
make -j${nproc} "${FLAGS[@]}"
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libp4est", :libp4est),
    LibraryProduct("libsc", :libsc)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8.1.0", julia_compat="1.6")

max_stars_repo_name: luoq/Yggdrasil

max_issues_repo_name: luoq/Yggdrasil

max_forks_repo_name: luoq/Yggdrasil


text:
insts = OrderedDict()
insts["minimal"] = [
    ((5, 2, 2, true, true, true),),
    ((5, 1, 2, false, true, true),),
    ((5, 1, 2, false, true, false),),
    ((5, 1, 2, true, false, false), :ExpPSD),
    ((5, 1, 2, true, false, false), :SOCExpPSD),
    ((:iris, 2, true, false, true),),
    ]
insts["fast"] = [
    ((10, 1, 10, true, false, false), :ExpPSD),
    ((10, 1, 10, true, false, false), :SOCExpPSD),
    ((100, 1, 250, true, true, true),),
    ((100, 2, 5, false, true, true),),
    ((200, 2, 20, true, true, true),),
    ((50, 3, 2, true, true, true),),
    ((50, 3, 2, false, false, true),),
    ((50, 3, 2, true, false, true), :ExpPSD),
    ((50, 3, 2, true, false, true), :SOCExpPSD),
    ((50, 3, 4, true, true, true),),
    ((500, 3, 14, true, true, true),),
    ((100, 8, 2, true, true, true),),
    ((250, 4, 4, true, true, true),),
    ((250, 4, 4, false, false, true),),
    ((200, 32, 2, true, false, true),),
    ((:iris, 4, false, false, false),),
    ((:iris, 5, true, false, false),),
    ((:iris, 6, true, true, true),),
    ((:cancer, 4, true, true, true),),
    ]
insts["various"] = [
    ((50, 2, 16, true, false, false),),
    ((50, 2, 16, true, false, false), :ExpPSD),
    ((50, 2, 16, true, false, false), :SOCExpPSD),
    ((50, 2, 16, false, false, false),),
    ((50, 2, 16, false, false, true),),
    ((50, 8, 4, true, false, false),),
    ((50, 8, 4, true, false, false), :ExpPSD),
    ((50, 8, 4, true, false, false), :SOCExpPSD),
    ((50, 8, 4, false, false, false),),
    ((50, 8, 4, false, false, true),),
    ((50, 32, 2, true, false, false),),
    ((50, 32, 2, true, false, false), :ExpPSD),
    ((50, 32, 2, true, false, false), :SOCExpPSD),
    ((:iris, 6, true, true, true),),
    ((:cancer, 4, true, true, true),),
    ]
return (DensityEstJuMP, insts)


text:
module Functional

# Functions

export range, repeat, first, last, lastindex, empty, reverse
export zip, unzip
export enumerate
export map, filter
export all, any
export take, takewhile
export skip, skipwhile
export reduce, foldright, foldleft
export sum, mul
export avg, min, max
export find, join, concat

include("functions.jl")

# Matching

export match, default, when, or

include("matching.jl")

end # module


text:
module ImageFeatures

# package code goes here
using Images, Distributions
using SparseArrays
import Random.seed!

include("core.jl")
include("const.jl")
include("corner.jl")
include("lbp.jl")
include("glcm.jl")
include("brief.jl")
include("orb.jl")
include("freak.jl")
include("brisk.jl")
include("houghtransform.jl")
include("hog.jl")

export Keypoint, Keypoints, Feature, Features, Params, BRIEF, ORB, FREAK, BRISK, HOG

export
    #Core
    create_descriptor,
    hamming_distance,
    match_keypoints,
    grade_matches,

    #Local Binary Patterns
    lbp,
    modified_lbp,
    direction_coded_lbp,
    lbp_original,
    lbp_uniform,
    lbp_rotation_invariant,
    multi_block_lbp,

    #Gray Level Co Occurence Matrix
    glcm,
    glcm_symmetric,
    glcm_norm,
    glcm_prop,
    max_prob,
    contrast,
    ASM,
    IDM,
    glcm_entropy,
    energy,
    dissimilarity,
    correlation,
    glcm_mean_ref,
    glcm_mean_neighbour,
    glcm_var_ref,
    glcm_var_neighbour,

    #Corners
    corner_orientations,

    #BRIEF
    random_uniform,
    random_coarse,
    gaussian,
    gaussian_local,
    center_sample,

    #Lines
    hough_transform_standard,

    #Circles
    hough_circle_gradient

"""
    desc, keypoints = create_descriptor(img, keypoints, params)
    desc, keypoints = create_descriptor(img, params)

Create a descriptor for each entry in `keypoints` from the image `img`. `params` specifies the parameters for any of several descriptors:
 - [`BRIEF`](@ref)
 - [`ORB`](@ref)
 - [`BRISK`](@ref)
 - [`FREAK`](@ref)
 - [`HOG`](@ref)

Some descriptors support discovery of the `keypoints` from [`fastcorners`](@ref).
"""
create_descriptor

end


text:
using YaoGraphRegister
using Test

@testset "YaoGraphRegister.jl" begin
    # Write your own tests here.
end

using LightGraphs, GraphPlot, Compose

g = SimpleGraph(10, 10)
gplot(g)

max_stars_repo_name: QuantumBFS/YaoGraphs

max_issues_repo_name: QuantumBFS/YaoGraphs

max_forks_repo_name: QuantumBFS/YaoGraphs


text:
using TensorOperations

"""

     interpolate( interp :: InterpolationType,
                  f  :: Array{Float64,1},
                  xe :: Array{Float64,1})

Computes interpolation of a 1D function via Hermite-tensor method.

Arguments:
 - interp : inteprolation type (Hermite or Radial)
 - f      : 1d array containing values at interpolation nodes (Chebyshev or Uniform)
 - xe     : vector of evaluation points

  © Anna Yurova, 2017

"""
function interpolate(interp :: InterpolationType,
                     f      :: Array{Float64,1},
		     xe     :: Array{Float64,1} )

   x = interp(xe)

   @tensoropt begin
      s[e] := x[e, c]*f[c]
   end

   s

end

"""
     interpolate( interp_x :: InterpolationType,
                  interp_y :: InterpolationType,
                  f  :: Array{Float64,2},
                  xe :: Array{Float64,1},
                  ye :: Array{Float64,1})

Computes interpolation of a 2D function via Hermite-tensor method.

Arguments:
 - interp : inteprolation type (Hermite or Radial)
 - f      : 2d array containing values at interpolation nodes (Chebyshev or Uniform)
 - xe,ye  : vector of evaluation points
"""
function interpolate(interp_x :: InterpolationType,
                     interp_y :: InterpolationType,
                     f        :: Array{Float64,2},
		     xe       :: Array{Float64,1},
		     ye       :: Array{Float64,1})

   x = interp_x(xe)
   y = interp_y(ye)

   @tensoropt begin
       s[e1,e2] := x[e1,c1]*y[e2,c2]*f[c1,c2]
   end

   s

end

"""
     interpolate( interp_x :: InterpolationType,
                  interp_y :: InterpolationType,
                  interp_z :: InterpolationType,
                  f  :: Array{Float64,3},
                  xe :: Array{Float64,1},
                  ye :: Array{Float64,1},
                  ze :: Array{Float64,1})

Computes interpolation of a 3D function via Hermite-tensor method.

Arguments:
 - interp_x : inteprolation type (Hermite or Radial)
 - interp_y : inteprolation type (Hermite or Radial)
 - interp_z : inteprolation type (Hermite or Radial)
 - f        : 3d array containing values at interpolation nodes (Chebyshev or Uniform)

Computing s = ( z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)

"""
function interpolate(interp_x :: InterpolationType,
                     interp_y :: InterpolationType,
                     interp_z :: InterpolationType,
                     f        :: Array{Float64,3},
		     xe       :: Array{Float64,1},
		     ye       :: Array{Float64,1},
		     ze       :: Array{Float64,1})

    x = interp_x(xe)
    y = interp_y(ye)
    z = interp_z(ze)

    @tensoropt begin
        s[e1,e2,e3] := x[e1,c1]*y[e2,c2]*z[e3,c3]*f[c1,c2,c3]
    end
    s
end

"""
     interpolate( interp_x :: InterpolationType,
                  interp_y :: InterpolationType,
                  interp_z :: InterpolationType,
                  interp_v :: InterpolationType,
                  f  :: Array{Float64,4},
                  xe :: Array{Float64,1},
                  ye :: Array{Float64,1},
                  ze :: Array{Float64,1},
                  ve :: Array{Float64,1})

Computes interpolation of a 4D function via Hermite-tensor method.

Arguments:
 - interp_x : inteprolation type (Hermite or Radial)
 - interp_y : inteprolation type (Hermite or Radial)
 - interp_z : inteprolation type (Hermite or Radial)
 - interp_v : inteprolation type (Hermite or Radial)
 - f        : 4d array containing values at interpolation nodes (Chebyshev or Uniform)

Computing s = (v ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)

"""
function interpolate(interp_x :: InterpolationType,
                     interp_y :: InterpolationType,
                     interp_z :: InterpolationType,
                     interp_v :: InterpolationType,
                     f        :: Array{Float64,4},
		     xe       :: Array{Float64,1},
		     ye       :: Array{Float64,1},
		     ze       :: Array{Float64,1},
		     ve       :: Array{Float64,1})

    x = interp_x(xe)
    y = interp_y(ye)
    z = interp_z(ze)
    v = interp_v(ve)

    @tensoropt begin
        s[e1,e2,e3,e4] := x[e1,c1]*y[e2,c2]*z[e3,c3]*v[e4,c4]*f[c1,c2,c3,c4]
    end
    s
end

"""
     interpolate( interp_x :: InterpolationType,
                  interp_y :: InterpolationType,
                  interp_z :: InterpolationType,
                  interp_v :: InterpolationType,
                  interp_w :: InterpolationType,
                  f  :: Array{Float64,4},
                  xe :: Array{Float64,1},
                  ye :: Array{Float64,1},
                  ze :: Array{Float64,1},
                  ve :: Array{Float64,1},
                  we :: Array{Float64,1})

Computes interpolation of a 4D function via Hermite-tensor method.

Arguments:
 - interp_x : interpolation type (Hermite or Radial)
 - interp_y : interpolation type (Hermite or Radial)
 - interp_z : interpolation type (Hermite or Radial)
 - interp_v : interpolation type (Hermite or Radial)
 - interp_w : interpolation type (Hermite or Radial)
 - f        : 5d arrray containing values at interpolation nodes (Chebyshev or Uniform)

Extract the number of collocation and evaluation points in each dimension
Computing s = (w ⊗ v ⊗ z ⊗ y ⊗ x) vec(f) (see Sec. 4.1 of the paper)

"""
function interpolate(interp_x :: InterpolationType,
                     interp_y :: InterpolationType,
                     interp_z :: InterpolationType,
                     interp_v :: InterpolationType,
                     interp_w :: InterpolationType,
                     f        :: Array{Float64,5},
		     xe       :: Array{Float64,1},
		     ye       :: Array{Float64,1},
		     ze       :: Array{Float64,1},
		     ve       :: Array{Float64,1},
		     we       :: Array{Float64,1})

    x = interp_x(xe)
    y = interp_y(ye)
    z = interp_z(ze)
    v = interp_v(ve)
    w = interp_w(we)

    @tensoropt begin
        s[e1,e2,e3,e4,e5] := x[e1,c1]*y[e2,c2]*z[e3,c3]*v[e4,c4]*w[e5,c5]*f[c1,c2,c3,c4,c5]
    end
    s

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AvailableProvidersList <: SwaggerModel
    countries::Any # spec type: Union{ Nothing, Vector{AvailableProvidersListCountry} } # spec name: countries

    function AvailableProvidersList(;countries=nothing)
        o = new()
        validate_property(AvailableProvidersList, Symbol("countries"), countries)
        setfield!(o, Symbol("countries"), countries)
        o
    end
end # type AvailableProvidersList

const _property_map_AvailableProvidersList = Dict{Symbol,Symbol}(Symbol("countries")=>Symbol("countries"))
const _property_types_AvailableProvidersList = Dict{Symbol,String}(Symbol("countries")=>"Vector{AvailableProvidersListCountry}")
Base.propertynames(::Type{ AvailableProvidersList }) = collect(keys(_property_map_AvailableProvidersList))
Swagger.property_type(::Type{ AvailableProvidersList }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AvailableProvidersList[name]))}
Swagger.field_name(::Type{ AvailableProvidersList }, property_name::Symbol) =  _property_map_AvailableProvidersList[property_name]

function check_required(o::AvailableProvidersList)
    (getproperty(o, Symbol("countries")) === nothing) && (return false)
    true
end

function validate_property(::Type{ AvailableProvidersList }, name::Symbol, val)
end


text:
module LibGEOS

using GEOS_jll
using GeoInterface
using CEnum

export Point,
    MultiPoint,
    LineString,
    MultiLineString,
    LinearRing,
    Polygon,
    MultiPolygon,
    GeometryCollection,
    parseWKT,
    geomFromWKT,
    geomToWKT,
    readgeom,
    writegeom,
    project,
    projectNormalized,
    interpolate,
    interpolateNormalized,
    buffer,
    bufferWithStyle,
    envelope,
    intersection,
    convexhull,
    difference,
    symmetricDifference,
    boundary,
    union,
    unaryUnion,
    pointOnSurface,
    centroid,
    node,
    polygonize,
    lineMerge,
    simplify,
    topologyPreserveSimplify,
    uniquePoints,
    sharedPaths,
    snap,
    delaunayTriangulation,
    delaunayTriangulationEdges,
    constrainedDelaunayTriangulation,
    disjoint,
    touches,
    intersects,
    crosses,
    within,
    overlaps,
    equals,
    equalsexact,
    covers,
    coveredby,
    prepareGeom,
    prepcontains,
    prepcontainsproperly,
    prepcoveredby,
    prepcovers,
    prepcrosses,
    prepdisjoint,
    prepintersects,
    prepoverlaps,
    preptouches,
    prepwithin,
    isEmpty,
    isSimple,
    isRing,
    hasZ,
    isClosed,
    isValid,
    interiorRings,
    exteriorRing,
    numGeometries,
    numPoints,
    startPoint,
    endPoint,
    area,
    geomLength,
    distance,
    hausdorffdistance,
    nearestPoints,
    getPrecision,
    setPrecision,
    getXMin,
    getYMin,
    getXMax,
    getYMax,
    minimumRotatedRectangle,
    getGeometry,
    getGeometries,
    STRtree,
    query

include("libgeos_api.jl")

mutable struct GEOSError <: Exception
    msg::String
end
Base.showerror(io::IO, err::GEOSError) = print(io, "GEOSError\n\t$(err.msg)")

function geosjl_errorhandler(message::Ptr{UInt8}, userdata)
    if unsafe_string(message) == "%s"
        throw(GEOSError(unsafe_string(Cstring(userdata))))
    else
        throw(GEOSError(unsafe_string(message)))
    end
end

mutable struct GEOSContext
    ptr::GEOSContextHandle_t

    function GEOSContext()
        context = new(GEOS_init_r())
        GEOSContext_setNoticeHandler_r(context.ptr, C_NULL)
        GEOSContext_setErrorHandler_r(
            context.ptr,
            @cfunction(geosjl_errorhandler, Ptr{Cvoid}, (Ptr{UInt8}, Ptr{Cvoid}))
        )
        finalizer(context -> (GEOS_finish_r(context.ptr); context.ptr = C_NULL), context)
        context
    end
end

mutable struct WKTReader
    ptr::Ptr{GEOSWKTReader}

    function WKTReader(context::GEOSContext)
        reader = new(GEOSWKTReader_create_r(context.ptr))
        finalizer(function (reader)
            GEOSWKTReader_destroy_r(context.ptr, reader.ptr)
            reader.ptr = C_NULL
        end, reader)
        reader
    end
end

mutable struct WKTWriter
    ptr::Ptr{GEOSWKTWriter}

    function WKTWriter(
        context::GEOSContext;
        trim::Bool = true,
        outputdim::Int = 3,
        roundingprecision::Int = -1,
    )
        writer = new(GEOSWKTWriter_create_r(context.ptr))
        GEOSWKTWriter_setTrim_r(context.ptr, writer.ptr, UInt8(trim))
        GEOSWKTWriter_setOutputDimension_r(context.ptr, writer.ptr, outputdim)
        GEOSWKTWriter_setRoundingPrecision_r(context.ptr, writer.ptr, roundingprecision)
        finalizer(function (writer)
            GEOSWKTWriter_destroy_r(context.ptr, writer.ptr)
            writer.ptr = C_NULL
        end, writer)
        writer
    end
end

mutable struct WKBReader
    ptr::Ptr{GEOSWKBReader}

    function WKBReader(context::GEOSContext)
        reader = new(GEOSWKBReader_create_r(context.ptr))
        finalizer(function (reader)
            GEOSWKBReader_destroy_r(context.ptr, reader.ptr)
            reader.ptr = C_NULL
        end, reader)
        reader
    end
end

mutable struct WKBWriter
    ptr::Ptr{GEOSWKBWriter}

    function WKBWriter(context::GEOSContext)
        writer = new(GEOSWKBWriter_create_r(context.ptr))
        finalizer(function (writer)
            GEOSWKBWriter_destroy_r(context.ptr, writer.ptr)
            writer.ptr = C_NULL
        end, writer)
        writer
    end
end

function __init__()
    global _context = GEOSContext()
end

include("geos_functions.jl")
include("geos_types.jl")
include("geos_operations.jl")
include("geo_interface.jl")
include("strtree.jl")
include("deprecated.jl")
end


text:
using jInv.Mesh
using jInv.LinearSolvers
using Base.Test
using KrylovMethods
using Multigrid

println("===  Example 2D DivSigGrad ====");

domain = [0.0, 1.0, 0.0, 1.0];
n      = [50,50];
Mr     = getRegularMesh(domain,n)
G      = getNodalGradientMatrix(Mr);
m      = spdiagm(exp.(randn(size(G,1))));
Ar     = G'*m*G;
Ar     = Ar + 1e-1*norm(Ar,1)*speye(size(Ar,2));
N      = size(Ar,2);
B      = Ar*rand(N,4);

levels      = 5;
numCores 	= 8;
maxIter     = 15;
relativeTol = 1e-2;
relaxType   = "SPAI";
relaxParam  = 1.0;
relaxPre 	= 2;
relaxPost   = 2;
cycleType   ='V';
coarseSolveType = "Julia";

MG = getMGparam(levels,numCores,maxIter,relativeTol,relaxType,relaxParam,relaxPre,relaxPost,cycleType,coarseSolveType);

sSAPCG   = getSA_AMGsolver(MG, "PCG",sym=1,out=1);
X,  = solveLinearSystem(Ar,B,sSAPCG);
@test vecnorm(Ar*X-B)/vecnorm(B) < sSAPCG.tol

clear!(sSAPCG);
sSABiCG   = copySolver(sSAPCG);
sSABiCG.Krylov = "BiCGSTAB"
X,  = solveLinearSystem(Ar,B,sSABiCG);
@test vecnorm(Ar*X-B)/vecnorm(B) < sSABiCG.tol


text:
println("Parent running.")
@async(begin sleep(1); println("This is the child process."); sleep(2); println("Child again.") end)
sleep(2)
println("This is the parent process again.")
sleep(2)
println("Parent again.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Functions for estimating a template CCF

Author: Eric Ford
Created: October 2020
"""

""" calc_normalized_ccfs( ccfs )
Normalizes each spectrum by fitting a line to the region outside the center
"""
function fit_ccf_normalizations( v_grid::A1, ccfs::A2; v_mid::Real, dv_min::Real = RvSpectMLBase.max_bc, dv_max::Real = Inf ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2} }
    @assert zero(dv_min) < dv_min < Inf
    @assert dv_min < dv_max
    idx = findall( x->dv_min < abs(x-v_mid) < dv_max , v_grid )
    @assert length(idx) >= 2
    map(col-> mean(ccfs[idx,col] ), 1:size(ccfs,2) )
end

function calc_normalized_ccfs( v_grid::A1, ccfs::A2; v_mid::Real, dv_min::Real = RvSpectMLBase.max_bc, dv_max::Real = Inf ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2} }
    norms = fit_ccf_normalizations(v_grid,ccfs,v_mid=v_mid,dv_min=dv_min, dv_max=dv_max)
    ccfs_out = copy(ccfs)
    for i in 1:size(ccfs,2)
        ccfs_out[:,i] ./= norms[i]
    end
    return ccfs_out
end

function calc_normalized_ccfs( v_grid::A1, ccfs::A2, ccf_vars::A3; v_mid::Real, dv_min::Real = RvSpectMLBase.max_bc, dv_max::Real = Inf ) where {T1<:Real, A1<:AbstractArray{T1,1}, T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,2}  }
    norms = fit_ccf_normalizations(v_grid,ccfs,v_mid=v_mid,dv_min=dv_min, dv_max=dv_max)
    ccfs_out = copy(ccfs)
    ccf_vars_out = copy(ccf_vars)
    for i in 1:size(ccfs,2)
        ccfs_out[:,i] ./= norms[i]
        ccf_vars_out[:,i] ./= norms[i]^2
    end
    return (ccfs=ccfs_out, ccf_vars=ccf_vars_out)
end

""" calc_normalized_ccfs( ccfs )
Normalizes each spectrum by its maximum value.
"""
function calc_normalized_ccfs( ccfs::A2 ) where {T2<:Real, A2<:AbstractArray{T2,2} }
    return ccfs ./maximum(ccfs,dims=1)
end

function calc_normalized_ccfs( ccfs::A2, ccf_vars::A3 ) where { T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,2}  }
    max_ccfs = maximum(ccfs,dims=1)
    ccfs_norm = ccfs ./ max_ccfs
    ccf_vars_norm = ccf_vars ./ max_ccfs
    return (ccfs=ccfs_norm, ccf_vars=ccf_vars_norm)
end

""" calc_ccf_template( ccfs, [ccf_vars] ; assume_normalized = false )
Calculates ccf template
Warning: uses maximum CCF for normalization, unless you normalize manually.
"""
function calc_ccf_template( ccfs::A2; assume_normalized::Bool = false ) where {T2<:Real, A2<:AbstractArray{T2,2} }
    ccfs_norm = assume_normalized ? ccfs : calc_normalized_ccfs(ccfs)
    ccf_template = mean(ccfs_norm,dims=2)
    est_ccf_vars_norm = var(ccfs_norm.-ccf_template,dims=1)
    ccf_template = vec(sum(ccfs_norm./est_ccf_vars_norm,dims=2) ./sum(1.0 ./est_ccf_vars_norm,dims=2))
    return ccf_template
end

function calc_ccf_template( ccfs::A2, ccf_vars::A3; assume_normalized::Bool = false  ) where { T2<:Real, A2<:AbstractArray{T2,2}, T3<:Real, A3<:AbstractArray{T3,2}  }
    (ccfs_norm, ccf_vars_norm) = assume_normalized ? (ccfs, ccf_vars) : calc_normalized_ccfs(ccfs, ccf_vars)
    ccf_template = vec(sum(ccfs./ccf_vars,dims=2) ./sum(1.0 ./ccf_vars,dims=2))
    return ccf_template
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs