text:
"""
`module  JAC.ManyElectron`
    ... a submodel of JAC that contains all structs and methods for many-electron CSF, basis and levels.
"""
module ManyElectron

    using  ..Basics, ..Defaults,  ..Radial
    export Configuration, ConfigurationR, AsfSettings, CsfR, Basis, Level, Multiplet,
           AbstractQedModel, NoneQed, QedPetersburg, QedSydney, LSjjSettings,
           AbstractConfigurationRestriction, RestrictNoElectronsTo, RestrictParity, RestrictToShellDoubles, RequestMinimumOccupation, RequestMaximumOccupation

    """
    `abstract type ManyElectron.AbstractQedModel`
        ... defines an abstract and a number of singleton types for dealing with QED estimates in the many-electron
            computations.

      + struct QedPetersburg
        ... to estimate the QED corrections due to the model Hamiltonian by Shabaev and coworkers (2013).

      + struct QedSydney
        ... to estimate the QED corrections by means of the radiative potential by Flambaum and Ginges (2004).

      + struct NoneQed
        ... No QED estimates are included into the computations.
    """
    abstract type  AbstractQedModel                      end
    struct         QedPetersburg  <:  AbstractQedModel   end
    struct         QedSydney      <:  AbstractQedModel   end
    struct         NoneQed        <:  AbstractQedModel   end

    """
    `abstract type ManyElectron.AbstractConfigurationRestriction`
        ... defines an abstract types for dealing with restrictions that need to be applied to a list of configurations.
            Typically, a loop through is made through all given restrictions and all configurations are tested to obey all
            these restrictions. Two contradicting restrictions, for instance RestrictParity(plus) & RestrictParity(minus),
            therefore leads zero configurations in all cases. It remains the reponsibility of the user to make sure that
            the given restrictions are consistent with what is to be achieved. The given set of restrictions can be easily
            extended if this need arises by the users.

        + RestrictNoElectronsTo(..)    ... to restrict the total number of electron in high subshells.
        + RestrictParity(..)           ... to restrict to configurations with a given parity.
        + RestrictToShellDoubles(..)   ... to allow only double occupations in high subshells.
        + RequestMinimumOccupation(..) ... to request a minimum occupation in a given set of shells.
        + RequestMaximumOccupation(..) ... to request a maximum occupation in a given set of shells.
    """
    abstract type  AbstractConfigurationRestriction      end

    """
    `struct  ManyElectron.RestrictNoElectronsTo  <: AbstractConfigurationRestriction`
        ... restrict the number of electron in all shells with principal quantum number n >= nmin or orbital angular momentum l >= lmin
            to a total of ne electrons.

        + ne            ::Int64     ... maximum number of (allowed) electrons in the specicied higher subshells.
        + nmin          ::Int64     ... principal quantum number nmin.
        + lmin          ::Int64     ... orbital angular momentum lmin.
    """
    struct   RestrictNoElectronsTo  <: AbstractConfigurationRestriction
        ne              ::Int64
        nmin            ::Int64
        lmin            ::Int64
    end

    function Base.string(res::RestrictNoElectronsTo)
        sa = "Restrict to configurations with a maximum of $(res.ne) electrons in shells with n >= $(res.nmin) & l >= $(res.lmin)."
        return( sa )
    end

    function Base.show(io::IO, res::RestrictNoElectronsTo)
        sa = string(res);       print(io, sa)
    end

    """
    `struct  ManyElectron.RestrictParity  <: AbstractConfigurationRestriction`
        ... restrict to configurations with a given parity.

        + parity        ::Basics.Parity   ... given parity.
    """
    struct   RestrictParity  <: AbstractConfigurationRestriction
        parity          ::Basics.Parity
    end

    function Base.string(res::RestrictParity)
        sa = "Restrict to configurations with parity $(res.parity)."
        return( sa )
    end

    function Base.show(io::IO, res::RestrictParity)
        sa = string(res);       print(io, sa)
    end

    """
    `struct  ManyElectron.RestrictToShellDoubles  <: AbstractConfigurationRestriction`
        ... restrict to a double electron occupation in all shells with principal quantum number n >= nmin or orbital angular momentum l >= lmin.

        + nmin          ::Int64     ... principal quantum number nmin.
        + lmin          ::Int64     ... orbital angular momentum lmin.
    """
    struct   RestrictToShellDoubles  <: AbstractConfigurationRestriction
        nmin            ::Int64
        lmin            ::Int64
    end

    function Base.string(res::RestrictToShellDoubles)
        sa = "Restrict to configurations with a double electron occupation in shells with n >= $(res.nmin) & l >= $(res.lmin)."
        return( sa )
    end

    function Base.show(io::IO, res::RestrictToShellDoubles)
        sa = string(res);       print(io, sa)
    end

    """
    `struct  ManyElectron.RequestMinimumOccupation  <: AbstractConfigurationRestriction`
        ... request a minimum occupation ne in the given (list of) shells.

        + ne            ::Int64           ... minimum electron occupation.
        + shells        ::Array{Shell,1}  ... list of shells.
    """
    struct   RequestMinimumOccupation  <: AbstractConfigurationRestriction
        ne              ::Int64
        shells          ::Array{Shell,1}
    end

    function Base.string(res::RequestMinimumOccupation)
        sa = "Request a minimum occupation of $(res.ne) electron in the (list of) shells $(res.shells)."
        return( sa )
    end

    function Base.show(io::IO, res::RequestMinimumOccupation)
        sa = string(res);       print(io, sa)
    end

    """
    `struct  ManyElectron.RequestMaximumOccupation  <: AbstractConfigurationRestriction`
        ... request a maximum occupation ne in the given (list of) shells.

        + ne            ::Int64           ... maximum electron occupation.
        + shells        ::Array{Shell,1}  ... list of shells.
    """
    struct   RequestMaximumOccupation  <: AbstractConfigurationRestriction
        ne              ::Int64
        shells          ::Array{Shell,1}
    end

    function Base.string(res::RequestMaximumOccupation)
        sa = "Request a maximum occupation of $(res.ne) electron in the (list of) shells $(res.shells)."
        return( sa )
    end

    function Base.show(io::IO, res::RequestMaximumOccupation)
        sa = string(res);       print(io, sa)
    end

    """
    `struct  ManyElectron.Configuration`
        ... defines a struct for a non-relativistic electron configuration that is fully speficied by its shell notations
            (such as "1s", "2s", "2p", ....) and the corresponding occupation numbers (>= 0). An electron configuration
            is independent of any order of the shells, and a zero occupation is assumed for all shells that do not appear as a
            key in the shell (dictionary). Therefore, the number of keys do not allow any conclusion about the underlying orbital
            space of any considered computation that include more than just a single configuration.

        + shells         ::Dict{Shell,Int64}   ... Dictionary that maps shells to their occupation.
        + NoElectrons    ::Int64               ... No. of electrons.
    """
    struct  Configuration
        shells           ::Dict{Shell,Int64}
        NoElectrons      ::Int64
    end

    """
    `ManyElectron.Configuration(sa::String)`
        ... constructor for a given configuration string, such as "[He]", "[Ne]", "[Ne] 3s 3p^6"  or "1s 2p^6 3s^2 3p".
    """
    function Configuration(sa::String)
        sax = strip(sa)
        shellOccList = split(sax, " ")
        ##x println("shellOccList = $shellOccList")
        NoElectrons = 0;    shellDict    = Dict{Shell,Int64}()
        for shocc in shellOccList
            if       shocc == " "    continue
            elseif   shocc  in [ "[He]", "[Ne]", "[Ar]", "[Kr]", "[Xe]" ]
                shellDict = Base.merge( shellDict, Dict( Shell("1s") => 2));                                       NoElectrons = NoElectrons + 2
                shocc == "[He]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("2s") => 2,  Shell("2p") => 6));                     NoElectrons = NoElectrons + 8
                shocc == "[Ne]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("3s") => 2,  Shell("3p") => 6));                     NoElectrons = NoElectrons + 8
                shocc == "[Ar]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("3d") => 10, Shell("4s") => 2, Shell("4p") => 6));   NoElectrons = NoElectrons + 18
                shocc == "[Kr]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("4d") => 10, Shell("5s") => 2, Shell("5p") => 6));   NoElectrons = NoElectrons + 18
                shocc == "[Xe]"  &&  continue
            else
                ia = findall(in("^"), shocc);
                if     length(ia) == 0     shellDict = Base.merge( shellDict, Dict( Shell(shocc) => 1));           NoElectrons = NoElectrons + 1
                else   occ = parse(Int64, shocc[ia[1]+1:end] );    sh = Shell( shocc[1:ia[1]-1] )
                       shellDict = Base.merge( shellDict, Dict( sh => occ));                                       NoElectrons = NoElectrons + occ
                end
            end
        end
        return( Configuration(shellDict, NoElectrons) )
    end

    # `Base.show(io::IO, conf::Configuration)`  ... prepares a proper printout of the non-relativistic configuration conf::Configuration.
    function Base.show(io::IO, conf::Configuration)
        sa = "Configuration: " * string(conf)
        print(io, sa)
    end

    # `Base.unique(confs::Array{Configuration,1})`  ... return a unique list of configurations.
    function Base.unique(confs::Array{Configuration,1})
        confList       = Configuration[]
        configurations = deepcopy(confs)
        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

        return( confList )
    end

    # `Base.string(conf::Configuration)`  ... provides a String notation for the variable conf::Configuration.
    function Base.string(conf::Configuration)
        wa = keys(conf.shells);   va = values(conf.shells)
        ##x wb = Defaults.getDefaults("ordered shell list: non-relativistic", 7)
        shells = Shell[];   for  k in wa    push!(shells, k)    end
        sortedShells = Base.sort( shells , lt=Base.isless)

        sa = ""
        for  k  in sortedShells
            occ = conf.shells[k]
            sa  = sa * string(k) * "^$occ "
        end

        return( sa )
    end

    """
    `Base.:(==)(confa::Configuration, confb::Configuration)`
        ... compares (recursively) two non-relativistic configurations and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(confa::Configuration, confb::Configuration)
        if   confa.NoElectrons  != confb.NoElectrons    return( false )    end
        allShells = merge(confa.shells, confb.shells)
        wk        = keys(allShells)
        ##x @show confa, confb, wk
        for  k in wk
            if  haskey(confa.shells, k)   &&   haskey(confb.shells, k)  &&
                                               confa.shells[k] != confb.shells[k]     return( false )    end
            if  !haskey(confa.shells, k)  &&   confb.shells[k] != 0                   return( false )    end
            if  !haskey(confb.shells, k)  &&   confa.shells[k] != 0                   return( false )    end
        end
        return( true )
    end

    """
    `struct  ManyElectron.ConfigurationR`
        ... defines a type for a relativistic electron configuration that is fully speficied by its subshell notations
            (such as "1s_1/2", "2s_1/2", "2p_3/2", ....) and the corresponding occupation numbers (>= 0). An electron
            configuration is independent of any order of the subshells, and a zero occupation is assumed for all shells that
            do not appears as a key. Therefore, the number of keys do not allow any conclusion about the underlying orbital
            space of any considered computation that include more than a single configuration.

        + subshells      ::Dict{Subshell,Int64}  ... Dictionary that maps subshells to their occupation.
        + NoElectrons    ::Int64                 ... No. of electrons.
    """
    struct  ConfigurationR
        subshells        ::Dict{Subshell,Int64}
        NoElectrons      ::Int64
    end

    # `Base.show(io::IO, conf::ConfigurationR)`  ... prepares a proper printout of the relativistic configuration conf::ConfigurationR.
    function Base.show(io::IO, conf::ConfigurationR)
        sa = "Configuration: " * string(conf)
        print(io, sa)
    end

    # `Base.string(conf::ConfigurationR)`  ... provides a String notation for the variable conf::ConfigurationR.
    function Base.string(conf::ConfigurationR)
        wa = keys(conf.subshells);   va = values(conf.subshells)
        ##x wb = Defaults.getDefaults("ordered subshell list: relativistic", 7)
        subshells = Subshell[];   for  k in wa    push!(subshells, k)    end
        sortedSubshells = Base.sort( subshells , lt=Base.isless)

        sa = ""
        for  k  in sortedSubshells
            occ = conf.subshells[k]
            sa  = sa * string(k) * "^$occ "
        end

        return( sa )
    end

    """
    `Base.:(==)(confa::ConfigurationR, confb::ConfigurationR)`
        ... compares (recursively) two relativistic configurations and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(confa::ConfigurationR, confb::ConfigurationR)
        if   confa.NoElectrons  != confb.NoElectrons    return( false )    end
        wk = keys(confa.subshells)
        for  k in wk
            if  !haskey(confb.subshells, k)  ||   confa.subshells[k] != confb.subshells[k]    return( false )    end
        end
        return( true )
    end

    """
    `abstract type ManyElectron.AbstractStartOrbitals`
        ... defines an abstract and a number of singleton types to determine the orbitals from which the SCF computations start.

        + struct StartFromHydrogenic   ... to start the SCF procedure from hydrogenic orbitals.
        + struct StartFromPrevious     ... to start the SCF procedure from previously calculated orbitals.
    """
    abstract type  AbstractStartOrbitals                          end
    struct     StartFromHydrogenic    <:  AbstractStartOrbitals   end

    export  AbstractStartOrbitals, StartFromHydrogenic, StartFromPrevious

    """
    `struct  ManyElectron.StartFromPrevious  <: ManyElectron.AbstractStartOrbitals`   ... to start from a set of previously calculated orbitals.

        + orbitals      ::Dict{Subshell, Orbital}  ... set of previously-calculated orbitals.
    """
    struct   StartFromPrevious  <: ManyElectron.AbstractStartOrbitals
        orbitals        ::Dict{Subshell, Orbital}
    end

    function Base.string(from::ManyElectron.StartFromPrevious)
        sa = " "
        for  (sh,v) in from.orbitals  sa = sa * "  " * string(sh)   end
        return( sa )
    end

    function Base.show(io::IO, from::ManyElectron.StartFromPrevious)
        sa = "Start from previously-calculated orbitals for shells:\n" * string(from);       print(io, sa, "\n")
    end

    """
    `struct  ManyElectron.LSjjSettings`  ... defines a type for the details and parameters for performing jj-LS expansions.

        + makeIt         ::Bool            ... True, if the jj-LS expansion is to be made and false otherwise.
        + minWeight      ::Float64         ... minimum weight with which a (relativistic) CSF must contribute
                                               to (at least one) level.
        + printWeight    ::Float64         ... minimum weight of a nonrelativistic CSF to be printed out in the
                                               final expansion.
        + levelSelection ::LevelSelection  ... Specifies the selected levels, if any.
    """
    struct LSjjSettings
        makeIt          ::Bool
        minWeight       ::Float64
        printWeight     ::Float64
        levelSelection  ::LevelSelection
    end

    """
    `LSjjSettings(makeIt::Bool)`  ... constructor for the default values of jj-LS transformations.
    """
    function LSjjSettings(makeIt::Bool)
        LSjjSettings(makeIt, 0.05, 0.1, LevelSelection() )
    end

    # `Base.show(io::IO, settings::LSjjSettings)`  ... prepares a proper printout of the variable settings::LSjjSettings.
    function Base.show(io::IO, settings::LSjjSettings)
        println(io, "makeIt:            $(settings.makeIt)  ")
        println(io, "minWeight:         $(settings.minWeight)  ")
        println(io, "printWeight:       $(settings.printWeight)  ")
        println(io, "levelSelection:    $(settings.levelSelection)  ")
    end

    """
    `struct  ManyElectron.AsfSettings`
        ... a struct for defining the settings for the atomic state functions, i.e. the self-consistent-field (SCF)
            and CI computations

        + generateScf          ::Bool                   ... True, if a SCF need to be generated, and false otherwise
                                                            (frozen orbitals).
        + eeInteraction        ::AbstractEeInteraction  ... Specify the e-e interaction to be included into the SCF
                                                            computations.
        + scField              ::AbstractScField        ... Specify the self-consistent field, for instance,
                                                            Basics.ALField(), etc.
        + startScfFrom         ::AbstractStartOrbitals  ... Specify the orbitals to start the SCF computations
        + maxIterationsScf     ::Int64                  ... maximum number of SCF iterations
        + accuracyScf          ::Float64                ... convergence criterion for the SCF field.
        + shellSequenceScf     ::Array{Subshell,1}      ... Sequence of subshells to be optimized.
        + frozenSubshells      ::Array{Subshell,1}      ... Sequence of subshells to be kept frozen.

    	+ eeInteractionCI      ::AbstractEeInteraction  ... Specify the e-e interaction to be included into the
                                                            CI computations.
    	+ qedModel             ::AbstractQedModel       ... model for estimating QED corrections {NoneQed(),
                                                            QedPetersburg(), QedSydney()}.
    	+ jjLS                 ::LSjjSettings           ... settings to control a jj-LS transformation of atomic
                                                            level, if requested.
        + levelSelectionCI     ::LevelSelection         ... Specifies the selected levels, if any.
    """
    struct  AsfSettings
        generateScf            ::Bool
        eeInteraction          ::AbstractEeInteraction
        scField                ::AbstractScField
        startScfFrom           ::AbstractStartOrbitals
        maxIterationsScf       ::Int64
        accuracyScf            ::Float64
        shellSequenceScf       ::Array{Subshell,1}
        frozenSubshells        ::Array{Subshell,1}
        #
    	eeInteractionCI        ::AbstractEeInteraction
    	qedModel               ::AbstractQedModel
    	jjLS                   ::LSjjSettings
    	levelSelectionCI       ::LevelSelection
     end

    """
    `ManyElectron.AsfSettings()`  ... constructor for setting the default values.
    """
    function AsfSettings()
    	AsfSettings(true, CoulombInteraction(), Basics.DFSField(), StartFromHydrogenic(), 24, 1.0e-6, Subshell[], Subshell[],
    	            CoulombInteraction(), NoneQed(), LSjjSettings(false), LevelSelection() )
    end

    """
    `ManyElectron.AsfSettings(settings::AsfSettings;`

                generateScf=..,       eeInteraction=..,       scField=..,            startScfFrom=..,
                maxIterationsScf=..,  accuracyScf=..,         shellSequenceScf=..,   frozenSubshells=..,
                eeInteractionCI=..,   qedModel=..,            jjLS=..,               levelSelection=..,
                printout::Bool=false)
        ... constructor for re-defining a settings::AsfSettings.
    """
    function AsfSettings(settings::AsfSettings;
        generateScf::Union{Nothing,Bool}=nothing,                       eeInteraction::Union{Nothing,AbstractEeInteraction}=nothing,
        scField::Union{Nothing,AbstractScField}=nothing,                startScfFrom::Union{Nothing,AbstractStartOrbitals}=nothing,
        maxIterationsScf::Union{Nothing,Int64}=nothing,                 accuracyScf::Union{Nothing,Float64}=nothing,
        shellSequenceScf::Union{Nothing,Array{Subshell,1}}=nothing,     frozenSubshells::Union{Nothing,Array{Subshell,1}}=nothing,
        eeInteractionCI::Union{Nothing,AbstractEeInteraction}=nothing,  qedModel::Union{Nothing,AbstractQedModel}=nothing,
        jjLS::Union{Nothing,LSjjSettings}=nothing,
        levelSelectionCI::Union{Nothing,LevelSelection}=nothing,        printout::Bool=false)

        if  generateScf         == nothing   generateScfx          = settings.generateScf           else   generateScfx          = generateScf          end
        if  eeInteraction       == nothing   eeInteractionx        = settings.eeInteraction         else   eeInteractionx        = eeInteraction        end
        if  scField             == nothing   scFieldx              = settings.scField               else   scFieldx              = scField              end
        if  startScfFrom        == nothing   startScfFromx         = settings.startScfFrom          else   startScfFromx         = startScfFrom         end
        if  maxIterationsScf    == nothing   maxIterationsScfx     = settings.maxIterationsScf      else   maxIterationsScfx     = maxIterationsScf     end
        if  accuracyScf         == nothing   accuracyScfx          = settings.accuracyScf           else   accuracyScf           = accuracyScf          end
        if  shellSequenceScf    == nothing   shellSequenceScfx     = settings.shellSequenceScf      else   shellSequenceScfx     = shellSequenceScf     end
        if  frozenSubshells     == nothing   frozenSubshellsx      = settings.frozenSubshells       else   frozenSubshellsx      = frozenSubshells      end
        if  eeInteractionCI     == nothing   eeInteractionCIx      = settings.eeInteractionCI       else   eeInteractionCIx      = eeInteractionCI      end
        if  qedModel            == nothing   qedModelx             = settings.qedModel              else   qedModelx             = qedModel             end
        if  jjLS                == nothing   jjLSx                 = settings.jjLS                  else   jjLSx                 = jjLS                 end
        if  levelSelectionCI    == nothing   levelSelectionCIx     = settings.levelSelectionCI      else   levelSelectionCIx     = levelSelectionCI     end

    	AsfSettings(generateScfx, eeInteractionx, scFieldx, startScfFromx, maxIterationsScfx, accuracyScfx,
    	            shellSequenceScfx, frozenSubshellsx, eeInteractionCIx, qedModelx, jjLSx, levelSelectionCIx)
    end

    # `Base.show(io::IO, settings::AsfSettings)`  ... prepares a proper printout of the settings::AsfSettings.
    function Base.show(io::IO, settings::AsfSettings)
    	  println(io, "generateScf:          $(settings.generateScf)  ")
    	  println(io, "eeInteraction:        $(settings.eeInteraction)  ")
    	  println(io, "scField:              $(settings.scField)  ")
    	  println(io, "startScfFrom:         $(settings.startScfFrom)  ")
    	  println(io, "maxIterationsScf:     $(settings.maxIterationsScf)  ")
    	  println(io, "accuracyScf:          $(settings.accuracyScf)  ")
    	  println(io, "shellSequenceScf:     $(settings.shellSequenceScf)  ")
    	  println(io, "frozenSubshells:      $(settings.frozenSubshells)  ")
    	  #
    	  println(io, "eeInteractionCI:      $(settings.eeInteractionCI)  ")
    	  println(io, "qedModel :            $(settings.qedModel)  ")
    	  println(io, "jjLS:                 $(settings.jjLS.makeIt)  ")
    	  println(io, "levelSelectionCI:     $(settings.levelSelectionCI)  ")
    end

    # `Base.string(settings::AsfSettings)`  ... provides a String notation for the variable settings::AsfSettings.
    function Base.string(settings::AsfSettings)
    	  error("Not yet implemented.")
    	  sa = "Asf settings: maximum No. of iterations = $(settings.maxIterationsScf), accuracy = (settings.accuracyScf)"
    	  return( sa )
    end

    """
    `Base.:(==)(seta::AsfSettings, setb::AsfSettings)`
        ... compares (recursively) two instances of AsfSettings and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(seta::AsfSettings, setb::AsfSettings)

    	if  seta.generateScf      !=  setb.generateScf                return( false )    end
    	if  seta.eeInteraction    !=  setb.eeInteraction              return( false )    end
    	if  seta.scField          !=  setb.scField                    return( false )    end
    	if  seta.startScfFrom     !=  setb.startScfFrom               return( false )    end
    	#   startOrbitals
    	if  seta.maxIterationsScf !=  setb.maxIterationsScf           return( false )    end
    	if  seta.accuracyScf      !=  setb.accuracyScf                return( false )    end
    	if  seta.shellSequenceScf !=  setb.shellSequenceScf           return( false )    end
    	if  seta.frozenSubshells  !=  setb.frozenSubshells            return( false )    end
    	if  seta.eeInteractionCI  !=  setb.eeInteractionCI            return( false )    end
    	#   qedModel
    	#   jjLS

    	  return( true )
    end

    """
    `struct  ManyElectron.CsfR`
        ... defines a type for a relativistic configuration state function (CSF) in terms of its orbitals (sequence of orbitals),
            their occupation as well as the angular momenta, seniority and the coupling of the corresponding antisymmetric subshell states.

        + useStandardSubshells ::Bool                 ... Determines whether the subshell list has to be obtained from some outer structure
                                                          given by the subshells field below.
        + J                    ::AngularJ64           ... Total angular momentum J.
        + parity               ::Parity               ... Total parity.
        + occupation           ::Array{Int64,1}       ... occupation of the orbitals with regard to the specified orbital list.
        + seniorityNr          ::Array{Int64,1}       ... list of seniority or Nr values of the antisymmetric subshell states
        + subshellJ            ::Array{AngularJ64,1}  ... list of J-values of the antisymmetric subshell states.
        + subshellX            ::Array{AngularJ64,1}  ... intermediate X-values in the coupling of the antisymmetric subshell states.
        + subshells            ::Array{Subshell,1}    ... Explicitly given subshell list if useStandardSubshells == false.
    """
    struct  CsfR
        useStandardSubshells   ::Bool
        J		               ::AngularJ64
        parity  	           ::Parity
        occupation	           ::Array{Int64,1}
        seniorityNr	           ::Array{Int64,1}
        subshellJ	           ::Array{AngularJ64,1}
        subshellX	           ::Array{AngularJ64,1}
        subshells	           ::Array{Subshell,1}
    end

    """
    `ManyElectron.CsfR(J::AngularJ64, parity::Parity)`
        ... simple constructor for given J and parity, and where standardOrbitals is set to false.
    """
    function CsfR(J::AngularJ64, parity::String)
    	  CsfR(false, J, parity, Int64[], Int64[], AngularJ64[], AngularJ64[], Subshell[] )
    end

    """
    `ManyElectron.CsfR(sa::String)`
        ... simple constructor for given closed-shell occupation (such as "[Ne]"), and where standardOrbitals is set to true.
    """
    function CsfR(sa::String)
    	  wa = subshellsFromClosedShellConfiguration(sa)

    	  occupation = Int64[];	  subshellJ = AngularJ64[];	subshellX  = AngularJ64[]
    	  subshells  = Subshell[];  seniorityNr = Int64[]

        for  sb  in  wa
    	      #x wb = Basics.SubshellQuantumNumbers(string(sb));   kappa = wb[2];   occ = wb[4] + 1
    	      occ = Basics.subshell_2j(sb) + 1
    	      occupation   = push!(occupation, occ)
    	      seniorityNr  = push!(seniorityNr,    0)
    	      subshellJ    = push!(subshellJ, AngularJ64(0) )
    	      subshellX    = push!(subshellX, AngularJ64(0) )
    	      subshells    = push!(subshells, sb )
    	  end
    	  J = AngularJ64(0);    parity = plus

    	  CsfR(false, J, parity, occupation, seniorityNr, subshellJ, subshellX, subshells)
    end

    # `Base.show(io::IO, csf::CsfR)`  ... prepares a proper printout of the variable orbital::CsfR.
    function Base.show(io::IO, csf::CsfR)
    	  print(io, string(csf) )
    end

    # `Base.string(csf::CsfR)`  ... provides a String notation for csf::CsfR.
    function Base.string(csf::CsfR)
    	  sa = "\n   CSF: "
        if    csf.useStandardSubshells

    	      for  i in 1:length(csf.occupation)
    		    sa = sa * Basics.subshellStateString(string(Defaults.GBL_STANDARD_SUBSHELL_LIST[i]), csf.occupation[i], csf.seniorityNr[i],
                                                  csf.subshellJ[i], csf.subshellX[i])
    		    sa = sa * ", "
    	      end
    	      sa = sa * ": J=" * string(csf.J) * string(csf.parity)

            else
    	      for  i in 1:length(csf.subshells)
    		    sa = sa * Basics.subshellStateString(string(csf.subshells[i]), csf.occupation[i], csf.seniorityNr[i],
                                                  csf.subshellJ[i], csf.subshellX[i])
    		    sa = sa * ", "
    	      end
    	      sa = sa * ": J=" * string(csf.J) * string(csf.parity)
    	  end
    	  return( sa )
    end

    """
    `Base.:(==)(csfa::CsfR, csfb::CsfR)`
        ... compares (recursively) two relativistic CSFs and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(csfa::CsfR, csfb::CsfR)
        # Stop the comparison with an error message if the CSF are not defined w.r.t a common subshell list.
        if   !(csfa.useStandardSubshells)   ||   !(csfb.useStandardSubshells)	## ||   length(csfa.occupation) != length(csfb.occupation)
            error("To determine the equivalence of two CSF requires that both are defined on a common subshell list.")
        end

        if  length(csfa.occupation) != length(csfb.occupation)      return( false )    end
        if  csfa.J  !=  csfb.J  ||  csfa.parity  !=  csfb.parity	return( false )    end
    	if  csfa.occupation      !=  csfb.occupation			    return( false )    end
        if  csfa.seniorityNr       !=  csfb.seniorityNr			        return( false )    end
        if  csfa.subshellJ       !=  csfb.subshellJ			        return( false )    end
        if  csfa.subshellX       !=  csfb.subshellX			        return( false )    end

    	  return( true )
    end

    """
    `ManyElectron.CsfRGrasp92(subshells::Array{Subshell,1}, coreSubshells::Array{Subshell,1}, sa::String, sb::String, sc::String)`
        ... to construct a CsfR from 3 Grasp-Strings and the given list of subshells and core subshells.
    """
    function CsfRGrasp92(subshells::Array{Subshell,1}, coreSubshells::Array{Subshell,1}, sa::String, sb::String, sc::String)
    	occupation = Int64[];	 seniorityNr = Int64[];	 subshellJ = AngularJ64[];    subshellX = AngularJ64[];    subshellsx = Subshell[]

    	# Define the occupation and quantum numbers of the closed shells
    	for i in 1:length(coreSubshells)
    	    push!(occupation, Basics.subshell_2j( coreSubshells[i] ) +1);
    	    push!(seniorityNr, 0);     push!(subshellJ, AngularJ64(0));	 push!(subshellX, AngularJ64(0))
    	end

    	scc = sc * "	  ";   scx = " "
    	i  = -1;   ishell = length(coreSubshells)
    	while true
    	    i = i + 1;    sax = sa[9i+1:9i+9];     sbx = sb[9i+1:9i+9];    scx = strip( scc[9i+6:9i+14] )

    	    sh = strip( sax[1:5] );    sh = Basics.subshellGrasp(sh);    occ = parse( sax[7:8] )
    	    search(sbx, ',') > 0    &&    error("stop a: missing seniorityNr; sb = $sb")   # Include seniorityNr more properly if it occurs
    	    subJx = parse( sbx )
    	    if      typeof(subJx) == Void     subJ = AngularJ64(0)
    	    elseif  typeof(subJx) == Int64    subJ = AngularJ64(subJx)
    	    elseif  typeof(subJx) == Expr     subJ = AngularJ64( subJx.args[2]// subJx.args[3] )
    	    else    error("stop b; type = $(typeof(subJx))")
    	    end

    	    if   length(scx) > 0  &&  scx[end:end] in ["+", "-"]    scx = scx[1:end-1]   end
    	    subXx = parse( scx )
    	    if      ishell == 0  &&  typeof(subXx)     == Void                                       subX = subJ
    	    elseif  ishell >= 1  &&  subshellX[ishell] == AngularJ64(0)  && typeof(subXx) == Void    subX = subJ
    	    elseif  ishell >= 1  &&  typeof(subXx)     == Void                                       subX = subshellX[ishell]
    	    elseif  typeof(subXx) == Void                     subX = subshellX[ishell-1]
    	    elseif  typeof(subXx) == Int64                    subX = AngularJ64(subXx)
    	    elseif  typeof(subXx) == Expr                     subX = AngularJ64( subXx.args[2]// subXx.args[3] )
    	    else    error("stop c")
    	    end

    	    # Fill the arrays due to the pre-defined sequence of subshells
    	    while  true
    		ishell = ishell + 1
    		if  subshells[ishell] == sh
    		    wa = ManyElectron.provideSubshellStates(sh, occ);    nu = -1
    		    for  a in wa
    			if   AngularJ64( a.Jsub2//2 ) == subJ	 nu = a.nu;    break   end
    		    end
    		    nu < 0    &&    error("stop d")
    		    push!(occupation, occ);    push!(seniorityNr, nu);	 push!(subshellJ, subJ);    push!(subshellX, subX)
    		    break
    		else
    		    push!(occupation, 0);      push!(seniorityNr, 0);	 push!(subshellJ, AngularJ64(0));    push!(subshellX, subshellX[ishell-1])
    		end
    	    end

    	    ishell > length(subshells)    &&	error("stop e")

    	    # Terminate if all subshells are read
    	    if  9i + 10 >  length(sa)	 break    end
    	end

      while  ishell < length(subshells)
    	    ishell = ishell + 1
          push!(occupation, 0);    push!(seniorityNr, 0);	 push!(subshellJ, AngularJ64(0));    push!(subshellX, subshellX[ishell-1])
      end

    	# Fill the remaining subshells

    	J = subshellX[end];    scx = strip(sc)
    	if	  scx[end:end] == "+"  parity = Parity("+")
    	elseif  scx[end:end] == "-"  parity = Parity("-")
    	else	  error("stop f")
    	end

    	wa = CsfR(true, J, parity, occupation, seniorityNr, subshellJ, subshellX, subshellsx)
    end

    """
    `ManyElectron.CsfRTransformToStandardSubshells(csf::CsfR, subshells::Array{Subshell,1})`
        ... constructor to re-define a CSF with regard to the given list of subshells. It incorporates 'empty' subshell state
            and their coupling (if necessary) but does not support a re-arrangement of the subshell order. This constructor
            is used to 'combine' individual CsfR into a common list (such as in a Basis) and, thus, use StandardSubshells = true
            is set. This method terminates with an error message if the given subshell list requires a re-coupling of the CSF.
    """
    function CsfRTransformToStandardSubshells(csf::CsfR, subshells::Array{Subshell,1})
    	J = csf.J;    parity = csf.parity;
    	occupation = Int64[];	 seniorityNr = Int64[];	  subshellJ = AngularJ64[];	subshellX = AngularJ64[]

    	ish = 0
    	for  subsh in subshells
    	    push!(subshells, subsh)
    	    found = false
    	    for  i in 1:length(csf.subshells)
    		if     subsh == csf.subshells[i]  &&   i >   ish    i = ish;  found = true;   break
    		elseif subsh == csf.subshells[i]  &&   i <=  ish    error("CSF requires recoupling; this is not supported.")
    		end
    	    end
    	    if  found
    		push!(occupation, csf.occupation[ish]);    push!(seniorityNr, csf.seniorityNr[ish])
    		push!(subshellJ,  csf.subshellJ[ish]);     push!(subshellX, csf.subshellX[ish])
    	    else
    		push!(occupation, 0);	 push!(seniorityNr, 0);	 push!(subshellJ,  0)
    		if    ish == 0    push!(subshellX, AngularJ64(0) )
    		else		  push!(subshellX, subshellX[end])
    		end
    	    end
    	end

    	CsfR(true, J, parity, occupation, seniorityNr, subshellJ, subshellX, Subshell[] )
    end

    """
    `struct  ManyElectron.Basis`
        ... defines a type for a relativistic atomic basis, including the (full) specification of the configuration space
            and radial orbitals.

        + isDefined      ::Bool                      ... Determines whether this atomic basis 'points' to some a well-defined basis
                                                         or to an empty instance, otherwise.
        + NoElectrons    ::Int64                     ... No. of electrons.
        + subshells      ::Array{Subshell,1}         ... Explicitly given subshell list for this basis.
        + csfs           ::Array{CsfR,1}             ... List of CSF.
        + coreSubshells  ::Array{Subshell,1}         ... List of (one-electron) core orbitals.
        + orbitals       ::Dict{Subshell, Orbital}   ... Dictionary of (one-electron) orbitals.
    """
    struct  Basis
        isDefined	       ::Bool
    	  NoElectrons	 ::Int64
        subshells	       ::Array{Subshell,1}
        csfs             ::Array{CsfR,1}
    	  coreSubshells	 ::Array{Subshell,1}
    	  orbitals	       ::Dict{Subshell, Orbital}
    end

    """
    `ManyElectron.Basis()`  ... constructor for an 'empty' instance::Basis with isDefined == false.
    """
    function Basis()
    	  Basis(false, 0, Subshell[], CsfR[], Subshell[], Dict{Subshell, Orbital}() )
    end

    """
    `ManyElectron.Basis(basis::Basis, subshells::Array{Subshell,1})`
        ... constructor to re-expand a given basis with regard to an extended subshell list; this expansion does not change
            the CSF basis but makes sure that the occupation and other quantum numbers are properly expanded.
            An error message is issued if the sequence of the subshells does not allow such a simple re-expansion.
            A newBasis::Basis is returned.
    """
    function Basis(basis::Basis, subshells::Array{Subshell,1})
        if basis.subshells == subshells   return( basis )   end
        # Check that the subshells allow such an re-expansion
        oldSubshells = basis.subshells
        newSubshells = Subshell[]
        for  subsh in subshells
            if subsh in oldSubshells   push!(newSubshells, subsh)   end
        end
        if  oldSubshells != newSubshells         error("Improper subshell order: $oldSubshells != $newSubshells")        end
        if  !basis.csfs[1].useStandardSubshells  error("Re-expansion of basis only if useStandardSubshells == true.")    end

        # Define the newCoreSubshells
        nx = length(basis.coreSubshells);   newCoreSubshells = Subshell[];   ny = 0
        for  ny = nx:-1:1
            if basis.coreSubshells[1:ny] == subshells[1:ny]     newCoreSubshells = copy(subshells[1:ny]);   break   end
        end
        ## x@show basis.coreSubshells, newCoreSubshells
        # Define the new CSF basis
        newCsfs = CsfR[]
        for  csf in basis.csfs
            noccupation = Int64[];   nseniorityNr = Int64[];   nsubshellJ = AngularJ64[];   nsubshellX = AngularJ64[]
            for  subsh  in  subshells
                if  subsh in oldSubshells
                    # Determine the index in oldSubshells
                    ii = 0;     for i = 1:length(oldSubshells)   if  subsh == oldSubshells[i]   ii = i; break   end     end
                    push!(noccupation,  csf.occupation[ii])
                    push!(nseniorityNr, csf.seniorityNr[ii])
                    push!(nsubshellJ,   csf.subshellJ[ii])
                    push!(nsubshellX,   csf.subshellX[ii])
                else
                    # Extend quantum numbers properly
                    push!(noccupation,  0)
                    push!(nseniorityNr, 0)
                    push!(nsubshellJ,   AngularJ64(0))
                    push!(nsubshellX,   nsubshellX[end])
                end
            end
            newcsf = CsfR(csf.useStandardSubshells, csf.J, csf.parity, noccupation, nseniorityNr, nsubshellJ, nsubshellX, csf.subshells)
            push!(newCsfs, newcsf)
        end

        Basis(basis.isDefined, basis.NoElectrons, subshells, newCsfs, newCoreSubshells, basis.orbitals)
    end

    """
    `ManyElectron.Basis("from Grasp2013", cslFilename::String, rwfFilename::String)`
        ... to construct an instance::Basis from the Grasp92/Grasp2013 .csl and .rwf files.
    """
    function Basis(sa::String, cslFilename::String, rwfFilename::String)
        if  sa == "from Grasp2013"
           isDefined = true
           #
           println("ManyElectron.Basis-aa: warning ... The standard grid is set; make an interactive request if not appropriate.")
           grid     = Radial.Grid(true)
           basis    = Basics.read("CSF list: Grasp92", cslFilename)
           orbitals = Basics.readOrbitalFileGrasp92(rwfFilename, grid)
        else  error("stop a")
        end

    	  Basis( isDefined, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals)
    end

    """
    `Base.show(io::IO, basis::Basis; full=false)`  ... prepares a proper printout of the variable basis::Basis.
    """
    function Base.show(io::IO, basis::Basis)
        print(io, string(basis), "\n")
    	  if  basis.isDefined
            println(io, "isDefined:          $(basis.isDefined)  ")
            println(io, "NoElectrons:        $(basis.NoElectrons)  ")
            println(io, "subshells:          $(basis.subshells)  ")
            println(io, "csfs:               $(basis.csfs)  ")
            println(io, "coreSubshells:      $(basis.coreSubshells)  ")
            println(io, "orbitals:           $(basis.orbitals)  ")
    	  end
    end

    # `Base.string(basis::Basis)`  ... provides a String notation for variable basis::Basis.
    function Base.string(basis::Basis)
    	if  basis.isDefined
    	    sa = "Atomic basis:  $(length(basis.csfs)) CSF defined for $(basis.NoElectrons) elecrons."
    	else
    	    sa = "Atomic basis:  Undefined."
    	end
    	return(sa)
    end

    """
    `struct  ManyElectron.Level`
        ... defines a type for an atomic level in terms of its quantum number, energy and with regard to a relativistic basis.

        + J            ::AngularJ64       ... Total angular momentum J.
        + M            ::AngularM64       ... Total projection M, only defined if a particular sublevel is referred to.
        + parity       ::Parity           ... Parity of the level.
        + index        ::Int64            ... Index of this level in its original multiplet, or 0.
        + energy       ::Float64          ... energy
        + relativeOcc  ::Float64          ... Relative occupation of this level, if involved in the evolution of a cascade.
        + hasStateRep  ::Bool             ... Determines whether this level 'point' to a physical representation of an state
                                              (i.e. an basis and corresponding mixing coefficient vector), or just to empty instances, otherwise.
        + basis        ::Basis            ... relativistic atomic basis
        + mc           ::Vector{Float64}  ... Vector of mixing coefficients w.r.t basis.
    """
    struct  Level
        J		       ::AngularJ64
        M		       ::AngularM64
        parity  	   ::Parity
        index  	       ::Int64
        energy  	   ::Float64
        relativeOcc    ::Float64
        hasStateRep    ::Bool
        basis	       ::Basis
        mc		       ::Vector{Float64}
    end

    """
    `ManyElectron.Level()`  ... constructor an empty level without useful data.
    """
    function Level()
        Level( AngularJ64(0),  AngularM64(0), Parity("+"), 0, 0., 0., false, Basis(), Vector{Float64}[] )
    end

    """
    `ManyElectron.Level(J::AngularJ64, parity::Parity, energy::Float64, relativeOcc::Float64)`
        ... constructor an atomic level without a state representation: hasStateRep == false.
    """
    function Level(J::AngularJ64, parity::Parity, energy::Float64, relativeOcc::Float64)
        Level(J, J, parity, 0, energy, relativeOcc, false, Basis(), Vector{Float64}[] )
    end

    """
    `ManyElectron.Level(level::Level, subshells::Array{Subshell,1})`
        ... constructor to re-expand a given level with regard to an extended subshell list; this expansion does not change
            the CSF basis nor the mixing coeffients but makes sure that the occupation and other quantum numbers are properly
            expanded. An error message is issued if the sequence of the subshells does not allow such a simple re-expansion.
            A newLevel::Level is returned.
    """
    function Level(level::Level, subshells::Array{Subshell,1})
        basis = Basis(level.basis, subshells)
        Level(level.J, level.M, level.parity, level.index, level.energy, level.relativeOcc, level.hasStateRep, basis, level.mc)
    end

    # `Base.show(io::IO, level::Level)`  ... prepares a proper printout of the variable level::Level.
    function Base.show(io::IO, level::Level)
        println(io, "Level: J = $(level.J), M = $(level.M), parity = $(level.parity), index = $(level.index) ")
        println(io, "energy:         $(level.energy)  ")
        println(io, "relativeOcc:    $(level.relativeOcc)  ")
        println(io, "hasStateRep:    $(level.hasStateRep)  ")
        println(io, "basis:           (level.basis)  ")
        println(io, "mc:             $(level.mc)  ")
    end

    """
    `struct  ManyElectron.Multiplet`
        ... defines a type for an ordered list of atomic levels; has only the default constructor.

    	  + name  	 ::String	         ... A name associated to the multiplet.
    	  + levels	 ::Array{Level,1}   ... A list of levels (pointers).
    """
    struct  Multiplet
        name       ::String
    	  levels  	 ::Array{Level,1}
    end

    """
     `ManyElectron.Multiplet()`  ... constructor for providing an 'empty' instance of this struct.
    """
    function Multiplet()
        Multiplet("", Level[] )
    end

    """
    `ManyElectron.Multiplet("from Ratip2012", cslFilename::String, rwfFilename::String, mixFilename::String)`
        ... to construct an instance::Multiplet from the Grasp92/Ratip2012 .csl, .rwf and .mix files.  Some consistency
            checks are made and the method terminates with an error message if the files don't fit together.

    + `("from Grasp2018", cslFilename::String, rwfFilename::String, mixFilename::String)`
        ... to do the same but for files from Grasp2018.
    """
    function Multiplet(sa::String, cslFilename::String, rwfFilename::String, mixFilename::String)
        if  sa == "from Ratip2012"
           name      = "from "*cslFilename[1:end-4]
           basis     = ManyElectron.Basis("from Grasp2013", cslFilename, rwfFilename)
           multiplet = Basics.readMixFileRelci(mixFilename, basis)
        else  error("Unrecognized keystrings:  $sa")
        end

        return( multiplet )
    end

    # `Base.show(io::IO, multiplet::Multiplet)`  ... prepares a proper printout of the variable multiplet::Multiplet.
    function Base.show(io::IO, multiplet::Multiplet)
        println(io, "name:        $(multiplet.name)  ")
        println(io, "levels:      $(multiplet.levels)  ")
    end

    # Functions/methods that are later added to the module ManyElectron
    function provideSubshellStates                                  end

end # module


text:
const CURLAUTH_ANY                           = -17
const CURLAUTH_ANYSAFE                       = -18
const CURLAUTH_BASIC                         = 1
const CURLAUTH_DIGEST                        = 2
const CURLAUTH_GSSNEGOTIATE                  = 4
const CURLAUTH_NTLM                          = 8
const CURLCLOSEPOLICY_CALLBACK               = 5
const CURLCLOSEPOLICY_LEAST_RECENTLY_USED    = 2
const CURLCLOSEPOLICY_LEAST_TRAFFIC          = 3
const CURLCLOSEPOLICY_OLDEST                 = 1
const CURLCLOSEPOLICY_SLOWEST                = 4
const CURLE_ABORTED_BY_CALLBACK              = 42
const CURLE_BAD_CALLING_ORDER                = 44
const CURLE_BAD_CONTENT_ENCODING             = 61
const CURLE_BAD_FUNCTION_ARGUMENT            = 43
const CURLE_BAD_PASSWORD_ENTERED             = 46
const CURLE_COULDNT_CONNECT                  = 7
const CURLE_COULDNT_RESOLVE_HOST             = 6
const CURLE_COULDNT_RESOLVE_PROXY            = 5
const CURLE_FAILED_INIT                      = 2
const CURLE_FILESIZE_EXCEEDED                = 63
const CURLE_FILE_COULDNT_READ_FILE           = 37
const CURLE_FTP_ACCESS_DENIED                = 9
const CURLE_FTP_BAD_DOWNLOAD_RESUME          = 36
const CURLE_FTP_CANT_GET_HOST                = 15
const CURLE_FTP_CANT_RECONNECT               = 16
const CURLE_FTP_COULDNT_GET_SIZE             = 32
const CURLE_FTP_COULDNT_RETR_FILE            = 19
const CURLE_FTP_COULDNT_SET_ASCII            = 29
const CURLE_FTP_COULDNT_SET_BINARY           = 17
const CURLE_FTP_COULDNT_STOR_FILE            = 25
const CURLE_FTP_COULDNT_USE_REST             = 31
const CURLE_FTP_PORT_FAILED                  = 30
const CURLE_FTP_QUOTE_ERROR                  = 21
const CURLE_FTP_SSL_FAILED                   = 64
const CURLE_FTP_USER_PASSWORD_INCORRECT      = 10
const CURLE_FTP_WEIRD_227_FORMAT             = 14
const CURLE_FTP_WEIRD_PASS_REPLY             = 11
const CURLE_FTP_WEIRD_PASV_REPLY             = 13
const CURLE_FTP_WEIRD_SERVER_REPLY           = 8
const CURLE_FTP_WEIRD_USER_REPLY             = 12
const CURLE_FTP_WRITE_ERROR                  = 20
const CURLE_FUNCTION_NOT_FOUND               = 41
const CURLE_GOT_NOTHING                      = 52
const CURLE_HTTP_NOT_FOUND                   = 22
const CURLE_HTTP_PORT_FAILED                 = 45
const CURLE_HTTP_POST_ERROR                  = 34
const CURLE_HTTP_RANGE_ERROR                 = 33
const CURLE_LDAP_CANNOT_BIND                 = 38
const CURLE_LDAP_INVALID_URL                 = 62
const CURLE_LDAP_SEARCH_FAILED               = 39
const CURLE_LIBRARY_NOT_FOUND                = 40
const CURLE_MALFORMAT_USER                   = 24
const CURLE_OBSOLETE                         = 50
const CURLE_OK                               = 0
const CURLE_OPERATION_TIMEOUTED              = 28
const CURLE_OUT_OF_MEMORY                    = 27
const CURLE_PARTIAL_FILE                     = 18
const CURLE_READ_ERROR                       = 26
const CURLE_RECV_ERROR                       = 56
const CURLE_SEND_ERROR                       = 55
const CURLE_SHARE_IN_USE                     = 57
const CURLE_SSH                              = 79
const CURLE_SSL_CACERT                       = 60
const CURLE_SSL_CERTPROBLEM                  = 58
const CURLE_SSL_CIPHER                       = 59
const CURLE_SSL_CONNECT_ERROR                = 35
const CURLE_SSL_ENGINE_NOTFOUND              = 53
const CURLE_SSL_ENGINE_SETFAILED             = 54
const CURLE_SSL_PEER_CERTIFICATE             = 51
const CURLE_TELNET_OPTION_SYNTAX             = 49
const CURLE_TOO_MANY_REDIRECTS               = 47
const CURLE_UNKNOWN_TELNET_OPTION            = 48
const CURLE_UNSUPPORTED_PROTOCOL             = 1
const CURLE_URL_MALFORMAT                    = 3
const CURLE_URL_MALFORMAT_USER               = 4
const CURLE_WRITE_ERROR                      = 23
const CURLFTPAUTH_DEFAULT                    = 0
const CURLFTPAUTH_SSL                        = 1
const CURLFTPAUTH_TLS                        = 2
const CURLFTPMETHOD_MULTICWD                 = 1
const CURLFTPMETHOD_NOCWD                    = 2
const CURLFTPMETHOD_SINGLECWD                = 3
const CURLFTPSSL_ALL                         = 3
const CURLFTPSSL_CONTROL                     = 2
const CURLFTPSSL_NONE                        = 0
const CURLFTPSSL_TRY                         = 1
const CURLINFO_CERTINFO                      = 4194338
const CURLINFO_CONNECT_TIME                  = 3145733
const CURLINFO_CONTENT_LENGTH_DOWNLOAD       = 3145743
const CURLINFO_CONTENT_LENGTH_UPLOAD         = 3145744
const CURLINFO_CONTENT_TYPE                  = 1048594
const CURLINFO_EFFECTIVE_URL                 = 1048577
const CURLINFO_FILETIME                      = 2097166
const CURLINFO_HEADER_OUT                    = 2
const CURLINFO_HEADER_SIZE                   = 2097163
const CURLINFO_HTTP_CODE                     = 2097154
const CURLINFO_NAMELOOKUP_TIME               = 3145732
const CURLINFO_PRETRANSFER_TIME              = 3145734
const CURLINFO_PRIVATE                       = 1048597
const CURLINFO_REDIRECT_COUNT                = 2097172
const CURLINFO_REDIRECT_TIME                 = 3145747
const CURLINFO_REDIRECT_URL                  = 1048607
const CURLINFO_REQUEST_SIZE                  = 2097164
const CURLINFO_SIZE_DOWNLOAD                 = 3145736
const CURLINFO_SIZE_UPLOAD                   = 3145735
const CURLINFO_SPEED_DOWNLOAD                = 3145737
const CURLINFO_SPEED_UPLOAD                  = 3145738
const CURLINFO_SSL_VERIFYRESULT              = 2097165
const CURLINFO_STARTTRANSFER_TIME            = 3145745
const CURLINFO_TOTAL_TIME                    = 3145731
const CURLMSG_DONE                           = 1
const CURLM_BAD_EASY_HANDLE                  = 2
const CURLM_BAD_HANDLE                       = 1
const CURLM_CALL_MULTI_PERFORM               = -1
const CURLM_INTERNAL_ERROR                   = 4
const CURLM_OK                               = 0
const CURLM_OUT_OF_MEMORY                    = 3
const CURLOPT_AUTOREFERER                    = 58
const CURLOPT_BINARYTRANSFER                 = 19914
const CURLOPT_BUFFERSIZE                     = 98
const CURLOPT_CAINFO                         = 10065
const CURLOPT_CAPATH                         = 10097
const CURLOPT_CERTINFO                       = 172
const CURLOPT_CLOSEPOLICY                    = 72
const CURLOPT_CONNECTTIMEOUT                 = 78
const CURLOPT_CONNECTTIMEOUT_MS              = 156
const CURLOPT_COOKIE                         = 10022
const CURLOPT_COOKIEFILE                     = 10031
const CURLOPT_COOKIEJAR                      = 10082
const CURLOPT_COOKIESESSION                  = 96
const CURLOPT_CRLF                           = 27
const CURLOPT_CUSTOMREQUEST                  = 10036
const CURLOPT_DNS_CACHE_TIMEOUT              = 92
const CURLOPT_DNS_USE_GLOBAL_CACHE           = 91
const CURLOPT_EGDSOCKET                      = 10077
const CURLOPT_ENCODING                       = 10102
const CURLOPT_FAILONERROR                    = 45
const CURLOPT_FILE                           = 10001
const CURLOPT_FILETIME                       = 69
const CURLOPT_FOLLOWLOCATION                 = 52
const CURLOPT_FORBID_REUSE                   = 75
const CURLOPT_FRESH_CONNECT                  = 74
const CURLOPT_FTPAPPEND                      = 50
const CURLOPT_FTPLISTONLY                    = 48
const CURLOPT_FTPPORT                        = 10017
const CURLOPT_FTPSSLAUTH                     = 129
const CURLOPT_FTP_CREATE_MISSING_DIRS        = 110
const CURLOPT_FTP_FILEMETHOD                 = 138
const CURLOPT_FTP_SKIP_PASV_IP               = 137
const CURLOPT_FTP_SSL                        = 119
const CURLOPT_FTP_USE_EPRT                   = 106
const CURLOPT_FTP_USE_EPSV                   = 85
const CURLOPT_HEADER                         = 42
const CURLOPT_HEADERFUNCTION                 = 20079
const CURLOPT_HTTP200ALIASES                 = 10104
const CURLOPT_HTTPAUTH                       = 107
const CURLOPT_HTTPGET                        = 80
const CURLOPT_HTTPHEADER                     = 10023
const CURLOPT_HTTPPROXYTUNNEL                = 61
const CURLOPT_HTTP_VERSION                   = 84
const CURLOPT_INFILE                         = 10009
const CURLOPT_INFILESIZE                     = 14
const CURLOPT_INTERFACE                      = 10062
const CURLOPT_IPRESOLVE                      = 113
const CURLOPT_KEYPASSWD                      = 10026
const CURLOPT_KRB4LEVEL                      = 10063
const CURLOPT_LOW_SPEED_LIMIT                = 19
const CURLOPT_LOW_SPEED_TIME                 = 20
const CURLOPT_MAXCONNECTS                    = 71
const CURLOPT_MAXREDIRS                      = 68
const CURLOPT_MAX_RECV_SPEED_LARGE           = 30146
const CURLOPT_MAX_SEND_SPEED_LARGE           = 30145
const CURLOPT_NETRC                          = 51
const CURLOPT_NOBODY                         = 44
const CURLOPT_NOPROGRESS                     = 43
const CURLOPT_NOSIGNAL                       = 99
const CURLOPT_PORT                           = 3
const CURLOPT_POST                           = 47
const CURLOPT_POSTFIELDS                     = 10015
const CURLOPT_POSTQUOTE                      = 10039
const CURLOPT_POSTREDIR                      = 161
const CURLOPT_PRIVATE                        = 10103
const CURLOPT_PROGRESSFUNCTION               = 20056
const CURLOPT_PROTOCOLS                      = 181
const CURLOPT_PROXY                          = 10004
const CURLOPT_PROXYAUTH                      = 111
const CURLOPT_PROXYPORT                      = 59
const CURLOPT_PROXYTYPE                      = 101
const CURLOPT_PROXYUSERPWD                   = 10006
const CURLOPT_PUT                            = 54
const CURLOPT_QUOTE                          = 10028
const CURLOPT_RANDOM_FILE                    = 10076
const CURLOPT_RANGE                          = 10007
const CURLOPT_READDATA                       = 10009
const CURLOPT_READFUNCTION                   = 20012
const CURLOPT_REDIR_PROTOCOLS                = 182
const CURLOPT_REFERER                        = 10016
const CURLOPT_RESUME_FROM                    = 21
const CURLOPT_RETURNTRANSFER                 = 19913
const CURLOPT_SSH_AUTH_TYPES                 = 151
const CURLOPT_SSH_HOST_PUBLIC_KEY_MD5        = 10162
const CURLOPT_SSH_PRIVATE_KEYFILE            = 10153
const CURLOPT_SSH_PUBLIC_KEYFILE             = 10152
const CURLOPT_SSLCERT                        = 10025
const CURLOPT_SSLCERTPASSWD                  = 10026
const CURLOPT_SSLCERTTYPE                    = 10086
const CURLOPT_SSLENGINE                      = 10089
const CURLOPT_SSLENGINE_DEFAULT              = 90
const CURLOPT_SSLKEY                         = 10087
const CURLOPT_SSLKEYPASSWD                   = 10026
const CURLOPT_SSLKEYTYPE                     = 10088
const CURLOPT_SSLVERSION                     = 32
const CURLOPT_SSL_CIPHER_LIST                = 10083
const CURLOPT_SSL_VERIFYHOST                 = 81
const CURLOPT_SSL_VERIFYPEER                 = 64
const CURLOPT_STDERR                         = 10037
const CURLOPT_TCP_NODELAY                    = 121
const CURLOPT_TIMECONDITION                  = 33
const CURLOPT_TIMEOUT                        = 13
const CURLOPT_TIMEOUT_MS                     = 155
const CURLOPT_TIMEVALUE                      = 34
const CURLOPT_TRANSFERTEXT                   = 53
const CURLOPT_UNRESTRICTED_AUTH              = 105
const CURLOPT_UPLOAD                         = 46
const CURLOPT_URL                            = 10002
const CURLOPT_USERAGENT                      = 10018
const CURLOPT_USERPWD                        = 10005
const CURLOPT_VERBOSE                        = 41
const CURLOPT_WRITEFUNCTION                  = 20011
const CURLOPT_WRITEHEADER                    = 10029
const CURLPROTO_ALL                          = -1
const CURLPROTO_DICT                         = 512
const CURLPROTO_FILE                         = 1024
const CURLPROTO_FTP                          = 4
const CURLPROTO_FTPS                         = 8
const CURLPROTO_HTTP                         = 1
const CURLPROTO_HTTPS                        = 2
const CURLPROTO_LDAP                         = 128
const CURLPROTO_LDAPS                        = 256
const CURLPROTO_SCP                          = 16
const CURLPROTO_SFTP                         = 32
const CURLPROTO_TELNET                       = 64
const CURLPROTO_TFTP                         = 2048
const CURLPROXY_HTTP                         = 0
const CURLPROXY_SOCKS4                       = 4
const CURLPROXY_SOCKS5                       = 5
const CURLSSH_AUTH_DEFAULT                   = -1
const CURLSSH_AUTH_HOST                      = 4
const CURLSSH_AUTH_KEYBOARD                  = 8
const CURLSSH_AUTH_NONE                      = 0
const CURLSSH_AUTH_PASSWORD                  = 2
const CURLSSH_AUTH_PUBLICKEY                 = 1
const CURLVERSION_NOW                        = 3
const CURL_HTTP_VERSION_1_0                  = 1
const CURL_HTTP_VERSION_1_1                  = 2
const CURL_HTTP_VERSION_NONE                 = 0
const CURL_IPRESOLVE_V4                      = 1
const CURL_IPRESOLVE_V6                      = 2
const CURL_IPRESOLVE_WHATEVER                = 0
const CURL_NETRC_IGNORED                     = 0
const CURL_NETRC_OPTIONAL                    = 1
const CURL_NETRC_REQUIRED                    = 2
const CURL_TIMECOND_IFMODSINCE               = 1
const CURL_TIMECOND_IFUNMODSINCE             = 2
const CURL_TIMECOND_LASTMOD                  = 3
const CURL_VERSION_IPV6                      = 1
const CURL_VERSION_KERBEROS4                 = 2
const CURL_VERSION_LIBZ                      = 8
const CURL_VERSION_SSL                       = 4

ext: jl


text:
# Load modules that are used in the CliMA project.
# These are general modules not necessarily specific
# to CliMA
using MPI
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

# Load modules specific to CliMA project
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Vtk

# Prognostic equations: ρ, (ρu), (ρv), (ρw), (ρe_tot), (ρq_tot)
# Even for the dry example shown here, we load the moist thermodynamics module
# and consider the dry equation set to be the same as the moist equations but
# with total specific humidity = 0.
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters: R_d, cp_d, grav, cv_d, MSLP, T_0

# For a three dimensional problem
const _nstate = 6
const _ρ, _U, _V, _W, _E, _QT = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E, QTid = _QT)
const statenames = ("ρ", "U", "V", "W", "E", "QT")

const _nviscstates = 6
const _τ11, _τ22, _τ33, _τ12, _τ13, _τ23 = 1:_nviscstates

const _ngradstates = 3
const _states_for_gradient_transform = (_ρ, _U, _V, _W)

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
  using Random
end

const γ_exact = 7 // 5
const μ_exact = 10
const xmin = 0
const ymin = 0
const zmin = 0
const xmax = 1500
const ymax = 1500
const zmax =  150
const xc   = xmax / 2
const yc   = ymax / 2
const zc   = zmax / 2

# -------------------------------------------------------------------------
# Preflux calculation: This function computes parameters required for the
# DG RHS (but not explicitly solved for as a prognostic variable)
# In the case of the rising_thermal_bubble example: the saturation
# adjusted temperature and pressure are such examples. Since we define
# the equation and its arguments here the user is afforded a lot of freedom
# around its behaviour.
# The preflux function interacts with the following
# Modules: NumericalFluxes.jl
# functions: wavespeed, cns_flux!, bcstate!
# -------------------------------------------------------------------------
@inline function preflux(Q,VF, aux, _...)
  γ::eltype(Q) = γ_exact
  gravity::eltype(Q) = grav
  R_gas::eltype(Q) = R_d
  @inbounds ρ, U, V, W, E, QT = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E], Q[_QT]
  ρinv = 1 / ρ
  x,y,z = aux[_a_x], aux[_a_y], aux[_a_z]
  u, v, w = ρinv * U, ρinv * V, ρinv * W
  e_int = (E - (U^2 + V^2+ W^2)/(2*ρ) - ρ * gravity * y) / ρ
  qt = QT / ρ
  # Establish the current thermodynamic state using the prognostic variables
  TS = PhaseEquil(e_int, qt, ρ)
  T = air_temperature(TS)
  P = air_pressure(TS) # Test with dry atmosphere
  (P, u, v, w, ρinv)
  # Preflux returns pressure, 3 velocity components, and 1/ρ
end

# -------------------------------------------------------------------------
# max eigenvalue
@inline function wavespeed(n, Q, aux, t, P, u, v, w, ρinv)
  γ::eltype(Q) = γ_exact
  @inbounds abs(n[1] * u + n[2] * v + n[3] * w) + sqrt(ρinv * γ * P)
end

# -------------------------------------------------------------------------
# ### Physical Flux (Required)
#md # Here, we define the physical flux function, i.e. the conservative form
#md # of the equations of motion for the prognostic variables ρ, U, V, W, E, QT
#md # $\frac{\partial Q}{\partial t} + \nabla \cdot \boldsymbol{F} = \boldsymbol {S}$
#md # $\boldsymbol{F}$ contains both the viscous and inviscid flux components
#md # and $\boldsymbol{S}$ contains source terms.
#md # Note that the preflux calculation is splatted at the end of the function call
#md # to cns_flux!
# -------------------------------------------------------------------------
cns_flux!(F, Q, VF, aux, t) = cns_flux!(F, Q, VF, aux, t, preflux(Q,VF, aux)...)
@inline function cns_flux!(F, Q, VF, aux, t, P, u, v, w, ρinv)
  @inbounds begin
    ρ, U, V, W, E, QT = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E], Q[_QT]
    # Inviscid contributions
    F[1, _ρ], F[2, _ρ], F[3, _ρ]    = U          , V          , W
    F[1, _U], F[2, _U], F[3, _U]    = u * U  + P , v * U      , w * U
    F[1, _V], F[2, _V], F[3, _V]    = u * V      , v * V + P  , w * V
    F[1, _W], F[2, _W], F[3, _W]    = u * W      , v * W      , w * W + P
    F[1, _E], F[2, _E], F[3, _E]    = u * (E + P), v * (E + P), w * (E + P)
    F[1, _QT], F[2, _QT], F[3, _QT] = u * QT  , v * QT     , w * QT
    # Stress tensor
    τ11, τ22, τ33 = VF[_τ11], VF[_τ22], VF[_τ33]
    τ12 = τ21 = VF[_τ12]
    τ13 = τ31 = VF[_τ13]
    τ23 = τ32 = VF[_τ23]
    # Viscous contributions
    F[1, _U] -= τ11; F[2, _U] -= τ12; F[3, _U] -= τ13
    F[1, _V] -= τ21; F[2, _V] -= τ22; F[3, _V] -= τ23
    F[1, _W] -= τ31; F[2, _W] -= τ32; F[3, _W] -= τ33
    # Energy dissipation
    F[1, _E] -= u * τ11 + v * τ12 + w * τ13
    F[2, _E] -= u * τ21 + v * τ22 + w * τ23
    F[3, _E] -= u * τ31 + v * τ32 + w * τ33
  end
