text:
# function linsolve_tridiagonal(left_diag::Vector, main_diag::Vector, right_diag::Vector, rhs::Vector)
function linsolve_tridiagonal(left_diag, main_diag, right_diag, rhs)
    n_row = length(main_diag)
    n = length(left_diag)
    if n+1 != n_row
        error("Invalid dimension for left diag")
    end
    if length(right_diag) != n
        error("Invalid dimension for right diag")
    end

    x = zeros(n_row)

    a = copy(main_diag)
    r = copy(rhs)

    # fw elimination
    @inbounds for ii = 2:n_row
        a[ii] -= left_diag[ii-1] * right_diag[ii-1] / a[ii-1]
        r[ii] -= left_diag[ii-1] * r[ii-1] / a[ii-1]
    end

    # bw elimination
    x[n_row] = r[n_row] / a[n_row]
    @inbounds for ii = (n_row-1):-1:1
        x[ii] = (r[ii] - x[ii+1] * right_diag[ii]) / a[ii]
    end

    return x
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""
`module  JAC.PhotoEmission`
    ... a submodel of JAC that contains all methods for computing Einstein coefficients, oscillator strength, etc. between
        some initial and final-state multiplets.    ***done***
"""
module PhotoEmission

    using Printf, ..AngularMomentum, ..Basics, ..Defaults, ..InteractionStrength, ..ManyElectron, ..Radial,
                  ..SpinAngular, ..TableStrings

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

        + multipoles              ::Array{EmMultipoles}     ... Specifies the (radiat. field) multipoles to be included.
        + gauges                  ::Array{UseGauge}         ... Gauges to be included into the computations.
        + calcAnisotropy          ::Bool                    ... True, if the anisotropy (structure) functions are to be
                                                                calculated and false otherwise
        + printBefore             ::Bool                    ... True, if all energies and lines are printed before comput.
        + lineSelection           ::LineSelection          ... Specifies the selected levels, if any.
        + photonEnergyShift       ::Float64                 ... An overall energy shift for all photon energies.
        + mimimumPhotonEnergy     ::Float64                 ... minimum transition energy for which (photon) transitions
                                                                are included into the computation.
        + maximumPhotonEnergy     ::Float64                 ... maximum transition energy for which (photon) transitions
                                                                are included.
    """
    struct Settings
        multipoles                ::Array{EmMultipole,1}
        gauges                    ::Array{UseGauge}
        calcAnisotropy            ::Bool
        printBefore               ::Bool
        lineSelection             ::LineSelection
        photonEnergyShift         ::Float64
        mimimumPhotonEnergy       ::Float64
        maximumPhotonEnergy       ::Float64
    end

    """
    `PhotoEmission.Settings()`  ... constructor for the default values of radiative line computations
    """
    function Settings()
        Settings(EmMultipole[E1], UseGauge[Basics.UseCoulomb], false, false, LineSelection(), 0., 0., 10000.)
    end

    """
    `PhotoEmission.Settings(set::PhotoEmission.Settings;`

            multipoles::=..,        gauges=..,          calcAnisotropy=..,          printBefore=..,
            line=..,         selectedLines=..,   photonEnergyShift=..,       mimimumPhotonEnergy=..,
            maximumPhotonEnergy=..)

        ... constructor for modifying the given PhotoEmission.Settings by 'overwriting' the previously selected parameters.
    """
    function Settings(set::PhotoEmission.Settings;
        multipoles::Union{Nothing,Array{EmMultipole,1}}=nothing,    gauges::Union{Nothing,Array{UseGauge}}=nothing,
        calcAnisotropy::Union{Nothing,Bool}=nothing,                printBefore::Union{Nothing,Bool}=nothing,
        lineSelection::Union{Nothing,LineSelection}=nothing,
        photonEnergyShift::Union{Nothing,Float64}=nothing,          mimimumPhotonEnergy::Union{Nothing,Float64}=nothing,
        maximumPhotonEnergy::Union{Nothing,Float64}=nothing)

        if  multipoles          == nothing   multipolesx          = set.multipoles              else  multipolesx          = multipoles            end
        if  gauges              == nothing   gaugesx              = set.gauges                  else  gaugesx              = gauges                end
        if  calcAnisotropy      == nothing   calcAnisotropyx      = set.calcAnisotropy          else  calcAnisotropyx      = calcAnisotropy        end
        if  printBefore         == nothing   printBeforex         = set.printBefore             else  printBeforex         = printBefore           end
        if  lineSelection       == nothing   lineSelectionx       = set.lineSelection           else  lineSelectionx       = lineSelection         end
        if  photonEnergyShift   == nothing   photonEnergyShiftx   = set.photonEnergyShift       else  photonEnergyShiftx   = photonEnergyShift     end
        if  mimimumPhotonEnergy == nothing   mimimumPhotonEnergyx = set.mimimumPhotonEnergy     else  mimimumPhotonEnergyx = mimimumPhotonEnergy   end
        if  maximumPhotonEnergy == nothing   maximumPhotonEnergyx = set.maximumPhotonEnergy     else  maximumPhotonEnergyx = maximumPhotonEnergy   end

        Settings( multipolesx, gaugesx, calcAnisotropyx, printBeforex, lineSelectionx,    ## selectLinesx, selectedLinesx,
                  photonEnergyShiftx, mimimumPhotonEnergyx, maximumPhotonEnergyx)
    end

    # `Base.show(io::IO, settings::PhotoEmission.Settings)`  ... prepares a proper printout of the variable settings::PhotoEmissionSettings.
    function Base.show(io::IO, settings::PhotoEmission.Settings)
        println(io, "multipoles:             $(settings.multipoles)  ")
        println(io, "gauges:                 $(settings.gauges)  ")
        println(io, "calcAnisotropy:         $(settings.calcAnisotropy)  ")
        println(io, "printBefore:            $(settings.printBefore)  ")
        println(io, "lineSelection:          $(settings.lineSelection)  ")
        println(io, "photonEnergyShift:      $(settings.photonEnergyShift)  ")
        println(io, "mimimumPhotonEnergy:    $(settings.mimimumPhotonEnergy)  ")
        println(io, "maximumPhotonEnergy:    $(settings.maximumPhotonEnergy)  ")
    end

    """
    `struct  PhotoEmission.Channel`
        ... defines a type for a single radiative emission/absorption channel that specifies the multipole, gauge and amplitude.

        + multipole         ::EmMultipole        ... Multipole of the photon emission/absorption.
        + gauge             ::EmGauge            ... Gauge for dealing with the (coupled) radiation field.
        + amplitude         ::Complex{Float64}   ... Amplitude of this multiple channel.
    """
    struct Channel
        multipole           ::EmMultipole
        gauge               ::EmGauge
        amplitude           ::Complex{Float64}
    end

    # `Base.show(io::IO, channel::PhotoEmission.Channel)`  ... prepares a proper printout of the variable channel::PhotoEmission.Channel.
    function Base.show(io::IO, channel::PhotoEmission.Channel)
        print(io, "PhotoEmission.Channel($(channel.multipole), $(channel.gauge), amp = $(channel.amplitude)) ")
    end

    """
    `struct  PhotoEmission.Line`
        ... defines a type for a radiative line that may include the definition of sublines and their corresponding amplitudes.

        + initialLevel   ::Level               ... initial-(state) level
        + finalLevel     ::Level               ... final-(state) level
        + omega          ::Float64             ... Transition frequency of this line; can be shifted w.r.t. the level energies.
        + photonRate     ::EmProperty          ... Total rate of this line.
        + angularBeta    ::EmProperty          ... Angular beta_2 coefficient.
        + hasSublines    ::Bool                ... Determines whether the sublines are defined in terms of their multipolarity, amplitude, or not.
        + channels       ::Array{PhotoEmission.Channel,1}  ... List of radiative (photon) channels
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        omega            ::Float64
        photonRate       ::EmProperty
        angularBeta      ::EmProperty
        hasSublines      ::Bool
        channels         ::Array{PhotoEmission.Channel,1}
    end

    """
    `PhotoEmission.Line(initialLevel::Level, finalLevel::Level, photonRate::Float64)`
        ... constructor an radiative line between a specified initial and final level.
    """
    function Line(initialLevel::Level, finalLevel::Level, omega::Float64, photonRate::EmProperty)
       Line(initialLevel, finalLevel, omega, photonRate, EmProperty(0., 0.), false, 0, PhotoEmission.Channel[])
    end

    # `Base.show(io::IO, line::PhotoEmissionLine)`  ... prepares a proper printout of the variable line::PhotoEmission.Line.
    function Base.show(io::IO, line::PhotoEmission.Line)
        println(io, "initialLevel:         $(line.initialLevel)  ")
        println(io, "finalLevel:           $(line.finalLevel)  ")
        println(io, "omega:                $(line.omega)  ")
        println(io, "photonRate:           $(line.photonRate)  ")
        println(io, "angularBeta:          $(line.angularBeta)  ")
        println(io, "hasSublines:          $(line.hasSublines)  ")
        println(io, "channels:             $(line.channels)  ")
    end

    """
    `PhotoEmission.amplitude(kind::String, Mp::EmMultipole, gauge::EmGauge, omega::Float64, finalLevel::Level, initialLevel::Level,
                             grid::Radial.Grid; display::Bool=false, printout::Bool=false)`
        ... to compute the kind = (absorption or emission) amplitude  <alpha_f J_f || O^(Mp, kind) || alpha_i J_i> for the
            interaction with  photon of multipolarity Mp and for the given transition energy and gauge. A value::ComplexF64 is
            returned. The amplitude value is printed to screen if display=true.
    """
    function amplitude(kind::String, Mp::EmMultipole, gauge::EmGauge, omega::Float64, finalLevel::Level, initialLevel::Level,
                       grid::Radial.Grid; display::Bool=false, printout::Bool=false)

        if      kind == "emission"
        #-------------------------
            if  initialLevel.basis.subshells == finalLevel.basis.subshells
                iLevel = initialLevel;   fLevel = finalLevel
            else
                subshells = Basics.merge(initialLevel.basis.subshells, finalLevel.basis.subshells)
                iLevel    = Level(initialLevel, subshells)
                fLevel    = Level(finalLevel, subshells)
            end

            nf = length(fLevel.basis.csfs);    ni = length(iLevel.basis.csfs)
            if  printout   printstyled("Compute radiative $(Mp) matrix of dimension $nf x $ni in the initial- and final-state bases " *
                                       "for the transition [$(iLevel.index)-$(fLevel.index)] ... ", color=:light_green)    end
            matrix = zeros(ComplexF64, nf, ni)
            #
            for  r = 1:nf
                if  fLevel.basis.csfs[r].J != fLevel.J      ||  fLevel.basis.csfs[r].parity  != fLevel.parity    continue    end
                for  s = 1:ni
                    if  iLevel.basis.csfs[s].J != iLevel.J  ||  iLevel.basis.csfs[s].parity  != iLevel.parity    continue    end
                    ##x wa = Basics.compute("angular coefficients: 1-p, Grasp92", 0, Mp.L, fLevel.basis.csfs[r], iLevel.basis.csfs[s])
                    ##x # Calculate the spin-angular coefficients
                    ##x if  Defaults.saRatip()
                    ##x     waR = Basics.compute("angular coefficients: 1-p, Grasp92", 0, Mp.L, finalLevel.basis.csfs[r], initialLevel.basis.csfs[s])
                    ##xx    wa  = waR
                    ##x end
                    ##x if  Defaults.saGG()
                    subshellList = fLevel.basis.subshells
                    opa = SpinAngular.OneParticleOperator(Mp.L, plus, true)
                    wa  = SpinAngular.computeCoefficients(opa, fLevel.basis.csfs[r], iLevel.basis.csfs[s], subshellList)
                    ##x     wa  = waG
                    ##x end
                    ##x if  Defaults.saRatip() && Defaults.saGG() && true
                    ##x     if  length(waR) != 0     println("\n>> Angular coeffients from GRASP/MCT   = $waR ")    end
                    ##x     if  length(waG) != 0     println(  ">> Angular coeffients from SpinAngular = $waG ")    end
                    ##x end
                    #
                    ##x @show finalLevel.basis  ## .csfs[r]
                    ##x @show initialLevel.basis  ## .csfs[s]
                    me = 0.
                    ##x println("Enter loop ... r=$r   s=$s  length(wa) = $(length(wa))")
                    for  coeff in wa
                        ##x ja = Basics.subshell_2j(finalLevel.basis.orbitals[coeff.a].subshell)
                        ##x jb = Basics.subshell_2j(initialLevel.basis.orbitals[coeff.b].subshell)
                        MbaCheng  = InteractionStrength.MbaEmissionCheng(Mp, gauge, omega, fLevel.basis.orbitals[coeff.a],
                                                                                iLevel.basis.orbitals[coeff.b], grid)
                        me = me + coeff.T * MbaCheng  ## sqrt( (ja + 1)/(jb + 1) ) * ##   0.707106781186548 *
                    end
                    ##x println("r = $r, s = $s, me = $me")
                    matrix[r,s] = me
                end
            end
            if  printout   printstyled("done. \n", color=:light_green)    end
            amplitude = transpose(fLevel.mc) * matrix * iLevel.mc
            #
            #
        elseif  kind == "absorption"
        #---------------------------
            iLevel = finalLevel;   fLevel = initialLevel
            amplitude = PhotoEmission.amplitude("emission", Mp, gauge, omega, fLevel, iLevel, grid, printout=printout)
            amplitude = conj(amplitude)
        else    error("stop a")
        end

        if  display
            println("    < level=$(finalLevel.index) [J=$(finalLevel.J)$(string(finalLevel.parity))] ||" *
                    " O^($Mp, $kind) ($omega a.u., $gauge) ||" *
                    " $(initialLevel.index) [$(initialLevel.J)$(string(initialLevel.parity))] >  = $amplitude  ")
        end

        return( amplitude )
    end

    """
    `PhotoEmission.computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid, settings::PhotoEmission.Settings;
                                output=true)`
        ... to compute the radiative transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{PhotoEmission.Lines} is returned.
    """
    function  computeLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid, settings::PhotoEmission.Settings; output=true)
        # Define a common subshell list for both multiplets
        subshellList = Basics.generate("subshells: ordered list for two bases", finalMultiplet.levels[1].basis, initialMultiplet.levels[1].basis)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=true)
        ##x Defaults.setDefaults("standard grid", grid)
        println("")
        printstyled("PhotoEmission.computeLines(): The computation of the transition amplitudes and properties starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        ##x println("basis = $(finalMultiplet.levels[1].basis) ")
        lines = PhotoEmission.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    PhotoEmission.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = PhotoEmission.Line[]
        for  line in lines
            newLine = PhotoEmission.computeAmplitudesProperties(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        PhotoEmission.displayRates(stdout, newLines)
        if  settings.calcAnisotropy    PhotoEmission.displayAnisotropies(stdout, newLines)    end
        PhotoEmission.displayLifetimes(stdout, newLines)
        #
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   PhotoEmission.displayRates(iostream, newLines)
                           PhotoEmission.displayAnisotropies(stdout, newLines)
                           PhotoEmission.displayLifetimes(iostream, newLines)
        end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `PhotoEmission.computeLinesCascade(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                       settings::PhotoEmission.Settings; output::Bool=true, printout::Bool=true)`
        ... to compute the radiative transition amplitudes and all properties as requested by the given settings. The computations
            and printout is adapted for larger cascade computations by including only lines with at least one channel and by sending
            all printout to a summary file only. A list of lines::Array{PhotoEmission.Lines} is returned.
    """
    function  computeLinesCascade(finalMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                                  settings::PhotoEmission.Settings; output=true, printout::Bool=true)
        # Define a common subshell list for both multiplets
        subshellList = Basics.generate("subshells: ordered list for two bases", finalMultiplet.levels[1].basis, initialMultiplet.levels[1].basis)
        Defaults.setDefaults("relativistic subshell list", subshellList; printout=false)
        lines = PhotoEmission.determineLines(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        # if  settings.printBefore    PhotoEmission.displayLines(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = PhotoEmission.Line[]
        for  (i,line)  in  enumerate(lines)
            if  rem(i,500) == 0    println("> Radiative line $i:")   end
            newLine = PhotoEmission.computeAmplitudesProperties(line, grid, settings, printout=printout)
            #
            # Don't add this line if it does not contribute to the decay
            wa = 0.
            for  ch in newLine.channels   wa = wa + abs(ch.amplitude)^2    end
            if   wa == 0.    continue    end
            push!( newLines, newLine)
        end
        # Print all results to a summary file, if requested
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   PhotoEmission.displayRates(iostream, newLines)    end
        #
        if    output    return( newLines )
        else            return( nothing )
        end
    end

    """
    `PhotoEmission.computeAmplitudesProperties(line::PhotoEmission.Line, grid::Radial.Grid, settings::Einstein.Settings; printout::Bool=true)`
        ... to compute all amplitudes and properties of the given line; a line::Einstein.Line is returned for which the amplitudes and
            properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::PhotoEmission.Line, grid::Radial.Grid, settings::PhotoEmission.Settings; printout::Bool=true)
        global JAC_counter
        newChannels = PhotoEmission.Channel[];    rateC = 0.;    rateB = 0.
        for channel in line.channels
            #
            amplitude = PhotoEmission.amplitude("emission", channel.multipole, channel.gauge, line.omega,
                                                line.finalLevel, line.initialLevel, grid, printout=printout)
            #
            push!( newChannels, PhotoEmission.Channel( channel.multipole, channel.gauge, amplitude) )
            if       channel.gauge == Basics.Coulomb     rateC = rateC + abs(amplitude)^2
            elseif   channel.gauge == Basics.Babushkin   rateB = rateB + abs(amplitude)^2
            elseif   channel.gauge == Basics.Magnetic    rateB = rateB + abs(amplitude)^2;   rateC = rateC + abs(amplitude)^2
            end
        end
        #
        # Calculate the photonrate and angular beta if requested
        wa = 8.0pi * Defaults.getDefaults("alpha") * line.omega / (Basics.twice(line.initialLevel.J) + 1) *
                                                      (Basics.twice(line.finalLevel.J) + 1)
        photonrate  = EmProperty(wa * rateC, wa * rateB)
        angularBeta = EmProperty(-9., -9.)
        line = PhotoEmission.Line( line.initialLevel, line.finalLevel, line.omega, photonrate, angularBeta, true, newChannels)
        return( line )
    end

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

    """
    `PhotoEmission.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PhotoEmission.Settings)`
        ... to determine a list of PhotoEmission Line's for transitions between the levels from the given initial- and
            final-state multiplets and by taking into account the particular selections and settings for this computation;
            an Array{PhotoEmission.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::PhotoEmission.Settings)
        lines = PhotoEmission.Line[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    omega = iLevel.energy - fLevel.energy   + settings.photonEnergyShift
                    if  omega <= settings.mimimumPhotonEnergy  ||  omega > settings.maximumPhotonEnergy    continue   end

                    channels = PhotoEmission.determineChannels(fLevel, iLevel, settings)
                    if   length(channels) == 0   continue   end
                    push!( lines, PhotoEmission.Line(iLevel, fLevel, omega, EmProperty(0., 0.), EmProperty(0., 0.), true, channels) )
                end
            end
        end
        return( lines )
    end

    """
    `PhotoEmission.displayAnisotropies(stream::IO, lines::Array{PhotoEmission.Line,1})`
        ... to list all energies and anisotropy parameters of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayAnisotropies(stream::IO, lines::Array{PhotoEmission.Line,1})
        nx = 153
        println(stream, " ")
        println(stream, "  Anisotropy (structure) functions:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        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(14, "Energy"   ; na=4);
        sb = sb * TableStrings.center(14,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(20, "Multipoles";   na=0);
        sa = sa * TableStrings.center(14, "f_2 (Coulomb)";   na=3);
        sa = sa * TableStrings.center(14, "f_2 (Babushkin)"; na=3);
        sa = sa * TableStrings.center(14, "f_4 (Coulomb)";   na=3);
        sa = sa * TableStrings.center(14, "f_4 (Babushkin)"; na=3);
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            f2Coulomb   = 0.0im;    f2Babushkin   = 0.0im;    f4Coulomb = 0.0im;    f4Babushkin = 0.0im;   mpList = EmMultipole[]
            normCoulomb = 0.;       normBabushkin = 0.
            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)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.omega)) * "    "
            # Now compute the anisotropy parameter from the amplitudes
            for  ch in line.channels
                if  !(ch.multipole  in  mpList)    push!(mpList, ch.multipole)    end
                for  chp  in line.channels
                    #
                    if  ch.gauge  in  [ EmGauge("Coulomb"), EmGauge("Magnetic")]
                        if  ch == chp    normCoulomb = normCoulomb + (abs(ch.amplitude)^2)    end
                        angL  = AngularJ64(ch.multipole.L);   p  = Basics.multipole_p(ch.multipole)
                        angLp = AngularJ64(ch.multipole.L);   pp = Basics.multipole_p(chp.multipole)
                        angJi = line.initialLevel.J;    angJf = line.finalLevel.J
                        #
                        f2Coulomb = f2Coulomb + 1.0im^(chp.multipole.L + pp - ch.multipole.L - p + 16)                         *
                                                AngularMomentum.phaseFactor([angJf, +1, angJi, +1, AngularJ64(3)])         *
                                                sqrt( (2ch.multipole.L+1) * (2chp.multipole.L+1) )                             *
                                    AngularMomentum.ClebschGordan(angL, AngularM64(1), angLp, AngularM64(-1), AngularJ64(2), AngularM64(0) ) *
                                                (1. + (-1.)^(ch.multipole.L + p + chp.multipole.L + pp - 2) )                  *
                                                AngularMomentum.Wigner_6j(angL, angLp, AngularJ64(2), angJi, angJi, angJf) *
                                                conj(ch.amplitude) * chp.amplitude
                        #
                        f4Coulomb = f4Coulomb + 1.0im^(chp.multipole.L + pp - ch.multipole.L - p + 16)                         *
                                                AngularMomentum.phaseFactor([angJf, +1, angJi, +1, AngularJ64(5)])         *
                                                sqrt( (2ch.multipole.L+1) * (2chp.multipole.L+1) )                             *
                                    AngularMomentum.ClebschGordan(angL, AngularM64(1), angLp, AngularM64(-1), AngularJ64(4), AngularM64(0) ) *
                                                (1. + (-1.)^(ch.multipole.L + p + chp.multipole.L + pp - 2) )                  *
                                                AngularMomentum.Wigner_6j(angL, angLp, AngularJ64(4), angJi, angJi, angJf) *
                                                conj(ch.amplitude) * chp.amplitude
                    end
                    #
                    if  ch.gauge  in  [ EmGauge("Babushkin"), EmGauge("Magnetic")]
                        if  ch == chp    normBabushkin = normBabushkin + (abs(ch.amplitude)^2)    end
                        angL  = AngularJ64(ch.multipole.L);   p  = Basics.multipole_p(ch.multipole)
                        angLp = AngularJ64(ch.multipole.L);   pp = Basics.multipole_p(chp.multipole)
                        angJi = line.initialLevel.J;    angJf = line.finalLevel.J
                        #
                        f2Babushkin = f2Babushkin + 1.0im^(chp.multipole.L + pp - ch.multipole.L - p + 16)                         *
                                                    AngularMomentum.phaseFactor([angJf, +1, angJi, +1, AngularJ64(3)])         *
                                                    sqrt( (2ch.multipole.L+1) * (2chp.multipole.L+1) )                             *
                                      AngularMomentum.ClebschGordan(angL, AngularM64(1), angLp, AngularM64(-1), AngularJ64(2), AngularM64(0) ) *
                                                    (1. + (-1.)^(ch.multipole.L + p + chp.multipole.L + pp - 2) )                  *
                                                    AngularMomentum.Wigner_6j(angL, angLp, AngularJ64(2), angJi, angJi, angJf) *
                                                    conj(ch.amplitude) * chp.amplitude
                        #
                        f4Babushkin = f4Babushkin + 1.0im^(chp.multipole.L + pp - ch.multipole.L - p + 16)                         *
                                                    AngularMomentum.phaseFactor([angJf, +1, angJi, +1, AngularJ64(5)])         *
                                                    sqrt( (2ch.multipole.L+1) * (2chp.multipole.L+1) )                             *
                                      AngularMomentum.ClebschGordan(angL, AngularM64(1), angLp, AngularM64(-1), AngularJ64(4), AngularM64(0) ) *
                                                    (1. + (-1.)^(ch.multipole.L + p + chp.multipole.L + pp - 2) )                  *
                                                    AngularMomentum.Wigner_6j(angL, angLp, AngularJ64(4), angJi, angJi, angJf) *
                                                    conj(ch.amplitude) * chp.amplitude
                    end
                end
            end
            f2Coulomb   = f2Coulomb   / normCoulomb   * sqrt(Basics.twice(line.initialLevel.J) + 1) / 2.
            f4Coulomb   = f4Coulomb   / normCoulomb   * sqrt(Basics.twice(line.initialLevel.J) + 1) / 2.
            f2Babushkin = f2Babushkin / normBabushkin * sqrt(Basics.twice(line.initialLevel.J) + 1) / 2.
            f4Babushkin = f4Babushkin / normBabushkin * sqrt(Basics.twice(line.initialLevel.J) + 1) / 2.
            #
            sa = sa * TableStrings.flushleft( 17, TableStrings.multipoleList(mpList);  na=3)
            sa = sa * TableStrings.flushright(15, @sprintf("%.8e", f2Coulomb.re);    na=3)
            sa = sa * TableStrings.flushright(15, @sprintf("%.8e", f2Babushkin.re);  na=3)
            sa = sa * TableStrings.flushright(15, @sprintf("%.8e", f4Coulomb.re);    na=3)
            sa = sa * TableStrings.flushright(15, @sprintf("%.8e", f4Babushkin.re);  na=3)
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `PhotoEmission.displayLifetimes(stream::IO, lines::Array{PhotoEmission.Line,1})`
        ... to list all lifetimes as derived from the selected lines. A neat table is printed but nothing is returned otherwise.
    """
    function  displayLifetimes(stream::IO, lines::Array{PhotoEmission.Line,1})
        # Determine all initial levels (and their level information) for printing the lifetimes
        ilevels = Int64[];   istr = String[]
        for  i = 1:length(lines)
            ii = lines[i].initialLevel.index;
            if  !(ii in ilevels)
               sa  = "  ";    sym = LevelSymmetry( lines[i].initialLevel.J, lines[i].initialLevel.parity )
               sa = sa * TableStrings.center(10, TableStrings.level(lines[i].initialLevel.index); na=2)
               sa = sa * TableStrings.center(10, string(sym); na=4)
               sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", lines[i].initialLevel.energy)) * "    "
               push!( ilevels, ii);    push!( istr, sa )
            end
        end
        # Determine the lifetime (in a.u.) of the selected initial levels
        irates = Basics.EmProperty[]
        for  ii in  ilevels
            waCoulomb = 0.;    waBabushkin = 0.
            for  i = 1:length(lines)
                if   lines[i].initialLevel.index == ii
                    waCoulomb   = waCoulomb   + lines[i].photonRate.Coulomb
                    waBabushkin = waBabushkin + lines[i].photonRate.Babushkin
                end
            end
            push!(irates, EmProperty( waCoulomb, waBabushkin) )
        end

        nx = 105
        println(stream, " ")
        println(stream, "  PhotoEmission lifetimes (as derived from these computations):")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                              sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                              sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(12, "Level energy"   ; na=3);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(12, "Used Gauge"    ; na=6);                     sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(26, "Lifetime"; na=3);
        sb = sb * TableStrings.center(26, "[a.u.]"*"          "*TableStrings.inUnits("time"); na=5)
        sa = sa * TableStrings.center(12, "Decay widths"; na=4);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  ii = 1:length(ilevels)
            sa = istr[ii]
            sa = sa * "Coulomb          " * @sprintf("%.6e",              1.0/irates[ii].Coulomb)     * "  "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("time: from atomic",   1.0/irates[ii].Coulomb) )   * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic",     irates[ii].Coulomb) )
            println(stream, sa)
            sa = repeat(" ", length(istr[ii]) )
            sa = sa * "Babushkin        " * @sprintf("%.6e",              1.0/irates[ii].Babushkin)   * "  "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("time: from atomic",   1.0/irates[ii].Babushkin) ) * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic",     irates[ii].Babushkin) )
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `PhotoEmission.displayLines(lines::Array{PhotoEmission.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{PhotoEmission.Line,1})
        nx = 95
        println(" ")
        println("  Selected radiative lines:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        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(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(30, "List of multipoles"; na=4);             sb = sb * TableStrings.hBlank(34)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "  ";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                           fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", line.omega)) * "    "
            mpGaugeList = Tuple{EmMultipole,EmGauge}[]
            for  i in 1:length(line.channels)
                push!( mpGaugeList, (line.channels[i].multipole, line.channels[i].gauge) )
            end
            sa = sa * TableStrings.multipoleGaugeTupels(50, mpGaugeList)
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        println(" ")
        #
        return( nothing )
    end

    """
    `PhotoEmission.displayRates(stream::IO, lines::Array{PhotoEmission.Line,1})`
        ... to list all results, energies, rates, etc. of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayRates(stream::IO, lines::Array{PhotoEmission.Line,1})
        nx = 145
        println(stream, " ")
        println(stream, "  Einstein coefficients, transition rates and oscillator strengths:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        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(12, "Energy"   ; na=4);
        sb = sb * TableStrings.center(12,TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center( 9, "Multipole"; na=0);                         sb = sb * TableStrings.hBlank(10)
        sa = sa * TableStrings.center(11, "Gauge"    ; na=4);                         sb = sb * TableStrings.hBlank(17)
        sa = sa * TableStrings.center(26, "A--Einstein--B"; na=3);
        sb = sb * TableStrings.center(26, TableStrings.inUnits("rate")*"          "*TableStrings.inUnits("rate"); na=2)
        sa = sa * TableStrings.center(11, "Osc. strength"    ; na=3);                 sb = sb * TableStrings.hBlank(17)
        sa = sa * TableStrings.center(12, "Decay widths"; na=4);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            for  ch in line.channels
                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)
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("energy: from atomic", line.omega)) * "    "
                sa = sa * TableStrings.center(9,  string(ch.multipole); na=4)
                sa = sa * TableStrings.flushleft(11, string(ch.gauge);  na=2)
                chRate =  8pi * Defaults.getDefaults("alpha") * line.omega / (Basics.twice(line.initialLevel.J) + 1) * (abs(ch.amplitude)^2) *
                                                                 (Basics.twice(line.finalLevel.J) + 1)
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to Einstein A",  line, chRate)) * "  "
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to Einstein B",  line, chRate)) * "    "
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to g_f",         line, chRate)) * "    "
                sa = sa * @sprintf("%.6e", Basics.recast("rate: radiative, to decay width", line, chRate)) * "      "
                ##x sa = sa * @sprintf("%.5e  %.5e", ch.amplitude.re, ch.amplitude.im) * "        "
                println(stream, sa)
            end
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module

ext: jl


text:
poisson(T, n) = sparse(Tridiagonal(fill(T(-1), n-1),
                        fill(T(2), n), fill(T(-1), n-1)))
poisson(n) = poisson(Float64, n)


text:
module aoc_14

export read_input

const mask_rx = r"mask\s*=\s*([10X]{36})"
const assign_rx = r"mem\[([0-9]*)\]\s=\s([0-9]*)"

abstract type Command end

struct Mask <: Command
    zero_mask::UInt64
    one_mask::UInt64
end

struct Assign <: Command
    address::UInt64
    value::UInt64
end

Mask(zero_mask::BitVector, one_mask::BitVector) = Mask(reverse(zero_mask).chunks[1], reverse(one_mask).chunks[1])
function Mask(m::AbstractString)
    Mask(convert(BitVector,[c .== '0' for c in m]), convert(BitVector,[c .== '1' for c in m]))
end
function Assign(address::AbstractString, value::AbstractString)
    Assign( Base.parse(Int64, address), Base.parse(Int64, value) )
end

function read_input(filename)
    lines = readlines(filename)
    commands = Vector{Command}();
    for line in lines
        if ( (m = match(assign_rx,line) ) != nothing)
            push!( commands, Assign(m.captures[1], m.captures[2]) )
        elseif ( ( m = match(mask_rx,line) ) != nothing)
            push!( commands, Mask(m.captures[1]) )
        else
            error("Cannot parse \"" * line * "\"")
        end
    end
    return commands
end

function parse(commands::Vector{Command})
    mask = Mask( falses( 36 ), falses( 36 ) )
    addresses = [assign.address for assign in commands if assign isa Assign]
    mem = zeros(UInt64, maximum(addresses) )
    for command in commands
        mask = evaluate!( command, mask, mem)
    end
    println(sum(mem))
    return mem
end

function evaluate!( command::Mask, mask::Mask, args...)
    return command
end

function evaluate!( assign::Assign, mask::Mask, mem::Vector{UInt64})
    mem[ assign.address ] = ~( ~( assign.value | mask.one_mask ) | mask.zero_mask )
    return mask
end

# Part 2

const thirtysixbitmask = UInt64(2^36-1)
const bitscan36 = 1 .<< (0:35)

getxs( command::Mask ) = ~command.zero_mask & ~command.one_mask & thirtysixbitmask
getxpos( mask::UInt64 ) = findall( mask .& bitscan36 .!= 0 )
getxpos( command::Mask ) = getxpos( getxs(command) )

function get_floating_memory( address::UInt64, mask::Mask)
    xs = getxs( mask )
    xpos = getxpos( xs )
    addresses = zeros(UInt64, 2^length(xpos))
    address = address | mask.one_mask
    address = address & ~xs
    for a = 1:length(addresses)
        bits = (UInt64(a) - 1) .& bitscan36[1:length(xpos)] .!= 0
        addresses[a] = address | sum(bits .<< (xpos .- 1))
    end
    return addresses
end

function parse2(commands::Vector{Command})
    mask = Mask( falses( 36 ), falses( 36 ) )
    # addresses = [assign.address for assign in commands if assign isa Assign]
    #mem = zeros(UInt64, 2^30 )
    mem = Dict{UInt64,UInt64}()
    for command in commands
        mask = evaluate2!( command, mask, mem)
    end
    #println(sum(mem))
    return mem
end

function evaluate2!( command::Mask, mask::Mask, args...)
    println( command )
    return command
end
function evaluate2!( command::Assign, mask::Mask, mem::Dict{UInt64,UInt64})
    println( command.address )
    addresses = get_floating_memory( command.address, mask)
    for address in addresses
        mem[ address ] = command.value
    end
    return mask
end

end


text:
module SnowballStemmer

using BinDeps
using Libdl

export stemmer_types, Stemmer, release
export stem

include("stemmer.jl")
end


text:
using Documenter, GridapMHD

makedocs(;
    modules=[GridapMHD],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jesusbonilla/GridapMHD.jl/blob/{commit}{path}#L{line}",
    sitename="GridapMHD.jl",
    authors="Jesús Bonilla, Large Scale Scientific Computing",
    assets=String[],
)

deploydocs(;
    repo="github.com/jesusbonilla/GridapMHD.jl",
)


text:
############################################################################################
#Define general functions that need to be dispatched for new particle kernels

"""
$(SIGNATURES)
Initiate particle with given initial distribution.

# Examples
```julia
```

"""
function initial(_rng::Random.AbstractRNG, kernel::ParticleKernel) end

"""
$(SIGNATURES)
Propagate particle forward given current particle.

# Examples
```julia
```

"""
function transition(
    _rng::Random.AbstractRNG, kernel::ParticleKernel, val::AbstractMatrix{P}, iter::Integer
) where {P} end

"""
$(SIGNATURES)
Calculate log transtion probability from particle given particle history.

# Examples
```julia
```

"""
function ℓtransition(
    _rng::Random.AbstractRNG,
    kernel::ParticleKernel,
    valₜ::Union{P,AbstractVector{P}},
    val::AbstractMatrix{P},
    iter::Integer,
) where {P} end

############################################################################################
#Include
include("markov.jl")
include("semimarkov.jl")

############################################################################################
"""
$(SIGNATURES)
Inplace initiate particle given kernel.

# Examples
```julia
```

"""
function initial!(
    _rng::Random.AbstractRNG, kernel::ParticleKernel, val::AbstractMatrix{P}, iter::Integer
) where {P}
    @inbounds for idx in Base.OneTo(size(val, 1))
        val[idx, iter] = initial(_rng, kernel)
    end
    return nothing
end

"""
$(SIGNATURES)
Inplace propagate particle forward given current particle and kernel.

# Examples
```julia
```

"""
function transition!(
    _rng::Random.AbstractRNG, kernel::ParticleKernel, val::AbstractMatrix{P}, iter::Integer
) where {P}
    #!NOTE: Transition to s_t ~ s_1:t-1, trajectories:t-1
    @inbounds for idx in Base.OneTo(size(val, 1))
        val[idx, iter] = transition(_rng, kernel, view(val, idx, :), iter)
    end
    return nothing
end

############################################################################################
"""
$(SIGNATURES)
Inplace calculate log transtion probability from current particle given trajectory.

# Examples
```julia
```

"""
function ℓtransition!(
    ℓcontainer::Vector{F},
    valₜ::Union{P,AbstractArray{P}},
    kernel::ParticleKernel,
    val::AbstractMatrix{P},
    iter::Integer,
) where {F<:AbstractFloat,P}
    @inbounds for idx in Base.OneTo(size(val, 1))
        ℓcontainer[idx] = ℓtransition(valₜ, kernel, view(val, idx, :), iter)
    end
    return ℓcontainer
end

############################################################################################

"""
$(SIGNATURES)
Evaluate data given particle trajectory.

# Examples
```julia
```

"""
function ℓevidence(
    dataₜ::D, kernel::ParticleKernel, val::AbstractVector{P}, iter::Integer
) where {D,P}
    return logpdf(kernel.evidence(val, iter), dataₜ)
end

"""
$(SIGNATURES)
Inplace logpdf evaluation of data given current particle trajectory.

# Examples
```julia
```

"""
function ℓevidence!(
    ℓcontainer::Vector{F},
    dataₜ::D,
    kernel::ParticleKernel,
    val::AbstractMatrix{P},
    iter::Integer,
) where {F<:AbstractFloat,D,P}
    @inbounds for idx in Base.OneTo(size(val, 1))
        #!NOTE: Need to evaluate e_t ~ s_1:t, e_1:t-1
        ℓcontainer[idx] = ℓevidence(dataₜ, kernel, view(val, idx, :), iter)
    end
    return ℓcontainer
end

############################################################################################
#export
export initial,
    transition, ℓtransition, initial!, transition!, ℓtransition!, ℓevidence, ℓevidence!


text:
type MIDIFile
    format::UInt16 # The format of the file. Can be 0, 1 or 2
    timedivision::Int16 # The time division of the track. Ticks per beat.
    tracks::Array{MIDITrack, 1} # An array of tracks

    MIDIFile() = new(0,96,MIDITrack[])
end

# Reads a file into a MIDIFile data type
function readMIDIfile(filename::AbstractString)
    f = open(filename)

    MIDIfile = MIDIFile()
    # Check that it's a valid MIDI file - first four bytes should spell MThd
    mthd = join(map(char, read(f, UInt8, 4)))
    if mthd != MTHD
        error("Not a valid MIDI file. Expected first 4 bytes to spell 'MThd', got $(mthd)")
    end

    # Skip the next four bytes - this is the header size, and it's always equal to 6.
    skip(f, 4)

    # Read the format code. 0 = single track, 1 = multiple tracks, 2 = multiple songs
    # Remember - MIDI files store data in big-endian format, which is why ntoh is used
    MIDIfile.format = ntoh(read(f, UInt16))

    # Get the number of tracks and time division
    numberoftracks = ntoh(read(f, UInt16))
    MIDIfile.timedivision = ntoh(read(f, Int16))
    MIDIfile.tracks = [readtrack(f) for x in 1:numberoftracks]
    close(f)

    MIDIfile
end

# Writes a MIDI file to the given filename
function writeMIDIfile(filename::AbstractString, data::MIDIFile)
    f = open(filename, "w")

    write(f, convert(Array{UInt8, 1}, MTHD)) # File identifier
    write(f, hton(convert(Uint32, 6))) # Header length
    write(f, hton(data.format))
    write(f, hton(convert(UInt16, length(data.tracks))))
    write(f, hton(data.timedivision))

    map(track->writetrack(f, track), data.tracks)

    close(f)
end

export MIDIFile, readMIDIfile, writeMIDIfile


text:
__precompile__()
module NPZ

# NPZ file format is described in
# https://github.com/numpy/numpy/blob/v1.7.0/numpy/lib/format.py

using ZipFile, Compat

@static if VERSION >=  v"0.7.0-DEV.2575"
    import Base.CodeUnits
else
    # CodeUnits not yet supported by Compat but not needed in julia 0.6...
    # codeunits function in Compat returns uintX instead of codeunits
    # therefore this 'stump' type should work
    abstract type CodeUnits{U, S} end
end

export npzread, npzwrite

const NPYMagic = UInt8[0x93, 'N', 'U', 'M', 'P', 'Y']
const ZIPMagic = UInt8['P', 'K', 3, 4]
const Version = UInt8[1, 0]

const MaxMagicLen = maximum(length.([NPYMagic, ZIPMagic]))

const TypeMaps = [
    ("b1", Bool),
    ("i1", Int8),
    ("i2", Int16),
    ("i4", Int32),
    ("i8", Int64),
    ("u1", UInt8),
    ("u2", UInt16),
    ("u4", UInt32),
    ("u8", UInt64),
    ("f2", Float16),
    ("f4", Float32),
    ("f8", Float64),
    ("c8", Complex{Float32}),
    ("c16", Complex{Float64}),
]
const Numpy2Julia = Dict{String, DataType}()
for (s,t) in TypeMaps
    Numpy2Julia[s] = t
end

const Julia2Numpy = Dict{DataType, String}()

@static if VERSION >= v"0.4.0"
    function __init__()
        for (s,t) in TypeMaps
            Julia2Numpy[t] = s
        end
    end
else
    for (s,t) in TypeMaps
        Julia2Numpy[t] = s
    end
end

# Julia2Numpy is a dictionary that uses Types as keys.
# This is problematic for precompilation because the
# hash of a Type changes everytime Julia is run.
# The hash of the keys when NPZ is precompiled will
# not be the same as when it is later run. This can
# be fixed by rehashing the Dict when the module is
# loaded.

readle(ios::IO, ::Type{T}) where T = ltoh(read(ios, T)) # ltoh is inverse of htol

function writecheck(io::IO, x::Any)
    n = write(io, x) # returns size in bytes
    n == sizeof(x) || error("short write") # sizeof is size in bytes
end

# Endianness only pertains to multi-byte things
writele(ios::IO, x::AbstractVector{UInt8}) = writecheck(ios, x)
writele(ios::IO, x::AbstractVector{CodeUnits{UInt8, <:Any}}) = writecheck(ios, x)
# codeunits returns vector of CodeUnits in 7+, uint in 6
writele(ios::IO, x::AbstractString) = writele(ios, codeunits(x))

writele(ios::IO, x::UInt16) = writecheck(ios, htol(x))

function parsechar(s::AbstractString, c::Char)
    firstchar = s[firstindex(s)]
    if  firstchar != c
        error("parsing header failed: expected character '$c', found '$firstchar'")
    end
    SubString(s, nextind(s, 1))
end

function parsestring(s::AbstractString)
    s = parsechar(s, '\'')
    parts = split(s, '\'', limit = 2)
    length(parts) != 2 && error("parsing header failed: malformed string")
    parts[1], parts[2]
end

function parsebool(s::AbstractString)
    if SubString(s, firstindex(s), thisind(s, 4)) == "True"
        return true, SubString(s, nextind(s, 4))
    elseif SubString(s, firstindex(s), thisind(s, 5)) == "False"
        return false, SubString(s, nextind(s, 5))
    end
    error("parsing header failed: excepted True or False")
end

function parseinteger(s::AbstractString)
    isdigit(s[firstindex(s)]) || error("parsing header failed: no digits")
    tail_idx = findfirst(c -> !isdigit(c), s)
    if tail_idx == nothing
        intstr = SubString(s, firstindex(s))
    else
        intstr = SubString(s, firstindex(s), prevind(s, tail_idx))
        if s[tail_idx] == 'L' # output of firstindex should be a valid code point
            tail_idx = nextind(s, tail_idx)
        end
    end
    n = parse(Int, intstr)
    return n, SubString(s, tail_idx)
end

function parsetuple(s::AbstractString)
    s = parsechar(s, '(')
    tup = Int[]
    while true
        s = strip(s)
        if s[firstindex(s)] == ')'
            break
        end
        n, s = parseinteger(s)
        push!(tup, n)
        s = strip(s)
        if s[firstindex(s)] == ')'
            break
        end
        s = parsechar(s, ',')
    end
    s = parsechar(s, ')')
    Tuple(tup), s
end

function parsedtype(s::AbstractString)
    dtype, s = parsestring(s)
    c = dtype[firstindex(s)]
    t = SubString(dtype, nextind(s, 1))
    if c == '<'
        toh = ltoh
    elseif c == '>'
        toh = ntoh
    elseif c == '|'
        toh = identity
    else
        error("parsing header failed: unsupported endian character $c")
    end
    if !haskey(Numpy2Julia, t)
        error("parsing header failed: unsupported type $t")
    end
    (toh, Numpy2Julia[t]), s
end

struct Header{T,N,F<:Function}
    descr::F
    fortran_order::Bool
    shape::NTuple{N,Int}
end

Header{T}(descr::F, fortran_order, shape::NTuple{N,Int}) where {T,N,F} = Header{T,N,F}(descr, fortran_order, shape)
Base.size(hdr::Header) = hdr.shape
Base.eltype(hdr::Header{T}) where T = T
Base.ndims(hdr::Header{T,N}) where {T,N} = N

function parseheader(s::AbstractString)
    s = parsechar(s, '{')

    dict = Dict{String,Any}()
    T = Any
    for _ in 1:3
        s = strip(s)
        key, s = parsestring(s)
        s = strip(s)
        s = parsechar(s, ':')
        s = strip(s)
        if key == "descr"
            (descr, T), s = parsedtype(s)
            dict[key] = descr
        elseif key == "fortran_order"
            dict[key], s = parsebool(s)
        elseif key == "shape"
            dict[key], s = parsetuple(s)
        else
            error("parsing header failed: bad dictionary key")
        end
        s = strip(s)
        if s[firstindex(s)] == '}'
            break
        end
        s = parsechar(s, ',')
    end
    s = strip(s)
    s = parsechar(s, '}')
    s = strip(s)
    if s != ""
        error("malformed header")
    end
    Header{T}(dict["descr"], dict["fortran_order"], dict["shape"])
end

function readheader(f::IO)
    @compat b = read!(f, Vector{UInt8}(undef, length(NPYMagic)))
    if b != NPYMagic
        error("not a numpy array file")
    end
    @compat b = read!(f, Vector{UInt8}(undef, length(Version)))

    # support for version 2 files
    if b[1] == 1
        hdrlen = UInt32(readle(f, UInt16))
    elseif b[1] == 2
        hdrlen = UInt32(readle(f, UInt32))
    else
        error("unsupported NPZ version")
    end

    @compat hdr = ascii(String(read!(f, Vector{UInt8}(undef, hdrlen))))
    parseheader(strip(hdr))
end

function _npzreadarray(f, hdr::Header{T}) where {T}
    toh = hdr.descr
    if hdr.fortran_order
        @compat x = map(toh, read!(f, Array{T}(undef, hdr.shape)))
    else
        @compat x = map(toh, read!(f, Array{T}(undef, reverse(hdr.shape))))
        if ndims(x) > 1
            x = permutedims(x, collect(ndims(x):-1:1))
        end
    end
    ndims(x) == 0 ? x[1] : x
end

function npzreadarray(f::IO)
    hdr = readheader(f)
    _npzreadarray(f, hdr)
end

function samestart(a::AbstractVector, b::AbstractVector)
    nb = length(b)
    length(a) >= nb && view(a, 1:nb) == b
end

function _maybetrimext(name::AbstractString)
    fname, ext = splitext(name)
    if ext == ".npy"
        name = fname
    end
    name
end

"""
    npzread(filename::AbstractString, [vars])

Read a variable or a collection of variables from `filename`.
The input needs to be either an `npy` or an `npz` file.
The optional argument `vars` is used only for `npz` files.
If it is specified, only the matching variables are read in from the file.

!!! note "Zero-dimensional arrays"
    Zero-dimensional arrays are stripped while being read in, and the values that they
    contain are returned. This is a notable difference from numpy, where
    numerical values are written out and read back in as zero-dimensional arrays.

# Examples

```julia
julia> npzwrite("temp.npz", x = ones(3), y = 3)

