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

# What is included in the file `FMI3_c.jl`?
# - default implementations for the `fmi3CallbackFunctions`
# - enum definition for `fmi2ComponentState` to protocol the current FMU-state
# - julia-implementaions of the functions inside the FMI-specification
# Any c-function `f(c::fmi2Component, args...)` in the spec is implemented as `f(c::fmi2Component, args...)`.
# Any c-function `f(args...)` without a leading `fmi2Component`-arguemnt is implented as `f(c_ptr, args...)` where `c_ptr` is a pointer to the c-function (inside the DLL).

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.1. Super State: FMU State Setable

Function that is called in the FMU, usually if an fmi3XXX function, does not behave as desired. If “logger” is called with “status = fmi3OK”, then the message is a pure information message.
instanceEnvironment - is the instance name of the model that calls this function.
category - is the category of the message. The meaning of “category” is defined by the modeling environment that generated the FMU. Depending on this modeling environment, none, some or all allowed values of “category” for this FMU are defined in the modelDescription.xml file via element “<fmiModelDescription><LogCategories>”, see section 2.4.5. Only messages are provided by function logger that have a category according to a call to fmi3SetDebugLogging (see below).
message - is provided in the same way and with the same format control as in function “printf” from the C standard library. [Typically, this function prints the message and stores it optionally in a log file.]
"""
# TODO error in the specification??
function fmi3CallbackLogger(instanceEnvironment::Ptr{Cvoid},
            category::Ptr{Cchar},
            status::Cuint,
            message::Ptr{Cchar})

    if message != C_NULL
        _message = unsafe_string(message)
    else
        _message = ""
    end

    if category != C_NULL
        _category = unsafe_string(category)
    else
        _category = "No category"
    end
    _status = fmi3StatusString(status)
    if status == Integer(fmi3StatusOK)
        @info "[$_status][$_category]: $_message"
    elseif status == Integer(fmi3StatusWarning)
        @warn "[$_status][$_category]: $_message"
    else
        @error "[$_status][$_category]: $_message"
    end

    nothing
end

"""
Source: FMISpec3.0, Version D5ef1c1: 4.2.2. State: Intermediate Update Mode

When a Co-Simulation FMU provides values for its output variables at intermediate points between two consecutive communication points, and is able to receive new values for input variables at these intermediate points, the Intermediate Update Callback function is called. This is typically required when the FMU uses a numerical solver to integrate the FMU's internal state between communication points in fmi3DoStep.
The callback function switches the FMU from Step Mode (see 4.2.1.) in the Intermediate Update Mode (see 4.2.2.) and returns to Step Mode afterwards. The parameters of this function are:

instanceEnvironment - is the instance name of the model that calls this function.

intermediateUpdateTime - is the internal value of the independent variable [typically simulation time] of the FMU at which the callback has been called for intermediate and final steps. If an event happens or an output Clock ticks, intermediateUpdateTime is the time of event or output Clock tick. In Co-Simulation, intermediateUpdateTime is restricted by the arguments to fmi3DoStep as follows:
currentCommunicationPoint ≤ intermediateUpdateTime ≤ (currentCommunicationPoint + communicationStepSize).
The FMU must not call the callback function fmi3CallbackIntermediateUpdate with an intermediateUpdateTime that is smaller than the intermediateUpdateTime given in a previous call of fmi3CallbackIntermediateUpdate with intermediateStepFinished == fmi3True.

If intermediateVariableSetRequested == fmi3True, the co-simulation algorithm may provide intermediate values for continuous input variables with intermediateUpdate = true by calling fmi3Set{VariableType}. The set of variables for which the co-simulation algorithm will provide intermediate values is declared through the requiredIntermediateVariables argument to fmi3InstantiateXXX. If a co-simulation algorithm does not provide a new value for any of the variables contained in the set it registered, the last value set remains.

If intermediateVariableGetAllowed == fmi3True, the co-simulation algorithm may collect intermediate output variables by calling fmi3Get{VariableType} for variables with intermediateUpdate = true. The set of variables for which the co-simulation algorithm can get values is supplied through the requiredIntermediateVariables argument to fmi3InstantiateXXX.

If intermediateStepFinished == fmi3False, the intermediate outputs of the FMU that the co-simulation algorithm inquires with fmi3Get{VariableType} resulting from tentative internal solver states and may still change for the same intermediateUpdateTime [e.g., if the solver deems the tentative state to cause a too high approximation error, it may go back in time and try to re-estimate the state using smaller internal time steps].
If intermediateStepFinished == fmi3True, intermediate outputs inquired by the co-simulation algorithm with fmi3Get{VariableType} correspond to accepted internal solver step.

When canReturnEarly == fmi3True the FMU signals to the co-simulation algorithm its ability to return early from the current fmi3DoStep.

earlyReturnRequested - If and only if canReturnEarly == fmi3True, the co-simulation algorithm may request the FMU to return early from fmi3DoStep by setting earlyReturnRequested == fmi3True.

earlyReturnTime is used to signal the FMU at which time to return early from the current fmi3DoStep, if the return value of earlyReturnRequested == fmi3True. If the earlyReturnTime is greater than the last signaled intermediateUpdateTime, the FMU may integrate up to the time instant earlyReturnTime.

If the ModelDescription has the "providesIntermediateUpdate" flag, the Intermediate update callback function is called. That flag is ignored in ModelExchange and ScheduledExecution.
"""
function fmi3CallbackIntermediateUpdate(instanceEnvironment::Ptr{Cvoid},
    intermediateUpdateTime::fmi3Float64,
    intermediateVariableSetRequested::fmi3Boolean,
    intermediateVariableGetAllowed::fmi3Boolean,
    intermediateStepFinished::fmi3Boolean,
    canReturnEarly::fmi3Boolean,
    earlyReturnRequested::Ptr{fmi3Boolean},
    earlyReturnTime::Ptr{fmi3Float64})
    @debug "To be implemented!"
end

"""
Source: FMISpec3.0, Version D5ef1c1: 5.2.2. State: Clock Activation Mode

A model partition of a Scheduled Execution FMU calls fmi3CallbackClockUpdate to signal that a triggered output Clock ticked or a new interval for a countdown Clock is available.
fmi3CallbackClockUpdate switches the FMU itself then into the Clock Update Mode (see 5.2.3.). The callback may be called from several model partitions.

instanceEnvironment - is the instance name of the model that calls this function.
"""
function fmi3CallbackClockUpdate(instanceEnvironment::Ptr{Cvoid})
    @debug "to be implemented!"
end

# this is a custom type to catch the internal state of the component
@enum fmi3InstanceState begin
    # ToDo
    fmi3InstanceStateToDo
end

"""
Source: FMISpec3.0, Version D5ef1c1:: 2.3.1. Super State: FMU State Setable

This function instantiates a Model Exchange FMU (see Section 3). It is allowed to call this function only if modelDescription.xml includes a <ModelExchange> element.
"""
function fmi3InstantiateModelExchange(cfunc::Ptr{Nothing},
        instanceName::fmi3String,
        fmuinstantiationToken::fmi3String,
        fmuResourceLocation::fmi3String,
        visible::fmi3Boolean,
        loggingOn::fmi3Boolean,
        instanceEnvironment::fmi3InstanceEnvironment,
        logMessage::Ptr{Cvoid})

    compAddr = ccall(cfunc,
        Ptr{Cvoid},
        (Cstring, Cstring, Cstring,
        Cuint, Cuint, Ptr{Cvoid}, Ptr{Cvoid}),
        instanceName, fmuinstantiationToken, fmuResourceLocation,
        visible, loggingOn, instanceEnvironment, logMessage)

    compAddr
end

"""
Source: FMISpec3.0, Version D5ef1c1:: 2.3.1. Super State: FMU State Setable

This function instantiates a Co-Simulation FMU (see Section 4). It is allowed to call this function only if modelDescription.xml includes a <CoSimulation> element.
"""
function fmi3InstantiateCoSimulation(cfunc::Ptr{Nothing},
    instanceName::fmi3String,
    instantiationToken::fmi3String,
    resourcePath::fmi3String,
    visible::fmi3Boolean,
    loggingOn::fmi3Boolean,
    eventModeUsed::fmi3Boolean,
    earlyReturnAllowed::fmi3Boolean,
    requiredIntermediateVariables::Array{fmi3ValueReference},
    nRequiredIntermediateVariables::Csize_t,
    instanceEnvironment::fmi3InstanceEnvironment,
    logMessage::Ptr{Cvoid},
    intermediateUpdate::Ptr{Cvoid})

    compAddr = ccall(cfunc,
        Ptr{Cvoid},
        (Cstring, Cstring, Cstring,
        Cint, Cint, Cint, Cint, Ptr{fmi3ValueReference}, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        instanceName, instantiationToken, resourcePath,
        visible, loggingOn, eventModeUsed, earlyReturnAllowed, requiredIntermediateVariables,
        nRequiredIntermediateVariables, instanceEnvironment, logMessage, intermediateUpdate)

    compAddr
end

# TODO not tested
"""
Source: FMISpec3.0, Version D5ef1c1:: 2.3.1. Super State: FMU State Setable

This function instantiates a Scheduled Execution FMU (see Section 4). It is allowed to call this function only if modelDescription.xml includes a <ScheduledExecution> element.
"""
function fmi3InstantiateScheduledExecution(cfunc::Ptr{Nothing},
    instanceName::fmi3String,
    instantiationToken::fmi3String,
    resourcePath::fmi3String,
    visible::fmi3Boolean,
    loggingOn::fmi3Boolean,
    instanceEnvironment::fmi3InstanceEnvironment,
    logMessage::Ptr{Cvoid},
    clockUpdate::Ptr{Cvoid},
    lockPreemption::Ptr{Cvoid},
    unlockPreemption::Ptr{Cvoid})
    @assert false "Not tested!"
    compAddr = ccall(cfunc,
        Ptr{Cvoid},
        (Cstring, Cstring, Cstring,
        Cint, Cint, Cint, Cint, Ptr{fmi3ValueReference}, Csize_t, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
        instanceName, instantiationToken, resourcePath,
        visible, loggingOn, eventModeUsed, earlyReturnAllowed, requiredIntermediateVariables,
        nRequiredIntermediateVariables, instanceEnvironment, logMessage, clockUpdate, lockPreemption, unlockPreemption)

    compAddr
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.1. Super State: FMU State Setable

Disposes the given instance, unloads the loaded model, and frees all the allocated memory and other resources that have been allocated by the functions of the FMU interface. If a NULL pointer is provided for argument instance, the function call is ignored (does not have an effect).
"""
function fmi3FreeInstance!(c::fmi3Instance)

    ind = findall(x->x==c, c.fmu.components)
    deleteat!(c.fmu.components, ind)
    ccall(c.fmu.cFreeInstance, Cvoid, (Ptr{Cvoid},), c.compAddr)

    nothing
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.4. Inquire Version Number of Header Files

