text:
"""
A `MappedDomain` represents the mapping of a domain.

The map of a domain `d` under the mapping `y=f(x)` consists of all points `f(x)`
with `x ∈ d`. The characteristic function of a mapped domain is defined in
terms of the inverse map `g = inverse(f)`:
```
x ∈ m ⟺ g(x) ∈ d
```
"""
abstract type AbstractMappedDomain{T} <: SimpleLazyDomain{T} end

superdomain(d::AbstractMappedDomain) = d.domain

const MappedVectorDomain{T} = AbstractMappedDomain{Vector{T}}

tointernalpoint(d::AbstractMappedDomain, x) = inverse_map(d, x)
toexternalpoint(d::AbstractMappedDomain, y) = forward_map(d, y)

canonicaldomain(d::AbstractMappedDomain) = canonicaldomain(superdomain(d))
mapfrom_canonical(d::AbstractMappedDomain) = forward_map(d) ∘ mapfrom_canonical(superdomain(d))
mapto_canonical(d::AbstractMappedDomain) = mapto_canonical(superdomain(d)) ∘ inverse_map(d)

# TODO: check whether the map alters the dimension
dimension(d::MappedVectorDomain) = dimension(superdomain(d))

# TODO: check whether the map affects these properties
isempty(d::AbstractMappedDomain) = isempty(superdomain(d))
isopenset(d::AbstractMappedDomain) = isopenset(superdomain(d))
isclosedset(d::AbstractMappedDomain) = isclosedset(superdomain(d))

## I/O functionality

show(io::IO, mime::MIME"text/plain", d::AbstractMappedDomain) = composite_show(io, mime, d)
Display.displaystencil(d::AbstractMappedDomain) =
    map_stencil_broadcast(forward_map(d), superdomain(d))
Display.object_parentheses(d::AbstractMappedDomain) =
    Display.object_parentheses(forward_map(d))
Display.stencil_parentheses(d::AbstractMappedDomain) =
    Display.stencil_parentheses(forward_map(d))

"A `MappedDomain` stores the inverse map of a mapped domain."
struct MappedDomain{T,F,D} <: AbstractMappedDomain{T}
    invmap  ::  F
    domain  ::  D
end

# In the constructor, we have to decide which T to use for the MappedDomain.
# - we don't know anything about invmap: deduce T from the given domain
MappedDomain(invmap, domain::Domain{T}) where {T} = MappedDomain{T}(invmap, domain)
# - if the map is a Map{T}, use that T for the MappedDomain
MappedDomain(invmap::Map{T}, domain::Domain{S}) where {S,T} = MappedDomain{T}(invmap, domain)
# - if T is given in the constructor, by all means we use that
MappedDomain{T}(invmap, domain::Domain) where {T} =
    MappedDomain{T,typeof(invmap),typeof(domain)}(invmap, domain)
# - in that case, if the map is a Map{S}, make sure that S matches T
MappedDomain{T}(invmap::Map{T}, domain::Domain) where {T} =
    MappedDomain{T,typeof(invmap),typeof(domain)}(invmap, domain)
MappedDomain{T}(invmap::Map{S}, domain::Domain) where {S,T} =
    MappedDomain{T}(convert(Map{T}, invmap), domain)

similardomain(d::MappedDomain, ::Type{T}) where {T} =
    MappedDomain{T}(d.invmap, d.domain)

forward_map(d::MappedDomain) = rightinverse(d.invmap)
forward_map(d::MappedDomain, x) = rightinverse(d.invmap, x)

inverse_map(d::MappedDomain) = d.invmap
inverse_map(d::MappedDomain, y) = d.invmap(y)

"Map a domain with the inverse of the given map"
map_domain(map, domain::Domain) = _map_domain(map, domain)

# Fallback: we don't know anything about map, just try to invert
_map_domain(map, domain) = mapped_domain(inverse(map), domain)
_map_domain(map::Map{T}, domain::Domain{T}) where {T} =
    mapped_domain(inverse(map), domain)
# If map is a Map{T}, then verify and if necessary update T
function _map_domain(map::Map, domain)
    U = codomaintype(map, eltype(domain))
    if U == Union{}
        error("incompatible types of $(map) and $(domain)")
    end
    mapped_domain(inverse(convert(Map{U}, map)), convert(Domain{U}, domain))
end

==(a::MappedDomain, b::MappedDomain) = (a.invmap == b.invmap) && (superdomain(a) == superdomain(b))

"Make a mapped domain with the given inverse map"
mapped_domain(invmap, domain::Domain) = _mapped_domain(invmap, domain)

# We face the same task as in the constructor: attempt to identify T
# Here, we are more flexible, and attempt to do more conversions. We assume
# that users invoking the MappedDomain constructor know what they are doing,
# but users invoking mapped_domain just expect it to work.

# - we don't know anything about invmap, just pass it on
_mapped_domain(invmap, domain) = MappedDomain(invmap, domain)
# - invmap is a Map{T}: its codomaintype should match the eltype of the domain
# -- first, update the numtype
_mapped_domain(invmap::Map, domain) =
    _mapped_domain(invmap, domain, promote_type(numtype(invmap),numtype(domain)))
_mapped_domain(invmap::Map{T}, domain::Domain{S}, ::Type{U}) where {S,T,U} =
    _mapped_domain2(convert_numtype(invmap,U), convert_numtype(domain,U))
# -- then, ensure the codomaintype of the map equals the element type of the domain
_mapped_domain2(invmap, domain) = _mapped_domain2(invmap, domain, codomaintype(invmap), eltype(domain))
# --- it's okay
_mapped_domain2(invmap, domain, ::Type{T}, ::Type{T}) where {T} =
    MappedDomain(invmap, domain)
# --- it's not okay: attempt to convert the map (triggers e.g. when combining a scalar
#        LinearMap with a vector domain)
_mapped_domain2(invmap, domain, ::Type{S}, ::Type{T}) where {S,T} =
    MappedDomain(convert(Map{T}, invmap), domain)

# TODO: deprecate this syntax, it is confusing
(∘)(domain::Domain, invmap::Function) = mapped_domain(invmap, domain)
(∘)(domain::Domain, invmap::AbstractMap) = mapped_domain(invmap, domain)

# Avoid nested mapping domains, construct a composite map instead
# This assumes that the map types can be combined using \circ
mapped_domain(invmap, d::MappedDomain) = mapped_domain(inverse_map(d) ∘ invmap, superdomain(d))

boundary(d::MappedDomain) = _boundary(d, superdomain(d), inverse_map(d))
_boundary(d::MappedDomain, superdomain, invmap) = MappedDomain(invmap, boundary(superdomain))
interior(d::MappedDomain) = _interior(d, superdomain(d), inverse_map(d))
_interior(d::MappedDomain, superdomain, invmap) = MappedDomain(invmap, interior(superdomain))
closure(d::MappedDomain) = _closure(d, superdomain(d), inverse_map(d))
_closure(d::MappedDomain, superdomain, invmap) = MappedDomain(invmap, closure(superdomain))

convert(::Type{MappedDomain}, d::Domain{T}) where {T} =
    MappedDomain{T}(mapto_canonical(d), canonicaldomain(d))
convert(::Type{MappedDomain{T}}, d::Domain) where {T} =
    MappedDomain{T}(mapto_canonical(d), canonicaldomain(d))

"A `ParametricDomain` stores the forward map of a mapped domain."
struct ParametricDomain{T,F,D} <: AbstractMappedDomain{T}
    fmap    ::  F
    domain  ::  D
end

ParametricDomain(fmap, domain::Domain) = ParametricDomain{codomaintype(fmap)}(fmap, domain)
ParametricDomain{T}(fmap, domain::Domain) where {T} =
    ParametricDomain{T,typeof(fmap),typeof(domain)}(fmap, domain)

similardomain(d::ParametricDomain, ::Type{T}) where {T} =
    ParametricDomain{T}(d.fmap, d.domain)

forward_map(d::ParametricDomain) = d.fmap
forward_map(d::ParametricDomain, x) = d.fmap(x)

inverse_map(d::ParametricDomain) = leftinverse(d.fmap)
inverse_map(d::ParametricDomain, y) = leftinverse(d.fmap, y)

"Return the domain that results from mapping the given domain."
parametric_domain(fmap, domain::Domain) = ParametricDomain(fmap, domain)

