text:
# This file holds all the functions related to processing simulation results.

# The prerequisite group that allows counting of personnel members.
# countPersonnelPrereqGroup = PrerequisiteGroup()
# addPrereq!( countPersonnelPrereqGroup, Prerequisite( :status, :active,
#     relation = ==, valType = Symbol ) )

export getActiveAtTime,
       getActiveAgesAtTime,
       getInFlux,
       getOutFlux,
       getLastChanges,
       getDatabaseAtTime

"""
```
getActiveAtTime( mpSim::ManpowerSimulation,
                 timePoint::T,
                 fields::Vector{String} = Vector{String}() )
    where T <: Real
```
This function retrieves a view of all the personnel in the simulation `mpSim`
who were active at time `timePoint`, and displays the ID key and all the columns
listed in `fields`.

If the given time is < 0 or larger than the current simulation time or the
length of the simulation (whichever is smaller), this function returns
`nothing`. Otherwise, it returns a `DataFrames.DataFrame` object.
"""
function getActiveAtTime( mpSim::ManpowerSimulation, timePoint::T,
    fields::Vector{String} = Vector{String}() ) where T <: Real

    if ( timePoint < 0 ) || ( timePoint > min( now( mpSim ), mpSim.simLength ) )
        return
    end  # if ( timePoint < 0 ) || ...

    tmpFields = unique( vcat( mpSim.idKey, fields ) )
    queryCmd = "SELECT '$(join( tmpFields, "', '" ))'
        FROM $(mpSim.personnelDBname)
        WHERE ( timeEntered <= $timePoint ) AND
            ( ( timeExited > $timePoint ) OR ( timeExited IS NULL ) )
        ORDER BY $(mpSim.idKey)"
    return SQLite.query( mpSim.simDB, queryCmd )

end  # getActiveAtTime( mpSim, timePoint, fields )

"""
```
getActiveAtTime( mpSim::ManpowerSimulation,
                 stateName::String,
                 timePoint::T,
                 fields::Vector{String} = Vector{String}() )
    where T <: Real
```
This function retrieves a view of all the personnel in the simulation `mpSim`
who had the state with name `stateName` at time `timePoint`, and displays the ID
key and all the columns listed in `fields`.

If the given time is < 0 or larger than the current simulation time or the
length of the simulation (whichever is smaller), this function returns
`nothing`. Otherwise, it returns a `DataFrames.DataFrame` object.
"""
function getActiveAtTime( mpSim::ManpowerSimulation, stateName::String,
    timePoint::T, fields::Vector{String} = Vector{String}() ) where T <: Real

    if ( timePoint < 0 ) || ( timePoint > min( now( mpSim ), mpSim.simLength ) )
        return
    end  # if ( timePoint < 0 ) || ...

    tmpFields = unique( vcat( mpSim.idKey, fields ) )
    queryCmd = "SELECT '$(join( tmpFields, "', '" ))'
        FROM $(mpSim.personnelDBname)
        WHERE timeEntered < 0"
    dummyResult = SQLite.query( mpSim.simDB, queryCmd )

    # Get the IDs of all active personnel members at that time.
    activeIDs = getActiveAtTime( mpSim, timePoint )[ 1 ]

    # Don't continue if there are no active personnel at that time.
    if isempty( activeIDs )
        return dummyResult
    end  # if isempty( activeIDs )

    # In that list, find all IDs of personnel that are in the required state at
    #   the given time.
    queryCmd = "SELECT $(mpSim.idKey), count( $(mpSim.idKey) ) appearances
	    FROM $(mpSim.transitionDBname)
	    WHERE ( $(mpSim.idKey) IN ('$(join( activeIDs, "', '" ))') ) AND
            ( timeIndex <= $timePoint ) AND
            ( ( endState IS '$stateName' ) OR ( startState IS '$stateName' ) )
	    GROUP BY $(mpSim.idKey)"
    activeIDs = SQLite.query( mpSim.simDB, queryCmd )
    activeIDs = activeIDs[ Symbol( mpSim.idKey ) ][
        activeIDs[ :appearances ] .== 1 ]

    # Don't continue if there are no active personnel with the given state at
    #   that time.
    if isempty( activeIDs )
        return dummyResult
    end  # if isempty( activeIDs )

    # Get the requested fields of the active IDs.
    queryCmd = "SELECT '$(join( tmpFields, "', '" ))'
        FROM $(mpSim.personnelDBname)
        WHERE $(mpSim.idKey) IN ('$(join( activeIDs, "', '" ))')"
    return SQLite.query( mpSim.simDB, queryCmd )

end  # getActiveAtTime( mpSim, stateName, timePoint, fields )

"""
```
getActiveAgesAtTime( mpSim::ManpowerSimulation,
                     timePoint::T )
    where T <: Real
```
This function returns the ages of all personnel members in the manpower
simulation `mpSim` who were active at time `timePoint`.

If the given time is < 0 or larger than the current simulation time or the
length of the simulation (whichever is smaller), the function returns `nothing`.
Otherwise, it returns a `Vector{float64}` with the ages.
"""
function getActiveAgesAtTime( mpSim::ManpowerSimulation, timePoint::T ) where T <: Real

    if ( timePoint < 0 ) || ( timePoint > min( now( mpSim ), mpSim.simLength ) )
        return
    end  # if ( timePoint < 0 ) || ...

    # Get all necessary information from the database.
    activePersonnel = getActiveAtTime( mpSim, timePoint,
        [ "ageAtRecruitment + $timePoint - timeEntered age" ] )

    # Compute the current age of all active personnel members.
    return activePersonnel[ :age ]

end  # getActiveAgesAtTime( mpSim, timePoint )

"""
```
getInFlux( mpSim::ManpowerSimulation,
           tBegin::T1,
           tEnd::T2 )
    where T1 <: Real
    where T2 <: Real
```
This function retrieves a view of all the personnel members in the simulation
`mpSim` who entered the organisation between `tBegin` (exclusive) and `tEnd`
(inclusive). People who entered and left the organisation in this interval are
not counted.

If the end time is <= 0, or the start time is larger than the current simulation
time or the length of the simulation (whichever is smaller), or if the end time
is smaller than the start time, this function returns `nothing`.
"""
function getInFlux( mpSim::ManpowerSimulation, tBegin::T1, tEnd::T2 ) where T1 <: Real where T2 <: Real

    if ( tEnd <= 0 ) || ( tBegin > min( now( mpSim ), mpSim.simLength ) ) ||
        ( tEnd <= tBegin )
        return
    end  # if ( tEnd <= 0 ) || ...

    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.personnelDBname)
        WHERE ( $tBegin < timeEntered ) AND ( timeEntered <= $tEnd ) AND
            ( ( timeExited > $tEnd ) OR ( timeExited IS NULL ) )"
    SQLite.query( mpSim.simDB, queryCmd )

end  # getInFlux( mpSim, tBegin, tEnd, fields )

"""
```
getInFlux( mpSim::ManpowerSimulation,
           stateName::String,
           tBegin::T1,
           tEnd::T2 )
    where T1 <: Real
    where T2 <: Real
```
This function retrieves a view of all the personnel members in the simulation
`mpSim` who entered the state with name `stateName` between `tBegin` (exclusive)
and `tEnd` (inclusive). People who entered and left the state (or the
organisation) in this interval are not counted. The view also includes the
source (start state) and the reason (transition) of the flux for each personnel
member.

If the end time is <= 0, or the start time is larger than the current simulation
time or the length of the simulation (whichever is smaller), or if the end time
is smaller than the start time, this function returns `nothing`.
"""
function getInFlux( mpSim::ManpowerSimulation, stateName::String, tBegin::T1,
    tEnd::T2 ) where T1 <: Real where T2 <: Real

    if ( tEnd <= 0 ) || ( tBegin > min( now( mpSim ), mpSim.simLength ) ) ||
        ( tEnd <= tBegin )
        return
    end  # if ( tEnd <= 0 ) || ...

    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.personnelDBname)
        WHERE timeEntered < 0"
    dummyResult = SQLite.query( mpSim.simDB, queryCmd )

    # Retrieve IDs of people leaving the organisation in time interval.
    queryCmd = "SELECT $(mpSim.idKey) FROM $(mpSim.personnelDBname)
        WHERE ( $tBegin < timeExited ) AND ( timeExited <= $tEnd )"
    fluxOutIDs = SQLite.query( mpSim.simDB, queryCmd )[ 1 ]

    # Retrieve IDs and amount of transitions starting or ending in the state.
    queryCmd = "SELECT $(mpSim.idKey), count($(mpSim.idKey)) numEntries
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) NOT IN ('$(join( fluxOutIDs, "', '" ))') ) AND
            ( ( startState IS '$stateName' ) OR ( endState IS '$stateName' ) )
        GROUP BY $(mpSim.idKey)"
    queryRes = SQLite.query( mpSim.simDB, queryCmd )
    nEntries, nCold = size( queryRes )
    numEntries = Dict{String, Vector{Int}}()

    for ii in 1:nEntries
        numEntries[ queryRes[ 1 ][ ii ] ] = [ queryRes[ 2 ][ ii ], 0 ]
    end  # for ii in 1:nEntries

    # Retrieve IDs and amount of transitions ending in the state.
    queryCmd = "SELECT $(mpSim.idKey), count($(mpSim.idKey)) numEntries
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) NOT IN ('$(join( fluxOutIDs, "', '" ))') ) AND
            ( endState IS '$stateName' )
        GROUP BY $(mpSim.idKey)"
    queryRes = SQLite.query( mpSim.simDB, queryCmd )
    nEntries, nCold = size( queryRes )

    for ii in 1:nEntries
        numEntries[ queryRes[ 1 ][ ii ] ][ 2 ] = queryRes[ 2 ][ ii ]
    end  # for ii in 1:nEntries

    # If most transitions are into the state, personnel member is in the state
    #   at the end of the time interval.
    fluxInIDs = collect( keys( numEntries ) )
    filter!( id -> ( numEntries[ id ][ 2 ] * 2 > numEntries[ id ][ 1 ] ),
        fluxInIDs )

    if isempty( fluxInIDs )
        return dummyResult
    end  # if isempty( fluxInIDs )

    # Get start state and transition name for previously identified people.
    queryCmd = "SELECT $(mpSim.idKey), startState, transition
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( endState IS '$stateName' ) AND
            ( $(mpSim.idKey) IN ('$(join( fluxInIDs, "', '" ))') )
        GROUP BY $(mpSim.idKey)"
    result = SQLite.query( mpSim.simDB, queryCmd )

    # If person has multiple transitions in/out of the state, set start state
    #   and transition to unknown.
    for ii in 1:length( fluxInIDs )
        id = result[ Symbol( mpSim.idKey ) ][ ii ]

        if numEntries[ id ][ 1 ] > 1
            result[ :startState ][ ii ] = "unknown"
            result[ :transition ][ ii ] = "unknown"
        end  # if numEntries[ id ][ 1 ] > 1
    end  # for ii in 1:length( fluxInIDs )

    return result

#=
    # Retrieve the list of people achieving the state in the time interval.
    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( endState IS '$stateName' )"
    fluxInIDs = SQLite.query( mpSim.simDB, queryCmd )[ 1 ]

    # If no personnel members entered the state, no need to continue.
    if isempty( fluxInIDs )
        return dummyResult
    end  # if isempty( fluxInIDs )

    # Retrieve the list of those personnel members leaving the state or the
    #   organisation in the time interval.
    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) IN ('$(join( fluxInIDs, "', '" ))') ) AND
            ( ( startState IS '$stateName' ) OR ( endState IS NULL ) )"
    fluxOutIDs = SQLite.query( mpSim.simDB, queryCmd )[ 1 ]
    filter!( id -> id ∉ fluxOutIDs, fluxInIDs )

    queryCmd = "SELECT $(mpSim.idKey), startState, transition
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) IN ('$(join( fluxInIDs, "', '" ))') ) AND
            ( endState IS '$stateName' )"

    return SQLite.query( mpSim.simDB, queryCmd )
=#
end  # getInFlux( mpSim, tBegin, tEnd, fields )

"""
```
getOutFlux( mpSim::ManpowerSimulation,
            tBegin::T1,
            tEnd::T2,
            fields::Vector{String} = Vector{String}() )
    where T1 <: Real
    where T2 <: Real
```
This function retrieves a view of all the personnel members in the simulation
`mpSim` who left the organisation between `tBegin` (inclusive) and `tEnd`
(exclusive). The attributes in `fields` will also be snown in the retrieved
view. People who entered and left the organisation in this interval are not
counted.

If the end time is <= 0, or the start time is larger than the current simulation
time or the length of the simulation (whichever is smaller), or if the end time
is smaller than the start time, this function returns `nothing`.
"""
function getOutFlux( mpSim::ManpowerSimulation, tBegin::T1, tEnd::T2,
    fields::Vector{String} = Vector{String}() ) where T1 <: Real where T2 <: Real

    if ( tEnd <= 0 ) || ( tBegin > min( now( mpSim ), mpSim.simLength ) ) ||
        ( tEnd <= tBegin )
        return
    end  # if ( tEnd <= 0 ) || ...

    tmpFields = unique( vcat( mpSim.idKey, fields ) )
    queryCmd = "SELECT '$(join( tmpFields, "', '" ))'
        FROM $(mpSim.personnelDBname)
        WHERE ( $tBegin < timeExited ) AND ( timeExited <= $tEnd ) AND
            ( timeEntered <= $tBegin )"
    SQLite.query( mpSim.simDB, queryCmd )

end  # getOutFlux( mpSim, tBegin, tEnd, fields )

"""
```
getOutFlux( mpSim::ManpowerSimulation,
            stateName::String,
            tBegin::T1,
            tEnd::T2 )
    where T1 <: Real
    where T2 <: Real
```
This function retrieves a view of all the personnel members in the simulation
`mpSim` who left the state with name `stateName` between `tBegin` (exclusive)
and `tEnd` (inclusive). People who entered and left the state (or the
organisation) in this interval are not counted. The view also includes the
target (end state) and the reason (transition) of the flux for each personnel
member.

If the end time is <= 0, or the start time is larger than the current simulation
time or the length of the simulation (whichever is smaller), or if the end time
is smaller than the start time, this function returns `nothing`.
"""
function getOutFlux( mpSim::ManpowerSimulation, stateName::String, tBegin::T1,
    tEnd::T2 ) where T1 <: Real where T2 <: Real

    if ( tEnd <= 0 ) || ( tBegin > min( now( mpSim ), mpSim.simLength ) ) ||
        ( tEnd <= tBegin )
        return
    end  # if ( tEnd <= 0 ) || ...

    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.personnelDBname)
        WHERE timeEntered < 0"
    dummyResult = SQLite.query( mpSim.simDB, queryCmd )

    # Get people who were in the state at the start of the time interval.
    inStateIDs = getActiveAtTime( mpSim, stateName, tBegin )[ 1 ]

    if isempty( inStateIDs )
        return dummyResult
    end  # if isempty( inStateIDs )

    # Get people who left the organisation in time interval.
    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.personnelDBname)
        WHERE ( $(mpSim.idKey) IN ('$(join( inStateIDs, "', '" ))') ) AND
            ( $tBegin < timeExited ) AND ( timeExited <= $tEnd )"
    fluxOutOfSystemIDs = SQLite.query( mpSim.simDB, queryCmd )[ 1 ]

    # Retrieve IDs and amount of transitions starting or ending in the state.
    queryCmd = "SELECT $(mpSim.idKey), count($(mpSim.idKey)) numEntries
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) IN ('$(join( inStateIDs, "', '" ))') ) AND
            ( ( startState IS '$stateName' ) OR ( endState IS '$stateName' ) )
        GROUP BY $(mpSim.idKey)"
    queryRes = SQLite.query( mpSim.simDB, queryCmd )
    nEntries, nCold = size( queryRes )
    numEntries = Dict{String, Vector{Int}}()

    for ii in 1:nEntries
        numEntries[ queryRes[ 1 ][ ii ] ] = [ queryRes[ 2 ][ ii ], 0 ]
    end  # for ii in 1:nEntries

    # Retrieve IDs and amount of transitions starting in the state.
    queryCmd = "SELECT $(mpSim.idKey), count($(mpSim.idKey)) numEntries
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) IN ('$(join( inStateIDs, "', '" ))') ) AND
            ( startState IS '$stateName' )
        GROUP BY $(mpSim.idKey)"
    queryRes = SQLite.query( mpSim.simDB, queryCmd )
    nEntries, nCold = size( queryRes )

    for ii in 1:nEntries
        numEntries[ queryRes[ 1 ][ ii ] ][ 2 ] = queryRes[ 2 ][ ii ]
    end  # for ii in 1:nEntries

    fluxOutIDs = collect( keys( numEntries ) )
    filter!( id -> numEntries[ id ][ 2 ] * 2 > numEntries[ id ][ 1 ],
        fluxOutIDs )
    fluxOutIDs = merge( fluxOutIDs, fluxOutOfSystemIDs )

    if isempty( fluxOutIDs )
        return dummyResult
    end  # if isempty( fluxInIDs )

    queryCmd = "SELECT $(mpSim.idKey), endState, transition
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) IN ('$(join( fluxOutIDs, "', '" ))') ) AND
            ( startState IN ('$stateName', 'active') )
        GROUP BY $(mpSim.idKey)"
    result = SQLite.query( mpSim.simDB, queryCmd )

    # If the person has multiple transition into/out of the state during time
    #   interval and hasn't left the system, set end state and transition to
    #   unknown.
    for ii in eachindex( fluxOutIDs )
        id = result[ Symbol( mpSim.idKey ) ][ ii ]

        if ( id ∉ fluxOutOfSystemIDs ) && ( numEntries[ id ][ 1 ] > 1 )
            result[ :endState ][ ii ] = "unknown"
            result[ :transition ][ ii ] = "unknown"
        end  # if id ∈ fluxOutOfSystemIDs
    end  # for ii in eachindex( fluxOutIDs )

    return result

    # Retrieve the list of those personnel members entering the state or the
    #   organisation in the time interval.
    if !isempty( fluxOutIDs )
        queryCmd = "SELECT $(mpSim.idKey)
            FROM $(mpSim.transitionDBname)
            WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
                ( $(mpSim.idKey) IN ('$(join( fluxOutIDs, "', '" ))') ) AND
                ( ( endState IS '$stateName' ) OR ( startState IS NULL ) )"
        fluxInIDs = SQLite.query( mpSim.simDB, queryCmd )[ 1 ]
        filter!( id -> id ∉ fluxInIDs, fluxOutIDs )
    end  # if !isempty( fluxOutIDs )

    # Get all personnel members who were in the state at the start of the time
    #   interval.
    inStateIDs = getActiveAtTime( mpSim, stateName, tBegin )[ 1 ]
    queryCmd = "SELECT $(mpSim.idKey)
        FROM $(mpSim.personnelDBname)
        WHERE ( $(mpSim.idKey) IN ('$(join( inStateIDs, "', '" ))') ) AND
            ( $tBegin < timeExited ) AND ( timeExited <= $tEnd )"
    fluxOutOfSystemIDs = SQLite.query( mpSim.simDB, queryCmd )[ 1 ]
    fluxOutIDs = merge( fluxOutIDs, fluxOutOfSystemIDs )

    if isempty( fluxOutIDs )
        return dummyResult
    end  # if isempty( fluxInIDs )

    queryCmd = "SELECT $(mpSim.idKey), endState, transition
        FROM $(mpSim.transitionDBname)
        WHERE ( $tBegin < timeIndex ) AND ( timeIndex <= $tEnd ) AND
            ( $(mpSim.idKey) IN ('$(join( fluxOutIDs, "', '" ))') ) AND
            ( startState IN ('$stateName', 'active') )"
    result2 = SQLite.query( mpSim.simDB, queryCmd )

    resMiss = filter( id -> id ∉ result[ 1 ], result2[ 1 ] )
    res2Miss = filter( id -> id ∉ result2[ 1 ], result[ 1 ] )

    print( "Time: $tBegin" )
    if isempty( resMiss ) && isempty( res2Miss )
        println( "\t\tNothing missing -- ", size( result ), " vs ", size( result2 ) )

        if length( result2[ 1 ] ) > length( unique( result2[ 1 ] ) )
            println( "Old way has duplicates." )
        end
    else
        if !isempty( resMiss )
            print( "\t\tMissing from new: ", resMiss )
        end

        if !isempty( res2Miss )
            print( "\t\tMissing from old: ", res2Miss )
        end

        println()
    end

    return SQLite.query( mpSim.simDB, queryCmd )

end  # getOutFlux( mpSim, tBegin, tEnd, fields )

"""
```
getLastChanges( mpSim::ManpowerSimulation,
                timePoint::T,
                fields::Vector{String} )
    where T <: Real
```
This function retrieves a view of the last changes in the simulation `mpSim` up
to time `timePoint` for all personnel members active at that time, and for the
fields in `fields`.

If the given time is < 0 or if the fields attribute is an empty vector, this
function returns `nothing`. Otherwise, it returns a `DataFrames.DataFrame`
object.
"""
function getLastChanges( mpSim::ManpowerSimulation, timePoint::T,
    fields::Vector{String} ) where T <: Real

    if ( timePoint < 0 ) || isempty( fields )
        return
    end  # if ( timePoint < 0 ) || ...

    persName = mpSim.personnelDBname
    histName = mpSim.historyDBname
    idKey = mpSim.idKey

    tmpFields = unique( fields )
    queryCmd = "SELECT $histName.$idKey, attribute, timeIndex, numValue,
        strValue FROM $histName
        INNER JOIN $persName ON $persName.$idKey == $histName.$idKey
        WHERE ( timeIndex <= $timePoint ) AND
            ( ( timeExited > $timePoint ) OR ( timeExited IS NULL ) ) AND
            ( attribute IN ('$(join( fields, "', '"))') )
        GROUP BY $histName.$idKey, attribute
        ORDER BY $histName.$idKey, attribute"
    return SQLite.query( mpSim.simDB, queryCmd )

end  # getLastChanges( mpSim, timePoint, fields )

"""
```
getDatabaseAtTime( mpSim::ManpowerSimulation,
                   timePoint::T,
                   fields::Vector{String} = Vector{String}() )
    where T <: Real
```
This function retrieves a view of all the personnel in the simulation `mpSim`
who were active at time `timePoint`, and displays the ID key and all the columns
listed in `fields`. The values in the database are those that the records had at
the given time of the simulation.

If the given time is < 0 or larger than the current simulation time or the
length of the simulation (whichever is smaller), or if the list of fields is
empty, this function returns `nothing`. Otherwise, it returns a
`DataFrames.DataFrame` object.
"""
function getDatabaseAtTime( mpSim::ManpowerSimulation, timePoint::T,
    fields::Vector{String} = Vector{String}() ) where T <: Real

    activeDB = getActiveAtTime( mpSim, timePoint, fields )

    # If no additional fields are requested, it's just the list of personnel
    #   members, and it can be returned.
    if isempty( fields )
        return activeDB
    end  # if isempty( fields )

    if activeDB === nothing
        return
    end  # if activeDB === nothing

    nPers = size( activeDB )[ 1 ]

    # If there aren't any active personnel members at that time, that's the
    #   database.
    if nPers == 0
        return activeDB
    end  # if nPers == 0

    changeList = getLastChanges( mpSim, timePoint, fields )
    nAttrs = floor( Int, size( changeList )[ 1 ] / nPers )
    attrsToChange = changeList[ 1:nAttrs, :attribute ]

    # XXX Thorough test needed to see if this works as advertised. If it does
    #   though... fantastic bit of coding.
    for ii in eachindex( attrsToChange )
        # The attribute to update.
        changedAttr = Symbol( attrsToChange[ ii ] )
        # Is it a numerical attribute?
        isNum = isa( changeList[ ii, :numValue ], Real )
        # Perform the update.
        activeDB[ changedAttr ] =
            changeList[ ii + (0:(nPers - 1)) * nAttrs,
            isNum ? :numValue : :strValue ]
        # This works because the view from the personnel database is sorted by
        #   personnel ID, and the view from the historic database is sorted by
        #   personnel ID, then by attribute.
        # In addition, any attribute that shows up in the historic database
        #   must show up for ALL personnel records as the starting value of the
        #   attribute must be stored for later recall.
    end  # for ii in eachindex( fields )

    return activeDB

end  # getDatabaseAtTime( mpSim, timePoint, fields )

export generateReport
"""
This function generates a basic simulation report of the simulation `mpSim` on
a timegrid with resolution `timeDelta`, and exports it to the Excel file
`fileName`. The `.xlsx` extension is added if necessary.
"""
function generateReport( mpSim::ManpowerSimulation, timeDelta::T,
    fileName::String ) where T <: Real

    tic()
    nRecords = countRecords( mpSim, timeDelta )
    nFluxIn = countFluxIn( mpSim, timeDelta )
    nFluxOut = countFluxOut( mpSim, timeDelta, true )

    workbook = Workbook()
    sheet = createSheet( workbook, "Personnel" )

    # General info
    sheet[ "A", 1 ] = "Simulation length"
    sheet[ "B", 1 ] = min( now( mpSim ), mpSim.simLength )
    sheet[ "A", 2 ] = "Report time resolution"
    sheet[ "B", 2 ] = timeDelta
    sheet[ "A", 3 ] = "Personnel cap"
    sheet[ "B", 3 ] = mpSim.personnelTarget
    sheet[ "A", 4 ] = "Report generated in"
    sheet[ "C", 4 ] = "s"

    # Data headers
    sheet[ "A", 6 ] = "Time"
    sheet[ "B", 6 ] = "Personnel"
    sheet[ "C", 6 ] = "Flux In"
    sheet[ "D", 6 ] = "Flux Out"
    sheet[ "E", 6 ] = "Net Flux"
    sheet[ "F", 6 ] = "Retired"
    sheet[ "G", 6 ] = "Resigned"

    # Put data in Excel file
    n = length( nRecords[ 1 ] )
    sheet[ "A", 7 ] = nRecords[ 1 ][ 1 ]
    sheet[ "B", 7 ] = nRecords[ 2 ][ 1 ]

    for ii in 2:n
        sheet[ "A", 6 + ii ] = nRecords[ 1 ][ ii ]
        sheet[ "B", 6 + ii ] = nRecords[ 2 ][ ii ]
        sheet[ "C", 6 + ii ] = nFluxIn[ 2 ][ ii - 1 ]
        sheet[ "D", 6 + ii ] = nFluxOut[ 2 ][ ii - 1 ]
        sheet[ "E", 6 + ii ] = sheet[ "C", 6 + ii ] - sheet[ "D", 6 + ii ]
        sheet[ "F", 6 + ii ] = haskey( nFluxOut[ 3 ], "retirement" ) ? nFluxOut[ 3 ][ "retirement" ][ ii - 1 ] : 0
        sheet[ "G", 6 + ii ] = haskey( nFluxOut[ 3 ], "attrition" ) ? nFluxOut[ 3 ][ "attrition" ][ ii - 1 ] : 0
    end  # for ii in 2:n

    # if n == length( nFluxIn[ 1 ] )
    #     sheet[ "A", 7 + n ] == nFluxIn[ 1 ][ n ]
    #     sheet[ "C", 7 + n ] == nFluxIn[ 2 ][ n ]
    #     sheet[ "D", 7 + n ] == nFluxOut[ 2 ][ n ]
    # end  # if n == length( nFluxIn[ 1 ] )

    # Wrap up.
    sheet[ "B", 4 ] = toq()  # Time it took to create the report.
    tmpFileName = endswith( fileName, ".xlsx" ) ? fileName : fileName * ".xlsx"
    write( tmpFileName, workbook )
    println( "Report created and saved to $tmpFileName." )

end  # generateReport( mpSim, timeDelta )


text:
function petk02()
#
# M. Petkovic testing polynomials, p118
#
    p = poly([1.0;1;-im;-im;-im;5*im;5*im;-5*im;-5*im]);
    z = [1 2; -im 3; 5*im 2; -5*im 2];

    p, PolyZeros(z)
end

ext: jl


text:
# This file is a part of TextSearch.jl

module TextSearch
    import Base: broadcastable
    import StatsBase: fit, predict
    using SimilaritySearch, InvertedFiles, LinearAlgebra, SparseArrays

    include("textconfig.jl")
    include("normalize.jl")
    include("tokenize.jl")
    include("bow.jl")
    include("vmodel.jl")
    include("emodel.jl")
    include("neardup.jl")
#    include("multi.jl")
end


text:
import DoubleFloats: Double64, ComplexDF64

@testset "LinearAlgebra" begin
	@testset "Constructor + basics" begin
		A = rand(ComplexF32, 12, 12)
		WS = HC.MatrixWorkspace(A)
		@test WS isa HC.MatrixWorkspace{Float64}
		@test_throws ArgumentError HC.update!(WS, rand(12, 13))

		WS.factorized[] = true
		B = randn(ComplexF64, 12, 12)
		HC.update!(WS, B)
		@test WS.A == B
		@test WS.A !== B
		@test WS.factorized[] == false

		HC.factorization!(WS, HC.LU_FACT)
		@test WS.fact[] == HC.LU_FACT
		WS.factorized[] = true
		HC.factorization!(WS, HC.QR_FACT)
		@test WS.fact[] == HC.QR_FACT
		@test WS.factorized[] == false

		HC.factorization!(WS, HC.LU_FACT)
		@test (@allocated HC.factorize!(WS)) == 0
		HC.factorization!(WS, HC.QR_FACT)
		@test WS.factorized[] == false
		HC.factorize!(WS) # first factorize resizes the work buffers
		WS.factorized[] = false
		@test (@allocated HC.factorize!(WS)) == 0

		# test indexing etc
		@test WS[2,1] == B[2, 1]
		@test WS[2] == B[2]
		WS[2,3] = 5.0
		@test WS[2,3] == 5.0
		WS[5] = 2.32
		@test WS[5] == 2.32
		@test size(WS) == (12,12)
		C = randn(ComplexF64, 12, 12)
		WS .= C
		@test WS.A == C
	end

	@testset "ldiv" begin
		A = randn(ComplexF64, 13, 13)
		b = randn(ComplexF64, 13)
		x = zeros(ComplexF64, 13)
		WS = HC.MatrixWorkspace(A)
		HC.update!(WS, A)
		ldiv!(x, WS, b)
		WS.factorized[] = false
		@test (@allocated ldiv!(x, WS, b)) == 0
		@test (lu(A) \ b) ≈ x atol=1e-12

		HC.factorization!(WS, HC.QR_FACT)
		x .= 0
		ldiv!(x, WS, b)
		@test (qr(A, Val(true)) \ b) ≈ x atol=1e-12
		WS.factorized[] = false
		@test (@allocated ldiv!(x, WS, b)) == 0

		# overdetermined
		A = randn(ComplexF64, 13, 8)
		x = randn(ComplexF64, 8)
		b = A*x
		WS = HC.MatrixWorkspace(A)
		HC.update!(WS, A)
		@test_throws ArgumentError HC.factorization!(WS, HC.LU_FACT)
		ldiv!(x, WS, b)
		WS.factorized[] = false
		x .= 0
		@test (@allocated ldiv!(x, WS, b)) == 0
		@test (qr(A, Val(true)) \ b) ≈ x atol=1e-12
	end

	@testset "(Mixed precision) residual" begin
		A = randn(ComplexF64, 13, 13)
		x = randn(ComplexF64, 13)
		b = A * x
		û = convert.(ComplexF64, ComplexDF64.(A) * ComplexDF64.(x) - ComplexDF64.(b))

		u = zeros(ComplexF64, 13)
		# fixed precision residual should be of order unit roundoff
		HC.residual!(u, A, x, b)
		@test maximum(abs, ComplexDF64.(u) - ComplexDF64.(û)) < 1e-14
		@test maximum(abs, ComplexDF64.(u) - ComplexDF64.(û)) > 1e-18

		HC.residual!(u, A, x, b, Complex{Double64})
		@test maximum(abs, ComplexDF64.(u) - ComplexDF64.(û)) < 1e-26
	end

	@testset "(Mixed precision) iterative refinement" begin
		A = randn(ComplexF64, 4, 4)
		b = randn(ComplexF64, 4)
		x̂ = zeros(ComplexF64, 4)
		x = copy(x̂)
		WS = HC.MatrixWorkspace(A)
		U, S, VT = svd(A)
		S[end] *= 1e-4
		A2 = U * diagm(0=>S) * VT'
		update!(WS, A2)
		ldiv!(x̂, WS, b)
		δx1 = HC.iterative_refinement_step!(x, WS, x̂, b, HC.InfNorm(), ComplexDF64)
		@test δx1 / norm(x, Inf) > 1e-14
		δx2 = HC.iterative_refinement_step!(WS, x, b, HC.InfNorm(), ComplexDF64)
		@test δx2/ norm(x, Inf) < 1e-16

		δx = HC.iterative_refinement_step!(x, WS, x̂, b)
		@test δx1 / norm(x, Inf) > 1e-14

		# overdetermined
		A = randn(ComplexF64, 4, 3)
		x̂ = zeros(ComplexF64, 3)
		x = copy(x̂)
		WS = HC.MatrixWorkspace(A)
		U, S, VT = svd(A)
		S[end] *= 1e-4
		A2 = U * diagm(0=>S) * VT'
		update!(WS, A2)
		b = A2 * randn(ComplexF64, 3)
		ldiv!(x̂, WS, b)
		δx1 = HC.iterative_refinement_step!(x, WS, x̂, b, HC.InfNorm(), ComplexDF64)
		@test δx1 / norm(x, Inf) > 1e-14
		δx2 = HC.iterative_refinement_step!(WS, x, b, HC.InfNorm(), ComplexDF64)
		@test δx2/ norm(x, Inf) < 1e-16

		δx = HC.iterative_refinement_step!(x, WS, x̂, b)
		@test δx1 / norm(x, Inf) > 1e-14
	end

	@testset "rcond" begin
		A = randn(ComplexF64, 13, 13)
		WS = HC.MatrixWorkspace(A)
		U, S, VT = svd(A)
		S[end] *= 1e-5
		A2 = U * diagm(0=>S) * VT'
		HC.update!(WS, A2)
		HC.factorization!(WS, HC.LU_FACT)
		HC.factorize!(WS)
		@test cond(A2, Inf) ≈ inv(HC.rcond(WS))
		HC.factorization!(WS, HC.QR_FACT)
		@test 0.5 ≤ cond(A2, Inf) / inv(HC.rcond(WS)) ≤ 1.5
		# @test (@allocated rcond(WS)) == 0

		A = randn(ComplexF64, 13, 10)
		WS = HC.MatrixWorkspace(A)
		U, S, VT = svd(A)
		S[end] *= 1e-8
		A2 = U * diagm(0=>S) * VT'
		HC.update!(WS, A2)
		HC.factorize!(WS)
		@test 1 / HC.rcond(WS) ≈ cond(UpperTriangular(WS.qr.R), Inf) rtol=1e-14
		# @test (@allocated rcond(WS)) == 0

		# meddle aroudn with the internals to check dead branch
		WS.fact[] = HC.LU_FACT
		@test isnan(HC.rcond(WS))
	end

	@testset "inf norm est" begin
		A = randn(ComplexF64, 5, 5)
		U, s, Vt = svd(A)
		s[end] *= 1e-6
		B = U*diagm(0=>s)*Vt'

		d = rand(5)*1000
		D = diagm(0=>d)
		g = rand(5) .* 1000
		G = diagm(0=>g)

		@test 0.1 ≤ opnorm(inv(B), Inf) / HC.inf_norm_est(lu(B)) ≤ 10
		@test 0.1 ≤ opnorm(inv(D)*inv(B), Inf) / HC.inf_norm_est(lu(B), nothing, d) ≤ 10
		@test 0.1 ≤ opnorm(inv(B)*G, Inf) / HC.inf_norm_est(lu(B), g, nothing) ≤ 10
		@test 0.1 ≤ opnorm(inv(D)*inv(B)*G, Inf) / HC.inf_norm_est(lu(B), g, d) ≤ 10

		WS = HC.MatrixWorkspace(B)
		@test 0.1 ≤ opnorm(inv(B), Inf) / HC.inf_norm_est(WS) ≤ 10
		@test 0.1 ≤ opnorm(inv(D)*inv(B), Inf) / HC.inf_norm_est(WS, nothing, d) ≤ 10
		@test 0.1 ≤ opnorm(inv(B)*G, Inf) / HC.inf_norm_est(WS, g, nothing) ≤ 10
		@test 0.1 ≤ opnorm(inv(D)*inv(B)*G, Inf) / HC.inf_norm_est(WS, g, d) ≤ 10
	end

	@testset "JacobianMonitor" begin
		A = randn(ComplexF64, 6, 6)
		b = randn(ComplexF64, 6)
		x̂ = similar(b)
		x = similar(b)
		JM = HC.JacobianMonitor(zeros(ComplexF64, 6, 6))
		@test JM isa HC.JacobianMonitor
		HC.jacobian(JM) .= A
		HC.updated!(JM)
		ldiv!(x̂, JM, b)
		ldiv!(x, jacobian(JM), b)
		@test x == x̂
	end
	@testset "Forward error" begin
		A = [1.81451+0.46473im -0.473651+0.813117im 0.130822-0.254704im -0.203277+0.397106im;
			 -2.04468-1.59228im 0.441157-1.39994im -0.0984804+1.42619im -0.818292-0.487898im;
			 0.35501+0.0235574im 1.32618+0.577264im 1.12561+0.43383im 0.579533+1.49484im;
			 0.144381+0.104175im 0.665384-1.43389im 0.22097+0.258016im 0.374103-0.635066im]
		b = [-0.14581-0.868847im, 0.0162454+0.91156im, -0.133769+0.510834im, -0.867022+0.752196im]
		x̂ = similar(b)
		x = similar(b)
		JM = HC.JacobianMonitor(A)
		HC.updated!(JM)
		ldiv!(x̂, JM, b)

		ferr = HC.forward_err!(x, JM, x̂, b, HC.InfNorm())
		@test HC.forward_err(JM) == ferr
		@test ferr < eps()*cond(A)
		@test ferr > 1e-12
		ferrD64 = HC.forward_err!(x, JM, x̂, b, HC.InfNorm(), ComplexDF64)
		@test 0.5 ≤ ferrD64 / ferr ≤ 2

		ferr = HC.forward_err!(JM, x̂, b, HC.InfNorm(), ComplexF64)
		@test ferr < eps()*cond(A)
		@test ferr > 1e-12

		D = diagm(0 => [10.0^(2i) for i in 1:4])
		JM = HC.JacobianMonitor(D*A)
		HC.updated!(JM)
		ldiv!(x̂, JM, D*b)
		conderr = HC.cond(D*A)*eps()
		sferr = HC.strong_forward_err!(x, JM, x̂, D*b, HC.InfNorm())
		@test HC.forward_err(JM) == sferr
		ferr = HC.forward_err!(x, JM, x̂, D*b, HC.InfNorm())
		@test ferr < sferr < conderr*1e-4 # check that error estimate is actually substantially less

		sferr2 = HC.strong_forward_err!(JM, x̂, D*b, HC.InfNorm())
		@test sferr == sferr2
		@test x̂ ≈ x

		DR = diagm(0=>[10.0^(2i-4) for i in 1:4])
		JM = HC.JacobianMonitor(A*DR)
		ldiv!(x̂, JM, b)

		WN = HC.WeightedNorm(HC.InfNorm(), x̂)
		WN .= abs.(x̂)

		wferr = HC.forward_err!(x, JM, x̂, b, WN)
		ferr = HC.forward_err!(x, JM, x̂, b, HC.InfNorm())

		strong_wferr = HC.strong_forward_err!(x, JM, x̂, b, WN)
		strong_ferr = HC.strong_forward_err!(x, JM, x̂, b, HC.InfNorm())
		@test ferr ≤ strong_ferr ≤ eps()*inv(HC.rcond(jacobian(JM)))
		@test wferr ≤ strong_wferr ≤ eps()*inv(HC.rcond(jacobian(JM)))

		# overdetermined
		A = randn(ComplexF64, 3, 2)
		x = randn(ComplexF64, 2)
		b = A * x
		x̂ = similar(x)

		JM = HC.JacobianMonitor(A)
		ldiv!(x̂, JM, b)

		ferr = HC.forward_err!(x, JM, x̂, b, HC.InfNorm())
		strong_ferr = HC.strong_forward_err!(x, JM, x̂, b, HC.InfNorm())
		@test ferr < 1e-15
		@test ferr == strong_ferr # test fallback to ferr for qr
	end

	@testset "cond" begin
		A = randn(ComplexF64, 13, 13)
		D = diagm(0=>[2.0^(i-5) for i in 1:13])
		JM = HC.JacobianMonitor(D*A)
		comp_cond = HC.cond!(JM)
		@test comp_cond === cond(JM)
		B = diagm(0=>inv.(JM.J.row_scaling)) * D * A
		@test abs(cond(JM) - cond(B, Inf)) / cond(JM) < 2

		A = randn(ComplexF64, 13, 12)
		JM = HC.JacobianMonitor(D*A)
		comp_cond = HC.cond!(JM)
		@test comp_cond === cond(JM)
	end

	@testset "Hermite Normal Form" begin
        A = [0 3 1 0; -2 2 -1 2; 1 -1 2 3; -3 3 3 2]
        H, U = HC.hnf(A)
        @test A * U == H

        D = zeros(Int, 1,1)
        D[1,1] = -2
        H, U = HC.hnf(D)
        @test H[1,1] == 2
        @test U[1,1] == -1

		# overflow A
		A = [2 4 4 2 -4 -3 -4 3 0 3;
		     -4 -2 -5 -5 1 0 -3 0 -1 -2;
			 -3 5 -3 1 4 3 -1 -4 -1 0;
			 1 3 3 4 -2 -3 -2 -5 -4 -3;
			 2 2 0 2 -4 4 3 -4 -2 -4;
			 -3 -3 1 1 -4 4 0 4 4 -4;
			 -3 -1 -1 -1 2 -4 -3 -4 4 -4;
			 0 1 -2 4 5 4 3 1 -5 2;
			 -5 -4 -5 3 1 5 0 -3 -3 -1;
			 2 -5 -3 -1 -1 5 -4 5 -3 1]
		@test_throws OverflowError HC.hnf(A)
		@test_throws OverflowError HC.hnf(A, Int128)
		H, U = HC.hnf(A, BigInt)
		@test A * U == H
    end
end


text:
# Setup loaded nodes entry
function update_equivalent_m_and_f!(data::Dict)
  ns = Int[]
  if :loaded_elements in keys(data)
    for t in data[:loaded_nodes]
      push!(ns, t[1])
    end
    for t in data[:loaded_elements]
      if t[1] in ns
        indx = findall(in(t[1]), ns)[1]
        data[:loaded_nodes][indx][2][1] += -2.8
        data[:loaded_nodes][indx][2][2] += -0.8
        t2 = t[1] + 1
        if !(t2 in ns)
          push!(ns, t2)
          push!(data[:loaded_nodes], (t2, [0.0  0.0]))
        end
        indx = findall(in(t2), ns)[1]
        data[:loaded_nodes][indx][2][1] += -1.2
        data[:loaded_nodes][indx][2][2] += 0.5333
      else
        force = -6.0
        moment = 3.0
        if !(t[1] in ns)
          push!(data[:loaded_nodes], (t[1], [force -moment]))
          push!(ns, t[1])
        else
          indx = findall(in(t[1]), ns)[1]
          data[:loaded_nodes][indx][2][1] += -force
          data[:loaded_nodes][indx][2][2] += -moment
        end
        t2 = t[1] + 1
        if !(t2 in ns)
          push!(data[:loaded_nodes], (t2, [force  moment]))
          push!(ns, t2)
        else
          indx = findall(in(t2), ns)[1]
          data[:loaded_nodes][indx][2][1] += force
          data[:loaded_nodes][indx][2][2] += moment
        end
      end
    end
  end
end

update_equivalent_m_and_f!(data::Dict)

data[:loaded_nodes]


text:
"""
    calc_nchk(n::Integer,k::Integer)

Calculates binomial coefficient: n choose k
"""
function calc_nchk(n::Integer,k::Integer)
    accum::Int = 1
    for i in 1:k
        accum = accum * (n-k+i) ÷ i
    end
    return accum
end

"""
    get_nchk(n::Integer,k::Integer)

Looks up binomial coefficient from a precomputed table: n choose k
"""
@inline function get_nchk(n,k)
    return binom_coeff[n+1,k+1]
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
mutable struct TTState
  currentunits :: Units
end

global function _traverseTree(e::Element, c::Channel, ::TTState)
  put!(c, e)
end
global function _traverseTree(u::Units, c::Channel, state::TTState)
  put!(c, u)
  state.currentunits = u
end
global function _traverseTree(g::Group, c::Channel, state::TTState)
  previousunits = state.currentunits
  if g.units != Units()
    put!(c, g.units)
    state.currentunits = g.units
  end
  for e in g.elements
    _traverseTree(e, c, state)
  end
  if g.units != Units()
    put!(c, previousunits)
    state.currentunits = previousunits
  end
end
function __traverseTree(e,c)
  _traverseTree(e,c,TTState(Units("mm")))
  close(c)
end

function traverseTree(e::Element)
  c = Channel{Element}(32)
  @async __traverseTree(e,c)
  return c
end


text:
export to_basictypes

"""
    to_basictypes(block::AbstractBlock)

convert gates to basic types

    * ChainBlock
    * PutBlock
    * PrimitiveBlock
"""
function to_basictypes end

to_basictypes(block::PrimitiveBlock) = block
function to_basictypes(block::AbstractBlock)
    throw(NotImplementedError(:to_basictypes, typeof(block)))
end

function to_basictypes(block::RepeatedBlock)
    chain(block.n, map(i -> put(block.n, i => content(block)), block.locs))
end

to_basictypes(block::CachedBlock) = content(block)
function to_basictypes(block::Subroutine{D,<:PrimitiveBlock}) where {D}
    put(nqudits(block), block.locs => content(block))
end
function to_basictypes(block::Subroutine)
    to_basictypes(map_address(content(block), AddressInfo(block.n, [block.locs...])))
end
function to_basictypes(block::Subroutine{D,<:Measure}) where {D}
    map_address(content(block), AddressInfo(nqudits(block), [block.locs...]))
end
to_basictypes(block::Daggered) = Daggered(block.content)
to_basictypes(block::Scale) = Scale(block.alpha, block.content)
to_basictypes(block::KronBlock) =
    chain(block.n, [put(block.n, i => block[i]) for i in block.locs])
to_basictypes(block::Union{Add,PutBlock,ChainBlock,ControlBlock}) = block


text:
function sinmf()

end


text:
#####
##### Nothing happens when your boundary condition is nothing
#####

fill_west_halo!(c,   ::Nothing, args...; kwargs...) = NoneEvent()
fill_east_halo!(c,   ::Nothing, args...; kwargs...) = NoneEvent()
fill_south_halo!(c,  ::Nothing, args...; kwargs...) = NoneEvent()
fill_north_halo!(c,  ::Nothing, args...; kwargs...) = NoneEvent()
fill_top_halo!(c,    ::Nothing, args...; kwargs...) = NoneEvent()
fill_bottom_halo!(c, ::Nothing, args...; kwargs...) = NoneEvent()

fill_west_and_east_halo!(c,  ::Nothing, ::Nothing, args...; kwargs...) = NoneEvent()
fill_south_and_north_halo!(c,::Nothing, ::Nothing, args...; kwargs...) = NoneEvent()
fill_bottom_and_top_halo!(c, ::Nothing, ::Nothing, args...; kwargs...) = NoneEvent()


text:
"""
A structure that holds all of the relevant information for constructing the model
"""
struct GLO{T1<:Real, T2<:Integer}
	kernel::Function  # kernel function
	n_kern_hyper::T2  # amount of hyperparameters for the kernel function
	n_dif::T2  # amount of times you are differenting the base kernel
	n_out::T2  # amount of scores you are jointly modelling
	x_obs::Vector{T1} # the observation times/phases
	y_obs::Vector{T1}  # the flattened, observed data
	noise::Vector{T1}  # the measurement noise at all observations
	normals::Vector{T1}  # the normalization of each section of y_obs
	a0::Matrix{T1}  # the meta kernel coefficients
	non_zero_hyper_inds::Vector{T2}  # the indices of the non-zero hyperparameters
	# The powers that each a0 coefficient
	# is taken to for each part of the matrix construction
	# used for constructing differentiated versions of the kernel
	coeff_orders::AbstractArray{T2,6}
	coeff_coeffs::AbstractArray{T2,4}
	covariance::Array{T1, 3}  # the measurement covariance at all observations
	has_covariance::Bool
	kernel_changes_with_output::Bool

	function GLO(
		kernel::Function,
		n_kern_hyper::T2,
		n_dif::T2,
		n_out::T2,
		x_obs::Vector{T1},
		y_obs::Vector{T1};
		noise::Vector{T1} = zeros(length(y_obs)),
		normals::Vector{T1} = ones(n_out),
		a0::Matrix{T1} = ones(n_out, n_dif),
		covariance::Array{T1, 3} = zeros(length(x_obs), n_out, n_out),
		kernel_changes_with_output::Bool=false
		) where {T1<:Real, T2<:Integer}

		@assert size(a0) == (n_out, n_dif)
		non_zero_hyper_inds = append!(findall(!iszero, collect(Iterators.flatten(a0))), collect(1:n_kern_hyper) .+ length(a0))
		coeff_orders, coeff_coeffs = coefficient_orders(n_out, n_dif, a=a0)
		has_covariance = (covariance != zeros(length(x_obs), n_out, n_out))

		return GLO(kernel, n_kern_hyper, n_dif, n_out, x_obs, y_obs, noise, normals, a0, non_zero_hyper_inds, coeff_orders, coeff_coeffs, covariance, has_covariance, kernel_changes_with_output)
	end
	function GLO(
		kernel::Function,  # kernel function
		n_kern_hyper::T2,  # amount of hyperparameters for the kernel function
		n_dif::T2,  # amount of times you are differenting the base kernel
		n_out::T2,  # amount of scores you are jointly modelling
		x_obs::Vector{T1}, # the observation times/phases
		y_obs::Vector{T1},  # the flattened, observed data
		noise::Vector{T1},  # the measurement noise at all observations
		normals::Vector{T1},  # the normalization of each section of y_obs
		a0::Matrix{T1},
		non_zero_hyper_inds::Vector{T2},
		coeff_orders::Array{T2,6},
		coeff_coeffs::Array{T2,4},
		covariance::Array{T1, 3},
		has_covariance::Bool,
		kernel_changes_with_output::Bool) where {T1<:Real, T2<:Integer}

		@assert isfinite(kernel(ones(n_kern_hyper), randn(), zeros(Int64, 2 + n_kern_hyper)))  # make sure the kernel is valid by testing a sample input
		@assert 0 < n_dif <= 3
		@assert 0 < n_out
		n_meas = length(x_obs)
		@assert (n_meas * n_out) == length(y_obs) == length(noise)
		@assert n_meas == size(covariance, 1)
		@assert length(normals) == n_out
		@assert size(a0) == (n_out, n_dif)
		@assert size(coeff_orders) == (n_out, n_out, n_dif, n_dif, n_out, n_dif)  # maybe unnecessary due to the fact that we construct it
		@assert size(coeff_coeffs) == (n_out, n_out, n_dif, n_dif)  # maybe unnecessary due to the fact that we construct it
		@assert length(non_zero_hyper_inds) == length(findall(!iszero, collect(Iterators.flatten(a0)))) + n_kern_hyper
		@assert n_out == size(covariance, 2) == size(covariance, 3)
		@assert (covariance != zeros(n_meas, n_out, n_out)) == has_covariance

		return new{typeof(x_obs[1]),typeof(n_kern_hyper)}(kernel, n_kern_hyper, n_dif, n_out, x_obs, y_obs, noise, normals, a0, non_zero_hyper_inds, coeff_orders, coeff_coeffs, covariance, has_covariance, kernel_changes_with_output)
	end
end

function normalize_problem_definition!(prob_def::GLO)
	n_obs = length(prob_def.x_obs)
	renorms = ones(prob_def.n_out)
	for i in 1:prob_def.n_out
		inds = (i:prob_def.n_out:length(prob_def.y_obs))
		prob_def.y_obs[inds] .-= mean(prob_def.y_obs[inds])
		renorms[i] = std(prob_def.y_obs[inds])
		# println(renorm)
		prob_def.normals[i] *= renorms[i]
		prob_def.y_obs[inds] /= renorms[i]
		prob_def.noise[inds] /= renorms[i]
	end
	if prob_def.has_covariance
		renorm_mat = renorms .* transpose(renorms)
		for i in 1:n_obs
			prob_def.covariance[i, :, :] ./= renorm_mat
		end
	end
end


text:
module GraphLasso
using DataFrames, Lasso, Graphs

export graphlasso

function graphlasso(S::Matrix{Float64}, α::Float64; tol::Float64=1e-5,
                    maxit::Int=1000, penalize_diag::Bool=true, verbose::Bool=false)
    p = size(S,1)
    adj = abs(S) .> α
    blocks = connected_components(sparse2adjacencylist(sparse(max(adj, eye(p)))))
    if(verbose)
        print(blocks)
    end
    W = zeros(p,p)
    Θ = zeros(p,p)
    for block = blocks
        W[block,block], Θ[block,block] = fit_block(S[block,block], α, tol, maxit, penalize_diag)
    end
    return W, Θ
end

function fit_block(S::Matrix{Float64}, α::Float64, tol::Float64,
                   maxit::Int, penalize_diag::Bool)
    p = size(S,1)
    W = copy(S)
    if penalize_diag
        W += α * eye(p)
    end
    if size(S,1) == 1
        return W, inv(W)
    end
    Θ = zeros(p,p)
    i = 0
    β = zeros(p-1,p)
    while i < maxit
        i += 1
        W_old = copy(W)
        for j = 1:p
            inds = collect(1:p)
            splice!(inds, j)
            W11 = W[inds,inds]
            sqrtW11 = sqrtm(W11)
            β[:,j] = fit(LassoPath, sqrtW11, sqrtW11 \ S[inds,j], λ=[α/(p-1)],
                         standardize=false, intercept=false).coefs
            W[inds,j] = W11 * β[:,j]
            W[j,inds] = W[inds,j]
        end
        if norm(W - W_old) < tol
            break
        end
    end
    if i == maxit
        warn("Maximum number of iterations reached, graphlasso failed to converge")
    end
    for j in 1:p
        inds = collect(1:p)
        splice!(inds, j)
        Θ[j,j] = 1/(W[j,j] - dot(W[inds,j],β[:,j]))
        Θ[inds,j] = -Θ[j,j] * β[:,j]
        Θ[j,inds] = Θ[inds,j]
    end
    return W, Θ
end
end


text:
#using ACEds.RPI: SparsePSHDegree, BasicPSH1pBasis

using ACE
using ACE: numz
#import ACE: standardevaluator, graphevaluator
using ACEds
using ACEds.RPI: RPIBasis, SparsePSHDegree, BasicPSH1pBasis, PIBasis, Rot3DCoeffs, _rpi_A2B_matrix
using Random, Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d, evaluate_ed
using JuLIP.MLIPs: combine
using StaticArrays

#---

@info("Basic test of RPIBasis construction and evaluation")
maxdeg = 10
N = 4
r0 = 1.0
rcut = 3.0
trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
D = SparsePSHDegree()
P1 = BasicPSH1pBasis(Pr; species = :X, D = D)

#---
pibasis = PIBasis(P1, N, D, maxdeg)
#rpibasis = RPIBasis(P1, N, D, maxdeg)
rpibasis = RPIBasis(P1, N, D, maxdeg,L=1)
#---

#---

using ReferenceFrameRotations: angle_to_dcm

Q = angle_to_dcm(.5,.4,-.3, :ZYX)

Nat = 10
Rs, Zs, z0 = rand_nhd(Nat, Pr, :X)

Rs_rot = deepcopy(Rs)
for i = 1:Nat
    Rs_rot[i] = Q*Rs_rot[i]
end

B = evaluate(rpibasis, Rs, Zs, z0)

#print(Q*real.(B[1]), "     ", B_rot[1],"\n")
N_basis = length(B)

@info("check for equivariance")

B_rot = evaluate(rpibasis, Rs_rot, Zs, z0)
#for i in 1:N_basis
#    print(norm(Q*B[i]-B_rot[i])<1E-5,"\n")
#    print(Q*B[i], "     ", B_rot[i],"\n")
#end
@test all([ norm(Q*B[i]-B_rot[i])<1E-5 for i in 1:length(B)])

@info("check for linear independence")

using LinearAlgebra: rank

N_eval = Integer(ceil(N_basis/3)) # may pass
N_eval = Integer(floor(N_basis/3)) # should not pass
N_eval = N_basis  # should pass

B_array = zeros(3*N_eval,N_basis)
atol = 1E-10
for k in 1:N_eval
    Rs, Zs, z0 = rand_nhd(Nat, Pr, :X)
    B = evaluate(rpibasis, Rs, Zs, z0)
    for i in 1:N_basis
        B_array[(3*(k-1)+1):(3*(k-1)+3), i] = B[i]
    end
end

print("Rank: ", rank(B_array,atol=atol),"\n")
@test rank(B_array,atol=atol) == N_basis


text:
"""
Run this 2D-demo with 4 procs : mpirun -n 4 julia demo_2D.jl
"""
module demo
using MPIStructuredMesh
using MPI # this is not necessary, only for demo display purpose

"""
MPI display is bad, this is just a display function for the demo
"""
function display_array(array, mesh)
    println("\n")
    nprocs = MPI.Comm_size(mesh.comm)
    for r in 0:nprocs-1
        if get_rank(mesh) == r
            println(mesh.coords) # print the mesh coordinates in the cartesian grid
            display(transpose(array)) # display array
            println("\n")
        end
        MPI.Barrier(mesh.comm)
    end
end

# Build mesh
n = (4, 3) # Number of elements in each space directions
ndims = (2, 2) # Number of processors in each space direction
noverlaps = (2,1) # Number of elements overlap (on both side) in each space direction
mesh = DistributedMesh(ndims, n, noverlaps)

# Build an array for demo. Only interior elements are initialized (they received a unique integer identifier)
array = zeros(Int, nelts_with_halo(mesh))
irange, jrange = local_indices(mesh)
nx, ny = get_nelts(mesh)
for i in 1:nx, j in 1:ny
    array[irange[i],jrange[j]] = (get_rank(mesh) + 1) * 100 + (j - 1) * nx + i
end

# Print the array before exchange
println("\nBefore update_halo!\n")
display_array(array, mesh)

# Update halo !
update_halo!(array, mesh)

# Print array after exchange
println("\nAfter update_halo!\n")
display_array(array, mesh)

# The end
finalize_mesh(mesh)

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule QPALM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("QPALM")
JLLWrappers.@generate_main_file("QPALM", UUID("ffff4089-a68e-5e40-8b55-c7b9cd3999ca"))
end  # module QPALM_jll


text:
module BayesInference

using LinearAlgebra
using SpecialFunctions

using Distributions
using PDMats
using StatsBase
using Dates

include("generalTools.jl")
include("conjugateUpdates.jl")
include("mcmcTools.jl")

end # module


text:
function ground(video, action, framenumber)
    background("white")
    sethue("blue")
    return framenumber
end

function ground_color(color_bg, color_pen, framenumber)
    background(color_bg)
    sethue(color_pen)
    return framenumber
end

function latex_title(args...)
    translate(0, -200)
    latex(L"E=mc^2", 20)
end

function circ(p=O, color="black")
    sethue(color)
    circle(p, 25, :fill)
    return p
end

function circ_ret_trans(p=O, color="black")
    sethue(color)
    circle(p, 25, :fill)
    return Transformation(p, 0.0)
end

function path!(points, pos, color)
    sethue(color)
    push!(points, pos)
    circle.(points, 2, :fill)
end

function rad(p1, p2, color)
    sethue(color)
    line(p1,p2, :stroke)
end

@testset "Dancing circles" begin
    p1 = Point(100,0)
    p2 = Point(100,80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(video, [
        BackgroundAction(1:25, ground, Rotation(0.0), Translation(Point(25, 25), Point(25, 25))),
        Action(latex_title),
        Action(Rel(-24:0), :red_ball, (args...)->circ(p1, "red"), Rotation(from_rot, to_rot)),
        Action(1:25, :blue_ball, (args...)->circ(p2, "blue"), Rotation(to_rot, from_rot, :red_ball)),
        Action(1:25, (video, args...)->path!(path_of_red, get_position(:red_ball), "red")),
        Action(:same, (video, args...)->path!(path_of_blue, pos(:blue_ball), "blue")),
        Action(1:25, (args...)->rad(pos(:red_ball), pos(:blue_ball), "black"))
    ], tempdirectory="images", pathname="dancing.gif")

    @test_reference "refs/dancing_circles_16.png" load("images/0000000016.png")
    @test isfile("dancing.gif")
    for i=1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
    rm("images/palette.bmp")
    rm("dancing.gif")
end

@testset "Dancing circles layered" begin
    p1 = Point(100,0)
    p2 = Point(100,80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(video, [
        Action(1:25, ground, Rotation(π/2, π/2, O), Translation(Point(25,25), Point(25,25)); in_global_layer=true),
        Action(1:25, latex_title),
        Action(1:25, :red_ball, (args...)->circ(p1, "red"), Rotation(to_rot)),
        Action(1:25, :blue_ball, (args...)->circ(p2, "blue"), Rotation(to_rot, from_rot, :red_ball)),
        Action(1:25, (video, args...)->path!(path_of_red, get_position(:red_ball), "red")),
        Action(1:25, (video, args...)->path!(path_of_blue, pos(:blue_ball), "blue")),
        Action(1:25, (args...)->rad(pos(:red_ball), pos(:blue_ball), "black"))
    ], tempdirectory="images", pathname="")

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i=1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Dancing circles layered return Transformation" begin
    p1 = Point(100,0)
    p2 = Point(100,80)
    from_rot = 0.0
    to_rot = 2π
    path_of_blue = Point[]
    path_of_red = Point[]

    video = Video(500, 500)
    javis(video, [
        Action(1:25, ground, Rotation(π/2, π/2, O), Translation(Point(25,25), Point(25,25)); in_global_layer=true),
        Action(1:25, latex_title),
        Action(1:25, :red_ball, (args...)->circ_ret_trans(p1, "red"), Rotation(to_rot)),
        Action(1:25, :blue_ball, (args...)->circ_ret_trans(p2, "blue"), Rotation(to_rot, from_rot, :red_ball)),
        Action(1:25, (video, args...)->path!(path_of_red, pos(:red_ball), "red")),
        Action(1:25, (video, args...)->path!(path_of_blue, pos(:blue_ball), "blue")),
        Action(1:25, (args...)->rad(pos(:red_ball), pos(:blue_ball), "black"))
    ], tempdirectory="images", pathname="")

    @test_reference "refs/dancing_circles_16_rot_trans.png" load("images/0000000016.png")
    for i=1:25
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Drawing grid" begin
    video = Video(500, 500)
    javis(
        video,
        [
            Action(1:40, ground),

            Action(1:10, draw_grid(direction = "BL", line_gap = 25)),
            Action(zero_lines(direction = "BL", line_thickness = 10)),

            Action(Rel(10), draw_grid(direction = "BR", line_gap = 25)),
            Action(zero_lines(direction = "BR", line_thickness = 10)),

            Action(Rel(10), draw_grid(direction = "TL", line_gap = 25)),
            Action(zero_lines(direction = "TL", line_thickness = 10)),

            Action(Rel(10), draw_grid(direction = "TR", line_gap = 25)),
            Action(zero_lines(direction = "TR", line_thickness = 10)),

        ],
        tempdirectory = "images", pathname=""
    )

    @test_reference "refs/grid_drawing_bl.png" load("images/0000000008.png")
    @test_reference "refs/grid_drawing_br.png" load("images/0000000018.png")
    @test_reference "refs/grid_drawing_tl.png" load("images/0000000028.png")
    @test_reference "refs/grid_drawing_tr.png" load("images/0000000038.png")
    for i=1:40
	    rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "Use returned angle" begin
    video = Video(500, 500)
    p = Point(100, 0)
    javis(video, [
        Action(1:10, ground),
        Action(:circ, (args...)->circ_ret_trans(p), Rotation(2π)),
        Action((args...)->line(Point(-200, 0), Point(-200, -10*ang(:circ)), :stroke))
    ], tempdirectory="images", pathname="")

    @test_reference "refs/circle_angle.png" load("images/0000000008.png")
    for i=1:10
	    rm("images/$(lpad(i, 10, "0")).png")
    end
end

astar(args...) = star(Point(-100,-100), 30)
acirc(args...) = circle(Point(100,100), 30)

@testset "morphing star2circle and back" begin
    video = Video(500, 500)
    javis(video, [
        BackgroundAction(1:20, :framenumber, (args...)->ground_color("white", "black", args[3])),
        Action(1:10, (args...)->circle(Point(-100, 0), val(:framenumber), :fill)),
        Action(1:10, morph(astar, acirc)),
        Action(11:20, morph(acirc, astar))
    ], tempdirectory="images", pathname="")

    @test_reference "refs/star2circle5.png" load("images/0000000005.png")
    @test_reference "refs/star2circle15.png" load("images/0000000015.png")
    for i=1:20
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_nicholas(args...)
    background("white")
    sethue("blue")
    setline(3)
end

function house_of_nicholas(;p1=O, width=100, color="black")
    # sethue(color)
    #     .p1
    # .p2   .p3
    #
    # .p4   .p5
    width2 = div(width, 2)
    p2 = p1 + Point(-width2, width2)
    p3 = p1 + Point(width2, width2)
    p4 = p2 + Point(0, width)
    p5 = p3 + Point(0, width)
    line(p4, p5, :stroke)
    line(p5, p2, :stroke)
    line(p2, p4, :stroke)
    line(p4, p3, :stroke)
    line(p3, p1, :stroke)
    line(p1, p2, :stroke)
    line(p2, p3, :stroke)
    setline(8)
    line(p3, p5, :stroke)
end

@testset "House of Nicholas line_width" begin
    demo = Video(500, 500)
    javis(demo, [
        BackgroundAction(1:50, ground_nicholas),
        Action((args...)->house_of_nicholas(); subactions = [
            SubAction(1:25, appear(:fade_line_width)),
            SubAction(26:50, disappear(:fade_line_width))
        ])
    ], tempdirectory="images", pathname="")

    @test_reference "refs/nicholas15.png" load("images/0000000015.png")
    @test_reference "refs/nicholas25.png" load("images/0000000025.png")
    @test_reference "refs/nicholas35.png" load("images/0000000035.png")
    for i=1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

function ground_opacity(args...)
    background("white")
    sethue("black")
    setopacity(0.5)
end

function square_opacity(p1, w)
    setopacity(1.0)
    sethue("red")
    rect(p1, w, w, :fill)
end

@testset "Circle/square appear opacity" begin
    demo = Video(500, 500)
    javis(demo, [
        BackgroundAction(1:50, ground_opacity),
        Action((args...)->circ(); subactions = [
            SubAction(1:25, appear(:fade)),
            SubAction(26:50, disappear(:fade))
        ]),
        Action(5:50, (args...)->square_opacity(Point(-100, 0), 60); subactions = [
            SubAction(1:15, appear(:fade)),
            SubAction(16:35, Translation(100, 50)),
            SubAction(36:40, disappear(:fade))
            # for global frames 46-50 it should still be disappeared
        ])
    ], tempdirectory="images", pathname="")

    @test_reference "refs/circlerSquare07opacity.png" load("images/0000000007.png")
    @test_reference "refs/circlerSquare25opacity.png" load("images/0000000025.png")
    @test_reference "refs/circlerSquare42opacity.png" load("images/0000000042.png")
    # test that the last frame is completely white
    @test sum(load("images/0000000050.png")) == RGB{Float64}(500*500, 500*500, 500*500)
    for i=1:50
        rm("images/$(lpad(i, 10, "0")).png")
    end
end

@testset "test default kwargs" begin
    video = Video(500, 500)
    pathname = javis(video, [
        Action(1:10, ground),
        Action(1:10, morph(astar, acirc))
    ])
    path, ext = splitext(pathname)
    @test ext == ".gif"
    @test isfile(pathname)
    rm(pathname)
end

@testset "test @error .mp3" begin
    video = Video(500, 500)
    @test_logs (:error,) javis(video, [
        Action(1:10, ground),
        Action(1:10, morph(astar, acirc))
    ]; pathname="test.mp3")
end


text:
using LinearAlgebra
using BenchmarkTools

dgg = zeros(4, 4)

dge = [
    0. 0 1 0;
    0 0 0 0;
    0 0 0 0;
    0 0 0 0
]

dee = [
    -1. 0 0 0;
    0 -1 0 0;
    0 0 -1 0;
    0 0 0 -1
]

Z = zeros(4, 4)

d  = [dgg Z; Z dee]
d′ = [Z dge; Z Z]

ρ = [
    1. 1 1 1 1 1 0 1;
    1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1;
    0 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 0
]
ρ₀ = similar(ρ); ρ₀ .= 0
ρ_tmp = similar(ρ); ρ_tmp .= 0

D = (d' * d)::Array{Float64, 2}
d′t = transpose(d)

function test1!(ρ₀, ρ, ρ_tmp, d, d′, d′t, D)
    mul!(ρ_tmp, ρ, d′t)
    mul!(ρ₀, d′, ρ_tmp)
    mul!(ρ₀, D, ρ, -1/2, 0)
    mul!(ρ₀, ρ, D, -1/2, 0)
    return nothing
end

function test2!(ρ₀, ρ, ρ_tmp, d, d′, d′t, D)
    BLAS.gemm!('N', 'T', 1., ρ,  d′, 1., ρ_tmp)
    BLAS.gemm!('N', 'N', 1., ρ₀, d′, 1., ρ_tmp)
    BLAS.gemm!('N', 'N', -1/2, D, ρ, 1., ρ₀)
    BLAS.gemm!('N', 'N', -1/2, ρ, D, 1., ρ₀)
    return nothing
end

@benchmark test1!($ρ₀, $ρ, $ρ_tmp, $d, $d′, $d′t, $D) evals=1
@benchmark test2!($ρ₀, $ρ, $ρ_tmp, $d, $d′, $d′t, $D) evals=1

test2!(ρ₀, ρ, ρ_tmp, d, d′, d′t, D); ρ₀

function test3!(t, d, ρ)

    @inbounds for ga in 1:4
        @inbounds for eb in 5:8
            t[ga, eb] = -ρ[ga, eb] / 2
            t[eb, ga] = t[ga, eb]
        end
    end

    @inbounds for ea in 5:8
        @inbounds for eb in 5:8
            t[ea, eb] = -ρ[ea, eb]
        end
    end

    @inbounds for ga in 1:4
        @inbounds for gb in 1:4
            @inbounds for ec′ in 5:8
                @inbounds for ec′′ in 5:8
                    t[ga, gb] += d[ec′, ga] * d[gb, ec′′] * ρ[ec′, ec′′]
                end
            end
        end
    end

    return nothing
end

@benchmark test3!(t_, $d, $ρ) setup=(t_ = zeros(8, 8)) evals=1

@time test3!(t, d, ρ)


text:
using BinaryBuilder

name = "FreeBSDLibc"
version = v"11.1"

# sources to build, such as mingw32, our patches, etc....
sources = [
    "https://download.freebsd.org/ftp/releases/amd64/11.1-RELEASE/base.txz" =>
    "62acaee7e7c9df66ee2c0c2d533d1da0ddf67d32833bc4b77d935ddd9fe27dab"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
sysroot="${prefix}/${target}/sys-root"

mkdir -p "${sysroot}"
mv usr/include "${sysroot}/"
mv usr/lib "${sysroot}/"
mv lib/* "${sysroot}/lib"
mkdir -p "${sysroot}/usr"
ln -sf "${sysroot}/include" "${sysroot}/usr/"
ln -sf "${sysroot}/lib" "${sysroot}/usr/"
ln -sf "libgcc_s.so.1" "${sysroot}/lib/libgcc_s.so"
ln -sf "libcxxrt.so.1" "${sysroot}/lib/libcxxrt.so"

# Many symlinks exist that point to `../../lib/libfoo.so`.
# We need them to point to just `libfoo.so`. :P
for f in $(find "${prefix}/${target}" -xtype l); do
	link_target="$(readlink "$f")"
	if [[ -n $(echo "${link_target}" | grep "^../../lib") ]]; then
		ln -vsf "${link_target#../../lib/}" "${f}"
	fi
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    FreeBSD(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libc", :libc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

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


text:
module InelasticScattering

    using DocStringExtensions
    using Parameters
    using OffsetArrays
    using Interpolations
    #import PhysicalConstants.CODATA2018:c_0, h, k_B

    include("raman_constants.jl")
    include("molecular_constructors.jl")
    include("inelastic_cross_section.jl")
    include("apply_lineshape.jl")

    "Speed of light in vacuum in `[cm/s]`"
    global const c = 2.99792458e10
    "Planck constant in `[erg ̇s]`"
    global const h = 6.62607015e-27
    "Boltzmann constant in `[erg/K]`"
    global const k_B = 1.380649e-16

    export compute_effective_coefficents!
    export compute_σ_Rayl_coeff!
    export compute_σ_Rayl_VibRaman_coeff_hires!
    export compute_energy_levels!
    export apply_lineshape!

    #export compute_σ_Raman_coeff!
    export compute_σ_VibRaman_coeff!, compute_σ_RoVibRaman_coeff!
end # module

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
using Test
using IncrementalInference

##

@testset "Basic test of belief prediction on alternate solveKey" begin

##

fg = initfg()

addVariable!(fg, :a, ContinuousScalar)
addVariable!(fg, :b, ContinuousScalar)

addFactor!(fg, [:a], Prior(Normal(10,1)), graphinit=false)
addFactor!(fg, [:a;:b], LinearRelative(Normal(10,1)), graphinit=false)

deleteVariableSolverData!(fg, :a, :default)
deleteVariableSolverData!(fg, :b, :default)

##

pts = sampleFactor(fg, :af1, 100)

IIF.setDefaultNodeData!(getVariable(fg, :a), 0, 100, 1, solveKey=:testSolveKey,
                        initialized=false, varType=ContinuousScalar(), dontmargin=false)
#

initManual!(fg, :a, pts, :testSolveKey)

@test !isInitialized(fg, :a)
@test isInitialized(fg, :a, :testSolveKey)

##

IIF.setDefaultNodeData!(getVariable(fg, :b), 0, 100, 1, solveKey=:testSolveKey,
                        initialized=false, varType=ContinuousScalar(), dontmargin=false)
#

@test !(:default in listSolveKeys(getVariable(fg, :a)))
@test !(:default in listSolveKeys(getVariable(fg, :b)))

@test (:testSolveKey in listSolveKeys(getVariable(fg, :a)))
@test (:testSolveKey in listSolveKeys(getVariable(fg, :b)))

##

doautoinit!(fg, :b, solveKey=:testSolveKey)

##

@test !(:default in listSolveKeys(getVariable(fg, :a)))
@test !(:default in listSolveKeys(getVariable(fg, :b)))

@test (:testSolveKey in listSolveKeys(getVariable(fg, :a)))
@test (:testSolveKey in listSolveKeys(getVariable(fg, :b)))

##

@test isapprox( calcPPE(fg, :a, solveKey=:testSolveKey).suggested[1], 10, atol=1)
@test isapprox( calcPPE(fg, :b, solveKey=:testSolveKey).suggested[1], 20, atol=1)

##

end

@testset "test solve with unique solveKey, see #1219" begin

##

fg = initfg()
getSolverParams(fg).graphinit=false

addVariable!(fg, :a, ContinuousScalar)
addVariable!(fg, :b, ContinuousScalar)
addVariable!(fg, :c, ContinuousScalar)
addVariable!(fg, :d, ContinuousScalar)
addVariable!(fg, :e, ContinuousScalar)

addFactor!(fg, [:a], Prior(Normal()))
addFactor!(fg, [:a;:b], LinearRelative(Normal(10, 1)))
addFactor!(fg, [:b;:c], LinearRelative(Normal(10, 1)))
addFactor!(fg, [:c;:d], LinearRelative(Normal(10, 1)))
addFactor!(fg, [:d;:e], LinearRelative(Normal(10, 1)))

getSolverParams(fg).graphinit=true

##

# getSolverParams(fg).limititers=30
solveTree!(fg, solveKey=:testSolveKey )

##

@test isapprox( calcPPE(fg, :a, solveKey=:testSolveKey).suggested[1], 0, atol=2)
@test isapprox( calcPPE(fg, :b, solveKey=:testSolveKey).suggested[1], 10, atol=2)
@test isapprox( calcPPE(fg, :c, solveKey=:testSolveKey).suggested[1], 20, atol=2)
@test isapprox( calcPPE(fg, :d, solveKey=:testSolveKey).suggested[1], 30, atol=2)
@test isapprox( calcPPE(fg, :e, solveKey=:testSolveKey).suggested[1], 40, atol=2)

##

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm,25cm)

# ##

# plotKDE(fg, ls(fg), solveKey=:testSolveKey)

##

end

#


text:
export Configuration;
export stack_depth, buffer_length, config_push!, config_pop!, unshift!, add_arc, is_unshiftable, set_reshiftable;
export has_head, del_arc;

export NONEXIST_TOKEN, UNKNOWN_TOKEN, NULL_TOKEN, ROOT_TOKEN

using DataStructures

const NONEXIST_TOKEN = -1
const UNKNOWN_TOKEN = "-UNKNOWN-"
const NULL_TOKEN = "-NULL-"
const LIMIT_TRANSITIONS_NUMBER = 400

mutable struct Configuration
  buffer::Vector{Integer}
  stack::Stack{Integer}
  sentence::Sentence
  unshiftable::Vector{Bool}
  tree::DependencyTree

  Configuration(
    buffer::Vector{Integer},
    stack::Stack{Integer},
    sentence::Sentence,
    unshiftable::Vector{Bool},
    tree::DependencyTree
  ) = new(buffer, stack, sentence, unshiftable, tree)

  function Configuration(sentence::Sentence)
    buffer = collect(1:sentence.length)
    unshiftable = zeros(Bool, sentence.length)
    stack = Stack{Integer}()
    push!(stack, 0)
    tree = DependencyTree(sentence)

    new(buffer, stack, sentence, unshiftable, tree)
  end
end

stack_depth(config::Configuration) = length(config.stack)

buffer_length(config::Configuration) = length(config.buffer)

function config_push!(config::Configuration)
  buffer_length(config) == 0 && return

  buffer_element = popfirst!(config.buffer)
  push!(config.stack, buffer_element)
end

function config_pop!(config::Configuration)
  stack_depth(config) == 0 && return

  pop!(config.stack)
end

function unshift!(config::Configuration)
  word_id = config_pop!(config)
  if word_id === nothing
    return
  end

  insert!(config.buffer, 1, word_id)
  config.unshiftable[word_id] = true
end

function add_arc(config::Configuration, head_id::Integer, word_id::Integer, label::String)
  set_arc(config.tree, word_id, head_id, label)
end

function is_unshiftable(config::Configuration, word_id::Integer)
  1 <= word_id <= length(config.unshiftable) ? config.unshiftable[word_id] : true
end

function set_reshiftable(config::Configuration, word_id::Integer)
  if 1 <= word_id <= length(config.unshiftable)
    config.unshiftable[word_id] = false
  end
end

function get_token(config::Configuration, index::Integer)
  index == 0 && return ROOT_TOKEN

  1 <= index <= config.sentence.length && return config.sentence.tokens[index].name

  NULL_TOKEN
end

function get_tag(config::Configuration, index::Integer)
  index == 0 && return ROOT_TOKEN

  1 <= index <= config.sentence.length && return config.sentence.pos_tags[index].name

  NULL_TOKEN
end

function get_label(config::Configuration, index::Integer)
  1 <= index <= config.sentence.length && return config.tree.nodes[index].label

  NULL_TOKEN
end

function get_buffer_element(config::Configuration, index::Integer)
  1 <= index <= buffer_length(config) && return config.buffer[index]

  NONEXIST_TOKEN
end

function get_stack_element(config::Configuration, index::Integer)
  stack_vector = collect(config.stack)
  1 <= index <= stack_depth(config) && return stack_vector[index]

  NONEXIST_TOKEN
end

function is_terminal(config::Configuration)
  stack_depth(config) == 1 && buffer_length(config) == 0
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
export BrickletAnalogInV3VoltageCallbackConfiguration
struct BrickletAnalogInV3VoltageCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletAnalogInV3Calibration
struct BrickletAnalogInV3Calibration
    offset::Integer
    multiplier::Integer
    divisor::Integer
end

export BrickletAnalogInV3SPITFPErrorCount
struct BrickletAnalogInV3SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletAnalogInV3Identity
struct BrickletAnalogInV3Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAnalogInV3
"""
Measures DC voltage between 0V and 42V
"""
mutable struct BrickletAnalogInV3 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAnalogInV3(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_analog_in_v3")
        deviceInternal = package.BrickletAnalogInV3(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_voltage
"""
    $(SIGNATURES)

Returns the measured voltage. The resolution is approximately 10mV to 1mV
depending on the oversampling configuration (:func:`Set Oversampling`).

If you want to get the value periodically, it is recommended to use the
:cb:`Voltage` callback. You can set the callback configuration
with :func:`Set Voltage Callback Configuration`.
"""
function get_voltage(device::BrickletAnalogInV3)
    return device.deviceInternal.get_voltage()
end

export set_voltage_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Voltage` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Voltage` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_voltage_callback_configuration(device::BrickletAnalogInV3, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_voltage_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_voltage_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Voltage Callback Configuration`.
"""
function get_voltage_callback_configuration(device::BrickletAnalogInV3)
    return device.deviceInternal.get_voltage_callback_configuration()
end

export set_oversampling
"""
    $(SIGNATURES)

Sets the oversampling between 32x and 16384x. The Bricklet
takes one 12bit sample every 17.5µs. Thus an oversampling
of 32x is equivalent to an integration time of 0.56ms and
a oversampling of 16384x is equivalent to an integration
time of 286ms.

The oversampling uses the moving average principle. A
new value is always calculated once per millisecond.

With increased oversampling the noise decreases. With decreased
oversampling the reaction time increases (changes in voltage will be
measured faster).
"""
function set_oversampling(device::BrickletAnalogInV3, oversampling)
    device.deviceInternal.set_oversampling(oversampling)
end

export get_oversampling
"""
    $(SIGNATURES)

Returns the oversampling value as set by :func:`Set Oversampling`.
"""
function get_oversampling(device::BrickletAnalogInV3)
    return device.deviceInternal.get_oversampling()
end

export set_calibration
"""
    $(SIGNATURES)

Sets a calibration for the measured voltage value.
The formula for the calibration is as follows::

 Calibrated Value = (Value + Offset) * Multiplier / Divisor

We recommend that you use the Brick Viewer to calibrate
the Bricklet. The calibration will be saved internally and only
has to be done once.
"""
function set_calibration(device::BrickletAnalogInV3, offset, multiplier, divisor)
    device.deviceInternal.set_calibration(offset, multiplier, divisor)
end

export get_calibration
"""
    $(SIGNATURES)

Returns the calibration as set by :func:`Set Calibration`.
"""
function get_calibration(device::BrickletAnalogInV3)
    return device.deviceInternal.get_calibration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletAnalogInV3)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletAnalogInV3, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletAnalogInV3)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletAnalogInV3, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletAnalogInV3, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletAnalogInV3, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletAnalogInV3)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletAnalogInV3)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletAnalogInV3)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletAnalogInV3, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletAnalogInV3)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAnalogInV3)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAnalogInV3, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
__precompile__(true)

module LearnBase

# Only reexport required functions by default
import StatsBase: nobs, fit, fit!, predict, params, params!
using Compat

import Base.issymmetric

"""
Baseclass for any kind of cost. Notable examples for
costs are `Loss` and `Penalty`.
"""
@compat abstract type Cost end

"""
Baseclass for all losses. A loss is some (possibly simplified)
function `L(features, targets, outputs)`, where `outputs` are the
result of some function `f(features)`.
"""
@compat abstract type Loss <: Cost end

"""
A loss is considered **supervised**, if all the information needed
to compute `L(features, targets, outputs)` are contained in
`targets` and `outputs`, and thus allows for the simplification
`L(targets, outputs)`.
"""
@compat abstract type SupervisedLoss <: Loss end

"""
A supervised loss, where the targets are in {-1, 1}, and which
can be simplified to `L(targets, outputs) = L(targets * outputs)`
is considered **margin-based**.
"""
@compat abstract type MarginLoss <: SupervisedLoss end

"""
A supervised loss that can be simplified to
`L(targets, outputs) = L(targets - outputs)` is considered
distance-based.
"""
@compat abstract type DistanceLoss <: SupervisedLoss end

"""
A loss is considered **unsupervised**, if all the information needed
to compute `L(features, targets, outputs)` are contained in
`features` and `outputs`, and thus allows for the simplification
`L(features, outputs)`.
"""
@compat abstract type UnsupervisedLoss <: Loss end

@compat abstract type Penalty <: Cost end

function scaled end

function value end
function value! end
function meanvalue end
function sumvalue end
function meanderiv end
function sumderiv end
function deriv end
function deriv2 end
function deriv! end
function value_deriv end
function addgrad! end
function value_grad end
function value_grad! end
function prox end
function prox! end

"Return the gradient of the learnable parameters w.r.t. some objective"
function grad end
"Do a backward pass, updating the gradients of learnable parameters and/or inputs"
function grad! end

function isminimizable end
function isdifferentiable end
function istwicedifferentiable end
function isconvex end
function isstrictlyconvex end
function isstronglyconvex end
function isnemitski end
function isunivfishercons end
function isfishercons end
function islipschitzcont end
function islocallylipschitzcont end
function islipschitzcont_deriv end # maybe overkill
function isclipable end
function ismarginbased end
function isclasscalibrated end
function isdistancebased end

# fallback to supersets
isstrictlyconvex(any) = isstronglyconvex(any)
isconvex(any) = isstrictlyconvex(any)
islocallylipschitzcont(any) = islipschitzcont(any)

"""
Anything that takes an input and performs some kind
of function to produce an output. For example a linear
prediction function.
"""
@compat abstract type Transformation end
@compat abstract type StochasticTransformation <: Transformation end
@compat abstract type Learnable <: Transformation end

function transform end
"Do a forward pass, and return the output"
function transform! end

"""
Baseclass for any prediction model that can be minimized.
This means that an object of a subclass contains all the
information needed to compute its own current loss.
"""
@compat abstract type Minimizable <: Learnable end

function update end
function update! end
function learn end
function learn! end

#    getobs(data, [idx], [obsdim])
#
# Return the observations corresponding to the observation-index
# `idx`. Note that `idx` can be of type `Int` or `AbstractVector`.
# Both options must be supported by a custom type.
#
# The returned observation(s) should be in the form intended to
# be passed as-is to some learning algorithm. There is no strict
# interface requirement on how this "actual data" must look like.
# Every author behind some custom data container can make this
# decision him-/herself. We do, however, expect it to be consistent
# for `idx` being an integer, as well as `idx` being an abstract
# vector, respectively.
#
# If it makes sense for the type of `data`, `obsdim` can be used
# to disptach on which dimension of `data` denotes the observations.
# See `?ObsDim`
#
# This function is implemented in MLDataPattern
function getobs end

#    getobs!(buffer, data, [idx], [obsdim])
#
# Inplace version of `getobs(data, idx, obsdim)`. If this method
# is defined for the type of `data`, then `buffer` should be used
# to store the result, instead of allocating a dedicated object.
#
# Implementing this function is optional. In the case no such
# method is provided for the type of `data`, then `buffer` will be
# *ignored* and the result of `getobs` returned. This could be
# because the type of `data` may not lend itself to the concept
# of `copy!`. Thus supporting a custom `getobs!(::MyType, ...)`
# is optional and not required.
#
# If it makes sense for the type of `data`, `obsdim` can be used
# to disptach on which dimension of `data` denotes the observations.
# See `?ObsDim`
#
# This function is implemented in MLDataPattern
function getobs! end

# --------------------------------------------------------------------

#    gettarget([f], observation)
#
# Use `f` (if provided) to extract the target from the single
# `observation` and return it. It is used internally by
# `targets` (only if `f` is provided) and by
# `eachtarget` (always) on each individual observation.
#
# Even though this function is not exported, it is intended to be
# extended by users to support their custom data storage types.
#
# This function is implemented in MLDataPattern
function gettarget end # not exported

#    gettargets(data, [idx], [obsdim])
#
# Return the targets corresponding to the observation-index `idx`.
# Note that `idx` can be of type `Int` or `AbstractVector`.
#
# Implementing this function for a custom type of `data` is
# optional. It is particularly useful if the targets in `data` can
# be provided without invoking `getobs`. For example if you have a
# remote data-source where the labels are part of some metadata
# that is locally available.
#
# If it makes sense for the type of `data`, `obsdim` can be used
# to disptach on which dimension of `data` denotes the observations.
# See `?ObsDim`
#
# This function is implemented in MLDataPattern
function gettargets end # not exported

#    targets([f], data, [obsdim])
#
# This function is implemented in MLDataPattern
function targets end

# --------------------------------------------------------------------

"""
    abstract DataView{TElem, TData} <: AbstractVector{TElem}

Baseclass for all vector-like views of some data structure.
This allow for example to see some design matrix as a vector of
individual observation-vectors instead of one matrix.

see `MLDataPattern.ObsView` and `MLDataPattern.BatchView` for examples.
"""
@compat abstract type DataView{TElem, TData} <: AbstractVector{TElem} end

"""
    abstract AbstractObsView{TElem, TData} <: DataView{TElem, TData}

Baseclass for all vector-like views of some data structure,
that views it as some form or vector of observations.

see `MLDataPattern.ObsView` for a concrete example.
"""
@compat abstract type AbstractObsView{TElem, TData} <: DataView{TElem, TData} end

"""
    abstract AbstractBatchView{TElem, TData} <: DataView{TElem, TData}

Baseclass for all vector-like views of some data structure,
that views it as some form or vector of equally sized batches.

see `MLDataPattern.BatchView` for a concrete example.
"""
@compat abstract type AbstractBatchView{TElem, TData} <: DataView{TElem, TData} end

# --------------------------------------------------------------------

"""
    abstract DataIterator{TElem,TData}

Baseclass for all types that iterate over a `data` source
in some manner. The total number of observations may or may
not be known or defined and in general there is no contract that
`getobs` or `nobs` has to be supported by the type of `data`.
Furthermore, `length` should be used to query how many elements
the iterator can provide, while `nobs` may return the underlying
true amount of observations available (if known).

see `MLDataPattern.RandomObs`, `MLDataPattern.RandomBatches`
"""
@compat abstract type DataIterator{TElem,TData} end

"""
    abstract ObsIterator{TElem,TData} <: DataIterator{TElem,TData}

Baseclass for all types that iterate over some data source
one observation at a time.

```julia
using MLDataPattern
@assert typeof(RandomObs(X)) <: ObsIterator

for x in RandomObs(X)
    # ...
end
```

see `MLDataPattern.RandomObs`
"""
@compat abstract type ObsIterator{TElem,TData} <: DataIterator{TElem,TData} end

"""
    abstract BatchIterator{TElem,TData} <: DataIterator{TElem,TData}

Baseclass for all types that iterate over of some data source one
batch at a time.

```julia
@assert typeof(RandomBatches(X, size=10)) <: BatchIterator

for x in RandomBatches(X, size=10)
    @assert nobs(x) == 10
    # ...
end
```

see `MLDataPattern.RandomBatches`
"""
@compat abstract type BatchIterator{TElem,TData} <: DataIterator{TElem,TData} end

# --------------------------------------------------------------------

# just for dispatch for those who care to
@compat const AbstractDataIterator{E,T}  = Union{DataIterator{E,T}, DataView{E,T}}
@compat const AbstractObsIterator{E,T}   = Union{ObsIterator{E,T},  AbstractObsView{E,T}}
@compat const AbstractBatchIterator{E,T} = Union{BatchIterator{E,T},AbstractBatchView{E,T}}

# --------------------------------------------------------------------

#    datasubset(data, [idx], [obsdim])
#
# Return a lazy subset of the observations in `data` that correspond
# to the given `idx`. No data should be copied except of the
# indices. Note that `idx` can be of type `Int` or `AbstractVector`.
# Both options must be supported by a custom type.
#
# If it makes sense for the type of `data`, `obsdim` can be used
# to disptach on which dimension of `data` denotes the observations.
# See `?ObsDim`
#
# This function is implemented in MLDataPattern
function datasubset end

"""
    default_obsdim(data)

The specify the default obsdim for a specific type of data.
Defaults to `ObsDim.Undefined()`
"""
function default_obsdim end

# just for dispatch for those who care to
"see `?ObsDim`"
@compat abstract type ObsDimension end

"""
    module ObsDim

Singleton types to define which dimension of some data structure
(e.g. some `Array`) denotes the observations.

- `ObsDim.First()`
- `ObsDim.Last()`
- `ObsDim.Contant(dim)`

Used for efficient dispatching
"""
module ObsDim
    using ..LearnBase.ObsDimension

    """
    Default value for most functions. Denotes that the concept of
    an observation dimension is not defined for the given data.
    """
    immutable Undefined <: ObsDimension end

    """
        ObsDim.Last <: ObsDimension

    Defines that the last dimension denotes the observations
    """
    immutable Last <: ObsDimension end

    """
        ObsDim.Constant{DIM} <: ObsDimension

    Defines that the dimension `DIM` denotes the observations
    """
    immutable Constant{DIM} <: ObsDimension end
    Constant(dim::Int) = Constant{dim}()

    """
        ObsDim.First <: ObsDimension

    Defines that the first dimension denotes the observations
    """
    const First = Constant{1}
end

Base.convert(::Type{ObsDimension}, dim) = throw(ArgumentError("Unknown way to specify a obsdim: $dim"))
Base.convert(::Type{ObsDimension}, dim::ObsDimension) = dim
Base.convert(::Type{ObsDimension}, ::Void) = ObsDim.Undefined()
Base.convert(::Type{ObsDimension}, dim::Int) = ObsDim.Constant(dim)
Base.convert(::Type{ObsDimension}, dim::String) = convert(ObsDimension, Symbol(lowercase(dim)))
Base.convert(::Type{ObsDimension}, dims::Tuple) = map(d->convert(ObsDimension, d), dims)
function Base.convert(::Type{ObsDimension}, dim::Symbol)
    if dim == :first || dim == :begin
        ObsDim.First()
    elseif dim == Symbol("end") || dim == :last
        ObsDim.Last()
    elseif dim == Symbol("nothing") || dim == :none || dim == :null || dim == :na || dim == :undefined
        ObsDim.Undefined()
    else
        throw(ArgumentError("Unknown way to specify a obsdim: $dim"))
    end
end

@deprecate obs_dim(dim) convert(ObsDimension, dim)

default_obsdim(data) = ObsDim.Undefined()
default_obsdim(A::AbstractArray) = ObsDim.Last()
default_obsdim(tup::Tuple) = map(default_obsdim, tup)

import Base: AbstractSet

"A continuous range (inclusive) between a lo and a hi"
immutable IntervalSet{T} <: AbstractSet{T}
    lo::T
    hi::T
end
function IntervalSet{A,B}(lo::A, hi::B)
    T = promote_type(A,B)
    IntervalSet{T}(convert(T,lo), convert(T,hi))
end

# numeric interval
randtype{T<:Number}(s::IntervalSet{T}) = Float64
Base.rand{T<:Number}(s::IntervalSet{T}, dims::Integer...) = rand(dims...) * (s.hi - s.lo) + s.lo
Base.in{T<:Number}(x::Number, s::IntervalSet{T}) = s.lo <= x <= s.hi
Base.length{T<:Number}(s::IntervalSet{T}) = 1

# vector of intervals
randtype{T<:AbstractVector}(s::IntervalSet{T}) = Vector{Float64}
Base.rand{T<:AbstractVector}(s::IntervalSet{T}) = Float64[rand() * (s.hi[i] - s.lo[i]) + s.lo[i] for i=1:length(s)]
Base.in{T<:AbstractVector}(x::AbstractVector, s::IntervalSet{T}) = all(i -> s.lo[i] <= x[i] <= s.hi[i], 1:length(s))
Base.length{T<:AbstractVector}(s::IntervalSet{T}) = length(s.lo)

"Set of discrete items"
immutable DiscreteSet{T<:AbstractArray} <: AbstractSet{T}
    items::T
end
randtype(s::DiscreteSet) = eltype(s.items)
Base.rand(s::DiscreteSet, dims::Integer...) = rand(s.items, dims...)
Base.in(x, s::DiscreteSet) = x in s.items
Base.length(s::DiscreteSet) = length(s.items)
Base.getindex(s::DiscreteSet, i::Int) = s.items[i]

# operations on arrays of sets
randtype{S<:AbstractSet,N}(sets::AbstractArray{S,N}) = Array{promote_type(map(randtype, sets)...), N}
Base.rand{S<:AbstractSet}(sets::AbstractArray{S}) = eltype(randtype(sets))[rand(s) for s in sets]
function Base.rand{S<:AbstractSet}(sets::AbstractArray{S}, dim1::Integer, dims::Integer...)
    A = Array{randtype(sets)}(dim1, dims...)
    for i in eachindex(A)
        A[i] = rand(sets)
    end
    A
end
function Base.in{S<:AbstractSet}(xs::AbstractArray, sets::AbstractArray{S})
    size(xs) == size(sets) && all(map(in, xs, sets))
end

"Groups several heterogenous sets. Used mainly for proper dispatch."
immutable TupleSet{T<:Tuple} <: AbstractSet{T}
    sets::T
end
TupleSet(sets::AbstractSet...) = TupleSet(sets)

# rand can return arrays or tuples, but defaults to arrays
randtype(sets::TupleSet, ::Type{Vector}) = Vector{promote_type(map(randtype, sets.sets)...)}
Base.rand(sets::TupleSet, ::Type{Vector}) = eltype(randtype(sets, Vector))[rand(s) for s in sets.sets]
randtype(sets::TupleSet, ::Type{Tuple}) = Tuple{map(randtype, sets.sets)...}
Base.rand(sets::TupleSet, ::Type{Tuple}) = map(rand, sets.sets)
function Base.rand{OT}(sets::TupleSet, ::Type{OT}, dim1::Integer, dims::Integer...)
    A = Array{randtype(sets, OT)}(dim1, dims...)
    for i in eachindex(A)
        A[i] = rand(sets, OT)
    end
    A
end
Base.length(sets::TupleSet) = sum(length(s) for s in sets.sets)
Base.start(sets::TupleSet) = start(sets.sets)
Base.done(sets::TupleSet, i) = done(sets.sets, i)
Base.next(sets::TupleSet, i) = next(sets.sets, i)

randtype(sets::TupleSet) = randtype(sets, Vector)
Base.rand(sets::TupleSet, dims::Integer...) = rand(sets, Vector, dims...)
Base.in(x, sets::TupleSet) = all(map(in, x, sets.sets))

"Returns an AbstractSet representing valid input values"
function inputdomain end

"Returns an AbstractSet representing valid output/target values"
function targetdomain end

export

    # Types
    Cost,
        Loss,
            SupervisedLoss,
                MarginLoss,
                DistanceLoss,
            UnsupervisedLoss,
        Penalty,

    Transformation,
        Learnable,
        StochasticTransformation,

    Minimizable,

    AbstractSet,
        IntervalSet,
        DiscreteSet,
        TupleSet,

    # Functions
    getobs,
    getobs!,

    scaled,

    learn,
    learn!,
    update,
    update!,
    transform,
    transform!,
    value,
    value!,
    meanvalue,
    sumvalue,
    meanderiv,
    sumderiv,
    deriv,
    deriv!,
    params,
    params!,
    grad,
    grad!,
    addgrad!,
    deriv2,
    value_deriv,
    value_deriv!,
    value_grad,
    value_grad!,
    prox,
    prox!,
    inputdomain,
    targetdomain,

    isminimizable,
    isdifferentiable,
    istwicedifferentiable,
    isconvex,
    isstrictlyconvex,
    isstronglyconvex,
    isnemitski,
    isunivfishercons,
    isfishercons,
    islipschitzcont,
    islocallylipschitzcont,
    islipschitzcont_deriv,
    isclipable,
    ismarginbased,
    isclasscalibrated,
    isdistancebased,

    # Base
    issymmetric,

    getobs,
    getobs!,
    datasubset,
    ObsDim,

    targets,

    AbstractDataIterator,
        AbstractObsIterator,
        AbstractBatchIterator,

    DataView,
        AbstractObsView,
        AbstractBatchView,

    DataIterator,
        ObsIterator,
        BatchIterator,

    # StatsBase
    fit,
    fit!,
    nobs

end # module