text:
# Based on:
# E. Cáceres, G.N. Gatica, F.A. Sequeira, A mixed virtual element method for
# the Brinkman problem, Math. Models Methods Appl. Sci. 27 (2017), no. 4,
# 707–743.

using jFEMTools
import Tensors
const jF = jFEMTools;
import FastGaussQuadrature
using LinearAlgebra

order = 2;
dim = 2;

mesh = unitSquareMesh2(RectangleCell, (3,3));

# Assamble local matrices
cell = jF.getcells(mesh)[1]
ci = 1

vertices = getverticescoords(mesh, ci)
tess = jF.get_2Dtesselation(vertices) #tesselate cell
centroid = cell_centroid(mesh, ci)
diameter = cell_diameter(mesh, ci)

# Init basis
Pe_ref_basis = jF.Monomials(dim-1,order,Tensors.Vec((0.5,)),1.0)
Pkm_basis = jF.Monomials(dim,order+1,centroid,diameter)
Pk_basis = jF.Monomials(dim,order,centroid,diameter)

points, weights = FastGaussQuadrature.gausslobatto(2*order-1)
# Shift interval from (-1,1) to (0,1)
weights *= 0.5
points = points .+ 1.0; points /= 2.0

# Assemble M_{mass,1} matrix
ϕdofs = jF.getnbasefunctions(Pe_ref_basis)

Mmass1ref = zeros(Float64,ϕdofs,ϕdofs)
for i in 1:ϕdofs, j = 1:ϕdofs, k = 1:size(points,1)
  Mmass1ref[i,j] += jF.value(Pe_ref_basis,i,Tensors.Vec((points[k],)))*jF.value(Pe_ref_basis,j,Tensors.Vec((points[k],)))*weights[k]
end
Mmass1ref_inv = inv(Mmass1ref)

# Assemble M_{edge} matrix
mdofs = jF.getnbasefunctions(Pkm_basis)
Medges = Vector{Matrix{Float64}}()
for e in 1:jF.getnedges(mesh,cell)
  Medge = zeros(Float64,ϕdofs,mdofs)
  edge_coords = jF.getverticescoords(mesh, jF.EdgeIndex(ci,e))
  he = norm(edge_coords[2] - edge_coords[1])
  for i in 1:ϕdofs, j = 1:mdofs, k = 1:size(points,1)
    mappedPoint = jF.map_unit_to_segment(Tensors.Vec((points[k],)),edge_coords)
    Medge[i,j] += jF.value(Pe_ref_basis,i,Tensors.Vec((points[k],)))*jF.value(Pkm_basis,j,mappedPoint)*weights[k]
  end
  push!(Medges,he*Medge)
end

# Assemble M_{mass,2}^{k+1}
Mmass2m = zeros(Float64,mdofs,mdofs)
for i in 1:mdofs, j = 1:mdofs
  αβ = Pkm_basis.indices[i] .+ Pkm_basis.indices[j]
  for e in 1:jF.getnedges(mesh,cell)
    edge_coords = jF.getverticescoords(mesh, jF.EdgeIndex(ci,e))
    normal = normal = jF.get_Normal(mesh, jF.EdgeIndex(ci, e))
    for k = 1:size(points,1)
      mappedPoint = jF.map_unit_to_segment(Tensors.Vec((points[k],)),edge_coords)
      Mmass2m[i,j] += diameter/(2+sum(αβ))*
                      dot(((1/diameter*(mappedPoint-centroid)).^αβ), normal)
    end
  end
end

#Assemble M_{mass,2} matrix
mkdofs = jF.getnbasefunctions(Pk_basis)
Mmass2 = view(Mmass2m,1:mkdofs,1:mkdofs)

#Assemble M_{grad}
Mgrad = zeros(Float64,mdofs-1,mdofs-1)
for i in 2:mdofs, j = 2:mdofs
  for e in 1:jF.getnedges(mesh,cell)
    edge_coords = jF.getverticescoords(mesh, jF.EdgeIndex(ci,e))
    normal = normal = jF.get_Normal(mesh, jF.EdgeIndex(ci, e))
    for k = 1:size(points,1)
      mappedPoint = jF.map_unit_to_segment(Tensors.Vec((points[k],)),edge_coords)
      # Product of δx derivatives
      α = Pkm_basis.indices[i] .- (1,0)
      β = Pkm_basis.indices[j] .- (1,0)
      if !(any(x -> x < 0, α) || any(x -> x < 0, β))
        Mgrad[i-1,j-1] += 1/(2+sum(α.+β))*((α[1]+β[1])/diameter)*
                        dot(((1/diameter*(mappedPoint-centroid)).^(α.+β)), normal)
      end
      # Product of δy derivatives
      α = Pkm_basis.indices[i] .- (0,1)
      β = Pkm_basis.indices[j] .- (0,1)
      if !(any(x -> x < 0, α) || any(x -> x < 0, β))
        Mgrad[i-1,j-1] += 1/(2+sum(α.+β))*((α[1]+β[1])/diameter)*
                        dot(((1/diameter*(mappedPoint-centroid)).^(α.+β)), normal)
      end
    end
  end
end

# Compute basis for G_k^{⟂}
gperpdofs = Int(order*(order+1)/2)
M0 = zeros(Float64, mdofs-1,2*mdofs)

for i in 2:mdofs, j = 1:mdofs
  for e in 1:jF.getnedges(mesh,cell)
    edge_coords = jF.getverticescoords(mesh, jF.EdgeIndex(ci,e))
    normal = normal = jF.get_Normal(mesh, jF.EdgeIndex(ci, e))
    for k = 1:size(points,1)
      mappedPoint = jF.map_unit_to_segment(Tensors.Vec((points[k],)),edge_coords)
      # δx m_{i+2} m_{j}
      α = Pkm_basis.indices[i] .- (1,0)
      β = Pkm_basis.indices[j]
      if !(any(x -> x < 0, α))
        M0[i-1,j] += (α[1]+1)/(2+sum(α.+β))*
                        dot(((1/diameter*(mappedPoint-centroid)).^(α.+β)), normal)
      end
      # δy m_{i+2} m_{j}
      α = Pkm_basis.indices[i] .- (0,1)
      β = Pkm_basis.indices[j]
      if !(any(x -> x < 0, α))
        M0[i-1,j+mdofs] += (α[2]+1)/(2+sum(α.+β))*
                        dot(((1/diameter*(mappedPoint-centroid)).^(α.+β)), normal)
      end
    end
  end
end

A0 = nullspace(M0)[:,1:gperpdofs]

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
workspace();
using BoxChopper;

box = [0 1; 0 1; 0 1]
r₀ = [0,0,0]
nout = [-1,-1,1]

nr₀ = nout⋅r₀
@code_warntype calc_vcbits(box, nout, nr₀)

ext: jl


text:
export ParametrizedMode
export scan_param!

struct ParametrizedMode{K,Kₑ,Kₘ,Kθ,Kθₛ,NT<:NamedTuple}
    θ::NT  # named tuple of parameters to scan; entries should be vectors of same length; must have θ.ω
    β::ArrComplexF{Kθ}  # β(θ)
    E::NTuple{Kₑ,ArrComplexF{Kθₛ}}  # E(θ)
    H::NTuple{Kₘ,ArrComplexF{Kθₛ}}  # H(θ)
end

function ParametrizedMode(θ::NamedTuple, mdl::Model{K,Kₑ,Kₘ}) where {K,Kₑ,Kₘ}
    haskey(θ, :ω) || @error "θ = $θ should have entry named ω."
    all(issorted.(values(θ))) || @error "Each entry of θ should be sorted in ascending order for interpolation."

    szθ = length.(values(θ))
    Kθ = length(szθ)

    szₛ = size(mdl.grid)  # size of shape dimension
    Kθₛ = Kθ + length(szₛ)

    β = ArrComplexF(undef, szθ)
    E = ntuple(k->ArrComplexF(undef, szθ..., szₛ...), Val(Kₑ))
    H = ntuple(k->ArrComplexF(undef, szθ..., szₛ...), Val(Kₘ))

    return ParametrizedMode{K,Kₑ,Kₘ,Kθ,Kθₛ,typeof(θ)}(θ,β,E,H)
end

Base.length(pm::ParametrizedMode) = length(pm.θ.ω)

function scan_param!(pm::ParametrizedMode{K,Kₑ,Kₘ,Kθ},
                     βguessₑ::Number,  # guess β for ending parameters (e.g., largest ω = shortest λ) in pm
                     mdl::Model{K,Kₑ,Kₘ};
                     update_model!::Any=(mdl,θ)->nothing,  # default: do not update model
                     nmode::Int=1  # mode number
                     ) where {K,Kₑ,Kₘ,Kθ}
    println("Begin parameter scan.")

    θ = pm.θ
    θkeys = keys(θ)
    θvals = values(θ)
    szθ = length.(θvals)
    Nθ = prod(szθ)

    CI = reverse(CartesianIndices(szθ))  # run from last corner of parameter space
    LI = reverse(LinearIndices(szθ))
    rng_shp = ntuple(k->Colon(), Val(K))  # represent all indices in shape dimensions
    t = @elapsed begin
        βguess = βguessₑ
        fₜref = ComplexF[]
        for χ = CI
            λ = LI[χ]  # linear index corresponding to χ
            println("\tScanning $λ out of $Nθ...")
            θval = t_ind(θvals, χ)
            θᵪ = NamedTuple{θkeys}(θval)
            ωᵪ = θᵪ.ω

            clear_objs!(mdl)
            update_model!(mdl, θᵪ)

            β, E, H, fₜref = calc_mode(mdl, ωᵪ, βguess; nmode, fₜref)

            pm.β[χ] = β  # filling from last corner of parameter space
            for k = 1:Kₑ
                view(pm.E[k], χ, rng_shp...) .= E[k]
            end
            for k = 1:Kₘ
                view(pm.H[k], χ, rng_shp...) .= H[k]
            end

            if λ ≠ Nθ  # still more parameters to examine, so update βguess and fₜref
                χnext = CI[λ+1]  # next Cartesian index
                ∆χ = χnext - χ
                sc = .!iszero.(∆χ.I)  # sc[k] = true if kth subscript will change
                if sum(sc) ≠ 1  # two or more subscripts will change
                    sc_max = findlast(sc)  # last dimension of changing subscript
                    ∆χclose = CartesianIndex(ntuple(k->(k==sc_max), Val(Kθ)))
                    χclose = χ + ∆χclose
                    βguess = pm.β[χclose]

                    # Below, use χclose.I... instead of χclose, because "iteration is
                    # deliberately unsupported for CartesianIndex."
                    Eguess = view.(pm.E, χclose.I..., rng_shp...)
                    Hguess = view.(pm.H, χclose.I..., rng_shp...)

                    fₜref = field2vec(Eguess, Hguess, mdl)
                else
                    βguess = β
                    # Don't update fₜref; use output of calc_mode().
                end
            end
        end
    end
    println("End parameter scan: $t seconds taken.")

    return nothing
end

include("dispersion.jl")


text:
module TestCore
using VoxelRayTracers
using Test

function collect_nonzero_voxels(itr)
    filter(collect(itr)) do hit
        hit.exit_time != hit.entry_time
    end
end
    #ray = (
    #    position = [0.01,-100, -100],
    #    velocity = [0.001, 1,1],
    #)
    #edgs = (-2:100.0, -50:50.0, sort!(randn(100)))
    #itr = @inferred eachtraversal(ray, edgs)
    #item, state = @inferred Nothing iterate(itr)
    #item, state = @inferred Nothing iterate(itr, state)
    #error()