end

# -------------------------------------------------------------------------
#md # Here we define a function to extract the velocity components from the
#md # prognostic equations (i.e. the momentum and density variables). This
#md # function is not required in general, but provides useful functionality
#md # in some cases.
# -------------------------------------------------------------------------
# Compute the velocity from the state
@inline function velocities!(vel, Q, _...)
  @inbounds begin
    # ordering should match states_for_gradient_transform
    ρ, U, V, W = Q[_ρ], Q[_U], Q[_V], Q[_W]
    ρinv = 1 / ρ
    vel[1], vel[2], vel[3] = ρinv * U, ρinv * V, ρinv * W
  end
end

# -------------------------------------------------------------------------
#md ### Auxiliary Function (Not required)
#md # In this example the auxiliary function is used to store the spatial
#md # coordinates. This may also be used to store variables for which gradients
#md # are needed, but are not available through teh prognostic variable
#md # calculations. (An example of this will follow - in the Smagorinsky model,
#md # where a local Richardson number via potential temperature gradient is required)
# -------------------------------------------------------------------------
const _nauxstate = 3
const _a_x, _a_y, _a_z, = 1:_nauxstate
@inline function auxiliary_state_initialization!(aux, x, y, z)
  @inbounds begin
    aux[_a_x] = x
    aux[_a_y] = y
    aux[_a_z] = z
  end
