text:
function wall_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, idx, configData, phi_i, phi_k)

    power::Float64 = configData["core"]["power"]::Float64
    limiter_flag::Float64 = configData["core"]["limiter_flag"]::Float64

    sum_delx_sqr = zero(Float64)
    sum_dely_sqr = zero(Float64)
    sum_delx_dely = zero(Float64)

    sum_delx_delf = zeros(Float64, 4)
    sum_dely_delf = zeros(Float64, 4)

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx::Float64 = ny
    ty::Float64 = -nx

    G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
    result = zeros(Float64,4)
    # phi_i = zeros(Float64,4)
    # phi_k = zeros(Float64,4)

    for itm in loc_globaldata[idx].xpos_conn

        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist^power

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr = sum_delx_sqr + dels*dels_weights
        sum_dely_sqr = sum_dely_sqr + deln*deln_weights
        sum_delx_dely = sum_delx_dely + dels*deln_weights

        qtilde_i = @. loc_globaldata[idx].q - 0.5*(delx * loc_globaldata[idx].dq[1] + dely * loc_globaldata[idx].dq[2])
        qtilde_k = @. globaldata_itm.q - 0.5*(delx * globaldata_itm.dq[1] + dely * globaldata_itm.dq[2])

        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), qtilde_k)
        # end

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata, idx, configData, phi_i)
            venkat_limiter(qtilde_k, loc_globaldata, idx, configData, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])
        end

        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), phi_i)
        #     println(IOContext(stdout, :compact => false), phi_k)
        #     println(IOContext(stdout, :compact => false), qtilde_k)
        # end

                #if limiter_flag == 2
        #    maxi = max_q_values(globaldata, idx)
        #    mini = min_q_values(globaldata, idx)
#
        #    for i in 1:4
        #        if qtilde_i[i] > maxi[i]
        #            qtilde_i[i] = maxi[i]
        #        end
        #        if qtilde_i[i] < mini[i]
        #            qtilde_i[i] = mini[i]
        #        end
        #        if qtilde_k[i] > maxi[i]
        #            qtilde_k[i] = maxi[i]
        #        end
        #        if qtilde_k[i] < mini[i]
        #            qtilde_k[i] = mini[i]
        #        end
        #    end
        #end
        # if idx == 76
        #     println("qtile_i ", qtilde_i)
        #     println("qtile_k ", qtilde_k)
        # end
        if idx == 100
            flag = 0
        else
            flag = 1
        end

        # if idx == 3
        #     println(IOContext(stdout, :compact => false), result)
        # end
        qtilde_to_primitive(result, qtilde_i, configData)
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), result)
        # end
        flux_quad_GxII(G_i, nx, ny, result[1], result[2], result[3], result[4], flag)
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), result)
        # end
        qtilde_to_primitive(result, qtilde_k, configData)
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), result)
        # end
        flux_quad_GxII(G_k, nx, ny, result[1], result[2], result[3], result[4], flag)
        # if idx == 100
        #     println("****")
        #     println(IOContext(stdout, :compact => false), phi_i)
        #     println(IOContext(stdout, :compact => false), phi_k)
        #     println(IOContext(stdout, :compact => false), qtilde_i)
        #     println(IOContext(stdout, :compact => false), qtilde_k)
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), result)
        #     println(IOContext(stdout, :compact => false), G_i)
        #     println(IOContext(stdout, :compact => false), G_k)
        #     # println(IOContext(stdout, :compact => false), dels_weights)
        #     # println(IOContext(stdout, :compact => false), deln_weights)
        #     println("****")
        # end
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), result)
        #     println(IOContext(stdout, :compact => false), G_i)
        #     println(IOContext(stdout, :compact => false), G_k)
        #     # println(IOContext(stdout, :compact => false), dels_weights)
        #     # println(IOContext(stdout, :compact => false), deln_weights)
        # end
        for i in 1:4
            sum_delx_delf[i] += (G_k[i] - G_i[i]) * dels_weights
            sum_dely_delf[i] += (G_k[i] - G_i[i]) * deln_weights
        end
    end

    det = @. sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det
    # if idx == 3
    #     println(IOContext(stdout, :compact => false), "===Gx===")
    #     # println(IOContext(stdout, :compact => false), sum_delx_sqr)
    #     # println(IOContext(stdout, :compact => false), sum_dely_sqr)
    #     # println(IOContext(stdout, :compact => false), sum_delx_dely)
    #     # println(IOContext(stdout, :compact => false), det)
    #     # println(IOContext(stdout, :compact => false), one_by_det)
    #     println(IOContext(stdout, :compact => false), sum_delx_delf)
    #     println(IOContext(stdout, :compact => false), sum_dely_delf)
    #     println(IOContext(stdout, :compact => false), G)
    #     println()
    # end
    return @. (sum_delx_delf*sum_dely_sqr - sum_dely_delf*sum_delx_dely)*one_by_det
end

function wall_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, configData, phi_i, phi_k)

    power::Float64 = configData["core"]["power"]::Float64
    limiter_flag::Float64 = configData["core"]["limiter_flag"]::Float64

    sum_delx_sqr = zero(Float64)
    sum_dely_sqr = zero(Float64)
    sum_delx_dely = zero(Float64)

    sum_delx_delf = zeros(Float64, 4)
    sum_dely_delf = zeros(Float64, 4)

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
    result = zeros(Float64,4)
    # phi_i = zeros(Float64,4)
    # phi_k = zeros(Float64,4)

    for itm in loc_globaldata[idx].xneg_conn
        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist^power

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr = sum_delx_sqr + dels*dels_weights
        sum_dely_sqr = sum_dely_sqr + deln*deln_weights

        sum_delx_dely = sum_delx_dely + dels*deln_weights

        qtilde_i = @. loc_globaldata[idx].q - 0.5*(delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
        qtilde_k = @. globaldata_itm.q - 0.5*(delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata, idx, configData, phi_i)
            venkat_limiter(qtilde_k, loc_globaldata, idx, configData, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])
        end

        #if limiter_flag == 2
        #    maxi = max_q_values(globaldata, idx)
        #    mini = min_q_values(globaldata, idx)
