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

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

    """
    `struct  PhotoDoubleIonization.Settings`  ... defines a type for the settings in estimating single-photon double-ionization cross sections.

        + multipoles              ::Array{EmMultipole}      ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}         ... Specifies the gauges to be included into the computations.
        + photonEnergies          ::Array{Float64,1}        ... List of photon energies.
        + green                   ::Array{GreenChannel,1}   ... Precalculated and user-specified Green function of the ion.
        + NoEnergySharings        ::Int64                   ... Number of energy sharings that are used in the computations for each line.
        + calcAnisotropy          ::Bool                    ... True, if the beta anisotropy parameters are to be calculated and false o/w.
        + printBefore             ::Bool                    ... True, if all energies and lines are printed before their evaluation.
        + maxKappa                ::Int64                   ... Maximum kappa value of partial waves to be included.
        + lineSelection           ::LineSelection           ... Specifies the selected levels, if any.
    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        photonEnergies            ::Array{Float64,1}
        green                     ::Array{GreenChannel,1}
        NoEnergySharings          ::Int64
        calcAnisotropy            ::Bool
        printBefore               ::Bool
        maxKappa                  ::Int64
        lineSelection             ::LineSelection
    end

    """
    `PhotoDoubleIonization()`  ... constructor for the default values of PhotoDoubleIonization line computations
    """
    function Settings()
        Settings(EmMultipole[E1], UseGauge[UseCoulomb, UseBabushkin], Float64[], GreenChannel[], 0, false, false, 0, LineSelection())
    end

    # `Base.show(io::IO, settings::PhotoDoubleIonization.Settings)`  ... prepares a proper printout of the variable settings::PhotoDoubleIonization.Settings.
    function Base.show(io::IO, settings::PhotoDoubleIonization.Settings)
        println(io, "multipoles:                   $(settings.multipoles)  ")
        println(io, "gauges:                       $(settings.gauges)  ")
        println(io, "photonEnergies:               $(settings.photonEnergies)  ")
        println(io, "green:                         (settings.green)  ")
        println(io, "NoEnergySharings:             $(settings.NoEnergySharings)  ")
        println(io, "calcAnisotropy:               $(settings.calcAnisotropy)  ")
        println(io, "printBefore:                  $(settings.printBefore)  ")
        println(io, "lineSelection:                $(settings.lineSelection)  ")
    end

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

        + multipole      ::EmMultipole          ... Multipole of the photon absorption.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + omega          ::Float64              ... photon energy
        + energy1        ::Float64              ... energy of the partial wave_1
        + kappa1         ::Int64                ... partial-wave of the free electron_1
        + phase1         ::Float64              ... phase of the partial wave_1
        + xSymmetry      ::LevelSymmetry        ... angular momentum and parity after coupling electron1.
        + energy2        ::Float64              ... energy of the partial wave_2
        + kappa2         ::Int64                ... partial-wave of the free electron_2
        + phase2         ::Float64              ... phase of the partial wave_2
        + tSymmetry      ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + amplitude      ::Complex{Float64}     ... Auger amplitude associated with the given channel.
    """
    struct  ReducedChannel
        multipole        ::EmMultipole
        gauge            ::EmGauge
        omega            ::Float64
        energy1          ::Float64
        kappa1           ::Int64
        phase1           ::Float64
        xSymmetry        ::LevelSymmetry
        energy2          ::Float64
        kappa2           ::Int64
        phase2           ::Float64
        tSymmetry        ::LevelSymmetry
        amplitude        ::Complex{Float64}
    end

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

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

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

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

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

        + initialLevel   ::Level            ... initial-(state) level
        + finalLevel     ::Level            ... final-(state) level
        + omega          ::Float64          ... photon energy of this line
        + totalCs        ::EmProperty       ... Total rate of this line.
        + sharings       ::Array{PhotoDoubleIonization.Sharing,1}  ... List of PhotoDoubleIonization sharings of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        omega            ::Float64
        totalCs          ::EmProperty
        sharings         ::Array{PhotoDoubleIonization.Sharing,1}
    end

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

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

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

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

    """
    `PhotoDoubleIonization.determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, omega::Float64, energy::Float64,
                                                        settings::PhotoDoubleIonization.Settings)`
        ... to determine a list of PhotoDoubleIonization Sharing's and Channel's for a transitions from the initial to final level and by taking into
            account the particular settings of for this computation; an Array{PhotoDoubleIonization.Sharing,1} is returned.
    """
    function determineSharingsAndChannels(finalLevel::Level, initialLevel::Level, omega::Float64, energy::Float64, settings::PhotoDoubleIonization.Settings)
        sharings = PhotoDoubleIonization.Sharing[];    eSharings = Basics.determineEnergySharings(energy, settings.NoEnergySharings)
        for  es in eSharings
            epsilon1  = es[1];    epsilon2 = es[2];    weight = es[3]
            channels  = PhotoDoubleIonization.ReducedChannel[];
            symi      = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
            for  mp in settings.multipoles
                for  gauge in settings.gauges
                    symList = AngularMomentum.allowedMultipoleSymmetries(symi, mp)
                    for  symt in symList
                        couplings = AngularMomentum.allowedDoubleKappaCouplingSequence(symf, symt, settings.maxKappa)
                        for  coupling  in  couplings
                            # Include further restrictions if appropriate
                            if     string(mp)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                                push!(channels, ReducedChannel(mp, Basics.Coulomb,   omega, epsilon1, coupling[1], 0., coupling[2],
                                                                                epsilon2, coupling[3], 0., symt, Complex(0.)) )
                            elseif string(mp)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                                push!(channels, ReducedChannel(mp, Basics.Babushkin, omega, epsilon1, coupling[1], 0., coupling[2],
                                                                                epsilon2, coupling[3], 0., symt, Complex(0.)) )
                            elseif string(mp)[1] == 'M'
                                push!(channels, ReducedChannel(mp, Basics.Magnetic,  omega, epsilon1, coupling[1], 0., coupling[2],
                                                                                epsilon2, coupling[3], 0., symt, Complex(0.)) )
                            end
                        end
                    end
                end
            end
            push!(sharings, PhotoDoubleIonization.Sharing(omega, epsilon1, epsilon2, weight, EmProperty(0., 0.), channels) )
        end
        return( sharings )
    end

    """
    `PhotoDoubleIonization.displayDifferentialCrossSections(stream::IO, lines::Array{PhotoDoubleIonization.Line,1}, settings::PhotoDoubleIonization.Settings)`
        ... to display all differential rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayDifferentialCrossSections(stream::IO, lines::Array{PhotoDoubleIonization.Line,1}, settings::PhotoDoubleIonization.Settings)
        #
        # First, print lines and sharings
        nx = 130
        println(stream, " ")
        println(stream, "  Energy-differential cross sections of selected photo-double ionization lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(38, "Energies (all in " * TableStrings.inUnits("energy") * ")"; na=3);
        sb = sb * TableStrings.flushleft(38, "  i -- f        epsilon_1     epsilon_2"; na=3)
        sa = sa * TableStrings.center(14, "Weight"; na=0);                            sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(34, "Cou -- diff. cross section -- Bab"; na=3)
        sb = sb * TableStrings.center(34, TableStrings.inUnits("cross section") * "        " *
                                          TableStrings.inUnits("cross section"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            for  (is, sharing)  in  enumerate(line.sharings)
                if  is == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon1))                        * "    "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon2))                        * "    "
                sb = sb * @sprintf("%.4e",                                              sharing.weight)                           * "       "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("cross section: from atomic", sharing.differentialCs.Coulomb))   * "   "
                sb = sb * @sprintf("%.5e", Defaults.convertUnits("cross section: from atomic", sharing.differentialCs.Babushkin)) * "   "
                println(stream,  sb )
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `PhotoDoubleIonization.displayLines(lines::Array{PhotoDoubleIonization.Line,1})`
        ... to display a list of lines, sharings and channels that have been selected due to the prior settings. A neat table
            of all selected transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{PhotoDoubleIonization.Line,1})
        #
        # First, print lines and sharings
        nx = 94
        println(" ")
        println("  Selected photo-double ionization lines & sharings:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(54, "Energies (all in )" * TableStrings.inUnits("energy") * ")"; na=5);
        sb = sb * TableStrings.flushleft(54, "  i -- f         omega        epsilon_1    epsilon_2"; na=5)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            for  (is, sharing)  in  enumerate(line.sharings)
                if  is == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.omega))    * "    "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon1)) * "   "
                sb = sb * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", sharing.epsilon2)) * "   "
                println( sb )
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        #
        # Second, print lines and channles
        nx = 150
        println(" ")
        println("  Selected photo-double ionization lines & channels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.flushleft(100, "Channels (all energies in " * TableStrings.inUnits("energy") * ")" ; na=5);
        sb = sb * TableStrings.flushleft(100, "Multipole  Gauge      omega       (epsilon, kappa)_1   -->    J^P_x   -->  (epsilon, kappa)_2   -->    J^P_t"; na=5)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            for (ic, ch) in enumerate(line.sharings[1].channels)
                if  ic == 1     sb = sa     else    sb = TableStrings.hBlank( length(sa) )    end
                skappa1 = "  " * string(ch.kappa1);         skappa2 = "  " * string(ch.kappa2);
                sxsym   = "      " * string(ch.xSymmetry);  stsym   = "      " * string(ch.tSymmetry)
                sb = sb * string(ch.multipole) * "         " * string(ch.gauge)[1:3] * "      "
                sb = sb * @sprintf("%.3e", Defaults.convertUnits("energy: from atomic", ch.omega)) * "     "
                sb = sb * " (" * @sprintf("%.3e", Defaults.convertUnits("energy: from atomic", ch.energy1)) * ", " * skappa1[end-2:end] * ")    -->"
                sb = sb * sxsym[end-8:end] * "   -->  "
                sb = sb * " (" * @sprintf("%.3e", Defaults.convertUnits("energy: from atomic", ch.energy2)) * ", " * skappa2[end-2:end] * ")    -->"
                sb = sb * stsym[end-8:end]
                println( sb )
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `PhotoDoubleIonization.displayTotalCrossSections(stream::IO, lines::Array{PhotoDoubleIonization.Line,1}, settings::PhotoDoubleIonization.Settings)`
        ... to display all total rates, etc. of the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayTotalCrossSections(stream::IO, lines::Array{PhotoDoubleIonization.Line,1}, settings::PhotoDoubleIonization.Settings)
        #
        # First, print lines and sharings
        nx = 103
        println(stream, " ")
        println(stream, "  Total rates of selected photo-double ionization lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=2);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(12, "f--Energy--i"; na=3)
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(12, "omega"     ; na=4)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(30, "Cou --   cross section   -- Bab"; na=4)
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section") * "        " *
                                          TableStrings.inUnits("cross section"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";      isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic", energy)) * "    "
            #
            sb = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.omega))                    * "         "
            sb = sb * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", line.totalCs.Coulomb))   * "     "
            sb = sb * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", line.totalCs.Babushkin)) * "   "
            println(stream,  sb )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
using Documenter, FeatherFiles

makedocs(
	modules = [FeatherFiles],
	sitename = "FeatherFiles.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/FeatherFiles.jl.git"
)

ext: jl


text:
mutable struct CuContext
    ptr::Ptr{Cvoid}
end

function CuContext(dev::Int)
    ref = Ref{Ptr{Cvoid}}()
    @apicall :cuCtxCreate (Ptr{Ptr{Cvoid}},Cuint,Cint) ref 0 dev
    ctx = CuContext(ref[])
    # finalizer(ctx, destroy)
    ctx
end

Base.:(==)(a::CuContext, b::CuContext) = a.ptr == b.ptr
Base.hash(ctx::CuContext, h::UInt) = hash(ctx.ptr, h)
Base.unsafe_convert(::Type{Ptr{Cvoid}}, ctx::CuContext) = ctx.ptr

function destroy(ctx::CuContext)
    setcontext(ctx)
    @apicall :cuCtxDestroy (Ptr{Cvoid},) ctx
end

function getcontext()
    ref = Ref{Ptr{Cvoid}}()
    @apicall :cuCtxGetCurrent (Ptr{Ptr{Cvoid}},) ref
    CuContext(ref[])
end

function setcontext(ctx::CuContext)
    @apicall :cuCtxSetCurrent (Ptr{Cvoid},) ctx
end
function setcontext(f::Function, ctx::CuContext)
    oldctx = getcontext()
    setcontext(ctx)
    f()
    setcontext(oldctx)
end

function synchronize()
    @apicall :cuCtxSynchronize ()
    empty!(ALLOCATED)
    empty!(CUDNN.ALLOCATED)
end


text:
module Tester

# stdlib dependencies
using Random
using Logging

# external dependencies
using JTAGManager
using ProgressMeter

export JTAG, led, validate

# Set start and stop addresses of DRAM.
dram_addr_start() = 0
dram_addr_end() = 0x03ff_ffff

# Write to the LED of the board
# LED's are at address 0x4000_0000
"""
    led(jtag::JTAG, val)