end

# -------------------------------------------------------------------------
#md ### Viscous fluxes.
#md # The viscous flux function compute_stresses computes the components of
#md # the velocity gradient tensor, and the corresponding strain rates to
#md # populate the viscous flux array VF. SijSij is calculated in addition
#md # to facilitate implementation of the constant coefficient Smagorinsky model
#md # (pending)
@inline function compute_stresses!(VF, grad_vel, _...)
  μ::eltype(VF) = μ_exact
  @inbounds begin
    dudx, dudy, dudz = grad_vel[1, 1], grad_vel[2, 1], grad_vel[3, 1]
    dvdx, dvdy, dvdz = grad_vel[1, 2], grad_vel[2, 2], grad_vel[3, 2]
    dwdx, dwdy, dwdz = grad_vel[1, 3], grad_vel[2, 3], grad_vel[3, 3]
    # strains
    ϵ11 = dudx
    ϵ22 = dvdy
    ϵ33 = dwdz
    ϵ12 = (dudy + dvdx) / 2
    ϵ13 = (dudz + dwdx) / 2
    ϵ23 = (dvdz + dwdy) / 2
    # --------------------------------------------
    # SMAGORINSKY COEFFICIENT COMPONENTS
    # --------------------------------------------
    SijSij = (ϵ11 + ϵ22 + ϵ33
              + 2.0 * ϵ12
              + 2.0 * ϵ13
              + 2.0 * ϵ23)
    # mod strain rate ϵ ---------------------------
    # deviatoric stresses
    VF[_τ11] = 2μ * (ϵ11 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ22] = 2μ * (ϵ22 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ33] = 2μ * (ϵ33 - (ϵ11 + ϵ22 + ϵ33) / 3)
    VF[_τ12] = 2μ * ϵ12
    VF[_τ13] = 2μ * ϵ13
    VF[_τ23] = 2μ * ϵ23

  end
