text:
"""
`module  JAC.AutoIonization`
    ... a submodel of JAC that contains all methods for computing Auger properties between some initial and final-state
        multiplets.
"""
module AutoIonization

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

    """
    `struct  Settings`  ... defines a type for the details and parameters of computing Auger lines.

        + calcAnisotropy          ::Bool                         ... True, if the intrinsic alpha_2,4 angular parameters are to be calculated,
                                                                     and false otherwise.
        + printBefore             ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + lineSelection           ::LineSelection                ... Specifies the selected levels, if any.
        + minAugerEnergy          ::Float64                      ... Minimum energy of free (Auger) electrons to be included.
        + maxAugerEnergy          ::Float64                      ... Maximum energy of free (Auger) electrons to be included.
        + maxKappa                ::Int64                        ... Maximum kappa value of partial waves to be included.
        + operator                ::AbstractEeInteraction        ... Auger operator that is to be used for evaluating the Auger amplitudes:
                                                                     allowed values are: CoulombInteraction(), BreitInteraction(), ...
    """
    struct Settings
        calcAnisotropy            ::Bool
        printBefore               ::Bool
        lineSelection             ::LineSelection
        minAugerEnergy            ::Float64
        maxAugerEnergy            ::Float64
        maxKappa                  ::Int64
        operator                  ::AbstractEeInteraction
    end

    """
    `AutoIonization.Settings()`  ... constructor for the default values of AutoIonization line computations
    """
    function Settings()
        Settings(false, false, LineSelection(), 0., 10e5, 100, CoulombInteraction())
    end

    """
    `AutoIonization.Settings(set::AutoIonization.Settings;`

            calcAnisotropy=..,      printBefore=..,
            minAugerEnergy=..,      maxAugerEnergy=..,          maxKappa=..,            operator=..)

        ... constructor for modifying the given AutoIonization.Settings by 'overwriting' the previously selected parameters.
    """
    function Settings(set::AutoIonization.Settings;
        calcAnisotropy::Union{Nothing,Bool}=nothing,            printBefore::Union{Nothing,Bool}=nothing,
        lineSelection::Union{Nothing,LineSelection}=nothing,
        minAugerEnergy::Union{Nothing,Float64}=nothing,         maxAugerEnergy::Union{Nothing,Float64}=nothing,
        maxKappa::Union{Nothing,Int64}=nothing,                 operator::Union{Nothing,String}=nothing)

        if  calcAnisotropy  == nothing   calcAnisotropyx  = set.calcAnisotropy    else  calcAnisotropyx  = calcAnisotropy   end
        if  printBefore     == nothing   printBeforex     = set.printBefore       else  printBeforex     = printBefore      end
        if  lineSelection   == nothing   lineSelectionx   = set.lineSelection     else  lineSelectionx   = lineSelection    end
        if  minAugerEnergy  == nothing   minAugerEnergyx  = set.minAugerEnergy    else  minAugerEnergyx  = minAugerEnergy   end
        if  maxAugerEnergy  == nothing   maxAugerEnergyx  = set.maxAugerEnergy    else  maxAugerEnergyx  = maxAugerEnergy   end
        if  maxKappa        == nothing   maxKappax        = set.maxKappa          else  maxKappax        = maxKappa         end
        if  operator        == nothing   operatorx        = set.operator          else  operatorx        = operator         end

        Settings( calcAnisotropyx, printBeforex, lineSelectionx, minAugerEnergyx, maxAugerEnergyx, maxKappax, operatorx)
    end

    # `Base.show(io::IO, settings::AutoIonization.Settings)`  ... prepares a proper printout of the variable settings::AutoIonization.Settings.
    function Base.show(io::IO, settings::AutoIonization.Settings)
        println(io, "calcAnisotropy:                $(settings.calcAnisotropy)  ")
        println(io, "printBefore:                   $(settings.printBefore)  ")
        println(io, "lineSelection:                 $(settings.lineSelection)  ")
        println(io, "minAugerEnergy:                $(settings.minAugerEnergy)  ")
        println(io, "maxAugerEnergy:                $(settings.maxAugerEnergy)  ")
        println(io, "maxKappa:                      $(settings.maxKappa)  ")
        println(io, "operator:                      $(settings.operator)  ")
    end

    """
    `struct  Channel`
        ... defines a type for a AutoIonization channel to help characterize a scattering (continuum) state of many
            electron-states with a single free electron.

        + kappa          ::Int64                ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... Auger amplitude associated with the given channel.
    """
    struct  Channel
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  Line`
        ... defines a type for a AutoIonization line that may include the definition of sublines and their
            corresponding amplitudes.

        + initialLevel   ::Level           ... initial-(state) level
        + finalLevel     ::Level           ... final-(state) level
        + electronEnergy ::Float64         ... Energy of the (incoming free) electron.
        + totalRate      ::Float64         ... Total rate of this line.
        + angularAlpha   ::Float64         ... Angular alpha_2 coefficient.
        + hasChannels    ::Bool
            ... Determines whether the individual scattering (sub-) channels are defined in terms of their free-electron energy, kappa
                and the total angular momentum/parity as well as the amplitude, or not.
        + channels       ::Array{AutoIonization.Channel,1}  ... List of AutoIonization channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        electronEnergy   ::Float64
        totalRate        ::Float64
        angularAlpha     ::Float64
        hasChannels      ::Bool
        channels         ::Array{AutoIonization.Channel,1}
    end

    """
    `AutoIonization.Line(initialLevel::Level, finalLevel::Level, totalRate::Float64)`
        ... constructor for an AutoIonization line between a specified initial and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level, totalRate::Float64)
        Line(initialLevel, finalLevel, 0., totalRate, 0., false, AutoIonization.Channel[])
    end

    # `Base.show(io::IO, line::AutoIonization.Line)`  ... prepares a proper printout of the variable line::AutoIonization.Line.
    function Base.show(io::IO, line::AutoIonization.Line)
        println(io, "initialLevel:           $(line.initialLevel)  ")
        println(io, "finalLevel:             $(line.finalLevel)  ")
        println(io, "electronEnergy:         $(line.electronEnergy)  ")
        println(io, "totalRate:              $(line.totalRate)  ")
        println(io, "angularAlpha:           $(line.angularAlpha)  ")
        println(io, "hasChannels:            $(line.hasChannels)  ")
        println(io, "channels:               $(line.channels)  ")
    end

    """
    `AutoIonization.amplitude(kind::String, channel::AutoIonization.Channel, continuumLevel::Level, initialLevel::Level,
                              grid::Radial.Grid; printout::Bool=true)`
        ... to compute the kind in  CoulombInteraction(), BreitInteraction(), CoulombBreit()   Auger amplitude
            <(alpha_f J_f, kappa) J_i || O^(Auger, kind) || alpha_i J_i>  due to the interelectronic interaction for the given
            final and initial level. A value::ComplexF64 is returned.
    """
    function amplitude(kind::AbstractEeInteraction, channel::AutoIonization.Channel, continuumLevel::Level, initialLevel::Level, grid::Radial.Grid;
                       printout::Bool=true)
        nt = length(continuumLevel.basis.csfs);    ni = length(initialLevel.basis.csfs);    partial = Subshell(9,channel.kappa)
        if  printout  printstyled("Compute ($kind) Auger matrix of dimension $nt x $ni in the continuum- and initial-state bases " *
                                  "for the transition [$(initialLevel.index)- ...] and for partial wave $(string(partial)[2:end]) ... ",
                                  color=:light_green)    end
        matrix = zeros(ComplexF64, nt, ni)
        #
        if      kind in [ CoulombInteraction(), BreitInteraction(), CoulombBreit()]        ## pure V^Coulomb interaction
        #--------------------------------------------------------------------------
            for  r = 1:nt
                for  s = 1:ni
                    if  initialLevel.basis.csfs[s].J != initialLevel.J  ||  initialLevel.basis.csfs[s].parity != initialLevel.parity      continue    end
                    ##x wa = compute("angular coefficients: e-e, Ratip2013", continuumLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                    # Calculate the spin-angular coefficients
                    if  Defaults.saRatip()
                        waR = compute("angular coefficients: e-e, Ratip2013", continuumLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                        wa  = waR
                    end
                    if  Defaults.saGG()
                        subshellList = continuumLevel.basis.subshells
                        opa  = SpinAngular.TwoParticleOperator(0, plus, true)
                        waG2 = SpinAngular.computeCoefficients(opa, continuumLevel.basis.csfs[r], initialLevel.basis.csfs[s], subshellList)
                        wa   = [1.0, waG2]
                    end
                    if  Defaults.saRatip() && Defaults.saGG() && true
                        if  length(waR[2]) != 0     println("\n>> Angular coeffients from Ratip2013   = $(waR[2]) ")    end
                        if  length(waG2)   != 0     println(  ">> Angular coeffients from SpinAngular = $waG2 ")        end
                    end
                    #
                    me = 0.
                    for  coeff in wa[2]
                        if   kind in [ CoulombInteraction(), CoulombBreit()]
                            me = me + coeff.V * InteractionStrength.XL_Coulomb(coeff.nu,
                                                    continuumLevel.basis.orbitals[coeff.a], continuumLevel.basis.orbitals[coeff.b],
                                                    initialLevel.basis.orbitals[coeff.c],   initialLevel.basis.orbitals[coeff.d], grid)   end
                        if   kind in [ BreitInteraction(), CoulombBreit()]
                            me = me + coeff.V * InteractionStrength.XL_Breit(coeff.nu,
                                                    continuumLevel.basis.orbitals[coeff.a], continuumLevel.basis.orbitals[coeff.b],
                                                    initialLevel.basis.orbitals[coeff.c],   initialLevel.basis.orbitals[coeff.d], grid)   end
                    end
                    matrix[r,s] = me
                end
            end
            if  printout  printstyled("done. \n", color=:light_green)    end
            amplitude = transpose(continuumLevel.mc) * matrix * initialLevel.mc
            amplitude = im^Basics.subshell_l(Subshell(101, channel.kappa)) * exp( -im*channel.phase ) * amplitude
            #
            #
         elseif  kind == "H-E"
        #--------------------
            iLevel = finalLevel;   fLevel = initialLevel
            amplitude = 0.;    error("stop a")
        else    error("stop b")
        end

        return( amplitude )
    end

    """
    `AutoIonization.channelAmplitude(kind::String, channel::AutoIonization.Channel, energy::Float64, finalLevel::Level,
                                     initialLevel::Level, grid::Radial.Grid)`
        ... to compute the kind = (CoulombInteraction(), BreitInteraction(), CoulombBreit())   Auger amplitude
            <(alpha_f J_f, kappa) J_i || O^(Auger, kind) || alpha_i J_i>  due to the interelectronic interaction for the given final and
            initial level. A newChannel::AutoIonization.Channel is returned.
    """
    function channelAmplitude(kind::String, channel::AutoIonization.Channel, energy::Float64, finalLevel::Level,
                              initialLevel::Level, grid::Radial.Grid)
        newiLevel = Basics.generateLevelWithSymmetryReducedBasis(initialLevel)
        newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newiLevel)
        newfLevel = Basics.generateLevelWithSymmetryReducedBasis(finalLevel)
        cOrbital, phase  = Continuum.generateOrbital(energy, Subshell(101, channel.kappa), newfLevel, grid, contSettings)
        newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newfLevel)
        newChannel = AutoIonization.Channel(channel.kappa, channel.symmetry, phase, 0.)
        amplitude = AutoIonization.amplitude(kind, channel, newcLevel, newiLevel, grid)

        newChannel = AutoIonization.Channel(newChannel.kappa, newChannel.symmetry, newChannel.phase, amplitude)
        return( newChannel )
    end

    """
    `AutoIonization.computeAmplitudesProperties(line::AutoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                                settings::AutoIonization.Settings; printout::Bool=true)`
        ... to compute all amplitudes and properties of the given line; a line::AutoIonization.Line is returned for which the amplitudes
            and properties are now evaluated.
    """
    function computeAmplitudesProperties(line::AutoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64, settings::AutoIonization.Settings;
                                         printout::Bool=true)
        newChannels = AutoIonization.Channel[];   contSettings = Continuum.Settings(false, nrContinuum);   rate = 0.
        for channel in line.channels
            newiLevel = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel, line.initialLevel.basis.subshells)
            newfLevel = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel, newiLevel.basis.subshells)
            newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newiLevel)
            cOrbital, phase  = Continuum.generateOrbitalForLevel(line.electronEnergy, Subshell(101, channel.kappa), newfLevel, nm, grid, contSettings)
            newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newfLevel)
            newChannel = AutoIonization.Channel(channel.kappa, channel.symmetry, phase, 0.)
            amplitude  = AutoIonization.amplitude(settings.operator, newChannel, newcLevel, newiLevel, grid, printout=printout)
            rate       = rate + conj(amplitude) * amplitude
            push!( newChannels, AutoIonization.Channel(newChannel.kappa, newChannel.symmetry, newChannel.phase, amplitude) )
        end
        totalRate = 2pi* rate;   angularAlpha = 0.
        newLine   = AutoIonization.Line(line.initialLevel, line.finalLevel, line.electronEnergy, totalRate, angularAlpha, true, newChannels)
        #
        if  settings.calcAnisotropy    angularAlpha = AutoIonization.computeIntrinsicAlpha(2, newLine)
            newLine   = AutoIonization.Line(line.initialLevel, line.finalLevel, line.electronEnergy, totalRate, real(angularAlpha), true, newChannels)
        end

        return( newLine )
    end

    """
    `AutoIonization.computeAmplitudesPropertiesPlasma(line::AutoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                                      settings::PlasmaShift.AugerSettings; printout::Bool=true)`
        ... to compute all amplitudes and properties of the given line but for the given plasma model; a line::AutoIonization.Line is returned
            for which the amplitudes and properties are now evaluated.
    """
    function computeAmplitudesPropertiesPlasma(line::AutoIonization.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                               settings::PlasmaShift.AugerSettings; printout::Bool=true)
        newChannels = AutoIonization.Channel[];   contSettings = Continuum.Settings(false, nrContinuum);   rate = 0.
        for channel in line.channels
            newiLevel = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel)
            newiLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newiLevel)
            newfLevel = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel)
            @warn "Adapt a proper continuum orbital for the plasma potential"
            cOrbital, phase  = Continuum.generateOrbitalForLevel(line.electronEnergy, Subshell(101, channel.kappa), newfLevel, nm, grid, contSettings)
            newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newfLevel)
            newChannel = AutoIonization.Channel(channel.kappa, channel.symmetry, phase, 0.)
            @warn "Adapt a proper Auger amplitude for the plasma e-e interaction"
            amplitude = 1.0
            # amplitude  = AutoIonization.amplitude(settings.operator, newChannel, newcLevel, newiLevel, grid)
            rate       = rate + conj(amplitude) * amplitude
            push!( newChannels, AutoIonization.Channel(newChannel.kappa, newChannel.symmetry, newChannel.phase, amplitude) )
        end
        totalRate = 2pi* rate;   angularAlpha = 0.
        newLine   = AutoIonization.Line(line.initialLevel, line.finalLevel, line.electronEnergy, totalRate, angularAlpha, true, newChannels)

        return( newLine )
    end

    """
    `AutoIonization.computeIntrinsicAlpha(k::Int64, line::AutoIonization.Line)`
        ... to compute the intrinsic alpha_k anisotropy parameter for the given line. A value::Float64 is returned.
    """
    function  computeIntrinsicAlpha(k::Int64, line::AutoIonization.Line)
        if  !line.hasChannels   error("No channels are defined for the given AutoIonization.line.")                   end
        wn = 0.;    for  channel in line.channels    wn = wn + conj(channel.amplitude) * channel.amplitude   end
        wa = 0.;    Ji = line.initialLevel.J;    Jf = line.finalLevel.J;
        for  cha  in line.channels
            j = AngularMomentum.kappa_j(cha.kappa);    l = AngularMomentum.kappa_l(cha.kappa)
            for  chp  in line.channels
                jp = AngularMomentum.kappa_j(chp.kappa);    lp = AngularMomentum.kappa_l(chp.kappa)
                wa = wa + sqrt( AngularMomentum.bracket([l, lp, j, jp]) ) *
                          AngularMomentum.ClebschGordan(l, AngularM64(0), lp, AngularM64(0), AngularJ64(k), AngularM64(0)) *
                          AngularMomentum.Wigner_6j(Ji, j, Jf, jp, Ji, AngularJ64(k)) *
                          AngularMomentum.Wigner_6j(l,  j, AngularJ64(1//2), jp, lp, AngularJ64(k)) *
                          cha.amplitude * conj(chp.amplitude)
            end
        end
        value = AngularMomentum.phaseFactor([Ji, +1, Jf, +1, AngularJ64(k), -1, AngularJ64(1//2)]) *
                sqrt(Basics.twice(Ji) + 1) * wa / wn

        if  false
            # Calculate the value given by M.H. Chen, PRA 47 (1993) 3733; this does not agree exactly so far.
            wa = 0.
            for  cha  in line.channels
                j = AngularMomentum.kappa_j(cha.kappa);    l = AngularMomentum.kappa_l(cha.kappa);    phase = cha.phase
                for  chp  in line.channels
                    jp = AngularMomentum.kappa_j(chp.kappa);    lp = AngularMomentum.kappa_l(chp.kappa);    phasep = chp.phase
                    wa = wa + AngularMomentum.phaseFactor([Ji, +1, Jf, -1, AngularJ64(1//2)]) * im^(Float64(l) - Float64(lp)) * cos(phase - phasep) *
                              sqrt( AngularMomentum.bracket([l, lp, j, jp, AngularJ64(k), Ji]) ) *
                              AngularMomentum.Wigner_3j(lp, l, AngularJ64(k), AngularJ64(0), AngularJ64(0), AngularJ64(0)) *
                              AngularMomentum.Wigner_6j(j, jp, AngularJ64(k), lp, l, AngularJ64(1//2)) *
                              AngularMomentum.Wigner_6j(Ji, Ji, AngularJ64(k), jp, j, Jf) *
                              cha.amplitude * conj(chp.amplitude)
                end
            end
            println("*** Comparison value = $value    Chen-value = $(wa/wn)")
        end

        return( value )
    end

    """
    `AutoIonization.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                 settings::AutoIonization.Settings; output=true, printout::Bool=true)`
        ... to compute the Auger transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{AutoIonization.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                           settings::AutoIonization.Settings; output=true, printout::Bool=true)
        println("")
        printstyled("AutoIonization.computeLines(): The computation of Auger rates and properties starts now ... \n", color=:light_green)
        printstyled("------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        lines = AutoIonization.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    AutoIonization.displayLines(lines)    end
        # Determine maximum energy and check for consistency of the grid
        maxEnergy = 0.;   for  line in lines   maxEnergy = max(maxEnergy, line.electronEnergy)   end
        nrContinuum = Continuum.gridConsistency(maxEnergy, grid)
        # Calculate all amplitudes and requested properties
        newLines = AutoIonization.Line[]
        for  line in lines
            newLine = AutoIonization.computeAmplitudesProperties(line, nm, grid, nrContinuum, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        AutoIonization.displayRates(stdout, newLines, settings)
        AutoIonization.displayLifetimes(stdout, newLines)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   AutoIonization.displayRates(iostream, newLines, settings);   AutoIonization.displayLifetimes(iostream, newLines)     end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `AutoIonization.computeLinesCascade(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                        settings::AutoIonization.Settings; output::Bool=true, printout::Bool=true)`
        ... to compute the Auger transition amplitudes and all properties as requested by the given settings. The computations
            and printout is adapted for large cascade computations by including only lines with at least one channel and by sending
            all printout to a summary file only. A list of lines::Array{AutoIonization.Lines} is returned.
    """
    function  computeLinesCascade(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                  settings::AutoIonization.Settings; output::Bool=true, printout::Bool=true)

        lines = AutoIonization.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        # if  settings.printBefore    AutoIonization.displayLines(lines)    end
        # Determine maximum energy and check for consistency of the grid
        maxEnergy = 0.;   for  line in lines   maxEnergy = max(maxEnergy, line.electronEnergy)   end
        nrContinuum = Continuum.gridConsistency(maxEnergy, grid)
        # Calculate all amplitudes and requested properties
        newLines = AutoIonization.Line[]
        for  (i,line)  in  enumerate(lines)
            if  rem(i,10) == 0    println("> Auger line $i:  ... calculated ")    end
            newLine = AutoIonization.computeAmplitudesProperties(line, nm, grid, nrContinuum, settings, printout=printout)
            push!( newLines, newLine)
        end
        # Print all results to a summary file, if requested
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   AutoIonization.displayRates(iostream, newLines, settings)     end
        #
        if    output    return( newLines )
        else            return( nothing )
        end
    end

    """
    `AutoIonization.computeLinesFromOrbitals(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                             settings::AutoIonization.Settings, contOrbitals::Dict{Subshell, Orbital}; output::Bool=true, printout::Bool=true)`
        ... to compute the Auger transition amplitudes and all properties as requested by the given settings but by using the given set of
            continuum orbitals. The computations and printout is adapted for large cascade computations by including only lines with at least
            one channel and by sending all printout to a summary file only. A list of lines::Array{AutoIonization.Lines} is returned.
    """
    function  computeLinesFromOrbitals(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                       settings::AutoIonization.Settings, contOrbitals::Dict{Subshell, Orbital}; output::Bool=true, printout::Bool=true)

        lines = AutoIonization.determineLines(finalMultiplet, initialMultiplet, settings)
        # Calculate all amplitudes and requested properties
        newLines = AutoIonization.Line[]
        for  (i,line)  in  enumerate(lines)
            if  rem(i,500) == 0    println("> Auger line $i:  ... calculated ")    end
            # Calculate the individual channels with the given orbitals
            newChannels = AutoIonization.Channel[];   rate = 0.
            for channel in line.channels
                newiLevel  = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel, line.initialLevel.basis.subshells)
                newfLevel  = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel, newiLevel.basis.subshells)
                sh         = Subshell(101, channel.kappa)
                if haskey(contOrbitals, sh)   cOrbital = contOrbitals[sh]      else    println(">>> skip Auger channel for $sh");   continue    end
                newiLevel  = Basics.generateLevelWithExtraSubshell(sh, newiLevel)
                newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newfLevel)
                newChannel = AutoIonization.Channel(channel.kappa, channel.symmetry, 0., 0.)
                amplitude  = AutoIonization.amplitude(settings.operator, newChannel, newcLevel, newiLevel, grid, printout=printout)
                rate       = rate + conj(amplitude) * amplitude
                push!( newChannels, AutoIonization.Channel(newChannel.kappa, newChannel.symmetry, newChannel.phase, amplitude) )
            end
            totalRate = 2pi* rate;   angularAlpha = 0.
            newLine   = AutoIonization.Line(line.initialLevel, line.finalLevel, line.electronEnergy, totalRate, angularAlpha, true, newChannels)
            ##
            ## if  settings.calcAnisotropy    angularAlpha = AutoIonization.computeIntrinsicAlpha(2, newLine)
            ##  newLine   = AutoIonization.Line(line.initialLevel, line.finalLevel, line.electronEnergy, totalRate, real(angularAlpha), true, newChannels)
            ## end
            push!( newLines, newLine)
        end
        # Print all results to a summary file, if requested
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   AutoIonization.displayRates(iostream, newLines, settings)     end
        #
        if    output    return( newLines )
        else            return( nothing )
        end
    end

    """
    `AutoIonization.computeLinesPlasma(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                       settings::PlasmaShift.AugerSettings; output=true)`
        ... to compute the Auger transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{AutoIonization.Lines} is returned.
    """
    function  computeLinesPlasma(finalMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid,
                                 settings::PlasmaShift.AugerSettings; output=true)
        println("")
        printstyled("AutoIonization.computeLinesPlasma(): The computation of Auger rates starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        augerSettings = AutoIonization.Settings(false, settings.printBefore, settings.selectLines, settings.selectedLines, 0., 1.0e6, 100, CoulombInteraction())
        lines = AutoIonization.determineLines(finalMultiplet, initialMultiplet, augerSettings)
        # Display all selected lines before the computations start
        if  settings.printBefore    AutoIonization.displayLines(lines)    end
        # Determine maximum energy and check for consistency of the grid
        maxEnergy = 0.;   for  line in lines   maxEnergy = max(maxEnergy, line.electronEnergy)   end
        nrContinuum = Continuum.gridConsistency(maxEnergy, grid)
        # Calculate all amplitudes and requested properties
        newLines = AutoIonization.Line[]
        for  line in lines
            newLine = AutoIonization.computeAmplitudesPropertiesPlasma(line, nm, grid, nrContinuum, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        AutoIonization.displayRates(stdout, newLines, augerSettings)
        AutoIonization.displayLifetimes(stdout, newLines)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   AutoIonization.displayRates(iostream, newLines, augerSettings);   AutoIonization.displayLifetimes(iostream, newLines)     end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `AutoIonization.determineChannels(finalLevel::Level, initialLevel::Level, settings::AutoIonization.Settings)`
        ... to determine a list of Auger Channel for a transitions from the initial to final level and by taking into account the particular
            settings of for this computation; an Array{AutoIonization.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::AutoIonization.Settings)
        channels  = AutoIonization.Channel[];
        symi      = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        kappaList = AngularMomentum.allowedKappaSymmetries(symi, symf)
        for  kappa in kappaList
            if  abs(kappa) > settings.maxKappa      continue    end
            push!(channels, AutoIonization.Channel(kappa, symi, 0., Complex(0.)) )
        end
        return( channels )
    end

    """
    `AutoIonization.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::AutoIonization.Settings)`
        ... to determine a list of AutoIonization.Line's for transitions between levels from the initial- and final-state multiplets, and
            by taking into account the particular selections and settings for this computation; an Array{AutoIonization.Line,1} is returned.
            Apart from the level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::AutoIonization.Settings)
        lines = AutoIonization.Line[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    energy = iLevel.energy - fLevel.energy
                    if   energy < 0.01                                                             continue   end
                    if   energy < settings.minAugerEnergy  ||  energy > settings.maxAugerEnergy    continue   end
                    channels = AutoIonization.determineChannels(fLevel, iLevel, settings)
                    push!( lines, AutoIonization.Line(iLevel, fLevel, energy, 0., 0., true, channels) )
                end
            end
        end
        return( lines )
    end

    """
    `AutoIonization.displayLines(lines::Array{AutoIonization.Line,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of all selected
            transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{AutoIonization.Line,1})
        nx = 150
        println(" ")
        println("  Selected Auger lines:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                                sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                                sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(14, "Energy e_A"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(37, "List of kappas and total symmetries"; na=4)
        sb = sb * TableStrings.flushleft(37, "partial (total J^P)                "; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy))  * "    "
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.electronEnergy))       * "   "
            kappaSymmetryList = Tuple{Int64,LevelSymmetry}[]
            for  i in 1:length(line.channels)
                push!( kappaSymmetryList, (line.channels[i].kappa, line.channels[i].symmetry) )
            end
            sa = sa * TableStrings.kappaSymmetryTupels(80, kappaSymmetryList)
            println( sa )
        end
        println("  ", TableStrings.hLine(nx), "\n")
        #
        return( nothing )
    end

    """
    `AutoIonization.displayLifetimes(stream::IO, lines::Array{AutoIonization.Line,1})`
        ... to list all lifetimes as associated with the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayLifetimes(stream::IO, lines::Array{AutoIonization.Line,1})
        nx = 104
        println(stream, " ")
        println(stream, "  Auger lifetimes, total rates and widths:")
        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( 8, "J^P";      na=4);                           sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(12, "Lifetime"; na=4);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("time"); na=4)
        sa = sa * TableStrings.center(14, "Total rate"; na=6);
        sb = sb * TableStrings.center(14,TableStrings.inUnits("rate"); na=4)
        sa = sa * TableStrings.center(42, "Widths"; na=2);
        sb = sb * TableStrings.center(42, "  Hartrees         Kaysers           eV"; na=2)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        notYetDone = trues(1000)
        for  line in lines
            if  notYetDone[line.initialLevel.index]
                notYetDone[line.initialLevel.index] = false
                totalRate = 0.
                for  ln in lines
                    if  ln.initialLevel.index == line.initialLevel.index   totalRate = totalRate + ln.totalRate    end
                end
                sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                sa = sa * TableStrings.center(10, TableStrings.level(line.initialLevel.index); na=2)
                sa = sa * TableStrings.center( 8, string(isym); na=4)
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("time: from atomic",  1/totalRate))            * "     "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("rate: from atomic",    totalRate))            * "      "
                sa = sa * @sprintf("%.6e", totalRate)                                                           * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic to Kayser",  totalRate))  * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic to eV",      totalRate))  * "    "
                println(stream, sa)
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `AutoIonization.displayRates(stream::IO, lines::Array{AutoIonization.Line,1}, settings::AutoIonization.Settings)`
        ... to list all results, energies, rates, etc. of the selected lines. A neat table is printed but nothing is returned
            otherwise.
    """
    function  displayRates(stream::IO, lines::Array{AutoIonization.Line,1}, settings::AutoIonization.Settings)
        nx = 106
        println(stream, " ")
        if  settings.calcAnisotropy    println(stream, "  Auger rates and intrinsic angular parameters: \n")
        else                           println(stream, "  Auger rates (without angular parameters): \n")        end
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "Energy"   ; na=2);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(14, "Electron energy"   ; na=2);
        sb = sb * TableStrings.center(14,TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(14, "Auger rate"; na=2);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("rate"); na=2)
        sa = sa * TableStrings.center(15, "alpha_2"; na=2);                           sb = sb * TableStrings.hBlank(18)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy))  * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy))       * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("rate: from atomic", line.totalRate))              * "    "
            sa = sa * TableStrings.flushright(13, @sprintf("%.4e", line.angularAlpha))            * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