@testset "1d" begin
    @testset "left to right" begin
        ray = (position=[0.0], velocity=[1.0])
        edges = (1:3,)
        itr = eachtraversal(ray, edges)
        hits = collect(itr)
        @test hits == [
            (voxelindex = CartesianIndex(1,), entry_time = 1.0, exit_time = 2.0),
            (voxelindex = CartesianIndex(2,), entry_time = 2.0, exit_time = 3.0),
           ]

        ray = (position=[0.0], velocity=[1.0])
        edges = ([1,2,3],)
        itr = eachtraversal(ray, edges)
        hits = collect(itr)
        @test hits == [
            (voxelindex = CartesianIndex(1,), entry_time = 1.0, exit_time = 2.0),
            (voxelindex = CartesianIndex(2,), entry_time = 2.0, exit_time = 3.0),
           ]
    end
    @testset "right to left" begin
        ray = (position=[100.0], velocity=[-1.0])
        edges = ([-10, -6.0, -2.0, 2.0], )
        itr = eachtraversal(ray, edges)
        hits = collect(itr)
        @test hits == [
            (voxelindex = CartesianIndex(3), entry_time = 98.0, exit_time = 102.0),
            (voxelindex = CartesianIndex(2), entry_time = 102.0, exit_time = 106.0),
            (voxelindex = CartesianIndex(1), entry_time = 106.0, exit_time = 110.0),
        ]
    end

    @testset "no hits" begin
        ray = (position=[3.0], velocity=[1.0])
        edges = (-3:2,)
        itr = @inferred eachtraversal(ray, edges)
        hits = @inferred collect(itr)
        @test isempty(hits)
    end

    @testset "start inside" begin
        ray = (position=[3.0], velocity=[1.0])
        edges = ([0, 2, 4, 6, 7],)
        itr = @inferred eachtraversal(ray, edges)
        hits = @inferred collect(itr)
        @test hits == [
            (voxelindex = CartesianIndex(2,), entry_time = 0.0, exit_time = 1.0),
            (voxelindex = CartesianIndex(3,), entry_time = 1.0, exit_time = 3.0),
            (voxelindex = CartesianIndex(4,), entry_time = 3.0, exit_time = 4.0),
        ]
    end

    @testset "start inside last voxel" begin
        ray = (position=[0], velocity=[1])
        edges = ([-2, 1],)
        hits = collect(eachtraversal(ray, edges))
        @test hits == [
            (voxelindex = CartesianIndex(1,), entry_time = 0.0, exit_time = 1.0),
        ]
    end
end
@testset "conversion" begin
    ray1 = (position=(1f0, 2f0), velocity=(1f0, 0f0))
    ray2 = (position=[1f0, 2f0], velocity=(1, 0))
    edges = (2f0:10f0, -3f0:3f0)
    @test eachtraversal(ray1, edges) == eachtraversal(ray2, edges)
    hit = first(eachtraversal(ray1, edges))
    @test hit.entry_time isa Float32
    @test hit.exit_time isa Float32
end

@testset "2d" begin

    @testset "default" begin
        pos = [0.0,-100]
        vel = [0.0,1]

        edges = (-2:5.0, [-10, -6.0, -2.0, 2.0])
        ray = (position=pos, velocity=vel)
        itr = eachtraversal(ray, edges)
        hits = @inferred collect(itr)
        @test eltype(hits) == typeof(first(hits))
        @test hits == [
            (voxelindex = CartesianIndex(3, 1), entry_time = 90.0, exit_time = 94.0),
            (voxelindex = CartesianIndex(3, 2), entry_time = 94.0, exit_time = 98.0),
            (voxelindex = CartesianIndex(3, 3), entry_time = 98.0, exit_time = 102.0),
        ]
    end

    @testset "diagonal" begin
        ray = (position=[0.0, 0.0], velocity=[1.0, 1.0])

        edges = ([1,3,4],1:4)
        itr = @inferred eachtraversal(ray, edges)
        @test collect_nonzero_voxels(itr) == [
         (voxelindex = CartesianIndex(1, 1), entry_time = 1.0, exit_time = 2.0),
         (voxelindex = CartesianIndex(1, 2), entry_time = 2.0, exit_time = 3.0),
         (voxelindex = CartesianIndex(2, 3), entry_time = 3.0, exit_time = 4.0),
        ]
        @test_broken collect(itr) == [
         (voxelindex = CartesianIndex(1, 1), entry_time = 1.0, exit_time = 2.0),
         (voxelindex = CartesianIndex(1, 2), entry_time = 2.0, exit_time = 3.0),
         (voxelindex = CartesianIndex(2, 3), entry_time = 3.0, exit_time = 4.0),
        ]
    end

    @testset "simultaneous entry on boundary wall and inner wall" begin
        ray = (position=[0.0, 5.0], velocity=[1.0, -2.0])
        edges = ([1, 10],[-10,3,4])
        itr = @inferred eachtraversal(ray, edges)
        ray = (position=[0.0, 5.0], velocity=[1.0, -2.0])
        edges = ([1, 10],[-10,3,4])
        itr = @inferred eachtraversal(ray, edges)
        @test collect_nonzero_voxels(itr) == [
            (voxelindex = CartesianIndex(1, 1), entry_time = 1.0, exit_time = 7.5)
        ]
        @test_broken collect(itr) == collect_nonzero_voxels(itr)
    end

    @testset "no hits, touch corner" begin
        ray = (position=[0,0], velocity=[1,1])
        edges = ([1,2], [-3,1])
        hits = eachtraversal(ray, edges)
        @test isempty(collect_nonzero_voxels(hits))
        @test_broken isempty(hits)
    end

    @testset "no hits, enter dim1 after leave dim2" begin
        ray = (position=[0,0], velocity=[1,2])
        edges = ([10, 11], [10,11])
        hits = eachtraversal(ray, edges)
        @test isempty(hits)
    end
end

@testset "3d" begin
    ray = (position=[10.0, 20, 30], velocity=[1.0, -1.0, -2])
    edges = ([11, 13, 14, 20], [5, 17, 18], [10, 25])
    itr = @inferred eachtraversal(ray, edges)
    @test collect_nonzero_voxels(itr) == [
        (voxelindex = CartesianIndex(1, 2, 1), entry_time = 2.5, exit_time = 3.0),
        (voxelindex = CartesianIndex(2, 1, 1), entry_time = 3.0, exit_time = 4.0),
        (voxelindex = CartesianIndex(3, 1, 1), entry_time = 4.0, exit_time = 10.0),
    ]
    @test_broken collect(itr) == collect_nonzero_voxels(itr)
end

@testset "inference" begin
    ray = (
        position = [0.01,-100, -100],
        velocity = [0.001, 1,1],
    )
    edgs = (-2:100.0, -50:50.0, sort!(randn(100)))
    itr = @inferred eachtraversal(ray, edgs)
    item, state = @inferred Nothing iterate(itr)
    item, state = @inferred Nothing iterate(itr, state)
end

@testset "eltype" begin
    for (P,V,E,T) in [
                      (Float32, Float32, Float32, Float32),
                      (Float64, Float64, Float64, Float64),
                      (Int64, Int64, Int64, Float64),
                      (Float64, Int64, Float32, Float64),
                   ]
        ray = (position=P[0,], velocity=V[1])
        edges = (E[1,2,3],)
        itr = @inferred eachtraversal(ray, edges)
        item, state = @inferred Nothing iterate(itr)
        item, state = @inferred Nothing iterate(itr, state)
        @test typeof(item.entry_time) == T
        @test typeof(item.exit_time) == T
    end
end

end#module


text:
using Test
using TDD

@testset "graph.jl" begin
        graph = [[1,2,3],
                [2,1],
                [3,1],
                [4,5],
                [5,4]]

        @test TDD.direct_neighbors(graph, 1) == [1,2,3]
        @test TDD.reachable(graph, 1) == [1,2,3]
        @test TDD.conn_comp(graph) == [[1,2,3],[4,5]]#2 connected sections

        #test if we get back self-connectivity even if the user doesn't specify it
        graph = [[2,3],
                [1],
                [1],
                [5],
                [4]]
        @test TDD.direct_neighbors(graph, 1) == [1,2,3]
        @test TDD.reachable(graph, 1) == [1,2,3]
        @test TDD.conn_comp(graph) == [[1,2,3],[4,5]]#2 connected sections

        graph = [[2,3],
                [1],
                [1],
                [5],
                [4],
                []]
        @test TDD.direct_neighbors(graph, 6) == [6]
        @test TDD.reachable(graph, 1) == [1,2,3]
        @test TDD.reachable(graph, 6) == [6]
        @test TDD.conn_comp(graph) == [[1,2,3],[4,5],[6]]#2 connected sections

        graph = [[1,2,3],
                [2,1,4],
                [3,1],
                [4,5],
                [5,4]]
        @test TDD.reachable(graph, 3) == [3,1,2,4,5]
        @test TDD.conn_comp(graph) == [[1,2,4,5,3]]#everything is connected

        graph = Bool[
                1 1 1 0 0
                1 1 0 0 0
                1 0 1 0 0
                0 0 0 1 1
                0 0 0 1 1]
        @test TDD.direct_neighbors(graph, 1) == [1,2,3]
        @test TDD.reachable(graph, 1) == [1,2,3]
        @test TDD.conn_comp(graph) == [[1,2,3],[4,5]]#2 connected sections

        graph = Bool[
                0 1 1 0 0
                1 0 0 0 0
                1 0 0 0 0
                0 0 0 0 1
                0 0 0 1 0]
        @test TDD.direct_neighbors(graph, 1) == [1,2,3]
        @test TDD.reachable(graph, 1) == [1,2,3]
        @test TDD.conn_comp(graph) == [[1,2,3],[4,5]]#2 connected sections

        graph = Bool[
                0 1 1 0 0 0
                1 0 0 0 0 0
                1 0 0 0 0 0
                0 0 0 0 1 0
                0 0 0 1 0 0
                0 0 0 0 0 0]
        @test TDD.direct_neighbors(graph, 6) == [6]
        @test TDD.reachable(graph, 1) == [1,2,3]
        @test TDD.reachable(graph, 6) == [6]
        @test TDD.conn_comp(graph) == [[1,2,3],[4,5],[6]]#2 connected sections
end


text:
# Density type that is useful across a variety of the sampling methods.
# NB: The density computed by f need not be normalized.

type Density
  f::Function
end

type DifferentiableDensity  # maybe <: DifferentiableFunction ?
  f::Function
  gradient::Function
end


text:
function smoothing_filter(N, filter_width)
    @assert N >= filter_width && isodd(filter_width)
    filter = zeros(Float32, N, N)
    half_width = Int((filter_width - 1) / 2)
    for i in 1:half_width
        filter[i, 1:half_width+i] .= 1 / (half_width+i)
        filter[end+1-i, end-(half_width+i-1):end] .= 1 / (half_width+i)
    end

    for i in half_width+1:N-half_width
        filter[i, i-half_width:i+half_width] .= 1 / filter_width
    end

    return filter
end


