text:
"""
`module  JAC.RayleighCompton`
    ... a submodel of JAC that contains all methods for computing elastic Rayleigh and inelastic Compton photon scattering cross sections.
"""
module RayleighCompton

    using Printf, ..AngularMomentum, ..AtomicState, ..Basics, ..Defaults, ..ManyElectron, ..Radial, ..TableStrings

    """
    `struct  RayleighCompton.Settings`  ... defines a type for the settings in estimating lastic Rayleigh and inelastic Compton
                                            photon scattering cross sections

        + multipoles              ::Array{EmMultipole}      ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge}         ... Specifies the gauges to be included into the computations.
        + photonEnergies          ::Array{Float64,1}        ... List of photon energies in units defined by the user.
        + green                   ::Array{GreenChannel,1}   ... Precalculated and user-specified Green function of the ion.
        + calcAngular             ::Bool                    ... Calculate the angular distribution of the fluorescence photon.
        + calcStokes              ::Bool                    ... Calculate the Stokes parameters of the Rayleigh-scattered photon.
        + printBefore             ::Bool                    ... True, if all energies and lines are printed before their evaluation.
        + incidentStokes          ::ExpStokes               ... Stokes parameters of the incident radiation.
        + solidAngles             ::Array{SolidAngle,1}     ... List of solid angles [(theta_1, pho_1), ...].
        + lineSelection           ::LineSelection           ... Specifies the selected levels, if any.

    """
    struct Settings
        multipoles                ::Array{EmMultipole}
        gauges                    ::Array{UseGauge}
        photonEnergies            ::Array{Float64,1}
        green                     ::Array{AtomicState.GreenChannel,1}
        calcAngular               ::Bool
        calcStokes                ::Bool
        printBefore               ::Bool
        incidentStokes            ::ExpStokes
        solidAngles               ::Array{SolidAngle,1}
        lineSelection             ::LineSelection
    end

    """
    `RayleighCompton.Settings()`  ... constructor for the default values of Rayleigh-Compton photon-scattering estimates.
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], GreenChannel[], 0., false, false, false, ExpStokes(), SolidAngle[], LineSelection() )
    end

    # `Base.show(io::IO, settings::RayleighCompton.Settings)`  ... prepares a proper printout of the variable settings::RayleighCompton.Settings.
    function Base.show(io::IO, settings::RayleighCompton.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "gauges:                   $(settings.gauges)  ")
        println(io, "photonEnergies:           $(settings.photonEnergies)  ")
        println(io, "green:                     (settings.green)  ")
        println(io, "calcAngular:              $(settings.calcAngular)  ")
        println(io, "calcStokes:               $(settings.calcStokes)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "incidentStokes:           $(settings.incidentStokes)  ")
        println(io, "solidAngles:              $(settings.solidAngles)  ")
        println(io, "lineSelection:            $(settings.lineSelection)  ")
    end

    """
    `struct  RayleighCompton.ReducedChannel`
        ... defines a type for a Rayleigh-Compton reduced channel/ U^(K, Compton) amplitude to help characterize the scattering of
            a single photon with well-defined multipolarities; in JAC, a Rayleigh-Compton channel is characterized by a reduced amplitude
            U^(K, Compton) (level_f, multipole_2, J_nu, omega, multipole_1, level_i).

        + K              ::AngularJ64           ... Rank K of the channel.
        + multipole1     ::EmMultipole          ... Multipole of the incoming photon.
        + multipole2     ::EmMultipole          ... Multipole of the outgoing, scattered photon.
        + gauge          ::EmGauge              ... Gauge for dealing with the (coupled) radiation field.
        + omega1         ::Float64              ... Photon frequency of multipole1.
        + omega2         ::Float64              ... Photon frequency of multipole2.
        + omega          ::Float64              ... Photon frequency in denominator.
        + Jsym           ::LevelSymmetry        ... Symmetry J_nu of the intermediate levels.
        + amplitude      ::Complex{Float64}     ... Two-photon reduced amplitude associated with the given channel.
   """
    struct  ReducedChannel
        K                ::AngularJ64
        multipole1       ::EmMultipole
        multipole2       ::EmMultipole
        gauge            ::EmGauge
        omega1           ::Float64
        omega2           ::Float64
        omega            ::Float64
        Jsym             ::LevelSymmetry
        amplitude        ::Complex{Float64}
    end

    """
    `struct  RayleighCompton.Line`  ... defines a type for a Rayleigh-Compton scattering line that may include the definition of channels.

        + initialLevel     ::Level              ... initial-(state) level
        + finalLevel       ::Level              ... final-(state) level
        + inOmega          ::Float64            ... Energy of the incoming photon.
        + outOmega         ::Float64            ... Energy of the outgoing, scattered photon.
        + crossSection     ::EmProperty         ... Cross section for this photoionization.
        + channels         ::Array{RayleighCompton.ReducedChannel,1}    ... List of (reduced) RayleighCompton.Channels of this line.
    """
    struct  Line
        initialLevel       ::Level
        finalLevel         ::Level
        inOmega            ::Float64
        outOmega           ::Float64
        crossSection       ::EmProperty
        channels           ::Array{RayleighCompton.ReducedChannel,1}
    end

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

    """
    `RayleighCompton.computeAmplitudesProperties(line::RayleighCompton.Line, grid::Radial.Grid, settings::RayleighCompton.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::RayleighCompton.Line is returned for which
            the amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesProperties(line::RayleighCompton.Line, grid::Radial.Grid, settings::RayleighCompton.Settings)
        newChannels = RayleighCompton.ReducedChannel[]
        for ch in line.channels
            amplitude = RayleighCompton.computeReducedAmplitude(ch.K, line.finalLevel, ch.multipole2, ch.omega2, ch.Jsym, ch.omega,
                                                                ch.multipole1, ch.omega1, line.initialLevel, ch.gauge, grid, settings.green)
            push!( newChannels, ReducedChannel(ch.K, ch.multipole1, ch.multipole2, ch.gauge, ch.omega1, ch.omega2, ch.omega, ch.Jsym, amplitude) )
        end
        cs_Cou       = 0.;   cs_Bab = 0.
        crossSection = EmProperty(cs_Cou, cs_Bab)
        newLine      = RayleighCompton.Line(line.initialLevel, line.finalLevel, line.inOmega, line.outOmega, crossSection, newChannels)
        return( newLine )
    end

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

        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `RayleighCompton.computeReducedAmplitude(K::AngularJ64, finalLevel::Level, multipole2::EmMultipole, omega2::Float64, Jsym::LevelSymmetry,
                                             omega::Float64, multipole1::EmMultipole, omega1::Float64, initialLevel::Level, gauge::EmGauge,
                                             grid::Radial.Grid, greenChannels::Array{AtomicState.GreenChannel,1})`
        ... to compute the reduced amplitudes U^(K, Compton) as well as all properties of the given line;
            a line::RayleighCompton.Line is returned for which the reduced amplitudes and properties are now evaluated.
    """
    function  computeReducedAmplitude(K::AngularJ64, finalLevel::Level, multipole2::EmMultipole, omega2::Float64, Jsym::LevelSymmetry,
                                      omega::Float64, multipole1::EmMultipole, omega1::Float64, initialLevel::Level, gauge::EmGauge,
                                      grid::Radial.Grid, greenChannels::Array{AtomicState.GreenChannel,1})
        amplitude = 0.0im
        return( amplitude )
    end

    """
    `RayleighCompton.determineChannels(finalLevel::Level, initialLevel::Level, inOmega::Float64, outOmega::Float64, settings::RayleighCompton.Settings)`
        ... to determine a list of Rayleigh-Compton ReducedChannelChannel's for a transitions from the initial to final level and by taking
            into account the particular settings of for this computation; an Array{RayleighCompton.ReducedChannelChannel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, inOmega::Float64, outOmega::Float64, settings::RayleighCompton.Settings)
        channels = RayleighCompton.ReducedChannel[];
        symi = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        for  mp1 in settings.multipoles
            for  mp2 in settings.multipoles
                symmetries  = AngularMomentum.allowedTotalSymmetries(symf, mp2, mp1, symi)
                Klist       = oplus(symf.J, symi.J)
                for  symn in symmetries
                    for  gauge in settings.gauges
                        # Include further restrictions if appropriate
                        if     string(mp1)[1] == 'E' || string(mp2)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                            for K in Klist  push!(channels, ReducedChannel(K, mp1, mp2, Basics.Coulomb,   inOmega, outOmega, inOmega, symn, 0.) )  end
                        elseif string(mp1)[1] == 'E' || string(mp2)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                            for K in Klist  push!(channels, ReducedChannel(K, mp1, mp2, Basics.Babushkin, inOmega, outOmega, inOmega, symn, 0.) )  end
                        elseif string(mp1)[1] == 'M' && string(mp2)[1] == 'M'
                            for K in Klist  push!(channels, ReducedChannel(K, mp1, mp2, Basics.Magnetic,  inOmega, outOmega, inOmega, symn, 0.) )  end
                        end
                    end
                end
            end
        end
        # Unify channels
        return( channels )
    end

    """
    `RayleighCompton.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::RayleighCompton.Settings)`
        ... to determine a list of RayleighCompton Line's for elastic scattering or Raman-type 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{RayleighCompton.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::RayleighCompton.Settings)
        lines = RayleighCompton.Line[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    for  omega in settings.photonEnergies
                        inOmega  = Defaults.convertUnits("energy: to atomic", omega)
                        outOmega = inOmega - (fLevel.energy - iLevel.energy);       if  outOmega <= 0.    continue   end
                        channels = RayleighCompton.determineChannels(fLevel, iLevel, inOmega, outOmega, settings)
                        if   length(channels) == 0   continue   end
                        push!( lines, RayleighCompton.Line(iLevel, fLevel, inOmega, outOmega, EmProperty(0., 0.), channels) )
                    end
                end
            end
        end
        return( lines )
    end

    """
    `RayleighCompton.displayAngular(stream::IO, lines::Array{RayleighCompton.Line,1}, settings::RayleighCompton.Settings)`
        ... to display a list of lines and the angle-differential cross sections for the selected solid angles.
            A neat table of all selected transitions, energies and solid angles is printed but nothing is returned otherwise.
    """
    function  displayAngular(stream::IO, lines::Array{RayleighCompton.Line,1}, settings::RayleighCompton.Settings)
        nx = 127
        println(stream, " ")
        println(stream, "  Angular distribution of the selected Rayleigh-Compton scattering lines at seleced solid angles: ... to be adapted !!")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "   ";   sb = "   ";   sc = "  "
        sa = sa * TableStrings.center(16, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(16, "i--J^P--f"; na=7);                         sb = sb * TableStrings.hBlank(23)
        sa = sa * TableStrings.center(20, "i -- Energy -- f"; na=3);
        sb = sb * TableStrings.center(20, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(24, "omega-in   omega-out"; na=2);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=1)
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=6)
        println(stream, sa[1:end]);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy)) * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.finalLevel.energy))   * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.inOmega))      * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.outOmega))     * "       "
            println(stream, sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `RayleighCompton.displayCrossSections(stream::IO, lines::Array{RayleighCompton.Line,1}, settings::RayleighCompton.Settings)`
        ... to display a list of lines and their cross sections as they were selected by the settings. A neat table of all selected
            transitions, energies and cross sections is printed but nothing is returned otherwise.
    """
    function  displayCrossSections(stream::IO, lines::Array{RayleighCompton.Line,1}, settings::RayleighCompton.Settings)
        nx = 127
        println(stream, " ")
        println(stream, "  Cross sections of the selected Rayleigh-Compton scattering lines:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "   ";   sb = "   ";   sc = "  "
        sa = sa * TableStrings.center(16, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(16, "i--J^P--f"; na=7);                         sb = sb * TableStrings.hBlank(23)
        sa = sa * TableStrings.center(20, "i -- Energy -- f"; na=3);
        sb = sb * TableStrings.center(20, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(24, "omega-in   omega-out"; na=2);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=1)
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=6)
        sa = sa * TableStrings.center(28, "Cou -- cross section -- Bab"; na=0);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("cross section"); na=3)
        sb = sb * TableStrings.center(10, TableStrings.inUnits("cross section"); na=4)
        println(stream, sa[1:end]);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy)) * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.finalLevel.energy))   * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.inOmega))      * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.outOmega))     * "       "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Coulomb))       * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Babushkin))     * "   "
            println(stream, sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `RayleighCompton.displayLines(stream::IO, lines::Array{RayleighCompton.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(stream::IO, lines::Array{RayleighCompton.Line,1})
        nx = 157
        println(stream, " ")
        println(stream, "  Selected Rayleigh-Compton scacttering lines and (reduced) channels:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "   ";   sb = "   ";   sc = "  "
        sa = sa * TableStrings.center(16, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(16, "i--J^P--f"; na=7);                         sb = sb * TableStrings.hBlank(23)
        sa = sa * TableStrings.center(20, "i -- Energy -- f"; na=3);
        sb = sb * TableStrings.center(20, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(24, "omega-in   omega-out"; na=0);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=1)
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(11, "No channels"; na=2);                       sb = sb * TableStrings.hBlank(11)
        sa = sa * TableStrings.flushleft(60, "List of multipoles & intermediate level symmetries";   na=1)
        sb = sb * TableStrings.flushleft(60, "(K-rank, multipole_1, Jsym, multipole_2, gauge), ..."; na=1)
        println(stream, sa[1:end]);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy)) * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.finalLevel.energy))   * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.inOmega))      * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.outOmega))     * "   "
            sc = "         " * string(length(line.channels)) * "      "
            sa = sa * sc[end-12:end]
            mpGaugeList = Tuple{AngularJ64, Basics.EmMultipole, LevelSymmetry, Basics.EmMultipole, Basics.EmGauge}[]
            for  channel in  line.channels
                push!( mpGaugeList, (channel.K, channel.multipole1, channel.Jsym, channel.multipole2, channel.gauge) )
            end
            wa = TableStrings.twoPhotonGaugeTupels(60, mpGaugeList)
            if  length(wa) > 0    sb = sa * wa[1];    println( sb )    end
            for  i = 1:length(wa)
                sb = TableStrings.hBlank( length(sa) );    sb = sb * wa[i];    println(stream, sb )
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `RayleighCompton.displayStokesParameters(stream::IO, lines::Array{RayleighCompton.Line,1}, settings::RayleighCompton.Settings)`
        ... to display a list of lines and the Stokes parameters at the selected solid angles.
            A neat table of all selected transitions, energies and solid angles is printed but nothing is returned otherwise.
    """
    function  displayStokesParameters(stream::IO, lines::Array{RayleighCompton.Line,1}, settings::RayleighCompton.Settings)
        nx = 127
        println(stream, " ")
        println(stream, "  StokesParameters of the selected Rayleigh-Compton scattering lines at seleced solid angles: ... to be adapted !!")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "   ";   sb = "   ";   sc = "  "
        sa = sa * TableStrings.center(16, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(16, "i--J^P--f"; na=7);                         sb = sb * TableStrings.hBlank(23)
        sa = sa * TableStrings.center(20, "i -- Energy -- f"; na=3);
        sb = sb * TableStrings.center(20, TableStrings.inUnits("energy"); na=5)
        sa = sa * TableStrings.center(24, "omega-in   omega-out"; na=2);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=1)
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=6)
        println(stream, sa[1:end]);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.initialLevel.energy)) * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.finalLevel.energy))   * "   "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.inOmega))      * " "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.outOmega))     * "       "
            println(stream, sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

end # module


text:
import LightGraphs.LinAlg: adjacency_matrix

"""
    adjacency_matrix(g[, T=Int; dir=:out])

