text:
mutable struct PriorParameters
    a::AbstractVector
    b::AbstractVector
    α::AbstractVector
    PriorParameters() = new()
    PriorParameters(a, b, α) = new(a, b, α)
end

mutable struct PoissionMixtureModel
    model::MixtureModel{Univariate,Discrete,Poisson,Float64}
    λs::AbstractMatrix
    ps::AbstractMatrix
    loglikelihoods::AbstractVector
    PoissionMixtureModel() = new()
end

function init_model(;K::Int = 3, a = 1.0, b = 0.01)
    α = init_α(K)
    λ = init_λ(a, b, K)
    p = init_p(α)
    prior = PriorParameters(ones(K) .* a, ones(K) .*b, α)
    prior, MixtureModel(Poisson, λ, p)
end

function generate_toy_data(
    N::Int = 10^4;
    true_λ::AbstractVector = [100.0, 50.0, 15.0],
    true_p::AbstractVector = [0.5, 0.2, 0.3],
)
    true_data = MixtureModel(Poisson, true_λ, true_p)
    toy_data = rand(true_data, N)
    true_data, toy_data
end

init_α(K) = rand(Uniform(0, 100), K)
init_λ(a, b, K) = rand(Gamma(a, 1 / b), K)
init_p(α) = rand(Dirichlet(α))

function sample_s(model, data)
    N = length(data)
    K = length(model.components)

    sₙ = Matrix{Int}(undef, N, K)
    η = Matrix{Float64}(undef, N, K)
    λ = [Distributions.params(comp)[1] for comp in model.components]
    p = probs(model)
    @inbounds for n = 1:N
        @simd for k = 1:K
            η[n, k] = data[n] * log(λ[k]) - λ[k] + log(p[k])
        end
        logsum_η = -logsumexp(η[n, :])
        @simd for k = 1:K
            η[n, k] = exp(η[n, k] + logsum_η)
        end
        sₙ[n, :] .= rand(Multinomial(1, η[n, :]))
    end
    sₙ
end

function sample_λ(model, prior, sₙ, toy_data)
    K = size(sₙ, 2)
    λ = Vector{Float64}(undef, K)
    @inbounds for k = 1:K
        â = sum(sₙ[:, k] .* toy_data) + prior.a[k]
        b̂ = sum(sₙ[:, k]) + prior.b[k]
        λ[k] = rand(Gamma(â, 1/b̂))
    end
    λ
end

function sample_p(model, prior, sₙ)
    N, K = size(sₙ)
    α̂ = Vector{Float64}(undef, K)
    @inbounds for k = 1:K
        α̂[k] = sum(sₙ[:, k]) + prior.α[k]
    end
    rand(Dirichlet(α̂))
end

function gibbs!(
    pmm::PoissionMixtureModel,
    prior::PriorParameters,
    toy_data::AbstractVector;
    maxiter::Int = 2000,
)
    model = pmm.model
    N = length(toy_data)
    K = length(model.components)

    pmm.λs = Matrix{Float64}(undef, maxiter, K)
    pmm.ps = Matrix{Float64}(undef, maxiter, K)
    pmm.loglikelihoods = Vector{Float64}(undef, maxiter)

    for k = 1:K
        pmm.λs[1, k] = model.components[k].λ
        pmm.ps[1, k] = probs(model)[k]
    end
    pmm.loglikelihoods[1] = sum(logpdf.(model, toy_data))

    @inbounds @simd for i = 1:maxiter
        sₙ = sample_s(model, toy_data)
        λ = sample_λ(model, prior, sₙ, toy_data)
        p = sample_p(model, prior, sₙ)

        pmm.λs[i, :] .= λ
        pmm.ps[i, :] .= p
        model = MixtureModel(Poisson, λ, p)
        pmm.loglikelihoods[i] = sum(logpdf.(model, toy_data))
        pmm.model = model
        fill!(sₙ, 0)
    end
    pmm
end

function update_s(prior, data)
    N = length(data)
    K = length(prior.α)

    η = Matrix{Float64}(undef, N, K)

    a = prior.a
    b = prior.b
    α = prior.α
    λ = a ./ b
    lnλ = digamma.(a) .- digamma.(b)
    lnp = [digamma(α[k]) + digamma(sum(α)) for k = 1:K]

    @inbounds for n = 1:N
        @simd for k = 1:K
            η[n, k] = data[n] * lnλ[k] - λ[k] + lnp[k]
        end
        logsum_η = -logsumexp(η[n, :])
        @simd for k = 1:K
            η[n, k] = exp(η[n, k] + logsum_η)
        end
    end
    η
end

function update_λ(prior, η, toy_data)
    K = size(η, 2)
    â = Vector{Float64}(undef, K)
    b̂ = Vector{Float64}(undef, K)
    a = prior.a
    b = prior.b
    @inbounds for k = 1:K
        â[k] = sum(η[:, k] .* toy_data) + a[k]
        b̂[k] = sum(η[:, k]) + b[k]
    end
    â, b̂
end

function update_p(prior, η)
    N, K = size(η)
    α̂ = Vector{Float64}(undef, K)
    α = prior.α
    @inbounds for k = 1:K
        α̂[k] = sum(η[:, k]) + α[k]
    end
    α̂
end

function vi!(prior::PriorParameters, toy_data::AbstractVector; maxiter::Int = 2000)
    N = length(toy_data)
    K = length(prior.α)

    #     elbo = Vector{Float64}(undef, N)
    #     elbo[1] =
    @inbounds for i = 1:maxiter
        η = update_s(prior, toy_data)
        â, b̂ = update_λ(prior, η, toy_data)
        α̂ = update_p(prior, η)

        prior.a = â
        prior.b = b̂
        prior.α = α̂
    end
end

max_stars_repo_name: SosUts/SuyamaBayes

max_issues_repo_name: SosUts/SuyamaBayes

max_forks_repo_name: SosUts/SuyamaBayes


text:
using FeatureSelectors
using Test
using RDatasets

@testset "FeatureSelectors.jl" begin
    @testset "Univariate feature selector" begin
        include("UnivariateFeatureSelector.jl")
    end
end

ext: jl


text:
## 00.jl
# 文字の逆順

str = "stressed"
println(reverse(str))


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

"""
    create_scaling_plan(application_source, scaling_instructions, scaling_plan_name)
    create_scaling_plan(application_source, scaling_instructions, scaling_plan_name, params::Dict{String,<:Any})

Creates a scaling plan.

# Arguments
- `application_source`: A CloudFormation stack or set of tags. You can create one scaling
  plan per application source. For more information, see ApplicationSource in the AWS Auto
  Scaling API Reference.
- `scaling_instructions`: The scaling instructions. For more information, see
  ScalingInstruction in the AWS Auto Scaling API Reference.
- `scaling_plan_name`: The name of the scaling plan. Names cannot contain vertical bars,
  colons, or forward slashes.

"""
create_scaling_plan(ApplicationSource, ScalingInstructions, ScalingPlanName; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("CreateScalingPlan", Dict{String, Any}("ApplicationSource"=>ApplicationSource, "ScalingInstructions"=>ScalingInstructions, "ScalingPlanName"=>ScalingPlanName); aws_config=aws_config)
create_scaling_plan(ApplicationSource, ScalingInstructions, ScalingPlanName, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("CreateScalingPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ApplicationSource"=>ApplicationSource, "ScalingInstructions"=>ScalingInstructions, "ScalingPlanName"=>ScalingPlanName), params)); aws_config=aws_config)

"""
    delete_scaling_plan(scaling_plan_name, scaling_plan_version)
    delete_scaling_plan(scaling_plan_name, scaling_plan_version, params::Dict{String,<:Any})

Deletes the specified scaling plan. Deleting a scaling plan deletes the underlying
ScalingInstruction for all of the scalable resources that are covered by the plan. If the
plan has launched resources or has scaling activities in progress, you must delete those
resources separately.

# Arguments
- `scaling_plan_name`: The name of the scaling plan.
- `scaling_plan_version`: The version number of the scaling plan. Currently, the only valid
  value is 1.

"""
delete_scaling_plan(ScalingPlanName, ScalingPlanVersion; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("DeleteScalingPlan", Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion); aws_config=aws_config)
delete_scaling_plan(ScalingPlanName, ScalingPlanVersion, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("DeleteScalingPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion), params)); aws_config=aws_config)

"""
    describe_scaling_plan_resources(scaling_plan_name, scaling_plan_version)
    describe_scaling_plan_resources(scaling_plan_name, scaling_plan_version, params::Dict{String,<:Any})

Describes the scalable resources in the specified scaling plan.

# Arguments
- `scaling_plan_name`: The name of the scaling plan.
- `scaling_plan_version`: The version number of the scaling plan. Currently, the only valid
  value is 1.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: The maximum number of scalable resources to return. The value must be
  between 1 and 50. The default value is 50.
- `"NextToken"`: The token for the next set of results.
"""
describe_scaling_plan_resources(ScalingPlanName, ScalingPlanVersion; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlanResources", Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion); aws_config=aws_config)
describe_scaling_plan_resources(ScalingPlanName, ScalingPlanVersion, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlanResources", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion), params)); aws_config=aws_config)

"""
    describe_scaling_plans()
    describe_scaling_plans(params::Dict{String,<:Any})

Describes one or more of your scaling plans.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ApplicationSources"`: The sources for the applications (up to 10). If you specify
  scaling plan names, you cannot specify application sources.
- `"MaxResults"`: The maximum number of scalable resources to return. This value can be
  between 1 and 50. The default value is 50.
- `"NextToken"`: The token for the next set of results.
- `"ScalingPlanNames"`: The names of the scaling plans (up to 10). If you specify
  application sources, you cannot specify scaling plan names.
- `"ScalingPlanVersion"`: The version number of the scaling plan. Currently, the only valid
  value is 1.  If you specify a scaling plan version, you must also specify a scaling plan
  name.
"""
describe_scaling_plans(; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlans"; aws_config=aws_config)
describe_scaling_plans(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("DescribeScalingPlans", params; aws_config=aws_config)

"""
    get_scaling_plan_resource_forecast_data(end_time, forecast_data_type, resource_id, scalable_dimension, scaling_plan_name, scaling_plan_version, service_namespace, start_time)
    get_scaling_plan_resource_forecast_data(end_time, forecast_data_type, resource_id, scalable_dimension, scaling_plan_name, scaling_plan_version, service_namespace, start_time, params::Dict{String,<:Any})

Retrieves the forecast data for a scalable resource. Capacity forecasts are represented as
predicted values, or data points, that are calculated using historical data points from a
specified CloudWatch load metric. Data points are available for up to 56 days.

# Arguments
- `end_time`: The exclusive end time of the time range for the forecast data to get. The
  maximum time duration between the start and end time is seven days.  Although this
  parameter can accept a date and time that is more than two days in the future, the
  availability of forecast data has limits. AWS Auto Scaling only issues forecasts for
  periods of two days in advance.
- `forecast_data_type`: The type of forecast data to get.    LoadForecast: The load metric
  forecast.     CapacityForecast: The capacity forecast.     ScheduledActionMinCapacity: The
  minimum capacity for each scheduled scaling action. This data is calculated as the larger
  of two values: the capacity forecast or the minimum capacity in the scaling instruction.
  ScheduledActionMaxCapacity: The maximum capacity for each scheduled scaling action. The
  calculation used is determined by the predictive scaling maximum capacity behavior setting
  in the scaling instruction.
- `resource_id`: The ID of the resource. This string consists of a prefix
  (autoScalingGroup) followed by the name of a specified Auto Scaling group (my-asg).
  Example: autoScalingGroup/my-asg.
- `scalable_dimension`: The scalable dimension for the resource. The only valid value is
  autoscaling:autoScalingGroup:DesiredCapacity.
- `scaling_plan_name`: The name of the scaling plan.
- `scaling_plan_version`: The version number of the scaling plan. Currently, the only valid
  value is 1.
- `service_namespace`: The namespace of the AWS service. The only valid value is
  autoscaling.
- `start_time`: The inclusive start time of the time range for the forecast data to get.
  The date and time can be at most 56 days before the current date and time.

"""
get_scaling_plan_resource_forecast_data(EndTime, ForecastDataType, ResourceId, ScalableDimension, ScalingPlanName, ScalingPlanVersion, ServiceNamespace, StartTime; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("GetScalingPlanResourceForecastData", Dict{String, Any}("EndTime"=>EndTime, "ForecastDataType"=>ForecastDataType, "ResourceId"=>ResourceId, "ScalableDimension"=>ScalableDimension, "ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion, "ServiceNamespace"=>ServiceNamespace, "StartTime"=>StartTime); aws_config=aws_config)
get_scaling_plan_resource_forecast_data(EndTime, ForecastDataType, ResourceId, ScalableDimension, ScalingPlanName, ScalingPlanVersion, ServiceNamespace, StartTime, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("GetScalingPlanResourceForecastData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EndTime"=>EndTime, "ForecastDataType"=>ForecastDataType, "ResourceId"=>ResourceId, "ScalableDimension"=>ScalableDimension, "ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion, "ServiceNamespace"=>ServiceNamespace, "StartTime"=>StartTime), params)); aws_config=aws_config)

"""
    update_scaling_plan(scaling_plan_name, scaling_plan_version)
    update_scaling_plan(scaling_plan_name, scaling_plan_version, params::Dict{String,<:Any})

Updates the specified scaling plan. You cannot update a scaling plan if it is in the
process of being created, updated, or deleted.

# Arguments
- `scaling_plan_name`: The name of the scaling plan.
- `scaling_plan_version`: The version number of the scaling plan. The only valid value is
  1. Currently, you cannot have multiple scaling plan versions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ApplicationSource"`: A CloudFormation stack or set of tags. For more information, see
  ApplicationSource in the AWS Auto Scaling API Reference.
- `"ScalingInstructions"`: The scaling instructions. For more information, see
  ScalingInstruction in the AWS Auto Scaling API Reference.
"""
update_scaling_plan(ScalingPlanName, ScalingPlanVersion; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("UpdateScalingPlan", Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion); aws_config=aws_config)
update_scaling_plan(ScalingPlanName, ScalingPlanVersion, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = auto_scaling_plans("UpdateScalingPlan", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ScalingPlanName"=>ScalingPlanName, "ScalingPlanVersion"=>ScalingPlanVersion), params)); aws_config=aws_config)


text:
module CreateTableRoles

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table

function up()
  create_table(:roles) do
    [
      primary_key()
      column(:name, :string, limit = 100)
    ]
  end

  add_index(:roles, :name)
end

function down()
  drop_table(:roles)
end

end


text:
# Tested in Julia 1.5.1
using DelimitedFiles
data = readdlm("./data/adventInput07.txt", '.')[:,1]
function bagsInDict(bagsList)
    bagsIn = Dict{String,String}()
    for b in bagsList
        if length(b) > 0
            nbags = b |> rstrip |> lstrip
            q = split(nbags, ' ')[1]
            bagsIn[nbags[length(q)+2:end]] = q
        end
    end
    bagsIn
end
function bagsOutDict(data)
    bagsOut = Dict{String,Dict}()
    for d in data
        pbag, mbag = split(d, r" contain ")
        bagsList = split(mbag, r"(bags|bag|,)")
        keyOut = split(pbag, r"(bags|bag|,)")[1] |> rstrip |> lstrip
        bagsOut[keyOut] = bagsInDict(bagsList)
    end
    bagsOut
end
# here all color bags with all sub bags...
bagColors = bagsOutDict(data)
# now what?...


text:
@testset "utils" begin
    data = LearningHorse.sample(1:100, 30)
    @test ones(size(data)...) == (data .<= 100)
    data, ncdata = LearningHorse.sample(1:100, 30, nc=true)
    @test ones(size(data)...) == (data .<= 100)
    @test ones(size(ncdata)) == (ncdata .<= 100)
end


text:
function implicit_keys()::Bool
    CODES = Dict("KEY" => 1)
    return "KEY" ∈ keys(CODES)
end

function explicit_keys()::Bool
    CODES = Dict("KEY" => 1)
    return "KEY" ∈ keys(CODES)
end

function dict_values()::Bool
    CODES = Dict("KEY" => 1)
    return 1 ∈ values(CODES)
end

function return_dict_index_str(key::String)::Int64
    CODES = Dict("KEY" => 1)
    return CODES[key]
end

function return_dict_index_int(key::Int64)::String
    CODES = Dict(1 => "one")
    return CODES[key]
end

if abspath(PROGRAM_FILE) == @__FILE__
    @assert(implicit_keys())
    @assert(explicit_keys())
    @assert(dict_values())
    @assert(return_dict_index_str("KEY") == 1)
    @assert(return_dict_index_int(1) == "one")
    println("OK")
end


text: println(gethostname())

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export HParallelotope,
       directions,
       offset,
       dim,
       base_vertex,
       extremal_vertices,
       center,
       genmat,
       generators,
       constraints_list

"""
    HParallelotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N}

Type that represents a parallelotope in constraint form.

### Fields

- `directions` -- square matrix such that each row is the direction of two parallel constraints
- `offset`     -- vector such that each element is the offset of the corresponding constraint

### Notes

Parallelotopes are centrally symmetric convex polytopes in ``\\mathbb{R}^n`` having
``2n`` pairwise parallel constraints. Every parallelotope is a zonotope.
As such, parallelotopes can be represented either in constraint form or in
generator form. The `HParallelotope` type represents parallelotopes in constraint form.

The constraint form of a parallelotope is described next. Let
``D ∈ \\mathbb{R}^{n × n}`` be a matrix and let ``c ∈ \\mathbb{R}^{2n}`` be
a vector. The parallelotope ``P ⊂ \\mathbb{R}^n`` generated by the directions
matrix ``D`` and the offset vector ``c`` is given by the set of points
``x ∈ \\mathbb{R}^`` such that:

```math
    D_i ⋅ x ≤ c_{i},\\text{  and  } -D_i ⋅ x ≤ c_{n+i}
```
for ``i = 1, …, n``. Here ``D_i`` represents the ``i``-th row of ``D`` and ``c_i`` the ``i``-th
component of ``c``.

For details on the notions given in these notes as well as applications of parallelotopes in
reachability analysis we refer to [1] and [2]. For conversions between set representations
we refer to [3].

### References

[1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for polynomial dynamical systems.*
    Formal Methods in System Design 50.1 (2017): 1-38.

[2] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Parallelotope bundles for polynomial reachability.*
    Proceedings of the 19th International Conference on Hybrid Systems: Computation and Control. ACM, 2016.

[3] Matthias Althoff, Olaf Stursberg, and Martin Buss. *Computing reachable sets of hybrid systems using
    a combination of zonotopes and polytopes.* Nonlinear analysis: hybrid systems 4.2 (2010): 233-249.
"""
struct HParallelotope{N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}} <: AbstractZonotope{N}
    directions::MN
    offset::VN

    # default constructor with dimension check
    function HParallelotope(D::MN, c::VN) where {N, VN<:AbstractVector{N}, MN<:AbstractMatrix{N}}
        @assert length(c) == 2*checksquare(D) "the length of the offset direction should be twice the size " *
            "of the directions matrix, but they are $(length(c)) and $(size(D)) dimensional respectively"
        return new{N, VN, MN}(D, c)
    end
end

isoperationtype(::Type{<:HParallelotope}) = false
isconvextype(::Type{<:HParallelotope}) = true

# =================
# Getter functions
# =================

"""
    directions(P::HParallelotope)

Return the directions matrix of a parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

A matrix where each row represents a direction of the parallelotope.
The negated directions `-D_i` are implicit (see [`HParallelotope`](@ref) for details).
"""
function directions(P::HParallelotope)
    return P.directions
end

"""
    offset(P::HParallelotope)

Return the offsets of a parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

A vector with the ``2n`` offsets of the parallelotope,
if ``n`` is the dimension of ``P``.
"""
function offset(P::HParallelotope)
    return P.offset
end

# ========================================
# Functions required by the set interface
# ========================================

"""
    dim(P::HParallelotope)

Return the dimension of a parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

An integer representing the dimension of the parallelotope
"""
function dim(P::HParallelotope)
    return size(P.directions, 1)
end

# ========================
# Type-specific functions
# ========================

"""
    base_vertex(P::HParallelotope)

Compute the base vertex of the given parallelotope in constraint representation.

### Input

- `P` -- parallelotope in constraint representation

### Output

The base vertex of ``P``.

### Algorithm

Intuitively, the base vertex is the point from which we get the relative
positions of all the other points.
The base vertex can be computed as the solution of the ``n``-dimensional linear
system ``D_i x = c_{n+i}`` for ``i = 1, \\ldots, n``, see [1, Section 3.2.1].

[1] Dreossi, Tommaso, Thao Dang, and Carla Piazza. *Reachability computation for
    polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38.
"""
function base_vertex(P::HParallelotope)
    D, c = P.directions, P.offset
    n = dim(P)
    v = to_negative_vector(view(c, (n+1):2n)) # converts to a dense vector as well
    return D \ v
end

"""
    extremal_vertices(P::HParallelotope{N, VN}) where {N, VN}

Compute the extremal vertices with respect to the base vertex of the
given parallelotope in constraint representation.

### Input

- `P` -- parallelotope in constraint representation

### Output

The list of vertices connected to the base vertex of ``P``.

### Algorithm

Let ``P`` be a parallelotope in constraint representation with directions
matrix ``D`` and offset vector ``c``. We denote the *extremal vertices* of ``P``
with respect to its base vertex ``q`` to be those vertices of ``P`` which have
an edge in common with ``q``.  The extremal vertices can be computed as the
solution of the ``n``-dimensional linear systems of equations ``D x = v_i``
where for each ``i = 1, \\ldots, n``,
``v_i = [-c_{n+1}, \\ldots, c_i, \\ldots, -c_{2n}]``.

We refer to [1, Section 3.2.1] for details.

[1] Tommaso Dreossi, Thao Dang, and Carla Piazza. *Reachability computation for
    polynomial dynamical systems.* Formal Methods in System Design 50.1 (2017): 1-38.
"""
function extremal_vertices(P::HParallelotope{N, VN}) where {N, VN}
    D, c = P.directions, P.offset
    n = dim(P)
    v = to_negative_vector(view(c, n+1:2n))
    vertices = Vector{VN}(undef, n)
    h = copy(v)
    @inbounds for i in 1:n
        h[i] = c[i]
        vertices[i] = D \ h
        h[i] = v[i]
    end
    return vertices
end

"""
    center(P::HParallelotope)

Return the center of a parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

The center of the parallelotope.

### Algorithm

Let ``P`` be a parallelotope with base vertex ``q`` and list of extremal
vertices with respect to ``q`` given by the set ``\\{v_i\\}`` for
``i = 1, \\ldots, n``. Then the center is located at

```math
    c = q + \\sum_{i=1}^n \\frac{v_i - q}{2} = q (1 - \\frac{2}) + \\frac{s}{2},
```
where ``s := \\sum_{i=1}^n v_i`` is the sum of extremal vertices.
"""
function center(P::HParallelotope)
    n = dim(P)
    q = base_vertex(P)
    E = extremal_vertices(P)
    s = sum(E)
    return q * (1 - n/2) + s/2
end

"""
    genmat(P::HParallelotope)

Return the generator matrix of a parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

A matrix where each column represents one generator of the parallelotope `P`.

### Algorithm

Let ``P`` be a parallelotope with base vertex ``q`` and list of extremal
vertices with respect to ``q`` given by the set ``\\{v_i\\}`` for
``i = 1, \\ldots, n``. Then, the ``i``-th generator of ``P``, represented as the
``i``-th column vector ``G[:, i]``, is given by:

```math
    G[:, i] = \\frac{v_i - q}{2}
```
for ``i = 1, \\ldots, n``.
"""
function genmat(P::HParallelotope)
    E = extremal_vertices(P)
    q = base_vertex(P)
    return 1/2 * hcat(E...) .- q/2
end

"""
    generators(P::HParallelotope)

Return an iterator over the generators of a parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

An iterator over the generators of `P`.
"""
function generators(P::HParallelotope)
    return generators_fallback(P)
end

"""
    constraints_list(P::HParallelotope{N, VN}) where {N, VN}

Return the list of constraints of the given parallelotope.

### Input

- `P` -- parallelotope in constraint representation

### Output

The list of constraints of `P`.
"""
function constraints_list(P::HParallelotope{N, VN}) where {N, VN}
    D, c = P.directions, P.offset
    n = dim(P)
    if isempty(D)
        return Vector{LinearConstraint{N, VN}}()
    end
    clist = Vector{LinearConstraint{N, VN}}(undef, 2n)
    @inbounds for i in 1:n
        clist[i] = LinearConstraint(D[i, :], c[i])
        clist[i+n] = LinearConstraint(-D[i, :], c[i+n])
    end
    return clist
end

"""
    rand(::Type{HParallelotope}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random parallelotope.

### Input

- `HParallelotope` -- type for dispatch
- `N`             -- (optional, default: `Float64`) numeric type
- `dim`           -- (optional, default: 2) dimension
- `rng`           -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`          -- (optional, default: `nothing`) seed for reseeding

### Output

A random parallelotope.

### Notes

All numbers are normally distributed with mean 0 and standard deviation 1.
"""
function rand(::Type{HParallelotope};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    D = randn(N, dim, dim)
    offset = randn(N, 2 * dim)
    return HParallelotope(D, offset)
end


text:
using Smg2s
using Test
using SparseArrays
using Plots

using LinearAlgebra
using Random

n = 1000
diag_l = -800
diag_u = -790
nbOne = 10

nilpvec = zeros(Int64, n-1)

cnt = 0

while cnt < n-1
    l₁ = rand(1:nbOne)
    l₀ = rand(1:nbOne)

    for i = cnt+1:l₁+cnt
        if i <= n-1
            nilpvec[i] = 1
        end
    end

    global cnt += l₀ + l₁

end

nilpMatrix = sparse(diagm(100=>nilpvec[1:n-100]))

#@info nilpvec
#nilp = Nilp(nilpMatrix, n)

nilp = Nilp(nbOne, 100, n)

## No hermitian case
function f1(idx::Integer, n = n)
    return cos((idx-1) * 2 * π / n) + 1 * rand() + 10 + 10 * sin((idx-1) * 2 * π / n) * im
end

spec1 = zeros(ComplexF64, n)
Spectrum!(spec1, f1, n)

Am1 = spzeros(ComplexF64, n, n)
initMat!(Am1, diag_l, diag_u, n; scale=0.1, sparsity=0.5)

@time genMat1 = nonherm(n, diag_l, diag_u, spec1, Am1, nilp)
@info "sparsity = " nnz(genMat1) / (n * n)
genspec1 = eigen!(Matrix(genMat1)).values

p1=scatter(real(spec1), imag(spec1), markercolor = :green, markern = 7, label="given spectrum")
p1=scatter!(real(genspec1), imag(genspec1), markercolor = :red, markern = 5, label="spectrum of generated matrix")

## Non symmetric 1
function f2(idx::Integer, n = n)
    rnd = rand()
    if idx % 2 == 0
        return 10 * cos((idx) * 2 * π / n) + 5 * rnd + 5 * sin((idx) * 2 * π / n) * im
    else
        return 10 * cos((idx+1) * 2 * π / n) + 5 * rnd - 5 * sin((idx+1) * 2 * π / n) * im
    end
end

spec2 = zeros(ComplexF64, n)
Spectrum!(spec2, f2, n)

Am2 = spzeros(Float64, n, n)
initMat!(Am2, diag_l, diag_u, n; scale=0.1, sparsity=0.5)
@time genMat2 = nonsym(n, diag_l, diag_u, spec2, Am2, nilp)

@info "sparsity = " nnz(genMat2) / (n * n)

genspec2 = eigen!(Matrix(genMat2)).values
p2=scatter(real(spec2), imag(spec2), markercolor = :green, markern = 7, label="given spectrum")
p2=scatter!(real(genspec2), imag(genspec2), markercolor = :red, markern = 5, label="spectrum of generated matrix")

## Non symmetric 2
function f3(idx::Integer, n = n)
    return 10 * cos(((idx-1) / 2) * 2 * π / n) + 5
end

spec3 = zeros(ComplexF64, n)
Spectrum!(spec3, f3, n)

Am3 = spzeros(Float64, n, n)
initMat!(Am3, diag_l, diag_u, n; scale=0.1, sparsity=.5)
@time genMat3 = nonsym(n, diag_l, diag_u, spec3, Am3, nilp)

@info "sparsity = " nnz(genMat3) / (n * n)

genspec3 = eigen!(Matrix(genMat3)).values
p3=scatter(real(spec3), imag(spec3), ylims=(-1,1), markercolor = :green, markern = 7, label="given spectrum")
p3=scatter!(real(genspec3), imag(genspec3), ylims=(-1,1), markercolor = :red, markern = 5, label="spectrum of generated matrix")

p4=spy(dropzeros(genMat3), legend = :none)

l = @layout [a  b;c d]
plot(p1, p2, p3, p4, layout = l)

#savefig("fig/example")


text:
using RobotZoo
using RobotDynamics
using StaticArrays
using BenchmarkTools

linmodel = RobotZoo.LinearModels.DoubleIntegrator(3)
A = RobotDynamics.get_A(linmodel)
B = RobotDynamics.get_B(linmodel)
@test A isa SizedMatrix{6,6}
@test B isa SizedMatrix{6,3}
@test !RobotDynamics.is_discrete(linmodel)

model = RobotZoo.DoubleIntegrator(3)

F1 = RobotDynamics.DynamicsJacobian(linmodel)
F2 = RobotDynamics.DynamicsJacobian(model)
z = KnotPoint(rand(model)..., 0.1)
b1 = @benchmark dynamics($linmodel, $z)
b2 = @benchmark dynamics($model, $z)
@test minimum(b2) < minimum(b1)  # standard model is faster since it avoids matrix mult

if !haskey(ENV, "CI")
    b1 = @benchmark jacobian!($F1, $linmodel, $z)
    b2 = @benchmark jacobian!($F2, $model, $z)
    @test minimum(b1) < minimum(b2)  # linear model Jacobian is faster
end


text:
using JuMP

function decision_variable(model::Model, S::States, d::DecisionNode, base_name::String="")
    # Create decision variables.
    dims = S[[d.I_j; d.j]]
    z_j = Array{VariableRef}(undef, dims...)
    for s in paths(dims)
        z_j[s...] = @variable(model, binary=true, base_name=base_name)
    end
    # Constraints to one decision per decision strategy.
    for s_I in paths(S[d.I_j])
        @constraint(model, sum(z_j[s_I..., s_j] for s_j in 1:S[d.j]) == 1)
    end
    return z_j
end

struct DecisionVariables
    D::Vector{DecisionNode}
    z::Vector{<:Array{VariableRef}}
end

"""
    DecisionVariables(model::Model, S::States, D::Vector{DecisionNode}; names::Bool=false, name::String="z")

Create decision variables and constraints.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `S::States`: States structure associated with the influence diagram.
- `D::Vector{DecisionNode}`: Vector containing decicion nodes.
- `names::Bool`: Use names or have JuMP variables be anonymous.
- `name::String`: Prefix for predefined decision variable naming convention.

# Examples
```julia
z = DecisionVariables(model, S, D)
```
"""
function DecisionVariables(model::Model, S::States, D::Vector{DecisionNode}; names::Bool=false, name::String="z")
    DecisionVariables(D, [decision_variable(model, S, d, (names ? "$(name)_$(d.j)$(s)" : "")) for d in D])
end

function is_forbidden(s::Path, forbidden_paths::Vector{ForbiddenPath})
    return !all(s[k]∉v for (k, v) in forbidden_paths)
end

function path_compatibility_variable(model::Model, z::DecisionVariables, base_name::String="")
    # Create a path compatiblity variable
    x = @variable(model, base_name=base_name)

    # Constraint on the lower and upper bounds.
    @constraint(model, 0 ≤ x ≤ 1.0)

    return x
end

struct PathCompatibilityVariables{N} <: AbstractDict{Path{N}, VariableRef}
    data::Dict{Path{N}, VariableRef}
end

Base.getindex(x_s::PathCompatibilityVariables, key) = getindex(x_s.data, key)
Base.get(x_s::PathCompatibilityVariables, key, default) = get(x_s.data, key, default)
Base.keys(x_s::PathCompatibilityVariables) = keys(x_s.data)
Base.values(x_s::PathCompatibilityVariables) = values(x_s.data)
Base.pairs(x_s::PathCompatibilityVariables) = pairs(x_s.data)
Base.iterate(x_s::PathCompatibilityVariables) = iterate(x_s.data)
Base.iterate(x_s::PathCompatibilityVariables, i) = iterate(x_s.data, i)

function decision_strategy_constraint(model::Model, S::States, d::DecisionNode, D::Vector{DecisionNode}, z::Array{VariableRef}, x_s::PathCompatibilityVariables)

    # states of nodes in information structure (s_j | s_I(j))
    dims = S[[d.I_j; d.j]]

    # Theoretical upper bound based on number of paths with information structure (s_j | s_I(j)) divided by number of possible decision strategies in other decision nodes
    other_decisions = map(d_n -> d_n.j, filter(d_n -> all(d_n.j != i for i in [d.I_j; d.j]), D))
    theoretical_ub = prod(S)/prod(dims)/ prod(S[other_decisions])

    # paths that have corresponding path compatibility variable
    existing_paths = keys(x_s)

    for s_j_s_Ij in paths(dims) # iterate through all information states and states of d
        # paths with (s_j | s_I(j)) information structure
        feasible_paths = filter(s -> s[[d.I_j; d.j]] == s_j_s_Ij, existing_paths)

        @constraint(model, sum(get(x_s, s, 0) for s in feasible_paths) ≤ z[s_j_s_Ij...] * min(length(feasible_paths), theoretical_ub))
    end
end

"""
    PathCompatibilityVariables(model::Model,
        z::DecisionVariables,
        S::States,
        P::AbstractPathProbability;
        names::Bool=false,
        name::String="x",
        forbidden_paths::Vector{ForbiddenPath}=ForbiddenPath[],
        fixed::Dict{Node, State}=Dict{Node, State}(),
        probability_cut::Bool=true)

Create path compatibility variables and constraints.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `z::DecisionVariables`: Decision variables from `DecisionVariables` function.
- `S::States`: States structure associated with the influence diagram.
- `P::AbstractPathProbability`: Path probabilities structure for which the function
  `P(s)` is defined and returns the path probabilities for path `s`.
- `names::Bool`: Use names or have JuMP variables be anonymous.
- `name::String`: Prefix for predefined decision variable naming convention.
- `forbidden_paths::Vector{ForbiddenPath}`: The forbidden subpath structures.
    Path compatibility variables will not be generated for paths that include
    forbidden subpaths.
- `fixed::Dict{Node, State}`: Path compatibility variable will not be generated
    for paths which do not include these fixed subpaths.
- `probability_cut` Includes probability cut constraint in the optimisation model.

# Examples
```julia
x_s = PathCompatibilityVariables(model, z, S, P; probability_cut = false)
```
"""
function PathCompatibilityVariables(model::Model,
    z::DecisionVariables,
    S::States,
    P::AbstractPathProbability;
    names::Bool=false,
    name::String="x",
    forbidden_paths::Vector{ForbiddenPath}=ForbiddenPath[],
    fixed::Dict{Node, State}=Dict{Node, State}(),
    probability_cut::Bool=true)

    if !isempty(forbidden_paths)
        @warn("Forbidden paths is still an experimental feature.")
    end

    # Create path compatibility variable for each effective path.
    N = length(S)
    variables_x_s = Dict{Path{N}, VariableRef}(
        s => path_compatibility_variable(model, z, (names ? "$(name)$(s)" : ""))
        for s in paths(S, fixed)
        if !iszero(P(s)) && !is_forbidden(s, forbidden_paths)
    )

    x_s = PathCompatibilityVariables{N}(variables_x_s)

    # Add decision strategy constraints for each decision node
    for (d, z_d) in zip(z.D, z.z)
        decision_strategy_constraint(model, S, d, z.D, z_d, x_s)
    end

    if probability_cut
        @constraint(model, sum(x * P(s) for (s, x) in x_s) == 1.0)
    end

    x_s
end

"""
    lazy_probability_cut(model::Model, x_s::PathCompatibilityVariables, P::AbstractPathProbability)

Adds a probability cut to the model as a lazy constraint.

# Examples
```julia
lazy_probability_cut(model, x_s, P)
```

!!! note
    Remember to set lazy constraints on in the solver parameters, unless your solver does this automatically. Note that Gurobi does this automatically.

"""
function lazy_probability_cut(model::Model, x_s::PathCompatibilityVariables, P::AbstractPathProbability)
    # August 2021: The current implementation of JuMP doesn't allow multiple callback functions of the same type (e.g. lazy)
    # (see https://github.com/jump-dev/JuMP.jl/issues/2642)
    # What this means is that if you come up with a new lazy cut, you must replace this
    # function with a more general function (see discussion and solution in https://github.com/gamma-opt/DecisionProgramming.jl/issues/20)

    function probability_cut(cb_data)
        xsum = sum(callback_value(cb_data, x) * P(s) for (s, x) in x_s)
        if !isapprox(xsum, 1.0)
            con = @build_constraint(sum(x * P(s) for (s, x) in x_s) == 1.0)
            MOI.submit(model, MOI.LazyConstraint(cb_data), con)
        end
    end
    MOI.set(model, MOI.LazyConstraintCallback(), probability_cut)
end

# --- Objective Functions ---

"""
    PositivePathUtility(S::States, U::AbstractPathUtility)

Positive affine transformation of path utility. Always evaluates positive values.

# Examples
```julia-repl
julia> U⁺ = PositivePathUtility(S, U)
julia> all(U⁺(s) > 0 for s in paths(S))
true
```
"""
struct PositivePathUtility <: AbstractPathUtility
    U::AbstractPathUtility
    min::Float64
    function PositivePathUtility(S::States, U::AbstractPathUtility)
        u_min = minimum(U(s) for s in paths(S))
        new(U, u_min)
    end
end

(U::PositivePathUtility)(s::Path) = U.U(s) - U.min + 1

"""
    NegativePathUtility(S::States, U::AbstractPathUtility)

Negative affine transformation of path utility. Always evaluates negative values.

# Examples
```julia-repl
julia> U⁻ = NegativePathUtility(S, U)
julia> all(U⁻(s) < 0 for s in paths(S))
true
```
"""
struct NegativePathUtility <: AbstractPathUtility
    U::AbstractPathUtility
    max::Float64
    function NegativePathUtility(S::States, U::AbstractPathUtility)
        u_max = maximum(U(s) for s in paths(S))
        new(U, u_max)
    end
end

(U::NegativePathUtility)(s::Path) = U.U(s) - U.max - 1

"""
    expected_value(model::Model,
        x_s::PathCompatibilityVariables,
        U::AbstractPathUtility,
        P::AbstractPathProbability;
        probability_scale_factor::Float64=1.0)

Create an expected value objective.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `x_s::PathCompatibilityVariables`: Path compatibility variables.
- `S::States`: States structure associated with the influence diagram.
- `P::AbstractPathProbability`: Path probabilities structure for which the function
  `P(s)` is defined and returns the path probabilities for path `s`.
- `probability_scale_factor::Float64`: Multiplies the path probabilities by this factor.

# Examples
```julia
EV = expected_value(model, x_s, U, P)
EV = expected_value(model, x_s, U, P; probability_scale_factor = 10.0)
```
"""
function expected_value(model::Model,
    x_s::PathCompatibilityVariables,
    U::AbstractPathUtility,
    P::AbstractPathProbability;
    probability_scale_factor::Float64=1.0)

    if probability_scale_factor ≤ 0
        throw(DomainError("The probability_scale_factor must be greater than 0."))
    end

    @expression(model, sum(P(s) * x * U(s) * probability_scale_factor for (s, x) in x_s))
end

"""
    conditional_value_at_risk(model::Model,
        x_s::PathCompatibilityVariables{N},
        U::AbstractPathUtility,
        P::AbstractPathProbability,
        α::Float64;
        probability_scale_factor::Float64=1.0) where N

Create a conditional value-at-risk (CVaR) objective.

# Arguments
- `model::Model`: JuMP model into which variables are added.
- `x_s::PathCompatibilityVariables`: Path compatibility variables.
- `S::States`: States structure associated with the influence diagram.
- `P::AbstractPathProbability`: Path probabilities structure for which the function
  `P(s)` is defined and returns the path probabilities for path `s`.
- `α::Float64`: Probability level at which conditional value-at-risk is optimised.
- `probability_scale_factor::Float64`: Adjusts conditional value at risk model to
   be compatible with the expected value expression if the probabilities were scaled there.

# Examples
```julia
α = 0.05  # Parameter such that 0 ≤ α ≤ 1
CVaR = conditional_value_at_risk(model, x_s, U, P, α)
CVaR = conditional_value_at_risk(model, x_s, U, P, α; probability_scale_factor = 10.0)
```
"""
function conditional_value_at_risk(model::Model,
    x_s::PathCompatibilityVariables{N},
    U::AbstractPathUtility,
    P::AbstractPathProbability,
    α::Float64;
    probability_scale_factor::Float64=1.0) where N

    if probability_scale_factor ≤ 0
        throw(DomainError("The probability_scale_factor must be greater than 0."))
    end
    if !(0 < α ≤ 1)
        throw(DomainError("α should be 0 < α ≤ 1"))
    end

    if !(probability_scale_factor == 1.0)
        @warn("The conditional value at risk is scaled by the probability_scale_factor. Make sure other terms of the objective function are also scaled.")
    end

    # Pre-computed parameters
    u = collect(Iterators.flatten(U(s) for s in keys(x_s)))
    u_sorted = sort(u)
    u_min = u_sorted[1]
    u_max = u_sorted[end]
    M = u_max - u_min
    u_diff = diff(u_sorted)
    ϵ = if isempty(u_diff) 0.0 else minimum(filter(!iszero, abs.(u_diff))) / 2 end

    # Variables and constraints
    η = @variable(model)
    @constraint(model, η ≥ u_min)
    @constraint(model, η ≤ u_max)
    ρ′_s = Dict{Path{N}, VariableRef}()
    for (s, x) in x_s
        u_s = U(s)
        λ = @variable(model, binary=true)
        λ′ = @variable(model, binary=true)
        ρ = @variable(model)
        ρ′ = @variable(model)
        @constraint(model, η - u_s ≤ M * λ)
        @constraint(model, η - u_s ≥ (M + ϵ) * λ - M)
        @constraint(model, η - u_s ≤ (M + ϵ) * λ′ - ϵ)
        @constraint(model, η - u_s ≥ M * (λ′ - 1))
        @constraint(model, 0 ≤ ρ)
        @constraint(model, 0 ≤ ρ′)
        @constraint(model, ρ ≤ λ * probability_scale_factor)
        @constraint(model, ρ′ ≤ λ′* probability_scale_factor)
        @constraint(model, ρ ≤ ρ′)
        @constraint(model, ρ′ ≤ x * P(s) * probability_scale_factor)
        @constraint(model, (x * P(s) - (1 - λ))* probability_scale_factor ≤ ρ)
        ρ′_s[s] = ρ′
    end
    @constraint(model, sum(values(ρ′_s)) == α * probability_scale_factor)

    # Return CVaR as an expression
    CVaR = @expression(model, sum(ρ_bar * U(s) for (s, ρ_bar) in ρ′_s) / α)

    return CVaR
end

# --- Construct decision strategy from JuMP variables ---

"""
    LocalDecisionStrategy(j::Node, z::Array{VariableRef})

Construct decision strategy from variable refs.
"""
function LocalDecisionStrategy(j::Node, z::Array{VariableRef})
    LocalDecisionStrategy(j, @. Int(round(value(z))))
end

"""
    DecisionStrategy(z::DecisionVariables)

Extract values for decision variables from solved decision model.

# Examples
```julia
Z = DecisionStrategy(z)
```
"""
function DecisionStrategy(z::DecisionVariables)
    DecisionStrategy(z.D, [LocalDecisionStrategy(d.j, v) for (d, v) in zip(z.D, z.z)])
end


text:
using Base: @propagate_inbounds

struct DenseTaggedImage{T, N, O <: Unsigned,AA <: AbstractArray} <: AbstractDenseTIFF{T, N}
    data::AA
    ifds::Vector{IFD{O}}

    function DenseTaggedImage(data::AbstractArray{T, N}, ifds::Vector{IFD{O}}) where {T, N, O}
        if N == 3
            @assert size(data, 3) == length(ifds)
        elseif N == 2
            @assert length(ifds) == 1
        end
        new{T, N, O, typeof(data)}(data, ifds)
    end
end

function DenseTaggedImage(data::AbstractArray{T, 2}, ifd::IFD{O}) where {T, O}
    DenseTaggedImage(data, IFD{O}[ifd])
end

DenseTaggedImage(data::AbstractArray{T, 3}) where {T} = DenseTaggedImage(data, _constructifd(data))
DenseTaggedImage(data::AbstractArray{T, 2}) where {T} = DenseTaggedImage(data, [_constructifd(data, UInt32)])

Base.size(t::DenseTaggedImage) = size(t.data)

@propagate_inbounds Base.getindex(img::DenseTaggedImage{T, N}, i::Vararg{Int, N}) where {T, N} = img.data[i...]
@propagate_inbounds Base.getindex(img::DenseTaggedImage, i...) = getindex(img.data, i...)
@propagate_inbounds Base.getindex(img::DenseTaggedImage{T, 3}, i::Int, j::Int, k::Int) where {T} = getindex(img.data, i, j, k)

Base.setindex!(img::DenseTaggedImage, i...) = setindex!(img.data, i...)

function Base.getindex(img::DenseTaggedImage{T, 3}, i::Colon, j::Colon, k) where {T}
    DenseTaggedImage(getindex(img.data, i, j, k), img.ifds[k])
end

"""
    offset(img)

Returns the type of the TIFF file offset. For most TIFFs, it should be UInt32,
while for BigTIFFs it will be UInt64.
"""
offset(::DenseTaggedImage{T, N, O, AA}) where {T, N, O, AA} = O

"""
    _constructifd(data, offset)

Generate a IFD with the minimal set of tags to describe `data`.
"""
function _constructifd(data::AbstractArray{T, 3}) where {T <: Colorant}
    offset = UInt32
    if sizeof(T) * length(data) >= 3.8*10^9
        @info "Array too large to fit in standard TIFF container, switching to BigTIFF"
        offset = UInt64
    end

    ifds = IFD{offset}[]

    for slice in 1:size(data, 3)
        push!(ifds, _constructifd(view(data, :, :, slice), offset))
    end

    ifds
end

function _constructifd(data::AbstractArray{T, 2}, ::Type{O}) where {T <: Colorant, O <: Unsigned}
    ifd = IFD(O)

    ifd[IMAGEWIDTH] = UInt32(size(data, 2))
    ifd[IMAGELENGTH] = UInt32(size(data, 1))
    n_samples = samplesperpixel(data)
    ifd[BITSPERSAMPLE] = fill(UInt16(bitspersample(data)), n_samples)
    ifd[PHOTOMETRIC] = interpretation(data)
    ifd[SAMPLESPERPIXEL] = UInt16(n_samples)
    ifd[SAMPLEFORMAT] = fill(UInt16(sampleformat(data)), n_samples)
    extra = extrasamples(data)
    if extra !== nothing
        ifd[EXTRASAMPLES] = extra
    end
    ifd
end

Base.write(io::IOStream, img::DenseTaggedImage) = write(Stream(format"TIFF", io, extract_filename(io)), img)

function Base.write(io::Stream, img::DenseTaggedImage)
    O = offset(img)
    tf = TiffFile(O)
    tf.io = io

    prev_ifd_record = write(tf) # record that will have be updated

    for (idx, ifd) in enumerate(img.ifds)
        data_pos = position(tf.io) # start of data

        write(tf, reinterpret(UInt8, permutedims(view(img.data, :, :, idx), [2, 1]))) # write data
        ifd_pos = position(tf.io)

        # update record of previous IFD to point to this new IFD
        seek(tf.io, prev_ifd_record)
        write(tf, O(ifd_pos))

        ifd = first(img.ifds)
        ifd[COMPRESSION] = COMPRESSION_NONE
        ifd[STRIPOFFSETS] = O(data_pos)
        ifd[STRIPBYTECOUNTS] = O(ifd_pos-data_pos)
        ifd[SOFTWARE] = "$(parentmodule(IFD)).jl v$PKGVERSION"
        sort!(ifd.tags)

        seek(tf.io, ifd_pos)
        prev_ifd_record = write(tf, ifd)
        seekend(tf.io)
    end
end


text:
@adjoint real(x::Complex) = real(x), r̄ -> (r̄ + zero(r̄)*im,)
@adjoint imag(x::Complex) = imag(x), ī -> (zero(ī) + ī*im,)

# The adjoint of the map z -> g*z is given by y -> g' * y.
# Therefore, for holomorphic functions (for which the differential is given by a complex multiplication),
# the gradient map is given by a multiplication with the conjugate of the derivative (in the holomorphic sense)
@adjoint log(x::Complex) = log(x), ȳ -> (ȳ/conj(x),)
@adjoint exp(x::Complex) = exp(x), ȳ -> (ȳ*conj(exp(x)),)
@adjoint sin(x::Complex) = sin(x), ȳ -> (ȳ*conj(cos(x)),)
@adjoint cos(x::Complex) = cos(x), ȳ -> (-ȳ*conj(sin(x)),)


text:
using Test
using ShuffleProofs: marshal_s_Gq, unmarshal, decode

s_Gq = "ModPGroup(safe-prime modulus=2*order+1. order bit-length = 511)::00000000020100000020636f6d2e766572696669636174756d2e61726974686d2e4d6f645047726f757000000000040100000041009a91c3b704e382e0c772fa7cf0e5d6363edc53d156e841555702c5b6f906574204bf49a551b695bed292e0218337c0861ee649d2fe4039174514fe2c23c10f6701000000404d48e1db8271c17063b97d3e7872eb1b1f6e29e8ab7420aaab8162db7c832ba1025fa4d2a8db4adf69497010c19be0430f7324e97f201c8ba28a7f1611e087b3010000004100300763b0150525252e4989f51e33c4e6462091152ef2291e45699374a3aa8acea714ff30260338bddbb48fc7446b273aaada90e3ee8326f388b582ea8a073502010000000400000001"

g = unmarshal(BigInt, decode(split(s_Gq, "::")[2]))

@test convert(String, marshal_s_Gq(g)) == s_Gq


text:
@testset "read-write maps" begin
    mktempdir() do tmpdir_ro
        mktempdir() do tmpdir_rw
            open(joinpath(tmpdir_ro, "myroinputfile"), "w") do io
                println(io, "Hello RO world!")
            end
            open(joinpath(tmpdir_rw, "myrwinputfile"), "w") do io
                println(io, "Hello RW world!")
            end
            x1 = read(joinpath(tmpdir_ro, "myroinputfile"), String)
            x2 = read(joinpath(tmpdir_rw, "myrwinputfile"), String)
            @test x1 == "Hello RO world!\n"
            @test x2 == "Hello RW world!\n"
            @test !ispath(joinpath(tmpdir_rw, "myrwoutputfile"))
            @test !isfile(joinpath(tmpdir_rw, "myrwoutputfile"))

            config = DockerConfig(;
                base_image = "julia:latest",
                read_only_maps = Dict(
                    "/home/myuser/workdir_ro" => tmpdir_ro,
                ),
                read_write_maps = Dict(
                    "/home/myuser/workdir_rw" => tmpdir_rw,
                ),
            )

            Utils.make_world_readable_recursively(tmpdir_ro)
            Utils.make_world_writeable_recursively(tmpdir_rw)

            with_container() do container
                code = """
                y1 = read("/home/myuser/workdir_ro/myroinputfile", String)
                y2 = read("/home/myuser/workdir_rw/myrwinputfile", String)
                if y1 != "Hello RO world!\n"
                    throw(ErrorException("File contents did not match"))
                end
                if y2 != "Hello RW world!\n"
                    throw(ErrorException("File contents did not match"))
                end
                if ispath("/home/myuser/workdir_rw/myrwoutputfile")
                    throw(ErrorException("The file exists, but it should not exist"))
                end
                if isfile("/home/myuser/workdir_rw/myrwoutputfile")
                    throw(ErrorException("The file exists, but it should not exist"))
                end
                open("/home/myuser/workdir_rw/myrwoutputfile", "w") do io
                    println(io, "I created and wrote this file inside the container.")
                end
                chmod("/home/myuser/workdir_rw/myrwoutputfile", 0o666)
                mktempdir() do tmpdir
                    open(joinpath(tmpdir, "sometemporaryfile"), "w") do io
                        println(io, "I wrote this file into a temporary directory.")
                    end
                end
                """
                @test success(container, config, `julia -e $(code)`)
            end

            @test ispath(joinpath(tmpdir_rw, "myrwoutputfile"))
            @test isfile(joinpath(tmpdir_rw, "myrwoutputfile"))
            z1 = read(joinpath(tmpdir_rw, "myrwoutputfile"), String)
            @test z1 == "I created and wrote this file inside the container.\n"
        end
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Vector products

function scltrip{T<:Number}(a::AbstractArray{T,1},
                            b::AbstractArray{T,1},
                            c::AbstractArray{T,1})
    dot(a, cross(b, c))
end

function vectrip{T<:Number}(a::AbstractArray{T,1},
                            b::AbstractArray{T,1},
                            c::AbstractArray{T,1})
    cross(a, cross(b, c))
end

a = [3, 4, 5]
b = [4, 3, 5]
c = [-5, -12, -13]

println("Test Vectors:")
println("   a = ", a)
println("   b = ", a)
println("   c = ", a)

println("\nVector Products:")
println("   a dot b = ", dot(a, b))
println("   a cross b = ", cross(a, b))
println("   a dot b cross c = ", scltrip(a, b, c))
println("   a cross b cross c = ", vectrip(a, b, c))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
Create a Julia package/project from
given template.

# Args

- `path`: path of the project to create

# Options

- `--template=<template name>`: template to use, default is the `project` template.
- `--name=<string>`: name of the package, default is the `basename` of `path`.
- `--uuid=<string>`: UUID of the package, default is generated using `uuid1`.
- `--authors=<string>`: authors list, default is empty.
- `--username=<string>`: username of `git`, default will use the system git config.

# Flags

- `-f,--force`: if `path` exists overwrite the content in `path`.
"""
@cast function create(path::String;
        template::String="project",
        name::String=basename(path),
        uuid::String="",
        authors::String="",
        username::String="",
        force::Bool=false
    )
    if !ispath(template)
        template_path = pkgdir(Internal, "templates", "configs", template * ".toml")
        isfile(template_path) || cmd_error("cannot find template $template")
    end
    t = from_toml(Blueprints.Template, template_path)

    # parse CLI inputs
    uuid = isempty(uuid) ? nothing : UUID(uuid)
    authors = isempty(authors) ? String[] : String.(split(authors, ','))
    username = isempty(username) ? nothing : username

    Blueprints.create(t, path; name, uuid, authors, username, force)
    return
end


text:
"""
     computeR1R2(μ)

Compute the non-dimensional distances of each body from the barycenter given the mass
parameter `μ` {NON}.
"""
function computeR1R2(μ)
    R₁ = -μ
    R₂ = 1-μ
    return R₁, R₂
end

"""
     computeR1R2(sys::System)

Compute the non-dimensional distances of each body from the barycenter given CR3BP system
`sys`.
"""
computeR1R2(sys::System) = computeR1R2(sys.μ)

"""
    computeR1R2(p::Array)

Compute the dimensional distances of each body from the barycenter given `p = [μ₁,μ₂,d]`
{km³/s², km³/s², km} which are the gravitational parameters of the first and second primary
bodies and the distance between them.
"""
function computeR1R2(p::Array)
    μ₁,μ₂,d = p
    R₁ = d*μ₂/(μ₁+μ₂)
    R₂ = d*μ₁/(μ₁+μ₂)
    return R₁, R₂
end

"""
    computer1r2(rv, μ)

Compute the non-dimensional position vectors of the particle from each body given the state
vector `rv = [r; v]` {NON; NON} and the mass parameter `μ` {NON}.
"""
function computer1r2(rv,μ)
    x,y,z = rv[1:3]
    r₁ = (x + μ)^2      + y^2 + z^2
    r₂ = (x - 1 + μ)^2  + y^2 + z^2
    return r₁, r₂
end

"""
    computer1r2(rv, sys::System)

Compute the non-dimensional position vectors of the particle from each body given the state
vector `rv = [r; v]` {NON; NON} and the CR3BP system `sys`.
"""
computer1r2(rv,sys::System) = computer1r2(rv,sys.μ)

"""
    computer1r2(rv, p::Array)

Compute the dimensional position vectors of the particle from each body given
`p = [μ₁,μ₂,d]` {km³/s², km³/s², km} which are the gravitational parameters of the first and
second primary bodies and the distance between them.
"""
function computer1r2(rv,p::Array)
    x,y,z = rv[1:3]
    R₁,R₂ = computeR1R2(p)
    r₁ = (x + R₁)^2 + y^2 + z^2
    r₂ = (x - R₂)^2 + y^2 + z^2
    return r₁, r₂
end

"""
    computeL1(μ; tol=1e-15)

Compute position vector of L1 in a normalized CR3BP given the mass parameter `μ` {NON}.
"""
function computeL1(μ;tol=1e-15)
    α = (μ/3 * (1 - μ))^(1/3)
    dα = 1
    count = 0
    while abs(dα) > tol
        α₀ = α
        α = (μ*(1 - α)^2 / (3 - 2*μ - α*(3 - μ - α)))^(1/3)
        dα = α - α₀
        count += 1
        count > 100 ? break : nothing
    end
    L1 = [1 - μ - α; 0; 0]
    return L1
end

"""
    computeL1(μ; tol=1e-15)

Compute position vector of L1 in a normalized CR3BP given the CR3BP system `sys`.
"""
computeL1(sys::System; tol=1e-15) = computeL1(sys.μ, tol=tol)

"""
    computeL1(p::Array;tol=1e-15)

Compute position vector of L1 in a non-normalized CR3BP given `p = [μ₁,μ₂,d]`
{km³/s², km³/s², km} which are the gravitational parameters of the first and second primary
bodies and the distance between them.
"""
function computeL1(p::Array;tol=1e-15)
    μ₁,μ₂,d = p
    μ = μ₂/(μ₁ + μ₂)
    L1 = computeL1(μ;tol=tol)*d
    return L1
end

"""
    computeL2(μ; tol=1e-15)

Compute position vector of L2 in a normalized CR3BP given the mass parameter `μ` {NON}.
"""
function computeL2(μ;tol=1e-15)
    β = (μ/3 * (1 - μ))^(1/3)
    dβ = 1
    count = 0
    while abs(dβ) > tol
        β₀ = β
        β = (μ*(1 + β)^2 / (3 - 2*μ + β*(3 - μ + β)))^(1/3)
        dβ = β - β₀
        count += 1
        count > 100 ? break : nothing
    end
    L2 = [1 - μ + β; 0; 0]
    return L2
end
computeL2(sys::System; tol=1e-15) = computeL2(sys.μ, tol=tol)

"""
    computeL2(p::Array;tol=1e-15)

Compute 3D L2 in a non-normalized CR3BP given p = [μ₁,μ₂,d], which are the
gravitational parameters of the first and second primary bodies [km³/s²]and the
distance between them [km].
"""
function computeL2(p::Array;tol=1e-15)
    μ₁,μ₂,d = p
    μ = μ₂/(μ₁ + μ₂)
    L2 = computeL2(μ;tol=tol)*d
    return L2
end

"""
    computeL3(μ;tol=1e-15)

Compute position vector of L3 in a normalized CR3BP given the mass parameter `μ` {NON}.
"""
function computeL3(μ;tol=1e-15)
    γ = -7*μ/12 + 1
    dγ = 1
    count = 0
    while abs(dγ) > tol
        γ₀ = γ
        γ = ((1 - μ)*(1 + γ)^2 / (1 + 2*μ + γ*(2 + μ + γ)))^(1/3)
        dγ = γ - γ₀
        count += 1
        count > 100 ? break : nothing
    end
    L3 = [-μ - γ; 0; 0]
    return L3
end

computeL3(sys::System; tol=1e-15) = computeL3(sys.μ, tol=tol)

"""
    computeL3(p::Array;tol=1e-15)

Compute 3D L3 in a non-normalized CR3BP given p = [μ₁,μ₂,d], which are the
gravitational parameters of the first and second primary bodies [km³/s²]and the
distance between them [km].
"""
function computeL3(p::Array;tol=1e-15)
    μ₁,μ₂,d = p
    μ = μ₂/(μ₁ + μ₂)
    L3 = computeL3(μ;tol=tol)*d
    return L3
end

"""
    computeL4(μ; tol=1e-15)

Compute position vector of L4 in a normalized CR3BP given the mass parameter `μ` {NON}.
"""
computeL4(μ;tol=1e-15) = [0.5-μ; √3/2; 0]
computeL4(sys::System; tol=1e-15) = computeL4(sys.μ, tol=tol)

"""
    computeL4(p::Array;tol=1e-15)

Compute 3D L4 in a non-normalized CR3BP given p = [μ₁,μ₂,d], which are the
gravitational parameters of the first and second primary bodies [km³/s²]and the
distance between them [km].
"""
function computeL4(p::Array;tol=1e-15)
    μ₁,μ₂,d = p
    μ = μ₂/(μ₁ + μ₂)
    L4 = computeL4(μ;tol=tol)*d
    return L4
end

"""
    computeL5(μ;tol=1e-15)

Compute position vector of L5 in a normalized CR3BP given the mass parameter `μ` {NON}.
"""
computeL5(μ;tol=1e-15) = [0.5-μ; -√3/2; 0]
computeL5(sys::System; tol=1e-15) = computeL5(sys.μ, tol=tol)

"""
    computeL5(p::Array;tol=1e-15)

Compute 3D L5 in a non-normalized CR3BP given p = [μ₁,μ₂,d], which are the
gravitational parameters of the first and second primary bodies [km³/s²]and the
distance between them [km].
"""
function computeL5(p::Array;tol=1e-15)
    μ₁,μ₂,d = p
    μ = μ₂/(μ₁ + μ₂)
    L5 = computeL5(μ;tol=tol)*d
    return L5
end

"""
    computeLpts(μ;tol=1e-15)

Compute 3D Lagrange points in a normalized CR3BP given μ, the system mass parameter.
"""
function computeLpts(μ; tol=1e-15)
    Lpts = [
        computeL1(μ,tol=tol),
        computeL2(μ,tol=tol),
        computeL3(μ,tol=tol),
        computeL4(μ,tol=tol),
        computeL5(μ,tol=tol),
        ]
    return Lpts
end

"""
    computeLpts(sys::System;tol=1e-15)

Compute 3D Lagrange points in a non-normalized CR3BP given system S.
"""
computeLpts(sys::System; tol=1e-15) = computeLpts(sys.μ, tol=tol)

"""
    computeLpts(p::Array;tol=1e-15)

Compute 3D Lagrange points in a non-normalized CR3BP given p = [μ₁,μ₂,d], which are the
gravitational parameters of the first and second primary bodies [km³/s²]and the
distance between them [km].
"""
function computeLpts(p::Array;tol=1e-15)
    μ₁,μ₂,d = p
    μ = μ₂/(μ₁ + μ₂)
    Lpts = computeLpts(μ, tol=tol)*d
    return Lpts
end

"""
    computeUeff(rv,μ)

Compute effective potential given normalized state `rv` {NON} and mass parameter `μ` {NON}.
"""
function computeUeff(rv,μ)
    x,y,z = rv[1:3]
    r₁,r₂ = computer1r2(rv,μ)
    Ueff = -(x^2 + y^2)/2 - (1-μ)/r₁ - μ/r₂;
    return Ueff
end

"""
    computeUeff(rv,sys::System)

Compute effective potential given state `rv = [r; v]`` {km; km/s} and CR3BP system `sys`.
"""
computeUeff(rv,sys::System) = computeUeff(rv,sys.μ)

"""
    computeUeff(rv,p::Array)

Compute effective potential given state `rv = [r; v]`` {km; km/s} and `p = [μ₁,μ₂,d]`
{km³/s²; km³/s²; km}, which are the gravitational parameters of the first and second primary
bodies and the distance between them.
"""
function computeUeff(rv,p::Array)
    x,y,z = rv[1:3]
    R₁,R₂ = computeR1R2(p)
    r₁,r₂ = computer1r2(rv,p)
    μ₁,μ₂,d = p
    ωₛ = sqrt((μ₁ + μ₂)/d^3)
    Ueff = -(x^2 + y^2)*ωₛ^2/2 - μ₁/r₁ - μ₂/r₂;
    return Ueff
end

"""
    computeC(rv,μ)

Compute Jacobi constant given normalized state rv {NON} and mass parameter {NON}
"""
function computeC(rv,μ)
    v = norm(rv[4:6])
    Ueff = computeUeff(rv,μ)
    C = -2*Ueff - v^2
    return C
end

"""
    computeC(rv,sys::System)

Compute Jacobi constant given normalized state rv {NON} and system
"""
computeC(rv,sys::System) = computeC(rv, sys.μ)

"""
    computeC(rv,p::Array)

Compute Jacobi constant given state rv = [r; v] {km; km/s} and p = [μ₁,μ₂,d],
which are the gravitational parameters of the first and second primary bodies
[km³/s²] and the distance between them [km].
"""
function computeC(rv,p::Array)
    μ₁,μ₂,d = p
    v = norm(rv[4:6])
    Ueff = computeUeff(rv,p)
    C = -2*Ueff - v^2
    return C
end

"""
    computeT(rv,sys::System) UNFINISHED (I don't like this name, "tisserand" is better)

Compute Tisserand parameter given orbital elements
"""
function computeT(a,e,i; aⱼ=7.783561990635208e8, ang_unit::Symbol=:deg)
    if ang_unit == :deg
    elseif ang_unit == :rad
        i = rad2deg(i)
    else
        error("ang_unit should be :rad or :deg")
    end

    T = aⱼ/a + s*sqrt(a/aⱼ*(1 - e^2)*cosd(i))
    return T
end

"""
    stability_index(Φ) UNFINISHED

Compute the stability index for a trajectory given its state transition matrix Φ
"""
function stability_index(Φ)
    return 1
end

"""
    rotx(θ)
"""
function rotx(θ)
    R = [1      0       0;
         0 cos(θ) -sin(θ);
         0 sin(θ)  cos(θ)]
    return R
end

"""
    rotxd(θ)
"""
function rotxd(θ)
    R = [1      0       0;
         0 cosd(θ) -sind(θ);
         0 sind(θ)  cosd(θ)]
    return R
end

"""
    roty(θ)
"""
function roty(θ)
    R = [cos(θ)  0  sin(θ);
              0  1       0;
        -sin(θ)  0  cos(θ)]
    return R
end

"""
    rotyd(θ)
"""
function rotyd(θ)
    R = [cosd(θ)  0  sind(θ);
               0  1        0;
        -sind(θ)  0  cosd(θ)]
    return R
end

"""
    rotz(θ)
"""
function rotz(θ)
    R = [cos(θ) -sin(θ) 0;
         sin(θ)  cos(θ) 0;
              0       0 1]
    return R
end

"""
    rotzd(θ)
"""
function rotzd(θ)
    R = [cosd(θ) -sind(θ) 0;
         sind(θ)  cosd(θ) 0;
               0        0 1]
    return R
end

"""
    rotlatlon(ϕ, λ; ang_unit=:deg)
"""
function rotlatlon(ϕ, λ; ang_unit=:deg)
    if ang_unit == :rad
        Φ = rotx(π/2 - ϕ)
        Λ = rotz(π/2 + λ)
    elseif ang_unit == :deg
        Φ = rotxd(90 - ϕ)
        Λ = rotzd(90 + λ)
    else
        error("ang_unit should be :deg or :rad")
    end
    R = Λ*Φ
    return R
end

"""
    date2mjd(ut1_date)
"""
function date2mjd(ut1_date)
    if size(ut1_date) == (3,)
        Y = ut1_date[1]
        M = ut1_date[2]
        D = ut1_date[3]
    elseif size(ut1_date) == (6,)
        Y = ut1_date[1]
        M = ut1_date[2]
        D = ut1_date[3] + ut1_date[4]/24 + ut1_date[5]/24/60 + ut1_date[6]/24/60/60
    else
        error("ut1_date should be [Y,M,D] or [Y,M,D,h,m,s]")
    end

    if M <= 2
        y = Y-1;
        m = M+12;
    else
        y = Y;
        m = M;
    end

    if Y < 1582
        B = -2 + ((y+4716/4)/4) - 1179;
    elseif Y > 1582
        B = y/400 - y/100 + y/4;
    else
        if M < 10
            B = -2 + ((y+4716/4)/4) - 1179;
        elseif M > 10
            B = y/400 - y/100 + y/4;
        else
            if D <= 4
                B = -2 + ((y+4716/4)/4) - 1179;
            elseif D >= 10
                B = y/400 - y/100 + y/4;
            else
                B = NaN;
            end
        end
    end

    ut1_mjd = 365*y - 679004 + floor(B) + floor(30.6001*(m+1)) + D;

end

"""
    mjd2gmst(ut1_mjd; ang_unit::Symbol=:rad)
"""
function mjd2gmst(ut1_mjd; ang_unit::Symbol=:rad)
    θ₀ = 280.4606 # [deg] Greenwich Mean Sidereal time at J2000 Epoch
    ωₑ = 360.9856473 # [deg/day] rotating rate of the Earth
    d = ut1_mjd - 51544.5 # Normalize by epoch (Jan. 1, 2000 12:00h)

    if ang_unit == :rad
        θ₀ = deg2rad(θ₀)
        ωₑ = deg2rad(ωₑ)
        θ = θ₀ + ωₑ*d
        θ = wrapto2pi(θ)
    elseif ang_unit == :deg
        θ = θ₀ + ωₑ*d
        θ = wrapto360(θ)
    else
        error("ang_unit should be :deg or :rad")
    end
    return θ
end

"""
    wrapto360(θ)
"""
function wrapto360(θ)
    while θ < 0;    θ += 360;   end
    while θ > 360;  θ -= 360;   end
    return θ
end

"""
    wrapto180(θ)
"""
function wrapto180(θ)
    while θ < -180;  θ += 360;   end
    while θ > 180; θ -= 360;   end
    return θ
end

"""
    wrapto2pi(θ)
"""
function wrapto2pi(θ)
    while θ < 0;    θ += 2π;    end
    while θ > 2π;   θ -= 2π;    end
    return θ
end

"""
    wraptopi(θ)
"""
function wraptopi(θ)
    while θ < -π;   θ += 2π;    end
    while θ > π;    θ -= 2π;    end
    return θ
end

"""
    date2str(date)
"""
function date2str(date)
    Y = Int(date[1])
    M = Int(date[2])
    D = Int(date[3])
    h = Int(date[4])
    m = Int(date[5])
    s = date[6]
    return string(Y, "-", M, "-", D, " ", h, ":", m, ":", s)
end

"""
    deserno_sphere(N_desired)

    Adapted from Lucas Bury's code implementing the algorithm in Deserno, M. (2004) How to
    Generate Equidistributed points on the Surface of a Sphere
"""
function deserno_sphere(N_desired)

area = 4*pi/N_desired
distance = sqrt(area)

M_theta = round(pi/distance);

d_theta = pi/M_theta

d_phi = area/d_theta;

N_new = 0;
xs = zeros(0)
ys = zeros(0)
zs = zeros(0)
for m in 0:(M_theta-1)

    theta = pi*(m+0.5)/M_theta;
    M_phi = round(2*pi*sin(theta)/d_phi); # not exact

    for n in 0:(M_phi-1)
        Phi = 2*pi*n/M_phi;

        N_new = N_new + 1;

        append!(xs, sin(theta)*cos(Phi))
        append!(ys, sin(theta)*sin(Phi))
        append!(zs, cos(theta))
    end
end

xyz_Sphere = copy(transpose(hcat(xs, ys, zs)))

return xyz_Sphere, N_new
end

"""
deserno_hemisphere(N_desired, normal)
"""
function deserno_hemisphere(N_desired, normal)
# -------------------------------------------------
### Obtaining spherical point cloud
# -------------------------------------------------
xyz_Sphere, N_sphere = deserno_sphere(N_desired*2);

# -------------------------------------------------
### Grab the desired hemisphere
# -------------------------------------------------
### Indices points with a 'z' value >= 0
indices = xyz_Sphere[3,:] .>= 0

### Use logical indexing to grab all points with a 'z' >= 0
xyz_hem = xyz_Sphere[:,indices]
N_new = size(xyz_hem, 2)

# -------------------------------------------------
### If new center is in the -z axis, just flip the signs on the existing hemisphere
# -------------------------------------------------
### In getSphere, the pattern is centered around the 'z' axis
oldNormal = [0, 0, 1]

### Turn input vector to unit vector
normal = normal ./ norm(normal)

### In case the new center is in the opposite direction of the old center, can
### save computational time with this (but also the other algorithm breaks)
if normal == -oldNormal
    xyz_unitHemisphere = zeros(size(xyz_hem))
    for kk in 1:N_new
        xyz_unitHemisphere[:,kk] = [xyz_hem[1,kk], xyz_hem[2,kk],-xyz_hem[3,kk]]
    end
    return xyz_unitHemisphere, N_new
end

# -------------------------------------------------
### Compute the rotation matrix to properly align the new hemisphere
# -------------------------------------------------
### Algorithm to compute a rotation matrix between two vectors
v = cross(oldNormal, normal)
skewSymmetricMat = [0 -v[3] v[2]; v[3] 0 -v[1]; -v[2] v[1] 0]
R_old2new = [1 0 0; 0 1 0; 0 0 1] + skewSymmetricMat + skewSymmetricMat*skewSymmetricMat*(1 / (1 + dot(oldNormal,normal)))

### Rotate old hemisphere to new position, centered about normal
xyz_unitHemisphere = zeros(size(xyz_hem))
for kk in 1:N_new
    xyz_unitHemisphere[:,kk] = R_old2new * xyz_hem[:,kk]
end

return xyz_unitHemisphere, N_new
end
"""
    Generates a ring of N vectors with origin c offset by α degrees about central vector r
"""
function spherical_ring(c::Array, r::Array, α; N=100)
    ϕ = 90 - α
    λ = LinRange(360/N,360,N)
    r_ring = [latlon2cart(ϕ,λ[i],norm(r),ang_unit=:deg) for i = 1:N]
    n = [1,0,0]
    r = r/norm(r)
    if r'*n > 0.9
        n[1:3] = [0,1,0]
    end
    p = cross(r,n)
    p = p/norm(p)
    q = cross(r,p)
    q = q/norm(q)
    A = [p q r]
    r_ring = [A*r_ring[i] + c for i = 1:N]
    return r_ring
end

"""
Calculate the ratio of Lagrange point distance from closest primary to distance between two primaries
"""
function gammaL(sys::System, Lpt::Int)
    μ = sys.μ

    # poly1 = [1, -1*(3-μ), (3-2*μ),     -μ,      2*μ,     -μ];
    # poly2 = [1,    (3-μ), (3-2*μ),     -μ,     -2*μ,     -μ];
    # poly3 = [1,    (2+μ), (1+2*μ), -(1-μ), -2*(1-μ), -(1-μ)];

    poly1 = [    -μ,      2*μ,     -μ, (3-2*μ), -1*(3-μ), 1];
    poly2 = [    -μ,     -2*μ,     -μ, (3-2*μ),    (3-μ), 1];
    poly3 = [-(1-μ), -2*(1-μ), -(1-μ), (1+2*μ),    (2+μ), 1];

    rt1 = roots(poly1)
    rt2 = roots(poly2)
    rt3 = roots(poly3)

    Γ = zeros(3)
    for i=1:5
            if isreal(rt1[i]) Γ[1]=rt1[i]; end
            if isreal(rt2[i]) Γ[2]=rt2[i]; end
            if isreal(rt3[i]) Γ[3]=rt3[i]; end
    end
    γ = Γ[Lpt];
end


text:
# This file is part of JuliaFMI.
# Licensed under MIT: https://github.com/AnHeuermann/JuliaFMI/blob/master/LICENSE.txt

# This file contains type definitions and constructores for those types.

"""
Declaration of FMI2 Types
"""

@enum fmuType begin
    modelExchange
    coSimulation
end

@enum NamingConvention begin
    flat
    structured
end
function NamingConvention(in::String)
    if in=="flat"
        return flat
    elseif in == "structured"
        return structured
    else
        error("Can not convert String \"$in\" to NamingConvention.")
    end
end

@enum ModelState begin
    modelUninstantiated
    modelInstantiated
    modelInitializationMode
    modelContinuousTimeMode
    modelEventMode
    modelTerminated
    modelError
    modelFatal
end

# FMI2 Errors
struct FMI2Warning <: Exception
    msg::String
end
FMI2Warning() = FMI2Warning("")
Base.showerror(io::IO, e::FMI2Warning) = print(io, "fmi2Warning", isempty(e.msg) ? "" : ": ", e.msg)

struct FMI2Discard <: Exception
    msg::String
end
FMI2Discard() = FMI2Discard("")
Base.showerror(io::IO, e::FMI2Discard) = print(io, "fmi2Discard", isempty(e.msg) ? "" : ": ", e.msg)

struct FMI2Error <: Exception
    msg::String
end
FMI2Error() = FMI2Error("")
Base.showerror(io::IO, e::FMI2Error) = print(io, "fmi2Error", isempty(e.msg) ? "" : ": ", e.msg)

struct FMI2Fatal <: Exception
    msg::String
end
FMI2Fatal() = FMI2Fatal("")
Base.showerror(io::IO, e::FMI2Fatal) = print(io, "fmi2Fatal", isempty(e.msg) ? "" : ": ", e.msg)

struct FMI2Pending <: Exception
    msg::String
end
FMI2Pending() = FMI2Pending("")
Base.showerror(io::IO, e::FMI2Pending) = print(io, "fmi2Pending", isempty(e.msg) ? "" : ": ", e.msg)

function fmiError(fmi2Status::Union{Unsigned, Integer}, message::String="")
    if fmi2Status == 1
        return FMI2Warning(message)
    elseif fmi2Status == 2
        return FMI2Discard(message)
    elseif fmi2Status == 3
        return FMI2Error(message)
    elseif fmi2Status ==  4
        return FMI2Fatal(message)
    elseif fmi2Status == 5
        return FMI2Pending(message)
    else
        return FMI2Error("Unknown error code")
    end
end

# Pointers to functions provided by the environment to be used by the FMU
struct CallbackFunctions
    callbackLogger::Ptr{Nothing}
    allocateMemory::Ptr{Nothing}
    freeMemory::Ptr{Nothing}
    stepFinished::Ptr{Nothing}

    componentEnvironmendt::Ptr{Nothing}
end

mutable struct fmi2ComponentEnvironment
    logFile::String     # if not empty location of file to write logger messages
                        # defaults to stderr ??
    numWarnings::Int
    numErrors::Int
    numFatals::Int
end

mutable struct RealVariable
    value::Float64
    valueReference::UInt
    name::String

    # attributes
    min::Float64
    max::Float64

    # Inner constructors
    RealVariable()=new()
    function RealVariable(value, valueReference, name)
        new(value, valueReference, name, -Inf64, Inf64)
    end
end

mutable struct IntVariable
    value::Int64
    valueReference::UInt
    name::String

    # attributes
    min::Int64
    max::Int64

    # Inner constructors
    IntVariable()=new()
    function IntVariable(value, valueReference, name)
        new(value, valueReference, name, -(2^63 - 1), 2^63 - 1)
    end
end

mutable struct BoolVariable
    value::Bool
    valueReference::UInt
    name::String

    # Inner Constructors
    BoolVariable() = new()
    function BoolVariable(value, valueReference, name)
        new(value, valueReference, name)
    end
end

mutable struct StringVariable
    value::String
    valueReference::UInt
    name::String

    # Inner constructors
    StringVariable() = new()
    function StringVariable(value, valueReference, name)
        new(value, valueReference, name)
    end
end

mutable struct EnumerationVariable
    value::UInt     # Holds a number representig an enumeration
    valueReference::UInt
    name::String

    # Inner constructors
    EnumerationVariable() = new()
    function EnumerationVariable(value, valueReference, name)
        new(value, valueReference, name)
    end
end

mutable struct ModelVariables
    reals::Array{RealVariable,1}
    ints::Array{IntVariable,1}
    bools::Array{BoolVariable,1}
    strings::Array{StringVariable,1}
    enumerations::Array{EnumerationVariable,1}

    function ModelVariables(n_reals, n_ints, n_bools, n_strings, n_enumerations)

        reals = Array{RealVariable}(undef, n_reals)
        ints = Array{IntVariable}(undef, n_ints)
        bools = Array{BoolVariable}(undef, n_bools)
        strings = Array{StringVariable}(undef, n_strings)
        enumerations = Array{EnumerationVariable}(undef, n_enumerations)

        new(reals, ints, bools, strings, enumerations)
    end
end

mutable struct SimulationData
    time::AbstractFloat
    lastStepTime::AbstractFloat
    modelVariables::ModelVariables
    eventIndicators::Array{Float64}

    SimulationData()=new()
    function SimulationData(n_reals, n_ints, n_bools, n_strings, n_enumerations,
        numberOfEventIndicators)

        modelVariables=ModelVariables(n_reals, n_ints, n_bools, n_strings,
            n_enumerations)
        eventIndicators = Array{Float64}(undef, numberOfEventIndicators)
        new(0, 0, modelVariables, eventIndicators)
    end
end

mutable struct ModelData
    numberOfStates::Int
    numberOfDerivatives::Int    # Is always numberOfStates

    numberOfReals::Int
    numberOfInts::Int
    numberOfBools::Int
    numberOfStrings::Int
    numberOfEnumerations::Int
    numberOfExterns::Int

    numberOfEventIndicators::Int

    ModelData() = new(0,0,0,0,0,0,0,0)
end

mutable struct ExperimentData
    startTime::AbstractFloat
    stopTime::AbstractFloat
    tolerance::AbstractFloat
    stepSize::AbstractFloat

    # Inner constructors
    ExperimentData() = new(0, 1, 1e-6, 1e-6/4)
    function ExperimentData(startTime, stopTime, tolerance, stepSize)
        if startTime >= stopTime
            error("ExperimentData not valid: startTime=$startTime >= stopTime=$stopTime")
        elseif tolerance <= 0
            error("ExperimentData not valid: tolerance=$tolerance not greater than zero")
        elseif stepSize <= 0
            error("ExperimentData not valid: stepSize=$stepSize not greater than zero")
        end
        new(startTime, stopTime, tolerance, stepSize)
    end
end

struct RealAttributes
    quantity::String
    unit::String            # TODO: make types for Units and functions
    displayUnit::String     #       for unit conversion
    relativeQuantity::Bool
    min::Real
    max::Real
    nominal::Real
    unbound::Bool

    # Inner constructor
    RealAttributes() = new()

    function RealAttributes(quantity, unit, displayUnit, relativeQuantity,
        min, max, nominal, unbound)

        if nominal <= 0
            error("Nominal > 0.0 required")
        elseif min > max
            error("Minimum is greater than maximum")
        end

        new(quantity, unit, displayUnit, relativeQuantity, min, max, nominal,
            unbound)
    end

    function RealAttributes(quantity, unit, displayUnit, min, max, nominal)

        if nominal <= 0
            error("Nominal > 0.0 required")
        elseif min > max
            error("Minimum is greater than maximum")
        end

        new(quantity, unit, displayUnit, false, min, max, nominal, false)
    end
end

struct IntegerAttributes
    quantity::String
    min::Int
    max::Int

    # Inner constructors
    IntegerAttributes()=new()
    IntegerAttributes(quantity, min, max)=new(quantity, min, max)
end

"""
Item of an enumeartion attribute
"""
struct EnumerationItemAttribute
    name::String
    "Unique number in the same enumeration"
    value::Int
    description::String

    # Constructor
    EnumerationItemAttribute() = new()
    EnumerationItemAttribute(name, value, description) = new(name, value, description)
end

"""
Attributes for simpleType enumeration
"""
struct EnumerationAttributes
    "Physical quantity of the variable"
    quantity::String
    "Items of an enumeration"
    items::Array{EnumerationItemAttribute}

    # Constructor
    EnumerationAttributes() = new()
    EnumerationAttributes(quantity, items) = new(quantity, items)
end

struct RealProperties
    declaredType::String
    variableAttributes::RealAttributes
    start::Float64
    derivative::UInt
    reinit::Bool

    # Inner constructors
    RealProperties() = new()
    RealProperties(declaredType, variableAttributes, start, derivative, reinit) = new(declaredType, variableAttributes, start, derivative, reinit)
    function RealProperties(declaredType, variableAttributes, start, derivative, reinit_in::String)
        if reinit_in == "true"
            reinit = true
        elseif reinit_in == "false"
            reinit = false
        else
            error("Could not parste input reinit=$reinit_in to Bool.")
        end

        return new(declaredType, variableAttributes, start, derivative, reinit)
    end
end

struct IntegerProperties
    declaredType::String
    variableAttributes::IntegerAttributes
    start::Int

    # Inner constructors
    IntegerProperies() = new()
    IntegerProperties(declaredType, variableAttributes, start) = new(declaredType, variableAttributes, start)
end

struct BooleanProperties
    declaredType::String
    start::Bool

    BooleanProperties() = new()
    BooleanProperties(declaredType, start) = new(declaredType, start)
end

struct StringProperties
    declaredType::String
    start::String

    StringProperties() = new()
    StringProperties(declaredType, start) = new(declaredType, start)
end

struct EnumerationProperties
    declaredType::String
    quantity::String
    min::Int
    max::Int
    start::Int

    EnumerationProperties() = new()
    EnumerationProperties(declaredType, quantity, min, max, start) = new(declaredType, quantity, min, max, start)
end

struct ScalarVariable
    name::String
    valueReference::Unsigned

    # Optional
    description::String
    causality::String           # TODO: Change to enumeration??
    variability::String         # TODO: Change to enumeration??
    initial::String             # TODO: Change to enumeration??
    canHandleMultipleSetPerTimelnstant::Bool

    # Type specific properties of ScalarVariable
    typeSpecificProperties::Union{RealProperties, IntegerProperties, BooleanProperties, StringProperties, EnumerationProperties}

    # Inner constructors
    function ScalarVariable(name, valueReference, description, typeSpecificProperties)
        if isempty(strip(name))
            error("ScalarVariable $name not valid: name can't be empty or only whitespace")
        elseif valueReference < 0
            error("ScalarVariable $name not valid: valueReference=$valueReference not unsigned")
        end
        new(name, Unsigned(valueReference), description, "local", "continous", "", false, typeSpecificProperties)
    end

    function ScalarVariable(name, valueReference, description, causality,
        variability, initial, canHandleMultipleSetPerTimelnstant, typeSpecificProperties)

        # Check name
        if isempty(strip(name))
            error("ScalarVariable $name not valid: name can't be empty or only whitespace")
        end

        # Check valueReference
        if valueReference < 0
            error("ScalarVariable $name not valid: valueReference=$valueReference not unsigned")
        end

        #check if causality and variability are correct
        if isempty(variability)
            variability = "continuous"
        elseif !in(variability, ["constant", "fixed", "tunable", "discrete", "continuous"])
            error("ScalarVariable $name not valid: variability has to be one of \"constant\", \"fixed\", \"tunable\", \"discrete\" or \"continuous\" but is \"$variability\"")
        end

        if isempty(causality)
            causality = "local"
        elseif !in(causality,["parameter", "calculatedParameter","input", "output", "local", "independent"])
          error("ScalarVariable $name not valid: causality has to be one of \"parameter\", \"calculatedParameter\", \"input\", \"output\", \"local\", \"independent\" but is \"$causality\"")
        elseif causality == "parameter"
            if (variability!="fixed" && variability!="tunable" && variability!="constant")
                error("ScalarVariable $name not valid: causality is \"parameter\", so variability has to be \"fixed\" or \"tunable\" but is \"$variability\"")
            end
            if isempty(initial)
                initial = "exact"
            elseif initial!="exact"
                error("ScalarVariable $name not valid: causality is \"parameter\", so initial has to be \"exact\" or empty but is \"$initial\"")
            end
        elseif causality=="calculatedParameter"
            if (variability!="fixed" && variability!="tunable")
                error("ScalarVariable $name not valid: causality is \"calculatedParameter\", so variability has to be \"fixed\" or \"tunable\" but is \"$variability\"")
            end
            if isempty(initial)
                initial = "calculated"
            elseif !in(initial, ["approx", "calculated"])
                error("ScalarVariable $name not valid: causality is \"calculatedParameter\", so initial has to be \"approx\", \"calculated\" or empty but is \"$initial\"")
            end
        elseif causality=="input"
            if !isempty(initial)
                error("ScalarVariable $name not valid: causality is \"input\", so initial has to be empty but is \"$initial\"")
            end
        elseif causality=="independent"
            if variability!="continuous"
                error("ScalarVariable $name not valid: causality is \"independent\", so variability has to be \"continuous\" but is \"$causality\"")
            end
        end

        new(name, Unsigned(valueReference), description, causality, variability, initial, canHandleMultipleSetPerTimelnstant, typeSpecificProperties)
    end
end

"""
Definition of a simple type according to fmi2SimpleType
"""
struct SimpleType
    "Unique name with respect to all other element of TypeDefinitions in ModelDescription"
    name::String
    description::String

    # ToDo: Add Bool and String to this, maybe as empty Attributes?
    attributes::Union{RealAttributes, IntegerAttributes, EnumerationAttributes}

    # Constructor
    SimpleType() = new()
    SimpleType(name, description, attributes)=new(name, description, attributes)
end

struct LogCategory
    name::String
    description::String

    LogCategory(name) = new(name, "")
    LogCategory(name, description) = new(name, description)
end

"""
Containing all informations from modelDescription.xml
"""
mutable struct ModelDescription
    # FMI model description
    fmiVersion::String
    modelName::String
    guid::String
    description::String
    author::String
    version::String
    copyright::String
    license::String
    generationTool::String
    generationDateAndTime::String
    variableNamingConvention::NamingConvention
    numberOfEventIndicators::Int

    # Model exchange
    isModelExchange::Bool
    # Co-Simulation
    isCoSimulation::Bool
    modelIdentifier::String

    # Unit definitions
    # Type definitions
    typeDefinitions::Union{Array{SimpleType}, Nothing}

    logCategories::Union{Array{LogCategory}, Nothing}

    # Default experiment
    defaultExperiment::ExperimentData

    # Vendor annotations

    # Model variables
    modelVariables::Array{ScalarVariable,1}

    # Model structure
    modelStructure

    # Constructor for uninitialized struct
    function ModelDescription()
        md = new()
        md.isModelExchange = false
        md.isCoSimulation = false
        return md
    end
end

mutable struct EventInfo
    newDiscreteStatesNeeded::Bool
    terminateSimulation::Bool
    nominalsOfContinuousStatesChanged::Bool
    valuesOfContinuousStatesChanged::Bool
    nextEventTimeDefined::Bool
    nextEventTime::Float64

    EventInfo() = new(true, true, true, true, true, -1.0)
end

"""
Functional Mockupt Unit (FMU) struct.
"""
mutable struct FMU
    modelName::String
    instanceName::String
    FMUPath::String                     # TODO: find better type for paths
    fmuResourceLocation::String         # is URI
    fmuGUID::String

    modelDescription::ModelDescription

    modelData::ModelData

    simulationData::SimulationData

    experimentData::ExperimentData

    eventInfo::EventInfo

    modelState::ModelState

    # CSV and log file
    csvFile::IOStream
    logFile::IOStream

    # Other stuff
    libHandle::Ptr{Nothing}
    libLoggerHandle::Ptr{Nothing}
    tmpFolder::String
    fmuType::fmuType
    fmiCallbackFunctions::CallbackFunctions
    fmi2Component::Ptr{Nothing}

    # Constructor
    FMU() = new()
end

max_stars_repo_name: JKRT/JuliaFMI

max_issues_repo_name: JKRT/JuliaFMI

max_forks_repo_name: JKRT/JuliaFMI


text:
using Parameters
using ForwardDiff
using NLsolve
using Devectorize

using Voigt
using Voigt.Unicode

using ConstLab

create_component_macros("mises_visc", (1, 1, 1))

const II = veye(6) ⊗ veye(6)
const Idev = eye(6,6) - 1/3 * II

@with_kw immutable ViscoPlastMS
    ₙεₚ::Vector{Float64} = zeros(6)
    ₙσ::Vector{Float64} = zeros(6)
    ₙκ::Float64 = 0.0
    ₙμ::Float64 = 0.0
    loading::Bool = false
end

@with_kw immutable ViscoPlastMP
    E::Float64
    ν::Float64
    Ee::Matrix{Float64} =   (G = E / 2(1 + ν);
                            K = E / 3(1 - 2ν);
                            Ee = 2 * G * Idev+ K * II;
                            Ee[4:6, 4:6] /= 2.0;
                            Ee)
    σy::Float64
    H::Float64
    n::Float64
    σc::Float64
    t★::Float64
end

function stress(ε, ∆t, matpar::ViscoPlastMP, matstat::ViscoPlastMS)
    @unpack matpar: E, Ee, ν, σy, H, n, σc, t★
    @unpack matstat: ₙεₚ, ₙσ, ₙκ, ₙμ

    σₜᵣ = Ee * (ε - ₙεₚ)
    σdevₜᵣ = dev(σₜᵣ)
    sₜᵣ = sqrt(3/2) * vnorm(σdevₜᵣ)
    Φ = sₜᵣ - σy - ₙκ

    if Φ < 0
        return σₜᵣ, ViscoPlastMS(ₙεₚ, σₜᵣ, ₙκ, ₙμ, false)
    else

        R!(x, fx) = compute_residual!(fx, x, sₜᵣ, ∆t, matpar, matstat)
        x0 = [sₜᵣ, ₙκ, ₙμ]

        res = nlsolve(R!, x0; iterations = 20, ftol=1e-5, method=:newton, autodiff = true)
        if !converged(res)
            error("Material did not converge")
        end
        X = res.zero::Vector{Float64}
        ConstLab.@unpack_mises_visc (s, κ, μ) =  X

        G = E / 2(1 + ν)
        @devec σ_dev = 1 ./ (1 + 3 .* G .* μ ./ s) .* σdevₜᵣ
        @devec σ = σₜᵣ - σdevₜᵣ + σ_dev
        @devec εₚ = ₙεₚ + 3 ./ (2 .* s) .* μ .* σ_dev

        return σ, ViscoPlastMS(εₚ, σ, κ, μ, true)
    end
end

function compute_residual!(fx, x, sₜᵣ, ∆t, matpar, matstat)
    @unpack matpar: E, ν, σy, H, n, σc, t★
    @unpack matstat: ₙεₚ, ₙσ, ₙκ, ₙμ

    ConstLab.@unpack_mises_visc (s, κ, μ) = x

    Φ = max(0.0, s - σy - κ)
    G = E / 2(1 + ν)
    R_s = s - sₜᵣ + 3G*μ
    R_κ = κ - ₙκ - H * μ
    Γ′ = (Φ / σc)^n
    Γ′′ =  n / σc * (Φ / σc)^(n-1)
    R_μ = Γ′ - t★ / ∆t * μ

    R_s′ = R_s/3G
    R_κ′ = R_κ/H
    R_μ′ = R_μ/Γ′′
    ConstLab.@pack_mises_visc fx = (R_s′, R_κ′, R_μ′)
    return
end

function compute_jacobian(x, ∆t, matpar, matstat)
    @unpack matpar: E, ν, σy, H, n, σc, t★

    ConstLab.@unpack_mises_visc (s, κ, μ) = x

    G = E / 2(1 + ν)
    Φ = max(zero(eltype(x)), s - σy - κ)
    Γ′′ =  n / σc * (Φ / σc)^(n-1)

    J = zeros(3,3)
    J[1,1] = 1/3G
    J[1,3] = 1
    J[2,2] = 1/H
    J[2,3] = -1
    J[3,1] = 1
    J[3,2] = -1
    J[3,3] = -1/Γ′′ *t★/∆t
    return J
end

function ATS(ε, ∆t, matpar::ViscoPlastMP, matstat::ViscoPlastMS)
    @unpack matpar: E, Ee, ν, σy, H, n, σc, t★
    @unpack matstat: ₙεₚ, ₙσ, ₙκ, ₙμ, loading

    if !loading
        return Ee
    else
        # ATS
        σ, μ, κ = ₙσ, ₙμ, ₙκ
        σdev = dev(σ)
        s = sqrt(3/2) * vnorm(σdev)
        Φ = s - σy - κ
        G = E / 2(1 + ν)

        σdevdev = σdev ⊗ σdev
        Γ′′ =  n / σc * (Φ / σc)^(n-1)
        h_at = 3 * G + H + 1 / Γ′′ * t★ / ∆t
        b = (3 * G * μ / s) / (1 + 3 * G * μ / s)
        @devec Q = Idev - (3./(2.*s.^2) .* (σdevdev))
        @devec ATS = Ee + (- 2 .* G .* b .* Q - 9 .* G.^2 ./ (h_at .* s.^2) .* (σdevdev))
        @devec ATS[4:6, 4:6] ./= 2.0
        return ATS
    end
end


text:
function getCoefficientsLobIII2(T=Float64)
    a = Array{Dec128}( [[0//1  0//1]
                        [1//1  0//1]] )
    b = Array{Dec128}( [1//2, 1//2] )
    c = Array{Dec128}( [0//1, 1//1] )

    CoefficientsRK(T, :LobIII2, 2, a, b, c)
end

function getCoefficientsLobIII3(T=Float64)
    a = Array{Dec128}( [[0      0    0 ]
                        [1//4  1//4  0 ]
                        [0     1     0 ]] )
    b = Array{Dec128}( [1//6,  2//3,  1//6] )
    c = Array{Dec128}( [0,     1//2,  1   ] )

    CoefficientsRK(T, :LobIII3, 4, a, b, c)
end

function getCoefficientsLobIII4(T=Float64)
    a = @dec128 [
         [      0             0             0     0 ]
         [ (5+√5)/60          1/6   (15-7*√5)/60  0 ]
         [ (5-√5)/60  (15+7*√5)/60          1/6   0 ]
         [      1/6      (5-√5)/12     (5+√5)/12  0 ]
        ]
    b = @dec128 [1/12,      5/12,      5/12,   1/12]
    c = @dec128 [0,    (5-√5)/10, (5+√5)/10,   1   ]

    CoefficientsRK(T, :LobIII4, 6, a, b, c)
end

function getCoefficientsLobIIIA2(T=Float64)
    a = Array{Dec128}( [[0     0   ]
                        [1//2  1//2]] )
    b = Array{Dec128}( [1//2, 1//2] )
    c = Array{Dec128}( [0,    1   ] )

    CoefficientsRK(T, :LobIIIA2, 2, a, b, c)
end

function getCoefficientsLobIIIA3(T=Float64)
    a = Array{Dec128}( [[0      0     0    ]
                        [5//24  1//3 -1//24]
                        [1//6   2//3  1//6 ]] )
    b = Array{Dec128}( [1//6,  2//3,  1//6] )
    c = Array{Dec128}( [0,     1//2,  1   ] )

    CoefficientsRK(T, :LobIIIA3, 4, a, b, c)
end

function getCoefficientsLobIIIA4(T=Float64)
    a = @dec128 [
         [0             0                0                0           ]
         [(11+√5)/120  (25-   √5)/120  (25-13*√5)/120  (-1+√5)/120]
         [(11-√5)/120  (25+13*√5)/120  (25+   √5)/120  (-1-√5)/120]
         [      1/12            5/12            5/12         1/12 ]
        ]
    b = @dec128 [1/12,      5/12,      5/12,   1/12]
    c = @dec128 [0,    (5-√5)/10, (5+√5)/10,   1   ]

    CoefficientsRK(T, :LobIIIA4, 6, a, b, c)
end

function getCoefficientsLobIIIB2(T=Float64)
    a = Array{Dec128}( [[1//2  0]
                        [1//2  0]] )
    b = Array{Dec128}( [1//2, 1//2] )
    c = Array{Dec128}( [1//2, 1//2] )
    # CORRECTED AN ERROR IN c: it was 0, 1 before

    CoefficientsRK(T, :LobIIIB2, 2, a, b, c)
end

function getCoefficientsLobIIIB3(T=Float64)
    a = Array{Dec128}( [[1//6 -1//6  0   ]
                        [1//6  1//3  0   ]
                        [1//6  5//6  0   ]] )
    b = Array{Dec128}( [1//6,  2//3,  1//6] )
    c = Array{Dec128}( [0,     1//2,  1//1] )

    CoefficientsRK(T, :LobIIIB3, 4, a, b, c)
end

function getCoefficientsLobIIIB4(T=Float64)
    a = @dec128 [
         [ 1/12  (-1-   √5)/24   (-1+   √5)/24    0 ]
         [ 1/12  (25+   √5)/120  (25-13*√5)/120   0 ]
         [ 1/12  (25+13*√5)/120  (25-   √5)/120   0 ]
         [ 1/12  (11-   √5)/24   (11+   √5)/24    0 ]
        ]
    b = @dec128 [1/12,       5/12,      5/12,   1/12]
    c = @dec128 [0,     (5-√5)/10, (5+√5)/10,   1   ]

    CoefficientsRK(T, :LobIIIB4, 6, a, b, c)
end

function getCoefficientsLobIIIC2(T=Float64)
    a = Array{Dec128}( [[1//2 -1//2]
                        [1//2  1//2]] )
    b = Array{Dec128}( [1//2, 1//2] )
    c = Array{Dec128}( [0,    1//1] )

    CoefficientsRK(T, :LobIIIC2, 2, a, b, c)
end

function getCoefficientsLobIIIC3(T=Float64)
    a = Array{Dec128}( [[1//6 -1//3   1//6 ]
                        [1//6  5//12 -1//12]
                        [1//6  2//3   1//6 ]] )
    b = Array{Dec128}( [1//6,  2//3,  1//6] )
    c = Array{Dec128}( [0,     1//2,  1   ] )

    CoefficientsRK(T, :LobIIIC3, 4, a, b, c)
end

function getCoefficientsLobIIIC4(T=Float64)
    a = @dec128 [
         [ 1/12        -√5/12         √5/12   -1/12 ]
         [ 1/12          1/4   (10-7*√5)/60   √5/60 ]
         [ 1/12  (10+7*√5)/60          1/4   -√5/60 ]
         [ 1/12          5/12          5/12    1/12 ]
        ]
    b = @dec128 [1/12,      5/12,      5/12,   1/12]
    c = @dec128 [0,    (5-√5)/10, (5+√5)/10,   1   ]

    CoefficientsRK(T, :LobIIIC4, 6, a, b, c)
end

function getCoefficientsLobIIIF2(T=Float64)
    a = Array{Dec128}( [[1//12 -1//12]
                        [7//12  5//12]] )
    b = Array{Dec128}( [1//2, 1//2] )
    c = Array{Dec128}( [0,    1   ] )

    CoefficientsRK(T, :LobIIIF2, 4, a, b, c)
end

function getCoefficientsLobIIIF3(T=Float64)
    a = Array{Dec128}( [[1//30  -1//15   1//30]
                        [5//24   1//3   -1//24]
                        [2//15  11//15   2//15]] )
    b = Array{Dec128}( [1//6,  2//3,  1//6] )
    c = Array{Dec128}( [0,     1//2,  1   ] )

    CoefficientsRK(T, :LobIIIF3, 6, a, b, c)
end

function getCoefficientsLobIIIF4(T=Float64)
    a = @dec128 [
         [  1/56                 -√5/56           √5/56   -1/56         ]
         [ 37/420+√5/120  5/24-   √5/210  5/24-47*√5/420  -1/210+√5/120 ]
         [ 37/420-√5/120  5/24+47*√5/420  5/24+   √5/210  -1/210-√5/120 ]
         [ 17/168         5/12-   √5/56   5/12+   √5/56   11/168        ]
        ]
    b = @dec128 [1/12,      5/12,      5/12,   1/12]
    c = @dec128 [0,    (5-√5)/10, (5+√5)/10,   1   ]

    CoefficientsRK(T, :LobIIIF4, 8, a, b, c)
end

function getCoefficientsLobIIID2(T=Float64)
    lob  = getCoefficientsLobIII2(Dec128)
    lobC = getCoefficientsLobIIIC2(Dec128)
    CoefficientsRK(T, :LobIIID2, lobC.o, (lob.a + lobC.a)/2, lobC.b, lobC.c)
end

function getCoefficientsLobIIID3(T=Float64)
    lob  = getCoefficientsLobIII3(Dec128)
    lobC = getCoefficientsLobIIIC3(Dec128)
    CoefficientsRK(T, :LobIIID3, lobC.o, (lob.a + lobC.a)/2, lobC.b, lobC.c)
end

function getCoefficientsLobIIID4(T=Float64)
    lob  = getCoefficientsLobIII4(Dec128)
    lobC = getCoefficientsLobIIIC4(Dec128)
    CoefficientsRK(T, :LobIIID4, lobC.o, (lob.a + lobC.a)/2, lobC.b, lobC.c)
end

function getCoefficientsLobIIIE2(T=Float64)
    lobA = getCoefficientsLobIIIA2(Dec128)
    lobB = getCoefficientsLobIIIB2(Dec128)
    CoefficientsRK(T, :LobIIIE2, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
end

function getCoefficientsLobIIIE3(T=Float64)
    lobA = getCoefficientsLobIIIA3(Dec128)
    lobB = getCoefficientsLobIIIB3(Dec128)
    CoefficientsRK(T, :LobIIIE3, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
end

function getCoefficientsLobIIIE4(T=Float64)
    lobA = getCoefficientsLobIIIA4(Dec128)
    lobB = getCoefficientsLobIIIB4(Dec128)
    CoefficientsRK(T, :LobIIIE4, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
end

function getCoefficientsLobIIIG2(T=Float64)
    symplecticize(getCoefficientsLobIIIF2(Dec128); name=:LobIIIG2, T=T)
end

function getCoefficientsLobIIIG3(T=Float64)
    symplecticize(getCoefficientsLobIIIF3(Dec128); name=:LobIIIG3, T=T)
end

function getCoefficientsLobIIIG4(T=Float64)
    symplecticize(getCoefficientsLobIIIF4(Dec128); name=:LobIIIG4, T=T)
end


text:
using MultiInterface
using Test

@testset "MultiInterface" begin
    include("basic.jl")
    # include("validation.jl")
    include("policy.jl")
    # include("wip.jl")
end