text:
"""
`module  PhotoExcitationAutoion`
    ... a submodel of JAC that contains all methods for computing photo-excitation-autoionization cross sections and rates.
"""
module PhotoExcitationAutoion

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

    """
    `struct  PhotoExcitationAutoion.Settings  <:  AbstractProcessSettings`
        ... defines a type for the details and parameters of computing photon-impact excitation-autoionization pathways
            |i(N)>  --> |m(N)>  --> |f(N-1)>.

        + multipoles              ::Array{Basics.EmMultipole,1} ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{Basics.UseGauge,1}    ... Specifies the gauges to be included into the computations.
        + calcPartialCs           ::Bool                        ... Calculate the partial excitation-autoionization cross sections
                                                                    sigma(i-e-f), and where the direct photoionization is neglected.
        + calcAngular             ::Bool                        ... Calculate the angular distribution of the resonantly emitted Auger electrons
                                                                    at selected solid angles.
        + calcFano                ::Bool                        ... Calculate the Fano parameters of individual cs resonances sigma(i-e-f).
        + printBefore             ::Bool                        ... True, if all energies and lines are printed before their evaluation.
        + incidentStokes          ::ExpStokes                   ... Stokes parameters of the incident radiation.
        + solidAngles             ::Array{SolidAngle,1}         ... List of solid angles [(theta_1, pho_1), ...].
        + electronEnergyShift     ::Float64                     ... An overall energy shift for all electron energies.
        + maxKappa                ::Int64                       ... Maximum kappa value of partial waves to be included for free electrons.
        + pathwaySelection        ::PathwaySelection            ... Specifies the selected levels/pathways, if any.
    """
    struct Settings  <:  AbstractProcessSettings
        multipoles                ::Array{Basics.EmMultipole,1}
        gauges                    ::Array{Basics.UseGauge,1}
        calcPartialCs             ::Bool
        calcAngular               ::Bool
        calcFano                  ::Bool
        printBefore               ::Bool
        incidentStokes            ::ExpStokes
        solidAngles               ::Array{SolidAngle,1}
        electronEnergyShift       ::Float64
        maxKappa                  ::Int64
        pathwaySelection          ::PathwaySelection
    end

    """
    `PhotoExcitationAutoion.Settings()`
        ... constructor for the default values of photon-impact excitation-autoionizaton settings.
    """
    function Settings()
        Settings( Basics.EmMultipole[], UseGauge[], false, false, false, false, Basics.ExpStokes(), SolidAngle[], 0., 0, PathwaySelection())
    end

    # `Base.show(io::IO, settings::PhotoExcitationAutoion.Settings)`
    # 	... prepares a proper printout of the variable settings::PhotoExcitationAutoion.Settings.
    function Base.show(io::IO, settings::PhotoExcitationAutoion.Settings)
        println(io, "multipoles:              $(settings.multipoles)  ")
        println(io, "gauges:                  $(settings.gauges)  ")
        println(io, "calcPartialCs:           $(settings.calcPartialCs)  ")
        println(io, "calcAngular:             $(settings.calcAngular)  ")
        println(io, "calcFano:                $(settings.calcFano)  ")
        println(io, "printBefore:             $(settings.printBefore)  ")
        println(io, "incidentStokes:          $(settings.incidentStokes)  ")
        println(io, "solidAngles:             $(settings.solidAngles)  ")
        println(io, "electronEnergyShift:     $(settings.electronEnergyShift)  ")
        println(io, "maxKappa:                $(settings.maxKappa)  ")
        println(io, "pathwaySelection:        $(settings.pathwaySelection)  ")
    end

    """
    `struct  PhotoExcitationAutoion.Pathway`
        ... defines a type for a photon-impact excitation pathway that may include the definition of different
            excitation and autoionization channels and their corresponding amplitudes.

        + initialLevel        ::Level           ... initial-(state) level
        + intermediateLevel   ::Level           ... intermediate-(state) level
        + finalLevel          ::Level           ... final-(state) level
        + excitEnergy         ::Float64         ... photon excitation energy of this pathway
        + electronEnergy      ::Float64         ... energy of the (finally outgoing, scattered) electron
        + partialCs           ::EmProperty      ... partial cross section sigma(i-e-f) of this pathway
        + qFano               ::EmProperty      ... Fano-q parameter of a resonance (i-e-f)
        + excitChannels       ::Array{PhotoEmission.Channel,1}      ... List of excitation channels of this pathway.
        + augerChannels       ::Array{AutoIonization.Channel,1}     ... List of Auger channels of this pathway.
        + photoChannels       ::Array{PhotoIonization.Channel,1}    ... List of photoionization channels of this pathway.
    """
    struct  Pathway
        initialLevel          ::Level
        intermediateLevel     ::Level
        finalLevel            ::Level
        excitEnergy           ::Float64
        electronEnergy        ::Float64
        partialCs             ::EmProperty
        qFano                 ::EmProperty
        excitChannels         ::Array{PhotoEmission.Channel,1}
        augerChannels         ::Array{AutoIonization.Channel,1}
        photoChannels         ::Array{PhotoIonization.Channel,1}
    end

    """
    `PhotoExcitationAutoion.Pathway()`
        ... 'empty' constructor for an photon-impact excitation-autoionization pathway between a specified initial, intermediate
            and final level.
    """
    function Pathway()
        Pathway(Level(), Level(), Level(), 0., 0., EmProperty(0., 0.), EmProperty(0., 0.), PhotoEmission.Channel[],
                AutoIonization.Channel[], PhotoIonization.Channel[] )
    end

    # `Base.show(io::IO, pathway::PhotoExcitationAutoion.Pathway)`
    #		... prepares a proper printout of the variable pathway::PhotoExcitationAutoion.Pathway.
    function Base.show(io::IO, pathway::PhotoExcitationAutoion.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "excitEnergy                 $(pathway.excitEnergy)  ")
        println(io, "electronEnergy              $(pathway.electronEnergy)  ")
        println(io, "partialCs:                  $(pathway.partialCs)  ")
        println(io, "qFano:                      $(pathway.qFano)  ")
        println(io, "excitChannels:              $(pathway.excitChannels)  ")
        println(io, "augerChannels:              $(pathway.augerChannels)  ")
        println(io, "photoChannels:              $(pathway.photoChannels)  ")
    end

    """
    `PhotoExcitationAutoion.computeAmplitudesProperties(pathway::PhotoExcitationAutoion.Pathway, nm::Nuclear.Model, grid::Radial.Grid,
                                                        nrContinuum::Int64, settings::PhotoExcitationAutoion.Settings)`
        ... to compute all amplitudes and properties of the given pathway; a line::PhotoExcitationAutoion.Pathway is returned for which
            the amplitudes and properties have now been evaluated.
    """
    function  computeAmplitudesProperties(pathway::PhotoExcitationAutoion.Pathway, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                          settings::PhotoExcitationAutoion.Settings)
        # Compute all excitation channels
        neweChannels = PhotoEmission.Channel[]
        for eChannel in pathway.excitChannels
            amplitude   = PhotoEmission.amplitude("absorption", eChannel.multipole, eChannel.gauge, pathway.excitEnergy,
                                                  pathway.intermediateLevel, pathway.initialLevel, grid)
            push!( neweChannels, PhotoEmission.Channel( eChannel.multipole, eChannel.gauge, amplitude))
        end
        # Compute all AutoIonization decay channels
        newaChannels = AutoIonization.Channel[];   contSettings = Continuum.Settings(false, nrContinuum)
        for aChannel in pathway.augerChannels
            newnLevel   = Basics.generateLevelWithSymmetryReducedBasis(pathway.intermediateLevel, pathway.intermediateLevel.basis.subshells)
            newnLevel   = Basics.generateLevelWithExtraSubshell(Subshell(101, aChannel.kappa), newnLevel)
            newfLevel   = Basics.generateLevelWithSymmetryReducedBasis(pathway.finalLevel, pathway.finalLevel.basis.subshells)
            cOrbital, phase  = Continuum.generateOrbitalForLevel(pathway.electronEnergy, Subshell(101, aChannel.kappa), newfLevel, nm, grid, contSettings)
            newcLevel   = Basics.generateLevelWithExtraElectron(cOrbital, aChannel.symmetry, newfLevel)
            newcChannel = AutoIonization.Channel( aChannel.kappa, aChannel.symmetry, phase, Complex(0.))
            amplitude = 1.0
            ## amplitude   = AutoIonization.amplitude("Coulomb", aChannel, newnLevel, newcLevel, grid)
            push!( newaChannels, AutoIonization.Channel( aChannel.kappa, aChannel.symmetry, phase, amplitude))
        end
        # Compute all photoionization channels
        newpChannels = PhotoIonization.Channel[]
        for pChannel in pathway.photoChannels
            ## amplitude   = PhotoIonization.amplitude("absorption", eChannel.multipole, eChannel.gauge, pathway.excitEnergy,
            ##                                         pathway.intermediateLevel, pathway.initialLevel, grid)
            amplitude = 1.0im
            push!( newpChannels, PhotoIonization.Channel( pChannel.multipole, pChannel.gauge, pChannel.kappa, pChannel.symmetry, pChannel.phase, amplitude))
        end
        #
        partialCs = EmProperty(-1., -1.)
        qFano     = EmProperty(-2., -2.)
        pathway = PhotoExcitationAutoion.Pathway( pathway.initialLevel, pathway.intermediateLevel, pathway.finalLevel, pathway.excitEnergy,
                                                  pathway.electronEnergy, partialCs, qFano, neweChannels, newaChannels, newpChannels)
        return( pathway )
    end

    """
    `PhotoExcitationAutoion.computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                            nm::Nuclear.Model, grid::Radial.Grid, settings::PhotoExcitation.Settings; output=true)`
        ... to compute the photo-excitation-autoionization amplitudes and all properties as requested by the given settings. A list of
            lines::Array{PhotoExcitationAutoion.Lines} is returned.
    """
    function  computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet, nm::Nuclear.Model,
                              grid::Radial.Grid, settings::PhotoExcitationAutoion.Settings; output=true)
        println("")
        printstyled("PhotoExcitationAutoion.computePathways(): The computation of photo-excitation-autoionization amplitudes starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        pathways = PhotoExcitationAutoion.determinePathways(finalMultiplet, intermediateMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    PhotoExcitationAutoion.displayPathways(stdout, pathways, settings)    end
        # Determine maximum (electron) energy and check for consistency of the grid
        maxEnergy = 0.;   for  pathway in pathways   maxEnergy = max(maxEnergy, pathway.electronEnergy)   end
        nrContinuum = Continuum.gridConsistency(maxEnergy, grid)
        # Calculate all amplitudes and requested properties
        newPathways = PhotoExcitationAutoion.Pathway[]
        for  pathway in pathways
            push!( newPathways, PhotoExcitationAutoion.computeAmplitudesProperties(pathway, nm, grid, nrContinuum, settings) )
        end
        # Print all results to screen
        if  settings.calcPartialCs      PhotoExcitationAutoion.displayPartialCs(stdout, newPathways)         end
        if  settings.calcAngular        PhotoExcitationAutoion.displayAngular(stdout, newPathways)           end
        if  settings.calcFano           PhotoExcitationAutoion.displayFanoParameters(stdout, newPathways)    end
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary
            if  settings.calcPartialCs      PhotoExcitationAutoion.displayPartialCs(iostream, newPathways)         end
            if  settings.calcAngular        PhotoExcitationAutoion.displayAngular(iostream, newPathways)           end
            if  settings.calcFano           PhotoExcitationAutoion.displayFanoParameters(iostream, newPathways)    end
        end
        #
        if    output    return( newPathways )
        else            return( nothing )
        end
    end

    """
    `PhotoExcitationAutoion.determinePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                              settings::PhotoExcitationAutoion.Settings)`
        ... to determine a list of photoexcitation-autoionization pathways between the levels from the given initial-, intermediate- and
            final-state multiplets and by taking into account the particular selections and settings for this computation; an
            Array{PhotoExcitationAutoion.Line,1} is returned. Apart from the level specification, all physical properties are set to zero
            during the initialization process.
    """
    function  determinePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                settings::PhotoExcitationAutoion.Settings)
        pathways = PhotoExcitationAutoion.Pathway[]
        for  iLevel  in  initialMultiplet.levels
            for  nLevel  in  intermediateMultiplet.levels
                for  fLevel  in  finalMultiplet.levels
                    if  Basics.selectLevelTriple(iLevel, nLevel, fLevel, settings.pathwaySelection)
                        eEnergy = nLevel.energy - iLevel.energy
                        aEnergy = nLevel.energy - fLevel.energy
                        pEnergy = fLevel.energy - iLevel.energy
                        ##x @show eEnergy, aEnergy
                        if  eEnergy < 0.   ||   aEnergy < 0    continue    end
                        rSettings = PhotoEmission.Settings( settings.multipoles, settings.gauges, false, false, LineSelection(), 0., 0., 0.)
                        eChannels = PhotoEmission.determineChannels(nLevel, iLevel, rSettings)
                        ##x @show length(eChannels)
                        aSettings = AutoIonization.Settings( false, false, LineSelection(), 0., 0., settings.maxKappa, CoulombInteraction())
                        aChannels = AutoIonization.determineChannels(fLevel, nLevel, aSettings)
                        ##x @show length(aChannels), settings.maxKappa, fLevel.energy, nLevel.energy
                        pSettings = PhotoIonization.Settings( settings.multipoles, settings.gauges, [pEnergy], false, false, false, false,
                                                              LineSelection(), ExpStokes())
                        pChannels = PhotoIonization.determineChannels(fLevel, iLevel, pSettings)
                        ##x @show length(pChannels), settings.maxKappa, fLevel.energy, iLevel.energy
                        push!( pathways, PhotoExcitationAutoion.Pathway(iLevel, nLevel, fLevel, eEnergy, aEnergy, EmProperty(0., 0.), EmProperty(0., 0.),
                                                                        eChannels, aChannels, pChannels) )
                    end
                end
            end
        end
        ##x @show length(pathways)
        return( pathways )
    end

    """
    `PhotoExcitationAutoion.displayPathways(stream::IO, pathways::Array{PhotoExcitationAutoion.Line,1}, settings::PhotoExcitationAutoion.Settings)`
        ... to display a list of pathways 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  displayPathways(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1}, settings::PhotoExcitationAutoion.Settings)
        nx = 158
        println(stream, " ")
        println(stream, "  Selected photo-excitation-autoionization pathways to be calculated:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "    ";   sb = "    "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=0);    sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(30, " Energies  " * TableStrings.inUnits("energy"); na=0);
        sb = sb * TableStrings.center(30, "  e--i         e--f    "; na=0)
        sa = sa * TableStrings.flushleft(57, "List of multipoles, gauges, kappas and total symmetries"; na=4)
        sb = sb * TableStrings.flushleft(57, "partial (multipole, gauge, total J^P)                  "; na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  pathway in pathways
            sa  = " ";    isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                          msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                          fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=2)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.excitEnergy))    * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronEnergy)) * "    "
            kappaMultipoleSymmetryList = Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry}[]
            for  ech in pathway.excitChannels
                for  ach in pathway.augerChannels
                    ##x eChannel = pathway.channels[i].excitationChannel;    aChannel = pathway.channels[i].augerChannel;
                    push!( kappaMultipoleSymmetryList, (ach.kappa, ech.multipole, ech.gauge, ach.symmetry) )
                end
            end
            wa = TableStrings.kappaMultipoleSymmetryTupels(85, kappaMultipoleSymmetryList)
            if      length(wa) == 0                       println(stream,  sa )
            elseif  length(wa) > 0    sb = sa * wa[1];    println(stream,  sb )                                 end
            for  i = 2:length(wa)   sb = TableStrings.hBlank( length(sa) ) * wa[i];    println(stream,  sb )    end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        println(stream, "\n  + Given solid angles:  $(settings.solidAngles) ")
        if settings.calcPartialCs  println(stream, "  + Calculate the partial cross sections sigma(i-e-f) at the given solid angles.")           end
        if settings.calcAngular    println(stream, "  + Calculate the angular distribution of the Auger electrons at the given solid angles.")   end
        if settings.calcFano       println(stream, "  + Calculate the Fano parameters of the cross section resonance sigma(i-e-f).")             end
        #
        return( nothing )
    end

    """
    `PhotoExcitationAutoion.displayAngular(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1})`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but
            nothing is returned otherwise.
    """
    function  displayAngular(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1})
        nx = 135
        println(stream, " ")
        println(stream, "  Angular parameters of electrons following the photo-excitation & autoionizationof resonances: ... to be adapted !!")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "    ";   sb = "    "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=0);    sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(30, " Energies  " * TableStrings.inUnits("energy"); na=0);
        sb = sb * TableStrings.center(30, "  e--i         e--f    "; na=0)
        sa = sa * TableStrings.center(10, "Multipoles"; na=3);                        sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(30, "Cou -- Cross sections -- Bab"; na=2);
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section")*"          "*
                                              TableStrings.inUnits("cross section"); na=2)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  pathway in pathways
            sa  = " ";     isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                           msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                           fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=2)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.excitEnergy))    * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronEnergy)) * "    "
            multipoles = EmMultipole[]
            for  ech in pathway.excitChannels
                multipoles = push!( multipoles, ech.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
            sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.partialCs.Coulomb))     * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.partialCs.Babushkin))   * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `PhotoExcitationAutoion.displayFanoParameters(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1})`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but
            nothing is returned otherwise.
    """
    function  displayFanoParameters(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1})
        nx = 135
        println(stream, " ")
        println(stream, "  Fano parameters of photo-excited Auger resonances: ... to be adapted !!")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "    ";   sb = "    "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=0);    sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(30, " Energies  " * TableStrings.inUnits("energy"); na=0);
        sb = sb * TableStrings.center(30, "  e--i         e--f    "; na=0)
        sa = sa * TableStrings.center(10, "Multipoles"; na=3);                        sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(30, "Cou -- Cross sections -- Bab"; na=2);
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section")*"          "*
                                          TableStrings.inUnits("cross section"); na=2)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  pathway in pathways
            sa  = " ";     isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                           msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                           fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=2)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.excitEnergy))    * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronEnergy)) * "    "
            multipoles = EmMultipole[]
            for  ech in pathway.excitChannels
                multipoles = push!( multipoles, ech.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
            sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.partialCs.Coulomb))     * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.partialCs.Babushkin))   * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `PhotoExcitationAutoion.displayPartialCs(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1})`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but
            nothing is returned otherwise.
    """
    function  displayPartialCs(stream::IO, pathways::Array{PhotoExcitationAutoion.Pathway,1})
        nx = 135
        println(stream, " ")
        println(stream, "  Partial photo-excitation & autoionization cross sections: ... to be adapted !!")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "    ";   sb = "    "
        sa = sa * TableStrings.center(23, "Levels"; na=2);            sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(23, "J^P symmetries"; na=0);    sb = sb * TableStrings.center(23, "i   --   e   --   f"; na=1);
        sa = sa * TableStrings.center(30, " Energies  " * TableStrings.inUnits("energy"); na=0);
        sb = sb * TableStrings.center(30, "  e--i         e--f    "; na=0)
        sa = sa * TableStrings.center(10, "Multipoles"; na=3);                        sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(30, "Cou -- Cross sections -- Bab"; na=2);
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section")*"          "*
                                          TableStrings.inUnits("cross section"); na=2)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  pathway in pathways
            sa  = " ";     isym = LevelSymmetry( pathway.initialLevel.J,      pathway.initialLevel.parity)
                           msym = LevelSymmetry( pathway.intermediateLevel.J, pathway.intermediateLevel.parity)
                           fsym = LevelSymmetry( pathway.finalLevel.J,        pathway.finalLevel.parity)
            sa = sa * TableStrings.center(23, TableStrings.levels_imf(pathway.initialLevel.index, pathway.intermediateLevel.index,
                                                                              pathway.finalLevel.index); na=2)
            sa = sa * TableStrings.center(23, TableStrings.symmetries_imf(isym, msym, fsym);  na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.excitEnergy))    * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", pathway.electronEnergy)) * "    "
            multipoles = EmMultipole[]
            for  ech in pathway.excitChannels
                multipoles = push!( multipoles, ech.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
            sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.partialCs.Coulomb))     * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", pathway.partialCs.Babushkin))   * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
using Test
using FlightMechanics

ψ, θ, ϕ = 0.5, 0.05, 0.2  # rad
p, q, r = 1.0, 0.2, 0.1  # rad/s

# Test regression
ear = body_angular_velocity_to_euler_angles_rates(p, q, r, θ, ϕ)
pqr = euler_angles_rates_to_body_angular_velocity(ear..., θ, ϕ)
@test isapprox(pqr, [p, q, r])

# Test against trimmer function
ψ_dot = 0.2
pqr1 = euler_angles_rates_to_body_angular_velocity(ψ_dot, 0.0, 0.0, θ, ϕ)
pqr2 = turn_rate_angular_velocity(ψ_dot, θ, ϕ)
@test isapprox(pqr1, pqr2)

# Test body_angular_velocity_to_quaternion_rates
quat = euler2quaternion(ψ, θ, ϕ)
qdot = body_angular_velocity_to_quaternion_rates(p, q, r, quat...)

Ω_inv = inv(
    [0 -p -q -r;
     p  0  r -q;
     q -r  0  p;
     r  q -p  0]
)
@test isapprox(Ω_inv * qdot * 2, quat)

# Test uvw_to_tasαβ
exp_tas = 100  # m/s
exp_α = 0.01  # rad
exp_β = 0.005  # rad
u, v, w = wind2body(tas, 0, 0, exp_α, exp_β)
tas, α, β = uvw_to_tasαβ(u, v, w)
@test isapprox(tas, exp_tas)
@test isapprox(α, exp_α)
@test isapprox(β, exp_β, atol=1e-7)

# Test tas_α_β_dot_from_uvw_dot
u = 50.0
v = 1.0
w = 5.0
u_dot = 5.0
v_dot = 0.5
w_dot = 2.0
tas_dot, α_dot, β_dot = tas_α_β_dot_from_uvw_dot(u, v, w, u_dot, v_dot, w_dot)

# Test alpha_dot against formula and tas_dot
tas = √(u*u + v*v + w*w)
@test isapprox(α_dot, (w_dot * u - w * u_dot) / (u*u + w*w))
@test isapprox(tas_dot, (u * u_dot + v * v_dot + w * w_dot) / tas)
# Test β_dot against Stevens implementation (function uses Morelli)
@test isapprox(β_dot, (v_dot * tas - v * tas_dot) / (tas * sqrt(u*u + w*w)))

# tas_α_β_dot_from_uvw_dot <--> uvw_dot_from_tas_α_β
tas, α, β = uvw_to_tasαβ(u, v, w)
u_dot_, v_dot_, w_dot_ = uvw_dot_from_tas_α_β(tas, α, β, tas_dot, α_dot, β_dot)
@test isapprox([u_dot, v_dot, w_dot], [u_dot_, v_dot_, w_dot_], atol=1e-4)

ext: jl


text:
export img_save

function img_save(img, path, filename)
    cur_path = pwd()
    img1 = Gray.(convert(Array{N0f16}, img))
    img1 = Gray.(convert.(Normed{UInt16,16}, img1))
    cd(path)
    FileIO.save(File{format"TIFF"}(filename), img1)
    cd(cur_path)
end


text:
# # Safe Feet
#
# Walk over water and air by replacing it with Dimaond Blocks
# Usage:
#

using PiCraft

function safeFeet()
   b = getTile() .+ (0, -1, 0)    #Get the Block coordinates on which the player is standing
   if getBlock(b) ∈  [PiCraft.AIR, PiCraft.WATER_STATIONARY, PiCraft.WATER_FLOWING, PiCraft.LAVA, PiCraft.LAVA_STATIONARY] #Check if the block is
       setBlock(b, PiCraft.DIAMOND_BLOCK)
   end
   sleep(0.05) #Pause as to not spam the API, results in a smoother experience
end

# Run this function in a loop, and press Ctrl-C to stop

# ```julia
#    t = now()
#    while (now() - t).value/1000 < time_period
#       safeFeet()
#    end
# ```

# ### Challenge:
# * Use a queue to remove blocks such that the maximum bridge length is 10.
# * Your script should automatically remove blocks which had too far behind.
# * Enter "\py bridge" for a demonstration


text:
#!/usr/bin/env julia
# Lists/Sequences

# Exercise 1
let
    done = false
    words = String[]

    while !done
        println("Gimme some words. I will stop on 'STOP': ")
        word = chomp(readline())
        if word == "STOP"
            done = true
        else
            push!(words, word)
        end
    end
    open("words.txt", "w") do f
        for word in words
            println(f, word)
        end
    end
end

# Exercise 2

try
    words = readlines("words.txt")
    open("words.txt", "w") do f
        for word in reverse(sort(words))
            println(f, word)
        end
    end
catch e
    msg = sprint(showerror, e)
    println(msg)
end

# Exercise 3

try
    accessions = readlines("data/ex5.acc")
    sort!(accessions)
    open("clean.acc", "w") do f
        println(f, accessions[1])
        for i in 2:length(accessions)
            if accessions[i] != accessions[i-1]
                println(f, accessions[i])
            end
        end
    end
catch e
    msg = sprint(showerror, e)
    println(msg)
end

# Exercise 4
try
    accessions = readlines("data/ex5.acc")
    sort!(accessions)
    open("clean.acc", "w") do f
        i = 2
        while i <= length(accessions)
            if accessions[i] == accessions[i-1]
                popat!(accessions, i)
            else
                i += 1
            end
        end
        for acc in accessions
            println(f, acc)
        end
    end
catch e
    msg = sprint(showerror, e)
    println(msg)
end

# Exercise 5
function ask_for_acc()
    println("Search in clean.acc: ")
    chomp(readline())
end

function linear_search_clean()
    accessions = readlines("clean.acc")
    search = ask_for_acc()
    while search != "STOP"
        if search in accessions
            println(search*" was found in the file.")
        else
            println(search*" was NOT found in the file.")
        end
        search = ask_for_acc()
    end
    println("Done.")
end

linear_search_clean()

# Exercise 6
function binary_search_clean()
    accessions = readlines("clean.acc")
    search = ask_for_acc()
    while search != "STOP"
        minimum = 1
        maximum = length(accessions)
        inlist = false

        while minimum != maximum && !inlist
            mid = convert(Int64, floor((minimum + maximum) / 2))
            if accessions[mid] < search
                minimum = mid + 1
            elseif accessions[mid] > search
                maximum = mid
            else
                inlist = true
            end
        end
        inlist ? println(search*" was found in the file.") : println(search*" was NOT found in the file.")
        println()
        search = ask_for_acc()
    end
    println("Done.")
end

binary_search_clean()

# Exercise 7.a
function print_avg(nums)
    nums = [parse(Int64, num) for num in nums]
    println(sum(nums)/length(nums))
end

if length(ARGS) == 2
    print_avg(ARGS)
else
    nums = readlines()
    if length(nums) > 0
        print_avg(nums)
    else
        println("Input was not given correctly.")
    end
end

# Exercise 7.b

function count_negative(f)
    negative_count = 0
    for ln = eachline(f)
        negative_count += count("-", ln)
    end
    println(negative_count)
    negative_count
end

if length(ARGS) == 0
    println("Enter filename:")
    filename = chomp(readline())
    open(count_negative, filename)
elseif length(ARGS) == 1
    filename = ARGS[1]
    open(count_negative, filename)
else
    println("You need to supply only one filename.")
end

# Exercise 8
using DelimitedFiles

function print_help()
    print("Arguments given incorrectly")
    print("The program must be used as such:")
    print("julia $PROGRAM_FILE -fx,y file")
    exit(1)
end

function cut()
    if length(ARGS) == 2
        filename = ARGS[2]
        if ARGS[1][begin:2] == "-f"
            cols = [parse(Int64, x) for x in split(ARGS[1][3:end], ",")]
        else
            print_help()
        end
        data = readdlm(filename, '\t', Any, '\n')
        open("$filename.cut", "w") do io
           writedlm(io, data[:, cols])
       end
    else
        print_help()
    end
end

cut()

# Exercise 9-10

using DelimitedFiles

function print_help()
    print("Arguments given incorrectly")
    print("The program must be used as such:")
    print("julia $PROGRAM_FILE file")
    exit(1)
end

function sum_cols()
    if length(ARGS) == 1
        try
            data = readdlm(ARGS[1], '\t', Float64, '\n')
            num_points, num_cols = size(data)
            avg = Vector{Float64}(undef, num_cols)
            for i in eachindex(avg)
                avg[i] = sum(view(data, :, i))/num_points
            end
            println(avg)
        catch e
            msg = sprint(showerror, e)
            println(msg)
            exit(1)
        end
    else
        print_help()
    end
end

sum_cols()


text:
module GPUifyLoops

if VERSION < v"1.1"
    @error "GPUifyLoops depends on Julia v1.1"
end

abstract type Device end
struct CPU <: Device end

abstract type GPU <: Device end
struct CUDA <: GPU end

#=
# Hopefully we can eventually support AMDGPUs through ROCm
struct ROCm <: GPU end
=#

export CPU, CUDA, Device

using StaticArrays
using Requires

export @setup, @loop, @synchronize
export @scratch, @shmem
export contextualize
export @unroll
export @launch

##
# contextualize
##
include("context.jl")

backend() = CPU()
# FIXME: Get backend from Context or have Context per backend
Cassette.overdub(ctx::Ctx, ::typeof(backend)) = CUDA()

macro launch(ex...)
    # destructure the `@launch` expression
    call = ex[end]
    kwargs = ex[2:end-1]

    device = ex[1]

    # destructure the kernel call
    if call.head != :call
        throw(ArgumentError("second argument to @launch should be a function call"))
    end

    f = call.args[1]
    args = call.args[2:end]

    quote
        $launch($(esc(device)), $(esc(f)), $(map(esc, args)...); $(map(esc, kwargs)...))
    end
end

"""
   launch(::Device, f, args..., kwargs...)