Return a sparse adjacency matrix for a graph, indexed by `[u, v]`
vertices. Non-zero values indicate an edge between `u` and `v`. Users may
override the default data type (`Int`) and specify an optional direction.

### Optional Arguments
`dir=:out`: `:in`, `:out`, or `:both` are currently supported.

### Implementation Notes
This function is optimized for speed and directly manipulates CSC sparse matrix fields.
"""
function adjacency_matrix(g::Graph, T::DataType=Float64; dir::Symbol=:out)
   n_v = nv(g)
   nz = ne(g)
   colpt = ones(Int64, n_v + 1)

   rowval = sizehint!(Vector{Int64}(), nz)
   weights   = sizehint!(Vector{T}(), nz)
    for j in 1:n_v  # this is by column, not by row.
       wgts = sizehint!(Vector{T}(), g.vertices[j].number_of_edges)
       dsts = sizehint!(Vector{Int64}(), g.vertices[j].number_of_edges)
        for e in g.vertices[j]
            push!(wgts,e.weight)
            push!(dsts,target_vertex(e,g.vertices[j]).id)
        end
        colpt[j + 1] = colpt[j] + length(dsts)

       dsts_indices = sortperm(dsts)
        append!(rowval, dsts[dsts_indices])
        append!(weights, wgts[dsts_indices])
    end
    return SparseMatrixCSC(n_v, n_v, colpt, rowval, weights)

end

ext: jl


text:
# using Ipaper
# using Test
# @subset(dt, y == 2)

# con = "x == 1 & y == 2 | z == 1"
# con_dt_transform(con; dname="df")
word = "[[:alnum:]_\\.\\[\\]]+"
pattern_lgl_left = "($word)( *)(?=[=<>≤≥\\!])" # logical operation
pattern_lgl_right = "(?<=[=<>≤≥])( *)($word)" # logical operation
pattern_op = r"([=<>≤≥\!]+)"

## Prepare for @subset
@testset "pattern key" begin
  # key = str_extract(x, pattern_lgl_left)
  @test str_extract_strip("x == 1", pattern_lgl_left) == "x"
  @test str_extract_strip("x != 1", pattern_lgl_left) == "x"
  @test str_extract_strip("x >= 1", pattern_lgl_left) == "x"
  @test str_extract_strip("x <= 1", pattern_lgl_left) == "x"
  @test str_extract_strip("x > 1", pattern_lgl_left) == "x"
  @test str_extract_strip("x < 1", pattern_lgl_left) == "x"
  # val = str_extract("x == 1", r"(?<=(==)|>=).*")
end

@testset "pattern val" begin
  # pattern_lgl_right = r"(?<===|\!=|>=|<=|>|≥|<|≤).*" # logical operation
  # key = str_extract(x, pattern_lgl_left)
  @test str_extract_strip("x == 1", pattern_lgl_right) == "1"
  @test str_extract_strip("x >= 1", pattern_lgl_right) == "1"
  @test str_extract_strip("x <= 1", pattern_lgl_right) == "1"
  @test str_extract_strip("x > 1", pattern_lgl_right) == "1"
  @test str_extract_strip("x < 1", pattern_lgl_right) == "1"
  @test str_extract_strip("x != 1", pattern_lgl_right) == "1"
end

@testset "pattern op" begin
  xs = [
    "x >= 1",
    "x == 1",
    "x != 1",
    "x > 1",
    "x < 1"]
  op = str_extract(xs, r"([=<>≤≥\!]+)")
  @test length(op) == length(xs)
end


text:
using BayesianSumProductNetworks
using StatsBase
using StatsFuns
using AxisArrays
using LinearAlgebra
using Random
using DelimitedFiles
using JLD2
using FileIO
using ArgParse
using Logging
using SparseArrays
using Dates
using HDF5
using Serialization

include("helper.jl")

s = ArgParseSettings()
@add_arg_table s begin
    "--alpha"
        help = "Concentration parameter of sum nodes"
        default = 1.0
        arg_type = Float64
    "--beta"
        help = "Concentration parameter of product nodes"
        default = 10.0
        arg_type = Float64
    "--gamma"
        help = "Concentration parameter of leaf distributions"
        default = 0.1
        arg_type = Float64
    "--Ksum"
        help = "Number of children under a sum node"
        default = 2
        arg_type = Int
    "--Kprod"
        help = "Number of children under a product node"
        default = 2
        arg_type = Int
    "--Kleaf"
        help = "Number of leaves in an atomic region"
        default = 2
        arg_type = Int
    "--J"
        help = "Number of partitions."
        default = 2
        arg_type = Int
    "--L"
        help = "Number of layers."
        default = 2
        arg_type = Int
    "--seed"
        help = "Random seed = seed * Int(first(dataset)) * chain."
        default = 42
        arg_type = Int
    "--chain"
        help = "Chain"
        default = 1
        arg_type = Int
    "--maxcategories"
        default = 20
        arg_type = Int
    "--numsamples"
        arg_type = Int
        default = 1000
        help = "Number of MCMC samples we want to store."
    "dataset"
        help = "Dataset name, assuming naming convention of 20 terrible datasets."
        required = true
        arg_type = String
    "datafolder"
        help = "Path to folder containing datasets."
        required = true
        arg_type = String
    "outputfolder"
        help = "Path to folder in which output files should be stored."
        required = true
        arg_type = String
    "burnin"
        help = "Number of burnin iterations."
        required = true
        arg_type = Int
    "iterations"
        help = "Number of iterations."
        required = true
        arg_type = Int
end

args = parse_args(s)

# Load the dataset
dataset = args["dataset"]
datafolder = args["datafolder"]

# results dir
resultsdir = args["outputfolder"]
if !isdir(resultsdir)
    mkpath(resultsdir)
end

@info "Processing $dataset in $datafolder"
@info "Storing results in $resultsdir"

donefile = joinpath(resultsdir, "done")
if isfile(donefile)
    @warn "Sampling already done."
    exit()
end

# check if results are locked and exit if this is the case
lockfile = joinpath(resultsdir, "filelock.lck")

# Output files.
modelfile = joinpath(resultsdir, "bspn_model.jld2")
burninfile = joinpath(resultsdir, "bspn_burnin.jld2")
llhfile = joinpath(resultsdir, "llhfile.jld2")
tracefile = joinpath(resultsdir, "tracefile.jld2")
iterfile = joinpath(resultsdir, "iteration")

# Check if we should continue the sampling...
if (isfile(modelfile) || isfile(burninfile)) && isfile(iterfile)
    d = now() - Dates.unix2datetime(stat(iterfile).mtime)
    if round(d, Minute) > Minute(60)
        # remove lock file
        if isfile(lockfile)
            rm(lockfile)
        end
    end
end

# Check if we should continue the sampling...
if isfile(lockfile)
    @info "MCMC sampling is locked, exiting."
    exit(-2)
else

    t = rand(1:15) + rand(1:15)
    sleep(t)

    if isfile(lockfile)
        @info "MCMC sampling is locked, exiting."
        exit(-2)
    end

    touch(lockfile)
end

# Write out the configuration of the run.
save(joinpath(resultsdir, "configuration.jld2"), Dict("configuration" => args))

# Start script.
open(joinpath(resultsdir, "log.txt"), "w+") do io
    logger = SimpleLogger(io)
    with_logger(logger) do
        if Threads.nthreads() == 1
            @info "Run export JULIA_NUM_THREADS=nproc before running julia to use multithreading."
        else
            @info "Using $(Threads.nthreads()) threads."

            # Test if all threads are available.
            thrsa = zeros(Threads.nthreads())
            Threads.@threads for t in 1:Threads.nthreads()
                thrsa[t] = Threads.threadid()
            end

            avThreads = length(unique(thrsa[thrsa .> 0]))
            @assert avThreads == Threads.nthreads() "Only $(avThreads) threads available."
        end

        @info "Loading dataset: $dataset stored under: $datafolder"

        x_train = map(Float64, h5read(joinpath(datafolder, dataset*"PP", "train.hdf5"), "train"))
        x_train = transpose(x_train)

        x_valid = map(Float64, h5read(joinpath(datafolder, dataset*"PP", "valid.hdf5"), "valid"))
        x_train = cat(x_train, transpose(x_valid), dims = 1)

        x_test = map(Float64, h5read(joinpath(datafolder, dataset*"PP", "test.hdf5"), "test"))
        x_test = transpose(x_test)

        # Set the number of iterations, burnin, etc.
        iterations = args["iterations"]
        burnin = args["burnin"]

        # ######### #
        # MCMC code #
        # ######### #
        ks = min(iterations, args["numsamples"])
        @assert isinteger(iterations / ks)
        thinning = Int(iterations / ks)

        (N_train, D) = size(x_train)
        (N_test, _) = size(x_test)

        @info "Processing dataset of size N=$N_train, D=$D."

        # Preprocessing
        dids = findall(map(d -> all(isinteger, filter(x -> !isnan(x), x_train[:,d])), 1:D))
        dids = filter(d -> length(filter(x -> !isnan(x), unique(x_train[:,d]))) < args["maxcategories"], 1:D)
        for d in dids
            ix = findall(!isnan, x_train[:,d])
            l = unique(x_train[ix,d])

            x_train[ix,d] = map(x -> findfirst(l .== x), x_train[ix,d])

            ix = findall(!isnan, x_test[:,d])
            x_test[ix,d] = map(x -> x ∈ l ? findfirst(l .== x) : length(l)+1, x_test[ix,d])

            if length(l) == 2
                # Make binary.
                x_train .-= 1
                x_test .-= 1
            end
        end

        priors = Vector{Vector{Distribution}}(undef, D)
        likelihoods = Vector{Vector{Distribution}}(undef, D)
        sstats = Vector{Vector{AbstractSufficientStats}}(undef, D)

        for d in 1:D
            priors[d] = Vector{Distribution}()
            likelihoods[d] = Vector{Distribution}()
            sstats[d] = Vector{Distribution}()
            ix = findall(!isnan, x_train[:,d])

            if d ∈ dids # is discrete ?
                if minimum(x_train[ix,d]) == 0
                    # Binary data
                    push!(priors[d], Beta(1/2, 1/2))
                    push!(likelihoods[d], Bernoulli())
                    push!(sstats[d], BetaSufficientStats())
                else
                    # Categorical data
                    l = unique(x_train[:,d])
                    push!(priors[d], Dirichlet(length(l), 0.1))
                    push!(likelihoods[d], Categorical(rand(last(priors[d]))))
                    push!(sstats[d], DirichletSufficientStats(length(l)))
                end

                # Poisson data
                push!(priors[d], Gamma())
                push!(likelihoods[d], Poisson())
                push!(sstats[d], GammaSufficienStats())
            else # continuous data
                if all((x) -> x >= 0, x_train[ix,d])
                    # Positive continuous data
                    push!(priors[d], Gamma())
                    push!(likelihoods[d], Exponential())
                    push!(sstats[d], GammaSufficienStats())
                end

                # Normal distributed data
                push!(priors[d], NormalInverseGamma(mean(x_train[ix,d]), 1.0/std(x_train[ix,d]), 2.0, 3.0))
                push!(likelihoods[d], Normal())
                push!(sstats[d], NormalInvGammaSufficientStats())
            end
        end

        # Set the random seed.
        rnd_seed = args["seed"]+args["chain"]
        rng = Random.seed!(rnd_seed)

        # Construct a second network which we use for sampling.
        @info "Constructing spn for sampling."
        spn = templateRegion(args["alpha"],
                             args["beta"],
                             args["gamma"],
                             priors,
                             likelihoods,
                             sstats,
                             N_train, D,
                             args["Ksum"], args["Kprod"],
                             args["J"], args["Kleaf"],
                             1, args["L"], root = true)

        @info "Starting MCMC for $iterations iterations."

        startIter = 1
        storePredictions = true
        storeTrace = false

        # Initialize SPN.
        if isfile(modelfile)
            # Resume.
            spn = open(modelfile, "r") do fin
                deserialize(fin)
            end
            r_ = read(iterfile, Int)
            startIter = 1 + r_
            @info "Continue MCMC."
        elseif isfile(burninfile)
            spn = open(burninfile, "r") do fin
                deserialize(fin)
            end
            @info "Continue burnin or start using previously converged chain, not storing trace."
            r_ = read(iterfile, Int)
            startIter = 1 + r_
        else
            @warn "Starting from random position, not storing predictions!"
            randomscopes!(spn, x_train)
        end

        for iteration in startIter:(burnin + iterations)
            dateprint = Dates.format(now(), "HH:MM:SS")

            storeTrace = iteration <= burnin
            storePredictions = iteration > burnin

            @info "[$(dateprint)] iteration = $(iteration) "

            ancestralsampling!(spn, x_train)
            gibbssamplescopes!(spn, x_train)

            if ((iteration-burnin) % thinning) == 0
                # Serialize model.
                fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                open(fname, "w") do fout
                    serialize(fout, spn)
                end
                if iteration > burnin
                    mv(fname, modelfile, force=true)
                else
                    mv(fname, burninfile, force=true)
                end

                if storeTrace
                    # Store parameters of root.
                    Nparam = length(spn.logweights[:])
                    θ = isfile(tracefile) ? load(tracefile, "weights") : zeros(0,Nparam)
                    fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                    save(fname, "weights", vcat(θ, spn.logweights[:]'))
                    mv(fname, tracefile, force=true)
                end

                if storePredictions
                    # Evaluate network.
                    lp_test = logpdf(spn, x_test)
                    lp_train = logpdf(spn, x_train)

                    b_ = isfile(llhfile)
                    llh_test = b_ ? load(llhfile, "test") : zeros(0,N_test)
                    llh_train = b_ ? load(llhfile, "train") : zeros(0,N_train)
                    fname = joinpath(resultsdir, "tmp_"*randstring()*".jld2")
                    save(fname, "train", vcat(llh_train, lp_train'), "test", vcat(llh_test, lp_test'))
                    mv(fname, llhfile, force=true)
                end

                write(iterfile, iteration)
            end
        end
        # Save sampling results to disk.
        @info "Finished MCMC"

        # Create done file and release lock.
        touch(donefile)
        rm(lockfile)

        tmp_files = filter(f -> isfile(joinpath(resultsdir,f)) && startswith(f, "tmp_"), readdir(resultsdir))
        for f in tmp_files
            rm(joinpath(resultsdir,f))
        end
    end
end

max_stars_repo_name: trappmartin/BayesianSumProductNetworks

max_issues_repo_name: trappmartin/BayesianSumProductNetworks

max_forks_repo_name: trappmartin/BayesianSumProductNetworks


text:
ztype(op::libgb.GrB_UnaryOp) = tojuliatype(ptrtogbtype[libgb.GxB_UnaryOp_ztype(op)])
xtype(op::libgb.GrB_UnaryOp) = tojuliatype(ptrtogbtype[libgb.GxB_UnaryOp_xtype(op)])
Base.show(io::IO, ::MIME"text/plain", u::libgb.GrB_UnaryOp) = gxbprint(io, u)

xtype(op::libgb.GrB_BinaryOp) = tojuliatype(ptrtogbtype[libgb.GxB_BinaryOp_xtype(op)])
ytype(op::libgb.GrB_BinaryOp) = tojuliatype(ptrtogbtype[libgb.GxB_BinaryOp_ytype(op)])
ztype(op::libgb.GrB_BinaryOp) = tojuliatype(ptrtogbtype[libgb.GxB_BinaryOp_ztype(op)])
Base.show(io::IO, ::MIME"text/plain", u::libgb.GrB_BinaryOp) = gxbprint(io, u)

operator(monoid::libgb.GrB_Monoid) = libgb.GxB_Monoid_operator(monoid)
xtype(monoid::libgb.GrB_Monoid) = xtype(operator(monoid))
ytype(monoid::libgb.GrB_Monoid) = ytype(operator(monoid))
ztype(monoid::libgb.GrB_Monoid) = ztype(operator(monoid))
Base.show(io::IO, ::MIME"text/plain", m::libgb.GrB_Monoid) = gxbprint(io, m)

multiplyop(rig::libgb.GrB_Semiring) = libgb.GxB_Semiring_multiply(rig)
addop(rig::libgb.GrB_Semiring) = libgb.GxB_Semiring_add(rig)
xtype(rig::libgb.GrB_Semiring) = xtype(multiplyop(rig))
ytype(rig::libgb.GrB_Semiring) = ytype(multiplyop(rig))
ztype(rig::libgb.GrB_Semiring) = ztype(addop(rig))
Base.show(io::IO, ::MIME"text/plain", s::libgb.GrB_Semiring) = gxbprint(io, s)


text:
function bisection_refine!(u_A, u_B, bsn_nfo, integ, basin_A, basin_B, tol)

    # shortcut functions
    function get_col(u0)
        a = get_color_point!(bsn_nfo, integ, u0)
        return iseven(a) ? Int(a/2) : Int((a-1)/2)
    end

    function get_col_next(u0)
        bsn_nfo.reinit_f!(integ,u0)
        bsn_nfo.iter_f!(integ)
        a=get_col(bsn_nfo.get_u(integ))
        return  a
    end

    # directing vector
    Δu = u_B - u_A

    # Midpoint method
    u_A_next=u_A
    u_B_next=u_B
    dist= norm(u_A_next-u_B_next)

    if dist < tol

        if get_col_next(u_A_next) ∉ basin_A
            # This basin has switch we try to perturb u_A_next until we get back to the initial basin.
            # We guarantee that the next iteration will stay also in the good basin
            eps=1e-8
            u_A_rnd = u_A_next + eps*randn(2)
            while get_col_next(u_A_rnd) ∉ basin_A
                eps += eps
                @show u_A_rnd = u_A_next + eps*randn(2)
            end
            u_A_next=u_A_rnd
        end

        if get_col_next(u_B_next) ∉ basin_B
            # This basin has switch we try to perturb u_B_next until we get back to the initial basin.
            eps=1e-8
            u_B_rnd = u_B_next + eps*randn(2)
            while get_col_next(u_B_rnd) ∉ basin_B
                eps += eps
                @show u_B_rnd = u_B_next + eps*randn(2)
            end
            u_B_next=u_B_rnd
        end

        return u_A_next,u_B_next

    else
        # refinement with midpoint method.
        while dist > tol

            mV= Δu/2

            if get_col(u_A_next+mV) ∈ basin_A
                u_A_next = u_A_next+mV
            elseif get_col(u_B_next-mV) ∈ basin_B
                u_B_next = u_B_next-mV
            else
                @show get_col(u_B_next-mV), get_col(u_A_next+mV)
                break
            end

            if get_col_next(u_A_next) ∉ basin_A
                # This basin has switch we try to perturb u_A_next until we get back to the initial basin.
                eps=1e-8
                u_A_rnd = u_A_next + eps*randn(2)
                while get_col_next(u_A_rnd) ∉ basin_A
                    eps += eps
                    @show u_A_rnd = u_A_next + eps*randn(2)
                end
                u_A_next=u_A_rnd
            end

            if get_col_next(u_B_next) ∉ basin_B
                # This basin has switch we try to perturb u_B_next until we get back to the initial basin.
                eps=1e-8
                u_B_rnd = u_B_next + eps*randn(2)
                while get_col_next(u_B_rnd) ∉ basin_B
                    eps += eps
                    @show u_B_rnd = u_B_next + eps*randn(2)
                end
                u_B_next=u_B_rnd
            end

            Δu = u_B_next - u_A_next
            dist= norm(u_A_next-u_B_next)
        end
    end

    return u_A_next,u_B_next
end

"""
    compute_saddle(bsn_nfo::BasinInfo, integ; max_iter=10)
