text:
# Lambert's Problem

using BoundaryValueDiffEq
using DiffEqBase, OrdinaryDiffEq, LinearAlgebra
using NLsolve
using Test

y0=[-4.7763169762853989E+06, -3.8386398704441520E+05, -5.3500183933132319E+06, -5528.612564911408, 1216.8442360202787, 4845.114446429901]
init_val = [-4.7763169762853989E+06,-3.8386398704441520E+05,-5.3500183933132319E+06,7.0526926403748598E+06,-7.9650476230388973E+05,-1.1911128863666430E+06]
J2=1.08262668E-3
req=6378137
myu=398600.4418E+9
t0=86400*2.3577475462484435E+04
t1=86400*2.3577522023524125E+04
tspan = (t0,t1)
# ODE solver
function orbital(dy,y,p,t)
  r2=(y[1]^2 + y[2]^2 + y[3]^2)
  r3=r2^(3/2)
  w=1+1.5J2*(req*req/r2)*(1-5y[3]*y[3]/r2)
  w2=1+1.5J2*(req*req/r2)*(3-5y[3]*y[3]/r2)
  dy[1] = y[4]
  dy[2] = y[5]
  dy[3] = y[6]
  dy[4] = -myu*y[1]*w/r3
  dy[5] = -myu*y[2]*w/r3
  dy[6] = -myu*y[3]*w2/r3
end

function bc!_generator(resid,sol,init_val)
  resid[1] = sol[1][1]   - init_val[1]
  resid[2] = sol[1][2]   - init_val[2]
  resid[3] = sol[1][3]   - init_val[3]
  resid[4] = sol[end][1] - init_val[4]
  resid[5] = sol[end][2] - init_val[5]
  resid[6] = sol[end][3] - init_val[6]
end
cur_bc! = (resid,sol,p,t) -> bc!_generator(resid,sol,init_val)
resid_f = Array{Float64}(undef, 6)

### Test the IVP Near the true solution
### Should be small
# prob = ODEProblem(f,y0,tspan)
# sol = solve(prob,DP5(),abstol=1e-13,reltol=1e-13)
# cur_bc!(resid_f,sol)

TestTol = 0.05

### Now use the BVP solver to get closer
bvp = BVProblem(orbital, cur_bc!, y0, tspan)
@time sol = solve(bvp, Shooting(DP5(),nlsolve=(f,u0) -> (res=NLsolve.nlsolve(f,u0,autodiff=:central,ftol=1e-13);(res.zero, res.f_converged))),force_dtmin=true,abstol=1e-13,reltol=1e-13)
cur_bc!(resid_f,sol,nothing,sol.t)
@test norm(resid_f, Inf) < TestTol

@time sol = solve(bvp, Shooting(DP5()),force_dtmin=true,abstol=1e-13,reltol=1e-13)
cur_bc!(resid_f,sol,nothing,sol.t)
@test norm(resid_f, Inf) < TestTol

@time sol = solve(bvp, Shooting(DP5(),nlsolve=(f,u0) -> (res=NLsolve.nlsolve(f,u0,autodiff=:forward,ftol=1e-13,xtol=1e-13);
(res.zero, res.f_converged))),force_dtmin=true,abstol=1e-13,reltol=1e-13)
cur_bc!(resid_f,sol,nothing,sol.t)
@test norm(resid_f, Inf) < TestTol


text:
using Tables, PrettyTables
import EzXML.parsehtml

function amlTable(x)
  io = IOBuffer()
  pretty_table(io, x, backend = :html, standalone = false)
  str = String(resize!(io.data, io.size))

  html = Base.findfirst("html/body/table",parsehtml(str))
  unlink!(html)
  return html
end

ext: jl


text:
using Turing
using Optim
using Test
using StatsBase
using NamedArrays
using ReverseDiff
using Random
using LinearAlgebra

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@testset "ModeEstimation.jl" begin
    @testset "MLE" begin
        Random.seed!(222)
        true_value = [0.0625031, 1.75]

        m1 = optimize(gdemo_default, MLE())
        m2 = optimize(gdemo_default, MLE(), NelderMead())
        m3 = optimize(gdemo_default, MLE(), Newton())
        m4 = optimize(gdemo_default, MLE(), true_value, Newton())
        m5 = optimize(gdemo_default, MLE(), true_value)

        @test all(isapprox.(m1.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m2.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m3.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m4.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m5.values.array - true_value, 0.0, atol=0.01))
    end

    @testset "MAP" begin
        Random.seed!(222)
        true_value = [49 / 54, 7 / 6]

        m1 = optimize(gdemo_default, MAP())
        m2 = optimize(gdemo_default, MAP(), NelderMead())
        m3 = optimize(gdemo_default, MAP(), Newton())
        m4 = optimize(gdemo_default, MAP(), true_value, Newton())
        m5 = optimize(gdemo_default, MAP(), true_value)

        @test all(isapprox.(m1.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m2.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m3.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m4.values.array - true_value, 0.0, atol=0.01))
        @test all(isapprox.(m5.values.array - true_value, 0.0, atol=0.01))
    end

    @testset "StatsBase integration" begin
        Random.seed!(54321)
        mle_est = optimize(gdemo_default, MLE())

        @test coefnames(mle_est) == ["s", "m"]

        diffs = coef(mle_est).array - [0.0625031; 1.75001]
        @test all(isapprox.(diffs, 0.0, atol=0.1))

        infomat = [0.003907027690416608 4.157954948417027e-7; 4.157954948417027e-7 0.03125155528962335]
        @test all(isapprox.(infomat - informationmatrix(mle_est), 0.0, atol=0.01))

        ctable = coeftable(mle_est)
        @test ctable isa StatsBase.CoefTable

        s = stderror(mle_est).array
        @test all(isapprox.(s - [0.06250415643292194, 0.17677963626053916], 0.0, atol=0.01))

        @test coefnames(mle_est) == params(mle_est)
        @test vcov(mle_est) == informationmatrix(mle_est)

        @test isapprox(loglikelihood(mle_est), -0.0652883561466624, atol=0.01)
    end

    @testset "Linear regression test" begin
        @model function regtest(x, y)
            beta ~ MvNormal(2,1)
            mu = x*beta
            y ~ MvNormal(mu, 1.0)
        end

        Random.seed!(987)
        true_beta = [1.0, -2.2]
        x = rand(40, 2)
        y = x*true_beta

        model = regtest(x, y)
        mle = optimize(model, MLE())

        vcmat = inv(x'x)
        vcmat_mle = informationmatrix(mle).array

        @test isapprox(mle.values.array, true_beta)
        @test isapprox(vcmat, vcmat_mle)
    end

    @testset "Dot tilde test" begin
        @model function dot_gdemo(x)
            s ~ InverseGamma(2,3)
            m ~ Normal(0, sqrt(s))

            (.~)(x, Normal(m, sqrt(s)))
        end

        model_dot = dot_gdemo([1.5, 2.0])

        mle1 = optimize(gdemo_default, MLE())
        mle2 = optimize(model_dot, MLE())

        map1 = optimize(gdemo_default, MAP())
        map2 = optimize(model_dot, MAP())

        @test isapprox(mle1.values.array, mle2.values.array)
        @test isapprox(map1.values.array, map2.values.array)
    end
end


text:
module ShaperOriginDesignLib

using Printf
using Unitful
using UnitfulUS

export SVG_NAMESPACE

SVG_NAMESPACE = "http://www.w3.org/2000/svg"

include("shaper_cut_attributes.jl")
include("svg_utils.jl")
include("custom_anchor.jl")

end

max_stars_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_issues_repo_name: MarkNahabedian/ShaperOriginDesignLib

max_forks_repo_name: MarkNahabedian/ShaperOriginDesignLib


text:
""" Média Ponderada
# Implemente um programa que lê quatro valores
# e calcule a média ponderada para pesos 1, 2, 3 e 4,
# respectivamente (multiplique cada nota pelo seu peso,
# some os produtos e divida o resultado pela soma dos pesos).
# Exemplo: 5*2 + 7*2 + 6*3 + 3*3
"""

println("Informe o primeiro produto")
avg = parse( Float16, readline() )
println("Informe o peso do primeiro produto")
weight = parse( Float16, readline() )
numerator = avg * weight
denominator = weight

println("Informe o segundo produto")
avg = parse( Float16, readline() )
println("Informe o peso do segundo produto")
weight = parse( Float16, readline() )
numerator += avg * weight
denominator += weight

println("Informe o terceiro produto")
avg = parse( Float16, readline() )
println("Informe o peso do terceiro produto")
weight = parse( Float16, readline() )
numerator += avg * weight
denominator += weight

println("Informe o quarto produto")
avg = parse( Float16, readline() )
println("Informe o peso do quarto produto")
weight = parse( Float16, readline() )
numerator += avg * weight
denominator += weight

avg = (numerator / denominator)
println( "A média ponderada = ",avg )

max_stars_repo_name: codinginbrazil/functional-programming

max_issues_repo_name: codinginbrazil/functional-programming

max_forks_repo_name: codinginbrazil/functional-programming


text:
using AdjacencyArrays, Test, SparseArrays

function tests()
	@testset "Empty" begin
		net = EmptyNetwork(10)
	    @test numnodes(net) == 10
	    @test numconnections(net) == 0
	    @test !isdirected(net)
	    @test connectivity(net, 1) == 0
		@test connectivities(net) == [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
	    @test connectivity(net) == 0
	    @test adjMat(net) == zeros(10, 10)
	    @test adjVet(net) == []
	end

	@testset "Global" begin
		# Undirected
		net = GlobalNetwork(10)
		@test numnodes(net) == 10
		@test numconnections(net) == 10 * 9 / 2
		@test !isdirected(net)
		@test connectivity(net, 1) == 9
		@test connectivities(net) == [9, 9, 9, 9, 9, 9, 9, 9, 9, 9]
		@test connectivity(net) == 9

		# Undirected
		net = GlobalNetwork(10, directed=true)
		@test numconnections(net) == 10 * 9
		@test isdirected(net)

		net = GlobalNetwork(4)
		@test adjVet(net) == [1, 2, 3, 4, 6, 7, 8, 9, 11, 12, 13, 14]
		@test adjMat(net) == BitArray([
			0 1 1 1;
			1 0 1 1;
			1 1 0 1;
			1 1 1 0
		])
	end

	@testset "Regular" begin
		# Undirected
		net = RegularNetwork(10, 4)
	    @test numnodes(net) == 10
	    @test numconnections(net) == 20
	    @test !isdirected(net)
	    @test connectivity(net, 1) == 4
		@test connectivities(net) == [4, 4, 4, 4, 4, 4, 4, 4, 4, 4]
	    @test connectivity(net) == 4

		# Directed
		net = RegularNetwork(10, 4, directed=true)
	    @test numconnections(net) == 40
	    @test isdirected(net)

	    net = RegularNetwork(5, 2)
	    @test adjMat(net) == [
	    	0 1 0 0 1;
	    	1 0 1 0 0;
	    	0 1 0 1 0;
	    	0 0 1 0 1;
	    	1 0 0 1 0
	    ]
	    @test adjMat(net, sparse=true) == sparse(
	    	[2, 5, 1, 3, 2, 4, 3, 5, 1, 4],
	    	[1, 1, 2, 2, 3, 3, 4, 4, 5, 5],
	    	1,
	    )
	    @test adjVet(net) == [1, 4, 5, 7, 11, 13, 17, 19, 20, 23]

	    @test adjMat(RegularNetwork(10, 2)) == adjMat(RegularNetwork(10, 2, directed=true))

		# Argument manip
		net = RegularNetwork(10, 5)
		@test connectivity(net) == 4
		@test numconnections(net) == 4 * 10 / 2
	    @test adjMat(net)  == adjMat(RegularNetwork(10, 4))
		net = RegularNetwork(10, -5)
		@test connectivity(net) == 0
		@test numconnections(net) == 0
	    @test adjMat(net) == adjMat(EmptyNetwork(10))
	    net = RegularNetwork(10, 20)
		@test connectivity(net) == 9
		@test adjMat(net) == adjMat(GlobalNetwork(10))
		@test numconnections(net) == sum(adjMat(net)) / 2
		net = RegularNetwork(11, 20)
		@test connectivity(net) == 10
		@test adjMat(net) == adjMat(GlobalNetwork(11))
		@test numconnections(net) == sum(adjMat(net)) / 2
	end

	@testset "ErdosRenyi p-initialized" begin
		# Undirected Network
		net = ErdosRenyiNetwork(500, 0.5, seed=1234)
		@test numnodes(net) == 500
		@test net.p == 0.5
		@test numconnections(net) ≈ (net.p * numnodes(net) * (numnodes(net) - 1) / 2) rtol=0.1
		@test !isdirected(net)
		@test net.seed == 1234

		@test AdjacencyArrays.issymmetric(adjMat(net))
		@test sum(adjMat(net)) == 2 * numconnections(net)
		@test length(adjVet(net)) == 2 * numconnections(net)

		# Directed Network
		net = ErdosRenyiNetwork(500, 0.5, seed=1234, directed=true)
		@test isdirected(net)
		@test numconnections(net) ≈ (net.p * numnodes(net) * (numnodes(net) - 1)) rtol=0.1
		@test sum(adjMat(net)) == numconnections(net)
		@test length(adjVet(net)) == numconnections(net)

		# Seed Behaviour
		net = ErdosRenyiNetwork(10, 0.5)
		@test adjMat(ErdosRenyiNetwork(10, 0.5, seed=net.seed)) == adjMat(net)

		# Argument manip
		@test ErdosRenyiNetwork(10, -1.).p == 0.0
		@test ErdosRenyiNetwork(10, 1.2).p == 1.0
		@test adjMat(ErdosRenyiNetwork(10, -1.)) == adjMat(EmptyNetwork(10))
		@test adjMat(ErdosRenyiNetwork(10, 1.2)) == adjMat(GlobalNetwork(10))
	end

	@testset "ErdosRenyi numConnections-initialized" begin
		# Undirected Network
		net = ErdosRenyiNetwork(10, 20, seed=1234)
		@test numnodes(net) == 10
		@test numconnections(net) == 20
		@test isnothing(net.p)
		@test !isdirected(net)
		@test net.seed == 1234

		@test AdjacencyArrays.issymmetric(adjMat(net))
		@test sum(adjMat(net)) == 2 * numconnections(net)
		@test length(adjVet(net)) == 2 * numconnections(net)

		# Directed Network
		net = ErdosRenyiNetwork(10, 20, seed=1234, directed=true)
		@test numconnections(net) == 20
		@test isnothing(net.p)
		@test isdirected(net)
		@test sum(adjMat(net)) == numconnections(net)
		@test length(adjVet(net)) == numconnections(net)

		# Seed Behaviour
		net = ErdosRenyiNetwork(10, 30)
		@test adjMat(ErdosRenyiNetwork(10, 30, seed=net.seed)) == adjMat(net)

		# Argument manip
		@test numconnections(ErdosRenyiNetwork(10, -1)) == 0
		@test numconnections(ErdosRenyiNetwork(10, 100)) == 45
		@test adjMat(ErdosRenyiNetwork(10, -1))  == adjMat(EmptyNetwork(10))
		@test adjMat(ErdosRenyiNetwork(10, 100)) == adjMat(GlobalNetwork(10))
	end

	@testset "WattsStrogatz β-initialized" begin
	    # Undirected Network
		net = WattsStrogatzNetwork(500, 10, 0.2, seed=1234)
		@test numnodes(net) == 500
		@test connectivity(net) == 10
		@test net.β == 0.2
		@test numconnections(net) == 10 * 500 / 2
		@test numshortcuts(net) ≈ (net.β * numconnections(net)) rtol=0.1
		@test !isdirected(net)
		@test net.seed == 1234

		@test AdjacencyArrays.issymmetric(adjMat(net))
		@test sum(adjMat(net)) == 2 * numconnections(net)
		@test length(adjVet(net)) == 2 * numconnections(net)

		# Directed Network
		net = WattsStrogatzNetwork(500, 10, 0.2, seed=1234, directed=true)
		@test numconnections(net) == 10 * 500
		@test numshortcuts(net) ≈ (net.β * numconnections(net)) rtol=0.1
		@test isdirected(net)
		@test sum(adjMat(net)) == numconnections(net)
		@test length(adjVet(net)) == numconnections(net)

	    # Seed Behaviour
		net = WattsStrogatzNetwork(20, 4, 0.2)
		@test adjMat(WattsStrogatzNetwork(20, 4, 0.2, seed=net.seed)) == adjMat(net)

		# Argument manip
		net = WattsStrogatzNetwork(10, -1, 0.2)
		@test connectivity(net) == 0
		@test numconnections(net) == 0
		@test numshortcuts(net) == 0
		@test isnothing(net.β)
		@test adjMat(net) == adjMat(EmptyNetwork(10))

		net = WattsStrogatzNetwork(10, 11, 0.2)
		@test connectivity(net) == 9
		@test numconnections(net) == 90 / 2
		@test numshortcuts(net) == 0
		@test isnothing(net.β)
		@test adjMat(net) == adjMat(GlobalNetwork(10))

		net = WattsStrogatzNetwork(10, 2, -0.5)
		@test net.β == 0.0
		@test numshortcuts(net) == 0
		@test adjMat(net) == adjMat(RegularNetwork(10, 2))

		net = WattsStrogatzNetwork(10, 2, 1.5)
		@test net.β == 1.0
		@test numshortcuts(net) == numconnections(net)
	end

	@testset "WattsStrogatz numShortcuts-initialized" begin
	    # Undirected Network
		net = WattsStrogatzNetwork(500, 10, 100, seed=1234)
		@test numnodes(net) == 500
		@test connectivity(net) == 10
		@test isnothing(net.β)
		@test numconnections(net) == 10 * 500 / 2
		@test numshortcuts(net) == 100
		@test !isdirected(net)
		@test net.seed == 1234

		@test AdjacencyArrays.issymmetric(adjMat(net))
		@test sum(adjMat(net)) == 2 * numconnections(net)
		@test length(adjVet(net)) == 2 * numconnections(net)

		# Directed Network
		net = WattsStrogatzNetwork(500, 10, 100, seed=1234, directed=true)
		@test numconnections(net) == 10 * 500
		@test numshortcuts(net) == 100
		@test isdirected(net)
		@test sum(adjMat(net)) == numconnections(net)
		@test length(adjVet(net)) == numconnections(net)

	    # Seed Behaviour
		net = WattsStrogatzNetwork(20, 4, 10)
		@test adjMat(WattsStrogatzNetwork(20, 4, 10, seed=net.seed)) == adjMat(net)

		# Argument manip
		net = WattsStrogatzNetwork(10, -1, 5)
		@test connectivity(net) == 0
		@test numconnections(net) == 0
		@test numshortcuts(net) == 0
		@test isnothing(net.β)
		@test adjMat(net) == adjMat(EmptyNetwork(10))

		net = WattsStrogatzNetwork(10, 11, 5)
		@test connectivity(net) == 9
		@test numconnections(net) == 90 / 2
		@test numshortcuts(net) == 0
		@test isnothing(net.β)
		@test adjMat(net) == adjMat(GlobalNetwork(10))

		net = WattsStrogatzNetwork(10, 2, -5)
		@test numshortcuts(net) == 0
		@test adjMat(net) == adjMat(RegularNetwork(10, 2))

		net = WattsStrogatzNetwork(10, 2, 40)
		@test numshortcuts(net) == numconnections(net)
	end

	@testset "Custom Networks" begin
		mat = [
			0 1 1 0;
			0 0 1 1;
			1 0 0 1;
			1 1 0 0
		]

	    @test adjMat(CustomNetwork(mat)) == mat
	    @test adjMat(CustomNetwork(adjMatToVet(mat), size(mat, 1))) == mat
	    @test adjMat(CustomNetwork(Float64.(mat))) == mat
	    @test adjMat(CustomNetwork(Bool.(mat))) == mat
	    @test adjMat(CustomNetwork(sparse(mat))) == mat

	    @test isdirected(CustomNetwork(mat))
	    @test !isdirected(CustomNetwork([0 1 0; 1 0 1; 0 1 0]))

	    @test_throws ArgumentError CustomNetwork([0 1 1; 1 0 1])
	    @test_throws ArgumentError CustomNetwork([0 2 1; 1 0 1; 3 2 0])
	    @test_throws ArgumentError CustomNetwork([0 0.5 1; 0.5 0 1; 0.75 0.2 0])
	    @test_throws ArgumentError CustomNetwork(mat, directed=false)
	end

	@testset "Trivial Networks" begin
		# N = 1
		trivialAdjMat = BitMatrix([0 for i in [1], j in [1]])

		@test adjMat(EmptyNetwork(1)) == trivialAdjMat
		@test adjMat(GlobalNetwork(1)) == trivialAdjMat
		@test adjMat(RegularNetwork(1, 2)) == trivialAdjMat
		@test adjMat(ErdosRenyiNetwork(1, 0.4)) == trivialAdjMat
		@test adjMat(ErdosRenyiNetwork(1, 4)) == trivialAdjMat
		@test adjMat(WattsStrogatzNetwork(1, 2, 0.4)) == trivialAdjMat
		@test adjMat(WattsStrogatzNetwork(1, 2, 4)) == trivialAdjMat

		@test adjVet(EmptyNetwork(1)) == []
		@test adjVet(GlobalNetwork(1)) == []
		@test adjVet(RegularNetwork(1, 2)) == []
		@test adjVet(ErdosRenyiNetwork(1, 0.4)) == []
		@test adjVet(ErdosRenyiNetwork(1, 4)) == []
		@test adjVet(WattsStrogatzNetwork(1, 2, 0.4)) == []
		@test adjVet(WattsStrogatzNetwork(1, 2, 4)) == []

		# N < 1
		@test_throws ArgumentError EmptyNetwork(0)
		@test_throws ArgumentError EmptyNetwork(-2)
		@test_throws ArgumentError GlobalNetwork(0)
		@test_throws ArgumentError GlobalNetwork(-2)
		@test_throws ArgumentError RegularNetwork(0, 2)
		@test_throws ArgumentError RegularNetwork(-2, 2)
		@test_throws ArgumentError ErdosRenyiNetwork(0, 0.4)
		@test_throws ArgumentError ErdosRenyiNetwork(-2, 0.4)
		@test_throws ArgumentError ErdosRenyiNetwork(0, 4)
		@test_throws ArgumentError ErdosRenyiNetwork(-2, 4)
		@test_throws ArgumentError WattsStrogatzNetwork(0, 2, 0.4)
		@test_throws ArgumentError WattsStrogatzNetwork(-2, 2, 0.4)
		@test_throws ArgumentError WattsStrogatzNetwork(0, 2, 4)
		@test_throws ArgumentError WattsStrogatzNetwork(-2, 2, 4)
	end

	@testset "Properties" begin
		# Connectivity
		net = GlobalNetwork(10)
	    @test_throws ArgumentError connectivity(net, 0)
	    @test_throws ArgumentError connectivity(net, 11)

	    net = RegularNetwork(10, 4, directed=true)
	    @test connectivity(net, 1, dir_behaviour=:in) == 4
	    @test connectivity(net, 1, dir_behaviour=:out) == 4
	    @test connectivity(net, 1, dir_behaviour=:total) == 2 * 4
	    @test connectivity(net, 1, dir_behaviour=:mean) == 4
	    @test connectivity(net, 1, dir_behaviour=:both) == (4, 4)

	    # Clustering Coefficient
		net = GlobalNetwork(10)
	    @test clusteringcoefficient(net, 1) == 1
	    @test clusteringcoefficient(net) == 1
	    @test clusteringcoefficients(net) == repeat([1], 10)

		net = RegularNetwork(20, 4)
	    @test clusteringcoefficient(net, 1) == 0.5
	    @test clusteringcoefficient(net) == 0.5
	    @test clusteringcoefficients(net) == repeat([0.5], 20)

		net = ErdosRenyiNetwork(300, 0.2, seed=1234)
	    @test clusteringcoefficient(net) ≈ 0.2 rtol=0.1

	    net = GlobalNetwork(10, directed=true)
	    @test clusteringcoefficient(net, 1) == 1
	    @test clusteringcoefficient(net) == 1
	    @test clusteringcoefficients(net) == repeat([1], 10)

		net = RegularNetwork(20, 4, directed=true)
	    @test clusteringcoefficient(net, 1) == 0.5
	    @test clusteringcoefficient(net) == 0.5
	    @test clusteringcoefficients(net) == repeat([0.5], 20)

	    # Shortest path length
	    @test shortestpath(GlobalNetwork(8), 1) == [0, 1, 1, 1, 1, 1, 1, 1]
	    @test shortestpath(GlobalNetwork(8)) == 1
	    @test shortestpath(RegularNetwork(8, 2), 1) == [0, 1, 2, 3, 4, 3, 2, 1]
	    @test shortestpath(RegularNetwork(8, 2), 1, 4) == 3
	    @test shortestpath(RegularNetwork(8, 2)) == 16 / 7

	    net = CustomNetwork([
	    	0 0 0 0 0;
	    	1 0 1 0 1;
	    	1 0 0 1 0;
	    	0 1 0 0 1;
	    	0 0 0 0 0
	    ])
	    @test shortestpath(net, 1) == [0, 1, 1, 2, Inf]
	    @test shortestpath(net, 1, 4) == 2
	    @test shortestpath(net) == 17 / (5 * 4)
	    @test averagepathlength(net) == shortestpath(net)

	    @test shortestpath(EmptyNetwork(10)) == 0

	    # Small-world-ness
	    net = WattsStrogatzNetwork(400, 4, 0.2; seed=1252, directed=true)
	    equivRan = equivalentRandomNetwork(net)
	    equivLat = equivalentLatticeNetwork(net)
	    @test connectivity(net, dir_behaviour=:total) == connectivity(equivRan, dir_behaviour=:total)
	    @test connectivity(net, dir_behaviour=:in) == connectivity(equivRan, dir_behaviour=:in)
	    @test clusteringcoefficient(net) < clusteringcoefficient(equivLat)
	end

	@testset "Operators" begin
		# Undirected netowrks
		@test allEdges(EmptyNetwork(5)) == []
	    @test allEdges(GlobalNetwork(4)) == [
	    	(1, 2),
	    	(1, 3),
	    	(1, 4),
	    	(2, 3),
	    	(2, 4),
	    	(3, 4)
	    ]
	    @test allEdges(RegularNetwork(5, 2)) == [
	    	(1, 2),
	    	(1, 5),
	    	(2, 3),
	    	(3, 4),
	    	(4, 5)
	    ]

	    # Directed netowrks
	    @test allEdges(GlobalNetwork(4, directed=true)) == [
	    	(1, 2),
	    	(1, 3),
	    	(1, 4),
	    	(2, 1),
	    	(2, 3),
	    	(2, 4),
	    	(3, 1),
	    	(3, 2),
	    	(3, 4),
	    	(4, 1),
	    	(4, 2),
	    	(4, 3)
	    ]
	    @test allEdges(RegularNetwork(5, 2, directed=true)) == [
	    	(1, 2),
	    	(1, 5),
	    	(2, 1),
	    	(2, 3),
	    	(3, 2),
	    	(3, 4),
	    	(4, 3),
	    	(4, 5),
	    	(5, 1),
	    	(5, 4)
	    ]

	    # Different first_index
	    @test allEdges(RegularNetwork(5, 2), first_index=0) == [
	    	(0, 1),
	    	(0, 4),
	    	(1, 2),
	    	(2, 3),
	    	(3, 4)
	    ]

	    @test_throws ArgumentError allEdges(GlobalNetwork(3), first_index=2)
	end

	@testset "Converters" begin
		net = RegularNetwork(10, 4)
		vet = adjVet(net)
		mat = adjMat(net)

		@test adjVetToMat(vet, 10) == mat
		@test adjMatToVet(mat) == vet
		@test adjMatToVet(Int.(mat)) == vet
		@test adjMatToVet(Bool.(mat)) == vet

		@test_throws ArgumentError adjMatToVet([1 0; 3 2])
		@test_throws AssertionError adjMatToVet([1 0 1; 0 0 1])
	end
end

tests()


text:
struct SARAH_prox_DNN_iterable{R<:Real,Tf,Tg,Tp}
    F::Union{Array{Tf},Tf}            # smooth term - array of f_i s
    g::Tg                   # nonsmooth part
    opt_params::Tp                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of ∇f_i
    μ::Maybe{Union{Array{R},R}}  # convexity moduli of the gradients
    γ::Maybe{R}             # stepsize (a single scalar)
    m::Maybe{Int}           # number of inner loop updates
    ꞵ::Maybe{R}             # a scalar momentum
    data::Tuple
    DNN_config!
end

mutable struct SARAH_prox_DNN_state{R<:Real,Tx}
    γ::R                    # stepsize
    m::Maybe{Int}           # number of inner loop updates
    ꞵ::R                    # momentum
    av::Tx                  # average of ∇f_i s in outer loop
    w::Tx                   # w_{t-1} in the inner loop
    w_plus::Tx              # w_{t} in the inner loop
    ind::Array{Int}         # running idx set (1 to N)

    # some extra placeholders
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx                # placeholder for gradients
end

function SARAH_prox_DNN_state(γ::R, m, ꞵ::R, av::Tx, w::Tx, w_plus::Tx, ind) where {R,Tx}
    return SARAH_prox_DNN_state{R,Tx}(γ, m, ꞵ, av, w, w_plus, ind, copy(av), copy(av))
end

batchmemaybe(x) = tuple(x)
batchmemaybe(x::Tuple) = x

function Base.iterate(iter::SARAH_prox_DNN_iterable{R}) where {R}
    N = iter.N
    m = iter.m
    batch = 1
    ω = (3*(N-batch))/(2*batch*(N-1))
    ind = collect(1:N) # a vector of 1 to N
    # updating the stepsize
    if iter.γ === nothing
        γ = (2*sqrt(ω * m)) / (4*sqrt(ω * m) + 1)
    else
        γ = iter.γ # provided γ
    end

    if iter.ꞵ === nothing
        if iter.L === nothing
            @warn "smoothness parameter absent"
            return nothing
        else
            L_M = maximum(iter.L)
            ꞵ = 1/(L_M * sqrt(ω * m))
        end
    else
        ꞵ = iter.ꞵ
    end

    # initializing the vectors
    x0 = iter.DNN_config!()
    # av = zero(x0)
    w_plus = zero(x0)
    # for i = 1:N
    #     ∇f = (Flux.gradient(() -> iter.F(i,iter.x0), params(iter.x0)))[iter.x0][:,1]
    #     # ∇f, ~ = gradient(iter.F[i], iter.x0)
    #     ∇f ./= N
    #     av .+= ∇f
    # end
    gs = Flux.gradient(iter.opt_params) do # sampled gradient
        iter.F(batchmemaybe(iter.data)...)
    end
    av = iter.DNN_config!(gs=gs)
    w = copy(x0)
    CIAOAlgorithms.prox!(w_plus, iter.g, w - γ * av, γ)
    w_plus .*= ꞵ
    w_plus .+= (1 - ꞵ) * w

    state = SARAH_prox_DNN_state(γ, m, ꞵ, av, w, w_plus, ind)
    return state, state
end

function Base.iterate(iter::SARAH_prox_DNN_iterable{R}, state::SARAH_prox_DNN_state{R}) where {R}
    # The inner cycle
    for i in rand(state.ind, state.m)
        temp_x = iter.data[1][:,i]
        temp_y = iter.data[2][:,i]
        # gradient!(state.temp, iter.F[i], state.w_plus)
        # state.temp .= (Flux.gradient(() -> iter.F(i,state.w_plus), params(state.w_plus)))[state.w_plus][:,1]
        iter.DNN_config!(gs=state.w_plus)
        gs = Flux.gradient(iter.opt_params) do # sampled gradient
            iter.F(temp_x,temp_y)
        end
        state.temp .= iter.DNN_config!(gs=gs)

        # gradient!(state.∇f_temp, iter.F[i], state.w)
        # state.∇f_temp .= (Flux.gradient(() -> iter.F(i,state.w), params(state.w)))[state.w][:,1]
        iter.DNN_config!(gs=state.w)
        gs = Flux.gradient(iter.opt_params) do # sampled gradient
            iter.F(temp_x,temp_y)
        end
        state.∇f_temp .= iter.DNN_config!(gs=gs)

        state.av .+= state.temp
        state.av .-= state.∇f_temp

        state.temp .= state.av
        state.temp .*= state.γ
        state.temp .*= -1
        state.temp .+= state.w_plus
        state.w .= state.w_plus
        CIAOAlgorithms.prox!(state.w_plus, iter.g, state.temp, state.γ)
        state.w_plus .*= state.ꞵ
        state.temp .= state.w
        state.temp .*= 1-state.ꞵ
        state.w_plus .+= state.temp
    end
    # full update
    state.w .= state.w_plus
    # state.av .= zero(state.w)  # for next iterate
    # for i = 1:iter.N
    #     # gradient!(state.∇f_temp, iter.F[i], state.w)
    #     state.∇f_temp .= (Flux.gradient(() -> iter.F(i,state.w), params(state.w)))[state.w][:,1]
    #     state.∇f_temp ./= iter.N
    #     state.av .+= state.∇f_temp
    # end
    iter.DNN_config!(gs=state.w)
    gs = Flux.gradient(iter.opt_params) do # sampled gradient
        iter.F(batchmemaybe(iter.data)...)
    end
    state.av .= iter.DNN_config!(gs=gs)

    state.temp .= state.av
    state.temp .*= state.γ
    state.temp .*= -1
    state.temp .+= state.w
    CIAOAlgorithms.prox!(state.w_plus, iter.g, state.temp, state.γ)
    state.w_plus .*= state.ꞵ
    state.temp .= state.w
    state.temp .*= 1-state.ꞵ
    state.w_plus .+= state.temp

    return state, state
end

solution(state::SARAH_prox_DNN_state) = state.w_plus


text:
"""
    VTKPolyhedron <: AbstractMeshCell

Represents a polyhedron cell in an unstructured grid.

Using `VTKPolyhedron` should be preferred to using a `MeshCell` with a cell type
`VTKCellTypes.VTK_POLYHEDRON`, since the latter cannot hold all the necessary
information to describe a polyhedron cell.

---

    VTKPolyhedron(connectivity, faces...)

Construct polyhedron cell from connectivity vector (see [`MeshCell`](@ref) for
details) and from a list of polyhedron faces.

# Example

Create a polyhedron with 8 points and 6 faces.
This can represent a cube if the 8 points are properly positioned.

```jldoctest
julia> cell = VTKPolyhedron(
           1:8,
           (1, 4, 3, 2),
           (1, 5, 8, 4),
           (5, 6, 7, 8),
           (6, 2, 3, 7),
           (1, 2, 6, 5),
           (3, 4, 8, 7),
       )
VTKPolyhedron{UnitRange{Int64}, NTuple{6, NTuple{4, Int64}}}(1:8, ((1, 4, 3, 2), (1, 5, 8, 4), (5, 6, 7, 8), (6, 2, 3, 7), (1, 2, 6, 5), (3, 4, 8, 7)))
```
"""
struct VTKPolyhedron{V <: Connectivity, Faces} <: AbstractMeshCell
    connectivity :: V
    faces :: Faces

    function VTKPolyhedron(connectivity, faces...)
        new{typeof(connectivity), typeof(faces)}(connectivity, faces)
    end
end

Base.eltype(::Type{<:VTKPolyhedron}) = VTKCellTypes.VTKCellType
cell_type(::VTKPolyhedron) = VTKCellTypes.VTK_POLYHEDRON

faces(cell::VTKPolyhedron) = cell.faces

function process_faces!(data, cell::VTKPolyhedron, offset)
    fs = faces(cell)
    num_faces = length(fs)
    num_values = 1 + sum(x -> length(x) + 1, fs)
    resize!(data, length(data) + num_values)
    @inbounds data[offset += 1] = num_faces
    @inbounds for f in fs
        data[offset += 1] = length(f)
        for idx in f
            data[offset += 1] = idx - 1  # switch to zero-based indexing
        end
    end
    num_values
end


text:
using SnpArrays
using Knockoffs
using DataFrames
using CSV

res = 0
for chr in 1:22
    chr == 10 && continue
    plinkfile = "/scratch/users/bbchu/ukb_SHAPEIT/knockoffs/ukb_gen_chr$(chr)_ibd1_res$(res)"
    xdata = SnpData(plinkfile)
    isknockoff = endswith.(xdata.snp_info[!, :snpid], ".k")
    original, knockoff = Int[], Int[]
    for i in 1:size(xdata)[2]
        isknockoff[i] ? push!(knockoff, i) : push!(original, i)
    end

    @time decorrelate_knockoffs(plinkfile, original, knockoff;
        outfile = "ukb_gen_chr$(chr)_ibd1_res$(res)_decorrelated",
        outdir = pwd(),
        r2_threshold = 0.95)
end


text:
if Base.VERSION >= v"0.7"
    using Pkg
end

println("Pkg.free(... all packages ...)")
Pkg.free(collect(Iterators.filter(x -> x != "LyapunovExponents",
                                  keys(Pkg.installed()))))
println("Pkg.resolve()")
Pkg.resolve()
println("Pkg.update()")
Pkg.update()


text:
module TableOperations

using Tables, SentinelArrays

struct TransformsRow{T, F} <: Tables.AbstractRow
    row::T
    funcs::F
end

getrow(r::TransformsRow) = getfield(r, :row)
getfuncs(r::TransformsRow) = getfield(r, :funcs)

Tables.getcolumn(row::TransformsRow, nm::Symbol) = (getfunc(row, getfuncs(row), nm))(Tables.getcolumn(getrow(row), nm))
Tables.getcolumn(row::TransformsRow, i::Int) = (getfunc(row, getfuncs(row), i))(Tables.getcolumn(getrow(row), i))
Tables.columnnames(row::TransformsRow) = Tables.columnnames(getrow(row))

struct Transforms{C, T, F}
    source::T
    funcs::F # NamedTuple of columnname=>transform function
end

Tables.columnnames(t::Transforms{true}) = Tables.columnnames(getfield(t, 1))
Tables.getcolumn(t::Transforms{true}, nm::Symbol) = Base.map(getfunc(t, getfield(t, 2), nm), Tables.getcolumn(getfield(t, 1), nm))
Tables.getcolumn(t::Transforms{true}, i::Int) = Base.map(getfunc(t, getfield(t, 2), i), Tables.getcolumn(getfield(t, 1), i))
# for backwards compat
Base.propertynames(t::Transforms{true}) = Tables.columnnames(t)
Base.getproperty(t::Transforms{true}, nm::Symbol) = Tables.getcolumn(t, nm)

"""
    Tables.transform(source, funcs) => Tables.Transforms
    source |> Tables.transform(funcs) => Tables.Transform

Given any Tables.jl-compatible source, apply a series of transformation functions, for the columns specified in `funcs`.
The tranform functions can be a NamedTuple or Dict mapping column name (`String` or `Symbol` or `Integer` index) to Function.
"""
function transform end

transform(funcs) = x->transform(x, funcs)
transform(; kw...) = transform(kw.data)
function transform(src::T, funcs::F) where {T, F}
    C = Tables.columnaccess(T)
    x = C ? Tables.columns(src) : Tables.rows(src)
    return Transforms{C, typeof(x), F}(x, funcs)
end

getfunc(row, nt::NamedTuple, nm::Symbol) = get(nt, nm, identity)
getfunc(row, d::Dict{String, <:Base.Callable}, nm::Symbol) = get(d, String(nm), identity)
getfunc(row, d::Dict{Symbol, <:Base.Callable}, nm::Symbol) = get(d, nm, identity)
getfunc(row, d::Dict{Int, <:Base.Callable}, nm::Symbol) = get(d, findfirst(isequal(nm), Tables.columnnames(row)), identity)

getfunc(row, nt::NamedTuple, i::Int) = get(nt, Tables.columnnames(row)[i], identity)
getfunc(row, d::Dict{String, <:Base.Callable}, i::Int) = get(d, String(Tables.columnnames(row)[i]), identity)
getfunc(row, d::Dict{Symbol, <:Base.Callable}, i::Int) = get(d, Tables.columnnames(row)[i], identity)
getfunc(row, d::Dict{Int, <:Base.Callable}, i::Int) = get(d, i, identity)

Tables.istable(::Type{<:Transforms}) = true
Tables.rowaccess(::Type{Transforms{C, T, F}}) where {C, T, F} = !C
Tables.rows(t::Transforms{false}) = t
Tables.columnaccess(::Type{Transforms{C, T, F}}) where {C, T, F} = C
Tables.columns(t::Transforms{true}) = t
# avoid relying on inference here and just let sinks figure things out
Tables.schema(t::Transforms) = nothing

Base.IteratorSize(::Type{Transforms{false, T, F}}) where {T, F} = Base.IteratorSize(T)
Base.length(t::Transforms{false}) = length(getfield(t, 1))
Base.eltype(t::Transforms{false, T, F}) where {T, F} = TransformsRow{eltype(getfield(t, 1)), F}

@inline function Base.iterate(t::Transforms{false}, st=())
    state = iterate(getfield(t, 1), st...)
    state === nothing && return nothing
    return TransformsRow(state[1], getfield(t, 2)), (state[2],)
end

# select
struct Select{T, columnaccess, names}
    source::T
end

"""
    Tables.select(source, columns...) => Tables.Select
    source |> Tables.select(columns...) => Tables.Select

Create a lazy wrapper that satisfies the Tables.jl interface and keeps only the columns given by the columns arguments, which can be `String`s, `Symbol`s, or `Integer`s
"""
function select end

select(names::Symbol...) = x->select(x, names...)
select(names::String...) = x->select(x, Base.map(Symbol, names)...)
select(inds::Integer...) = x->select(x, Base.map(Int, inds)...)

function select(x::T, names...) where {T}
    colaccess = Tables.columnaccess(T)
    r = colaccess ? Tables.columns(x) : Tables.rows(x)
    return Select{typeof(r), colaccess, names}(r)
end

Tables.istable(::Type{<:Select}) = true

Base.@pure function typesubset(::Tables.Schema{names, types}, nms::NTuple{N, Symbol}) where {names, types, N}
    return Tuple{Any[Tables.columntype(names, types, nm) for nm in nms]...}
end

Base.@pure function typesubset(::Tables.Schema{names, types}, inds::NTuple{N, Int}) where {names, types, N}
    return Tuple{Any[fieldtype(types, i) for i in inds]...}
end

typesubset(::Tables.Schema{names, types}, ::Tuple{}) where {names, types} = Tuple{}

namesubset(::Tables.Schema{names, types}, nms::NTuple{N, Symbol}) where {names, types, N} = nms
Base.@pure namesubset(::Tables.Schema{names, T}, inds::NTuple{N, Int}) where {names, T, N} = ntuple(i -> names[inds[i]], N)
namesubset(::Tables.Schema{names, types}, ::Tuple{}) where {names, types} = ()
namesubset(names, nms::NTuple{N, Symbol}) where {N} = nms
namesubset(names, inds::NTuple{N, Int}) where {N} = ntuple(i -> names[inds[i]], N)
namesubset(names, ::Tuple{}) = ()

function Tables.schema(s::Select{T, columnaccess, names}) where {T, columnaccess, names}
    sch = Tables.schema(getfield(s, 1))
    sch === nothing && return nothing
    return Tables.Schema(namesubset(sch, names), typesubset(sch, names))
end

# columns: make Select property-accessible
Tables.getcolumn(s::Select{T, true, names}, nm::Symbol) where {T, names} = Tables.getcolumn(getfield(s, 1), nm)
Tables.getcolumn(s::Select{T, true, names}, i::Int) where {T, names} = Tables.getcolumn(getfield(s, 1), i)
Tables.columnnames(s::Select{T, true, names}) where {T, names} = namesubset(Tables.columnnames(getfield(s, 1)), names)
Tables.columnaccess(::Type{Select{T, C, names}}) where {T, C, names} = C
Tables.columns(s::Select{T, true, names}) where {T, names} = s
# for backwards compat
Base.propertynames(s::Select{T, true, names}) where {T, names} = Tables.columnnames(s)
Base.getproperty(s::Select{T, true, names}, nm::Symbol) where {T, names} = Tables.getcolumn(s, nm)

# rows: implement Iterator interface
Base.IteratorSize(::Type{Select{T, false, names}}) where {T, names} = Base.IteratorSize(T)
Base.length(s::Select{T, false, names}) where {T, names} = length(getfield(s, 1))
Base.IteratorEltype(::Type{Select{T, false, names}}) where {T, names} = Base.IteratorEltype(T)
Base.eltype(s::Select{T, false, names}) where {T, names} = SelectRow{eltype(getfield(s, 1)), names}
Tables.rowaccess(::Type{Select{T, columnaccess, names}}) where {T, columnaccess, names} = !columnaccess
Tables.rows(s::Select{T, false, names}) where {T, names} = s

# we need to iterate a "row view" in case the underlying source has unknown schema
# to ensure each iterated row only has `names` Tables.columnnames
struct SelectRow{T, names} <: Tables.AbstractRow
    row::T
end

Tables.getcolumn(row::SelectRow, nm::Symbol) = Tables.getcolumn(getfield(row, 1), nm)
Tables.getcolumn(row::SelectRow{T, names}, i::Int) where {T, names} = Tables.getcolumn(getfield(row, 1), names[i])
Tables.getcolumn(row::SelectRow, ::Type{T}, i::Int, nm::Symbol) where {T} = Tables.getcolumn(getfield(row, 1), T, Tables.columnindex(Tables.columnnames(getfield(row, 1)), nm), nm)

getprops(row, nms::NTuple{N, Symbol}) where {N} = nms
getprops(row, inds::NTuple{N, Int}) where {N} = ntuple(i->Tables.columnnames(getfield(row, 1))[inds[i]], N)
getprops(row, ::Tuple{}) = ()

Tables.columnnames(row::SelectRow{T, names}) where {T, names} = getprops(row, names)

@inline function Base.iterate(s::Select{T, false, names}) where {T, names}
    state = iterate(getfield(s, 1))
    state === nothing && return nothing
    row, st = state
    return SelectRow{typeof(row), names}(row), st
end

@inline function Base.iterate(s::Select{T, false, names}, st) where {T, names}
    state = iterate(getfield(s, 1), st)
    state === nothing && return nothing
    row, st = state
    return SelectRow{typeof(row), names}(row), st
end

# filter
struct Filter{F, T}
    f::F
    x::T
end

"""
    Tables.filter(f, source) => Tables.Filter
    source |> Tables.filter(f) => Tables.Filter

Create a lazy wrapper that satisfies the Tables.jl interface and keeps the rows where `f(row)` is true.
"""
function filter end

function filter(f::F, x) where {F <: Base.Callable}
    r = Tables.rows(x)
    return Filter{F, typeof(r)}(f, r)
end
filter(f::Base.Callable) = x->filter(f, x)

Tables.isrowtable(::Type{<:Filter}) = true
Tables.schema(f::Filter) = Tables.schema(f.x)

Base.IteratorSize(::Type{Filter{F, T}}) where {F, T} = Base.SizeUnknown()
Base.IteratorEltype(::Type{Filter{F, T}}) where {F, T} = Base.IteratorEltype(T)
Base.eltype(f::Filter) = eltype(f.x)

 @inline function Base.iterate(f::Filter)
    state = iterate(f.x)
    state === nothing && return nothing
    while !f.f(state[1])
        state = iterate(f.x, state[2])
        state === nothing && return nothing
    end
    return state
end

 @inline function Base.iterate(f::Filter, st)
    state = iterate(f.x, st)
    state === nothing && return nothing
    while !f.f(state[1])
        state = iterate(f.x, state[2])
        state === nothing && return nothing
    end
    return state
end

# map
struct Map{F, T}
    func::F
    source::T
end

"""
    Tables.map(f, source) => Tables.Map
    source |> Tables.map(f) => Tables.Map

Create a lazy wrapper that satisfies the Tables.jl interface and will apply the function `f(row)` to each
row in the input table source. Note that `f` must take and produce a valid Tables.jl `Row` object.
"""
function map end

function map(f::F, x::T) where {F <: Base.Callable, T}
    r = Tables.rows(x)
    return Map{F, typeof(r)}(f, r)
end
map(f::Base.Callable) = x->map(f, x)

Tables.isrowtable(::Type{<:Map}) = true
Tables.schema(m::Map) = nothing

Base.IteratorSize(::Type{Map{T, F}}) where {T, F} = Base.IteratorSize(T)
Base.length(m::Map) = length(m.source)
Base.IteratorEltype(::Type{<:Map}) = Base.EltypeUnknown()

@inline function Base.iterate(m::Map)
    state = iterate(m.source)
    state === nothing && return nothing
    return m.func(state[1]), state[2]
end

@inline function Base.iterate(m::Map, st)
    state = iterate(m.source, st)
    state === nothing && return nothing
    return m.func(state[1]), state[2]
end

# joinpartitions
struct JoinedPartitions{S} <: Tables.AbstractColumns
    schema::S
    x::Vector{ChainedVector}
    lookup::Dict{Symbol, ChainedVector}
end

Tables.istable(::Type{<:JoinedPartitions}) = true
Tables.columnaccess(::Type{<:JoinedPartitions}) = true
Tables.columns(x::JoinedPartitions) = Tables.CopiedColumns(x)
Tables.columnnames(x::JoinedPartitions) = Tables.schema(x).names
Tables.getcolumn(x::JoinedPartitions, i::Int) = getfield(x, :x)[i]
Tables.getcolumn(x::JoinedPartitions, nm::Symbol) = getfield(x, :lookup)[nm]
Tables.schema(x::JoinedPartitions) = getfield(x, :schema)

"""
    TableOperations.joinpartitions(x) => TableOperations.JoinedPartitions
    x |> TableOperations.joinpartitions() => TableOperations.JoinedPartitions

Take an input `x` that implements `Tables.partitions` and "join" the partitions into
a single, "long" table. Each column is lazily appended using `SentinelArrays.ChainedVector`
so each partition's column is a single chain, and all partitions together are treated as a
single column. This can be helpful for "materializing" a partitioned input if single-column
operations are desired. No copy of the input data is made to avoid excessive memory allocations.

The returned object, `TableOperations.JoinedPartitions`, satisfies itself the `Tables.columns`
interface, so access to individual columns is supported via `x.col1`, `x[1]`, or
`Tables.getcolumn(x, :col1)`, in addition to the normal Tables.jl compatibility with sink
functions, like `df = DataFrame(TableOperations.joinpartitions(x))`.
"""
function joinpartitions(x)
    schema = Ref{Tables.Schema}()
    joined = ChainedVector[]
    N = 0
    for partition in Tables.partitions(x)
        cols = Tables.columns(partition)
        if isempty(joined)
            schema[] = Tables.schema(cols)
            N = length(schema[].names)
            foreach(i -> push!(joined, ChainedVector([Tables.getcolumn(cols, i)])), 1:N)
        else
            foreach(i -> append!(joined[i], Tables.getcolumn(cols, i)), 1:N)
        end
    end
    return JoinedPartitions(schema[], joined, Dict{Symbol, ChainedVector}(nm => col for (nm, col) in zip(schema[].names, joined)))
end

joinpartitions() = x -> joinpartitions(x)

struct NarrowTypes{T}
    x::T
    schema::Tables.Schema
end

narrow_arr(x) = mapreduce(typeof, promote_type, x)
narrow_types(t) = NarrowTypes(t, Tables.Schema(Tables.columnnames(t), [narrow_arr(getproperty(t, nm)) for nm in Tables.columnnames(t)]))

Tables.getcolumn(nt::NarrowTypes, nm::Symbol) = Vector{getproperty(nt.schema.types, nm)}(Tables.getcolumn(getfield(nt, 1), nm))
Tables.getcolumn(nt::NarrowTypes, i::Int) = Vector{nt.schema.types[i]}(Tables.getcolumn(getfield(nt, 1), i))

Tables.columnnames(nt::NarrowTypes) = Tables.columnnames(getfield(nt, 1)) # or nt.sch.names?
Tables.schema(nt::NarrowTypes) = nt.schema

Tables.istable(::Type{<:NarrowTypes}) = true

end # module


text:
function set_wally(nodes, bdcon, wally, cellcenter, cellxmax, cellymax, icell)
	# serch wall
	swith_wall = zeros(4)   # x-, x+, y-, y+を検索する
							# 1:壁
	for i in 1:4
		if Int(bdcon[i][1]) == 2 || Int(bdcon[i][1]) == 7
			swith_wall[i] = 1
		end
	end

	# cal number of ponits
	nop = 0
	if swith_wall[1] == 1
		nop = nop + (cellymax-1)
	end
	if swith_wall[2] == 1
		nop = nop + (cellymax-1)
	end
	if swith_wall[3] == 1
		nop = nop + (cellxmax-1)
	end
	if swith_wall[4] == 1
		nop = nop + (cellxmax-1)
	end

	wallpoint = zeros(nop, 2) # (格子数, 座標)
	ite = 1

	# x-
	if swith_wall[1] == 1
		for j in 1+icell:cellymax-icell
			wallpoint[ite,1] = nodes[icell+1,j,1]
			wallpoint[ite,2] = nodes[icell+1,j,2]
			ite = ite + 1
		end
	end
	# x+
	if swith_wall[2] == 1
		for j in 1+icell:cellymax-icell
			wallpoint[ite,1] = nodes[cellxmax-icell+1,j,1]
			wallpoint[ite,2] = nodes[cellxmax-icell+1,j,2]
			ite = ite + 1
		end
	end
	# y-
	if swith_wall[3] == 1
		for i in 2:cellxmax-1
			wallpoint[ite,1] = nodes[i,icell+1,1]
			wallpoint[ite,2] = nodes[i,icell+1,2]
			ite = ite + 1
		end
	end
	# y+
	if swith_wall[4] == 1
		for i in 2:cellxmax-1
			wallpoint[ite,1] = nodes[i,cellymax-icell+1,1]
			wallpoint[ite,2] = nodes[i,cellymax-icell+1,2]
			ite = ite + 1
		end
	end

	# 距離格納用
	distance_list = zeros(nop,2)
	distance = 0

	# 最大値を2つ抽出
	for j in 1+icell:cellymax-icell
		for i in 1+icell:cellxmax-icell
			for n in 1:nop
				distance = ((cellcenter[i,j,1]-wallpoint[n,1])^2 +
							(cellcenter[i,j,2]-wallpoint[n,2])^2
							)^0.5
				distance_list[n,1] = n
				distance_list[n,2] = distance
			end

			qdistance_list = copy(distance_list)
			qdistance_list = quicksort(qdistance_list,1,nop,2,2)

			# wall_point
			point1 = Int(qdistance_list[1,1])
			point2 = Int(qdistance_list[2,1])

			dis_bottom = ((wallpoint[point1, 1] - wallpoint[point2, 1])^2 + (wallpoint[point1, 2] - wallpoint[point2, 2])^2)^0.5

			# O型格子の関係で同じ点が二つあるため
			if dis_bottom == 0.0
				point1 = Int(qdistance_list[1,1])
				point2 = Int(qdistance_list[3,1])

				dis_bottom = ((wallpoint[point1, 1] - wallpoint[point2, 1])^2 + (wallpoint[point1, 2] - wallpoint[point2, 2])^2)^0.5
			end

			dis1 = distance_list[point1, 2]
			dis2 = distance_list[point2, 2]

			# ヘロンの公式から距離を算出
			s = (dis_bottom + dis1 + dis2)*0.5
			surface = (s*(s-dis_bottom)*(s-dis1)*(s-dis2))^0.5
			wally[i,j] = surface/dis_bottom * 2
		end
	end

	# 仮想セルへの代入
	for i in 1+icell:cellxmax-icell
		for j in 1:icell
			wally[i,j] = wally[i,icell+1]
			wally[i,cellymax-j+1] = wally[i,cellymax-icell]
		end
	end
	for j in 1+icell:cellymax-icell
		for i in 1:icell
			wally[i,j] = wally[icell+1,j]
			wally[cellxmax-i+1,j] = wally[cellxmax-icell,j]
		end
	end

	return wally, swith_wall
end

function cal_yplus(yplus, Qbase, wally, swith_wall, mu, cellxmax, cellymax, vecAx, vecAy, volume, icell)
	tvAxx = 0.0
    tvAxy = 0.0
    Axx = 0.0
    Axy = 0.0
    tvAyx = 0.0
    tvAyy = 0.0
    Ayx = 0.0
    Ayy = 0.0

	if swith_wall[1] == 1
		for j in 1+icell:cellymax-icell
			# 壁に並行な方向の速度を求めるためx-の壁面に対しvecAyを使用
			for i in 1:cellxmax
				tvAyx = 0.5*(vecAy[i,j-1,1]+vecAy[i,j,1])
				tvAyy = 0.5*(vecAy[i,j-1,2]+vecAy[i,j,2])
				Ayx = tvAyx / (tvAyx^2+tvAyy^2)^0.5
				Ayy = tvAyy / (tvAyx^2+tvAyy^2)^0.5
				nu = mu[i,j] / Qbase[i,j,1]
				ubar = abs(Ayx*Qbase[i,j,2] + Ayy*Qbase[i,j,3])
				yplus[i,j], up = cal_wall_val_spalding(ubar, wally[i,j], nu)
			end
		end
	end
	if swith_wall[2] == 1
		for j in 1+:cellymax-icell
			# 壁に並行な方向の速度を求めるためx+の壁面に対しvecAyを使用
			for i in 1:cellxmax
				tvAyx = 0.5*(vecAy[i,j-1,1]+vecAy[i,j,1])
				tvAyy = 0.5*(vecAy[i,j-1,2]+vecAy[i,j,2])
				Ayx = tvAyx / (tvAyx^2+tvAyy^2)^0.5
				Ayy = tvAyy / (tvAyx^2+tvAyy^2)^0.5
				nu = mu[i,j] / Qbase[i,j,1]
				ubar = abs(Ayx*Qbase[i,j,2] + Ayy*Qbase[i,j,3])
				yplus[i,j], up = cal_wall_val_spalding(ubar, wally[i,j], nu)
			end
		end
	end
	if swith_wall[3] == 1
		for j in 1:cellymax
			for i in 1+icell:cellxmax-icell
			# 壁に並行な方向の速度を求めるためy-の壁面に対しvecAxを使用

				tvAxx = 0.5*(vecAx[i-1,j,1]+vecAx[i,j,1])
				tvAxy = 0.5*(vecAx[i-1,j,2]+vecAx[i,j,2])
				Axx = tvAxx / (tvAxx^2+tvAxy^2)^0.5
				Axy = tvAxy / (tvAxx^2+tvAxy^2)^0.5
				nu = mu[i,j] / Qbase[i,j,1]
				ubar = abs(Axx*Qbase[i,j,2] + Axy*Qbase[i,j,3])
				yplus[i,j], up = cal_wall_val_spalding(ubar, wally[i,j], nu)

				if yplus[i,j] > 500
					break
				end
			end
		end
	end
	if swith_wall[4] == 1
		for j in 1:cellymax
			# 壁に並行な方向の速度を求めるためy-の壁面に対しvecAxを使用
			for i in 1+icell:cellxmax-icell
				tvAxx = 0.5*(vecAx[i-1,j,1]+vecAx[i,j,1])
				tvAxy = 0.5*(vecAx[i-1,j,2]+vecAx[i,j,2])
				Axx = tvAxx / (tvAxx^2+tvAxy^2)^0.5
				Axy = tvAxy / (tvAxx^2+tvAxy^2)^0.5
				nu = mu[i,j] / Qbase[i,j,1]
				ubar = abs(Axx*Qbase[i,j,2] + Axy*Qbase[i,j,3])
				yplus[i,j], up = cal_wall_val_spalding(ubar, wally[i,j], nu)
			end
		end
	end

	return yplus
end

function cal_wall_val_spalding(u,y,nu)
	# spalding 定数
	k = 0.4
	B = 5.5

	# ニュートン法
	up = 200.0
	for i in 1:100
		old = up

		up = up - spalding(up,u,y,nu,k,B)/spalding_dash(up,u,y,nu,k,B)

		if abs(up-old)/up < 10^(-4)
			break
		end
	end

	yplus = u*y/nu/up
	return yplus, up
end

function spalding(up,u,y,nu,k,B)
	t = k*up
	F = - 1/up* u*y/nu + up + exp(-k*B)*(exp(t)-1-t-t^2/2-t^3/6)
	return F
end

function spalding_dash(up,u,y,nu,k,B)
	t = k*up
	Fdash = 1/up^2* u*y/nu + 1 + exp(-k*B)*(k*exp(t)-k-k^2*up-k^3*up^2/2)
	return Fdash
end

# test
#=
uu = 6.64037328
yy = 0.06299976860
nunu = 1.406107133e-5
println(cal_wall_val_spalding(uu,yy,nunu))
throw(UndefVarError(:x))
=#

function wallf_Van_Driest(yplus)
	f = 1 - exp(-yplus/26)
	return f
end


text:
f(x) = cos(x)
g(x) = -sin(x)
h(x) = -cos(x)

x = 2.0

x = optimize_newton(x, f, g, h)

@test x ≈ π   rtol = 1e-4


text:
export meantoecc, benchmark_kepler

function benchmark_kepler(times::Int)
    r = [8.59072560e+02, -4.13720368e+03, 5.29556871e+03]
    v = [7.37289205e+00, 2.08223573e+00, 4.39999794e-01]
    mu = 3.986004418e5
    el = elements(r, v, mu)
    worst = -Inf
    best = Inf
    total = 0.0
    meantoecc(pi/2, el[2])
    for i = 1:times
        gc_enable(false)
        t = @elapsed meantoecc(pi/2, el[2])
        gc_enable(true)
        if t > worst
            worst = t
        end
        if t < best
            best = t
        end
        total += t
    end
    println("[$(total/times),$best,$worst]")
end

function period(a::Float64, mu::Float64)
    return 2pi*sqrt(a^3/mu)
end

function newton(x0::Float64, func::Function, derivative::Function, maxiter::Int=50, tol::Float64=sqrt(eps()))
    p0 = x0
    for i = 1:maxiter
        p = p0 - func(p0)/derivative(p0)
        if abs(p - p0) < tol
            return p
        end
        p0 = p
    end
    error("Not converged.")
end

function meantoecc(M::Float64, ecc::Float64)
    kepler(E::Float64) = E - ecc*sin(E) - M
    kepler_der(E::Float64) = 1 - ecc*cos(E)
    return newton(M, kepler, kepler_der)
end

max_issues_repo_name: OpenAstrodynamics/benchmarks

max_forks_repo_name: OpenAstrodynamics/benchmarks


text:
"""
    struct Configuration{<:AbstractOrbital}

Represents a configuration -- a set of orbitals and their associated occupation number.
Furthermore, each orbital can be in one of the following _states_: `:open`, `:closed` or
`:inactive`.
"""
struct Configuration{O<:AbstractOrbital}
    orbitals::Vector{O}
    occupancy::Vector{Int}
    states::Vector{Symbol}
    function Configuration(
        orbitals::Vector{O},
        occupancy::Vector{Int},
        states::Vector{Symbol}=[:open for o in orbitals]) where {O<:AbstractOrbital}
        length(orbitals) == length(occupancy) ||
            throw(ArgumentError("Need to specify occupation numbers for all orbitals"))
        length(states) ≤ length(orbitals) ||
            throw(ArgumentError("Cannot provide more states than orbitals"))

        length(orbitals) == length(unique(orbitals)) ||
            throw(ArgumentError("Not all orbitals are unique"))

        if length(states) < length(orbitals)
            append!(states, repeat([:open], length(orbitals) - length(states)))
        end
        sd = setdiff(states, [:open, :closed, :inactive])
        isempty(sd) || throw(ArgumentError("Unknown orbital states $(sd)"))

        for i in eachindex(orbitals)
            occ = occupancy[i]
            occ < 1 && throw(ArgumentError("Invalid occupancy $(occ)"))
            orb = orbitals[i]
            degen_orb = degeneracy(orb)
            if occ > degen_orb
                if O <: Orbital || O <: SpinOrbital
                    throw(ArgumentError("Higher occupancy than possible for $(orb) with degeneracy $(degen_orb)"))
                elseif O <: RelativisticOrbital
                    orb_conj = flip_j(orb)
                    degen_orb_conj = degeneracy(orb_conj)
                    if occ == degen_orb + degen_orb_conj && orb_conj ∉ orbitals
                        occupancy[i] = degen_orb
                        push!(orbitals, orb_conj)
                        push!(occupancy, degen_orb_conj)
                        push!(states, states[i])
                    elseif orb_conj ∈ orbitals
                        throw(ArgumentError("Higher occupancy than possible for $(orb) with degeneracy $(degen_orb)"))
                    else
                        throw(ArgumentError("Can only specify higher occupancy for $(orb) if completely filling the $(orb),$(orb_conj) subshell (for total degeneracy of $(degen_orb+degen_orb_conj))"))
                    end
                end
            end
        end

        for i in eachindex(orbitals)
            states[i] == :closed && occupancy[i] < degeneracy(orbitals[i]) &&
                throw(ArgumentError("Can only close filled orbitals"))
        end

        p = sortperm(orbitals)

        new{O}(orbitals[p], occupancy[p], states[p])
    end
end

function Configuration(orbs::Vector{Tuple{O,Int,Symbol}}) where {O<:AbstractOrbital}
    orbitals = Vector{O}()
    occupancy = Vector{Int}()
    states = Vector{Symbol}()
    for (orb,occ,state) in orbs
        push!(orbitals, orb)
        push!(occupancy, occ)
        push!(states, state)
    end
    Configuration(orbitals, occupancy, states)
end

Configuration(orbital::O, occupancy::Int, state::Symbol=:open) where {O<:AbstractOrbital} =
    Configuration([orbital], [occupancy], [state])

Configuration{O}() where {O<:AbstractOrbital} =
    Configuration(O[], Int[], Symbol[])

const RelativisticConfiguration{N} = Configuration{RelativisticOrbital{N}}

issimilar(a::Configuration{<:O}, b::Configuration{<:O}) where {O<:AbstractOrbital} =
    a.orbitals == b.orbitals && a.occupancy == b.occupancy

Base.:(==)(a::Configuration{<:O}, b::Configuration{<:O}) where {O<:AbstractOrbital} =
    issimilar(a, b) && a.states == b.states

"""
    fill(configuration)

Ensure all orbitals are at their maximum occupancy.
"""
function Base.fill(config::Configuration)
    map(config) do (orb,occ,state)
        orb,degeneracy(orb),state
    end |> Configuration
end

function Base.close(config::Configuration)
    map(config) do (orb,occ,state)
        orb,occ,:closed
    end |> Configuration
end

function write_orbitals(io::IO, config::Configuration)
    for (i,(orb,occ,state)) in enumerate(config)
        i > 1 && write(io, " ")
        write(io, "$(orb)")
        occ > 1 && write(io, to_superscript(occ))
        state == :closed && write(io, "ᶜ")
        state == :inactive && write(io, "ⁱ")
    end
end

function Base.show(io::IO, config::Configuration{O}) where O
    nc = length(config)
    if nc == 0
        write(io, "∅")
        return
    end
    noble_core_name = get_noble_core_name(config)
    core_config = core(config)
    ncc = length(core_config)
    if !isnothing(noble_core_name)
        write(io, "[$(noble_core_name)]ᶜ")
        ngc = length(get_noble_gas(O, noble_core_name))
        core_config = core(config)
        if ncc > ngc
            write(io, ' ')
            write_orbitals(io, core_config[ngc+1:end])
        end
        nc > ncc && write(io, ' ')
    elseif ncc > 0
        write_orbitals(io, core_config)
    end
    write_orbitals(io, peel(config))
end

"""
    get_noble_core_name(config::Configuration)

Returns the name of the noble gas with the most electrons whose configuration still forms
the first part of the closed part of `config`, or `nothing` if no such element is found.

```jldoctest
julia> AtomicLevels.get_noble_core_name(c"[He] 2s2")
He

julia> AtomicLevels.get_noble_core_name(c"1s2c 2s2c 2p6c 3s2c")
Ne

julia> AtomicLevels.get_noble_core_name(c"1s2") === nothing
true
```
"""
function get_noble_core_name(config::Configuration{O}) where O
    nc = length(config)
    nc == 0 && return nothing
    core_config = core(config)
    ncc = length(core_config)
    if length(core_config) > 0
        for gas in Iterators.reverse(noble_gases)
            gas_cfg = get_noble_gas(O, gas)
            ngc = length(gas_cfg)
            if ncc ≥ ngc && issimilar(core_config[1:length(gas_cfg)], gas_cfg)
                return gas
            end
        end
    end
    return nothing
end

function state_sym(state::AbstractString)
    if state == "c"
        :closed
    elseif state == "i"
        :inactive
    else
        :open
    end
end

function core_configuration(::Type{O}, element::AbstractString, state::AbstractString) where {O<:AbstractOrbital}
    element ∉ keys(noble_gases_configurations[O]) && throw(ArgumentError("Unknown noble gas $(element)"))
    state = state_sym(state == "" ? "c" : state) # By default, we'd like cores to be frozen
    core_config = noble_gases_configurations[O][element]
    Configuration(core_config.orbitals, core_config.occupancy,
                  [state for o in core_config.orbitals])
end

function parse_orbital(::Type{O}, orb_str) where {O<:AbstractOrbital}
    m = match(r"^(([0-9]+|.)([a-z]|\[[0-9]+\])[-]{0,1})([0-9]*)([*ci]{0,1})$", orb_str)
    orbital_from_string(O, m[1]) , (m[4] == "") ? 1 : parse(Int, m[4]), state_sym(m[5])
end

function Base.parse(::Type{Configuration{O}}, conf_str::AbstractString) where {O<:AbstractOrbital}
    isempty(conf_str) && return Configuration{O}()
    orbs = split(conf_str, r"[\. ]")
    core_m = match(r"\[([a-zA-Z]+)\]([*ci]{0,1})", first(orbs))
    if core_m != nothing
        core_config = core_configuration(O, core_m[1], core_m[2])
        if length(orbs) > 1
            peel_config = Configuration(parse_orbital.(Ref(O), orbs[2:end]))
            Configuration(vcat(core_config.orbitals, peel_config.orbitals),
                          vcat(core_config.occupancy, peel_config.occupancy),
                          vcat(core_config.states, peel_config.states))
        else
            core_config
        end
    else
        Configuration(parse_orbital.(Ref(O), orbs))
    end
end

"""
    @c_str -> Configuration{Orbital}

Construct a [`Configuration`](@ref), representing a non-relativistic configuration, out of
a string.

```jldoctest
julia> c"1s2 2s"
1s² 2s

julia> c"[Kr] 4d10 5s2 4f2"
[Kr]ᶜ 4d¹⁰ 4f² 5s²
```
"""
macro c_str(conf_str)
    parse(Configuration{Orbital}, conf_str)
end

"""
    @rc_str -> Configuration{RelativisticOrbital}

Construct a [`Configuration`](@ref) representing a relativistic configuration out of a
string.

```jldoctest
julia> rc"[Ne] 3s 3p- 3p"
[Ne]ᶜ 3s 3p⁻ 3p

julia> rc"[Ne] 3s 3p-2 3p4"
[Ne]ᶜ 3s 3p⁻² 3p⁴
```
"""
macro rc_str(conf_str)
    parse(Configuration{RelativisticOrbital}, conf_str)
end

Base.getindex(conf::Configuration{O}, i::Integer) where O =
    (conf.orbitals[i], conf.occupancy[i], conf.states[i])
Base.getindex(conf::Configuration{O}, i::Union{<:UnitRange{<:Integer},<:AbstractVector{<:Integer}}) where O =
    Configuration([conf[ii] for ii in i])

Base.iterate(conf::Configuration{O}, (el, i)=(length(conf)>0 ? conf[1] : nothing,1)) where O =
    i > length(conf) ? nothing : (el, (conf[i==length(conf) ? i : i+1],i+1))

Base.length(conf::Configuration) = length(conf.orbitals)
Base.lastindex(conf::Configuration) = length(conf)
Base.eltype(conf::Configuration{O}) where O = (O,Int,Symbol)

function Base.isless(a::Configuration{<:O}, b::Configuration{<:O}) where {O<:AbstractOrbital}
    l = min(length(a),length(b))
    # If they are equal up to orbital l, designate the shorter config
    # as the smaller one.
    a[1:l] == b[1:l] && return length(a) == l
    norm_occ = (orb,w) -> 2w ≥ degeneracy(orb) ? degeneracy(orb) - w : w
    for ((orba,occa,statea),(orbb,occb,stateb)) in zip(a[1:l],b[1:l])
        if orba < orbb
            return true
        elseif orba == orbb
            # This is slightly arbitrary, but we designate the orbital
            # with occupation closest to a filled shell as the smaller
            # one.
            norm_occ(orba,occa) < norm_occ(orbb,occb) && return true
        else
            return false
        end
    end
    false
end

"""
    num_electrons(c::Configuration) -> Int

Return the number of electrons in the configuration.

```jldoctest
julia> num_electrons(c"1s2")
2

julia> num_electrons(rc"[Kr] 5s2 5p-2 5p2")
42
```
"""
num_electrons(c::Configuration) = sum(c.occupancy)

"""
    num_electrons(c::Configuration, o::AbstractOrbital) -> Int

Returns the number of electrons on orbital `o` in configuration `c`. If `o` is not part of
the configuration, returns `0`.

```jldoctest
julia> num_electrons(c"1s 2s2", o"2s")
2

julia> num_electrons(rc"[Rn] Qf-5 Pf3", ro"Qf-")
5

julia> num_electrons(c"[Ne]", o"3s")
0
```
"""
function num_electrons(c::Configuration, o::AbstractOrbital)
    idx = findfirst(isequal(o), c.orbitals)
    isnothing(idx) && return 0
    c.occupancy[idx]
end

Base.in(orb::O, conf::Configuration{O}) where {O<:AbstractOrbital} =
    orb ∈ conf.orbitals

Base.filter(f::Function, conf::Configuration) =
    conf[filter(j -> f(conf[j]...), eachindex(conf.orbitals))]

core(conf::Configuration) = filter((orb,occ,state) -> state == :closed, conf)
peel(conf::Configuration) = filter((orb,occ,state) -> state != :closed, conf)
inactive(conf::Configuration) = filter((orb,occ,state) -> state == :inactive, conf)
active(conf::Configuration) = filter((orb,occ,state) -> state != :inactive, peel(conf))
bound(conf::Configuration) = filter((orb,occ,state) -> isbound(orb), conf)
continuum(conf::Configuration) = filter((orb,occ,state) -> !isbound(orb), peel(conf))

parity(conf::Configuration) = mapreduce(o -> parity(o[1])^o[2], *, conf)
Base.count(conf::Configuration) = mapreduce(o -> o[2], +, conf)

function Base.replace(conf::Configuration{O₁}, orbs::Pair{O₂,O₃}) where {O<:AbstractOrbital,O₁<:O,O₂<:O,O₃<:O}
    src,dest = orbs
    orbitals = promote_type(O₁,O₂,O₃)[]
    append!(orbitals, conf.orbitals)
    occupancy = copy(conf.occupancy)
    states = copy(conf.states)

    i = findfirst(isequal(src), orbitals)
    isnothing(i) && throw(ArgumentError("$(src) not present in $(conf)"))

    j = findfirst(isequal(dest), orbitals)
    if isnothing(j)
        push!(orbitals, dest)
        push!(occupancy, 1)
        push!(states, :open)
    else
        occupancy[j] == degeneracy(dest) &&
            throw(ArgumentError("$(dest) already maximally occupied in $(conf)"))
        occupancy[j] += 1
    end

    occupancy[i] -= 1
    if occupancy[i] == 0
        deleteat!(orbitals, i)
        deleteat!(occupancy, i)
        deleteat!(states, i)
    end

    Configuration(orbitals, occupancy, states)
end

"""
    -(configuration::Configuration, orbital::AbstractOrbital[, n=1])

Remove `n` electrons in the orbital `orbital` from the configuration
`configuration`. If the orbital had previously been `:closed` or
`:inactive`, it will now be `:open`.
"""
function Base.:(-)(configuration::Configuration{O₁}, orbital::O₂, n::Int=1) where {O<:AbstractOrbital,O₁<:O,O₂<:O}
    orbitals = promote_type(O₁,O₂)[]
    append!(orbitals, configuration.orbitals)
    occupancy = copy(configuration.occupancy)
    states = copy(configuration.states)

    i = findfirst(isequal(orbital), orbitals)
    isnothing(i) && throw(ArgumentError("$(orbital) not present in $(configuration)"))

    occupancy[i] ≥ n ||
        throw(ArgumentError("Trying to remove $(n) electrons from orbital $(orbital) with occupancy $(occupancy[i])"))

    occupancy[i] -= n
    states[i] = :open
    if occupancy[i] == 0
        deleteat!(orbitals, i)
        deleteat!(occupancy, i)
        deleteat!(states, i)
    end

    Configuration(orbitals, occupancy, states)
end

function Base.:(+)(a::Configuration{O₁}, b::Configuration{O₂}) where {O<:AbstractOrbital,O₁<:O,O₂<:O}
    orbitals = promote_type(O₁,O₂)[]
    append!(orbitals, a.orbitals)
    occupancy = copy(a.occupancy)
    states = copy(a.states)
 for (orb,occ,state) in b
        i = findfirst(isequal(orb), orbitals)
        if isnothing(i)
            push!(orbitals, orb)
            push!(occupancy, occ)
            push!(states, state)
        else
            occupancy[i] += occ
            states[i] == state || throw(ArgumentError("Incompatible states for $(orb): $(states[i]) and $state"))
        end
    end
    Configuration(orbitals, occupancy, states)
end

"""
    delete!(c::Configuration, o::AbstractOrbital)

Remove the entire subshell corresponding to orbital `o` from configuration `c`.

```jldoctest
julia> delete!(c"[Ar] 4s2 3d10 4p2", o"4s")
[Ar]ᶜ 3d¹⁰ 4p²
```
"""
function Base.delete!(c::Configuration{O}, o::O) where O <: AbstractOrbital
    idx = findfirst(isequal(o), [co[1] for co in c])
    idx === nothing && return c
    deleteat!(c.orbitals, idx)
    deleteat!(c.occupancy, idx)
    deleteat!(c.states, idx)
    return c
end

"""
    ⊗(::Union{Configuration, Vector{Configuration}}, ::Union{Configuration, Vector{Configuration}})

Given two collections of `Configuration`s, it creates an array of `Configuration`s with all
possible juxtapositions of configurations from each collection.

# Examples

```jldoctest
julia> [c"1s", c"2s"] ⊗ [c"2p-", c"2p"]
4-element Array{Configuration{RelativisticOrbital},1}:
 1s 2p⁻
 1s 2p
 2s 2p⁻
 2s 2p

julia> c"1s" ⊗ [c"2s2", c"2s 2p-"]
2-element Array{Configuration{RelativisticOrbital},1}:
 1s 2s²
 1s 2s 2p⁻
```
"""
⊗(a::Vector{<:Configuration}, b::Vector{<:Configuration}) =
    [x+y for x in a for y in b]
⊗(a::Union{<:Configuration,Vector{<:Configuration}}, b::Configuration) =
    a ⊗ [b]
⊗(a::Configuration, b::Vector{<:Configuration}) =
    [a] ⊗ b

"""
    rconfigurations_from_orbital(n, ℓ, occupancy)

Generate all `Configuration`s with relativistic orbitals corresponding to the
non-relativistic orbital with `n` and `ℓ` quantum numbers, with given occupancy.

# Examples

```jldoctest
julia> rconfigurations_from_orbital(3, 1, 2)
3-element Array{Configuration{RelativisticOrbital},1}:
 3p⁻²
 3p⁻ 3p
 3p²
```
"""
function rconfigurations_from_orbital(n::N, ℓ::Int, occupancy::Int) where {N<:MQ}
    n isa Integer && ℓ + 1 > n && throw(ArgumentError("ℓ=$ℓ too high for given n=$n"))
    occupancy > 2*(2ℓ + 1) && throw(ArgumentError("occupancy=$occupancy too high for given ℓ=$ℓ"))

    degeneracy_ℓm, degeneracy_ℓp = 2ℓ, 2ℓ + 2 # degeneracies of nℓ- and nℓ orbitals
    nlow_min = max(occupancy - degeneracy_ℓp, 0)
    nlow_max = min(degeneracy_ℓm, occupancy)
    confs = RelativisticConfiguration[]
    for nlow = nlow_max:-1:nlow_min
        nhigh = occupancy - nlow
        conf = if nlow == 0
            Configuration([RelativisticOrbital(n, ℓ, ℓ + 1//2)], [nhigh])
        elseif nhigh == 0
            Configuration([RelativisticOrbital(n, ℓ, ℓ - 1//2)], [nlow])
        else
            Configuration(
                [RelativisticOrbital(n, ℓ, ℓ - 1//2),
                 RelativisticOrbital(n, ℓ, ℓ + 1//2)],
                [nlow, nhigh]
            )
        end
        push!(confs, conf)
    end
    return confs
end

"""
    rconfigurations_from_orbital(orbital::Orbital, occupancy)

Generate all `Configuration`s with relativistic orbitals corresponding to the
non-relativistic version of the `orbital` with a given occupancy.

# Examples

```jldoctest
julia> rconfigurations_from_orbital(o"3p", 2)
3-element Array{Configuration{RelativisticOrbital},1}:
 3p⁻²
 3p⁻ 3p
 3p²
```
"""
function rconfigurations_from_orbital(orbital::Orbital, occupation::Integer)
    rconfigurations_from_orbital(orbital.n, orbital.ℓ, occupation)
end

function rconfigurations_from_nrstring(orb_str::AbstractString)
    m = match(r"^([0-9]+|.)([a-z]+)([0-9]+)?$", orb_str)
    isnothing(m) && throw(ArgumentError("Invalid orbital string: $(orb_str)"))
    n = parse_orbital_n(m)
    ℓ = parse_orbital_ℓ(m)
    occupancy = isnothing(m[3]) ? 1 : parse(Int, m[3])
    return rconfigurations_from_orbital(n, ℓ, occupancy)
end

"""
    @rcs_str -> Vector{Configuration{RelativisticOrbital}}

Construct a `Vector` of all `Configuration`s corresponding to the non-relativistic `nℓ`
orbital with the given occupancy from the input string.

The string is assumed to have the following syntax: `\$(n)\$(ℓ)\$(occupancy)`, where `n`
and `occupancy` are integers, and `ℓ` is in spectroscopic notation.

# Examples

```jldoctest
julia> rcs"3p2"
3-element Array{Configuration{Int64},1}:
 3p⁻²
 3p⁻ 3p
 3p²
```
"""
macro rcs_str(s)
    rconfigurations_from_nrstring(s)
end

"""
    spin_configurations(configuration)

Generate all possible configurations of spin-orbitals from
`configuration`, i.e. all permissible values for the quantum numbers
`n`, `ℓ`, `mℓ`, `ms` for each electron. Example:

    spin_configuration(c"1s2") -> 1s₀α 1s₀β
"""
function spin_configurations(c::Configuration{O}) where {O<:Orbital}
    states = Dict{Orbital,Symbol}()
    orbitals = map(c) do (orb,occ,state)
        states[orb] = state
        sorbs = spin_orbitals(orb)
        comb = collect(combinations(sorbs, occ))
        O == Orbital ? Vector{Vector{SpinOrbital}}(comb) : comb
    end
    map(allchoices(orbitals)) do choice
        c = vcat(choice...)
        s = [states[orb.orb] for orb in c]
        Configuration(c, ones(Int,length(c)), s)
    end
end

"""
    spin_configurations(configurations)

For each configuration in `configurations`, generate all possible
configurations of spin-orbitals.
"""
spin_configurations(cs::Vector{<:Configuration}) =
    sort(vcat(map(spin_configurations, cs)...))

function Base.show(io::IO, c::Configuration{<:SpinOrbital})
    orbitals = Dict{Orbital, Vector{<:SpinOrbital}}()
    core_orbitals = sort(unique(map(o -> o.orb, core(c).orbitals)))
    core_cfg = Configuration(core_orbitals, ones(Int, length(core_orbitals))) |> fill |> close

    if !isempty(core_cfg)
        show(io, core_cfg)
        write(io, " ")
    end
    for orb in peel(c).orbitals
        orbitals[orb.orb] = push!(get(orbitals, orb.orb, SpinOrbital[]), orb)
    end
    map(sort(collect(keys(orbitals)))) do orb
        ℓ = orb.ℓ
        g = degeneracy(orb)
        sub_shell = orbitals[orb]
        if length(sub_shell) == g
            format("{1:s}{2:s}", orb, to_superscript(g))
        else
            map(mℓrange(orb)) do mℓ
                mℓshell = findall(o -> o.mℓ == mℓ, sub_shell)
                if length(mℓshell) == 2
                    format("{1:s}{2:s}{3:s}", orb, to_subscript(mℓ), to_superscript(2))
                elseif length(mℓshell) == 1
                    string(sub_shell[mℓshell[1]])
                else
                    ""
                end
            end |> so -> join(filter(s -> !isempty(s), so), " ")
        end
    end |> so -> write(io, join(so, " "))
end

"""
    substitutions(src::Configuration{<:SpinOrbital}, dst::Configuration{<:SpinOrbital})

Find all orbital substitutions going from spin-configuration `src` to
configuration `dst`.
"""
function substitutions(src::Configuration{A}, dst::Configuration{B}) where {A<:SpinOrbital,B<:SpinOrbital}
    src == dst && return []
    # This is only valid for spin-configurations, since occupation
    # numbers are not dealt with.
    num_electrons(src) == num_electrons(dst) ||
        throw(ArgumentError("Configurations not of same amount of electrons"))
    r = Vector{Pair{A,B}}()
    missing = A[]
    same = Int[]
    for (i,(orb,occ,state)) in enumerate(src)
        j = findfirst(isequal(orb), dst.orbitals)
        if j === nothing
            push!(missing, orb)
        else
            push!(same, j)
        end
    end
    new = [j for j ∈ 1:num_electrons(dst)
           if j ∉ same]
    [mo => dst.orbitals[j] for (mo,j) in zip(missing, new)]
end

# We need to declare noble_gases first, with empty entries for Orbital and RelativisticOrbital
# since parse(Configuration, ...) uses it.
const noble_gases_configurations = Dict(
    O => Dict{String,Configuration{<:O}}()
    for O in [Orbital, RelativisticOrbital]
)
const noble_gases_configuration_strings = [
    "He" => "1s2",
    "Ne" => "[He] 2s2 2p6",
    "Ar" => "[Ne] 3s2 3p6",
    "Kr" => "[Ar] 3d10 4s2 4p6",
    "Xe" => "[Kr] 4d10 5s2 5p6",
    "Rn" => "[Xe] 4f14 5d10 6s2 6p6",
]
const noble_gases = [gas for (gas, _) in noble_gases_configuration_strings]
for (gas, configuration) in noble_gases_configuration_strings, O in [Orbital, RelativisticOrbital]
    noble_gases_configurations[O][gas] = parse(Configuration{O}, configuration)
end

# This construct is needed since when showing configurations, they
# will be specialized on the Orbital parameterization, which we cannot
# index noble_gases with.
for O in [Orbital, RelativisticOrbital]
    @eval get_noble_gas(::Type{<:$O}, k) = noble_gases_configurations[$O][k]
end

export Configuration, @c_str, @rc_str,
    num_electrons, core, peel, active, inactive, bound, continuum, parity, ⊗, @rcs_str,
    spin_configurations, substitutions