end
# -------------------------------------------------------------------------

# -------------------------------------------------------------------------
# generic bc for 2d , 3d
@inline function bcstate!(QP, VFP, auxP, nM, QM, VFM, auxM, bctype, t, PM, uM, vM, wM, ρinvM)
  @inbounds begin
    x, y, z = auxM[_a_x], auxM[_a_y], auxM[_a_z]
    ρM, UM, VM, WM, EM, QTM = QM[_ρ], QM[_U], QM[_V], QM[_W], QM[_E], QM[_QT]
    UnM = nM[1] * UM + nM[2] * VM + nM[3] * WM
    QP[_U] = UM - 2 * nM[1] * UnM
    QP[_V] = VM - 2 * nM[2] * UnM
    QP[_W] = WM - 2 * nM[3] * UnM
    QP[_ρ] = ρM
    QP[_E] = EM
    QP[_QT] = QTM
    VFP .= VFM
    # To calculate PP, uP, vP, wP, ρinvP we use the preflux function
    nothing
    #preflux(QP, auxP, t)
    # Required return from this function is either nothing or preflux with plus state as arguments
  end
end
# -------------------------------------------------------------------------

@inline stresses_boundary_penalty!(VF, _...) = VF.=0

@inline function stresses_penalty!(VF, nM, velM, QM, aM, velP, QP, aP, t)
  @inbounds begin
    n_Δvel = similar(VF, Size(3, 3))
    for j = 1:3, i = 1:3
      n_Δvel[i, j] = nM[i] * (velP[j] - velM[j]) / 2
    end
    compute_stresses!(VF, n_Δvel)
  end
