text:
"""
`module  JAC.Cascade
    ... a submodel of JAC that contains all methods to set-up and process (simple) atomic cascade computations.
"""
module Cascade

    using Printf, ..AutoIonization, ..Basics, ..Defaults, ..DecayYield, ..Radial, ..ManyElectron, ..Nuclear, ..PhotoEmission,
                  ..PhotoIonization, ..TableStrings

    """
    `abstract type Cascade.AbstractApproach`
        ... defines an abstract and a number of singleton types for the computational approach/model that is applied to
            generate and evaluate a cascade.

      + struct AverageSCA
        ... to evaluate the level structure and transitions of all involved levels in single-configuration approach but
            without configuration interaction and from just a single orbital set.

      + struct SCA
        ... to evaluate the level structure and transitions of all involved levels in single-configuration approach but
            by calculating all fine-structure resolved transitions.
    """
    abstract type  AbstractApproach                   end
    struct         AverageSCA  <:  AbstractApproach   end
    struct         SCA         <:  AbstractApproach   end
    struct         UserMCA     <:  AbstractApproach   end

    """
    `struct  Cascade.Block`
        ... defines a type for an individual block of configurations that are treatet together within the cascade. Such an block is given
            by a list of configurations that may occur as initial- or final-state configurations in some step of the canscade and that are
            treated together in a single multiplet to allow for configuration interaction but to avoid 'double counting' of individual levels.

        + NoElectrons     ::Int64                     ... Number of electrons in this block.
        + confs           ::Array{Configuration,1}    ... List of one or several configurations that define the multiplet.
        + hasMultiplet    ::Bool                      ... true if the (level representation in the) multiplet has already been computed and
                                                          false otherwise.
        + multiplet       ::Multiplet                 ... Multiplet of the this block.
    """
    struct  Block
        NoElectrons       ::Int64
        confs             ::Array{Configuration,1}
        hasMultiplet      ::Bool
        multiplet         ::Multiplet
    end

    """
    `Cascade.Block()`  ... constructor for an 'empty' instance of a Cascade.Block.
    """
    function Block()
        Block( )
    end

    # `Base.show(io::IO, block::Cascade.Block)`  ... prepares a proper printout of the variable block::Cascade.Block.
    function Base.show(io::IO, block::Cascade.Block)
        println(io, "NoElectrons:        $(block.NoElectrons)  ")
        println(io, "confs :             $(block.confs )  ")
        println(io, "hasMultiplet:       $(block.hasMultiplet)  ")
        println(io, "multiplet:          $(block.multiplet)  ")
    end

    """
    `struct  Cascade.Step`
        ... defines a type for an individual step of an excitation and/or decay cascade. Such an individual step is given by a well-defined
            process, such as Auger, PhotoEmission, or others and two lists of initial- and final-state configuration that are (each) treated
            together in a multiplet to allow for configuration interaction but to avoid 'double counting' of individual levels.

        + process          ::JBasics.AtomicProcess         ... Atomic process that 'acts' in this step of the cascade.
        + settings         ::Union{PhotoEmission.Settings, AutoIonization.Settings}
                                                       ... Settings for this step of the cascade.
        + initialConfs     ::Array{Configuration,1}    ... List of one or several configurations that define the initial-state multiplet.
        + finalConfs       ::Array{Configuration,1}    ... List of one or several configurations that define the final-state multiplet.
        + initialMultiplet ::Multiplet                 ... Multiplet of the initial-state levels of this step of the cascade.
        + finalMultiplet   ::Multiplet                 ... Multiplet of the final-state levels of this step of the cascade.
    """
    struct  Step
        process            ::Basics.AtomicProcess
        settings           ::Union{PhotoEmission.Settings, AutoIonization.Settings}
        initialConfs       ::Array{Configuration,1}
        finalConfs         ::Array{Configuration,1}
        initialMultiplet   ::Multiplet
        finalMultiplet     ::Multiplet
    end

    """
    `Cascade.Step()`  ... constructor for an 'empty' instance of a Cascade.Step.
    """
    function Step()
        Step( Basics.NoProcess, PhotoEmission.Settings, Configuration[], Configuration[], Multiplet(), Multiplet())
    end

    # `Base.show(io::IO, step::Cascade.Step)`  ... prepares a proper printout of the variable step::Cascade.Step.
    function Base.show(io::IO, step::Cascade.Step)
        println(io, "process:                $(step.process)  ")
        println(io, "settings:               $(step.settings)  ")
        println(io, "initialConfs:           $(step.initialConfs)  ")
        println(io, "finalConfs:             $(step.finalConfs)  ")
        println(io, "initialMultiplet :      $(step.initialMultiplet )  ")
        println(io, "finalMultiplet:         $(step.finalMultiplet)  ")
    end

    """
    `struct  Cascade.Computation`
        ... defines a type for a cascade computation, i.e. for the computation of a whole excitation and/or decay cascade. The data
            from this computation can be modified, adapted and refined to the practical needs before the actual computations are
            carried out. Initially, this struct contains the physical metadata about the cascade to be calculated but gets enlarged
            in course of the computation to keep also wave functions, level multiplets, etc.

        + name               ::String                         ... A name for the cascade
        + nuclearModel       ::Nuclear.Model                  ... Model, charge and parameters of the nucleus.
        + grid               ::Radial.Grid                    ... The radial grid to be used for the computation.
        + asfSettings        ::AsfSettings                    ... Provides the settings for the SCF process.
        + approach           ::Cascade.AbstractApproach       ... Computational approach/model that is applied to generate and evaluate the
                                                                  cascade; possible approaches are: {'single-configuration', ...}
        + processes          ::Array{Basics.AtomicProcess,1}   ... List of the atomic processes that are supported and should be included into the
                                                                  cascade.
        + initialConfs       ::Array{Configuration,1}         ... List of one or several configurations that contain the level(s) from which the
                                                                  cascade starts.
        + initialLevels      ::Array{Tuple{Int64,Float64},1}  ... List of one or several (tupels of) levels together with their relative population
                                                                  from which the cascade starts.
        + maxElectronLoss    ::Int64                          ... (Maximum) Number of electrons in which the initial- and final-state
                                                                  configurations can differ from each other; this also determines the maximal steps
                                                                  of any particular decay path.
        + NoShakeDisplacements ::Int64                        ... Maximum number of electron displacements due to shake-up or shake-down processes
                                                                  in any individual step of the cascade.
        + shakeFromShells ::Array{Shell,1}                    ... List of shells from which shake transitions may occur.
        + shakeToShells   ::Array{Shell,1}                    ... List of shells into which shake transitions may occur.
        + steps           ::Array{Cascade.Step,1}             ... List of individual steps between well-defined atomic multiplets that are
                                                                  included into the cascade.
    """
    struct  Computation
        name                 ::String
        nuclearModel         ::Nuclear.Model
        grid                 ::Radial.Grid
        asfSettings          ::AsfSettings
        approach             ::Union{Cascade.AbstractApproach}
        processes            ::Array{Basics.AtomicProcess,1}
        initialConfs         ::Array{Configuration,1}
        initialLevels        ::Array{Tuple{Int64,Float64},1}
        maxElectronLoss      ::Int64
        NoShakeDisplacements ::Int64
        shakeFromShells      ::Array{Shell,1}
        shakeToShells        ::Array{Shell,1}
        steps                ::Array{Cascade.Step,1}
    end

    """
    `Cascade.Computation()`  ... constructor for an 'empty' instance of a Cascade.Computation.
    """
    function Computation()
        Computation("",  Nuclear.Model(0.), Radial.Grid(), AsfSettings(), averageSCA, [Auger],
                    Configuration[], [(0, 0.)], 0, 0, Shell[], Shell[], CascadeComputationStep[] )
    end

    # `Base.show(io::IO, computation::Cascade.Computation)`  ... prepares a proper printout of the variable computation::Cascade.Computation.
    function Base.show(io::IO, computation::Cascade.Computation)
        println(io, "name:                     $(computation.name)  ")
        println(io, "nuclearModel:             $(computation.nuclearModel)  ")
        println(io, "grid:                     $(computation.grid)  ")
        println(io, "asfSettings:                computation.asfSettings  ")
        println(io, "approach:                 $(computation.approach)  ")
        println(io, "processes:                $(computation.processes)  ")
        println(io, "initialConfs:             $(computation.initialConfs)  ")
        println(io, "initialLevels:            $(computation.initialLevels)  ")
        println(io, "maxElectronLoss:          $(computation.maxElectronLoss)  ")
        println(io, "NoShakeDisplacements:     $(computation.NoShakeDisplacements)  ")
        println(io, "shakeFromShells:          $(computation.shakeFromShells)  ")
        println(io, "shakeToShells:            $(computation.shakeToShells)  ")
        println(io, "steps:                    $(computation.steps)  ")
    end

    """
    `struct  Cascade.Data`  ... defines a type for an atomic cascade, i.e. lists of radiative, Auger and photoionization lines.

        + name           ::String                               ... A name for the cascade.
        + linesR         ::Array{PhotoEmission.Line,1}          ... List of radiative lines.
        + linesA         ::Array{AutoIonization.Line,1}         ... List of Auger lines.
        + linesP         ::Array{PhotoIonization.Line,1}        ... List of photoionization lines.
    """
    struct  Data
        name             ::String
        linesR           ::Array{PhotoEmission.Line,1}
        linesA           ::Array{AutoIonization.Line,1}
        linesP           ::Array{PhotoIonization.Line,1}
    end

    """
    `Cascade.Data()`  ... (simple) constructor for cascade data.
    """
    function Data()
        Data("", Array{PhotoEmission.Line,1}[], Array{AutoIonization.Line,1}[], Array{PhotoIonization.Line,1}[] )
    end

    # `Base.show(io::IO, data::Cascade.Data)`  ... prepares a proper printout of the variable data::Cascade.Data.
    function Base.show(io::IO, data::Cascade.Data)
        println(io, "name:                    $(data.name)  ")
        println(io, "linesR:                  $(data.linesR)  ")
        println(io, "linesA:                  $(data.linesA)  ")
        println(io, "linesP:                  $(data.linesP)  ")
    end

    #==
    """
    `@enum   Cascade.Property`  ... defines a enumeration for the various properties that can be obtained from the simulation of cascade data.

        + IonDist               ... simulate the 'ion distribution' as it is found after all cascade processes are completed.
        + FinalDist             ... simulate the 'final-level distribution' as it is found after all cascade processes are completed.
        + DecayPathes           ... determine the major 'decay pathes' of the cascade.
        + ElectronIntensity     ... simulate the electron-line intensities as function of electron energy.
        + PhotonIntensity       ... simulate  the photon-line intensities as function of electron energy.
        + ElectronCoincidence   ... simulate electron-coincidence spectra.
    """
    @enum   Property    IonDist  FinalDist  DecayPathes  ElectronIntensity  PhotonIntensity  ElectronCoincidence

    """
    `@enum   Cascade.SimulationMethod`  ... defines a enumeration for the various methods that can be used to run the simulation of cascade data.

        + ProbPropagation     ... to propagate the (occupation) probabilites of the levels until no further changes occur.
        + MonteCarlo          ... to simulate the cascade decay by a Monte-Carlo approach of possible pathes (not yet considered).
        + RateEquations       ... to solve the cascade by a set of rate equations (not yet considered).
    """
    @enum   SimulationMethod    ProbPropagation   MonteCarlo  RateEquations   ==#

    """
    `abstract type  Cascade.Property`
        ... defines an abstract and various singleton types for the different properties that can be obtained from the simulation of cascade data.

        + struct IonDist               ... simulate the 'ion distribution' as it is found after all cascade processes are completed.
        + struct FinalDist             ... simulate the 'final-level distribution' as it is found after all cascade processes are completed.
        + struct DecayPathes           ... determine the major 'decay pathes' of the cascade.
        + struct ElectronIntensity     ... simulate the electron-line intensities as function of electron energy.
        + struct PhotonIntensity       ... simulate  the photon-line intensities as function of electron energy.
        + struct ElectronCoincidence   ... simulate electron-coincidence spectra.
    """
    abstract type  Property                      end
    struct   IonDist              <:  Property   end
    struct   FinalDist            <:  Property   end
    struct   DecayPathes          <:  Property   end
    struct   ElectronIntensity    <:  Property   end
    struct   PhotonIntensity      <:  Property   end
    struct   ElectronCoincidence  <:  Property   end

    """
    abstract type  Cascade.SimulationMethod`
        ... defines a enumeration for the various methods that can be used to run the simulation of cascade data.

        + struct ProbPropagation     ... to propagate the (occupation) probabilites of the levels until no further changes occur.
        + struct MonteCarlo          ... to simulate the cascade decay by a Monte-Carlo approach of possible pathes (not yet considered).
        + struct RateEquations       ... to solve the cascade by a set of rate equations (not yet considered).
    """
    abstract type  SimulationMethod                  end
    struct   ProbPropagation  <:  SimulationMethod   end
    struct   MonteCarlo       <:  SimulationMethod   end
    struct   RateEquations    <:  SimulationMethod   end

    """
    `struct  Cascade.SimulationSettings`  ... defines settings for performing the simulation of some cascade (data).

        + minElectronEnergy ::Float64     ... Minimum electron energy for the simulation of electron spectra.
        + maxElectronEnergy ::Float64     ... Maximum electron energy for the simulation of electron spectra.
        + minPhotonEnergy   ::Float64     ... Minimum photon energy for the simulation of electron spectra.
        + maxPhotonEnergy   ::Float64     ... Maximum photon energy for the simulation of electron spectra..
    """
    struct  SimulationSettings
        minElectronEnergy   ::Float64
        maxElectronEnergy   ::Float64
        minPhotonEnergy     ::Float64
        maxPhotonEnergy     ::Float64
    end

    """
    `Cascade.SimulationSettings()`  ... constructor for an 'empty' instance of a Cascade.Block.
    """
    function SimulationSettings()
        SimulationSettings(0., 1.0e6,  0., 1.0e6 )
    end

    # `Base.show(io::IO, settings::SimulationSettings)`  ... prepares a proper printout of the variable settings::SimulationSettings.
    function Base.show(io::IO, settings::SimulationSettings)
        println(io, "minElectronEnergy:        $(settings.minElectronEnergy)  ")
        println(io, "maxElectronEnergy:        $(settings.maxElectronEnergy)  ")
        println(io, "minPhotonEnergy:          $(settings.minPhotonEnergy)  ")
        println(io, "maxPhotonEnergy:          $(settings.maxPhotonEnergy)  ")
    end

    """
    `struct  Cascade.Simulation`  ... defines a simulation on some given cascade (data).

        + properties      ::Array{Cascade.Property,1}   ... Properties that are considered in this simulation of the cascade (data).
        + method          ::Cascade.SimulationMethod    ... Method that is used in the cascade simulation; cf. Cascade.SimulationMethod.
        + settings        ::Cascade.SimulationSettings  ... Settings for performing these simulations.
    """
    struct  Simulation
        properties        ::Array{Cascade.Property,1}
        method            ::Cascade.SimulationMethod
        settings          ::Cascade.SimulationSettings
    end

    # `Base.show(io::IO, simulation::Cascade.Simulation)`  ... prepares a proper printout of the variable simulation::Cascade.Simulation.
    function Base.show(io::IO, simulation::Cascade.Simulation)
        println(io, "properties:        $(simulation.properties)  ")
        println(io, "method:            $(simulation.method)  ")
        println(io, "settings:          $(simulation.settings)  ")
    end

    """
    `struct  Cascade.LineIndex`  ... defines a line index with regard to the various lineLists of data::Cascade.Data.

        + process      ::Basics.AtomicProcess    ... refers to the particular lineList of cascade (data).
        + index        ::Int64                       ... index of the corresponding line.
    """
    struct  LineIndex
        process        ::Basics.AtomicProcess
        index          ::Int64
    end

    # `Base.show(io::IO, index::Cascade.LineIndex)`  ... prepares a proper printout of the variable index::Cascade.LineIndex.
    function Base.show(io::IO, index::Cascade.LineIndex)
        println(io, "process:        $(index.process)  ")
        println(io, "index:          $(index.index)  ")
    end

    """
    `mutable struct  Cascade.Level`  ... defines a level specification for dealing with cascade transitions.

        + energy       ::Float64                     ... energy of the level.
        + J            ::AngularJ64                  ... total angular momentum of the level
        + parity       ::Basics.Parity               ... total parity of the level
        + NoElectrons  ::Int64                       ... total number of electrons of the ion to which this level belongs.
        + relativeOcc  ::Float64                     ... relative occupation
        + parents      ::Array{Cascade.LineIndex,1}  ... list of parent lines that (may) populate the level.
        + daugthers    ::Array{Cascade.LineIndex,1}  ... list of daugther lines that (may) de-populate the level.
    """
    mutable struct  Level
        energy         ::Float64
        J              ::AngularJ64
        parity         ::Basics.Parity
        NoElectrons    ::Int64
        relativeOcc    ::Float64
        parents        ::Array{Cascade.LineIndex,1}
        daugthers      ::Array{Cascade.LineIndex,1}
    end

    # `Base.show(io::IO, level::Cascade.Level)`  ... prepares a proper printout of the variable level::Cascade.Level.
    function Base.show(io::IO, level::Cascade.Level)
        println(io, "energy:        $(level.energy)  ")
        println(io, "J:             $(level.J)  ")
        println(io, "parity:        $(level.parity)  ")
        println(io, "NoElectrons:   $(level.NoElectrons)  ")
        println(io, "relativeOcc:   $(level.relativeOcc)  ")
        println(io, "parents:       $(level.parents)  ")
        println(io, "daugthers:     $(level.daugthers)  ")
    end

    """
    `Cascade.computeDecayYieldOutcome(outcome::DecayYield.Outcome, linesR::Array{PhotoEmission.Line,1},
                                      linesA::Array{AutoIonization.Line,1}, settings::DecayYield.Settings)`
        ... to compute the flourescence and Auger yields for a single outcome as specified by its level; an
            outcome::DecayYield.Outcome is returned in which are physical parameters are now specified.
    """
    function computeDecayYieldOutcome(outcome::DecayYield.Outcome, linesR::Array{PhotoEmission.Line,1},
                                      linesA::Array{AutoIonization.Line,1}, settings::DecayYield.Settings)
        # Identify the level key of the given level also in the lists of radiative and Auger lines
        level = outcome.level;    levelKey = LevelKey( LevelSymmetry(level.J, level.parity), level.index, level.energy, 0.)
        similarKey = LevelKey();  rateR = 0.;    rateA = 0.;   NoPhotonLines = 0;   NoAugerLines = 0
        ##x println("** levelKey = $levelKey")
        for  line in linesR
            compareKey = LevelKey( LevelSymmetry(line.initialLevel.J, line.initialLevel.parity), line.initialLevel.index, line.initialLevel.energy, 0.)
            if   Basics.isSimilar(levelKey, compareKey, 1.0e-3)    println("** compareKey = $compareKey");   similarKey = deepcopy(compareKey)    end
        end
        if   similarKey == LevelKey()    error("No similar level found !")   end

        for  line in linesR
            ##x println("** Radiative: similarKey =$similarKey     $(line.initialLevel)")
            if  similarKey == LevelKey( LevelSymmetry(line.initialLevel.J, line.initialLevel.parity), line.initialLevel.index, line.initialLevel.energy, 0.)
                rateR = rateR + line.photonRate;   NoPhotonLines = NoPhotonLines + 1
                ##x println("NoPhotonLines = $NoPhotonLines")
            end
        end
        for  line in linesA
            ##x println("** Auger: similarKey =$similarKey     $(line.initialLevel)")
            if  similarKey == LevelKey( LevelSymmetry(line.initialLevel.J, line.initialLevel.parity), line.initialLevel.index, line.initialLevel.energy, 0.)
                rateA = rateA + line.totalRate;   NoAugerLines = NoAugerLines + 1
            end
        end

        omegaR = rateR / (rateR + rateA);   omegaA = rateA / (rateR + rateA)
        newOutcome = DecayYield.Outcome(level, NoPhotonLines, NoAugerLines, rateR, rateA, omegaR, omegaA)
        return( newOutcome )
    end

    """
    `Cascade.computeSteps(comp::Cascade.Computation, stepList::Array{Cascade.Step,1})`
        ... computes in turn all the requested transition amplitudes and PhotoEmission.Line's, AutoIonization.Line's,
            etc. for all pre-specified decay steps of the cascade. When compared with the standard atomic process computations,
            however, the amount of output is largely reduced. A set of  data::Cascade.Data  is returned.
    """
    function computeSteps(comp::Cascade.Computation, stepList::Array{Cascade.Step,1})
        linesA = AutoIonization.Line[];    linesR = PhotoEmission.Line[];    linesP = PhotoIonization.Line[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        nt = 0;   st = 0
        for  step  in  stepList
            st = st + 1
            nc = length(step.initialMultiplet.levels) * length(step.finalMultiplet.levels)
            println("\n  $st) Perform $(string(step.process)) amplitude computations for up to $nc lines (without selection rules): ")
            if  printSummary   println(iostream, "\n* $st) Perform $(string(step.process)) amplitude computations for " *
                                                 "up to $nc lines (without selection rules): ")   end
            if      step.process == Basics.Auger
                newLines = AutoIonization.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.nuclearModel, comp.grid,
                                                              step.settings, output=true, printout=false)
                append!(linesA, newLines);    nt = length(linesA)
            elseif  step.process == Basics.Radiative
                newLines = PhotoEmission.computeLinesCascade(step.finalMultiplet, step.initialMultiplet, comp.grid,
                                                             step.settings, output=true, printout=false)
                append!(linesR, newLines);    nt = length(linesR)
            else   error("Unsupported atomic process for cascade computations.")
            end
            println("     Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, giving now rise " *
                    "to a total of $nt $(string(step.process)) lines." )
            if  printSummary   println(iostream, "\n*    Step $st:: A total of $(length(newLines)) $(string(step.process)) lines are calculated, " *
                                                 "giving now rise to a total of $nt $(string(step.process)) lines." )   end
        end
        #
        data = Cascade.Data(comp.name, linesR, linesA, linesP)
    end

    """
    `Cascade.determineSteps(comp::Cascade.Computation, blockList::Array{Cascade.Block,1})`
        ... determines all steps::Cascade.Step that need to be computed for this cascade. It cycles through the given processes and
            distinguished between the different cascade approaches. It also checks that the averaged energies of the configuration
            allows such a step energetically. A stepList::Array{Cascade.Step,1} is returned
    """
    function determineSteps(comp::Cascade.Computation, blockList::Array{Cascade.Block,1})
        stepList = Cascade.Step[]
        ##x println("comp.approach = $(comp.approach)")
        if  comp.approach  in  [Cascade.AverageSCA(), Cascade.SCA()]
            for  a = 1:length(blockList)
                for  b = 1:length(blockList)
                    minEn = 100000.;   maxEn = -100000.
                    for  p = 1:length(blockList[a].multiplet.levels),  q = 1:length(blockList[b].multiplet.levels)
                        minEn = min(minEn, blockList[a].multiplet.levels[p].energy - blockList[b].multiplet.levels[q].energy)
                        maxEn = max(maxEn, blockList[a].multiplet.levels[p].energy - blockList[b].multiplet.levels[q].energy)
                    end
                    for  process  in  comp.processes
                        if      process == Basics.Radiative
                            if  a == b   ||   minEn < 0.    continue   end
                            if  blockList[a].NoElectrons == blockList[b].NoElectrons
                                settings = PhotoEmission.Settings([E1], [UseBabushkin], false, false, false, Tuple{Int64,Int64}[], 0., 0., 1.0e6)
                                push!( stepList, Cascade.Step(process, settings, blockList[a].confs, blockList[b].confs,
                                                              blockList[a].multiplet, blockList[b].multiplet) )
                            end
                        elseif  process == Basics.Auger
                            if  a == b   ||   minEn < 0.    continue   end
                            if  blockList[a].NoElectrons == blockList[b].NoElectrons + 1
                                settings = AutoIonization.Settings(false, false, false, Tuple{Int64,Int64}[], 0., 1.0e6, 2, "Coulomb")
                                push!( stepList, Cascade.Step(process, settings, blockList[a].confs, blockList[b].confs,
                                                              blockList[a].multiplet, blockList[b].multiplet) )
                            end
                        else    error("stop a")
                        end
                    end
                end
            end
            #
        else  error("Unsupported cascade approach.")
        end
        return( stepList )
    end

    """
    `Cascade.displayBlocks(stream::IO, blockList::Array{Cascade.Block,1})`
        ... group & display the blocks of the cascade with same No. of electrons; this blocks are displayed with the
            minimum and maximum energy of each multiplet.
    """
    function displayBlocks(stream::IO, blockList::Array{Cascade.Block,1})
        #
        println(stream, "\n* Configuration 'blocks' (multiplets) in the given cascade model: \n")
        println(stream, "  ", TableStrings.hLine(134))
        println(stream, "      No.   Configurations                                                                       " *
                        "      Range of total energies " * TableStrings.inUnits("energy") )
        println(stream, "  ", TableStrings.hLine(134))
        i = 0
        for  block  in  blockList
            i = i + 1;
            sa = "   " * TableStrings.flushright( 6, string(i); na=2)
            sb = " ";         for conf  in blockList[i].confs   sb = sb * string(conf) * ", "    end
            en = Float64[];   for level in  block.multiplet.levels    push!(en, level.energy)    end
            minEn = minimum(en);   minEn = Defaults.convertUnits("energy: from atomic", minEn)
            maxEn = maximum(en);   maxEn = Defaults.convertUnits("energy: from atomic", maxEn)
            sa = sa * TableStrings.flushleft(90, sb[1:end-2]; na=2)
            sa = sa * TableStrings.flushleft(30, string( round(minEn)) * " ... " * string( round(maxEn)); na=2)
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(134))

        return( nothing )
    end

    """
    `Cascade.displayInitialLevels(stream::IO, multiplet::Multiplet, initialLevels::Array{Tuple{Int64,Float64},1})`
        ... display the calculated initial levels to screen together with their given relative occupation.
    """
    function displayInitialLevels(stream::IO, multiplet::Multiplet, initialLevels::Array{Tuple{Int64,Float64},1})
        println(stream, " ")
        println(stream, "* Initial levels of the given cascade, relative to the lowest, and their given occupation:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(64))
        println(stream, "    Level  J Parity          Energy " * TableStrings.inUnits("energy") * "         rel. occupation ")
        println(stream, "  ", TableStrings.hLine(64))
        for  i = 1:length(multiplet.levels)
            lev = multiplet.levels[i]
            en  = lev.energy - multiplet.levels[1].energy;    en_requested = Defaults.convertUnits("energy: from atomic", en)
            wx = 0.
            for  ilevel in initialLevels
                if  i in ilevel   wx = ilevel[2];    break    end
            end
            sb  = "          "   * string(wx)
            sc  = "   "  * TableStrings.level(i) * "     " * string(LevelSymmetry(lev.J, lev.parity)) * "     "
            @printf(stream, "%s %.15e %s %s", sc, en_requested, sb, "\n")
        end
        println(stream, "  ", TableStrings.hLine(64))
        return( nothing )
    end

    """
    `Cascade.displaySteps(stream::IO, steps::Array{Cascade.Step,1})`
        ... displays all predefined steps in a neat table and supports to delete individual steps from the list.
    """
    function displaySteps(stream::IO, steps::Array{Cascade.Step,1})
        println(stream, " ")
        println(stream, "* Steps that are defined for the current cascade due to the given approach:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(170))
        sa = "  "
        sa = sa * TableStrings.center( 9, "Step-No"; na=2)
        sa = sa * TableStrings.flushleft(11, "Process"; na=1)
        sa = sa * TableStrings.flushleft(55, "Initial:  No CSF, configuration(s)"; na=4)
        sa = sa * TableStrings.flushleft(55, "Final:  No CSF, configuration(s)"; na=4)
        sa = sa * TableStrings.flushleft(40, "Energies from ... to in " * TableStrings.inUnits("energy"); na=4)
        println(stream, sa)
        println(stream, "  ", TableStrings.hLine(170))
        #
        for  i = 1:length(steps)
            sa = " " * TableStrings.flushright( 7, string(i); na=5)
            sa = sa  * TableStrings.flushleft( 11, string(steps[i].process); na=1)
            sb = "";   for conf in steps[i].initialConfs   sb = sb * string(conf) * ", "    end
            sa = sa  * TableStrings.flushright( 5, string( length(steps[i].initialMultiplet.levels[1].basis.csfs) )*", "; na=0)
            sa = sa  * TableStrings.flushleft( 50, sb[1:end-2]; na=4)
            sb = "";   for conf in steps[i].finalConfs   sb = sb * string(conf) * ", "    end
            sa = sa  * TableStrings.flushright( 5, string( length(steps[i].finalMultiplet.levels[1].basis.csfs) )*", "; na=0)
            sa = sa  * TableStrings.flushleft( 50, sb[1:end-2]; na=4)
            minEn = 1000.;   maxEn = -1000.;
            for  p = 1:length(steps[i].initialMultiplet.levels),  q = 1:length(steps[i].finalMultiplet.levels)
                minEn = min(minEn, steps[i].initialMultiplet.levels[p].energy - steps[i].finalMultiplet.levels[q].energy)
                maxEn = max(maxEn, steps[i].initialMultiplet.levels[p].energy - steps[i].finalMultiplet.levels[q].energy)
            end
            minEn = Defaults.convertUnits("energy: from atomic", minEn);   maxEn = Defaults.convertUnits("energy: from atomic", maxEn)
            sa = sa * string( round(minEn)) * " ... " * string( round(maxEn))
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(170))
    end

    """
    `Cascade.generateBlocks(comp::Cascade.Computation, confs::Array{Configuration,1})`
        ... generate all block::Cascade.Block that need to be computed for this cascade and compute the corresponding multiplets.
            The different cascade approches follow different strategies in defining and computing these blocks.
            A blockList::Array{Cascade.Block,1} is returned.
    """
    function generateBlocks(comp::Cascade.Computation, confs::Array{Configuration,1}, initalOrbitals::Dict{Subshell, Orbital})
        blockList = Cascade.Block[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        if    comp.approach == AverageSCA()
            println("\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println("    + orbitals from the initial multiplet are applied throughout; ")
            println("    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; ")
            println("    + only E1 dipole transitions are applied in all radiative decay stets; ")
            println("    + for each decay step, a (single) set of continuum orbitals with an configuration averaged energy is applied " *
                           "for all transitions of the same step. \n")
            if  printSummary
            println(iostream, "\n* In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println(iostream, "    + orbitals from the initial multiplet are applied throughout; ")
            println(iostream, "    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; ")
            println(iostream, "    + only E1 dipole transitions are applied in all radiative decay stets; ")
            println(iostream, "    + for each decay step, a (single) set of continuum orbitals with an configuration averaged energy is applied " *
                              "for all transitions of the same step. \n")
            end
            #
            for  confa  in confs
                print("  Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")
                if  printSummary   println(iostream, "\n*  Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")   end
                multiplet = perform("computation: mutiplet from orbitals, no CI, CSF diagonal", [confa],  initalOrbitals,
                                    comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
                ##x if  printSummary   println(iostream, "* ... and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")   end
            end
        elseif    comp.approach == SCA()
            println("\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println("    + orbitals are generated independently for each multiplet (block); ")
            println("    + configuration interaction is included for each block; ")
            println("    + only E1 dipole transitions are applied in all radiative decay stets; \n")
            if  printSummary
            println(iostream, "\n* In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println(iostream, "    + orbitals are generated independently for each multiplet (block); ")
            println(iostream, "    + configuration interaction is included for each block; ")
            println(iostream, "    + only E1 dipole transitions are applied in all radiative decay stets; \n")
            end
            #
            i = 0
            for  confa  in confs
                ## i = i + 1;    if   i in [1,2, 4,5,6,7,8,9,10,11,12,13,14]  ||  i > 15   println("  Block $i omitted.");    continue    end
                ## i = i + 1;    if   i < 11  ||  i > 11   println("  Block $i omitted.");    continue    end
                print("  Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")
                if  printSummary   print(iostream, "* Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")   end
                basis     = perform("computation: SCF", [confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                multiplet = perform("computation: CI",  basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
                ##x if  printSummary   println(iostream, "and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")   end
            end
        else  error("Unsupported cascade approach.")
        end

        return( blockList )
    end

    """
    `Cascade.generateConfigurationList(initialConfs::Array{Configuration,1}, further::Int64, NoShake::Int64)`
        ... generates all possible (decay) configurations with up to further holes and with NoShake displacements. First, all
            configuratons are generated for which the hole is either moved 'outwards' or is moved and a second 'outer' hole is
            created; this step is repated further + 2 times. From the generated list, only those configurations are kept with
            up to further holes, when compared with the initial configuration. A confList::Array{Configuration,1} is returned.
    """
    function generateConfigurationList(initialConfs::Array{Configuration,1}, further::Int64, NoShake::Int64)
        confList = copy(initialConfs);    cList = copy(initialConfs);   initialNoElectrons = initialConfs[1].NoElectrons
        # First, move and generate new 'outer' hole without displacements
        for  fur = 1:further+1
            newConfList = Configuration[]
            for conf  in cList
                holeList = Basics.determineHoleShells(conf)
                for  holeShell in holeList
                    wa = generateConfigurationsWith1OuterHole(conf,  holeShell);   append!(newConfList, wa)
                    wa = generateConfigurationsWith2OuterHoles(conf, holeShell);   append!(newConfList, wa)
                end
            end
            if  length(newConfList) > 0    newConfList = Basics.excludeDoubles(newConfList)    end
            cList = newConfList
            append!(confList, newConfList)
        end
        # Make sure that only configurations with up to further holes are returned
        newConfList = Configuration[]
        for   conf in confList
            if  conf.NoElectrons + further >= initialNoElectrons   push!(newConfList, conf)    end
        end
        # Add further shake-displacements if appropriate
        newConfList = Basics.excludeDoubles(newConfList)
        return( newConfList )
    end

    """
    `Cascade.generateConfigurationsWith1OuterHole(conf,  holeShell)`
        ... generates all possible (decay) configurations where the hole in holeShell is moved 'outwards'.
            A confList::Array{Configuration,1} is returned.
    """
    function generateConfigurationsWith1OuterHole(conf::Configuration,  holeShell::Shell)
         shList = Basics.generate("shells: ordered list for NR configurations", [conf]);   i0 = 0
         for  i = 1:length(shList)
             if   holeShell == shList[i]    i0 = i;    break    end
         end
         if  i0 == 0   error("stop a")   end
         #
         # Now move the hole 'outwards'
         confList = Configuration[]
         for  i = i0+1:length(shList)
             if  haskey(conf.shells, shList[i])  &&  conf.shells[ shList[i] ] >= 1
                 newshells = copy( conf.shells )
                 newshells[ shList[i] ] = newshells[ shList[i] ] - 1
                 newshells[ holeShell ] = newshells[ holeShell ] + 1
                 push!(confList, Configuration( newshells, conf.NoElectrons ) )
             end
         end
         return( confList )
    end

    """
    `Cascade.generateConfigurationsWith2OuterHoles(conf,  holeShell)`
        ... generates all possible (decay) configurations where the hole in holeShell is moved 'outwards'.
            A confList::Array{Configuration,1} is returned.
    """
    function generateConfigurationsWith2OuterHoles(conf::Configuration,  holeShell::Shell)
         shList = Basics.generate("shells: ordered list for NR configurations", [conf]);   i0 = 0
         for  i = 1:length(shList)
             if   holeShell == shList[i]    i0 = i;    break    end
         end
         if  i0 == 0   error("stop a")   end
         #
         # Now move the hole 'outwards'
         confList = Configuration[]
         for  i = i0+1:length(shList)
             if  haskey(conf.shells, shList[i])  &&  conf.shells[ shList[i] ] >= 2
                 newshells = copy( conf.shells )
                 newshells[ shList[i] ] = newshells[ shList[i] ] - 2
                 newshells[ holeShell ] = newshells[ holeShell ] + 1
                 push!(confList, Configuration( newshells, conf.NoElectrons - 1 ) )
             end
             #
             for  j = i0+1:length(shList)
                 if  i != j   &&   haskey(conf.shells, shList[i])  &&  conf.shells[ shList[i] ] >= 1   &&
                                   haskey(conf.shells, shList[j])  &&  conf.shells[ shList[j] ] >= 1
                     newshells = copy( conf.shells )
                     newshells[ shList[i] ] = newshells[ shList[i] ] - 1
                     newshells[ shList[j] ] = newshells[ shList[j] ] - 1
                     newshells[ holeShell ] = newshells[ holeShell ] + 1
                     push!(confList, Configuration( newshells, conf.NoElectrons - 1 ) )
                 end
             end
         end
         return( confList )
    end

    """
    `Cascade.groupDisplayConfigurationList(Z::Float64, confs::Array{Configuration,1})`
        ... group & display the configuration list into sublists with the same No. of electrons; this lists are displayed together
            with an estimated total energy.
    """
    function groupDisplayConfigurationList(Z::Float64, confs::Array{Configuration,1})
        minNoElectrons = 1000;   maxNoElectrons = 0
        for  conf in confs
            minNoElectrons = min(minNoElectrons, conf.NoElectrons)
            maxNoElectrons = max(maxNoElectrons, conf.NoElectrons)
        end
        #
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        println("\n  Electron configuration used in the cascade:")
        if  printSummary   println(iostream, "\n* Electron configuration used in the cascade:")    end
        confList = Configuration[];   nc = 0
        for  n = maxNoElectrons:-1:minNoElectrons
            println("\n    Configuration(s) with $n electrons:")
            if  printSummary   println(iostream, "\n    Configuration(s) with $n electrons:")      end
            for  conf in confs
                if n == conf.NoElectrons
                    nc = nc + 1
                    push!(confList, conf )
                    wa = Basics.provide("binding energy", round(Int64, Z), conf);    wa = Defaults.convertUnits("energy: from atomic", wa)
                    sa = "   av. BE = "  * string( round(-wa) ) * "  " * TableStrings.inUnits("energy")
                    println("      " * string(conf) * sa * "      ($nc)" )
                    if  printSummary   println(iostream, "      " * string(conf) * sa * "      ($nc)")      end
                end
            end
        end

        println("\n  A total of $nc configuration have been defined for this cascade, and selected configurations could be " *
                "removed here:  [currently not supported] \n")
        if  printSummary   println(iostream, "\n* A total of $nc configuration have been defined for this cascade, and selected " *
                                             "configurations could be removed here:  [currently not supported] \n")      end
        return( confList )
    end

    """
    `Cascade.modifySteps(stepList::Array{Cascade.Step,1})`
        ... allows the user to modify the steps, for instance, by deleting selected steps of the cascade or by modifying the settings of
            one or several steps. A newStepList::Array{Cascade.Step,1} for which the transition data are eventually computed.
    """
    function modifySteps(stepList::Array{Cascade.Step,1})
        #
        newStepList = Cascade.Step[]
        #
        println("\n  Here, modify the individual steps explicitly in the code, if needed, ...... and just do it !!")
        #
        #  Delete individual steps from stepList
        #  if  i in [1,2,5, ...] modify the particular settings, etc.
        for  i = 1:length(stepList)
            step = stepList[i]
            #
            if  i in []
                println("  Modify step $i :")
                newStep = Cascade.Step(step.process, step.settings, step.initialConfs, step.finalConfs, step.initialMultiplet, step.initialMultiplet)
                push!(newStepList, newStep)
            else
                push!(newStepList, step)
            end
        end
        #
        # wa = [1,2,3]
        # delete from list
        #
        println("\n  A total of $(length(newStepList)) steps are still defined in the cascade.")
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   println(iostream, "\n* A total of $(length(newStepList)) steps are still defined in the cascade.")    end

        return( newStepList )
    end

    """
    `Cascade.simulateLevelDistribution(simulation::Cascade.Simulation, data::Cascade.Data)`
        ... sorts all levels as given by data and propagates their (occupation) probability until no further changes occur. For this
            propagation, it runs through all levels and propagates the probabilty until no level probability changes anymore. The final
            level distribution is then used to derive the ion distribution or the level distribution, if appropriate. Nothing is returned.
    """
    function simulateLevelDistribution(simulation::Cascade.Simulation, data::Cascade.Data)
        levels = Cascade.extractLevels(data)
        Cascade.displayLevelTree(data.name, levels, data)
        Cascade.propagateProbability!(levels, data)
        println("aaa")
        if   Cascade.IonDist()   in simulation.properties    Cascade.displayIonDistribution(data.name, levels)     end
        if   Cascade.FinalDist() in simulation.properties    Cascade.displayLevelDistribution(data.name, levels)   end

        return( nothing )
    end

    """
    `Cascade.displayIonDistribution(sc::String, levels::Array{Cascade.Level,1})`
        ... displays the (current or final) ion distribution in a neat table. Nothing is returned.
    """
    function displayIonDistribution(sc::String, levels::Array{Cascade.Level,1})
        minElectrons = 1000;   maxElectrons = 0
        for  level in levels   minElectrons = min(minElectrons, level.NoElectrons);   maxElectrons = max(maxElectrons, level.NoElectrons)   end
        println(" ")
        println("  (Final) Ion distribution for the cascade:  $sc   mine = $minElectrons maxe = $maxElectrons ")
        println(" ")
        println("  ", TableStrings.hLine(31))
        sa = "  "
        sa = sa * TableStrings.center(14, "No. electrons"; na=4)
        sa = sa * TableStrings.center(10,"Rel. occ.";      na=2)
        println(sa)
        println("  ", TableStrings.hLine(31))
        for n = maxElectrons:-1:minElectrons
            sa = "             " * string(n);   sa = sa[end-10:end];   prob = 0.
            for  level in levels    if  n == level.NoElectrons   prob = prob + level.relativeOcc    end    end
            sa = sa * "         " * @sprintf("%.5e", prob)
            println(sa)
        end
        println("  ", TableStrings.hLine(31))

        return( nothing )
    end

    """
    `Cascade.displayLevelDistribution(sc::String, levels::Array{Cascade.Level,1})`
        ... displays the (current or final) level distribution in a neat table. Only those levels with a non-zero
            occupation are displayed here. Nothing is returned.
    """
    function displayLevelDistribution(sc::String, levels::Array{Cascade.Level,1})
        minElectrons = 1000;   maxElectrons = 0;   energies = zeros(length(levels))
        for  i = 1:length(levels)
            minElectrons = min(minElectrons, levels[i].NoElectrons);   maxElectrons = max(maxElectrons, levels[i].NoElectrons)
            energies[i]  = levels[i].energy
        end
        enIndices = sortperm(energies, rev=true)
        # Now printout the results
        println(" ")
        println("  (Final) Level distribution for the cascade:  $sc")
        println(" ")
        println("  ", TableStrings.hLine(57))
        sa = "  "
        sa = sa * TableStrings.center(14, "No. electrons"; na=2)
        sa = sa * TableStrings.center( 6, "J^P"          ; na=2);
        sa = sa * TableStrings.center(16, "Energy " * TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(10, "Rel. occ.";                                    na=2)
        println(sa)
        println("  ", TableStrings.hLine(57))
        for n = maxElectrons:-1:minElectrons
            sa = "            " * string(n);   sa = sa[end-10:end]
            for  en in enIndices
                if  n == levels[en].NoElectrons  ##    &&  levels[en].relativeOcc > 0
                    sb = sa * "       " * string( LevelSymmetry(levels[en].J, levels[en].parity) )     * "    "
                    sb = sb * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", levels[en].energy))  * "      "
                    sb = sb * @sprintf("%.5e", levels[en].relativeOcc)
                    sa = "           "
                    println(sb)
                end
            end
        end
        println("  ", TableStrings.hLine(57))

        return( nothing )
    end

    """
    `Cascade.displayLevelTree(sc::String, levels::Array{Cascade.Level,1}, data::Cascade.Data)`
        ... displays all defined levels  in a neat table, together with their No. of electrons, symmetry, level energy,
            current (relative) population as well as analogue information about their parents and daugther levels. This
            enables one to recognize (and perhaps later add) missing parent and daughter levels. Nothing is returned.
    """
    function displayLevelTree(sc::String, levels::Array{Cascade.Level,1}, data::Cascade.Data)
        minElectrons = 1000;   maxElectrons = 0;   energies = zeros(length(levels))
        for  i = 1:length(levels)
            minElectrons = min(minElectrons, levels[i].NoElectrons);   maxElectrons = max(maxElectrons, levels[i].NoElectrons)
            energies[i]  = levels[i].energy
        end
        enIndices = sortperm(energies, rev=true)
        # Now printout the results
        println(" ")
        println("  Level tree of this cascade:  $sc")
        println(" ")
        println("  ", TableStrings.hLine(175))
        sa = "  "
        sa = sa * TableStrings.center( 6, "No. e-"; na=2)
        sa = sa * TableStrings.center( 6, "J^P"          ; na=2);
        sa = sa * TableStrings.center(16, "Energy " * TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(10, "Rel. occ.";                                    na=5)
        sb = "Parents P(A: No_e, sym, energy) and Daughters D(R: No_e, sym, energy);  all energies in " * TableStrings.inUnits("energy")
        sa = sa * TableStrings.flushleft(100, sb; na=2)
        println(sa)
        println("  ", TableStrings.hLine(175))
        for n = maxElectrons:-1:minElectrons
            sa = "            " * string(n);   sa = sa[end-5:end]
            for  en in enIndices
                if  n == levels[en].NoElectrons
                    sb = sa * "      " * string( LevelSymmetry(levels[en].J, levels[en].parity) )      * "   "
                    sb = sb * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", levels[en].energy))  * "      "
                    sb = sb * @sprintf("%.5e", levels[en].relativeOcc)                                 * "    "
                    pProcessSymmetryEnergyList = Tuple{AtomicProcess,Int64,LevelSymmetry,Float64}[]
                    dProcessSymmetryEnergyList = Tuple{Basics.AtomicProcess,Int64,LevelSymmetry,Float64}[]
                    for  p in levels[en].parents
                        idx = p.index
                        if      p.process == Basics.Auger         lev = data.linesA[idx].initialLevel
                        elseif  p.process == Basics.Radiative    lev = data.linesR[idx].initialLevel
                        elseif  p.process == Basics.Photo         lev = data.linesP[idx].initialLevel
                        else    error("stop a")    end
                        push!( pProcessSymmetryEnergyList, (p.process, lev.basis.NoElectrons, LevelSymmetry(lev.J, lev.parity), lev.energy) )
                    end
                    for  d in levels[en].daugthers
                        idx = d.index
                        if      d.process == Basics.Auger         lev = data.linesA[idx].finalLevel
                        elseif  d.process == Basics.Radiative    lev = data.linesR[idx].finalLevel
                        elseif  d.process == Basics.Photo         lev = data.linesP[idx].finalLevel
                        else    error("stop b")    end
                        push!( dProcessSymmetryEnergyList, (d.process, lev.basis.NoElectrons, LevelSymmetry(lev.J, lev.parity), lev.energy) )
                    end
                    wa = TableStrings.processSymmetryEnergyTupels(120, pProcessSymmetryEnergyList, "P")
                    if  length(wa) > 0    sc = sb * wa[1];    println( sc )    else    println( sb )   end
                    for  i = 2:length(wa)
                        sc = TableStrings.hBlank( length(sb) ) * wa[i];    println( sc )
                    end
                    wa = TableStrings.processSymmetryEnergyTupels(120, dProcessSymmetryEnergyList, "D")
                    for  i = 1:length(wa)
                        sc = TableStrings.hBlank( length(sb) ) * wa[i];    println( sc )
                    end
                    sa = "      "
                end
            end
        end
        println("  ", TableStrings.hLine(175))

        return( nothing )
    end

    """
    `Cascade.propagateProbability!(levels::Array{Cascade.Level,1}, data::Cascade.Data)`
        ... propagates the relative level occupation through the levels of the cascade until no further change occur in the
            relative level occupation. The argument levels is modified during the propagation, but nothing is returned otherwise.
    """
    function propagateProbability!(levels::Array{Cascade.Level,1}, data::Cascade.Data)
        n = 0
        println("\n  Probability propagation through $(length(levels)) levels of the cascade:")
        while true
            n = n + 1;    totalProbability = 0.
            print("    $n-th round ... ")
            for  level in levels
                if   level.relativeOcc > 0.   && length(level.daugthers) > 0
                    # A level with relative occupation > 0 has still 'daugther' levels; collect all decay rates for this level
                    prob  = level.relativeOcc;   totalProbability = totalProbability + prob;   rates = zeros(length(level.daugthers))
                    level.relativeOcc = 0.
                    for  i = 1:length(level.daugthers)
                        idx = level.daugthers[i].index
                        if      level.daugthers[i].process == Basics.Radiative    rate[i] = data.lineR[idx].photonRate.Babushkin
                        elseif  level.daugthers[i].process == Basics.Auger         rate[i] = data.lineA[idx].totalRate
                        else    error("stop a; process = $(level.daugthers[i].process) ")
                        end
                    end
                    totalRate = sum(rates)
                    # Shift the relative occupation to the 'daugther' levels due to the different decay pathes
                    for  i = 1:length(level.daugthers)
                        idx = level.daugthers[i].index
                        if      level.daugthers[i].process == Basics.Radiative    line = data.lineR[idx]
                        elseif  level.daugthers[i].process == Basics.Auger         line = data.lineA[idx]
                        else    error("stop b; process = $(level.daugthers[i].process) ")
                        end
                        level = Cascade.Level( line.finalLevel.energy, line.finalLevel.J, line.finalLevel.parity,
                                               line.finalLevel.basis.NoElectrons, 0., Cascade.LineIndex[], Cascade.LineIndex[] )
                        kk    = Cascade.findLevelIndex(level, levels)
                        levels[kk].relativeOcc = levels[kk].relativeOcc + prob * rates[i] / totalRate
                    end
                end
            end
            println("has propagated a total of $totalProbability level occupation.")
            # Cycle once more if the relative occupation has still changed
            if  totalProbability == 0.    break    end
        end

        return( nothing )
    end

    """
    `Cascade.findLevelIndex(level::Cascade.Level, levels::Array{Cascade.Level,1})`
        ... find the index of the given level within the given list of levels; an idx::Int64 is returned and an error message is
            issued if the level is not found in the list.
    """
    function findLevelIndex(level::Cascade.Level, levels::Array{Cascade.Level,1})
        for  k = 1:length(levels)
            if  level.energy == levels[k].energy  &&   level.J == levels[k].J   &&   level.parity == levels[k].parity   &&
                level.NoElectrons == levels[k].NoElectrons
                kk = k;   return( kk )
            end
        end
        error("findLevelIndex():  No index was found;\n   level = $(level) ")
    end

    """
    `Cascade.extractLevels(data::Cascade.Data)`
        ... extracts and sorts all levels from the given cascade data into a new levelList::Array{Cascade.Level,1} to simplify the
            propagation of the probabilities. In this list, every level of the overall cascade just occurs just once, together
            with its parent lines (which may populate the level) and the daugther lines (to which the pobability may decay).
            A levelList::Array{Cascade.Level,1} is returned.
    """
    function extractLevels(data::Cascade.Data)
        levels = Cascade.Level[]
        print("\n  Extract, sort and unify the list of levels of the cascade ... ")
        for  i = 1:length(data.linesR)
            line = data.linesR[i]
            iLevel = Cascade.Level( line.initialLevel.energy, line.initialLevel.J, line.initialLevel.parity, line.initialLevel.basis.NoElectrons,
                                    line.initialLevel.relativeOcc, Cascade.LineIndex[], [ Cascade.LineIndex(Basics.Radiative, i)] )
            Cascade.pushLevels!(levels, iLevel)
            fLevel = Cascade.Level( line.finalLevel.energy, line.finalLevel.J, line.finalLevel.parity, line.finalLevel.basis.NoElectrons,
                                    line.finalLevel.relativeOcc, [ Cascade.LineIndex(Basics.Radiative, i)], Cascade.LineIndex[] )
            Cascade.pushLevels!(levels, fLevel)
        end

        for  i = 1:length(data.linesA)
            line = data.linesA[i]
            iLevel = Cascade.Level( line.initialLevel.energy, line.initialLevel.J, line.initialLevel.parity, line.initialLevel.basis.NoElectrons,
                                    line.initialLevel.relativeOcc, Cascade.LineIndex[], [ Cascade.LineIndex(Basics.Auger, i)] )
            Cascade.pushLevels!(levels, iLevel)
            fLevel = Cascade.Level( line.finalLevel.energy, line.finalLevel.J, line.finalLevel.parity, line.finalLevel.basis.NoElectrons,
                                    line.finalLevel.relativeOcc, [ Cascade.LineIndex(Basics.Auger, i)], Cascade.LineIndex[] )
            Cascade.pushLevels!(levels, fLevel)
        end

        for  i = 1:length(data.linesP)
            line = data.linesP[i]
            iLevel = Cascade.Level( line.initialLevel.energy, line.initialLevel.J, line.initialLevel.parity, line.initialLevel.basis.NoElectrons,
                                    line.initialLevel.relativeOcc, Cascade.LineIndex[], [ Cascade.LineIndex(Basics.Photo, i)] )
            Cascade.pushLevels!(levels, iLevel)
            fLevel = Cascade.Level( line.finalLevel.energy, line.finalLevel.J, line.finalLevel.parity, line.finalLevel.basis.NoElectrons,
                                    line.finalLevel.relativeOcc, [ Cascade.LineIndex(Basics.Photo, i)], Cascade.LineIndex[] )
            Cascade.pushLevels!(levels, fLevel)
        end

        println("a total of $(length(levels)) levels were found.")
        return( levels )
    end

    """
    `Cascade.pushLevels!(levels::Array{Cascade.Level,1}, newLevel::Cascade.Level)`
        ... push's the information of newLevel of levels. This is the standard 'push!(levels, newLevel)' if newLevel is not yet
            including in levels, and the proper modification of the parent and daugther lines of this level otherwise. The argument
            levels::Array{Cascade.Level,1} is modified and nothing is returned otherwise.
    """
    function pushLevels!(levels::Array{Cascade.Level,1}, newLevel::Cascade.Level)
        for  i = 1:length(levels)
            if  newLevel.energy == levels[i].energy  &&  newLevel.J == levels[i].J  &&  newLevel.parity == levels[i].parity
                append!(levels[i].parents,   newLevel.parents)
                append!(levels[i].daugthers, newLevel.daugthers)
                return( nothing )
            end
        end
        push!( levels, newLevel)
        ##x info("... one level added, n = $(length(levels)) ")
        return( nothing )
    end

end # module


text:
module NIDAQ

include("DAQmx/DAQmx.jl")

using .DAQmx
using PrettyTables, OrderedCollections, InteractiveUtils, OffsetArrays
using DataStructures, Distributed
using DelimitedFiles, Dates, HDF5
import Base.write, Base.read!, Base.push!, Base.show

export DAQmx
export DAQDevice, DefaultDev, DAQTask, TaskChannel
export AnalogIn, AnalogOut, DigitalIn, DigitalOut, CounterIn, CounterOut
export AI, AO, DI, DO, CI, CO

export lsdev, lschan
export start, stop, clear, isrunning
export taskalloc, read!
export rate, maxrate
export record!
export pulses, pulses!, pulsetiming, pulsetiming!

include("types.jl")
include("constructors.jl")
include("util.jl")
include("task.jl")
include("properties.jl")
include("analog.jl")
#include("digital.jl")
include("counter.jl")

const VERSION_SUPPORT = v"20.1.0"
global CLIENT_VERSION

const event_notify_c = Ref{Ptr{Cvoid}}(0)
const done_notify_c = Ref{Ptr{Cvoid}}(0)

function __init__()
    event_notify_c[] = @cfunction(event_notify, Cint,
                                  (TaskHandle, Cint, Cuint, Ptr{Cvoid}))
    done_notify_c[] = @cfunction(done_notify, Cint,
                                 (TaskHandle, Cint, Ptr{Cvoid}))
end

end

ext: jl


text:
function bariers_count(r::Robot)
    ans = 0
    move_while_can(r, (West, Sud)) # west sud
    dir = Ost
    ans += process_row(r, dir, Nord)
    while !isborder(r, Nord)     # идем змейкой
        move!(r, Nord)
        dir = inverse(dir)
        ans += process_row(r, dir, Nord)
    end
    move_while_can(r, (West, Nord))

    dir = Sud
    ans += process_row(r, dir, Ost)

    while !isborder(r, Ost)
        move!(r, Ost)
        dir = inverse(dir)
        ans += process_row(r, dir, Ost)
    end

    return ans - 2
end

function process_row(r::Robot, side, to_check_side) # checking on the side for barriers
    count = 0
    flag = false
    while !isborder(r, side) || bypass(r, side)
        if isborder(r, to_check_side)
            if !flag
                flag = true   #меняем флаг на тру и фолс если почитали перегородку то1 , пока есть перегородка будет тру,есил что то поенятеся все обнулится и появится новая перегородка и снова каунт +1
                count += 1
            end
        else
            flag = false
        end
        if !isborder(r, side)
            move!(r, side)
        end
    end
    return count
end

function bypass(r::Robot, dir)        # проверяем можно ли пройти в ту сторну если можно то возвращаем тру в ином случае фолс, всегда вовращаемся
    bypass_dir = left(dir)
    inverse_bypass_dir = inverse(bypass_dir)
    num_steps = 0
while !isborder(r, bypass_dir) && isborder(r, dir) # while no border side of passing
    move!(r, bypass_dir)
    num_steps += 1                                # counting the steps of pass
end
if isborder(r, dir) && isborder(r,bypass_dir)
    moves!(r, inverse_bypass_dir, num_steps)     # move in inverse side of passing
    return false
end
move!(r, dir)
moves!(r, inverse_bypass_dir, num_steps)        # in any case we should come back
return true
end

# move while it is possible for us
function move_while_can(r::Robot, sides)
    moving = true
    # it is possible to move
    while moving
        moving = false
        for side in sides
            if !isborder(r, side)
                move!(r, side)
                moving = true
            end
        end
    end
    end

###-------------------------------------------------------------------------
function moves!(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
    end
end
function moves!(r::Robot, side::HorizonSide, num_steps::Int)
    for _ in 1:num_steps
        move!(r, side)
    end
end

function putmarkers!(r::Robot, side::HorizonSide)
    while isborder(r, side) == false
        move!(r, side)
        putmarker!(r)
    end

function move_markers!(r::Robot, side::HorizonSide)
    while ismarker(r) == true
        move!(r, side)
    end
end

left(side::HorizonSide) = HorizonSide(mod(Int(side) + 1, 4))

inverse(side::HorizonSide) = HorizonSide(mod(Int(side) + 2, 4))

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
"""
Fundamental computation for CTMC
"""

"""
unifstep!(tr, P, poi, range, weight, x, y)

Compute the probability vector using the uniformized CTMC.

y = exp(tr(Q)*t) * x

where Q is unifomed by P = I - Q/qv. In the computation, Poisson p.m.f. with mean qv*t is used.

Parameters:
- tr: transpose operator
- P: The uniformed matrix
- poi: Poisson p.m.f.
- range: domain of Poisson p.m.f
- weight: The normalizing constant for Poisson p.m.f.
- x: Array (inout). x may be changed after executing
- y: Array (out). y should be zero before executing
Return value: nothing
"""

function unifstep!(tr::Symbol, P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv,
    x::ArrayT, y::ArrayT)::Nothing where {Ti,Tv,ArrayT<:AbstractArray{Tv}}
    _unifstep!(Val(tr), P, poi, range, weight, x, y)
end

@origin (poi => left) function _unifstep!(::Val{:T}, P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv,
    x::ArrayT, y::ArrayT)::Nothing where {Ti,Tv,ArrayT<:AbstractArray{Tv}}
    left, right = range
    Pdash = P'
    @axpy(poi[left], x, y)
    @inbounds for i = left+1:right
        x .= Pdash * x
        @axpy(poi[i], x, y)
    end
    @scal(1/weight, y)
    nothing
end

@origin (poi => left) function _unifstep!(::Val{:N}, P::AbstractMatrix{Tv},
    poi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv,
    x::ArrayT, y::ArrayT)::Nothing where {Ti,Tv,ArrayT<:AbstractArray{Tv}}
    left, right = range
    @axpy(poi[left], x, y)
    @inbounds for i = left+1:right
        x .= P * x
        @axpy(poi[i], x, y)
    end
    @scal(1/weight, y)
    nothing
end

"""
cunifstep!(tr, P, poi, cpoi, range, weight, qv_weight, x, y)

Compute the probability vector and cmulative one using the uniformized CTMC.

y = exp(tr(Q)*t) * x
cy = cy + int_0^t exp(tr(Q)*u) * x du

where Q is unifomed by P = I - Q/qv. In the computation, Poisson p.m.f. with mean qv*t is used.

Parameters:
- tr: transpose operator
- P: The uniformed matrix
- poi: Poisson p.m.f.
- cpoi: Poisson c.c.d.f.
- range: domain of Poisson p.m.f
- weight: The normalizing constant for Poisson p.m.f.
- qv_weight: The normalizing constant for Poisson c.c.d.f.
- x: Array (inout). x may be changed after executing
- y: Array (out). y should be zero before executing
- cy: Array (inout). cy should be zero before executing
Return value: nothing
"""

function cunifstep!(tr::Symbol, P::AbstractMatrix{Tv},
    poi::Vector{Tv}, cpoi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::ArrayT, y::ArrayT, cy::ArrayT)::Nothing where {Ti,Tv,ArrayT<:AbstractArray{Tv}}
    _cunifstep!(Val(tr), P, poi, cpoi, range, weight, qv_weight, x, y, cy)
end

@origin (poi => left, cpoi => left) function _cunifstep!(::Val{:T}, P::AbstractMatrix{Tv},
    poi::Vector{Tv}, cpoi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::ArrayT, y::ArrayT, cy::ArrayT)::Nothing where {Ti,Tv,ArrayT<:AbstractArray{Tv}}
    left, right = range
    Pdash = P'
    @axpy(poi[left], x, y)
    @axpy(cpoi[left], x, cy)
    @inbounds for i = left+1:right
        x .= Pdash * x
        @axpy(poi[i], x, y)
        @axpy(cpoi[i], x, cy)
    end
    @scal(1/weight, y)
    @scal(1/qv_weight, cy)
    nothing
end

@origin (poi => left, cpoi => left) function _cunifstep!(::Val{:N}, P::AbstractMatrix{Tv},
    poi::Vector{Tv}, cpoi::Vector{Tv}, range::Tuple{Ti,Ti}, weight::Tv, qv_weight::Tv,
    x::ArrayT, y::ArrayT, cy::ArrayT)::Nothing where {Ti,Tv,ArrayT<:AbstractArray{Tv}}
    left, right = range
    @axpy(poi[left], x, y)
    @axpy(cpoi[left], x, cy)
    @inbounds for i = left+1:right
        x .= P * x
        @axpy(poi[i], x, y)
        @axpy(cpoi[i], x, cy)
    end
    @scal(1/weight, y)
    @scal(1/qv_weight, cy)
    nothing
end


text:
module testModule

using Distributed
using GenGlobal
using StatsFuns
using SharedArrays

@GenGlobal globalx g_ysame2 g_ydiff2 g_sa
@GenGlobal myglobalx

export pplus, globf, update_shared!, compute_shared, remote_mapreduce, get_y2, myglobalx

pplus(x...) = broadcast(+, x...)

"each worker in `pids` returns `f()`, and results are reduced using `R()`"
remote_mapreduce(pids::AbstractVector{<:Integer}, f::Function, R::Function) = mapreduce((p::Int) -> remotecall_fetch(f, p), R, pids )

"each worker returns `f()`, and results are reduced using `R()`"
remote_mapreduce(f::Function, R::Function) = remote_mapreduce(workers(), f, R)

"set sa[i,j] = i ∀ i,j"
function update_shared!(sa::SharedArray{T,2}, j::Integer) where {T}
    # This is the inner function
    j ∈ 1:size(sa,2) || throw(DomainError())
    @views sa_vw = sa[:,j]
    for i in eachindex(sa_vw)
        sa_vw[i] = i
    end
end

# This is the "outer-wrapper" for above & ensures type-stability
function update_shared!(j::Integer)
    global g_sa
    update_shared!(g_sa, j)
end

"return `logsumexp(sa[:,j])` and `myid()`... as well as adding these to `ysame2` and `ydiff2`"
function compute_shared(sa::SharedArray{T,2}, j::Integer, ysame2::Matrix{T}, ydiff2::Matrix{T}) where {T}
    @views sa_vw = sa[:,j]
    ysame = logsumexp(sa_vw)
    ydiff = myid()
    ysame2 .+= ysame
    ydiff2 .+= ydiff
    return ysame, ydiff
end

function compute_shared(j::Integer)
    global g_sa
    global g_ysame2
    global g_ydiff2
    compute_shared(g_sa, j, g_ysame2, g_ydiff2)
end

"return global y2s as tuple from a worker"
function get_y2()
    # could also be done as get_y2() = (get_g_ysame2(), get_g_ydiff2(), )
    global g_ysame2
    global g_ydiff2
    return (g_ysame2, g_ydiff2,)
end

end # module end


text:
struct Where
    fn::Any
end

where(args...; kwargs...) = Where(args...; kwargs...)

struct Across
    innames::Any
    fn::Any
    outnames::Any

    Across(innames, fn::AbstractVector; names = nothing) = new(innames, fn, names)
    Across(innames, fn::Tuple; names = nothing) = new(innames, fn, names)
    Across(innames, fn::NamedTuple; names = nothing) = new(innames, fn, names)
    Across(innames, fn; names = "{col}") = new(innames, (fn,), names)

end

across(args...; kwargs...) = Across(args...; kwargs...)


text:
module Utils

export batched_outer, batched_kronecker, add_dim, gram_matrix

using Flux

"""
# Parameters

- `A`: array of shape (`m`, `num_batch`)
- `B`: array of shape (`n`, `num_batch`)

# Return

- `C`: array of shape (`m`, `n`, `num_batch`), where `C[:, :, i] = A[:, i] * B[:, i]'`

"""
function batched_outer(A, B)
    @assert ndims(A) == 2
    @assert ndims(B) == 2
    size(A, 2) == size(B, 2) || throw(DimensionMismatch("batch size mismatch"))
    num_batch = size(A, 2)
    A′ = reshape(A, :, 1, num_batch)
    B′ = reshape(B, 1, :, num_batch)
    Flux.batched_mul(A′, B′)
end

"""
mat1: (m, n, batch)
mat2: (j, k, batch)
"""
function batched_kronecker(mat1, mat2)
    @assert ndims(mat1) == 3
    @assert ndims(mat2) == 3
    @assert size(mat1, 3) == size(mat2, 3)
    num_batch = size(mat1, 3)
    mat1_flatten = reshape(mat1, :, 1, num_batch)
    mat2_flatten = reshape(mat2, 1, :, num_batch)
    curried_reshape(new_shape...) = (m)-> reshape(m, new_shape...)
    curried_permdims(perms...) = (m) -> permutedims(m, perms)
    Flux.batched_mul(
        mat1_flatten,
        mat2_flatten,
        ) |>
    curried_reshape(size(mat1)[1:end-1]..., size(mat2)[1:end-1]..., num_batch) |>
    curried_permdims(3, 1, 4, 2, 5) |>
    curried_reshape(size(mat1, 1)*size(mat2, 1), size(mat1, 2)*size(mat2, 2), num_batch)
end

"""
https://stackoverflow.com/a/59849291/3908797
"""
add_dim(x::AbstractArray{T, N}) where {T, N} = reshape(x, Val(N+1))

"""
- `A`: (d, m, batch)
- `B`: (d, n, batch)

return C: (m, n, batch)
"""
gram_matrix(A, B) = Flux.batched_mul(permutedims(A, (2, 1, 3)), B)

"""
`labels`: (m, batch)
"""
function onehot(labels)
    num_batch = size(labels, 2)
    encoded = []
    for i in 1:num_batch
        uniq_label = unique(labels[:, i])
        encoder(l) = Flux.onehot(l, uniq_label)
        push!(
            encoded,
            vcat(encoder.(labels[:, i])...)
        )
    end
    hcat(encoded...)
end

end


text:
function update_alpha!(s::State, c::Constants, d::Data, sb_ibp::Bool)
  newShape = shape(c.alpha_prior) + c.K

  if sb_ibp
    newRate = rate(c.alpha_prior) - sum(log.(s.v))
  else
    newRate = rate(c.alpha_prior) - sum(log.(s.v)) / c.K
  end

  post = Gamma(newShape, 1.0 / newRate)
  s.alpha = rand(post)
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using Plots

function circle(x)
    θ = x[1];
    r = 1/2 + 1/2 * (2*x[2]/(1+x[2]^2));
    return [1-r*cos(θ), 1-r*sin(θ)];
end


text:
module PHcalc
using Optim

export Neutral,Acid, System, α,pHsolve, minimise, pHfast

mutable struct Acid
	ka::Vector{Float64}
	conc::Float64
	charge::Vector{Int64}
	pKA::Vector{Float64}
	_ka::Vector{Float64}
	function Acid(ka, conc, charge::Int64)
			if typeof(ka) <: Number
				ka=[ka]
			end
			_ka=copy(ka)
			sort!(ka, rev=:true)
			prepend!(_ka,1.)
			charge=charge:-1:charge-(length(ka)) |> collect
			# charge=LinRange(charge,-length(ka), length(ka))
			new(ka, conc,charge , -log10.(ka),_ka )
		end
end

mutable struct Neutral
	charge::Union{Vector{Int64},Int64}
	conc::Float64
end

struct System
	species
	function System(species...)
		new(species)
	end
end

function α(self::Acid,pH) # solo para un valor de pH, se podria pensar para un array de pHs.
		h3o=10.0^(-pH)
		power=length(self._ka)-1:-1:0
		h3o_power=(x->h3o^x).(power)
		Ka_prod=cumprod(self._ka)
		h3o_Ka=h3o_power.*(Ka_prod)
		h3o_Ka./sum(h3o_Ka)
	end

function α(self::Neutral,pH) # Devielve 1 porque no tiene disociacion, devuelvo en array para que se parezca a la otra implementación.
	[1]
end

function pHsolve(sys)
	function minimise(pH)
		h3o=10.0^(-pH)
     	oh = (10.0^(-14))/h3o
	 	x = (h3o - oh)
		for specie in sys.species
			x+=(specie.conc.*specie.charge.*α(specie,pH))|> sum
		end
	 	abs(x)
	end

	function phguess()
		phv=1:.1:14
		minimise.(phv) |> argmin |> x -> getindex(phv, x)
	end
	guess=phguess()
	optimize(x->minimise(first(x)), [guess], BFGS()).minimizer |> first
end

function minimise(sys::System,pH)
	h3o=10.0^(-pH)
     	oh = (10.0^(-14))/h3o
	 x = (h3o - oh)
	for specie in sys.species
		x+=(specie.conc.*specie.charge.*α(specie,pH))|> sum
	end
	 abs(x)
end

function pHfast(sys, precision=.01)
	function minimise(pH)
		h3o=10.0^(-pH)
     	oh = (10.0^(-14))/h3o
	 	x = (h3o - oh)
		for specie in sys.species
			x+=(specie.conc.*specie.charge.*α(specie,pH))|> sum
		end
	 	abs(x)
	end
	phv=1:precision:14
	minimise.(phv) |> argmin |> x -> getindex(phv,x)

end

end


text:
module Exceptions

import Genie
import HTTP

export ExceptionalResponse, RuntimeException, InternalServerException, NotFoundException, FileExistsException

"""
    struct ExceptionalResponse <: Exception

A type of exception which wraps an HTTP Response object.
The thrown exception will propagate until it is caught up the app stack or ultimately by Genie
and the wrapped response is sent to the client.

### Example
If the user is not authenticated, an `ExceptionalResponse` is thrown - if the exception is not caught
in the app's stack, Genie will catch it and return the wrapped `Response` object, forcing an HTTP redirect to the login page.

```julia
isauthenticated() || throw(ExceptionalResponse(redirect(:show_login)))
```
"""
struct ExceptionalResponse <: Exception
  response::HTTP.Response
end

Base.show(io::IO, ex::ExceptionalResponse) = print(io, "ExceptionalResponseException: $(ex.response.status) - $(Dict(ex.response.headers))")

###

"""
    RuntimeException

Represents an unexpected and unhandled runtime exceptions. An error event will be logged and the
exception will be sent to the client, depending on the environment
(the error stack is dumped by default in dev mode or an error message is displayed in production).

It allows defining custom error message and info, as well as an error code, in addition to the exception object.

# Arguments
- `message::String`
- `info::String`
- `code::Int`
- `ex::Union{Nothing,Exception}`
"""
struct RuntimeException <: Exception
  message::String
  info::String
  code::Int
  ex::Union{Nothing,Exception}
end

"""
    RuntimeException(message::String, code::Int)

`RuntimeException` constructor using `message` and `code`.
"""
RuntimeException(message::String, code::Int) = RuntimeException(message, "", code, nothing)

"""
    RuntimeException(message::String, info::String, code::Int)

`RuntimeException` constructor using `message`, `info` and `code`.
"""
RuntimeException(message::String, info::String, code::Int) = RuntimeException(message, info, code, nothing)

"""
    Base.show(io::IO, ex::RuntimeException)

Custom printing of `RuntimeException`
"""
Base.show(io::IO, ex::RuntimeException) = print(io, ex.message)

###

"""
    struct InternalServerException <: Exception

Dedicated exception type for server side exceptions. Results in a 500 error by default.

# Arguments
- `message::String`
- `info::String`
- `code::Int`
"""
struct InternalServerException <: Exception
  message::String
  info::String
  code::Int
end

"""
    InternalServerException(message::String)

External `InternalServerException` constructor accepting a custome message.
"""
InternalServerException(message::String) = InternalServerException(message, "", 500)

"""
    InternalServerException()

External `InternalServerException` using default values.
"""
InternalServerException() = InternalServerException("Internal Server Error")

###

"""
    struct NotFoundException <: Exception

Specialized exception representing a not found resources. Results in a 404 response being sent to the client.

# Arguments
- `message::String`
- `info::String`
- `code::Int`
- `resource::String`
"""
struct NotFoundException <: Exception
  message::String
  info::String
  code::Int
  resource::String
end

Base.show(io::IO, ex::NotFoundException) = print(io, ex.message)

"""
    NotFoundException(resource::String)

External constructor allowing to pass the name of the not found resource.
"""
NotFoundException(resource::String) = NotFoundException("not found", "", 404, resource)

NotFoundException(resource::String, info::String) = NotFoundException("$resource can not be found", info, 404, resource)

"""
    NotFoundException()

External constructor using default arguments.
"""
NotFoundException() = NotFoundException("")

###

"""
    struct FileExistsException <: Exception

Custom exception type for signaling that the requested file already exists.
"""
struct FileExistsException <: Exception
  path::String
end

"""
    Base.show(io::IO, ex::FileExistsException)

Custom printing for `FileExistsException`
"""
Base.show(io::IO, ex::FileExistsException) = print(io, "FileExistsException: $(ex.path)")

end


text:
module MyLibPrecomp

greet() = print("Hello World!")

add(n::Int, m::Int) = n + m

end # module

max_stars_repo_name: RevoltBI/Julia-Meetup

max_issues_repo_name: RevoltBI/Julia-Meetup

max_forks_repo_name: RevoltBI/Julia-Meetup


text:
using ContactStructures

plotlyjs()
plot_along_param_curve(t -> (2*cos(2π*t),2*sin(2π*t),0), standard_structure, 0.0625;  midpoint_marker_color=:red, ε=0.35)
plot!(xlims=(-3, 3), ylims=(-3,3), zlims=(-3,3),legend=false)


text:
export
ruleSPGaussianWeightedMeanPrecisionOutNPP,
ruleVBGaussianWeightedMeanPrecisionOut

ruleSPGaussianWeightedMeanPrecisionOutNPP(  msg_out::Nothing,
                                            msg_weighted_mean::Message{PointMass, V},
                                            msg_prec::Message{PointMass}) where V<:VariateType =
    Message(V, GaussianWeightedMeanPrecision, xi=deepcopy(msg_weighted_mean.dist.params[:m]), w=deepcopy(msg_prec.dist.params[:m]))

ruleVBGaussianWeightedMeanPrecisionOut( dist_out::Any,
                                        dist_weighted_mean::ProbabilityDistribution{V},
                                        dist_prec::ProbabilityDistribution) where V<:VariateType =
    Message(V, GaussianWeightedMeanPrecision, xi=unsafeMean(dist_weighted_mean), w=unsafeMean(dist_prec))


text:
struct PlaceHolder end

import Base

□ = PlaceHolder()

Base.show(io::IO, ::PlaceHolder) = print(io, "□")


text:
module Distributions

using Stats

export                                  # types
    CauchitLink,
    CloglogLink,
    Distribution,
    UnivariateDistribution,
    MultivariateDistribution,
    MatrixDistribution,
    NonMatrixDistribution,
    DiscreteDistribution,
    ContinuousDistribution,
    DiscreteUnivariateDistribution,
    DiscreteMultivariateDistribution,
    DiscreteMatrixDistribution,
    ContinuousUnivariateDistribution,
    ContinuousMultivariateDistribution,
    ContinuousMatrixDistribution,
    Arcsine,
    Bernoulli,
    Beta,
    BetaPrime,
    Binomial,
    Categorical,
    Cauchy,
    Chisq,
    Dirichlet,
    DiscreteUniform,
    DoubleExponential,
    EmpiricalDistribution,
    Erlang,
    Exponential,
    FDist,
    Gamma,
    Geometric,
    Gumbel,
    HyperGeometric,
    IdentityLink,
    InverseLink,
    InverseWishart,
    InvertedGamma,
    Laplace,
    Levy,
    Link,
    Logistic,
    LogitLink,
    LogLink,
    logNormal,
    MixtureModel,
    Multinomial,
    MultivariateNormal,
    NegativeBinomial,
    NoncentralBeta,
    NoncentralChisq,
    NoncentralF,
    NoncentralT,
    Normal,
    Pareto,
    Poisson,
    ProbitLink,
    Rayleigh,
    TDist,
    Triangular,
    TruncatedNormal,
    Uniform,
    Weibull,
    Wishart,
                                        # methods
    binaryentropy, # entropy of distribution in bits
    canonicallink, # canonical link function for a distribution
    ccdf,          # complementary cdf, i.e. 1 - cdf
    cdf,           # cumulative distribution function
    cf,            # characteristic function
    cgf,           # cumulant generating function
    cquantile,     # complementary quantile (i.e. using prob in right hand tail)
    cumulant,      # cumulants of distribution
    deviance,      # deviance of fitted and observed responses
    devresid,      # vector of squared deviance residuals
    entropy,       # entropy of distribution in nats
    fit,           # fit a distribution to data
    freecumulant,  # free cumulants of distribution
    insupport,     # predicate, is x in the support of the distribution?
    invlogccdf,    # complementary quantile based on log probability
    invlogcdf,     # quantile based on log probability
    kde,           # Kernel density estimator
    kurtosis,      # kurtosis of the distribution
    linkfun,       # link function mapping mu to eta, the linear predictor
    linkinv,       # inverse link mapping eta to mu
    logccdf,       # ccdf returning log-probability
    logcdf,        # cdf returning log-probability
    logpdf,        # log probability density
    logpdf!,       # evaluate log pdf to provided storage
    logpmf,        # log probability mass
    logpmf!,       # evaluate log pmf to provided storage
    mean,          # mean of distribution
    median,        # median of distribution
    mgf,           # moment generating function
    modes,         # mode(s) of distribution as vector
    moment,        # moments of distribution
    mueta,         # derivative of inverse link function
    mustart,       # starting values of mean vector in GLMs
    pdf,           # probability density function (ContinuousDistribution)
    pmf,           # probability mass function (DiscreteDistribution)
    quantile,      # inverse of cdf (defined for p in (0,1))
    rand,          # random sampler
    rand!,         # replacement random sampler
    sample,        # another random sampler - not sure why this is here
    skewness,      # skewness of the distribution
    sprand,        # random sampler for sparse matrices
    std,           # standard deviation of distribution
    valideta,      # validity check on linear predictor
    validmu,       # validity check on mean vector
    var            # variance of distribution

import Base.mean, Base.median, Base.quantile
import Base.rand, Base.std, Base.var, Base.cor, Base.cov
import Base.show, Base.sprand

include("drawtable.jl")
include("tvpack.jl")

abstract Distribution
abstract UnivariateDistribution             <: Distribution
abstract MultivariateDistribution           <: Distribution
abstract MatrixDistribution                 <: Distribution

abstract DiscreteUnivariateDistribution     <: UnivariateDistribution
abstract ContinuousUnivariateDistribution   <: UnivariateDistribution

abstract DiscreteMultivariateDistribution   <: MultivariateDistribution
abstract ContinuousMultivariateDistribution <: MultivariateDistribution

abstract ContinuousMatrixDistribution       <: MatrixDistribution
abstract DiscreteMatrixDistribution         <: MatrixDistribution

typealias NonMatrixDistribution Union(UnivariateDistribution, MultivariateDistribution)
typealias DiscreteDistribution Union(DiscreteUnivariateDistribution, DiscreteMultivariateDistribution)
typealias ContinuousDistribution Union(ContinuousUnivariateDistribution, ContinuousMultivariateDistribution)

include("fallbacks.jl")
include("rmath.jl")

# TODO: Move these two methods into Stats
xlogx(x::Real) = x == 0.0 ? 0.0 : x * log(x)
xlogxdmu(x::Real, mu::Real) = x == 0.0 ? 0.0 : x * log(x / mu)

# Univariate distributions
include(joinpath("univariate", "alpha.jl"))
include(joinpath("univariate", "arcsine.jl"))
include(joinpath("univariate", "bernoulli.jl"))
include(joinpath("univariate", "beta.jl"))
include(joinpath("univariate", "betaprime.jl"))
include(joinpath("univariate", "binomial.jl"))
include(joinpath("univariate", "categorical.jl"))
include(joinpath("univariate", "cauchy.jl"))
include(joinpath("univariate", "chi.jl"))
include(joinpath("univariate", "chisq.jl"))
include(joinpath("univariate", "discreteuniform.jl"))
include(joinpath("univariate", "exponential.jl"))
include(joinpath("univariate", "fdist.jl"))
include(joinpath("univariate", "gamma.jl"))
include(joinpath("univariate", "erlang.jl"))
include(joinpath("univariate", "geometric.jl"))
include(joinpath("univariate", "gumbel.jl"))
include(joinpath("univariate", "hypergeometric.jl"))
include(joinpath("univariate", "invertedgamma.jl"))
include(joinpath("univariate", "laplace.jl"))
include(joinpath("univariate", "logistic.jl"))
include(joinpath("univariate", "lognormal.jl"))
include(joinpath("univariate", "negativebinomial.jl"))
include(joinpath("univariate", "noncentralbeta.jl"))
include(joinpath("univariate", "noncentralchisq.jl"))
include(joinpath("univariate", "noncentralf.jl"))
include(joinpath("univariate", "noncentralt.jl"))
include(joinpath("univariate", "normal.jl"))
include(joinpath("univariate", "pareto.jl"))
include(joinpath("univariate", "poisson.jl"))
include(joinpath("univariate", "rayleigh.jl"))
include(joinpath("univariate", "tdist.jl"))
include(joinpath("univariate", "triangular.jl"))
include(joinpath("univariate", "uniform.jl"))
include(joinpath("univariate", "weibull.jl"))

# Multivariate distributions
include(joinpath("multivariate", "dirichlet.jl"))
include(joinpath("multivariate", "multinomial.jl"))
include(joinpath("multivariate", "multivariatenormal.jl"))

# Matrix distributions
include(joinpath("matrix", "inversewishart.jl"))
include(joinpath("matrix", "wishart.jl"))

# Truncated distributions
include("truncate.jl")
include(joinpath("univariate", "truncated", "normal.jl"))

# Mixture distributions
include("mixturemodel.jl")

# Sample w/ and w/o replacement
include("sample.jl")

# Link functions for GLM's
include("glmtools.jl")

# REPL representations
include("show.jl")

# Fit distributions to data
include("fit.jl")

# Kernel density estimators
include("kde.jl")

# Expectations, entropy, KL divergence
include("functionals.jl")

end # module


text:
using Distributed
#addprocs(5)

using LinearAlgebra
using GeoPhyInv
using Test
using Misfits
using BenchmarkTools
using Profile
using LinearMaps

"""
* `files` are in folder
* `gfiles` are files in the generated folder
"""
function run_test(files,folder,gfiles=[])
	for file in files
		fp = joinpath(folder, string(file, ".jl"))
		println(" ***********************************************")
		println(" *********** $(fp) ***************")
		println(" ***********************************************")
		include(fp)
	end
	for file in gfiles
		fpg = joinpath("generated",folder,string(file, ".jl"))
		println(" ***********************************************")
		println(" ******** $(fpg) *********")
		println(" ***********************************************")
		include(fpg)
	end
end

run_test(["convert", "param_adj"],"media", ["doc"])
run_test([],"ageom", ["doc"])
run_test([],"srcwav", ["doc"])
run_test([],"data", ["doc"])
run_test(["accuracy","backprop", "rho_projection"],"fdtd", ["doc","create_snaps","reuse_expt"])
run_test(["gradient_accuracy", "born_map"],"fwi", ["doc", "pizza", "born_tutorial"])
run_test(["interp_tests"],"Interpolation", [])
run_test(["testscript_RandomEigenfns", "adj_state_expt", "adj_state", "testdAdx"], "Poisson", ["doc","forw", "test_born"])


text:
using Distributions
using Base.Test

## Fisher's noncentral hypergeometric distribution
ns = 80
nf = 60
n = 100

# http://en.wikipedia.org/wiki/Fisher's_noncentral_hypergeometric_distribution
ω = 10.0
d = FisherNoncentralHypergeometric(ns, nf, n, ω)

@test_approx_eq_eps mean(d) 71.95759 1e-5
@test mode(d) == 72
@test quantile(d, 0.05) == 68
@test quantile(d, 0.95) == 75
@test_approx_eq_eps logpdf(d, 75) -2.600957 1e-6
@test_approx_eq_eps pdf(d, 75) 0.07420252 1e-7
@test_approx_eq_eps cdf(d, 75) 0.9516117 1e-7

ω = 0.361
d = FisherNoncentralHypergeometric(ns, nf, n, ω)

@test_approx_eq_eps mean(d) 50.42299 1e-5
@test mode(d) == 50
@test quantile(d, 0.05) == 47
@test quantile(d, 0.95) == 54
@test_approx_eq_eps logpdf(d, 50) -1.820071 1e-6
@test_approx_eq_eps pdf(d, 50) 0.1620142 1e-7
@test_approx_eq_eps cdf(d, 50) 0.5203835 1e-7

# distributions are both equal to the (central) hypergeometric distribution when the odds ratio is 1.
ω = 1.0
d = FisherNoncentralHypergeometric(ns, nf, n, ω)

ref = Hypergeometric(ns,nf,n)
@test_approx_eq_eps mean(d) mean(ref) 1e-5
@test_approx_eq_eps logpdf(d, 51) logpdf(ref, 51) 1e-7
@test_approx_eq_eps pdf(d, 51) pdf(ref, 51) 1e-7
@test_approx_eq_eps cdf(d, 51) cdf(ref, 51) 1e-7
@test_approx_eq_eps quantile(d, 0.05) quantile(ref, 0.05) 1e-7
@test_approx_eq_eps quantile(d, 0.95) quantile(ref, 0.95) 1e-7
@test mode(d) == mode(ref)

## Wallenius' noncentral hypergeometric distribution
ns = 80
nf = 60
n = 100

# http://en.wikipedia.org/wiki/Fisher's_noncentral_hypergeometric_distribution
ω = 10.0
d = WalleniusNoncentralHypergeometric(ns, nf, n, ω)

@test_approx_eq_eps mean(d) 78.82945 1e-5
@test mode(d) == 80
@test quantile(d, 0.05) == 77
@test quantile(d, 0.95) == 80
@test_approx_eq_eps logpdf(d, 75) -4.750073 1e-6
@test_approx_eq_eps pdf(d, 75) 0.00865106 1e-7
@test_approx_eq_eps cdf(d, 75) 0.01133378 1e-7

ω = 0.361
d = WalleniusNoncentralHypergeometric(ns, nf, n, ω)

@test_approx_eq_eps mean(d) 45.75323 1e-5
@test mode(d) == 45
@test quantile(d, 0.05) == 42
@test quantile(d, 0.95) == 49
@test_approx_eq_eps logpdf(d, 50) -3.624478 1e-6
@test_approx_eq_eps pdf(d, 50) 0.026663 1e-7
@test_approx_eq_eps cdf(d, 50) 0.983674 1e-7

# distributions are both equal to the (central) hypergeometric distribution when the odds ratio is 1.
ω = 1.0
d = WalleniusNoncentralHypergeometric(ns, nf, n, ω)

ref = Hypergeometric(ns,nf,n)
@test_approx_eq_eps mean(d) mean(ref) 1e-5
@test_approx_eq_eps logpdf(d, 51) logpdf(ref, 51) 1e-7
@test_approx_eq_eps pdf(d, 51) pdf(ref, 51) 1e-7
@test_approx_eq_eps cdf(d, 51) cdf(ref, 51) 1e-7
@test_approx_eq_eps quantile(d, 0.05) quantile(ref, 0.05) 1e-7
@test_approx_eq_eps quantile(d, 0.95) quantile(ref, 0.95) 1e-7
@test mode(d) == mode(ref)