text:
"""
`module  JAC.PhotoExcitation`
    ... a submodel of JAC that contains all methods for computing photo-excitation properties between some initial
        and final-state multiplets.
"""
module PhotoExcitation

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

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

        + multipoles              ::Array{EmMultipole,1}         ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge,1}            ... Specifies the gauges to be included into the computations.
        + calcForStokes           ::Bool                         ... True, if the excitation cross sections are to be calculated (and false otherwise)
                                                                     for given Stokes parameter of the incident plane-wave photons.
        + calcPhotonDm            ::Bool                         ... True, if the photon density matrix of a subsequently emitted fluorescence photon is to be
                                                                     calculated and false otherwise.
        + calcTensors             ::Bool                         ... True, if the statistical tensors of the excited atom are to be calculated and false otherwise.
        + printBeforeComputation  ::Bool                         ... True, if all energies and lines are printed before their evaluation.
        + selectLines             ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines           ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
        + 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.
        + stokes                  ::ExpStokes                    ... Stokes parameters of the incident radiation.
    """
    struct Settings
        multipoles                ::Array{EmMultipole,1}
        gauges                    ::Array{UseGauge,1}
        calcForStokes             ::Bool
        calcPhotonDm              ::Bool
        calcTensors               ::Bool
        printBeforeComputation    ::Bool
        selectLines               ::Bool
        selectedLines             ::Array{Tuple{Int64,Int64},1}
        photonEnergyShift         ::Float64
        mimimumPhotonEnergy       ::Float64
        maximumPhotonEnergy       ::Float64
        stokes                    ::ExpStokes
    end

    """
    `PhotoExcitation.Settings()`  ... 'empty' constructor for the default values of photo-excitation line computations
    """
    function Settings()
        Settings(EmMultipole[], UseGauge[], false, false, false, false, false, Array{Tuple{Int64,Int64},1}[], 0., 0., 0.)
    end

    # `Base.show(io::IO, settings::PhotoExcitation.Settings)`  ... prepares a proper printout of the variable settings::PhotoExcitation.Settings.
    function Base.show(io::IO, settings::PhotoExcitation.Settings)
        println(io, "multipoles:               $(settings.multipoles)  ")
        println(io, "use-gauges:               $(settings.gauges)  ")
        println(io, "calcForStokes:            $(settings.calcForStokes)  ")
        println(io, "calcPhotonDm:             $(settings.calcPhotonDm)  ")
        println(io, "calcTensor:               $(settings.calcTensor)  ")
        println(io, "printBeforeComputation:   $(settings.printBeforeComputation)  ")
        println(io, "selectLines:              $(settings.selectLines)  ")
        println(io, "selectedLines:            $(settings.selectedLines)  ")
        println(io, "photonEnergyShift:        $(settings.photonEnergyShift)  ")
        println(io, "mimimumPhotonEnergy:      $(settings.mimimumPhotonEnergy)  ")
        println(io, "maximumPhotonEnergy:      $(settings.maximumPhotonEnergy)  ")
        println(io, "stokes:                   $(settings.stokes)  ")
    end

    """
    `struct  PhotoExcitation.Line`
        ... defines a type for a photo-excitation 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.
        + crossSection   ::EmProperty                  ... Total cross section of this line.
        + staTensor      ::Array{TensorComp,1}         ... Array of statistical tensor components rho_kq
        + hasSublines    ::Bool                        ... Determines whether the individual sublines are defined in terms of their
                                                           multipolarity, amplitude, or not; cf. PhotoEmission.Channel
        + channels       ::Array{PhotoEmission.Channel,1}  ... List of radiative (photon) channels
    """
    struct  Line
        initialLevel     ::Level
        finalLevel       ::Level
        omega            ::Float64
        crossSection     ::EmProperty
        staTensor        ::Array{TensorComp,1}
        hasSublines      ::Bool
        channels         ::Array{PhotoEmission.Channel,1}
    end

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

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

    """
    `PhotoExcitation.computeAmplitudesProperties(line::PhotoExcitation.Line, grid::Radial.Grid, settings::PhotoExcitation.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::PhotoExcitation.Line is returned for which
            the amplitudes and properties have now been evaluated.
    """
    function  computeAmplitudesProperties(line::PhotoExcitation.Line, grid::Radial.Grid, settings::PhotoExcitation.Settings)
        newChannels = PhotoEmission.Channel[]
        for  channel  in  line.channels
            ##x println("channel = $channel ")
            amplitude = PhotoEmission.amplitude("absorption", channel.multipole, channel.gauge, line.omega, line.finalLevel, line.initialLevel, grid)
            push!( newChannels, PhotoEmission.Channel( channel.multipole, channel.gauge, amplitude) )
        end
        # Calculate the crossSection and stastistical tensors if requested
        crossSection = EmProperty(-1., -1.)    # Cross sections for linear-polarized plane-wave photons
        staTensor    = TensorComp[]
        if  settings.calcTensors
            push!(staTensor, TensorComp(0, 0, 1.))
        end
        line = PhotoExcitation.Line( line.initialLevel, line.finalLevel, line.omega, crossSection, staTensor, true, newChannels)
        return( line )
    end

    """
    `PhotoExcitation.computeCrossSection(line::PhotoExcitation.Line, stokes::ExpStokes)`
        ... to compute the excitation cross section for the excitation of unpolarized atoms by plane-wave photons, whose polarization
            is described by the given (experimental) Stokes parameters. A cs::EmProperty is returned.
    """
    function  computeCrossSection(line::PhotoExcitation.Line, stokes::ExpStokes)
        csCoulomb = -2.;    csBabushkin = -2.
        cs = EmProperty( csCoulomb, csBabushkin )
        return( cs )
    end

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

    """
    `PhotoExcitation.computeMatrix_obsolete(Mp::EmMultipole, gauge::EmGauge, omega::Float64, finalBasis::Basis, initialBasis::Basis,
                                            grid::Radial.Grid, settings::PhotoExcitation.Settings)`
        ... to compute the transition matrix (O^Mp_rs) = (<finalCSF_r|| O^Mp (omega; gauge) ||initialCSF_s>) of the Mp multiplet field
            for the given transition energy and gauge, and between the CSF_r from the finalBasis and the CSF_s from the initialBasis.
            A (non-quadratic) matrix::Array{Float64,2} with dimensions [length(finalBasis.csfs) x length(initialBasis.csfs)] is returned.
            Note that this transition matrix is specific to transitions of the given multipolarity and omega.
    """
    function computeMatrix_obsolete(Mp::EmMultipole, gauge::EmGauge, omega::Float64, finalBasis::Basis, initialBasis::Basis,
                                    grid::Radial.Grid, settings::PhotoExcitation.Settings)
        nf = length(finalBasis.csfs);    ni = length(initialBasis.csfs)

        print("Compute radiative T^L matrix of dimension $nf x $ni for the given initial- and final-state bases ...")
        matrix = zeros(Float64, nf, ni)
        for  r = 1:nf
            for  s = 1:ni
                wa = compute("angular coefficients: 1-p, Ratip2013", Mp.L, finalBasis.csfs[r], initialBasis.csfs[s])
                me = 0.
                for  coeff in wa
                    jj = Basics.subshell_2j(finalBasis.orbitals[coeff.a].subshell)
                    me = me + coeff.T * sqrt( jj + 1) * InteractionStrength.multipole(Mp, gauge, omega, finalBasis.orbitals[coeff.a],
                                                                                initialBasis.orbitals[coeff.b], grid)
                end
                matrix[r,s] = me
            end
        end
        println("   ... done.")
        return( matrix )
    end

    """
    `PhotoExcitation.computeStatisticalTensor(k::Int64, q::Int64, line::PhotoExcitation.Line, stokes::ExpStokes)`
        ... to compute the statistical tensor (component) rho_{k,q} of the final level for the excitation of unpolarized atoms by
            plane-wave photons, whose polarization is described by the given (experimental) Stokes parameters.
            A rho_kq::EmProperty is returned.
    """
    function  computeStatisticalTensor(k::Int64, q::Int64, line::PhotoExcitation.Line, stokes::ExpStokes)
        if  k == q  ||  k == -q   rhoCoulomb = -3.;    rhoBabushkin = -3.    else   rhoCoulomb = 0.;    rhoBabushkin = 0.    end
        rho = EmProperty( rhoCoulomb, rhoBabushkin )
        return( rho )
    end

    """
    `PhotoExcitation.determineChannels(finalLevel::Level, initialLevel::Level, settings::PhotoExcitation.Settings)`
        ... to determine a list of PhotoExcitation.Channel for a transitions from the initial to final level and by taking
            into account the particular settings of for this computation; an Array{PhotoExcitation.Channel,1} is returned.
    """
    function determineChannels(finalLevel::Level, initialLevel::Level, settings::PhotoExcitation.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
        ##x println("PhotoExcitation.determineChannels-aa: channels = $channels ")
        return( channels )
    end

    """
    `PhotoExcitation.determineLines(finalMultiplet::Multiplet, initialMultiplet::Multiplet, settings::PhotoExcitation.Settings)`
        ... to determine a list of photo-excitation 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{PhotoExcitation.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::PhotoExcitation.Settings)
        if    settings.selectLines    selectLines   = true;   selectedLines = Basics.determine("selected lines", settings.selectedLines)
        else                          selectLines   = false
        end

        lines = PhotoExcitation.Line[]
        for  i = 1:length(initialMultiplet.levels)
            for  f = 1:length(finalMultiplet.levels)
                if  selectLines  &&  !(haskey(selectedLines, (i,f)) )    continue   end
                ##x println("PhotoExcitation.determineLines-aa: angular i = $i, f = $f")
                omega    = abs( finalMultiplet.levels[f].energy - initialMultiplet.levels[i].energy) + settings.photonEnergyShift
                if   omega == 0.  ||  omega < settings.mimimumPhotonEnergy  ||  omega > settings.maximumPhotonEnergy    continue   end

                channels = PhotoExcitation.determineChannels(finalMultiplet.levels[f], initialMultiplet.levels[i], settings)
                if   length(channels) == 0   continue   end
                push!( lines, PhotoExcitation.Line(initialMultiplet.levels[i], finalMultiplet.levels[f], omega, EmProperty(0., 0.),
                                                   TensorComp[], true, channels) )
            end
        end
        return( lines )
    end

    """
    `PhotoExcitation.displayLines(lines::Array{PhotoExcitation.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{PhotoExcitation.Line,1})
        println(" ")
        println("  Selected photo-excitation lines:")
        println(" ")
        println("  ", TableStrings.hLine(120))
        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(120))
        #
        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{Basics.EmMultipole,Basics.EmGauge}[]
            for  i in 1:length(line.channels)
                push!( mpGaugeList, (line.channels[i].multipole, line.channels[i].gauge) )
            end
            ##x println("PhotoExcitation.diplayLines-ad: mpGaugeList = ", mpGaugeList)
            sa = sa * TableStrings.multipoleGaugeTupels(50, mpGaugeList)
            println( sa )
        end
        println("  ", TableStrings.hLine(120))
        #
        return( nothing )
    end

    """
    `PhotoExcitation.displayCrossSections(stream::IO, lines::Array{PhotoExcitation.Line,1}, settings::PhotoExcitation.Settings)`
        ... to list all results, energies, cross sections, etc. of the selected lines. A neat table is printed but nothing
            is returned otherwise.
    """
    function  displayCrossSections(stream::IO, lines::Array{PhotoExcitation.Line,1}, settings::PhotoExcitation.Settings)
        println(stream, " ")
        println(stream, "  Photo-excitation cross sections for (completely) linearly-polarized plane-wave photons:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(105))
        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=3);
        sb = sb * TableStrings.center(14,TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center(10, "Multipoles"; na=3);                        sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(30, "Cou -- Cross section -- Bab"; na=2);
        sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section")*"          "*
                                              TableStrings.inUnits("cross section"); na=2)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(105))
        #
        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)) * "    "
            multipoles = EmMultipole[]
            for  ch in line.channels
                multipoles = push!( multipoles, ch.multipole)
            end
            multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
            sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Coulomb))     * "    "
            sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", line.crossSection.Babushkin))   * "    "
            ##x sa = sa * @sprintf("%.6e", line.anisotropy.Coulomb)      * "    "
            ##x sa = sa * @sprintf("%.6e", line.anisotropy.Babushkin)    * "    "
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(105))
        #
        #
        if  settings.calcForStokes
            stokes = settings.stokes
            println(stream, " ")
            println(stream, "  Photo-excitation cross sections for incident plane-wave photons with given $stokes:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(105))
            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=3);
            sb = sb * TableStrings.center(14,TableStrings.inUnits("energy"); na=3)
            sa = sa * TableStrings.center(10, "Multipoles"; na=3);                        sb = sb * TableStrings.hBlank(14)
            sa = sa * TableStrings.center(30, "Cou -- Cross section -- Bab"; na=2);
            sb = sb * TableStrings.center(30, TableStrings.inUnits("cross section")*"          "*
                                                  TableStrings.inUnits("cross section"); na=2)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(105))
            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)) * "    "
                multipoles = EmMultipole[]
                for  ch in line.channels
                    multipoles = push!( multipoles, ch.multipole)
                end
                multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
                sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
                crossSection = PhotoExcitation.computeCrossSection(line, stokes)
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", crossSection.Coulomb))     * "    "
                sa = sa * @sprintf("%.6e", Defaults.convertUnits("cross section: from atomic", crossSection.Babushkin))   * "    "
                println(stream, sa)
            end
            println(stream, "  ", TableStrings.hLine(105))
        end
        #
        #
        if  settings.calcTensors
            stokes = settings.stokes
            println(stream, " ")
            println(stream, "  Statistical tensors rho_kq  and alignment parameters A_kq for the excitation by incident plane-wave photons")
            println(stream, "  with given $stokes:")
            println(stream, " ")
            println(stream, "  ", TableStrings.hLine(145))
            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=3);
            sb = sb * TableStrings.center(14,TableStrings.inUnits("energy"); na=3)
            sa = sa * TableStrings.center(10, "Multipoles"; na=1);                        sb = sb * TableStrings.hBlank(12)
            sa = sa * TableStrings.center(12, "     k    q    "; na=0);                   sb = sb * TableStrings.hBlank(14)
            sa = sa * TableStrings.center(26, "Cou -- rho_kq -- Bab"; na=5);              sb = sb * TableStrings.hBlank(14)
            sa = sa * TableStrings.center(26, "Cou --  A_kq  -- Bab"; na=2);              sb = sb * TableStrings.hBlank(14)
            println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(145))
            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)) * "    "
                multipoles = EmMultipole[]
                for  ch in line.channels
                    multipoles = push!( multipoles, ch.multipole)
                end
                multipoles = unique(multipoles);   mpString = TableStrings.multipoleList(multipoles) * "          "
                sa = sa * TableStrings.flushleft(11, mpString[1:10];  na=3)
                println(stream, sa);   sc = TableStrings.hBlank( length(sa)-1 )
                tc00   = PhotoExcitation.computeStatisticalTensor(0, 0, line, stokes)
                for  k = 2:2
                    for  q = -k:k
                        tc = PhotoExcitation.computeStatisticalTensor(k, q, line, stokes)
                        if  abs(tc.Coulomb) == abs(tc.Babushkin) == 0.    continue    end
                        #
                        aCoulomb = tc.Coulomb / tc00.Coulomb;    aBabushkin = tc.Babushkin / tc00.Babushkin
                        sd = sc * TableStrings.level(k) * TableStrings.level(q) * "   "
                        sd = sd * @sprintf("%.6e", tc.Coulomb)     * "  "
                        sd = sd * @sprintf("%.6e", tc.Babushkin)   * "    "
                        sd = sd * @sprintf("%.6e", aCoulomb)       * "  "
                        sd = sd * @sprintf("%.6e", aBabushkin)     * "    "
                        println(stream, sd)
                    end
                end
            end
            println(stream, "  ", TableStrings.hLine(145))
        end
        #
        #
        if  settings.calcTensors   &&   settings.calcPhotonDm
            stokes = settings.stokes
            println(stream, " ")
            println(stream, "  Reduced density matrix of the fluorescence photons at selected solid angles for the prior excitation by incident plane-wave photons")
            println(stream, "  with given $stokes:")
            println(stream, " ")
            println(stream, "  Not yet implemented !!")
        end
        #
        return( nothing )
    end

