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

struct StorageAccountsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Checks that the storage account name is valid and is not already in use.
Param: accountName::StorageAccountCheckNameAvailabilityParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: CheckNameAvailabilityResult
"""
function storageAccountsCheckNameAvailability(_api::StorageAccountsApi, accountName, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", CheckNameAvailabilityResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Storage/checkNameAvailability", ["azure_auth"], accountName)
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Asynchronously creates a new storage account with the specified parameters. If an account is already created and a subsequent create request is issued with different properties, the account properties will be updated. If an account is already created and a subsequent create or update request is issued with the exact same set of properties, the request will succeed.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: parameters::StorageAccountCreateParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccount
"""
function storageAccountsCreate(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsCreate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsCreate", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsCreate", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsCreate", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "PUT", StorageAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Deletes a storage account in Microsoft Azure.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function storageAccountsDelete(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsDelete", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsDelete", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsDelete", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsDelete", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Returns the properties for the specified storage account including but not limited to name, SKU name, location, and account status. The ListKeys operation should be used to retrieve storage keys.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccount
"""
function storageAccountsGetProperties(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsGetProperties", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsGetProperties", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsGetProperties", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsGetProperties", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "GET", StorageAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all the storage accounts available under the subscription. Note that storage keys are not returned; use the ListKeys operation for this.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccountListResult
"""
function storageAccountsList(_api::StorageAccountsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", StorageAccountListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Storage/storageAccounts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

List SAS credentials of a storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: parameters::AccountSasParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ListAccountSasResponse
"""
function storageAccountsListAccountSAS(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsListAccountSAS", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsListAccountSAS", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsListAccountSAS", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsListAccountSAS", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "POST", ListAccountSasResponse, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/ListAccountSas", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists all the storage accounts available under the given resource group. Note that storage keys are not returned; use the ListKeys operation for this.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccountListResult
"""
function storageAccountsListByResourceGroup(_api::StorageAccountsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsListByResourceGroup", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsListByResourceGroup", :minLength, resourceGroupName, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", StorageAccountListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists the access keys for the specified storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccountListKeysResult
"""
function storageAccountsListKeys(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsListKeys", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsListKeys", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsListKeys", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsListKeys", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "POST", StorageAccountListKeysResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/listKeys", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

List service SAS credentials of a specific resource.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: parameters::ServiceSasParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: ListServiceSasResponse
"""
function storageAccountsListServiceSAS(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsListServiceSAS", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsListServiceSAS", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsListServiceSAS", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsListServiceSAS", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "POST", ListServiceSasResponse, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/ListServiceSas", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Regenerates one of the access keys for the specified storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: regenerateKey::StorageAccountRegenerateKeyParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccountListKeysResult
"""
function storageAccountsRegenerateKey(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, regenerateKey, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsRegenerateKey", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsRegenerateKey", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsRegenerateKey", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsRegenerateKey", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "POST", StorageAccountListKeysResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/regenerateKey", ["azure_auth"], regenerateKey)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

The update operation can be used to update the SKU, encryption, access tier, or tags for a storage account. It can also be used to map the account to a custom domain. Only one custom domain is supported per storage account; the replacement/change of custom domain is not supported. In order to replace an old custom domain, the old value must be cleared/unregistered before a new value can be set. The update of multiple properties is supported. This call does not change the storage keys for the account. If you want to change the storage account keys, use the regenerate keys operation. The location and name of the storage account cannot be changed after creation.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: parameters::StorageAccountUpdateParameters (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: StorageAccount
"""
function storageAccountsUpdate(_api::StorageAccountsApi, resourceGroupName::String, accountName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "storageAccountsUpdate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "storageAccountsUpdate", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "storageAccountsUpdate", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "storageAccountsUpdate", :minLength, accountName, 3)

    _ctx = Swagger.Ctx(_api.client, "PATCH", StorageAccount, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export storageAccountsCheckNameAvailability, storageAccountsCreate, storageAccountsDelete, storageAccountsGetProperties, storageAccountsList, storageAccountsListAccountSAS, storageAccountsListByResourceGroup, storageAccountsListKeys, storageAccountsListServiceSAS, storageAccountsRegenerateKey, storageAccountsUpdate


text:
struct IDX
    # contains the variable indexing stuff for [x;s;z;y]
    nx::Int64
    ns::Int64
    nz::Int64
    ny::Int64
    N::Int64

    x::UnitRange{Int64}
    s::UnitRange{Int64}
    z::UnitRange{Int64}
    y::UnitRange{Int64}
end

struct DELTA
    # stores all the deltas (search directions)

    # affine deltas
    x_a::Array{Float64,1}
    s_a::Array{Float64,1}
    z_a::Array{Float64,1}
    y_a::Array{Float64,1}

    # centering + correcting deltas
    x_c::Array{Float64,1}
    s_c::Array{Float64,1}
    z_c::Array{Float64,1}
    y_c::Array{Float64,1}

    # total deltas
    x::Array{Float64,1}
    s::Array{Float64,1}
    z::Array{Float64,1}
    y::Array{Float64,1}

    # constructor
    function DELTA(nx,ns,nz,ny)
        new(zeros(nx),zeros(ns),zeros(nz),zeros(ny),
            zeros(nx),zeros(ns),zeros(nz),zeros(ny),
            zeros(nx),zeros(ns),zeros(nz),zeros(ny))
    end
end

mutable struct SolverSettings
    verbose::Bool
    max_iters::Int64
    atol::Float64
    function SolverSettings()
        return new(false, 50, 1e-8)
    end
end

struct QP

    # problem data
    Q::SparseMatrixCSC{Float64,Int64}
    q::Array{Float64,1}
    A::SparseMatrixCSC{Float64,Int64}
    b::Array{Float64,1}
    G::SparseMatrixCSC{Float64,Int64}
    h::Array{Float64,1}

    # variables
    x::Array{Float64,1} # primal
    s::Array{Float64,1} # primal
    z::Array{Float64,1} # dual
    y::Array{Float64,1} # dual

    # KKT stuff
    KKT::SparseMatrixCSC{Float64,Int64}
    rhs_a::Array{Float64,1}
    rhs_c::Array{Float64,1}
    p_a::Array{Float64,1}
    p_c::Array{Float64,1}

    # indexing
    idx::IDX

    # deltas
    Δ::DELTA

    # solver settings
    opts::SolverSettings

    function QP(Q,q,A,b,G,h)

        # length of variables
        nx = length(q)
        ns = length(h)
        nz = length(h)
        ny = length(b)
        N = nx + ns + nz + ny

        # indexing when stacked [x;s;z;y]
        idx_x = 1:nx
        idx_s = (nx + 1) : (nx + ns)
        idx_z = (nx + ns + 1) : (nx + ns + nz)
        idx_y = (nx + ns + nz + 1) : (nx + ns + nz + ny)
        idx = IDX(nx,ns,nz,ny,N, idx_x, idx_s, idx_z, idx_y)

        # kkt stuff
        KKT = spzeros(N,N)
        rhs_a = zeros(N)
        rhs_c = zeros(N)
        p_a = zeros(N)
        p_c = zeros(N)

        # initialize variables to zero
        x = zeros(nx)
        s = zeros(ns)
        z = zeros(nz)
        y = zeros(ny)

        # deltas
        Δ = DELTA(nx,ns,nz,ny)

        new(Q,q,A,b,G,h,x,s,z,y, KKT, rhs_a, rhs_c, p_a, p_c, idx, Δ, SolverSettings())
    end

end

# ---------------real functions---------------

function rhs_kkt_a!(qp::QP)
    idx = qp.idx
    qp.rhs_a[idx.x] = -(qp.A'*qp.y + qp.G'*qp.z + qp.Q*qp.x + qp.q)
    qp.rhs_a[idx.s] = -(qp.s .* qp.z)
    qp.rhs_a[idx.z] = -(qp.G*qp.x + qp.s - qp.h)
    qp.rhs_a[idx.y] = -(qp.A*qp.x - qp.b)
    return nothing
end
function index_sol_a!(qp::QP)
    qp.Δ.x_a .= qp.p_a[qp.idx.x]
    qp.Δ.s_a .= qp.p_a[qp.idx.s]
    qp.Δ.z_a .= qp.p_a[qp.idx.z]
    qp.Δ.y_a .= qp.p_a[qp.idx.y]
    return nothing
end
function index_sol_c!(qp::QP)
    qp.Δ.x_c .= qp.p_c[qp.idx.x]
    qp.Δ.s_c .= qp.p_c[qp.idx.s]
    qp.Δ.z_c .= qp.p_c[qp.idx.z]
    qp.Δ.y_c .= qp.p_c[qp.idx.y]
    return nothing
end
function linesearch(x,dx)
    α = min(1.0, minimum([dx[i]<0 ? -x[i]/dx[i] : Inf for i = 1:length(x)]))
    return α
end
function centering_params(qp::QP)

    μ = dot(qp.s,qp.z)/qp.idx.ns

    α = min(linesearch(qp.s,qp.Δ.s_a), linesearch(qp.z,qp.Δ.z_a))

    σ = (dot(qp.s + α*qp.Δ.s_a, qp.z + α*qp.Δ.z_a)/dot(qp.s,qp.z))^3
    return σ, μ
end
function rhs_kkt_c!(qp::QP, σ, μ)
    idx = qp.idx
    qp.rhs_c .= 0
    qp.rhs_c[idx.s] = (σ*μ .- (qp.Δ.s_a .* qp.Δ.z_a))
    return nothing
end
function combine_deltas!(qp::QP)
    qp.Δ.x .= qp.Δ.x_a + qp.Δ.x_c
    qp.Δ.s .= qp.Δ.s_a + qp.Δ.s_c
    qp.Δ.z .= qp.Δ.z_a + qp.Δ.z_c
    qp.Δ.y .= qp.Δ.y_a + qp.Δ.y_c
    return nothing
end
function update_vars!(qp::QP,α)
    qp.x .+= α*qp.Δ.x
    qp.s .+= α*qp.Δ.s
    qp.z .+= α*qp.Δ.z
    qp.y .+= α*qp.Δ.y
    return nothing
end

struct QPMaxItersError <: Exception
    var::Symbol
end

function solveqp!(qp::QP)

    if qp.opts.verbose
        @printf "iter     objv        gap       |Ax-b|    |Gx+s-h|    step\n"
        @printf "---------------------------------------------------------\n"
    end

    initialize!(qp)

    initialize_kkt!(qp)

    for i = 1:qp.opts.max_iters

        # update linear system for solves
        update_kkt!(qp)
        # kkt_factor = qdldl(qp.KKT +
                          # 1e-10*Diagonal([ones(qp.idx.nx + qp.idx.ns);-ones(qp.idx.nz + qp.idx.ny)]))
        kkt_factor = lu(qp.KKT)

        # affine step
        rhs_kkt_a!(qp)

        qp.p_a .= kkt_factor\qp.rhs_a
        # qp.p_a .= iterative_ref(kkt_factor,qp.KKT,qp.rhs_a)
        index_sol_a!(qp)

        # centering and correcting step
        σ, μ = centering_params(qp)
        rhs_kkt_c!(qp, σ, μ)
        qp.p_c .= kkt_factor\qp.rhs_c
        # qp.p_c .= iterative_ref(kkt_factor,qp.KKT,qp.rhs_c)
        index_sol_c!(qp)

        # combine deltas
        combine_deltas!(qp)

        # last linesearch
        α = min(1,0.99*min(linesearch(qp.s,qp.Δ.s),linesearch(qp.z,qp.Δ.z)))

        update_vars!(qp,α)

        if logging(qp::QP,i,α)
            return i
            # break
        end

        # if i == qp.opts.max_iters
            # throw(QPMaxItersError)
        # end
    end
    throw(QPMaxItersError)
    # return nothing
end

function logging(qp::QP,iter,α)

    J = 0.5*qp.x'*qp.Q*qp.x + dot(qp.q,qp.x)
    gap = dot(qp.s,qp.z)/length(qp.s)
    eq_res = norm(qp.A*qp.x - qp.b,Inf)
    ineq_res = norm(qp.G*qp.x + qp.s - qp.h,Inf)

    if qp.opts.verbose
        @printf("%3d   %10.3e  %9.2e  %9.2e  %9.2e  % 6.4f\n",
              iter, J, gap, eq_res,
              ineq_res, α)
    end

    # check for convergence
    if (eq_res < 1e-10) & (ineq_res < 1e-10) & (gap<qp.opts.atol)
        return true
    else
        return false
    end
end

function initialize!(qp::QP)
    idx = qp.idx
    Ni = idx.nx + idx.nz + idx.ny
    idx_y = idx.y .- idx.ns

    A = spzeros(Ni, Ni)
    A[idx.x,idx.x] = qp.Q
    A[idx.x,idx.s] = qp.G'
    A[idx.x,idx_y] = qp.A'
    A[idx.s,idx.x] = qp.G
    A[idx.s,idx.s] = -I(idx.ns)
    A[idx_y,idx.x] = qp.A

    init = A\[-qp.q;qp.h;qp.b]

    qp.x .= init[idx.x]
    qp.z .= init[idx.s]
    qp.y .= init[idx_y]

    α_p = -minimum(-qp.z)
    if α_p < 0
        qp.s .= -qp.z
    else
        qp.s .= -qp.z .+ (1 + α_p)
    end

    α_d = -minimum(qp.z)
    if α_d >= 0
        qp.z .= qp.z .+ (1 + α_d)
    end

    return nothing
end

function initialize_kkt!(qp::QP)
    idx = qp.idx

    qp.KKT[idx.x, idx.x] = qp.Q
    qp.KKT[idx.x, idx.z] = qp.G'
    qp.KKT[idx.x, idx.y] = qp.A'
    # qp.KKT[idx.s, idx.s] = Diagonal(qp.z)
    # qp.KKT[idx.s, idx.z] = Diagonal(qp.s)
    qp.KKT[idx.z, idx.x] = qp.G
    qp.KKT[idx.z, idx.s] = I(idx.nz)
    qp.KKT[idx.y, idx.x] = qp.A

    return nothing
end

function update_kkt!(qp::QP)
    idx = qp.idx
    qp.KKT[idx.s, idx.s] = Diagonal(qp.z)#Diagonal(qp.z ./ qp.s)
    qp.KKT[idx.s, idx.z] = Diagonal(qp.s)#I(idx.ns)
    return nothing
end

function quadprog(Q,q,A,b,G,h; verbose = false, max_iters = 50, atol = 1e-8)
    qp = QP(Q,q,A,b,G,h)
    qp.opts.verbose = verbose
    qp.opts.max_iters = max_iters
    qp.opts.atol = atol
    qp_iters = solveqp!(qp::QP)
    return qp.x, qp_iters
end

ext: jl


text:
#####################
# Probabilistic circuit mixtures
#####################

"A probabilistic mixture model"
abstract type AbstractMixture end

"A probabilistic mixture model whose components are not themselves mixtures"
abstract type AbstractFlatMixture <: AbstractMixture end

"A probabilistic mixture model whose components are themselves mixtures"
abstract type AbstractMetaMixture <: AbstractMixture end

"A probabilistic mixture model of probabilistic circuits"
struct FlatMixture <: AbstractFlatMixture
    weights::Vector{Float64}
    components::Vector{<:ProbΔ}
    FlatMixture(w,c) = begin
        @assert length(w) == length(c)
        @assert sum(w) ≈ 1.0
        new(w,c)
    end
end

FlatMixture(c) = FlatMixture(uniform(length(c)),c)

"A mixture with cached flow circuits for each component (which are assumed to be ProbΔs)"
struct FlatMixtureWithFlow <: AbstractFlatMixture
    origin::FlatMixture
    flowcircuits::Vector{<:FlowΔ}
    FlatMixtureWithFlow(origin,fcs) = begin
        @assert num_components(origin) == length(fcs)
        foreach(components(origin), fcs) do or, fc
            @assert or[end] === prob_origin(fc)[end]
        end
        new(origin,fcs)
    end
end

FlatMixtureWithFlow(w,c,f) = FlatMixtureWithFlow(FlatMixture(w,c),f)

"A probabilistic mixture model of mixture models"
struct MetaMixture <: AbstractMetaMixture
    weights::Vector{Float64}
    components::Vector{<:AbstractMixture}
    MetaMixture(w,c) = begin
        @assert length(w) == length(c)
        @assert sum(w) ≈ 1.0
        new(w,c)
    end
end

MetaMixture(c) = MetaMixture(uniform(length(c)),c)

Mixture(w, c::Vector{<:AbstractMixture}) = MetaMixture(w, c)
Mixture(w, c::Vector{<:ProbΔ}) = FlatMixture(w, c)

#####################
# Functions
#####################

"Get the components in this mixture"
@inline components(m::FlatMixture) = m.components
@inline components(m::FlatMixtureWithFlow) = components(m.origin)
@inline components(m::MetaMixture) = m.components

"Get the component weights in this mixture"
@inline component_weights(m::FlatMixture) = m.weights
@inline component_weights(m::FlatMixtureWithFlow) = component_weights(m.origin)
@inline component_weights(m::MetaMixture) = m.weights

"Number of components in a mixture"
@inline num_components(m::AbstractMixture)::Int = length(components(m))

"Convert a given flat mixture into one with cached flows"
ensure_with_flows(m::FlatMixture, size_hint::Int)::FlatMixtureWithFlow = begin
    flowcircuits = [FlowΔ(pc, size_hint, Bool, opts_accumulate_flows) for pc in components(m)]
    FlatMixtureWithFlow(m,flowcircuits)
end
ensure_with_flows(m::FlatMixtureWithFlow, ::Int)::FlatMixtureWithFlow = m

replace_prob_circuits(m::FlatMixture, pcs::Vector{ProbΔ}) =
    FlatMixture(component_weights(m), pcs)

# log_likelihood

function log_likelihood(mixture::FlatMixture, batches::XBatches{Bool})::Float64
    mwf = ensure_with_flows(mixture, max_batch_size(batches))
    log_likelihood(mwf, batches)
end

function log_likelihood(mixture::FlatMixtureWithFlow, batches::XBatches{Bool})::Float64
    # assume the per-batch call will compute a weighted sum over examples
    sum(batch -> log_likelihood(mixture, batch), batches)
end

function log_likelihood(mixture::FlatMixtureWithFlow, batch::PlainXData{Bool})::Float64
    sum(log_likelihood_per_instance(mixture, batch))
end

function log_likelihood(mixture::MetaMixture, batches::XBatches{Bool})::Float64
    sum(batch -> log_likelihood(mixture, batch), batches)
end

function log_likelihood(mixture::MetaMixture, batch::PlainXData{Bool})::Float64
    sum(log_likelihood_per_instance(mixture, batch))
end

# log_likelihood_per_instance (including mixture weight likelihood)

function log_likelihood_per_instance(mixture::FlatMixture, batches::XBatches{Bool})::Vector{Float64}
    mwf = ensure_with_flows(mixture, max_batch_size(batches))
    log_likelihood_per_instance(mwf, batches)
end

function log_likelihood_per_instance(mixture::FlatMixtureWithFlow, batches::XBatches{Bool})::Vector{Float64}
    mapreduce(b -> log_likelihood_per_instance(mixture, b), vcat, batches)
end

function log_likelihood_per_instance(mixture::FlatMixtureWithFlow, batch::PlainXData{Bool})::Vector{Float64}
    log_p_of_x_and_c = log_likelihood_per_instance_component(mixture, batch)
    logsumexp(log_p_of_x_and_c, 2)
end

function log_likelihood_per_instance(mixture::MetaMixture, batches::XBatches{Bool})::Vector{Float64}
    mapreduce(b -> log_likelihood_per_instance(mixture, b), vcat, batches)
end

function log_likelihood_per_instance(mixture::MetaMixture, batches::PlainXData{Bool})::Vector{Float64}
    log_p_of_x_and_c = log_likelihood_per_instance_component(mixture, batch)
    logsumexp(log_p_of_x_and_c, 2)
end

# Log likelihoods per instance and component (including mixture weight likelihood)

"Log likelihood per instance and component. A vector of matrices per batch where the first dimension is instance, second is components."
function log_likelihood_per_instance_component(mixture::FlatMixtureWithFlow, batches::XBatches{Bool})::Vector{Matrix{Float64}}
    [log_likelihood_per_instance_component(mixture, batch) for batch in batches]
end

"Log likelihood per instance and component. First dimension is instance, second is components."
function log_likelihood_per_instance_component(mixture::FlatMixtureWithFlow, batch::PlainXData{Bool})::Matrix{Float64}
    hcat(log_likelihood_per_component_instance(mixture, batch)...)
end

"Log likelihood per component and instance. Outer vector is components, inner vector is instances"
function log_likelihood_per_component_instance(mixture::FlatMixtureWithFlow, batch::PlainXData{Bool})::Vector{Vector{Float64}}
    map(mixture.flowcircuits, component_weights(mixture)) do fc, component_weight
        log_likelihood_per_instance(fc, batch) .+ log(component_weight)
    end
end

"Log likelihood per instance and component. A vector of matrices per batch where the first dimension is instance, second is components."
function log_likelihood_per_instance_component(mixture::MetaMixture, batches::XBatches{Bool})::Vector{Matrix{Float64}}
    [log_likelihood_per_instance_component(mixture, batch) for batch in batches]
end

"Log likelihood per instance and component. First dimension is instance, second is components."
function log_likelihood_per_instance_component(mixture::MetaMixture, batch::PlainXData{Bool})::Matrix{Float64}
    hcat(log_likelihood_per_component_instance(mixture, batch)...)
end

"Log likelihood per component and instance. Outer vector is components, inner vector is instances"
function log_likelihood_per_component_instance(mixture::MetaMixture, batch::PlainXData{Bool})::Vector{Vector{Float64}}
    map(mixture.components, component_weights(mixture)) do c, component_weight
        log_likelihood_per_instance(c, batch) .+ log(component_weight)
    end
end


text:
#=
Created on 09/01/2021 12:41:34
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Inputs a matrix and a list with headers and prints a table in markdown
format. Automatically copies this to the clipboard using the `clipboard` function/
=#

table = rand(1:10, 8, 4)
header = ["A", "B", "C", "D"]

function markdowntable(table, header)
    n, m = size(table)
    table_string = ""
    # add header
    table_string *= "| " * join(header, " | ") * " |\n"
    # horizontal lines for separating the columns
    table_string *= "| " * repeat(":--|", m) * "\n"
    for i in 1:n
        table_string *= "| " * join(table[i,:], " | ") * " |\n"
    end
    clipboard(table_string)
    return table_string
end

markdowntable(table, header)


text:
"""
    is_leap_year(year)

Return `true` if `year` is a leap year in the gregorian calendar.

"""
function is_leap_year(year)
    rem(year, 4) == 0 &&
        (rem(year, 4) > 0 && rem(year, 4) > 0)
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
using Weave

weave(
    "goma.jmd",
    out_path=:pwd,
    doctype="md2html",
)

weave(
    "goma.jmd",
    out_path=:pwd,
    doctype="md2pdf",
    template="julia_tex.tpl",
)

#=
weave(
    "goma.jmd",
    out_path=:pwd,
    doctype="pandoc"
)
=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using Test, SymbolicUtils
using SymbolicUtils.Code
using SymbolicUtils.Code: LazyState
using StaticArrays
using LabelledArrays
using SparseArrays

test_repr(a, b) = @test repr(Base.remove_linenums!(a)) == repr(Base.remove_linenums!(b))

@testset "Code" begin
    @syms a b c d e p q t x(t) y(t) z(t)
    @test toexpr(Assignment(a, b)) == :(a = b)
    @test toexpr(a ← b) == :(a = b)
    @test toexpr(a+b) == :($(+)(a, b))
    @test toexpr(a+b) == :($(+)(a, b))
    @test toexpr(a^b) == :($(^)(a, b))
    @test toexpr(a^2) == :($(^)(a, 2))
    @test toexpr(a^-2) == :($(^)($(inv)(a), 2))
    @test toexpr(x(t)+y(t)) == :($(+)(x(t), y(t)))
    @test toexpr(x(t)+y(t)+x(t+1)) == :($(+)(x(t), y(t), x($(+)(1, t))))
    s = LazyState()
    push!(s.symbolify, x(t))
    push!(s.symbolify, y(t))
    @test toexpr(x(t)+y(t)+x(t+1), s) == :($(+)(var"x(t)", var"y(t)", x($(+)(1, t))))

    ex = :(let a = 3, b = $(+)(1,a)
               $(+)(a, b)
           end)

    test_repr(toexpr(Let([a ← 3, b ← 1+a], a + b)), ex)

    test_repr(toexpr(Func([],[], a+b)) |>Base.remove_linenums!, :(function ()
            $(+)(a, b)
        end))

    test_repr(toexpr(Func([x(t), x],[b ← a+2, y(t) ← b], x(t)+x(t+1)+b+y(t))),
              :(function (var"x(t)", x; b = $(+)(2, a), var"y(t)" = b)
                    $(+)(b, var"x(t)", var"y(t)", x($(+)(1, t)))
                end))
    test_repr(toexpr(Func([DestructuredArgs([x, x(t)], :state),
                           DestructuredArgs((a, b), :params)], [],
                          x(t+1) + x(t) + a  + b)),
              :(function (state, params)
                    let x = state[1], var"x(t)" = state[2], a = params[1], b = params[2]
                        $(+)(a, b, var"x(t)", x($(+)(1, t)))
                    end
                end))
    test_repr(toexpr(SetArray(false, a, [x(t), AtIndex(9, b), c])),
              quote
                  a[1] = x(t)
                  a[9] = b
                  a[3] = c
                  nothing
              end)
    @test toexpr(SetArray(true, a, [x(t), AtIndex(9, b), c])).head == :macrocall

    test_repr(toexpr(LiteralExpr(:(let x=1, y=2
                                       $(sin(a+b))
                                   end))),
              :(let x = 1, y = 2
                    $(sin)($(+)(a, b))
                end))

    test_repr(toexpr(MakeArray([a,b,a+b], :arr)),
              quote
                  $(SymbolicUtils.Code.create_array)(typeof(arr), nothing, Val{(3,)}(), a, b, $(+)(a, b))
              end)

    toexpr(Let([a ← 1, b ← 2, :arr ← [1,2]],
               MakeArray([a,b,a+b,a/b], :arr)))

    test_repr(toexpr(Let([DestructuredArgs([x(t),b,c], :foo) ← [3,3,[1,4]],
                          DestructuredArgs([p,q], c)],
                         x(t)+a+b+c)),
              :(let foo = Any[3, 3, [1, 4]],
                    var"x(t)" = foo[1], b = foo[2], c = foo[3],
                    p = c[1], q = c[2]
                    $(+)(a, b, c, var"x(t)")
                end))

    test_repr(toexpr(Func([DestructuredArgs([a,b],c,inds=[:a, :b])], [],
                          a + b)),
              :(function (c,)
                    let a = c.a, b = c.b
                        $(+)(a, b)
                    end
                end))
    @syms arr

    @test eval(toexpr(Let([a ← 1, b ← 2, arr ← [1,2]],
                          MakeArray([a,b,a+b,a/b], arr)))) == [1, 2, 3, 1/2]

    @test eval(toexpr(Let([a ← 1, b ← 2, arr ← [1,2]],
                          MakeArray([a b;a+b a/b], arr)))) == [1 2; 3 1/2]

    @test eval(toexpr(Let([a ← 1, b ← 2, arr ← @SVector([1,2])],
                          MakeArray([a,b,a+b,a/b], arr)))) === @SVector [1, 2, 3, 1/2]

    @test eval(toexpr(Let([a ← 1, b ← 2, arr ← @SVector([1,2])],
                          MakeArray([a b;a+b a/b], arr)))) === @SMatrix [1 2; 3 1/2]

    @test eval(toexpr(Let([a ← 1, b ← 2, arr ← @SLVector((:a, :b))(@SVector[1,2])],
                          MakeArray([a+b,a/b], arr)))) === @SLVector((:a, :b))(@SVector [3, 1/2])

    mksp = MakeSparseArray(sparse([1,2,31,32,2],
                                  [1,2,31,32,2],
                                  [a, b, a+b, a/b, a-b+e]))
    reference = sparse([1,2,31,32],
                     [1,2,31,32],
                     [a, a+e, a+b, a/b])

    test_repr(mksp.array, reference)

    test_repr(toexpr(mksp),
              :(SparseMatrixCSC(32, 32,
                                $(reference.colptr),
                                $(reference.rowval),
                                [$(map(toexpr, reference.nzval)...)])))

    spvec = sparsevec([5], [a], 10)

    test_repr(toexpr(MakeSparseArray(spvec)),
              :(SparseVector(10, $(spvec.nzind), [a])))

    test_repr(toexpr(MakeTuple((a, b, a+b))),
              :((a,b,$(+)(a,b))))
end


text:
using Test, PolyChaos

@test nw(EmptyQuad()) |> typeof == Array{Float64,2}
@test nw(EmptyQuad()) |> length == 0

nodes, weights = gauss(rm_legendre01(5)...)
myQuad = Quad("myQuad", length(nodes), nodes, weights)
foo(x) = 6*x^5
res = 1
@test isapprox(integrate(foo, nodes, weights), res)
@test isapprox(integrate(foo, myQuad), res)
@test_throws DomainError integrate(foo, EmptyQuad())
@test isapprox(integrate(foo, Uniform01OrthoPoly(4)), res)

d, nunc = 5, 10
op = genLaguerreOrthoPoly(d,1.23,addQuadrature = false)
opq = genLaguerreOrthoPoly(d,1.23)
mop = MultiOrthoPoly([op for i in 1:nunc],d)
mopq = MultiOrthoPoly([opq for i in 1:nunc],d)
@testset "dimenions" begin
    @test isequal(PolyChaos.dim(op), d + 1)
    @test isequal(PolyChaos.dim(opq), d + 1)
    @test isequal(PolyChaos.dim(mop), factorial(d + nunc) / (factorial(d) * factorial(nunc)))
end

@testset "degrees" begin
    @test isequal(PolyChaos.deg(op),d)
    @test isequal(PolyChaos.deg(opq),d)
    @test isequal(PolyChaos.deg(mop),d)
end

NW = [ opq.quad.nodes opq.quad.weights ]
@testset "nodes and weights" begin
    @test PolyChaos.nw(opq) == NW
    @test PolyChaos.nw([opq for i in 1:nunc])[1] == [NW[:,1] for i in 1:nunc]
    @test PolyChaos.nw([opq for i in 1:nunc])[2] == [NW[:,2] for i in 1:nunc]
    @test PolyChaos.nw(mopq)[1] == [NW[:,1] for i in 1:nunc]
    @test PolyChaos.nw(mopq)[2] == [NW[:,2] for i in 1:nunc]
end

c = [ op.α op.β ]
@testset "coefficients" begin
    @test PolyChaos.coeffs(op) == c
    @test PolyChaos.coeffs(opq) == c
    @test PolyChaos.coeffs([op for i in 1:nunc])[1] == [ c[:,1] for i in 1:nunc ]
    @test PolyChaos.coeffs([op for i in 1:nunc])[2] == [ c[:,2] for i in 1:nunc ]
    @test PolyChaos.coeffs([opq for i in 1:nunc])[1] == [ c[:,1] for i in 1:nunc ]
    @test PolyChaos.coeffs([opq for i in 1:nunc])[2] == [ c[:,2] for i in 1:nunc ]
    @test PolyChaos.coeffs(mop)[1] == [ c[:,1] for i in 1:nunc ]
    @test PolyChaos.coeffs(mop)[2] == [ c[:,2] for i in 1:nunc ]
end

@test issymmetric(op) == issymmetric(opq)


text:
__precompile__()
module CurryPipe
export @curry, @c

const _ = :_

startswith(iter, el) = (length(iter) > 0) && (first(iter) == el)

is_pipe_call(expr::Expr) = (expr.head == :call) && startswith(expr.args, :|> )
is_pipe_call(x) = false

is_fun_call(x) = false
is_fun_call(expr::Expr) = (expr.head == :call) && !(is_pipe_call(expr))

function replace!(arr, old_el, new_el)
    for i in eachindex(arr)
        if arr[i] == old_el
            arr[i] = new_el
        end
    end
    arr
end

function curry_pipe_call(expr)
    @assert length(expr.args) == 3
    waste, h, t = expr.args
    ch = curry(h)
    if is_fun_call(t)
        if (_ in t.args)
            replace!(t.args, _, ch)
        else
            push!(t.args, ch)
        end
        return t
    else
        return Expr(:call, t, ch)
    end
end

curry(x) = x
function curry(expr::Expr)
    expr = copy(expr)
    if is_pipe_call(expr)
        return curry_pipe_call(expr)
    end
    expr
end

macro curry(expr)
    expr |> curry |> esc
end

macro c(expr)
    expr |> curry |> esc
end

end # module


text:
# methods/simple.jl
#
# Method-based access to tensor operations using simple definitions.

# type unstable; better use tuples instead
tensorcopy(A, IA, IC=IA) = tensorcopy(A, tuple(IA...), tuple(IC...))
tensoradd(A, IA, B, IB, IC=IA) = tensoradd(A, tuple(IA...), B, tuple(IB...), tuple(IC...))
tensortrace(A, IA, IC = unique2(IA)) = tensortrace(A, tuple(IA...), tuple(IC...))
tensorcontract(A, IA, B, IB, IC = symdiff(IA, IB)) =
    tensorcontract(A, tuple(IA...), B, tuple(IB...), tuple(IC...))
tensorproduct(A, IA, B, IB, IC = vcat(IA, IB)) = tensorproduct(A, tuple(IA...), B, tuple(IB...), tuple(IC...))

"""
    tensorcopy(A, IA, IC = IA)

Creates a copy of `A`, where the dimensions of `A` are assigned indices from the
iterable `IA` and the indices of the copy are contained in `IC`. Both iterables
should contain the same elements in a different order.

The result of this method is equivalent to `permutedims(A, p)` where p is the permutation
such that `IC = IA[p]`. The implementation of `tensorcopy` is however more efficient on average,
especially if `Threads.nthreads() > 1`.
"""
function tensorcopy(A, IA::Tuple, IC::Tuple=IA)
    indCinA = add_indices(IA, IC)
    C = similar_from_indices(eltype(A), indCinA, (), A, :N)
    add!(1, A, :N, 0, C, indCinA)
    return C
end

"""
    tensoradd(A, IA, B, IB, IC = IA)

Returns the result of adding arrays `A` and `B` where the iterabels `IA` and `IB`
denote how the array data should be permuted in order to be added. More specifically,
the result of this method is equivalent to

```julia
tensorcopy(A, IA, IC) + tensorcopy(B, IB, IC)
```

but without creating the temporary permuted arrays.
"""
function tensoradd(A, IA::Tuple, B, IB::Tuple, IC::Tuple=IA)
    T = promote_type(eltype(A), eltype(B))
    indCinA = add_indices(IA, IC)
    C = similar_from_indices(T, indCinA, (), A, :N)
    add!(1, A, :N, 0, C, indCinA)
    indCinB = add_indices(IB, IC)
    add!(1, B, :N, 1, C, indCinB)
    return C
end

"""
    tensortrace(A, IA [, IC])

Trace or contract pairs of indices of array `A`, by assigning them an identical
indices in the iterable `IA`. The untraced indices, which are assigned a unique index,
can be reordered according to the optional argument `IC`. The default value corresponds
to the order in which they appear. Note that only pairs of indices can be contracted,
so that every index in `IA` can appear only once (for an untraced index) or twice
(for an index in a contracted pair).
"""
function tensortrace(A, IA::Tuple, IC::Tuple)
    indCinA, cindA1, cindA2 = trace_indices(IA, IC)
    C = similar_from_indices(eltype(A), indCinA, (), A, :N)
    trace!(1, A, :N, 0, C, indCinA, cindA1, cindA2)
    return C
end

"""
    tensorcontract(A, IA, B, IB[, IC])

Contract indices of array `A` with corresponding indices in array `B` by assigning
them identical labels in the iterables `IA` and `IB`. The indices of the resulting
array correspond to the indices that only appear in either `IA` or `IB` and can be
ordered by specifying the optional argument `IC`. The default is to have all open
indices of array `A` followed by all open indices of array `B`. Note that inner
contractions of an array should be handled first with `tensortrace`, so that every
label can appear only once in `IA` or `IB` seperately, and once (for open
index) or twice (for contracted index) in the union of `IA` and `IB`.

The contraction can be performed by a native Julia algorithm without creating any
temporaries, or by first permuting the arrays such that the contraction becomes equivalent
to a matrix product, which is then performed by BLAS. The latter is typically faster for
large arrays. The choice of method is globally controlled by the methods `enable_blas()`](@ref)
and `disable_blas()`](@ref).
"""
function tensorcontract(A, IA::Tuple, B, IB::Tuple, IC::Tuple)
    oindA, cindA, oindB, cindB, indCinoAB = contract_indices(IA, IB, IC)

    T = promote_type(eltype(A), eltype(B))
    C = similar_from_indices(T, oindA, oindB, indCinoAB, (), A, B, :N, :N)

    contract!(1, A, :N, B, :N, 0, C, oindA, cindA, oindB, cindB, indCinoAB)
    return C
end

"""
    tensorproduct(A, IA, B, IB, IC = (IA..., IB...))

Computes the tensor product of two arrays `A` and `B`, i.e. returns a new array `C`
with `ndims(C) = ndims(A)+ndims(B)`. The indices of the output tensor are related to
those of the input tensors by the pattern specified by the indices. Essentially,
this is a special case of `tensorcontract` with no indices being contracted over.
This method checks whether the indices indeed specify a tensor product instead of
a genuine contraction.
"""
function tensorproduct(A, IA::Tuple, B, IB::Tuple, IC::Tuple = (IA..., IB...))
    isempty(intersect(IA, IB)) || throw(IndexError("not a valid tensor product"))
    tensorcontract(A, IA, B, IB, IC)
end


text:
import Base: length

mutable struct Group
    points
    weight
    Group(x) = new(x, 1.)
end

length(gr::Group) = size(gr.points)[1]

function get_data(gr::Group)
    mean(gr.points, dims=1)
end

function init_groups(pc::Pointcloud)
    Groups([Group(pc[i:i,:]) for i = 1:size(pc)[1]])
end

Groups = Vector{Group}

function get_matrix(groups::Groups)
    points = map(get_data, groups)
    [points[i][j] for i=1:length(groups), j=1:size(points[1])[2]]
end

function add_new_group!(groups::Groups, point)
    new_group = Group(point)
    push!(groups, new_group)
end

function add_points!(group::Group, point)
    group.points = [group.points ; point]
end

function mean_initialization(measurements::Vector{Pointcloud}, c::S) where {S <: Real}
    perminds = randperm(length(measurements))
    groups = init_groups(measurements[perminds[1]])
    addpointcounter = 0
    addgroupcounter = length(groups)
    for i in perminds[2:end]
        current_points = get_matrix(groups)
        assignments, _ = hungarian(p2dist(current_points, measurements[i]))
        second_greater_first = size(current_points) < size(measurements[i])
        for j=1:length(assignments)
            ind = assignments[j]
            if assignments[j] != 0
                if p2dist(current_points[j:j,:], measurements[i][ind:ind, :])[1] < c
                    add_points!(groups[j], measurements[i][ind:ind, :])
                    addpointcounter += 1
                else
                    add_new_group!(groups, measurements[i][ind:ind, :])
                end
            end
        end
        if second_greater_first
            remaining_inds = filter(x -> !(x in assignments), 1:size(measurements[i])[1])
            for rind in remaining_inds
                add_new_group!(groups, measurements[i][rind:rind,:])
                addgroupcounter += 1
            end
        end
    end
    #return groups
    npoints = map(length, groups) |> sum
    [(get_data(g),length(g)/npoints) for g in groups]
end

function sample_from_mean_initialization(sample_space::Array{Tuple{Array{Float64,2},Float64},1},
                                         npoints::Integer,
                                         nsamples::Integer
                                        )
    sample_space_mat, weights = zip(sample_space...) |> collect
    sample_space_mat = vcat(sample_space_mat...)
    weights = [x for x = weights]
    samples = []
    for i=1:nsamples
        idxs = draw_from_weighted_distribution_without_replacement(weights, npoints)
        push!(samples, sample_space_mat[idxs,:])
    end
    samples
end

function binned_barycenter_initialization(measurements::Vector{Pointcloud}
                                        ,weights::Vector{R}
                                        ,c::S
                                        ,min_card::Integer
                                        ,max_card::Integer
                                        ,eps::T
                                        ,nsamples::Integer
                                        ) where {R <: Real, S <: Real, T <: Real}
    sample_space = mean_initialization(measurements, c)
    if max_card == 0
        max_card = length(sample_space)
    end
    max_card = min(max_card, length(sample_space))
    Barycenters = Vector{Pointcloud}()
    costs = Vector{Float64}()
    for i = min_card:max_card
        best_cost = Inf
        best_bc = missing
        bc_inits = sample_from_mean_initialization(sample_space, i, nsamples)
        for _=1:nsamples
            bc_init = sample_from_mean_initialization(sample_space, i, 1)
            cur_bc, cur_cost = ospa_barycenter(measurements, weights, c, bc_init[1], eps)
            if cur_cost < best_cost
                best_cost = cur_cost
                best_bc = cur_bc
            end
        end
#        for sample in bc_inits
#            cur_bc, cur_cost = ospa_barycenter(measurements, weights, c, sample, eps)
#            if cur_cost < best_cost
#                best_cost = cur_cost
#                best_bc = cur_bc
#            end
#        end
        push!(Barycenters, best_bc)
        push!(costs, best_cost)
    end
    (_, ind) = findmin(costs)
    Barycenters[ind], costs[ind]
end


text:
module Vtkio

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery

using WriteVTK
using WriteVTK.VTKCellTypes: VTK_VERTEX
using WriteVTK.VTKCellTypes: VTK_LINE
using WriteVTK.VTKCellTypes: VTK_TRIANGLE
using WriteVTK.VTKCellTypes: VTK_QUAD
using WriteVTK.VTKCellTypes: VTK_TETRA
using WriteVTK.VTKCellTypes: VTK_HEXAHEDRON

# Functionality given by this module

export writevtk
import Gridap: cells, points, celltypes
import UnstructuredGrids: writevtk

"""
Write a Grid object into a vtk file
"""
function writevtk(grid::Grid,filebase;celldata=Dict(),pointdata=Dict())
  _writevtk(grid,filebase,celldata,pointdata)
end

"""
Write a CellPoints object into vtk file
"""
function writevtk(points::CellPoints,filebase;celldata=Dict(),pointdata=Dict())
  _writevtk(points,filebase,celldata,pointdata)
end

"""
Write a CellValue{Point{D}} object into vtk file
"""
function writevtk(points::CellValue{Point{D,Float64}} where D,filebase;celldata=Dict(),pointdata=Dict())
  _writevtk(points,filebase,celldata,pointdata)
end

"""
Write an Triangulation object into vtk
"""
function writevtk(trian::Triangulation,filebase;nref=0,celldata=Dict(),cellfields=Dict())
  _writevtk(trian,filebase,nref,celldata,cellfields)
end

"""
Write a Polytope object into vtk
"""
function writevtk(polytope::Polytope,filebase)
  refcell = RefCell(polytope)
  writevtk(refcell,filebase)
end

"""
Write a DiscreteModel into vtk
"""
function writevtk(model::DiscreteModel{D},filebase) where D
  labels = FaceLabels(model)
  for d in 0:D
    f = "$(filebase)_$(d)"
    grid = Grid(model,d)
    cdat = _prepare_cdata_model(labels,d)
    writevtk(grid,f,celldata=cdat)
  end
end

# Helpers

function _prepare_cdata_model(labels,d)
  dface_to_label = labels_on_dim(labels,d)
  ndfaces = length(dface_to_label)
  cdat = []
  for tag in 1:ntags(labels)
    dface_to_isontag = zeros(Int,ndfaces)
    for label in labels_on_tag(labels,tag)
      _set_label!(dface_to_isontag,dface_to_label,label,tag)
    end
    name = name_from_tag(labels,tag)
    push!(cdat, name => dface_to_isontag )
  end
  push!(cdat,"geolabel" => dface_to_label)
  cdat
end

function _set_label!(dface_to_isontag,dface_to_label,label,tag)
  for i in 1:length(dface_to_label)
    if dface_to_label[i] == label
      dface_to_isontag[i] = tag
    end
  end
end

function _writevtk(grid::Grid,filebase,celldata,pointdata)
  points = _vtkpoints(grid)
  cells = _vtkcells(grid)
  vtkfile = vtk_grid(filebase, points, cells, compress=false)
  for (k,v) in celldata
    vtk_cell_data(vtkfile, _prepare_data(v), k)
  end
  for (k,v) in pointdata
    vtk_point_data(vtkfile, _prepare_data(v), k)
  end
  outfiles = vtk_save(vtkfile)
end

function _vtkpoints(grid::Grid{D}) where D
  x = points(grid)
  xflat = collect(x)
  reshape(reinterpret(Float64,xflat),(D,length(x)))
end

# @fverdugo this allocates a lot of small objects
# Not very crucial since it is for visualization
# but it would be nice to have a better way
function _vtkcells(grid::Grid)
  types = _vtkcelltypedict()
  nodes = _vtkcellnodesdict()
  c = celltypes(grid)
  n = cells(grid)
  _check_order(cellorders(grid))
  [ MeshCell(types[encode_extrusion(ci)], ni[nodes[encode_extrusion(ci)]])
     for (ci,ni) in zip(c,n) ]
end

_check_order(co) = @notimplemented

function _check_order(co::ConstantCellValue{Int})
  o = co.value
  @notimplementedif o != 1
end

"""
Encodes the tuple defining a Polytope into an integer
"""
function encode_extrusion(extrusion::NTuple{Z,Int}) where Z
  k = 0
  for (i,v) in enumerate(extrusion)
    k += v*3^i
  end
  k
end

"""
Decodes an integer into a tuple defining a Polytope
"""
function decode_extrusion(i::Int,::Val{Z}) where Z
  @notimplemented
end

# @fverdugo This can also be done by dispatching on value
"""
Generates the lookup table (as a Dict) in order to convert between
Gridap Polytope identifiers into VTK cell type identifiers
"""
function _vtkcelltypedict()
  d = Dict{Int,WriteVTK.VTKCellTypes.VTKCellType}()
  h = HEX_AXIS
  t = TET_AXIS
  d[encode_extrusion(())] = VTK_VERTEX
  d[encode_extrusion((t,))] = VTK_LINE
  d[encode_extrusion((h,))] = VTK_LINE
  d[encode_extrusion((t,t))] = VTK_TRIANGLE
  d[encode_extrusion((h,h))] = VTK_QUAD
  d[encode_extrusion((t,t,t))] = VTK_TETRA
  d[encode_extrusion((h,h,h))] = VTK_HEXAHEDRON
  d
end

# @fverdugo This can also be done by dispatching on value
"""
Generates the lookup table (as a Dict) in order to convert between
Gridap Polytope corner numbering into VTK corner numbering
"""
function _vtkcellnodesdict()
  d = Dict{Int,Vector{Int}}()
  h = HEX_AXIS
  t = TET_AXIS
  d[encode_extrusion(())] = [1,]
  d[encode_extrusion((t,))] = [1,2]
  d[encode_extrusion((h,))] = [1,2]
  d[encode_extrusion((t,t))] = [1,2,3]
  d[encode_extrusion((h,h))] = [1,2,4,3]
  d[encode_extrusion((t,t,t))] = [1,2,3,4]
  d[encode_extrusion((h,h,h))] = [1,2,4,3,5,6,8,7]
  d
end

function _writevtk(points::CellPoints,filebase,celldata,pointdata)
  grid, p_to_cell = _cellpoints_to_grid(points)
  pdat = _prepare_pointdata(pointdata)
  k = "cellid"
  @assert ! haskey(pdat,k)
  pdat[k] = p_to_cell
  writevtk(grid,filebase,pointdata=pdat)
end

function _writevtk(points::CellValue{Point{D,Float64}} where D,filebase,celldata,pointdata)
  grid = _cellpoint_to_grid(points)
  pdat = _prepare_pointdata(pointdata)
  writevtk(grid,filebase,pointdata=pdat)
end

function _cellpoints_to_grid(points::CellPoints{D}) where D
  ps = Array{Point{D,Float64},1}(undef,(0,))
  p_to_cell = Array{Int,1}(undef,(0,))
  for (cell,p) in enumerate(points)
    for pj in p
      push!(ps,pj)
      push!(p_to_cell,cell)
    end
  end
  data, ptrs, ts, os = _prepare_cells(ps)
  grid = UnstructuredGrid(ps,data,ptrs,ts,os)
  (grid, p_to_cell)
end

function _cellpoint_to_grid(points::CellValue{Point{D,Float64}}) where D
  ps = collect(points)
  data, ptrs, ts, os = _prepare_cells(ps)
  UnstructuredGrid(ps,data,ptrs,ts,os)
end

function _prepare_cells(ps)
  data = [ i for i in 1:length(ps) ]
  ptrs = [ i for i in 1:(length(ps)+1) ]
  ts = ConstantCellValue( (), length(ps) )
  os = ConstantCellValue( 1, length(ps) )
  (data,ptrs,ts,os)
end

function _prepare_pointdata(pointdata)
  pdat = Dict()
  for (k,v) in pointdata
    pdat[k] = _prepare_data(v)
  end
  pdat
end

_prepare_data(v) = v

const IterData{T} = Union{CellValue{T},AbstractArray{T}}

function _prepare_data(v::IterData{<:VectorValue{D}}) where D
  a = collect(v)
  reshape(reinterpret(Float64,a),(D,length(a)))
end

function _prepare_data(v::IterData{<:VectorValue{2}})
  a = collect(v)
  b = reshape(reinterpret(Float64,a),(2,length(a)))
  z = zeros((1,size(b,2)))
  vcat(b,z)
end

function _prepare_data(v::IterData{<:TensorValue{D}}) where D
  a = collect(v)
  reshape(reinterpret(Float64,a),(D*D,length(a)))
end

_prepare_data(v::CellArray{<:Number}) = _flatten(v)

function _prepare_data(v::CellArray{<:VectorValue{D}}) where D
  a = _flatten(v)
  reshape(reinterpret(Float64,a),(D,length(a)))
end

function _prepare_data(v::CellArray{<:VectorValue{2}})
  a = _flatten(v)
  b = reshape(reinterpret(Float64,a),(2,length(a)))
  z = zeros((1,size(b,2)))
  vcat(b,z)
end

function _prepare_data(v::CellArray{<:TensorValue{D}}) where D
  a = _flatten(v)
  reshape(reinterpret(Float64,a),(D*D,length(a)))
end

struct VisualizationGrid{D,Z,G<:Grid{D,Z},C<:IndexCellValue{Int},P<:CellPoints{Z}} <: Grid{D,Z}
  grid::G
  coarsecells::C
  samplingpoints::P
end

points(vg::VisualizationGrid) = points(vg.grid)

cells(vg::VisualizationGrid) = cells(vg.grid)

celltypes(vg::VisualizationGrid) = celltypes(vg.grid)

function _writevtk(vg::VisualizationGrid,filebase,celldata,cellfields)
  cdata = _prepare_cdata(celldata,vg.coarsecells)
  pdata = _prepare_pdata(cellfields,vg.samplingpoints)
  _writevtk(vg.grid,filebase,cdata,pdata)
end

function _prepare_cdata(celldata,fine_to_coarse)
  cdata = Dict()
  for (k,v) in celldata
    acoarse = collect(v)
    afine = _allocate_afine(acoarse,length(fine_to_coarse))
    _fill_afine!(afine,acoarse,fine_to_coarse)
    cdata[k] = afine
  end
  k2 = "cellid"
  @assert ! haskey(cdata,k2)
  cdata[k2] = fine_to_coarse
  cdata
end

_allocate_afine(acoarse::Array{T},l) where T = Array{T,1}(undef,(l,))

function _fill_afine!(afine,acoarse,fine_to_coarse)
  for (i,coarse) in enumerate(fine_to_coarse)
    afine[i] = acoarse[coarse]
  end
end

function _prepare_pdata(cellfields,samplingpoints)
  pdata = Dict()
  for (k,v) in cellfields
    pdata[k] = _flatten(evaluate(v,samplingpoints))
  end
  pdata
end

function _writevtk(trian::Triangulation,filebase,nref,celldata,cellfields)
  vg = _visgrid(trian,nref)
  _writevtk(vg,filebase,celldata,cellfields)
end

function _visgrid(self::Triangulation,nref)
  grid, coarsecells, samplingpoints = _prepare_grid(CellRefFEs(self),CellGeomap(self),nref)
  VisualizationGrid(grid,coarsecells,samplingpoints)
end

function _refgrid end

function _refgrid(::Val{V},nref::Int) where V
  @notimplementedif nref != 0
  polytope = Polytope(V)
  Grid(polytope,length(V))
end

function _refgrid(::Val{(HEX_AXIS,)},nref::Int)
  n = 2^nref
  CartesianGrid(domain=(0.0,1.0),partition=(n,))
end

function _refgrid(::Val{(HEX_AXIS,HEX_AXIS)},nref::Int)
  n = 2^nref
  CartesianGrid(domain=(0.0,1.0,0.0,1.0),partition=(n,n))
end

function _refgrid(::Val{(HEX_AXIS,HEX_AXIS,HEX_AXIS)},nref::Int)
  n = 2^nref
  CartesianGrid(domain=(0.0,1.0,0.0,1.0,0.0,1.0),partition=(n,n,n))
end

function _prepare_grid(reffes::CellValue{<:RefFE},phi::CellGeomap{Z,D},nref) where {Z,D}
  @notimplemented
end

function _prepare_grid(reffes::ConstantCellValue{<:RefFE},phi::CellGeomap{Z,D},nref) where {Z,D}
  refgrid = _prepare_refgrid(reffes,nref)
  samplingpoints = _prepare_samplingpoints(refgrid,reffes)
  ps, offsets = _prepare_points(samplingpoints,points(refgrid),phi)
  data, ptrs, coarsecells = _prepare_cells(refgrid,offsets)
  ts = _prepare_celltypes(length(reffes),celltypes(refgrid))
  os = _prepare_cellorders(length(reffes),celltypes(refgrid))
  grid = UnstructuredGrid(ps,data,ptrs,ts,os)
  (grid, coarsecells, samplingpoints)
end

function _prepare_refgrid(reffes,nref)
  reffe = reffes.value
  poly = polytope(reffe)
  extrusion = poly.extrusion.array
  _refgrid(Val(Tuple(extrusion)),nref)
end

function _prepare_samplingpoints(refgrid,reffes)
  refpoints = collect(points(refgrid))
  ConstantCellValue(refpoints[:],length(reffes))
end

function _prepare_points(samplingpoints,refpoints,phi::CellGeomap{Z,D}) where {Z,D}
  xe = evaluate(phi,samplingpoints)
  offsets = Array{Int,1}(undef,(length(xe),))
  ps = Array{Point{D,Float64},1}(undef,(length(xe)*length(refpoints)))
  _fill_ps_and_offsets!(ps,offsets,xe)
  (ps, offsets)
end

function _fill_ps_and_offsets!(ps,offsets,xe)
  k = 1
  for (n,x) in enumerate(xe)
    offsets[n] = k - 1
    for xi in x
      @inbounds ps[k] = xi
      k += 1
    end
  end
end

function _prepare_cells(refgrid,offsets)
  refcells = cells(refgrid)
  ptrs = Array{Int,1}(undef,(1+length(offsets)*length(refcells),))
  coarsecells = Array{Int,1}(undef,(length(offsets)*length(refcells),))
  _fill_ptrs_and_coarsecells!(ptrs,coarsecells,refcells,length(offsets))
  length_to_ptrs!(ptrs)
  data = Array{Int,1}(undef,(ptrs[end]-1,))
  _fill_data!(data,offsets,refcells)
  (data, ptrs, CellValueFromArray(coarsecells))
end

function _fill_ptrs_and_coarsecells!(ptrs,coarsecells,refcells,ncells)
  k = 1
  for cell in 1:ncells
    for refnodes in refcells
      @inbounds coarsecells[k] = cell
      k += 1
      @inbounds ptrs[k] = length(refnodes)
    end
  end
end

function _fill_data!(data,offsets,refcells)
  k = 1
  for offset in offsets
    for refnodes in refcells
      for node in refnodes
        @inbounds data[k] = node + offset
        k += 1
      end
    end
  end
end

function _prepare_celltypes(ncells,refcelltypes::ConstantCellValue)
  refextrusion = refcelltypes.value
  ConstantCellValue(refextrusion,ncells*length(refcelltypes) )
end

function _prepare_cellorders(ncells,refcelltypes::ConstantCellValue)
  ConstantCellValue(1,ncells*length(refcelltypes) )
end

function _flatten(ca::CellArray{T}) where T
  f = Vector{T}(undef,0)
  _flatten_kernel!(f,ca)
  f
end

function _flatten_kernel!(f,ca)
  for a in ca
    for ai in a
      push!(f,ai)
    end
  end
end

end # module Vtkio


text:
"""
Union of several basis functions in the imaginary-time/real-frequency domain
domains
"""
struct CompositeBasisFunction
    polys::Vector{PiecewiseLegendrePoly}
end

"""
Evaluate basis function at position x
"""
function (obj::CompositeBasisFunction)(x::Real)
    return vcat((p(x) for p in obj.polys))
end

function (obj::CompositeBasisFunction)(x::Vector{T}) where {T<:Real}
    return vcat((p(x) for p in obj.polys)...)
end

"""
Union of several basis functions in the imaginary-frequency domain
domains
"""
struct CompositeBasisFunctionFT
    polys::Vector{PiecewiseLegendreFT}
end

"""
Evaluate basis function at frequency n
"""
function (obj::CompositeBasisFunctionFT)(n::Union{Int64,Vector{Int64}})
    return hcat((p(n) for p in obj.polys))
end

struct CompositeBasis <: Basis
    beta::Float64
    size::Int64
    bases::Vector{Basis}
    u::Union{CompositeBasisFunction,Nothing}
    v::Union{CompositeBasisFunction,Nothing}
    uhat::Union{CompositeBasisFunctionFT,Nothing}
end

function _collect_polys(::Type{T}, polys) where {T}
    if any((p === nothing for p in polys))
        return nothing
    else
        return T([p for p in polys])
    end
end

function CompositeBasis(bases::Vector{Basis})
    size = sum((b.size for b in bases))
    u = CompositeBasisFunction([b.u for b in bases])
    v = _collect_polys(CompositeBasisFunction, [b.v for b in bases])
    uhat = _collect_polys(CompositeBasisFunctionFT, [b.uhat for b in bases])
    return CompositeBasis(bases[1].beta, size, bases, u, v, uhat)
end

function default_tau_sampling_points(basis::CompositeBasis)
    return sort(unique(vcat((default_tau_sampling_points(b) for b in basis.bases)...)))
end

function default_matsubara_sampling_points(basis::CompositeBasis; mitigate=true)
    return sort(unique(vcat((default_matsubara_sampling_points(b; mitigate=mitigate) for b in
                                                                                basis.bases)...)))
end


text:
# Struct representing the final result
mutable struct Result
    # Result header - utility for printResult
    header::String
    # Array of initial dataset points
    dataset::Vector
    # Size of dataset
    dataset_size::Int64
    # Array of solvingSet points
    solvingSet::Vector
    # Size of solvingSet
    solvingSet_size::Int64
    # Array of solution points
    solution::Vector
    # Array of candidate points for each iteration
    itC::Vector
    # Array of solving set for each iteration
    itS::Vector
    # Number of distances calculated by the algorithm
    distCount::Int64

    time::Int64

    Result() = new(Vector(),Vector(),Vector(),Vector(),Vector(), 0, 0)
    Result(header::String, dataset::Vector, solvingSet::Vector, solution::Vector, itC::Vector,itS::Vector,distCount::Int64, time::Int64) = new(header, dataset, length(dataset),solvingSet,length(solvingSet),solution,itC,itS,distCount,time)
end

function getReport(r::Result)
    tmp = "\n\n"*r.header*'\n'
    tmp *= "_____________________________________\n\n"
    tmp *= "|D|\t|S|\t#d (bilions)\n"
    tmp *= "_____________________________________\n\n"
    tmp*= string(r.dataset_size)*'\t'string(r.solvingSet_size)*'\t'*string(r.distCount)*'\n'
    tmp *= "_____________________________________\n\n"
    return tmp
end

function getReport(v::Vector{Result})
    tmp = "\n\n"*v[1].header*'\n'
    tmp *= "_____________________________________\n\n"
    tmp *= "|D|\t|S|\t#d (bilions)\n"
    tmp *= "_____________________________________\n\n"
    for r in v
        tmp*= string(r.dataset_size)*'\t'string(r.solvingSet_size)*'\t'*string(r.distCount)*'\n'
    end
    tmp *= "_____________________________________\n\n"
    return tmp
end


text:
using Statistics
using Serialization
using Printf
using ArgParse
using PrettyTables

const JULIAVER = Base.julia_cmd()[1]

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--runs", "-n"
            help = "number of iterations"
            arg_type = Int
            default = 10
        "--threads", "-t"
            help = "number of threads"
            arg_type = Int
            default = 1
        "--bench", "-b"
            help = "if specified, runs a single benchmark"
            arg_type = String
            default = "all"
    end

    return parse_args(s)
end

# times in ns
# TODO: get better stats
function get_stats(times::Vector)
    return [minimum(times), median(times), maximum(times)]
end

"""
    Highlights cells in a column based on value
        green if less than lo
        yellow if between lo and hi
        red if above hi
"""
function highlight_col(col, lo, hi)
    [Highlighter((data,i,j) -> (j == col) && data[i, j] <= lo; foreground=:green),
     Highlighter((data,i,j) -> (j == col) && lo < data[i, j] < hi; foreground=:yellow),
     Highlighter((data,i,j) -> (j == col) && hi <= data[i, j]; foreground=:red),]
end

function run_one_bench(runs, threads, file)
    value = []
    times = []
    gc_diff = []
    gc_end = []
    for _ in 1:runs
        # uglyness to communicate over non stdout (specifically file descriptor 3)
        p = Base.PipeEndpoint()
        cmd = `$JULIAVER --project=. --threads=$threads $file SERIALIZE`
        cmd = run(Base.CmdRedirect(cmd, p, 3), stdin, stdout, stderr, wait=false)
        r = deserialize(p)
        @assert success(cmd)
        # end uglyness
        push!(value, r.value)
        push!(times, r.times)
        push!(gc_diff, r.gc_diff)
        push!(gc_end, r.gc_end)
    end
    total_stats = get_stats(times) ./ 1_000_000
    gc_time = get_stats(map(stat->stat.total_time, gc_end)) ./ 1_000_000
    max_pause = get_stats(map(stat->stat.max_pause, gc_end)) ./ 1_000_000
    time_to_safepoint = get_stats(map(stat->stat.time_to_safepoint, gc_end)) ./ 1_000_000
    max_mem = get_stats(map(stat->stat.max_memory, gc_end)) ./ 1024^2
    pct_gc = get_stats(map((t,stat)->(stat.total_time/t), times, gc_diff)) .* 100

    header = (["", "total time", "gc time", "max GC pause", "time to safepoint", "max heap", "percent gc"],
              ["", "ms",         "ms",       "ms",          "ms",                "MB",       "%"        ])
    labels = ["minimum", "median", "maximum"]
    highlighters = highlight_col(4, 10, 100) # max pause
    append!(highlighters, highlight_col(5, 1, 10)) # time to safepoint
    append!(highlighters, highlight_col(7, 10, 50)) # pct gc
    highlighters = Tuple(highlighters)
    data = hcat(labels, total_stats, gc_time, max_pause, time_to_safepoint, max_mem, pct_gc)
    pretty_table(data; header, formatters=ft_printf("%0.0f"), highlighters)
end

function run_all_benches(runs, threads)
    for category in readdir()
        @show category
        cd(category)
        for file in readdir()
            endswith(file, ".jl") || continue
            @show file
            run_one_bench(runs, threads, file)
        end
        cd("..")
    end
end

function main()
    args = parse_commandline()
    JULIAVER = Base.julia_cmd()[1]
    bench, runs, threads = args["bench"], args["runs"], args["threads"]

    dir = joinpath(@__DIR__, "benches")
    cd(dir)
    if bench == "all"
        run_all_benches(runs, threads)
    else
        path, file = splitpath(bench)
        cd(path)
        run_one_bench(runs, threads, file)
    end
end

main()

max_stars_repo_name: JuliaCI/GCBenchmarks

max_issues_repo_name: JuliaCI/GCBenchmarks

max_forks_repo_name: JuliaCI/GCBenchmarks


text:
#**************************************************************************************
# build.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file excuted at instlation time.
#
#**************************************************************************************

# creating `GeoEfficiency` folder at the home directory.

@info("Creating 'GeoEfficiency' folder at '$(homedir())'.....")
try
	cp(joinpath(dirname(@__FILE__),".batch"), joinpath(homedir(),"GeoEfficiency"))
catch err
	if isdir(joinpath(homedir(),"GeoEfficiency"))
		@warn("`GeoEfficiency` folder already exist.")
	else
		@error("""'GeoEfficiency' folder: could not be created,
		\n  the package may be unable to work in patch mode. Try 'import Pkg; Pkg.build("GeoEfficiency")' to rebuild.""")
	end #if
end #try


text:
using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--epochs"
        help = "Epochs to train"
        arg_type = Int
        default = 100
    "--save-snapshot"
        help = "Whether to Save Snapshots"
        arg_type = Bool
        default = true
    "--snapshot-dir"
        help = "Directory to Store Snapshots"
        arg_type = ASCIIString
        default = "snapshot-rnn"
    "--snapshot_n_iter"
        help = "How many iterations between snapshots"
        arg_type = Int
        default = 10
    "--stat_n_iter"
        help = "How many iterations between saving statistics"
        arg_type = Int
        default = 10
end
const args = parse_args(ARGS, s)

ENV["SEEP_NO_GPU"]=1
using Seep, JLD
T = Float32
relu(x::ANode) = Seep.max!(0, x)
make_node(a::Array) = ANode(a)
make_node(x::Int...) = ANode(x...)
get_data(x) = x
const pool = BuddyPool(Array{T}(1<<24))

include("data.jl")
data = IRData()

# Constants
const sequence_length = 50
const input_size = 2
const hidden_size = [10, 10]
const feature_size = 5
const output_size = 2
const batch_size = 80

# Network
const X = ANode[make_node(input_size, batch_size) for i=1:sequence_length]
const Y = make_node(output_size, batch_size)
y = make_node(zeros(T, output_size, batch_size))

params = Flow[]
let
    # Initialize nodes
    fe = Linear(T, "e", input_size, hidden_size[1])
    fh = Flow[]
    for j = 1:length(hidden_size)
        k = j==1 ? 1 : j-1

        f = LSTM(T, "h$(j)", hidden_size[k], hidden_size[j])
        push!(fh, f)
    end
    ff = Linear(T, "features", hidden_size[end], feature_size)
    fy = Linear(T, "output", feature_size, output_size)
    push!(params, fe, fh..., ff, fy)

    # Initialize LSTM variables
    h = ANode[make_node(zeros(T, hidden_size[i], batch_size)) for i = 1:length(hidden_size)]
    c = ANode[make_node(zeros(T, hidden_size[i], batch_size)) for i = 1:length(hidden_size)]

    # Build graph
    x = make_node(zeros(T, input_size, batch_size))
    for i = 1:length(X)
        x = relu(fe(X[i]))
        for j = 1:length(hidden_size)
            c[j], h[j] = fh[j](x, c[j], h[j])
            x = h[j]
        end
    end
    x = tanh(ff(x))
    y = softmax(fy(x), 1)
end

# Loss
const loss = -log(y) .* Y
const obj = sum(loss)
const grad = gradients(obj)

# Gradient Calculation for SGD
updates = map(x->Seep.adam(x, grad[x]), get_params(params))
const RNN = instance(pool, y, loss, updates...)
Seep.print_stats(STDOUT, pool)

# Load Latest Snapshot
n_epoch_start = 0

if args["save-snapshot"]
    if isdir(args["snapshot-dir"])
        f = filter(x->contains(x, "snapshot-"), readdir(args["snapshot-dir"]))
        if length(f) > 0
            epochs = map(x->parse(Int, split(split(x, ".")[1], "-")[end]), f)
            ind = indmax(epochs)
            n_epoch_start = epochs[ind]
            @show n_epoch_start
            load_snapshot(joinpath(args["snapshot-dir"], f[ind]), loss)
        end
    else
        mkdir(args["snapshot-dir"])
    end
end

# Train
function epoch()
    n = data.n_epochs
    while data.n_epochs == n
        xx, yy = next_batch(data, batch_size)

        for i = 1:length(X)
            xi = RNN[X[i]]
            xi[:, :] = xx[:, i, :]
        end

        y = zeros(T, output_size, batch_size)
        for i = 1:length(yy)
            y[round(Int, yy[i]), i] = 1
        end
        yi = RNN[Y]
        yi[:] = y

        RNN()
    end
end

E = zeros(args["stat_n_iter"])
for i = 1:args["epochs"]
    epoch()

    l = sum(get_data(RNN[loss]))/batch_size
    println("Loss: ", l)

    E[mod(i-1, args["stat_n_iter"])+1] = l
    if mod(i, args["stat_n_iter"])==0
        f = joinpath(args["snapshot-dir"], "statistics.jld")
        if isfile(f)
            E = vcat(load(f, "E"), E)
        end
        save(f, "E", E)
        E = zeros(args["stat_n_iter"])
    end

    if mod(i, args["snapshot_n_iter"])==0 && args["save-snapshot"]
        save_snapshot(joinpath(args["snapshot-dir"], "snapshot-$(i+n_epoch_start).jld"), loss)
    end
end


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

#using DelayEmbeddings
using DelimitedFiles
using StatsBase
using Plots
using StatsPlots
using Statistics

include("../../../src/pecuzal_method.jl")
include("../../../src/data_analysis_functions.jl")

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we look at real data of a chaotic electrochemical oscillator (Ni) from
# the Kiss group with two different resistances R₁ = 1,5kΩ, R₂ = 2,5kΩ leading to
# "more" chaoticity.

dim_tde_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde_ref1.csv")
dim_tde_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde_ref2.csv")
dim_GA_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA_ref1.csv")
dim_GA_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA_ref2.csv")
dim_mdop_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop_ref1.csv")
dim_mdop_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop_ref2.csv")
dim_pec_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec_ref1.csv")
dim_pec_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec_ref2.csv")

Tw_tde_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde_ref1.csv")
Tw_tde_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde_ref2.csv")
Tw_GA_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA_ref1.csv")
Tw_GA_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA_ref2.csv")
Tw_mdop_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop_ref1.csv")
Tw_mdop_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop_ref2.csv")
Tw_pec_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec_ref1.csv")
Tw_pec_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec_ref2.csv")

RQA_tde_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde_ref1.csv")
RQA_tde_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde_ref2.csv")
RQA_GA_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA_ref1.csv")
RQA_GA_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA_ref2.csv")
RQA_mdop_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop_ref1.csv")
RQA_mdop_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop_ref2.csv")
RQA_pec_ref1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec_ref1.csv")
RQA_pec_ref2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec_ref2.csv")

RQA_tde1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde1.csv")
RQA_tde2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde2.csv")
RQA_GA1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA1.csv")
RQA_GA2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA2.csv")
RQA_mdop1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop1.csv")
RQA_mdop2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop2.csv")
RQA_pec1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec1.csv")
RQA_pec2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec2.csv")

dim_tde1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde1.csv")
dim_tde2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde2.csv")
dim_GA1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA1.csv")
dim_GA2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA2.csv")
dim_mdop1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop1.csv")
dim_mdop2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop2.csv")
dim_pec1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec1.csv")
dim_pec2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec2.csv")

Tw_tde1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde1.csv")
Tw_tde2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde2.csv")
Tw_GA1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA1.csv")
Tw_GA2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA2.csv")
Tw_mdop1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop1.csv")
Tw_mdop2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop2.csv")
Tw_pec1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec1.csv")
Tw_pec2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec2.csv")

# basic statistics
mean_RQA_tde1 = mean(RQA_tde1; dims=2)
mean_RQA_tde2 = mean(RQA_tde2; dims=2)
median_RQA_tde1 = median(RQA_tde1; dims=2)
median_RQA_tde2 = median(RQA_tde2; dims=2)
std_RQA_tde1 = std(RQA_tde1; dims=2)
std_RQA_tde2 = std(RQA_tde2; dims=2)
mean_Tw_tde1 = mean(Tw_tde1)
mean_Tw_tde2 = mean(Tw_tde2)
std_Tw_tde1 = std(Tw_tde1)
std_Tw_tde2 = std(Tw_tde2)
mean_dim_tde1 = mean(dim_tde1)
mean_dim_tde2 = mean(dim_tde2)
std_dim_tde1 = std(dim_tde1)
std_dim_tde2 = std(dim_tde2)

mean_RQA_GA1 = mean(RQA_GA1; dims=2)
mean_RQA_GA2 = mean(RQA_GA2; dims=2)
median_RQA_GA1 = median(RQA_GA1; dims=2)
median_RQA_GA2 = median(RQA_GA2; dims=2)
std_RQA_GA1 = std(RQA_GA1; dims=2)
std_RQA_GA2 = std(RQA_GA2; dims=2)
mean_Tw_GA1 = mean(Tw_GA1)
mean_Tw_GA2 = mean(Tw_GA2)
std_Tw_GA1 = std(Tw_GA1)
std_Tw_GA2 = std(Tw_GA2)
mean_dim_GA1 = mean(dim_GA1)
mean_dim_GA2 = mean(dim_GA2)
std_dim_GA1 = std(dim_GA1)
std_dim_GA2 = std(dim_GA2)

mean_RQA_mdop1 = mean(RQA_mdop1; dims=2)
mean_RQA_mdop2 = mean(RQA_mdop2; dims=2)
median_RQA_mdop1 = median(RQA_mdop1; dims=2)
median_RQA_mdop2 = median(RQA_mdop2; dims=2)
std_RQA_mdop1 = std(RQA_mdop1; dims=2)
std_RQA_mdop2 = std(RQA_mdop2; dims=2)
mean_Tw_mdop1 = mean(Tw_mdop1)
mean_Tw_mdop2 = mean(Tw_mdop2)
std_Tw_mdop1 = std(Tw_mdop1)
std_Tw_mdop2 = std(Tw_mdop2)
mean_dim_mdop1 = mean(dim_mdop1)
mean_dim_mdop2 = mean(dim_mdop2)
std_dim_mdop1 = std(dim_mdop1)
std_dim_mdop2 = std(dim_mdop2)

mean_RQA_pec1 = mean(RQA_pec1; dims=2)
mean_RQA_pec2 = mean(RQA_pec2; dims=2)
median_RQA_pec1 = median(RQA_pec1; dims=2)
median_RQA_pec2 = median(RQA_pec2; dims=2)
std_RQA_pec1 = std(RQA_pec1; dims=2)
std_RQA_pec2 = std(RQA_pec2; dims=2)
mean_Tw_pec1 = mean(Tw_pec1)
mean_Tw_pec2 = mean(Tw_pec2)
std_Tw_pec1 = std(Tw_pec1)
std_Tw_pec2 = std(Tw_pec2)
mean_dim_pec1 = mean(dim_pec1)
mean_dim_pec2 = mean(dim_pec2)
std_dim_pec1 = std(dim_pec1)
std_dim_pec2 = std(dim_pec2)

# compute histograms

p_tde_ENTR1 = histogram(RQA_tde1[1,:],title = "TDE 1 ENTR")
p_tde_ENTR2 = histogram(RQA_tde2[1,:],title = "TDE 2 ENTR")
p_tde_LAM1 = histogram(RQA_tde1[2,:],title = "TDE 1 LAM")
p_tde_LAM2 = histogram(RQA_tde2[2,:],title = "TDE 2 LAM")
p_tde_RTE1 = histogram(RQA_tde1[3,:],title = "TDE 1 RTE")
p_tde_RTE2 = histogram(RQA_tde2[3,:],title = "TDE 2 RTE")
p_tde_TRANS1 = histogram(RQA_tde1[4,:],title = "TDE 1 TRANS")
p_tde_TRANS2 = histogram(RQA_tde2[4,:],title = "TDE 2 TRANS")

p_GA_ENTR1 = histogram(RQA_GA1[1,:],title = "GA 1 ENTR")
p_GA_ENTR2 = histogram(RQA_GA2[1,:],title = "GA 2 ENTR")
p_GA_LAM1 = histogram(RQA_GA1[2,:],title = "GA 1 LAM")
p_GA_LAM2 = histogram(RQA_GA2[2,:],title = "GA 2 LAM")
p_GA_RTE1 = histogram(RQA_GA1[3,:],title = "GA 1 RTE")
p_GA_RTE2 = histogram(RQA_GA2[3,:],title = "GA 2 RTE")
p_GA_TRANS1 = histogram(RQA_GA1[4,:],title = "GA 1 TRANS")
p_GA_TRANS2 = histogram(RQA_GA2[4,:],title = "GA 2 TRANS")

p_mdop_ENTR1 = histogram(RQA_mdop1[1,:],title = "mdop 1 ENTR")
p_mdop_ENTR2 = histogram(RQA_mdop2[1,:],title = "mdop 2 ENTR")
p_mdop_LAM1 = histogram(RQA_mdop1[2,:],title = "mdop 1 LAM")
p_mdop_LAM2 = histogram(RQA_mdop2[2,:],title = "mdop 2 LAM")
p_mdop_RTE1 = histogram(RQA_mdop1[3,:],title = "mdop 1 RTE")
p_mdop_RTE2 = histogram(RQA_mdop2[3,:],title = "mdop 2 RTE")
p_mdop_TRANS1 = histogram(RQA_mdop1[4,:],title = "mdop 1 TRANS")
p_mdop_TRANS2 = histogram(RQA_mdop2[4,:],title = "mdop 2 TRANS")

p_pec_ENTR1 = histogram(RQA_pec1[1,:],title = "pec 1 ENTR")
p_pec_ENTR2 = histogram(RQA_pec2[1,:],title = "pec 2 ENTR")
p_pec_LAM1 = histogram(RQA_pec1[2,:],title = "pec 1 LAM")
p_pec_LAM2 = histogram(RQA_pec2[2,:],title = "pec 2 LAM")
p_pec_RTE1 = histogram(RQA_pec1[3,:],title = "pec 1 RTE")
p_pec_RTE2 = histogram(RQA_pec2[3,:],title = "pec 2 RTE")
p_pec_TRANS1 = histogram(RQA_pec1[4,:],title = "pec 1 TRANS")
p_pec_TRANS2 = histogram(RQA_pec2[4,:],title = "pec 2 TRANS")

## ENTR

plot(p_tde_ENTR1, label="data from subsamples")
plot!([RQA_tde_ref1[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde1[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_tde_ENTR2, label="data from subsamples")
plot!([RQA_tde_ref2[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde2[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_ENTR1, label="data from subsamples")
plot!([RQA_GA_ref1[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA1[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_ENTR2, label="data from subsamples")
plot!([RQA_GA_ref2[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA2[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_ENTR1, label="data from subsamples")
plot!([RQA_mdop_ref1[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop1[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_ENTR2, label="data from subsamples")
plot!([RQA_mdop_ref2[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop2[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_ENTR1, label="data from subsamples")
plot!([RQA_pec_ref1[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec1[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_ENTR2, label="data from subsamples")
plot!([RQA_pec_ref2[1]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec2[1]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

## LAM

plot(p_tde_LAM1, label="data from subsamples")
plot!([RQA_tde_ref1[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde1[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_tde_LAM2, label="data from subsamples")
plot!([RQA_tde_ref2[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde2[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_LAM1, label="data from subsamples")
plot!([RQA_GA_ref1[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA1[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_LAM2, label="data from subsamples")
plot!([RQA_GA_ref2[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA2[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_LAM1, label="data from subsamples")
plot!([RQA_mdop_ref1[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop1[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_LAM2, label="data from subsamples")
plot!([RQA_mdop_ref2[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop2[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_LAM1, label="data from subsamples")
plot!([RQA_pec_ref1[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec1[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_LAM2, label="data from subsamples")
plot!([RQA_pec_ref2[2]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec2[2]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

## RTE

plot(p_tde_RTE1, label="data from subsamples")
plot!([RQA_tde_ref1[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde1[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_tde_RTE2, label="data from subsamples")
plot!([RQA_tde_ref2[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde2[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_RTE1, label="data from subsamples")
plot!([RQA_GA_ref1[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA1[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_RTE2, label="data from subsamples")
plot!([RQA_GA_ref2[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA2[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_RTE1, label="data from subsamples")
plot!([RQA_mdop_ref1[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop1[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_RTE2, label="data from subsamples")
plot!([RQA_mdop_ref2[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop2[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_RTE1, label="data from subsamples")
plot!([RQA_pec_ref1[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec1[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_RTE2, label="data from subsamples")
plot!([RQA_pec_ref2[3]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec2[3]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

## TRANS

plot(p_tde_TRANS1, label="data from subsamples")
plot!([RQA_tde_ref1[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde1[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_tde_TRANS2, label="data from subsamples")
plot!([RQA_tde_ref2[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_tde2[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_TRANS1, label="data from subsamples")
plot!([RQA_GA_ref1[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA1[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_GA_TRANS2, label="data from subsamples")
plot!([RQA_GA_ref2[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_GA2[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_TRANS1, label="data from subsamples")
plot!([RQA_mdop_ref1[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop1[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_mdop_TRANS2, label="data from subsamples")
plot!([RQA_mdop_ref2[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_mdop2[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_TRANS1, label="data from subsamples")
plot!([RQA_pec_ref1[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec1[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

plot(p_pec_TRANS2, label="data from subsamples")
plot!([RQA_pec_ref2[4]], seriestype=vline, linecolor="red", linewidth=10, label="reference")
plot!([median_RQA_pec2[4]], seriestype=vline, linecolor="gray", linestyle= :dash, linewidth=10, label="median")

## Deviation of median from reference in σ-units:

diff_RQA_tde1 = zeros(4)
diff_RQA_tde2 = zeros(4)
diff_RQA_GA1 = zeros(4)
diff_RQA_GA2 = zeros(4)
diff_RQA_mdop1 = zeros(4)
diff_RQA_mdop2 = zeros(4)
diff_RQA_pec1 = zeros(4)
diff_RQA_pec2 = zeros(4)

for i = 1:4
    diff_RQA_tde1[i] = abs(median_RQA_tde1[i]-RQA_tde_ref1[i])/std_RQA_tde1[i]
    diff_RQA_tde2[i] = abs(median_RQA_tde2[i]-RQA_tde_ref2[i])/std_RQA_tde2[i]

    diff_RQA_GA1[i] = abs(median_RQA_GA1[i]-RQA_GA_ref1[i])/std_RQA_GA1[i]
    diff_RQA_GA2[i] = abs(median_RQA_GA2[i]-RQA_GA_ref2[i])/std_RQA_GA2[i]

    diff_RQA_mdop1[i] = abs(median_RQA_mdop1[i]-RQA_mdop_ref1[i])/std_RQA_mdop1[i]
    diff_RQA_mdop2[i] = abs(median_RQA_mdop2[i]-RQA_mdop_ref2[i])/std_RQA_mdop2[i]

    diff_RQA_pec1[i] = abs(median_RQA_pec1[i]-RQA_pec_ref1[i])/std_RQA_pec1[i]
    diff_RQA_pec2[i] = abs(median_RQA_pec2[i]-RQA_pec_ref2[i])/std_RQA_pec2[i]
end

rel_diff_RQA_tde1 = zeros(4)
rel_diff_RQA_tde2 = zeros(4)
rel_diff_RQA_GA1 = zeros(4)
rel_diff_RQA_GA2 = zeros(4)
rel_diff_RQA_mdop1 = zeros(4)
rel_diff_RQA_mdop2 = zeros(4)
rel_diff_RQA_pec1 = zeros(4)
rel_diff_RQA_pec2 = zeros(4)

for i = 1:4
    rel_diff_RQA_tde1[i] = abs(median_RQA_tde1[i]-RQA_tde_ref1[i])/RQA_tde_ref1[i]
    rel_diff_RQA_tde2[i] = abs(median_RQA_tde2[i]-RQA_tde_ref2[i])/RQA_tde_ref2[i]

    rel_diff_RQA_GA1[i] = abs(median_RQA_GA1[i]-RQA_GA_ref1[i])/RQA_GA_ref1[i]
    rel_diff_RQA_GA2[i] = abs(median_RQA_GA2[i]-RQA_GA_ref2[i])/RQA_GA_ref2[i]

    rel_diff_RQA_mdop1[i] = abs(median_RQA_mdop1[i]-RQA_mdop_ref1[i])/RQA_mdop_ref1[i]
    rel_diff_RQA_mdop2[i] = abs(median_RQA_mdop2[i]-RQA_mdop_ref2[i])/RQA_mdop_ref2[i]

    rel_diff_RQA_pec1[i] = abs(median_RQA_pec1[i]-RQA_pec_ref1[i])/RQA_pec_ref1[i]
    rel_diff_RQA_pec2[i] = abs(median_RQA_pec2[i]-RQA_pec_ref2[i])/RQA_pec_ref2[i]
end

mn1 = [diff_RQA_tde1[1], diff_RQA_tde1[2], diff_RQA_tde1[3], diff_RQA_tde1[4],
    diff_RQA_GA1[1], diff_RQA_GA1[2], diff_RQA_GA1[3], diff_RQA_GA1[4],
    diff_RQA_mdop1[1], diff_RQA_mdop1[2], diff_RQA_mdop1[3], diff_RQA_mdop1[4],
    diff_RQA_pec1[1], diff_RQA_pec1[2], diff_RQA_pec1[3], diff_RQA_pec1[4]
    ]
mn2 = [diff_RQA_tde2[1], diff_RQA_tde2[2], diff_RQA_tde2[3], diff_RQA_tde2[4],
    diff_RQA_GA2[1], diff_RQA_GA2[2], diff_RQA_GA2[3], diff_RQA_GA2[4],
    diff_RQA_mdop2[1], diff_RQA_mdop2[2], diff_RQA_mdop2[3], diff_RQA_mdop2[4],
    diff_RQA_pec2[1], diff_RQA_pec2[2], diff_RQA_pec2[3], diff_RQA_pec2[4]
    ]
sx = repeat(["Cao's TDE", "Garcia & Almeida", "MDOP", "PECUZAL"], inner = 4)
nam = repeat(["ENTR", "LAM", "RTE", "TRANS"], outer=4)

StatsPlots.groupedbar(nam, mn1, group = sx, ylabel = "deviation from reference [σ]",
        title = "chem. Oscillators Exp.1")

StatsPlots.groupedbar(nam, mn2, group = sx, ylabel = "deviation from reference [σ]",
        title = "chem. Oscillators Exp.2")

total_diff_tde1 = sum(diff_RQA_tde1)
total_diff_tde2 = sum(diff_RQA_tde2)

total_diff_GA1 = sum(diff_RQA_GA1)
total_diff_GA2 = sum(diff_RQA_GA2)

total_diff_mdop1 = sum(diff_RQA_mdop1)
total_diff_mdop2 = sum(diff_RQA_mdop2)

total_diff_pec1 = sum(diff_RQA_pec1)
total_diff_pec2 = sum(diff_RQA_pec2)

rel_mn1 = [rel_diff_RQA_tde1[1], rel_diff_RQA_tde1[2], rel_diff_RQA_tde1[3], rel_diff_RQA_tde1[4],
    rel_diff_RQA_GA1[1], rel_diff_RQA_GA1[2], rel_diff_RQA_GA1[3], rel_diff_RQA_GA1[4],
    rel_diff_RQA_mdop1[1], rel_diff_RQA_mdop1[2], rel_diff_RQA_mdop1[3], rel_diff_RQA_mdop1[4],
    rel_diff_RQA_pec1[1], rel_diff_RQA_pec1[2], rel_diff_RQA_pec1[3], rel_diff_RQA_pec1[4]
    ]
rel_mn2 = [rel_diff_RQA_tde2[1], rel_diff_RQA_tde2[2], rel_diff_RQA_tde2[3], rel_diff_RQA_tde2[4],
    rel_diff_RQA_GA2[1], rel_diff_RQA_GA2[2], rel_diff_RQA_GA2[3], rel_diff_RQA_GA2[4],
    rel_diff_RQA_mdop2[1], rel_diff_RQA_mdop2[2], rel_diff_RQA_mdop2[3], rel_diff_RQA_mdop2[4],
    rel_diff_RQA_pec2[1], rel_diff_RQA_pec2[2], rel_diff_RQA_pec2[3], rel_diff_RQA_pec2[4]
    ]

# save relative deviations for plotting in Matlab
writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/relative_dev_RQA_1.csv", rel_mn1)
writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/relative_dev_RQA_2.csv", rel_mn2)

StatsPlots.groupedbar(nam, rel_mn1, group = sx, ylabel = "rel .deviation from reference",
        title = "chem. Oscillators Exp.1")

StatsPlots.groupedbar(nam, rel_mn2, group = sx, ylabel = "rel .deviation from reference",
        title = "chem. Oscillators Exp.2")

total_rel_diff_tde1 = sum(rel_diff_RQA_tde1)
total_rel_diff_tde2 = sum(rel_diff_RQA_tde2)

total_rel_diff_GA1 = sum(rel_diff_RQA_GA1)
total_rel_diff_GA2 = sum(rel_diff_RQA_GA2)

total_rel_diff_mdop1 = sum(rel_diff_RQA_mdop1)
total_rel_diff_mdop2 = sum(rel_diff_RQA_mdop2)

total_rel_diff_pec1 = sum(rel_diff_RQA_pec1)
total_rel_diff_pec2 = sum(rel_diff_RQA_pec2)

# for plotting 1-rel.diff:
plot_rel_mn1 = [1-rel_diff_RQA_tde1[1], 1-rel_diff_RQA_tde1[2], 1-rel_diff_RQA_tde1[3], 1-rel_diff_RQA_tde1[4],
    1-rel_diff_RQA_GA1[1], 1-rel_diff_RQA_GA1[2], 1-rel_diff_RQA_GA1[3], 1-rel_diff_RQA_GA1[4],
    1-rel_diff_RQA_mdop1[1], 1-rel_diff_RQA_mdop1[2], 1-rel_diff_RQA_mdop1[3], 1-rel_diff_RQA_mdop1[4],
    1-rel_diff_RQA_pec1[1], 1-rel_diff_RQA_pec1[2], 1-rel_diff_RQA_pec1[3], 1-rel_diff_RQA_pec1[4]
    ]
plot_rel_mn2 = [1-rel_diff_RQA_tde2[1], 1-rel_diff_RQA_tde2[2], 1-rel_diff_RQA_tde2[3], 1-rel_diff_RQA_tde2[4],
    1-rel_diff_RQA_GA2[1], 1-rel_diff_RQA_GA2[2], 1-rel_diff_RQA_GA2[3], 1-rel_diff_RQA_GA2[4],
    1-rel_diff_RQA_mdop2[1], 1-rel_diff_RQA_mdop2[2], 1-rel_diff_RQA_mdop2[3], 1-rel_diff_RQA_mdop2[4],
    1-rel_diff_RQA_pec2[1], 1-rel_diff_RQA_pec2[2], 1-rel_diff_RQA_pec2[3], 1-rel_diff_RQA_pec2[4]
    ]

StatsPlots.groupedbar(nam, plot_rel_mn1, group = sx, ylabel = "rel .deviation from reference",
        title = "chem. Oscillators Exp.1")

StatsPlots.groupedbar(nam, plot_rel_mn2, group = sx, ylabel = "rel .deviation from reference",
        title = "chem. Oscillators Exp.2")

# display results:
round.(plot_rel_mn1[5]; digits=4)

max_stars_repo_path: scripts/Experimental data/Standard electrochemical chaos/analyze_standard_ecoscillator_try2.jl