Write `val` to the LEDs on the DE10 Board.
"""
led(jtag::JTAG, val) = write(jtag, 0x4000_0000, [val])

# Struct for tracking what has been written and where.
mutable struct AddressSlot
    address :: Int
    value :: UInt8
    iswritten :: Bool
    isread :: Bool
end
AddressSlot(address = 0) = AddressSlot(address, 0x00, false, false)

iswritten(a::AddressSlot) = a.iswritten
isread(a::AddressSlot) = a.isread

setvalue(a::AddressSlot, value) = (a.value = value; a.iswritten = true; nothing)
getvalue(a::AddressSlot) = a.value

getaddress(a::AddressSlot) = a.address
Base.read(a::AddressSlot) = (a.isread = true; a.value)

"""
    validate(jtag::JTAG, ntests; kwargs...)

Test the DRAM for the device referenced by `jtag`. Preform approximately
`ntests` number of reads and writes.

Key Word Agruments
------------------
* `writesize_max = 2^14` - The maximum number of bytes to write in any single
    write action.
* `writesize_min = 2` - The minimum number of bytes to write in any single write
    action.
* `readsize_max = 2^14` - The maximum number of bytes to read in any single read
    action.
* `readsize_min = 1` - The minimum number of bytes to read in any single read
    action.
* `start_address = 0` - The first memory mapped address of the DRAM. If the
    SDRAM is mapped to another address, change this.
    * `end_address = 0x3ff_ffff` - The final memory mapped address of the DRAM.