This function returns fmi3Version of the fmi3Functions.h header file which was used to compile the functions of the FMU. This function call is allowed always and in all interface types.

The standard header file as documented in this specification has version "3.0-beta.2", so this function returns "3.0-beta.2".
"""
function fmi3GetVersion(cfunc::Ptr{Nothing})

    fmi3Version = ccall(cfunc,
                        Cstring,
                        ())

    unsafe_string(fmi3Version)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.1. Super State: FMU State Setable

The function controls debug logging that is output via the logger function callback. If loggingOn = fmi3True, debug logging is enabled, otherwise it is switched off.
"""
function fmi3SetDebugLogging(c::fmi3Instance, logginOn::fmi3Boolean, nCategories::Unsigned, categories::Ptr{Nothing})
    status = ccall(c.fmu.cSetDebugLogging,
                   Cuint,
                   (Ptr{Nothing}, Cint, Csize_t, Ptr{Nothing}),
                   c.compAddr, logginOn, nCategories, categories)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

Informs the FMU to enter Initialization Mode. Before calling this function, all variables with attribute <Datatype initial = "exact" or "approx"> can be set with the “fmi3SetXXX” functions (the ScalarVariable attributes are defined in the Model Description File, see section 2.4.7). Setting other variables is not allowed.
Also sets the simulation start and stop time.
"""
function fmi3EnterInitializationMode(c::fmi3Instance,toleranceDefined::fmi3Boolean,
    tolerance::fmi3Float64,
    startTime::fmi3Float64,
    stopTimeDefined::fmi3Boolean,
    stopTime::fmi3Float64)
    ccall(c.fmu.cEnterInitializationMode,
          Cuint,
          (Ptr{Nothing}, fmi3Boolean, fmi3Float64, fmi3Float64, fmi3Boolean, fmi3Float64),
          c.compAddr, toleranceDefined, tolerance, startTime, stopTimeDefined, stopTime)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.3. State: Initialization Mode

Informs the FMU to exit Initialization Mode.
"""
function fmi3ExitInitializationMode(c::fmi3Instance)
    ccall(c.fmu.cExitInitializationMode,
          Cuint,
          (Ptr{Nothing},),
          c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.4. Super State: Initialized

Informs the FMU that the simulation run is terminated.
"""
function fmi3Terminate(c::fmi3Instance)
    ccall(c.fmu.cTerminate, Cuint, (Ptr{Nothing},), c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.1. Super State: FMU State Setable

Is called by the environment to reset the FMU after a simulation run. The FMU goes into the same state as if fmi3InstantiateXXX would have been called.
"""
function fmi3Reset(c::fmi3Instance)
    ccall(c.fmu.cReset, Cuint, (Ptr{Nothing},), c.compAddr)
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetFloat32!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Float32}, nvalue::Csize_t)
    ccall(c.fmu.cGetFloat32,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float32}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetFloat32(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Float32}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetFloat32,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float32}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetFloat64!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Float64}, nvalue::Csize_t)
    ccall(c.fmu.cGetFloat64,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetFloat64(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Float64}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetFloat64,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetInt8!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int8}, nvalue::Csize_t)
    ccall(c.fmu.cGetInt8,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int8}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetInt8(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int8}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetInt8,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int8}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetUInt8!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt8}, nvalue::Csize_t)
    ccall(c.fmu.cGetUInt8,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt8}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetUInt8(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt8}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetUInt8,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt8}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetInt16!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int16}, nvalue::Csize_t)
    ccall(c.fmu.cGetInt16,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int16}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetInt16(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int16}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetInt16,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int16}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetUInt16!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt16}, nvalue::Csize_t)
    ccall(c.fmu.cGetUInt16,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt16}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetUInt16(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt16}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetUInt16,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt16}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetInt32!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int32}, nvalue::Csize_t)
    ccall(c.fmu.cGetInt32,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int32}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetInt32(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int32}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetInt32,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int32}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetUInt32!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt32}, nvalue::Csize_t)
    ccall(c.fmu.cGetUInt32,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt32}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetUInt32(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt32}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetUInt32,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt32}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetInt64!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int64}, nvalue::Csize_t)
    ccall(c.fmu.cGetInt64,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int64}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetInt64(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Int64}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetInt64,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int64}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

# TODO test, no variable in FMUs
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetUInt64!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt64}, nvalue::Csize_t)
    ccall(c.fmu.cGetUInt64,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt64}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetUInt64(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3UInt64}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetUInt64,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt64}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetBoolean!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Boolean}, nvalue::Csize_t)
    status = ccall(c.fmu.cGetBoolean,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Boolean}, Csize_t),
          c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetBoolean(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Boolean}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetBoolean,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Boolean}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetString!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Vector{Ptr{Cchar}}, nvalue::Csize_t)
    status = ccall(c.fmu.cGetString,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t,  Ptr{Cchar}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetString(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Union{Array{Ptr{Cchar}}, Array{Ptr{UInt8}}}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetString,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{Cchar}, Csize_t),
                c.compAddr, vr, nvr, value, nvalue)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValues - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetBinary!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, valueSizes::Array{Csize_t}, value::Array{fmi3Binary}, nvalue::Csize_t)
    ccall(c.fmu.cGetBinary,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{Csize_t}, Ptr{fmi3Binary}, Csize_t),
                c.compAddr, vr, nvr, valueSizes, value, nvalue)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetBinary(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, valueSizes::Array{Csize_t}, value::Array{fmi3Binary}, nvalue::Csize_t)
    status = ccall(c.fmu.cSetBinary,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{Csize_t}, Ptr{fmi3Binary}, Csize_t),
                c.compAddr, vr, nvr, valueSizes, value, nvalue)
    status
end

# TODO, Clocks not implemented so far thus not tested
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3GetClock!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Clock})
    status = ccall(c.fmu.cGetClock,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t,  Ptr{fmi3Clock}),
                c.compAddr, vr, nvr, value)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.2. Getting and Setting Variable Values

Functions to get and set values of variables idetified by their valueReference.

nValue - is different from nvr if the value reference represents an array and therefore are more values tied to a single value reference.
"""
function fmi3SetClock(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, value::Array{fmi3Clock})
    status = ccall(c.fmu.cSetClock,
                Cuint,
                (Ptr{Nothing},Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Clock}),
                c.compAddr, vr, nvr, value)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

fmi3GetFMUstate makes a copy of the internal FMU state and returns a pointer to this copy
"""
function fmi3GetFMUState!(c::fmi3Instance, FMUstate::Ref{fmi3FMUState})
    status = ccall(c.fmu.cGetFMUState,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3FMUState}),
                c.compAddr, FMUstate)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

fmi3SetFMUstate copies the content of the previously copied FMUstate back and uses it as actual new FMU state.
"""
function fmi3SetFMUState(c::fmi3Instance, FMUstate::fmi3FMUState)
    status = ccall(c.fmu.cSetFMUState,
                Cuint,
                (Ptr{Nothing}, fmi3FMUState),
                c.compAddr, FMUstate)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

fmi3FreeFMUstate frees all memory and other resources allocated with the fmi3GetFMUstate call for this FMUstate.
"""
function fmi3FreeFMUState!(c::fmi3Instance, FMUstate::Ref{fmi3FMUState})
    status = ccall(c.fmu.cFreeFMUState,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3FMUState}),
                c.compAddr, FMUstate)
    status
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

fmi3SerializedFMUstateSize returns the size of the byte vector which is needed to store FMUstate in it.
"""
function fmi3SerializedFMUStateSize!(c::fmi3Instance, FMUstate::fmi3FMUState, size::Ref{Csize_t})
    status = ccall(c.fmu.cSerializedFMUStateSize,
                Cuint,
                (Ptr{Nothing}, Ptr{Cvoid}, Ptr{Csize_t}),
                c.compAddr, FMUstate, size)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

fmi3SerializeFMUstate serializes the data which is referenced by pointer FMUstate and copies this data in to the byte vector serializedState of length size
"""
function fmi3SerializeFMUState!(c::fmi3Instance, FMUstate::fmi3FMUState, serialzedState::Array{fmi3Byte}, size::Csize_t)
    status = ccall(c.fmu.cSerializeFMUState,
                Cuint,
                (Ptr{Nothing}, Ptr{Cvoid}, Ptr{Cchar}, Csize_t),
                c.compAddr, FMUstate, serialzedState, size)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.6.4. Getting and Setting the Complete FMU State

fmi3DeSerializeFMUstate deserializes the byte vector serializedState of length size, constructs a copy of the FMU state and returns FMUstate, the pointer to this copy.
"""
function fmi3DeSerializeFMUState!(c::fmi3Instance, serialzedState::Array{fmi3Byte}, size::Csize_t, FMUstate::Ref{fmi3FMUState})
    status = ccall(c.fmu.cDeSerializeFMUState,
                Cuint,
                (Ptr{Nothing}, Ptr{Cchar}, Csize_t, Ptr{fmi3FMUState}),
                c.compAddr, serialzedState, size, FMUstate)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9. Clocks

fmi3SetIntervalDecimal sets the interval until the next clock tick
"""
# TODO Clocks and dependencies functions
function fmi3SetIntervalDecimal(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, intervals::Array{fmi3Float64})
    @assert false "Not tested"
    status = ccall(c.fmu.cSetIntervalDecimal,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}),
                c.compAddr, vr, nvr, intervals)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9. Clocks

fmi3SetIntervalFraction sets the interval until the next clock tick
Only allowed if the attribute 'supportsFraction' is set.
"""
function fmi3SetIntervalFraction(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, intervalCounters::Array{fmi3UInt64}, resolutions::Array{fmi3UInt64})
    @assert false "Not tested"
    status = ccall(c.fmu.cSetIntervalFraction,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt64}, Ptr{fmi3UInt64}),
                c.compAddr, vr, nvr, intervalCounters, resolutions)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9. Clocks

fmi3GetIntervalDecimal retrieves the interval until the next clock tick.

For input Clocks it is allowed to call this function to query the next activation interval.
For changing aperiodic Clock, this function must be called in every Event Mode where this clock was activated.
For countdown aperiodic Clock, this function must be called in every Event Mode.
Clock intervals are computed in fmi3UpdateDiscreteStates (at the latest), therefore, this function should be called after fmi3UpdateDiscreteStates.
For information about fmi3IntervalQualifiers, call ?fmi3IntervalQualifier
"""
function fmi3GetIntervalDecimal!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, intervals::Array{fmi3Float64}, qualifiers::fmi3IntervalQualifier)
    @assert false "Not tested"
    status = ccall(c.fmu.cGetIntervalDecimal,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}, Ptr{fmi3IntervalQualifier}),
                c.compAddr, vr, nvr, intervals, qualifiers)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9. Clocks

fmi3GetIntervalFraction retrieves the interval until the next clock tick.

For input Clocks it is allowed to call this function to query the next activation interval.
For changing aperiodic Clock, this function must be called in every Event Mode where this clock was activated.
For countdown aperiodic Clock, this function must be called in every Event Mode.
Clock intervals are computed in fmi3UpdateDiscreteStates (at the latest), therefore, this function should be called after fmi3UpdateDiscreteStates.
For information about fmi3IntervalQualifiers, call ?fmi3IntervalQualifier
"""
function fmi3GetIntervalFraction!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, intervalCounters::Array{fmi3UInt64}, resolutions::Array{fmi3UInt64}, qualifiers::fmi3IntervalQualifier)
    @assert false "Not tested"
    status = ccall(c.fmu.cGetIntervalFraction,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt64}, Ptr{fmi3UInt64}, Ptr{fmi3IntervalQualifier}),
                c.compAddr, vr, nvr, intervalCounters, resolutions, qualifiers)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9. Clocks

fmi3GetShiftDecimal retrieves the delay to the first Clock tick from the FMU.
"""
function fmi3GetShiftDecimal!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, shifts::Array{fmi3Float64})
    @assert false "Not tested"
    status = ccall(c.fmu.cGetShiftDecimal,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}),
                c.compAddr, vr, nvr, shifts)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.9. Clocks

fmi3GetShiftFraction retrieves the delay to the first Clock tick from the FMU.
"""
function fmi3GetShiftFraction!(c::fmi3Instance, vr::Array{fmi3ValueReference}, nvr::Csize_t, shiftCounters::Array{fmi3UInt64}, resolutions::Array{fmi3UInt64})
    @assert false "Not tested"
    status = ccall(c.fmu.cGetShiftFraction,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3UInt64}, Ptr{fmi3UInt64}),
                c.compAddr, vr, nvr, shiftCounters, resolutions)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 5.2.2. State: Clock Activation Mode

During Clock Activation Mode (see 5.2.2.) after fmi3ActivateModelPartition has been called for a calculated, tunable or changing Clock the FMU provides the information on when the Clock will tick again, i.e. when the corresponding model partition has to be scheduled the next time.

Each fmi3ActivateModelPartition call is associated with the computation of an exposed model partition of the FMU and therefore to an input Clock.
"""
function fmi3ActivateModelPartition(c::fmi3Instance, vr::fmi3ValueReference, activationTime::Array{fmi3Float64})
    @assert false "Not tested"
    status = ccall(c.fmu.cActivateModelPartition,
                Cuint,
                (Ptr{Nothing}, fmi3ValueReference, Ptr{fmi3Float64}),
                c.compAddr, vr, activationTime)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.10. Dependencies of Variables

The number of dependencies of a given variable, which may change if structural parameters are changed, can be retrieved by calling fmi3GetNumberOfVariableDependencies.

This information can only be retrieved if the 'providesPerElementDependencies' tag in the ModelDescription is set.
"""
# TODO not tested
function fmi3GetNumberOfVariableDependencies!(c::fmi3Instance, vr::fmi3ValueReference, nvr::Ref{Csize_t})
    @assert false "Not tested"
    status = ccall(c.fmu.cGetNumberOfVariableDependencies,
                Cuint,
                (Ptr{Nothing}, fmi3ValueReference, Ptr{Csize_t}),
                c.compAddr, vr, nvr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.10. Dependencies of Variables

The actual dependencies (of type dependenciesKind) can be retrieved by calling the function fmi3GetVariableDependencies:

dependent - specifies the valueReference of the variable for which the dependencies should be returned.

nDependencies - specifies the number of dependencies that the calling environment allocated space for in the result buffers, and should correspond to value obtained by calling fmi3GetNumberOfVariableDependencies.

elementIndicesOfDependent - must point to a buffer of size_t values of size nDependencies allocated by the calling environment.
It is filled in by this function with the element index of the dependent variable that dependency information is provided for. The element indices start with 1. Using the element index 0 means all elements of the variable. (Note: If an array has more than one dimension the indices are serialized in the same order as defined for values in Section 2.2.6.1.)

independents -  must point to a buffer of fmi3ValueReference values of size nDependencies allocated by the calling environment.
It is filled in by this function with the value reference of the independent variable that this dependency entry is dependent upon.

elementIndicesIndependents - must point to a buffer of size_t values of size nDependencies allocated by the calling environment.
It is filled in by this function with the element index of the independent variable that this dependency entry is dependent upon. The element indices start with 1. Using the element index 0 means all elements of the variable. (Note: If an array has more than one dimension the indices are serialized in the same order as defined for values in Section 2.2.6.1.)

dependencyKinds - must point to a buffer of dependenciesKind values of size nDependencies allocated by the calling environment.
It is filled in by this function with the enumeration value describing the dependency of this dependency entry.
For more information about dependenciesKinds, call ?fmi3DependencyKind

If this function is called before the fmi3ExitInitializationMode call, it returns the initial dependencies. If this function is called after the fmi3ExitInitializationMode call, it returns the runtime dependencies.
The retrieved dependency information of one variable becomes invalid as soon as a structural parameter linked to the variable or to any of its depending variables are set. As a consequence, if you change structural parameters affecting B or A, the dependency of B becomes invalid. The dependency information must change only if structural parameters are changed.

This information can only be retrieved if the 'providesPerElementDependencies' tag in the ModelDescription is set.
"""
function fmi3GetVariableDependencies!(c::fmi3Instance, vr::fmi3ValueReference, elementIndiceOfDependents::Array{Csize_t}, independents::Array{fmi3ValueReference},  elementIndiceOfInpendents::Array{Csize_t}, dependencyKind::Array{fmi3DependencyKind}, ndependencies::Csize_t)
    @assert false "Not tested"
    status = ccall(c.fmu.cGetVariableDependencies,
                Cuint,
                (Ptr{Nothing}, fmi3ValueReference, Ptr{Csize_t}, Ptr{fmi3ValueReference}, Ptr{Csize_t}, Ptr{fmi3DependencyKind}, Csize_t),
                c.compAddr, vr, elementIndiceOfDependents, independents, elementIndiceOfInpendents, dependencyKind, ndependencies)
end

# TODO not tested
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

This function computes the directional derivatives v_{sensitivity} = J ⋅ v_{seed} of an FMU.

unknowns - contains value references to the unknowns.

nUnknowns - contains the length of argument unknowns.

knowns - contains value references of the knowns.

nKnowns - contains the length of argument knowns.

seed - contains the components of the seed vector.

nSeed - contains the length of seed.

sensitivity - contains the components of the sensitivity vector.

nSensitivity - contains the length of sensitivity.

This function can only be called if the 'ProvidesDirectionalDerivatives' tag in the ModelDescription is set.
"""
function fmi3GetDirectionalDerivative!(c::fmi3Instance,
                                       unknowns::Array{fmi3ValueReference},
                                       nUnknowns::Csize_t,
                                       knowns::Array{fmi3ValueReference},
                                       nKnowns::Csize_t,
                                       seed::Array{fmi3Float64},
                                       nSeed::Csize_t,
                                       sensitivity::Array{fmi3Float64},
                                       nSensitivity::Csize_t)
    @assert fmi3ProvidesDirectionalDerivative(c.fmu) ["fmi3GetDirectionalDerivative!(...): This FMU does not support build-in directional derivatives!"]
    ccall(c.fmu.cGetDirectionalDerivative,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}, Csize_t, Ptr{fmi3Float64}, Csize_t),
          c.compAddr, unknowns, nUnknowns, knowns, nKnowns, seed, nSeed, sensitivity, nSensitivity)

end

# TODO not tested
"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.11. Getting Partial Derivatives

This function computes the adjoint derivatives v^T_{sensitivity}= v^T_{seed} ⋅ J of an FMU.

unknowns - contains value references to the unknowns.

nUnknowns - contains the length of argument unknowns.

knowns - contains value references of the knowns.

nKnowns - contains the length of argument knowns.

seed - contains the components of the seed vector.

nSeed - contains the length of seed.

sensitivity - contains the components of the sensitivity vector.

nSensitivity - contains the length of sensitivity.

This function can only be called if the 'ProvidesAdjointDerivatives' tag in the ModelDescription is set.
"""
function fmi3GetAdjointDerivative!(c::fmi3Instance,
                                       unknowns::Array{fmi3ValueReference},
                                       nUnknowns::Csize_t,
                                       knowns::Array{fmi3ValueReference},
                                       nKnowns::Csize_t,
                                       seed::Array{fmi3Float64},
                                       nSeed::Csize_t,
                                       sensitivity::Array{fmi3Float64},
                                       nSensitivity::Csize_t)
    @assert fmi3ProvidesAdjointDerivatives(c.fmu) ["fmi3GetAdjointDerivative!(...): This FMU does not support build-in adjoint derivatives!"]
    ccall(c.fmu.cGetAdjointDerivative,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Float64}, Csize_t, Ptr{fmi3Float64}, Csize_t),
          c.compAddr, unknowns, nUnknowns, knowns, nKnowns, seed, nSeed, sensitivity, nSensitivity)

end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.2.12. Getting Derivatives of Continuous Outputs

Retrieves the n-th derivative of output values.

valueReferences - is a vector of value references that define the variables whose derivatives shall be retrieved. If multiple derivatives of a variable shall be retrieved, list the value reference multiple times.

nValueReferences - is the dimension of the arguments valueReferences and orders.

orders - contains the orders of the respective derivative (1 means the first derivative, 2 means the second derivative, …, 0 is not allowed).
If multiple derivatives of a variable shall be retrieved, provide a list of them in the orders array, corresponding to a multiply occurring value reference in the valueReferences array.
The highest order of derivatives retrievable can be determined by the 'maxOutputDerivativeOrder' tag in the ModelDescription.

values - is a vector with the values of the derivatives. The order of the values elements is derived from a twofold serialization: the outer level corresponds to the combination of a value reference (e.g., valueReferences[k]) and order (e.g., orders[k]), and the inner level to the serialization of variables as defined in Section 2.2.6.1. The inner level does not exist for scalar variables.

nValues - is the size of the argument values. nValues only equals nValueReferences if all corresponding output variables are scalar variables.
"""
function fmi3GetOutputDerivatives!(c::fmi3Instance,  vr::Array{fmi3ValueReference}, nValueReferences::Csize_t, order::Array{fmi3Int32}, values::Array{fmi3Float64}, nValues::Csize_t)
    status = ccall(c.fmu.cGetOutputDerivatives,
                Cuint,
                (Ptr{Nothing}, Ptr{fmi3ValueReference}, Csize_t, Ptr{fmi3Int32}, Ptr{fmi3Float64}, Csize_t),
                c.compAddr, vr, nValueReferences, order, values, nValues)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

If the importer needs to change structural parameters, it must move the FMU into Configuration Mode using fmi3EnterConfigurationMode.
"""
function fmi3EnterConfigurationMode(c::fmi3Instance)
    ccall(c.fmu.cEnterConfigurationMode,
            Cuint,
            (Ptr{Nothing},),
            c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.6. State: Configuration Mode

Exits the Configuration Mode and returns to state Instantiated.
"""
function fmi3ExitConfigurationMode(c::fmi3Instance)
    ccall(c.fmu.cExitConfigurationMode,
          Cuint,
          (Ptr{Nothing},),
          c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

This function returns the number of continuous states.
This function can only be called in Model Exchange.

fmi3GetNumberOfContinuousStates must be called after a structural parameter is changed. As long as no structural parameters changed, the number of states is given in the modelDescription.xml, alleviating the need to call this function.
"""
function fmi3GetNumberOfContinuousStates!(c::fmi3Instance, nContinuousStates::Ref{Csize_t})
    ccall(c.fmu.cGetNumberOfContinuousStates,
            Cuint,
            (Ptr{Nothing}, Ptr{Csize_t}),
            c.compAddr, nContinuousStates)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.2. State: Instantiated

This function returns the number of event indicators.
This function can only be called in Model Exchange.

fmi3GetNumberOfEventIndicators must be called after a structural parameter is changed. As long as no structural parameters changed, the number of states is given in the modelDescription.xml, alleviating the need to call this function.
"""
function fmi3GetNumberOfEventIndicators!(c::fmi3Instance, nEventIndicators::Ref{Csize_t})
    ccall(c.fmu.cGetNumberOfEventIndicators,
            Cuint,
            (Ptr{Nothing}, Ptr{Csize_t}),
            c.compAddr, nEventIndicators)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.3. State: Initialization Mode

Return the states at the current time instant.

This function must be called if fmi3UpdateDiscreteStates returned with valuesOfContinuousStatesChanged == fmi3True. Not allowed in Co-Simulation and Scheduled Execution.
"""
function fmi3GetContinuousStates!(c::fmi3Instance, nominals::Array{fmi3Float64}, nContinuousStates::Csize_t)
    ccall(c.fmu.cGetContinuousStates,
            Cuint,
            (Ptr{Nothing}, Ptr{fmi3Float64}, Csize_t),
            c.compAddr, nominals, nContinuousStates)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.3. State: Initialization Mode

Return the nominal values of the continuous states.

If fmi3UpdateDiscreteStates returned with nominalsOfContinuousStatesChanged == fmi3True, then at least one nominal value of the states has changed and can be inquired with fmi3GetNominalsOfContinuousStates.
Not allowed in Co-Simulation and Scheduled Execution.
"""
function fmi3GetNominalsOfContinuousStates!(c::fmi3Instance, x_nominal::Array{fmi3Float64}, nx::Csize_t)
    status = ccall(c.fmu.cGetNominalsOfContinuousStates,
                    Cuint,
                    (Ptr{Nothing}, Ptr{fmi3Float64}, Csize_t),
                    c.compAddr, x_nominal, nx)
end

# TODO not testable not supported by FMU
"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.3. State: Initialization Mode

This function is called to trigger the evaluation of fdisc to compute the current values of discrete states from previous values.
The FMU signals the support of fmi3EvaluateDiscreteStates via the capability flag providesEvaluateDiscreteStates.
"""
function fmi3EvaluateDiscreteStates(c::fmi3Instance)
    ccall(c.fmu.cEvaluateDiscreteStates,
            Cuint,
            (Ptr{Nothing},),
            c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.5. State: Event Mode

This function is called to signal a converged solution at the current super-dense time instant. fmi3UpdateDiscreteStates must be called at least once per super-dense time instant.
"""
function fmi3UpdateDiscreteStates(c::fmi3Instance, discreteStatesNeedUpdate::Ref{fmi3Boolean}, terminateSimulation::Ref{fmi3Boolean},
                                    nominalsOfContinuousStatesChanged::Ref{fmi3Boolean}, valuesOfContinuousStatesChanged::Ref{fmi3Boolean},
                                    nextEventTimeDefined::Ref{fmi3Boolean}, nextEventTime::Ref{fmi3Float64})
    ccall(c.fmu.cUpdateDiscreteStates,
            Cuint,
            (Ptr{Nothing}, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}, Ptr{fmi3Float64}),
            c.compAddr, discreteStatesNeedUpdate, terminateSimulation, nominalsOfContinuousStatesChanged, valuesOfContinuousStatesChanged, nextEventTimeDefined, nextEventTime)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.5. State: Event Mode

The model enters Continuous-Time Mode and all discrete-time equations become inactive and all relations are “frozen”.
This function has to be called when changing from Event Mode (after the global event iteration in Event Mode over all involved FMUs and other models has converged) into Continuous-Time Mode.
"""
function fmi3EnterContinuousTimeMode(c::fmi3Instance)
    ccall(c.fmu.cEnterContinuousTimeMode,
          Cuint,
          (Ptr{Nothing},),
          c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 2.3.5. State: Event Mode

This function must be called to change from Event Mode into Step Mode in Co-Simulation (see 4.2.).
"""
function fmi3EnterStepMode(c::fmi3Instance)
    ccall(c.fmu.cEnterStepMode,
          Cuint,
          (Ptr{Nothing},),
          c.compAddr)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

Set a new time instant and re-initialize caching of variables that depend on time, provided the newly provided time value is different to the previously set time value (variables that depend solely on constants or parameters need not to be newly computed in the sequel, but the previously computed values can be reused).
"""
function fmi3SetTime(c::fmi3Instance, time::fmi3Float64)
    c.t = time
    ccall(c.fmu.cSetTime,
          Cuint,
          (Ptr{Nothing}, fmi3Float64),
          c.compAddr, time)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

Set a new (continuous) state vector and re-initialize caching of variables that depend on the states. Argument nx is the length of vector x and is provided for checking purposes
"""
function fmi3SetContinuousStates(c::fmi3Instance,
                                 x::Array{fmi3Float64},
                                 nx::Csize_t)
    ccall(c.fmu.cSetContinuousStates,
         Cuint,
         (Ptr{Nothing}, Ptr{fmi3Float64}, Csize_t),
         c.compAddr, x, nx)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

Compute first-oder state derivatives at the current time instant and for the current states.
"""
function fmi3GetContinuousStateDerivatives(c::fmi3Instance,
                            derivatives::Array{fmi3Float64},
                            nx::Csize_t)
    ccall(c.fmu.cGetContinuousStateDerivatives,
          Cuint,
          (Ptr{Nothing}, Ptr{fmi3Float64}, Csize_t),
          c.compAddr, derivatives, nx)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

Compute event indicators at the current time instant and for the current states. EventIndicators signal Events by their sign change.
"""
function fmi3GetEventIndicators!(c::fmi3Instance, eventIndicators::Array{fmi3Float64}, ni::Csize_t)
    status = ccall(c.fmu.cGetEventIndicators,
                    Cuint,
                    (Ptr{Nothing}, Ptr{fmi3Float64}, Csize_t),
                    c.compAddr, eventIndicators, ni)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

This function must be called by the environment after every completed step of the integrator provided the capability flag needsCompletedIntegratorStep = true.
If enterEventMode == fmi3True, the event mode must be entered
If terminateSimulation == fmi3True, the simulation shall be terminated
"""
function fmi3CompletedIntegratorStep!(c::fmi3Instance,
                                      noSetFMUStatePriorToCurrentPoint::fmi3Boolean,
                                      enterEventMode::Ref{fmi3Boolean},
                                      terminateSimulation::Ref{fmi3Boolean})
    ccall(c.fmu.cCompletedIntegratorStep,
          Cuint,
          (Ptr{Nothing}, fmi3Boolean, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}),
          c.compAddr, noSetFMUStatePriorToCurrentPoint, enterEventMode, terminateSimulation)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.2.1. State: Continuous-Time Mode

The model enters Event Mode from the Continuous-Time Mode in ModelExchange oder Step Mode in CoSimulation and discrete-time equations may become active (and relations are not “frozen”).
"""
function fmi3EnterEventMode(c::fmi3Instance, stepEvent::fmi3Boolean, stateEvent::fmi3Boolean, rootsFound::Array{fmi3Int32}, nEventIndicators::Csize_t, timeEvent::fmi3Boolean)
    ccall(c.fmu.cEnterEventMode,
          Cuint,
          (Ptr{Nothing},fmi3Boolean, fmi3Boolean, Ptr{fmi3Int32}, Csize_t, fmi3Boolean),
          c.compAddr, stepEvent, stateEvent, rootsFound, nEventIndicators, timeEvent)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 4.2.1. State: Step Mode

The computation of a time step is started.
"""
function fmi3DoStep!(c::fmi3Instance, currentCommunicationPoint::fmi3Float64, communicationStepSize::fmi3Float64, noSetFMUStatePriorToCurrentPoint::fmi3Boolean,
                    eventEncountered::Ref{fmi3Boolean}, terminateSimulation::Ref{fmi3Boolean}, earlyReturn::Ref{fmi3Boolean}, lastSuccessfulTime::Ref{fmi3Float64})
    @assert c.fmu.cDoStep != C_NULL ["fmi3DoStep(...): This FMU does not support fmi3DoStep, probably it's a ME-FMU with no CS-support?"]

    ccall(c.fmu.cDoStep, Cuint,
          (Ptr{Nothing}, fmi3Float64, fmi3Float64, fmi3Boolean, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}, Ptr{fmi3Boolean}, Ptr{fmi3Float64}),
          c.compAddr, currentCommunicationPoint, communicationStepSize, noSetFMUStatePriorToCurrentPoint, eventEncountered, terminateSimulation, earlyReturn, lastSuccessfulTime)
end


text:
abstract type AbstractExcitation end

λ₀(excitation::AbstractExcitation) = c₀ / frequency(excitation)
ω₀(excitation::AbstractExcitation) = 2π * frequency(excitation)
k₀(excitation::AbstractExcitation) = 2π / λ₀(excitation)

@with_kw struct Harmonic <: AbstractExcitation
    freq
    amp = 1.0
end
frequency(signal::Harmonic) = signal.freq
function (signal::Harmonic)(t)
    @unpack freq, amp = signal
    ω₀ = 2π * freq
    amp * exp(-im * (ω₀ * t))
end
@with_kw struct Gaussian <: AbstractExcitation
    freq
    df
    t0 = 2 * inv(df)
    amp = 1.0
end
frequency(signal::Gaussian) = signal.freq
function (signal::Gaussian)(t)
    @unpack freq, df, t0, amp = signal
    ω₀ = 2π * freq
    σₜ = 1 / (2π * df)
    amp * exp(-0.5((t - t0) * inv(σₜ))^2) * exp(-im * (ω₀ * t))
end

ext: jl


text:
"""
    Plastic(E, ν, σ_y, H, r, κ_∞, α_∞)

Plasticity with von Mises yield surface and mixed non-linear kinematic + non-linear isotropic hardening. Both hardening laws are of saturation type.
# Arguments
- `E::Float64`: Young's modulus
- `ν::Float64`: Poisson's ratio
- `σ_y:` yield limit
- `H:` hardening modulus
- `r:` coupling parameter between isotropic and kinematic hardening
- `κ_∞:` saturation stress for isotropic hardening
- `α_∞:` saturation stress for kinematic hardening
"""
struct Plastic <: AbstractMaterial
    E::Float64 # Young's modulus
    ν::Float64 # Poisson's ratio
    σ_y::Float64 # yield limit
    H::Float64 # hardening modulus
    r::Float64 # coupling parameter between isotropic and kinmatic hardening
    κ_∞::Float64 # saturation stress for isotropic hardening
    α_∞::Float64 # saturation stress for kinematic hardening

    # precomputed elastic tangent tensor and hardening moduli
    Eᵉ::SymmetricTensor{4,3,Float64,36}
    H_κ::Float64
    H_α::Float64

    function Plastic(E, ν, σ_y, H, r, κ_∞, α_∞)
        Eᵉ = elastic_tangent_3D(E, ν)
        H_κ = -r*H
        H_α = -2/3*(1-r)*H

        return new(E, ν, σ_y, H, r, κ_∞, α_∞, Eᵉ, H_κ, H_α)
    end
end

# keyword argument constructor
Plastic(;E, ν, σ_y, H, r, κ_∞, α_∞) = Plastic(E, ν, σ_y, H, r, κ_∞, α_∞)

struct PlasticState{dim,T,M} <: AbstractMaterialState
    εᵖ::SymmetricTensor{2,dim,T,M}
    κ::T
    α::SymmetricTensor{2,dim,T,M}
    μ::T
end

Base.zero(::Type{PlasticState{dim,T,M}}) where {dim,T,M} = PlasticState(zero(SymmetricTensor{2,dim,T,M}), zero(T), zero(SymmetricTensor{2,dim,T,M}), zero(T))
initial_material_state(::Plastic) = zero(PlasticState{3,Float64,6})

function get_cache(m::Plastic)
    state = initial_material_state(m)
    # it doesn't actually matter which state and strain step we use here,
    # f is overwritten in the constitutive driver before being used.
    f(r_vector, x_vector) = vector_residual!(((x)->MaterialModels.residuals(x, m, state, zero(SymmetricTensor{2,3}))), r_vector, x_vector, m)
    v_cache = Vector{Float64}(undef, get_n_scalar_equations(m))
    cache = NLsolve.OnceDifferentiable(f, v_cache, v_cache; autodiff = :forward)
    return cache
end

# constitutive driver operates in 3D, so these can always be 3D
# TODO: Should Residuals have a Type Parameter N for the number of scalar equations?
struct ResidualsPlastic{T}
    σ::SymmetricTensor{2,3,T,6}
    κ::T
    α::SymmetricTensor{2,3,T,6}
    μ::T
end

# doesn't allow other dimesion or unsymmetric tensors, so this is a constant
get_n_scalar_equations(::Plastic) = 14

# Residuals(::Plastic) = ResidualsPlastic(zero(SymmetricTensor{2,3}), zero(Float64), zero(SymmetricTensor{2,3}), zero(Float64))
Tensors.get_base(::Type{Plastic}) = ResidualsPlastic # needed for frommandel

# specified version, fewer allocations, faster than generic fallback
function Tensors.tomandel!(v::Vector{T}, r::ResidualsPlastic{T}) where T
    M=6
    # TODO check vector length
    tomandel!(view(v, 1:M), r.σ)
    v[M+1] = r.κ
    tomandel!(view(v, M+2:2M+1), r.α)
    v[2M+2] = r.μ
    return v
end

function Tensors.frommandel(::Type{ResidualsPlastic}, v::Vector{T}) where T
    σ = frommandel(SymmetricTensor{2,3}, view(v, 1:6))
    κ = v[7]
    α = frommandel(SymmetricTensor{2,3}, view(v, 8:13))
    μ = v[14]
    return ResidualsPlastic{T}(σ, κ, α, μ)
end

"""
    material_response(m::Plastic, ε::SymmetricTensor{2,3,T,6}, state::PlasticState{3}; <keyword arguments>)

Return the stress tensor, stress tangent and the new `MaterialState` for the given strain ε and previous material state `state`.

Plastic free energy:
```math
\\Psi^\\text{p} = \\frac{1}{2} \\, r \\, H \\, k^2
+ \\frac{1}{2}\\left( 1-r \\right) \\, H \\, \\left[ \\sqrt{\\frac{2}{3}} \\left| \\text{dev} \\left(\\mathbf{a} \\right) \\right| \\right]^2
```
Von Mises yield function:
```math
\\Phi = \\sqrt{\\frac{3}{2}} \\left| \\text{dev} \\left( \\boldsymbol{\\sigma} - \\boldsymbol{\\alpha} \\right) \\right| - \\sigma_y - \\kappa
```
An associative flow rule and non-associative hardening rules are used. The evolution equations for the hardening variables are:
```math
\\begin{aligned}
\\dot{k} &= -\\lambda \\left( 1 - \\frac{\\kappa}{\\kappa_{\\infty}} \\right) \\\\
\\dot{\\mathbf{a}} &= -\\lambda \\left( \\frac{\\partial\\Phi}{\\partial\\boldsymbol{\\sigma}}
+ \\frac{3}{2\\alpha_{\\infty}} \\, \\text{dev} \\left( \\boldsymbol{\\alpha} \\right) \\right) \\,.
\\end{aligned}
```
# Keyword arguments
- `cache`: Cache for the iterative solver, used by NLsolve.jl. It is strongly recommended to pre-allocate the cache for repeated calls to `material_response`. See [`get_cache`](@ref).
- `options::Dict{Symbol, Any}`: Solver options for the non-linear solver. Under the key `:nlsolve_params` keyword arguments for `nlsolve` can be handed over.
See [NLsolve documentation](https://github.com/JuliaNLSolvers/NLsolve.jl#common-options). By default the Newton solver will be used.
"""
function material_response(m::Plastic, ε::SymmetricTensor{2,3,T,6}, state::PlasticState{3},
    Δt=nothing; cache=get_cache(m), options::Dict{Symbol, Any} = Dict{Symbol, Any}()) where T

    σ_trial = m.Eᵉ ⊡ (ε - state.εᵖ)
    Φ = sqrt(3/2)*norm(dev(σ_trial-state.α)) - m.σ_y - state.κ

    if Φ <= 0
        return σ_trial, m.Eᵉ, state
    else
        # set the current residual function that depends only on the variables
        # cache.f = (r_vector, x_vector) -> vector_residual!(((x)->residuals(x,m,state,Δε)), r_vector, x_vector, m)
        f(r_vector, x_vector) = vector_residual!(((x)->residuals(x,m,state,ε)), r_vector, x_vector, m)
        update_cache!(cache, f)
        # initial guess
        x0 = ResidualsPlastic(σ_trial, state.κ, state.α, state.μ)
        # convert initial guess to vector
        tomandel!(cache.x_f, x0)
        # solve for variables x
        nlsolve_options = get(options, :nlsolve_params, Dict{Symbol, Any}(:method=>:newton))
        haskey(nlsolve_options, :method) || merge!(nlsolve_options, Dict{Symbol, Any}(:method=>:newton)) # set newton if the user did not supply another method
        result = NLsolve.nlsolve(cache, cache.x_f; nlsolve_options...)
        if result.f_converged
            x = frommandel(ResidualsPlastic, result.zero::Vector{T})
            εᵖ = state.εᵖ + x.μ*(3/(2*sqrt(3/2)*norm(dev(x.σ-x.α)))*dev(x.σ-x.α))
            dRdx = cache.DF
            inv_J_σσ = frommandel(SymmetricTensor{4,3}, inv(dRdx))
            ∂σ∂ε = inv_J_σσ ⊡ m.Eᵉ
            return x.σ, ∂σ∂ε, PlasticState(εᵖ, x.κ, x.α, x.μ)
        else
            error("Material model not converged. Could not find material state.")
        end
    end
end

function residuals(vars::ResidualsPlastic, m::Plastic, material_state::PlasticState{3}, ε)
    σ = vars.σ; κ = vars.κ; α = vars.α; μ = vars.μ

    σ_trial = m.Eᵉ ⊡ (ε - material_state.εᵖ)
    ν = 3/(2*sqrt(3/2)*norm(dev(σ-α)))*dev(σ-α)
    Rσ = σ - σ_trial + m.Eᵉ ⊡ (μ*ν) # R_σ(εᵖ, α, μ)
    Rκ = κ - material_state.κ - μ*m.H_κ*(-1 +κ/m.κ_∞) # R_κ(κ, μ)
    Rα = α - material_state.α - μ*m.H_α * (-ν + 3/(2*m.α_∞)*dev(α)) # R_α(σ, α, μ)
    Rμ = sqrt(3/2)*norm(dev(σ-α)) - m.σ_y - κ
    return ResidualsPlastic(Rσ, Rκ, Rα, Rμ)
end


text:
using Glob
using Literate
using Random

config = Dict(
    :binder_root_url =>
        "https://mybinder.org/v2/gh/maxmouchet/HMMBase.jl/master?filepath=",
    :repo_root_url => "https://github.com/maxmouchet/HMMBase.jl/blob/master",
)

function literate_documenter(inputdir, outputdir)
    Random.seed!(2019)
    map(glob("*.jl", inputdir)) do file
        Literate.markdown(file, outputdir, documenter = true; config...)
    end
end

function literate_notebooks(inputdir, outputdir)
    Random.seed!(2019)
    map(glob("*.jl", inputdir)) do file
        Literate.notebook(file, outputdir; config...)
    end
end

literate_documenter("examples/", "docs/src/examples/")
literate_notebooks("examples/", "examples/")


text:
using Test
using ClimateMachine.SurfaceFluxes.UniversalFunctions
using ClimateMachine.SurfaceFluxes.UniversalFunctions: b_m, a_m, a_h, b_h, Pr_0
using CLIMAParameters
using CLIMAParameters.SurfaceFluxes.UniversalFunctions
using CLIMAParameters.Planet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

# TODO: Right now, we test these functions for
# type stability and correctness in the asymptotic
# limit. We may want to extend correctness tests.

@testset "UniversalFunctions" begin
    @testset "Type stability" begin
        FT = Float32
        ζ = FT(-2):FT(0.01):FT(200)
        for L in (-FT(10), FT(10))
            args = (param_set, L)
            for uf in (Gryanik(args...), Grachev(args...), Businger(args...))
                for transport in (MomentumTransport(), HeatTransport())
                    ϕ = phi.(uf, ζ, transport)
                    @test eltype(ϕ) == FT
                    ψ = psi.(uf, ζ, transport)
                    @test eltype(ψ) == FT
                end
            end
        end
    end
    @testset "Conversions" begin
        FT = Float32
        ζ = FT(10)
        L = FT(10)
        args = (param_set, L)

        uf = Gryanik(args...)
        @test Businger(uf) isa Businger
        @test Grachev(uf) isa Grachev

        uf = Grachev(args...)
        @test Businger(uf) isa Businger
        @test Gryanik(uf) isa Gryanik

        uf = Businger(args...)
        @test Grachev(uf) isa Grachev
        @test Gryanik(uf) isa Gryanik
    end
    @testset "Asymptotic range" begin
        FT = Float32

        ϕ_h_ζ∞(uf::Grachev) = 1 + FT(b_h(uf))
        ϕ_m_ζ∞(uf::Grachev, ζ) = FT(a_m(uf)) / FT(b_m(uf)) * ζ^FT(1 / 3)

        ϕ_h_ζ∞(uf::Gryanik) = FT(Pr_0(uf)) * (1 + FT(a_h(uf) / b_h(uf)))
        ϕ_m_ζ∞(uf::Gryanik, ζ) = FT(a_m(uf) / b_m(uf)^FT(2 / 3)) * ζ^FT(1 / 3)

        for L in (-FT(10), FT(10))
            args = (param_set, L)
            for uf in (Grachev(args...), Gryanik(args...))
                for ζ in FT(10) .^ (4, 6, 8, 10)
                    ϕ_h = phi(uf, ζ, HeatTransport())
                    @test isapprox(ϕ_h, ϕ_h_ζ∞(uf))
                end
                for ζ in FT(10) .^ (8, 9, 10)
                    ϕ_m = phi(uf, ζ, MomentumTransport())
                    @test isapprox(ϕ_m, ϕ_m_ζ∞(uf, ζ))
                end
            end
        end

    end

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using ..Config: GDConfig

mutable struct GDStats{T <: Unsigned} <: Statistics
    T::Type{T}
    base::BaseStats
    network::GDNetStats

    hash_round_changes::Vector{Float32}
    deviation_round_changes::Vector{Float32}

    GDStats{T}(config::Configuration, num_weights::Int) where T <: Unsigned = begin
        gdfile_length = ceil(Int, num_weights / config.chunksize)
        # call hash function with fake data to find the hash size
        hash_size = sizeof(config.fingerprint([0x00]))
        lsbsize = sizeof(T) * 8 - config.msbsize
        deviation_size = lsbsize * config.chunksize / 8
        basis_size = config.msbsize * config.chunksize / 8

        return new(
            T,
            BaseStats(
                config.base.num_comm_rounds,
                config.base.fraction_clients,
                config.base.num_total_clients,
                num_weights
            ),
            GDNetStats(gdfile_length, basis_size, deviation_size, hash_size),
            Vector{Float32}(undef, 0),
            Vector{Float32}(undef, 0)
        )
    end
end

function update_stats!(stats::GDStats)
    req_hashes = stats.base.req_data.hashes
    res_hashes = map(gdfile -> gdfile.hashes, stats.base.res_data)
    push!(stats.hash_round_changes, compute_round_changes(req_hashes, res_hashes))

    req_dev = stats.base.req_data.deviations
    res_dev = map(gdfile -> gdfile.deviations, stats.base.res_data)
    push!(stats.deviation_round_changes, compute_round_changes(req_dev, res_dev))
end


text:
using SteadyState
using Test
using QuantumOptics, SparseArrays
using IterativeSolvers

@testset "SteadyState.jl" begin
    b = GenericBasis(15)
    H = SparseOperator(b,sprand(ComplexF64,b.shape[],b.shape[],0.1))
    H = H + dagger(H)
    J = [SparseOperator(b,sprand(ComplexF64,b.shape[],b.shape[],0.1)) for i in 1:3]

    tol = 1e-9

    # Reference
    r = steadystate.eigenvector(H, J; tol=tol)

    # Dense
    r1 = steadystate_bicg(DenseOperator(H), DenseOperator.(J), 4; tol=tol)
    @test tr(r1) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r1.data .- r1.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r1.data .- r.data')) ≈ 0.0 atol=tol

    #Sparse
    #=
    r2 = steadystate_bicg(H, J, 4; tol=tol)
    @test tr(r2) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r2.data .- r2.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r2.data .- r.data')) ≈ 0.0 atol=tol
    =#

    r0 = SparseOperator(b)
    r0.data[1,1] = one(ComplexF64)

    # Dense
    r3 = steadystate_iterative!(deepcopy(DenseOperator(r0)), DenseOperator(H), DenseOperator.(J), :bicgstabl!, 2; tol=tol)
    @test tr(r3) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r3.data .- r3.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r3.data .- r.data')) ≈ 0.0 atol=tol

    #Sparse
    #=
    r4 = steadystate_iterative!(deepcopy(r0), H, J, :bicgstabl!, 2; tol=tol)
    @test tr(r4) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r4.data .- r4.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r4.data .- r.data')) ≈ 0.0 atol=tol
    =#

    # Dense
    r5 = steadystate_iterative!(deepcopy(DenseOperator(r0)), DenseOperator(H), DenseOperator.(J), :gmres!; tol=tol)
    @test tr(r5) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r5.data .- r5.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r5.data .- r.data')) ≈ 0.0 atol=tol

    # Dense
    r6 = steadystate_iterative!(deepcopy(DenseOperator(r0)), DenseOperator(H), DenseOperator.(J), :idrs!; s=4, tol=tol)
    @test tr(r6) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r6.data .- r6.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r6.data .- r.data')) ≈ 0.0 atol=tol

    # Sparse
    r7 = steadystate_iterative!(deepcopy(r0), H, J, :idrs!; s=4, tol=tol)
    @test tr(r7) ≈ one(ComplexF64) atol=tol
    @test sum(abs2.(r7.data .- r7.data')) ≈ 0.0 atol=tol
    @test sum(abs2.(r7.data .- r.data')) ≈ 0.0 atol=tol
end


text:
const TOP_INSTITUTION_FIELDS = [:name, :fisyr, :netassets, :charitytype,
  :categorynamefiltered]

const TOP_INSTITUTION_METHODS = [:rank, :quantile, :cumulative]

const TOP_INSTITUTIONS_YEAR = 2015
const TOP_INSTITUTIONS_QUANTILE = 0.95
const TOP_INSTITUTIONS_CUMULATIVE = 0.95
const TOP_INSTITUTIONS_RANK = 500
const TOP_INSTITUTIONS_PREFIX = "top_"

processmissing(v::AbstractVector) = collect(skipmissing(v))

#writes out the institutions in the subdataframe
function writeTopInstitutions(subdf::SubDataFrame, label::String;
  outputPath::String = TABLE_PATH,
  topInstitutionFields::Vector{Symbol} = TOP_INSTITUTION_FIELDS)::Nothing

  outPath = "$outputPath\\$(label).csv"
  uCSV.write(outPath, subdf[topInstitutionFields])

  return nothing
end

getTopInstitutionsRank(subdf::SubDataFrame,
  topInstitutionsRank::Int)::SubDataFrame = view(subdf, 1:topInstitutionsRank, :)

function getTopInstitutionsQuantile(subdf::SubDataFrame,
  topInstitutionsQuantile::Float64)
  v::Vector{Float64} = subdf[:netassets]
  qFunc::Function = StatsBase.ecdf(v)
  dfquantile::Vector{Float64} = (qFunc).(subdf[:netassets])

  return view(subdf, dfquantile .> topInstitutionsQuantile, :)
end

function getTopInstitutionsCumulative(subdf::SubDataFrame,
  topInstitutionsCumulative::Float64)

  v::Vector{Float64} = subdf[:netassets]
  cumv::Vector{Float64} = similar(v) #will hold the running total
  cumv[1] = v[1]
  for i::Int ∈ 2:length(v)
    cumv[i] = cumv[i-1] + v[i]
  end
  cumv .= cumv ./ cumv[end] #normalize

  return view(subdf, cumv .< topInstitutionsCumulative, :)
end

#gets the top institutions given a method and SORTED subdataframe
function getTopInstitutions(subdf::SubDataFrame, topMethod::Symbol;
  topInstitutionsYear::Int = TOP_INSTITUTIONS_YEAR,
  topInstitutionsQuantile::Float64 = TOP_INSTITUTIONS_QUANTILE,
  topInstitutionsRank::Int = TOP_INSTITUTIONS_RANK,
  topInstitutionsCumulative::Float64 = TOP_INSTITUTIONS_CUMULATIVE,
  topInstitutionsPrefix::String = TOP_INSTITUTIONS_PREFIX)

  if topMethod == :rank #this assumes we want the top x ranked asset holders
    subdf = getTopInstitutionsRank(subdf, topInstitutionsRank)

  elseif topMethod == :quantile #this assumes we want the top  x percentile asset holders
    subdf = getTopInstitutionsQuantile(subdf, topInstitutionsQuantile)

  elseif topMethod == :cumulative #this assumes we want the top holders which account for x% of assets
    subdf = getTopInstitutionsCumulative(subdf, topInstitutionsCumulative)
  elseif topMethod ≠ :full
      @warn "Top institutions method not detected, assuming full file"
  end

  label::String =
    "$(topInstitutionsPrefix)_$(topInstitutionsYear)_$(topMethod)"
  writeTopInstitutions(subdf, label)

  return nothing
end

#takes in an unsorted data object and gets the top institutions by assets
function getTopInstitutions(data::NCCSData, topMethod::Symbol;
    topInstitutionsYear::Int = TOP_INSTITUTIONS_YEAR)
  sort!(data.df, :netassets, rev=true)
  subdf::SubDataFrame =
    view(data.df, (data.df[:fisyr] .== topInstitutionsYear) .& (
    (!ismissing).(data.df[:netassets])), :)

  getTopInstitutions(subdf, topMethod)

  return nothing
end

#convenience method for running multiple methods at once
function getTopInstitutions(data::NCCSData;
  topMethods::Vector{Symbol} = TOP_INSTITUTION_METHODS,
  topInstitutionsYear::Int = TOP_INSTITUTIONS_YEAR,
  sorted::Bool = false)

  (!sorted) && sort!(data.df, [:fisyr, :adjnetassets], (:fisyr, order(:adjnetassets, rev=true)))
  subdf::SubDataFrame =
    view(data.df, (data.df[:fisyr] .== topInstitutionsYear) .& (
    (!ismissing).(data.df[:netassets])),:)

  #Execute the tasks
  for s::Symbol ∈ topMethods
    getTopInstitutions(subdf, s)
  end

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# CWT (continuous wavelet transform)
# cwt(Y::AbstractVector, ::ContWave)

@doc """
     wave = cwt(Y::AbstractArray{T,N}, c::CWT{W, S, WT}, daughters, rfftPlan =
             plan_rfft([1]), fftPlan = plan_fft([1])) where {N, T<:Real,
                                                             S<:Real,
                                                             U<:Number,
                                                             W<:WaveletBoundary,
                                                             WT<:Union{<:Morlet,
                                                                       <:Paul}}

  return the continuous wavelet transform along the first axis with averaging.
  `wave`, is (signalLength)×(nscales+1)×(previous dimensions), of type T of
  Y. averagingLength defines the number of octaves (powers of 2) that are
  replaced by an averaging function. This has form averagingType, which can be
  one of `Father()` or `Dirac()`- in the `Father()` case, it uses the same form
  as for the wavelets, while the `Dirac` uses a constant window. J1 is the
  total number of scales; default (when J1=NaN, or is negative) is just under
  the maximum possible number, i.e. the log base 2 of the length of the signal,
  times the number of wavelets per octave. If you have sampling information,
  you will need to scale wave by δt^(1/2).
  The default assumption is that the sampling rate is 1kHz.

  """
function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N, T}
    @assert typeof(N)<:Integer
    # vectors behave a bit strangely, so we reshape them
    if N==1
        Y= reshape(Y,(length(Y), 1))
    end
    n1 = size(Y, 1);

    _, nScales, _ = getNWavelets(n1, cWav)
    #construct time series to analyze, pad if necessary
    x = reflect(Y, boundaryType(cWav)()) #this function is defined below

    # check if the plans we were given are dummies or not
    x̂, fromPlan = prepSignalAndPlans(x, cWav, fftPlans)
    # If the vector isn't long enough to actually have any other scales, just
    # return the averaging
    if nScales <= 0 || size(daughters,2) == 1
        daughters = daughters[:,1:1]
        nScales = 1
    end

    if isAnalytic(cWav.waveType)
        OutType = ensureComplex(T)
    else
        OutType = T
    end

    wave = zeros(OutType, size(x)..., nScales);  # result array;
    # faster if we put the example index on the outside loop through all scales
    # and compute transform
    if isAnalytic(cWav.waveType)
        if eltype(x) <: Real
            analyticTransformReal!(wave, daughters, x̂, fromPlan, cWav.averagingType)
        else
            analyticTransformComplex!(wave, daughters, x̂, fromPlan, cWav.averagingType)
        end
    else
        otherwiseTransform!(wave, daughters, x̂, fromPlan, cWav.averagingType)
    end
    wave = permutedims(wave, [1, ndims(wave), (2:(ndims(wave)-1))...])
    ax = axes(wave)
    wave = wave[1:n1, ax[2:end]...]
    if N==1
        wave = dropdims(wave, dims=3)
    end

    return wave
end
function ensureComplex(T)
    if T <: Real
        return Complex{T}
    else
        return T
    end
end
# there are 4 cases to deal with
#       input Type | Real | Complex
#       analytic?  |----------------
#             yes  | both | fft
#              no  | rfft | fft
#  Analytic on Real input
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy, N, true}, fftPlans) where {T <: Real, W,S,WaTy,N}
    # analytic wavelets that are being applied on real inputs
    if fftPlans isa Tuple{<:AbstractFFTs.Plan{<:Real}, <:AbstractFFTs.Plan{<:Complex}}
        # they handed us the right kind of thing, so no need to make new ones
        x̂ = fftPlans[1] * x
        fromPlan = fftPlans[2]
    else
        toPlan = plan_rfft(x,1)
        x̂ = toPlan * x
        fromPlan = plan_fft(x,1)
    end
    return x̂, fromPlan
end

#  Non-analytic on Real input
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy, N, false}, fftPlans) where {T <: Real, W,S,WaTy,N}
    # real wavelets that are being applied on real inputs
    if fftPlans isa AbstractFFTs.Plan{<:Real}
        # they handed us the right kind of thing, so no need to make new ones
        x̂ = fftPlans * x
        fromPlan = fftPlans
    else
        fromPlan = plan_rfft(x,1)
        x̂ = fromPlan * x
    end
    return x̂, fromPlan
end
# complex input
function prepSignalAndPlans(x::AbstractArray{T}, cWav, fftPlans) where {T <: Complex, W,S,WaTy,N}
    # real wavelets that are being applied on real inputs
    if fftPlans isa AbstractFFTs.Plan{<:Complex}
        # they handed us the right kind of thing, so no need to make new ones
        x̂ = fftPlans * x
        fromPlan = fftPlans
    else
        fromPlan = plan_fft(x,1)
        x̂ = fromPlan * x
    end
    return x̂, fromPlan
end

# analytic on real data with an averaging function
function analyticTransformReal!(wave, daughters, x̂, fftPlan, averagingType::Union{Father, Dirac})
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    isSourceEven = mod(size(wave,1)+1,2)
    # the averaging function isn't analytic, so we need to do both positive and
    # negative frequencies
    @views tmpWave = x̂ .* daughters[:,1]
    @views wave[(n1+1):end, outer..., 1] = reverse(conj.(tmpWave[2:end-isSourceEven, outer...]), dims=1)
    @views wave[1:n1, outer..., 1] = tmpWave
    @views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # averaging
    for j in 2:size(daughters,2)
        @views wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

# analytic on complex data with an averaging function
function analyticTransformComplex!(wave, daughters, x̂, fftPlan, averagingType::Union{Father, Dirac})
    outer = axes(x̂)[2:end]
    n1 = size(daughters, 1)
    isSourceEven = mod(size(wave,1)+1,2)
    # the averaging function isn't analytic, so we need to do both positive and
    # negative frequencies
    @views positiveFreqs = x̂[1:n1, outer...] .* daughters[:,1]
    @views negativeFreqs = x̂[(n1-isSourceEven+1):end, outer...] .* reverse(conj.(daughters[2:end,1]))
    @views wave[(n1-isSourceEven+1):end, outer..., 1] = negativeFreqs
    @views wave[1:n1, outer..., 1] = positiveFreqs
    @views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # averaging
    for j in 2:size(daughters,2)
        @views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:,j]
        @views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function analyticTransformComplex!(wave, daughters, x̂, fftPlan, averagingType)
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    for j in 1:size(daughters,2)
        @views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:,j]
        @views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

# analytic on real data without an averaging function
function analyticTransformReal!(wave, daughters, x̂, fftPlan, averagingType::NoAve)
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    # the no averaging version
    for j in 1:size(daughters,2)
        wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function otherwiseTransform!(wave::AbstractArray{<:Real}, daughters, x̂, fromPlan, averagingType)
    # real wavelets on real data: that just makes sense
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    for j in 1:size(daughters, 2)
        @views tmp = x̂ .* daughters[:,j]
        @views wave[:, outer..., j] = fromPlan \ tmp  # wavelet transform
    end
end

# if it isn't analytic, the output is complex only if the input is complex
function otherwiseTransform!(wave::AbstractArray{<:Complex}, daughters, x̂, fromPlan, averagingType)
    # applying a real transform to complex data is maybe a bit odd, but you do you
    outer = axes(x̂)[2:end]
    n1 = size(daughters, 1)
    isSourceEven = mod(size(fromPlan, 1) + 1, 2)
    for j in 1:size(daughters, 2)
        @views wave[1:n1, outer..., j] = @views x̂[1:n1,outer...] .* daughters[:,j]
        @views wave[n1-isSourceEven+1:end, outer..., j] =  x̂[n1-isSourceEven+1:end,outer...] .* reverse(conj.(daughters[2:end,j]))
        @views wave[:, outer..., j] = fromPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function reflect(Y, bt)
    n1 = size(Y, 1)
    if typeof(bt) <: ZPBoundary
        base2 = ceil(Int, log2(n1));   # power of 2 nearest to N
        x = cat(Y, zeros(2^(base2)-n1, size(Y)[2:end]...), dims=1)
    elseif typeof(bt) <: SymBoundary
        x = cat(Y, reverse(Y,dims = 1), dims = 1)
    else
        x = Y
    end
    return x
end

function cwt(Y::AbstractArray{T}, c::CWT{W}; varArgs...) where {T<:Number, S<:Real, V<: Real,
                                                                                W<:WaveletBoundary}
    daughters,ω = computeWavelets(size(Y, 1), c; varArgs...)
    return cwt(Y, c, daughters)
end

"""
period,scale, coi = caveats(Y::AbstractArray{T}, c::CWT{W}; J1::S=NaN) where {T<:Real, S<:Real, W<:WaveletBoundary}

returns the period, the scales, and the cone of influence for the given wavelet transform. If you have sampling information, you will need to scale the vector scale appropriately by 1/δt, and the actual transform by δt^(1/2).
"""
function caveats(n1, c::CWT{W}; J1::Int64=-1, dt::S=1/1000, s0::V=NaN) where {S<:Real, W<:WaveletBoundary, V <: Real}
    # don't alter scaling with sampling information if it doesn't exists
    fλ = (4*π) / (c.σ[1] + sqrt(2 + c.σ[1]^2))
    if isnan(dt) || (dt<0)
        dt = 1
    end

    # smallest resolvable scale
    if isnan(s0) || (s0<0)
        s0 = 2 * dt / fλ
    end
    sj = s0 * 2.0.^(collect(0:J1)./c.Q)
    # Fourier equivalent frequencies
    freqs = 1 ./ (fλ .* sj)

    # Determines the cone-of-influence. Note that it is returned as a function
    # of time in Fourier periods. Uses triangualr Bartlett window with
    # non-zero end-points.
    coi = (n1 / 2 .- abs.(collect(0:n1-1) .- (n1 - 1) ./ 2))
    coi = (fλ * dt / sqrt(2)).*coi

    n1 = length(Y);
    # J1 is the total number of elements
    if isnan(J1) || (J1<0)
        J1=floor(Int,(log2(n1))*c.Q);
    end
    #....construct time series to analyze, pad if necessary
    if boundaryType(c) == ZPBoundary
        base2 = round(Int,log(n1)/log(2));   # power of 2 nearest to N
        n = length(Y)+2^(base2+1)-n1
    elseif boundaryType(c) == PerBoundary
        n = length(Y)*2
    end
    ω = [0:floor(Int, n/2); -floor(Int,n/2)+1:-1]*2π
    period = c.fourierFactor*2 .^((0:J1)/c.Q)
    scale = [1E-5; 1:((n1+1)/2-1); reverse((1:(n1/2-1)),dims=1); 1E-5]
    coi = c.coi*scale  # COI [Sec.3g]
    return sj, freqs, period, scale, coi
end
cwt(Y::AbstractArray{T}, w::ContWaveClass; varargs...) where {T<:Number, S<:Real, V<:Real} = cwt(Y,CWT(w); varargs...)
caveats(Y::AbstractArray{T}, w::ContWaveClass; J1::S=NaN) where {T<: Number, S<: Real} = caveats(Y,CWT(w),J1=J1)
cwt(Y::AbstractArray{T}) where T<:Real = cwt(Y,Morlet())
caveats(Y::AbstractArray{T}) where T<:Real = caveats(Y,Morlet())

abstract type InverseType end
struct DualFrames <: InverseType end
struct NaiveDelta <: InverseType end
struct PenroseDelta <: InverseType end

"""
    icwt(res::AbstractArray, cWav::CWT, inverseStyle::InverseType=PenroseDelta())
Compute the inverse wavelet transform using one of three dual frames. The default uses delta functions with weights chosen via a least squares method, the `PenroseDelta()` below. This is chosen as a default because the Morlet wavelets tend to fail catastrophically using the canonical dual frame (the `dualFrames()` type).

    icwt(res::AbstractArray, cWav::CWT, inverseStyle::PenroseDelta)
return the inverse continuous wavelet transform, computed using the simple dual frame ``β_jδ_{ji}``, where ``β_j`` is chosen to solve the least squares problem ``\\|Ŵβ-1\\|_2^2``, where ``Ŵ`` is the Fourier domain representation of the `cWav` wavelets. In both this case and `NaiveDelta()`, the fourier transform of ``δ`` is the constant function, thus this least squares problem.

    icwt(res::AbstractArray, cWav::CWT, inverseStyle::NaiveDelta)
return the inverse continuous wavelet transform, computed using the simple dual frame ``β_jδ_{ji}``, where ``β_j`` is chosen to negate the scale factor ``(^1/_s)^{^1/_p}``. Generally less accurate than choosing the weights using `PenroseDelta`. This is the method discussed in Torrence and Compo.

    icwt(res::AbstractArray, cWav::CWT, inverseStyle::dualFrames)
return the inverse continuous wavelet transform, computed using the canonical dual frame ``\\tilde{\\widehat{ψ}} = \\frac{ψ̂_n(ω)}{∑_n\\|ψ̂_n(ω)\\|^2}``. The algorithm is to compute the cwt again, but using the canonical dual frame; consiquentially, it is the most computationally intensive of the three algorithms, and typically the best behaved. Will be numerically unstable if the high frequencies of all of the wavelets are too small however, and tends to fail spectacularly in this case.

"""
function icwt(res::AbstractArray, cWav::CWT, inverseStyle::PenroseDelta)
    Ŵ = computeWavelets(size(res,1), cWav)[1]
    β = computeDualWeights(Ŵ, cWav)
    testDualCoverage(β, Ŵ)
    compXRecon = sum(res .* β, dims=2)
    imagXRecon = irfft(im*rfft(imag.(compXRecon),1), size(compXRecon,1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
    return imagXRecon + real.(compXRecon)
end

function icwt(res::AbstractArray, cWav::CWT, inverseStyle::NaiveDelta)
    Ŵ = computeWavelets(size(res,1), cWav)[1]
    β = computeNaiveDualWeights(Ŵ, cWav, size(res,1))
    testDualCoverage(β, Ŵ)
    compXRecon = sum(res .* β, dims=2)
    imagXRecon = irfft(im*rfft(imag.(compXRecon),1), size(compXRecon,1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
    return imagXRecon + real.(compXRecon)
end

function icwt(res::AbstractArray, cWav::CWT, inverseStyle::DualFrames)
    Ŵ = computeWavelets(size(res,1), cWav)[1]
    canonDualFrames = explicitConstruction(Ŵ)
    testDualCoverage(canonDualFrames)
    convolved = cwt(res, cWav, canonDualFrames)
    ax = axes(convolved)
    @views xRecon = sum(convolved[:,i,i,ax[4:end]...] for i=1:size(Ŵ,2))
    return xRecon
end

icwt(Y::AbstractArray, w::ContWaveClass; varargs...) where {S<:Real, T<:Real, V<:Real} = icwt(Y,CWT(w), PenroseDelta(); varargs...)
icwt(Y::AbstractArray; varargs...) = icwt(Y,Morlet(), PenroseDelta(); varargs...)

# CWT (continuous wavelet transform directly) TODO: direct if sufficiently small


text:
"""
Based on the hao algorithm

"""
function isinside(polygon, point)
    k = 0
    f = 0
    u1 = 0
    v1 = 0
    u2 = 0
    v2 = 0
    currentP = null
    nextP = null

    x = point[0]
    y = point[1]

    numContours = length(polygon)
    for i = 1:numContours
        ii = 0
        contourLen = length(polygon[i]) - 1
        contour = polygon[i]
        currentP = contour[1]
        v1 = currentP[0] - x
        u1 = currentP[1] - y

        for ii=1:contourLen
              nextP = contour[ii + 1]

              v2 = nextP[1] - y

              if (v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)
                  currentP = nextP
                  v1 = v2
                  u1 = currentP[0] - x
                  continue
              end

              u2 = nextP[0] - p[0]

              if (v2 > 0 && v1 <= 0)
                  f = (u1 * v2) - (u2 * v1)
                  if (f > 0)
                      k = k + 1
                  elseif (f === 0)
                      return 0
                  end
              elseif (v1 > 0 && v2 <= 0)
                  f = (u1 * v2) - (u2 * v1)
                  if (f < 0)
                      k = k + 1
                  elseif (f === 0)
                      return 0
                  end
              elseif (v2 === 0 && v1 < 0)
                  f = (u1 * v2) - (u2 * v1)
                  if (f === 0)
                      return 0
                  end
              elseif (v1 === 0 && v2 < 0)
                  f = u1 * v2 - u2 * v1
                  if (f === 0)
                      return 0
                  end
              elseif (v1 === 0 && v2 === 0)
                  if (u2 <= 0 && u1 >= 0)
                      return 0
                  elseif (u1 <= 0 && u2 >= 0)
                      return 0
              end
              currentP = nextP
              v1 = v2
              u1 = u2
          end
    end

    if k % 2 == 0
        return false
    end
    return true

end


text:
using Documenter
using NewsvendorModel
using Distributions

makedocs(
    sitename = "NewsvendorModel.jl",
    format = Documenter.HTML(),
    modules = [NewsvendorModel],
    doctest = false
)

deploydocs(
    repo = "github.com/frankhuettner/NewsvendorModel.jl.git",
    devbranch = "main",
)


text:
# Test dissipator definition
@test diss(q('x'))(q('z')) ≈ q('x') * q('z') * q('x')' - (1/2)*(q('x')'*q('x')*q('z') + q('z')*q('x')'*q('x'))
@test diss(o('x'))(o('n')) ≈ o('x') * o('n') * o('x')' - (1/2)*(o('x')'*o('x')*o('n') + o('n')*o('x')'*o('x'))
@test diss(o('y'))(o('n')) ≈ o('y') * o('n') * o('y')' - (1/2)*(o('y')'*o('y')*o('n') + o('n')*o('y')'*o('y'))
@test unsuperket(sdiss(o('y')) * superket(o('n'))) ≈ o('y') * o('n') * o('y')' - (1/2)*(o('y')'*o('y')*o('n') + o('n')*o('y')'*o('y'))

# Test innovation definition
@test inn(q('x'))(q('z')) == q('x') * q('z') +  q('z') * q('x')' - tr((q('x') + q('x')')*q('z'))*q('z')
@test inn(o('x'))(o('n')) == o('x') * o('n') +  o('n') * o('x')' - tr((o('x') + o('x')')*o('n'))*o('n')
@test inn(o('y'))(o('n')) == o('y') * o('n') +  o('n') * o('y')' - tr((o('y') + o('y')')*o('n'))*o('n')