text:
"""
`module  JAC.Hfs`
    ... a submodel of JAC that contains all methods for computing HFS A and B coefficients, hyperfine sublevel representations, etc.
"""
module Hfs

    using Printf, ..AngularMomentum, ..Basics,  ..Defaults, ..InteractionStrength, ..ManyElectron, ..Radial, ..Nuclear, ..TableStrings

    """
    `struct  Hfs.InteractionMatrix`  ... defines a type for storing the T^1 and T^2 interaction matrices for a given basis.

        + calcT1   ::Bool               ... true, if the matrixT1 has been calculated and false otherwise.
        + calcT2   ::Bool               ... true, if the matrixT2 has been calculated and false otherwise.
        + matrixT1 ::Array{Float64,2}   ... T1 interaction matrix
        + matrixT2 ::Array{Float64,2}   ... T2 interaction matrix

    """
    struct InteractionMatrix
        calcT1     ::Bool
        calcT2     ::Bool
        matrixT1   ::Array{Float64,2}
        matrixT2   ::Array
    end

    """
    `Hfs.InteractionMatrix()`  ... constructor for an `empty` instance of InteractionMatrix.
    """
    function InteractionMatrix()
        InteractionMatrix(false, false, zeros(2,2), zeros(2,2))
    end

    # `Base.show(io::IO, im::Hfs.InteractionMatrix)`  ... prepares a proper printout of the variable InteractionMatrix.
    function Base.show(io::IO, im::Hfs.InteractionMatrix)
        println(io, "calcT1:           $(im.calcT1)  ")
        println(io, "calcT1:           $(im.calcT1)  ")
        println(io, "matrixT1:         $(im.matrixT1)  ")
        println(io, "matrixT2:         $(im.matrixT2)  ")
    end

    """
    `struct  Hfs.IJF_Vector`  ... defines a type for a IJF-coupled basis vector, here based on an ASF.

        + I        ::AngularJ64   ... Nuclear spin I.
        + F        ::AngularJ64   ... Total angular momentum F
        + levelJ   ::Level        ... Atomic level with well-defined total (electronic) angular momentum J

    """
    struct IJF_Vector
        I          ::AngularJ64
        F          ::AngularJ64
        levelJ     ::Level
    end

    """
    `Hfs.IJF_Vector()`  ... constructor for an `empty` instance of IJF_Vector.
    """
    function IJF_Vector()
        IJF_Vector(AngularJ64(0), AngularJ64(0), Level())
    end

    # `Base.show(io::IO, ijfVector::Hfs.IJF_Vector)`  ... prepares a proper printout of the variable ijfVector.
    function Base.show(io::IO, ijfVector::Hfs.IJF_Vector)
        println(io, "I:           $(ijfVector.I)  ")
        println(io, "F:           $(ijfVector.F)  ")
        println(io, "levelJ:      $(ijfVector.levelJ)  ")
    end

   """
    `struct  Hfs.IJF_Basis`  ... defines a type for a IJF-coupled basis.

        + vectors   ::Array{IJF_Vector,1}   ... List of IJF_Vectors that form the basis.
        + basisJ    ::Basis                 ... Electronic basis that allows access to the electronic orbitals and CSF basis.
    """
    struct IJF_Basis
        vectors     ::Array{IJF_Vector,1}
        basisJ      ::Basis
    end

    """
   `Hfs.IJF_Basis()`  ... constructor for an `empty` instance of IJF_Basis.
    """
    function IJF_Basis()
        IJF_Basis(IJF_Vector[], Basis())
    end

    # `Base.show(io::IO, ijfBasis::Hfs.IJF_Basis)`  ... prepares a proper printout of the variable ijfBasis.
    function Base.show(io::IO, ijfBasis::Hfs.IJF_Basis)
        println(io, "levelFs:           $(ijfBasis.levelFs)  ")
        println(io, "basisJ:            $(ijfBasis.basisJ)  ")
    end

    """
    `struct  Hfs.IJF_Level`  ... defines a type for a IJF-coupled level.

        + I              ::AngularJ64       ... Nuclear spin I.
        + F              ::AngularJ64       ... Total angular momentum F.
        + M              ::AngularM64       ... Total projection M, only defined if a particular sublevel is referred to.
        + parity         ::Parity           ... Parity of the level which corresponds to the electronic system.
        + energy         ::Float64          ... energy
        + hasStateRep    ::Bool             ... Determines whether this IJF-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          ::IJF_Basis        ... basis for this level
        + mc             ::Vector{Float64}  ... Vector of mixing coefficients w.r.t basis.
    """
    struct IJF_Level
        I                ::AngularJ64
        F                ::AngularJ64
        M                ::AngularM64
        parity           ::Parity
        energy           ::Float64
        hasStateRep      ::Bool
        basis            ::IJF_Basis
        mc               ::Vector{Float64}
    end

    """
    `Hfs.IJF_Level()`  ... constructor for an `empty` instance of IJF_Level.
    """
    function IJF_Level()
        IJF_Level(AngularJ64(0), AngularJ64(0), AngularJ64(0), AngularM64(0), Basics.plus, 0., false, IJF_Basis(), Float64[])
    end

    # `Base.show(io::IO, ijfLevel::Hfs.IJF_Level)`  ... prepares a proper printout of the variable ijfLevel::Hfs.IJF_Level.
    function Base.show(io::IO, ijfLevel::Hfs.IJF_Level)
        println(io, "I:              $(ijfLevel.I)  ")
        println(io, "J:              $(ijfLevel.J)  ")
        println(io, "F:              $(ijfLevel.F)  ")
        println(io, "M:              $(ijfLevel.M)  ")
        println(io, "parity:         $(ijfLevel.parity)  ")
        println(io, "energy:         $(ijfLevel.energy)  ")
        println(io, "hasStateRep:    $(ijfLevel.hasStateRep)  ")
        println(io, "basis:          $(ijfLevel.basis)  ")
        println(io, "mc:             $(ijfLevel.mc)  ")
    end

    """
    `struct  Hfs.IJF_Multiplet`  ... defines a type for a multiplet of IJF-coupled levels.

        + name     ::String                ... A name associated to the multiplet.
        + levelFs  ::Array{IJF_Level,1}    ... List of IJF-coupled levels.

    """
    struct IJF_Multiplet
        name       ::String
        levelFs    ::Array{IJF_Level,1}
    end

    """
    `Hfs.IJF_Multiplet()`  ... constructor for an `empty` instance of Hfs.IJF_Multiplet.
    """
    function IJF_Multiplet()
        IJF_Multiplet("", IJF_Level[])
    end

    # `Base.show(io::IO, ijfMultiplet::Hfs.IJF_Multiplet)`  ... prepares a proper printout of the variable ijfMultiplet::Hfs.IJF_Multiplet.
    function Base.show(io::IO, ijfMultiplet::Hfs.IJF_Multiplet)
        println(io, "name:           $(ijfMultiplet.name)  ")
        println(io, "levelFs:        $(ijfMultiplet.levelFs)  ")
    end

    """
    `struct  Hfs.Outcome`
        ... defines a type to keep the outcome of a HFS computation, such as the HFS A and B coefficients as well
            other results.

        + Jlevel                    ::Level            ... Atomic level to which the outcome refers to.
        + AIoverMu                  ::Float64          ... HFS A * I / mu value.
        + BoverQ                    ::Float64          ... HFS B / Q value
        + amplitudeT1               ::Complex{Float64} ... T1 amplitude
        + amplitudeT2               ::Complex{Float64} ... T2 amplitude
        + nuclearI                  ::AngularJ64       ... nuclear spin
        + hasFlevels                ::Bool             ... True if the F-values and F-energies are specified explicitly.
        + Flevels                   ::Bool             ... Array to keep F-values and F-energies.
    """
    struct Outcome
        Jlevel                      ::Level
        AIoverMu                    ::Float64
        BoverQ                      ::Float64
        amplitudeT1                 ::Complex{Float64}
        amplitudeT2                 ::Complex{Float64}
        nuclearI                    ::AngularJ64
        hasFlevels                  ::Bool
        Flevels                     ::Bool
    end

    """
    `Hfs.Outcome()`  ... constructor for an `empty` instance of Hfs.Outcome for the computation of HFS properties.
    """
    function Outcome()
        Outcome(Level(), 0., 0., 0., 0., AngularJ64(0), false, false)
    end

    # `Base.show(io::IO, outcome::Hfs.Outcome)`  ... prepares a proper printout of the variable outcome::Hfs.Outcome.
    function Base.show(io::IO, outcome::Hfs.Outcome)
        println(io, "Jlevel:                    $(outcome.Jlevel)  ")
        println(io, "AIoverMu:                  $(outcome.AIoverMu)  ")
        println(io, "BoverQ:                    $(outcome.BoverQ)  ")
        println(io, "amplitudeT1:               $(outcome.amplitudeT1)  ")
        println(io, "amplitudeT2:               $(outcome.amplitudeT2)  ")
        println(io, "nuclearI:                  $(outcome.nuclearI)  ")
        println(io, "hasFlevels:                $(outcome.hasFlevels)  ")
        println(io, "Flevels:                   $(outcome.Flevels)  ")
    end

    """
    `struct  Settings`  ... defines a type for the details and parameters of computing HFS A and B coefficients.

        + calcT1                    ::Bool             ... True if T1-amplitudes (HFS A values) need to be calculated, and false otherwise.
        + calcT2                    ::Bool             ... True if T2-amplitudes (HFS B values) need to be calculated, and false otherwise.
        + calcNondiagonal           ::Bool             ... True if also (non-)diagonal hyperfine amplitudes are to be calculated and
                                                           printed, and false otherwise.
        + calcIJFexpansion          ::Bool             ... True if the selected atomic levels are to be represented in a IJF-coupled basis,
                                                           and false otherwise.
        + printBefore    ::Bool                        ... True if a list of selected levels is printed before the actual computations start.
        + printDeltaEF              ::Bool             ... True if the energy shift of E_F (with regard to E_J) is to be printed, and false otherwise.
        + levelSelection            ::LevelSelection   ... Specifies the selected levels, if any.
    """
    struct Settings
        calcT1                      ::Bool
        calcT2                      ::Bool
        calcNondiagonal             ::Bool
        calcIJFexpansion            ::Bool
        printBefore                 ::Bool
        printDeltaEF                ::Bool
        levelSelection              ::LevelSelection
    end

    """
    `Hfs.Settings(; calcT1::Bool=true,` calcT2::Bool=false, calcNondiagonal::Bool=false, calcIJFexpansion::Bool=false,
                        printBefore::Bool=false, printDeltaEF::Bool=false,
                        levelSelection::LevelSelection=LevelSelection())
        ... keyword constructor to overwrite selected value of Einstein line computations.
    """
    function Settings(; calcT1::Bool=true, calcT2::Bool=false, calcNondiagonal::Bool=false, calcIJFexpansion::Bool=false,
                        printBefore::Bool=false, printDeltaEF::Bool=false,
                        levelSelection::LevelSelection=LevelSelection())
        Settings(calcT1, calcT2, calcNondiagonal, calcIJFexpansion, printBefore, printDeltaEF, levelSelection)
    end

    # `Base.show(io::IO, settings::Hfs.Settings)`  ... prepares a proper printout of the variable settings::Hfs.Settings.
    function Base.show(io::IO, settings::Hfs.Settings)
        println(io, "calcT1:                   $(settings.calcT1)  ")
        println(io, "calcT2:                   $(settings.calcT2)  ")
        println(io, "calcNondiagonal:          $(settings.calcNondiagonal)  ")
        println(io, "calcIJFexpansion:         $(settings.calcIJFexpansion)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "printDeltaEF:             $(settings.printDeltaEF)  ")
        println(io, "levelSelection:           $(settings.levelSelection)  ")
    end

    """
    `Hfs.amplitude(kind::String, rLevel::Level, sLevel::Level, grid::Radial.Grid; printout::Bool=true)`
        ... to compute either the  T^(1) or T^(2) hyperfine amplitude <alpha_r J_r || T^(n)) || alpha_s J_s>
            for a given pair of levels. A value::ComplexF64 is returned.
    """
    function  amplitude(kind::String, rLevel::Level, sLevel::Level, grid::Radial.Grid; printout::Bool=true)
        #
        if  rLevel.parity != sLevel.parity   return( ComplexF64(0.) )   end
        nr = length(rLevel.basis.csfs);    ns = length(sLevel.basis.csfs);    matrix = zeros(ComplexF64, nr, ns)
        if  printout   printstyled("Compute hyperfine $(kind[1:5]) matrix of dimension $nr x $ns ... \n", color=:light_green)   end
        #
        for  r = 1:nr
            for  s = 1:ns
                me = 0.
                if  rLevel.basis.csfs[r].parity  != rLevel.parity    ||  sLevel.basis.csfs[s].parity  != sLevel.parity  ||
                    rLevel.parity != sLevel.parity    continue
                end
                #
                if      kind == "T^(1) amplitude"
                #--------------------------------
                    wa = compute("angular coefficients: 1-p, Grasp92", 0, 1, rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.hfs_t1(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], grid)
                        #
                        ## println("**  <$(coeff.a) || t1 || $(coeff.b)>  = $(coeff.T * tamp)   = $(coeff.T) * $tamp" )
                        me = me + coeff.T * tamp
                    end
                #
                elseif  kind == "T^(2) amplitude"
                #--------------------------------
                    wa = compute("angular coefficients: 1-p, Grasp92", 0, 2, rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.hfs_t2(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], grid)
                        #
                        ## println("**  <$(coeff.a) || t2 || $(coeff.b)>  = $(coeff.T * tamp)   = $(coeff.T) * $tamp" )
                        me = me + coeff.T * tamp
                    end
                #
                else    error("stop a")
                end
                #
                matrix[r,s] = me
                println("**                                                                    " *
                        "                      <$r || $kind || $s> = $me" )
            end
        end
        if  printout   printstyled("done.\n", color=:light_green)   end
        amplitude = transpose(rLevel.mc) * matrix * sLevel.mc
        #
        return( amplitude )
    end

    """
    `Hfs.computeAmplitudesProperties(outcome::Hfs.Outcome, grid::Radial.Grid, settings::Hfs.Settings, im::Hfs.InteractionMatrix)
        ... to compute all amplitudes and properties of for a given level; an outcome::Hfs.Outcome is returned for which the
            amplitudes and properties are now evaluated explicitly.
    """
    function  computeAmplitudesProperties(outcome::Hfs.Outcome, grid::Radial.Grid, settings::Hfs.Settings, im::Hfs.InteractionMatrix)
        AIoverMu = 0.0;   BoverQ = 0.0;   amplitudeT1 = 0.0;   amplitudeT2 = 0.0;    J = AngularMomentum.oneJ(outcome.Jlevel.J)
        if  settings.calcT1  &&  outcome.Jlevel.J != AngularJ64(0)
            if  im.calcT1   amplitudeT1 = transpose(outcome.Jlevel.mc) * im.matrixT1 * outcome.Jlevel.mc
            else            amplitudeT1 = Hfs.amplitude("T^(1) amplitude", outcome.Jlevel, outcome.Jlevel, grid)
            end

            AIoverMu = amplitudeT1 / sqrt(J * (J+1))
        end
        #
        if  settings.calcT2  &&  outcome.Jlevel.J != AngularJ64(0)
            if  im.calcT2   amplitudeT2 = transpose(outcome.Jlevel.mc) * im.matrixT2 * outcome.Jlevel.mc
            else            amplitudeT2 = Hfs.amplitude("T^(2) amplitude", outcome.Jlevel, outcome.Jlevel, grid)
            end

            BoverQ   = 2 * amplitudeT2 * sqrt( (2J-1) / ((J+1)*(2J+3)) )  ## * sqrt(J)
        end
        newOutcome = Hfs.Outcome( outcome.Jlevel, AIoverMu, BoverQ, amplitudeT1, amplitudeT2,
                                  outcome.nuclearI, outcome.hasFlevels, outcome.Flevels)
        return( newOutcome )
    end

    """
    `Hfs.computeHyperfineMultiplet(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::Hfs.Settings; output=true)`
        ... to compute a hyperfine multiplet, i.e. a representation of hyperfine levels within a hyperfine-coupled basis as defined by the
            given (electronic) multiplet; a hfsMultiplet::IJF_Multiplet is returned.
    """
    function computeHyperfineMultiplet(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::Hfs.Settings; output=true)
        println("")
        printstyled("Hfs.computeHyperfineMultiplet(): The computation of the hyperfine multiplet starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------ \n", color=:light_green)
        #
        hfsBasis     = Hfs.defineHyperfineBasis(multiplet, nm, settings)
        hfsMultiplet = Hfs.computeHyperfineRepresentation(hfsBasis, nm, grid, settings)
        # Print all results to screen
        hfsMultiplet = Basics.sortByEnergy(hfsMultiplet)
        Basics.tabulate("multiplet: energies",                                      hfsMultiplet)
        Basics.tabulate("multiplet: energy of each level relative to lowest level", hfsMultiplet)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary
            Basics.tabulate("multiplet: energies",                                      hfsMultiplet, stream=iostream)
            Basics.tabulate("multiplet: energy of each level relative to lowest level", hfsMultiplet, stream=iostream)
        end
        #
        if    output    return( hfsMultiplet )
        else            return( nothing )
        end
    end

    """
    `Hfs.computeHyperfineRepresentation(hfsBasis::IJF_Basis, nm::Nuclear.Model, grid::Radial.Grid, settings::Hfs.Settings)`
        ... to set-up and diagonalized the Hamiltonian matrix of H^(DFB) + H^(hfs) within the atomic hyperfine (IJF-coupled) basis;
            a hfsMultiplet::IJF_Multiplet is returned.
    """
    function computeHyperfineRepresentation(hfsBasis::IJF_Basis, nm::Nuclear.Model, grid::Radial.Grid, settings::Hfs.Settings)
        n = length(hfsBasis.vectors);   matrix = zeros(n,n)
        for  r = 1:n
            for  s = 1:n
                matrix[r,s] = 0.
                if  r == s    matrix[r,s] = matrix[r,s] + hfsBasis.vectors[r].levelJ.energy                  end
                if  hfsBasis.vectors[r].F              !=  hfsBasis.vectors[s].F                continue     end
                if  hfsBasis.vectors[r].levelJ.parity  !=  hfsBasis.vectors[s].levelJ.parity    continue     end
                # Now add the hyperfine matrix element for K = 1 (magnetic) and K = 2 (electric) contributions
                spinI = nm.spinI;   Jr = hfsBasis.vectors[r].levelJ.J;   Js = hfsBasis.vectors[s].levelJ.J;   F = hfsBasis.vectors[r].F
                Ix    = AngularMomentum.oneJ(spinI)
                #
                wa = AngularMomentum.phaseFactor([spinI, +1, Jr, +1, F])  *
                     AngularMomentum.Wigner_6j(spinI, Jr, F, Js, spinI, AngularJ64(1))
                wb = Hfs.amplitude("T^(1) amplitude", hfsBasis.vectors[r].levelJ, hfsBasis.vectors[s].levelJ, grid::Radial.Grid)
                wc = nm.mu * sqrt( (Ix+1)/Ix )
                matrix[r,s] = matrix[r,s] + wa * wb * wc     * 1.0e-6 ## fudge-factor to keep HFS interaction small
                #
                if  spinI  in [AngularJ64(0), AngularJ64(1//2)]                                 continue     end
                wa = AngularMomentum.phaseFactor([spinI, +1, Jr, +1, F]) *
                     AngularMomentum.Wigner_6j(spinI, Jr, F, Js, spinI, AngularJ64(2))
                wb = Hfs.amplitude("T^(2) amplitude", hfsBasis.vectors[r].levelJ, hfsBasis.vectors[s].levelJ, grid::Radial.Grid)
                wc = nm.Q / 2. * sqrt( (Ix+1)*(2Ix+3)/ (Ix*(2Ix-1)) )
                matrix[r,s] = matrix[r,s] + wa * wb * wc     * 1.0e-6 ## fudge-factor to keep HFS interaction small
            end
        end
        #
        # Diagonalize the matrix and set-up the representation
        eigen    = Basics.diagonalize("matrix: LinearAlgebra", matrix)
        levelFs  = Hfs.IJF_Level[]
        for  ev = 1:length(eigen.values)
            # Construct the eigenvector with regard to the given basis (not w.r.t the symmetry block)
            evector   = eigen.vectors[ev];    en = eigen.values[ev]
            parity    = Basics.plus;    F = AngularJ64(0);     MF = AngularM64(0)
            for  r = 1:length(hfsBasis.vectors)
                if  abs(evector[r]) > 1.0e-6
                    parity = hfsBasis.vectors[r].levelJ.parity
                    F      = hfsBasis.vectors[r].F;     MF = AngularM64(hfsBasis.vectors[r].F);   break    end
            end
            newlevelF = Hfs.IJF_Level(nm.spinI, F, MF, parity, en, true, hfsBasis, evector)
            push!( levelFs, newlevelF)
        end
        hfsMultiplet = Hfs.IJF_Multiplet("hyperfine", levelFs)

        return( hfsMultiplet )
    end

    """
    `Hfs.computeInteractionMatrix(basis::Basis, grid::Radial.Grid, settings::Hfs.Settings)`
        ... to compute the T^1 and/or T^2 interaction matrices for the given basis, i.e. (<csf_r || T^(n)) || csf_s>).
            An im::Hfs.InteractionMatrix is returned.
    """
    function  computeInteractionMatrix(basis::Basis, grid::Radial.Grid, settings::Hfs.Settings)
        #
        ncsf = length(basis.csfs);    matrixT1 = zeros(ncsf,ncsf);    matrixT2 = zeros(ncsf,ncsf)
        #
        if  settings.calcT1
            calcT1 = true;    matrixT1 = zeros(ncsf,ncsf)
            for  r = 1:ncsf
                for  s = 1:ncsf
                    if  basis.csfs[r].parity  != basis.csfs[s].parity   continue    end
                    #
                    wa = compute("angular coefficients: 1-p, Grasp92", 0, 1, basis.csfs[r], basis.csfs[s])
                    for  coeff in wa
                        ja   = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                        jb   = Basics.subshell_2j(basis.orbitals[coeff.b].subshell)
                        tamp = InteractionStrength.hfs_t1(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid)
                        matrixT1[r,s] = matrixT1[r,s] + coeff.T * tamp
                    end
                end
            end
        else
            calcT1 = false;    matrixT1 = zeros(2,2)
        end
        #
        if  settings.calcT2
            calcT2 = true;    matrixT2 = zeros(ncsf,ncsf)
            for  r = 1:ncsf
                for  s = 1:ncsf
                    if  basis.csfs[r].parity  != basis.csfs[s].parity   continue    end
                    #
                    wa = compute("angular coefficients: 1-p, Grasp92", 0, 2, basis.csfs[r], basis.csfs[s])
                    for  coeff in wa
                        ja   = Basics.subshell_2j(basis.orbitals[coeff.a].subshell)
                        jb   = Basics.subshell_2j(basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.hfs_t2(basis.orbitals[coeff.a], basis.orbitals[coeff.b], grid)
                        matrixT2[r,s] = matrixT2[r,s] + coeff.T * tamp
                    end
                end
            end
        else
            calcT2 = false;    matrixT2 = zeros(2,2)
        end
        #
        im = Hfs.InteractionMatrix(calcT1, calcT2, matrixT1, matrixT2)
        #
        return( im )
    end

    """
    `Hfs.computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::Hfs.Settings; output=true)`
        ... to compute (as selected) the HFS A and B parameters as well as hyperfine energy splittings for the levels
            of the given multiplet and as specified by the given settings. The results are printed in neat tables to
            screen and, if requested, an arrays{Hfs.Outcome,1} with all the results are returned.
    """
    function computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::Hfs.Settings; output=true)
        println("")
        printstyled("Hfs.computeOutcomes(): The computation of the Hyperfine amplitudes and parameters starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------------ \n", color=:light_green)
        println("")
        outcomes = Hfs.determineOutcomes(multiplet, settings)
        # Display all selected levels before the computations start
        if  settings.printBefore    Hfs.displayOutcomes(outcomes)    end
        # Calculate all amplitudes and requested properties
        im = computeInteractionMatrix(multiplet.levels[1].basis, grid, settings)
        newOutcomes = Hfs.Outcome[]
        for  outcome in outcomes
            newOutcome = Hfs.computeAmplitudesProperties(outcome, grid, settings, im)
            push!( newOutcomes, newOutcome)
        end
        # Print all results to screen
        Hfs.displayResults(stdout, newOutcomes, nm, settings)
        # Compute and display the non-diagonal hyperfine amplitudes, if requested
        if  settings.calcNondiagonal    Hfs.displayNondiagonal(stdout, multiplet, grid, settings)   end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary
            Hfs.displayResults(iostream, newOutcomes, nm, settings)
            if  settings.calcNondiagonal    Hfs.displayNondiagonal(iostream, multiplet, grid, settings)   end
        end
        #
        if    output    return( newOutcomes )
        else            return( nothing )
        end
    end

    """
    `Hfs.defineHyperfineBasis(multiplet::Multiplet, nm::Nuclear.Model, settings::Hfs.Settings)`
        ... to define/set-up an atomic hyperfine (IJF-coupled) basis for the given electronic multipet;
            a hfsBasis::IJF_Basis is returned.
    """
    function  defineHyperfineBasis(multiplet::Multiplet, nm::Nuclear.Model, settings::Hfs.Settings)
        function  display_ijfVector(i::Int64, vector::Hfs.IJF_Vector)
            si = string(i);   ni = length(si);    sa = repeat(" ", 5);    sa = sa[1:5-ni] * si * ")  "
            sa = sa * "[" * string( LevelSymmetry(vector.levelJ.J, vector.levelJ.parity) ) * "] " * string(vector.F) * repeat(" ", 4)
            return( sa )
        end

        vectors = IJF_Vector[]
        #
        for i = 1:length(multiplet.levels)
            Flist = Basics.oplus(nm.spinI, multiplet.levels[i].J)
            for  F in Flist
                push!(vectors,  IJF_Vector(nm.spinI, F, multiplet.levels[i]) )
            end
        end
        #
        hfsBasis = IJF_Basis(vectors, multiplet.levels[1].basis)
        #
        println(" ")
        println("  Construction of a atomic hyperfine (IJF-coupled) basis of dimension $(length(hfsBasis.vectors)), based on " *
                "$(length(multiplet.levels)) electronic levels and nuclear spin $(hfsBasis.vectors[1].I).")
        println("  Basis vectors are defined as [J^P] F: \n")
        sa = "  "
        for  k = 1:length(hfsBasis.vectors)
            if  length(sa) > 100    println(sa);   sa = "  "   end
            sa = sa * display_ijfVector(k, hfsBasis.vectors[k])
        end
        if   length(sa) > 5    println(sa)   end
        println(" ")

        return( hfsBasis )
    end

    """
    `Hfs.determineOutcomes(multiplet::Multiplet, settings::Hfs.Settings)`
        ... to determine a list of Outcomes's for the computation of HFS A- and B-parameters for the given multiplet.
            It takes into account the particular selections and settings. An Array{Hfs.Outcome,1} is returned. Apart from the
            level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineOutcomes(multiplet::Multiplet, settings::Hfs.Settings)
        outcomes = Hfs.Outcome[]
        for  level  in  multiplet.levels
            if  Basics.selectLevel(level, settings.levelSelection)
                push!( outcomes, Hfs.Outcome(level, 0., 0., 0., 0., AngularJ64(0), false, false) )
            end
        end
        return( outcomes )
    end

    """
    `Hfs.displayNondiagonal(stream::IO, multiplet::Multiplet, grid::Radial.Grid, settings::Hfs.Settings)`
        ... to compute and display all non-diagonal hyperfine amplitudes for the selected levels. A small neat table of
            all (pairwise) hyperfine amplitudes is printed but nothing is returned otherwise.
    """
    function  displayNondiagonal(stream::IO, multiplet::Multiplet, grid::Radial.Grid, settings::Hfs.Settings)
        # Determine pairs to be calculated
        pairs = Tuple{Int64,Int64}[]
        for  (f, fLevel)  in  enumerate(multiplet.levels)
            for  (i, iLevel)  in  enumerate(multiplet.levels)
                if  Basics.selectLevel(fLevel, settings.levelSelection)  &&   Basics.selectLevel(iLevel, settings.levelSelection)
                    push!( pairs, (f,i) )
                end
            end
        end
        #
        nx = 107
        println(stream, " ")
        println(stream, "  Selected (non-) diagonal hyperfine amplitudes:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  "
        sa = sa * TableStrings.center(10, "Level_f"; na=2)
        sa = sa * TableStrings.center(10, "Level_i"; na=2)
        sa = sa * TableStrings.center(10, "J^P_f";   na=3)
        sa = sa * TableStrings.center(57, "T1   --   Amplitudes   --   T2"; na=4);
        sa = sa * TableStrings.center(10, "J^P_f";   na=4)
        println(stream, sa);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  (f,i) in pairs
            sa   = "  ";
            sa   = sa * TableStrings.center(10, string(f); na=2)
            sa   = sa * TableStrings.center(10, string(i); na=2)
            symf = LevelSymmetry( multiplet.levels[f].J, multiplet.levels[f].parity)
            symi = LevelSymmetry( multiplet.levels[i].J, multiplet.levels[i].parity)
            sa   = sa * TableStrings.center(10, string(symf); na=4)
            T1   = Hfs.amplitude("T^(1) amplitude", multiplet.levels[f], multiplet.levels[i], grid, printout=false)
            T2   = Hfs.amplitude("T^(2) amplitude", multiplet.levels[f], multiplet.levels[i], grid, printout=false)
            sa   = sa * @sprintf("%.5e %s %.5e", T1.re, "  ", T1.im) * "    "
            sa   = sa * @sprintf("%.5e %s %.5e", T2.re, "  ", T2.im) * "    "
            sa   = sa * TableStrings.center(10, string(symi); na=4)
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `Hfs.displayOutcomes(outcomes::Array{Hfs.Outcome,1})`
        ... to display a list of levels that have been selected for the computations A small neat table of all
            selected levels and their energies is printed but nothing is returned otherwise.
    """
    function  displayOutcomes(outcomes::Array{Hfs.Outcome,1})
        nx = 43
        println(" ")
        println("  Selected HFS levels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.Jlevel.J, outcome.Jlevel.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.Jlevel.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.Jlevel.energy)) * "    "
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        println(" ")
        #
        return( nothing )
    end

    """
    `Hfs.displayResults(stream::IO, outcomes::Array{Hfs.Outcome,1}, nm::Nuclear.Model, settings::Hfs.Settings)`
        ... to display the energies, A- and B-values, Delta E_F energy shifts, etc. for the selected levels. All nuclear
            parameters are taken from the nuclear model. A neat table is printed but nothing is returned otherwise.
    """
    function  displayResults(stream::IO, outcomes::Array{Hfs.Outcome,1}, nm::Nuclear.Model, settings::Hfs.Settings)
        nx = 117
        println(stream, " ")
        println(stream, "  HFS parameters and amplitudes:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(31, "A/mu [mu_N]  --  HFS  --  B/Q [barn]"; na=4);
        sb = sb * TableStrings.center(31, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(32, "T1 -- Amplitudes -- T2"    ; na=4);        sb = sb * TableStrings.hBlank(36)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.Jlevel.J, outcome.Jlevel.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.Jlevel.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            energy = outcome.Jlevel.energy
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))           * "      "
            wa = Defaults.convertUnits("energy: from atomic", outcome.AIoverMu) / nm.spinI.num * nm.spinI.den
            wa = wa / (2 * 1836.15267 ) * Defaults.getDefaults("alpha")  # take 2m_p and alpha into account
            sa = sa * @sprintf("%.8e", wa)                                                   * "    "
            wa = Defaults.convertUnits("energy: from atomic", outcome.BoverQ)
            wa = wa / Defaults.convertUnits("cross section: from atomic to barn", 1.0)  # take Q [barn] into account
            sa = sa * @sprintf("%.8e", wa)                                                   * "    "
            sa = sa * @sprintf("%.8e %s %.8e", outcome.amplitudeT1.re, "  ", outcome.amplitudeT2.re) * "    "
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        #
        if !settings.printDeltaEF   return( nothing )   end
        #
        # Printout the Delta E_F energy shifts of the hyperfine levels |alpha F> with regard to the (electronic) levels |alpha J>
        nx = 90
        println(stream, " ")
        println(stream, "  HFS Delta E_F energy shifts with regard to the (electronic) level energies E_J:")
        println(stream, " ")
        println(stream, "    Nuclear spin I:                        $(nm.spinI) ")
        println(stream, "    Nuclear magnetic-dipole moment    mu = $(nm.mu)    ")
        println(stream, "    Nuclear electric-quadrupole moment Q = $(nm.Q)     ")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);                            sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(10, "F^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Delta E_F"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(14, "C factor"; na=4);
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.Jlevel.J, outcome.Jlevel.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.Jlevel.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            energy = outcome.Jlevel.energy
            #
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))           * "    "
            J     = AngularMomentum.oneJ(outcome.Jlevel.J);   spinI = AngularMomentum.oneJ(nm.spinI)
            Flist = Basics.oplus(nm.spinI, outcome.Jlevel.J)
            first = true
            for  Fang in Flist
                Fsym    = LevelSymmetry(Fang, outcome.Jlevel.parity)
                F       = AngularMomentum.oneJ(Fang)
                Cfactor = F*(F+1) - J*(J+1) - spinI*(spinI+1)
                energy  = outcome.AIoverMu * nm.mu / spinI * Cfactor / 2.
                if  abs(outcome.BoverQ) > 1.0e-10
                    energy  = energy +  outcome.BoverQ * nm.Q * 3/4 * (Cfactor*(Cfactor+1) - spinI*(spinI+1)*J*(J+1) ) /
                                        ( 2spinI*(2spinI-1)*J*(2J-1) )
                end
                sb = TableStrings.center(10, string(Fsym); na=2)
                sb = sb * TableStrings.flushright(16, @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))) * "    "
                sb = sb * TableStrings.flushright(12, @sprintf("%.5e", Cfactor))
                #
                if   first    println(stream,  sa*sb );   first = false
                else          println(stream,  TableStrings.hBlank( length(sa) ) * sb )
                end
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
println("hello world")
#> hello world