"""
function validate(
        jtag :: JTAG,
        ntests :: Integer;
        # kwargs!
        writesize_max = 2 ^ 14,
        writesize_min = 2,
        readsize_max = 2^14,
        readsize_min = 1,
        start_address = dram_addr_start(),
        end_address = dram_addr_end()
    )

    # Keep an array tracking all writes that have occurred.
    dram = [AddressSlot(i) for i in start_address:end_address]

    # Track number of reads and writes.
    write_actions = 0
    read_actions = 0
    failed_read_actions = 0

    nwrites = 0
    nreads = 0

    # Set up progress bar
    p = Progress(ntests, 1)

    # Get address range.
    # Make sure to only take even addresses because the Intel SDRAM controller
    # doesn't have low-byte/high-byte enables, so we must always write 2 bytes
    # to an even address.
    address_range = start_address:2:end_address

    while write_actions < ntests || read_actions < ntests
        # Perform a random number of writes of random size.
        for _ in 1:rand(1:5)
            # Pick a random starting address
            base = rand(address_range)

            # Pick a length for the write
            writesize = rand(writesize_min:2:writesize_max)

            # Perform a quick length check to avoid going off the end of the
            # dram.
            if base + writesize > end_address
                writesize = end_address - base
            end

            # Generate the writes - first locally into the local copy of
            # DRAM, then to the actual jtag device.
            for addr in base:base + writesize - 1
                setvalue(dram[addr], rand(UInt8))
            end

            # Create a vector from the values we just wrote and write to jtag.
            data = [getvalue(dram[addr]) for addr in base:base + writesize - 1]
            write(jtag, base, data)

            # Record statistics
            write_actions += 1
            nwrites += writesize
        end

        # Perform a random number of reads
        for _ in 1:rand(1:5)
            # Get a random starting address from addresses written.
            base = getaddress(rand(filter(iswritten, dram)))

            # Inspect DRAM, make sure we only read memory that has been written
            # to.
            local readsize
            for outer readsize in readsize_min:rand(readsize_min:readsize_max)
                if !iswritten(dram[base + readsize - 1])
                    readsize -= 1
                    break
                end
            end

            # Perform a read
            data = read(jtag, base, readsize)
            expected = [read(dram[addr]) for addr in base:base + readsize -1]

            if data != expected
                @error """
                Invalid read of size $readsize starting at address $base.
                """
                failed_read_actions += 1
            end

            # Record stats
            read_actions += 1
            nreads += readsize
            next!(p)
        end
    end

    # Compute the precent of addresses read and written.
    addresses_written = count(iswritten, dram)
    addresses_read = count(isread, dram)

    @info """
    Write actions: $write_actions
    Number of writes: $nwrites
    Unique Addresses Written: $addresses_written
    Write Coverage: $(addresses_written / length(dram))

    Read actions: $read_actions
    Number of reads: $nreads
    Unique Addresses Read: $addresses_read
    Read Coverage: $(addresses_read / length(dram))
    """

    if failed_read_actions == 0
        @info "All operations successful"
    else
        @error "There were $(failed_read_actions) failed read operations"
    end
end

end


text:
#=
This file takes the code from each file in the examples folder and
creates a markdown file of the same name as the julia file then
puts each example from the julia file into a code block and adds
a short html div below with the interactive output.
=#
using PlotlyJS

# used in examples
using Distributions, HTTP, DataFrames, RDatasets, Colors, CSV, JSON
using Random, Dates, LinearAlgebra, DelimitedFiles

const THIS_DIR = dirname(@__FILE__)

# Walk through each example in a file and get the markdown from `single_example`
function single_example_file(filename::String)
    base_fn = split(filename, ".")[1]
    start_example = "```@example $(base_fn)"
    end_example = "```"
    # Open a file to write to
    open(joinpath(THIS_DIR, "src", "examples", "$(base_fn).md"), "w") do outfile

        write_example(ex) = println(outfile, start_example, "\n", ex, "\n", end_example, "\n")

        fn_h1 = titlecase(replace(base_fn, "_" => " "))
        println(outfile, "# $(fn_h1)\n")

        # Read lines from a files
        fulltext = open(
            f->read(f, String),
            joinpath(THIS_DIR, "..", "examples", filename),
            "r"
        )
        all_lines = split(fulltext, "\n")
        l = 1
        regex = r"^function ([^_].+?)\("
        regex_end = r"^end$"

        # find preamble
        if base_fn == "subplots"  # special case
            preamble = "using PlotlyJS, Dates\ninclude(\"../../../examples/line_scatter.jl\")"
            write_example(preamble)
        else
            first_line = findfirst(x -> match(regex, x) !== nothing, all_lines)
            if first_line !== nothing
                preamble = strip(join(all_lines[1:first_line-1], "\n"))
                write_example(preamble)
            end
        end

        while true
            # Find next function name (break if none)
            l = findnext(x -> match(regex, x) !== nothing, all_lines, l+1)
            if l == 0 || l === nothing
                break
            end
            # find corresponding end for this function
            end_l = findnext(x -> match(regex_end, x) !== nothing, all_lines, l+1)

            # Pull out function text
            func_block = join(all_lines[l:end_l], "\n")
            fun_name = match(regex, all_lines[l])[1]
            # println("adding $fun_name")
            an_ex = string(func_block, "\n", fun_name, "()")
            write_example(an_ex)
            l = end_l
        end
    end  # do outfile

    return nothing
end

function main()
    # Read all file names in
    if length(ARGS) == 0
        all_file_names = readdir(joinpath(THIS_DIR, "..", "examples"))
    else
        all_file_names = [endswith(i, ".jl") ? i : "$(i).jl" for i in ARGS]
    end
    all_julia_files = filter(x -> endswith(x, ".jl"), all_file_names)

    foreach(single_example_file, all_julia_files)
end


text:
module UnicdeTest
using GenericRegex
using Base.Test
include("parse/unicode.jl")
end


text:
using Logging, Random
using JuMP, Gurobi
using DecisionProgramming

Random.seed!(13)

const N = 4
const c_k = rand(N)
const b = 0.03
fortification(k, a) = [c_k[k], 0][a]

@info("Creating the influence diagram.")
diagram = InfluenceDiagram()

add_node!(diagram, ChanceNode("L", [], ["high", "low"]))

for i in 1:N
    add_node!(diagram, ChanceNode("R$i", ["L"], ["high", "low"]))
    add_node!(diagram, DecisionNode("A$i", ["R$i"], ["yes", "no"]))
end

add_node!(diagram, ChanceNode("F", ["L", ["A$i" for i in 1:N]...], ["failure", "success"]))

add_node!(diagram, ValueNode("T", ["F", ["A$i" for i in 1:N]...]))

generate_arcs!(diagram)

X_L = [rand(), 0]
X_L[2] = 1.0 - X_L[1]
add_probabilities!(diagram, "L", X_L)

for i in 1:N
    x_R, y_R = rand(2)
    X_R = ProbabilityMatrix(diagram, "R$i")
    X_R["high", "high"] = max(x_R, 1-x_R)
    X_R["high", "low"] = 1 - max(x_R, 1-x_R)
    X_R["low", "low"] = max(y_R, 1-y_R)
    X_R["low", "high"] = 1-max(y_R, 1-y_R)
    add_probabilities!(diagram, "R$i", X_R)
end

X_F = ProbabilityMatrix(diagram, "F")
x_F, y_F = rand(2)
for s in paths([State(2) for i in 1:N])
    denominator = exp(b * sum(fortification(k, a) for (k, a) in enumerate(s)))
    s1 = [s...]
    X_F[1, s1..., 1] = max(x_F, 1-x_F) / denominator
    X_F[1, s..., 2] = 1.0 - X_F[1, s..., 1]
    X_F[2, s..., 1] = min(y_F, 1-y_F) / denominator
    X_F[2, s..., 2] = 1.0 - X_F[2, s..., 1]
end
add_probabilities!(diagram, "F", X_F)

Y_T = UtilityMatrix(diagram, "T")
for s in paths([State(2) for i in 1:N])
    cost = sum(-fortification(k, a) for (k, a) in enumerate(s))
    Y_T[1, s...] = 0 + cost
    Y_T[2, s...] = 100 + cost
end
add_utilities!(diagram, "T", Y_T)

generate_diagram!(diagram, positive_path_utility=true)

model = Model()
z = DecisionVariables(model, diagram)
x_s = PathCompatibilityVariables(model, diagram, z, probability_cut = false)
EV = expected_value(model, diagram, x_s)
@objective(model, Max, EV)

@info("Starting the optimization process.")
optimizer = optimizer_with_attributes(
    () -> Gurobi.Optimizer(Gurobi.Env()),
    "IntFeasTol"      => 1e-9,
)
set_optimizer(model, optimizer)
optimize!(model)

@info("Extracting results.")
Z = DecisionStrategy(z)
U_distribution = UtilityDistribution(diagram, Z)
S_probabilities = StateProbabilities(diagram, Z)

@info("Printing decision strategy:")
print_decision_strategy(diagram, Z, S_probabilities)

@info("State probabilities:")
print_state_probabilities(diagram, S_probabilities, ["L"])
print_state_probabilities(diagram, S_probabilities, [["R$i" for i in 1:N]...])
print_state_probabilities(diagram, S_probabilities, [["A$i" for i in 1:N]...])
print_state_probabilities(diagram, S_probabilities, ["F"])

@info("Printing utility distribution.")
print_utility_distribution(U_distribution)

@info("Printing statistics")
print_statistics(U_distribution)


text:
# The entry format is (cudaname, julianame, kernelcode)
# With single name entries cudaname=julianame and kernelcode=name(xi,yi).
# I commented out functions if I don't know the Julia equivalent.

binary_ops = [
    ("add","+","xi+yi"),
    ("sub","-","xi-yi"),
    ("mul","*","xi*yi"),
    ("div","/","xi/yi"),
    ("pow","^","pow(xi,yi)"),
    ("max","max","(xi>yi?xi:yi)"),
    ("min","min","(xi<yi?xi:yi)"),
    ("eq","==","xi==yi"),
    ("ne","!=","xi!=yi"),
    ("gt",">","xi>yi"),
    ("ge",">=","xi>=yi"),
    ("lt","<","xi<yi"),
    ("le","<=","xi<=yi"),
    # "hypot",
    # "rhypot",
    # "atan2",
    # "frexp",
    # "ldexp",
    # "scalbn",
    # "scalbln",
    # "fmod",
    # "remainder",
    # "mod",
    # "fdim",
    ("invxback","invxback","(-xi*yi*yi)"),
    ("reluback","reluback","(yi>0?xi:0)"),
    ("eluback","eluback","(yi>0?xi:xi*(1+yi))"),
    ("seluback","seluback","(yi>0?1.0507009873554805*xi:xi*(1.7580993408473773+yi))"),
    ("sigmback","sigmback","(xi*yi*(1-yi))"),
    ("tanhback","tanhback","(xi*(1-yi*yi))"),
    # ("rpow","rpow","pow(yi,xi)"),   # need this for Array.^Scalar -> cuda bug #108 switching to CUDA.jl for pow
]

unary_ops_with_int_degree = [
    # cuda does not define negative degrees for bessel, fix it here:
    ("jn", "besselj", "(d>=0 ? jn(d,xi) : d%2==0 ? jn(-d,xi) : -jn(-d,xi))"),
    ("yn", "bessely", "(d>=0 ? yn(d,xi) : d%2==0 ? yn(-d,xi) : -yn(-d,xi))"), # bessely not defined for negative x!
]

# The following list comes from the NVIDIA math docs with some extras.
# http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#mathematical-functions-appendix
# The entry format is (cudaname, julianame, kernelcode)
# With single name entries cudaname=julianame and kernelcode=name(xi).
# I commented out functions if I don't know the Julia equivalent.
unary_ops = [
("abs2", "abs2", "(xi*xi)"),
("abs", "abs", "(xi<0?-xi:xi)"),
"acos",
"acosh",
"asin",
"asinh",
"atan",
"atanh",
"cbrt",
"ceil",
"cos",
"cosh",
"cospi",
# ("cyl_bessel_i0", "besseli0"),  # besseli0,i1 is not defined in SpecialFunctions
# ("cyl_bessel_i1", "besseli1"),
"erf",
"erfc",
"erfcinv",
"erfcx",
"erfinv",
"exp",
"exp10",
"exp2",
"expm1",
"floor",
# "ilogb",
("invx", "invx", "1/xi"),
("j0", "besselj0"),
("j1", "besselj1"),
("gamma_impl", "gamma"),
("lgamma", "lgamma"), # TODO: SpecialFunctions 0.8: lgamma(x::Real)` is deprecated, use `(logabsgamma(x))[1]` instead. Other alternative is loggamma, throws a DomainError if gamma(x) is negative.
("digamma_impl", "digamma"),
("trigamma_impl", "trigamma"),
#
# "llrint",
# "llround",
"log",
"log10",
"log1p",
"log2",
# "logb",
# "lrint",
# "lround",
# "nearbyint",
("neg", "-", "-xi"),
# "normcdf",
# "normcdfinv",
("one", "one", "1"),
# "rcbrt",
("relu", "relu", "(xi>0?xi:0)"),
("elu", "elu", "(xi>0?xi:exp(xi)-1)"),
("selu", "selu", "1.0507009873554805*(xi>0?xi:1.6732632423543778*(exp(xi)-1))"),
# "rint",
"round",
# "rsqrt",
("sigm", "sigm", "(xi>=0?1/(1+exp(-xi)):(exp(xi)/(1+exp(xi))))"),
("sign", "sign", "(xi>0?1:xi<0?-1:0)"),
"sin",
"sinh",
"sinpi",
"sqrt",
"tan",
"tanh",
# "tgamma",
"trunc",
("y0","bessely0"),
("y1","bessely1"),
("zero", "zero", "0"),
]

reduction_ops = [
# The entry format is (cudaname, julianame, merge, item, init)
# ai is the accumulator, xi is the array element
("sum","sum","ai+xi","xi","0"),
("prod","prod","ai*xi","xi","1"),
("maximum","maximum","(ai>xi ? ai : xi)","xi","(-INFINITY)"),
("minimum","minimum","(ai<xi ? ai : xi)","xi","INFINITY"),
("sumabs","sumabs","ai+xi","(xi<0 ? -xi : xi)","0"),
("sumabs2","sumabs2","ai+xi","(xi*xi)","0"),
("maxabs","maxabs","(ai>xi ? ai : xi)","(xi<0 ? -xi : xi)","0"),
("minabs","minabs","(ai<xi ? ai : xi)","(xi<0 ? -xi : xi)","INFINITY"),
("countnz","countnz","ai+xi","(xi!=0)","0"),
]


text:
include("gausslegendre.jl")
include("getdigit_n_xc.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n_xc
using .GaussLegendre
using .MyPrintf

const DIGIT = 17
const MAXN = 1000

function gauss_gl(xc::BigFloat, n::Int64)
    res = GaussLegendre.gauleg(x::BigFloat -> exp(- x * x), BigFloat(0), xc, n)

    return BigFloat(4) * res * res
end

function main()
    xc = GetDigit_n_xc.getxc(gauss_gl, MAXN, DIGIT, true)
    n = GetDigit_n_xc.getn_linear(gauss_gl, xc, MAXN, DIGIT, true)
    precision = GetDigit_n_xc.getdigit(gauss_gl, xc, n, DIGIT, true)
    precision = GetDigit_n_xc.getdigit_2(gauss_gl, xc, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, gauss_gl(xc, n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("xc = %.14f, n = %d, 計算%d桁で、%d桁求まりました\n", xc, n, precision, DIGIT)
end

main()


text:
### Imports ###

import Base: size, setindex!, getindex, broadcast!, Broadcast.fill!, Broadcast.fill

### Exports ###

export ZeroArray, RingBuffer, AnyNode
export fasttanh, canthread

### Types ###

# Useful type for writing/reading inactive containers efficiently
mutable struct ZeroArray{T, N} <: DenseArray{T, N}
  size::NTuple{N, Int}
end
ZeroArray(T::Type, size::NTuple{N,Int}) where N = ZeroArray{T,N}(size)
ZeroArray(arr::DenseArray{T,N}) where {T,N} = ZeroArray{T,N}(size(arr))

mutable struct RingBuffer{T}
  buf::Matrix{T}
  pos::Int
end
RingBuffer(T, width, length) = RingBuffer(zeros(T, width, length), 1)

# TODO: LockedArray for locking individual arrays

mutable struct AnyNode{T} <: AbstractNode
  item::T
end

### Methods ###

Base.size(za::ZeroArray) = za.size
# TODO: Check for OOB situation?
Base.setindex!(za::ZeroArray{T,N}, value, idx...) where {T,N} = zero(T)
Base.getindex(za::ZeroArray{T,N}, idx...) where {T,N} = zero(T)
# TODO: More broadcast overrides
Base.broadcast!(f, za::ZeroArray{T,N}, args...) where {T,N} = za
Base.Broadcast.fill!(za::ZeroArray{T,N}, value) where {T,N} = za
Base.Broadcast.fill(za::ZeroArray{T,N}, value::T) where {T,N} =
  (value == zero(T) ? za : error("Cannot fill a ZeroArray with a non-zero element"))

function rotate!(rb::RingBuffer)
  if rb.pos < size(rb.buf, 2)
    rb.pos += 1
  else
    rb.pos = 1
  end
end
function _rbindex(rb::RingBuffer, idx)
  @assert 0 <= idx <= size(rb.buf, 2)-1 "Invalid RingBuffer index: $idx"
  npos = rb.pos - idx
  if npos < 1
    npos += size(rb.buf, 2)
  end
  npos
end
Base.getindex(rb::RingBuffer, idx::Int) = getindex(rb.buf, Colon(), _rbindex(rb, idx))
Base.getindex(rb::RingBuffer, ::Colon) = getindex(rb.buf, Colon(), rb.pos)
function Base.getindex(rb::RingBuffer{T}, idx::Vector{Int}) where T
  @assert length(idx) == size(rb.buf, 1) "Invalid RingBuffer indexing vector size: $(length(idx))"
  result = zeros(T, size(rb.buf, 1))
  for i in axes(rb.buf, 1)
    result[i] = rb.buf[i, _rbindex(rb, idx[i])]
  end
  result
end
Base.setindex!(rb::RingBuffer, value, idx::Int) = (rb.buf[:, _rbindex(rb, idx)] .= value)
Base.setindex!(rb::RingBuffer, value, ::Colon) = (rb.buf[:, rb.pos] .= value)
function Base.setindex!(rb::RingBuffer{T}, value::Vector{T}, idx::Vector{Int}) where T
  @assert length(idx) == size(rb.buf, 1) "Invalid RingBuffer indexing vector size: $(length(idx))"
  for i in axes(rb.buf, 1)
    rb.buf[i, _rbindex(rb, idx[i])] = value[i]
  end
end
Base.lastindex(rb::RingBuffer) = size(rb.buf, 2)
clear!(rb::RingBuffer{T}) where T = fill!(rb.buf, zero(T))

# Returns a resized array, with as many existing entries copied as possible
function resize_array(arr::Array{T,N} where {T,N}, new_size::NTuple{N,Int} where N, fill_value=zero(T))
  old_size = size(arr)
  new_arr = fill(fill_value, new_size)
  shared_region = CartesianRange(min.(old_size, new_size))
  copy!(new_arr, shared_region, arr, shared_region)
end

### Miscellaneous Utilities ###

fasttanh(x) = x * ( 27 + x * x ) / ( 27 + 9 * x * x )

canthread() = SAMANTHA_THREADS && !Threads.in_threaded_loop.x && Threads.nthreads() > 1


text:
# MTS Algorithms
# MTS: Multiple Trajectory Search for Large Scale Global Optimization,
# By [Lin-Yu Tseng and Chun Chen, 2008](https://sci2s.ugr.es/sites/default/files/files/TematicWebSites/EAMHCO/contributionsCEC08/tseng08mts.pdf)

# MTS Implementation

# Algs
struct MTSAlg <: AbstractOptimizer end
struct LS1Alg <: AbstractOptimizer end

# Options
@with_kw struct MTSOptions
    M = 100
    maxiter=200
    search_range_tol=1e-15
    n_foreground = 80
    n_local_search = 200
    n_local_search_test = 3
    n_local_search_best = 300
    BONUS1 = 10
    BONUS2 = 1
    a_min = 0.4
    a_max = 0.5
    b_min = 0.1
    b_max = 0.3
    c_min = 0
    c_max = 1
    # Fixed parameters
    REDUCE_SEARCH_RANGE_FACTOR = 2.5
    SEARCH_RANGE_DEGERATE_FACTOR = 2
    Y1_INCR = 0.1
    Y1_DECR = 0.1
    X2_INCR = 0.2
end

@with_kw struct LS1Options
    M = 100
    maxiter=200
    search_range_tol=1e-15
    # Fixed parameters
    REDUCE_SEARCH_RANGE_FACTOR = 2.5
    SEARCH_RANGE_DEGERATE_FACTOR = 2
    Y1_INCR = 0.1
    Y1_DECR = 0.1
    X2_INCR = 0.2
    # Dummy parameters
    BONUS1 = 10
    BONUS2 = 1
end

# Workspaces
@params mutable struct MTSWorkspace <: Workspace
    model::VecModel
    x0::AbstractVector
    x::AbstractVector
    options::MTSOptions
    enable::BitVector
    improve::BitVector
    search_range::AbstractVector
    # Volitale variables
    optimal_x::AbstractVector
    optimal_ind::Int
    optimal_val::Real
end

@params mutable struct LS1Workspace <: Workspace
    model::VecModel
    x0::AbstractVector
    x::AbstractVector
    options::LS1Options
    enable::BitVector
    improve::BitVector
    search_range::AbstractVector
    # Volitale variables
    optimal_x::AbstractVector
    optimal_ind::Int
    optimal_val::Real
end

# Workspace constructors
function MTSWorkspace(model::VecModel, x0::AbstractVector, options::MTSOptions; kwargs...)
    @unpack box_min, box_max = model
    M = options.M
    # Initialize improve and serch range
    enable = trues(M)
    improve =  trues(M)
    search_range = [(box_max-box_min) ./ 2 for _ in 1:M]
    MTSWorkspace(model, x0, copy(x0), options, enable, improve, search_range, x0[1], -1, Inf)
end

function LS1Workspace(model::VecModel, x0::AbstractVector, options::LS1Options; kwargs...)
    @unpack box_min, box_max = model
    M = options.M
    # Initialize improve and serch range
    enable = trues(M)
    improve =  trues(M)
    search_range = [(box_max-box_min) ./ 2 for _ in 1:M]
    LS1Workspace(model, x0, copy(x0), options, enable, improve, search_range, x0[1], -1, Inf)
end

# Exposed workspace constructors
function Workspace(model::VecModel, optimizer::LS1Alg, x0::AbstractVector; options::LS1Options=LS1Options(), kwargs...,)
    @assert length(x0) > 0 && x0[1] isa AbstractVector
    if length(model.ineq_constraints) > 0 || length(model.eq_constraints) > 0
        @warn "LS1 does not support (in)equality constraints. Your input would be ignored. "
    end
    return LS1Workspace(model, x0, options)
end

# LS1 Workspace constructor without x0 (use method in paper to initialize)
function Workspace(model::VecModel, optimizer::LS1Alg; options::LS1Options=LS1Options(), kwargs...)
    x0 = initialize_x(model, options)
    return Workspace(model, optimizer, x0; options=options)
end

@params struct MTSResult <: AbstractResult
    minimum
    minimizer
end

# Tool functions
function initialize_x(model::VecModel, options::Union{MTSOptions, LS1Options})
    @unpack box_min, box_max = model
    @unpack M = options
    n_vars = getdim(model)[2]
    SOA = build_SOA(M, n_vars, M)
    x0 = Array{Real,2}(undef, M, n_vars)
    for i in 1:M
        for j in 1:n_vars
            # To be confirmed: At the 5th line of "Multi Trajectory Search" algorithm in paper, I do think (u_i, l_i) shoule be (u_j, l_j)
            x0[i, j] = box_min[j] + (box_max[j]-box_min[j])*(SOA[i, j]/(M-1))
        end
    end
    [x0[i, :] for i in 1:size(x0, 1)]
end

function reduce_search_range(search_range, k, n_vars, box_min, box_max, search_range_tol, REDUCE_SEARCH_RANGE_FACTOR)
    search_range[k] ./= 2
    for i in 1:n_vars
        if search_range[k][i] < search_range_tol
            search_range[k][i] = (box_max[i] - box_min[i]) / REDUCE_SEARCH_RANGE_FACTOR
        end
    end
end

function get_buffered_clamp_and_evaluate!()
    buffer = Dict()
    function buffered_clamp_and_evaluate!(model::AbstractModel, x::AbstractVector)
        if !haskey(buffer, (model, x))
            buffer[((model, x))] = clamp_and_evaluate!(model::AbstractModel, x::AbstractVector)
        end
        return buffer[((model, x))]
    end
    return buffered_clamp_and_evaluate!
end

# Subalgorithms
function _localsearch1(workspace, k)
    @unpack model, options = workspace
    @unpack x, improve, search_range = workspace
    @unpack BONUS1, BONUS2, search_range_tol = options
    @unpack SEARCH_RANGE_DEGERATE_FACTOR, REDUCE_SEARCH_RANGE_FACTOR = options
    @unpack box_min, box_max = model
    n_vars = getdim(model)[2]
    grade = 0
    # search_range in paper is one-dimensional. Expand it to multidimensional.
    if improve[k] == false
        reduce_search_range(search_range, k, n_vars, box_min, box_max, search_range_tol, REDUCE_SEARCH_RANGE_FACTOR)
    end
    improve[k] = false
    buffered_clamp_and_evaluate! = get_buffered_clamp_and_evaluate!()
    for i in 1:n_vars
        # Original value
        _xk = copy(x[k])
        xk_val = buffered_clamp_and_evaluate!(model, _xk)
        update_xki = true
        _xk[i] -= search_range[k][i]
        # Value after update
        _xk_val = buffered_clamp_and_evaluate!(model, _xk)
        # Better than current best solution
        if _xk_val < workspace.optimal_val
            grade += BONUS1
            workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
        end
        # Value stays the same
        if _xk_val == xk_val
            # Restore
            update_xki = false
        else
            # Value degenerates
            if _xk_val > xk_val
                # Restore x_k
                _xk = copy(x[k])
                _xk[i] += search_range[k][i] / SEARCH_RANGE_DEGERATE_FACTOR
                # Current value
                _xk_val = buffered_clamp_and_evaluate!(model, _xk)
                if _xk_val < workspace.optimal_val
                    grade += BONUS1
                    workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
                end
                if _xk_val >= xk_val
                    # Restore
                    update_xki = false
                else
                    grade += BONUS2
                    improve[k] = true
                end
            else
                grade += BONUS2
                improve[k] = true
            end
        end
        if update_xki
            x[k][i] = _xk[i]
        end
    end
    return grade
end

function _localsearch2(workspace, k)
    @unpack model, options = workspace
    @unpack x, improve, search_range = workspace
    @unpack BONUS1, BONUS2, search_range_tol = options
    @unpack SEARCH_RANGE_DEGERATE_FACTOR, REDUCE_SEARCH_RANGE_FACTOR = options
    @unpack box_min, box_max = model
    n_vars = getdim(model)[2]
    grade = 0
    # search_range in paper is one-dimensional. Expand it to multidimensional.
    if improve[k] == false
        reduce_search_range(search_range, k, n_vars, box_min, box_max, search_range_tol, REDUCE_SEARCH_RANGE_FACTOR)
    end
    improve[k] = false
    D = zeros(Int, n_vars)
    r = zeros(Int, n_vars)
    buffered_clamp_and_evaluate! = get_buffered_clamp_and_evaluate!()
    for i in 1:n_vars
        # Original value
        xk_val = buffered_clamp_and_evaluate!(model, x[k])
        update_xk = true
        D .= rand.(Ref([-1, 1]))
        r .= rand.(Ref([0, 1, 2, 3]))
        # Value after update
        _xk = copy(x[k])
        _xk .= [r[_i] == 0 ? _xk[_i]-search_range[k][i]*D[_i] : _xk[_i] for _i in 1:length(r)]
        _xk_val = buffered_clamp_and_evaluate!(model, _xk)
        if _xk_val < workspace.optimal_val
            grade += BONUS1
            workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
        end
        # Value stays the same
        if _xk_val == xk_val
            # Restore
            update_xk = false
        else
            # Value degenerates
            if _xk_val > xk_val
                # Restore x_k
                _xk = copy(x[k])
                _xk .= [r[_i] == 0 ? _xk[_i]+(search_range[k][i]*D[_i] / SEARCH_RANGE_DEGERATE_FACTOR) : _xk[_i] for _i in 1:length(r)]
                _xk_val = buffered_clamp_and_evaluate!(model, _xk)
                if _xk_val < workspace.optimal_val
                    grade += BONUS1
                    workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
                end
                if _xk_val >= xk_val
                    update_xk = false
                else
                    grade += BONUS2
                    improve[k] = true
                end
            else
                grade += BONUS2
                improve[k] = true
            end
        end
        if update_xk
            x[k] = _xk
        end
    end
    return grade
end

function _localsearch3(workspace, k)
    @unpack model, options = workspace
    @unpack x = workspace
    @unpack BONUS1, BONUS2 = options
    @unpack Y1_INCR, Y1_DECR, X2_INCR = options
    @unpack a_min, a_max, b_min, b_max, c_min, c_max = options
    n_vars = getdim(model)[2]
    grade = 0
    _xk = copy(x[k])
    update_xk = false
    buffered_clamp_and_evaluate! = get_buffered_clamp_and_evaluate!()
    for i in 1:n_vars
        # Original value
        _xk_val = buffered_clamp_and_evaluate!(model, _xk)
        update_xk = true
        # Heuristic search
        _xk_x1, _xk_y1, _xk_x2  = copy(_xk), copy(_xk), copy(_xk)
        _xk_x1[i] += Y1_INCR
        _xk_y1[i] -= Y1_DECR
        _xk_x2[i] += X2_INCR
        _xk_x1_val, _xk_y1_val, _xk_x2_val = buffered_clamp_and_evaluate!(model, _xk_x1), buffered_clamp_and_evaluate!(model, _xk_y1), buffered_clamp_and_evaluate!(model, _xk_x2)
        if _xk_x1_val < workspace.optimal_val
            grade += BONUS1
            workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
        end
        if _xk_y1_val < workspace.optimal_val
            grade += BONUS1
            workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
        end
        if _xk_x2_val < workspace.optimal_val
            grade += BONUS1
            workspace.optimal_x, workspace.optimal_ind, workspace.optimal_val = copy(_xk), k, _xk_val
        end
        D1, D2, D3 = _xk_val - _xk_x1_val, _xk_val - _xk_y1_val, _xk_val - _xk_x2_val
        grade += ((D1>0) + (D2>0) + (D3>0)) * BONUS2
        if update_xk == false
            update_xk = (D1>0) || (D2>0) || (D3>0)
        end
        a = (rand() * (a_max - a_min)) + a_min
        b = (rand() * (b_max - b_min)) + b_min
        c = (rand() * (c_max - c_max)) + c_min
        _xk[i] += a*(D1 - D2) + b*(D3 - 2*D1) + c
    end
    if update_xk
        # Update x[k]
        x[k] = _xk
    end
    return grade
end

const LOCAL_SEARCH_METHODS = [_localsearch1, _localsearch2, _localsearch3]

function mts(workspace::MTSWorkspace)
    # Multiple Trajectory Search
    @unpack options, enable = workspace
    @unpack M, n_foreground, n_local_search, n_local_search_test, n_local_search_best = options
    grade_x = [-Inf for _ in 1:M]
    for i in 1:M
        if !enable[i]
            continue
        end
        LS_testgrades = [0 for _ in 1:length(LOCAL_SEARCH_METHODS)]
        LS_testgrades[1] += _localsearch1(workspace, i)
        LS_testgrades[2] += _localsearch2(workspace, i)
        LS_testgrades[3] += _localsearch3(workspace, i)
        max_ind = findmax(LS_testgrades)[2]
        if max_ind == 1
            _best_local_search = _localsearch1
        elseif max_ind == 2
            _best_local_search = _localsearch2
        else
            _best_local_search = _localsearch3
        end
        grade_x[i] = 0
        for _ in 1:n_local_search_best
            grade_x[i] += _best_local_search(workspace, i)
        end
    end
    for _ in 1:n_local_search
        _localsearch1(workspace, workspace.optimal_ind)
    end
    enable_ind = reverse(sortperm(grade_x))[begin:n_foreground]
    enable[:] .= false
    enable[enable_ind] .= true
end

function optimize!(workspace::MTSWorkspace)
    options = workspace.options
    for iter in 1:options.maxiter
        #if debugging[] && iter % 50 == 0
        #        println("Iter ", iter, " max iter: ", options.maxiter)
        #end
        mts(workspace)
    end
    MTSResult(workspace.optimal_val, workspace.optimal_x)
end

# Build a SOA (Simulated Orthogonal Array). Refers to section 2 of the paper posted above.
function build_SOA(m, k, q)
    SOA = Array{Int,2}(undef, m, k)
    for c in 1:k
        q_perm = randperm(q) .- 1
        p = 0
        for r in randperm(m)
            p = (p%q)+1
            SOA[r, c] = q_perm[p]
        end
    end
    SOA
end

# mts Workspace constructor with x0
function Workspace(model::VecModel, optimizer::MTSAlg, x0::AbstractVector; options::MTSOptions=MTSOptions(), kwargs...,)
    @assert length(x0) > 0 && x0[1] isa AbstractVector
    if length(model.ineq_constraints) > 0 || length(model.eq_constraints) > 0
        @warn "MTS does not support (in)equality constraints. Your input would be ignored. "
    end
    return MTSWorkspace(model, x0, options)
end

# mts Workspace constructor without x0 (use method in paper to initialize)
function Workspace(model::VecModel, optimizer::MTSAlg; options::MTSOptions=MTSOptions(), kwargs...)
    x0 = initialize_x(model, options)
    return Workspace(model, optimizer, x0; options=options)
end

# Export localsearch1 independently
function localsearch1(workspace::Union{MTSWorkspace, LS1Workspace})
    M = workspace.options.M
    for i in 1:M
        _localsearch1(workspace, i)
    end
end

# Export LS1 independently
function optimize!(workspace::LS1Workspace)
    options = workspace.options
    for iter in 1:options.maxiter
        #if debugging[] && iter % 50 == 0
        #        println("Iter ", iter, " max iter: ", options.maxiter)
        #end
        localsearch1(workspace)
    end
    MTSResult(workspace.optimal_val, workspace.optimal_x)
end


text:
type GDA
    n_class::Integer
    class_mean::Matrix
    class_cov::Matrix
    class_priors::Vector
    class_::Vector
end

function GDA(;
            class_mean=zeros(2,2),
            class_cov=zeros(2,2),
            class_priors=zeros(2),
            n_class = 2,
            class_ = zeros(2))
    return GDA(n_class, class_mean, class_cov, class_priors,class_)
end

function train!(model::GDA, X::Matrix, y::Vector)
    n_feature = size(X,2)
    n_sample = size(X,1)
    model.class_ = sort(unique(y))
    model.n_class = length(unique(y))
    model.class_mean = zeros(n_feature, model.n_class)
    model.class_cov = cov(X)
    model.class_priors = zeros(model.n_class)
    for (i,class) in enumerate(model.class_)
        X_temp = X[vec(y.==class),:]
        model.class_mean[:,i] = mean(X_temp,1)
        model.class_priors[i] = size(X_temp,1)/n_sample
    end
end

function predict(model::GDA,
                 X::Matrix)
    n_sample = size(X,1)
    res = zeros(n_sample)
    for i = 1:n_sample
        res[i] = predict(model, vec(X[i,:]))
    end
    return res
end

function predict(model::GDA,
                 X::Vector)
    temp = zeros(model.n_class)
    for (i, class) in enumerate(model.n_class)
        p = log_pdf(model, X, i)
        temp[i] = p[1]
    end
    res = indmax(exp(temp - log(sum(exp(temp)))))
    @show res
    class = model.class_[res]
    return class
end

function log_pdf(model::GDA,
                 X::Vector, class::Integer)
    t = size(X,1)
    c = -0.5 * (X - model.class_mean[:, class])' * pinv(model.class_cov) * (X- model.class_mean[:, class])
    d = -0.5 * log(det(model.class_cov))
    p = -t/2*log(2*pi)
    pp = log(model.class_priors[class])
    return c + d + p + pp
end

# therer is some problem to be fixed
function test_GDA()
    X_train, X_test, y_train, y_test = make_cla()
    model = GDA()
    train!(model,X_train, y_train)
    predictions = predict(model,X_test)
    print("classification accuracy", accuracy(y_test, predictions))
end


text:
"""    ftransform(f, s::Shape)