Launch a kernel on the GPU. `kwargs` are passed to `@cuda`
`kwargs` can be any of the compilation and runtime arguments
normally passed to `@cuda`.
"""
launch(::CPU, f, args...; kwargs...) = f(args...)

"""
    launch_config(::F, maxthreads, args...; kwargs...)

Calculate a valid launch configuration based on the typeof(F), the
maximum number of threads, the functions arguments and the particular
launch configuration passed to the call.

Return a NamedTuple that has `blocks`, `threads`, `shmem`, and `stream`.
All arguments are optional, but blocks and threads is recommended.
"""
function launch_config(@nospecialize(f), maxthreads, args...; kwargs...)
    return kwargs
end

function split_kwargs(kwargs)
    compiler_kws = [:minthreads, :maxthreads, :blocks_per_sm, :maxregs]
    call_kws     = [:blocks, :threads, :shmem, :stream]
    compiler_kwargs = []
    call_kwargs = []
    for kwarg in kwargs
        key, val = kwarg
        if isa(key, Symbol)
            if key in compiler_kws
                push!(compiler_kwargs, kwarg)
            elseif key in call_kws
                push!(call_kwargs, kwarg)
            else
                throw(ArgumentError("unknown keyword argument '$key'"))
            end
        else
            throw(ArgumentError("non-symbolic keyword '$key'"))
        end
    end
    return compiler_kwargs, call_kwargs
end

@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative

    function launch(::CUDA, f::F, args...; kwargs...) where F
        compiler_kwargs, call_kwargs = split_kwargs(kwargs)
        GC.@preserve args begin
            kernel_args = map(cudaconvert, args)
            kernel_tt = Tuple{map(Core.Typeof, kernel_args)...}
            kernel = cufunction(contextualize(f), kernel_tt; compiler_kwargs...)

            maxthreads = CUDAnative.maxthreads(kernel)
            config = launch_config(f, maxthreads, args...; call_kwargs...)

            kernel(kernel_args...; config...)
        end
        return nothing
    end
end

isdevice(::CPU) = false
isdevice(::Device) = true
isdevice() = isdevice(backend())

sync(::CPU) = nothing
sync() = sync(backend())

@init @require CUDAnative="be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
    using .CUDAnative
    sync(::CUDA) = CUDAnative.sync_threads()
end

@deprecate iscpu(::Val{:GPU}) isdevice()
@deprecate iscpu(::Val{:CPU}) !isdevice()
@deprecate sync(::Val{:GPU}) sync()
@deprecate sync(::Val{:CPU}) sync()

"""
    @syncronize

