text:
"""
`module  JAC.PhotoRecombination`
    ... a submodel of JAC that contains all methods for computing photo-recomnbination, i.e. radiative recombination
        and radiative electron capture properties between some initial and final-state multiplets.
"""
module PhotoRecombination

    using Printf, ..AngularMomentum, ..Basics, ..Continuum, ..ManyElectron, ..PhotoEmission, ..Radial, ..Nuclear, ..TableStrings
    ##x global JAC_counter = 0

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

        + multipoles              ::Array{EmMultipole}           ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}              ... Specifies the gauges to be included into the computations.
        + electronEnergies        ::Array{Float64,1}             ... List of electron energies.
        + ionEnergies             ::Array{Float64,1}             ... List of ion energies.
        + useIonEnergies          ::Bool                         ... Make use of ion energies in [MeV/u] to obtain the electron energies.
        + calcAnisotropy          ::Bool                         ... True, if the overall anisotropy is to be calculated.
        + calcTensors             ::Bool                         ... True, if the statistical tensors are to be calculated and false otherwise.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        electronEnergies          ::Array{Float64,1}
        ionEnergies               ::Array{Float64,1}
        useIonEnergies            ::Bool
        calcAnisotropy            ::Bool
        calcTensors               ::Bool
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
    end

    """
    `PhotoRecombination.Settings()`  ... constructor for the default values of photo recombination line computations
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], Float64[], Float64[], false, false, false, false, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::PhotoRecombination.Settings)`
    #		... prepares a proper printout of the variable settings::PhotoRecombination.Settings.
    function Base.show(io::IO, settings::PhotoRecombination.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "electronEnergies:         $(settings.electronEnergies)  ")
        println(io, "ionEnergies:              $(settings.ionEnergies)  ")
        println(io, "useIonEnergies:           $(settings.useIonEnergies)  ")
        println(io, "calcAnisotropy:           $(settings.calcAnisotropy)  ")
        println(io, "calcTensors:              $(settings.calcTensors)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
    end

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

        + multipole      ::EmMultipole          ... Multipole of the photon emission/absorption.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + kappa          ::Int64                ... partial-wave of the free electron
        + symmetry       ::LevelSymmetry        ... total angular momentum and parity of the scattering state
        + phase          ::Float64              ... phase of the partial wave
        + amplitude      ::Complex{Float64}     ... Rec amplitude associated with the given channel.
    """
    struct  Channel
        multipole        ::EmMultipole
        gauge            ::EmGauge
        kappa            ::Int64
        symmetry         ::LevelSymmetry
        phase            ::Float64
        amplitude        ::Complex{Float64}
    end

    """
    `struct  PhotoRecombination.Line`
        ... defines a type for a Photorecombination line that may include the definition of channels.

        + initialLevel   ::Level                  ... initial-(state) level
        + finalLevel     ::Level                  ... final-(state) level
        + electronEnergy ::Float64                ... Energy of the (incoming free) electron.
        + photonEnergy   ::Float64                ... Energy of the emitted photon.
        + crossSection   ::EmProperty             ... Cross section for this electron capture.
        + hasChannels    ::Bool                   ... Determines whether the individual (sub-) channels are defined in terms of their
                                                      free-electron energy, kappa, multipole, etc., or not.
        + channels       ::Array{PhotoRecombination.Channel,1}    ... List of photorecombination channels of this line.
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        electronEnergy   ::Float64
        photonEnergy     ::Float64
        crossSection     ::EmProperty
        hasChannels      ::Bool
        channels         ::Array{PhotoRecombination.Channel,1}
    end

    """
    `PhotoRecombination.Line()`
        ... constructor for an `empty instance` of a photorecombination line between a specified initial
            and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level)
        Line(initialLevel::Level, finalLevel::Level, 0., 0., EmProperty(0., 0.), false, Channel[])
    end

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

    """
    `PhotoRecombination.amplitude(kind::String, channel::PhotoRecombination.Channel, energy::Float64, finalLevel::Level,
                                  continuumLevel::Level, grid::Radial.Grid)`
        ... to compute the kind = (photorecombination) amplitude
            < alpha_f J_f || O^(photorecombination) || (alpha_i J_i, epsilon kappa) J_t>  due to the electron-photon
            interaction for the given final and continuum level, the partial wave of the outgoing electron as well as
            the given multipole and gauge. A value::ComplexF64 is returned.
    """
    function amplitude(kind::String, channel::PhotoRecombination.Channel, energy::Float64, finalLevel::Level,
                       continuumLevel::Level, grid::Radial.Grid)
        if      kind in [ "photorecombination"]
        #-----------------------------------
            amplitude = PhotoEmission.amplitude("emission", channel.multipole, channel.gauge, energy, finalLevel, continuumLevel, grid)
            amplitude = im^Basics.subshell_l(Subshell(101, channel.kappa)) * exp( im*channel.phase ) * amplitude
        else    error("stop b")
        end

        return( amplitude )
    end

    """
    `PhotoRecombination.computeAmplitudesProperties(line::PhotoRecombination.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                                    settings::PhotoRecombination.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::PhotoRecombination.Line is returned for
            which the amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::PhotoRecombination.Line, nm::Nuclear.Model, grid::Radial.Grid, nrContinuum::Int64,
                                          settings::PhotoRecombination.Settings)
        newChannels = PhotoRecombination.Channel[];;   contSettings = Continuum.Settings(false, nrContinuum);    csC = 0.;    csB = 0.
        for channel in line.channels
            newfLevel = Basics.generateLevelWithSymmetryReducedBasis(line.finalLevel)
            newfLevel = Basics.generateLevelWithExtraSubshell(Subshell(101, channel.kappa), newfLevel)
            newiLevel = Basics.generateLevelWithSymmetryReducedBasis(line.initialLevel)
            cOrbital, phase  = Continuum.generateOrbitalForLevel(line.electronEnergy, Subshell(101, channel.kappa), newiLevel, nm, grid, contSettings)
            newcLevel  = Basics.generateLevelWithExtraElectron(cOrbital, channel.symmetry, newiLevel)
            newChannel = PhotoRecombination.Channel(channel.multipole, channel.gauge, channel.kappa, channel.symmetry, phase, 0.)
            amplitude  = PhotoRecombination.amplitude("photorecombination", channel, line.photonEnergy, newfLevel, newcLevel, grid)
            push!( newChannels, PhotoRecombination.Channel(newChannel.multipole, newChannel.gauge, newChannel.kappa, newChannel.symmetry,
                                                           newChannel.phase, amplitude) )
            if       channel.gauge == Basics.Coulomb     csC = csC + abs(amplitude)^2
            elseif   channel.gauge == Basics.Babushkin   csB = csB + abs(amplitude)^2
            elseif   channel.gauge == Basics.Magnetic    csB = csB + abs(amplitude)^2;   csC = csC + abs(amplitude)^2
            end
        end
        Ji2 = AngularMomentum.twoJ(line.initialLevel.J)
        csFactor     = 8 * pi^3 * Defaults.getDefaults("alpha")^3 * line.photonEnergy / (Ji2 + 1)
        crossSection = EmProperty(csFactor * csC, csFactor * csB)
        newLine = PhotoRecombination.Line( line.initialLevel, line.finalLevel, line.electronEnergy, line.photonEnergy, crossSection, true, newChannels)
        return( newLine )
    end

    """
    `PhotoRecombination.computeAnisotropyParameter(nu::Int64, gauge::EmGauge, line::PhotoRecombination.Line)`
        ... to compute the anisotropy parameter of the emitted photons for the photorecombination of an initially unpolarized ion.
            A value::ComplexF64 is returned.
    """
    function  computeAnisotropyParameter(nu::Int64, gauge::EmGauge, line::PhotoRecombination.Line)
        if  !line.hasChannels   error("No channels are defined for the given PhotoRecombination.line.")         end
        wa = 0.0im;    Ji = line.initialLevel.J;    Jf = line.finalLevel.J
        wn = 0.;
        for  ch in line.channels
            if  gauge != ch.gauge  &&  gauge != Basics.Magnetic   continue    end
            wn = wn + conj(ch.amplitude) * ch.amplitude
        end

        for  cha  in line.channels
            if  gauge != cha.gauge  &&  gauge != Basics.Magnetic   continue    end
            J = cha.symmetry.J;    L = cha.multipole.L;    if  cha.multipole.electric   p = 1   else    p = 0   end
            j = AngularMomentum.kappa_j(cha.kappa);    l = AngularMomentum.kappa_l(cha.kappa)
            #
            for  chp  in line.channels
                if  gauge != chp.gauge  &&  gauge != Basics.Magnetic   continue    end
                Jp = chp.symmetry.J;    Lp = chp.multipole.L;    if  chp.multipole.electric   pp = 1   else    pp = 0   end
                jp = AngularMomentum.kappa_j(chp.kappa);     lp = AngularMomentum.kappa_l(chp.kappa)
                #
                if  1 + (-1)^(L + p + Lp + pp - nu) == 0    continue    end
                wa = wa + 1.0im^(L + p - Lp - pp) * AngularMomentum.phaseFactor([Ji, -1, AngularJ64(1//2), -1, Jf]) *
                                 sqrt( AngularMomentum.bracket([AngularJ64(L), AngularJ64(Lp), l, lp, j, jp, J, Jp]) ) *
                                 AngularMomentum.ClebschGordan( l, AngularM64(0), lp, AngularM64(0),  AngularJ64(nu),  AngularM64(0)) *
                                 AngularMomentum.ClebschGordan( AngularJ64(L), AngularM64(1), AngularJ64(Lp), AngularM64(-1),
                                                                    AngularJ64(nu), AngularM64(0)) *
                                 AngularMomentum.Wigner_6j(J, Jp, AngularJ64(nu), AngularJ64(Lp), AngularJ64(L), Jf) *
                                 AngularMomentum.Wigner_6j(J, Jp, AngularJ64(nu), jp, j, Ji) *
                                 AngularMomentum.Wigner_6j(j, jp, AngularJ64(nu), lp, l, AngularJ64(1//2)) *
                                 conj(cha.amplitude) * chp.amplitude
            end
        end

        value = - 0.5 * wa / wn
        return( value )
    end

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

    """
    `PhotoRecombination.determineChannels(finalLevel::Level, initialLevel::Level, settings::PhotoRecombination.Settings)`
        ... to determine a list of RecChannel for a transitions from the initial to final level and by taking into account
            the particular settings of for this computation; an Array{PhotoRecombination.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::PhotoRecombination.Settings)
        channels = PhotoRecombination.Channel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        for  mp in settings.multipoles
            for  gauge in settings.gauges
                symList = AngularMomentum.allowedMultipoleSymmetries(symf, mp)
                for  symt in symList
                    kappaList = AngularMomentum.allowedKappaSymmetries(symi, symt)
                    for  kappa in kappaList
                        # Include further restrictions if appropriate
                        if     string(mp)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                            push!(channels, PhotoRecombination.Channel(mp, Basics.Coulomb,   kappa, symt, 0., Complex(0.)) )
                        elseif string(mp)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                            push!(channels, PhotoRecombination.Channel(mp, Basics.Babushkin, kappa, symt, 0., Complex(0.)) )
                        elseif string(mp)[1] == 'M'
                            push!(channels, PhotoRecombination.Channel(mp, Basics.Magnetic,  kappa, symt, 0., Complex(0.)) )
                        end
                    end
                end
            end
        end
        return( channels )
    end

    """
    `PhotoRecombination.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PhotoRecombination.Settings)`
        ... to determine a list of PhotoRecombination.Line's for transitions between levels from the initial- and final-state
            multiplets, and by taking into account the particular selections and settings for this computation;
            an Array{PhotoRecombination.Line,1} is returned. Apart from the level specification, all physical properties are set
            to zero during the initialization process.
    """
    function  determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PhotoRecombination.Settings)
        if    settings.selectLines    selectLines   = true;
                       selectedLines = Basics.determineSelectedLines(settings.selectedLines, initialMultiplet, finalMultiplet)
        else                          selectLines   = false
        end

        lines = PhotoRecombination.Line[]
        for  i = 1:length(initialMultiplet.levels)
            for  f = 1:length(finalMultiplet.levels)
                if  selectLines  &&  !((i,f) in selectedLines )    continue   end
                ##x println("PhotoRecombination.determineLines-aa:  i = $i, f = $f")
                #
                if  settings.useIonEnergies
                    error("Not yet implemented: Given ion energies.")
                else
                    for  en in settings.electronEnergies
                        if  en < 0    continue   end
                        # Electron energies are still in 'pre-defined' units; convert to Hartree
                        en_au = Defaults.convertUnits("energy: to atomic", en)
                        omega = en_au + initialMultiplet.levels[i].energy - finalMultiplet.levels[f].energy

                        channels = PhotoRecombination.determineChannels(finalMultiplet.levels[f], initialMultiplet.levels[i], settings)
                        push!( lines, PhotoRecombination.Line(initialMultiplet.levels[i], finalMultiplet.levels[f], en_au, omega,
                                                              EmProperty(0., 0.), true, channels) )
                    end
                end
            end
        end
        return( lines )
    end

    """
    `PhotoRecombination.displayLines(lines::Array{PhotoRecombination.Line,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of all selected
            transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines(lines::Array{PhotoRecombination.Line,1})
        println(" ")
        println("  Selected photorecombination lines:")
        println(" ")
        println("  ", TableStrings.hLine(175))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                                sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                                sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(10, "Energy_if"; na=2);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(12, "Energy e_r"; na=1);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(10, "omega"; na=5);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=4)
        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(sa);    println(sb);    println("  ", TableStrings.hLine(175))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))              * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
            kappaMultipoleSymmetryList = Tuple{Int64,EmMultipole,EmGauge,LevelSymmetry}[]
            for  i in 1:length(line.channels)
                push!( kappaMultipoleSymmetryList, (line.channels[i].kappa, line.channels[i].multipole, line.channels[i].gauge,
                                                    line.channels[i].symmetry) )
            end
            wa = TableStrings.kappaMultipoleSymmetryTupels(85, kappaMultipoleSymmetryList)
            sb = sa * wa[1];    println( sb )
            for  i = 2:length(wa)
                sb = TableStrings.hBlank( length(sa) ) * wa[i];    println( sb )
            end
        end
        println("  ", TableStrings.hLine(175), "\n")
        #
        return( nothing )
    end

    """
    `PhotoRecombination.displayResults(stream::IO, lines::Array{PhotoRecombination.Line,1}, settings::PhotoRecombination.Settings)`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayResults(stream::IO, lines::Array{PhotoRecombination.Line,1}, settings::PhotoRecombination.Settings)
        println(stream, " ")
        println(stream, "  Photorecombination cross sections for beta^2 * gamma^2 = 1:")
        println(stream, "  ")
        println(stream, "  ", TableStrings.hLine(133))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"   ; na=0);                       sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(16, "i--J^P--f"   ; na=3);                       sb = sb * TableStrings.hBlank(19)
        sa = sa * TableStrings.center(12, "i--Energy--f"; na=4)
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(12, "omega"     ; na=4)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(12, "Energy e_p"; na=3)
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(10, "Multipoles"; na=1);                              sb = sb * TableStrings.hBlank(13)
        sa = sa * TableStrings.center(30, "Cou -- Cross section -- Bab"; na=3)
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section") * "          " *
                                              TableStrings.inUnits("cross section"); na=3)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(133))
        #
        for  line in lines
            sa  = " ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                          fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(17, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(17, TableStrings.symmetries_if(isym, fsym); na=3)
            en = line.initialLevel.energy - line.finalLevel.energy
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                  * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "    "
            multipoles = EmMultipole[]
            for  ch in line.channels
                multipoles = push!( multipoles, ch.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
            sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=2)
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Coulomb))     * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Babushkin))   * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(133))
        #
        #
        if  settings.calcAnisotropy
            println(stream, " ")
            println(stream, "  Anisotropy angular parameters beta_nu of the emitted photons for initially unpolarized ions:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(133))
            sa = "  ";   sb = "  "
            sa = sa * TableStrings.center(18, "i-level-f"   ; na=0);                       sb = sb * TableStrings.hBlank(20)
            sa = sa * TableStrings.center(16, "i--J^P--f"   ; na=3);                       sb = sb * TableStrings.hBlank(19)
            sa = sa * TableStrings.center(12, "i--Energy--f"; na=4)
            sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
            sa = sa * TableStrings.center(12, "omega"     ; na=4)
            sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
            sa = sa * TableStrings.center(12, "Energy e_p"; na=3)
            sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=3)
            sa = sa * TableStrings.center(10, "Multipoles"; na=5);                         sb = sb * TableStrings.hBlank(15)
            sa = sa * TableStrings.flushleft(57, "beta's"; na=4)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(133))
            #
            for  line in lines
                sa  = " ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                              fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
                sa = sa * TableStrings.center(17, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
                sa = sa * TableStrings.center(17, TableStrings.symmetries_if(isym, fsym); na=3)
                en = line.initialLevel.energy - line.finalLevel.energy
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", en))                  * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.photonEnergy))   * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.electronEnergy)) * "    "
                multipoles = EmMultipole[]
                for  ch in line.channels
                    multipoles = push!( multipoles, ch.multipole)
                end
                multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
                sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
                be = PhotoRecombination.computeAnisotropyParameter(1, Basics.Coulomb,   line);   sa = sa * @sprintf("%.4e", be.re) * " (1, C);  "
                be = PhotoRecombination.computeAnisotropyParameter(1, Basics.Babushkin, line);   sa = sa * @sprintf("%.4e", be.re) * " (1, B);  "
                println(stream, sa);   sa = TableStrings.hBlank(102)
                be = PhotoRecombination.computeAnisotropyParameter(2, Basics.Coulomb,   line);   sa = sa * @sprintf("%.4e", be.re) * " (2, C);  "
                be = PhotoRecombination.computeAnisotropyParameter(2, Basics.Babushkin, line);   sa = sa * @sprintf("%.4e", be.re) * " (2, B);  "
                println(stream, sa);   sa = TableStrings.hBlank(102)
                be = PhotoRecombination.computeAnisotropyParameter(3, Basics.Coulomb,   line);   sa = sa * @sprintf("%.4e", be.re) * " (3, C);  "
                be = PhotoRecombination.computeAnisotropyParameter(3, Basics.Babushkin, line);   sa = sa * @sprintf("%.4e", be.re) * " (3, B);  "
                println(stream, sa);   sa = TableStrings.hBlank(102)
                be = PhotoRecombination.computeAnisotropyParameter(4, Basics.Coulomb,   line);   sa = sa * @sprintf("%.4e", be.re) * " (4, C);  "
                be = PhotoRecombination.computeAnisotropyParameter(4, Basics.Babushkin, line);   sa = sa * @sprintf("%.4e", be.re) * " (4, B);  "
                println(stream, sa);   sa = TableStrings.hBlank(102)
            end
            println(stream, "  ", TableStrings.hLine(133))
        end
        #
        #
        if  settings.calcTensors
            println(stream, " ")
            println(stream, "  Reduced statistical tensors of the recombined ion ... not yet implemented !!")
            println(stream, " ")
        end
        #
        return( nothing )
    end

end # module


text:
using ExprOptimization, ExprRules
using Test, Random

let
    grammar = @grammar begin
        R = R + R
        R = |(1:2)
    end

    function loss(node::RuleNode, grammar::Grammar)
        Core.eval(node, grammar)
    end

    Random.seed!(0)
    p = CrossEntropy(10, 5, 4, 5)
    res = optimize(p, grammar, :R, loss)
    @test res.expr == 1
    @test Core.eval(res.tree, grammar) == 1
    @test res.loss == 1

    iter = ExpressionIterator(grammar, 2, :R)
    pop = collect(iter)

    losses = Vector{Float64}(undef,length(pop))
    (best_tree, best_loss) = CrossEntropys.evaluate!(p, loss, grammar, pop, losses, pop[1], Inf)
    @test Core.eval(best_tree, grammar) == 1
    @test best_loss == 1

    p = CrossEntropy(10, 5, 4, 5; track_method=CrossEntropys.TopKTracking(3))
    res = optimize(p, grammar, :R, loss)
    res.alg_result[:top_k]
end

ext: jl


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

struct PacketCapturesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Create and start a packet capture on the specified VM.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: packetCaptureName::String (required)
Param: parameters::PacketCapture (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PacketCaptureResult
"""
function _swaggerinternal_packetCapturesCreate(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", PacketCaptureResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/packetCaptures/{packetCaptureName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "packetCaptureName", packetCaptureName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function packetCapturesCreate(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesCreate(_api, resourceGroupName, networkWatcherName, packetCaptureName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function packetCapturesCreate(_api::PacketCapturesApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesCreate(_api, resourceGroupName, networkWatcherName, packetCaptureName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified packet capture session.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: packetCaptureName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_packetCapturesDelete(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/packetCaptures/{packetCaptureName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "packetCaptureName", packetCaptureName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function packetCapturesDelete(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesDelete(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function packetCapturesDelete(_api::PacketCapturesApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesDelete(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets a packet capture session by name.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: packetCaptureName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PacketCaptureResult
"""
function _swaggerinternal_packetCapturesGet(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PacketCaptureResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/packetCaptures/{packetCaptureName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "packetCaptureName", packetCaptureName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function packetCapturesGet(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesGet(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function packetCapturesGet(_api::PacketCapturesApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesGet(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Query the status of a running packet capture session.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: packetCaptureName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PacketCaptureQueryStatusResult
"""
function _swaggerinternal_packetCapturesGetStatus(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", PacketCaptureQueryStatusResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/packetCaptures/{packetCaptureName}/queryStatus", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "packetCaptureName", packetCaptureName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function packetCapturesGetStatus(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesGetStatus(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function packetCapturesGetStatus(_api::PacketCapturesApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesGetStatus(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Lists all packet capture sessions within the specified resource group.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PacketCaptureListResult
"""
function _swaggerinternal_packetCapturesList(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", PacketCaptureListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/packetCaptures", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function packetCapturesList(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesList(_api, resourceGroupName, networkWatcherName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function packetCapturesList(_api::PacketCapturesApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesList(_api, resourceGroupName, networkWatcherName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Stops a specified packet capture session.
Param: resourceGroupName::String (required)
Param: networkWatcherName::String (required)
Param: packetCaptureName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_packetCapturesStop(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkWatchers/{networkWatcherName}/packetCaptures/{packetCaptureName}/stop", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkWatcherName", networkWatcherName)  # type String
    Swagger.set_param(_ctx.path, "packetCaptureName", packetCaptureName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function packetCapturesStop(_api::PacketCapturesApi, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesStop(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function packetCapturesStop(_api::PacketCapturesApi, response_stream::Channel, resourceGroupName::String, networkWatcherName::String, packetCaptureName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_packetCapturesStop(_api, resourceGroupName, networkWatcherName, packetCaptureName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export packetCapturesCreate, packetCapturesDelete, packetCapturesGet, packetCapturesGetStatus, packetCapturesList, packetCapturesStop


text:
using Glob: glob
using Random: shuffle
using Distributed: procs, addprocs
using Distributed: @everywhere

if length(procs()) < Sys.CPU_THREADS-2
    addprocs(Sys.CPU_THREADS-2) ;println("$(length(procs())) cores running.") ;end
@everywhere include("gru_api(seq_to_vec).jl") ; @everywhere include("utils.jl")

const learning_rate = .01
const hm_epochs     = 50

const in_size       = 52
const layers        = [52]
const out_size      = 52

create_model_definitions(layers)

const model = Model(mk_layers(in_size, layers, out_size)...)

const data = [[[randn(1, in_size) for i in 1:rand(1:16)], softmax(randn(1, out_size))] for _ in 1:100] # const data = import_data(glob("class*.txt")[1])

for i in 1:hm_epochs
    print("epoch: $i ")
    train!(model, shuffle(data), learning_rate)
    end
end

save_model(model)

print("Hit enter to continue..")
chomp(readline())


text:
###########
# READING #
###########
function read_field(name::String, default, info, map)
    x = get(info, name, default)
    if default === nothing
        x === nothing && return nothing
    else
        x == default && return default
    end
    x isa String || pkgerror("Expected field `$name` to be a String.")
    return map(x)
end

read_pinned(::Nothing) = false
read_pinned(pinned::Bool) = pinned
read_pinned(::Any) = pkgerror("Expected field `pinned` to be a Boolean.")

function safe_SHA1(sha::String)
    try sha = SHA1(sha)
    catch err
        err isa ArgumentError || rethrow()
        pkgerror("Could not parse `git-tree-sha1` field as a SHA.")
    end
    return sha
end

function safe_uuid(uuid::String)::UUID
    try uuid = UUID(uuid)
    catch err
        err isa ArgumentError || rethrow()
        pkgerror("Could not parse `uuid` field as a UUID.")
    end
    return uuid
end

function safe_bool(bool::String)
    try bool = parse(Bool, bool)
    catch err
        err isa ArgumentError || rethrow()
        pkgerror("Could not parse `pinned` field as a Bool.")
    end
    return bool
end

# note: returns raw version *not* parsed version
function safe_version(version::String)::VersionNumber
    try version = VersionNumber(version)
    catch err
        err isa ArgumentError || rethrow()
        pkgerror("Could not parse `version` as a `VersionNumber`.")
    end
    return version
end

# turn /-paths into \-paths on Windows
function safe_path(path::String)
    if Sys.iswindows() && !isabspath(path)
        path = joinpath(split(path, "/")...)
    end
    return path
end

read_deps(::Nothing) = Dict{String, UUID}()
read_deps(deps) = pkgerror("Expected `deps` field to be either a list or a table.")
function read_deps(deps::AbstractVector)
    ret = String[]
    for dep in deps
        dep isa String || pkgerror("Expected `dep` entry to be a String.")
        push!(ret, dep)
    end
    return ret
end
function read_deps(raw::Dict{String, Any})::Dict{String,UUID}
    deps = Dict{String,UUID}()
    for (name, uuid) in raw
        deps[name] = safe_uuid(uuid)
    end
    return deps
end

struct Stage1
    uuid::UUID
    entry::PackageEntry
    deps::Union{Vector{String}, Dict{String,UUID}}
end

normalize_deps(name, uuid, deps, manifest) = deps
function normalize_deps(name, uuid, deps::Vector{String}, manifest::Dict{String,Vector{Stage1}})
    if length(deps) != length(unique(deps))
        pkgerror("Duplicate entry in `$name=$uuid`'s `deps` field.")
    end
    final = Dict{String,UUID}()
    for dep in deps
        infos = get(manifest, dep, nothing)
        if infos === nothing
            pkgerror("`$name=$uuid` depends on `$dep`, ",
                     "but no such entry exists in the manifest.")
        end
        # should have used dict format instead of vector format
        length(infos) == 1 || pkgerror("Invalid manifest format. ",
                                       "`$name=$uuid`'s dependency on `$dep` is ambiguous.")
        final[dep] = infos[1].uuid
    end
    return final
end

function validate_manifest(julia_version::Union{Nothing,VersionNumber}, manifest_format::VersionNumber, stage1::Dict{String,Vector{Stage1}}, other::Dict{String, Any})
    # expand vector format deps
    for (name, infos) in stage1, info in infos
        info.entry.deps = normalize_deps(name, info.uuid, info.deps, stage1)
    end
    # invariant: all dependencies are now normalized to Dict{String,UUID}
    deps = Dict{UUID, PackageEntry}()
    for (name, infos) in stage1, info in infos
        deps[info.uuid] = info.entry
    end
    # now just verify the graph structure
    for (entry_uuid, entry) in deps, (name, uuid) in entry.deps
        dep_entry = get(deps, uuid, nothing)
        if dep_entry === nothing
            pkgerror("`$(entry.name)=$(entry_uuid)` depends on `$name=$uuid`, ",
                     "but no such entry exists in the manifest.")
        end
        if dep_entry.name != name
            pkgerror("`$(entry.name)=$(entry_uuid)` depends on `$name=$uuid`, ",
                     "but entry with UUID `$uuid` has name `$(dep_entry.name)`.")
        end
    end
    return Manifest(; julia_version, manifest_format, deps, other)
end

function Manifest(raw::Dict, f_or_io::Union{String, IO})::Manifest
    julia_version = haskey(raw, "julia_version") ? VersionNumber(raw["julia_version"]) : nothing
    manifest_format = VersionNumber(raw["manifest_format"])
    if !in(manifest_format.major, 1:2)
        if f_or_io isa IO
            @warn "Unknown Manifest.toml format version detected in streamed manifest. Unexpected behavior may occur" manifest_format
        else
            @warn "Unknown Manifest.toml format version detected in file `$(f_or_io)`. Unexpected behavior may occur" manifest_format maxlog = 1 _id = Symbol(f_or_io)
        end
    end
    stage1 = Dict{String,Vector{Stage1}}()
    if haskey(raw, "deps") # deps field doesn't exist if there are no deps
        for (name, infos) in raw["deps"], info in infos
            entry = PackageEntry()
            entry.name = name
            uuid = nothing
            deps = nothing
            try
                entry.pinned      = read_pinned(get(info, "pinned", nothing))
                uuid              = read_field("uuid",          nothing, info, safe_uuid)::UUID
                entry.version     = read_field("version",       nothing, info, safe_version)
                entry.path        = read_field("path",          nothing, info, safe_path)
                entry.repo.source = read_field("repo-url",      nothing, info, identity)
                entry.repo.rev    = read_field("repo-rev",      nothing, info, identity)
                entry.repo.subdir = read_field("repo-subdir",   nothing, info, identity)
                entry.tree_hash   = read_field("git-tree-sha1", nothing, info, safe_SHA1)
                entry.uuid        = uuid
                deps = read_deps(get(info::Dict, "deps", nothing))
            catch
                # TODO: Should probably not unconditionally log something
                @debug "Could not parse manifest entry for `$name`" f_or_io
                rethrow()
            end
            entry.other = info::Union{Dict,Nothing}
            stage1[name] = push!(get(stage1, name, Stage1[]), Stage1(uuid, entry, deps))
        end
        # by this point, all the fields of the `PackageEntry`s have been type casted
        # but we have *not* verified the _graph_ structure of the manifest
    end
    other = Dict{String, Any}()
    for (k, v) in raw
        if k in ("julia_version", "deps", "manifest_format")
            continue
        end
        other[k] = v
    end
    return validate_manifest(julia_version, manifest_format, stage1, other)
end

function read_manifest(f_or_io::Union{String, IO})
    raw = try
        if f_or_io isa IO
            TOML.parse(read(f_or_io, String))
        else
            isfile(f_or_io) ? parse_toml(f_or_io) : return Manifest()
        end
    catch e
        if e isa TOML.ParserError
            pkgerror("Could not parse manifest: ", sprint(showerror, e))
        end
        rethrow()
    end
    if Base.is_v1_format_manifest(raw)
        raw = convert_v1_format_manifest(raw)
    end
    return Manifest(raw, f_or_io)
end

function convert_v1_format_manifest(old_raw_manifest::Dict)
    new_raw_manifest = Dict{String, Any}(
            "deps" => old_raw_manifest,
            "manifest_format" => "1.0.0" # must be a string here to match raw dict
            # don't set julia_version as it is unknown in old manifests
        )
    return new_raw_manifest
end

###########
# WRITING #
###########
function destructure(manifest::Manifest)::Dict
    function entry!(entry, key, value; default=nothing)
        if value == default
            delete!(entry, key)
        else
            entry[key] = value
        end
    end

    unique_name = Dict{String,Bool}()
    for (uuid, entry) in manifest
        unique_name[entry.name] = !haskey(unique_name, entry.name)
    end

    # maintain the format of the manifest when writing
    if manifest.manifest_format.major == 1
        raw = Dict{String,Vector{Dict{String,Any}}}()
    elseif manifest.manifest_format.major == 2
        raw = Dict{String,Any}()
        if !isnothing(manifest.julia_version) # don't write julia_version if nothing
            raw["julia_version"] = manifest.julia_version
        end
        raw["manifest_format"] = string(manifest.manifest_format.major, ".", manifest.manifest_format.minor)
        raw["deps"] = Dict{String,Vector{Dict{String,Any}}}()
        for (k, v) in manifest.other
            raw[k] = v
        end
    end

    for (uuid, entry) in manifest
        new_entry = something(entry.other, Dict{String,Any}())
        new_entry["uuid"] = string(uuid)
        entry!(new_entry, "version", entry.version)
        entry!(new_entry, "git-tree-sha1", entry.tree_hash)
        entry!(new_entry, "pinned", entry.pinned; default=false)
        path = entry.path
        if path !== nothing && Sys.iswindows() && !isabspath(path)
            path = join(splitpath(path), "/")
        end
        entry!(new_entry, "path", path)
        repo_source = entry.repo.source
        if repo_source !== nothing && Sys.iswindows() && !isabspath(repo_source) && !isurl(repo_source)
            repo_source = join(splitpath(repo_source), "/")
        end
        entry!(new_entry, "repo-url", repo_source)
        entry!(new_entry, "repo-rev", entry.repo.rev)
        entry!(new_entry, "repo-subdir", entry.repo.subdir)
        if isempty(entry.deps)
            delete!(new_entry, "deps")
        else
            if all(dep -> unique_name[first(dep)], entry.deps)
                new_entry["deps"] = sort(collect(keys(entry.deps)))
            else
                new_entry["deps"] = Dict{String,String}()
                for (name, uuid) in entry.deps
                    new_entry["deps"][name] = string(uuid)
                end
            end
        end
        if manifest.manifest_format.major == 1
            push!(get!(raw, entry.name, Dict{String,Any}[]), new_entry)
        elseif manifest.manifest_format.major == 2
            push!(get!(raw["deps"], entry.name, Dict{String,Any}[]), new_entry)
        end
    end
    return raw
end

function write_manifest(env::EnvCache)
    mkpath(dirname(env.manifest_file))
    write_manifest(env.manifest, env.manifest_file)
end
function write_manifest(manifest::Manifest, manifest_file::AbstractString)
    if manifest.manifest_format.major == 1
        @warn """The active manifest file at `$(manifest_file)` has an old format that is being maintained.
            To update to the new format, which is supported by Julia versions ≥ 1.6.2, run `Pkg.upgrade_manifest()` which will upgrade the format without re-resolving.
            To then record the julia version re-resolve with `Pkg.resolve()` and if there are resolve conflicts consider `Pkg.update()`.""" maxlog = 1 _id = Symbol(manifest_file)
    end
    return write_manifest(destructure(manifest), manifest_file)
end
function write_manifest(io::IO, manifest::Manifest)
    return write_manifest(io, destructure(manifest))
end
function write_manifest(io::IO, raw_manifest::Dict)
    print(io, "# This file is machine-generated - editing it directly is not advised\n\n")
    TOML.print(io, raw_manifest, sorted=true) do x
        (typeof(x) in [String, Nothing, UUID, SHA1, VersionNumber]) && return string(x)
        error("unhandled type `$(typeof(x))`")
    end
    return nothing
end
function write_manifest(raw_manifest::Dict, manifest_file::AbstractString)
    str = sprint(write_manifest, raw_manifest)
    write(manifest_file, str)
end

############
# METADATA #
############

function check_warn_manifest_julia_version_compat(manifest::Manifest, manifest_file::String)
    isempty(manifest.deps) && return
    if manifest.manifest_format < v"2"
        @warn """The active manifest file is an older format with no julia version entry. Dependencies may have \
        been resolved with a different julia version.""" maxlog = 1 _file = manifest_file _line = 0 _module = nothing
        return
    end
    v = manifest.julia_version
    if v === nothing
        @warn """The active manifest file is missing a julia version entry. Dependencies may have \
        been resolved with a different julia version.""" maxlog = 1 _file = manifest_file _line = 0 _module = nothing
        return
    end
    if Base.thisminor(v) != Base.thisminor(VERSION)
        @warn """The active manifest file has dependencies that were resolved with a different julia \
        version ($(manifest.julia_version)). Unexpected behavior may occur.""" maxlog = 1 _file = manifest_file _line = 0 _module = nothing
    end
end


text:
using Clang.Generators
using MacroTools
using GEOS_jll
using JuliaFormatter: format

"Functions that return a Cstring are wrapped in unsafe_string to return a String"
function rewrite(ex::Expr)
    if @capture(ex, function fname_(fargs__)
        @ccall lib_.cname_(cargs__)::rettype_
    end)

        # bind the ccall such that we can easily wrap it
        cc = :(@ccall $lib.$cname($(cargs...))::$rettype)

        cc′ = if rettype == :Cstring
            :(unsafe_string($cc))
        else
            cc
        end

        # stitch the modified function expression back together
        f = :(function $fname($(fargs...))
            $cc′
        end) |> prettify

        return f
    end
    return ex
end

function rewrite!(dag::ExprDAG)
    for n in dag.nodes
        map!(rewrite, n.exprs, n.exprs)
    end
end

cd(@__DIR__)

include_dir = joinpath(GEOS_jll.artifact_dir, "include") |> normpath

# wrapper generator options
options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")

header_dir = include_dir
headers = [joinpath(header_dir, "geos_c.h")]

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx, BUILDSTAGE_NO_PRINTING)
rewrite!(ctx.dag)
build!(ctx, BUILDSTAGE_PRINTING_ONLY)

# run JuliaFormatter on the whole package
format(joinpath(@__DIR__, ".."))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Torch"
version = v"1.4.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/dhairyagandhi96/Torch.jl.git", "85bd08d39e7fba29ec4a643f60dd006ed8be8ede"),
    ArchiveSource("https://download.pytorch.org/libtorch/cu101/libtorch-cxx11-abi-shared-with-deps-1.4.0.zip", "f214bfde532877aa5d4e0803e51a28fa8edd97b6a44b6615f75a70352b6b542e"),
    ArchiveSource("https://github.com/JuliaGPU/CUDABuilder/releases/download/v0.3.0/CUDNN+CUDA10.1.v7.6.5.x86_64-linux-gnu.tar.gz", "79de5b5085a33bc144b87028e998a1d295a15c3424d6d45b25defe500f616974", unpack_target = "cudnn"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir

mv cudnn $prefix
mv libtorch/share/* $prefix/share/
mv libtorch/lib/* $prefix/lib/
rm -r libtorch/lib
rm -r libtorch/share
mv libtorch/* $prefix
rm -r libtorch

mkdir -p /usr/local/cuda/lib64
cd /usr/local/cuda/lib64
ln -s ${prefix}/cuda/lib64/libcudart.so libcudart.so
ln -s ${prefix}/cuda/lib64/libnvToolsExt.so libnvToolsExt.so

cd $WORKSPACE/srcdir/Torch.jl/build
mkdir build && cd build
cmake -DCMAKE_PREFIX_PATH=$prefix -DTorch_DIR=$prefix/share/cmake/Torch -DCUDA_TOOLKIT_ROOT_DIR=$prefix/cuda ..
cmake --build .

mkdir -p "${libdir}"
cp -r $WORKSPACE/srcdir/Torch.jl/build/build/*.${dlext} "${libdir}"
rm -rf $prefix/cuda
install_license ${WORKSPACE}/srcdir/Torch.jl/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc", cxxstring_abi = "cxx11"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libdoeye_caml", :libdoeye_caml, dont_dlopen = true),
    LibraryProduct("libtorch", :libtorch, dont_dlopen = true),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="CUDA_full_jll", version=v"10.1.243")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll")),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7.1.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module Day16

import Base.show
import ..data_dir # from parent module
input = read(joinpath(data_dir, "day16"), String)

export part1, part2

## Functions for defining packets from inputs

# To avoid integer overflows, we convert each hex to a 4-character bitstring
# individually, rather than converting hex to bits all at once
function bits(hex::AbstractString)
    if length(hex) == 1
        return (last(bitstring(parse(Int64, hex; base = 16)), 4))
    end
    return join(bits.(split(hex, "")))
end

decimal(bitstring::AbstractString) = parse(Int64, bitstring; base = 2)
version(packet_bitstring::AbstractString) = packet_bitstring[1:3] |> decimal
type_id(packet_bitstring::AbstractString) = packet_bitstring[4:6] |> decimal

## Packet definitions

abstract type Packet end

struct LiteralPacket <: Packet
    decoded::Int64
    length::Int64
    version::Int64
    type_id::Int64
end

struct OperatorPacket <: Packet
    subpackets::Vector{Packet}
    length::Int64
    version::Int64
    type_id::Int64
end

const type_id_function = Dict(
    # Need to make it so that all functions work on vectors
    0 => sum,
    1 => prod,
    2 => minimum,
    3 => maximum,
    # omit 4, which is reserved for literals
    # the following only ever contain two packets
    5 => x -> Int(>(x...)),
    6 => x -> Int(<(x...)),
    7 => x -> Int(==(x...)),
)

## Packet show functions

const type_id_description = Dict(
    0 => "+",
    1 => "×",
    2 => "min",
    3 => "max",
    5 => ">",
    6 => "<",
    7 => "=",
)

function show(io::IO, literal_packet::LiteralPacket; indent = 0)
    print(io, repeat("  ", indent))
    println(io, "$(literal_packet.decoded) <v$(literal_packet.version)>")
end

function show(io::IO, operator_packet::OperatorPacket; indent = 0)
    print(io, repeat("  ", indent))
    println(io, "$(type_id_description[operator_packet.type_id]) <v$(operator_packet.version)>")
    for subpacket in operator_packet.subpackets
        show(io, subpacket; indent = indent + 1)
    end
end

## Packet extractions

function extract_next_packet(packet_bitstring::AbstractString)
    if type_id(packet_bitstring) == 4
        return extract_literal_packet(packet_bitstring)
    else
        return extract_operator_packet(packet_bitstring)
    end
end

function extract_literal_packet(packet_bitstring::AbstractString)
    last_bit = 11 # header + one group
    while packet_bitstring[last_bit-4] != '0'
        last_bit += 5
    end
    content = packet_bitstring[7:last_bit]
    groups = [content[i:min(i + 4, end)] for i = 1:5:length(content)]
    decoded = [group[2:end] for group in groups] |> join |> decimal
    packet = LiteralPacket(
        decoded,
        last_bit,
        version(packet_bitstring[1:last_bit]),
        4
    )
    return packet, packet_bitstring[last_bit+1:end]
end

function extract_operator_packet(packet_bitstring::AbstractString)
    length_type_id = packet_bitstring[7]
    if length_type_id == '0'
        header_length = 3 + 3 + 1 + 15
        content_start = header_length + 1
        next_subpackets_length = decimal(packet_bitstring[8:content_start-1]) # next 15 bits

        subpackets = Vector{Packet}()
        remainder = packet_bitstring[content_start:end]
        while true
            packet, remainder = extract_next_packet(remainder)
            push!(subpackets, packet)
            subpackets_length = sum(subpacket.length for subpacket in subpackets)
            if subpackets_length == next_subpackets_length
                break
            end
        end
    else # length_type_id == '1'
        header_length = 3 + 3 + 1 + 11
        content_start = header_length + 1
        number_subpackets = decimal(packet_bitstring[8:content_start-1]) # next 11 bits

        subpackets = Vector{Packet}()
        remainder = packet_bitstring[content_start:end]
        while true
            packet, remainder = extract_next_packet(remainder)
            push!(subpackets, packet)
            if length(subpackets) == number_subpackets
                break
            end
        end
        subpackets_length = sum(subpacket.length for subpacket in subpackets)
    end

    packet = OperatorPacket(
        subpackets,
        header_length + subpackets_length,
        version(packet_bitstring),
        type_id(packet_bitstring)
    )

    return packet, remainder
end

# Functions for part 1

sum_versions(packet::LiteralPacket) = packet.version
sum_versions(packet::OperatorPacket) = packet.version + sum(sum_versions(subpacket) for subpacket in packet.subpackets)

function part1(input = input)
    packet_bitstring = bits(input)
    outermost_packet, _ = extract_next_packet(packet_bitstring) # _ is remainder
    sum_versions(outermost_packet)
end

# Functions for part 2

decode(packet::LiteralPacket) = packet.decoded

function decode(packet::OperatorPacket)
    operation = type_id_function[packet.type_id]
    operation([decode(subpacket) for subpacket in packet.subpackets])
end

function part2(input = input)
    packet_bitstring = bits(input)
    outermost_packet, _ = extract_next_packet(packet_bitstring) # _ is remainder
    decode(outermost_packet)
end

end #module


text:
using SortFilters
using Test
using Random
using Statistics

len = 100
Random.seed!(0)
T = Int
x = rand(T(-10):T(10), len)
window = 9

mx_ref = similar(x)
mn_ref = similar(x)
med_ref = similar(x, Float64)
p20_ref = similar(x, Float64)
p25_p75_ref = similar(x, NTuple{2,Float64})

for i = 1:len
    chunk = x[max(1, i-window+1) : i]
    mx_ref[i] = maximum(chunk)
    mn_ref[i] = minimum(chunk)
    med_ref[i] = median(chunk)
    p20_ref[i] = quantile(chunk, 0.20)
    p25_p75_ref[i] = quantile(chunk, (0.25, 0.75))
end

@testset "movsort" begin

mx = movsort(x, window, 1.0)
mx_ = movsort(x, window, 100)

mn = movsort(x, window, 0.0)
mn_ = movsort(x, window, 0)

med = movsort(x, window, 0.5)
med_ = movsort(x, window, 50)

p25_p75 = movsort(x, window, (0.25, 0.75))
p25_p75_ = movsort(x, window, (25, 75))

@test mx == mx_ == mx_ref
@test mn == mn_ == mn_ref
@test med[window:end] == med_[window:end] == med_ref[window:end]
@test p25_p75[window:end] == p25_p75_[window:end] == p25_p75_ref[window:end]

p20 = movsort(x, window, 0.20)
@test p20[window:end] == p20_ref[window:end] # should match nearest index quantiles

end

@testset "movsort!" begin

mx 		 = similar(mx_ref, Float64)
mx_ 	 = similar(mx_ref)
mn 		 = similar(mn_ref, Float64)
mn_ 	 = similar(mn_ref)
med 	 = similar(med_ref, Float64)
med_ 	 = similar(med_ref)
p25_p75  = similar(p25_p75_ref, NTuple{2,Float64})
p25_p75_ = similar(p25_p75_ref)
p20      = similar(p20_ref, Float64)

movsort!(mx, x, window, 1.0)
movsort!(mx_, x, window, 100)

movsort!(mn, x, window, 0.0)
movsort!(mn_, x, window, 0)

movsort!(med, x, window, 0.5)
movsort!(med_, x, window, 50)

movsort!(p25_p75, x, window, (0.25, 0.75))
movsort!(p25_p75_, x, window, (25, 75))

@test mx == mx_ == mx_ref
@test mn == mn_ == mn_ref
@test med[window:end] == med_[window:end] == med_ref[window:end]
@test p25_p75[window:end] == p25_p75_[window:end] == p25_p75_ref[window:end]

movsort!(p20, x, window, 0.20)
@test p20[window:end] == p20_ref[window:end] # should not match nearest index quantiles

end


text:
# XpressMathProgInterface
# Standardized MILP interface

export XpressSolver

type XpressMathProgModel <: AbstractLinearQuadraticModel
    inner::Model

    lazycb
    cutcb
    heuristiccb
    infocb

    options
end
function XpressMathProgModel(;options...)

   m = XpressMathProgModel(Model(), nothing, nothing, nothing, nothing, options)

   for (name,value) in options
       setparam!(m.inner, XPRS_CONTROLS_DICT[name], value)
   end
   return m
end

immutable XpressSolver <: AbstractMathProgSolver
    options
end
XpressSolver(;kwargs...) = XpressSolver(kwargs)
LinearQuadraticModel(s::XpressSolver) = XpressMathProgModel(;s.options...)
ConicModel(s::XpressSolver) = LPQPtoConicBridge(LinearQuadraticModel(s))

supportedcones(m::XpressSolver) = [:Free]#,:Zero,:NonNeg,:NonPos,:SOC,:SOCRotated]

loadproblem!(m::XpressMathProgModel, filename::AbstractString) = read_model(m.inner, filename)

function updatemodel!(m::XpressMathProgModel)
end
# = Base.warn_once("Model update not necessary for Xpress.")
function update_model!(m::XpressMathProgModel)
end
# = Base.warn_once("Model update not necessary for Xpress.")

function loadproblem!(m::XpressMathProgModel, A, collb, colub, obj, rowlb, rowub, sense)
  # throw away old model
  #env = m.inner.env
  #m.inner.finalize_env = false
  free_model(m.inner)

  m.inner = Model()

  for (name,value) in m.options
       setparam!(m.inner, XPRS_CONTROLS_DICT[name], value)
  end

  add_cvars!(m.inner, float(obj), float(collb), float(colub))

  neginf = typemin(eltype(rowlb))
  posinf = typemax(eltype(rowub))

  # check if we have any range constraints
  # to properly support these, we will need to keep track of the
  # slack variables automatically added by gurobi.
  rangeconstrs = any((rowlb .!= rowub) & (rowlb .> neginf) & (rowub .< posinf))
  if rangeconstrs
      warn("Julia Xpress interface doesn't properly support range (two-sided) constraints.")
      add_rangeconstrs!(m.inner, float(A), float(rowlb), float(rowub))
  else
      b = Array(Float64,length(rowlb))
      senses = Array(Cchar,length(rowlb))
      for i in 1:length(rowlb)
          if rowlb[i] == rowub[i]
              senses[i] = XPRS_EQ#'='
              b[i] = rowlb[i]
          elseif rowlb[i] > neginf
              senses[i] = XPRS_GEQ#'>'
              b[i] = rowlb[i]
          else
              @assert rowub[i] < posinf
              senses[i] = XPRS_LEQ#'<'
              b[i] = rowub[i]
          end
      end
      add_constrs!(m.inner, float(A), senses, b)
  end

  setsense!(m,sense)
end

writeproblem(m::XpressMathProgModel, filename::AbstractString) = write_model(m.inner, filename)

getvarLB(m::XpressMathProgModel)     = get_lb(m.inner)
setvarLB!(m::XpressMathProgModel, l) = set_lb!(m.inner,l)

getvarUB(m::XpressMathProgModel)     = get_ub(m.inner)
setvarUB!(m::XpressMathProgModel, u) = set_ub!(m.inner,u)

function getconstrLB(m::XpressMathProgModel)
    sense = get_rowtype(m.inner)
    ret   = get_rhs(m.inner)
    for i = 1:num_constrs(m.inner)
        if sense[i] == XPRS_GEQ || sense[i] == XPRS_EQ
            # Do nothing
        else
            # LEQ constraint, so LB is -Inf
            ret[i] = -Inf
        end
     end
     return ret
end

function getconstrUB(m::XpressMathProgModel)
    sense = get_rowtype(m.inner)
    ret   = get_rhs(m.inner)
    for i = 1:num_constrs(m.inner)
        if sense[i] == XPRS_LEQ || sense[i] == XPRS_EQ
            # Do nothing
        else
            # GEQ constraint, so UB is +Inf
            ret[i] = +Inf
        end
    end
    return ret
end

setconstrLB!(m::XpressMathProgModel, lb) = set_constrLB!(m.inner,lb)
setconstrUB!(m::XpressMathProgModel, ub) = set_constrUB!(m.inner,ub)

getobj(m::XpressMathProgModel)     = get_obj(m.inner)
setobj!(m::XpressMathProgModel, c) = set_obj!(m.inner,c)

function addvar!(m::XpressMathProgModel, constridx, constrcoef, l, u, objcoef)

    add_var!(m.inner, length(constridx), constridx, float(constrcoef), float(objcoef), float(l), float(u))
end
function addvar!(m::XpressMathProgModel, l, u, objcoef)

    add_cvar!(m.inner, float(objcoef), float(l), float(u))
end

function addconstr!(m::XpressMathProgModel, varidx, coef, lb, ub)

    if lb == -Inf
        # <= constraint
        add_constr!(m.inner, varidx, coef, XPRS_LEQ, ub)
    elseif ub == +Inf
        # >= constraint
        add_constr!(m.inner, varidx, coef, XPRS_GEQ, lb)
    elseif lb == ub
        # == constraint
        add_constr!(m.inner, varidx, coef, XPRS_EQ, lb)
    else
        # Range constraint
        error("Adding range constraints not supported yet.")
    end
end

getconstrmatrix(m::XpressMathProgModel) = get_constrmatrix(m.inner)

function setsense!(m::XpressMathProgModel, sense)
  if sense == :Min
    set_sense!(m.inner, :minimize)
  elseif sense == :Max
    set_sense!(m.inner, :maximize)
  else
    error("Unrecognized objective sense $sense")
  end
end

function getsense(m::XpressMathProgModel)
  if model_sense(m.inner) == :maximize
    return :Max
  else
    return :Min
  end
end

numvar(m::XpressMathProgModel)    = num_vars(m.inner)
numconstr(m::XpressMathProgModel) = num_constrs(m.inner)
numlinconstr(m::XpressMathProgModel) = num_constrs(m.inner) - num_qconstrs(m.inner)
numquadconstr(m::XpressMathProgModel) = num_qconstrs(m.inner)

function optimize!(m::XpressMathProgModel)
    # set callbacks if present
    #if m.lazycb != nothing || m.cutcb != nothing || m.heuristiccb != nothing || m.infocb != nothing
    #    setmathprogcallback!(m)
    #end
    #if m.lazycb != nothing
    #  setparam!(m.inner, "LazyConstraints", 1)
    #end
    optimize(m.inner)
end

#fix status
# add mip outputs
function status(m::XpressMathProgModel)
  is_mip(m.inner) ? s = get_mip_status(m.inner) : s = get_lp_status(m.inner)

  if s == :optimal
    return :Optimal
  elseif s == :infeasible
    return :Infeasible
  elseif s == :unbounded
    return :Unbounded
  elseif s == :cutoff
    return :CutOff
  elseif s == :unfinished
    return :Unfinished
  elseif s == :cutoff_in_dual
    return :CutOff_In_Dual
  elseif s == :unsolved
    return :Unsolved
  elseif s == :nonconvex
    return :Nonconvex
  #mip
  elseif s == :mip_optimal
    return :Optimal
  elseif s == :mip_infeasible
    return :Infeasible
  elseif s == :mip_unbounded
    return :Unbounded
  elseif s == :mip_not_loaded
    return :Unfinished
  elseif s == :mip_lp_not_optimal
    return :Unfinished
  elseif s == :mip_lp_optimal
    return :Unfinished
  elseif s == :mip_no_sol_found
    return :Unfinished
  elseif s == :mip_suboptimal
    return :Unfinished
#old
  elseif s == :iteration_limit || s == :node_limit || s == :time_limit || s == :solution_limit
    return :UserLimit
  elseif s == :numeric
    return :Numeric
  elseif s == :suboptimal
    return :Suboptimal # not very useful status
  elseif s == :interrupted # ended by user?
    return :UserLimit
  else
    error("Unrecognized solution status: $s")
  end
end

getobjval(m::XpressMathProgModel)   = get_objval(m.inner)
getobjbound(m::XpressMathProgModel) = get_bestbound(m.inner) #get_objbound(m.inner)
getsolution(m::XpressMathProgModel) = get_solution(m.inner)

function getconstrsolution(m::XpressMathProgModel)
    rows = num_constrs(m.inner)
    sense = get_rowtype(m.inner)
    rhs   = get_rhs(m.inner)
    slack = get_slack(m.inner)
    ret   = zeros(rows)
    for i = 1:rows
        if sense[i] == XPRS_EQ
            # Must be equal to RHS if feasible
            ret[i] = rhs[i]
        else
            # Slack variable is non-negative for <= constraints, and non-positive for >= constraints
            ret[i] = rhs[i] - slack[i]
        end
    end
    return ret
end

function getreducedcosts(m::XpressMathProgModel)
    if is_qcp(m.inner) #&& get_int_param(m.inner.env, "QCPDual") == 0
        return fill(NaN, num_vars(m.inner))
    else
        return get_reducedcost(m.inner)
    end
end

function getconstrduals(m::XpressMathProgModel)
    if is_qcp(m.inner) #&& get_int_param(m.inner.env, "QCPDual") == 0
        return fill(NaN, num_constrs(m.inner))
    else
        return get_dual(m.inner)
    end
end
#=
function getquadconstrduals(m::XpressMathProgModel)
    if is_qcp(m.inner) #&& get_int_param(m.inner.env, "QCPDual") == 0
        return fill(NaN, num_qconstrs(m.inner))
    else
        return get_dual(m.inner)
    end
end
=#

getinfeasibilityray(m::XpressMathProgModel) = getdualray(m.inner)
getunboundedray(m::XpressMathProgModel) = getprimalray(m.inner)

getbasis(m::XpressMathProgModel) = get_basis(m.inner)

getrawsolver(m::XpressMathProgModel) = m.inner

const var_type_map = Dict(
  'C' => :Cont,
  'B' => :Bin,
  'I' => :Int,
  'S' => :SemiCont,
  'R' => :SemiInt
)

const rev_var_type_map = Dict(
  :Cont => 'C',
  :Bin => 'B',
  :Int => 'I',
  :SemiCont => 'S',
  :SemiInt => 'R'
)

function setvartype!(m::XpressMathProgModel, vartype::Vector{Symbol})
    # do this to make sure we deal with new columns
    ind = collect( 1:num_vars(m.inner) )

    nvartype = map(x->Cchar(rev_var_type_map[x]), vartype)

    chgcoltypes!(m.inner, ind, nvartype )
end

function getvartype(m::XpressMathProgModel)
    ret = get_coltype(m.inner)
    map(x->var_type_map[x], ret)
end
#=
function setwarmstart!(m::XpressMathProgModel, v)
    for j = 1:length(v)
        if isnan(v[j])
            v[j] = 1e101  # GRB_UNDEFINED
        end
    end
    set_dblattrarray!(m.inner, "Start", 1, num_vars(m.inner), v)
end
=#
addsos1!(m::XpressMathProgModel, idx, weight) = add_sos!(m.inner, :SOS1, idx, weight)
addsos2!(m::XpressMathProgModel, idx, weight) = add_sos!(m.inner, :SOS2, idx, weight)

# QCQP

function setquadobj!(m::XpressMathProgModel, rowidx, colidx, quadval)
    delq!(m.inner)

    # xpress only accept one input per matrix slot
    k = ones(Bool,length(rowidx)) #replicates holder (they are falses)
    for i in 1:length(rowidx), j in (i+1):length(rowidx)
      if rowidx[i] == rowidx[j] && colidx[i] == colidx[j]
        quadval[i] += quadval[j]
        quadval[j] = 0.0
        rowidx[j] = -1
        colidx[j] = -1
        k[j] = false
      end
    end
    scaledvals = similar(quadval)
    for i in 1:length(rowidx)
      if rowidx[i] == colidx[i]
        # rescale from matrix format to "terms" format
        scaledvals[i] = quadval[i] #/ 2
      else
        scaledvals[i] = quadval[i]
      end
    end
    add_qpterms!(m.inner, rowidx[k], colidx[k], scaledvals[k])

end

function addquadconstr!(m::XpressMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)

    add_qconstr!(m.inner, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
end

#getsolvetime(m::XpressMathProgModel) = get_runtime(m.inner)
getnodecount(m::XpressMathProgModel) = get_node_count(m.inner)

# Callbacks

#=
setlazycallback!(m::XpressMathProgModel,f) = (m.lazycb = f)
setcutcallback!(m::XpressMathProgModel,f) = (m.cutcb = f)
setheuristiccallback!(m::XpressMathProgModel,f) = (m.heuristiccb = f)
setinfocallback!(m::XpressMathProgModel,f) = (m.infocb = f)

type GurobiCallbackData <: MathProgCallbackData
    cbdata::CallbackData
    state::Symbol
    where::Cint
    sol::Vector{Float64}  # Used for heuristic callbacks
#    model::XpressMathProgModel # not needed?
end

function cbgetmipsolution(d::GurobiCallbackData)
    @assert d.state == :MIPSol
    return cbget_mipsol_sol(d.cbdata, d.where)
end

function cbgetmipsolution(d::GurobiCallbackData,output)
    @assert d.state == :MIPSol
    return cbget_mipsol_sol(d.cbdata, d.where, output)
end

function cbgetlpsolution(d::GurobiCallbackData)
    @assert d.state == :MIPNode
    return cbget_mipnode_rel(d.cbdata, d.where)
end

function cbgetlpsolution(d::GurobiCallbackData, output)
    @assert d.state == :MIPNode
    return cbget_mipnode_rel(d.cbdata, d.where, output)
end

# TODO: macro for these getters?
function cbgetobj(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_objbst(d.cbdata, d.where)
    elseif d.state == :MIPInfo
        return cbget_mip_objbst(d.cbdata, d.where)
    elseif d.state == :MIPSol
        error("Gurobi does not implement cbgetobj when state == MIPSol")
        # https://groups.google.com/forum/#!topic/gurobi/Az_X6Ag-y6k
        # https://github.com/JuliaOpt/MathProgBase.jl/issues/23
        # A complex workaround is possible but hasn't been implemented.
        # return cbget_mipsol_objbst(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

function cbgetbestbound(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_objbnd(d.cbdata, d.where)
    elseif d.state == :MIPSol
        return cbget_mipsol_objbnd(d.cbdata, d.where)
    elseif d.state == :MIPInfo
        return cbget_mip_objbnd(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

function cbgetexplorednodes(d::GurobiCallbackData)
    if d.state == :MIPNode
        return cbget_mipnode_nodcnt(d.cbdata, d.where)
    elseif d.state == :MIPSol
        return cbget_mipsol_nodcnt(d.cbdata, d.where)
    elseif d.state == :MIPInfo
        return cbget_mip_nodcnt(d.cbdata, d.where)
    else
        error("Unrecognized callback state $(d.state)")
    end
end

# returns :MIPNode :MIPSol :Other
cbgetstate(d::GurobiCallbackData) = d.state

function cbaddcut!(d::GurobiCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode
    cbcut(d.cbdata, convert(Vector{Cint}, varidx), float(varcoef), sense, float(rhs))
end

function cbaddlazy!(d::GurobiCallbackData,varidx,varcoef,sense,rhs)
    @assert d.state == :MIPNode || d.state == :MIPSol
    cblazy(d.cbdata, convert(Vector{Cint}, varidx), float(varcoef), sense, float(rhs))
end

function cbaddsolution!(d::GurobiCallbackData)
    # Gurobi doesn't support adding solutions on MIPSol.
    # TODO: support this anyway
    @assert d.state == :MIPNode
    cbsolution(d.cbdata, d.sol)
    # "Wipe" solution back to GRB_UNDEFINIED
    for i in 1:length(d.sol)
        d.sol[i] = 1e101  # GRB_UNDEFINED
    end
end

function cbsetsolutionvalue!(d::GurobiCallbackData,varidx,value)
    d.sol[varidx] = value
end

# breaking abstraction, define our low-level callback to eliminatate
# a level of indirection

function mastercallback(ptr_model::Ptr{Void}, cbdata::Ptr{Void}, where::Cint, userdata::Ptr{Void})

    model = unsafe_pointer_to_objref(userdata)::XpressMathProgModel
    grbrawcb = CallbackData(cbdata,model.inner)
    if where == CB_MIPSOL
        state = :MIPSol
        grbcb = GurobiCallbackData(grbrawcb, state, where, [0.0])
        if model.lazycb != nothing
            ret = model.lazycb(grbcb)
            if ret == :Exit
                terminate(model.inner)
            end
        end
    elseif where == CB_MIPNODE
        state = :MIPNode
        # skip callback if node is reported to be cut off or infeasible --
        # nothing to do.
        # TODO: users may want this information
        status = cbget_mipnode_status(grbrawcb, where)
        if status != 2
            return convert(Cint,0)
        end
        grbcb = GurobiCallbackData(grbrawcb, state, where, [0.0])
        if model.cutcb != nothing
            ret = model.cutcb(grbcb)
            if ret == :Exit
                terminate(model.inner)
            end
        end
        if model.heuristiccb != nothing
            grbcb.sol = fill(1e101, numvar(model))  # GRB_UNDEFINED
            ret = model.heuristiccb(grbcb)
            if ret == :Exit
                terminate(model.inner)
            end
        end
        if model.lazycb != nothing
            ret = model.lazycb(grbcb)
            if ret == :Exit
                terminate(model.inner)
            end
        end
    elseif where == CB_MIP
        state = :MIPInfo
        grbcb = GurobiCallbackData(grbrawcb, state, where, [0.0])
        if model.infocb != nothing
            ret = model.infocb(grbcb)
            if ret == :Exit
                terminate(model.inner)
            end
        end
    end
    return convert(Cint,0)
end

# User callback function should be of the form:
# callback(cbdata::MathProgCallbackData)
# return :Exit to indicate an error

function setmathprogcallback!(model::XpressMathProgModel)

    @windows_only WORD_SIZE == 64 || error("Callbacks not currently supported on Win32. Use 64-bit Julia with 64-bit Gurobi.")
    grbcallback = cfunction(mastercallback, Cint, (Ptr{Void}, Ptr{Void}, Cint, Ptr{Void}))
    ret = @grb_ccall(setcallbackfunc, Cint, (Ptr{Void}, Ptr{Void}, Any), model.inner.ptr_model, grbcallback, model)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end
=#


text:
@testset "GKL - svdsolve full" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T, (n,n))
            alg = GKL(orth = orth, krylovdim = 2*n, maxiter = 1, tol = precision(T))
            S, lvecs, rvecs, info = @constinferred svdsolve(wrapop(A), wrapvec2(A[:,1]), n, :LR, alg)

            @test S ≈ svdvals(A)

            U = hcat(unwrapvec2.(lvecs)...)
            V = hcat(unwrapvec.(rvecs)...)
            @test U'*U ≈ I
            @test V'*V ≈ I
            @test A*V ≈ U*Diagonal(S)
        end
    end
end

@testset "GKL - svdsolve iteratively" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for orth in (cgs2, mgs2, cgsr, mgsr)
            A = rand(T, (2*N,N))
            v = rand(T, (2*N,))
            n₁ = div(n, 2)
            alg = GKL(orth = orth, krylovdim = n, maxiter = 10, tol = precision(T), eager = true)
            S, lvecs, rvecs, info = @constinferred svdsolve(wrapop(A), wrapvec2(v), n₁, :LR, alg)

            l = info.converged
            @test S[1:l] ≈ svdvals(A)[1:l]

            U = hcat(unwrapvec2.(lvecs)...)
            V = hcat(unwrapvec.(rvecs)...)
            @test U[:,1:l]'*U[:,1:l] ≈ I
            @test V[:,1:l]'*V[:,1:l] ≈ I

            R = hcat(unwrapvec2.(info.residual)...)
            @test A' * U ≈ V * Diagonal(S)
            @test A * V ≈ U * Diagonal(S) + R
        end
    end
end


text:
sess = Session()

@testset "fem_to_gauss_points" begin
m = 10
n = 10
h = 0.1
u = rand((m+1)*(n+1))
u0 = fem_to_gauss_points(u, m, n, h)
@test run(sess, fem_to_gauss_points(constant(u), m, n, h)) ≈ u0
end


text:
# We wish to extend operations to identically named methods dispatched
# on `Machine`s and `NodalMachine`s. For example, we have from the model API
#
# `predict(model::M, fitresult, X) where M<:Supervised`
#
# but want also want
#
# `predict(machine::Machine, X)` where `X` is data
#
# and "networks.jl" requires us to define
#
# `predict(machine::NodalMachine, X)` where `X` is data
#
# and we would like the syntactic sugar (for `X` a node):
#
# `predict(machine::NodalMachine, X::Node)=node(predict, machine, X)`
#
# (If an operation has zero arguments, we cannot achieve the last
# desire because of ambiguity with the preceding one.)

## TODO: need to add checks on the arguments of
## predict(::AbstractMachine, ) and transform(::AbstractMachine, )

for operation in (:predict, :predict_mean, :predict_mode, :predict_median,
                  :transform, :inverse_transform)
    ex = quote
        function $(operation)(machine::AbstractMachine{M}, args...) where M
            if isdefined(machine, :fitresult) || M <: Static
                return $(operation)(machine.model, machine.fitresult, args...)
            else
                throw(error("$machine has not been trained."))
            end
        end
        $(operation)(machine::Machine; rows=:) =
            $(operation)(machine, selectrows(machine.args[1], rows))
        $(operation)(machine::NodalMachine, args::AbstractNode...) =
            node($(operation), machine, args...)
    end
    eval(ex)
end

# the zero argument special cases:
"""
    fitted_params(mach)

Return the learned parameters for a machine `mach` that has been
`fit!`, for example the coefficients in a linear model.

This is a named tuple and human-readable if possible.

If `mach` is a machine for a composite model, then the returned value
has keys `machines` and `fitted_params_given_machine`, whose
corresponding values are a vector of (nodal) machines appearing in the
underlying learning network, and a dictionary of reports keyed on
those machines.

```julia
using MLJ
X, y = @load_crabs;
pipe = @pipeline MyPipe(
    std = Standardizer(),
    clf = @load LinearBinaryClassifier pkg=GLM
)
mach = machine(MyPipe(), X, y) |> fit!
fp = fitted_params(mach)
machs = fp.machines
2-element Array{Any,1}:
 NodalMachine{LinearBinaryClassifier{LogitLink}} @ 1…57
 NodalMachine{Standardizer} @ 7…33

fp.fitted_params_given_machine[machs[1]]
(coef = [121.05433477939319, 1.5863921128182814,
         61.0770377473622, -233.42699281787324, 72.74253591435117],
 intercept = 10.384459260848505,)
```

"""
function fitted_params(machine::AbstractMachine)
    if isdefined(machine, :fitresult)
        return fitted_params(machine.model, machine.fitresult)
    else
        throw(error("$machine has not been trained."))
    end
end


text:
module _Sine

using Brainstorm._Math._MPFR: add!, precision!

const guard_digits = 5
const bits_per_digit = log2(10)

pi(digits::T) where {T<:Integer} = begin
    setprecision(8) do
        target_prec = ceil(T, bits_per_digit * (digits + guard_digits))
        _pi(target_prec)
    end
end

# http://numbers.computation.free.fr/Constants/Pi/iterativePi.html
_pi(target_prec) = begin
    prec = 20
    setprecision(prec)
    α = big(355/113)

    while prec < target_prec
        prec = min(3*prec, target_prec)
        precision!(α, prec); setprecision(prec)

        add!(α, sin(α))
    end
    α
end

end


text:
"""
`Rot3DCoeffsEquiv: ` storing recursively precomputed coefficients for a
rotation-invariant basis.
"""

using LinearAlgebra: dot

struct Rot3DCoeffsEquiv{T,L}<: R3DC{T}
   vals::Vector{Dict}
   cg::ClebschGordan{T}
end

struct MRangeEq{N, T2}
   ll::SVector{N, Int}
   cartrg::T2
   L::Int
end

Base.length(mr::MRangeEq) = sum(_->1, _mrange(mr.ll, mr.L))

_mrange(ll, L) = MRangeEq(ll, Iterators.Stateful(
               filter((x) -> abs(sum(x))<= L, Tuple.(CartesianIndices(ntuple(i -> -ll[i]:ll[i], length(ll)))))
                     ),L)

function Base.iterate(mr::MRangeEq, args...)
	while true
		if isempty(mr.cartrg)
   		return nothing
		end
		mpre = popfirst!(mr.cartrg)
		return SVector(mpre), nothing
	end
end

function get0(L::Int,T=Float64)
	if L== 0
		return T(0)
	elseif L == 1
		return @SArray zeros(Complex{T},3,3)
		#SMatrix{3, 3, Complex{T}, 9}(0, 0, 0, 0, 0, 0, 0, 0, 0)
	else
		#For L== 2 @SArray zeros(Complex{T},3,3,3)
		ErrorException("Only types for L in {0,1} implemented")
	end
end

dicttype(A::Rot3DCoeffsEquiv, N::Integer) = dicttype(A::Rot3DCoeffsEquiv,Val(N))

dicttype(A::Rot3DCoeffsEquiv{T,L},::Val{N}) where {T,L,N} =
   Dict{Tuple{SVector{N,Int}, SVector{N,Int}, SVector{N,Int}}, typeof(get0(L,T))}

Rot3DCoeffsEquiv(L,T=Float64) = Rot3DCoeffsEquiv{T,L}(Dict[], ClebschGordan(T))

function get_vals(A::Rot3DCoeffsEquiv, valN::Val{N}) where {N}
	if length(A.vals) < N
		for n = length(A.vals)+1:N
			push!(A.vals, dicttype(A,n)())
		end
	end
   	return A.vals[N]::dicttype(A,valN)
end

function (A::Rot3DCoeffsEquiv{T,L})(ll::StaticVector{N},
							mm::StaticVector{N},
							kk::StaticVector{N}) where {T,L, N}
   if       abs(sum(mm)) > L ||
			abs(sum(kk)) > L ||
			!all(abs.(mm) .<= ll) ||
			!all(abs.(kk) .<= ll)
	  return get0(L,T)
   end
   vals = get_vals(A, Val(N))  # this should infer the type!
   key = _key(ll, mm, kk)
   if haskey(vals, key)
	  val  = vals[key]
   else
	  val = _compute_val(A, key...)
	  vals[key] = val
   end
   return val
end

function (A::Rot3DCoeffsEquiv{T,0})(ll::StaticVector{1},
							mm::StaticVector{1},
							kk::StaticVector{1}) where {T}
   if ll[1] == mm[1] == kk[1] == 0
	  return T(8 * pi^2)
   else
	  return get0(0,T)
   end
end

function (A::Rot3DCoeffsEquiv{T,1})(ll::StaticVector{1},
							mm::StaticVector{1},
							kk::StaticVector{1}) where {T}
   if ll[1] == 1 && abs(mm[1]) <= 1 && abs(kk[1]) <= 1
	  return  rmatrices[(mm[1],kk[1])]
   else
	  return get0(1,T)
   end
end

#function get0val(m,k,L)

#end
rmatrices = Dict(
  (-1,-1) => SMatrix{3, 3, ComplexF64, 9}(1/6, 1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
  (-1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
  (-1,1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, -1im/6, 0, -1im/6, 1/6, 0, 0, 0, 0),
  (0,-1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
  (0,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, 0, 0, 1/3),
  (0,1) => SMatrix{3, 3, ComplexF64, 9}(0, 0, -1/(3*sqrt(2)), 0, 0, -1im/(3*sqrt(2)), 0, 0, 0),
  (1,-1) => SMatrix{3, 3, ComplexF64, 9}(-1/6, 1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0),
  (1,0) => SMatrix{3, 3, ComplexF64, 9}(0, 0, 0, 0, 0, 0, -1/(3*sqrt(2)), 1im/(3*sqrt(2)), 0),
  (1,1) => SMatrix{3, 3, ComplexF64, 9}(1/6, -1im/6, 0, 1im/6, 1/6, 0, 0, 0, 0)
  )

function (A::Rot3DCoeffsEquiv{T,L})(ll::StaticVector{1},
							 mm::StaticVector{1},
							 kk::StaticVector{1}) where {T,L}
	ErrorException("Not implemented for L = " + string(L))
end

function _compute_val(A::Rot3DCoeffsEquiv{T,L}, ll::StaticVector{N},
                                          mm::StaticVector{N},
                                          kk::StaticVector{N}) where {T,L,N}

     val = get0(L,T)
     llp = ll[1:N-2]
     mmp = mm[1:N-2]
     kkp = kk[1:N-2]
     for j = abs(ll[N-1]-ll[N]):(ll[N-1]+ll[N])
        if abs(kk[N-1]+kk[N]) > j || abs(mm[N-1]+mm[N]) > j
           continue
        end
  		cgk = try
  			A.cg(ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
  		catch
  			@show (ll[N-1], kk[N-1], ll[N], kk[N], j, kk[N-1]+kk[N])
  			get0(L,T)
  		end
  		cgm = A.cg(ll[N-1], mm[N-1], ll[N], mm[N], j, mm[N-1]+mm[N])
  		if cgk * cgm  != 0
  			val += cgk * cgm * A( SVector(llp..., j),
  								       SVector(mmp..., mm[N-1]+mm[N]),
  								       SVector(kkp..., kk[N-1]+kk[N]) )
  		end
     end
     return val
 end

function vec3_symm_basis(A::Rot3DCoeffsEquiv{T,1},
						 nn::SVector{N, TN},
						 ll::SVector{N, Int}) where {T, N, TN}
	Ure, Mre = re_basis(A, ll)
	G = _gramian(nn, ll, Ure, Mre)
   	S = svd(G)
   	rk = rank(G; rtol =  1e-7)
	Urpe = S.U[:, 1:rk]'
	Utemp = Diagonal(sqrt.(S.S[1:rk])) * Urpe * Ure
	Utilde = zeros(SArray{Tuple{3},Complex{T},1,3}, rk, length(Mre))
	for α in 1:rk
		for (imu, mu) in enumerate(Mre)
			for k in 1:3
				for (i, mm) in enumerate(Mre)
					Utilde[α,imu] += Utemp[α,i] * A(ll,mu,mm)[:,k]
				end
			end
		end
	end

	return Utilde, Mre
end

function re_basis(A::Rot3DCoeffsEquiv{T}, ll::SVector; ordered=false) where {T}
	GGre, Mre = _gramianERot(A, ll)
	S = svd(GGre)
	rk = rank(Diagonal(S.S))
	return Diagonal(sqrt.(S.S[1:rk])) * S.U[:, 1:rk]', Mre
end

function _gramianERot(A::Rot3DCoeffsEquiv{T,1}, ll::SVector{N}) where{T,N}
	Mre = collect(_mrange(ll, 1))
  	len = length(Mre)*3
	GGre = zeros(Complex{T}, len, len)

	for (i1, (k1,mm1)) in enumerate(Iterators.product(1:3, Mre)), (i2, (k2,mm2)) in enumerate(Iterators.product(1:3, Mre))
		for mu1 in Mre, mu2 in Mre
			#show(typeof(dot(A(ll, mu2, mm2)[:,k2],A(ll, mu1, mm1)[:,k1])))
			GGre[i1, i2] += dot(A(ll, mu2, mm2)[:,k2],A(ll, mu1, mm1)[:,k1])
      	end
	end
   	return GGre, Mre
	#@assert all(abs.(imag(GG)) .<= imag_tol)
	#return real.(GG)
end


text:
# Bayesian logistic regression posterior abstract type
using StatsFuns: logistic

abstract SteinLogisticRegressionPosterior <: SteinPosterior

# each row of betas is assumed to be a beta sample
function loglikelihood(d::SteinLogisticRegressionPosterior,
                       beta::Array{Float64, 1};
                       idx=1:length(d.y))
    X = d.X[idx,:]
    y = d.y[idx]
    # condprobs is a N x b matrix
    condprobs = logistic(y .* (X * beta))
    sum(log(condprobs))
end

function gradloglikelihood(d::SteinLogisticRegressionPosterior,
                           beta::Array{Float64,1};
                           idx=1:length(d.y))
    X = d.X[idx,:]
    y = d.y[idx]
    # logitscores if passed to logit would be conditional probs
    logitscores = vec(y .* (X * beta))
    yX = y .* X

    gradients = logistic(-logitscores) .* yX
    vec(sum(gradients, 1))
end

function diaghessian(d::SteinLogisticRegressionPosterior,
                     beta::Array{Float64,1};
                     idx=1:length(d.y))
    X = d.X[idx,:]
    y = d.y[idx]
    # logitscores if passed to logit would be conditional probs
    logitscores = vec(y .* (X * beta))
    obsweights = logistic(-logitscores) .* logistic(logitscores)

    dhessian = sum(broadcast(*, X .^ 2, obsweights), 1)
    priordhessian = priordiaghessian(d, beta)

    vec(dhessian) + priordhessian
end


text:
using PolyChaos, Test

deg = 5
op = OrthoPoly("genlaguerre",deg,Dict(:shape=>1.23))
opq = OrthoPolyQ(op)

n = 10
X = [ exp(3 + rand()), exp.(3 .+ rand(n)) ]
@testset "evaluate univariate basis at point(s)" begin
    for x in X
        @test evaluate(collect(0:deg),x,op.α,op.β) == evaluate(collect(0:deg),x,op) == evaluate(collect(0:deg),x,opq)
        @test evaluate(x,op) == evaluate(x,opq)
        for d in 0:deg
            @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op) == evaluate(d,x,opq)
        end
    end
end

ops = [  OrthoPoly("genlaguerre",deg,Dict(:shape=>1.23)),
        OrthoPoly("genhermite",deg,Dict(:mu=>3.4)),
        OrthoPoly("legendre",deg),
        OrthoPoly("jacobi",deg,Dict(:shape_a=>4.3,:shape_b=>10.))
     ]
opqs = map(OrthoPolyQ, ops)
mop = MultiOrthoPoly(ops,deg)
mopq = MultiOrthoPoly(opqs,deg)

nop = length(ops)
X = [ exp.(3 .+ rand(nop)) , exp.(3 .+ rand(n,nop)) ]

@testset "evaluate multivariate basis at point(s)" begin
    for x in X, op in [mop, mopq]
        @test evaluate(op.ind,x,coeffs(op)...) == evaluate(op.ind,x,op) == evaluate(x,op)
        for d in eachrow(op.ind[1:3,:])
            @test evaluate(collect(d),x,coeffs(op)...) == evaluate(collect(d),x,op)
        end
    end
end