ext: jl

max_stars_repo_name: lavirlifiliol/hello-world

max_issues_repo_name: lavirlifiliol/hello-world

max_forks_repo_name: lavirlifiliol/hello-world


text:
firstpass(str) = match(r"^\s*(\w+)\s+(-?\d+)\s+(-?\d+)\s*$", str)

secondpass(::Nothing, str) = Symbol(strip(str)), nothing
function secondpass(m, _)
    poi, _x, _y = m.captures
    x = parse(Float64, _x)
    y = parse(Float64, _y)
    return Symbol(poi), Space(x, y)
end

parsepoi(str) = secondpass(firstpass(str), str)

function parsepois(row)
    poi_names = Symbol[]
    expected_locations = Dict{Symbol, Space}()
    for poi in split(row, ',')
        poi_name, xy = parsepoi(poi)
        push!(poi_names, poi_name)
        if !isnothing(xy)
            expected_locations[poi_name] = xy
        end
    end
    (; poi_names, expected_locations)
end

function _adjust_expected(pois, expected_locations)
    e = @LArray Matrix(hcat(values(expected_locations)...)) (; Dict(k => (2i-1:2i) for (i,k) in enumerate(keys(expected_locations)))...)
    r = Matrix(hcat((only(space(pois[k])) for k in keys(expected_locations))...))
    R, μ = getR(e, r)
    Dict(k => Space(R * e[k] + μ) for k in keys(expected_locations))