boundary(d::ParametricDomain) = _boundary(d, superdomain(d), forward_map(d))
_boundary(d::ParametricDomain, superdomain, fmap) = ParametricDomain(fmap, boundary(superdomain))
interior(d::ParametricDomain) = _interior(d, superdomain(d), forward_map(d))
_interior(d::ParametricDomain, superdomain, fmap) = ParametricDomain(fmap, interior(superdomain))
closure(d::ParametricDomain) = _closure(d, superdomain(d), forward_map(d))
_closure(d::ParametricDomain, superdomain, fmap) = ParametricDomain(fmap, closure(superdomain))


text:
include("image_geom.jl")
include("ellipse_im.jl")

export ir_dump
export embed
export image_geom
export MIRT_image_geom
export ellipse_im

ext: jl


text:
##############################################################################
#
# LSA
#
##############################################################################

lsa(dtm::DocumentTermMatrix) = svd(tf_idf(dtm))
lsa(crps::Corpus) = svd(tf_idf(DocumentTermMatrix(crps)))


text:
################################################################################
#                    Poly - the not so primitive, primitive                    #
################################################################################

"""
Special method for polys so we don't fall back to atomic meshes, which are much more
complex and slower to draw than standard paths with single color.
"""
function draw_plot(scene::Scene, screen::CairoScreen, poly::Poly)
    # dispatch on input arguments to poly to use smarter drawing methods than
    # meshes if possible
    draw_poly(scene, screen, poly, to_value.(poly.input_args)...)
end

"""
Fallback method for args without special treatment.
"""
function draw_poly(scene::Scene, screen::CairoScreen, poly, args...)
    draw_poly_as_mesh(scene, screen, poly)
end

function draw_poly_as_mesh(scene, screen, poly)
    draw_plot(scene, screen, poly.plots[1])
    draw_plot(scene, screen, poly.plots[2])
end

function draw_poly(scene::Scene, screen::CairoScreen, poly, points::Vector{<:Point2})

    # in the rare case of per-vertex colors redirect to mesh drawing
    if poly.color[] isa Array
        draw_poly_as_mesh(scene, screen, poly)
        return
    end

    model = poly.model[]
    points = project_position.(Ref(scene), points, Ref(model))
    Cairo.move_to(screen.context, points[1]...)
    for p in points[2:end]
        Cairo.line_to(screen.context, p...)
    end
    Cairo.close_path(screen.context)
    Cairo.set_source_rgba(screen.context, rgbatuple(to_color(poly.color[]))...)
    Cairo.fill_preserve(screen.context)
    Cairo.set_source_rgba(screen.context, rgbatuple(to_color(poly.strokecolor[]))...)
    Cairo.set_line_width(screen.context, poly.strokewidth[])
    Cairo.stroke(screen.context)
end

function draw_poly(scene::Scene, screen::CairoScreen, poly, points_list::Vector{<:Vector{<:Point2}})
    for points in points_list
        draw_poly(scene, screen, poly, points)
    end
end

draw_poly(scene::Scene, screen::CairoScreen, poly, rect::Rect2D) = draw_poly(scene, screen, poly, [rect])

function draw_poly(scene::Scene, screen::CairoScreen, poly, rects::Vector{<:Rect2D})
    model = poly.model[]
    projected_rects = project_rect.(Ref(scene), rects, Ref(model))

    color = poly.color[]
    if color isa AbstractArray{<:Number}
        color = numbers_to_colors(color, poly)
    elseif color isa String
        # string is erroneously broadcasted as chars otherwise
        color = to_color(color)
    end
    strokecolor = poly.strokecolor[]
    if strokecolor isa AbstractArray{<:Number}
        strokecolor = numbers_to_colors(strokecolor, poly)
    elseif strokecolor isa String
        # string is erroneously broadcasted as chars otherwise
        strokecolor = to_color(strokecolor)
    end

    broadcast_foreach(projected_rects, color, strokecolor, poly.strokewidth[]) do r, c, sc, sw
        Cairo.rectangle(screen.context, origin(r)..., widths(r)...)
        Cairo.set_source_rgba(screen.context, rgbatuple(to_color(c))...)
        Cairo.fill_preserve(screen.context)
        Cairo.set_source_rgba(screen.context, rgbatuple(to_color(sc))...)
        Cairo.set_line_width(screen.context, sw)
        Cairo.stroke(screen.context)
    end
end


text:
using Plots, RDatasets, DataFrames, StatPlots, Query

gr()

iris = dataset("datasets", "iris")

function dfmean(df::DataFrame)
    return [mean(i) for i in df.columns]
end

function dfcov(df)
    cols = df.columns
    ncols = size(cols)[1]
    ndatas = size(cols[1])[1]
    array = zeros(ncols, ncols)
    mean = dfmean(df)
    for i in 1:ncols
        for j in 1:ncols
            array[i, j] = 1/ndatas * (cols[i] .- mean[i])' * (cols[j] .- mean[j])
        end
    end
    return array
end

function whitening(df)
    groups = Array(df[:Species])
    petals = df[[:PetalLength, :PetalWidth]]
    colnames = names(petals)

    values = Array(petals)
    mu = reshape(dfmean(petals), (1, 2))

    coviris = dfcov(petals)
    S = -swap(eigvecs(coviris))

    Lambda = S' * coviris * S
    L = zeros(size(Lambda))
    for i in 1:size(L)[1]
        L[i, i] = sqrt(Lambda[i, i])
    end
    u = (values .- mu) * S * inv(L)

    tf_df = DataFrame()
    tf_df[colnames[1]] = u[:, 1]
    tf_df[colnames[2]] = u[:, 2]
    tf_df[:Species] = groups
    return tf_df
end

function swap(array)
    temp = zeros(2, 2)
    temp[:, 1] = array[:, 2]
    temp[:, 2] = array[:, 1]

    return temp
end

@df iris scatter(:PetalLength, :PetalWidth, group=:Species,
                 m=(0.5, [:+ :h :star7], 12), bg=RGB(.2,.2,.2))

scatter!(xlabel="Petal Length", ylabel="Petal Width", title="No Whitening",
         xlims=(1, 7), ylims=(-1, 4))
savefig("../figures/no_whitening.png")

whitened = whitening(iris)

@df whitened scatter(:PetalLength, :PetalWidth, group=:Species,
                         m=(0.5, [:+ :h :star7], 12), bg=RGB(.2,.2,.2))
scatter!(xlabel="Petal Length", ylabel="Petal Width", title="Whitened",
         xlims=(-3, 3), ylims=(-3, 3))
savefig("../figures/whitened.png")


text:
    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
            basis     = Basics.performSCF(computation.configs, nModel, computation.grid, computation.asfSettings)
            multiplet = Basics.performCI( 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
            for settings  in computation.propertySettings
                if      typeof(settings) == Einstein.Settings
                    outcome = Einstein.computeLines(multiplet,        computation.grid, settings)
                    if output    results = Base.merge( results, Dict("Einstein lines:" => outcome) )                  end
                    #
                elseif  typeof(settings) == Hfs.Settings    && settings.calcIJFexpansion
                    outcome = Hfs.computeHyperfineMultiplet(multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("IJF multiplet:" => outcome) )                   end
                    #
                elseif  typeof(settings) == Hfs.Settings
                    outcome = Hfs.computeOutcomes(multiplet, nModel,  computation.grid, settings)
                    if output    results = Base.merge( results, Dict("HFS outcomes:" => outcome) )                    end
                    #
                elseif  typeof(settings) == LandeZeeman.Settings
                    outcome = LandeZeeman.computeOutcomes(multiplet, nModel,  computation.grid, settings)
                    if output    results = Base.merge( results, Dict("Zeeman parameter outcomes:" => outcome) )       end
                    #
                elseif  typeof(settings) == IsotopeShift.Settings
                    outcome = IsotopeShift.computeOutcomes(multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("Isotope parameter outcomes:" => outcome) )      end
                    #
                elseif  typeof(settings) == AlphaVariation.Settings
                    outcome = AlphaVariation.computeOutcomes(multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("alpha variation parameter outcomes:" => outcome) )      end
                    #
                elseif  typeof(settings) == FormFactor.Settings
                    outcome = FormFactor.computeOutcomes(multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("Form factor outcomes:" => outcome) )            end
                    #
                elseif  typeof(settings) == DecayYield.Settings
                    outcome = DecayYield.computeOutcomes(computation.configs, computation.asfSettings,
                                                            multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("Fluorescence and AutoIonization yield outcomes:" => outcome) )   end
                    #
                elseif  typeof(settings) == MultipolePolarizibility.Settings
                    outcome = MultipolePolarizibility.computeOutcomes(multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("Polarizibility outcomes:" => outcome) )         end
                    #
                elseif  typeof(settings) == ReducedDensityMatrix.Settings
                    outcome = ReducedDensityMatrix.computeOutcomes(multiplet, nModel, computation.grid, settings)
                    if output    results = Base.merge( results, Dict("RDM outcomes:" => outcome) )         end
                    #
                elseif  typeof(settings) == PlasmaShift.Settings
                    outcome = PlasmaShift.computeOutcomes(multiplet, nModel, computation.grid, computation.asfSettings, settings)
                    if output    results = Base.merge( results, Dict("Plasma shift outcomes:" => outcome) )            end
                end
            end

        # Evaluate processes that need special SCF and CI procedures
        elseif  typeof(computation.processSettings)  in [PlasmaShift.AugerSettings, PlasmaShift.PhotoSettings]
            pSettings        = computation.processSettings
            plasmaSettings   = PlasmaShift.Settings(pSettings.plasmaModel, pSettings.lambdaDebye, pSettings.ionSphereR0, pSettings.NoBoundElectrons)
            initialBasis     = Basics.performSCF(computation.initialConfigs, nModel, computation.grid, computation.initialAsfSettings)
            initialMultiplet = perform("computation: CI for plasma",  initialBasis, nModel, computation.grid, computation.initialAsfSettings, plasmaSettings)
            finalBasis       = Basics.performSCF(computation.finalConfigs, nModel, computation.grid, computation.finalAsfSettings)
            finalMultiplet   = perform("computation: CI for plasma",  finalBasis, nModel, computation.grid, computation.finalAsfSettings, plasmaSettings)
            #
            if      typeof(computation.processSettings)  == PlasmaShift.AugerSettings
                outcome = AutoIonization.computeLinesPlasma(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("AutoIonization lines in plasma:" => outcome) )         end
            elseif  typeof(computation.processSettings)  == PlasmaShift.PhotoSettings
                outcome = 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     = Basics.performSCF(computation.initialConfigs, nModel, computation.grid, computation.initialAsfSettings)
            initialMultiplet = Basics.performCI( initialBasis, nModel, computation.grid, computation.initialAsfSettings)
            LSjj.expandLevelsIntoLS(initialMultiplet, computation.initialAsfSettings.jjLS)
            finalBasis       = Basics.performSCF(computation.finalConfigs, nModel, computation.grid, computation.finalAsfSettings)
            finalMultiplet   = Basics.performCI(   finalBasis, nModel, computation.grid, computation.finalAsfSettings)
            LSjj.expandLevelsIntoLS(finalMultiplet, computation.finalAsfSettings.jjLS)
            #
            if typeof(computation.processSettings) in [PhotoExcitationFluores.Settings, PhotoExcitationAutoion.Settings, PhotoIonizationFluores.Settings,
                                                       PhotoIonizationAutoion.Settings, ImpactExcitationAutoion.Settings, Dielectronic.Settings]
                intermediateBasis     = Basics.performSCF(computation.intermediateConfigs, nModel, computation.grid, computation.intermediateAsfSettings)
                intermediateMultiplet = Basics.performCI( intermediateBasis, nModel, computation.grid, computation.intermediateAsfSettings)
            end
            #
            if      typeof(computation.processSettings) == AutoIonization.Settings
                outcome = AutoIonization.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("AutoIonization lines:" => outcome) )                  end
            elseif  typeof(computation.processSettings) == RayleighCompton.Settings
                outcome = RayleighCompton.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Rayleigh-Compton lines:" => outcome) )                end
            elseif  typeof(computation.processSettings) == DoubleAutoIonization.Settings
                outcome = DoubleAutoIonization.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Double-Auger lines:" => outcome) )                    end
            elseif  typeof(computation.processSettings) == Dielectronic.Settings
                outcome = Dielectronic.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, nModel,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("dielectronic recombination pathways:" => outcome) )   end
            elseif  typeof(computation.processSettings) == MultiPhotonIonization.Settings
                outcome = MultiPhotonDeExcitation.computeLines(computation.processSettings.process,
                                                               finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("multi-photon-de-excitation lines:" => outcome) )      end
            elseif  typeof(computation.processSettings) == PhotoIonization.Settings
                outcome = PhotoIonization.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photoionization lines:" => outcome) )                 end
            elseif  typeof(computation.processSettings) == PhotoDoubleIonization.Settings
                outcome = PhotoDoubleIonization.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Single-photon double-ionization lines:" => outcome) ) end
            elseif  typeof(computation.processSettings) == PhotoExcitation.Settings
                outcome = PhotoExcitation.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-excitation lines:" => outcome) )                end
            elseif  typeof(computation.processSettings) == PhotoExcitationAutoion.Settings
                outcome = PhotoExcitationAutoion.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, nModel,
                                                                 computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-excitation-autoionization pathways:" => outcome) )      end
            elseif  typeof(computation.processSettings) == PhotoExcitationFluores.Settings
                outcome = PhotoExcitationFluores.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                 computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-excitation-fluorescence pathways:" => outcome) )        end
            elseif  typeof(computation.processSettings) == PhotoEmission.Settings
                outcome = PhotoEmission.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("radiative lines:" => outcome) )                       end
            elseif  typeof(computation.processSettings) == RadiativeAuger.Settings
                outcome = RadiativeAuger.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("radiative Auger sharings:" => outcome) )              end
            elseif  typeof(computation.processSettings) == PhotoRecombination.Settings
                outcome = PhotoRecombination.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo recombination lines:" => outcome) )             end
            elseif  typeof(computation.processSettings) == ImpactExcitation.Settings
                outcome = ImpactExcitation.computeLines(finalMultiplet, initialMultiplet, nModel, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("impact-excitation lines:" => outcome) )               end
                #
                #
            elseif  typeof(computation.processSettings) == PairA1P()
                outcome = PairAnnihilation1Photon.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("pair-annihilation-1-photon lines:" => outcome) )       end
            elseif  typeof(computation.processSettings) == Coulex()
                outcome = CoulombExcitation.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Coulomb excitation lines:" => outcome) )               end
            elseif  typeof(computation.processSettings) == Coulion()
                outcome = CoulombIonization.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("Coulomb ionization lines:" => outcome) )               end
            elseif  typeof(computation.processSettings) == PhotoIonAuto()
                outcome = PhotoIonizationAutoion.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-ionization-autoionization pathways:" => outcome) )      end
            elseif  typeof(computation.processSettings) == PhotoIonFluor()
                outcome = PhotoIonizationFluores.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("photo-ionizatiton-fluorescence pathways:" => outcome) )        end
            elseif  typeof(computation.processSettings) == ImpactExcAuto()
                outcome = ImpactExcitationAutoion.computePathways(finalMultiplet, intermediateMultiplet, initialMultiplet,
                                                                    computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("impact-excitation-autoionization pathways:" => outcome) )     end
            elseif  typeof(computation.processSettings) == MultiPI()
                outcome = MultiPhotonIonization.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("multi-photon single ionization:" => outcome) )        end
            elseif  typeof(computation.processSettings) == MultiPDI()
                outcome = MultiPhotonDoubleIon.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("multi-photon double ionization:" => outcome) )        end
            elseif  typeof(computation.processSettings) == InternalConv()
                outcome = InternalConversion.computeLines(finalMultiplet, initialMultiplet, computation.grid, computation.processSettings)
                if output    results = Base.merge( results, Dict("internal conversion lines:" => 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::Bool=true, outputToFile::Bool=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, outputToFile::Bool=true)
        Cascade.perform(comp.scheme, comp::Cascade.Computation, output=output, outputToFile=outputToFile)
    end

    """
    `Basics.perform(comp::Cascade.Simulation)`
        ... 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.Simulation; 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.Simulation; output::Bool=false)
        Cascade.perform(comp::Cascade.Simulation, output=output)
    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")
        return( Basics.performCI(configs::Array{Configuration,1}, initalOrbitals::Dict{Subshell, Orbital}, nuclearModel::Nuclear.Model,
                                 grid::Radial.Grid, settings::AsfSettings; printout=printout) )
    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")
        return( Basics.performCI(basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                                 settings::AsfSettings, plasmaSettings::PlasmaShift.Settings; printout=printout) )
    end

    """
    `Basics.performSCF(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.performSCF(configs::Array{Configuration,1}, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings;
                               printout::Bool=true)
        if  printout    println("\n... in performSCF ...")    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(">> include ", 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

        wa = Bsplines.generatePrimitives(grid)

        # Generate start orbitals
        if  typeof(settings.startScfFrom) == StartFromHydrogenic
            if  printout   println("Start SCF process with hydrogenic orbitals.")   end
            # Generate start orbitals for the SCF field by using B-splines
            ##x orbitals  = Bsplines.generateOrbitalsHydrogenic(waL, nsL, waS, nsS, nuclearModel, subshellList)
            orbitals  = Bsplines.generateOrbitalsHydrogenic(wa, nuclearModel, subshellList; printout=printout)
        elseif  typeof(settings.startScfFrom) == StartFromPrevious
            if  printout   println("Start SCF process from given list of orbitals.energy")    end
            # Taking starting orbitals for the given dictionary; non-relativistic orbitals with a proper nuclear charge
            # are adapted if no orbital is found
            orbitals = Dict{Subshell, Orbital}()
            for  subsh in subshellList
                if  haskey(settings.startScfFrom.orbitals, subsh)
                    orb = settings.startScfFrom.orbitals[subsh]
                    orbitals = Base.merge( orbitals, Dict( subsh => orb) )
                else
                    println("Start orbitals do not contain an Orbital for subshell $subsh ")
                    orb      = HydrogenicIon.radialOrbital(subsh, nuclearModel.Z, grid)
                    orb      = Orbital(orb.subshell, orb.isBound, true, orb.energy, orb.P, orb.Q, orb.Pprime, orb.Qprime, Radial.Grid())
                    orbitals = Base.merge( orbitals, Dict( subsh => orb) )
                end
            end
        else  error("stop a")
        end

        # Perform the SCF computations for the orbitals due to the given settings
        basis    = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, orbitals)
        if   settings.scField in [Basics.DFSField(), Basics.HSField()]
            basis = Bsplines.solveSelfConsistentMeanField(wa, nuclearModel, basis, settings; printout=printout)
        elseif   settings.scField in [Basics.NuclearField()]  && settings.startScfFrom == StartFromHydrogenic()
            # Return the basis as already generated.
        elseif   settings.scField in [Basics.ALField()]
            basis = Bsplines.solveSelfConsistentALField(wa, nuclearModel, basis, settings; printout=printout)
        elseif   typeof(settings.scField) == Basics.EOLField
            basis = Bsplines.solveSelfConsistentEOLField(wa, nuclearModel, basis, settings; printout=printout)
        else  error("stop b")
        end

        return( basis )
    end

    """
    `Basics.performSCF(basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid, frozenShells::Array{Shell,1},
                       settings::AtomicState.RasSettings; printout::Bool=true)`
        ... to generate an atomic basis and to compute the self-consistent field (SCF) for this basis due to the given settings;
            all subshells due to the list of frozenShells are kept fixed. A newBasis::Basis is returned.
    """
    function Basics.performSCF(basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid, frozenShells::Array{Shell,1},
                               settings::AtomicState.RasSettings; printout::Bool=true)
        if  printout    println("\n... in performSCF['for RAS step'] ...")    end

        # Determine the list of frozen subshells
        frozenSubshells = Subshell[]
        for  sh in frozenShells
            if  sh.l == 0    push!( frozenSubshells, Subshell(sh.n, -1))
            else             push!( frozenSubshells, Subshell(sh.n, sh.l));    push!( frozenSubshells, Subshell(sh.n, -sh.l -1))
            end
        end

        wa = Bsplines.generatePrimitives(grid)
        # Specify the AsfSettings for the given RAS step
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),
                                  settings.maxIterationsScf, settings.accuracyScf, Subshell[], frozenSubshells,
                                  CoulombInteraction(), NoneQed(), LSjjSettings(false), settings.levelSelectionCI )

        basis = Bsplines.solveSelfConsistentMeanField(wa, nuclearModel, basis, asfSettings; printout=printout)

        return( basis )
    end

    """
    `Basics.performCI(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.performCI(configs::Array{Configuration,1}, initalOrbitals::Dict{Subshell, Orbital}, nuclearModel::Nuclear.Model,
                              grid::Radial.Grid, settings::AsfSettings; printout::Bool=true)
        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(">> include ", 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 = Nuclear.nuclearPotential(nuclearModel, grid)

        keep = true
        InteractionStrength.XL_Coulomb_reset_storage(keep)
        #
        n = length(csfList);    matrix = zeros(Float64, n, n)
        for  r = 1:n
            ##x wa = compute("angular coefficients: e-e, Ratip2013", basis.csfs[r], basis.csfs[r])
            # Calculate the spin-angular coefficients
            if  Defaults.saRatip()
                waR = compute("angular coefficients: e-e, Ratip2013", basis.csfs[r], basis.csfs[r])
                wa  = waR
            end
            if  Defaults.saGG()
                subshellList = basis.subshells
                opa  = SpinAngular.OneParticleOperator(0, plus, true)
                waG1 = SpinAngular.computeCoefficients(opa, basis.csfs[r], basis.csfs[r], subshellList)
                opa  = SpinAngular.TwoParticleOperator(0, plus, true)
                waG2 = SpinAngular.computeCoefficients(opa, basis.csfs[r], basis.csfs[r], subshellList)
                wa   = [waG1, waG2]
            end
            if  Defaults.saRatip() && Defaults.saGG() && true
                if  length(waR[1]) != 0     println(  ">> Angular coeffients from Ratip2013   = $(waR[1]) ")    end
                if  length(waG1)   != 0     println("\n>> Angular coeffients from SpinAngular = $waG1 ")        end
                if  length(waR[2]) != 0     println(  ">> Angular coeffients from Ratip2013   = $(waR[2]) ")    end
                if  length(waG2)   != 0     println("\n>> Angular coeffients from SpinAngular = $waG2 ")        end
            end
            #
            me = 0.
            for  coeff in wa[1]
                jj = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                me = me + coeff.T * sqrt( jj + 1) * RadialIntegrals.GrantIab(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid, potential)
            end

            for  coeff in wa[2]
                if  settings.eeInteractionCI in [CoulombInteraction(), CoulombBreit()]
                    me = me + coeff.V * InteractionStrength.XL_Coulomb(coeff.nu, basis.orbitals[coeff.a], basis.orbitals[coeff.b],
                                                                                basis.orbitals[coeff.c], basis.orbitals[coeff.d], grid, keep=keep)   end
                if  settings.eeInteractionCI in [BreitInteraction(), CoulombBreit()]
                    me = me + coeff.V * 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: LinearAlgebra", 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)

        levelNos = Int64[]
        for (ilev, level) in  enumerate(mp.levels)       push!(levelNos, ilev)   end

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

        return( mp )
    end

    """
    `Basics.performCI(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.performCI(basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid, settings::AsfSettings; printout::Bool=true)

        # Determine the J^P symmetry blocks
        symmetries = Dict{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
            # Skip the symmetry block if it not selected
            ##x @show  sym, settings.levelSelectionCI,  Basics.selectSymmetry(sym, settings.levelSelectionCI)
            if  !Basics.selectSymmetry(sym, settings.levelSelectionCI)     continue    end
            matrix = compute("matrix: CI, J^P symmetry", sym, basis, nuclearModel, grid, settings; printout=printout)
            #
            ##x if  size(matrix,1) == 2   matrix[1,2] = matrix[2,1] = -1.0 * matrix[2,1]     end
            eigen  = Basics.diagonalize("matrix: LinearAlgebra", matrix)
            ##x eigen  = Basics.diagonalize("matrix: Julia, eigfact", matrix)
            ##x @show matrix, eigen
            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)

        # Determine the level list to be printed out
        levelNos = Int64[]
        for (ilev, level) in  enumerate(mp.levels)
            if  Basics.selectLevel(level, settings.levelSelectionCI)    push!(levelNos, ilev)    end
        end

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

        return( mp )
    end

    """
    `Basics.performCI(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.performCI(basis::Basis, nuclearModel::Nuclear.Model, grid::Radial.Grid,
                              settings::AsfSettings, plasmaSettings::PlasmaShift.Settings; printout::Bool=true)

        # Determine the J^P symmetry blocks
        symmetries = Dict{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: LinearAlgebra", 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)

        levelNos = Int64[]
        for (ilev, level) in  enumerate(mp.levels)       push!(levelNos, ilev)   end

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

        return( mp )
    end


text:
module Scale

using Colors
using Compose
using DataStructures
using Gadfly
using Showoff
using IndirectArrays
using CategoricalArrays
using Printf
using Base.Iterators, Dates

import Gadfly: element_aesthetics, isconcrete, concrete_length, discretize_make_ia,
    aes2str, valid_aesthetics, Maybe
import Distributions: Distribution

include("color_misc.jl")

iscategorical(scales::Dict{Symbol, Gadfly.ScaleElement}, var::Symbol) =
        haskey(scales, var) && isa(scales[var], DiscreteScale)

function apply_scales(scales, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for scale in scales
        apply_scale(scale, aess, datas...)
    end

    for (aes, data) in zip(aess, datas)
        aes.titles = data.titles
    end
end

function apply_scales(scales, datas::Gadfly.Data...)
    aess = Gadfly.Aesthetics[Gadfly.Aesthetics() for _ in datas]
    apply_scales(scales, aess, datas...)
    aess
end

# Transformations on continuous scales
struct ContinuousScaleTransform
    f::Function     # transform function
    finv::Function  # f's inverse

    # A function taking one or more values and returning an array of
    # strings.
    label::Function
end

identity_formatter(xs::AbstractArray, format=:auto) = showoff(xs, format)
const identity_transform = ContinuousScaleTransform(identity, identity, identity_formatter)

log10_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("10%s", x) for x in showoff(xs, format)]
const log10_transform = ContinuousScaleTransform(log10, x -> 10^x, log10_formatter)

log2_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("2%s", x) for x in showoff(xs, format)]
const log2_transform = ContinuousScaleTransform(log2, x -> 2^x, log2_formatter)

log_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("e%s", x) for x in showoff(xs, format)]
const log_transform = ContinuousScaleTransform(log, exp, log_formatter)

asinh_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("sinh(%s)", x) for x in showoff(xs, format)]
const asinh_transform = ContinuousScaleTransform(asinh, sinh, asinh_formatter)

sqrt_formatter(xs::AbstractArray, format=:plain) =
        [@sprintf("%s2", x) for x in showoff(xs, format)]
const sqrt_transform = ContinuousScaleTransform(sqrt, x -> x^2, sqrt_formatter)

# Continuous scale maps data on a continuous scale simple by calling
# `convert(Float64, ...)`.
struct ContinuousScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}
    trans::ContinuousScaleTransform
    minvalue
    maxvalue
    minticks
    maxticks
    labels::Union{(Nothing), Function}
    format
    scalable

    function ContinuousScale(vars, trans,
                             minvalue, maxvalue, minticks, maxticks,
                             labels, format, scalable)
        minvalue != nothing && maxvalue != nothing && minvalue > maxvalue &&
                error("Cannot construct a ContinuousScale with minvalue > maxvalue")
        new(vars, trans, minvalue, maxvalue, minticks, maxticks, labels, format, scalable)
    end
end

function ContinuousScale(vars, trans;
                         labels=nothing,
                         minvalue=nothing, maxvalue=nothing,
                         minticks=2, maxticks=10, format=nothing, scalable=true)
    ContinuousScale(vars, trans, minvalue, maxvalue, minticks, maxticks, labels, format, scalable)
end

function make_labeler(scale::ContinuousScale)
    if scale.labels != nothing
        xs -> [scale.labels(x) for x in xs]
    elseif scale.format == nothing
        scale.trans.label
    else
        xs -> scale.trans.label(xs, scale.format)
    end
end

const x_vars = [:x, :xmin, :xmax, :xintercept, :xviewmin, :xviewmax, :xend]
const y_vars = [:y, :ymin, :ymax, :yintercept, :intercept, :middle, :upper_fence, :lower_fence,
                :upper_hinge, :lower_hinge, :yviewmin, :yviewmax, :yend]

element_aesthetics(scale::ContinuousScale) = scale.vars

function continuous_scale_partial(vars::Vector{Symbol}, trans::ContinuousScaleTransform)
    function(; minvalue=nothing, maxvalue=nothing, labels=nothing, format=nothing,
             minticks=2, maxticks=10, scalable=true)
        ContinuousScale(vars, trans, minvalue=minvalue, maxvalue=maxvalue,
                        labels=labels, format=format, minticks=minticks,
                        maxticks=maxticks, scalable=scalable)
    end
end

xy_continuous_docstr(var, aess) = """
    $(var)_continuous[(; minvalue=nothing, maxvalue=nothing, labels=nothing,
                       format=nothing, minticks=2, maxticks=10, scalable=true)]

Map $(aess) to $var positions in Cartesian coordinates, which are presumed to
be numerical, using an identity transform.  `minvalue` and `maxvalue` set soft
lower and upper bounds.  (Use [`Coord.cartesian`](@ref) to enforce a hard
bound.)  `labels` is a function which maps a coordinate value to a string
label.  `format` is one of `:plain`, `:scientific`, `:engineering`, or `:auto`.
Set `scalable` to false to prevent zooming on this axis.  See also
[`$(var)_log10`](@ref), [`$(var)_log2`](@ref), [`$(var)_log`](@ref),
[`$(var)_asinh`](@ref), and [`$(var)_sqrt`](@ref) for alternatives to the
identity transform.
"""

# can be put on two lines with julia 0.7
@doc xy_continuous_docstr("x", aes2str(element_aesthetics(x_continuous()))) const x_continuous = continuous_scale_partial(x_vars, identity_transform)

@doc xy_continuous_docstr("y", aes2str(element_aesthetics(y_continuous()))) const y_continuous = continuous_scale_partial(y_vars, identity_transform)

xy_fun_docstr(var,fun) = """
    $(var)_$(fun)[(; minvalue=nothing, maxvalue=nothing, labels=nothing,
                       format=nothing, minticks=2, maxticks=10, scalable=true)]

Similar to [`Scale.$(var)_continuous`](@ref), except that the aesthetics are
`$(fun)` transformed and the `labels` function inputs transformed values.
"""

for xy in [:x,:y], fun in [:log10, :log2, :log, :asinh, :sqrt]
  @eval @doc $(xy_fun_docstr(xy,fun)) const $(Symbol(xy,'_',fun)) = continuous_scale_partial($(Symbol(xy,"_vars")), $(Symbol(fun,"_transform")))
end

"""
    size_continuous[(; minvalue=nothing, maxvalue=nothing, labels=nothing,
                     format=nothing, minticks=2, maxticks=10, scalable=true)]

To be updated in Gadfly 2.0. Now try out the new [`Scale.size_radius`](@ref) and [`Scale.size_area`](@ref).
"""
const size_continuous = continuous_scale_partial([:size], identity_transform)

"""
    slope_continuous[(; minvalue=nothing, maxvalue=nothing, labels=nothing,
                     format=nothing, minticks=2, maxticks=10, scalable=true)]
"""
const slope_continuous = continuous_scale_partial([:slope], identity_transform)

"""
    alpha_continuous[(; minvalue=0.0, maxvalue=1.0, labels=nothing,
                     format=nothing, minticks=2, maxticks=10, scalable=true)]

Rescale the data values between `minvalue` and `maxvalue` to opacity (alpha) values between 0 and 1.
"""
alpha_continuous(; minvalue=0.0, maxvalue=1.0, labels=nothing, format=nothing, minticks=2, maxticks=10, scalable=true) =
     ContinuousScale([:alpha], identity_transform, minvalue=minvalue, maxvalue=maxvalue,
               labels=labels, format=format, minticks=minticks, maxticks=maxticks, scalable=scalable)

function apply_scale(scale::ContinuousScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            vals = getfield(data, var)
            if vals isa CategoricalArray
                throw(ArgumentError("continuous scale for $var aesthetic when stored as a CategoricalArray. Consider using a discrete scale or convert data to an Array."))
            end
            vals === nothing && continue

            if eltype(vals) <: Measure
                setfield!(aes, var, vals)
                continue
            end

            # special case for function arrays bound to :y
            # pass the function values through and wait for the scale to
            # be reapplied by Stat.func
            if var == :y && eltype(vals) <: Function
                aes.y = vals
                continue
            end

            # special case for Distribution values bound to :x or :y. wait for
            # scale to be re-applied by Stat.qq
            if in(var, [:x, :y]) && eltype(vals) <: Distribution
                setfield!(aes, var, vals)
                continue
            end

            T = Any
            for d in vals
                if isconcrete(d)
                    T = typeof(scale.trans.f(d))
                    break
                end
            end

#            T <: Measure && (T = Measure)

            ds = any(ismissing, vals) ? Array{Union{Missing,T}}(undef,length(vals)) :
                    Array{T}(undef,length(vals))
            apply_scale_typed!(ds, vals, scale)

            if var == :alpha
                ds = (vals.-scale.minvalue)./(scale.maxvalue-scale.minvalue)
             end

            if var == :xviewmin || var == :xviewmax ||
               var == :yviewmin || var == :yviewmax
                setfield!(aes, var, ds[1])
            else
                setfield!(aes, var, ds)
            end

            if var in x_vars
                label_var = :x_label
            elseif var in y_vars
                label_var = :y_label
            else
                label_var = Symbol(var, "_label")
            end

            if in(label_var, valid_aesthetics)
                setfield!(aes, label_var, make_labeler(scale))
            end
        end

        if scale.minvalue != nothing
            if scale.vars === x_vars
                aes.xviewmin = scale.trans.f(scale.minvalue)
            elseif scale.vars === y_vars
                aes.yviewmin = scale.trans.f(scale.minvalue)
            end
        end

        if scale.maxvalue != nothing
            if scale.vars === x_vars
                aes.xviewmax = scale.trans.f(scale.maxvalue)
            elseif scale.vars === y_vars
                aes.yviewmax = scale.trans.f(scale.maxvalue)
            end
        end
    end
end

function apply_scale_typed!(ds, field, scale::ContinuousScale)
    for i in 1:length(field)
        d = field[i]
        ds[i] = isconcrete(d) ? scale.trans.f(d) : d
    end
end

function discretize(values, levels=nothing, order=nothing, preserve_order=true)
    if levels == nothing
        if preserve_order
            levels = OrderedSet()
            for value in values
                push!(levels, value)
            end
            da = discretize_make_ia(values, collect(eltype(values), levels))
        else
            da = discretize_make_ia(values)
        end
    else
        da = discretize_make_ia(values, levels)
    end

    if order != nothing
        return discretize_make_ia(da, da.values[order])
    else
        return da
    end
end

struct DiscreteScaleTransform
    f::Function
end

struct DiscreteScale <: Gadfly.ScaleElement
    vars::Vector{Symbol}
    labels::Union{(Nothing), Function}
    levels::Union{(Nothing), AbstractVector}
    order::Union{(Nothing), AbstractVector}
end
DiscreteScale(vals; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale(vals, labels, levels, order)

const discrete = DiscreteScale

element_aesthetics(scale::DiscreteScale) = scale.vars

xy_discrete_docstr(var, aess) = """
    $(var)_discrete[(; labels=nothing, levels=nothing, order=nothing)]

Map $(aess), which are presumed to be categorical, to Cartesian coordinates.
Unlike [`Scale.x_continuous`](@ref), each unique $var value will be mapped to
equally spaced positions, regardless of value.

By default continuous scales are applied to numerical data. If data consists of
numbers specifying categories, explicitly adding `Scale.$(var)_discrete` is the
easiest way to get that data to plot appropriately.

`labels` is either a function which maps a coordinate value to a string label,
or a vector of strings of the same length as the number of unique values in the
aesthetic.  `levels` gives values for the scale.  Order will be respected and
anything in the data that's not respresented in `levels` will be set to
`missing`.  `order` is a vector of integers giving a permutation of the levels
default order.

See also [`group_discrete`](@ref), [`shape_discrete`](@ref),
[`size_discrete`](@ref), [`linestyle_discrete`](@ref), and [`alpha_discrete`](@ref).
"""

@doc xy_discrete_docstr("x", aes2str(element_aesthetics(x_discrete()))) x_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale(x_vars, labels=labels, levels=levels, order=order)

@doc xy_discrete_docstr("y", aes2str(element_aesthetics(y_discrete()))) y_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale(y_vars, labels=labels, levels=levels, order=order)

type_discrete_docstr(aes) = """
    $(aes)_discrete[(; labels=nothing, levels=nothing, order=nothing)]

Similar to [`Scale.x_discrete`](@ref), except applied to the `$aes` aesthetic.
"""

@doc type_discrete_docstr("group") group_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:group], labels=labels, levels=levels, order=order)

@doc type_discrete_docstr("shape") shape_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:shape], labels=labels, levels=levels, order=order)

@doc type_discrete_docstr("size") size_discrete(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:size], labels=labels, levels=levels, order=order)

"""
    alpha_discrete[(; labels=nothing, levels=nothing, order=nothing)]

Similar to [`Scale.x_discrete`](@ref), except applied to the `alpha` aesthetic. The alpha palette
is set by `Theme(alphas=[])`.
"""
alpha_discrete(; labels=nothing, levels=nothing, order=nothing) =
            DiscreteScale([:alpha], labels=labels, levels=levels, order=order)

function apply_scale(scale::DiscreteScale, aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)

    d = []
    for data in datas
        for var in scale.vars
            datavar = getfield(data, var)
            (datavar===nothing || eltype(datavar)<:Measure) && continue
            append!(d, skipmissing(datavar))
        end
    end
    levelset = unique(d)

    scale_levels = (scale.levels==nothing) ? [levelset...] : scale.levels

    for (aes, data) in zip(aess, datas)
        for var in scale.vars
            label_var = Symbol(var, "_label")
            datavar = getfield(data, var)
            datavar===nothing && continue

            if eltype(datavar) <: Measure
                setfield!(aes, var, datavar)
                continue
            end

            disc_data = discretize(datavar, scale_levels, scale.order)
            setfield!(aes, var, discretize_make_ia(disc_data.index))

            # The leveler for discrete scales is a closure over the discretized data.
            if scale.labels === nothing
                function default_labeler(xs)
                    lvls = filter(!ismissing, disc_data.values)
                    vals = Any[1 <= x <= length(lvls) ? lvls[x] : "" for x in xs]
                    if all([isa(val, AbstractFloat) for val in vals])
                        return showoff(vals)
                    else
                        return [string(val) for val in vals]
                    end
                end
                labeler = default_labeler
            else
                function explicit_labeler(xs)
                    lvls = filter(!ismissing, disc_data.values)
                    return [string(scale.labels(lvls[x])) for x in xs]
                end
                labeler = explicit_labeler
            end

            in(label_var, valid_aesthetics) && setfield!(aes, label_var, labeler)
        end
    end
end

# Error message if a plot has both a discrete and continuous color variable
const colorclash = """There appears to be mixing of a continuous variable and a discrete variable in the `color` aesthetic.
This can happen when the color aesthetic is assigned in more than 1 layer (which is ok to do),
or you have forgotten to assign it. Either:
    1) make sure that the `color` variable in all layers is the same (i.e. discrete or continuous), or
    2) add `Scale.color_continuous` or `Scale.color_discrete` to the plot, or
    3) in a layer where `color` is a discrete variable, instead assign e.g. `color=[colorant"red"]`, or a vector of colorants."""

struct NoneColorScale <: Gadfly.ScaleElement
end

"""
    color_none

Suppress the default color scale that some statistics impose by setting the
`color` aesthetic to `nothing`. (To be deprecated)
"""

function color_none()
    @warn """`Scale.color_none` to be deprecated. Instead use e.g. `plot(..., Geom.contour, color=[colorant"black"])`"""
    return NoneColorScale()
end

element_aesthetics(scale::NoneColorScale) = [:color]

function apply_scale(scale::NoneColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for aes in aess
        aes.color = nothing
    end
end

struct IdentityColorScale <: Gadfly.ScaleElement
end

"""
    color_identity
"""
const color_identity = IdentityColorScale

element_aesthetics(scale::IdentityColorScale) = [:color]

function apply_scale(scale::IdentityColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        data.color === nothing && continue
        aes.color = discretize_make_ia(data.color)
        aes.color_key_colors = Dict()
    end
end

struct DiscreteColorScale <: Gadfly.ScaleElement
    f::Function
    levels::Union{(Nothing), AbstractVector}
    order::Union{(Nothing), AbstractVector}
    preserve_order::Bool
end
DiscreteColorScale(f; levels=nothing, order=nothing, preserve_order=true) =
        DiscreteColorScale(f, levels, order, preserve_order)

element_aesthetics(scale::DiscreteColorScale) = [:color]

default_discrete_colors(n) = convert(Vector{Color},
        distinguishable_colors(n,
                               [LCHab(70, 60, 240)],
                               transform=c -> deuteranopic(c, 0.5),
                               lchoices=Float64[65, 70, 75, 80],
                               cchoices=Float64[0, 50, 60, 70],
                               hchoices=range(0, stop=330, length=24)))

"""
    color_discrete_hue[(f; levels=nothing, order=nothing, preserve_order=true)]

Create a discrete color scale that maps the categorical levels in the `color`
aesthetic to `Color`s.  `f` is a function that produces a vector of colors.
`levels` gives values for the scale.  Order will be respected and anything in
the data that's not represented in `levels` will be set to missing.  `order` is
a vector of integers giving a permutation of the levels default order.  If
`preserve_order` is `true` orders levels as they appear in the data.

Either input `Stat.color_discrete_hue` as an argument to `plot`, or set
`discrete_color_scale` in a [Theme](@ref Themes).

# Examples

```
julia> x = Scale.color_discrete_hue()
Gadfly.Scale.DiscreteColorScale(Gadfly.Scale.default_discrete_colors, nothing, nothing, true)

julia> x.f(3)
3-element Array{ColorTypes.Color,1}:
 LCHab{Float32}(70.0,60.0,240.0)
 LCHab{Float32}(80.0,70.0,100.435)
 LCHab{Float32}(65.8994,62.2146,353.998)
```
"""
color_discrete_hue(f=default_discrete_colors;
                   levels=nothing, order=nothing, preserve_order=true) =
        DiscreteColorScale(f, levels=levels, order=order, preserve_order=preserve_order)

@deprecate discrete_color_hue(; levels=nothing, order=nothing, preserve_order=true) color_discrete_hue(; levels=levels, order=order, preserve_order=preserve_order)

const color_discrete = color_discrete_hue   ### WHY HAVE THIS ALIAS?

@deprecate discrete_color(; levels=nothing, order=nothing, preserve_order=true) color_discrete(; levels=levels, order=order, preserve_order=preserve_order)

"""
    color_discrete_manual(colors...; levels=nothing, order=nothing)

Similar to [`color_discrete_hue`](@ref) except that colors for each level are
specified directly instead of being computed by a function.
"""
color_discrete_manual(colors::AbstractString...; levels=nothing, order=nothing) =
        color_discrete_manual(Gadfly.parse_colorant(colors)...; levels=levels, order=order)

function color_discrete_manual(colors::Color...; levels=nothing, order=nothing)
    cs = [colors...]
    f = n -> distinguishable_colors(n, cs)
    DiscreteColorScale(f, levels=levels, order=order)
end

@deprecate discrete_color_manual(colors...; levels=nothing, order=nothing) color_discrete_manual(colors...; levels=levels, order=order)

function apply_scale(scale::DiscreteColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    levelset = OrderedSet()
    for (aes, data) in zip(aess, datas)
        (data.color === nothing || eltype(data.color)<:Colorant) && continue
        for d in data.color
            # Remove missing values
            # FixMe! The handling of missing values shouldn't be this scattered across the source
            ismissing(d) || push!(levelset, d)
        end
    end

    if scale.levels == nothing
        scale_levels = [levelset...]
        scale.preserve_order || sort!(scale_levels)
    else
        scale_levels = scale.levels
    end
    scale.order == nothing || permute!(scale_levels, scale.order)
    colors = convert(Vector{RGB{Float32}}, scale.f(length(scale_levels)))

    color_map = Dict([(color, string(label))
                              for (label, color) in zip(scale_levels, colors)])
    labeler(xs) = [color_map[x] for x in xs]

    for (aes, data) in zip(aess, datas)
        data.color === nothing && continue
        if eltype(data.color)<:Colorant
            aes.color = data.color
            continue
        end
        # Remove missing values
        # FixMe! The handling of missing values shouldn't be this scattered across the source
        ds = discretize([c for c in data.color if !ismissing(c)], scale_levels)

        colorvals = colors[ds.index]

        colored_ds = discretize_make_ia(colorvals, colors)
        aes.color = colored_ds

        aes.color_label = labeler
        aes.color_key_colors = OrderedDict()
        for (i, c) in enumerate(colors)
            aes.color_key_colors[c] = i
        end
    end
end

struct ContinuousColorScale <: Gadfly.ScaleElement
    # A function of the form f(p) where 0 <= p <= 1, that returns a color.
    f::Function
    trans::ContinuousScaleTransform

    minvalue
    maxvalue
end
ContinuousColorScale(f, trans=identity_transform; minvalue=nothing, maxvalue=nothing) =
        ContinuousColorScale(f, trans, minvalue, maxvalue)

function continuous_color_scale_partial(trans::ContinuousScaleTransform)
    lch_diverge2 = function(l0=30, l1=100, c=40, h0=260, h1=10, hmid=20, power=1.5)
        lspan = l1 - l0
        hspan1 = hmid - h0
        hspan0 = h1 - hmid
        function(r)
            r2 = 2r - 1
            return LCHab(min(80, l1 - lspan * abs(r2)^power), max(10, c * abs(r2)),
                         (1-r)*h0 + r * h1)
        end
    end

    function(; minvalue=nothing, maxvalue=nothing, colormap=lch_diverge2())
        ContinuousColorScale(colormap, trans, minvalue=minvalue, maxvalue=maxvalue)
    end
end

element_aesthetics(::ContinuousColorScale) = [:color]

"""
    color_continuous[(; minvalue=nothing, maxvalue=nothing, colormap)]

Create a continuous color scale by mapping
$(aes2str(element_aesthetics(color_continuous()))) to a `Color`.  `minvalue`
and `maxvalue` specify the data values corresponding to the bottom and top of
the color scale.  `colormap` is a function defined on the interval from 0 to 1
that returns a `Color`. See also [`lab_gradient`](@ref).

Either input `Stat.color_continuous` as an argument to `plot`, or
set `continuous_color_scale` in a [Theme](@ref Themes).

See also [`color_log10`](@ref), [`color_log2`](@ref), [`color_log`](@ref),
[`color_asinh`](@ref), and [`color_sqrt`](@ref).
"""
const color_continuous = continuous_color_scale_partial(identity_transform)

color_fun(fun) = """
    color_$(fun)[(; minvalue=nothing, maxvalue=nothing, colormap)]

Similar to [`Scale.color_continuous`](@ref), except that `color` is $(fun) transformed.
"""

for fun in [:log10, :log2, :log, :asinh, :sqrt]
  @eval @doc $(color_fun(fun)) const $(Symbol("color_",fun)) = continuous_color_scale_partial($(Symbol(fun,"_transform")))
end

const color_continuous_gradient = color_continuous  ### WHY HAVE THIS ALIAS?

@deprecate continuous_color_gradient(;minvalue=nothing, maxvalue=nothing) color_continuous_gradient(;minvalue=minvalue, maxvalue=maxvalue)

@deprecate continuous_color(;minvalue=nothing, maxvalue=nothing) color_continuous(;minvalue=nothing, maxvalue=nothing)

function apply_scale(scale::ContinuousColorScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)

    for (aes, data) in zip(aess, datas)
        eltype(data.color)<:Colorant && (aes.color = data.color)
    end

    cdata = skipmissing(
        Iterators.flatten(i.color for i in datas if !(i.color == nothing || eltype(i.color)<:Colorant)))
    if !isempty(cdata)
        typeof(first(cdata))<:Union{Number, Date} || error(colorclash)
        cmin, cmax = extrema(cdata)
    else
        return
    end

    strict_span = false
    scale.minvalue != nothing && scale.maxvalue != nothing && (strict_span=true)
    scale.minvalue != nothing && (cmin=scale.minvalue)
    scale.maxvalue != nothing && (cmax=scale.maxvalue)

    cmin, cmax = promote(cmin, cmax)

    cmin = scale.trans.f(cmin)
    cmax = scale.trans.f(cmax)

    ticks, viewmin, viewmax = Gadfly.optimize_ticks(cmin, cmax, strict_span=strict_span)
    ticks[1] == 0 && cmin >= 1 && !strict_span && (ticks[1] = 1)

    cmin, cmax = ticks[1], ticks[end]
    cspan = cmax != cmin ? cmax - cmin : 1.0

    for (aes, data) in zip(aess, datas)
        (data.color===nothing  || eltype(data.color)<:Colorant) && continue

        aes.color = Array{RGB{Float32}}(undef, length(data.color))
        apply_scale_typed!(aes.color, data.color, scale, cmin, cspan)

        color_key_colors = Dict{Color, Float64}()
        color_key_labels = Dict{Color, AbstractString}()

        tick_labels = scale.trans.label(ticks)
        for (i, j, label) in zip(ticks, ticks[2:end], tick_labels[1:end-1])
            r = (i - cmin) / cspan
            c = scale.f(r)
            color_key_colors[c] = r
            color_key_labels[c] = label
        end
        c = scale.f((ticks[end] - cmin) / cspan)
        color_key_colors[c] = (ticks[end] - cmin) / cspan
        color_key_labels[c] = tick_labels[end]

        labeler(xs) = [get(color_key_labels, x, "") for x in xs]

        aes.color_function = scale.f
        aes.color_label = labeler
        aes.color_key_colors = color_key_colors
        aes.color_key_continuous = true
    end
end

function apply_scale_typed!(ds, field, scale::ContinuousColorScale,
                            cmin, cspan)
    for (i, d) in enumerate(field)
        if isconcrete(d)
            ds[i] = convert(RGB{Float32},
                        scale.f((scale.trans.f(d) - cmin) / cspan))
        else
            ds[i] = missing
        end
    end
end

# Label scale is always discrete, hence we call it 'label' rather
# 'label_discrete'.
struct LabelScale <: Gadfly.ScaleElement
end

function apply_scale(scale::LabelScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        data.label === nothing && continue
        aes.label = discretize(data.label)
    end
end

element_aesthetics(::LabelScale) = [:label]

"""
   label
"""
const label = LabelScale

"""
    xgroup[(; labels=nothing, levels=nothing, order=nothing)]

A discrete scale for use with [`Geom.subplot_grid`](@ref Gadfly.Geom.subplot_grid).
"""
xgroup(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:xgroup], labels=labels, levels=levels, order=order)

"""
    ygroup[(; labels=nothing, levels=nothing, order=nothing)]

A discrete scale for use with [`Geom.subplot_grid`](@ref Gadfly.Geom.subplot_grid).
"""
ygroup(; labels=nothing, levels=nothing, order=nothing) =
        DiscreteScale([:ygroup], labels=labels, levels=levels, order=order)

# Catchall scale for when no transformation of the data is necessary
struct IdentityScale <: Gadfly.ScaleElement
    var::Symbol
end

element_aesthetics(scale::IdentityScale) = [scale.var]

function apply_scale(scale::IdentityScale,
                     aess::Vector{Gadfly.Aesthetics}, datas::Gadfly.Data...)
    for (aes, data) in zip(aess, datas)
        getfield(data, scale.var) === nothing && continue
        setfield!(aes, scale.var, getfield(data, scale.var))
    end
end

### could these be defined as `const z_func = ` ?
"""
    x_func()
"""
x_func() = IdentityScale(:x)

"""
    z_func()
"""
z_func() = IdentityScale(:z)

"""
    y_func()
"""
y_func() = IdentityScale(:y)

"""
    x_distribution()
"""
x_distribution() = IdentityScale(:x)

"""
    y_distribution()
"""
y_distribution() = IdentityScale(:y)

"""
    shape_identity()
"""
shape_identity() = IdentityScale(:shape)

"""
    size_identity()
"""
size_identity() = IdentityScale(:size)

linestyle_identity() = IdentityScale(:linestyle)

include("scale/scales.jl")

end # module Scale


text:
# module ExponentialFamilies
using Distributions

# A parameterized exponential family
struct ExponentialFamily{P,X}
    logh :: Function # :: X -> Real
    logg :: Function # :: P -> Real
    η    :: Function # :: P -> StaticVector{N, Real}
    t    :: Function # :: X -> StaticVector{N, Real}
end

# A specific instance of a distribution belonging to an exponential family
struct ExpFamDist{P,X}
    fam :: ExponentialFamily{P,X}
    θ :: P
end

# Find the exponential family distribution instance for a given distribution.
# This needs to be reworked, since many distributions are not from exponential families
efam(d::ExpFamDist{P,X} where {P,X}) = d

Distributions.logdensity(d::ExpFamDist{P,X}, x :: X) where {P,X} = d.fam.logh(x) + d.fam.logg(d.θ) + sum(d.fam.η(d.θ) .* d.fam.t(x))

# Exponential families are closed under iid
function iid(n::Integer, d::ExpFamDist{P,X} ) where {P,X}
    fam = d.fam
    logh(xs) = sum(fam.logh.(xs))
    logg(θ) = n * fam.logg(θ)
    η(θ) = fam.η(θ)
    t(xs) = sum(fam.t.(xs))

    newfam = ExponentialFamily{P,Vector{X}}(logh, logg, η, t)
    ExpFamDist{P,Vector{X}}(newfam,d.θ)
end

iid(n, d ) = iid(n,efam(d))

################################
# Bernoulli example

logit(p) = log(p/(1-p))

function efam(d::Bernoulli{P}) where {P <: Real}
    logh(x) = 0
    logg(p::P) = -log(1 - p)
    η(p::P) = logit(p)
    t(x) = Float64(x)

    fam = ExponentialFamily{P,Bool}(logh, logg, η, t)

    ExpFamDist(fam, d.p)
end

efam(Bernoulli(0.2))

binom = iid(1000,Bernoulli(0.2))

using BenchmarkTools

@btime logdensity(binom, repeat([true],1000))

@btime sum(logdensity.(Bernoulli(0.2),repeat([true],1000)))

################################
# Normal example

const sqrt2π = sqrt(2*π)

const invsqrt2π = 1/sqrt(2*π)

function efam(d::Normal{P}) where {P <: Real}
    logh(x) = invsqrt2π
    logg(θ) = 0.5 * (θ.μ/θ.σ)^2 + log(θ.σ)
    η(θ) = (θ.μ/θ.σ^2, -1/(2*θ.σ^2))
    t(x) = (x,x^2)

    fam = ExponentialFamily{NamedTuple{(:μ, :σ),Tuple{P,P}}, Float64}(logh, logg, η, t)

    ExpFamDist(fam, (μ=d.μ,σ=d.σ))
end

const xs = 3 .+ 4 .* randn(1000);

@btime logdensity(iid(1000,Normal(3.1,4.3)), xs)

@btime sum(logdensity.(Normal(3,4),xs))

# end # module


text:
using DrWatson
quickactivate(@__DIR__)
@show projectdir()
##
using CovidSurvey
using CSV
using DataFrames
using Dates
# using RCall
using Turing
using StatsPlots
using Plots.PlotMeasures
using Random
## ============================================================================

## ===========================================================================
# load data and model
# ============================================================================
fdir_p = projectdir("reports/2021-06-07/")
fdir_np = projectdir("reports/variable-delays/")
# fname = "cases=2020-06-01_chains=10_hospit=true_model=parametric-cases_observ=2021-01-01_preds=CF,CC,CR,CS_sero=true_steps=4000_warmup=1000"
# fname = "cases=2020-06-01_chains=10_hospit=true_model=parametric-hospit_observ=2021-01-01_preds=CF,CC,CR,CS_sero=true_steps=4000_warmup=1000"
fname_p = "cases=2020-06-01_chains=10_hospit=true_model=parametric-deaths_observ=2021-01-01_preds=CF,CC,CR,CS_sero=true_steps=4000_warmup=1000"
fname_np = "cases=2020-06-01_chains=10_hospit=true_model=deaths_observ=2021-02-06_sero=true_steps=4000_warmup=1000"

## read chain
chain = read(fdir*"CHAIN_"*fname*".jls", Chains)
chain = chain[1:10:end,:,:]
# chain = chain[1:10:end,:,[1,2,3,4,6,7,8,9,10]]
diagnostics = gelmandiag(chain)
"/home/and/dev/CovidSurvey/reports/2021-06-07/Rt_"
Plots.scalefontsizes(1.3)
pgfplotsx()
p = let
    rtp = load(fdir*"Rt_"*fname*".csv") |> DataFrame
    p = plot(legend=:bottomright, ylabel="Rt")
    plot!(rtp.date, rtp.mq, ribbon=(rtp.mq - rtp.llq, rtp.uuq - rtp.mq), lw=2, label="with predictors (contact data)")
    plot!(rt.date, rt.mq, ribbon=(rt.mq - rt.llq, rt.uuq - rt.mq), lw=2, label="without predictors")
    # xlims!((t[1],t[end]))
end
savefig(projectdir("figures/compare_Rt.png"))

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey