text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using Libdl
using ZipFile

include("FMI2_md.jl")

"""
Source: FMISpec2.0.2[p.19]: 2.1.5 Creation, Destruction and Logging of FMU Instances

The mutable struct representing a FMU and all it instantiated instances in the FMI 2.0.2 Standard.
Also contains the paths to the FMU and ZIP folder as well als all the FMI 2.0.2 function pointers
"""
mutable struct FMU2
    modelName::fmi2String
    instanceName::fmi2String
    fmuResourceLocation::fmi2String
    fmuGUID::fmi2String

    modelDescription::fmi2ModelDescription

    fmuType::fmi2Type
    callbackFunctions::fmi2CallbackFunctions
    components::Array{fmi2Component}

    t::Real         # current time
    next_t::Real    # next time

    # paths of ziped and unziped FMU folders
    fmu2Path::fmi2String
    zipPath::fmi2String

    # c-functions
    cInstantiate::Ptr{Cvoid}
    cGetTypesPlatform::Ptr{Cvoid}
    cGetVersion::Ptr{Cvoid}
    cFreeInstance::Ptr{Cvoid}
    cSetDebugLogging::Ptr{Cvoid}
    cSetupExperiment::Ptr{Cvoid}
    cEnterInitializationMode::Ptr{Cvoid}
    cExitInitializationMode::Ptr{Cvoid}
    cTerminate::Ptr{Cvoid}
    cReset::Ptr{Cvoid}
    cGetReal::Ptr{Cvoid}
    cSetReal::Ptr{Cvoid}
    cGetInteger::Ptr{Cvoid}
    cSetInteger::Ptr{Cvoid}
    cGetBoolean::Ptr{Cvoid}
    cSetBoolean::Ptr{Cvoid}
    cGetString::Ptr{Cvoid}
    cSetString::Ptr{Cvoid}
    cGetFMUstate::Ptr{Cvoid}
    cSetFMUstate::Ptr{Cvoid}
    cFreeFMUstate::Ptr{Cvoid}
    cSerializedFMUstateSize::Ptr{Cvoid}
    cSerializeFMUstate::Ptr{Cvoid}
    cDeSerializeFMUstate::Ptr{Cvoid}
    cGetDirectionalDerivative::Ptr{Cvoid}

    # Co Simulation function calls
    cSetRealInputDerivatives::Ptr{Cvoid}
    cGetRealOutputDerivatives::Ptr{Cvoid}
    cDoStep::Ptr{Cvoid}
    cCancelStep::Ptr{Cvoid}
    cGetStatus::Ptr{Cvoid}
    cGetRealStatus::Ptr{Cvoid}
    cGetIntegerStatus::Ptr{Cvoid}
    cGetBooleanStatus::Ptr{Cvoid}
    cGetStringStatus::Ptr{Cvoid}

    # Model Exchange function calls
    cEnterContinuousTimeMode::Ptr{Cvoid}
    cGetContinuousStates::Ptr{Cvoid}
    cGetDerivatives::Ptr{Cvoid}
    cSetTime::Ptr{Cvoid}
    cSetContinuousStates::Ptr{Cvoid}
    cCompletedIntegratorStep::Ptr{Cvoid}
    cEnterEventMode::Ptr{Cvoid}
    cNewDiscreteStates::Ptr{Cvoid}
    cGetEventIndicators::Ptr{Cvoid}
    cGetNominalsOfContinuousStates::Ptr{Cvoid}

    # c-libraries
    libHandle::Ptr{Nothing}

    # sensitivity logging for backward pass (under development)
    sensitivities

    # Constructor
    FMU2() = new()
end

""" struct to handle FMU simulation data / results """
mutable struct fmi2SimulationResult
    valueReferences::Array{fmi2ValueReference}
    dataPoints::Array
    fmu::FMU2

    fmi2SimulationResult() = new()
end

""" collects all data points for variable save index ´i´ inside a fmi2SimulationResult ´sd´ """
function fmi2SimulationResultGetValuesAtIndex(sd::fmi2SimulationResult, i)
    collect(dataPoint[i] for dataPoint in sd.dataPoints)
end

""" collects all time data points inside a fmi2SimulationResult ´sd´ """
function fmi2SimulationResultGetTime(sd::fmi2SimulationResult)
    fmi2SimulationResultGetValuesAtIndex(sd, 1)
end

""" collects all data points for variable with value reference ´tvr´ inside a fmi2SimulationResult ´sd´ """
function fmi2SimulationResultGetValues(sd::fmi2SimulationResult, tvr::fmi2ValueReference)
    @assert tvr != nothing ["fmi2SimulationResultGetValues(...): value referrnce is nothing!"]
    @assert length(sd.dataPoints) > 0 ["fmi2SimulationResultGetValues(...): simulation results are empty!"]

    numVars = length(sd.dataPoints[1])-1
    for i in 1:numVars
        vr = sd.valueReferences[i]
        if vr == tvr
            return fmi2SimulationResultGetValuesAtIndex(sd, i+1)
        end
    end

    nothing
end

""" collects all data points for variable with value name ´s´ inside a fmi2SimulationResult ´sd´ """
function fmi2SimulationResultGetValues(sd::fmi2SimulationResult, s::String)
    fmi2SimulationResultGetValues(sd, fmi2String2ValueReference(sd.fmu, s))
end

"""
Returns the ValueReference coresponding to the variable name
"""
function fmi2String2ValueReference(fmu2::FMU2, name::String)
    reference = nothing
    if haskey(fmu2.modelDescription.stringValueReferences, name)
        reference = fmu2.modelDescription.stringValueReferences[name]
    else
        @warn "no variable with this name found"
    end
    reference
end

"""
Returns an array of ValueReferences coresponding to the variable names
"""
function fmi2String2ValueReference(fmu2::FMU2, names::Array{String})
    vr = Array{fmi2ValueReference}(undef,0)
    for string in names
        reference = fmi2String2ValueReference(fmu2, string)
        if reference == nothing
            @error "valueReference not found"
        else
            push!(vr, reference)
        end
    end
    vr
end
"""
Returns an array of variable names matching a fmi2ValueReference
"""
function fmi2ValueReference2String(fmu2::FMU2, reference::fmi2ValueReference)
    variables = [k for (k,v) in fmu2.modelDescription.stringValueReferences if v == reference]
end

"""
create a copy of the .fmu file as a .zip folder and unzips it
returns the paths to the ziped and ziped folders
"""
function fmi2Unzip(pathTofmu::String)
    # set paths for fmu handling
    zipname = splitext(pathTofmu)
    zipPath = string(zipname[1], ".zip")
    unzipedfmuPath = string(zipname[1], "/")

    if !isdir(unzipedfmuPath)
        cp(pathTofmu, zipPath; force = true)
        mkdir(unzipedfmuPath)

        fileFullPath = isabspath(zipPath) ?  zipPath : joinpath(pwd(),zipPath)
        basePath = dirname(fileFullPath)
        outPath = (unzipedfmuPath == "" ? basePath : (isabspath(unzipedfmuPath) ? unzipedfmuPath : joinpath(pwd(),unzipedfmuPath)))
        isdir(outPath) ? "" : mkdir(outPath)
        zarchive = ZipFile.Reader(fileFullPath)
        for f in zarchive.files
            fullFilePath = joinpath(outPath,f.name)
            if (endswith(f.name,"/") || endswith(f.name,"\\"))
                #mkdir(fullFilePath)
                mkpath(fullFilePath)
            else
                mkpath(dirname(fullFilePath))
                write(fullFilePath, read(f))
            end
        end
        close(zarchive)

        # InfoZIP.unzip(zipPath, unzipedfmuPath)
    end

    @assert isdir(unzipedfmuPath) ["Error:FMU not loaded"]

    (unzipedfmuPath, zipPath)
end

"""
Checks with dlsym for available function in library.
Prints an info text and returns C_NULL if not (soft-check).
"""
function dlsym_opt(libHandle, symbol)
    addr = dlsym(libHandle, symbol; throw_error=false)
    if addr == C_NULL
        @info "This FMU does not support optional function '$symbol'."
    end
    addr
end

"""
sets the properties of the fmu by reading the modelDescription.xml
retrieve all the pointers of binary functions for later @userplot

returns the instance of the fmu struct
"""
function fmi2Load(pathTofmu2::String)
    # Create uninitialized FMU
    fmu_2 = FMU2()
    fmu_2.components = []

    # set paths for fmu handling
    (fmu_2.fmu2Path, fmu_2.zipPath) = fmi2Unzip(pathTofmu2)

    # set paths for modelExchangeScripting and binary
    tmpName = splitpath(fmu_2.fmu2Path)
    fmuName = tmpName[length(tmpName)]
    pathToModelDescription = joinpath(fmu_2.fmu2Path, "modelDescription.xml")

    directoryBinary = ""
    pathToBinary = ""

    if Sys.iswindows()
        directories = ["binaries/win64", "binaries/x86_64-windows"]
        for directory in directories
            directoryBinary = joinpath(fmu_2.fmu2Path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$fmuName.dll")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is Windows, but can't valid find FMU binary."
    elseif Sys.islinux()
        directories = ["binaries/linux64", "binaries/x86_64-linux"]
        for directory in directories
            directoryBinary = joinpath(fmu_2.fmu2Path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$fmuName.so")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is Linux, but can't find valid FMU binary."
    elseif Sys.isapple()
        directories = ["binaries/darwin64", "binaries/x86_64-darwin"]
        for directory in directories
            directoryBinary = joinpath(fmu_2.fmu2Path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$fmuName.dylib")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is macOS, but can't find valid FMU binary."
    else
        @assert false "Unknown target platform."
    end

    # parse modelDescription.xml
    fmu_2.modelDescription = fmi2readModelDescription(pathToModelDescription)
    fmu_2.modelName = fmu_2.modelDescription.modelName
    fmu_2.instanceName = fmu_2.modelDescription.modelName

    lastDirectory = pwd()
    cd(directoryBinary)

    # set FMU binary handler
    fmu_2.libHandle = dlopen(pathToBinary)

    cd(lastDirectory)

    # set fmu properties
    fmu_2.fmuGUID = fmu_2.modelDescription.guid

    if fmu_2.modelDescription.isCoSimulation == fmi2True &&
        fmu_2.modelDescription.isModelExchange == fmi2True
        @info "fmi2Load: FMU supports both CS and ME, using CS as default if nothing specified."
    end

    if fmu_2.modelDescription.isCoSimulation == fmi2True
        fmu_2.fmuType = fmi2CoSimulation::fmi2Type
    elseif fmu_2.modelDescription.isModelExchange == fmi2True
        fmu_2.fmuType = fmi2ModelExchange::fmi2Type
    else
        error(unknownFMUType)
    end

    tmpResourceLocation = string("file:/", fmu_2.fmu2Path)
    tmpResourceLocation = joinpath(tmpResourceLocation, "resources")
    fmu_2.fmuResourceLocation = replace(tmpResourceLocation, "\\" => "/")
    @info "FMU ressource location: $(fmu_2.fmuResourceLocation)"

    # retrieve functions
    fmu_2.cInstantiate                  = dlsym(fmu_2.libHandle, :fmi2Instantiate)
    fmu_2.cGetTypesPlatform             = dlsym(fmu_2.libHandle, :fmi2GetTypesPlatform)
    fmu_2.cGetVersion                   = dlsym(fmu_2.libHandle, :fmi2GetVersion)
    fmu_2.cFreeInstance                 = dlsym(fmu_2.libHandle, :fmi2FreeInstance)
    fmu_2.cSetDebugLogging              = dlsym(fmu_2.libHandle, :fmi2SetDebugLogging)
    fmu_2.cSetupExperiment              = dlsym(fmu_2.libHandle, :fmi2SetupExperiment)
    fmu_2.cEnterInitializationMode      = dlsym(fmu_2.libHandle, :fmi2EnterInitializationMode)
    fmu_2.cExitInitializationMode       = dlsym(fmu_2.libHandle, :fmi2ExitInitializationMode)
    fmu_2.cTerminate                    = dlsym(fmu_2.libHandle, :fmi2Terminate)
    fmu_2.cReset                        = dlsym(fmu_2.libHandle, :fmi2Reset)
    fmu_2.cGetReal                      = dlsym(fmu_2.libHandle, :fmi2GetReal)
    fmu_2.cSetReal                      = dlsym(fmu_2.libHandle, :fmi2SetReal)
    fmu_2.cGetInteger                   = dlsym(fmu_2.libHandle, :fmi2GetInteger)
    fmu_2.cSetInteger                   = dlsym(fmu_2.libHandle, :fmi2SetInteger)
    fmu_2.cGetBoolean                   = dlsym(fmu_2.libHandle, :fmi2GetBoolean)
    fmu_2.cSetBoolean                   = dlsym(fmu_2.libHandle, :fmi2SetBoolean)

    fmu_2.cGetString                    = dlsym_opt(fmu_2.libHandle, :fmi2GetString)
    fmu_2.cSetString                    = dlsym_opt(fmu_2.libHandle, :fmi2SetString)
    fmu_2.cGetFMUstate                  = dlsym_opt(fmu_2.libHandle, :fmi2GetFMUstate)
    fmu_2.cSetFMUstate                  = dlsym_opt(fmu_2.libHandle, :fmi2SetFMUstate)
    fmu_2.cFreeFMUstate                 = dlsym_opt(fmu_2.libHandle, :fmi2FreeFMUstate)
    fmu_2.cSerializedFMUstateSize       = dlsym_opt(fmu_2.libHandle, :fmi2SerializedFMUstateSize)
    fmu_2.cSerializeFMUstate            = dlsym_opt(fmu_2.libHandle, :fmi2SerializeFMUstate)
    fmu_2.cDeSerializeFMUstate          = dlsym_opt(fmu_2.libHandle, :fmi2DeSerializeFMUstate)
    fmu_2.cGetDirectionalDerivative     = dlsym_opt(fmu_2.libHandle, :fmi2GetDirectionalDerivative)

    # CS specific function calls
    if fmu_2.modelDescription.isCoSimulation == fmi2True
        fmu_2.cSetRealInputDerivatives      = dlsym(fmu_2.libHandle, :fmi2SetRealInputDerivatives)
        fmu_2.cGetRealOutputDerivatives     = dlsym(fmu_2.libHandle, :fmi2GetRealOutputDerivatives)
        fmu_2.cDoStep                       = dlsym(fmu_2.libHandle, :fmi2DoStep)
        fmu_2.cCancelStep                   = dlsym(fmu_2.libHandle, :fmi2CancelStep)
        fmu_2.cGetStatus                    = dlsym(fmu_2.libHandle, :fmi2GetStatus)
        fmu_2.cGetRealStatus                = dlsym(fmu_2.libHandle, :fmi2GetRealStatus)
        fmu_2.cGetIntegerStatus             = dlsym(fmu_2.libHandle, :fmi2GetIntegerStatus)
        fmu_2.cGetBooleanStatus             = dlsym(fmu_2.libHandle, :fmi2GetBooleanStatus)
        fmu_2.cGetStringStatus              = dlsym(fmu_2.libHandle, :fmi2GetStringStatus)
    end

    # ME specific function calls
    if fmu_2.modelDescription.isModelExchange == fmi2True
        fmu_2.cEnterContinuousTimeMode      = dlsym(fmu_2.libHandle, :fmi2EnterContinuousTimeMode)
        fmu_2.cGetContinuousStates          = dlsym(fmu_2.libHandle, :fmi2GetContinuousStates)
        fmu_2.cGetDerivatives               = dlsym(fmu_2.libHandle, :fmi2GetDerivatives)
        fmu_2.cSetTime                      = dlsym(fmu_2.libHandle, :fmi2SetTime)
        fmu_2.cSetContinuousStates          = dlsym(fmu_2.libHandle, :fmi2SetContinuousStates)
        fmu_2.cCompletedIntegratorStep      = dlsym(fmu_2.libHandle, :fmi2CompletedIntegratorStep)
        fmu_2.cEnterEventMode               = dlsym(fmu_2.libHandle, :fmi2EnterEventMode)
        fmu_2.cNewDiscreteStates            = dlsym(fmu_2.libHandle, :fmi2NewDiscreteStates)
        fmu_2.cGetEventIndicators           = dlsym(fmu_2.libHandle, :fmi2GetEventIndicators)
        fmu_2.cGetNominalsOfContinuousStates= dlsym(fmu_2.libHandle, :fmi2GetNominalsOfContinuousStates)
    end

    fmu_2
end

"""
Starts a simulation of the fmu instance for the matching fmu type, if both types are available the CoSimulation Simulation is started
"""
function fmi2Simulate(fmu2::FMU2, dt::Real, t_start::Real = 0.0, t_stop::Real = 1.0, recordValues::Union{Array{fmi2ValueReference}, Array{String}} = [])

    if fmu2.fmuType == fmi2CoSimulation::fmi2Type
        fmi2SimulateCS(fmu2, dt, t_start, t_stop, recordValues)
    elseif fmu2.fmuType == fmi2ModelExchange::fmi2Type
        fmi2SimulateME(fmu2, dt, t_start, t_stop)
    else
        error(unknownFMUType)
    end
end

"""
Starts a simulation of the CoSimulation FMU instance
"""
function fmi2SimulateCS(fmu2::FMU2, dt::Real, t_start::Real, t_stop::Real, recordValues::Array{fmi2ValueReference} = [], setup=true)

    if setup
        fmi2SetupExperiment(fmu2, t_start, t_stop)
        fmiEnterInitializationMode(fmu2)
        fmiExitInitializationMode(fmu2)
    end

    t = t_start

    sd = nothing
    record = length(recordValues) > 0

    numDigits = length(string(round(Integer, 1/dt)))

    if record
        sd = fmi2SimulationResult()
        sd.valueReferences = recordValues
        sd.dataPoints = []
        sd.fmu = fmu2

        values = fmi2GetReal(fmu2, sd.valueReferences)
        push!(sd.dataPoints, (t, values...))

        while t < t_stop
            fmiDoStep(fmu2, t, dt)
            t = round(t + dt, digits=numDigits)

            values = fmi2GetReal(fmu2, sd.valueReferences)
            push!(sd.dataPoints, (t, values...))
        end
    else
        while t < t_stop
            fmiDoStep(fmu2, t, dt)
            t = round(t + dt, digits=numDigits)
        end
    end

    sd
end

"""
Starts a simulation of the CoSimulation FMU instance
"""
function fmi2SimulateCS(fmu2::FMU2, dt::Real, t_start::Real, t_stop::Real, recordValues::Array{String}, setup=true)
    vr = fmi2String2ValueReference(fmu2, recordValues)
    fmi2SimulateCS(fmu2, dt, t_start, t_stop, vr, setup)
end

"""
Unload a fmu under the FMI-Standard 2.0.2

Free the allocated memory, close the binaries and remove temporary zip and unziped fmu fmi2ModelDescription
"""
function fmi2Unload(fmu2::FMU2, cleanUp::Bool = true)
    fmi2FreeInstance!(fmu2)

    dlclose(fmu2.libHandle)

    if cleanUp
        try
            rm(fmu2.fmu2Path; recursive = true, force = true)
            rm(fmu2.zipPath; recursive = true, force = true)
        catch e
            @warn "Cannot delete unpacked data on disc."
        end
    end
end

# Comfort functions for fmi2 functions

"""
Returns a string representing the header file used to compile the fmi2 functions.function

Returns "default" by default

For more information call ?fmi2GetTypesPlatform
"""
function fmi2GetTypesPlatform(fmu2::FMU2)
    fmi2GetTypesPlatform(fmu2.cGetTypesPlatform)
end

"""
Returns the version of the FMI Standard used in this FMU

For more information call ?fmi2GetVersion
"""
function fmi2GetVersion(fmu2::FMU2)
    fmi2GetVersion(fmu2.cGetVersion)
end

"""
Create a new instance of the given fmu2, adds a logger if logginOn == true

Returns the instance of the fmu2

For more information call ?fmi2Instantiate
"""
function fmi2Instantiate!(fmu2::FMU2; visible::Bool = false, loggingOn::Bool = false)

    ptrLogger = @cfunction(cbLogger, Cvoid, (Ptr{Cvoid}, Ptr{Cchar}, Cuint, Ptr{Cchar}, Ptr{Cchar}))
    ptrAllocateMemory = @cfunction(cbAllocateMemory, Ptr{Cvoid}, (Csize_t, Csize_t))
    ptrFreeMemory = @cfunction(cbFreeMemory, Cvoid, (Ptr{Cvoid},))
    ptrStepFinished = C_NULL
    fmu2.callbackFunctions = fmi2CallbackFunctions(ptrLogger, ptrAllocateMemory, ptrFreeMemory, ptrStepFinished, C_NULL)

    compAddr = fmi2Instantiate(fmu2.cInstantiate, fmu2.instanceName, fmu2.fmuType, fmu2.fmuGUID, fmu2.fmuResourceLocation, fmu2.callbackFunctions, fmi2Boolean(visible), fmi2Boolean(loggingOn))

    component = fmi2Component(compAddr, fmu2)
    push!(fmu2.components, component)

    component
end

"""
Free the allocated memory used for the looger and fmu2 instance and destroy the instance

For more information call ?fmi2FreeInstance
"""
function fmi2FreeInstance!(fmu2::FMU2)
    for component in fmu2.components
        fmi2FreeInstance(component)
    end

    fmu2.components = []
end

"""
Set the DebugLogger for the FMU

For more information call ?fmi2SetDebugLogging
"""
function fmi2SetDebugLogging(fmu2::FMU2)
    fmi2SetDebugLogging(fmu2.components[end], fmi2False, Unsigned(0), C_NULL)
end

"""
Set the start and end time for a simulation, can only be called after fmi2Instantiate and before fmi2EnterInitializationMode

For more information call ?fmi2SetupExperiment
"""
function fmi2SetupExperiment(fmu2::FMU2,
    toleranceDefined::Bool,
                tolerance::Real,
                startTime::Real,
                stopTimeDefined::Bool,
                stopTime::Real)
    fmu2.t = startTime
    fmi2SetupExperiment(fmu2.components[end], fmi2Boolean(toleranceDefined), fmi2Real(tolerance),
                            fmi2Real(startTime), fmi2Boolean(stopTimeDefined), fmi2Real(stopTime))
end

"""
Setup the simulation but without defining all of the parameters

For more information call ?fmi2SetupExperiment
"""
function fmi2SetupExperiment(fmu2::FMU2, startTime::Real = 0.0, stopTime::Real = startTime; tolerance::Real = 0.0)

    toleranceDefined = (tolerance > 0.0)
    stopTimeDefined = (stopTime > startTime)

    fmi2SetupExperiment(fmu2, toleranceDefined, tolerance, startTime, stopTimeDefined, stopTime)
end

"""
FMU enters Initialization mode

For more information call ?fmi2EnterInitializationMode
"""
function fmi2EnterInitializationMode(fmu2::FMU2)
    fmi2EnterInitializationMode(fmu2.components[end])
end

"""
FMU exits Initialization mode

For more information call ?fmi2ExitInitializationMode
"""
function fmi2ExitInitializationMode(fmu2::FMU2)
    fmi2ExitInitializationMode(fmu2.components[end])
end

"""
Informs FMU that simulation run is terminated

For more information call ?fmi2Terminate
"""
function fmi2Terminate(fmu2::FMU2)
    fmi2Terminate(fmu2.components[end])
end

"""
Reset FMU

For more information call ?fmi2Reset
"""
function fmi2Reset(fmu2::FMU2)
    fmi2Reset(fmu2.components[end])
end

"""
Get the values of an array of fmi2Real variables

For more information call ?fmi2GetReal
"""
function fmi2GetReal(fmu2::FMU2, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    values = zeros(fmi2Real, nvr)
    fmi2GetReal!(fmu2.components[end], vr, nvr, values)
    values
end
"""
Get the value of a fmi2Real variable

For more information call ?fmi2GetReal
"""
function fmi2GetReal(fmu2::FMU2, vr::fmi2ValueReference)
    values = fmi2GetReal(fmu2, [vr])
    values[1]
end
"""
Get the values of an array of fmi2Real variables by variable name

For more information call ?fmi2GetReal
"""
function fmi2GetReal(fmu2::FMU2, vr_string::Array{String})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2GetReal(fmu2, vr)
    end
end
"""
Get the value of a fmi2Real variable by variable name

For more information call ?fmi2GetReal
"""
function fmi2GetReal(fmu2::FMU2, vr_string::String)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2GetReal(fmu2, vr)
    end
end

"""
Get the values of an array of fmi2Real variables

For more information call ?fmi2GetReal
"""
function fmi2GetReal!(fmu2::FMU2, vr::Array{fmi2ValueReference}, values::Array{<:Real})
    vars = zeros(fmi2Real, length(values))
    if length(values) != length(vr)
        @error "Number of value references and in place array doesn't match"
    else
        fmi2GetReal!(fmu2.components[end], vr, Csize_t(length(values)), vars)
    end
    values[:] = vars
end
"""
Get the values of an array of fmi2Real variables by variable name

For more information call ?fmi2GetReal
"""
function fmi2GetReal!(fmu2::FMU2, vr_string::Array{String}, values::Array{<:Real})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    fmi2GetReal!(fmu2, vr, values)
end
"""
Set the values of an array of fmi2Real variables

For more information call ?fmi2SetReal
"""
function fmi2SetReal(fmu2::FMU2, vr::Array{fmi2ValueReference}, value::Array{<:Real})
    nvr = Csize_t(length(vr))
    fmi2SetReal(fmu2.components[end], vr, nvr, Array{fmi2Real}(value))
end

"""
Set the value of a fmi2Real variable

For more information call ?fmi2SetReal
"""
function fmi2SetReal(fmu2::FMU2, vr::fmi2ValueReference, value::Real)
    fmi2SetReal(fmu2, [vr], [value])
end
"""
Set the values of an array of fmi2Real variables by variable name

For more information call ?fmi2SetReal
"""
function fmi2SetReal(fmu2::FMU2, vr_string::Array{String}, value::Array{<:Real})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2SetReal(fmu2, vr, value)
    end
end
"""
Set the value of a fmi2Real variable by variable name

For more information call ?fmi2SetReal
"""
function fmi2SetReal(fmu2::FMU2, vr_string::String, value::Real)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2SetReal(fmu2, vr, value)
    end
end
"""
Get the values of an array of fmi2Integer variables

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(fmu2::FMU2, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    values = zeros(fmi2Integer, nvr)

    fmi2GetInteger!(fmu2.components[end], vr, nvr, values)
    values
end

"""
Get the value of a fmi2Integer variable

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(fmu2::FMU2, vr::fmi2ValueReference)
    values = fmi2GetInteger(fmu2, [vr])
    values[1]
end
"""
Get the values of an array of fmi2Integer variables by variable name

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(fmu2::FMU2, vr_string::Array{String})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2GetInteger(fmu2, vr)
    end
end
"""
Get the value of a fmi2Integer variable by variable name

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger(fmu2::FMU2, vr_string::String)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2GetInteger(fmu2, vr)
    end
end
"""
Get the values of an array of fmi2Integer variables

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger!(fmu2::FMU2, vr::Array{fmi2ValueReference}, values::Array{<:Integer})
    vars = zeros(fmi2Integer, length(values))
    if length(values) != length(vr)
        @error "Number of value references and in place array doesn't match"
    else
        fmi2GetInteger!(fmu2.components[end], vr, Csize_t(length(values)), vars)
    end
    values[:] = vars
end
"""
Get the values of an array of fmi2Integer variables by variable name

For more information call ?fmi2GetInteger
"""
function fmi2GetInteger!(fmu2::FMU2, vr_string::Array{String}, values::Array{<:Integer})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    fmi2GetInteger!(fmu2, vr, values)
end
"""
Set the values of an array of fmi2Integer variables

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(fmu2::FMU2, vr::Array{fmi2ValueReference},value::Array{<:Integer})
    nvr = Csize_t(length(vr))
    fmi2SetInteger(fmu2.components[end], vr, nvr, Array{fmi2Integer}(value))
end
"""
Get the value of a fmi2Integer variable

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(fmu2::FMU2, vr::fmi2ValueReference, value::Integer)
    fmi2SetInteger(fmu2,[vr], [value])
end
"""
Set the values of an array of fmi2Integer variables by variable name

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(fmu2::FMU2, vr_string::Array{String}, value::Array{<:Integer})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2SetInteger(fmu2, vr, value)
    end
end
"""
Set the value of a fmi2Integer variable by variable name

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(fmu2::FMU2, vr_string::String, value::Integer)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2SetInteger(fmu2, vr, value)
    end
end
"""
Get the values of an array of fmi2Boolean variables

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(fmu2::FMU2, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    values = zeros(fmi2Boolean, nvr)
    fmi2GetBoolean!(fmu2.components[end], vr, nvr, values)
    values
end
"""
Get the value of a fmi2Boolean variable

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(fmu2::FMU2, vr::fmi2ValueReference)
    values = fmi2GetBoolean(fmu2, [vr])
    values[1]
end
"""
Get the values of an array of fmi2Boolean variables by variable name

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(fmu2::FMU2, vr_string::Array{String})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2GetBoolean(fmu2, vr)
    end
end
"""
Get the value of a fmi2Boolean variable by variable name

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean(fmu2::FMU2, vr_string::String)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2GetBoolean(fmu2, vr)
    end
end
"""
Get the values of an array of fmi2Boolean variables

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean!(fmu2::FMU2, vr::Array{fmi2ValueReference}, values::Array{Bool})
    vars = zeros(fmi2Boolean, length(values))
    if length(values) != length(vr)
        @error "Number of value references and in place array doesn't match"
    else
        fmi2GetBoolean!(fmu2.components[end], vr, Csize_t(length(values)),vars)
    end
    values[:]  = vars
end
"""
Get the values of an array of fmi2Boolean variables by variable name

For more information call ?fmi2GetBoolean
"""
function fmi2GetBoolean!(fmu2::FMU2, vr_string::Array{String}, values::Array{Bool})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    fmi2GetBoolean!(fmu2, vr, values)
end
"""
Set the values of an array of fmi2Boolean variables

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(fmu2::FMU2, vr::Array{fmi2ValueReference}, value::Array{Bool})
    nvr = Csize_t(length(vr))
    fmi2SetBoolean(fmu2.components[end], vr, nvr, Array{fmi2Boolean}(value))
end
"""
Set the value of a fmi2Boolean variable

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(fmu2::FMU2, vr::fmi2ValueReference, value::Bool)
    fmi2SetBoolean(fmu2, [vr], [value])
end
"""
Set the values of an array of fmi2Boolean variables by variable name

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(fmu2::FMU2, vr_string::Array{String}, value::Array{Bool})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2SetBoolean(fmu2, vr, value)
    end
end
"""
Set the value of a fmi2Boolean variable by variable name

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(fmu2::FMU2, vr_string::String, value::Bool)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2SetBoolean(fmu2, vr, value)
    end
end
"""
Get the values of an array of fmi2String variables

For more information call ?fmi2GetString
"""
function fmi2GetString(fmu2::FMU2, vr::Array{fmi2ValueReference})
    nvr = Csize_t(length(vr))
    value = Vector{Ptr{Cchar}}(undef, nvr)
    fmi2GetString!(fmu2.components[end], vr, nvr, value)
    unsafe_string.(value)
end
"""
Get the value of a fmi2String variable

For more information call ?fmi2GetString
"""
function fmi2GetString(fmu2::FMU2, vr::fmi2ValueReference)
    values = fmi2GetString(fmu2, [vr])
    values[1]
end
"""
Get the values of an array of fmi2String variables by variable name

For more information call ?fmi2GetString
"""
function fmi2GetString(fmu2::FMU2, vr_string::Array{String})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2GetString(fmu2, vr)
    end
end
"""
Get the value of a fmi2String variable by variable name

For more information call ?fmi2GetString
"""
function fmi2GetString(fmu2::FMU2, vr_string::String)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2GetString(fmu2, vr)
    end
end
"""
Get the values of an array of fmi2String variables

For more information call ?fmi2GetString
"""
function fmi2GetString!(fmu2::FMU2, vr::Array{fmi2ValueReference}, values::Array{String})
    vars = Vector{Ptr{Cchar}}(undef, length(vr))
    if length(values) != length(vr)
        @error "Number of value references and in place array doesn't match"
    else
        fmi2GetString!(fmu2.components[end], vr, Csize_t(length(values)), vars)
        values[:] = unsafe_string.(vars)
    end

end
"""
Get the values of an array of fmi2String variables by variable name

For more information call ?fmi2GetString
"""
function fmi2GetString!(fmu2::FMU2, vr_string::Array{String}, values::Array{String})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    fmi2GetString!(fmu2, vr, values)
end
"""
Set the values of an array of fmi2String variables

For more information call ?fmi2SetString
"""
function fmi2SetString(fmu2::FMU2, vr::Array{fmi2ValueReference}, value::Array{String})
    nvr = Csize_t(length(vr))
    ptrs = pointer.(value)
    fmi2SetString(fmu2.components[end], vr, nvr, ptrs)
end
"""
Set the values of a fmi2String variable

For more information call ?fmi2SetString
"""
function fmi2SetString(fmu2::FMU2, vr::fmi2ValueReference, value::String)
    fmi2SetString(fmu2, [vr], [value])
end
"""
Set the values of an array of fmi2String variables by variable name

For more information call ?fmi2SetString
"""
function fmi2SetString(fmu2::FMU2, vr_string::Array{String}, value::Array{String})
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if length(vr) == 0
        @error "no valueReferences could be converted"
    else
        fmi2SetString(fmu2, vr, value)
    end
end
"""
Set the value of a fmi2String variable by variable name

For more information call ?fmi2SetString
"""
function fmi2SetString(fmu2::FMU2, vr_string::String, value::String)
    vr = fmi2String2ValueReference(fmu2, vr_string)
    if vr == nothing
        @error "valueReference not found"
    else
        fmi2SetString(fmu2, vr, value)
    end
end
"""
Get the pointer to the current FMU state

For more information call ?fmi2GetFMUstate
"""
function fmi2GetFMUstate(fmu2::FMU2)
    state = fmi2FMUstate()
    stateRef = Ref(state)
    fmi2GetFMUstate(fmu2.components[end], stateRef)
    state = stateRef[]
    state
end

"""
Set the FMU state to the given fmi2FMUstate

For more information call ?fmi2SetFMUstate
"""
function fmi2SetFMUstate(fmu2::FMU2, state::fmi2FMUstate)
    fmi2SetFMUstate(fmu2.components[end], state)
end

"""
Free the allocated memory for the FMU state

For more information call ?fmi2FreeFMUstate
"""
function fmi2FreeFMUstate(fmu2::FMU2, state::fmi2FMUstate)
    stateRef = Ref(state)
    fmi2FreeFMUstate(fmu2.components[end], stateRef)
    state = stateRef[]
end

"""
Returns the size of a byte vector the FMU can be stored in

For more information call ?fmi2SerzializedFMUstateSize
"""
function fmi2SerializedFMUstateSize(fmu2::FMU2, state::fmi2FMUstate)
    size = 0
    sizeRef = Ref(Csize_t(size))
    fmi2SerializedFMUstateSize(fmu2.components[end], state, sizeRef)
    size = sizeRef[]
end

"""
Serialize the data in the FMU state pointer

For more information call ?fmi2SerzializeFMUstate
"""
function fmi2SerializeFMUstate(fmu2::FMU2, state::fmi2FMUstate)
    size = fmi2SerializedFMUstateSize(fmu2, state)
    serializedState = Array{fmi2Byte}(undef, size)
    fmi2SerializeFMUstate(fmu2.components[end], state, serializedState, size)
    serializedState
end

"""
Deserialize the data in the serializedState fmi2Byte field

For more information call ?fmi2DeSerzializeFMUstate
"""
function fmi2DeSerializeFMUstate(fmu2::FMU2, serializedState::Array{fmi2Byte})
    size = length(serializedState)
    state = fmi2FMUstate()
    stateRef = Ref(state)
    fmi2DeSerializeFMUstate(fmu2.components[end], serializedState, Csize_t(size), stateRef)
    state = stateRef[]
end

"""
Computes directional derivatives

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative(fmu2::FMU2,
                                      vUnknown_ref::Array{fmi2ValueReference},
                                      vKnown_ref::Array{fmi2ValueReference},
                                      dvKnown::Array{fmi2Real} = Array{fmi2Real}([]))

    nKnown = Csize_t(length(vKnown_ref))
    nUnknown = Csize_t(length(vUnknown_ref))

    if length(dvKnown) < nKnown
        dvKnown = ones(fmi2Real, nKnown)
    end

    dvUnknown = zeros(fmi2Real, nUnknown)

    fmi2GetDirectionalDerivative!(fmu2.components[end], vUnknown_ref, nUnknown, vKnown_ref, nKnown, dvKnown, dvUnknown)

    dvUnknown
end
"""
Computes directional derivatives

For more information call ?fmi2GetDirectionalDerivatives
"""
function fmi2GetDirectionalDerivative(fmu::FMU2,
                                      vUnknown_ref::fmi2ValueReference,
                                      vKnown_ref::fmi2ValueReference,
                                      dvKnown::fmi2Real = 1.0,
                                      dvUnknown::fmi2Real = 1.0)

    fmi2GetDirectionalDerivative(fmu, [vUnknown_ref], [vKnown_ref], [dvKnown])[1]
end

# CoSimulation specific functions
"""
Sets the n-th time derivative of real input variables.

For more information call ?fmi2SetRealInputDerivatives
"""
function fmi2SetRealInputDerivatives(fmu2::FMU2, vr::fmi2ValueReference, nvr::Cint, order::Integer, value::Real)
    fmi2SetRealInputDerivatives(fmu2.components[end], vr, nvr, fmi2Integer(order), fmi2Real(value))
end

"""
Retrieves the n-th derivative of output values.

For more information call ?fmi2GetRealOutputDerivatives
"""
function fmi2GetRealOutputDerivatives(fmu2::FMU2, vr::fmi2ValueReference, nvr::Cint, order::Integer, value::Real)
    fmi2GetRealOutputDerivatives(fmu2.components[end], vr, nvr, fmi2Integer(order), fmi2Real(value))
end

"""
The computation of a time step is started.

For more information call ?fmi2DoStep
"""
function fmi2DoStep(fmu2::FMU2, currentCommunicationPoint::Real, communicationStepSize::Real, noSetFMUStatePriorToCurrentPoint::Bool = true)
    fmi2DoStep(fmu2.components[end], fmi2Real(currentCommunicationPoint), fmi2Real(communicationStepSize), fmi2Boolean(noSetFMUStatePriorToCurrentPoint))
end

function fmi2DoStep(fmu2::FMU2, communicationStepSize::Real)
    fmi2DoStep(fmu2.components[end], fmi2Real(fmu2.t), fmi2Real(communicationStepSize), fmi2True)
    fmu2.t += communicationStepSize
end

"""
Cancels a DoStep that didn't finish correctly

For more information call ?fmi2CancelStep
"""
function fmi2CancelStep(fmu2::FMU2)
    fmi2CancelStep(fmu2.components[end])
end

"""
Informs Master of fmi2Status

For more information call ?fmi2GetStatus
"""
function fmi2GetStatus(fmu2::FMU2, s::fmi2StatusKind, value::fmi2Status)
    fmi2GetStatus(fmu2.components[end], s, value)
end

"""
Informs Master of fmi2Status of fmi2Real

For more information call ?fmi2GetRealStatus
"""
function fmi2GetRealStatus(fmu2::FMU2, s::fmi2StatusKind, value::Real)
    fmi2GetRealStatus(fmu2.components[end], s, fmi2Real(value))
end

"""
Informs Master of fmi2Status of fmi2Integer

For more information call ?fmi2GetIntegerStatus
"""
function fmi2GetIntegerStatus(fmu2::FMU2, s::fmi2StatusKind, value::Integer)
    fmi2GetIntegerStatus(fmu2.components[end], s, fmi2Integer(value))
end

"""
Informs Master of fmi2Status of fmi2Boolean

For more information call ?fmi2GetBooleanStatus
"""
function fmi2GetBooleanStatus(fmu2::FMU2, s::fmi2StatusKind, value::Bool)
    fmi2GetBooleanStatus(fmu2.components[end], s, fmiBoolean(value))
end

"""
Informs Master of fmi2Status of fmi2String

For more information call ?fmi2GetStringStatus
"""
function fmi2GetStringStatus(fmu2::FMU2, s::fmi2StatusKind, value::String)
    fmi2GetStringStatus(fmu2.components[end], s, fmiString(value))
end
# Model Exchange specific functions
"""
Set a new time instant

For more information call ?fmi2SetTime
"""
function fmi2SetTime(fmu2::FMU2, time::Real)
    fmi2SetTime(fmu2.components[end], fmi2Real(time))
end

"""
Set a new (continuous) state vector

For more information call ?fmi2SetContinuousStates
"""
function fmi2SetContinuousStates(fmu2::FMU2, x::Union{Array{Float32}, Array{Float64}})
    nx = Csize_t(length(x))
    fmi2SetContinuousStates(fmu2.components[end], Array{fmi2Real}(x), nx)
end

"""
The model enters Event Mode

For more information call ?fmi2EnterEventMode
"""
function fmi2EnterEventMode(fmu2::FMU2)
    fmi2EnterEventMode(fmu2.components[end])
end
"""
Returns the next discrete states

For more information call ?fmi2NewDiscreteStates
"""
function fmi2NewDiscreteStates(fmu2::FMU2)
    eventInfo = fmi2EventInfo()
    fmi2NewDiscreteStates(fmu2.components[end], eventInfo)
    eventInfo
end

"""
The model enters Continuous-Time Mode

For more information call ?fmi2EnterContinuousTimeMode
"""
function fmi2EnterContinuousTimeMode(fmu2::FMU2)
    fmi2EnterContinuousTimeMode(fmu2.components[end])
end

"""
This function must be called by the environment after every completed step

For more information call ?fmi2CompletedIntegratorStep
"""
function fmi2CompletedIntegratorStep(fmu2::FMU2,
                                     noSetFMUStatePriorToCurrentPoint::fmi2Boolean)
    enterEventMode = fmi2Boolean(false)
    terminateSimulation = fmi2Boolean(false)
    status = fmi2CompletedIntegratorStep!(fmu2.components[end],
                                         noSetFMUStatePriorToCurrentPoint,
                                         enterEventMode,
                                         terminateSimulation)
    (status, enterEventMode, terminateSimulation)
end

"""
Compute state derivatives at the current time instant and for the current states.

For more information call ?fmi2GetDerivatives
"""
function  fmi2GetDerivatives(fmu2::FMU2)
    nx = Csize_t(fmu2.modelDescription.numberOfContinuousStates)
    derivatives = zeros(fmi2Real, nx)
    fmi2GetDerivatives(fmu2.components[end], derivatives, nx)
    derivatives
end
"""
Returns the event indicators of the FMU

For more information call ?fmi2GetEventIndicators
"""
function fmi2GetEventIndicators(fmu2::FMU2)
    ni = Csize_t(fmu2.modelDescription.numberOfEventIndicators)
    eventIndicators = zeros(fmi2Real, ni)
    fmi2GetEventIndicators(fmu2.components[end], eventIndicators, ni)
    eventIndicators
end
"""
Return the new (continuous) state vector x

For more information call ?fmi2GetContinuousStates
"""
function fmi2GetContinuousStates(fmu2::FMU2)
    nx = Csize_t(fmu2.modelDescription.numberOfContinuousStates)
    x = zeros(fmi2Real, nx)
    fmi2GetContinuousStates(fmu2.components[end], x, nx)
    x
end

"""
Return the new (continuous) state vector x

For more information call ?fmi2GetNominalsOfContinuousStates
"""
function fmi2GetNominalsOfContinuousStates(fmu2::FMU2)
    nx = Csize_t(fmu2.modelDescription.numberOfContinuousStates)
    x = zeros(fmi2Real, nx)
    fmi2GetNominalsOfContinuousStates(fmu2.components[end], x, nx)
    x
end


text:
## Utility function

using Base: check_reducedims

# Determine if there are any true values in a BitArray in a given
# range. We use this for reductions with skipmissing=false along the first
# dimension.
function _any(B::BitArray, istart::Int, iend::Int)
    chunks = B.chunks
    startidx, startbit = Base.get_chunks_id(istart)
    endidx, endbit = Base.get_chunks_id(iend)
    startidx == endidx && return chunks[startidx] >> startbit << (63-endbit+startbit) != 0

    chunks[startidx] >> startbit != 0 && return true
    for i = startidx+1:endidx-1
        chunks[i] != 0 && return true
    end
    chunks[endidx] << (63-endbit) != 0 && return true
    return false
end

# Counts the number of ones in a given range. We use this for counting
# the values for mean and var with skipmissing=false along the first
# dimension.
function _count(B::BitArray, istart::Int, iend::Int)
    chunks = B.chunks
    startidx, startbit = Base.get_chunks_id(istart)
    endidx, endbit = Base.get_chunks_id(iend)
    startidx == endidx && return count_ones(chunks[startidx] >> startbit << (63-endbit+startbit))

    n = 0
    n += count_ones(chunks[startidx] >> startbit)
    for i = startidx+1:endidx-1
        n += count_ones(chunks[i])
    end
    n += count_ones(chunks[endidx] << (63-endbit))
    return n
end

## missing-preserving
@generated function _mapreducedim!(f::SafeMapFuns, op::SafeReduceFuns,
                                   R::DataArray, A::DataArray{T,N} where {T}) where {N}
    quote
        data = A.data
        na = A.na

        lsiz = check_reducedims(R, data)
        isempty(data) && return R

        if lsiz > 16
            # use mapreduce_impl, which is probably better tuned to achieve higher performance
            nslices = div(length(A), lsiz)
            ibase = 0
            extr = daextract(R)
            for i = 1:nslices
                if _any(na, ibase+1, ibase+lsiz)
                    unsafe_setmissing!(R, extr, i)
                else
                    v = Base.mapreduce_impl(f, op, data, ibase+1, ibase+lsiz)
                    @inbounds unsafe_dasetindex!(R, extr, v, i)
                end
                ibase += lsiz
            end
        else
            @nextract $N sizeR d->size(R,d)
            na_chunks = A.na.chunks

            new_data = R.data

            # If reducing to a DataArray, skip strides with missings.
            # In my benchmarks, it is actually faster to compute a new missing
            # array and BitArray it than to operate on the BitArray
            # directly.
            new_na = fill(false, size(new_data))

            @nexprs 1 d->(state_0 = state_{$N} = 1)
            @nexprs $N d->(skip_d = sizeR_d == 1)
            k = 1
            @nloops($N, i, A,
                d->(state_{d-1} = state_d),
                d->(skip_d || (state_d = state_0)), begin
                    @inbounds vna = new_na[state_0] | Base.unsafe_bitgetindex(na_chunks, k)
                    if vna
                        @inbounds new_na[state_0] = true
                    else
                        @inbounds x = data[k]
                        v = f(x)
                        @inbounds v0 = new_data[state_0]
                        nv = op(v0, v)
                        @inbounds new_data[state_0] = nv
                    end

                    state_0 += 1
                    k += 1
                end)

            R.na = BitArray(new_na)
        end
        return R
    end
end

## missing-preserving to array
@generated function _mapreducedim!(f::SafeMapFuns, op::SafeReduceFuns,
                                   R::AbstractArray, A::DataArray{T,N} where {T}) where {N}
    quote
        data = A.data
        na   = A.na

        lsiz = check_reducedims(R, data)
        isempty(data) && return R

        if lsiz > 16
            # use mapreduce_impl, which is probably better tuned to achieve higher performance
            nslices = div(length(A), lsiz)
            ibase = 0
            extr = daextract(R)
            for i = 1:nslices
                if _any(na, ibase+1, ibase+lsiz)
                    throw(MissingException("array contains missing values but output array does not support them"))
                else
                    v = Base.mapreduce_impl(f, op, data, ibase+1, ibase+lsiz)
                    @inbounds unsafe_dasetindex!(R, extr, v, i)
                end
                ibase += lsiz
            end
        else
            @nextract $N sizeR d->size(R,d)

            # If reducing to a non-DataArray, throw an error at the start on missing
            any(ismissing, A) && throw(MissingException("array contains missing values: pass skipmissing=true to skip them"))
            @nloops $N i data d->(j_d = sizeR_d==1 ? 1 : i_d) begin
                @inbounds x = (@nref $N data i)
                v = f(x)
                @inbounds v0 = (@nref $N R j)
                nv = op(v0, v)
                @inbounds (@nref $N R j) = nv
            end
        end
        return R
    end
end
_mapreducedim!(f, op, R, A) = Base._mapreducedim!(f, op, R, A)

## missing-skipping
_getdata(A) = A
_getdata(A::DataArray) = A.data

# mapreduce across a dimension. If specified, C contains the number of
# non-missing values reduced into each element of R.
@generated function _mapreducedim_skipmissing_impl!(f, op, R::AbstractArray,
                                                        C::Union{Array{Int}, Nothing},
                                                        A::DataArray{T,N} where {T}) where {N}
    quote

        data = A.data
        na = A.na
        na_chunks = na.chunks
        new_data = _getdata(R)

        C === nothing || size(R) == size(C) || throw(DimensionMismatch("R and C must have same size"))
        lsiz = check_reducedims(new_data, data)
        isempty(data) && return R

        if lsiz > 16
            # keep the accumulator as a local variable when reducing along the first dimension
            nslices = div(length(A), lsiz)
            ibase = 0
            for i = 1:nslices
                # TODO: use pairwise impl for sum
                @inbounds v = new_data[i]
                @inbounds C !== nothing && (C[i] = lsiz - _count(na, ibase+1, ibase+lsiz))
                for k = ibase+1:ibase+lsiz
                    @inbounds Base.unsafe_bitgetindex(na_chunks, k) && continue
                    @inbounds x = data[k]
                    v = convert(typeof(v), op(f(x), v))::typeof(v)
                end
                @inbounds new_data[i] = v
                ibase += lsiz
            end
        else
            # general implementation
            @nextract $N sizeR d->size(new_data,d)
            @nexprs 1 d->(state_0 = state_{$N} = 1)
            @nexprs $N d->(skip_d = sizeR_d == 1)
            k = 1
            C !== nothing && fill!(C, div(length(A), length(R)))
            @nloops($N, i, A,
                d->(state_{d-1} = state_d),
                d->(skip_d || (state_d = state_0)), begin
                    @inbounds xna = Base.unsafe_bitgetindex(na_chunks, k)
                    if xna
                        C !== nothing && @inbounds C[state_0] -= 1
                    else
                        @inbounds x = data[k]
                        v = f(x)
                        @inbounds v0 = new_data[state_0]
                        nv = op(v0, v)
                        @inbounds new_data[state_0] = nv
                    end

                    state_0 += 1
                    k += 1
            end)
        end
        return R
    end
end

_mapreducedim_skipmissing!(f, op, R::AbstractArray, A::DataArray) =
    _mapreducedim_skipmissing_impl!(f, op, R, nothing, A)

# for MinFun/MaxFun, min or max is missing if all values along a dimension are missing
function _mapreducedim_skipmissing!(f, op::Union{typeof(min), typeof(max)}, R::DataArray, A::DataArray)
    R.na = BitArray(all!(fill(true, size(R)), A.na))
    _mapreducedim_skipmissing_impl!(f, op, R, nothing, A)
end
function _mapreducedim_skipmissing!(f, op::Union{typeof(min), typeof(max)}, R::AbstractArray, A::DataArray)
    if any(all!(fill(true, size(R)), A.na))
        throw(MissingException("some dimensions of the array only contain missing values"))
    end
    _mapreducedim_skipmissing_impl!(f, op, R, nothing, A)
end

## general reducedim interface

for op in (+, *, &, |, min, max)
    @eval begin
        function Base.initarray!(a::DataArray{T}, op::typeof($op), init::Bool) where T
            if init
                Base.initarray!(a.data, op, true)
                Base.fill!(a.na, false)
            end
            a
        end
    end
end

function Base.reducedim_initarray(A::DataArray, region, v0, ::Type{R}) where R
    rd = length.(Base.reduced_indices(A.data, region))
    DataArray(fill!(similar(A.data, R, rd), v0), falses(rd))
end
function Base.reducedim_initarray0(A::DataArray, region, v0, ::Type{R}) where R
    rd = length.(Base.reduced_indices0(A,region))
    DataArray(fill!(similar(A.data, R, rd), v0), falses(rd))
end

function Base.mapreducedim!(f::Function, op, R::AbstractArray, A::DataArray;
                            skipmissing::Bool=false, skipna::Bool=false)
    if skipna && !skipmissing
        Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
        skipmissing = true
    end
    skipmissing ? _mapreducedim_skipmissing!(f, op, R, A) : _mapreducedim!(f, op, R, A)
end
function Base.mapreducedim!(f, op, R::AbstractArray, A::DataArray;
                            skipmissing::Bool=false, skipna::Bool=false)
    if skipn && !skipmissing
        Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
        skipmissing = true
    end
    skipmissing ? _mapreducedim_skipna!(f, op, R, A) : _mapreducedim!(f, op, R, A)
end
Base.reducedim!(op, R::DataArray{RT}, A::AbstractArray;
                skipmissing::Bool=false, skipna::Bool=false) where {RT} =
    Base.mapreducedim!(identity, op, R, A, zero(RT); skipna=skipna)

Base.mapreducedim(f, op, A::DataArray, region, v0;
                  skipmissing::Bool=false, skipna::Bool=false) =
    Base.mapreducedim!(f, op, Base.reducedim_initarray(A, region, v0), A;
                       skipmissing=skipmissing, skipna=skipna)
Base.mapreducedim(f, op, A::DataArray{T}, region;
                  skipmissing::Bool=false, skipna::Bool=false) where {T} =
    Base.mapreducedim!(f, op, Base.reducedim_init(f, op, A, region), A;
                       skipmissing=skipmissing, skipna=skipna)

Base.reducedim(op, A::DataArray, region, v0; skipna::Bool=false) =
    Base.mapreducedim(identity, op, A, region, v0; skipmissing=skipmissing, skipna=skipna)
Base.reducedim(op, A::DataArray, region; skipna::Bool=false) =
    Base.mapreducedim(identity, op, A, region; skipmissing=skipmissing, skipna=skipna)

## usual reductions

for (basfn, Op) in [(:sum, +), (:prod, *),
                    (:maximum, max), (:minimum, min),
                    (:all, &), (:any, |)]
    fname = Expr(:., :Base, Base.Meta.quot(basfn))
    fname! = Expr(:., :Base, Base.Meta.quot(Symbol(string(basfn, '!'))))
    @eval begin
        $(fname!)(f::Union{Function,$(supertype(typeof(abs)))}, r::AbstractArray, A::DataArray;
                  init::Bool=true, skipmissing::Bool=false, skipna::Bool=false) =
            Base.mapreducedim!(f, $(Op), Base.initarray!(r, $(Op), init), A;
                               skipmissing=skipmissing, skipna=skipna)
        $(fname!)(r::AbstractArray, A::DataArray;
                  init::Bool=true, skipmissing::Bool=false, skipna::Bool=false) =
            $(fname!)(identity, r, A; init=init, skipmissing=skipmissing, skipna=skipna)

        $(fname)(f::Union{Function,$(supertype(typeof(abs)))}, A::DataArray, region;
                 skipmissing::Bool=false, skipna::Bool=false) =
            Base.mapreducedim(f, $(Op), A, region; skipmissing=skipmissing, skipna=skipna)
        $(fname)(A::DataArray, region; skipmissing::Bool=false, skipna::Bool=false) =
            $(fname)(identity, A, region; skipmissing=skipmissing, skipna=skipna)
    end
end

## mean

function Base.mean!(R::AbstractArray{T}, A::DataArray;
                    skipmissing::Bool=false, skipna::Bool=false, init::Bool=true) where {T}
    init && fill!(R, 0)
    if skipna || skipmissing
        C = Array{Int}(size(R))
        _mapreducedim_skipmissing_impl!(identity, +, R, C, A)
        broadcast!(/, R, R, C)
    else
        sum!(R, A; skipmissing=false)
        broadcast!(/, R, R, convert(T, length(A)/length(R)))
        R
    end
end

Base.mean(A::DataArray{T}, region; skipmissing::Bool=false, skipna::Bool=false) where {T} =
    mean!(Base.reducedim_initarray(A, region, zero(Base.momenttype(T))), A;
          skipmissing=skipmissing, skipna=skipna, init=false)

## var

struct MapReduceDim2ArgHelperFun{F,T}
    f::F
    val::T
end
(f::MapReduceDim2ArgHelperFun)(x) = f.f(x, f.val)

# A version of _mapreducedim! that accepts an array S of the same size
# as R, the elements of which are passed as a second argument to f.
@generated function _mapreducedim_2arg!(f, op, R::DataArray,
                                                    A::DataArray{T,N} where {T},
                                                    S::AbstractArray) where {N}
    quote
        data = A.data
        na = A.na
        Sextr = daextract(S)

        lsiz = check_reducedims(R, data)
        size(R) == size(S) || throw(DimensionMismatch("R and S must have same size"))
        isempty(data) && return R

        if lsiz > 16
            # use mapreduce_impl, which is probably better tuned to achieve higher performance
            nslices = div(length(A), lsiz)
            ibase = 0
            extr = daextract(R)
            for i = 1:nslices
                if unsafe_ismissing(S, Sextr, i) || _any(na, ibase+1, ibase+lsiz)
                    unsafe_setmissing!(R, extr, i)
                else
                    @inbounds s = unsafe_getindex_notmissing(S, Sextr, i)
                    v = Base.mapreduce_impl(MapReduceDim2ArgHelperFun(f, s), op, data, ibase+1, ibase+lsiz)
                    @inbounds unsafe_dasetindex!(R, extr, v, i)
                end
                ibase += lsiz
            end
        else
            @nextract $N sizeR d->size(R,d)
            na_chunks = A.na.chunks
            new_data = R.data
            new_na = isa(S, DataArray) ? Array(S.na) : fill(false, size(S))

            @nexprs 1 d->(state_0 = state_{$N} = 1)
            @nexprs $N d->(skip_d = sizeR_d == 1)
            k = 1
            @nloops($N, i, A,
                d->(state_{d-1} = state_d),
                d->(skip_d || (state_d = state_0)), begin
                    @inbounds vna = new_na[state_0] | Base.unsafe_bitgetindex(na_chunks, k)
                    if vna
                        @inbounds new_na[state_0] = true
                    else
                        @inbounds s = unsafe_getindex_notmissing(S, Sextr, state_0)
                        @inbounds x = data[k]
                        v = f(x, s)
                        @inbounds v0 = new_data[state_0]
                        nv = op(v0, v)
                        @inbounds new_data[state_0] = nv
                    end

                    state_0 += 1
                    k += 1
                end)

            R.na = BitArray(new_na)
        end
        return R
    end
end

# A version of _mapreducedim_skipmissing! that accepts an array S of the same size
# as R, the elements of which are passed as a second argument to f.
@generated function _mapreducedim_skipmissing_2arg!(f, op, R::AbstractArray,
                                                        C::Union{Array{Int}, Nothing},
                                                        A::DataArray{T,N} where {T}, S::AbstractArray) where {N}
    quote
        data = A.data
        na = A.na
        na_chunks = na.chunks
        new_data = _getdata(R)
        Sextr = daextract(S)

        lsiz = check_reducedims(new_data, data)
        C === nothing || size(R) == size(C) || throw(DimensionMismatch("R and C must have same size"))
        size(R) == size(S) || throw(DimensionMismatch("R and S must have same size"))
        isempty(data) && return R
        @nextract $N sizeR d->size(new_data,d)
        sizA1 = size(data, 1)

        # If there are any missings in S, assume these will produce missings in R
        if isa(S, DataArray)
            copy!(R.na, S.na)
        end

        if lsiz > 16
            # keep the accumulator as a local variable when reducing along the first dimension
            nslices = div(length(A), lsiz)
            ibase = 0
            for i = 1:nslices
                @inbounds v = new_data[i]
                !isa(C, Nothing) && (C[i] = lsiz - _count(na, ibase+1, ibase+lsiz))

                # If S[i] is missing, skip this iteration
                @inbounds sna = unsafe_ismissing(S, Sextr, i)
                if !sna
                    @inbounds s = unsafe_getindex_notmissing(S, Sextr, i)
                    # TODO: use pairwise impl for sum
                    for k = ibase+1:ibase+lsiz
                        @inbounds Base.unsafe_bitgetindex(na_chunks, k) && continue
                        @inbounds x = data[k]
                        v = convert(typeof(v), op(f(x, s), v))::typeof(v)
                    end

                    @inbounds new_data[i] = v
                end

                ibase += lsiz
            end
        else
            # general implementation
            @nexprs 1 d->(state_0 = state_{$N} = 1)
            @nexprs $N d->(skip_d = sizeR_d == 1)
            k = 1
            !isa(C, Nothing) && fill!(C, div(length(A), length(R)))
            @nloops($N, i, A,
                d->(state_{d-1} = state_d),
                d->(skip_d || (state_d = state_0)), begin
                    @inbounds xna = Base.unsafe_bitgetindex(na_chunks, k) | unsafe_ismissing(S, Sextr, state_0)
                    if xna
                        !isa(C, Nothing) && @inbounds C[state_0] -= 1
                    else
                        @inbounds s = unsafe_getindex_notmissing(S, Sextr, state_0)
                        @inbounds x = data[k]
                        v = f(x, s)
                        @inbounds v0 = new_data[state_0]
                        nv = op(v0, v)
                        @inbounds new_data[state_0] = nv
                    end

                    state_0 += 1
                    k += 1
            end)
        end
        return R
    end
end

struct Abs2MinusFun end
(::Abs2MinusFun)(x, m) = abs2(x - m)

function Base.varm!(R::AbstractArray, A::DataArray, m::AbstractArray; corrected::Bool=true,
                    skipmissing::Bool=false, skipna::Bool=false, init::Bool=true)
    if skipna && !skipmissing
            Base.depwarn("skipna=$skipna is deprecated, use skipmissing=$skipna instead", :mapreduce)
            skipmissing = true
    end
    if isempty(A)
        fill!(R, convert(eltype(R), NaN))
    else
        init && fill!(R, zero(eltype(R)))
        if skipmissing
            C = Array{Int}(size(R))

            # Compute R = abs2(A-m)
            _mapreducedim_skipmissing_2arg!(Abs2MinusFun(), +, R, C, A, m)

            # Divide by number of non-missing values
            if corrected
                for i = 1:length(C)
                    @inbounds C[i] = max(C[i] - 1, 0)
                end
            end
            broadcast!(/, R, R, C)
        else
            # Compute R = abs2(A-m)
            _mapreducedim_2arg!(Abs2MinusFun(), +, R, A, m)

            # Divide by number of values
            broadcast!(/, R, R, div(length(A), length(R)) - corrected)
        end
    end
end

Base.varm(A::DataArray{T}, m::AbstractArray, region; corrected::Bool=true,
          skipmissing::Bool=false, skipna::Bool=false) where {T} =
    Base.varm!(Base.reducedim_initarray(A, region, zero(Base.momenttype(T))), A, m;
               corrected=corrected, skipmissing=skipmissing, skipna=skipna, init=false)

function Base.var(A::DataArray{T}, region::Union{Integer, AbstractArray, Tuple};
                  corrected::Bool=true, mean=nothing,
                  skipmissing::Bool=false, skipna::Bool=false) where T
    if mean == 0
        Base.varm(A, Base.reducedim_initarray(A, region, zero(Base.momenttype(T))), region;
                  corrected=corrected, skipmissing=skipmissing, skipna=skipna)
    elseif mean == nothing
        if skipna || skipmissing
            # Can reduce mean into ordinary array
            m = zeros(Base.momenttype(T), length.(Base.reduced_indices(A, region)))
            Base.varm(A, Base.mean!(m, A; skipmissing=skipmissing, skipna=skipna), region;
                      corrected=corrected, skipmissing=skipmissing, skipna=skipna)
        else
            Base.varm(A, Base.mean(A, region; skipmissing=skipmissing, skipna=skipna), region;
                      corrected=corrected, skipmissing=skipmissing, skipna=skipna)
        end
    elseif isa(mean, AbstractArray)
        Base.varm(A, mean::AbstractArray, region;
                  corrected=corrected, skipmissing=skipmissing, skipna=skipna)
    else
        throw(ErrorException("invalid value of mean"))
    end
end

ext: jl


text:
"""
Prints a list of logical combos, variables defined, number of feasible outcomes,
number of commands, as well as the last command.

discover(x::LogicalCombo)
discover()

"""
function discover(x::LogicalCombo)
   nameset = names(Main)

   logicalcombonames = String[]
   variables         = String[]
   nfeasiblelist     = Int16[]
   ncommands         = Int16[]
   lastcommand       = String[]

   for v in nameset
     # global logicalcombonames, variables, nfeasible, ncommands
     x = getfield(Main, v)
     if isa(x, LogicalCombo)
       push!(logicalcombonames, string(v))
       push!(variables, join(keys(x),", "))
       push!(nfeasiblelist, nfeasible(x))
       push!(ncommands, length(x.commands))
       (length(x.commands) > 0 ? push!(lastcommand, last(x.commands)) : push!(lastcommand, ""))
     end
   end

   (length(logicalcombonames)==0) && (println("No LogicalCombos found in Main"); return)

   txtout = "| Name | Variables | #feasible | #commands | Last Command | \n"
   txtout *= "| :---: | :---: | :---: | :---: | :---: | \n"

   for i in 1:length(logicalcombonames)
       txtout *= "| $(logicalcombonames[i]) | $(variables[i]) | " *
             "$(nfeasiblelist[i]) | $(ncommands[i]) | $(lastcommand[i]) | \n"
   end

   printmarkdown(txtout)
end

discover() = discover(LogicalCombo())


text:
#!/usr/bin/env julia
################################################################################
#
#  Copyright (c) 2015 Wojciech Migda
#  All rights reserved
#  Distributed under the terms of the Apache 2.0 license
#
################################################################################
#
#  Filename: Worker.jl
#
#  Decription:
#       Generate MFCC feature vectors
#
#  Authors:
#       Wojciech Migda
#
################################################################################
#
#  History:
#  --------
#  Date         Who  Ticket     Description
#  ----------   ---  ---------  ------------------------------------------------
#  2015-09-20   wm              Initial version
#
################################################################################

module Worker

export reduce

using MFCC

const THIS_DIR = dirname(Base.source_path())

@windows? (
begin
    ENV["PATH"] = "$(ENV["PATH"]);$(THIS_DIR);$(THIS_DIR)/../external/bin"
    #println(ENV["PATH"])
    const LIB_MP3_DECODER = "libmp3decoder"
end
:
begin
    const LIB_MP3_DECODER = find_library(["libmp3decoder"], [THIS_DIR, "$(THIS_DIR)/../external/bin"])
end
)

type Mp3Params
    rate::Uint32
    channels::Uint32
    encoding::Uint32
end

function mp3decoder!(ifname::String, osig::AbstractVector{Int16})
    # we expect 10 secs of samples, int16, stereo at 44100 Hz
    mp3params = Mp3Params(uint32(0), uint32(0), uint32(0))

    const NREAD =
        ccall((:mp3decoder, LIB_MP3_DECODER),
            Csize_t, (Ptr{Uint8}, Ptr{Void}, Csize_t, Ptr{Mp3Params}),
            ifname, osig, size(osig, 1) * sizeof(osig[1]), &mp3params)
    @assert mp3params.encoding == 0xD0 # MPG123_ENC_SIGNED_16

    return NREAD, mp3params
end

function downsample{_iTp, _oTp<:FloatingPoint}(
    # code_native(downsample, (Int, Bool, Vector{Int}, Vector{Float64}))
    ratio::Int,
    tomono::Bool,
    ivec::AbstractVector{_iTp},
    ovec::AbstractVector{_oTp})

    const STRIDE::Int = ratio * (tomono ? 2 : 1)
    const ILEN::Int = size(ivec, 1)

    resize!(ovec, div(ILEN, STRIDE))

    opos::Int = 1
    ipos::Int = 1
    @simd for ipos = 1:STRIDE:ILEN - STRIDE + 1
        @inbounds ovec[opos] = 0
        @simd for sub = 0:STRIDE - 1
            @inbounds ovec[opos] += ivec[ipos + sub]
        end
        @inbounds ovec[opos] /= STRIDE

        opos += 1
    end
end

function gen_mfcc{_Tp}(MP3_DIR::String, TRAIN_X::AbstractVector{_Tp}, start::Int;
                       _winlen=0.025, _winstep=0.01, _numcep=13,
                       _nfilt=26, _nfft=512, _preemph=0.97, _ceplifter=22, _appendEnergy=true)

    const DATA_DIR = "$(THIS_DIR)/../../data"
    const TRAIN_MP3 = MP3_DIR

    const NSAMP = 10 * 44100
    const NCHAN = 2
    full_sig = Array(Int16, NSAMP * NCHAN)
    sig = similar(full_sig, Float64)

    vvX = Vector{Float64}[]

    for irow = 1:size(TRAIN_X, 1)
        const NREAD, MP3PARAMS = mp3decoder!("$(TRAIN_MP3)/$(TRAIN_X[irow])", full_sig)
        @inbounds full_sig[1 + NREAD / 2:end] = 0

        const DOWN_RATIO = 4
        downsample(DOWN_RATIO, true, full_sig, sig)

        mfcc_feat = mfcc(sig, samplerate=MP3PARAMS.rate / float(DOWN_RATIO),
                         winlen=_winlen, winstep=_winstep, numcep=_numcep, nfilt=_nfilt, nfft=_nfft, preemph=_preemph, ceplifter=_ceplifter, appendEnergy=_appendEnergy)

        push!(vvX, vec(mfcc_feat))
    end

    return (vvX, start)
end

function reduce(
    lhs::(Vector{Vector{Float64}}, Int),
    rhs::(Vector{Vector{Float64}}, Int)
    )
    return [lhs, rhs]
end

function reduce(
    lhs::Vector{(Vector{Vector{Float64}}, Int)},
    rhs::(Vector{Vector{Float64}}, Int)
    )
    push!(lhs, rhs)
    return lhs
end

function reduce(
    lhs::Vector{(Vector{Vector{Float64}}, Int)},
    rhs::Vector{(Vector{Vector{Float64}}, Int)}
    )
    append!(lhs, rhs)
    return lhs
end

end # module Worker


text:
# SINGLE NODE/CORE VERSION

import Pkg; Pkg.instantiate()

include("/users/yh31/scratch/projects/gigaword_64k/src/gigaword_64k.jl")

const path_xin_dir = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga5/data/xin_eng"

const path_intermed_data = "/users/yh31/scratch/datasets/entity_linking/raw_data/gigaword/giga_topk/intermediate_files"
const path_intermed_xin = joinpath(path_intermed_data, "xin")

# const path_of_tree = path_xin_dir
# const path_output = path_intermed_xin
# const n_items = 8

process_part_of_tree(path_xin_dir, path_intermed_xin, 4)


text:
@comment """
# AbstractDataFrame
"""

"""

An abstract type for which all concrete types expose a database-like
interface.

### Common methods

An AbstractDataFrame is a two-dimensional table with Symbols for
column names. An AbstractDataFrame is also similar to an Associative
type in that it allows indexing by a key (the columns).

The following are normally implemented for AbstractDataFrames:

* `describe(d)` : summarize columns
* `dump(d)` : show structure
* `hcat(d1, d2)` : horizontal concatenation
* `vcat(d1, d2)` : vertical concatenation
* `names(d)` : columns names
* `names!(d, vals)` : set columns names
* `rename!(d, args)` : rename columns names based on keyword arguments
* `eltypes(d)` : `eltype` of each column
* `length(d)` : number of columns
* `size(d)` : (nrows, ncols)
* `head(d, n = 5)` : first `n` rows
* `tail(d, n = 5)` : last `n` rows
* `convert(Array, d)` : convert to an array
* `DataArray(d)` : convert to a DataArray
* `complete_cases(d)` : indexes of complete cases (rows with no NA's)
* `complete_cases!(d)` : remove rows with NA's
* `nonunique(d)` : indexes of duplicate rows
* `unique!(d)` : remove duplicate rows
* `similar(d)` : a DataFrame with similar columns as `d`

### Indexing

Table columns are accessed (`getindex`) by a single index that can be
a symbol identifier, an integer, or a vector of each. If a single
column is selected, just the column object is returned. If multiple
columns are selected, some AbstractDataFrame is returned.

```julia
d[:colA]
d[3]
d[[:colA, :colB]]
d[[1:3; 5]]
```

Rows and columns can be indexed like a `Matrix` with the added feature
of indexing columns by name.

```julia
d[1:3, :colA]
d[3,3]
d[3,:]
d[3,[:colA, :colB]]
d[:, [:colA, :colB]]
d[[1:3; 5], :]
```

`setindex` works similarly.
"""
abstract AbstractDataFrame

##############################################################################
##
## Interface (not final)
##
##############################################################################

# index(df) => AbstractIndex
# nrow(df) => Int
# ncol(df) => Int
# getindex(...)
# setindex!(...) exclusive of methods that add new columns

##############################################################################
##
## Basic properties of a DataFrame
##
##############################################################################

immutable Cols{T <: AbstractDataFrame}
    df::T
end
Base.start(::Cols) = 1
Base.done(itr::Cols, st) = st > length(itr.df)
Base.next(itr::Cols, st) = (itr.df[st], st + 1)
Base.length(itr::Cols) = length(itr.df)

# N.B. where stored as a vector, 'columns(x) = x.vector' is a bit cheaper
columns{T <: AbstractDataFrame}(df::T) = Cols{T}(df)

Base.names(df::AbstractDataFrame) = names(index(df))
_names(df::AbstractDataFrame) = _names(index(df))

"""
Set column names

```julia
names!(df::AbstractDataFrame, vals)
```

### Arguments

* `df` : the AbstractDataFrame
* `vals` : column names, normally a Vector{Symbol} the same length as
  the number of columns in `df`

### Result

* `::AbstractDataFrame` : the updated result

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
names!(df, [:a, :b, :c])
```

"""
function names!(df::AbstractDataFrame, vals)
    names!(index(df), vals)
    return df
end

function rename!(df::AbstractDataFrame, args...)
    rename!(index(df), args...)
    return df
end
rename!(f::Function, df::AbstractDataFrame) = rename!(df, f)

rename(df::AbstractDataFrame, args...) = rename!(copy(df), args...)
rename(f::Function, df::AbstractDataFrame) = rename(df, f)

"""
Rename columns

```julia
rename!(df::AbstractDataFrame, from::Symbol, to::Symbol)
rename!(df::AbstractDataFrame, d::Associative)
rename!(f::Function, df::AbstractDataFrame)
rename(df::AbstractDataFrame, from::Symbol, to::Symbol)
rename(f::Function, df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame
* `d` : an Associative type that maps the original name to a new name
* `f` : a function that has the old column name (a symbol) as input
  and new column name (a symbol) as output

### Result

* `::AbstractDataFrame` : the updated result

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
rename(x -> symbol(uppercase(string(x))), df)
rename(df, @compat(Dict(:i=>:A, :x=>:X)))
rename(df, :y, :Y)
rename!(df, @compat(Dict(:i=>:A, :x=>:X)))
```

"""
(rename!, rename)

"""
Column elemental types

```julia
eltypes(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame

### Result

* `::Vector{Type}` : the elemental type of each column

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
eltypes(df)
```

"""
function eltypes(df::AbstractDataFrame)
    ncols = size(df, 2)
    res = Array(Type, ncols)
    for j in 1:ncols
        res[j] = eltype(df[j])
    end
    return res
end

Base.size(df::AbstractDataFrame) = (nrow(df), ncol(df))
function Base.size(df::AbstractDataFrame, i::Integer)
    if i == 1
        nrow(df)
    elseif i == 2
        ncol(df)
    else
        throw(ArgumentError("DataFrames have only two dimensions"))
    end
end

Base.length(df::AbstractDataFrame) = ncol(df)
Base.endof(df::AbstractDataFrame) = ncol(df)

Base.ndims(::AbstractDataFrame) = 2

##############################################################################
##
## Similar
##
##############################################################################

Base.similar(df::AbstractDataFrame, dims::Int) =
    DataFrame([similar(x, dims) for x in columns(df)], copy(index(df)))

nas{T}(dv::AbstractArray{T}, dims::Union(Int, @compat Tuple{Vararg{Int}})) =   # TODO move to datavector.jl?
    DataArray(Array(T, dims), trues(dims))

nas{T,R}(dv::PooledDataArray{T,R}, dims::Union(Int, @compat Tuple{Vararg{Int}})) =
    PooledDataArray(DataArrays.RefArray(zeros(R, dims)), dv.pool)

nas(df::AbstractDataFrame, dims::Int) =
    DataFrame(Any[nas(x, dims) for x in columns(df)], copy(index(df)))

##############################################################################
##
## Equality
##
##############################################################################

function Base.isequal(df1::AbstractDataFrame, df2::AbstractDataFrame)
    size(df1, 2) == size(df2, 2) || return false
    isequal(index(df1), index(df2)) || return false
    for idx in 1:size(df1, 2)
        isequal(df1[idx], df2[idx]) || return false
    end
    return true
end

function Base.(:(==))(df1::AbstractDataFrame, df2::AbstractDataFrame)
    size(df1, 2) == size(df2, 2) || return false
    isequal(index(df1), index(df2)) || return false
    eq = true
    for idx in 1:size(df1, 2)
        coleq = df1[idx] == df2[idx]
        # coleq could be NA
        !isequal(coleq, false) || return false
        eq &= coleq
    end
    return eq
end

##############################################################################
##
## Associative methods
##
##############################################################################

Base.haskey(df::AbstractDataFrame, key::Any) = haskey(index(df), key)
Base.get(df::AbstractDataFrame, key::Any, default::Any) = haskey(df, key) ? df[key] : default
Base.isempty(df::AbstractDataFrame) = ncol(df) == 0

##############################################################################
##
## Description
##
##############################################################################

DataArrays.head(df::AbstractDataFrame, r::Int) = df[1:min(r,nrow(df)), :]
DataArrays.head(df::AbstractDataFrame) = head(df, 6)
DataArrays.tail(df::AbstractDataFrame, r::Int) = df[max(1,nrow(df)-r+1):nrow(df), :]
DataArrays.tail(df::AbstractDataFrame) = tail(df, 6)

"""
Show the first or last part of an AbstractDataFrame

```julia
head(df::AbstractDataFrame, r::Int = 6)
tail(df::AbstractDataFrame, r::Int = 6)
```

### Arguments

* `df` : the AbstractDataFrame
* `r` : the number of rows to show

### Result

* `::AbstractDataFrame` : the first or last part of `df`

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
head(df)
tail(df)
```

"""
(head, tail)

# get the structure of a DF
"""
Show the structure of an AbstractDataFrame, in a tree-like format

```julia
dump(df::AbstractDataFrame, n::Int = 5)
dump(io::IO, df::AbstractDataFrame, n::Int = 5)
```

### Arguments

* `df` : the AbstractDataFrame
* `n` : the number of levels to show
* `io` : optional output descriptor

### Result

* nothing

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
str(df)
```

"""
function Base.dump(io::IO, df::AbstractDataFrame, n::Int, indent)
    println(io, typeof(df), "  $(nrow(df)) observations of $(ncol(df)) variables")
    if n > 0
        for (name, col) in eachcol(df)
            print(io, indent, "  ", name, ": ")
            dump(io, col, n - 1, string(indent, "  "))
        end
    end
end

function Base.dump(io::IO, dv::AbstractDataVector, n::Int, indent)
    println(io, typeof(dv), "(", length(dv), ") ", dv[1:min(4, end)])
end

# summarize the columns of a DF
# TODO: clever layout in rows
"""
Summarize the columns of an AbstractDataFrame

```julia
describe(df::AbstractDataFrame)
describe(io, df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame
* `io` : optional output descriptor

### Result

* nothing

### Details

If the column's base type derives from Number, compute the minimum, first
quantile, median, mean, third quantile, and maximum. NA's are filtered and
reported separately.

For boolean columns, report trues, falses, and NAs.

For other types, show column characteristics and number of NAs.

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
describe(df)
```

"""
StatsBase.describe(df::AbstractDataFrame) = describe(STDOUT, df)
function StatsBase.describe(io, df::AbstractDataFrame)
    for (name, col) in eachcol(df)
        println(io, name)
        describe(io, col)
        println(io, )
    end
end
StatsBase.describe(dv::AbstractArray) = describe(STDOUT, dv)
function StatsBase.describe{T<:Number}(io, dv::AbstractArray{T})
    if all(isna(dv))
        println(io, " * All NA * ")
        return
    end
    filtered = float(dropna(dv))
    qs = quantile(filtered, [0, .25, .5, .75, 1])
    statNames = ["Min", "1st Qu.", "Median", "Mean", "3rd Qu.", "Max"]
    statVals = [qs[1:3]; mean(filtered); qs[4:5]]
    for i = 1:6
        println(io, string(rpad(statNames[i], 8, " "), " ", string(statVals[i])))
    end
    nas = sum(isna(dv))
    println(io, "NAs      $nas")
    println(io, "NA%      $(round(nas*100/length(dv), 2))%")
    return
end
function StatsBase.describe{T}(io, dv::AbstractArray{T})
    ispooled = isa(dv, PooledDataVector) ? "Pooled " : ""
    # if nothing else, just give the length and element type and NA count
    println(io, "Length  $(length(dv))")
    println(io, "Type    $(ispooled)$(string(eltype(dv)))")
    println(io, "NAs     $(sum(isna(dv)))")
    println(io, "NA%     $(round(sum(isna(dv))*100/length(dv), 2))%")
    println(io, "Unique  $(length(unique(dv)))")
    return
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

"""
Indexes of complete cases (rows without NA's)

```julia
complete_cases(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame

### Result

* `::Vector{Bool}` : indexes of complete cases

See also `complete_cases!`.

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
df[[1,4,5], :x] = NA
df[[9,10], :y] = NA
complete_cases(df)
```

"""
function complete_cases(df::AbstractDataFrame)
    ## Returns a Vector{Bool} of indexes of complete cases (rows with no NA's).
    res = !isna(df[1])
    for i in 2:ncol(df)
        res &= !isna(df[i])
    end
    res
end

"""
Delete rows with NA's.

```julia
complete_cases!(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame

### Result

* `::AbstractDataFrame` : the updated version

See also `complete_cases`.

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
df[[1,4,5], :x] = NA
df[[9,10], :y] = NA
complete_cases!(df)
```

"""
complete_cases!(df::AbstractDataFrame) = deleterows!(df, find(!complete_cases(df)))

function Base.convert(::Type{Array}, df::AbstractDataFrame)
    convert(Matrix, df)
end
function Base.convert(::Type{Matrix}, df::AbstractDataFrame)
    T = reduce(typejoin, eltypes(df))
    convert(Matrix{T}, df)
end
function Base.convert{T}(::Type{Array{T}}, df::AbstractDataFrame)
    convert(Matrix{T}, df)
end
function Base.convert{T}(::Type{Matrix{T}}, df::AbstractDataFrame)
    n, p = size(df)
    res = Array(T, n, p)
    idx = 1
    for col in columns(df)
        anyna(col) && error("DataFrame contains NAs")
        copy!(res, idx, data(col))
        idx += n
    end
    return res
end

function Base.convert(::Type{DataArray}, df::AbstractDataFrame)
    convert(DataMatrix, df)
end
function Base.convert(::Type{DataMatrix}, df::AbstractDataFrame)
    T = reduce(typejoin, eltypes(df))
    convert(DataMatrix{T}, df)
end
function Base.convert{T}(::Type{DataArray{T}}, df::AbstractDataFrame)
    convert(DataMatrix{T}, df)
end
function Base.convert{T}(::Type{DataMatrix{T}}, df::AbstractDataFrame)
    n, p = size(df)
    res = DataArray(T, n, p)
    idx = 1
    for col in columns(df)
        copy!(res, idx, col)
        idx += n
    end
    return res
end

"""
Indexes of complete cases (rows without NA's)

```julia
nonunique(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame

### Result

* `::Vector{Bool}` : indicates whether the row is a duplicate of some
  prior row

See also `unique` and `unique!`.

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
df = vcat(df, df)
nonunique(df)
```

"""
function nonunique(df::AbstractDataFrame)
    res = fill(false, nrow(df))
    rows = Set{DataFrameRow}()
    for i in 1:nrow(df)
        arow = DataFrameRow(df, i)
        if in(arow, rows)
            res[i] = true
        else
            push!(rows, arow)
        end
    end
    res
end

unique!(df::AbstractDataFrame) = deleterows!(df, find(nonunique(df)))

# Unique rows of an AbstractDataFrame.
Base.unique(df::AbstractDataFrame) = df[!nonunique(df), :]

"""
Delete duplicate rows

```julia
unique(df::AbstractDataFrame)
unique!(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame

### Result

* `::AbstractDataFrame` : the updated version

See also `nonunique` and `unique`.

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
df = vcat(df, df)
unique(df)   # doesn't modify df
unique!(df)  # modifies df
```

"""
(unique, unique!)

function nonuniquekey(df::AbstractDataFrame)
    # Here's another (probably a lot faster) way to do `nonunique`
    # by grouping on all columns. It will fail if columns cannot be
    # made into PooledDataVector's.
    gd = groupby(df, _names(df))
    idx = [1:length(gd.idx)][gd.idx][gd.starts]
    res = fill(true, nrow(df))
    res[idx] = false
    res
end

# Count the number of missing values in every column of an AbstractDataFrame.
function colmissing(df::AbstractDataFrame) # -> Vector{Int}
    nrows, ncols = size(df)
    missing = zeros(Int, ncols)
    for j in 1:ncols
        missing[j] = countna(df[j])
    end
    return missing
end

function without(df::AbstractDataFrame, icols::Vector{Int})
    newcols = _setdiff(1:ncol(df), icols)
    df[newcols]
end
without(df::AbstractDataFrame, i::Int) = without(df, [i])
without(df::AbstractDataFrame, c::Any) = without(df, index(df)[c])

##############################################################################
##
## Hcat / vcat
##
##############################################################################

# hcat's first argument must be an AbstractDataFrame
# Trailing arguments (currently) may also be DataVectors, Vectors, or scalars.

# hcat! is defined in dataframes/dataframes.jl
# Its first argument (currently) must be a DataFrame.

# catch-all to cover cases where indexing returns a DataFrame and copy doesn't
Base.hcat(df::AbstractDataFrame, x) = hcat!(df[:, :], x)

Base.hcat(df::AbstractDataFrame, x, y...) = hcat!(hcat(df, x), y...)

# vcat only accepts DataFrames. Finds union of columns, maintaining order
# of first df. Missing data becomes NAs.

Base.vcat(df::AbstractDataFrame) = df

Base.vcat(dfs::AbstractDataFrame...) = vcat(collect(dfs))

Base.vcat(dfs::Vector{None}) = dfs
function Base.vcat{T<:AbstractDataFrame}(dfs::Vector{T})
    coltyps, colnams, similars = _colinfo(dfs)

    res = DataFrame()
    Nrow = sum(nrow, dfs)
    for j in 1:length(colnams)
        colnam = colnams[j]
        col = similar(similars[j], coltyps[j], Nrow)

        i = 1
        for df in dfs
            if haskey(df, colnam)
                copy!(col, i, df[colnam])
            end
            i += size(df, 1)
        end

        res[colnam] = col
    end
    res
end

_isnullable(::AbstractArray) = false
_isnullable(::AbstractDataArray) = true
const EMPTY_DATA = DataArray(Void, 0)

function _colinfo{T<:AbstractDataFrame}(dfs::Vector{T})
    df1 = dfs[1]
    colindex = copy(index(df1))
    coltyps = eltypes(df1)
    similars = collect(columns(df1))
    nonnull_ct = Int[_isnullable(c) for c in columns(df1)]

    for i in 2:length(dfs)
        df = dfs[i]
        for j in 1:size(df, 2)
            col = df[j]
            cn, ct = _names(df)[j], eltype(col)
            if haskey(colindex, cn)
                idx = colindex[cn]

                oldtyp = coltyps[idx]
                if !(ct <: oldtyp)
                    coltyps[idx] = promote_type(oldtyp, ct)
                end
                nonnull_ct[idx] += !_isnullable(col)
            else # new column
                push!(colindex, cn)
                push!(coltyps, ct)
                push!(similars, col)
                push!(nonnull_ct, !_isnullable(col))
            end
        end
    end

    for j in 1:length(colindex)
        if nonnull_ct[j] < length(dfs) && !_isnullable(similars[j])
            similars[j] = EMPTY_DATA
        end
    end
    colnams = _names(colindex)

    coltyps, colnams, similars
end

##############################################################################
##
## Hashing
##
## Make sure this agrees with is_equals()
##
##############################################################################

function Base.hash(df::AbstractDataFrame)
    h = hash(size(df)) + 1
    for i in 1:size(df, 2)
        h = hash(df[i], h)
    end
    return @compat UInt(h)
end

## Documentation for methods defined elsewhere

"""
Number of rows or columns in an AbstractDataFrame

```julia
nrow(df::AbstractDataFrame)
ncol(df::AbstractDataFrame)
```

### Arguments

* `df` : the AbstractDataFrame

### Result

* `::AbstractDataFrame` : the updated version

See also `size`.

NOTE: these functions may be depreciated for `size`.

### Examples

```julia
df = DataFrame(i = 1:10, x = rand(10), y = rand(["a", "b", "c"], 10))
size(df)
nrow(df)
ncol(df)
```

""";
[:nrow, :ncol]


text:
include("buffer.jl")

import Random
import DataStructures
import MacroTools: @forward
import StatsBase

abstract type AbstractReplay end

proc_state(er::AbstractReplay, x) = identity(x)
Base.keys(er::AbstractReplay) = 1:length(er)

function Base.iterate(asr::AbstractReplay)

    state = 1
    result = asr[state]
    state += 1
    (result, state)
end

function Base.iterate(asr::AbstractReplay, state::Integer)

    if state > length(asr)
        return nothing
    end

    result = asr[state]
    state += 1

    (result, state)
end

mutable struct ExperienceReplay{CB<:CircularBuffer} <: AbstractReplay
    buffer::CB
end

const ExperienceReplayDef{TPL, TPS} = ExperienceReplay{CircularBuffer{TPL, TPS, (:s, :a, :sp, :r, :t)}}

ExperienceReplay(size, types, shapes, column_names) = begin
    cb = CircularBuffer(size, types, shapes, column_names)
    ExperienceReplay(cb)
end

ExperienceReplayDef(size, obs_size, obs_type) =
    ExperienceReplay(size,
                     (obs_type, Int, obs_type, Float32, Bool),
                     (obs_size, 1, obs_size, 1, 1),
                     (:s, :a, :sp, :r, :t))

Base.length(er::ExperienceReplay) = length(er.buffer)
Base.getindex(er::ExperienceReplay, idx) = er.buffer[idx]
Base.view(er::ExperienceReplay, idx) = @view er.buffer[idx]

Base.push!(er::ExperienceReplay, experience) = push!(er.buffer, experience)

StatsBase.sample(er::ExperienceReplay, batch_size) = StatsBase.sample(Random.GLOBAL_RNG, er, batch_size)

function StatsBase.sample(rng::Random.AbstractRNG, er::ExperienceReplay, batch_size)
    idx = rand(rng, 1:length(er), batch_size)
    return er[idx]
end

mutable struct DynaExperienceReplay{CB<:CircularBuffer} <: AbstractReplay
    buffer::CB
end

DynaExperienceReplay(size, types, shapes, column_names) = begin
    cb = CircularBuffer(size, types, shapes, column_names)
    DynaExperienceReplay(cb)
end

Base.length(er::DynaExperienceReplay) = length(er.buffer)
Base.getindex(er::DynaExperienceReplay, idx) = er.buffer[idx]
Base.view(er::DynaExperienceReplay, idx) = @view er.buffer[idx]

Base.push!(er::DynaExperienceReplay, experience) = push!(er.buffer, experience)

StatsBase.sample(er::DynaExperienceReplay, batch_size) = StatsBase.sample(Random.GLOBAL_RNG, er, batch_size)

function StatsBase.sample(rng::Random.AbstractRNG, er::DynaExperienceReplay, batch_size)
    idxes = rand(rng, 1:length(er), batch_size)
    [er[idx] for idx in idxes]
end


text:
using Random, Test
using LatexLH, ModelObjectsLH, ModelParams, CollegeEntry

ce = CollegeEntry;

function constructor_test()
    @testset "Constructors" begin
        J = 8; nc = 3;
        switches = make_entry_switches_oneloc(J, nc);
        @test validate_es(switches)
    end
end

function access_test(switches)
    @testset "Access routines" begin
        # println("\n--------------------------")
        # println(switches);
        objId = ObjectId(:entryDecision);
        st = ce.make_test_symbol_table();
        e = init_entry_decision(objId, switches, st);
        # println(e);

        J = n_types(e);
        nc = n_colleges(e);
        nl = n_locations(e);

        @test 0.0 < min_entry_prob(e) < max_entry_prob(e) < 1.0
        @test entry_pref_scale(e) > 0.0

        @test isa(CollegeEntry.value_local(e), Float64)

        capacityV = capacities(e);
        if any(capacityV .< 1e6)
            @test CollegeEntry.limited_capacity(e)
        else
            @test !CollegeEntry.limited_capacity(e)
        end

        typeMass_jlM = type_mass_jl(e);
        @test size(typeMass_jlM) == (J, nl)
        @test all(typeMass_jlM .>= 0.0)
        for j = 1 : J
            @test isapprox(typeMass_jlM[j,:], type_mass_jl(e, j))
            @test isapprox(sum(typeMass_jlM[j,:]),  type_mass_j(e, j))
        end
        @test type_mass_jl(e, J, nl) == typeMass_jlM[J, nl]

        capacity_clM = capacities(e);
        @test size(capacity_clM) == (nc, nl)
        for ic = 1 : nc
            @test isapprox(capacity_clM[ic,:], CollegeEntry.capacity(e, ic))
        end

        if nl == 1
            @test CollegeEntry.value_local(e) == 0.0
        end
    end
end

function subset_switches_test(switches)
    @testset "Subset switches" begin
        idxV = 2 : 2 : n_types(switches);
        CollegeEntry.subset_types!(switches, idxV);
        @test validate_es(switches)
        @test n_types(switches) == length(idxV)
	end
end

# Test `entry_probs` which has no notion of locations
function entry_test(switches, prefShocks :: Bool)
    rng = MersenneTwister(49);
    @testset "Entry probs: $switches" begin
        F1 = Float64;
        st = ce.make_test_symbol_table();
        e = init_entry_decision(ObjectId(:entry), switches, st);
        # println("\n------------")
        # println(e);
        # println("Preference shocks: $prefShocks");

        J = n_types(e);
        nc = n_colleges(e);
        nl = n_locations(e);
        vWork_jV, vCollege_jcM = CollegeEntry.values_for_test(rng, J, nc, nl);
        admitV = [1, 3];
        rejectV = [2];
        prob_jcM, eVal_jV = entry_probs(e, vWork_jV, vCollege_jcM, admitV;
            prefShocks = prefShocks);
        @test all(prob_jcM .>= 0.0)
        @test all(prob_jcM .<= 1.0)
        @test size(prob_jcM) == (J, nc)
        @test size(eVal_jV) == (J,)
        @test all(prob_jcM[:, rejectV] .== 0.0)

        # One person at a time
        for j = 1 : J
            prob_cV, eVal =
                entry_probs(e, vWork_jV[j], vCollege_jcM[j,:], admitV;
                    prefShocks = prefShocks);
            @test isapprox(prob_cV, prob_jcM[j,:])
            @test isapprox(eVal, eVal_jV[j])
        end

        # Increasing a value should increase probability
        # Not for limited capacities, though
        # Also not without pref shocks
        if !limited_capacity(e)  &&  prefShocks
            idx = admitV[end];
            otherAdmitV = admitV[1 : (end-1)];
            vCollege_jcM[:, idx] .+= 0.1;
            prob2_jcM, eVal2_jV = entry_probs(e, vWork_jV, vCollege_jcM, admitV;
                prefShocks = prefShocks);
            @test all(prob2_jcM[:, idx] .> prob_jcM[:, idx])
            @test all(prob2_jcM[:, otherAdmitV] .< prob_jcM[:, otherAdmitV])
            @test all(eVal2_jV .> eVal_jV)
        end

        # Empty admission set
        prob_jcM, eVal_jV = entry_probs(e, vWork_jV, vCollege_jcM, [];
            prefShocks = prefShocks);
        @test size(prob_jcM) == (J, nc)
        @test all(prob_jcM .== 0.0)
        # if !isa(e, EntryTwoStep)
            @test isapprox(eVal_jV, vWork_jV)
        # end
    end
end

# Check that small preference shocks give about the same answer as no preference shocks
function small_pref_entry_test(switches)
    rng = MersenneTwister(49);
    @testset "Entry probs" begin
        F1 = Float64;
        ce.set_pref_scale!(switches, 0.001);
        st = ce.make_test_symbol_table();
        e = init_entry_decision(ObjectId(:entry), switches, st);
        # println("\n------------")
        # println(e);

        J = n_types(e);
        nc = n_colleges(e);
        nl = n_locations(e);
        vWork_jV, vCollege_jcM = CollegeEntry.values_for_test(rng, J, nc, nl);
        admitV = [1, 3];
        rejectV = [2];

        prob_jcM, eVal_jV = entry_probs(e, vWork_jV, vCollege_jcM, admitV;
            prefShocks = true);
        @test all(prob_jcM .>= 0.0)
        @test all(prob_jcM .<= 1.0)
        @test size(prob_jcM) == (J, nc)
        @test size(eVal_jV) == (J,)
        @test all(prob_jcM[:, rejectV] .== 0.0)

        prob2_jcM, eVal2_jV = entry_probs(e, vWork_jV, vCollege_jcM, admitV;
            prefShocks = false);
        @test isapprox(prob_jcM, prob2_jcM, atol = 0.01)
        @test isapprox(eVal_jV, eVal2_jV, rtol = 0.01)
    end
end

function sim_entry_test(switches)
    @testset "Simulate Entry probs" begin
        rng = MersenneTwister(123);
        F1 = Float64;
        st = ce.make_test_symbol_table();
        e = init_entry_decision(ObjectId(:entry), switches, st);
        # println("\n------------")
        # println(e);

        J = n_types(e);
        nc = n_colleges(e);
        nl = n_locations(e);
        vCollege_clM = test_value_cl(rng, nc, nl);
        vWork = 1.0 .+ sum(vCollege_clM) / length(vCollege_clM);
        avail_clM = rand(rng, Bool, nc, nl);
        # Make sure one college is available
        avail_clM[1,1] = true;

        prob_clV, eVal = entry_probs(e, vWork, vec(vCollege_clM), vec(avail_clM));
        prob_clM = reshape(prob_clV, nc, nl);
        @test !any(prob_clM[.!avail_clM] .> 0.0)

        # Tests expected value by simulation
        nSim = Int(1e5);
        prob2_clM, eVal2 = CollegeEntry.sim_entry_probs(e,
            vWork, vCollege_clM, avail_clM,
            nSim, rng);
        @test !any(prob2_clM[.!avail_clM] .> 0.0)

        @test isapprox(eVal, eVal2, atol = 1e-4)
        # @show eVal, eVal2;
        @test isapprox(prob_clM, prob2_clM, atol = 0.02)
    end
end

function scale_entry_probs_test()
    rng = MersenneTwister(123);
    @testset "Scale entry probs" begin
        J = 30; nl = 4; nc = 5;
        entryProb_jlcM = make_test_entry_probs(rng, J, nc, nl);
        minEntryProb = 0.03;
        maxEntryProb = 0.8;
        scale_entry_probs!(entryProb_jlcM, minEntryProb, maxEntryProb);
        entryProb_jlM = sum(entryProb_jlcM, dims = 3);
        @test all(entryProb_jlM .< maxEntryProb + 0.01)
        @test all(entryProb_jlcM .> 0.0)

        # With probs in range, nothing should change
        entryProb_jlcM = fill(minEntryProb + 0.01, J, nl, nc);
        entryProb2_jlcM = copy(entryProb_jlcM);
        scale_entry_probs!(entryProb_jlcM, minEntryProb, maxEntryProb);
        @test isapprox(entryProb_jlcM, entryProb2_jlcM)
	end
end

@testset "All" begin
    constructor_test()
    scale_entry_probs_test();

    J = 8; nc = 3;
    for switches in test_entry_switches(J, nc)
        access_test(switches);
        subset_switches_test(switches);
        for prefShocks ∈ (true, false)
            entry_test(switches, prefShocks);
            small_pref_entry_test(switches);
        end
        sim_entry_test(switches);
    end
end

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


text:
abstract type SDEIntegratorCache{DT,D,M} <: IntegratorCache{DT,D} end
abstract type PSDEIntegratorCache{DT,D,M} <: IntegratorCache{DT,D} end