end

function getR(e, r)
    e .-= mean(e, dims = 2)
    μr = mean(r, dims = 2)
    r .-= μr
    F = svd(e * r')
    R = F.V * F.U'
    return R, μr
end

adjust_expected(pois, expected_locations) = length(expected_locations) ≤ 1 ?  Dict(k => only(space(pois[k])) for k in keys(expected_locations)) : _adjust_expected(pois, expected_locations)

function flipy!(pois)
    ymax = maximum(maximum(last, space(v)) for v in values(pois))
    for v in values(pois), row in LazyRows(v.xyt)
        x2, y = row.xy
        row.xy = Space(x2, ymax - y)
    end
end


text:
# Simplicial Complex

"""Simplicial Complex Type

Simplex parameter type `P` must have implemented total order function `isless`.
"""
mutable struct SimplicialComplex{P} <: AbstractComplex
    cells::Dict{Int,Vector{Simplex{P}}}   # cells per dimension
    # order::Dict{Pair{Int,Int},Int}      # order of cells
end
Base.show(io::IO, cplx::SimplicialComplex) = print(io, "SimplicialComplex($(size(cplx)))")
Base.copy(cplx::SimplicialComplex) = SimplicialComplex(deepcopy(cplx.cells))

# ---------------
# Private Methods
# ---------------

"""Add a simplex to the complex (as is) and return a complex size, a dimension of simplex and its index in it

This function **doesn't** add missing faces of the simplex to the complex. Use `addsimplex!` function for instead.
"""
function addsimplex(cplx::SimplicialComplex{P}, splx::Simplex{P}) where {P}
    d = dim(splx)
    d < 0 && return (0,d,0)
    !haskey(cplx.cells, d) && setindex!(cplx.cells, Simplex{P}[], d)
    i = length(cplx.cells[d])+1
    splx[:index] = i
    push!(cplx.cells[d], splx)
    # cplx.order[d=>i] = length(cplx.order)+1
    return splx
end

"""Add a simplex to the complex and all of its faces recursivly and return a complex size, a dimension of simplex and its index in it"""
function addsimplex!(cplx::SimplicialComplex{P}, splx::Simplex{P}) where {P}
    # cache already
    added = Set{Simplex{P}}()

    toprocess = Simplex{P}[]
    ret = Simplex{P}[]
    addedidxs = NTuple{3,Int}[]

    # add simplex to complex
    s = addsimplex(cplx, splx)
    push!(addedidxs, (sum(size(cplx)), dim(s), s[:index]))
    push!(ret, s)
    for f in faces(splx) # add simplex faces for processing
        push!(toprocess, f)
    end

    # process simplex faces and add missing simplexes
    while(length(toprocess) > 0)
        tmp = pop!(toprocess) # processing simplex faces
        d = dim(tmp)

        cplx[tmp, d] <= size(cplx, d) && continue # skip if already in complex
        tmp in added && continue # skip if already processed

        s = addsimplex(cplx, tmp) # add simples to the complex
        push!(addedidxs, (sum(size(cplx)), dim(s), s[:index]))
        push!(added, tmp) # mark as processed
        push!(ret, s)

        dim(tmp) == 0 && continue # do not create faces for 0-simplexes
        for f in faces(tmp) # add simplex faces for processing
            push!(toprocess, f)
        end
    end

    return ret
end

#
# Public Interface
#

celltype(cplx::SimplicialComplex{P}) where {P} = Simplex{P}

function cells(cplx::SimplicialComplex)
    CCT = valtype(cplx.cells)
    length(cplx.cells) == 0 && return CCT[] # no cell in complex
    dims = maximum(keys(cplx.cells))
    CCT[haskey(cplx.cells, d) ? cplx.cells[d] : CCT() for d in 0:dims]
end

cells(cplx::SimplicialComplex{P}, d::Int) where {P} = get(cplx.cells, d,  nothing)

dim(cplx::SimplicialComplex) = length(size(cplx))-1

function boundary(cplx::SimplicialComplex, idx::Int, d::Int, ::Type{R}) where {R}
    ch = Chain(d-1, R)
    if d == 0
        setdim!(ch, -1)
        return ch
    end
    pos = one(R)
    neg = -one(R)

    splx = cplx[idx, d]
    splx === nothing && return ch

    sgn = true
    for face in faces(splx)
        fidx = cplx[face, d-1]
        push!(ch, (sgn ? pos : neg), fidx)
        sgn = !sgn
    end

    return ch
end

function coboundary(cplx::SimplicialComplex, idx::Int, d::Int, ::Type{R}) where {R}

    d == 0 && return Chain(d+1, R)

    cbd = Dict{Int,Chain{R}}()
    for i in 1:size(cplx,d-1)
        cbd[i] = Chain(d+1, R)
    end

    for i in 1:size(cplx,d)
        bd = boundary(cplx, i, d, R)
        for (c, el) in bd
            push!(cbd[el], c, i)
        end
    end

    return cbd[idx]
end

Base.push!(cplx::SimplicialComplex{P}, splx::Simplex{P}; recursive=false) where {P} =
    recursive ? addsimplex!(cplx, splx) : [addsimplex(cplx, splx)]

#
# Constructors
#

SimplicialComplex(::Type{P}) where P = SimplicialComplex(Dict{Int,Vector{Simplex{P}}}())
(::Type{SimplicialComplex{P}})() where P = SimplicialComplex(P)

function SimplicialComplex(splxs::Simplex{P}...) where {P}
    cplx = SimplicialComplex(P)
    added = Set{Simplex}()
    toprocess = Vector{Simplex}()
    for splx in splxs
        push!(toprocess, splx) # add simplex to process its faces
        while(length(toprocess) > 0)
            tmp = pop!(toprocess) # processing simplex faces
            tmp in added && continue # skip if already processed
            addsimplex(cplx, tmp) # add simples to complex
            push!(added, tmp) # mark as processed
            for f in faces(tmp) # add simplex faces for processing
                push!(toprocess, f)
            end
        end
    end

    return cplx
end

SimplicialComplex(splxs::Vector{Simplex}) = SimplicialComplex(splxs...)

#
# Miscellaneous
#

function Base.read(io::IO, ::Type{SimplicialComplex{P}}) where {P}
    splxs = Simplex[]
    while !eof(io)
        l = chomp(readline(io))
        si = map(e->parse(P, e), split(l, ' '))
        push!(splxs, Simplex(si...))
    end
    return SimplicialComplex(splxs)
end

function Base.write(io::IO, cplx::SimplicialComplex)
    zsplxs = Dict{Int,Int}()
    for s in cplx.cells[0]
        zsplxs[s.vs[]] = s.idx
        write(io, "$(s.idx)")
        write(io, 0x0A)
    end
    for d in 1:dim(cplx)
        for s in cplx.cells[d]
            for (j,i) in enumerate(s.vs)
                write(io, "$(zsplxs[i])")
                if j>d
                    write(io, 0x0A)
                else
                    write(io, ' ')
                end
            end
        end
    end
    return
end


text:
# FeynmanTheory.jl

"""
frohlichalpha(ε_Inf,ε_S,freq,m_eff)

    Calculates the Frohlich alpha parameter, for a given dielectric constant,
    frequency (f) of phonon in Hertz, and effective mass (in units of the
    bare electron mass).

    See Feynman 1955:
    http://dx.doi.org/10.1103/PhysRev.97.660

"""
function frohlichalpha(ε_optic,ε_static,freq,m_eff)
    ω=freq*2*pi #frequency to angular velocity
    # Note: we need to add a 4*pi factor to the permitivity of freespace.
    # This gives numeric agreement with literature values.  This is required as
    # the contemporary 1950s and 1960s literature implicitly used atomic units,
    # where the electric constant ^-1 has this factor baked in, k_e=1/(4πϵ_0).
    α= 0.5/(4*π*ε_0) *              # Units: m/F
       (1/ε_optic - 1/ε_static) *   # Units: none
       (q^2/(hbar*ω)) *             # Units: F
       sqrt(2*me*m_eff*ω/ħ)         # Units: 1/m
    return (α)
end

#####
# Athermal (Feynman 1955) model
# Set up equations for the polaron free energy, which we will variationally improve upon

# Integrand of (31) in Feynman I (Feynman 1955, Physical Review, "Slow electrons...")
fF(τ,v,w)=(w^2 * τ + (v^2-w^2)/v*(1-exp(-v*τ)))^-0.5 * exp(-τ)
# (31) in Feynman I
AF(v,w,α)=π^(-0.5) * α*v * quadgk(τ->fF(τ,v,w),0,Inf)[1]
# (33) in Feynman I
F(v,w,α)=(3/(4*v))*(v-w)^2-AF(v,w,α)

# Let's wrap the Feynman athermal variation approximation in a simple function
"""
    feynmanvw(α; v=7.0, w=6.0)

    Calculate v and w variational polaron parameters,
    for the supplied α Frohlich coupling.
    This version uses the original athermal action (Feynman 1955).
	Returns v,w.
"""
function feynmanvw(α; v=7.0, w=6.0) # v,w defaults
    initial=[v,w]
    # Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
    lower=[0.1,0.1]
    upper=[1000.0,1000.0]

    myf(x) = F(x[1],x[2],α) # Wraps the function so just the two variational params are exposed

    # Now updated to use Optim.jl > 0.15.0 call signature (Julia >0.6 only)
    res=optimize(OnceDifferentiable(myf, initial; autodiff = :forward),
                 lower, upper, initial, Fminbox( BFGS() ) )
    # specify Optim.jl optimizer. This is doing all the work.

    v,w=Optim.minimizer(res)

    return v,w
end

# Hellwarth et al. 1999 PRB - Part IV; T-dep of the Feynman variation parameter

# Originally a Friday afternoon of hacking to try and implement the T-dep electron-phonon coupling from the above PRB
# Which was unusually successful! And more or less reproduced Table III

# In Julia we have 'Multiple dispatch', so let's just construct the free
# energies (temperature-dependent) with the same name as before, but withthe
# thermodynamic beta where required

# Define Osaka's free-energies (Hellwarth1999 version) as Julia functions
# Equation numbers follow above Hellwarth et al. 1999 PRB
# 62b
A(v,w,β)=3/β*( log(v/w) - 1/2*log(2*π*β) - log(sinh(v*β/2)/sinh(w*β/2)))
# 62d
Y(x,v,β)=1/(1-exp(-v*β))*(1+exp(-v*β)-exp(-v*x)-exp(v*(x-β)))
# 62c integrand
f(x,v,w,β)=(exp(β-x)+exp(x))/(w^2*x*(1-x/β)+Y(x,v,β)*(v^2-w^2)/v)^(1/2)
# 62c
B(v,w,β,α) = α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->f(x,v,w,β),0,β/2)[1]
# 62e
C(v,w,β)=3/4*(v^2-w^2)/v * (coth(v*β/2)-2/(v*β))
# 62a
F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β))