end
# -------------------------------------------------------------------------

# --------------DEFINE SOURCES HERE -------------------------------#
#  TODO: Make sure that the source values are not being over-written
# ------------------------------------------------------------------
"""
The function source! collects all the individual source terms
associated with a given problem. We do not define sources here,
rather we only call those source terms which are necessary based
on the governing equations.
by terms defined elsewhere
"""
@inline function source!(S,Q,aux,t)

  # Initialise the final block source term
  S .= 0

  # Typically these sources are imported from modules
  @inbounds begin
  #source_squircle_sponge!(S,Q,aux,t)
  source_geopot!(S, Q, aux, t)
  #source_radiation!(S,Q,aux,t)
  #source_ls_subsidence!(S,Q,aux,t)
  end
end

"""
Rayleigh sponge function: Linear damping / relaxation to specified
reference values. In the current implementation we relax velocities
at the boundaries to a still atmosphere.
"""
@inline function source_squircle_sponge!(S,Q,aux,t)
  gravity::eltype(Q) = grav
  α = 1.0
  U, V, W = Q[_U], Q[_V], Q[_W]
  x, y, z = aux[_a_x], aux[_a_y], aux[_a_z]
  rp = (x^4 + y^4 + z^4)^(1/4)
  rsponge = 0.85 * xmax # Sponge damper extent
  @inbounds begin
    if rp > rsponge
      S[_U] += -α * sinpi(1/2 * (rp-rsponge)/rsponge) ^ 4 * U
      S[_V] += -α * sinpi(1/2 * (rp-rsponge)/rsponge) ^ 4 * V
      S[_W] += -α * sinpi(1/2 * (rp-rsponge)/rsponge) ^ 4 * W
    elseif rp > 2 * rsponge
      S[_U] -= U
      S[_V] -= V
      S[_W] -= W
    end
  end
end

"""
Geopotential source term. Gravity forcing applied to the vertical
momentum equation
"""
@inline function source_geopot!(S,Q,aux,t)
  gravity::eltype(Q) = grav
  @inbounds begin
    ρ, U, V, W, E  = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]
    S[_V] += - ρ * gravity
  end
end

"""
Large scale subsidence common to several atmospheric observational
campaigns. In the absence of a GCM to drive the flow we may need to
specify a large scale forcing function.
"""
@inline function source_ls_subsidence!(S,Q,aux,t)
  @inbounds begin
    nothing
  end
end

# ------------------------------------------------------------------
# -------------END DEF SOURCES-------------------------------------#

# initial condition
"""
User-specified. Required.
This function specifies the initial conditions for the Rising Thermal
Bubble driver.
"""
function rising_thermal_bubble!(dim, Q, t, x, y, z, _...)
  DFloat                = eltype(Q)
  γ::DFloat             = γ_exact
  # can override default gas constants
  # to moist values later in the driver
  R_gas::DFloat         = R_d
  c_p::DFloat           = cp_d
  c_v::DFloat           = cv_d
  p0::DFloat            = MSLP
  gravity::DFloat       = grav
  # initialise with dry domain
  q_tot::DFloat         = 0
  q_liq::DFloat         = 0
  q_ice::DFloat         = 0
  # perturbation parameters for rising bubble
  r                     = sqrt((x - xc)^2 + (y - yc)^2)
  rc::DFloat            = 300
  θ_ref::DFloat         = 300
  θ_c::DFloat           = 5.0
  Δθ::DFloat            = 0.0
  if r <= rc
    Δθ = θ_c * (1 + cospi(r/rc))/2
  end
  θ                     = θ_ref + Δθ # potential temperature
  π_exner               = 1.0 - gravity / (c_p * θ) * y # exner pressure
  ρ                     = p0 / (R_gas * θ) * (π_exner)^ (c_v / R_gas) # density

  P                     = p0 * (R_gas * (ρ * θ) / p0) ^(c_p/c_v) # pressure (absolute)
  T                     = P / (ρ * R_gas) # temperature
  U, V, W               = 0.0 , 0.0 , 0.0  # momentum components
  # energy definitions
  e_kin                 = (U^2 + V^2 + W^2) / (2*ρ)/ ρ
  e_pot                 = gravity * y
  e_int                 = c_v * (T - T_0) #internal_energy(T, q_tot, q_liq, q_ice)
  E                     = ρ * (e_int + e_kin + e_pot)  #* total_energy(e_kin, e_pot, T, q_tot, q_liq, q_ice)
  @inbounds Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E], Q[_QT]= ρ, U, V, W, E, ρ * q_tot