#
        #    for i in 1:4
        #        if qtilde_i[i] > maxi[i]
        #            qtilde_i[i] = maxi[i]
        #        end
        #        if qtilde_i[i] < mini[i]
        #            qtilde_i[i] = mini[i]
        #        end
        #        if qtilde_k[i] > maxi[i]
        #            qtilde_k[i] = maxi[i]
        #        end
        #        if qtilde_k[i] < mini[i]
        #            qtilde_k[i] = mini[i]
        #        end
        #    end
        #end
        qtilde_to_primitive(result, qtilde_i, configData)
        flux_quad_GxI(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, configData)
        flux_quad_GxI(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            sum_delx_delf[i] += (G_k[i] - G_i[i]) * dels_weights
            sum_dely_delf[i] += (G_k[i] - G_i[i]) * deln_weights
        end

    end
    det = @. sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det
    return @. (sum_delx_delf*sum_dely_sqr - sum_dely_delf*sum_delx_dely)*one_by_det
end

function wall_dGy_neg(loc_globaldata, loc_ghost_holder, dist_length, idx, configData, phi_i, phi_k)

    power::Float64 = configData["core"]["power"]::Float64
    limiter_flag::Float64 = configData["core"]["limiter_flag"]::Float64

    sum_delx_sqr = zero(Float64)
    sum_dely_sqr = zero(Float64)
    sum_delx_dely = zero(Float64)

    sum_delx_delf = zeros(Float64, 4)
    sum_dely_delf = zeros(Float64, 4)

    x_i = loc_globaldata[idx].x
    y_i = loc_globaldata[idx].y

    nx = loc_globaldata[idx].nx
    ny = loc_globaldata[idx].ny

    tx = ny
    ty = -nx

    G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
    result = zeros(Float64,4)
    # phi_i = zeros(Float64,4)
    # phi_k = zeros(Float64,4)

    for itm in loc_globaldata[idx].yneg_conn

        if itm <= dist_length
            globaldata_itm = loc_globaldata[itm]
        else
            globaldata_itm = loc_ghost_holder[1][itm]
        end
        x_k = globaldata_itm.x
        y_k = globaldata_itm.y

        delx = x_k - x_i
        dely = y_k - y_i

        dels = delx*tx + dely*ty
        deln = delx*nx + dely*ny

        dist = hypot(dels, deln)
        weights = dist^power

        dels_weights = dels*weights
        deln_weights = deln*weights

        sum_delx_sqr = sum_delx_sqr + dels*dels_weights
        sum_dely_sqr = sum_dely_sqr + deln*deln_weights

        sum_delx_dely = sum_delx_dely + dels*deln_weights

        qtilde_i = @. loc_globaldata[idx].q - 0.5*(delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
        qtilde_k = @. globaldata_itm.q - 0.5*(delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])

        if limiter_flag == 1
            venkat_limiter(qtilde_i, loc_globaldata, idx, configData, phi_i)
            venkat_limiter(qtilde_k, loc_globaldata, idx, configData, phi_k)
            @. qtilde_i = loc_globaldata[idx].q - 0.5 * phi_i * (delx*loc_globaldata[idx].dq[1] + dely*loc_globaldata[idx].dq[2])
            @. qtilde_k = globaldata_itm.q - 0.5 * phi_k * (delx*globaldata_itm.dq[1] + dely*globaldata_itm.dq[2])
        end
                #if limiter_flag == 2
        #    maxi = max_q_values(globaldata, idx)
        #    mini = min_q_values(globaldata, idx)
#
        #    for i in 1:4
        #        if qtilde_i[i] > maxi[i]
        #            qtilde_i[i] = maxi[i]
        #        end
        #        if qtilde_i[i] < mini[i]
        #            qtilde_i[i] = mini[i]
        #        end
        #        if qtilde_k[i] > maxi[i]
        #            qtilde_k[i] = maxi[i]
        #        end
        #        if qtilde_k[i] < mini[i]
        #            qtilde_k[i] = mini[i]
        #        end
        #    end
        #end
        qtilde_to_primitive(result, qtilde_i, configData)
        flux_Gyn(G_i, nx, ny, result[1], result[2], result[3], result[4])

        qtilde_to_primitive(result, qtilde_k, configData)
        flux_Gyn(G_k, nx, ny, result[1], result[2], result[3], result[4])

        for i in 1:4
            sum_delx_delf[i] += (G_k[i] - G_i[i]) * dels_weights
            sum_dely_delf[i] += (G_k[i] - G_i[i]) * deln_weights
        end
        # if idx == 3
        #     println(IOContext(stdout, :compact => false), itm)
        #     println(IOContext(stdout, :compact => false), result)
        #     println(IOContext(stdout, :compact => false), G_i)
        #     println(IOContext(stdout, :compact => false), G_k)
        # end
    end

    det = @. sum_delx_sqr*sum_dely_sqr - sum_delx_dely*sum_delx_dely
    one_by_det = 1.0 / det
    # if idx == 3
    #     println(IOContext(stdout, :compact => false), "===Gx===")
    #     println(IOContext(stdout, :compact => false), sum_delx_delf)
    #     println(IOContext(stdout, :compact => false), sum_dely_delf)
    #     # println(IOContext(stdout, :compact => false), G)
    #     # println()
    # end
    return @. (sum_dely_delf*sum_delx_sqr - sum_delx_delf*sum_delx_dely)*one_by_det
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
using Test
using Polynomials
using Jacobi

include("testjacobi.jl")
include("testgauss.jl")

ext: jl


text:
function wall_locations(grid::Grid{Rect})::Array{Int32,3}
    nrows, ncols = grid.nrows, grid.ncols
    walls = zeros(Int32, (nrows, ncols, 4))
    directions = [north, south, east, west]
    for row in 1:nrows
        for col in 1:ncols
            cell = grid[row, col]
            for i in 1:4
                d = directions[i]
                if !is_connected(cell, d(grid, cell))
                    walls[row, col, i] = 1
                end
            end
        end
    end
    return walls
end


text:
module Banners

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export banner

register_normal_element("q__banner", context = @__MODULE__)

function banner(content::String = "",
                args...;
                buttons::Vector{String} = String[],
                icon::Union{String,Nothing} = nothing,
                kwargs...)

  q__banner(args...; kw(kwargs)...) do
    string(
      (icon !== nothing ? Genie.Renderer.Html.template(()->icon, Symbol("v-slot:avatar")) : ""),
      content,
      (! isempty(buttons) ? Genie.Renderer.Html.template(()->join(buttons, "\n"), Symbol("v-slot:action")) : "")
    )
  end
end

end


text:
using BioDataExpression
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@test 1 == 1

dbfile = joinpath(pwd(),"test", "data","database.sqlite");

@test isfile(dbfile)

data = BioData( dbfile );

@test length( nonzeros(data.data) ) == 97784

Normalize(data, 100)

@test length( nonzeros(data.raw) ) == 97784

@test length( nonzeros(data.data) ) == 77660

@test corspearman( log(toPlot( data.data[1:100,1])+1),log10(toPlot(data.raw[1:100,1])+1) ) > 0.9

res = zersos( size(data.data)[2] )
for i=1:size(data.data)[2]
  res[i] = corspearman( log(toPlot( data.data[1:100,i])+1),log10(toPlot(data.raw[1:100,i])+1) ) > 0.8
end

@test length(find(res .== 0)) == 0

Zscore(data)

@test length( nonzeros(data.zscored) ) == 77660

@test corspearman( log(toPlot( data.data[1:100,1])+1),log10(toPlot(data.zscored[1:100,1])+1) ) > 0.9

res = zersos( size(data.data)[2] )
for i=1:size(data.data)[2]
  res[i] = corspearman( log(toPlot( data.data[1:100,i])+1),log10(toPlot(data.zscored[1:100,i])+1) ) > 0.8
end

@test length(find(res .== 0)) == 0


text:
@testset "910.smallest-range-ii.jl" begin
    @test smallest_rangeII([1], 0) == 0
    @test smallest_rangeII([0, 10], 2) == 6
    @test smallest_rangeII([1, 3, 6], 3) == 3
    @test smallest_rangeII([1, 1, 3, 4, 10, 12, 13, 14, 26, 27, 28, 42, 42, 48, 50, 50, 53, 58, 60, 61, 63, 66, 72, 73, 75, 75, 80, 86, 90, 93], 48) == 82
end


text:
# tools for model building

import Hypatia
import Hypatia.Cones: smat_to_svec!, svec_to_smat!, svec_side, svec_length

const rt2 = sqrt(2)

function svec(mat::AbstractMatrix{T}) where {T}
    vec = zeros(T, svec_length(size(mat, 1)))
    return smat_to_svec!(vec, mat, rt2)
end

function svec(mat::AbstractMatrix{Complex{T}}) where {T}
    vec = zeros(T, svec_length(ComplexF64, size(mat, 1)))
    return smat_to_svec!(vec, mat, rt2)
end

function smat(::Type{Float64}, vec::AbstractVector{T}) where {T}
    side = svec_side(length(vec))
    mat = zeros(T, side, side)
    return svec_to_smat!(mat, vec, rt2)
end

function smat(::Type{ComplexF64}, vec::AbstractVector{T}) where {T}
    side = svec_side(ComplexF64, length(vec))
    mat = zeros(Complex{T}, side, side)
    return svec_to_smat!(mat, vec, rt2)
end

function scale_svec(
    R::Type{<:Union{Float64, ComplexF64}},
    vec::AbstractVector{T},
    scal::Float64,
) where {T}
    incr = (R == Float64 ? 1 : 2)
    svec = copy(vec)
    k = 1
    for j in 1:svec_side(R, length(vec))
        for _ in 1:(incr * (j - 1))
            # scale off-diagonal
            svec[k] *= scal
            k += 1
        end
        k += 1
    end
    @assert k == length(vec) + 1
    return svec
end


text:
using OrdinaryDiffEq, StaticArrays, Test

function time_derivative(du,u,p,t)
  du[1] = -t
end
function time_derivative_static(u,p,t)
  SVector(-t)
end
function time_derivative_analytic(u0,p,t)
  u0 .- t.^2 ./ 2
end

const CACHE_TEST_ALGS = [Euler(),Midpoint(),RK4(),SSPRK22(),SSPRK33(),SSPRK53(),
  SSPRK63(),SSPRK73(),SSPRK83(),SSPRK43(),SSPRK432(),SSPRK932(),SSPRK54(),SSPRK104(),CarpenterKennedy2N54(),
  BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),Feagin14(),TanYam7(),
  Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9(),OwrenZen3(),OwrenZen4(),OwrenZen5()]

tspan = (0.0,1.0)

for (ff_time_derivative, u0) in
  ((ODEFunction(time_derivative, analytic=time_derivative_analytic), [1.0]),
   (ODEFunction(time_derivative_static, analytic=time_derivative_analytic), SVector(1.0)))

  @info "StaticArrays?: $(u0 isa StaticArray)"

  prob = ODEProblem(ff_time_derivative,u0,tspan)

  for _autodiff in (true, false)
    @info "autodiff=$(_autodiff)"

    @show Rosenbrock23
    sol = solve(prob,Rosenbrock32(autodiff=_autodiff),reltol=1e-9,abstol=1e-9)
    @test sol.errors[:final] < 1e-5*10^(!_autodiff)
    sol = solve(prob,Rosenbrock23(autodiff=_autodiff),reltol=1e-9,abstol=1e-9)
    @test sol.errors[:final] < 1e-10*10_000^(!_autodiff)

    @show Rodas4
    sol = solve(prob,Rodas4(autodiff=_autodiff),reltol=1e-9,abstol=1e-9)
    @test sol.errors[:final] < 1e-10*100_000^(!_autodiff)

    @show Rodas5
    sol = solve(prob,Rodas5(autodiff=_autodiff),reltol=1e-9,abstol=1e-9)
    @test sol.errors[:final] < 1e-10*1_000_000_000^(!_autodiff)

    @show Veldd4
    sol = solve(prob,Veldd4(autodiff=_autodiff),reltol=1e-9,abstol=1e-9)
    @test sol.errors[:final] < 1e-10*100_000^(!_autodiff)

    @show KenCarp3
    sol = solve(prob,KenCarp3(autodiff=_autodiff),reltol=1e-12,abstol=1e-12)
    @test length(sol) > 2
    @test sol.errors[:final] < 1e-10

    @show KenCarp4
    sol = solve(prob,KenCarp4(autodiff=_autodiff),reltol=1e-12,abstol=1e-12)
    @test length(sol) > 2
    @test sol.errors[:final] < 1e-10

    @show KenCarp47
    sol = solve(prob,KenCarp47(autodiff=_autodiff),reltol=1e-12,abstol=1e-12)
    @test length(sol) > 2
    @test sol.errors[:final] < 1e-10

    @show KenCarp5
    sol = solve(prob,KenCarp5(autodiff=_autodiff),reltol=1e-12,abstol=1e-12)
    @test length(sol) > 2
    @test sol.errors[:final] < 1e-10

    @show KenCarp58
    sol = solve(prob,KenCarp58(autodiff=_autodiff),reltol=1e-12,abstol=1e-12)
    @test length(sol) > 2
    @test sol.errors[:final] < 1e-10

    @show TRBDF2
    sol = solve(prob,TRBDF2(autodiff=_autodiff),reltol=1e-9,abstol=1e-9)
    @test sol.errors[:final] < 1e-10

    @show ImplicitEuler
    sol = solve(prob,ImplicitEuler(autodiff=_autodiff),dt=1/10)
    @test sol.errors[:final] < 1e-1

    @show Trapezoid
    sol = solve(prob,Trapezoid(autodiff=_autodiff),dt=1/10)
    @test sol.errors[:final] < 1e-12
  end

  @show Euler
  sol = solve(prob,Euler(),dt=1/100)
  @test sol.errors[:final] < 6e-3

  for alg in CACHE_TEST_ALGS
    @show alg
    sol = solve(prob,alg,dt=1/10)
    if !(typeof(alg) <: Euler)
      @test sol.errors[:final] < 4e-14
    end
  end
end


text:
function printsequence(io::IO, source::Enumerable)
    T = eltype(source)
    rows = Base.IteratorSize(source) == Base.HasLength() ? length(source) : "?"

    print(io, "$(rows)-element query result")

    max_element_to_show = 10

    i = 1
    foo = iterate(source)
    while foo!==nothing
        v, s = foo
        println(io)
        if i==max_element_to_show+1
            print(io, "... with ")
            if Base.IteratorSize(source)!=Base.HasLength()
                print(io, " more elements")
            else
                extra_rows = length(source) - max_element_to_show
                print(io, "$extra_rows more $(extra_rows==1 ? "element" : "elements")")
            end
            break
        else
            print(io, " ")
            show(IOContext(io, :compact => true), v)
        end
        i += 1

        foo = iterate(source, s)
    end
end

function Base.show(io::IO, source::Enumerable)
    if eltype(source) <: NamedTuple
        TableShowUtils.printtable(io, source, "query result")
    else
        printsequence(io, source)
    end
end

function Base.show(io::IO, ::MIME"text/html", source::Enumerable)
    if eltype(source) <: NamedTuple
        TableShowUtils.printHTMLtable(io, source)
    else
        error("Cannot write this Enumerable as text/html.")
    end
end

function Base.Multimedia.showable(::MIME"text/html", source::Enumerable)
    return eltype(source) <: NamedTuple
end

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::Enumerable)
    if eltype(source) <: NamedTuple
        TableShowUtils.printdataresource(io, source)
    else
        error("Cannot write this Enumerable as 'application/vnd.dataresource+json'.")
    end
end

function Base.Multimedia.showable(::MIME"application/vnd.dataresource+json", source::Enumerable)
    return eltype(source) <: NamedTuple
end


text:
include("operators.jl")

### Terms
# - A Term is something that, given a state, returns a named tuple (E, hams) with an energy
#   and a list of RFO (for each kpoint).
# - Each term must overload
#     `ene_ops(t, ψ, occ; kwargs...)` -> (E::Real, ops::Vector{RealFourierOperator}).
# - Note that terms are allowed to hold on to references to ψ (eg Fock term),
#   so ψ should not mutated after ene_ops

# The Hamiltonian is defined as half of the gradient of the energy
# with respect to the density matrix sum_n fn |ψn><ψn|.
# In particular, dE/dψn = 2 fn |Hψn> (plus weighting for kpoint sampling)
abstract type Term end

### Builders are objects X that store the term parameters, and produce a
# XTerm <: Term when instantiated with a `basis`

"""
A term with a constant zero energy.
"""
struct NoopTerm <: Term
    basis::PlaneWaveBasis
end
function ene_ops(term::NoopTerm, ψ, occ; kwargs...)
    (E=zero(eltype(term.basis)), ops=[NoopOperator(term.basis, kpoint)
                                      for kpoint in term.basis.kpoints])
end

include("Hamiltonian.jl")

# breaks_symmetries on a term builder answers true if this term breaks
# the symmetries of the lattice/atoms (in which case kpoint reduction
# is invalid)
breaks_symmetries(term_builder::Any) = false

include("kinetic.jl")

include("local.jl")
breaks_symmetries(term_builder::ExternalFromReal) = true
breaks_symmetries(term_builder::ExternalFromFourier) = true

include("nonlocal.jl")
include("hartree.jl")
include("power_nonlinearity.jl")
include("xc.jl")
include("ewald.jl")
include("psp_correction.jl")
include("entropy.jl")

include("magnetic.jl")
breaks_symmetries(term_builder::Magnetic) = true

include("anyonic.jl")
breaks_symmetries(term_builder::Anyonic) = true

# forces computes either nothing or an array forces[el][at][α]
compute_forces(term::Term, ψ, occ; kwargs...) = nothing  # by default, no force

