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

FMU_NUM_STATES = 0
# FMU_NUM_DERIVATIVES = FMU_NUM_STATES
FMU_NUM_OUTPUTS = 0
FMU_NUM_INPUTS = 0
FMU_NUM_EVENTS = 0
FMU_NUM_PARAMETERS = 0

FMU_FCT_INIT = function() return ([],[],[],[],[]) end
FMU_FCT_EVALUATE = function(t, x, ẋ, u, p) return ([], [], []) end
FMU_FCT_OUTPUT = function(t, x, ẋ, u, p) return [] end
FMU_FCT_EVENT = function(t, x, ẋ, u, p) return [] end

##############

function logInfo(_component::fmi2Component, status::fmi2Status, message)
    component = dereferenceComponent(_component)
    if component != nothing
        if component.loggingOn == fmi2True
            ccall(component.callbackFunctions.logger, Cvoid, (fmi2ComponentEnvironment, fmi2String, fmi2Status, fmi2String, fmi2String), component.callbackFunctions.componentEnvironment, component.instanceName, status, "info", message * "\n")
        end
    end
end

function logWarning(_component::fmi2Component, status::fmi2Status, message)
    component = dereferenceComponent(_component)
    if component != nothing
        ccall(component.callbackFunctions.logger, Cvoid, (fmi2ComponentEnvironment, fmi2String, fmi2Status, fmi2String, fmi2String), component.callbackFunctions.componentEnvironment, component.instanceName, status, "warning", message * "\n")
    end
end

function logError(_component::fmi2Component, status::fmi2Status, message)
    component = dereferenceComponent(_component)
    if component != nothing
        ccall(component.callbackFunctions.logger, Cvoid, (fmi2ComponentEnvironment, fmi2String, fmi2Status, fmi2String, fmi2String), component.callbackFunctions.componentEnvironment, component.instanceName, status, "error", message * "\n")
    end
end

function dereferenceComponent(addr::fmi2Component)
    global FMIBUILD_FMU
    for component in FMIBUILD_FMU.components
        if addr == component.compAddr
            return component
        end
    end

    @warn "Unknown fmi2Component at $(addr)."
    return nothing
end

function reset(_component::fmi2Component)
    component = dereferenceComponent(_component)

    (component.t, x, ẋ, u, p) = FMU_FCT_INIT()
    component.z = FMU_FCT_EVENT(component.t, x, ẋ, u, p)
    y = FMU_FCT_OUTPUT(component.t, x, ẋ, u, p)

    applyValues(component.compAddr, x, ẋ, u, y, p)
end

function evaluate(_component::fmi2Component)
    component = dereferenceComponent(_component)

    x, ẋ, u, y, p = extractValues(_component)

    x, ẋ, p = FMU_FCT_EVALUATE(component.t, x, ẋ, u, p)
    y       = FMU_FCT_OUTPUT(  component.t, x, ẋ, u, p)
    z_new   = FMU_FCT_EVENT(   component.t, x, ẋ, u, p)

    applyValues(_component, x, ẋ, u, y, p)

    # event triggering

    if (component.z_prev != nothing) && (sign.(z_new) != sign.(component.z_prev))
        component.continuousStatesChanged = fmi2True
        component.z_prev = z_new
        component.z      = z_new
    else
        component.z_prev = component.z
        component.z      = z_new
    end
end

function applyValues(_component::fmi2Component, x, ẋ, u, y, p)
    component = dereferenceComponent(_component)

    for i in 1:length(FMIBUILD_FMU.modelDescription.stateValueReferences)
        vr = FMIBUILD_FMU.modelDescription.stateValueReferences[i]
        component.realValues[vr] = x[i]
    end

    for i in 1:length(FMIBUILD_FMU.modelDescription.derivativeValueReferences)
        vr = FMIBUILD_FMU.modelDescription.derivativeValueReferences[i]
        component.realValues[vr] = ẋ[i]
    end

    for i in 1:length(FMIBUILD_FMU.modelDescription.inputValueReferences)
        vr = FMIBUILD_FMU.modelDescription.inputValueReferences[i]
        component.realValues[vr] = u[i]
    end

    for i in 1:length(FMIBUILD_FMU.modelDescription.outputValueReferences)
        vr = FMIBUILD_FMU.modelDescription.outputValueReferences[i]
        component.realValues[vr] = y[i]
    end

    for i in 1:length(FMIBUILD_FMU.modelDescription.parameterValueReferences)
        vr = FMIBUILD_FMU.modelDescription.parameterValueReferences[i]
        component.realValues[vr] = p[i]
    end

    nothing
end

function extractValues(_component::fmi2Component)
    component = dereferenceComponent(_component)

    xs = collect(component.realValues[vr] for vr in component.fmu.modelDescription.stateValueReferences)
    ẋs = collect(component.realValues[vr] for vr in component.fmu.modelDescription.derivativeValueReferences)
    us = collect(component.realValues[vr] for vr in component.fmu.modelDescription.inputValueReferences)
    ys = collect(component.realValues[vr] for vr in component.fmu.modelDescription.outputValueReferences)
    ps = collect(component.realValues[vr] for vr in component.fmu.modelDescription.parameterValueReferences)

    return xs, ẋs, us, ys, ps
end

const STRING_TYPES_PLATFORM_DEFAULT = "default"
const STRING_VERSION_2_0 = "2.0"

# 2.1.4
function simple_fmi2GetTypesPlatform()
    global STRING_TYPES_PLATFORM_DEFAULT
    return pointer(STRING_TYPES_PLATFORM_DEFAULT)
end

function simple_fmi2GetVersion()
    global STRING_VERSION_2_0
    return pointer(STRING_VERSION_2_0)
end

# 2.1.5
function simple_fmi2Instantiate(instanceName::fmi2String,
                                         fmuType::fmi2Type,
                                         fmuGUID::fmi2String,
                                         fmuResourceLocation::fmi2String,
                                         functions::Ptr{fmi2CallbackFunctions},
                                         visible::fmi2Boolean,
                                         loggingOn::fmi2Boolean)

    global FMIBUILD_FMU

    component = FMU2Component()
    component.loggingOn = loggingOn
    component.callbackFunctions = unsafe_load(functions)
    component.instanceName = unsafe_string(instanceName)

    component.compAddr = pointer_from_objref(component)
    component.fmu = FMIBUILD_FMU
    push!(FMIBUILD_FMU.components, component)

    reset(component.compAddr)

    logInfo(component.compAddr, fmi2StatusOK, "fmi2Instantiate: New fmi2Component at $(component.compAddr)")

    return component.compAddr
end

function simple_fmi2FreeInstance(_component::fmi2Component)
    component = dereferenceComponent(_component)

    if component != nothing
        logInfo(component.compAddr, fmi2StatusOK, "fmi2FreeInstance")

        global FMIBUILD_FMU
        for i in 1:length(FMIBUILD_FMU.components)
            if FMIBUILD_FMU.components[i].compAddr == component.compAddr
                deleteat!(FMIBUILD_FMU.components, i)
                break
            end
        end
        logInfo(component.compAddr, fmi2StatusOK, "fmi2FreeInstance: Freed fmi2Component at $(component.compAddr).")
    end

    nothing
end

function simple_fmi2SetDebugLogging(_component::fmi2Component, loggingOn::fmi2Boolean, nCategories::Csize_t, categories::Ptr{fmi2String})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetDebugLogging")

    component.loggingOn = loggingOn
    # ToDo: categories (arguments)

    return fmi2StatusOK
end

function simple_fmi2SetupExperiment(_component::fmi2Component, toleranceDefined::fmi2Boolean, tolerance::fmi2Real, startTime::fmi2Real, stopTimeDefined::fmi2Boolean, stopTime::fmi2Real)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetupExperiment")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2EnterInitializationMode(_component::fmi2Component)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2EnterInitializationMode")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2ExitInitializationMode(_component::fmi2Component)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2ExitInitializationMode")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2Terminate(_component::fmi2Component)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2Terminate")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2Reset(_component::fmi2Component)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2Reset")

    reset(_component)

    return fmi2StatusOK
end

function simple_fmi2GetReal(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2Real})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetReal")

    value = unsafe_wrap(Array{fmi2Real}, _value, nvr)
    vr = unsafe_wrap(Array{fmi2ValueReference}, _vr, nvr)

    global FMU_NUM_STATES, FMU_NUM_OUTPUTS, FMU_NUM_INPUTS, FMU_NUM_PARAMETERS

    for i in 1:nvr
        valueRef = vr[i]

        try
            value[i] = component.realValues[valueRef]
        catch e
            logError(component.compAddr, fmi2StatusError, "fmi2SetReal: Unknown value reference $(valueRef).\n")
            return fmi2StatusError
        end
    end

    return fmi2StatusOK
end

function simple_fmi2GetInteger(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2Integer})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetInteger")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2GetBoolean(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2Boolean})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetBoolean")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2GetString(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2String})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetString")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2SetReal(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2Real})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetReal")

    value = unsafe_wrap(Array{fmi2Real}, _value, nvr)
    vr = unsafe_wrap(Array{fmi2ValueReference}, _vr, nvr)

    for i in 1:nvr
        valueRef = vr[i]

        try
            component.realValues[valueRef] = value[i]
        catch e
            logError(component.compAddr, fmi2StatusError, "fmi2SetReal: Unknown value reference $(valueRef).\n")
            return fmi2StatusError
        end
    end

    evaluate(_component)

    return fmi2StatusOK
end

function simple_fmi2SetInteger(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2Integer})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetInteger")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2SetBoolean(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2Boolean})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetBoolean")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2SetString(_component::fmi2Component, _vr::Ptr{fmi2ValueReference}, nvr::Csize_t, _value::Ptr{fmi2String})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetString")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2SetTime(_component::fmi2Component, time::fmi2Real)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetTime")

    component.t = time

    evaluate(_component)

    return fmi2StatusOK
end

function simple_fmi2SetContinuousStates(_component::fmi2Component, _x::Ptr{fmi2Real}, nx::Csize_t)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2SetContinuousStates")

    if nx != length(component.fmu.modelDescription.stateValueReferences)
        logWarning(component.compAddr, fmi2StatusWarning, "fmi2SetContinuousStates: Model has $(length(component.fmu.modelDescription.stateValueReferences)) states, but `nx`=$(nx).")
    end

    x = unsafe_wrap(Array{fmi2Real}, _x, nx)

    for i in 1:nx
        vr = component.fmu.modelDescription.stateValueReferences[i]
        component.realValues[vr] = x[i]
    end

    evaluate(_component)

    return fmi2StatusOK
end

function simple_fmi2EnterEventMode(_component::fmi2Component)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2EnterEventMode")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2NewDiscreteStates(_component::fmi2Component, _fmi2eventInfo::Ptr{fmi2EventInfo})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2NewDiscreteStates")

    # ToDo: This is not efficient (copy struct and overwrite), direct memory access would be much nicer!
    eventInfo = unsafe_load(_fmi2eventInfo)
    eventInfo.newDiscreteStatesNeeded = fmi2False
    eventInfo.terminateSimulation = fmi2False
    eventInfo.nominalsOfContinuousStatesChanged = fmi2False
    eventInfo.valuesOfContinuousStatesChanged = component.continuousStatesChanged
    eventInfo.nextEventTimeDefined = fmi2False
    eventInfo.nextEventTime = 0.0
    unsafe_store!(_fmi2eventInfo, eventInfo);

    component.continuousStatesChanged = fmi2False

    return fmi2StatusOK
end

function simple_fmi2EnterContinuousTimeMode(_component::fmi2Component)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2EnterContinuousTimeMode")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2CompletedIntegratorStep(_component::fmi2Component, noSetFMUStatePriorToCurrentPoint::fmi2Boolean, enterEventMode::Ptr{fmi2Boolean}, terminateSimulation::Ptr{fmi2Boolean})
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2CompletedIntegratorStep")

    # ToDo

    return fmi2StatusOK
end

function simple_fmi2GetDerivatives(_component::fmi2Component, _derivatives::Ptr{fmi2Real}, nx::Csize_t)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetDerivatives")

    if nx != length(component.fmu.modelDescription.derivativeValueReferences)
        logWarning(component.compAddr, fmi2StatusWarning, "fmi2GetDerivatives: Model has $(length(component.fmu.modelDescription.derivativeValueReferences)) states, but `nx`=$(nx).")
    end

    derivatives = unsafe_wrap(Array{fmi2Real}, _derivatives, nx)

    for i in 1:nx
        vr = component.fmu.modelDescription.derivativeValueReferences[i]
        derivatives[i] = component.realValues[vr]
    end

    return fmi2StatusOK
end

function simple_fmi2GetEventIndicators(_component::fmi2Component, _eventIndicators::Ptr{fmi2Real}, ni::Csize_t)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetEventIndicators")

    if ni != length(component.fmu.modelDescription.numberOfEventIndicators)
        logWarning(component.compAddr, fmi2StatusWarning, "fmi2GetEventIndicators: Model has $(length(component.eventIndicators)) states, but `ni`=$(ni).")
    end

    eventIndicators = unsafe_wrap(Array{fmi2Real}, _eventIndicators, ni)

    for i in 1:ni
        eventIndicators[i] = component.z[i]
    end

    return fmi2StatusOK
end

function simple_fmi2GetContinuousStates(_component::fmi2Component, _x::Ptr{fmi2Real}, nx::Csize_t)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetContinuousStates")

    if nx != length(component.fmu.modelDescription.stateValueReferences)
        logWarning(component.compAddr, fmi2StatusWarning, "fmi2GetContinuousStates: Model has $(length(component.fmu.modelDescription.stateValueReferences)) states, but `nx`=$(nx).")
    end

    x = unsafe_wrap(Array{fmi2Real}, _x, nx)

    for i in 1:nx
        vr = component.fmu.modelDescription.stateValueReferences[i]
        x[i] = component.realValues[vr]
    end

    return fmi2StatusOK
end

