text:
function fornberg_weights!(c::AbstractMatrix{T}, x, ξ) where T
    nd,m = size(c) .- 1
    c = OffsetArray(c, 0:nd, 0:m)
    x = OffsetArray(x, 0:nd)
    length(x) == nd+1 ||
        throw(DimensionMismatch("Number of grid points does not match requested number of derivative coefficients"))

    c₁ = one(T)
    c₄ = x[0] - ξ

    c .= zero(T)
    c[0,0] = 1

    for i = 1:nd
        mn = min(i,m)
        c₂ = one(T)
        c₅ = c₄
        c₄ = x[i] - ξ
        for j = 0:i-1
            c₃ = x[i] - x[j]
            c₂ *= c₃
            for k = mn:-1:1
                c[i,k] = c₁*(k*c[i-1,k-1] - c₅*c[i-1,k])/c₂
            end
            c[i,0] = -c₁*c₅*c[i-1,0]/c₂
            for k = mn:-1:1
                c[j,k] = (c₄*c[j,k] - k*c[j,k-1])/c₃
            end
            c[j,0] = c₄*c[j,0]/c₃
        end
        c₁ = c₂
    end
    c
end

function fornberg_weights(x, ξ::T, m::Integer) where T
    c = zeros(promote_type(eltype(x),T), length(x), m+1)
    fornberg_weights!(c, x, ξ)
    c
end

function fornberg_all_weights!(c::AbstractArray{T,3}, x, ξ) where T
    nd,m = size(c,1) .- 1,size(c,3) .- 1
    c = OffsetArray(c, 0:nd, 0:nd, 0:m)
    x = OffsetArray(x, 0:nd)
    length(x) == nd+1 ||
        throw(DimensionMismatch("Number of grid points does not match requested number of derivative coefficients"))

    c .= zero(T)
    c[0,0,0] = 1

    c₁ = one(T)
    c₄ = x[0] - ξ

    for i = 1:nd
        mn = min(i,m)
        c₂ = one(T)
        c₅ = c₄
        c₄ = x[i] - ξ
        for j = 0:i-1
            c₃ = x[i] - x[j]
            c₂ *= c₃
            c[j,i,0] = c₄*c[j,i-1,0]/c₃
            for k = 1:mn
                c[j,i,k] = (c₄*c[j,i-1,k] - k*c[j,i-1,k-1])/c₃
            end
        end
        c[i,i,0] = -c₁*c₅*c[i-1,i-1,0]/c₂

        for k = 1:mn
            c[i,i,k] = c₁*(k*c[i-1,i-1,k-1] - c₅*c[i-1,i-1,k])/c₂
        end
        c₁ = c₂
    end
    c
end

function fornberg_all_weights(x, ξ::T, m::Integer) where T
    c = zeros(promote_type(eltype(x),T), length(x), length(x), m+1)
    fornberg_all_weights!(c, x, ξ)
    c
end

function fornberg_test_driver(x, ξ)
    pnum(n::AbstractFloat) = printfmt("{1:12.4f}", n)
    pnum(n) = print(format(n, width=12))
    pnum(n::Rational) =
        print(format(isinteger(n) ? Int(n) : n, width=12))

    M = min(4,length(x)-1)
    C = fornberg_all_weights(x, ξ, M)
    N = size(C,2)-1
    print("j ")
    for j = 0:N
        pnum(j)
    end
    println()
    print("x ")
    foreach(pnum, x)
    println()
    println()
    for (ki,k) = enumerate(0:M)
        printfmtln("k = {1:2d}", k)
        for i = k:N
            ii = i+1
            printfmt("{1:2d}", i)
            foreach(pnum, view(C,1:ii,ii,ki))
            println()
        end
        println()
    end
    println()
    C = fornberg_weights(x, ξ, M)
    for (ki,k) = enumerate(0:M)
        printfmt("{1:2d}", k)
        foreach(pnum, view(C,:,ki))
        println(stdout)
    end
end


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

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

