text:
module StaticLint
using SymbolServer
using CSTParser
using CSTParser: isidentifier
using CSTParser: Scope, Binding, EXPR, PUNCTUATION, IDENTIFIER, KEYWORD, OPERATOR
using CSTParser: Call, UnaryOpCall, BinaryOpCall, WhereOpCall, Import, Using, Export, TopLevel, ModuleH, BareModule, Quote, Quotenode, MacroName, MacroCall, Macro, x_Str, FileH, Parameters
using CSTParser: setparent!, setscope!
# to be removed after CSTParser change
kindof(x::EXPR) = x.kind
kindof(t::CSTParser.Tokens.AbstractToken) = t.kind
valof(x::EXPR) = x.val
typof(x::EXPR) = x.typ
parentof(x::EXPR) = x.parent
parentof(s::Scope) = s.parent
scopeof(x::EXPR) = x.scope
bindingof(x::EXPR) = x.binding
refof(x::EXPR) = x.ref

function setref!(x::EXPR, r)
    x.ref = r
    return x
end

const NoReference = Binding("0NoReference", nothing, nothing, [], nothing)

mutable struct State
    file
    scope::Scope
    delayed::Bool
    ignorewherescope::Bool
    quoted::Bool
    urefs::Vector{EXPR}
    server
end

function (state::State)(x)
    delayed = state.delayed # store states
    isquoted = state.quoted

    if quoted(x)
        state.quoted = true
    elseif state.quoted && unquoted(x)
        state.quoted = false
    end
    # imports
    if typof(x) === Using || typof(x) === Import
        resolve_import(x, state)
    elseif typof(x) === Export # Allow delayed resolution
        state.delayed = true
    end

    #bindings
    add_binding(x, state)
    mark_globals(x, state)

    #macros
    handle_macro(x, state)

    # scope
    clear_scope(x)
    s0 = state.scope
    if typof(x) === FileH
        setscope!(x, state.scope)
    elseif scopeof(x) isa Scope
        if CSTParser.defines_function(x) || CSTParser.defines_macro(x)
            state.delayed = true # Allow delayed resolution
        end
        scopeof(x) != s0 && setparent!(scopeof(x), s0)
        state.scope = scopeof(x)
        if typof(x) === ModuleH # Add default modules to a new module
            state.scope.modules = Dict{String,Any}()
            state.scope.modules["Base"] = getsymbolserver(state.server)["Base"]
            state.scope.modules["Core"] = getsymbolserver(state.server)["Core"]
        elseif typof(x) === BareModule
            state.scope.modules = Dict{String,Any}()
            state.scope.modules["Core"] = getsymbolserver(state.server)["Core"]
        end
        if (typof(x) === CSTParser.ModuleH || typof(x) === CSTParser.BareModule) && bindingof(x) !== nothing # Add reference to out of scope binding (i.e. itself)
            state.scope.names[bindingof(x).name] = bindingof(x)
        elseif typof(x) === CSTParser.Flatten && typof(x.args[1]) === CSTParser.Generator && x.args[1].args isa Vector{EXPR} && length(x.args[1].args) > 0 && typof(x.args[1].args[1]) === CSTParser.Generator
            setscope!(x.args[1].args[1], nothing)
        end
    end
    followinclude(x, state) # follow include
    if state.quoted
        if isidentifier(x) && !hasref(x)
            setref!(x, NoReference)
        end
    elseif typof(x) === CSTParser.Quotenode && length(x.args) == 2 && kindof(x.args[1]) === CSTParser.Tokens.COLON && typof(x.args[2]) === CSTParser.IDENTIFIER
        setref!(x.args[2], NoReference)
    elseif (isidentifier(x) && !hasref(x)) || resolvable_macroname(x) || typof(x) === x_Str || (typof(x) === BinaryOpCall && kindof(x.args[2]) === CSTParser.Tokens.DOT)
        resolved = resolve_ref(x, state.scope, state)
        if !resolved && (state.delayed || isglobal(valof(x), state.scope))
            push!(state.urefs, x)
        end
    end
    # traverse across children (evaluation order)
    if typof(x) === CSTParser.BinaryOpCall && (CSTParser.is_assignment(x) && !CSTParser.is_func_call(x.args[1]) || typof(x.args[2]) === CSTParser.Tokens.DECLARATION) && !(CSTParser.is_assignment(x) && typof(x.args[1]) === CSTParser.Curly)
        state(x.args[3])
        state(x.args[2])
        state(x.args[1])
    elseif typof(x) === CSTParser.WhereOpCall
        @inbounds for i = 3:length(x.args)
            state(x.args[i])
        end
        state(x.args[1])
        state(x.args[2])
    elseif typof(x) === CSTParser.Generator
        @inbounds for i = 2:length(x.args)
            state(x.args[i])
        end
        state(x.args[1])
    elseif typof(x) === CSTParser.Flatten && x.args !== nothing && length(x.args) === 1 && x.args[1].args !== nothing && length(x.args[1]) >= 3 && length(x.args[1].args[1]) >= 3
        for i = 3:length(x.args[1].args[1].args)
            state(x.args[1].args[1].args[i])
        end
        for i = 3:length(x.args[1].args)
            state(x.args[1].args[i])
        end
        state(x.args[1].args[1].args[1])
    elseif x.args !== nothing
        @inbounds for i in 1:length(x.args)
            state(x.args[i])
        end
    end
    checks(x, state.server)

    # return to previous states
    state.scope != s0 && (state.scope = s0)
    state.delayed = delayed
    state.quoted = isquoted
    return state.scope
end

function add_binding(x, state)
    scope = state.scope
    if bindingof(x) isa Binding
        # check for global marker
        if isglobal(bindingof(x).name, scope)
            scope = _get_global_scope(state.scope)
        end

        if typof(x) === Macro
            scope.names[string("@", bindingof(x).name)] = bindingof(x)
            mn = CSTParser.get_name(x)
            if typof(mn) === IDENTIFIER
                setref!(mn, bindingof(x))
            end
        else
            if haskey(scope.names, bindingof(x).name)
                bindingof(x).overwrites = scope.names[bindingof(x).name]
            end
            scope.names[bindingof(x).name] = bindingof(x)
        end
        infer_type(bindingof(x), scope, state)
    elseif bindingof(x) isa SymbolServer.SymStore
        scope.names[valof(x)] = bindingof(x)
    end
end

isglobal(name, scope) = haskey(scope.names, "#globals") && name in scope.names["#globals"].refs

function mark_globals(x, state)
    if typof(x) === CSTParser.Global
        if !haskey(state.scope.names, "#globals")
            state.scope.names["#globals"] = CSTParser.Binding("#globals", nothing, nothing, String[], nothing)
        end
        if x.args isa Vector{EXPR}
            for i = 2:length(x.args)
                if typof(x.args[i]) === CSTParser.IDENTIFIER && !haskey(state.scope.names, valof(x.args[i]))
                    push!(state.scope.names["#globals"].refs, valof(x.args[i]))
                end
            end
        end
    end
end

function followinclude(x, state::State)
    if typof(x) === Call && typof(x.args[1]) === IDENTIFIER && valof(x.args[1]) == "include"
        path = get_path(x)
        if isempty(path)
        elseif hasfile(state.server, path)
        elseif canloadfile(state.server, path)
            loadfile(state.server, path)
        elseif hasfile(state.server, joinpath(dirname(getpath(state.file)), path))
            path = joinpath(dirname(getpath(state.file)), path)
        elseif canloadfile(state.server, joinpath(dirname(getpath(state.file)), path))
            path = joinpath(dirname(getpath(state.file)), path)
            loadfile(state.server, path,)
        else
            path = ""
        end
        if !isempty(path)
            oldfile = state.file
            state.file = getfile(state.server, path)
            setroot(state.file, getroot(oldfile))
            setscope!(getcst(state.file), nothing)
            state(getcst(state.file))
            state.file = oldfile
        else
            # (printstyled(">>>>Can't follow include", color = :red);printstyled(" $(Expr(x)) from $(dirname(state.path))\n"))
            # error handling for broken `include` here
        end
    end
end

include("server.jl")
include("imports.jl")
include("references.jl")
include("macros.jl")
include("checks.jl")
include("type_inf.jl")
include("utils.jl")
end


text:
"""
    u(xy::AbstractVector{T}, λ::arb, k::Int; boundary = nothing, notransform::Bool = false)
    u(xy::AbstractVector{T}, λ::arb, ks::UnitRange{Int}; boundary = nothing, notransform::Bool = false)
    u(xy::AbstractVector{T}, λ::arb; boundary = nothing, notransform::Bool = false)

Evaluate the eigenfunction.

The point can be given in either Cartesian or Polar coordinates. By
default a coordinate transformation is applied to switch from the
standard coordinate system to that used by `u`. If `notransform` is
true then this coordinate transformation is not performed, this
assumes that they already given in the coordinate system used by `u`.

If `k::Int` is given then evaluate only that basis function and don't
multiply with any coefficient. If `ks::UnitRange{Int}` is given then
return a vector with the values of the corresponding basis functions.

If `boundary` is set to an integer then the point is assumed to come
from the corresponding boundary of the domain, some eigenfunctions are
then identically equal to zero.

See also: [`coordinate_transform`](@ref)
"""
function (u::AbstractPlanarEigenfunction) end

(u::AbstractPlanarEigenfunction)(
    xy::AbstractVector{T},
    λ::arb,
    k::Integer;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}} =
    u(xy, λ, k:k, boundary = boundary, notransform = notransform)

function (u::AbstractPlanarEigenfunction)(
    xy::AbstractVector{T},
    λ::arb;
    boundary = nothing,
    notransform::Bool = false,
) where {T<:Union{arb,arb_series}}
    # TODO: This computes all terms even if some terms become NaN very
    # early on. Could potentially be optimized.
    coeffs = coefficients(u)
    isempty(coeffs) && return zero(λ)
    return sum(coeffs .* u(xy, λ, 1:length(coeffs), boundary = boundary))
end

ext: jl


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges

const MOIU = MOI.Utilities
MOIU.@model(ECOSModelData,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.ExponentialCone),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction,),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))
const optimizer = MOIU.CachingOptimizer(ECOSModelData{Float64}(), ECOS.Optimizer(verbose=false))

# SOC2 requires 1e-4
const config = MOIT.TestConfig(atol=1e-4, rtol=1e-4)

@testset "Unit" begin
    MOIT.unittest(MOIB.SplitInterval{Float64}(optimizer), config,
                  [# Quadratic functions are not supported
                   "solve_qcp_edge_cases", "solve_qp_edge_cases",
                   # Integer and ZeroOne sets are not supported
                   "solve_integer_edge_cases", "solve_objbound_edge_cases"])
end

@testset "Continuous linear problems" begin
    MOIT.contlineartest(MOIB.SplitInterval{Float64}(optimizer), config)
end

@testset "Continuous conic problems" begin
    exclude = ["sdp", "rootdet", "logdet"]
    MOIT.contconictest(MOIB.GeoMean{Float64}(MOIB.RSOC{Float64}(optimizer)),
                       config, exclude)
end

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "ECOS"
end

@testset "supports_allocate_load" begin
    @test MOIU.supports_allocate_load(optimizer.optimizer, false)
    @test !MOIU.supports_allocate_load(optimizer.optimizer, true)
end


text:
module Boltzmann

export RBM,
       Units,
       BaseRBM,
       DBN,
       fit,
       transform,
       generate,
       components,
       save_params,
       load_params

include("core.jl")

end


text:
export TR_ls_step_computation
"""
Compute the direction used in the one dimensionnal trust region scheme.
Return either the Newton direction or the boundary of the trust region if the
Newton direction cannot be used. Inputs:
  - h, value of h(t).
  - g, value of h'(t).
  - dN, value of the Newton direction.
  - Δn, lower bound of the trust region.
  - Δp, upper bound of the trust region.
"""
function TR_ls_step_computation(h :: Float64, g :: Float64, dN :: Float64,
                                Δn :: Float64, Δp ::Float64)

  if h > 0.0
    if g > 0.0
      d = max(Δn, dN)
    else
      d = min(dN, Δp)
    end
  else
    if g > 0.0
      d = Δn
    else
      d = Δp
    end
  end

  return d