function read_OD(file, meta_data_file)
    # Read file per line
    s = open(file) do f
        readlines(f)
    end
    # Look for OD header
    OD_header = findall(x-> occursin("OD", x), s)[1]

    # Find empty lines
    empty_lines = findall(x->x=="", s)

    # Read OD data
    if length(empty_lines)==1
        df = CSV.read(file, skipto=OD_header+2)
    else
        ind = findfirst(x->x>OD_header, empty_lines)
        df = CSV.read(file, datarow=OD_header+2, limit=empty_lines[ind+1]-empty_lines[ind]-1, header=false)
    end
    # Rename time and temperature
    rename!(df, Dict("Column1"=>"time_[s]", "Column2"=>"temp_[C]"))
    # Make tidy dataframe
    df  = stack(df, 3:ncol(df))

    # Rename columns
    rename!(df, Dict("variable"=>"Column", "value"=>"OD"))

    # Drop missing data
    dropmissing!(df)

    # Transform time to seconds from start
    df[!, Symbol("time_[s]")] = map(x->Dates.value(Dates.Second(x - Dates.Time(0))), df[!, Symbol("time_[s]")])

    xf = XLSX.readxlsx(meta_data_file)
    meta_df = DataFrame()
    for x in XLSX.sheetnames(xf)
        arr = convert(Array{String},xf[x]["A1:H6"])
        flat_arr = vcat([arr[i, :] for i in 1:size(arr, 1)]...)
        meta_df[!, Symbol(x)]= flat_arr
    end
    meta_df[! , :Column] = [Symbol("Column$i") for i in 3:50]
    return join(meta_df, df, on = :Column, kind=:outer)