Calls `sync_threads()` on the GPU and nothing on the CPU.
"""
macro synchronize()
    :($sync())
end

"""
    @loop for i in (A; B)
        # body
    end

Take a `for i in (A; B)` expression and on the CPU lowers it to:

```julia
for i in A
    # body
end
```

and on the GPU:
```julia
for i in B
    if !(i in A)
        continue
    end
    # body
end
```
"""
macro loop(expr)
    if expr.head != :for
        error("Syntax error: @loop needs a for loop")
    end

    induction = expr.args[1]
    body = expr.args[2]

    if induction.head != :(=)
        error("Syntax error: @loop needs a induction variable")
    end

    rhs = induction.args[2]
    if rhs.head == :block
        @assert length(rhs.args) == 3
        # rhs[2] is a linenode
        cpuidx = rhs.args[1]
        gpuidx = rhs.args[3]

        rhs = Expr(:if, :(!$isdevice()), cpuidx, gpuidx)
        induction.args[2] = rhs

        # use cpuidx calculation to check bounds of on GPU.
        bounds_chk = quote
            if $isdevice() && !($gpuidx in $cpuidx)
                continue
            end
        end

        pushfirst!(body.args, bounds_chk)
    end

    return esc(Expr(:for, induction, body))
end

###
# Scratch and shared-memory
###
include("scratch.jl")
include("shmem.jl")

###
# Loopinfo
# - `@unroll`
###
include("loopinfo.jl")
using .LoopInfo

end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

doc = """UnitCommitment.jl Benchmark Runner

Usage:
  run.jl [-s ARG]... [-m ARG]... [-c ARG]... [-f ARG]... [options]

Examples:

    1. Benchmark all solvers, methods and formulations:

        julia run.jl

    2. Benchmark formulations "default" and "ArrCon200" using Gurobi:

        julia run.jl -s gurobi -f default -f ArrCon2000

    3. Benchmark a few test cases, using all solvers, methods and formulations:

        julia run.jl -c or-lib/20_0_1_w -c matpower/case1888rte/2017-02-01

    4. Solve 4 test cases in parallel, with 2 threads available per worker:

        JULIA_NUM_THREADS=2 julia --procs 4 run.jl

Options:
  -h --help             Show this screen.
  -s --solver=ARG       Mixed-integer linear solver (e.g. gurobi)
  -c --case=ARG         Unit commitment test case (e.g. or-lib/20_0_1_w)
  -m --method=ARG       Solution method (e.g. default)
  -f --formulation=ARG  Formulation (e.g. ArrCon2000)
  --time-limit=ARG      Time limit in seconds [default: 3600]
  --gap=ARG             Relative MIP gap tolerance [default: 0.001]
  --trials=ARG          Number of trials [default: 5]