# Can now evaluate, e.g.
# F(v,w,β,α)=F(7.2,6.5,1.0,1.0)
# BUT - this is just the objective function! (The temperature-dependent free-energy.) Not the optimised parameters.
# Also there's a scary numeric integration (quadgk) buried within...

# Unbundled version of the finite temperature one
"""
    feynmanvw(α, βred; v=7.1, w=6.5, verbose::Bool=false)

    Calculate v and w variational polaron parameters,
    for the supplied α Frohlich coupling and βred reduced thermodynamic temperature.
    This uses the Osaka finite temperature action.
    Returns v,w.
"""
function feynmanvw(α, βred; v=7.1, w=6.5, verbose::Bool=false) # v,w defaults
    # Initial v,w to use
    initial=[v,w]
    # Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
    lower=[0.1,0.1]
    upper=[100.0,100.0]

    myf(x) = F(x[1],x[2],βred,α) # Wraps the function so just the two variational params are exposed

    # Now updated to use Optim > 0.15.0 call signature (Julia >0.6 only)
    res=optimize(OnceDifferentiable(myf, initial; autodiff = :forward),
                 lower, upper, initial, Fminbox( BFGS() ) )
    # specify Optim.jl optimizer. This is doing all the work.

    if Optim.converged(res) == false
        print("\tWARNING: Failed to converge to v,w soln? : ",Optim.converged(res) )
    end

    if verbose
        println()
        show(res)
    end

    v,w=Optim.minimizer(res)
    return v,w