julia> npzread("temp.npz") # Reads all variables
Dict{String,Any} with 2 entries:
  "x" => [1.0, 1.0, 1.0]
  "y" => 3

julia> npzread("temp.npz", ["x"]) # Reads only "x"
Dict{String,Array{Float64,1}} with 1 entry:
  "x" => [1.0, 1.0, 1.0]
```
"""
function npzread(filename::AbstractString, vars...)
    # Detect if the file is a numpy npy array file or a npz/zip file.
    f = open(filename)
    @compat b = read!(f, Vector{UInt8}(undef, MaxMagicLen))

    if samestart(b, ZIPMagic)
        fz = ZipFile.Reader(filename)
        data = npzread(fz, vars...)
        close(fz)
    elseif samestart(b, NPYMagic)
        seekstart(f)
        data = npzreadarray(f)
    else
        close(f)
        error("not a NPY or NPZ/Zip file: $filename")
    end
    close(f)
    return data
end

function npzread(dir::ZipFile.Reader,
    vars = map(f -> _maybetrimext(f.name), dir.files))

    Dict(_maybetrimext(f.name) => npzreadarray(f)
        for f in dir.files
            if f.name in vars || _maybetrimext(f.name) in vars)
end

"""
    readheader(filename, [vars...])

Return a header or a collection of headers corresponding to each variable contained in `filename`.
The header contains information about the `eltype` and `size` of the array that may be extracted using
the corresponding accessor functions.
"""
function readheader(filename::AbstractString, vars...)
    # Detect if the file is a numpy npy array file or a npz/zip file.
    f = open(filename)
    @compat b = read!(f, Vector{UInt8}(undef, MaxMagicLen))

    if samestart(b, ZIPMagic)
        fz = ZipFile.Reader(filename)
        data = readheader(fz, vars...)
        close(fz)
    elseif samestart(b, NPYMagic)
        seekstart(f)
        data = readheader(f)
    else
        close(f)
        error("not a NPY or NPZ/Zip file: $filename")
    end

    close(f)
    return data
end
function readheader(dir::ZipFile.Reader,
    vars = map(f -> _maybetrimext(f.name), dir.files))

    Dict(_maybetrimext(f.name) => readheader(f)
        for f in dir.files
            if f.name in vars || _maybetrimext(f.name) in vars)
end

function npzwritearray(
    f::IO, x::AbstractArray{UInt8}, T::DataType, shape)

    if !haskey(Julia2Numpy, T)
        error("unsupported type $T")
    end
    writele(f, NPYMagic)
    writele(f, Version)

    descr =  (ENDIAN_BOM == 0x01020304 ? ">" : "<") * Julia2Numpy[T]
    dict = "{'descr': '$descr', 'fortran_order': True, 'shape': $(Tuple(shape)), }"

    # The dictionary is padded with enough whitespace so that
    # the array data is 16-byte aligned
    n = length(NPYMagic)+length(Version)+2+length(dict)
    pad = (div(n+16-1, 16)*16) - n
    if pad > 0
        dict *= " "^(pad-1) * "\n"
    end

    writele(f, UInt16(length(dict)))
    writele(f, dict)
    if write(f, x) != length(x)
        error("short write")
    end
end

function npzwritearray(f::IO, x::AbstractArray)
    npzwritearray(f, reinterpret(UInt8, vec(x)), eltype(x), size(x))
end

function npzwritearray(f::IO, x::Number)
    npzwritearray(f, reinterpret(UInt8, [x]), typeof(x), ())
end

"""
    npzwrite(filename::AbstractString, x)

Write the variable `x` to the `npy` file `filename`.
Unlike `numpy`, the extension `.npy` is not appened to `filename`.

!!! warn "Warning"
    Any existing file with the same name will be overwritten.

# Examples

```julia
julia> npzwrite("abc.npy", zeros(3))