text:
#========================================================================================#
#	Laboratory 1
#
# Introduction to Julia, variables and data.
#
# Author: Niall Palfreyman, 22/12/2021
#========================================================================================#
[
	Activity(
		"""
		Hi! In this laboratory we use activities and questions that introduce us to the basics
		of the Julia language. Feel free to experiment at any time by entering your own
		commands at the Julia prompt.

		I will interpret anything you enter as an argument of the function reply() as an
		answer to the current question. Try this now - at the Julia prompt, enter the
		following code:

		city = "Ingolstadt"

		Then enter your answer as "reply(city)":
		""",
		"We're still getting acquainted - simply follow the above instructions precisely. :-)",
		x -> (x=="Ingolstadt")
	),
	Activity(
		"""
		You can use the Julia REPL like a calculator. Arithmetic operations have their usual
		priority, so tell me the result you get if you enter "(3+4)*5":
		""",
		"35",
		x -> x==35
	),
	Activity(
		"""
		A very useful feature of the Julia REPL is that you can always retrieve the result of
		your previous calculation by entering "ans" at the Julia prompt. Press up-arrow a few
		times until you see your previous calculation "(3+4)*5". Enter this calculation, and
		then enter "ans" at the Julia prompt. What value is returned?
		""",
		"35",
		x -> x==35
	),
	Activity(
		"""
		Remove the brackets: \"3+4*5\". First calculate this value, then retrieve it again
		by entering the variable \"ans\". Finally, show me the value using \"reply(ans)\":
		""",
		"23",
		x -> x==23
	),
	Activity(
		"What answer do you get if you remove the multiplication sign: \"5(3+4)\"?",
		"35",
		x -> x==35
	),
	Activity(
		"What do you get if you use floating-point division \"1 + 1/7\"?",
		"$(1+1/7)",
		x -> abs(x-(1+1/7)) < 0.01
	),
	Activity(
		"Or you can use exact division to produce a rational number \"a = 1 + 1//7\":",
		"$(1+1//7)",
		x -> x == 1+1//7
	),
	Activity(
		"You can then convert this to a floating-point number: \"Float64(a)\"",
		"$(1+1/7)",
		x -> abs(x-(1+1/7)) < 0.01
	),
	Activity(
		"""
		Let's try working with Vectors. A Vector is a linearly ordered list in square
		brackets. They can contain all kinds of elements. Enter this definition of the
		Vector b: \"b = [3,5,"4",[1,2]]\". Then use length() to find the length of b:
		""",
		"length(b)",
		x -> x == 4
	),
	Activity(
		"Use b[2] to access the 2nd element of b, and use the function isodd() to see if it is odd:",
		"isodd(b[2])",
		x -> x == true
	),
	Activity(
		"What is the type of (typeof) the 4th element of b?",
		"Vector{Int64}",
		x -> (x <: Vector)
	),
	Activity(
		"Change the 3rd element of b to my name \"Niall\" and find the length of that 3rd element",
		"b[3] = \"Niall\"",
		x -> x == 5
	),
	Activity(
		"""
		Ranges are closely related to Vectors. Investigate the expression \"2:4\", then use
		collect() to convert it into a Vector. Finally, use b[2:4] to extract the last three
		elements of b:
		""",
		"[5,\"4\",[1,2]",
		x -> x == [5,"Niall",[1,2]]
	),
	Activity(
		"""
		We can also work with BigInt numbers: \"c = factorial(big(100))\". Find out how
		many digits this number has by using digits() to convert it into a Vector.
		""",
		"158",
		x -> x == 158
	),
	Activity(
		"""
		Functions are centrally important in Julia, and the language offers us three
		different ways of defining them. Use this format to define your own squaring
		function, then use it to compute the square of -253:
		function sq1( x)
			x * x
		end
		""",
		"64009",
		x -> x == 64009
	),
	Activity(
		"""
		Now use this shorter format to define a new squaring function, and use it to
		compute the square of the imaginary number 5im:
		sq2( x) = x * x
		""",
		"-25",
		x -> x == -25
	),
	Activity(
		"""
		Finally, functions are so important in Julia that we often want to create an
		anonymous function quickly in the middle of a calculation. Use the following
		anonymous notation to compute the square of the complex number 1+im:
		(x -> x*x)(1+im)
		""",
		"2im",
		x -> x == 2im
	),
	Activity(
		"""
		Use tuples to assign values in parallel to the variables x and y:
		(x,y) = (2,3)
		Then give me (as a string) the code to swap the values of x and y:
		""",
		"(x,y) = (y,x)",
		x -> (r=replace(x," "=>""); r == "(x,y)=(y,x)" || r == "(y,x)=(x,y)")
	),
	Activity(
		"""
		You can use conditionals in the usual way, for example:

			function f(x)
				if x ≥ 0
					x
				else
					-x
				end
			end

		Now use the shorter notation "x ≥ 0 ? x : -x", and assign g equal to an anonymous function
		that calculates the factorial of its argument. That is, g(n) returns 1 if n ≤ 1, and
		returns n * g(n-1) otherwise. Use g to calculate 9!:
		""",
		"n -> (n ≤ 1 ? 1 : n * g(n-1))",
		x -> (x == 362880)
	),
	Activity(
		"""
		We use double quotation marks to denote Strings. Find out the typeof() "Niall":
		""",
		"typeof(\"Niall\")",
		x -> (x == String)
	),
	Activity(
		"""
		We can concatenate two or more Strings into a single String using the operator *.
		Define the two strings:

		name1 = \"Niall\"
		name2 = \"Palfreyman\"

		What result do you get if you concatenate these two Strings?
		""",
		"name1 * name2",
		x -> (x == "NiallPalfreyman")
	),
	Activity(
		"""
		I prefer my name with a space between the first and second names. We can do this
		using the call join([name1,name2],\" \"). What String results from this?
		""",
		"\"Niall Palfreyman\"",
		x -> (x == "Niall Palfreyman")
	),
	Activity(
		"""
		We can pattern text more simply using STRING INTERPOLATION. Enter this code:

		\"\$name1 \$name2\"

		The dollar sign (\$) interpolates the values of name1 and name2 into the string.
		Enter the following code and then see if you can work out the return value of
		the function println():

		x = π/2; println("The sine of \$x is \$(sin(x)).")

		(Note: You can enter the symbol π by typing "\\pi" and then pressing Tab.
		Alternatively, you can replace this symbol by the named constant "pi".)
		""",
		"\"typeof(ans)\"",
		x -> (x == Nothing)
	),
	Activity(
		"""
		Julia offers us many different functions for manipulating strings. Let's start
		by creating a String that we can play around with:

		a_string = "We are one and we are many, and from all the lands on Earth we come."

		Now ask whether the substring "Earth" occurs in the string a_quote: What is the
		return value of the function call occursin("Earth",a_string)?
		""",
		"",
		x -> x == true
	),
	Activity(
		"""
		What value does the function call startswith("Earth",a_string) return?
		""",
		"",
		x -> x == false
	),
	Activity(
		"""
		What value does the function call endswith("Earth",a_string) return?
		""",
		"",
		x -> x == false
	),
	Activity(
		"""
		Try out the following function calls: lowercase(a_string), uppercase(a_string)
		and titlecase(a_string).
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Try out the call: replace(a_string, "lands" => "countries"). Then find out the
		type of the second argument of this function call:
		""",
		"typeof(\"lands\" => \"countries\")",
		x -> (x <: Pair)
	),
	Activity(
		"""
		The call split(a_string, " ") divides a_string into tokens separated by the
		space symbol " ". How many tokens arise if we split a_string using the
		separator "an"?
		""",
		"split(a_string,\"an\")",
		x -> x == 5
	),
	Activity(
		"""
		We can convert numbers to strings. For example, what result do you get if
		you enter string(2.718)?
		""",
		"",
		x -> x == "2.718"
	),
	Activity(
		"""
		And converting back again, what result do you get if you enter
		parse(Float64,"2.718")?
		""",
		"parse() parses argument 2 according to the type of argument 1.",
		x -> x == 2.718
	),
]


text:
abstract AbstractPathState

# modified from http://stackoverflow.com/questions/25678112/insert-item-into-a-sorted-list-with-julia-with-and-without-duplicates
# returns true if insert succeeded, false if it was a duplicate
_insert_and_dedup!(v::Vector{Int}, x::Int) = isempty(splice!(v, searchsorted(v,x), x))

"""A type representing a single edge between two vertices of a graph."""
typealias Edge Pair{Int,Int}

"""Return source of an edge."""
src(e::Edge) = e.first
"""Return destination of an edge."""
dst(e::Edge) = e.second

 is_ordered(e::Edge) = src(e) <= dst(e)
@deprecate rev(e::Edge) reverse(e)

==(e1::Edge, e2::Edge) = (e1.first == e2.first && e1.second == e2.second)

"""A type representing an undirected graph."""
type Graph
    vertices::UnitRange{Int}
    ne::Int
    fadjlist::Vector{Vector{Int}} # [src]: (dst, dst, dst)
end

"""A type representing a directed graph."""
type DiGraph
    vertices::UnitRange{Int}
    ne::Int
    fadjlist::Vector{Vector{Int}} # [src]: (dst, dst, dst)
    badjlist::Vector{Vector{Int}} # [dst]: (src, src, src)
end

typealias SimpleGraph Union{Graph, DiGraph}

"""Return the vertices of a graph."""
vertices(g::SimpleGraph) = g.vertices

"""Return an iterator to the edges of a graph.
The returned iterator is valid for one pass over the edges, and is invalidated by changes to `g`.
"""
edges(g::SimpleGraph) = EdgeIter(g)

"""Returns the forward adjacency list of a graph.

The Array, where each vertex the Array of destinations for each of the edges eminating from that vertex.
This is equivalent to:

    fadj = [Vector{Int}() for _ in vertices(g)]
    for e in edges(g)
        push!(fadj[src(e)], dst(e))
    end
    fadj

For most graphs types this is pre-calculated.

The optional second argument take the `v`th vertex adjacency list, that is:

    fadj(g, v::Int) == fadj(g)[v]

NOTE: returns a reference, not a copy. Do not modify result.
"""
fadj(g::SimpleGraph) = g.fadjlist
fadj(g::SimpleGraph, v::Int) = g.fadjlist[v]

"""Returns true if all of the vertices and edges of `g` are contained in `h`."""
function issubset{T<:SimpleGraph}(g::T, h::T)
    (gmin, gmax) = extrema(vertices(g))
    (hmin, hmax) = extrema(vertices(h))
    return (hmin <= gmin <= gmax <= hmax) && issubset(edges(g), edges(h))
end

"""Add `n` new vertices to the graph `g`. Returns true if all vertices
were added successfully, false otherwise."""
function add_vertices!(g::SimpleGraph, n::Integer)
    added = true
    for i = 1:n
        added &= add_vertex!(g)
    end
    return added
end

"""Return true if the graph `g` has an edge from `u` to `v`."""
has_edge(g::SimpleGraph, u::Int, v::Int) = has_edge(g, Edge(u, v))

"""
    in_edges(g, v)

Returns an Array of the edges in `g` that arrive at vertex `v`.
`v=dst(e)` for each returned edge `e`.
"""
in_edges(g::SimpleGraph, v::Int) = [Edge(x,v) for x in badj(g, v)]

"""
    out_edges(g, v)

Returns an Array of the edges in `g` that depart from vertex `v`.
`v = src(e)` for each returned edge `e`.
"""
out_edges(g::SimpleGraph, v::Int) = [Edge(v,x) for x in fadj(g,v)]

"""Return true if `v` is a vertex of `g`."""
has_vertex(g::SimpleGraph, v::Int) = v in vertices(g)

"""
    nv(g)

The number of vertices in `g`.
"""
nv(g::SimpleGraph) = length(vertices(g))
"""
    ne(g)

The number of edges in `g`.
"""
ne(g::SimpleGraph) = g.ne

"""
    add_edge!(g, u, v)

Add a new edge to `g` from `u` to `v`.
Will return false if add fails (e.g., if vertices are not in the graph); true otherwise.
"""
add_edge!(g::SimpleGraph, u::Int, v::Int) = add_edge!(g, Edge(u, v))

"""
    rem_edge!(g, u, v)

Remove the edge from `u` to `v`.

Returns false if edge removal fails (e.g., if edge does not exist); true otherwise.
"""
rem_edge!(g::SimpleGraph, u::Int, v::Int) = rem_edge!(g, Edge(u, v))

"""
    rem_vertex!(g, v)

Remove the vertex `v` from graph `g`.
This operation has to be performed carefully if one keeps external data structures indexed by
edges or vertices in the graph, since internally the removal is performed swapping the vertices `v`  and `n=nv(g)`,
and removing the vertex `n` from the graph.
After removal the vertices in the ` g` will be indexed by 1:n-1.
This is an O(k^2) operation, where `k` is the max of the degrees of vertices `v` and `n`.
Returns false if removal fails (e.g., if vertex is not in the graph); true otherwise.
"""
function rem_vertex!(g::SimpleGraph, v::Int)
    v in vertices(g) || return false
    n = nv(g)

    edgs = in_edges(g, v)
    for e in edgs
        rem_edge!(g, e)
    end
    neigs = copy(in_neighbors(g, n))
    for i in neigs
        rem_edge!(g, Edge(i, n))
    end
    if v != n
        for i in neigs
            add_edge!(g, Edge(i, v))
        end
    end

    if is_directed(g)
        edgs = out_edges(g, v)
        for e in edgs
            rem_edge!(g, e)
        end
        neigs = copy(out_neighbors(g, n))
        for i in neigs
            rem_edge!(g, Edge(n, i))
        end
        if v != n
            for i in neigs
                add_edge!(g, Edge(v, i))
            end
        end
    end

    g.vertices = 1:n-1
    pop!(g.fadjlist)
    if is_directed(g)
        pop!(g.badjlist)
    end
    return true