Maps a function `f(x, y, z)` to another `g(λ, μ, ν) * J(λ, μ, ν)`, where
`g(λ, μ, ν)` is basically `f(x(λ, μ, ν), y(λ, μ, ν), z(λ, μ, ν))` within the
volume of the shape `s` but under a change of variables to a rectangular
domain, and `J(λ, μ, ν)` is the Jacobian determinant of the transformation. The
limits of the domain are given by `domain(s)`
"""
function ftransform(f, s::Shape)
    function g(x, y, z)
        p = Point(x, y, z)
        r = f(x, y, z)
        return ifelse(p in s, r, zero(r))
    end
    return g
end

"""    ptransform(s::Shape)

Returns a function that maps a point `(λ, μ, ν)` on the domain given by
`domain(s)` to a tuple `(j, x, y, z)`, where `p = (x, y, z)` corresponds to the
cartesian coordinates of a point inside `s`, and `j` is the is the Jacobian
determinant of the transformation `(x, y, z) ↦ (λ, μ, ν)` evaluated on `p`.
"""
function ptransform end

for S in (:Cube, :Cylinder, :Ellipsoid, :EllipticCylinder, :HollowCylinder,
          :Parallelepiped, :RectangularPyramid, :Ring, :Sphere, :SphericalCap,
          :SquarePyramid, :TriangularToroid, :Torus, :TSP)

    T = Symbol(S, :PT)

    @eval begin
        ptransform(s::$S) = $T(s)
        ftransform{F}(f::F, s::$S) = FunctionTransformation(f, $T(s))
    end
end

function (T::FunctionTransformation)(λ, μ, ν)
    j, x, y, z = T.t(λ, μ, ν)
    return j * T.f(x, y, z)
end

@inline function (T::SpherePT)(λ, θ, φ)
    sθ, cθ = sincos(θ)
    sφ, cφ = sincos(φ)
    r = λ * T.s.r
    rsθ = r * sθ
    j = T.s.r * r * rsθ
    x = rsθ * cφ
    y = rsθ * sφ
    z = r * cθ
    return j, x, y, z
end

@inline function (T::SphericalCapPT)(λ, φ, ν)
    sφ, cφ = sincos(φ)
    κ = 1 - ν
    μ = κ * (T.a - κ)
    k = λ * √μ
    j = T.w * λ * μ
    x = k * cφ
    y = k * sφ
    z = T.s.c * ν
    return j, x, y, z
end

@inline function (T::EllipsoidPT)(λ, θ, φ)
    sθ, cθ = sincos(θ)
    sφ, cφ = sincos(φ)
    λsθ = λ * sθ
    j = T.w * λ * λsθ
    x = λsθ * T.s.a * cφ
    y = λsθ * T.s.b * sφ
    z = λ * T.s.c * cθ
    return j, x, y, z
end

@inline function (T::CylinderPT)(λ, φ, ν)
    sφ, cφ = sincos(φ)
    ρ = λ * T.s.r
    j = T.w * ρ
    x = ρ * cφ
    y = ρ * sφ
    z = ν * T.s.c
    return j, x, y, z
end

@inline function (T::HollowCylinderPT)(λ, φ, ν)
    sφ, cφ = sincos(φ)
    ρ = λ * T.a + T.s.r
    j = T.w * ρ
    x = ρ * cφ
    y = ρ * sφ
    z = ν * T.s.c
    return j, x, y, z
end

@inline function (T::TriangularToroidPT){N}(λ, φ, ν::N)
    sφ, cφ = sincos(φ)
    κ = (1 - ν) * N(0.5)
    ρ = λ * κ * T.s.b + T.s.r
    j = T.w * κ * ρ
    x = ρ * cφ
    y = ρ * sφ
    z = ν * T.s.c
    return j, x, y, z
end

@inline function (T::EllipticCylinderPT)(λ, φ, ν)
    sφ, cφ = sincos(φ)
    j = T.w * λ
    x = λ * T.s.a * cφ
    y = λ * T.s.b * sφ
    z = ν * T.s.c
    return j, x, y, z
end

@inline function (T::CubePT)(λ, μ, ν)
    j = T.w
    x = T.s.a * λ
    y = T.s.a * μ
    z = T.s.a * ν
    return j, x, y, z
end

@inline function (T::ParallelepipedPT)(λ, μ, ν)
    j = T.w
    x = T.s.a * λ
    y = T.s.b * μ
    z = T.s.c * ν
    return j, x, y, z
end

@inline function (T::SquarePyramidPT){N}(λ, μ, ν::N)
    κ = (1 - ν) * N(0.5)
    κa = κ * T.s.a
    j = T.w * κ^2
    x = κa * λ
    y = κa * μ
    z = T.s.b * ν
    return j, x, y, z
end

@inline function (T::RectangularPyramidPT){N}(λ, μ, ν::N)
    κ = (1 - ν) * N(0.5)
    j = T.w * κ^2
    x = κ * T.s.a * λ
    y = κ * T.s.b * μ
    z = T.s.c * ν
    return j, x, y, z
end

@inline function (T::TSPPT){N}(λ, μ, ν::N)
    κ = 1 - T.s.r * (1 + ν) * N(0.5)
    κa = κ * T.s.a
    j = T.w * κ^2
    x = κa * λ
    y = κa * μ
    z = T.a * ν
    return j, x, y, z
end

@inline function (T::RingPT)(λ, θ, φ)
    sθ, cθ = sincos(θ)
    sφ, cφ = sincos(φ)
    aλ = T.s.a * λ
    ρ = aλ * cθ + T.s.R
    j = T.s.b * aλ * ρ
    x = ρ * cφ
    y = ρ * sφ
    z = T.s.b * λ * sθ
    return j, x, y, z
end

@inline function (T::TorusPT)(λ, θ, φ)
    sθ, cθ = sincos(θ)
    sφ, cφ = sincos(φ)
    rλ = T.s.r * λ
    ρ = rλ * cθ + T.s.R
    j = T.s.r * rλ * ρ
    x = ρ * cφ
    y = ρ * sφ
    z = rλ * sθ
    return j, x, y, z
end

### Variables domains
domain(::Sphere            ) = (( 0.0,  0.0,  0.0), (1.0,  1π,  2π))
domain(::Ellipsoid         ) = (( 0.0,  0.0,  0.0), (1.0,  1π,  2π))
domain(::Cylinder          ) = (( 0.0,  -1π, -1.0), (1.0,  1π, 1.0))
domain(::HollowCylinder    ) = (( 0.0,  -1π, -1.0), (1.0,  1π, 1.0))
domain(::EllipticCylinder  ) = (( 0.0,  -1π, -1.0), (1.0,  1π, 1.0))
domain(::TriangularToroid  ) = ((-1.0,  -1π, -1.0), (1.0,  1π, 1.0))
domain(::SphericalCap      ) = (( 0.0,  0.0, -1.0), (1.0,  2π, 1.0))
domain(::Cube              ) = ((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0))
domain(::Parallelepiped    ) = ((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0))
domain(::RectangularPyramid) = ((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0))
domain(::SquarePyramid     ) = ((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0))
domain(::TSP               ) = ((-1.0, -1.0, -1.0), (1.0, 1.0, 1.0))
domain(::Ring              ) = (( 0.0,  -1π,  -1π), (1.0,  1π,  1π))
domain(::Torus             ) = (( 0.0,  -1π,  -1π), (1.0,  1π,  1π))


text:
# Tests that the weighting scheme works for binary_classification

@testset "weightsTest.jl" begin
    try
        if isfile( joinpath(ENV["LIGHTGBM_PATH"],"examples/binary_classification/binary.test") )
            binary_test = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/binary_classification/binary.test"), '\t');
            binary_train = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/binary_classification/binary.train"), '\t');
            binary_test_weight = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/binary_classification/binary.test.weight"), '\t')[:,1];
            binary_train_weight = readdlm(joinpath(ENV["LIGHTGBM_PATH"] , "examples/binary_classification/binary.train.weight"), '\t')[:,1];
        else
            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/binary_classification/binary.test");
            work=String(res.body);
            binary_test =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false));

            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/binary_classification/binary.train");
            work=String(res.body);
            binary_train =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false));

            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/binary_classification/binary.test.weight");
            work=String(res.body);
            binary_test_weight =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false))[:,1];

            res = HTTP.get("https://raw.githubusercontent.com/microsoft/LightGBM/v2.3.1/examples/binary_classification/binary.train.weight");
            work=String(res.body);
            binary_train_weight =convert(Matrix,CSV.read(IOBuffer(work),delim='\t',header=false))[:,1];
        end

        X_train = binary_train[:, 2:end]
        y_train = binary_train[:, 1]
        X_test = binary_test[:, 2:end]
        y_test = binary_test[:, 1]

        # Test binary estimator.
        estimator = LightGBM.LGBMBinary(num_iterations = 20,
                                        learning_rate = .1,
                                        early_stopping_round = 1,
                                        feature_fraction = .8,
                                        bagging_fraction = .9,
                                        bagging_freq = 1,
                                        num_leaves = 1000,
                                        metric = ["auc", "binary_logloss"],
                                        is_training_metric = true,
                                        max_bin = 255,
                                        min_sum_hessian_in_leaf = 0.,
                                        min_data_in_leaf = 1);

        # Test fitting.
        LightGBM.fit(estimator, X_train, y_train, weights=binary_train_weight[:,1]);
        LightGBM.fit(estimator, X_train, y_train, (X_test, y_test), weights=binary_train_weight[:,1]);

        @test true
    catch
        @test false
    end
end


text:
# An implementation of the German variant of the Beaufort cipher in Julia. The only difference is that the key and plainext have switched roles

function germanBeaufort()
    #Get plaintext and key
    println("Please enter the text: ")
    text = readline()
    println("Please enter the keyword: ")
    keyword = readline()
    text = uppercase(text)
    keyword = uppercase(keyword)
    key = ""

    # The keyword must be repeated until it matches the length of the text
    for i = firstindex(text):lastindex(text)
        #pos = ((c - 'A') % length(keyword)) + 1
        pos = i
        while pos > length(keyword)
            pos -= length(keyword)
        end
        key = string(key, keyword[pos])
    end

    # Now encrypt. Works the same as Vigenere, except subtracts rather than adding distance.
    index = 1
    for char in text
        distance = key[index] - 'A'
        newchar = char - distance
        if(newchar > 'Z')
            newchar -= 26
        end
        if(newchar < 'A')
            newchar += 26
        end
        index += 1
        print(newchar)
    end
end
germanBeaufort()


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

relaxed_tols = (default_tol_relax = 100,)
insts = OrderedDict()
insts["minimal"] = [
    ((2, 2), nothing, relaxed_tols),
    ]
insts["various"] = [
    ((2, 2), nothing, relaxed_tols),
    ((1, 4), nothing, relaxed_tols),
    ((4, 1), nothing, relaxed_tols),
    ]
return (EntanglementAssisted, insts)


text:
using BioSequences
using FASTX
using Test
using DataFrames
using CSV
#tes Alt parsing and repr update saving writing to fasta
using st16SseqJuliaTools

primers = "/home/gediminas/working/primers/primers.fasta"
db_file = "/home/gediminas/Desktop/tm/GRCh37_latest_rna.fna"

#collect primer  st16SseqJuliaTools
primer_inf = Array{Tuple{String,String},1}()
reader = FASTA.Reader(open(primers, "r"))
record = FASTA.Record()

ctp = 0
while !eof(reader)
    global ctp += 1
    read!(reader, record)
    id = identifier(record)
    seq = sequence(record)
    push!(primer_inf,(id,seq))
    if ctp == 100
        break
    end
end
close(reader)

data = DataFrame(
 ID=Array{String,1}(),
 SEQ=Array{String,1}(),
 K=Array{Int64,1}(),
 K2=Array{Int64,1}(),
 N=Array{Int64,1}(),
 NMBBIND=Array{Int64,1}(),
 TAIL=Array{Float64,1}(),
 E=Array{Float64,1}(),
 KMER=Array{String,1}(),
 KMER_COUNT=Array{Int64,1}(),
 KMER_MSIMAT=Array{Int64,1}(),
 TM_DNA_DNA = Array{Float64,1}(),
 TM_DNA_RNA_HS = Array{Float64,1}(),
 TM_DNA_RNA_LS = Array{Float64,1}(),
 TM_DNA_DNA_KMER = Array{Float64,1}(),
)

ctn = 0
for (k,k2,n) in [(15,3,1),(15,2,1),(15,0,1),(15,3,2),(15,2,2),(15,0,2),(15,3,3),(15,2,3),(15,0,3),(6,3,1),(6,2,1),(6,0,1),(6,3,2),(6,2,2),(6,0,2),(6,3,3),(6,2,3),(6,0,3)]
    kmers_db = KmerData(db_file,k)
    sm = get_significance_model(k,k2,n,kmers_db)
    global ctn += 1
    println("Working with this set", (k,k2,n) )
    for (id,seq) in primer_inf
        parts = []
        push!(parts,id)
        push!(parts,seq)
        push!(parts,k)
        push!(parts,k2)
        push!(parts,n)
        rescal = collect_matches_data(seq,k,k2,n,kmers_db,sm)
        for r in rescal
            push!(parts,r)
        end
        push!(data,parts)
    end
    CSV.write("REZ_$ctn.csv",data)
end
println(data)


text:
format_iteration(ndigits::Integer, ralign::Bool=true) = generate_formatter("%$(ralign ? "" : "-")$(ndigits)d")

format_percentage(precision::Integer=2, hundredalign::Bool=true, ralign::Bool=true) =
  generate_formatter("%$(ralign ? "" : "-")$(hundredalign ? 4+precision : 3+precision).$(precision)f")


text:
module DateUtils

using Missings, Compat.Dates

export yyyymmdd2yyyy, yyyymmdd2yyyymm, yyyymmdd2date, yyyymm2date, dduuuyyyy2date, uuuddyyyy2date,
  datetimestring2datetime, date2yyyymmdd, date2yyyymm, date2yyyy, datetimestring2yyyymmdd,
  date2yyyyqq, yearqtr2date, date2yyyyww, yearmonth2date

# package code goes here
include("utils.jl")

end # module


text:
Pkg.add("IJulia")
Pkg.add("PyCall")
Pkg.add("PyPlot")

max_stars_repo_name: jessewei/docker-tensorflow-watson

max_issues_repo_name: jessewei/docker-tensorflow-watson

max_forks_repo_name: jessewei/docker-tensorflow-watson


text:
f(x::Int64, y::Float64) = 5x + y

sig = Signature(first(methods(f)))
sig_f = Signature(f)

@testset "Signatures" begin
    @test argnames(sig) == ["x", "y"]
    @test argtypes(sig) == ["Int64", "Float64"]
    @test argnames(sig_f) == argnames(sig)
    @test argtypes(sig_f) == argtypes(sig)
end


text:
using ArgParse

include("convert.jl")

function main()
  s = ArgParseSettings()
  @add_arg_table s begin
    "--categoric", "-c"
      help = "optional: keep categoric features as strings."
      action = :store_true
    "--missing", "-m"
      help = "optional: keep rows with missing values."
      action = :store_true
  end
  parsed_args = parse_args(s)
  keepcat = parsed_args["categoric"]
  keepna = parsed_args["missing"]

  processalldirs("classification", keepcat, keepna)
  processalldirs("regression", keepcat, keepna)
end

main()

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
module Static

import IfElse: ifelse
using Base: @propagate_inbounds, Slice, AbstractCartesianIndex

export StaticInt, StaticFloat64, StaticSymbol, True, False, StaticBool, NDIndex
export dynamic, is_static, known, static

@static if VERSION >= v"1.7.0-DEV.421"
    using Base: @aggressive_constprop
else
    macro aggressive_constprop(ex)
        ex
    end
end

include("int.jl")
include("bool.jl")
include("float.jl")
include("symbol.jl")
include("ndindex.jl")
include("tuples.jl")

"""
    known(::Type{T})