end


text:
export databuf_t, databuf_desc_t, databuf_data, databuf_create, databuf_clear, databuf_attach,
        databuf_detach, check_databuf, databuf_wait_filled, databuf_wait_free, databuf_set_filled,
        databuf_set_free

abstract type HashpipeDatabuf end

"""
    databuf_t <: HashpipeDatabuf

Contain data pertaining to a Hashpipe databuffer.
"""
struct databuf_t <: HashpipeDatabuf
    data_type::NTuple{64, UInt8}
    header_size::Int # May need to change to Csize_t
    block_size::Int # May need to change to Csize_t
    n_block::Cint
    shmid::Cint
    semid::Cint
end
# From C code:
# These typedefs are used to declare pointers to a pipeline thread's init and
# run functions.
# typedef int (* initfunc_t)(thread_args_t *);
# typedef void * (* runfunc_t)(thread_args_t *);

# This typedefs are used to declare pointers to a pipline thread's data buffer
# create function.
# typedef databuf_t * (* databuf_createfunc_t)(int, int);

# typedef struct {
#   databuf_createfunc_t create;
# } databuf_desc_t;

"""
    struct databuf_desc_t

Contain the pointer to the databuf create function.
"""
struct databuf_desc_t
    """
    C-compatible pointer to Julia databuf create function.
    Generate using @cfunction macro.
    """
    create::Ptr{Cvoid}