function simple_fmi2GetNominalsOfContinuousStates(_component::fmi2Component, _x_nominal::Ptr{fmi2Real}, nx::Csize_t)
    component = dereferenceComponent(_component)
    logInfo(component.compAddr, fmi2StatusOK, "fmi2GetNominalsOfContinuousStates")

    if nx != length(component.fmu.modelDescription.stateValueReferences)
        logWarning(component.compAddr, fmi2StatusWarning, "fmi2GetNominalsOfContinuousStates: Model has $(length(component.fmu.modelDescription.stateValueReferences)) states, but `nx`=$(nx).")
    end

    x_nominal = unsafe_wrap(Array{fmi2Real}, _x_nominal, nx)

    for i in 1:nx
        vr = component.fmu.modelDescription.stateValueReferences[i]
        x_nominal[i] = component.realValues[vr]
    end

    return fmi2StatusOK
end

"""
    initializationFct           # () -> (t, x, ẋ, u, p)
    evaluationFct               # (t, x, ẋ, u, p) -> (x, ẋ, p)
    outputFct                   # (t, x, ẋ, u, p) -> y
    eventFct                    # (t, x, ẋ, u, p) -> e
"""
function fmi2CreateSimple(;
    initializationFct=nothing,
    evaluationFct=nothing,
    outputFct=nothing,
    eventFct=nothing)

    global FMIBUILD_FMU

    global FMU_FCT_INIT
    global FMU_FCT_EVALUATE
    global FMU_FCT_OUTPUT
    global FMU_FCT_EVENT

    global FMU_NUM_STATES
    global FMU_NUM_OUTPUTS
    global FMU_NUM_INPUTS
    global FMU_NUM_EVENTS
    global FMU_NUM_PARAMETERS

    FMIBUILD_FMU = fmi2Create()

    FMU_FCT_INIT = initializationFct
    FMU_FCT_EVALUATE = evaluationFct
    FMU_FCT_OUTPUT = outputFct
    FMU_FCT_EVENT = eventFct

    t, x, ẋ, u, p = FMU_FCT_INIT()
    y = FMU_FCT_OUTPUT(t, x, ẋ, u, p)
    e = FMU_FCT_EVENT(t, x, ẋ, u, p)

    FMU_NUM_STATES = length(x)
    FMU_NUM_OUTPUTS = length(y)
    FMU_NUM_INPUTS = length(u)
    FMU_NUM_EVENTS = length(e)
    FMU_NUM_PARAMETERS  = length(p)

    fmi2SetFctGetVersion(FMIBUILD_FMU,                     simple_fmi2GetVersion)
    fmi2SetFctGetTypesPlatform(FMIBUILD_FMU,               simple_fmi2GetTypesPlatform)
    fmi2SetFctInstantiate(FMIBUILD_FMU,                    simple_fmi2Instantiate)
    #FMIBUILD_FMU.cInstantiate = @cfunction(simple_fmi2Instantiate, fmi2Component, (fmi2String, fmi2Type, fmi2String, fmi2String, Ptr{fmi2CallbackFunctions}, fmi2Boolean, fmi2Boolean))
    fmi2SetFctFreeInstance(FMIBUILD_FMU,                   simple_fmi2FreeInstance)
    fmi2SetFctSetDebugLogging(FMIBUILD_FMU,                simple_fmi2SetDebugLogging)
    fmi2SetFctSetupExperiment(FMIBUILD_FMU,                simple_fmi2SetupExperiment)
    fmi2SetEnterInitializationMode(FMIBUILD_FMU,           simple_fmi2EnterInitializationMode)
    fmi2SetFctExitInitializationMode(FMIBUILD_FMU,         simple_fmi2ExitInitializationMode)
    fmi2SetFctTerminate(FMIBUILD_FMU,                      simple_fmi2Terminate)
    fmi2SetFctReset(FMIBUILD_FMU,                          simple_fmi2Reset)
    fmi2SetFctGetReal(FMIBUILD_FMU,                        simple_fmi2GetReal)
    fmi2SetFctGetInteger(FMIBUILD_FMU,                     simple_fmi2GetInteger)
    fmi2SetFctGetBoolean(FMIBUILD_FMU,                     simple_fmi2GetBoolean)
    fmi2SetFctGetString(FMIBUILD_FMU,                      simple_fmi2GetString)
    fmi2SetFctSetReal(FMIBUILD_FMU,                        simple_fmi2SetReal)
    fmi2SetFctSetInteger(FMIBUILD_FMU,                     simple_fmi2SetInteger)
    fmi2SetFctSetBoolean(FMIBUILD_FMU,                     simple_fmi2SetBoolean)
    fmi2SetFctSetString(FMIBUILD_FMU,                      simple_fmi2SetString)
    fmi2SetFctSetTime(FMIBUILD_FMU,                        simple_fmi2SetTime)
    fmi2SetFctSetContinuousStates(FMIBUILD_FMU,            simple_fmi2SetContinuousStates)
    fmi2SetFctEnterEventMode(FMIBUILD_FMU,                 simple_fmi2EnterEventMode)
    fmi2SetFctNewDiscreteStates(FMIBUILD_FMU,              simple_fmi2NewDiscreteStates)
    fmi2SetFctEnterContinuousTimeMode(FMIBUILD_FMU,        simple_fmi2EnterContinuousTimeMode)
    fmi2SetFctCompletedIntegratorStep(FMIBUILD_FMU,        simple_fmi2CompletedIntegratorStep)
    fmi2SetFctGetDerivatives(FMIBUILD_FMU,                 simple_fmi2GetDerivatives)
    fmi2SetFctGetEventIndicators(FMIBUILD_FMU,             simple_fmi2GetEventIndicators)
    fmi2SetFctGetContinuousStates(FMIBUILD_FMU,            simple_fmi2GetContinuousStates)
    fmi2SetFctGetNominalsOfContinuousStates(FMIBUILD_FMU,  simple_fmi2GetNominalsOfContinuousStates)

    return FMIBUILD_FMU
end


text:
module Garch

using Random, Statistics, Optim, Parameters

export Std, GJR, fit!, simulate!, predict!

@with_kw mutable struct Std
    @deftype Float64
    σ = 0.01
    @assert σ > 0.0
    μ₁ = rand()
    @assert μ₁ > 0.0 && μ₁ < 1.0
    μ₂ = rand()
    @assert μ₂ > 0.0 && μ₂ < 1.0
end

@with_kw mutable struct GJR
    @deftype Float64
    σ = 0.01
    @assert σ > 0.0
    μ₁ = rand()
    @assert μ₁ > 0.0 && μ₁ < 1.0
    μ₂ = rand()
    @assert μ₂ > 0.0 && μ₂ < 1.0
    μ₃ = rand()
    @assert μ₃ > 0.0 && μ₃ < 1.0
    isdown::Bool = true
end

function nll(z, σ, rt)
    μ₁, μ₂ = exp.(-exp.(-z))
    ll = 0.0
    n = length(rt)
    σ² = σ * σ
    h = σ²
    μ = 1.0 - μ₂
    for t ∈ 2:n
        h = σ² + μ₁ * (μ₂ * h + μ * rt[t-1] * rt[t-1] - σ²)
        if h > 0.0
            ll += 0.5 * (log(h) + rt[t] * rt[t] / h)
        end
    end
    return ll
end

function nll(z, σ, rt, J)
    μ₁, μ₂, μ₃ = exp.(-exp.(-z))
    ll = 0.0
    n = length(rt)
    σ² = σ * σ
    h = σ²
    for t ∈ 2:n
        h = σ² + μ₁ * ((1.0 - μ₂ + 2.0 * (μ₂ - μ₃) * J[t-1]) * rt[t-1] * rt[t-1] + μ₃ * h - σ²)
        if h > 0.0
            ll += 0.5 * (log(h) + rt[t] * rt[t] / h)
        end
    end
    return ll
end

function fit!(model::GJR, x::Vector{Float64}; time_limit = 60.0, solver = NelderMead(), g_tol = 1.0e-6)
    @unpack σ, μ₁, μ₂, μ₃, isdown = model
    σ = std(x)
    J = 1.0 .* (isdown ? (x .< 0.0) : (x .> 0.0))
    p = [-log(-log(μ₁)), -log(-log(μ₂)), -log(-log(μ₃))]
    p0 = 4.0 .* rand(3)
    res = optimize(p -> nll(p, σ, x, J), p0, solver, Optim.Options(time_limit = time_limit, g_tol = g_tol))
    μ₁, μ₂ = exp.(-exp.(-res.minimizer))
    @pack! model = σ, μ₁, μ₂, μ₃, isdown
    return res
end

function fit!(model::Std, x::Vector{Float64}; time_limit = 60.0, solver = NelderMead(), g_tol = 1.0e-6)
    @unpack σ, μ₁, μ₂ = model
    σ = std(x)
    p = [-log(-log(μ₁)), -log(-log(μ₂))]
    p0 = 4.0 .* rand(2)
    res = optimize(p -> nll(p, σ, x), p0, solver, Optim.Options(time_limit = time_limit, g_tol = g_tol))
    μ₁, μ₂ = exp.(-exp.(-res.minimizer))
    @pack! model = σ, μ₁, μ₂
    return res
end

function simulate(model::Std, pathlen::Integer; r0 = 0.0)
    r = fill(0.0, pathlen + 1)
    r[1] = r0
    @unpack σ, μ₁, μ₂ = model
    σ² = σ * σ
    h = σ²
    μ = 1.0 - μ₂
    for t ∈ 1:pathlen
        h = max(0.0, σ² + μ₁ * (μ₂ * h + μ * r[t] * r[t] - σ²))
        r[t+1] = √h * randn()
    end
    return r[2:end]
end

function simulate(model::GJR, pathlen::Integer; r0 = 0.0)
    r = fill(0.0, pathlen + 1)
    r[1] = r0
    @unpack σ, μ₁, μ₂, μ₃, isdown = model
    σ² = σ * σ
    h = σ²
    for t ∈ 1:pathlen
        J = isdown ? (r[t] < 0.0) : (r[t] > 0.0)
        h = max(σ² + μ₁ * ((1.0 - μ₂ + 2.0 * (μ₂ - μ₃) * J) * r[t] * r[t] + μ₃ * h - σ²), 0.0)
        r[t+1] = √h * randn()
    end
    return r[2:end]
end

function predict(model::Union{GJR,Std}, paths::Integer, pathlen::Integer)
    vol = fill(0.0, paths)
    vol .= [std(simulate(model, pathlen, r0 = 0.0)) for i ∈ 1:paths]
    return vol
end

end

ext: jl


text:
export RandomVertexCover

struct RandomVertexCover end

"""
    vertex_cover(g, RandomVertexCover())

Find a set of vertices such that every edge in `g` has some vertex in the set as
atleast one of its end point.

### Implementation Notes
Performs [Approximate Minimum Vertex Cover](https://en.wikipedia.org/wiki/Vertex_cover#Approximate_evaluation) once.
Returns a vector of vertices representing the vertices in the Vertex Cover.

### Performance
Runtime: O(|V|+|E|)
Memory: O(|E|)
Approximation Factor: 2
"""
function vertex_cover(
    g::AbstractGraph{T},
    alg::RandomVertexCover
    ) where T <: Integer

    nvg::T = nv(g)
    cover = Vector{T}()
    in_cover = falses(nvg)

    @inbounds for e in shuffle(collect(edges(g)))
        if !(in_cover[e.src] || in_cover[e.dst])
            in_cover[e.src] = in_cover[e.dst] = true
            push!(cover, e.src)
            push!(cover, e.dst)
        end
    end

    return [v for v in vertices(g) if in_cover[v]]
end


text:
function gaussjacobi(n::Int32, a, b)
#GAUSS-JACOBI QUADRATURE NODES AND WEIGHTS

    if ( a == 0 && b == 0 )
        x = gausslegendre(n)
    elseif ( a == -0.5 && b == -0.5 )
        x = gausschebyshev(n,1)
    elseif ( a == 0.5 && b == 0.5)
        x = gausschebyshev(n,2)
    elseif ( a == -0.5 && b == 0.5)
        x = gausschebyshev(n,3)
    elseif ( a == 0.5 && b == -0.5 )
        x = gausschebyshev(n,4)
    elseif ( n == 0 )
        x = (Float64[],Float64[])
    elseif ( n == 1 )
        x = ([(b-a)/(a+b+2)], [2^(a+b+1)*beta(a+1, b+1)])
    elseif ( n <= 100 )
        x = JacobiRec(n, a, b)
    elseif ( n > 100 )
        x = JacobiAsy(n, a, b)
    else
        error("1st argument must be a positive integer.")
    end
    return x
end

function JacobiRec(n::Int32, a::Float64, b::Float64)
#Compute nodes and weights using recurrrence relation.
    x1 = HalfRec(n, a, b, 1)
    x2 = HalfRec(n, b, a, 0)
    x = vcat( -flipud(x2[1]) , x1[1] )
    ders = vcat( flipud(x2[2]) , x1[2] )
    w = 1./((1-x.^2).*ders.^2)
    w = 2^(a+b+1)*gamma(2+a) * gamma(2+b) / (gamma(2+a+b)*(a+1)*(b+1)) * w ./ sum(w)
    return x, w
end

function HalfRec(n::Int32, a::Float64, b::Float64, flag::Int32)
#HALFREC   Jacobi polynomial recurrence relation.
    # Asymptotic formula - only valid for positive x.
    r = (flag==1) ? [ceil(n/2):-1:1] : [floor(n/2):-1:1]
    C = (2*r+a-.5)*pi/(2*n+a+b+1)
    x = cos( C + 1/(2*n+a+b+1)^2 * ((.25-a^2)*cot(.5*C) - (.25-b^2)*tan(.5*C)) )
    dx = 1.0
    counter = 0
    # Loop until convergence:
    while ( norm(dx,Inf) > sqrt(eps(Float64))/1000 && counter < 10 )
        counter = counter + 1
        P = innerJacobiRec(n, x, a, b)
        dx = -P[1]./P[2]
        x += dx
    end
    # Once more for derivatives:
    P = innerJacobiRec(n, x, a, b)
    return x, P[2]
end

function innerJacobiRec(n::Int32, x::Array{Float64}, a::Float64, b::Float64)
# EVALUATE JACOBI POLYNOMIALS AND ITS DERIVATIVE USING THREE-TERM RECURRENCE.
    #P, Pm1, PP, PPm1 = [.5*(a-b+(a+b+2)*x)], ones(x), .5*(a+b+2)*ones(x), zeros(x)
    N = length(x)
    P = Array(Float64,N); Pm1 = Array(Float64,N); PP = Array(Float64,N); PPm1 = Array(Float64,N);
    for j = 1:N
        P[j] = .5*(a-b+(a+b+2)*x[j]); Pm1[j] = 1.0; PP[j] = .5*(a+b+2); PPm1[j] = 0.0
        for k = 1:n-1
            A = 2*(k + 1)*(k + a+b + 1)*(2*k + a+b)
            B = (2*k + a+b + 1)*(a^2 - b^2)
            C = prod(2*k + a+b + [0:2])
            D = 2*(k + a)*(k + b)*(2*k + a+b + 2)
            Pm1[j], P[j] = P[j], ( (B+C*x[j])*P[j] - D*Pm1[j] ) / A
            PPm1[j], PP[j] = PP[j], ( (B+C*x[j])*PP[j] + C*Pm1[j] - D*PPm1[j] ) / A
        end
    end
    return P, PP
end

function weightsConstant( n::Int32, a::Float64, b::Float64 )
    # Compute the constant for weights:
    M = min(20, n-1)
    C = 1.0
    p::Float64 = -a*b/n
    for m = 1:M
        C += p
        p *= -(m+a)*(m+b)/(m+1)/(n-m)
        if ( abs(p / C) < eps(Float64)/100 ) break end
    end
    return 2^(a+b+1)*C
end

function JacobiAsy(n::Int32, a::Float64, b::Float64)
#ASY   Compute nodes and weights using asymptotic formulae.

    # Determine switch between interior and boundary regions:
    nbdy = 10
    bdyidx1 = n-(nbdy-1):n
    bdyidx2 = nbdy:-1:1

    # Interior formula:
    x = asy1(n, a, b, nbdy)
    w = x[2]; x = x[1]

    # Boundary formula (right):
    xbdy = boundary(n, a, b, nbdy);
    x[bdyidx1] = xbdy[1]
    w[bdyidx1] = xbdy[2]

    # Boundary formula (left):
    if ( !(a == b) )
        xbdy = boundary(n, b, a, nbdy)
    end
    x[bdyidx2] = -xbdy[1]
    w[bdyidx2] = xbdy[2]

    w *= weightsConstant(n, a, b)
    return x[:], w[:]
end

function asy1(n::Int32, a::Float64, b::Float64, nbdy)
# Algorithm for computing nodes and weights in the interior.

    # Approximate roots via asymptotic formula: (Gatteschi and Pittaluga, 1985)
    K = (2*[n:-1:1]+a-.5)*pi/(2*n+a+b+1)
    tt = K + 1/(2*n+a+b+1)^2*((.25-a^2)*cot(.5*K)-(.25-b^2)*tan(.5*K))

    # First half (x > 0):
    t = tt[tt .<= pi/2]'
    mint = t[end-nbdy+1]
    idx = 1:max(findfirst(float64(t .< mint))-1, 1)

    dt = 1.0; counter = 0
    # Newton iteration
    while ( norm(dt,Inf) > sqrt(eps(Float64))/100 && counter < 10 )
        vals = feval_asy1(n, a, b, t, idx, 0)         # Evaluate
        dt = vals[1]./vals[2]                                 # Newton update
        t += dt                                     # Next iterate
        counter += 1
        dt = dt[idx]
    end
    vals = feval_asy1(n, a, b, t, idx, 0)      # Once more for luck
    t += vals[1]./vals[2]                                 # Newton update.

    # Store:
    x = cos(t)
    w = 1./vals[2].^2

    # Second half (x < 0):
    tmp = a; a = b; b = tmp
    t = pi - tt[1:(n-length(x))]'
    mint = t[nbdy]
    idx = max(findfirst(float64(t .> mint)), 1):length(t)

    dt = 1.0; counter = 0;
    # Newton iteration
    while ( norm(dt,Inf) > sqrt(eps(Float64))/100 && counter < 10 )
        vals = feval_asy1(n, a, b, t, idx, 0)  # Evaluate.
        dt = vals[1]./vals[2]                                # Newton update.
        t += dt                                     # Next iterate.
        counter += 1
        dt = dt[idx]
    end
    vals = feval_asy1(n, a, b, t, idx, 0)     # Once more for luck.
    t += vals[1]./vals[2]                                 # Newton update.

    # Store:
    x = [-cos(t) x];
    w = [1./vals[2].^2 w];
    return x, w
end

function feval_asy1(n::Int32, a::Float64, b::Float64, t, idx, flag)
# Evaluate the interior asymptotic formula at x = cos(t).

    # Number of terms in the expansion:
    M = 20

    # Some often used vectors/matrices:
    onesT = ones(1,length(t)); onesM = ones(M); MM = ([0:M-1]')';

    # The sine and cosine terms:
    alpha = (.5*(2*n+a+b+1+MM))*onesT .* (onesM*t) - .5*(a+.5)*pi
    cosA = cos(alpha); sinA = sin(alpha)

    if ( flag == 1 )
        k = ( idx[1] == 1 ) ? [length(t):-1:1]' : [1:length(t)]'
        ta = float64(float32(t))
        tb = t - ta; hi = n*ta; lo = n*tb + (a+b+1)*.5*t
        pia = float64(float32(pi))
        pib = -8.742278000372485e-08; #pib = pi - pia
        dh = (hi - (k-.25)*pia) + lo - .5*a*pia - (k-.25+.5*a)*pib
        tmp = 0; sgn = 1; fact = 1; DH = dh; dh2 = dh.*dh
        for j = 0:20
            dc = sgn*DH/fact
            tmp = tmp + dc
            sgn = -sgn
            fact = fact*(2*j+3)*(2*j+2)
            DH = DH.*dh2
            if ( norm(dc,Inf) < eps(Float64)/2000 )
                break
            end
        end
        tmp = -tmp; tmp[1] = -tmp[1]
        tmp = sign(cosA[1,2]*tmp[2])*tmp
        cosA[1,:] = tmp
    end

    sinT = onesM*sin(t)
    cosT = onesM*cos(t)
    cosA2 = cosA.*cosT + sinA.*sinT
    sinA2 = sinA.*cosT - cosA.*sinT

    one = ones(t)
    sinT = vcat( one , cumprod(onesM[2:end]*(.5*csc(.5*t))))
    cosT = .5*sec(.5*t)

    j = transpose([0:M-2])
    vec = (.5+a+j).*(.5-a+j)./(j+1)./(2*n+a+b+j+2)
    P1 = [1 cumprod(vec,2)]
    P1[3:4:end] = -P1[3:4:end]
    P1[4:4:end] = -P1[4:4:end]
    P2 = eye(M)
    for l = 0:M-1
        j = transpose([0:(M-l-2)])
        vec = (.5+b+j).*(.5-b+j)./(j+1)./(2*n+a+b+j+l+2)
        P2[l+1+(1:length(j)),l+1] = cumprod(vec,2)
    end
    PHI = repmat(P1,M,1).*P2

    j = transpose([0:M-2])
    vec = (.5+a+j).*(.5-a+j)./(j+1)./(2*(n-1)+a+b+j+2)
    P1 = [1 cumprod(vec,2)]
    P1[3:4:end] = -P1[3:4:end]
    P1[4:4:end] = -P1[4:4:end]
    P2 = eye(M)
    for l = 0:M-1
        j = transpose([0:(M-l-2)])
        vec = (.5+b+j).*(.5-b+j)./(j+1)./(2*(n-1)+a+b+j+l+2)
        P2[l+1+(1:length(j)),l+1] = cumprod(vec,2)
    end
    PHI2 = repmat(P1,M,1).*P2

    S = 0; S2 = 0;
    SC = sinT
    for m = 0:M-1
        l = [0:2:m]
        phi = PHI[m+1,l+1]
        dS1 = phi*SC[l+1,:].*cosA[m+1,:]
        phi2 = PHI2[m+1,l+1]
        dS12 = phi2*SC[l+1,:].*cosA2[m+1,:]
        l = [1:2:m]
        phi = PHI[m+1,l+1]
        dS2 = phi*SC[l+1,:].*sinA[m+1,:]
        phi2 = PHI2[m+1,l+1]
        dS22 = phi2*SC[l+1,:].*sinA2[m+1,:]
        if m > 10 && norm(dS1[idx] + dS2[idx], Inf) < eps(Float64)/100 break end
        S = S + dS1 + dS2
        S2 = S2 + dS12 + dS22
        SC[1:m+1,:] = broadcast(.*,SC[1:m+1,:],cosT)
    end

    # Constant out the front:
    dsa = .5*(a^2)/n; dsb = .5*(b^2)/n; dsab = .25*(a+b)^2/n
    ds = dsa + dsb - dsab; s = ds; j = 1
    dsold = ds # to fix a = -b bug.
    while ( (abs(ds/s) + dsold) > eps(Float64)/10 )
        dsold = abs(ds/s)
        j += 1
        tmp = -(j-1)/(j+1)/n
        dsa = tmp*dsa*a
        dsb = tmp*dsb*b
        dsab = .5*tmp*dsab*(a+b)
        ds = dsa + dsb - dsab
        s = s + ds
    end
    p2 = exp(s)*sqrt(2*pi)*sqrt((n+a)*(n+b)/(2*n+a+b))/(2*n+a+b+1)
    g = [1, 1/12, 1/288, -139/51840, -571/2488320, 163879/209018880,
         5246819/75246796800, -534703531/902961561600,
         -4483131259/86684309913600, 432261921612371/514904800886784000]
    f(g,z) = sum(g.*[1,cumprod(ones(9)./z)])
    C = p2*(f(g,n+a)*f(g,n+b)/f(g,2*n+a+b))*2/pi
    C2 = C*(a+b+2*n).*(a+b+1+2*n)./(4*(a+n).*(b+n))

    vals = C*S
    S2 = C2*S2

    # Use relation for derivative:
    ders = (n*(a-b-(2*n+a+b)*cos(t)).*vals + 2*(n+a)*(n+b)*S2)/(2*n+a+b)./sin(t)
    denom = 1./real(sin(t/2).^(a+.5).*cos(t/2).^(b+.5))
    vals = vals.*denom
    ders = ders.*denom

    return (vals, ders)
end

function boundary(n::Int32, a::Float64, b::Float64, npts)
# Algorithm for computing nodes and weights near the boundary.

    # Use Newton iterations to find the first few Bessel roots:
    smallK = min(30, npts);
    jk = besselRoots(a, min(npts, smallK))
    # Use asy formula for larger ones (See NIST 10.21.19, Olver 1974 p247)
    if ( npts > smallK )
        mu = 4*a^2;
        a8 = 8*([length(jk)+1:npts]'+.5*a-.25)*pi
        jk2 = .125*a8-(mu-1)./a8 - 4*(mu-1)*(7*mu-31)/3./a8.^3 -
          32*(mu-1)*(83*mu.^2-983*mu+3779)/15./a8.^5 -
          64*(mu-1)*(6949*mu^3-153855*mu^2+1585743*mu-6277237)/105./a8.^7
        jk = [jk ; jk2]
    end
    jk = real(jk[1:npts])

    # Approximate roots via asymptotic formula: (see Olver 1974)
    phik = jk/(n + .5*(a + b + 1))
    x = cos( phik + ((a^2-.25)*(1-phik.*cot(phik))./(8*phik) - .25*(a^2-b^2)*tan(.5*phik))/(n + .5*(a + b + 1))^2 )

    dx = 1.0; counter = 0;
    # Newton iteration:
    while ( norm(dx,Inf) > sqrt(eps(Float64))/200 && counter < 10)
        vals = innerJacobiRec(n, x, a, b)   # Evaluate via asymptotic formula.
        dx = -vals[1]./vals[2]                   # Newton update.
        x += dx                             # Next iterate.
        counter += 1
    end
    vals = innerJacobiRec(n, x, a, b);     # Evaluate via asymptotic formula.
    dx = -vals[1]./vals[2]                        # Newton update
    x += dx

    # flip:
    x = x[npts:-1:1]
    ders = vals[2]; vals = vals[1]
    ders = ders[npts:-1:1]

    # Revert to x-space:
    w = 1./((1-x.^2).*ders.^2)
    return x, w
end

function besselRoots(nu::Float64, m::Int32)
# BESSELROOTS(NU, M) returns the first M roots of besselj(nu, x).
# Find an approximation:
    jk = Array(Float64, m)
    m1 = 3
    if ( nu == 0 )
        xs = 2.404825557695773
    elseif ( nu > 0 )
        xs = nu + 1.8557*nu^(1/3)
    else
        nu1 = nu + 1;
        # See Piessens 1984:
        xs = 2*sqrt(nu+1)*(1 + nu1/4 - 7*nu1^2/96 + 49*nu1^3/1152 - 8363*nu1/276480);
        m1 = min(max(2*ceil(abs(log10(nu1))), 3), m);
    end

    jk[1] = besselNewton(nu, xs);
    if ( m > 1 )
        jk[2] = besselNewton(nu, jk[1]+.9*pi)
    end
    if ( m > 2)
        for k = 3:m1
            jk[k] = besselNewton(nu, jk[k-1]+.99*pi)
        end
        for k = m1+1:m
            jk[k] = besselNewton(nu, jk[k-1]+pi);
        end
    end
    return jk
end

function besselNewton(nu::Float64, x::Float64)
# BESSELNEWTON(NU, JK)   Find roots of Bessel function using Newton iteration.
    dx = 1.0
    counter = 0
    while ( (dx > sqrt(eps(Float64)/1000)) && counter < 20 )
        u = besselj(nu, x)
        du = besselj(nu-1, x) - nu./x*u
        dx = u./du
        x -= dx
        counter += 1
    end
    return x
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module NetworkManagementClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

# export models
export convert, Access
export convert, AuthenticationMethod
export convert, AzureAsyncOperationResult
export convert, CloudError
export convert, CloudErrorBody
export convert, Error
export convert, ErrorDetails
export convert, IPAllocationMethod
export convert, IPVersion
export convert, ManagedServiceIdentity
export convert, ManagedServiceIdentityUserAssignedIdentities
export convert, ProvisioningState
export convert, Resource
export convert, SubResource
export convert, TagsObject

# export operations
export convert

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
on_windows = @windows ? true : false
on_osx = @osx ? true : false

oldwdir = pwd()

nemo_pkgdir = Pkg.dir("Nemo")
nemo_wdir = Pkg.dir("Nemo", "deps")
nemo_vdir = Pkg.dir("Nemo", "local")

pkgdir = Pkg.dir("Hecke")
wdir = Pkg.dir("Hecke", "deps")
vdir = Pkg.dir("Hecke", "local")

if !ispath(Pkg.dir("Hecke", "local"))
    mkdir(Pkg.dir("Hecke", "local"))
end
if !ispath(Pkg.dir("Hecke", "local", "lib"))
    mkdir(Pkg.dir("Hecke", "local", "lib"))
end

function download_dll(url_string, location_string)
   try
      run(`curl -o $(location_string) -L $(url_string)`)
   catch
      download(url_string, location_string)
   end
end

LDFLAGS = "-Wl,-rpath,$nemo_vdir/lib -Wl,-rpath,\$\$ORIGIN/../share/julia/site/v$(VERSION.major).$(VERSION.minor)/Nemo/local/lib -Wl,-rpath,$vdir/lib -Wl,-rpath,\$\$ORIGIN/../share/julia/site/v$(VERSION.major).$(VERSION.minor)/Hecke/local/lib"
DLCFLAGS = "-fPIC -fno-common"

cd(wdir)

if on_windows
   if Int == Int64
      download_dll("http://www.mathematik.uni-kl.de/~thofmann/hecke/bin/libhecke.dll", joinpath(vdir, "lib", "libhecke.dll"))
   else
      println("There is no libhecke for 32 bit Windows")
   end
else
   cd("$wdir/hecke")
   withenv("LD_LIBRARY_PATH"=>"$vdir/lib:$nemo_vdir/lib", "LDFLAGS"=>LDFLAGS) do
      run(`./configure --prefix=$vdir --disable-static --enable-shared --with-mpir=$nemo_vdir --with-mpfr=$nemo_vdir --with-flint=$nemo_vdir`)
      run(`make -j4`)
      run(`make install`)
   end
end

push!(Libdl.DL_LOAD_PATH, Pkg.dir("Hecke", "local", "lib"))

cd(oldwdir)


text:
dir = @__DIR__
include(joinpath(dir,"test_function_optimisation.jl"))
include(joinpath(dir,"test_function_optimisation_noise.jl"))
include(joinpath(dir,"test_function_optimisation_multi_obj.jl"))
include(joinpath(dir,"test_function_optimisation_fast_options.jl"))
include(joinpath(dir,"test_function_optimisation_discontinuous.jl"))
include(joinpath(dir,"test_function_optimisation_constraint.jl"))
include(joinpath(dir,"test_function_optimisation_categorical.jl"))


text:
using HypothesisTests, Base.Test
using HypothesisTests: default_tail

# www.uni-siegen.de/phil/sozialwissenschaften/soziologie/mitarbeiter/ludwig-mayerhofer/statistik/statistik_downloads/statistik_ii_7.pdf
u5 = [620, 5350, 7220]
u250 = [3580, 4180, 5690]
u2500 = [3600, 3820, 3840, 3850, 4160, 5300, 6900, 7120, 9370]
more = [3920, 4520, 4760, 8560, 10350]
t = HypothesisTests.KruskalWallisTest(u5, u250, u2500, more)

@test t.n_i == [length(u5), length(u250), length(u2500), length(more)]
@test t.df == 3
@test t.R_i == [31, 25, 88, 66]
@test t.H ≈ 1.5803174603174597
@test t.tie_adjustment == 1
@test pvalue(t) ≈ 0.6638608922384397
@test default_tail(t) == :right
show(IOBuffer(), t)

# http://www.brightstat.com/index.php?option=com_content&task=view&id=41&Itemid=1&limit=1&limitstart=2
city1 = [68, 93, 123, 83, 108, 122]
city2 = [119, 116, 101, 103, 113, 84]
city3 = [70, 68, 54, 73, 81, 68]
city4 = [61, 54, 59, 67, 59, 70]
t = HypothesisTests.KruskalWallisTest(city1, city2, city3, city4)

@test t.n_i == [length(city1), length(city2), length(city3), length(city4)]
@test t.df == 3
@test t.R_i == [104, 113, 53, 30]
@test t.H ≈ 16.028783253379856
@test t.tie_adjustment ≈ 0.9969565217391304
@test pvalue(t) ≈ 0.0011186794961869423
show(IOBuffer(), t)

# example with non-integer rank sum
t1 = [1.2,1.9,2.1]
t2 = [3.4,1.9,5.6]
t3 = [1.3,4.4,9.9]
t = KruskalWallisTest(t1, t2, t3)
@test t.n_i == [length(t1), length(t2), length(t3)]
@test t.df == 2
@test t.R_i ==  [9.5, 17.5, 18.0]
@test t.H ≈ 2.039215686274513
@test t.tie_adjustment ≈ 0.9916666666666667
@test pvalue(t) ≈ 0.3607363776845705
show(IOBuffer(), t)


text:
using SparseArrays

export boundary1
function boundary1(mesh::DECMesh)
    δ = spzeros(length(mesh.tri.edges), length(mesh.tri.vertices))
    edges = mesh.tri.edges
    for i = 1:length(edges)
        e = edges[i]
        δ[i, e[1]] = -1
        δ[i, e[2]] = 1
    end
    return δ
end

export boundary2d
function boundary2d(mesh::DECMesh)
    polyns = mesh.dtri.dvertices
    edges = mesh.dtri.dedges
    δ = spzeros(length(polyns), length(edges))
    for i = 1:length(polyns)
        p = polyns[i]
        e = vcat([[p[i - 1], p[i]] for i = 2:length(p)], [[p[end], p[begin]]])
        for j = 1:length(e)
            if e[j] in edges
                k = findfirst(isequal(e[j]), edges)
                δ[i, k] = 1
            elseif reverse!(e[j]) in edges
                k = findfirst(isequal(e[j]), edges)
                δ[i, k] = -1
            end
        end
    end
    return δ
end

export hodgestar1
function hodgestar1(mesh::DECMesh)
    vertices = mesh.tri.vertices
    edges = mesh.tri.edges
    dtriangles = mesh.dtri.dtriangles
    dedges = mesh.dtri.dedges
    star = spzeros(length(edges), length(edges))
    for i = 1:length(edges)
        edge = edges[i]
        dedge = dedges[i]
        star[i, i] = longitude(dtriangles[dedge]...) / longitude(vertices[edge]...)
    end
    return star
end

export hodgestar2
function hodgestar2(mesh::DECMesh)
    dtriangles = mesh.dtri.dtriangles
    dvertices = mesh.dtri.dvertices
    star = spzeros(length(mesh.tri.vertices), length(mesh.tri.vertices))
    for i = 1:length(mesh.tri.vertices)
        star[i, i] = 1 / area(dtriangles[dvertices[i]])
    end
    return star
end

export laplacian
function laplacian(mesh::DECMesh)
    d₀ = boundary1(mesh)
    d₁ = boundary2d(mesh)
    *₁ = hodgestar1(mesh)
    *₂ = hodgestar2(mesh)

    Δf =  *₂  * d₁ *  *₁  * d₀
end


text:
function count_questions(group_lines)
    counts = zeros(Int, (26))
    n_answers = 1
    for c in group_lines
        int_c = Int(c)
        if int_c == 10
            n_answers += 1
        else
            counts[int_c - 96] += 1
        end
    end
    return sum(counts .== n_answers)
end

function run(s)
    # Your code here
    return sum(count_questions(group_lines) for group_lines in split(s, "\n\n"))
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
module Intersection

export intersect

using IntervalSets

# A, B is disjoint sets
function IntervalSets.intersect(A::AbstractArray{<:AbstractInterval}, B::AbstractArray{<:AbstractInterval})
    # check order
    @assert all(A[i].right ≤ A[i + 1].left for i in 1:(length(A) - 1))
    @assert all(B[i].right ≤ B[i + 1].left for i in 1:(length(B) - 1))

    return filter(I -> !isempty(I), reshape(A, (1, length(A))) .∩ B)
end

# I is an interval, and A is disjoint sets
function IntervalSets.intersect(I::AbstractInterval, A::AbstractArray{<:AbstractInterval})
    # check order
    @assert all(A[i].right ≤ A[i + 1].left for i in 1:(length(A) - 1))

    return filter(I -> !isempty(I), Ref(I) .∩ A)
end
IntervalSets.intersect(A::AbstractArray{<:AbstractInterval}, I::AbstractInterval) =  I ∩ A

end # module Intersection


text:
function toposort(data::Dict{T,Set{T}}) where T
    data = copy(data)
    for (k, v) in data
        delete!(v, k)
    end
    extraitems = setdiff(reduce(∪, values(data)), keys(data))
    for item in extraitems
        data[item] = Set{T}()
    end
    rst = Vector{T}()
    while true
        ordered = Set(item for (item, dep) in data if isempty(dep))
        if isempty(ordered) break end
        append!(rst, ordered)
        data = Dict{T,Set{T}}(item => setdiff(dep, ordered) for (item, dep) in data if item ∉ ordered)
    end
    @assert isempty(data) "a cyclic dependency exists amongst $(keys(data))"
    return rst
end

data = Dict{String,Set{String}}(
    "des_system_lib" => Set(split("std synopsys std_cell_lib des_system_lib dw02 dw01 ramlib ieee")),
    "dw01" =>           Set(split("ieee dw01 dware gtech")),
    "dw02" =>           Set(split("ieee dw02 dware")),
    "dw03" =>           Set(split("std synopsys dware dw03 dw02 dw01 ieee gtech")),
    "dw04" =>           Set(split("dw04 ieee dw01 dware gtech")),
    "dw05" =>           Set(split("dw05 ieee dware")),
    "dw06" =>           Set(split("dw06 ieee dware")),
    "dw07" =>           Set(split("ieee dware")),
    "dware" =>          Set(split("ieee dware")),
    "gtech" =>          Set(split("ieee gtech")),
    "ramlib" =>         Set(split("std ieee")),
    "std_cell_lib" =>   Set(split("ieee std_cell_lib")),
    "synopsys" =>       Set(),
    )

println("# Topologically sorted:\n - ", join(toposort(data), "\n - "))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ForwardDiff, UnPack, LinearAlgebra
using DiffEqBase: DiffCache, get_tmp, dualcache

mutable struct RALΛ{L <: Function, LC}
    Λ::L
    cache::LC
end

function RALΛ(Λ::Function, z::C1, matrix_type::DataType, dims::Tuple{Int, Int}) where {C1 <: AbstractVector{<: Number}}
    cache = matrix_type(undef, 0, 0) # Create empty matrix first, just to check if Λ is in place or not
    if applicable(Λ, cache, z)
        cache = matrix_type(undef, dims)
        Λnew = function _Λ_ip(cache::LCN, z::C1N) where {LCN <: DiffCache, C1N <: AbstractVector{<: Number}}
            Λ(get_tmp(cache, z), z)
            return get_tmp(cache, z)
        end
        return RALΛ(Λnew, dualcache(cache, Val{length(z)}))
    else
        function _Λ_oop(cache::LCN, z::C1N) where {LCN <: Nothing, C1N <: AbstractVector{<: Number}}
            return Λ(z)
        end
        return RALΛ(Λnew, nothing)
    end
end

function RALΛ(Λin::LC, z::C1) where {LC <: AbstractMatrix{<: Number}, C1 <: AbstractVector{<: Number}}
    Λ(cache::LCN, z::C1N) where {LCN <: AbstractMatrix{<: Number}, C1N <: AbstractVector{<: Number}} = cache
    return RALΛ{Function, LC}(Λ, Λin)
end

function (ralλ::RALΛ)(z::C1) where {C1 <: AbstractVector{<: Number}}
    return ralλ.Λ(ralλ.cache, z)
end

mutable struct RALΣ{S <: Function, SC}
    Σ::S
    cache::SC
end

function RALΣ(Σ::Function, z::C1, matrix_type::DataType, dims::Tuple{Int, Int}) where {C1 <: AbstractVector{<: Number}}
    cache = matrix_type(undef, 0, 0)
    if applicable(Σ, cache, z)
        cache = matrix_type(undef, dims)
        Σnew = function _Σ_ip(cache::SCN, z::C1N) where {SCN <: DiffCache, C1N <: AbstractVector{<: Number}}
            du = get_tmp(cache, z)
            Σ(du, z)
            return du
        end
        return RALΣ(Σnew, dualcache(cache, Val{length(z)}))
    else
        Σnew = function _Σ_oop(cache::SCN, z::C1N) where {SCN <: Nothing, C1N <: AbstractVector{<: Number}}
            return Σ(z)
        end
        return RALΣ(Σnew, nothing)
    end
end

function RALΣ(Σin::SC, z::C1) where {SC <: AbstractMatrix{<: Number}, C1 <: AbstractVector{<: Number}}
    Σ(cache::SCN, z::C1N) where {SCN <: AbstractMatrix{<: Number}, C1N <: AbstractVector{<: Number}} = cache
    return RALΣ{Function, SC}(Σ, Σin)
end

function (ralσ::RALΣ)(z::C1) where {C1 <: AbstractVector{<: Number}}
    return ralσ.Σ(ralσ.cache, z)
end

mutable struct RALNonlinearSystem{M <: Function, L <: RALΛ, S <: RALΣ, X <: Function, V <: Function,
                                  VC1 <: AbstractVector{<: Number}, VC2 <: AbstractVector{<: Number}, VC3 <: AbstractVector{<: Number}}
    μ::M         # Functions
    Λ::L         # no type assertion for L b/c it can be Function or Matrix of zeros
    Σ::S         # no type assertion for S b/c it can be Function or constant Matrix
    ξ::X
    𝒱::V
    μ_sss::VC1    # Stochastic steady state values, for caching
    ξ_sss::VC2
    𝒱_sss::VC3
    inplace::NamedTuple{(:μ, :ξ, :𝒱), NTuple{3, Bool}}
end

function RALNonlinearSystem(μ::M, Λ::L, Σ::S, ξ::X, 𝒱::V, μ_sss::VC1, ξ_sss::VC2, 𝒱_sss::VC3,
                            z::C1, y::C1, Ψ::C2, Γ₅::JC5, Γ₆::JC6) where {M <: Function, L <: RALΛ, S <: RALΣ, X <: Function, V <: Function,
                                                                          VC1 <: AbstractVector{<: Number}, VC2 <: AbstractVector{<: Number},
                                                                          VC3 <: AbstractVector{<: Number},
                                                                          C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number},
                                                                          JC5 <: AbstractMatrix{<: Number}, JC6 <: AbstractMatrix{<: Number}}

    inplace = (μ = applicable(μ, μ_sss, z, y), ξ = applicable(ξ, ξ_sss, z, y), 𝒱 = applicable(𝒱, 𝒱_sss, z, Ψ, Γ₅, Γ₆))

    return RALNonlinearSystem{M, L, S, X, V, VC1, VC2, VC3}(μ, Λ, Σ, ξ, 𝒱, μ_sss, ξ_sss, 𝒱_sss, inplace)
end

function update!(m::RALNonlinearSystem, z::C1, y::C1, Ψ::C2,
                 Γ₅::JC5, Γ₆::JC6) where {C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number},
                                          JC5 <: AbstractMatrix{<: Number}, JC6 <: AbstractMatrix{<: Number}}

    if m.inplace[:μ]
        m.μ(m.μ_sss, z, y)
    else
        m.μ_sss .= m.μ(z, y)
    end

    if m.inplace[:ξ]
        m.ξ(m.ξ_sss, z, y)
    else
        m.ξ_sss .= m.ξ(z, y)
    end

    if m.inplace[:𝒱]
        m.𝒱(m.𝒱_sss, z, Ψ, Γ₅, Γ₆)
    else
        m.𝒱_sss .= m.𝒱(z, Ψ, Γ₅, Γ₆)
    end

    m
end

mutable struct RALLinearizedSystem{Mz <: Function, My <: Function, Xz <: Function, Xy <: Function, J <: Function,
                                   JC1 <: AbstractMatrix{<: Number}, JC2 <: AbstractMatrix{<: Number},
                                   JC3 <: AbstractMatrix{<: Number}, JC4 <: AbstractMatrix{<: Number},
                                   JC5 <: AbstractMatrix{<: Number}, JC6 <: AbstractMatrix{<: Number},
                                   JC7 <: AbstractMatrix{<: Number}}
    μz::Mz     # Functions
    μy::My
    ξz::Xz
    ξy::Xy
    J𝒱::J
    Γ₁::JC1    # Jacobians, for caching
    Γ₂::JC2
    Γ₃::JC3
    Γ₄::JC4
    Γ₅::JC5
    Γ₆::JC6
    JV::JC7
    inplace::NamedTuple{(:μz, :μy, :ξz, :ξy, :J𝒱), NTuple{5, Bool}}
end

function RALLinearizedSystem(μz::Mz, μy::My, ξz::Xz, ξy::Xy, J𝒱::J,
                             Γ₁::JC1, Γ₂::JC2, Γ₃::JC3, Γ₄::JC4, Γ₅::JC5, Γ₆::JC6,
                             JV::JC7, z::C1, y::C1, Ψ::C2,
                             μ_sss::VC1, ξ_sss::VC2, 𝒱_sss::VC3) where {Mz <: Function, My <: Function, Xz <: Function,
                                                                        Xy <: Function, J <: Function,
                                                                        JC1 <: AbstractMatrix{<: Number}, JC2 <: AbstractMatrix{<: Number},
                                                                        JC3 <: AbstractMatrix{<: Number}, JC4 <: AbstractMatrix{<: Number},
                                                                        JC5 <: AbstractMatrix{<: Number}, JC6 <: AbstractMatrix{<: Number},
                                                                        JC7 <: AbstractMatrix{<: Number},
                                                                        C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number},
                                                                        VC1 <: AbstractVector{<: Number}, VC2 <: AbstractVector{<: Number},
                                                                        VC3 <: AbstractVector{<: Number},}

    inplace = (μz = applicable(μz, Γ₁, z, y, μ_sss), μy = applicable(μy, Γ₂, z, y, μ_sss), ξz = applicable(ξz, Γ₃, z, y, ξ_sss),
               ξy = applicable(ξy, Γ₄, z, y, ξ_sss), J𝒱 = applicable(J𝒱, JV, z, Ψ, Γ₅, Γ₆, 𝒱_sss))

    return RALLinearizedSystem(μz, μy, ξz, ξy, J𝒱, Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆, JV, inplace)
end

function update!(m::RALLinearizedSystem, z::C1, y::C1, Ψ::C2,
                 μ_sss::VC1, ξ_sss::VC2, 𝒱_sss::VC3) where {C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number},
                                                            VC1 <: AbstractVector{<: Number}, VC2 <: AbstractVector{<: Number},
                                                            VC3 <: AbstractVector{<: Number}}

    if m.inplace[:μz]
        m.μz(m.Γ₁, z, y, μ_sss)
    else
        m.μz(m.Γ₁, z, y)
    end

    if m.inplace[:μy]
        m.μy(m.Γ₂, z, y, μ_sss)
    else
        m.μy(m.Γ₂, z, y)
    end

    if m.inplace[:ξz]
        m.ξz(m.Γ₃, z, y, ξ_sss)
    else
        m.ξz(m.Γ₃, z, y)
    end

    if m.inplace[:ξy]
        m.ξy(m.Γ₄, z, y, ξ_sss)
    else
        m.ξy(m.Γ₄, z, y)
    end

    if m.inplace[:J𝒱]
        m.J𝒱(m.JV, z, Ψ, m.Γ₅, m.Γ₆, 𝒱_sss)
    else
        m.J𝒱(m.JV, z, Ψ, m.Γ₅, m.Γ₆)
    end

    m
end

abstract type AbstractRiskAdjustedLinearization end

"""
    RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, 𝒱, Nz, Ny, Nε)

Creates a first-order perturbation around the stochastic steady state ``(z, y)`` of a discrete-time dynamic model.

(TODO: Move more of the formality to documentation, and make this shorter and concise, w/out explanation of matrix equations)
The affine approximation of the model is
``math
\\begin{aligned}
    \\mathbb{E}[z_{t + 1}] & = \\mu(z, y) + \\Gamma_1(z_t - z) + \\Gamma_2(y_t - y)\\\\
    0                      & = \\xi(z, y) + \\Gamma_3(z_t - z) + \\Gamma_4(y_t - y) + \\Gamma_5 \\mathbb{E}_t z_{t + 1} + \\Gamma_6 \\mathbb{E}_t y_{t + 1} + \\mathscr{V}(z) + J\\mathscr{V}(z)(z_t  - z),
\\end{aligned}
``

where ``\\Gamma_1, \\Gamma_2`` are the Jacobians of ``\\mu`` with respect to ``z_t`` and ``y_t``, respectively;
``\\Gamma_3, \\Gamma_4`` are the Jacobians of ``\\xi`` with respect to ``z_t`` and ``y_t``, respectively;
``\\Gamma_5, \\Gamma_6`` are constant matrices; ``\\mathscr{V}(z)`` is the model's entropy;
``J\\mathscr{V}(z)`` is the Jacobian of the entropy;

and the state variables ``z_t`` and jump variables ``y_t`` follow
``math
\\begin{aligned}
    z_{t + 1} & = z + \\Gamma_1(z_t - z) + \\Gamma_2(y_t - y) + (I_{n_z} - \\Lambda(z_t) \\Psi)^{-1}\\Sigma(z_t)\\varepsilon_{t + 1},\\\\
    y_t       & = y + \\Psi(z_t - z)
\\end{aligned}
``

The unknowns ``(z, y, \\Psi)`` solve the system of equations
``math
\\begin{aligned}
0 & = \\mu(z, y) - z,\\\\
0 & = \\xi(z, y) + \\Gamma_5 z + \\Gamma_6 y + \\mathscr{V}(z),\\\\
0 & = \\Gamma_3 + \\Gamma_4 \\Psi + (\\Gamma_5 + \\Gamma_6 \\Psi)(\\Gamma_1 + \\Gamma_2 \\Psi) + J\\mathscr{V}(z).
\\end{aligned}
``
(TODO: Move the nonlinear model statement to documentation)
The true nonlinear equations defining model are assumed to take the form

``math
\\begin{aligned}
    z_{t + 1} & = \\mu(z_t, y_t) + \\Lambda(z_t)(y_{t + 1} - \\mathbb{E}_t y_{t + 1}) + \\Sigma(z_t) \\varepsilon_{t + 1},\\\\
    0 & = \\log\\mathbb{E}_t[\\exp(\\xi(z_t, y_t) + \\Gamma_5 z_{t + 1} + \\Gamma_6 y_{t + 1})].
\\end{aligned}
``

The vectors ``z_t\\in \\mathbb{R}^{n_z}`` and ``y_t \\in \\mathbb{R}^{n_y}`` are the state and jump variables, respectively.
The first vector equation comprise the model's expectational equations, which are typically
the first-order conditions for the jump variables from agents' optimization problem.
The second vector equation comprise the transition equations of the state variables. The exogenous shocks
``\\varepsilon\\in\\mathbb{R}^{n_\\varepsilon}`` form a martingale difference sequence whose distribution
is described by the differentiable, conditional cumulant generating function (ccgf)

``math
\\begin{aligned}
\\kappa[\\alpha(z_t) \\mid z_t] = \\log\\mathbb{E}_t[\\exp(\\alpha(z_t)' \\varepsilon_{t + 1})],\\quad \text{ for any differentiable map }\\alpha::\\mathbb{R}^{n_z}\\rightarrow\\mathbb{R}^{n_\\varepsilon}.
\\end{aligned}
``

The functions
``math
\\begin{aligned}
\\xi:\\mathbb{R}^{2n_y + 2n_z}\\rightarrow \\mathbb{R}^{n_y},& \\quad \\mu:\\mathbb{R}^{n_y + n_z}\\rightarrow \\mathbb{R}^{n_z},\\\\
\\Lambda::\\mathbb{R}^{n_z} \\rightarrow \\mathbb{R}^{n_z \\times n_y}, & \\quad \\Sigma::\\mathbb{R}^{n_z}\\rightarrow \\mathbb{R}^{n_z\\times n_\\varepsilon}
\\end{aligned}
are differentiable. The first two functions characterize the effects of time ``t`` variables on the expectational and
state transition equations. The function ``\\Lambda`` characterizes heteroskedastic endogenous risk that depends on
innovations in jump variables while the function ``\\Sigma`` characterizes exogenous risk.

Refer to Lopz et al. (2018) "Risk-Adjusted Linearizations of Dynamic Equilibrium Models" for details.
"""
mutable struct RiskAdjustedLinearization{A <: RALNonlinearSystem, B <: RALLinearizedSystem,
                                         C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number}} <: AbstractRiskAdjustedLinearization
    nonlinear::A
    linearization::B
    z::C1        # Coefficients
    y::C1
    Ψ::C2
    Nz::Int      # Dimensions
    Ny::Int
    Nε::Int
end
# TODO
# 1.UPDATE THE PRINTING, maybe just write out "risk-adjusted linearization with dimensions ()"
#
# 2. Test update! functions for the various blocks as well as access functions for RiskAdjustedLinearization
#
# 3. Check inplace inference is correct, check construction of each block plus main block
#=
TODO: Finish this once the final struct is completed
# A series of lower level constructors
function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, 𝒱::V, μz::Mz, μy::My, ξz::Xz, ξy::Xy, J𝒱::J,
                                   μ_sss::AbstractVector{T}, ξ_sss::AbstractVector{T}, 𝒱_sss::AbstractVector{T},
                                   Γ₁::AbstractMatrix{T}, Γ₂::AbstractMatrix{T}, Γ₃::AbstractMatrix{T}
                                   Γ₄::AbstractMatrix{T}, Γ₅::AbstractMatrix{T}, Γ₆::AbstractMatrix{T},
                                   JV::AbstractMatrix{T}, z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1) where {T <: Number, M <: Function, L,
                                                        S, X <: Function, V <: Function,
                                                        Mz <: Function, My <: Function, Xz <: Function,
                                                        Xy <: Function, J <: Function}

    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        Nε = size(Σ(z), 2)
    end

    return RiskAdjustedLinearization{T, M, L, S, X, V, J}(μ, Λ, Σ, ξ, 𝒱, J𝒱, μ_sss, ξ_sss, 𝒱_sss,
                                                          Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆,
                                                          JV, z, y, Ψ, Nz, Ny, Nε)
end

function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, 𝒱::V, μz::Mz, μy::My, ξz::Xz, ξy::Xy, J𝒱::J,
                                   Γ₁::AbstractMatrix{T}, Γ₂::AbstractMatrix{T}, Γ₃::AbstractMatrix{T}
                                   Γ₄::AbstractMatrix{T}, Γ₅::AbstractMatrix{T}, Γ₆::AbstractMatrix{T},
                                   JV::AbstractMatrix{T}, z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1) where {T <: Number, M <: Function, L,
                                                        S, X <: Function, V <: Function,
                                                        Mz <: Function, My <: Function, Xz <: Function,
                                                        Xy <: Function, J <: Function}
    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        Nε = size(Σ(z), 2)
    end

    # Cache stochastic steady state vectors
    μ_sss, ξ_sss, 𝒱_sss = _cache_sss_vectors(z, y)

    return RiskAdjustedLinearization{T, M, L, S, X, V, J}(μ, Λ, Σ, ξ, 𝒱, J𝒱, μ_sss, ξ_sss, 𝒱_sss,
                                                          Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆,
                                                          JV, z, y, Ψ, Nz, Ny, Nε)
end

function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, 𝒱::V, μz::Mz, μy::My, ξz::Xz, ξy::Xy, J𝒱::J,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1) where {T <: Number, M <: Function, L,
                                                        S, X <: Function, V <: Function,
                                                        Mz <: Function, My <: Function, Xz <: Function,
                                                        Xy <: Function, J <: Function}
    # Get dimensions
    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        Nε = size(Σ(z), 2)
    end

    # Cache stochastic steady state vectors
    μ_sss, ξ_sss, 𝒱_sss = _cache_sss_vectors(z, y)

    # Cache stochastic steady state Jacobians
    Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆, JV = _cache_jacobians(Ψ, Nz, Ny)

    return RiskAdjustedLinearization{T, M, L, S, X, V, J}(μ, Λ, Σ, ξ, 𝒱, J𝒱, μ_sss, ξ_sss, 𝒱_sss,
                                                          Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆,
                                                          JV, z, y, Ψ, Nz, Ny, Nε)
end
=#
function RiskAdjustedLinearization(nonlinear::A, linearization::B, z::C1, y::C1, Ψ::C2,
                                   Nz::Int, Ny::Int, Nε::Int;
                                   check_inputs::Bool = true) where {A <: RALNonlinearSystem, B <: RALLinearizedSystem,
                                                                     C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number}}

    # Make sure inputs are well-formed
    if check_inputs
        _check_inputs(nonlinear, linearization, z, y, Ψ)
    end

    return RiskAdjustedLinearization{A, B, C1, C2}(nonlinear, linearization, z, y, Ψ, Nz, Ny, Nε)
end

# Constructor that uses ForwardDiff to calculate Jacobian functions
# NOTE THAT here we pass in the ccgf, rather than 𝒱
function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, Γ₅::JC5, Γ₆::JC6, ccgf::CF,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nz::Int, Ny::Int, Nε::Int; sss_vector_type::DataType = Vector{T},
                                   jacobian_type::DataType = Matrix{T}) where {T <: Number, M <: Function, L <: RALΛ, S <: RALΣ,
                                                                               X <: Function,
                                                                               JC5 <: AbstractMatrix{<: Number},
                                                                               JC6 <: AbstractMatrix{<: Number},
                                                                               CF <: Function}

    # Cache stochastic steady state vectors
    μ_sss, ξ_sss, 𝒱_sss = _cache_sss_vectors(z, y)

    # Cache stochastic steady state Jacobians
    Γ₁, Γ₂, Γ₃, Γ₄, JV = _cache_jacobians(Ψ, Nz, Ny, jacobian_type)

    # Use cached Jacobians to create Jacobian functions for μ, ξ
    if applicable(μ, z, y) # Check if μ is in place or not
        μz = (F, z, y) -> ForwardDiff.jacobian!(F, x -> μ(x, y), z) # not in place
        μy = (F, z, y) -> ForwardDiff.jacobian!(F, x -> μ(z, x), y)
    else # in place
        μz = (F, z, y, μ_sss) -> ForwardDiff.jacobian!(F, (G, x) -> μ(G, x, y), μ_sss, z)
        μy = (F, z, y, μ_sss) -> ForwardDiff.jacobian!(F, (G, x) -> μ(G, z, x), μ_sss, y)
    end

    if applicable(ξ, z, y) # Check if ξ is in place or not
        ξz = (F, z, y) -> ForwardDiff.jacobian!(F, x -> ξ(x, y), z) # not in place
        ξy = (F, z, y) -> ForwardDiff.jacobian!(F, x -> ξ(z, x), y)
    else # in place
        ξz = (F, z, y, ξ_sss) -> ForwardDiff.jacobian!(F, (G, x) -> ξ(G, x, y), ξ_sss, z)
        ξy = (F, z, y, ξ_sss) -> ForwardDiff.jacobian!(F, (G, x) -> ξ(G, z, x), ξ_sss, y)
    end

    # Create 𝒱 and its Jacobian J𝒱
    if applicable(ccgf, Γ₅, z) # Check if ccgf is in place or not
        𝒱 = function _𝒱(F, z, Ψ, Γ₅, Γ₆)
            F .= ccgf((Γ₅ + Γ₆ * Ψ) * ((I - Λ(z) * Ψ) \ Σ(z)), z)
        end
    else # in place
        𝒱 = (F, z, Ψ, Γ₅, Γ₆) -> ccgf(F, (Γ₅ + Γ₆ * Ψ) * ((I - Λ(z) * Ψ) \ Σ(z)), z)
    end
    J𝒱 = function _J𝒱(F, z, Ψ, Γ₅, Γ₆, 𝒱_sss)
        ForwardDiff.jacobian!(F, (G, x) -> 𝒱(G, x, Ψ, Γ₅, Γ₆), 𝒱_sss, z)
    end

    # Form underlying RAL blocks
    nonlinear_system  = RALNonlinearSystem(μ, Λ, Σ, ξ, 𝒱, μ_sss, ξ_sss, 𝒱_sss, z, y, Ψ, Γ₅, Γ₆)
    linearized_system = RALLinearizedSystem(μz, μy, ξz, ξy, J𝒱, Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆, JV, z, y, Ψ, μ_sss, ξ_sss, 𝒱_sss)

    return RiskAdjustedLinearization(nonlinear_system, linearized_system, z, y, Ψ, Nz, Ny, Nε)
end

function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, Γ₅::JC5, Γ₆::JC6, ccgf::CF,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int; sss_vector_type::DataType = Vector{T}, sss_matrix_type::DataType = Matrix{T},
                                   jacobian_type::DataType = Matrix{T}) where {T <: Number, M <: Function, L <: Function, S <: Function,
                                                                               X <: Function,
                                                                               JC5 <: AbstractMatrix{<: Number},
                                                                               JC6 <: AbstractMatrix{<: Number},
                                                                               CF <: Function}
    # Get dimensions
    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        error("Nε cannot be negative")
    end

    # Create wrappers enabling caching for Λ and Σ
    Λ = RALΛ(Λ, z, sss_matrix_type, (Nz, Ny))
    Σ = RALΣ(Σ, z, sss_matrix_type, (Nz, Nε))

    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nz, Ny, Nε, sss_vector_type = sss_vector_type,
                                     jacobian_type = jacobian_type)
end

function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, Γ₅::JC5, Γ₆::JC6, ccgf::CF,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1; sss_vector_type::DataType = Vector{T}, sss_matrix_type::DataType = Matrix{T},
                                   jacobian_type::DataType = Matrix{T}) where {T <: Number, M <: Function, L <: AbstractMatrix{<: Number}, S <: Function,
                                                                               X <: Function,
                                                                               JC5 <: AbstractMatrix{<: Number},
                                                                               JC6 <: AbstractMatrix{<: Number},
                                                                               CF <: Function}
    # Get dimensions
    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        error("Nε cannot be negative")
    end

    # Create wrappers enabling caching for Λ and Σ
    Λ = RALΛ(Λ, z)
    Σ = RALΣ(Σ, z, sss_matrix_type, (Nz, Nε))

    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nz, Ny, Nε, sss_vector_type = sss_vector_type,
                                     jacobian_type = jacobian_type)
end

function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, Γ₅::JC5, Γ₆::JC6, ccgf::CF,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1; sss_vector_type::DataType = Vector{T}, sss_matrix_type::DataType = Matrix{T},
                                   jacobian_type::DataType = Matrix{T}) where {T <: Number, M <: Function, L <: Function, S <: AbstractMatrix{<: Number},
                                                                               X <: Function,
                                                                               JC5 <: AbstractMatrix{<: Number},
                                                                               JC6 <: AbstractMatrix{<: Number},
                                                                               CF <: Function}
    # Get dimensions
    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        error("Nε cannot be negative")
    end

    # Create wrappers enabling caching for Λ and Σ
    Λ = RALΛ(Λ, z, sss_matrix_type, (Nz, Ny))
    Σ = RALΣ(Σ, z)

    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nz, Ny, Nε, sss_vector_type = sss_vector_type,
                                     jacobian_type = jacobian_type)
end

function RiskAdjustedLinearization(μ::M, Λ::L, Σ::S, ξ::X, Γ₅::JC5, Γ₆::JC6, ccgf::CF,
                                   z::AbstractVector{T}, y::AbstractVector{T}, Ψ::AbstractMatrix{T},
                                   Nε::Int = -1; sss_vector_type::DataType = Vector{T}, sss_matrix_type::DataType = Matrix{T},
                                   jacobian_type::DataType = Matrix{T}) where {T <: Number, M <: Function,
                                                                               L <: AbstractMatrix{<: Number}, S <: AbstractMatrix{<: Number},
                                                                               X <: Function,
                                                                               JC5 <: AbstractMatrix{<: Number},
                                                                               JC6 <: AbstractMatrix{<: Number},
                                                                               CF <: Function}
    # Get dimensions
    Nz = length(z)
    Ny = length(y)
    if Nε < 0
        error("Nε cannot be negative")
    end

    # Create wrappers enabling caching for Λ and Σ
    Λ = RALΛ(Λ, z)
    Σ = RALΣ(Σ, z)

    return RiskAdjustedLinearization(μ, Λ, Σ, ξ, Γ₅, Γ₆, ccgf, z, y, Ψ, Nz, Ny, Nε, sss_vector_type = sss_vector_type,
                                     jacobian_type = jacobian_type)
end

function _cache_jacobians(Ψ::AbstractMatrix{T}, Nz::Int, Ny::Int, mat_type::DataType) where {T <: Number}

    Γ₁ = mat_type(undef, Nz, Nz)
    Γ₂ = mat_type(undef, Nz, Ny)
    Γ₃ = similar(Ψ)
    Γ₄ = mat_type(undef, Ny, Ny)
    JV = similar(Ψ)

    return Γ₁, Γ₂, Γ₃, Γ₄, JV
end

function _cache_sss_vectors(z::AbstractVector{T}, y::AbstractVector{T}) where {T <: Number, L, S}

    μ_sss = similar(z)
    ξ_sss = similar(y)
    𝒱_sss = similar(y)

   return μ_sss, ξ_sss, 𝒱_sss
end

function _check_inputs(nonlinear::A, linearization::B, z::C1, y::C1, Ψ::C2) where {A <: RALNonlinearSystem, B <: RALLinearizedSystem,
                                                                                C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number}}

    # Get contents of nonlinear and linearization blocks
    @unpack μ, ξ, 𝒱, μ_sss, ξ_sss, 𝒱_sss = nonlinear
    @unpack μz, μy, ξz, ξy, J𝒱, Γ₁, Γ₂, Γ₃, Γ₄, Γ₅, Γ₆, JV = linearization

    @assert applicable(μ, z, y) ||
        applicable(μ, μ_sss, z, y) "The function μ must take either the form " *
        "μ(z, y) or the in-place equivalent μ(F, z, y)"

    @assert applicable(ξ, z, y) ||
        applicable(ξ, ξ_sss, z, y) "The function μ must take either the form " *
        "ξ(z, y) or the in-place equivalent ξ(F, z, y)"

    @assert applicable(𝒱, z, Ψ, Γ₅, Γ₆) ||
        applicable(𝒱, y, z, Ψ, Γ₅, Γ₆) "The function 𝒱 must take either the form " *
        "𝒱(z, Ψ, Γ₅, Γ₆) or the in-place equivalent 𝒱(F, z, Ψ, Γ₅, Γ₆)"

    @assert applicable(μz, Γ₁, z, y) ||
        applicable(μz, Γ₁, z, y, μ_sss) "The function μz must take either the form " *
        "μz(F, z, y) or μz(F, z, y, μ_sss)"

    @assert applicable(μy, Γ₂, z, y) ||
        applicable(μy, Γ₂, z, y, μ_sss) "The function μy must take either the form " *
        "μy(F, z, y) or μy(F, z, y, μ_sss)"

    @assert applicable(ξz, Γ₃, z, y) ||
        applicable(ξz, Γ₃, z, y, ξ_sss) "The function ξz must take either the form " *
        "ξz(F, z, y) or ξz(F, z, y, ξ_sss)"

    @assert applicable(ξy, Γ₄, z, y) ||
        applicable(ξy, Γ₄, z, y, ξ_sss) "The function ξy must take either the form " *
        "ξy(F, z, y) or ξy(F, z, y, ξ_sss)"

    @assert applicable(J𝒱, z, Ψ, Γ₅, Γ₆) ||
        applicable(J𝒱, JV, z, Ψ, Γ₅, Γ₆, 𝒱_sss) "The function J𝒱 must take either the form " *
        "J𝒱(F, z, Ψ, Γ₅, Γ₆) or J𝒱(F, z, Ψ, Γ₅, Γ₆, 𝒱_sss)"
end

## Methods for using RiskAdjustedLinearization
@inline Γ₁(m::RiskAdjustedLinearization) = m.linearization.Γ₁
@inline Γ₂(m::RiskAdjustedLinearization) = m.linearization.Γ₂
@inline Γ₃(m::RiskAdjustedLinearization) = m.linearization.Γ₃
@inline Γ₄(m::RiskAdjustedLinearization) = m.linearization.Γ₄
@inline Γ₅(m::RiskAdjustedLinearization) = m.linearization.Γ₅
@inline Γ₆(m::RiskAdjustedLinearization) = m.linearization.Γ₆
@inline JV(m::RiskAdjustedLinearization) = m.linearization.JV
@inline getvalues(m::RiskAdjustedLinearization) = (m.z, m.y, m.Ψ)
@inline getvecvalues(m::RiskAdjustedLinearization) = vcat(m.z, m.y, vec(m.Ψ))
@inline nonlinear_system(m::RiskAdjustedLinearization) = m.nonlinear
@inline linearized_system(m::RiskAdjustedLinearization) = m.linearization

function update!(m::RiskAdjustedLinearization, z::C1, y::C1, Ψ::C2;
                 update_cache::Bool = true) where {C1 <: AbstractVector{<: Number}, C2 <: AbstractMatrix{<: Number}}

    # Update values of the affine approximation
    m.z .= z
    m.y .= y
    m.Ψ .= Ψ

    # Update the cached vectors and Jacobians
    if update_cache
        update!(nonlinear_system(m), m.z, m.y, m.Ψ, Γ₅(m), Γ₆(m))

        update!(linearized_system(m), m.z, m.y, m.Ψ, m.nonlinear.μ_sss, m.nonlinear.ξ_sss, m.nonlinear.𝒱_sss)
    end

    m
end

function Base.show(io::IO, m::RiskAdjustedLinearization)
    @printf io "Risk-Adjusted Linearization of a Dynamic Economic Model\n"
    @printf io "No. of state variables:              %i\n" m.Nz
    @printf io "No. of jump variables:               %i\n" m.Ny
    @printf io "No. of exogenous shocks:             %i\n" m.Nε
end


text:
module TestPath

using Modia
using Modia.Measurements
@usingModiaPlot

const ptp_path = PTP_path(["angle1", "angle2", "angle3"],
                          positions = [0.0 2.0 3.0;  # angle1=0.0, angle2=2.0, angle3=3.0
                                       0.5 3.0 4.0;
                                       0.8 1.5 0.3;
                                       0.2 1.5 0.8],
                          startTime = (0.1/60)u"minute",
                          v_max = 2*ones(3),
                          a_max = 3*ones(3))
angles = zeros(3)
getPosition!(ptp_path, 0.5, angles)   # angles = [0.12, 2.24, 3.24]
path = getPath(ptp_path)
printResultInfo(path)

plot(path, [("angle1", "angle2", "angle3"),
            ("der(angle1)", "der(angle2)", "der(angle3)"),
            ("der2(angle1)", "der2(angle2)", "der2(angle3)")], figure=1)

plotPath(ptp_path,plot,onlyPositions=false, heading="ptp_path", figure=2)

nom(v) = measurement(v)
err(v) = v ± (0.1*v)
const ptp_path2 = PTP_path{Measurement{Float64}}(
                          ["angle1", "angle2", "angle3"],
                           positions = [nom(0.0) nom(2.0) nom(3.0);  # angle1=0.0, angle2=2.0, angle3=3.0
                                        err(0.5) err(3.0) err(4.0);
                                        err(0.8) err(1.5) err(0.3);
                                        err(0.2) err(1.5) err(0.8)],
                           startTime = nom(0.1),
                           v_max = fill(measurement(2.0),3),
                           a_max = fill(measurement(3.0),3))
angles = zeros(Measurement{Float64}, 3)
getPosition!(ptp_path2, nom(0.5), angles)   # angles = [0.12, 2.24, 3.24]
path2 = getPath(ptp_path2)
printResultInfo(path2)

plot(path2, [("angle1", "angle2", "angle3"),
             ("der(angle1)", "der(angle2)", "der(angle3)"),
             ("der2(angle1)", "der2(angle2)", "der2(angle3)")], figure=3)

plotPath(ptp_path2,plot,onlyPositions=false, heading="ptp_path2", figure=4)
end


text:
using JLD2, Test

fn = joinpath(mktempdir(), "test.jld")
f = jldopen(fn, "w")
write(f, "/test_group_1/x1", 1)
@test !isempty(read(f, "test_group_1"))
write(f, "test_group_1/x2", 2)
write(f, "/test_group_2/x1", 3)
write(f, "test_group_3/x1", 4)
read(f, "test_group_3")["contained_group/x1"] = 5
JLD2.Group(f, "test_group_4")["x1"] = 6
JLD2.Group(read(f, "test_group_4"), "contained_group")["x1"] = 7
@test isempty(JLD2.Group(f, "empty_group"))
JLD2.Group(f, "group_with_one_group_child/empty_group")
@test !isempty(read(f, "group_with_one_group_child"))
@test_throws ArgumentError JLD2.Group(f, "test_group_4")
@test_throws ArgumentError write(f, "test_group_4/contained_group", 0)
@test !haskey(f, "path_that_does_not_exist")
@test get(f, "path_that_does_not_exist", 2.0) == 2.0
@test get!(f, "path_created_by_get!", "x") == "x"
@test f["path_created_by_get!"] == "x"
@test get(JLD2.Group(f, "test_group_5"), "path_that_does_not_exist", 3.0) == 3.0
@test get!(f["test_group_5"], "path_created_by_get!", 4.0) == 4.0
@test f["test_group_5"]["path_created_by_get!"] == 4.0

# Test that the default element constructor is called only when necessary
# First, we'll access a group which *does* exist, so the default constructor is not called:
@test get(() -> error("this should not be called because the key exists"), f["test_group_5"], "path_created_by_get!") == 4.0

# Then we'll access a few keys that do not exist, and verify that the default
# constructor was called:
default_called = false
get(f, "path_that_does_not_exist") do
    global default_called
    default_called = true
end
@test default_called

default_called = false
get!(JLD2.Group(f, "test_group_6"), "path_created_by_get!") do
    global default_called
    default_called = true
end
@test default_called
close(f)

f = jldopen(fn, "r")
@test read(f, "/test_group_1/x1") == 1
@test read(f, "/test_group_1/x2") == 2
@test read(f, "/test_group_2/x1") == 3
@test read(f, "/test_group_3/x1") == 4
@test read(f, "/test_group_3/contained_group/x1") == 5
@test_throws KeyError read(f, "/test_group_3/contained_group/x2")
@test_throws KeyError read(f, "/test_group_3/contained_group/x2/oops")

@test read(f, "test_group_1/x1") == 1
@test read(f, "test_group_1/x2") == 2
@test read(f, "test_group_2/x1") == 3
@test read(f, "test_group_3/x1") == 4
@test read(f, "test_group_3/contained_group/x1") == 5

@test read(f, "test_group_1")["x1"]  == 1
@test read(f, "test_group_1")["x2"]  == 2
@test read(f, "test_group_2")["x1"]  == 3
@test read(f, "test_group_3")["x1"]  == 4
@test read(f, "test_group_3")["contained_group"]["x1"]  == 5

@test read(f, "test_group_1")["/test_group_2/x1"]  == 3
@test read(f, "test_group_3/contained_group")["x1"]  == 5
@test read(f, "test_group_3")["contained_group/x1"]  == 5

@test !isempty(read(f, "group_with_one_group_child"))
@test !isempty(read(f, "test_group_1"))
@test isempty(read(f, "empty_group"))
@test isempty(read(f, "group_with_one_group_child/empty_group"))

@test get(f, "path_that_does_not_exist", 1.0) == 1.0
@test get(f["test_group_5"], "path_that_does_not_exist", 1.0) == 1.0
@test f["test_group_5"]["path_created_by_get!"] == 4.0

# Throws an ArgumentError because the element doesn't exist and cannot
# be created for a read-only file:
@test_throws ArgumentError get!(f, "path_that_does_not_exist", 1.0)
@test_throws ArgumentError get!(f["test_group_5"], "path_that_does_not_exist", 1.0)
# Does not throw an ArgumentError because the element does exist:
@test get!(f["test_group_1"], "x1", 2) == 1

# Make sure printing doesn't error
show(IOBuffer(), f)
close(f)


text:
# Here starts the MCPhyloTree module
module MCPhyloTree

# load necessary packages
using Statistics
using DataStructures
using RecipesBase
using Random

### Get Files

# Type
include("Node.jl")

# Basic Functionality
include("Basics/Tree_Basics.jl")
include("Basics/Tree_Traversal.jl")
include("Basics/Tree_Pruning.jl")
include("Basics/Tree_Search.jl")
include("Basics/Represenations.jl")

# Tree Building
include("Building/Tree2Matrix.jl")
include("Building/Tree_Building.jl")
include("Building/Tree_Clustering.jl")
include("Building/Tree_Consensus.jl")
include("Building/Tree_Ladderizing.jl")
include("Building/ParseNewick.jl")

# Tree Distance
include("Distance/Tree_Distance.jl")

# Tree Moves
include("Moves/NNI.jl")
include("Moves/SPR.jl")
include("Moves/EdgeLength.jl")
include("Moves/Randomize.jl")
include("Moves/Rerooting.jl")

# utils
include("utils.jl")

# plotting recipes
include("Plotting/tree_plot.jl")
include("Plotting/ascii.jl")

# Export Type
export AbstractNode, GeneralNode, Node

export add_child!,
    delete_node!,
    remove_child!,
    insert_node!,
    tree_height,
    tree_length,
    node_height,
    node_age,
    node_depth,
    node_distance,
    path_length,
    set_binary!,
    random_node,
    set_branchlength_vector!,
    find_lca,
    check_binary,
    check_leafsets,
    prune_tree,
    prune_tree!,
    find_num,
    find_binary,
    find_root,
    post_order,
    pre_order,
    level_order,
    newick,
    find_by_name,
    print_ascii

# getters
export get_branchlength_vector,
    get_leaves, get_mother, get_path, get_sister, get_bipartitions

# Building
export to_covariance,
    to_covariance_ultra,
    to_df,
    from_df,
    leave_incidence_matrix,
    to_distance_matrix,
    ParseNewick,
    parsing_newick_string,
    upgma,
    neighbor_joining,
    majority_consensus_tree,
    loose_consensus_tree,
    greedy_consensus_tree,
    ladderize_tree,
    ladderize_tree!,
    create_tree_from_leaves,
    cov2tree

# Distance
export RF, RF_weighted, BHV_bounds

# Moves
export slide,
    slide!,
    swing,
    swing!,
    change_edge_length!,
    NNI,
    NNI!,
    randomize,
    randomize!,
    reroot,
    SPR,
    SPR!

end


text:
# Sectors corresponding to irreducible representations of compact groups
#==============================================================================#
# Irreps of Abelian groups
#------------------------------------------------------------------------------#
abstract type AbelianIrrep <: Sector end

Base.@pure fusiontype(::Type{<:AbelianIrrep}) = Abelian
Base.@pure braidingtype(::Type{<:AbelianIrrep}) = Bosonic

Nsymbol(a::G, b::G, c::G) where {G<:AbelianIrrep} = c == first(a ⊗ b)
Fsymbol(a::G, b::G, c::G, d::G, e::G, f::G) where {G<:AbelianIrrep} =
    Int(Nsymbol(a,b,e)*Nsymbol(e,c,d)*Nsymbol(b,c,f)*Nsymbol(a,f,d))
frobeniusschur(a::AbelianIrrep) = 1
Bsymbol(a::G, b::G, c::G) where {G<:AbelianIrrep} = Int(Nsymbol(a, b, c))
Rsymbol(a::G, b::G, c::G) where {G<:AbelianIrrep} = Int(Nsymbol(a, b, c))

# ZNIrrep: irreps of Z_N are labelled by integers mod N; do we ever want N > 127?
struct ZNIrrep{N} <: AbelianIrrep
    n::Int8
    function ZNIrrep{N}(n::Integer) where {N}
        new{N}(mod(n, N))
    end
end
Base.one(::Type{ZNIrrep{N}}) where {N} =ZNIrrep{N}(0)
Base.conj(c::ZNIrrep{N}) where {N} = ZNIrrep{N}(-c.n)
⊗(c1::ZNIrrep{N}, c2::ZNIrrep{N}) where {N} = (ZNIrrep{N}(c1.charge+c2.charge),)

Base.convert(Z::Type{<:ZNIrrep}, n::Real) = Z(convert(Int, n))

const ℤ₂ = ZNIrrep{2}
const ℤ₃ = ZNIrrep{3}
const ℤ₄ = ZNIrrep{4}
const Parity = ZNIrrep{2}
Base.show(io::IO, ::Type{ZNIrrep{2}}) = print(io, "ℤ₂")
Base.show(io::IO, ::Type{ZNIrrep{3}}) = print(io, "ℤ₃")
Base.show(io::IO, ::Type{ZNIrrep{4}}) = print(io, "ℤ₄")
Base.show(io::IO, ::Type{ZNIrrep{4}}) = print(io, "ℤ₄")
Base.show(io::IO, c::ZNIrrep{2}) = get(io, :compact, false) ? print(io, c.n) : print(io, "ℤ₂(", c.n, ")")
Base.show(io::IO, c::ZNIrrep{3}) = get(io, :compact, false) ? print(io, c.n) : print(io, "ℤ₃(", c.n, ")")
Base.show(io::IO, c::ZNIrrep{4}) = get(io, :compact, false) ? print(io, c.n) : print(io, "ℤ₄(", c.n, ")")
Base.show(io::IO, c::ZNIrrep{N}) where {N} = get(io, :compact, false) ? print(io, c.n) :
    print(io, "ZNIrrep{", N, "}(" , c.n, ")")

# U1Charge: irreps of U1 are labelled by integers
struct U1Irrep <: AbelianIrrep
    charge::Int
end
Base.one(::Type{U1Irrep}) = U1Irrep(0)
Base.conj(c::U1Irrep) = U1Irrep(-c.charge)
⊗(c1::U1Irrep, c2::U1Irrep) = (U1Irrep(c1.charge+c2.charge),)

Base.convert(::Type{U1Irrep}, c::Real) = U1Irrep(convert(Int, c))

const U₁ = U1Irrep
Base.show(io::IO, ::Type{U1Irrep}) = print(io, "U₁")
Base.show(io::IO, c::U1Irrep) = get(io, :compact, false) ? print(io, c.charge) : print(io, "U₁(", c.charge, ")")

# NOTE: FractionalU1Charge?

# Nob-abelian groups
#------------------------------------------------------------------------------#
# SU2Charge: irreps of SU2 are labelled by half integers j, internally we use the integer dimension 2j+1 instead
import WignerSymbols

struct SU2IrrepException <: Exception end
Base.show(io::IO, ::SU2IrrepException) = print(io, "Irreps of (bosonic or fermionic) `SU₂` should be labelled by non-negative half integers, i.e. elements of `Rational{Int}` with denominator 1 or 2")

struct SU2Irrep <: Sector
    dim::Int
    # Let constructor take the actual half integer value j
    SU2Irrep(j::Int) = j >= 0 ? new(2*j+1) : throw(SU2IrrepException)
    function SU2Irrep(j::Rational{Int})
        if j.den == 2
            new(j.num+1)
        elseif j.den == 1
            new(2*j.num+1)
        else
            throw(SU2IrrepException)
        end
    end
end
_getj(s::SU2Irrep) = (s.dim-1)//2

Base.one(::Type{SU2Irrep}) = SU2Irrep(0)
Base.conj(s::SU2Irrep) = s
⊗(s1::SU2Irrep, s2::SU2Irrep) = SectorSet{SU2Irrep}( abs(_getj(s1)-_getj(s2)):(_getj(s1)+_getj(s2)) )

Base.convert(::Type{SU2Irrep}, s::Real) = SU2Irrep(convert(Int, 2*s)//2)

dim(s::SU2Irrep) = s.dim

Base.@pure fusiontype(::Type{SU2Irrep}) = SimpleNonAbelian
Base.@pure braidingtype(::Type{SU2Irrep}) = Bosonic

Nsymbol(sa::SU2Irrep, sb::SU2Irrep, sc::SU2Irrep) = WignerSymbols.δ(_getj(sa), _getj(sb), _getj(sc))
Fsymbol(s1::SU2Irrep, s2::SU2Irrep, s3::SU2Irrep, s4::SU2Irrep, s5::SU2Irrep, s6::SU2Irrep) =
    WignerSymbols.racahW(map(_getj,(s1,s2,s4,s3,s5,s6))...)*sqrt(dim(s5)*dim(s6))
function Rsymbol(sa::SU2Irrep, sb::SU2Irrep, sc::SU2Irrep)
    Nsymbol(sa, sb, sc) || return 0
    iseven(convert(Int, _getj(sa)+_getj(sb)-_getj(sc))) ? 1.0 : -1.0
end

const SU₂ = SU2Irrep
Base.show(io::IO, ::Type{SU2Irrep}) = print(io, "SU₂")
Base.show(io::IO, s::SU2Irrep) = get(io, :compact, false) ? print(io, _getj(s)) : print(io, "SU₂(", _getj(s), ")")


text:
include("bson_tests.jl")
include("mongodb_tests.jl")


text:
using Pkg

Pkg.activate(".")
Pkg.test("Fjage")


text:
#export Simplex

# U: the dimension of the universe
# D: the dimension of the manifold
# N: the number of vertices
# T: the type of the coordinates
# C: the complimentary dimension (should always be U-D)
struct Simplex{U,D,C,N,T}
    vertices::SVector{N,SVector{U,T}}
    tangents::SVector{D,SVector{U,T}}
    normals::SVector{C,SVector{U,T}}
    volume::T
end

"""
    coordtype(simplex)

Return coordinate type used by simplex.
"""
coordtype(::Type{Simplex{U,D,C,N,T}}) where {U,D,C,N,T} = T
coordtype(p::Simplex) = coordtype(typeof(p))

"""
    volume(simplex)

Return the volume of the simplex.
"""
volume(p::Simplex) = p.volume

"""
A tuple of points, aka an interval behaves trivially like a chart
"""
volume(x::Tuple{T,T}) where {T} = norm(x[2]-x[1])

"""
    dimension(simplex)

Return the manifold dimension of the simplex.
"""
dimension(::Type{Simplex{U,D,C,N,T}}) where {U,D,C,N,T} = D

"""
    dimension(simplex)

Return the manifold dimension of the simplex.
"""
dimension(p::Simplex) = dimension(typeof(p))

"""
    length(simplex)

Returns the number of vertices (equals dimension + 1)
"""
Base.length(p::Simplex) = dimension(typeof(p))+1

"""
  universedimension(p)

Return the dimension of the universe in which `p` is embedded.
"""
universedimension(::Type{Simplex{U,D,C,N,T}}) where {U,D,C,N,T} = U
universedimension(p::Simplex) = universedimension(typeof(p))

"""
    getindex(simplex, I)
    simplex[I]

Get the vertices at index I (scalar or array) defining the simplex
"""
getindex(p::Simplex, I::Union{Number,SVector,Array}) = p.vertices[I]

"""
    simplex(vertices)
    simplex(v1, v2, ...)
    simplex(vertices, Val{D})

Build a D-dimensional simplex. The vertices can be passed in
an array (static or dynamic), or supplied separately. If the
length of the array is not part of its type, the speed of the
construction can be improved by supplying an extra Val{D}
argument. In case it is not clear from the context whether
the vertex array is dynamically or statically sized, use the
third form as it will not incur notable performance hits.

Note that D is the dimension of the simplex, i.e. the number
of vertices supplied minus one.
"""
@generated function simplex(vertices::SVector{D1,P}) where {D1,P}
    U = length(P)
    D = D1 - 1
    C = U-D
    T = eltype(P)
    xp1 =:(())
    for i in 1:D
        push!(xp1.args, :(vertices[$i]-vertices[end]))
    end
    xp2 = :(SVector{$D,P}($xp1))
    quote
        tangents = $xp2
        normals, volume = _normals(tangents, Val{$C})
        Simplex(vertices, tangents, normals, volume)
    end
end

simplex(vertices...) = simplex(SVector((vertices...,)))

@generated function simplex(vertices, ::Type{Val{D}}) where D
    P = eltype(vertices)
    D1 = D + 1
    xp = :(())
    for i in 1:D1
        push!(xp.args, :(vertices[$i]))
    end
    :(simplex(SVector{$D1,$P}($xp)))
end

normal(s::Simplex{3,2,1,3,T}) where {T} = normalize(cross(s[1]-s[3],s[2]-s[3]))

function _normals(tangents, ::Type{Val{1}})
    PT = eltype(tangents)
    D  = length(tangents)
    T  = eltype(PT)

    n = zeros(T,D+1)
    b = Array{T}(undef,D,D)

    for i in 1:D+1
        fill!(b, zero(T))
        for j in 1:D
            for k in 1:i-1
                b[k,j] = tangents[j][k]
            end
            for k in i:D
                b[k,j] = tangents[j][k+1]
            end
        end
        n[i] = (-1)^(i-1) * det(b)
    end

    n *= (-1)^D / norm(n)
    normals = SVector{1,PT}([PT(n)])

    metric = T[dot(tangents[i], tangents[j]) for i in 1:D, j in 1:D]
    volume = sqrt(abs(det(metric))) /  factorial(D)

    return normals, volume

end

function _normals(tangents::SVector{2,SVector{3,T}}, ::Type{Val{1}}) where {T}

    n = tangents[1] × tangents[2]
    l = norm(n)

    P = SVector{3,T}
    SVector{1,P}(n/l), 0.5*l
end

function _normals(tangents, ::Type{Val{C}}) where C
    PT = eltype(tangents)
    D  = length(tangents)
    U = length(PT)
    T  = eltype(PT)

    metric = T[dot(tangents[i], tangents[j]) for i in 1:D, j in 1:D]
    volume = sqrt(abs(det(metric))) / factorial(D)

    # Fix this. This function needs to become gneerated
    normals = SVector{C,PT}(PT[zero(PT) for i in 1:C])

    return normals, volume
end

"""
    barytocart(simplex, uv)

Returns the point in the simplex with barycentric coordinates uv
"""
function barytocart(mani::Simplex, u)
    r = last(mani.vertices)
    for i in 1 : dimension(mani)
        ti = mani.tangents[i]
        ui = u[i]
        #r += mani.tangents[i] * u[i]
        r += ti * ui
    end
    return r
end

"""
    carttobary(simplex, point) -> barycoords

Compute the barycentric coordinates on 'simplex' of 'point'.
"""
function carttobary(p::Simplex{U,D,C,N,T}, cart) where {U,D,C,N,T}

    G = [dot(p.tangents[i], p.tangents[j]) for i in 1:D, j in 1:D]
    #w = [dot(p.tangents[i], cart - p.vertices[end]) for i in 1:D]

    o = p.vertices[end]
    w = [sum(t[j]*(cart[j]-o[j]) for j in 1:length(cart)) for t in p.tangents]

    u = G \ w

    return SVector{D}(u)
end

function edges(s::Simplex{3,3})
    C = [SVector{2}(c) for c in combinations(@SVector[1,2,3,4],2)]
    T = eltype(eltype(s.vertices))
    P = Simplex{3,1,2,2,T}
    Edges = Vector{P}(undef, length(C))
    # @show eltype(Edges)
    # @show eltype(s.vertices)
    for c in C
        q = relorientation(c, @SVector[1,2,3,4])
        Q = abs(q)
        Edges[Q] = q > 0 ?
            simplex(s.vertices[c[1]], s.vertices[c[2]]) :
            simplex(s.vertices[c[2]], s.vertices[c[1]])
    end
    return Edges
end

function edges(s::Simplex{3,2})
    return [
        simplex(s.vertices[2], s.vertices[3]),
        simplex(s.vertices[3], s.vertices[1]),
        simplex(s.vertices[1], s.vertices[2])]
end

"""
    ReferenceSimplex{Dimension, CoordType, NumVertices}

This domain is defined to bootstrap the quadrature generation strategy. The generic definition
of numquads on a chart pulls back to the domain. For a limit set of reference domains, explicit
quadrature rules are defined. The weights and points are then pushed forward to the configuaration
space element over which integration is desired.

For more details see the implementation in quadpoints.jl
"""
struct ReferenceSimplex{D,T,N}
    # N: number of defining points
    # D: dimension of the simplex
    # T: type of the coordinates
    simplex::Simplex{D,D,0,N,T}
end

function ReferenceSimplex{D,T,N}() where {D,T,N}
    P = SVector{D,T}[]
    for i in 1:D
        a = zeros(T,D)
        a[i] = 1
        p = SVector{D,T}(a)
        push!(P,p)
    end

    o = zero(SVector{D,T})
    push!(P,o)
    ReferenceSimplex{D,T,N}(simplex(P...))
end

barytocart(ch::ReferenceSimplex, u) = barytocart(ch.simplex, u)
carttobary(ch::ReferenceSimplex, p) = carttobary(ch.simplex, p)

domain(ch::Simplex{U,D,C,N,T}) where {U,D,C,T,N} = ReferenceSimplex{D,T,N}()
neighborhood(ch::ReferenceSimplex, u) = u


text:
# my second attempt at writing a dCRAB alg
using OCToolbox
using Optim
using QuantumInformation

n_SI = 10
# generate frequency components
ω = rand(n_SI)

# we could write lots of these functions for any choice of basis
# there's already a library we could use I think
function gen_basis_fns(ω)
    f(A, B, t) = A * cos(2π * ω * t/T) + B * sin(2π * ω * t/T)
end

# sim stuff
K = 2
N = 10
T = 2
Δt = T/N

coeffs = zeros(N, K)

# functional that you want to minimise
function test(controls)
    controls = reshape(controls, (K, N))
    U = stack_props(controls, 0 * sz)
    C2(ρ, U * Ψ)
end

function full_propagator(c_vec, H0)
    exp((H0 + c_vec[1] .* sx .+ c_vec[2] .* sy) .* (-Δt*π*im))
end

function stack_props(c_mat,H0)
    init = I(2)
    for i in 1:size(c_mat)[2]
        init = full_propagator(c_mat[:,i], H0) * init
    end
    init
end

# we need some way to turn the continuous functions into something
# piecewise for now
function discretise(fn_list, coeffs, T, Δt)
    # fn_list contains a list of the functions
    discrete_time = range(0, (N-1) * Δt, step = Δt)
    envelope = zeros(length(discrete_time))
    for (i, fn) in enumerate(fn_list)
        envelope += fn.(coeffs[i, :]..., discrete_time)
    end
    envelope
end

fn_list = []
coeffs .* 0
for iteration in 1:n_SI
    @show iteration
    append!(fn_list, [gen_basis_fns(ω[iteration])])

    function fn_to_minimize(x)
        coeffs[iteration:iteration, :] = x
        control1 = discretise(fn_list, coeffs[1:iteration, :], T, Δt)
        # we use the same controls for both x and y because we are lazy
        controls = transpose(hcat(control1, control1))
        test(controls)
    end
    o = optimize(fn_to_minimize, [0.5, 0.5], NelderMead(), Optim.Options(show_trace = true))
    coeffs[iteration:iteration, :] = o.minimizer
end

coeffs

control1 = discretise(fn_list, coeffs, T, Δt)
controls = transpose(hcat(control1, control1))
test(controls)


text:
module Perfusion

include("utils.jl")
export @extract, interquartile_mean, percent_error, apply_mask, crop

using DICOM
include("load_dicom.jl")
export load_vfa_dicom, load_dce_dicom

using LinearAlgebra: norm
using LsqFit
using Statistics: mean
include("relaxation.jl")
export signal_to_concentration,
    concentration_to_signal,
    signal_to_R1,
    R1_to_signal,
    concentration_to_R1,
    R1_to_concentration
export fit_relaxation, fit_relaxation_nls, fit_relaxation_despot, fit_relaxation_novifast

using SpecialFunctions: gamma
include("aif.jl")
export aif_biexponential, aif_fritzhansen, aif_weinmann
export aif_orton1, aif_orton2, aif_orton3
export aif_georgiou, aif_parker

using LsqFit
using NamedTupleTools: select
using NumericalIntegration: cumul_integrate, TrapezoidalFast
using Statistics: mean, quantile, std
include("model_fitting.jl")
export model_tofts, model_exchange, model_filtration, model_uptake
export fit_model
export fit_tofts_nls, fit_tofts_lls
export fit_extendedtofts_nls, fit_extendedtofts_lls
export fit_uptake_nls, fit_uptake_lls
export fit_exchange_nls, fit_exchange_lls
export fit_filtration_nls, fit_filtration_lls
export fit_rrm_nls, fit_rrm_lls, fit_crrm_nls, fit_crrm_lls
export fit_errm_lls, fit_cerrm_lls
export fit_rrift, fit_rrift_with_crrm, fit_rrift_with_cerrm

using CancerImagingArchive
include("invivo_demo.jl")

end # module


text:
abstract type Macroscopic_quantity end

@with_kw mutable struct Macroquant{T, K} <: Macroscopic_quantity
    # In one spatial dimension T::Vector{Float32,64} <: Array{Float32,64, 1}
    # In two spatial dimension T::Matrix{Float32,64} <: Array{Float32.64, 2}
    # The macro with_kw adds a specific constructor operation to the struct.
    height::T
    velocity::K
    pressure::T
    energy::T
end

@with_kw mutable struct Forces{T} <: Macroscopic_quantity
    slip::T
    h∇p::T
    bathymetry::T
    thermal::T
end

@with_kw mutable struct Twovector{T} <: Macroscopic_quantity
    x::T
    y::T
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
# activate the test environment if needed
using Pkg
if ! ("PackageCompiler" ∈ keys(Pkg.project().dependencies))
    using TestEnv; TestEnv.activate()
end
@info "Loading packages ..."
using Dierckx, StaticArrays, LinearAlgebra, Parameters, NLsolve, DocStringExtensions, Sundials, KiteUtils, KitePodModels, Plots, AtmosphericModels
using PackageCompiler

@info "Creating sysimage ..."
push!(LOAD_PATH,joinpath(pwd(),"src"))

PackageCompiler.create_sysimage(
    [:Dierckx, :StaticArrays, :Parameters, :NLsolve, :DocStringExtensions, :Sundials, :KiteUtils, :KitePodModels, :Plots, :AtmosphericModels];
    sysimage_path="kps-image_tmp.so",
    precompile_execution_file=joinpath("test", "test_for_precompile.jl")
)