Returns the known value corresponding to a static type `T`. If `T` is not a static type then
`nothing` is returned.

See also: [`static`](@ref), [`is_static`](@ref)
"""
known
@aggressive_constprop known(x) = known(typeof(x))
known(::Type{T}) where {T} = nothing
known(::Type{StaticInt{N}}) where {N} = N::Int
known(::Type{StaticFloat64{N}}) where {N} = N::Float64
known(::Type{StaticSymbol{S}}) where {S} = S::Symbol
known(::Type{Val{V}}) where {V} = V
known(::Type{True}) = true
known(::Type{False}) = false
known(::Type{NDIndex{N,I}}) where {N,I} = known(I)
_get_known(::Type{T}, dim::StaticInt{D}) where {T,D} = known(_get_tuple(T, dim))
function known(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}}
    return eachop(_get_known, nstatic(Val(N)), T)
end

"""
    static(x)

Returns a static form of `x`. If `x` is already in a static form then `x` is returned. If
there is no static alternative for `x` then an error is thrown.

See also: [`is_static`](@ref), [`known`](@ref)

```julia
julia> using Static

julia> static(1)
static(1)

julia> static(true)
True()

julia> static(:x)
static(:x)

```
"""
static
@aggressive_constprop static(x::X) where {X} = ifelse(is_static(X), identity, _no_static_type)(x)
@aggressive_constprop static(x::Int) = StaticInt(x)
@aggressive_constprop static(x::Union{Int8,UInt8,Int16,UInt16}) = StaticInt(x % Int)
@static if sizeof(Int) == 8
    @aggressive_constprop static(x::Union{Int32,UInt32}) = StaticInt(x % Int)
    @aggressive_constprop static(x::UInt64) = StaticInt(Int(x))
else
    @aggressive_constprop static(x::UInt32) = StaticInt(Int(x))
end
@aggressive_constprop static(x::Float64) = StaticFloat64(x)
@aggressive_constprop static(x::Bool) = StaticBool(x)
@aggressive_constprop static(x::Symbol) = StaticSymbol(x)
@aggressive_constprop static(x::Tuple{Vararg{Any}}) = map(static, x)
@generated static(::Val{V}) where {V} = static(V)
function _no_static_type(@nospecialize(x))
    error("There is no static alternative for type $(typeof(x)).")
end
static(x::CartesianIndex) = NDIndex(static(Tuple(x)))

"""
    is_static(::Type{T}) -> StaticBool