end

ext: jl


text:
test_str = "939
7,13,x,x,59,x,31,19"

test_list = split(test_str, '\n')

list = readlines("day13_input.txt")

function departure_time(list)
    earliest = parse(Int, list[1])
    buses = map(x-> parse(Int, x), split(replace(list[2], ",x" => ""), ','))
    next_departures = Dict(map(x->(div(earliest, x) + 1 ) * x => x, buses))
    departure = minimum(keys(next_departures))
    return (departure - earliest) * next_departures[departure]
end

@show departure_time(test_list)
@show departure_time(list)
@show departure_time(test_list) == 295
@show departure_time(list) == 246

function time_stamp(line)
    buses = map(x-> parse(Int, x), split(replace(line, "x" => "-1"), ','))
    imax = argmax(buses .+ Array(0:-1:1-length(buses)))
    bmax = buses[imax]
    n = 1
    t = n*bmax + 1 - imax
    len = length(buses)
    while !all(≤(0), map(i->mod(t+i-1, buses[i]), 1:len))
        n += 1
        t = n*bmax + 1 - imax
    end
    return t
end

@show time_stamp(test_list[2])

@show time_stamp("17,x,13,19") == 3417
@show time_stamp("67,7,59,61") == 754018
@show time_stamp("67,x,7,59,61") == 779210
@show time_stamp("67,7,x,59,61") == 1261476
@show time_stamp("1789,37,47,1889") == 1202161486
#@show time_stamp(list[2]) # Too slow

# Implementing algorithm for solving linear Diophantine equations as
# explained in https://www.math.uwaterloo.ca/~wgilbert/Research/GilbertPathria.pdf
# but there is some bug. It works on all test examples, but not on the problem list.

function time_stamp_diophantine(line)
    buses = map(x-> parse(Int, x), split(replace(line, "x" => "0"), ','))
    running = [b for b=buses if b > 0]
    delays = [1-i for i=1:length(buses) if buses[i] > 0]
    offset = (mod.(running+delays,running))[2:end]
    len=length(running)
    T = fill(0, len, len)
    for i=1:len
        T[i,i] = 1
    end
    A = fill(0, len-1, len)
    for i=1:len-1
        A[i,1] = running[1]
        A[i,i+1] = - running[i+1]
    end
    R = deepcopy(A')
    for i=2:len
        while R[i,i-1] ≠ 0
            if abs(R[i-1,i-1]) > abs(R[i, i-1])
                a = div(R[i-1,i-1], R[i, i-1])
                if mod(R[i-1,i-1], R[i, i-1]) == 0
                    a -= 1
                end
                R[i-1,:] .-= a * R[i, :]
                T[i-1,:] .-= a * T[i,:]
            else
                a = div(R[i, i-1], R[i-1,i-1])
                R[i,:] .-= a * R[i-1,:]
                T[i,:] .-= a * T[i-1,:]
            end
        end
    end
    k = BigInt.((R[1:end-1,:]')\offset)
    k = [k; BigInt(floor(-sum([T[i,1]*k[i]/T[len,1] for i=1:len-1])))]
    n1 = sum([T[i,1]*k[i] for i=1:len])
#    while !all(>(0), [sum([T[i,j]*k[i] for i=1:len]) for j=1:len])
    while n1 < 0
        k[len] += sign(T[len,1])
        n1 = sum([T[i,1]*k[i] for i=1:len])
    end

    t = n1*running[1]
    return running, offset, len, A, T, R, k, n1, t
end

@show time_stamp_diophantine(test_list[2])[end] == 1068781
@show time_stamp_diophantine("7,13,x,x,59,x,31,19")[end] == 1068781
@show time_stamp_diophantine("17,x,13,19")[end] == 3417
@show time_stamp_diophantine("67,7,59,61")[end]  == 754018
@show time_stamp_diophantine("67,x,7,59,61")[end]  == 779210
@show time_stamp_diophantine("67,7,x,59,61")[end]  == 1261476
@show time_stamp_diophantine("1789,37,47,1889")[end]  == 1202161486
@show time_stamp_diophantine(list[2])[end] == 939490236001473

function test_result(line)
    running, offset, len, A, T, R, k, n1, t = time_stamp_diophantine(line)

    check = all(>(0), [sum([T[i,j]*k[i] for i=1:len]) for j=1:len])
    timings = (T'k) .* running
    check_offset = timings[2:end] - fill(timings[1], len-1)
    difference = all(==(0), check_offset + offset)
    return check, difference, T*A'==R
end

# Using the algorithm for the Chinese Remainder Theorem

function time_stamp_CRG(line)
    buses = map(x-> parse(Int, x), split(replace(line, "x" => "0"), ','))
    running = [b for b=buses if b > 0]
    offset = [mod(b+1-i,b) for (i,b) in enumerate(buses) if b>0]
    N = prod(running)
    d = div.(N,running)
    result = sum(offset .* d .* invmod.(d,running))
    return mod(result, N)
end

# There is an option to use Mods.jl package, which already has an implemention
# of the Chinese Remainder Theorem, as I learned from other Humans of Julia on
# Discord, but it is more obscure and slower, I prefer the direct implementation
# above.

using Mods

function time_stamp_ModsCRT(line)
    buses = map(x-> parse(Int, x), split(replace(line, "x" => "0"), ','))
    running = [b for b=buses if b > 0]
    x_delays = [i-1 for i=1:length(buses) if buses[i] > 0]
    mod_prep = Mod[Mod(b,d) for (d,b) in zip(running, x_delays)]
    rem = Mods.CRT(mod_prep...)
    return modulus(rem) - rem.val
end

@show time_stamp_CRG(test_list[2]) == 1068781
@show time_stamp_CRG(list[2]) == 939490236001473
@show time_stamp_ModsCRT(test_list[2]) == 1068781
@show time_stamp_ModsCRT(list[2]) == 939490236001473

#=
```julia
julia> @btime time_stamp_ModsCRT(list[2])
  23.110 μs (45 allocations: 5.47 KiB)
939490236001473

julia> @btime time_stamp_CRG(list[2])
  15.109 μs (24 allocations: 5.28 KiB)
939490236001473
```
 =#

nothing

#=

If `x=buses` (with `x=-1` if bus = 'x'), time `t` is given by the condition
that `t+i-1` is an integer multiple of `x[i]`, i.e. `mod(t+i-1,x) == 0` or
`mod(t, x) = 1-i` or `mod(t,x) = x-i+1`

In case `t` is not so big, we can just search for this condition
by stepping up `t` by multiples of `x[1]`. That was my first approach,
which works fast for the test set. But for the real list, that is not
computationally feasible. So we do differently.

The conditions

    `t + i - 1 = nᵢ xᵢ,   ∀ i=1:length(x) s.th. xᵢ > 0`

are equivalent to
    `t = n₁x₁ = n₂x₂ - 1 = ⋯ = nᵢxᵢ - i + 1,     ∀ i, with xᵢ > 0`

And we know that `x₁ > 0`. Hence we just iterate over `i` to find

    `For i = 1:
        n₁¹ = 1
    For i = 2:
        n₂² such that n₂²x₂ - 1 = jn₁¹x₁
        n₁² = j₂²n₁¹
    For i = 3:
        n₃³ such that n₃³x₂ - 2 = jn₁²x₁
        n₁² = jn₁²
    `n₁ such that t = n₁ x₁`
    `n₂ such that n₂x₂ - 1 = `

```julia
n = 1
for i = 2:length(x)
    j = 1
    while mod(j*n*buses[i]-i+1, buses[1]) > 0
        j += 1
    end
    n *= j
end
```

But this is slow in some cases. So we notice the problem can also be written as
    `t = 0 (mod x₁); t = x₂-1 (mod x₂); … t = xᵢ+1-i (mod xᵢ); …`

which is in the context of the Chinese Remainder Problem (see e.g.
[Chinese Remainder Theorem](https://en.wikipedia.org/wiki/Chinese_remainder_theorem))
and [Chinese Remainder Theorem with algorithm](https://shainer.github.io/crypto/math/2017/10/22/chinese-remainder-theorem.html))

An algorithm, in python, from the latter homepage, is the following

def ChineseRemainderGauss(x, a):
    """
    Return t such that mod(t,xᵢ) = aᵢ, ∀i=1,…,n, if 0 ≤ aᵢ< xᵢ
    and the xᵢ are pairwise coprime.
    """
    from math import prod
    P = prod(x)
    result = 0

    for i in range(len(n)):
        ai = a[i]
        xi = x[i]
        bi = P // xi

        result += ai * bi * invmod(bi, xi)

    return result % P

See also https://shainer.github.io/crypto/math/2017/10/22/chinese-remainder-theorem.html

 =#


text:
using Bigleaf
using Documenter, Latexify, DataFrames

# allow plot to work without display
# https://discourse.julialang.org/t/generation-of-documentation-fails-qt-qpa-xcb-could-not-connect-to-display/60988/2
ENV["GKSwstype"] = "100"

DocMeta.setdocmeta!(Bigleaf, :DocTestSetup, :(using Bigleaf, Latexify, DataFrames);
    recursive=true, warn=false)
doctest(Bigleaf, manual = false)

makedocs(;
    # modules=[Bigleaf], # uncomment to show warnings on non-included docstrings
    authors="Thomas Wutzler <twutz@bgc-jena.mpg.de>, Jürgen Knauer <Juergen.Knauer@csiro.au> and contributors",
    repo="https://github.com/bgctw/Bigleaf.jl/blob/{commit}{path}#{line}",
    sitename="Bigleaf.jl",
    doctestfilters=[r".*Info.*"],
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/Bigleaf.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        #"Meteorological variables" => "metorological_variables.md",
        #"Unit conversions" => "unit_conversions.md",
        "Walkthrough" => "walkthrough.md",
        hide("metorological_variables.md"),
        hide("stability_correction.md"),
        hide("surface_roughness.md"),
        hide("boundary_layer_conductance.md"),
        hide("surface_conductance.md"),
        hide("aerodynamic_conductance.md"),
        hide("evapotranspiration.md"),
        hide("surface_conductance.md"),
        hide("global_radiation.md"),
        hide("unit_conversions.md"),
        hide("BigleafConstants.md"),
        hide("filter_data.md"),
        #"Index" => "autodocs.md",
        ],
)

deploydocs(;
    repo="github.com/bgctw/Bigleaf.jl",
    devbranch="main",
)


text:
using MATLAB
include("../Cell-Video/utils/data.jl")
include("../Cell-Video/utils/measure.jl")
using .Data

function run_KTHSE(; σ::Float64 = 1.0)
    X, M = gen_zoomed_data(15*4, 10, σ=1, zoom=2)
    mX = mxarray(Array{Matrix{Float64}, 1}([hcat(x...)' for x in X]))
    path = mat"RunWithoutImages($mX, $σ)"
    return acc2fβ(calc_path_accuracy_point(path, M, length.(X)))
end

function run_KTHSE(X, M; σ::Float64 = 1.0, fβ = true)
    mX = mxarray(Array{Matrix{Float64}, 1}([hcat(x...)' for x in X]))
    path = mat"RunWithoutImages($mX, $σ)"

    if size(path, 2) == 1
        # convert Array{, 1} to Array{, 2}
        if fβ
            return acc2fβ(calc_path_accuracy_point(reshape(path, :, 1), M, length.(X)))
        else
            return calc_path_accuracy_point(reshape(path, :, 1), M, length.(X))
        end
    elseif size(path, 2) == 0
        # no cell
        if fβ
            return acc2fβ(calc_path_accuracy_point(-ones(size(path, 1), 1), M, length.(X)))
        else
            return calc_path_accuracy_point(-ones(size(path, 1), 1), M, length.(X))
        end
    else
        if fβ
            return acc2fβ(calc_path_accuracy_point(path, M, length.(X)))
        else
            return calc_path_accuracy_point(path, M, length.(X))
        end
    end
end

function cum_run_KTHSE(X, M; σ::Float64 = 1.0)
    f = length(X)
    acc = ones(f)
    for t = 2:f
        mX = mxarray(Array{Matrix{Float64}, 1}([hcat(x...)' for x in X[1:t]]))
        path = mat"RunWithoutImages($mX, $σ)"
        if size(path, 2) == 1
            acc[t] = acc2fβ(calc_path_accuracy(reshape(path, :, 1), M[1:t-1], length.(X[1:t]))...)
        elseif size(path, 2) == 0
            acc[t] = acc2fβ(calc_path_accuracy(-ones(size(path, 1), 1), M[1:t-1], length.(X[1:t]))...)
        else
            acc[t] = acc2fβ(calc_path_accuracy(path, M[1:t-1], length.(X[1:t]))...)
        end
    end
    return acc
end

max_stars_repo_name: szcf-weiya/BaxterAlgorithms

max_issues_repo_name: szcf-weiya/BaxterAlgorithms

max_forks_repo_name: szcf-weiya/BaxterAlgorithms


text:
# This file was formerly a part of Julia. License is MIT: https://julialang.org/license

include("../perfutil.jl")

function cat2d_perf(n, iter)
    a = rand(n,n)
    b = rand(n,n)
    for i=1:iter
        c = [a b; b a]
    end
end

function cat2d_perf2(n, iter)
    a = rand(n,n)
    b = rand(n,n)
    for i=1:iter
        c = Matrix{Float64}(uninitialized, 2n,2n)
        c[  1:n,  1:n] = a
        c[  1:n,  n+1:end] = b
        c[n+1:end,1:n] = b
        c[n+1:end,n+1:end] = a
    end
end

function hcat_perf(n, iter)
    a = rand(n,n)
    b = rand(n,n)
    for i=1:iter
        c = [a b b a]
    end
end

function hcat_perf2(n, iter)
    a = rand(n,n)
    b = rand(n,n)
    for i=1:iter
        c = Matrix{Float64}(uninitialized, n, 4n)
        c[:,    1:  n] = a
        c[:,  n+1: 2n] = b
        c[:, 2n+1: 3n] = b
        c[:, 3n+1:end] = a
    end
end

function vcat_perf(n, iter)
    a = rand(n,n)
    b = rand(n,n)
    for i = 1:iter
        c = [a; b; b; a]
    end
end

function vcat_perf2(n, iter)
    a = rand(n,n)
    b = rand(n,n)
    for i=1:iter
        c = Matrix{Float64}(uninitialized, 4n, n)
        c[   1: n, :] = a
        c[ n+1:2n, :] = b
        c[2n+1:3n, :] = b
        c[3n+1:4n, :] = a
    end
end

function catnd_perf(n, iter)
    a = rand(1,n,n,1)
    b = rand(1,n,n)
    for i = 1:iter
        c = cat(3, a, b, b, a)
    end
end

function catnd_perf2(n, iter)
    a = rand(1,n,n,1)
    b = rand(1,n,n)
    for i = 1:iter
        c = Array{Float64}(uninitialized, 1, n, 4n, 1)
        c[1,:,   1: n,1] = a
        c[1,:, n+1:2n,1] = b
        c[1,:,2n+1:3n,1] = b
        c[1,:,3n+1:4n,1] = a
    end
end

problemsizes = [(5, 20000, "small"), (500, 2, "large")]
testdata = [(cat2d_perf,  "hvcat",        "horizontal/vertical matrix concatenation", problemsizes),
            (cat2d_perf2, "hvcat_setind", "horizontal/vertical matrix concatenation using setindex", problemsizes),
            (hcat_perf,   "hcat",         "horizontal matrix concatenation", problemsizes),
            (hcat_perf2,  "hcat_setind",  "horizontal matrix concatenation using setindex", problemsizes),
            (vcat_perf,   "vcat",         "vertical matrix concatenation", problemsizes),
            (vcat_perf2,  "vcat_setind",  "vertical matrix concatenation using setindex", problemsizes),
            (catnd_perf,  "catnd",        "N-dimensional matrix concatenation", problemsizes),
            (catnd_perf2, "catnd_setind", "N-dimensional matrix concatenation using setindex", problemsizes)]
include("../perfgeneric.jl")

max_stars_repo_name: JuliaTagBot/Benchmarks

max_issues_repo_name: JuliaArchive/Benchmarks

max_forks_repo_name: JuliaArchive/Benchmarks


text:
using DistributionVectors
using Test, Random

@testset "DistributionVector" begin
    include("distributionvector.jl")
end

@testset "autocor" begin
    include("autocor.jl")
end

@testset "sumdistributionvector" begin
    include("sumdistributionvector.jl")
end

const tests = [
    "normal",
    "lognormal",
]

printstyled("Running Distribution tests:\n", color=:blue)

for t in tests
    @testset "Test $t" begin
        Random.seed!(345679)
        include("$t.jl")
    end
end

# print method ambiguities
# println("Potentially stale exports: ")
# display(Test.detect_ambiguities(DistributionVectors))
# println()


text:
struct DBFGS{T1,T2,T3} <: QuasiNewton{T1}
    approx::T1
    theta::T2
    P::T3
end
DBFGS(approx) = DBFGS(approx, 0.2, nothing)
DBFGS(; inverse = true, theta = 0.2) = DBFGS(inverse ? Inverse() : Direct(), theta, nothing)
hasprecon(::DBFGS{<:Any,<:Any,<:Nothing}) = NoPrecon()
hasprecon(::DBFGS{<:Any,<:Any,<:Any}) = HasPrecon()
summary(dbfgs::DBFGS{Inverse}) = "Inverse Damped BFGS"
summary(dbfgs::DBFGS{Direct}) = "Direct Damped BFGS"

function update!(scheme::DBFGS{<:Direct,<:Any,<:Any}, B, s, y)
    # We could write this as
    #     B .+= (y*y')/dot(s, y) - (B*s)*(s'*B)/(s'*B*s)
    #     B .+= (y*y')/dot(s, y) - b*b'/dot(s, b)
    # where b = B*s
    # But instead, we split up the calculations. First calculate the denominator
    # in the first term
    σ = dot(s, y)
    ρ = inv(σ) # scalar
    # Then calculate the vector b
    b = B * s # vector temporary
    sb = dot(s, b)
    if σ ≥ scheme.theta * sb
        θ = 1.0
        # Calculate one vector divided by dot(s, b)
        ρbb = inv(sb) * b
        # And calculate
        B .+= (inv(σ) * y) * y' .- ρbb * b'
    else
        θ = 0.8 * sb / (sb - σ)
        r = y * θ + (1 - θ) * b
        # Calculate one vector divided by dot(s, b)
        ρbb = inv(dot(s, b)) * b
        # And calculate
        B .+= (inv(dot(s, r)) * r) * r' .- ρbb * b'
    end
end
function update(scheme::DBFGS{<:Direct,<:Any}, B, s, y)
    # As above, but out of place

    σ = dot(s, y)
    b = B * s
    sb = dot(s, b)
    if σ ≥ scheme.theta * sb
        θ = 1.0
        # Calculate one vector divided by dot(s, b)
        ρbb = inv(sb) * b
        # And calculate
        return B = B .+ (inv(σ) * y) * y' .- ρbb * b'
    else
        θ = 0.8 * sb / (sb - σ)
        r = y * θ + (1 - θ) * b
        # Calculate one vector divided by dot(s, b)
        ρbb = inv(dot(s, b)) * b
        # And calculate
        return B = B .+ (inv(dot(s, r)) * r) * r' .- ρbb * b'
    end
end
function update(scheme::DBFGS{<:Inverse,<:Any}, H, s, y)
    σ = dot(s, y)
    ρ = inv(σ)
    #   if isfinite(ρ)
    C = (I - ρ * s * y')
    H = C * H * C' + ρ * s * s'
    #   end
    H
end
function update!(scheme::DBFGS{<:Inverse,<:Any}, H, s, y)
    σ = dot(s, y)
    ρ = inv(σ)

    if isfinite(ρ)
        Hy = H * y
        H .= H .+ ((σ + y' * Hy) .* ρ^2) * (s * s')
        Hys = Hy * s'
        Hys .= Hys .+ Hys'
        H .= H .- Hys .* ρ
    end
    H
end

function update!(scheme::DBFGS{<:Inverse,<:Any}, A::UniformScaling, s, y)
    update(scheme, A, s, y)
end
function update!(scheme::DBFGS{<:Direct,<:Any}, A::UniformScaling, s, y)
    update(scheme, A, s, y)
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
using SafeTestsets

@safetestset "Harmonic Oscillator                                                             " begin include("harmonic_oscillator_tests.jl") end
@safetestset "Kubo Oscillator                                                                 " begin include("kubo_oscillator_tests.jl") end
@safetestset "Lorenz Attractor                                                                " begin include("lorenz_attractor_tests.jl") end
@safetestset "Lotka-Volterra 2D                                                               " begin include("lotka_volterra_2d_tests.jl") end
@safetestset "Lotka-Volterra 2D with singular Lagrangian                                      " begin include("lotka_volterra_2d_singular_tests.jl") end
@safetestset "Lotka-Volterra 2D with symmetric Lagrangian                                     " begin include("lotka_volterra_2d_symmetric_tests.jl") end
@safetestset "Lotka-Volterra 2D with symmetric Lagrangian with gauge terms                    " begin include("lotka_volterra_2d_gauge_tests.jl") end
@safetestset "Lotka-Volterra 3D                                                               " begin include("lotka_volterra_3d_tests.jl") end
@safetestset "Lotka-Volterra 4D                                                               " begin include("lotka_volterra_4d_tests.jl") end
@safetestset "Lotka-Volterra 4D (Lagrangian)                                                  " begin include("lotka_volterra_4d_lagrangian_tests.jl") end
@safetestset "Point Vortices                                                                  " begin include("point_vortices_tests.jl") end


text:
using Test, ReferenceTests, BSON
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil

ENV["USE_GPU"] = false
ENV["DO_VIZ"]  = false
ENV["DO_SAVE"] = false
ENV["DO_SAVE_VIZ"] = false
ENV["NX"] = 64
ENV["USE_RETURN"] = true

# # Unit tests
# @testset "Unit-tests" begin
#     # some unit tests
# end

# Reference test using ReferenceTests.jl
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2) for (v1,v2) in zip(values(d1), values(d2))])

## 1D tests
@reset_parallel_stencil()
include("../scripts/diff_1D_lin.jl")
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d1d1 = Dict(:X=> xc[inds], :H=>H[inds])

@reset_parallel_stencil()
include("../scripts/diff_1D_linstep.jl")
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d1d2 = Dict(:X=> xc[inds], :H=>H[inds])

@reset_parallel_stencil()
include("../scripts/diff_1D_nonlin.jl")
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d1d3 = Dict(:X=> xc[inds], :H=>H[inds])

@reset_parallel_stencil()
include("../scripts/diff_1D_react.jl")
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d1d4 = Dict(:X=> xc[inds], :H=>H[inds])

@testset "Reference-tests diffusion 1D" begin
    @test_reference "reftest-files/test_diff_1D_lin.bson" d1d1 by=comp
    @test_reference "reftest-files/test_diff_1D_linstep.bson" d1d2 by=comp
    @test_reference "reftest-files/test_diff_1D_nonlin.bson" d1d3 by=comp
    @test_reference "reftest-files/test_diff_1D_react.bson" d1d4 by=comp
end

@reset_parallel_stencil()


text:
"""
`module Experimental`

TODO: write documentation
"""
module Experimental

using JuLIP: AbstractAtoms, dofs, set_dofs!, atomdofs
using JuLIP.Potentials
using ForwardDiff

export newton!, constrained_bond_newton!

function newton(x, g, h; maxsteps=10, tol=1e-5)
    for i in 1:maxsteps
        gi = g(x)
        norm(gi, Inf) < tol && break
        x = x - h(x) \ gi
        @show i, norm(gi, Inf)
    end
    return x
end

function newton!(a::AbstractAtoms; maxsteps=20, tol=1e-10)
    x = newton(dofs(a),
               x -> gradient(a, x),
               x -> hessian(a, x))
    set_dofs!(a, x)
    return x
end

function constrained_bond_newton!(a::AbstractAtoms, i::Integer, j::Integer, bondlength::Float64;
                                  maxsteps=20, tol=1e-10)
    I1 = atomdofs(a, i)
    I2 = atomdofs(a, j)
    I1I2 = [I1; I2]

    # bondlength of target bond
    blen(x) = norm(x[I2] - x[I1])

    # constraint function
    c(x) = blen(x) - bondlength

    # gradient of constraint
    function dc(x)
        r = zeros(length(x))
        r[I1] = (x[I1]-x[I2])/blen(x)
        r[I2] = (x[I2]-x[I1])/blen(x)
        return r
    end

    # define some auxiliary index arrays that start at 1 to
    # allow the ForwardDiff hessian to be done only on relevent dofs, x[I1I2]
    _I1 = 1:length(I1)
    _I2 = length(I1)+1:length(I1)+length(I2)
    _blen(x) = norm(x[_I2] - x[_I1])
    _c(x) = _blen(x) - bondlength

    function ddc(x)
        s = spzeros(length(x), length(x))
        s[I1I2, I1I2] = ForwardDiff.hessian(_c, x[I1I2])
        return s
    end

    # Define Lagrangian L(x, λ) = E - λC and its gradient and hessian
    L(x, λ) = energy(a, x) - λ*c(x)
    L(z) = L(z[1:end-1], z[end])

    dL(x, λ) = [gradient(a, x) - λ * dc(x); -c(x)]
    dL(z) = dL(z[1:end-1], z[end])

    ddL(x, λ) = [(hessian(a, x) - λ * ddc(x)) -dc(x); -dc(x)' 0.]
    ddL(z) = ddL(z[1:end-1], z[end])

    # Use Newton scheme to find saddles L where ∇L = 0
    z = [dofs(a); -10.0*c(dofs(a))]
    z = newton(z, dL, ddL, maxsteps=maxsteps, tol=tol)
    set_dofs!(a, z[1:end-1])
    return z
end

end


text:
module KernelsTests

using CellwiseValues
using TensorValues

k = NumberKernelFromFunction(-)

test_number_kernel(k,1,4,3)

test_number_kernel(k,-1.0,1.0)

k = NumberKernelFromFunction(sum)

test_number_kernel(k,7,[1,2,4])

k = ArrayKernelFromBroadcastedFunction(*)

test_array_kernel(k,[3,4,3],[1,2,3],[3,2,1])

test_array_kernel(k,[6,4,2],2,[3,2,1])

test_array_kernel(k,[6,4,2],[3,2,1],2)

test_array_kernel(k,[6,4,2],[3,2,1],2,1)

u = rand(Int,2,3,1)
v = rand(Int,1,3,4)
w = u .* v
test_array_kernel(k,w,u,v)
test_array_kernel(k,w,u,v,1)

k = ArrayKernelFromBroadcastedFunction(+)

v1 = VectorValue(2,3)
v2 = VectorValue(3,2)
v3 = VectorValue(1,2)
u = [v1,v2,v3]
v = v1
w = u .+ v
test_array_kernel(k,w,u,v)

w = broadcast(+,u,v,0)
test_array_kernel(k,w,u,v,0)

end # module