The algorithm the saddle that lies in  a boundary of the basin of attraction of a dynamical system. When the boundary is fractal, this
set is known as the chaotic saddle and is responsible for the transient dynamics of the system. The saddle is computed with the Saddle-Straddle
algorithm. It is necessary to define two `generalized basin`, that is we must separate the basin into two sets (see also keyword arguments).

[H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 2012]

## Arguments
* `bsn_nfo` : structure that holds the information of the basin as well as the map function. This structure is set when the basin is first computed with `basin_stroboscopic_map` or `basin_poincare_map`.
* `integ` : the matrix containing the information of the basin.
* `bas_A` : vector with the indices of the attractors that will represent the generalized basin A
* `bas_B` : vector with the indices of the attractors that will represent the generalized basin B. Notice that `bas_A ∪ bas_B = [1:N]` and `bas_A ∩ bas_B = ∅`

## Keyword arguments
* `N` : number of points of the saddle to compute
* `init_tol`: length of the straddle segment

## Example
```
# compute 1000 points between generalized A = [1] and basin B = [2,3]
sa,sb = compute_saddle(bsn, integ_df, [1], [2,3],1000)
# sa is the "left" set and sb the right "set"
```

"""
function compute_saddle(integ, bsn_nfo::BasinInfo, bas_A, bas_B; N=100, init_tol = 1e-6)

    # shortcut functions
    function get_col(u0)
        a = get_color_point!(bsn_nfo, integ, u0)
        return iseven(a) ? Int(a/2) : Int((a-1)/2)
    end

    Na = bsn_nfo.Na
    # basic check
    if (Set(bas_A) ∪ Set(bas_B)) != Set(collect(1:Na))
        @error "Generalized basins are not well defined"
        return nothing,nothing
    end

    if length(Set(bas_A) ∩ Set(bas_B)) > 0
        @error "Generalized basins are not well defined"
        return nothing, nothing
    end

    tol = init_tol
    xg = bsn_nfo.xg; yg = bsn_nfo.yg; # aliases

    # Set initial condition near attractors of the selected basins.
    u_A = u_B = [0., 0.]
    v = bsn_nfo.attractors[bas_A[1]*2]
    u_A = v[1]
    v = bsn_nfo.attractors[bas_B[1]*2]
    u_B = v[1]

    saddle_series_A = Array{typeof(u_A),1}(undef, N)
    saddle_series_B = Array{typeof(u_A),1}(undef, N)
    u_A_it = u_B_it = [0.,0.]
    Ttr= 20 # skip 20 points
    switch_bas = 0 # keep track of switching between basin if it happens too much the basin is problematic
    k=1;
    while k<= N

        # Initial bisection
        u_A_r,u_B_r = bisection_refine!(u_A, u_B, bsn_nfo, integ, bas_A, bas_B, tol)

        dist = norm(u_A_r-u_B_r)
        bsn_nfo.reinit_f!(integ,u_A_r)
        bsn_nfo.iter_f!(integ)
        u_A_it = deepcopy(bsn_nfo.get_u(integ)) # for some reason I have to deepcopy this vector
        ca=get_col(u_A_it)

        bsn_nfo.reinit_f!(integ,u_B_r)
        bsn_nfo.iter_f!(integ)
        u_B_it = deepcopy(bsn_nfo.get_u(integ))
        cb=get_col(u_B_it)

        if ca ∉ bas_A || cb ∉ bas_B
            # one of the basins switched attractor, we retry refinement but with bigger tolerance.
            println("Basin have switched this is unusual, basin A:", ca, " basin B:  ", cb)
            switch_bas += 1
            tol = tol*2
            Ttr= 20; # skip 20 points
        else
            if Ttr <= 0 &&  dist < tol
                saddle_series_A[k]=u_A
                saddle_series_B[k]=u_B
                k +=1
                switch_bas = 0
            else
                Ttr -= 1;
            end
            u_A=u_A_it;
            u_B=u_B_it;
            #@show u_A,u_B
            tol = init_tol
        end

        if switch_bas > 20
            @warn "Problem with the basin: it may be a riddled boundary or the initial resolution of the basin is too low."
            break
        end
    end

    return saddle_series_A, saddle_series_B
end


text:
# This file was generated, do not modify it. # hide
using BenchmarkTools

function array_work(f, μ, σ, x, y)
    @inbounds for i in eachindex(x)
        y[i] = logdensity(f(μ[i], σ[i]), x[i])
    end
end

time_normal(f) = @belapsed $array_work($f, $μ, $σ, $x, $y)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
#FIXME: PyObject not found when pyscf only
import PySCF: scf_rhf, mymp2, short_mp2 #, pyscf

# Workaround: repeat the imports.
using PyCall: pyimport
pyscf = pyimport("pyscf")
mp = pyimport("pyscf.mp")   # Had to import mp alone ??!

@testset "MP₂: H₂O" begin

   mol = pyscf.gto.Mole()
   mol.build(verbose = 0,
             basis = "sto-3g",
             atom = """8 0  0.    0
                       1 0 -0.757 0.587
                       1 0  0.757 0.587""",
            )

   eri = mol.intor("cint2e_sph", aosym="s8")
   rhf      = scf_rhf(mol,verbose=false)
   e_rhf_me, e, C = rhf.Emol, rhf.e, rhf.C
   mf = pyscf.scf.RHF(mol) # set the model
   mf.kernel()             # run RHF (populate mf with data)

   e_mp2_me1   = mymp2(e,C,eri,mol,alg=:ao2mo_noddy) # Naive loops ao2mo with my MP₂
   e_mp2_me2   = mymp2(e,C,eri,mol,alg=:ao2mo_smart) # Smart loops ao2mo with my MP₂
   e_mp2_me3   = short_mp2(e,C,eri,mol)              # einsum      ao2mo with my MP₂
   e_mp2_pyscf = mp.MP2(mf).kernel()[1]              # pure PySCF MP₂ solution

   @test isapprox(e_mp2_me1, e_mp2_pyscf, atol=1e-7)
   @test isapprox(e_mp2_me2, e_mp2_pyscf, atol=1e-7)
   @test isapprox(e_mp2_me3, e_mp2_pyscf, atol=1e-7)
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

MyGame.Example.eval(quote

FlatBuffers.@STRUCT struct Vec3
    x::Float32
    y::Float32
    z::Float32
    test1::Float64
    test2::Color
    test3::Test
end
FlatBuffers.@ALIGN(Vec3, 8)

Vec3(buf::AbstractVector{UInt8}) = FlatBuffers.read(Vec3, buf)
Vec3(io::IO) = FlatBuffers.deserialize(io, Vec3)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
taxonomy_path(msa::MSA) = path(msa) * "_taxonomy"

function load_with_taxonomy(msa::MSA)
    if !isfile(taxonomy_path(msa))
        df = load(msa)
        tax = get_taxonomies(df.uniprot_id)
        dftax = DataFrame(uniprot_id = df.uniprot_id, taxonomy = tax)
        CSV.write(taxonomy_path(msa), dftax)
    end
    return CSV.read(taxonomy_path(msa), DataFrame)
end

function get_taxonomies(uniprot_ids::AbstractVector{String})
    idx = Dict(id => k for (k, id) in enumerate(uniprot_ids))
    tax = Vector{Union{Missing,String}}(missing, length(uniprot_ids))
    open(datadep"uniprot/uniprot.txt.gz") do io
        gzip = GzipDecompressorStream(io)
        for line in eachline(gzip)
            rows = split(line, "\t")
            id = rows[2]
            if haskey(idx, id)
                @assert length(rows) == 22
                tax[idx[id]] = rows[10]
            end
        end
    end
    return tax
end

#
# function uniprot_headers()
# 	headers = [
# 		(:uniprot_acc, String), (:uniprot_id, String),
# 		(:seq_version, Int), (:crc64, String), (:md5, String),
# 		(:description, String), (:evidence, Int), (:length, Int),
# 		(:species, String), (:taxonomy, String),
# 		(:is_fragment, Int), (:sequence, String),
# 		(:updated, String),	(:created, String),
# 		(:ncbi_taxid, Int),	(:ref_proteome, Int),
# 		(:complete_proteome, Int), (:treefam_acc, String),
# 		(:rp15, Int), (:rp35, Int), (:rp55, Int), (:rp75, Int)
# 	]
# end
#
# function uniprot_taxonomy_table()
# 	uniprot_headers = uniprot_csv_headers()
# 	uniprot_csv() do csv
# 		df = csv |> select(:uniprot_id, :taxonomy) |> DataFrame
# 	end
# 	return df
# end
#
# function uniprot_csv(f)
# 	io = open(datadep"uniprot/uniprot.txt.gz")
# 	try
# 		gzip = GzipDecompressorStream(io)
# 		csv = CSV.File(gzip; header=first.(headers), types=last.(headers))
# 		f(csv)
# 	finally
# 		close(io)
# 	end
# end


text:
import UncertainData: InterpolateAndBin

"""
    InterpolateBinTest(test::CausalityTest, interpolate_bin_resampling::InterpolateAndBin)
    InterpolateBinTest(test::CausalityTest, interpolate_bin_resampling::InterpolateAndBin, n::Int)