end

"""Return the number of edges which start at vertex `v`."""
indegree(g::SimpleGraph, v::Int) = length(badj(g,v))
"""Return the number of edges which end at vertex `v`."""
outdegree(g::SimpleGraph, v::Int) = length(fadj(g,v))

indegree(g::SimpleGraph, v::AbstractArray{Int,1} = vertices(g)) = [indegree(g,x) for x in v]
outdegree(g::SimpleGraph, v::AbstractArray{Int,1} = vertices(g)) = [outdegree(g,x) for x in v]
degree(g::SimpleGraph, v::AbstractArray{Int,1} = vertices(g)) = [degree(g,x) for x in v]

"Return the maxium `outdegree` of vertices in `g`."
Δout(g) = noallocextreme(outdegree,(>), typemin(Int), g)
"Return the minimum `outdegree` of vertices in `g`."
δout(g) = noallocextreme(outdegree,(<), typemax(Int), g)
"Return the maximum `indegree` of vertices in `g`."
δin(g)  = noallocextreme(indegree,(<), typemax(Int), g)
"Return the minimum `indegree` of vertices in `g`."
Δin(g)  = noallocextreme(indegree,(>), typemin(Int), g)
"Return the minimum `degree` of vertices in `g`."
δ(g)    = noallocextreme(degree,(<), typemax(Int), g)
"Return the maximum `degree` of vertices in `g`."
Δ(g)    = noallocextreme(degree,(>), typemin(Int), g)

"computes the extreme value of `[f(g,i) for i=i:nv(g)]` without gathering them all"
function noallocextreme(f, comparison, initial, g)
    value = initial
    for i in 1:nv(g)
        funci = f(g, i)
        if comparison(funci, value)
            value = funci
        end
    end
    return value
end

"""
    degree_histogram(g)

Returns a `StatsBase.Histogram` of the degrees of vertices in `g`.
"""
degree_histogram(g::SimpleGraph) = fit(Histogram, degree(g))

"""Returns a list of all neighbors connected to vertex `v` by an incoming edge.

NOTE: returns a reference, not a copy. Do not modify result.
"""
in_neighbors(g::SimpleGraph, v::Int) = badj(g,v)
"""Returns a list of all neighbors connected to vertex `v` by an outgoing edge.

NOTE: returns a reference, not a copy. Do not modify result.
"""
out_neighbors(g::SimpleGraph, v::Int) = fadj(g,v)

"""Returns a list of all neighbors of vertex `v` in `g`.

For DiGraphs, this is equivalent to `out_neighbors(g, v)`.

NOTE: returns a reference, not a copy. Do not modify result.
"""
neighbors(g::SimpleGraph, v::Int) = out_neighbors(g, v)

"Returns the neighbors common to vertices `u` and `v` in `g`."
common_neighbors(g::SimpleGraph, u::Int, v::Int) = intersect(neighbors(g,u), neighbors(g,v))

"Returns true if `g` has any self loops."
has_self_loop(g::SimpleGraph) = any(v->has_edge(g, v, v), vertices(g))


text:
"""Restarts a calibration process."""

# Import modules to all processes
using ArgParse
using Glob
using CalibrateEDMF
using CalibrateEDMF.DistributionUtils
using CalibrateEDMF.Pipeline
using CalibrateEDMF.HelperFuncs
# Import EKP modules
using EnsembleKalmanProcesses
using EnsembleKalmanProcesses.ParameterDistributions
using JLD2

s = ArgParseSettings()
@add_arg_table s begin
    "--output_dir"
    help = "Results output directory"
    arg_type = String
    "--job_id"
    help = "Job identifier"
    arg_type = String
    default = "12345"
end
parsed_args = parse_args(ARGS, s)
# Recover inputs for restart
outdir_path = parsed_args["output_dir"]
include(joinpath(outdir_path, "config.jl"))
ekobjs = glob("ekobj_iter_*.jld2", outdir_path)
iters = @. parse(Int64, getfield(match(r"(?<=ekobj_iter_)(\d+)", basename(ekobjs)), :match))
last_iteration = maximum(iters)

priors = deserialize_prior(load(joinpath(outdir_path, "prior.jld2")))
ekobj = load(ekobj_path(outdir_path, last_iteration))["ekp"]
config = get_config()
restart_calibration(ekobj, priors, last_iteration, config, outdir_path, job_id = parsed_args["job_id"])