end # module


text:
"""
ODBC API Function Definitions
By Jacob Quinn, 2016
In general, the ODBC functions are implemented to mirror the C header files (sql.h,sqlext.h,sqltypes.h,sqlucode.h)
A few liberties are taken in utliizing standard Julia functions and idioms
Format:
  * function name
  * URL reference
  * short function description
  * valid const definitions
  * relevant notes
  * working and tested status
  * function definition code

Contents
 * Macros and Utility Functions
 * Handle Functions
 * Connection Functions
 * Resultset Metadata Functions
 * Query Functions
 * Resultset Retrieval Functions
 * DBMS Meta Functions
 * Error Handling and Diagnostics
"""
module API

using WeakRefStrings
import Compat: @__MODULE__, Sys

if VERSION < v"0.7.0-DEV.2575"
    const Dates = Base.Dates
else
    import Dates
end

include("types.jl")

#### Macros and Utility Functions ####

"""
# MAXFETCHSIZE sets the default rowset fetch size
# used in retrieving resultset blocks from queries
"""
const MAXFETCHSIZE = 65535

# success codes
const SQL_SUCCESS           = Int16(0)
const SQL_SUCCESS_WITH_INFO = Int16(1)

# error codes
const SQL_ERROR             = Int16(-1)
const SQL_INVALID_HANDLE    = Int16(-2)

# status codes
const SQL_STILL_EXECUTING   = Int16(2)
const SQL_NO_DATA           = Int16(100)

const RETURN_VALUES = Dict(SQL_ERROR   => "SQL_ERROR",
                           SQL_NO_DATA => "SQL_NO_DATA",
                           SQL_INVALID_HANDLE  => "SQL_INVALID_HANDLE",
                           SQL_STILL_EXECUTING => "SQL_STILL_EXECUTING")