A causality test where the data is interpolated and then binned before the test is applied.

`n` controls the number of independent resamplings over which the test
is performed (if not provided, `n` is set to `1` by default).

## Examples

Record `N` points from the built-in `ar1_unidir` system, collect the 1st and 2nd
variables as `X` and `Y` and add some uncertainties to both the indices and the values.

```julia
N = 300
sys = ar1_unidir(c_xy = 0.8)
X, Y = example_uncertain_indexvalue_datasets(sys, N, (1, 2),
    d_xval = Uniform(0.001, 0.05), d_yval = Uniform(0.001, 0.05));
```

Say we want to interpolate both `X` and `Y` to some very fine grid
using linear interpolation, then bin to a coarser grid and summarise
each bin using the mean of the values in each bin.

```julia
# Define interpolation grid over the range of available index values
tmin = max(minimum(mean.(X.indices)), minimum(mean.(Y.indices)))
tmax = max(maximum(mean.(X.indices)), maximum(mean.(Y.indices)))
intp_grid = tmin:0.01:tmax

# Define binning grid
left_bin_edges = tmin:5:tmax

# Define the InterpolateAndBin instance
intp_bin = InterpolateAndBin(mean, left_bin_edges, Linear(), intp_grid, Flat(OnGrid()))
```
"""
struct InterpolateBinTest{
        CT <: CausalityTest,
        IB <: InterpolateAndBin{L} where L} <: PreprocessedDataCausalityTest{CT}
    test::CT
    interpolate_bin_resampling::IB
    n::Int
end

function InterpolateBinTest(test, interpolate_bin_resampling)
    InterpolateBinTest(test, interpolate_bin_resampling, 1)
end;

export InterpolateBinTest


text:
function param_count(model)
    ps = Flux.params(model)
    res = 0
    for p in keys(ps.params.dict)
        res += prod(size(p))
    end
    res
end

struct Flatten end

function (l::Flatten)(input)
    reshape(input, :, size(input, 4))
end

struct Sigmoid end
function (l::Sigmoid)(input)
    σ.(input)
end

struct ReLU end

function (l::ReLU)(input)
    relu.(input)
end

struct LeakyReLU end
function (l::LeakyReLU)(input)
    leakyrelu.(input)
end

struct Tanh end
function (l::Tanh)(input)
    tanh.(input)
end

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
module FDiffUtils

using ForwardDiff
using ForwardDiff: Dual
using Parameters: @with_kw
using StaticArrays: SVector

function vderiv2(f, x::Number = 0.0)
    return ForwardDiff.jacobian(
        y -> ForwardDiff.jacobian(z -> f(z[1]), y)[:, 1],
        SVector(x),
    )[:, 1]
end

function deriv2(f, x::Number = 0.0)
    return ForwardDiff.derivative(
        a -> ForwardDiff.derivative(f, a),
        x)
end

function deriv3(f, x::Number = 0.0)
    return ForwardDiff.derivative(
        a -> ForwardDiff.derivative(
            b -> ForwardDiff.derivative(f, b),
            a),
        x)
end

nondual(x::Dual) = nondual(x.value)
nondual(x::AbstractArray) = nondual.(x)
nondual(x) = x

_all(f, x::AbstractArray) = all(f, x)
_all(f, x::Dual) = _all(f, x.value) && all(_all.(Ref(f), x.partials))
_all(f, x) = f(x)
_any(f, x) = ! _all(!f, x)

all_isfinite(x) = _all(isfinite, x)
any_isnan(x) = _any(isnan, x)

@with_kw struct NotFinite <: Exception
    value
    expr = nothing
end

function Base.showerror(io::IO, err::NotFinite)
    if err.expr != nothing
        show(io, expr)
    else
        print(io, "Value")
    end
    println(io, " contains non-finite value(s):")
    show(IOContext(io, :compact => false), err.value)
    println(io)
end

macro assert_finite(expr)
    quote
        x = $(esc(expr))
        if ! all_isfinite(x)
            error(NotFinite(x, $(QuoteNode(expr))))
        end
        x
    end
end

end  # module


text:
module NASAPrecipitation

## Modules Used
using Logging
using NCDatasets
using Printf
using Statistics

## Reexporting exported functions within these modules
using Reexport
@reexport using Dates
@reexport using GeoRegions

import Base: download, show, read

## Exporting the following functions:
export
        IMERGFinalHH, IMERGLateHH, IMERGEarlyHH,
        IMERGFinalDY, IMERGLateDY, IMERGEarlyDY,
        IMERGMonthly,

        TRMM3Hourly, TRMM3HourlyNRT,
        TRMMDaily,   TRMMDailyNRT,
        TRMMMonthly,

        getIMERGlsm, getTRMMlsm, addNPDGeoRegions,

        download, read, npdfnc

## Abstract types
"""
    NASAPrecipitationDataset{ST<:AbstractString, DT<:TimeType}
        npdID :: ST
        lname :: ST
        doi   :: ST
        dtbeg :: DT
        dtend :: DT
        sroot :: ST
        hroot :: ST
        fpref :: ST
        fsuff :: ST
    end

Abstract supertype for NASA Precipitation datasets on NASA OPeNDAP Servers.

Fields:
* `npdID` : ID for the `NASAPrecipitationDataset`, used in determining containing folders and filenames of the NetCDF
* `lname` : The name describing the `NASAPrecipitationDataset`, used mostly in Logging
* `doi`   : The DOI identifier, to be saved into the NetCDF
* `dtbeg` : The start date (Y,M,D) of our download / analysis
* `dtend` : The end date (Y,M,D) of our download / analysis
* `sroot` : The directory in which to save our downloads and analysis files to
* `hroot` : The URL of the NASA's EOSDIS OPeNDAP server for which this dataset is stored
* `fpref` : The prefix component of the NetCDF files to be downloaded
* `fsuff` : The suffix component of the NetCDF files to be downloaded

Of these fields, only `dtbeg`, `dtend` and `sroot` are user-defined.  All other fields are predetermined depending on the type of NASA Precipitation Dataset called.
"""
abstract type NASAPrecipitationDataset end

"""
    IMERGDataset <: NASAPrecipitationDataset

Abstract supertype for GPM IMERG datasets on NASA OPeNDAP Servers, a subType of `NASAPrecipitationDataset`.
"""
abstract type IMERGDataset <: NASAPrecipitationDataset end

"""
    TRMMDataset <: NASAPrecipitationDataset

Abstract supertype for TRMM TMPA datasets on NASA OPeNDAP Servers, a subType of `NASAPrecipitationDataset`.
"""
abstract type TRMMDataset  <: NASAPrecipitationDataset end

modulelog() = "$(now()) - NASAPrecipitation.jl"

function __init__()
    @info "$(modulelog()) - Checking to see if GeoRegions required by NASAPrecipitation.jl have been added to the list of available GeoRegions"
    disable_logging(Logging.Warn)
	if !isGeoRegion("IMERG",throw=false) ||
	    !isGeoRegion("TRMM",throw=false) ||
	    !isGeoRegion("TRMMLSM",throw=false)
        disable_logging(Logging.Debug)
        @info "$(modulelog()) - At least one of the required three GeoRegions (IMERG, TRMM, TRMMLSM) has not been added, proceeding to add them again ..."
        disable_logging(Logging.Warn)
	    addGeoRegions(joinpath(@__DIR__,"NPDGeoRegions.txt"))
    else
        disable_logging(Logging.Debug)
        @info "$(modulelog()) - All of the required three GeoRegions (IMERG, TRMM, TRMMLSM) have been added"
        disable_logging(Logging.Warn)
	end
    disable_logging(Logging.Debug)
end

## Including Relevant Files

include("IMERG/halfhourly.jl")
include("IMERG/daily.jl")
include("IMERG/monthly.jl")
include("IMERG/landseamask.jl")

include("TRMM/3hourly.jl")
include("TRMM/daily.jl")
include("TRMM/monthly.jl")
include("TRMM/landseamask.jl")

include("downloads.jl")
include("save.jl")
include("read.jl")
include("filesystem.jl")
include("backend.jl")

end

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text:
using SafeTestsets

# Blocks
@safetestset "Blocks: math" begin include("Blocks/math.jl") end
@safetestset "Blocks: nonlinear" begin include("Blocks/nonlinear.jl") end
@safetestset "Blocks: continuous" begin include("Blocks/continuous.jl") end
@safetestset "Blocks: sources" begin include("Blocks/sources.jl") end

# Electrical
@safetestset "Analog Circuits" begin include("Electrical/analog.jl") end
#@safetestset "Digital Circuits" begin include("Electrical/digital.jl") end
@safetestset "RC Circuit Demo" begin include("demo.jl") end

# Thermal
@safetestset "Thermal Circuits" begin include("Thermal/thermal.jl") end
@safetestset "Thermal Demo" begin include("Thermal/demo.jl") end

# Magnetic
@safetestset "Magnetic" begin include("Magnetic/magnetic.jl") end

# Mechanical
@safetestset "Mechanical" begin include("Mechanical/rotational.jl") end


text:
function to_rna(dna)
    trans = ""
    for n in dna
        # print(n)
        if n == ""
            trans = trans
        elseif n == 'C'
            trans *= 'G'
        elseif n == 'G'
            trans *= 'C'
        elseif n == 'T'
            trans *= 'A'
        elseif n == 'A'
            trans *= 'U'
        else
            throw(ErrorException("bad input"))
        end
    end

    return trans
end

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
################################################################################
#  Copyright 2021, Tom Van Acker                                               #
################################################################################
# StochasticPowerModels.jl                                                     #
# An extention package of PowerModels.jl for Stochastic (Optimal) Power Flow   #
# See http://github.com/timmyfaraday/StochasticPowerModels.jl                  #
################################################################################

# bus voltage
"variable: `ve[i]` for `i` in `bus`es"
function variable_bus_voltage_expectation(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, aux_fix::Bool=false)
    ve = _PM.var(pm, nw)[:ve] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :bus)], base_name="$(nw)_ve",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :bus, i), "ve_start", 1.0)
    )

    if bounded
        for (i, bus) in _PM.ref(pm, nw, :bus)
            JuMP.set_lower_bound(ve[i], 0.0)
        end
    end

    if aux_fix
        JuMP.fix.(ve, 1.0; force=true)
    end

    report && _PM.sol_component_value(pm, nw, :bus, :ve, _PM.ids(pm, nw, :bus), ve)
end

"variable: `vv[i]` for `i` in `bus`es"
function variable_bus_voltage_variance(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, aux_fix::Bool=false)
    vv = _PM.var(pm, nw)[:vv] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :bus)], base_name="$(nw)_vv",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :bus, i), "vv_start", 1.0)
    )

    if bounded
        for (i, bus) in _PM.ref(pm, nw, :bus)
            JuMP.set_lower_bound(vv[i],  0.0)
        end
    end

    if aux_fix
        JuMP.fix.(vv, 1.0; force=true)
    end

    report && _PM.sol_component_value(pm, nw, :bus, :vv, _PM.ids(pm, nw, :bus), vv)
end

"variable: `vs[i]` for `i` in `bus`es"
function variable_bus_voltage_squared(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, aux_fix::Bool=false)
    vs = _PM.var(pm, nw)[:vs] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :bus)], base_name="$(nw)_vs",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :bus, i), "vs_start", 1.0)
    )

    if bounded
        for (i, bus) in _PM.ref(pm, nw, :bus)
            JuMP.set_lower_bound(vs[i], -2.0 * bus["vmax"]^2)
            JuMP.set_upper_bound(vs[i],  2.0 * bus["vmax"]^2)
        end
    end

    if aux_fix
        JuMP.fix.(vs, 1.0; force=true)
    end

    report && _PM.sol_component_value(pm, nw, :bus, :vs, _PM.ids(pm, nw, :bus), vs)
end

# load current
"variable: `crd[j]` for `j` in `load`"
function variable_load_current_real(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    crd = _PM.var(pm, nw)[:crd] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :load)], base_name="$(nw)_crd",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :load, i), "crd_start")
    )

    report && _PM.sol_component_value(pm, nw, :load, :crd, _PM.ids(pm, nw, :load), crd)
end

"variable: `cid[j]` for `j` in `load`"
function variable_load_current_imaginary(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    cid = _PM.var(pm, nw)[:cid] = JuMP.@variable(pm.model,
        [i in _PM.ids(pm, nw, :load)], base_name="$(nw)_cid",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :load, i), "cid_start")
    )

    report && _PM.sol_component_value(pm, nw, :load, :cid, _PM.ids(pm, nw, :load), cid)
end

# branch current
"variable: `cr[l,i,j]` for `(l,i,j)` in `arcs`"
function expression_variable_branch_current_real(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    cr = _PM.var(pm, nw)[:cr] = Dict()

    bus = _PM.ref(pm, nw, :bus)
    branch = _PM.ref(pm, nw, :branch)

    for (l,i,j) in _PM.ref(pm, nw, :arcs_from)
        b = branch[l]
        tm = b["tap"]
        tr, ti = _PM.calc_branch_t(b)
        g_sh_fr, b_sh_fr = b["g_fr"], b["b_fr"]
        g_sh_to, b_sh_to = b["g_to"], b["b_to"]

        vr_fr = _PM.var(pm, nw, :vr, i)
        vi_fr = _PM.var(pm, nw, :vi, i)

        vr_to = _PM.var(pm, nw, :vr, j)
        vi_to = _PM.var(pm, nw, :vi, j)

        csr_fr = _PM.var(pm, nw, :csr, l)
        csi_fr = _PM.var(pm, nw, :csi, l)

        cr[(l,i,j)] = (tr * csr_fr - ti * csi_fr + g_sh_fr * vr_fr - b_sh_fr * vi_fr) / tm^2
        cr[(l,j,i)] = -csr_fr + g_sh_to * vr_to - b_sh_to * vi_to

        # ub = Inf
        # if haskey(b, "rate_a")
        #     rate_fr = b["rate_a"]*b["tap"]
        #     rate_to = b["rate_a"]
        #     ub = max(rate_fr/bus[i]["vmin"], rate_to/bus[j]["vmin"])
        # end
        # if haskey(b, "c_rating_a")
        #     ub = b["c_rating_a"]
        # end

        # if !isinf(ub)
        #     JuMP.@constraint(pm.model, cr[(l,i,j)] >= -ub)
        #     JuMP.@constraint(pm.model, cr[(l,i,j)] <= ub)

        #     JuMP.@constraint(pm.model, cr[(l,j,i)] >= -ub)
        #     JuMP.@constraint(pm.model, cr[(l,j,i)] <= ub)
        # end
    end

    report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :branch, :cr_fr, :cr_to, _PM.ref(pm, nw, :arcs_from), _PM.ref(pm, nw, :arcs_to), cr)
end

"variable: `ci[l,i,j]` for `(l,i,j)` in `arcs`"
function expression_variable_branch_current_imaginary(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    ci = _PM.var(pm, nw)[:ci] = Dict()

    bus = _PM.ref(pm, nw, :bus)
    branch = _PM.ref(pm, nw, :branch)

    for (l,i,j) in _PM.ref(pm, nw, :arcs_from)
        b = branch[l]
        tm = b["tap"]
        tr, ti = _PM.calc_branch_t(b)
        g_sh_fr, b_sh_fr = b["g_fr"], b["b_fr"]
        g_sh_to, b_sh_to = b["g_to"], b["b_to"]

        vr_fr = _PM.var(pm, nw, :vr, i)
        vi_fr = _PM.var(pm, nw, :vi, i)

        vr_to = _PM.var(pm, nw, :vr, j)
        vi_to = _PM.var(pm, nw, :vi, j)

        csr_fr = _PM.var(pm, nw, :csr, l)
        csi_fr = _PM.var(pm, nw, :csi, l)

        ci[(l,i,j)] = (tr * csi_fr + ti * csr_fr + g_sh_fr * vi_fr + b_sh_fr * vr_fr) / tm^2
        ci[(l,j,i)] = -csi_fr + g_sh_to * vi_to + b_sh_to * vr_to

        # ub = Inf
        # if haskey(b, "rate_a")
        #     rate_fr = b["rate_a"]*b["tap"]
        #     rate_to = b["rate_a"]
        #     ub = max(rate_fr/bus[i]["vmin"], rate_to/bus[j]["vmin"])
        # end
        # if haskey(b, "c_rating_a")
        #     ub = b["c_rating_a"]
        # end

        # if !isinf(ub)
        #     JuMP.@constraint(pm.model, ci[(l,i,j)] >= -ub)
        #     JuMP.@constraint(pm.model, ci[(l,i,j)] <= ub)

        #     JuMP.@constraint(pm.model, ci[(l,j,i)] >= -ub)
        #     JuMP.@constraint(pm.model, ci[(l,j,i)] <= ub)
        # end
    end

    report && _IM.sol_component_value_edge(pm, _PM.pm_it_sym, nw, :branch, :ci_fr, :ci_to, _PM.ref(pm, nw, :arcs_from), _PM.ref(pm, nw, :arcs_to), ci)
end

"variable: `cse[l,i,j]` for `(l,i,j)` in `arcs_from`"
function variable_branch_series_current_expectation(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, aux_fix::Bool=false, report::Bool=true)
    cse = _PM.var(pm, nw)[:cse] = JuMP.@variable(pm.model,
        [l in _PM.ids(pm, nw, :branch)], base_name="$(nw)_cse",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "cse_start", 0.0)
    )

    if bounded
        for l in _PM.ids(pm, nw, :branch)
            JuMP.set_lower_bound(cse[l], 0.0)
        end
    end

    if aux_fix
        JuMP.fix.(cse, 0.0; force=true)
    end

    report && _PM.sol_component_value(pm, nw, :branch, :cse, _PM.ids(pm, nw, :branch), cse)
end

"variable: `csv[l,i,j]` for `(l,i,j)` in `arcs_from`"
function variable_branch_series_current_variance(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, aux_fix::Bool=false, report::Bool=true)
    csv = _PM.var(pm, nw)[:csv] = JuMP.@variable(pm.model,
        [l in _PM.ids(pm, nw, :branch)], base_name="$(nw)_csv",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "csv_start", 0.0)
    )

    if bounded
        for l in _PM.ids(pm, nw, :branch)
                JuMP.set_lower_bound(csv[l], 0.0)
        end
    end

    if aux_fix
        JuMP.fix.(csv, 0.0; force=true)
    end

    report && _PM.sol_component_value(pm, nw, :branch, :csv, _PM.ids(pm, nw, :branch), csv)
end

"variable: `css[l,i,j]` for `(l,i,j)` in `arcs_from`"
function variable_branch_series_current_squared(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, aux_fix::Bool=false, report::Bool=true)
    css = _PM.var(pm, nw)[:css] = JuMP.@variable(pm.model,
        [l in _PM.ids(pm, nw, :branch)], base_name="$(nw)_css",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "css_start", 0.0)
    )

    if bounded
        bus = _PM.ref(pm, nw, :bus)
        branch = _PM.ref(pm, nw, :branch)

        for (l,i,j) in _PM.ref(pm, nw, :arcs_from)
            b = branch[l]
            ub = Inf
            if haskey(b, "rate_a")
                rate = b["rate_a"] * b["tap"]
                y_fr = abs(b["g_fr"] + im * b["b_fr"])
                y_to = abs(b["g_to"] + im * b["b_to"])
                shunt_current = max(y_fr * bus[i]["vmax"]^2, y_to * bus[j]["vmax"]^2)
                series_current = max(rate / bus[i]["vmin"], rate / bus[j]["vmin"])
                ub = series_current + shunt_current
            end
            if haskey(b, "c_rating_a")
                total_current = b["c_rating_a"]
                y_fr = abs(b["g_fr"] + im * b["b_fr"])
                y_to = abs(b["g_to"] + im * b["b_to"])
                shunt_current = max(y_fr * bus[i]["vmax"]^2, y_to * bus[j]["vmax"]^2)
                ub = total_current + shunt_current
            end

            if !isinf(ub)
                JuMP.set_lower_bound(css[l], -2.0 * ub^2)
                JuMP.set_upper_bound(css[l],  2.0 * ub^2)
            end
        end
    end

    if aux_fix
        JuMP.fix.(css, 0.0; force=true)
    end

    report && _PM.sol_component_value(pm, nw, :branch, :css, _PM.ids(pm, nw, :branch), css)
end

# branch voltage drop
"variable: `vbdr[l,i,j]` for `(l,i,j)` in `arcs_from`"
function variable_branch_voltage_drop_real(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    vbdr = _PM.var(pm, nw)[:vbdr] = JuMP.@variable(pm.model,
        [l in _PM.ids(pm, nw, :branch)], base_name="$(nw)_vbdr",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "vbdr_start", 0.01)
    )

    report && _PM.sol_component_value(pm, nw, :branch, :vbdr, _PM.ids(pm, nw, :branch), vbdr)
end

"variable: `vbdi[l,i,j]` for `(l,i,j)` in `arcs_from`"
function variable_branch_voltage_drop_img(pm::AbstractPowerModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)

    vbdi = _PM.var(pm, nw)[:vbdi] = JuMP.@variable(pm.model,
        [l in _PM.ids(pm, nw, :branch)], base_name="$(nw)_vbdi",
        start = _PM.comp_start_value(_PM.ref(pm, nw, :branch, l), "vbdi_start", 0.01)
    )

    report && _PM.sol_component_value(pm, nw, :branch, :vbdi, _PM.ids(pm, nw, :branch), vbdi)
end

"expression: voltage drop real"
function expression_branch_voltage_drop_real(pm::AbstractPowerModel, b::Int; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    if !haskey(_PM.var(pm, nw), :vbdr)
         vbdr=_PM.var(pm, nw)[:vbdr] = Dict()
    end

        branch = _PM.ref(pm, nw, :branch, b)
        f_bus = branch["f_bus"]
        t_bus = branch["t_bus"]

        #vbdr[b] = (vr_fr-vr_to)

        expression_branch_voltage_drop_real(pm, nw, b, f_bus, t_bus)
end

"expression: voltage drop img"
function expression_branch_voltage_drop_img(pm::AbstractPowerModel, b::Int; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true)
    if !haskey(_PM.var(pm, nw), :vbdi)
        vbdi=_PM.var(pm, nw)[:vbdi] = Dict()
    end
        #branch = _PMs.ref(pm, nw, :branch)
    branch = _PM.ref(pm, nw, :branch, b)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]

    expression_branch_voltage_drop_img(pm, nw, b, f_bus, t_bus)
end


text:
# vary transition prob, mcts_depth, and number of nodes in network
cond1 = make_condition_dict(num_nets=500,mcts_d=[10],n_bnds=[8,35],seed=5,role=:trusted)
cond2 = make_condition_dict(num_nets=100,mcts_d=[1],n_bnds=[8,35],seed=4)
cond3 = make_condition_dict(num_nets=100,mcts_d=[2],n_bnds=[8,35],seed=3)
cond4 = make_condition_dict(num_nets=100,mcts_d=[3],n_bnds=[8,35],seed=2)
cond5 = make_condition_dict(num_nets=100,mcts_d=[4,7],n_bnds=[8,35],seed=1)

exp_dict = Dict(:name=>"mturk", :conditions=>Dict(1=>cond1, 2=>cond2, 3=> cond3,4=> cond4,5=> cond5),
                :xQ=>Dict(:inpts=>[:tprob,:N],:epocs=>1000,:cmp=>[2,3,4,5],:nn_loc=>"nn_logs"))

fldr = "logs"
for x in keys(exp_dict[:conditions])
    fname = make_useful_fname(exp_dict,x)
    exp_dict[:conditions][x][:fname] = fname
    exp_dict[:conditions][x][:fldr] = fldr
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using Test
using Flux: onehotbatch, mse, crossentropy, logitcrossentropy,
            σ, bce_loss, logitbce_loss, flatten,
            xlogx, xlogy

const ϵ = 1e-7

@testset "xlogx & xlogy" begin
  @test iszero(xlogx(0))
  @test isnan(xlogx(NaN))
  @test xlogx(2) ≈ 2.0 * log(2.0)
  @inferred xlogx(2)
  @inferred xlogx(0)

  @test iszero(xlogy(0, 1))
  @test isnan(xlogy(NaN, 1))
  @test isnan(xlogy(1, NaN))
  @test isnan(xlogy(NaN, NaN))
  @test xlogy(2, 3) ≈ 2.0 * log(3.0)
  @inferred xlogy(2, 3)
  @inferred xlogy(0, 1)
end

@testset "losses" begin
  # First, regression-style y's
  y = [1, 1, 0, 0]
  ŷ = [.9, .1, .1, .9]

  @testset "mse" begin
    @test mse(ŷ, y) ≈ (.1^2 + .9^2)/2
  end

  @testset "mae" begin
    @test Flux.mae(ŷ, y) ≈ 1/2
  end

  @testset "huber_loss" begin
    @test Flux.huber_loss(ŷ, y) ≈ 0.20500000000000002
  end

  y = [123.0,456.0,789.0]
  ŷ = [345.0,332.0,789.0]
  @testset "msle" begin
    @test Flux.msle(ŷ, y) ≈ 0.38813985859136585
  end

  # Now onehot y's
  y = onehotbatch([1, 1, 0, 0], 0:1)
  ŷ = [.1 .9; .9 .1; .9 .1; .1 .9]'
  v = log(.1 / .9)
  logŷ = [v 0.0; 0.0 v; 0.0 v; v 0.0]'
  lossvalue = 1.203972804325936

  @testset "crossentropy" begin
    @test crossentropy([0.1,0.0,0.9], [0.1,0.0,0.9]) ≈ crossentropy([0.1,0.9], [0.1,0.9])
    @test crossentropy(ŷ, y) ≈ lossvalue
  end

  @testset "logitcrossentropy" begin
    @test logitcrossentropy(logŷ, y) ≈ lossvalue
  end

  logŷ, y = randn(3), rand(3)
  @testset "bce_loss" begin
    @test bce_loss(σ.(logŷ), y; ϵ=0) ≈ mean(-y.*log.(σ.(logŷ)) - (1 .- y).*log.(1 .- σ.(logŷ)))
    @test bce_loss(σ.(logŷ), y) ≈ mean(-y.*log.(σ.(logŷ) .+ eps.(σ.(logŷ))) - (1 .- y).*log.(1 .- σ.(logŷ) .+ eps.(σ.(logŷ))))
  end

  @testset "logitbce_loss" begin
    @test logitbce_loss(logŷ, y) ≈ bce_loss(σ.(logŷ), y; ϵ=0)
  end

  y = [1 2 3]
  ŷ = [4.0 5.0 6.0]
  @testset "kldivergence" begin
    @test Flux.kldivergence([0.1,0.0,0.9], [0.1,0.0,0.9]) ≈ Flux.kldivergence([0.1,0.9], [0.1,0.9])
    @test Flux.kldivergence(ŷ, y) ≈ -1.7661057888493457
    @test Flux.kldivergence(y, y) ≈ 0
  end

  y = [1 2 3 4]
  ŷ = [5.0 6.0 7.0 8.0]
  @testset "hinge_loss" begin
    @test Flux.hinge_loss(ŷ, y) ≈ 0
    @test Flux.hinge_loss(y, 0.5 .* y) ≈ 0.125
  end

  @testset "squared_hinge_loss" begin
    @test Flux.squared_hinge_loss(ŷ, y) ≈ 0
    @test Flux.squared_hinge_loss(y, 0.5 .* y) ≈ 0.0625
  end

  y = [0.1 0.2 0.3]
  ŷ = [0.4 0.5 0.6]
  @testset "poisson_loss" begin
    @test Flux.poisson_loss(ŷ, y) ≈ 0.6278353988097339
    @test Flux.poisson_loss(y, y) ≈ 0.5044459776946685
  end

  y = [1.0 0.5 0.3 2.4]
  ŷ = [0 1.4 0.5 1.2]
  @testset "dice_coeff_loss" begin
    @test Flux.dice_coeff_loss(ŷ, y) ≈ 0.2799999999999999
    @test Flux.dice_coeff_loss(y, y) ≈ 0.0
  end

  @testset "tversky_loss" begin
    @test Flux.tversky_loss(ŷ, y) ≈ -0.06772009029345383
    @test Flux.tversky_loss(ŷ, y, β=0.8) ≈ -0.09490740740740744
    @test Flux.tversky_loss(y, y) ≈ -0.5576923076923075
  end

  @testset "no spurious promotions" begin
    for T in (Float32, Float64)
      y = rand(T, 2)
      ŷ = rand(T, 2)
      for f in (mse, crossentropy, logitcrossentropy, Flux.kldivergence, Flux.hinge_loss, Flux.poisson_loss,
              Flux.mae, Flux.huber_loss, Flux.msle, Flux.squared_hinge_loss, Flux.dice_coeff_loss, Flux.tversky_loss)
        fwd, back = Flux.pullback(f, ŷ, y)
        @test fwd isa T
        @test eltype(back(one(T))[1]) == T
      end
    end
  end
end

@testset "helpers" begin
  @testset "flatten" begin
    x = randn(Float32, 10, 10, 3, 2)
    @test size(flatten(x)) == (300, 2)
  end
end


text:
##############################################################################
#
# InvertedGamma distribution, cf Bayesian Theory (Barnardo & Smith) p 119
#
#  Note that B&S parametrize in terms of shape/rate, but this is uses
#  shape/scale so that it is consistent with the implementation of Gamma
#
##############################################################################

immutable InvertedGamma <: ContinuousUnivariateDistribution
    shape::Float64 # location
    scale::Float64 # scale
    function InvertedGamma(sh::Real, sc::Real)
    	if sh > 0.0 && sc > 0.0
    		new(float64(sh), float64(sc))
    	else
    		error("Both shape and scale must be positive")
    	end
    end
end

cdf(d::InvertedGamma, x::Real) = 1.0 - cdf(Gamma(d.shape, d.scale), 1.0 / x)

function entropy(d::InvertedGamma)
    a = d.shape
    b = d.scale
    return a + log(b) + lgamma(a) - (1.0 + a) * digamma(a)
end

insupport(d::InvertedGamma, x::Number) = isreal(x) && isfinite(x) && 0 <= x

function kurtosis(d::InvertedGamma)
    a = d.shape
    if a > 4.0
        return (30.0 * a - 66.0) / ((a - 3.0) * (a - 4.0))
    else
        error("kurtosis not defined unless shape > 4")
    end
end

function mean(d::InvertedGamma)
	if d.shape > 1.0
		return (1.0 / d.scale) / (d.shape - 1.0)
	else
		error("Expectation only defined if shape > 1")
	end
end

function mgf(d::InvertedGamma, t::Real)
    a, b = d.shape, d.scale
    return (2 * (-b * t)^(a / 2)) / gamma(a) *
           besselk(a, sqrt(-4.0 * b * t))
end

function cf(d::InvertedGamma, t::Real)
    a, b = d.shape, d.scale
    return (2 * (-im * b * t)^(a / 2)) / gamma(a) *
           besselk(a, sqrt(-4.0 * im * b * t))
end

modes(d::InvertedGamma) = [(1.0 / d.scale) / (d.shape + 1.0)]

pdf(d::InvertedGamma, x::Real) = exp(logpdf(d, x))

function logpdf(d::InvertedGamma, x::Real)
    return -(d.shape * log(d.scale)) - lgamma(d.shape) -
           ((d.shape + 1.0) * log(x)) - 1.0/ (d.scale * x)
end

function quantile(d::InvertedGamma, p::Real)
    return 1.0 / quantile(Gamma(d.shape, d.scale), 1.0 - p)
end

rand(d::InvertedGamma) = 1.0 / rand(Gamma(d.shape, d.scale))

function rand!(d::InvertedGamma, A::Array{Float64})
    A = rand!(Gamma(d.shape, d.scale), A)
    for i in 1:length(A)
    	A[i] = 1.0 / A[i]
    end
    return A
end

function skewness(d::InvertedGamma)
    a = d.shape
    if a > 3.0
        return (4.0 * sqrt(a - 2.0)) / (a - 3.0)
    else
        error("skewness not defined unless shape > 3")
    end
end

function var(d::InvertedGamma)
	if d.shape > 2.0
		(1.0 / d.scale)^2 / ((d.shape - 1.0)^2 * (d.shape - 2.0))
	else
		error("Variance only defined if shape > 2")
	end
end


text:
"""
$(TYPEDSIGNATURES)

```jldoctest
julia> random_design((Uniform(2, 3), DiscreteUniform(-1, 5), Uniform(5, 10)), 10)
10×3 Array{Real,2}:
 2.04922   5  8.85741
 2.25659  -1  6.57617
 2.86526   4  5.41422
 2.81201  -1  5.40944
 2.92412   1  6.22902
 2.20051   5  8.80749
 2.69459   3  6.34626
 2.28902   2  8.72759
 2.95173   0  5.42517
 2.35163   4  5.89413

```
"""
function random_design(distributions::Tuple, size::Int)
    permutedims(foldl((x, y) -> hcat(collect(x),
                                     collect(y)),
                      (rand(d) for d in distributions) for i = 1:size))
end


text:
export triangulate

"""
    triangulate(poly::Polygon{N,Point{2,T}}) where {N,T}

Return a Vec of the `N`-2 triangles that partition the 2D `Polygon`.

Generated using fan triangulation, which assumes the `Polygon` is convex.
"""
function triangulate(p::Polygon{N, Point{2,T}}) where {N,T}
    return Vec{N-2, Triangle{Point{2,T}}}(
                ntuple(i->Triangle(p[1], p[i+1], p[i+2]), Val(N-2))
               )
end

"""
    triangulate(quad::Quadrilateral{Point{3,T}}, ::Val{N}) where {N,T}

Return a Vec of the (2`N`+2)^2 triangles that approximately partition the `Quadrilateral`.
"""
@generated function triangulate(quad::Quadrilateral{Point{3,T}}, ::Val{N}) where {T,N}
    # N is the number of divisions of each edge
    N1 = N + 1
    inv_N1 = T(1)/N1
    if N === 0
        exprs = [:(Triangle(quad[1], quad[2], quad[3])),
                 :(Triangle(quad[3], quad[4], quad[1]))
                ]
    else
        exprs = Expr[]
        for j ∈ 0:N
            s₀ = inv_N1*j
            s₁ = inv_N1*(j + 1)
            for i ∈ 0:N
                r₀ = inv_N1*i
                r₁ = inv_N1*(i + 1)
                push!(exprs, :(Triangle(quad($r₀, $s₀), quad($r₁, $s₀), quad($r₀, $s₁))))
                push!(exprs, :(Triangle(quad($r₀, $s₁), quad($r₁, $s₀), quad($r₁, $s₁))))
            end
        end
    end
    return quote
        return Vec{$(2N1^2), Triangle{Point{3,$T}}}(tuple($(exprs...)))
    end
end

"""
    triangulate(tri6::QuadraticTriangle{Point{Dim,T}}, ::Val{N}) where {Dim,T,N}

Return a Vec of (2`N`+2)^2 triangles that approximately partition the `QuadraticTriangle`.
"""
@generated function triangulate(tri6::QuadraticTriangle{Point{Dim,T}}, ::Val{N}) where {Dim,T,N}
    # N is the number of divisions of each edge
    N1 = N + 1
    inv_N1 = T(1)/N1
    if N === 0
        exprs = [:(Triangle(tri6[1], tri6[2], tri6[3]))]
    else
        exprs = Expr[]
        for s ∈ 1:N
            s₋₁ = inv_N1*(s - 1)
            s₀  = inv_N1*s
            s₁  = inv_N1*(s + 1)
            for r ∈ 0:N-s
                r₀ = inv_N1*r
                r₁ = inv_N1*(r + 1)
                push!(exprs, :(Triangle(tri6($r₀, $s₀), tri6($r₁, $s₀ ), tri6($r₀, $s₁))))
                push!(exprs, :(Triangle(tri6($r₀, $s₀), tri6($r₁, $s₋₁), tri6($r₁, $s₀))))
            end
        end
        j = N1*N + 1
        s₀ = zero(T)
        s₁ = inv_N1
        for r ∈ 0:N
            r₀ = inv_N1*r
            r₁ = inv_N1*(r + 1)
            push!(exprs, :(Triangle(tri6($r₀, $s₀), tri6($r₁, $s₀), tri6($r₀, $s₁))))
        end
    end
    return quote
        return Vec{$(N1^2), Triangle{Point{$Dim,$T}}}(tuple($(exprs...)))
    end
end

"""
    triangulate(quad::QuadraticQuadrilateral{Point{Dim,T}}, ::Val{N}) where {Dim,T,N}

Return a Vec of (2`N`+2)^2 triangles that approximately partition the `QuadraticQuadrilateral`.
"""
@generated function triangulate(quad8::QuadraticQuadrilateral{Point{Dim,T}},
                                ::Val{N}) where {Dim,T,N}
    # N is the number of divisions of each edge
    N1 = N + 1
    inv_N1 = T(1)/N1
    if N === 0
        exprs = [:(Triangle(quad8[1], quad8[2], quad8[3])),
                 :(Triangle(quad8[3], quad8[4], quad8[1]))
                ]
    else
        exprs = Expr[]
        for j ∈ 0:N
            s₀ = inv_N1*j
            s₁ = inv_N1*(j + 1)
            for i ∈ 0:N
                r₀ = inv_N1*i
                r₁ = inv_N1*(i + 1)
                push!(exprs, :(Triangle(quad8($r₀, $s₀), quad8($r₁, $s₀), quad8($r₀, $s₁))))
                push!(exprs, :(Triangle(quad8($r₀, $s₁), quad8($r₁, $s₀), quad8($r₁, $s₁))))
            end
        end
    end
    return quote
        return Vec{$(2N1^2), Triangle{Point{3,$T}}}(tuple($(exprs...)))
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using Gadfly

n = 10
plot(x=rand(n), y=rand(n),
     Guide.xticks(ticks=[0.0, 0.1, 0.9, 1.0]),
     Guide.yticks(ticks=[0.4, 0.5, 0.6]))


text:
module BART

  using Distributed
  using Distributions
  using LinearAlgebra
  using MCMCChains
  using StatsBase
  using SpecialFunctions

  import DecisionTree

  include("trees.jl")
  include("models.jl")
  include("proposals.jl")
  include("predict.jl")
  include("fit.jl")

  export
    BartModel,
    Chains,
    Hypers,
    Opts,
    TrainData

  export
    depth,
    fit,
    ptd,
    predict,
    update

end


text:
module JSONRPCEndpoints

import JSON, UUIDs

export JSONRPCEndpoint, send_notification, send_request, send_success_response, send_error_response

mutable struct JSONRPCEndpoint
    pipe_in
    pipe_out

    out_msg_queue::Channel{Any}
    in_msg_queue::Channel{Any}

    outstanding_requests::Dict{String, Channel{Any}}

    err_handler::Union{Nothing,Function}

    function JSONRPCEndpoint(pipe_in, pipe_out, err_handler=nothing)
        return new(pipe_in, pipe_out, Channel{Any}(Inf), Channel{Any}(Inf), Dict{String, Channel{Any}}(), err_handler)
    end
end

function write_transport_layer(stream, response)
    response_utf8 = transcode(UInt8, response)
    n = length(response_utf8)
    write(stream, "Content-Length: $n\r\n\r\n")
    write(stream, response_utf8)
end

function read_transport_layer(stream)
    header_dict = Dict{String,String}()
    line = chomp(readline(stream))
    # Check whether the socket was closed
    if line == ""
        return nothing
    end
    while length(line) > 0
        h_parts = split(line, ":")
        header_dict[chomp(h_parts[1])] = chomp(h_parts[2])
        line = chomp(readline(stream))
    end
    message_length = parse(Int, header_dict["Content-Length"])
    message_str = String(read(stream, message_length))
    return message_str
end

function Base.run(x::JSONRPCEndpoint)
    @async try
        for msg in x.out_msg_queue
            write_transport_layer(x.pipe_out, msg)
        end
    catch err
        bt = catch_backtrace()
        if x.err_handler!==nothing
            x.err_handler(err, bt)
        else
            Base.display_error(stderr, err, bt)
        end
    end

    @async try
        while true
            message = read_transport_layer(x.pipe_in)

            if message===nothing
                break
            end

            message_dict = JSON.parse(message)

            if haskey(message_dict, "method")
                put!(x.in_msg_queue, message_dict)
            else
                # This must be a response
                id_of_request = message_dict["id"]

                channel_for_response = x.outstanding_requests[id_of_request]
                put!(channel_for_response, message_dict)
            end
        end
    catch err
        bt = catch_backtrace()
        if x.err_handler!==nothing
            x.err_handler(err, bt)
        else
            Base.display_error(stderr, err, bt)
        end
    end
end

function send_notification(x::JSONRPCEndpoint, method::AbstractString, params)
    message = Dict("jsonrpc"=>"2.0", "method"=>method, "params"=>params)

    message_json = JSON.json(message)

    put!(x.out_msg_queue, message_json)
end

function send_request(x::JSONRPCEndpoint, method::AbstractString, params)
    id = string(UUIDs.uuid4())
    message = Dict("jsonrpc"=>"2.0", "method"=>method, "params"=>params, "id"=>id)

    response_channel = Channel{Any}(1)
    x.outstanding_requests[id] = response_channel

    message_json = JSON.json(message)

    put!(x.out_msg_queue, message_json)

    response = take!(response_channel)

    if haskey(response, "result")
        return response["result"]
    else
        @info "We recieved an error response: $response"
        # error("Throw a better error here.")
    end
end

function get_next_message(endpoint::JSONRPCEndpoint)
    msg = take!(endpoint.in_msg_queue)

    return msg
end

function send_success_response(endpoint, original_request, result)
    response = Dict("jsonrpc"=>"2.0", "id"=>original_request["id"], "result"=>result)

    response_json = JSON.json(response)

    put!(endpoint.out_msg_queue, response_json)
end

function send_error_response(endpoint, original_request, code, message, data)
    response = Dict("jsonrpc"=>"2.0", "id"=>original_request["id"], "error"=>Dict("code"=>code, "message"=>message, "data"=>data))

    response_json = JSON.json(response)

    put!(endpoint.out_msg_queue, response_json)
end

end


text:
module EasyTranspose

export ᵀ

struct ᵀ end
Base.:(*)(arr::AbstractVecOrMat,::typeof(ᵀ)) = permutedims(arr)

end # module


text:
using Base.Test
using MathProgBase
using CPLEX

include("constants.jl")
include("low_level_api.jl")
@testset "LP" begin
    include("lp_01.jl")
    include("lp_02.jl")
    include("lp_03.jl")
end
include("mip_01.jl")
@testset "QP" begin
    include("qp_01.jl")
    include("qp_02.jl")
end
include("qcqp_01.jl")
include("env.jl")
include("sos.jl")
include("problemtype.jl")
include("miqcp.jl")
@testset "MathProgBase" begin
    include("mathprog.jl")
end

include("MOIWrapper.jl")