"""
Compute the forces of an obtained SCF solution. Returns the forces wrt. the fractional
lattice vectors. To get cartesian forces use [`compute_forces_cart`](@ref).
Returns a list of lists of forces
`[[force for atom in positions] for (element, positions) in atoms]`
which has the same structure as the `atoms` object passed to the underlying [`Model`](@ref).
"""
@timing function compute_forces(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    if !any(iszero(kpt.coordinate) for kpt in basis.kpoints)
        @warn "Forces for shifted k-Grids not tested"
    end

    # TODO optimize allocs here
    T = eltype(basis)
    forces = [zeros(Vec3{T}, length(positions)) for (element, positions) in basis.model.atoms]
    for term in basis.terms
        f_term = compute_forces(term, ψ, occ; kwargs...)
        if !isnothing(f_term)
            forces += f_term
        end
    end
    forces
end

"""
Compute the cartesian forces of an obtained SCF solution in Hartree / Bohr.
Returns a list of lists of forces
`[[force for atom in positions] for (element, positions) in atoms]`
which has the same structure as the `atoms` object passed to the underlying [`Model`](@ref).
"""
function compute_forces_cart(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    lattice = basis.model.lattice
    forces = compute_forces(basis::PlaneWaveBasis, ψ, occ; kwargs...)
    [[lattice \ f for f in forces_for_element] for forces_for_element in forces]
end

function compute_forces(scfres)
    compute_forces(scfres.basis, scfres.ψ, scfres.occupation; ρ=scfres.ρ, ρspin=scfres.ρspin)
end
function compute_forces_cart(scfres)
    compute_forces_cart(scfres.basis, scfres.ψ, scfres.occupation;
                        ρ=scfres.ρ, ρspin=scfres.ρspin)
end

@doc raw"""
    compute_kernel(basis::PlaneWaveBasis; kwargs...)

Computes a matrix representation of the full response kernel
(derivative of potential with respect to density) in real space.
For non-spin-polarized calculations the matrix dimension is
`prod(basis.fft_size)` × `prod(basis.fft_size)` and
for collinear spin-polarized cases it is
`2prod(basis.fft_size)` × `2prod(basis.fft_size)`.
In this case the matrix has effectively 4 blocks, which are:
```math
\left(\begin{array}{cc}
    K_{α, \text{tot}} & K_{α, \text{spin}}\\
    K_{β, \text{tot}} & K_{β, \text{spin}}
\end{array}\right)
```
i.e. corresponding to a mapping
``(ρ_\text{tot}, ρ_\text{spin})^T = (ρ_α + ρ_β, ρ_α - ρ_β)^T ↦ (V_α, V_β)^T``.
"""
@timing function compute_kernel(basis::PlaneWaveBasis{T}; kwargs...) where {T}
    n_spin = basis.model.n_spin_components
    kernel = zeros(T, n_spin * prod(basis.fft_size), n_spin * prod(basis.fft_size))
    for term in basis.terms
        isnothing(term) && continue
        kernel .+= compute_kernel(term; kwargs...)
    end
    kernel
end
compute_kernel(::Term; kwargs...) = nothing  # By default no kernel

"""
    apply_kernel(basis::PlaneWaveBasis, dρ, dρspin=nothing; kwargs...)

Computes the potential response to a perturbation `(dρ, dρspin)` in real space.
Returns the array `[dV_α, dV_β]` for collinear spin-polarized
calculations, else the array [dV_{tot}].
"""
@timing function apply_kernel(basis::PlaneWaveBasis, dρ, dρspin=nothing;
                              RPA=false, kwargs...)
    n_spin = basis.model.n_spin_components
    (n_spin == 1) && @assert isnothing(dρspin)
    (n_spin == 2) && @assert !isnothing(dρspin)
    @assert 1 ≤ n_spin ≤ 2

    dV = [RealFourierArray(basis) for _ in 1:n_spin]
    for term in basis.terms
        # Skip XC term if RPA is selected
        RPA && term isa TermXc && continue

        dV_term = apply_kernel(term, dρ, dρspin; kwargs...)
        if !isnothing(dV_term)
            for σ in 1:n_spin
                dV[σ].real .+= dV_term[σ].real
            end
        end
    end
    dV
end
apply_kernel(::Term, dρ, dρspin; kwargs...) = nothing  # by default, no kernel


text:
using LinearSolve, LinearAlgebra, SparseArrays
using CUDA
using Test

n = 8
A = Matrix(I,n,n)
b = ones(n)
A1 = A/1; b1 = rand(n); x1 = zero(b)
A2 = A/2; b2 = rand(n); x2 = zero(b)

prob1 = LinearProblem(A1, b1; u0=x1)
prob2 = LinearProblem(A2, b2; u0=x2)

cache_kwargs = (;verbose=true, abstol=1e-8, reltol=1e-8, maxiter=30,)

function test_interface(alg, prob1, prob2)
    A1 = prob1.A; b1 = prob1.b; x1 = prob1.u0
    A2 = prob2.A; b2 = prob2.b; x2 = prob2.u0

    y = solve(prob1, alg; cache_kwargs...)
    @test A1 *  y  ≈ b1

    cache = SciMLBase.init(prob1,alg; cache_kwargs...) # initialize cache
    y = solve(cache)
    @test A1 *  y  ≈ b1

    cache = LinearSolve.set_A(cache,copy(A2))
    y = solve(cache)
    @test A2 *  y  ≈ b1

    cache = LinearSolve.set_b(cache,b2)
    y = solve(cache)
    @test A2 *  y  ≈ b2

    return
end

test_interface(GPUOffloadFactorization(), prob1, prob2)

A1 = prob1.A; b1 = prob1.b; x1 = prob1.u0
y = solve(prob1)
@test A1 *  y  ≈ b1


text:
"""
    advection_matrix(domain, n)

Assemble discrete advection matrix.
"""
function advection_matrix end

function advection_matrix(domain::ClosedIntervalDomain, n)
    Δx = 1 // n * (domain.right - domain.left)
    inds = [-1, 1]
    stencil = [-1 // 2, 1 // 2] / Δx
    diags = [i => fill(s, n + 1 - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    C = spdiagm(diags...)
    C[1, 1] = 2stencil[1]
    C[1, 2] = 2stencil[2]
    C[end, end-1] = 2stencil[1]
    C[end, end] = 2stencil[2]
    C
end

function advection_matrix(domain::PeriodicIntervalDomain, n)
    Δx = 1 // n * (domain.right - domain.left)
    inds = [-1, 1]
    stencil = [-1 // 2, 1 // 2] / Δx
    diags = [i => fill(s, n - abs(i)) for (i, s) ∈ zip(inds, stencil)]
    C = spdiagm(diags...)
    C[1, end] = stencil[1]
    C[end, 1] = stencil[2]
    C
end


text:
struct EnumerableOrderby{T,S,KS<:Function,TKS} <: Enumerable
    source::S
    keySelector::KS
    descending::Bool
end

Base.IteratorSize(::Type{EnumerableOrderby{T,S,KS,TKS}}) where {T,S,KS,TKS} = haslength(S)

Base.eltype(::Type{EnumerableOrderby{T,S,KS,TKS}}) where {T,S,KS,TKS} = T

Base.length(iter::EnumerableOrderby{T,S,KS,TKS}) where {T,S,KS,TKS} = length(iter.source)

function orderby(source::Enumerable, f::Function, f_expr::Expr)
    T = eltype(source)
    TKS = Base._return_type(f, Tuple{T,})

    KS = typeof(f)

    return EnumerableOrderby{T,typeof(source), KS,TKS}(source, f, false)
end

function orderby_descending(source::Enumerable, f::Function, f_expr::Expr)
    T = eltype(source)
    TKS = Base._return_type(f, Tuple{T,})

    KS = typeof(f)

    return EnumerableOrderby{T,typeof(source),KS,TKS}(source, f, true)
end

function Base.iterate(iter::EnumerableOrderby{T,S,KS,TKS}) where {T,S,KS,TKS}
    rows = (Base.IteratorSize(typeof(iter)) isa Base.HasLength || Base.IteratorSize(typeof(iter)) isa Base.HasShape) ? length(iter) : 0

    elements = Array{T}(undef, rows)

    if Base.IteratorSize(typeof(iter)) isa Base.HasLength || Base.IteratorSize(typeof(iter)) isa Base.HasShape
        for i in enumerate(iter.source)
            elements[i[1]] = i[2]
        end
    else
        for i in iter.source
            push!(elements, i)
        end
    end

    if length(elements)==0
        return nothing
    end

    sort!(elements, by=iter.keySelector, rev=iter.descending)

    return elements[1], (elements, 2)
end

function Base.iterate(iter::EnumerableOrderby{T,S,KS,TKS}, state) where {T,S,KS,TKS}
    if state[2]>length(state[1])
        return nothing
    else
        return state[1][state[2]], (state[1], state[2]+1)
    end
end

struct EnumerableThenBy{T,S,KS<:Function,TKS} <: Enumerable
    source::S
    keySelector::KS
    descending::Bool
end

Base.eltype(::Type{EnumerableThenBy{T,S,KS,TKS}}) where {T,S,KS,TKS} = T

Base.length(iter::EnumerableThenBy{T,S,KS,TKS}) where {T,S,KS,TKS} = length(iter.source)

function thenby(source::Enumerable, f::Function, f_expr::Expr)
    T = eltype(source)
    TKS = Base._return_type(f, Tuple{T,})
    KS = typeof(f)
    return EnumerableThenBy{T,typeof(source),KS,TKS}(source, f, false)
end

function thenby_descending(source::Enumerable, f::Function, f_expr::Expr)
    T = eltype(source)
    TKS = Base._return_type(f, Tuple{T,})
    KS = typeof(f)
    return EnumerableThenBy{T,typeof(source),KS,TKS}(source, f, true)
end

function Base.iterate(iter::EnumerableThenBy{T,S,KS,TKS}) where {T,S,KS,TKS}
    # Find start of ordering sequence
    source = iter.source
    keySelectors = [source.keySelector,iter.keySelector]
    directions = [source.descending, iter.descending]
    while !isa(source, EnumerableOrderby)
        source = source.source
        insert!(keySelectors,1,source.keySelector)
        insert!(directions,1,source.descending)
    end
    keySelector = element->[i(element) for i in keySelectors]

    lt = (t1,t2) -> begin
        n1, n2 = length(t1), length(t2)
        for i = 1:min(n1, n2)
            a, b = t1[i], t2[i]
            descending = directions[i]
            if !isequal(a, b)
                return descending ? !isless(a, b) : isless(a, b)
            end
        end
        return n1 < n2
    end

    rows = (Base.IteratorSize(typeof(iter)) isa Base.HasLength || Base.IteratorSize(typeof(iter)) isa Base.HasShape) ? length(iter) : 0

    elements = Array{T}(undef, rows)

    if (Base.IteratorSize(typeof(iter)) isa Base.HasLength || Base.IteratorSize(typeof(iter)) isa Base.HasShape)
        for i in enumerate(iter.source)
            elements[i[1]] = i[2]
        end
    else
        for i in iter.source
            push!(elements, i)
        end
    end

    if length(elements)==0
        return nothing
    end

    sort!(elements, by=keySelector, lt=lt)

    return elements[1], (elements, 2)
end

function Base.iterate(iter::EnumerableThenBy{T,S,KS,TKS}, state) where {T,S,KS,TKS}
    if state[2]>length(state[1])
        return nothing
    else
        return state[1][state[2]], (state[1], state[2]+1)
    end
end


text:
"""
        set!(widget, prop, value)

Sets the `value` of the Gtk property `prop` of `widget`.

# Examples

```julia
julia> set!(window, :title, "My window")

julia> set!(entry, :text, "Foo bar")
```
"""
function set!(widget::GtkWidget, prop::SString, value::T) where {T}
    set_gtk_property!(widget, prop, value)
end

"""
        set!(widget, prop, values::NamedTuple)

Sets a Gtk property of `widget` for each key-value pair in `values`.

The name of the final property will be `prop-key` where key is a key of `values`.

For example, instead of setting the margins of a `widget` each time, like this:

```julia
set!(mywidget, :margin_top, 10)
set!(mywidget, :margin_bottom, 20)
set!(mywidget, :margin_left, 15)
set!(mywidget, :margin_right, 10)
```

We can simple do:

```julia
set!(mywidget, :margin, (
    top = 10,
    bottom = 20,
    left = 15,
    right = 10
))
```
"""
function set!(widget::GtkWidget, prop::SString, values::NamedTuple)
    for (sprop, value) in pairs(values)
        set!(widget, "$prop-$sprop", value)
    end
end

"""
        set!(widget, props::NamedTuple)

Sets a Gtk property of `widget` for each key-value pair in `props`.

# Examples

```julia
julia> set!(button, (label = "Click-me"))

julia> set!(color_btn, (color = colorant"#f1a"))
```
"""
function set!(widget::GtkWidget, props::NamedTuple)
    for (prop, value) in props
        set!(widget, prop, value)
    end
end

"""
        set!(widget, props::Dict)

Sets a Gtk property of `widget` for each key-value pair in `props`.
"""
function set!(widget::GtkWidget, props::Dict)
    for (prop, value) in props
        set!(widget, prop, value)
    end
end

"""
        set!(widget; props...)

Sets a Gtk property of `widget` for each key-value pair in `props`.
"""
function set!(widget::GtkWidget; props...)
    if length(props) > 0
        set!(widget, Dict(pairs(props)))
    end
end

"""
        getprop(widget, key)

Returns the Gtk property `key` from `widget`.

# Examples

```julia
julia> getprop(window, :title)
"My window"
```
"""
function getprop(widget::GtkWidget, key::SString)
    return get_gtk_property(widget, key)
end
"""
        getprop(widget, key, T)

Returns the Gtk property `key` with type `T` from `widget`.

# Examples

```julia
julia> getprop(window, :title, String)
"My window"
```
"""
function getprop(widget::GtkWidget, key::SString, ::Type{T}) where T
    return get_gtk_property(widget, key, T)
end

"""
        onevent(callback, event, widget)

Adds a event listener for `widget`. The `callback` is called when the `event` is fired.

# Examples

```julia
julia> onevent(:destroy, window) do
            println("The window has been closed")
       end

```
"""
function onevent(callback::Function, event::SString, widget::GtkWidget)
    return signal_connect(widget, event) do args...
        try
            if applicable(callback, args...)
                callback(args...)
            else
                callback()
            end
        catch e
            @error "Error in event '$event' callback." exception=e
            Base.show_backtrace(stderr, catch_backtrace())
        end
    end
end

function onevent(callback::Function, event::SString, widgets::Vector{<:GtkWidget})
    return [onevent(callback, event, widget) for widget in widgets]
end

function set!(widget::GtkWidget, event::SString, callback::Function)
    if startswith(string(event), "on")
        onevent(callback, String(string(event)[3:end]), widget)
    else
        set_gtk_property!(widget, event, callback)
    end
end

"""
        offevent(widget, eventid)

Disconnect the event of id `eventid` from `widget`.

# Examples

```julia
julia> id = onevent(:destroy, window) do
            println("The window has been closed")
       end

julia> destroy(window)
The window has been closed

julia> offevent(window, id)

julia> show(window)
GtkWindow...

julia> destroy(window)

```
"""
function offevent(widget, id)
    signal_handler_disconnect(widget, id)
end

macro onevent(event, widget, callback)
    args = esc.([callback, event, widget])
    :( onevent($(args...)) )
end

function add!(parent::GtkWidget, child::GtkWidget)
    push!(parent, child)
end

function add!(parent::GtkWidget, children::GtkWidget...)
    for child in children
        add!(parent, child)
    end
end

function add!(parent::GtkWidget, children::Union{Tuple, Vector})
    add!(parent, children...)
end

function add!(parent::GtkWidget, children::Function)
    add!(parent, children())
end

function value(slider::Union{GtkScale, GtkSpinButton})
    return GAccessor.value(GAccessor.adjustment(slider))
end

function setvalue!(slider::Union{GtkScale, GtkSpinButton}, value::Real)
    adj = GAccessor.adjustment(slider)
    set_gtk_property!(adj, :value, value)
end

function value(w::Union{Gtk.GtkSwitch, GtkToggleButton, GtkCheckButton})
    return getprop(w, :active, Bool)
end

function setvalue!(w::Union{Gtk.GtkSwitch, GtkToggleButton, GtkCheckButton}, value::Bool)
    set!(w, :active, value)
end

function value(entry::GtkEntry)
    return bytestring(GAccessor.text(entry))
end

function value(cbtn::GtkColorButton)
    return convert(RGBA, getprop(cbtn, :rgba, Gtk.GdkRGBA))
end

function setvalue!(cbtn::GtkColorButton, color::Colorant)
    set!(cbtn, rgba = convert(GdkRGBA, color))
end


text:
#=
							Multiples of 3 and 5

	If we list all the natural numbers below 10 that are multiples of 3 or 5, we
	get 3, 5, 6 and 9. The sum of these multiples is 23.

	Find the sum of all the multiples of 3 or 5 below 1000.

							Answer: 233168
=#

function multiples(a::Int, b::Int, limit::Int)
    array = []

    for i in 2:limit-1
        if 0 == i%a || 0 == i%b
            push!(array, i)
        end
    end

    return array
end

println(sum(multiples(3, 5, 1000)))

max_stars_repo_name: Fazendaaa/project-euler

max_issues_repo_name: Fazendaaa/project-euler

max_forks_repo_name: Fazendaaa/project-euler


text:
function post(alg::BFFPSV18{N, ST}, ivp::IVP{<:AbstractContinuousSystem}, tspan;
              Δt0::TimeInterval=zeroI, kwargs...) where {N, ST}

    @unpack δ, approx_model, vars, block_indices,
            row_blocks, column_blocks = alg

    if haskey(kwargs, :NSTEPS)
        NSTEPS = kwargs[:NSTEPS]
        T = NSTEPS * δ
    else
        # get time horizon from the time span imposing that it is of the form (0, T)
        T = _get_T(tspan, check_zero=true, check_positive=true)
        NSTEPS = ceil(Int, T / δ)
    end

    # normalize system to canonical form
    # x' = Ax, x in X, or
    # x' = Ax + u, x in X, u in U
    ivp_norm = _normalize(ivp)

    # homogenize the initial-value problem
    if haskey(kwargs, :homogenize) && kwargs[:homogenize] == true
        ivp_norm = homogenize(ivp_norm)
    end

    # discretize system
    ivp_discr = discretize(ivp_norm, δ, approx_model)
    Φ = state_matrix(ivp_discr)
    Ω0 = initial_state(ivp_discr)
    X = stateset(ivp_discr)

    # true <=> there is no input, i.e. the system is of the form x' = Ax, x ∈ X
    got_homogeneous = !hasinput(ivp_discr)

    # decompose the initial states into a cartesian product
    # TODO add option to do the lazy decomposition
    Xhat0 = _decompose(Ω0, column_blocks, ST)
    Φ = state_matrix(ivp_discr)
    X = stateset(ivp_discr) # invariant

    # force using sparse type for the matrix exponential
    if alg.sparse
        Φ = SparseArrays.sparse(Φ)
    end

    # variables that are actually computed
    vars = reduce(vcat, alg.row_blocks)

    # preallocate output flowpipe
    CP = CartesianProductArray{N, ST}
    F = Vector{SparseReachSet{N, CP, length(vars)}}(undef, NSTEPS)

    # option to use array views
    viewval = Val(alg.view)

    if got_homogeneous
        reach_homog_BFFPSV18!(F, Xhat0, Φ, NSTEPS, δ, X, ST,
                              vars, block_indices,
                              row_blocks, column_blocks, Δt0, viewval)

    else
        U = inputset(ivp_discr)
        @assert isa(U, LazySet) "expected input of type `<:LazySet`, but got $(typeof(U))"
        reach_inhomog_BFFPSV18!(F, Xhat0, Φ, NSTEPS, δ, X, U, ST,
                                vars, block_indices,
                                row_blocks, column_blocks, Δt0, viewval)
    end
    return Flowpipe(F)
end


text:
# ```@meta
# CurrentModule = AtomicMiscellany
# ```
#
# # [Notes: hydrogenic energies](@id notes-hydrogenic)

using Pkg #src
Pkg.activate(joinpath(@__DIR__, "..", "..")) #src
using AtomicMiscellany: AtomicMiscellany, hydrogenic_energy, UniformShellNucleus, NRElectron, DiracElectron, α #hide
using Plots, Formatting #hide

# ## Point nucleus
#
# The energies for point nuclei (represented with [`PointNucleus`](@ref)) are calculated using analytic formulae.
#
# For a non-relativisic particles ([`NRParticle`](@ref)) the energy is given by
#
# ```math
# E_n = - \frac{(Z\alpha)^2}{2 n^2} mc^2
# ```
#
# and for relativistic particles ([`DiracParticle`](@ref)) by
#
# ```math
# m c^2 \left\{\left[1 + \left(
# \frac{Z\alpha}{n - |\kappa| + \sqrt{\kappa^2 - (Z\alpha)^2}}
# \right)^2 \right]^{-1/2} - 1 \right\}
# ```
#
# The only thing to note is that for the relativistic case, we set the ``E = 0`` for ``Z = 0``, so that the convention for the energy scale for the relativistic and non-relativistic cases would be the same.

# ## Uniform spherical shell
#
# There are no analytical expressions for the energies for finite nuclei. However, in the case of the uniform spherical shell ([`UniformShellNucleus`](@ref); aka. cut-off Coulomb, top-slice model), the potential is simple and the equations solvable in the ranges ``[0, R]`` and ``[R, \infty)`` separately.
#
# Following the expressions in [`Andrae2000`](@ref), we implement a root-finding based algorithm for determining the energies, in both the relativistic and non-relativistic case.

## We use this function below to showcase the performance of the implementation of hydrogenic_energies.
function plot_fnc_energies(p; n,
    Rs = [0.5, 1e-1, 5e-2, 2e-2, 1e-2, 7e-3, 3e-3, 1e-3, 1e-4],
    qnumbers...
)
    Z = range(1, isa(p, AtomicMiscellany.DiracParticle) ? α^-1 : 150, length=1001)
    emin = hydrogenic_energy.(p, Z; n = n, qnumbers...)
    emax = hydrogenic_energy.(p, Z; n = n + 1, qnumbers...)
    p1 = plot(hcat(Z,Z), hcat(emin,emax), c=:black, ls=:dash, legend=:bottomleft, label=false)
    p2 = plot(hcat(Z,Z), hcat(emin ./ abs.(emin), emax ./ abs.(emin)), label=false, c=:black, ls=:dash)
    for (i, R) = enumerate(Rs)
        es = map(Z) do Z
            try
                hydrogenic_energy(p, UniformShellNucleus(R), Z; n = n, qnumbers...)
            catch e
                return NaN
            end
        end
        plot!(p1, Z, es, label="R=$R", c = i)
        plot!(p2, Z, es./ abs.(emin), label=false, c = i)
    end
    ylabel!(p1, "E"); ylabel!(p2, "E / E(PNC)")
    xlabel!(p2, "Nuclear charge Z")
    plot(p1, p2, layout=(2,1), size=(800, 600))
    vline!([α^-1 α^-1], label=false, c=:black, ls=:dot)
    title!("n = $n, $(qnumbers...)", subplot=1)
end
nothing # hide

# ### Non-relativisic energies
#
# Calculates the matching function ``f`` derived from [`Andrae2000`](@ref), equations (247-248, 256-258). ``f`` is derived by multiplying the equation through with the denominators to remove any infinities these denominators might bring.
#
# The function `AtomicMiscellany.f_andrae_nonrel_topslice` implements the following function, which is then fed into a root-finder:
#
# ```math
# \begin{aligned}
# f(E, r, Z, R, \ell) &=
# \beta j_\ell(x_{\rm{in}}) U(a, b, x_{\rm{out}})
# \\ &+ 2 \beta a j_\ell(x_{\rm{in}}) U(a+1, b+1, x_{\rm{out}})
# \\ &- \alpha j_{\ell+1}(x_{\rm{in}}) U(a, b, x_{\rm{out}}) ≡ 0
# \end{aligned}
# \\[1em]
# \alpha_{n\ell} = \sqrt{2(E + Z/R)}, \quad
# x_{\rm{in}} = \alpha_{n\ell} r, \quad
# \beta_{n\ell} = \sqrt{-2E}, \quad
# x_{\rm{out}} = 2 \beta_{n\ell} r,
# \\[1em]
# a = \frac{b}{2} - \nu_{n\ell}, \quad
# b = 2(\ell + 1), \quad
# \nu_{n\ell} = Z/\beta_{n\ell}
# ```

## We can use this function to visualize f:
function plot_nonrel_fnc_f(; R, Z, ℓ)
    emin = hydrogenic_energy(NRElectron, Z; n = 1)
    ## Excluding the first point, since at E = -Z/R, f is zero, which causes problems
    es = range(max(1.1emin, -Z/R), 0, length=502)[2:end-1]
    ys = (E -> AtomicMiscellany.f_andrae_nonrel_topslice(Z=Z, r=R, R=R, ℓ=ℓ, E = E)).(es)
    any(iszero.(ys)) && @warn "Zeros in ys"
    plot(es, abs.(ys), yaxis=:log10, legend=false, c = (y -> y>0 ? 2 : 1).(ys))
    absy_sorted = filter(!iszero, sort(abs.(ys)))
    ylims!(absy_sorted[1]/2, absy_sorted[490])
    e0s = map(ℓ .+ (1:6)) do n
        try
            hydrogenic_energy(NRElectron, UniformShellNucleus(R), Z; n = n, ℓ = ℓ, verbose=false)
        catch
            NaN
        end
    end
    xlabel!("E"); ylabel!("log |f|")
    title!("Z=$(fmt(".2f", Z)), ℓ=$ℓ, R=$(fmt(".2e", R))")
    vline!(e0s, c=:gray)
    vline!([hydrogenic_energy(NRElectron, Z; n = n) for n = 1:3], c=:black, ls=:dash)
    plot!(size=(800, 400))
end
nothing # hide

# Here is ``f(E)`` for ``Z = 1`` and ``\ell = 0``, for two very different nuclear sizes.
# ``R = 10^{-4}~\mathrm{a.u.}`` is approximately physical, whereas ``R = 2.5~\mathrm{a.u.}`` is unphysically huge. However, the latter is useful for showing how the graph of the function shifts as we make the nucleus bigger and bigger.
#
# The black dashed vertical lines indicate the hydrogenic energies for point nuclei, which can be determined using explicit formulae. As we need to specify a bracketing interval for the root-finder, we use those energies.
# This seems like a reasonable choice, as we know that the FNC energy will always we higher than the corresponding PNC energy, but unlikely to reach the PNC energy of the 2s, unless the parameters are very extreme.

plot(
    plot_nonrel_fnc_f(R = 1e-4, Z = 1, ℓ = 0),
    plot_nonrel_fnc_f(R = 2.5, Z = 1, ℓ = 0),
    layout = (2, 1), size=(800, 600),
)

# To illustrate the potential problems with the bracketing intervals, let's see what happens (for a huge nucleus) when we increase ``Z`` a lot. Eventually, we will hit the `2s` energy. However, this should not really happen in any physically relevant cases.

a = @animate for Z = range(1, 550, length=101)
    plot_nonrel_fnc_f(R = 1e-2, Z = Z, ℓ = 0)
end
gif(a, fps=10)

# To showcase the algorithm an its limits, here are the `1s` energies for various nuclear sizes (`R` values).
# The black dashed lines once again indicate the bracketing interval of PNC energies.

plot_fnc_energies(NRElectron, n = 1, ℓ = 0)

# Now, if we do the same thing for ``n = 2``, we start seeing weird cases where

plot_fnc_energies(NRElectron, n = 2, ℓ = 0)

# This can easily be understood by looking at the behaviour of ``f`` as you increase ``Z`` for a large nucleus (`R = 0.1`):

a = @animate for Z = range(25, 100, length=101)
    plot_nonrel_fnc_f(R = 0.1, Z = Z, ℓ = 0)
end
gif(a, fps=10)

# At one point, the zero corresponding to the `1s` energy enters the bracketing interval of the `2s`. Root-finder, of course, can not distiguish that, so it just returns the `1s` energy. Again though, this is only an issue if the value of ``R`` or ``Z`` gets extremely large.

# We see the same happening for higher ``\ell`` values too. One thing to note is that we also get in trouble for the lowest eigenvalue, because the ``E = -Z/R``, below which ``f(E)`` becomes complex, enters the ``n=2`` bracketing interval.

plot_fnc_energies(NRElectron, n = 2, ℓ = 1)

#-

a = @animate for Z = range(25, 100, length=101)
    plot_nonrel_fnc_f(R = 0.1, Z = Z, ℓ = 1)
end
gif(a, fps=10)

# However, again, for reasonable ``R`` and ``Z`` values, everything behaves reasonably.
# Also, correctly, we do not get any zeros at lower ``n`` values for higher ``\ell`` values:

plot(
    plot_nonrel_fnc_f(R = 1e-4, Z = 1, ℓ = 0),
    plot_nonrel_fnc_f(R = 1e-4, Z = 1, ℓ = 1),
    plot_nonrel_fnc_f(R = 1e-4, Z = 1, ℓ = 2),
    layout = (3, 1), size=(800, 900),
)

# ### Relativistic energies
#
# The relativistic energy is determined in essentially the same way as in the non-relativistic case.
# The relevant equation in [`Andrae2000`](@ref) is (284), setting up the matching condition between the inner and outer parts of the orbital:
#
# ```math
# L^{\mathrm{in}}_{n\kappa}(R) - \frac{\gamma}{R} = - T^{\mathrm{out}}_{n\kappa}(R)
# ```
#
# ``L^{\mathrm{in}}_{n\kappa}(r)`` is the logarithmic derivative of the ``P`` component of the inner part of the wavefunction (i.e. for ``r \leq R``).
# According to equations (273) and (274), it can be written as
#
# ```math
# L^{\mathrm{in}}_{n\kappa}(r) = \frac{\ell + 1}{r} - T^{\mathrm{in}}_{n\kappa}(r)
# = \frac{\ell + 1}{r} - \alpha_{n\kappa} \frac{j_{\ell+1}(x_{\mathrm{in}})}{j_{\ell}(x_{\mathrm{in}})}
# \\[1em]
# \alpha_{n\kappa} = \sqrt{2 (E + Z/R) + [(E + Z/R)/c]^2},
# \quad
# x_{\mathrm{in}} = \alpha_{n\kappa} r
# ```
#
# ``T^{\mathrm{out}}_{n\kappa}(r)``, which is related to the outer part of the orbital, is given in equation (283)
#
# ```math
# T^{\mathrm{out}}_{n\kappa}(r) = \beta_{n\kappa}
# + 2 \beta_{n\kappa} \frac{
#     a U(a+1, b+1, x_{\mathrm{out}}) + A \cdot (a+1) U(a+2, b+1, x_{\mathrm{out}})
# }{
#     U(a, b, x_{\mathrm{out}}) + A \cdot U(a + 1, b, x_{\mathrm{out}})
# }
# \\[1em]
# x_{\mathrm{out}} = 2 \beta_{n\kappa} r, \quad
# \beta_{n\kappa} = C^+ C^- / c, \quad
# C^{\pm} = \sqrt{c^2 \pm W_{n\kappa}}
# \\[1em]
# W_{n\kappa} = E + c^2, \quad
# v_{n\kappa} = \frac{1}{2} + \frac{Z W_{n\kappa}}{\beta_{n\kappa} c^2}, \quad
# A = \frac{Z}{\beta_{n\kappa}} + \kappa
# \\[1em]
# a = b/2 - v_{n\kappa}, \quad
# b = 2\gamma + 1, \quad
# \gamma = \sqrt{\kappa^2 - (Z/c)^2}
# ```
#
# We can get rid of the denominators, which may cause some unwanted singularities, and write down the following condition
#
# ```math
# \begin{aligned}
# f(E) & =
# \left(
#     \frac{\ell + 1 - \gamma}{R} + \beta_{n\kappa}
# \right) (U_{00} + A U_{10}) J_0 \\
# & - \alpha_{n\kappa} (U_{00} + A U_{10}) J_1 \\
# & + 2 \beta_{n\kappa} [a U_{11} + A (a + 1) U_{21}] ≡ 0
# \end{aligned}
# \\[1em]
# J_i = j_{\ell + 1}(x_{\mathrm{in}}), \quad
# U_{ij} = U(a + i, b + j, x_{\mathrm{out}})
# ```
#
# The special functions ``j_\ell(x)`` and ``U(a, b, x)`` are the spherical bessel and irregular confluent hypergeometric functions, respectively.
# They are evaluated using the corresponding functions from GSL ([`sf_bessel_jl`](https://www.gnu.org/software/gsl/doc/html/specfunc.html?highlight=sf_bessel_jl#c.gsl_sf_bessel_jl), [`sf_hyperg_U`](https://www.gnu.org/software/gsl/doc/html/specfunc.html?highlight=sf_hyperg_u#c.gsl_sf_hyperg_U)), via [GSL.jl](https://github.com/JuliaMath/GSL.jl).

## The following function can be used to visualize ``f``.
function plot_rel_fnc_f(; R, Z, κ, energies = nothing, plotfn=plot)
    es = if isnothing(energies)
        emin = hydrogenic_energy(DiracElectron, Z; n = 1, κ = κ)
        ## Excluding the first point, since at E = -Z/R, f is zero, which causes problems
        es = range(max(1.1emin, -Z/R), 0, length=502)[2:end-1]
    else
        energies
    end
    ys = (E -> AtomicMiscellany.f_andrae_dirac_topslice(Z=Z, R=R, κ=κ, E = E)).(es)
    any(iszero.(ys)) && @warn "Zeros in ys"
    plotfn(es, abs.(ys), yaxis=:log10, legend=false, c = (y -> y>0 ? 2 : 1).(ys))
    if isnothing(energies)
        absy_sorted = filter(!iszero, sort(abs.(ys)))
        ylims!(absy_sorted[1]/2, absy_sorted[490])
    end
    e0s = map(AtomicMiscellany.κ2ℓ(κ) .+ (1:6)) do n
        try
            hydrogenic_energy(DiracElectron, UniformShellNucleus(R), Z; n = n, κ = κ, verbose=false)
        catch e
            @warn "Unable to find root for n=$n" error=(e, catch_backtrace())
            NaN
        end
    end
    xlabel!("E"); ylabel!("log |f|")
    title!("Z=$(fmt(".2f", Z)), κ=$κ, R=$(fmt(".2e", R))")
    vline!(e0s, c=:gray)
    vline!([hydrogenic_energy(DiracElectron, Z; n = n, κ = κ) for n = 1:6], c=:black, ls=:dash)
    plot!(size=(800, 400))
end
nothing # hide

# The general behavior of ``f`` is very similar to the non-relativistic case.
# However, the expressions given in [`Andrae2000`](@ref) require us to evaluate ``\gamma = \sqrt{\kappa^2 - (Z/c)^2}``, so the expressions would become complex for ``Z > α^{-1}``. For that reason we restrict ourselves to ``Z \leq α^{-1}``.

plot(
    plot_rel_fnc_f(R = 1e-4, Z = 1, κ = -1),
    plot_rel_fnc_f(R = 5e-1, Z = 1, κ = -1),
    layout = (2, 1), size=(800, 600),
)

#-

plot_fnc_energies(DiracElectron, n = 1, κ = -1)

# As in the non-relativistic case, we run into issues with the lower energy zeros coming into the bracketing interval once the parameters become extreme enough (as we're limited in ``Z``, this primarily means ``R`` in this case).

plot_fnc_energies(DiracElectron, n = 2, κ = -1)
#-
plot_fnc_energies(DiracElectron, n = 3, κ = -1)
#-
plot_fnc_energies(DiracElectron, n = 2, κ = 1)

# We do run into trouble if we try to find the zeros for higher angular momenta:

plot_fnc_energies(DiracElectron, n = 2, κ = -2)
#-
plot_fnc_energies(DiracElectron, n = 3, κ = 2)

# This, turns out, is due to the limited resolution of 64-bit floats. We can zoom in on around the `1s` zero. We see that ``f(E)`` evaluates to the same value for many ``E`` values, implying a rounding problem. This, in turn, leads to the root-finder not being happy, since the bracketing interval does not encompass any positive ``f(E)`` values.

let R = 1e-3, Z = 1, κ = -2,
    e0 = hydrogenic_energy(DiracElectron, Z; n = 2, κ = κ),
    (emin, emax) = e0 .+ (-1, 1) .* 0.000_000_000_05
    plot_rel_fnc_f(R=R, Z=Z, κ=κ, energies = range(emin, emax, length=151), plotfn=scatter)
    xlims!(emin, emax)
end

# There are workarounds (modify the mathematical expressions for better numerical stability, use `BigFloat`s), but for now this is simply a limitation of the current implementation.
#
# Also, the non-relativistic version should also run into this issue at high enough ``\ell``.


text:
# abstracttensor.jl
#
# Abstract Tensor type
#----------------------
"""
    abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end

Abstract supertype of all tensor maps, i.e. linear maps between tensor products
of vector spaces of type `S<:IndexSpace`. An `AbstractTensorMap` maps from
an input space of type `ProductSpace{S,N₂}` to an output space of type
`ProductSpace{S,N₁}`.
"""
abstract type AbstractTensorMap{S<:IndexSpace, N₁, N₂} end
"""
    AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{T,S,N,0}

Abstract supertype of all tensors, i.e. elements in the tensor product space
of type `ProductSpace{S,N}`, built from elementary spaces of type `S<:IndexSpace`.

An `AbstractTensor{S,N}` is actually a special case `AbstractTensorMap{S,N,0}`,
i.e. a tensor map with only a non-trivial output space.
"""
const AbstractTensor{S<:IndexSpace, N} = AbstractTensorMap{S, N, 0}

# tensor characteristics
Base.@pure Base.eltype(T::Type{<:AbstractTensorMap}) = eltype(storagetype(T))
Base.@pure similarstoragetype(TT::Type{<:AbstractTensorMap}, ::Type{T}) where {T} =
    Core.Compiler.return_type(similar, Tuple{storagetype(TT), Type{T}})

storagetype(t::AbstractTensorMap) = storagetype(typeof(t))
similarstoragetype(t::AbstractTensorMap, T) = similarstoragetype(typeof(t), T)
Base.eltype(t::AbstractTensorMap) = eltype(typeof(t))
spacetype(t::AbstractTensorMap) = spacetype(typeof(t))
sectortype(t::AbstractTensorMap) = sectortype(typeof(t))
field(t::AbstractTensorMap) = field(typeof(t))
numout(t::AbstractTensorMap) = numout(typeof(t))
numin(t::AbstractTensorMap) = numin(typeof(t))
numind(t::AbstractTensorMap) = numind(typeof(t))

Base.@pure spacetype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = S
Base.@pure sectortype(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = sectortype(S)
Base.@pure field(::Type{<:AbstractTensorMap{S}}) where {S<:IndexSpace} = field(S)
Base.@pure numout(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁, N₂} = N₁
Base.@pure numin(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁, N₂} = N₂
Base.@pure numind(::Type{<:AbstractTensorMap{<:IndexSpace,N₁,N₂}}) where {N₁, N₂} = N₁ + N₂

const order = numind

# tensormap implementation should provide codomain(t) and domain(t)
codomain(t::AbstractTensorMap, i) = codomain(t)[i]
domain(t::AbstractTensorMap, i) = domain(t)[i]
source(t::AbstractTensorMap) = domain(t) # categorical terminology
target(t::AbstractTensorMap) = codomain(t) # categorical terminology
space(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂} =
    i <= N₁ ? codomain(t,i) : dual(domain(t, i-N₁))

space(t::AbstractTensor) = codomain(t)
space(t::AbstractTensor, i) = space(t)[i]

# some index manipulation utilities
codomainind(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}) where {N₁,N₂} =
    ntuple(n->n, StaticLength(N₁))
domainind(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}) where {N₁,N₂} =
    ntuple(n->N₁+n, StaticLength(N₂))

adjointtensorindex(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂} =
    ifelse(i<=N₁, N₂+i, i-N₁)

adjointtensorindices(t::AbstractTensorMap, I::IndexTuple) =
    map(i->adjointtensorindex(t, i), I)

# NOTE: do we still need this
tensor2spaceindex(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂} =
    ifelse(i<=N₁, i, 2N₁+N₂+1-i)
space2tensorindex(t::AbstractTensorMap{<:IndexSpace,N₁,N₂}, i) where {N₁,N₂} =
    ifelse(i<=N₁, i, 2N₁+N₂+1-i)

# Defining vector spaces:
#------------------------
const TensorSpace{S<:IndexSpace} = Union{S, ProductSpace{S}}
const TensorMapSpace{S<:IndexSpace, N₁, N₂} = Pair{ProductSpace{S,N₂},ProductSpace{S,N₁}}

# Little unicode hack to define TensorMapSpace
→(dom::ProductSpace{S}, codom::ProductSpace{S}) where {S<:IndexSpace} = dom => codom
→(dom::S, codom::ProductSpace{S}) where {S<:IndexSpace} = ProductSpace(dom) => codom
→(dom::ProductSpace{S}, codom::S) where {S<:IndexSpace} = dom => ProductSpace(codom)
→(dom::S, codom::S) where {S<:IndexSpace} = ProductSpace(dom) => ProductSpace(codom)

←(codom::ProductSpace{S}, dom::ProductSpace{S}) where {S<:IndexSpace} = dom => codom
←(codom::S, dom::ProductSpace{S}) where {S<:IndexSpace} = dom => ProductSpace(codom)
←(codom::ProductSpace{S}, dom::S) where {S<:IndexSpace} = ProductSpace(dom) => codom
←(codom::S, dom::S) where {S<:IndexSpace} = ProductSpace(dom) => ProductSpace(codom)

# NOTE: do we still need this
function blocksectors(codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₂}) where {S,N₁,N₂}
    sectortype(S) === Trivial && return (Trivial(),)
    return intersect(blocksectors(codom), blocksectors(dom))
end

# Conversion to Array:
#----------------------
# probably not optimized for speed, only for checking purposes
function Base.convert(::Type{Array}, t::AbstractTensorMap{S,N₁,N₂}) where {S,N₁,N₂}
    G = sectortype(t)
    if G === Trivial
        convert(Array, t[])
    elseif FusionStyle(G) isa Abelian || FusionStyle(G) isa SimpleNonAbelian
        # TODO: Frobenius-Schur indicators!, and fermions!
        cod = codomain(t)
        dom = domain(t)
        A = fill(zero(eltype(t)), (dims(cod)..., dims(dom)...))
        for (f1,f2) in fusiontrees(t)
            F1 = convert(Array, f1)
            F2 = convert(Array, f2)
            for i = 1:N₁
                if isdual(cod[i])
                    a = f1.uncoupled[i]
                    Z = sqrt(dim(a))*permutedims(conj(reshape(fusiontensor(a,dual(a),one(a)), (dim(a),dim(a)))),(2,1))
                    indF = ntuple(k->(k == i ? -i : k), StaticLength(N₁)+StaticLength(1))
                    indout = ntuple(identity, StaticLength(N₁)+StaticLength(1))
                    F1 = TensorOperations.tensorcontract(Z,(i,-i), F1, indF, indout)
                end
            end
            for i = 1:N₂
                if isdual(dom[i])
                    a = f2.uncoupled[i]
                    Z = sqrt(dim(a))*permutedims(conj(reshape(fusiontensor(a,dual(a),one(a)), (dim(a),dim(a)))),(2,1))
                    indF = ntuple(k->(k == i ? -i : k), StaticLength(N₂)+StaticLength(1))
                    indout = ntuple(identity, StaticLength(N₂)+StaticLength(1))
                    F2 = TensorOperations.tensorcontract(Z,(i,-i), F2, indF, indout)
                end
            end
            sz1 = size(F1)
            sz2 = size(F2)
            d1 = TupleTools.front(sz1)
            d2 = TupleTools.front(sz2)
            F = reshape(reshape(F1, TupleTools.prod(d1), sz1[end])*reshape(F2, TupleTools.prod(d2), sz2[end])', (d1...,d2...))
            Aslice = StridedView(A)[axes(cod, f1.uncoupled)..., axes(dom, f2.uncoupled)...]
            axpy!(1, StridedView(_kron(convert(Array, t[f1,f2]), F)), Aslice)
        end
        return A
    else
        # TODO: implement DegenerateNonAbelian case
        throw(MethodError(convert, (Array, t)))
    end
end
# TODO: Reverse conversion


text:
"""
    InputAugmentation(analyzer, n, [std=1, rng=GLOBAL_RNG])
    InputAugmentation(analyzer, n, distribution, [rng=GLOBAL_RNG])

A wrapper around analyzers that augments the input with `n` samples of additive noise sampled from `distribution`.
This input augmentation is then averaged to return an `Explanation`.
"""
struct InputAugmentation{A<:AbstractXAIMethod,D<:Distribution,R<:AbstractRNG} <:
       AbstractXAIMethod
    analyzer::A
    n::Integer
    distribution::D
    rng::R
end
function InputAugmentation(analyzer, n, distr, rng=GLOBAL_RNG)
    return InputAugmentation(analyzer, n, distr, rng)
end
function InputAugmentation(analyzer, n, σ::Real=0.1f0, args...)
    return InputAugmentation(analyzer, n, Normal(0.0f0, Float32(σ)^2), args...)
end

function (aug::InputAugmentation)(input, ns::AbstractNeuronSelector)
    # Regular forward pass of model
    output = aug.analyzer.model(input)
    output_indices = ns(output)

    # Call regular analyzer on augmented batch
    augmented_input = add_noise(augment_batch_dim(input, aug.n), aug.distribution, aug.rng)
    augmented_indices = augment_indices(output_indices, aug.n)
    augmented_expl = aug.analyzer(augmented_input, AugmentationSelector(augmented_indices))

    # Average explanation
    return Explanation(
        reduce_augmentation(augmented_expl.attribution, aug.n),
        output,
        output_indices,
        augmented_expl.analyzer,
        Nothing,
    )
end

function add_noise(A::AbstractArray{T}, distr::Distribution, rng::AbstractRNG) where {T}
    return A + T.(rand(rng, distr, size(A)))
end

"""
    augment_batch_dim(input, n)

Repeat each sample in input batch n-times along batch dimension.
This turns arrays of size `(..., B)` into arrays of size `(..., B*n)`.

## Example
```julia-repl
julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4

julia> augment_batch_dim(A, 3)
2×6 Matrix{Int64}:
 1  1  1  2  2  2
 3  3  3  4  4  4
```
"""
function augment_batch_dim(input::AbstractArray{T,N}, n) where {T,N}
    return repeat(input; inner=(ntuple(Returns(1), N - 1)..., n))
end

"""
    reduce_augmentation(augmented_input, n)

Reduce augmented input batch by averaging the explanation for each augmented sample.
"""
function reduce_augmentation(input::AbstractArray{T,N}, n) where {T<:AbstractFloat,N}
    # Allocate output array
    in_size = size(input)
    in_size[end] % n != 0 &&
        throw(ArgumentError("Can't reduce augmented batch size of $(in_size[end]) by $n"))
    out_size = (in_size[1:(end - 1)]..., div(in_size[end], n))
    out = similar(input, eltype(input), out_size)

    axs = axes(input, N)
    inds_before_N = ntuple(Returns(:), N - 1)
    for (i, ax) in enumerate(first(axs):n:last(axs))
        view(out, inds_before_N..., i) .=
            sum(view(input, inds_before_N..., ax:(ax + n - 1)); dims=N) / n
    end
    return out
end
"""
    augment_indices(indices, n)

Strip batch indices and return inidices for batch augmented by n samples.

## Example
```julia-repl
julia> inds = [CartesianIndex(5,1), CartesianIndex(3,2)]
2-element Vector{CartesianIndex{2}}:
 CartesianIndex(5, 1)
 CartesianIndex(3, 2)

julia> augment_indices(inds, 3)
6-element Vector{CartesianIndex{2}}:
 CartesianIndex(5, 1)
 CartesianIndex(5, 2)
 CartesianIndex(5, 3)
 CartesianIndex(3, 4)
 CartesianIndex(3, 5)
 CartesianIndex(3, 6)
```
"""
function augment_indices(inds::Vector{CartesianIndex{N}}, n) where {N}
    indices_wo_batch = [i.I[1:(end - 1)] for i in inds]
    return map(enumerate(repeat(indices_wo_batch; inner=n))) do (i, idx)
        CartesianIndex{N}(idx..., i)
    end
end