text:
# Autogenerated wrapper script for cilantro_jll for i686-w64-mingw32-cxx11
export libcilantro

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("cilantro")
JLLWrappers.@declare_library_product(libcilantro, "libcilantro.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcilantro,
        "bin\\libcilantro.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# MIT License

# Copyright (c) 2021 Dr. William A. Pisani

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
"""
    impropers(bonds::Array{Array{Int64,N} where N,1}) -> Array{NTuple{4,Int64},1}

Determines the impropers of a molecule or molecular system given the bonds.
"""
function impropers(bonds::Array{Array{Int64,N} where N,1})
    println("Determining impropers...")
    count_atomid1s = countmap(atomid[1] for atomid in bonds)
    count_atomid2s = countmap(atomid[2] for atomid in bonds)
    count_atomids = mergewith(+,count_atomid1s,count_atomid2s)

    # Get atomids with 3 or more bonds, those atoms have impropers
    atomids = [key for (key,value) in count_atomids if value >= 3]

    improper_array = Tuple{Int64,Int64,Int64,Int64}[]

    # Find atom ids connected to improper atoms
    for atomid in atomids
        atomid_connections = []
        for bond in bonds
            if atomid == bond[1]
                push!(atomid_connections,bond[2])
            elseif atomid == bond[2]
                push!(atomid_connections,bond[1])
            end
        end

        # Now that connections are found, get combinations of impropers
        if length(atomid_connections) == 3 # If only three bonds
            imp = (atomid_connections[1],atomid,atomid_connections[2],atomid_connections[3])
            push!(improper_array,imp)

        elseif length(atomid_connections) == 4 # If four bonds
            imp = (atomid_connections[1],atomid,atomid_connections[2],atomid_connections[3]) # 1-2-3
            push!(improper_array,imp)

            imp = (atomid_connections[2],atomid,atomid_connections[3],atomid_connections[4]) # 2-3-4
            push!(improper_array,imp)

            imp = (atomid_connections[1],atomid,atomid_connections[2],atomid_connections[4]) # 1-2-4
            push!(improper_array,imp)

            imp = (atomid_connections[1],atomid,atomid_connections[3],atomid_connections[4]) # 1-3-4
            push!(improper_array,imp)

        end
    end
    println("Impropers determined!")
    return improper_array
end

max_stars_repo_name: wapisani/CompMatSciTools

max_issues_repo_name: wapisani/CompMatSciTools

max_forks_repo_name: wapisani/CompMatSciTools


text:
"""
`module  JAC.BascisPerform`
    ... a submodel of JAC that contains methods that support tasks related to spectroscopic computation.
"""
module BascisPerform

    using Printf, JAC, ..AngularMomentum, ..Atomic, ..Basics, ..Bsplines, ..Cascade, ..Continuum, ..Defaults, ..Einstein, ..ManyElectron,
                  ..Nuclear, ..PlasmaShift, ..Radial

    export perform

    """
    `Basics.perform(computation::Atomic.Computation)`
        ... to perform the computation as prescribed by comp. All relevant intermediate and final results are printed to screen (stdout).
            Nothing is returned.

    `Basics.perform(computation::Atomic.Computation; output=true)`
        ... to perform the same but to return the complete output in a dictionary; the particular output depends on the type and
            specifications of the computations but can easily accessed by the keys of this dictionary.
    """
    function Basics.perform(computation::Atomic.Computation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel = computation.nuclearModel

        # Distinguish between the computation of level energies and properties and the simulation of atomic processes
        if   length(computation.configs) != 0
            ##x if  computation.asfSettings.breitCI
            ##x     @warn("perform: Test the Breit interaction; asfSettings = $(computation.asfSettings) ")
            ##x     wx = ManyElectron.Multiplet("from Ratip2012", "../test/testBreit/FeX-small-csl.inp",
            ##x                                     "../test/testBreit/FeX-small-scf.out","../test/testBreit/FeX-small-relci.mix-coulomb")
            ##x     basis     = wx.levels[1].basis
            ##x     multiplet = perform("computation: CI", basis, nModel, computation.grid, computation.asfSettings)
            ##x     error("stop after Breit test")
            ##x end =#
            basis     = perform("computation: SCF", computation.configs, nModel, computation.grid, computation.asfSettings)
            multiplet = perform("computation: CI", basis, nModel, computation.grid, computation.asfSettings)
            LSjj.expandLevelsIntoLS(multiplet, computation.asfSettings.jjLS)
            #
            if output    results = Base.merge( results, Dict("multiplet:" => multiplet) )
                        results = Base.merge( results, Dict("grid:"      => computation.grid) )  end

            # Now compute all requested properties
            if  JAC.EinsteinX      in computation.properties
                outcome = JAC.Einstein.computeLines(multiplet,        computation.grid, computation.einsteinSettings)
                if output    results = Base.merge( results, Dict("Einstein lines:" => outcome) )                  end
            end
            if  JAC.HFS            in computation.properties       && computation.hfsSettings.calcIJFexpansion
                outcome = JAC.Hfs.computeHyperfineMultiplet(multiplet, nModel, computation.grid, computation.hfsSettings)
                if output    results = Base.merge( results, Dict("IJF multiplet:" => outcome) )                   end
            end
            if  JAC.HFS            in computation.properties
                outcome = JAC.Hfs.computeOutcomes(multiplet, nModel,  computation.grid, computation.hfsSettings)
                if output    results = Base.merge( results, Dict("HFS outcomes:" => outcome) )                    end
            end
            if  JAC.LandeJ         in computation.properties
                outcome = JAC.LandeZeeman.computeOutcomes(multiplet, nModel,  computation.grid, computation.zeemanSettings)
                if output    results = Base.merge( results, Dict("Zeeman parameter outcomes:" => outcome) )       end
            end
            if  JAC.Isotope        in computation.properties
                outcome = JAC.IsotopeShift.computeOutcomes(multiplet, nModel, computation.grid, computation.isotopeSettings)
                if output    results = Base.merge( results, Dict("Isotope parameter outcomes:" => outcome) )      end
            end
            if  JAC.AlphaX         in computation.properties
                outcome = JAC.AlphaVariation.computeOutcomes(multiplet, nModel, computation.grid, computation.alphaSettings)
                if output    results = Base.merge( results, Dict("alpha variation parameter outcomes:" => outcome) )      end
            end
            if  JAC.FormF          in computation.properties
                outcome = JAC.FormFactor.computeOutcomes(multiplet, nModel, computation.grid, computation.formSettings)
                if output    results = Base.merge( results, Dict("Form factor outcomes:" => outcome) )            end
            end
            if  JAC.Yields         in computation.properties
                outcome = JAC.DecayYield.computeOutcomes(computation.configs, computation.asfSettings,
                                                        multiplet, nModel, computation.grid, computation.yieldSettings)
                if output    results = Base.merge( results, Dict("Fluorescence and AutoIonization yield outcomes:" => outcome) )   end
            end
            if  JAC.Green          in computation.properties
                outcome = JAC.GreenFunction.computeRepresentation(computation.configs, multiplet, nModel, computation.grid, computation.greenSettings)
                if output    results = Base.merge( results, Dict("Green function outcomes:" => outcome) )         end
            end
            if  JAC.Polarity       in computation.properties
                outcome = JAC.MultipolePolarizibility.computeOutcomes(multiplet, nModel, computation.grid, computation.polaritySettings)
                if output    results = Base.merge( results, Dict("Polarizibility outcomes:" => outcome) )         end
            end
            if  JAC.Plasma         in computation.properties
                outcome = JAC.PlasmaShift.computeOutcomes(multiplet, nModel, computation.grid, computation.asfSettings, computation.plasmaSettings)
                if output    results = Base.merge( results, Dict("Plasma shift outcomes:" => outcome) )            end
            end

        # Evaluate processes that need special SCF and CI procedures
        elseif  computation.process in [AugerInPlasma, PhotoInPlasma]
            pSettings        = computation.processSettings
            plasmaSettings   = PlasmaShift.Settings(pSettings.plasmaModel, pSettings.lambdaDebye, pSettings.ionSphereR0, pSettings.NoBoundElectrons)
            initialBasis     = perform("computation: SCF", computation.initialConfigs, nModel, computation.grid,
                                                        computation.initialAsfSettings)
            initialMultiplet = perform("computation: CI for plasma",  initialBasis, nModel, computation.grid, computation.initialAsfSettings,
                                                                    plasmaSettings)
            finalBasis       = perform("computation: SCF", computation.finalConfigs, nModel, computation.grid, computation.finalAsfSettings)
            finalMultiplet   = perform("computation: CI for plasma",  finalBasis, nModel, computation.grid, computation.finalAsfSettings,
                                                                    plasmaSettings)
            #
            if      computation.process == JAC.AugerInPlasma
                outcome = JAC.AutoIonization.computeLinesPlasma(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("AutoIonization lines in plasma:" => outcome) )         end
            elseif  computation.process == JAC.PhotoInPlasma
                outcome = JAC.PhotoIonization.computeLinesPlasma(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Photo lines in plasma:" => outcome) )                  end
            else
                error("stop a")
            end

        else
            initialBasis     = perform("computation: SCF", computation.initialConfigs, nModel, computation.grid, computation.initialAsfSettings)
            initialMultiplet = perform("computation: CI",  initialBasis, nModel, computation.grid, computation.initialAsfSettings)
            finalBasis       = perform("computation: SCF", computation.finalConfigs, nModel, computation.grid, computation.finalAsfSettings)
            finalMultiplet   = perform("computation: CI",  finalBasis, nModel, computation.grid, computation.finalAsfSettings)
            #
            if computation.process in [PhotoExcFluor, PhotoExcAuto, PhotoIonFluor, PhotoIonAuto, ImpactExcAuto, Dierec]
                intermediateBasis     = perform("computation: SCF", computation.intermediateConfigs, nModel, computation.grid,
                                                                    computation.intermediateAsfSettings)
                intermediateMultiplet = perform("computation: CI",  intermediateBasis, nModel, computation.grid,
                                                                    computation.intermediateAsfSettings)
            end
            #
            if      computation.process == JAC.Auger
                ##x println(" ")
                ##x for  i = 1:3   println("  perform-WARNING: Code modified to read in multiplets from Grasp computations and for testing integrals " *
                ##x                        "and Auger amplitudes !!")
                ##x end
                ##x initialMultiplet = JAC.ManyElectron.Multiplet("from Grasp2013", "TestAuger/belike-resonance-a-csl.inp",
                ##x                                               "TestAuger/belike-resonance-a-scf.out","TestAuger/belike-resonance-a-relci.mix")
                ##x finalMultiplet   = JAC.ManyElectron.Multiplet("from Grasp2013", "TestAuger/belike-ground-a-csl.inp",
                ##x                                               "TestAuger/belike-ground-a-scf.out","TestAuger/belike-ground-a-relci.mix")
                outcome = JAC.AutoIonization.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("AutoIonization lines:" => outcome) )                  end
            elseif  computation.process == JAC.Radiative
                ##x println(" ")
                ##x for  i = 1:3   println("  perform-WARNING: Code modified to read in multiplets from Grasp computations and for testing integrals " *
                ##x                        "and radiative amplitudes !!")
                ##x end
                ##x finalMultiplet   = JAC.ManyElectron.Multiplet("from Grasp2013", "TestRaditive/helike-resonance-a-csl.inp",
                ##x                                               "TestRaditive/helike-resonance-a-scf.out","TestRaditive/helike-resonance-a-relci.mix")
                ##x initialMultiplet = JAC.ManyElectron.Multiplet("from Grasp2013", "TestRaditive/helike-resonance-a-csl.inp",
                ##x                                               "TestRaditive/helike-resonance-a-scf.out","TestRaditive/helike-resonance-a-relci.mix")
                outcome = JAC.PhotoEmission.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("radiative lines:" => outcome) )              end
            elseif  computation.process == JAC.Photo
                outcome = JAC.PhotoIonization.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photoionization lines:" => outcome) )        end
            elseif  computation.process == JAC.Rec
                outcome = JAC.PhotoRecombination.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo recombination lines:" => outcome) )    end
            elseif  computation.process == JAC.Eimex
                outcome = JAC.ImpactExcitation.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("impact-excitation lines:" => outcome) )      end
            elseif  computation.process == JAC.PhotoExc
                outcome = JAC.PhotoExcitation.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-excitation lines:" => outcome) )       end
            elseif  computation.process == JAC.PairA1P
                outcome = JAC.PairAnnihilation1Photon.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("pair-annihilation-1-photon lines:" => outcome) )       end
            elseif  computation.process == JAC.MultiPhotonDE
                outcome = JAC.MultiPhotonDeExcitation.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("multi-photon-de-excitation lines:" => outcome) )       end
            elseif  computation.process == JAC.RAuger
                outcome = JAC.RadiativeAuger.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("radiative Auger sharings:" => outcome) )               end
            elseif  computation.process == JAC.Coulex
                outcome = JAC.CoulombExcitation.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Coulomb excitation lines:" => outcome) )               end
            elseif  computation.process == JAC.Coulion
                outcome = JAC.CoulombIonization.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Coulomb ionization lines:" => outcome) )               end
            elseif  computation.process == JAC.PhotoExcAuto
                outcome = JAC.PhotoExcitationAutoion.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, nModel,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-excitation-autoionization pathways:" => outcome) )      end
            elseif  computation.process == JAC.PhotoExcFluor
                outcome = JAC.PhotoExcitationFluores.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-excitation-fluorescence pathways:" => outcome) )        end
            elseif  computation.process == JAC.PhotoIonAuto
                outcome = JAC.PhotoIonizationAutoion.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-ionization-autoionization pathways:" => outcome) )      end
            elseif  computation.process == JAC.PhotoIonFluor
                outcome = JAC.PhotoIonizationFluores.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-ionizatiton-fluorescence pathways:" => outcome) )        end
            elseif  computation.process == JAC.ImpactExcAuto
                outcome = JAC.ImpactExcitationAutoion.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("impact-excitation-autoionization pathways:" => outcome) )     end
            elseif  computation.process == JAC.Compton
                outcome = JAC.RayleighCompton.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Rayleigh-Compton lines:" => outcome) )        end
            elseif  computation.process == JAC.MultiPI
                outcome = JAC.MultiPhotonIonization.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("multi-photon single ionization:" => outcome) )        end
            elseif  computation.process == JAC.MultiPDI
                outcome = JAC.MultiPhotonDoubleIon.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("multi-photon double ionization:" => outcome) )        end
            elseif  computation.process == JAC.InternalConv
                outcome = JAC.InternalConversion.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("internal conversion lines:" => outcome) )        end
            elseif  computation.process == JAC.Dierec
                outcome = JAC.Dielectronic.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, nModel,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("dielectronic recombination pathways:" => outcome) )           end
            else
                error("stop b")
            end
        end

        Defaults.warn(PrintWarnings)
        Defaults.warn(ResetWarnings)
        return( results )
    end

    """
    `Basics.perform(comp::Cascade.Computation)`
        ... to set-up and perform a cascade computation that starts from a given set of initial configurations and proceeds via
            various steps until a given number of electrons has been removed or the decay stops at some stable levels with regard
            to the given atomic processes. The results of all individual steps are printed to screen but nothing is returned
            otherwise.

    `Basics.perform(comp::Cascade.Computation; output=true)`
        ... to perform the same but to return the complete output in a dictionary;  the particular output depends on the type
            and specifications of the cascade but can easily accessed by the keys of this dictionary.
    """
    function Basics.perform(comp::Cascade.Computation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        # Perform the SCF and CI computation for the intial-state multiplet and print them out with their relative occupation
        basis     = perform("computation: SCF", comp.initialConfs, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
        multiplet = perform("computation: CI", basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
        multiplet = Multiplet("initial states", multiplet.levels)
        JAC.Cascade.displayInitialLevels(stdout, multiplet, comp.initialLevels)
        if  printSummary   JAC.Cascade.displayInitialLevels(iostream, multiplet, comp.initialLevels)    end
        if output    results = Base.merge( results, Dict("initial multiplet:" => multiplet) )           end
        #
        # Generate subsequent cascade configurations as well as display and group them together
        wa = JAC.Cascade.generateConfigurationList(comp.initialConfs, comp.maxElectronLoss, comp.NoShakeDisplacements)
        wb = JAC.Cascade.groupDisplayConfigurationList(comp.nuclearModel.Z, wa)
        #
        # Determine first all configuration 'blocks' and from them the individual steps of the cascade
        wc = JAC.Cascade.generateBlocks(comp, wb, basis.orbitals)
        JAC.Cascade.displayBlocks(stdout, wc)
        if  printSummary   JAC.Cascade.displayBlocks(iostream, wc)    end
        # Determine, modify and compute the transition data for all steps, ie. the PhotoEmission.Line's, the AutoIonization.Line's, etc.
        wd = JAC.Cascade.determineSteps(comp, wc)
        JAC.Cascade.displaySteps(stdout, wd)
        if  printSummary   JAC.Cascade.displaySteps(iostream, wd)    end
        we   = JAC.Cascade.modifySteps(wd)
        data = JAC.Cascade.computeSteps(comp, we)
        if output
            results = Base.merge( results, Dict("cascade data:" => data) )
            #
            #  Write out the result to file to later continue with simulations on the cascade data
            filename = "zzz-Cascade-" * string(now())[1:13] * ".jld"
            println("\nWrite all results to disk; use:\n   save($filename, results) \n   using JLD " *
                    "\n   results = load($filename)    ... to load the results back from file.")
            if  printSummary   println(iostream, "\nWrite all results to disk; use:\n   save($filename, results) \n   using JLD " *
                                                "\n   results = load($filename)    ... to load the results back from file." )      end
            save(filename, results)
        end
        ## return( results )
        return( data )
    end

    """
    `Basics.perform(simulation::Cascade.Simulation, data::Cascade.Data)`
        ... to simulate a cascade decay (and excitation) by using the given data::Cascade.Data. Different computational methods and
            different properties of the ionic system, such as the ion distribution or final-level distribution can be derived and
            displayed from these simulations. Of course, the details of these simulations strongly depend on the atomic processes
            and data that have been generated before by performing a computation::Cascade.Computation. The results of all individual
            steps are printed to screen but nothing is returned otherwise.

    `Basics.perform(simulation::Cascade.Simulation, data::Cascade.Data; output=true)`
        ... to perform the same but to return the complete output in a dictionary; the particular output depends on the method and
            specifications of the cascade but can easily accessed by the keys of this dictionary.
    """
    function Basics.perform(simulation::Cascade.Simulation, data::Cascade.Data; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        #
        # Distinguish between the different computational methods for running the simulations
        if  simulation.method == Cascade.ProbPropagation()
            if   JAC.Cascade.FinalDist()  in  simulation.properties   ||   JAC.Cascade.IonDist()  in  simulation.properties
                wa = Cascade.simulateLevelDistribution(simulation, data)
            end
            if   JAC.Cascade.ElectronIntensity()    in simulation.properties    ||   JAC.Cascade.PhotonIntensity()  in simulation.properties   ||
                JAC.Cascade.ElectronCoincidence()  in simulation.properties
                error("stop a: Not yet implemented")
            end
        else  error("stop b")
        end

        if output
            wx = 0.
            results = Base.merge( results, Dict("ion distribution:" => wx) )
            #  Write out the result to file to later continue with simulations on the cascade data
            filename = "xxx-simulation-" * string(now())[1:13] * ".jld"
            println("\nWrite all results to disk; use:\n   save($filename, results) \n   using JLD " *
                    "\n   results = load($filename)    ... to load the results back from file.")
            save(filename, results)
        end
        return( results )
    end

    """
    `Basics.perform(computation::Atomic.CasComputation)`
        ... to perform the computation ... . Nothing is returned.  **Not yet implemented !**

    `Basics.perform(computation::Atomic.CasComputation; output=true)`
        ... to perform the same ....  **Not yet implemented !**
    """
    function Basics.perform(computation::Atomic.CasComputation; output::Bool=false)
        error("Not yet implemented")
    end

    """
    `Basics.perform("computation: SCF", configs::Array{Configuration,1}, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings;
                                        printout::Bool=true)`
        ... to generate an atomic basis and to compute the self-consistent field (SCF) for this basis due to the given settings;
            a basis::Basis is returned.
    """
    function Basics.perform(sa::String, configs::Array{Configuration,1}, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings;
                    printout::Bool=true)
        !(sa == "computation: SCF")   &&   error("Unsupported keystring = $sa")
        if  printout    println("\n... in perform('computation: SCF', ...")    end

        # Generate a list of relativistic configurations and determine an ordered list of subshells for these configurations
        relconfList = ConfigurationR[]
        for  conf in configs
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        if  printout    for  i = 1:length(relconfList)    println("perform-aa: ", relconfList[i])    end   end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=printout)

        # Generate the relativistic CSF's for the given subshell list
        csfList = CsfR[]
        for  relconf in relconfList
            newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
            append!( csfList, newCsfs)
        end

        # Determine the number of electrons and the list of coreSubshells
        NoElectrons      = sum( csfList[1].occupation )
        coreSubshellList = Subshell[]
        for  k in 1:length(subshellList)
            mocc = Basics.subshell_2j(subshellList[k]) + 1;    is_filled = true
            for  csf in csfList
                if  csf.occupation[k] != mocc    is_filled = false;    break   end
            end
            if   is_filled    push!( coreSubshellList, subshellList[k])    end
        end

        ##x waL = JAC.Bsplines.generatePrimitives(7, 7, grid);    nsL = length(waL.bsplines) - 3
        ##x waS = JAC.Bsplines.generatePrimitives(8, 7, grid);    nsS = length(waS.bsplines) - 3
        wa = Bsplines.generatePrimitives(grid)

        # Generate start orbitals
        if  settings.startScf == "hydrogenic"
            # Generate start orbitals for the SCF field by using B-splines
            ##x orbitals  = JAC.Bsplines.generateOrbitalsHydrogenic(waL, nsL, waS, nsS, nuclearModel, subshellList)
            orbitals  = JAC.Bsplines.generateOrbitalsHydrogenic(wa, nuclearModel, subshellList; printout=printout)
        elseif  settings.startScf == "fromNRorbitals"
            # Generate starting orbitals for this csfList by adapting non-relativistic orbitals with a proper nuclear charge
            orbitals = Dict{Subshell, Orbital}()
            for  subsh in subshellList
                orb      = JAC.HydrogenicIon.radialOrbital(subsh, nuclearModel.Z, grid)
                orbitals = Base.merge( orbitals, Dict( subsh => orb) )
            end
        else  error("stop a")
        end

        ##x error("stop ** for test reasons ** ")

        # Perform the SCF computations for the orbitals due to the given settings
        basis    = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, orbitals)
        if   settings.methodScf in ["meanDFS", "meanHS"]
            ##x newBasis = JAC.Bsplines.solveSelfConsistentFieldMean(waL, nsL, waS, nsS, nuclearModel, basis, settings)
            newBasis = JAC.Bsplines.solveSelfConsistentFieldMean(wa, nuclearModel, basis, settings; printout=printout)
        elseif   settings.methodScf in ["AL", "OL"]
            ##x newBasis = JAC.Bsplines.solveSelfConsistentFieldOptimized(waL, nsL, waS, nsS, nuclearModel, basis, settings)
            newBasis = JAC.Bsplines.solveSelfConsistentFieldOptimized(wa, nuclearModel, basis, settings; printout=printout)
        else  error("stop b")
        end

        return( newBasis )
    end

    """
    `Basics.perform("computation: mutiplet from orbitals, no CI, CSF diagonal", configs::Array{Configuration,1},
                    initalOrbitals::Dict{Subshell, Orbital}, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings; printout::Bool=true)`
        ... to generate from the given initial orbitals a multiplet of single-CSF levels by just using the diagonal
            part of the Hamiltonian matrix; a multiplet::Multiplet is returned.
    """
    function Basics.perform(sa::String, configs::Array{Configuration,1}, initalOrbitals::Dict{Subshell, Orbital}, nuclearModel::Nuclear.Model,
                    grid::Radial.Grid, settings::AsfSettings; printout::Bool=true)
        !(sa == "computation: mutiplet from orbitals, no CI, CSF diagonal")   &&   error("Unsupported keystring = $sa")
        if  printout    println("\n... in perform('computation: mutiplet from orbitals, no CI, CSF diagonal', ...")    end

        # Generate a list of relativistic configurations and determine an ordered list of subshells for these configurations
        relconfList = ConfigurationR[]
        for  conf in configs
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        if  printout    for  i = 1:length(relconfList)    println("perform-aa: ", relconfList[i])    end   end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=printout)

        # Generate the relativistic CSF's for the given subshell list
        csfList = CsfR[]
        for  relconf in relconfList
            newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
            append!( csfList, newCsfs)
        end

        # Determine the number of electrons and the list of coreSubshells
        NoElectrons      = sum( csfList[1].occupation )
        coreSubshellList = Subshell[]
        for  k in 1:length(subshellList)
            mocc = Basics.subshell_2j(subshellList[k]) + 1;    is_filled = true
            for  csf in csfList
                if  csf.occupation[k] != mocc    is_filled = false;    break   end
            end
            if   is_filled    push!( coreSubshellList, subshellList[k])    end
        end

        # Set-up a basis for calculating the Hamiltonian matrix
        basis = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, initalOrbitals)

        # Calculate the diagonal part of the Hamiltonian matrix and define a multiplet for these approximate ASf;
        # Generate first an effective nuclear charge Z(r) on the given grid
        potential = JAC.Nuclear.nuclearPotential(nuclearModel, grid)

        n = length(csfList);    matrix = zeros(Float64, n, n)
        for  r = 1:n
            wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[r], basis.csfs[r])
            me = 0.
            for  coeff in wa[1]
                jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                me = me + coeff.T * sqrt( jj + 1) * JAC.RadialIntegrals.GrantIab(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
            end

            for  coeff in wa[2]
                if  settings.coulombCI
                    me = me + coeff.V * JAC.InteractionStrength.XL_Coulomb(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)   end
                if  settings.breitCI
                    me = me + coeff.V * JAC.InteractionStrength.XL_Breit(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid)     end
            end
            matrix[r,r] = me
        end
        #
        eigen  = Basics.diagonalize("matrix: Julia, eigfact", matrix)
        levels = Level[]
        for  ev = 1:length(eigen.values)
            vector = eigen.vectors[ev];   ns = 0
            for  r = 1:length(vector)
                if      (vector[r])^2 > 0.99    ns = r;   break
                elseif  (vector[r])^2 > 0.01    error("stop a; unexpected mixing coefficieint; vector = $vector ")
                end
            end
            J = csfList[ns].J
            newlevel = Level( J, AngularM64(J.num//J.den), csfList[ns].parity, 0, eigen.values[ev], 0., true, basis, vector )
            push!( levels, newlevel)
        end
        mp = Multiplet("noName", levels)
        mp = Basics.sortByEnergy(mp)

        # Display all level energies and energy splittings
        if  printout
            Basics.tabulate("multiplet: energies", mp)
            Basics.tabulate("multiplet: energy relative to immediately lower level",    mp)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", mp)
        end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary
            Basics.tabulate("multiplet: energies", mp, stream=iostream)
            Basics.tabulate("multiplet: energy relative to immediately lower level",    mp, stream=iostream)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", mp, stream=iostream)
        end

        return( mp )
    end

    """
    `Basics.perform("computation: CI", basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings; printout::Bool=true)`
        ... to  set-up and diagonalize from the (SCF) basis the configuration-interaction matrix and to derive and display the
            level structure of the corresponding multiplet due to the given settings; a multiplet::Multiplet is returned.
    """
    function Basics.perform(sa::String, basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings; printout::Bool=true)
        !(sa == "computation: CI")   &&   error("Unsupported keystring = $sa")

        # Determine the J^P symmetry blocks
        symmetries = Dict{JAC.LevelSymmetry,Int64}()
        for  csf in basis.csfs
            sym = LevelSymmetry(csf.J, csf.parity)
            if     haskey(symmetries, sym)    symmetries[sym] = symmetries[sym] + 1
            else   symmetries = Base.merge( symmetries, Dict( sym => 1 ) )
            end
        end

        # Test the total number of CSF
        NoCsf = 0
        for (sym,v) in symmetries   NoCsf = NoCsf + v   end

        # Calculate for each symmetry block the corresponding CI matrix, diagonalize it and append a Multiplet for this block
        multiplets = Multiplet[]
        for  (sym,v) in  symmetries
            matrix = compute("matrix: CI, J^P symmetry", sym, basis, nuclearModel, grid, settings; printout=printout)
            eigen  = Basics.diagonalize("matrix: Julia, eigfact", matrix)
            levels = Level[]
            for  ev = 1:length(eigen.values)
                # Construct the eigenvector with regard to the given basis (not w.r.t the symmetry block)
                evSym = eigen.vectors[ev];    vector = zeros( length(basis.csfs) );   ns = 0
                for  r = 1:length(basis.csfs)
                    if LevelSymmetry(basis.csfs[r].J, basis.csfs[r].parity) == sym    ns = ns + 1;   vector[r] = evSym[ns]   end
                end
                newlevel = Level( sym.J, AngularM64(sym.J.num//sym.J.den), sym.parity, 0, eigen.values[ev], 0., true, basis, vector )
                push!( levels, newlevel)
            end
            wa = Multiplet("noName", levels)
            push!( multiplets, wa)
        end

        # Merge all multiplets into a single one
        mp = Basics.merge(multiplets)
        mp = Basics.sortByEnergy(mp)

        # Display all level energies and energy splittings
        if  printout
            Basics.tabulate("multiplet: energies", mp)
            Basics.tabulate("multiplet: energy relative to immediately lower level",    mp)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", mp)
        end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary
            Basics.tabulate("multiplet: energies", mp, stream=iostream)
            Basics.tabulate("multiplet: energy relative to immediately lower level",    mp, stream=iostream)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", mp, stream=iostream)
        end

        return( mp )
    end

    """
    `Basics.perform("computation: CI for plasma", basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
        settings::AsfSettings, plasmaSettings::PlasmaShift.Settings; printout::Bool=true)`
        ... to  set-up and diagonalize from the given (SCF) basis the configuration-interaction matrix and to derive and
            display the level structure of the corresponding multiplet due to the given settings. Here, the CI matrix
            includes the modifications of the Hamiltonian due to the given plasmaSettings; a multiplet::Multiplet is returned.
    """
    function Basics.perform(sa::String, basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                    settings::AsfSettings, plasmaSettings::PlasmaShift.Settings; printout::Bool=true)
        !(sa == "computation: CI for plasma")   &&   error("Unsupported keystring = $sa")

        # Determine the J^P symmetry blocks
        symmetries = Dict{JAC.LevelSymmetry,Int64}()
        for  csf in basis.csfs
            sym = LevelSymmetry(csf.J, csf.parity)
            if     haskey(symmetries, sym)    symmetries[sym] = symmetries[sym] + 1
            else   symmetries = Base.merge( symmetries, Dict( sym => 1 ) )
            end
        end

        # Test the total number of CSF
        NoCsf = 0
        for (sym,v) in symmetries   NoCsf = NoCsf + v   end

        # Calculate for each symmetry block the corresponding CI matrix for the given plasma model and parameters as specified by
        # plasmaSettings, diagonalize it and append a Multiplet for this block
        multiplets = Multiplet[]
        for  (sym,v) in  symmetries
            matrix = compute("matrix: CI for plasma, J^P symmetry", sym, basis, nuclearModel, grid, settings, plasmaSettings; printout=printout)
            eigen  = Basics.diagonalize("matrix: Julia, eigfact", matrix)
            levels = Level[]
            for  ev = 1:length(eigen.values)
                # Construct the eigenvector with regard to the given basis (not w.r.t the symmetry block)
                evSym = eigen.vectors[ev];    vector = zeros( length(basis.csfs) );   ns = 0
                for  r = 1:length(basis.csfs)
                    if LevelSymmetry(basis.csfs[r].J, basis.csfs[r].parity) == sym    ns = ns + 1;   vector[r] = evSym[ns]   end
                end
                newlevel = Level( sym.J, AngularM64(sym.J.num//sym.J.den), sym.parity, 0, eigen.values[ev], 0., true, basis, vector )
                push!( levels, newlevel)
            end
            wa = Multiplet("noName", levels)
            push!( multiplets, wa)
        end

        # Merge all multiplets into a single one
        mp = Basics.merge(multiplets)
        mp = Basics.sortByEnergy(mp)

        # Display all level energies and energy splittings
        if  printout
            Basics.tabulate("multiplet: energies", mp)
            Basics.tabulate("multiplet: energy relative to immediately lower level",    mp)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", mp)
        end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary
            Basics.tabulate("multiplet: energies", mp, stream=iostream)
            Basics.tabulate("multiplet: energy relative to immediately lower level",    mp, stream=iostream)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", mp, stream=iostream)
        end

        return( mp )
    end

end # module


text:
using CosmoMMF
using Test

@testset "CosmoMMF.jl" begin
    # Write your tests here.

end


text:
export AbstractCorrector,
    AbstractCorrectorCache,
    CorrectorResult,
    cache,
    correct!,
    NewtonCorrector

# Interface
"""
Each corrector has to be a subtype of `AbstractCorrector`.
"""
abstract type AbstractCorrector end

"""
Each corrector has to construct a cache via `[cache](@ref)` which has
to be a subtype of `AbstractCorrectorCache`.
"""
abstract type AbstractCorrectorCache end

module CorrectorReturnCode
    import ..NewtonReturnCode

    @enum codes begin
        converged
        terminated
        terminated_no_approximate
        maximal_iterations
    end

    function code(c::NewtonReturnCode.codes)
        if c == NewtonReturnCode.converged
            converged
        elseif c == NewtonReturnCode.terminated
            terminated
        elseif c == NewtonReturnCode.terminated_no_approximate
            terminated_no_approximate
        else # maximal_iterations
            maximal_iterations
        end
    end
end

"""
    CorrectorResult{T}

Structure holding information about a `correct!` step. The fields are
* `converged::Bool` Indicating whether the correction was successfull
* `res::T` The residual of the final result.
* `iters::Int` The number of iterations used.
"""
struct CorrectorResult{T}
    return_code::CorrectorReturnCode.codes
    accuracy::T
    iters::Int
    ω₀::Float64
    ω::Float64
    norm_Δx₀::T
end

Base.show(io::IO, ::MIME"application/prs.juno.inline", r::CorrectorResult) = r
Base.show(io::IO, result::CorrectorResult) = print_fieldnames(io, result)

"""
    isconverged(result::Result)

Returns whether the correction was successfull.
"""
isconverged(result::CorrectorResult) = result.return_code == CorrectorReturnCode.converged

"""
    cache(::AbstractCorrector, ::HomotopyWithCache{M, N}, x, t)::AbstractCorrectorCache

Construct a cache to avoid allocations.
"""
cache(c::AbstractCorrector, args...) = throw(MethodError(cache, tuple(c, args...)))

"""
    correct!(xnext, ::AbstractCorrector, ::AbstractCorrectorCache, H::HomotopyWithCache, x, t, tol)::Result

Perform a correction step such that in the end `H(xnext,t) < tol`.
Returns a [`Result`](@ref).
"""
function correct! end

# Newton Corrector
"""
    NewtonCorrector(;simplified_last_step=true)

An ordinary Newton's method. If `simplified_last_step` is `true`, then for the last iteration
the previously Jacobian will be used. This uses an LU-factorization for square systems
and a QR-factorization for overdetermined.
"""
struct NewtonCorrector <: AbstractCorrector
    simplified_last_step::Bool
end

NewtonCorrector(;simplified_last_step=true) = NewtonCorrector(simplified_last_step)

struct NewtonCorrectorCache{FH<:FixedHomotopy, NC<:AbstractNewtonCache} <: AbstractCorrectorCache
    F::FH
    C::NC
end

function cache(::NewtonCorrector, H::HomotopyWithCache, x, t)
    F = FixedHomotopy(H, t)
    C = newton_cache(F, x)

    NewtonCorrectorCache(F, C)
end

@inline function correct!(out, alg::NewtonCorrector, cache::NewtonCorrectorCache,
                  H::HomotopyWithCache, x, t, norm, jacobian::Jacobian, tol::Float64,
                  maxiters::Int, step_size_model;
                  precision::PrecisionOption=PRECISION_FIXED_64,
                  update_jacobian_infos::Bool=false,
                  use_qr::Bool=false)
    cache.F.t = t
    result = newton!(out, cache.F, x, norm, cache.C, jacobian, tol, 1, maxiters,
        alg.simplified_last_step, update_jacobian_infos, use_qr, step_size_model.ω,
        step_size_model.expected_Δx₀, precision)
    CorrectorResult(result)
end

function CorrectorResult(R::NewtonResult)
    CorrectorResult(CorrectorReturnCode.code(R.return_code),
                    R.accuracy, R.iters, R.ω₀, R.ω, R.norm_Δx₀,)
end


text:
using SigmaRidgeRegression
using LinearAlgebra
using Test

tmp_block = BlockCovarianceDesign([IdentityCovarianceDesign(), IdentityCovarianceDesign(missing)], missing)

id = IdentityCovarianceDesign()
groups = GroupedFeatures([200;200])

@test set_groups(id, 400) == set_groups(id, groups)

instantiated_block = set_groups(tmp_block, groups)

bla = simulate_rotated_design(instantiated_block, 20)
@test size(bla) == (20,400)

instantiated_block.blocks[1]

spectrum(instantiated_block)
cov1 = SigmaRidgeRegression.UniformScalingCovarianceDesign(p=100, scaling=2.5)


text:
@testset "Problems" begin
    @testset "General Problem" begin
        F = equations(katsura(5))
        G = equations(cyclic(6))

        P1 = StartTargetInput(G, F)
        (PP1, start1) = problem_startsolutions(
            P1,
            [rand(ComplexF64, 6), rand(ComplexF64, 6)],
        )
        @test PP1 isa Problem{AffineTracking}
        @test length(start1) == 2

        P2 = TargetSystemInput(G)
        (PP2, start2) = problem_startsolutions(P2)
        @test PP2 isa Problem{AffineTracking}
        @test length(start2) == 720

        @test homotopy(PP2) isa AbstractHomotopy
        @test homvars(PP2) == nothing

        @polyvar x y z
        F = [x^2 + y^2 + 2 * z^2, x + y + 3z]
        P, _ = problem_startsolutions(TargetSystemInput(F), nothing, homvar = y)
        @test homvars(P) == (2,)

        P, _ = problem_startsolutions(TargetSystemInput(F), nothing)
        @test homvars(P) === nothing

        @test_throws ArgumentError problem_startsolutions(
            StartTargetInput([x^2 + y^2 + z^2, x^4 + y^4 + z^3], [x^3 + z^3, y^3 - z^3]),
            [rand(ComplexF64, 3)],
        )

        P, _ = problem_startsolutions(
            StartTargetInput([x^2 + y^2 + z^2, x^4 + y^4 + z^4], [x^3 + z^3, y^3 - z^3]),
            [rand(ComplexF64, 3)],
        )
        @test P isa Problem{ProjectiveTracking}
        @test homvars(P) === nothing

        P, _ = problem_startsolutions(
            StartTargetInput([x^2 + y^2 + z^2, x^4 + y^4 + z^2], [x^3 + z^3, y^3 - z^2]),
            [rand(ComplexF64, 3)],
        )
        @test P isa Problem{AffineTracking}

        # not supported right now
        @test_throws MethodError problem_startsolutions(
            StartTargetInput(
                SPSystem([x^2 + y^2 + z^2, x^4 + y^4 + z^2]),
                SPSystem([x^3 + z^3, y^3 - z^2]),
            ),
            [rand(ComplexF64, 3)],
        )

        P, starts = problem_startsolutions(HC.input_startsolutions(
            [x^2 + y^2 + 1, x^4 + y^4 + 1],
            [x^3 + 1, y^3 - 1],
            rand(ComplexF64, 3),
        )...)
        @test P isa Problem{AffineTracking}
        @test length(starts) == 1

        P, _ = problem_startsolutions(
            StartTargetInput([x^2 + y^2 + z^2, x^4 + y^4 + z^4], [x^3 + z^3, y^3 - z^3]),
            [rand(ComplexF64, 3)],
            homvar = z,
        )
        @test homvars(P) === (3,)

        F = FPSystem(homogenize(equations(cyclic(6))))
        P, startvals = problem_startsolutions(TargetSystemInput(F))
        @test P isa Problem{ProjectiveTracking}
        @test homvars(P) == nothing

        F = homogenize(equations(cyclic(6)))
        P, _ = problem_startsolutions(TargetSystemInput(F))
        @test P isa Problem{ProjectiveTracking}

        F = FPSystem(homogenize(equations(cyclic(6))))
        P, startvals = problem_startsolutions(TargetSystemInput(F), homvar = 5)
        @test P isa Problem{ProjectiveTracking}
        @test homvars(P) == (5,)
        @test length(startvals) == 720
    end

    @testset "PolyhedralProblem" begin
        @polyvar x y z
        F = [x^2 + y^2 + 2 * z^2, x + y + 3z]
        P, starts = problem_startsolutions(F; start_system = :polyhedral)
        @test P isa HC.PolyhedralProblem
        @test starts isa HC.PolyhedralStartSolutionsIterator

        @polyvar u v w
        F = [x^2 + y^2 + 2 * z^2, x + y + 3z] ∘ [u, v, w]
        P, starts = problem_startsolutions(F; start_system = :polyhedral)
        @test P isa HC.PolyhedralProblem
        @test starts isa HC.PolyhedralStartSolutionsIterator
    end

    @testset "multi-homogeneous" begin
        @polyvar z[1:6, 1:3]
        F = let
            p = [1, 1, 0]
            α = randn(5)
            a = randn(9)
            ## define the system of polynomials
            f = [z[i, :] ⋅ z[i, :] for i = 2:5]
            g = [z[i, :] ⋅ z[i+1, :] for i = 1:5]
            h = sum(a[i] .* (z[i, :] × z[i+1, :]) for i = 1:3) +
                sum(a[i+4] .* z[i, :] for i = 2:5)
            F′ = [f .- 1; g .- cos.(α); h .- p]
            ## assign values to z₁ and z₆
            [subs(f, z[1, :] => [1, 0, 0], z[6, :] => [1, 0, 0]) for f in F′]
        end
        group1 = [[z[2, :]; z[4, :]], [z[3, :]; z[5, :]]]
        group2 = [z[2, :], z[4, :], z[3, :], z[5, :]]

        @test bezout_number(F) == 1024
        @test bezout_number(F; variable_groups = group1) == 320
        @test bezout_number(F; variable_groups = group2) == 576

        _, starts = problem_startsolutions(F; variable_groups = group1)
        @test length(starts) == 320
        @test length(collect(starts)) == 320

        _, starts = problem_startsolutions(F; variable_groups = group2)
        @test length(starts) == 576
        @test length(collect(starts)) == 576
    end

    @testset "Overdetermined" begin
        @polyvar x y z
        prob, starts = problem_startsolutions([
            x - 2,
            y^2 + 3 * z^2,
            z^3 + x^3,
            z + x^2 + 3,
        ])
        @test prob isa HC.OverdeterminedProblem
        @test starts isa HC.TotalDegreeSolutionIterator
        @test starts.degrees == [3, 2, 2]

        # overdetermined, homogeneous
        prob, starts = problem_startsolutions([x - 2y, y^2 + 3 * x^2])
        @test prob isa HC.OverdeterminedProblem
        @test starts isa HC.TotalDegreeSolutionIterator
        @test starts.degrees == [2]

        prob, starts = problem_startsolutions([x - 2y, y^2 + 3 * x^2, x^3 + y^3])
        @test prob isa HC.OverdeterminedProblem
        @test starts isa HC.TotalDegreeSolutionIterator
        @test starts.degrees == [3]

        # Overdetermined, polyhedral
        prob, starts = problem_startsolutions(
            [x - 2, y^2 + 3 * z^2, z^3 + x^3, z + x^2 + 3];
            start_system = :polyhedral,
        )
        @test prob isa HC.OverdeterminedProblem
        @test prob.problem isa HC.PolyhedralProblem
        @test starts isa HC.PolyhedralStartSolutionsIterator

        # polyhedral needs to expand Composition
        @polyvar u v w
        prob, starts = problem_startsolutions(
            [x - 2, y^2 + 3 * z^2, z^3 + x^3, z + x^2 + 3] ∘ [u, v, w];
            start_system = :polyhedral,
        )
        @test prob isa HC.OverdeterminedProblem
        @test prob.problem isa HC.PolyhedralProblem
        @test starts isa HC.PolyhedralStartSolutionsIterator

        # Abstract Systems
        prob, starts = problem_startsolutions(FPSystem([
            x - 2z,
            y^2 + 3 * z^2,
            z^3 + x^3,
            z + x,
        ]))
        @test prob isa HC.OverdeterminedProblem{HC.ProjectiveTracking}
        @test starts isa HC.TotalDegreeSolutionIterator
        @test starts.degrees == [1, 2, 3]
    end

    @testset "Invalid input" begin
        @polyvar x y z
        @test_throws ArgumentError problem_startsolutions([x - 2y + 2, 0])
        @test_throws ArgumentError problem_startsolutions(
            [x - 2, y - 2],
            [x - 2, y - 2, y + 2],
            [[2, -3]],
        )
        @test_throws ArgumentError problem_startsolutions(
            [x - 2z, y^2 + 3z, z^3 + x^3],
            homvar = z,
        )

        # constant term
        @test_throws ArgumentError problem_startsolutions([subs(x + 2, x => 2), y^2 + 3x])
        @polyvar u v
        @test_throws ArgumentError problem_startsolutions([subs(x + 2, x => 2), y^2 + 3x] ∘
                                                          [u, v])
        @test_throws ArgumentError problem_startsolutions([x - 2z, y^2 + 3z], homvar = z)
        # test numerical homogeneous check fails
        @polyvar x y z
        G = FPSystem([x - 2z, y^2 + 3z])
        @test_throws ArgumentError problem_startsolutions(G, homvar = 3)
    end
end


text:
function test_element_wise_layer(backend::Backend, op::ElementWiseFunctorType, jl_op::Function, T, eps)
  println("-- Testing ElementWiseLayer{$op} on $(typeof(backend)){$T}...")

  tensor_dim = abs(rand(Int)) % 6 + 1
  dims = tuple((abs(rand(Int, tensor_dim)) % 8 + 1)...)
  println("    > $dims")

  NArg = get_num_args(op)
  inputs = [rand(T, dims) for i = 1:NArg]
  input_blobs = Blob[make_blob(backend, x) for x in inputs]
  diff_blobs = Blob[make_blob(backend, x) for x in inputs]

  layer = ElementWiseLayer(bottoms=Array(Symbol, NArg), tops=[:result], operation=op)
  state = setup(backend, layer, input_blobs, diff_blobs)

  forward(backend, state, input_blobs)
  expected_res = jl_op(inputs...)
  got_res = similar(inputs[1])
  copy!(got_res, state.blobs[1])
  @test all(abs(got_res - expected_res) .< eps)

  top_diff = rand(T, size(inputs[1]))
  copy!(state.blobs_diff[1], top_diff)

  backward(backend, state, input_blobs, diff_blobs)
  got_grads = map(diff_blobs) do blob
    arr = zeros(T, size(blob))
    copy!(arr, blob)
    arr
  end

  if jl_op == (+)
    for i = 1:length(got_grads)
      @test all(abs(got_grads[i] - top_diff) .< eps)
    end
  elseif jl_op == (-)
    @test all(abs(got_grads[1] - top_diff) .< eps)
    @test all(abs(got_grads[2] + top_diff) .< eps)
  elseif jl_op == (.*)
    @test all(abs(got_grads[1] - top_diff.*inputs[2]) .< eps)
    @test all(abs(got_grads[2] - top_diff.*inputs[1]) .< eps)
  elseif jl_op == (./)
    @test all(abs(got_grads[1] - top_diff./inputs[2]) .< eps)
    #@test all(abs(got_grads[2] + top_diff.*inputs[1]./(inputs[2].*inputs[2])) .< eps)
    @test all(abs(got_grads[2] + top_diff.*got_res./inputs[2]) .< eps)
  else
    error("Unknown operation $jl_op")
  end

  shutdown(backend, state)
end

function test_element_wise_layer(backend::Backend, T, eps)
  for (functor,jl_op) in ((ElementWiseFunctors.Add(), (+)),
                          (ElementWiseFunctors.Subtract(),(-)),
                          (ElementWiseFunctors.Multiply(), (.*)),
                          (ElementWiseFunctors.Divide(), (./)))
    test_element_wise_layer(backend, functor, jl_op, T, eps)
  end
end

function test_element_wise_layer(backend::Backend)
  test_element_wise_layer(backend, Float32, 1e-1)
  test_element_wise_layer(backend, Float64, 1e-5)
end

if test_cpu
  test_element_wise_layer(backend_cpu)
end
if test_gpu
  test_element_wise_layer(backend_gpu)
end


text:
#this file holds functions specific to acquiring the NCCS data

function downloadNCCSCSV(url::String, fileName::String, outPath::String,
  outName::String = fileName; validate::Bool=true,
  refreshData::Bool=true, maxAttempts::Int = 12)::Nothing

  if !refreshData && isfile("$outPath\\$(outName).csv.gz")# if the file exists
    if validate #validate if desired
    try #validate the file
        validateSingleTypeCSV("$outPath\\$outName.csv.gz")
        println("Download Validated: $outName")
        return nothing
      catch err
        #do nothing and do not return anything
      end
    else
      return nothing
    end
  end

  attempt::Int = 0
  retry::Bool = true
  while attempt < maxAttempts && retry #loops and tries to download all the files
    attempt+=1
    retry = false
    #download the file
    oStream::GZipStream = GZip.open("$outPath\\$outName.csv.gz", "w")
    #b::IOBuffer = IOBuffer()
    #weird error if using the transcode codec

    #acquire the data
    HTTP.open(io::IO->
      while !eof(io)
        write(oStream, readavailable(io))
      end,
      "GET", "$url/$fileName", retries=4)#, connection_limit=16,
      #pipeline_limit=32, readtimeout=120)

    #bString::String = String(take!(b))
    #TODO: Might need to uncomment #bString = replace(bString, "\0"=>"")
    #write(oStream, bString)
    GZip.close(oStream)

    #check if its a valid file
    try
      validateSingleTypeCSV("$outPath\\$outName.csv.gz")
      println("Download Validated: $outName")
    catch err
      @warn "DL Error with $outName. $(
        attempt<maxAttempts ? "Retrying ($attempt/$maxAttempts)" : "Max attempts Reached for file $outName.csv.")
        \nErrorMessage:\n$err\nurl:$url/$fileName"
      retry = true
    end

  end

  return nothing
end

#vectorized version (useful for parallelization)
function downloadNCCSCSVs(urls::Vector{String}, fileNames::Vector{String},
  outPath::String, outNames::Vector{String}, mode::Symbol = :serial,
  ;delay::Float64 = 0.0,
  validate::Bool = true, refreshData::Bool=true)::Nothing

  #short-circuit for parallel operation
  if mode == :parallel
    return downloadNCCSCSVsPar(urls, fileNames, outPath, outNames,
    delay=delay, validate=validate, refreshData=refreshData)
  end

  for i::Int ∈ 1:length(urls)

    if delay ≥ .001 #allows for a delay to avoid slamming the site
      sleep(delay)
    end

    downloadNCCSCSV(urls[i], fileNames[i], outPath, outNames[i],
      validate=validate, refreshData=refreshData)
  end

  return nothing
end

#useful for parallelization
function downloadNCCSCSVsPar(urls::Vector{String}, fileNames::Vector{String},
  outPath::String, outNames::Vector{String};
  delay::Float64 = 0.0, validate::Bool = true, refreshData::Bool=true)::Nothing

  numFiles::Int = length(urls)
  pids = workers()
  np::Int = length(pids)

  #divide up the work
  assignments::Vector{Vector{Int}} = ((i::Int)->Vector{Int}()).(1:np)
  for i::Int ∈ 1:numFiles
    push!(assignments[i%np+1], i) #assigns a file to a pid
  end

  #set up a vector for pushing and receiving assignments
  futures::Vector{Future} = Vector{Future}()
  sizehint!(futures, np)

  for i::Int ∈ 1:np
    if length(assignments[i]) > 0
      push!(futures,
        (@spawnat pids[i] downloadNCCSCSVs(urls[assignments[i]],
          fileNames[assignments[i]], outPath,
          outNames[assignments[i]], delay=delay, validate=validate,
          refreshData=refreshData)))
    end
  end

  #get the work
  (fetch).(futures)

  return nothing
end

function readNCCSCSVsIndex(;inPath::String = DATA_PATH, filesIndexName::String = FILES_INDEX_NAME)::DataFrame
  return CSV.read("$inPath\\$filesIndexName.csv", categorical=false)
  #return load("$inPath\\$filesIndexName.csv") |> DataFrame
end

function downloadNCCSCSVs(outPath::String; mode::Symbol=:parallel,
  delay::Float64=0.0, validate::Bool = true, refreshData::Bool=true)

  #read in the metadata
  fileDF::DataFrame = readNCCSCSVsIndex()
  numFiles::Int = size(fileDF,1)

  #build the urls
  urls::Vector{String} = ((urlStem::String, urlYear::Int)->
    "$urlStem/$urlYear").(fileDF[:urlStems], fileDF[:urlYears])
  fileNames::Vector{String} =  fileDF[:fileNames]
  outNames::Vector{String} =  fileDF[:outNames]

  return downloadNCCSCSVs(urls, fileNames, outPath, outNames, mode, delay=delay,
    validate=validate, refreshData=refreshData)

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text: # concept here is to show how to implement a custom encode_f using the case of the highest occupied and lowest unoccupied atomic orbitals and pulling their energies from the features table