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

using Base.Filesystem: mktempdir

include("FMI2_c.jl")
include("FMI2_comp.jl")
include("FMI2_md.jl")

"""
ToDo
"""
@enum fmi2Dependency begin
    fmi2DependencyUnknown
    fmi2DependencyDependent
    fmi2DependencyIndependent
    fmi2DependencyFixed
end

"""
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 <: FMU
    modelName::fmi2String
    instanceName::fmi2String
    fmuResourceLocation::fmi2String

    modelDescription::fmi2ModelDescription

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

    # paths of ziped and unziped FMU folders
    path::String
    zipPath::String

    # 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}

    # START: experimental section (to FMIFlux.jl)
    dependencies::Matrix{fmi2Dependency}

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

    x       # current state
    next_x  # next state

    dx      # current state derivative
    simulationResult

    jac_dxy_x::Matrix{fmi2Real}
    jac_dxy_u::Matrix{fmi2Real}
    jac_x::Array{fmi2Real}
    jac_t::fmi2Real
    # END: experimental section

    # Constructor
    FMU2() = new()
end

# wrapper functions on the model description
function fmi2GetModelName(fmu::FMU2)
    fmi2GetModelName(fmu.modelDescription)
end
function fmi2GetGUID(fmu::FMU2)
    fmi2GetGUID(fmu.modelDescription)
end
function fmi2GetGenerationTool(fmu::FMU2)
    fmi2GetGenerationTool(fmu.modelDescription)
end
function fmi2GetGenerationDateAndTime(fmu::FMU2)
    fmi2GetGenerationDateAndTime(fmu.modelDescription)
end
function fmi2GetVariableNamingConvention(fmu::FMU2)
    fmi2GetVariableNamingConvention(fmu.modelDescription)
end
function fmi2GetNumberOfEventIndicators(fmu::FMU2)
    fmi2GetNumberOfEventIndicators(fmu.modelDescription)
end

function fmi2CanGetSetState(fmu::FMU2)
    fmi2CanGetSetState(fmu.modelDescription)
end
function fmi2CanSerializeFMUstate(fmu::FMU2)
    fmi2CanSerializeFMUstate(fmu.modelDescription)
end
function fmi2ProvidesDirectionalDerivative(fmu::FMU2)
    fmi2ProvidesDirectionalDerivative(fmu.modelDescription)
end
function fmi2IsCoSimulation(fmu::FMU2)
    fmi2IsCoSimulation(fmu.modelDescription)
end
function fmi2IsModelExchange(fmu::FMU2)
    fmi2IsModelExchange(fmu.modelDescription)
end

# TODO: add to documentation and export list, if not already done

"""
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 reference 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 an array of ValueReferences coresponding to the variable names.
"""
function fmi2String2ValueReference(md::fmi2ModelDescription, names::Array{String})
    vr = Array{fmi2ValueReference}(undef,0)
    for name in names
        reference = fmi2String2ValueReference(md, name)
        if reference == nothing
            @warn "Value reference for variable '$name' not found, skipping."
        else
            push!(vr, reference)
        end
    end
    vr
end

"""
Returns the ValueReference coresponding to the variable name.
"""
function fmi2String2ValueReference(md::fmi2ModelDescription, name::String)
    reference = nothing
    if haskey(md.stringValueReferences, name)
        reference = md.stringValueReferences[name]
    else
        @warn "No variable named '$name' found."
    end
    reference
end

function fmi2String2ValueReference(fmu::FMU2, name::Union{String, Array{String}})
    fmi2String2ValueReference(fmu.modelDescription, name)
end

"""
Returns an array of variable names matching a fmi2ValueReference.
"""
function fmi2ValueReference2String(md::fmi2ModelDescription, reference::fmi2ValueReference)
    [k for (k,v) in md.stringValueReferences if v == reference]
end
function fmi2ValueReference2String(md::fmi2ModelDescription, reference::Int64)
    fmi2ValueReference2String(md, fmi2ValueReference(reference))
end

function fmi2ValueReference2String(fmu::FMU2, reference::Union{fmi2ValueReference, Int64})
    fmi2ValueReference2String(fmu.modelDescription, reference)
end

"""
Create a copy of the .fmu file as a .zip folder and unzips it.
Returns the paths to the zipped and unzipped folders.

Via optional argument ```unpackPath```, a path to unpack the FMU can be specified (default: system temporary directory).
"""
function fmi2Unzip(pathToFMU::String; unpackPath=nothing)

    fileNameExt = basename(pathToFMU)
    (fileName, fileExt) = splitext(fileNameExt)

    if unpackPath == nothing
        # cleanup=true leads to issues with automatic testing on linux server.
        unpackPath = mktempdir(; prefix="fmijl_", cleanup=false)
    end

    zipPath = joinpath(unpackPath, fileName * ".zip")
    unzippedPath = joinpath(unpackPath, fileName)

    # only copy ZIP if not already there
    if !isfile(zipPath)
        cp(pathToFMU, zipPath; force=true)
    end

    @assert isfile(zipPath) ["fmi2Unzip(...): ZIP-Archive couldn't be copied to `$zipPath`."]

    zipAbsPath = isabspath(zipPath) ?  zipPath : joinpath(pwd(), zipPath)
    unzippedAbsPath = isabspath(unzippedPath) ? unzippedPath : joinpath(pwd(), unzippedPath)

    @assert isfile(zipAbsPath) ["fmi2Unzip(...): Can't deploy ZIP-Archive at `$(zipAbsPath)`."]

    numFiles = 0

    # only unzip if not already done
    if !isdir(unzippedAbsPath)
        mkpath(unzippedAbsPath)

        zarchive = ZipFile.Reader(zipAbsPath)
        for f in zarchive.files
            fileAbsPath = normpath(joinpath(unzippedAbsPath, f.name))

            if endswith(f.name,"/") || endswith(f.name,"\\")
                mkpath(fileAbsPath) # mkdir(fileAbsPath)

                @assert isdir(fileAbsPath) ["fmi2Unzip(...): Can't create directory `$(f.name)` at `$(fileAbsPath)`."]
            else
                # create directory if not forced by zip file folder
                mkpath(dirname(fileAbsPath))

                numBytes = write(fileAbsPath, read(f))

                if numBytes == 0
                    @info "fmi2Unzip(...): Written file `$(f.name)`, but file is empty."
                end

                @assert isfile(fileAbsPath) ["fmi2Unzip(...): Can't unzip file `$(f.name)` at `$(fileAbsPath)`."]
                numFiles += 1
            end
        end
        close(zarchive)
    end

    @assert isdir(unzippedAbsPath) ["fmi2Unzip(...): ZIP-Archive couldn't be unzipped at `$(unzippedPath)`."]
    @info "fmi2Unzip(...): Successfully unzipped $numFiles files at `$unzippedAbsPath`."

    (unzippedAbsPath, zipAbsPath)
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 == nothing
        @info "This FMU does not support optional function '$symbol'."
        addr = Ptr{Cvoid}(C_NULL)
    end
    addr
end

"""
Sets the properties of the fmu by reading the modelDescription.xml.
Retrieves all the pointers of binary functions.

Returns the instance of the FMU struct.

Via optional argument ```unpackPath```, a path to unpack the FMU can be specified (default: system temporary directory).
"""
function fmi2Load(pathToFMU::String; unpackPath=nothing)
    # Create uninitialized FMU
    fmu = FMU2()
    fmu.components = []

    pathToFMU = normpath(pathToFMU)

    # set paths for fmu handling
    (fmu.path, fmu.zipPath) = fmi2Unzip(pathToFMU; unpackPath=unpackPath)

    # set paths for modelExchangeScripting and binary
    tmpName = splitpath(fmu.path)
    pathToModelDescription = joinpath(fmu.path, "modelDescription.xml")

    # parse modelDescription.xml
    fmu.modelDescription = fmi2ReadModelDescription(pathToModelDescription)
    fmu.modelName = fmu.modelDescription.modelName
    fmu.instanceName = fmu.modelDescription.modelName
    fmuName = fmi2GetModelIdentifier(fmu.modelDescription) # tmpName[length(tmpName)]

    directoryBinary = ""
    pathToBinary = ""

    if Sys.iswindows()
        directories = [joinpath("binaries", "win64"), joinpath("binaries","x86_64-windows")]
        for directory in directories
            directoryBinary = joinpath(fmu.path, directory)
            if isdir(directoryBinary)
                pathToBinary = joinpath(directoryBinary, "$(fmuName).dll")
                break
            end
        end
        @assert isfile(pathToBinary) "Target platform is Windows, but can't find valid FMU binary at `$(pathToBinary)` for path `$(fmu.path)`."
    elseif Sys.islinux()
        directories = [joinpath("binaries", "linux64"), joinpath("binaries", "x86_64-linux")]
        for directory in directories
            directoryBinary = joinpath(fmu.path, 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 at `$(pathToBinary)` for path `$(fmu.path)`."
    elseif Sys.isapple()
        directories = [joinpath("binaries", "darwin64"), joinpath("binaries", "x86_64-darwin")]
        for directory in directories
            directoryBinary = joinpath(fmu.path, 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 at `$(pathToBinary)` for path `$(fmu.path)`."
    else
        @assert false "Unsupported target platform. Supporting Windows64, Linux64 and Mac64."
    end

    lastDirectory = pwd()
    cd(directoryBinary)

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

    cd(lastDirectory)

    if fmi2IsCoSimulation(fmu)
        fmu.type = fmi2CoSimulation::fmi2Type
    elseif fmi2IsModelExchange(fmu)
        fmu.type = fmi2ModelExchange::fmi2Type
    else
        error(unknownFMUType)
    end

    if fmi2IsCoSimulation(fmu) && fmi2IsModelExchange(fmu)
        @info "fmi2Load(...): FMU supports both CS and ME, using CS as default if nothing specified."
    end

    # make URI ressource location
    tmpResourceLocation = string("file:///", fmu.path)
    tmpResourceLocation = joinpath(tmpResourceLocation, "resources")
    fmu.fmuResourceLocation = replace(tmpResourceLocation, "\\" => "/") # URIs.escapeuri(tmpResourceLocation)

    @info "fmi2Load(...): FMU resources location is `$(fmu.fmuResourceLocation)`"

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

    fmu.cGetString                    = dlsym_opt(fmu.libHandle, :fmi2GetString)
    fmu.cSetString                    = dlsym_opt(fmu.libHandle, :fmi2SetString)

    if fmi2CanGetSetState(fmu)
        fmu.cGetFMUstate                  = dlsym_opt(fmu.libHandle, :fmi2GetFMUstate)
        fmu.cSetFMUstate                  = dlsym_opt(fmu.libHandle, :fmi2SetFMUstate)
        fmu.cFreeFMUstate                 = dlsym_opt(fmu.libHandle, :fmi2FreeFMUstate)
    end

    if fmi2CanSerializeFMUstate(fmu)
        fmu.cSerializedFMUstateSize       = dlsym_opt(fmu.libHandle, :fmi2SerializedFMUstateSize)
        fmu.cSerializeFMUstate            = dlsym_opt(fmu.libHandle, :fmi2SerializeFMUstate)
        fmu.cDeSerializeFMUstate          = dlsym_opt(fmu.libHandle, :fmi2DeSerializeFMUstate)
    end

    if fmi2ProvidesDirectionalDerivative(fmu)
        fmu.cGetDirectionalDerivative     = dlsym_opt(fmu.libHandle, :fmi2GetDirectionalDerivative)
    end

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

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

    # initialize further variables
    fmu.jac_x = zeros(Float64, fmu.modelDescription.numberOfContinuousStates)
    fmu.jac_t = -1.0
    fmu.jac_dxy_x = zeros(fmi2Real,0,0)
    fmu.jac_dxy_u = zeros(fmi2Real,0,0)

    # dependency matrix
    # fmu.dependencies

    fmu
end

"""
Returns how a variable depends on another variable based on the model description.
"""
function fmi2VariableDependsOnVariable(fmu::FMU2, vr1::fmi2ValueReference, vr2::fmi2ValueReference)
    i1 = fmu.modelDescription.valueReferenceIndicies[vr1]
    i2 = fmu.modelDescription.valueReferenceIndicies[vr2]
    return fmi2GetDependencies(fmu)[i1, i2]
end

"""
Returns the FMU's dependency-matrix for fast look-ups on dependencies between value references.

Entries are from type fmi2Dependency, possible values are "Dependent", "Independent" and "Unknown".
"""
function fmi2GetDependencies(fmu::FMU2)
    if !isdefined(fmu, :dependencies)
        dim = length(fmu.modelDescription.valueReferences)
        @info "fmi2GetDependencies: Started building dependency matrix $(dim) x $(dim) ..."

        if fmi2DependenciesSupported(fmu.modelDescription)
            fmu.dependencies = fill(fmi2DependencyIndependent::fmi2Dependency, dim, dim)

            for i in 1:dim
                modelVariable = fmi2ModelVariablesForValueReference(fmu.modelDescription, fmu.modelDescription.valueReferences[i])[1]

                if modelVariable.dependencies != nothing
                    indicies = collect(fmu.modelDescription.valueReferenceIndicies[fmu.modelDescription.modelVariables[dependency].valueReference] for dependency in modelVariable.dependencies)
                    dependenciesKind = modelVariable.dependenciesKind

                    k = 1
                    for j in 1:dim
                        if j in indicies
                            if dependenciesKind[k] == "fixed"
                                fmu.dependencies[i,j] = fmi2DependencyFixed::fmi2Dependency
                            elseif dependenciesKind[k] == "dependent"
                                fmu.dependencies[i,j] = fmi2DependencyDependent::fmi2Dependency
                            else
                                @warn "Unknown dependency kind for index ($i, $j) = `$(dependenciesKind[k])`."
                            end
                            k += 1
                        end
                    end
                end
            end
        else
            fmu.dependencies = fill(fmi2DependencyUnknown::fmi2Dependency, dim, dim)
        end

        @info "fmi2GetDependencies: Building dependency matrix $(dim) x $(dim) finished."
    end

    fmu.dependencies
end

function fmi2PrintDependencies(fmu::FMU2)
    dep = fmi2GetDependencies(fmu)
    ni, nj = size(dep)

    for i in 1:ni
        str = ""
        for j in 1:nj
            str = "$(str) $(Integer(dep[i,j]))"
        end
        println(str)
    end
end

"""
Prints FMU related information.
"""
function fmi2Info(fmu::FMU2)
    println("#################### Begin information for FMU ####################")

    println("\tModel name:\t\t\t$(fmi2GetModelName(fmu))")
    println("\tFMI-Version:\t\t\t$(fmi2GetVersion(fmu))")
    println("\tGUID:\t\t\t\t$(fmi2GetGUID(fmu))")
    println("\tGeneration tool:\t\t$(fmi2GetGenerationTool(fmu))")
    println("\tGeneration time:\t\t$(fmi2GetGenerationDateAndTime(fmu))")
    println("\tVar. naming conv.:\t\t$(fmi2GetVariableNamingConvention(fmu))")
    println("\tEvent indicators:\t\t$(fmi2GetNumberOfEventIndicators(fmu))")

    println("\tInputs:\t\t\t\t$(length(fmu.modelDescription.inputValueReferences))")
    for vr in fmu.modelDescription.inputValueReferences
        println("\t\t$(vr) $(fmi2ValueReference2String(fmu, vr))")
    end

    println("\tOutputs:\t\t\t$(length(fmu.modelDescription.outputValueReferences))")
    for vr in fmu.modelDescription.outputValueReferences
        println("\t\t$(vr) $(fmi2ValueReference2String(fmu, vr))")
    end

    println("\tStates:\t\t\t\t$(length(fmu.modelDescription.stateValueReferences))")
    for vr in fmu.modelDescription.stateValueReferences
        println("\t\t$(vr) $(fmi2ValueReference2String(fmu, vr))")
    end

    println("\tSupports Co-Simulation:\t\t$(fmi2IsCoSimulation(fmu))")
    if fmi2IsCoSimulation(fmu)
        println("\t\tModel identifier:\t$(fmu.modelDescription.CSmodelIdentifier)")
        println("\t\tGet/Set State:\t\t$(fmu.modelDescription.CScanGetAndSetFMUstate)")
        println("\t\tSerialize State:\t$(fmu.modelDescription.CScanSerializeFMUstate)")
        println("\t\tDir. Derivatives:\t$(fmu.modelDescription.CSprovidesDirectionalDerivative)")

        println("\t\tVar. com. steps:\t$(fmu.modelDescription.CScanHandleVariableCommunicationStepSize)")
        println("\t\tInput interpol.:\t$(fmu.modelDescription.CScanInterpolateInputs)")
        println("\t\tMax order out. der.:\t$(fmu.modelDescription.CSmaxOutputDerivativeOrder)")
    end

    println("\tSupports Model-Exchange:\t$(fmi2IsModelExchange(fmu))")
    if fmi2IsModelExchange(fmu)
        println("\t\tModel identifier:\t$(fmu.modelDescription.MEmodelIdentifier)")
        println("\t\tGet/Set State:\t\t$(fmu.modelDescription.MEcanGetAndSetFMUstate)")
        println("\t\tSerialize State:\t$(fmu.modelDescription.MEcanSerializeFMUstate)")
        println("\t\tDir. Derivatives:\t$(fmu.modelDescription.MEprovidesDirectionalDerivative)")
    end

    println("##################### End information for FMU #####################")
end

"""
Starts a simulation of the fmu instance for the matching fmu type. If both types are available, CS is preferred over ME.
"""
function fmi2Simulate(fmu::FMU2, t_start::Real = 0.0, t_stop::Real = 1.0;
                      recordValues::fmi2ValueReferenceFormat = nothing, saveat=[], setup=true)
    fmi2Simulate(fmu.components[end], t_start, t_stop;
                 recordValues=recordValues, saveat=saveat, setup=setup)
end

"""
Starts a simulation of a FMU in CS-mode.
"""
function fmi2SimulateCS(fmu::FMU2, t_start::Real, t_stop::Real;
                        recordValues::fmi2ValueReferenceFormat = nothing, saveat=[], setup=true)
    fmi2SimulateCS(fmu.components[end], t_start, t_stop;
                   recordValues=recordValues, saveat=saveat, setup=setup)
end

"""
Starts a simulation of a FMU in ME-mode.
"""
function fmi2SimulateME(fmu::FMU2, t_start::Real, t_stop::Real;
                        recordValues::fmi2ValueReferenceFormat = nothing, saveat = [], setup = true, solver = nothing)
    fmi2SimulateME(fmu.components[end], t_start, t_stop;
                   recordValues=recordValues, saveat=saveat, setup=setup, solver=solver)
end

"""
Unload a FMU.

Free the allocated memory, close the binaries and remove temporary zip and unziped FMU model description.
"""
function fmi2Unload(fmu::FMU2, cleanUp::Bool = true)

    while length(fmu.components) > 0
        fmi2FreeInstance!(fmu.components[end])
    end

    dlclose(fmu.libHandle)

    # the components are removed from the component list via call to fmi2FreeInstance!
    @assert length(fmu.components) == 0 "fmi2Unload(...): Failure during deleting components, $(length(fmu.components)) remaining in stack."

    if cleanUp
        try
            rm(fmu.path; recursive = true, force = true)
            rm(fmu.zipPath; recursive = true, force = true)
        catch e
            @warn "Cannot delete unpacked data on disc. Maybe some files are opened in another application."
        end
    end
end

"""
TODO: FMI specification reference.

Returns a string representing the header file used to compile the FMU.

Returns "default" by default.

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

"""
Returns the number of states of the FMU.
"""
function fmi2GetNumberOfStates(fmu::FMU2)
    length(fmu.modelDescription.stateValueReferences)
end

"""
TODO: FMI specification reference.

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

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

"""
TODO: FMI specification reference.

Create a new instance of the given fmu, adds a logger if logginOn == true.

Returns the instance of a new FMU component.

For more information call ?fmi2Instantiate
"""
function fmi2Instantiate!(fmu::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
    fmu.callbackFunctions = fmi2CallbackFunctions(ptrLogger, ptrAllocateMemory, ptrFreeMemory, ptrStepFinished, C_NULL)

    compAddr = fmi2Instantiate(fmu.cInstantiate, fmu.instanceName, fmu.type, fmu.modelDescription.guid, fmu.fmuResourceLocation, fmu.callbackFunctions, fmi2Boolean(visible), fmi2Boolean(loggingOn))

    if compAddr == Ptr{Cvoid}(C_NULL)
        @error "fmi2Instantiate!(...): Instantiation failed!"
        return nothing
    end

    component = fmi2Component(compAddr, fmu)
    push!(fmu.components, component)
    component
end

"""
TODO: FMI specification reference.

Free the allocated memory used for the logger and fmu2 instance and destroy the instance.

For more information call ?fmi2FreeInstance
"""
function fmi2FreeInstance!(fmu::FMU2)
    fmi2FreeInstance!(fmu.components[end])
    pop!(fmu.components)
    nothing
end

"""
TODO: FMI specification reference.

Sets debug logging for the FMU.

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

"""
TODO: FMI specification reference.

Setup the simulation.

For more information call ?fmi2SetupExperiment
"""
function fmi2SetupExperiment(fmu::FMU2, startTime::Real = 0.0, stopTime::Real = startTime; tolerance::Real = 0.0)
    fmi2SetupExperiment(fmu.components[end], startTime, stopTime; tolerance=tolerance)
end

"""
TODO: FMI specification reference.

FMU enters Initialization mode.

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

"""
TODO: FMI specification reference.

FMU exits Initialization mode.

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

"""
TODO: FMI specification reference.

Informs FMU that simulation run is terminated.

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

"""
TODO: FMI specification reference.

Resets FMU.

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

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Real variables.

For more information call ?fmi2GetReal!
"""
function fmi2GetReal(fmu::FMU2, vr::fmi2ValueReferenceFormat)
    fmi2GetReal(fmu.components[end], vr)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Real variables.

For more information call ?fmi2GetReal!
"""
function fmi2GetReal!(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{<:Real}, <:Real})
    fmi2GetReal!(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2Real variables.

For more information call ?fmi2SetReal
"""
function fmi2SetReal(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{<:Real}, <:Real})
    fmi2SetReal(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Integer variables.

For more information call ?fmi2GetInteger!
"""
function fmi2GetInteger(fmu::FMU2, vr::fmi2ValueReferenceFormat)
    fmi2GetInteger(fmu.components[end], vr)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Integer variables.

For more information call ?fmi2GetInteger!
"""
function fmi2GetInteger!(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{<:Integer}, <:Integer})
    fmi2GetInteger!(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2Integer variables.

For more information call ?fmi2SetInteger
"""
function fmi2SetInteger(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{<:Integer}, <:Integer})
    fmi2SetInteger(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Boolean variables.

For more information call ?fmi2GetBoolean!
"""
function fmi2GetBoolean(fmu::FMU2, vr::fmi2ValueReferenceFormat)
    fmi2GetBoolean(fmu.components[end], vr)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2Boolean variables.

For more information call ?fmi2GetBoolean!
"""
function fmi2GetBoolean!(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{Bool}, Bool, Array{fmi2Boolean}})
    fmi2GetBoolean!(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2Boolean variables.

For more information call ?fmi2SetBoolean
"""
function fmi2SetBoolean(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{Bool}, Bool})
    fmi2SetBoolean(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2String variables.

For more information call ?fmi2GetString!
"""
function fmi2GetString(fmu::FMU2, vr::fmi2ValueReferenceFormat)
    fmi2GetString(fmu.components[end], vr)
end

"""
TODO: FMI specification reference.

Get the values of an array of fmi2String variables.

For more information call ?fmi2GetString!
"""
function fmi2GetString!(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{String}, String})
    fmi2GetString!(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Set the values of an array of fmi2String variables.

For more information call ?fmi2SetString
"""
function fmi2SetString(fmu::FMU2, vr::fmi2ValueReferenceFormat, values::Union{Array{String}, String})
    fmi2SetString(fmu.components[end], vr, values)
end

"""
TODO: FMI specification reference.

Get the pointer to the current FMU state.

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

"""
TODO: FMI specification reference.

Set the FMU to the given fmi2FMUstate.

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

"""
TODO: FMI specification reference.

Free the allocated memory for the FMU state.

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

"""
TODO: FMI specification reference.

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

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

"""
TODO: FMI specification reference.

Serialize the data in the FMU state pointer.

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

"""
TODO: FMI specification reference.

Deserialize the data in the serializedState fmi2Byte field.

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

"""
TODO: FMI specification reference.

Retrieves directional derivatives.

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

    fmi2GetDirectionalDerivative(fmu.components[end], vUnknown_ref, vKnown_ref, dvKnown)
end

"""
TODO: FMI specification reference.

Retrieves directional derivatives.

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

    fmi2GetDirectionalDerivative(fmu.components[end], vUnknown_ref, vKnown_ref, dvKnown)
end

"""
TODO: FMI specification reference.

Retrieves directional derivatives in-place.

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

    fmi2GetDirectionalDerivative!(fmu.components[end], vUnknown_ref, vKnown_ref, dvUnknown, dvKnown)
end

"""
This function approximates the directional derivative by sampling corresponding values (1st order, central differences).
"""
function fmi2SampleDirectionalDerivative(fmu::FMU2,
                                       vUnknown_ref::Array{fmi2ValueReference},
                                       vKnown_ref::Array{fmi2ValueReference},
                                       steps::Array{fmi2Real} = ones(fmi2Real, length(vKnown_ref)).*1e-5)
    fmi2SampleDirectionalDerivative(fmu.components[end], vUnknown_ref, vKnown_ref, steps)
end
function fmi2SampleDirectionalDerivative!(fmu::FMU2,
                                          vUnknown_ref::Array{fmi2ValueReference},
                                          vKnown_ref::Array{fmi2ValueReference},
                                          dvUnknown::AbstractArray,
                                          steps::Array{fmi2Real} = ones(fmi2Real, length(vKnown_ref)).*1e-5)
    fmi2SampleDirectionalDerivative!(fmu.components[end], vUnknown_ref, vKnown_ref, dvUnknown, steps)
end
function fmi2SampleDirectionalDerivative(fmu::FMU2,
                                       vUnknown_ref::fmi2ValueReference,
                                       vKnown_ref::fmi2ValueReference,
                                       steps::fmi2Real = 1e-5)
    fmi2SampleDirectionalDerivative(fmu.components[end], [vUnknown_ref], [vKnown_ref], [steps])[1]
end

"""
TODO: FMI specification reference.

Sets the n-th time derivative of real input variables.

vr defines the value references of the variables
the array order specifies the corresponding order of derivation of the variables

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

"""
TODO: FMI specification reference.

Retrieves the n-th derivative of output values.

vr defines the value references of the variables
the array order specifies the corresponding order of derivation of the variables

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

"""
TODO: FMI specification reference.

The computation of a time step is started.

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

"""
TODO: FMI specification reference.

Cancels a DoStep operation.

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

"""
TODO: FMI specification reference.

Informs Master of fmi2Status.

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

"""
TODO: FMI specification reference.

Informs Master of fmi2Status of fmi2Real.

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

"""
TODO: FMI specification reference.

Informs Master of fmi2Status of fmi2Integer.

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

"""
TODO: FMI specification reference.

Informs Master of fmi2Status of fmi2Boolean.

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

"""
TODO: FMI specification reference.

Informs Master of fmi2Status of fmi2String.

For more information call ?fmi2GetStringStatus
"""
function fmi2GetStringStatus(fmu::FMU2, s::fmi2StatusKind, value::String)
    fmi2GetStringStatus(fmu.components[end], s, fmiString(value))
end

"""
TODO: FMI specification reference.

Set independent variable time and reinitialize chaching of variables that depend on time.

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

"""
TODO: FMI specification reference.

Set a new (continuous) state vector and reinitialize chaching of variables that depend on states.

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

"""
TODO: FMI specification reference.

The model enters Event Mode.

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

"""
TODO: FMI specification reference.

Increment the super dense time in event mode.

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

"""
TODO: FMI specification reference.

The model enters Continuous-Time Mode.

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

"""
TODO: FMI specification reference.

This function must be called by the environment after every completed step
If enterEventMode == fmi2True, the event mode must be entered
If terminateSimulation == fmi2True, the simulation shall be terminated

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

"""
TODO: FMI specification reference.

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

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

"""
TODO: FMI specification reference.

Returns the event indicators of the FMU.

For more information call ?fmi2GetEventIndicators
"""
function fmi2GetEventIndicators(fmu::FMU2)
    ni = Csize_t(fmu.modelDescription.numberOfEventIndicators)
    eventIndicators = zeros(fmi2Real, ni)
    fmi2GetEventIndicators(fmu.components[end], eventIndicators, ni)
    eventIndicators
end

"""
TODO: FMI specification reference.

Return the new (continuous) state vector x.

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

"""
TODO: FMI specification reference.

Return the nominal values of the continuous states.

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

"""
Returns the start/default value for a given value reference.

TODO: Add this command in the documentation.
"""
function fmi2GetStartValue(fmu::FMU2, vr::fmi2ValueReferenceFormat)
    fmi2GetStartValue(fmu.components[end], vr)
end


text:
function solve_for_velocities!(state::MechanismState, momentum::Momentum, fixedJoints::Vector{GenericJoint{Float64}})
    solver = Gurobi.GurobiSolver()
    Gurobi.setparameters!(solver, Silent = true)
    model = Model(solver = solver)
    nv = num_velocities(state)
    A = momentum_matrix(state)

    momentum = transform(state, momentum, A.frame)

    @variable(model, v[1 : nv])
    @constraint(model, Array(A.linear) * v .== Array(momentum.linear))
    @constraint(model, Array(A.angular) * v .== Array(momentum.angular))
    for joint in fixedJoints
        @constraint(model, v[velocity_range(state, joint)] .== 0)
    end
    @objective(model, Min, sum(v[i]^2 for i = 1 : nv))

    status = solve(model)
    set_velocity!(state, getvalue(v))
end

function initialize_state!(state::MechanismState, xd::Float64, positionControlledJoints)
    setdirty!(state)
    com = center_of_mass(state)
    centroidalFrame = CartesianFrame3D()
    comdot = FreeVector3D(centroidalFrame, SVector(xd, 0.0, 0.0))
    linearMomentum = mass(state.mechanism) * comdot
    initialMomentum = Momentum{Float64}(centroidalFrame, zero(SVector{3}), linearMomentum.v)
    centroidalToWorld = Transform3D(centroidalFrame, root_frame(state.mechanism), com.v)
    initialMomentum = transform(initialMomentum, centroidalToWorld)
    solve_for_velocities!(state, initialMomentum, positionControlledJoints)
end

ext: jl

max_stars_repo_name: tkoolen/VariableHeightMomentumControl

max_issues_repo_name: tkoolen/VariableHeightMomentumControl

max_forks_repo_name: tkoolen/VariableHeightMomentumControl


text:
module GaussianRandomFields

# dependencies
using SpecialFunctions, Plots, FastGaussQuadrature, Arpack, FFTW

using LinearAlgebra, Statistics

using Base.Cartesian

# import statements
import Plots: plot, plot!, surface, heatmap, contour, contourf

# export statements
export IsotropicCovarianceStructure, AnisotropicCovarianceStructure, CovarianceFunction, apply # from covariance_functions.jl

export SeparableCovarianceFunction # from separable.jl

export Matern # from matern.jl

export Exponential # from exponential.jl

export SquaredExponential, Gaussian # from squaredexponential.jl

export AnisotropicExponential # from anisotropic_exponential.jl

export GaussianRandomField, sample, randdim # from gaussian_random_fields.jl

export Cholesky # from cholesky.jl

export Spectral # from spectral.jl

export KarhunenLoeve, rel_error # from karhunen_loeve.jl

export CirculantEmbedding # from circulant_embedding.jl

export GaussLegendre, EOLE, Simpson, Midpoint, Trapezoidal, AbstractEigenSolver, EigsSolver, EigenSolver, compute # from quadrature.jl

export plot, plot!, heatmap, surface, contour, contourf, plot_eigenvalues, plot_eigenfunction, plot_covariance_matrix # from plots.jl

export star, Lshape # from data/

# include statements
include("covariance_functions/covariance_functions.jl")

include("covariance_functions/matern.jl")

include("covariance_functions/exponential.jl")

include("covariance_functions/squared_exponential.jl")

include("covariance_functions/anisotropic_exponential.jl")

include("gaussian_random_fields.jl")

include("fem.jl")

include("generators/cholesky.jl")

include("generators/spectral.jl")

include("generators/quadrature.jl")

include("generators/karhunen_loeve.jl")

include("generators/circulant_embedding.jl")

include("analytic.jl")

include("separable.jl")

include("../data/star.jl")

include("../data/Lshape.jl")

include("plot.jl")

end # module


text:
# generator.jl
# generator.jl uses Clang to parse zdict.h, zstd.h, zstd_errors.h and generates LibZstd.jl

# Clang.jl 0.14 (40e3b903)
# As of 2021-09-18 0.14 is currently unreleased
# See https://github.com/JuliaInterop/Clang.jl/milestone/1

# Configuration files:
# generator.toml is used to configure the Clang generator
# prologue.jl adds some definitions to beginning LibZstd.jl

# Directions:
# 1. Change to this directory e.g. `import CodecZstd; cd(joinpath(dirname(pathof(CodecZstd)), "libzstd"))`
# 2. Activate the project in this directory e.g. `using Pkg; Pkg.activate(".")`
# 2a. You may need to dev Clang to get 0.14 or greater: `Pkg.dev("Clang")`
# 3. Include this file: `include("generator.jl")``
# Alternatively, from the shell, `cd libzstd; julia --project=. generator.jl`

#=
# Example Code Block to Run Generator
import CodecZstd; cd(joinpath(dirname(pathof(CodecZstd)), "libzstd"))
using Pkg; Pkg.activate(".")
# Pkg.dev("Clang")
include("generator.jl")
=#

using Clang.Generators
using Zstd_jll

cd(@__DIR__)

include_dir = joinpath(Zstd_jll.artifact_dir, "include") |> normpath
options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")
push!(args, "-DZSTD_STATIC_LINKING_ONLY")

headers = [joinpath(include_dir, header) for header in readdir(include_dir) if endswith(header, ".h")]

# create context
ctx = create_context(headers, args, options)

# build without printing so we can do custom rewriting
build!(ctx, BUILDSTAGE_NO_PRINTING)

# custom rewriter
function rewrite!(e::Expr)
    if e.head == :struct && length(e.args) > 1 && (e.args[2] == :ZSTD_inBuffer_s || e.args[2] == :ZSTD_outBuffer_s)
        # Make Zstd_inBuffer and Zstd_outBuffer mutable
        e.args[1] = true
    end
    e
end

function rewrite!(dag::ExprDAG)
    for node in get_nodes(dag)
        if node.id == :ZSTD_WINDOWLOG_MAX
            # Unsupported macro
            node.exprs[1].args[1] = :(const ZSTD_WINDOWLOG_MAX = sizeof(Csize_t) == 4 ? ZSTD_WINDOWLOG_MAX_32 : ZSTD_WINDOWLOG_MAX_64)
        elseif node.id == :ZSTD_CHAINLOG_MAX
            # Unsupported macro
            node.exprs[1].args[1] = :(const ZSTD_CHAINLOG_MAX = sizeof(Csize_t) == 4 ? ZSTD_CHAINLOG_MAX_32 : ZSTD_CHAINLOG_MAX_64)
        elseif node.id == :ZSTD_LIB_VERSION
            # Unsupported macro
            node.exprs[1].args[1] = :(const ZSTD_LIB_VERSION = VersionNumber(ZSTD_VERSION_MAJOR, ZSTD_VERSION_MINOR, ZSTD_VERSION_RELEASE))
        elseif node isa ExprNode{Generators.MacroUnsupported}
            # Usually just export symbols for ZSTDLIB_API, ZDICTLIB_API, ZSTDERRORLIB_API
            node.exprs[1].args[1] = :(const $(node.id) = nothing)
        end
        for expr in get_exprs(node)
            rewrite!(expr)
        end
    end
end

rewrite!(ctx.dag)

# print
build!(ctx, BUILDSTAGE_PRINTING_ONLY)


text:
using Printf

function monteπ(n)
    s = count(rand() ^ 2 + rand() ^ 2 < 1 for _ in 1:n)
    return 4s / n
end

for n in 10 .^ (3:8)
    p = monteπ(n)
    println("$(lpad(n, 9)): π ≈ $(lpad(p, 10)), pct.err = ", @sprintf("%2.5f%%", abs(p - π) / π))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@doc raw"""
    GeneralizedGrassmann{n,k,𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultEmbeddingType}

The generalized Grassmann manifold $\operatorname{Gr}(n,k,B)$ consists of all subspaces
spanned by $k$ linear independent vectors $𝔽^n$, where $𝔽  ∈ \{ℝ, ℂ\}$ is either the real- (or complex-) valued vectors.
This yields all $k$-dimensional subspaces of $ℝ^n$ for the real-valued case and all $2k$-dimensional subspaces
of $ℂ^n$ for the second.

The manifold can be represented as

````math
\operatorname{Gr}(n, k, B) := \bigl\{ \operatorname{span}(p)\ \big|\ p ∈ 𝔽^{n × k}, p^\mathrm{H}Bp = I_k\},
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate (or Hermitian) transpose and
$I_k$ is the $k × k$ identity matrix. This means, that the columns of $p$
form an unitary basis of the subspace with respect to the scaled inner product, that is a
point on $\operatorname{Gr}(n,k,B)$, and hence the subspace can actually be represented by
a whole equivalence class of representers. For $B=I_n$ this simplifies to the [`Grassmann`](@ref) manifold.

The tangent space at a point (subspace) $p$ is given by

````math
T_x\mathrm{Gr}(n,k,B) = \bigl\{
X ∈ 𝔽^{n × k} :
X^{\mathrm{H}}Bp + p^{\mathrm{H}}BX = 0_{k} \bigr\},
````

where $0_{k}$ denotes the $k × k$ zero matrix.

Note that a point $p ∈ \operatorname{Gr}(n,k,B)$ might be represented by
different matrices (i.e. matrices with $B$-unitary column vectors that span
the same subspace). Different representations of $p$ also lead to different
representation matrices for the tangent space $T_p\mathrm{Gr}(n,k,B)$

The manifold is named after
[Hermann G. Graßmann](https://en.wikipedia.org/wiki/Hermann_Grassmann) (1809-1877).

# Constructor

    GeneralizedGrassmann(n, k, B=I_n, field=ℝ)

Generate the (real-valued) Generalized Grassmann manifold of $n\times k$ dimensional
orthonormal matrices with scalar product `B`.
"""
struct GeneralizedGrassmann{n,k,𝔽,TB<:AbstractMatrix} <:
       AbstractEmbeddedManifold{𝔽,DefaultEmbeddingType}
    B::TB
end

function GeneralizedGrassmann(
    n::Int,
    k::Int,
    B::AbstractMatrix=Matrix{Float64}(I, n, n),
    field::AbstractNumbers=ℝ,
)
    return GeneralizedGrassmann{n,k,field,typeof(B)}(B)
end

@doc raw"""
    check_point(M::GeneralizedGrassmann{n,k,𝔽}, p)

Check whether `p` is representing a point on the [`GeneralizedGrassmann`](@ref) `M`, i.e. its
a `n`-by-`k` matrix of unitary column vectors with respect to the B inner prudct and
of correct `eltype` with respect to `𝔽`.
"""
function check_point(M::GeneralizedGrassmann{n,k,𝔽}, p; kwargs...) where {n,k,𝔽}
    mpv = invoke(
        check_point,
        Tuple{typeof(get_embedding(M)),typeof(p)},
        get_embedding(M),
        p;
        kwargs...,
    )
    mpv === nothing || return mpv
    c = p' * M.B * p
    if !isapprox(c, one(c); kwargs...)
        return DomainError(
            norm(c - one(c)),
            "The point $(p) does not lie on $(M), because x'Bx is not the unit matrix.",
        )
    end
    return nothing
end

@doc raw"""
    check_vector(M::GeneralizedGrassmann{n,k,𝔽}, p, X; kwargs...)

Check whether `X` is a tangent vector in the tangent space of `p` on
the [`GeneralizedGrassmann`](@ref) `M`, i.e. that `X` is of size and type as well as that

````math
    p^{\mathrm{H}}BX + \overline{X^{\mathrm{H}}Bp} = 0_k,
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transpose or Hermitian,
$\overline{\cdot}$ the (elementwise) complex conjugate, and $0_k$ denotes the $k × k$ zero natrix.
"""
function check_vector(M::GeneralizedGrassmann{n,k,𝔽}, p, X; kwargs...) where {n,k,𝔽}
    mpv = invoke(
        check_vector,
        Tuple{typeof(get_embedding(M)),typeof(p),typeof(X)},
        get_embedding(M),
        p,
        X;
        kwargs...,
    )
    mpv === nothing || return mpv
    if !isapprox(p' * M.B * X, -conj(X' * M.B * p); kwargs...)
        return DomainError(
            norm(p' * M.B * X + conj(X' * M.B * p)),
            "The matrix $(X) does not lie in the tangent space of $(p) on $(M), since x'Bv + v'Bx is not the zero matrix.",
        )
    end
    return nothing
end

function decorated_manifold(M::GeneralizedGrassmann{N,K,𝔽}) where {N,K,𝔽}
    return Euclidean(N, K; field=𝔽)
end

@doc raw"""
    distance(M::GeneralizedGrassmann, p, q)

Compute the Riemannian distance on [`GeneralizedGrassmann`](@ref)
manifold `M`$= \mathrm{Gr}(n,k,B)$.

Let $USV = p^\mathrm{H}Bq$ denote the SVD decomposition of
$p^\mathrm{H}Bq$, where $\cdot^{\mathrm{H}}$ denotes the complex
conjugate transposed or Hermitian. Then the distance is given by
````math
d_{\mathrm{Gr}(n,k,B)}(p,q) = \operatorname{norm}(\operatorname{Re}(b)).
````
where

````math
b_{i}=\begin{cases}
0 & \text{if} \; S_i ≥ 1\\
\arccos(S_i) & \, \text{if} \; S_i<1.
\end{cases}
````
"""
function distance(M::GeneralizedGrassmann, p, q)
    p ≈ q && return zero(real(eltype(p)))
    a = svd(p' * M.B * q).S
    return sqrt(sum(x -> abs2(acos(clamp(x, -1, 1))), a))
end

@doc raw"""
    exp(M::GeneralizedGrassmann, p, X)

Compute the exponential map on the [`GeneralizedGrassmann`](@ref) `M`$= \mathrm{Gr}(n,k,B)$
starting in `p` with tangent vector (direction) `X`. Let $X^{\mathrm{H}}BX = USV$ denote the
SVD decomposition of $X^{\mathrm{H}}BX$. Then the exponential map is written using

````math
\exp_p X = p V\cos(S)V^\mathrm{H} + U\sin(S)V^\mathrm{H},
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian and the
cosine and sine are applied element wise to the diagonal entries of $S$.
"""
exp(::GeneralizedGrassmann, ::Any...)

function exp!(M::GeneralizedGrassmann, q, p, X)
    norm(M, p, X) ≈ 0 && return copyto!(q, p)
    d = svd(X' * M.B * X)
    V = d.Vt
    S = abs.(sqrt.(d.S))
    U = X * (V / Diagonal(S))
    z = p * V * Diagonal(cos.(S)) * V + U * Diagonal(sin.(S)) * V
    return copyto!(q, project(M, z))
end

@doc raw"""
    injectivity_radius(M::GeneralizedGrassmann)
    injectivity_radius(M::GeneralizedGrassmann, p)

Return the injectivity radius on the [`GeneralizedGrassmann`](@ref) `M`,
which is $\frac{π}{2}$.
"""
injectivity_radius(::GeneralizedGrassmann) = π / 2
injectivity_radius(::GeneralizedGrassmann, ::ExponentialRetraction) = π / 2
injectivity_radius(::GeneralizedGrassmann, ::Any) = π / 2
injectivity_radius(::GeneralizedGrassmann, ::Any, ::ExponentialRetraction) = π / 2
eval(
    quote
        @invoke_maker 1 AbstractManifold injectivity_radius(
            M::GeneralizedGrassmann,
            rm::AbstractRetractionMethod,
        )
    end,
)

@doc raw"""
    inner(M::GeneralizedGrassmann, p, X, Y)

Compute the inner product for two tangent vectors `X`, `Y` from the tangent space
of `p` on the [`GeneralizedGrassmann`](@ref) manifold `M`. The formula reads

````math
g_p(X,Y) = \operatorname{tr}(X^{\mathrm{H}}BY),
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
"""
inner(M::GeneralizedGrassmann{n,k}, p, X, Y) where {n,k} = dot(X, M.B * Y)

function Base.isapprox(M::GeneralizedGrassmann, p, X, Y; kwargs...)
    return isapprox(sqrt(inner(M, p, zero_vector(M, p), X - Y)), 0; kwargs...)
end
function Base.isapprox(M::GeneralizedGrassmann, p, q; kwargs...)
    return isapprox(distance(M, p, q), 0.0; kwargs...)
end

@doc raw"""
    log(M::GeneralizedGrassmann, p, q)

Compute the logarithmic map on the [`GeneralizedGrassmann`](@ref) `M`$ = \mathcal M=\mathrm{Gr}(n,k,B)$,
i.e. the tangent vector `X` whose corresponding [`geodesic`](@ref) starting from `p`
reaches `q` after time 1 on `M`. The formula reads

````math
\log_p q = V\cdot \operatorname{atan}(S) \cdot U^\mathrm{H},
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian.
The matrices $U$ and $V$ are the unitary matrices, and $S$ is the diagonal matrix
containing the singular values of the SVD-decomposition

````math
USV = (q^\mathrm{H}Bp)^{-1} ( q^\mathrm{H} - q^\mathrm{H}Bpp^\mathrm{H}).
````

In this formula the $\operatorname{atan}$ is meant elementwise.
"""
log(::GeneralizedGrassmann, ::Any...)

function log!(M::GeneralizedGrassmann{n,k}, X, p, q) where {n,k}
    z = q' * M.B * p
    At = q' - z * p'
    Bt = z \ At
    d = svd(Bt')
    return X .= view(d.U, :, 1:k) * Diagonal(atan.(view(d.S, 1:k))) * view(d.Vt, 1:k, :)
end

@doc raw"""
    manifold_dimension(M::GeneralizedGrassmann)

Return the dimension of the [`GeneralizedGrassmann(n,k,𝔽)`](@ref) manifold `M`, i.e.

````math
\dim \operatorname{Gr}(n,k,B) = k(n-k) \dim_ℝ 𝔽,
````

where $\dim_ℝ 𝔽$ is the [`real_dimension`](@ref) of `𝔽`.
"""
function manifold_dimension(::GeneralizedGrassmann{n,k,𝔽}) where {n,k,𝔽}
    return k * (n - k) * real_dimension(𝔽)
end

"""
    mean(
        M::GeneralizedGrassmann,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolationWithinRadius(π/4);
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::AbstractManifold, args...)) of `x` using
[`GeodesicInterpolationWithinRadius`](@ref).
"""
mean(::GeneralizedGrassmann{n,k} where {n,k}, ::Any...)

function Statistics.mean!(
    M::GeneralizedGrassmann{n,k},
    p,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
) where {n,k}
    return mean!(M, p, x, w, GeodesicInterpolationWithinRadius(π / 4); kwargs...)
end

@doc raw"""
    project(M::GeneralizedGrassmann, p)

Project `p` from the embedding onto the [`GeneralizedGrassmann`](@ref) `M`, i.e. compute `q`
as the polar decomposition of $p$ such that $q^{\mathrm{H}}Bq$ is the identity,
where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transpose.
"""
project(::GeneralizedGrassmann, ::Any)

function project!(M::GeneralizedGrassmann, q, p)
    s = svd(p)
    e = eigen(s.U' * M.B * s.U)
    qsinv = e.vectors * Diagonal(1 ./ sqrt.(e.values))
    q .= s.U * qsinv * e.vectors' * s.V'
    return q
end

@doc raw"""
    project(M::GeneralizedGrassmann, p, X)

Project the `n`-by-`k` `X` onto the tangent space of `p` on the
[`GeneralizedGrassmann`](@ref) `M`, which is computed by

````math
\operatorname{proj_p}(X) = X - pp^{\mathrm{H}}B^\mathrm{T}X,
````

where $\cdot^{\mathrm{H}}$ denotes the complex conjugate transposed or Hermitian
and $\cdot^{\mathrm{T}}$ the transpose.
"""
project(::GeneralizedGrassmann, ::Any, ::Any)

function project!(M::GeneralizedGrassmann, Y, p, X)
    A = p' * M.B' * X
    copyto!(Y, X - p * Hermitian((A + A') / 2))
    return Y
end

@doc raw"""
    representation_size(M::GeneralizedGrassmann{n,k})

Return the represenation size or matrix dimension of a point on the [`GeneralizedGrassmann`](@ref)
`M`, i.e. $(n,k)$ for both the real-valued and the complex value case.
"""
@generated representation_size(::GeneralizedGrassmann{n,k}) where {n,k} = (n, k)

@doc raw"""
    retract(M::GeneralizedGrassmann, p, X, ::PolarRetraction)

Compute the SVD-based retraction [`PolarRetraction`](@ref) on the
[`GeneralizedGrassmann`](@ref) `M`, by
[`project`](@ref project(M::GeneralizedGrassmann, p))ing $p + X$ onto `M`.
"""
retract(::GeneralizedGrassmann, ::Any, ::Any, ::PolarRetraction)

function retract!(M::GeneralizedGrassmann, q, p, X, ::PolarRetraction)
    project!(M, q, p + X)
    return q
end
function retract!(M::GeneralizedGrassmann, q, p, X, ::ProjectionRetraction)
    project!(M, q, p + X)
    return q
end

function Base.show(io::IO, M::GeneralizedGrassmann{n,k,𝔽}) where {n,k,𝔽}
    return print(io, "GeneralizedGrassmann($(n), $(k), $(M.B), $(𝔽))")
end

@doc raw"""
    vector_transport_to(M::GeneralizedGrassmann, p, X, q, ::ProjectionTransport)

Compute the vector transport of the tangent vector `X` at `p` to `q`,
using the [`project`](@ref project(::GeneralizedGrassmann, ::Any...))
of `X` to `q`.
"""
vector_transport_to(::GeneralizedGrassmann, ::Any, ::Any, ::Any, ::ProjectionTransport)

function vector_transport_to!(M::GeneralizedGrassmann, Y, p, X, q, ::ProjectionTransport)
    project!(M, Y, q, X)
    return Y
end

@doc raw"""
    zero_vector(M::GeneralizedGrassmann, p)

Return the zero tangent vector from the tangent space at `p` on the
[`GeneralizedGrassmann`](@ref) `M`, which is given by a zero matrix the same size as `p`.
"""
zero_vector(::GeneralizedGrassmann, ::Any...)

zero_vector!(::GeneralizedGrassmann, X, p) = fill!(X, 0)


text:
using GeometryBasics: Point, Rect3D, Sphere, GeometryPrimitive, origin, widths, radius
using LinearAlgebra: dot, norm
using Random: rand

"""
`Position` : A point in 3-D.  Ultimately, derived from StaticArray.
"""
const Position = Point{3,Float64}

"""
The MonteCarlo uses the shapes defined in GeometryBasics basics as the foundation for its
sample construction mechanisms.  However, GeometryBasics basics does not provide all the
necessary methods.  Three additional methods are

    isinside(r::Shape, pos::Position)

Is `pos` strictly inside `r`?

    intersection(r::Shape, pos0::Particle, pos1::Particle)::Float64

Return a number `f` which represent the fraction of the distance from `pos0` to `pos1` that
first intersects the `Shape` `r`.  The intersection point will equal `pos0 .+ f*(pos1 .- pos0)`.
If `f` is between 0.0 and 1.0 then the intersection is on the interval between `pos0` and `pos1`.
If the ray from `pos0` towards `pos2` does not intersect `r` then this function returns Inf64.
"""
const RectangularShape = Rect3D{Float64}

isinside(rr::RectangularShape, pos::AbstractArray{Float64}) =
    all(pos .> minimum(rr)) && all(pos .< maximum(rr))

function intersection(
    rr::RectangularShape,
    pos1::AbstractArray{Float64},
    pos2::AbstractArray{Float64},
)::Float64
    _between(a, b, c) = (a > b) && (a < c)
    t = Inf64
    corner1, corner2 = minimum(rr), maximum(rr)
    for i in eachindex(pos1)
        j, k = i % 3 + 1, (i + 1) % 3 + 1
        if pos2[i] != pos1[i]
            u = (corner1[i] - pos1[i]) / (pos2[i] - pos1[i])
            if (u > 0.0) &&
               (u <= t) && #
               _between(pos1[j] + u * (pos2[j] - pos1[j]), corner1[j], corner2[j]) && #
               _between(pos1[k] + u * (pos2[k] - pos1[k]), corner1[k], corner2[k])
                t = u
            end
            u = (corner2[i] - pos1[i]) / (pos2[i] - pos1[i])
            if (u > 0.0) &&
               (u <= t) && #
               _between(pos1[j] + u * (pos2[j] - pos1[j]), corner1[j], corner2[j]) && #
               _between(pos1[k] + u * (pos2[k] - pos1[k]), corner1[k], corner2[k])
                t = u
            end
        end
    end
    return t
end

const SphericalShape = Sphere{Float64}

isinside(sr::SphericalShape, pos::AbstractArray{Float64}) =
    norm(pos .- origin(sr)) < radius(sr)

function intersection(
    sr::SphericalShape,
    pos0::AbstractArray{Float64},
    pos1::AbstractArray{Float64},
)::Float64
    d, m = pos1 .- pos0, pos0 .- origin(sr)
    ma2, b = -2.0 * dot(d, d), 2.0 * dot(m, d)
    f = b^2 + ma2 * 2.0 * (dot(m, m) - radius(sr)^2)
    if f >= 0.0
        up, un = (b + sqrt(f)) / ma2, (b - sqrt(f)) / ma2
        return min(up < 0.0 ? Inf64 : up, un < 0.0 ? Inf64 : un)
    end
    return Inf64
end

"""
    random_point_inside(shape)

Generate a randomized point that is guaranteed to be in the interior of the shape.
"""
function random_point_inside(shape)::Position
    res = Position(origin(shape) .+ rand(Position) .* widths(shape))
    while !isinside(shape, res)
        res = Position(origin(shape) .+ rand(Position) .* widths(shape))
    end
    return res
end

"""
Particle represents a type that may be simulated using a transport Monte Carlo.  It must provide
these methods:

    position(el::Particle)::Position
    previous(el::Particle)::Position
    energy(el::Particle)::Float64

The position of the current and previous elastic scatter locations which are stored in that Particle type.

    T(prev::Position, curr::Position, energy::Energy) where {T <: Particle }
    T(el::T, 𝜆::Float64, 𝜃::Float64, 𝜑::Float64, ΔE::Float64) where {T <: Particle }

Two constructors: One to create a defined Particle and the other to create a new Particle based off
another which is translated by `λ` at a scattering angle (`θ`, `ϕ`) which energy change of `ΔE`

    transport(pc::T, mat::Material)::NTuple{4, Float64} where {T <: Particle }

A function that generates the values of ( `λ`, `θ`, `ϕ`, `ΔE`) for the specified `Particle` in the specified `Material`.
"""
abstract type Particle end

struct Electron <: Particle
    previous::Position
    current::Position
    energy::Float64 # eV

    """
        Electron(prev::Position, curr::Position, energy::Float64)
        Electron(el::Electron, 𝜆::Float64, 𝜃::Float64, 𝜑::Float64, ΔE::Float64)::Electron

    Create a new `Electron` from this one in which the new `Electron` is a distance `𝜆` from the
    first along a trajectory that is `𝜃` and `𝜑` off the current trajectory.
    """
    Electron(prev::AbstractArray{Float64}, curr::AbstractArray{Float64}, energy::Float64) =
        new(prev, curr, energy)

    function Electron(el::Electron, 𝜆::Float64, 𝜃::Float64, 𝜑::Float64, ΔE::Float64)
        (u, v, w) = LinearAlgebra.normalize(position(el) .- previous(el))
        sc =
            1.0 - abs(w) > 1.0e-8 ? #
            Position( #
                u * cos(𝜃) + sin(𝜃) * (u * w * cos(𝜑) - v * sin(𝜑)) / sqrt(1.0 - w^2), #
                v * cos(𝜃) + sin(𝜃) * (v * w * cos(𝜑) + u * sin(𝜑)) / sqrt(1.0 - w^2), #
                w * cos(𝜃) - sqrt(1.0 - w^2) * sin(𝜃) * cos(𝜑), #
            ) :
            Position( #
                sign(w) * sin(𝜃) * cos(𝜑), #
                sign(w) * sin(𝜃) * sin(𝜑), #
                sign(w) * cos(𝜃),
            )
        return new(position(el), position(el) .+ 𝜆 * sc, el.energy + ΔE)
    end
end

Base.show(io::IO, el::Electron) = print(io, "Electron[$(position(el)), $(energy(el)) eV]")
Base.position(el::Particle) = el.current
previous(el::Particle) = el.previous
energy(el::Particle) = el.energy

"""
    transport(pc::Electron, mat::Material, ecx=Liljequist1989, bethe=JoyLuo)::NTuple{4, Float64}

The default function defining elastic scattering and energy loss for an Electron.

Returns ( `λ`, `θ`, `ϕ`, `ΔE`) where `λ` is the mean path length, `θ` is the elastic scatter angle, `ϕ` is the azimuthal elastic scatter
angle and `ΔE` is the energy loss for transport over the distance `λ`.
"""
function transport(
    pc::Electron,
    mat::Material,
    ecx::Type{<:ElasticScatteringCrossSection} = Liljequist1989,
    bethe::Type{<:BetheEnergyLoss} = JoyLuo,
)::NTuple{4,Float64}
    (𝜆′, θ′, ϕ′) = rand(ecx, mat, pc.energy)
    return (𝜆′, θ′, ϕ′, 𝜆′ * dEds(bethe, pc.energy, mat))
end

"""
    pathlength(el::Particle)

Length of the line segment represented by `el`.
"""
pathlength(el::Particle) = norm(position(el) .- previous(el))

intersection(r, p::Particle) = intersection(r, previous(p), position(p))

"""
    Region

A `Region` combines a geometric primative and a `Material` (with `:Density` property) and may fully contain zero or more child `Region`s.
"""
struct Region
    shape::GeometryPrimitive{3,Float64}
    material::Material
    parent::Union{Nothing,Region}
    children::Vector{Region}
    name::String

    function Region(
        sh::T,
        mat::Material,
        parent::Union{Nothing,Region},
        name::Union{Nothing,String} = nothing,
        ntests = 1000,
    ) where {T}
        @assert mat[:Density] > 0.0
        name = something(
            name,
            isnothing(parent) ? "Root" : "$(parent.name)[$(length(parent.children)+1)]",
        )
        res = new(sh, mat, parent, Region[], name)
        if !isnothing(parent)
            @assert all(
                _ -> isinside(parent.shape, random_point_inside(sh)),
                Base.OneTo(ntests),
            ) "The child $sh is not fully contained within the parent $(parent.shape)."
            @assert all(
                ch -> all(
                    _ -> !isinside(ch.shape, random_point_inside(sh)),
                    Base.OneTo(ntests),
                ),
                parent.children,
            ) "The child $sh overlaps a child of the parent shape."
            push!(parent.children, res)
        else
        end
        return res
    end
end

Base.show(io::IO, reg::Region) = print(
    io,
    "Region[$(reg.name), $(reg.shape), $(reg.material), $(length(reg.children)) children]",
)

"""
    childmost_region(reg::Region, pos::Position)::Region

Find the inner-most `Region` within `reg` containing the point `pos`.
"""
function childmost_region(reg::Region, pos::AbstractArray{Float64})::Region
    res = findfirst(ch -> isinside(ch.shape, pos), reg.children)
    return !isnothing(res) ? childmost_region(reg.children[res], pos) : reg
end

"""
    take_step(p::T, reg::Region, 𝜆::Float64, 𝜃::Float64, 𝜑::Float64)::Tuple{T, Region, Bool} where { T<: Particle}

Returns a `Tuple` containing a new `Particle` and the child-most `Region` in which the new `Particle` is found based
on a scatter event consisting a translation of up to `𝜆` mean-free path along a new direction given relative
to the current direction of `p` via the scatter angles `𝜃` and `𝜑`.

Returns the updated `Particle` reflecting the last trajectory step and the Region for the next step.
"""
function take_step(
    p::T,
    reg::Region,
    𝜆::Float64,
    𝜃::Float64,
    𝜑::Float64,
    ΔE::Float64,
    ϵ::Float64 = 1.0e-12,
)::Tuple{T,Region,Bool} where {T<:Particle}
    newP, nextReg = T(p, 𝜆, 𝜃, 𝜑, ΔE), reg
    t = min(
        intersection(reg.shape, newP), # Leave this Region?
        (intersection(ch.shape, newP) for ch in reg.children)..., # Enter a new child Region?
    )
    scatter = t > 1.0
    if !scatter # Enter new region
        newP = T(p, (t + ϵ) * 𝜆, 𝜃, 𝜑, (t + ϵ) * ΔE)
        nextReg = childmost_region(isnothing(reg.parent) ? reg : reg.parent, position(newP))
    end
    return (newP, nextReg, scatter)
end

"""
trajectory(p::T, reg::Region, scf::Function=transport; minE::Float64=50.0) where {T <: Particle}
trajectory(eval::Function, p::T, reg::Region, scf::Function, terminate::Function) where { T <: Particle }

Run a single particle trajectory from `p` to `minE` or until the particle exits `reg`.

  * `eval(part::T, region::Region)` a function evaluated at each scattering point
  * `p` defines the initial position, direction and energy of the particle (often created with `gun(T, ...)`)
  * `reg` The outer-most region for the trajectory (usually created with `chamber()`)
  * `scf` A function from (<:Particle, Material) -> ( λ, θ, ϕ, ΔE ) that implements the transport dynamics
  * `minE` Stopping criterion
  * `terminate` a function taking `T` and `Region` that returns false except on the last step (like `terminate = (pc,r)->pc.energy < 50.0`)
"""
function trajectory(
    eval::Function,
    p::T,
    reg::Region,
    scf::Function,
    terminate::Function,
) where {T<:Particle}
    (pc, nextr) = (p, childmost_region(reg, position(p)))
    θ, ϕ = 0.0, 0.0
    while (!terminate(pc, reg)) && isinside(reg.shape, position(pc))
        prevr = nextr
        (λ, θₙ, ϕₙ, ΔZ) = scf(pc, nextr.material)
        (pc, nextr, scatter) = take_step(pc, nextr, λ, θ, ϕ, ΔZ)
        (θ, ϕ) = scatter ? (θₙ, ϕₙ) : (0.0, 0.0)
        eval(pc, prevr)
    end
end

function trajectory(
    eval::Function,
    p::T,
    reg::Region,
    scf::Function = (t::T, mat::Material) -> transport(t, mat);
    minE::Float64 = 50.0,
) where {T<:Particle}
    term(pc::T, _::Region) = pc.energy < minE
    trajectory(eval, p, reg, scf, term)
end


text:
abstract AbstractMaterial

abstract AbstractMaterialStatus

stiffness() = error("Not implemented")
create_matstat() = error("Not implemented")
stress() = error("Not implemented")
get_kalpha() = error("Not implemented")
get_kalphas() = error("Not implemented")

@lazymod LinearIsotropicMod "linear_isotropic.jl"


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

JLLWrappers.@generate_main_file_header("licensecheck")
JLLWrappers.@generate_main_file("licensecheck", UUID("4ecb348a-8b88-51ea-b912-4c460483ee91"))
end  # module licensecheck_jll