text:
"""
    SubsetEstimates

BLB linear mixed model estimates from one subset.
"""
struct SubsetEstimates{T <: LinearAlgebra.BlasReal}
    # blb parameter
    n_boots::Int64
    # dimension of parameter estimates, for pre-allocation
    p::Int64
    q::Int64
    # blb results
    βs::Matrix{T}
    Σs::Array{T}
    σ²s::Vector{T}
end

"""
    SubsetEstimates(n_boots, p, q)

Constructor for `SubsetEstimates`.
"""
function SubsetEstimates(n_boots::Int64, p::Int64, q::Int64)
    # vector of vector/matrices approach
    # βs  = [Vector{Float64}(undef,  p) for _ in 1:n_boots] #Vector{Vector{Float64}}()
    # Σs  = [Matrix{Float64}(undef, q, q) for _ in 1:n_boots] #Vector{Matrix{Float64}}()

    # multi-dimensional array approach
    βs  = Matrix{Float64}(undef, n_boots, p)
    Σs  = Array{Float64}(undef, q, q, n_boots) # 3d array
    σ²s = Vector{Float64}(undef, n_boots)
    SubsetEstimates(n_boots, p, q, βs, Σs, σ²s)
end

"""
    blbEstimates

BLB linear mixed model estimates, which contains blb parameters and a vector of `SubsetEstimates`.
"""
struct blbEstimates{T <: LinearAlgebra.BlasReal}
    # blb parameters
    n_subsets::Int64
    subset_size::Int64
    n_boots::Int64
    fenames::Vector{String}
    renames::Vector{String}
    # subset estimates from all subsets
    all_estimates::Vector{SubsetEstimates{T}}
    runtime::Vector{Float64}
    # estimation method
    method::Symbol
end

"""
    save_bootstrap_result!(subset_estimates, β, Σ, σ²)

Save the result from one bootstrap iteration to subset_estimates

## Positional arguments
* `subset_estimates`: an object of the SubsetEstimates type
* `i`: update the ith element of βs, Σs and σ²s
* `β`: parameter estimates for fixed effects
* `Σ`: the covariance matrix of variance components
* `σ²`: estimate of error variance
"""
function save_bootstrap_result!(
    subset_estimates::SubsetEstimates{T},
    i::Int64,
    β::Vector{T},
    Σ::Matrix{T},
    σ²::T
    ) where T <: BlasReal
    # using the vector of vectors/matrices approach
    # copyto!(subset_estimates.βs[i], β)
    # copyto!(subset_estimates.Σs[i], Σ)

    # using the multi-dimensional array approach
    subset_estimates.βs[i, :] = β
    subset_estimates.Σs[:, :, i] = Σ
    subset_estimates.σ²s[i] = σ²
end

"""
    update_w!(m, w)

Update the weight vector `m.w` using `w`.
"""
function update_w!(
    m::Union{blblmmModel{T}, WSVarLmmModel{T}},
    w::Vector{Int64}
    ) where T <: BlasReal
    if typeof(m) == blblmmModel{T}
        copyto!(m.w, w)
    else
        copyto!(m.obswts, w)
    end
end

"""
    blb_one_subset(rng, m; N, subset_size, n_boots, method, solver, verbose, nonparametric_boot)

Performs Bag of Little Bootstraps on a subset.

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `m`: an object of type blblmmModel or WSVarLmmModel

# Keyword arguments
- `N`: number of individuals/clusters in the full dataset.
- `subset_size`: number of individuals/clusters in each subset.
- `n_boots`: number of bootstrap iterations. Default to 1000.
- `method`: a symbol, either `:ML` or `:WiSER`.
- `solver`: solver for the optimization problem.
- `verbose`: Bool, whether to print bootstrap progress.
- `nonparametric_boot`: Bool, whether to use nonparametric bootstrap. For WiSER models, only nonparametric bootstrap is supported.

# Values
- `subset_estimates`: an object of type `SubsetEstimates`
"""
function blb_one_subset(
    rng::Random.AbstractRNG,
    m::Union{blblmmModel{T}, WSVarLmmModel{T}};
    N::Int64,
    subset_size::Int64,
    n_boots::Int = 1000,
    method::Symbol,
    solver,
    verbose::Bool,
    nonparametric_boot::Bool
    ) where T <: BlasReal

    if method == :ML
        MixedModelsBLB.init_ls!(m)
        MixedModelsBLB.fit!(m, solver)
    else
        WiSER.init_ls!(m)
        WiSER.fit!(m, solver, verbose = verbose)
    end

    # Initalize an instance of SubsetEstimates type for storing results
    subset_estimates = SubsetEstimates(n_boots, m.p, m.q)

    # construct the simulator
    if nonparametric_boot
        verbose && print("Using Non-parametric Bootstrap\n")
        simulator = NonparametricBootSimulator(m, N = N, subset_size = subset_size)
    else
        verbose && print("Using Parametric Bootstrap\n")
        simulator = ParametricBootSimulator(m, N = N, subset_size = subset_size)
    end

    # Bootstrapping
    @inbounds for k = 1:n_boots
        if verbose
            flush(stdout)
            print("Bootstrap iteration ", k, "\n")
        end

        if !nonparametric_boot
            # Parametric bootstrapping. Updates m.data[i].y for all i
            simulate!(rng, m, simulator)
        end

        # Get weights by drawing N i.i.d. samples from multinomial
        Random.rand!(rng, simulator.mult_dist, simulator.ns)

        # Update weights in blblmmModel
        update_w!(m, simulator.ns)

        # Fit model on the bootstrap sample
        if method == :ML
            MixedModelsBLB.fit!(m, solver)
        else
            WiSER.fit!(m, solver, verbose = verbose)
        end

        # Save estimates
        if method == :ML
            save_bootstrap_result!(subset_estimates, k, m.β, m.Σ, m.σ²[1])
        else
            save_bootstrap_result!(subset_estimates, k, m.β, m.Σγ, exp(m.τ[1]))
        end

        if !nonparametric_boot
            # Reset model parameter to subset estimates because
            # using the bootstrap estimates from each iteration may be unstable.
            copyto!(m.β, simulator.β_subset)
            copyto!(m.Σ, simulator.Σ_subset)
            copyto!(m.σ², simulator.σ²_subset)
        end
    end
    return subset_estimates
end
blb_one_subset(m::Union{blblmmModel, WSVarLmmModel}; N::Int64, subset_size::Int64, n_boots::Int = 1000,
                method::Symbol, solver, verbose::Bool, nonparametric_boot::Bool)  =
    blb_one_subset(Random.GLOBAL_RNG, m; N = N, subset_size = subset_size, n_boots = n_boots, method = method,
                    solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)

"""
    blblmmobs(datatable)

Constructor for the type `blblmmObs`

# Positional arguments
- `data_obs`: a table object that is compatible with Tables.jl
- `feformula`: the formula for the fixed effects
- `reformula`: the formula for the random effects
"""
function blblmmobs(data_obs, feformula::FormulaTerm, reformula::FormulaTerm)
    y, X = StatsModels.modelcols(feformula, data_obs)
    Z = StatsModels.modelmatrix(reformula, data_obs)
    return blblmmObs(y, X, Z)
end
blblmmobs(feformula::FormulaTerm, reformula::FormulaTerm) = data_obs -> blblmmobs(data_obs, feformula, reformula)

"""
    wsvarlmmobs(datatable)

Constructor for the type `WSVarLmmObs`

# Positional arguments
- `data_obs`: a table object that is compatible with `Tables.jl`
- `feformula`: formula for the fixed effects
- `reformula`: formula for the random effects
- `wsvarformula`: formula for the fixed effects of the within-subject variance
"""
function wsvarlmmobs(data_obs, feformula::FormulaTerm, reformula::FormulaTerm, wsvarformula::FormulaTerm)
    y, X = StatsModels.modelcols(feformula, data_obs)
    Z = StatsModels.modelmatrix(reformula, data_obs)
    W = StatsModels.modelmatrix(wsvarformula, data_obs)
    return WSVarLmmObs(y, X, Z, W)
end
wsvarlmmobs(feformula::FormulaTerm, reformula::FormulaTerm, wsvarformula::FormulaTerm) =
    data_obs -> wsvarlmmobs(data_obs, feformula, reformula, wsvarformula)

"""
    count_levels(data_columns, cat_names)

Count the number of levels of each categorical variable.

# Positional arguments
- `data_columns`: an object of the AbstractColumns type
- `cat_names`: a character vector of the names of the categorical variables

# Values
- `cat_levels`: a dictionary that contains the number of levels of each categorical variable.
"""
function count_levels(data_columns, cat_names::Vector{String})
    cat_levels = Dict{String, Int64}()
    @inbounds for cat_name in cat_names
        cat_levels[cat_name] = length(countmap(Tables.getcolumn(data_columns, Symbol(cat_name))))
    end
    return cat_levels
end
count_levels(cat_names::Vector{String}) = data_columns -> count_levels(data_columns, cat_names)

"""
    subsetting!(rng, subset_id, data_columns, id_name, unique_id, cat_names, cat_levels)

Draw a subset from the full dataset. The IDs of the subset is stored in `subset_id`

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `subset_id`: a vector for storing the IDs of the subset
- `data_columns`: an object of the AbstractColumns type, or a DataFrame
- `unique_id`: a vector of the unique ID in the full data set
- `cat_names`: a character vector of the names of the categorical variables
- `cat_levels`: a dictionary that contains the number of levels of each categorical variable
"""
function subsetting!(
    rng::Random.AbstractRNG,
    subset_id::Vector,
    data_columns,
    id_name::Symbol,
    unique_id::Vector,
    cat_names::Vector{String},
    cat_levels::Dict
    )
    if length(cat_names) == 0
        sample!(rng, unique_id, subset_id; replace = false)
    else
        # there are categorical variables
        good_subset = false
        while !good_subset
            # Sample from the full dataset
            sample!(rng, unique_id, subset_id; replace = false)
            # subset_indices = LinearIndices(id)[findall(in(blb_id_unique), id)]
            cat_levels_subset = data_columns |>
                TableOperations.filter(x -> Tables.getcolumn(x, Symbol(id_name)) .∈ Ref(Set(subset_id))) |>
                Tables.columns |>
                count_levels(cat_names)
            # If the subset levels do not match the full dataset levels,
            # skip the current iteration and take another subset
            if cat_levels_subset == cat_levels
                good_subset = true
            end
        end
    end
    sort!(subset_id)
end
subsetting!(subset_id::Vector, data_columns, id_name::Symbol, unique_id::Vector, cat_names::Vector{String}, cat_levels::Dict) =
    subsetting!(Random.GLOBAL_RNG, subset_id, data_columns, id_name, unique_id, cat_names, cat_levels)

"""
    x_in_y(x, sorted_y, k)

Test whether x is in a sorted vector y or not.

# Positional arguments
- `x`: a value of type T
- `y`: a sorted vector of type T
"""
function x_in_y(x::T, sorted_y::Vector{T}, k::Int) where T
    # test whether x is in sorted_y
    index = searchsortedfirst(sorted_y, x)
    return index ≤ k && sorted_y[index] == x
end

"""
    blb_full_data(rng, datatable; feformula, reformula, wsvarformula, id_name, cat_names, subset_size, n_subsets, n_boots, method, solver, verbose,  nonparametric_boot)

Performs Bag of Little Bootstraps on the full dataset

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `datatable`: a data table type that is compatible with Tables.jl.

# Keyword arguments
- `feformula`: model formula for the fixed effects.
- `reformula`: model formula for the random effects.
- `wsvarformula`: model formula for the fixed effects of the within-subject variance. For linear mixed models, it should be `@formula(y ~ 1)`. Only need to be specified when `method = :WiSER`.
- `id_name`: name of the cluster identifier variable. String.
- `cat_names`: a vector of the names of the categorical variables.
- `subset_size`: number of clusters in the subset.
- `n_subsets`: number of subsets.
- `n_boots`: number of bootstrap iterations. Default to 1000
- `solver`: solver for the optimization problem.
- `method`: fitting the model by maximum-likelihood (:ML) or GEE (:WiSER)
- `verbose`: Bool, whether to print bootstrap progress (percentage completion)
- `nonparametric_boot`: Bool, whether to use Nonparametric bootstrap

# Values
- `result`: an object of the blbEstimates type
"""
function blb_full_data(
    rng::Random.AbstractRNG,
    datatable;
    feformula::FormulaTerm,
    reformula::FormulaTerm,
    wsvarformula::Union{FormulaTerm, Nothing} = nothing,
    id_name::String,
    cat_names::Vector{String} = Vector{String}(),
    subset_size::Int,
    n_subsets::Int = 10,
    n_boots::Int = 1000,
    method::Symbol = :ML,
    solver = Ipopt.IpoptSolver(print_level=0, mehrotra_algorithm = "yes", warm_start_init_point = "yes"),
    verbose::Bool = false,
    nonparametric_boot::Bool = true
    )

    method ∉ [:ML, :WiSER] && error("`method` must be :ML or :WiSER.")
    if method == :WiSER && nonparametric_boot == false
        error("When method = :WiSER, use nonparametric bootstrap only.")
    end

    # Create Tables.Columns type for subsequent processing
    datatable_cols = Tables.columns(datatable) # This step does not allocate memory
    # Get the unique ids, which will be used for subsetting
    typeof(id_name) <: String && (id_name = Symbol(id_name))
    unique_id = unique(Tables.getcolumn(datatable_cols, id_name))
    # Same performance as: unique(Tables.columntable(TableOperations.select(datatable, id_name))[id_name])
    N = length(unique_id) # number of individuals/clusters in the full dataset
    if N < subset_size
        error(string("The subset size should not be bigger than the total number of clusters. \n",
                        "Total number of clusters = ", N, "\n",
                        "Subset size = ", subset_size, "\n"))
    end

    # apply df-wide schema
    feformula = apply_schema(feformula, schema(feformula, datatable))
    reformula = apply_schema(reformula, schema(reformula, datatable))
    if method == :WiSER
        wsvarformula = apply_schema(wsvarformula, schema(wsvarformula, datatable))
    end

    #fenames = coefnames(feformula)[2]
    if typeof(StatsModels.coefnames(feformula)[2]) == String
        fenames = [StatsModels.coefnames(feformula)[2]]
    else
        fenames = StatsModels.coefnames(feformula)[2]
    end

    if typeof(StatsModels.coefnames(reformula)[2]) == String
        renames = [StatsModels.coefnames(reformula)[2]]
    else
        renames = StatsModels.coefnames(reformula)[2]
    end

    if method == :WiSER
        if typeof(StatsModels.coefnames(wsvarformula)[2]) == String
            wsvarnames = [StatsModels.coefnames(wsvarformula)[2]]
        else
            wsvarnames = StatsModels.coefnames(wsvarformula)[2]
        end
    end

    # By chance, some factors of a categorical variable may not show up in a subset.
    # To make sure this does not happen, we
    # 1. count the number of levels of each categorical variable in the full data;
    # 2. for each sampled subset, we check whether the number of levels match.
    # If they do not match, the subset is resampled.
    if length(cat_names) > 0
        cat_levels = count_levels(datatable_cols, cat_names)
    else
        cat_levels = Dict()
    end

    # Initialize an array to store the unique IDs for the subset
    subset_id = Vector{eltype(unique_id)}(undef, subset_size)
    # Initialize a vector of SubsetEstimates for storing estimates from subsets
    all_estimates = Vector{SubsetEstimates{Float64}}(undef, n_subsets)
    # Initialize a vector of obsvec
    if method == :ML
        obsvec = Vector{blblmmObs{Float64}}(undef, subset_size)
    else
        obsvec = Vector{WSVarLmmObs{Float64}}(undef, subset_size)
        clswts = ones(subset_size)
        respname = StatsModels.coefnames(feformula)[1]
    end
    # Initialize a vector for storing the runtime
    # Currently, runtime doesn't make sense when there are multiple workers
    runtime = Vector{Float64}(undef, n_subsets)
    # do grouping once, and collect
    datatable_grouped = datatable |> @groupby(_.id) |> collect

    if Distributed.nworkers() > 1
        # Using multi-processing
        futures = Vector{Future}(undef, n_subsets)
        # wks_schedule assigns each subset to a worker
        wks_schedule = Vector{Int}(undef, n_subsets)
        # nwks = Distributed.nworkers()
        wk = 2 # note that the worker number starts from 2
        wk_max = maximum(workers())
        @inbounds for i = 1:n_subsets
            wks_schedule[i] = wk
            wk == wk_max ? wk = 2 : wk += 1
        end
        flush(stdout)
        print("wks_schedule = ", wks_schedule, "\n")

        @inbounds for j = 1:n_subsets
            time0 = time_ns()
            # Take a subset
            subsetting!(rng, subset_id, datatable_cols, id_name, unique_id, cat_names, cat_levels)
            # Construct blblmmObs objects
            # if use_groupby
            if method == :ML
                obsvec = datatable_grouped |> @filter(x_in_y(key(_), subset_id, subset_size)) |>
                        @map(blblmmobs(_, feformula, reformula)) |> collect |> Array{blblmmObs{Float64}, 1}
                # Construct the blblmmModel type
                m = blblmmModel(obsvec, fenames, renames, N)
            else
                obsvec = datatable_grouped |> @filter(x_in_y(key(_), subset_id, subset_size)) |>
                        @map(wsvarlmmobs(_, feformula, reformula, wsvarformula)) |> collect |> Array{WSVarLmmObs{Float64}, 1}
                m = WSVarLmmModel(obsvec, obswts = clswts, respname = respname, meannames = fenames,
                                    renames = renames, wsvarnames = wsvarnames,
                                    meanformula = feformula, reformula = reformula, wsvarformula = wsvarformula,
                                    ids = subset_id)
            end

            # Process this subset on worker "wks_schedule[j]"
            futures[j] = remotecall(blb_one_subset, wks_schedule[j], rng, m;
                                    N = N, subset_size = subset_size, n_boots = n_boots, method = method,
                                    solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)
            # A remote call returns a Future to its result. Remote calls return immediately;
            # the process that made the call proceeds to its next operation while the remote call happens somewhere else.
            # You can wait for a remote call to finish by calling wait on the returned Future,
            # and you can obtain the full value of the result using fetch.
            runtime[j] = (time_ns() - time0) / 1e9
        end
        @inbounds for j = 1:n_subsets
            # Fetch results from workers
            all_estimates[j] = fetch(futures[j])
        end
    else
        # Not using multi-processing
        @inbounds for j = 1:n_subsets
            time0 = time_ns()
            # Take a subset
            subsetting!(rng, subset_id, datatable_cols, id_name, unique_id, cat_names, cat_levels)
            # Construct blblmmObs objects
            if method == :ML
                obsvec = datatable_grouped |> @filter(x_in_y(key(_), subset_id, subset_size)) |>
                        @map(blblmmobs(_, feformula, reformula)) |> collect |> Array{blblmmObs{Float64}, 1}
                # Construct the blblmmModel type
                m = blblmmModel(obsvec, fenames, renames, N)
            else
                obsvec = datatable_grouped |> @filter(x_in_y(key(_), subset_id, subset_size)) |>
                        @map(wsvarlmmobs(_, feformula, reformula, wsvarformula)) |> collect |> Array{WSVarLmmObs{Float64}, 1}
                m = WSVarLmmModel(obsvec, obswts = clswts, respname = respname, meannames = fenames,
                                    renames = renames, wsvarnames = wsvarnames,
                                    meanformula = feformula, reformula = reformula, wsvarformula = wsvarformula,
                                    ids = subset_id)
            end
            all_estimates[j] = blb_one_subset(rng, m; N = N, subset_size = subset_size, n_boots = n_boots, method = method,
                                              solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)
            runtime[j] = (time_ns() - time0) / 1e9
        end
    end
    # Create a blbEstimates instance for storing results from all subsets
    result = blbEstimates{Float64}(n_subsets, subset_size, n_boots, fenames, renames, all_estimates, runtime, method)
    return result
end

blb_full_data(datatable; feformula::FormulaTerm, reformula::FormulaTerm, wsvarformula::Union{FormulaTerm, Nothing} = nothing, id_name::String,
                cat_names::Vector{String} = Vector{String}(), subset_size::Int, n_subsets::Int = 10, n_boots::Int = 200,
                method::Symbol = :ML, solver = Ipopt.IpoptSolver(), verbose::Bool = false, nonparametric_boot::Bool = true) =
    blb_full_data(Random.GLOBAL_RNG, datatable; feformula = feformula, reformula = reformula, wsvarformula = wsvarformula,
                id_name = id_name, cat_names = cat_names, subset_size = subset_size, n_subsets = n_subsets, n_boots = n_boots,
                method = method, solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)

"""
    blb_db(rng, con, table_name; feformula, reformula, id_name, cat_names, subset_size, n_subsets, n_boots, solver, verbose,  nonparametric_boot)

Performs Bag of Little Bootstraps on databases.

# Positional arguments
- `rng`: random number generator. Default to the global rng.
- `con`: an object of type `MySQL.Connection` created by the function `DBInterface.connect`.
- `table_name`: table name for the longitudinal data.

# Keyword arguments
- `feformula`: model formula for the fixed effects.
- `reformula`: model formula for the random effects.
- `id_name`: name of the cluster identifier variable. String.
- `cat_names`: a vector of the names of the categorical variables.
- `subset_size`: number of clusters in the subset.
- `n_subsets`: number of subsets.
- `n_boots`: number of bootstrap iterations. Default to 1000
- `solver`: solver for the optimization problem.
- `verbose`: Bool, whether to print bootstrap progress (percentage completion)
- `nonparametric_boot`: Bool, whether to use Nonparametric bootstrap

# Values
- `result`: an object of the blbEstimates type
"""
function blb_db(
    rng::Random.AbstractRNG,
    con,
    table_name::String;
    feformula::FormulaTerm,
    reformula::FormulaTerm,
    id_name::String,
    cat_names::Vector{String} = Vector{String}(),
    subset_size::Int,
    n_subsets::Int = 10,
    n_boots::Int = 1000,
    solver = Ipopt.IpoptSolver(print_level=0, mehrotra_algorithm = "yes", warm_start_init_point = "yes"),
    verbose::Bool = false,
    # use_threads::Bool = false,
    nonparametric_boot::Bool = true
    )

    # initialize two empty vectors for storage
    fenames = String[]
    renames = String[]

    query = string("SELECT DISTINCT(", id_name, ") FROM ", table_name, ";")
    unique_id = DataFrame(DBInterface.execute(con,  query))[:, id_name]
    N = length(unique_id)
    if N < subset_size
        error(string("The subset size should not be bigger than the total number of clusters. \n",
                        "Total number of clusters = ", N, "\n",
                        "Subset size = ", subset_size, "\n"))
    end

    # Initialize an array to store the unique IDs for the subset
    subset_id = Vector{eltype(unique_id)}(undef, subset_size)
    # Initialize a vector of SubsetEstimates for storing estimates from subsets
    all_estimates = Vector{SubsetEstimates{Float64}}(undef, n_subsets)
    # Initialize a vector of the blblmmObs objects
    obsvec = Vector{blblmmObs{Float64}}(undef, subset_size)
    # Initialize a vector for storing the runtime
    # Currently, runtime doesn't make sense when there are multiple workers
    runtime = Vector{Float64}(undef, n_subsets)

    if Distributed.nworkers() > 1
        # Using multi-processing
        futures = Vector{Future}(undef, n_subsets)
        # wks_schedule assigns each subset to a worker
        wks_schedule = Vector{Int}(undef, n_subsets)
        # nwks = Distributed.nworkers()
        wk = 2 # note that the worker number starts from 2
        wk_max = maximum(workers())
        @inbounds for i = 1:n_subsets
            wks_schedule[i] = wk
            wk == wk_max ? wk = 2 : wk += 1
        end
        flush(stdout)
        print("wks_schedule = ", wks_schedule, "\n")

        @inbounds for j = 1:n_subsets
            time0 = time_ns()

            # Take a subset
            sample!(unique_id, subset_id, replace = false, ordered = true)

            # create a table to store subset id. this makes the WHERE step below easier
            subset_table = string("subset", j)
            query = string("DROP TABLE IF EXISTS ", subset_table, ";")
            DBInterface.execute(con,  query)
            query = string("CREATE TABLE ", subset_table, " (id INT, INDEX(id));")
            DBInterface.execute(con,  query)
            for i in 1:subset_size
                query = string("INSERT INTO ", subset_table, " VALUES (", subset_id[i], ");")
                DBInterface.execute(con, query)
            end
            # DBInterface.execute(con, "SELECT * FROM subset1") |> DataFrame

            # filter data using subset_id
            query = string("SELECT * FROM ", table_name, " WHERE id IN (SELECT id FROM ", subset_table, ");")
            datatable = DBInterface.execute(con,  query) |> DataFrame
            # drop the temp table
            query = string("DROP TABLE IF EXISTS ", subset_table, ";")
            DBInterface.execute(con,  query)

            # apply schema
            feformula = apply_schema(feformula, schema(feformula, datatable))
            reformula = apply_schema(reformula, schema(reformula, datatable))

            if j == 1
                # do this once on the first subset
                if typeof(coefnames(feformula)[2]) == String
                    puch!(fenames, coefnames(feformula)[2])
                else
                    temp = coefnames(feformula)[2]
                    for k in temp
                        push!(fenames, k)
                    end
                end
                if typeof(coefnames(reformula)[2]) == String
                    puch!(renames, coefnames(reformula)[2])
                else
                    temp = coefnames(reformula)[2]
                    for k in temp
                        push!(renames, k)
                    end
                end
            end

            # group by id and construct obsvec
            datatable_grouped = datatable |> @groupby(_.id) |> collect
            obsvec = datatable_grouped|> @map(MixedModelsBLB.blblmmobs(_, feformula, reformula)) |> collect |>
                        Array{blblmmObs{Float64}, 1}

            # Construct the blblmmModel type
            m = blblmmModel(obsvec, fenames, renames, N)

            # Process this subset on worker "wks_schedule[j]"
            futures[j] = remotecall(blb_one_subset, wks_schedule[j], rng, m;
                                    n_boots = n_boots, method = :ML, solver = solver, verbose = verbose,
                                    nonparametric_boot = nonparametric_boot)
            # A remote call returns a Future to its result. Remote calls return immediately;
            # the process that made the call proceeds to its next operation while the remote call happens somewhere else.
            # You can wait for a remote call to finish by calling wait on the returned Future,
            # and you can obtain the full value of the result using fetch.
            runtime[j] = (time_ns() - time0) / 1e9
        end
        @inbounds for j = 1:n_subsets
            # Fetch results from workers
            all_estimates[j] = fetch(futures[j])
        end
    else
        # Not using multi-processing
        @inbounds for j = 1:n_subsets
            time0 = time_ns()

            # Take a subset
            sample!(unique_id, subset_id, replace = false, ordered = true)

            # create a table to store subset id. this makes the WHERE step below easier
            subset_table = string("subset", j)
            query = string("DROP TABLE IF EXISTS ", subset_table, ";")
            DBInterface.execute(con,  query)
            query = string("CREATE TABLE ", subset_table, " (id INT, INDEX(id));")
            DBInterface.execute(con,  query)
            for i in 1:subset_size
                query = string("INSERT INTO ", subset_table, " VALUES (", subset_id[i], ");")
                DBInterface.execute(con, query)
            end

            # filter data using subset_id
            query = string("SELECT * FROM ", table_name, " WHERE id IN (SELECT id FROM ", subset_table, ");")
            datatable = DBInterface.execute(con,  query) |> DataFrame
            # drop the temp table
            query = string("DROP TABLE IF EXISTS ", subset_table, ";")
            DBInterface.execute(con,  query)

            # apply schema
            feformula = apply_schema(feformula, schema(feformula, datatable))
            reformula = apply_schema(reformula, schema(reformula, datatable))

            if j == 1
                # do this once on the first subset
                if typeof(coefnames(feformula)[2]) == String
                    puch!(fenames, coefnames(feformula)[2])
                else
                    temp = coefnames(feformula)[2]
                    for k in temp
                        push!(fenames, k)
                    end
                end
                if typeof(coefnames(reformula)[2]) == String
                    puch!(renames, coefnames(reformula)[2])
                else
                    temp = coefnames(reformula)[2]
                    for k in temp
                        push!(renames, k)
                    end
                end
            end

            # group by id and construct obsvec
            datatable_grouped = datatable |> @groupby(_.id) |> collect
            obsvec = datatable_grouped|> @map(MixedModelsBLB.blblmmobs(_, feformula, reformula)) |> collect |>
                        Array{blblmmObs{Float64}, 1}

            # Construct the blblmmModel type
            m = blblmmModel(obsvec, fenames, renames, N)
            all_estimates[j] = blb_one_subset(rng, m; n_boots = n_boots, method = :ML, solver = solver, verbose = verbose, nonparametric_boot = nonparametric_boot)
            runtime[j] = (time_ns() - time0) / 1e9
        end
    end

    # Create a blbEstimates instance for storing results from all subsets
    result = blbEstimates{Float64}(n_subsets, subset_size, n_boots, fenames, renames, all_estimates, runtime)
    return result