end

"""
    databuf_data(p_databuf::Ptr{databuf_t}, block_id::Int)

Return the pointer to the associated databuffer's block's data.
"""
function databuf_data(p_databuf::Ptr{databuf_t}, block_id::Int)
    p_data::Ptr{UInt8} = ccall((:hashpipe_databuf_data, libhashpipe),
                            Ptr{UInt8}, (Ptr{status_t}, Int8), p_databuf, block_id)
    return p_data
end

"""
    databuf_create(instance_id::Int, db_id::Int, header_size::Int, block_size::Int, n_block::Int)

Create a databuffer with given parameters.
"""
function databuf_create(instance_id::Int, db_id::Int,
            header_size::Int, block_size::Int, n_block::Int)
    p_databuf::Ptr{databuf_t} =
            ccall((:hashpipe_databuf_create, libhashpipe),
                Ptr{databuf_t},
                (Int8, Int8, Int, Int, Int),
                instance_id, db_id, header_size, block_size, n_block)
    return p_databuf
end

"""
    databuf_clear(p_databuf::Ptr{databuf_t})

Clear the data in a databuf.
"""
function databuf_clear(p_databuf::Ptr{databuf_t})
    ccall((:hashpipe_databuf_clear, libhashpipe),
            Cvoid, (Ptr{status_t},), p_databuf)
    return nothing
end

"""
    databuf_attach(instance_id::Int, db_id::Int)

Attach a databuf to a Hashpipe instance.
"""
function databuf_attach(instance_id::Int, db_id::Int)
    p_databuf::Ptr{databuf_t} = ccall((:hashpipe_databuf_attach, libhashpipe),
                    Ptr{databuf_t}, (Int8, Int8), instance_id, db_id)
    return p_databuf
end

"""
    databuf_detach(instance_id::Int, db_id::Int)

Detach a databuf from a Hashpipe instance.
"""
function databuf_detach(p_databuf::Ptr{databuf_t})
    error::Int = ccall((:hashpipe_databuf_detach, libhashpipe),
                    Int, (Ptr{databuf_t},), p_databuf)
    return error
end

"""
    check_databuf(instance_id=0, db_id=1)

Display databuf information of a given databuf of a Hashpipe instance.
"""
function check_databuf(instance_id::Int = 0, db_id::Int = 1)
    p_databuf = databuf_attach(instance_id, db_id)
    if p_databuf == C_NULL
        println("Error attaching to databuf $db_id (may not exist).")
        return nothing
    end
    println("--- Databuf $db_id Stats ---")
    display(p_databuf)
    return nothing
end

"""
    databuf_block_status(p_databuf::Ptr{databuf_t}, block_id::Int)

Return the status of the selected data block of the given databuffer.
"""
function databuf_block_status(p_databuf::Ptr{databuf_t}, block_id::Int)
    block_status::Int = ccall((:hashpipe_databuf_block_status, libhashpipe),
                    Int, (Ptr{databuf_t}, Int), p_databuf, block_id)
    return block_status
end

"""
    databuf_total_status(p_databuf::Ptr{databuf_t})

Return the total lock status for the given databuffer.
"""
function databuf_total_status(p_databuf::Ptr{databuf_t})
    total_status::UInt64 = ccall((:hashpipe_databuf_total_status, libhashpipe),
                    Int, (Ptr{databuf_t},), p_databuf)
    return total_status
end

"""
    databuf_total_mask(p_databuf::Ptr{databuf_t})
"""
function databuf_total_mask(p_databuf::Ptr{databuf_t})
    total_mask::UInt64 = ccall((:hashpipe_databuf_total_mask, libhashpipe),
                    UInt64, (Ptr{databuf_t},), p_databuf)
    return total_mask
end

# Databuf locking functions
# Each block in the buffer can be marked as free or filled.
# The "wait" functions block (i.e. sleep) until the specified state happens.
# The "busywait" functions busy-wait (i.e. do NOT sleep)
# until the specified state happens.  The "set" functions
# put the buffer in the specified state, returning error if
# it is already in that state.
"""
    databuf_wait_filled(p_databuf::Ptr{databuf_t}, block_id::Int)

Wait for the given block of data to be filled.

See also: [`databuf_wait_free`](@ref), [`databuf_set_filled`](@ref), [`databuf_set_free`](@ref)
"""
function databuf_wait_filled(p_databuf::Ptr{databuf_t}, block_id::Int)
    error::Int = ccall((:hashpipe_databuf_wait_filled, libhashpipe),
                    Int, (Ptr{databuf_t}, Int), p_databuf, block_id)
    return error
end

"""
    databuf_wait_free(p_databuf::Ptr{databuf_t}, block_id::Int)

Wait for the given block of data to be freed.

See also: [`databuf_wait_filled`](@ref), [`databuf_set_filled`](@ref), [`databuf_set_free`](@ref)
"""
function databuf_wait_free(p_databuf::Ptr{databuf_t}, block_id::Int)
    error::Int = ccall((:hashpipe_databuf_wait_free, libhashpipe),
                    Int, (Ptr{databuf_t}, Int), p_databuf, block_id)
    return error
end

"""
    databuf_set_filled(p_databuf::Ptr{databuf_t}, block_id::Int)

Set the given block of data as filled. Return an error if the block is already filled.

See also: [`databuf_wait_filled`](@ref), [`databuf_wait_free`](@ref), [`databuf_set_free`](@ref)
"""
function databuf_set_filled(p_databuf::Ptr{databuf_t}, block_id::Int)
    error::Int = ccall((:hashpipe_databuf_set_filled, libhashpipe),
                    Int, (Ptr{databuf_t}, Int), p_databuf, block_id)
    return error
end

"""
    databuf_set_free(p_databuf::Ptr{databuf_t}, block_id::Int)

Set the given block of data as free. Return an error if the block is already free.

See also: [`databuf_wait_filled`](@ref), [`databuf_wait_free`](@ref), [`databuf_set_filled`](@ref)
"""
function databuf_set_free(p_databuf::Ptr{databuf_t}, block_id::Int)
    error::Int = ccall((:hashpipe_databuf_set_free, libhashpipe),
                    Int, (Ptr{databuf_t}, Int), p_databuf, block_id)
    return error
end

"""
    Base.display(d::Hashpipe.databuf_t)

Display a Hashpipe databuffer nicely in REPL.
"""
function Base.display(d::Hashpipe.databuf_t)
    # Convert Ntuple to array and strip 0s before converting to string
    data_type_string = String(filter(x->x!=0x00, collect(d.data_type)))
    println("Data Type: $(data_type_string)")
    println("Header Size: $(d.header_size)")
    println("Num Blocks: $(d.n_block)")
    println("Block Size: $(d.block_size)")
    println("shmid: $(d.shmid)")
    println("semid: $(d.semid)")
    return nothing
end

"""
    Base.display(p::Ptr{databuf_t})

Display Hashpipe databuffer from pointer nicely in REPL.
"""
function Base.display(p::Ptr{databuf_t})
    databuf::databuf_t = unsafe_wrap(Array, p, 1)[]
    display(databuf)
    return nothing
end


text:
lines = readlines(open("input.txt"))
#lines = readlines(open("example.txt"))
data = cat([[x=='#' for x in line] for line in lines]..., dims=2)'

state = zeros(Bool, 403, 5, 5)
state[202, :, :] = copy(data)

function neighbours(state, d, i, j)
    nneigh = 0
    # up
    if i == 1
        nneigh += state[d-1, 2, 3]
    elseif (i == 4) && (j == 3)
        nneigh += sum(state[d+1, 5, 1:5])
    else
        nneigh += state[d, i-1, j]
    end
    # down
    if i == 5
        nneigh += state[d-1, 4, 3]
    elseif (i == 2) && (j == 3)
        nneigh += sum(state[d+1, 1, 1:5])
    else
        nneigh += state[d, i+1, j]
    end
    # left
    if j == 1
        nneigh += state[d-1, 3, 2]
    elseif (j == 4) && (i == 3)
        nneigh += sum(state[d+1, 1:5, 5])
    else
        nneigh += state[d, i, j-1]
    end
    # right
    if j == 5
        nneigh += state[d-1, 3, 4]
    elseif (j == 2) && (i == 3)
        nneigh += sum(state[d+1, 1:5, 1])
    else
        nneigh += state[d, i, j+1]
    end
    nneigh
end

function step!(next, state)
    for d=2:size(state, 1)-1
        for i=1:5
            for j=1:5
                if (i==3) && (j==3)
                    continue
                end
                nneigh = neighbours(state, d, i, j)
                if state[d, i, j]
                    next[d, i, j] = nneigh == 1
                else
                    next[d, i, j] = (nneigh == 1) || (nneigh == 2)
                end
            end
        end
    end
end

function run(state, nsteps)
    state = copy(state)
    next = copy(state)
    for i=1:nsteps
        step!(next, state)
        state, next = next, state
    end
    state
end

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
### Resampling: INPUT state_matrix, transition_matrix, C, T, number_zones, p_drive, p_dest, datamatrix, distance_matrix_km and OUTPUT state_matrix, transition_matrix

function resampling(state_matrix, transition_matrix, C, number_zones, p_drive, p_dest, datamatrix, distance_matrix_km)

        ## Sampling an entire day with imporved intitial value.
        for t=1:T

            # Sampling if a car drives
            for i=1:C
                origin = state_matrix[i,t];
                RndVar = rand();
                driving_probability = p_drive[origin,t];
                if(RndVar <= driving_probability)
                    drive = 1;
                else
                    drive = 0;
                    transition_matrix[i,t,2]= origin;
                end
                transition_matrix[i,t,1] = drive;
            end

            # Sampling travel destinations
            for i=1:C
                drive = transition_matrix[i,t,1];
                if (drive == 1)
                    RndVar = rand();
                    range_up = 0;
                    range_low = 0;
                    destination = 0;
                    origin = state_matrix[i,t];
                    distribution = p_dest[origin, :, t];
                    if(sum(distribution) == 0)
                        destination = origin;
                    else
                        for j=1:number_zones
                            range_up = range_up + distribution[j];
                            if(range_low < RndVar <= range_up)
                                destination = j;
                                break;
                            end
                            range_low = range_up;
                        end
                    end
                    transition_matrix[i,t,2]= destination;
                end
            end

            # sampling travel times and distances
            for i=1:C
                drive = transition_matrix[i,t,1];
                if(drive == 1)
                    origin = state_matrix[i,t,1];
                    destination = round.(Int, transition_matrix[i,t,2]);
                    if(origin == destination) # this is the case that a car drives to a destination within the origin zone
                        transition_matrix[i,t,3] = 5 * 60; # assuming arbitrarily a driving time of 5 minutes given in seconds
                        transition_matrix[i,t,4] = 1; # assuming arbitraliy a driving distance of 1 km given in km.
                    else
                        # sampling travel times
                        mean = datamatrix[origin, destination, t, 1];
                        std_deviation = datamatrix[origin, destination, t, 2];
                        if(std_deviation == 0) # catches the case of missing std deviation data where it is set to 0 in the original datasets
                            std_deviation = 0.1 * mean;
                        end
                        travel_time = rand(Truncated(Normal(mean, std_deviation),0.9 * mean, 1.1* mean)); # travel time chosen from truncated normal distribution
                        transition_matrix[i,t,3] = travel_time;

                        # sampling travel distances
                        mean = distance_matrix_km[origin, destination]; # constant zone distances as mean
                        std_deviation = 0.1 * mean; # arbitrary chosen standard deviation as 10% of the mean value, to create some variation.
                        travel_distance = rand(Truncated(Normal(mean, std_deviation),0.9 * mean, 1.1 * mean)); # travel distance chosen from truncated normal distribution.
                        transition_matrix[i,t,4] = travel_distance;
                    end
                end
            end

            # state update
            if(t<T)
                state_matrix[:,t+1] = transition_matrix[:,t,2];
            end
            println(trunc(Int,t/0.23),"% of the dataset is simulated")
        end

        return state_matrix, transition_matrix

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
correct = ( pi/4 + atan(3.0) ) / 2
for endpt in [neither, both, left, right]
    x, w = legendre(40, endpt)
    if endpt in [left, both]
        @test x[1] == -1
    end
    if endpt in [right, both]
        @test x[end] == 1
    end
    Q = sum(w .* 1 ./ ( 1 .+ (2x .- 1).^2))
    @test Q ≈ correct
end


text:
#Test DataRS manipulations
#-------------------------------------------------------------------------------

using MDDatasets
using Test

#==Input data
===============================================================================#
#get_ydata: Generates sine data from input parameters
get_ydata(t, tbit, vdd, trise) = sin(2pi*t/tbit)*(trise/tbit)+vdd

t = DataF1((1:2)*1e-9)

#Construct DataRS object from parametric sweeps:
data = fill(DataRS, PSweep("tbit", [1, 3] * 1e-9)) do tbit
	#Need to explicitly specify DataRS{DataF1} as type for inner-most sweep:
	fill(DataRS{DataF1}, PSweep("VDD", 0.9 * [0.9, 1])) do vdd
		trise = 0.1*tbit
		return get_ydata(t, tbit, vdd, trise)
	end
end

#==Tests
===============================================================================#
@testset "Construction of DataRS objects" begin
	#TODO: write tests to verify structure
	@info("Visually confirm structure of \"data\":")
	@show data

	@test data.sweep.id == "tbit" #should be outer-most loop
	@test data.elem[1].sweep.id == "VDD" #should be inner-most loop
end

@testset "Access parameter data" begin
	#Extract value of tbit (first sweep):
	data_tbit = parameter(data, "tbit")

	#Get sweep values:
	sw_tbit = data.sweep.v
	sw_vdd = data.elem[1].sweep.v

	for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
		@test data_tbit.elem[itbit].elem[ivdd] == sw_tbit[itbit]
	end
end

@testset "Broacasting operations on DataRS" begin
	#Get sweep values:
	sw_tbit = data.sweep.v
	sw_vdd = data.elem[1].sweep.v

	_v = 0*data+5 #Same dimensionality as data
	@testset "_v = 0*data+5" begin
		v5 = zeros(length(t)) .+ 5

		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _v.elem[itbit].elem[ivdd]
			@test _elem.x == t.x
			@test _elem.y == v5
		end
	end

	_vshift = data + _v
	@testset "_vshift = data + _v" begin
		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _vshift.elem[itbit].elem[ivdd]
			@test _elem.x == t.x
			@test _elem.y == data.elem[itbit].elem[ivdd].y .+ 5
		end
	end

	_m1 = maximum(data)
	@testset "_m1 = maximum(data)" begin
		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _m1.elem[itbit].elem[ivdd]
			@test typeof(_elem) == Float64
			@test _elem == maximum(data.elem[itbit].elem[ivdd].y)
		end
	end

	_delta = data - _m1
	@testset "_delta = data - _m1" begin
		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _delta.elem[itbit].elem[ivdd]
			@test _elem.x == t.x
			_y = data.elem[itbit].elem[ivdd].y
			@test _elem.y == _y .- maximum(_y)
		end
	end

	_m2 = maximum(_m1) #Collapse inner-most dimension again
	@testset "_m2 = maximum(_m1)" begin
		for itbit in 1:length(sw_tbit)
			_elem = _m2.elem[itbit]
			@test typeof(_elem) == Float64
			@test _elem == maximum(_m1.elem[itbit])
		end
	end

	_m3 = maximum(_m2) #Collapse inner-most dimension again
	@testset "_m3 = maximum(_m2)" begin
		_elem = _m3
		m2val = [_m2.elem[i] for i in 1:length(_m2)] #collect values in m2
		@test typeof(_elem) == Float64
		@test _elem == maximum(m2val)
	end

end

:Test_Complete


text:
"""
    geometric_factors(x, y, Dr, Ds)
    geometric_factors(x, y, z, Dr, Ds, Dt, Filters=(I,I,I))

Compute metrics of mappings between "real" elements and reference elements,
outward pointing normals on faces of every elements, and Jacobian.

x,y,z are arrays of coordinates, and Dr, Ds, Dt are nodal differentiation matrices
Filters = tuple of filtering matrices (e.g., to reduce degree in r,s, and t for GCL)

Geometric terms in 3D are constructed to ensure satisfaction of free-stream preservation
using the curl-based construction from 'Metric identities and the DG-SEM on curvilinear
meshes' (Kopriva 2006).
"""

function geometric_factors(x, y, Dr, Ds)
    xr, xs = Dr*x, Ds*x
    yr, ys = Dr*y, Ds*y

    J = @. -xs * yr + xr * ys
    rxJ, sxJ =  ys, -yr
    ryJ, syJ = -xs,  xr

    return rxJ, sxJ, ryJ, syJ, J
end

function geometric_factors(x, y, z, Dr, Ds, Dt, Filters=(I,I,I))

    xr, xs, xt = Dr*x, Ds*x, Dt*x
    yr, ys, yt = Dr*y, Ds*y, Dt*y
    zr, zs, zt = Dr*z, Ds*z, Dt*z

    Fr, Fs, Ft = (Dr*y).*z, (Ds*y).*z, (Dt*y).*z
    Fr,Fs,Ft = ((A,x)->A*x).(Filters,(Fr,Fs,Ft))
    rxJ = Dt*(Fs) - Ds*(Ft)
    sxJ = Dr*(Ft) - Dt*(Fr)
    txJ = Ds*(Fr) - Dr*(Fs)

    Fr, Fs, Ft = (Dr*x).*z, (Ds*x).*z, (Dt*x).*z
    Fr,Fs,Ft = ((A,x)->A*x).(Filters,(Fr,Fs,Ft))
    ryJ = -(Dt*(Fs) - Ds*(Ft))
    syJ = -(Dr*(Ft) - Dt*(Fr))
    tyJ = -(Ds*(Fr) - Dr*(Fs))

    Fr, Fs, Ft = (Dr*y).*x, (Ds*y).*x, (Dt*y).*x
    Fr,Fs,Ft = ((A,x)->A*x).(Filters,(Fr,Fs,Ft))
    rzJ = -(Dt*(Fs) - Ds*(Ft))
    szJ = -(Dr*(Ft) - Dt*(Fr))
    tzJ = -(Ds*(Fr) - Dr*(Fs))

    J = @. xr * (ys * zt - zs * yt) - yr * (xs * zt - zs * xt) + zr * (xs * yt - ys * xt)

    return rxJ, sxJ, txJ, ryJ, syJ, tyJ, rzJ, szJ, tzJ, J
end

"""
    estimate_h(rd::RefElemData, md::MeshData)
    estimate_h(e, rd::RefElemData, md::MeshData) # e = element index

Estimates the mesh size via min size_of_domain * |J|/|sJ|, since |J| = O(hᵈ) and |sJ| = O(hᵈ⁻¹).
"""
function estimate_h(rd::RefElemData{DIM}, md::MeshData{DIM}) where {DIM}
    hmin = Inf
    for e in 1:md.num_elements
        h_e = estimate_h(e, rd, md)
        hmin = min(hmin, h_e)
    end
    return hmin * compute_domain_size(rd, md)^(1/DIM)
end

function estimate_h(e, rd::RefElemData, md::MeshData)
    sJ_e = reshape(view(md.sJ, :, e), rd.Nfq ÷ rd.Nfaces, rd.Nfaces)
    sJ_face = zero(eltype(md.sJ))
    for f in 1:rd.Nfaces
        sJ_face = max(sJ_face, minimum(view(sJ_e, :, f)) / face_scaling(rd, f))
    end
    h_e = minimum(view(md.J, :, e)) / sJ_face
    return h_e
end

face_scaling(rd, f) = 1.0
face_scaling(rd::RefElemData{2, Tri}, f) = f == 3 ? sqrt(2) : 1.0 # sJ incorporates length of long triangle edge
face_scaling(rd::RefElemData{3, Tet}, f) = f == 2 ? sqrt(3) : 1.0 # sJ incorporates area of larger triangle face
compute_domain_size(rd::RefElemData, md::MeshData) = sum(rd.M * md.J)


text:
using Test
using CliOptions

include("testutils.jl")

@testset "AbstractOptionGroup()" begin
    @testset "show(::OptionGroup); $(join(v[1],','))" for v in [
        (["-a"], "OptionGroup(Option(:a))"),
        (["-a", "-b"], "OptionGroup(Option(:a),Option(:b))"),
    ]
        names, expected_repr = v
        result = CliOptions.ParseResult()
        option = OptionGroup([Option(name) for name in names]...)
        @test repr(option) == expected_repr
    end

    @testset "show(::MutexGroup); $(join(v[1],','))" for v in [
        (["-a"], "MutexGroup(Option(:a))"),
        (["-a", "-b"], "MutexGroup(Option(:a),Option(:b))"),
    ]
        names, expected_repr = v
        result = CliOptions.ParseResult()
        option = MutexGroup([Option(name) for name in names]...)
        @test repr(option) == expected_repr
    end

    option0 = OptionGroup()
    option1 = OptionGroup(Option("-a"))
    option2 = OptionGroup(Option("-a"), Option("-b"))
    @testset "length(::OptionGroup)" begin
        @test length(option0) == 0
        @test length(option1) == 1
        @test length(option2) == 2
    end
    @testset "iterate(::OptionGroup)" begin
        @test [o.names[1] for o in option0] == String[]
        @test [o.names[1] for o in option1] == String["-a"]
        @test [o.names[1] for o in option2] == String["-a", "-b"]
    end

    option0 = MutexGroup()
    option1 = MutexGroup(Option("-a"))
    option2 = MutexGroup(Option("-a"), Option("-b"))
    @testset "length(::MutexGroup)" begin
        @test length(option0) == 0
        @test length(option1) == 1
        @test length(option2) == 2
    end
    @testset "iterate(::MutexGroup)" begin
        @test [o.names[1] for o in option0] == String[]
        @test [o.names[1] for o in option1] == String["-a"]
        @test [o.names[1] for o in option2] == String["-a", "-b"]
    end
end