function _swaggerinternal_adminConversationsArchive(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.archive", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Archive a public or private channel.
Param: channel_id::String (required)
Return: Dict{String, Any}
"""
function adminConversationsArchive(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsArchive(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsArchive(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsArchive(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsConvertToPrivate(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.convertToPrivate", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Convert a public channel to a private channel.
Param: channel_id::String (required)
Return: Dict{String, Any}
"""
function adminConversationsConvertToPrivate(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsConvertToPrivate(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsConvertToPrivate(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsConvertToPrivate(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsCreate(_api::AdminConversationsApi, name::String, is_private::Bool; description=nothing, org_wide=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.create", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_param(_ctx.form, "description", description)  # type String
    Swagger.set_param(_ctx.form, "is_private", is_private)  # type Bool
    Swagger.set_param(_ctx.form, "org_wide", org_wide)  # type Bool
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Create a public or private channel-based conversation.
Param: name::String (required)
Param: is_private::Bool (required)
Param: description::String
Param: org_wide::Bool
Param: team_id::String
Return: Dict{String, Any}
"""
function adminConversationsCreate(_api::AdminConversationsApi, name::String, is_private::Bool; description=nothing, org_wide=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsCreate(_api, name, is_private; description=description, org_wide=org_wide, team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsCreate(_api::AdminConversationsApi, response_stream::Channel, name::String, is_private::Bool; description=nothing, org_wide=nothing, team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsCreate(_api, name, is_private; description=description, org_wide=org_wide, team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsDelete(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.delete", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Delete a public or private channel.
Param: channel_id::String (required)
Return: Dict{String, Any}
"""
function adminConversationsDelete(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsDelete(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsDelete(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsDelete(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsDisconnectShared(_api::AdminConversationsApi, channel_id::String; leaving_team_ids=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.disconnectShared", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_param(_ctx.form, "leaving_team_ids", leaving_team_ids)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Disconnect a connected channel from one or more workspaces.
Param: channel_id::String (required)
Param: leaving_team_ids::String
Return: Dict{String, Any}
"""
function adminConversationsDisconnectShared(_api::AdminConversationsApi, channel_id::String; leaving_team_ids=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsDisconnectShared(_api, channel_id; leaving_team_ids=leaving_team_ids, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsDisconnectShared(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; leaving_team_ids=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsDisconnectShared(_api, channel_id; leaving_team_ids=leaving_team_ids, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsGetConversationPrefs(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.conversations.getConversationPrefs", ["slackAuth"])
    Swagger.set_param(_ctx.query, "channel_id", channel_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Get conversation preferences for a public or private channel.
Param: channel_id::String (required)
Return: Dict{String, Any}
"""
function adminConversationsGetConversationPrefs(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsGetConversationPrefs(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsGetConversationPrefs(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsGetConversationPrefs(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsGetTeams(_api::AdminConversationsApi, channel_id::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.conversations.getTeams", ["slackAuth"])
    Swagger.set_param(_ctx.query, "channel_id", channel_id)  # type String
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Get all the workspaces a given public or private channel is connected to within this Enterprise org.
Param: channel_id::String (required)
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function adminConversationsGetTeams(_api::AdminConversationsApi, channel_id::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsGetTeams(_api, channel_id; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsGetTeams(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsGetTeams(_api, channel_id; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsInvite(_api::AdminConversationsApi, user_ids::String, channel_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.invite", ["slackAuth"])
    Swagger.set_param(_ctx.form, "user_ids", user_ids)  # type String
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Invite a user to a public or private channel.
Param: user_ids::String (required)
Param: channel_id::String (required)
Return: Dict{String, Any}
"""
function adminConversationsInvite(_api::AdminConversationsApi, user_ids::String, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsInvite(_api, user_ids, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsInvite(_api::AdminConversationsApi, response_stream::Channel, user_ids::String, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsInvite(_api, user_ids, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsRename(_api::AdminConversationsApi, channel_id::String, name::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.rename", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Rename a public or private channel.
Param: channel_id::String (required)
Param: name::String (required)
Return: Dict{String, Any}
"""
function adminConversationsRename(_api::AdminConversationsApi, channel_id::String, name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsRename(_api, channel_id, name; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsRename(_api::AdminConversationsApi, response_stream::Channel, channel_id::String, name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsRename(_api, channel_id, name; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsSearch(_api::AdminConversationsApi; team_ids=nothing, query=nothing, limit=nothing, cursor=nothing, search_channel_types=nothing, sort=nothing, sort_dir=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.conversations.search", ["slackAuth"])
    Swagger.set_param(_ctx.query, "team_ids", team_ids)  # type String
    Swagger.set_param(_ctx.query, "query", query)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "search_channel_types", search_channel_types)  # type String
    Swagger.set_param(_ctx.query, "sort", sort)  # type String
    Swagger.set_param(_ctx.query, "sort_dir", sort_dir)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Search for public or private channels in an Enterprise organization.
Param: team_ids::String
Param: query::String
Param: limit::Int32
Param: cursor::String
Param: search_channel_types::String
Param: sort::String
Param: sort_dir::String
Return: Dict{String, Any}
"""
function adminConversationsSearch(_api::AdminConversationsApi; team_ids=nothing, query=nothing, limit=nothing, cursor=nothing, search_channel_types=nothing, sort=nothing, sort_dir=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsSearch(_api; team_ids=team_ids, query=query, limit=limit, cursor=cursor, search_channel_types=search_channel_types, sort=sort, sort_dir=sort_dir, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsSearch(_api::AdminConversationsApi, response_stream::Channel; team_ids=nothing, query=nothing, limit=nothing, cursor=nothing, search_channel_types=nothing, sort=nothing, sort_dir=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsSearch(_api; team_ids=team_ids, query=query, limit=limit, cursor=cursor, search_channel_types=search_channel_types, sort=sort, sort_dir=sort_dir, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsSetConversationPrefs(_api::AdminConversationsApi, channel_id::String, prefs::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.setConversationPrefs", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_param(_ctx.form, "prefs", prefs)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set the posting permissions for a public or private channel.
Param: channel_id::String (required)
Param: prefs::String (required)
Return: Dict{String, Any}
"""
function adminConversationsSetConversationPrefs(_api::AdminConversationsApi, channel_id::String, prefs::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsSetConversationPrefs(_api, channel_id, prefs; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsSetConversationPrefs(_api::AdminConversationsApi, response_stream::Channel, channel_id::String, prefs::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsSetConversationPrefs(_api, channel_id, prefs; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsSetTeams(_api::AdminConversationsApi, channel_id::String; team_id=nothing, target_team_ids=nothing, org_channel=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.setTeams", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "target_team_ids", target_team_ids)  # type String
    Swagger.set_param(_ctx.form, "org_channel", org_channel)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set the workspaces in an Enterprise grid org that connect to a public or private channel.
Param: channel_id::String (required)
Param: team_id::String
Param: target_team_ids::String
Param: org_channel::Bool
Return: Dict{String, Any}
"""
function adminConversationsSetTeams(_api::AdminConversationsApi, channel_id::String; team_id=nothing, target_team_ids=nothing, org_channel=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsSetTeams(_api, channel_id; team_id=team_id, target_team_ids=target_team_ids, org_channel=org_channel, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsSetTeams(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; team_id=nothing, target_team_ids=nothing, org_channel=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsSetTeams(_api, channel_id; team_id=team_id, target_team_ids=target_team_ids, org_channel=org_channel, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminConversationsUnarchive(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.conversations.unarchive", ["slackAuth"])
    Swagger.set_param(_ctx.form, "channel_id", channel_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Unarchive a public or private channel.
Param: channel_id::String (required)
Return: Dict{String, Any}
"""
function adminConversationsUnarchive(_api::AdminConversationsApi, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsUnarchive(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminConversationsUnarchive(_api::AdminConversationsApi, response_stream::Channel, channel_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminConversationsUnarchive(_api, channel_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminConversationsArchive, adminConversationsConvertToPrivate, adminConversationsCreate, adminConversationsDelete, adminConversationsDisconnectShared, adminConversationsGetConversationPrefs, adminConversationsGetTeams, adminConversationsInvite, adminConversationsRename, adminConversationsSearch, adminConversationsSetConversationPrefs, adminConversationsSetTeams, adminConversationsUnarchive

ext: jl


text:
# Define logistic loss

# Define data for logistic regression problem
struct SparseLogitData{T <: Real} <: AbstractDataset{T}
    # Store X as sparse matrix
    X::AbstractSparseMatrix{T, Int}
    y::Vector{T}
    y_pred::Vector{T}
    # Cache values of exp(y[i] * ypred[i])
    cache_σ::Vector{T}
end
function SparseLogitData(X::AbstractSparseMatrix{T, Int}, y::Vector{T}) where T
    @assert size(X, 1) == size(y, 1)
    return SparseLogitData(X, y, zeros(T, size(y, 1)), zeros(T, size(y, 1)))
end
function SparseLogitData(libsmv_file::String; scale_data=false)
    dataset = parse_libsvm(libsmv_file)
    # Convert to sparse matrix
    X = to_sparse(dataset)
    if scale_data
        scale!(NormalScaler(), X)
    end
    y = copy(dataset.labels)
    format_label!(y)
    return SparseLogitData(X, y, zeros(size(y, 1)),
                           zeros(size(y, 1)))
end

ndata(dat::SparseLogitData) = length(dat.y)
nfeatures(data::SparseLogitData) = size(data.X, 2)

function predict!(data::SparseLogitData{T}, x::AbstractVector{T}, intercept=NaN) where T
    mul!(data.y_pred, data.X, x)
    if isfinite(intercept)
        data.y_pred .+= intercept
    end
    # Update σ
    @inbounds for j in 1:ndata(data)
        data.cache_σ[j] = expit(-data.y_pred[j] * data.y[j])
    end
end

"""
Compute logistic loss of given vector parameter `ω ∈  R^p`.

## Complexity
O(n)

"""
function loss(ω::AbstractVector{T}, data::SparseLogitData{T}) where T
    res = zero(T)
    n = ndata(data)
    # Could not take into account sparsity pattern here
    @inbounds for i in 1:n
        res += loss(Logit(), data.y_pred[i], data.y[i])
    end
    return one(T) / n * res
end

"""
Compute gradient of logistic loss for given vector parameter `ω ∈  R^p`.

## Complexity
O(n * p)

"""
function gradient!(grad::AbstractVector{T}, ω::AbstractVector{T},
                   data::SparseLogitData{T}, fit_intercept::Bool=false) where T
    # Sanity check
    @assert length(ω) == length(grad) == (nfeatures(data) + fit_intercept)
    invn = -one(T) / ndata(data)
    rowsv = SparseArrays.rowvals(data.X)
    xvals = SparseArrays.nonzeros(data.X)
    @inbounds for ncol in 1:nfeatures(data)
        acc = zero(T)
        for j ∈  SparseArrays.nzrange(data.X, ncol)
            i = rowsv[j]
            vals = xvals[j]
            acc += invn * data.y[i] * data.cache_σ[i] * vals
        end
        grad[ncol] += acc
    end

    if fit_intercept
        for j in 1:ndata(data)
            @inbounds grad[end] += invn * data.y[j] * data.cache_σ[j]
        end
    end

    return nothing
end

"""
Compute Hessian of logistic loss for given vector parameter `ω ∈  R^p`.

## Complexity
O(n * p * (p-1) /2)

"""
function hessian!(hess::AbstractVector{T}, ω::AbstractVector{T},
                  data::SparseLogitData{T}, fit_intercept::Bool=false) where T
    error("Currently  not implemented")
end

"""
Compute Hessian-vector product of logistic loss for given vector
parameter `ω ∈  R^p`.

## Complexity
O(n * 2 * p)

"""
function hessvec!(hessvec::AbstractVector{T}, ω::AbstractVector{T},
                  vec::AbstractVector{T}, data::SparseLogitData{T},
                  fit_intercept::Bool=false) where T
    p = nfeatures(data)
    n = ndata(data)
    @assert length(ω) == length(vec) == length(hessvec) == (p + fit_intercept)
    rowsv = rowvals(data.X)
    xvals = nonzeros(data.X)
    invn = one(T) / n

    acc = zeros(n)
    @inbounds for ncol in 1:p
        for j in nzrange(data.X, ncol)
            i = rowsv[j]
            vals = xvals[j]
            cst = invn * data.cache_σ[i] * (1.0 - data.cache_σ[i])
            acc[i] += vec[ncol] * vals * cst
        end
    end
    @inbounds for ncol in 1:p
        @simd for j in nzrange(data.X, ncol)
            i = rowsv[j]
            vals = xvals[j]
            hessvec[ncol] += acc[i] * vals
        end
    end

    if fit_intercept
        # TODO: simplify computation of Hv for intercept
        @inbounds for ncol in 1:p
            @simd for j in nzrange(data.X, ncol)
                i = rowsv[j]
                vals = xvals[j]
                σz = data.cache_σ[i]
                cst = invn * σz * (T(1) - σz)
                hessvec[ncol] += cst * vals * vec[end]
            end
        end
        @inbounds for i in 1:n
            σz = data.cache_σ[i]
            cst = invn * σz * (T(1) - σz)
            hessvec[end] += acc[i] + cst * vec[end]
        end
    end
    return nothing
end

"""
Compute diagonal of Hessian for given vector parameter `ω ∈  R^p`.

## Complexity
O(n * p)

"""
function diaghess!(diagh::AbstractVector{T}, ω::AbstractVector{T},
                   data::SparseLogitData{T}, fit_intercept::Bool=false) where T
    # Sanity check
    n = ndata(data)
    p = nfeatures(data)
    @assert length(ω) == length(diagh) == (p + fit_intercept)
    invn = one(T) / n
    rowsv = rowvals(data.X)
    xvals = nonzeros(data.X)
    @inbounds for ncol in 1:p
        acc = zero(T)
        for j in nzrange(data.X, ncol)
            i = rowsv[j]
            vals = xvals[j]
            σz = data.cache_σ[i]
            acc += invn * σz * (one(T) - σz) * vals^2
        end
        diagh[ncol] = acc
    end

    if fit_intercept
        for j in 1:n
            σz = data.cache_σ[i]
            diagh[end] += invn * σz * (one(T) - σz)
        end
    end
    return nothing
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL
using LinearAlgebra

#include("")
store = []
scaling = .85
V,(VV,EV,FV,CV) = Lar.cuboid([0.25,0.25,0.25],true,[-0.25,-0.25,-0.25])
mybox = (V,CV,FV,EV)
for k=1:10
	size = rand()*scaling
	scale = Lar.s(size,size,size)
	transl = Lar.t(rand(3)...)
	alpha = 2*pi*rand()
	rx = Lar.r(alpha,0,0); ry = Lar.r(0,alpha,0); rz = Lar.r(0,0,alpha)
	rot = rx * ry * rz
	str = Lar.Struct([ transl, scale, rot, mybox ])
	push!(store, Lar.struct2lar(str))
end

str = Lar.Struct(store)
V,CV,FV,EV = Lar.struct2lar(str)
# V = Plasm.normalize3D(V) TODO:  solve MethodError bug

GL.VIEW([ GL.GLPol(V,CV, GL.COLORS[1]) ]);

#preparazione dati
cop_EV = Lar.coboundary_0(EV::Lar.Cells);#stessi dati ma come matrice sparsa
cop_EW = convert(Lar.ChainOp, cop_EV);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');#cobordo 1

FV = Lar.compute_FV( cop_EW, cop_FE )

model = (convert(Lar.Points,W'), FV)

sp_idx = Lar.spaceindex(model)

fs_num = size(cop_FE, 1)

rV = Array{Float64,2}(undef,0,3)

using SparseArrays

rEV = SparseArrays.spzeros(Int8,0,0)

rFE = SparseArrays.spzeros(Int8,0,0)

function pippo(rV,rEV,rFE)
	for sigma in 1:fs_num
		#print(sigma, "/", fs_num, "\r")
		nV, nEV, nFE = Lar.Arrangement.frag_face(W, cop_EW, cop_FE, sp_idx, sigma)
            nV = convert(Lar.Points, nV)
		a,b,c = Lar.skel_merge( rV,rEV,rFE,  nV,nEV,nFE )
		rV=a;  rEV=b;  rFE=c
	end
	return rV, rEV, rFE
end

rV, rEV, rFE = pippo(rV,rEV,rFE)

function merge_vertices(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp, err=1e-4)
    vertsnum = size(V, 1)
    edgenum = size(EV, 1)
    facenum = size(FE, 1)
    newverts = zeros(Int, vertsnum)
    # KDTree constructor needs an explicit array of Float64
    V = Array{Float64,2}(V)
    W = convert(Lar.Points, LinearAlgebra.transpose(V))
    kdtree = KDTree(W)
	# remove vertices congruent to a single representative
    todelete = []
    i = 1
    for vi in 1:vertsnum
        if !(vi in todelete)
            nearvs = Lar.inrange(kdtree, V[vi, :], err)
            newverts[nearvs] .= i
            nearvs = setdiff(nearvs, vi)
            todelete = union(todelete, nearvs)
            i = i + 1
        end
    end
    nV = V[setdiff(collect(1:vertsnum), todelete), :]

    # translate edges to take congruence into account
    edges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    oedges = Array{Tuple{Int, Int}, 1}(undef, edgenum)
    for ei in 1:edgenum
        v1, v2 = EV[ei, :].nzind
        edges[ei] = Tuple{Int, Int}(sort([newverts[v1], newverts[v2]]))
        oedges[ei] = Tuple{Int, Int}(sort([v1, v2]))
    end
    nedges = union(edges)
    # remove edges of zero length
    nedges = filter(t->t[1]!=t[2], nedges)
    nedgenum = length(nedges)
    nEV = spzeros(Int8, nedgenum, size(nV, 1))

    etuple2idx = Dict{Tuple{Int, Int}, Int}()
    for ei in 1:nedgenum
    	begin
        	nEV[ei, collect(nedges[ei])] .= 1
        	nEV
        end
        etuple2idx[nedges[ei]] = ei
    end
    for e in 1:nedgenum
    	v1,v2 = findnz(nEV[e,:])[1]
    	nEV[e,v1] = -1; nEV[e,v2] = 1
    end

    # compute new faces to take congruence into account
    faces = [[
        map(x->newverts[x], FE[fi, ei] > 0 ? oedges[ei] : reverse(oedges[ei]))
        for ei in FE[fi, :].nzind
    ] for fi in 1:facenum]

    visited = []
    function filter_fn(face)

        verts = []
        map(e->verts = union(verts, collect(e)), face)
        verts = Set(verts)

        if !(verts in visited)
            push!(visited, verts)
            return true
        end
        return false
    end

    nfaces = filter(filter_fn, faces)

    nfacenum = length(nfaces)
    nFE = spzeros(Int8, nfacenum, size(nEV, 1))

    for fi in 1:nfacenum
        for edge in nfaces[fi]
            ei = etuple2idx[Tuple{Int, Int}(sort(collect(edge)))]
            nFE[fi, ei] = sign(edge[2] - edge[1])
        end
    end

    return Lar.Points(nV), nEV, nFE
end

rV, rEV, rFE = Lar.Arrangement.merge_vertices(rV, rEV, rFE)

function check_single_non_zero_CSC(A::SparseMatrixCSC)
	sleep_time = 10 #seconds
	temp = 0
	result = []
	for i in range(1,length=size(A,2))
		push!(result, temp)
		temp = 0
		for j in range(1, length=size(A,1))
	        temp+=abs(A[j,i])
		end
	#println(temp)#show column sums
	end
	count = 0
	for k in result
		if isodd(k)
			count+=1
		end
	end
	if count!=0
		println("//////////////////////////////////////////////////////////////////////////////////////////////////////////")
		println("CHECKING SINGLE NON ZERO ELEMENT")
		println("WARNING, THE MATRIX CONTAINS $count COLUMN WITH A SINGLE NON ZERO ELEMENT")
		println("the execution will sleep for $sleep_time seconds")
		println("//////////////////////////////////////////////////////////////////////////////////////////////////////////")
		sleep(sleep_time)
	end
end

bicon_comps = Lar.Arrangement.biconnected_components(rEV)

check_single_non_zero_CSC(rFE)

rV2, rEV2, rFE2, rCF2 = Lar.Arrangement.spatial_arrangement_2(rV, rEV, rFE)

triangulated_faces = Lar.triangulate2D(rV2, [rEV2, rFE2])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
V3 = convert(Lar.Points, rV2')
GL.VIEW( GL.GLExplode(V3,FVs,1.5,1.5,1.5,99,1) );

EVs = Lar.FV2EVs(rEV2, rFE2) # polygonal face fragments
GL.VIEW( GL.GLExplode(V3,EVs,1.5,1.5,1.5,99,1) );


text:
# Interface

###########################################################################################################
##### Methods common to all: AbstractLinearQuadraticModel, AbstractConicModel, AbstractNonlinearModel #####
##### see: http://mathprogbasejl.readthedocs.io/en/latest/solverinterface.html                        #####
###########################################################################################################

getobjgap(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl!")

getrawsolver(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl!")

getsolvetime(m::SCIPMathProgModel) = ccall((:SCIPgetSolvingTime, libcsip), Cdouble, (Ptr{Void},), _getInternalSCIP(m))

getsense(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl!")

numvar(m::SCIPMathProgModel) = _getNumVars(m)

"The number of proper constraints, excluding those from lazy callbacks."
numconstr(m::SCIPMathProgModel) = _getNumConss(m)

freemodel!(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl!")

# TODO: mapping for :SemiCont, :SemiInt
const vartypemap = Dict{Symbol, Cint}(
  :Cont => 3,
  :Bin => 0,
  :Int => 1
)
const revvartypemap = Dict{Cint, Symbol}(
  3 => :Cont,
  0 => :Bin,
  1 => :Int
)

function getvartype(m::SCIPMathProgModel)
    nvars = _getNumVars(m)
    types = zeros(nvars)
    for idx = one(Cint):nvars
        types[idx] = _getVarType(m, idx - one(Cint))
    end
    return map(x -> revvartypemap[x], types)
end

function setvartype!(m::SCIPMathProgModel, vartype::Vector{Symbol})
    nvars = Cint(length(vartype))
    scipvartypes = map(vt -> vartypemap[vt], vartype)
    for idx = one(Cint):nvars
        _chgVarType(m, idx - one(Cint), scipvartypes[idx])
    end
end

optimize!(m::SCIPMathProgModel) = _solve(m)

function status(m::SCIPMathProgModel)
    statusmap = [:Optimal,
                 :Infeasible,
                 :Unbounded,
                 :InfeasibleOrUnbounded,
                 :UserLimit, # node limit
                 :UserLimit, # time limit
                 :UserLimit, # memory limit
                 :UserLimit, # user limit
                 :Unknown    # TODO: find good value
                 ]
    stat = _getStatus(m)
    return statusmap[stat + 1]
end

getobjbound(m::SCIPMathProgModel) = _getObjBound(m)

getobjval(m::SCIPMathProgModel) = _getObjValue(m)

function getsolution(m::SCIPMathProgModel)
    nvars = _getNumVars(m)
    values = zeros(nvars)
    _getVarValues(m, values)
    values
end

function setsense!(m::SCIPMathProgModel, sense)
    if sense == :Max
        _setSenseMaximize(m)
    else
        _setSenseMinimize(m)
    end
end

###########################################################################
##### Methods specific to AbstractLinearQuadraticModel                #####
##### see: http://mathprogbasejl.readthedocs.io/en/latest/lpqcqp.html #####
###########################################################################

loadproblem!(m::SCIPMathProgModel, filename::AbstractString) =
    error("Not implemented for SCIP.jl")

function loadproblem!(m::SCIPMathProgModel, A, varlb, varub, obj, rowlb, rowub, sense)
    # TODO: clean old model?

    nrows, ncols = size(A)
    nvars = Cint(ncols)
    varindices = collect(zero(Cint):nvars - one(Cint))

    for v in 1:ncols
        # TODO: define enum for vartype?
        _addVar(m, float(varlb[v]), float(varub[v]),
                Cint(3), Ptr{Cint}(C_NULL))
    end
    for c in 1:nrows
        # TODO: care about sparse matrices
        denserow = float(collect(A[c, :]))
        _addLinCons(m, nvars, varindices, denserow,
                    float(rowlb[c]), float(rowub[c]), Ptr{Cint}(C_NULL))
    end

    _setObj(m, nvars, varindices, float(obj))

    if sense == :Max
        _setSenseMaximize(m)
    else
        _setSenseMinimize(m)
    end
end

writeproblem(m::SCIPMathProgModel, filename::AbstractString) =
    error("Not implemented for SCIP.jl")

getvarLB(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

function setvarLB!(m::SCIPMathProgModel, lb)
    indices = collect(zero(Cint):Cint(numvar(m) - 1))
    _chgVarLB(m, Cint(numvar(m)), indices, float(lb))
end

getvarUB(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

function setvarUB!(m::SCIPMathProgModel, ub)
    indices = collect(zero(Cint):Cint(numvar(m) - 1))
    _chgVarUB(m, Cint(numvar(m)), indices, float(ub))
end

getconstrLB(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

# setconstrLB!(m::SCIPMathProgModel, lb) = error("Not implemented for SCIP.jl")

getconstrUB(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

# setconstrUB!(m::SCIPMathProgModel, ub) = error("Not implemented for SCIP.jl")

getobj(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

function setobj!(m::SCIPMathProgModel, obj)
    indices = collect(zero(Cint):Cint(numvar(m) - 1))
    _setObj(m, Cint(numvar(m)), indices, float(obj))
end

getconstrmatrix(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

addvar!(m::SCIPMathProgModel, constridx, constrcoef, l, u, objcoef) =
    error("Not implemented for SCIP.jl")

addvar!(m::SCIPMathProgModel, l, u, objcoef) = error("Not implemented for SCIP.jl")

function addconstr!(m::SCIPMathProgModel, varidx, coef, lb, ub)
    _addLinCons(m, Cint(length(varidx)), Vector{Cint}(varidx - 1), float(coef),
                float(lb), float(ub), Ptr{Cint}(C_NULL))
end

numlinconstr(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getconstrsolution(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getreducedcosts(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getconstrduals(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getinfeasibilityray(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getbasis(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getunboundedray(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getsimplexiter(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getbarrieriter(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

function setwarmstart!(m::SCIPMathProgModel, v)
    # does not support incomplete solutions (with NaN)
    _setInitialSolution(m, float(v))
end

##########################################################################
##### Methods specific to Integer Programming                        #####
##########################################################################

getnodecount(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

function addsos1!(m::SCIPMathProgModel, idx, weight)
    nidx = Cint(length(idx))
    cidx = convert(Vector{Cint}, idx - 1)
    _addSOS1(m, nidx, cidx, weight, Ptr{Cint}(C_NULL))
end

function addsos2!(m::SCIPMathProgModel, idx, weight)
    nidx = Cint(length(idx))
    cidx = convert(Vector{Cint}, idx - 1)
    _addSOS2(m, nidx, cidx, weight, Ptr{Cint}(C_NULL))
end

##########################################################################
##### Methods specific to Quadratic Programming                      #####
##########################################################################

numquadconstr(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

setquadobj!{T<:Real}(m::SCIPMathProgModel, Q::Array{T, 2}) =
    error("Not implemented for SCIP.jl")

setquadobj!(m::SCIPMathProgModel, rowidx, colidx, quadval) =
    error("Not implemented for SCIP.jl")

setquadobjterms!(m::SCIPMathProgModel, rowidx, colidx, quadval) =
    _setQuadObj(m, Cint(0), Array{Cint}(0), Array{Cdouble}(0),
                 Cint(length(rowidx)), convert(Vector{Cint}, rowidx - 1),
                 convert(Vector{Cint}, colidx - 1), quadval)

function addquadconstr!(m::SCIPMathProgModel, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
    clhs = -Inf
    crhs =  Inf
    if sense == '<'
        crhs = rhs
    elseif sense == '>'
        clhs = rhs
    else
        @assert sense == '='
        clhs = rhs
        crhs = rhs
    end
    _addQuadCons(m, Cint(length(linearidx)), convert(Vector{Cint}, linearidx - 1),
                 linearval, Cint(length(quadrowidx)), convert(Vector{Cint}, quadrowidx - 1),
                 convert(Vector{Cint}, quadcolidx - 1), quadval, clhs, crhs, Ptr{Cint}(C_NULL))
end

getquadconstrsolution(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getquadconstrduals(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getquadinfeasibilityray(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

getquadconstrRHS(m::SCIPMathProgModel) = error("Not implemented for SCIP.jl")

setquadconstrRHS!(m::SCIPMathProgModel, lb) = error("Not implemented for SCIP.jl")

##########################################################################
##### Methods specific to MIP Callbacks                              #####
##########################################################################

# use a different type for heuristic callback and multiple dispatch to implements
# the methods that they share
abstract SCIPCallbackData <: MathProgCallbackData

type SCIPLazyCallbackData <: SCIPCallbackData
    model::SCIPMathProgModel
    csip_lazydata::Ptr{Void}
end

# this is the function that should fit the CSIP_LAZYCALLBACK signature
function lazycb_wrapper(csip_model::Ptr{Void}, csip_lazydata::Ptr{Void},
                        userdata::Ptr{Void})
    # m, f = unsafe_pointer_to_objref(userdata)::(SCIPMathProgModel, Function)
    # WTF: TypeError: typeassert: expected Type{T}, got Tuple{DataType,DataType}
    m, f = unsafe_pointer_to_objref(userdata)
    d = SCIPLazyCallbackData(m, csip_lazydata)
    ret = f(d)
    ret == :Exit && _interrupt(m)

    return convert(Cint, 0) # CSIP_RETCODE_OK
end

function setlazycallback!(m::SCIPMathProgModel, f)
    # f is function(d::SCIPLazyCallbackData)

    cbfunction = cfunction(lazycb_wrapper, Cint, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    m.lazy_userdata = (m, f)

    _addLazyCallback(m, cbfunction, m.lazy_userdata)
end

# if we are called from a lazy callback, we check whether the LP relaxation is integral
function cbgetstate(d::SCIPLazyCallbackData)
    context = _lazyGetContext(d.csip_lazydata)
    mapping = [:MIPNode, :MIPSol, :Other]
    mapping[context + 1]
end

function cbgetlpsolution(d::SCIPLazyCallbackData, output)
    _lazyGetVarValues(d.csip_lazydata, output)
end

#TODO: what should we do if there is no best solution?
cbgetmipsolution(d::SCIPLazyCallbackData, output) = cbgetlpsolution(d, output)

function cbgetlpsolution(d::SCIPCallbackData)
    output = Array(Float64, _getNumVars(m))
    cbgetlpsolution(d, output)
end

#TODO: what should we do if there is no best solution?
cbgetmipsolution(d::SCIPLazyCallbackData) = cbgetlpsolution(d)

function cbaddlazy!(d::SCIPLazyCallbackData, varidx, varcoef, sense, rhs)
    clhs = -Inf
    crhs =  Inf
    if sense == '<'
        crhs = rhs
    elseif sense == '>'
        clhs = rhs
    else
        @assert sense == '='
        clhs = rhs
        crhs = rhs
    end
    _lazyAddLinCons(d.csip_lazydata, convert(Cint, length(varidx)),
                  convert(Vector{Cint}, varidx - 1), varcoef, clhs, crhs,
                  convert(Cint, 0))
end

# function cbaddlazylocal!(d::SCIPLazyCallbackData, varidx, varcoef, sense, rhs)
# end

type SCIPHeurCallbackData <: SCIPCallbackData
    model::SCIPMathProgModel
    csip_heurdata::Ptr{Void}
    sol::Vector{Float64}
end

# this is the function that should fit the CSIP_HEURCALLBACK signature
function heurcb_wrapper(csip_model::Ptr{Void}, csip_heurdata::Ptr{Void},
                        userdata::Ptr{Void})
    # m, f = unsafe_pointer_to_objref(userdata)::(SCIPMathProgModel, Function)
    # WTF: TypeError: typeassert: expected Type{T}, got Tuple{DataType,DataType}
    m, f = unsafe_pointer_to_objref(userdata)
    d = SCIPHeurCallbackData(m, csip_heurdata, fill(NaN, numvar(m)))

    ret = f(d)
    ret == :Exit && _interrupt(m)

    return convert(Cint, 0) # CSIP_RETCODE_OK
end

function setheuristiccallback!(m::SCIPMathProgModel, f)
    # f is function(d::SCIPHeurCallbackData)

    cbfunction = cfunction(heurcb_wrapper, Cint, (Ptr{Void}, Ptr{Void}, Ptr{Void}))
    m.heur_userdata = (m, f)

    _addHeuristicCallback(m, cbfunction, m.heur_userdata)
end

# TODO: detect :MIPSol like with lazy constraints?
cbgetstate(d::SCIPHeurCallbackData) = :MIPNode

function cbgetlpsolution(d::SCIPHeurCallbackData, output)
    _heurGetVarValues(d.csip_heurdata, output)
end

function cbaddsolution!(d::SCIPHeurCallbackData)
    # check for unspecified values (NaN)
    complete = (findfirst(v -> v === NaN, d.sol) == 0)

    if complete
        # add solution that was filled from cbsetsolutionvalue
        _heurAddSolution(d.csip_heurdata, d.sol)
    else
        warn("Incomplete solutions not supported, skipping candidate.")
    end

    # reset solution vector for the next solution
    d.sol = fill(NaN, numvar(d.model))
end

function cbsetsolutionvalue!(d::SCIPHeurCallbackData, varidx, value)
    d.sol[varidx] = value
end

##########################################################################
##### Methods specific to AbstractConicModel                         #####
##### see: http://mathprogbasejl.readthedocs.io/en/latest/conic.html #####
##########################################################################

########################################################################
##### Methods specific to AbstractNonlinear                        #####
##### see: http://mathprogbasejl.readthedocs.io/en/latest/nlp.html #####
########################################################################
const julia_op_to_csip_op = Dict{Symbol, Int}(
   :- => 9,
   :/ => 11,
   :sqrt => 13,
   :^ => 14,
   :exp => 17,
   :log => 18,
   :+ => 64,
   :* => 65
)

#inspired in ReverseDiffSparse's conversion.jl
immutable CSIPNodeData
    nodetype::Int
    childids::Array{Int}
end

# for backward compatibility with 0.4 (from JuMP/macros.jl)
function comparison_to_call(ex)
    if Meta.isexpr(ex,:comparison) && length(ex.args) == 3
        return Expr(:call, ex.args[2], ex.args[1], ex.args[3])
    else
        return ex
    end
end

# parse constraint expression into something easier for CSIP
function constr_expr_to_nodedata(ex::Expr)
    values = Float64[]
    csipex = CSIPNodeData[]

    # for compatibility with 0.4
    ex = comparison_to_call(ex)

    # a ex.args is [side, comp, expr, comp, side] or [expr, comp, side]
    if length(ex.args) == 5 # two sided constraint
        expr_to_csip(ex.args[3], values, csipex)
    else
        @assert length(ex.args) == 3
        expr_to_csip(ex.args[2], values, csipex)
    end
    return csipex, values
end

# parse objective expression into something easier for CSIP
function obj_expr_to_nodedata(ex::Expr)
    values = Float64[]
    csipex = CSIPNodeData[]
    expr_to_csip(ex, values, csipex)
    return csipex, values
end

# store expression tree as an array of CSIPNodeData
# Each element of a CSIPNodeData consists of the operator identifier
# and an array with the position of the operator's children.
# For each node in the expression tree, first process its children
# (creating an array with their position) and afterwards add create a new
# node with the node's operator and the position of its children.
function expr_to_csip(ex::Expr, values, csipex)
    if Meta.isexpr(ex,:call) # functional operator
        # create node
        childids = []
        op = ex.args[1]
        nchildren = length(ex.args)
        for c in 2:nchildren  # the first actual children is in positin 2
            pos = expr_to_csip(ex.args[c], values, csipex)
            push!(childids, pos)
        end
        push!(csipex, CSIPNodeData(julia_op_to_csip_op[op], childids))

    elseif Meta.isexpr(ex, :ref) # variable
        @assert ex.args[1] == :x
        push!(csipex, CSIPNodeData(1, [ex.args[2]]))

    else # not supported I guess
        error("Unrecognized expression $ex: $(ex.head), probably not supported")
    end
    return length(csipex)
end

# parsing values
function expr_to_csip(ex::Number, values, csipex)
    valueidx = length(values)+1
    push!(values,ex)
    push!(csipex, CSIPNodeData(2, [valueidx]))
    return length(csipex)
end

# the AbstractNLPEvaluator contains the constraints and objective.
# One can get different data from it (with initialize)
# SCIP doesn't need most of the data, just the epression
function loadproblem!(m::SCIPMathProgModel, numVars, numConstr,
                      l, u, lb, ub, sense, d::AbstractNLPEvaluator)
    # add variables
    for v in 1:numVars
        _addVar(m, float(l[v]), float(u[v]), Cint(3), Ptr{Cint}(C_NULL))
    end

    # we want to get the expressions
    initialize(d, [:ExprGraph])

    # add constraints
    for c in 1:numConstr
        csipex, values = constr_expr_to_nodedata(constr_expr(d, c))
        children = Cint[]
        beg = Cint[1]
        ops = Cint[]
        for node in csipex
            ops = [ops; node.nodetype]
            children = [children; node.childids]
            push!(beg, beg[end]+length(node.childids))
        end
        _addNonLinCons(m, Cint(length(ops)),
                       convert(Vector{Cint},ops),
                       convert(Vector{Cint},children - 1),
                       convert(Vector{Cint},beg - 1),
                       values, float(lb[c]), float(ub[c]),
                       Ptr{Cint}(C_NULL))
    end

    # add objective
    csipex, values = obj_expr_to_nodedata(obj_expr(d))
    children = Cint[]
    beg = Cint[1]
    ops = Cint[]
    for node in csipex
        ops = [ops; node.nodetype]
        children = [children; node.childids]
        push!(beg, beg[end]+length(node.childids))
    end
    _setNonlinearObj(m, Cint(length(ops)),
                     convert(Vector{Cint},ops),
                     convert(Vector{Cint},children - 1),
                     convert(Vector{Cint},beg - 1),
                     values)

    # set sense
    if sense == :Max
        _setSenseMaximize(m)
    else
        _setSenseMinimize(m)
    end
end


text:
package java.io;

import java.io.InputStream;

public class HighInputStream extends InputStream {
    public HighInputStream() {}
    public int read() { return ({H})42; }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using ControlSystems
"""
Example for designing an LQG speed controller for an electrical DC motor.
"""

# Constants
Ke                      = 0.006156           # electromotive force constant in V/rpm
Kt                      = 0.0728             # Torque constant (Nm/A)
J                       = 2.8*700e-6;        # Inertia of motor plus load (kgm^2)
Rel                     = 0.11;              # DC motor resistance (Ω)
L                       = 34e-6;             # DC motor inductance (H)

# helper to convert a scalar into an 1x1 Matrix
mat(sc) = reshape([sc],1,1)

# Create an instance of a DC motor model in state space form.
# Ke: electromotive force constant in V/rpm
# Kt: torque constant (Nm/A)
# L: Inductance; R: Resistance, J: Inertia; b: viscous friction coefficient
function motor(Ke, Kt, L, R, J, b=1e-3)
    Ke1 = Ke*60.0/2π   # change unit to rad/s
    A = [-b/J     Kt/J
         -Ke1/L   -R/L];
    B = [0
        1/L];
    C = [1     0];
    D = 0;
    sys = ss(A, B, C, D);
end

p60 = motor(Ke, Kt, L, Rel, J)
f1 = stepplot(p60, 1)
f2 = bodeplot(p60)

# LQR control
Q = [1.     0;
     0      1]

R = 20.
K = lqr(p60.A, p60.B, Q, R)
# needs to be modified if Nbar is not a scalar
Nbar = 1. ./ (p60.D - (p60.C - p60.D*K) * inv(p60.A - p60.B*K) * p60.B)

# Kalman filter based observer
Vd = [10.     0   # covariance of the speed estimation
     0       10]; # covariance of the current estimation
Vn = 0.04;       # covariance for the speed measurement (radians/s)^2
Ka = kalman(p60, Vd, mat(Vn))
C = ControlSystems.observer_controller(p60, K, Ka)

Gcl = let (A,B,C,D) = ssdata(p60)
     Acl = [A-B*K B*K; zero(A) A-Ka*C]
     Bcl = [B * Nbar; zero(B)]
     Ccl = [C zero(C)]
     ss(Acl, Bcl, Ccl, 0)
end

T = Gcl
S = 1-T

# 1000 logarithmically spaced values from -3 to 3
f3 = bodeplot([Gcl, S, T], exp10.(range(-3, stop=3, length=1000)))
f4 = stepplot(Gcl, 1, label="Closed loop system using LQG")
Plots.plot(f1, f2, f3, f4, layout=(2,2), size=(800, 600))


text:
module TestDescribe

using Test
using RLEVectors

@testset begin

    x = RLEVector([4, 5, 6], [3, 6, 9])

    # nrun
    @test nrun(x) == 3

    # length
    @test length(x) == 9

    # size
    @test size(x) == (9,)
    @test size(x, 1) == 9
    @test size(x, 2) == (9, 1)

    # isempty
    @test isempty(x) == false
    @test isempty(RLEVector(Int[], Int[])) == true

    # ==
    @test x == x

    # isequal
    @test isequal(x, x)

    # eltype, endtype
    @test endtype(x) == Int64
    @test eltype(x) == Int64

end # testset

end # module


text:
using OffsetArrays
using SpecialFunctions: gamma, erfcx
using Polynomials
using PyPlot
using Printf

T = Float64
α = 1 / parse(T, "2")
β = one(T)
a = OffsetArray{T}(undef, 0:4)
b = OffsetArray{T}(undef, 0:4)
a[0] = 1 / gamma(β)
b[0] = zero(T)
let
pow = one(T)
for k = 1:4
    pow = -pow
    a[k] = pow / gamma(β+k*α)
    b[k] = pow * sinpi(k*α-β) * gamma(1+k*α-β) / π
end
end

A_43 = [ 1 0 -a[0]    0     0
      0 1 -a[1] -a[0]    0
      0 0 -a[2] -a[1] -a[0]
      1 0 -b[0] -b[1] -b[2]
      0 1    0  -b[0] -b[1] ]

rhs_43 = [ a[1], a[2], a[3], 0, 0 ]
pq = A_43 \ rhs_43

p = OffsetArray{T}(undef,0:2)
q = OffsetArray{T}(undef,0:3)
p[0] = a[0]
p[1:2] = pq[1:2]
q[0] = one(T)
q[1:3] = pq[3:5]

P_43 = Polynomial(p[0:2])
Q_43 = Polynomial(q[0:3])
R_43(x) = P_43(x) / Q_43(x)

A_52 = [ 1 0 -a[0]    0     0
         0 1 -a[1] -a[0]    0
         0 0 -a[2] -a[1] -a[0]
         0 0 -a[3] -a[2] -a[1]
         0 1    0  -b[0] -b[1] ]

rhs_52 = [ a[1], a[2], a[3], a[4], 0 ]
pq = A_52 \ rhs_52

p = OffsetArray{T}(undef,0:3)
q = OffsetArray{T}(undef,0:3)
p[0] = a[0]
p[1:2] = pq[1:2]
p[3] = zero(T)
q[0] = one(T)
q[1:3] = pq[3:5]

P_52 = Polynomial(p[0:2])
Q_52 = Polynomial(q[0:3])
R_52(x) = P_52(x) / Q_52(x)

figure(1)
x = range(0, 10, length=201)
Eα = erfcx.(x)
plot(x, R_43.(x)-Eα, x, R_52.(x)-Eα)
grid(true)
legend((L"$m=4, n=3$", L"m=5, n=2"))
xlabel(L"$x$", fontsize=12)
savefig("pade_43_52.pdf")

for k = 1:3
    @printf("%4d& %8.5f& %8.5f& %8.5f& %8.5f\\\\\n", k-1,
            P_43.coeffs[k], Q_43.coeffs[k], P_52.coeffs[k], Q_52.coeffs[k])
end
@printf("%4d& %8.5f& %8.5f& %8.5f& %8.5f\n", 3,
        0.0, Q_43.coeffs[4], 0.0, Q_52.coeffs[4])


text:
module Elsa

using Requires

# function __init__()
#     # @require Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("diagonalizers/pardiso.jl")
#     # @require KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" include("diagonalizers/krylovkit.jl")
#     @require Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" include("diagonalizers/arpack.jl")
#     # @require ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" include("diagonalizers/arnoldimethod.jl")
#     # @require IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" include("diagonalizers/iterativesolvers.jl")
# end

using StaticArrays, NearestNeighbors, SparseArrays, LinearAlgebra, OffsetArrays,
      ProgressMeter, LinearMaps, Random

using SparseArrays: getcolptr, AbstractSparseMatrix

export sublat, bravais, lattice, dims, sites, supercell, unitcell,
       hopping, onsite, onsiteselector, hoppingselector, onsite!, hopping!,
       hamiltonian, parametric, bloch, bloch!, optimize!, similarmatrix,
       flatten, wrap, transform!, combine,
       spectrum, bandstructure, marchingmesh, defaultmethod, bands, vertices,
       energies, states,
       momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM

export LatticePresets, RegionPresets, HamiltonianPresets

export LinearAlgebraPackage, ArpackPackage, KrylovKitPackage

export @SMatrix, @SVector, SMatrix, SVector

export ishermitian, I

const NameType = Symbol
const nametype = Symbol

const TOOMANYITERS = 10^8

include("iterators.jl")
include("presets.jl")
include("lattice.jl")
include("model.jl")
include("hamiltonian.jl")
include("parametric.jl")
include("mesh.jl")
include("diagonalizer.jl")
include("bandstructure.jl")
include("KPM.jl")
include("convert.jl")
include("tools.jl")

end

max_issues_repo_name: pablosanjose/QBox

max_forks_repo_name: pablosanjose/QBox


text:
using Cards
using Test

# write your own tests here
hand = Hand([2♣, 3♣, 6♣, 7♣, 8♣, 9♣, 2♢, 3♢, 4♢, 7♢, 10♢, J♢, A♢, 4♡, 5♡, 6♡, 7♡, Q♡, K♡, A♡, 4♠, 6♠, 9♠, K♠, A♠])

@testset "Card in Hand" begin
    @test 2♣ in hand
    @test !(4♣ in hand)
    @test !(A♣ in hand)
    @test A♠ in hand
end

@testset "Suit intersetion with Hand" begin
    @test ♡ ∩ hand == Hand([4♡, 5♡, 6♡, 7♡, Q♡, K♡, A♡])
    @test ♠ ∩ hand == Hand([4♠, 6♠, 9♠, K♠, A♠])
    @test length(♣ ∩ hand) == 6
end


text:
for (f_jl, f_np) in libdevice_1
    isdefined(Base, f_jl) || continue
    @eval Base.$f_jl(a::AbstractJaxScalar{T}) where T = numpy.$f_np(a)
end

for (f_jl, f_np) in libdevice_2
    isdefined(Base, f_jl) || continue
    @eval Base.$f_jl(a::AbstractJaxScalar{T}, b::Number) where T = numpy.$f_np(a, b)
    @eval Base.$f_jl(a::Number, b::AbstractJaxScalar{T}) where T = numpy.$f_np(a, b)
    #@eval Base.$f_jl(a::AbstractJaxScalar{T}, b::AbstractJaxScalar{T2}) where {T,T2} = numpy.$f_np(a, b)
end

for op in overridenbfuncs
    bop = Symbol(".", op)
    @eval begin
        $op(a::AbstractJaxArray{T,0}, b) where T = $bop(a, b)
        $op(a, b::AbstractJaxArray{T,0}) where T = $bop(a, b)
        $op(a::AbstractJaxArray{T1,0},
                b::AbstractJaxArray{T2,0}) where {T1,T2} = $bop(a, b)

        $op(a::AbstractJaxArray{T,0}, b::Number) where T = $bop(a, b)
        $op(a::Number, b::AbstractJaxArray{T,0}) where T = $bop(a, b)
    end
end


text:
using Plots
using MOONS.Grids
FT = Float64
a, b, n = FT(0.0), FT(1.0), 10
n_βs = 4
R = 1 .+ (10 .^ range(FT(-3), stop=FT(1), length=n_βs))
y_buff = 2

p = plot()
for (y,β) in enumerate(R)
  c = Coordinates(a,b,n; warpfun=Roberts_both, args=(β,))
  plot!(c.n.h, zeros(length(c.n.h)) .+ y,
    markershape=:circle, label="β=$β",
    ylims=(0,n_βs+y_buff), yaxis=nothing)
  plot!(layout=(n_βs,1), title="Roberts stretching (both)")
end
savefig("Roberts_both.svg")

p = plot()
for (y,β) in enumerate(R)
  c = Coordinates(a,b,n; warpfun=Roberts_left, args=(β,))
  plot!(c.n.h, zeros(length(c.n.h)) .+ y,
    markershape=:circle, label="β=$β",
    ylims=(0,n_βs+y_buff), yaxis=nothing)
  plot!(layout=(n_βs,1), title="Roberts stretching (left)")
end
savefig("Roberts_left.svg")

p = plot()
for (y,β) in enumerate(R)
  c = Coordinates(a,b,n; warpfun=Roberts_right, args=(β,))
  plot!(c.n.h, zeros(length(c.n.h)) .+ y,
    markershape=:circle, label="β=$β",
    ylims=(0,n_βs+y_buff), yaxis=nothing)
  plot!(layout=(n_βs,1), title="Roberts stretching (right)")
end
savefig("Roberts_right.svg")


text:
# Used for Dijisktra's algorithm
module PriorityQueue

import Base:
    ∈, length, minimum,
    sizehint!, push!, take!, insert!

export RankedQueue
export update!

"""
    parent(i)

Return the index of the parent of node `i`.
"""
parent(i) = i÷2

"""
    left(i)

Return the index of the left child of node `i`.
"""
left(i)   = 2*i
"""
    right(i)

Return the index of the right child of node `i`.
"""
right(i)  = 2*i+1

"""
    struct RankedQueue{T <: Real, S <: Any}
        rank :: Array{T, 1}
        data :: Array{S, 1}
    end

Maintains a priority queue of `data`.
Each datum has a `rank` that determines it's priority in the queue.
`rank` and `data` are sorted in ascending order.
"""
struct RankedQueue{T <: Real, S <: Any}
    rank :: Array{T, 1}
    data :: Array{S, 1}
end

∈(x::S, q::RankedQueue{T,S}) where {T <: Real, S <: Any} = x ∈ q.data

function sizehint!(q::RankedQueue, n)
    sizehint!(q.rank, n)
    sizehint!(q.data, n)
end

"""
    rotateup!(q::RankedQueue, i)

Modify the ranked queue by pushing up the node at index `i` until the priority is sorted again.
"""
function rotateup!(q::RankedQueue, i)
    i == 1 && return i

    while i > 1 && (p=parent(i); q.rank[i] < q.rank[p])
        q.rank[i], q.rank[p] = q.rank[p], q.rank[i]
        q.data[i], q.data[p] = q.data[p], q.data[i]
        i = p
    end

    return i
end

"""
    rotateup!(q::RankedQueue)

Modify the ranked queue by pushing up the last element until the priority is sorted.
"""
rotateup!(q::RankedQueue) = rotateup!(q, length(q.rank))

"""
    rotatedown!(q::RankedQueue, i)

Modify the ranked queue by pushing down the node at index `i` until the priority is sorted.
"""
function rotatedown!(q::RankedQueue, i)
    left(i) > length(q) && return i

    child(i) = (right(i) > length(q) || q.rank[left(i)] < q.rank[right(i)]) ? left(i) : right(i)
    while left(i) <= length(q) && (c=child(i); q.rank[i] > q.rank[c])
        q.rank[i], q.rank[c] = q.rank[c], q.rank[i]
        q.data[i], q.data[c] = q.data[c], q.data[i]
        i = c
    end

    return i
end

"""
    rotatedown!(q::RankedQueue)

Modify the ranked queue by pushing down the root until the priority is sorted.
"""
rotatedown!(q::RankedQueue) = rotatedown!(q, 1)

function RankedQueue(X::Tuple{S, T}...) where {T <: Real, S <: Any}
    q = RankedQueue{T,S}([x[2] for x in X], [x[1] for x in X])
    i = length(X) ÷ 2
    while i >= 1
        rotatedown!(q, i)
        i -= 1
    end

    return q
end

length(q::RankedQueue) = length(q.rank)

# external methods
# XXX: there be dragons - we don't check for duplicated data being passed to us

minimum(q::RankedQueue) = (data=q.data[1], rank=q.rank[1])

"""
    insert!(q::RankedQueue{T}, data::S, rank::T) where {T <: Real, S <: Any}

Push a new element `data` with priority `rank` onto the ranked queue `q`.
Rotates the queue until priority is sorted in ascending order.
"""
function insert!(q::RankedQueue{T}, data::S, rank::T) where {T <: Real, S <: Any}
    push!(q.rank, rank)
    push!(q.data, data)

    rotateup!(q)
end

"""
    take!(q::RankedQueue)

Pop off the element with element with lowest rank/highest priority.
"""
function take!(q::RankedQueue)
    r, q.rank[1] = q.rank[1], q.rank[end]
    pop!(q.rank)

    d, q.data[1] = q.data[1], q.data[end]
    pop!(q.data)

    rotatedown!(q)

    return (data=d, rank=r)
end

"""
    update!(q::RankedQueue{T, S}, data::S, new::T) where {T <: Real, S <: Any}

Change the priority of element `data` to rank `new`.
Will panic if `data` is not contained in queue `q`.
"""
function update!(q::RankedQueue{T, S}, data::S, new::T) where {T <: Real, S <: Any}
    (i = findfirst(d -> d == data, q.data)) == nothing && panic("attempting to update a non-existent data value")

    old = q.rank[i]
    q.rank[i] = new

    return if new > old
               rotatedown!(q, i)
           elseif new < old
               rotateup!(q, i)
           else
               i
           end
end

function test()
    Q = RankedQueue((0,1), (1,2), (2,10), (3,23), (4,0))
    @show Q
    insert!(Q, 5, -1)
    @show Q

    nothing
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
using Klio
using Sockets

Klio.settings.server_host = IPv4(0)
Klio.settings.expl_sqlite_file = joinpath(kliodir, "db/expl.sqlite")

Klio.run()


text:
using Base.Test
import WaterData

@testset "Constant verification" begin
    @test isapprox(WaterData.R_h2o, WaterData.R / WaterData.h2o_molar_mass, rtol=0.01)
end


text:
# Correlation Functions
# #+HTML: Source
# #+HTML:

# [[file:~/.julia/dev/MatrixProductStates/README.org::*Correlation%20Functions][Correlation Functions:1]]
"""
    two_point_correlator((i, op_i)::Pair{Int, Matrix}, (j, op_j)::Pair{Int, Matrix}, L)

Create an MPO on `L` sites (with bond dimension 1) representing identity operators everywhere except
sites `i` and `j` where `op_i` and `op_j` are inserted instead. ie.

    𝟙 ⊗ 𝟙 ⊗ ... ⊗ op_i ⊗ 𝟙 ⊗ ... ⊗ op_j ⊗ 𝟙 ⊗ ... ⊗ 𝟙

example: spin-spin correlation function

we can construct ⟨σᶻᵢσᶻⱼ⟩ on a 12 site lattice as
    σᶻ = [1 0; 0 -1]
    two_point_correlator(i=>σᶻ, j=>σᶻ, 12)
"""
function two_point_correlator((i, op_i), (j, op_j), L)
    d = size(op_i)[1]
    @assert (size(op_i) == (d, d)) && (size(op_j) == (d, d))
    @assert i in 1:L
    @assert j in 1:L
    id = diagm(0 => ones(Complex{Float64}, d))

    op_i_tnsr = reshape(convert(Matrix{Complex{Float64}}, op_i), 1, 1, d, d)
    op_j_tnsr = reshape(convert(Matrix{Complex{Float64}}, op_j), 1, 1, d, d)
    id_tnsr   = reshape(id, 1, 1, d, d)

    tensors = map(1:L) do l
        O_tnsr = (l == i ? op_i_tnsr :
                  l == j ? op_j_tnsr :
                  id_tnsr)
    end
    MPO{L,Complex{Float64}}(tensors)
end
# Correlation Functions:1 ends here


text:
function initialize_frequency_estimator!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{C, O, IC, DC, PSY.KauraPLL, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    F <: PSY.Filter,
}
    Vr_filter = get_inner_vars(dynamic_device)[Vr_filter_var]
    Vi_filter = get_inner_vars(dynamic_device)[Vi_filter_var]

    #Get parameters
    pll_control = PSY.get_freq_estimator(dynamic_device)
    kp_pll = PSY.get_kp_pll(pll_control)
    ki_pll = PSY.get_ki_pll(pll_control)

    #Get initial guess
    θ0_pll = angle(Vr_filter + Vi_filter * 1im)
    Vpll_d0 = Vr_filter
    Vpll_q0 = 0.0
    ϵ_pll0 = 0.0

    function f!(out, x)
        vpll_d = x[1]
        vpll_q = x[2]
        ϵ_pll = x[3]
        θ_pll = x[4]

        V_dq_pll = ri_dq(θ_pll + pi / 2) * [Vr_filter; Vi_filter]

        out[1] = (V_dq_pll[d] - vpll_d)
        out[2] = (V_dq_pll[q] - vpll_q)
        out[3] = atan(vpll_q / vpll_d)
        out[4] = (kp_pll * atan(vpll_q / vpll_d) + ki_pll * ϵ_pll)
    end

    x0 = [Vpll_d0, Vpll_q0, ϵ_pll0, θ0_pll]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in PLL failed")
    else
        sol_x0 = sol.zero

        #Obtain indices for component w/r to device
        local_ix = get_local_state_ix(dynamic_device, PSY.KauraPLL)

        #Update guess of PLL states
        pll_states = @view device_states[local_ix]
        pll_states[1] = sol_x0[1]
        pll_states[2] = sol_x0[2]
        pll_states[3] = sol_x0[3]
        pll_states[4] = sol_x0[4]

        #Update guess of frequency estimator
        get_inner_vars(dynamic_device)[ω_freq_estimator_var] = PSY.get_ω_ref(dynamic_device)
        get_inner_vars(dynamic_device)[θ_freq_estimator_var] = sol_x0[4]
    end
end

function initialize_frequency_estimator!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{C, O, IC, DC, PSY.ReducedOrderPLL, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    F <: PSY.Filter,
}
    Vr_filter = get_inner_vars(dynamic_device)[Vr_filter_var]
    Vi_filter = get_inner_vars(dynamic_device)[Vi_filter_var]

    #Get parameters
    pll_control = PSY.get_freq_estimator(dynamic_device)
    kp_pll = PSY.get_kp_pll(pll_control)
    ki_pll = PSY.get_ki_pll(pll_control)

    #Get initial guess
    θ0_pll = angle(Vr_filter + Vi_filter * 1im)
    Vpll_q0 = 0.0
    ϵ_pll0 = 0.0

    function f!(out, x)
        vpll_q = x[1]
        ϵ_pll = x[2]
        θ_pll = x[3]

        V_dq_pll = ri_dq(θ_pll + pi / 2) * [Vr_filter; Vi_filter]

        out[1] = V_dq_pll[q] - vpll_q
        out[2] = vpll_q
        out[3] = kp_pll * vpll_q + ki_pll * ϵ_pll
    end

    x0 = [Vpll_q0, ϵ_pll0, θ0_pll]
    sol = NLsolve.nlsolve(f!, x0)
    if !NLsolve.converged(sol)
        @warn("Initialization in PLL failed")
    else
        sol_x0 = sol.zero

        #Obtain indices for component w/r to device
        local_ix = get_local_state_ix(dynamic_device, PSY.ReducedOrderPLL)

        #Update guess of PLL states
        pll_states = @view device_states[local_ix]
        pll_states[1] = sol_x0[1]
        pll_states[2] = sol_x0[2]
        pll_states[3] = sol_x0[3]

        #Update guess of frequency estimator
        get_inner_vars(dynamic_device)[ω_freq_estimator_var] = PSY.get_ω_ref(dynamic_device)
        get_inner_vars(dynamic_device)[θ_freq_estimator_var] = sol_x0[3]
    end
end

function initialize_frequency_estimator!(
    device_states,
    static::PSY.StaticInjection,
    dynamic_device::PSY.DynamicInverter{C, O, IC, DC, PSY.FixedFrequency, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    F <: PSY.Filter,
}
    #Get parameters
    pll_control = PSY.get_freq_estimator(dynamic_device)
    frequency = PSY.get_frequency(pll_control)

    #Update guess of frequency estimator
    get_inner_vars(dynamic_device)[ω_freq_estimator_var] = frequency
end


text:
### LombScargle.jl ---  Perform Lomb-Scargle periodogram
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is MIT "Expat".
#
### Code:

__precompile__()

module LombScargle

using Measurements

export lombscargle, power, freq, freqpower, findmaxpower, findmaxfreq,
prob, probinv, fap, fapinv

# This is similar to Periodogram type of DSP.Periodograms module, but for
# unevenly spaced frequencies.
immutable Periodogram{T<:AbstractFloat}
    power::AbstractVector{T}
    freq::AbstractVector{T}
    times::AbstractVector{T}
    norm::AbstractString
end

include("extirpolation.jl")
include("utils.jl")

# Original algorithm that doesn't take into account uncertainties and doesn't
# fit the mean of the signal.  This is implemented following the recipe by
# * Townsend, R. H. D. 2010, ApJS, 191, 247 (URL:
#   http://dx.doi.org/10.1088/0067-0049/191/2/247,
#   Bibcode: http://adsabs.harvard.edu/abs/2010ApJS..191..247T)
function _lombscargle_orig{R1<:Real,R2<:Real,R3<:Real}(times::AbstractVector{R1},
                                                       signal::AbstractVector{R2},
                                                       freqs::AbstractVector{R3},
                                                       center_data::Bool)
    P_type = promote_type(float(R1), float(R2))
    P = Vector{P_type}(freqs)
    N = length(signal)
    nil = zero(P_type)
    # If "center_data" keyword is true, subtract the mean from each point.
    if center_data
        X = signal - mean(signal)
    else
        X = signal
    end
    XX = X⋅X

    @inbounds for n in eachindex(freqs)
        ω = 2pi*freqs[n]
        XC = XS = CC = CS = nil
        for j in eachindex(times)
            ωt = ω*times[j]
            C = cos(ωt)
            S = sin(ωt)
            XC += X[j]*C
            XS += X[j]*S
            CC += C*C
            CS += C*S
        end
        SS      = N - CC
        τ       = 0.5*atan2(2CS, CC - SS)/ω
        c_τ     = cos(ω*τ)
        s_τ     = sin(ω*τ)
        c_τ2    = c_τ*c_τ
        s_τ2    = s_τ*s_τ
        cs_τ_CS = 2c_τ*s_τ*CS
        P[n] = (abs2(c_τ*XC + s_τ*XS)/(c_τ2*CC + cs_τ_CS + s_τ2*SS) +
                abs2(c_τ*XS - s_τ*XC)/(c_τ2*SS - cs_τ_CS + s_τ2*CC))/XX
    end
    return P
end

# Generalised Lomb-Scargle algorithm: this takes into account uncertainties and
# fit the mean of the signal.  This is implemented following the recipe by
# * Zechmeister, M., Kürster, M. 2009, A&A, 496, 577  (URL:
#   http://dx.doi.org/10.1051/0004-6361:200811296,
#   Bibcode: http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)
# In addition, some tricks suggested by
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
#   http://dx.doi.org/10.1086/167197,
#   Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# to make computation faster are adopted.
function _generalised_lombscargle{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::AbstractVector{R1},
                                                                       signal::AbstractVector{R2},
                                                                       w::AbstractVector{R3},
                                                                       freqs::AbstractVector{R4},
                                                                       center_data::Bool,
                                                                       fit_mean::Bool)
    P_type = promote_type(float(R1), float(R2), float(R3), float(R4))
    P = Vector{P_type}(freqs)
    nil = zero(P_type)
    # If "center_data" keyword is true, subtract the mean from each point.
    if center_data || fit_mean
        y = signal - mean(signal)
    else
        y = signal
    end
    YY = w⋅(y.^2)
    if fit_mean
        Y  = w⋅y
    end

    @inbounds for n in eachindex(freqs)
        ω = 2pi*freqs[n]

        # Find τ for current angular frequency
        C = S = CS = CC = nil
        for i in eachindex(times)
            ωt  = ω*times[i]
            W   = w[i]
            c   = cos(ωt)
            s   = sin(ωt)
            CS += W*c*s
            CC += W*c*c
            if fit_mean
                C  += W*c
                S  += W*s
            end
        end
        if fit_mean
            CS -= C*S
            SS  = 1.0 - CC - S*S
            CC -= C*C
        else
            SS  = 1.0 - CC
        end
        τ   = 0.5*atan2(2CS, CC - SS)/ω

        # Now we can compute the power
        YC_τ = YS_τ = CC_τ = nil
        for i in eachindex(times)
            ωt    = ω*(times[i] - τ)
            W     = w[i]
            c     = cos(ωt)
            s     = sin(ωt)
            YC_τ += W*y[i]*c
            YS_τ += W*y[i]*s
            CC_τ += W*c*c
        end
        if fit_mean
            # "C_τ" and "S_τ" are computed following equation (7) of Press &
            # Rybicki, this formula simply comes from angle difference
            # trigonometric identities.
            cos_ωτ = cos(ω*τ)
            sin_ωτ = sin(ω*τ)
            C_τ    = C*cos_ωτ + S*sin_ωτ
            S_τ    = S*cos_ωτ - C*sin_ωτ
            YC_τ  -= Y*C_τ
            YS_τ  -= Y*S_τ
            SS_τ   = 1.0 - CC_τ - S_τ*S_τ
            CC_τ  -= C_τ*C_τ
            YY    -= Y*Y
        else
            SS_τ  = 1.0 - CC_τ
        end
        P[n] = (abs2(YC_τ)/CC_τ + abs2(YS_τ)/SS_τ)/YY
    end
    return P
end

# Fast, but approximate, method to compute the Lomb-Scargle periodogram for
# evenly spaced data.  See
# * Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
#   http://dx.doi.org/10.1086/167197,
#   Bibcode: http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
# This is adapted from Astropy implementation of the method.  See
# `lombscargle_fast' function.
function _press_rybicki{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::Range{R1},
                                                             signal::AbstractVector{R2},
                                                             w::AbstractVector{R3},
                                                             freqs::Range{R4},
                                                             center_data::Bool,
                                                             fit_mean::Bool,
                                                             oversampling::Int,
                                                             Mfft::Int)
    P_type = promote_type(float(R1), float(R2))
    P = Vector{P_type}(freqs)
    nil = zero(P_type)
    # If "center_data" keyword is true, subtract the mean from each point.
    if center_data || fit_mean
        y = signal - w⋅signal
    else
        y = signal
    end

    df = step(freqs)
    N  = length(freqs)
    f0 = minimum(freqs)
    #---------------------------------------------------------------------------
    # 1. compute functions of the time-shift tau at each frequency
    Ch, Sh = trig_sum(times, w .* y, df, N, f0, 1, oversampling, Mfft)
    C2, S2 = trig_sum(times, w,      df, N, f0, 2, oversampling, Mfft)
    if fit_mean
        C, S = trig_sum(times, w, df, N, f0, 1, oversampling, Mfft)
        tan_2ωτ = (S2 .- 2.0 * S .* C) ./ (C2 .- (C .* C .- S .* S))
    else
        tan_2ωτ = S2 ./ C2
    end
    # This is what we're computing below; the straightforward way is slower and
    # less stable, so we use trig identities instead
    #
    # ωτ = 0.5 * atan(tan_2ωτ)
    # S2w, C2w = sin(2 * ωτ), cos(2 * ωτ)
    # Sw, Cw = sin(ωτ), cos(ωτ)
    C2w = 1./(sqrt(1.0 + tan_2ωτ .* tan_2ωτ))
    S2w = tan_2ωτ .* C2w
    Cw = sqrt(0.5 * (1.0 + C2w))
    Sw = sqrt(0.5) .* sign(S2w) .* sqrt(1.0 - C2w)
    #---------------------------------------------------------------------------
    # 2. Compute the periodogram, following Zechmeister & Kurster
    #    and using tricks from Press & Rybicki.
    YY = w⋅(y.^2)
    YC = Ch .* Cw .+ Sh .* Sw
    YS = Sh .* Cw .- Ch .* Sw
    CC = 0.5 * (1.0 + C2 .* C2w .+ S2 .* S2w)
    SS = 0.5 * (1.0 - C2 .* C2w .- S2 .* S2w)
    if fit_mean
        CC -= (C .* Cw .+ S .* Sw) .^ 2
        SS -= (S .* Cw .- C .* Sw) .^ 2
    end
    P = (YC .* YC ./ CC .+ YS .* YS ./ SS)/YY
    return P
end

# This is the switch to select the appropriate function to run
function _lombscargle{R1<:Real,R2<:Real,R3<:Real,R4<:Real}(times::AbstractVector{R1},
                                                           signal::AbstractVector{R2},
                                                           floatrange::Bool,
                                                           with_errors::Bool,
                                                           w::AbstractVector{R3}=ones(signal)/length(signal);
                                                           normalization::AbstractString="standard",
                                                           noise_level::Real=1.0,
                                                           center_data::Bool=true,
                                                           fit_mean::Bool=true,
                                                           oversampling::Integer=5,
                                                           Mfft::Integer=4,
                                                           samples_per_peak::Integer=5,
                                                           nyquist_factor::Integer=5,
                                                           minimum_frequency::Real=NaN,
                                                           maximum_frequency::Real=NaN,
                                                           frequencies::AbstractVector{R4}=
                                                           autofrequency(times,
                                                                         samples_per_peak=samples_per_peak,
                                                                         nyquist_factor=nyquist_factor,
                                                                         minimum_frequency=minimum_frequency,
                                                                         maximum_frequency=maximum_frequency),
                                                           fast::Bool=(length(frequencies) > 200))
    @assert length(times) == length(signal) == length(w)
    # By default, we will use the fast algorithm only if there are more than 200
    # frequencies.  In any case, times must have been passed as a Range.
    if floatrange && fast
        P = _press_rybicki(times, signal, w, frequencies, center_data,
                              fit_mean, oversampling, Mfft)
    else
        if fit_mean || with_errors
            P = _generalised_lombscargle(times, signal, w, frequencies,
                                            center_data, fit_mean)
        else
            P = _lombscargle_orig(times, signal, frequencies, center_data)
        end
    end

    N = length(signal)
    # Normalize periodogram
    if normalization == "standard"
    elseif normalization == "model"
        P = P./(1.0 - P)
    elseif normalization == "log"
        P = -log(1.0 - P)
    elseif normalization == "psd"
        P *= 0.5*N*(w⋅signal.^2)
    elseif normalization == "Scargle"
        P /= noise_level
    elseif normalization == "HorneBaliunas"
        P *= 0.5*(N - 1.0)
    elseif normalization == "Cumming"
        P *= 0.5*(N - 3.0)/(1.0 - maximum(P))
    else
        error("normalization \"", normalization, "\" not supported")
    end

    return Periodogram(P, frequencies, times, normalization)
end

# No uncertainties
lombscargle{R1<:Real,R2<:Real}(times::Range{R1},
                               signal::AbstractVector{R2};
                               kwargs...) =
                                   _lombscargle(times,
                                                signal,
                                                true,
                                                false;
                                                kwargs...)

lombscargle{R1<:Real,R2<:Real}(times::AbstractVector{R1},
                               signal::AbstractVector{R2};
                               kwargs...) =
                                   _lombscargle(times,
                                                signal,
                                                false,
                                                false;
                                                kwargs...)

# Uncertainties provided
function _lombscargle_with_errors{R1<:Real,R2<:Real,R3<:Real}(times::AbstractVector{R1},
                                                              signal::AbstractVector{R2},
                                                              errors::AbstractVector{R3},
                                                              floatrange::Bool;
                                                              kwargs...)
    # Compute weights vector
    w = 1.0./errors.^2
    w /= sum(w)
    return _lombscargle(times, signal, floatrange, true, w; kwargs...)
end

function lombscargle{R1<:Real,R2<:Real,R3<:Real}(times::Range{R1},
                                                 signal::AbstractVector{R2},
                                                 errors::AbstractVector{R3};
                                                 kwargs...)
    return _lombscargle_with_errors(times, signal, errors, true; kwargs...)
end

function lombscargle{R1<:Real,R2<:Real,R3<:Real}(times::AbstractVector{R1},
                                                 signal::AbstractVector{R2},
                                                 errors::AbstractVector{R3};
                                                 kwargs...)
    return _lombscargle_with_errors(times, signal, errors, false; kwargs...)
end

# Uncertainties provided via Measurement type
function lombscargle{T<:Real,F<:AbstractFloat}(times::Range{T},
                                               signal::AbstractVector{Measurement{F}};
                                               kwargs...)
    return _lombscargle_with_errors(times,
                                    value(signal),
                                    uncertainty(signal),
                                    true;
                                    kwargs...)
end

function lombscargle{T<:Real,F<:AbstractFloat}(times::AbstractVector{T},
                                               signal::AbstractVector{Measurement{F}};
                                               kwargs...)
    return _lombscargle_with_errors(times,
                                    value(signal),
                                    uncertainty(signal),
                                    false;
                                    kwargs...)
end

"""
    lombscargle(times::AbstractVector{Real}, signal::AbstractVector{Real},
                errors::AbstractVector{Real}=ones(signal);
                normalization::AbstractString="standard",
                noise_level::Real=1.0,
                center_data::Bool=true,
                fit_mean::Bool=true,
                fast::Bool=true,
                oversampling::Integer=5,
                Mfft::Integer=4,
                samples_per_peak::Integer=5,
                nyquist_factor::Integer=5,
                minimum_frequency::Real=NaN,
                maximum_frequency::Real=NaN,
                frequencies::AbstractVector{Real}=
                autofrequency(times,
                              samples_per_peak=samples_per_peak,
                              nyquist_factor=nyquist_factor,
                              minimum_frequency=minimum_frequency,
                              maximum_frequency=maximum_frequency))

Compute the Lomb-Scargle periodogram of the `signal` vector, observed at
`times`.  You can also specify the uncertainties for each signal point with in
`errors` argument.  All these vectors must have the same length.

Optional keywords arguments are:

* `normalization`: how to normalize the periodogram.  Valid choices are:
  `"standard"`, `"model"`, `"log"`, `"psd"`, `"Scargle"`, `"HorneBaliunas"`,
  `"Cumming"`
* `noise_level`: the noise level used to normalize the periodogram when
  `normalization` is set to `"Scargle"`
* `fit_mean`: if `true`, fit for the mean of the signal using the Generalised
  Lomb-Scargle algorithm (see Zechmeister & Kürster paper below).  If this is
  `false` and no uncertainty on the signal is provided, the original algorithm
  by Lomb and Scargle will be employed (see Townsend paper below)
* `center_data`: if `true`, subtract the mean of `signal` from `signal` itself
  before performing the periodogram.  This is especially important if `fit_mean`
  is `false`
* `frequencies`: the frequecy grid (not angular frequencies) at which the
  periodogram will be computed, as a vector.  If not provided, it is
  automatically determined with `LombScargle.autofrequency` function, which see.
  See below for other available keywords that can be used to affect the
  frequency grid without directly setting `frequencies`
* `fast`: whether to use the fast method by Press & Rybicki, overriding the
  default choice.  In any case, `times` must be a `Range` object in order to use
  this method
* `oversampling`: oversampling the frequency factor for the approximation;
  roughly the number of time samples across the highest-frequency sinusoid.
  This parameter contains the tradeoff between accuracy and speed.  Used only
  when the fast method is employed
* `Mfft`: the number of adjacent points to use in the FFT approximation.  Used
  only when the fast method is employed

In addition, you can use all optional keyword arguments of
`LombScargle.autofrequency` function in order to tune the `frequencies` vector
without calling the function:

* `samples_per_peak`: the approximate number of desired samples across the
  typical peak
* `nyquist_factor`: the multiple of the average Nyquist frequency used to choose
  the maximum frequency if `maximum_frequency` is not provided
* `minimum_frequency`: if specified, then use this minimum frequency rather than
  one chosen based on the size of the baseline
* `maximum_frequency`: if specified, then use this maximum frequency rather than
  one chosen based on the average Nyquist frequency

If the signal has uncertainties, the `signal` vector can also be a vector of
`Measurement` objects (from
[`Measurements.jl`](https://github.com/giordano/Measurements.jl) package), in
which case you don’t need to pass a separate `errors` vector for the
uncertainties of the signal.  See `Measurements.jl` manual at
http://measurementsjl.readthedocs.io/ for details on how to create a vector of
`Measurement` objects.

The algorithm used here are reported in the following papers:

* Press, W. H., Rybicki, G. B. 1989, ApJ, 338, 277 (URL:
  http://dx.doi.org/10.1086/167197, Bibcode:
  http://adsabs.harvard.edu/abs/1989ApJ...338..277P)
* Townsend, R. H. D. 2010, ApJS, 191, 247 (URL:
  http://dx.doi.org/10.1088/0067-0049/191/2/247,
  Bibcode: http://adsabs.harvard.edu/abs/2010ApJS..191..247T)
* Zechmeister, M., Kürster, M. 2009, A&A, 496, 577  (URL:
  http://dx.doi.org/10.1051/0004-6361:200811296,
  Bibcode: http://adsabs.harvard.edu/abs/2009A%26A...496..577Z)
"""
lombscargle

end # module


text:
# Autogenerated wrapper script for libsqm_jll for x86_64-apple-darwin-libgfortran5
export libsqm

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("libsqm")
JLLWrappers.@declare_library_product(libsqm, "@rpath/libsqm.dylib")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsqm,
        "lib/libsqm.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    Send(buf, dest::Integer, tag::Integer, comm::Comm)

Perform a blocking send from the buffer `buf` to MPI rank `dest` of communicator `comm`
using the message tag `tag`.

# External links
$(_doc_external("MPI_Send"))
"""
function Send(buf::Buffer, dest::Integer, tag::Integer, comm::Comm)
    # int MPI_Send(const void* buf, int count, MPI_Datatype datatype, int dest,
    #              int tag, MPI_Comm comm)
    @mpichk ccall((:MPI_Send, libmpi), Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm),
          buf.data, buf.count, buf.datatype, dest, tag, comm)
    return nothing
end
Send(arr::Union{Ref,AbstractArray}, dest::Integer, tag::Integer, comm::Comm) =
    Send(Buffer(arr), dest, tag, comm)

"""
    Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T

Complete a blocking send of `obj` to MPI rank `dest` of communicator `comm`
using with the message tag `tag`.
"""
function Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T
    buf = Ref{T}(obj)
    Send(buf, dest, tag, comm)
end

"""
    send(obj, dest::Integer, tag::Integer, comm::Comm)

Complete a blocking send using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.
"""
function send(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Send(buf, dest, tag, comm)
end

"""
    Isend(data, dest::Integer, tag::Integer, comm::Comm)

Starts a nonblocking send of `data` to MPI rank `dest` of communicator `comm` using with
the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which [`Buffer_send`](@ref) is defined.

Returns the [`Request`](@ref) object for the nonblocking send.

# External links
$(_doc_external("MPI_Isend"))
"""
function Isend(buf::Buffer, dest::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Isend(const void* buf, int count, MPI_Datatype datatype, int dest,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Isend, libmpi), Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf.data, buf.count, buf.datatype, dest, tag, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return req
end
Isend(data, dest::Integer, tag::Integer, comm::Comm) =
    Isend(Buffer_send(data), dest, tag, comm)

"""
    isend(obj, dest::Integer, tag::Integer, comm::Comm)

Starts a nonblocking send of using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.

Returns the commication `Request` for the nonblocking send.
"""
function isend(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Isend(buf, dest, tag, comm)
end

"""
    data = Recv!(recvbuf, src::Integer, tag::Integer, comm::Comm, nothing)
    data, status = Recv!(recvbuf, src::Integer, tag::Integer, comm::Comm, MPI.Status)

Completes a blocking receive into the buffer `recvbuf` from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

`recvbuf` can be a [`Buffer`](@ref), or any object for which `Buffer(recvbuf)` is defined.

Returns the [`Status`](@ref) of the receive.

# See also
- [`Recv`](@ref)
- [`recv`](@ref)

# External links
$(_doc_external("MPI_Recv"))
"""
function Recv!(recvbuf::Buffer, src::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status},Nothing})
    # int MPI_Recv(void* buf, int count, MPI_Datatype datatype, int source,
    #              int tag, MPI_Comm comm, MPI_Status *status)
    @mpichk ccall((:MPI_Recv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{Status}),
                  recvbuf.data, recvbuf.count, recvbuf.datatype, src, tag, comm, something(status, MPI_STATUS_IGNORE))
    return recvbuf.data
end
Recv!(recvbuf, src::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status},Nothing}) =
    Recv!(Buffer(recvbuf), src, tag, comm, status)
function Recv!(recvbuf, src::Integer, tag::Integer, comm::Comm, ::Type{Status})
    status = Ref{Status}()
    data = Recv!(recvbuf, src, tag, comm, status)
    return data, status[]
end

"""
    data = Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm, nothing)
    data, status = Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm, MPI.Status)

Completes a blocking receive of a single object of type `T` from MPI rank `src` of
communicator `comm` using with the message tag `tag`.

Returns a tuple of the object of type `T` and the [`Status`](@ref) of the receive.

# See also
- [`Recv!`](@ref)
- [`recv`](@ref)

# External links
$(_doc_external("MPI_Recv"))
"""
function Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status}, Nothing}) where T
    data = Recv!(Ref{T}(), src, tag, comm, status)
    return data[]
end
function Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm, ::Type{Status}) where T
    status = Ref{Status}()
    val = Recv(T, src, tag, comm, status)
    return val, status[]
end

"""
    recv(src::Integer, tag::Integer, comm::Comm)

Completes a blocking receive of a serialized object from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

Returns the deserialized object and the [`Status`](@ref) of the receive.
"""
function recv(src::Integer, tag::Integer, comm::Comm, status::Union{Ref{Status}, Nothing})
    stat = Probe(src, tag, comm)
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    return MPI.deserialize(buf)
end
function recv(src::Integer, tag::Integer, comm::Comm, ::Type{Status})
    status = Ref{Status}()
    val = recv(src, tag, comm, status)
    return val, status[]
end

"""
    Irecv!(data, src::Integer, tag::Integer, comm::Comm)

Starts a nonblocking receive into the buffer `data` from MPI rank `src` of communicator
`comm` using with the message tag `tag`.

`data` can be a [`Buffer`](@ref), or any object for which `Buffer(data)` is defined.

Returns the [`Request`](@ref) for the nonblocking receive.

# External links
$(_doc_external("MPI_Irecv"))
"""
function Irecv!(buf::Buffer, src::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Irecv(void* buf, int count, MPI_Datatype datatype, int source,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Irecv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf.data, buf.count, buf.datatype, src, tag, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return req
end
Irecv!(data, src::Integer, tag::Integer, comm::Comm) =
    Irecv!(Buffer(data), src, tag, comm)

function irecv(src::Integer, tag::Integer, comm::Comm)
    (flag, stat) = Iprobe(src, tag, comm)
    if !flag
        return (false, nothing, nothing)
    end
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (true, MPI.deserialize(buf), stat)
end

"""
    data = Sendrecv!(
            sendbuf, dest::Integer, sendtag::Integer,
            recvbuf, source::Integer, recvtag::Integer,
            comm::Comm, nothing)
    data, status = Sendrecv!(
            sendbuf, dest::Integer, sendtag::Integer,
            recvbuf, source::Integer, recvtag::Integer,
            comm::Comm, MPI.Status)

Complete a blocking send-receive operation over the MPI communicator `comm`. Send
`sendcount` elements of type `sendtype` from `sendbuf` to the MPI rank `dest` using message
tag `tag`, and receive `recvcount` elements of type `recvtype` from MPI rank `source` into
the buffer `recvbuf` using message tag `tag`. Return a [`Status`](@ref) object.

If not provided, `sendtype`/`recvtype` and `sendcount`/`recvcount` are derived from the
element type and length of `sendbuf`/`recvbuf`, respectively.

# External links
$(_doc_external("MPI_Sendrecv"))
"""
function Sendrecv!(sendbuf::Buffer, dest::Integer, sendtag::Integer,
                   recvbuf::Buffer, source::Integer, recvtag::Integer,
                   comm::Comm, status::Union{Ref{Status}, Nothing})
    # int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, int dest,   int sendtag,
    #                        void *recvbuf, int recvcount, MPI_Datatype recvtype, int source, int recvtag,
    #                    MPI_Comm comm, MPI_Status *status)
    @mpichk ccall((:MPI_Sendrecv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPI_Comm, Ptr{Status}),
                  sendbuf.data, sendbuf.count, sendbuf.datatype, dest, sendtag,
                  recvbuf.data, recvbuf.count, recvbuf.datatype, source, recvtag,
                  comm, status)
    return recvbuf.data
end
Sendrecv!(sendbuf, dest::Integer, sendtag::Integer, recvbuf, source::Integer, recvtag::Integer, comm::Comm, status::Union{Ref{Status}, Nothing}) =
    Sendrecv!(Buffer(sendbuf), dest, sendtag, Buffer(recvbuf), source, recvtag, comm, status)

function Sendrecv!(sendbuf, dest::Integer, sendtag::Integer, recvbuf, source::Integer, recvtag::Integer, comm::Comm, ::Type{Status})
    status = Ref{Status}()
    data = Sendrecv!(sendbuf, dest, sendtag, recvbuf, source, recvtag, comm, status)
    return data, status[]
end