end

"""
    confint(subset_ests, level)

Calculate confidence intervals using estimates from one subset.

# Positional arguments
- `subset_ests`: an object of type `SubsetEstimates`
- `level`: confidence level, usually set to 0.95
"""
function confint(subset_ests::SubsetEstimates, level::Real)
    # ci_βs = Matrix{Float64}(undef, subset_ests.p, 2) # p-by-2 matrix
    # for i in 1:subset_ests.p
    #     ci_βs[i, :] = StatsBase.percentile(view(subset_ests.βs, :, i), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
    # end
    # q = subset_ests.q
    # ci_Σs = Matrix{Float64}(undef, ◺(q), 2)
    # k = 1
    # # For Σ, we get the CI for the diagonals first, then the upper off-diagonals
    # @inbounds for i in 1:q
    #     ci_Σs[k, :] = StatsBase.percentile(view(subset_ests.Σs, i, i, :), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
    #     k += 1
    # end
    # @inbounds for i in 1:q, j in (i+1):q
    #     ci_Σs[k, :] = StatsBase.percentile(view(subset_ests.Σs, i, j, :), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
    #     k += 1
    # end
    # ci_σ²s = reshape(StatsBase.percentile(subset_ests.σ²s, 100 * [(1 - level) / 2, 1 - (1-level) / 2]), 1, 2)
    # return ci_βs, ci_Σs, ci_σ²s
    ci_βs = Matrix{Float64}(undef, subset_ests.p, 2) # p-by-2 matrix
    for i in 1:subset_ests.p
        if all(y -> isnan(y), view(subset_ests.βs, :, i))
            ci_βs[i, :] .= NaN
        else
            ci_βs[i, :] = StatsBase.percentile(filter(y -> !isnan(y), view(subset_ests.βs, :, i)), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
        end
    end
    q = subset_ests.q
    ci_Σs = Matrix{Float64}(undef, ◺(q), 2)
    k = 1
    # For Σ, we get the CI for the diagonals first, then the upper off-diagonals
    @inbounds for i in 1:q
        if all(y -> isnan(y), view(subset_ests.Σs, i, i, :))
            ci_Σs[k, :] .= NaN
        else
            ci_Σs[k, :] = StatsBase.percentile(filter(y -> !isnan(y), view(subset_ests.Σs, i, i, :)), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
        end
        k += 1
    end
    @inbounds for i in 1:q, j in (i+1):q
        if all(y -> isnan(y), view(subset_ests.Σs, i, j, :))
            ci_Σs[k, :] .= NaN
        else
            ci_Σs[k, :] = StatsBase.percentile(filter(y -> !isnan(y), view(subset_ests.Σs, i, j, :)), 100 * [(1 - level) / 2, 1 - (1-level) / 2])
        end
        k += 1
    end
    if all(y -> isnan(y), subset_ests.σ²s)
        ci_σ²s = [NaN, NaN]
    else
        ci_σ²s = reshape(StatsBase.percentile(filter(y -> !isnan(y), subset_ests.σ²s), 100 * [(1 - level) / 2, 1 - (1-level) / 2]), 1, 2)
    end
    return ci_βs, ci_Σs, ci_σ²s
end

"""
    confint(blb_ests, level)

Calculate confidence intervals using estimates from all subsets.

# Positional arguments
- `blb_ests`: an object of type `blbEstimates`
- `level`: confidence level, usually set to 0.95
"""
function confint(blb_ests::blbEstimates, level::Real)
    # initialize arrays for storing the CIs from each subset
    cis_βs = Array{Float64}(undef, blb_ests.all_estimates[1].p, 2, blb_ests.n_subsets)
    cis_Σs = Array{Float64}(undef, ◺(blb_ests.all_estimates[1].q), 2, blb_ests.n_subsets)
    cis_σ²s = Array{Float64}(undef, 1, 2, blb_ests.n_subsets)
    @inbounds for i in 1:blb_ests.n_subsets
        cis_βs[:, :, i], cis_Σs[:, :, i], cis_σ²s[:, :, i] = confint(blb_ests.all_estimates[i], level)
    end

    if any(y -> isnan(y), cis_βs)
        ci_β = fill(NaN, (blb_ests.all_estimates[1].p, 2))
    else
        ci_β  = mean(cis_βs, dims = 3)[:, :, 1]
    end

    if any(y -> isnan(y), cis_Σs)
        ci_Σ = fill(NaN, (◺(blb_ests.all_estimates[1].q), 2))
    else
        ci_Σ  = mean(cis_Σs, dims = 3)[:, :, 1]
    end

    if any(y -> isnan(y), cis_σ²s)
        ci_σ² = fill(NaN, 1, 2)
    else
        ci_σ² = mean(cis_σ²s, dims = 3)[:, :, 1]
    end

    return ci_β, ci_Σ, ci_σ²
end
confint(blb_ests::blbEstimates) = confint(blb_ests, 0.95)

"""
    fixef(blb_ests)

Calculate BLB fixed effect estimates, which are averages of fixed effect estimates from all subsets.

# Positional arguments
- `blb_ests`: an object of type `blbEstimates`
"""
# returns fixed effect estimates
function fixef(blb_ests::blbEstimates)
    means_βs = Matrix{Float64}(undef, blb_ests.n_subsets, blb_ests.all_estimates[1].p) # n-by-p matrix
    for i in 1:blb_ests.n_subsets
        means_βs[i, :] = mean(blb_ests.all_estimates[i].βs, dims = 1)
    end
    mean_β = mean(means_βs, dims = 1)
    return mean_β
end

"""
    vc(blb_ests)

Calculate BLB variance components estimates, which are averages of variance component estimates from all subsets.

# Positional arguments
- `blb_ests`: an object of type `blbEstimates`
"""
# returns variance components estimates
function vc(blb_ests::blbEstimates)
    means_Σs = Array{Float64}(undef, blb_ests.all_estimates[1].q, blb_ests.all_estimates[1].q, blb_ests.n_subsets)
    # print(means_Σs, "\n")
    means_σ²s = Vector{Float64}(undef, blb_ests.n_subsets)
    for i in 1:blb_ests.n_subsets
        # print(mean(blb_ests.all_estimates[i].Σs, dims = 3), "\n")
        means_Σs[:,:,i] = mean(blb_ests.all_estimates[i].Σs, dims = 3)
        # view(means_Σs, :, :, i) .= mean(blb_ests.all_estimates[i].Σs, dims = 3)
        means_σ²s[i] = mean(blb_ests.all_estimates[i].σ²s)
    end
    mean_Σ = mean(means_Σs, dims = 3)[:, :, 1]
    mean_σ² = mean(means_σ²s)
    return mean_Σ, mean_σ²
end

function StatsBase.coeftable(ests::Vector, ci::Matrix, varnames::Vector{String})
    CoefTable(
        hcat(reshape(ests, :, 1), ci[:, 1], ci[:, 2]),
        ["Estimate", "CI Lower", "CI Upper"],
        varnames
        # 4 # pvalcol
    )
end

function vectorize_Σ(Σ::Matrix)
    q  = size(Σ, 1)
    v = Vector{Float64}(undef, ◺(q))
    idx = 1
    # Extract the diagonals first
    @inbounds for i in 1:q
        v[idx] = Σ[i, i]
        idx += 1
    end
    # Extract the upper off-diagonals
    @inbounds for i in 1:q, j in (i+1):q
        v[idx] = Σ[i, j]
        idx += 1
    end
    return v
end

function vc_names(blb_ests::blbEstimates)
    q = length(blb_ests.renames)
    re_names = Vector{String}(undef, ◺(q))
    idx = 1
    # Extract the diagonals' names first
    @inbounds for i in 1:q
        re_names[idx] = blb_ests.renames[i]
        idx += 1
    end
    # Extract the upper off-diagonals' names
    @inbounds for i in 1:q, j in (i+1):q
        re_names[idx] = string(blb_ests.renames[i], " : ", blb_ests.renames[j])
        idx += 1
    end
    return vcat(re_names, "Residual")
end

function Base.show(io::IO, blb_ests::blbEstimates)
    println("Bag of Little Boostrap (BLB) for linear mixed models.")
    println("Method: ", blb_ests.method)
    println("Number of subsets: ", blb_ests.n_subsets)
    println("Number of grouping factors per subset: ", blb_ests.subset_size)
    println("Number of bootstrap samples per subset: ", blb_ests.n_boots)
    println("Confidence interval level: 95%")
    println(io)

    # calculate all CIs and mean estimates
    ci_β, ci_Σ, ci_σ² = confint(blb_ests)
    mean_β = fixef(blb_ests)
    mean_Σ, mean_σ² = vc(blb_ests)

    println("Variance Components parameters")
    show(io, StatsBase.coeftable(vcat(vectorize_Σ(mean_Σ), mean_σ²), vcat(ci_Σ, ci_σ²), vc_names(blb_ests)))
    println(io)
    println(io)

    println("Fixed-effect parameters")
    show(io, StatsBase.coeftable(vec(mean_β), ci_β, blb_ests.fenames))
end


text:
using GeneticAlgorithms, Distributed

module testga

    using GeneticAlgorithms

    mutable struct TestMonster <: Entity
        genes
        fitness

        TestMonster(num) = new(num, nothing)
    end

    # -------

    function create_entity(num)
        TestMonster(num)
    end

    function fitness(entity)
        println("score: ", entity.genes)
        entity.genes
    end

    function isless(lhs::TestMonster, rhs::TestMonster)
        abs(lhs.fitness) > abs(rhs.fitness)
    end

    function group_entities(population)
        if population[1].fitness >= 16
            return
        end

        freeze(population[1])

        for i in 1:length(population)
            produce([1, i])
        end
    end

    function crossover(entities)
        TestMonster(max(entities[1].genes, entities[2].genes))
    end

    function mutate(entity)
        if rand(Float64) > 0.5
            entity.genes += 1
        end
    end

end

# -------

function test_serial()
    runga(testga; initial_pop_size = 8)
end

# -------

function test_parallel(; nprocs_to_add = 2)
    addprocs(nprocs_to_add)
    @everywhere include(Pkg.dir("GeneticAlgorithms", "test", "testga.jl"))

    println("nprocs: $(nprocs())")

    runga(testga; initial_pop_size = 8)
end

ext: jl


text:
"""
Calculate ethane critical properties
T --- temperature in Kelvin
p --- pressure in bar
ro --- density in mol/l
"""
function ethaneCritical()
	sub=newCritical()
	sub.T=305.322
	sub.p=4.8722*10.0	#MPa to bar
	sub.ro=206.18/30.06904	#mol/l #322 kg/m^3
	return sub
end

"""
Calculate ethane properties at triple point
T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function ethaneTriple()
	sub=newTriple()
	sub.T=90.368
	sub.p=1.14*10^-5	#Pa to bar
#	sub.lRo=999.793/18.015268	#kg/m^3 to mol/l
#	sub.vRo=0.00485458/18.015268	#kg/m^3s to mol/l
	return sub
end

"""
calculate water properties at saturation line
Wagner W. The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use // J. Phys. Chem. Ref. Data. 1999. Vol. 31, № 2. P. 387.
10.1063/1.1461829
T --- temperature in Kelvin
p --- pressure in bar
vRo --- vapor density in mol/l
lRo --- liquid density in mol/l
"""
function ethaneSaturate(T)
a=[-6.48647577, 1.47010078, -1.66261122, 3.57898378, -4.79105705]
b=[1.56138026, -0.381552776, 0.0785372040, 0.0370315089]
c=[-1.89879145, -3.65459262, 0.850562745, 0.363965487, -1.50005943, -2.26690389]
sub=newSaturate()
crit=ethaneCritical()
if (T<crit.T)
	v=1.0-T/crit.T
	sub.p=crit.p*exp(crit.T/T*(a[1]*v+a[2]*v^1.5+a[3]*v^2.5+a[4]*v^3.5+a[5]*v^4.0))
	sub.lRo=crit.ro*exp(b[1]*v^0.329+b[2]*v^(4.0/6.0)+b[3]*v^(8.0/6.0)+b[4]*v^(19.0/6.0))
	sub.vRo=crit.ro*exp(crit.T/T*(c[1]*v^(0.346)+c[2]*v^(5.0/6.0)+c[3]*v+c[4]*v^2.0+c[5]*v^3.0+c[6]*v^5))
	sub.T=T
else
	sub.p=-1
	sub.lRo=-1
	sub.vRo=-1
	sub.T=-1
end
return sub
end

"""
calculate ethane properties
A Reference Equation of State for the Thermodynamic Properties
of Ethane for Temperatures from the Melting Line to 675 K
and Pressures up to 900 MPa
D. Bücker  and W. Wagner
10.1063/1.1859286
T --- temperature, Kelvin
p --- pressure, bar
ro --- density, mol/l
z ---
s --- entropy,
s0 --- ideal part of entropy,
ds --- nonideal part derived k_B, dimensionless
g --- Gibbs energy
g0 --- ideal part of Gibbs energy
dg --- nonideal part of Gibbs energy
f ---  Helmholtz free energy
f0 --- ideal part
df --- nonideal part, dimensionless
u --- internal energy
du --- nonideal part
h --- entalpy
dh --- nonideal part
w --- speed of sound,
cp ---
cv ---
"""
function ethaneSinglephase(T,ro)

sub=newSingle()
crit=ethaneCritical()
delta=ro/crit.ro
tau=crit.T/T
R=8.3144598
mm=0.03006904

nOne=5
nTwo=34
nThree=5

c=[0.00, 0.00, 0.00, 0.00, 0.00,
   1.00, 1.00, 1.00, 1.00, 1.00,	#10
   1.00, 1.00, 1.00, 1.00, 1.00,
   2.00, 2.00, 2.00, 2.00, 2.00,	#20
   2.00, 2.00, 2.00, 3.00, 3.00,
   3.00, 3.00, 3.00, 3.00, 3.00,	#30
   3.00, 3.00, 3.00, 3.00, 4.00,
   4.00, 4.00, 4.00, 4.00, 0.00,	#40
   0.00, 0.00, 0.00, 0.00]

d=[1.00, 1.00, 2.00, 2.00, 4.00,
   1.00, 1.00, 2.00, 2.00, 3.00,	#10
   6.00, 6.00, 7.00, 9.00, 10.00,
   2.00, 4.00, 4.00, 5.00, 5.00,	#20
   6.00, 8.00, 9.00, 2.00, 3.00,
   3.00, 3.00, 4.00, 4.00, 5.00,	#30
   5.00, 6.00, 11.00, 14.00, 3.00,
   3.00, 4.00, 8.00, 10.00, 1.00,	#40
   1.00, 3.00, 3.00, 2.00 ]

t=[ 0.25, 1.00, 0.25, 0.75, 0.75,
	2.00, 4.25, 0.75, 2.25, 3.00,	#10
	1.00, 1.25, 2.75, 1.00, 2.00,
	2.50, 5.50, 7.00, 0.50, 5.50,	#20
	2.50, 4.00, 2.00, 10.0, 16.0,
	18.0, 20.0, 14.0, 18.0, 12.0,	#30
	19.0, 7.00, 15.0, 9.00, 26.0,
	28.0, 28.0, 22.0, 13.0, 0.00,	#40
	3.00, 3.00, 0.00, 3.00]

n=[0.83440745735241, -0.14287360607171e1,
0.34430242210927, -0.42096677920265,
0.12094500886549e-1, -0.57976201597341,
-0.33127037870838e-1, -0.11751654894130,
-0.11160957833067, 0.62181592654406e-1,
0.98481795434443e-1, -0.98268582682358e-1,
-0.23977831007049e-3, 0.69885663328821e-3,
0.19665987803305e-4, -0.14586152207928e-1,
0.46354100536781e-1, 0.60764622180645e-2,
-0.26447330147828e-2, -0.42931872689904e-1,
0.29987786517263e-2, 0.52919335175010e-2,
-0.10383897798198e-2, -0.54260348214694e-1,
-0.21959362918493, 0.35362456650354,
-0.12477390173714, 0.18425693591517,
-0.16192256436754, -0.82770876149064e-1,
0.50160758096437e-1, 0.93614326336655e-2,
-0.27839186242864e-3, 0.23560274071481e-4,
0.39238329738527e-2, -0.76488325813618e-3,
-0.49944304440730e-2, 0.18593386407186e-2,
-0.61404353331199e-3,-0.23312179367924e-2,
0.29301047908760e-2, -0.26912472842883e-3,
0.18413834111814e3, -0.10397127984854e2]

al=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#20
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#30
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 15.0,	#40
	15.0, 15.0, 20.0, 20.0]

bt=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#20
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#30
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 150.0,	#40
	150.0, 150.0, 275.0, 400.0]

gm=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#20
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#30
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 1.05,	#40
	1.05, 1.05, 1.22, 1.16]

ep=[0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#10
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#20
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 0.0,	#30
	0.0, 0.0, 0.0, 0.0, 0.0,
	0.0, 0.0, 0.0, 0.0, 1.0,	#40
	1.0, 1.0, 1.0, 1.0]

phidelta=0.0
phir=0.0
phitau=0.0
phidelta2=0.0
phitau2=0.0
dphideltadtau=0.0

i=0
for ni=1:nOne
	i=i+1
	phidelta=phidelta+n[i]*d[i]*delta^(d[i]-1.0)*tau^t[i]
	phir=phir+n[i]*delta^d[i]*tau^t[i]
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1)
	phidelta2=phidelta2+n[i]*d[i]*(d[i]-1.0)*delta^(d[i]-2)*tau^t[i]
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2.0)
	dphideltadtau=dphideltadtau+n[i]*d[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)
end
for ni=1:nTwo
	i=i+1
	phidelta=phidelta+n[i]*exp(-delta^c[i])*(delta^(d[i]-1.0)*tau^t[i]*(d[i]-c[i]*delta^c[i]))
	phir=phir+n[i]*delta^(d[i])*tau^t[i]*exp(-delta^c[i])
	phitau=phitau+n[i]*t[i]*delta^d[i]*tau^(t[i]-1.0)*exp(-delta^c[i]);
	phidelta2=phidelta2+n[i]*exp(-delta^c[i])*(delta^(d[i]-2)*tau^(t[i])*((d[i]-c[i]*delta^c[i])*(d[i]-1.0-c[i]*delta^c[i])-c[i]^2*delta^c[i]))
	phitau2=phitau2+n[i]*t[i]*(t[i]-1.0)*delta^d[i]*tau^(t[i]-2)*exp(-delta^c[i])
	dphideltadtau=dphideltadtau+n[i]*t[i]*delta^(d[i]-1.0)*tau^(t[i]-1.0)*(d[i]-c[i]*delta^c[i])*exp(-delta^c[i])
end
for pn=1:nThree
	i=i+1
	phidelta=phidelta+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta-2.0*al[i]*(delta-ep[i]))
	phir=phir+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)
	phitau=phitau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(t[i]/tau-2*bt[i]*(tau-gm[i]))
	phidelta2=phidelta2+n[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(-2*al[i]*delta^d[i]+4*al[i]^2*delta^d[i]*(delta-ep[i])^2-4*d[i]*al[i]*delta^(d[i]-1)*(delta-ep[i])+d[i]*(d[i]-1)*delta^(d[i]-2))
	phitau2=phitau2+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*((t[i]/tau-2*bt[i]*(tau-gm[i]))^2-t[i]/tau^2-2*bt[i])
	dphideltadtau=dphideltadtau+n[i]*delta^d[i]*tau^t[i]*exp(-al[i]*(delta-ep[i])^2-bt[i]*(tau-gm[i])^2)*(d[i]/delta -2*al[i]*(delta-ep[i]))*(t[i]/tau-2*bt[i]*(tau-gm[i]))
end

sub.T=T
sub.ro=ro
sub.p=(1.0+phidelta*delta)*ro*R*T/100.0	#bar
sub.z=1.0+phidelta*delta	#
#ideal part
n0=[9.212802589,
-4.682248550,
3.003039265,
1.117433359,
3.467773215,
6.941944640,
5.970850948,
0.0]
gamma0=[0.0, 0.0, 0.0, 1.4091052332, 4.0099170712, 6.5967098342, 13.9798102659,1.0]

#entropy
phi0=log(delta)+n0[1]+n0[2]*tau+n0[3]*log(tau)+n0[4]*log(1.0-exp(-gamma0[4]*tau))+n0[5]*log(1.0-exp(-gamma0[5]*tau))+n0[6]*log(1.0-exp(-gamma0[6]*tau))+n0[7]*log(1.0-exp(-gamma0[7]*tau))+n0[8]*log(1.0-exp(-gamma0[8]*tau))
phi0tau=n0[2]+n0[3]/tau+n0[4]*gamma0[4]*(1.0/(1.0-exp(-gamma0[4]*tau))-1.0)+n0[5]*gamma0[5]*(1.0/(1.0-exp(-gamma0[5]*tau))-1.0)+n0[6]*gamma0[6]*(1.0/(1.0-exp(-gamma0[6]*tau))-1.0)+n0[7]*gamma0[7]*(1.0/(1.0-exp(-gamma0[7]*tau))-1.0)+n0[8]*gamma0[8]*(1.0/(1.0-exp(-gamma0[8]*tau))-1.0)
phi0delta=1/delta
phi0delta2=-1/delta^2
phi0tau2=-n0[3]/tau^2-n0[4]*gamma0[4]^2*exp(-gamma0[4]*tau)*(1.0-exp(-gamma0[4]*tau))^(-2)-n0[5]*gamma0[5]^2*exp(-gamma0[5]*tau)*(1.0-exp(-gamma0[5]*tau))^(-2)-n0[6]*gamma0[6]^2*exp(-gamma0[6]*tau)*(1.0-exp(-gamma0[6]*tau))^(-2)-n0[7]*gamma0[7]^2*exp(-gamma0[7]*tau)*(1.0-exp(-gamma0[7]*tau))^(-2)-n0[8]*gamma0[8]^2*exp(-gamma0[8]*tau)*(1.0-exp(-gamma0[8]*tau))^(-2)
phi0taudelta=0

#Entropy
sub.s=(tau*(phi0tau+phitau)-phi0-phir)*R
sub.s0=(tau*(phi0tau)-phi0)*R
sub.ds=(tau*(phitau)-phir)	#devidet to R

#out.entrfrac=(tau*(phi0tau+phitau)-phi0-phir)/(tau*(phi0tau)-phi0)

dtaudtau=-tau/T
dphidt=dtaudtau*dphideltadtau

#sub.dpdt=(delta*cur_temp*dphidt+delta*phidelta+1)*R*T/100;
#phi0
#phir
#phidelta
#phidelta*delta
#(1+phi0+phir+phidelta*delta)
sub.g=(1+phi0+phir+phidelta*delta)*R*T	#R*cur_temp*
sub.g0=(1+phi0)*R*T
sub.dg=(phir+phidelta*delta)	#deriver RT
#
sub.f=(phi0+phir)*R*T
sub.f0=phi0*R*T
sub.df=phir

#Energy
sub.u=tau*(phi0tau+phitau)*R*T
sub.du=tau*(phitau)
#
sub.h=(1+tau*(phi0tau+phitau)+delta*phidelta)*R*T
sub.dh=1+tau*(phitau)+delta*phidelta;
#speed of sound
sub.w=sqrt(R/mm*T*(1+2*delta*phidelta+delta^2*phidelta2-(1+delta*phidelta-delta*tau*dphideltadtau)^2/(tau^2*(phi0tau2+phitau2))))	#0.018 перевод в килограмы
#
sub.cv=-tau^2*(phi0tau2+phitau2)*R
sub.cp=(-tau^2*(phi0tau2+phitau2)+(1+delta*phidelta-delta*tau*dphideltadtau)^2/(1+2*delta*phidelta+delta^2*phidelta2))*R
sub.dpdt=R/100.0*T*(1.0+2*delta*phidelta+delta^2*phidelta2)
#debug
# println("phi0 $(phi0)")
# println("phi0delta $(phi0delta)")
# println("phi0delta2 $(phi0delta2)")
# println("phi0tau $(phi0tau)")
# println("phi0tau2 $(phi0tau2)")
# println("phi0taudelta $(phi0taudelta)")
#
# println("phir $(phir)")
# println("phidelta $(phidelta)")
# println("phidelta2 $(phidelta2)")
# println("phitau $(phitau)")
# println("phitau2 $(phitau2)")
# println("dphideltadtau $(dphideltadtau)")
#enddebug
return sub
end


text:
"""
	grdsample(cmd0::String="", arg1=[], kwargs...)

Reads a grid file and interpolates it to create a new grid file with either: a
different registration; or a new grid-spacing or number of nodes, and perhaps
also a new sub-region

Full option list at [`grdsample`](http://gmt.soest.hawaii.edu/doc/latest/grdsample.html)

Parameters
----------

- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdsample(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdsample.html#g)
- **I** : **inc** : -- Str or Number --

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/grdsample.html#i)
- $(GMT.opt_R)
- **T** : **toggle** : -- Bool --

    Toggle the node registration for the output grid so as to become the opposite of the input grid
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdsample.html#t)
- $(GMT.opt_V)
- $(GMT.opt_f)
- $(GMT.opt_n)
- $(GMT.opt_r)
- $(GMT.opt_x)
"""
function grdsample(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("grdsample", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_n(cmd, d)
	cmd, = parse_r(cmd, d)
	cmd, = parse_x(cmd, d)

    cmd = add_opt(cmd, 'G', d, [:G :outgrid])
    cmd = add_opt(cmd, 'I', d, [:I :inc])
	cmd = add_opt(cmd, 'T', d, [:T :toggle])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end
	return common_grd(d, cmd, got_fname, 1, "grdsample", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdsample(arg1=[], cmd0::String=""; kw...) = grdsample(cmd0, arg1; kw...)


text:
const TIMEOUT = 1000

"""
This function executes a classical rewriting
algorithm on a Julia expression `ex`. Classical rewriting
applies rule in order with a fixed point iteration:

This algorithm heavily relies on
[RuntimeGeneratedFunctions.jl](https://github.com/SciML/RuntimeGeneratedFunctions.jl)
and the [MatchCore](https://github.com/SciML/RuntimeGeneratedFunctions.jl)
pattern matcher.
**NOTE**: this does not involve the use of [`EGraphs.EGraph`](@ref) or
equality saturation ([`EGraphs.saturate!`](@ref)).
When using `rewrite`, be aware of infinite loops:
Since rules are matched in order in every iteration, it
is possible that commonly used symbolic rules such as commutativity
or associativity of operators may cause this algorithm to
have a cycling computation instantly. This algorithm
detects cycling computation by keeping an history of hashes,
and instantly returns when a cycle is detected.

This algorithm is suitable for simple, deterministic symbolic rewrites.
For more advanced use cases, where it is needed to apply multiple
rewrites at the same time, or it is known that rules are causing loops,
please use [`EGraphs.EGraph`](@ref) and
equality saturation ([`EGraphs.saturate!`](@ref)).
"""
function rewrite(ex, theory::Theory;
        __source__=LineNumberNode(0),
        order=:inner,                   # evaluation order
        m::Module=@__MODULE__,
		timeout::Int=TIMEOUT
    )
    ex=cleanast(ex)

	if !(theory isa Function)
		theory = compile_theory(theory, m; __source__=__source__)
	end

    # n = iteration count. useful to protect against ∞ loops
    # let's use a closure :)
    n = 0
    countit = () -> begin
        n += 1
        n >= timeout ? error("max reduction iterations exceeded") : nothing
    end

    step = x -> cleanast(theory(x, m))
    norm_step = x -> normalize_nocycle(step, x; callback=countit)

    # evaluation order: outer = suitable for symbolic maths
    # inner = suitable for semantics
    walk = if order == :inner
        (x, y) -> df_walk!(x,y; skip_call=true)
    elseif order == :outer
        (x, y) -> bf_walk!(x,y; skip_call=true)
    else
        error(`unknown evaluation order $order`)
    end

    normalize_nocycle(x -> walk(norm_step, x), ex; )
end

macro rewrite(ex, theory, order)
	t = gettheory(theory, __module__)
    rewrite(ex, t; order=order, __source__=__source__, m=__module__) |> quot
end
macro rewrite(ex, theory) :(@rewrite $ex $theory outer) end

# escapes the expression instead of returning it.
macro esc_rewrite(ex, theory, order)
	t = gettheory(theory, __module__)
    rewrite(ex, t; order=order, __source__=__source__, m=__module__) |> esc
end
macro esc_rewrite(ex, theory) :(@ret_reduce $ex $theory outer) end

macro rewriter(te, order)
	t = gettheory(te, __module__)
	quote (ex) -> rewrite(ex, $t;
			order=$order, __source__=$__source__, m=$__module__)
	end
end


text:
"""
    prune(x::GaussianMixture, T::Real, U::Real, J_max::Integer)

Prune the posterior Gaussian-Mixture

Arguments:
- `x::GaussianMixture` Posterior state distribution. [Gaussian Mixture]
- `T::Real` Truncation threshold.  Drop distributions with weight less than T
- `U::Real` Merging threshold.  Merge if (μ_1-μ_2)^T * Σ^-1 * (μ_1-μ_2) < U
- `J_max::Integer` Maximum number of features
"""
function prune(x::GaussianMixture, T::Real, U::Real, J_max::Integer)
	l = 0
	J_k = x.N
	n = length(x.μ[1])

	# prune indices based on truncation threshold
	I = collect(1:J_k)
	I = I[x.w.>T]

	# setup output vectors
	w_new = zeros(J_k)
	μ_new = [zeros(n) for i=1:J_k]
	Σ_new = [zeros(n,n) for i=1:J_k]

	# perform pruning
	while length(I) != 0
		l = l+1
		j = I[argmax(x.w[I])]

		# Merge close features
		L = []
		for i in I
			delta = x.μ[i]-x.μ[j]
			if dot(delta, inv(x.Σ[i])*delta) < U
				push!(L,i)
			end
		end

		# determine merged parameters
		w_tilde = sum(x.w[L])
		μ_tilde = 1/w_tilde * sum(x.w[L].*x.μ[L])
		Σ_tilde = zeros(n,n)
		for i in L
			tmp = x.Σ[i] + (μ_tilde-x.μ[i])*(μ_tilde-x.μ[i])'
			Σ_tilde = Σ_tilde + x.w[i]*tmp
		end
		Σ_tilde = 1/w_tilde * Σ_tilde

		# store values
		w_new[l] = w_tilde
		μ_new[l] = μ_tilde
		Σ_new[l] = Σ_tilde

		# prune merged features
		I = setdiff(I,L)
	end

	if l > J_max 	# only keep the J_max with highest weights
		idxs = sort!([1:length(w_new);], by=i->(w_new[i]), rev=true)
		idxs = idxs[1:J_max]
		w_new = w_new[idxs]
		μ_new = μ_new[idxs]
		Σ_new = Σ_new[idxs]
	else		# drop zeros at end
		w_new = w_new[1:l]
		μ_new = μ_new[1:l]
		Σ_new = Σ_new[1:l]
	end

	x_new = GaussianMixture(w_new,μ_new,Σ_new)
	return x_new
end


text:
########################################################################
# Codes for the paper:
# An Anthropocene Map of Genetic Diversity
#
# Andreia Miraldo, Sen Li, Michael K. Borregaard, Alexander Floréz-Rodriguéz, Shyam Gopalakrishnan, Mirneza Risvanovic, Zhiheng Wang, Carsten Rahbek, Katharine A. Marske & David Nogués-Bravo
#
# Submitted to Science, 2016
# Code in this file by Michael K. Borregaard
#########################################################################

#Load the necessary library
using DataFrames, DataFramesMeta, StatsBase

"""
Calculate the GD (Genetic Diversity) value

**Parameters**
* 'dats' : A master matrix data frame generated by the 'create_master_matrices()' function
"""
function calcGD(dats)
    meanpaspec = by(dats, [:cell, :species], x -> mean(x[:num_per_bp]))
    GD = by(meanpaspec, :cell, x -> mean(x[:x1]))
    names!(GD, [:site, :GDval])
    GD, meanpaspec
end

"""
Calculate the weighted (by total # of sequence pairs) GD (Genetic Diversity) value

**Parameters**
* 'dats' : A master matrix data frame generated by the 'create_master_matrices()' function
"""
function calcGD_weighted(dats)
    meanpaspec = by(dats, [:cell, :species], x -> DataFrame(PIhat = mean(x[:num_per_bp]), numpairs = size(x,1)))
    GD = by(meanpaspec, :cell, x -> mean(x[:PIhat], WeightVec(x[:numpairs])))
    names!(GD, [:site, :GDval])
    GD, meanpaspec
end

"""
Calculate the total number of base pairs in a cell

**Parameters**
* 'dat' :A master matrix data frame generated by the 'create_master_matrices()' function
"""
function totalbasepairs(dat)
    seqs = unique!(dat[[:cell, :species, :seq2, :length_seq2]])
    seq1s = by(dat, [:cell, :species], df->df[[ :seq1, :length_seq1]][1,:])
    names!(seq1s, names(seqs))
    seqs = vcat(seqs, seq1s)
    by(seqs, [:cell], df-> sum(dropna(df[:length_seq2])))
end

"""
Calculate the total number of species in a cell

**Parameters**
* 'dat' :A master matrix data frame generated by the 'create_master_matrices()' function
"""
function totalspecies(dat)
    by(dat, :cell) do df
        sps = unique(df[:species])
        size(sps,1)
    end
end

"""
Calculate the richness and basepairs in a cell

**Parameters**
* 'dat' :A master matrix data frame generated by the 'create_master_matrices()' function
"""
function calc_cellvalues(dats)
    ret, full = calcGD(dats)
    ret[:richness] = by(dats, :cell, totalspecies)[:x1]
    ret[:basepairs] = by(dats, :cell, totalbasepairs)[:x1]
    ret, full
end

function ready_data(filename::AbstractString)
    ret = readtable(filename)
    ret = @where(ret, !isna(:seq1))
    ret = @where(ret, :overlap .> 0.5)
    ret
end

max_stars_repo_path: Miraldo et al. 2016/GD_summary_functions.jl

max_stars_repo_name: mkborregaard/codes-for-published-papers

max_issues_repo_name: DalavanCloud/codes-for-published-papers

max_forks_repo_name: DalavanCloud/codes-for-published-papers


text:
function boundary_condition_L_elas(c)
    g = [0.0, 0.0]
    return g
end
function rhs_L_elas(c)
    f = [0.0, 0.0]
    if c[1] > 0.9
        f[1] = -1.0
    end
    return f
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
function characElt(verts)
    n, n_sides = size(verts, 1), size(verts, 1)
    area_components = verts[:, 1] .* [verts[2:end, 2]..., verts[1, 2]] .-
                     [verts[2:end, 1]..., verts[1, 1]] .* verts[:, 2]
    area = 0.5 * abs(sum(area_components))
    centroid = sum((verts .+ vcat(verts[2:end, :], verts[1, :]')) .*
                    repeat(area_components, 1, 2), dims = 1) / (6 * area)
    diameter = 0 # compute the diameter by looking at every pair of vertices
    for i = 1:(n_sides - 1)
        for j = (i + 1):n_sides
            diameter = max(diameter, norm(verts[i, :] .- verts[j, :]))
        end
    end
    return area, centroid, diameter
end
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-------------------------
"""
   IK, JK, SK, IF, SF = assembKM_vemKsource_elas2D(pv, cells, rhs,
                                                   boundary_condition,
                                                   meshboundary)

Assemble stiffness and source terms of virtual element approach
"""
function assembKM_vemKsource_elas2D(pv, cells, rhs, boundary_condition,
                                    meshboundary)
    # material parameters
    lm = 0.5769
    mu = 0.3846
    # type of input (for overloaded differentiation calls)
    typein = typeof(pv[1])
    n_dofs, n_polys = size(pv, 1), 3 # method has 1 degree of freedom per vertex
    # stiffness coefficients and source coefficients
    IK, JK, IF = Int64[], Int64[], Int64[]
    SK, SF = if typein == Float64
        Float64[], Float64[]
    else
        Any[], Any[]
    end
    # a utility function for wrapping around a vector
    mod_wrap(x, a) = mod(x - 1, a) + 1
    # (3 * 3) matrix D of elasticity over constant strain fields
    Di = [ (2.0 * mu + lm)              lm         0.0;
                       lm   (2.0 * mu + lm)        0.0;
                      0.0              0.0    4.0 * mu]
    for el_id = 1:length(cells)
        #println(cells[el_id])
        v = cells[el_id]      # global indices of the vertices of el
        vc = pv[v, :]   # na * 2 array of vertices coordinates
        na = length(v)
        # start computing the geometric information
        vol, bar, diameter = characElt(vc)
        D = vol * Di
        # basis functions for R = [(1,0), (0,1)  , (-(y-bary, x-barx)) ]
        # basis functions for C = [(x-barx, 0), (0,y-bary) ,(y-bary, x-barx)]
        # (2n*3) Matrix WC of components of basis functions over constant
        # strain fields
        WC = zeros(typein, 2 * na, 3)
        # (2n*3) Matrix WR of components of basis functions over rigid-body
        # motions
        WR = zeros(typein, 2 * na, 3)
        # (2n*3) matrix NC: (NC)_{2i - 1,l} = (c_k(x_i))_0, (NC)_{2i,l} =
        # (c_k(x_i))_1
        NC = zeros(typein, 2 * na, 3)
        # (2n*3) matrix NR: (NR)_{2i - 1,l} = (r_k(x_i))_0, (NR)_{2i,l} =
        # (r_k(x_i))_1
        NR = zeros(typein, 2 * na, 3)
        for i = 1:na
            # coordinates of vertices i-1, i, i+1
            ci = vc[i, :] # this vertex and its neighbours
            cim = vc[mod_wrap(i - 1, na), :]
            cip = vc[mod_wrap(i + 1, na), :]
            # normal vector at vertex (multiplied by its length)
            nor = [cip[2] - cim[2], - (cip[1] - cim[1])]

            q1 = 0.5 * nor[1] / vol
            q2 = 0.5 * nor[2] / vol

            WC[2 * i - 1, 1]   = q1
            WC[2 * i - 1, 2]   = 0.0
            WC[2 * i - 1, 3]   = 0.5 * q2

            WC[2 * i, 1] = 0.0
            WC[2 * i, 2] = q2
            WC[2 * i, 3] = 0.5 * q1

            WR[2 * i - 1, 1]   = 1.0 / na
            WR[2 * i - 1, 2]   = 0.0
            WR[2 * i - 1, 3]   = 0.5 * q2

            WR[2 * i, 1] = 0.0
            WR[2 * i, 2] = 1.0 / na
            WR[2 * i, 3] = -0.5 * q1

            NC[2 * i - 1, 1] = ci[1] - bar[1]
            NC[2 * i - 1, 2] = 0.0
            NC[2 * i - 1, 3] = ci[2] - bar[2]

            NC[2 * i, 1] = 0.0
            NC[2 * i, 2] = ci[2] - bar[2]
            NC[2 * i, 3] = ci[1] - bar[1]

            NR[2 * i - 1, 1] = 1.0
            NR[2 * i - 1, 2] = 0.0
            NR[2 * i - 1, 3] = ci[2] - bar[2]

            NR[2 * i, 1] = 0.0
            NR[2 * i, 2] = 1.0
            NR[2 * i, 3] = - (ci[1] - bar[1])
        end
        # projection matrix over polynomial displacements
        PP = WC * NC' + WR * NR'
        # projection matrix over higher-order functions
        QP = eye(2 * na) - PP
        # polynomial block in the local stiffness matrix
        PROJ = WC * D * WC'
        # higher-order functions block in the local stiffness matrix
        alphaE = tr(D) / tr(NC' * NC)
        STAB = alphaE * QP * QP'
        # local stiffness matrix
        Ke = PROJ + STAB
        # copy local to global
        # K[vert_ids, vert_ids] += local_stiffness
        # F[vert_ids] .+= rhs(centroid) * area / n_sides
        for il = 1:na
            ig = v[il]
            for jl = 1:na
                jg = v[jl]
                push!(IK, 2 * ig - 1)
                push!(IK, 2 * ig - 1)
                push!(IK, 2 * ig)
                push!(IK, 2 * ig)

                push!(JK, 2 * jg - 1)
                push!(JK, 2 * jg)
                push!(JK, 2 * jg - 1)
                push!(JK, 2 * jg)

                push!(SK, Ke[2 * il - 1, 2 * jl - 1])
                push!(SK, Ke[2 * il - 1, 2 * jl])
                push!(SK, Ke[2 * il, 2 * jl - 1])
                push!(SK, Ke[2 * il, 2 * jl])
            end
        end
        # update force vector
        f = rhs(bar)
        for il = 1:na
            ig = v[il]
            push!(IF, 2 * ig - 1)
            push!(IF, 2 * ig)
            push!(SF, vol / na * f[1])
            push!(SF, vol / na * f[2])
        end
    end
    # add boundary contribution in the source term
    nb = length(meshboundary)
    boundary_vals = [boundary_condition(pv[meshboundary[k], :]) for k = 1:nb]
    # boundary dofs
    Vmeshboundary = vcat(2 * meshboundary .- 1, 2 * meshboundary)
    Jinmesh = indexin(JK, Vmeshboundary)
    Jboundary = findall(Jinmesh .!= nothing)
    for l ∈ Jboundary
        push!(IF, IK[l])
        if mod(Jinmesh[l], 2) == 1
            push!(SF, - SK[l] * boundary_vals[(Jinmesh[l] + 1) ÷ 2][1])
        else
            push!(SF, - SK[l] * boundary_vals[Jinmesh[l] ÷ 2][2])
        end
    end
    return IK, JK, SK, IF, SF
end
#-------------------------------------------------------------------------------
"""
    u, p, t, meshboundary = vem_elas2D(filename::String =
                                                         "squarepolmesh_coarse",
                                       nc::Int64 = 1_00;
                                       resolution::Int64 = 400)

N.B.
* nc = 100, 1_000 or 10_000
* adapted code from the article / matlab's code  "The virtual element method in 50 lines of matlab". See https://arxiv.org/pdf/1604.06021.pdf
"""
function vem_elas2D(filename::String = "Lpolmesh", nc::Int64 = 1_00;
                    resolution::Int64 = 400, visible::Bool = false)
    # computes the virtual element solution of the Poisson problem on
    # the specified mesh
    # load the mesh  + pv       : vertices of the cells
    #                + cells    : polygonal cells (⚠ orientation is crucial)
    #                + (pb, tb) : restricted delaunay mesh
    mesh_filename = "$(@__DIR__)/../test/data/$(filename)_$(nc).jld2"
    println("read file $(mesh_filename)")
    JLD2.@load(mesh_filename, ps, Iv, pv, bbelem, elem, pb, tb)
    cells, bcells = bbelem, bbelem
    # boundary points / rhs
    rhs, boundary_condition, meshboundary = if occursin("square", mesh_filename)
        rhs_sqr_elas, boundary_condition_sqr_elas, []
    else
        EPS = 1e-4
        rhs_L_elas, boundary_condition_L_elas, findall(pv[:, 2] .> (1 - EPS))
    end
    # call assemble function
    IK, JK, SK, IF, SF = assembKM_vemKsource_elas2D(pv, cells, rhs,
                                                    boundary_condition,
                                                    meshboundary)
    K, F, internal_dofs, u = solve_elas2D(IK, JK, SK, IF, SF, pv, meshboundary,
                                          boundary_condition)
    # plot
    if resolution > 0
        Iv = visible ? Iv : 1:length(cells)
        plotsolution(u, pv, cells, bcells, Iv, ps,
                     mesh_filename = mesh_filename, resolution = resolution)
    end
    return nothing
end
#-------------------------------------------------------------------------------
function solve_elas2D(IK, JK, SK, IF, SF, pv, meshboundary, boundary_condition)
    n_dofs = 2 * size(pv, 1)
    u = zeros(n_dofs) # degrees of freedom of the virtual element solution
    K = sparse(IK, JK, SK, n_dofs, n_dofs)
    F = Vector(sparsevec(IF, SF, n_dofs))
    nb = length(meshboundary)
    boundary_vals = [boundary_condition(pv[meshboundary[k], :]) for k = 1:nb]
    Vboundary_vals = [boundary_vals[k][l] for l = 1:2, k = 1:nb][:]
    # vertices which aren’t on the boundary
    Vmeshboundary = vcat(2 * meshboundary .- 1, 2 * meshboundary)
    internal_dofs =  setdiff(1:n_dofs, Vmeshboundary)
    # solve
    u[internal_dofs] = K[internal_dofs, internal_dofs] \ F[internal_dofs]
    u[Vmeshboundary] .= Vboundary_vals # set the boundary values
    return K, F, internal_dofs, u
end


text:
"""
David F. Anderson, Masanori Koyama; An asymptotic relationship between coupling
methods for stochastically modeled population processes. IMA J Numer Anal 2015;
35 (4): 1757-1778. doi: 10.1093/imanum/dru044
"""
SplitCoupledJumpProblem(prob::DiffEqBase.AbstractJumpProblem,prob_control::DiffEqBase.AbstractJumpProblem,aggregator::AbstractAggregatorAlgorithm,coupling_map::Vector{Tuple{Int,Int}})= JumpProblem(cat_problems(prob.prob,prob_control.prob),aggregator,build_split_jumps(prob,prob_control,coupling_map)...)

# make new problem by joining initial_data
function cat_problems(prob::DiscreteProblem,prob_control::DiscreteProblem)
  u0_coupled = CoupledArray(prob.u0,prob_control.u0,true)
  DiscreteProblem(u0_coupled,prob.tspan,prob.p)
end

function cat_problems(prob::DiffEqBase.AbstractODEProblem,prob_control::DiffEqBase.AbstractODEProblem)
  l = length(prob.u0) # add l_c = length(prob_control.u0)
  new_f = function (du,u,p,t)
    prob.f(@view(du[1:l]),u.u,p,t)
    prob_control.f(@view(du[l+1:2*l]),u.u_control,p,t)
  end
  u0_coupled = CoupledArray(prob.u0,prob_control.u0,true)
  ODEProblem(new_f,u0_coupled,prob.tspan,prob.p)
end

function cat_problems(prob::DiscreteProblem,prob_control::DiffEqBase.AbstractODEProblem)
  l = length(prob.u0) # add l_c = length(prob_control.u0)
  if !(typeof(prob.f) <: typeof(DiffEqBase.DISCRETE_INPLACE_DEFAULT))
    @warn("Coupling to DiscreteProblem with nontrivial f. Note that, unless scale_by_time=true, the meaning of f will change when using an ODE/SDE/DDE/DAE solver.")
  end
  new_f = function (du,u,p,t)
    prob.f(@view(du[1:l]),u.u,p,t)
    prob_control.f(@view(du[l+1:2*l]),u.u_control,p,t)
  end
  u0_coupled = CoupledArray(prob.u0,prob_control.u0,true)
  ODEProblem(new_f,u0_coupled,prob.tspan,prob.p)
end

function cat_problems(prob::DiffEqBase.AbstractSDEProblem,prob_control::DiffEqBase.AbstractSDEProblem)
  l = length(prob.u0)
  new_f = function (du,u,p,t)
    prob.f(@view(du[1:l]),u.u,p,t)
    prob_control.f(@view(du[l+1:2*l]),u.u_control,p,t)
  end
  new_g = function (du,u,p,t)
    prob.g(@view(du[1:l]),u.u,p,t)
    prob_control.g(@view(du[l+1:2*l]),u.u_control,p,t)
  end
  u0_coupled = CoupledArray(prob.u0,prob_control.u0,true)
  SDEProblem(new_f,new_g,u0_coupled,prob.tspan,prob.p)
end

function cat_problems(prob::DiffEqBase.AbstractSDEProblem,prob_control::DiffEqBase.AbstractODEProblem)
  l = length(prob.u0)
  new_f = function (du,u,p,t)
    prob.f(@view(du[1:l]),u.u,p,t)
    prob_control.f(@view(du[l+1:2*l]),u.u_control,p,t)
  end
  new_g = function (du,u,p,t)
    prob.g(@view(du[1:l]),u.u,p,t)
    for i in l+1:2*l
      du[i] = 0.
    end
  end
  u0_coupled = CoupledArray(prob.u0,prob_control.u0,true)
  SDEProblem(new_f,new_g,u0_coupled,prob.tspan,prob.p)
end

function cat_problems(prob::DiffEqBase.AbstractSDEProblem,prob_control::DiscreteProblem)
  l = length(prob.u0)
  if !(typeof(prob_control.f) <: typeof(DiffEqBase.DISCRETE_INPLACE_DEFAULT))
    @warn("Coupling to DiscreteProblem with nontrivial f. Note that, unless scale_by_time=true, the meaning of f will change when using an ODE/SDE/DDE/DAE solver.")
  end
  new_f = function (du,u,p,t)
    prob.f(@view(du[1:l]),u.u,p,t)
    prob_control.f(@view(du[l+1:2*l]),u.u_control,p,t)
  end
  new_g = function (du,u,p,t)
    prob.g(@view(du[1:l]),u.u,p,t)
    for i in l+1:2*l
      du[i] = 0.
    end
  end
  u0_coupled = CoupledArray(prob.u0,prob_control.u0,true)
  SDEProblem(new_f,new_g,u0_coupled,prob.tspan)
end

cat_problems(prob_control::DiffEqBase.AbstractODEProblem,prob::DiscreteProblem) = cat_problems(prob,prob_control)
cat_problems(prob_control::DiscreteProblem,prob::DiffEqBase.AbstractSDEProblem) = cat_problems(prob,prob_control)
cat_problems(prob_control::DiffEqBase.AbstractODEProblem,prob::DiffEqBase.AbstractSDEProblem) = cat_problems(prob,prob_control)

# this only depends on the jumps in prob, not prob.prob
function build_split_jumps(prob::DiffEqBase.AbstractJumpProblem,prob_control::DiffEqBase.AbstractJumpProblem,coupling_map::Vector{Tuple{Int,Int}})
  num_jumps = length(prob.discrete_jump_aggregation.rates)
  num_jumps_control = length(prob_control.discrete_jump_aggregation.rates)
  jumps = []
  # overallocates, will fix later
  uncoupled = deleteat!(Vector(1:num_jumps),[c[1] for c in coupling_map ])
  uncoupled_control = deleteat!(Vector(1:num_jumps_control),[c[2] for c in coupling_map ])
  for c in uncoupled   # make uncoupled jumps in prob
    new_rate = prob.discrete_jump_aggregation.rates[c]
    new_affect! = prob.discrete_jump_aggregation.affects![c]
    push!(jumps,ConstantRateJump(new_rate,new_affect!))
  end
  for c in uncoupled_control  # make uncoupled jumps in prob_control
    rate = prob_control.discrete_jump_aggregation.rates[c]
    new_rate = (u,p,t)->rate(u.u_control,p,t)
    affect! = prob_control.discrete_jump_aggregation.affects![c]
    new_affect! = function (integrator)
        flip_u!(integrator.u)
        affect!(integrator)
        flip_u!(integrator.u)
    end
    push!(jumps,ConstantRateJump(new_rate,new_affect!))
  end

  for c in coupling_map # make coupled jumps. 3 new jumps for each old one
    rate = prob.discrete_jump_aggregation.rates[c[1]]
    rate_control = prob_control.discrete_jump_aggregation.rates[c[2]]
    affect! = prob.discrete_jump_aggregation.affects![c[1]]
    affect_control! = prob_control.discrete_jump_aggregation.affects![c[2]]
    # shared jump
    new_affect! = function (integrator)
        affect!(integrator)
        flip_u!(integrator.u)
        affect_control!(integrator)
        flip_u!(integrator.u)
    end
    new_rate = (u,p,t)->min(rate(u.u,p,t),rate_control(u.u_control,p,t))
    push!(jumps,ConstantRateJump(new_rate,new_affect!))
     # only prob
    new_affect! = affect!
    new_rate = (u,p,t)->rate(u.u,p,t)-min(rate(u.u,p,t),rate_control(u.u_control,p,t))
    push!(jumps,ConstantRateJump(new_rate,new_affect!))
    # only prob_control
    new_affect! = function (integrator)
        flip_u!(integrator.u)
        affect!(integrator)
        flip_u!(integrator.u)
    end
    new_rate = (u,p,t)->rate_control(u.u_control,p,t)-min(rate(u.u,p,t),rate_control(u.u_control,p,t))
    push!(jumps,ConstantRateJump(new_rate,new_affect!))
  end
  jumps
end


text:
############################################################################
### Type
############################################################################

struct RAMMatrices
    A_ind
    S_ind
    F_ind
    M_ind
    parameters
    colnames
    constants
    size_F
end

############################################################################
### Constructor
############################################################################

function RAMMatrices(;A, S, F, M = nothing, parameters, colnames)
    A_indices = get_parameter_indices(parameters, A)
    S_indices = get_parameter_indices(parameters, S)
    isnothing(M) ? M_indices = nothing : M_indices = get_parameter_indices(parameters, M)
    F_indices = findall([any(isone.(col)) for col in eachcol(F)])
    constants = get_RAMConstants(A, S, M)
    return RAMMatrices(A_indices, S_indices, F_indices, M_indices, parameters, colnames, constants, size(F))
end

############################################################################
### Constants
############################################################################

struct RAMConstant
    matrix
    index
    value
end

import Base.==

function ==(c1::RAMConstant, c2::RAMConstant)
    res = ( (c1.matrix == c2.matrix) & (c1.index == c2.index) &
            (c1.value == c2.value)
    )
    return res
end

function get_RAMConstants(A, S, M)

    constants = Vector{RAMConstant}()

    for index in CartesianIndices(A)
        if (A[index] isa Number) && !iszero(A[index])
            push!(constants, RAMConstant(:A, index, A[index]))
        end
    end

    for index in CartesianIndices(S)
        if (S[index] isa Number) && !iszero(S[index])
            push!(constants, RAMConstant(:S, index, S[index]))
        end
    end

    if !isnothing(M)
        for index in CartesianIndices(M)
            if (M[index] isa Number) && !iszero(M[index])
                push!(constants, RAMConstant(:M, index, M[index]))
            end
        end
    end

    return constants

end

function set_RAMConstant!(A, S, M, rc::RAMConstant)
    if rc.matrix == :A
        A[rc.index] = rc.value
    elseif rc.matrix == :S
        S[rc.index] = rc.value
        S[rc.index[2], rc.index[1]] = rc.value
    elseif rc.matrix == :M
        M[rd.index] = rc.value
    end
end

function set_RAMConstants!(A, S, M, rc_vec::Vector{RAMConstant})
    for rc in rc_vec set_RAMConstant!(A, S, M, rc) end
end

############################################################################
### get RAMMatrices from parameter table
############################################################################

function RAMMatrices(partable::ParameterTable; par_id = nothing)

    if isnothing(par_id)
        parameters, n_par, par_positions = get_par_npar_identifier(partable)
    else
        parameters, n_par, par_positions =
        par_id[:parameters], par_id[:n_par], par_id[:par_positions]
    end

    n_observed = size(partable.variables[:observed_vars], 1)
    n_latent = size(partable.variables[:latent_vars], 1)
    n_node = n_observed + n_latent

    # F indices
    if length(partable.variables[:sorted_vars]) != 0
        F_ind = findall(x -> x ∈ partable.variables[:observed_vars], partable.variables[:sorted_vars])
    else
        F_ind = 1:n_observed
    end

    # indices of the colnames
    if length(partable.variables[:sorted_vars]) != 0
        positions = Dict(zip(partable.variables[:sorted_vars], collect(1:n_observed+n_latent)))
        colnames = copy(partable.variables[:sorted_vars])
    else
        positions = Dict(zip([partable.variables[:observed_vars]; partable.variables[:latent_vars]], collect(1:n_observed+n_latent)))
        colnames = [partable.variables[:observed_vars]; partable.variables[:latent_vars]]
    end

    # fill Matrices
    # known_labels = Dict{Symbol, Int64}()

    A_ind = Vector{Vector{Int64}}(undef, n_par)
    for i in 1:length(A_ind) A_ind[i] = Vector{Int64}() end
    S_ind = Vector{Vector{Int64}}(undef, n_par); S_ind .= [Vector{Int64}()]
    for i in 1:length(S_ind) S_ind[i] = Vector{Int64}() end

    # is there a meanstructure?
    if any(partable.columns[:from] .== Symbol("1"))
        M_ind = Vector{Vector{Int64}}(undef, n_par)
        for i in 1:length(M_ind) M_ind[i] = Vector{Int64}() end
    else
        M_ind = nothing
    end

    # handel constants
    constants = Vector{RAMConstant}()

    for i in 1:length(partable)

        from, parameter_type, to, free, value_fixed, identifier = partable[i]

        row_ind = positions[to]
        if from != Symbol("1") col_ind = positions[from] end

        if !free
            if (parameter_type == :→) & (from == Symbol("1"))
                push!(constants, RAMConstant(:M, row_ind, value_fixed))
            elseif (parameter_type == :→)
                push!(constants, RAMConstant(:A, CartesianIndex(row_ind, col_ind), value_fixed))
            else
                push!(constants, RAMConstant(:S, CartesianIndex(row_ind, col_ind), value_fixed))
            end
        else
            par_ind = par_positions[identifier]
            if (parameter_type == :→) && (from == Symbol("1"))
                push!(M_ind[par_ind], row_ind)
            elseif parameter_type == :→
                push!(A_ind[par_ind], (row_ind + (col_ind-1)*n_node))
            else
                push!(S_ind[par_ind], row_ind + (col_ind-1)*n_node)
                if row_ind != col_ind
                    push!(S_ind[par_ind], col_ind + (row_ind-1)*n_node)
                end
            end
        end

    end

    return RAMMatrices(A_ind, S_ind, F_ind, M_ind, parameters, colnames, constants, (n_observed, n_node))
end

############################################################################
### get parameter table from RAMMatrices
############################################################################

function ParameterTable(ram_matrices::RAMMatrices)

    partable = ParameterTable(nothing)

    colnames = ram_matrices.colnames
    position_names = Dict{Int64, Symbol}(1:length(colnames) .=> colnames)

    # observed and latent variables
    names_obs = colnames[ram_matrices.F_ind]
    names_lat = colnames[findall(x -> !(x ∈ ram_matrices.F_ind), 1:length(colnames))]

    partable.variables = Dict(
        :sorted_vars => Vector{Symbol}(),
        :observed_vars => names_obs,
        :latent_vars => names_lat
    )

    # constants
    for c in ram_matrices.constants
        push!(partable, get_partable_row(c, position_names))
    end

    # parameters
    for (i, par) in enumerate(ram_matrices.parameters)
        push_partable_rows!(
            partable, position_names,
            par, i,
            ram_matrices.A_ind,
            ram_matrices.S_ind,
            ram_matrices.M_ind,
            ram_matrices.size_F[2])
    end

    return partable
end

############################################################################
### get RAMMatrices from EnsembleParameterTable
############################################################################

function RAMMatrices(partable::EnsembleParameterTable)

    ram_matrices = Dict{Symbol, RAMMatrices}()

    parameters, n_par, par_positions = get_par_npar_identifier(partable)
    par_id = Dict(:parameters => parameters, :n_par => n_par, :par_positions => par_positions)

    for key in keys(partable.tables)
        ram_mat = RAMMatrices(partable.tables[key]; par_id = par_id)
        push!(ram_matrices, key => ram_mat)
    end

    return ram_matrices
end

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, ram_matrices::RAMMatrices)
    print_type_name(io, ram_matrices)
    print_field_types(io, ram_matrices)
end

############################################################################
### Additional Functions
############################################################################

function get_par_npar_identifier(partable::ParameterTable)
    parameters = unique(partable.columns[:identifier])
    filter!(x -> x != :const, parameters)
    n_par = length(parameters)
    par_positions = Dict(parameters .=> 1:n_par)
    return parameters, n_par, par_positions
end

function get_par_npar_identifier(partable::EnsembleParameterTable)

    parameters = Vector{Symbol}()
    for key in keys(partable.tables)
        append!(parameters, partable.tables[key].columns[:identifier])
    end
    parameters = unique(parameters)
    filter!(x -> x != :const, parameters)

    n_par = length(parameters)

    par_positions = Dict(parameters .=> 1:n_par)

    return parameters, n_par, par_positions
end

function get_partable_row(c::RAMConstant, position_names)
    # variable names
    from = position_names[c.index[2]]
    to = position_names[c.index[1]]
    # parameter type
    if c.matrix == :A
        parameter_type = :→
    elseif c.matrix == :S
        parameter_type = :↔
    elseif c.matrix == :M
        parameter_type = :→
    end
    free = false
    value_fixed = c.value
    start = 0.0
    estimate = 0.0
    identifier = :const
    return Dict(
        :from => from,
        :parameter_type => parameter_type,
        :to => to,
        :free => free,
        :value_fixed => value_fixed,
        :start => start,
        :estimate => estimate,
        :identifier => identifier)
end

function cartesian_is_known(index, known_indices)
    known = false
    for k_in in known_indices
        if (index == k_in) | ((index[1] == k_in[2]) & (index[2] == k_in[1]))
            known = true
        end
    end
    return known
end

cartesian_is_known(index, known_indices::Nothing) = false

function get_partable_row(par, position_names, index, matrix, n_nod, known_indices)

    # variable names
    if matrix == :M
        from = Symbol("1")
        to = position_names[index]
    else
        index = linear2cartesian(index, (n_nod, n_nod))

        if (matrix == :S) & (cartesian_is_known(index, known_indices))
            return nothing
        elseif matrix == :S
            push!(known_indices, index)
        end

        from = position_names[index[2]]
        to = position_names[index[1]]
    end

    # parameter type
    if matrix == :A
        parameter_type = :→
    elseif matrix == :S
        parameter_type = :↔
    elseif matrix == :M
        parameter_type = :→
    end

    free = true
    value_fixed = 0.0
    start = 0.0
    estimate = 0.0
    identifier = par

    return Dict(
        :from => from,
        :parameter_type => parameter_type,
        :to => to,
        :free => free,
        :value_fixed => value_fixed,
        :start => start,
        :estimate => estimate,
        :identifier => identifier)
end

function push_partable_rows!(partable, position_names, par, i, A_ind, S_ind, M_ind, n_nod)
    A_ind = A_ind[i]
    S_ind = S_ind[i]
    isnothing(M_ind) || (M_ind = M_ind[i])

    for ind in A_ind
        push!(partable, get_partable_row(par, position_names, ind, :A, n_nod, nothing))
    end

    known_indices = Vector{CartesianIndex}()
    for ind in S_ind
        push!(partable, get_partable_row(par, position_names, ind, :S, n_nod, known_indices))
    end

    if !isnothing(M_ind)
        for ind in M_ind
            push!(partable, get_partable_row(par, position_names, ind, :M, n_nod, nothing))
        end
    end

    return nothing

end

function ==(mat1::RAMMatrices, mat2::RAMMatrices)
    res = ( (mat1.A_ind == mat2.A_ind) & (mat1.S_ind == mat2.S_ind) &
            (mat1.F_ind == mat2.F_ind) & (mat1.M_ind == mat2.M_ind) &
            (mat1.parameters == mat2.parameters) &
            (mat1.colnames == mat2.colnames) & (mat1.size_F == mat2.size_F) &
            (mat1.constants == mat2.constants)
    )
    return res
end

function get_group(d::Dict, group)
    return d[group]
end


text:
"""
	rose(cmd0::String="", arg1=nothing; kwargs...)

Reads (length,azimuth) pairs and plot a windrose diagram (polar histograms).

Full option list at [`psrose`]($(GMTdoc)rose.html)

Parameters
----------

- $(GMT.opt_J)
- **A** | **sector** :: [Type => Str | Number]

	Gives the sector width in degrees for sector and rose diagram.
	($(GMTdoc)rose.html#a)
- $(GMT.opt_B)
- **C** | **color** :: [Type => Str | GMTcpt]

	Give a CPT. The mid x-value for each bar is used to look-up the bar color.
	($(GMTdoc)rose.html#c)
- **E** | **vectors** :: [Type => Str]

	Plot vectors showing the principal directions given in the mode_file file.
	($(GMTdoc)rose.html#e)
- **D** | **shift** :: [Type => Bool]

	Shift sectors so that they are centered on the bin interval (e.g., first sector is centered on 0 degrees).
	($(GMTdoc)rose.html#d)
- **F** | **no_scale** :: [Type => Bool]

	Do not draw the scale length bar [Default plots scale in lower right corner].
	($(GMTdoc)rose.html#f)
- **G** | **fill** :: [Type => Str | Number]

	Selects shade, color or pattern for filling the sectors [Default is no fill].
	($(GMTdoc)rose.html#g)
- **I** | **inquire** :: [Type => Bool]

	Inquire. Computes statistics needed to specify a useful -R. No plot is generated.
	($(GMTdoc)rose.html#i)
- **L** | **labels** :: [Type => Str | Number]

	Specify labels for the 0, 90, 180, and 270 degree marks.
	($(GMTdoc)rose.html#l)
- **M** | **vector_params** :: [Type => Str]

	Used with -C to modify vector parameters.
	($(GMTdoc)rose.html#m)
- $(GMT.opt_P)
- **Q** | **alpha** :: [Type => Str | []]

	Sets the confidence level used to determine if the mean resultant is significant.
	($(GMTdoc)rose.html#q)
- $(GMT.opt_R)
- **S** | **norm** | **normalize** :: [Type => Bool]

	Specifies radius of plotted circle (append a unit from c|i|p).
	($(GMTdoc)rose.html#s)
- **T** | **orientation** :: [Type => Bool]

	Specifies that the input data are orientation data (i.e., have a 180 degree ambiguity)
	instead of true 0-360 degree directions [Default].
	($(GMTdoc)rose.html#t)
- **W** | **pen** :: [Type => Str | Tuple]

	Set pen attributes for sector outline or rose plot. [Default is no outline].
	($(GMTdoc)rose.html#w)
- **Z** | **scale** :: [Type => Str]

	Multiply the data radii by scale.
	($(GMTdoc)rose.html#z)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_w)
- $(GMT.opt_swap_xy)
"""
function rose(cmd0::String="", arg1=nothing; first=true, kwargs...)

    gmt_proggy = (IamModern[1]) ? "rose "  : "psrose "
	length(kwargs) == 0 && return monolitic(gmt_proggy, cmd0, arg1)

	arg2 = nothing		# May be needed if GMTcpt type is sent in via C
	N_args = (arg1 === nothing) ? 0 : 1

	d, K, O = init_module(first, kwargs...)		# Also checks if the user wants ONLY the HELP mode

	# If inquire, no plotting so do it and return
	cmd = add_opt(d, "", 'I', [:I :inquire])
	if (cmd != "")
		cmd = add_opt(d, cmd, 'A', [:A :sector])
		if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end
		return gmt("psrose " * cmd, arg1)
	end

	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c")
	cmd, = parse_common_opts(d, cmd, [:UVXY :c :e :p :t :w :params], first)
	cmd  = parse_these_opts(cmd, d, [[:D :shift], [:F :no_scale], [:L :labels], [:M :vector_params],
	                                 [:Q :alpha], [:S :norm :normalize], [:T :orientation], [:Z :scale]])
	cmd = add_opt(d, cmd, 'A', [:A :sector], (width="", rose="_+r"))

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, opt_R)
	if (isa(arg1, Array{<:GMTdataset,1}))  arg1 = arg1[1].data  end	# WHY I HAVE TO DO THIS?

	cmd = add_opt(d, cmd, 'E', [:E :vectors])
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, CPTaliases, 'C', N_args, arg1, arg2)
	cmd = add_opt_fill(cmd, d, [:G :fill], 'G')
	cmd *= opt_pen(d, 'W', [:W :pen])

	return finish_PS_module(d, gmt_proggy * cmd, "", K, O, true, arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
rose!(cmd0::String="", arg1=nothing; kw...) = rose(cmd0, arg1; first=false, kw...)
rose(arg1,  cmd0::String=""; kw...)  = rose(cmd0, arg1; first=true, kw...)
rose!(arg1, cmd0::String=""; kw...) = rose(cmd0, arg1; first=false, kw...)

const psrose  = rose 			# Alias
const psrose! = rose!			# Alias


text:
# This file was generated, do not modify it. # hide
X = (x1=rand(100), x2=rand(100), x3=rand(100))
y = 2X.x1 - X.x2 + 0.05 * randn(100);

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
const woods = let res_init=zeros(7), jac_init=zeros(7,4), x_init=[-3.0, -1.0, -3.0, -1.0]

    function res(x, r)
        r[1] = 10.0 * (x[2] - x[1]^2 )
        r[2] = 1.0 - x[1]
        r[3] = sqrt(90.0) * (x[4] - x[3]^2 )
        r[4] = 1.0 - x[3]
        r[5] = sqrt(9.9) * (x[2] + x[4] - 2.0)
        t = sqrt(0.2)
        r[6] = t * (x[2] - 1.0)
        r[7] = t * (x[4] - 1.0)
        return r
    end

    function jac(x, jac)
        jac[:] .= 0.0
        jac[1, 1] = -20x[1]
        jac[1, 2] = 10.0
        jac[2, 1] = -1.0
        jac[3, 4] = sqrt(90.0)
        jac[3, 3] = -2x[3] * jac[3, 4]
        jac[4, 3] = -1.0
        jac[5, 2] = sqrt(9.9)
        jac[5, 4] = jac[5, 2]
        jac[6, 2] = sqrt(0.2)
        jac[7, 4] = jac[6, 2]
        return jac
    end

    f(;scale=1, verbose=false, print_steps=false) = testone("woods", res, jac,
                                                   res_init, jac_init, x_init;
                                                   scale=scale, verbose=verbose,
                                                   print_steps=print_steps)
end


text:
using Test
using JTAGManager

using Sockets

# Create a little server to connect to.
@async begin
    storage = zeros(UInt8, 10)
    server = listen(2540)
    while true

        socket = accept(server)
        @async while isopen(socket)
            message = readline(socket)

            # Check if writing. Store written data to "storage"
            if startswith(message, "jtag_write")
                # Get the start address and the data.
                split_message = split(message, " ")

                # Add 1 to start address to correct for base 0 addressing.
                start_address = parse(Int, split_message[2]) + 1
                payload = parse.(UInt8, split_message[3:end])

                # Construct writing range.
                write_range = start_address:(start_address + length(payload) - 1)
                storage[write_range] .= payload

            # Check if a read is being performed. Return the requested
            # data.
            elseif startswith(message, "jtag_read")
                split_message = split(message, " ")

                # Add 1 to start address to correct for base 0 addressing.
                start_address = parse(Int, split_message[2]) + 1
                num_bytes = parse(Int, split_message[3])

                read_range = start_address:(start_address + num_bytes - 1)
                return_message = join(storage[read_range], " ")

                println(socket, return_message)
            end
        end
    end
end

@testset "Testing Writer" begin
    # Wait for a little bit for the above task to get scheduled.
    sleep(1)

    # Open up default JTAG.
    jtag = JTAG()
    @test jtag.isopen == false

    # Test 1
    write(jtag, 0, [1,2,3])
    @test jtag.isopen == true
    sleep(0.5)
    data = read(jtag, 0, 3)
    @test data == [1,2,3]

    # Test 2
    write(jtag, 2, [1,2,3,4,5])
    sleep(0.5)
    data = read(jtag, 0, 5)
    @test data == [1,2,1,2,3]

    # Test 3 - single integer version
    write(jtag, 9, 10)
    sleep(0.5)
    data = read(jtag, 9, 1)
    @test data == [10]

    close(jtag)
    @test jtag.isopen == false
end


text:
struct BisectIndexBottleneckSplitter{F}
    f::F
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::BisectIndexBottleneckSplitter, args...; kwargs...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)
        f = oracle_stripe(SparseHint(), method.f, A, args...; kwargs...)

        #=
            search returns the largest j′ such that f(j, j′) <= c, returns
            max(j, j′_lo) - 1 if no such j′ can be found
        =#
        @inline function search(j, j′_lo, j′_hi, k, c)
            j′_lo = max(j, j′_lo)
            let j′
                while j′_lo <= j′_hi
                    j′ = fld2(j′_lo + j′_hi)
                    if f(j, j′, k) <= c
                        j′_lo = j′ + 1
                    else
                        j′_hi = j′ - 1
                    end
                end
                return j′_hi
            end
        end

        spl_lo = ones(Int, K + 1)
        spl_lo[K + 1] = n + 1

        spl_hi = fill(n + 1, K + 1)
        spl_hi[1] = 1

        spl = undefs(Int, K + 1)
        spl[1] = 1
        spl[K + 1] = n + 1

        c_lo, c_hi = bound_stripe(A, K, args..., f) ./ 1

        for k = 1:K
            j′_hi = spl_hi[k + 1]
            j′_lo = max(spl[k], spl_lo[k + 1])
            while j′_lo <= j′_hi
                j′ = fld2(j′_lo + j′_hi)
                c = f(spl[k], j′, k)
                if c_lo <= c < c_hi
                    chk = true
                    spl[k + 1] = j′
                    for k′ = k + 1 : K - 1
                        spl[k′ + 1] = search(spl[k′], spl_lo[k′ + 1], spl_hi[k′ + 1], k′, c)
                        if spl[k′ + 1] < spl[k′]
                            chk = false
                            spl[k′ + 1 : K] .= spl[k′]
                            break
                        end
                    end
                    if chk && f(spl[K], spl[K + 1], K) <= c
                        c_hi = c
                        j′_hi = j′ - 1
                        spl_hi .= spl
                    else
                        c_lo = c
                        j′_lo = j′ + 1
                        spl_lo .= spl
                    end
                elseif c >= c_hi
                    j′_hi = j′ - 1
                else
                    j′_lo = j′ + 1
                end
            end

            if j′_hi < spl[k]
                break
            end

            spl[k + 1] = j′_hi
        end

        return SplitPartition(K, spl_hi)
    end
end

struct FlipBisectIndexBottleneckSplitter{F}
    f::F
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::FlipBisectIndexBottleneckSplitter, args...; kwargs...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)
        f = oracle_stripe(SparseHint(), method.f, A, args...; kwargs...)

        #=
            search returns the smallest j′ such that f(j, j′) <= c, returns
            j′_hi + 1 if no such j′ can be found
        =#
        @inline function search(j, j′_lo, j′_hi, k, c)
            j′_lo = max(j, j′_lo)
            let j′
                while j′_lo <= j′_hi
                    j′ = fld2(j′_lo + j′_hi)
                    if f(j, j′, k) <= c
                        j′_hi = j′ - 1
                    else
                        j′_lo = j′ + 1
                    end
                end
                return j′_lo
            end
        end

        spl_lo = ones(Int, K + 1)
        spl_lo[K + 1] = n + 1

        spl_hi = fill(n + 1, K + 1)
        spl_hi[1] = 1

        spl = undefs(Int, K + 1)
        spl[1] = 1
        spl[K + 1] = n + 1

        c_lo, c_hi = bound_stripe(A, K, args..., f) ./ 1

        for k = 1:K
            j′_hi = spl_hi[k + 1]
            j′_lo = max(spl[k], spl_lo[k + 1])
            while j′_lo <= j′_hi
                j′ = fld2(j′_lo + j′_hi)
                c = f(spl[k], j′, k)
                if c_lo <= c < c_hi
                    chk = true
                    spl[k + 1] = j′
                    for k′ = k + 1 : K - 1
                        spl[k′ + 1] = search(spl[k′], spl_lo[k′ + 1], spl_hi[k′ + 1], k′, c)
                        if spl[k′ + 1] > n + 1
                            chk = false
                            spl[k′ + 1 : K] .= n + 1
                            break
                        end
                    end
                    if chk && f(spl[K], spl[K + 1], K) <= c
                        c_hi = c
                        j′_lo = j′ + 1
                        spl_lo .= spl
                    else
                        c_lo = c
                        j′_hi = j′ - 1
                        spl_hi .= spl
                    end
                elseif c >= c_hi
                    j′_lo = j′ + 1
                else
                    j′_hi = j′ - 1
                end
            end

            if j′_lo > n + 1
                break
            end

            spl[k + 1] = j′_lo
        end

        return SplitPartition(K, spl_lo)
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ONNX_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ONNX")
JLLWrappers.@generate_main_file("ONNX", UUID("4df3ec0b-c10b-5542-9400-16acb2a6b670"))
end  # module ONNX_jll


text:
module Data

include("types.jl")
export NumericalTable

include("storage/storage.jl")

include("manipulation/manipulation.jl")

include("manager/manager.jl")

end  # module


text:
using Conda
Conda.add("scipy")