macro odbc(func,args,vals...)
    if Sys.iswindows()
        esc(quote
            ret = ccall( ($func, odbc_dm), stdcall, SQLRETURN, $args, $(vals...))
        end)
    else
        esc(quote
            ret = ccall( ($func, odbc_dm),          SQLRETURN, $args, $(vals...))
        end)
    end
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712400(v=vs.85).aspx"
function SQLDrivers(env::Ptr{Void},
                    dir,
                    driver_desc::Ptr{SQLWCHAR},
                    desclen,
                    desc_length::Ref{SQLSMALLINT},
                    driver_attr::Ptr{SQLWCHAR},
                    attrlen,
                    attr_length::Ref{SQLSMALLINT})
    @odbc(:SQLDriversW,
                (Ptr{Void}, SQLUSMALLINT, Ptr{SQLWCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}, Ptr{SQLWCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}),
                env, dir, driver_desc, desclen, desc_length, driver_attr, attrlen, attr_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711004(v=vs.85).aspx"
function SQLDataSources(env::Ptr{Void},
                        dir,
                        dsn_desc::Ptr{SQLWCHAR},
                        desclen,
                        desc_length::Ref{SQLSMALLINT},
                        dsn_attr::Ptr{SQLWCHAR},
                        attrlen,
                        attr_length::Ref{SQLSMALLINT})
    @odbc(:SQLDataSourcesW,
                (Ptr{Void}, SQLUSMALLINT, Ptr{SQLWCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}, Ptr{SQLWCHAR}, SQLSMALLINT, Ref{SQLSMALLINT}),
                env, dir, dsn_desc, desclen, desc_length, dsn_attr, attrlen, attr_length)
end

#### Handle Functions ####

# SQLAllocHandle
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712455(v=vs.85).aspx"
# Description: allocates an environment, connection, statement, or descriptor handle
# Valid handle types
const SQL_HANDLE_ENV  = Int16(1)
const SQL_HANDLE_DBC  = Int16(2)
const SQL_HANDLE_STMT = Int16(3)
const SQL_HANDLE_DESC = Int16(4)
const SQL_NULL_HANDLE = C_NULL

#Status: Tested on Windows, Linux, Mac 32/64-bit
function SQLAllocHandle(handletype::SQLSMALLINT, parenthandle::Ptr{Void}, handle::Ref{Ptr{Void}})
    @odbc(:SQLAllocHandle,
                (SQLSMALLINT, Ptr{Void}, Ref{Ptr{Void}}),
                handletype, parenthandle, handle)
end

# SQLFreeHandle
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710123(v=vs.85).aspx"
# Description: frees resources associated with a specific environment, connection, statement, or descriptor handle
# See SQLAllocHandle for valid handle types
# Status: Tested on Windows, Linux, Mac 32/64-bit
function SQLFreeHandle(handletype::SQLSMALLINT,handle::Ptr{Void})
    @odbc(:SQLFreeHandle,
                (SQLSMALLINT, Ptr{Void}),
                handletype, handle)
end

# SQLSetEnvAttr
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709285(v=vs.85).aspx"
# Description: sets attributes that govern aspects of environments
# Valid attributes; valid values for attribute are indented
const SQL_ATTR_CONNECTION_POOLING = 201
const SQL_CP_OFF = UInt(0)
const SQL_CP_ONE_PER_DRIVER = UInt(1)
const SQL_CP_ONE_PER_HENV = UInt(2)
const SQL_CP_DEFAULT = SQL_CP_OFF
const SQL_ATTR_CP_MATCH = 202
const SQL_CP_RELAXED_MATCH = UInt(1)
const SQL_CP_STRICT_MATCH = UInt(0)
const SQL_ATTR_ODBC_VERSION = 200
const SQL_OV_ODBC2 = 2
const SQL_OV_ODBC3 = 3
const SQL_ATTR_OUTPUT_NTS = 10001
const SQL_TRUE = 1
const SQL_FALSE = 0

#Status: Tested on Windows, Linux, Mac 32/64-bit
function SQLSetEnvAttr(env_handle::Ptr{Void}, attribute::Int, value::T) where {T<:Union{Int,UInt}}
    @odbc(:SQLSetEnvAttr,
                (Ptr{Void}, Int, T, Int), env_handle, attribute, value, 0)
end

# SQLGetEnvAttr
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709276(v=vs.85).aspx"
# Description: returns the current setting of an environment attribute
# Valid attributes: See SQLSetEnvAttr
# Status:
function SQLGetEnvAttr(env::Ptr{Void},attribute::Int,value::Array{Int,1},bytes_returned::Array{Int,1})
    @odbc(:SQLGetEnvAttr,
                (Ptr{Void}, Int, Ptr{Int}, Int, Ptr{Int}),
                env, attribute, value, 0, bytes_returned)
end

# SQLSetConnectAttr
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713605(v=vs.85).aspx"
# Description: sets attributes that govern aspects of connections.
# Valid attributes
const SQL_ATTR_ACCESS_MODE = 101
const SQL_MODE_READ_ONLY = UInt(1)
const SQL_MODE_READ_WRITE = UInt(0)
#const SQL_ATTR_ASYNC_DBC_EVENT
#pointer
#const SQL_ATTR_ASYNC_DBC_FUNCTIONS_ENABLE
#const SQL_ASYNC_DBC_ENABLE_ON = UInt()
#const SQL_ASYNC_DBC_ENABLE_OFF = UInt()
#const SQL_ATTR_ASYNC_DBC_PCALLBACK
#pointer
#const SQL_ATTR_ASYNC_DBC_PCONTEXT
#pointer
const SQL_ATTR_ASYNC_ENABLE = 4
const SQL_ASYNC_ENABLE_OFF = UInt(0)
const SQL_ASYNC_ENABLE_ON = UInt(1)
const SQL_ATTR_AUTOCOMMIT = 102
const SQL_AUTOCOMMIT_OFF = UInt(0)
const SQL_AUTOCOMMIT_ON = UInt(1)
const SQL_ATTR_CONNECTION_TIMEOUT = 113
#uint of how long you want the connection timeout
const SQL_ATTR_CURRENT_CATALOG = 109
#string/Ptr{UInt8} of default database to use
#const SQL_ATTR_DBC_INFO_TOKEN
#pointer
const SQL_ATTR_ENLIST_IN_DTC = 1207
#pointer: Pass a DTC OLE transaction object that specifies the transaction to export to
# SQL Server, or SQL_DTC_DONE to end the connection's DTC association.
const SQL_ATTR_LOGIN_TIMEOUT = 103
#uint of how long you want the login timeout
const SQL_ATTR_METADATA_ID = 10014
#SQL_TRUE, SQL_FALSE
const SQL_ATTR_ODBC_CURSORS = 110
const SQL_CUR_USE_IF_NEEDED = UInt(0)
const SQL_CUR_USE_ODBC = UInt(1)
const SQL_CUR_USE_DRIVER = UInt(2)
const SQL_ATTR_PACKET_SIZE = 112
#uint for network packet size
const SQL_ATTR_QUIET_MODE = 111
#window handle pointer
const SQL_ATTR_TRACE = 104
const SQL_OPT_TRACE_OFF = UInt(0)
const SQL_OPT_TRACE_ON = UInt(1)
const SQL_ATTR_TRACEFILE = 105
#A null-terminated character string containing the name of the trace file.
const SQL_ATTR_TRANSLATE_LIB = 106
# A null-terminated character string containing the name of a library containing the functions SQLDriverToDataSource and
# SQLDataSourceToDriver that the driver accesses to perform tasks such as character set translation.
const SQL_ATTR_TRANSLATE_OPTION = 107
#A 32-bit flag value that is passed to the translation DLL.
const SQL_ATTR_TXN_ISOLATION = 108
#A 32-bit bitmask that sets the transaction isolation level for the current connection.

#Valid value_length
const SQL_IS_POINTER = -4
const SQL_IS_INTEGER = -6
const SQL_IS_UINTEGER = -5
const SQL_NTS = -3

#length of string or binary stream
#Status:
function SQLSetConnectAttr(dbc::Ptr{Void},attribute::Int,value::UInt,value_length::Int)
    @odbc(:SQLSetConnectAttrW,
                (Ptr{Void},Int,UInt,Int),
                dbc,attribute,value,value_length)
end

function SQLSetConnectAttr(dbc::Ptr{Void},attribute::Int,value::Array{Int},value_length::Int)
    @odbc(:SQLSetConnectAttrW,
                (Ptr{Void},Int,Ptr{Int},Int),
                dbc,attribute,value,value_length)
end

#SQLGetConnectAttr
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710297(v=vs.85).aspx"
#Description: returns the current setting of a connection attribute.
#Valid attributes: see SQLSetConnectAttr in addition to those below
const SQL_ATTR_AUTO_IPD = 10001
#SQL_TRUE, SQL_FALSE
const SQL_ATTR_CONNECTION_DEAD = 1209
const SQL_CD_TRUE = 1
const SQL_CD_FALSE = 0
#Status:
function SQLGetConnectAttr(dbc::Ptr{Void},attribute::Int,value::Array{T,N},bytes_returned::Array{Int,1}) where {T,N}
    @odbc(:SQLGetConnectAttrW,
                (Ptr{Void},Int,Ptr{T},Int,Ptr{Int}),
                dbc,attribute,value,sizeof(T)*N,bytes_returned)
end

#SQLSetStmtAttr
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712631(v=vs.85).aspx"
#Description: sets attributes related to a statement.
#Valid attributes
const SQL_ATTR_ROW_STATUS_PTR = 25
const SQL_ATTR_ROWS_FETCHED_PTR  = 26
const SQL_ATTR_ROW_ARRAY_SIZE = 27
const SQL_ATTR_CURSOR_TYPE = 6
const SQL_ATTR_CURSOR_SCROLLABLE = -1
const SQL_NONSCROLLABLE = 0
const SQL_SCROLLABLE = 1
const SQL_CURSOR_DYNAMIC = UInt(2)
const SQL_CURSOR_STATIC = UInt(3)
#this sets the rowset size for ExtendedFetch and FetchScroll
#Valid value_length: See SQLSetConnectAttr; SQL_IS_POINTER, SQL_IS_INTEGER, SQL_IS_UINTEGER, SQL_NTS
#Status:
function SQLSetStmtAttr(stmt::Ptr{Void},attribute,value::Ref{SQLLEN},value_length)
    @odbc(:SQLSetStmtAttrW,
                (Ptr{Void},SQLINTEGER,Ref{SQLLEN},SQLINTEGER),
                stmt,attribute,value,value_length)
end

function SQLSetStmtAttr(stmt::Ptr{Void},attribute,value,value_length)
    @odbc(:SQLSetStmtAttrW,
                (Ptr{Void},SQLINTEGER,SQLULEN,SQLINTEGER),
                stmt,attribute,value,value_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715438(v=vs.85).aspx"
function SQLGetStmtAttr(stmt::Ptr{Void},attribute::Int,value::Array{T,N},bytes_returned::Array{Int,1}) where {T,N}
    @odbc(:SQLGetStmtAttrW,
                (Ptr{Void},Int,Ptr{T},Int,Ptr{Int}),
                stmt,attribute,value,sizeof(T)*N,bytes_returned)
end

#SQLFreeStmt
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709284(v=vs.85).aspx"
#Description: stops processing associated with a specific statement,
# closes any open cursors associated with the statement,
# discards pending results, or, optionally,
# frees all resources associated with the statement handle.
#Valid param
const SQL_CLOSE = UInt16(0)
const SQL_RESET_PARAMS = UInt16(3)
const SQL_UNBIND = UInt16(2)

#Status:
function SQLFreeStmt(stmt::Ptr{Void},param::UInt16)
    @odbc(:SQLFreeStmt,
                (Ptr{Void},UInt16),
                stmt, param)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713560(v=vs.85).aspx"
function SQLSetDescField(desc::Ptr{Void},i::Int16,field_id::Int16,value::Array{T,N}) where {T,N}
    @odbc(:SQLSetDescFieldW,
                (Ptr{Void},Int16,Int16,Ptr{T},Int),
                desc,i,field_id,value,length(value))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716370(v=vs.85).aspx"
function SQLGetDescField(desc::Ptr{Void},i::Int16,attribute::Int16,value::Array{T,N},bytes_returned::Array{Int,1}) where {T,N}
    @odbc(:SQLGetDescFieldW,
                (Ptr{Void},Int16,Int16,Ptr{T},Int,Ptr{Int}),
                desc,i,attribute,value,sizeof(T)*N,bytes_returned)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710921(v=vs.85).aspx"
function SQLGetDescRec(desc::Ptr{Void},i::Int16,name::Array{UInt8,1},name_length::Array{Int16,1},type_ptr::Array{Int16,1},subtype_ptr::Array{Int16,1},length_ptr::Array{Int,1},precision_ptr::Array{Int16,1},scale_ptr::Array{Int16,1},nullable_ptr::Array{Int16,1},)
    @odbc(:SQLGetDescRecW,
                (Ptr{Void},Int16,Ptr{UInt8},Int16,Ptr{Int16},Ptr{Int16},Ptr{Int16},Ptr{Int},Ptr{Int16},Ptr{Int16},Ptr{Int16}),
                desc,i,name,length(name),name_length,type_ptr,subtype_ptr,length_ptr,precision_ptr,scale_ptr,nullable_ptr)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715378(v=vs.85).aspx"
function SQLCopyDesc(source_desc::Ptr{Void},dest_desc::Ptr{Void})
    @odbc(:SQLCopyDesc,
                (Ptr{Void},Ptr{Void}),
                source_desc,dest_desc)
end

### Connection Functions ###
# SQLConnect
#"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711810(v=vs.85).aspx"
# Description: establishes connections to a driver and a data source
# Status:
function SQLConnect(dbc::Ptr{Void},dsn,username,password)
    @odbc(:SQLConnectW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                dbc,transcode(SQLWCHAR,dsn),length(transcode(SQLWCHAR,dsn)),transcode(SQLWCHAR,username),length(transcode(SQLWCHAR,username)),transcode(SQLWCHAR,password),length(transcode(SQLWCHAR,password)))
end

#SQLDriverConnect
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715433(v=vs.85).aspx"
#Description:
#Valid driver_prompt
const SQL_DRIVER_COMPLETE = UInt16(1)
const SQL_DRIVER_COMPLETE_REQUIRED = UInt16(3)
const SQL_DRIVER_NOPROMPT = UInt16(0)
const SQL_DRIVER_PROMPT = UInt16(2)
#Status:
function SQLDriverConnect(dbc::Ptr{Void},window_handle::Ptr{Void},conn_string,out_conn::Ptr{SQLWCHAR},out_len,out_buff::Ref{Int16},driver_prompt)
    @odbc(:SQLDriverConnectW,
                (Ptr{Void},Ptr{Void},Ptr{SQLWCHAR},SQLSMALLINT,Ptr{SQLWCHAR},SQLSMALLINT,Ptr{SQLSMALLINT},SQLUSMALLINT),
                dbc,window_handle,transcode(SQLWCHAR,conn_string),length(transcode(SQLWCHAR,conn_string)),out_conn,out_len,out_buff,driver_prompt)
end
#SQLBrowseConnect
 "http://msdn.microsoft.com/en-us/library/windows/desktop/ms714565(v=vs.85).aspx"
 #Description: supports an iterative method of discovering and enumerating the attributes and attribute values required to connect to a data source
 #Status:
function SQLBrowseConnect(dbc::Ptr{Void},instring::AbstractString,outstring::Array{SQLWCHAR,1},indicator::Array{Int16,1})
    @odbc(:SQLBrowseConnectW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{Int16}),
                dbc,transcode(SQLWCHAR,instring),length(transcode(SQLWCHAR,instring)),transcode(SQLWCHAR,outstring),length(transcode(SQLWCHAR,outstring)),indicator)
end
#SQLDisconnect
 "http://msdn.microsoft.com/en-us/library/windows/desktop/ms713946(v=vs.85).aspx"
 #Description: closes the connection associated with a specific connection handle
 #Status:
function SQLDisconnect(dbc::Ptr{Void})
    @odbc(:SQLDisconnect,
                (Ptr{Void},),
                dbc)
end
#SQLGetFunctions
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709291(v=vs.85).aspx"
#Descriptions:
#Valid functionid

#supported will be SQL_TRUE or SQL_FALSE
#Status:
function SQLGetFunctions(dbc::Ptr{Void},functionid::UInt16,supported::Array{UInt16,1})
    @odbc(:SQLGetFunctions,
                (Ptr{Void},UInt16,Ptr{UInt16}),
                dbc,functionid,supported)
end

#SQLGetInfo
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711681(v=vs.85).aspx"
#Description:
#Status:
function SQLGetInfo(dbc::Ptr{Void},attribute::Int,value::Array{T,N},bytes_returned::Array{Int,1}) where {T,N}
    @odbc(:SQLGetInfoW,
                (Ptr{Void},Int,Ptr{T},Int,Ptr{Int}),
                dbc,attribute,value,sizeof(T)*N,bytes_returned)
end

#### Query Functions ####
#SQLNativeSql
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714575(v=vs.85).aspx"
#Description: returns the SQL string as modified by the driver
#Status:
function SQLNativeSql(dbc::Ptr{Void},query_string::AbstractString,output_string::Array{SQLWCHAR,1},length_ind::Array{Int,1})
    @odbc(:SQLNativeSql,
                (Ptr{Void},Ptr{SQLWCHAR},Int,Ptr{SQLWCHAR},Int,Ptr{Int}),
                dbc,transcode(SQLWCHAR,query_string),length(transcode(SQLWCHAR,query_string)),output_string,length(output_string),length_ind)
end

#SQLGetTypeInfo
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714632(v=vs.85).aspx"
#Description:
#valid sqltype
#const SQL_ALL_TYPES =
#Status:
function SQLGetTypeInfo(stmt::Ptr{Void},sqltype::Int16)
    @odbc(:SQLGetTypeInfo,
                (Ptr{Void},Int16),
                stmt,sqltype)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713824(v=vs.85).aspx"
function SQLPutData(stmt::Ptr{Void},data::Array{T},data_length::Int) where {T}
    @odbc(:SQLPutData,
                (Ptr{Void},Ptr{T},Int),
                stmt,data,data_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710926(v=vs.85).aspx"
function SQLPrepare(stmt::Ptr{Void},query_string::AbstractString)
    @odbc(:SQLPrepareW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,query_string),length(transcode(SQLWCHAR,query_string)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713584(v=vs.85).aspx"
function SQLExecute(stmt::Ptr{Void})
    @odbc(:SQLExecute,
                (Ptr{Void},),
                stmt)
end

#SQLExecDirect
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713611(v=vs.85).aspx"
#Description: executes a preparable statement
#Status:
function SQLExecDirect(stmt::Ptr{Void},query::AbstractString)
    @odbc(:SQLExecDirectW,
                (Ptr{Void},Ptr{SQLWCHAR},Int),
                stmt,transcode(SQLWCHAR,query),length(transcode(SQLWCHAR,query)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714112(v=vs.85).aspx"
function SQLCancel(stmt::Ptr{Void})
    @odbc(:SQLCancel,
                (Ptr{Void},),
                stmt)
end

#### Resultset Metadata Functions ####
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715393(v=vs.85).aspx"
function SQLNumResultCols(stmt::Ptr{Void},cols::Ref{Int16})
    @odbc(:SQLNumResultCols,
                (Ptr{Void},Ref{Int16}),
                stmt, cols)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711835(v=vs.85).aspx"
function SQLRowCount(stmt::Ptr{Void},rows::Ref{Int})
    @odbc(:SQLRowCount,
                (Ptr{Void},Ref{Int}),
                stmt, rows)
end

# "http://msdn.microsoft.com/en-us/library/windows/desktop/ms713558(v=vs.85).aspx"
# function SQLColAttribute(stmt::Ptr{Void},x::Int,)
# = @odbc(:SQLColAttributeW,
#                 (Ptr{Void},UInt16,UInt16,Ptr,Int16,Ptr{Int16},Ptr{Int}),
#                 stmt,x,)
# end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716289(v=vs.85).aspx"
function SQLDescribeCol(stmt,x,nm::Ptr{SQLWCHAR},nmlen,len::Ref{Int16},dt::Ref{Int16},cs::Ref{SQLULEN},dd::Ref{Int16},nul::Ref{Int16})
    @odbc(:SQLDescribeColW,
                (Ptr{Void},SQLUSMALLINT,Ptr{SQLWCHAR},SQLSMALLINT,Ref{SQLSMALLINT},Ref{SQLSMALLINT},Ref{SQLULEN},Ref{SQLSMALLINT},Ref{SQLSMALLINT}),
                stmt,x,nm,nmlen,len,dt,cs,dd,nul)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710188(v=vs.85).aspx"
function SQLDescribeParam(stmt::Ptr{Void},x::Int,sqltype::Array{Int16,1},column_size::Array{Int,1},decimal_digits::Array{Int16,1},nullable::Array{Int16,1})
    @odbc(:SQLDescribeParam,
                (Ptr{Void},UInt16,Ptr{Int16},Ptr{Int},Ptr{Int16},Ptr{Int16}),
                stmt,x,sqltype,column_size,decimal_digits,nullable)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712366(v=vs.85).aspx"
function SQLParamData(stmt::Ptr{Void},ptr_buffer::Array{Ptr{Void},1})
    @odbc(:SQLParamData,
                (Ptr{Void},Ptr{Void}),
                stmt,ptr_buffer)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715409(v=vs.85).aspx"
function SQLNumParams(stmt::Ptr{Void},param_count::Array{Int16,1})
    @odbc(:SQLNumParams,
                (Ptr{Void},Ptr{Int16}),
                stmt,param_count)
end

#### Resultset Retrieval Functions ####
#SQLBindParameter
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710963(v=vs.85).aspx"
#Description:
#valid iotype
const SQL_PARAM_INPUT = Int16(1)
const SQL_PARAM_OUTPUT = Int16(4)
const SQL_PARAM_INPUT_OUTPUT = Int16(2)
#const SQL_PARAM_INPUT_OUTPUT_STREAM = Int16()
#const SQL_PARAM_OUTPUT_STREAM = Int16()
#Status:
function SQLBindParameter(stmt::Ptr{Void},x::Int,iotype::Int16,ctype::Int16,sqltype::Int16,column_size::Int,decimal_digits::Int,param_value,param_size::Int,len::Ref{SQLLEN})
    @odbc(:SQLBindParameter,
                (Ptr{Void},UInt16,Int16,Int16,Int16,UInt,Int16,Ptr{Void},Int,Ptr{SQLLEN}),
                stmt,x,iotype,ctype,sqltype,column_size,decimal_digits,param_value,param_size,len)
end
SQLSetParam = SQLBindParameter

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711010(v=vs.85).aspx"
function SQLBindCols(stmt::Ptr{Void},x,ctype,mem,jlsize,indicator::Vector{SQLLEN})
    @odbc(:SQLBindCol,
                (Ptr{Void},SQLUSMALLINT,SQLSMALLINT,Ptr{Void},SQLLEN,Ptr{SQLLEN}),
                stmt,x,ctype,mem,jlsize,indicator)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711707(v=vs.85).aspx"
function SQLSetCursorName(stmt::Ptr{Void},cursor::AbstractString)
    @odbc(:SQLSetCursorNameW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,cursor),length(transcode(SQLWCHAR,cursor)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716209(v=vs.85).aspx"
function SQLGetCursorName(stmt::Ptr{Void},cursor::Array{UInt8,1},cursor_length::Array{Int16,1})
    @odbc(:SQLGetCursorNameW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{Int16}),
                stmt,transcode(SQLWCHAR,cursor),length(transcode(SQLWCHAR,cursor)),cursor_length)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715441(v=vs.85).aspx"
function SQLGetData(stmt::Ptr{Void},i,ctype,mem,jlsize,indicator::Ref{SQLLEN})
    @odbc(:SQLGetData,
                (Ptr{Void},SQLUSMALLINT,SQLSMALLINT,Ptr{Void},SQLLEN,Ptr{SQLLEN}),
                stmt,i,ctype,mem,jlsize,indicator)
end

#SQLFetchScroll
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714682(v=vs.85).aspx"
#Description:
#valid fetch_orientation
const SQL_FETCH_NEXT = Int16(1)
const SQL_FETCH_PRIOR = Int16(4)
const SQL_FETCH_FIRST = Int16(2)
const SQL_FETCH_LAST = Int16(3)
const SQL_FETCH_ABSOLUTE = Int16(5)
const SQL_FETCH_RELATIVE = Int16(6)
const SQL_FETCH_BOOKMARK = Int16(8)
#Status:
function SQLFetchScroll(stmt::Ptr{Void},fetch_orientation::Int16,fetch_offset::Int)
    @odbc(:SQLFetchScroll,
                (Ptr{Void},Int16,Int),
                stmt,fetch_orientation,fetch_offset)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713591(v=vs.85).aspx"
function SQLExtendedFetch(stmt::Ptr{Void},fetch_orientation::UInt16,fetch_offset::Int,row_count_ptr::Array{Int,1},row_status_array::Array{Int16,1})
    @odbc(:SQLExtendedFetch,
                (Ptr{Void},UInt16,Int,Ptr{Int},Ptr{Int16}),
                stmt,fetch_orientation,fetch_offset,row_count_ptr,row_status_array)
end

#SQLSetPos
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713507(v=vs.85).aspx"
#Description:
#valid operation
const SQL_POSITION = UInt16(0) #SQLSetPos
const SQL_REFRESH = UInt16(1) #SQLSetPos
const SQL_UPDATE = UInt16(2) #SQLSetPos
const SQL_DELETE = UInt16(3) #SQLSetPos
#valid lock_type
const SQL_LOCK_NO_CHANGE = UInt16(0) #SQLSetPos
const SQL_LOCK_EXCLUSIVE = UInt16(1) #SQLSetPos
const SQL_LOCK_UNLOCK = UInt16(2) #SQLSetPos
#Status
function SQLSetPos(stmt::Ptr{Void},rownumber::T,operation::UInt16,lock_type::UInt16) where {T}
    @odbc(:SQLSetPos,
                (Ptr{Void},T,UInt16,UInt16),
                stmt,rownumber,operation,lock_type)
end #T can be Uint64 or UInt16 it seems

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714673(v=vs.85).aspx"
function SQLMoreResults(stmt::Ptr{Void})
    @odbc(:SQLMoreResults,
                (Ptr{Void},),
                stmt)
end

#SQLEndTran
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716544(v=vs.85).aspx"
#Description:
#valid completion_type
const SQL_COMMIT = Int16(0) #SQLEndTran
const SQL_ROLLBACK = Int16(1) #SQLEndTran
#Status:
function SQLEndTran(handletype::Int16,handle::Ptr{Void},completion_type::Int16)
    @odbc(:SQLEndTran,
                (Int16,Ptr{Void},Int16),
                handletype,handle,completion_type)
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709301(v=vs.85).aspx"
function SQLCloseCursor(stmt::Ptr{Void})
    @odbc(:SQLCloseCursor,
                (Ptr{Void},),
                stmt)
end

#SQLBulkOperations
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms712471(v=vs.85).aspx"
#Description:
#valid operation
const SQL_ADD = UInt16(4) #SQLBulkOperations
const SQL_UPDATE_BY_BOOKMARK = UInt16(5) #SQLBulkOperations
const SQL_DELETE_BY_BOOKMARK = UInt16(6) #SQLBulkOperations
const SQL_FETCH_BY_BOOKMARK = UInt16(7) #SQLBulkOperations
#Status:
function SQLBulkOperations(stmt::Ptr{Void},operation::UInt16)
    @odbc(:SQLBulkOperations,
                (Ptr{Void},UInt16),
                stmt,operation)
end

#### DBMS Meta Functions ####
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711683(v=vs.85).aspx"
function SQLColumns(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,column::AbstractString)
    @odbc(:SQLColumnsW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),transcode(SQLWCHAR,column),length(transcode(SQLWCHAR,column)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716336(v=vs.85).aspx"
function SQLColumnPrivileges(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,column::AbstractString)
    @odbc(:SQLColumnPrivilegesW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),transcode(SQLWCHAR,column),length(transcode(SQLWCHAR,column)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms709315(v=vs.85).aspx"
function SQLForeignKeys(stmt::Ptr{Void},pkcatalog::AbstractString,pkschema::AbstractString,pktable::AbstractString,fkcatalog::AbstractString,fkschema::AbstractString,fktable::AbstractString)
    @odbc(:SQLForeignKeysW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,pkcatalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,pkschema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,pktable)),transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,fkcatalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,fkschema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,fktable)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711005(v=vs.85).aspx"
function SQLPrimaryKeys(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString)
    @odbc(:SQLPrimaryKeysW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711701(v=vs.85).aspx"
function SQLProcedureColumns(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,proc::AbstractString,column::AbstractString)
    @odbc(:SQLProcedureColumnsW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),proc,length(proc),transcode(SQLWCHAR,column),length(transcode(SQLWCHAR,column)))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms715368(v=vs.85).aspx"
function SQLProcedures(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,proc::AbstractString)
    @odbc(:SQLProceduresW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),proc,length(proc))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711831(v=vs.85).aspx"
function SQLTables(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,table_type::AbstractString)
    @odbc(:SQLTablesW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),table_type,length(table_type))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms713565(v=vs.85).aspx"
function SQLTablePrivileges(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString)
    @odbc(:SQLTablePrivilegesW,
                (Ptr{Void},Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)))
end

#SQLStatistics
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms711022(v=vs.85).aspx"
#Description:
#valid unique
const SQL_INDEX_ALL = UInt16(1)
const SQL_INDEX_CLUSTERED = UInt16(1)
const SQL_INDEX_HASHED = UInt16(2)
const SQL_INDEX_OTHER = UInt16(3)
const SQL_INDEX_UNIQUE = UInt16(0)
#valid reserved
const SQL_ENSURE = UInt16(1)
const SQL_QUICK = UInt16(0)
#Status:
function SQLStatistics(stmt::Ptr{Void},catalog::AbstractString,schema::AbstractString,table::AbstractString,unique::UInt16,reserved::UInt16)
    @odbc(:SQLStatisticsW,
                (Ptr{Void},Ptr{UInt8},Int16,Ptr{UInt8},Int16,Ptr{UInt8},Int16,UInt16,UInt16),
                stmt,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),unique,reserved)
end

#SQLSpecialColumns
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms714602(v=vs.85).aspx"
#Description:
#valid id_type
const SQL_BEST_ROWID        = Int16(1) #SQLSpecialColumns
const SQL_ROWVER            = Int16(2) #SQLSpecialColumns
#valid scope
const SQL_SCOPE_CURROW      = Int16(0) #SQLSpecialColumns
const SQL_SCOPE_SESSION     = Int16(2) #SQLSpecialColumns
const SQL_SCOPE_TRANSACTION = Int16(1) #SQLSpecialColumns
#valid nullable
const SQL_NO_NULLS          = Int16(0) #SQLSpecialColumns
const SQL_NULLABLE          = Int16(1) #SQLSpecialColumns
#const SQL_NULLABLE_UNKNOWN = Int16() #SQLSpecialColumns
#Status:
function SQLSpecialColumns(stmt::Ptr{Void},id_type::Int16,catalog::AbstractString,schema::AbstractString,table::AbstractString,scope::Int16,nullable::Int16)
    @odbc(:SQLSpecialColumnsW,
                (Ptr{Void},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Ptr{SQLWCHAR},Int16,Int16,Int16),
                stmt,id_type,transcode(SQLWCHAR,catalog),length(transcode(SQLWCHAR,catalog)),transcode(SQLWCHAR,schema),length(transcode(SQLWCHAR,schema)),transcode(SQLWCHAR,table),length(transcode(SQLWCHAR,table)),scope,nullable)
end

#### Error Handling Functions ####
#TODO: add consts
"http://msdn.microsoft.com/en-us/library/windows/desktop/ms710181(v=vs.85).aspx"
function SQLGetDiagField(handletype::Int16,handle::Ptr{Void},i::Int16,diag_id::Int16,diag_info::Array{SQLWCHAR,1},buffer_length::Int16,diag_length::Array{Int16,1})
    @odbc(:SQLGetDiagFieldW,
                (Int16,Ptr{Void},Int16,Int16,Ptr{SQLWCHAR},Int16,Ptr{Int16}),
                handletype,handle,i,diag_id,transcode(SQLWCHAR,diag_info),buffer_length,transcode(SQLWCHAR,msg_length))
end

"http://msdn.microsoft.com/en-us/library/windows/desktop/ms716256(v=vs.85).aspx"
function SQLGetDiagRec(handletype,handle,i,state::Ptr{SQLWCHAR},native::Ref{SQLINTEGER},error_msg,errlen,msg_length)
    @odbc(:SQLGetDiagRecW,
                (SQLSMALLINT,Ptr{Void},SQLSMALLINT,Ptr{SQLWCHAR},Ref{SQLINTEGER},Ptr{SQLWCHAR},SQLSMALLINT,Ref{SQLSMALLINT}),
                handletype,handle,i,state,native,error_msg,errlen,msg_length)
end

end # module

ext: jl


text: map(x->print(length(x)),split(replace("Now I need a drink, alcoholic of course, after the heavy lectures involving quantum mechanics.",r",|\.", "")))


text:
using Base.Test

tests = [
    "candlesticks",
    "levels",
    "momentum",
    "movingaverages",
    "utilities",
    "volatility",
    "volume",
]

@testset "MarketTechnicals" begin
    println("Running tests:")

    for test ∈ tests
        println("\t* $test ...")
        include("$test.jl")
    end
end


text:
# GPU CuArray Blueprint
function gpu_gen_code_replica!(
    code_replica::CuArray{ComplexF32},
    system::AbstractGNSSSystem,
    code_frequency,
    sampling_frequency,
    start_code_phase::AbstractFloat,
    start_sample::Integer,
    num_samples::Integer,
    early_late_sample_shift,
    prn::Integer
)
    # code_replica] = GNSSSignals.get_code(
    #         system,
    #         code_frequency .* (start_sample:start_sample + num_samples + 2*early_late_sample_shift) ./ sampling_frequency .+ start_code_phase,
    #         prn
    #     )
    idxs = start_sample:start_sample - 1 + num_samples + 2*early_late_sample_shift
    phases = code_frequency .* (0:num_samples - 1 + 2 * early_late_sample_shift) ./ sampling_frequency .+ start_code_phase
    code_length = get_code_length(system) * get_secondary_code_length(system)
    @inbounds @views code_replica[idxs] .= system.codes[2 .+ mod.(floor.(Int, phases), code_length), prn]
end

# CPU Array Blueprint
function gpu_gen_code_replica!(
    code_replica::Array,
    system::AbstractGNSSSystem,
    code_frequency,
    sampling_frequency,
    start_code_phase::AbstractFloat,
    start_sample::Integer,
    num_samples::Integer,
    early_late_sample_shift,
    prn::Integer
)
    # code_replica = GNSSSignals.get_code(
    #         system,
    #         code_frequency .* (start_sample:start_sample + num_samples + 2*early_late_sample_shift) ./ sampling_frequency .+ start_code_phase,
    #         prn
    #     )
    idxs = start_sample:start_sample - 1 + num_samples + 2*early_late_sample_shift
    phases = code_frequency .* (0:num_samples - 1 + 2 * early_late_sample_shift) ./ sampling_frequency .+ start_code_phase
    code_length = get_code_length(system) * get_secondary_code_length(system)
    @inbounds @views code_replica[idxs] .= system.codes[2 .+ mod.(floor.(Int, phases), code_length), prn]
end


text:
using MonteCarloObservable
using Test, Statistics
import HDF5

@testset "All Tests" begin

    include("observable.jl")
    include("lightobservable.jl")

end


text:
""""""
function preprocess(txt::AbstractString, lang::Language)
    txt = replace(txt, r"[.:,;?!()\[\]=*/+-]" => " ")
    tokens = tokenize(lang, txt)
    # TODO Stop word removal
    stemmer = Stemmer(isocode(lang))
    tokens = stem(stemmer, tokens)
    return join(tokens, " ")
end
preprocess(txt::AbstractString) = preprocess(txt, LanguageDetector()(txt)[1])

""""""
function TextAnalysis.tf_idf(crps::Corpus)
    update_lexicon!(crps)
    dtm = DocumentTermMatrix(crps)
    return tf_idf(dtm)
end
TextAnalysis.tf_idf(txts::Vector{<:AbstractString}) =
    tf_idf((Corpus(StringDocument.(txts))))

""""""
TextAnalysis.lsa(M::SparseMatrixCSC) = svd(Array(M))


text:
using GLMakie, LibFaceDetection, FileIO
path = download("https://thumbs.dreamstime.com/z/many-faces-2754451.jpg", "faces.jpg")
img = load(path)
faces = detect_faces(img)
scene = image(rotr90(img), scale_plot=false)
for face in faces
    linesegments!(face.boundingbox)
    scatter!(face.points)
end
display(scene)


text:
################################################################################
#
#  Unit groups of orders in étale algebras
#
################################################################################

function unit_group(O::AlgAssAbsOrd)
  @assert iscommutative(O)
  mU = get_attribute!(O, :unit_group) do
    if ismaximal(O)
      U, mU = _unit_group_maximal(O)
    else
      U, mU = _unit_group_non_maximal(O)
    end
    return mU
  end::Map
  return domain(mU), mU
end

function _unit_group_maximal(O::AlgAssAbsOrd)
  A = algebra(O)
  fields_and_maps = as_number_fields(A)
  unit_groups = Tuple{GrpAbFinGen, MapUnitGrp{NfOrd}}[ unit_group(maximal_order(field)) for (field, map) in fields_and_maps ]
  G = unit_groups[1][1]
  for i = 2:length(unit_groups)
    G = direct_product(G, unit_groups[i][1], task = :none)::GrpAbFinGen
  end
  S, StoG = snf(G)

  local disc_exp
  let StoG = StoG, unit_groups = unit_groups, fields_and_maps = fields_and_maps
    function disc_exp(x::GrpAbFinGenElem)
      g = StoG(x)
      v = zero(O)
      offset = 1
      for i = 1:length(fields_and_maps)
        K, AtoK = fields_and_maps[i]
        U, mU = unit_groups[i]
        u = U(sub(g.coeff, 1:1, offset:(offset + ngens(U) - 1)))
        v += O(AtoK\elem_in_nf(mU(u), copy = false))
        offset += ngens(U)
      end
      return v
    end
  end

  local disc_log
  let fields_and_maps = fields_and_maps, unit_groups = unit_groups, StoG = StoG, G = G
    function disc_log(x::AlgAssAbsOrdElem)
      g = zero_matrix(FlintZZ, 1, 0)
      for i = 1:length(fields_and_maps)
        K, AtoK = fields_and_maps[i]
        U, mU = unit_groups[i]
        OK = codomain(mU)
        y = OK(AtoK(elem_in_algebra(x, copy = false)))
        @assert isunit(y)
        u = mU\y
        g = hcat(g, u.coeff)
      end
      return StoG\G(g)
    end
  end

  StoO = MapUnitGrp{typeof(O)}()
  StoO.header = MapHeader(S, O, disc_exp, disc_log)
  return S, StoO
end

# Returns the group (O_A/F)^\times/(O/F)^\times and a map from this group to
# O_A/F where F is the conductor of O in the maximal order O_A.
function OO_mod_F_mod_O_mod_F(O::AlgAssAbsOrd)
  OA = maximal_order(algebra(O))
  F = conductor(O, OA, :left)
  FOA = F*OA
  Q1, toQ1 = quo(OA, FOA)
  H1, H1toQ1 = unit_group(Q1)
  Q2, toQ2 = quo(O, F)
  H2, H2toQ2 = unit_group(Q2)
  H2inH1, _ = sub(H1, [ H1toQ1\(toQ1(OA(elem_in_algebra(toQ2\(H2toQ2(H2[i])), copy = false)))) for i = 1:ngens(H2) ])
  H, toH = quo(H1, H2inH1)
  S, StoH = snf(H)
  local _disc_log
  let toH = toH, StoH = StoH, Q1 = Q1
    function _disc_log(x)
      @assert parent(x) === Q1
      s = StoH\(toH(H1toQ1\x))
      return s.coeff
    end
  end

  StoQ1 = GrpAbFinGenToAbsOrdQuoRingMultMap(S, Q1, [ H1toQ1(toH\(StoH(S[i]))) for i = 1:ngens(S) ], _disc_log)
  return S, StoQ1, toQ1
end
# for _unit_group_non_maximal see NfOrd/PicardGroup.jl

# Given an order O in an étale algebra, determine the
# O^+ = {x \in O | x_v > 0 for v in rlpl}
function unit_group_positive(O::AlgAssAbsOrd, rlpl)
  U, mU = unit_group(O)
  co = components(Field, algebra(O))
  n = sum(length.(rlpl))
  A = abelian_group(fill(2, n))
  t = Vector{elem_type(A)}()
  for j in 1:ngens(U)
    u = elem_in_algebra(mU(gen(U, j)))
    imu = Int[]
    for i in 1:length(co)
      uinK = co[i][2](u)
      for r in rlpl[i]
        @assert number_field(r) === co[i][1]
        if ispositive(uinK, r)
          push!(imu, 0)
        else
          push!(imu, 1)
        end
      end
    end
    push!(t, A(imu))
  end
  h = hom(U, A, t)
  K, KtoU = kernel(h)
  S, StoK = snf(K)
  StoU = StoK * KtoU

  exp = function(s)
    @assert parent(s) === S
    return mU(StoU(s))
  end

  log = function(uu)
    @assert parent(uu) === O
    kk = mU\uu
    fl, k = haspreimage(StoU, kk)
    if !fl
      error("Element not positive at described places")
    end
    return k
  end

  StoO = MapUnitGrp{typeof(O)}()
  StoO.header = MapHeader(S, O, exp, log)
  return S, StoO
end

# Determine the kernel of
# O^\times -> (O/g*oo)^
# O must be maximal order in an étale algebra
# mU is the unit group (map) of O
function _unit_group_one_units(O::AlgAssAbsOrd, g, real_places, mU)
  A = algebra(O)
  U = domain(mU)
  Q, mQ = quo(O, g)
  Quni, mQuni = unit_group(Q)
  c = components(Field, A)
  # At the moment I don't need the "exponent" map, so I will do the obvious
  # thing to deal with the infinite places. If one ever wants the exponent
  # map, one would have to patch together
  # infinite_primes_map(OK, p, m)
  # for the different components of A.
  npl = sum(length.(real_places))
  T = abelian_group(fill(2, npl))
  D, (QuniToD, TtoD) = direct_product(Quni, T, task = :sum)
  @assert domain(QuniToD) === Quni
  @assert domain(TtoD) === T
  imgs = elem_type(D)[]
  for a in gens(U)
    u = mU(a)
    img_inf = Int[]
    for i in 1:length(c)
      K, KtoC = c[i]
      for r in real_places[i]
        if sign(KtoC(elem_in_algebra(u)), r) == 1
          push!(img_inf, 0)
        else
          push!(img_inf, 1)
        end
      end
    end
    u_img_inf = T(img_inf)
    push!(imgs, QuniToD((mQuni\(mQ(u)))) + TtoD(u_img_inf))
  end
  h = hom(U, D, imgs)
  K, KtoU = kernel(h)
  S, StoK = snf(K)
  return S, (StoK * KtoU)
end


text:
## mateusz recommendation Manifolds.jl ProductRepr, see AMP.jl#41

using StaticArrays, Manifolds, NearestNeighbors, Distances

M = SpecialEuclidean(2)
N = 100
# convert point to coordinates
function coords(p)
    return SA[p.parts[1][1], p.parts[1][2], acos(p.parts[2][1,1])]
end
# reverse of `coords`
function uncoords(p)
    α = p[3]
    return ProductRepr((SA[p[1], p[2]]), SA[cos(α) -sin(α); sin(α) cos(α)])
end
# some random points to make a tree from
pts = [uncoords(@SVector randn(3)) for _ in 1:N]

# The variant in ManifoldML doesn't support `ProductRepr` currently.
struct SE2Distance{TM<:Manifold} <: Distances.Metric
    manifold::TM
end
function (dist::SE2Distance)(a, b)
    return distance(dist.manifold, uncoords(a), uncoords(b))
end

dist = SE2Distance(M)
# making a tree
vector_elem = coords.(pts)
balltree = BallTree(vector_elem, dist)

# point_matrix = reduce(hcat, map(a -> coords(a), pts))
# balltree = BallTree(point_matrix, dist)
# finding nearest neighbors
k = 3
idxs, dists = knn(balltree, coords(pts[2]), k)

## try make ManifoldFakeArray

mutable struct ManifoldFakeArray{T,N} <: Base.AbstractArray{T,N}
    element::T
end

ManifoldFakeArray{N}(el::T) where {T,N} = ManifoldFakeArray{T,N}(el)

struct ManifoldFakeArrayDistance{TM<:Manifold} <: Distances.Metric
    manifold::TM
end

function (dist::ManifoldFakeArrayDistance{TM})(a, b) where TM
    return distance(dist.manifold, a.element, b.element)
end

import Base: size
Base.size(el::Union{<:ManifoldFakeArray{T,N},Type{<:ManifoldFakeArray{T,N}}}) where {T,N} = N
Base.length(el::Union{<:ManifoldFakeArray{T,N},Type{<:ManifoldFakeArray{T,N}}}) where {T,N} = N

# Base.tail(el::ManifoldFakeArray{T,N}) where {T,N} = T[]

function Base.show(io, mime, mdfa::ManifoldFakeArray{T,N}) where {T,N}
    println(io, "ManifoldFakeArray{T,N}")
    println(io, "  N: $N")
    # println(io, "  T: $T")
    # println(io, "    .element::T: $(mdfa.element)")
end

function Base.show(io, mime, mdfa::AbstractArray{ <: ManifoldFakeArray{T,N}}) where {T,N}
    for i in 1:length(mdfa)
        show(io,mime, mdfa[i])
    end
end

Base.Multimedia.display(mdfa::ManifoldFakeArray{T,N}) where {T,N} = display(T)
Base.Multimedia.display(mdfa::AbstractArray{<:ManifoldFakeArray{T,N}}) where {T,N} = display.(mdfa)

faPts = map(x->ManifoldFakeArray{3}(x), pts)

distFA = ManifoldFakeArrayDistance(M)

distFA(faPts[1], faPts[2])

## real test

balltree = BallTree(faPts, distFA)

## build BallTree for KDE density estimation

using NearestNeighbors: BallTree, RMSDeviation
using Colors
using PyPlot
using PyCall
@pyimport matplotlib.patches as patch

using StaticArrays
using LinearAlgebra

using BenchmarkTools

##

allpts = [SVector(randn(2)...) for i in 1:100];

tree = BallTree(allpts, RMSDeviation(); leafsize = 50)

tree = BallTree(allpts, RMSDeviation(); leafsize = 25)

tree = BallTree(allpts, RMSDeviation(); leafsize = 13)

tree = BallTree(allpts, RMSDeviation(); leafsize = 7)

tree = BallTree(rand(2,100), RMSDeviation(); leafsize = 1)

## Make own distance

# using TransformUtils
using Manifolds

# carfully only include a single definition for use with NearestNeighbors.jl
using Distances: Metric

struct ManiRMSDeviation{M <: ManifoldsBase.AbstractManifold} <: Metric where M
  manifold::M
end

##

mr = ManiRMSDeviation(Euclidean(2))

##

distance(mr.manifold, 1,2)

##

distance(mr.manifold, SVector(1,2),SVector(2,3))

##

M_se2 = SpecialEuclidean(2)

sI = SMatrix{3,3,Float64}(diagm(ones(3)))

T1 = Manifolds.hat(M_se2, sI, SVector(0,0,pi/2))

T2 = Manifolds.hat(M_se2, sI, SVector(1,0,0))

distance(M_se2, T1, T2)

##

md = ManiRMSDeviation{}()

rd = RMSDeviation()

isconcretetype(rd)

# (dist::MeanSqDeviation)(a, b) = sqeuclidean(a, b) / length(a)
# msd(a, b) = MeanSqDeviation()(a, b)

# (dist::RMSDeviation)(a, b) = sqrt(MeanSqDeviation()(a, b))
# rmsd(a, b) = RMSDeviation()(a, b)

@code_warntype rmsd(1,1)

sqeuclidean(a, b) / length(a)

(::ManiRMSDeviation{T})(a, b) where T = sqrt(MeanSqDeviation()(a, b))

##

import NearestNeighbors.HyperSphere

##

# Adds a sphere to an axis
function add_sphere(ax, hs::HyperSphere, col)
    ell = patch.Circle(hs.center, radius = hs.r, facecolor="none", edgecolor=col)
    ax.add_artist(ell)
end

# Skip non leaf nodes
offset = tree.tree_data.n_internal_nodes + 1
nleafs = tree.tree_data.n_leafs

# Range of leaf nodes
index_range = offset: offset + nleafs - 1

# Generate some nice colors
cols = distinguishable_colors(length(index_range), RGB(0,0,0))

# Create figure
cfig = figure()
ax = cfig.add_subplot(1,1,1)
ax.set_aspect("equal")
axis((-.25,1.25,-.25,1.25))

for (i, idx) = enumerate(index_range)
  col = cols[i]
  # Get the indices of the leaf nodes into the tree data
  range = NearestNeighbors.get_leaf_range(tree.tree_data, idx)
  d = tree.data[range]

  # Plot the points in the hyper spehre
  plot(getindex.(d, 1), getindex.(d, 2), "*", color = (col.r, col.g, col.b))

  # And the hypersphere itself
  sphere = tree.hyper_spheres[idx]
  add_sphere(ax, sphere, (col.r, col.g, col.b))
end

title("Leaf nodes with their corresponding points")

##

using Distances

MeanSqDeviation()([1;2;3;4], [2;3;4;5])

#


text:
using Test, Agents, Random, LightGraphs, SQLite, DataFrames

mutable struct Agent0 <: AbstractAgent
  id::Int
end

mutable struct Agent1 <: AbstractAgent
  id::Int
  pos::Tuple{Int,Int}
end

mutable struct Agent2 <: AbstractAgent
  id::Int
  weight::Float64
end

mutable struct Agent3 <: AbstractAgent
  id::Int
  pos::Tuple{Int,Int}
  weight::Float64
end

mutable struct Agent4 <: AbstractAgent
  id::Int
  pos::Tuple{Int,Int}
  p::Int
end

mutable struct Agent5 <: AbstractAgent
  id::Int
  pos::Int
  weight::Float64
end

mutable struct Agent6 <: AbstractAgent
  id::Int
  pos::NTuple{2, Float64}
  vel::NTuple{2, Float64}
  weight::Float64
end

#
mutable struct Agent7 <: AbstractAgent
  id::Int
  pos::Int
  f1::Bool
  f2::Int
end

Agent7(id, pos; f1, f2) = Agent7(id, pos, f1, f2)

@testset "Agents.jl Tests" begin

include("api_tests.jl")
include("model_access.jl")
include("space_test.jl")
include("interaction_tests.jl")
include("data_collector_test.jl")
include("CA_test.jl")
include("continuousSpace_tests.jl")
include("collisions_tests.jl")

end


text:
# Modified from Graphs.jl.
# Workaround for edge_index not being O(1) on SimpleGraph.
# edge_index was only needed to test for cycles, so
# this implementation skips that check.

# Graphs.jl is licensed under the MIT License:
#
# Copyright (c) 2012: John Myles White and other contributors.
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

mutable struct TopologicalSortVisitor
    vertices::Vector{Int}
    parents::Vector{Int}

    function TopologicalSortVisitor(n::Int)
        vs = Int[]
        sizehint!(vs, n)
        return new(vs, zeros(Int, n))
    end
end

function depth_first_visit_impl!(
    adjlist,
    offsets,
    vertex_stack,                   # an (initialized) stack of vertex
    index_stack,                    # stack with out edge indices
    vertexcolormap::Vector{Int},    # an (initialized) color-map to indicate status of vertices
    visitor::TopologicalSortVisitor,
)
    while !isempty(vertex_stack)
        u = pop!(vertex_stack)
        out_idx = pop!(index_stack)
        #uegs = out_edges(u,graph)
        len_uegs = offsets[u+1] - offsets[u]
        found_new_vertex = false

        while out_idx <= len_uegs && !found_new_vertex
            v = adjlist[offsets[u]+out_idx-1]
            out_idx += 1
            v_color = vertexcolormap[v]

            if v_color == 0
                found_new_vertex = true
                vertexcolormap[v] = 1
                push!(vertex_stack, u)
                push!(index_stack, out_idx)
                visitor.parents[v] = u

                push!(vertex_stack, v)
                push!(index_stack, 1)
            end
        end

        if !found_new_vertex
            close_vertex!(visitor, u)
            vertexcolormap[u] = 2
        end
    end
end

function traverse_graph(
    adjlist,
    offsets,
    s,
    visitor::TopologicalSortVisitor,
    vertexcolormap,
    vertex_stack,
    index_stack,
)
    vertexcolormap[s] = 1

    resize!(vertex_stack, 1)
    vertex_stack[1] = s
    resize!(index_stack, 1)
    index_stack[1] = 1

    return depth_first_visit_impl!(
        adjlist,
        offsets,
        vertex_stack,
        index_stack,
        vertexcolormap,
        visitor,
    )
end

function close_vertex!(visitor::TopologicalSortVisitor, v::Int)
    return push!(visitor.vertices, v)
end

function reverse_topological_sort_by_dfs(
    adjlist,
    offsets,
    num_vertices,
    cmap::Vector{Int},
    vertex_stack::Vector{Int} = Int[],
    index_stack::Vector{Int} = Int[],
)
    @assert length(cmap) == num_vertices
    fill!(cmap, 0)
    visitor = TopologicalSortVisitor(num_vertices)

    for s in 1:num_vertices
        if cmap[s] == 0
            traverse_graph(
                adjlist,
                offsets,
                s,
                visitor,
                cmap,
                vertex_stack,
                index_stack,
            )
        end
    end

    return visitor.vertices, visitor.parents
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SteamTurbineGov1 <: TurbineGov
        R::Float64
        T1::Float64
        valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        T2::Float64
        T3::Float64
        D_T::Float64
        DB_h::Float64
        DB_l::Float64
        T_rate::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes.StateType}
        internal::InfrastructureSystemsInternal
    end

Steam Turbine-Governor. This model considers both TGOV1 or TGOV1DU in PSS/E.

# Arguments
- `R::Float64`: Droop parameter, validation range: `(0, 0.1)`, action if invalid: `warn`
- `T1::Float64`: Governor time constant, validation range: `(eps(), 0.5)`, action if invalid: `error`
- `valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Valve position limits
- `T2::Float64`: Lead Lag Lead Time constant , validation range: `(0, nothing)`, action if invalid: `warn`
- `T3::Float64`: Lead Lag Lag Time constant , validation range: `(eps(), 10)`, action if invalid: `error`
- `D_T::Float64`: Turbine Damping, validation range: `(0, 0.5)`, action if invalid: `warn`
- `DB_h::Float64`: Deadband for overspeed, validation range: `(0, nothing)`, action if invalid: `warn`
- `DB_l::Float64`: Deadband for underspeed, validation range: `(nothing, 0)`, action if invalid: `warn`
- `T_rate::Float64`: Turbine Rate (MW). If zero, generator base is used., validation range: `(0, nothing)`, action if invalid: `warn`
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the SteamTurbineGov1 model are:
	x_g1: Valve Opening,
	x_g2: Lead-lag state
- `n_states::Int`: TGOV1 has 2 states
- `states_types::Vector{StateTypes.StateType}`: TGOV1 has 2 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SteamTurbineGov1 <: TurbineGov
    "Droop parameter"
    R::Float64
    "Governor time constant"
    T1::Float64
    "Valve position limits"
    valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Lead Lag Lead Time constant "
    T2::Float64
    "Lead Lag Lag Time constant "
    T3::Float64
    "Turbine Damping"
    D_T::Float64
    "Deadband for overspeed"
    DB_h::Float64
    "Deadband for underspeed"
    DB_l::Float64
    "Turbine Rate (MW). If zero, generator base is used."
    T_rate::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the SteamTurbineGov1 model are:
	x_g1: Valve Opening,
	x_g2: Lead-lag state"
    states::Vector{Symbol}
    "TGOV1 has 2 states"
    n_states::Int
    "TGOV1 has 2 differential states"
    states_types::Vector{StateTypes.StateType}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SteamTurbineGov1(R, T1, valve_position_limits, T2, T3, D_T, DB_h, DB_l, T_rate, P_ref=1.0, ext=Dict{String, Any}(), )
    SteamTurbineGov1(R, T1, valve_position_limits, T2, T3, D_T, DB_h, DB_l, T_rate, P_ref, ext, [:x_g1, :x_g2], 2, [StateTypes.Differential, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function SteamTurbineGov1(; R, T1, valve_position_limits, T2, T3, D_T, DB_h, DB_l, T_rate, P_ref=1.0, ext=Dict{String, Any}(), states=[:x_g1, :x_g2], n_states=2, states_types=[StateTypes.Differential, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    SteamTurbineGov1(R, T1, valve_position_limits, T2, T3, D_T, DB_h, DB_l, T_rate, P_ref, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function SteamTurbineGov1(::Nothing)
    SteamTurbineGov1(;
        R=0,
        T1=0,
        valve_position_limits=(min=0.0, max=0.0),
        T2=0,
        T3=0,
        D_T=0,
        DB_h=0,
        DB_l=0,
        T_rate=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`SteamTurbineGov1`](@ref) `R`."""
get_R(value::SteamTurbineGov1) = value.R
"""Get [`SteamTurbineGov1`](@ref) `T1`."""
get_T1(value::SteamTurbineGov1) = value.T1
"""Get [`SteamTurbineGov1`](@ref) `valve_position_limits`."""
get_valve_position_limits(value::SteamTurbineGov1) = value.valve_position_limits
"""Get [`SteamTurbineGov1`](@ref) `T2`."""
get_T2(value::SteamTurbineGov1) = value.T2
"""Get [`SteamTurbineGov1`](@ref) `T3`."""
get_T3(value::SteamTurbineGov1) = value.T3
"""Get [`SteamTurbineGov1`](@ref) `D_T`."""
get_D_T(value::SteamTurbineGov1) = value.D_T
"""Get [`SteamTurbineGov1`](@ref) `DB_h`."""
get_DB_h(value::SteamTurbineGov1) = value.DB_h
"""Get [`SteamTurbineGov1`](@ref) `DB_l`."""
get_DB_l(value::SteamTurbineGov1) = value.DB_l
"""Get [`SteamTurbineGov1`](@ref) `T_rate`."""
get_T_rate(value::SteamTurbineGov1) = value.T_rate
"""Get [`SteamTurbineGov1`](@ref) `P_ref`."""
get_P_ref(value::SteamTurbineGov1) = value.P_ref
"""Get [`SteamTurbineGov1`](@ref) `ext`."""
get_ext(value::SteamTurbineGov1) = value.ext
"""Get [`SteamTurbineGov1`](@ref) `states`."""
get_states(value::SteamTurbineGov1) = value.states
"""Get [`SteamTurbineGov1`](@ref) `n_states`."""
get_n_states(value::SteamTurbineGov1) = value.n_states
"""Get [`SteamTurbineGov1`](@ref) `states_types`."""
get_states_types(value::SteamTurbineGov1) = value.states_types
"""Get [`SteamTurbineGov1`](@ref) `internal`."""
get_internal(value::SteamTurbineGov1) = value.internal

"""Set [`SteamTurbineGov1`](@ref) `R`."""
set_R!(value::SteamTurbineGov1, val) = value.R = val
"""Set [`SteamTurbineGov1`](@ref) `T1`."""
set_T1!(value::SteamTurbineGov1, val) = value.T1 = val
"""Set [`SteamTurbineGov1`](@ref) `valve_position_limits`."""
set_valve_position_limits!(value::SteamTurbineGov1, val) = value.valve_position_limits = val
"""Set [`SteamTurbineGov1`](@ref) `T2`."""
set_T2!(value::SteamTurbineGov1, val) = value.T2 = val
"""Set [`SteamTurbineGov1`](@ref) `T3`."""
set_T3!(value::SteamTurbineGov1, val) = value.T3 = val
"""Set [`SteamTurbineGov1`](@ref) `D_T`."""
set_D_T!(value::SteamTurbineGov1, val) = value.D_T = val
"""Set [`SteamTurbineGov1`](@ref) `DB_h`."""
set_DB_h!(value::SteamTurbineGov1, val) = value.DB_h = val
"""Set [`SteamTurbineGov1`](@ref) `DB_l`."""
set_DB_l!(value::SteamTurbineGov1, val) = value.DB_l = val
"""Set [`SteamTurbineGov1`](@ref) `T_rate`."""
set_T_rate!(value::SteamTurbineGov1, val) = value.T_rate = val
"""Set [`SteamTurbineGov1`](@ref) `P_ref`."""
set_P_ref!(value::SteamTurbineGov1, val) = value.P_ref = val
"""Set [`SteamTurbineGov1`](@ref) `ext`."""
set_ext!(value::SteamTurbineGov1, val) = value.ext = val
"""Set [`SteamTurbineGov1`](@ref) `states_types`."""
set_states_types!(value::SteamTurbineGov1, val) = value.states_types = val


text: if "foo" end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using ACE: State, DState, ForwardDiff
using StaticArrays
# using NamedTupleTools
using JuLIP: AtomicNumber

function housholderreflection(rr0::SVector{3})
   r0 = norm(rr0)
   I3x3 = SMatrix{3,3}(1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
   if r0 == 0
      return I3x3
   end
   v = SVector(0, 0, 1) - rr0/r0
   v̂ = v / (norm(v) + eps(eltype(rr0))^2)
   return I3x3 - 2 * v̂ * v̂'
end

function pullback_housholderreflection(rr0::SVector{3})
   H = housholderreflection(rr0)
   dH = ForwardDiff.jacobian(rr -> housholderreflection(rr)[:], rr0)
   dHt = SMatrix{3,9}(dH)'
   return H, g -> SMatrix{3,3}(dHt * SVector{3}(g))
end

_xyz2rθz(ss) = SVector(sqrt(ss[1]^2 + ss[2]^2),
                       atan(ss[2], ss[1]),
                       ss[3])

# function _rθz2xyz_d(ss)
#    s, θ, z = _xyz2rθz(ss)
#    sθ, cθ = sincos(θ)
#    return SA[ cθ -s*sθ 0;
#               sθ  s*cθ 0;
#                0    0  1 ]
# end

@doc raw"""
This implements the jacobi matrix
```math
   J = \frac{\partial (r, \theta, z)}{\partial (x, y, z)}
```
"""
function _xyz2rθz_d(ss)
   s, θ, z = _xyz2rθz(ss)
   sθ, cθ = sincos(θ)
   return SA[ cθ    sθ 0;
              -sθ/s cθ/s 0;
               0    0 1 ]
end

# this is just a helper function which should not be used for
# performance critical code
function eucl2cyl(rr0::SVector, rr::SVector)
   H = housholderreflection(rr0)
   ss = H * rr
   r, θ, z = _xyz2rθz(ss)
   return r, θ, z
end

"""
transforms euclidean to cylindral bond environment

### Input
* `rr0` : central bond
* `Rs` : list of environment atom positions, relative to rr0/2
* `Zs` : associated list of atomic species

### Output
* An ACE configuration where each state represents either the bond in terms
of just its length; or an atom from the environment in terms of its
cylindrical coordinates (r, θ, z) and species mu, with origin at `rr0/2`.
"""
function eucl2cyl(rrij::SVector, Zi, Zj,
                  Rs::AbstractVector{<: SVector},
                  Zs::AbstractVector{<: AtomicNumber})
   @assert length(Rs) == length(Zs)
   H = housholderreflection(rrij)
   rij = norm(rrij)

   function _eucl2cyl(rr, mu)
      ss = H * rr
      r, θ, z = _xyz2rθz(ss)
      return State( mu = mu,
                    r = r, θ = θ, z = z,
                    rij = rij,
                    mui = Zi,
                    muj = Zj,
                    be = :env )
   end

   Y0 = State( mu = AtomicNumber(0), r = 0.0, θ = 0.0, z = 0.0,
               rij = rij, mui = Zi, muj = Zj, be = :bond )
   cfg = Vector{typeof(Y0)}(undef, length(Rs)+1)
   cfg[1] = Y0
   for i = 1:length(Rs)
      cfg[i+1] = _eucl2cyl(Rs[i], Zs[i])
   end

   return cfg
end

"""
ssj = H * rrj
cj = xyz2rθz(ssj)
∂_rrj(f) = ∂_cj(f) * ∂_ssj(cj) * ∂_rrj(ssj)
∇_rrj(f) = ∂_rrj(ssj)' * ∂_ssj(cj)' * ∇_cj(f)
∇_rrj(f) = ∂_rr0(ssj)' * ∂_ssj(cj)' * ∇_cj(f)
"""
function rrule_eucl2cyl(rr0::SVector, Zi, Zj,
                        Rs::AbstractVector{<: SVector},
                        Zs::AbstractVector{<: AtomicNumber},
                        g_cyl::AbstractMatrix{<: DState})
   lenB = size(g_cyl, 1)
   lenR = length(Rs)
   @assert size(g_cyl, 2) == lenR + 1
   H = housholderreflection(rr0)
   H, pbH = pullback_housholderreflection(rr0)
   r̂0 = rr0 / norm(rr0) # ∇f(r0) = f'(r0) * r̂0

   g_Rs = zeros(SVector{3, Float64}, lenB, lenR)
   # g_cyl[:, 1] = derivative w.r.t. rr0 only
   g_rr0 = [ g_cyl[n, 1].rij * r̂0 for n = 1:lenB ]

   for j = 1:lenR
      rrj = Rs[j]
      ss = H * rrj
      s, θ, z = _xyz2rθz(ss)
      J = _xyz2rθz_d(ss)
      for n = 1:lenB
         gj = g_cyl[n, j+1]
         gj1 = J' * SVector(gj.r, gj.θ, gj.z)
         g_rr0[n] += gj.rij * r̂0 + pbH(gj1)' * rrj
         g_Rs[n, j] = H' * gj1
      end
   end

   return g_rr0, g_Rs
end

function rrule_eucl2cyl(rr0::SVector, Zi, Zj,
                        Rs::AbstractVector{<: SVector},
                        Zs::AbstractVector{<: AtomicNumber},
                        g_cyl::AbstractVector{<: DState})
   lenR = length(Rs)
   @assert length(g_cyl) == lenR + 1
   H = housholderreflection(rr0)
   H, pbH = pullback_housholderreflection(rr0)
   r̂0 = rr0 / norm(rr0) # ∇f(r0) = f'(r0) * r̂0

   g_Rs = zeros(SVector{3, Float64}, lenR)
   # deriv. of first element w.r.t. rr0 only
   g_rr0 = g_cyl[1].rij * r̂0

   for j = 1:lenR
      rrj = Rs[j]
      ss = H * rrj
      s, θ, z = _xyz2rθz(ss)
      J = _xyz2rθz_d(ss)
      gj = g_cyl[j+1]
      gj1 = J' * SVector(gj.r, gj.θ, gj.z)
      g_rr0 += gj.rij * r̂0 + pbH(gj1)' * rrj
      g_Rs[j] = H' * gj1
   end

   return g_rr0, g_Rs
end

# monkey-patch JuLIP since we aren't really updating it anymore
Base.isapprox(a::AtomicNumber, b::AtomicNumber) = (a == b)

"""
For testing only. Generates a cylindrical bond environment, and returns it
both in euclidean and cylindrical coordinates. The origin for both is the
mid-point of the bond.
"""
function rand_env(r0cut, rcut, zcut; Nenv = 10, species = [:Al, :Ti])
   species = AtomicNumber.(species)
   r0 = 1 + rand()
   rr0 = randn(SVector{3, Float64})
   rr0 = r0 * (rr0/norm(rr0))
   H = ACEbonds.housholderreflection(rr0)
   Rs = SVector{3, Float64}[]
   Zs = AtomicNumber[]
   Xcyl = []
   Zi, Zj = rand(species), rand(species)
   for _ = 1:Nenv
      z = (rand() - 0.5) * 2 * (zcut + r0cut)
      r = rand() * rcut
      θ = rand() * 2 * π - π
      sθ, cθ = sincos(θ)
      rr = H' * SVector(cθ * r, sθ * r, z)
      mu = rand(species)
      push!(Zs, mu)
      push!(Rs, rr)
      push!(Xcyl, State(mu = mu, r = r, θ=θ, z = z, rij = r0, mui = Zi, muj = Zj, be = :env))
   end
   return rr0, Zi, Zj, Rs, Zs, Xcyl
end


text:
#Pkg.update()
#Pkg.add("CSV")
#Pkg.add("Distributions")
#Pkg.add("JSON")
#Pkg.add("FileIO")
#Pkg.add("DataFrames")
using CSV, Distributions, JSON, FileIO, DataFrames

### Import source files
include("src/processgeodata.jl");
include("src/createdatamatrix.jl");
include("src/createpdrive.jl");
include("src/createpdestin.jl");
include("src/initializestates.jl");
include("src/solveinitialvalueproblem.jl");
include("src/resampling.jl");
include("src/averagedrivingtime.jl");
include("src/correctparameters.jl");
include("src/saveresults.jl");
include("src/saveparameters.jl");
include("src/createresultsdirectory.jl");

### Provide the path to the Uber datasets and the path to the folder where you want to save your results
path_to_cities =  # "***Insert path to Uber Movement travel time data ***" ;
path_to_results_folder = # **Insert path to directory where sampling results should be saved **;
city_list = readdir(path_to_cities); #getting the names of directories

############### Insert a list of particular cities that you want to sample here ##############
#city_list = [" ***Insert a city here*** "]
##############################################################################################

############################### Define these parameters manually #############################
e_drive = 0.5; # An exponent that determines the dynamics of driving activity. Can be changed.
e_dest = 2; # An exponent that determines the dynamics of destination popularity. Can be changed.
p_min = 0.1; # A parameter for p_drive. Can be changed.
p_max = 0.9; # A parameter for p_drive. Can be changed.
cars_per_zone = 1000; # the number of cars to be sampled per each city zone
T = 24; # The datasets we analyze have an hourly time resolution
##############################################################################################

for city_index = 1:length(city_list)

    ### Set the path to city files that we want to simulate
    city = city_list[city_index];
    path_to_data = string(path_to_cities, city);
    dataset_list = readdir(path_to_data);
    path_to_json_data = string(path_to_data,"/",dataset_list[length(dataset_list)])
    csv_dataset_list = dataset_list[1:(length(dataset_list)-1)];

    ### create a directory for saving the results of the city
    path_to_results = createresultsdirectory(path_to_results_folder, city)

    ### process the GEOJSON data: INPUT is julia_file and OUTPUT is distance_matrix_km , number_zones
    distance_matrix_km, number_zones = processgeodata(path_to_json_data, path_to_data, csv_dataset_list, path_to_results);

    ### determine the number of cars that are simulated for the current city.
    C = number_zones * cars_per_zone;
    C = trunc(Int, C); # transforms C from type Float64 to Integer. Necessary for subsequent matrix definitions.
    println("Simulation starts for the city of ", city, " with a vehicle fleet of C = ", C);

    ### initialize values for counting
    A_drive = 0;
    count_dataset = 0;

    ### Starts the simulation for the data sets of the current city
    for dataset_index = 1:length(csv_dataset_list)

        ### Set the path to dataset that we want to simulate
        data_set = csv_dataset_list[dataset_index];
        path_to_csv_data = string(path_to_data,"/",data_set);
        println("Simulating ", data_set)

        ### Processing dataset: INPUT is path_to_csv_data, number_zones and OUTPUT is datamatrix
        datamatrix = createdatamatrix(path_to_csv_data, number_zones);

        ### Creating driving activity: INPUT is datamatrix, p_min, p_max, number_zones and OUTPUT is p_drive
        p_drive = createpdrive(datamatrix, distance_matrix_km, number_zones);

        ### Creating destination popularity: INPUT is datamatrix,  number_zones and OUTPUT is p_dest
        p_dest = createpdestin(datamatrix, number_zones);

        ### Initializing state and state transitions: INPUT is C and OUTPUT state_matrix, transition_matrix
        state_matrix, transition_matrix = initializestates(C);

        ### Initial Value Problem: INPUT is state_matrix, p_drive, p_dest, C and OUTPUT initial_state
        initial_state = solveinitialvalueproblem(state_matrix, transition_matrix, p_drive, p_dest, C, number_zones)
        state_matrix[:,1] = initial_state;

        ### Resampling: INPUT state_matrix, transition_matrix, C, number_zones, p_drive, p_dest, datamatrix, distance_matrix_km and OUTPUT state_matrix, transition_matrix
        state_matrix, transition_matrix = resampling(state_matrix, transition_matrix, C, number_zones, p_drive, p_dest, datamatrix, distance_matrix_km);

        ### Calculating average driving time: INPUT C, A_drive[N], transition_matrix and OUTPUT A_drive
        A_drive = averagedrivingtime(C, A_drive, transition_matrix)
        count_dataset += 1; # increment the number of datasets contained in this average

        ### saving results
        saveresults(number_zones, state_matrix, transition_matrix, path_to_results, data_set, C);

    end

    ### Average A_drive for all sampled datasets of the current city
    A_drive  = A_drive/count_dataset;

    ### save parameters
    saveparameters(path_to_results, T, number_zones, cars_per_zone, C, e_drive, p_min, p_max, e_dest, A_drive);

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps