text:
using DynamicGrids, Setfield, Test, LinearAlgebra, StaticArrays, OffsetArrays

using DynamicGrids.Neighborhoods

import DynamicGrids.Neighborhoods: _window, hoodsize, radius

init = [0 0 0 1 1 1
        1 0 1 1 0 1
        0 1 1 1 1 1
        0 1 0 0 1 0
        0 0 0 0 1 1
        0 1 0 1 1 0]

win1 = [0 0 0
        0 1 0
        0 0 0]
win2 = [1 1 1
        1 0 1
        1 1 1]
win3 = [1 1 1
        0 0 1
        0 0 1]

@testset "Moore" begin
    moore = Moore{1}(init[1:3, 1:3])
    @test _window(moore) == init[1:3, 1:3]
    @test hoodsize(moore) == 3
    @test moore[1] == 0
    @test length(moore) == 8
    @test eltype(moore) == Int
    @test neighbors(moore) isa Tuple
    @test collect(neighbors(moore)) == [0, 1, 0, 0, 1, 0, 1, 1]
    @test sum(moore) == sum(neighbors(moore)) == 4
    @test offsets(moore) ==
        ((-1, -1), (0, -1), (1, -1), (-1, 0), (1, 0), (-1, 1), (0, 1), (1, 1))

    moore1 = @set moore._window = win1
    @test moore1._window == win1
    moore2 = DynamicGrids.setwindow(moore, win2)
    @test moore2._window == win2

    @test sum(Moore{1}(win1)) == 0
    @test sum(Moore{1}(win2)) == 8
    @test sum(Moore{1}(win3)) == 5
end

@testset "Window" begin
    @test Window{1}() == Window(1) == Window(zeros(3, 3))
    window = Window{1}(init[1:3, 1:3])
    @test _window(window) == init[1:3, 1:3]
    @test hoodsize(window) == 3
    @test window[1] == 0
    @test window[2] == 1
    @test length(window) == 9
    @test eltype(window) == Int
    @test neighbors(window) isa Array
    @test neighbors(window) == _window(window)
    @test sum(window) == sum(neighbors(window)) == 4
    @test offsets(window) == ((-1, -1), (0, -1), (1, -1), (-1, 0), (0, 0),
                              (1, 0), (-1, 1), (0, 1), (1, 1))

    window1 = @set window._window = win1
    @test window1._window == win1
    window2 = DynamicGrids.setwindow(window, win2)
    @test window2._window == win2

    @test sum(Window{1}(win1)) == 1
    @test sum(Window{1}(win2)) == 8
    @test sum(Window{1}(win3)) == 5
end

@testset "VonNeumann" begin
    vonneumann = VonNeumann(1, init[1:3, 1:3])
    @test offsets(vonneumann) == ((0, -1), (-1, 0), (1, 0), (0, 1))
    @test _window(vonneumann) == init[1:3, 1:3]
    @test hoodsize(vonneumann) == 3
    @test vonneumann[1] == 1
    @test vonneumann[2] == 0
    @test length(vonneumann) == 4
    @test eltype(vonneumann) == Int
    @test neighbors(vonneumann) isa Tuple
    @test collect(neighbors(vonneumann)) == [1, 0, 1, 1]
    @test sum(neighbors(vonneumann)) == 3
    vonneumann2 = VonNeumann(2)
    @test offsets(vonneumann2) ==
       ((0, -2), (-1, -1), (0, -1), (1, -1),
         (-2 , 0), (-1, 0), (1, 0), (2, 0),
         (-1, 1), (0, 1), (1, 1), (0, 2))

    vonneumann1 = @set vonneumann._window = win2
    @test vonneumann1._window == win2
    vonneumann2 = DynamicGrids.setwindow(vonneumann, win3)
    @test vonneumann2._window == win3

    @test sum(VonNeumann(1, win1)) == 0
    @test sum(VonNeumann(1, win2)) == 4
    @test sum(VonNeumann(1, win3)) == 2
end

@testset "Positional" begin
    win = [0 1 0 0 1
           0 0 1 0 0
           0 0 0 1 1
           0 0 1 0 1
           1 0 1 0 1]
    custom1 = Positional(((-1,-1), (2,-2), (2,2), (-1,2), (0,0)), win)
    custom2 = Positional{((-1,-1), (0,-1), (1,-1), (2,-1), (0,0))}(win)
    layered = LayeredPositional(
        (Positional((-1,1), (-2,2)), Positional((1,2), (2,2))), win)

    @test neighbors(custom1) isa Tuple
    @test collect(neighbors(custom1)) == [0, 1, 1, 0, 0]

    @test sum(custom1) == 2
    @test sum(custom2) == 0
    @test sum(layered) == (1, 2)
    @test offsets(layered) == (((-1, 1), (-2, 2)), ((1, 2), (2, 2)))

    layered1 = @set layered._window = 2win
    @test layered1._window == 2win
    layered2 = DynamicGrids.setwindow(layered, 3win)
    @test layered2._window == 3win

    win = reshape(2:10, 3, 3)
    hood = Positional(((-1, -1), (1, 1)), win)
    @test neighbors(hood) == (2, 10)
    @test offsets(hood) == ((-1, -1), (1, 1))
    @test positions(hood, (2, 2)) == ((1, 1), (3, 3))
end

@testset "LayeredPositional" begin
    lhood = LayeredPositional(
        Positional(((-1, -1), (1, 1)), ), Positional(((-2, -2), (2, 2)), )
    )
    @test offsets(lhood) == (((-1, -1), (1, 1)), ((-2, -2), (2, 2)))
    @test collect.(collect(positions(lhood, (1, 1)))) ==
        [[(0, 0), (2, 2)],
         [(-1, -1), (3, 3)]]

    win = reshape(1:25, 5, 5)
    lhood_win = DynamicGrids.setwindow(lhood, win)
    @test lhood_win._window == lhood_win.layers[1]._window ===
          lhood_win.layers[2]._window === win
      lhood_win.layers[2]._window
    @test map(radius, lhood_win.layers) == (2, 2)
    @test neighbors(lhood_win) == ((7, 19), (1, 25))
end

@testset "Kernel" begin
    win = reshape(1:9, 3, 3)
    @testset "Window" begin
        mat = zeros(3, 3)
        @test Kernel(mat) == Kernel(Window(1), mat)
        @test_throws ArgumentError Kernel(Window(2), mat)
        k = Kernel(Window{1,2,9,typeof(win)}(win), SMatrix{3,3}(reshape(1:9, 3, 3)))
        @test kernelproduct(k) == sum((1:9).^2)
        @test neighbors(k) == reshape(1:9, 3, 3)
        @test offsets(k) == ((-1, -1), (0, -1), (1, -1), (-1, 0), (0, 0),
                             (1, 0), (-1, 1), (0, 1), (1, 1))
        @test positions(k, (2, 2)) == ((1, 1), (2, 1), (3, 1), (1, 2),
                                       (2, 2), (3, 2), (1, 3), (2, 3), (3, 3))
    end
    @testset "Moore" begin
        k = Kernel(Moore{1,2,8,typeof(win)}(win), (1:4..., 6:9...))
        @test kernelproduct(k) == sum((1:4).^2) + sum((6:9).^2)
    end
    @testset "Positional" begin
        off = ((0,-1),(-1,0),(1,0),(0,1))
        hood = Positional{off,1,2,4,typeof(win)}(win)
        k = Kernel(hood, 1:4)
        @test kernelproduct(k) == 1 * 2 + 2 * 4 + 3 * 6 + 4 * 8
    end
end

@testset "neighbors works on rule" begin
    rule = Life(;neighborhood=Moore{1}([0 1 1; 0 0 0; 1 1 1]))
    @test sum(neighbors(rule)) == 5
end

@testset "readwindow" begin
    grid1 = [0, 1, 2, 3, 4, 0]
    grid2 = [
         0  0  0  0  0  0
         0  1  2  3  4  0
         0  5  6  7  8  0
         0  9 10 11 12  0
         0  0  0  0  0  0
    ]
    @test DynamicGrids.readwindow(Moore{1,1}(), grid1, (2,)) == [0, 1, 2]
    @test DynamicGrids.readwindow(Moore{1,1}(), OffsetArray(grid1, (0:5)), (2,)) == [1, 2, 3]
    @test_throws DimensionMismatch DynamicGrids.readwindow(Moore{1,1}(), grid2, (2,))
    @test_throws DimensionMismatch DynamicGrids.readwindow(Moore{1,2}(), grid2, (2,))
    @test DynamicGrids.readwindow(Moore{1,2}(), grid2, (2, 2)) == [0 0 0; 0 1 2; 0 5 6]
    @test DynamicGrids.readwindow(Moore{1,2}(), OffsetArray(grid2, (0:4, 0:5)), (2, 2)) == [1 2 3; 5 6 7; 9 10 11]
end

@testset "pad/unpad axes" begin
    A = zeros(6, 7)
    @test pad_axes(A, 2) == (-1:8, -1:9)
    @test pad_axes(A, Moore(3)) == (-2:9, -2:10)
    @test unpad_axes(A, 2) == (3:4, 3:5)
    @test unpad_axes(A, VonNeumann(1)) == (2:5, 2:6)
end


text:
"""
Genie utilities for working over HTTP.
"""
module HTTPUtils
using DocStringExtensionsMock

import HTTP

"""
$TYPEDSIGNATURES

Converts a `HTTP.Request` to a `Dict`.
"""
function Base.Dict(req::HTTP.Request) :: Dict{String,String}
  result = Dict{String,String}()
  for (k,v) in Dict(req.headers)
    result[lowercase(string(k))] = lowercase(string(v))
  end

  result
end

end

ext: jl


text:
module SimpleCG

include("basic_types.jl")
export Vector3, Ray3, Canvas, ColorRGB,
    Black, White, Red, Green, Blue,
    imshow

include("config.jl")
export LightDirection, LightColor

include("materials.jl")
export Material, sample,
    CheckerMaterial, PhongMaterial

include("geometries.jl")
export Geometry, IntersectResult, hit,
    UnionGeometry, GeometryWithMaterial,
    Sphere, Plane

include("cameras.jl")
export Camera, right, up, generate_ray,
    PerspectiveCamera

include("lights.jl")
export Light, Lights, LightSample, InShadow, in_shadow,
    DirectionalLight, PointLight, SpotLight

include("render.jl")
export render!

end


text:
export getSensMatVec

function getSensMatVec(v::Vector,m::Vector,pFor::FWIparam)

    # extract pointers
    M    			= pFor.Mesh
    omega 			= pFor.omega
    gamma 			= pFor.gamma
    Q     			= pFor.Sources
    P     			= pFor.Receivers
    Ainv   			= pFor.ForwardSolver
    nsrc 			= size(Q,2)
    nrec 			= size(P,2)
	batchSize 		= pFor.forwardSolveBatchSize;
	select  		= pFor.sourceSelection;

	if length(select) > 0
		nsrc = length(select);
	end

	if batchSize > nsrc
		batchSize = nsrc;
	end

	numBatches 	= ceil(Int64,nsrc/batchSize);

	# derivative of mass matrix
	An2cc = getNodalAverageMatrix(M);
	## ALL AT ONCE CODE
	H = spzeros(ComplexF64,prod(M.n.+1),prod(M.n.+1));
	gamma_t = An2cc'*gamma;
	m_t = An2cc'*m;
	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		H = GetHelmholtzOperator(M,m_t,omega, gamma_t, true,useSommerfeldBC);
		Ainv.helmParam = HelmholtzParam(M,gamma_t,m_t,omega,true,useSommerfeldBC);
		H = H + GetHelmholtzShiftOP(m_t, omega,Ainv.shift[1]);
		H = sparse(H');
		# H is actually shifted laplacian now...
	elseif isa(Ainv,JuliaSolver)
		H = GetHelmholtzOperator(M,m_t,omega, gamma_t, true,useSommerfeldBC);
	end
	Jv = zeros(ComplexF64,nrec,nsrc);
	t = ((1.0.-1im*vec(gamma_t./omega)).*(An2cc'*v));
	for k_batch = 1:numBatches
		batchIdxs = (k_batch-1)*batchSize + 1 : min(k_batch*batchSize,nsrc);
		if pFor.useFilesForFields
			filename = getFieldsFileName(omega);
			file     = matopen(filename);
			U        = read(file,string("Ubatch_",k_batch));
			close(file);
		else
			U = pFor.Fields[:,batchIdxs];
		end
		U = t.*U;
		U,Ainv = solveLinearSystem(H,U,Ainv,0);
		Jv[:,batchIdxs] = omega^2*(P'*U);
	end

    return vec(Jv)
end


text:
using PyCall
using Requests
using domainUtils

# set-up
url="http://lp-sandbox.ekstep.org:8080/taxonomy-service/v2/analytics/domain/map"
baseURL = "http://lp-sandbox.ekstep.org:8080/taxonomy-service/v2/content/"

# connect to neo4j via pycall

# connect to cassandra
@pyimport cassandra.cluster as cluster
REcluster=cluster.Cluster()

# connect to neo4j
@pyimport neo4jrestclient.client as n4c
gdb = n4c.GraphDatabase("http://localhost:7474/db/data/", username="neo4j", password="soma")
n = gdb.nodes.create(color="Red", widht=16, height=32)

get("http://httpbin.org/get"; query = Dict("title" => "page1"))
url="http://localhost:7474/user/neo4j/password"
resp = get(url,Dict("neo4j" => "soma"))

# using pycall and py2neo

@pyimport py2neo

max_stars_repo_name: harshavardhanc/sunbird-analytics

max_issues_repo_name: harshavardhanc/sunbird-analytics

max_forks_repo_name: harshavardhanc/sunbird-analytics


text:
using ImageFiltering
using Images
using Colors
using Luxor
using ColorTypes: ARGB32

function drawimagematrix(m)
    d = Drawing(500, 500, "/home/src/Projects/experiments/test.png")
    origin()
    circle(O, 170, :clip)
    w, h = size(m)
    t = Tiler(500, 500, w, h)
    mi = convertmatrixtocolors(m)
    for (pos, n) in t
        c = mi[t.currentrow, t.currentcol]
        setcolor(c)
        box(pos, t.tilewidth + 1, t.tileheight + 1, :fill)
    end
    sethue("black")
    setline(1)
    clipreset()
    circle(O, 170, :stroke)
    setline(2)
    radius = 15
    electrode_locations = [
        O,
        Point(-70, 0),
        Point(70, 0),
        Point(-140, 0),
        Point(140, 0),
        Point(0, 70),
        Point(-50, 70),
        Point(50, 70),
        Point(0, -70),
        Point(-50, -70),
        Point(50, -70),
        Point(115, -80),
        Point(-115, -80),
        Point(115, 80),
        Point(-115, 80),
        Point(40, -135),
        Point(-40, -135),
        Point(-190, -10),
        Point(190, -10),
        Point(-40, 135),
        Point(40, 135),
    ]

    electrode_names = [
        "Cz",
        "C3",
        "C4",
        "T3",
        "T4",
        "Pz",
        "P3",
        "P4",
        "Fz",
        "F3",
        "F4",
        "F8",
        "F7",
        "T6",
        "T5",
        "Fp2",
        "Fp1",
        "A1",
        "A2",
        "O1",
        "O2",
    ]
    electrode.(electrode_locations, "white", "black", :fill, radius, electrode_names)

    finish()
    return d
end

function convertmatrixtocolors(m)
    return convert.(Colors.RGBA, m)
end

function electrode(
    p = O,
    fill_color = "white",
    outline_color = "black",
    action = :fill,
    radius = 25,
    circ_text = "",
)
    sethue(fill_color)
    circle(p, radius, :fill)
    sethue(outline_color)
    circle(p, radius, :stroke)
    text(circ_text, p, valign = :middle, halign = :center)
end

test = @imagematrix begin
    background("white")
    t = Tiler(200, 200, 200, 200)
    for (pos, n) in t
        sethue(HSB(360 * noise(t.currentrow * 0.009, t.currentcol * 0.009), 1, 1))
        box(pos, 1, 1, :fillstroke)
    end
end 150 150

img = drawimagematrix(test)

#=

        Uncomment to demonstrate creating synthetic data.
        Can be an end point for actual EEG data in the future

=#

# I = [ARGB32(1, 1, rand(), 1) for _ = 1:30, _ = 1:30]
# for val in 1:length(I)
    # if rand([true, false])
        # I[val] = ARGB32(1, 0, 1, 1)
    # end
# end

# img = imfilter(I, Kernel.gaussian(2))

# drawimagematrix(img)

max_stars_repo_name: TheCedarPrince/NeuriViz

max_issues_repo_name: TheCedarPrince/NeuriViz

max_forks_repo_name: TheCedarPrince/NeuriViz


text:
#-----------------------------------------------------------------------# General
const Tup = Union{Tuple, NamedTuple}
const XY = Union{AbstractVector, Tup, Pair}
const VectorOb = Union{AbstractVector, Tup}
const TwoThings{T,S} = Union{Tuple{T,S}, Pair{T,S}, NamedTuple{names, Tuple{T,S}}} where names

_dot(x::AbstractVector, y::AbstractVector) = dot(x, y)
@generated function _dot(x::VectorOb, y::VectorOb)
    n = length(fieldnames(x))
    quote
        out = zero(promote_type(typeof(x[1]), typeof(y[1])))
        Base.Cartesian.@nexprs $n i-> (out += x[i] * y[i])
        out
    end
end

const ϵ = 1e-7

#-----------------------------------------------------------------------# BiasVec
"""
    BiasVec(x)

Lightweight wrapper of a vector which adds a "bias" term at the end.

# Example

    BiasVec(rand(5))
"""
struct BiasVec{T, A <: VectorOb} <: AbstractVector{T}
    x::A
    bias::T
end
BiasVec(x::AbstractVector{T}) where {T} = BiasVec(x, one(T))
BiasVec(x::Tup) = BiasVec(x, one(typeof(first(x))))

Base.length(v::BiasVec) = length(v.x) + 1
Base.size(v::BiasVec) = (length(v), )
Base.getindex(v::BiasVec, i::Int) = i > length(v.x) ? v.bias : v.x[i]
Base.IndexStyle(::Type{<:BiasVec}) = IndexLinear()

#-----------------------------------------------------------------------# fit!
@deprecate fit!(o::OnlineStat{VectorOb}, x::AbstractMatrix) fit!(o, eachrow(x))
@deprecate fit!(o::OnlineStat{XY}, x::AbstractMatrix) fit!(o, eachrow(x))
@deprecate fit!(o::OnlineStat{XY}, xy::Tuple{AbstractMatrix, AbstractVector}) fit!(o, zip(eachrow(xy[1]), xy[2]))

# #-----------------------------------------------------------------------# sparkline
# const ticks = ['▁','▂','▃','▄','▅','▆','▇']

# sparkline(x) = sparkline(stdout, x)
# sparkline(io::IO, x) = print(io, sparkstring(x))

# function sparkstring(x)
#     out = ""
#     if !isempty(x)
#         a, b = extrema(x)
#         for xi in x
#             i = round(Int, 7 * (xi-a) / (b-a))
#             out *= i == 0 ? ' ' : ticks[i]
#         end
#     end
#     out
# end


text:
# FASTQ Record
# ============

mutable struct Record
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # indexes
    identifier::UnitRange{Int}
    description::UnitRange{Int}
    sequence::UnitRange{Int}
    quality::UnitRange{Int}
end

"""
    FASTQ.Record()

Create an unfilled FASTQ record.
"""
function Record()
    return Record(UInt8[], 1:0, 1:0, 1:0, 1:0, 1:0)
end

"""
    FASTQ.Record(data::Vector{UInt8})

Create a FASTQ record object from `data`.

This function verifies and indexes fields for accessors.
Note that the ownership of `data` is transferred to a new record object.
"""
function Record(data::Vector{UInt8})
    return convert(Record, data)
end

function Base.convert(::Type{Record}, data::Vector{UInt8})
    record = Record(data, 1:0, 1:0, 1:0, 1:0, 1:0)
    index!(record)
    return record
end

"""
    FASTQ.Record(str::AbstractString)

Create a FASTQ record object from `str`.

This function verifies and indexes fields for accessors.
"""
function Record(str::AbstractString)
    return Record(Vector{UInt8}(str))
end

"""
    FASTQ.Record(identifier, sequence, quality; offset=33)

Create a FASTQ record from `identifier`, `sequence` and `quality`.
"""
function Record(identifier::AbstractString, sequence, quality::Vector; offset=33)
    return Record(identifier, nothing, sequence, quality, offset=offset)
end

"""
    FASTQ.Record(identifier, description, sequence, quality; offset=33)

Create a FASTQ record from `identifier`, `description`, `sequence` and `quality`.
"""
function Record(identifier::AbstractString, description::Union{AbstractString,Nothing}, sequence, quality::Vector; offset=33)
    if length(sequence) != length(quality)
        throw(ArgumentError("the length of sequence doesn't match the length of quality"))
    end
    buf = IOBuffer()
    print(buf, '@', identifier)
    if description != nothing
        print(buf, ' ', description)
    end
    print(buf, '\n')
    print(buf, sequence, '\n')
    print(buf, "+\n")
    ascii_quality = convert(Vector{UInt8}, quality .+ offset)
    write(buf, ascii_quality, '\n')
    return Record(take!(buf))
end

function Base.convert(::Type{Record}, str::AbstractString)
    return convert(Record, convert(Vector{UInt8}, str))
end

function Base.copy(record::Record)
    return Record(
        record.data[record.filled],
        record.filled,
        record.identifier,
        record.description,
        record.sequence,
        record.quality)
end

function Base.write(io::IO, record::Record)
    return unsafe_write(io, pointer(record.data, first(record.filled)), length(record.filled))
end

function Base.print(io::IO, record::Record)
    write(io, record)
    return nothing
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   identifier: ", hasidentifier(record) ? identifier(record) : "<missing>")
        println(io, "  description: ", hasdescription(record) ? description(record) : "<missing>")
        println(io, "     sequence: ", hassequence(record) ? sequence(String, record) : "<missing>")
          print(io, "      quality: ", hasquality(record) ? quality(record) : "<missing>")
    else
        print(io, " <not filled>")
    end
end

function initialize!(record::Record)
    record.filled = 1:0
    record.identifier = 1:0
    record.description = 1:0
    record.sequence = 1:0
    record.quality = 1:0
    return record
end

function BioCore.isfilled(record::Record)
    return !isempty(record.filled)
end

# Accessor functions
# ------------------

"""
    identifier(record::Record)::String

Get the sequence identifier of `record`.
"""
function identifier(record::Record)::String
    checkfilled(record)
    if !hasidentifier(record)
        missingerror(:identifier)
    end
    return String(record.data[record.identifier])
end

"""
    hasidentifier(record::Record)

Checks whether or not the `record` has an identifier.
"""
function hasidentifier(record::Record)
    return isfilled(record)
end

"""
    description(record::Record)::String

Get the description of `record`.
"""
function description(record::Record)::String
    checkfilled(record)
    if !hasdescription(record)
        missingerror(:description)
    end
    return String(record.data[record.description])
end

"""
    hasdescription(record::Record)

Checks whether or not the `record` has a description.
"""
function hasdescription(record)
    return isfilled(record) && record.description != 1:0
end

"""
    sequence(::Type{S}, record::Record, [part::UnitRange{Int}])::S

Get the sequence of `record`.

`S` can be either a subtype of `BioSequences.Sequence` or `String`.
If `part` argument is given, it returns the specified part of the sequence.
"""
function sequence(::Type{S}, record::Record, part::UnitRange{Int}=1:lastindex(record.sequence))::S where S <: BioSequences.Sequence
    checkfilled(record)
    if !hassequence(record)
        missingerror(:sequence)
    end
    seqpart = record.sequence[part]
    return S(record.data, first(seqpart), last(seqpart))
end

function sequence(::Type{String}, record::Record, part::UnitRange{Int}=1:lastindex(record.sequence))::String
    checkfilled(record)
    if !hassequence(record)
        missingerror(:sequence)
    end
    return String(record.data[record.sequence[part]])
end

"""
    sequence(record::Record, [part::UnitRange{Int}])::BioSequences.DNASequence

Get the sequence of `record`.
"""
function sequence(record::Record, part::UnitRange{Int}=1:lastindex(record.sequence))::BioSequences.DNASequence
    checkfilled(record)
    return sequence(BioSequences.DNASequence, record, part)
end

"""
    hassequence(record::Record)

Checks whether or not a sequence record contains a sequence.
"""
function hassequence(record::Record)
    # zero-length sequence may exist
    return isfilled(record)
end

"""
    quality(record::Record, [offset::Integer=33, [part::UnitRange]])::Vector{UInt8}

Get the base quality of `record`.
"""
function quality(record::Record, offset::Integer=33, part::UnitRange{Int}=1:lastindex(record.quality))::Vector{UInt8}
    checkfilled(record)
    quality = record.data[record.quality[part]]
    for i in 1:lastindex(part)
        # TODO: Checked arithmetic?
        @inbounds quality[i] -= offset
    end
    return quality
end

"""
    quality(record::Record, encoding_name::Symbol, [part::UnitRange])::Vector{UInt8}

Get the base quality of `record` by decoding with `encoding_name`.

The `encoding_name` can be either `:sanger`, `:solexa`, `:illumina13`, `:illumina15`, or `:illumina18`.
"""
function quality(record::Record, encoding_name::Symbol, part::UnitRange{Int}=1:lastindex(record.quality))::Vector{UInt8}
    checkfilled(record)
    encoding = (
        encoding_name == :sanger     ?     SANGER_QUAL_ENCODING :
        encoding_name == :solexa     ?     SOLEXA_QUAL_ENCODING :
        encoding_name == :illumina13 ? ILLUMINA13_QUAL_ENCODING :
        encoding_name == :illumina15 ? ILLUMINA15_QUAL_ENCODING :
        encoding_name == :illumina18 ? ILLUMINA18_QUAL_ENCODING :
        throw(ArgumentError("quality encoding ':$(encoding_name)' is not supported")))
    quality = Vector{UInt8}(undef, length(part))
    if !isempty(part)
        qpart = record.quality[part]
        check_quality_string(encoding, record.data, first(qpart), last(qpart))
        decode_quality_string!(encoding, record.data, quality, first(qpart), last(qpart))
    end
    return quality
end

"""
    hasquality(record::Record)

Check whether the given FASTQ `record` has a quality string.
"""
function hasquality(record::Record)
    return isfilled(record)
end

function BioCore.seqname(record::Record)
    return identifier(record)
end

function BioCore.hasseqname(record::Record)
    return hasidentifier(record)
end

function BioCore.sequence(record::Record)
    return sequence(record)
end

function BioCore.sequence(::Type{S}, record::Record) where S <: BioSequences.Sequence
    return sequence(S, record)
end

function BioCore.hassequence(record::Record)
    return hassequence(record)
end

# Helper functions
# ----------------

function checkfilled(record)
    if !isfilled(record)
        throw(ArgumentError("unfilled FASTQ record"))
    end
end


text:
"""
    logistic2_bidir(u₀, c_xy, c_yx, r₁, r₂, σ_xy, σ_yx)

Equations of motion for a bidirectional logistic model for the chaotic
population dynamics of two interacting species. This system is from [1],
and is given by

```math
\\begin{align}
x(t+1) &= r_1 f_{yx}^{t}(1 - f_{yx}^{t}) \\
y(t+1) &= r_2 f_{xy}^{t}(1 - f_{xy}^{t}) \\
f_{xy}^t &= \\dfrac{y(t) + c_{xy}(x(t) + \\sigma_{xy} \\xi_{xy}^t )}{1 + c_{xy} (1 + \\sigma_{xy} )} \\
f_{yx}^t &= \\dfrac{x(t) + c_{yx}(y(t) + \\sigma_{yx} \\xi_{yx}^t )}{1 + c_{yx} (1 + \\sigma_{yx} )},
\\end{align}
```

where the coupling strength ``c_{xy}`` controls how strongly species ``x`` influences species
``y``, and vice versa for ``c_{yx}``. To simulate time-varying influence of unobserved
processes, we use the dynamical noise terms ``\\xi_{xy}^t`` and ``\\xi_{yx}^t``, drawn from a
uniform distribution with support on ``[0, 1]``. If ``\\sigma_{xy} > 0``, then the influence
of ``x`` on ``y`` is masked by dynamical noise equivalent to ``\\sigma_{xy} \\xi_{xy}^{t}`` at
the ``t``-th iteration of the map, and vice versa for ``\\sigma_{yx}``.

## References

1. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
function eom_logistic2_bidir(dx, x, p, n)

    # c_xy is the coupling from x to y
    # c_yx is the coupling from y to x
    # σ_yx is the dynamical noise from y to x
    # σ_xy is the dynamical noise from y to x
    c_xy, c_yx, r₁, r₂, σ_xy, σ_yx = (p...,)

    ξ₁ = rand() # random number from flat distribution on [0, 1]
    ξ₂ = rand() # random number from flat distribution on [0, 1]
    x, y = x[1], x[2]

    f_xy = (y +  c_xy*(x + σ_xy*ξ₁) ) / (1 + c_xy*(1+σ_xy))
    f_yx = (x +  c_yx*(y + σ_yx*ξ₂) ) / (1 + c_yx*(1+σ_yx))

    dx[1] = r₁ * (f_yx) * (1 - f_yx)
    dx[2] = r₂ * (f_xy) * (1 - f_xy)
    return
end

function logistic2_bidir(u₀, c_xy, c_yx, r₁, r₂, σ_xy, σ_yx)
    p = [c_xy, c_yx, r₁, r₂, σ_xy, σ_yx]
    DiscreteDynamicalSystem(eom_logistic2_bidir, u₀, p)
end

"""
    logistic2_bidir(;u₀ = rand(2), c_xy = 0.1, c_yx = 0.1,
        r₁ = 3.78, r₂ = 3.66, σ_xy = 0.05, σ_yx = 0.05)

A bidirectional logistic model for the chaotic population dynamics of two interacting
species [1].

## Equations of motion

The equations of motion are

```math
\\begin{align}
x(t+1) &= r_1 f_{yx}^{t}(1 - f_{yx}^{t}) \\\\
y(t+1) &= r_2 f_{xy}^{t}(1 - f_{xy}^{t}) \\\\
f_{xy}^t &= \\dfrac{y(t) + c_{xy}(x(t) + \\sigma_{xy} \\xi_{xy}^t )}{1 + c_{xy} (1 + \\sigma_{xy} )} \\\\
f_{yx}^t &= \\dfrac{x(t) + c_{yx}(y(t) + \\sigma_{yx} \\xi_{yx}^t )}{1 + c_{yx} (1 + \\sigma_{yx} )},
\\end{align}
```

where the coupling strength ``c_{xy}`` controls how strongly species ``x`` influences species
``y``, and vice versa for ``c_{yx}``. To simulate time-varying influence of unobserved
processes, we use the dynamical noise terms ``\\xi_{xy}^t`` and ``\\xi_{yx}^t``, drawn from a
uniform distribution with support on ``[0, 1]``. If ``\\sigma_{xy} > 0``, then the influence
of ``x`` on ``y`` is masked by dynamical noise equivalent to ``\\sigma_{xy} \\xi_{xy}^{t}`` at
the ``t``-th iteration of the map, and vice versa for ``\\sigma_{yx}``.

## References

1. Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation
    using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212.
"""
logistic2_bidir(;u₀ = rand(2), c_xy = 0.1, c_yx = 0.1,
    r₁ = 3.78, r₂ = 3.66, σ_xy = 0.05, σ_yx = 0.05) =
    logistic2_bidir(u₀, c_xy, c_yx, r₁, r₂, σ_xy, σ_yx)


text:
IMGUI_DEMO_BUTTON_STATES = Dict(
    "basic_button" => OnOffButtons.State("Button"),
    "basic_button2" => Buttons.State("Button2"),
)

IMGUI_DEMO_RADIO_BUTTON_GROUP_STATES =
    RadioButtonGroups.State("radio button group",
                            [RadioButtons.State("radio a"),
                             RadioButtons.State("radio b"),
                             RadioButtons.State("radio c")])

IMGUI_DEMO_COLOR_BUTTON_STATES = [
    ColorButtons.State("click##1"),
    ColorButtons.State("click##2"),
    ColorButtons.State("click##3"),
    ColorButtons.State("click##4"),
    ColorButtons.State("click##5"),
    ColorButtons.State("click##6"),
    ColorButtons.State("click##7"),
]


text:
module problem110

description = """
Diophantine reciprocals II

In the following equation x, y, and n are positive integers.
1/x + 1/y = 1/n

It can be verified that when n = 1260 there are 113 distinct solutions and this is the least value of n for which the total number of distinct solutions exceeds one hundred.
What is the least value of n for which the number of distinct solutions exceeds four million?
"""

# From the forums in problem 108, I found that the number of solutions to N = ((number of divisors of N^2) + 1)/2
# Number of divisors can be efficiently calculated by taking the product of (power+1) for the powers of the prime factors of n

function nextprime(n)
  m = n+1
  while !isprime(m)
    m += 1
  end
  return m
end

# solution_count(N) = (sigma0(N^2) + 1) / 2
# N has be factored into product of primes already
# Sigma0(N^2) is inlined by doubling the power of each prime
function solution_count(d :: Dict)
  return div(prod([2*v + 1 for v in values(d)]) + 1, 2)
end

# Seach for the least N solving 1/x + 1/y = 1/N with at least the given number of solutions
# Iteratively search using numbers of the form 2^p2 * 3^p3 * 5^p5 ..., where p2 >= p3 >= p5 etc..
# Stop searching once increasing the power on the 2 term doesn't reduce the minimum solution
function least_n_with_num_solutions(numSolutions :: Integer)
  leastN = -1
  p2 = 1
  while true
    n = solution_search(numSolutions, [2 => p2])
    if leastN < 0 || n < leastN
      leastN = n
    elseif n > leastN
      break
    end
    p2 += 1
  end
  return leastN
end

# Workhorse function, searches for value with > numSolutions
# value is kept in factorized form in the Dict d, but is returned as BigInt
function solution_search(numSolutions :: Integer, d :: Dict)

  # base case: D already has enough solutions, so return the value as a BigInt
  if solution_count(d) >= numSolutions
    return prod([BigInt(prime)^power for (prime, power) in d])
  end

  # Recusive case: Find the next prime to add to d, and check with all powers up to the power of the current max prime in d.
  leastN = -1
  maxPrime = maximum(keys(d))
  p = nextprime(maxPrime)
  for exponent = 1 : d[maxPrime]
    d[p] = exponent
    n = solution_search(numSolutions, d)
    if leastN < 0 || (n > 0 && n < leastN)
      leastN = n
    end
  end

  # Cleanup the Dict on the way out
  delete!(d, p)
  return leastN
end

using Base.Test

@test solution_count(factor(1260)) == 113
@test least_n_with_num_solutions(113) == 1260
@test least_n_with_num_solutions(1000) == 180180 #from problem 108

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
fname = :ARCMA57_0
c = Combi(hessian_sparse,PDataMA57_0,solve_modelARCDiag,preprocessMA57_0,decreaseFact,Tparam())
include("Template.jl")


text:
function pointmin(_I :: Array{Float64,2})
    if isempty(_I)
        return false
    end
    I = copy(_I)
    Fx=zeros(size(I))
    Fy=zeros(size(I))
    Fz=zeros(size(I))

    J=zeros( tuple([size(I,i)+2 for i = 1:ndims(I)]...) )
    # J[:,:]=maximum(I[:])
    J .= maximum(I[:])

    J[2:end-1,2:end-1]=I
    Ne = [-1 -1; -1  0; -1  1; 0 -1; 0  1; 1 -1;  1  0; 1  1]
    for i=1:maximum(size(Ne))
        In=J[2+Ne[i,1]:end-1+Ne[i,1],2+Ne[i,2]:end-1+Ne[i,2]]
        check = In.<I #TODO Check this is what is intended (In<I)
        I[check]= copy(In[check])
        D=Ne[i,:]
        D=D./sqrt(sum(D.^2))
        Fx[check] .= D[1]
        Fy[check] .= D[2]
    end

    return (Fx,Fy)
end

function pointmin(_I :: Array{Float64,3})
    I = copy(_I)
    Fx=zeros(size(I))
    Fy=zeros(size(I))
    Fz=zeros(size(I))

    J=zeros( tuple([size(I,i)+2 for i = 1:ndims(I)]...) )
    J[:,:]=maximum(I[:])

    J[2:end-1,2:end-1,2:end-1]=I
    Ne=[-1 -1 -1; -1 -1  0; -1 -1  1; -1  0 -1; -1  0  0; -1  0  1; -1  1 -1; -1  1  0; -1  1  1;
    0 -1 -1;  0 -1  0;  0 -1  1;  0  0 -1;            0  0  1;  0  1 -1;  0  1  0;  0  1  1;
        1 -1 -1;  1 -1  0;  1 -1  1;  1  0 -1;  1  0  0;  1  0  1;  1  1 -1;  1  1  0;  1  1  1];

    for i=1:maximum(size(Ne))
        In = J[2+Ne[i,1]:end-1+Ne[i,1],2+Ne[i,2]:end-1+Ne[i,2],2+Ne[i,3]:end-1+Ne[i,3]]
        check = In.<I
        I[check]= In[check]
        D=Ne[i,:]
        D=D./sqrt(sum(D.^2))
        Fx[check]= D[1]
        Fy[check]= D[2]
        Fz[check]= D[3]
    end

    return (Fx,Fy,Fz)
end


text:
const F_NODE = 0
const C_NODE = 1
const U_NODE = 2

struct RS
end

function remove_diag!(a)
	n = size(a, 1)
	for i = 1:n
		for j in nzrange(a, i)
			if a.rowval[j] == i
				a.nzval[j] = 0
			end
   		end
	end
	dropzeros!(a)
end

function (::RS)(S)
	remove_diag!(S)
	RS_CF_splitting(S, adjoint(S))
end

function RS_CF_splitting(S::SparseMatrixCSC, T::SparseMatrixCSC)

	n = size(S,1)

	lambda = zeros(Int, n)

	Tp = T.colptr
	Tj = T.rowval
	Sp = S.colptr
	Sj = S.rowval

	interval_ptr = zeros(Int, n+1)
	interval_count = zeros(Int, n+1)
	index_to_node = zeros(Int,n)
	node_to_index = zeros(Int,n)

    for i = 1:n
		# compute lambda[i] - the number of nodes strongly coupled to node i
		lambda[i] = Sp[i+1] - Sp[i]
        interval_count[lambda[i] + 1] += 1
    end

	# initial interval_ptr
	@views accumulate!(+, interval_ptr[2:end], interval_count[1:end-1])

	# sort the nodes by the number of nodes strongly coupled to them:
	#   `index_to_node[idx]` is the node with `idx`th least number of nodes coupled to it
	#   `node_to_index[idx]` is the position of the `idx`th node in `index_to_node`
	# linear time and allocation-free equivalent to:
	#   sortperm!(index_to_node, lambda)
	#   node_to_index[index_to_node] .= 1:n
	interval_count .= 0 # temporarily zeroed, goes back to its original at end of loop
    for i = 1:n
        lambda_i = lambda[i] + 1
        interval_count[lambda_i] += 1
        index = interval_ptr[lambda_i] + interval_count[lambda_i]
        index_to_node[index] = i
        node_to_index[i]     = index
    end
	splitting = fill(U_NODE, n)

    # all nodes which no other nodes are strongly coupled to become F nodes
    for i = 1:n
		if lambda[i] == 0
            splitting[i] = F_NODE
		end
    end

	# i = index_to_node[top_index] can either refer to an F node or to the U node with the
	#	highest lambda[i].

	# index_to_node[interval_ptr[i]+1 : interval_ptr[i+1]] includes the set of U nodes with
	#	i-1 nodes strongly coupled to them, and other "inactive" F and C nodes.

	# C nodes are always in index_to_node[top_index:n]. So at the end of the last
	#	non-empty interval_ptr[i]+1 : interval_ptr[i+1] will be all the C nodes together
	#	with some other inactive F nodes.

	# when lambda_k+1 > lambda_i, i.e. lambda_k == lambda_i,  where lambda_k+1 = lambda_i+1
	#	is the new highest lambda[i], the line: `interval_ptr[lambda_k+1] = new_pos - 1`
	#	pushes the all the inactive C and F points to the end of the next now-non-empty
	#	interval.
	for top_index = n:-1:1
		i = index_to_node[top_index]
		lambda_i = lambda[i] + 1
		interval_count[lambda_i] -= 1

		splitting[i] == F_NODE && continue

		@assert splitting[i] == U_NODE
		splitting[i] = C_NODE
		for j in nzrange(S, i)
			row = S.rowval[j]
			if splitting[row] == U_NODE
				splitting[row] = F_NODE

				# increment lambda for all U nodes that node `row` is strongly coupled to
				for k in nzrange(T, row)
					rowk = T.rowval[k]

					if splitting[rowk] == U_NODE
						# to ensure `intervalcount` is inbounds
						lambda[rowk] >= n - 1 && continue

						# move rowk to the end of its current interval
						lambda_k = lambda[rowk] + 1
						old_pos  = node_to_index[rowk]
						new_pos  = interval_ptr[lambda_k] + interval_count[lambda_k]

						swap_node = index_to_node[new_pos]
						(index_to_node[old_pos], index_to_node[new_pos]) = (swap_node, rowk)
						node_to_index[rowk] = new_pos
						node_to_index[swap_node] = old_pos

						# increment lambda[rowk]
						lambda[rowk] += 1

						# update intervals
						interval_count[lambda_k]   -= 1
						interval_count[lambda_k+1] += 1
						interval_ptr[lambda_k+1]    = new_pos - 1
					end
				end
			end
		end

		# decrement lambda for all U nodes that node i is strongly coupled to
		for j in nzrange(T, i)
			row = T.rowval[j]
			if splitting[row] == U_NODE
				# to ensure `intervalcount` is inbounds
				lambda[row] == 0 && continue

				# move row to the beginning of its current interval
				lambda_j = lambda[row] + 1
				old_pos  = node_to_index[row]
				new_pos  = interval_ptr[lambda_j] + 1

				swap_node = index_to_node[new_pos]
				(index_to_node[old_pos], index_to_node[new_pos]) = (swap_node, row)
				node_to_index[row] = new_pos
				node_to_index[swap_node] = old_pos

				# decrement lambda[row]
				lambda[row] -= 1

				# update intervals
				interval_count[lambda_j]   -= 1
				interval_count[lambda_j-1] += 1
				interval_ptr[lambda_j]     += 1
			end
		end
	end
	splitting
end


text:
export ft_setdtr

function ft_setdtr(ft_handle::Culong)
  ft_status = ccall((:FT_SetDtr, d2xx),Cuint,(Culong,),ft_handle)
  checkstatus(ft_status)
  return nothing
end


text:
mutable struct DataP2G{N,T<:Real} <: AbstractData{N,T}
	pos::SVector{N,T}
	idx::Int64
	hsml::T
    mass::T
    field::T
	idx_ngbs::Vector{Int64} #ngb list for leaf node
end

DataP2G{N,T}() where {N,T} = DataP2G{N,T}(zero(SVector{N,T}),0,0,0,0,[])

function assign_additional_node_data!(n::DataP2G, old::DataP2G, new::DataP2G)
end

const MAX_DEPTH = Int(1000)

#return AMR structure (0 = leaf, 1=non-leaf node)
function get_AMRgrid(Tint::DataType, tree::Node{N,T,D}; max_depth::Int=MAX_DEPTH) where {N,T,D}
	gridAMR = Tint[]
	root_node_length = tree.length[1]
	get_AMRgrid_recursive!(gridAMR, tree, max_depth, root_node_length)
	return gridAMR
end

#default type is Int8 (to be compatible with RADMC3D)
function get_AMRgrid(tree::Node{N,T,D}; max_depth::Int=MAX_DEPTH) where {N,T,D}
	get_AMRgrid(Int8, tree, max_depth=max_depth)
end

#root_node_length is only used to figure out the node's depth
function get_AMRgrid_recursive!(gridAMR::Vector{Tint}, node::Node{N,T,D}, max_depth::Int, root_node_length::T) where {N,T,D, Tint<:Integer}
    if isLeaf(node)
        #println("in a leaf node")
		push!(gridAMR, 0)  #0 = leaf
    else
        #println("This is a node... ")
        #always open the node until we find a leaf
		depth = log2(round(root_node_length / node.length[1]))
		if depth < max_depth
			push!(gridAMR, 1)  #0 = leaf
	    	@inbounds for i in 1:2^N
	        	#println("open this node")
				get_AMRgrid_recursive!(gridAMR, node.child[i], max_depth, root_node_length)
	    	end
		else
			#don't open this node
			push!(gridAMR, 0)
		end
    end
end

#return grid/node volume following the order of the tree
function get_AMRgrid_volumes(tree::Node{N,T,D}; max_depth::Int=MAX_DEPTH) where {N,T,D}
	grid_volumes = T[]
	root_node_length = tree.length[1]
	get_AMRgrid_volumes_recursive!(grid_volumes, tree, max_depth, root_node_length)
	return grid_volumes
end

function get_AMRgrid_volumes_recursive!(volumearray, node::Node{N,T,D}, max_depth, root_node_length::T) where {N,T,D}
    if isLeaf(node)
		push!(volumearray, prod(node.length))
    else
		depth = log2(round(root_node_length / node.length[1]))
		if depth < max_depth
	    	@inbounds for i in 1:2^N
				get_AMRgrid_volumes_recursive!(volumearray, node.child[i], max_depth, root_node_length)
	    	end
		else
			push!(volumearray, prod(node.length))
		end
    end
end

#main function to map the particle information to an AMR grid

function map_particle_to_AMRgrid!(treeAMR::Node{N,T,D}, field::Vector{T}, volume::Vector{T}, X::Vector{SVector{N,T}}, hsml::Vector{T},
	boxsizes::SVector{N,T}, treeNgb::Union{Node{N,T,D}, Nothing}, serial::Bool) where {N,T,D}
	if serial
		map_particle_to_AMRgrid_recursive!(treeAMR, field, volume, X, hsml, boxsizes, treeNgb)
	else
		#4-layer unrolled (4^4 tasks)
		if isLeaf(treeAMR)
			map_particle_to_AMRgrid_recursive!(treeAMR, field, volume, X, hsml, boxsizes, treeNgb)
			return
		end
		@sync for i in 1:2^N
			if isLeaf(treeAMR.child[i])
				map_particle_to_AMRgrid_recursive!(treeAMR.child[i], field, volume, X, hsml, boxsizes, treeNgb)
			else
				for j in 1:2^N
					if isLeaf(treeAMR.child[i].child[j])
						map_particle_to_AMRgrid_recursive!(treeAMR.child[i].child[j], field, volume, X, hsml, boxsizes, treeNgb)
					else
						for k in 1:2^N
							if isLeaf(treeAMR.child[i].child[j].child[k])
								map_particle_to_AMRgrid_recursive!(treeAMR.child[i].child[j].child[k], field, volume, X, hsml, boxsizes, treeNgb)
							else
								for l in 1:2^N
									#@show i,j,k
									Threads.@spawn map_particle_to_AMRgrid_recursive!(treeAMR.child[i].child[j].child[k].child[l], field, volume, X, hsml, boxsizes, treeNgb)
								end
							end
						end
					end
				end
			end
		end

		#parallel version (1-layer unrolled)
		#@sync for i in 1:2^N
		#	Threads.@spawn map_particle_to_AMRgrid_recursive!(treeAMR.child[i], field, volume, X, hsml, boxsizes, treeNgb)
		#end
	end
end

function map_particle_to_AMRgrid_SPH!(treeAMR, field, volume, X, hsml, boxsizes;
	treeNgb::Union{Node{N,T,D}, Nothing}=nothing, serial::Bool=false) where {N,T,D}
	map_particle_to_AMRgrid!(treeAMR, field, volume, X, hsml, boxsizes, treeNgb, serial)
end
function map_particle_to_AMRgrid_MFM!(treeAMR, field, X, hsml, boxsizes;
	treeNgb::Union{Node{N,T,D}, Nothing}=nothing, serial::Bool=false) where {N,T,D}
	volume = eltype(field)[]
	map_particle_to_AMRgrid!(treeAMR, field, volume, X, hsml, boxsizes, treeNgb, serial)
end
function map_particle_to_AMRgrid_NGP!(treeAMR::Node{N,T,D}, field; serial::Bool=false) where {N,T,D}
	map_particle_to_AMRgrid!(treeAMR, field, T[], SVector{N,T}[], T[], SVector{N}(zeros(T,N)), nothing, serial)
end

#function map_particle_to_AMRgrid_thread_2nd!(tree::Node{N,T,D}, field::Vector{T}, volume::Vector{T},
#X::Vector{SVector{N,T}}, hsml::Vector{T}, node::Node{N,T,D}, boxsizes::SVector{N,T}) where {N,T,D}
#	@sync for i in 1:2^N, j in 1:2^N
#		Threads.@spawn map_particle_to_AMRgrid_recursive!(tree, field, volume, X, hsml, node.child[i].child[j], boxsizes)
#	end
#end

#recursively walk the tree
function map_particle_to_AMRgrid_loopP_recursive!(node::Node{N,T,D}, f_j::T, V_j::T, h_j::T, X_j::SVector{N,T}, boxsizes::SVector{N,T}) where {N,T,D}

	dx = nearest.(X_j - node.center, boxsizes)

	if !isLeaf(node)
		#@show "node"
		if all(abs.(dx) .< 0.5*node.length .+ h_j) #skip the node if no overlap
			@inbounds for i in 1:2^N
	        	map_particle_to_AMRgrid_loopP_recursive!(node.child[i], f_j, V_j, h_j, X_j, boxsizes)
			end
		end
	else
		#node can be used (either a leaf or a deep enough node)
		#println("leaf")
		dr = norm(dx)
		h_inv = 1 / h_j
		if dr < h_j
			Wij = kernel_cubic(dr*h_inv) * h_inv^N
			if node.n == nothing
				node.n = D()
			end
			node.n.field += f_j * (V_j * Wij)
		end
    end
end

#recursively walk the tree
function map_particle_to_AMRgrid_recursive!(node::Node{N,T,D}, field::Vector{T}, volume::Vector{T}, X::Vector{SVector{N,T}}, hsml::Vector{T},
	boxsizes::SVector{N,T}, treeNgb::Union{Node{N,T,D}, Nothing} ) where {N,T,D}

	if !isLeaf(node)
		@inbounds for i in 1:2^N
	        map_particle_to_AMRgrid_recursive!(node.child[i], field, volume, X, hsml, boxsizes, treeNgb)
		end
	else
	#if isLeaf(node)
		if hsml==T[] #NGP
			if node.n == nothing
				node.n = D()
			end
			if node.p !== nothing
				node.n.field = field[node.p.idx]
			else
				node.n.field = zero(T)
			end
			return
		end

		if treeNgb == nothing
			#already have a ngb list so we don't have to search for ngbs again
			node.n.field = zero(T)
		else
			if node.n == nothing
				node.n = D()
			end
        	node.n.idx_ngbs = get_scatter_ngb_tree(node.center, treeNgb, boxsizes)
		end

		if volume == T[]
			node.n.field = kernel_average_MFM(field, node.center, node.n.idx_ngbs, X, hsml, boxsizes)
		else
			node.n.field = kernel_average_SPH(field, node.center, node.n.idx_ngbs, X, hsml, volume, boxsizes)
		end
    end
end

function kernel_average_SPH(field::Vector{T}, X0::SVector{N,T}, idx_ngbs::Vector{Int64}, X::Vector{SVector{N,T}},
							hsml::Vector{T}, volume::Vector{T}, boxsizes::SVector{N,T}) where {N,T}
	res = zero(T)
	for k in eachindex(idx_ngbs)
		j = idx_ngbs[k]
		dx = nearest.(X[j] - X0, boxsizes)
		dr = norm(dx)
		Wij = kernel_cubic(dr/hsml[j]) / hsml[j]^N
		res += field[j] * (volume[j] * Wij)
		#res += 1.0  #debug
	end
	return res
end

function kernel_average_MFM(field::Vector{T}, X0::SVector{N,T}, idx_ngbs::Vector{Int64}, X::Vector{SVector{N,T}},
							hsml::Vector{T}, boxsizes::SVector{N,T}) where {N,T}
	res = zero(T)
	sigma = zero(T)
	for k in eachindex(idx_ngbs)
		j = idx_ngbs[k]
		dx = nearest.(X[j] - X0, boxsizes)
		dr = norm(dx)
		Wij = kernel_cubic(dr/hsml[j]) / hsml[j]^N
		res += field[j] * Wij
		sigma += Wij
	end
	return length(idx_ngbs) > 0 ? (res / sigma) : zero(T)
end

########## return the AMR field stored in the tree as an array (following the order of the tree)
function get_AMRfield(tree::Node{N,T,D}; max_depth::Int=MAX_DEPTH) where {N,T,D}
	fieldAMR = T[]
	root_node_length = tree.length[1]
	get_AMRfield_recursive!(fieldAMR, tree, max_depth, root_node_length)
	return fieldAMR
end

function get_AMRfield_recursive!(fieldAMR::Vector{T}, node::Node{N,T,D}, max_depth::Int, root_node_length::T) where {N,T,D, Tint<:Integer}
    if isLeaf(node)
		if node.n == nothing
			node.n = D()
		end
		push!(fieldAMR, node.n.field)
    else
		depth = log2(round(root_node_length / node.length[1]))
		if depth < max_depth
	    	@inbounds for i in 1:2^N
				get_AMRfield_recursive!(fieldAMR, node.child[i], max_depth, root_node_length)
	    	end
		else
			push!(fieldAMR, node.n.field)
		end
    end
end

########## project the 3D AMR field to a 2D image
#function project_AMRgrid_to_image(nx, ny, dimx, dimy, tree::Node{N,T,D}, rootcenter::SVector{N,T}, boxsizes::SVector{N,T}) where {N,T,D}
#	image = zeros(T, nx, ny)
#	project_AMRgrid_to_image_recursive!(image, nx, ny, dimx, dimy, tree, rootcenter, boxsizes);
#	return image
#end

function project_AMRgrid_to_image(nx, ny, dimx, dimy, tree::Node{N,T,D}, rootcenter::SVector{N,T}, boxsizes::SVector{N,T}; serial::Bool=false) where {N,T,D}
	image = zeros(T, nx, ny)
	if serial
		project_AMRgrid_to_image_recursive!(image, nx, ny, dimx, dimy, tree, rootcenter, boxsizes);
		return image
	end

	image_thread = [zeros(T, nx, ny) for i in 1:nthreads()]; #each thread has its own image
	if isLeaf(tree)
		project_AMRgrid_to_image_recursive!(image_thread[threadid()], nx, ny, dimx, dimy, tree, rootcenter, boxsizes)
	else
		@sync for i in 1:2^N
			if isLeaf(tree.child[i])
				project_AMRgrid_to_image_recursive!(image_thread[threadid()], nx, ny, dimx, dimy, tree.child[i], rootcenter, boxsizes)
			else
				for j in 1:2^N
					if isLeaf(tree.child[i].child[j])
						project_AMRgrid_to_image_recursive!(image_thread[threadid()], nx, ny, dimx, dimy, tree.child[i].child[j], rootcenter, boxsizes)
					else
						for k in 1:2^N
							Threads.@spawn project_AMRgrid_to_image_recursive!(image_thread[threadid()], nx, ny, dimx, dimy, tree.child[i].child[j].child[k], rootcenter, boxsizes)
						end
					end
				end
			end
		end
	end
	image .= sum(image_thread)
	return image
end

function project_AMRgrid_to_image_recursive!(image, nx, ny, dimx, dimy, node::Node{N,T,D}, rootcenter::SVector{N,T}, boxsizes::SVector{N,T}) where {N,T,D}
	@assert ispow2(nx) && ispow2(ny)

    if !isLeaf(node)
		@inbounds for i in 1:2^N
			project_AMRgrid_to_image_recursive!(image, nx, ny, dimx, dimy, node.child[i], rootcenter, boxsizes)
		end
	else
		#the lower left corner of the tree
		xmin = rootcenter[dimx] - 0.5 * boxsizes[dimx]
		ymin = rootcenter[dimy] - 0.5 * boxsizes[dimy]
		#tree cell boundaries
		x1 = node.center[dimx] - xmin - 0.5 * node.length[dimx]
		x2 = node.center[dimx] - xmin + 0.5 * node.length[dimx]
		y1 = node.center[dimy] - ymin - 0.5 * node.length[dimy]
		y2 = node.center[dimy] - ymin + 0.5 * node.length[dimy]
		#fraction of area when pixel is larger than node
		f_area = (node.length[dimx]*node.length[dimy]) / (boxsizes[dimx]*boxsizes[dimy]/(nx*ny))
		fac = min(f_area, 1.0)
		if N==3
			dimlos = filter!(x->(x!=dimx)&&(x!=dimy), [1,2,3])[1]
			fac *= node.length[dimlos]
		end
		dx = boxsizes[dimx] / nx
		dy = boxsizes[dimy] / ny
		dtol = 10
		#get rid of the round-off error with round()
		ix1r = round(x1 / dx, digits=dtol)
		ix2r = round(x2 / dx, digits=dtol)
		iy1r = round(y1 / dy, digits=dtol)
		iy2r = round(y2 / dy, digits=dtol)
		#convert into integer index
		ix1 = Int64(floor(ix1r) + 1)  #+1 becuase julia is one-based
		ix2 = Int64(floor(ix2r) + 1)
		iy1 = Int64(floor(iy1r) + 1)
		iy2 = Int64(floor(iy2r) + 1)
		tol = 10.0^(-dtol)
		@assert ix2 >= ix2r && iy2 >= iy2r
		if isinteger(ix2r) ix2 -= 1 end
		if isinteger(iy2r) iy2 -= 1 end
		if ix1>ix2 ix2=ix1 end #happens when node length < pixel size
		if iy1>iy2 iy2=iy1 end
		@assert ix2>=ix1 && iy2>=iy1

		if ix2 > nx || ix1 < 1
			@show ix1, ix2, ix1r, ix2r, x1, x2, dx, boxsizes[dimx], depth
		end
		if iy2 > ny || iy1 < 1
			@show iy1, iy2, y1, y2, dy, boxsizes[dimy], depth
		end

		#image[ix1:ix2,iy1:iy2] .+= node.n.field * fac
		#avoid using sub-array to reduce allocation
		for j in iy1:iy2, i in ix1:ix2
			image[i,j] += node.n.field * fac
		end
    end
end


text:
import Base: showerror

const AF_SUCCESS = 0
const AF_ERR_NO_MEM = 1
const AF_ERR_DRIVER = 2
const AF_ERR_RUNTIME = 3
const AF_ERR_INVALID_ARRAY = 4
const AF_ERR_ARG  = 5
const AF_ERR_SIZE = 6
const AF_ERR_TYPE = 7
const AF_ERR_DIFF_TYPE = 8
const AF_ERR_BATCH = 9
const AF_ERR_DEVICE = 10
const AF_ERR_NOT_SUPPORTED = 11
const AF_ERR_NOT_CONFIGURED  = 12
const AF_ERR_NONFREE = 13
const AF_ERR_NO_DBL = 14
const AF_ERR_NO_GFX = 15
const AF_ERR_LOAD_LIB = 16
const AF_ERR_LOAD_SYM  = 17
const AF_ERR_ARR_BKND_MISMATCH = 18
const AF_ERR_INTERNAL = 19
const AF_ERR_UNKNOWN = 20

export  AF_SUCCESS,
        AF_ERR_NO_MEM,
        AF_ERR_DRIVER,
        AF_ERR_RUNTIME,
        AF_ERR_INVALID_ARRAY,
        AF_ERR_ARG ,
        AF_ERR_SIZE,
        AF_ERR_TYPE,
        AF_ERR_DIFF_TYPE,
        AF_ERR_BATCH,
        AF_ERR_DEVICE ,
        AF_ERR_NOT_SUPPORTED ,
        AF_ERR_NOT_CONFIGURED  ,
        AF_ERR_NONFREE ,
        AF_ERR_NO_DBL ,
        AF_ERR_NO_GFX ,
        AF_ERR_LOAD_LIB ,
        AF_ERR_LOAD_SYM  ,
        AF_ERR_ARR_BKND_MISMATCH ,
        AF_ERR_INTERNAL ,
        AF_ERR_UNKNOWN

function throwAFerror(err::Integer)
    str = ccall((:af_err_to_string, af_lib),
                Cstring, (Cint, ), err)
    throw(ErrorException("ArrayFire Error ($err) : $(unsafe_string(str))"))
end


text:
root = RectanglePacker(Rect2(0,0,1024,1024))
push!(root, Rect2(0,0,20,20))
push!(root, Rect2(0,0,20,20))
push!(root, [Rect2(0,0, rand(5:50), rand(5:50)) for i=1:20])

# function get_rectangles(packer::Nothing, rectangles=Rect2i[])
#     return rectangles
# end
#
# function get_rectangles(packer, rectangles)
#     push!(rectangles, packer.area)
#     get_rectangles(packer.children.left, rectangles)
#     get_rectangles(packer.children.right, rectangles)
#     return rectangles
# end
# function get_rectangles(packer)
#     rectangles = Rect2i[]
#     get_rectangles(packer.children.left, rectangles)
#     get_rectangles(packer.children.right, rectangles)
#     return rectangles
# end
#
# rectangles = get_rectangles(root)
# using Makie
# linesegments(Rectf(0, 0, 1024, 1024))
# poly!(rectangles, color = (:red, 0.1), strokewidth=1, strokecolor=:black)
#


text:
@testset "1551.minimum-operations-to-make-array-equal.jl" begin
    @test min_operations_1551(3) == 2
    @test min_operations_1551(6) == 9
    @test min_operations_1551(58) == 841
    @test min_operations_1551(79) == 1560
end


text:
@testset "Supercell angles" begin
@testset "fcc" begin
    @test all(MkCell.cell_angles(fcc1*cellfcc)   .≈ (90.0, 120.0, 60.0))
    @test all(MkCell.cell_angles(fcc4*cellfcc)   .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(fcc7*cellfcc)   .≈ (99.59406822686046, 80.40593177313954, 80.40593177313954))
    @test all(MkCell.cell_angles(fcc13*cellfcc)  .≈ (72.54239687627792, 107.4576031237221, 72.54239687627792))
    @test all(MkCell.cell_angles(fcc23*cellfcc)  .≈ (90.0, 85.90395624184767, 107.9752838098089))
    @test all(MkCell.cell_angles(fcc26*cellfcc)  .≈ (85.90395624184767, 94.09604375815232, 94.09604375815232))
    @test all(MkCell.cell_angles(fcc32*cellfcc)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(fcc38*cellfcc)  .≈ (93.1847385367204, 86.8152614632796, 86.8152614632796))
    @test all(MkCell.cell_angles(fcc50*cellfcc)  .≈ (97.83747976308429, 82.16252023691571, 82.16252023691571))
    @test all(MkCell.cell_angles(fcc70*cellfcc)  .≈ (90.0, 87.79577249602796, 79.68680947509941))
    @test all(MkCell.cell_angles(fcc92*cellfcc)  .≈ (88.30279727450392, 88.08978682829007, 95.09758601873176))
    @test all(MkCell.cell_angles(fcc96*cellfcc)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(fcc100*cellfcc) .≈ (93.2771683126184, 86.7228316873816, 79.83575137827653))
    @test all(MkCell.cell_angles(fcc108*cellfcc) .≈ (90.0, 90.0, 90.0))
end;
@testset "cubic" begin
    @test all(MkCell.cell_angles(cubic1*cellcubic)   .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(cubic3*cellcubic)   .≈ (90.0, 90.0, 60.00000000000001))
    @test all(MkCell.cell_angles(cubic7*cellcubic)   .≈ (104.96321743330712, 75.03678256669288, 66.42182152179818))
    @test all(MkCell.cell_angles(cubic8*cellcubic)   .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(cubic13*cellcubic)  .≈ (90.0, 79.48026510934139, 99.59406822686046))
    @test all(MkCell.cell_angles(cubic14*cellcubic)  .≈ (99.59406822686046, 99.59406822686046, 99.59406822686046))
    @test all(MkCell.cell_angles(cubic24*cellcubic)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(cubic27*cellcubic)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(cubic33*cellcubic)  .≈ (90.0, 95.73917047726677, 90.0))
    @test all(MkCell.cell_angles(cubic38*cellcubic)  .≈ (85.20308428929593, 95.21590857045413, 75.47084566659821))
    @test all(MkCell.cell_angles(cubic52*cellcubic)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(cubic56*cellcubic)  .≈ (97.44883282786358, 98.21321070173819, 97.44883282786358))
    @test all(MkCell.cell_angles(cubic64*cellcubic)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(cubic73*cellcubic)  .≈ (90.0, 80.12503860217872, 93.1847385367204))
    @test all(MkCell.cell_angles(cubic76*cellcubic)  .≈ (93.1847385367204, 93.1847385367204, 93.1847385367204))
    @test all(MkCell.cell_angles(cubic93*cellcubic)  .≈ (92.666603476751, 92.86957619589583, 95.61383424102878))
    @test all(MkCell.cell_angles(cubic100*cellcubic) .≈ (90.0, 90.0, 90.0))
end;
@testset "bcc" begin
    @test all(MkCell.cell_angles(bcc1*cellbcc)   .≈ (70.52877936550931, 109.47122063449069, 70.52877936550931))
    @test all(MkCell.cell_angles(bcc2*cellbcc)   .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(bcc6*cellbcc)   .≈ (90.0, 90.0, 115.2394018206789))
    @test all(MkCell.cell_angles(bcc9*cellbcc)   .≈ (84.78409142954588, 95.21590857045413, 84.78409142954588))
    @test all(MkCell.cell_angles(bcc16*cellbcc)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(bcc21*cellbcc)  .≈ (86.98303869018402, 95.88878512679155, 95.88878512679155))
    @test all(MkCell.cell_angles(bcc25*cellbcc)  .≈ (90.0, 85.06265611930155, 104.96321743330712))
    @test all(MkCell.cell_angles(bcc28*cellbcc)  .≈ (80.40593177313954, 99.59406822686046, 80.40593177313954))
    @test all(MkCell.cell_angles(bcc35*cellbcc)  .≈ (87.87744863785665, 92.12255136214335, 87.87744863785665))
    @test all(MkCell.cell_angles(bcc48*cellbcc)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(bcc53*cellbcc)  .≈ (96.05074601591176, 83.86309610848119, 73.63722782068534))
    @test all(MkCell.cell_angles(bcc54*cellbcc)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(bcc65*cellbcc)  .≈ (77.91853011529311, 102.08146988470688, 77.91853011529311))
    @test all(MkCell.cell_angles(bcc84*cellbcc)  .≈ (92.24963404707006, 87.75036595292995, 84.78409142954588))
    @test all(MkCell.cell_angles(bcc91*cellbcc)  .≈ (88.87648134229958, 91.12351865770042, 88.87648134229958))
    @test all(MkCell.cell_angles(bcc103*cellbcc) .≈ (92.06928050584378, 93.99037534653964, 94.25094462500175))
    @test all(MkCell.cell_angles(bcc107*cellbcc) .≈ (90.0, 94.09604375815232, 86.00962465346038))
    @test all(MkCell.cell_angles(bcc112*cellbcc) .≈ (94.86140539579802, 85.13859460420197, 85.13859460420197))
    @test all(MkCell.cell_angles(bcc128*cellbcc) .≈ (90.0, 90.0, 90.0))
end;
@testset "hcp" begin
    @test all(MkCell.cell_angles(hcp1*cellhcp)   .≈ (90.0, 90.0, 60.0))
    @test all(MkCell.cell_angles(hcp3*cellhcp)   .≈ (90.0, 90.0, 60.0))
    @test all(MkCell.cell_angles(hcp6*cellhcp)   .≈ (90.0, 90.0, 97.7493663782984))
    @test all(MkCell.cell_angles(hcp8*cellhcp)   .≈ (86.62771331656612, 90.0, 90.0))
    @test all(MkCell.cell_angles(hcp11*cellhcp)  .≈ (93.10882512596606, 81.58250966338895, 76.22172246000032))
    @test all(MkCell.cell_angles(hcp14*cellhcp)  .≈ (95.7391704772668, 80.43938360731833, 78.01745450750637))
    @test all(MkCell.cell_angles(hcp16*cellhcp)  .≈ (90.0, 90.0, 98.21321070173819))
    @test all(MkCell.cell_angles(hcp21*cellhcp)  .≈ (85.05457028211667, 94.94542971788334, 85.05457028211667))
    @test all(MkCell.cell_angles(hcp27*cellhcp)  .≈ (99.04109614857515, 80.95890385142484, 80.95890385142484))
    @test all(MkCell.cell_angles(hcp37*cellhcp)  .≈ (87.38866773903177, 96.46147790595835, 87.99219492277635))
    @test all(MkCell.cell_angles(hcp39*cellhcp)  .≈ (93.26480184633287, 80.51388629606356, 88.73996568504259))
    @test all(MkCell.cell_angles(hcp42*cellhcp)  .≈ (91.26003431495744, 91.26003431495741, 82.05063659026476))
    @test all(MkCell.cell_angles(hcp53*cellhcp)  .≈ (97.09038962820946, 88.0060223885153, 85.52821741779152))
    @test all(MkCell.cell_angles(hcp62*cellhcp)  .≈ (93.27084136698927, 97.84269934104985, 95.41194366545838))
    @test all(MkCell.cell_angles(hcp65*cellhcp)  .≈ (81.55960115382774, 88.0060223885153, 93.63706568080396))
    @test all(MkCell.cell_angles(hcp69*cellhcp)  .≈ (90.0, 90.0, 94.30661909550145))
    @test all(MkCell.cell_angles(hcp76*cellhcp)  .≈ (89.57829310179544, 91.72387995423274, 94.12652542280888))
    @test all(MkCell.cell_angles(hcp80*cellhcp)  .≈ (86.62771331656612, 90.0, 90.0))
    @test all(MkCell.cell_angles(hcp88*cellhcp)  .≈ (90.80700949307113, 90.74639151793876, 90.74639151793875))
    @test all(MkCell.cell_angles(hcp90*cellhcp)  .≈ (90.0, 90.0, 90.0))
    @test all(MkCell.cell_angles(hcp100*cellhcp) .≈ (87.97513170272266, 90.0, 90.0))
end;

end;


text:
import Base.convert

convert(::Type{Int},    v::SubString{String}) = parse(Int, v)

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
function cost_matrices_with_dummy(truck_cost_mtx, drone_cost_mtx)
    Ct = [
        truck_cost_mtx          truck_cost_mtx[:, 1];
        truck_cost_mtx[1, :]'    0.0
    ]

    Cd = [
        drone_cost_mtx          drone_cost_mtx[:, 1];
        drone_cost_mtx[1, :]'    0.0
    ]

    return Ct, Cd
end

function _cost_matrices_with_dummy(x, y, speed_of_truck, speed_of_drone)
    n_nodes = length(x)
    @assert length(x) == length(y)

    dist = zeros(Float64, n_nodes, n_nodes)
    for i in 1:n_nodes
        for j in 1:n_nodes
            dist[i, j] = (x[i]-x[j])^2 + (y[i]-y[j])^2 |> sqrt
        end
    end

    Ct = speed_of_truck .* dist
    Cd = speed_of_drone .* dist

    return Ct, Cd
end

function cost_matrices_with_dummy(x, y, speed_of_truck, speed_of_drone)
    xx = copy(x)
    yy = copy(y)
    push!(xx, x[1])
    push!(yy, y[1])
    return _cost_matrices_with_dummy(xx, yy, speed_of_truck, speed_of_drone)
end

function travel_cost(path::Vector{Int}, C::Matrix{T}) where T
    # @show path
    sum = zero(T)
    for i in 1:length(path)-1
        sum += C[path[i], path[i+1]]
    end
    return sum
end

function objective_value(truck_route, drone_route, Ct, Cd)
    combined_nodes = intersect(truck_route, drone_route)
    obj_val = 0.0
    for i in 1:length(combined_nodes)-1
        j1 = combined_nodes[i]
        j2 = combined_nodes[i+1]

        t_idx1 = findfirst(x -> x == j1, truck_route)
        t_idx2 = findfirst(x -> x == j2, truck_route)
        t_cost = travel_cost(truck_route[t_idx1:t_idx2], Ct)

        d_idx1 = findfirst(x -> x == j1, drone_route)
        d_idx2 = findfirst(x -> x == j2, drone_route)
        d_cost = travel_cost(drone_route[d_idx1:d_idx2], Cd)

        obj_val += max(t_cost, d_cost)
    end
    return obj_val
end


text:
module Atum

using Reexport
using StaticArrays
@reexport using Bennu

using KernelAbstractions
using KernelAbstractions.Extras: @unroll
using CUDAKernels
using UnPack

include("balancelaw.jl")
include("numericalfluxes.jl")
include("dgsem.jl")
include("esdgsem.jl")
include("odesolvers.jl")

include("balancelaws/advection.jl")
include("balancelaws/shallow_water.jl")
include("balancelaws/euler.jl")
include("balancelaws/euler_gravity.jl")

end # module


text:
# We have to use the original rules as they may be modified
# elsewhere, like in an Interact.jl interface. So here we build them
# into the SimData object after running _modifyrules on them
function _updaterules(rules::Tuple, sd::SimData)
    newrs = ModelParameters.setparent(
        ruleset(sd),
        _modifyrules(_proc_setup(proc(sd), ModelParameters.stripparams(rules)), sd)
    )
    @set sd.ruleset = newrs
end

_modifyrules(rules::Tuple, simdata) =
    (modifyrule(rules[1], simdata), _modifyrules(tail(rules), simdata)...)
_modifyrules(rules::Tuple{}, simdata) = ()

# The default is to return a rule unchanged
modifyrule(rule, simdata) = rule


text:
export DCTOp

mutable struct DCTOp{T} <: AbstractLinearOperator{T}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod :: Function
  tprod :: Nothing
  ctprod :: Function
  nprod :: Int
  ntprod :: Int
  nctprod :: Int
  plan
  iplan
  dcttype::Int64
end

"""
  DCTOp(T::Type, shape::Tuple, dcttype=2)

returns a `DCTOp <: AbstractLinearOperator` which performs a DCT on a given input array.

# Arguments:
* `T::Type`       - type of the array to transform
* `shape::Tuple`  - size of the array to transform
* `dcttype`       - type of DCT (currently `2` and `4` are supported)
"""
function DCTOp(T::Type, shape::Tuple, dcttype=2)
  if dcttype == 2
    plan = plan_dct(zeros(T,shape))
    iplan = plan_idct(zeros(T,shape))
    return DCTOp{T}(prod(shape), prod(shape), true, false
            , x->vec(plan*reshape(x,shape))
            , nothing
            , y->vec(iplan*reshape(y,shape))
            , 0, 0, 0
            , plan
            , iplan
            , dcttype)
  elseif dcttype == 4
    factor = sqrt(1.0/(prod(shape)* 2^length(shape)) )
    plan = FFTW.plan_r2r(zeros(T,shape),FFTW.REDFT11)
    # return DCTOp{T}(prod(shape), prod(shape), true, false
    #         , x->vec((plan*reshape(x,shape)).*factor)
    #         , x->vec((plan*reshape(x,shape)).*factor)
    #         , nothing
    #         , 0, 0, 0
    #         , plan
    #         , plan
    #         , dcttype)
    return LinearOperator{T}(prod(shape), prod(shape), true, false
            , x->T.( vec((plan*reshape(x,shape)).*factor) )
            , x->T.( vec((plan*reshape(x,shape)).*factor) )
            , nothing )
  else
    error("DCT type $(dcttype) not supported")
  end
end

function Base.copy(S::DCTOp)
  return DCTOp(eltype(S), size(S.plan), S.dcttype)
end


text:
# This file is a part of InverseFunctions.jl, licensed under the MIT License (MIT).

"""
    inverse(f)

Returns the inverse of a function `f`.

The following conditions must be satisfied:

* `inverse(f) ∘ f` must be equivalent to `identity`.
* `inverse(f)(f(x)) ≈ x`
* `inverse(inverse(f))` must be equivalent (ideally identical) to `f`.

`inverse` supports mapped/broadcasted functions (via `Base.Fix1`) and (on
Julia >=v1.6) function composition.

Example:

```julia
foo(x) = inv(exp(-x) + 1)
inv_foo(y) = log(y / (1 - y))

InverseFunctions.inverse(::typeof(foo)) = inv_foo
InverseFunctions.inverse(::typeof(inv_foo)) = foo

x = 4.2
@assert inverse(foo)(foo(x)) ≈ x
@assert inverse(inverse(foo)) == foo

X = rand(10)
broadcasted_foo = Base.Fix1(broadcast, foo)
Y = broadcasted_foo(X)
@assert inverse(broadcasted_foo)(Y) ≈ X

# Requires Julia >= v1.6:
bar = log ∘ foo
@assert inverse(bar)(bar(x)) ≈ x
```
"""
function inverse end
export inverse

inverse(::typeof(inverse)) = inverse

@static if VERSION >= v"1.6"
    inverse(f::Base.ComposedFunction) = Base.ComposedFunction(inverse(f.inner), inverse(f.outer))
end

inverse(mapped_f::Base.Fix1{<:Union{typeof(map),typeof(broadcast)}}) = Base.Fix1(mapped_f.f, inverse(mapped_f.x))

inverse(::typeof(identity)) = identity
inverse(::typeof(inv)) = inv
inverse(::typeof(adjoint)) = adjoint
inverse(::typeof(transpose)) = transpose

inverse(::typeof(exp)) = log
inverse(::typeof(log)) = exp

inverse(::typeof(exp2)) = log2
inverse(::typeof(log2)) = exp2

inverse(::typeof(exp10)) = log10
inverse(::typeof(log10)) = exp10

inverse(::typeof(expm1)) = log1p
inverse(::typeof(log1p)) = expm1


text:
export SeisSrc

"""
SeisSrc: container for descriptions of a seismic source process

    S = SeisSrc()

Initialize an empty SeisSrc object. Fields can be initialized at creation with keywords; for example, S = SeisSrc(m0 = 1.6e22).

| Field   | Type                | Meaning                                 |
|:----    |:-----               |:--------                                |
| id      | String              | source process ID                       |
| eid     | String              | event ID (note, generally :id != :eid)  |
| m0      | Float64             | scalar seismic moment                   |
| mt      | Array{Float64,1}    | seismic moment tensor                   |
| dm      | Array{Float64,1}    | seismic moment tensor misfit            |
| npol    | Int64               | number of polarities in focal mechanism |
| gap     | Float64             | max azimuthal gap in focal mechanism    |
| pax     | Array{Float64,2}    | principal axes                          |
| planes  | Array{Float64,2}    | nodal planes                            |
| src     | String              | data source string (filename or URL)    |
| st      | SourceTime          | source-time subfield                    |
| misc    | Dict{String,Any}    | dictionary of non-essential information |
| notes   | Array{String,1}     | notes and automated logging             |
See also: EQLoc, EQMag
"""
mutable struct SeisSrc
  id      ::String
  eid     ::String
  m0      ::Float64
  mt      ::Array{Float64,1}
  dm      ::Array{Float64,1}
  npol    ::Int64
  gap     ::Float64
  pax     ::Array{Float64,2}
  planes  ::Array{Float64,2}
  src     ::String
  st      ::SourceTime
  misc    ::Dict{String,Any}
  notes   ::Array{String,1}

  function SeisSrc(
                    id      ::String,
                    eid     ::String,
                    m0      ::Float64,
                    mt      ::Array{Float64,1},
                    dm      ::Array{Float64,1},
                    npol    ::Int64,
                    gap     ::Float64,
                    pax     ::Array{Float64,2},
                    planes  ::Array{Float64,2},
                    src     ::String,
                    st      ::SourceTime,
                    misc    ::Dict{String,Any},
                    notes   ::Array{String,1},
                    )
    return new(id, eid, m0, mt, dm, npol, gap, pax, planes, src, st, misc, notes)
  end
end
SeisSrc(;
          id     ::String             = "",
          eid     ::String            = "",
          m0      ::Float64           = zero(Float64),
          mt      ::Array{Float64,1}  = Float64[],
          dm      ::Array{Float64,1}  = Float64[],
          npol    ::Int64             = zero(Int64),
          gap     ::Float64           = zero(Float64),
          pax     ::Array{Float64,2}  = Array{Float64, 2}(undef, 0, 0),
          planes  ::Array{Float64,2}  = Array{Float64, 2}(undef, 0, 0),
          src     ::String            = "",
          st      ::SourceTime        = SourceTime(),
          misc    ::Dict{String,Any}  = Dict{String,Any}(),
          notes   ::Array{String,1}   = String[],
          ) = SeisSrc(id, eid, m0, mt, dm, npol, gap, pax, planes, src, st, misc, notes)

function isempty(S::SeisSrc)
  q::Bool = min((getfield(S, :m0) == zero(Float64)),
                (getfield(S, :gap) == zero(Float64)),
                (getfield(S, :npol) == zero(Int64)))
  for f in (:id, :eid, :misc, :mt, :dm, :notes, :pax, :planes, :src, :st)
    q = min(q, isempty(getfield(S, f)))
  end
  return q
end

function isequal(S::SeisSrc, U::SeisSrc)
  q::Bool = true
  for f in fieldnames(SeisSrc)
    if f != :notes
      q = min(q, getfield(S,f) == getfield(U,f))
    end
  end
  return q
end
==(S::SeisSrc, U::SeisSrc) = isequal(S, U)

sizeof(S::SeisSrc) = 136 +
  sum([sizeof(getfield(S, f)) for f in (:id, :eid, :mt, :dm, :pax, :planes, :src, :st, :misc, :notes)])

function show(io::IO, S::SeisSrc)
  if get(io, :compact, false) == false
    for f in fieldnames(SeisSrc)
      fn = lpad(uppercase(String(f)), 6, " ")
      if f == :misc
        println(io, fn, ": ", length(getfield(S, :misc)), " items")
      elseif f == :notes
        println(io, fn, ": ", length(getfield(S, :notes)), " entries")
      elseif f == :id || f == :src
        println(io, fn, ": ", getfield(S, f))
      else
        println(io, fn, ": ", repr(getfield(S, f), context=:compact=>true))
      end
    end
  else
    c = :compact => true

    # Order of preference: mt, np, pax
    mech_str = string("m₀ = ", repr(getfield(S, :m0), context=c), "; S = ",
                      repr(getfield(S, :mt), context=c), "; ",
                      "NP = ", repr(getfield(S, :planes), context=c), "; ",
                      "PAX = ", repr(getfield(S, :pax), context=c))

    L = length(mech_str)
    L_max = displaysize(io)[2]
    if L > L_max
      print(io, mech_str[1:L_max-1], "…")
    else
      print(io, mech_str)
    end
  end
  return nothing
end
show(S::SeisSrc) = show(stdout, S)

# SeisSrc
function write(io::IO, S::SeisSrc)
  id  = codeunits(getfield(S, :id))
  eid = codeunits(getfield(S, :eid))
  src = codeunits(getfield(S, :src))

  # Write begins ------------------------------------------------------
  write(io, Int64(length(id)))
  write(io, id)                                             # id
  write(io, Int64(length(eid)))
  write(io, eid)                                            # eid
  write(io, S.m0)                                           # m0
  write(io, Int64(length(S.mt)))
  write(io, S.mt)                                           # mt
  write(io, Int64(length(S.dm)))
  write(io, S.dm)                                           # dm
  write(io, S.npol)                                         # npol
  write(io, S.gap)                                          # gap
  r, c = size(S.pax)
  write(io, Int64(r), Int64(c)),
  write(io, S.pax)                                          # pax
  r, c = size(S.planes)
  write(io, Int64(r), Int64(c)),
  write(io, S.planes)                                       # planes
  write(io, Int64(length(src)))
  write(io, src)                                            # src
  write(io, getfield(S, :st))                               # st
  write_misc(io, getfield(S, :misc))                        # misc
  write_string_vec(io, getfield(S, :notes))                 # notes

  # Write ends --------------------------------------------------------
  return nothing
end

function read(io::IO, ::Type{SeisSrc})
  u = getfield(BUF, :buf)
  return SeisSrc( String(read(io, read(io, Int64))),                    # :id
                  String(read(io, read(io, Int64))),                    # :eid
                  read(io, Float64),                                    # :m0
                  read!(io, Array{Float64, 1}(undef, read(io, Int64))), # :mt
                  read!(io, Array{Float64, 1}(undef, read(io, Int64))), # :dm
                  read(io, Int64),                                      # :npol
                  read(io, Float64),                                    # :gap
                  read!(io, Array{Float64, 2}(undef,
                    read(io, Int64), read(io, Int64))),                 # :pax
                  read!(io, Array{Float64, 2}(undef,
                    read(io, Int64), read(io, Int64))),                 # :planes
                  String(read(io, read(io, Int64))),                    # :src
                  read(io, SourceTime),                                 # :st
                  read_misc(io, u),                                     # :misc
                  read_string_vec(io, u) )                              # :notes
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 60941eaa-1aea-11eb-1277-97b991548781
begin
    using PlutoVista
    using VoronoiFVM
	using PlutoUI
end

# ╔═╡ 5e13b3db-570c-4159-939a-7e2268f0a102
md"""
# Bernoulli function test

We test the implementation of the Bernoulli function in VoronoiFVM against the evaluation
with BigFloat. This allows to optimize thresholds for switching between evaluation expressions.
"""

# ╔═╡ 20ce5908-7fe7-4e8f-baf8-c858d2d12ce5
TableOfContents(title="")

# ╔═╡ 1774defe-1939-45f9-95e5-6a9fff70f3f2
md"""
## The candidates
"""

# ╔═╡ b47b781b-ec11-485a-9de6-061ad0957f46
"""
Benchmark implementation  using BigFloat
"""
function B_Big(x)
	bx=BigFloat(x)
	Float64(bx/(exp(bx)-one(bx)))
end

# ╔═╡ 01dd7870-45a7-42e5-ad26-7933bdec0c60
"""
Naive implementation using FLoat64
"""
B_Float64(x)=x/(exp(x)-1)

# ╔═╡ 36158a4c-4a4b-44c8-9226-448713304335
"""
Approximation for large positive and negative arguments
"""
B_large(x)= x>0 ? 0.0 : -x

# ╔═╡ 54de85ff-8d61-4d15-b7c6-939ce0ea10fc
"""
Taylor scheme approximation for small  positive and negative arguments
"""
B_taylor(x)=VoronoiFVM.bernoulli_horner(x)

# ╔═╡ 8e77f24b-b3fd-4069-82ba-799dc6aa602a
"""
Implementation  in VoronoiFVM
"""
B_vfvm(x)=VoronoiFVM.fbernoulli(x)

# ╔═╡ 984bd562-271b-4e7e-b289-95d95e92ec2b
md"""
This implementation uses `B_taylor` if  `|exp(x)-1|<` $(VoronoiFVM.bernoulli_small_threshold) and
`B_large` if `|x|>` $(VoronoiFVM.bernoulli_large_threshold).
"""

# ╔═╡ c3fd0ff2-7111-4165-ad93-d6d7257301fa
md"""
## Approximation for small x

For small values of x, a Taylor approximation implementd using a Horner scheme is utilized, as the exponential expression runs into errors in the vicinity of zero and fails to evaluate at zero.. As as long as its error is large than that of the Taylor approximation calculated with the Taylor scheme, we should use the later one.
"""

# ╔═╡ 56ff3f5c-6fe9-4d44-a5ae-449c42efca62
smallX=collect(-0.5:1.0e-4+1.0e-8:0.5);

# ╔═╡ 6e7c197b-8ad2-4b9d-a0bc-40a48db32387
let
	p=PlutoVistaPlot(resolution=(500,200),legend=:ct,xlabel="x",ylabel="error")
    plot!(p,smallX,abs.(B_Big.(smallX).-B_Float64.(smallX)),yscale=:log,label="|B_Big(x)-B_Float64(x)|")
    plot!(p,smallX,abs.(B_Big.(smallX).- B_taylor.(smallX)),label="|B_Big(x)-B_taylor(x)|")
end

# ╔═╡ 65aecf78-f88b-4399-abac-717c3c62a285
md"""
## Approximation for large x

Here, an important aspect is to prevent floating point overflow for large x.
"""

# ╔═╡ 26cdb920-291a-4b54-963f-fd9bd610662f
largeX=-100:1.00001e-3:100;

# ╔═╡ 8633969b-33cb-486a-8731-2ec7dcb881d7
let
	p=PlutoVistaPlot(resolution=(500,200),legend=:lt)
    plot!(p,largeX,abs.(B_Big.(largeX).-B_Float64.(largeX)) ,yscale=:log,label="|B_Big(x)-B_Float64(x)|")
    plot!(p,largeX,abs.(B_Big.(largeX).-B_large.(largeX)),label="|B_Big(x)-B_large(x)|")
end

# ╔═╡ 5a293797-beb9-493e-af12-d978c50d6148
md"""
## Test of the Implementation in VoronoiFVM
"""

# ╔═╡ ed8f172e-1d74-4514-b1e6-815ec9c87ae5
maxerror(X)=maximum(abs.(B_Big.(X).-B_vfvm.(X)));

# ╔═╡ feb21ce6-0ddc-45fb-90f4-1e46261a9110
plot(smallX,abs.(B_Big.(smallX).-B_vfvm.(smallX)),label="|B_Big(x)-B_vfvm(x)|",legend=:lt,yscale=:log,resolution=(500,200),color=:red,xlabel="x",ylabel="error", title="Max error for small x: $(maxerror(smallX))" )

# ╔═╡ 91d3b907-9053-4467-a8ab-be9c5597741a
plot(largeX,abs.(B_Big.(largeX).-B_vfvm.(largeX)),label="|B_Big(x)-B_vfvm(x)|",legend=:lb,yscale=:log,resolution=(500,200),color=:red,xlabel="x",ylabel="error", title="Max error for large x: $(maxerror(largeX))" )

# ╔═╡ 12f268d2-baa6-4c1b-bed5-e9df53b469fc
html"<hr>"

# ╔═╡ 80263ff2-a7c5-42c1-bd3e-ee7799622f76
md"""
begin
    using Pkg
    Pkg.activate(".testenv")
    Pkg.add("Revise")
    using Revise
    Pkg.develop("VoronoiFVM")
    Pkg.add(["PlutoVista", "PlutoUI"])
end
""";

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PlutoVista = "646e1f28-b900-46d7-9d87-d554eb38a413"
VoronoiFVM = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"

[compat]
PlutoUI = "~0.7.16"
PlutoVista = "~0.8.7"
VoronoiFVM = "~0.13.2"
"""

# ╔═╡ 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-rc2"
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 = "d9352737cef8525944bf9ef34392d756321cbd54"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.38"

[[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 = "3533f5a691e60601fe60c90d8bc47a27aa2907ec"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.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 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.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 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[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 = "1e8e50f054057f23e908fbd6935766dca6293cc2"
uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
version = "0.8.7"

[[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 = "2db648b6712831ecb333eae76dbfd1c156ca13bb"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.2"

[[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", "Preferences", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "63777916efbcb0ab6173d09a658fb7f2783de485"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.21"

[[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 = "925ba2f11df005d894b113292d32fca9afe3f8c8"
uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
version = "0.3.9"

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

[[deps.HypertextLiteral]]
git-tree-sha1 = "5efcf53d798efede8fee5b2c8b09284be359bf24"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.2"

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

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

[[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.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "f0c6489b12d28fb4c2103073ec7452f3423bd308"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.1"

[[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", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "6193c3815f13ba1b78a51ce391db8be016ae9214"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.4"

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

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

[[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 = "d911b6a12ba974dabe2291c6d450094a7226b372"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.1"

[[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", "DocStringExtensions", "GridVisualize", "UUIDs"]
git-tree-sha1 = "b99d4e38e7dba4535cee937e0444aed5912245d0"
uuid = "646e1f28-b900-46d7-9d87-d554eb38a413"
version = "0.8.7"

[[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 = ["SHA", "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 = "c944fa4adbb47be43376359811c0a14757bdc8a8"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.20.0"

[[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 = "999444c932d3bd32e7cb5a349e219bec9cfa2bd5"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "1.17.1"

[[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 = "f0bccf98e16759818ffc5d97ac3ebf87eb950150"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.8.1"

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

[[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 = "65f1d72aa575d3d6c348fcfc0d0d905fdbd7f5c3"
uuid = "82b139dc-5afc-11e9-35da-9b9bdfd336f3"
version = "0.13.2"

[[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:
# ╠═60941eaa-1aea-11eb-1277-97b991548781
# ╟─5e13b3db-570c-4159-939a-7e2268f0a102
# ╟─20ce5908-7fe7-4e8f-baf8-c858d2d12ce5
# ╟─1774defe-1939-45f9-95e5-6a9fff70f3f2
# ╠═b47b781b-ec11-485a-9de6-061ad0957f46
# ╠═01dd7870-45a7-42e5-ad26-7933bdec0c60
# ╠═36158a4c-4a4b-44c8-9226-448713304335
# ╠═54de85ff-8d61-4d15-b7c6-939ce0ea10fc
# ╠═8e77f24b-b3fd-4069-82ba-799dc6aa602a
# ╟─984bd562-271b-4e7e-b289-95d95e92ec2b
# ╟─c3fd0ff2-7111-4165-ad93-d6d7257301fa
# ╠═56ff3f5c-6fe9-4d44-a5ae-449c42efca62
# ╠═6e7c197b-8ad2-4b9d-a0bc-40a48db32387
# ╠═65aecf78-f88b-4399-abac-717c3c62a285
# ╠═26cdb920-291a-4b54-963f-fd9bd610662f
# ╠═8633969b-33cb-486a-8731-2ec7dcb881d7
# ╟─5a293797-beb9-493e-af12-d978c50d6148
# ╠═ed8f172e-1d74-4514-b1e6-815ec9c87ae5
# ╟─feb21ce6-0ddc-45fb-90f4-1e46261a9110
# ╟─91d3b907-9053-4467-a8ab-be9c5597741a
# ╟─12f268d2-baa6-4c1b-bed5-e9df53b469fc
# ╠═80263ff2-a7c5-42c1-bd3e-ee7799622f76
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002