text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: resource_groups_tagging_api
using AWS.Compat
using AWS.UUIDs
"""
    DescribeReportCreation()

Describes the status of the StartReportCreation operation.  You can call this operation only from the organization's master account and from the us-east-1 Region.

"""

describe_report_creation(; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("DescribeReportCreation"; aws_config=aws_config)
describe_report_creation(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("DescribeReportCreation", args; aws_config=aws_config)

"""
    GetComplianceSummary()

Returns a table that shows counts of resources that are noncompliant with their tag policies. For more information on tag policies, see Tag Policies in the AWS Organizations User Guide.  You can call this operation only from the organization's master account and from the us-east-1 Region.

# Optional Parameters
- `GroupBy`: A list of attributes to group the counts of noncompliant resources by. If supplied, the counts are sorted by those attributes.
- `MaxResults`: A limit that restricts the number of results that are returned per page.
- `PaginationToken`: A string that indicates that additional data is available. Leave this value empty for your initial request. If the response includes a PaginationToken, use that string for this value to request an additional page of data.
- `RegionFilters`: A list of Regions to limit the output by. If you use this parameter, the count of returned noncompliant resources includes only resources in the specified Regions.
- `ResourceTypeFilters`: The constraints on the resources that you want returned. The format of each resource type is service[:resourceType]. For example, specifying a resource type of ec2 returns all Amazon EC2 resources (which includes EC2 instances). Specifying a resource type of ec2:instance returns only EC2 instances.  The string for each service name and resource type is the same as that embedded in a resource's Amazon Resource Name (ARN). Consult the AWS General Reference for the following:   For a list of service name strings, see AWS Service Namespaces.   For resource type strings, see Example ARNs.   For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.   You can specify multiple resource types by using an array. The array can include up to 100 items. Note that the length constraint requirement applies to each resource type filter.
- `TagKeyFilters`: A list of tag keys to limit the output by. If you use this parameter, the count of returned noncompliant resources includes only resources that have the specified tag keys.
- `TargetIdFilters`: The target identifiers (usually, specific account IDs) to limit the output by. If you use this parameter, the count of returned noncompliant resources includes only resources with the specified target IDs.
"""

get_compliance_summary(; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetComplianceSummary"; aws_config=aws_config)
get_compliance_summary(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetComplianceSummary", args; aws_config=aws_config)

"""
    GetResources()

Returns all the tagged or previously tagged resources that are located in the specified Region for the AWS account. Depending on what information you want returned, you can also specify the following:    Filters that specify what tags and resource types you want returned. The response includes all tags that are associated with the requested resources.   Information about compliance with the account's effective tag policy. For more information on tag policies, see Tag Policies in the AWS Organizations User Guide.     You can check the PaginationToken response parameter to determine if a query is complete. Queries occasionally return fewer results on a page than allowed. The PaginationToken response parameter value is null only when there are no more results to display.

# Optional Parameters
- `ExcludeCompliantResources`: Specifies whether to exclude resources that are compliant with the tag policy. Set this to true if you are interested in retrieving information on noncompliant resources only. You can use this parameter only if the IncludeComplianceDetails parameter is also set to true.
- `IncludeComplianceDetails`: Specifies whether to include details regarding the compliance with the effective tag policy. Set this to true to determine whether resources are compliant with the tag policy and to get details.
- `PaginationToken`: A string that indicates that additional data is available. Leave this value empty for your initial request. If the response includes a PaginationToken, use that string for this value to request an additional page of data.
- `ResourceTypeFilters`: The constraints on the resources that you want returned. The format of each resource type is service[:resourceType]. For example, specifying a resource type of ec2 returns all Amazon EC2 resources (which includes EC2 instances). Specifying a resource type of ec2:instance returns only EC2 instances.  The string for each service name and resource type is the same as that embedded in a resource's Amazon Resource Name (ARN). Consult the AWS General Reference for the following:   For a list of service name strings, see AWS Service Namespaces.   For resource type strings, see Example ARNs.   For more information about ARNs, see Amazon Resource Names (ARNs) and AWS Service Namespaces.   You can specify multiple resource types by using an array. The array can include up to 100 items. Note that the length constraint requirement applies to each resource type filter.
- `ResourcesPerPage`: A limit that restricts the number of resources returned by GetResources in paginated output. You can set ResourcesPerPage to a minimum of 1 item and the maximum of 100 items.
- `TagFilters`: A list of TagFilters (keys and values). Each TagFilter specified must contain a key with values as optional. A request can include up to 50 keys, and each key can include up to 20 values.  Note the following when deciding how to use TagFilters:   If you do specify a TagFilter, the response returns only those resources that are currently associated with the specified tag.    If you don't specify a TagFilter, the response includes all resources that were ever associated with tags. Resources that currently don't have associated tags are shown with an empty tag set, like this: \"Tags\": [].   If you specify more than one filter in a single request, the response returns only those resources that satisfy all specified filters.   If you specify a filter that contains more than one value for a key, the response returns resources that match any of the specified values for that key.   If you don't specify any values for a key, the response returns resources that are tagged with that key irrespective of the value. For example, for filters: filter1 = {key1, {value1}}, filter2 = {key2, {value2,value3,value4}} , filter3 = {key3}:   GetResources( {filter1} ) returns resources tagged with key1=value1   GetResources( {filter2} ) returns resources tagged with key2=value2 or key2=value3 or key2=value4   GetResources( {filter3} ) returns resources tagged with any tag containing key3 as its tag key, irrespective of its value   GetResources( {filter1,filter2,filter3} ) returns resources tagged with ( key1=value1) and ( key2=value2 or key2=value3 or key2=value4) and (key3, irrespective of the value)
- `TagsPerPage`: AWS recommends using ResourcesPerPage instead of this parameter. A limit that restricts the number of tags (key and value pairs) returned by GetResources in paginated output. A resource with no tags is counted as having one tag (one key and value pair).  GetResources does not split a resource and its associated tags across pages. If the specified TagsPerPage would cause such a break, a PaginationToken is returned in place of the affected resource and its tags. Use that token in another request to get the remaining data. For example, if you specify a TagsPerPage of 100 and the account has 22 resources with 10 tags each (meaning that each resource has 10 key and value pairs), the output will consist of three pages. The first page displays the first 10 resources, each with its 10 tags. The second page displays the next 10 resources, each with its 10 tags. The third page displays the remaining 2 resources, each with its 10 tags. You can set TagsPerPage to a minimum of 100 items and the maximum of 500 items.
"""

get_resources(; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetResources"; aws_config=aws_config)
get_resources(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetResources", args; aws_config=aws_config)

"""
    GetTagKeys()

Returns all tag keys in the specified Region for the AWS account.

# Optional Parameters
- `PaginationToken`: A string that indicates that additional data is available. Leave this value empty for your initial request. If the response includes a PaginationToken, use that string for this value to request an additional page of data.
"""

get_tag_keys(; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetTagKeys"; aws_config=aws_config)
get_tag_keys(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetTagKeys", args; aws_config=aws_config)

"""
    GetTagValues()

Returns all tag values for the specified key in the specified Region for the AWS account.

# Required Parameters
- `Key`: The key for which you want to list all existing values in the specified Region for the AWS account.

# Optional Parameters
- `PaginationToken`: A string that indicates that additional data is available. Leave this value empty for your initial request. If the response includes a PaginationToken, use that string for this value to request an additional page of data.
"""

get_tag_values(Key; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetTagValues", Dict{String, Any}("Key"=>Key); aws_config=aws_config)
get_tag_values(Key, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("GetTagValues", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Key"=>Key), args)); aws_config=aws_config)

"""
    StartReportCreation()

Generates a report that lists all tagged resources in accounts across your organization and tells whether each resource is compliant with the effective tag policy. Compliance data is refreshed daily.  The generated report is saved to the following location:  s3://example-bucket/AwsTagPolicies/o-exampleorgid/YYYY-MM-ddTHH:mm:ssZ/report.csv  You can call this operation only from the organization's master account and from the us-east-1 Region.

# Required Parameters
- `S3Bucket`: The name of the Amazon S3 bucket where the report will be stored; for example:  awsexamplebucket  For more information on S3 bucket requirements, including an example bucket policy, see the example S3 bucket policy on this page.

"""

start_report_creation(S3Bucket; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("StartReportCreation", Dict{String, Any}("S3Bucket"=>S3Bucket); aws_config=aws_config)
start_report_creation(S3Bucket, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("StartReportCreation", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("S3Bucket"=>S3Bucket), args)); aws_config=aws_config)

"""
    TagResources()

Applies one or more tags to the specified resources. Note the following:   Not all resources can have tags. For a list of services that support tagging, see this list.   Each resource can have up to 50 tags. For other limits, see Tag Naming and Usage Conventions in the AWS General Reference.    You can only tag resources that are located in the specified Region for the AWS account.   To add tags to a resource, you need the necessary permissions for the service that the resource belongs to as well as permissions for adding tags. For more information, see this list.    Do not store personally identifiable information (PII) or other confidential or sensitive information in tags. We use tags to provide you with billing and administration services. Tags are not intended to be used for private or sensitive data.

# Required Parameters
- `ResourceARNList`: A list of ARNs. An ARN (Amazon Resource Name) uniquely identifies a resource. For more information, see Amazon Resource Names (ARNs) and AWS Service Namespaces in the AWS General Reference.
- `Tags`: The tags that you want to add to the specified resources. A tag consists of a key and a value that you define.

"""

tag_resources(ResourceARNList, Tags; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("TagResources", Dict{String, Any}("ResourceARNList"=>ResourceARNList, "Tags"=>Tags); aws_config=aws_config)
tag_resources(ResourceARNList, Tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("TagResources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARNList"=>ResourceARNList, "Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResources()

Removes the specified tags from the specified resources. When you specify a tag key, the action removes both that key and its associated value. The operation succeeds even if you attempt to remove tags from a resource that were already removed. Note the following:   To remove tags from a resource, you need the necessary permissions for the service that the resource belongs to as well as permissions for removing tags. For more information, see this list.   You can only tag resources that are located in the specified Region for the AWS account.

# Required Parameters
- `ResourceARNList`: A list of ARNs. An ARN (Amazon Resource Name) uniquely identifies a resource. For more information, see Amazon Resource Names (ARNs) and AWS Service Namespaces in the AWS General Reference.
- `TagKeys`: A list of the tag keys that you want to remove from the specified resources.

"""

untag_resources(ResourceARNList, TagKeys; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("UntagResources", Dict{String, Any}("ResourceARNList"=>ResourceARNList, "TagKeys"=>TagKeys); aws_config=aws_config)
untag_resources(ResourceARNList, TagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = resource_groups_tagging_api("UntagResources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ResourceARNList"=>ResourceARNList, "TagKeys"=>TagKeys), args)); aws_config=aws_config)


text:
# # My Example
#
# Here we will try to create a literate example that is also code

# First we import our module
using TrialDocs, Plots

# Now we use the functions in the module
# Specifically the test_args_kw function

test_args_kw(1, "tomorrow", ones(2,2), keyword = "my keyword")

# ## Math
# We can test some math
# ```math
#   \frac{1}{2}
# ```
# ## Plotting
# Now we can test plotting
x = randn(3)
scatter(x, grid = true, gridstyle = :dash, gridalpha = 0.25, framestyle = :box, label = false, title = "simple graph", xlabel = "horizontal axis", ylabel = "vertical axis")

# !!! tip
#     maybe?

ext: jl

max_stars_repo_name: sandreza/TrialDocs

max_issues_repo_name: sandreza/TrialDocs

max_forks_repo_name: sandreza/TrialDocs


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ProductPartitioner(p₁, p₂)

A method for partitioning spatial objects using the product of two
partitioners `p₁` and `p₂`.
"""
struct ProductPartitioner{P1<:AbstractPartitioner,
                          P2<:AbstractPartitioner} <: AbstractPartitioner
  p₁::P1
  p₂::P2
end

# general case
function partition(object::AbstractSpatialObject,
                   partitioner::ProductPartitioner)
  # individual partition results
  s₁ = subsets(partition(object, partitioner.p₁))
  s₂ = subsets(partition(object, partitioner.p₂))

  # label-based representation
  l₁ = Vector{Int}(undef, npoints(object))
  l₂ = Vector{Int}(undef, npoints(object))
  for (i, inds) in enumerate(s₁)
    l₁[inds] .= i
  end
  for (i, inds) in enumerate(s₂)
    l₂[inds] .= i
  end

  # product of labels
  counter = 0
  d = Dict{Tuple{Int,Int},Int}()
  l = Vector{Int}(undef, npoints(object))
  for i in 1:npoints(object)
    t = (l₁[i], l₂[i])
    if t ∉ keys(d)
      counter += 1
      d[t] = counter
    end
    l[i] = d[t]
  end

  # return partition using label predicate
  pred(i, j) = l[i] == l[j]
  partition(object, PredicatePartitioner(pred))
end

# predicate partitioner
function partition(object::AbstractSpatialObject,
                   partitioner::ProductPartitioner{P1,P2}) where {P1<:AbstractPredicatePartitioner,
                                                                  P2<:AbstractPredicatePartitioner}
  pred(i, j) = partitioner.p₁(i, j) * partitioner.p₂(i, j)
  partition(object, PredicatePartitioner(pred))
end

# spatial predicate partitioner
function partition(object::AbstractSpatialObject,
                   partitioner::ProductPartitioner{P1,P2}) where {P1<:AbstractSpatialPredicatePartitioner,
                                                                  P2<:AbstractSpatialPredicatePartitioner}
  pred(x, y) = partitioner.p₁(x, y) * partitioner.p₂(x, y)
  partition(object, SpatialPredicatePartitioner(pred))
end

Base.:*(p₁::AbstractPartitioner, p₂::AbstractPartitioner) =
  ProductPartitioner(p₁, p₂)


text:
"""
    Newton{M}

`Newton` dispatches to the Newton direction for `LineSearch` algorithm and uses the Hessian
in `TrustRegion` models.

# Fields
- `modify::M`: the Hessian modification method to enforce positive definiteness
"""
struct Newton{M}
    modify::M
end

Newton() = Newton(MultipleIdentity())

order(::Newton) = SecondOrder()

function approx_hessian(n::Newton, state)
    ∇²f = state.∇²f
    isposdef(∇²f) && return ∇²f
    return n.modify(∇²f)
end

function direction(n::Newton, state)
    ∇f = state.∇f
    ∇²f = approx_hessian(n::Newton, state)
    return -inv(∇²f) * ∇f
end

"""
    MultipleIdentity{B, S}

Modify the Hessian by adding some multiples, τ, of the identity matrix, as described in
Nocedal and Wright's Numerical Optimization.

# Fields
- `β::B`: constant to compute appropriate τ
- `scale::S`: scale factor to expand τ
"""
struct MultipleIdentity{B, S}
    β::B
    scale::S

    function MultipleIdentity(β, ρ)
        β > 0 || error("β must be positive")
        ρ > 1 || error("ρ must be larger than 1")
        new{typeof(β), typeof(ρ)}(β, ρ)
    end
end

"""
    MultipleIdentity(; beta=1e-3, scale=2.0)

Initiate `MultipleIdentity` method to modify the Hessian.
"""
MultipleIdentity(; beta=1e-3, scale=2.0) = MultipleIdentity(beta, scale)

function (mi::MultipleIdentity)(∇²f::Matrix{T}) where T
    β = T(mi.β)
    scale = T(mi.scale)
    τ = -minimum(diag(∇²f)) + β
    h = similar(∇²f)
    while true
        h .= ∇²f + τ*I
        isposdef(h) && return h
        τ = max(scale*τ, β)
    end
end


text:
using Base.Cartesian
using StaticArrays

export Lanczos, Lanczos4OpenCV

abstract type AbstractLanczos <: InterpolationType end

"""
    Lanczos{N}(a=4)

Lanczos resampling via a kernel with scale parameter `a` and support over `N` neighbors.

This form of interpolation is merely the discrete convolution of the samples with a Lanczos kernel of size `a`. The size is directly related to how "far" the interpolation will reach for information, and has `O(N^2)` impact on runtime. An alternative implementation matching `lanczos4` from OpenCV is available as Lanczos4OpenCV.
"""
struct Lanczos{N} <: AbstractLanczos
    a::Int

    function Lanczos{N}(a) where N
        N < a && @warn "Using a smaller support than scale for Lanczos window. Proceed with caution."
        new{N}(a)
    end
end

Lanczos(a=4) = Lanczos{a}(a)

"""
    LanczosInterpolation
"""
struct LanczosInterpolation{T,N,IT <: DimSpec{AbstractLanczos},A <: AbstractArray{T,N},P <: Tuple{Vararg{AbstractArray,N}}} <: AbstractInterpolation{T,N,IT}
    coefs::A
    parentaxes::P
    it::IT
end

@generated degree(::Lanczos{N}) where {N} = :($N)

getknots(itp::LanczosInterpolation) = axes(itp)
coefficients(itp::LanczosInterpolation) = itp.coefs
itpflag(itp::LanczosInterpolation) = itp.it

size(itp::LanczosInterpolation) = map(length, itp.parentaxes)
axes(itp::LanczosInterpolation) = itp.parentaxes
lbounds(itp::LanczosInterpolation) = map(first, itp.parentaxes)
ubounds(itp::LanczosInterpolation) = map(last, itp.parentaxes)

function interpolate(A::AbstractArray{T}, it::AbstractLanczos) where T
    Apad = copy_with_padding(float(T), A, it)
    return LanczosInterpolation(Apad, axes(A), it)
end

@inline function (itp::LanczosInterpolation{T,N})(x::Vararg{<:Number,N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
    itp.coefs[wis...]
end

function weightedindex_parts(fs, it::AbstractLanczos, ax::AbstractUnitRange{<:Integer}, x)
    pos, δx = positions(it, ax, x)
    (position = pos, coefs = fmap(fs, it, δx))
end

function positions(it::AbstractLanczos, ax, x)
    xf = floorbounds(x, ax)
    δx = x - xf
    fast_trunc(Int, xf) - degree(it) + 1, δx
end

function value_weights(it::Lanczos, δx::S) where S
    N = degree(it)
    # short-circuit if integral
    isinteger(δx) && return ntuple(i->convert(float(S), i == N - δx), Val(2N))

    # LUTs
    #it.a === N === 4 && return _lanczos4(δx)

    cs = ntuple(i -> lanczos(N - i + δx, it.a, N), Val(2N))
    sum_cs = sum(cs)
    normed_cs = ntuple(i -> cs[i] / sum_cs, Val(length(cs)))
    return normed_cs
end

function padded_axis(ax::AbstractUnitRange, it::AbstractLanczos)
    N = degree(it)
    return first(ax) - N + 1:last(ax) + N
end

# precise implementations for fast evaluation of common kernels

"""
    lanczos(x, a, n=a)

Implementation of the [Lanczos kernel](https://en.wikipedia.org/wiki/Lanczos_resampling)
"""
lanczos(x::T, a::Integer, n=a) where {T} = abs(x) < n ? T(sinc(x) * sinc(x / a)) : zero(T)

"""
    Lanczos4OpenCV()

Alternative implementation of Lanczos resampling using algorithm `lanczos4` function of OpenCV:
https://github.com/opencv/opencv/blob/de15636724967faf62c2d1bce26f4335e4b359e5/modules/imgproc/src/resize.cpp#L917-L946
"""
struct Lanczos4OpenCV <: AbstractLanczos end

degree(::Lanczos4OpenCV) = 4

value_weights(::Lanczos4OpenCV, δx::S) where S = ifelse(isinteger(δx),ntuple(i->convert(float(S), i == 4 - δx), Val(8)) ,_lanczos4_opencv(δx))

# s45 = sqrt(2)/2
const s45 = 0.70710678118654752440084436210485
"""
    l4_2d_cs is a lookup table that could be generated by
    x = (0:7)*45*5
    l4_2d_cs = [cosd.(x) sind.(x)]'
"""
const l4_2d_cs = SA[1 0; -s45 -s45; 0 1; s45 -s45; -1 0; s45 s45; 0 -1; -s45 s45]

function _lanczos4_opencv(δx)
    p_4 = π / 4
    y0 = -(δx + 3) * p_4
    s0, c0 = sincos(y0)
    cs = ntuple(8) do i
        y = (δx + 4 - i) * p_4
        # Numerator is the sin subtraction identity
        # It is equivalent to the following
        # f(δx,i) = sin( π/4*( 5*(i-1)-δx-3 ) )
        (l4_2d_cs[i, 1] * s0 + l4_2d_cs[i, 2] * c0) / y^2
    end
    sum_cs = sum(cs)
    normed_cs = ntuple(i -> cs[i] / sum_cs, Val(8))
    return normed_cs
end


text:
"""
    struct CertificateMonomials <: MOI.AbstractConstraintAttribute end

A constraint attribute for the monomials indexing the
[`GramMatrixAttribute`](@ref) and [`MomentMatrixAttribute`](@ref) certificates.
"""
struct CertificateMonomials <: MOI.AbstractConstraintAttribute end

"""
    GramMatrixAttribute(N)
    GramMatrixAttribute()

A constraint attribute for the [`GramMatrix`](@ref) of a constraint, that is,
the positive semidefinite matrix `Q` indexed by the monomials in the vector `X`
such that ``X^\\top Q X`` is the sum-of-squares certificate of the constraint.
"""
struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute
    N::Int
end
GramMatrixAttribute() = GramMatrixAttribute(1)

"""
    MomentMatrixAttribute(N)
    MomentMatrixAttribute()

A constraint attribute fot the `MomentMatrix` of a constraint.
"""
struct MomentMatrixAttribute <: MOI.AbstractConstraintAttribute
    N::Int
end
MomentMatrixAttribute() = MomentMatrixAttribute(1)

"""
    LagrangianMultipliers(N)
    LagrangianMultipliers()

A constraint attribute fot the `LagrangianMultipliers` associated to the
inequalities of the domain of a constraint. There is one multiplier per
inequality and no multiplier for equalities as the equalities are handled by
reducing the polynomials over the ideal they generate instead of explicitely
creating multipliers.
"""
struct LagrangianMultipliers <: MOI.AbstractConstraintAttribute
    N::Int
end
LagrangianMultipliers() = LagrangianMultipliers(1)

# Needs to declare it set by optimize that it is not queried in the Caching
# optimize, even of `CertificateMonomials` which is set befor optimize.
function MOI.is_set_by_optimize(::Union{CertificateMonomials,
                                        GramMatrixAttribute,
                                        MomentMatrixAttribute,
                                        LagrangianMultipliers})
    return true
end

# If a variable is bridged, the `VectorOfVariables`-in-`SOSPolynomialSet` is
# bridged by `MOI.Bridges.Constraint.VectorFunctionizeBridge` and it has
# to pass the constraint to the SOS bridge.
function MOI.Bridges.Constraint.invariant_under_function_conversion(::Union{
    CertificateMonomials, GramMatrixAttribute,
    MomentMatrixAttribute, LagrangianMultipliers})
    return true
end

# This is type piracy but we tolerate it.
const ObjectWithoutIndex = GramMatrix{<:MOI.Utilities.ObjectWithoutIndex}
const ObjectOrTupleWithoutIndex = Union{ObjectWithoutIndex, Tuple{Vararg{ObjectWithoutIndex}}}
const ObjectOrTupleOrArrayWithoutIndex = Union{ObjectOrTupleWithoutIndex, AbstractArray{<:ObjectOrTupleWithoutIndex}}
MOI.Utilities.map_indices(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x
MOI.Utilities.substitute_variables(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x


text:
context("Primitives") do
    @fact Pyramid(Point(0,0,0),1,2) --> Pyramid{Int}(Point{3,Int}((0,0,0)),1,2)
    @fact Particle(Point(0,0,0), Vec(1,1,1)) --> Particle{3,Int}(Point{3,Int}((0,0,0)),Vec{3,Int}((1,1,1)))
    @fact HyperCube(Vec(0,0,0), 1) --> HyperCube{3,Int}(Vec{3,Int}(0,0,0), 1)

    centered_rect = centered(HyperCube)
    @fact centered_rect --> HyperCube{3,Float32}(Vec3f0(-0.5), 1f0)
    centered_rect = centered(HyperCube{2})
    @fact centered_rect --> HyperCube{2,Float32}(Vec2f0(-0.5f0), 1f0)

    centered_rect = centered(HyperCube{2, Float64})
    @fact centered_rect --> HyperCube{2,Float64}(Vec(-0.5,-0.5), 1.)

    centered_rect = centered(HyperCube{3, Float32})
    @fact centered_rect --> HyperCube{3,Float32}(Vec3f0(-0.5), 1f0)

end

context("Simplex") do
    s = Simplex(Point(1,2,3))
    @fact s --> Simplex{1,Point{3,Int}}((Point{3,Int}((1,2,3)),))
    s = Simplex(Point(1,2), Point(2,3))
    @fact s --> Simplex{2,Point{2,Int}}((Point{2,Int}((1,2)),Point{2,Int}((2,3))))
end


text:
using Documenter
using JuSwalbe
doctest(JuSwalbe; manual = false)

makedocs(
    sitename="JuSwalbe.jl",
    modules=[JuSwalbe],
    pages=["index.md", "manual.md", "reference.md", "devnotes.md"])

deploydocs(
    repo = "github.com/Zitzeronion/JuSwalbe.jl.git",
)

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
#=
The following examples are taken from:
W.H. Enright and H. Hayashi, The evaluation of numerical software for delay
differential equations.
=#

# Problem A1
@doc raw"""
    prob_dde_DDETST_A1

Delay differential equation model of blood production, given by

```math
u'(t) = \frac{0.2 u(t - 14)}{1 + u(t - 14)^{10}} - 0.1 u(t)
```

for ``t \in [0, 500]`` and history function ``\phi(t) = 0.5`` for ``t \leq 0``.

# References

Mackey, M. C. and Glass, L. (1977). Oscillation and chaos in physiological control systems,
Science (197), pp. 287-289.
"""
prob_dde_DDETST_A1

function f_dde_DDETST_A1(u, h, p, t)
  z = h(p, t - 14)

  0.2 * z / (1 + z^10) - 0.1 * u
end

function h_dde_DDETST_A1(p, t)
  t ≤ 0 || error("history function is only implemented for t ≤ 0")

  0.5
end

const prob_dde_DDETST_A1 =
  DDEProblem(f_dde_DDETST_A1, h_dde_DDETST_A1, (0.0, 500.0); constant_lags = [14])

# Problem A2
@doc raw"""
    prob_dde_DDETST_A2

Delay differential equation model of chronic granulocytic leukemia, given by

```math
u_1'(t) = \frac{1.1}{1 + \sqrt{10} u_1(t - 20)^{5/4}} - \frac{10 u_1(t)}{1 + 40 u_2(t)},
```
```math
u_2'(t) = \frac{100 u_1(t)}{1 + 40 u_2(t)} - 2.43 u_2(t),
```

for ``t \in [0, 100]`` and history function

```math
\phi_1(t) = 1.05767027/3,
```
```math
\phi_2(t) = 1.030713491/3,
```

for ``t \leq 0``.

# References

Wheldon, T., Kirk, J. and Finlay, H. (1974). Cyclical granulopoiesis in chronic granulocytic
leukemia: A simulation study., Blood (43), pp. 379-387.
"""
prob_dde_DDETST_A2

function f_dde_DDETST_A2!(du, u, h, p, t)
  z = 10 * u[1] / (1 + 40 * u[2])

  du[1] = 1.1 / (1 + sqrt(10) * (h(p, t - 20; idxs = 1)^(5/4))) - z
  du[2] = 10 * z - 2.43 * u[2]

  nothing
end

function h_dde_DDETST_A2(p, t; idxs::Union{Nothing,Int} = nothing)
  t ≤ 0 || error("history function is only implemented for t ≤ 0")

  if idxs === nothing
    [1.05767027/3, 1.030713491/3]
  elseif idxs == 1
    1.05767027/3
  elseif idxs == 2
    1.030713491/3
  else
    error("delay differential equation consists of two components")
  end
end

const prob_dde_DDETST_A2 =
  DDEProblem(f_dde_DDETST_A2!, h_dde_DDETST_A2, (0.0, 100.0); constant_lags = [20])

# Problem B2
@doc raw"""
    prob_dde_DDETST_B1

Delay differential equation

```math
u'(t) = 1 - u(\exp(1 - 1/t))
```

for ``t \in [0.1, 10]`` with history function ``\phi(t) = \log t`` for ``t \in (0, 0.1]``.

# Solution

The analytical solution for ``t \in [0.1, 10]`` is

```math
u(t) = \log t.
```

# References

Neves, K. W. (1975). Automatic integration of functional differential equations: An
approach, ACM Trans. Math. Soft. (1), pp. 357-368.
"""
prob_dde_DDETST_B1

function f_dde_DDETST_B1(u, h, p, t)
  1 - h(p, exp(1 - 1/t))
end

function h_dde_DDETST_B1(p, t)
  0 < t ≤ 0.1 || error("history function is only implemented for t ∈ (0, 0.1]")

  log(t)
end

function fanalytic_dde_DDETST_B1(u₀, ::typeof(h_dde_DDETST_B1), p, t)
  0.1 ≤ t ≤ 10 && u₀ == log(0.1) ||
    error("analytical solution is only implemented for t ∈ [0.1, 10] and u(0) = log 0.1")

  log(t)
end

const prob_dde_DDETST_B1 =
  DDEProblem(DDEFunction(f_dde_DDETST_B1; analytic = fanalytic_dde_DDETST_B1),
             h_dde_DDETST_B1, (0.1, 10.0);
             dependent_lags = ((u, p, t) -> t - exp(1 - 1/t),))

# Problem B2
@doc raw"""
    prob_dde_DDETST_B2

Delay differential equation

```math
u'(t) = - 1 - u(t) + 2 [u(t / 2) < 0]
```

for ``t \in [0, 2 \log 66]`` with history function ``\phi(0) = 1``.

# Solution

The analytical solution for ``t \in [0, 2 \log 66]`` is

```math
u(t) = \begin{cases}
  2 \exp(-t) - 1 & \text{if } t \in [0, 2 \log 2], \\
  1 - 6 \exp(-t) & \text{if  } t \in (2 \log 2, 2 \log 6], \\
  66 \exp(-t) - 1 & \text{if } t \in (2 \log 6, 2 \log 66].
\end{cases}
```

# References

Neves, K. W. and Thompson, S. (1992). Solution of systems of functional differential
equations with state dependent delays, Technical Report TR-92-009, Computer Science,
Radford University.
"""
prob_dde_DDETST_B2

function f_dde_DDETST_B2(u, h, p, t)
  if h(p, t/2) < 0
    1 - u
  else
    -1 - u
  end
end

function h_dde_DDETST_B2(p, t)
  t == 0 || error("history function is only implemented for t = 0")

  1.0
end

function fanalytic_dde_DDETST_B2(u₀, ::typeof(h_dde_DDETST_B2), p, t)
  0 ≤ t ≤ 2 * log(66) && u₀ == 1 ||
    error("analytical solution is only implemented for t ∈ [0, 2 log(66)] and u(0) = 1")

  if t ≤ 2 * log(2)
    2 * exp(-t) - 1
  elseif t ≤ 2 * log(6)
    1 - 6 * exp(-t)
  else
    66 * exp(-t) - 1
  end
end

const prob_dde_DDETST_B2 =
  DDEProblem(DDEFunction(f_dde_DDETST_B2; analytic = fanalytic_dde_DDETST_B2),
             h_dde_DDETST_B2, (0.0, 2 * log(66));
             dependent_lags = ((u, p, t) -> t / 2,))

# Problem C1
@doc raw"""
    prob_dde_DDETST_C1

Delay differential equation

```math
u'(t) = - 2 u(t - 1 - |u(t)|) (1 - u(t)^2)
```

for ``t \in [0, 30]`` with history function ``\phi(t) = 0.5`` for ``t \leq 0``.

# References

Paul, C. A. H. (1994). A test set of functional differential equations, Technical Report
249, The Department of Mathematics, The University of Manchester, Manchester, England.
"""
prob_dde_DDETST_C1

f_dde_DDETST_C1(u, h, p, t) = - 2 * h(p, t - 1 - abs(u)) * (1 - u^2)

function h_dde_DDETST_C1(p, t)
  t ≤ 0 || error("history function is only implemented for t ≤ 0")

  0.5
end

const prob_dde_DDETST_C1 =
  DDEProblem(f_dde_DDETST_C1, h_dde_DDETST_C1, (0.0, 30.0);
             dependent_lags = ((u, p, t) -> 1 + abs(u),))

# Problem C2
@doc raw"""
    prob_dde_DDETST_C2

Delay differential equation

```math
u_1'(t) = - 2 u_1(t - u_2(t)),
```
```math
u_₂'(t) = \frac{|u_1(t - u_2(t))| - |u_1(t)|}{1 + |u_1(t - u_2(t))|},
```

for ``t \in [0, 40]`` with history function

```math
\phi_1(t) = 1,
```
```math
\phi_2(t) = 0.5,
```

for ``t \leq 0``.

# References

Paul, C. A. H. (1994). A test set of functional differential equations, Technical Report
249, The Department of Mathematics, The University of Manchester, Manchester, England.
"""
prob_dde_DDETST_C2

function f_dde_DDETST_C2!(du, u, h, p, t)
  z = h(p, t - u[2]; idxs = 1)
  absz = abs(z)

  du[1] = -2 * z
  du[2] = (absz - abs(u[1])) / (1 + absz)

  nothing
end

function h_dde_DDETST_C2(p, t; idxs::Union{Nothing,Int} = nothing)
  t ≤ 0 || error("history function is only implemented for t ≤ 0")

  if idxs === nothing
    [1.0, 0.5]
  elseif idxs == 1
    1.0
  elseif idxs == 2
    0.5
  else
    error("delay differential equation consists of two components")
  end
end

const prob_dde_DDETST_C2 =
  DDEProblem(f_dde_DDETST_C2!, h_dde_DDETST_C2, (0.0, 30.0);
             dependent_lags = ((u, p, t) -> u[2],))

# Problem C3
@doc raw"""
    prob_dde_DDETST_C3

Delay differential equation model of hematopoiesis, given by

```math
u_1'(t) = \hat{s}_0 u_2(t - T_1) - \gamma u_1(t) - Q,
```
```math
u_2'(t) = f(u_1(t)) - k u_2(t),
```
```math
u_3'(t) = 1 - \frac{Q \exp(\gamma u_3(t))}{\hat{s}_0 u_2(t - T_1 - u_3(t))},
```

for ``t \in [0, 300]`` with history function ``\phi_1(0) = 3.325``, ``\phi_3(0) = 120``, and

```math
\phi_2(t) = \begin{cases}
  10 & \text{if } t \in [- T_1, 0],\\
  9.5 & \text{if } t < - T_1,
\end{cases}
```

where ``f(y) = a / (1 + K y^r)``, ``\hat{s}_0 = 0.0031``, ``T_1 = 6``, ``\gamma = 0.001``,
``Q = 0.0275``, ``k = 2.8``, ``a = 6570``, ``K = 0.0382``, and ``r = 6.96``.

# References

Mahaffy, J. M., Belair, J. and Mackey, M. C. (1996). Hematopoietic model with moving
boundary condition and state dependent delay, Private communication.
"""
const prob_dde_DDETST_C3 =
  let s₀ = 0.0031, T₁ = 6, γ = 0.001, Q = 0.0275, k = 2.8, a = 6570, K = 0.0382, r = 6.96
    global function f_dde_DDETST_C3!(du, u, h, p, t)
      du[1] = s₀ * h(p, t - T₁; idxs = 2) - γ * u[1] - Q
      du[2] = a / (1 + K * u[1]^r) - k * u[2]
      du[3] = 1 - Q * exp(γ * u[3]) / (s₀ * h(p, t - T₁ - u[3]; idxs = 2))

      nothing
    end

    global function h_dde_DDETST_C3(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
        [3.325, t < - T₁ ? 9.5 : 10.0, 120.0]
      elseif idxs == 1
        3.325
      elseif idxs == 2
        t < - T₁ ? 9.5 : 10.0
      elseif idxs == 3
        120.0
      else
        error("delay differential equation consists of three components")
      end
    end

    DDEProblem(f_dde_DDETST_C3!, h_dde_DDETST_C3, (0.0, 300.0);
               constant_lags = [T₁], dependent_lags = ((u, p, t) -> T₁ + u[3],))
  end

# Problem C4
@doc raw"""
    prob_dde_DDETST_C4

Delay differential equation model of hematopoiesis, given by the same delay differential
equation as [`prob_dde_DDETST_C3`](@ref)

```math
u_1'(t) = \hat{s}_0 u_2(t - T_1) - \gamma u_1(t) - Q,
```
```math
u_2'(t) = f(u_1(t)) - k u_2(t),
```
```math
u_3'(t) = 1 - \frac{Q \exp(\gamma u_3(t))}{\hat{s}_0 u_2(t - T_1 - u_3(t))},
```

for ``t \in [0, 100]`` with history function
``\phi_1(0) = 3.5``, ``\phi_3(0) = 50``, and ``\phi_2(t) = 10`` for ``t \leq 0``, where
``f(y) = a / (1 + K y^r)``, ``\hat{s}_0 = 0.00372``, ``T_1 = 3``, ``\gamma = 0.1``,
``Q = 0.00178``, ``k = 6.65``, ``a = 15600``, ``K = 0.0382``, and ``r = 6.96``.

# References

Mahaffy, J. M., Belair, J. and Mackey, M. C. (1996). Hematopoietic model with moving
boundary condition and state dependent delay, Private communication.
"""
const prob_dde_DDETST_C4 =
  let s₀ = 0.00372, T₁ = 3, γ = 0.01, Q = 0.00178, k = 6.65, a = 15600, K = 0.0382, r = 6.96
    global function f_dde_DDETST_C4!(du, u, h, p, t)
      du[1] = s₀ * h(p, t - T₁; idxs = 2) - γ * u[1] - Q
      du[2] = a / (1 + K * u[1]^r) - k * u[2]
      du[3] = 1 - Q * exp(γ * u[3]) / (s₀ * h(p, t - T₁ - u[3]; idxs = 2))

      nothing
    end

    global function h_dde_DDETST_C4(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
        [3.5, 10.0, 50.0]
      elseif idxs == 1
        3.5
      elseif idxs == 2
        10.0
      elseif idxs == 3
        50.0
      else
        error("delay differential equation consists of three components")
      end
    end

    DDEProblem(f_dde_DDETST_C4!, [3.5, 10.0, 50.0], h_dde_DDETST_C4, (0.0, 100.0);
               constant_lags = [T₁], dependent_lags = ((u, p, t) -> T₁ + u[3],))
  end

# Problem D1
@doc raw"""
    prob_dde_DDETST_D1

Delay differential equation

```math
u_1'(t) = u_2(t),
```
```math
u_2'(t) = - u_2(\exp(1 - u_2(t))) u_2(t)^2 \exp(1 - u_2(t)),
```

for ``t \in [0.1, 5]`` with history function

```math
\phi_1(t) = \log t,
```
```math
\phi_2(t) = 1 / t,
```

for ``t \in (0, 0.1]``.

# Solution

The analytical solution for ``t \in [0.1, 5]`` is

```math
u_1(t) = \log t,
```
```math
u_2(t) = 1 / t.
```

# References

Neves, K. W. (1975). Automatic integration of functional differential equations: An
approach, ACM Trans. Math. Soft. (1), pp. 357-368.
"""
prob_dde_DDETST_D1

function f_dde_DDETST_D1!(du, u, h, p, t)
  s = exp(1 - u[2])

  du[1] = u[2]
  du[2] = -h(p, s; idxs = 2) * u[2]^2 * s

  nothing
end

function h_dde_DDETST_D1(p, t; idxs::Union{Nothing,Int} = nothing)
  0 < t ≤ 0.1 || error("history function is only implemented for 0 < t ≤ 0.1")

  if idxs === nothing
    [log(t), 1 / t]
  elseif idxs == 1
    log(t)
  elseif idxs == 2
    1 / t
  else
    error("delay differential equation consists of two components")
  end
end

function fanalytic_dde_DDETST_D1(u₀, ::typeof(h_dde_DDETST_D1), p, t)
  0.1 ≤ t ≤ 5 && u₀ == [log(0.1), 10] ||
    error("analytical solution is only implemented for t ∈ [0.1, 5] and u(0.1) = [log 0.1, 10]")

  [log(t), 1 / t]
end

const prob_dde_DDETST_D1 =
  DDEProblem(DDEFunction(f_dde_DDETST_D1!; analytic = fanalytic_dde_DDETST_D1),
             h_dde_DDETST_D1, (0.1, 5.0);
             dependent_lags = ((u, p, t) -> t - exp(1 - u[2]),))

# Problem D2
@doc raw"""
    prob_dde_DDETST_D2

Delay differential equation model of antigen antibody dynamics with fading memory, given by

```math
u_1'(t) = - r_1 u_1(t) u_2(t) + r_2 u_3(t),
```
```math
u_2'(t) = - r_1 u_1(t) u_2(t) + \alpha r_1 u_1(t - u_4(t)) u_2(t - u_4(t)),
```
```math
u_3'(t) = r_1 u_1(t) u_2(t) - r_2 u_3(t),
```
```math
u_4'(t) = 1 + \frac{3 \delta - u_1(t) u_2(t) - u_3(t)}{u_1(t - u_4(t)) u_2(t - u_4(t)) + u_3(t - u_4(t))} \exp(\delta u_4(t)),
```

for ``t \in [0, 40]`` with history function

```math
\phi_1(t) = 5,
```
```math
\phi_2(t) = 0.1,
```
```math
\phi_3(t) = 0,
```
```math
\phi_4(t) = 0,
```

for ``t \leq 0``, where ``r_1 = 0.02``, ``r_2 = 0.005``, ``\alpha = 3``, and ``\delta = 0.01``.

# References

Gatica, J. and Waltman, P. (1982). A threshold model of antigen antibody dynamics with
fading memory, in Lakshmikantham (ed.), Nonlinear phenomena in mathematical science,
Academic Press, New York, pp. 425-439.
"""
const prob_dde_DDETST_D2 =
  let r₁ = 0.02, r₂ = 0.005, α = 3, δ = 0.01
    global function f_dde_DDETST_D2!(du, u, h, p, t)
      u₁u₂ = u[1] * u[2]
      r₁u₁u₂ = r₁ * u₁u₂
      r₂u₃ = r₂ * u[3]

      v = h(p, t - u[4])
      v₁v₂ = v[1] * v[2]

      du[1] = -r₁u₁u₂ + r₂u₃
      du[2] = -r₁u₁u₂ + α * r₁ * v₁v₂
      du[3] = -du[1]
      du[4] = 1 + (3 * δ - u₁u₂ - u[3]) * exp(δ * u[4]) / (v₁v₂ + v[3])

      nothing
    end

    global function h_dde_DDETST_D2(p, t)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      [5.0, 0.1, 0.0, 0.0]
    end

    DDEProblem(f_dde_DDETST_D2!, h_dde_DDETST_D2, (0.0, 40.0);
               dependent_lags = ((u, p, t) -> u[4],))
  end

# Problem E1
@doc raw"""
    prob_dde_DDETST_E1

Delay differential equation model of a food-limited population, given by

```math
u(t) = r u(t) (1 - u(t - 1) - c u'(t - 1))
```

for ``t \in [0, 40]`` with history function ``\phi(t) = 2 + t`` for ``t \leq 0``, where
``r = \pi / \sqrt{3} + 1/20`` and ``c = \sqrt{3} / (2 \pi) - 1 / 25``.

# References

Kuang, Y. and Feldstein, A. (1991). Boundedness of solutions of a nonlinear nonautonomous
neutral delay equation, J. Math. Anal. Appl. (156), pp. 293-304.
"""
const prob_dde_DDETST_E1 =
  let r = π / sqrt(3) + 1 / 20, c = sqrt(3) / (2 * π) - 1 / 25
    global f_dde_DDETST_E1(u, h, p, t) = r * u * (1 - h(p, t - 1) - c * h(p, t - 1, Val{1}))

    global function h_dde_DDETST_E1(p, t)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      2 + t
    end

    global function h_dde_DDETST_E1(p, t, ::Type{Val{1}})
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      1.0
    end

    DDEProblem(f_dde_DDETST_E1, h_dde_DDETST_E1, (0.0, 40.0);
               constant_lags = [1], neutral = true)
  end

# Problem E2
@doc raw"""
    prob_dde_DDETST_E2

Delay differential equation model of a logistic Gauss-type predator-prey system, given by

```math
u_1'(t) = u_1(t) (1 - u_1(t - \tau) - \rho u_1'(t - \tau)) - \frac{u_2(t) u_1(t)^2}{u_1(t)^2 + 1},
```
```math
u_2'(t) = u_2(t) \left(\frac{u_1(t)^2}{u_1(t)^2 + 1} - \alpha\right),
```

for ``t \in [0, 2]`` with history function

```math
\phi_1(t) = 0.33 - t / 10,
```
```math
\phi_2(t) = 2.22 + t / 10,
```

for ``t \leq 0``, where ``\alpha = 0.1``, ``\rho = 2.9``, and ``\tau = 0.42``.

# References

Kuang, Y. (1991). On neutral delay logistics Gauss-type predator-prey systems, Dyn. Stab.
Systems (6), pp. 173-189.
"""
const prob_dde_DDETST_E2 =
  let α = 0.1, ρ = 2.9, τ = 0.42
    global function f_dde_DDETST_E2!(du, u, h, p, t)
      v = u[1]^2
      z = u[2] * v / (v + 1)

      du[1] = u[1] * (1 - h(p, t - τ; idxs = 1) - ρ * h(p, t - 1, Val{1}; idxs = 1)) - z
      du[2] = z - α * u[2]

      nothing
    end

    global function h_dde_DDETST_E2(p, t; idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
        [0.33 - 0.1 * t, 2.22 + 0.1 * t]
      elseif idxs == 1
        0.33 - 0.1 * t
      elseif idxs == 2
        2.22 + 0.1 * t
      else
        error("delay differential equation consists of two components")
      end
    end

    global function h_dde_DDETST_E2(p, t, ::Type{Val{1}};
                                    idxs::Union{Nothing,Int} = nothing)
      t ≤ 0 || error("history function is only implemented for t ≤ 0")

      if idxs === nothing
        [-0.1, -0.1]
      elseif idxs == 1 || idxs == 2
        -0.1
      else
        error("delay differential equation consists of two components")
      end
    end

    DDEProblem(f_dde_DDETST_E2!, h_dde_DDETST_E2, (0.0, 2.0);
               constant_lags = [τ], neutral = true)
  end

# Problem F1
@doc raw"""
    prob_dde_DDETST_F1

Delay differential equation

```math
u'(t) = 2 \cos(2t) u(t / 2)^{2 \cos t} + \log(u'(t / 2)) - \log(2 \cos t) - \sin t
```

for ``t \in [0, 1]`` with history function ``\phi(0) = 1`` and ``\phi'(0) = 2``.

# Solution

The analytical solution for ``t \in [0, 1]`` is

```math
u(t) = \exp(\sin(2t)).
```

# References

Jackiewicz, Z. (1981). One step methods for the numerical solution of Volterra functional
differential equations of neutral type, Applicable Anal. (12), pp. 1-11.
"""
prob_dde_DDETST_F1

function f_dde_DDETST_F1(u, h, p, t)
  s, c = sincos(t)
  c2 = 2 * c

  to2 = t / 2
  h0 = h(p, to2)
  h1 = h(p, to2, Val{1})

  2 * cos(2 * t) * h0^c2 + log(h1) - log(c2) - s
end

function h_dde_DDETST_F1(p, t)
  iszero(t) || error("history function is only implemented for t = 0")

  1.0
end

function h_dde_DDETST_F1(p, t, ::Type{Val{1}})
  iszero(t) || error("history function is only implemented for t = 0")

  2.0
end

function fanalytic_dde_DDETST_F1(u₀, ::typeof(h_dde_DDETST_F1), p, t)
  0 ≤ t ≤ 1 && u₀ == 1 ||
    error("analytical solution is only implemented for t ∈ [0, 1] and u(0) = 1")

  exp(sin(2 * t))
end

const prob_dde_DDETST_F1 =
  DDEProblem(DDEFunction(f_dde_DDETST_F1; analytic = fanalytic_dde_DDETST_F1),
             h_dde_DDETST_F1, (0.0, 0.1);
             dependent_lags = ((u, p, t) -> t / 2,), neutral = true)

# Problem F2
@doc raw"""
    prob_dde_DDETST_F2

Delay differential equation

```math
u'(t) = u'(2t - 0.5)
```

for ``t \in [0.25, 0.499]`` with history function ``\phi(t) = \exp(-t^2)`` and
``\phi'(t) = -2t \exp(-t^2)`` for ``t \leq 0.25``.

# Solution

The analytical solution for ``t \in [0.25, 0.499]`` is

```math
u(t) = u_i(t) = \exp(-4^i t^2 + B_i t + C_i) / 2^i + K_i
```

if ``t \in [x_i, x_{i + 1}]``, where

```math
x_i = (1 - 2^{-i}) / 2,
```
```math
B_i = 2 (4^{i-1} + B_{i-1}),
```
```math
C_i = - 4^{i-2} - B_{i-1} / 2 + C_{i-1},
```
```math
K_i = - \exp(-4^i x_i^2 + B_i x_i + C_i) / 2^i + u_{i-1}(x_i),
```

and ``B_0 = C_0 = K_0 = 0``.

# References

Neves, K. W. and Thompson, S. (1992). Solution of systems of functional differential
equations with state dependent delays, Technical Report TR-92-009, Computer Science,
Radford University.
"""
prob_dde_DDETST_F2

f_dde_DDETST_F2(u, h, p, t) = h(p, 2 * t - 0.5, Val{1})

function h_dde_DDETST_F2(p, t)
  t ≤ 0.25 || error("history function is only implemented for t ≤ 0.25")

  exp(-t^2)
end

function h_dde_DDETST_F2(p, t, ::Type{Val{1}})
  t ≤ 0.25 || error("history function is only implemented for t ≤ 0.25")

  -2 * t * exp(-t^2)
end

function fanalytic_dde_DDETST_F2(u₀, ::typeof(h_dde_DDETST_F2), p, t)
  0.25 ≤ t ≤ 0.499 && u₀ == exp(-0.25^2) ||
    error("analytical solution is only implemented for t ∈ [0.25, 0.499] and u(0.25) = exp(-0.0625)")

  i = 0
  x = 0.0
  B = C = K = 0.0

  while true
    i += 1
    y = (1 - 2^(-i)) / 2

    if y ≥ t
      return exp(- 4^(i - 1) * t^2 + B * t + C) / 2^(i - 1) + K
    else
      K += exp(- 4^(i - 1) * x^2 + B * x + C) / 2^(i-1)
      C += - 4^(i - 2) - B / 2
      B = 2 * (4^(i - 1) + B)
      K -= exp(- 4^i * x^2 + B * x + C) / 2^i
      x = y
    end
  end
end

const prob_dde_DDETST_F2 =
  DDEProblem(DDEFunction(f_dde_DDETST_F2; analytic = fanalytic_dde_DDETST_F2),
             h_dde_DDETST_F2, (0.25, 0.499);
             dependent_lags = ((u, p, t) -> 1/2 - t,), neutral = true)

# Problem F3
@doc raw"""
    prob_dde_DDETST_F3

Delay differential equation

```math
u'(t) = \exp(-u(t)) + L_3 \left[\sin(u'(\alpha(t))) - \sin\left(\frac{1}{3 + \alpha(t)}\right)\right]
```

for ``t \in [0, 10]`` with history function ``\phi(0) = \log 3`` and ``\phi'(0) = 1 / 3``,
where ``\alpha(t) = 0.5 t (1 - \cos(2 \pi t))`` and ``L_3 = 0.2``.

# Solution

The analytical solution for ``t \in [0, 10]`` is

```math
u(t) = \log(t + 3).
```
"""
prob_dde_DDETST_F3

let L₃ = 0.2
  global function f_dde_DDETST_F3(u, h, p, t)
    α = 0.5 * t * (1 - cos(2 * π * t))

    exp(-u) + L₃ * (sin(h(p, α, Val{1})) - sin(inv(3 + α)))
  end
end

function h_dde_DDETST_F345(p, t)
  iszero(t) || error("history function is only implemented for t = 0")

  log(3)
end

function h_dde_DDETST_F345(p, t, ::Type{Val{1}})
  iszero(t) || error("history function is only implemented for t = 0")

  inv(3)
end

function fanalytic_dde_DDETST_F345(u₀, ::typeof(h_dde_DDETST_F345), p, t)
  0 ≤ t ≤ 10 && u₀ == log(3) ||
    error("analytical solution is only implemented for t ∈ [0, 10] and u(0) = log 3")

  log(t + 3)
end

const prob_dde_DDETST_F3 =
  DDEProblem(DDEFunction(f_dde_DDETST_F3; analytic = fanalytic_dde_DDETST_F345),
             h_dde_DDETST_F345, (0.0, 10.0);
             dependent_lags = ((u, p, t) -> 0.5 * t * (1 + cos(2 * π * t)),),
             neutral = true)

# Problem F4
"""
    prob_dde_DDETST_F4

Same delay differential equation as [`prob_dde_DDETST_F3`](@ref) with ``L_3 = 0.4``.
"""
prob_dde_DDETST_F4

let L₃ = 0.4
  global function f_dde_DDETST_F4(u, h, p, t)
    α = 0.5 * t * (1 - cos(2 * π * t))

    exp(-u) + L₃ * (sin(h(p, α, Val{1})) - sin(inv(3 + α)))
  end
end

const prob_dde_DDETST_F4 =
  remake(prob_dde_DDETST_F3;
         f = DDEFunction(f_dde_DDETST_F4; analytic = fanalytic_dde_DDETST_F345))

# Problem F5
"""
    prob_dde_DDETST_F5

Same delay differential equation as [`prob_dde_DDETST_F3`](@ref) with ``L_3 = 0.6``.
"""
prob_dde_DDETST_F5

let L₃ = 0.6
  global function f_dde_DDETST_F5(u, h, p, t)
    α = 0.5 * t * (1 - cos(2 * π * t))

    exp(-u) + L₃ * (sin(h(p, α, Val{1})) - sin(inv(3 + α)))
  end
end

const prob_dde_DDETST_F5 =
  remake(prob_dde_DDETST_F3;
         f = DDEFunction(f_dde_DDETST_F5; analytic = fanalytic_dde_DDETST_F345))

# Problem G1
@doc raw"""
    prob_dde_DDETST_G1

Delay differential equation

```math
u'(t) = - u'(t - u(t)^2 / 4)
```

for ``t \in [0, 1]`` with history function ``\phi(t) = 1 - t`` for ``t \leq 0`` and
``\phi'(t) = -1`` for ``t < 0``.

# Solution

The analytical solution for ``t \in [0, 1]`` is

```math
u(t) = t + 1.
```

# References

El'sgol'ts, L. E. and Norkin, S. B. (1973). Introduction to the Theory and Application of
Differential Equations with Deviating Arguments, Academic Press, New York, p. 44.
"""
prob_dde_DDETST_G1

f_dde_DDETST_G1(u, h, p, t) = - h(p, t - u^2 / 4, Val{1})

function h_dde_DDETST_G1(p, t)
  t ≤ zero(t) || error("history function is only implemented for t ≤ 0")

  1 - t
end

function h_dde_DDETST_G1(p, t, ::Type{Val{1}})
  t < zero(t) || error("history function is only implemented for t < 0")

  - 1.0
end

function fanalytic_dde_DDETST_G1(u₀, ::typeof(h_dde_DDETST_G1), p, t)
  0 ≤ t ≤ 1 && u₀ == 1 ||
    error("analytical solution is only implemented for t ∈ [0, 1] and u(0) = 1")

  t + 1
end

const prob_dde_DDETST_G1 =
  DDEProblem(DDEFunction(f_dde_DDETST_G1; analytic = fanalytic_dde_DDETST_G1),
             h_dde_DDETST_G1, (0.0, 1.0);
             dependent_lags = ((u, p, t) -> u^2 / 4,), neutral = true)

# Problem G2
@doc raw"""
    prob_dde_DDETST_G2

Delay differential equation

```math
u'(t) = - u'(u(t) - 2)
```

for ``t \in [0, 1]`` with history function ``\phi(t) = 1 - t`` for ``t \leq 0`` and
``\phi'(t) = -1`` for ``t < 0``.

# Solution

The analytical solution for ``t \in [0, 1]`` is

```math
u(t) = t + 1.
```

El'sgol'ts, L. E. and Norkin, S. B. (1973). Introduction to the Theory and Application of
Differential Equations with Deviating Arguments, Academic Press, New York, pp. 44-45.
"""
prob_dde_DDETST_G2

f_dde_DDETST_G2(u, h, p, t) = - h(p, u - 2, Val{1})

function h_dde_DDETST_G2(p, t)
  t ≤ zero(t) || error("history function is only implemented for t ≤ 0")

  1 - t
end

function h_dde_DDETST_G2(p, t, ::Type{Val{1}})
  t < zero(t) || error("history function is only implemented for t < 0")

  - 1.0
end

function fanalytic_dde_DDETST_G2(u₀, ::typeof(h_dde_DDETST_G2), p, t)
  0 ≤ t ≤ 1 && u₀ == 1.0 ||
    error("analytical solution is only implemented for t ∈ [0, 1] and u(0) = 1")

  t + 1
end

const prob_dde_DDETST_G2 =
  DDEProblem(DDEFunction(f_dde_DDETST_G2; analytic = fanalytic_dde_DDETST_G2),
             h_dde_DDETST_G2, (0.0, 1.0);
             dependent_lags = ((u, p, t) -> t + 2 -  u,), neutral = true)

# Problem H1
@doc raw"""
    prob_dde_DDETST_H1

Delay differential equation

```math
u'(t) = - \frac{4 t u(t)^2}{4 + \log(\cos(2t))^2} + \tan(2t) + 0.5 \arctan\left(u'\left(\frac{t u(t)^2}{1 + u(t)^2}\right)\right)
```

for ``t \in [0, 0.225 \pi]`` with history function ``\phi(0) = 0`` and ``\phi'(0) = 0``.

# Solution

The analytical solution for ``t \in [0, 0.225 \pi]`` is

```math
u(t) = - \log(\cos(2t)) / 2.
```

# References

Castleton, R. N. and Grimm, L. J. (1973). A first order method for differential equations of
neutral type, Math. Comput. (27), pp. 571-577.
"""
prob_dde_DDETST_H1

function f_dde_DDETST_H1(u, h, p, t)
  v = u^2
  z = t * v

  - 4 * z / (4 + log(cos(2 * t))^2) + tan(2 * t) + 0.5 * atan(h(p, z / (1 + v), Val{1}))
end

function h_dde_DDETST_H1(p, t)
  iszero(t) || error("history function is only implemented for t = 0")

  0.0
end

function h_dde_DDETST_H1(p, t, ::Type{Val{1}})
  iszero(t) || error("history function is only implemented for t = 0")

  0.0
end

function fanalytic_dde_DDETST_H1(u₀, ::typeof(h_dde_DDETST_H1), p, t)
  0 ≤ t ≤ 0.225 * π && iszero(u₀) ||
    error("analytical solution is only implemented for t ∈ [0, 0.225π] and u(0) = 0")

  - log(cos(2 * t)) / 2
end

const prob_dde_DDETST_H1 =
  DDEProblem(DDEFunction(f_dde_DDETST_H1; analytic = fanalytic_dde_DDETST_H1),
             h_dde_DDETST_H1, (0.0, 0.225 * π);
             dependent_lags = ((u, p, t) -> t / (1 + u^2),), neutral = true)

# Problem H2
@doc raw"""
    prob_dde_DDETST_H2

Delay differential equation

```math
u'(t) = \cos(t) (1 + u(t u(t)^2)) + L_3 u(t) u'(t u(t)^2) + (1 - L_3) \sin(t) \cos(t \sin(t)^2) - \sin(t + t \sin(t)^2)
```

for ``t \in [0, \pi]`` with history function ``\phi(0) = 0`` and ``\phi'(0) = 1``, where
``L_3 = 0.1``.

# Solution

The analytical solution for ``t \in [0, \pi]`` is

```math
u(t) = \sin(t).
```

# References

Hayashi, H. (1996). Numerical solution of retarded and neutral delay differential equations
using continuous Runge-Kutta methods, PhD thesis, Department of Computer Science, University
of Toronto, Toronto, Canada.
"""
prob_dde_DDETST_H2

let L₃ = 0.1
  global function f_dde_DDETST_H2(u, h, p, t)
    s, c = sincos(t)
    a = t * s^2

    z = t * u^2

    c * (1 + h(p, z)) + L₃ * u * h(p, z, Val{1}) + (1 - L₃) * s * cos(a) - sin(t + a)
  end
end

function h_dde_DDETST_H234(p, t)
  iszero(t) || error("history function is only implemented for t = 0")

  0.0
end

function h_dde_DDETST_H234(p, t, ::Type{Val{1}})
  iszero(t) || error("history function is only implemented for t = 0")

  1.0
end

function fanalytic_dde_DDETST_H234(u₀, ::typeof(h_dde_DDETST_H234), p, t)
  0 ≤ t ≤ π && iszero(u₀) ||
    error("analytical solution is only implemented for t ∈ [0, π] and u(0) = 0")

  sin(t)
end

const prob_dde_DDETST_H2 =
  DDEProblem(DDEFunction(f_dde_DDETST_H2; analytic = fanalytic_dde_DDETST_H234),
             h_dde_DDETST_H234, (0.0, π);
             dependent_lags = ((u, p, t) -> t * (1 - u^2),), neutral = true)

# Problem H3
"""
    prob_dde_DDETST_H3

Same delay differential equation as [`prob_dde_DDETST_H2`](@ref) with ``L_3 = 0.3``.

# References

Hayashi, H. (1996). Numerical solution of retarded and neutral delay differential equations
using continuous Runge-Kutta methods, PhD thesis, Department of Computer Science, University
of Toronto, Toronto, Canada.
"""
prob_dde_DDETST_H3

let L₃ = 0.3
  global function f_dde_DDETST_H3(u, h, p, t)
    s, c = sincos(t)
    a = t * s^2

    z = t * u^2

    c * (1 + h(p, z)) + L₃ * u * h(p, z, Val{1}) + (1 - L₃) * s * cos(a) - sin(t + a)
  end
end

const prob_dde_DDETST_H3 =
  remake(prob_dde_DDETST_H2;
         f = DDEFunction(f_dde_DDETST_H3; analytic = fanalytic_dde_DDETST_H234))

# Problem H4
"""
    prob_dde_DDETST_H4

Same delay differential equation as [`prob_dde_DDETST_H2`](@ref) with ``L_3 = 0.5``.

# References

Hayashi, H. (1996). Numerical solution of retarded and neutral delay differential equations
using continuous Runge-Kutta methods, PhD thesis, Department of Computer Science, University
of Toronto, Toronto, Canada.
"""
prob_dde_DDETST_H4

let L₃ = 0.5
  global function f_dde_DDETST_H4(u, h, p, t)
    s, c = sincos(t)
    a = t * s^2

    z = t * u^2

    c * (1 + h(p, z)) + L₃ * u * h(p, z, Val{1}) + (1 - L₃) * s * cos(a) - sin(t + a)
  end
end

const prob_dde_DDETST_H4 =
  remake(prob_dde_DDETST_H2;
         f = DDEFunction(f_dde_DDETST_H4; analytic = fanalytic_dde_DDETST_H234))


text:
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:CRR},
    ::Val{:Pendulum},
    type::AbstractString;
    save_dir = nothing,
    seed = 123,
)
    rng = StableRNG(seed)
    inner_env = PendulumEnv(T = Float32, rng = rng)
    A = action_space(inner_env)
    low = A.left
    high = A.right
    ns = length(state(inner_env))
    na = 1

    trajectory_num = 10000
    dataset_size = 10000
    batch_size = 128

    env = ActionTransformedEnv(
        inner_env;
        action_mapping = x -> low + (x[1] + 1) * 0.5 * (high - low),
    )
    init = glorot_uniform(rng)

    create_policy_net() = GaussianNetwork(
        pre = Chain(
            Dense(ns, 64, relu),
            Dense(64, 64, relu),
        ),
        μ = Chain(Dense(64, na, init = init)),
        logσ = Chain(Dense(64, na, init = init)),
    )

    create_q_net() = Chain(
        Dense(ns + na, 64, relu; init = init),
        Dense(64, 64, relu; init = init),
        Dense(64, 1; init = init),
    )

    agent = Agent(
        policy = OfflinePolicy(
            learner = CRRLearner(
                approximator = ActorCritic(
                    actor = create_policy_net() |> cpu,
                    critic = create_q_net() |> cpu,
                    optimizer = ADAM(3e-3),
                ),
                target_approximator = ActorCritic(
                    actor = create_policy_net() |> cpu,
                    critic = create_q_net() |> cpu,
                    optimizer = ADAM(3e-3),
                ),
                γ = 0.99f0,
                batch_size = batch_size,
                policy_improvement_mode = :exp,
                ratio_upper_bound = 20.0f0,
                β = 1.0f0,
                advantage_estimator = :mean,
                m = 4,
                update_freq = 1,
                continuous = true,
            ),
            dataset = gen_JuliaRL_dataset(:SAC, :Pendulum, type; dataset_size = dataset_size),
            continuous = true,
            batch_size = batch_size,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
            action = Vector{Float32} => (na,),
        ),
    )

    stop_condition = StopAfterStep(trajectory_num, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "CRR <-> Pendulum ($type dataset)")
end

using Plots
pyplot() #hide
ex = E`JuliaRL_CRR_Pendulum(medium)`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_CRR_Pendulum_medium.png") #hide

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module SSLTestSubsetting
using SimpleSDMLayers
using Test

temp = SimpleSDMPredictor(WorldClim, BioClim, 1)
@test size(temp) == (1080, 2160)

coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
l1 = clip(temp; coords...)
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
tempfile = tempname()
geotiff(tempfile, l2)
l3 = geotiff(SimpleSDMPredictor, tempfile)

@test size(l1) == size(l2)
@test size(l1) == size(l3)
@test l1.grid == l2.grid
@test l1.grid == l3.grid

for l in (l1, l2, l3)
    @test size(l) == (330, 570)

    @test l.left == coords.left
    @test l.right == coords.right
    @test l.bottom == coords.bottom
    @test l.top == coords.top

    @test stride(l) == stride(temp)
    @test longitudes(l)[1] == -144.91666666666666
    @test longitudes(l)[end] == -50.083333333333336
    @test latitudes(l)[1] == 20.083333333333332
    @test latitudes(l)[end] == 74.91666666666667
end

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Date/DateTime Ranges

# Override default step; otherwise it would be Millisecond(1)
Base.colon(start::T, stop::T) where {T<:DateTime} = StepRange(start, Day(1), stop)
Base.colon(start::T, stop::T) where {T<:Date}     = StepRange(start, Day(1), stop)
Base.colon(start::T, stop::T) where {T<:Time}     = StepRange(start, Second(1), stop)

Base.range(start::DateTime, len::Integer)  = range(start, Day(1), len)
Base.range(start::Date, len::Integer)      = range(start, Day(1), len)

(::Type{StepRange{<:Dates.DatePeriod,<:Real}})(start, step, stop) =
    throw(ArgumentError("must specify step as a Period when constructing Dates ranges"))

# Given a start and end date, how many steps/periods are in between
guess(a::DateTime, b::DateTime, c) = floor(Int64, (Int128(value(b)) - Int128(value(a))) / toms(c))
guess(a::Date, b::Date, c) = Int64(div(value(b - a), days(c)))
len(a::Time, b::Time, c) = Int64(div(value(b - a), tons(c)))
function len(a, b, c)
    lo, hi, st = min(a, b), max(a, b), abs(c)
    i = guess(a, b, c) - 1
    while lo + st * i <= hi
        i += 1
    end
    return i - 1
end
Base.length(r::StepRange{<:TimeType}) = isempty(r) ? Int64(0) : len(r.start, r.stop, r.step) + 1
# Period ranges hook into Int64 overflow detection
Base.length(r::StepRange{<:Period}) = length(StepRange(value(r.start), value(r.step), value(r.stop)))

# Used to calculate the last valid date in the range given the start, stop, and step
# last = stop - steprem(start, stop, step)
Base.steprem(a::T, b::T, c) where {T<:TimeType} = b - (a + c * len(a, b, c))

import Base.in
function in(x::T, r::StepRange{T}) where T<:TimeType
    n = len(first(r), x, step(r)) + 1
    n >= 1 && n <= length(r) && r[n] == x
end

Base.start(r::StepRange{<:TimeType}) = 0
Base.next(r::StepRange{<:TimeType}, i::Int) = (r.start + r.step*i, i + 1)
Base.done(r::StepRange{<:TimeType,<:Period}, i::Integer) = length(r) <= i

+(x::Period, r::Range{<:TimeType}) = (x + first(r)):step(r):(x + last(r))
+(r::Range{<:TimeType}, x::Period) = x + r
-(r::Range{<:TimeType}, x::Period) = (first(r)-x):step(r):(last(r)-x)

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
# poisson processes
#
struct InhomogeneousPoisson end
struct InhomogeneousPoissonPiecewiseDecreasing end

"""
next_event_time(rate::Function, max_rate::Float64, rng::AbstractRNG)::Float64

Generate a new event from an inhomogeneous poisson process with rate Lambda(t).
Based on (Ogata’s Modified Thinning Algorithm: Ogata,  1981,  p.25,  Algorithm  2)
see also https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf

# Arguments
- `rate`: rate(dt) has to be defined outside (e..g t-> rate(t+t0,args))
- `max_rate`: maximal rate in near future (has to be evaluated externally)
- `rng`: random number generator

API - output
* returns the next event time

"""
function next_time(alg::InhomogeneousPoisson, rng::AbstractRNG, rate::Function, max_rate::Float64)::Float64
    dt = 0.0
    theta = 1.0 / max_rate
    while true
        # generate next event from bounding homogeneous Poisson process with max_rate
        dt += randexp(rng) * theta
        # accept next event with probability rate(t)/rate_max [Thinning algorithm]
        if rand(rng) < rate(dt) / max_rate
            return dt
        end
    end
end

"""
Generate a new event from an inhomogeneous poisson process with rate Lambda(t) under
the assumption that rate(dt) is monotonically decreasing.
Based on (Ogata’s Modified Thinning Algorithm: Ogata,  1981,  p.25,  Algorithm  3)
see also https://www.math.fsu.edu/~ychen/research/Thinning%20algorithm.pdf

# Arguments
- `rate`: rate(dt) has to be defined outside (e..g t-> rate(t+t0,args))
- `rng`: random number generator

API - output
* returns the next event time

"""
function next_time(alg::InhomogeneousPoissonPiecewiseDecreasing, rng::AbstractRNG, rate::Function)::Float64
    dt = 0.0
    while true
        # future rate can only be smaller than current rate
        max_rate = rate(dt)
        # generate next event from bounding homogeneous Poisson process with max_rate
        dt += randexp(rng) / max_rate
        # accept next event with probability rate(t)/rate_max [Thinning algorithm]
        if rand(rng) < rate(dt) / max_rate
            return dt
        end
    end
end

"""
Generate a new event id from a collection of inhomogeneous poisson processes with
rates Lambda(t).
# Arguments
- `rates`: rates(dt); Float -> [Float]
- `max_rate`: maximal rate in near future (has to be evaluated externally)
- `rng`: random number generator

API - output
* returns the next event id
"""
function next_event(alg::InhomogeneousPoisson, rng::AbstractRNG, rates::Array{Float64}, max_rate::Float64)::Int
    next_index = 1
    cumulated_rates = cumsum(rates)
    sum_rate = cumulated_rates[end]

    theta = rand(rng) * sum_rate
    id = 1
    # catch lower-bound case that cannot be reached by binary search
    if theta >= cumulated_rates[1]
        # binary search
        index_l = 1
        index_r = length(cumulated_rates)
        while index_l < index_r - 1
            # index_m = floor(Int,(index_l+index_r)/2)
            index_m = fld(index_l + index_r, 2)
            if cumulated_rates[index_m] < theta
                index_l = index_m
            else
                index_r = index_m
            end
        end
        id = index_r
    end

    return id
end

"""
Generate a new event from a collection of inhomogeneous poisson processes with
rates Lambda(t).
# Arguments
- `rates`: rates(dt); Float -> [Float]
- `max_rate`: maximal rate in near future (has to be evaluated externally)
- `rng`: random number generator

API - output
* returns the next event time and event id as tuple (dt, id)
"""
function next(alg::InhomogeneousPoisson, rng::AbstractRNG, rates::Function, max_rate::Float64)::Tuple{Float64,Int}
    rate(t) = sum(rates(t))
    dt = next_time(alg, rng, rate, max_rate)
    id = next_event(alg, rng, rates(dt), max_rate)
    return dt, id
end


text:
# global axes
const globalX = [1., 0., 0.]
const globalY = [0., 1., 0.]
const globalZ = [0., 0., 1.]

# expands elemental start/end nodes to global DOF indices
function dofExpander(element::Element, nodes::Vector{Node})
    element.dofIndex = vcat(nodes[element.nodeIndex[1]].globalIndex, nodes[element.nodeIndex[2]].globalIndex)
end

# dofExpander for a vector of elements
function dofExpander(elements::Vector{Element}, nodes::Vector{Node})
    for element in elements
        element.dofIndex = vcat(nodes[element.nodeIndex[1]].globalIndex, nodes[element.nodeIndex[2]].globalIndex)
    end
end

# global dof index of each node
function nodeGlobalIndex(nodes::Vector{Node})
    n = length(nodes)
    nNodalDOFS = length(nodes[1].DOFS)
    for i = 1:n
        nodes[i].globalIndex = i * nNodalDOFS - (nNodalDOFS-1) .+ collect(0:nNodalDOFS-1)
    end
end

```
Local coordinate system of element
```
function lcs(element::Element, Ψ; tol = 0.001)

    # local x vector
    xvec = normalize(element.posEnd .- element.posStart)

    if norm(cross(xvec, globalY)) < tol
        CYx = xvec[2] #cosine to global Y axis
        xvec = CYx * globalY
        yvec = -CYx * globalX * cos(Ψ) + sin(Ψ) * globalZ
        zvec = CYx * globalX * sin(Ψ) + cos(Ψ) * globalZ
    else
        zbar = normalize(cross(xvec, [0, 1, 0]))
        ybar = normalize(cross(zbar, xvec))

        yvec = cos(Ψ) * ybar + sin(Ψ) * zbar
        zvec = -sin(Ψ) * ybar + cos(Ψ) * zbar
    end

    return Vec3.([xvec, yvec, zvec])
end

# creates the elemental stiffness matrix in global coordinate system
function k_elemental(element::Element, dims::Int; return_k = false, tol = 1e-3)
    type = element.type
    if type == :truss
        if dims == 3
            Cx, Cy, Cz = (element.posEnd .- element.posStart) ./ element.length

            # rotation matrix
            R = [Cx Cy Cz 0 0 0; 0 0 0 Cx Cy Cz]
            element.R = R

            # stiffness matrix
            K = element.E * element.A / element.length .* [1 -1; -1 1]

            element.k = Symmetric(R' * K * R)
        elseif dims == 2
            l, m = (element.posEnd .- element.posStart) ./ element.length

            # rotation matrix
            R = [l m 0 0; 0 0 l m]
            element.R = R

            # stiffness matrix
            K = element.E * element.A / element.length .* [1 -1; -1 1]

            element.k = Symmetric(R' * K * R)
        else
            error("dims must be 2 or 3")
        end
    elseif type == :frame
        if dims == 3

            ##
            # deprecated
            ##
            #rotation matrix
            # CXx, CYx, CZx = (element.posEnd .- element.posStart) ./ element.length
            # D = sqrt(CXx^2 + CYx^2)
            # CXy = -CYx / D
            # CYy = CXx / D
            # CZy = 0
            # CXz = -CXx * CZx / D
            # CYz = -CYx * CZx / D
            # CZz = D

            # if abs(element.posStart[1]-element.posEnd[1]) < tol && abs(element.posStart[2]-element.posEnd[2]) < tol
            #     if element.posEnd[3] > element.posStart[3]
            #         Λ = [0 0 1; 0 1 0; -1 0 0]
            #     else
            #         Λ = [0 0 -1; 0 1 0; 1 0 0]
            #     end
            # else
            #     Λ = [CXx CYx CZx; CXy CYy CZy; CXz CYz CZz]
            # end

            ##
            # New method provides precise pitch angle control w/r/t local x axis
            xvec = normalize(element.posEnd .- element.posStart) # local x vector
            CXx, CYx, CZx = xvec # local x cosines

            if norm(cross(xvec, globalY)) < tol #special case for horizontal members aligned with global Y
                Λ = [0. CYx 0.;
                    -CYx*cos(element.Ψ) 0 sin(element.Ψ);
                    CYx*sin(element.Ψ) 0 cos(element.Ψ)]
            else # all other
                b1 = (-CXx * CYx * cos(element.Ψ) - CZx * sin(element.Ψ)) / sqrt(CXx^2 + CZx^2)
                b2 = sqrt(CXx^2 + CZx^2) * cos(element.Ψ)
                b3 = (-CYx * CZx * cos(element.Ψ) + CXx * sin(element.Ψ)) / sqrt(CXx^2 + CZx^2)

                c1 = (CXx * CYx * sin(element.Ψ) - CZx * cos(element.Ψ)) / sqrt(CXx^2 + CZx^2)
                c2 = -sqrt(CXx^2 + CZx^2) * sin(element.Ψ)
                c3 = (CYx * CZx * sin(element.Ψ) + CXx * cos(element.Ψ)) / sqrt(CXx^2 + CZx^2)

                Λ = [CXx CYx CZx;
                    b1 b2 b3;
                    c1 c2 c3]
            end

            R = [Λ zeros(3,9); zeros(3,3) Λ zeros(3,6); zeros(3,6) Λ zeros(3,3); zeros(3,9) Λ]
            element.R = R

            # Faulty method from Ferreira "MATLAB Codes for Finite Element Analysis"
            #stiffness matrix
            #nodal stiffnesses
            # k1 = element.E * element.A / element.length
            # k2 = 12 * element.E * element.Iz / element.length^3
            # k3 = 6 * element.E * element.Iz / element.length^2
            # k4 = 4 * element.E * element.Iz / element.length
            # k5 = 2 * element.E * element.Iz / element.length
            # k6 = 12 * element.E * element.Iy / element.length^3
            # k7 = 6 * element.E * element.Iy / element.length^2
            # k8 = 4 * element.E * element.Iy / element.length
            # k9 = 2 * element.E * element.Iy / element.length
            # k10 = element.G * element.J / element.length

            # #stiffness blocks
            # a = [k1 0 0; 0 k2 0; 0 0 k6]
            # b = [0 0 0; 0 0 k3; 0 -k7 0]
            # c = [k10 0 0; 0 k8 0; 0 0 k4]
            # d = [-k10 0 0; 0 k9 0; 0 0 k5]

            # element.k = Symmetric(R' * [a b -a b; b' c b d; -a' b' a -b; b' d' -b' c] * R)

            E = element.E
            length = element.length
            A = element.A
            Iz = element.Iz
            Iy = element.Iy
            G = element.G
            J = element.J
            k = E / length^3 * [
                A*length^2 0 0 0 0 0 -A*length^2 0 0 0 0 0;
                0 12Iz 0 0 0 6length*Iz 0 -12Iz 0 0 0 6length*Iz;
                0 0 12Iy 0 -6length*Iy 0 0 0 -12Iy 0 -6length*Iy 0;
                0 0 0 G*J*length^2/E 0 0 0 0 0 -G*J*length^2/E 0 0;
                0 0 -6length*Iy 0 4length^2*Iy 0 0 0 6length*Iy 0 2length^2*Iy 0;
                0 6length*Iz 0 0 0 4length^2*Iz 0 -6length*Iz 0 0 0 2length^2*Iz;
                -A*length^2 0 0 0 0 0 A*length^2 0 0 0 0 0;
                0 -12Iz 0 0 0 -6length*Iz 0 12Iz 0 0 0 -6length*Iz;
                0 0 -12Iy 0 6length*Iy 0 0 0 12Iy 0 6length*Iy 0;
                0 0 0 -G*J*length^2/E 0 0 0 0 0 G*J*length^2/E 0 0;
                0 0 -6length*Iy 0 2length^2*Iy 0 0 0 6length*Iy 0 4length^2*Iy 0;
                0 6length*Iz 0 0 0 2length^2*Iz 0 -6length*Iz 0 0 0 4length^2*Iz
                ]

            element.k = Symmetric(R' * k * R)
        elseif dims == 2
            l, m = (element.posEnd .- element.posStart) ./ element.length

            # rotation matrix
            R = [l m 0 0 0 0;
                -m l 0 0 0 0;
                0 0 1 0 0 0;
                0 0 0 l m 0;
                0 0 0 -m l 0;
                0 0 0 0 0 1]
            element.R = R

            # stiffness matrix
            k = zeros(6,6)

            p = [element.A / element.length, 0, 0, -element.A/element.length, 0, 0]
            v = [0, 12 * element.Iz/element.length^3, 6 * element.Iz/element.length^2, 0, -12element.Iz/element.length^3, 6element.Iz/element.length^2]
            m = [0, 6element.Iz/element.length^2, 4element.Iz/element.length, 0, -6element.Iz/element.length^2, 2element.Iz/element.length]

            k[1,:] = p
            k[2,:] = v
            k[3,:] = m
            k[4,:] = -p
            k[5,:] = -v
            k[6,:] = [0, 6element.Iz/element.length^2, 2element.Iz/element.length, 0, -6element.Iz/element.length^2, 4element.Iz/element.length]

            k .*= element.E

            element.k = Symmetric(R' * k * R)
        else
            error("dims must be 2 or 3")
        end
    else
        error("type must be :truss or :frame")
    end

    if return_k
        return element.k
    end
end

function K(elements::Array{Element}, nDOFS::Int)
    K = spzeros(nDOFS, nDOFS)

    for element in elements
        K[element.dofIndex, element.dofIndex] .+= element.k
    end

    return K
end

#scales displacement s.t. maximum displacement is 30% of the mean element length
function autoScaleFactor(structure::Structure)
    lengths = [e.length for e in structure.elements]
    if structure.dims == 2
        displacements = vcat([structure.U[n.globalIndex[1:2]] for n in structure.nodes]...)
    else
        displacements = vcat([structure.U[n.globalIndex[1:3]] for n in structure.nodes]...)
    end
    factor = 0.30 / maximum(abs.(displacements)) * mean(lengths)
    return round(factor)
end

function postProcess(structure::Structure; scaleFactor = :auto)

    ### create displaced nodes and elements (scaled)
    if scaleFactor == :auto
        scaleFactor = autoScaleFactor(structure)
    end

    positions = [node.position for node in structure.nodes]

    if structure.dims == 2
        displacements = vcat([structure.U[n.globalIndex[1:2]] for n in structure.nodes])
    else
        displacements = vcat([structure.U[n.globalIndex[1:3]] for n in structure.nodes])
    end

    newpositions = [positions[i] .+ scaleFactor .* displacements[i] for i = 1:length(structure.nodes)]

    structure.displacedNodes = Node.(newpositions)
    structure.displacedElements = [Element(structure.displacedNodes, element.nodeIndex, element.type) for element in structure.elements]

    ### Internal forces in local CS
    for element in structure.elements
        element.localForce = element.R * element.k * structure.U[element.dofIndex]
        element.globalForce = element.k * structure.U[element.dofIndex]
        element.axialForce = axialLoad(element)
    end

    ### Nodal displacements and reactions
    for node in structure.nodes
        node.disp = structure.U[node.globalIndex]
        if any(node.DOFS .== false)
            node.reaction = structure.reactions[node.globalIndex]
        end
    end
end

function reactions(structure::Structure)
    if isdefined(structure, :U) == false
        error("Perform analysis before solving reactions.")
    end

    # index of boundary condition DOFs
    fixedDOFS = findall(structure.DOFS .== false)

    # solve for loads at reactions
    reactions = structure.K[fixedDOFS, :] * structure.U

    # initialize
    structure.reactions = zeros(structure.nDOFS)

    # update
    structure.reactions[fixedDOFS] = reactions
end

function addNodeElements(elements::Vector{Element}, nodes::Vector{Node})
    for i = 1:length(elements)
        j, k = elements[i].nodeIndex
        push!(nodes[j].elements, (i, -1))
        push!(nodes[k].elements, (i, 1))
    end

    if isdefined(nodes[1], :elements)
        for i = 1:length(nodes)
            nodes[i].elements = unique(nodes[i].elements)
        end
    end
end

function addNodeLoads(loads::Vector{Load}, nodes::Vector{Node})
   for i = 1:length(loads)
    nodes[loads[i].index].load = loads[i].load
   end
end

function analyze(structure::Structure; forceK = false, SF = :auto)
    if !isdefined(structure, :F)
        error("No loads defined.")
    end

    # associated elements to each node
    addNodeElements(structure.elements, structure.nodes)
    # associated load vectors to individual nodes
    addNodeLoads(structure.loads, structure.nodes)

    # prevents rebuilding of stiffness matrix if structure stays the same
    # setting forceK = true rebuilds the stiffness matrix
    if isdefined(structure, :K) && !forceK
        #displacement
        U = structure.K[structure.freeDOFS, structure.freeDOFS] \ structure.F[structure.freeDOFS]

        #compliance
        structure.compliance = U' * structure.F[structure.freeDOFS]
        return
    end
    #create nodal DOF indices
    nodeGlobalIndex(structure.nodes)

    #create elemental DOF indices
    dofExpander(structure.elements, structure.nodes)

    #create elemental stiffness matrices
    k_elemental.(structure.elements, structure.dims)

    #create global stiffness matrix
    structure.K = K(structure.elements, structure.nDOFS)

    #displacement
    U = structure.K[structure.freeDOFS, structure.freeDOFS] \ structure.F[structure.freeDOFS]

    #compliance
    structure.compliance = U' * structure.F[structure.freeDOFS]

    #long form displacement
    structure.U = zeros(structure.nDOFS)
    structure.U[structure.freeDOFS] = U

    #reactions
    reactions(structure)

    #post processing
    postProcess(structure; scaleFactor = SF)
end

function axialLoad(element::Element)
    if isdefined(element, :globalForce) == false
        error("No internal loads defined; run analysis and/or postprocessing")
    end

    dims = length(element.posStart)

    if element.type == :truss
        if dims == 2
            if dot(element.posEnd - element.posStart, element.globalForce[3:4]) < 0
                return - norm(element.globalForce[3:4])
            else
                return norm(element.globalForce[3:4])
            end
        else
            if dot(element.posEnd - element.posStart, element.globalForce[4:6]) < 0
                return - norm(element.globalForce[4:6])
            else
                return norm(element.globalForce[4:6])
            end
        end
    else
        if dims == 2
            if dot(element.posEnd - element.posStart, element.globalForce[4:5]) < 0
                return - norm(element.globalForce[4:5])
            else
                return norm(element.globalForce[4:5])
            end
        else
            if dot(element.posEnd - element.posStart, element.globalForce[7:9]) < 0
                return - norm(element.globalForce[7:9])
            else
                return norm(element.globalForce[7:9])
            end
        end
    end
end

```
Typical DOF settings:
dims : 2 or 3
type : :truss or :frame
fixity : :x/y/zfree, x/y/zfixed, :free, :fixed, :pinfixed
```
function dofMaker(dims::Int64, type::Symbol, fixity::Symbol)
    if dims == 2
        if type == :truss
            if fixity == :free
                return [true, true]
            elseif fixity == :fixed
                return [false, false]
            elseif fixity == :xfree
                return [true, false]
            elseif fixity == :yfree
                return [false, true]
            else
                error("unknown fixity, use :free, :fixed, :xfree, :yfree")
            end
        elseif type == :frame
            if fixity == :free
                return [true, true, true]
            elseif fixity == :pinfixed
                return [false, false, true]
            elseif fixity == :fixed
                return [false, false, false]
            elseif fixity == :xfree
                return [true, false, true]
            elseif fixity == :yfree
                return [false, true, false]
            else
                error("unknown fixity, use :free, :pinfixed, :fixed, :xfree, :yfree")
            end
        else
            error("type must be :truss or :frame")
        end
    elseif dims == 3
        if type == :truss
            if fixity == :free
                return [true, true, true]
            elseif fixity == :fixed
                return [false, false, false]
            elseif fixity == :zfixed
                return [true, true, false]
            elseif fixity == :xfixed
                return [false, true, true]
            elseif fixity == :yfixed
                return [true, false, true]
            elseif fixity == :xfree
                return [true, false, false]
            elseif fixity == :yfree
                return [false, true, false]
            elseif fixity == :zfree
                return [false, false, true]
            else
                error("unknown fixity, use :free, :fixed, :x/y/zfixed, :x/y/zfree")
            end
        elseif type == :frame
            if fixity == :free
                return [true, true, true, true, true, true]
            elseif fixity == :pinfixed
                return [false, false, false, true, true, true]
            elseif fixity == :fixed
                return [false, false, false, false, false, false]
            elseif fixity == :zfixed
                return [true, true, false, true, true, true]
            elseif fixity == :xfixed
                return  [false, true, true, true, true, true]
            elseif fixity == :yfixed
                return [true, false, true, true, true, true]
            elseif fixity == :xfree
                return [true, false, false, true, true, true]
            elseif fixity == :yfree
                return [false, true, false, true, true, true]
            elseif fixity == :zfree
                return [false, false, true, true, true, true]
            else
                error("unknown fixity, use :free, :pinfixed, :fixed, :x/y/zfixed, :x/y/zfree")
            end
        else
            error("type must be :truss or :frame")
        end
    else
        error("dims must be 2 or 3")
    end
end

function trussSizer(P, k, L, E, ρ, σ, minA)

    Astress = P/σ

    if sign(P) < 0
        D = (64 * abs(P) * (k * L)^2 / pi^3 / E / (1-ρ^4))^(1/4)
        Abuckling = pi/4 * D^2 *(1-ρ^2)
        if !isnothing(minA)
            A = max(Abuckling, Astress, minA)
        else
            A = max(Abuckling, Astress)
        end
    else
        if !isnothing(minA)
            A = max(Astress, minA)
        else
            A = Astress
        end
    end

    return A
end

# pseudo size is an approximated optimization function for sizing member areas. It performs a linear elastic analysis using uniform section properties, then resizes members to meet stress and buckling limitiations
function pseudoSize!(structure::Structure,
        maxStress;
        minArea = nothing,
        dDratio = 0.9,
        k = 1,
        maxIters = 1)

    # perform analysis

    analyze(structure)

    iter = 0
    while iter < maxIters
        for element in structure.elements
            element.A = trussSizer(element.axialForce, k, element.length, element.E, dDratio, maxStress, minArea)
        end
        analyze(structure; forceK = true)
        iter += 1
    end
end

function pseudoSize(structure::Structure,
    maxStress;
    minArea = nothing,
    dDratio = 0.9,
    k = 1,
    maxIters = 1)

    # perform analysis
    newstructure = deepcopy(structure)
    analyze(newstructure)

    iter = 0
    while iter < maxIters
        if storeAreas
            areas = Vector{Float64}()
        end
        for element in newstructure.elements
            element.A = trussSizer(element.axialForce, k, element.length, element.E, dDratio, maxStress, minArea)
            if storeAreas
                push!(areas, element.A)
            end
        end
        analyze(newstructure; forceK = true)

        iter += 1
    end

    return newstructure
end

function structureMass(structure::Structure, ρ)
    volume = 0.0
    for element in structure.elements
        volume += element.A * element.length
    end
    return volume * ρ
end

max_stars_repo_name: keithjlee/Asap

max_issues_repo_name: keithjlee/Asap

max_forks_repo_name: keithjlee/Asap


text:
"""
    genr_box_diag(data::Dict)
This will return the matrix of box size information.
"""
function genr_box_diag(data::Dict)
    box_diag = copy_array(data["box_info"])
    box_diag = mean(box_diag, dims=1)[1, :, :]
    box_diag = box_diag[:, 2] - box_diag[:, 1]
    box_diag = diag(box_diag)
    return box_diag
end

"""
    genr_box_vol(data::Dict)
This will return the volume of box in unit of your chossen.
"""
function genr_box_vol(data::Dict)
    box_vol = copy_array(data["box_info"])
    box_vol = mean(box_vol, dims=1)[1, :, :]
    box_vol = box_vol[:, 2] - box_vol[:, 1]
    return prod(box_vol)
end

"""
    genr_box_inv(data::Dict)
This will return the inverse matrix of box size information.
"""
function genr_box_inv(data::Dict)
    box_diag = genr_box_diag(data)
    return inv(box_diag)
end

"""
      PBC!(data::Dict)
This will add a variable "coord_corr" to `data` which is created by `read_dump`. `coord_corr` has handled the PBC issues.
"""
function PBC!(data::Dict)
    coord = data["coord"]
    box_diag = genr_box_diag(data)
    box_inv = genr_box_inv(data)
    coord_corr = zeros(size(coord))
    for step = 1 : data["num_steps"]
        coord_corr[step, :, :] = coord[step, :, :] * box_inv
    end
    coord_diff = coord_corr[2:end, :, :] .- coord_corr[1:end-1, :, :]
    coord_diff .-= round.(coord_diff)
    coord_corr[2:end, :, :] .= 0
    coord_corr[1, :, :] *= box_diag
    for step = 2 : data["num_steps"]
        coord_corr[step, :, :] = coord_corr[step-1, :, :] .+ coord_diff[step-1, :, :] * box_diag
    end
    data["coord_corr"] = coord_corr;

    println("\"coord_corr\"\thas been added to `data`")
end


text:
module DiffEqOperators

import Base: +, -, *, /, \, size, getindex, setindex!, Matrix, convert
using DiffEqBase, StaticArrays, LinearAlgebra
import LinearAlgebra: mul!, ldiv!, lmul!, rmul!, axpy!, opnorm, factorize, I
import DiffEqBase: AbstractDiffEqLinearOperator, update_coefficients!, is_constant
using SparseArrays

abstract type AbstractDerivativeOperator{T} <: AbstractDiffEqLinearOperator{T} end
abstract type AbstractDiffEqCompositeOperator{T} <: AbstractDiffEqLinearOperator{T} end

### Common default methods for the operators
include("common_defaults.jl")

### Basic Operators
include("basic_operators.jl")

### Derivative Operators
include("derivative_operators/fornberg.jl")
include("derivative_operators/upwind_operator.jl")
include("derivative_operators/derivative_irreg_operator.jl")
include("derivative_operators/derivative_operator.jl")
include("derivative_operators/abstract_operator_functions.jl")
include("derivative_operators/boundary_operators.jl")

### Composite Operators
include("composite_operators.jl")

# The (u,p,t) and (du,u,p,t) interface
for T in [DiffEqScalar, DiffEqArrayOperator, FactorizedDiffEqArrayOperator, DiffEqIdentity,
  DiffEqScaledOperator, DiffEqOperatorCombination, DiffEqOperatorComposition]
  (L::T)(u,p,t) = (update_coefficients!(L,u,p,t); L * u)
  (L::T)(du,u,p,t) = (update_coefficients!(L,u,p,t); mul!(du,L,u))
end

export DiffEqScalar, DiffEqArrayOperator, DiffEqIdentity, getops
export AbstractDerivativeOperator, DerivativeOperator, UpwindOperator, FiniteDifference
end # module


text:
for (gels, gesv, getrs, getri, elty) in
    ((:dgels_,:dgesv_,:dgetrs_,:dgetri_,:Float64),
     (:sgels_,:sgesv_,:sgetrs_,:sgetri_,:Float32),
     (:zgels_,:zgesv_,:zgetrs_,:zgetri_,:ComplexF64),
     (:cgels_,:cgesv_,:cgetrs_,:cgetri_,:ComplexF32))
    @eval begin
        #      SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,INFO)
        # *     .. Scalar Arguments ..
        #       CHARACTER          TRANS
        #       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS
        function mygels!(res::AbstractMatrix{$elty}, A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty}, work::Vector{$elty})
            trans = 'N'
            Base.require_one_based_indexing(A, B)
            LAPACK.chktrans(trans)
            chkstride1(A, B)
            btrn  = trans == 'T'
            m, n  = size(A)
            if size(B,1) != (btrn ? n : m)
                throw(DimensionMismatch("matrix A has dimensions ($m,$n), transposed: $btrn, but leading dimension of B is $(size(B,1))"))
            end
            info  = Ref{BlasInt}()
            lwork = BlasInt(length(work))
            ccall((@blasfunc($gels), liblapack), Cvoid,
                      ( Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
                      (btrn ? 'T' : 'N'), m, n, size(B,2), A, max(1,stride(A,2)),
                      B, max(1,stride(B,2)), work, lwork, info, 1)
            LAPACK.chklapackerror(info[])
            k   = min(m, n)
            res .= @view(A[1:k, 1:k])
            m < n ? tril!(res) : triu!(res)
            res
        end

        function calc_gels_working_size(A::AbstractMatrix{$elty}, B::AbstractVecOrMat{$elty})
            trans = 'N'
            Base.require_one_based_indexing(A, B)
            LAPACK.chktrans(trans)
            chkstride1(A, B)
            btrn  = trans == 'T'
            m, n  = size(A)
            if size(B,1) != (btrn ? n : m)
                throw(DimensionMismatch("matrix A has dimensions ($m,$n), transposed: $btrn, but leading dimension of B is $(size(B,1))"))
            end
            info  = Ref{BlasInt}()
            work  = Vector{$elty}(undef, 1)
            lwork = BlasInt(-1)
            ccall((@blasfunc($gels), liblapack), Cvoid,
                      ( Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                       Ptr{$elty}, Ref{BlasInt}, Ptr{BlasInt}, Clong),
                      (btrn ? 'T' : 'N'), m, n, size(B,2), A, max(1,stride(A,2)),
                      B, max(1,stride(B,2)), work, lwork, info, 1)
            LAPACK.chklapackerror(info[])
            BlasInt(real(work[1]))
        end
    end
end

"""
    mygels!(trans, A, B) -> (F)
Solves the linear equation `A * X = B`, `transpose(A) * X = B`, or `adjoint(A) * X = B` using
a QR or LQ factorization. Modifies the matrix/vector `B` in place with the
solution. `A` is overwritten with its `QR` or `LQ` factorization. `trans`
may be one of `N` (no modification), `T` (transpose), or `C` (conjugate
transpose). `gels!` searches for the minimum norm/least squares solution.
`A` may be under or over determined. The solution is returned in `B`.
"""
mygels!(trans::AbstractChar, A::AbstractMatrix, B::AbstractVecOrMat, work::Vector)

"""
    mygels!(trans, A, B) -> (working_size)
Calculates working size for mygels!
"""
calc_gels_working_size(trans::AbstractChar, A::AbstractMatrix, B::AbstractVecOrMat)