Returns `True` if `T` is a static type.

See also: [`static`](@ref), [`known`](@ref)
"""
is_static(@nospecialize(x)) = is_static(typeof(x))
is_static(@nospecialize(x::Type{<:StaticInt})) = True()
is_static(@nospecialize(x::Type{<:StaticBool})) = True()
is_static(@nospecialize(x::Type{<:StaticSymbol})) = True()
is_static(@nospecialize(x::Type{<:Val})) = True()
is_static(@nospecialize(x::Type{<:StaticFloat64})) = True()
is_static(x::Type{T}) where {T} = False()

@aggressive_constprop _tuple_static(::Type{T}, i) where {T} = is_static(_get_tuple(T, i))
function is_static(::Type{T}) where {N,T<:Tuple{Vararg{Any,N}}}
    if all(eachop(_tuple_static, nstatic(Val(N)), T))
        return True()
    else
        return False()
    end
end

"""
    dynamic(x)

Returns the "dynamic" or non-static form of `x`.
"""
dynamic(x::X) where {X} = _dynamic(is_static(X), x)
_dynamic(::True, x::X) where {X} = known(X)
_dynamic(::False, x::X) where {X} = x
@aggressive_constprop dynamic(x::Tuple) = map(dynamic, x)
dynamic(x::NDIndex) = CartesianIndex(dynamic(Tuple(x)))

end


text:
""" Basic tests for the `Hsq` struct and associated function. """

include("../src/LDScore.jl")

using Main
using Main.LDScore
using Test

eps = 1.0e-6 # accuracy of approximatly equal; the lower the more accurate

approx(x, y, eps) = abs(x - y) <= eps

function vector_approx(x, y, eps)
    if size(x) != size(y)
        return false
    end
    for (i, j) in enumerate(x)
        if !approx(x[i], y[i], eps)
            return false
        end
    end
    return true
end

# Test to ensure that the basic structure of `Hsq` is intact
@testset "Hsq struct" begin
    # column vectors d = 6 * 1
    y = [[0.2061] [0.2601] [4.3514] [6.1703] [5.0221] [2.418]]'
    x = [[6.1388] [6.4785] [8.88] [4.8064] [3.6219] [3.7244]]'
    w = [[6.1388] [6.4785] [8.88] [4.8064] [3.6219] [3.7244]]'
    N = [[592497.0] [591333.0] [696881.0] [617319.0] [687344.0] [685637.0]]'
    M = [[1173569.0]]
    n_blocks = 2
    intercept = nothing
    slow = false
    step1_ii = [[true] [true] [true] [true] [true] [true]]'
    old_weights = false

    hsq = LDScore.Hsq(
        y, x, w, N, M, n_blocks = n_blocks, intercept = intercept, slow = slow,
        step1_ii = step1_ii, old_weights = old_weights,
    )

    hsq.old_weights == false
    LDScore.ld_score_regression(
        hsq, y, x, w, N, M, n_blocks, intercept, slow, step1_ii, old_weights,
    )

    @test typeof(hsq) == LDScore.Hsq
    @test hsq.__null_intercept__ == 1
    @test hsq.y == y
    @test hsq.x == x
    @test hsq.w == w
    @test hsq.N == N
    @test hsq.M == M
    @test hsq.n_blocks == n_blocks
    @test hsq.intercept == intercept
    @test hsq.slow == slow
    @test hsq.step1_ii == step1_ii
end

@testset "Hsq functions" begin
    χ² = ones((4, 1)) .* 4
    ld = ones((4, 1))
    w_ld = ones((4, 1))
    N = 9 .* ones((4, 1))
    M = [[7]]
    hsq = LDScore.Hsq(
        χ², ld, w_ld, N, M, n_blocks = 3, intercept = 1,
        slow = false, old_weights = false, step1_ii = nothing,
    )

    @testset "weights" begin
        x = [[0.09570313] [0.09570313] [0.09570313] [0.09570313]]'
        a = LDScore.hsq_weights(ld, w_ld, N, M, 1)
        b = LDScore.hsq_weights(ld, w_ld, N, M, 2)
        @test vector_approx(a, x, eps)
        @test vector_approx(a, b, eps)

        x = [[0.5] [0.5] [0.5] [0.5]]'
        a = LDScore.hsq_weights(ld, w_ld, N, M, 0)
        b = LDScore.hsq_weights(ld, w_ld, N, M, -1)
        @test vector_approx(a, x, eps)
        @test vector_approx(a, b, eps)
    end

    @testset "aggregate" begin
        χ² = ones((10, 1)) .* 3 ./ 2
        ld = ones((10, 1)) .* 100
        N = ones((10, 1)) .* 100000
        M = 1e7
        agg = LDScore.aggregate(hsq, χ², ld, N, M)
        approx(agg, 0.5, eps)
        agg = LDScore.aggregate(hsq, χ², ld, N, M; intercept = 1.5)
        approx(agg, 0, eps)
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule TauDEM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("TauDEM")
JLLWrappers.@generate_main_file("TauDEM", UUID("f79d145f-20fd-5c1a-ae67-8074679d1f1f"))
end  # module TauDEM_jll


text:
struct Level{Tprice, Tvolume}
    price::Tprice
    volume::Tvolume
end

import Base: isless
isless(l1::Level, l2::Level) = isless(l1.price, l2.price)

function isunique(v_lev::Array{Level{Tprice,Tvolume},1}) where {Tprice,Tvolume}
    length(unique([lv.price for lv in v_lev])) == length(v_lev)
end


text:
@testset "Rosenbrock" begin

    function test_result(result::Evolutionary.EvolutionaryOptimizationResults, N::Int, tol::Float64)
        fitness = minimum(result)
        extremum = Evolutionary.minimizer(result)
        @test sum(abs2,extremum.-ones(N)) ≈ 0.0 atol=tol
        @test fitness ≈ 0.0 atol=tol
    end

    # Objective function
    rosenbrock(x::AbstractVector) = (1.0 - x[1])^2 + 100 * (x[2] - x[1]^2)^2

    # Parameters
    N = 2

    # Testing: (15,3/(+,)100)-ES
    settings = [
        :isotropic=>(gaussian, gaussian, IsotropicStrategy(N)),
        :isotropic=>(cauchy, gaussian, IsotropicStrategy(N)),
        :anisotropic=>(gaussian, gaussian, AnisotropicStrategy(N))
    ]
    selections = [:plus, :comma]
    @testset "ES settings" for (sn,ss) in settings, sel in selections
        m = ES(
            initStrategy = ss[3],
            recombination = average, srecombination = average,
            mutation = ss[1], smutation = ss[2],
            μ = 15, ρ = 3, λ = 100, selection=sel
        )
        result = Evolutionary.optimize( rosenbrock, (() -> rand(N)), m, Evolutionary.Options(iterations=1000, successive_f_tol=25))
        println("$(summary(m)):$(sn) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
        test_result(result, N, sel == :plus ? 0.1 : 0.5)
    end
    m = ES( initStrategy = IsotropicStrategy(N),
                recombination = average, srecombination = average,
                mutation = gaussian, smutation = gaussian,
                μ = 15, ρ = 3, λ = 100)
    result = Evolutionary.optimize(rosenbrock, BoxConstraints(0.0, 0.5, N), (() -> rand(N)), m)
    println("$(summary(m)) [box] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.minimizer(result) ≈ [0.5, 0.25] atol=1e-1

    # Testing: CMA-ES
    result = Evolutionary.optimize(rosenbrock, (() -> rand(Float32,N)), CMAES(mu = 5, lambda = 100, weights=zeros(Float32,100)))
    println("(5/5,100)-CMA-ES => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-2)
    result = Evolutionary.optimize(rosenbrock, fill(0.0, N), fill(0.5, N), (() -> rand(N)), CMAES(mu = 5, lambda = 100))
    println("(5/5,100)-CMA-ES [box] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.minimizer(result) ≈ [0.5, 0.25] atol=1e-5

    con_c!(x) = [sum(x)]
    c = PenaltyConstraints(100.0, fill(0.0, 5N), Float64[], [1.0], [1.0], con_c!)
    result = Evolutionary.optimize(rosenbrock, c, (() -> rand(5N)), CMAES(mu = 40, lambda = 100))
    println("(5/5,100)-CMA-ES [penalty] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.minimizer(result) |> sum ≈ 1.0 atol=1e-1

    c = PenaltyConstraints(100.0, fill(0.0, 2N), fill(0.5, 2N), [1.0], [1.0], con_c!)
    result = Evolutionary.optimize(rosenbrock, c, (() -> rand(2N)), CMAES(mu = 8, lambda = 100))
    println("(5/5,100)-CMA-ES [penalty] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.minimizer(result) |> sum ≈ 1.0 atol=1e-1
    @test all(0.0 <= x+0.01 && x-0.01 <= 0.5 for x in abs.(Evolutionary.minimizer(result)))

    # Testing: GA
    m = GA(
        populationSize = 100,
        ɛ = 0.1,
        selection = rouletteinv,
        crossover = intermediate(0.25),
        mutation = domainrange(fill(0.5,N))
    )
    result = Evolutionary.optimize(rosenbrock, (() -> rand(N)), m)
    println("GA(p=100,x=0.8,μ=0.1,ɛ=0.1) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-1)
    result = Evolutionary.optimize(rosenbrock, BoxConstraints(0.0, 0.5, N), (() -> rand(N)), m)
    println("GA(p=100,x=0.8,μ=0.1,ɛ=0.1)[box] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.minimizer(result) ≈ [0.5, 0.25] atol=1e-1

    # Testing: DE
    result = Evolutionary.optimize(rosenbrock, (() -> rand(N)), DE(populationSize = 100))
    println("DE/rand/1/bin(F=1.0,Cr=0.5) => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    test_result(result, N, 1e-1)
    result = Evolutionary.optimize(rosenbrock, BoxConstraints(0.0, 0.5, N), (() -> rand(Float32,N)), DE(populationSize = 100))
    println("DE/rand/1/bin(F=1.0,Cr=0.5)[box] => F: $(minimum(result)), C: $(Evolutionary.iterations(result))")
    @test Evolutionary.minimizer(result) ≈ [0.5, 0.25] atol=1e-1

end


text:
# Autogenerated wrapper script for flux_core_jll for x86_64-linux-gnu
export flux, libflux_core

using libsodium_jll
using ZeroMQ_jll
using czmq_jll
using Jansson_jll
using Libuuid_jll
using Lz4_jll
using Hwloc_jll
using SQLite_jll
using Lua_jll
JLLWrappers.@generate_wrapper_header("flux_core")
JLLWrappers.@declare_executable_product(flux)
JLLWrappers.@declare_library_product(libflux_core, "libflux-core.so.2")
function __init__()
    JLLWrappers.@generate_init_header(libsodium_jll, ZeroMQ_jll, czmq_jll, Jansson_jll, Libuuid_jll, Lz4_jll, Hwloc_jll, SQLite_jll, Lua_jll)
    JLLWrappers.@init_executable_product(
        flux,
        "bin/flux",
    )

    JLLWrappers.@init_library_product(
        libflux_core,
        "lib/libflux-core.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Autogenerated wrapper script for ViennaRNA_jll for x86_64-w64-mingw32-cxx11
export AnalyseDists, AnalyseSeqs, Kinfold, Kinfold, RNA2Dfold, RNALalifold, RNALfold, RNAPKplex, RNAaliduplex, RNAalifold, RNAcofold, RNAdistance, RNAdos, RNAduplex, RNAeval, RNAfold, RNAforester, RNAheat, RNAinverse, RNAlocmin, RNAmultifold, RNApaln, RNAparconv, RNApdist, RNAplex, RNAplfold, RNAplot, RNApvmin, RNAsnoop, RNAsubopt, RNAup, b2ct, ct2db, libRNA, popt

using MPFR_jll
using CompilerSupportLibraries_jll
using GSL_jll
JLLWrappers.@generate_wrapper_header("ViennaRNA")
JLLWrappers.@declare_executable_product(AnalyseDists)
JLLWrappers.@declare_executable_product(AnalyseSeqs)
JLLWrappers.@declare_executable_product(Kinfold)
JLLWrappers.@declare_executable_product(Kinfold)
JLLWrappers.@declare_executable_product(RNA2Dfold)
JLLWrappers.@declare_executable_product(RNALalifold)
JLLWrappers.@declare_executable_product(RNALfold)
JLLWrappers.@declare_executable_product(RNAPKplex)
JLLWrappers.@declare_executable_product(RNAaliduplex)
JLLWrappers.@declare_executable_product(RNAalifold)
JLLWrappers.@declare_executable_product(RNAcofold)
JLLWrappers.@declare_executable_product(RNAdistance)
JLLWrappers.@declare_executable_product(RNAdos)
JLLWrappers.@declare_executable_product(RNAduplex)
JLLWrappers.@declare_executable_product(RNAeval)
JLLWrappers.@declare_executable_product(RNAfold)
JLLWrappers.@declare_executable_product(RNAforester)
JLLWrappers.@declare_executable_product(RNAheat)
JLLWrappers.@declare_executable_product(RNAinverse)
JLLWrappers.@declare_executable_product(RNAlocmin)
JLLWrappers.@declare_executable_product(RNAmultifold)
JLLWrappers.@declare_executable_product(RNApaln)
JLLWrappers.@declare_executable_product(RNAparconv)
JLLWrappers.@declare_executable_product(RNApdist)
JLLWrappers.@declare_executable_product(RNAplex)
JLLWrappers.@declare_executable_product(RNAplfold)
JLLWrappers.@declare_executable_product(RNAplot)
JLLWrappers.@declare_executable_product(RNApvmin)
JLLWrappers.@declare_executable_product(RNAsnoop)
JLLWrappers.@declare_executable_product(RNAsubopt)
JLLWrappers.@declare_executable_product(RNAup)
JLLWrappers.@declare_executable_product(b2ct)
JLLWrappers.@declare_executable_product(ct2db)
JLLWrappers.@declare_library_product(libRNA, "libRNA.dll")
JLLWrappers.@declare_executable_product(popt)
function __init__()
    JLLWrappers.@generate_init_header(MPFR_jll, CompilerSupportLibraries_jll, GSL_jll)
    JLLWrappers.@init_executable_product(
        AnalyseDists,
        "bin\\AnalyseDists.exe",
    )

    JLLWrappers.@init_executable_product(
        AnalyseSeqs,
        "bin\\AnalyseSeqs.exe",
    )

    JLLWrappers.@init_executable_product(
        Kinfold,
        "bin\\kinwalker.exe",
    )

    JLLWrappers.@init_executable_product(
        Kinfold,
        "bin\\Kinfold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNA2Dfold,
        "bin\\RNA2Dfold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNALalifold,
        "bin\\RNALalifold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNALfold,
        "bin\\RNALfold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAPKplex,
        "bin\\RNAPKplex.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAaliduplex,
        "bin\\RNAaliduplex.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAalifold,
        "bin\\RNAalifold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAcofold,
        "bin\\RNAcofold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAdistance,
        "bin\\RNAdistance.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAdos,
        "bin\\RNAdos.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAduplex,
        "bin\\RNAduplex.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAeval,
        "bin\\RNAeval.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAfold,
        "bin\\RNAfold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAforester,
        "bin\\RNAforester.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAheat,
        "bin\\RNAheat.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAinverse,
        "bin\\RNAinverse.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAlocmin,
        "bin\\RNAlocmin.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAmultifold,
        "bin\\RNAmultifold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNApaln,
        "bin\\RNApaln.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAparconv,
        "bin\\RNAparconv.exe",
    )

    JLLWrappers.@init_executable_product(
        RNApdist,
        "bin\\RNApdist.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAplex,
        "bin\\RNAplex.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAplfold,
        "bin\\RNAplfold.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAplot,
        "bin\\RNAplot.exe",
    )

    JLLWrappers.@init_executable_product(
        RNApvmin,
        "bin\\RNApvmin.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAsnoop,
        "bin\\RNAsnoop.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAsubopt,
        "bin\\RNAsubopt.exe",
    )

    JLLWrappers.@init_executable_product(
        RNAup,
        "bin\\RNAup.exe",
    )

    JLLWrappers.@init_executable_product(
        b2ct,
        "bin\\b2ct.exe",
    )

    JLLWrappers.@init_executable_product(
        ct2db,
        "bin\\ct2db.exe",
    )

    JLLWrappers.@init_library_product(
        libRNA,
        "bin\\libRNA.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        popt,
        "bin\\popt.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Mock2
using Test

@testset "Mock2.jl" begin
  @test  addone(1)==2 ;
  @test  addtwo(2)==4;
  @test  my_f(0)==1;
  @test  my_f_der(1)==-3;
end


text:
using CellularAutomata

const gens = 3

blinker = [[0, 0, 0, 0, 0] [0, 0, 1, 0, 0] [0, 0, 1, 0, 0] [0, 0, 1, 0, 0] [0, 0, 0, 0, 0]]

ca = CellularAutomaton(Life((3, (2,3))), blinker, gens)

@test ca.evolution == cat([0 0 0 0 0; 0 0 0 0 0; 0 1 1 1 0; 0 0 0 0 0; 0 0 0 0 0],
[0 0 0 0 0; 0 0 1 0 0; 0 0 1 0 0; 0 0 1 0 0; 0 0 0 0 0],
[0 0 0 0 0; 0 0 0 0 0; 0 1 1 1 0; 0 0 0 0 0; 0 0 0 0 0], dims=3)


text:
struct Symel
    symbol
    rrep
end

struct Op
    optype
    order
    nprimes
end

struct Chartable
    name
    symels
    irreps
    characters
end

struct PG
    str
    family
    n
    subfamily
end

struct Rep
    irrep::Bool
    d::Int
    p::Int
    v::Int
    h::Int
    i::Int
    l::Int
    name
end

function Base.:(==)(A::Symel, B::Symel)
    if sum(A.rrep .- B.rrep) < tol
        return true
    else
        return false
    end
end


text:
# TODO this should be reimplemented to work for nd to nd stuff like reshaped arrays
struct Reshape{I,O,S} <: IndexTransform{I,O}
    size::Size{S}

    Reshape(s::Size{S}) where {S} = new{1,known_length(S),S}(s)
end

function (r::Reshape{1,O,S})(i::CanonicalInt) where {I,O,S}
    _lin2sub(NOnes(static(O)), getfield(s, :size), i)
end
@generated function _lin2sub(o::O, s::S, i::I) where {O,S,I}
    out = Expr(:block, Expr(:meta, :inline))
    t = Expr(:tuple)
    iprev = :(i - 1)
    N = length(S.parameters)
    for i in 1:N
        if i === N
            push!(t.args, :($iprev + getfield(o, $i)))
        else
            len = gensym()
            inext = gensym()
            push!(out.args, :($len = getfield(s, $i)))
            push!(out.args, :($inext = div($iprev, $len)))
            push!(t.args, :($iprev - $len * $inext + getfield(o, $i)))
            iprev = inext
        end
    end
    push!(out.args, :(NDIndex($(t))))
    out
end