text:
# Elliptical Frechet and Marshal Olkim copulas generators

## Elliptical copulas

"""
    Gaussian_cop

Gaussian copula

Fields:
   - Σ - the correlation matrix must be symmetric, positively defined and with ones on diagonal

Constructor

    Gaussian_cop(Σ::Matrix{T}) where T <: Real

```jldoctest

julia> Gaussian_cop([1. 0.5; 0.5 1.])
Gaussian_cop([1.0 0.5; 0.5 1.0])

```
"""
struct Gaussian_cop{T} <: Copula{T}
  Σ::Matrix{T}
  n::Int
  function(::Type{Gaussian_cop})(Σ::Matrix{T}) where T <: Real
    Σ ≈ transpose(Σ) || throw(DomainError("Σ matrix not symmetric"))
    isposdef(Σ) || throw(DomainError("Σ matrix not positivelly defined"))
    prod(diag(Σ)) ≈ 1.0 || throw(DomainError("Σ matrix do not have ones on diagonal"))
    new{T}(Σ,size(Σ,1))
  end
end

function simulate_copula!(U, copula::Gaussian_cop{T}; rng = Random.GLOBAL_RNG) where T
  Σ = copula.Σ
  z = transpose(rand(rng, MvNormal(Σ),size(U,1)))
  for i in 1:size(Σ, 1)
    d = Normal(0, sqrt.(Σ[i,i]))
    z[:,i] = cdf.(d, z[:,i])
  end
  U .= Array(z)
  return nothing
end

"""
    Student_cop

t-Student copula

fields
  - Σ::Matrix{Real} - the correlation matrix must be symmetric, positively defined and with ones on diagonal
  - ν::Int - the parameter n.o. degrees of freedom we require ν > 0

Constructor:

    Student_cop(Σ::Matrix{Real}, ν::Int)

```jldoctest

julia> Student_cop([1. 0.5; 0.5 1.], 4)
Student_cop([1.0 0.5; 0.5 1.0], 4)

```
"""
struct Student_cop{T} <: Copula{T}
  Σ::Matrix{T}
  ν::Int
  n::Int
  function(::Type{Student_cop})(Σ::Matrix{T}, ν::Int) where T <: Real
    Σ ≈ transpose(Σ) || throw(DomainError("Σ matrix not symmetric"))
    isposdef(Σ) || throw(DomainError("Σ matrix not positivelly defined"))
    prod(diag(Σ)) ≈ 1.0 || throw(DomainError("Σ matrix do not have ones on diagonal"))
    ν > 0 || throw(DomainError("ν lower or equal zero"))
    new{T}(Σ, ν,size(Σ,1))
  end
end

function simulate_copula!(U, copula::Student_cop{T}; rng = Random.GLOBAL_RNG) where T
  Σ = copula.Σ
  ν = copula.ν
  z = transpose(rand(rng, MvNormal(Σ),size(U,1)))
  V = rand(rng, Chisq(ν), size(z, 1))
  V = T.(V)
  for i in 1:size(Σ, 1)
    x = z[:,i].*sqrt.(ν./V)./sqrt(Σ[i,i])
    z[:,i] = cdf.(TDist(ν), x)
  end
  U .= Array(z)
  return nothing
end

"""
    Frechet_cop

The Frechet copula

Fileds:
  - n - number of marginals
  - α - the parameter of the maximal copula
  - β - the parameter of the minimal copula

Constructor

    Frechet_cop(n::Int, α::Real)

The one parameter Frechet copula is a combination of maximal copula with  weight α
and independent copula with  weight 1-α.

Constructor

    Frechet_cop(n::Int, α::Real, β::Real)

The two parameters Frechet copula C = α C{max} + β C{min} + (1- α - β) C{⟂}, supported
only for n = 2.

```jldoctest
julia> Frechet_cop(4, 0.5)
Frechet_cop(4, 0.5, 0.0)

julia> Frechet_cop(2, 0.5, 0.3)
Frechet_cop(2, 0.5, 0.3)
```
"""
struct Frechet_cop{T} <: Copula{T}
  n::Int
  α::T
  β::T
  function(::Type{Frechet_cop})(n::Int, α::T) where T <: Real
    0 <= α <= 1 || throw(DomainError("generaton not supported for α ∉ [0,1]"))
    n > 1 || throw(DomainError("n must be greater than 1"))
    new{T}(n, α, 0.)
  end
  function(::Type{Frechet_cop})(n::Int, α::T, β::T) where T <: Real
    0 <= α <= 1 || throw(DomainError("generaton not supported for α ∉ [0,1]"))
    0 <= β <= 1 || throw(DomainError("generaton not supported for β ∉ [0,1]"))
    n == 2 || throw(AssertionError("two parameters Frechet copula supported only for n = 2"))
    0 <= α+β <= 1 || throw(DomainError("α+β must be in range [0,1]"))
    new{T}(n, α, β)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Frechet_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the Frechet copula - Frechet_cop
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> f = Frechet_cop(3, 0.5)
Frechet_cop(3, 0.5, 0.0)

julia> u = zeros(1,3)
1×3 Array{Real,2}:
 0.0  0.0  0.0

julia> Random.seed!(43);

julia> simulate_copula!(u,f)

julia> u
1×3 Array{Real,2}:
 0.180975  0.775377  0.888934
```
"""
function simulate_copula!(U, copula::Frechet_cop{T}; rng = Random.GLOBAL_RNG) where T
  n = copula.n
  α = copula.α
  β = copula.β
  size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
  if (β > 0) & (n == 2)
    for j in 1:size(U,1)
      u_el = rand(rng, T, n)
      frechet_el2!(u_el, α, β, rand(rng, T))
      U[j,:] = u_el
    end
  else
    for j in 1:size(U,1)
      u_el = rand(rng, T, n)
      frechet_el!(u_el, α, rand(rng))
      U[j,:] = u_el
    end
  end
end

"""
  frechet(t::Int, n::Int, α::Real; rng::AbstractRNG)

Given n-variate random data u ∈ R^{t, n}
Returns t realization of n variate data generated from one parameter Frechet_cop(n, α).

```jldoctest
julia> Random.seed!(43);

julia> u = rand(10, 2);

julia> frechet(0.5, u)
10×2 Array{Real,2}:
 0.180975   0.661781
 0.0742681  0.0742681
 0.888934   0.125437
 0.0950087  0.0950087
 0.130474   0.130474
 0.912603   0.740184
 0.828727   0.00463791
 0.400537   0.0288987
 0.521601   0.521601
 0.955881   0.851275
```
"""
function frechet(α::T, u; rng) where T
  for j in 1:size(u, 1)
    v = rand(rng, T)
    el = u[j, :]
    frechet_el!(el, α, v)
    u[j,:] = el
  end
  u
end

"""
  frechet_el!(u::Vector{Real}, α::Real, v::Real)

Given n-variate random vector changes it to such modeled by the two parameters Frechet_cop(n, α, β).
v is the random number form [0,1].
"""
function frechet_el!(u, α, v)
  if (α >= v)
    for i in 1:length(u)-1
      u[i] = u[end]
    end
  end
end

"""
  function frechet_el2!(u::Vector{Real}, α::Real, β::Real, v::Real)

Given bivariate random vector changes it to such modeled by the two parameters Frechet_cop(n, α, β).
v is the random number form [0,1]

"""
function frechet_el2!(u, α, β, v)
  if (α >= v)
    u[1] = u[2]
  elseif (α < v <= α+β)
    u[1] = 1-u[2]
  end
end

### Marshall - Olkin familly

"""
    Marshall_Olkin_cop

Fields:
  - n::Int - number of marginals
  - λ::Vector{Real} - vector of non-negative parameters λₛ, i.e.:
      λ = [λ₁, λ₂, ..., λₙ, λ₁₂, λ₁₃, ..., λ₁ₙ, λ₂₃, ..., λₙ₋₁ₙ, λ₁₂₃, ..., λ₁₂...ₙ]
      and n = ceil(Int, log(2, length(λ)-1)).

Constructor

    Marshall_Olkin_cop(λ)

length(λ) ≧ 3 is required

```jldoctest
julia> Marshall_Olkin_cop([0.5, 0.5, 0.6])
Marshall_Olkin_cop(2, [0.5, 0.5, 0.6])

julia> Marshall_Olkin_cop([0.5, 0.5, 0.6, 0.7, 0.7, 0.7, 0.8])
Marshall_Olkin_cop(3, [0.5, 0.5, 0.6, 0.7, 0.7, 0.7, 0.8])
```
"""
struct Marshall_Olkin_cop{T} <: Copula{T}
  n::Int
  λ::Vector{T}
  function(::Type{Marshall_Olkin_cop})(λ::Vector{T}) where T <: Real
    minimum(λ) >= 0 || throw(AssertionError("all parameters must by >= 0 "))
    length(λ) >= 3 || throw(AssertionError("not supported for length(λ) < 3"))
    n = floor(Int, log(2, length(λ)+1))
    new{T}(n, λ)
  end
end

"""
    simulate_copula!(U::Matrix{Real}, copula::Marshall_Olkin_cop; rng::AbstractRNG = Random.GLOBAL_RNG)

Given the preallocated output U, Returns size(U,1) realizations from the Marshall  Olkin copula - Marshall_Olkin_cop
N.o. marginals is size(U,2), requires size(U,2) == copula.n

```jldoctest
julia> u = zeros(1,2)
1×2 Array{Float64,2}:
 0.0  0.0

julia> cop = Marshall_Olkin_cop([1.,2.,3.])
Marshall_Olkin_cop(2, [1.0, 2.0, 3.0])

julia> Random.seed!(43);

julia> simulate_copula!(u,cop)

julia> u
1×2 Array{Float64,2}:
 0.854724  0.821831
```
"""
function simulate_copula!(U, copula::Marshall_Olkin_cop{T}; rng = Random.GLOBAL_RNG) where T
  λ = copula.λ
  n = copula.n
  size(U, 2) == n || throw(AssertionError("n.o. margins in pre allocated output and copula not equal"))
  s = collect(combinations(1:n))
  for j in 1:size(U,1)
    u = rand(rng, T, 2^n-1)
    U[j,:] = mocopula_el(u, n, λ, s)
  end
end

"""
  mocopula(u::Matrix{Real}, n::Int, λ::Vector{Real})

  Returns: t x n Matrix{Float}, t realizations of n-variate data generated from Marshall-Olkin
  copula with parameter vector λ of non-negative elements λₛ, given [0,1]ᵗˣˡ ∋ u, where
  l = 2ⁿ-1

```jldoctest

  julia> mocopula([0.2 0.3 0.4; 0.3 0.4 0.6; 0.4 0.5 0.7], 2, [1., 1.5, 2.])
  3×2 Array{Float64,2}:
   0.252982  0.201189
   0.464758  0.409039
   0.585662  0.5357

```
"""
function mocopula(u, n, λ)
  T = eltype(u)
  t = size(u,1)
  U = zeros(T, t, n)
  s = collect(combinations(1:n))
  for j in 1:t
      U[j,:] = mocopula_el(u[j,:], n, λ, s)
  end
  U
end

"""
    mocopula_el(u::Vector{Real}, n::Int, λ::Vector{Real}, s::Vector{Vector{Int}})

```jldoctest

julia> mocopula_el([0.1, 0.2, 0.3], 2, [1., 2., 3.], s)
2-element Array{Float64,1}:
 0.20082988502465082
 0.1344421423967149
```
"""
function mocopula_el(u, n, λ, s)
  T = eltype(u)
  l = length(u)
  U = zeros(T, n)
  for i in 1:n
    inds = findall([i in s[k] for k in 1:l])
    x = minimum([-log(u[k])./(λ[k]) for k in inds])
    Λ = sum(λ[inds])
    U[i] = exp.(-Λ*x)
  end
    U
end


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

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

This operation cancels a job. Jobs can be cancelled only when they are in the WAITING state.

# Arguments
- `job_id`: The unique identifier for a job.

"""
cancel_job(JobId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/jobs/$(JobId)"; aws_config=aws_config)
cancel_job(JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/jobs/$(JobId)", params; aws_config=aws_config)

"""
    create_data_set(asset_type, description, name)
    create_data_set(asset_type, description, name, params::Dict{String,<:Any})

This operation creates a data set.

# Arguments
- `asset_type`: The type of file your data is stored in. Currently, the supported asset
  type is S3_SNAPSHOT.
- `description`: A description for the data set. This value can be up to 16,348 characters
  long.
- `name`: The name of the data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Tags"`: A data set tag is an optional label that you can assign to a data set when you
  create it. Each tag consists of a key and an optional value, both of which you define. When
  you use tagging, you can also use tag-based access control in IAM policies to control
  access to these data sets and revisions.
"""
create_data_set(AssetType, Description, Name; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets", Dict{String, Any}("AssetType"=>AssetType, "Description"=>Description, "Name"=>Name); aws_config=aws_config)
create_data_set(AssetType, Description, Name, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AssetType"=>AssetType, "Description"=>Description, "Name"=>Name), params)); aws_config=aws_config)

"""
    create_job(details, type)
    create_job(details, type, params::Dict{String,<:Any})

This operation creates a job.

# Arguments
- `details`: The details for the CreateJob request.
- `type`: The type of job to be created.

"""
create_job(Details, Type; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/v1/jobs", Dict{String, Any}("Details"=>Details, "Type"=>Type); aws_config=aws_config)
create_job(Details, Type, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/v1/jobs", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Details"=>Details, "Type"=>Type), params)); aws_config=aws_config)

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

This operation creates a revision for a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Comment"`: An optional comment about the revision.
- `"Tags"`: A revision tag is an optional label that you can assign to a revision when you
  create it. Each tag consists of a key and an optional value, both of which you define. When
  you use tagging, you can also use tag-based access control in IAM policies to control
  access to these data sets and revisions.
"""
create_revision(DataSetId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets/$(DataSetId)/revisions"; aws_config=aws_config)
create_revision(DataSetId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/v1/data-sets/$(DataSetId)/revisions", params; aws_config=aws_config)

"""
    delete_asset(asset_id, data_set_id, revision_id)
    delete_asset(asset_id, data_set_id, revision_id, params::Dict{String,<:Any})

This operation deletes an asset.

# Arguments
- `asset_id`: The unique identifier for an asset.
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
delete_asset(AssetId, DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)"; aws_config=aws_config)
delete_asset(AssetId, DataSetId, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", params; aws_config=aws_config)

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

This operation deletes a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

"""
delete_data_set(DataSetId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)"; aws_config=aws_config)
delete_data_set(DataSetId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)", params; aws_config=aws_config)

"""
    delete_revision(data_set_id, revision_id)
    delete_revision(data_set_id, revision_id, params::Dict{String,<:Any})

This operation deletes a revision.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
delete_revision(DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)"; aws_config=aws_config)
delete_revision(DataSetId, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)", params; aws_config=aws_config)

"""
    get_asset(asset_id, data_set_id, revision_id)
    get_asset(asset_id, data_set_id, revision_id, params::Dict{String,<:Any})

This operation returns information about an asset.

# Arguments
- `asset_id`: The unique identifier for an asset.
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
get_asset(AssetId, DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)"; aws_config=aws_config)
get_asset(AssetId, DataSetId, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", params; aws_config=aws_config)

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

This operation returns information about a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

"""
get_data_set(DataSetId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)"; aws_config=aws_config)
get_data_set(DataSetId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)", params; aws_config=aws_config)

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

This operation returns information about a job.

# Arguments
- `job_id`: The unique identifier for a job.

"""
get_job(JobId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs/$(JobId)"; aws_config=aws_config)
get_job(JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs/$(JobId)", params; aws_config=aws_config)

"""
    get_revision(data_set_id, revision_id)
    get_revision(data_set_id, revision_id, params::Dict{String,<:Any})

This operation returns information about a revision.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

"""
get_revision(DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)"; aws_config=aws_config)
get_revision(DataSetId, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)", params; aws_config=aws_config)

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

This operation lists a data set's revisions sorted by CreatedAt in descending order.

# Arguments
- `data_set_id`: The unique identifier for a data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
"""
list_data_set_revisions(DataSetId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions"; aws_config=aws_config)
list_data_set_revisions(DataSetId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions", params; aws_config=aws_config)

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

This operation lists your data sets. When listing by origin OWNED, results are sorted by
CreatedAt in descending order. When listing by origin ENTITLED, there is no order and the
maxResults parameter is ignored.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
- `"origin"`: A property that defines the data set as OWNED by the account (for providers)
  or ENTITLED to the account (for subscribers).
"""
list_data_sets(; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets"; aws_config=aws_config)
list_data_sets(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets", params; aws_config=aws_config)

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

This operation lists your jobs sorted by CreatedAt in descending order.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"dataSetId"`: The unique identifier for a data set.
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
- `"revisionId"`: The unique identifier for a revision.
"""
list_jobs(; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs"; aws_config=aws_config)
list_jobs(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/jobs", params; aws_config=aws_config)

"""
    list_revision_assets(data_set_id, revision_id)
    list_revision_assets(data_set_id, revision_id, params::Dict{String,<:Any})

This operation lists a revision's assets sorted alphabetically in descending order.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of results returned by a single call.
- `"nextToken"`: The token value retrieved from a previous call to access the next page of
  results.
"""
list_revision_assets(DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets"; aws_config=aws_config)
list_revision_assets(DataSetId, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets", params; aws_config=aws_config)

"""
    list_tags_for_resource(resource-arn)
    list_tags_for_resource(resource-arn, params::Dict{String,<:Any})

This operation lists the tags on the resource.

# Arguments
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.

"""
list_tags_for_resource(resource_arn; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("GET", "/tags/$(resource-arn)", params; aws_config=aws_config)

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

This operation starts a job.

# Arguments
- `job_id`: The unique identifier for a job.

"""
start_job(JobId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/jobs/$(JobId)"; aws_config=aws_config)
start_job(JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/jobs/$(JobId)", params; aws_config=aws_config)

"""
    tag_resource(resource-arn, tags)
    tag_resource(resource-arn, tags, params::Dict{String,<:Any})

This operation tags a resource.

# Arguments
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.
- `tags`: A label that consists of a customer-defined key and an optional value.

"""
tag_resource(resource_arn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource-arn, tag_keys)
    untag_resource(resource-arn, tag_keys, params::Dict{String,<:Any})

This operation removes one or more tags from a resource.

# Arguments
- `resource-arn`: An Amazon Resource Name (ARN) that uniquely identifies an AWS resource.
- `tag_keys`: The key tags.

"""
untag_resource(resource_arn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_asset(asset_id, data_set_id, name, revision_id)
    update_asset(asset_id, data_set_id, name, revision_id, params::Dict{String,<:Any})

This operation updates an asset.

# Arguments
- `asset_id`: The unique identifier for an asset.
- `data_set_id`: The unique identifier for a data set.
- `name`: The name of the asset. When importing from Amazon S3, the S3 object key is used
  as the asset name. When exporting to Amazon S3, the asset name is used as default target S3
  object key.
- `revision_id`: The unique identifier for a revision.

"""
update_asset(AssetId, DataSetId, Name, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", Dict{String, Any}("Name"=>Name); aws_config=aws_config)
update_asset(AssetId, DataSetId, Name, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)/assets/$(AssetId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name), params)); aws_config=aws_config)

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

This operation updates a data set.

# Arguments
- `data_set_id`: The unique identifier for a data set.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Description"`: The description for the data set.
- `"Name"`: The name of the data set.
"""
update_data_set(DataSetId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)"; aws_config=aws_config)
update_data_set(DataSetId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)", params; aws_config=aws_config)

"""
    update_revision(data_set_id, revision_id)
    update_revision(data_set_id, revision_id, params::Dict{String,<:Any})

This operation updates a revision.

# Arguments
- `data_set_id`: The unique identifier for a data set.
- `revision_id`: The unique identifier for a revision.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Comment"`: An optional comment about the revision.
- `"Finalized"`: Finalizing a revision tells AWS Data Exchange that your changes to the
  assets in the revision are complete. After it's in this read-only state, you can publish
  the revision to your products.
"""
update_revision(DataSetId, RevisionId; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)"; aws_config=aws_config)
update_revision(DataSetId, RevisionId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = dataexchange("PATCH", "/v1/data-sets/$(DataSetId)/revisions/$(RevisionId)", params; aws_config=aws_config)

ext: jl


text:
export iauGc2gd
"""
Transform geocentric coordinates to geodetic using the specified
reference ellipsoid.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical transformation.

Given:
   n       int        ellipsoid identifier (Note 1)
   xyz     double[3]  geocentric vector (Note 2)

Returned:
   elong   double     longitude (radians, east +ve, Note 3)
   phi     double     latitude (geodetic, radians, Note 3)
   height  double     height above ellipsoid (geodetic, Notes 2,3)

Returned (function value):
         int         status:  0 = OK
                              -1 = illegal identifier (Note 3)
                              -2 = internal error (Note 3)

Notes:

   1. The identifier n is a number that specifies the choice of
   reference ellipsoid.  The following are supported:

      n    ellipsoid

      1     WGS84
      2     GRS80
      3     WGS72

   The n value has no significance outside the SOFA software.  For
   convenience, symbols WGS84 etc. are defined in sofam.h.

   2. The geocentric vector (xyz, given) and height (height, returned)
   are in meters.

   3. An error status -1 means that the identifier n is illegal.  An
   error status -2 is theoretically impossible.  In all error cases,
   all three results are set to -1e9.

   4. The inverse transformation is performed in the function iauGd2gc.

Called:
   iauEform     Earth reference ellipsoids
   iauGc2gde    geocentric to geodetic transformation, general

This revision:  2013 September 1

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauGc2gd( int n, double xyz[3],
#                double *elong, double *phi, double *height )

function iauGc2gd(n::Int, xyz::AbstractVector{<:Real})
   ref_elong  = Ref{Float64}(0.0)
   ref_phi    = Ref{Float64}(0.0)
   ref_height = Ref{Float64}(0.0)

   status = ccall((:iauGc2gd, libsofa_c), Cint,
         (Cint, Ptr{Cdouble},
         Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
         convert(Int32, n), xyz,
         ref_elong, ref_phi, ref_height)

   return status, ref_elong[], ref_phi[], ref_height[]
 end


text:
function log_normal_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[])

    closure = OrderedDict()
    closure_exp = OrderedDict() # here it does not play a role
    N = sys.N

    if isempty(binary_vars)
        isbernoulli = false
    else
        sys = bernoulli_moment_eqs(sys, binary_vars)
        isbernoulli = true
    end

    if sys isa CentralMomentEquations
        M = copy(sys.M)
        μ = central_to_raw_moments(N, sys.m_order)
        μ_symbolic = define_μ(N, sys.q_order)
    else
        μ = copy(sys.μ)
        μ_symbolic = copy(μ)
    end

    Σ = Dict()
    for j in 1:N
        for k in j:N
            eⱼ = sys.iter_1[j]
            eₖ = sys.iter_1[k]
            if sys isa CentralMomentEquations
                Σ[(j,k)] = 1. + M[eⱼ .+ eₖ] / μ[eⱼ] / μ[eₖ]
            else
                Σ[(j,k)] = μ[eⱼ .+ eₖ] / μ[eⱼ] / μ[eₖ]
            end
        end
    end

    unique_iter_q = unique(sort(i) for i in sys.iter_q)
    sub = Dict()

    for i in unique_iter_q
        term = prod([μ[sys.iter_1[j]]^i[j] for j in 1:N])
        for j in 1:N
            for k in j+1:N
                term *= Σ[j,k]^(i[j]*i[k])
            end
            term *= Σ[j,j]^(i[j]*(i[j]-1)/2)
        end
        μ[i] = simplify(term)
        closure[μ_symbolic[i]] = μ[i]

        perms = collect(multiset_permutations(i, length(i)))[2:end]

        for iter_perm in perms

            iter_perm_ind = sortperm(sortperm(iter_perm))
            for r in sys.iter_all
                sub[μ_symbolic[r]] = μ_symbolic[r[iter_perm_ind]]
            end

            iter_perm = Tuple(iter_perm)
            μ[iter_perm] = substitute(μ[i], sub)
            closure[μ_symbolic[iter_perm]] = μ[iter_perm]
        end

    end

    if typeof(sys) == CentralMomentEquations
        # construct the corresponding truncated expressions of higher order
        # central moments from the obtained log-normal raw moment expressions
        raw_to_central = raw_to_central_moments(N, sys.q_order, μ, bernoulli=isbernoulli)
        central_to_raw = central_to_raw_moments(N, sys.q_order)
        closure_M = OrderedDict()
        for i in sys.iter_q
            closure_exp[M[i]] = raw_to_central[i]
            closure_M[M[i]] = closure[μ_symbolic[i]]-(central_to_raw[i]-M[i])
            closure_M[M[i]] = simplify(closure_M[M[i]])
        end
        closure = closure_M
    else
        closure_exp = closure
    end

    close_eqs(sys, closure_exp, closure, false)

end


text:
D2ds = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_4^3", [0.0 1.0 0.0; -1.0 0.0 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmad_2", [0.0 1.0 0.0; 1.0 0.0 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [0.0 -1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D3ds = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000003 0.8660254037844384 0.0; -0.8660254037844384 -0.5000000000000003 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmad_2", [1.0 0.0 -0.0; 0.0 -1.0 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_3", [-0.4999999999999998 0.8660254037844392 0.0; 0.8660254037844392 0.49999999999999933 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [-0.5000000000000009 -0.8660254037844383 0.0; -0.8660254037844383 0.5000000000000008 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.5000000000000002 -0.8660254037844387 0.0; -0.8660254037844387 0.5000000000000009 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [-0.4999999999999991 0.8660254037844394 0.0; 0.8660254037844394 0.4999999999999998 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D4ds = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_4^3", [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_8^1", [0.7071067811865476 -0.7071067811865475 0.0; 0.7071067811865475 0.7071067811865476 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_8^3", [-0.7071067811865474 -0.7071067811865477 -0.0; 0.7071067811865477 -0.7071067811865474 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_8^5", [-0.7071067811865479 0.7071067811865471 0.0; -0.7071067811865471 -0.7071067811865479 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_8^7", [0.707106781186547 0.707106781186548 0.0; -0.707106781186548 0.707106781186547 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmad_2", [0.7071067811865475 -0.7071067811865477 -0.0; -0.7071067811865477 -0.7071067811865475 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_3", [0.7071067811865479 0.7071067811865472 0.0; 0.7071067811865472 -0.7071067811865479 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_4", [-0.7071067811865479 0.7071067811865481 0.0; 0.7071067811865481 0.707106781186547 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [-0.7071067811865483 -0.707106781186547 0.0; -0.707106781186547 0.7071067811865481 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(4)", [-1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D5ds = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_5^1", [0.30901699437494745 -0.9510565162951535 0.0; 0.9510565162951535 0.30901699437494745 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^2", [-0.8090169943749473 -0.5877852522924731 0.0; 0.5877852522924731 -0.8090169943749473 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^3", [-0.8090169943749475 0.587785252292473 0.0; -0.587785252292473 -0.8090169943749475 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_5^4", [0.30901699437494734 0.9510565162951535 0.0; -0.9510565162951535 0.30901699437494734 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_10^1", [0.8090169943749475 -0.5877852522924731 0.0; 0.5877852522924731 0.8090169943749475 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_10^3", [-0.3090169943749474 -0.9510565162951536 -0.0; 0.9510565162951536 -0.3090169943749474 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_10^7", [-0.30901699437494756 0.9510565162951538 0.0; -0.9510565162951538 -0.30901699437494756 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_10^9", [0.8090169943749477 0.5877852522924735 0.0; -0.5877852522924735 0.8090169943749477 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmad_2", [0.30901699437494745 -0.9510565162951536 -0.0; -0.9510565162951536 -0.30901699437494745 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_3", [1.0 -0.0 -0.0; -0.0 -1.0 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_4", [0.30901699437494745 0.9510565162951536 0.0; 0.9510565162951536 -0.30901699437494745 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_5", [-0.8090169943749472 0.5877852522924732 0.0; 0.5877852522924732 0.8090169943749473 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [-0.8090169943749472 -0.587785252292473 0.0; -0.587785252292473 0.8090169943749475 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.8090169943749473 0.5877852522924734 0.0; 0.5877852522924734 0.8090169943749481 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [0.3090169943749479 -0.951056516295154 0.0; -0.951056516295154 -0.3090169943749471 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(4)", [0.30901699437494834 0.9510565162951539 0.0; 0.9510565162951539 -0.30901699437494745 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(5)", [-0.8090169943749475 -0.5877852522924731 0.0; -0.5877852522924731 0.8090169943749481 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D6ds = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1]),
    Symmetry.Symel("C_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844388 0.0; 0.8660254037844388 -0.4999999999999998 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_3^2", [-0.5000000000000006 0.8660254037844385 0.0; -0.8660254037844385 -0.5000000000000006 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("S_12^1", [0.8660254037844387 -0.49999999999999994 0.0; 0.49999999999999994 0.8660254037844387 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_12^5", [-0.8660254037844385 -0.5000000000000006 -0.0; 0.5000000000000006 -0.8660254037844385 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_12^7", [-0.866025403784439 0.49999999999999944 0.0; -0.49999999999999944 -0.866025403784439 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_4^3", [0.0 1.0 0.0; -1.0 0.0 -0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("S_12^11", [0.8660254037844383 0.500000000000001 0.0; -0.500000000000001 0.8660254037844383 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("sigmad_2", [0.0 -1.0 -0.0; -1.0 0.0 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_3", [0.8660254037844384 -0.5000000000000004 -0.0; -0.5000000000000004 -0.8660254037844384 -0.0; -0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_4", [0.8660254037844389 0.4999999999999995 0.0; 0.4999999999999995 -0.8660254037844393 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_5", [0.0 1.0000000000000002 0.0; 1.0000000000000002 0.0 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_6", [-0.8660254037844384 0.5000000000000012 0.0; 0.5000000000000012 0.8660254037844379 -0.0; 0.0 -0.0 1.0]),
    Symmetry.Symel("sigmad_1", [-0.8660254037844393 -0.49999999999999883 0.0; -0.49999999999999883 0.8660254037844393 0.0; 0.0 0.0 1.0]),
    Symmetry.Symel("C_2(2)", [-0.4999999999999998 0.8660254037844388 0.0; 0.8660254037844388 0.4999999999999998 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(3)", [-0.5000000000000004 -0.8660254037844385 0.0; -0.8660254037844385 0.5000000000000007 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(4)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(5)", [-0.4999999999999989 0.8660254037844395 0.0; 0.8660254037844395 0.49999999999999933 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(6)", [-0.5000000000000013 -0.8660254037844378 0.0; -0.8660254037844378 0.5000000000000009 0.0; 0.0 0.0 -1.0]),
    Symmetry.Symel("C_2(1)", [1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])]

D2dcn = ["A1", "A2", "B1", "B2", "E"]
D2dct = [1.0  1.0  1.0  1.0  1.0;
         1.0  1.0  1.0 -1.0 -1.0;
         1.0 -1.0  1.0  1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0;
         2.0  0.0 -2.0  0.0  0.0]
D3dcn = ["A1g", "A2g", "Eg", "A1u", "A2u", "Eu"]
D3dct = [1.0  1.0  1.0  1.0  1.0  1.0;
         1.0  1.0 -1.0  1.0  1.0 -1.0;
         2.0 -1.0  0.0  2.0 -1.0  0.0;
         1.0  1.0  1.0 -1.0 -1.0 -1.0;
         1.0  1.0 -1.0 -1.0 -1.0  1.0;
         2.0 -1.0  0.0 -2.0  1.0  0.0]
D4dcn = ["A1", "A2", "B1", "B2", "E1", "E2", "E3"]
D4dct = [1.0  1.0  1.0  1.0  1.0  1.0  1.0;
         1.0  1.0  1.0  1.0  1.0 -1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0  1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0;
         2.0  rt2  0.0 -rt2 -2.0  0.0  0.0;
         2.0  0.0 -2.0  0.0  2.0  0.0  0.0;
         2.0 -rt2  0.0  rt2 -2.0  0.0  0.0]
D5dcn = ["A1g", "A2g", "E1g", "E2g", "A1u", "A2u", "E1u", "E2u"]
D5dct = [1.0  1.0       1.0        1.0  1.0  1.0       1.0        1.0;
         1.0  1.0       1.0       -1.0  1.0  1.0       1.0       -1.0;
         2.0  twocos72  twocos144  0.0  2.0  twocos72  twocos144  0.0;
         2.0  twocos144 twocos72   0.0  2.0  twocos144 twocos72   0.0;
         1.0  1.0       1.0        1.0 -1.0 -1.0      -1.0       -1.0;
         1.0  1.0       1.0       -1.0 -1.0 -1.0      -1.0        1.0;
         2.0  twocos72  twocos144  0.0 -2.0 -twocos72 -twocos144  0.0;
         2.0  twocos144 twocos72   0.0 -2.0 -twocos144 -twocos72   0.0]
D6dcn = ["A1", "A2", "B1", "B2", "E1", "E2", "E3", "E4", "E5"]
D6dct = [1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0;
         1.0  1.0  1.0  1.0  1.0  1.0  1.0 -1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0  1.0 -1.0;
         1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0;
         2.0  rt3  1.0  0.0 -1.0 -rt3 -2.0  0.0  0.0;
         2.0  1.0 -1.0 -2.0 -1.0  1.0  2.0  0.0  0.0;
         2.0  0.0 -2.0  0.0  2.0  0.0 -2.0  0.0  0.0;
         2.0 -1.0 -1.0  2.0 -1.0 -1.0  2.0  0.0  0.0;
         2.0 -rt3  1.0  0.0 -1.0  rt3 -2.0  0.0  0.0]


text:
### functions to analyze stellar clusters
###
### HR diagram metric using Delaunay tessalation to estimate density in the HRD
###
## The tesselation algorithm works with coordinates in range [1,2]
#
function metricHRD(hrdi::Df , label)

    QA = []
    QP = []

    for ilab in label
        nhr = size(ilab)
        points = []

        for n in 1:nhr[1]
            ## Fixed 27.02.2020!!
            push!(points,[hrdi.data[7,ilab[n]],-hrdi.data[6,ilab[n]]])
        end

        ## println("Voronoi tesselation ...")
        ## min_pts is the minimum of pts
        if length(points) > 4
            per, area = voronoi(points)
        else
            per  = [1e10,1e10,1e10,1e10]
            area = [1e10,1e10,1e10,1e10]
        end

        p = sort(per ./ nhr[1])
        a = sort(area ./ nhr[1])

        n10 = max(1, trunc(Int, nhr[1] / 10))
        start10 = max(1 , 1*n10)
        end10   = min(nhr[1] , 9*n10)

        q1a  = mean(a[start10:end10])
        q1p  = mean(p[start10:end10])
        dq1a = std(a[start10:end10])
        dq1p = std(p[start10:end10])

        qa  = -log(max(1e-15, q1a))
        dqa = abs(dq1a/qa)
        qp  = -log(max(1e-15, q1p))
        dqp = abs(dq1p / qp)

        push!(QA, (qa,dqa))
        push!(QP, (qp,dqp))

    end
    return(QA,QP)
end

##  proj :
##  spatial2d : y,z (id 2,3)
##  spatial3d : x,y,z (id 1,2,3)
##  velocity  : vrad,vdec (id 4,5)
##  HRD       : Diagram HR (id 6,7)
##
## labels : labels of the clusters##
##
## !!! works with dfcart !!! (cartesian and normalized)

function metric(s::Df, labels ,proj = "spatial2d", APERTURE = 1.0 , MAXAPERTURE = 15.0, NBOOTSTRAP = 50)

    if proj == "spatial2d"
        ycenter = mean(s.data[2,:])
        zcenter = mean(s.data[3,:])
    elseif proj == "spatial3d"
        xcenter = mean(s.data[1,:])
        ycenter = mean(s.data[2,:])
        zcenter = mean(s.data[3,:])
    elseif proj == "velocity"
        vxcenter = mean(s.data[4,:])
        vycenter = mean(s.data[5,:])
    elseif proj == "HRD"
        res = metricHRD(s, labels)
        return(res)
    end

    aper2 = APERTURE * APERTURE
    Q = []
    for ilab in labels
        slab = s.data[:,ilab]

        if proj == "spatial2d"
            yy = slab[2,:]
            zz = slab[3,:]
            ycenter = mean(yy)
            zcenter = mean(zz)
            angle_out = 2π * rand(Float64, NBOOTSTRAP)
            rad_out   = (MAXAPERTURE - 2APERTURE) * rand(Float64, NBOOTSTRAP)
            angle_in  = 2π * rand(Float64, NBOOTSTRAP)
            rad_in    = APERTURE * randexp(Float64, NBOOTSTRAP)
        elseif proj == "spatial3d"
            xx = slab[1,:]
            yy = slab[2,:]
            zz = slab[3,:]
            xcenter = mean(xx)
            ycenter = mean(yy)
            zcenter = mean(zz)
            angle_out = 2π * rand(Float64, NBOOTSTRAP)
            phi_out   = π * rand(Float64, NBOOTSTRAP)
            rad_out   = (MAXAPERTURE - 2APERTURE) * rand(Float64, NBOOTSTRAP)
            angle_in  = 2π * rand(Float64, NBOOTSTRAP)
            phi_in    = π * rand(Float64, NBOOTSTRAP)
            rad_in    = APERTURE * randexp(Float64, NBOOTSTRAP)
        elseif proj == "velocity"
            vx = slab[4,:]
            vy = slab[5,:]
            vxcenter = mean(vx)
            vycenter = mean(vy)
            angle_out = 2π * rand(Float64, NBOOTSTRAP)
            rad_out   = (MAXAPERTURE - 2APERTURE) * rand(Float64, NBOOTSTRAP)
            angle_in  = 2π * rand(Float64, NBOOTSTRAP)
            rad_in    = APERTURE * rand(Float64, NBOOTSTRAP)
        end

        qc = zeros(NBOOTSTRAP)

        for k in 1:NBOOTSTRAP
            if proj == "spatial2d"
                yout = ycenter + rad_out[k] * cos(angle_out[k])
                zout = zcenter + rad_out[k] * sin(angle_out[k])
                radii_out = (slab[2,:] .- yout) .* (slab[2,:] .- yout) .+ (slab[3,:] .- zout) .* (slab[3,:] .- zout)
                yin = ycenter + rad_in[k] * cos(angle_in[k])
                zin = zcenter + rad_in[k] * sin(angle_in[k])
                radii_in = (slab[2,:] .- yin) .* (slab[2,:] .- yin) .+ (slab[3,:] .- zin) .* (slab[3,:] .- zin)
            elseif proj == "spatial3d"
                xout = xcenter + rad_out[k] * cos(angle_out[k]) *  cos(phi_out[k])
                yout = ycenter + rad_out[k] * sin(angle_out[k]) *  cos(phi_out[k])
                zout = zcenter + rad_out[k] * sin(phi_out[k])
                radii_out = (slab[1,:] .- xout) .* (slab[1,:] .- xout) .+ (slab[2,:] .- yout) .* (slab[2,:] .- yout).+ (slab[3,:] .- zout) .* (slab[3,:] .- zout)
                xin = xcenter + rad_in[k] * cos(angle_in[k]) *  cos(phi_in[k])
                yin = ycenter + rad_in[k] * sin(angle_in[k]) *  cos(phi_in[k])
                zin = zcenter + rad_in[k] * sin(phi_in[k])
                radii_in = (slab[1,:] .- xin) .* (slab[1,:] .- xin) .+ (slab[2,:] .- yin) .* (slab[2,:] .- yin).+ (slab[3,:] .- zin) .* (slab[3,:] .- zin)
            elseif proj == "velocity"
                vxout = vxcenter + rad_out[k] * cos(angle_out[k])
                vyout = vycenter + rad_out[k] * sin(angle_out[k])
                radii_out = (slab[4,:] .- vxout) .* (slab[4,:] .- vxout) .+ (slab[5,:] .- vyout) .* (slab[5,:] .- vyout)
                vxin = vxcenter + rad_in[k] * cos(angle_in[k])
                vyin = vycenter + rad_in[k] * sin(angle_in[k])
                radii_in = (slab[4,:] .- vxin) .* (slab[4,:] .- vxin) .+ (slab[5,:] .- vyin) .* (slab[5,:] .- vyin)
            end

            isdout = radii_out .< aper2
            nout   = length(radii_out[isdout])
            isin   = radii_in .< aper2
            nin    = length(radii_in[isin])

            qc[k] = log(max(1,nin)) - log(max(1,nout))

        end

        Qc    = mean(qc)
        Q_std = std(qc)
        #println("Q: $Qc -- Q_std : $Q_std")
        push!(Q,(Qc,Q_std))
    end
    return(Q)
end

### DBSCAN
function clusters(data , epsilon, leaf , minneigh, mincluster)
    eps = epsilon
    leafsize = leaf
    min_neighbors = minneigh
    min_cluster_size = mincluster

    res = dbscan(data , eps , leafsize = leaf, min_neighbors = minneigh, min_cluster_size=mincluster)

    label = Vector{Vector{Int}}()

    for cl in res
        indx = cl.core_indices
        append!(indx, cl.boundary_indices)
        push!(label,indx)
    end
    return(label)
end

## find clusters from dbscan with a metric
##
## compute the dbscan clusters, the metric and the Metropolis acceptance
##
## df: Gaia struct (data.jl), normally a weighted cartersian df (dfcartnorm)
## dfcart: Gaia struct in cartesian w/o any transformation (dfcart)
##
## !! *qq* : composite metric. The weight could be adjusted.
## !! *qstar* : number of stars in the highest qq
## m: dbscan parameters.
##
## WARNING: the default values for APERTURE
## Check the metric(2) function used...
## Note 2.3.2020: the HRD metric is discarded and the weight scaling is down 6

### find_clusters with meta parameters.
###
function find_clusters2(df::GaiaClustering.Df, dfcart::GaiaClustering.Df , m::GaiaClustering.modelfull,
    param::GaiaClustering.meta, verbose= false)
    let
        labels = clusters(df.data , m.eps , 20, m.min_nei, m.min_cl)
        nsol= length(labels)
        if nsol == 0 || nsol > param.clustermax
            if verbose println("### Warning... $nsol clusters found with find_clusters2, returns 0 ") end
            return(0, 0)
        end

    ### metrics of the clusters
        q2d = metric2(dfcart, labels, "spatial2d" , param.aperture2d, param.maxaperture2d, param.nboot)
        q3d = metric2(dfcart, labels, "spatial3d" , param.aperture3d, param.maxaperture3d, param.nboot)
        qv  = metric2(dfcart, labels, "velocity" , param.aperturev, param.maxaperturev, param.nboot)
        qp, qa = metric2(dfcart, labels, "HRD" )

        nlab = []
        for ilab in labels
            push!(nlab,length(ilab))
        end

    #### metric for the number of stars in the cluster
        qn = []
        for nl in nlab
            push!(qn,log10(nl))
        end

        qc = 0.
        qstar = 0
        for i in 1:length(nlab)
            k1 = q2d[i][1]
            k1bis = q3d[i][1]
            k2 = qv[i][1]
            k3 = qa[i][1]
            k4 = qn[i]
            ############### Composite metric ###
            qq = (2k1 + k1bis + 3k2 + k3 + k4) / 8.0
            # qq = (3k1 + k1bis + 3k2 + k4) / 8.0
            ###############
            if qq > qc
                qc = qq
                qstar = nlab[i]
            end
        end
        return(qc, qstar)
    end
end
## label from the dbscan labels with maximum stars
function find_cluster_label(labels)
    let
    println("## Selecting best cluster based on Nmax...")
    i = 1 ; nmax = 0 ; ilabel = 0
        for ilab in labels
            nlab = length(ilab)
            if nlab > nmax
                ilabel = i
                nmax= nlab
            end
            i += 1
        end
    return(ilabel, nmax)
    end
end

## label with higher Qc is selected
function find_cluster_label2(labels, df::GaiaClustering.Df, dfcart::GaiaClustering.Df ,
    m::GaiaClustering.meta)
    let

        if length(labels) == 0
            println("## no clusters in the DBSCAN selection...")
            return(0, 0, 0)
        end
        ### metrics of the clusters
        q2d = metric2(dfcart, labels, "spatial2d" , m.aperture2d, m.maxaperture2d, m.nboot)
        q3d = metric2(dfcart, labels, "spatial3d" , m.aperture3d, m.maxaperture3d, m.nboot)     #### Added
        qv  = metric2(dfcart, labels, "velocity" , m.aperturev, m.maxaperturev, m.nboot)
        qp, qa = metric2(dfcart, labels, "HRD" )

        nlab = []
        for ilab in labels
            push!(nlab,length(ilab))
        end

        #### metric for the number of stars in the cluster
        qn = []
        for nl in nlab
            push!(qn,log10(nl))
        end

        qc= []
        for i in 1:length(nlab)
            k1 = q2d[i][1]
            k1bis = q3d[i][1]
            k2 = qv[i][1]
            k3 = qa[i][1]
            k4 = qn[i]
            ############### Composite metric ###
            qq = (2k1 + k1bis + 3k2 + k3 + k4) / 8.0
            # qq = (3k1 + k1bis + 3k2 + k4) / 8.0
            ###############
            push!(qc,qq)
        end

        if m.labels == "Qc"
            println("## Selecting best cluster based on Qc")
            println("## Qc: $qc")
            println("## Qn: $nlab")
            bestlabel= findmax(qc)[2]
        elseif m.labels == "Qn"
            println("## Selecting best cluster based on Qn")
            println("## Qc: $qc")
            println("## Qn: $nlab")
            bestlabel= findmax(nlab)[2]
        elseif m.labels == "QcQn"
            println("## Selecting best cluster based on Qc and Qn")
            println("## Qc: $qc")
            println("## Qn: $nlab")
            qcqn= get_qcqn.(qc, nlab, (m , ))
            println("## QcQn: $qcqn")
            bestlabel= findmax(qcqn)[2]
        elseif m.labels == "QcQnhigh"
            println("## Selecting best cluster based on Qc and maxQn. Make sure that is what you want...")
            println("## Qc: $qc")
            println("## Qn: $nlab")
            mtemp= m
            mtemp.minQn= 30
            qcqn= get_qcqn.(qc, nlab, (mtemp , ))
            println("## QcQn: $qcqn")
            bestlabel= findmax(qcqn)[2]
        else
            println("## No rules for selecting best cluster, label 1 set by default...")
            bestlabel= 1
        end

        return(bestlabel, nlab[bestlabel], qc[bestlabel])
    end
end
## function to project (Qc,Qn) values on the best value with QcQn methd
## qn = n
## favors qc for small qn otherwise qc
function sigmoid(x, c1, c2)
    f= 1/(1+exp(-c1*(x-c2)))
    return(f)
end
function get_qcqn(qc, qn,m::GaiaClustering.meta)
    k= qc*sigmoid(qn, 0.1, m.minQn)*(1-sigmoid(qn,0.1,m.maxQn))
    # k= qc*sigmoid(qn, 0.1, m.minQn)
    return(k)
end
## compute the properties of the cluster with indices indx
function get_properties_SC(indx, df::GaiaClustering.Df, dfcart::GaiaClustering.Df)::SCproperties
    nstars   = length(indx)
    distance = mean(df.data[3,indx])
    l        = mean(df.data[1,indx])
    b        = mean(df.data[2,indx])
    ra       = mean(df.raw[1,indx])
    dec      = mean(df.raw[2,indx])
    vl       = mean(df.data[4,indx])
    vb       = mean(df.data[5,indx])
    xdisp    = std(dfcart.data[1,indx])
    ydisp    = std(dfcart.data[2,indx])
    zdisp    = std(dfcart.data[3,indx])
    vldisp   = std(df.data[4,indx])
    vbdisp   = std(df.data[5,indx])
    parallax = mean(df.raw[5,indx])
    pmra     = mean(df.raw[6,indx])
    pmdec    = mean(df.raw[7,indx])
    pml      = mean(df.raw[8,indx])
    pmb      = mean(df.raw[9,indx])

    vrad     = 0.
    vraddisp = 0.
    indvrad = isnotnan(df.raw[13,indx])

    if length(indvrad) >  1
        vrad = mean(df.raw[13,indx[indvrad]])
        vraddisp = std(df.raw[13,indx[indvrad]])
    elseif length(indvrad) == 1
        vrad     = mean(df.raw[13,indx[indvrad]])
        vraddisp = 0.
    end

    sc = SCproperties(nstars , distance, ra , dec , l , b , parallax, pmra , pmdec , pml, pmb, vl , vb, vrad, xdisp ,
        ydisp , zdisp , vldisp, vbdisp, vraddisp)
    return(sc)

end

## compute the properties of the cluster with indices indx
## new extended version of SCproperties.
## mean -> median
function get_properties_SC2(indx, df::GaiaClustering.Df, dfcart::GaiaClustering.Df)::SCproperties2
    nstars   = length(indx)
    distance = median(df.data[3,indx])
    l        = median(df.data[1,indx])
    b        = median(df.data[2,indx])
    ra       = median(df.raw[1,indx])
    dec      = median(df.raw[2,indx])
    vl       = median(df.data[4,indx])
    vb       = median(df.data[5,indx])
    xdisp    = std(dfcart.data[1,indx])
    ydisp    = std(dfcart.data[2,indx])
    zdisp    = std(dfcart.data[3,indx])
    vldisp   = std(df.data[4,indx])
    vbdisp   = std(df.data[5,indx])
    parallax = median(df.raw[5,indx])
    pmra     = median(df.raw[6,indx])
    pmdec    = median(df.raw[7,indx])
    pml      = median(df.raw[8,indx])
    pmb      = median(df.raw[9,indx])

    vrad     = 0.
    vraddisp = 0.
    indvrad = isnotnan(df.raw[13,indx])

    if length(indvrad) >  1
        vrad = mean(df.raw[13,indx[indvrad]])
        vraddisp = std(df.raw[13,indx[indvrad]])
    elseif length(indvrad) == 1
        vrad     = mean(df.raw[13,indx[indvrad]])
        vraddisp = 0.
    end

    doff= sqrt(median(dfcart.data[2,indx])^2+median(dfcart.data[3,indx])^2)
    doffdeg= atand(doff/distance)

    edgeratio1, edgeratio2= edge_ratio(dfcart, indx)

    sc = SCproperties2()
    sc.nstars= nstars
    sc.distance= distance
    sc.ra= ra
    sc.dec= dec
    sc.l= l
    sc.b= b
    sc.parallax= parallax
    sc.pmra= pmra
    sc.pmdec= pmdec
    sc.pml= pml
    sc.pmb= pmb
    sc.vl= vl
    sc.vb= vb
    sc.vrad= vrad
    sc.xdisp= xdisp
    sc.ydisp= ydisp
    sc.zdisp= zdisp
    sc.vldisp= vldisp
    sc.vbdisp= vbdisp
    sc.vraddisp= vraddisp
    sc.offdeg= doffdeg
    sc.edgratg= edgeratio1
    sc.edgratm= edgeratio2

    return(sc)
end
## some of the parameters are changed in find_clusters..
function metric2(s::GaiaClustering.Df, labels ,proj = "spatial2d", APERTURE = 1.0 ,
        MAXAPERTURE = 15.0, NBOOTSTRAP = 50 , MAXDISP2D = 10.,  MAXDISP3D = 30. , MAXDISPVEL = 4.)

    ### probabilities of the dispersion..
    p2d  = Normal(0., MAXDISP2D)
    p3d  = Normal(0., MAXDISP3D)
    pvel = Normal(0., MAXDISPVEL)

    if proj == "spatial2d"
        ycenter = mean(s.data[2,:])
        zcenter = mean(s.data[3,:])
    elseif proj == "spatial3d"
        xcenter = mean(s.data[1,:])
        ycenter = mean(s.data[2,:])
        zcenter = mean(s.data[3,:])
    elseif proj == "velocity"
        vxcenter = mean(s.data[4,:])
        vycenter = mean(s.data[5,:])
    elseif proj == "HRD"
        res = metricHRD(s, labels)
        return(res)
    end

    aper2 = APERTURE * APERTURE
    Q = []
    for ilab in labels
        slab = s.data[:,ilab]

        if proj == "spatial2d"
            yy = slab[2,:]
            zz = slab[3,:]
            ycenter = mean(yy)
            zcenter = mean(zz)
            angle_out = 2π * rand(Float64, NBOOTSTRAP)
            rad_out   = (MAXAPERTURE - 2APERTURE) * rand(Float64, NBOOTSTRAP)
            angle_in  = 2π * rand(Float64, NBOOTSTRAP)
            rad_in    = APERTURE * randexp(Float64, NBOOTSTRAP)
            ydisp = std(slab[2,:])
            zdisp = std(slab[3,:])
            prob2d = pdf(p2d,sqrt(ydisp^2+zdisp^2)) / pdf(p2d,0.)
        elseif proj == "spatial3d"
            xx = slab[1,:]
            yy = slab[2,:]
            zz = slab[3,:]
            xcenter = mean(xx)
            ycenter = mean(yy)
            zcenter = mean(zz)
            angle_out = 2π * rand(Float64, NBOOTSTRAP)
            phi_out   = π * rand(Float64, NBOOTSTRAP)
            rad_out   = (MAXAPERTURE - 2APERTURE) * rand(Float64, NBOOTSTRAP)
            angle_in  = 2π * rand(Float64, NBOOTSTRAP)
            phi_in    = π * rand(Float64, NBOOTSTRAP)
            rad_in    = APERTURE * randexp(Float64, NBOOTSTRAP)
            xdisp = std(slab[1,:])
            ydisp = std(slab[2,:])
            zdisp = std(slab[3,:])
            prob3d = pdf(p3d,sqrt(xdisp^2+ydisp^2+zdisp^2)) / pdf(p3d,0.)
        elseif proj == "velocity"
            vx = slab[4,:]
            vy = slab[5,:]
            vxcenter = mean(vx)
            vycenter = mean(vy)
            angle_out = 2π * rand(Float64, NBOOTSTRAP)
            rad_out   = (MAXAPERTURE - 2APERTURE) * rand(Float64, NBOOTSTRAP)
            angle_in  = 2π * rand(Float64, NBOOTSTRAP)
            rad_in    = APERTURE * rand(Float64, NBOOTSTRAP)
            vxdisp = std(slab[4,:])
            vydisp = std(slab[5,:])
            probvel = pdf(pvel,sqrt(vxdisp^2+vydisp^2)) / pdf(pvel,0.)
        end

        qc = zeros(NBOOTSTRAP)

        for k in 1:NBOOTSTRAP
            if proj == "spatial2d"
                prob =  prob2d
                yout = ycenter + rad_out[k] * cos(angle_out[k])
                zout = zcenter + rad_out[k] * sin(angle_out[k])
                radii_out = (slab[2,:] .- yout) .* (slab[2,:] .- yout) .+ (slab[3,:] .- zout) .* (slab[3,:] .- zout)
                yin = ycenter + rad_in[k] * cos(angle_in[k])
                zin = zcenter + rad_in[k] * sin(angle_in[k])
                radii_in = (slab[2,:] .- yin) .* (slab[2,:] .- yin) .+ (slab[3,:] .- zin) .* (slab[3,:] .- zin)
            elseif proj == "spatial3d"
                prob = prob3d
                xout = xcenter + rad_out[k] * cos(angle_out[k]) *  cos(phi_out[k])
                yout = ycenter + rad_out[k] * sin(angle_out[k]) *  cos(phi_out[k])
                zout = zcenter + rad_out[k] * sin(phi_out[k])
                radii_out = (slab[1,:] .- xout) .* (slab[1,:] .- xout) .+ (slab[2,:] .- yout) .* (slab[2,:] .- yout).+ (slab[3,:] .- zout) .* (slab[3,:] .- zout)
                xin = xcenter + rad_in[k] * cos(angle_in[k]) *  cos(phi_in[k])
                yin = ycenter + rad_in[k] * sin(angle_in[k]) *  cos(phi_in[k])
                zin = zcenter + rad_in[k] * sin(phi_in[k])
                radii_in = (slab[1,:] .- xin) .* (slab[1,:] .- xin) .+ (slab[2,:] .- yin) .* (slab[2,:] .- yin).+ (slab[3,:] .- zin) .* (slab[3,:] .- zin)
            elseif proj == "velocity"
                prob = probvel
                vxout = vxcenter + rad_out[k] * cos(angle_out[k])
                vyout = vycenter + rad_out[k] * sin(angle_out[k])
                radii_out = (slab[4,:] .- vxout) .* (slab[4,:] .- vxout) .+ (slab[5,:] .- vyout) .* (slab[5,:] .- vyout)
                vxin = vxcenter + rad_in[k] * cos(angle_in[k])
                vyin = vycenter + rad_in[k] * sin(angle_in[k])
                radii_in = (slab[4,:] .- vxin) .* (slab[4,:] .- vxin) .+ (slab[5,:] .- vyin) .* (slab[5,:] .- vyin)
            end

            isdout = radii_out .< aper2
            nout   = length(radii_out[isdout])
            isin   = radii_in .< aper2
            nin    = length(radii_in[isin])

            qc[k] = prob * (log(max(1,nin)) - log(max(1,nout)))

        end

        Qc    = mean(qc)
        Q_std = std(qc)
        #println("Q: $Qc -- Q_std : $Q_std")
        push!(Q,(Qc,Q_std))
    end
    return(Q)
end

### get the metrics
### still a list of label lists...
##
function get_metrics(labels, dfcart::GaiaClustering.Df , m::meta)
    let
        #### metrics of the clusters
        q2d = metric2(dfcart, labels, "spatial2d" , m.aperture2d, m.maxaperture2d, m.nboot)
        q3d = metric2(dfcart, labels, "spatial3d" , m.aperture3d, m.maxaperture3d, m.nboot)
        qv  = metric2(dfcart, labels, "velocity" ,  m.aperturev, m.maxaperturev, m.nboot)
        qp, qa = metric2(dfcart, labels, "HRD" )

        nlab = []
        for ilab in labels
            push!(nlab,length(ilab))
        end

        #### metric for the number of stars in the cluster
        qn = []
        for nl in nlab
            push!(qn,log10(nl))
        end

        qc = 0.
        qstar = 0
        for i in 1:length(nlab)
            k1 = q2d[i][1]
            k1bis = q3d[i][1]
            k2 = qv[i][1]
            k3 = qa[i][1]
            k4 = qn[i]
            ############### Composite metric ###
            qq = (2k1 + k1bis + 3k2 + k3 + k4) / 8.0
            # qq = (3k1 + k1bis + 3k2 + k4) / 8.0
            ###############
            if qq > qc
                qc = qq
                qstar = nlab[i]
            end
        end
        return(qc, qstar)
    end
end

#### cycle extraction
#### namefile: prefix for plotfile
function cycle_extraction(df::GaiaClustering.Df, dfcart::GaiaClustering.Df, m::GaiaClustering.meta)
    let
        println("############### cycle_extraction #########")

        votname= m.votname
        cyclerun= true ; cycle= 1 ; FLAG= 0

        sclist= [] ; mcmclist= [] ; perflist= [] ; chainlist= []

        println("##")
        while cyclerun
            FLAG= -1
            tstart= now()
            println("###############")
            print("## "); println(blue("starting cycle $cycle ..."))
            @printf("## starting time: %s \n",tstart)
            ## extraction one cycle.. MCMC optimization
            mc , iter, FLAGmcmc= abc_mcmc_dbscan_full2(dfcart, m)
            println("## ABC/MCMC flag: $FLAGmcmc")
            nchain= length(mc.qc)
            println("## $iter iterations performed...")
            println("## $nchain chains")

            if FLAGmcmc== -1 || nchain > m.minchainreached
                println("## optimization completed..")
                println("## analyzing solutions...")
                plot_dbscanfull_mcmc(m.plotdir, "$votname.$cycle", mc , false)

                ## get the cluster and plot it
                println("## extracting the cluster using DBSCAN/WEIGHTING with:")
                res= extraction_mcmc(mc, m.votname)
                eps= res.epsm[1]
                min_nei= trunc(Int,res.mneim[1] + 0.5)
                min_cl= trunc(Int,res.mclm[1] + 0.5)
                w3d= res.w3dm[1]
                wvel= res.wvelm[1]
                whrd= res.whrdm[1]

                mres = GaiaClustering.modelfull(eps,min_nei,min_cl,w3d,wvel,whrd)
                dfcartnorm = getDfcartnorm(dfcart, mres)
                labels = clusters(dfcartnorm.data ,eps  , 20, min_nei, min_cl)
                labelmax , nmax, qc = find_cluster_label2(labels, df, dfcart, m)
                println("## label $labelmax written to oc...")
                export_df("$votname.$cycle", m.ocdir, df , dfcart , labels , labelmax)

                ## Principal components
                pc, pcres= compute_PC(df, dfcart, labels, labelmax)

                edgeratio1, edgeratio2= edge_ratio(dfcart, labels[labelmax])
                scproperties = get_properties_SC2(labels[labelmax] , df, dfcart)
                scdf= convertStruct2Df(scproperties)
                insertcols!(scdf, 1, :votname => votname)
                s=size(scdf)
                insertcols!(scdf, 2, :cycle => cycle)
                insertcols!(scdf, 3, :pc3 => pcres[3])
                insertcols!(scdf, 3, :pc2 => pcres[2])
                insertcols!(scdf, 3, :pc1 => pcres[1])

                insertcols!(res, 2,  :cycle => cycle)
                push!(sclist, scdf)
                push!(mcmclist, res)

                ## create DF chain
                dfchain= create_DFchain(mc, votname, cycle)
                push!(chainlist,dfchain)

                println("###")
                println("### solution label: $labelmax")
                print("### "); println(red(@sprintf("PC1: %3.1f , PC2: %3.1f , PC3: %3.1f", pcres[1], pcres[2], pcres[3])))
                println("### Offdeg: $(scproperties.offdeg)")
                println("### Edge ratio: $(scproperties.edgratm)")
                println("### N stars: $nmax")
                println("### Qc: $qc")
                println("###")

                k= score_cycle(qc, nmax, nchain, iter)
                @printf("## score cycle %d: %3.3f \n",cycle, k)

                extraplot= DataFrame(cycle=cycle, score_cycle=k, qc=qc, votname=votname, pc1=pcres[1],pc2=pcres[2], pc3=pcres[3])

                plot_cluster2(m.plotdir, "$votname.$cycle", labels[labelmax], scproperties,
                    dfcart , false, extraplot)

                jump= 50  # how many stars to jump in the plot
                plot_rawdata(m.plotdir, "$votname.$cycle", labels[labelmax], scproperties,
                    dfcart , pc, jump, false, extraplot)

                println("###")
                println("### subtracting BEST solution from Df...")
                dfnew, dfcartnew= remove_stars(df, dfcart, labels[labelmax])
                df= dfnew
                dfcart= dfcartnew

                ########################### STOP conditions #########
                FLAG= 0
                if nmax < m.minstarstop
                    FLAG= FLAG | (1<<0)
                    println("### extraction stopped at cycle $cycle")
                    println("### nmax too low...")
                    cyclerun= false
                end
                if cycle == m.cyclemax
                    FLAG= FLAG | (1<<1)
                    println("### extraction stopped at cycle $cycle")
                    println("### cyclemax reached...")
                    cyclerun= false
                end
                if qc < m.qcminstop
                    FLAG= FLAG | (1<<2)
                    println("### extraction stopped at cycle $cycle")
                    println("### Qc too low...")
                    cyclerun= false
                end
                if w3d/wvel < m.wratiominstop || wvel/w3d < m.wratiominstop
                    FLAG= FLAG | (1<<3)
                    println("### extraction stopped at cycle $cycle")
                    println("### weight ratio too low...")
                    cyclerun= false
                end
                if FLAGmcmc == 3 && nchain > m.minchainreached
                    FLAG= FLAG | (1<<4)
                    println("## extraction stopped at cycle $cycle")
                    println("## chain iteration not performed completely but sufficient to keep...")
                    cyclerun= false
                end
                #####################################################
                ## Time
                tend= now()
                duration= Dates.value(tend-tstart) / (1000*1)
                nstar= size(df.raw)[2]
                timeperiterstar= duration / (iter*nstar)
                timeperchainstar= duration / (nchain*nstar)
                @printf("## \n")
                @printf("## Time: \n")
                @printf("## duration per cycle %3.3f sec \n", duration)
                @printf("## duration per iteration*star %3.3e sec \n", timeperiterstar)
                @printf("## duration per chain*star %3.3e sec \n", timeperchainstar)
                @printf("##\n")

                ## log the results of performances
                dfout= DataFrame(votname=votname, cycle=cycle, nstar=nstar, qc=qc, nmax=nmax, nchain=nchain, iter=iter,
                scorecycle=k, duration=duration, timeperiterstar=timeperiterstar ,
                timeperchainstar= timeperchainstar )
                push!(perflist, dfout)

                cycle += 1
            else
                println("## nothing found, stopped...")
                FLAG= 0
                cyclerun= false
            end
        end
        if cycle >= 2
            save_cycle(sclist, mcmclist, perflist, chainlist, m)
        end
        return(cycle-1, FLAG)
    end
end

function score_cycle(qc, qn, nchain, iter)
    k= log10(qc*qn*nchain /iter)
    return(k)
end

function remove_stars(df::GaiaClustering.Df, dfcart::GaiaClustering.Df, idx)
    s0=size(df.data)

    diff= setdiff(1:s0[2],idx)
    dfdata= df.data[:,setdiff(1:end,idx)]
    dfraw= df.raw[:,setdiff(1:end,idx)]
    dferr= df.err[:,setdiff(1:end,idx)]
    dfsid= df.sourceid[:,setdiff(1:end,idx)]

    dfcartdata= dfcart.data[:,setdiff(1:end,idx)]
    dfcartraw= dfcart.raw[:,setdiff(1:end,idx)]
    dfcarterr= dfcart.err[:,setdiff(1:end,idx)]
    dfcartsid= dfcart.sourceid[:,setdiff(1:end,idx)]

    s=size(dfdata)

    dfnew= GaiaClustering.Df(s[2],dfdata,dfraw,dferr,dfsid)
    dfcartnew= GaiaClustering.Df(s[2],dfcartdata,dfcartraw,dfcarterr,dfcartsid)

    nrem= length(idx)
    println("### $nrem stars removed")
    return(dfnew, dfcartnew)
end

#### cycle extraction
#### namefile: prefix for plotfile
### adding boolean optim for optimization or not
function cycle_extraction_optim(df::GaiaClustering.Df, dfcart::GaiaClustering.Df, m::GaiaClustering.meta, optim=true)
    let
        println("############### cycle_extraction #########")
        println("## maximum cycle: $(m.cyclemax)")
        println("## Q metric on DBSCAN solutions: $(m.labels)")
        if optim
            println("## DBSCAN/weighting optimization. It may take time...")
        else
            println("## NO optimization, the weightings and DBSCAN parameters are fixed.")
            println("## w3d : $(m.w3d)")
            println("## wvel : $(m.wvel)")
            println("## whrd : $(m.whrd)")
            println("## ϵ : $(m.eps)")
            println("## min_cluster : $(m.mcl)")
            println("## min_neighbor : $(m.mnei)")
        end
        votname= basename(m.votname)       # get rid of the leading part only useful to read data
        cyclerun= true ; cycle= 1 ; FLAG= 0

        sclist= [] ; mcmclist= [] ; perflist= [] ; chainlist= []

        println("##")
        while cyclerun
            FLAG= -1
            tstart= now()
            println("###############")
            print("## "); println(blue("starting cycle $cycle ..."))
            @printf("## starting time: %s \n",tstart)

            if optim
                ## extraction one cycle.. MCMC optimization
                mc , iter, FLAGmcmc= abc_mcmc_dbscan_full2(dfcart, m)
                println("## ABC/MCMC flag: $FLAGmcmc")
                nchain= length(mc.qc)
                println("## $iter iterations performed...")
                println("## $nchain chains")
            else
                println("## setting the weightings/DBSCAN")
                eps= m.eps
                min_nei= m.mnei
                min_cl= m.mcl
                w3d= m.w3d
                wvel= m.wvel
                whrd= m.whrd

                mres = GaiaClustering.modelfull(eps,min_nei,min_cl,w3d,wvel,whrd)
                dfcartnorm = getDfcartnorm(dfcart, mres)
                labels = clusters(dfcartnorm.data ,eps  , 20, min_nei, min_cl)

                if length(labels) == 0
                    FLAGmcmc= 0   ## to force stop even w/o optimization
                    nchain= 0
                    println("### no solution from DBSCAN...")
                else
                    FLAGmcmc= -1
                    nchain=  m.minchainreached+1
                end
            end

            if FLAGmcmc== -1 || nchain > m.minchainreached
                if optim
                    println("## optimization completed..")
                    println("## analyzing solutions...")
                    plot_dbscanfull_mcmc(m.plotdir, "$votname.$cycle", mc , false)

                    ## get the cluster and plot it
                    println("## extracting the cluster using DBSCAN/WEIGHTING with:")
                    res= extraction_mcmc(mc, m.votname)
                    eps= res.epsm[1]
                    min_nei= trunc(Int,res.mneim[1] + 0.5)
                    min_cl= trunc(Int,res.mclm[1] + 0.5)
                    w3d= res.w3dm[1]
                    wvel= res.wvelm[1]
                    whrd= res.whrdm[1]

                    mres = GaiaClustering.modelfull(eps,min_nei,min_cl,w3d,wvel,whrd)
                    dfcartnorm = getDfcartnorm(dfcart, mres)
                    labels = clusters(dfcartnorm.data ,eps  , 20, min_nei, min_cl)
                end

                labelmax , nmax, qc = find_cluster_label2(labels, df, dfcart, m)
                println("## label $labelmax written as an oc solution...")
                export_df("$votname.$cycle", m.ocdir, df , dfcart , labels , labelmax)

                ## Principal components
                pc, pcres= compute_PC(df, dfcart, labels, labelmax)

                edgeratio1, edgeratio2= edge_ratio(dfcart, labels[labelmax])
                scproperties = get_properties_SC2(labels[labelmax] , df, dfcart)
                scdf= convertStruct2Df(scproperties)
                insertcols!(scdf, 1, :votname => votname)
                # s=size(scdf)
                insertcols!(scdf, 2, :uuid => string(m.uuid))
                insertcols!(scdf, 3, :cycle => cycle)
                insertcols!(scdf, 4, :pc3 => pcres[3])
                insertcols!(scdf, 4, :pc2 => pcres[2])
                insertcols!(scdf, 4, :pc1 => pcres[1])

                ## Xg, Yg, Zg median Galactic position
                dg = df.data[3,labels[labelmax]]
                lg = df.data[1,labels[labelmax]]
                bg = df.data[2,labels[labelmax]]
                xg= [] ; yg= [] ; zg= []
                for i in 1:size(dg)[1]
                    x1 = galXYZ(lg[i],bg[i],dg[i])
                    push!(xg,x1[1]) ; push!(yg,x1[2]) ; push!(zg,x1[3])
                end
                Xgm= median(xg) ; Ygm= median(yg); Zgm= median(zg)

                if optim
                    insertcols!(res, 2,  :cycle => cycle)
                    push!(mcmclist, res)
                    ## create DF chain
                    dfchain= create_DFchain(mc, votname, cycle)
                    push!(chainlist,dfchain)
                end

                push!(sclist, scdf)

                println("###")
                println("### solution label: $labelmax")
                print("### "); println(red(@sprintf("PC1: %3.1f , PC2: %3.1f , PC3: %3.1f", pcres[1], pcres[2], pcres[3])))
                println("### Offdeg: $(scproperties.offdeg)")
                println("### Edge ratio: $(scproperties.edgratm)")
                println("### N stars: $nmax")
                println("### Qc: $qc")
                println("###")

                ## score only makes sense for optimization, otherwise set to 1.
                if optim k= score_cycle(qc, nmax, nchain, iter) else k= 1. end
                @printf("## score cycle %d: %3.3f \n",cycle, k)

                extraplot= DataFrame(cycle=cycle, score_cycle=k, qc=qc, votname=votname, pc1=pcres[1],pc2=pcres[2], pc3=pcres[3], xg=Xgm, yg=Ygm,zg=Zgm, uuid= m.uuid)

                plot_cluster2(m.plotdir, "$votname.$cycle", labels[labelmax], scproperties,
                    dfcart , false, extraplot)

                jump= 50  # how many stars to jump in the plot
                plot_rawdata(m.plotdir, "$votname.$cycle", labels[labelmax], scproperties,
                    dfcart , pc, jump, false, extraplot)

                println("###")
                println("### subtracting BEST solution from Df...")
                dfnew, dfcartnew= remove_stars(df, dfcart, labels[labelmax])
                df= dfnew
                dfcart= dfcartnew

                ########################### STOP conditions #########
                FLAG= 0
                if nmax < m.minstarstop
                    FLAG= FLAG | (1<<0)
                    println("### extraction stopped at cycle $cycle")
                    println("### nmax too low...")
                    cyclerun= false
                end
                if cycle == m.cyclemax
                    FLAG= FLAG | (1<<1)
                    println("### extraction stopped at cycle $cycle")
                    println("### cyclemax reached...")
                    cyclerun= false
                end
                if qc < m.qcminstop
                    FLAG= FLAG | (1<<2)
                    println("### extraction stopped at cycle $cycle")
                    println("### Qc too low...")
                    cyclerun= false
                end
                if w3d/wvel < m.wratiominstop || wvel/w3d < m.wratiominstop
                    FLAG= FLAG | (1<<3)
                    println("### extraction stopped at cycle $cycle")
                    println("### weight ratio too low...")
                    cyclerun= false
                end
                if FLAGmcmc == 3 && nchain > m.minchainreached
                    FLAG= FLAG | (1<<4)
                    println("## extraction stopped at cycle $cycle")
                    println("## chain iteration not performed completely but sufficient to keep...")
                    cyclerun= false
                end
                #####################################################
                ## Time
                ##
                tend= now()
                duration= Dates.value(tend-tstart) / (1000*1)
                nstar= size(df.raw)[2]
                if optim timeperiterstar= duration / (iter*nstar) end
                if optim timeperchainstar= duration / (nchain*nstar) end
                @printf("## \n")
                @printf("## Time: \n")
                @printf("## duration per cycle %3.3f sec \n", duration)
                if optim @printf("## duration per iteration*star %3.3e sec \n", timeperiterstar) end
                if optim @printf("## duration per chain*star %3.3e sec \n", timeperchainstar) end
                @printf("##\n")

                if optim
                    ## log the results of performances
                    dfout= DataFrame(votname=votname, cycle=cycle, nstar=nstar, qc=qc, nmax=nmax, nchain=nchain, iter=iter,
                    scorecycle=k, duration=duration, timeperiterstar=timeperiterstar ,
                    timeperchainstar= timeperchainstar )
                    push!(perflist, dfout)
                end

                cycle += 1
            else
                println("## nothing found, stopped...")
                FLAG= 0
                cyclerun= false
            end
        end
        if cycle >= 2
            save_cycle_optim(sclist, mcmclist, perflist, chainlist, m, optim)
        end
        return(cycle-1, FLAG)
    end
end
##
### save results cycle in csv
###
function save_cycle(sc, mcmc, perf, chain,  m::GaiaClustering.meta)
    filesc= @sprintf("%s.sc.csv", m.prefile)
    filemcmc= @sprintf("%s.mcmc.csv", m.prefile)
    fileperf= @sprintf("%s.perf.csv", m.prefile)
    filechain= @sprintf("%s.chain.csv", m.prefile)
    votname= m.votname

    ncycle= length(sc)

    for i in 1:ncycle
        if !isfile(filesc)
            CSV.write(filesc,sc[i],delim=';')
            println("## $filesc created...")
            CSV.write(filemcmc,mcmc[i],delim=';')
            println("## $filemcmc created...")
            CSV.write(fileperf,perf[i],delim=';')
            println("## $fileperf created...")
            CSV.write(filechain,chain[i],delim=';')
            println("## $filechain created...")
        else
            res= CSV.File(filesc, delim=";") |> DataFrame
            append!(res,sc[i]) ; CSV.write(filesc,res,delim=';')
            res= CSV.File(filemcmc, delim=";") |> DataFrame
            append!(res,mcmc[i]) ; CSV.write(filemcmc,res,delim=';')
            res= CSV.File(fileperf, delim=";") |> DataFrame
            append!(res,perf[i]) ; CSV.write(fileperf,res,delim=';')
            # res= CSV.File(filechain, delim=";") |> DataFrame
            CSV.write(filechain,chain[i],delim=';', append=true)
        end
    end
end
##
### save results cycle in csv with the optim option
###
function save_cycle_optim(sc, mcmc, perf, chain,  m::GaiaClustering.meta,optim)
    filesc= @sprintf("%s.sc.csv", m.prefile)
    if optim
        filemcmc= @sprintf("%s.mcmc.csv", m.prefile)
        fileperf= @sprintf("%s.perf.csv", m.prefile)
        filechain= @sprintf("%s.chain.csv", m.prefile)
    end
    votname= m.votname
    ncycle= length(sc)

    for i in 1:ncycle
        if !isfile(filesc)
            CSV.write(filesc,sc[i],delim=';')
            println("## $filesc created...")
            if optim
                CSV.write(filemcmc,mcmc[i],delim=';')
                println("## $filemcmc created...")
                CSV.write(fileperf,perf[i],delim=';')
                println("## $fileperf created...")
                CSV.write(filechain,chain[i],delim=';')
                println("## $filechain created...")
            end
        else
            println("## cycle $i results appended to $filesc ...")
            res= CSV.File(filesc, delim=";") |> DataFrame
            append!(res,sc[i]) ; CSV.write(filesc,res,delim=';')
            if optim
                if !isfile(filemcmc)
                    CSV.write(filemcmc,mcmc[i],delim=';')
                    println("## $filemcmc created...")
                    CSV.write(fileperf,perf[i],delim=';')
                    println("## $fileperf created...")
                    CSV.write(filechain,chain[i],delim=';')
                    println("## $filechain created...")
                else
                    res= CSV.File(filemcmc, delim=";") |> DataFrame
                    append!(res,mcmc[i]) ; CSV.write(filemcmc,res,delim=';')
                    res= CSV.File(fileperf, delim=";") |> DataFrame
                    append!(res,perf[i]) ; CSV.write(fileperf,res,delim=';')
                    # res= CSV.File(filechain, delim=";") |> DataFrame
                    CSV.write(filechain,chain[i],delim=';', append=true)
                end
            end
        end
    end
end
## compute the closeness of the cluster solution to the edge of the data.
## return a ratio in the range [0., 1]. 1 means at the edge and 0 in the center
function edge_ratio(dfcart::GaiaClustering.Df, ind)
    r2= dfcart.data[2,ind] .* dfcart.data[2,ind] .+ dfcart.data[3,ind] .* dfcart.data[3,ind]
    minX= minimum(dfcart.data[1,ind])
    maxX= maximum(dfcart.data[1,ind])
    xg= median(dfcart.data[2,ind]) ; yg= median(dfcart.data[3,ind]) ; dg= median(dfcart.data[1,ind])

    indx= (dfcart.data[1,:] .<= maxX) .& (dfcart.data[1,:] .>= minX)
    rtot2=  dfcart.data[2,indx] .* dfcart.data[2,indx]  .+ dfcart.data[3,indx] .* dfcart.data[3,indx]

    ratio= sqrt((xg^2+yg^2) / maximum(rtot2))

    alpha= maximum(sqrt.(rtot2)) / maxX
    rg= alpha*dg
    ratio_2= sqrt(maximum(r2)) / rg

    return(ratio, ratio_2)
end

## compute principal components of a salution
## compute cumulated ratio for the first 3 PCs and the 3 first PCs of the normalized data
##
function compute_PC(df::GaiaClustering.Df, dfcart::GaiaClustering.Df, labels, labelmax)
        print("### Computing principal components... \n")
        s=size(labels[labelmax])
        data= zeros(8,s[1])

        X= dfcart.data[1, labels[labelmax]]
        Y= dfcart.data[2, labels[labelmax]]
        Z= dfcart.data[3, labels[labelmax]]
        vl= df.data[4,labels[labelmax]]
        vb= df.data[5,labels[labelmax]]
        gbar= df.raw[10,labels[labelmax]]
        rp= df.raw[11,labels[labelmax]]
        bp= df.raw[12,labels[labelmax]]

        data[1,:]= X
        data[2,:]= Y
        data[3,:]= Z
        data[4,:]= vl
        data[5,:]= vb
        data[6,:]= gbar
        data[7,:]= gbar .- rp
        data[8,:]= bp .- gbar

        # d=Array(data')
        dt= StatsBase.fit(ZScoreTransform, data, dims=2)
        d2= StatsBase.transform(dt, data)
        M = fit(PCA, d2)
        Yt = MultivariateStats.transform(M, d2)

        totvar= tvar(M)
        pvs= principalvars(M)
        ratioac= accumulate(+, pvs ./ totvar)

        if length(ratioac) >= 3
            pcres= 100 .* ratioac[1:3]
        else
            pcres= [100,100,100]
        end

        return(Yt, pcres)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
using Test, AFK

_eltype = Float32
_test_atype = (AFK.gpu() >= 0 ? AFK.KnetArray{_eltype} : Array{_eltype})

_INPUTdim = 10; _OUTPUTdim = 2;

@time include("primitive.jl")


text:
# ----- Outer Constructors -------------------------------------------------- #

# The following provides an outer constructor whose argument signature matches
# that of the inner constructor provided in typedefs.jl: constructs a NullableArray
# from an Array of values and an Array{Bool} mask.
function NullableArray{T, N}(A::AbstractArray{T, N},
                             m::Array{Bool, N}) # -> NullableArray{T, N}
    return NullableArray{T, N}(A, m)
end

# TODO: Uncomment this doc entry when Base Julia can parse it correctly.
# @doc """
# Allow users to construct a quasi-uninitialized `NullableArray` object by
# specifing:
#
# * `T`: The type of its elements.
# * `dims`: The size of the resulting `NullableArray`.
#
# NOTE: The `values` field will be truly uninitialized, but the `isnull` field
# will be initialized to `true` everywhere, making every entry of a new
# `NullableArray` a null value by default.
# """ ->
function NullableArray{T}(::Type{T}, dims::Dims) # -> NullableArray{T, N}
    return NullableArray(Array(T, dims), fill(true, dims))
end

# Constructs an empty NullableArray of type parameter T and number of dimensions
# equal to the number of arguments given in 'dims...', where the latter are
# dimension lengths.
function NullableArray(T::Type, dims::Int...) # -> NullableArray
    return NullableArray(T, dims)
end

@compat (::Type{NullableArray{T}}){T}(dims::Dims) = NullableArray(T, dims)
@compat (::Type{NullableArray{T}}){T}(dims::Int...) = NullableArray(T, dims)

# The following method constructs a NullableArray from an Array{Any} argument
# 'A' that contains some placeholder of type 'T' for null values.
#
# e.g.: julia> NullableArray([1, nothing, 2], Int, Void)
#       3-element NullableArrays.NullableArray{Int64,1}:
#       Nullable(1)
#       Nullable{Int64}()
#       Nullable(2)
#
#       julia> NullableArray([1, "notdefined", 2], Int, ASCIIString)
#       3-element NullableArrays.NullableArray{Int64,1}:
#       Nullable(1)
#       Nullable{Int64}()
#       Nullable(2)
#
# TODO: think about dispatching on T = Any in method above to call
# the following method passing 'T=Void' for pseudo-literal
# NullableArray construction
function NullableArray{T, U}(A::AbstractArray,
                             ::Type{T}, ::Type{U}) # -> NullableArray{T, N}
    res = NullableArray(T, size(A))
    for i in 1:length(A)
        if !isa(A[i], U)
            @inbounds setindex!(res, A[i], i)
        end
    end
    return res
end

# The following method constructs a NullableArray from an Array{Any} argument
# `A` that contains some placeholder value `na` for null values.
#
# e.g.: julia> NullableArray(Any[1, "na", 2], Int, "na")
#       3-element NullableArrays.NullableArray{Int64,1}:
#       Nullable(1)
#       Nullable{Int64}()
#       Nullable(2)
#
function NullableArray{T}(A::AbstractArray,
                             ::Type{T},
                             na::Any;
                             conversion::Base.Callable=Base.convert) # -> NullableArray{T, N}
    res = NullableArray(T, size(A))
    for i in 1:length(A)
        if !isequal(A[i], na)
            @inbounds setindex!(res, A[i], i)
        end
    end
    return res
end

# The following method allows for the construction of zero-element
# NullableArrays by calling the parametrized type on zero arguments.
@compat (::Type{NullableArray{T, N}}){T, N}() = NullableArray(T, ntuple(i->0, N))

# ----- Conversion to NullableArrays ---------------------------------------- #
# Also provides constructors from arrays via the fallback mechanism.

#----- Conversion from arrays (of non-Nullables) -----------------------------#
function Base.convert{S, T, N}(::Type{NullableArray{T, N}},
                               A::AbstractArray{S, N}) # -> NullableArray{T, N}
    NullableArray{T, N}(convert(Array{T, N}, A), fill(false, size(A)))
end

function Base.convert{S, T, N}(::Type{NullableArray{T}},
                               A::AbstractArray{S, N}) # -> NullableArray{T, N}
    convert(NullableArray{T, N}, A)
end

function Base.convert{T, N}(::Type{NullableArray},
                            A::AbstractArray{T, N}) # -> NullableArray{T, N}
    convert(NullableArray{T, N}, A)
end

#----- Conversion from arrays of Nullables -----------------------------------#
function Base.convert{S, T, N}(::Type{NullableArray{T, N}},
                               A::AbstractArray{Nullable{S}, N}) # -> NullableArray{T, N}
   out = NullableArray{T, N}(Array(T, size(A)), Array(Bool, size(A)))
   for i = 1:length(A)
       if !(out.isnull[i] = isnull(A[i]))
           out.values[i] = A[i].value
       end
   end
   out
end

#----- Conversion from NullableArrays of a different type --------------------#
Base.convert{T, N}(::Type{NullableArray}, X::NullableArray{T,N}) = X

function Base.convert{S, T, N}(::Type{NullableArray{T}},
                               A::AbstractArray{Nullable{S}, N}) # -> NullableArray{T, N}
    convert(NullableArray{T, N}, A)
end

function Base.convert{T, N}(::Type{NullableArray},
                            A::AbstractArray{Nullable{T}, N}) # -> NullableArray{T, N}
    convert(NullableArray{T, N}, A)
end

function Base.convert{S, T, N}(::Type{NullableArray{T, N}},
                               A::NullableArray{S, N}) # -> NullableArray{T, N}
    NullableArray(convert(Array{T, N}, A.values), A.isnull)
end


text:
#AnalogMTK: Base functionnality and types
#-------------------------------------------------------------------------------

#==
===============================================================================#
struct Ideal
end

#Last line


text:
function ensure_downloaded_files()
    global download_cache
    download_cache::String
    if isempty(readdir(download_cache))
        force_download_files(download_cache)
    end
    return download_cache::String
end

function force_download_files(dest_directory::AbstractString)
    let url = "https://www.hcup-us.ahrq.gov/toolssoftware/ccs/AppendixASingleDX.txt"
        @info "Downloading files for CCS9"
        sleep(3) # Sleep for a bit, so that we can be compliant with any rate limits
        src = Downloads.download(url)
        dest_filename = "ccs9.txt"
        dest = joinpath(
            dest_directory,
            dest_filename,
        )
        mv(src, dest; force = true)
    end
    let url = "https://www.hcup-us.ahrq.gov/toolssoftware/ccsr/DXCCSR_v2022-1.zip"
        @info "Downloading files for CCS10"
        sleep(3) # Sleep for a bit, so that we can be compliant with any rate limits
        downloaded_filename = Downloads.download(url)
        dest_filename = "ccs10.csv"
        dest = joinpath(
            dest_directory,
            dest_filename,
        )
        _extract_zip_file(downloaded_filename) do extraction_output_dir
            src = joinpath(extraction_output_dir, "DXCCSR_v2022-1.CSV")
            mv(src, dest)
        end
    end
    return nothing
end

function _extract_zip_file(
        f::F,
        zip_file::AbstractString,
    ) where {F <: Function}
    mktempdir() do extraction_output_dir
        cd(extraction_output_dir) do
            extraction_cmd = `$(Pkg.PlatformEngines.exe7z())`
            push!(extraction_cmd.exec, "x")
            push!(extraction_cmd.exec, "$(zip_file)")
            push!(extraction_cmd.exec, "-o$(extraction_output_dir)")
            run(extraction_cmd)
            f(extraction_output_dir)
        end
    end
    return nothing
end


text:
using Test
using Transpiler
const T = Transpiler
const root_dir = dirname(dirname(pathof(Transpiler)))
is_project_file(x) = any(endswith(x,y) for y in (".py",".pxd","pyx"))

@testset "SCAMPy" begin
    foreign = joinpath(root_dir, "..", "CliMA", "SCAMPyRepos", "SCAMPyOriginal", "SCAMPy")
    example_root=joinpath(root_dir,"examples","SCAMPy")

    paths = Dict(
        :foreign=>foreign,
        :local_original=>joinpath(example_root, "original"),
        :example_root=>example_root,
        )

    project = T.configure_project(;
        paths=paths,
        file_filter=is_project_file,
        folder_filter=x->!occursin("__pycache__", x),
        skip_import=true,
        language_to=Julia(),
        language_from=Cython()
        )

    T.transpile_py_2_jl(project)
end


text:
using GR

function points_from_image(img, npts)
  w, h = size(img)
  xpts = Float64[]
  ypts = Float64[]
  cols = Int[]
  for i in 1:npts
    x = rand() * w
    y = rand() * h
    c = img[Int64(floor(x)) + 1, Int64(floor(h-y)) + 1]
    r = ( c        & 0xff) / 255.0
    g = ((c >> 8 ) & 0xff) / 255.0
    b = ((c >> 16) & 0xff) / 255.0
    if 0.2989 * r + 0.5870 * g + 0.1140 * b > 0.8
      if rand() < 0.1
        push!(xpts, x)
        push!(ypts, y)
        push!(cols, 1)
      end
      continue
    end
    push!(xpts, x)
    push!(ypts, y)
    push!(cols, inqcolorfromrgb(r, g, b))
  end
  xpts, ypts, cols
end

w, h, img = readimage("julia_logo.png")
x, y, cols = points_from_image(img, 30000)

settransparency(0.5)
setmarkersize(0.5)
scatter(x, y, size=(900,600))

n, tri = delaunay(x, y)
for i in 1:n
  if all(cols[tri[i,:]] .!= 1)
    settransparency(0.5)
    setfillcolorind(cols[tri[i,1]])
  else
    settransparency(0.1)
    setfillcolorind(1)
  end
  fillarea(x[tri[i,:]], y[tri[i,:]])
end
updatews()