"""

using Distributed
using Pkg
Pkg.activate(".")
@everywhere using Pkg
@everywhere Pkg.activate(".")

using DocOpt
args = docopt(doc)

@everywhere using UnitCommitment
@everywhere UnitCommitment._setup_logger()

using UnitCommitment
using Gurobi
using Logging
using JuMP

import UnitCommitment:
    ArrCon2000,
    CarArr2006,
    DamKucRajAta2016,
    Formulation,
    Gar1962,
    KnuOstWat2018,
    MorLatRam2013,
    PanGua2016,
    XavQiuWanThi2019

# Benchmark test cases
# -----------------------------------------------------------------------------
cases = [
    "pglib-uc/ca/2014-09-01_reserves_0",
    "pglib-uc/ca/2014-09-01_reserves_1",
    "pglib-uc/ca/2015-03-01_reserves_0",
    "pglib-uc/ca/2015-06-01_reserves_0",
    "pglib-uc/ca/Scenario400_reserves_1",
    "pglib-uc/ferc/2015-01-01_lw",
    "pglib-uc/ferc/2015-05-01_lw",
    "pglib-uc/ferc/2015-07-01_hw",
    "pglib-uc/ferc/2015-10-01_lw",
    "pglib-uc/ferc/2015-12-01_lw",
    "pglib-uc/rts_gmlc/2020-04-03",
    "pglib-uc/rts_gmlc/2020-09-20",
    "pglib-uc/rts_gmlc/2020-10-27",
    "pglib-uc/rts_gmlc/2020-11-25",
    "pglib-uc/rts_gmlc/2020-12-23",
    "or-lib/20_0_1_w",
    "or-lib/20_0_5_w",
    "or-lib/50_0_2_w",
    "or-lib/75_0_2_w",
    "or-lib/100_0_1_w",
    "or-lib/100_0_4_w",
    "or-lib/100_0_5_w",
    "or-lib/200_0_3_w",
    "or-lib/200_0_7_w",
    "or-lib/200_0_9_w",
    "tejada19/UC_24h_290g",
    "tejada19/UC_24h_623g",
    "tejada19/UC_24h_959g",
    "tejada19/UC_24h_1577g",
    "tejada19/UC_24h_1888g",
    "tejada19/UC_168h_72g",
    "tejada19/UC_168h_86g",
    "tejada19/UC_168h_130g",
    "tejada19/UC_168h_131g",
    "tejada19/UC_168h_199g",
    "matpower/case1888rte/2017-02-01",
    "matpower/case1951rte/2017-02-01",
    "matpower/case2848rte/2017-02-01",
    "matpower/case3012wp/2017-02-01",
    "matpower/case3375wp/2017-02-01",
    "matpower/case6468rte/2017-02-01",
    "matpower/case6515rte/2017-02-01",
]

# Formulations
# -----------------------------------------------------------------------------
formulations = Dict(
    "default" => Formulation(),
    "ArrCon2000" => Formulation(ramping = ArrCon2000.Ramping()),
    "CarArr2006" => Formulation(pwl_costs = CarArr2006.PwlCosts()),
    "DamKucRajAta2016" => Formulation(ramping = DamKucRajAta2016.Ramping()),
    "Gar1962" => Formulation(pwl_costs = Gar1962.PwlCosts()),
    "KnuOstWat2018" => Formulation(pwl_costs = KnuOstWat2018.PwlCosts()),
    "MorLatRam2013" => Formulation(ramping = MorLatRam2013.Ramping()),
    "PanGua2016" => Formulation(ramping = PanGua2016.Ramping()),
)

# Solution methods
# -----------------------------------------------------------------------------
const gap_limit = parse(Float64, args["--gap"])
const time_limit = parse(Float64, args["--time-limit"])
methods = Dict(
    "default" => XavQiuWanThi2019.Method(
        time_limit = time_limit,
        gap_limit = gap_limit,
    ),
)

# MIP solvers
# -----------------------------------------------------------------------------
optimizers = Dict(
    "gurobi" => optimizer_with_attributes(
        Gurobi.Optimizer,
        "Threads" => Threads.nthreads(),
    ),
)

# Parse command line arguments
# -----------------------------------------------------------------------------
if !isempty(args["--case"])
    cases = args["--case"]
end
if !isempty(args["--formulation"])
    formulations = filter(p -> p.first in args["--formulation"], formulations)
end
if !isempty(args["--method"])
    methods = filter(p -> p.first in args["--method"], methods)
end
if !isempty(args["--solver"])
    optimizers = filter(p -> p.first in args["--solver"], optimizers)
end
const ntrials = parse(Int, args["--trials"])

# Print benchmark settings
# -----------------------------------------------------------------------------
function printlist(d::Dict)
    for key in keys(d)
        @info "  - $key"
    end
end

function printlist(d::Vector)
    for key in d
        @info "  - $key"
    end
end

@info "Computational environment:"
@info "  - CPU: $(Sys.cpu_info()[1].model)"
@info "  - Logical CPU cores: $(length(Sys.cpu_info()))"
@info "  - System memory: $(round(Sys.total_memory() / 2^30, digits=2)) GiB"
@info "  - Available workers: $(nworkers())"
@info "  - Available threads per worker: $(Threads.nthreads())"

@info "Parameters:"
@info "  - Number of trials: $ntrials"
@info "  - Time limit (s): $time_limit"
@info "  - Relative MIP gap tolerance: $gap_limit"

@info "Solvers:"
printlist(optimizers)

@info "Methods:"
printlist(methods)

@info "Formulations:"
printlist(formulations)

@info "Cases:"
printlist(cases)

# Run benchmarks
# -----------------------------------------------------------------------------
UnitCommitment._run_benchmarks(
    cases = cases,
    formulations = formulations,
    methods = methods,
    optimizers = optimizers,
    trials = 1:ntrials,
)


text:
module api

include("types.jl")

const paths = is_apple() ? String["/System/Library/Frameworks/OpenCL.framework"] : String[]

const libopencl = Libdl.find_library(["libOpenCL", "OpenCL"], paths)
@assert libopencl != ""

function _ocl_func(func, ret_type, arg_types)
    local args_in = Symbol[Symbol("arg$i::$T")
                           for (i, T) in enumerate(arg_types.args)]
    esc(quote
        function $func($(args_in...))
            ccall(($(string(func)), libopencl),
                   $ret_type,
                   $arg_types,
                   $(args_in...))
        end
    end)
end

macro ocl_func(func, ret_type, arg_types)
    _ocl_func(func, ret_type, arg_types)
end

const CL_callback  = Ptr{Void}

abstract type CL_user_data_tag end
const CL_user_data = Ptr{CL_user_data_tag}

Base.cconvert{T}(::Type{Ptr{CL_user_data_tag}}, obj::T) = Ref{T}(obj)
Base.unsafe_convert{T}(::Type{Ptr{CL_user_data_tag}}, ref::Ref{T}) =
    Ptr{CL_user_data_tag}(isbits(T) ? pointer_from_objref(ref) : pointer_from_objref(ref[]))

Base.cconvert(::Type{Ptr{CL_user_data_tag}}, ptr::Ptr) = ptr
Base.unsafe_convert(::Type{Ptr{CL_user_data_tag}}, ptr::Ptr) = Ptr{CL_user_data_tag}(ptr)

include("api/opencl_1.0.0.jl")
include("api/opencl_1.1.0.jl")
include("api/opencl_1.2.0.jl")
include("api/opencl_2.0.0.jl")

function parse_version(version_string)
    mg = match(r"^OpenCL ([0-9]+)\.([0-9]+) .*$", version_string)
    if mg === nothing
        error("Non conforming version string: $(ver)")
    end
    return VersionNumber(parse(Int, mg.captures[1]),
                                 parse(Int, mg.captures[2]))
end

end


text:
# The simple
struct Newton{T1, Tlin} <: QuasiNewton{T1}
   approx::T1
   linsolve::Tlin
end
hasprecon(::Newton) = NoPrecon()
# struct DefaultSequence end
DefaultNewtonLinsolve(B::Number, g) = B\g
function DefaultNewtonLinsolve(B, g)
  B\g
end
function DefaultNewtonLinsolve(d, B, g)
  d .= (B\g)
end
Newton(;approx=Direct(), linsolve=DefaultNewtonLinsolve) = Newton(approx, linsolve)
summary(::Newton{<:Direct, <:typeof(DefaultNewtonLinsolve)}) = "Newton's method with default linsolve"
summary(::Newton{<:Direct, Any}) = "Newton's method with user supplied linsolve"


text:
#
# Uniform scores test statistic
#
# Mardia, K. V. (1972). A Multi-Sample Uniform Scores Test on a Circle
# and Its Parametric Competitor. Journal of the Royal Statistical
# Society. Series B (Methodological), 34(1), 102–113.
#
using StatsBase
immutable UniformScores <: PairwiseStatistic
    groups::Vector{Int32}
    ningroup::Vector{Int32}

    function UniformScores{T<:Integer}(groups::Vector{T})
        ext = extrema(groups)
        ext[1] == 1 || throw(ArgumentError("group indices must start at 1"))
        ningroup = zeros(Int32, ext[2])
        for i = 1:length(groups)
            ningroup[groups[i]] += 1
        end
        new(groups, ningroup)
    end
end
Base.eltype{T<:Real}(::UniformScores, X::AbstractArray{Complex{T}}) = T

allocwork{T<:Real}(t::UniformScores, X::AbstractVecOrMat{Complex{T}}) =
    (zeros(T, size(X)), zeros(T, ntrials(X)), zeros(Int, ntrials(X)), [cis(2π*i/size(X, 1)) for i = 1:size(X, 1)], zeros(Complex{T}, length(t.ningroup)))

# Single input matrix
using Base.Sort: DEFAULT_UNSTABLE, Perm, Forward
function computestat!{T<:Real}(t::UniformScores, out::AbstractMatrix{T},
                               work::Tuple{Matrix{T}, Vector{T}, Vector{Int}, Vector{Complex{T}}, Vector{Complex{T}}}, X::AbstractVecOrMat{Complex{T}})
    chkinput(out, X)
    angles, anglediff, ranks, circranks, d = work
    (size(X, 1) == size(angles, 1) == length(anglediff) == length(ranks) == length(circranks) == length(t.groups) &&
        size(X, 2) == size(angles, 2)) && length(d) == length(t.ningroup) || throw(ArgumentError("invalid work object"))
    for i = 1:length(X)
        @inbounds angles[i] = angle(X[i])
    end

    for k = 1:size(X, 2), j = 1:k-1
        @inbounds @simd for i = 1:size(X, 1)
            anglediff[i] = mod2pi(angles[i, j] - angles[i, k])
            ranks[i] = i
        end
        # sortperm!, but don't allocate memory
        sort!(ranks, DEFAULT_UNSTABLE, Perm(Forward, anglediff))
        fill!(d, zero(T))
        @simd for i = 1:size(X, 1)
            @inbounds d[t.groups[ranks[i]]] += circranks[i]
        end
        Wr = zero(T)
        @simd for i = 1:length(d)
            @inbounds Wr += abs2(d[i])/t.ningroup[i]
        end
        out[j, k] = 2*Wr
    end
    out
end


text:
@testset "indexing" begin
    @testset "normal use case" begin
        mktempdir() do path
            Checkpoints.config("TestPkg.bar", path)
            TestPkg.bar([1,2,3])

            index = index_checkpoint_files(path)
            @test Tables.istable(Tables.columntable(index))
            @test length(index) == 1
            entry = only(index)
            @test tags(entry) == (:date=>"2017-01-01",)
            @test entry.date == "2017-01-01"
            @test prefixes(entry) == ("TestPkg",)
            @test checkpoint_name(entry) == "bar"
            @test ==(
                checkpoint_path(entry),
                Path(joinpath(path, "date=2017-01-01", "TestPkg", "bar.jlso"))
            )
        end
    end

    @testset "files not saved by Checkpoints.jl" begin
        mktempdir(SystemPath) do path
            Checkpoints.config("TestPkg.bar", path)
            TestPkg.bar([1,2,3])
            other_file = joinpath(path, "date=2021-01-01", "other_file.txt")
            mkpath(dirname(other_file))
            write(other_file, 1)
            index = index_files(path)
            @test length(index) == 2
            @test other_file == only(checkpoint_path(entry) for entry in index if entry.date == "2021-01-01")
        end
    end

    @testset "clashing tags" begin
        mktempdir() do path
            Checkpoints.config("TestPkg.bar", path)
            with_checkpoint_tags(:tags=>"clash", :a=>"a1", :b=>"b1", :a=>"a2") do
                TestPkg.bar([1,2,3])
            end

            index = index_checkpoint_files(path)
            @test_throws ErrorException Tables.istable(Tables.columntable(index))
            @test length(index) == 1
            entry = only(index)
            @test ==(
                tags(entry),
                (:tags=>"clash", :a=>"a1", :b=>"b1", :a=>"a2", :date=>"2017-01-01",)
            )
            @test entry.date == "2017-01-01"

            @test_throws ErrorException entry.a
            @test_throws ErrorException entry.tags
        end
    end
end


text:
using JSON
include("plot.jl")
include("evo/data.jl")

function get_individual(filename::String)
    ind_dict = JSON.parse(open(filename, "r"))
    genes = Array{Array{Float64}}(ind_dict["genes"])
    fitness = Array{Float64}(ind_dict["fitness"])
    NeurodevoInd(genes, fitness)
end

function classify_performance(m::Model, X::Array{Float64}, Y::Array{Int64},
                              dataset::String, seed::Int64)
    total_time = 0
    total_memory = 0
    nout = length(unique(Y))
    df = DataFrame(seed=Int[], dataset=String[], epoch=Int[], index=Int[], label=Int[],
                   guess=Int[], certainty=Float64[], correct=Int[])
    for epoch in 1:m.cfg["n_epochs"]
        labels = zeros(Int64, length(Y))
        correct = 0
        for i in eachindex(Y)
            outputs = zeros(nout)
            for s in 1:m.cfg["n_step"]
                outputs, t, bytes, gctime, memallocs = @timed Neurodevo.step!(m, X[:, i])
            end
            labels[i] = argmax(outputs)
            if labels[i] == Y[i]
                correct += 1
                Neurodevo.reward!(m, [1.0])
            end
            push!(df, [seed, dataset, epoch, i, Y[i], labels[i],
                       outputs[labels[i]], correct])
        end
    end
    df
end

function get_performance(ind::NeurodevoInd)
    cfg = cgp_cfg(Config("cfg/evo.yaml"), ind.genes[1])
    c = cgp_controller(cfg, ind.genes[2:end])
    m = Model(cfg, c)

    X, Y = ind.func(ind.seed)
    nin = size(X, 1)
    nout = length(unique(Y))
    nsample = length(Y)
    layered_init!(m, nin, nout; nreward=1, nhidden=cfg["n_hidden"])
    df = classify_performance(m, X, Y, String(split(repr(ind.func), "_")[2]),
                              ind.seed)
    df[:step] = df[:index] + (df[:epoch] .- 1) * nsample
    df[:accuracy] = df[:correct] ./ df[:index]
    df
end

function plot_performance(res::DataFrame; title="Training accuracy",
                          filename="training.pdf")
    plt = plot(res, x=:index, y=:accuracy, color=:epoch,
               Geom.smooth,
               Guide.xlabel("Step"),
               Guide.ylabel("Accuracy"),
               Scale.color_discrete,
               Guide.title(title))
    draw(PDF(filename, 8inch, 6inch), plt)
end

function plot_all_data(ind::NeurodevoInd, id::String)
    ind.func = get_iris
    res = get_performance(ind)
    plot_performance(res; title="Iris training",
                     filename=string(id, "_iris_training.pdf"))
    # ind.func = get_diabetes
    # res = get_performance(ind)
    # plot_performance(res; title="Diabetes training",
    #                  filename=string(id, "_diabetes_training.pdf"))
    # ind.func = get_glass
    # res = get_performance(ind)
    # plot_performance(res; title="Glass training",
    #                  filename=string(id, "_glass_training.pdf"))
end

function plot_all()
    ind = get_individual("gens/1/1000/0100.dna")
    plot_all_data(ind, "1_1000")
    # ind = get_individual("gens/1/2000/0100.dna")
    # plot_all_data(ind, "1_2000")
    # ind = get_individual("gens/89462/1000/0100.dna")
    # plot_all_data(ind, "89462")
end


text:
export InDims, OutDims
const AxisDescriptorAll = Union{AxisDescriptor,String,Type{T},CubeAxis,Function} where T<:CubeAxis
using ..Cubes.Axes: get_descriptor, ByFunction
using ...YAXArrays: YAXDefaults
using DataFrames: DataFrame
using YAXArrayBase: yaxcreate

wrapWorkArray(::Type{Array},a,axes) = a
wrapWorkArray(T,a,axes) = yaxcreate(T,a,map(axname,axes),map(i->i.values,axes),nothing)
# import DataFrames
# function wrapWorkArray(t::AsDataFrame,a,cablabaxes)
#   colnames = map(Symbol,cablabaxes[2].values)
#   df = DataFrames.DataFrame(a,colnames)
#   if t.dimcol
#     df[Symbol(axname(cablabaxes[1]))]=collect(cablabaxes[1].values)
#   end
#   df
# end

abstract type ProcFilter end
struct AllMissing <: ProcFilter end
struct NValid     <: ProcFilter
  n::Int
end
struct AnyMissing <: ProcFilter end
struct AnyOcean   <: ProcFilter end
struct NoFilter   <: ProcFilter end
struct StdZero    <: ProcFilter end
struct UserFilter{F} <: ProcFilter
  f::F
end

checkskip(::NoFilter,x)         = false
checkskip(::AllMissing,x::AbstractArray)   = all(ismissing,x)
checkskip(::AllMissing,df::DataFrame)  = any(map(i->all(ismissing,getindex(df,i)),names(df)))
checkskip(::AnyMissing,x::AbstractArray)   = any(ismissing,x)
checkskip(::AnyMissing,df::DataFrame)  = any(map(i->any(ismissing,getindex(df,i)),names(df)))
checkskip(nv::NValid,x::AbstractArray)     = count(!ismissing,x) <= nv.n
checkskip(uf::UserFilter,x) = uf.f(x)
checkskip(::StdZero,x)      = all(i->i==x[1],x)
docheck(pf::ProcFilter,x)::Bool = checkskip(pf,x)
docheck(pf::Tuple,x)        = reduce(|,map(i->docheck(i,x),pf))

getprocfilter(f::Function) = (UserFilter(f),)
getprocfilter(pf::ProcFilter) = (pf,)
getprocfilter(pf::NTuple{N,<:ProcFilter}) where N = pf

"""
    InDims(axisdesc;...)

Creates a description of an Input Data Cube for cube operations. Takes a single
  or a Vector/Tuple of axes as first argument. Axes can be specified by their
  name (String), through an Axis type, or by passing a concrete axis.

### Keyword arguments

* `artype` how shall the array be represented in the inner function. Defaults to `AsArray`, alternatives are `AsDataFrame` or `AsAxisArray`
* `filter` define some filter to skip the computation, e.g. when all values are missing. Defaults to
    `AllMissing()`, possible values are `AnyMissing()`, `AnyOcean()`, `StdZero()`, `NValid(n)`
    (for at least n non-missing elements). It is also possible to provide a custom one-argument function
    that takes the array and returns `true` if the compuation shall be skipped and `false` otherwise.
"""
mutable struct InDims
  axisdesc::Tuple
  artype
  procfilter::Tuple
end
function InDims(axisdesc::AxisDescriptorAll...; artype=Array, filter = AllMissing())
  descs = map(get_descriptor,axisdesc)
  any(i->isa(i,ByFunction),descs) && error("Input cubes can not be specified through a function")
  InDims(descs,artype,getprocfilter(filter))
end

struct OutDims
  axisdesc
  backend::Symbol
  backendargs
  update::Bool
  artype
  chunksize::Any
  outtype::Union{Int,DataType}
end
"""
    OutDims(axisdesc;...)

Creates a description of an Output Data Cube for cube operations. Takes a single
  or a Vector/Tuple of axes as first argument. Axes can be specified by their
  name (String), through an Axis type, or by passing a concrete axis.

- `axisdesc`: List of input axis names
- `backend` : specifies the dataset backend to write data to, must be either :auto or a key in `YAXArrayBase.backendlist`
- `update` : specifies wether the function operates inplace or if an output is returned
- `artype` : specifies the Array type inside the inner function that is mapped over
- `chunksize`: A Dict specifying the chunksizes for the output dimensions of the cube, or `:input` to copy chunksizes from input cube axes or `:max` to not chunk the inner dimensions
- `outtype`: force the output type to a specific type, defaults to `Any` which means that the element type of the first input cube is used
"""
function OutDims(axisdesc...;
           backend=:auto,
           update=false,
           artype=Array,
           chunksize=YAXDefaults.chunksize[],
           outtype=1,
           backendargs...)
  descs = get_descriptor.(axisdesc)
  OutDims(descs,backend,backendargs,update,artype,chunksize,outtype)
end

registerDATFunction(a...;kwargs...)=@warn("Registration does not exist anymore, ignoring....")


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 5323b5ac-0c51-4952-8f27-4f489c1f2d1a
begin
    import Pkg
    Pkg.activate() # disable Pluto's package management
end

# ╔═╡ d8aeb9f4-990c-11eb-0dbb-078cda4928d5
begin
    using GSL # provides the special functions
    using GR
    using PlutoUI
end

# ╔═╡ 6127fb6a-22fc-4f35-8228-cc06ffebf9c0
begin
    a0 = 5.291772109217e-11

    # The unitless radial coordinate
    ρ(r, n) = 2r/(n*a0)
end

# ╔═╡ 7aa443bd-b706-41e2-bdc8-0d8c10c99b02
# The θ and ϕ dependence
function Yml(m::Int, l::Int, θ::Real, ϕ::Real)
    (-1.0)^m*sf_legendre_Plm(l, abs(m), cos(θ))*ℯ^(im*m*ϕ)
end

# ╔═╡ 778c4fa2-68c0-494f-ae22-477189a27f2b
# The Radial dependence
function R(n::Int, l::Int, ρ::Real)
    sf_laguerre_n(n-l-1, 2*l+1, ρ)*ℯ^(-ρ/2)*ρ^l
end

# ╔═╡ cf4d7a50-afea-4f10-b08d-0f0150b7da7a
# A normalization: This is dependent on the choice of polynomial representation
function norm(n::Int, l::Int)
    sqrt((2/n)^3 * factorial(n-l-1)/(2n*factorial(n+l)))
end

# ╔═╡ d3ff17f0-3dd6-4ffe-a045-5b573aec02a5
# Generates an Orbital Funtion of (r, θ, ϕ) for a specificied n, l, m.
function Orbital(n::Int, l::Int, m::Int)
    # we make sure l and m are within proper bounds
    if l > n || abs(m) > l
        throw(DomainError())
    end
    Ψ(ρ, θ, ϕ) = norm(n, l) * R(n, l, ρ) * abs(Yml(m, l, θ, ϕ))
    Ψ
end

# ╔═╡ 94158386-447f-40cd-a321-8b333053d558
function CarttoSph(x::Array, y::Array, z::Array)
    r = sqrt.(x.^2+y.^2+z.^2)
    θ = acos.(z./r)
    ϕ = atan.(y./x)
    r, θ, ϕ
end

# ╔═╡ ad5ffb61-b2cc-499a-87a4-ad71833de766
function calculate_electronic_density(n, l, m)
    N = 50
    r = 1e-9
    s = LinRange(-r, r, N)
    x, y, z = meshgrid(s, s, s)

    r, θ, ϕ = CarttoSph(x, y, z)
    Ψ = Orbital(n, l, m)

    Ψv = zeros(Float32, N, N, N)
    for i in 1:N
        for j in 1:N
            for k in 1:N
                Ψv[i, j, k] = abs(Ψ(ρ(r[i, j, k], n), θ[i, j, k], ϕ[i, j, k]))
            end
        end
    end

    (Ψv.-minimum(Ψv))./(maximum(Ψv).-minimum(Ψv))
end

# ╔═╡ a4c51114-4bf9-4d40-b123-1629a3be54e1
@bind Ψ Select(["2,0,0", "3,0,0", "2,1,0", "3,1,0", "3,1,1", "2,1,1", "3,2,0", "3,2,1", "3,2,2", "4,0,0", "4, 1,0", "4,2,0", "4,2,1", "4,2,2", "4,3,0", "4,3,1", "4,3,2", "4,3,3"], default="3,2,0")

# ╔═╡ 11cfffdf-91d6-4909-bb59-e8d12a95e5e8
@bind ϕ Slider(0:360, default=30)

# ╔═╡ 9899f86c-2fba-4d72-abab-acbef1dc6988
@bind iso Slider(0.1:0.05:0.8, default=0.2)

# ╔═╡ 16697f21-f879-48d1-ad41-f39e1dde5830
begin
    n, m, l = parse.(Int, split(Ψ, ','))
    Ψv = calculate_electronic_density(n, m, l)
    isosurface(Ψv, isovalue=iso, rotation=ϕ)
end

# ╔═╡ Cell order:
# ╠═5323b5ac-0c51-4952-8f27-4f489c1f2d1a
# ╠═d8aeb9f4-990c-11eb-0dbb-078cda4928d5
# ╠═6127fb6a-22fc-4f35-8228-cc06ffebf9c0
# ╠═7aa443bd-b706-41e2-bdc8-0d8c10c99b02
# ╠═778c4fa2-68c0-494f-ae22-477189a27f2b
# ╠═cf4d7a50-afea-4f10-b08d-0f0150b7da7a
# ╠═d3ff17f0-3dd6-4ffe-a045-5b573aec02a5
# ╠═94158386-447f-40cd-a321-8b333053d558
# ╠═ad5ffb61-b2cc-499a-87a4-ad71833de766
# ╠═a4c51114-4bf9-4d40-b123-1629a3be54e1
# ╠═11cfffdf-91d6-4909-bb59-e8d12a95e5e8
# ╠═9899f86c-2fba-4d72-abab-acbef1dc6988
# ╠═16697f21-f879-48d1-ad41-f39e1dde5830


text:
# ------------------------------------------------------------
# Functions and types for interacting with interior filters

immutable InteriorFilter
	van_moment::Int64
	support::DaubSupport
	filter::Vector{Float64}

	function InteriorFilter(p, support, filter)
		if p >= 0
			new(p, support, filter)
		else
			throw(AssertionError("Not a valid interior filter"))
		end
	end
end

function Base.show(io::IO, IF::InteriorFilter)
	S = support(IF)
	println(io, "Filter for Daubechies ", van_moment(IF), " scaling function on [", left(S), ", ", right(S), "]:")

	show(io, IF.filter)
end

@doc """
	ifilter(p::Int)

Internal Daubechies filter with `p` vanishing moments with the **normal** filters.

	ifilter(p::Int, true)

Internal Daubechies filter with `p` vanishing moments with the symmlet filters.
This is the default.
"""->
function ifilter(p::Integer, symmlet::Bool=true)
	p < 1 && throw(DomainError())

	if symmlet && 1 < p <= 8
		return InteriorFilter(p, DaubSupport(-p+1,p), INTERIOR_FILTERS[p])
	else
		return InteriorFilter(p, DaubSupport(0,2*p-1), wavelet(WT.Daubechies{p}()).qmf)
	end
end

for func in [:support, :van_moment]
	@eval begin
		function $func(C::InteriorFilter)
			C.$func
		end
	end
end

# Returns a copy; otherwise e.g. scale!(coef(C), 2) will modify C.filter
coef(C::InteriorFilter) = copy(C.filter)

function Base.length(C::InteriorFilter)
	length(coef(C))
end

# ------------------------------------------------------------
# Functions and types for interacting with boundary filters

immutable BoundaryFilter
	side::Char
	van_moment::Int64
	support::DaubSupport
	filter::Array{Vector{Float64}}

	function BoundaryFilter(side, p, S, F)
		if (side == 'L' || side == 'R') && (2 <= p <= 8)
			new(side, p, S, F)
		else
			throw(AssertionError("Not a valid boundary filter"))
		end
	end
end

left(B::BoundaryFilter) = left(B.support)
right(B::BoundaryFilter) = right(B.support)
side(B::BoundaryFilter) = B.side

function Base.show(io::IO, B::BoundaryFilter)
	S = support(B)
	p = van_moment(B)

	side = (B.side == 'L' ? "left" : "right")
	print(io, "Filters for ", side, " Daubechies ", p, " scaling function on [", left(S), ", ", right(S), "]:")

	for k in 0:p-1
		print(io, "\nk = ", k, ": ")
		show(io, bfilter(B,k))
	end
end

@doc """
	integers(B::BoundaryFilter)

The non-zero integers in the support of the boundary scaling function with filter `B`.
"""->
function integers(B::BoundaryFilter)
	if side(B) == 'L'
		return right(support(B)):-1:1
	else side(B) == 'R'
		# The explicit step ensures type stability
		return left(support(B)):1:-1
	end
end

@doc """
	bfilter(p::Int, boundary::Char) -> BoundaryFilter

Return the boundary filters for the scaling functions with `p` vanishing moments.
`N` can be between 2 and 8.

`boundary` is either `'L'` or `'R'`.
"""->
function bfilter(p::Integer, boundary::Char)
	2 <= p <= 8 || throw(AssertionError())

	if boundary == 'L'
		 return BoundaryFilter('L', p, DaubSupport(0, 2*p-1), LEFT_SCALING_FILTERS[p])
	elseif boundary == 'R'
		 return BoundaryFilter('R', p, DaubSupport(-2*p+1, 0), RIGHT_SCALING_FILTERS[p])
	else
		error("Boundary must be 'L' or 'R'")
	end
end

@doc """
	bfilter(BoundaryFilter, k::Int) -> Vector

Return the boundary filter for the `k`'th scaling function (0 <= `k` < the number of vanishing moments).
"""->
function bfilter(B::BoundaryFilter, k::Int)
	0 <= k < van_moment(B) || throw(DomainError())
	return B.filter[k+1]
end

@doc """
	van_moment(F::BoundaryFilter) -> Integer

