text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mq

using AWS.Compat
using AWS.UUIDs
"""
    CreateBroker()

Creates a broker. Note: This API is asynchronous.

# Optional Parameters
- `authenticationStrategy`: The authentication strategy used to secure the broker.
- `autoMinorVersionUpgrade`: Required. Enables automatic upgrades to new minor versions for brokers, as Apache releases the versions. The automatic upgrades occur during the maintenance window of the broker or after a manual broker reboot.
- `brokerName`: Required. The name of the broker. This value must be unique in your AWS account, 1-50 characters long, must contain only letters, numbers, dashes, and underscores, and must not contain whitespaces, brackets, wildcard characters, or special characters.
- `configuration`: A list of information about the configuration.
- `creatorRequestId`: The unique ID that the requester receives for the created broker. Amazon MQ passes your ID with the API action. Note: We recommend using a Universally Unique Identifier (UUID) for the creatorRequestId. You may omit the creatorRequestId if your application doesn't require idempotency.
- `deploymentMode`: Required. The deployment mode of the broker.
- `encryptionOptions`: Encryption options for the broker.
- `engineType`: Required. The type of broker engine. Note: Currently, Amazon MQ supports only ACTIVEMQ.
- `engineVersion`: Required. The version of the broker engine. For a list of supported engine versions, see https://docs.aws.amazon.com/amazon-mq/latest/developer-guide/broker-engine.html
- `hostInstanceType`: Required. The broker's instance type.
- `ldapServerMetadata`: The metadata of the LDAP server used to authenticate and authorize connections to the broker.
- `logs`: Enables Amazon CloudWatch logging for brokers.
- `maintenanceWindowStartTime`: The parameters that determine the WeeklyStartTime.
- `publiclyAccessible`: Required. Enables connections from applications outside of the VPC that hosts the broker's subnets.
- `securityGroups`: The list of security groups (1 minimum, 5 maximum) that authorizes connections to brokers.
- `storageType`: The broker's storage type.
- `subnetIds`: The list of groups (2 maximum) that define which subnets and IP ranges the broker can use from different Availability Zones. A SINGLE_INSTANCE deployment requires one subnet (for example, the default subnet). An ACTIVE_STANDBY_MULTI_AZ deployment requires two subnets.
- `tags`: Create tags when creating the broker.
- `users`: Required. The list of ActiveMQ users (persons or applications) who can access queues and topics. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.
"""

create_broker(; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/brokers", Dict{String, Any}("creatorRequestId"=>string(uuid4())); aws_config=aws_config)
create_broker(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/brokers", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("creatorRequestId"=>string(uuid4())), args)); aws_config=aws_config)

"""
    CreateConfiguration()

Creates a new configuration for the specified configuration name. Amazon MQ uses the default configuration (the engine type and version).

# Optional Parameters
- `authenticationStrategy`: The authentication strategy associated with the configuration.
- `engineType`: Required. The type of broker engine. Note: Currently, Amazon MQ supports only ACTIVEMQ.
- `engineVersion`: Required. The version of the broker engine. For a list of supported engine versions, see https://docs.aws.amazon.com/amazon-mq/latest/developer-guide/broker-engine.html
- `name`: Required. The name of the configuration. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 1-150 characters long.
- `tags`: Create tags when creating the configuration.
"""

create_configuration(; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/configurations"; aws_config=aws_config)
create_configuration(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/configurations", args; aws_config=aws_config)

"""
    CreateTags()

Add a tag to a resource.

# Required Parameters
- `resource-arn`: The Amazon Resource Name (ARN) of the resource tag.

# Optional Parameters
- `tags`: The key-value pair for the resource tag.
"""

create_tags(resource_arn; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/tags/$(resource-arn)"; aws_config=aws_config)
create_tags(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    CreateUser()

Creates an ActiveMQ user.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.

# Optional Parameters
- `consoleAccess`: Enables access to the the ActiveMQ Web Console for the ActiveMQ user.
- `groups`: The list of groups (20 maximum) to which the ActiveMQ user belongs. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.
- `password`: Required. The password of the user. This value must be at least 12 characters long, must contain at least 4 unique characters, and must not contain commas.
"""

create_user(broker_id, username; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/brokers/$(broker-id)/users/$(username)"; aws_config=aws_config)
create_user(broker_id, username, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/brokers/$(broker-id)/users/$(username)", args; aws_config=aws_config)

"""
    DeleteBroker()

Deletes a broker. Note: This API is asynchronous.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

"""

delete_broker(broker_id; aws_config::AWSConfig=global_aws_config()) = mq("DELETE", "/v1/brokers/$(broker-id)"; aws_config=aws_config)
delete_broker(broker_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("DELETE", "/v1/brokers/$(broker-id)", args; aws_config=aws_config)

"""
    DeleteTags()

Removes a tag from a resource.

# Required Parameters
- `resource-arn`: The Amazon Resource Name (ARN) of the resource tag.
- `tagKeys`: An array of tag keys to delete

"""

delete_tags(resource_arn, tagKeys; aws_config::AWSConfig=global_aws_config()) = mq("DELETE", "/v1/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
delete_tags(resource_arn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("DELETE", "/v1/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    DeleteUser()

Deletes an ActiveMQ user.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.

"""

delete_user(broker_id, username; aws_config::AWSConfig=global_aws_config()) = mq("DELETE", "/v1/brokers/$(broker-id)/users/$(username)"; aws_config=aws_config)
delete_user(broker_id, username, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("DELETE", "/v1/brokers/$(broker-id)/users/$(username)", args; aws_config=aws_config)

"""
    DescribeBroker()

Returns information about the specified broker.

# Required Parameters
- `broker-id`: The name of the broker. This value must be unique in your AWS account, 1-50 characters long, must contain only letters, numbers, dashes, and underscores, and must not contain whitespaces, brackets, wildcard characters, or special characters.

"""

describe_broker(broker_id; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers/$(broker-id)"; aws_config=aws_config)
describe_broker(broker_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers/$(broker-id)", args; aws_config=aws_config)

"""
    DescribeBrokerEngineTypes()

Describe available engine types and versions.

# Optional Parameters
- `engineType`: Filter response by engine type.
- `maxResults`: The maximum number of engine types that Amazon MQ can return per page (20 by default). This value must be an integer from 5 to 100.
- `nextToken`: The token that specifies the next page of results Amazon MQ should return. To request the first page, leave nextToken empty.
"""

describe_broker_engine_types(; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/broker-engine-types"; aws_config=aws_config)
describe_broker_engine_types(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/broker-engine-types", args; aws_config=aws_config)

"""
    DescribeBrokerInstanceOptions()

Describe available broker instance options.

# Optional Parameters
- `engineType`: Filter response by engine type.
- `hostInstanceType`: Filter response by host instance type.
- `maxResults`: The maximum number of instance options that Amazon MQ can return per page (20 by default). This value must be an integer from 5 to 100.
- `nextToken`: The token that specifies the next page of results Amazon MQ should return. To request the first page, leave nextToken empty.
- `storageType`: Filter response by storage type.
"""

describe_broker_instance_options(; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/broker-instance-options"; aws_config=aws_config)
describe_broker_instance_options(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/broker-instance-options", args; aws_config=aws_config)

"""
    DescribeConfiguration()

Returns information about the specified configuration.

# Required Parameters
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.

"""

describe_configuration(configuration_id; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations/$(configuration-id)"; aws_config=aws_config)
describe_configuration(configuration_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations/$(configuration-id)", args; aws_config=aws_config)

"""
    DescribeConfigurationRevision()

Returns the specified configuration revision for the specified configuration.

# Required Parameters
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.
- `configuration-revision`: The revision of the configuration.

"""

describe_configuration_revision(configuration_id, configuration_revision; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations/$(configuration-id)/revisions/$(configuration-revision)"; aws_config=aws_config)
describe_configuration_revision(configuration_id, configuration_revision, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations/$(configuration-id)/revisions/$(configuration-revision)", args; aws_config=aws_config)

"""
    DescribeUser()

Returns information about an ActiveMQ user.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.

"""

describe_user(broker_id, username; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers/$(broker-id)/users/$(username)"; aws_config=aws_config)
describe_user(broker_id, username, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers/$(broker-id)/users/$(username)", args; aws_config=aws_config)

"""
    ListBrokers()

Returns a list of all brokers.

# Optional Parameters
- `maxResults`: The maximum number of brokers that Amazon MQ can return per page (20 by default). This value must be an integer from 5 to 100.
- `nextToken`: The token that specifies the next page of results Amazon MQ should return. To request the first page, leave nextToken empty.
"""

list_brokers(; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers"; aws_config=aws_config)
list_brokers(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers", args; aws_config=aws_config)

"""
    ListConfigurationRevisions()

Returns a list of all revisions for the specified configuration.

# Required Parameters
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.

# Optional Parameters
- `maxResults`: The maximum number of configurations that Amazon MQ can return per page (20 by default). This value must be an integer from 5 to 100.
- `nextToken`: The token that specifies the next page of results Amazon MQ should return. To request the first page, leave nextToken empty.
"""

list_configuration_revisions(configuration_id; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations/$(configuration-id)/revisions"; aws_config=aws_config)
list_configuration_revisions(configuration_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations/$(configuration-id)/revisions", args; aws_config=aws_config)

"""
    ListConfigurations()

Returns a list of all configurations.

# Optional Parameters
- `maxResults`: The maximum number of configurations that Amazon MQ can return per page (20 by default). This value must be an integer from 5 to 100.
- `nextToken`: The token that specifies the next page of results Amazon MQ should return. To request the first page, leave nextToken empty.
"""

list_configurations(; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations"; aws_config=aws_config)
list_configurations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/configurations", args; aws_config=aws_config)

"""
    ListTags()

Lists tags for a resource.

# Required Parameters
- `resource-arn`: The Amazon Resource Name (ARN) of the resource tag.

"""

list_tags(resource_arn; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/tags/$(resource-arn)"; aws_config=aws_config)
list_tags(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    ListUsers()

Returns a list of all ActiveMQ users.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

# Optional Parameters
- `maxResults`: The maximum number of ActiveMQ users that can be returned per page (20 by default). This value must be an integer from 5 to 100.
- `nextToken`: The token that specifies the next page of results Amazon MQ should return. To request the first page, leave nextToken empty.
"""

list_users(broker_id; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers/$(broker-id)/users"; aws_config=aws_config)
list_users(broker_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("GET", "/v1/brokers/$(broker-id)/users", args; aws_config=aws_config)

"""
    RebootBroker()

Reboots a broker. Note: This API is asynchronous.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

"""

reboot_broker(broker_id; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/brokers/$(broker-id)/reboot"; aws_config=aws_config)
reboot_broker(broker_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("POST", "/v1/brokers/$(broker-id)/reboot", args; aws_config=aws_config)

"""
    UpdateBroker()

Adds a pending configuration change to a broker.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

# Optional Parameters
- `authenticationStrategy`: The authentication strategy used to secure the broker.
- `autoMinorVersionUpgrade`: Enables automatic upgrades to new minor versions for brokers, as Apache releases the versions. The automatic upgrades occur during the maintenance window of the broker or after a manual broker reboot.
- `configuration`: A list of information about the configuration.
- `engineVersion`: The version of the broker engine. For a list of supported engine versions, see https://docs.aws.amazon.com/amazon-mq/latest/developer-guide/broker-engine.html
- `hostInstanceType`: The host instance type of the broker to upgrade to. For a list of supported instance types, see https://docs.aws.amazon.com/amazon-mq/latest/developer-guide//broker.html#broker-instance-types
- `ldapServerMetadata`: The metadata of the LDAP server used to authenticate and authorize connections to the broker.
- `logs`: Enables Amazon CloudWatch logging for brokers.
- `securityGroups`: The list of security groups (1 minimum, 5 maximum) that authorizes connections to brokers.
"""

update_broker(broker_id; aws_config::AWSConfig=global_aws_config()) = mq("PUT", "/v1/brokers/$(broker-id)"; aws_config=aws_config)
update_broker(broker_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("PUT", "/v1/brokers/$(broker-id)", args; aws_config=aws_config)

"""
    UpdateConfiguration()

Updates the specified configuration.

# Required Parameters
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.

# Optional Parameters
- `data`: Required. The base64-encoded XML configuration.
- `description`: The description of the configuration.
"""

update_configuration(configuration_id; aws_config::AWSConfig=global_aws_config()) = mq("PUT", "/v1/configurations/$(configuration-id)"; aws_config=aws_config)
update_configuration(configuration_id, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("PUT", "/v1/configurations/$(configuration-id)", args; aws_config=aws_config)

"""
    UpdateUser()

Updates the information for an ActiveMQ user.

# Required Parameters
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: Required. The username of the ActiveMQ user. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.

# Optional Parameters
- `consoleAccess`: Enables access to the the ActiveMQ Web Console for the ActiveMQ user.
- `groups`: The list of groups (20 maximum) to which the ActiveMQ user belongs. This value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100 characters long.
- `password`: The password of the user. This value must be at least 12 characters long, must contain at least 4 unique characters, and must not contain commas.
"""

update_user(broker_id, username; aws_config::AWSConfig=global_aws_config()) = mq("PUT", "/v1/brokers/$(broker-id)/users/$(username)"; aws_config=aws_config)
update_user(broker_id, username, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mq("PUT", "/v1/brokers/$(broker-id)/users/$(username)", args; aws_config=aws_config)


text:
abstract type Weight end

Base.show(io::IO, w::Weight) = print(io, name(w))

function Base.:(==)(o1::Weight, o2::Weight)
    typeof(o1) == typeof(o2) || return false
    all(x -> getfield(o1, x) == getfield(o2, x), fieldnames(typeof(o1)))
end
Base.copy(w::Weight) = deepcopy(w)

#-----------------------------------------------------------------------# EqualWeight
"""
    EqualWeight()

Equally weighted observations.

``γ(t) = 1 / t``
"""
struct EqualWeight <: Weight end
(::EqualWeight)(n) = inv(n)

#-----------------------------------------------------------------------# ExponentialWeight
"""
    ExponentialWeight(λ::Float64)
    ExponentialWeight(lookback::Int)

Exponentially weighted observations.  Each weight is `λ = 2 / (lookback + 1)`.

`ExponentialWeight` does not satisfy the usual assumption that `γ(1) == 1`.  Therefore, some
statistics have an implicit starting value.

```
# E.g. Mean has an implicit starting value of 0.
o = Mean(weight=ExponentialWeight(.1))
fit!(o, 10)
value(o) == 1
```

``γ(t) = λ``
"""
struct ExponentialWeight <: Weight
    λ::Float64
    ExponentialWeight(λ::Real = .1) = new(λ)
    ExponentialWeight(lookback::Integer) = new(2 / (lookback + 1))
end
(w::ExponentialWeight)(n) = w.λ
Base.show(io::IO, w::ExponentialWeight) = print(io, name(w) * "(λ = $(w.λ))")

#-----------------------------------------------------------------------# LearningRate
"""
    LearningRate(r = .6)

Slowly decreasing weight.  Satisfies the standard stochastic approximation assumption
``∑ γ(t) = ∞, ∑ γ(t)^2 < ∞`` if ``r ∈ (.5, 1]``.

``γ(t) = inv(t ^ r)``
"""
struct LearningRate <: Weight
    r::Float64
    LearningRate(r = .6) = new(r)
end
(w::LearningRate)(n) = inv(n ^ w.r)
Base.show(io::IO, w::LearningRate) = print(io, name(w) * "(r = $(w.r))")

#-----------------------------------------------------------------------# LearningRate2
"""
    LearningRate2(c = .5)

Slowly decreasing weight.

``γ(t) = inv(1 + c * (t - 1))``
"""
struct LearningRate2 <: Weight
    c::Float64
    LearningRate2(c = .5) = new(c)
end
(w::LearningRate2)(n) = 1 / (1 + w.c * (n - 1))
Base.show(io::IO, w::LearningRate2) = print(io, name(w) * "(c = $(w.c))")

#-----------------------------------------------------------------------# HarmonicWeight
"""
    HarmonicWeight(a = 10.0)

Weight determined by harmonic series.

``γ(t) = a / (a + t - 1)``
"""
struct HarmonicWeight <: Weight
    a::Float64
    HarmonicWeight(a = 10.0) = new(a)
end
(w::HarmonicWeight)(n) = w.a / (w.a + n - 1)
Base.show(io::IO, w::HarmonicWeight) = print(io, name(w) * "(a = $(w.a))")

#-----------------------------------------------------------------------# McclainWeight
"""
    McclainWeight(α = .1)

Weight which decreases into a constant.

``γ(t) = γ(t-1) / (1 + γ(t-1) - α)``
"""
mutable struct McclainWeight <: Weight
    α::Float64
    last::Float64
    McclainWeight(α = .1) = new(α, 1.0)
end
(w::McclainWeight)(n) = n == 1 ? 1.0 : (w.last = w.last / (1 + w.last - w.α))
Base.show(io::IO, w::McclainWeight) = print(io, name(w) * "(α = $(w.α))")

ext: jl


text:
# License for this file: MIT (expat)
# Copyright 2020, DLR Institute of System Dynamics and Control
# Developer: Martin Otter, DLR-SR
#
# This file is part of module ModiaResult
#
# Utility functions that are usually not directly called.

using  Unitful
import OrderedCollections
import DataFrames
import Measurements
import MonteCarloMeasurements
import Pkg
isinstalled(pkg::String) = any(x -> x.name == pkg && x.is_direct_dep, values(Pkg.dependencies()))

const AvailableModiaPlotPackages = ["GLMakie", "WGLMakie", "CairoMakie", "PyPlot", "NoPlot", "SilentNoPlot"]
const ModiaPlotPackagesStack = String[]

"""
    usePlotPackage(plotPackage::String)

Define the ModiaPlot package that shall be used by command `ModiaResult.@usingModiaPlot`.
If a ModiaPlot package is already defined, save it on an internal stack
(can be reactivated with `usePreviousPlotPackage()`.

Possible values for `plotPackage`:
- `"GLMakie"`
- `"WGLMakie"`
- `"CairoMakie"`
- `"PyPlot"`
- `"NoPlot"`
- `"SilentNoPlot"`

# Example

```julia
import ModiaResult

ModiaResult.usePlotPackage("GLMakie")

module MyTest
    ModiaResult.@usingModiaPlot

    t = range(0.0, stop=10.0, length=100)
    result = Dict{String,Any}("time" => t, "phi" => sin.(t))

    plot(result, "phi")  # use GLMakie for the rendering
end
```
"""
function usePlotPackage(plotPackage::String; pushPreviousOnStack=true)::Bool
    success = true
    if plotPackage == "NoPlot" || plotPackage == "SilentNoPlot"
        if  pushPreviousOnStack && haskey(ENV, "MODIA_PLOT")
            push!(ModiaPlotPackagesStack, ENV["MODIA_PLOT"])
        end
        if plotPackage == "NoPlot"
            ENV["MODIA_PLOT"] = "NoPlot"
        else
            ENV["MODIA_PLOT"] = "SilentNoPlot"
        end
    else
        plotPackageName = "ModiaPlot_" * plotPackage
        if plotPackage in AvailableModiaPlotPackages
            # Check that plotPackage is defined in current environment
            if isinstalled(plotPackageName)
                if pushPreviousOnStack && haskey(ENV, "MODIA_PLOT")
                    push!(ModiaPlotPackagesStack, ENV["MODIA_PLOT"])
                end
                ENV["MODIA_PLOT"] = plotPackage
            else
                @warn "... usePlotPackage(\"$plotPackage\"): Call ignored, since package $plotPackageName is not in your current environment"
                success = false
            end
        else
            @warn "\n... usePlotPackage(\"$plotPackage\"): Call ignored, since argument not in $AvailableModiaPlotPackages."
            success = false
        end
    end
    return success
end

"""
    usePreviousPlotPackage()

Pop the last saved ModiaPlot package from an internal stack
and call `usePlotPackage(<popped ModiaPlot package>)`.
"""
function usePreviousPlotPackage()::Bool
    if length(ModiaPlotPackagesStack) > 0
        plotPackage = pop!(ModiaPlotPackagesStack)
        success = usePlotPackage(plotPackage, pushPreviousOnStack=false)
    else
        @warn "usePreviousPlotPackage(): Call ignored, because nothing saved."
        success = false
    end
    return success
end

"""
    currentPlotPackage()

Return the name of the plot package as a string that was
defined with [`usePlotPackage`](@ref).
For example, the function may return "GLMakie", "PyPlot" or "NoPlot" or
or "", if no PlotPackage is defined.
"""
currentPlotPackage() = haskey(ENV, "MODIA_PLOT") ? ENV["MODIA_PLOT"] : ""

"""
    @usingModiaPlot()

Execute `using XXX`, where `XXX` is the ModiaPlot package that was
activated with `usePlotPackage(plotPackage)`.
"""
macro usingModiaPlot()
    if haskey(ENV, "MODIA_PLOT")
        ModiaPlotPackage = ENV["MODIA_PLOT"]
        if !(ModiaPlotPackage in AvailableModiaPlotPackages)
            @warn "ENV[\"MODIA_PLOT\"] = \"$ModiaPlotPackage\" is not supported!. Using \"NoPlot\"."
            @goto USE_NO_PLOT
        elseif ModiaPlotPackage == "NoPlot"
            @goto USE_NO_PLOT
        elseif ModiaPlotPackage == "SilentNoPlot"
            expr = :( import ModiaResult.SilentNoPlot: plot, showFigure, saveFigure, closeFigure, closeAllFigures )
            return esc( expr )
        else
            ModiaPlotPackage = Symbol("ModiaPlot_" * ModiaPlotPackage)
            expr = :(using $ModiaPlotPackage)
            println("$expr")
            return esc( :(using $ModiaPlotPackage) )
        end

    else
        @warn "No plot package activated. Using \"NoPlot\"."
        @goto USE_NO_PLOT
    end

    @label USE_NO_PLOT
    expr = :( import ModiaResult.NoPlot: plot, showFigure, saveFigure, closeFigure, closeAllFigures )
    println("$expr")
    return esc( expr )
end

const signalTypeToString = ["Independent", "Continuous", "Clocked"]

"""
    resultInfo(result)

Return information about the result as DataFrames.DataFrame object
with columns:

```julia
name::String, unit::String, nTime::String, signalType::String, valueSize::String, eltype::String
```
"""
function resultInfo(result)
    if isnothing(result)
        @info "The call of showInfo(result) is ignored, since the argument is nothing."
        return
    end

    resultInfoTable = DataFrames.DataFrame(name=String[], unit=String[], nTime=String[], signalType=String[], valueSize=String[], eltype=String[])

    timeSigName = timeSignalName(result)
    for name in signalNames(result)
        (signalType, nTime, valueSize, elType, sigUnit) = signalInfo(result, name)
        if isnothing(elType)
            sigUnit2    = "???"
            nTime2      = "???"
            signalType2 = "???"
            valueSize2  = "???"
            elType2     = "???"

        else
            sigUnit2    = string(sigUnit)
            nTime2      = name==timeSigName && !hasOneTimeSignal(result) ? "---" : string(nTime)
            signalType2 = signalTypeToString[Int(signalType)]
            valueSize2  = string(valueSize)
            elType2     = string(elType)
        end

        push!(resultInfoTable, [name, sigUnit2, nTime2, signalType2, valueSize2, elType2] )
    end

    return resultInfoTable
end

"""
    printResultInfo(result)

Print info about result.

# Example
```julia
using ModiaResult
using Unitful
ModiaResult.@usingModiaPlot

t = range(0.0, stop=10.0, length=100)
result = OrderedDict{String,Any}("time"=> t*u"s", "phi" => sin.(t)*u"rad")
printResultInfo(result)

# Gives output:
 # │ name  unit  nTime  signalType  valueSize  eltype
───┼───────────────────────────────────────────────────
 1 │ time  s     100    Independent ()         Float64
 2 │ phi   rad   100    Continuous  ()         Float64
```

"""
function printResultInfo(result)::Nothing
    resultInfoTable = resultInfo(result)
    show(stdout, resultInfoTable, summary=false, rowlabel=Symbol("#"), allcols=true, eltypes=false, truncate=50)
    println(stdout)
    return nothing
end

"""
    ResultDict(args...; defaultHeading="", hasOneTimeSignal=true)

Return a new ResultDict dictionary (is based on OrderedCollections.OrderedDict).

- A key of the dictionary is a String. Key `"time"` characterizes the
  independent variable.

- A value of the dictionary is a tuple `(timeSignal, signal, signalType)`
  where `timeSignal::Vector{AbstractVector}`,
  `signal::Vector{AbstractVector}` and
  `signalType::ModiaResult.SignalType.
  `signal[i][j]` is the signalValue at time instant `timeSignal[i][j]` in segment `i`.

# Example

```julia
using ModiaResult

time0 = [0.0, 7.0]
t     = ([time0], [time0], ModiaResult.Independent)

time1 = 0.0 : 0.1  : 2.0
time2 = 3.0 : 0.01 : 3.5
time3 = 5.0 : 0.1  : 7.0
sigA1 = sin.(time1)u"m"
sigA2 = cos.(time2)u"m"
sigA3 = sin.(time3)u"m"
sigA  = ([time1, time2, time3],
         [sigA1, sigA2, sigA3 ],
         ModiaResult.SignalType)
sigB  = ([time2], [sin.(time2)], ModiaResult.SignalType)
sigC  = ([time3], [sin.(time3)], ModiaResult.Clocked)

result = ModiaResult.ResultDict("time" => t,
                                "sigA" => sigA,
                                "sigB" => sigB,
                                "sigC" => sigC,
                                defaultHeading = "Three test signals",
                                hasOneTimeSignal = false)
showInfo(result)
```
"""
struct ResultDict    <: AbstractDict{String,Tuple{Any,Any,SignalType}}
    dict::OrderedCollections.OrderedDict{String,Tuple{Any,Any,SignalType}}
    defaultHeading::String
    hasOneTimeSignal::Bool

    ResultDict(args...; defaultHeading="", hasOneTimeSignal=true) =
        new(OrderedCollections.OrderedDict{String,Tuple{Any,Any,SignalType}}(args...),
            defaultHeading, hasOneTimeSignal)
end

        #new(OrderedCollections.OrderedDict{String,Tuple{Vector{AbstractVector},
        #                                            Vector{AbstractVector},
        #                                            ModiaResult.SignalType}}(args...),

# Overload AbstractDict methods
Base.haskey(result::ResultDict, key) = Base.haskey(result.dict, key)

Base.get(result::ResultDict, key, default)     = Base.get(result.dict, key, default)
Base.get(f::Function, result::ResultDict, key) = Base.get(f, result.dict, key)

Base.get!(result::ResultDict, key, default)     = Base.get!(result.dict, key, default)
Base.get!(f::Function, result::ResultDict, key) = Base.get!(f, result.dict, key)

Base.getkey(result::ResultDict, key, default) = Base.getkey(result.dict, key, default)

Base.delete!(result::ResultDict, key) = Base.delete!(result.dict, key)

Base.keys(result::ResultDict) = Base.keys(result.dict)

Base.pop!(result::ResultDict, key)          = Base.pop!(result.dict, key)
Base.pop!(result::ResultDict, key, default) = Base.pop!(result.dict, key, default)

Base.setindex!(result::ResultDict, value, key...) = Base.setindex!(result.dict, value, key...)

Base.values(result::ResultDict) = Base.values(result.dict)

"""
    signalLength(signal)

Return the total number of values of `signal::Vector{AbstractVector}`.
If signal[i] is nothing or missing, a length of zero is returned.
"""
function signalLength(signal::AbstractVector)
    for s in signal
        if ismissing(s) || isnothing(s)
            return 0
        end
    end
    return sum( length(s) for s in signal )
end

"""
    hasSameSegments(signal1, signal2)

Return true, if the lengths of the segments in `signal1` and in `signal2` are the same.
"""
function hasSameSegments(signal1::Vector{AbstractVector}, signal2::Vector{AbstractVector})
    if length(signal1) != length(signal2)
        return false
    end

    for i = 1:length(signal1)
        if length(signal1[i]) != length(signal2[i])
            return false
        end
    end

    return true
end

"""
    hasDimensionMismatch(signal, timeSignal, timeSignalName)

Print a warning message if signalLength(signal) != signalLength(timeSignal)
and return true. Otherwise, return false
"""
function hasDimensionMismatch(signal, signalName, timeSignal, timeSignalName::AbstractString)
    if signalLength(signal) != signalLength(timeSignal)
        lensignal = signalLength(signal)
        lentime   = signalLength(timeSignal)
        @warn "signalLength of \"$signalName\" = $lensignal but signalLength of \"$timeSignalName\" = $lentime"
        return true
    end
    return false
end

"""
    (sigType, nTime, sigSize, sigElType, sigUnit) = signalInfo(result, name)

Return information about a signal, given the `name` of the signal in `result`:

- `sigType::SignalType`: Ìnterpolation type of signal.

- `nTime::Int`: Number of signal time points.

- `sigSize`: size(signal[1][1])

- `sigElType`: ustrip( eltype(signal[1][1]) ), that is the element type of the signal without unit.

- `sigUnit`: Unit of signal

If `name` is defined, but no signal is available (= nothing, missing or zero length),
return `nTime=0` and `nothing` for `sigSize, sigElType, sigUnit`.
"""
function signalInfo(result, name::AbstractString)
    (timeSignal, signal, sigType) = rawSignal(result,name)
    if ismissing(signal) || isnothing(signal) || signalLength(signal) == 0 ||
       hasDimensionMismatch(signal, name, timeSignal, timeSignalName(result))
        return (sigType, 0, nothing, nothing, nothing)
    end

    value     = signal[1][1]
    valueSize = size(value)
    valueUnit = unit(value[1])

    if typeof(value) <: MonteCarloMeasurements.Particles
        elTypeAsString = string(typeof(ustrip.(value[1])))
        nparticles     = length(value)
        valueElType    = "MonteCarloMeasurements.Particles{" * elTypeAsString * ",$nparticles}"
    elseif typeof(value) <: MonteCarloMeasurements.StaticParticles
        elTypeAsString = string(typeof(ustrip.(value[1])))
        nparticles     = length(value)
        valueElType    = "MonteCarloMeasurements.StaticParticles{" * elTypeAsString * ",$nparticles}"
    else
        valueElType = typeof( ustrip.(value) )
    end
    nTime = signalLength(timeSignal)
    return (sigType, nTime, valueSize, valueElType, valueUnit)
end

# Default implementation of getSignalDetails
function getSignalDetails(result, name::AbstractString)
    sigPresent = false
    if hasSignal(result, name)
        (timeSig, sig2, sigType) = rawSignal(result, name)
        timeSigName = timeSignalName(result)
        if !( isnothing(sig2) || ismissing(sig2) || signalLength(sig2) == 0 ||
              hasDimensionMismatch(sig2, name, timeSig, timeSigName) )
            sigPresent = true
            value      = sig2[1][1]
            if ndims(value) == 0
                sig            = sig2
                arrayName      = name
                arrayIndices   = ()
                nScalarSignals = 1
            else
                arrayName      = name
                arrayIndices   = Tuple(1:Int(ni) for ni in size(value))
                nScalarSignals = length(value)
                sig = Vector{Matrix{eltype(value)}}(undef, length(sig2))
                for segment = 1:length(sig2)
                    sig[segment] = zeros(eltype(value), length(sig2[segment]), nScalarSignals)
                    siga  = sig[segment]
                    sig2a = sig2[segment]
                    for (i, value_i) in enumerate(sig2a)
                        for j in 1:nScalarSignals
                            siga[i,j] = sig2a[i][j]
                        end
                    end
                end
            end
        end

    else
        # Handle signal arrays, such as a.b.c[3] or a.b.c[2:3, 1:5, 3]
        if name[end] == ']'
            i = findlast('[', name)
            if i >= 2
                arrayName = name[1:i-1]
                indices   = name[i+1:end-1]
                if hasSignal(result, arrayName)
                    (timeSig, sig2, sigType) = rawSignal(result, arrayName)
                    timeSigName = timeSignalName(result)
                    if !( isnothing(sig2) || ismissing(sig2) || signalLength(sig2) == 0 ||
                          hasDimensionMismatch(sig2, arrayName, timeSig, timeSigName) )
                        sigPresent = true
                        value = sig2[1][1]

                        # Determine indices as tuple
                        arrayIndices = eval( Meta.parse( "(" * indices * ",)" ) )

                        # Determine number of signals
                        #nScalarSignals = sum( length(indexRange) for indexRange in arrayIndices )

                        # Extract sub-matrix
                        sig = Vector{Any}(undef,length(sig2))
                        for segment = 1:length(sig2)
                            sig2a = sig2[segment]
                            sig[segment] = [getindex(sig2a[i], arrayIndices...) for i in eachindex(sig2a)]
                        end

                        # Determine number of signals
                        nScalarSignals = length(sig[1][1])

                        # "flatten" array to matrix
                        eltypeValue = eltype(value)
                        if !(eltypeValue <: Number)
                            @warn "eltype($name) = $eltypeValue and this is not <: Number!"
                            return (nothing, nothing, nothing, nothing, name, (), 0)
                        end
                        for segment = 1:length(sig2)
                            sig[segment] = zeros(eltypeValue, length(sig2[segment]), nScalarSignals)
                            siga  = sig[segment]
                            sig2a = sig2[segment]
                            for (i, value_i) in enumerate(sig2a)
                                for j in 1:nScalarSignals
                                    siga[i,j] = sig2a[i][j]
                                end
                            end
                        end
                    end
                end
            end
        end
    end

    if sigPresent
        return (sig, timeSig, timeSigName, sigType, arrayName, arrayIndices, nScalarSignals)
    else
        return (nothing, nothing, nothing, nothing, name, (), 0)
    end
end

"""
    (signal, timeSignal, timeSignalName, signalType, arrayName, arrayIndices, nScalarSignals) =
         getSignalWithWarning(result, name)

Call getSignal(result,name) and print a warning message if `signal == nothing`
"""
function getSignalDetailsWithWarning(result,name::AbstractString)
    (sig, timeSig, timeSigName, sigType, arrayName, arrayIndices, nScalarSignals) = getSignalDetails(result,name)
    if isnothing(sig)
        @warn "\"$name\" is not correct or is not defined or has no values."
    end
    return (sig, timeSig, timeSigName, sigType, arrayName, arrayIndices, nScalarSignals)
end

appendUnit2(name, unit) = unit == "" ? name : string(name, " [", unit, "]")

function appendUnit(name, value)
    if typeof(value) <: MonteCarloMeasurements.StaticParticles ||
       typeof(value) <: MonteCarloMeasurements.Particles
        appendUnit2(name, string(unit(value.particles[1])))
    else
        appendUnit2(name, string(unit(value)))
    end
end

"""
    (xsig, xsigLegend, ysig, ysigLegend, ysigType) = getPlotSignal(result, ysigName; xsigName=nothing)

Given the result data structure `result` and a variable `ysigName::AbstractString` with
or without array range indices (for example `ysigName = "a.b.c[2,3:5]"`) and an optional
variable name `xsigName::AbstractString` for the x-axis, return

- `xsig::Vector{T1<:Real}`: The vector of the x-axis signal without a unit. Segments are concatenated
  and separated by NaN.

- `xsigLegend::AbstractString`: The legend of the x-axis consisting of the x-axis name
  and its unit (if available).

- `ysig::Vector{T2}` or `::Matrix{T2}`: the y-axis signal, either as a vector or as a matrix
  of values without units depending on the given name. For example, if `ysigName = "a.b.c[2,3:5]"`, then
  `ysig` consists of a matrix with three columns corresponding to the variable values of
  `"a.b.c[2,3]", "a.b.c[2,4]", "a.b.c[2,5]"` with the (potential) units are stripped off.
  Segments are concatenated and separated by NaN.

- `ysigLegend::Vector{AbstractString}`: The legend of the y-axis as a vector
  of strings, where `ysigLegend[1]` is the legend for `ysig`, if `ysig` is a vector,
  and `ysigLegend[i]` is the legend for the i-th column of `ysig`, if `ysig` is a matrix.
  For example, if variable `"a.b.c"` has unit `m/s`, then `ysigName = "a.b.c[2,3:5]"` results in
  `ysigLegend = ["a.b.c[2,3] [m/s]", "a.b.c[2,3] [m/s]", "a.b.c[2,5] [m/s]"]`.

- `ysigType::`[`SignalType`](@ref): The signal type of `ysig` (either `ModiaResult.Continuous`
  or `ModiaResult.Clocked`).

If `ysigName` is not valid, or no signal values are available, the function returns
`(nothing, nothing, nothing, nothing, nothing)`, and prints a warning message.
"""
function getPlotSignal(result, ysigName::AbstractString; xsigName=nothing)
    (ysig, xsig, timeSigName, ysigType, ysigArrayName, ysigArrayIndices, nysigScalarSignals) = getSignalDetailsWithWarning(result, ysigName)

    # Check y-axis signal and time signal
    if isnothing(ysig) || isnothing(xsig) || isnothing(timeSigName)  || signalLength(ysig) == 0
        @goto ERROR
    end

    # Get xSigName or check xSigName
    if isnothing(xsigName)
        xsigName = timeSigName
    elseif xsigName != timeSigName
        (xsig, xsigTime, xsigTimeName, xsigType, xsigArrayName, xsigArrayIndices, nxsigScalarSignals) = getSignalDetailsWithWarning(result, xsigName)
        if isnothing(xsig) || isnothing(xsigTime) || isnothing(xsigTimeName) || signalLength(xsig) == 0
            @goto ERROR
        elseif !hasSameSegments(ysig, xsig)
            @warn "\"$xsigName\" (= x-axis) and \"$ysigName\" (= y-axis) have not the same time signal vector."
            @goto ERROR
        end
    end

    # Check x-axis signal
    xsigValue = first(first(xsig))
    if length(xsigValue) != 1
        @warn "\"$xsigName\" does not characterize a scalar variable as needed for the x-axis."
        @goto ERROR
    elseif !( typeof(xsigValue) <: Number )
        @warn "\"$xsigName\" has no Number type values, but values of type " * string(typeof(xsigValue)) * "."
        @goto ERROR
    elseif typeof(xsigValue) <: Measurements.Measurement
        @warn "\"$xsigName\" is a Measurements.Measurement type and this is not (yet) supported for the x-axis."
        @goto ERROR
    elseif typeof(xsigValue) <: MonteCarloMeasurements.StaticParticles
        @warn "\"$xsigName\" is a MonteCarloMeasurements.StaticParticles type and this is not supported for the x-axis."
        @goto ERROR
    elseif typeof(xsigValue) <: MonteCarloMeasurements.Particles
        @warn "\"$xsigName\" is a MonteCarloMeasurements.Particles type and this is not supported for the x-axis."
        @goto ERROR
    end

    # Build xsigLegend
    xsigLegend = appendUnit(xsigName, xsigValue)

    # Get one segment of the y-axis and check it
    ysegment1 = first(ysig)
    if !( typeof(ysegment1) <: AbstractVector || typeof(ysegment1) <: AbstractMatrix )
        @error "Bug in function: typeof of an y-axis segment is neither a vector nor a Matrix, but " * string(typeof(ysegment1))
    elseif !(eltype(ysegment1) <: Number)
        @warn "\"$ysigName\" has no Number values but values of type " * string(eltype(ysegment1))
        @goto ERROR
    end

    # Build ysigLegend
    value = ysegment1[1]
    if ysigArrayIndices == ()
        # ysigName is a scalar variable
        ysigLegend = [appendUnit(ysigName, value)]

    else
        # ysigName is an array variable
        ysigLegend = [ysigArrayName * "[" for i = 1:nysigScalarSignals]
        i = 1
        ySizeLength = Int[]
        for j1 in eachindex(ysigArrayIndices)
            push!(ySizeLength, length(ysigArrayIndices[j1]))
            i = 1
            if j1 == 1
                for j2 in 1:div(nysigScalarSignals, ySizeLength[1])
                    for j3 in ysigArrayIndices[1]
                        ysigLegend[i] *= string(j3)
                        i += 1
                    end
                end
            else
                ncum = prod( ySizeLength[1:j1-1] )
                for j2 in ysigArrayIndices[j1]
                    for j3 = 1:ncum
                        ysigLegend[i] *= "," * string(j2)
                        i += 1
                    end
                end
            end
        end

        for i = 1:nysigScalarSignals
            ysigLegend[i] *= appendUnit("]", ysegment1[1,i])
        end
    end

    #xsig2 = Vector{Any}(undef, length(xsig))
    #ysig2 = Vector{Any}(undef, length(ysig))
    #for i = 1:length(xsig)
    #    xsig2[i] = collect(ustrip.(xsig[i]))
    #    ysig2[i] = collect(ustrip.(ysig[i]))
    #end

    xsig2 = collect(ustrip.(first(xsig)))
    ysig2 = collect(ustrip.(first(ysig)))
    if length(xsig) > 1
        xNaN = convert(eltype(xsig2), NaN)
        if ndims(ysig2) == 1
            yNaN = convert(eltype(ysig2), NaN)
        else
            yNaN = fill(convert(eltype(ysig2), NaN), 1, size(ysig2,2))
        end

        for i = 2:length(xsig)
            xsig2 = vcat(xsig2, xNaN, collect(ustrip.(xsig[i])))
            ysig2 = vcat(ysig2, yNaN, collect(ustrip.(ysig[i])))
        end
    end

    return (xsig2, xsigLegend, ysig2, ysigLegend, ysigType)

    @label ERROR
    return (nothing, nothing, nothing, nothing, nothing)
end

"""
    getHeading(result, heading)

Return `heading` if no empty string. Otherwise, return `defaultHeading(result)`.
"""
getHeading(result, heading::AbstractString) = heading != "" ? heading : defaultHeading(result)

"""
    prepend!(prefix, signalLegend)

Add `prefix` string in front of every element of the `signalLegend` string-Vector.
"""
function prepend!(prefix::AbstractString, signalLegend::Vector{AbstractString})
   for i in eachindex(signalLegend)
      signalLegend[i] = prefix*signalLegend[i]
   end
   return signalLegend
end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
module TimeSeriesA2P1

#Use this to print
#=
using Weave
codeName = "TimeSeriesA2P1"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, CTIO, CTReg, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]

const DEFAULT_BURN = 1_000
const NUM_SIMS = 10_000
const DEFAULT_LAG = 5

BLAS.set_num_threads(1)

mutable struct MCSpec
  T::Int
  ρ::Float64
  burn::Int
end

#generates the VAR data
function getVARData!(MC::MCSpec,
    y::Vector{Float64} = Vector{Float64}(MC.T),
    z::Vector{Float64} = Vector{Float64}(MC.T))::NTuple{2, Vector{Float64}}

  #pregenerate RVs
  εy::Vector{Float64} = rand(Normal(), MC.T+MC.burn)
  εz::Vector{Float64} = rand(Normal(), MC.T+MC.burn)

  #first do the burn
  yburn::Float64 = 0.0
  zburn::Float64 = 0.0

  @fastmath for i::Int ∈ 1:(MC.burn+1)
    yburn = MC.ρ*yburn + εy[i]
    zburn = MC.ρ*zburn + εz[i]
  end

  y[1] = yburn
  z[1] = zburn

  @fastmath for i::Int ∈ 2:MC.T
    y[i] = MC.ρ*y[i-1] + εy[i+MC.burn]
    z[i] = MC.ρ*z[i-1] + εz[i+MC.burn]
  end

  return (y,z)
end

mutable struct VARReg
  N::Int
  MC::MCSpec
  lag::Int
  βs::Matrix{Float64}
  SEs::Matrix{Float64}
  Ts::Matrix{Float64}
end

function VARReg(N::Int, MC::MCSpec; lag::Int=DEFAULT_LAG, verbose::Bool = false,
    βs::Matrix{Float64} = Matrix{Float64}(N,2),
    SEs::Matrix{Float64} = Matrix{Float64}(N,2),
    Ts::Matrix{Float64} = Matrix{Float64}(N,2))::VARReg

  #preallocate regression components
  X::Matrix{Float64} = ones(MC.T,2)
  Y::Vector{Float64} = Vector{Float64}(MC.T)
  XNames::Vector{Symbol} = [:intercept, :X]
  YName::Symbol = :Y

  #container for holding the Z data
  Z::Vector{Float64} = Vector{Float64}(MC.T)

  for i::Int ∈ 1:N
    getVARData!(MC, Y, Z)
    X[:, 2] .= Z
    reg::CTLM = CTLM(X, Y, XNames=XNames, YName=YName)
    βs[i,:] .= reg.β
    SEs[i,:] .= (diag(getNeweyWest!(reg, lag))).^0.5
    Ts[i,:] .= βs[i,:] ./ SEs[i,:]

    if verbose
      println("NW: $(SEs[i,:]), NWSlow: $((diag(getNeweyWestSlow(reg, lag))).^0.5)")
    end

  end

  return VARReg(N, MC, lag, βs, SEs, Ts)
end

#convenience method for the above
VARReg(;N::Int=NUM_SIMS, T::Int=50, ρ::Float64=1.0, burn::Int = DEFAULT_BURN, lag=DEFAULT_LAG,
    βs::Matrix{Float64} = Matrix{Float64}(N,2),
    SEs::Matrix{Float64} = Matrix{Float64}(N,2),
    Ts::Matrix{Float64} = Matrix{Float64}(N,2))::VARReg =
  VARReg(N, MCSpec(T, ρ, burn), lag=lag, βs=βs, SEs=SEs, Ts=Ts)

function visualizeA2P1(VRs::Vector{VARReg})

  #convenience containers
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  #convenience method for number formatting
  n2s(x::T where T<:Real) = num2Str(x, 3, Ints=true)

  push!(plotNames, "Betas")
  push!(plots, Vector{Plot}())
  for VR::VARReg ∈ VRs
    push!(plots[end],
      plot(x=VR.βs[:,2], Geom.histogram(bincount=20, density=true, position=:dodge),
        Guide.xlabel("β"),Guide.ylabel("Density"),
        Guide.title("ρ=$(n2s(VR.MC.ρ)), T=$(n2s(VR.MC.T)), " *
          "μ(β)=$(n2s(mean(VR.βs[:,2]))), σ(β)=$(n2s(std(VR.βs[:,2])))")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "SEs")
  push!(plots,Vector{Plot}())
  for VR::VARReg ∈ VRs
    push!(plots[end],
      plot(x=VR.SEs[:,2], Geom.histogram(bincount=20, density=true, position=:dodge),
        Guide.xlabel("SE"),Guide.ylabel("Density"),
        Guide.title("ρ=$(n2s(VR.MC.ρ)), T=$(n2s(VR.MC.T)), " *
          "μ(SE)=$(n2s(mean(VR.SEs[:,2]))), σ(SE)=$(n2s(std(VR.SEs[:,2])))")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "Ts")
  push!(plots,Vector{Plot}())
  for VR::VARReg ∈ VRs
    push!(plots[end],
      plot(x=VR.Ts[:,2], Geom.histogram(bincount=20, density=true, position=:dodge),
        Guide.xlabel("t-stat (Newey-West Lag=5)"),Guide.ylabel("Density"),
        Guide.title("ρ=$(n2s(VR.MC.ρ)), T=$(n2s(VR.MC.T)), " *
          "μ(t)=$(n2s(mean(VR.Ts[:,2]))), σ(t)=$(n2s(std(VR.Ts[:,2])))")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
  end

end

#########################getNeweyWest ##########OLS Only
function getNeweyWest!(X::Matrix{Float64}, xqr::CTQR, ε::Vector{Float64}, lag::Int,
  Σ::Matrix{Float64} = Matrix{Float64}(xqr.K, xqr.K))::Matrix{Float64}

  #pre-allocate for the spectral matrix

  Rv::Matrix{Float64} = Matrix{Float64}(xqr.K, xqr.K) #pre-allocate working matrix
  RRInv::Matrix{Float64} = BLAS.gemm('N', 'T', xqr.RInv, xqr.RInv) #this is equivelent to [X'X]^-1

  #need to multiply through by the error
  Xe::Matrix{Float64} = X .* ε
  ST::Matrix{Float64} = BLAS.gemm('T','N',1.0/xqr.N, Xe, Xe)
  for v::Int ∈ 1:lag
    #overwrites Rv with (1/N)R'R
    BLAS.gemm!('T', 'N', 1.0/xqr.N, view(Xe, (v+1):(xqr.N), :),view(Xe, 1:(xqr.N-v), :), 0.0, Rv)
    #Rv .= view(Xe, (v+1):(xqr.N), :)' * view(Xe, 1:(xqr.N-v), :) .* (1.0/xqr.N)
    ST .+= (lag+1.-v)/(lag+1.) .* (Rv .+ Rv')
  end

  #this is [X'X]^-1S=[R'R]^-1S
  RRInvS::Matrix{Float64} = BLAS.gemm('N', 'N', RRInv, ST)

  #finally we have T[X'X]^-1S[X'X]^-1
  BLAS.gemm!('N','N',Float64(xqr.N), RRInvS, RRInv, 0.0, Σ)
  return Σ
end

function getNeweyWest!(lin::CTLM, lag::Int,
  Σ::Matrix{Float64} = Matrix{Float64}(lin.K, lin.K))::Matrix{Float64}

  return getNeweyWest!(lin.X, lin.xqr, lin.ε, lag, Σ)
end

#helper function in case the input requires a single argument function
getNeweyWest!(lag::Int) = (lin::CTLM)-> getNeweyWest!(lin, lag)

#test function to work
function testSalil()::Void
  Y::Vector{Float64} = SALIL_Y
  X::Matrix{Float64} = ones(length(SALIL_Z),2)
  X[:,2] .= SALIL_Z
  reg::CTLM = CTLM(X, Y)

  println("βs: $(reg.β)")
  println("σs: $((diag(getNeweyWest!(reg, 5))).^0.5)")

  return nothing
end

function A2P1Script(refreshData::Bool = false)::Void
  @time VRs::Vector{VARReg} = [VARReg(N=10_000, T=50, ρ=1.0),
      VARReg(N=10_000, T=50, ρ=0.9),
      VARReg(N=10_000, T=50, ρ=0.5),
      VARReg(N=10_000, T=1000, ρ=1.0),
      VARReg(N=10_000, T=1000, ρ=0.9),
      VARReg(N=10_000, T=1000, ρ=0.5)]

  visualizeA2P1(VRs)

  return nothing
end

#uncomment to run
@time begin
  A2P1Script()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using ArgParse
using CUDAapi
using Dates
using LinearAlgebra
using Logging
using MPI
using Printf
using Requires

using ..Atmos
using ..VTK
using ..ColumnwiseLUSolver
using ..Diagnostics
using ..GenericCallbacks
using ..ODESolvers
using ..Mesh.Grids: EveryDirection, VerticalDirection, HorizontalDirection
using ..MoistThermodynamics
using ..MPIStateArrays
using ..DGmethods: vars_state, vars_aux, update_aux!, update_aux_diffusive!
using ..TicToc
using ..VariableTemplates

Base.@kwdef mutable struct CLIMA_Settings
    disable_gpu::Bool = false
    mpi_knows_cuda::Bool = false
    show_updates::Bool = true
    update_interval::Int = 60
    enable_diagnostics::Bool = true
    diagnostics_interval::Int = 10000
    enable_vtk::Bool = false
    vtk_interval::Int = 10000
    log_level::String = "INFO"
    output_dir::String = "output"
    integration_testing::Bool = false
    array_type
end

const Settings = CLIMA_Settings(array_type = Array)

array_type() = Settings.array_type

include("Configurations.jl")

const cuarray_pkgid = Base.PkgId(Base.UUID("3a865a2d-5b23-5a0f-bc46-62713ec82fae"), "CuArrays")

@init begin
    if get(ENV, "CLIMA_GPU", "") != "false" && CUDAapi.has_cuda_gpu()
        CuArrays = Base.require(cuarray_pkgid)
    end
end

function _init_array(::Type{Array})
    return nothing
end

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    using .CuArrays, .CuArrays.CUDAdrv, .CuArrays.CUDAnative
    function _init_array(::Type{CuArray})
        comm = MPI.COMM_WORLD
        # allocate GPUs among MPI ranks
        local_comm = MPI.Comm_split_type(comm, MPI.MPI_COMM_TYPE_SHARED,  MPI.Comm_rank(comm))
        # we intentionally oversubscribe GPUs for testing: may want to disable this for production
        CUDAnative.device!(MPI.Comm_rank(local_comm) % length(devices()))
        CuArrays.allowscalar(false)
        return nothing
    end
end

function gpu_allowscalar(val)
    if haskey(Base.loaded_modules, CLIMA.cuarray_pkgid)
        Base.loaded_modules[CLIMA.cuarray_pkgid].allowscalar(val)
    end
    return nothing
end

"""
    parse_commandline()
"""
function parse_commandline()
    exc_handler = isinteractive() ? ArgParse.debug_handler : ArgParse.default_handler
    s = ArgParseSettings(exc_handler=exc_handler)

    @add_arg_table s begin
        "--disable-gpu"
            help = "do not use the GPU"
            action = :store_true
        "--mpi-knows-cuda"
            help = "MPI is CUDA-enabled"
            action = :store_true
        "--no-show-updates"
            help = "do not show simulation updates"
            action = :store_true
        "--update-interval"
            help = "interval in seconds for showing simulation updates"
            arg_type = Int
            default = 60
        "--disable-diagnostics"
            help = "disable the collection of diagnostics to <output-dir>"
            action = :store_true
        "--diagnostics-interval"
            help = "interval in simulation steps for gathering diagnostics"
            arg_type = Int
            default = 10000
        "--enable-vtk"
            help = "output VTK to <output-dir> every <vtk-interval> simulation steps"
            action = :store_true
        "--vtk-interval"
            help = "interval in simulation steps for VTK output"
            arg_type = Int
            default = 10000
        "--log-level"
            help = "set the log level to one of debug/info/warn/error"
            arg_type = String
            default = "info"
        "--output-dir"
            help = "directory for output data"
            arg_type = String
            default = "output"
        "--integration-testing"
            help = "enable integration testing"
            action = :store_true
    end

    return parse_args(s)
end

"""
    CLIMA.init(; disable_gpu=false)

Initialize MPI, allocate GPUs among MPI ranks if using GPUs, parse command
line arguments for CLIMA, and return a Dict of any additional command line
arguments.
"""
function init(; disable_gpu=false)
    # initialize MPI
    if !MPI.Initialized()
        MPI.Init()
    end

    # set up timing mechanism
    tictoc()

    # parse command line arguments
    try
        parsed_args = parse_commandline()
        Settings.disable_gpu = disable_gpu || parsed_args["disable-gpu"]
        Settings.mpi_knows_cuda = parsed_args["mpi-knows-cuda"]
        Settings.show_updates = !parsed_args["no-show-updates"]
        Settings.update_interval = parsed_args["update-interval"]
        Settings.enable_diagnostics = !parsed_args["disable-diagnostics"]
        Settings.diagnostics_interval = parsed_args["diagnostics-interval"]
        Settings.enable_vtk = parsed_args["enable-vtk"]
        Settings.vtk_interval = parsed_args["vtk-interval"]
        Settings.output_dir = parsed_args["output-dir"]
        Settings.integration_testing = parsed_args["integration-testing"]
        Settings.log_level = uppercase(parsed_args["log-level"])
    catch
        Settings.disable_gpu = disable_gpu
    end

    # set up the array type appropriately depending on whether we're using GPUs
    if !Settings.disable_gpu && get(ENV, "CLIMA_GPU", "") != "false" && CUDAapi.has_cuda_gpu()
        atyp = CuArrays.CuArray
    else
        atyp = Array
    end
    _init_array(atyp)
    Settings.array_type = atyp

    # create the output directory if needed
    if Settings.enable_diagnostics || Settings.enable_vtk
        mkpath(Settings.output_dir)
    end

    # set up logging
    loglevel = Settings.log_level == "DEBUG" ? Logging.Debug :
        Settings.log_level == "WARN"  ? Logging.Warn  :
        Settings.log_level == "ERROR" ? Logging.Error : Logging.Info
    # TODO: write a better MPI logging back-end and also integrate Dlog for large scale
    logger_stream = MPI.Comm_rank(MPI.COMM_WORLD) == 0 ? stderr : devnull
    global_logger(ConsoleLogger(logger_stream, loglevel))

    return nothing
end

"""
    CLIMA.SolverConfiguration

Parameters needed by `CLIMA.solve!()` to run a simulation.
"""
struct SolverConfiguration{FT}
    name::String
    mpicomm::MPI.Comm
    dg::DGModel
    Q::MPIStateArray
    t0::FT
    timeend::FT
    dt::FT
    forcecpu::Bool
    numberofsteps::Int
    init_args
    solver
end

"""
    CLIMA.setup_solver(t0, timeend, driver_config)

Set up the DG model per the specified driver configuration and set up the ODE solver.
"""
function setup_solver(t0::FT, timeend::FT,
                      driver_config::DriverConfiguration,
                      init_args...;
                      forcecpu=false,
                      ode_solver_type=nothing,
                      ode_dt=nothing,
                      modeldata=nothing,
                      Courant_number=0.4
                     ) where {FT<:AbstractFloat}
    @tic setup_solver

    bl = driver_config.bl
    grid = driver_config.grid
    numfluxnondiff = driver_config.numfluxnondiff
    numfluxdiff = driver_config.numfluxdiff
    gradnumflux = driver_config.gradnumflux

    # create DG model, initialize ODE state
    dg = DGModel(bl, grid, numfluxnondiff, numfluxdiff, gradnumflux,
                 modeldata=modeldata)
    @info @sprintf("Initializing %s", driver_config.name)
    Q = init_ode_state(dg, FT(0), init_args...; forcecpu=forcecpu)

    # TODO: using `update_aux!()` to apply filters to the state variables and
    # `update_aux_diffusive!()` to calculate the vertical component of velocity
    # using the integral kernels -- these should have their own interface
    update_aux!(dg, bl, Q, FT(0))
    update_aux_diffusive!(dg, bl, Q, FT(0))

    # if solver has been specified, use it
    if ode_solver_type !== nothing
        solver_type = ode_solver_type
    else
        solver_type = driver_config.solver_type
    end

    # create the linear model for IMEX solvers
    linmodel = nothing
    if isa(solver_type, ExplicitSolverType)
        dtmodel = bl
    else # solver_type === IMEXSolverType
        linmodel = solver_type.linear_model(bl)
        dtmodel = linmodel
    end

    # initial Δt specified or computed
    if ode_dt !== nothing
        dt = ode_dt
    else
        dt = calculate_dt(grid, dtmodel, Courant_number)
    end
    numberofsteps = convert(Int, cld(timeend, dt))
    dt = timeend / numberofsteps

    # create the solver
    if isa(solver_type, ExplicitSolverType)
        solver = solver_type.solver_method(dg, Q; dt=dt, t0=t0)
    else # solver_type === IMEXSolverType
        vdg = DGModel(linmodel, grid, numfluxnondiff, numfluxdiff, gradnumflux,
                      auxstate=dg.auxstate, direction=VerticalDirection())

        solver = solver_type.solver_method(dg, vdg, solver_type.linear_solver(), Q;
                                           dt=dt, t0=t0)
    end

    @toc setup_solver

    return SolverConfiguration(driver_config.name, driver_config.mpicomm, dg, Q,
                               t0, timeend, dt, forcecpu, numberofsteps,
                               init_args, solver)
end

"""
    CLIMA.invoke!(solver_config)

Run the simulation.
"""
function invoke!(solver_config::SolverConfiguration;
                 user_callbacks=(),
                 check_euclidean_distance=false,
                 adjustfinalstep=false
                )
    mpicomm = solver_config.mpicomm
    dg = solver_config.dg
    bl = dg.balancelaw
    Q = solver_config.Q
    FT = eltype(Q)
    timeend = solver_config.timeend
    forcecpu = solver_config.forcecpu
    init_args = solver_config.init_args
    solver = solver_config.solver

    # set up callbacks
    callbacks = ()
    if Settings.show_updates
        # set up the information callback
        upd_starttime = Ref(now())
        cbinfo = GenericCallbacks.EveryXWallTimeSeconds(Settings.update_interval, mpicomm) do (init=false)
            if init
                upd_starttime[] = now()
            else
                runtime = Dates.format(convert(Dates.DateTime,
                                               Dates.now()-upd_starttime[]),
                                       Dates.dateformat"HH:MM:SS")
                energy = norm(solver_config.Q)
                @info @sprintf("""Update
                               simtime = %.16e
                               runtime = %s
                               norm(Q) = %.16e""",
                               ODESolvers.gettime(solver),
                               runtime,
                               energy)
            end
            nothing
        end
        callbacks = (callbacks..., cbinfo)
    end
    if Settings.enable_diagnostics
        # set up diagnostics to be collected via callback
        cbdiagnostics = GenericCallbacks.EveryXSimulationSteps(Settings.diagnostics_interval) do (init=false)
            if init
                dia_starttime = replace(string(now()), ":" => ".")
                Diagnostics.init(mpicomm, dg, Q, dia_starttime, Settings.output_dir)
            end
            currtime = ODESolvers.gettime(solver)
            @info @sprintf("""Diagnostics
                           collecting at %s""",
                           string(currtime))
            Diagnostics.collect(currtime)
            nothing
        end
        callbacks = (callbacks..., cbdiagnostics)
    end
    if Settings.enable_vtk
        # set up VTK output callback
        step = [0]
        cbvtk = GenericCallbacks.EveryXSimulationSteps(Settings.vtk_interval) do (init=false)
            vprefix = @sprintf("%s_mpirank%04d_step%04d", solver_config.name,
                               MPI.Comm_rank(mpicomm), step[1])
            outprefix = joinpath(Settings.output_dir, vprefix)
            statenames = flattenednames(vars_state(bl, FT))
            auxnames = flattenednames(vars_aux(bl, FT))
            writevtk(outprefix, Q, dg, statenames, dg.auxstate, auxnames)
            # Generate the pvtu file for these vtk files
            if MPI.Comm_rank(mpicomm) == 0
                # name of the pvtu file
                pprefix = @sprintf("%s_step%04d", solver_config.name, step[1])
                pvtuprefix = joinpath(Settings.output_dir, pprefix)
                # name of each of the ranks vtk files
                prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
                    @sprintf("%s_mpirank%04d_step%04d", solver_config.name, i-1, step[1])
                end
                writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...))
            end
            step[1] += 1
            nothing
        end
        callbacks = (callbacks..., cbvtk)
    end
    callbacks = (callbacks..., user_callbacks...)

    # initial condition norm
    eng0 = norm(Q)
    @info @sprintf("""Starting %s
                   dt              = %.5e
                   timeend         = %.5e
                   number of steps = %d
                   norm(Q)         = %.16e""",
                   solver_config.name,
                   solver_config.dt,
                   solver_config.timeend,
                   solver_config.numberofsteps,
                   eng0)

    # run the simulation
    @tic solve!
    solve!(Q, solver; timeend=timeend, callbacks=callbacks, adjustfinalstep=adjustfinalstep)
    @toc solve!

    engf = norm(solver_config.Q)

    @info @sprintf("""Finished
                   norm(Q)            = %.16e
                   norm(Q) / norm(Q₀) = %.16e
                   norm(Q) - norm(Q₀) = %.16e""",
                   engf,
                   engf/eng0,
                   engf-eng0)

    if check_euclidean_distance
        Qe = init_ode_state(dg, timeend, init_args...; forcecpu=forcecpu)
        engfe = norm(Qe)
        errf = euclidean_distance(solver_config.Q, Qe)
        @info @sprintf("""Euclidean distance
                       norm(Q - Qe)            = %.16e
                       norm(Q - Qe) / norm(Qe) = %.16e""",
                       errf,
                       errf/engfe)
    end

    return engf / eng0
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
# doc.jl - infrastructure

"""
    process(::Array{AppliSales.Order, 1}; path="./test_invoices.sqlite")

    process(::Array{UnpaidInvoice, 1}, ::Array{AppliGeneralLedger.BankStatement, 1}; path="./test_invoices.sqlite")

- Creates UnpaidInvoice's from AppliSale.Order's, archive them, and creates AppliGeneralLedger.Entry's for the general ledger.
- Creates PaidInvoices's from UnpaidInvoices by using AppliGeneralLedger.BankStatement's, and creates AppliGeneralLedger.Entry's for the general ledger.

# Example

```
julia> using AppliSales

julia> using AppliGeneralLedger

julia> using AppliInvoicing

julia> const PATH_CSV = "./bank.csv"

julia> orders = AppliSales.process()

julia> journal_entries_1 = AppliInvoicing.process(orders)

julia> stms = AppliInvoicing.read_bank_statements(PATH_CSV)

julia> unpaid_invoices = retrieve_unpaid_invoices()

julia> journal_entries_2 = AppliInvoicing.process(unpaid_invoices, stms)

julia> cmd = `rm test_invoicing.sqlite`

julia> run(cmd)
```
"""
function process end

"""
    read_bank_statements(path::String)

Retrieves bank statements from a CSV-file.

# Example

```
julia> const PATH_CSV = "./bank.csv"

julia> stms = AppliInvoicing.read_bank_statements(PATH_CSV)
```
"""
function read_bank_statements end

"""
    retrieve_unpaid_invoices(path="./test_invoicing.sqlite")::Array{UnpaidInvoice, 1}

Retrieves UnpaidInvoice's from a SQLite.jl database.

# Example

```
julia> using AppliSales

julia> using AppliInvoicing

julia> orders = AppliSales.process()

julia> AppliInvoicing.process(orders)

julia> unpaid_invoices = retrieve_unpaid_invoices()

julia> cmd = `rm test_invoicing.sqlite`

julia> run(cmd)
```
"""
function retrieve_unpaid_invoices end

"""
    retrieve_paid_invoices(path="./test_invoicing.sqlite")::Array{PaidInvoice, 1}

Retrieves PaidInvoice's from a SQLite.jl database.

# Example

```
TODO
```
"""
function retrieve_paid_invoices end


text:
# Model not supporting Interval
MOIU.@model(SimpleModel,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.RotatedSecondOrderCone, MOI.GeometricMeanCone,
             MOI.PositiveSemidefiniteConeTriangle, MOI.ExponentialCone),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction,))

function test_noc(bridgedmock, F, S, n)
    @test MOI.get(bridgedmock, MOI.NumberOfConstraints{F, S}()) == n
    @test length(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{F, S}())) == n
    @test ((F, S) in MOI.get(bridgedmock, MOI.ListOfConstraints())) == !iszero(n)
end

# Test deletion of bridge
function test_delete_bridge(m::MOIB.AbstractBridgeOptimizer, ci::MOI.ConstraintIndex{F, S}, nvars::Int, nocs::Tuple) where {F, S}
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    test_noc(m, F, S, 1)
    for noc in nocs
        test_noc(m, noc...)
    end
    @test MOI.is_valid(m, ci)
    MOI.delete(m, ci)
    @test_throws MOI.InvalidIndex{typeof(ci)} MOI.delete(m, ci)
    try
        MOI.delete(m, ci)
    catch err
        @test err.index == ci
    end
    @test !MOI.is_valid(m, ci)
    @test isempty(m.bridges)
    test_noc(m, F, S, 0)
    # As the bridge has been removed, if the constraints it has created where not removed, it wouldn't be there to decrease this counter anymore
    @test MOI.get(m, MOI.NumberOfVariables()) == nvars
    for noc in nocs
        test_noc(m, noc...)
    end
end

@testset "BridgeOptimizer" begin
    mock = MOIU.MockOptimizer(SimpleModel{Float64}())
    bridgedmock = MOIB.SplitInterval{Float64}(mock)

    @testset "Issue #453" begin
        MOI.empty!(bridgedmock)
        MOIU.loadfromstring!(bridgedmock, """
            variables: x
            maxobjective: 3.0x
            c: 2.0x in Interval(1.0, 4.0)
            d: x in LessThan(1.5)
        """)
        x = MOI.get(bridgedmock, MOI.VariableIndex, "x")
        @test isa(x, MOI.VariableIndex)
        c1 = MOI.get(bridgedmock, MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}, "c")
        @test isa(c1, MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}})
        c2 = MOI.get(bridgedmock, MOI.ConstraintIndex, "c")
        @test c1 == c2
        d1 = MOI.get(bridgedmock, MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}, "d")
        @test isa(d1, MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}})
        d2 = MOI.get(bridgedmock, MOI.ConstraintIndex, "d")
        @test d1 == d2
    end

    MOI.empty!(bridgedmock)

    @testset "Name test" begin
        MOIT.nametest(bridgedmock)
    end

    @testset "Copy test" begin
        MOIT.failcopytestc(bridgedmock)
        MOIT.failcopytestia(bridgedmock)
        MOIT.failcopytestva(bridgedmock)
        MOIT.failcopytestca(bridgedmock)
        MOIT.copytest(bridgedmock, SimpleModel{Float64}())
    end

    @testset "Custom test" begin
        model = MOIB.SplitInterval{Int}(SimpleModel{Int}())
        @test !MOIB.supports_bridging_constraint(model, MOI.VectorAffineFunction{Float64}, MOI.Interval{Float64})

        x, y = MOI.add_variables(model, 2)
        @test MOI.get(model, MOI.NumberOfVariables()) == 2

        f1 = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(3, x)], 7)
        c1 = MOI.add_constraint(model, f1, MOI.Interval(-1, 1))

        @test MOI.get(model, MOI.ListOfConstraints()) == [(MOI.ScalarAffineFunction{Int},MOI.Interval{Int})]
        test_noc(model, MOI.ScalarAffineFunction{Int}, MOI.GreaterThan{Int}, 0)
        test_noc(model, MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}, 1)
        @test (@inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Int},MOI.Interval{Int}}())) == [c1]

        f2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2, -1], [x, y]), 2)
        c2 = MOI.add_constraint(model, f1, MOI.GreaterThan(-2))

        @test MOI.get(model, MOI.ListOfConstraints()) == [(MOI.ScalarAffineFunction{Int},MOI.GreaterThan{Int}), (MOI.ScalarAffineFunction{Int},MOI.Interval{Int})]
        test_noc(model, MOI.ScalarAffineFunction{Int}, MOI.GreaterThan{Int}, 1)
        test_noc(model, MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}, 1)
        @test (@inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Int},MOI.Interval{Int}}())) == [c1]
        @test (@inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Int},MOI.GreaterThan{Int}}())) == [c2]

        @test MOI.is_valid(model, c2)
        MOI.delete(model, c2)

        @test MOI.get(model, MOI.ListOfConstraints()) == [(MOI.ScalarAffineFunction{Int},MOI.Interval{Int})]
        test_noc(model, MOI.ScalarAffineFunction{Int}, MOI.GreaterThan{Int}, 0)
        test_noc(model, MOI.ScalarAffineFunction{Int}, MOI.Interval{Int}, 1)
        @test (@inferred MOI.get(model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Int},MOI.Interval{Int}}())) == [c1]
    end

    @testset "Continuous Linear" begin
        MOIT.contlineartest(bridgedmock, MOIT.TestConfig(solve=false))
    end
end

# Model not supporting RotatedSecondOrderCone
MOIU.@model(NoRSOCModel,
            (),
            (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
            (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives, MOI.SecondOrderCone,
             MOI.ExponentialCone, MOI.PositiveSemidefiniteConeTriangle),
            (),
            (MOI.SingleVariable,),
            (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction),
            (MOI.VectorOfVariables,),
            (MOI.VectorAffineFunction, MOI.VectorQuadraticFunction))

@testset "LazyBridgeOptimizer" begin
    mock = MOIU.MockOptimizer(NoRSOCModel{Float64}())
    bridgedmock = MOIB.LazyBridgeOptimizer(mock, Model{Float64}())
    MOIB.add_bridge(bridgedmock, MOIB.SplitIntervalBridge{Float64})
    MOIB.add_bridge(bridgedmock, MOIB.RSOCtoPSDBridge{Float64})
    MOIB.add_bridge(bridgedmock, MOIB.SOCtoPSDBridge{Float64})
    MOIB.add_bridge(bridgedmock, MOIB.RSOCBridge{Float64})

    @testset "Name test" begin
        MOIT.nametest(bridgedmock)
    end

    @testset "Copy test" begin
        MOIT.failcopytestc(bridgedmock)
        MOIT.failcopytestia(bridgedmock)
        MOIT.failcopytestva(bridgedmock)
        MOIT.failcopytestca(bridgedmock)
        MOIT.copytest(bridgedmock, NoRSOCModel{Float64}())
    end

    # Test that RSOCtoPSD is used instead of RSOC+SOCtoPSD as it is a shortest path.
    @testset "Bridge selection" begin
        MOI.empty!(bridgedmock)
        @test !(MOI.supports_constraint(bridgedmock, MOI.VectorAffineFunction{Float64}, MOI.LogDetConeTriangle))
        x = MOI.add_variables(bridgedmock, 3)
        c = MOI.add_constraint(bridgedmock, MOI.VectorOfVariables(x), MOI.RotatedSecondOrderCone(3))
        @test MOIB.bridge(bridgedmock, c) isa MOIB.RSOCtoPSDBridge
        @test bridgedmock.dist[(MathOptInterface.VectorOfVariables, MathOptInterface.RotatedSecondOrderCone)] == 1
    end

    @testset "Supports" begin
        fullbridgedmock = MOIB.fullbridgeoptimizer(mock, Float64)
        for F in [MOI.SingleVariable, MOI.ScalarAffineFunction{Float64},
                  MOI.ScalarQuadraticFunction{Float64}]
            @test MOI.supports_constraint(fullbridgedmock, F,
                                          MOI.Interval{Float64})
        end
        for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
                  MOI.VectorQuadraticFunction{Float64}]
            @test MOI.supports_constraint(fullbridgedmock, F,
                                          MOI.PositiveSemidefiniteConeSquare)
            @test MOI.supports_constraint(fullbridgedmock, F,
                                          MOI.GeometricMeanCone)
        end
        for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}]
            # The bridges in this for loop do not support yet
            # VectorQuadraticFunction. See TODO's for the reason.
            # TODO: Missing vcat for quadratic for supporting quadratic.
            @test MOI.supports_constraint(fullbridgedmock, F,
                                          MOI.RotatedSecondOrderCone)
            # TODO: Det bridges need to use MOIU.operate to support quadratic.
            @test MOI.supports_constraint(fullbridgedmock, F,
                                          MOI.LogDetConeTriangle)
            @test MOI.supports_constraint(fullbridgedmock, F,
                                          MOI.RootDetConeTriangle)
        end
    end

    @testset "Combining two briges" begin
        fullbridgedmock = MOIB.fullbridgeoptimizer(mock, Float64)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1, 1, 0, 1, √2])
        config = MOIT.TestConfig()
        MOIT.rootdett1vtest(fullbridgedmock, config)
        MOIT.rootdett1ftest(fullbridgedmock, config)
        # Dual is not yet implemented for RootDet and GeoMean bridges
        ci = first(MOI.get(fullbridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle}()))
        test_delete_bridge(fullbridgedmock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                    (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
                                                    (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)))

    end

    @testset "Continuous Linear" begin
        MOIT.contlineartest(bridgedmock, MOIT.TestConfig(solve=false))
    end

    @testset "Continuous Conic" begin
        MOIT.contconictest(MOIB.fullbridgeoptimizer(mock, Float64), MOIT.TestConfig(solve=false), ["logdets", "rootdets", "psds"])
    end
end

@testset "Bridge tests" begin
    mock = MOIU.MockOptimizer(SimpleModel{Float64}())
    config = MOIT.TestConfig()

    @testset "Interval" begin
        bridgedmock = MOIB.SplitInterval{Float64}(mock)
        MOIT.basic_constraint_tests(bridgedmock, config,
                                    include=[(MOI.SingleVariable,
                                              MOI.Interval{Float64}),
                                             (MOI.ScalarAffineFunction{Float64},
                                              MOI.Interval{Float64}),
                                             (MOI.ScalarQuadraticFunction{Float64},
                                              MOI.Interval{Float64})])
        MOIU.set_mock_optimize!(mock,
             (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [5.0, 5.0],
                  (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0],
                  (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [-1]),
             (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [2.5, 2.5],
                  (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [1],
                  (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64})    => [0]),
             (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0]),
             (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [6.0, 6.0]))
        MOIT.linear10test(bridgedmock, config)
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}}()))
        newf = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([1.0, -1.0], MOI.get(bridgedmock, MOI.ListOfVariableIndices())), 0.0)
        MOI.set(bridgedmock, MOI.ConstraintFunction(), ci, newf)
        @test MOI.get(bridgedmock, MOI.ConstraintFunction(), ci) ≈ newf
        test_delete_bridge(bridgedmock, ci, 2, ((MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
                                                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64},    0)))
   end

    @testset "RSOC" begin
        bridgedmock = MOIB.RSOC{Float64}(mock)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2, 0.5, 1.0],
                              (MOI.SingleVariable,                MOI.EqualTo{Float64}) => [-√2, -1/√2],
                              (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)  => [[3/2, 1/2, -1.0, -1.0]])
        MOIT.rotatedsoc1vtest(bridgedmock, config)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2],
                              (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)  => [[3/2, 1/2, -1.0, -1.0]])
        MOIT.rotatedsoc1ftest(bridgedmock, config)
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}()))
        test_delete_bridge(bridgedmock, ci, 2, ((MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone, 0),))
    end

    @testset "QuadtoSOC" begin
        bridgedmock = MOIB.QuadtoSOC{Float64}(mock)
        @testset "Error for non-convex quadratic constraints" begin
            x = MOI.add_variable(bridgedmock)
            @test_throws ErrorException begin
                MOI.add_constraint(bridgedmock,
                                   MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                                               [MOI.ScalarQuadraticTerm(1.0, x, x)],
                                                               0.0),
                                   MOI.GreaterThan(0.0))
            end
            @test_throws ErrorException begin
                MOI.add_constraint(bridgedmock,
                                   MOI.ScalarQuadraticFunction(MOI.ScalarAffineTerm{Float64}[],
                                                               [MOI.ScalarQuadraticTerm(-1.0, x, x)],
                                                               0.0),
                                   MOI.LessThan(0.0))
            end
        end
        @testset "Quadratic constraints with 2 variables" begin
            MOIU.set_mock_optimize!(mock,
                (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
                    MOI.Success,
                    (MOI.FeasiblePoint, [0.5, 0.5])
                ),
                (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock,
                    MOI.Success,
                    (MOI.FeasiblePoint, [0.5, (√13 - 1)/4])
                )
            )
            MOIT.solve_qcp_edge_cases(bridgedmock, config)
            ci = first(MOI.get(mock,
                               MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64},
                                                           MOI.RotatedSecondOrderCone}()))
            x, y = MOI.get(mock, MOI.ListOfVariableIndices())
            # The matrix is
            # 2 1
            # 1 2
            # for which the cholesky factorization is U' * U with U =
            # √2 √2/2
            #  . √3/√2
            expected = MOI.VectorAffineFunction{Float64}(MOI.VectorAffineTerm.([3, 3, 4],
                                                                               MOI.ScalarAffineTerm.([√2, √2/2, √3/√2],
                                                                                [x, y, y])),
                                                         [1.0, 1.0, 0.0, 0.0])
            @test MOI.get(mock, MOI.ConstraintFunction(), ci) ≈ expected
        end
        @testset "QCP tests" begin
            MOIU.set_mock_optimize!(mock,
                (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/2, 7/4], MOI.FeasiblePoint))
            MOIT.qcp1test(bridgedmock, config)
            MOIU.set_mock_optimize!(mock,
                (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FeasiblePoint))
            MOIT.qcp2test(bridgedmock, config)
            MOIU.set_mock_optimize!(mock,
                (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [√2], MOI.FeasiblePoint))
            MOIT.qcp3test(bridgedmock, config)
            @testset "Bridge deletion" begin
                ci = first(MOI.get(bridgedmock,
                                   MOI.ListOfConstraintIndices{MOI.ScalarQuadraticFunction{Float64},
                                                               MOI.LessThan{Float64}}()))
                test_delete_bridge(bridgedmock, ci, 1, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),))
            end
        end
    end

    @testset "SquarePSD" begin
        bridgedmock = MOIB.SquarePSD{Float64}(mock)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
                              (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [2, 2])
        MOIT.psds0vtest(bridgedmock, config)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, ones(4),
                              (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[1, -1, 1]],
                              (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})                 => [2, 2])
        MOIT.psds0ftest(bridgedmock, config)
        ci = first(MOI.get(bridgedmock,
                           MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64},
                                                       MOI.PositiveSemidefiniteConeSquare}()))
        test_delete_bridge(bridgedmock, ci, 4,
                           ((MOI.VectorAffineFunction{Float64},
                             MOI.PositiveSemidefiniteConeTriangle, 0),
                            (MOI.ScalarAffineFunction{Float64},
                             MOI.EqualTo{Float64}, 1)))
    end

    @testset "GeoMean" begin
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [ones(4); 2; √2; √2])
        bridgedmock = MOIB.GeoMean{Float64}(mock)
        MOIT.geomean1vtest(bridgedmock, config)
        MOIT.geomean1ftest(bridgedmock, config)
        # Dual is not yet implemented for GeoMean bridge
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone}()))
        test_delete_bridge(bridgedmock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64},      1)))
    end

    @testset "SOCtoPSD" begin
        bridgedmock = MOIB.SOCtoPSD{Float64}(mock)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1/√2, 1/√2],
                              (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[√2/2, -1/2, √2/4, -1/2, √2/4, √2/4]],
                              (MOI.VectorAffineFunction{Float64}, MOI.Zeros)                            => [[-√2]])
        MOIT.soc1vtest(bridgedmock, config)
        MOIT.soc1ftest(bridgedmock, config)
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone}()))
        test_delete_bridge(bridgedmock, ci, 3, ((MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0),))
    end

    @testset "RSOCtoPSD" begin
        bridgedmock = MOIB.RSOCtoPSD{Float64}(mock)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2, 0.5, 1.0],
                              (MOI.SingleVariable,                MOI.EqualTo{Float64})       => [-√2, -1/√2],
                              (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[√2, -1/2, √2/8, -1/2, √2/8, √2/8]])
        MOIT.rotatedsoc1vtest(bridgedmock, config)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2],
                              (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle) => [[√2, -1/2, √2/8, -1/2, √2/8, √2/8]])
        MOIT.rotatedsoc1ftest(bridgedmock, config)
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone}()))
        test_delete_bridge(bridgedmock, ci, 2, ((MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0),))
    end

    @testset "LogDet" begin
        bridgedmock = MOIB.LogDet{Float64}(mock)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 1, 1, 0, 1, 0, 0])
        MOIT.logdett1vtest(bridgedmock, config)
        MOIT.logdett1ftest(bridgedmock, config)
        # Dual is not yet implemented for LogDet bridge
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.LogDetConeTriangle}()))
        test_delete_bridge(bridgedmock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone, 0), (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)))
    end

    @testset "RootDet" begin
        bridgedmock = MOIB.RootDet{Float64}(mock)
        mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1, 1, 0, 1])
        MOIT.rootdett1vtest(bridgedmock, config)
        MOIT.rootdett1ftest(bridgedmock, config)
        # Dual is not yet implemented for RootDet bridge
        ci = first(MOI.get(bridgedmock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle}()))
        test_delete_bridge(bridgedmock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                                (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
                                                (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)))
    end
end

# We need to test this in a module at the top level because it can't be defined
# in a testset. If it runs without error, then we're okay.
module TestExternalBridge
    using MathOptInterface

    struct StrictlyGreaterThan <: MathOptInterface.AbstractScalarSet end
    struct StrictlyGreaterBridge{T} <: MathOptInterface.Bridges.AbstractBridge end

    function StrictlyGreaterBridge(
            model,
            func::MathOptInterface.SingleVariable,
            set::StrictlyGreaterThan)
        return StrictlyGreaterBridge{Float64}()
    end

    function MathOptInterface.supports_constraint(
            ::Type{StrictlyGreaterBridge{T}},
            ::Type{MathOptInterface.SingleVariable},
            ::Type{StrictlyGreaterThan}) where {T}
        return true
    end

    function MathOptInterface.Bridges.added_constraint_types(
            ::Type{StrictlyGreaterBridge{T}},
            ::Type{MathOptInterface.SingleVariable},
            ::Type{StrictlyGreaterThan}) where {T}
        return [(
            MathOptInterface.SingleVariable,
            MathOptInterface.GreaterThan{T}
        )]
    end

    MathOptInterface.Bridges.@bridge(StrictlyGreater, StrictlyGreaterBridge,
        (StrictlyGreaterThan, ),
        (),
        (),
        (),
        (MathOptInterface.SingleVariable, ),
        (),
        (),
        ()
    )
end

@testset "@bridge with external components" begin
    model = SimpleModel{Float64}();
    @test MOI.supports_constraint(model, MOI.SingleVariable, MOI.GreaterThan{Float64})
    @test !MOI.supports_constraint(model, MOI.SingleVariable, TestExternalBridge.StrictlyGreaterThan)

    bridge = TestExternalBridge.StrictlyGreater{Float64}(model);
    @test MOI.supports_constraint(bridge, MOI.SingleVariable, MOI.GreaterThan{Float64})
    @test MOI.supports_constraint(bridge, MOI.SingleVariable, TestExternalBridge.StrictlyGreaterThan)
end


text:
n = [3.36, 3.692, 2.05, 22.851, 0.389, 0.351, 0.919, 0.813, 0.055, 6.6];
d =  [1.136, 13.6192, 12.05, 2.851, 9, 10.351, 0.111919, 10.813, 1.055, 1.6];
leftIdx = 1:4
rightIdx = 5:length(n)
nl = n[leftIdx]
nr = n[rightIdx]
dl = d[leftIdx]
dr = d[rightIdx]

w = copy(d)
wl = w[leftIdx]
wr = w[rightIdx]

include(joinpath("testing", "customVariance_fns.jl"))

zzr = CustomVariance()
merge!(cv,zzr)

cv = CustomVariance(n, d)
cv.meanratio
cv.m2
cv.var
cvl = CustomVariance(nl, dl)
cvr = CustomVariance(nr, dr)

merge!(cvl,cvr)
@code_warntype merge!(cvl, cvr)

@code_warntype unmerge!(cvl, cvr)
unmerge!(cvl,cvr)
cvl
CustomVariance(nl,dl)
# var want 53.3375356025
# m2 want 160.0126068077387

unmerge!(cvl,CustomVariance(nl, dl))
cvl
CustomVariance(nr,dr)

sum(n) / sum(d)
sum(d .* (n ./ d - sum(n) / sum(d)).^2)
sum(nl) / sum(dl)
sum(nr) / sum(dr)
0

CustomVariance()

cvl = CustomVariance(nl, dl)
merge!(cvl,CustomVariance())
for i in [5,6]
    addSinglePoint!(cvl, n[i], d[i])
    removeSinglePoint!(cvl, n[i], d[i])
end
cvl

CustomVariance(n[1:6],d[1:6])


text:
using ECCodes
using PyPlot

withmessages("gfs.t18z.pgrb2.0p50.f015") do msg
    if (msg["typeOfLevel"]=="heightAboveGround")&&
        (msg["shortName"]=="10v")&&
        (msg["level"]=="10")&&
        (msg["dataDate"]=="20170709")&&
        (msg["dataTime"]=="1800")&&
        (msg["forecastTime"]=="15")
        global lats,lons,vals
        lats,lons,vals=ECCodes.data(msg)
    end
end

PyPlot.tricontourf(lons,lats,vals,100)