end

function run(mpicomm, dim, Ne, N, timeend, DFloat, dt)

  ArrayType = Array
  # CuArray option (TODO merge new master)

  brickrange = (range(DFloat(xmin), length=Ne[1]+1, DFloat(xmax)),
                range(DFloat(ymin), length=Ne[2]+1, DFloat(ymax)),
                range(DFloat(zmin), length=Ne[3]+1, DFloat(zmax)))

  # User defined periodicity in the topl assignment
  # brickrange defines the domain extents
  topl = BrickTopology(mpicomm, brickrange, periodicity=(false,false,false))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DFloat,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N)

  numflux!(x...) = NumericalFluxes.rusanov!(x..., cns_flux!, wavespeed, preflux)
  numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(x..., cns_flux!, bcstate!, wavespeed, preflux)

  # spacedisc = data needed for evaluating the right-hand side function
  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = _nstate,
                           flux! = cns_flux!,
                           numerical_flux! = numflux!,
                           numerical_boundary_flux! = numbcflux!,
                           number_gradient_states = _ngradstates,
                           states_for_gradient_transform =
                            _states_for_gradient_transform,
                           number_viscous_states = _nviscstates,
                           gradient_transform! = velocities!,
                           viscous_transform! = compute_stresses!,
                           viscous_penalty! = stresses_penalty!,
                           viscous_boundary_penalty! = stresses_boundary_penalty!,
                           auxiliary_state_length = _nauxstate,
                           auxiliary_state_initialization! =
                           auxiliary_state_initialization!,
                           source! = source!)

  # This is a actual state/function that lives on the grid
  initialcondition(Q, x...) = rising_thermal_bubble!(Val(dim), Q, DFloat(0), x...)
  Q = MPIStateArray(spacedisc, initialcondition)

  lsrk = LowStorageRungeKutta(spacedisc, Q; dt = dt, t0 = 0)

  eng0 = norm(Q)
  @info @sprintf """Starting
  norm(Q₀) = %.16e""" eng0

  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(5, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      energy = norm(Q)
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     energy)
    end
  end

  step = [0]
  mkpath("vtk")
  cbvtk = GenericCallbacks.EveryXSimulationSteps(10) do (init=false)
    outprefix = @sprintf("vtk/cns_%dD_mpirank%04d_step%04d", dim,
                         MPI.Comm_rank(mpicomm), step[1])
    @debug "doing VTK output" outprefix
    writevtk(outprefix, Q, spacedisc, statenames)
    step[1] += 1
    nothing
  end

  # solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, ))
  solve!(Q, lsrk; timeend=timeend, callbacks=(cbinfo, cbvtk))

  # Print some end of the simulation information
  engf = norm(Q)
  if integration_testing
    Qe = MPIStateArray(spacedisc,
                       (Q, x...) -> initialcondition!(Val(dim), Q,
                                                      DFloat(timeend), x...))
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf/eng0 engf-eng0 errf errf / engfe
  else
    @info @sprintf """Finished
    norm(Q)            = %.16e
    norm(Q) / norm(Q₀) = %.16e
    norm(Q) - norm(Q₀) = %.16e""" engf engf/eng0 engf-eng0
  end
  integration_testing ? errf : (engf / eng0)
end

using Test
let
  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())
  mpicomm = MPI.COMM_WORLD
  if MPI.Comm_rank(mpicomm) == 0
    ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
    loglevel = ll == "DEBUG" ? Logging.Debug :
    ll == "WARN"  ? Logging.Warn  :
    ll == "ERROR" ? Logging.Error : Logging.Info
    global_logger(ConsoleLogger(stderr, loglevel))
  else
    global_logger(NullLogger())
  end
    # User defined number of elements
    # User defined timestep estimate
    # User defined simulation end time
    # User defined polynomial order
    numelem = (5,5,1)
    dt = 1e-2
    timeend = 10
    polynomialorder = 5
    for DFloat in (Float64,) #Float32)
      for dim = 3:3
        engf_eng0 = run(mpicomm, dim, numelem[1:dim], polynomialorder, timeend,
                        DFloat, dt)
      end
    end
  end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
@testset "base interface" begin

############################################################################################

@testset "method caching" begin

import InteractiveUtils: _dump_function

# #17057 fallout
@noinline post17057_child(i) = sink(i)
function post17057_parent(arr::Ptr{Int64})
    i = post17057_child(0)
    unsafe_store!(arr, i, i)
end

# bug: default module activation segfaulted on NULL child function if cached=false
params = Base.CodegenParams(cached=false)
if VERSION >= v"1.1.0-DEV.762"
    _dump_function(post17057_parent, Tuple{Ptr{Int64}},
                   #=native=#false, #=wrapper=#false, #=strip=#false,
                   #=dump_module=#true, #=syntax=#:att, #=optimize=#false, :none,
                   params)
else
    _dump_function(post17057_parent, Tuple{Ptr{Int64}},
                   #=native=#false, #=wrapper=#false, #=strip=#false,
                   #=dump_module=#true, #=syntax=#:att, #=optimize=#false,
                   params)
end

end

############################################################################################

end


text:
## The macro is adapated from Suppressor.jl
## https://github.com/JuliaIO/Suppressor.jl/blob/master/src/Suppressor.jl

const stdout_in_R = R"JuliaCall:::stdout_display"

macro capture_out1(block)
    quote
        if ccall(:jl_generating_output, Cint, ()) == 0
            original_stdout = Main.STDOUT
            out_rd, out_wr = redirect_stdout()
            out_reader = @async read(out_rd, String)
        end

        out = try
            Core.eval(Main, Meta.parse($(string(block))))
        finally
            if ccall(:jl_generating_output, Cint, ()) == 0
                redirect_stdout(original_stdout)
                close(out_wr)
            end
        end

        if ccall(:jl_generating_output, Cint, ()) == 0
            stdout_in_R(fetch(out_reader))
        end

        out
    end
end

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
const FWD_SUFFIX = "_fwd"
const BKWD_SUFFIX = "_bkwd"

function split_rxns(metnet::MetNet, iders;
        get_fwd_ider::Function = (rxn) -> string(rxn, FWD_SUFFIX),
        get_bkwd_ider::Function = (rxn) -> string(rxn, BKWD_SUFFIX),
        on_split::Function = (new_model, fwd_idx, bkwd_idx) -> nothing
    )

    M, N = size(metnet)
    idxs = rxnindex.([metnet], iders)
    newM, newN = M, N + length(idxs)

    new_model = expanded_model(metnet, newM, newN)

    for (i, fwd_idx) in enumerate(idxs)

        # backup
        orig_rxn = metnet.rxns[fwd_idx]
        orig_lb = metnet.lb[fwd_idx]

        # transform to forward reaction
        new_model.rxns[fwd_idx] = get_fwd_ider(orig_rxn)
        new_model.lb[fwd_idx] = 0.0

        # add backward reaction
        bkwd_idx = findempty(new_model, :rxns)
        new_model.rxns[bkwd_idx] = get_bkwd_ider(orig_rxn)
        new_model.S[:,bkwd_idx] .= -new_model.S[:,fwd_idx]
        new_model.c[bkwd_idx] = 0.0
        new_model.lb[bkwd_idx] = 0.0
        new_model.ub[bkwd_idx] = abs(orig_lb)

        on_split(new_model, fwd_idx, bkwd_idx)
    end

    return new_model

end

"""
    Returns a new MetNet with no reversible reactions.
"""
function split_revs(metnet::MetNet;
        get_fwd_ider::Function = (rxn) -> string(rxn, FWD_SUFFIX),
        get_bkwd_ider::Function = (rxn) -> string(rxn, BKWD_SUFFIX),
        on_rev::Function = (new_model, fwd_idx, bkwd_idx) -> nothing
    ) # TODO Add tests

    return split_rxns(metnet, revs(metnet);
        get_fwd_ider, get_bkwd_ider, on_split = on_rev
    )
end

# function split_revs(metnet::MetNet;
#         get_fwd_ider::Function = (rxn) -> string(rxn, FWD_SUFFIX),
#         get_bkwd_ider::Function = (rxn) -> string(rxn, BKWD_SUFFIX),
#         on_rev!::Function = (new_model, fwd_idx, bkwd_idx) -> nothing
#     ) # TODO Add tests

#     M, N = size(metnet)
#     revs = revs(metnet)
#     newM, newN = M, N + length(revs)

#     new_model = expanded_model(metnet, newM, newN)

#     # I only need to update the rev reactions
#     _check_and_push!(v::AbstractVector, idx, val) = checkbounds(Bool, v, idx) && push!(v, val)

#     for (i, fwd_idx) in enumerate(findall(revs))

#         # backup
#         orig_rxn = metnet.rxns[fwd_idx]
#         orig_lb = metnet.lb[fwd_idx]

#         # transform to forward reaction
#         new_model.rxns[fwd_idx] = get_fwd_ider(orig_rxn)
#         new_model.lb[fwd_idx] = 0.0

#         # add backward reaction
#         bkwd_idx = findempty(new_model, :rxns)
#         new_model.rxns[bkwd_idx] = get_bkwd_ider(orig_rxn)
#         new_model.S[:,bkwd_idx] .= -new_model.S[:,fwd_idx]
#         new_model.c[bkwd_idx] = 0.0
#         new_model.lb[bkwd_idx] = 0.0
#         new_model.ub[bkwd_idx] = abs(orig_lb)

#         on_rev!(new_model, fwd_idx, bkwd_idx)
#     end

#     return new_model
# end


text:
# ---
# title: 1026. Maximum Difference Between Node and Ancestor
# id: problem1026
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/maximum-difference-between-node-and-ancestor/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, find the maximum value `V` for which there
# exist **different** nodes `A` and `B` where `V = |A.val - B.val|` and `A` is
# an ancestor of `B`.
#
# A node `A` is an ancestor of `B` if either: any child of `A` is equal to `B`,
# or any child of `A` is an ancestor of `B`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/09/tmp-tree.jpg)
#
#
#
#     Input: root = [8,3,10,1,6,null,14,null,null,4,7,13]
#     Output: 7
#     Explanation: We have various ancestor-node differences, some of which are given below :
#     |8 - 3| = 5
#     |3 - 7| = 4
#     |8 - 1| = 7
#     |10 - 13| = 3
#     Among all possible differences, the maximum value of 7 is obtained by |8 - 1| = 7.
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/09/tmp-tree-1.jpg)
#
#
#
#     Input: root = [1,null,2,null,0,3]
#     Output: 3
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[2, 5000]`.
#   * `0 <= Node.val <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=

In this file we simulate a simple rocket landing

Units in kg, m, s

=#
include("..\\src\\rocket\\rocket-setup.jl")
include("..\\src\\dynamics\\trajectory-setup.jl")
include("..\\src\\objective\\LQR_objective.jl")
include("..\\src\\constraints\\constraintManager.jl")
include("..\\src\\auglag\\auglag-core.jl")
include("..\\src\\solver\\AL-Primal-Main-Solver.jl")
include("..\\src\\solver\\QP-AffineEquality-Solver.jl")
include("..\\src\\other_utils\\parsePrimalDual.jl")

include("..\\src\\results\\trajectoryParsing.jl")
include("..\\src\\results\\plotTrajectory.jl")
include("..\\src\\results\\plotConstraintViolation.jl")
include("..\\src\\results\\plotObjective.jl")

# Based on the Falcon 9
# 549,054 kg (Mass)
# 282 s (Specific Impulse)
# select "y" as the vertical
mass = 1
isp = 1
grav = [0; -9.81]
deltaTime = 1
rocket = rocket_simple(mass, isp, grav, deltaTime)

# in m
# The Karman Line (100 km)
rocketStart = [2.0; 20.0; 0.0; 0.0]
rocketEnd = [0.0; 0.0; 0.0; 0.0]#[-5.0; 0.0; 0.0; 0.0]

uHover = mass * grav

# Number of time steps to discretize the trajectory
NSteps = 40
# Initialize the trajectory with a line
initTraj = initializeTraj(rocketStart, rocketEnd, uHover, uHover, NSteps)

# Use a Linear Quadratic Regulator as the cost function
lqrQMat = 10 * Diagonal(I, size(rocketStart, 1))
lqrRMat = 2 * Diagonal(I, Int64(size(rocketStart, 1) / 2))
costFun = makeLQR_TrajReferenced(lqrQMat, lqrRMat, NSteps, initTraj)

ADyn, BDyn = rocketDynamicsFull(rocket, rocketStart, rocketEnd, NSteps)
dynConstraint = AL_AffineEquality(ADyn, BDyn)
lambdaInit = zeros(size(BDyn))

cMRocket = constraintManager([dynConstraint], [lambdaInit])

penaltyStart = 1.0

# Test that the evaluations work
println("\n--------------------------------------------")
println(" Testing constraints are inputted correctly ")
println("--------------------------------------------")
println("Starting constraint violation: ")
println([evalConstraints(cMRocket, initTraj, penaltyStart)])
println("Starting gradient of constraints: ")
println(evalGradConstraints(cMRocket, initTraj, penaltyStart))
println("Starting hessian of constraints: ")
println(evalHessConstraints(cMRocket, initTraj))

# Equiped with the constraint term and the objective term, I now build the
# Augmented Lagrangian
alRocket = augLag(costFun, cMRocket, penaltyStart)

# Test that the evaluations work
println("--------------------------------------------")
println("           Testing AL is Functional         ")
println("--------------------------------------------")
println("Evaluating augmented lagrangian: ")
println([evalAL(alRocket, initTraj)])
println("Evaluating gradient of augmented lagrangian: ")
println(evalGradAL(alRocket, initTraj))
println("Evaluating hessian of augmented lagrangian: ")
println(size(evalHessAl(alRocket, initTraj)))

# Next we select resonable solver parameters
currSolveParams = solverParams(0.1, 0.5,
                                6, 2,
                                10^-4,
                                10, 10^6,
                                2.5, 2, 0.2, 0.2, 0.4)
solParamPrint(currSolveParams)
println()
println()

println("--------------------------------------------")
println("             Beginning Solve                ")
println("--------------------------------------------")

# # Solve the Trajectory Optimization problem
# trajStates, resArr = ALPrimalNewtonMain(initTraj, alRocket, currSolveParams)
#
#
# if true
#     # Get the parsed list of trajectories
#     nDim = size(grav, 1)
#     pltTraj, pltCV, pltObj, plts, pltv, pltu = batchPlot(trajStates, nDim)
# end

costQ = costFun.lqr_qp.QR_Full
costP = - (initTraj' * costQ)'
trajLambdaSolved = solveQP_AffineEq(costQ, costP, ADyn, BDyn)
trajSolved = parsePrimalDualVec(trajLambdaSolved, size(initTraj, 1))

if true
    # Get the parsed list of trajectories
    nDim = size(grav, 1)
    pltTraj, pltCV, pltObj, plts, pltv, pltu = batchPlot(
                                    [initTraj, trajSolved.primals], nDim)
end

if true
    header = "DirectLanding"
    saveBulk(pltTraj, pltCV, pltObj, plts, pltv, pltu, header)
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
import SimpleCache
ENV["JULIA_DEBUG"] = Main;

"Run a Cmd object, returning the stdout & stderr contents plus the exit code"
function execute(cmd::Cmd)
  out = Pipe()
  err = Pipe()

  process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
  close(out.in)
  close(err.in)

  (
    stdout = String(read(out)),
    stderr = String(read(err)),
    code = process.exitcode
  )
end

# Which of the hashes SHOULD change when which value changes is encoded in
# dependency_graph
dependency_graph = Dict(
    "VALUE_A" => [1, 5],
    "VALUE_B" => [2],
    "VALUE_C" => [3],
    "VALUE_D" => [4],
    "VALUE_E" => [5],
)

Base.replace(s::String, oldnews::Pair...) = foldl(replace, oldnews, init=s)

# calc_hashes replaces VALUE_A and VALUE_B in original.jl according to replacedict, then runs this file and report the calculated hashes
function calc_hashes(prefix_newlines::Bool, replacedict)
    src = read("original.jl", String)
    if prefix_newlines
        src = "\n\n" * src
    end
    src = replace(src, replacedict...)
    write("generated.jl", src)
    output = execute(`julia generated.jl`)
    @assert output.code == 0

    ar = split(output.stdout, "\n")
    ar = [el for el in ar if el != ""]
    hashes = parse.(UInt64, ar)
    hashes
end

function gen_default_dict()
    default_dict = Dict()
    counter = 0
    for el in dependency_graph
        counter += 1
        default_dict[el[1]] = counter
    end
    default_dict
end

h1 = calc_hashes(false, gen_default_dict())
h2 = calc_hashes(false, gen_default_dict())
@assert(h1 == h2 , "Instable across restarts!")
h3 = calc_hashes(true, gen_default_dict())
@assert(h1 == h3 , "Instable againes linechanges!")

for el in dependency_graph
    dict = gen_default_dict()
    dict[el[1]] += 100
    hr = calc_hashes(false, dict)
    for i in 1:length(hr)
        if i in el[2]
            if i == 5 # The 5th test (function pointer(x)) is known broken
                continue
            end
            @assert(hr[i] != h1[i], string(i) * "-th Hash stayed identical, but should have changed when " * el[1] * " changed")
        else
            @assert(hr[i] == h1[i], string(i) * "-th Hash changed, but should have stayed identical when " * el[1] * " changed")
        end
    end
end


text:
# # Power flow optimization
# The data for example is taken from [MATPOWER](http://www.pserc.cornell.edu/matpower/) website. MATPOWER is Matlab package for solving power flow and optimal power flow problems.

using Convex, SCS
using Test
using MAT   #Pkg.add("MAT")
aux(str) = joinpath(@__DIR__, "aux_files", str) # path to auxiliary files

TOL = 1e-2;
input = matopen(aux("Data.mat"))
varnames = names(input)
Data = read(input, "inj", "Y");

n = size(Data[2], 1);
Y = Data[2];
inj = Data[1];
W = ComplexVariable(n, n);
objective = real(sum(diag(W)));
c1 = Constraint[];
for i = 2:n
    push!(c1, sum(W[i,:] .* (Y[i,:]')) == inj[i]);
end
c2 = W in :SDP
c3 = real(W[1,1]) == 1.06^2;
push!(c1, c2)
push!(c1, c3)
p = maximize(objective, c1);
solve!(p, () -> SCS.Optimizer(verbose = 0))
p.optval
#15.125857662600703
evaluate(objective)
#15.1258578588357

output = matopen("Res.mat")
names(output)
outputData = read(output, "Wres");
Wres = outputData
real_diff = real(evaluate(W)) - real(Wres);
imag_diff = imag(evaluate(W)) - imag(Wres);
@test real_diff ≈ zeros(n, n) atol = TOL
@test imag_diff ≈ zeros(n, n) atol = TOL

real_diff = real(evaluate(W)) - (real(evaluate(W)))';
imag_sum = imag(evaluate(W)) + (imag(evaluate(W)))';
@test real_diff ≈ zeros(n, n) atol = TOL
@test imag_diff ≈ zeros(n, n) atol = TOL


text:
# Implement code block 2.8 in McElreath's Statistical Rethinking, 2nd edition.

using Random
using Distributions
using Plots

Random.seed!(17)

num_samples = 2000

probas = [0.5]

successes = 6
trials = 9

for i = 2:num_samples

    proba_new = rand(Normal(probas[end], 0.1))

    if proba_new < 0; proba_new = abs(proba_new); end
    if proba_new > 1; proba_new = 2 - proba_new; end

    q0 = pdf(Binomial(trials, probas[end]), successes)
    q1 = pdf(Binomial(trials, proba_new), successes)

    push!(probas, rand(Uniform(0, 1)) < q1 / q0 ? proba_new : probas[end])
end

# println(probas)

display(histogram(probas, bins=:40))


text:
module TestLoading

using Test
using MLJModels

function isloaded(name::String, pkg::String)
    (name, pkg) in map(localmodels()) do m
        (m.name, m.package_name)
    end
end

@load AdaBoostStumpClassifier pkg=DecisionTree verbosity=0
@test isloaded("AdaBoostStumpClassifier", "DecisionTree")

# built-ins load fine:
@load Standardizer verbosity=0

# load one version of a RidgeRegressor:
@test !isloaded("RidgeRegressor", "MultivariateStats")
@load RidgeRegressor pkg=MultivariateStats verbosity=0
@test isloaded("RidgeRegressor", "MultivariateStats")

# error if ambiguous:
@test_throws ArgumentError @load RidgeRegressor

# error if not in project:
@test !isloaded("KMeans", "Clustering")
@test_throws ArgumentError @load KMeans pkg=Clustering verbosity=0

# use add option:
@load KMeans pkg=Clustering verbosity=0 add=true
@test isloaded("KMeans", "Clustering")

# deprecated methods:
@test_throws Exception load("model", pkg = "pkg")
@test_throws Exception load(models()[1])

module FooBar
using MLJModels
function regressor()
    Regressor = @load LinearRegressor pkg=MultivariateStats verbosity=0
    return Regressor()
end
end
using .FooBar

@testset "@load from within a function within a module" begin
    model = FooBar.regressor()
    @test isdefined(model, :bias)
end

end # module

true