text:
"""
`module  JAC.DoubleAutoIonization`
    ... a submodel of JAC that contains all methods for computing double Auger and autoionization amplitudes and rates.
"""
module DoubleAutoIonization

    using Printf, JAC, ..AngularMomentum, ..AtomicState, ..AutoIonization, ..Basics, ..ManyElectron, ..TableStrings

    """
    `struct  DoubleAutoIonization.Settings`  ... defines a type for the settings in estimating double-Auger and autoionization rates.

        + green                   ::Array{GreenChannel,1}   ... Precalculated and user-specified Green function of the ion.
        + NoEnergySharings        ::Int64                   ... Number of energy sharings that are used in the computations for each line.
        + printBefore             ::Bool                    ... True, if all energies and lines are printed before their evaluation.
        + 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(), ...
        + lineSelection           ::LineSelection           ... Specifies the selected levels, if any.
    """
    struct Settings
        green                     ::Array{GreenChannel,1}
        NoEnergySharings          ::Int64
        printBefore               ::Bool
        minAugerEnergy            ::Float64
        maxAugerEnergy            ::Float64
        maxKappa                  ::Int64
        operator                  ::AbstractEeInteraction
        lineSelection             ::LineSelection
    end

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

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

    """
    `struct  DoubleAutoIonization.ReducedChannel`
        ... defines a type for a DoubleAutoIonization (reduced) channel to help characterize a scattering (continuum) state of many
            electron-states with two free electrons.

        + energy1        ::Float64              ... energy of the partial wave_1
        + kappa1         ::Int64                ... partial-wave of the free electron_1
        + phase1         ::Float64              ... phase of the partial wave_1
        + xSymmetry      ::LevelSymmetry        ... angular momentum and parity after coupling electron1.
        + energy2        ::Float64              ... energy of the partial wave_2
        + kappa2         ::Int64                ... partial-wave of the free electron_2
        + phase2         ::Float64              ... phase of the partial wave_2
        + amplitude      ::Complex{Float64}     ... Auger amplitude associated with the given channel.
    """
    struct  ReducedChannel
        energy1          ::Float64
        kappa1           ::Int64
        phase1           ::Float64
        xSymmetry        ::LevelSymmetry
        energy2          ::Float64
        kappa2           ::Int64
        phase2           ::Float64
        amplitude        ::Complex{Float64}
    end

    # `Base.show(io::IO, channel::DoubleAutoIonization.ReducedChannel)`
    #       ... prepares a proper printout of the variable cannel::DoubleAutoIonization.ReducedChannel.
    function Base.show(io::IO, channel::DoubleAutoIonization.ReducedChannel)
        ## sa = "reduced DA channel for J^P=(channel.tSymmetry)  with " *
        ##      "kappa1=$(channel.kappa1), energy1=$(channel.energy1), phase1=$(channel.phase1), " *
        ##      "kappa2=$(channel.kappa2), energy2=$(channel.energy2), phase2=$(channel.phase2), amp=$(channel.amplitude) "
        ## println(io, sa)
        println(io, "multipole:                    $(channel.multipole)  ")
        println(io, "gauge:                        $(channel.gauge)  ")
        println(io, "omega:                        $(channel.omega)  ")
        println(io, "energy1:                      $(channel.energy1)  ")
        println(io, "kappa1:                       $(channel.kappa1)  ")
        println(io, "phase1:                       $(channel.phase1)  ")
        println(io, "xSymmetry:                    $(channel.xSymmetry)  ")
        println(io, "energy2:                      $(channel.energy2)  ")
        println(io, "kappa2:                       $(channel.kappa2)  ")
        println(io, "phase2:                       $(channel.phase2)  ")
        println(io, "amplitude:                    $(channel.amplitude)  ")
    end

    """
    `struct  DoubleAutoIonization.Sharing`
        ... defines a type for a DoubleAutoIonization sharing to help characterize energy sharing between the two emitted electrons.

        + epsilon1         ::Float64         ... Energy of (free) electron 1.
        + epsilon2         ::Float64         ... Energy of (free) electron 2.
        + weight           ::Float64         ... Gauss-Lengendre weight of this sharing for energy-integrated quantities.
        + differentialRate ::EmProperty      ... differential cross section of this energy sharing.
        + channels         ::Array{DoubleAutoIonization.ReducedChannel,1}  ... List of DoubleAutoIonization (reduced) channels of this line.
    """
    struct  Sharing
        epsilon1           ::Float64
        epsilon2           ::Float64
        weight             ::Float64
        differentialRate   ::EmProperty
        channels           ::Array{DoubleAutoIonization.ReducedChannel,1}
    end

    # `Base.show(io::IO, sharing::DoubleAutoIonization.Sharing)`  ... prepares a proper printout of the variable sharing::DoubleAutoIonization.Sharing.
    function Base.show(io::IO, sharing::DoubleAutoIonization.Sharing)
        println(io, "epsilon1:               $(sharing.epsilon1)  ")
        println(io, "epsilon2:               $(sharing.epsilon2)  ")
        println(io, "differentialCs:         $(sharing.differentialCs)  ")
        println(io, "channels:               $(sharing.channels)  ")
    end

    """
    `struct  DoubleAutoIonization.Line`  ... defines a type for a double Auger line that includes sharings and their reduced amplitudes.

        + initialLevel   ::Level            ... initial-(state) level
        + finalLevel     ::Level            ... final-(state) level
        + totalRate      ::EmProperty       ... Total rate of this line.
        + sharings       ::Array{DoubleAutoIonization.Sharing,1}  ... List of DoubleAutoIonization sharings of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        totalRate        ::EmProperty
        sharings         ::Array{DoubleAutoIonization.Sharing,1}
    end

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

    """
    `DoubleAutoIonization.computeAmplitudesProperties(line::DoubleAutoIonization.Line, grid::Radial.Grid, settings::DoubleAutoIonization.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::DoubleAutoIonization.Line is returned for which the amplitudes
            and properties have now been evaluated.
    """
    function  computeAmplitudesProperties(line::DoubleAutoIonization.Line, grid::Radial.Grid, settings::DoubleAutoIonization.Settings)
        newSharings = DoubleAutoIonization.Sharing[]
        for sharing in line.sharings
            newChannels = DoubleAutoIonization.ReducedChannel[]
            for ch in sharing.channels
                # Generate a continuum orbital
                phase     = 1.3
                amplitude = 1.0im
                push!( newChannels, ReducedChannel(ch.energy1, ch.kappa1, phase, ch.xSymmetry, ch.energy2, ch.kappa2, phase, amplitude) )
            end
            # Calculate the differential rate
            diffRate = EmProperty(-1., -1.)
            push!( newSharings, DoubleAutoIonization.Sharing(sharing.epsilon1, sharing.epsilon2, sharing.weight, diffRate, newChannels) )
        end
        # Calculate the total rate
        totalRate = EmProperty(-1., -1.)
        line = DoubleAutoIonization.Line(line.initialLevel, line.finalLevel, totalRate, newSharings)
        return( line )
    end

    """
    `DoubleAutoIonization.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                       settings::DoubleAutoIonization.Settings; output=true)`
        ... to compute the double-Auger transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{DoubleAutoIonization.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                           settings::DoubleAutoIonization.Settings; output=true)
        println("")
        printstyled("DoubleAutoIonization.computeLines(): The computation of double-Auger rates starts now ... \n", color=:light_green)
        printstyled("----------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        lines = DoubleAutoIonization.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    DoubleAutoIonization.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = DoubleAutoIonization.Line[]
        for  line in lines
            newLine = DoubleAutoIonization.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        DoubleAutoIonization.displayTotalRates(stdout, lines, settings)
        DoubleAutoIonization.displayDifferentialRates(stdout, lines, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   DoubleAutoIonization.displayTotalRates(iostream, lines, settings)
                           DoubleAutoIonization.displayDifferentialRates(iostream, lines, settings)     end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `DoubleAutoIonization.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::DoubleAutoIonization.Settings)`
        ... to determine a list of DoubleAutoIonization.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{DoubleAutoIonization.Line,1} is
            returned. Apart from the level specification and sharing, all physical properties are set to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::DoubleAutoIonization.Settings)
        lines = DoubleAutoIonization.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 < settings.minAugerEnergy  ||  energy > settings.maxAugerEnergy    continue   end
                    sharings = DoubleAutoIonization.determineSharingsAndChannels(fLevel, iLevel, energy, settings)
                    push!( lines, DoubleAutoIonization.Line(iLevel, fLevel, EmProperty(0., 0.,), sharings) )
                end
            end
        end
        return( lines )
    end

    """
    `DoubleAutoIonization.determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, energy::Float64, settings::DoubleAutoIonization.Settings)`
        ... to determine a list of DoubleAutoIonization Sharing's and Channel's for a transitions from the initial to final level and by taking into
            account the particular settings of for this computation; an Array{DoubleAutoIonization.Sharing,1} is returned.
    """
    function determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, energy::Float64, settings::DoubleAutoIonization.Settings)
        sharings = DoubleAutoIonization.Sharing[];    eSharings = Basics.determineEnergySharings(energy, settings.NoEnergySharings)
        for  es in eSharings
            epsilon1  = es[1];    epsilon2 = es[2];    weight = es[3]
            channels  = DoubleAutoIonization.ReducedChannel[];
            symi      = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
            couplings = AngularMomentum.allowedDoubleKappaCouplingSequence(symf, symi, settings.maxKappa)
            for  coupling  in  couplings
                push!(channels, ReducedChannel(epsilon1, coupling[1], 0., coupling[2], epsilon2, coupling[3], 0., Complex(0.)) )
            end
            push!(sharings, DoubleAutoIonization.Sharing(epsilon1, epsilon2, weight, EmProperty(0., 0.), channels) )
        end
        return( sharings )
    end

    """
    `DoubleAutoIonization.displayDifferentialRates(stream::IO, lines::Array{DoubleAutoIonization.Line,1}, settings::DoubleAutoIonization.Settings)`
        ... to display all differential rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayDifferentialRates(stream::IO, lines::Array{DoubleAutoIonization.Line,1}, settings::DoubleAutoIonization.Settings)
        #
        # First, print lines and sharings
        nx = 128
        println(stream, " ")
        println(stream, "  Energy-differential rates of selected double-Auger lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(38, "Energies (all in " * TableStrings.inUnits("energy") * ")"; na=3);
        sb = sb * TableStrings.flushleft(38, "  i -- f        epsilon_1     epsilon_2"; na=3)
        sa = sa * TableStrings.center(14, "Weight"; na=0);                            sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(34, "Cou --  diff. rate  -- Bab"; na=3)
        sb = sb * TableStrings.center(34, TableStrings.inUnits("rate") * "        " *
                                          TableStrings.inUnits("rate"); na=3)
        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)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            for  (is, sharing)  in  enumerate(line.sharings)
                if  is == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon1))                 * "    "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon2))                 * "    "
                sb = sb * @sprintf("%.4e",                                              sharing.weight)                    * "       "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", sharing.differentialRate.Coulomb))   * "   "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", sharing.differentialRate.Babushkin)) * "   "
                println(stream,  sb )
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `DoubleAutoIonization.displayLines(lines::Array{DoubleAutoIonization.Line,1})`
        ... to display a list of lines, sharings 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{DoubleAutoIonization.Line,1})
        #
        # First, print lines and sharings
        nx = 82
        println(" ")
        println("  Selected double-Auger lines & sharings:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(54, "Energies (all in )" * TableStrings.inUnits("energy") * ")"; na=5);
        sb = sb * TableStrings.flushleft(54, "   i -- f       epsilon_1     epsilon_2"; na=5)
        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)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            for  (is, sharing)  in  enumerate(line.sharings)
                if  is == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon1)) * "    "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon2)) * "    "
                println( sb )
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        #
        # Second, print lines and channles
        nx = 123
        println(" ")
        println("  Selected photo-double ionization lines & channels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(77, "Channels (all energies in " * TableStrings.inUnits("energy") * ")" ; na=5);
        sb = sb * TableStrings.flushleft(77, "(epsilon, kappa)_1   -->    J^P_x  -->   (epsilon, kappa)_2   -->  J^P_t = J^P_i"; na=5)
        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)
            for (ic, ch) in enumerate(line.sharings[1].channels)
                if  ic == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                skappa1 = "  " * string(ch.kappa1);    skappa2 = "  " * string(ch.kappa2);    sxsym = "      " * string(ch.xSymmetry)
                sb = sb * " (" * @sprintf("%.3e", Defaults.convertUnits("energy: from atomic", ch.energy1)) * ", " * skappa1[end-2:end] * ")    -->"
                sb = sb * sxsym[end-8:end] * "  -->   "
                sb = sb * " (" * @sprintf("%.3e", Defaults.convertUnits("energy: from atomic", ch.energy2)) * ", " * skappa2[end-2:end] * ")    -->     "
                sb = sb * string(isym)
                println( sb )
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `DoubleAutoIonization.displayTotalRates(stream::IO, lines::Array{DoubleAutoIonization.Line,1}, settings::DoubleAutoIonization.Settings)`
        ... to display all total rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayTotalRates(stream::IO, lines::Array{DoubleAutoIonization.Line,1}, settings::DoubleAutoIonization.Settings)
        #
        # First, print lines and sharings
        nx = 88
        println(stream, " ")
        println(stream, "  Total rates of selected double-Auger lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=2);                         sb = sb * TableStrings.hBlank(21)
        sa = sa * TableStrings.center(12, "i--Energy--f"; na=3)
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(30, "Cou --   total rate   -- Bab"; na=4)
        sb = sb * TableStrings.center(30, TableStrings.inUnits("rate") * "          " *
                                          TableStrings.inUnits("rate"); na=3)
        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)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))                 * "      "
            #
            sb = sa * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", line.totalRate.Coulomb))   * "     "
            sb = sb * @sprintf("%.5e", Defaults.convertUnits("rate: from atomic", line.totalRate.Babushkin)) * "   "
            println(stream,  sb )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
using Documenter, ESDL

makedocs(
    modules  = [ESDL],
    clean    = true,
    doctest  = false,
    format   = Documenter.HTML(),
    sitename = "ESDL.jl",
    authors  = "Fabian Gans",
    pages    = [
        "Home" => "index.md",
        "Examples" => Any[
            "examples/time_mean.md",
        ],
        "API Documentation" => Any[
            "Index" => "API/api_index.md",
            "Accessing the Data Cube" => "API/cube_access.md",
            "Analisys" => "API/analysis.md",
            "Loading and Saving Results" => "API/iotools.md"
        ],
        "Other functions" => "lib/misc.md"
        ]
)

deploydocs(
    repo   = "github.com/esa-esdl/ESDL.jl.git",
    push_preview = true
)

ext: jl


text:
module ROCKernels

import AMDGPU
import AMDGPU: rocfunction, HSAAgent, HSAQueue, HSASignal, HSAStatusSignal, Mem
import SpecialFunctions
import StaticArrays
import StaticArrays: MArray
import Cassette
import Adapt
import KernelAbstractions

export ROCDevice

const FREE_QUEUES = HSAQueue[]
const QUEUES = HSAQueue[]
const QUEUE_GC_THRESHOLD = Ref{Int}(16)

# This code is loaded after an `@init` step
if haskey(ENV, "KERNELABSTRACTIONS_QUEUES_GC_THRESHOLD")
    global QUEUE_GC_THRESHOLD[] = parse(Int, ENV["KERNELABSTRACTIONS_QUEUES_GC_THRESHOLD"])
end

## Stream GC
# Simplistic stream gc design in which when we have a total number
# of streams bigger than a threshold, we start scanning the streams
# and add them back to the freelist if all work on them has completed.
# Alternative designs:
# - Enqueue a host function on the stream that adds the stream back to the freelist
# - Attach a finalizer to events that adds the stream back to the freelist
# Possible improvements
# - Add a background task that occasionally scans all streams
# - Add a hysterisis by checking a "since last scanned" timestamp
# - Add locking
const QUEUE_GC_LOCK = Threads.ReentrantLock()
function next_queue()
    lock(QUEUE_GC_LOCK) do
        if !isempty(FREE_QUEUES)
            return pop!(FREE_QUEUES)
        end

        if length(QUEUES) > QUEUE_GC_THRESHOLD[]
            for queue in QUEUES
                #if AMDGPU.queued_packets(queue) == 0
                    push!(FREE_QUEUES, queue)
                #end
            end
        end

        if !isempty(FREE_QUEUES)
            return pop!(FREE_QUEUES)
        end

        # FIXME: Which agent?
        queue = HSAQueue()
        push!(QUEUES, queue)
        return queue
    end
end

import KernelAbstractions: Event, CPUEvent, NoneEvent, MultiEvent, CPU, GPU, isdone, failed

struct ROCDevice <: GPU end

struct ROCEvent{T<:Union{AMDGPU.HSA.Signal,HSAStatusSignal}} <: Event
    event::T
end
ROCEvent(signal::HSASignal) = ROCEvent(signal.signal[])

cpuwait(ev::ROCEvent{AMDGPU.HSA.Signal}) = wait(HSASignal(ev.event))
cpuwait(ev::ROCEvent{HSAStatusSignal}) = wait(ev.event)
gpuwait(ev::ROCEvent{AMDGPU.HSA.Signal}) = AMDGPU.device_signal_wait(ev.event, 0)
gpuwait(ev::ROCEvent{HSAStatusSignal}) = wait(ev.event)

Adapt.adapt_storage(::AMDGPU.Adaptor, ev::ROCEvent{HSAStatusSignal}) =
    ROCEvent(ev.event.signal.signal[])

failed(::ROCEvent) = false # FIXME
isdone(ev::ROCEvent) = true # FIXME

function Event(::ROCDevice)
    queue = AMDGPU.get_default_queue()
    event = AMDGPU.barrier_and!(queue, AMDGPU.active_kernels(queue))
    ROCEvent(event.signal)
end

import Base: wait

wait(ev::ROCEvent, progress=nothing) = wait(CPU(), ev, progress)

function wait(::CPU, ev::ROCEvent, progress=nothing)
    if progress === nothing
        cpuwait(ev)
    else
        # FIXME: One-shot wait
        while !cpuwait(ev)
            progress()
            # do we need to `yield` here?
        end
    end
end

wait(::ROCDevice, ev::ROCEvent, progress=nothing, queue=nothing) = gpuwait(ev)
wait(::ROCDevice, ev::NoneEvent, progress=nothing, queue=nothing) = nothing

function wait(::ROCDevice, ev::MultiEvent, progress=nothing, queue=AMDGPU.get_default_queue())
    dependencies = collect(ev.events)
    rocdeps = map(d->d.event isa HSAStatusSignal ? d.event.signal.signal[] : d.event, filter(d->d isa ROCEvent, dependencies))
    wait.(rocdeps)
    #isempty(rocdeps) || wait(AMDGPU.barrier_and!(queue, rocdeps))
    for event in filter(d->!(d isa ROCEvent), dependencies)
        wait(ROCDevice(), event, progress)
    end
end

# TODO: Create a hostcall that waits on a specified Task
wait(::ROCDevice, ev::CPUEvent, progress=nothing, queue=nothing) = error("GPU->CPU wait not implemented")

function KernelAbstractions.async_copy!(::ROCDevice, A, B; dependencies=nothing, progress=nothing)
    # TODO: Register A and B
    wait(ROCDevice(), MultiEvent(dependencies), progress, AMDGPU.get_default_queue())
    GC.@preserve A B begin
        destptr = pointer(A)
        srcptr  = pointer(B)
        N       = length(A)
        # TODO: memory_copy_multiple, or custom copy kernel
        AMDGPU.HSA.memory_copy(destptr, srcptr, N) |> AMDGPU.check
    end

    return ROCEvent(HSASignal(0))
end

import KernelAbstractions: Kernel, StaticSize, DynamicSize, partition, blocks, workitems, launch_config

###
# Kernel launch
###
function launch_config(kernel::Kernel{ROCDevice}, ndrange, workgroupsize)
    if ndrange isa Integer
        ndrange = (ndrange,)
    end
    if workgroupsize isa Integer
        workgroupsize = (workgroupsize, )
    end

    # partition checked that the ndrange's agreed
    if KernelAbstractions.ndrange(kernel) <: StaticSize
        ndrange = nothing
    end

    iterspace, dynamic = if KernelAbstractions.workgroupsize(kernel) <: DynamicSize &&
        workgroupsize === nothing
        # use ndrange as preliminary workgroupsize for autotuning
        partition(kernel, ndrange, ndrange)
    else
        partition(kernel, ndrange, workgroupsize)
    end

    return ndrange, workgroupsize, iterspace, dynamic
end

function threads_to_workgroupsize(threads, ndrange)
    total = 1
    return map(ndrange) do n
        x = min(div(threads, total), n)
        total *= x
        return x
    end
end

function (obj::Kernel{ROCDevice})(args...; ndrange=nothing, dependencies=nothing, workgroupsize=nothing, progress=nothing)

    ndrange, workgroupsize, iterspace, dynamic = launch_config(obj, ndrange, workgroupsize)
    # this might not be the final context, since we may tune the workgroupsize
    ctx = mkcontext(obj, ndrange, iterspace)
    #= TODO: Autotuning
    kernel = AMDGPU.@roc launch=false name=String(nameof(obj.f)) Cassette.overdub(ctx, obj.f, args...)

    # figure out the optimal workgroupsize automatically
    if KernelAbstractions.workgroupsize(obj) <: DynamicSize && workgroupsize === nothing
        config = AMDGPU.launch_configuration(kernel.fun; max_threads=prod(ndrange))
        workgroupsize = threads_to_workgroupsize(config.threads, ndrange)
        iterspace, dynamic = partition(obj, ndrange, workgroupsize)
        ctx = mkcontext(obj, ndrange, iterspace)
    end

    # If the kernel is statically sized we can tell the compiler about that
    if KernelAbstractions.workgroupsize(obj) <: StaticSize
        maxthreads = prod(KernelAbstractions.get(KernelAbstractions.workgroupsize(obj)))
    else
        maxthreads = nothing
    end
    =#

    nblocks = length(blocks(iterspace))
    threads = length(workitems(iterspace))

    if nblocks == 0
        return MultiEvent(dependencies)
    end

    queue = next_queue()
    wait(ROCDevice(), MultiEvent(dependencies), progress, queue)

    # Launch kernel
    event = AMDGPU.@roc(groupsize=threads, gridsize=nblocks*threads, queue=queue,
                        name=String(nameof(obj.f)), # TODO: maxthreads=maxthreads,
                        Cassette.overdub(ctx, obj.f, args...))

    return ROCEvent(event.event)
end

Cassette.@context ROCCtx

import KernelAbstractions: CompilerMetadata, CompilerPass, DynamicCheck, LinearIndices
import KernelAbstractions: __index_Local_Linear, __index_Group_Linear, __index_Global_Linear, __index_Local_Cartesian, __index_Group_Cartesian, __index_Global_Cartesian, __validindex, __print
import KernelAbstractions: mkcontext, expand, __iterspace, __ndrange, __dynamic_checkbounds

function mkcontext(kernel::Kernel{ROCDevice}, _ndrange, iterspace)
    metadata = CompilerMetadata{KernelAbstractions.ndrange(kernel), DynamicCheck}(_ndrange, iterspace)
    Cassette.disablehooks(ROCCtx(pass = CompilerPass, metadata=metadata))
end
function mkcontext(kernel::Kernel{ROCDevice}, I, _ndrange, iterspace, ::Dynamic) where Dynamic
    metadata = CompilerMetadata{KernelAbstractions.ndrange(kernel), Dynamic}(I, _ndrange, iterspace)
    Cassette.disablehooks(ROCCtx(pass = CompilerPass, metadata=metadata))
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__index_Local_Linear))
    return AMDGPU.threadIdx().x
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__index_Group_Linear))
    return AMDGPU.blockIdx().x
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__index_Global_Linear))
    I =  @inbounds expand(__iterspace(ctx.metadata), AMDGPU.blockIdx().x, AMDGPU.threadIdx().x)
    # TODO: This is unfortunate, can we get the linear index cheaper
    @inbounds LinearIndices(__ndrange(ctx.metadata))[I]
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__index_Local_Cartesian))
    @inbounds workitems(__iterspace(ctx.metadata))[AMDGPU.threadIdx().x]
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__index_Group_Cartesian))
    @inbounds blocks(__iterspace(ctx.metadata))[AMDGPU.blockIdx().x]
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__index_Global_Cartesian))
    return @inbounds expand(__iterspace(ctx.metadata), AMDGPU.blockIdx().x, AMDGPU.threadIdx().x)
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__validindex))
    if __dynamic_checkbounds(ctx.metadata)
        I = @inbounds expand(__iterspace(ctx.metadata), AMDGPU.blockIdx().x, AMDGPU.threadIdx().x)
        return I in __ndrange(ctx.metadata)
    else
        return true
    end
end

import KernelAbstractions: groupsize, __groupsize, __workitems_iterspace, add_float_contract, sub_float_contract, mul_float_contract

KernelAbstractions.generate_overdubs(@__MODULE__, ROCCtx)

###
# AMDGPU specific method rewrites
###

@inline Cassette.overdub(::ROCCtx, ::typeof(^), x::Float64, y::Float64) = AMDGPU.pow(x, y)
@inline Cassette.overdub(::ROCCtx, ::typeof(^), x::Float32, y::Float32) = AMDGPU.pow(x, y)
@inline Cassette.overdub(::ROCCtx, ::typeof(^), x::Float64, y::Int32)   = AMDGPU.pow(x, y)
@inline Cassette.overdub(::ROCCtx, ::typeof(^), x::Float32, y::Int32)   = AMDGPU.pow(x, y)
#@inline Cassette.overdub(::ROCCtx, ::typeof(^), x::Union{Float32, Float64}, y::Int64) = AMDGPU.pow(x, y)

# device lib functions
# TODO: Query these from AMDGPU.MATH_INTRINSICS
const rocfuns = (:cos, :cospi, :sin, :sinpi, :tan,
          :acos, :asin, :atan,
          :cosh, :sinh, :tanh,
          :acosh, :asinh, :atanh,
          :log, :log10, :log1p, :log2,
          :exp, :exp2, :exp10, :expm1, :ldexp,
          # :isfinite, :isinf, :isnan, :signbit,
          :abs,
          :sqrt, :cbrt,
          :ceil, :floor,)
for f in rocfuns
    @eval function Cassette.overdub(ctx::ROCCtx, ::typeof(Base.$f), x::Union{Float32, Float64})
        @Base._inline_meta
        return AMDGPU.$f(x)
    end
end

@inline Cassette.overdub(::ROCCtx, ::typeof(sincos), x::Union{Float32, Float64}) = (AMDGPU.sin(x), AMDGPU.cos(x)) # TODO: AMDGPU.sincos
# TODO: @inline Cassette.overdub(::ROCCtx, ::typeof(exp), x::Union{ComplexF32, ComplexF64}) = AMDGPU.exp(x)
# TODO: SpecialFunctions intrinsics

@inline Cassette.overdub(::ROCCtx, ::typeof(SpecialFunctions.gamma), x::Union{Float32, Float64}) = AMDGPU.tgamma(x)
@inline Cassette.overdub(::ROCCtx, ::typeof(SpecialFunctions.erf), x::Union{Float32, Float64}) = AMDGPU.erf(x)
@inline Cassette.overdub(::ROCCtx, ::typeof(SpecialFunctions.erfc), x::Union{Float32, Float64}) = AMDGPU.erfc(x)

import KernelAbstractions: ConstAdaptor, SharedMemory, Scratchpad, __synchronize, __size

###
# GPU implementation of shared memory
###
@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(SharedMemory), ::Type{T}, ::Val{Dims}, ::Val{Id}) where {T, Dims, Id}
    ptr = AMDGPU.alloc_special(Val(Id), T, Val(AMDGPU.AS.Local), Val(prod(Dims)))
    AMDGPU.ROCDeviceArray(Dims, ptr)
end

###
# GPU implementation of scratch memory
# - private memory for each workitem
###

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(Scratchpad), ::Type{T}, ::Val{Dims}) where {T, Dims}
    MArray{__size(Dims), T}(undef)
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__synchronize))
    AMDGPU.sync_workgroup()
end

@inline function Cassette.overdub(ctx::ROCCtx, ::typeof(__print), args...)
    for arg in args
        AMDGPU.@rocprintf("%s", arg)
    end
end

###
# GPU implementation of `@Const`
###
#=
struct ConstROCDeviceArray{T,N,A} <: AbstractArray{T,N}
    shape::Dims{N}
    ptr::Core.LLVMPtr{T,A}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    ConstROCDeviceArray{T,N,A}(shape::Dims{N}, ptr::Core.LLVMPtr{T,A}) where {T,A,N} = new(shape,ptr)
end

Adapt.adapt_storage(to::ConstAdaptor, a::AMDGPU.ROCDeviceArray{T,N,A}) where {T,N,A} = ConstROCDeviceArray{T, N, A}(a.shape, a.ptr)

Base.pointer(a::ConstROCDeviceArray) = a.ptr
Base.pointer(a::ConstROCDeviceArray, i::Integer) =
    pointer(a) + (i - 1) * Base.elsize(a)

Base.elsize(::Type{<:ConstROCDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::ConstROCDeviceArray) = g.shape
Base.length(g::ConstROCDeviceArray) = prod(g.shape)
Base.IndexStyle(::Type{<:ConstROCDeviceArray}) = Base.IndexLinear()

Base.unsafe_convert(::Type{Core.LLVMPtr{T,A}}, a::ConstROCDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

@inline function Base.getindex(A::ConstROCDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = Base.datatype_alignment(T)
    AMDGPU.unsafe_cached_load(pointer(A), index, Val(align))::T
end

@inline function Base.unsafe_view(arr::ConstROCDeviceArray{T, 1, A}, I::Vararg{Base.ViewIndex,1}) where {T, A}
    ptr = pointer(arr) + (I[1].start-1)*sizeof(T)
    len = I[1].stop - I[1].start + 1
    return ConstROCDeviceArray{T,1,A}(len, ptr)
end
=#

KernelAbstractions.argconvert(::Kernel{ROCDevice}, arg) = AMDGPU.rocconvert(arg)

end


text:
using DelimitedFiles
using WriteVTK
using MeshCore: AbsShapeDesc, SHAPE_DESC, P1, L2, T3, Q4, T4, H8, T6, Q8
using MeshCore: ShapeColl, shapedesc, nshapes, IncRel, nrelations
using MeshCore: VecAttrib, attribute
using LinearAlgebra: norm

include("import.jl")
include("export.jl")


text:
module SimpleParser

export  Parser,
        look_ahead_type,
        next_token,
        match,
        Lexer,
        start,
        next,
        done,
        show,
        Token,
        NUMBER,
        STRING,
        UNKNOWN,
        EOF,
        LPAREN,
        RPAREN,
        LBRACE,
        RBRACE,
        COMMA,
        EQUAL,
        SEMICOLON

import Base: ==, start, next, done, show, match, iteratorsize

include("lexer.jl")
include("parser.jl")

end


text:
module Server

import BufferedStreams
import Sockets
import ThreadPools

include("parsers/Http1_1.jl")
include("Utils.jl")
include("Writer.jl")

struct Error
    code::Int16
    message::String
end

struct Request
    target::String
    method::String
    headers::Dict{String, String}
    body::String
end

SERVER = nothing

"""
    _handle_connection(buffered_stream::BufferedStreams.BufferedInputStream) :: Tuple{Error, Request}

Forward BufferedInputStream to Http1.1 parser
If error occurs, return Error struct and empty Request struct
"""
function _handle_connection(buffered_stream::BufferedStreams.BufferedInputStream) :: Tuple{Error, Request}
    body::String = ""
    error::Error = Error(0, "")
    headers::Dict{String, String} = Dict()
    method::String = ""
    target::String = ""
    request::Request = Request(target, method, headers, body)

    try
        target, method, headers, body = ParserHttp1_1.parse_request(buffered_stream)
        request = Request(target, method, headers, body)
    catch e
        if method!="" && e isa EOFError
            request = Request(target, method, headers, body)
        elseif contains(string(typeof(e)), "TimeoutError")
            error = Error(408, "Request Timeout")
        elseif contains(string(typeof(e)), "InternalError")
            error = Error(500, sprint(showerror, e))
        else
            error = Error(401, "Bad Request")
        end
    end

    return error, request
end

"""
    close()

Stop listening on Socket
"""
function close()
    Sockets.close(SERVER)
end

"""
    start(f::Function, host::Union{Sockets.IPAddr, String}, port::Integer) :: Request

Listen to incoming requests on separate threads
Write response (execute corresponding Julia function) to Socket, before closing it
"""
function start(f::Function, host::Union{Sockets.IPAddr, String}, port::Integer) :: Request
    error::Error = Error(0, "")
    request::Request = Request("", "", Dict(), "")
    response::String = ""

    global SERVER
    SERVER = Sockets.listen(parse(Sockets.IPAddr, host), port)  # TODO increase backlog=511

    while isopen(SERVER)
        ThreadPools.@bthreads for x in 1:1
            socket = Sockets.accept(SERVER)
            bytes::Vector{UInt8} = readavailable(socket)

            # TODO: Find better way for PUT/POST requests with blank line before body content
            # As simply `bytesavailable` will return 0
            if ( bytes[1:4] == UInt8[0x50, 0x4f, 0x53, 0x54] || bytes[1:3] == UInt8[0x50, 0x55, 0x54] )
                if ( bytes[length(bytes)-3:length(bytes)] == UInt8[0x0d, 0x0a, 0x0d, 0x0a] )
                    bytes = append!(bytes, readavailable(socket))
                end
            end
            buffered_stream = BufferedStreams.BufferedInputStream(bytes)

            try
                error, request = _handle_connection(buffered_stream)
                if error.message != ""
                    response = Writer.write(error.code, request.headers; body=error.message)
                else
                    response = f(request)
                end
            catch e
                error = Error(500, sprint(showerror, e))
                @error "Error handling request" exception=(e, stacktrace(catch_backtrace()))
                if e isa Base.IOError && e.code == -54
                    @warn "Connection reset by peer (ECONNRESET)"
                else
                    @error exception=(e, stacktrace(catch_backtrace()))
                end
            finally
                if isopen(socket)
                    write(socket, response)
                    BufferedStreams.close(buffered_stream)
                    Sockets.close(socket)
                else
                    println("stream closed ...")
                end
            end
        end
    end

    return request
end

end


text:
using Combinatorics
using DataFrames
using DataStructures
using Dates
using LightGraphs
using LinearAlgebra
using Random
using StaticArrays
using Statistics
using TeaFiles
using Test
using TimeDag
using TimeZones

using TimeDag: Node
using TimeDag: IntersectAlignment, LeftAlignment, UnionAlignment
using TimeDag: duplicate, evaluate_until!, start_at

#! format: off

@testset "TimeDag.jl" begin
    # Perform these tests first, because they do naughty things to the global identity map
    # which would e.g. break nodes defined in common.jl.
    @testset "identity map" begin include("identity_map.jl") end

    include("common.jl")

    @testset "block" begin include("block.jl") end
    @testset "constants" begin include("constants.jl") end
    @testset "maybe" begin include("maybe.jl") end
    @testset "output_type" begin include("output_type.jl") end

    @testset "ops" begin
        @testset "align" begin include("ops/align.jl") end
        @testset "array" begin include("ops/array.jl") end
        @testset "common" begin include("ops/common.jl") end
        @testset "conditional" begin include("ops/conditional.jl") end
        @testset "history" begin include("ops/history.jl") end
        @testset "lagging" begin include("ops/lagging.jl") end
        @testset "random" begin include("ops/random.jl") end
        @testset "simple" begin include("ops/simple.jl") end
        @testset "sources" begin include("ops/sources.jl") end
        @testset "window" begin include("ops/window.jl") end
    end
end

#! format: on


text:
"""
    DimensionError <: Exception
    DimensionError(name, dim_expected, dim_found)

    Error for unexpected dimension.
    Output: "DimensionError: Input `name` should have length `dim_expected` not `dim_found`"

    ## From NLPModels.jl
"""
struct DimensionError <: Exception
    name::Union{Symbol,String}
    dim_expected::Int
    dim_found::Int
end

function Base.showerror(io::IO, e::DimensionError)
    print(
        io,
        "DimensionError: Input $(e.name) should have length $(e.dim_expected) not $(e.dim_found)",
    )
end

# https://groups.google.com/forum/?fromgroups=#!topic/julia-users/b6RbQ2amKzg
# https://github.com/JuliaSmoothOptimizers/NLPModels.jl
"""
    @lencheck n x y z …

    Check that arrays `x`, `y`, `z`, etc. have a prescribed length `n`.
"""
macro lencheck(l, vars...)
    exprs = Expr[]
    for var in vars
        varname = string(var)
        push!(exprs, :(
            if length($(esc(var))) != $(esc(l))
                throw(DimensionError($varname, $(esc(l)), length($(esc(var)))))
            end
        ))
    end
    Expr(:block, exprs...)
end

"""
    ContractError(abstract_type, behavior) <: Exception

    Occurs when a concrete type of `abstract_type` has not implemented a specified `behavior`
    Output: "SubtypingError: `behavior` is a requirement of `abstract_type` contract"
"""
struct ContractError <: Exception
    abstract_type::Type
    behavior::Symbol
end

function Base.showerror(io::IO, e::ContractError)
    print(
        io,
        "ContractError: $(e.behavior) is a requirement of $(e.abstract_type) contract!",
    )
end

"""
    @contract abstract_type :behavior

    Declare symbol `:behavior` is a hard-contract requirement of `abstract_type`.
"""
macro contract(abstract_type, behavior)
    return :(throw(ContractError($(esc(abstract_type)), $(esc(behavior)))))
end

"""
    AccessError(type) <: Exception

    Output: "AccessError: `:type` restricts field access
"""
struct AccessError <: Exception
    type::Type
end

function Base.showerror(io::IO, e::AccessError)
    print(io, "AccessError: $(e.type) restricts field access")
end

"""
@restrict T::Type

Restricts field access of type `T` by overloading `Base.getproperty` as:

```julia
Base.getproperty(T, s) = @restrict typeof(T) s
```

then accessing fields of `T` throws an `AccessError`.
"""
macro restrict(type)
    return :(throw(AccessError($(esc(type)))))
end

# """
# _parse_type(table::Dict{Symbol, Any}, T::Type, instance)

# Stores symbol, value pairs consisting of field symbols of `T` and the value
# stored at `instance<:T`.
# """
# function _parse_type!(table::Dict{Symbol, Any}, T::Type, instance)
#     if instance isa T
#         for field in fieldnames(T)
#             table[field] = getfield(m, field)
#         end
#     end
#     return nothing
# end


text:
"""
**Distance to a primary producer**

This function measures, for every species, its shortest path to a primary
producer using matrix exponentiation. A primary producer has a value of 1,
a primary consumer a value of 2, and so forth.

"""
function distance_to_producer(L::Array{Int64, 2})

    # We identify producers
    is_producer = vec(sum(L, dims = 2) .== 0)

    # Producers have a distance of 1
    d = zeros(Int64, length(is_producer))
    d[is_producer] .= 1

    # We work on a copy of the matrix with no self-loops
    K = copy(L)
    for i in eachindex(d)
        K[i ,i] = 0
    end

    # We loop as long as there are species with unknown distance
    i = 1
    while (i < length(d))|(sum(d) == 0)
        connected_at_length = (K^i * is_producer) .> 0
        d[(d .== 0) .* (connected_at_length)] .= i+1
        i = i+1
    end
    return d
end

"""
**Trophic rank**

Based on the average distance of preys to primary producers. Specifically, the
rank is defined as the average of the distance of preys to primary producers
(recursively). Primary producers always have a trophic rank of 1.

"""
function trophic_rank(L::Array{Int64, 2})
    # Average of positve elements, 0 otherwise
    nonzeromean = (x) -> maximum(x) == 0 ? 0 : mean(x[x.>0])
    d = distance_to_producer(L)
    dL = L .* d'
    TL = zeros(length(d))
    for i in eachindex(d)
        TL[i] = nonzeromean(dL[i,:]) + 1
    end
    return TL
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
function +(a,b)
    println(a,"+",b)
    return a-b
end

i=8+9

function emmm(x::Int=9)
    println(x+8)
end

emmm()
emmm(1)

max_stars_repo_name: Rratic/JuliaGuide

max_issues_repo_name: Rratic/JuliaGuide

max_forks_repo_name: Rratic/JuliaGuide


text:
# this file contains the functions to evaluate the right hand side of the weak form in the pdf

# Organization:
# The code is organized into 3 levels of functions.  High level functions
# only call mid level functions, passing the Mesh, Equation, and SBP objects
# to them.  At this point, there shouldn't be any need for more high level
# functions.  High level functions should use AbstractMesh and
# AbstractEulerData as the types of the mesh and equation objects, to signify
# that they are implementation independent.
#
# The mid level functions either call SBP functions or have loops over
# the arrays stored in the Equation, calling low level functions
# to calculate quantities at each node.  Use the ArrayViews package to pass
# portions of a larger array to the low level functions.  These functions
# should use AbstractMesh and EulerData as the types of the mesh and equation
# objects, to signify they are independent of the mesh implementation but
# do interact with the fields of the equation object.
#
# Low level functions calculate quantities for a single node.  Although
# these functions will likely take in an ArrayView, their argument
# types should be AbstractArrays{T, N}, specifying T according to
# what type the original array is (from either the Equation or Mesh object)
# They should also take the eqn.params object as an argument, which has
# all of the static parameters as the the equation object, so it can be
# used to dispatch to the correct method based on the static paremters.

# For high and mid level functions, the proper order for arguments
# is (mesh, sbp, eqn, opts).
# For low level functions, the order of arguments should be:  eqn.params first,
# then all input (read-only) arguments, then any output arguments.

# Rules:
# 1. function that take a composite type with abstract fields *cannot* use those
# fields directly (they can pass those fields to other functions however
# This leads directly to a two level structure for the code: high level function
# that take in composite types and low level function that take in arrays and
# perform calculations on them.

# The reason for this is that the compiler does not compile new version of the
#   function based on the types of the fields of a composite type. Passing the
#   fields of the type to other functions fixes this problem because the fields
#   are now arguments, so the compiler can specialize the code.

# 2.  Arrays should not be returned from functions.  The caller should allocate
# and array and pass it into the function.

# This allows reusing the same array during a loop (rather than
#   allocating a new array).

@doc """
### EulerEquationMod General Description
This module is organized into 3 levels of functions: high, middle, and low.

The high level functions take the mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerEquation,
and opts (options dictionary), They do not know the types their
arguments are paramaterized on. There is only one method for each high level function.
All they do is call mid level functions.

Mid level functions take the same arguments as high level functions but know
the types they are paramaterized on and do the right thing for all parameters.
There is only one method for each mid level function.  Mid level functions
typically pass pieces of the large arrays stored in the eqn and mesh objects
to low level functions that actually do computation.

Low level functions usually take in vectors of whatever values they need, do
a computation, and store the result in an argument vector.  These function
are usually passed an ArrayView, so writing to the argument vector writes
directly to the large array stored in the eqn object.  Thus there is no
loss of efficiency by using low level functions.
The Params type is also passed to low level function, which has all the
type parameters as the EulerEquation object, so it can be used for dispatch.
"""

# Rules for paramaterization:
# Tmsh = mesh data type
# Tsbp = SBP operator data type
# Tsol = equation, res_vec, q_vec data type

import PDESolver.evalResidual

@doc """
### EulerEquationMod.evalResidual

  This function drives the evaluation of the EulerEquations.
  It is agnostic to the dimension of the equation. and the types the arguments
  are paramaterized on.

  The function calls only high level functions, all of which take the same
  four arguments.  Mid level function also take the same arguments.

  The input/output variables are eqn.q and eqn.res, respectively.
  eqn.q_vec and eqn.res_vec exist for reusable storage *outside* the residual
  evaluation.  They should never be used inside the residual evaluation.

  The function array1DTo3D takes q_vec and puts it into eqn.q
  The function array3DTo1D takes eqn.res and puts it into res_vec

  Arguments:
    * mesh  : a mesh object
    * sbp   : SBP operator object
    * eqn   : an EulerData object
    * opts  : options dictionary

  The optional time argument is used for unsteady equations

"""->
# this function is what the timestepper calls
# high level function
function evalResidual(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData,
                     opts::Dict, t=0.0)

  time = eqn.params.time
  eqn.params.t = t  # record t to params
  myrank = mesh.myrank

  #println("\nentered evalResidual")
  time.t_send += @elapsed if opts["parallel_type"] == 1
    setParallelData(eqn.shared_data, opts["parallel_data"])
    startSolutionExchange(mesh, sbp, eqn, opts)
  end

  time.t_dataprep += @elapsed dataPrep(mesh, sbp, eqn, opts)
#  println("dataPrep @time printed above")

  time.t_volume += @elapsed if opts["addVolumeIntegrals"]
    evalVolumeIntegrals(mesh, sbp, eqn, opts)
    #println("volume integral @time printed above")
  end

  if opts["use_GLS"]
    GLS(mesh,sbp,eqn)
  end

  time.t_bndry += @elapsed if opts["addBoundaryIntegrals"]
    evalBoundaryIntegrals(mesh, sbp, eqn, opts)
    #println("boundary integral @time printed above")
  end

  time.t_stab += @elapsed if opts["addStabilization"]
    addStabilization(mesh, sbp, eqn, opts)
    #println("stabilizing @time printed above")
  end

  time.t_face += @elapsed if mesh.isDG && opts["addFaceIntegrals"]
    evalFaceIntegrals(mesh, sbp, eqn, opts)
    #println("face integral @time printed above")
  end

  time.t_sharedface += @elapsed if mesh.commsize > 1
    evalSharedFaceIntegrals(mesh, sbp, eqn, opts)
    #println("evalSharedFaceIntegrals @time printed above")
  end

  time.t_shock += @elapsed if opts["addShockCapturing"]
    evalShockCapturing(mesh, sbp, eqn, opts)
    #println("shock capturing @time printed above")
  end

  time.t_source += @elapsed evalSourceTerm(mesh, sbp, eqn, opts)
  #println("source integral @time printed above")

  # apply inverse mass matrix to eqn.res, necessary for CN
  if opts["use_Minv"]
    applyMassMatrixInverse3D(mesh, sbp, eqn, opts, eqn.res)
  end

  return nothing
end  # end evalResidual

@doc """
### EulerEquationMod.init

  This function performs any operations that need to be performed before
  the first residual evaluation.
  Any operations that need to be performed at the beginning of *each*
  residual evaluation should go in dataPrep

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * pmesh (equal to `mesh` by default)

  **Keyword Arguments**

   * init_mesh: set the boundary condition functors, default true,
                should false if `eqn`is nested inside another physics
                module's [`AbstractSolutionData`](@ref) object.
"""
# high level functions
function init(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
              eqn::AbstractEulerData{Tsol, Tres}, opts, pmesh=mesh; init_mesh=true) where {Tmsh, Tsol, Tres}

  # println("\nInitializing Euler module")

  # get BC functors
  if init_mesh
    getBCFunctors(mesh, sbp, eqn, opts)
    getBCFunctors(pmesh, sbp, eqn, opts)
  end

  getSRCFunctors(mesh, sbp, eqn, opts)
  if mesh.isDG
    getFluxFunctors(mesh, sbp, eqn, opts)
    getFaceElementFunctors(mesh, sbp, eqn, opts)
    getFluxFunctors_diff(mesh, sbp, eqn, opts)
  end

  if opts["use_staggered_grid"]
    mesh2 = mesh.mesh2
    sbp2 = mesh.sbp2

    if init_mesh
      getBCFunctors(mesh2, sbp2, eqn, opts)
    end

    getSRCFunctors(mesh2, sbp2, eqn, opts)
    if mesh2.isDG
      getFluxFunctors(mesh2, sbp2, eqn, opts)
      getFaceElementFunctors(mesh2, sbp2, eqn, opts)
    end
  end

  init_revm(mesh, sbp, eqn, opts, pmesh)
  init_revq(mesh, sbp, eqn, opts, pmesh)

  return nothing
end

function init_revm(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
              eqn::AbstractEulerData{Tsol, Tres}, opts, pmesh=mesh) where {Tmsh, Tsol, Tres}

  # Get functors for the boundary conditions
  getBCFunctors_revm(mesh, sbp, eqn, opts)
  getBCFunctors_revm(pmesh, sbp, eqn, opts)

  if mesh.isDG
    # For element interface flux
    getFluxFunctors_revm(mesh, sbp, eqn, opts)
  end

  return nothing
end

function init_revq(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
              eqn::AbstractEulerData{Tsol, Tres}, opts, pmesh=mesh) where {Tmsh, Tsol, Tres}

  # Get functors for the boundary conditions
  getBCFunctors_revq(mesh, sbp, eqn, opts)
  getBCFunctors_revq(pmesh, sbp, eqn, opts)

  if mesh.isDG
    # For element interface flux
    getFluxFunctors_revq(mesh, sbp, eqn, opts)
  end

  return nothing
end

function majorIterationCallback(itr::Integer,
                                mesh::AbstractMesh{Tmsh},
                                sbp::AbstractOperator,
                                eqn::EulerData{Tsol, Tres, Tdim}, opts, f::IO) where {Tmsh, Tsol, Tres, Tdim}

#  println("Performing major Iteration Callback")

  myrank = mesh.myrank
  output_freq = opts["output_freq"]::Int

#  println("eqn.q = \n", eqn.q)

  if opts["write_vis"] && (((itr % opts["output_freq"])) == 0 || itr == 1)
    vals = real(eqn.q_vec)  # remove unneded imaginary part

    #TODO: testing
    if opts["addShockCapturing"]
      writeShockSensorField(mesh, sbp, eqn, opts, getShockSensor(eqn.shock_capturing))
    end

    saveSolutionToMesh(mesh, vals)
    fname = string("solution_", itr)
    writeVisFiles(mesh, fname)

    if opts["write_vorticity_vis"]
      # write vorticity field
      new_field = vec(getVorticity(mesh, sbp, eqn, opts))
      saveSolutionToMesh(mesh, new_field)
      fname = string("solution_vorticity_", itr)
      writeVisFiles(mesh, fname)
    end

    vals = real(eqn.res_vec)  # remove unneded imaginary part
    saveSolutionToMesh(mesh, vals)
    fname = string("residual_", itr)
    writeVisFiles(mesh, fname)

#=
    # DEBUGGING: write error to file
    q_exact = zeros(eqn.q_vec)
    ex_func = ICDict[opts["IC_name"]]
    ex_func(mesh, sbp, eqn, opts, q_exact)
    q_err = real(eqn.q_vec) - q_exact
    saveSolutionToMesh(mesh, q_err)
    fname = string("error_", itr)
    writeVisFiles(mesh, fname)
=#
  end

  if opts["callback_write_qvec"] && (itr % opts["output_freq"] == 0)
    writeSolutionFiles(mesh, sbp, eqn, opts, "callback_q_vec_$itr")
  end

  if opts["write_entropy"]
    if mesh.isDG
      # undo multiplication by inverse mass matrix
      res_vec_orig = eqn.M.*copy(eqn.res_vec)
      res_orig = reshape(res_vec_orig, mesh.numDofPerNode,
                         mesh.numNodesPerElement, mesh.numEl)
    end

    @mpi_master f = eqn.file_dict[opts["write_entropy_fname"]]

    if(itr % opts["write_entropy_freq"] == 0)
      # calculate the entropy norm
      val = real(calcEntropyIntegral(mesh, sbp, eqn, opts, eqn.q_vec))

      # compute w^T * res_vec
      val2 = real(contractResEntropyVars(mesh, sbp, eqn, opts, eqn.q_vec, res_vec_orig))
      @mpi_master println(f, itr, " ", eqn.params.t, " ",  val, " ", val2)
    end

    @mpi_master if (itr % output_freq) == 0
      flush(f)
    end
  end  # end if write_entropy

  if opts["write_integralq"]
    integralq_vals = integrateQ(mesh, sbp, eqn, opts, eqn.q_vec)
    @mpi_master begin
      f = eqn.file_dict[opts["write_integralq_fname"]]
      print(f, itr, " ", eqn.params.t)
      for i=1:length(integralq_vals)
        print(f, " ", integralq_vals[i])
      end
      print(f, "\n")

      if (itr % output_freq) == 0
        flush(f)
      end
    end
  end  # end if write_integralq

  if opts["write_enstrophy"]
    @mpi_master f = eqn.file_dict[opts["write_enstrophy_fname"]]

    if (itr % opts["write_enstrophy_freq"]) == 0
      enstrophy = real(calcEnstrophy(mesh, sbp, eqn, opts, eqn.q))
      @mpi_master println(f, itr, " ", eqn.params.t, " ", enstrophy)
    end

    @mpi_master if (itr % opts["output_freq"]) == 0
      flush(f)
    end
  end

  if opts["write_kinetic_energy"]
    @mpi_master f = eqn.file_dict[opts["write_kinetic_energy_fname"]]

    if (itr % opts["write_kinetic_energy_freq"]) == 0
      kinetic_energy = real(calcKineticEnergy(mesh, sbp, eqn, opts, eqn.q_vec))
      @mpi_master println(f, itr, " ", eqn.params.t, " ", kinetic_energy)
    end

    @mpi_master if (itr % opts["output_freq"]) == 0
      flush(f)
    end
  end

  if opts["write_kinetic_energydt"]
    @mpi_master f = eqn.file_dict[opts["write_kinetic_energydt_fname"]]

    if (itr % opts["write_kinetic_energydt_freq"]) == 0
      kinetic_energydt = real(calcKineticEnergydt(mesh, sbp, eqn, opts, eqn.q_vec, eqn.res_vec))
      @mpi_master println(f, itr, " ", eqn.params.t, " ", kinetic_energydt)
    end

    @mpi_master if (itr % opts["output_freq"]) == 0
      flush(f)
    end
  end

  if opts["write_drag"]
    #=
    @mpi_master f = eqn.file_dict[opts["write_drag_fname"]]

    if (itr % opts["write_drag_freq"]) == 0
      objective = EulerEquationMod.createObjectiveFunctionalData(mesh, sbp, eqn, opts)
      drag = real(evalFunctional(mesh, sbp, eqn, opts, objective))
      @mpi_master println(f, itr, " ", drag)
    end

    @mpi_master if (itr % opts["output_freq"]) == 0
      flush(f)
    end
    =#
  end

  return nothing

end

@doc """
### EulerEquationMod.dataPrep

  This function calculates all the quantities that are stored over the entire
  mesh, and performs a few sanity check such as ensuring density and pressure
  are positive.

  This is a high level function
"""
# high level function
function dataPrep(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                   eqn::AbstractEulerData{Tsol, Tres}, opts) where {Tmsh, Tsol, Tres}
# gather up all the data needed to do vectorized operatinos on the mesh
# calculates all mesh wide quantities in eqn
# this is almost the exact list of everything we *shouldn't* be storing, but
# rather recalculating on the fly

  # apply filtering to input
  if eqn.params.use_filter
    applyFilter(mesh, sbp, eqn, eqn.q, opts)
  end

  # zero out res
  fill!(eqn.res, 0.0)
  fill!(eqn.res_edge, 0.0)

  getAuxVars(mesh, eqn)
#  println("  getAuxVars @time printed above")

  if opts["check_density"]
    checkDensity(eqn, mesh)
    # println("  checkDensity @time printed above")
  end

  if opts["check_pressure"]
    checkPressure(eqn, mesh)
    # println("  checkPressure @time printed above")
  end

  # calculate fluxes

  if opts["use_staggered_grid"]
    aux_vars = zeros(Tres, 1, mesh.mesh2.numNodesPerElement)
    for i=1:mesh.numEl
      qs = ro_sview(eqn.q, :, :, i)
      qf = sview(eqn.q_flux, :, :, i)
      interpolateElementStaggered(eqn.params, mesh, qs, aux_vars, qf)
    end
  end

  #TODO: consider setting precompute_volume_flux to false for type 2 volume
  #      integrals
  if opts["precompute_volume_flux"]
    getEulerFlux(mesh, sbp,  eqn, opts)
    # println("  getEulerFlux @time printed above")
  end

  if mesh.isDG
    if opts["precompute_q_face"] && opts["addFaceIntegrals"]
      interpolateFace(mesh, sbp, eqn, opts, eqn.q, eqn.q_face)
      # println("  interpolateFace @time printed above")
    end

    if opts["precompute_face_flux"] && opts["addFaceIntegrals"]
      calcFaceFlux(mesh, sbp, eqn, eqn.flux_func, mesh.interfaces, eqn.flux_face)
      #  println("  interpolateFace @time printed above")
    end
    if opts["precompute_q_bndry"] && opts["addBoundaryIntegrals"]
      interpolateBoundary(mesh, sbp, eqn, opts, eqn.q, eqn.q_bndry, eqn.aux_vars_bndry)
      # println("  interpolateFace @time printed above")
    end
  end

  if opts["precompute_boundary_flux"] && opts["addBoundaryIntegrals"]
    fill!(eqn.bndryflux, 0.0)
    getBCFluxes(mesh, sbp, eqn, opts)
    # println("  getBCFluxes @time printed above")
  end

  # is this needed for anything besides edge stabilization?
  if eqn.params.use_edgestab
    stabscale(mesh, sbp, eqn)
  end
#  println("  stabscale @time printed above")

  return nothing
end # end function dataPrep

#------------------------------------------------------------------------------
# sanity check functions
#------------------------------------------------------------------------------
@doc """
### EulerEquationMod.checkDensity

  This function checks that the density is positive. If not, an error is thrown.
  Because density is stored in the eqn.q array, performing this check takes
  very little time.

  Arguments:
    * EulerData
    * mesh

  This is a mid level function.
"""->
# mid level function
function checkDensity(eqn::EulerData{Tsol}, mesh) where Tsol
# check that density is positive

(ndof, nnodes, numel) = size(eqn.q)
q_cons = zeros(Tsol, ndof)  # conservative variables
for i=1:numel
  for j=1:nnodes
    convertToConservative(eqn.params, sview(eqn.q, :, j, i), q_cons)

    if real(q_cons[1]) <= 0.0
      println(STDERR, "Negative density at element ", i, ", node ", j)
      println(STDERR, "Coordinates = ", mesh.coords[:, j, i])
      println(STDERR, "q = ", q_cons)
      error("Negative density detected")
    end

    #=
    if real(q_cons[1]) < 0.0
      println("q_conservative = ", q_cons)
    end
    @assert( real(q_cons[1]) > 0.0, "element $i, node $j. Density < 0")
    =#
  end
end

return nothing

end

@doc """
### EulerEquationMod.checkPressure

  This function checks that pressure is positive.  If not, an error is thrown.
  Because pressure is stored in the auxiliary variables array, this check
  takes very little time.

  Arguments:
    * EulerData
    * mesh

  This is a mid level function
"""->
function checkPressure(eqn::EulerData, mesh)
# check that density is positive

(ndof, nnodes, numel) = size(eqn.q)

for i=1:numel
  for j=1:nnodes
    q = sview(eqn.q, :, j, i)
    aux_vars = sview(eqn.aux_vars,:, j, i)
    press = @getPressure(aux_vars)
    if real(press) <= 0.0
      println(STDERR, "Negative pressure at element ", i, ", node ", j)
      println(STDERR, "Coordinates = ", mesh.coords[:, j, i])
      println(STDERR, "q = ", q)
      println(STDERR, "press = ", press)
      error("Negative pressure detected")
    end
#    @assert( real(press) > 0.0, "element $i, node $j, q = $q, press = $press")
  end
end

return nothing

end

#------------------------------------------------------------------------------
# functions that evaluate terms in the weak form
#------------------------------------------------------------------------------
"""
### EulerEquationMod.evalVolumeIntegrals

  This function evaluates the volume integrals of the Euler equations by
  calling the appropriate SBP function on the Euler flux stored in
  eqn.flux_parametric.
  This function knows the dimension of the equation and does the right
  thing for all dimensions.  eqn.res is updated with the result

  This is a mid level function.
"""
function evalVolumeIntegrals(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                             opts) where {Tmsh,  Tsol, Tres, Tdim}

  integral_type = opts["volume_integral_type"]
  if integral_type == 1  # regular volume integrals

    if opts["precompute_volume_flux"]

      if opts["Q_transpose"] == true
        for i=1:Tdim
          weakdifferentiate!(sbp, i, sview(eqn.flux_parametric, :, :, :, i), eqn.res, trans=true)
        end
      else
        for i=1:Tdim
          weakdifferentiate!(sbp, i, sview(eqn.flux_parametric, :, :, :, i), eqn.res, SummationByParts.Subtract(), trans=false)
        end
      end  # end if Q_transpose
    else  # not precomputing the volume flux
      calcVolumeIntegrals_nopre(mesh, sbp, eqn, opts)
    end  # end if precompute_volume _flux

  elseif integral_type == 2  # entropy stable formulation
    calcVolumeIntegralsSplitForm(mesh, sbp, eqn, opts, eqn.volume_flux_func)
  else
    throw(ErrorException("Unsupported volume integral type = $integral_type"))
  end

  # artificialViscosity(mesh, sbp, eqn)

end  # end evalVolumeIntegrals

"""
  This function evaluates the boundary integrals in the Euler equations by
  calling the appropriate SBP function on eqn.bndryflux, which must be populated
  before calling this function.  eqn.res is updated with the result

  This is a mid level function

"""
function evalBoundaryIntegrals(mesh::AbstractMesh{Tmsh},
       sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh, Tsol, Tres, Tdim}

  #TODO: remove conditional
  if mesh.isDG
    if opts["precompute_boundary_flux"]
      boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, eqn.bndryflux, eqn.res, SummationByParts.Subtract())
    # else do nothing
    else
      # when precompute_boundary_flux == false, this fuunction does the
      # integration too, updating res
      getBCFluxes(mesh, sbp, eqn, opts)
    end

  else
    boundaryintegrate!(mesh.sbpface, mesh.bndryfaces, eqn.bndryflux, eqn.res, SummationByParts.Subtract())
  end

  return nothing

end  # end evalBoundaryIntegrals

@doc """
### EulerEquationMod.addStabilization

  This function add whatever form of stabilization opts specifies by calling
  the appropriate function.  Some arrays may need to be populated before
  calling this function, depending on the type of stabilization used.

  This is a mid level function
"""->
# mid level function
function addStabilization(mesh::AbstractMesh{Tmsh},
             sbp::AbstractOperator, eqn::EulerData{Tsol}, opts) where {Tmsh,  Tsol}

#  println("==== start of addStabilization ====")

  # ----- Edge Stabilization -----#
  if eqn.params.use_edgestab
#    println("applying edge stabilization")
    if opts["use_edge_res"]
      edgestabilize!(mesh, sbp, eqn, mesh.interfaces, eqn.q, mesh.coords,
                     mesh.dxidx, mesh.jac, eqn.edgestab_alpha, eqn.stabscale,
                     eqn.res, eqn.res_edge)
    else
      edgestabilize!(sbp, mesh.interfaces, mesh, eqn.q, mesh.coords, mesh.dxidx,
                     mesh.jac, eqn.edgestab_alpha, eqn.stabscale, eqn.res)
    end
  end

  # ----- Filtering -----
  if eqn.params.use_res_filter
#    println("applying residual filter")
    applyFilter(mesh, sbp, eqn, eqn.res, opts, trans=true)
  end

  # ----- Artificial Dissipation -----
  if eqn.params.use_dissipation
#    println("applying artificial dissipation")
    applyDissipation(mesh, sbp, eqn, opts, eqn.q)
  end

  if opts["use_GLS2"]
     applyGLS3(mesh, sbp, eqn, opts)
#    test_GLS(mesh, sbp, eqn, opts)
  end

  if opts["use_lps"]
    applyLPStab(mesh, sbp, eqn, opts)
  end
#  println("==== end of addStabilization ====")

  return nothing

end

@doc """
### EulerEquationMod.evalFaceIntegrals

  This function evaluates the face integrals in a DG formulation and
  updates the residual.  The array eqn.flux_face must already be populated
  with the face flux.

  **Inputs**:

   * mesh: an AbstractDGMesh
   * sbp: an SBP operator
   * eqn: an EulerData object
   * opts: the options dictonary

  **Outputs**:

    none

"""->
function evalFaceIntegrals(mesh::AbstractDGMesh{Tmsh},
               sbp::AbstractOperator, eqn::EulerData{Tsol}, opts) where {Tmsh, Tsol}

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1
#    println("calculating regular face integrals")
    if opts["precompute_face_flux"]
      interiorfaceintegrate!(mesh.sbpface, mesh.interfaces, eqn.flux_face,
                             eqn.res, SummationByParts.Subtract())
    else
      calcFaceIntegral_nopre(mesh, sbp, eqn, opts, eqn.flux_func, mesh.interfaces)
    end

  elseif face_integral_type == 2
#    println("calculating ESS face integrals")
    if opts["use_staggered_grid"]
      getFaceElementIntegral(mesh, mesh.mesh2, sbp, mesh.sbp2, eqn,
                             eqn.face_element_integral_func,
                             eqn.flux_func, mesh.mesh2.sbpface, mesh.interfaces)

    else
      getFaceElementIntegral(mesh, sbp, eqn, eqn.face_element_integral_func,
                             eqn.flux_func, mesh.sbpface, mesh.interfaces)
    end

  else
    throw(ErrorException("Unsupported face integral type = $face_integral_type"))
  end

  # do some output here?
  return nothing

end

#=
@doc """
### EulerEquationMod.sendParallelData

  This function interpolates the data into the send buffer and post
  the Isends and Irecvs.  It does not wait for them to finish

  Inputs:
    mesh
    sbp
    eqn
    opts
"""->
function sendParallelData(mesh::AbstractDGMesh, sbp, eqn, opts)

  for i=1:mesh.npeers
    # interpolate
    mesh.send_waited[i] = getSendData(mesh, opts, eqn.q, mesh.bndries_local[i], eqn.q_face_send[i], mesh.send_reqs[i], mesh.send_waited[i])
  end

  exchangeFaceData(mesh, opts, eqn.q_face_send, eqn.q_face_recv)

  return nothing
end
=#
@doc """
### EulerEquationMod.evalSharedFaceIntegrals

  This function does the computation that needs the parallel
  communication to have finished already, namely the face integrals
  for the shared faces

  **Inputs**:

   * mesh
   * sbp
   * eqn
   * opts
"""->
function evalSharedFaceIntegrals(mesh::AbstractDGMesh, sbp, eqn, opts)

  face_integral_type = opts["face_integral_type"]
  if face_integral_type == 1

    if opts["parallel_data"] == PARALLEL_DATA_FACE
      finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceIntegrals)
    elseif opts["parallel_data"] == PARALLEL_DATA_ELEMENT

      finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceIntegrals_element)
#      calcSharedFaceIntegrals_element(mesh, sbp, eqn, opts, eqn.flux_func)
    else
      throw(ErrorException("unsupported parallel data type"))
    end

  elseif face_integral_type == 2

      finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, calcSharedFaceElementIntegrals_element)
#    getSharedFaceElementIntegrals_element(mesh, sbp, eqn, opts, eqn.face_element_integral_func,  eqn.flux_func)
  else
    throw(ErrorException("unsupported face integral type = $face_integral_type"))
  end

  return nothing
end

"""
### EulerEquationMod.evalSourceTerm

  This function performs all the actions necessary to update eqn.res
  with the source term.  The source term is stored in eqn.src_func.  It is
  an abstract field, so it cannot be accessed (performantly) directly, so
  it is passed to an inner function.

  **Inputs**:

   * mesh : Abstract mesh type
   * sbp  : Summation-by-parts operator
   * eqn  : Euler equation object
   * opts : options dictonary

  Outputs: none

  Aliasing restrictions: none

"""
function evalSourceTerm(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  # placeholder for multiple source term functionality (similar to how
  # boundary conditions are done)
  if opts["use_src_term"]
    applySourceTerm(mesh, sbp, eqn, opts, eqn.src_func)
  end

  return nothing
end  # end function

"""
  This function adds the shock capturing terms (if requested) to the residual.

  **Inputs**:

   * mesh : Abstract mesh type
   * sbp  : Summation-by-parts operator
   * eqn  : Euler equation object
   * opts : options dictonary

  Outputs: none

  Aliasing restrictions: none
"""
function evalShockCapturing(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                     opts) where {Tmsh, Tsol, Tres, Tdim}

  # currently there is only one choice for the shock capturing scheme.
  # If there are more, need something more sophisiticated to choose.
  # Perhaps add an abstract typed field to eqn containing the structs

  capture = eqn.shock_capturing
  applyShockCapturing(mesh, sbp, eqn, opts, capture)

  return nothing
end

#TODO: put this in Utils?
@doc """
### EulerEquationMod.applyMassMatrixInverse3D

  This function applies the 3D inverse mass matrix to an array.
    The array passed in should always be eqn.res

  **Inputs**:

   * mesh: mesh object, needed for numEl and numDofPerNode fields
   * sbp: sbp object, needed for numnodes field
   * eqn: equation object, needed for Minv3D field
   * opts
   * arr: the 3D array to have the 3D mass matrix inverse applied to it

"""->
function applyMassMatrixInverse3D(mesh, sbp, eqn, opts, arr)

  for i = 1:mesh.numEl
    for j = 1:sbp.numnodes
      for k = 1:mesh.numDofPerNode
        arr[k, j, i] = eqn.Minv3D[k, j, i] * arr[k, j, i]
      end
    end
  end

  return arr
end


text:
using Dates
using Printf

import Base: show, +, -

"Clock that displays 24 hour clock that rolls over properly"
struct Clock{T<:Integer}
    h::T
    min::T

    function Clock(h::T, min::T) where {T}
        h += fld(min, 60)
        h %= 24
        min %= 60
        if h < 0
            h += 24
        end
        if min < 0
            min += 60
        end
        new{T}(h, min)
    end
end

show(io::IO, c::Clock) = show(io, @sprintf("%02d:%02d", c.h, c.min))

+(c::Clock, min::Dates.Minute) = Clock(c.h, c.min + min.value)
-(c::Clock, min::Dates.Minute) = Clock(c.h, c.min - min.value)


text:
suite["ancillaries"] = BenchmarkGroup(["identity", "field", "vertical", "left", "right"])
Ls   = buildLs(A, H; env_maxdim=χ, maxdim=D, cutoff=0.0)
Rs   = buildRs(A, H; env_maxdim=χ, maxdim=D, cutoff=0.0)
dummy = cuITensor(1.0)
suite["ancillaries"]["identity"] = @benchmarkable PEPS.makeAncillaryIs($A, $Ls[$col-1], $Rs[$col+1], $col)
vH = PEPS.getDirectional(vcat(H[:, col]...), PEPS.Vertical)
suite["ancillaries"]["vertical"] = @benchmarkable PEPS.makeAncillaryVs($A, $Ls[$col-1], $Rs[$col+1], $vH, $col)
fH = PEPS.getDirectional(vcat(H[:, col]...), PEPS.Field)
suite["ancillaries"]["field"]    = @benchmarkable PEPS.makeAncillaryFs($A, $Ls[$col-1], $Rs[$col+1], $fH, $col)
lH = PEPS.getDirectional(vcat(H[:, col-1]...), PEPS.Horizontal)
suite["ancillaries"]["left"]     = @benchmarkable PEPS.makeAncillarySide($A, $Ls[$col-1], $Rs[$col+1], $lH, $col, :left)
rH = PEPS.getDirectional(vcat(H[:, col]...), PEPS.Horizontal)
suite["ancillaries"]["right"]    = @benchmarkable PEPS.makeAncillarySide($A, $Rs[$col+1], $Ls[$col-1], $rH, $col, :right)


text:
import  NLsolve

using TimerOutputs
using Test
using Modia: SimulationOptions, emptyResult!, derivatives!, timeEventCondition!, affectTimeEvent!, outputs!,  reinitEventHandler, propagateEvaluateAndInstantiate!, eventIteration!, get_xe, hasParticles, resizeLinearEquations!, updateEquationInfo!

using Modia
using DataFrames

#***************************************************************************************
#***************************************************************************************
#***************************************************************************************

printstyled("!!override!!  function Modia.init!  -- VECTOR 'SIM8.jl'\n",color=:light_yellow)

#NEU
function Modia.init!(m::SimulationModel{FloatType,TimeType})::Bool where {FloatType,TimeType}
    emptyResult!(m)
    eh = m.eventHandler
    reinitEventHandler(eh, m.options.stopTime, m.options.logEvents)
    eh.firstInitialOfAllSegments = true

	# Apply updates from merge Map and propagate/instantiate/evaluate the resulting evaluatedParameters
    if length(m.options.merge) > 0
        m.parameters = mergeModels(m.parameters, m.options.merge)
        m.evaluatedParameters = propagateEvaluateAndInstantiate!(FloatType, m.unitless, m.modelModule, m.parameters, m.equationInfo, m.previous_dict, m.previous, m.pre_dict, m.pre, m.hold_dict, m.hold)
        if isnothing(m.evaluatedParameters)
            return false
        end

        # Resize linear equation systems if dimensions of vector valued tearing variables changed
        resizeLinearEquations!(m, m.options.log)

        # Resize state vector memory
        m.x_start = updateEquationInfo!(m.equationInfo, FloatType)
        nx = length(m.x_start)
        resize!(m.x_init, nx)
        resize!(m.der_x, nx)
        eqInfo = m.equationInfo
        m.x_vec = [zeros(FloatType, eqInfo.x_info[i].length) for i in eqInfo.nxFixedLength+1:length(eqInfo.x_info)]
    end

    # Initialize auxiliary arrays for event iteration
    m.x_init .= 0
    m.der_x  .= 0

    # Log parameters
    if m.options.logParameters
        parameters = m.parameters
        @showModel parameters
    end
    if m.options.logEvaluatedParameters
        evaluatedParameters = m.evaluatedParameters
        @showModel evaluatedParameters
    end

    if m.options.logStates
        # List init/start values
        x_table = DataFrames.DataFrame(state=String[], init=Any[], unit=String[], nominal=Float64[])
        for xe_info in m.equationInfo.x_info
            xe_init = get_xe(m.x_start, xe_info)
            if hasParticles(xe_init)
                xe_init = string(minimum(xe_init)) * " .. " * string(maximum(xe_init))
            end
            push!(x_table, (xe_info.x_name, xe_init, xe_info.unit, xe_info.nominal))
        end
        show(stdout, x_table; allrows=true, allcols=true, rowlabel = Symbol("#"), summary=false, eltypes=false)
        println("\n")
    end

    # Initialize model, linearEquations and compute and store all variables at the initial time
    if m.options.log
        println("      Initialization at time = ", m.options.startTime, " s")
    end

    # Perform initial event iteration
    m.nGetDerivatives = 0
    m.nf = 0
    m.isInitial   = true
    eh.initial    = true
#    m.storeResult = true
#    m.getDerivatives!(m.der_x, m.x_start, m, startTime)
#    Base.invokelatest(m.getDerivatives!, m.der_x, m.x_start, m, startTime)
    for i in eachindex(m.x_init)
        m.x_init[i] = deepcopy(m.x_start[i])
    end

##------NEUNEU------------------------------------------------------------------------------

ListeEIN = Dict{String, Any}()
ListeAUS = Dict{String, Any}()

for info in m.equationInfo.x_info
    ListeEIN[info.x_name] = info
    ListeAUS[info.der_x_name] = info
end

startValue = Vector{Float64}()
inputIndex = Vector{Any}()

logINITvalues = true

function report(txt::String, color::Symbol)
    if logINITvalues
        printstyled(txt, color=color)
    end
end

report("\nINPUT values\n", :blue)
for ein in EIN
    if haskey(ListeEIN, ein[1])
        rec = ListeEIN[ein[1]]
        if rec.scalar
            report("\t STATE scalar : $(ein[1])\n", :green)
            push!(inputIndex, rec.startIndex)
            push!(startValue, ein[2])
        else
            report("\t STATE vector : $(ein[1]) \n", :light_green)
            for i = 1:rec.length
                push!(inputIndex, rec.startIndex+i-1)
                push!(startValue, ein[2][i])   # keine Prüfung der Länge
            end
        end
    elseif haskey(m.parameters, Symbol(ein[1]))
        report("\t         PARA : $(ein[1])\n", :cyan)
        push!(startValue, ein[2])
        push!(inputIndex, Symbol(ein[1]))
    else
        printstyled("\tnicht gefunden: $(ein[1])\n", color=:red)
    end
end

if logINITvalues
    @show length(inputIndex)
    @show inputIndex
    @show startValue
end

targetIndex = Vector{Int64}()
targetValue = Vector{Float64}()

#TODO
#TODO Überprüfung der Länge bei vec
#TODO

#INFO  INdex negativ für OUTPUT target

report("\nTARGET values\n", :blue)
for aus in AUS
    if haskey(ListeAUS, aus[1])
        rec = ListeAUS[aus[1]]
        if rec.scalar
            report("\t STATE scalar : $(aus[1])\n", :green)
            push!(targetIndex, rec.startIndex)
            push!(targetValue, aus[2])
        else
            report("\t STATE vector : $(aus[1])\n", :light_green)
            for i = 1:rec.length
                push!(targetIndex, rec.startIndex+i-1)
                push!(targetValue, aus[2][i])   # keine Prüfung der Länge
            end
        end
    elseif haskey(m.result_info, aus[1])
        report("\t       RESULT : $(aus[1])\n", :cyan)
        push!(targetIndex, -m.result_info[aus[1]].index)
        push!(targetValue, aus[2])
    else
        printstyled("\tnicht gefunden: $(aus[1])\n", color=:red)
    end

end

if logINITvalues
    @show length(targetIndex)
    @show targetIndex
    @show targetValue
end

function fromSolver2State!(m, Xsolver, Xstate)

    for (i,w) in enumerate(inputIndex)
        if typeof(w) == Symbol
            m.evaluatedParameters[w] = Xsolver[i]
        else
            Xstate[w] = Xsolver[i]
        end
    end
end

function fromState2Solver!(m, Xstate,resi)

    allVars= pop!(m.init_vars)
    for (i,t) in enumerate(targetIndex)
        if t>0 # state
            resi[i] = Xstate[t] - targetValue[i]
        else
            V = allVars[-t]  # Zielwert
            resi[i] = V - targetValue[i]
        end
    end
end

#------------------------------

        TSTART = m.options.startTime

        solver0 = Vector{Float64}(startValue)
        state0 = vec(deepcopy(m.x_init))
        fromSolver2State!(m, solver0, state0)
        resid = similar(state0)
        derivatives!(resid,state0,m,TSTART)

        Lu = []
        Ldu = []
        stateU = Array{Float64,1}(state0)
        statedU =Array{Float64,1}(zeros(length(m.x_init)))

        function F!(residual,solverU)
            fromSolver2State!(m, solverU, stateU)
            derivatives!(statedU,stateU,m, TSTART)
            fromState2Solver!(m, statedU, residual)
        end

        m.storeInit = true
        RES = NLsolve.nlsolve(F!,solver0, iterations=100, xtol = 0.0, ftol=1e-8, method=:newton, show_trace=false)
        if !NLsolve.converged(RES)
            @show RES
            error("!NLsolve.converged(RES)")
        end

        if logINITvalues
            @show RES
        else
            printstyled("Initialization: Function Calls (f): $(RES.f_calls) \n",color=:light_blue)
        end

        fromSolver2State!(m, RES.zero,stateU)
        m.storeInit = false

        m.nf = 0    # reset
        m.x_init = copy(stateU)

##-Ende -----NEUNEU--------------------------------------------------------------------------

    eventIteration!(m, m.x_init, m.options.startTime)
    m.success     = false   # is set to true at the first outputs! call.
    eh.initial    = false
    m.isInitial   = false
    m.storeResult = false
    eh.afterSimulationStart = true
    return true
end

#***************************************************************************************
#***************************************************************************************
#***************************************************************************************

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
import Base.==
import Base.!=

const null = Nullable{Union{}}()

# C# spec section 7.10.9

=={T}(a::Nullable{T},b::Nullable{Union{}}) = isnull(a)
=={T}(a::Nullable{Union{}},b::Nullable{T}) = isnull(b)
!={T}(a::Nullable{T},b::Nullable{Union{}}) = !isnull(a)
!={T}(a::Nullable{Union{}},b::Nullable{T}) = !isnull(b)

# Strings

for op in (:lowercase,:uppercase,:reverse,:ucfirst,:lcfirst,:chop,:chomp)
    @eval begin
        import Base.$(op)
        function $op{T<:AbstractString}(x::Nullable{T})
            if isnull(x)
                return Nullable{T}()
            else
                return Nullable($op(get(x)))
            end
        end
    end
end

import Base.getindex
function getindex{T<:AbstractString}(s::Nullable{T},i)
    if isnull(s)
        return Nullable{T}()
    else
        return Nullable(get(s)[i])
    end
end

import Base.endof
function endof{T<:AbstractString}(s::Nullable{T})
    if isnull(s)
        # TODO Decide whether this makes sense?
        return 0
    else
        return endof(get(s))
    end
end

import Base.length
function length{T<:AbstractString}(s::Nullable{T})
    if isnull(s)
        return Nullable{Int}()
    else
        return Nullable{Int}(length(get(s)))
    end
end

# C# spec section 7.3.7

for op in (:+, :-, :!, :~)
    @eval begin
        import Base.$(op)
        $op{T<:Number}(x::Nullable{T}) = isnull(x) ? Nullable{T}() : Nullable($op(get(x)))
    end
end

for op in (:+, :-, :*, :/, :%, :&, :|, :^, :<<, :>>)
    @eval begin
        import Base.$(op)
        $op{T<:Number}(a::Nullable{T},b::Nullable{T}) = isnull(a) || isnull(b) ? Nullable{T}() : Nullable($op(get(a), get(b)))
        $op{T1<:Number,T2<:Number}(x::Nullable{T1},y::T2) = isnull(x) ? Nullable{promote_type(T1,T2)}() : Nullable{promote_type(T1,T2)}($op(get(x), y))
        $op{T1<:Number,T2<:Number}(x::T1,y::Nullable{T2}) = isnull(y) ? Nullable{promote_type(T1,T2)}() : Nullable{promote_type(T1,T2)}($op(x, get(y)))
    end
end

=={T1,T2}(a::Nullable{T1},b::Nullable{T2}) = isnull(a) && isnull(b) ? true : !isnull(a) && !isnull(b) ? get(a)==get(b) : false
=={T1,T2}(a::Nullable{T1},b::T2) = isnull(a) ? false : get(a)==b
=={T1,T2}(a::T1,b::Nullable{T2}) = isnull(b) ? false : a==get(b)

!={T1,T2}(a::Nullable{T1},b::Nullable{T2}) = isnull(a) && isnull(b) ? false : !isnull(a) && !isnull(b) ? get(a)!=get(b) : true
!={T1,T2}(a::Nullable{T1},b::T2) = isnull(a) ? true : get(a)!=b
!={T1,T2}(a::T1,b::Nullable{T2}) = isnull(b) ? true : a!=get(b)

for op in (:<,:>,:<=,:>=)
    @eval begin
        import Base.$(op)
        $op{T<:Number}(a::Nullable{T},b::Nullable{T}) = isnull(a) || isnull(b) ? false : $op(get(a), get(b))
        $op{T1<:Number,T2<:Number}(x::Nullable{T1},y::T2) = isnull(x) ? false : $op(get(x), y)
        $op{T1<:Number,T2<:Number}(x::T1,y::Nullable{T2}) = isnull(y) ? false : $op(x, get(y))
    end
end

# C# spec 7.11.4
function (&)(x::Nullable{Bool},y::Nullable{Bool})
    if isnull(x)
        if isnull(y) || get(y)==true
            return Nullable{Bool}()
        else
            return Nullable(false)
        end
    elseif get(x)==true
        return y
    else
        return Nullable(false)
    end
end

(&)(x::Bool,y::Nullable{Bool}) = x ? y : Nullable(false)
(&)(x::Nullable{Bool},y::Bool) = y ? x : Nullable(false)

function (|)(x::Nullable{Bool},y::Nullable{Bool})
    if isnull(x)
        if isnull(y) || !get(y)
            return Nullable{Bool}()
        else
            return Nullable(true)
        end
    elseif get(x)
        return Nullable(true)
    else
        return y
    end
end

(|)(x::Bool,y::Nullable{Bool}) = x ? Nullable(true) : y
(|)(x::Nullable{Bool},y::Bool) = y ? Nullable(true) : x