text:
    export generate

    """
    `Basics.generate(representation::AtomicState.Representation)`
        ... to generate an atomic representation as specified by the representation.repType::AbstractRepresentationType.
            All relevant intermediate and final results are printed to screen (stdout). Nothing is returned.

    `Basics.perform(representation::AtomicState.Representation; output=true)`
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and
            specifications of the representation but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(representation::AtomicState.Representation; output::Bool=false)
        results = Basics.generate(representation.repType, representation; output=output)

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

    """
    `Basics.generate(repType::AtomicState.MeanFieldBasis, representation::AtomicState.Representation)`
        ... to generate a mean-field basis (representation) for a set of reference configurations; all relevant intermediate
            and final results are printed to screen (stdout). Nothing is returned.

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

        # The asfSettings only define the SCF part and are partly derived from the MeanFieldSettings
        asfSettings = AsfSettings(AsfSettings(); scField=repType.settings.scField)

        basis      = Basics.performSCF(rep.refConfigs, nModel, rep.grid, asfSettings; printout=true)
        if output    results = Base.merge( results, Dict("mean-field basis" => basis) )          end

        return( results )
    end

    """
    `Basics.generate(repType::AtomicState.CiExpansion, representation::AtomicState.Representation)`
        ... to generate a configuration-interaction expansion for a single level symmetry and based on a set of reference configurations
            and a number of pre-specified steps. All relevant intermediate and final results are printed to screen (stdout).
            Nothing is returned.

    `Basics.generate(repType::AtomicState.CiExpansion, representation::AtomicState.Representation; output=true)`
        ... to generate 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.generate(repType::AtomicState.CiExpansion, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel    = rep.nuclearModel
        orbitals  = repType.applyOrbitals

        # Generate a list of relativistic configurations and  CSF's for the given subshell list
        relconfList = ConfigurationR[]
        for  conf in rep.refConfigs
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        ##x if  false    for  i = 1:length(relconfList)    println(">> include ", relconfList[i])    end   end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        ##x Defaults.setDefaults("relativistic subshell list", subshellList; printout=printout)
        csfList = CsfR[]
        for  relconf in relconfList
            newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList)
            append!( csfList, newCsfs)
        end

        # Generate all level symmetries that need to be taken into account
        if     length(repType.settings.levelSelectionCI.symmetries) != 0    symmetries = repType.settings.levelSelectionCI.symmetries
        else   symmetries = LevelSymmetry[]
            for csf in csfList
                if      LevelSymmetry(csf.J, csf.parity)  in  symmetries
                else    push!( symmetries, LevelSymmetry(csf.J, csf.parity) )
                end
            end
        end
        println("*** Level symmetries = $symmetries ")

        # The asfSettings only define the CI part and are partly derived from the CiSettings
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),    0, 0., Subshell[], Subshell[],
                                  repType.settings.eeInteractionCI, NoneQed(), LSjjSettings(false),
                                  repType.settings.levelSelectionCI)

        basis      = Basics.generateBasis(rep.refConfigs, symmetries, repType.excitations)
        basis      = Basis( true, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals )
        # Test that all orbitals are available
        for  subshell in  basis.subshells
            if   !(haskey(orbitals, subshell))  error("No orbital available for subshell $subshell ")     end
        end
        multiplet  = Basics.performCI(basis,  nModel, rep.grid, asfSettings; printout=true)
        if output    results = Base.merge( results, Dict("CI multiplet" => Multiplet("CI multiplet:", multiplet.levels)) )              end

        return( results )
    end

    """
    `Basics.generate(repType::AtomicState.RasExpansion, representation::AtomicState.Representation)`
        ... to generate a restricted active-space expansion for a single level symmetry and based on a set of reference configurations
            and a number of pre-specified steps. All relevant intermediate and final results are printed to screen (stdout).
            Nothing is returned.

    `Basics.generate(repType::AtomicState.RasExpansion, representation::AtomicState.Representation; output=true)`
        ... to generate 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.generate(repType::AtomicState.RasExpansion, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel   = rep.nuclearModel
        # First perform a SCF+CI computations for the reference configurations below to generate a spectrum of start orbitals
        asfSettings    = AsfSettings()  ## Use default settings to define a first multiplet from the reference configurations;
                                        ## all further details are specified for each step
        priorBasis     = Basics.performSCF(rep.refConfigs, nModel, rep.grid, asfSettings; printout=true)
        priorMultiplet = Basics.performCI(priorBasis, nModel, rep.grid, asfSettings; printout=true)
        nuclearPot     = Nuclear.nuclearPotential(nModel, rep.grid)
        ## electronicPot  = Basics.compute("radial potential: Dirac-Fock-Slater", rep.grid, priorMultiplet.levels[1].basis)
        ## meanPot        = Basics.add(nuclearPot, electronicPot)
        subshellList   = Basics.extractRelativisticSubshellList(rep)             ## extract all subshells that occur in the RAS computation
        startOrbitals  = Basics.generateOrbitalsForPotential(rep.grid, nuclearPot, subshellList)  ## generate a spectrum of sufficient size
        if output    results = Base.merge( results, Dict("reference multiplet" => Multiplet("Reference multiplet:", priorMultiplet.levels) ) )  end

        # The asfSettings only define the CI part of the RAS steps and partly derived from the RasSettings
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),    0, 0., Subshell[], Subshell[],
                                  repType.settings.eeInteractionCI, NoneQed(), LSjjSettings(true), repType.settings.levelSelectionCI )

        # Now, cycle over all steps of the RasExpansion
        for (istep, step)  in  enumerate(repType.steps)
            println("")
            printstyled("++ Compute the orbitals, orbitals and multiplet for step $istep ... \n", color=:light_green)
            printstyled("--------------------------------------------------------------      \n", color=:light_green)
            basis      = Basics.generateBasis(rep.refConfigs, [repType.symmetry], step)
            orbitals   = Basics.generateOrbitalsForBasis(basis, step.frozenShells, priorMultiplet.levels[1].basis, startOrbitals)
            basis      = Basis( true, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals )
            basis      = Basics.performSCF(basis, nModel, rep.grid, step.frozenShells, repType.settings; printout=true)
            multiplet  = Basics.performCI(basis,  nModel, rep.grid, asfSettings; printout=true)
            if output    results = Base.merge( results, Dict("step"*string(istep) => Multiplet("Multiplet:", multiplet.levels)) )              end
            priorMultiplet = multiplet
        end

        return( results )
    end

    """
    `Basics.generate(repType::AtomicState.GreenExpansion, representation::AtomicState.Representation)`
        ... to generate a Green (function) expansion for a given approach and excitation scheme of the electron,
            based on a set of reference configurations, a list of level symmetries as well as for given settings.
            All relevant intermediate and final results are printed to screen (stdout). Nothing is returned.

    `Basics.generate(repType::AtomicState.GreenExpansion, representation::AtomicState.Representation; output=true)`
        ... to generate the same but to return the complete output in a dictionary; the particular output depends on the type and
            specifications of the representation but can easily accessed by the keys of this dictionary.
    """
    function Basics.generate(repType::AtomicState.GreenExpansion, rep::AtomicState.Representation; output::Bool=false)
        if  output    results = Dict{String, Any}()    else    results = nothing    end
        nModel    = rep.nuclearModel
        settings  = repType.settings
        # First perform a SCF+CI computations for the reference configurations below to generate a spectrum of start orbitals
        asfSettings   = AsfSettings()  ## Use default settings to define a first multiplet from the reference configurations;
                                       ## all further details are specified for each step
        refBasis      = Basics.performSCF(rep.refConfigs, nModel, rep.grid, asfSettings; printout=true)
        refMultiplet  = Basics.performCI(refBasis, nModel, rep.grid, asfSettings; printout=true)
        nuclearPot    = Nuclear.nuclearPotential(nModel, rep.grid)
        electronicPot = Basics.compute("radial potential: Dirac-Fock-Slater", rep.grid, refBasis)
        meanPot       = Basics.add(nuclearPot, electronicPot)

        println("")
        printstyled("Compute an approximate Green function expansion ... \n", color=:light_green)
        printstyled("--------------------------------------------------- \n", color=:light_green)

        channels = AtomicState.GreenChannel[]

        # Generate all (non-relativistic) configurations from the bound configurations due to the given excitation scheme
        confList = Basics.generateConfigurationsForExcitationScheme(rep.refConfigs, repType.excitationScheme, settings.nMax, settings.lValues)
        # Print (if required) information about the generated configuration list
        if  settings.printBefore    Basics.display(stdout, confList)    end

        # Generate shell list abd a full single-electron spectrum for this potential
        subshellList = Basics.extractRelativisticSubshellList(confList)                      ## extract all subshells that occur in confList
        orbitals     = Basics.generateOrbitalsForPotential(rep.grid, meanPot, subshellList)  ## generate a spectrum of sufficient size
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)
        Basics.display(stdout, orbitals, rep.grid)

        # The asfSettings only define the CI part of the Green channels and are partly derived from the GreenSettings
        asfSettings = AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(),    0, 0., Subshell[], Subshell[],
                                  CoulombInteraction(), NoneQed(), LSjjSettings(false), settings.levelSelection )

        # Cycle over all selected level symmetries to generate the requested channels
        for  levelSymmetry  in  repType.levelSymmetries
            basis      = Basics.generateBasis(confList, [levelSymmetry])
            basis      = Basis( true, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals )
            multiplet  = Basics.computeMultipletForGreenApproach(repType.approach, basis, nModel, rep.grid, asfSettings, settings; printout=true)
            push!( channels, GreenChannel( levelSymmetry, multiplet) )
        end

        # Print all results to screen
        Basics.display(stdout, channels)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    Basics.display(iostream, channels)         end

        if output    results = Base.merge( results, Dict("Green channels" => channels) )   end
        return( results )
    end

    """
    `Basics.generate("condensed multiplet: by single weight", multiplet::Multiplet)`
        ... to condense/reduce the number of CSF in the basis of the given multiplet due to a single 'weight';
            a multiplet::Multiplet is returned.  **Not yet implemented !**
    """
    function Basics.generate(sa::String, multiplet::Multiplet)
        !(sa == "condensed multiplet: by single weight")   &&   error("Unsupported keystring = $sa")
        error("Not yet implemented !")
    end

    """
    `Basics.generate("configuration list: NR, from basis", basis::Basis)`
        ... to (re-) generate the list of NR configurations from the given basis; a confList::Array{Configuration,1} is returned.
    """
    function Basics.generate(sa::String, basis::Basis)
        !(sa == "configuration list: NR, from basis")   &&   error("Unsupported keystring = $sa")
        confList    = Configuration[]
        NoElectrons = sum( basis.csfs[1].occupation )
        for  csf in basis.csfs
            shellList = Dict{Shell,Int64}[]
            for  i = 1:length(csf.occupation)
                n = csf.subshell[i].n;    l = Basics.subshell_l( csf.subshell[i] );    sh = Shell(n, l)
                if     haskey(shellList, sh)    shellList(sh) = shellList(sh) + csf.occupation[i]
                else   shellList = merge( shellList, Dict(sh => csf.occupation[i]) )
                end
            end
            confNew = Configuration( shellList, NoElectrons )
            push!( confList, confNew )
        end
        confListNew = ConfigurationsExcludeDoubles(confList)

        return( confListNew )
    end

    """
    `Basics.generate("configuration list: NR, single-configuration", refConf::Configuration, NoExcitations::Int64, fromShells::Array{Shell,1},`
                     toShells::Array{Shell,1})
        ... to generate a non-relativistic configuration list, including the given reference configuration (refConf) and with
            all configurations that differ by NoExcitations from the fromShells into the toShells; an Array{Configuration,1}
            is returned.
    """
    function Basics.generate(sa::String, refConf::Configuration, NoExcitations::Int64, fromShells::Array{Shell,1}, toShells::Array{Shell,1})
        !(sa == "configuration list: NR, single-configuration")   &&   error("Unsupported keystring = $sa")
        confList = [refConf]
        # First prepare a proper reference configuration that also includes all fromShells and toShells with zero occupation
        shellDict = deepcopy(refConf.shells)
        allShells = deepcopy(fromShells)
        for  sh in toShells
            if     sh in allShells
            else   append!(allShells, [sh] )
            end
        end

        for  sh in allShells
        if  !haskey(shellDict, sh )    shellDict = merge( shellDict, Dict( sh => 0 ) )    end
        end

        rconf = Configuration(shellDict, refConf.NoElectrons)
        println("aa", rconf)

        # Add single excitations
        if      NoExcitations == 0
        elseif  NoExcitations == 1
            for  fromsh in fromShells
                for  tosh in toShells
                    confnew = deepcopy(rconf)
                    confnew.shells[fromsh] = confnew.shells[fromsh] - 1
                    confnew.shells[tosh]   = confnew.shells[tosh]   + 1
                    println("bb", confnew)

                    # Check that this is valid configuration and append, if appropriate
                    add = true
                    for (k,v) in  confnew.shells
                        if  v <= -1  ||  Parity(rconf) != Parity(confnew)    add = false    end
                    end
                    if  add
                        for (k,v) in confnew.shells
                        if  v == 0   delete!(confnew.shells, k)   end
                        end
                        push!(confList, confnew);   println("fromsh = $fromsh  tosh = $tosh")
                    end
                end
            end

        # Add double excitations
        elseif  NoExcitations == 2
            for  fromsha in fromShells
            for  fromshb in fromShells
                for  tosha in toShells
                for  toshb in toShells
                    confnew = deepcopy(rconf)
                    confnew.shells[fromsha] = confnew.shells[fromsha] - 1
                    confnew.shells[fromshb] = confnew.shells[fromshb] - 1
                    confnew.shells[tosha]   = confnew.shells[tosha]   + 1
                    confnew.shells[toshb]   = confnew.shells[toshb]   + 1
                    println("bb", confnew)

                    # Check that this is valid configuration and append, if appropriate
                    add = true
                    for (k,v) in  confnew.shells
                        if  v <= -1  ||  Parity(rconf) != Parity(confnew)    add = false    end
                    end
                    if  add
                        for (k,v) in confnew.shells
                        if  v == 0   delete!(confnew.shells, k)   end
                        end
                        push!(confList, confnew)
                    end
                end
                end
            end
            end

        else
            error("stop a")
        end

        # Exclude all configurations that appear twice (and more) in the list
        confList = Basics.ConfigurationsExcludeDoubles(confList)

        return( confList )
    end

    """
    `Basics.generate("configuration list: relativistic", conf::Configuration)`
        ... to split/decompose a non-relativistic configuration into an list of relativistic ConfigurationR[]. The proper
            occupuation of the relativistic subshells is taken into account.
    """
    function Basics.generate(sa::String, conf::Configuration)
        subshellList = Subshell[]
        confList     = ConfigurationR[]

        !(sa == "configuration list: relativistic")  &&   error("Unsupported keystring = $sa")

        initialize = true;    NoElectrons = 0
        for (k, occ)  in conf.shells
            NoElectrons     = NoElectrons + occ
            wa              = Basics.shellSplitOccupation(k, occ)
            subshellListNew = Dict{Subshell,Int64}[]

            if  initialize
                subshellListNew = wa;    initialize = false
            else
                for  s in 1:length(subshellList)
                    for  a in 1:length(wa)
                        wb = Base.merge( subshellList[s], wa[a] )
                        push!(subshellListNew, wb)
                    end
                end
            end
            subshellList = deepcopy(subshellListNew)
        end

        for subsh in subshellList
        wa = ConfigurationR(subsh, NoElectrons)
        push!(confList, wa)
        end

        return( confList )
    end

    """
    `Basics.generate("CSF list: from single ConfigurationR", conf::ConfigurationR, subshellList::Array{Subshell,1})`
        ... to construct from a given (relativistic) configuration all possible CSF with regard to the subshell order as specified
            by subshellList; a list::Array{CsfR,1} is returned.
    """
    function Basics.generate(sa::String, conf::ConfigurationR, subshellList::Array{Subshell,1})
        parity  = Basics.determineParity(conf)
        csfList = CsfR[];   useStandardSubshells = true;    first = true;    previousCsfs = CsfR[]
        # subhshellList = Subshell[];
        for  subsh in subshellList
            if   subsh in keys(conf.subshells)    occ = conf.subshells[subsh]    else    occ = 0    end
            if   first
                stateList   = ManyElectron.provideSubshellStates(subsh, occ)
                currentCsfs = CsfR[]
                for  state in stateList
                    push!( currentCsfs, CsfR( true, AngularJ64(state.Jsub2//2), parity, [state.occ], [state.nu],
                                                [AngularJ64(state.Jsub2//2)], [AngularJ64(state.Jsub2//2)], Subshell[]) )
                end
                previousCsfs = copy(currentCsfs)
                first        = false
            else
                # Now support also all couplings of the subshell states with the CSFs that were built-up so far
                stateList   = ManyElectron.provideSubshellStates(subsh, occ)
                currentCsfs = CsfR[]
                for  csf in  previousCsfs
                    for  state in stateList
                        occupation = deepcopy(csf.occupation);    seniorityNr = deepcopy(csf.seniorityNr);
                        subshellJ  = deepcopy(csf.subshellJ);     subshells = deepcopy(csf.subshells)
                        push!(occupation, state.occ);   push!(seniorityNr, state.nu);   push!(subshellJ, AngularJ64(state.Jsub2//2) )
                        push!(subshells, subsh)
                        newXList = oplus( csf.subshellX[end], AngularJ64(state.Jsub2//2) )
                        for  newX in newXList
                            subshellX = deepcopy(csf.subshellX);   push!(subshellX, newX)
                            push!( currentCsfs, CsfR( true, subshellX[end], parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[]) )
                        end
                    end
                end
                previousCsfs = copy(currentCsfs)
            end
        end

        return( previousCsfs )
    end

    """
    `Basics.generate("shells: ordered list for NR configurations", confs::Array{Configuration,1})`
        ... to generate for confs, i.e. all the given (non-relativistic) configurations, a common and ordered shell list;
            a list::Array{Shell,1} is returned.
    """
    function Basics.generate(sa::String, confs::Array{Configuration,1})
        shells = Shell[]

        !(sa == "shells: ordered list for NR configurations")  &&   error("Unsupported keystring = $sa")

        wa = Defaults.getDefaults("ordered shell list: non-relativistic", 7)
        for  a in wa
            for  cf in 1:length(confs)
                ks = keys(confs[cf].shells)
                if  a in ks   push!(shells, a);    break    end
            end
        end

        return( shells )
    end

    """
    `Basics.generate("subshells: ordered list for relativistic configurations", confs::Array{ConfigurationR,1})`
        ... to generate for confs, i.e. all the given (relativistic) configurations, common and ordered subshell list;
            a list::Array{Subshell,1} is returned.
    """
    function Basics.generate(sa::String, confs::Array{ConfigurationR,1})
        subshells = Subshell[]

        !(sa == "subshells: ordered list for relativistic configurations")  &&   error("Unsupported keystring = $sa")

        for  conf in confs
            for  subsh in keys(conf.subshells)      push!(subshells, subsh)     end
        end
        subshells = Base.unique(subshells)
        subshells = Base.sort( subshells, lt=Base.isless)
        ## Do include 'empty' subshells into the subshell list if they are specified by the given configurations
        ## if  a in ks   &&   confs[cf].subshells[a]  !=  0     push!(subshells, a);    break    end

        #== wa = Defaults.getDefaults("ordered subshell list: relativistic", 7)
        for  a in wa
            for  cf in 1:length(confs)
                ks = keys(confs[cf].subshells)
                ## Do include 'empty' subshells into the subshell list if they are specified by the given configurations
                ## if  a in ks   &&   confs[cf].subshells[a]  !=  0     push!(subshells, a);    break    end
                if  a in ks     push!(subshells, a);    break    end
            end
        end ==#

        return( subshells )
    end

    """
    `Basics.generate("subshells: ordered list for two bases", basisA::Basis,  basisB::Basis)`
        ... to generate common and ordered subshell list for the two basis A and B; a list::Array{Subshell,1} is returned.
    """
    function Basics.generate(sa::String, basisA::Basis,  basisB::Basis)
        function areEqual(nn::Int64, sha::Array{Subshell,1}, shb::Array{Subshell,1})
            # Determines whether the first nn subshells are equal in sha and shb (true) or not (false)
            for  i = 1:nx
                if    sha[i] != shb[i]    return( false )   end
            end
            return( true )
        end

        subshells = Subshell[]

        !(sa == "subshells: ordered list for two bases")  &&   error("Unsupported keystring = $sa")

        nx = min(length(basisA.subshells), length(basisB.subshells))
        if  areEqual(nx, basisA.subshells, basisB.subshells)
            # If subshell order is equal, any subshell order is accepted for those subshells that occur in both basis
            for  i = 1:nx
                if    basisA.subshells[i] != basisB.subshells[i]   error("Inconsistent subshells of two bases.")
                else  push!( subshells, basisA.subshells[i])
                end
            end
            #
            if       length(basisA.subshells) > nx
                for  i = nx+1:length(basisA.subshells)    push!(subshells, basisA.subshells[i])    end
            elseif   length(basisB.subshells) > nx
                for  i = nx+1:length(basisB.subshells)    push!(subshells, basisB.subshells[i])    end
            end
        else
            println("basisA.subshells = $(basisA.subshells)")
            println("basisB.subshells = $(basisB.subshells)")
            # If subshell order is NOT equal, all subshell in basisA and basisB must follow standard order
            standardList = Defaults.getDefaults("ordered subshell list: relativistic", 7)
            na = 0;  nb = 0
            for  sh in standardList
                if  sh in basisA.subshells  ||   sh in basisB.subshells   push!( subshells, sh)    end
            end
            nn = 0;
            for sh in basisA.subshells
                if    !(sh in subshells) push!( subshells, sh)
                else  wb = findall(x->x==sh, subshells);    if wb[1] <= nn   error("stop a")   else   nn = wb[1]  end
                end
            end
            nn = 0;
            for sh in basisB.subshells
                if    !(sh in subshells) push!( subshells, sh)
                else  wb = findall(x->x==sh, subshells);    if wb[1] <= nn   error("stop a")   else   nn = wb[1]  end
                end
            end
            println("*** extended subshells from two basis = $subshells")
        end

        return( subshells )
    end

    """
    `Basics.generate("single-electron spectrum: STO", N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0,
                     beta_0::Float64=1.1)`
        ... to generate a complete one-electron spectrum with N positive and N negative states, and by using even-tempered Slater-type
            orbitals (STO) with parameters ``\alpha_i = \alpha_0 \beta_0^i``; a spectrum::SingleElecSpectrum is returned where just
            N0 positive and N_0 negative are kept for later use.  **Not yet implemented !**

    `Basics.generate("single-electron spectrum: STO, positive", N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0,
                     beta_0::Float64=1.1)`
        ... to generate the same but to return only the N_0 positive states.  **Not yet implemented !**
    """
    function Basics.generate(sa::String, N::Int64, potential::Radial.Potential, grid::Radial.Grid; N_0::Int64=30, alpha_0::Float64=1.0, beta_0::Float64=1.1)
        !(sa == "single-electron spectrum: STO"  &&  sa == "single-electron spectrum: STO, positive")   &&   error("Unsupported keystring = $sa")
        error("Not yet implemented !")

        #= * define N normalized states in stoplus::Array{Vector{Float64},1} and stominus::Array{Vector{Float64},1}
            on the given grid and with given parameters
        * calculate one-electron (NxN) Hamiltonian matrix with potential h_pq and overlap S_pq, p,q = 1..N
        * solve generalized eigenvalue problem h c_vec = epsilon S c_vec
        * set spectrum::SingleElecSpectrum with just N0 functions
        * check normalization and orthogonality of the functions   =#

        return( nothing )
    end

    """
    `Basics.generateBasis(confList::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1})`
        ... generates the CSF basis for the given configuration list and level symmetries; a basis::Basis is returned but
            without a valid representation of the radial orbitals.
    """
    function Basics.generateBasis(confList::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1})
        #
        relconfList = ConfigurationR[]
        for  conf in confList
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)

        # 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)
            for  csf in newCsfs     if  LevelSymmetry(csf.J, csf.parity) in symmetries   push!(csfList, csf)   end   end
        end
        #
        if length(csfList) == 0     error("There are no CSF with $symmetries in the given configuration list.")      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

        basis = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, Dict{Subshell, Orbital}())

        println("Construct a basis with $(length(basis.csfs)) CSF for J^P = $symmetries with $(length(basis.subshells)) subshells: " *
                "$(basis.subshells[1])  $(basis.subshells[2]) ...  $(basis.subshells[end-1])  $(basis.subshells[end])")
        return( basis )
    end

    """
    `Basics.generateBasis(refConfigs::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1}, step::AtomicState.RasStep)`
        ... generates the CSF basis for the given reference configurations and single, double, ... excitations of
            electrons from the corresponding fromShells --> toShells. A basis::Basis is returned but without a valid
            representation of the radial orbitals
    """
    function Basics.generateBasis(refConfigs::Array{Configuration,1}, symmetries::Array{LevelSymmetry,1}, step::AtomicState.RasStep)
        # Single excitations
        if      step.seFrom == Shell[]  ||  step.seTo == Shell[]    confSingles = Configuration[]
        else    confSingles = Basics.generateConfigurations(refConfigs, step.seFrom, step.seTo)
        end
        # Double excitations
        if      step.deFrom == Shell[]  ||  step.deTo == Shell[]    confDoubles = Configuration[]
        else    confDoubles = Basics.generateConfigurations(refConfigs, step.deFrom, step.deTo)
                confDoubles = Basics.generateConfigurations(confDoubles,     step.deFrom, step.deTo)
        end
        # Triple excitations
        if      step.teFrom == Shell[]  ||  step.teTo == Shell[]    confTriples = Configuration[]
        else    confTriples = Basics.generateConfigurations(refConfigs,      step.teFrom, step.teTo)
                confTriples = Basics.generateConfigurations(confTriples,     step.teFrom, step.teTo)
                confTriples = Basics.generateConfigurations(confTriples,     step.teFrom, step.teTo)
        end
        # Quadruple excitations
        if      step.teFrom == Shell[]  ||  step.teTo == Shell[]    confQuadruples = Configuration[]
        else    confQuadruples = Basics.generateConfigurations(refConfigs, step.qeFrom, step.qeTo)
                confQuadruples = Basics.generateConfigurations(confQuadruples,  step.qeFrom, step.qeTo)
                confQuadruples = Basics.generateConfigurations(confQuadruples,  step.qeFrom, step.qeTo)
                confQuadruples = Basics.generateConfigurations(confQuadruples,  step.qeFrom, step.qeTo)
        end

        # Now get a unique set of configurations and generate the relativistic configurations
        # configurations = Base.unique((==), configurations) does not work here
        confList       = Configuration[]
        configurations = deepcopy(refConfigs);      append!(configurations, confSingles);    append!(configurations, confDoubles)
                                                    append!(configurations, confTriples);    append!(configurations, confQuadruples)
        for  confa in configurations
            addTo = true
            for confb in confList    if   confa == confb    addTo = false;    break     end     end
            if  addTo    push!(confList, confa)     end
        end
        for  i = 1:length(confList)    println(">> include ", confList[i])    end
        #
        relconfList = ConfigurationR[]
        for  conf in confList
            wa = Basics.generate("configuration list: relativistic", conf)
            append!( relconfList, wa)
        end
        subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)

        # 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)
            for  csf in newCsfs     if  LevelSymmetry(csf.J, csf.parity) in symmetries   push!(csfList, csf)   end   end
        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

        basis = Basis(true, NoElectrons, subshellList, csfList, coreSubshellList, Dict{Subshell, Orbital}())

        println("Construct a basis with $(length(basis.csfs)) CSF for J^P = $symmetries with $(length(basis.subshells)) subshells: " *
                "$(basis.subshells[1])  $(basis.subshells[2]) ...  $(basis.subshells[end-1])  $(basis.subshells[end])")
        return( basis )
    end

    """
    `Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1})`
        ... generates all nonrelativistic configurations due to the (single) excitation of an electron fromShells --> toShells
            for all given reference configurations. A confList::Array{Configuration,1} is returned.
    """
    function Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1})
        confList = Configuration[]
        for  config in refConfigs
            for  fromShell in fromShells
                for  toShell in toShells
                    # Cycle through all shells of config to built-up a new configuration with a single excitation if fromShell occurs
                    if  haskey(config.shells, fromShell )
                        newShells = Dict{Shell,Int64}();   addShell = true;   addConfiguration = true
                        for (k,v) in config.shells
                            if        k == fromShell == toShell       newShells = Base.merge( newShells, Dict( k => v))
                                      addShell = false
                            elseif    k == fromShell    if  v-1 > 0   newShells = Base.merge( newShells, Dict( fromShell => v-1))  end
                            elseif    k == toShell
                                if    v+1 <= 2*(2*toShell.l + 1)      newShells = Base.merge( newShells, Dict( toShell => v+1))
                                      addShell = false
                                else  addConfiguration = false
                                end
                            else      newShells = Base.merge( newShells, Dict( k => v))
                            end
                            #
                            if  addShell                              newShells = Base.merge( newShells, Dict( toShell => 1))      end
                        end
                        if  addConfiguration   push!(confList, Configuration(newShells, config.NoElectrons ))          end
                    else    # No configuration to add if fromShell does not occur in the reference configuration
                    end
                end
            end
        end

        confList = Base.unique(confList)

        return( confList )
    end

    """
    `Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1},
                                   noex::Int64; restrictions::Array{AbstractConfigurationRestriction,1}=AbstractConfigurationRestriction[])`
        ... generates all nonrelativistic configurations with excitation of up to noex electrons fromShells --> toShells
            for all given reference configurations. Moreover, the list of restrictions (if any is given) is finally applied
            to restrict the configurations to a given set of limitations. It remains the reponsibility of the user to make sure that
            the given restrictions do not contradict each other and are consistent with what is to be achieved. The given set of
            restrictions can be easily extended if this need arises by the users. A confList::Array{Configuration,1} is returned.
    """
    function Basics.generateConfigurations(refConfigs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1},
                                           noex::Int64; restrictions::Array{AbstractConfigurationRestriction,1}=AbstractConfigurationRestriction[])
        if  noex > 0
            newConfigs = Basics.generateConfigurations(refConfigs, fromShells, toShells)
            newConfigs = Base.unique(newConfigs)
            for ne = 2:noex
                newConfigs = Basics.generateConfigurations(newConfigs, fromShells, toShells)
                newConfigs = Base.unique(newConfigs)
            end
        end

        # Now apply in turn all given restrictions, if any, and append if no restriction is violated
        confList = Configuration[]
        for  conf  in  newConfigs
            addConf = true
            for res in restrictions
                if  Basics.isViolated(conf,res)     addConf = false;    break   end
            end
            if  addConf     push!(confList, conf)   end
        end

        return( confList )
    end

    """
    `Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.NoExcitationScheme,
                                                      nMax::Int64, lValues::Array{Int64,1})`
        ... generates a list of non-relativistic configurations for the given (reference) confs and the excitation scheme.
            All orbitals in standard order are considered up to maxShell = (n_max, l_max).
    """
    function Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.NoExcitationScheme,
                                                            nMax::Int64, lValues::Array{Int64,1})
        @warn(  "No excitations are included if exScheme::NoExcitationScheme. ")
        println("No excitations are included if exScheme::NoExcitationScheme. ")
        #
        newConfList = confs
        return( newConfList )
    end

    """
    `Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.DeExciteSingleElectron,
                                                      nMax::Int64, lValues::Array{Int64,1})`
        ... generates a list of non-relativistic configurations for the given (reference) confs and the excitation scheme.
            All orbitals in standard order are considered up to maxShell = (n_max, l_max).
    """
    function Basics.generateConfigurationsForExcitationScheme(confs::Array{Configuration,1}, exScheme::Basics.DeExciteSingleElectron,
                                                              nMax::Int64, lValues::Array{Int64,1})
        confList  = Configuration[];    NoElectrons = confs[1].NoElectrons
        shellList = Basics.generateShellList(confs, nMax, lValues)

        # Create all configurations with single excitations from the given list of configurations
        for  conf in confs
            # Take one electron fromShell and add toShell
            for  (fromShell,occ)  in  conf.shells
                for  toShell  in  shellList
                    newShells = deepcopy( conf.shells )
                    if      fromShell == toShell
                    elseif  haskey(conf.shells, toShell )   fromOcc = occ;   toOcc = conf.shells[toShell]
                            if  fromOcc - 1 < 0                        println("..");    continue    end
                            if  toOcc   + 1 > 2*(2*toShell.l + 1)      println("..");    continue    end
                            newShells[fromShell] = newShells[fromShell] - 1
                            newShells[toShell]   = newShells[toShell] + 1
                    else    fromOcc = occ
                            if  fromOcc - 1 < 0                        println("..");    continue    end
                            newShells[fromShell] = newShells[fromShell] - 1
                            newShells = Base.merge( newShells, Dict( toShell => 1))
                    end
                    # Add a new configuration
                    push!( confList, Configuration( newShells, NoElectrons))
                    if  true  println(">> Generate $(Configuration( newShells, NoElectrons))")   end
                end
            end
        end

        nbefore     = length(confList)
        newConfList = Configuration[]
        for  conf  in  confList
            if  conf in newConfList    continue;    else    push!( newConfList,  conf)    end
        end
        nafter      = length(newConfList)
        println(">> Number of generated configurations for $exScheme is: $nbefore (before) and $nafter (after).")

        return( newConfList )
    end

    """
    `Basics.generateConfigurationsWithElectronCapture(confs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1}, noex::Int64)`
        ... generates a list of non-relativistic configurations for the given (reference) confs and with one additional (cpatured)
            electron. All (doubly) excited configurations with upto :NoExcitations displacements of electrons fromShells into toShells
            and 'one' additional electron in the toShells are taken into account. The may result in large configuration lists even for
            a moderate number of fromShell and/or toShells.
    """
    function Basics.generateConfigurationsWithElectronCapture(confs::Array{Configuration,1}, fromShells::Array{Shell,1}, toShells::Array{Shell,1},
                                                              noex::Int64)
        newConfList = Configuration[];     NoElectrons = confs[1].NoElectrons + 1
        confList    = Basics.generateConfigurations(confs, fromShells, toShells, noex)
        # Now add one (captured) electron from the toShells to all configurations in confList
        for  conf in confList
            # Take one electron toShells and `add' it to conf
            for  toShell  in  toShells
                newShells = deepcopy( conf.shells )
                if      haskey(conf.shells, toShell )  &&  conf.shells[toShell]  + 1 > 2*(2*toShell.l + 1)     continue
                elseif  haskey(conf.shells, toShell )  newShells[toShell] = newShells[toShell] + 1
                        push!( newConfList, Configuration( newShells, NoElectrons))
                else    newShells = Base.merge( newShells, Dict( toShell => 1))
                        push!( newConfList, Configuration( newShells, NoElectrons))
                end
                if  false  println(">> Generate $(Configuration( newShells, NoElectrons)) with electron capture.")   end
            end
        end
        newConfList = unique(newConfList)

        return( newConfList )
    end

    """
    `Basics.generateLevelWithExtraElectron(newOrbital, symt::LevelSymmetry, level::Level)`
        ... generates a (new) level with one extra electron in subshell newOrbital.subshell and with overall symmetry symt. The function
            assumes that all CSF in the basis of level have the same symmetry as the level itself, that the new subshell is not yet part
            of the basis and that CSF with the total symmetry can be constructed; it terminates with an error if one of this assumptions
            is violated. A newLevel::Level with the same number of CSF and the same representation (eigenvector) as the given level is
            returned. From a physics viewpoint, newLevel refers to a antisymmetrized product state of [level x orbital(sh, energy)] J^P
            (symt) and with an total energy = level.energy + newOrbital.energy.
    """
    function Basics.generateLevelWithExtraElectron(newOrbital, symt::LevelSymmetry, level::Level)
        basis = level.basis;    newSubshells = copy(basis.subshells);    newCsfs = CsfR[];   J = level.J;   parity = level.parity

        push!(newSubshells, newOrbital.subshell)
        Defaults.setDefaults("relativistic subshell list", newSubshells; printout=false)
        newOrbitals = Base.merge( copy(basis.orbitals), Dict( newOrbital.subshell => newOrbital ))

        for  i = 1:length(basis.csfs)
            stateList   = ManyElectron.provideSubshellStates(newOrbital.subshell, 1)
            if  basis.csfs[i].J != J   ||   basis.csfs[i].parity != parity    error("Improper symmetry of CSF.")             end
            if  length(stateList) != 1                                        error("Improper number of subshell states.")   end
            if  AngularMomentum.triangularDelta(J, Basics.subshell_j(newOrbital.subshell), symt.J) != 1
                                                                            error("Improper coupling of (new) subshell & total J.")   end
            substate   = stateList[1]
            occupation = copy(basis.csfs[i].occupation);    push!(occupation, substate.occ)
            seniorityNr  = copy(basis.csfs[i].seniorityNr);     push!(seniorityNr,  substate.nu)
            subshellJ  = copy(basis.csfs[i].subshellJ);     push!(subshellJ,  AngularJ64(substate.Jsub2//2) )
            subshellX  = copy(basis.csfs[i].subshellX);     push!(subshellX,  symt.J)
            push!(newCsfs, CsfR(true, symt.J, symt.parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[] ) )
        end

        newBasis = Basis(true, basis.NoElectrons+1, newSubshells, newCsfs, copy(basis.coreSubshells), newOrbitals)
        newLevel = Level(symt.J, AngularM64(symt.J), symt.parity, level.index, level.energy+newOrbital.energy, level.relativeOcc,
                        true, newBasis, level.mc)
        return( newLevel )
    end

    """
    `Basics.generateLevelWithExtraSubshell(sh::Subshell, level::Level)`
        ... generates a (new) level with one extra subshell sh but with the same overall symmetry as before. The function assumes that the
            new subshell is not yet part of the basis; it terminates with an error if one of this assumptions is violated.
            A newLevel::Level with the same symmetry, number of CSF, energy and the same representation (eigenvector) as the given
            level is returned.
    """
    function Basics.generateLevelWithExtraSubshell(sh::Subshell, level::Level)
        basis = level.basis;    newSubshells = copy(basis.subshells);    newCsfs = CsfR[];   J = level.J;   parity = level.parity
        newOrbitals = copy(basis.orbitals)

        push!(newSubshells, sh)
        newOrbitals = Base.merge( newOrbitals, Dict( sh => Orbital(sh, -10000.) ))
        Defaults.setDefaults("relativistic subshell list", newSubshells; printout=false)

        for  i = 1:length(basis.csfs)
            stateList   = ManyElectron.provideSubshellStates(sh, 0)
            if  length(stateList) != 1                                        error("Improper number of subshell states.")   end
            substate   = stateList[1]
            occupation = copy(basis.csfs[i].occupation);    push!(occupation, substate.occ)
            seniorityNr  = copy(basis.csfs[i].seniorityNr);     push!(seniorityNr,  substate.nu)
            subshellJ  = copy(basis.csfs[i].subshellJ);     push!(subshellJ,  AngularJ64(substate.Jsub2//2) )
            subshellX  = copy(basis.csfs[i].subshellX);     push!(subshellX,  subshellX[end])
            push!(newCsfs, CsfR(true, basis.csfs[i].J, basis.csfs[i].parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[] ) )
        end

        newBasis = Basis(true, basis.NoElectrons, newSubshells, newCsfs, copy(basis.coreSubshells), newOrbitals)
        newLevel = Level(level.J, level.M, level.parity, level.index, level.energy, level.relativeOcc, true, newBasis, level.mc)
        return( newLevel )
    end

    """
    `Basics.generateLevelWithSymmetryReducedBasis(level::Level, subshells::Array{Subshell,1})`
        ... generates a level with a new basis and representation that only includes the CSF with the same symmetry as the level
            itself. It also ensures that the new basis is based on the given subshells so that it can be later used with levels
            which have one or several subshells more (at it often appears in ionization processes). Both, the underlying CSF basis
            and the eigenvectors are adopted accordingly. The procedures only supports the 'addition' of subshells but it stops
            if (1) the individual sequence differs or if the basis of level contains more subshells than the given subshell-list.
            A (new) level::levelR is returned.
    """
    function Basics.generateLevelWithSymmetryReducedBasis(level::Level, subshells::Array{Subshell,1})
        testBasis = level.basis;    newCsfs = CsfR[];    newMc = Float64[]

        # Decide of whether the basis need to be extended by some subshells
        nt = length(testBasis.subshells)
        if      testBasis.subshells == subshells      basis = testBasis;    nLevel = level
        elseif  nt >= length(subshells)               error("stop a")
        elseif  nt + 1 == length(subshells)  &&  testBasis.subshells == subshells[1:nt]
                nLevel = Basics.generateLevelWithExtraSubshell(subshells[end], level)
                basis  = nLevel.basis
        else
            nLevel  = level
            for  ns = nt+1:length(subshells)
                nLevel = Basics.generateLevelWithExtraSubshell(subshells[ns], nLevel)
            end
            basis  = nLevel.basis
        end

        for  i = 1:length(basis.csfs)
            if  basis.csfs[i].J == nLevel.J  &&   basis.csfs[i].parity == nLevel.parity        &&    abs(nLevel.mc[i]) > 1.0e-7
                push!(newCsfs, deepcopy(basis.csfs[i]));    push!(newMc, deepcopy(nLevel.mc[i]))
            end
        end
        newBasis = Basis(true, basis.NoElectrons, deepcopy(basis.subshells), newCsfs, deepcopy(basis.coreSubshells), deepcopy(basis.orbitals))
        newLevel = Level(nLevel.J, nLevel.M, nLevel.parity, nLevel.index, nLevel.energy, nLevel.relativeOcc, true, newBasis, newMc)
        return( newLevel )
    end

    """
    `Basics.generateOrbitalsForPotential(grid::Radial.Grid, meanPot::Radial.Potential, subshellList::Array{Subshell,1})`
        ... generates a set of (start) orbitals from the given potential and for all the subshells in subshellList.
            A set of orbitals::Dict{Subshell, Orbital} is returned.
    """
    function  Basics.generateOrbitalsForPotential(grid::Radial.Grid, meanPot::Radial.Potential, subshellList::Array{Subshell,1})
        orbitals = Dict{Subshell, Orbital}()

        # Determine kappaMin and kappaMax
        kappaMin = 0;   kappaMax = 0
        for  subsh in subshellList
            if       subsh.kappa < kappaMin    kappaMin = subsh.kappa
            elseif   subsh.kappa > kappaMax    kappaMax = subsh.kappa
            end
        end

        # Generate the primitives for a B-spline basis
        wa = Bsplines.generatePrimitives(grid)

        # Now cycle through all kappa symmetries
        for  kappa = kappaMin:kappaMax
            # Determine all requested subshells of symmetry kappa ... to compute them together
            shList     = Subshell[];    for  subsh in subshellList    if kappa == subsh.kappa    push!( shList, subsh)    end    end
            shOrbitals = Bsplines.generateOrbitalsForPotential(wa, kappa, shList, meanPot; printout=false)
            for  sh in shList
                orbitals = Base.merge( orbitals, Dict( sh => shOrbitals[sh]))
            end
        end

        return( orbitals )
    end

    """
    `Basics.generateOrbitalsForBasis(basis::Basis, frozenShells::Array{Shell,1}, priorBasis::Basis, startOrbitals::Dict{Subshell, Orbital})`
        ... generates a dict of (relativistic) orbitals as specificed by subshells of the basis. These orbitals are taken from priorBasis
            if contained in frozen-shells, and from spectrum (start orbitals) otherwises. An error message is issued if some
            requested orbital is not found. A orbList::Dict{Subshell, Orbital} is returned.
    """
    function Basics.generateOrbitalsForBasis(basis::Basis, frozenShells::Array{Shell,1}, priorBasis::Basis,
                                             startOrbitals::Dict{Subshell, Orbital})
        orbitals = Dict{Subshell, Orbital}()
        for  subsh in basis.subshells
            if  subsh in frozenShells    &&   haskey(priorBasis.orbitals, subsh)
                orbitals = Base.merge( orbitals, Dict( subsh => priorBasis.orbitals[subsh]) )
            elseif  subsh in frozenShells     error("Frozen orbital $subsh not found in prior basis.")
            elseif  haskey(priorBasis.orbitals, subsh)
                println(">> Start orbital $subsh is taken from prior basis")
                orbitals = Base.merge( orbitals, Dict( subsh => priorBasis.orbitals[subsh]) )
            else
                println(">> Start orbital $subsh is taken from hydrogenic orbitals")
                orbitals = Base.merge( orbitals, Dict( subsh => startOrbitals[subsh]) )
            end
        end
        return( orbitals )
    end

    """
    `Basics.generateOrbitalSuperposition(a::Orbital, b::Orbital, cx::Float64, grid::Radial.Grid)`
        ... generates a superposition  a + cx * b of two given orbitals; the function just takes the linear combination
            of the large and small components as well as the energy. The function assumes that both orbitals are defined on the same
            grid. A re-normalized newOrbital::Orbital is returned for which the derivatives are not defined.
            This function is used to accelerate the convergence (hopefully).
    """
    function Basics.generateOrbitalSuperposition(a::Orbital, b::Orbital, cx::Float64, grid::Radial.Grid)
        if  a.subshell != b.subshell  ||  !(a.isBound)  ||  !(b.isBound)   error("stop a")     end
        mtp  = max(size(a.P, 1), size(b.P, 1));     newP = zeros(mtp);  newQ = zeros(mtp);   newEnergy = (a.energy + cx * b.energy) / (1 + cx)
        mtpa = size(a.P, 1);    newP[1:mtpa] = a.P[1:mtpa];    newQ[1:mtpa] = a.Q[1:mtpa]
        mtpb = size(b.P, 1);    newP[1:mtpb] = newP[1:mtpb] + cx * b.P[1:mtpb]
                                newQ[1:mtpb] = newQ[1:mtpb] + cx * b.Q[1:mtpb]
        newOrbital = Orbital(a.subshell, a.isBound, a.useStandardGrid, newEnergy, newP, newQ, Float64[], Float64[], a.grid)
        norm       = RadialIntegrals.overlap(newOrbital, newOrbital, grid)
        newP       = newP / sqrt(norm);      newQ = newQ / sqrt(norm)
        newOrbital = Orbital(a.subshell, a.isBound, a.useStandardGrid, newEnergy, newP, newQ, Float64[], Float64[], a.grid)

        return( newOrbital )
    end

    """
    `Basics.generateShellList(confs::Array{Configuration,1}, nMax::Int64, lValues::Array{Int64,1})`
        ... generates a list of (non-relativistic) shells in standard order that contains all shells from the given list of
            configurations as well as the shells up to the principal quantum number nMax and the list of orbital quantum
            numbers lValues. A shellList::Array{Shell,1} is returned.
    """
    function Basics.generateShellList(confs::Array{Configuration,1}, nMax::Int64, lValues::Array{Int64,1})
        shellList = Shell[]
        # Add all shells from the given configurations
        for  conf in confs
            for (k,v) in conf.shells    push!( shellList, k)    end
        end
        # Add all shells for the given quantum numbers
        for  n = 1:nMax
            for  ll in  lValues
                if  n >= ll + 1     push!( shellList, Shell(n,ll))      end
            end
        end

        # Now bring the shells in standard order
        newShellList = Shell[]
        for  n = 1:nMax
            for  l = 0:30
                if  n >= l + 1  && Shell(n,l) in shellList    push!( newShellList, Shell(n,l))     end
            end
        end
        println(">> From configurations generated shell list $newShellList ")

        return( newShellList )
    end


text:
"""
    AbstractHMM{F<:VariateForm}

A custom HMM type must at-least implement the following interface:
```julia
struct CustomHMM{F,T} <: AbstractHMM{F}
    a::AbstractVector{T}               # Initial state distribution
    A::AbstractMatrix{T}               # Transition matrix
    B::AbstractVector{Distribution{F}} # Observations distributions
    # Optional, custom, fields ....
end
```
"""
abstract type AbstractHMM{F<:VariateForm} end

"""
    HMM([a, ]A, B) -> HMM

Build an HMM with transition matrix `A` and observation distributions `B`.
If the initial state distribution `a` is not specified, a uniform distribution is assumed.

Observations distributions can be of different types (for example `Normal` and `Exponential`),
but they must be of the same dimension.

**Arguments**
- `a::AbstractVector{T}`: initial probabilities vector.
- `A::AbstractMatrix{T}`: transition matrix.
- `B::AbstractVector{<:Distribution{F}}`: observations distributions.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
```
"""
struct HMM{F,T} <: AbstractHMM{F}
    a::Vector{T}
    A::Matrix{T}
    B::Vector{Distribution{F}}
    HMM{F,T}(a, A, B) where {F,T} = assert_hmm(a, A, B) && new(a, A, B)
end

HMM(a::AbstractVector{T}, A::AbstractMatrix{T}, B::AbstractVector{<:Distribution{F}}) where {F,T} = HMM{F,T}(a, A, B)
HMM(A::AbstractMatrix{T}, B::AbstractVector{<:Distribution{F}}) where {F,T} = HMM{F,T}(ones(size(A)[1])/size(A)[1], A, B)

"""
    assert_hmm(a, A, B)

Throw an `ArgumentError` if the initial state distribution and the transition matrix rows does not sum to 1,
and if the observation distributions do not have the same dimensions.
"""
function assert_hmm(a::AbstractVector,
                    A::AbstractMatrix,
                    B::AbstractVector{<:Distribution})
    @argcheck isprobvec(a)
    @argcheck istransmat(A)
    @argcheck all(length.(B) .== length(B[1])) ArgumentError("All distributions must have the same dimensions")
    @argcheck length(a) == size(A,1) == length(B)
    return true
end

"""
    issquare(A) -> Bool

Return true if `A` is a square matrix.
"""
issquare(A::AbstractMatrix) = size(A,1) == size(A,2)

"""
    istransmat(A) -> Bool

Return true if `A` is square and its rows sums to 1.
"""
istransmat(A::AbstractMatrix) = issquare(A) && all([isprobvec(A[i,:]) for i in 1:size(A,1)])

==(h1::AbstractHMM, h2::AbstractHMM) = (h1.a == h2.a) && (h1.A == h2.A) && (h1.B == h2.B)

"""
    rand([rng, ]hmm, T; init, seq) -> Array | (Vector, Array)

Sample a trajectory of `T` timesteps from `hmm`.

**Keyword Arguments**
- `init::Integer = rand(Categorical(hmm.a))`: initial state.
- `seq::Bool = false`: whether to return the hidden state sequence or not.

**Output**
- `Vector{Int}` (if `seq == true`): hidden state sequence.
- `Vector{Float64}` (for `Univariate` HMMs): observations (`T`).
- `Matrix{Float64}` (for `Multivariate` HMMs): observations (`T x dim(obs)`).

**Examples**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, 1000) # or
z, y = rand(hmm, 1000, seq = true)
size(y) # (1000,)
```

```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [MvNormal(ones(2)), MvNormal(ones(2))])
y = rand(hmm, 1000) # or
z, y = rand(hmm, 1000, seq = true)
size(y) # (1000, 2)
```
"""
function rand(rng::AbstractRNG, hmm::AbstractHMM, T::Integer; init = rand(rng, Categorical(hmm.a)), seq = false)
    z = Vector{Int}(undef, T)
    (T >= 1) && (z[1] = init)
    for t in 2:T
        z[t] = rand(rng, Categorical(hmm.A[z[t-1],:]))
    end
    y = rand(rng, hmm, z)
    seq ? (z, y) : y
end

"""
    rand([rng, ]hmm, z) -> Array

Sample observations from `hmm` according to trajectory `z`.

**Output**
- `Vector{Float64}` (for `Univariate` HMMs): observations (`T`).
- `Matrix{Float64}` (for `Multivariate` HMMs): observations (`T x dim(obs)`).

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, [1, 1, 2, 2, 1])
```
"""
function rand(rng::AbstractRNG, hmm::AbstractHMM{Univariate}, z::AbstractVector{<:Integer})
    y = Vector{Float64}(undef, length(z))
    for t in eachindex(z)
        y[t] = rand(rng, hmm.B[z[t]])
    end
    y
end

function rand(rng::AbstractRNG, hmm::AbstractHMM{Multivariate}, z::AbstractVector{<:Integer})
    y = Matrix{Float64}(undef, length(z), size(hmm, 1))
    for t in eachindex(z)
        y[t,:] = rand(rng, hmm.B[z[t]])
    end
    y
end

rand(hmm::AbstractHMM, T::Integer; kwargs...) = rand(GLOBAL_RNG, hmm, T; kwargs...)

rand(hmm::AbstractHMM, z::AbstractVector{<:Integer}) = rand(GLOBAL_RNG, hmm, z)

"""
    size(hmm, [dim]) -> Int | Tuple

Return the number of states in `hmm` and the dimension of the observations.

**Example**
```jldoctest
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
size(hmm)
# output
(2, 1)
```
"""
size(hmm::AbstractHMM, dim=:) = (length(hmm.B), length(hmm.B[1]))[dim]

"""

    copy(hmm) -> HMM

Return a copy of `hmm`.
"""
copy(hmm::HMM) = HMM(copy(hmm.a), copy(hmm.A), copy(hmm.B))

"""
    permute(hmm, perm) -> HMM

Permute the states of `hmm` according to `perm`.

**Arguments**

- `perm::Vector{<:Integer}`: permutation of the states.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.8 0.2; 0.1 0.9], [Normal(0,1), Normal(10,1)])
hmm = permute(hmm, [2, 1])
hmm.A # [0.9 0.1; 0.2 0.8]
hmm.B # [Normal(10,1), Normal(0,1)]
```
"""
function permute(hmm::AbstractHMM, perm::Vector{<:Integer})
    a = hmm.a[perm]
    B = hmm.B[perm]
    A = zeros(size(hmm.A))
    for i in 1:size(A,1), j in 1:size(A,2)
        A[i,j] = hmm.A[perm[i],perm[j]]
    end
    HMM(a, A, B)
end

"""
    nparams(hmm) -> Int

Return the number of _free_ parameters in `hmm`, without counting the
observation distributions parameters.

**Example**
```jldoctest
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
nparams(hmm)
# output
3
```
"""
function nparams(hmm::AbstractHMM)
    (length(hmm.a) - 1) + (length(hmm.A) - size(hmm.A, 1))
end

"""
    statdists(hmm) -> Vector{Vector}

Return the stationnary distribution(s) of `hmm`.
That is, the eigenvectors of transpose(hmm.A) with eigenvalues 1.
"""
function statdists(hmm::AbstractHMM)
    eig = eigen(collect(transpose(hmm.A)))
    dists = []
    for (i, val) in enumerate(eig.values)
        if val == 1.0
            dist = eig.vectors[:,i]
            dist /= sum(dist)
            push!(dists, dist)
        end
    end
    dists
end

ext: jl


text:
module Qutilities

using LinearAlgebra: Diagonal, eigvals, Hermitian, svdvals, tr

export
    binent,
    concurrence,
    concurrence_lb,
    formation,
    mutinf,
    negativity,
    ptrace,
    ptranspose,
    purity,
    sigma_x,
    sigma_y,
    sigma_z,
    spinflip,
    S_renyi,
    S_vn

# All logarithms are base 2.
const LOG = log2

"""
    hermitize(rho::AbstractMatrix)

Make `rho` Hermitian, but carefully.
"""
function hermitize(rho::AbstractMatrix)
    size(rho, 1) == size(rho, 2) || throw(DomainError(size(rho), "Only square matrices are supported."))

    err = maximum(abs.(rho - rho'))
    if err > 1e-12
        @warn "Matrix is not Hermitian: $(err)"
    end

    # Make the diagonal strictly real.
    rho_H = copy(rho)
    for i in 1:size(rho_H, 1)
        rho_H[i, i] = real(rho_H[i, i])
    end

    Hermitian(rho_H)
end

hermitize(rho::Diagonal) = rho

"""
    nonneg(x::Real)

`x` if `x` is non-negative; zero if `x` is negative.
"""
nonneg(x::Real) = x < 0 ? zero(x) : x

"""
    shannon(xs::AbstractVector)

Shannon entropy of `xs`.
"""
shannon(xs::AbstractVector) = -sum([x * LOG(x) for x in xs if x > 0])

# Single-qubit Pauli matrices.
const sigma_x = [[0.0, 1.0] [1.0, 0.0]]
const sigma_y = [[0.0, im] [-im, 0.0]]
const sigma_z = [[1.0, 0.0] [0.0, -1.0]]

"""
    ptrace(rho::AbstractMatrix{T}, dims, which::Int)

Partial trace of `rho` along dimension `which` from `dims`.
"""
function ptrace(rho::AbstractMatrix{T}, dims, which::Int) where {T}
    size(rho) == (prod(dims), prod(dims)) || throw(DomainError(size(rho), "Only square matrices are supported."))

    size_before = prod(dims[1:(which-1)])
    size_at = dims[which]
    size_after = prod(dims[(which+1):end])

    result = zeros(T, size_before*size_after, size_before*size_after)
    for i1 in 1:size_before
        for j1 in 1:size_before
            for k in 1:size_at
                for i2 in 1:size_after
                    for j2 in 1:size_after
                        row1 = size_after*(i1-1) + i2
                        col1 = size_after*(j1-1) + j2
                        row2 = size_at*size_after*(i1-1) + size_after*(k-1) + i2
                        col2 = size_at*size_after*(j1-1) + size_after*(k-1) + j2
                        result[row1, col1] += rho[row2, col2]
                    end
                end
            end
        end
    end
    result
end

"""
    ptrace(rho::AbstractMatrix, which::Int=2)

Partial trace of `rho` along dimension `which`.

`rho` is split into halves and the trace is over the second half by default.
"""
function ptrace(rho::AbstractMatrix, which::Int=2)
    size(rho, 1) % 2 == 0 || throw(DomainError(size(rho, 1), "Matrix size must be even."))

    s = div(size(rho, 1), 2)
    ptrace(rho, (s, s), which)
end

"""
    ptranspose(rho::AbstractMatrix, dims, which::Int)

Partial transpose of `rho` along dimension `which` from `dims`.
"""
function ptranspose(rho::AbstractMatrix, dims, which::Int)
    size(rho) == (prod(dims), prod(dims)) || throw(DomainError(size(rho), "Only square matrices are supported."))

    size_before = prod(dims[1:(which-1)])
    size_at = dims[which]
    size_after = prod(dims[(which+1):end])

    result = similar(rho)
    for i1 in 1:size_before
        for j1 in 1:size_before
            for i2 in 1:size_at
                for j2 in 1:size_at
                    for i3 in 1:size_after
                        for j3 in 1:size_after
                            row1 = size_at*size_after*(i1-1) + size_after*(i2-1) + i3
                            col1 = size_at*size_after*(j1-1) + size_after*(j2-1) + j3
                            row2 = size_at*size_after*(i1-1) + size_after*(j2-1) + i3
                            col2 = size_at*size_after*(j1-1) + size_after*(i2-1) + j3
                            result[row1, col1] = rho[row2, col2]
                        end
                    end
                end
            end
        end
    end
    result
end

"""
    ptranspose(rho::AbstractMatrix, which::Int=2)

Partial transpose of `rho` along dimension `which`.

`rho` is split into halves and the transpose is over the second half by
default.
"""
function ptranspose(rho::AbstractMatrix, which::Int=2)
    size(rho, 1) % 2 == 0 || throw(DomainError(size(rho, 1), "Matrix size must be even."))

    s = div(size(rho, 1), 2)
    ptranspose(rho, (s, s), which)
end

"""
    binent(x::Real)

Binary entropy of `x`.
"""
binent(x::Real) = shannon([x, one(x) - x])

"""
    purity(rho::AbstractMatrix)

Purity of `rho`.
"""
purity(rho::AbstractMatrix) = rho^2 |> tr |> real

"""
    S_vn(rho::AbstractMatrix)

Von Neumann entropy of `rho`.
"""
S_vn(rho::AbstractMatrix) = rho |> hermitize |> eigvals |> shannon

"""
    S_renyi(rho::AbstractMatrix, alpha::Real=2)

Order `alpha` Rényi entropy of `rho`.

The `alpha` parameter may have any value on the interval `[0, Inf]` except 1.
It defaults to 2.
"""
function S_renyi(rho::AbstractMatrix, alpha::Real=2)
    E = rho |> hermitize |> eigvals
    alpha == Inf && return E |> maximum |> LOG |> -
    LOG(sum(E.^alpha)) / (1 - alpha)
end

"""
    mutinf(rho::AbstractMatrix, S::Function=S_vn)

Mutual information of `rho` using the entropy function `S`.
"""
function mutinf(rho::AbstractMatrix, S::Function=S_vn)
    S(ptrace(rho, 1)) + S(ptrace(rho, 2)) - S(rho)
end

"""
    spinflip(rho::AbstractMatrix)

Wootters spin-flip operation for two qubits in the mixed state `rho` in the
standard basis.

Ref: Wootters, W. K. (1998). Entanglement of formation of an arbitrary state of
two qubits. Physical Review Letters, 80(10), 2245.
"""
function spinflip(rho::AbstractMatrix)
    size(rho) == (4, 4) || throw(DomainError(size(rho), "Matrix must be 4 by 4."))

    Y = kron(sigma_y, sigma_y)
    Y * conj(rho) * Y
end

"""
    concurrence(rho::AbstractMatrix)

Concurrence of a mixed state `rho` for two qubits in the standard basis.

Ref: Wootters, W. K. (1998). Entanglement of formation of an arbitrary state of
two qubits. Physical Review Letters, 80(10), 2245.
"""
function concurrence(rho::AbstractMatrix)
    size(rho) == (4, 4) || throw(DomainError(size(rho), "Matrix must be 4 by 4."))

    E = rho * spinflip(rho) |> eigvals
    if any(imag(E) .> 1e-15)
        @warn "Complex eigenvalues: $(maximum(imag(E)))"
    end
    if any(real(E) .< -1e-12)
        @warn "Negative eigenvalues: $(minimum(real(E)))"
    end
    F = sort(sqrt.(nonneg.(real(E))))
    nonneg(F[end] - sum(F[1:(end-1)]))
end

"""
    concurrence_lb(rho::AbstractMatrix)

Lower bound on the concurrence for two qubits in the mixed state `rho` in the
standard basis.

Ref: Mintert, F., & Buchleitner, A. (2007). Observable entanglement measure for
mixed quantum states. Physical Review Letters, 98(14), 140505.
"""
function concurrence_lb(rho::AbstractMatrix)
    size(rho) == (4, 4) || throw(DomainError(size(rho), "Matrix must be 4 by 4."))

    2.0*(purity(rho) - purity(ptrace(rho))) |> nonneg |> sqrt
end

"""
    formation(C::Real)

Entanglement of formation for two qubits with concurrence `C`.

Ref: Wootters, W. K. (1998). Entanglement of formation of an arbitrary state of
two qubits. Physical Review Letters, 80(10), 2245.
"""
formation(C::Real) = binent(0.5 * (1.0 + sqrt(1.0 - C^2)))

"""
    formation(rho::AbstractMatrix)

Entanglement of formation for two qubits in the mixed state `rho` in the
standard basis.

Ref: Wootters, W. K. (1998). Entanglement of formation of an arbitrary state of
two qubits. Physical Review Letters, 80(10), 2245.
"""
formation(rho::AbstractMatrix) = rho |> concurrence |> formation

"""
    negativity(rho::AbstractMatrix)

Logarithmic negativity for a symmetric bipartition of `rho`.

Ref: Plenio, M. B. (2005). Logarithmic negativity: A full entanglement monotone
that is not convex. Physical Review Letters, 95(9), 090503.
"""
negativity(rho::AbstractMatrix) = rho |> ptranspose |> svdvals |> sum |> LOG

end


text:
export IntervalDomain, ProductDomain, ⊗, CircleDomain

abstract type AbstractDomain{T,N} end

struct VarDomainPairing
  variables
  domain::AbstractDomain
end
Base.:∈(variable::ModelingToolkit.Operation,domain::AbstractDomain) = VarDomainPairing(variable,domain)
Base.:∈(variables::NTuple{N,ModelingToolkit.Operation},domain::AbstractDomain) where N = VarDomainPairing(variables,domain)

## Specific Domains

struct IntervalDomain{T} <: AbstractDomain{T,1}
  lower::T
  upper::T
end

struct ProductDomain{D,T,N} <: AbstractDomain{T,N}
  domains::D
end
⊗(args::AbstractDomain{T}...) where T = ProductDomain{typeof(args),T,length(args)}(args)

struct CircleDomain <: AbstractDomain{Float64,2}
  polar::Bool
  CircleDomain(polar=false) = new(polar)
end


text:
using MAT
using Datetime
using JSON

include("config.jl")
include("utils.jl")

jobs = Dict[]
resultStreams = Dict()
working_on = Dict{Int, RemoteRef}()
idle = Int[]

function runCluster(fileName, outFolder = "results"; configName = "cluster_configs/all.json", loadTime=-1)
    try
        mkdir(outFolder)
    end

    dt = now()
    out = "$outFolder/$(year(dt))-$(month(dt))-$(day(dt))_$(hour(dt)):$(minute(dt)):$(second(dt))"
    mkdir(out)

    config = matread("data/$(fileName).mat")

    jobs, resultStreams = parseClusterConfig(configName, out)

    # Setup everything
    nc, ng, c_ocl = determineWorkload()

    addprocs(nc + ng)
    idle = workers()
    require("job.jl")

    count = 1
    for id in workers()
        if count <= nc
            if c_ocl
                @at_proc id CPU_OCL = true
            else
                @at_proc id CPU_OCL = false
            end
            @at_proc id GPU_OCL = false
        else
            @at_proc id CPU_OCL = false
            @at_proc id GPU_OCL = true
        end
        count += 1
    end

    println("Starting simulations on $(nworkers()) workers")

    try
        running = true
        while running
            for (w, rref) in working_on
                if isready(rref)
                    rref = pop!(working_on, w) # remove from dict
                    try
                        fN, params, result = fetch(rref) # fetch Remote Reference
                        writedlm_row(resultStreams[fN], [params, result...])
                    catch e
                        warn("Could not write $(fetch(rref)) because of $e")
                    end
                    push!(idle, w) # mark worker as idle
                end
            end

            while (length(idle) > 0) && (length(jobs) > 0)
                work_item = pop!(jobs) # get a work item
                w = pop!(idle) # get an idle worker
                rref = @spawnat w runProcess(config, work_item, out, loadTime)
                working_on[w] = rref
            end

            println("Running jobs: $(length(working_on))")
            println("Jobs left: $(length(jobs))")

            if (length(jobs) == 0) && (length(working_on) == 0)
                running = false
            else
		sleep(30.0)
	    	# timedwait(1.0, pollint=0.5) do
		#	rrefs = values(working_on)
		#	isempty(rrefs) || any(isready, rrefs)
		# end
            end
        end
    catch e
        println(e)
    finally
        map(close, values(resultStreams))
    end
end

function parseClusterConfig(filePath, outputPath)
    config = JSON.parse(open(filePath, "r"))
    vals = (String, Dict)[]
    outputStreams = Dict()

    for (fileName,  subconfig) in config
        if haskey(subconfig, "name")
            name = subconfig["name"]
            parameterNames, subvals =
                if name == "punch_local"
                    parsePunchLocal(subconfig)
                elseif name == "punch_random"
                    parsePunchRandom(subconfig)
                elseif name == "gaussian_blur"
                    parseGaussianBlur(subconfig)
                elseif name == "global"
                    parseGlobal(subconfig)
                else
                    warn("Can't handle $name")
                    Dict[]
                end
            times = if haskey(subconfig, "times")
                        subconfig["times"]
                    else
                        1
                    end
            subvals = [(fileName, v) for v in subvals]
            for i in 1:times
                append!(vals, subvals)
            end
            out = open("$outputPath/$fileName", "w")
            writedlm_row(out, ["name", parameterNames, "time", "timeToStable", "valueOfInterest", "stable", "meanMField", "meanAField", "structM", "structA", "structF", "structW", "structD", "fileName"])
            outputStreams[fileName] = out
        else
            warn("$fileName has no disturbance assign to it.")
        end
    end
    return (vals, outputStreams)
end

function parsePunchLocal(config)
    x = config["x"]
    y = config["y"]
    alpha = config["alpha"]
    beta = config["beta"]

    ranges = map(torange, (x, y, alpha, beta))
    (["x", "y", "alpha", "beta"], createDisturbance(:punch_local, ranges))
end

function parsePunchRandom(config)
    alpha = config["alpha"]
    beta = config["beta"]
    ranges = map(torange, (alpha, beta))

    (["alpha", "beta"], createDisturbance(:punch_random, ranges))
end

function parseGaussianBlur(config)
    sigma = config["sigma"]
    (["sigma"], createDisturbance1(:gaussian_blur, torange(sigma)))
end

function parseGlobal(config)
    mu = config["mu"]
    sig = config["sig"]

    ranges = map(torange, (mu, sig))
    (["mu", "sig"], createDisturbance(:global, ranges))
end


text:
@safetestset "Router tests" begin

  @safetestset "Basic routing" begin
    using Genie, Genie.Router

    route("/hello") do
      "Hello"
    end

  end;

end;


text:
function crawl(::Val{:ref}, args...)
    return :(dereference($(crawl.(args)...)))
end
function deparse(::Val{:index}, b::JSNode, idx)::JSString
    return jsstring(deparse(b), "[", deparse(idx), "]")
end

"""
Generate a JSNode should be inserted when a value is dereferences.

This function is called with the value that is to be dereferenced and not the
expression itself that is being dereferenced (i.e., after the crawl phase).

This is used mainly to allow WebIO to customize behavior when dereferencing
Observables.
"""
function dereference(args...)
    return JSAST(:index, args...)
end


text:
# @testset "integrate" begin

#     # Test basic properties of the prior.
#     let f = GP(ZeroMean(), EQ(), GPC())
#         v = ∫(StandardNormal(), f)
#         @test mean(v)(1) == 0.0
#         @test mean_vector(v) == [0.0]
#         @test kernel(v)(1, 1) == 1 / sqrt(3)
#         @test kernel(v, f)(1, 5.0) == kernel(f, v)(5.0, 1)
#         @test kernel(v, f)(1, 0.0) == 1 / sqrt(2)
#     end

#     # Test properties of the posterior when the integral is observed.
#     let f = GP(ZeroMean(), EQ(), GPC())
#         v = ∫(StandardNormal(), f)
#         v′, f′ = (v | (v←[1])), (f | (v←[1]))
#         x = collect(range(-10.0, stop=10.0, length=100))
#         g(x) = exp(-0.5 * x^2) / sqrt(2π)
#         @test abs(quadgk(x->mean(f′)(x) * g(x), -10.0, 10.0)[1] - 1.0) < 1e-6
#     end

#     # Test properties of the posterior when the function is observed.
#     let rng = MersenneTwister(123456), f = GP(ZeroMean(), EQ(), GPC())
#         v = ∫(StandardNormal(), f)
#         x = collect(range(-10.0, stop=10.0, length=100))
#         fs = sample(rng, f(x))
#         v′, μf′ = mean(v | (f(x)←fs))(1), mean(f | (f(x)←fs))
#         g(x) = exp(-0.5 * x^2) / sqrt(2π)
#         @test abs(quadgk(x->μf′(x) * g(x), -10.0, 10.0)[1] - v′) < 1e-6
#     end
# end


text:
"""
    search(command::String, logicset::LogicalCombo; verbose=true)

Searches for a possible match among a LogicalCombo in which the wildcard term
is true. Search requires the use of a wildcard. In the event that a wildcard
is missing, search will insert a {{i}} to the left of the command.*{{i+1}}* can
be used to search for relationships between the ith column and another
column.

### Examples
```julia
julia> myset = logicalparse("v1, v2, v3 ∈ 1:10")
v1, v2, v3 ∈ 1:10        feasible outcomes 1000 ✓        :6 6 10

julia> myset = logicalparse("{{i}} >> {{i+1}}", myset)
{{i}} >> {{i+1}}
>>> v1 >> v2
>>> v2 >> v3
         feasible outcomes 56 ✓          :10 7 3

julia> search("{{i}} == 4", myset)
Checking: v1 == 4
Checking: v2 == 4
Checking: v3 == 4

:v1 is a not match with 0 feasible combinations out of 56.
:v2 is a possible match with 10 feasible combinations out of 56.
:v3 is a possible match with 6 feasible combinations out of 56.

julia> search("== 4", myset, verbose=false) == search("{{i}} == 4", myset, verbose=false)
true

julia> search("{{i}} > {{!i}}", myset)
Checking: v1 > v2
Checking: v1 > v3
Checking: v2 > v1
Checking: v2 > v3
Checking: v3 > v1
Checking: v3 > v2

:v1 is a match with 56 feasible combinations out of 56.
:v2 is a not match with 0 feasible combinations out of 56.
:v3 is a not match with 0 feasible combinations out of 56.
```
"""
function search(command::String, logicset::LogicalCombo; verbose=true)

    clearcounter()

    rowsin = nfeasible(logicset)

    if rowsin == 0
        println("No outcomes feasible outcomes with set!")
        return missing
    end

    !occursin(raw"{{", command) &&
      !occursin(r"^[><=|!^&]", command) &&
      throw("No syntax match. Enter form '{{i}} = ...' or '= ...'")

    !occursin(raw"{{", command) &&
      occursin(r"^[><=|!^&]", command) &&
      (command = "{{i}} " * command)

    colcheck = operatorspawn(command,
      logicset,
      returnlogical=true,
      prefix="Checking: ",
      verbose=verbose)

    colcount = Any[]

    for i in 1:size(logicset, 2)
      select = colcheck[1][:, colcheck[2] .== logicset.keys[i]]
      (size(select,2) == 0) && push!(colcount, missing)
      (size(select,2) > 0 ) && push!(colcount, [all(select[j,:]) for j in 1:size(select,1)] |> sum)
    end

    # colcount = [sum(colcheck[:,i]) for i in 1:size(colcheck,2)]

    if verbose
        println()
        for i in 1:size(logicset,2)
          if colcount[i] === missing
              println(":$(logicset.keys[i]) could not be evaluated at $command.")
              continue
          end
          (colcount[i] == rowsin) && print(":$(logicset.keys[i]) is a match")
          (colcount[i] == 0) && print(":$(logicset.keys[i]) is a not match")
          (colcount[i] > 0) && (colcount[i] < rowsin) &&
            print(":$(logicset.keys[i]) is a possible match")
          println(" with $(colcount[i]) feasible combinations out of $rowsin.")
        end
    end

    return colcount ./ rowsin
end


text:
function part1(input)
    binaries = map(line -> collect(line) .== '1', readlines(input))
    majority = sum(binaries) .> length(binaries) / 2
    γ = decode(majority)
    ϵ = decode(.!majority)
    return γ * ϵ
end

function part2(input)
    binaries = map(line -> collect(line) .== '1', readlines(input))
    oxygen = decode(find_common(copy(binaries), 1, false))
    CO2 = decode(find_common(copy(binaries), 1, true))
    return oxygen * CO2
end

decode(x) = sum(x .<< (reverse(eachindex(x)) .- 1))

function find_common(binaries, n, choose_uncommon)
    length(binaries) == 1 && return only(binaries)
    common = sum(getindex.(binaries, n)) >= length(binaries) / 2
    return find_common(filter(x -> x[n] == common ⊻ choose_uncommon, binaries),
                       n + 1, choose_uncommon)
end


text:
using LogarithmicGreensFunctions
using Base.Test

using PyPlot

symset!(A,i,j,v) = (A[i,j] = A[j,i] = v)

# One interval

n = 100
H = full(Tridiagonal(-ones(n-1),zeros(n),-ones(n-1)))
symset!(H,n,1,-1)
d = [minimum(abs(i-j+s*n) for s = -1:1) for i = 1:n, j = 1:n]

figure()
for (i,z) in enumerate((3,im,0.1-0.1*im))
    semilogy(vec(d), vec(abs.(inv(H - z*I))), "C$(i-1)o", ms=2,mew=0);
    semilogy(0:n÷2, exp.(-greens(domain(-2,2),z)*(0:n÷2)), "C$(i-1)--")
end
show()

# Two intervals

H = full(Tridiagonal(-ones(n-1),repmat([-1,1],n÷2),-ones(n-1)))
symset!(H,n,1,-1)
d = [minimum(abs(i-j+s*n) for s = -1:1) for i = 1:n, j = 1:n]

figure()
for (i,z) in enumerate((0,0.5,3,2+0.1im))
    semilogy(vec(d), vec(abs.(inv(H - z*I))), "C$(i-1)o", ms=2,mew=0);
    semilogy(0:n÷2, exp.(-greens(domain(-sqrt(5),-1,1,sqrt(5)),z)*(0:n÷2)), "C$(i-1)--")
end
show()


text:
push!(LOAD_PATH, string(pwd(),"/src"))

using Base.Test
using ModifiedStokesSolver
import AnalyticDomains

@testset "BlockCorrections" begin
    # Params
    alpha = 10.0
    # Discretize
    numpanels = 20
    panelorder = 16
    curve = AnalyticDomains.starfish(n_arms = 3, amplitude=0.2)
    dcurve = discretize(curve, numpanels, panelorder)
    N = dcurve.numpoints

    # Assemble complete double layer matrix
    Dmat = doublelayer_matrix(dcurve, alpha)

    # Assemble block correction matrix
    BCmat = ModifiedStokesSolver.correction_block_matrix(dcurve, alpha)

    tovec(x) = reshape(x, :)
    fromvec(x) = reshape(x, 2, :)
    matvec(x) = tovec(doublelayer_self(dcurve, fromvec(x), alpha, zerolimit=false)) + BCmat*x

    x = rand(2*N)
    v1 = Dmat*x
    v2 = matvec(x)
    relerr = norm(v1-v2, Inf) / norm(v1, Inf)

    @test relerr < 1e-14
end


text:
# test of transforms

using NumericExtensions
using Base.Test

## linear transforms

x = rand(5)
X = rand(5, 4)

# scalar
a = 2.5
@test_approx_eq transform(a, x) a * x
@test_approx_eq transform(a, X) a * X

y = copy(x); transform!(a, y)
@test_approx_eq y a * x
Y = copy(X); transform!(a, Y)
@test_approx_eq Y a * X

# vector
a = rand(5) + 1.0
@test_approx_eq transform(a, x) a .* x
@test_approx_eq transform(a, X) a .* X

y = copy(x); transform!(a, y)
@test_approx_eq y a .* x
Y = copy(X); transform!(a, Y)
@test_approx_eq Y a .* X

# matrix
a = rand(3, 5) + 1.0

@test_approx_eq transform(a, x) a * x
@test_approx_eq transform(a, X) a * X

# transpose
a = rand(5, 3) + 1.0

@test_approx_eq transform(Transpose(a), x) a'x
@test_approx_eq transform(Transpose(a), X) a'X

## affine transforms

b = randn(5)

# scalar
a = 2.5
@test_approx_eq transform(AffineTransform(a), x) a * x
@test_approx_eq transform(AffineTransform(a), X) a * X
@test_approx_eq transform(AffineTransform(a, b), x) a * x .+ b
@test_approx_eq transform(AffineTransform(a, b), X) a * X .+ b

y = copy(x); transform!(AffineTransform(a), y)
@test_approx_eq y a * x
Y = copy(X); transform!(AffineTransform(a), Y)
@test_approx_eq Y a * X
y = copy(x); transform!(AffineTransform(a, b), y)
@test_approx_eq y a * x .+ b
Y = copy(X); transform!(AffineTransform(a, b), Y)
@test_approx_eq Y a * X .+ b

# vector
a = rand(5) + 1.0
@test_approx_eq transform(AffineTransform(a), x) a .* x
@test_approx_eq transform(AffineTransform(a), X) a .* X
@test_approx_eq transform(AffineTransform(a, b), x) a .* x .+ b
@test_approx_eq transform(AffineTransform(a, b), X) a .* X .+ b

y = copy(x); transform!(AffineTransform(a), y)
@test_approx_eq y a .* x
Y = copy(X); transform!(AffineTransform(a), Y)
@test_approx_eq Y a .* X
y = copy(x); transform!(AffineTransform(a, b), y)
@test_approx_eq y a .* x .+ b
Y = copy(X); transform!(AffineTransform(a, b), Y)
@test_approx_eq Y a .* X .+ b

# matrix
a = rand(3, 5)
b = rand(3)
@test_approx_eq transform(AffineTransform(a), x) a * x
@test_approx_eq transform(AffineTransform(a), X) a * X
@test_approx_eq transform(AffineTransform(a, b), x) a * x .+ b
@test_approx_eq transform(AffineTransform(a, b), X) a * X .+ b

# transpose
a = rand(5, 3)
b = rand(3)
@test_approx_eq transform(AffineTransform(Transpose(a)), x) a'x
@test_approx_eq transform(AffineTransform(Transpose(a)), X) a'X
@test_approx_eq transform(AffineTransform(Transpose(a), b), x) a'x .+ b
@test_approx_eq transform(AffineTransform(Transpose(a), b), X) a'X .+ b


text:
function visualise_pulses(init_ctrl_arr, final_ctrl_arr)
        fig_init = plot(vec(init_ctrl_arr), linetype=:steppre, label="initial guess", lc=:red, xlabel="time")
        fig_final = plot(vec(final_ctrl_arr), linetype=:steppre, label="optimal pulses", lc=:blue, xlabel="time")

        fig = plot(fig_init, fig_final, layout = (2, 1))
        savefig(fig, "fig.png")
end


text:
# Borrowed from Zygote tests
function ngradient(f, xs::AbstractArray...)
    grads = zero.(xs)
    for (x, Δ) in zip(xs, grads), i = 1:length(x)
        δ = sqrt(eps())
        tmp = x[i]
        x[i] = tmp - δ / 2
        y1 = f(xs...)
        x[i] = tmp + δ / 2
        y2 = f(xs...)
        x[i] = tmp
        Δ[i] = (y2 - y1) / δ
    end
    return grads
end

gradcheck(f, xs...) =
    all(isapprox.(ngradient(f, xs...), gradient(f, xs...), rtol = 1e-5, atol = 1e-5))

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)


text:
module test_julia_abstractarraymath_eachrow_eachcol_eachslice

using Test

@test collect(eachrow([1, 2, 3])) == [[1], [2], [3]]
@test collect(eachrow([1 2 3; ])) == [[1, 2, 3]]
@test collect(eachrow([])) == []
@test collect(eachrow([]')) == [[]]

@test collect(eachcol([1, 2, 3])) == [[1, 2, 3]]
@test collect(eachcol([1 2 3; ])) == [[1], [2], [3]]
@test collect(eachcol([])) == [[]]
@test collect(eachcol([]')) == []

@test collect(eachslice([1 2 3;], dims=1)) == [[1, 2, 3]]
@test collect(eachslice([1 2 3;], dims=2)) == [[1], [2], [3]]

end # module test_julia_abstractarraymath_eachrow_eachcol_eachslice

module test_julia_abstractarraymath_repeat_circshift

using Test

@test repeat([1; 2; 3], 1, 5) == [1 1 1 1 1;
                                  2 2 2 2 2;
                                  3 3 3 3 3;]

@test repeat([1 2 3;], 5, 1) == [1 2 3;
                                 1 2 3;
                                 1 2 3;
                                 1 2 3;
                                 1 2 3;]

@test circshift([1, 2, 3, 4, 5], 3) ==
                      [3, 4, 5, 1, 2]

end # module test_julia_abstractarraymath_repeat_circshift


text:
using DataFrames, TableOperations, Tables, Random
using TreeParzen
using Recommenders:
    Movielens1M,
    load_dataset,
    ratio_split,
    BPR,
    evaluate_u2i,
    MeanPrecision,
    MeanRecall,
    MeanNDCG

function one_write(filepath, params, result)
    n_epochs = convert(Int, params[:n_epochs])
    n_negatives = convert(Int, params[:n_negatives])
    learning_rate = params[:learning_rate]
    dimension = convert(Int, 2^params[:log2_dimension])
    reg_coeff = params[:reg_coeff]

    ndcg10 = result[:ndcg10]
    precision10 = result[:precision10]
    recall10 = result[:recall10]

    str = "$n_epochs,$n_negatives,$learning_rate,$dimension,$reg_coeff,$ndcg10,$precision10,$recall10\n"

    open(filepath, "a") do io
        write(io, str)
    end
end

function search(filepath)
    ml1M = Movielens1M()
    download(ml1M)
    rating, _, _ = load_dataset(ml1M)

    Random.seed!(1234)
    train_valid_table, test_table = ratio_split(rating, 0.8)
    train_table, valid_table = ratio_split(train_valid_table, 0.8)

    prec10 = MeanPrecision(10)
    recall10 = MeanRecall(10)
    ndcg10 = MeanNDCG(10)
    metrics = [prec10, recall10, ndcg10]

    function invert_output(params)
        @info params
        n_epochs = convert(Int, params[:n_epochs])
        n_negatives = convert(Int, params[:n_negatives])
        learning_rate = params[:learning_rate]
        dimension = convert(Int, 2^params[:log2_dimension])
        reg_coeff = params[:reg_coeff]

        model = BPR(dimension, reg_coeff)

        result = evaluate_u2i(
            model,
            train_table,
            valid_table,
            metrics,
            10,
            col_item = :movieid,
            n_epochs = n_epochs,
            n_negatives = n_negatives,
            learning_rate = learning_rate,
            drop_history = true,
            verbose = 16,
        )
        @info result
        one_write(filepath, params, result)
        return -result[:ndcg10]
    end

    open(filepath, "w") do io
        header = "#n_epochs,n_negatives,learning_rate,dimension,reg_coeff,ndcg10,precision10,recall10\n"
        write(io, header)
    end

    @info "coarse search"
    best_result = Inf
    best_params = nothing
    for lr in [1e-3, 3e-3, 1e-2], n_neg in [4, 8, 16]
        params = Dict(
            :n_epochs => 128,
            :n_negatives => n_neg,
            :learning_rate => lr,
            :log2_dimension => 6,
            :reg_coeff => 1e-20,
        )
        result = invert_output(params)
        if result < best_result
            best_result = result
            best_params = params
        end
    end

    for reg_coeff in [1e-3, 3e-3, 1e-2]
        params = copy(best_params)
        params[:reg_coeff] = reg_coeff
        result = invert_output(params)
        if result < best_result
            best_result = result
            best_params = params
        end
    end

    @info "refinement"
    space = Dict(
        :n_epochs =>
            HP.Choice(:n_epochs, [best_params[:n_epochs], best_params[:n_epochs] * 2]),
        :n_negatives => HP.QuantUniform(
            :n_negatives,
            best_params[:n_negatives] - 2.0,
            best_params[:n_negatives] + 2.0,
            1.0,
        ),
        :learning_rate => HP.Uniform(
            :learning_rate,
            best_params[:learning_rate] * 0.5,
            best_params[:learning_rate] * 2.0,
        ),
        :log2_dimension => HP.Choice(:log2_dimension, [6, 7, 8]),
        :reg_coeff => HP.Uniform(
            :reg_coeff,
            best_params[:reg_coeff] * 0.5,
            best_params[:reg_coeff] * 2.0,
        ),
    )

    best = fmin(invert_output, space, 20, logging_interval = -1)
    @info best

    @info "Evaluate best model."

    best_model = BPR(convert(Int, 2^best[:log2_dimension]), best[:reg_coeff])
    result = evaluate_u2i(
        best_model,
        train_valid_table,
        test_table,
        metrics,
        10,
        col_item = :movieid,
        n_epochs = convert(Int, best[:n_epochs]),
        n_negatives = convert(Int, best[:n_negatives]),
        learning_rate = best[:learning_rate],
        drop_history = true,
        verbose = 16,
    )
    @info result

end

search("./log_param_search_bpr_ml1M.txt")


text:
# Adaptive Metropolis Random Walk
# https://github.com/TuringLang/AdvancedMH.jl/blob/master/src/mh-core.jl

# NOTE:
# - Is there a more computationally efficient way?
# - Can we use bijectors to transform continuous parameters to unconstrained space?

import Pkg; Pkg.activate(".")
include("amh.jl")
using StatsPlots
using MCMCChains
import Random

Random.seed!(0)

trilindex(x) = x[tril(ones(Int, size(x))) .== 1]
Base.range(x::AbstractArray) = [minimum(x), maximum(x)]
# x = reshape(collect(1:9), 3, 3)
# trilindex(x)

K = 6
meantrue = randn(K) # zeros(K)
covtrue = let
  M = randn(K,K)
  M*M'
end

m = DensityModel(x -> logpdf(MvNormal(meantrue, covtrue), x))
# p = RWMH(MvNormal(K, .1))
p = RWAM(K)
c = sample(m, p, Int(60000), chain_type=Chains);
c = c[30001:10:end];

plot(vec(get(c, :lp)[1].data), label=nothing)
meanx = mean(c).nt[:mean]

X = c.value.data[:, 1:end-1, 1]
covtrue
cov(X)
scatter(trilindex(covtrue), trilindex(cov(X)), label=nothing, color=:blue, alpha=.6)
plot!(range(covtrue), range(covtrue), ls=:dot, label=nothing)
corrplot(X[:, 1:3], grid=false)

scatter(meantrue, meanx, label=nothing, color=:blue, alpha=.6)
plot!(range(meantrue), range(meantrue), ls=:dot, label=nothing)

# boxplot(X, label=nothing)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
# Neutron mass [MeV/c²]
const mₙ = 939.565
# Proton mass [MeV/c²]
const mₚ = 939.272
# Reduced np mass [MeV/c²]
const m_np= 2*(mₙ*mₚ)/(mₙ+mₚ) /197
const m_pn = m_np
# Pion± mass [MeV/c²]
const m_πpm = 139.570
# Pion0 mass [MeV/c²]
const m_π0 = 134.977
const m_π = (2m_πpm + m_π0)/3

max_stars_repo_name: ErlendLima/Scattering

max_issues_repo_name: ErlendLima/Scattering

max_forks_repo_name: ErlendLima/Scattering


text:
function RangeRateTroposphericDelayModifier(arg0::DiscreteTroposphericModel, arg1::jboolean)
    return RangeRateTroposphericDelayModifier((DiscreteTroposphericModel, jboolean), arg0, arg1)
end

function get_parameters_drivers(obj::RangeRateTroposphericDelayModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::RangeRateTroposphericDelayModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end

function range_rate_error_tropospheric_model(obj::RangeRateTroposphericDelayModifier, arg0::GroundStation, arg1::FieldSpacecraftState, arg2::Vector{RealFieldElement})
    return jcall(obj, "rangeRateErrorTroposphericModel", RealFieldElement, (GroundStation, FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function range_rate_error_tropospheric_model(obj::RangeRateTroposphericDelayModifier, arg0::GroundStation, arg1::SpacecraftState)
    return jcall(obj, "rangeRateErrorTroposphericModel", jdouble, (GroundStation, SpacecraftState), arg0, arg1)
end