julia> npzread("abc.npy")
3-element Array{Float64,1}:
 0.0
 0.0
 0.0
```
"""
function npzwrite(filename::AbstractString, x)
    open(filename, "w") do f
        npzwritearray(f, x)
    end
end

"""
    npzwrite(filename::AbstractString, vars::Dict{<:AbstractString})
    npzwrite(filename::AbstractString, args...; kwargs...)

In the first form, write the variables in `vars` to an `npz` file named `filename`.

In the second form, collect the variables in `args` and `kwargs` and write them all
to `filename`. The variables in `args` are saved with names `arr_0`, `arr_1`
and so on, whereas the ones in `kwargs` are saved with the specified names.

Unlike `numpy`, the extension `.npz` is not appened to `filename`.

!!! warn "Warning"
    Any existing file with the same name will be overwritten.

# Examples

```julia
julia> npzwrite("temp.npz", Dict("x" => ones(3), "y" => 3))

julia> npzread("temp.npz")
Dict{String,Any} with 2 entries:
  "x" => [1.0, 1.0, 1.0]
  "y" => 3

julia> npzwrite("temp.npz", ones(2,2), x = ones(3), y = 3)

julia> npzread("temp.npz")
Dict{String,Any} with 3 entries:
  "arr_0" => [1.0 1.0; 1.0 1.0]
  "x"     => [1.0, 1.0, 1.0]
  "y"     => 3
```
"""
function npzwrite(filename::AbstractString, vars::Dict{<:AbstractString})
    dir = ZipFile.Writer(filename)

    if length(vars) == 0
        @warn "no data to be written to $filename. It might not be possible to read the file correctly."
    end

    for (k, v) in vars
        f = ZipFile.addfile(dir, k * ".npy")
        npzwritearray(f, v)
        close(f)
    end

    close(dir)
end

function npzwrite(filename::AbstractString, args...; kwargs...)
    dkwargs = Dict(string(k) => v for (k,v) in kwargs)
    dargs = Dict("arr_"*string(i-1) => v for (i,v) in enumerate(args))

    d = convert(Dict{String,Any}, merge(dargs, dkwargs))

    npzwrite(filename, d)
end

end # module


text:
export initialize

function initialize(model;groups=[],expr="",parallel=false)
  var=match(r"^.*(?==)",expr).match
  value=match(r"(?<==).*",expr).match

  for name in [groups;]
    p=model.groups[name][:processor]
    model.initialization*="@spawnat $(p+1) @eval $var=[Float64($value) for i=1:N]\n"
  end
end


text:
function on_date(obj::UnnormalizedSphericalHarmonicsProvider, arg0::AbsoluteDate)
    return jcall(obj, "onDate", UnnormalizedSphericalHarmonicsProvider_UnnormalizedSphericalHarmonics, (AbsoluteDate,), arg0)
end


text:
##ORTHOGONALITY HELPER
function consolidate_srcs(con_idxs::Dict{Integer,Vector{Integer}}, m::ICA_PWM_Model, obs_array::AbstractMatrix{<:Integer}, obs_lengths::AbstractVector{<:Integer}, bg_scores::AbstractMatrix{<:AbstractFloat}, contour::AbstractFloat,  models::AbstractVector{<:Model_Record}; iterates::Integer=length(m.sources)*2, remote=false)
    new_log_Li=-Inf;  iterate = 1
    T,O = size(obs_array); T=T-1; S = length(m.sources)
    new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
    a, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true)

    while new_log_Li <= contour && iterate <= iterates #until we produce a model more likely than the lh contour or exceed iterates
        new_sources=deepcopy(m.sources); new_mix=deepcopy(m.mix_matrix)
        clean=Vector{Bool}(trues(O))

        for host_src in filter(!in(vcat(values(con_idxs)...)), keys(con_idxs)) #copy mix information to the source to be consolidated on as host
            for cons_src in con_idxs[host_src]
                new_mix[:,host_src].=[new_mix[o,host_src] || new_mix[o,cons_src] for o in 1:size(new_mix,1)]
            end
        end

        remote ? (merger_m = deserialize(rand(models).path)) : (merger_m = remotecall_fetch(deserialize, 1, rand(models).path)) #randomly select a model to merge
        used_m_srcs=Vector{Int64}()

        for src in unique(vcat(values(con_idxs)...)) #replace all non-host sources with sources from a merger model unlike the one being removed
            distvec=[pwm_distance(m.sources[src][1],m_src[1]) for m_src in merger_m.sources]
            m_src=findmax(distvec)[2]
            while m_src in used_m_srcs
                distvec[m_src]=0.; m_src=findmax(distvec)[2]
                length(used_m_srcs)==length(merger_m.sources) && break;break
            end

            clean[new_mix[:,src]].=false #mark dirty any obs that start with the source
            new_sources[src]=merger_m.sources[m_src]
            new_mix[:,src]=merger_m.mix_matrix[:,m_src]
            clean[new_mix[:,src]].=false #mark dirty any obs that end with the source
            push!(used_m_srcs, m_src)
        end

        if consolidate_check(new_sources)[1] #if the new sources pass the consolidate check
            new_log_Li, cache = IPM_likelihood(new_sources, obs_array, obs_lengths, bg_scores, new_mix, true, true, cache, clean) #assess likelihood
        end

        iterate += 1
    end

    return ICA_PWM_Model("candidate","consolidated $(m.origin)", new_sources, m.source_length_limits, new_mix, new_log_Li)
end

function consolidate_check(sources::AbstractVector{<:Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer}}; thresh=CONSOLIDATE_THRESH, revcomp=REVCOMP)
    pass=true
    lengthδmat=[size(src1[1],1) - size(src2[1],1) for src1 in sources, src2 in sources]
    cons_idxs=Dict{Integer,Vector{Integer}}()
    for src1 in 1:length(sources), src2 in src1+1:length(sources)
        if lengthδmat[src1,src2]==0
            revcomp ? (info_condition = (pwm_distance(sources[src1][1],sources[src2][1]) < thresh || pwm_distance(sources[src1][1],revcomp_pwm(sources[src2][1])) < thresh)) : (info_condition = (pwm_distance(sources[src1][1],sources[src2][1]) < thresh))
            if info_condition
                if !in(src1,keys(cons_idxs))
                    cons_idxs[src1]=[src2]; pass=false
                else
                    push!(cons_idxs[src1], src2)
                end
            end
        end
    end

    return pass, cons_idxs
end

                function pwm_distance(pwm1,pwm2)
                    minwml=min(size(pwm1,1),size(pwm2,1))
                    return sum([euclidean(exp.(pwm1[pos,:]), exp.(pwm2[pos,:])) for pos in 1:minwml])/minwml
                end


text:
# Automatically-generated headers for cusolver
#
# DO NOT EDIT THIS FILE DIRECTLY.

# Automatically generated using Clang.jl

const CUSOLVER_VER_MAJOR = 10
const CUSOLVER_VER_MINOR = 5
const CUSOLVER_VER_PATCH = 0
const CUSOLVER_VER_BUILD = 218
const CUSOLVER_VERSION = CUSOLVER_VER_MAJOR * 1000 + CUSOLVER_VER_MINOR * 100 + CUSOLVER_VER_PATCH
const cusolverDnContext = Cvoid
const cusolverDnHandle_t = Ptr{cusolverDnContext}
const syevjInfo = Cvoid
const syevjInfo_t = Ptr{syevjInfo}
const gesvdjInfo = Cvoid
const gesvdjInfo_t = Ptr{gesvdjInfo}
const cusolverDnIRSParams = Cvoid
const cusolverDnIRSParams_t = Ptr{cusolverDnIRSParams}
const cusolverDnIRSInfos = Cvoid
const cusolverDnIRSInfos_t = Ptr{cusolverDnIRSInfos}
const cusolverDnParams = Cvoid
const cusolverDnParams_t = Ptr{cusolverDnParams}

@cenum cusolverDnFunction_t::UInt32 begin
    CUSOLVERDN_GETRF = 0
end

const cusolver_int_t = Cint

@cenum cusolverStatus_t::UInt32 begin
    CUSOLVER_STATUS_SUCCESS = 0
    CUSOLVER_STATUS_NOT_INITIALIZED = 1
    CUSOLVER_STATUS_ALLOC_FAILED = 2
    CUSOLVER_STATUS_INVALID_VALUE = 3
    CUSOLVER_STATUS_ARCH_MISMATCH = 4
    CUSOLVER_STATUS_MAPPING_ERROR = 5
    CUSOLVER_STATUS_EXECUTION_FAILED = 6
    CUSOLVER_STATUS_INTERNAL_ERROR = 7
    CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED = 8
    CUSOLVER_STATUS_NOT_SUPPORTED = 9
    CUSOLVER_STATUS_ZERO_PIVOT = 10
    CUSOLVER_STATUS_INVALID_LICENSE = 11
    CUSOLVER_STATUS_IRS_PARAMS_NOT_INITIALIZED = 12
    CUSOLVER_STATUS_IRS_PARAMS_INVALID = 13
    CUSOLVER_STATUS_IRS_PARAMS_INVALID_PREC = 14
    CUSOLVER_STATUS_IRS_PARAMS_INVALID_REFINE = 15
    CUSOLVER_STATUS_IRS_PARAMS_INVALID_MAXITER = 16
    CUSOLVER_STATUS_IRS_INTERNAL_ERROR = 20
    CUSOLVER_STATUS_IRS_NOT_SUPPORTED = 21
    CUSOLVER_STATUS_IRS_OUT_OF_RANGE = 22
    CUSOLVER_STATUS_IRS_NRHS_NOT_SUPPORTED_FOR_REFINE_GMRES = 23
    CUSOLVER_STATUS_IRS_INFOS_NOT_INITIALIZED = 25
    CUSOLVER_STATUS_IRS_INFOS_NOT_DESTROYED = 26
    CUSOLVER_STATUS_IRS_MATRIX_SINGULAR = 30
    CUSOLVER_STATUS_INVALID_WORKSPACE = 31
end

@cenum cusolverEigType_t::UInt32 begin
    CUSOLVER_EIG_TYPE_1 = 1
    CUSOLVER_EIG_TYPE_2 = 2
    CUSOLVER_EIG_TYPE_3 = 3
end

@cenum cusolverEigMode_t::UInt32 begin
    CUSOLVER_EIG_MODE_NOVECTOR = 0
    CUSOLVER_EIG_MODE_VECTOR = 1
end

@cenum cusolverEigRange_t::UInt32 begin
    CUSOLVER_EIG_RANGE_ALL = 1001
    CUSOLVER_EIG_RANGE_I = 1002
    CUSOLVER_EIG_RANGE_V = 1003
end

@cenum cusolverNorm_t::UInt32 begin
    CUSOLVER_INF_NORM = 104
    CUSOLVER_MAX_NORM = 105
    CUSOLVER_ONE_NORM = 106
    CUSOLVER_FRO_NORM = 107
end

@cenum cusolverIRSRefinement_t::UInt32 begin
    CUSOLVER_IRS_REFINE_NOT_SET = 1100
    CUSOLVER_IRS_REFINE_NONE = 1101
    CUSOLVER_IRS_REFINE_CLASSICAL = 1102
    CUSOLVER_IRS_REFINE_CLASSICAL_GMRES = 1103
    CUSOLVER_IRS_REFINE_GMRES = 1104
    CUSOLVER_IRS_REFINE_GMRES_GMRES = 1105
    CUSOLVER_IRS_REFINE_GMRES_NOPCOND = 1106
    CUSOLVER_PREC_DD = 1150
    CUSOLVER_PREC_SS = 1151
    CUSOLVER_PREC_SHT = 1152
end

@cenum cusolverPrecType_t::UInt32 begin
    CUSOLVER_R_8I = 1201
    CUSOLVER_R_8U = 1202
    CUSOLVER_R_64F = 1203
    CUSOLVER_R_32F = 1204
    CUSOLVER_R_16F = 1205
    CUSOLVER_R_16BF = 1206
    CUSOLVER_R_TF32 = 1207
    CUSOLVER_R_AP = 1208
    CUSOLVER_C_8I = 1211
    CUSOLVER_C_8U = 1212
    CUSOLVER_C_64F = 1213
    CUSOLVER_C_32F = 1214
    CUSOLVER_C_16F = 1215
    CUSOLVER_C_16BF = 1216
    CUSOLVER_C_TF32 = 1217
    CUSOLVER_C_AP = 1218
end

@cenum cusolverAlgMode_t::UInt32 begin
    CUSOLVER_ALG_0 = 0
    CUSOLVER_ALG_1 = 1
end

@cenum cusolverStorevMode_t::UInt32 begin
    CUBLAS_STOREV_COLUMNWISE = 0
    CUBLAS_STOREV_ROWWISE = 1
end

@cenum cusolverDirectMode_t::UInt32 begin
    CUBLAS_DIRECT_FORWARD = 0
    CUBLAS_DIRECT_BACKWARD = 1
end

const cusolverSpContext = Cvoid
const cusolverSpHandle_t = Ptr{cusolverSpContext}
const csrqrInfo = Cvoid
const csrqrInfo_t = Ptr{csrqrInfo}


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of helpers to print tables.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "@pt" begin
    # To get the output of the macro @pt, we must redirect the stdout.
    old_stdout = stdout
    in, out    = redirect_stdout()

    # The configurations must not interfere with the printings of any type
    # except for the `Dict`.
    @ptconf sortkeys = true

    # Test 1
    # --------------------------------------------------------------------------

    expected = """
.......................
: 1 :     2 :   3 : 4 :
:...:.......:.....:...:
: 1 : false : 1.0 : 1 :
: 2 :  true : 2.0 : 2 :
: 3 : false : 3.0 : 3 :
: 4 :  true : 4.0 : 4 :
: 5 : false : 5.0 : 5 :
: 6 :  true : 6.0 : 6 :
:...:.......:.....:...:
"""

    @ptconf tf = tf_ascii_dots
    @pt :header = ["1","2","3","4"] data

    result = String(readavailable(in))
    @test result == expected

    # Test 2
    # --------------------------------------------------------------------------

    expected = """
╭────────┬────────┬────────┬────────╮
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│   1    │ false  │  1.0   │   1    │
│   2    │  true  │  2.0   │   2    │
│   3    │ false  │  3.0   │   3    │
│   4    │  true  │  4.0   │   4    │
│   5    │ false  │  5.0   │   5    │
│   6    │  true  │  6.0   │   6    │
╰────────┴────────┴────────┴────────╯
"""

    @ptconf tf = tf_unicode_rounded
    @ptconf alignment = :c
    @pt data

    result = String(readavailable(in))
    @test result == expected

    # Test 3
    # --------------------------------------------------------------------------

    expected = """
.----------.----------.
| Column 1 | Column 2 |
|   Sub. 1 |   Sub. 2 |
:----------+----------:
|        1 |        2 |
|        3 |        4 |
|        5 |        6 |
'----------'----------'
"""

    header = ["Column 1" "Column 2"; "Sub. 1" "Sub. 2"]

    @ptconfclean
    @ptconf tf = tf_ascii_rounded
    @pt :header = header [1 2; 3 4; 5 6]

    result = String(readavailable(in))
    @test result == expected

    # Test 4
    # --------------------------------------------------------------------------

    expected = """
.-------.---------------------.
| Keys  | Values              |
| Int64 | String              |
:-------+---------------------:
| 1     | São José dos Campos |
| 2     | SP                  |
| 3     | Brasil              |
'-------'---------------------'
"""

    @ptconf sortkeys = true
    @ptconf alignment = :l
    @pt d = Dict(Int64(1) => "São José dos Campos",
                 Int64(2) => "SP",
                 Int64(3) => "Brasil")

    result = String(readavailable(in))
    @test result == expected

    # Restore the original stdout.
    close(in)
    close(out)
    redirect_stdout(old_stdout)

    # The expression after `@pt` must be evaluated. Hence, `d` must be a dict
    # with the defined elements.
    @test d isa Dict
    @test d[1] == "São José dos Campos"
    @test d[2] == "SP"
    @test d[3] == "Brasil"

    # Clean configurations to avoid test failures if this is run in the same
    # section.
    @ptconfclean
end


text:
using SimilaritySearch

function main(n, m, dim)
    db = MatrixDatabase(rand(Float32, dim, n))
    queries = MatrixDatabase(rand(Float32, dim, m))
    seq = ExhaustiveSearch(; db)
    k = 10
    @time searchbatch(seq, queries, k; parallel=Threads.nthreads() > 1)
end

@info "warming"
main(100, 10, 3)
@info "large benchmark"
I, D = main(100_000, 1000, 3)
size(I), size(D)

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
using QuadPencils
using Test

@testset "QuadPencils.jl" begin
    # Write your tests here.
end


text:
THIS_DIR = dirname(@__FILE__)
CSV_PATH = joinpath(THIS_DIR, "../../data/hr/hr.csv")

TREASURER_DEPT = "TREASURER"
TREASURER_HEAD = "CITY TREASURER"
ACC_HEAD = "DIR OF ACCOUNTING"
ACCS = ["AUDITOR IV", "ACCOUNTANT IV", "ACCOUNTANT III", "STAFF ASST", "ACCOUNTANT I"]

csv = readcsv(CSV_PATH, header=true)[1]
N = size(csv, 1)

dept_name_set = Set{String}()

tr_head = nothing
trs = []
acc_head = nothing
accs = []

for eid = 1:N
    department_name = csv[eid, 1]
    push!(dept_name_set, department_name)
    if department_name == TREASURER_DEPT
        position = csv[eid, 4]
        if position == TREASURER_HEAD
            tr_head = eid
        elseif position == ACC_HEAD
            acc_head = eid
            push!(trs, eid)
        elseif position in ACCS
            push!(accs, eid)
        else
            push!(trs, eid)
        end
    end
end

M = length(dept_name_set)

dept_data = 1:M
emp_data = 1:N

dept_name_data = String[]
dept_employee_data = Vector{Int}[]

emp_name_data = String[]
emp_position_data = String[]
emp_salary_data = Int[]
emp_department_data = Int[]
emp_manager_data = Nullable{Int}[]
emp_subordinate_data = Vector{Int}[]

did_by_name = Dict{String, Int}()

for eid = 1:N
    department_name, name, surname, position, salary = csv[eid, :]
    if !(department_name in keys(did_by_name))
        did = length(did_by_name)+1
        push!(dept_name_data, department_name)
        push!(dept_employee_data, Int[])
        did_by_name[department_name] = did
    end
    did = did_by_name[department_name]
    push!(dept_employee_data[did], eid)
    push!(emp_name_data, name * " " * surname[1:1])
    push!(emp_position_data, position)
    push!(emp_salary_data, salary)
    push!(emp_department_data, did)
    if department_name == TREASURER_DEPT
        if position == TREASURER_HEAD
            push!(emp_manager_data, nothing)
            push!(emp_subordinate_data, trs)
        elseif position == ACC_HEAD
            push!(emp_manager_data, tr_head)
            push!(emp_subordinate_data, accs)
        elseif position in ACCS
            push!(emp_manager_data, acc_head)
            push!(emp_subordinate_data, Int[])
        else
            push!(emp_manager_data, tr_head)
            push!(emp_subordinate_data, Int[])
        end
    else
        push!(emp_manager_data, nothing)
        push!(emp_subordinate_data, Int[])
    end
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
const GLSLScalarTypes = Union{Float32, Int32, UInt32}

"""
Statically sized uniform buffer.
Supports push!, but with fixed memory, so it will error after reaching
it's preallocated length.
"""
struct UniformBuffer{T, N}
    buffer::Buffer{T}
    offsets::NTuple{N, Int}
    elementsize::Int
    length::Int
end

"""
    Pre allocates an empty buffer with `max_batch_size` size
    which can be used to store multiple uniform blocks of type T
"""
function UniformBuffer(::Type{T}, max_batch_size = 1024, mode = GL_STATIC_DRAW) where T
    offsets, elementsize = std140_offsets(T)
    buffer = Buffer{T}(
        max_batch_size,
        elementsize * max_batch_size,
        GL_UNIFORM_BUFFER, mode
    )
    UniformBuffer(buffer, offsets, elementsize, 0)
end
"""
    Creates an Uniform buffer with the contents of `data`
"""
function UniformBuffer{T}(data::T, mode = GL_STATIC_DRAW) where T
    buffer = UniformBuffer(T, 1, mode)
    push!(buffer, data)
    buffer
end

"""
Returns the alignment of the `Type` of T as assumed in https://khronos.org/registry/OpenGL/specs/gl/glspec45.core.pdf#page=159,
returning a tuple with the first element being the 'base' alignment, and the second the total size inside memory.
"""
function glsl_alignment_size(T)
    function ceil4(i)
        while i%4 != 0
            i += 1
        end
        return i
    end
    T <: Bool && return sizeof(Int32), sizeof(Int32)
    N = sizeof(T)
    T <: GLSLScalarTypes && return N, N
    T <: Function && return sizeof(Vec4f0), sizeof(Vec4f0) # sizeof(EmptyStruct) padded to Vec4f0
    ET = eltype(T)
    N  = sizeof(ET)
    if T <: Matrix
        nrows, ncols = size(T)
        ncols        = ceil4(ncols)
        return div(ncols, 4) * N, ncols * nrows * N
    end
    if T <: Vector
        return ceil4(length(T)) * N, length(T) * N
    end
    @error "Struct $T not supported yet. Please help by implementing all rules from https://khronos.org/registry/OpenGL/specs/gl/glspec45.core.pdf#page=159"
end

function std140_offsets(::Type{T}) where T
    elementsize = 0
    offsets = if T <: GLSLScalarTypes
        elementsize = sizeof(T)
        (0,)
    else
        offset = 0
        offsets = ntuple(nfields(T)) do i
            ft = fieldtype(T, i)
            alignement, sz = glsl_alignment_size(ft)
            if offset % alignement != 0
                offset = (div(offset, alignement) + 1) * alignement
            end
            of = offset
            offset += sz
            of
        end
        elementsize = offset
        offsets
    end
    offsets, elementsize
end

Base.convert(::Type{UniformBuffer}, x) = UniformBuffer(x)
Base.convert(::Type{UniformBuffer}, x::UniformBuffer) = x
Base.eltype(::UniformBuffer{T}) where T = T

function Base.setindex!(buffer::UniformBuffer{T}, element::T, idx::Integer) where T
    if idx > length(buffer.buffer)
        throw(BoundsError(buffer, idx))
    end
    buff = buffer.buffer
    glBindBuffer(buff.buffertype, buff.id)
    dptr = Ptr{UInt8}(glMapBuffer(buff.buffertype, GL_WRITE_ONLY))
    for (offset, ptr, size) in iterate_fields(buffer, element, idx)
        unsafe_copy!(dptr + offset, ptr, size)
    end
    glUnmapBuffer(buff.buffertype)
    bind(buff, 0)
    element
end

function Base.push!(buffer::UniformBuffer{T}, element::T) where T
    buffer.length += 1
    buffer[buffer.length] = element
    buffer
end

function assert_blocksize(buffer::UniformBuffer, program, blockname::String)
    block_index = glGetUniformBlockIndex(program, blockname)
    blocksize_ref = Ref{GLint}(0)
    glGetActiveUniformBlockiv(
        program, block_index,
        GL_UNIFORM_BLOCK_DATA_SIZE, blocksize_ref
    )
    blocksize = blocksize_ref[]
    @assert buffer.elementsize * length(buffer.buffer) == blocksize
end

_getfield(x::GLSLScalarTypes, i) = x
_getfield(x, i) = getfield(x, i)

function iterate_fields(buffer::UniformBuffer{T, N}, x, index) where {T, N}
    offset = buffer.elementsize * (index - 1)
    x_ref = isimmutable(x) ? Ref(x) : x
    base_ptr = Ptr{UInt8}(pointer_from_objref(x_ref))
    ntuple(Val{N}) do i
        offset + buffer.offsets[i], base_ptr + fieldoffset(T, i), sizeof(fieldtype(T, i))
    end
end

extract_val(::Val{X}) where X = X

# function Base.setindex!{T <: Composable, N, TF}(x::UniformBuffer{T, N}, val::TF, field::Type{<: Field})
#     index = extract_val(FieldTraits.fieldindex(T, field))
#     if index === 0
#         throw(BoundsError(x, field))
#     end
#     val_conv = convert(fieldtype(T, index), val)
#     val_ref = if isbits(val)
#         Base.Ref(val)
#     elseif isimmutable(val)
#         error("Struct $TF contains pointers and can't be transferred to GPU")
#     else
#         pointer_from_objref(val)
#     end
#     buff = x.buffer
#     bind(buff) do
#         BufferSubData(buff.buffertype, x.offsets[index], sizeof(val_conv), val_ref)
#     end
#     x
# end

# function Base.getindex{T <: Composable, N}(x::UniformBuffer{T, N}, field::Type{<: Field})
#     index = extract_val(FieldTraits.fieldindex(T, field)[1])
#     if index == 0
#         throw(BoundsError(x, field))
#     end
#     ET = fieldtype(T, index)
#     val_ref = Ref{ET}()
#     bind(x.buffer) do
#         glGetBufferSubData(x.buffer.buffertype, x.offsets[index], sizeof(ET), val_ref)
#     end
#     val_ref[]
# end


text:
using LinearAlgebra
using IterativeSolvers
# using CUDA

# include("./utils/types.jl")
# using .Utils

dim = 819
light = 64

# CUDA.@profile begin
A = rand(dim, dim)
B = rand(light, light)
# end

function power_iteration(A, num_simulations::Int64)
    lambda, v = powm(A, tol=1e-3, maxiter=num_simulations)
    return v
end

# warm up
@time power_iteration(B, 1)

@time power_iteration(A, 16)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

#=
@testset "2d nonlinear elasticity: test nonhomogeneous boundary conditions and stress calculation" begin

    # field problem
    block = Problem(Elasticity, "BLOCK", 2)
    block.properties.formulation = :plane_stress
    block.properties.finite_strain = true
    block.properties.geometric_stiffness = true

    nodes = Dict{Int, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])

    element = Element(Quad4, [1, 2, 3, 4])
    update!(element, "geometry", nodes)
    update!(element, "youngs modulus", 288.0)
    update!(element, "poissons ratio", 1/3)

    plastic_parameters = Dict{Any, Any}("type" => JuliaFEM.ideal_plasticity!,
                                        "yield_surface" => Val{:von_mises},
                                        "params" => Dict("yield_stress" => 175.0))
    to_integ_points = Dict()
    map(x-> to_integ_points[x] = plastic_parameters, get_connectivity(element))
    update!(element, "plasticity", to_integ_points)
    push!(block, element)

    # boundary conditions
    bc = Problem(Dirichlet, "bc", 2, "displacement")
    bel1 = Element(Seg2, [1, 2])
    bel2 = Element(Seg2, [3, 4])
    bel3 = Element(Seg2, [4, 1])
    update!([bel1, bel2, bel3], "geometry", nodes)
    update!(bel1, "displacement 2", 0.0)
    update!(bel2, "displacement 2", 0.5)
    update!(bel3, "displacement 1", 0.0)
    push!(bc, bel1, bel2, bel3)

    solver = NonlinearSolver("solve block problem")
    solver.time = 1.0
    push!(solver, block, bc)
    solver()

    # from code aster
    eps_expected = [-2.08333312468287E-01, 6.25000000000000E-01, 0.0]
    sig_expected = [ 4.50685020821470E-06, 4.62857140373777E+02, 0.0]
    u3_expected =  [-2.36237356855269E-01, 5.00000000000000E-01]

    u3 = reshape(block.assembly.u, 2, 4)[:, 3]
    @info("u3 = $u3")
    @test isapprox(u3, u3_expected, atol=1.0e-5)
end

=#


text:
const REPO_URL = "https://github.com/TuringLang/TuringTutorials"
const CLONED_DIR = joinpath(REPO_DIR, "ClonedTuringTutorials")

"""
    clean_weave_cache()

On the one hand, we need `cache = :all` to have quick builds.
On the other hand, we don't need cache files committed to the repo which break the build.
Therefore, this method manually cleans the cache just to be sure.
"""
function clean_weave_cache()
    for (root, dirs, files) in walkdir(pkgdir(TuringTutorials); onerror=x->())
        if "cache" in dirs
            cache_dir = joinpath(root, "cache")
            rm(cache_dir; force=true, recursive=true)
        end
    end
end

"""
    clone_tutorials_output()

Ensure that `$CLONED_DIR` exists and contains the latest commit from the output branch for `$REPO_URL`.
"""
function clone_tutorials_output()
    branch = "artifacts"
    args = [
        "clone",
        "--depth=1",
        "--branch=$branch"
    ]
    if isdir(CLONED_DIR)
        try
            cd(CLONED_DIR) do
                run(`git checkout $branch`)
                run(`git pull --ff --allow-unrelated-histories`)
            end
        catch
            rm(CLONED_DIR; recursive=true, force=true)
            run(`git $args $REPO_URL $CLONED_DIR`)
            cd(CLONED_DIR) do
                run(`git checkout $branch`)
            end
        end
    else
        run(`git $args $REPO_URL $CLONED_DIR`)
    end
end

"""
    download_artifacts()

Explicitly copy all the updated tutorials from the artifacts branch.
This allows orphaning the artifacts branch on each deploy to ease debugging, cleanup old and
unused files, and have a smaller branch.
"""
function download_artifacts()
    if !isdir(CLONED_DIR)
        clone_tutorials_output()
    end
    T = tutorials()
    for tutorial in T
        for dir in ["html", "markdown", "notebook", "script"]
            from_dir = joinpath(CLONED_DIR, dir, tutorial)
            to_dir = joinpath(REPO_DIR, dir, tutorial)
            mkpath(to_dir)
            # from_dir is missing for new/renamed tutorials.
            if isdir(from_dir)
                cp(from_dir, to_dir; force=true)
            end
        end
    end
end

function file_changed(old_dir, new_dir, file)
    old_path = joinpath(old_dir, file)
    new_path = joinpath(new_dir, file)
    old = read(old_path, String)
    new = isfile(new_path) ? read(new_path, String) : ""
    return old != new
end

"""
    any_changes(tutorial::String)

Return whether there are any changes for the local source files, such as `.jmd` and `Manifest.toml`,
compared to the files in `$CLONED_DIR`.
"""
function any_changes(tutorial::String)
    old_dir = joinpath(CLONED_DIR, "tutorials", tutorial)
    new_dir = joinpath(REPO_DIR, "tutorials", tutorial)
    if isdir(old_dir)
        files = readdir(old_dir)
        files = filter(!=(WEAVE_LOG_FILE), files)
        return any(file_changed.(old_dir, new_dir, files))
    else
        # A newly added tutorial.
        return true
    end
end

"""
    changed_tutorials()

Return the tutorials which have changed compared to the output branch at $REPO_URL.
"""
function changed_tutorials()
    clone_tutorials_output()
    T = tutorials()
    changed = filter(any_changes, T)
    n = length(changed)
    println("Found changes for the tutorials $changed ($(n)/$(length(T)))")
    if length(T) == 0
        changed = first(T)
        println("Running the first tutorial to be able to verify that the CI job works.")
    end
    changed
end

max_stars_repo_name: staticfloat/TuringTutorials

max_issues_repo_name: staticfloat/TuringTutorials

max_forks_repo_name: staticfloat/TuringTutorials


text:
Base.@kwdef mutable struct Agent
    handle::Int
    earliest_departure::Int
    latest_arrival::Int
    initial_position::Tuple{Int,Int}
    initial_direction::Int
    target_position::Tuple{Int,Int}
end

function Agent(pyagent::Py)
    handle = pyconvert(Int, pyagent.handle) + 1
    earliest_departure = pyconvert(Int, pyagent.earliest_departure) + 1
    latest_arrival = pyconvert(Int, pyagent.latest_arrival) + 1
    initial_position = pyconvert(Tuple{Int,Int}, pyagent.initial_position) .+ 1
    initial_direction = pyconvert(Int, pyagent.initial_direction) + 1
    target_position = pyconvert(Tuple{Int,Int}, pyagent.target) .+ 1

    return Agent(;
        handle=handle,
        earliest_departure=earliest_departure,
        latest_arrival=latest_arrival,
        initial_position=initial_position,
        initial_direction=initial_direction,
        target_position=target_position,
    )
end

function initial_label(agent::Agent)
    return (agent.initial_position..., agent.initial_direction, DEPARTURE)
end

function target_label(agent::Agent)
    return (agent.target_position..., NO_DIRECTION, ARRIVAL)
end

function station_positions(agents::AbstractVector{Agent})
    initial_positions = unique(agent.initial_position for agent in agents)
    target_positions = unique(agent.target_position for agent in agents)
    return union(initial_positions, target_positions)
end