Return the number of vanishing moments of the boundary scaling functions defined by `F`.
"""->
function van_moment(F::BoundaryFilter)
	return F.van_moment
end

@doc """
	support(B::BoundaryFilter)

Union of the supports of the boundary scaling functions defined by the filters `B`.
"""->
function support(B::BoundaryFilter)
	B.support
end

@doc """
	support(B::BoundaryFilter, k)

Support of the `k`'th boundary scaling function defined by the filters `B`.
"""->
function support(B::BoundaryFilter, k::Integer)
	0 <= k < (vm = van_moment(B)) || throw(DomainError())
	if B.side == 'L'
		return DaubSupport(0, vm+k)
	elseif B.side == 'R'
		return DaubSupport(-vm-k, 0)
	end
end

# ------------------------------------------------------------
# Boundary low pass filter coefficients for Daubechies wavelets
# http://www.pacm.princeton.edu/~ingrid/publications/54.txt

const INTERIOR_FILTERS = Dict{Int, Vector{Float64}}(
2 => [ 0.482962913145 ; 0.836516303738 ; 0.224143868042 ; -0.129409522551 ]
,
3 => [ .332670552950 ; .806891509311 ; .459877502118 ; -.135011020010 ; -.085441273882 ; .035226291882 ]
,
4 => [ 0.045570345896 ; -0.0178247014417 ; -0.140317624179 ; 0.421234534204 ; 1.13665824341 ; 0.703739068656 ; -0.0419109651251 ; -0.107148901418 ] / sqrt2
,
5 => [ 0.0276321529578 ; -0.0298424998687 ; -0.247951362613 ; 0.0234789231361 ; 0.89658164838 ; 1.02305296689 ; 0.281990696854 ; -0.0553441861166 ; 0.0417468644215 ; 0.0386547959548 ] / sqrt2
,
6 => [ -0.0110318675094 ; 0.00249992209279 ; 0.06325056266 ; -0.0297837512985  ; -0.102724969862 ; 0.477904371333 ; 1.11389278393 ; 0.694457972958 ; -0.0683231215866  ; -0.166863215412 ; 0.00493661237185 ; 0.0217847003266 ] / sqrt2
,
7 => [ 0.014521394762 ; 0.00567134268574 ; -0.152463871896 ; -0.198056706807 ; 0.408183939725 ; 1.08578270981 ; 0.758162601964  ; 0.0246656594886  ; -0.070078291222 ; 0.0960147679355 ; 0.043155452582 ; -0.0178704316511 ; -0.0014812259146 ; 0.0037926585342 ] / sqrt2
,
8 => [ 0.00267279339281 ; -0.000428394300246 ; -0.0211456865284 ; 0.00538638875377 ; 0.0694904659113 ; -0.0384935212634 ; -0.0734625087609 ; 0.515398670374 ; 1.09910663054 ; 0.68074534719 ; -0.0866536154058 ; -0.202648655286 ; 0.0107586117505 ; 0.0448236230437 ; -0.000766690896228 ; -0.0047834585115 ] / sqrt2
)

const LEFT_SCALING_FILTERS = Dict{Int, Array{Vector}}(
2 => Any[
[ 0.6033325119E+00 ; 0.6908955318E+00 ; -0.3983129977E+00 ]
,
[ 0.3751746045E-01 ; 0.4573276599E+00 ; 0.8500881025E+00 ; 0.2238203570E+00 ; -0.1292227434E+00 ]
]
,
3 => Any[
[ 0.3888997639E+00 ; -0.8820782813E-01 ; -0.8478413085E+00 ; 0.3494874367E+00 ]
,
[ -0.6211483178E+00 ; 0.5225273932E+00 ; -0.2000080031E+00 ; 0.3378673486E+00 ; -0.3997707705E+00 ; 0.1648201297E+00 ]
,
[ -0.9587863831E-02 ; 0.3712255319E-03 ; 0.3260097101E+00 ; 0.8016481645E+00 ; 0.4720552620E+00 ; -0.1400420809E+00 ; -0.8542510010E-01 ; 0.3521962365E-01 ]
]
,
4 => Any[
[ 0.9097539258E+00 ; 0.4041658894E+00 ; 0.8904031866E-01 ; -0.1198419201E-01 ; -0.3042908414E-01 ]
,
[ -0.2728514077E+00 ; 0.5090815232E+00 ; 0.6236424434E+00 ; 0.4628400863E+00 ; 0.2467476417E+00 ; -0.1766953329E-01 ; -0.4517364549E-01 ]
,
[ 0.1261179286E+00 ; -0.2308557268E+00 ; -0.5279923525E-01 ; 0.2192651713E+00 ; 0.4634807211E+00 ; 0.7001197140E+00 ; 0.4120325790E+00 ; -0.2622276250E-01 ; -0.6704069413E-01 ]
,
[ -0.2907980427E-01 ; 0.5992807229E-01 ; 0.6176427778E-02 ; -0.4021099904E-01 ; -0.3952587013E-01 ; -0.5259906257E-01 ; 0.3289494480E+00 ; 0.7966378967E+00 ; 0.4901130336E+00 ; -0.2943287768E-01 ; -0.7524762313E-01 ]
]
,
5 => Any[
[ 0.9302490657E+00 ; 0.3488878121E+00 ; 0.1098578445E+00 ; 0.2701025958E-01 ; 0.7897329608E-02 ; 0.7300852033E-02 ]
,
[ -0.3099501188E+00 ; 0.5987214008E+00 ; 0.6208411495E+00 ; 0.3720164550E+00 ; 0.1451060379E+00 ; 0.6263262144E-02 ; 0.1687072021E-01 ; 0.1562115518E-01 ]
,
[ 0.1138401479E+00 ; -0.3641831811E+00; 0.1228916167E-01 ; 0.4117590210E+00 ; 0.5954973448E+00 ; 0.5457944280E+00 ; 0.1756416185E+00 ; -0.1751671276E-01 ; 0.2376581050E-01 ; 0.2200554624E-01 ]
,
[ -0.4098816981E-01 ; 0.1431055982E+00 ; -0.6777711666E-01 ; -0.1740980466E+00 ; -0.9559132536E-01 ; 0.1512981220E+00 ; 0.6357421629E+00 ; 0.6823643465E+00 ; 0.1963091141E+00 ; -0.3301996982E-01 ; 0.2833562039E-01 ; 0.2623688365E-01 ]
,
[ 0.7965029532E-02 ; -0.2796122002E-01 ; 0.1706238600E-01 ; 0.2928655475E-01 ; 0.2857806991E-02 ; -0.4686793277E-01 ; -0.1702677471E+00 ; 0.3028147873E-01 ; 0.6351733467E+00 ; 0.7207789094E+00 ; 0.1995469939E+00 ; -0.3864256419E-01 ; 0.2947288315E-01 ; 0.2728991267E-01 ]
]
,
6 => Any[
[ 0.9231184460E+00 ; 0.3781064509E+00 ; 0.6815972396E-01 ; -0.1063392935E-01 ; -0.9724819798E-02 ; 0.1297887318E-02 ; 0.5723757426E-02 ]
,
[ -0.2868924125E+00 ; 0.5940469387E+00 ; 0.6609604336E+00 ; 0.3412988002E+00 ; 0.9716877076E-01 ; -0.1860855577E-01 ; -0.4013700006E-01 ; 0.1537307481E-02 ; 0.6783960390E-02 ]
,
[ 0.1472870932E+00 ; -0.3468945123E+00 ; 0.6221043773E-01 ; 0.4938022108E+00 ; 0.5849105895E+00 ; 0.4598788833E+00 ; 0.2274224529E+00 ; -0.2788953680E-01 ; -0.6367943024E-01 ; 0.2177711568E-02 ; 0.9609989671E-02 ]
,
[ -0.8324138021E-01 ; 0.2035630921E+00 ; -0.8289131955E-01 ; -0.1939904105E+00 ; -0.1908457968E-01 ; 0.2636403540E+00 ; 0.5123094150E+00 ; 0.6538041187E+00 ; 0.3695601186E+00 ; -0.3993953251E-01 ; -0.9479455728E-01 ; 0.2986558645E-02 ; 0.1317933842E-01 ]
,
[ 0.3164853364E-01 ; -0.8214471961E-01 ; 0.4327551059E-01 ; 0.6830653932E-01 ; -0.4840048148E-02 ; -0.7320834090E-01 ; -0.5786918156E-01 ; 0.2055855208E-01 ; 0.4001944206E+00 ; 0.7644983569E+00 ; 0.4653802512E+00 ; -0.4686421163E-01 ; -0.1137136979E+00 ; 0.3413297317E-02 ; 0.1506248690E-01 ]
,
[ -0.5550804971E-02 ; 0.1471087926E-01 ; -0.8859117912E-02 ; -0.1158923365E-01 ; 0.1569832259E-02 ; 0.1157111171E-01 ; 0.7445526684E-02 ; 0.2902493903E-01 ; -0.2999940975E-01 ; -0.6576340468E-01 ; 0.3438258786E+00 ; 0.7867130346E+00 ; 0.4895844920E+00 ; -0.4825468670E-01 ; -0.1177993251E+00 ; 0.3488475665E-02 ; 0.1539424027E-01 ]
]
,
7 => Any[
[ 0.9426568713E+00 ; 0.2691475418E+00 ; 0.1685897688E+00 ; 0.9379442214E-01 ; 0.4001600171E-01 ; 0.1057520816E-01 ; -0.1804517915E-02 ; 0.4620417362E-02 ]
,
[ -0.3237338194E+00 ; 0.6746802222E+00 ; 0.5985769136E+00 ; 0.2331427899E+00 ; 0.7615188825E-02 ; 0.8787238208E-01 ; 0.1327822780E+00 ; -0.4208296665E-01 ; -0.4763766998E-02 ; 0.1219755973E-01 ]
,
[ -0.2713311808E-01 ; -0.2899110477E+00 ; 0.6661404639E-01 ; 0.7251175151E+00 ; 0.5750729858E+00 ; 0.4362186445E-01 ; -0.1857980605E+00 ; 0.1201045634E+00 ; 0.5011061230E-01 ; -0.2991531277E-01 ; -0.1573393605E-02 ; 0.4028652635E-02 ]
,
[ 0.2428297386E-01 ; -0.1012534619E-01 ; -0.1589124403E+00 ; -0.1400872185E+00 ; 0.3332116140E+00 ; 0.7053065273E+00 ; 0.5793528972E+00 ; 0.3896614469E-02 ; -0.5955471765E-01 ; 0.7479188863E-01 ; 0.3081374376E-01 ; -0.1340206604E-01 ; -0.1047350551E-02 ; 0.2681726647E-02 ]
,
[ 0.2396415492E-02 ; -0.2199441967E-01 ; 0.3831273498E-01 ; 0.8416826933E-02 ; -0.1404553727E+00 ; -0.9975512873E-01 ; 0.2671255526E+00 ; 0.7725679637E+00 ; 0.5400947777E+00 ; 0.1395136806E-01 ; -0.4947167749E-01 ; 0.6809969285E-01 ; 0.3044230228E-01 ; -0.1259630803E-01 ; -0.1045026849E-02 ; 0.2675776840E-02 ]
,
[ -0.9493922980E-03 ; 0.8445740060E-02 ; -0.8462283108E-02 ; -0.3100522289E-02 ; 0.2269052507E-01 ; -0.9675176425E-02 ; -0.1013545517E+00 ; -0.1405275956E+00 ; 0.2877816730E+00 ; 0.7686169120E+00 ; 0.5358841401E+00 ; 0.1743190445E-01 ; -0.4949055461E-01 ; 0.6785067390E-01 ; 0.3050565786E-01 ; -0.1262873733E-01 ; -0.1047102354E-02 ; 0.2681091142E-02 ]
,
[ 0.1364120169E-03 ; -0.1153014522E-02 ; 0.1050196147E-02 ; 0.6396328807E-03 ; -0.1940525665E-02 ; 0.1962825019E-02 ; 0.9431646581E-02 ; 0.4027633417E-02 ; -0.1076854888E+00 ; -0.1401281648E+00 ; 0.2886799721E+00 ; 0.7677614604E+00 ; 0.5360881498E+00 ; 0.1744975917E-01 ; -0.4954996630E-01 ; 0.6789074612E-01 ; 0.3051526119E-01 ; -0.1263602514E-01 ; -0.1047379021E-02 ; 0.2681799545E-02 ]
]
,
8 => Any[
[ 0.9281136260E+00 ; 0.3668618286E+00 ; 0.6266997447E-01 ; -0.7647264163E-02 ; -0.5088004810E-02 ; 0.1171898539E-02 ; 0.1656055473E-02 ; -0.1929333000E-03 ; -0.1203653707E-02 ]
,
[ -0.2926634432E+00 ; 0.6373078980E+00 ; 0.6452435440E+00 ; 0.2965161471E+00 ; 0.5906254343E-01 ; -0.1366965641E-01 ; -0.1323318522E-01 ; 0.2213708145E-02 ; 0.8306912600E-02 ; -0.1892199826E-03 ; -0.1180561737E-02 ]
,
[ 0.1537394280E+00 ; -0.3933877381E+00 ; 0.1368645026E+00 ; 0.6064470308E+00 ; 0.5685338138E+00 ; 0.3166230149E+00 ; 0.9637901837E-01 ; -0.2215743488E-01 ; -0.4105959919E-01 ; 0.3062891878E-02 ; 0.1191510029E-01 ; -0.2473506634E-03 ; -0.1543244666E-02 ]
,
[ -0.9143117079E-01 ; 0.2571933236E+00 ; -0.1858137201E+00 ; -0.2318000501E+00 ; 0.1598859051E+00 ; 0.4871009258E+00 ; 0.5734443542E+00 ; 0.4450618899E+00 ; 0.2033090498E+00 ; -0.3331081673E-01 ; -0.6859620275E-01 ; 0.4406096842E-02 ; 0.1763488599E-01 ; -0.3388224811E-03 ; -0.2113946167E-02 ]
,
[ 0.5709470789E-01 ; -0.1665850443E+00 ; 0.1454047216E+00 ; 0.1015538648E+00 ; -0.1362852159E+00 ; -0.1769345243E+00 ; 0.8596479652E-02 ; 0.3062987789E+00 ; 0.5431752460E+00 ; 0.6140955155E+00 ; 0.3284966144E+00 ; -0.4723534035E-01 ; -0.1039762360E+00 ; 0.6053819361E-02 ; 0.2473853830E-01 ; -0.4480360998E-03 ; -0.2795340476E-02 ]
,
[ -0.2771885178E-01 ; 0.8389262409E-01 ; -0.8147463473E-01 ; -0.3866640174E-01 ; 0.7464550842E-01 ; 0.6920099668E-01 ; -0.2231653624E-01 ; -0.1020129060E+00 ; -0.6275852804E-01 ; 0.8683887629E-01 ; 0.4540741128E+00 ; 0.7328223558E+00 ; 0.4319437621E+00 ; -0.5734179466E-01 ; -0.1314896433E+00 ; 0.7195965493E-02 ; 0.2979946085E-01 ; -0.5190323961E-03 ; -0.3238293222E-02 ]
,
[ 0.8529185218E-02 ; -0.2624362023E-01 ; 0.2697489980E-01 ; 0.9748173003E-02 ; -0.2431726904E-01 ; -0.1901319060E-01 ; 0.8695522469E-02 ; 0.2729552258E-01 ; 0.1231155208E-01 ; 0.9187597846E-02 ; -0.4745631970E-01 ; -0.2924026466E-01 ; 0.3831732042E+00 ; 0.7723764353E+00 ; 0.4743837775E+00 ; -0.6082976218E-01 ; -0.1418145038E+00 ; 0.7565491348E-02 ; 0.3149214931E-01 ; -0.5401012058E-03 ; -0.3369743560E-02 ]
,
[ -0.1148179452E-02 ; 0.3567347972E-02 ; -0.3784785333E-02 ; -0.1105680377E-02 ; 0.3387180654E-02 ; 0.2395725348E-02 ; -0.1281174092E-02 ; -0.3407989030E-02 ; -0.1283416571E-02 ; -0.9841020551E-02 ; 0.6341196709E-02 ; 0.4619539219E-01 ; -0.2943416082E-01 ; -0.5094613096E-01 ; 0.3655480429E+00 ; 0.7770499808E+00 ; 0.4810695095E+00 ; -0.6126039684E-01 ; -0.1432436897E+00 ; 0.7606483137E-02 ; 0.3168981524E-01 ; -0.5420982042E-03 ; -0.3382203026E-02 ]
]
)

const RIGHT_SCALING_FILTERS = Dict{Int, Array{Vector}}(
2 => Any[
[ 0.8705087534E+00 ; 0.4348969980E+00 ; 0.2303890438E+00 ]
,
[ -0.1942334074E+00 ; 0.1901514184E+00 ; 0.3749553316E+00 ; 0.7675566693E+00 ; 0.4431490496E+00 ]
]
,
3 => Any[
[ 0.9096849943E+00 ; 0.3823606559E+00 ; 0.1509872153E+00 ; 0.5896101069E-01 ]
,
[ -0.2904078511E+00 ; 0.4189992290E+00 ; 0.4969643721E+00 ; 0.4907578307E+00 ; 0.4643627674E+00 ; 0.1914505442E+00 ]
,
[ 0.8183541840E-01 ; -0.1587582156E+00 ; -0.9124735623E-01 ; 0.6042558204E-03 ; 0.7702933610E-01 ; 0.5200601778E+00 ; 0.7642591993E+00 ; 0.3150938230E+00 ]
]
,
4 => Any[
[ 0.9154705188E+00 ; 0.3919142810E+00 ; 0.5947771124E-01 ; -0.2519180851E-01 ; 0.6437934569E-01 ]
,
[ -0.2191626469E+00 ; 0.4488001781E+00 ; 0.7540005084E+00 ; 0.3937758157E+00 ; -0.1581338944E+00 ; -0.1614201190E-01 ; 0.4126840881E-01 ]
,
[ 0.1290078289E-01 ; -0.1390716006E+00 ; 0.2921367950E-01 ; 0.4606168537E+00 ; 0.8164119742E+00 ; 0.2986473346E+00 ; -0.1027663536E+00 ; -0.1257488211E-01 ; 0.3214874197E-01 ]
,
[ -0.6775603652E-02 ; 0.1913244129E-01 ; -0.1770918425E-01 ; -0.6765916174E-01 ; -0.3023588481E-01 ; 0.4977920821E+00 ; 0.8039495996E+00 ; 0.2977111011E+00 ; -0.9910804055E-01 ; -0.1259895190E-01 ; 0.3221027840E-01 ]
]
,
5 => Any[
[ 0.6629994791E+00 ; 0.5590853114E+00 ; -0.3638895712E+00 ; -0.3190292420E+00 ; -0.8577516141E-01 ; 0.7938922949E-01 ]
,
[ -0.1615976602E+00 ; 0.7429922432E+00 ; 0.5771521705E+00 ; 0.2773360273E+00 ; 0.1530007177E-01 ; -0.1063942878E+00 ; -0.1216769629E-01 ; 0.1126647052E-01 ]
,
[ 0.2576681838E+00 ; -0.3995181391E-01 ; 0.1290658921E-01 ; 0.4163708782E+00 ; 0.6469374504E+00 ; 0.5608421202E+00 ; 0.4995700940E-02 ; -0.1569482750E+00 ; -0.2009485730E-01 ; 0.1860648984E-01 ]
,
[ -0.8955158953E-01 ; 0.3895616532E-01 ; 0.4751157158E-02 ; -0.6596188737E-01 ; 0.1520881337E-01 ; 0.2466644601E+00 ; 0.7199998123E+00 ; 0.6131240289E+00 ; 0.1537319456E-01 ; -0.1721947543E+00 ; -0.2083858871E-01 ; 0.1929513523E-01 ]
,
[ 0.1546583496E-01 ; -0.6559247482E-02 ; 0.2728093541E-02 ; 0.1830645747E-01 ; 0.1708296448E-01 ; 0.2156026058E-01 ; -0.3692524465E-01 ; 0.2045334663E+00 ; 0.7233635779E+00 ; 0.6327590443E+00 ; 0.1649999072E-01 ; -0.1751674471E+00 ; -0.2109311508E-01 ; 0.1953080958E-01 ]
]
,
6 => Any[
[ 0.9236675275E+00 ; 0.3793343885E+00 ; 0.4986459778E-01 ; -0.1688947163E-01 ; -0.3311806998E-02 ; 0.2804629640E-02 ; -0.1237553624E-01 ]
,
[ -0.2831955072E+00 ; 0.6209232596E+00 ; 0.6494435367E+00 ; 0.3279016732E+00 ; 0.4613465423E-01  ; -0.2201105042E-01 ; 0.4772996268E-01 ; 0.1811083865E-02 ; -0.7992103956E-02 ]
,
[ 0.1415889950E+00 ; -0.3264117354E+00 ; 0.5898718517E-01 ; 0.5662575494E+00 ; 0.6511436245E+00 ; 0.3434709841E+00 ; -0.6897565912E-01 ; -0.2139852264E-01 ; 0.4514540543E-01 ; 0.1807019829E-02 ; -0.7974169836E-02 ]
,
[ -0.4080929146E-01 ; 0.1251977353E+00 ; -0.9208945297E-01 ; -0.1815441604E+00 ; 0.7492893632E-01 ; 0.4921890166E+00 ; 0.7598885402E+00 ; 0.3298535864E+00 ; -0.6885203537E-01 ; -0.2056224689E-01 ; 0.4357038095E-01 ; 0.1729458762E-02 ; -0.7631901801E-02 ]
,
[ 0.1473338385E-01 ; -0.3879928306E-01 ; 0.3222685441E-01 ; 0.4181351746E-01 ; -0.4470160402E-01 ; -0.1143592801E+00 ; -0.3079264597E-01 ; 0.4953303341E+00 ; 0.7831063037E+00 ; 0.3365261769E+00 ; -0.7161554542E-01 ; -0.2096645980E-01 ; 0.4445995152E-01 ; 0.1762244193E-02 ; -0.7776580122E-02 ]
,
[ -0.2279778462E-02 ; 0.5934267127E-02 ; -0.4653920620E-02 ; -0.4697326198E-02 ; 0.7905498252E-02 ; 0.1478426761E-01 ; 0.6356609850E-03 ; -0.1184137622E+00 ; -0.4719524338E-01 ; 0.4914056494E+00 ; 0.7873827096E+00 ; 0.3378777169E+00 ; -0.7256793624E-01 ; -0.2105690525E-01 ; 0.4471373463E-01 ; 0.1767573251E-02 ; -0.7800096641E-02 ]
]
,
7 => Any[
[ 0.9091449027E+00 ; 0.4122956702E+00 ; 0.8463996964E-02 ; -0.4985469571E-01 ; -0.2885243693E-01 ; -0.8468475604E-02 ; 0.9290995299E-03 ; 0.2377975298E-02 ]
,
[ -0.2699308042E+00 ; 0.6287804833E+00 ; 0.6425124921E+00 ; 0.3295936378E+00 ; 0.9680040842E-01 ; -0.1505090528E-02 ; -0.2576996713E-01 ; -0.1663793375E-01 ; 0.7886583178E-03 ; 0.2019348750E-02 ]
,
[ 0.2157402688E+00 ; -0.3660961628E+00 ; 0.9979560410E-01 ; 0.5343306371E+00 ; 0.5651956498E+00 ; 0.3984651316E+00 ; 0.2006608171E+00 ; 0.3790996947E-01 ; -0.5279896883E-01 ; -0.3509616877E-01 ; 0.1600620645E-02 ; 0.4098367097E-02 ]
,
[ -0.1219999250E+00 ; 0.2497494925E+00 ; -0.1489636838E+00 ; -0.2314941689E+00 ; 0.4306729538E-01 ; 0.3349988297E+00 ; 0.5170449138E+00 ; 0.5314908191E+00 ; 0.3984867666E+00 ; 0.1134176909E+00 ; -0.8664578083E-01 ; -0.6206801690E-01 ; 0.2555164772E-02 ; 0.6542464174E-02 ]
,
[ 0.6643773700E-01 ; -0.1425076588E+00 ; 0.1150002072E+00 ; 0.1110261320E+00 ; -0.5456405522E-01 ; -0.1403852745E+00 ; -0.9002913124E-01 ; 0.9927107272E-01 ; 0.3252695570E+00 ; 0.5955739150E+00 ; 0.6177120462E+00 ; 0.2102093627E+00 ; -0.1207077310E+00 ; -0.9029360085E-01 ; 0.3498469159E-02 ; 0.8957782052E-02 ]
,
[ -0.2472698999E-01 ; 0.5451727819E-01 ; -0.4913268141E-01 ; -0.3543504434E-01 ; 0.2981657384E-01 ; 0.5443081267E-01 ; 0.3080352023E-01 ; -0.2311127518E-01 ; -0.2663146891E-01 ; -0.3986350788E-01 ; 0.1012109980E+00 ; 0.5654368499E+00 ; 0.7419995864E+00 ; 0.2730680885E+00 ; -0.1373359727E+00 ; -0.1050857718E+00 ; 0.3942767529E-02 ; 0.1009540190E-01 ]
,
[ 0.4060965940E-02 ; -0.9082349871E-02 ; 0.8675066430E-02 ; 0.5263347825E-02 ; -0.5426164014E-02 ; -0.8281419894E-02 ; -0.2928559727E-02 ; 0.7407888982E-02 ; 0.1392100926E-02 ; 0.2865053765E-01 ; 0.5527738563E-01 ; -0.5250033174E-01 ; 0.2337385282E-01 ; 0.5390029057E+00 ; 0.7666684946E+00 ; 0.2878132266E+00 ; -0.1399692009E+00 ; -0.1077108034E+00 ; 0.4008607630E-02 ; 0.1026398458E-01 ]
]
,
8 => Any[
[ 0.9280961080E+00 ; 0.3685482772E+00 ; 0.5169019500E-01 ; -0.1072856966E-01 ; -0.3547264354E-02 ; 0.1881724221E-02 ; 0.2253426546E-03 ; -0.4074549423E-03 ; 0.2542121693E-02 ]
,
[ -0.2880715337E+00 ; 0.6422316271E+00 ; 0.6540393284E+00 ; 0.2730170577E+00 ; 0.4300012686E-01 ; -0.1532507378E-01 ; -0.3095082211E-02 ; 0.3081063250E-02 ; -0.1192117285E-01 ; -0.2510321213E-03 ; 0.1566213637E-02 ]
,
[ 0.1554699662E+00 ; -0.3928508188E+00 ; 0.1655006354E+00 ; 0.6186900270E+00 ; 0.5757733437E+00 ; 0.2765573214E+00 ; 0.4215137090E-01 ; -0.2022141803E-01 ; 0.3170527734E-01 ; 0.2969565406E-02 ; -0.1133537854E-01 ; -0.2472629557E-03 ; 0.1542697449E-02 ]
,
[ -0.9273812077E-01 ; 0.2525480154E+00 ; -0.1903081753E+00 ; -0.2325263183E+00 ; 0.1951123754E+00 ; 0.5747752998E+00 ; 0.6031889613E+00 ; 0.3197987067E+00 ; -0.2318110629E-01 ; -0.2335276594E-01 ; 0.3904642048E-01 ; 0.3358324305E-02 ; -0.1296942483E-01 ; -0.2744734004E-03 ; 0.1712466040E-02 ]
,
[ 0.4096472123E-01 ; -0.1240767851E+00 ; 0.1278114403E+00 ; 0.7138489316E-01 ; -0.1584871307E+00 ; -0.1727866355E+00 ; 0.1143640280E+00 ; 0.5129116028E+00 ; 0.7137757767E+00 ; 0.3414317954E+00 ; -0.3881329228E-01 ; -0.2526829049E-01 ; 0.4424255303E-01 ; 0.3575822034E-02 ; -0.1394488775E-01 ; -0.2875903358E-03 ; 0.1794303867E-02 ]
,
[ -0.1790191282E-01 ; 0.5295161530E-01 ; -0.5704547538E-01 ; -0.1894676577E-01 ; 0.6968256270E-01 ; 0.4547075753E-01 ; -0.7337937479E-01 ; -0.1444748191E+00 ; -0.1409191494E-01 ; 0.4948704363E+00 ; 0.7618895951E+00 ; 0.3591100053E+00 ; -0.4810103061E-01 ; -0.2674149513E-01 ; 0.4783092078E-01 ; 0.3754819526E-02 ; -0.1470868404E-01 ; -0.2997260590E-03 ; 0.1870019815E-02 ]
,
[ 0.4842573875E-02 ; -0.1435818114E-01 ; 0.1564773617E-01 ; 0.3635704163E-02 ; -0.1807616145E-01 ; -0.7951972463E-02 ; 0.2177252927E-01 ; 0.3159516732E-01 ; -0.5393073080E-02 ; -0.1454356948E+00 ; -0.5484912621E-01 ; 0.4834662634E+00 ; 0.7753953496E+00 ; 0.3639514649E+00 ; -0.5141318290E-01 ; -0.2717238099E-01 ; 0.4899043649E-01 ; 0.3804055936E-02 ; -0.1492982098E-01 ; -0.3026846009E-03 ; 0.1888478444E-02 ]
,
[ -0.5886912723E-03 ; 0.1750384584E-02 ; -0.1931250902E-02 ; -0.3621632497E-03 ; 0.2085884179E-02 ; 0.6409242844E-03 ; -0.2736289413E-02 ; -0.3368917714E-02 ; 0.1038171472E-02 ; 0.3193936333E-01 ; 0.6823471193E-02 ; -0.1435014917E+00 ; -0.6098305316E-01 ; 0.4814477634E+00 ; 0.7771077176E+00 ; 0.3644295324E+00 ; -0.5192289139E-01 ; -0.2721779258E-01 ; 0.4913249792E-01 ; 0.3808651449E-02 ; -0.1495173150E-01 ; -0.3029170562E-03 ; 0.1889928755E-02 ]
]
)


text:
@testset "Continuity" begin
    @test Continuity(Vector{Int}) == Discrete
    @test Continuity(UnitRange{Int}) == Continuous
    @test Continuity(1:10) == Continuous
end


text:
#=convert3dto2d
function that flattens a 3d array to a 2d dataframe with colnames
and adds an additional column to indicate which slice of the lost
dimension each row came from, and another additional column to indicate
which row of the lost dimension
inputs: array3d: array to be flattened,
    ind_dim: dimension that gets the column indicator
    y_dim: dimension to be vertical,
    x_dim: dimension to be horizonal
    colnames: set of column names
    indnames: names to indicate levels of flattened dimension
Called by main model for writing model outputs to csv
checks that the lengths of testnames, DRNAMES are appropriate and returns a dataframe
Patricia Levi July 2018
=#
function convert3dto2d(array3d,ind_dim,x_dim,y_dim,colnames,indnames)
    if length(indnames) != size(array3d,ind_dim)
        error("indnames must have same length as ind_dim")
    end
    if length(colnames) != size(array3d,x_dim)+2
        error("colnames must be one longer than x_dim")
    end

    array2d = Array{}
    for j in 1:size(array3d,ind_dim)
        newdf = array3d[j,:,:]
        if j>1
            # print("try to append to array2d \n")
            array2d = vcat(array2d,hcat(newdf,fill(indnames[j],size(array3d,2)),[1:size(array3d,2)...]))
        else
            array2d = hcat(newdf,fill(indnames[j],size(array3d,2)),[1:size(array3d,2)...])
            # print("array2d initialized\n")
        end
    end

    df2d = DataFrame(array2d)
    names!(df2d,[Symbol("$input") for input in colnames])

    return df2d
end

# TESTING
# in this example, ind_dim = 1, y_dim = 2, x_dim = 3
# test = reshape([1:(12*5*2)...],(2,12,5))
# testnames = ["o1","o2","o3","o4","o5","DR_IND"]
# DRNAMES = ["dr1","dr2"]
# xx = convert3dto2d(test,1,3,2,testnames,DRNAMES)
# CSV.write(string(default_fol,"/test2.csv"),xx)

# desired output is a 12x6 dataframe with above colnames
# for j in 1:size(test)[1] #iterate over each generator
#     newdf = test[j,:,:]
#     # dftest[j] = reshape(test[j,:,:], size(test)[2]*size(test)[3])
#     if j>1
#         dftest = vcat(dftest,hcat(newdf,fill(DRNAMES[j],size(test,2))))
#     else
#         dftest = hcat(newdf,fill(DRNAMES[j],size(test,2)))
#     end
# end
# dfend = DataFrame(dftest)
# names!(dfend,[Symbol("$input") for input in testnames])

max_stars_repo_name: pajalevi/Stoch-UC-model-for-DR

max_issues_repo_name: pajalevi/Stoch-UC-model-for-DR

max_forks_repo_name: pajalevi/Stoch-UC-model-for-DR


text:
#= Finding the Greatest Common Divisor of 2 numbers using the Euclidean
Formula to lessen the time complexity.
=#

## Function

function GCD(a, b)
    if (a == 0)
        return b
    end
    if (b == 0)
        return a
    end
    if (a < b)
        a, b = b, a
    end
    ans = a % b
    while (ans != 0)
        a = b
        b = ans
        ans = a % b
    end
    return b
end

## Input

a = readline()
a = parse(Int64, a)
b = readline()
b = parse(Int64, b)

## Calling the function

GCD(a, b)

#=
Sample Test case:
Input:
    a = 1001  b = 2352
Output:
    7
Time complexity: O( log (min(a,b)) )
=#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct LogAnalyticsOutput <: SwaggerModel
    output::Any # spec type: Union{ Nothing, String } # spec name: output

    function LogAnalyticsOutput(;output=nothing)
        o = new()
        validate_property(LogAnalyticsOutput, Symbol("output"), output)
        setfield!(o, Symbol("output"), output)
        o
    end
end # type LogAnalyticsOutput

const _property_map_LogAnalyticsOutput = Dict{Symbol,Symbol}(Symbol("output")=>Symbol("output"))
const _property_types_LogAnalyticsOutput = Dict{Symbol,String}(Symbol("output")=>"String")
Base.propertynames(::Type{ LogAnalyticsOutput }) = collect(keys(_property_map_LogAnalyticsOutput))
Swagger.property_type(::Type{ LogAnalyticsOutput }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_LogAnalyticsOutput[name]))}
Swagger.field_name(::Type{ LogAnalyticsOutput }, property_name::Symbol) =  _property_map_LogAnalyticsOutput[property_name]

function check_required(o::LogAnalyticsOutput)
    true
end

function validate_property(::Type{ LogAnalyticsOutput }, name::Symbol, val)
end


text:
using ITensors
using Random

include(joinpath(ITensors.examples_dir(), "src", "electronk.jl"))
include(joinpath(ITensors.examples_dir(), "src", "hubbard.jl"))

function main(;
  Nx::Int=6,
  Ny::Int=3,
  U::Float64=4.0,
  t::Float64=1.0,
  maxdim::Int=3000,
  conserve_ky=true,
  use_splitblocks=true,
  seed=1234,
)
  Random.seed!(seed)
  @show Threads.nthreads()
  @show ITensors.using_threaded_blocksparse()

  N = Nx * Ny

  sweeps = Sweeps(10)
  maxdims = min.([100, 200, 400, 800, 2000, 3000, maxdim], maxdim)
  setmaxdim!(sweeps, maxdims...)
  setcutoff!(sweeps, 1e-6)
  setnoise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  @show sweeps

  sites = siteinds("ElecK", N; conserve_qns=true, conserve_ky=conserve_ky, modulus_ky=Ny)

  ampo = hubbard(; Nx=Nx, Ny=Ny, t=t, U=U, ky=true)
  H = MPO(ampo, sites)

  # This step makes the MPO more sparse.
  # It generally improves DMRG performance
  # at large bond dimensions but makes DMRG slower at
  # small bond dimensions.
  if use_splitblocks
    H = splitblocks(linkinds, H)
  end

  # Number of structural nonzero elements in a bulk
  # Hamiltonian MPO tensor
  @show nnz(H[end ÷ 2])
  @show nnzblocks(H[end ÷ 2])

  # Create start state
  state = Vector{String}(undef, N)
  for i in 1:N
    x = (i - 1) ÷ Ny
    y = (i - 1) % Ny
    if x % 2 == 0
      if y % 2 == 0
        state[i] = "Up"
      else
        state[i] = "Dn"
      end
    else
      if y % 2 == 0
        state[i] = "Dn"
      else
        state[i] = "Up"
      end
    end
  end

  psi0 = randomMPS(sites, state, 10)

  energy, psi = @time dmrg(H, psi0, sweeps; svd_alg="divide_and_conquer")
  @show Nx, Ny
  @show t, U
  @show flux(psi)
  @show maxlinkdim(psi)
  @show energy
  return energy, H, psi
end


text:
module NWNTools

    include("objects.jl")
    include("intersections.jl")
    include("graph_tools.jl")
    include("circuit_tools.jl")
    include("randomization.jl")
    include("stat_tools.jl")
    include("plotting.jl")

end


text:
using Colors, Test

@testset "Colormatch" begin
    @test colormatch(500) === colormatch(CIE1931_CMF, 500)

    cmfs = (CIE1931_CMF, CIE1964_CMF, CIE1931J_CMF, CIE1931JV_CMF, CIE2006_2_CMF, CIE2006_10_CMF)
    @testset "$cmf" for cmf in cmfs
        @test colormatch(cmf, 350.0) === XYZ(0.0, 0.0, 0.0)
        @test colormatch(cmf, 850.0) === XYZ(0.0, 0.0, 0.0)
        xyz450 = colormatch(cmf, 450.0)
        @test xyz450.y < xyz450.x < xyz450.z
        xyz550 = colormatch(cmf, 550.0)
        @test xyz550.z < xyz550.x < xyz550.y
        xyz600 = colormatch(cmf, 600.0)
        @test xyz600.z < xyz600.y < xyz600.x
        xyz551 = colormatch(cmf, 551.0)
        @test colormatch(cmf, 550.5) ≈ (xyz550 + colormatch(cmf, 551)) / 2
    end
end


text:
## Convective adjustment

𝒟_train  = LearnConvection.Data.data(train, problem; D=16, N=N);

x = 𝒟_train.x[1300]

plott(x) = plot(x, collect(1:16), legend=false)
p(x) = plot!(x, collect(1:16), legend=false)

plott(x)

y = copy(x)
y[12]+=0.1
y[13]+=0.05
y[3]+=0.1
y[4]+=0.1
y[7]+=0.1

myplot = plott(y)
p(y)

png(myplot, "myyy.png")

r = copy(y)
convective_adjust!(r)

myplot2 = plott(r)
png(myplot2, "myyy3.png")

function convective_adjust!(x)
    # remove negative gradients from temperature profile
    for i in length(x)-2:-1:2
        if x[i] > x[i+1]
            if x[i-1] > x[i]; x[i] = x[i+1]
            else; x[i] = (x[i-1]+x[i+1])/2
            end
        end
    end
end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection