text:
# This file holds the functions to configure a manpower simulation from an
#   SQLite configuration database.

export configureSimFromDatabase

"""
```
configureSimFromDatabase( mpSim::ManpowerSimulation,
                          dbName::String,
                          overWriteResults:Bool = true;
                          configName::String = "config" )
```
This function configures the manpower simulation `mpSim` form the configuration
table with name `configName` in the SQLite database with filename `dbName`. If
the flag `overWriteResults` is set to `true`, the results of the previous
simulation, if any, will be wiped. If the database is missing the appropriate
extension, `.sqlite`, it will be appended to the name.

The function returns a `Bool`, indicating whether the configuration was
succesful or if the database was corrupted.
"""
function configureSimFromDatabase( mpSim::ManpowerSimulation, dbName::String,
    overwriteResults::Bool = true;
    configName::String = "config" )::Bool

    tmpDBname = dbName * ( endswith( dbName, ".sqlite" ) ? "" : ".sqlite" )

    # Check if the file exists.
    if !ispath( tmpDBname )
        error( "File '$tmpDBname' does not exist." )
    end  # if !ispath( tmpDBname )

    configDB = SQLite.DB( tmpDBname )

    # check if a configuration table exists.
    if configName ∉ SQLite.tables( configDB )[ :name ]
        error( "Database `$tmpDBname` does not contain a manpower simulation configuration." )
    end  # if configName ∉ SQLite.tables( configDB )[ :name ]

    # Get the general parameters
    try
        readGeneralParsFromDatabase( configDB, mpSim, configName )
        readAttributesFromDatabase( configDB, mpSim, configName )
        readAttritionFromDatabase( configDB, mpSim, configName )
        readStatesFromDatabase( configDB, mpSim, configName )
        readTransTypesFromDatabase( configDB, mpSim, configName )
        readTransitionsFromDatabase( configDB, mpSim, configName )
        readRecruitmentFromDatabase( configDB, mpSim, configName )
        readRetirementFromDatabase( configDB, mpSim, configName )

        # Current sim time.
        queryCmd = "SELECT realPar FROM $configName WHERE parName IS 'Sim time'"
        currentTime = SQLite.query( configDB, queryCmd )[ 1, 1 ]

        # Wipe simulation results when requested, or advance the simulation to
        #   time in config file.
        if overwriteResults
            resetSimulation( mpSim )
        elseif currentTime < mpSim.simLength
            run( mpSim.sim, currentTime )
        else
            run( mpSim.sim )
        end  # if currentTime < simDB.simLength

        mpSim.isInitialised = true
        mpSim.isVirgin = now( mpSim ) == 0.0
    catch
        warn( "Configuration not properly read due to database corruption.\n" *
        "Simulation is not well configured and should not be executed." )
        return false
    end  # try

    return true

end  # configureSimFromDatabase( mpSim, dbName, configName )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

"""
```
readGeneralParsFromDatabase( configDB::SQLite.DB,
                             mpSim::ManpowerSimulation,
                             configName::String )
```
This function reads the general simulation parameters from the SQLite database
`configDB`, looking in the table with name `configName`, and uses these
parameters to configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the funuction will
issue warnings, or throw an error depending on the severity.
"""
function readGeneralParsFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'General'"
    generalPars = SQLite.query( configDB, queryCmd )

    # Config file name.
    index = findfirst( generalPars[ :parName ], "Config file" )
    mpSim.parFileName = generalPars[ :strPar ][ index ]

    # Database name.
    index = findfirst( generalPars[ :parName ], "DB name" )
    mpSim.dbName = generalPars[ :strPar ][ index ]
    mpSim.simDB = SQLite.DB( mpSim.dbName )

    # This line ensures that foreign key logic works.
    SQLite.execute!( mpSim.simDB, "PRAGMA foreign_keys = ON" )

    # Catalogue name.
    index = findfirst( generalPars[ :parName ], "Catalogue name" )
    mpSim.catFileName = generalPars[ :strPar ][ index ]

    # Simulation name.
    index = findfirst( generalPars[ :parName ], "Sim name" )
    mpSim.simName = generalPars[ :strPar ][ index ]
    mpSim.personnelDBname = "Personnel_" * mpSim.simName
    mpSim.historyDBname = "History_" * mpSim.simName
    mpSim.transitionDBname = "Transitions_" * mpSim.simName

    # ID key.
    index = findfirst( generalPars[ :parName ], "ID key" )
    mpSim.idKey = generalPars[ :strPar ][ index ]

    # Personnel target.
    index = findfirst( generalPars[ :parName ], "Personnel target" )
    setPersonnelCap( mpSim, generalPars[ :intPar ][ index ] )

    # Start date.
    index = findfirst( generalPars[ :parName ], "Start date" )
    setSimStartDate( mpSim, Date( generalPars[ :strPar ][ index ] ) )

    # Sim length.
    index = findfirst( generalPars[ :parName ], "Sim length" )
    setSimulationLength( mpSim, generalPars[ :realPar ][ index ] )

    # Number of DB commits.
    index = findfirst( generalPars[ :parName ], "DB commits" )
    setDatabaseCommitTime( mpSim, mpSim.simLength /
        generalPars[ :intPar ][ index ] )

    return

end  # readGeneralParsFromDatabase( configDB, mpSim, configName )

"""
```
readAttributesFromDatabase( configDB::SQLite.DB,
                            mpSim::ManpowerSimulation,
                            configName::String )
```
This function reads the personnel attributes from the SQLite database
`configDB`, looking in the table with name `configName`, and uses these
parameters to configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the funuction will
issue warnings, or throw an error depending on the severity.
"""
function readAttributesFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    clearAttributes!( mpSim )
    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Attribute'"
    attributes = SQLite.query( configDB, queryCmd )
    nAttrs = length( attributes[ :parName ] )

    for ii in 1:nAttrs
        newAttr = PersonnelAttribute( strip( attributes[ :parName ][ ii ] ) )

        # Is the attribute fixed?
        setFixed!( newAttr, parse( Bool, attributes[ :boolPar ][ ii ] ) )

        # Read the parameters stored as a string.
        attrPars = split( attributes[ :strPar ][ ii ], ";" )

        # Read the possible values the attribute can take.
        if attrPars[ 1 ] != "[]"
            possibleVals = map( val -> String( strip( val ) ),
                split( attrPars[ 1 ][ 2:(end - 1) ], "," ) )
            setPossibleValues!( newAttr, possibleVals )
        end  # if if attrPars[ 1 ] != "[]"

        # Read the initial values of the attribute with their respective
        #   generation probabilities.
        if attrPars[ 2 ] != "[]"
            initVals = split( attrPars[ 2 ][ 2:(end - 1) ], "," )
            initValDict = Dict{String, Float64}()

            for valPair in initVals
                valPair = split( valPair, ":" )
                initValDict[ strip( valPair[ 1 ] ) ] = parse( Float64,
                    valPair[ 2 ] )
            end  # for valPair in initVals

            setAttrValues!( newAttr, initValDict )
        end  # if attrPars[ 2 ] != "[]"

        addAttribute!( mpSim, newAttr )
    end  # for ii in 1:nAttrs

    return

end  # readAttributesFromDatabase( configDB, mpSim, configName )

"""
```
readAttritionFromDatabase( configDB::SQLite.DB,
                           mpSim::ManpowerSimulation,
                           configName::String )
```
This function reads the retirement parameters from the SQLite database
`configDB`, looking in the table with name `configName`, and uses these
parameters to configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the funuction will
issue warnings, or throw an error depending on the severity.
"""
function readAttritionFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    clearAttritionSchemes!( mpSim )
    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Attrition'"
    attrList = SQLite.query( configDB, queryCmd )
    nSchemes = length( attrList[ :parName ] )
    isDefaultFound = false

    # If no attrition is in the database, set the attrition scheme to nothing
    #   and finish.
    if nSchemes == 0
        setAttrition( mpSim )
        return
    end  # if length( retirement ) == 0

    # Process attrition schemes.
    for ii in 1:nSchemes
        attrName = strip( attrList[ :parName ][ ii ] )
        newAttrScheme = Attrition( attrName )

        # Read attrition period.
        setAttritionPeriod( newAttrScheme, attrList[ :realPar ][ ii ] )

        # Read attrition curve.
        attrCurve = split( attrList[ :strPar ][ ii ], ";" )
        attrCurveDict = Dict{Float64, Float64}()

        for attrPair in attrCurve
            attrPair = split( attrPair, ":" )
            attrCurveDict[ parse( Float64, attrPair[ 1 ] ) ] = parse( Float64,
                attrPair[ 2 ] )
        end  # for attrPair in attrCurve

        setAttritionCurve( newAttrScheme, attrCurveDict )

        if attrName == "default"
            setAttrition( mpSim, newAttrScheme )
            isDefaultFound = true
        else
            addAttritionScheme!( newAttrScheme, mpSim )
        end  # if attrName == "default"
    end  # for ii in 1:nSchemes

    # Add a default attrition rate of 0% if there's no default in the database.
    if !isDefaultFound
        setAttrition( mpSim )
    end

    return

end  # readAttritionFromDatabase( configDB, mpSim, configName )

"""
```
readStatesFromDatabase( configDB::SQLite.DB,
                        mpSim::ManpowerSimulation,
                        configName::String )
```
This function reads the personnel states from the SQLite database `configDB`,
looking in the table with name `configName`, and uses these parameters to
configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the function will
issue warnings, or throw an error depending on the severity.
"""
function readStatesFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    clearStates!( mpSim )
    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'State'"
    states = SQLite.query( configDB, queryCmd )
    nStates = length( states[ :parName ] )

    for ii in 1:nStates
        newState = State( strip( states[ :parName ][ ii ] ) )
        setInitial!( newState, parse( Bool, states[ :boolPar ][ ii ] ) )
        setStateTarget!( newState, states[ :intPar ][ ii ] )
        setStateRetirementAge!( newState, states[ :realPar ][ ii ] )

        # Read all the parameters stored as string.
        statePars = split( states[ :strPar ][ ii ], ";" )

        # Read the state requirements.
        if statePars[ 1 ] != "[]"
            stateReqs = split( statePars[ 1 ][ 2:(end - 1) ], "," )

            for stateReq in stateReqs
                stateReq = split( stateReq, ":" )
                addRequirement!( newState, String( strip( stateReq[ 1 ] ) ),
                    map( val -> String( strip( val ) ),
                    split( stateReq[ 2 ], "//" ) ) )
            end  # for stateReq in stateReqs
        end  # if statePars[ 1 ] != "[]"

        # Read the state attrition regime.
        stateAttr = split( statePars[ 2 ], ":" )

        if length( stateAttr ) == 1
            schemeName = strip( stateAttr[ 1 ] )
            attrScheme = schemeName == "default" ?
                mpSim.defaultAttritionScheme :
                mpSim.attritionSchemes[ schemeName ]
            setStateAttritionScheme!( newState, attrScheme )
        elseif stateAttr[ 1 ] == "Attrition"
            # This one is necessary because generated attrition schemes start
            #   with "Attrition:".
            schemeName = join( stateAttr, ":" )
            attrScheme = mpSim.attritionSchemes[ schemeName ]
            setStateAttritionScheme!( newState, attrScheme )
        else
            setStateAttritionScheme!( newState,
                parse( Float64, stateAttr[ 2 ] ),
                parse( Float64, stateAttr[ 1 ] ), mpSim )
        end  # if length( stateAttr ) == 1

        addState!( mpSim, newState, newState.isInitial )
    end  # for ii in 1:nAttrs

    return

end  # readStatesFromDatabase( configDB, mpSim, configName )

function readCompoundstatesFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    clearCompoundStates!( mpSim )
    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Compound state'"
    compStates = SQLite.query( configDB, queryCmd )
    nCompStates = length( compStates[ :parName ] )

    for ii in 1:nCompStates
        isStateList = parse( Bool, compStates[ :boolPar ][ ii ] )
        stateName = compStates[ :parName ][ ii ]
        newCompState = isStateList ? CompoundState( stateName ) :
            State( stateName )
        setStateTarget!( newCompState, compStates[ :intPar ][ ii ] )
        strParProc = String.( split( compStates[ :strPar ][ ii ], ";" ) )

        if isStateList
            addStateToCompound!( newCompState, strParProc... )
            addCompoundState!( mpSim, newCompState )
        else
            strParProc = map( ii -> String.( split( strParProc[ ii ], ":" ) ),
                eachindex( strParProc ) )
            foreach( attrPair -> addRequirement!( newCompState, attrPair[ 1 ],
                attrPair[ 2 ] ), strParProc )
            mpSim.compoundStates[ stateName ] = newCompState
        end  # if isStateList
    end  # for ii in 1:nCompStates

    return

end  # readCompoundstatesFromDatabase( configDB, mpSim, configName )

"""
```
readTransTypesFromDatabase( configDB::SQLite.DB,
                            mpSim::ManpowerSimulation,
                            configName::String )
```
This function reads the transition types from the SQLite database `configDB`,
looking in the table with name `configName`, and uses these parameters to
configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the function will
issue warnings, or throw an error depending on the severity.
"""
function readTransTypesFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    clearTransitions!( mpSim )
    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Transition type'"
    transTypes = SQLite.query( configDB, queryCmd )
    nTypes = length( transTypes[ :parName ] )
    foreach( 1:nTypes ) do ii
        setTransTypePriority!( mpSim, transTypes[ :parName ][ ii ],
            transTypes[ :intPar ][ ii ] )
    end  # foreach( 1:nTypes ) do ii

    return

end  # readTransTypesFromDatabase( configDB, mpSim, configName )

"""
```
readTransitionsFromDatabase( configDB::SQLite.DB,
                             mpSim::ManpowerSimulation,
                             configName::String )
```
This function reads the state transitions from the SQLite database `configDB`,
looking in the table with name `configName`, and uses these parameters to
configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the function will
issue warnings, or throw an error depending on the severity.
"""
function readTransitionsFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Transition'"
    transitions = SQLite.query( configDB, queryCmd )
    nTrans = length( transitions[ :parName ] )
    stateList = collect( keys( merge( mpSim.initStateList,
        mpSim.otherStateList ) ) )

    for ii in 1:nTrans
        transPars = split( transitions[ :strPar ][ ii ], ";" )
        newTrans = Transition( strip( transitions[ :parName ][ ii ] ),
            mpSim.stateList[ strip( transPars[ 1 ] ) ],
            mpSim.stateList[ strip( transPars[ 2 ] ) ] )
        setSchedule( newTrans, parse( Float64, transPars[ 3 ] ),
            parse( Float64, transPars[ 4 ] ) )

        # Read the extra conditions on the transition.
        if transPars[ 5 ] != "[]"
            conds = split( transPars[ 5 ][ 2:(end - 1) ], "," )

            for cond in conds
                condParts = split( cond, ":" )
                condVal = tryparse( Float64, condParts[ 3 ] )

                if condVal.hasvalue
                    condVal = condVal.value / 12.0
                else
                    condVal = replace( strip( String( condParts[ 3 ] ) ), "//",
                        "," )
                end  # if condVal.hasvalue

                newCond = processCondition( strip( String( condParts[ 1 ] ) ),
                    strip( String( condParts[ 2 ] ) ), condVal )
                addCondition!( newTrans, newCond[ 1 ] )
            end  # for cond in conds
        end  # if transPars[ 5 ] != "[]"

        # Read the extra changes of the transition.
        if transPars[ 6 ] != "[]"
            changes = split( transPars[ 6 ][ 2:(end - 1) ], "," )

            for change in changes
                change = split( change, ":" )
                addAttributeChange!( newTrans, String( change[ 1 ] ),
                    String( change[ 2 ] ) )
            end  # for change in changes
        end  # if transPars[ 6 ] != "[]"

        # Read the final parameters.
        probList = map( prob -> parse( Float64, prob ),
            split( transPars[ 7 ][ 2:(end - 1) ], "," ) )
        setTransProbabilities( newTrans, probList )
        setMaxAttempts( newTrans, parse( Int, transPars[ 8 ] ) )
        setFireAfterFail( newTrans, parse( Bool, transPars[ 9 ] ) )
        setMaxFlux( newTrans, parse( Int, transPars[ 10 ] ) )
        setHasPriority( newTrans, parse( Bool, transPars[ 11 ] ) )
        addTransition!( mpSim, newTrans )
    end  # for ii in 1:nTrans

    return

end  # readTransitionsFromDatabase( configDB, mpSim, configName )

"""
```
readRecruitmentFromDatabase( configDB::SQLite.DB,
                             mpSim::ManpowerSimulation,
                             configName::String )
```
This function reads the recruitment schemes from the SQLite database
`configDB`, looking in the table with name `configName`, and uses these
parameters to configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the function will
issue warnings, or throw an error depending on the severity.
"""
function readRecruitmentFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    clearRecruitmentSchemes!( mpSim )
    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Recruitment'"
    recruitment = SQLite.query( configDB, queryCmd )
    numSchemes = length( recruitment[ :parName ] )

    # If no recruitment schemes are in the database, do nothing.
    if numSchemes == 0
        return
    end  # if numSchemes == 0

    for ii in 1:numSchemes
        recScheme = Recruitment( recruitment[ :parName ][ ii ], 1 )
        recPars = split( recruitment[ :strPar ][ ii ], ";" )
        setRecruitmentSchedule( recScheme, parse( Float64, recPars[ 1 ] ),
            parse( Float64, recPars[ 2 ] ) )
        setRecruitState( recScheme, String( strip( recPars[ 3 ] ) ) )

        # Read the distribution of the number of people to recruit.
        if parse( Bool, recruitment[ :boolPar ][ ii ] )
            setRecruitmentLimits( recScheme, parse( Int, recPars[ 4 ] ),
                parse( Int, recPars[ 5 ] ) )
        else
            distType = Symbol( strip( recPars[ 4 ] ) )
            distNodes = split( recPars[ 5 ][ 2:(end - 1) ], "," )
            distNodesDict = Dict{Int, Float64}()

            for node in distNodes
                node = split( node, ":" )
                distNodesDict[ parse( Float64, node[ 1 ] ) ] =
                    parse( Float64, node[ 2 ] )
            end  # for node in distNodes

            setRecruitmentDistribution( recScheme, distNodesDict, distType )
        end  # if isAdaptive

        ageDistType = Symbol( strip( recPars[ 6 ] ) )
        ageDistNodes = split( recPars[ 7 ][ 2:(end - 1) ], "," )
        ageDistNodesDict = Dict{Float64, Float64}()

        for node in ageDistNodes
            node = split( node, ":" )
            ageDistNodesDict[ parse( Float64, node[ 1 ] ) ] =
                parse( Float64, node[ 2 ] )
        end  # for node in distNodes

        setAgeDistribution( recScheme, ageDistNodesDict, ageDistType )
        addRecruitmentScheme!( mpSim, recScheme )
    end  # for ii in 1:numSchemes

    return

end  # readRecruitmentFromDatabase( configDB, mpSim, configName )

"""
```
readRetirementFromDatabase( configDB::SQLite.DB,
                            mpSim::ManpowerSimulation,
                            configName::String )
```
This function reads the retirement parameters from the SQLite database
`configDB`, looking in the table with name `configName`, and uses these
parameters to configure the manpower simulation `mpSim`.

This function returns `nothing`. If information is missing, the funuction will
issue warnings, or throw an error depending on the severity.
"""
function readRetirementFromDatabase( configDB::SQLite.DB,
    mpSim::ManpowerSimulation, configName::String )::Void

    queryCmd = "SELECT * FROM $configName
        WHERE parType IS 'Retirement' AND parName IS 'Retirement'"
    retirement = SQLite.query( configDB, queryCmd )

    # If no retirement is in the database, set the retirement scheme to nothing
    #   and finish.
    if length( retirement ) == 0
        setRetirement( mpSim )
        return
    end  # if length( retirement ) == 0

    retScheme = Retirement( isEither = parse( Bool,
        retirement[ :boolPar ][ 1 ] ) )

    retPars = split( retirement[ :strPar ][ 1 ], ";" )
    setCareerLength( retScheme, parse( Float64, retPars[ 1 ] ) )
    setRetirementAge( retScheme, parse( Float64, retPars[ 2 ] ) )
    setRetirementSchedule( retScheme, parse( Float64, retPars[ 3 ] ),
        parse( Float64, retPars[ 4 ] ) )
    setRetirement( mpSim, retScheme )

    return

end  # readRetirementFromDatabase( configDB, mpSim, configName )


text:
for c in [',', ' ', ':', '.', '!', '?', '\n', '\'', '’']
    DICASCII[c] = missing
end
DICASCII['É'] = 'e';
DICASCII['Â'] = 'a';
DICASCII['À'] = 'a';

function nettoie(s)
    sn = Array{UInt8, 1}()
    for c in s
        cn = DICASCII[c]
        if !ismissing(cn)
            push!(sn, cn)
        end
    end
    sn
end

function control(AMOTS, lir, sn)
    if length(sn) == 0
        # println(lir)
	return (true, lir)
    end
    for i in 1:length(AMOTS)
        if length(lir) > 0 && i == lir[end]
	    continue
	end
        if sn[1] in AMOTS[i]
            nlir = [lir; i]
	    ret = control(AMOTS, nlir, sn[2:end])
	    if ret[1]
	        return ret
            end
        end
    end
    return false
end

control(AMOTS::Array{Array{UInt8,1},1},
        sn::Array{UInt8,1}) = control(AMOTS, Array{UInt8, 1}(), sn)
control(amots::String, s::String) = control(liste_mots(amots), nettoie(s))

ext: jl


text:
using Test
using JSONQ
using JSON3

@testset "JSONQ" begin
    maths_test_dataset_file = joinpath(@__DIR__, "sample_math_test_results.json")
    maths_test_dataset_url = "https://data.cityofnewyork.us/api/views/7yig-nj52/rows.json?accessType=DOWNLOAD"
    isfile(maths_test_dataset_file) || download(maths_test_dataset_url, maths_test_dataset_file)
    maths_test_dataset = JSON3.read(read(maths_test_dataset_file, String))

    ctx = JSONQ.Ctx()
    JSONQ.compile(ctx, ".meta.view.attributionLink"; name="attr")
    JSONQ.compile(ctx, ".meta.view.columns[1].name"; name="col1name")

    @test JSONQ.execute(ctx, maths_test_dataset, "attr") == "http://schools.nyc.gov/NR/rdonlyres/3E439D63-A14E-4B6E-B019-B634A8915D71/0/DistrictMathResults20062012Public.xlsx"
    @test JSONQ.execute(ctx, maths_test_dataset, "col1name") == "sid"
end


text:
@doc """
    merge!(S::SeisData, U::SeisData[, prune_only=true])

Merge channels of two SeisData structures.

    merge!(S::SeisData[, prune_only=true])

"Flatten" a SeisData structure by merging channels with identical properties.

If `prune_only=true`, the only action taken is deletion of empty and duplicate
channels; `merge!(S, U, prune_only=true)` is identical to an in-place `S+U`.
""" merge!
function merge!(S::Y; v::Integer=KW.v, purge_only::Bool=false) where Y<:GphysData
  # Required preprocessing
  prune!(S)

  # Initialize variables
  ID    = getfield(S, :id)
  NAME  = getfield(S, :name)
  LOC   = getfield(S, :loc)
  FS    = getfield(S, :fs)
  GAIN  = getfield(S, :gain)
  RESP  = getfield(S, :resp)
  UNITS = getfield(S, :units)
  SRC   = getfield(S, :src)
  MISC  = getfield(S, :misc)
  NOTES = getfield(S, :notes)
  T     = getfield(S, :t)
  X     = getfield(S, :x)

  UID = unique(getfield(S, :id))
  cnt = 0
  note_head = string(SeisIO.timestamp(), " ¦ ")
  to_delete = Array{Int64,1}(undef, 0)
  while cnt < length(UID)
    cnt = cnt + 1
    id = getindex(UID, cnt)
    GRP = findall(S.id.==id)
    SUBGRPS = get_subgroups(LOC, FS, RESP, UNITS, GRP)

    for subgrp in SUBGRPS
      dup_check!(subgrp, to_delete, T, X)
      (purge_only == true) && continue

      N     = length(subgrp)
      i1    = getindex(subgrp, 1)
      fs    = getindex(FS, i1)
      if fs == 0.0
        Ω = merge_non_ts!(S, subgrp)
        append!(to_delete, subgrp[subgrp.!=Ω])
        continue
      end
      Δ     = round(Int64, sμ/fs)
      W, Ω  = get_merge_w(Δ, subgrp, T, X)
      rest  = subgrp[subgrp.!=Ω]

      # GAIN, MISC, NAME, NOTES, SRC ==========================================
      gain  = getindex(GAIN, Ω)
      notes = getindex(NOTES, Ω)
      misc  = getindex(MISC, Ω)
      for i in rest
        scalefac = gain / getindex(GAIN, i)
        if scalefac != 1.0
          rmul!(getindex(X, i), scalefac)
          GAIN[i] = gain            # added 2020-12-03, in case merge breaks
        end
        if getindex(SRC, i) != getindex(SRC, Ω)
          push!(notes, string(note_head, "+source ¦ ", getindex(SRC, i)))
        end
        if getindex(NAME, i) != getindex(NAME, Ω)
          push!(notes, string(note_head, "alt name ¦ ", getindex(NAME, i)))
        end
        append!(notes, getindex(NOTES, i))
        merge!(misc, getindex(MISC, i))
      end
      sort!(notes)

      # Extra fields ==========================================================
      merge_ext!(S, Ω, rest)

      # T,X ===================================================================
      T_Ω, X_Ω = segment_merge(Δ, Ω, W, X)
      setindex!(X, X_Ω, Ω)
      setindex!(T, T_Ω, Ω)

      # Document it
      proc_note!(S, Ω, "merge!", string("combined channels ",
                 repr(subgrp), " (N = ", N, ") into :t, :x"))
      append!(to_delete, rest)
    end
    # Done with SUBGRPS
  end
  deleteat!(S, to_delete)
  sort!(S)
  return nothing
end

# The following must be the same (or unset) to merge:
# :loc
# :resp
# :fs         cannot be unset; not a nullable array
# :units

# The following are easily dealt with:
# :gain       can be translated with no trouble
# :misc       can call merge! method for dictionaries
# :notes      append and sort
# :name       not important, can log extra names to :notes
# :src        use most recent, log extras to :notes

function merge!(C::T1, D::T2) where {T1<:GphysChannel, T2<:GphysChannel}

  # Identical structures or empty D
  ((C == D) || isempty(D.x) || isempty(D.t)) && return nothing

  # empty C
  if (isempty(C.x) || isempty(C.t))
    ff = (T1 == T2) ? fieldnames(T1) : SeisIO.datafields
    for f in ff
      setfield!(C, f, deepcopy(getfield(D, f)))
    end
    return nothing
  end

  # partial match or exit
  m = cmatch_p!(C, D)
  (m == false) && (@warn("Critical field mismatch! Not merged!"); return nothing)

  # auto-convert T2 to T1
  if T2 != T1
    D = convert(T1, D)
  end

  note_head = string(SeisIO.timestamp(), " ¦ ")
  ttest = (endtime(C.t, C.fs) > endtime(D.t, D.fs))
  ω = ttest ? C : D
  α = ttest ? D : C

  # at this point, (:fs, :gain, :loc, :resp, :units) are known to match
  # this is stricter than merge!, which allows gain mismatches within S

  # MISC, NAME, NOTES, SRC ================================================
  if ω.src != α.src
    push!(ω.notes, string(note_head, "+source ¦ ", α.src))
  end
  if ω.name != α.name
    push!(ω.notes, string(note_head, "alt name ¦ ", α.name))
  end

  append!(C.notes, D.notes)
  sort!(C.notes)
  if ttest
    merge!(C.misc, D.misc)
  else
    C.misc = merge(D.misc, C.misc)
  end

  # Extra fields ==========================================================
  merge_ext!(ω, α)

  # T,X ===================================================================
  if C.fs == 0.0
    merge_non_ts!(C, D)
    return nothing
  end

  # Proceed only for time-series data
  subgrp = ttest ? [1, 2] : [2, 1]
  T = [ω.t, α.t]
  X = Array{FloatArray, 1}(undef, 2)
  X[1] = ω.x
  X[2] = α.x
  dup_check!(subgrp, Int64[], T, X)
  N = length(subgrp)

  # At this point, N == 1 means the same data at the same times; nothing to do
  (N == 1) && return nothing

  # on to time windows
  fs          = C.fs
  Δ           = round(Int64, sμ/fs)
  W, Ω        = get_merge_w(Δ, subgrp, T, X)
  TΩ, XΩ      = segment_merge(Δ, Ω, W, X)
  C.t = TΩ
  C.x = XΩ
  proc_note!(C, "merge!", string("combined data from a structure of type $T2"))
  return nothing
end


text:
"""
    struct PSpace1D{TR,TI<:Integer,TA<:AbstractArray} <: AbstractPhysicalSpace
        x0::TR
        x1::TR
        nx::TI
        x::TA
        dx::TA
    end

1D physical space with structured mesh

"""
struct PSpace1D{TR,TI<:Integer,TA<:AbstractArray} <: AbstractPhysicalSpace1D
    x0::TR
    x1::TR
    nx::TI
    x::TA
    dx::TA
end

function PSpace1D(X0::TR, X1::TR, NX::TI, NG = 0::Integer) where {TR,TI}
    TX = ifelse(TR == Float32, Float32, Float64)

    δ = (X1 - X0) / NX
    x = OffsetArray{TX}(undef, 1-NG:NX+NG)
    dx = similar(x)

    # uniform mesh
    for i in eachindex(x)
        x[i] = X0 + (i - 0.5) * δ
        dx[i] = δ
    end

    return PSpace1D{TR,TI,typeof(x)}(X0, X1, NX, x, dx)
end

PSpace1D() = PSpace1D(0, 1, 100)
PSpace1D(X0::T, X1::T) where {T} = PSpace1D(X0, X1, 100)

"""
    struct PSpace2D{TR<:Real,TI<:Integer,TA<:AbstractMatrix{<:Real},TB<:AbstractArray{<:Real,4}} <: AbstractPhysicalSpace2D
        x0::TR
        x1::TR
        nx::TI
        y0::TR
        y1::TR
        ny::TI
        x::TA
        y::TA
        dx::TA
        dy::TA
        vertices::TB
    end

2D Physical space with structured mesh

"""
struct PSpace2D{
    TR<:Real,
    TI<:Integer,
    TA<:AbstractMatrix{<:Real},
    TB<:AbstractArray{<:Real,4},
} <: AbstractPhysicalSpace2D
    x0::TR
    x1::TR
    nx::TI
    y0::TR
    y1::TR
    ny::TI
    x::TA
    y::TA
    dx::TA
    dy::TA
    vertices::TB
end

function PSpace2D(
    X0::TR,
    X1::TR,
    NX::TI,
    Y0::TR,
    Y1::TR,
    NY::TI,
    NGX = 0::Integer,
    NGY = 0::Integer,
) where {TR,TI}
    TX = ifelse(TR == Float32, Float32, Float64)

    δx = (X1 - X0) / NX
    δy = (Y1 - Y0) / NY
    x = OffsetArray{TX}(undef, 1-NGX:NX+NGX, 1-NGY:NY+NGY)
    y = OffsetArray{TX}(undef, 1-NGX:NX+NGX, 1-NGY:NY+NGY)
    dx = OffsetArray{TX}(undef, 1-NGX:NX+NGX, 1-NGY:NY+NGY)
    dy = OffsetArray{TX}(undef, 1-NGX:NX+NGX, 1-NGY:NY+NGY)

    for j in axes(x, 2)
        for i in axes(x, 1)
            x[i, j] = X0 + (i - 0.5) * δx
            y[i, j] = Y0 + (j - 0.5) * δy
            dx[i, j] = δx
            dy[i, j] = δy
        end
    end

    vertices = OffsetArray{TX}(undef, 1-NGX:NX+NGX, 1-NGY:NY+NGY, 4, 2)
    for j in axes(vertices, 2), i in axes(vertices, 1)
        vertices[i, j, 1, 1] = x[i, j] - 0.5 * dx[i, j]
        vertices[i, j, 2, 1] = x[i, j] + 0.5 * dx[i, j]
        vertices[i, j, 3, 1] = x[i, j] + 0.5 * dx[i, j]
        vertices[i, j, 4, 1] = x[i, j] - 0.5 * dx[i, j]

        vertices[i, j, 1, 2] = y[i, j] - 0.5 * dy[i, j]
        vertices[i, j, 2, 2] = y[i, j] - 0.5 * dy[i, j]
        vertices[i, j, 3, 2] = y[i, j] + 0.5 * dy[i, j]
        vertices[i, j, 4, 2] = y[i, j] + 0.5 * dy[i, j]
    end

    return PSpace2D{TR,TI,typeof(x),typeof(vertices)}(
        X0,
        X1,
        NX,
        Y0,
        Y1,
        NY,
        x,
        y,
        dx,
        dy,
        vertices,
    )
end

PSpace2D() = PSpace2D(0, 1, 45, 0, 1, 45)
PSpace2D(X0::T, X1::T, Y0::T, Y1::T) where {T} = PSpace2D(X0, X1, 45, Y0, Y1, 45)

"""

2D Circular space in polar coordinates

"""
struct CSpace2D{
    TR<:Real,
    TI<:Integer,
    TA<:AbstractMatrix{<:Real},
    TB<:AbstractArray{<:Real,4},
} <: AbstractPhysicalSpace2D
    r0::TR
    r1::TR
    nr::TI
    θ0::TR
    θ1::TR
    nθ::TI
    r::TA
    θ::TA
    x::TA
    y::TA
    dr::TA
    dθ::TA
    darc::TA
    vertices::TB
end

function CSpace2D(
    R0::TR,
    R1::TR,
    NR::TI,
    θ0,
    θ1,
    Nθ::TI,
    NGR = 0::TI,
    NGθ = 0::TI,
) where {TR,TI<:Integer}
    TX = ifelse(TR == Float32, Float32, Float64)

    δr = (R1 - R0) / NR
    δθ = (θ1 - θ0) / Nθ
    r = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)
    θ = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)
    x = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)
    y = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)
    dr = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)
    dθ = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)
    darc = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ)

    for j in axes(r, 2)
        for i in axes(r, 1)
            r[i, j] = R0 + (i - 0.5) * δr
            θ[i, j] = θ0 + (j - 0.5) * δθ
            x[i, j] = r[i, j] * cos(θ[i, j])
            y[i, j] = r[i, j] * sin(θ[i, j])
            dr[i, j] = δr
            dθ[i, j] = δθ
            darc[i, j] = r[i, j] * dθ[i, j]
        end
    end

    vertices = OffsetArray{TX}(undef, 1-NGR:NR+NGR, 1-NGθ:Nθ+NGθ, 4, 2)
    for j in axes(vertices, 2), i in axes(vertices, 1)
        vertices[i, j, 1, 1] = (r[i, j] - 0.5 * dr[i, j]) * cos(θ[i, j] - 0.5 * dθ[i, j])
        vertices[i, j, 2, 1] = (r[i, j] + 0.5 * dr[i, j]) * cos(θ[i, j] - 0.5 * dθ[i, j])
        vertices[i, j, 3, 1] = (r[i, j] + 0.5 * dr[i, j]) * cos(θ[i, j] + 0.5 * dθ[i, j])
        vertices[i, j, 4, 1] = (r[i, j] - 0.5 * dr[i, j]) * cos(θ[i, j] + 0.5 * dθ[i, j])

        vertices[i, j, 1, 2] = (r[i, j] - 0.5 * dr[i, j]) * sin(θ[i, j] - 0.5 * dθ[i, j])
        vertices[i, j, 2, 2] = (r[i, j] + 0.5 * dr[i, j]) * sin(θ[i, j] - 0.5 * dθ[i, j])
        vertices[i, j, 3, 2] = (r[i, j] + 0.5 * dr[i, j]) * sin(θ[i, j] + 0.5 * dθ[i, j])
        vertices[i, j, 4, 2] = (r[i, j] - 0.5 * dr[i, j]) * sin(θ[i, j] + 0.5 * dθ[i, j])
    end

    return CSpace2D{TR,TI,typeof(x),typeof(vertices)}(
        R0,
        R1,
        NR,
        TX(θ0),
        TX(θ1),
        Nθ,
        r,
        θ,
        x,
        y,
        dr,
        dθ,
        darc,
        vertices,
    )
end

"""
    uniform_mesh(x0::Real, xnum::Int, dx::Real)

Generate uniform mesh
"""
function uniform_mesh(x0, xnum::T, dx) where {T<:Int}
    points = zeros(xnum)
    for i = 1:xnum
        points[i] = x0 + (i - 0.5) * dx
    end

    return points
end

"""
    ndgrid(v::AbstractVector)
    ndgrid(v1::AbstractVector{T}, v2::AbstractVector{T}) where T
    ndgrid(vs::AbstractVector{T}...) where T

Equivalent N-dimensional mesh generator as matlab
"""
ndgrid(v::AbstractVector) = copy(v)

function ndgrid(v1::AbstractVector{T}, v2::AbstractVector{T}) where {T}
    m, n = length(v1), length(v2)
    v1 = reshape(v1, m, 1)
    v2 = reshape(v2, 1, n)

    return (repeat(v1, 1, n), repeat(v2, m, 1))
end

function ndgrid(vs::AbstractVector{T}...) where {T}
    ndgrid_fill(a, v, s, snext) = begin
        for j = 1:length(a)
            a[j] = v[div(rem(j - 1, snext), s)+1]
        end
    end

    n = length(vs)
    sz = map(length, vs)
    out = ntuple(i -> Array{T}(undef, sz), n)
    s = 1
    for i = 1:n
        a = out[i]::Array
        v = vs[i]
        snext = s * size(a, i)
        ndgrid_fill(a, v, s, snext)
        s = snext
    end

    return out
end

"""
    2D: meshgrid(x::AbstractArray{<:Real,1}, y::AbstractArray{<:Real,1})
    3D: meshgrid(x::AbstractArray{<:Real,1}, y::AbstractArray{<:Real,1}, z::AbstractArray{<:Real,1})

Equivalent structured mesh generator as matlab
"""
meshgrid(v::AbstractVector{T}) where {T} = meshgrid(v, v)

function meshgrid(x::T, y::T) where {T<:AbstractVector{<:Real}}
    X = [i for j in y, i in x]
    Y = [j for j in y, i in x]

    return X, Y
end

function meshgrid(x::T, y::T, z::T) where {T<:AbstractVector{<:Real}}
    X = [i for k in z, j in y, i in x]
    Y = [j for k in z, j in y, i in x]
    Z = [k for k in z, j in y, i in x]

    return X, Y, Z
end

"""
    find_idx(x::AbstractArray{<:Real,1}, p::Real; mode = :nonuniform::Symbol)

Find the location index of a point in mesh

- @arg x: center locations of 1D mesh
- @arg p: point location
- @arg mode: choose uniform / nonuniform formulations

"""
function find_idx(
    x::T,
    p::Real;
    mode = :nonuniform::Symbol,
) where {T<:AbstractArray{<:Real,1}}

    if mode == :uniform
        dx = x[2] - x[1]
        return Int(ceil((p - x[1] + 0.5 * dx) / dx)) # point location
    else
        return argmin(abs.(x .- p)) # center location
    end

end

# ------------------------------------------------------------
# Extended Base.show()
# ------------------------------------------------------------

function Base.show(io::IO, ps::PSpace1D{TR,TI,TA}) where {TR,TI,TA}
    print(
        io,
        "PhysicalSpace1D{$TR,$TI,$TA}\n",
        "domain: ($(ps.x0),$(ps.x1))\n",
        "resolution: $(ps.nx)\n",
        "ghost: $(1-firstindex(ps.x))\n",
    )
end

function Base.show(io::IO, ps::PSpace2D{TR,TI,TA}) where {TR,TI,TA}
    print(
        io,
        "PhysicalSpace2D{$TR,$TI,$TA}\n",
        "domain: ($(ps.x0),$(ps.x1)) × ($(ps.y0),$(ps.y1))\n",
        "resolution: $(ps.nx) × $(ps.ny)\n",
        "ghost in x: $(1-firstindex(ps.x[:, 1]))\n",
        "ghost in y: $(1-firstindex(ps.y[1, :]))\n",
    )
end

function Base.show(io::IO, ps::CSpace2D{TR,TI,TA}) where {TR,TI,TA}
    print(
        io,
        "CircularSpace2D{$TR,$TI,$TA}\n",
        "domain: ($(ps.r0),$(ps.r1)) × ($(ps.θ0),$(ps.θ1))\n",
        "resolution: $(ps.nr) × $(ps.nθ)\n",
        "ghost in r: $(1-firstindex(ps.x[:, 1]))\n",
        "ghost in θ: $(1-firstindex(ps.y[1, :]))\n",
    )
end


text:
## -----
## BigFloat numbers are used; The MutableArithmetics.jl package was
## mimicked below, though not included directly, so rounding modes
## could be specified.  This implementation below still allocates much
## more than is desirable.

# for more generic uses
Big(a;precision=53) = BigFloat(a; precision=precision)

#import Base: precision
precision(x::Float64) = Int(53)
precision(x::BigFloat) = Int(Base.MPFR.precision(x))

## ---- bigfloat
## We add some methods from libmpfr to give rounding
## hints. (cf. MutableArithmetics.jl for inspiration)

# This caused errors which is too bad
# as it means many more allocations because we must create new
# BigFloat values every time the precision is changes
set_prec!(x::BigFloat,n::Int) = ccall(
    (:mpfr_set_prec, :libmpfr),
    Cvoid,
    (Ref{BigFloat}, Cint, Base.MPFR.MPFRRoundingMode),
    x,Cint(n), Base.MPFR.ROUNDING_MODE[]
)

# this is copy!(y, x)
# we can not use `set_prec` to set precision though
set!(y::BigFloat, x, rnd=Base.MPFR.ROUNDING_MODE[]) = ccall(
    (:mpfr_set, :libmpfr),
    Int32,
    (Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
    y, x, rnd
)
zero!(a::BigFloat) = set!(a, 0)

swap!(y::BigFloat,x::BigFloat) = ccall(
    (:mpfr_swap, :libmpfr),
    Cvoid,
    (Ref{BigFloat}, Ref{BigFloat}),
    y, x
)

# in place +,-,*./
function radd!(a::BigFloat, b::BigFloat, rnd=Base.MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_add, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
        a,a,b, rnd
    )
    a
end

function rsub!(a::BigFloat, b::BigFloat, rnd=Base.MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_sub, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
        a,a,b, rnd
    )
    a
end

function LinearAlgebra.rmul!(a::BigFloat, b::Number, rnd=Base.MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_mul, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
        a,a,b, rnd
    )
    a
end

function LinearAlgebra.rdiv!(a::BigFloat, b::BigFloat, rnd=Base.MPFR.ROUNDING_MODE[])
    ccall(
        (:mpfr_div, :libmpfr),
        Int32,
        (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Base.MPFR.MPFRRoundingMode),
        a,a,b, rnd
    )
    a
end

# a <- a*b+c
function muladd!(a::T,b,c) where {T <: BigFloat}
    rmul!(a,b)
    radd!(a,c)
end

# nonallocating online minimum and maximum
function min!(x::T,a::T) where {T <: BigFloat}
    a < x && set!(x, a)
    nothing
end

function max!(x::T,a::T) where {T <: BigFloat}
    a > x && set!(x, a)
    nothing
end

epsᴸ(L::Int) = setprecision(() -> eps(BigFloat), L)

# evalpoly with interval arithmetic to precision L
function evalpolyᴸ(x, p, L)
    Σ = ball(0, L)
    evalpolyᴸ!(Σ, x, p, L)
    return Σ
end

function evalpolyᴸ!(Σ, x, p, L)
    ϵ = epsᴸ(L)
    𝐱 = ball(x, ϵ)
    𝐩ᵢ = ball(L)

    w,d,u = BigFloat(;precision=L),BigFloat(;precision=L),BigFloat(;precision=L)

    evalpolyᴸ!(Σ, 𝐱, p, L, 𝐩ᵢ, w, d, u)
end

function evalpolyᴸ!(Σ, 𝐱,  p, L, # storage, fat x, p, L
                    𝐩ᵢ, w, d, u)  # storage

    ϵ = epsᴸ(L)
    ball!(𝐩ᵢ, p[end], ϵ)
    setprecision(L) do
        copy!(Σ, 𝐩ᵢ)
        for i ∈ length(p)-1:-1:1
            # Σ = muladd(Σ, x, p[i])
            rmul!(Σ, 𝐱, w, d, u)
            ball!(𝐩ᵢ, p[i], ϵ)
            radd!(Σ, 𝐩ᵢ)
        end
        Σ
    end
end


text:
# This file provides an introductory tour of Onda.jl by generating, storing,
# and loading a toy Onda dataset. Run lines in the REPL to inspect output at
# each step! Tests are littered throughout to demonstrate functionality in a
# concrete manner, and so that we can ensure examples stay updated as the
# package evolves.
#
# NOTE: You should read https://github.com/beacon-biosignals/README.md#the-onda-format-specification
# before and/or alongside the completion of this tour; it explains the
# purpose/structure of the format.

using Onda, Legolas, Arrow, TimeSpans, DataFrames, Dates, UUIDs, Test

#####
##### generate some mock data
#####
# Let's kick off the tour by generating some mock data to play with in subsequent sections!
#
# Onda is primarily concerned with manipulating 3 interrelated entities. Paraphrasing from the
# Onda specification, these entities are:
#
# - "signals": A signal is the digitized output of a process, comprised of metadata (e.g. LPCM encoding,
#   channel information, sample data path/format information, etc.) and associated multi-channel sample
#   data.
#
# - "recordings": A recording is a collection of one or more signals recorded simultaneously over some
#   time period.
#
# - "annotations": An annotation is a a piece of (meta)data associated with a specific time span within
#   a specific recording.
#
# Signals and annotations are serialized as Arrow tables, while each sample data file is serialized to
# the file format specified by its corresponding signal's metadata. A "recording" is simply the collection
# of signals and annotations that share a common `recording` field.
#
# Below, we generate a bunch of signals/annotations across 10 recordings, writing the corresponding
# Arrow tables and sample data files to a temporary directory.

saws(info, duration) = [(j + i) % 100 * info.sample_resolution_in_unit for
                        i in 1:channel_count(info), j in 1:sample_count(info, duration)]

root = mktempdir()

signals = Signal[]
signals_recordings = [uuid4() for _ in 1:2]
for recording in signals_recordings
    for (kind, channels) in ("eeg" => ["fp1", "f3", "c3", "p3",
                                       "f7", "t3", "t5", "o1",
                                       "fz", "cz", "pz",
                                       "fp2", "f4", "c4", "p4",
                                       "f8", "t4", "t6", "o2"],
                             "ecg" => ["avl", "avr"],
                             "spo2" => ["spo2"])
        file_format = rand(("lpcm", "lpcm.zst"))
        file_path = joinpath(root, string(recording, "_", kind, ".", file_format))
        Onda.log("generating $file_path...")
        info = SamplesInfo(; kind=kind, channels=channels,
                           sample_unit="microvolt",
                           sample_resolution_in_unit=rand((0.25, 1)),
                           sample_offset_in_unit=rand((-1, 0, 1)),
                           sample_type=rand((Float32, Int16, Int32)),
                           sample_rate=rand((128, 256, 143.5)))
        data = saws(info, Minute(rand(1:10)))
        samples = Samples(data, info, false)
        start = Second(rand(0:30))
        signal = store(file_path, file_format, samples, recording, start)
        push!(signals, signal)
    end
end
path_to_signals = joinpath(root, "test.onda.signal.arrow")
Onda.write_signals(path_to_signals, signals)
Onda.log("wrote out $path_to_signals")

annotations = Annotation[]
sources = (uuid4(), uuid4(), uuid4())
annotations_recordings = vcat(signals_recordings[1:end-1], uuid4()) # overlapping but not equal to signals_recordings
for recording in annotations_recordings
    for i in 1:rand(3:10)
        start = Second(rand(0:60))
        annotation = Annotation(; recording=recording, id=uuid4(), span=TimeSpan(start, start + Second(rand(1:30))),
                                rating=rand(1:100), quality=rand(("good", "bad")), source=rand(sources))
        push!(annotations, annotation)
    end
end
path_to_annotations = joinpath(root, "test.onda.annotation.arrow")
Onda.write_annotations(path_to_annotations, annotations)
Onda.log("wrote out $path_to_annotations")

#####
##### basic Onda + DataFrames patterns
#####
# Since signals and annotations are represented tabularly, any package
# that supports the Tables.jl interface can be used to interact with
# them. Here, we show how you can use DataFrames.jl to perform a variety
# of common operations.
#
# If you're going to be working with Onda frequently, then it's probably
# worthwhile to become fluent in Julia's Tables.jl/DataFrames.jl
# ecosystem. This tour will give you a solid head start!
#
# Note that most of these operations are only shown here on a single table
# to avoid redundancy, but these examples are generally applicable to both
# signals and annotations tables.

# Read Onda Arrow files into `DataFrame`s:
signals = DataFrame(Legolas.read(path_to_signals))
annotations = DataFrame(Legolas.read(path_to_annotations))

# Get all signals from a given recording:
target = rand(signals.recording)
view(signals, findall(==(target), signals.recording), :)

# One of the consumer/producer-friendly properties of Onda is that signals
# and annotations are both represented in flat tables, enabling you to easily
# impose whatever indexing structure is most convenient for your use case.
#
# For example, if you wish to primarily access signals by recording and kind,
# you can easily create a structure with that index via `groupby`:
target = rand(signals.recording)
grouped = groupby(signals, Cols(:recording, :kind))
grouped[(target, "eeg")]

# Group/index signals + annotations by recording together:
target = rand(signals.recording)
dict = Legolas.gather(:recording, signals, annotations)
dict[target]

# Count number of signals in each recording:
combine(groupby(signals, :recording), nrow)

# Grab the longest signal in each recording:
combine(s -> s[argmax(duration.(s.span)), :], groupby(signals, :recording))

# Grab all multichannel signals greater than 5 minutes long:
filter(s -> length(s.channels) > 1 && duration(s.span) > Minute(5), signals)

# Load all sample data for a given recording:
target = rand(signals.recording)
transform(view(signals, findall(==(target), signals.recording), :),
          AsTable(:) => ByRow(load) => :samples)

# `mmap` sample data for a given LPCM signal:
target = first(s for s in eachrow(signals) if s.file_format == "lpcm")
Onda.mmap(target)

# Delete all sample data for a given recording (uncomment the
# inline-commented section to actual delete filtered signals'
# sample data!):
target = rand(signals.recording)
signals_copy = copy(signals) # we're gonna keep using `signals` afterwards, so let's work with a copy
filter!(s -> s.recording != target #=|| (rm(s.file_path); false)=#, signals_copy)

# Merge overlapping annotations of the same `quality` in the same recording.
# `merged` is an annotations table with a custom column of merged ids:
merged = DataFrame(mapreduce(merge_overlapping_annotations, vcat, groupby(annotations, :quality)))
m = rand(eachrow(merged)) # let's get the original annotation(s) from this merged annotation
view(annotations, findall(in(m.from), annotations.id), :)

# Load all the annotated segments that fall within a given signal's timespan:
within_signal(ann, sig) = ann.recording == sig.recording && TimeSpans.contains(sig.span, ann.span)
sig = first(sig for sig in eachrow(signals) if any(within_signal(ann, sig) for ann in eachrow(annotations)))
transform(filter(ann -> within_signal(ann, sig), annotations),
          :span => ByRow(span -> load(sig, translate(span, -start(sig.span)))) => :samples)

# In the above, we called `load(sig, span)` for each `span`. This invocation attempts to load
# *only* the sample data corresponding to `span`, which can be very efficient if the sample data
# file format + storage system supports random access and the full sample data file is very large.
# However, if random access isn't supported, or the sample data file is relatively small, or the
# requested set of `span`s heavily overlap, this approach may be less efficient than simply loading
# the whole file upfront. Here we demonstrate the latter as an alternative (note: in the future, we
# want to support an optimal batch loader):
samples = load(sig)
transform(filter(ann -> within_signal(ann, sig), annotations),
          :span => ByRow(span -> view(samples, :, translate(span, -start(sig.span)))) => :samples)

#####
##### working with `Samples`
#####
# A `Samples` struct wraps a matrix of interleaved LPCM-encoded (or decoded) sample data,
# along with a `SamplesInfo` instance that allows this matrix to be encoded/decoded.
# In this matrix, the rows correspond to channels and the columns correspond to timesteps.

# Let's grab a `Samples` instance for one of our mock EEG signals:
eeg_signal = signals[findfirst(==("eeg"), signals.kind), :]
eeg = load(eeg_signal)

# # Here are some basic functions for examining `Samples` instances:
@test eeg isa Samples && !eeg.encoded
@test sample_count(eeg) == sample_count(eeg_signal, duration(eeg)) == index_from_time(eeg.info.sample_rate, duration(eeg)) - 1
@test channel_count(eeg) == channel_count(eeg_signal) == length(eeg.info.channels)
@test channel(eeg, "f3") == channel(eeg_signal, "f3") == findfirst(==("f3"), eeg.info.channels)
@test channel(eeg, 2) == channel(eeg_signal, 2) == eeg.info.channels[2]
@test duration(eeg) == duration(eeg_signal.span)

# Here are some basic indexing examples using `getindex` and `view` wherein
# channel names and sample-rate-agnostic `TimeSpan`s are employed as indices:
span = TimeSpan(Second(3), Second(9))
span_range = index_from_time(eeg.info.sample_rate, span)
@test eeg[:, span].data == view(eeg, :, span_range).data
@test eeg["f3", :].data == view(eeg, channel(eeg, "f3"), :).data
@test eeg["f3", 1:10].data == view(eeg, channel(eeg, "f3"), 1:10).data
@test eeg["f3", span].data == view(eeg, channel(eeg, "f3"), span_range).data
rows = ["f3", "c3", "p3"]
@test eeg[rows, 1:10].data == view(eeg, channel.(Ref(eeg), rows), 1:10).data
rows = ["c3", 4, "f3"]
@test eeg[rows, span].data == view(eeg, channel.(Ref(eeg), rows), span_range).data

# One can also index rows by regular expressions. For example, to match all the
# channels which have an `f`:
f_channels = ["fp1", "f3","f7", "fz", "fp2", "f4", "f8"]
@test eeg[r"f", span].data == view(eeg, channel.(Ref(eeg), f_channels), span_range).data

# Onda overloads the necessary Arrow.jl machinery to enable individual sample data
# segments (specifically, `Samples` and `SamplesInfo` values) to be (de)serialized
# to/from Arrow for storage or IPC purposes; see below for an example. Note that if
# you wanted to use Arrow as a storage format for whole sample data files w/ Onda,
# it'd make more sense to create an `AbstractLPCMFormat` subtype for your Arrow <-> LPCM
# mapping (an example of this can be seen in `examples/flac.jl`).
x = (a=[eeg], b=[eeg.info])
y = Arrow.Table(Arrow.tobuffer(x))
@test x.a == y.a
@test x.b == y.b

# Note that `Samples` is not an `AbstractArray` subtype; the special indexing
# behavior above is only defined for convenient data manipulation. It is fine
# to access the sample data matrix directly via the `data` field if you need
# to manipulate the matrix directly or pass it to downstream computations.


text:
# test InertialPose3

using RoME, Distributions
using JLD

@load joinpath(dirname(@__FILE__),"preintstationarydata.jld")
# DATA[1]

N = 100
fg = initfg()

initCov = 0.001*eye(15)
[initCov[i,i] = 0.0001^2 for i in 4:6];
[initCov[i,i] = 0.0002^2 for i in 10:15];
odoCov = deepcopy(initCov)

println("Adding PriorInertialPose3 to graph...")
v1 = addNode!(fg, :x1,  0.1*randn(15,N),  N=N)
initPosePrior = PriorInertialPose3( MvNormal( zeros(15), initCov) )
f1  = addFactor!(fg, [v1], initPosePrior)

n = 1
v2 = addNode!(fg, :x2, dims=15,  N=N)
noise = MvNormal(zeros(15),(DATA[n][3]+DATA[n][3]')*0.5 )
inerodo = InertialPose3(noise,DATA[n][1],DATA[n][2])
f2  = addFactor!(fg, [v1;v2], inerodo )

initializeNode!(fg, :x2, N=N)

# drawDensityMesh(fg, :x2)

plotKDE(fg, :x2, dims=[1;2]  , title="x,y")
plotKDE(fg, :x2, dims=[2;3]  , title="y,z")
plotKDE(fg, :x2, dims=[6;3]  , title="ψ,z")
plotKDE(fg, :x2, dims=[4;5]  , title=",")
plotKDE(fg, :x2, dims=[7;8]  , title=",")
plotKDE(fg, :x2, dims=[9;10] , title=",")
plotKDE(fg, :x2, dims=[11;12], title=",")
plotKDE(fg, :x2, dims=[13;14], title=",")
plotKDE(fg, :x2, dims=[15]   , title=",")

DATA[n][3][1,1]

res = zeros(15)
idx = 1
# pos, so3, vel, bw, ba
meas = vectoarr2([0;0;9.81/2; 0;0;0; 0;0;9.81; 0;0;0; 0;0;0])
wIPi = zeros(15,1)
wIPj = zeros(15,1)

inerodo(res, idx, (meas,), wIPi, wIPj)
# function (ip3::InertialPose3)(
#         res::Vector{Float64},
#         idx::Int,
#         meas::Tuple,
#         wIPi::Array{Float64,2},
#         wIPj::Array{Float64,2}  )

using Optim

ggo = (x) -> inerodo(res, idx, (meas,), wIPi, vectoarr2(x))

ggos = (x) -> ggo([x[1:5]...,0.0,x[6:end]...])

ret = optimize(ggo, zeros(15))

using Gadfly

obj = (x) -> ggo([0,0.0,0, 0,0,x, 0,0,0, 0,0,0, 0,0,0])
plot(obj, -2, 2)

ran = linspace(-2,2,100)
ggoxy = (x,y) -> ggo([0,0,0.0, 0,0,x, y,0,0, 0,0,0, 0,0,0])
plot(z=ggoxy, x=ran, y=ran, Geom.contour)

@show res
using NLsolve

gg = (res, x) -> inerodo(res, idx, (meas,), wIPi, vectoarr2(x))

ret = nlsolve(gg, wIPj[:])

@show ret.zero

using KernelDensityEstimate
p = kde!( getSample(inerodo, 100)[1] )
plotKDE(p, dims=[1;2]  )
plotKDE(p, dims=[6;3]  )
plotKDE(p, dims=[4;5]  )
plotKDE(p, dims=[7;8]  )
plotKDE(p, dims=[9;10] )
plotKDE(p, dims=[11;12])
plotKDE(p, dims=[13;14])
plotKDE(p, dims=[15]   )

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl", 22)

using Gadfly

spy(DATA[n][3])

norm(DATA[n][3]-DATA[n][3]')

tree = wipeBuildNewTree!(fg)

inferOverTree!(fg, tree)

plotKDE(fg, :x1, dims=[4;5])

writeGraphPdf(fg)
run(`evince fg.pdf`)

#


text:
function configs_to_dataframe(configs::Dict; path::String, write::Bool=false)
    config_df = init_configs_dataframe()
    progress_bar = ProgressMeter.Progress(length(keys(configs)), 0.2, "Converting configurations...")
    for config_key in keys(configs)
        config_row = extract_config(configs, config_key)
        push!(config_df, config_row)
        ProgressMeter.next!(progress_bar)
    end
    if write
        write_dataframe(df=config_df, filename="configs.feather", path=path)
    end
    return config_df
end

function init_configs_dataframe()
    return DataFrames.DataFrame(
        Config=String[],
        Nodes=Int64[],
        N0=Int64[],
        K=Int64[],
        AgentCount=Int64[],
        InitialTrash=Int64[],
        InitialLitterers=Int64[],
        GarbageCollection=Int64[],
        Iterations=Int64[]
    )
end

function extract_config(configs::Dict, config_key::Any)
    return (
        Config="config_" * lpad(config_key, 2, "0"),
        Nodes=configs[config_key].n_nodes,
        N0=configs[config_key].n0,
        K=configs[config_key].k,
        AgentCount=configs[config_key].agentcount,
        InitialTrash=configs[config_key].initial_trash,
        InitialLitterers=configs[config_key].initial_litterers,
        GarbageCollection=configs[config_key].garbage_collection,
        Iterations=configs[config_key].n_iter
    )
end

function names_to_camelcase!(df::DataFrames.DataFrame)
    new_colnames = [to_camelcase(name) for name in names(df)]
    rename!(df, Symbol.(new_colnames))
    return df
end

function to_camelcase(varname::String)
	varname_split_titlecase = titlecase.(split(deepcopy(varname), "_"))
	varname_modified = join(varname_split_titlecase, "")
	return varname_modified
end


text:
using Revise
using LinearAlgebra
using Plots

includet("src/InternalArray.jl")
includet("src/FixedStepMethods.jl")
includet("src/MixedStepMethods.jl")
includet("src/Interpolation.jl")
includet("src/Arenstorf.jl")
IArray = InternalArray

const closing_err = 5e-2
const stop_time = 40
const mstep = 1e-5

x, time = @timed fixed_step_ode_solver(
    Arenstorf.derivative,
    0,
    [Arenstorf.p0; Arenstorf.v0],
    (t, x) ->
        t > 20 && norm(Arenstorf.p0 .- last(x)[1:2]) < closing_err || t > stop_time,
    step = mstep,
    maxit = 10^7,
    method = euler_step,
)

p = Arenstorf.init_plot()
inds = [i for i in 1:length(x)÷200:length(x)]
if last(inds) != length(x)
    push!(inds, length(x))
end
plot!(
    IArray(x[inds], 1),
    IArray(x[inds], 2),
    label = "satellite",
    title = "Arenstorf Orbit - Euler's Method",
)
savefig("media/arenstorf/euler.png")

begin
    print("\033[32;1marenstorf>\033[0m \n")
    print("\tMethod: Euler's\n")
    print("\tStep: $(mstep)\n")
    print("\tTarget Error (closing condition): $(closing_err)\n")
    println()
    numit = length(x)-1
    print("\tNumber of iterations: $(numit)\n")
    print("\tRunning time: $(time)s\n")
    print("\tTime (solution): $(numit*mstep)\n")
    print("\tNumber of evaluations: $(numit)\n")
    print("\tOrbit closed with target tolerance: $(numit*mstep < stop_time)\n")
    print("\tDistance from last point to start point: ",
          "$(norm(last(x)[1:2] .- Arenstorf.p0))\n")
    n = length(x)
    last5x = [p[1] for p in x[n-4:n]]
    last5y = [p[2] for p in x[n-4:n]]
    closingx = Interpolation.hermite(last5y, last5x, 0)
    print("\tDistance from interpolated closing point to start point: ",
          "$(abs(closingx - Arenstorf.p0[1]))\n")
end

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: March 2022
# Copyright: Georgia Institute of Technology, 2022

export _wsave, sef_plot_configs

_wsave(s, fig::Figure; dpi::Int = 200) = fig.savefig(s, bbox_inches = "tight", dpi = dpi)

function sef_plot_configs(; fontsize = 12)
    set_style("whitegrid")
    rc("font", family = "serif", size = fontsize)
    font_prop = matplotlib.font_manager.FontProperties(
        family = "serif",
        style = "normal",
        size = fontsize,
    )
    sfmt = matplotlib.ticker.ScalarFormatter(useMathText = true)
    sfmt.set_powerlimits((0, 0))
    matplotlib.use("Agg")

    return font_prop, sfmt
end


text:
"""
$(SIGNATURES)

Implements a sparse regression, given an `AbstractOptimizer` or `AbstractSubspaceOptimizer`.
`X` denotes the coefficient matrix, `A` the design matrix and `Y` the matrix of observed or target values.
`X` can be derived via `init(opt, A, Y)`.
`maxiter` indicates the maximum iterations for each call of the optimizer, `abstol` the absolute tolerance of
the difference between iterations in the 2 norm. If the optimizer is called with a `Vector` of thresholds, each `maxiter` indicates
the maximum iterations for each threshold.

If `progress` is set to `true`, a progressbar will be available. `progress_outer` and `progress_offset` are used to compute the initial offset of the
progressbar.

If used with a `Vector` of thresholds, the functions `f` with signature `f(X, A, Y)` and `g` with signature `g(x, threshold) = G(f(X, A, Y))` with the arguments given as stated above can be passed in. These are
used for finding the pareto-optimal solution to the sparse regression.
"""
function sparse_regression!(X, A, Y, opt::AbstractOptimizer{T};
    maxiter::Int = maximum(size(A)),
    abstol = eps(eltype(T)), progress::Bool = false,
    f::Function = F(opt),
    g::Function = G(opt),
    progress_outer::Int = 1, progress_offset::Int = 0, kwargs...) where T

    # Closure for the pareto function
    fg(x, A, y, lambda) = (g∘f)(x, A, y, lambda)

    # TODO Tmp Result
    X_tmp = similar(X)
    X_tmp .= X

    λ = get_threshold(opt)
    λ = issorted(λ) ? λ : sort(λ)
    λs = zeros(eltype(λ), size(Y,2))

    if progress
        progress =  init_progress(opt, X, A, Y, length(λ)*progress_outer, progress_offset)
    else
        progress = nothing
    end

    obj = zero(eltype(X))
    objtmp = zero(eltype(X))
    sparsity = zero(eltype(X))
    sparsitytmp = zero(eltype(X))

    @views for (i,λi) in enumerate(λ)
        init!(X_tmp, opt, A, Y)

        opt(X_tmp, A, Y, λi, maxiter = maxiter, abstol = abstol, f = f, g = g)

        # Increasing the threshold makes no sense
        all(iszero(X_tmp)) && break

        for j in 1:size(Y, 2)
            if fg(X_tmp[:, j], A, Y[:, j], λi) < fg(X[:, j], A, Y[:, j], λi)
                λs[j] =λi
                X[:, j] .= X_tmp[:, j]
            end
        end

        if !isnothing(progress)
            sparsity, obj = f(X, A, Y, λi)
            sparsitytmp, objtmp = f(X_tmp, A, Y, λi)

            ProgressMeter.next!(
            progress;
            showvalues = [
                (:Threshold, λi), (Symbol("Best Objective"), obj), (Symbol("Best Sparsity"), sparsity),
                (Symbol("Current Objective"), objtmp), (Symbol("Current Sparsity"), sparsitytmp)
            ]
            )
        end
    end

    for i in 1:size(Y, 2)
        @views clip_by_threshold!(X[:, i], λs[i])
    end

    return λs
end


text:
############################################################################################
#					             CONSTRAINT LIST										   #
############################################################################################
"""
	AbstractConstraintSet

Stores constraint error and Jacobian values, correctly accounting for the error state if
necessary.

# Interface
- `get_convals(::AbstractConstraintSet)::Vector{<:ConVal}` where the size of the Jacobians
	match the full state dimension
- `get_errvals(::AbstractConstraintSet)::Vector{<:ConVal}` where the size of the Jacobians
	match the error state dimension
- must have field `c_max::Vector{<:AbstractFloat}` of length `length(get_convals(conSet))`

# Methods
Once the previous interface is defined, the following methods are defined
- `Base.iterate`: iterates over `get_convals(conSet)`
- `Base.length`: number of independent constraints
- `evaluate!(conSet, Z::Traj)`: evaluate the constraints over the entire trajectory `Z`
- `jacobian!(conSet, Z::Traj)`: evaluate the constraint Jacobians over the entire trajectory `Z`
- `error_expansion!(conSet, model, G)`: evaluate the Jacobians for the error state using the
	state error Jacobian `G`
- `max_violation(conSet)`: return the maximum constraint violation
- `findmax_violation(conSet)`: return details about the location of the maximum
	constraint violation in the trajectory
"""
abstract type AbstractConstraintSet end

"""
	ConstraintList

Stores the set of constraints included in a trajectory optimization problem. Includes a list
of both the constraint types [`AbstractConstraint`](@ref) as well as the knot points at which
the constraint is applied. Each constraint is assumed to apply to a contiguous set of knot points.

A `ConstraintList` supports iteration and indexing over the `AbstractConstraint`s, and
iteration of both the constraints and the indices of the knot points at which they apply
via `zip(cons::ConstraintList)`.

Constraints are added via the [`add_constraint!`](@ref) method, which verifies that the constraint
dimension is consistent with the state and control dimensions of the problem.

The total number of constraints at each knot point can be queried using the
[`num_constraints`](@ref) method.

The constraint list can also be sorted to separate `StageConstraint`s and `CoupledConstraint`s
via the `sort!` method.

A constraint list can be queried if it has a `DynamicsConstraint` via
`has_dynamics_constraint(::ConstraintList)`.

# Constructor
	ConstraintList(n::Int, m::Int, N::Int)
"""
struct ConstraintList <: AbstractConstraintSet
	n::Int
	m::Int
	constraints::Vector{AbstractConstraint}
	inds::Vector{UnitRange{Int}}
	p::Vector{Int}
	function ConstraintList(n::Int, m::Int, N::Int)
		constraints = AbstractConstraint[]
		inds = UnitRange{Int}[]
		p = zeros(Int,N)
		new(n, m, constraints, inds, p)
	end
end

"""
	add_constraint!(cons::ConstraintList, con::AbstractConstraint, inds::UnitRange, [idx])

Add constraint `cons` to `ConstraintList` `cons` for knot points given by `inds`.

Use `idx` to determine the location of the constraint in the constraint list.
`idx=-1` (default) adds the constraint at the end of the list.

# Example
Here is an example of adding a goal and control limit constraint for a cartpole swing-up.
```julia
# Dimensions of our problem
n,m,N = 4,1,51    # 51 knot points

# Create our list of constraints
cons = ConstraintList(n,m,N)

# Create the goal constraint
xf = [0,π,0,0]
goalcon = GoalConstraint(xf)
add_constraint!(cons, goalcon, N)  # add to the last time step

# Create control limits
ubnd = 3
bnd = BoundConstraint(n,m, u_min=-ubnd, u_max=ubnd, idx=1)  # make it the first constraint
add_constraint!(cons, bnd, 1:N-1)  # add to all but the last time step

# Indexing
cons[1] === bnd                            # (true)
cons[2] === goal                           # (true)
allcons = [con for con in cons]
cons_and_inds = [(con,ind) in zip(cons)]
cons_and_inds[1] == (bnd,1:n-1)            # (true)
```
"""
function add_constraint!(cons::ConstraintList, con::AbstractConstraint, inds::UnitRange{Int}, idx=-1)
	@assert check_dims(con, cons.n, cons.m) "New constaint not consistent with n=$(cons.n) and m=$(cons.m)"
	@assert inds[end] <= length(cons.p) "Invalid inds, inds[end] must be less than number of knotpoints, $(length(cons.p))"
	if isempty(cons)
		idx = -1
	end
	if idx == -1
		push!(cons.constraints, con)
		push!(cons.inds, inds)
	elseif 0 < idx <= length(cons)
		insert!(cons.constraints, idx, con)
		insert!(cons.inds, idx, inds)
	else
		throw(ArgumentError("cannot insert constraint at index=$idx. Length = $(length(cons))"))
	end
	num_constraints!(cons)
	@assert length(cons.constraints) == length(cons.inds)
end

@inline add_constraint!(cons::ConstraintList, con::AbstractConstraint, k::Int, idx=-1) =
	add_constraint!(cons, con, k:k, idx)

# Iteration
Base.iterate(cons::ConstraintList) = length(cons) == 0 ? nothing : (cons[1], 1)
Base.iterate(cons::ConstraintList, i::Int) = i < length(cons) ? (cons[i+1], i+1) : nothing
@inline Base.length(cons::ConstraintList) = length(cons.constraints)
Base.IteratorSize(::ConstraintList) = Base.HasLength()
Base.IteratorEltype(::ConstraintList) = Base.HasEltype()
Base.eltype(::ConstraintList) = AbstractConstraint
Base.firstindex(::ConstraintList) = 1
Base.lastindex(cons::ConstraintList) = length(cons.constraints)

Base.zip(cons::ConstraintList) = zip(cons.inds, cons.constraints)

@inline Base.getindex(cons::ConstraintList, i::Int) = cons.constraints[i]

for method in (:deepcopy, :copy)
	@eval function Base.$method(cons::ConstraintList)
		cons2 = ConstraintList(cons.n, cons.m, length(cons.p))
		for i in eachindex(cons.constraints)
			con_ = $(method == :deepcopy ? :(copy(cons.constraints[i])) : :(cons.constraints[i]))
			add_constraint!(cons2, con_, copy(cons.inds[i]))
		end
		return cons2
	end
end

"""
	num_constraints(::ConstraintList)
	num_constraints(::Problem)
	num_constraints(::TrajOptNLP)

Return a vector of length `N` constaining the total number of constraint values at each
knot point.
"""
@inline num_constraints(cons::ConstraintList) = cons.p

function num_constraints!(cons::ConstraintList)
	cons.p .*= 0
	for i = 1:length(cons)
		p = length(cons[i])
		for k in cons.inds[i]
			cons.p[k] += p
		end
	end
end

function change_dimension(cons::ConstraintList, n::Int, m::Int, ix=1:n, iu=1:m)
	new_list = ConstraintList(n, m, length(cons.p))
	for (i,con) in enumerate(cons)
		new_con = change_dimension(con, n, m, ix, iu)
		add_constraint!(new_list, new_con, cons.inds[i])
	end
	return new_list
end

# sort the constraint list by stage < coupled, preserving ordering
function Base.sort!(cons::ConstraintList; rev::Bool=false)
	lt(con1,con2) = false
	lt(con1::StageConstraint, con2::CoupledConstraint) = true
	inds = sortperm(cons.constraints, alg=MergeSort, lt=lt, rev=rev)
	permute!(cons.inds, inds)
	permute!(cons.constraints, inds)
	return cons
end

function has_dynamics_constraint(conSet::ConstraintList)
	for con in conSet
		if con isa DynamicsConstraint
			return true
		end
	end
	return false
end

############################################################################################
#                                 NLP FUNCTIONS
############################################################################################
"""
	primal_bounds!(zL, zU, cons::ConstraintList; remove=true)

Get the lower and upper bounds on the primal variables imposed by the constraints in `cons`,
where `zL` and `zU` are vectors of length `NN`, where `NN` is the total number of primal
variables in the problem. Returns the modified lower bound `zL` and upper bound `zU`.

If any of the bound constraints are redundant, the strictest bound is returned.

If `remove = true`, these constraints will be removed from `cons`.
"""
function primal_bounds!(zL, zU, cons::ConstraintList, remove::Bool=true)
	NN = length(zL)
	n,m = cons.n, cons.m
	isequal = NN % (n+m) == 0
	N = isequal ? Int(NN / (n+m)) : Int((NN+m)/(n+m))
	rm_inds = Int[]
	for (j,(inds,con)) in enumerate(zip(cons))
		is_bound = false
		for (i,k) in enumerate(inds)
			off = (k-1)*(n+m)
			if !isequal && k == N  # don't allow the indexing to go out of bounds at the last time step
				zind = off .+ (1:n)
			else
				zind = off .+ (1:n+m)
			end
			is_bound |= primal_bounds!(view(zL, zind), view(zU, zind), con)
		end
		is_bound && push!(rm_inds,j)
	end
	if remove
		deleteat!(cons.constraints, rm_inds)
		deleteat!(cons.inds, rm_inds)
	end
	return zL, zU
end

"""
	JacobianStructure

Contains information describing the structure (sparsity) of the constraint Jacobian.
"""
struct JacobianStructure
	NN::Int  # primal variables
	P::Int   # dual variables (constraints)
	nD::Int  # non-zero entries in constraint Jacobian
	cinds::Vector{Vector{UnitRange{Int}}}  # indices into the constraint vector
	zinds::Vector{Matrix{UnitRange{Int}}}  # indices into primal vector
	linds::Vector{Matrix{UnitRange{Int}}}  # indices into Jacobian vector
end

function JacobianStructure(cons::ConstraintList, structure=:by_knotpoint)
	n,m = cons.n, cons.m
    N = length(cons.p)
	isequal = false
	if has_dynamics_constraint(cons)
		isequal = integration(cons[end]) <: Implicit
	end
	NN = num_vars(cons.n, cons.m, N, isequal)
	P = sum(num_constraints(cons))
    numcon = length(cons.constraints)

	zinds_ = gen_zinds(n, m, N, isequal)
	cinds = [[1:0 for j in eachindex(cons.inds[i])] for i in 1:length(cons)]
	zinds = [[1:0 for j in eachindex(cons.inds[i]), k in widths(cons[i],n,m)] for i in 1:length(cons)]
	linds = [[1:0 for j in eachindex(cons.inds[i]), k in widths(cons[i],n,m)] for i in 1:length(cons)]

    # Dynamics and general constraints
    idx = 0
	nD = 0
	if structure == :by_constraint
	    for (i,con) in enumerate(cons.constraints)
			len = length(con)
			zs = get_inds(con)
			for (j,k) in enumerate(cons.inds[i])
				inds = idx .+ (1:len)
				cinds[i][j] = inds
				idx += len
				off_z = (k-1)*(n+m)
				for l in length(zs)
					zi = zs[l] .+ off_z       # indices of depedent primal variables
					blk_len = len*length(zi)  # number of elements in the Jacobian block
					zinds[i][j,l] = zi
					linds[i][j,l] = nD .+ (1:blk_len)
					nD += blk_len
				end
	        end
	    end
	elseif structure == :by_knotpoint
		for k = 1:N
			for (i,con) in enumerate(cons)
				inds = cons.inds[i]
				len = length(con)
				zs = get_inds(con,n,m)
				if k in inds
					j = k - inds[1] + 1
					cinds[i][j] = idx .+ (1:len)
					idx += len

					off_z = (k-1)*(n+m)
					for l in 1:length(zs)
						zi = zs[l] .+ off_z       # indices of depedent primal variables
						blk_len = len*length(zi)  # number of elements in the Jacobian block
						zinds[i][j,l] = zi
						linds[i][j,l] = nD .+ (1:blk_len)
						nD += blk_len
					end
				end
			end
		end
	end
	return JacobianStructure(NN, P, nD, cinds, zinds, linds)
end

function jacobian_structure!(D::AbstractMatrix, jac::JacobianStructure)
	ncons = length(jac.cinds)
	for i = 1:ncons
		cinds = jac.cinds[i]
		zinds = jac.zinds[i]
		linds = jac.linds[i]
		for j = 1:length(cinds)
			for k = 1:size(linds,2)
				D[cinds[j], zinds[j,k]] = linds[j,k]
			end
		end
	end
	return D
end

@inline jacobian_structure(jac::JacobianStructure) =
	jacobian_structure!(spzeros(Int,jac.P, jac.NN), jac)

function gen_convals(Dv::AbstractVector, d::AbstractVector, cons::ConstraintList, jac=JacobianStructure(cons))
	n,m = cons.n, cons.m  # TODO: account for state diff size
	C = map(enumerate(zip(cons))) do (i,(inds,con))
		len = length(con)
		ws = widths(con, n, m)
		[reshape(view(Dv, jac.linds[i][j,l]), len, ws[l])
			for j = 1:length(inds), l = 1:length(ws)]
	end
	c = map(enumerate(zip(cons))) do (i,(inds,con))
		[view(d, jac.cinds[i][j]) for j = 1:length(inds)]
	end
	return C,c
end

function gen_convals(D::AbstractMatrix, d::AbstractVector, cons::ConstraintList, jac=JacobianStructure(cons))
	n,m = cons.n, cons.m  # TODO: account for state diff size
	ncons = length(cons)
	C = map(enumerate(zip(cons))) do (i,(inds,con))
		len = length(con)
		ws = widths(con, n, m)
		[view(D, jac.cinds[i][j], jac.zinds[i][j,l])
			for j = 1:length(inds), l = 1:length(ws)]
	end

	c = map(enumerate(zip(cons))) do (i,(inds,con))
		[view(d, jac.cinds[i][j]) for j = 1:length(inds)]
	end
    return C,c
end

"""
	gen_con_inds(cons::ConstraintList, structure::Symbol)

Generate the indices into the concatenated constraint vector for each constraint.
Determines the bandedness of the Jacobian
"""
function gen_con_inds(conSet::ConstraintList, structure=:by_knotpoint)
    return cinds, linds
end

# function jacobian_structure(cons::ConstraintList)
# 	isequal = false
# 	if has_dynamics_constraint(cons)
# 		isequal = integration(cons[end]) <: Implicit
# 	end
# 	N = length(cons.p)
# 	NN = num_vars(cons.n, cons.m, N, isequal)
# 	P = sum(num_constraints(cons))
# 	D = spzeros(NN,P)
#
# 	idx = 0
# 	for k = 1:N
# 		for (i,conval) in enumerate(get_convals(conSet))
# 			p = length(conval.con)
# 			if k in conval.inds
# 				for (l,w) in enumerate(widths(conval.con))
# 					blk_len = p*w
# 					inds = reshape(idx .+ (1:blk_len), p, w)
# 					j = _index(conval,k)
# 					# linds[i][j] = inds
# 					conval.jac[j,l] .= inds
# 					idx += blk_len
# 				end
# 			end
# 		end
# 	end
# end


text:
# ------------ Execution pipelines ------------ #

# This abstract type represents a "compiler execution context".
# This context is allowed to look at certain calls, and re-interpret them using dispatch.
# This allows the expression of tracing as a tracked side-effect.

abstract type ExecutionContext{J, K, T} end

# These are "soft" interfaces, not all of these methods apply to every subtype of ExecutionContext.

increment!(ctx::T, w::F) where {T <: ExecutionContext, F <: AbstractFloat} = ctx.weight += w

visit!(ctx::T, addr) where T <: ExecutionContext = visit!(ctx.visited, addr)

get_prev(ctx::T, addr) where T <: ExecutionContext = get_sub(ctx.prev, addr)

function add_choice!(ctx::T, addr, score, v) where T <: ExecutionContext
    ctx.score += score
    set_sub!(ctx.tr, addr, Choice(score, v))
end

function add_value!(ctx::T, addr, score, v) where T <: ExecutionContext
    ctx.score += score
    set_sub!(ctx.map, addr, Value(v))
end

function add_call!(ctx::T, addr, cs::C) where {T <: ExecutionContext, C <: CallSite}
    ctx.score += get_score(cs)
    set_sub!(ctx.tr, addr, cs)
end

# ----------- Control compiler passes with options ------------ #

abstract type CompilationOptions end

struct DefaultCompilationOptions{H, S, Sp, AA} <: CompilationOptions end

const NoStatic = DefaultCompilationOptions{false, false, false, false}
const DefaultPipeline = DefaultCompilationOptions{false, true, false, false}
const StaticWithLints = DefaultCompilationOptions{true, true, false, false}
const SpecializerPipeline = DefaultCompilationOptions{true, true, true, false}
const AutomaticAddressingPipeline = DefaultCompilationOptions{true, true, false, true}

extract_options(::Type{DefaultCompilationOptions{H, S, Sp, AA}}) where {H, S, Sp, AA} = (H = H, S = S, Sp = Sp, AA = AA)

# ------------ includes ------------ #

# Note - includes must come before the dynamo definitions because of world age issues.

include("pipelines/gradient_store.jl")
include("pipelines/contexts.jl")

const DiffAware = Union{UpdateContext, RegenerateContext}
const DoesNotCare = Union{SimulateContext, GenerateContext, AssessContext, ProposeContext, ForwardModeContext, BackpropagationContext}
const AllContexts = Union{SimulateContext, GenerateContext, AssessContext, ProposeContext, ForwardModeContext, BackpropagationContext}

include("pipelines/default_pipeline.jl")
include("pipelines/gfi/generate.jl")
include("pipelines/gfi/simulate.jl")
include("pipelines/gfi/propose.jl")
include("pipelines/gfi/assess.jl")
include("pipelines/gfi/update.jl")
include("pipelines/gfi/regenerate.jl")
include("pipelines/gfi/backpropagate.jl")
include("pipelines/gfi/forwardmode.jl")

# Staging for "doesn't care about incremental diff" contexts.
# These basically deal with straightforward execution - they don't require us to think about incremental computing.
@dynamo function (sx::ExecutionContext{J, K, T})(a...) where {J, K, T <: DoesNotCare}
    ir = IR(a...)
    ir == nothing && return
    ir = staged_pipeline(ir, J, K, T)
    ir
end

# Diff aware contexts are a bit tricky.
# Sneaky invoke hackz. Allows you to pass in arguments with types which don't match the dispatch table. Dangerous - but required for incremental computation on incremental computation unaware code.
@dynamo function (sx::ExecutionContext{J, K, T})(f, ::Type{S}, args...) where {J, S <: Tuple, K, T <: DiffAware}
    ir = IR(f, S.parameters...)
    ir == nothing && return
    ir = staged_pipeline(ir, J, K, UpdateContext)
    ir
end

# Each context has an associated function.
get_interface(::Type{SimulateContext}) = simulate
get_interface(::Type{GenerateContext}) = generate
get_interface(::Type{ProposeContext}) = propose
get_interface(::Type{AssessContext}) = assess
get_interface(::Type{RegenerateContext}) = regenerate
get_interface(::Type{UpdateContext}) = update

# Fixes for Base. Unfortunately we have to copy definitions for two different branches in the type hierarchy (required by the way we've split dynamos).
function (sx::ExecutionContext{J, K, T})(::typeof(Core._apply_iterate), f, c::typeof(trace), args...) where {J, K, T <: DiffAware}
    flt = flatten(args)
    addr, rest = flt[1], flt[2 : end]
    ret, cl = get_interface(T)(rest...)
    add_call!(sx, addr, cl)
    ret
end

function (sx::ExecutionContext{J, K, T})(::typeof(Core._apply_iterate), f, c::typeof(trace), args...) where {J, K, T <: DoesNotCare}
    flt = flatten(args)
    addr, rest = flt[1], flt[2 : end]
    ret, cl = get_interface(T)(rest...)
    add_call!(sx, addr, cl)
    ret
end

function (sx::ExecutionContext{J, K, T})(::typeof(Base.collect), generator::Base.Generator) where {J, K, T <: DiffAware}
    map(generator.iter) do i
        δ = Δ(i, NoChange()) # DiffAware contexts require this Diff arg.
        sx(generator.f, tupletype(δ), δ)
    end
end

function (sx::ExecutionContext{J, K, T})(::typeof(Base.collect), generator::Base.Generator) where {J, K, T <: DoesNotCare}
    map(generator.iter) do i
        sx(generator.f, i)
    end
end

function (sx::ExecutionContext{J, K, T})(::typeof(Base.map), fn, iter) where {J, K, T <: DiffAware}
    map(iter) do i
        δ = Δ(i, NoChange()) # DiffAware contexts require this Diff arg.
        sx(fn, tupletype(δ), δ)
    end
end

function (sx::ExecutionContext{J, K, T})(::typeof(Base.map), fn, iter) where {J, K, T <: DoesNotCare}
    map(iter) do i
        sx(fn, i)
    end
end

function (sx::ExecutionContext{J, K, T})(::typeof(Base.filter), fn, iter) where {J, K, T <: DiffAware}
    filter(iter) do i
        δ = Δ(i, NoChange()) # DiffAware contexts require this Diff arg.
        sx(fn, tupletype(δ), δ)
    end
end

function (sx::ExecutionContext{J, K, T})(::typeof(Base.filter), fn, iter) where {J, K, T <: DoesNotCare}
    filter(iter) do i
        sx(fn, i)
    end
end

# ------------ Documentation ------------ #


text:
"""
    latexarray{T}(arr::AbstractArray{T, 2})
Create a LaTeX array environment using [`latexraw`](@ref).

# Examples
```julia
arr = [1 2; 3 4]
latexarray(arr)
```
```math
"\\begin{equation}\n\\left[\n\\begin{array}{cc}\n1 & 2\\\\ \n3 & 4\\\\ \n\\end{array}\n\\right]\n\\end{equation}\n"
```
"""
function latexarray(arr::AbstractArray; adjustment::Symbol=:c, transpose=false, double_linebreak=false,
    starred=false, kwargs...)
    transpose && (arr = permutedims(arr))
    rows = first(size(arr))
    columns = length(size(arr)) > 1 ? size(arr)[2] : 1

    eol = double_linebreak ? " \\\\\\\\\n" : " \\\\\n"

    str = "\\begin{equation$(starred ? "*" : "")}\n"
    str *= "\\left[\n"
    str *= "\\begin{array}{" * "$(adjustment)"^columns * "}\n"

    arr = latexraw(arr; kwargs...)
    for i=1:rows, j=1:columns
        str *= arr[i,j]
        j==columns ? (str *= eol) : (str *= " & ")
    end

    str *= "\\end{array}\n"
    str *= "\\right]\n"
    str *= "\\end{equation$(starred ? "*" : "")}\n"
    latexstr = LaTeXString(str)
    COPY_TO_CLIPBOARD && clipboard(latexstr)
    return latexstr
end

latexarray(args::AbstractArray...; kwargs...) = latexarray(hcat(args...); kwargs...)
latexarray(arg::AbstractDict; kwargs...) = latexarray(collect(keys(arg)), collect(values(arg)); kwargs...)
latexarray(arg::Tuple...; kwargs...) = latexarray([collect(i) for i in arg]...; kwargs...)

function latexarray(arg::Tuple; kwargs...)
    if first(arg) isa Tuple || first(arg) isa AbstractArray
        return latexarray([collect(i) for i in arg]...; kwargs...)
    end
    return latexarray(collect(arg); kwargs...)
end


text:
## The Equation class
@fenicsclass Equation

#compare(form1::Form, form2::Form) = Equation(form1.pyobject[:__eq__](form2.pyobject) )
#Equation(form1::Form, form2::Form) = Equation(dolfin.Equation(form1.pyobject, form2.pyobject, 23) )
export Equation

@fenicsclass Vector
Vector(comm::MPI_Comm, n::Int) = Vector(dolfin.Vector(comm.pyobject, n))
Vector(n::Int) = Vector(mpi_comm_world(), n)
array(v::Vector) = fenicscall(v, :array)::Array{Float64}
export array

## The Function class
# (still problems with the import/export clash with Base.Function or Core.Function)
#@fenicsclass CoefficientFunction
@fenicsclass Function
Function(V::FunctionSpace) = Function(dolfin.Function(V.pyobject))

vector(f::Function) = Vector(fenicscall(f, :vector))
array(f::Function) = array(vector(f))
interpolate(ex::Expression, V::FunctionSpace) = Function(dolfin.interpolate(ex.pyobject, V.pyobject))

export Function, vector, interpolate

## LinearVariationalProblem and Solver
@fenicsclass LinearVariationalProblem
LinearVariationalProblem(a::Form, L::Form, u::Function) =
  LinearVariationalProblem(dolfin.LinearVariationalProblem(a.pyobject, L.pyobject, u.pyobject))
# TODO: add boundary conditions (class is also missing)
# TODO: add parameters

@fenicsclass LinearVariationalSolver
LinearVariationalSolver(problem::LinearVariationalProblem) =
  LinearVariationalSolver(dolfin.LinearVariationalSolver(problem.pyobject))

solve(solver::LinearVariationalSolver) = fenicscall(solver, :solve)
solve(problem::LinearVariationalProblem) = solve(LinearVariationalSolver(problem))
export
  LinearVariationalProblem,
  LinearVariationalSolver, solve


text:
@testset "Predicates" begin
    # ispalindromic
    @test  ispalindromic(dna"")
    @test !ispalindromic(dna"A")
    @test !ispalindromic(dna"C")
    @test  ispalindromic(dna"AT")
    @test  ispalindromic(dna"CG")
    @test !ispalindromic(dna"AC")
    @test !ispalindromic(dna"TT")
    @test  ispalindromic(dna"ANT")
    @test  ispalindromic(dna"ACGT")
    @test !ispalindromic(dna"ACNT")

    @test  ispalindromic(DNAKmer("ACGT"))
    @test !ispalindromic(DNAKmer("CACG"))

    @test  ispalindromic(rna"")
    @test !ispalindromic(rna"A")
    @test !ispalindromic(rna"C")
    @test  ispalindromic(rna"AU")
    @test  ispalindromic(rna"CG")
    @test !ispalindromic(rna"AC")
    @test !ispalindromic(rna"UU")
    @test  ispalindromic(rna"ANU")
    @test  ispalindromic(rna"ACGU")
    @test !ispalindromic(rna"ACNU")

    @test  ispalindromic(RNAKmer("ACGU"))
    @test !ispalindromic(RNAKmer("CACG"))

    @test_throws Exception ispalindromic(aa"PQ")

    # hasambiguity
    @test !hasambiguity(dna"")
    @test !hasambiguity(dna"A")
    @test  hasambiguity(dna"N")
    @test !hasambiguity(dna"ACGT")
    @test  hasambiguity(dna"ANGT")

    @test !hasambiguity(DNAKmer("ACGT"))

    @test !hasambiguity(rna"")
    @test !hasambiguity(rna"A")
    @test  hasambiguity(rna"N")
    @test !hasambiguity(rna"ACGU")
    @test  hasambiguity(rna"ANGU")

    @test !hasambiguity(RNAKmer("ACGU"))

    @test !hasambiguity(aa"")
    @test !hasambiguity(aa"A")
    @test !hasambiguity(aa"P")
    @test  hasambiguity(aa"B")
    @test  hasambiguity(aa"X")
    @test !hasambiguity(aa"ARNDCQEGHILKMFPSTWYVOU")
    @test  hasambiguity(aa"ARXDCQEGHILKMFPSTWYVOU")

    # isrepetitive
    @test  isrepetitive(dna"")
    @test !isrepetitive(dna"", 1)
    @test  isrepetitive(dna"A")
    @test  isrepetitive(dna"A", 1)
    @test  isrepetitive(dna"AAA")
    @test !isrepetitive(dna"ACGT", 2)
    @test  isrepetitive(dna"AAGT", 2)
    @test  isrepetitive(dna"ACCG", 2)
    @test  isrepetitive(dna"ACGG", 2)
    @test !isrepetitive(dna"ACGTCCGT", 3)
    @test  isrepetitive(dna"ACGCCCGT", 3)

    @test !isrepetitive(DNAKmer("ACGT"), 2)
    @test  isrepetitive(DNAKmer("ACCT"), 2)
    @test !isrepetitive(DNAKmer("ACCT"), 3)
    @test  isrepetitive(DNAKmer("ACCCT"), 2)
    @test  isrepetitive(DNAKmer("ACCCT"), 3)

    @test  isrepetitive(rna"")
    @test !isrepetitive(rna"", 1)
    @test  isrepetitive(rna"A")
    @test  isrepetitive(rna"A", 1)
    @test  isrepetitive(rna"AAA")
    @test !isrepetitive(rna"ACGU", 2)
    @test  isrepetitive(rna"AAGU", 2)
    @test  isrepetitive(rna"ACCG", 2)
    @test  isrepetitive(rna"ACGG", 2)
    @test !isrepetitive(rna"ACGUCCGU", 3)
    @test  isrepetitive(rna"ACGCCCGU", 3)

    @test !isrepetitive(RNAKmer("ACGU"), 2)
    @test  isrepetitive(RNAKmer("ACCU"), 2)
    @test !isrepetitive(RNAKmer("ACCU"), 3)
    @test  isrepetitive(RNAKmer("ACCCU"), 2)
    @test  isrepetitive(RNAKmer("ACCCU"), 3)

    @test  isrepetitive(aa"")
    @test !isrepetitive(aa"PGQQ")
    @test  isrepetitive(aa"PGQQ", 2)
    @test !isrepetitive(aa"PPQQ", 3)
    @test  isrepetitive(aa"PPPQQ", 3)
end


text:
include(Pkg.dir("PseudoGenomes", "src","pseudogenome-snp.jl"))
include(Pkg.dir("PseudoGenomes", "src","pseudogenome-deletions.jl"))

mm10path =joinpath(homedir(),"Dropbox/Data/GENOME_DATA/mm10full/mm10.fa")
mm10_masked_fasta = joinpath( tempdir(), "mm10.castmaskdeletion.fasta")

SVs_filename = Pkg.dir( "PseudoGenomes","testdata","28strains.REL-1410-SV.sdp.SAMPLE.tab.gz")
mask_fasta_file_with_strain_SV_deletions(mm10path,mm10_masked_fasta,SVs_filename,"CAST_EiJ")

mm10_masked_and_snp_subst_fasta = joinpath( tempdir(), "mm10.castmaskdeletion.snpsubst.fasta")
substitute_snps_in_fasta_file(mm10_masked_fasta,vcf_path,mm10_masked_and_snp_subst_fasta)


text:
function compute_jacobian(obj::ValueAndJacobianFunction, arg0::Vector{jdouble})
    return jcall(obj, "computeJacobian", RealMatrix, (Vector{jdouble},), arg0)
end

function compute_value(obj::ValueAndJacobianFunction, arg0::Vector{jdouble})
    return jcall(obj, "computeValue", RealVector, (Vector{jdouble},), arg0)
end


text:
function createFakeTiles(room::Maple.Room, x::Int, y::Int, width::Int, height::Int, material::Char='3'; blendIn::Bool=false)::Array{Char, 2}
    tx, ty = floor(Int, x / 8) + 1, floor(Int, y / 8) + 1
    tw, th = floor(Int, width / 8), floor(Int, height / 8)

    fakeTiles = fill('0', (th + 2, tw + 2))

    if blendIn
        fakeTiles[1:end, 1:end] = get(room.fgTiles.data, (ty - 1:ty + th, tx - 1:tx + tw), '0')
    end

    fakeTiles[2:end - 1, 2:end - 1] .= material

    return fakeTiles
end

function getObjectTiles(room::Maple.Room, x::Int, y::Int, width::Int, height::Int)::Maple.ObjectTiles
    tx, ty = floor(Int, x / 8) + 1, floor(Int, y / 8) + 1
    tw, th = floor(Int, width / 8), floor(Int, height / 8)

    return Maple.ObjectTiles(get(room.objTiles.data, (ty - 1:ty + th, tx - 1:tx + tw), -1))
end

function materialTileTypeKey(entity::Maple.Entity)
    if entity.name == "exitBlock" || entity.name == "conditionBlock"
        return "tileType"

    else
        return "tiletype"
    end
end

function validTileEntityTiles(layer::String="fgTiles", ignore::Array{Char, 1}=[Maple.tile_fg_names["Air"], Maple.tile_fg_names["Template"]])
    return filter(m -> !(m in ignore), validTiles(layer))
end

function tiletypeEditingOptions()
    validTiles = Ahorn.validTiles("fgTiles")
    tileNames = Ahorn.tileNames("fgTiles")

    return Dict{String, String}(
        tileNames[mat] => string(mat)
        for mat in validTiles if mat != Maple.tile_fg_names["Air"]
    )
end

# Not the most efficient, but renders correctly
function drawTileEntity(ctx::Cairo.CairoContext, room::Maple.Room, entity::Maple.Entity; material=nothing, alpha=nothing, blendIn::Bool=false)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    blendIn = get(entity.data, "blendin", blendIn)

    if material === nothing
        key = materialTileTypeKey(entity)
        tile = get(entity.data, key, "3")
        material = isa(tile, Number) ? string(tile) : tile
    end

    # Don't draw air versions, even though they shouldn't exist
    if material[1] in validTileEntityTiles()
        fakeTiles = createFakeTiles(room, x, y, width, height, material[1], blendIn=blendIn)
        objTiles = getObjectTiles(room, x, y, width, height)

        drawFakeTiles(ctx, room, fakeTiles, objTiles, true, x, y, alpha=alpha, clipEdges=true)
    end
end

function drawFakeTiles(ctx::Cairo.CairoContext, room::Maple.Room, tiles::Array{Char, 2}, objTiles::Maple.ObjectTiles, fg, x, y; alpha=nothing, clipEdges=false)
    fakeDr = DrawableRoom(
        loadedState.map,
        Maple.Room(
            name="$(room.name)-$x-$y",
            fgTiles=Maple.Tiles(fg ? tiles : Matrix{Char}(undef, 0, 0)),
            bgTiles=Maple.Tiles(!fg ? tiles : Matrix{Char}(undef, 0, 0)),
            objTiles=objTiles
        ),

        TileStates(),
        TileStates(),

        nothing,
        Layer[],

        colors.background_room_fill
    )

    Cairo.save(ctx)

    if clipEdges
        height, width = (size(tiles) .- 2) .* 8
        rectangle(ctx, x, y, width, height)
        clip(ctx)

        # Offset the drawing since we trimmed away the border
        x -= 8
        y -= 8
    end

    translate(ctx, x, y)
    drawTiles(ctx, fakeDr, fg, alpha=alpha, useObjectTiles=true)

    Cairo.restore(ctx)
end

function tileEntityFinalizer(entity::Maple.Entity)
    key = materialTileTypeKey(entity)
    defaultTile = string(Maple.tile_fg_names["Snow"])
    tile = string(get(persistence, "brushes_material_fgTiles", defaultTile))
    tile = tile[1] in Ahorn.validTileEntityTiles() ? tile : defaultTile
    entity.data[key] = tile
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
###############################################################################
#
# # ListofStates tutorial : 1/2
#
# We illustrate here the use of ListofStates in dealing with a warm start
# procedure.
#
# ListofStates can also prove the user history over the iteration process.
#
###############################################################################

using NLPModels, Stopping, Test

# Random search in [0,1] of the global minimum for unconstrained optimization
function algo_rand(stp::NLPStopping)
  x0 = stp.current_state.x
  n = length(x0)
  OK = start!(stp)

  while !OK
    x = rand(n)
    OK = update_and_stop!(stp, x = x, fx = obj(nlp, x), gx = grad(nlp, x))
  end

  return stp
end

include("../test-stopping/rosenbrock.jl")
x0 = 1.5 * ones(6)
nlp = ADNLPModel(rosenbrock, x0, zeros(6), ones(6))

state = NLPAtX(x0)
stop_lstt = NLPStopping(
  nlp,
  state,
  list = ListofStates(state),
  max_iter = 10,
  optimality_check = optim_check_bounded,
)
algo_rand(stop_lstt)
print(stop_lstt.listofstates, print_sym = [:fx, :x])
@test length(stop_lstt.listofstates.list) == 12

#Note the difference if the length of the ListofStates is limited
reinit!(stop_lstt, rstate = true, x = x0)
stop_lstt.listofstates = ListofStates(state, n = 5)
algo_rand(stop_lstt)
print(stop_lstt.listofstates, print_sym = [:fx, :x])
@test length(stop_lstt.listofstates.list) == 5

#Take the best out of 5:
bestfx, best = findmax([stop_lstt.listofstates[i].fx for i = 1:length(stop_lstt.listofstates)])
best_state = copy(stop_lstt.listofstates[best])
reinit!(stop_lstt)
stop_lstt.current_state = best_state
stop_lstt.listofstates = ListofStates(best_state, n = 5)
algo_rand(stop_lstt)
print(stop_lstt.listofstates, print_sym = [:fx, :x])
@test length(stop_lstt.listofstates.list) == 5