end


text:
module problem122

description = """
Efficient exponentiation

The most naive way of computing n15 requires fourteen multiplications:
n × n × ... × n = n15

But using a "binary" method you can compute it in six multiplications:

n × n = n^2
n^2 × n^2 = n^4
n^4 × n^4 = n^8
n^8 × n^4 = n^12
n^12 × n^2 = n^14
n^14 × n = n^15

However it is yet possible to compute it in only five multiplications:

n × n = n^2
n^2 × n = n^3
n^3 × n^3 = n^6
n^6 × n^6 = n^12
n^12 × n^3 = n^15

We shall define m(k) to be the minimum number of multiplications to compute nk; for example m(15) = 5.
For 1 ≤ k ≤ 200, find ∑ m(k).
"""

# Calculates the minimum number of multiplications to compute n^k
# given the known powers, and a maximum cutoff
function muls(k :: Int, known = [1], maxp = k-1)
  if (k % 2) == 0
    return muls(div(k,2)) + 1 # for even powers, solution is always to square the half power

  elseif known[end] == k
    return length(known)-1  # solution found

  elseif length(known) >= maxp # no solution for given known path
    return -1

  elseif (maxp - length(known)) < floor(log2(k / known[end])) # no possible solution for given knowns
    return -1
  end

  # search for solutions by summing items in the knowns path
  for i = length(known):-1:1
    for j = i:-1:1
      n = known[i] + known[j]
      if n <= k && n > known[end]
        push!(known, n)
        result = muls(k, known, maxp)
        pop!(known)
        if result == -1
          break
        elseif result == length(known)
          return result
        elseif result < maxp
          maxp = result
        end
      end
    end
  end

  return maxp
end

function solve(maxk = 200)
  total = 0
  for i = 1:maxk
    total += muls(i)
  end
  return total
end

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:
function mitosis1(which_cells::Array{CartesianIndex{2},1})

    cell_size = size(which_cells)[1]

    centre_of_mass = (Tuple(sum(which_cells)) ./ cell_size)

    xx = 0.
    xy = 0.
    yx = 0.
    yy = 0.

    for i in Tuple.(which_cells)

        xx += (i[1] - centre_of_mass[1])^2
        xy += (i[1] - centre_of_mass[1])*(i[2] - centre_of_mass[2])
        yx += (i[2] - centre_of_mass[2])*(i[1] - centre_of_mass[1])
        yy += (i[2] - centre_of_mass[2])^2
    end

    T = (xx + yy)/cell_size
    D = (xx*yy - xy*yx)/(cell_size^2)

    L1 = 0.5*(T + sqrt(T^2 - 4*D +  1e-10))
    L2 = 0.5*(T - sqrt(T^2 - 4*D + 1e-10))

    if yx != 0
        normal_dir= (-yx/cell_size,max(L1,L2) - yy/cell_size)

    else
        normal_dir = (0,-1)
    end

    hyp_const = -sum(normal_dir .* centre_of_mass)

    left::Array{CartesianIndex{2},1} = []
    right::Array{CartesianIndex{2},1} = []

    for i in which_cells
        if sum(Tuple(i) .* normal_dir) + hyp_const >= 0
            push!(left,i)
        else
            push!(right,i)
        end
    end

    return left,right
end

function cell_mitosis_and_apoptosis!(environment::Array{Float64,3},all_cells::Array{Float64,1},n_div::Dict{Float64,Integer},
    cell_config::CellularConfig)

    new_cells::Array{Float64,1} = []
    counter = 1

    @inbounds for cell in all_cells

        cell_size = 0.
        which_cells::Array{CartesianIndex{2},1} = []

        @inbounds for j in 1:size(environment,2)
            @inbounds for i in 1:size(environment,1)
                if environment[i,j,1] == cell
                    cell_size += 1.
                    push!(which_cells,CartesianIndex(i, j))
                end
            end
        end

        if cell_size >= cell_config.cell_size_upper_limit
            left,right = mitosis1(which_cells)
            new_id = maximum(all_cells) + counter
            environment[left,1] .= new_id
            push!(new_cells,new_id)
            counter += 1
            n_div[cell] += 1
        end
    end

    append!(all_cells,new_cells)

    for cell in new_cells
        n_div[cell] = 0
    end

end

max_stars_repo_name: harrybooth/FractalDevelopment

max_issues_repo_name: harrybooth/FractalDevelopment

max_forks_repo_name: harrybooth/FractalDevelopment


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of colors.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Default color" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data),
                    context = :color => true)

    @test result == expected

    header = ([1, 2, 3, 4],
              [5, 6, 7, 8])

    expected = """
┌───┬───────┬─────┬───┐
│\e[1m 1 \e[0m│\e[1m     2 \e[0m│\e[1m   3 \e[0m│\e[1m 4 \e[0m│
│\e[90m 5 \e[0m│\e[90m     6 \e[0m│\e[90m   7 \e[0m│\e[90m 8 \e[0m│
├───┼───────┼─────┼───┤
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""

    result = sprint((io)->pretty_table(io, data; header = header),
                    context = :color => true)

    @test result == expected
end

@testset "Row number" begin
    expected = """
┌─────┬────────┬────────┬────────┬────────┐
│\e[1m Row \e[0m│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├─────┼────────┼────────┼────────┼────────┤
│   1 │      1 │  false │    1.0 │      1 │
│   2 │      2 │   true │    2.0 │      2 │
│   3 │      3 │  false │    3.0 │      3 │
│   4 │      4 │   true │    4.0 │      4 │
│   5 │      5 │  false │    5.0 │      5 │
│   6 │      6 │   true │    6.0 │      6 │
└─────┴────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data, show_row_number = true),
                    context = :color => true)

    @test result == expected
end

@testset "Row name" begin
    expected = """
┌──────┬────────┬────────┬────────┬────────┐
│\e[1m Name \e[0m│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├──────┼────────┼────────┼────────┼────────┤
│\e[1m    A \e[0m│      1 │  false │    1.0 │      1 │
│\e[1m    B \e[0m│      2 │   true │    2.0 │      2 │
│\e[1m    C \e[0m│      3 │  false │    3.0 │      3 │
│\e[1m    D \e[0m│      4 │   true │    4.0 │      4 │
│\e[1m    E \e[0m│      5 │  false │    5.0 │      5 │
│\e[1m    F \e[0m│      6 │   true │    6.0 │      6 │
└──────┴────────┴────────┴────────┴────────┘
"""

    row_names = ['A'+i for i = 0:5]
    result = sprint((io)->pretty_table(io, data,
                                       row_names = row_names,
                                       row_name_column_title = "Name"),
                    context = :color => true)

    @test result == expected
end

@testset "Border color" begin
   expected = """
\e[33m┌\e[0m\e[33m────────\e[0m\e[33m┬\e[0m\e[33m────────\e[0m\e[33m┬\e[0m\e[33m────────\e[0m\e[33m┬\e[0m\e[33m────────\e[0m\e[33m┐\e[0m
\e[33m│\e[0m\e[1m Col. 1 \e[0m\e[33m│\e[0m\e[1m Col. 2 \e[0m\e[33m│\e[0m\e[1m Col. 3 \e[0m\e[33m│\e[0m\e[1m Col. 4 \e[0m\e[33m│\e[0m
\e[33m├\e[0m\e[33m────────\e[0m\e[33m┼\e[0m\e[33m────────\e[0m\e[33m┼\e[0m\e[33m────────\e[0m\e[33m┼\e[0m\e[33m────────\e[0m\e[33m┤\e[0m
\e[33m│\e[0m      1 \e[33m│\e[0m  false \e[33m│\e[0m    1.0 \e[33m│\e[0m      1 \e[33m│\e[0m
\e[33m│\e[0m      2 \e[33m│\e[0m   true \e[33m│\e[0m    2.0 \e[33m│\e[0m      2 \e[33m│\e[0m
\e[33m│\e[0m      3 \e[33m│\e[0m  false \e[33m│\e[0m    3.0 \e[33m│\e[0m      3 \e[33m│\e[0m
\e[33m│\e[0m      4 \e[33m│\e[0m   true \e[33m│\e[0m    4.0 \e[33m│\e[0m      4 \e[33m│\e[0m
\e[33m│\e[0m      5 \e[33m│\e[0m  false \e[33m│\e[0m    5.0 \e[33m│\e[0m      5 \e[33m│\e[0m
\e[33m│\e[0m      6 \e[33m│\e[0m   true \e[33m│\e[0m    6.0 \e[33m│\e[0m      6 \e[33m│\e[0m
\e[33m└\e[0m\e[33m────────\e[0m\e[33m┴\e[0m\e[33m────────\e[0m\e[33m┴\e[0m\e[33m────────\e[0m\e[33m┴\e[0m\e[33m────────\e[0m\e[33m┘\e[0m
"""

    result = sprint((io)->pretty_table(io, data, border_crayon = crayon"yellow"),
                    context = :color => true)

    @test result == expected
end

@testset "Highlighters" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │\e[33;1m    3.0 \e[0m│      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    hl = Highlighter((data, i, j)-> i == 3 && j == 3;
                     bold = true,
                     foreground = :yellow)

    result = sprint((io)->pretty_table(io, data, highlighters = hl),
                    context = :color => true)

    @test result == expected

    hl = Highlighter((data, i, j)-> i == 3 && j == 3,
                     crayon"yellow bold")

    result = sprint((io)->pretty_table(io, data, highlighters = hl),
                    context = :color => true)

    @test result == expected

    hl2 = Highlighter((data, i, j)-> i == 3 && j == 2,
                      crayon"blue bold")

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │\e[34;1m  false \e[0m│\e[33;1m    3.0 \e[0m│      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data, highlighters = (hl, hl2)),
                    context = :color => true)

    @test result == expected

    hl3 = Highlighter((data, i, j)-> data[i,j] isa AbstractFloat && data[i,j] == 3,
                      (h, data, i, j)->crayon"yellow bold")

    result = sprint((io)->pretty_table(io, data, highlighters = (hl3, hl2)),
                    context = :color => true)

    @test result == expected
end

@testset "Highlighters with table cropping and filters" begin
    matrix = [1:1:100 1:1:100 1:1:100]

    hl = Highlighter((data, i, j)-> i == 100 && j == 2,
                        crayon"yellow")

    expected = """
┌────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│
├────────┼────────┼────────┤
│      1 │      1 │      1 │
│      2 │      2 │      2 │
│      3 │      3 │      3 │
│      4 │      4 │      4 │
│   ⋮    │   ⋮    │   ⋮    │
│     98 │     98 │     98 │
│     99 │     99 │     99 │
│    100 │\e[33m    100 \e[0m│    100 │
└────────┴────────┴────────┘
\e[36m             93 rows omitted\e[0m
"""

    result = sprint((io)->pretty_table(io, matrix,
                                       crop = :both,
                                       display_size = (15,-1),
                                       highlighters = (hl,),
                                       vcrop_mode = :middle),
                    context = :color => true)

    @test result == expected

    expected = """
┌────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│
├────────┼────────┼────────┤
│      2 │      2 │      2 │
│      4 │      4 │      4 │
│      6 │      6 │      6 │
│      8 │      8 │      8 │
│   ⋮    │   ⋮    │   ⋮    │
│     96 │     96 │     96 │
│     98 │     98 │     98 │
│    100 │\e[33m    100 \e[0m│    100 │
└────────┴────────┴────────┘
\e[36m             43 rows omitted\e[0m
"""

    result = sprint((io)->pretty_table(io, matrix,
                                       crop = :both,
                                       display_size = (15,-1),
                                       filters_row = ((data, i)->i % 2 == 0,),
                                       highlighters = (hl,),
                                       vcrop_mode = :middle),
                    context = :color => true)

    @test result == expected

    expected = """
┌────────┬────────┐
│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│
├────────┼────────┤
│      2 │      2 │
│      4 │      4 │
│      6 │      6 │
│      8 │      8 │
│   ⋮    │   ⋮    │
│     96 │     96 │
│     98 │     98 │
│\e[33m    100 \e[0m│    100 │
└────────┴────────┘
\e[36m    43 rows omitted\e[0m
"""

    result = sprint((io)->pretty_table(io, matrix,
                                       crop = :both,
                                       display_size = (15,-1),
                                       filters_col = ((data, i)->i != 1,),
                                       filters_row = ((data, i)->i % 2 == 0,),
                                       highlighters = (hl,),
                                       vcrop_mode = :middle),
                    context = :color => true)

    @test result == expected
end

@testset "Pre-defined highlighters" begin

    # hl_cell
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │\e[33m    3.0 \e[0m│      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data;
                                       highlighters = hl_cell(3, 3, crayon"yellow")),
                    context = :color => true)

    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │\e[33m    2.0 \e[0m│      2 │
│      3 │  false │\e[33m    3.0 \e[0m│      3 │
│      4 │   true │    4.0 │\e[33m      4 \e[0m│
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data;
                                       highlighters = hl_cell([(2, 3), (3, 3), (4, 4)],
                                                              crayon"yellow")),
                    context = :color => true)

    @test result == expected

    # hl_col
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │\e[33m  false \e[0m│    1.0 │      1 │
│      2 │\e[33m   true \e[0m│    2.0 │      2 │
│      3 │\e[33m  false \e[0m│    3.0 │      3 │
│      4 │\e[33m   true \e[0m│    4.0 │      4 │
│      5 │\e[33m  false \e[0m│    5.0 │      5 │
│      6 │\e[33m   true \e[0m│    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data;
                                       highlighters = hl_col(2, crayon"yellow")),
                    context = :color => true)

    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │\e[33m  false \e[0m│    1.0 │\e[33m      1 \e[0m│
│      2 │\e[33m   true \e[0m│    2.0 │\e[33m      2 \e[0m│
│      3 │\e[33m  false \e[0m│    3.0 │\e[33m      3 \e[0m│
│      4 │\e[33m   true \e[0m│    4.0 │\e[33m      4 \e[0m│
│      5 │\e[33m  false \e[0m│    5.0 │\e[33m      5 \e[0m│
│      6 │\e[33m   true \e[0m│    6.0 │\e[33m      6 \e[0m│
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data;
                                       highlighters = hl_col([2,4], crayon"yellow")),
                    context = :color => true)

    @test result == expected

    # hl_row
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│\e[33m      2 \e[0m│\e[33m   true \e[0m│\e[33m    2.0 \e[0m│\e[33m      2 \e[0m│
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data;
                                       highlighters = hl_row(2, crayon"yellow")),
                    context = :color => true)

    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│\e[33m      2 \e[0m│\e[33m   true \e[0m│\e[33m    2.0 \e[0m│\e[33m      2 \e[0m│
│      3 │  false │    3.0 │      3 │
│\e[33m      4 \e[0m│\e[33m   true \e[0m│\e[33m    4.0 \e[0m│\e[33m      4 \e[0m│
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data;
                                       highlighters = hl_row([2,4], crayon"yellow")),
                    context = :color => true)

    @test result == expected

    # hl_lt
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│\e[31;1m      1 \e[0m│\e[31;1m  false \e[0m│\e[31;1m    1.0 \e[0m│\e[31;1m      1 \e[0m│
│\e[31;1m      2 \e[0m│\e[31;1m   true \e[0m│\e[31;1m    2.0 \e[0m│\e[31;1m      2 \e[0m│
│      3 │\e[31;1m  false \e[0m│    3.0 │      3 │
│      4 │\e[31;1m   true \e[0m│    4.0 │      4 │
│      5 │\e[31;1m  false \e[0m│    5.0 │      5 │
│      6 │\e[31;1m   true \e[0m│    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data; highlighters = hl_lt(3)),
                    context = :color => true)

    @test result == expected

    # hl_leq
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│\e[31;1m      1 \e[0m│\e[31;1m  false \e[0m│\e[31;1m    1.0 \e[0m│\e[31;1m      1 \e[0m│
│\e[31;1m      2 \e[0m│\e[31;1m   true \e[0m│\e[31;1m    2.0 \e[0m│\e[31;1m      2 \e[0m│
│\e[31;1m      3 \e[0m│\e[31;1m  false \e[0m│\e[31;1m    3.0 \e[0m│\e[31;1m      3 \e[0m│
│      4 │\e[31;1m   true \e[0m│    4.0 │      4 │
│      5 │\e[31;1m  false \e[0m│    5.0 │      5 │
│      6 │\e[31;1m   true \e[0m│    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data; highlighters = hl_leq(3)),
                    context = :color => true)

    @test result == expected

    # hl_gt
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│\e[34;1m      4 \e[0m│   true │\e[34;1m    4.0 \e[0m│\e[34;1m      4 \e[0m│
│\e[34;1m      5 \e[0m│  false │\e[34;1m    5.0 \e[0m│\e[34;1m      5 \e[0m│
│\e[34;1m      6 \e[0m│   true │\e[34;1m    6.0 \e[0m│\e[34;1m      6 \e[0m│
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data; highlighters = hl_gt(3)),
                    context = :color => true)

    @test result == expected

    # hl_geq
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│\e[34;1m      3 \e[0m│  false │\e[34;1m    3.0 \e[0m│\e[34;1m      3 \e[0m│
│\e[34;1m      4 \e[0m│   true │\e[34;1m    4.0 \e[0m│\e[34;1m      4 \e[0m│
│\e[34;1m      5 \e[0m│  false │\e[34;1m    5.0 \e[0m│\e[34;1m      5 \e[0m│
│\e[34;1m      6 \e[0m│   true │\e[34;1m    6.0 \e[0m│\e[34;1m      6 \e[0m│
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data; highlighters = hl_geq(3)),
                    context = :color => true)

    @test result == expected

    # hl_value
    # ==========================================================================

    expected = """
┌────────┬────────┬────────┬────────┐
│\e[1m Col. 1 \e[0m│\e[1m Col. 2 \e[0m│\e[1m Col. 3 \e[0m│\e[1m Col. 4 \e[0m│
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│\e[33;1m      3 \e[0m│  false │\e[33;1m    3.0 \e[0m│\e[33;1m      3 \e[0m│
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = sprint((io)->pretty_table(io, data; highlighters = hl_value(3)),
                    context = :color => true)

    @test result == expected
end

@testset "Markdown" begin
    # With linebreaks
    # ==========================================================================

    a = md"""
       # Header

       This is a paragraph.

       ```julia
       function test()
           return 1
       end
       ```
       """;

    data = [1 a
            2 a]

    expected = """
┌────────┬────────────────────────┐
│\e[1m Col. 1 \e[0m│\e[1m                 Col. 2 \e[0m│
├────────┼────────────────────────┤
│      1 │ \e[1m  Header\e[22m\e[0m               │
│        │ \e[1m  ≡≡≡≡≡≡≡≡\e[22m\e[0m             │
│        │ \e[0m                       │
│        │   This is a paragraph.\e[0m │
│        │ \e[0m                       │
│        │ \e[36m  function test()\e[39m\e[0m      │
│        │ \e[36m      return 1\e[39m\e[0m         │
│        │ \e[36m  end\e[39m\e[0m                  │
├────────┼────────────────────────┤
│      2 │ \e[1m  Header\e[22m\e[0m               │
│        │ \e[1m  ≡≡≡≡≡≡≡≡\e[22m\e[0m             │
│        │ \e[0m                       │
│        │   This is a paragraph.\e[0m │
│        │ \e[0m                       │
│        │ \e[36m  function test()\e[39m\e[0m      │
│        │ \e[36m      return 1\e[39m\e[0m         │
│        │ \e[36m  end\e[39m\e[0m                  │
└────────┴────────────────────────┘
"""

    result = sprint((io)->pretty_table(io, data,
                                       hlines = :all,
                                       linebreaks = true),
                    context = :color => true)

    @test result == expected

    # Without linebreaks
    # ==========================================================================

    a = md"""
    **bold**
    *italics*
    """

    data = [1 a]

    expected = """
┌────────┬────────────────┐
│\e[1m Col. 1 \e[0m│\e[1m         Col. 2 \e[0m│
├────────┼────────────────┤
│      1 │   \e[1mbold\e[22m \e[4mitalics\e[24m │
└────────┴────────────────┘
"""

    result = sprint((io)->pretty_table(io, data,
                                       hlines = :all,
                                       linebreaks = false),
                    context = :color => true)

    @test result == expected
end


text:
using Plots, SNN

SNN.IF_curve(SNN.IF, 0:0.1:20) |> display
SNN.IF_curve(SNN.IZ, 0:0.1:20) |> display
SNN.IF_curve(SNN.HH, 1e-5*(0:0.01:4)) |> display


text:
#
# time utility functions
#
function default_dtmod_helper(stenciltype::Symbol, Q::Int, D::Int, mindel::T, velmax::T, dtmod::T, dtrec::T, alpha::T, ndim::Int) where T<:Real
    a = stencilcoeffs(:Fornberg, Q, D, ndim, dtmod, velmax, mindel, T)
    c = courantnumber(a, ndim)
    default_dtmod_helper_helper(mindel, velmax, c, dtmod, dtrec, alpha)
end
function default_dtmod_helper_helper(mindel::T, velmax::T, c::T, dtmod::T, dtrec::T, alpha::T) where T<:Real
    dtmod_auto = alpha*c*mindel/velmax
    if dtmod < 0
        dtmod = dtmod_auto
    end
    if 0.99*dtmod > dtmod_auto
        throw(ArgumentError("dtmod=$(dtmod), dtmod_auto=$(dtmod_auto), mindel=$(mindel), velmax=$(velmax)"))
    end
    T(dtrec/ceil(dtrec/min(dtmod,dtrec) - .1)) # tricky
end

"""
# 2D
    dtmod = WaveFD.default_dtmod(stenciltype, Q, D, dz, dx, velmax, dtmod, dtrec, alpha)
    dtmod = WaveFD.default_dtmod(alpha, dz, dx, velmax, dtrec, dtmod) # assumes courant=1 .. should go away once we become less lazy about doing analysis
"""
default_dtmod(stenciltype::Symbol, Q::Int, D::Int, dz::T, dx::T, velmax::T, dtmod::T, dtrec::T, alpha::T) where {T<:AbstractFloat} = default_dtmod_helper(stenciltype, Q, D, min(dz,dx), velmax, dtmod, dtrec, alpha, 2)
default_dtmod(alpha::T, dz::T, dx::T, velmax::T, dtrec::T, dtmod::T) where {T<:AbstractFloat} = default_dtmod_helper_helper(min(dz,dx), velmax, one(T), dtmod, dtrec, alpha)

"""
# 3D
    dtmod = WaveFD.default_dtmod(stenciltype, Q, D, dz, dy, dx, velmax, dtmod, dtrec, alpha)
    dtmod = WaveFD.default_dtmod(alpha, dz, dy, dx, velmax, dtrec, dtmod) # assumes courant=1, lazy option, rather than doing proper analysis

where `dtmod::AbstractFloat` is a stable finite-difference time step determined by the C-F-L condition (https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition).
A modified Courant number is used for the C-F-L condition, and is determined by the choice of finite-difference stencil and a scaling factor `alpha`.  Making `alpha` smaller should reduce
numerical dispersion, but at the expense of a smaller `dtmod`.  We have found that setting `alpha=0.25` produces results with small dispersion.  For simplicity, we compute `dtmod`
such that `dtrec` is a scalar multiple of `dtmod`.

## paramters
* `proptype` choose betweeen `:O28_C`,`:O28_Julia`,`:O24_C`, or `:O24_Julia`.  The Courant number depends on the finite difference stencil.  In this case, either 4th or 8th order
* `stenciltype` choose between `:Fornberg` and `:Nihei`.  `:Fornberg` are standard coefficients derived from Taylor. `:Nihei` is an attempt at producing optimal coefficients but only work for constant earth models.
* `Q` is the accuracy order of the stencil, choose between `4` and `8`
* `D` is the derivative order of the stencil, choose between `1` and `2`. Note that `1` assumes staggered, and `2` assumes centered.
* `dz`,`dy`,`dx` are the grid spacing
* `velmax` expected maximum velocity in the earth model
* `dtmod` desired modeling time interval.  If `dtmod<0`, then `dtmod` will be computed according to C-F-L and the recording dt.  Otherwise, we will attempt to use the user specified `dtmod` and throw
an error if it violates C-F-L condition.  However, we will, if needed, reduce `dtmod` so that `dtrec` is a scalar multiple of `dtmod`.
* `dtrec` recording time interval
* `alpha` Courant number multiplier

## Example
```julia
dtmod=-1.0
dz,dx,velmax,dtrec,alpha=0.01,0.01,5.0,0.004,0.25
Q,D=8,1
dtmod = WaveFD.default_dtmod(:Fornberg, Q, D, dz, dx, velmax, dtmod, dtrec, alpha)
```

## Notes
This method is not exported, because it is not the expectation that the user will need to call this method unless they want/need to by-pass the logic in the
Jot Wavefield propagation operators.
"""
default_dtmod(stenciltype::Symbol, Q::Int, D::Int, dz::T, dy::T, dx::T, velmax::T, dtmod::T, dtrec::T, alpha::T) where {T<:AbstractFloat} = default_dtmod_helper(stenciltype, Q, D, min(dz,dy,dx), velmax, dtmod, dtrec, alpha, 3)
default_dtmod(alpha::T, dz::T, dy::T, dx::T, velmax::T, dtrec::T, dtmod::T) where {T<:AbstractFloat} = default_dtmod_helper_helper(min(dz,dy,dx), velmax, one(T), dtmod, dtrec, alpha)

"""
    it0, ntmod = WaveFD.default_ntmod(dtrec,dtmod,t0,ntrec)

`it0` the time index corresponding to zero-time, and `ntmod` the total number of time points, including negative time.
`dtrec`,`dtmod` and `ntrec` are `Real`, while `t0` can either be `Real`, `Array{Array{Real}}` or `DArray{Array{Real}}`.  If `t0` is
of type `Array{Real}`, then `it0` is determined using `minimum(t0)`.
"""
function default_ntmod(dtrec::Real, dtmod::Real, t0::Real, ntrec::Int)
    ntmod = round(Int,dtrec/dtmod)*(ntrec - 1) + 1
    it0 = max(round(Int64, -t0/dtmod) + 1, 1)
    ntmod += it0-1
    it0, ntmod
end
default_ntmod(dtrec::Real, dtmod::Real, t0::Array, ntrec::Int) = default_ntmod(dtrec, dtmod, mapreduce(i->minimum(t0[i]), min, 1:length(t0)), ntrec)

function default_ntmod(dtrec::Real, dtmod::Real, t0::DArray, ntrec::Int)
    pids = procs(t0)
    npids = length(pids)
    r = Array{Future}(undef, npids)
    @sync for (i,pid) in enumerate(pids)
        r[i] = @spawnat pid begin
            t0_local = localpart(t0)
            t0min_local = minimum(t0_local[1])
            for j = 2:length(t0_local)
                t0min_local = minimum(t0_local[j]) < t0min_local ? minimum(t0_local[j]) : t0min_local
            end
            t0min_local
        end
    end
    t0min = fetch(r[1])
    for i = 2:npids
        t0test = fetch(r[i])
        t0min = t0test < t0min ? t0test : t0min
    end
    default_ntmod(dtrec, dtmod, t0min, ntrec)
end

"""
    ntmod = default_ntmod(dtrec,dtmod,ntrec)

`ntmod` is the total number of time points.  It is assuming that the first time index corresponds to zero-time.
"""
default_ntmod(dtrec::Real, dtmod::Real, ntrec::Int) = default_ntmod(dtrec, dtmod, 0.0, ntrec)[2]

#
# time interpolation codes
#

mutable struct TimeInterp{T<:AbstractFloat,C<:Language}
    h::Array{Array{T,1},1}
    nthreads::Int64
end

"""
    h = WaveFD.interpfilters(dtmod, dtrec [, mode=0, impl=WaveFD.LangC, nthreads=Sys.CPU_THREADS])

Build an 8 point sinc filter, mapping between `dtmod` and `dtrec` sampling.  The optional parameters are:

* `mode::Int` determines if amplitude is preserved in the forward (mode=0) or adjoint operation (mode=1)
* `impl::WaveFD.Language` can be set to either `WaveFD.LangJulia()` or `WaveFD.LangC()` to determine which code-path to follow
* `nthreads`  If `impl=WaveFD.LangC()`, then OMP is used for threading the 2D and 3D arrays where time is assumed to be along the fast dimension

# Notes
It is assumed that `dtmod<dtrec`:

* forward operator (see `WaveFD.interpforward!`) - interpolates from `dtmod` to `dtrec`, preserves amplitude if `mode=0`.  This is the default behaviour.
* adjoint operator (see `WaveFD.interpadjoint!`) - interpolates from `dtrec` to `dtmod`, preserves amplitude if `mode=1`
"""
function interpfilters(dtmod::T, dtrec::T, mode::Int=0, impl::Language=LangC(), nthreads=Sys.CPU_THREADS) where T<:Real
    @assert dtmod <= dtrec
    n = round(Int, dtrec/dtmod)
    if abs(n*dtmod - dtrec)/dtrec > 1e-6
        throw(ArgumentError("abs(n*dtmod-dtrec)/dtrec > 1e-6, n=$(n), dtrec=$(dtrec), dtmod=$(dtmod), abs(n*dtmod-dtrec)/dtrec=$(abs(n*dtmod-dtrec)/n)"))
    end

    d = dtmod/dtrec*[0:n-1;]
    fmax = T(0.066 + 0.265*log(8.0))
    j = [1:8;]

    a = @. sinc((j-1)*fmax)
    A = [a[1] a[2] a[3] a[4] a[5] a[6] a[7] a[8] ;
         a[2] a[1] a[2] a[3] a[4] a[5] a[6] a[7] ;
         a[3] a[2] a[1] a[2] a[3] a[4] a[5] a[6] ;
         a[4] a[3] a[2] a[1] a[2] a[3] a[4] a[5] ;
         a[5] a[4] a[3] a[2] a[1] a[2] a[3] a[4] ;
         a[6] a[5] a[4] a[3] a[2] a[1] a[2] a[3] ;
         a[7] a[6] a[5] a[4] a[3] a[2] a[1] a[2] ;
         a[8] a[7] a[6] a[5] a[4] a[3] a[2] a[1]]

    h = Array{T}[]
    for i = 1:n
        c = @. sinc((4-j+d[i])*fmax)
        push!(h, A\c)
        if mode == 0
            h[i] ./= n # denom is to make sure that the amplitude in the forward is correct
        end
    end

    TimeInterp{T,typeof(impl)}(h,nthreads)
end

"""
    WaveFD.interpadjoint!(h, m, d)

Interpolate from `d::Array{T,N}` (coarse) to `m::Array{T,N}` (fine) using the sinc filter coefficients in `h::Array{Array{T,1},1}`.  `h` is built using
`WaveFD.interpfilters`.  For example:

    WaveFD.interpadjoint!(WaveFD.interpfilters(.001,.004), m, d)

Note that we support, `N=1`, `N=2` or `N=3`. If `N=2` or `N=3`, then interpolation is done along the fast dimension.  By default, `interpadjoint!` does not preserve
amplitude (see `WaveFD.interpfilters`).
"""
function interpadjoint!(H::TimeInterp{T,LangJulia}, m::StridedArray{T,1}, d::StridedArray{T,1}) where T
    if length(m) == length(d)
        m .= d
        return nothing
    end
    h = H.h
    nh = length(h)
    nd = length(d)
    nm = length(m)
    fill!(m,0.0)

    interpadjoint_helper!(m,d,h,nm,nh,1,4,8) # i=1
    interpadjoint_helper!(m,d,h,nm,nh,2,3,8) # i=2
    interpadjoint_helper!(m,d,h,nm,nh,3,2,8) # i=3
    for i = 4:nd-4
        interpadjoint_helper!(m,d,h,nm,nh,i,1,8)
    end
    interpadjoint_helper!(m,d,h,nm,nh,nd-3,1,7) # i=nd-3
    interpadjoint_helper!(m,d,h,nm,nh,nd-2,1,6) # i=nd-2
    interpadjoint_helper!(m,d,h,nm,nh,nd-1,1,5) # i=nd-1
    nothing
end

@inline function interpadjoint_helper!(m,d,h,nm,nh,i,lo,hi)
    j = (i-1)*nh
    im4 = i-4
    for k=1:nh
        jk = j+k
        jk > nm && continue
        for l = lo:hi
            m[jk] += h[k][l] * d[im4+l]
        end
    end
end

"""
    WaveFD.interpforward!(h, d, m)

Interpolate from `m::Array{T,N}` (fine) to `d::Array{T,N}` (coarse) using the sinc filter coefficients in `h::Array{Array{T,1},1}`.  `h` is built using
`WaveFD.interpfilters`.  For example:

    WaveFD.interpadjoint!(WaveFD.interpfilters(.001,.004), m, d)

Note that we support, `N=1`, `N=2` or `N=3`. If `N=2` or `N=3`, then interpolation is done along the fast dimension.  By default, `interpforward!` preserves
amplitude (see `WaveFD.interpfilters!`).
"""
function interpforward!(H::TimeInterp{T,LangJulia}, d::StridedArray{T,1}, m::StridedArray{T,1}) where T
    if length(d) == length(m)
        d .= m
        return nothing
    end
    h = H.h
    nh = length(h)
    nd = length(d)
    nm = length(m)
    fill!(d,0.0)
    interpforward_helper!(d,m,h,nh,1,4,8) # i=1
    interpforward_helper!(d,m,h,nh,2,3,8) # i=2
    interpforward_helper!(d,m,h,nh,3,2,8) # i=3
    for i = 4:nd-4
        interpforward_helper!(d,m,h,nh,i,1,8)
    end
    interpforward_helper!(d,m,h,nh,nd-3,1,7) # i=nd-3
    interpforward_helper!(d,m,h,nh,nd-2,1,6) # i=nd-2
    interpforward_helper!(d,m,h,nh,nd-1,1,5) # i=nd-1
    nothing
end

@inline function interpforward_helper!(d,m,h,nh,i,lo,hi)
    j = (i-1)*nh
    im4 = i-4
    for k = 1:nh
        jk = j + k
        jk > length(m) && continue
        for l = lo:hi
            d[im4+l] += h[k][l]*m[jk]
        end
    end
end

function interpforward!(H::TimeInterp{T,LangJulia}, d::StridedArray{T,2}, m::StridedArray{T,2}) where T
    for i2 = 1:size(m,2)
        d_trace = @view d[:,i2]
        m_trace = @view m[:,i2]
        interpforward!(H, d_trace, m_trace)
    end
    nothing
end

function interpforward!(H::TimeInterp{T,LangJulia}, d::StridedArray{T,3}, m::StridedArray{T,3}) where T
    for i3 = 1:size(m,3), i2 = 1:size(m,2)
        d_trace = @view d[:,i2,i3]
        m_trace = @view m[:,i2,i3]
        interpforward!(H, d_trace, m_trace)
    end
    nothing
end

function interpadjoint!(H::TimeInterp{T,LangJulia}, m::StridedArray{T,2}, d::StridedArray{T,2}) where T
    for i2 = 1:size(m,2)
        d_trace = @view d[:,i2]
        m_trace = @view m[:,i2]
        interpadjoint!(H, m_trace, d_trace)
    end
    nothing
end

function interpadjoint!(H::TimeInterp{T,LangJulia}, m::StridedArray{T,3}, d::StridedArray{T,3}) where T
    for i3 = 1:size(m,3), i2=1:size(m,2)
        d_trace = @view d[:,i2,i3]
        m_trace = @view m[:,i2,i3]
        interpadjoint!(H, m_trace, d_trace)
    end
    nothing
end

#
# fast-path (OpenMP) codes for time interpolation (should be able to delete these once Julia gets shared memory parallel support)
#
function interpadjoint!(H::TimeInterp{Float32,LangC}, m::StridedArray{Float32,1}, d::StridedArray{Float32,1})
    if length(m) == length(d)
        m .= d
        return nothing
    end
    ccall(
        (:interpadjoint_1d_float, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cfloat}}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t,     Csize_t,   Csize_t),
         H.h,              m,           d,           length(H.h), length(m), length(d)
    )
end

function interpadjoint!(H::TimeInterp{Float64,LangC}, m::StridedArray{Float64,1}, d::StridedArray{Float64,1})
    if length(m) == length(d)
        m .= d
        return nothing
    end
    ccall(
        (:interpadjoint_1d_double, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cdouble}, Csize_t,     Csize_t,   Csize_t),
         H.h,               m,            d,            length(H.h), length(m), length(d)
    )
end

function interpadjoint!(H::TimeInterp{Float32,LangC}, m::StridedArray{Float32,2}, d::StridedArray{Float32,2})
    if length(m) == length(d)
        m .= d
        return nothing
    end
    ccall(
        (:interpadjoint_nd_float, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cfloat}}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t,     Csize_t,   Csize_t,   Csize_t,   Csize_t),
         H.h,              m,           d,           length(H.h), size(m,1), size(d,1), size(d,2), H.nthreads
    )
    nothing
end

function interpadjoint!(H::TimeInterp{Float64,LangC}, m::StridedArray{Float64,2}, d::StridedArray{Float64,2})
    if length(m) == length(d)
        m .= d
        return nothing
    end
    ccall(
        (:interpadjoint_nd_double, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cdouble}, Csize_t,     Csize_t,   Csize_t,   Csize_t,   Csize_t),
         H.h,               m,            d,            length(H.h), size(m,1), size(d,1), size(d,2), H.nthreads
    )
    nothing
end

function interpadjoint!(H::TimeInterp{Float32,LangC}, m::StridedArray{Float32,3}, d::StridedArray{Float32,3})
    if length(m) == length(d)
        m .= d
        return nothing
    end
    ccall(
        (:interpadjoint_nd_float, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cfloat}}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t,     Csize_t,   Csize_t,   Csize_t,             Csize_t),
         H.h,              m,           d,           length(H.h), size(m,1), size(d,1), size(d,2)*size(d,3), H.nthreads
    )
    nothing
end

function interpadjoint!(H::TimeInterp{Float64,LangC}, m::StridedArray{Float64,3}, d::StridedArray{Float64,3})
    if length(m) == length(d)
        m .= d
        return nothing
    end
    ccall(
        (:interpadjoint_nd_double, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cdouble}, Csize_t,     Csize_t,   Csize_t,   Csize_t,             Csize_t),
         H.h,               m,            d,            length(H.h), size(m,1), size(d,1), size(d,2)*size(d,3), H.nthreads
    )
    nothing
end

function interpforward!(H::TimeInterp{Float32,LangC}, d::StridedArray{Float32,1}, m::StridedArray{Float32,1})
    if length(d) == length(m)
        d .= m
        return nothing
    end
    ccall(
        (:interpforward_1d_float, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cfloat}}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t,     Csize_t,   Csize_t),
         H.h,              d,           m,           length(H.h), length(d), length(m)
    )
    nothing
end

function interpforward!(H::TimeInterp{Float64,LangC}, d::StridedArray{Float64,1}, m::StridedArray{Float64,1})
    if length(d) == length(m)
        d .= m
        return nothing
    end
    ccall(
        (:interpforward_1d_double, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cdouble}, Csize_t,     Csize_t,   Csize_t),
         H.h,               d,            m,            length(H.h), length(d), length(m)
    )
    nothing
end

function interpforward!(H::TimeInterp{Float32,LangC}, d::StridedArray{Float32,2}, m::StridedArray{Float32,2})
    if length(d) == length(m)
        d .= m
        return nothing
    end
    ccall(
        (:interpforward_nd_float, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cfloat}}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t,      Csize_t,   Csize_t,   Csize_t,   Csize_t),
         H.h,              d,           m,           length(H.h),  size(d,1), size(m,1), size(m,2), H.nthreads
    )
    nothing
end

function interpforward!(H::TimeInterp{Float64,LangC}, d::StridedArray{Float64,2}, m::StridedArray{Float64,2})
    if length(d) == length(m)
        d .= m
        return nothing
    end
    ccall(
        (:interpforward_nd_double, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cdouble}, Csize_t,     Csize_t,   Csize_t,   Csize_t,   Csize_t),
         H.h,               d,            m,            length(H.h), size(d,1), size(m,1), size(m,2), H.nthreads
    )
    nothing
end

function interpforward!(H::TimeInterp{Float32,LangC}, d::StridedArray{Float32,3}, m::StridedArray{Float32,3})
    if length(d) == length(m)
        d .= m
        return nothing
    end
    ccall(
        (:interpforward_nd_float, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cfloat}}, Ptr{Cfloat}, Ptr{Cfloat}, Csize_t,      Csize_t,   Csize_t,   Csize_t,             Csize_t),
         H.h,               d,           m,           length(H.h), size(d,1), size(m,1), size(m,2)*size(m,3), H.nthreads
    )
    nothing
end

function interpforward!(H::TimeInterp{Float64,LangC}, d::StridedArray{Float64,3}, m::StridedArray{Float64,3})
    if length(d) == length(m)
        d .= m
        return nothing
    end
    ccall(
        (:interpforward_nd_double, libspacetime),
        Ptr{Cvoid},
        (Ptr{Ptr{Cdouble}}, Ptr{Cdouble}, Ptr{Cdouble}, Csize_t,     Csize_t,   Csize_t,   Csize_t,             Csize_t),
         H.h,               d,            m,            length(H.h), size(d,1), size(m,1), size(m,2)*size(m,3), H.nthreads
    )
    nothing
end

#
# space interpolation codes
#
#
# Graham Hicks' Kaiser windowed sinc function method, Geophysics 67(1), P. 156-166
# paper suggests r=4, with b=6.31 for kmax=pi/2 and with b=4.14 for kmax=2pi/3
#
# For our problem we are using 5 grid points per wave-length --
# k_max = w_max / c_min = 2 pi fmax / c_min
# k_nyq = pi / del = pi / (c_min / (5fmax)) = 5 fmax * pi / c_min
#
# k_max   2 pi fmax / c_min   2   4        1   5
# ----- = ----------------- = - = --    <  - = --
# k_nyq   5 pi fmax / c_min   5   10       2   10
#
# Hicks suggests b=6.42 for kmax=1/2q
kaiser(b::T, r::Int64, x::T) where {T<:Union{Float32,Float64}} = abs(x) > r ? zero(T) : T(besseli(0,b*sqrt(1-(x/r)^2))/besseli(0,b))

function hicks!(d::Matrix, b::Real, r::Integer, alphaz::Real, alphax::Real, fz::Vector, fx::Vector)
    l = 2*r
    @inbounds for i = 1:l
        a = i - r
        z = a - alphaz
        x = a - alphax
        fz[i] = kaiser(b,r,z)*sinc(z)
        fx[i] = kaiser(b,r,x)*sinc(x)
    end
    @inbounds for j = 1:l, i = 1:l
        d[i,j] = fz[i]*fx[j]
    end
end

function hicks(b::Real, r::Integer, alphaz::Real, alphax::Real)
    l = 2*r
    d = Array{eltype(b)}(undef, l, l)
    fz = Array{eltype(b)}(undef, l)
    fx = Array{eltype(b)}(undef, l)
    hicks!(d, b, r, alphaz, alphax, fz, fx)
    d
end

function hicks!(d::Array{T,3}, b::Real, r::Integer, alphaz::Real, alphay::Real, alphax::Real, fz::Vector, fy::Vector, fx::Vector) where T<:Real
    l = 2*r
    @inbounds for i = 1:l
        a = i - r
        z = a - alphaz
        y = a - alphay
        x = a - alphax
        fz[i] = kaiser(b,r,z)*sinc(z)
        fy[i] = kaiser(b,r,y)*sinc(y)
        fx[i] = kaiser(b,r,x)*sinc(x)
    end
    @inbounds for k = 1:l, j = 1:l, i = 1:l
        d[i,j,k] = fz[i]*fy[j]*fx[k]
    end
end

function hicks(b::Real, r::Integer, alphaz::Real, alphay::Real, alphax::Real)
    l = 2*r
    d = Array{eltype(b)}(undef, l, l, l)
    fz = Array{eltype(b)}(undef, l)
    fy = Array{eltype(b)}(undef, l)
    fx = Array{eltype(b)}(undef, l)
    hicks!(d, b, r, alphaz, alphay, alphax, fz, fy, fx)
    d
end

function hickscoeffs(dz::T, dx::T, z0::Float64, x0::Float64, nz::Int64, nx::Int64, z::Array{Float64,1}, x::Array{Float64,1}, fs_index=1, zstagger=0.0, xstagger=0.0) where T
    hicks_b = T(6.42)
    hicks_r = 4
    hicks_l = 2*hicks_r
    hicks_fz = Array{T}(undef, hicks_l)
    hicks_fx = Array{T}(undef, hicks_l)
    iz = Array{Array{Int64,1}}(undef, length(z))
    ix = Array{Array{Int64,1}}(undef, length(z))
    c = Array{Array{T,2}}(undef, length(z))
    zerotol = Base.eps(T)*1e3
    for i = 1:length(z)
        iz_mid = floor(Int, (z[i] - z0)/dz) + 1
        ix_mid = floor(Int, (x[i] - x0)/dx) + 1
        alphaz = T((z[i] - (z0 + (iz_mid - 1 + zstagger)*dz))/dz)
        alphax = T((x[i] - (x0 + (ix_mid - 1 + xstagger)*dx))/dx)

        ongrid = abs(alphaz) < zerotol && abs(alphax) < zerotol
        if ongrid == true
            iz[i] = [iz_mid]
            ix[i] = [ix_mid]
            c[i] = ones(T,1,1)
        else
            iz[i] = collect(iz_mid-hicks_r+1:iz_mid+hicks_r)
            ix[i] = collect(ix_mid-hicks_r+1:ix_mid+hicks_r)
            c[i] = Array{T}(undef, hicks_l, hicks_l)
            hicks!(c[i], hicks_b, hicks_r, alphaz, alphax, hicks_fz, hicks_fx)

            # reflect source terms above the free-surface
            # we assume that if there is no free-surface, then
            # the model would have being padded accordingly, making
            # the below if statement always evaluate to false
            if iz[i][1] == fs_index - 3
                c[i][5,:] -= c[i][3,:]
                c[i][6,:] -= c[i][2,:]
                c[i][7,:] -= c[i][1,:]
            end
            if iz[i][1] == fs_index - 2
                c[i][4,:] -= c[i][2,:]
                c[i][5,:] -= c[i][1,:]
            end
            if iz[i][1] == fs_index - 1
                c[i][3,:] -= c[i][1,:]
            end

            # cut the coefficients beyond the edge of the domain:
            iz_in = findall(z->(fs_index+1)<=z<=nz,iz[i]) # exclude the free-surface since its coefficients should always be zero
            ix_in = findall(x->1<=x<=nx,ix[i])
            c[i] = c[i][iz_in,ix_in]
            iz[i] = iz[i][iz_in]
            ix[i] = ix[i][ix_in]
        end
    end
    iz, ix, c
end

function hickscoeffs(dz::T, dy::T, dx::T, z0::Float64, y0::Float64, x0::Float64, nz::Int64, ny::Int64, nx::Int64, z::Array{Float64,1}, y::Array{Float64,1}, x::Array{Float64,1}, fs_index=1, zstagger=0.0, ystagger=0.0, xstagger=0.0) where T
    hicks_b = T(6.42)
    hicks_r = 4
    hicks_l = 2*hicks_r
    hicks_fz = Array{T}(undef, hicks_l)
    hicks_fy = Array{T}(undef, hicks_l)
    hicks_fx = Array{T}(undef, hicks_l)
    iz = Array{Array{Int64,1}}(undef, length(z))
    iy = Array{Array{Int64,1}}(undef, length(z))
    ix = Array{Array{Int64,1}}(undef, length(z))
    c = Array{Array{T,3}}(undef, length(z))
    zerotol = Base.eps(T)*1e3
    for i = 1:length(z)
        iz_mid = floor(Int, (z[i] - z0)/dz) + 1
        iy_mid = floor(Int, (y[i] - y0)/dy) + 1
        ix_mid = floor(Int, (x[i] - x0)/dx) + 1
        alphaz = T((z[i] - (z0 + (iz_mid - 1 + zstagger)*dz))/dz)
        alphay = T((y[i] - (y0 + (iy_mid - 1 + ystagger)*dy))/dy)
        alphax = T((x[i] - (x0 + (ix_mid - 1 + xstagger)*dx))/dx)

        ongrid = abs(alphaz) < zerotol && abs(alphay) < zerotol && abs(alphax) < zerotol
        if ongrid == true
            iz[i] = [iz_mid]
            iy[i] = [iy_mid]
            ix[i] = [ix_mid]
            c[i] = ones(T,1,1,1)
        else
            iz[i] = collect(iz_mid-hicks_r+1:iz_mid+hicks_r)
            iy[i] = collect(iy_mid-hicks_r+1:iy_mid+hicks_r)
            ix[i] = collect(ix_mid-hicks_r+1:ix_mid+hicks_r)
            c[i] = Array{T}(undef, hicks_l, hicks_l, hicks_l)
            hicks!(c[i], hicks_b, hicks_r, alphaz, alphay, alphax, hicks_fz, hicks_fy, hicks_fx)

            # reflect source terms above the free-surface
            # we assume that if there is no free-surface, then
            # the model would have being padded accordingly, making
            # the below if statement always evaluate to false
            if iz[i][1] == fs_index-3
                c[i][5,:,:] -= c[i][3,:,:]
                c[i][6,:,:] -= c[i][2,:,:]
                c[i][7,:,:] -= c[i][1,:,:]
            end
            if iz[i][1] == fs_index-2
                c[i][4,:,:] -= c[i][2,:,:]
                c[i][5,:,:] -= c[i][1,:,:]
            end
            if iz[i][1] == fs_index-1
                c[i][3,:,:] -= c[i][1,:,:]
            end

            # cut the coefficients beyond the edge of the domain:
            iz_in = findall(z->(fs_index+1)<=z<=nz,iz[i]) # exclude the free-surface since its coefficients should always be zero
            iy_in = findall(y->1<=y<=ny,iy[i])
            ix_in = findall(x->1<=x<=nx,ix[i])
            c[i] = c[i][iz_in,iy_in,ix_in]
            iz[i] = iz[i][iz_in]
            iy[i] = iy[i][iy_in]
            ix[i] = ix[i][ix_in]
        end
    end
    iz, iy, ix, c
end

function linearcoeffs(dz::T, dx::T, z0::Float64, x0::Float64, nz::Int64, nx::Int64, z::Array{Float64,1}, x::Array{Float64,1}, fs_index=1) where T
    iz = Array{Array{Int64,1}}(undef, length(z))
    ix = Array{Array{Int64,1}}(undef, length(z))
    c = Array{Array{T,2}}(undef, length(z))
    for i = 1:length(z)
        kx = floor(Int64, (x[i] - x0) / dx) + 1
        kz = floor(Int64, (z[i] - z0) / dz) + 1

        rx = (x[i] - (x0 + dx*(kx-1))) / dx
        rz = (z[i] - (z0 + dz*(kz-1))) / dz

        c[i] = Array{T}(undef, 2,2)
        c[i][1,1] = (1.0 - rx) * (1.0 - rz)
        c[i][2,1] = (1.0 - rx) *        rz
        c[i][1,2] =        rx  * (1.0 - rz)
        c[i][2,2] =        rx  *        rz

        iz[i] = collect(kz:kz+1)
        ix[i] = collect(kx:kx+1)

        # cut the coefficients beyond the edge of the domain:
        iz_in = findall(z->(fs_index+1)<=z<=nz,iz[i]) # exclude the free-surface since its coefficients should always be zero
        ix_in = findall(x->1<=x<=nx,ix[i])
        c[i] = c[i][iz_in,ix_in]
        iz[i] = iz[i][iz_in]
        ix[i] = ix[i][ix_in]
    end
    iz, ix, c
end

function linearcoeffs(dz::T, dy::T, dx::T, z0::Float64, y0::Float64, x0::Float64, nz::Int64, ny::Int64, nx::Int64, z::Array{Float64,1}, y::Array{Float64,1}, x::Array{Float64,1}, fs_index=1) where T
    iz = Array{Array{Int64,1}}(undef, length(z))
    iy = Array{Array{Int64,1}}(undef, length(z))
    ix = Array{Array{Int64,1}}(undef, length(z))
    c = Array{Array{T,3}}(undef, length(z))
    for i = 1:length(z)
        kx = floor(Int64, (x[i] - x0) / dx) + 1
        ky = floor(Int64, (y[i] - y0) / dy) + 1
        kz = floor(Int64, (z[i] - z0) / dz) + 1

        rx = (x[i] - (x0 + dx*(kx-1))) / dx
        ry = (y[i] - (y0 + dy*(ky-1))) / dy
        rz = (z[i] - (z0 + dz*(kz-1))) / dz

        c[i] = Array{T}(undef, 2,2,2)
        c[i][1,1,1] = (1.0 - rx) * (1.0 - ry) * (1.0 - rz)
        c[i][2,1,1] = (1.0 - rx) * (1.0 - ry) *        rz
        c[i][1,2,1] = (1.0 - rx) *        ry  * (1.0 - rz)
        c[i][2,2,1] = (1.0 - rx) *        ry  *        rz
        c[i][1,1,2] =        rx  * (1.0 - ry) * (1.0 - rz)
        c[i][2,1,2] =        rx  * (1.0 - ry) *        rz
        c[i][1,2,2] =        rx  *        ry  * (1.0 - rz)
        c[i][2,2,2] =        rx  *        ry  *        rz

        iz[i] = collect(kz:kz+1)
        iy[i] = collect(ky:ky+1)
        ix[i] = collect(kx:kx+1)

        # cut the coefficients beyond the edge of the domain:
        iz_in = findall(z->(fs_index+1)<=z<=nz,iz[i]) # exclude the free-surface since its coefficients should always be zero
        iy_in = findall(y->1<=y<=ny,iy[i])
        ix_in = findall(x->1<=x<=nx,ix[i])
        c[i] = c[i][iz_in,iy_in,ix_in]
        iz[i] = iz[i][iz_in]
        iy[i] = iy[i][iy_in]
        ix[i] = ix[i][ix_in]
    end
    iz, iy, ix, c
end

function ongridcoeffs(dz::T, dx::T, z0::Float64, x0::Float64, nz::Int64, nx::Int64, z::Array{Float64,1}, x::Array{Float64,1}) where T
    iz = Array{Int64}(undef, length(z))
    ix = Array{Int64}(undef, length(z))
    c  = Array{T}(undef, length(z))
    for i = 1:length(z)
        iz[i] = round(Int64, (z[i] - z0)/dz) + 1
        ix[i] = round(Int64, (x[i] - x0)/dx) + 1
        c[i]  = one(T)
    end
    iz, ix, c
end

function ongridcoeffs(dz::T, dy::T, dx::T, z0::Float64, y0::Float64, x0::Float64, nz::Int64, ny::Int64, nx::Int64, z::Array{Float64,1}, y::Array{Float64,1}, x::Array{Float64,1}) where T
    iz = Array{Int64}(undef, length(z))
    iy = Array{Int64}(undef, length(z))
    ix = Array{Int64}(undef, length(z))
    c  = Array{T}(undef, length(z))
    for i = 1:length(z)
        iz[i] = round(Int64, (z[i] - z0)/dz) + 1
        iy[i] = round(Int64, (y[i] - y0)/dy) + 1
        ix[i] = round(Int64, (x[i] - x0)/dx) + 1
        c[i]  = one(T)
    end
    iz, iy, ix, c
end

function sourceblock_range(nb, irng)
    b,rm = divrem(length(irng),nb)
    blocks = Vector{UnitRange{Int}}(undef, nb)
    i1 = irng[1]
    for iblock = 1:nb
        i2 = i1 + b - 1
        if iblock <= rm
            i2 += 1
        end
        blocks[iblock] = i1:i2
        i1 = i2+1
    end
    blocks
end

function get_blocks_for_source_point(i, blocks)
    iblock_min = findfirst(block->in(i[1],block), blocks)
    iblock_max = findfirst(block->in(i[end],block), blocks)
    iblock_min:iblock_max
end

struct SourcePoint{T,N}
    c::Array{T,N}
    r::NTuple{N,Vector{Int}}
    index::Int
end

struct SourceBlock{T,N}
    block_range::NTuple{N,UnitRange{Int}}
    source_points::Vector{SourcePoint{T,N}}
end

Base.range(block::SourceBlock, i) = block.block_range[i]

function source_blocking_range(i)
    a = mapreduce(minimum, min, i)::Int
    b = mapreduce(maximum, max, i)::Int
    a:b
end

"""
    sourceblocks = source_blocking(nz, nx, nbz, nbx, iz, ix, c)

Blocks the axes with `nbz*nbx` blocks, and partitions `iz`, `ix` and `c` into
sourceblocks::Vector{SourceBlock}.

# Inputs
* nz,nx is the finite-difference model size
* iz,ix are the injection points
* c are the injection filters (e.g. Hick's coefficients)
* nbz,nbx are the number of blocks in each dimension

see also WaveFD.injectdata!
"""
function source_blocking(nz, nx, nbz, nbx, iz, ix, c::Vector{<:AbstractMatrix{T}}) where {T}
    izrng = source_blocking_range(iz)
    ixrng = source_blocking_range(ix)
    zblocks = sourceblock_range(nbz, izrng) # TODO change name
    xblocks = sourceblock_range(nbx, ixrng)

    nblocks = nbz*nbx

    source_blocks = [SourceBlock((zblock,xblock), SourcePoint{T,2}[]) for zblock in zblocks, xblock in xblocks]
    for kshot = 1:length(iz)
        source_point = SourcePoint(c[kshot], (iz[kshot],ix[kshot]), kshot)

        jz_blocks = get_blocks_for_source_point(iz[kshot], zblocks)
        jx_blocks = get_blocks_for_source_point(ix[kshot], xblocks)

        for jx_block in jx_blocks, jz_block in jz_blocks
            push!(source_blocks[jz_block,jx_block].source_points, source_point)
        end
    end
    source_blocks
end

function source_blocking(nz, ny, nx, nbz, nby, nbx, iz, iy, ix, c::Vector{<:AbstractArray{T,3}}) where {T}
    izrng = source_blocking_range(iz)
    iyrng = source_blocking_range(iy)
    ixrng = source_blocking_range(ix)
    zblocks = sourceblock_range(nbz, izrng)
    yblocks = sourceblock_range(nby, iyrng)
    xblocks = sourceblock_range(nbx, ixrng)

    nblocks = nbz*nbx

    source_blocks = [SourceBlock((zblock,yblock,xblock), SourcePoint{T,3}[]) for zblock in zblocks, yblock in yblocks, xblock in xblocks]
    for kshot = 1:length(iz)
        source_point = SourcePoint(c[kshot], (iz[kshot],iy[kshot],ix[kshot]), kshot)

        jz_blocks = get_blocks_for_source_point(iz[kshot], zblocks)
        jy_blocks = get_blocks_for_source_point(iy[kshot], yblocks)
        jx_blocks = get_blocks_for_source_point(ix[kshot], xblocks)

        for jx_block in jx_blocks, jy_block in jy_blocks, jz_block in jz_blocks
            push!(source_blocks[jz_block,jy_block,jx_block].source_points, source_point)
        end
    end
    source_blocks
end

"""
  injectdata!(field, blks, data, it)

Inject data from data[it] into field. `blks` is computed using `WaveFD.source_blocking`.

The general work-flow is:
```
nz,nx=size(field)
nthreads=20
blks = WaveFD.source_blocking(nz,nx,nbz,nbx,iz,ix,c)
for it = 1:ntrec # time loop
    WaveFD.injectdata!(field, blks, data, it)
    ...
end
```
"""
function injectdata!(field::AbstractArray, blocks::AbstractArray{<:SourceBlock}, data::AbstractArray, it::Int)
    @sync for block in blocks
        Threads.@spawn injectdata_block!(field, block, data, it)
    end
end

function injectdata_block!(field::AbstractArray{T,2}, block::SourceBlock, data::AbstractArray{T}, it::Int) where {T}
    @inbounds for source_point in block.source_points
        iz,ix,i,c = source_point.r[1],source_point.r[2],source_point.index,source_point.c
        for jz = 1:length(iz)
            kz = iz[jz]
            kz ∈ range(block, 1) || continue
            for jx = 1:length(ix)
                kx = ix[jx]
                kx ∈ range(block, 2) || continue

                field[kz,kx] += c[jz,jx]*data[it,i]
            end
        end
    end
end

function injectdata_block!(field::AbstractArray{T,3}, block::SourceBlock, data::AbstractArray{T}, it::Int) where {T}
    @inbounds for source_point in block.source_points
        iz,iy,ix,i,c = source_point.r[1],source_point.r[2],source_point.r[3],source_point.index,source_point.c
        for jz = 1:length(iz)
            kz = iz[jz]
            kz ∈ range(block, 1) || continue
            for jy = 1:length(iy)
                ky = iy[jy]
                ky ∈ range(block, 2) || continue
                for jx = 1:length(ix)
                    kx = ix[jx]
                    kx ∈ range(block, 3) || continue

                    field[kz,ky,kx] += c[jz,jy,jx]*data[it,i]
                end
            end
        end
    end
end

"""
  injectdata!(field, data, it, iz, ix, c[, bz=10])
"""
function injectdata!(field::AbstractArray{T,2}, data::AbstractArray{T,2}, it::Integer, iz::AbstractVector{Vector{Int64}}, ix::AbstractVector{Vector{Int64}}, c::AbstractVector{Array{C,2}}, nbz::Integer=10, nbx::Integer=10) where {T,C}
    blocks = WaveFD.source_blocking(size(field)..., nbz, nbx, iz, ix, c)
    injectdata!(field, blocks, data, it)
end

"""
  injectdata!(field, data, it, iz, iy, ix, c[,bz=10])
"""
function injectdata!(field::AbstractArray{T,3}, data::AbstractArray{T,2}, it::Integer, iz::AbstractVector{Vector{Int64}}, iy::AbstractVector{Vector{Int64}}, ix::AbstractVector{Vector{Int64}}, c::AbstractVector{Array{C,3}}, nbz::Integer=1, nby::Integer=10, nbx::Integer=10) where {T,C}
    blocks = WaveFD.source_blocking(size(field)..., nbz, nby, nbx, iz, iy, ix, c)
    injectdata!(field, blocks, data, it)
end

function extractdata!(data::AbstractArray{C,2}, field::AbstractArray{C,2}, it::Integer, iz::AbstractArray{Array{Int64,1},1}, ix::AbstractArray{Array{Int64,1},1}, c::AbstractArray{Array{T,2},1}) where {T,C}
    for i = 1:length(iz)
        nc_z, nc_x = size(c[i])
        for ihx=1:nc_x
            @simd ivdep for ihz=1:nc_z
                @inbounds begin
                    data[it,i] += c[i][ihz,ihx]*field[iz[i][ihz],ix[i][ihx]]
                end
            end
        end
    end
end

function extractdata!(data::AbstractArray{C,2}, field::AbstractArray{C,3}, it::Integer, iz::AbstractArray{Array{Int64,1},1}, iy::AbstractArray{Array{Int64,1},1}, ix::AbstractArray{Array{Int64,1},1}, c::AbstractArray{Array{T,3},1}) where {T,C}
    Threads.@threads for i = 1:length(iz)
        nc_z, nc_y, nc_x = size(c[i])
        for ihx=1:nc_x, ihy=1:nc_y, ihz=1:nc_z
            @inbounds begin
                data[it,i] += c[i][ihz,ihy,ihx]*field[iz[i][ihz],iy[i][ihy],ix[i][ihx]]
            end
        end
    end
end

#
# injection/extraction methods specialized to on-grid source/receivers:
#
function allongrid(dz::T, dx::T, z0::Float64, x0::Float64, z::Array{Float64,1}, x::Array{Float64,1}; rtol=1e-6) where T
    for i = 1:length(z)
        if isapprox(z0 + round(Int, (z[i] - z0)/dz)*dz, z[i], rtol=rtol) == false
            return false
        end
        if isapprox(x0 + round(Int, (x[i] - x0)/dx)*dx, x[i], rtol=rtol) == false
            return false
        end
    end
    true
end

function allongrid(dz::T, dy::T, dx::T, z0::Float64, y0::Float64, x0::Float64, z::Array{Float64,1}, y::Array{Float64,1}, x::Array{Float64,1}; rtol=1e-6) where T
    for i = 1:length(z)
        if isapprox(z0 + round(Int, (z[i] - z0)/dz)*dz, z[i], rtol=rtol) == false
            return false
        end
        if isapprox(y0 + round(Int, (y[i] - y0)/dy)*dy, y[i], rtol=rtol) == false
            return false
        end
        if isapprox(x0 + round(Int, (x[i] - x0)/dx)*dx, x[i], rtol=rtol) == false
            return false
        end
    end
    true
end

function injectdata!(field::Array{T,2}, data::Array{T,2}, it::Integer, iz::Array{Int64,1}, ix::Array{Int64,1}, c::Array{T,1}, language::Language=LangC(), nthreads::Int=Sys.CPU_THREADS) where T
    if isa(language, LangC)
        if T == Float32
            ccall(
                (:injectdata_2d_ongrid_float, libspacetime),
                Cvoid,
                (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clong}, Ptr{Clong}, Csize_t, Csize_t,      Csize_t,   Csize_t,       Csize_t,            Csize_t),
                 field,       data,        c,           iz,         ix,         it,      size(data,1), length(c), size(field,1), size(field,2),      nthreads
            )
        else
            ccall(
                (:injectdata_2d_ongrid_double, libspacetime),
                Cvoid,
                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong}, Csize_t, Csize_t,      Csize_t,   Csize_t,       Csize_t,       Csize_t),
                 field,       data,          c,            iz,         ix,         it,      size(data,1), length(c), size(field,1), size(field,2), nthreads
            )
        end
        return
    end

    @inbounds for i = 1:length(iz)
        field[iz[i],ix[i]] += c[i]*data[it,i]
    end
end

function injectdata!(field::Array{T,3}, data::Array{T,2}, it::Integer, iz::Array{Int64,1}, iy::Array{Int64,1}, ix::Array{Int64,1}, c::Array{T,1}, language::Language=LangC(), nthreads::Int=Sys.CPU_THREADS) where T
    if isa(language, LangC)
        if T == Float32
            ccall(
                (:injectdata_3d_ongrid_float, libspacetime),
                Cvoid,
                (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Csize_t, Csize_t,      Csize_t,   Csize_t,       Csize_t,       Csize_t,       Csize_t),
                 field,       data,        c,           iz,         iy,         ix,         it,      size(data,1), length(c), size(field,1), size(field,2), size(field,3), nthreads
            )
        else
            ccall(
                (:injectdata_3d_ongrid_double, libspacetime),
                Cvoid,
                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Csize_t, Csize_t,      Csize_t,   Csize_t,       Csize_t,       Csize_t,       Csize_t),
                 field,       data,          c,            iz,         iy,         ix,         it,      size(data,1), length(c), size(field,1), size(field,2), size(field,3), nthreads
            )
        end
        return
    end

    @inbounds for i = 1:length(iz)
        field[iz[i],iy[i],ix[i]] += c[i]*data[it,i]
    end
end

function extractdata!(data::Array{T,2}, field::Array{T,2}, it::Integer, iz::Array{Int64,1}, ix::Array{Int64,1}, c::Array{T,1}, language::Language=LangC(), nthreads::Int=Sys.CPU_THREADS) where T
    if isa(language, LangC)
        if T == Float32
            ccall(
                (:extractdata_2d_ongrid_float, libspacetime),
                Cvoid,
                (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clong}, Ptr{Clong}, Clong, Clong,        Clong,     Clong,         Clong,         Clong),
                 data,        field,       c,           iz,         ix,         it,    size(data,1), length(c), size(field,1), size(field,2), nthreads
            )
        else
            ccall(
                (:extractdata_2d_ongrid_double, libspacetime),
                Cvoid,
                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong}, Clong, Clong,        Clong,     Clong,         Clong,         Clong),
                 data,         field,        c,            iz,         ix,         it,    size(data,1), length(c), size(field,1), size(field,2), nthreads
            )
        end
        return
    end
    @inbounds for i = 1:length(iz)
        data[it,i] += c[i]*field[iz[i],ix[i]]
    end
end

function extractdata!(data::Array{T,2}, field::Array{T,3}, it::Integer, iz::Array{Int64,1}, iy::Array{Int64,1}, ix::Array{Int64,1}, c::Array{T,1}, language::Language=LangC(), nthreads::Int=Sys.CPU_THREADS) where T
    if isa(language, LangC)
        if T == Float32
            ccall(
                (:extractdata_3d_ongrid_float, libspacetime),
                Cvoid,
                (Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Clong, Clong,        Clong,     Clong,         Clong,         Clong,         Clong),
                 data,        field,       c,           iz,         iy,         ix,         it,    size(data,1), length(c), size(field,1), size(field,2), size(field,3), nthreads
            )
        else
            ccall(
                (:extractdata_3d_ongrid_double, libspacetime),
                Cvoid,
                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Clong, Clong,        Clong,     Clong,         Clong,         Clong,         Clong),
                 data,         field,        c,            iz,         iy,         ix,         it,    size(data,1), length(c), size(field,1), size(field,2), size(field,3), nthreads
            )
        end
        return
    end
    @inbounds for i = 1:length(iz)
        data[it,i] += c[i]*field[iz[i],iy[i],ix[i]]
    end
end

#
# taper methods
#
function taper_region(nz, nx, iz, ix, iradius, iflat; interior=false)
    region = falses(nz, nx)
    addz = interior ? iflat : iflat+iradius
    addx = interior ? iflat : iflat+iradius
    for is = 1:length(ix)
        ilbz,iubz = clamp(minimum(iz[is])-addz,1,nz),clamp(maximum(iz[is])+addz,1,nz)
        ilbx,iubx = clamp(minimum(ix[is])-addx,1,nx),clamp(maximum(ix[is])+addx,1,nx)
        region[ilbz:iubz,ilbx:iubx] = true
    end

    region
end
function taper_region(nz, ny, nx, iz, iy, ix, iradius, iflat; interior=false)
    region = falses(nz, ny, nx)
    addz = interior ? iflat : iflat+iradius
    addy = interior ? iflat : iflat+iradius
    addx = interior ? iflat : iflat+iradius
    @inbounds for is = 1:length(ix)
        ilbz,iubz = clamp(minimum(iz[is])-addz,1,nz),clamp(maximum(iz[is])+addz,1,nz)
        ilby,iuby = clamp(minimum(iy[is])-addy,1,ny),clamp(maximum(iy[is])+addy,1,ny)
        ilbx,iubx = clamp(minimum(ix[is])-addx,1,nx),clamp(maximum(ix[is])+addx,1,nx)
        region[ilbz:iubz,ilby:iuby,ilbx:iubx] = true
    end

    region
end

function taper_region(taper_int::AbstractArray{Bool,2}, taper_ext)
    region_tpr = SVector{2,Float32}[]
    nz, nx = size(taper_int)
    @inbounds for jx = 1:nx, jz = 1:nz
        if !taper_int[jz,jx] && taper_ext[jz,jx]
            push!(region_tpr, SVector{2,Float32}(Float32(jz),Float32(jx)))
        end
    end
    region_tpr
end
function taper_region(taper_int::AbstractArray{Bool,3}, taper_ext)
    region_tpr = SVector{3,Float32}[]
    nz, ny, nx = size(taper_int)
    @inbounds for jx = 1:nx, jy = 1:ny, jz = 1:nz
        if !taper_int[jz,jy,jx] && taper_ext[jz,jy,jx]
            push!(region_tpr, SVector{3,Float32}(Float32(jz),Float32(jy),Float32(jx)))
        end
    end
    region_tpr
end

function taper_edge(region::AbstractArray{Bool,2})
    nz,nx = size(region)
    edge = SVector{2,Float32}[]
    @inbounds for ix = 2:nx-1, iz = 2:nz-1
        if region[iz,ix]
            if !region[iz+1,ix]
                push!(edge, SVector{2,Float32}(iz, ix))
                continue
            end
            if !region[iz-1,ix]
                push!(edge, SVector{2,Float32}(iz, ix))
                continue
            end
            if !region[iz,ix+1]
                push!(edge, SVector{2,Float32}(iz, ix))
                continue
            end
            if !region[iz,ix-1]
                push!(edge, SVector{2,Float32}(iz, ix))
            end
        end
    end
    KDTree(edge)
end
function taper_edge(region::AbstractArray{Bool,3})
    nz,ny,nx = size(region)
    edge = SVector{3,Float32}[]
    @inbounds for ix = 2:nx-1, iy = 2:ny-1, iz = 2:nz-1
        if region[iz,iy,ix]
            if !region[iz+1,iy,ix]
                push!(edge, SVector{3,Float32}(iz, iy, ix))
                continue
            end
            if !region[iz-1,iy,ix]
                push!(edge, SVector{3,Float32}(iz, iy, ix))
                continue
            end
            if !region[iz,iy+1,ix]
                push!(edge, SVector{3,Float32}(iz, iy, ix))
                continue
            end
            if !region[iz,iy-1,ix]
                push!(edge, SVector{3,Float32}(iz, iy, ix))
                continue
            end
            if !region[iz,iy,ix+1]
                push!(edge, SVector{3,Float32}(iz, iy, ix))
                continue
            end
            if !region[iz,iy,ix-1]
                push!(edge, SVector{3,Float32}(iz, iy, ix))
            end
        end
    end
    KDTree(edge)
end

function taper_build(t::Type{T}, nz, nx, region_tpr, region_int_edge, region_ext_edge) where T
    tpr = zeros(T, nz, nx)
    dint = knn(region_int_edge, region_tpr, 1)[2]
    dext = knn(region_ext_edge, region_tpr, 1)[2]
    for (i,tpr_pt) in enumerate(region_tpr)
        jz = Int(tpr_pt[1])
        jx = Int(tpr_pt[2])
        tpr[jz,jx] = dint[i][1]/(dint[i][1]+dext[i][1])
    end
    tpr
end
function taper_build(t::Type{T}, nz, ny, nx, region_tpr, region_int_edge, region_ext_edge) where T
    tpr = zeros(T, nz, ny, nx)
    dint = knn(region_int_edge, region_tpr, 1)[2]
    dext = knn(region_ext_edge, region_tpr, 1)[2]
    for (i,tpr_pt) in enumerate(region_tpr)
        jz = Int(tpr_pt[1])
        jy = Int(tpr_pt[2])
        jx = Int(tpr_pt[3])
        tpr[jz,jy,jx] = dint[i][1]/(dint[i][1]+dext[i][1])
    end
    tpr
end

function taper_etatozero(η, iz, ix, iradius, iflat)
    T = eltype(η)
    nz,nx = size(η)

    region_int = taper_region(nz, nx, iz, ix, iradius, iflat, interior=true)
    region_ext = taper_region(nz, nx, iz, ix, iradius, iflat, interior=false)

    region_int_edge = taper_edge(region_int)
    region_ext_edge = taper_edge(region_ext)
    region_tpr = taper_region(region_int, region_ext)

    tpr = taper_build(T, nz, nx, region_tpr, region_int_edge, region_ext_edge)

    for jx = 1:nx, jz = 1:nz
        if region_int[jz,jx] == 1
            η[jz,jx] = 0.0
            continue
        end
        if region_ext[jz,jx] == 1
            η[jz,jx] *= tpr[jz,jx]
        end
    end
    nothing
end
function taper_etatozero(η, iz, iy, ix, iradius, iflat)
    T = eltype(η)
    nz,ny,nx = size(η)

    region_int = taper_region(nz, ny, nx, iz, iy, ix, iradius, iflat, interior=true)
    region_ext = taper_region(nz, ny, nx, iz, iy, ix, iradius, iflat, interior=false)

    region_int_edge = taper_edge(region_int)
    region_ext_edge = taper_edge(region_ext)
    region_tpr = taper_region(region_int, region_ext)

    tpr = taper_build(T, nz, ny, nx, region_tpr, region_int_edge, region_ext_edge)

    for jx = 1:nx, jy = 1:ny, jz = 1:nz
        if region_int[jz,jy,jx] == 1
            η[jz,jy,jx] = 0.0
            continue
        end
        if region_ext[jz,jy,jx] == 1
            η[jz,jy,jx] *= tpr[jz,jy,jx]
        end
    end
    nothing
end

function taper_epsilon2delta(δ, ϵ, iz, ix, iradius, iflat)
    T = eltype(δ)
    nz,nx = size(δ)

    region_int = taper_region(nz, nx, iz, ix, iradius, iflat, interior=true)
    region_ext = taper_region(nz, nx, iz, ix, iradius, iflat, interior=false)

    region_int_edge = taper_edge(region_int)
    region_ext_edge = taper_edge(region_ext)
    region_tpr = taper_region(region_int, region_ext)

    tpr = taper_build(T, nz, nx, region_tpr, region_int_edge, region_ext_edge)

    for jx = 1:nx, jz = 1:nz
        if region_int[jz,jx]
            δ[jz,jx] = ϵ[jz,jx]
            continue
        end
        if region_ext[jz,jx]
            δ[jz,jx] = (1-tpr[jz,jx])*ϵ[jz,jx] + tpr[jz,jx]*δ[jz,jx]
        end
    end
    nothing
end
function taper_epsilon2delta(δ, ϵ, iz, iy, ix, iradius, iflat)
    T = eltype(δ)
    nz,ny,nx = size(δ)

    region_int = taper_region(nz, ny, nx, iz, iy, ix, iradius, iflat, interior=true)
    region_ext = taper_region(nz, ny, nx, iz, iy, ix, iradius, iflat, interior=false)

    region_int_edge = taper_edge(region_int)
    region_ext_edge = taper_edge(region_ext)
    region_tpr = taper_region(region_int, region_ext)

    tpr = taper_build(T, nz, ny, nx, region_tpr, region_int_edge, region_ext_edge)

    for jx = 1:nx, jy = 1:ny, jz = 1:nz
        if region_int[jz,jy,jx]
            δ[jz,jy,jx] = ϵ[jz,jy,jx]
            continue
        end
        if region_ext[jz,jy,jx]
            δ[jz,jy,jx] = (1-tpr[jz,jy,jx])*ϵ[jz,jy,jx] + tpr[jz,jy,jx]*δ[jz,jy,jx]
        end
    end
    nothing
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
using Random

@testset "Random Number Generation" begin
    model = ABM(Agent2)
    @test model.rng == Random.default_rng()
    rng = StableRNG(42)
    rng0 = StableRNG(42)

    model = ABM(Agent1, GridSpace((3,3)); rng)
    agent = Agent1(1, (1,1))
    add_agent_pos!(agent, model)
    agent = Agent1(2, (1,1))
    add_agent_single!(agent, model)
    # Test that model rng pool was used
    @test model.rng ≠ rng0
    @test agent.pos == (2,1)
    seed!(model, 42)
    @test model.rng == rng0

    model = ABM(Agent2; rng = RandomDevice())
    @test_throws MethodError seed!(model, 64)
end

@testset "sample!" begin
    rng = StableRNG(42)
    model = ABM(Agent2)
    for i in 1:20
        add_agent!(model, rand(model.rng) / rand(model.rng))
    end
    allweights = [i.weight for i in values(model.agents)]
    mean_weights = sum(allweights) / length(allweights)
    sample!(model, 12, :weight)
    @test nagents(model) == 12
    allweights = [i.weight for i in values(model.agents)]
    mean_weights_new = sum(allweights) / length(allweights)
    @test mean_weights_new > mean_weights
    sample!(model, 40, :weight)
    @test nagents(model) == 40
    allweights = [i.weight for i in values(model.agents)]
    mean_weights_new = sum(allweights) / length(allweights)
    @test mean_weights_new > mean_weights

    Random.seed!(6459)
    model2 = ABM(Agent3, GridSpace((10, 10)))
    for i in 1:20
        add_agent_single!(Agent3(i, (1, 1), rand(model2.rng) / rand(model2.rng)), model2)
    end
    @test sample!(model2, 10) === nothing
    @test sample!(model2, 10, :weight) === nothing
    allweights = [i.weight for i in values(model2.agents)]
    mean_weights = sum(allweights) / length(allweights)
    sample!(model2, 12, :weight)
    @test nagents(model2) == 12
    allweights = [i.weight for i in values(model2.agents)]
    mean_weights_new = sum(allweights) / length(allweights)
    @test mean_weights_new > mean_weights

    sample!(model2, 40, :weight)
    @test nagents(model2) == 40

    Random.seed!(6459)
    #Guarantee all starting weights are unique
    model3 = ABM(Agent2)
    while true
        for i in 1:20
            add_agent!(model3, rand(model3.rng) / rand(model3.rng))
        end
        allweights = [i.weight for i in values(model3.agents)]
        allunique(allweights) && break
    end
    # Cannot draw 50 samples out of a pool of 20 without replacement
    @test_throws ErrorException sample!(model3, 50, :weight; replace = false)
    sample!(model3, 15, :weight; replace = false)
    allweights = [i.weight for i in values(model3.agents)]
    @test allunique(allweights)
    model3 = ABM(Agent2)
    while true
        for i in 1:20
            add_agent!(model3, rand(model3.rng) / rand(model3.rng))
        end
        allweights = [i.weight for i in values(model3.agents)]
        allunique(allweights) && break
    end
    sample!(model3, 100, :weight; replace = true)
    allweights = [i.weight for i in values(model3.agents)]
    @test !allunique(allweights)
end

@testset "random agent" begin
    space = GridSpace((10, 10))
    model = ABM(Union{Daisy,Land}, space; warn = false)
    fill_space!(Daisy, model, "black")
    add_agent!(Land(999, (1, 1), 999), model)

    a = random_agent(model)
    @test typeof(a) <: Union{Daisy,Land}

    c1(a) = a isa Land
    a = random_agent(model, c1)
    @test a.id == 999

    c2(a) = a isa Float64
    a = random_agent(model, c2)
    @test isnothing(a)
end