text:
"""
    RegressionBasedDID <: DiffinDiffsEstimator

Estimation procedure for regression-based difference-in-differences.

A `StatsSpec` for this procedure accepts the following arguments:

| Key | Type restriction | Default value | Description |
|:---|:---|:---|:---|
| `data` | | | A `Tables.jl`-compatible data table |
| `tr` | `DynamicTreatment{SharpDesign}` | | Treatment specification |
| `pr` | `TrendOrUnspecifiedPR{Unconditional,Exact}` | | Parallel trend assumption |
| `yterm` | `AbstractTerm` | | A term for outcome variable |
| `treatname` | `Symbol` | | Column name for the variable representing treatment time |
| `subset` | `Union{BitVector,Nothing}` | `nothing` | Rows from `data` to be used for estimation |
| `weightname` | `Union{Symbol,Nothing}` | `nothing` | Column name of the sample weight variable |
| `vce` | `Vcov.CovarianceEstimator` | `Vcov.CovarianceEstimator` | Variance-covariance estimator |
| `treatintterms` | `TermSet` | `TermSet()` | Terms interacted with the treatment indicators |
| `xterms` | `TermSet` | `TermSet()` | Terms for covariates and fixed effects |
| `contrasts` | `Union{Dict{Symbol,Any},Nothing}` | `nothing` | Contrast coding to be processed by `StatsModels.jl` |
| `drop_singletons` | `Bool` | `true` | Drop singleton observations for fixed effects |
| `nfethreads` | `Int` | `Threads.nthreads()` | Number of threads to be used for solving fixed effects |
| `fetol` | `Float64` | `1e-8` | Tolerance level for the fixed effect solver |
| `femaxiter` | `Int` | `10000` | Maximum number of iterations allowed for the fixed effect solver |
| `cohortinteracted` | `Bool` | `true` | Interact treatment indicators by treatment time |
| `solvelsweights` | `Bool` | `false` | Solve the cell-level least-square weights with default cell partition |
| `lswtnames` | Iterable of `Symbol`s | `tuple()` | Column names from `treatcells` defining the cell partition used for solving least-square weights |
"""
const RegressionBasedDID = DiffinDiffsEstimator{:RegressionBasedDID,
    Tuple{CheckData, GroupTreatintterms, GroupXterms, GroupContrasts,
    CheckVcov, CheckVars, GroupSample,
    ParseFEterms, GroupFEterms, MakeFEs, CheckFEs, MakeWeights, MakeFESolver,
    MakeYXCols, MakeTreatCols, SolveLeastSquares, EstVcov, SolveLeastSquaresWeights}}

"""
    Reg <: DiffinDiffsEstimator

Alias for [`RegressionBasedDID`](@ref).
"""
const Reg = RegressionBasedDID

function valid_didargs(d::Type{Reg}, ::DynamicTreatment{SharpDesign},
        ::TrendOrUnspecifiedPR{Unconditional,Exact}, args::Dict{Symbol,Any})
    name = get(args, :name, "")::String
    treatintterms = haskey(args, :treatintterms) ? args[:treatintterms] : TermSet()
    xterms = haskey(args, :xterms) ? args[:xterms] : TermSet()
    solvelsweights = haskey(args, :lswtnames) || get(args, :solvelsweights, false)::Bool
    ntargs = (data=args[:data],
        tr=args[:tr]::DynamicTreatment{SharpDesign},
        pr=args[:pr]::TrendOrUnspecifiedPR{Unconditional,Exact},
        yterm=args[:yterm]::AbstractTerm,
        treatname=args[:treatname]::Symbol,
        subset=get(args, :subset, nothing)::Union{BitVector,Nothing},
        weightname=get(args, :weightname, nothing)::Union{Symbol,Nothing},
        vce=get(args, :vce, Vcov.RobustCovariance())::Vcov.CovarianceEstimator,
        treatintterms=treatintterms::TermSet,
        xterms=xterms::TermSet,
        contrasts=get(args, :contrasts, nothing)::Union{Dict{Symbol,Any},Nothing},
        drop_singletons=get(args, :drop_singletons, true)::Bool,
        nfethreads=get(args, :nfethreads, Threads.nthreads())::Int,
        fetol=get(args, :fetol, 1e-8)::Float64,
        femaxiter=get(args, :femaxiter, 10000)::Int,
        cohortinteracted=get(args, :cohortinteracted, true)::Bool,
        solvelsweights=solvelsweights::Bool,
        lswtnames=get(args, :lswtnames, ()))
    return name, d, ntargs
end

"""
    RegressionBasedDIDResult{TR,CohortInteracted,Haslsweights} <: DIDResult{TR}

Estimation results from regression-based difference-in-differences.

# Fields
- `coef::Vector{Float64}`: coefficient estimates.
- `vcov::Matrix{Float64}`: variance-covariance matrix for the estimates.
- `vce::CovarianceEstimator`: variance-covariance estiamtor.
- `tr::TR`: treatment specification.
- `pr::AbstractParallel`: parallel trend assumption.
- `treatweights::Vector{Float64}`: total sample weights from observations for which the corresponding treatment indicator takes one.
- `treatcounts::Vector{Int}`: total number of observations for which the corresponding treatment indicator takes one.
- `esample::BitVector`: indicator for the rows from `data` involved in estimation.
- `nobs::Int`: number of observations involved in estimation.
- `dof_residual::Int`: residual degree of freedom.
- `F::Float64`: F-statistic for overall significance of regression model.
- `p::Float64`: p-value corresponding to the F-statistic.
- `yname::String`: name of the outcome variable.
- `coefnames::Vector{String}`: coefficient names.
- `coefinds::Dict{String, Int}`: a map from `coefnames` to integer indices for retrieving estimates by name.
- `treatcells::VecColumnTable`: a tabular description of cells where a treatment indicator takes one.
- `treatname::Symbol`: column name for the variable representing treatment time.
- `yxterms::Dict{AbstractTerm, AbstractTerm}`: a map from all specified terms to concrete terms.
- `yterm::AbstractTerm`: the specified term for outcome variable.
- `xterms::Vector{AbstractTerm}`: the specified terms for covariates and fixed effects.
- `contrasts::Union{Dict{Symbol, Any}, Nothing}`: contrast coding to be processed by `StatsModels.jl`.
- `weightname::Union{Symbol, Nothing}`: column name of the sample weight variable.
- `fenames::Vector{String}`: names of the fixed effects.
- `nfeiterations::Union{Int, Nothing}`: number of iterations for the fixed effect solver to reach convergence.
- `feconverged::Union{Bool, Nothing}`: whether the fixed effect solver has converged.
- `nfesingledropped::Int`: number of singleton observations for fixed effects that have been dropped.
- `lsweights::Union{TableIndexedMatrix, Nothing}`: cell-level least-square weights.
- `cellymeans::Union{Vector{Float64}, Nothing}`: cell-level averages of the outcome variable.
- `cellweights::Union{Vector{Float64}, Nothing}`: total sample weights for each cell.
- `cellcounts::Union{Vector{Int}, Nothing}`: number of observations for each cell.
"""
struct RegressionBasedDIDResult{TR,CohortInteracted,Haslsweights} <: DIDResult{TR}
    coef::Vector{Float64}
    vcov::Matrix{Float64}
    vce::CovarianceEstimator
    tr::TR
    pr::AbstractParallel
    treatweights::Vector{Float64}
    treatcounts::Vector{Int}
    esample::BitVector
    nobs::Int
    dof_residual::Int
    F::Float64
    p::Float64
    yname::String
    coefnames::Vector{String}
    coefinds::Dict{String, Int}
    treatcells::VecColumnTable
    treatname::Symbol
    yxterms::Dict{AbstractTerm, AbstractTerm}
    yterm::AbstractTerm
    xterms::Vector{AbstractTerm}
    contrasts::Union{Dict{Symbol, Any}, Nothing}
    weightname::Union{Symbol, Nothing}
    fenames::Vector{String}
    nfeiterations::Union{Int, Nothing}
    feconverged::Union{Bool, Nothing}
    nfesingledropped::Int
    lsweights::Union{TableIndexedMatrix{Float64, Matrix{Float64}, VecColumnTable, VecColumnTable}, Nothing}
    cellymeans::Union{Vector{Float64}, Nothing}
    cellweights::Union{Vector{Float64}, Nothing}
    cellcounts::Union{Vector{Int}, Nothing}
end

function result(::Type{Reg}, @nospecialize(nt::NamedTuple))
    yterm = nt.yxterms[nt.yterm]
    yname = coefnames(yterm)
    cnames = _treatnames(nt.treatcells)
    cnames = append!(cnames, coefnames.(nt.xterms))[nt.basiscols]
    coefinds = Dict(cnames .=> 1:length(cnames))
    didresult = RegressionBasedDIDResult{typeof(nt.tr),
        nt.cohortinteracted, nt.lsweights!==nothing}(
        nt.coef, nt.vcov_mat, nt.vce, nt.tr, nt.pr, nt.treatweights, nt.treatcounts,
        nt.esample, sum(nt.esample), nt.dof_resid, nt.F, nt.p,
        yname, cnames, coefinds, nt.treatcells, nt.treatname, nt.yxterms,
        yterm, nt.xterms, nt.contrasts, nt.weightname,
        nt.fenames, nt.nfeiterations, nt.feconverged, nt.nsingle,
        nt.lsweights, nt.cellymeans, nt.cellweights, nt.cellcounts)
    return merge(nt, (result=didresult,))
end

"""
    has_fe(r::RegressionBasedDIDResult)

Test whether any fixed effect is involved in regression.
"""
has_fe(r::RegressionBasedDIDResult) = !isempty(r.fenames)

_top_info(r::RegressionBasedDIDResult) = (
    "Number of obs" => nobs(r),
    "Degrees of freedom" => nobs(r) - dof_residual(r),
    "F-statistic" => TestStat(r.F),
    "p-value" => PValue(r.p)
)

_nunique(t, s::Symbol) = length(unique(getproperty(t, s)))

_excluded_rel_str(tr::DynamicTreatment) =
    isempty(tr.exc) ? "none" : join(string.(tr.exc), " ")

_treat_info(r::RegressionBasedDIDResult{<:DynamicTreatment, true}) = (
    "Number of cohorts" => _nunique(r.treatcells, r.treatname),
    "Interactions within cohorts" => length(columnnames(r.treatcells)) - 2,
    "Relative time periods" => _nunique(r.treatcells, :rel),
    "Excluded periods" => NoQuote(_excluded_rel_str(r.tr))
)

_treat_info(r::RegressionBasedDIDResult{<:DynamicTreatment, false}) = (
    "Relative time periods" => _nunique(r.treatcells, :rel),
    "Excluded periods" => NoQuote(_excluded_rel_str(r.tr))
)

_treat_spec(r::RegressionBasedDIDResult{DynamicTreatment{SharpDesign}, true}) =
    "Cohort-interacted sharp dynamic specification"

_treat_spec(r::RegressionBasedDIDResult{DynamicTreatment{SharpDesign}, false}) =
    "Sharp dynamic specification"

_fe_info(r::RegressionBasedDIDResult) = (
    "Converged" => r.feconverged,
    "Singletons dropped" => r.nfesingledropped
)

show(io::IO, ::RegressionBasedDIDResult) = print(io, "Regression-based DID result")

function show(io::IO, ::MIME"text/plain", r::RegressionBasedDIDResult;
        totalwidth::Int=70, interwidth::Int=4+mod(totalwidth,2))
    halfwidth = div(totalwidth-interwidth, 2)
    top_info = _top_info(r)
    fe_info = has_fe(r) ? _fe_info(r) : ()
    tr_info = _treat_info(r)
    blocks = (top_info, tr_info, fe_info)
    fes = has_fe(r) ? join(string.(r.fenames), " ") : "none"
    fetitle = string("Fixed effects: ", fes)
    blocktitles = ("Summary of results: Regression-based DID",
        _treat_spec(r), fetitle[1:min(totalwidth,length(fetitle))])

    for (ib, b) in enumerate(blocks)
        println(io, repeat('─', totalwidth))
        println(io, blocktitles[ib])
        if !isempty(b)
            # Print line between block title and block content
            println(io, repeat('─', totalwidth))
            for (i, e) in enumerate(b)
                print(io, e[1], ':')
                print(io, lpad(e[2], halfwidth - length(e[1]) - 1))
                print(io, isodd(i) ? repeat(' ', interwidth) : '\n')
            end
        end
    end
    print(io, repeat('─', totalwidth))
end

"""
    AggregatedRegDIDResult{TR,Haslsweights,P<:RegressionBasedDIDResult,I} <: AggregatedDIDResult{TR,P}

Estimation results aggregated from a [`RegressionBasedDIDResult`](@ref).
See also [`agg`](@ref).

# Fields
- `parent::P`: the [`RegressionBasedDIDResult`](@ref) from which the results are generated.
- `inds::I`: indices of the coefficient estimates from `parent` used to generate the results.
- `coef::Vector{Float64}`: coefficient estimates.
- `vcov::Matrix{Float64}`: variance-covariance matrix for the estimates.
- `coefweights::Matrix{Float64}`: coefficient weights used to aggregate the coefficient estimates from `parent`.
- `treatweights::Vector{Float64}`: sum of `treatweights` from `parent` over combined `treatcells`.
- `treatcounts::Vector{Int}`: sum of `treatcounts` from `parent` over combined `treatcells`.
- `coefnames::Vector{String}`: coefficient names.
- `coefinds::Dict{String, Int}`: a map from `coefnames` to integer indices for retrieving estimates by name.
- `treatcells::VecColumnTable`: cells combined from the `treatcells` from `parent`.
- `lsweights::Union{TableIndexedMatrix, Nothing}`: cell-level least-square weights.
- `cellymeans::Union{Vector{Float64}, Nothing}`: cell-level averages of the outcome variable.
- `cellweights::Union{Vector{Float64}, Nothing}`: total sample weights for each cell.
- `cellcounts::Union{Vector{Int}, Nothing}`: number of observations for each cell.
"""
struct AggregatedRegDIDResult{TR,Haslsweights,P<:RegressionBasedDIDResult,I} <: AggregatedDIDResult{TR,P}
    parent::P
    inds::I
    coef::Vector{Float64}
    vcov::Matrix{Float64}
    coefweights::Matrix{Float64}
    treatweights::Vector{Float64}
    treatcounts::Vector{Int}
    coefnames::Vector{String}
    coefinds::Dict{String, Int}
    treatcells::VecColumnTable
    lsweights::Union{TableIndexedMatrix{Float64, Matrix{Float64}, VecColumnTable, VecColumnTable}, Nothing}
    cellymeans::Union{Vector{Float64}, Nothing}
    cellweights::Union{Vector{Float64}, Nothing}
    cellcounts::Union{Vector{Int}, Nothing}
end

"""
    agg(r::RegressionBasedDIDResult{<:DynamicTreatment}, names=nothing; kwargs...)

Aggregate coefficient estimates from `r` by values taken by
the columns from `r.treatcells` indexed by `names`
with weights proportional to `treatweights` within each relative time.

# Keywords
- `bys=nothing`: columnwise transformations over `r.treatcells` before grouping by `names`.
- `subset=nothing`: subset of treatment coefficients used for aggregation.
"""
function agg(r::RegressionBasedDIDResult{<:DynamicTreatment}, names=nothing;
        bys=nothing, subset=nothing)
    inds = subset === nothing ? Colon() : _parse_subset(r, subset, false)
    bycells = view(treatcells(r), inds)
    isempty(bycells) && throw(ArgumentError(
        "No coefficient left for aggregation after taking the subset"))
    _parse_bycells!(getfield(bycells, :columns), bycells, bys)
    names === nothing || (bycells = subcolumns(bycells, names, nomissing=false))

    tcells, rows = cellrows(bycells, findcell(bycells))
    ncell = length(rows)
    pcf = view(treatcoef(r), inds)
    cweights = zeros(length(pcf), ncell)
    ptreatweights = view(r.treatweights, inds)
    ptreatcounts = view(r.treatcounts, inds)
    # Ensure the weights for each relative time always sum up to one
    rels = view(r.treatcells.rel, inds)
    for (i, rs) in enumerate(rows)
        if length(rs) > 1
            relgroups = _groupfind(view(rels, rs))
            for ids in values(relgroups)
                if length(ids) > 1
                    cwts = view(ptreatweights, view(rs, ids))
                    cweights[view(rs, ids), i] .= cwts ./ sum(cwts)
                else
                    cweights[rs[ids[1]], i] = 1.0
                end
            end
        else
            cweights[rs[1], i] = 1.0
        end
    end
    cf = cweights' * pcf
    v = cweights' * view(treatvcov(r), inds, inds) * cweights
    treatweights = [sum(ptreatweights[rows[i]]) for i in 1:ncell]
    treatcounts = [sum(ptreatcounts[rows[i]]) for i in 1:ncell]
    cnames = _treatnames(tcells)
    coefinds = Dict(cnames .=> keys(cnames))
    if r.lsweights === nothing
        lswt = nothing
    else
        lswtmat = view(r.lsweights.m, :, inds) * cweights
        lswt = TableIndexedMatrix(lswtmat, r.lsweights.r, tcells)
    end
    return AggregatedRegDIDResult{typeof(r.tr), lswt!==nothing, typeof(r), typeof(inds)}(
        r, inds, cf, v, cweights, treatweights, treatcounts, cnames, coefinds, tcells,
        lswt, r.cellymeans, r.cellweights, r.cellcounts)
end

vce(r::AggregatedRegDIDResult) = vce(parent(r))
treatment(r::AggregatedRegDIDResult) = treatment(parent(r))
nobs(r::AggregatedRegDIDResult) = nobs(parent(r))
outcomename(r::AggregatedRegDIDResult) = outcomename(parent(r))
weights(r::AggregatedRegDIDResult) = weights(parent(r))
treatnames(r::AggregatedRegDIDResult) = coefnames(r)
dof_residual(r::AggregatedRegDIDResult) = dof_residual(parent(r))

"""
    RegDIDResultOrAgg{TR,Haslsweights}

Union type of [`RegressionBasedDIDResult`](@ref) and [`AggregatedRegDIDResult`](@ref).
"""
const RegDIDResultOrAgg{TR,Haslsweights} =
    Union{RegressionBasedDIDResult{TR,<:Any,Haslsweights},
    AggregatedRegDIDResult{TR,Haslsweights}}

"""
    has_lsweights(r::RegDIDResultOrAgg)

Test whether `r` contains computed least-sqaure weights (`r.lsweights!==nothing`).
"""
has_lsweights(::RegDIDResultOrAgg{TR,H}) where {TR,H} = H

"""
    ContrastResult{T,M,R,C} <: AbstractMatrix{T}

Matrix type that holds least-square weights obtained from one or more
[`RegDIDResultOrAgg`](@ref)s computed over the same set of cells
and cell-level averages.
See also [`contrast`](@ref).

The least-square weights are stored in a `Matrix` that can be retrieved
with property name `:m`,
where the weights for each treatment coefficient
are stored columnwise starting from the second column and
the first column contains the cell-level averages of outcome variable.
The indices for cells can be accessed with property name `:r`;
and indices for identifying the coefficients can be accessed with property name `:c`.
The [`RegDIDResultOrAgg`](@ref)s used to generate the `ContrastResult`
can be accessed by calling `parent`.
"""
struct ContrastResult{T,M,R,C} <: AbstractMatrix{T}
    rs::Vector{RegDIDResultOrAgg}
    m::TableIndexedMatrix{T,M,R,C}
    function ContrastResult(rs::Vector{RegDIDResultOrAgg},
            m::TableIndexedMatrix{T,M,R,C}) where {T,M,R,C}
        cnames = columnnames(m.c)
        cnames[1]==:iresult && cnames[2]==:icoef && cnames[3]==:name || throw(ArgumentError(
            "Table paired with column indices has unaccepted column names"))
        return new{T,M,R,C}(rs, m)
    end
end

_getmat(cr::ContrastResult) = getfield(cr, :m)
Base.size(cr::ContrastResult) = size(_getmat(cr))
Base.getindex(cr::ContrastResult, i) = _getmat(cr)[i]
Base.getindex(cr::ContrastResult, i, j) = _getmat(cr)[i,j]
Base.IndexStyle(::Type{<:ContrastResult{T,M}}) where {T,M} = IndexStyle(M)
Base.getproperty(cr::ContrastResult, n::Symbol) = getproperty(_getmat(cr), n)
Base.parent(cr::ContrastResult) = getfield(cr, :rs)

"""
    contrast(r1::RegDIDResultOrAgg, rs::RegDIDResultOrAgg...; kwargs)

Construct a [`ContrastResult`](@ref) by collecting the computed least-square weights
from each of the [`RegDIDResultOrAgg`](@ref).

# Keywords
- `subset=nothing`: indices for cells to be included (rows in output).
- `coefs=nothing`: indices for coefficients from each result to be included (columns in output).
"""
function contrast(r1::RegDIDResultOrAgg, rs::RegDIDResultOrAgg...;
        subset=nothing, coefs=nothing)
    has_lsweights(r1) && all(r->has_lsweights(r), rs) || throw(ArgumentError(
        "Results must contain computed least-sqaure weights"))
    ri = r1.lsweights.r
    rinds = subset === nothing ? Colon() : _parse_subset(ri, subset)
    # Make a copy to avoid accidentally overwriting cells for DIDResult
    ri = deepcopy(view(ri, rinds))
    dcoefs = coefs === nothing ? NamedTuple() : IdDict{Int,Any}(coefs)
    rs = RegDIDResultOrAgg[r1, rs...]
    m = Vector{Array}(undef, length(rs)+1)
    m[1] = view(r1.cellymeans, rinds)
    iresult = [0]
    icoef = [0]
    names = ["cellymeans"]
    for (i, r) in enumerate(rs)
        view(r.lsweights.r, rinds) == ri || throw(ArgumentError(
            "Cells for least-square weights comparisons must be identical across the inputs"))
        cinds = get(dcoefs, i, nothing)
        cinds = cinds === nothing ? (1:ntreatcoef(r)) : _parse_subset(treatcells(r), cinds)
        m[i+1] = view(r.lsweights.m, rinds, cinds)
        ic = view(1:ntreatcoef(r), cinds)
        push!(iresult, (i for k in 1:length(ic))...)
        append!(icoef, ic)
        append!(names, view(treatnames(r), cinds))
    end
    m = hcat(m...)
    ci = VecColumnTable((iresult=iresult, icoef=icoef, name=names))
    return ContrastResult(rs, TableIndexedMatrix(m, ri, ci))
end

function Base.:(==)(x::ContrastResult, y::ContrastResult)
    # Assume no missing
    x.m == y.m || return false
    x.r == y.r || return false
    x.c == y.c || return false
    return parent(x) == parent(y)
end

function Base.sort!(cr::ContrastResult; @nospecialize(kwargs...))
    p = sortperm(cr.r; kwargs...)
    @inbounds for col in cr.r
        col .= col[p]
    end
    @inbounds cr.m .= cr.m[p,:]
    return cr
end

function Base.view(cr::ContrastResult, subset)
    inds = _parse_subset(cr.r, subset)
    r = view(cr.r, inds)
    m = view(cr.m, inds, :)
    return ContrastResult(parent(cr), TableIndexedMatrix(m, r, cr.c))
end

function _checklengthmatch(v, name::String, N::Int)
    length(v) == N || throw(ArgumentError(
        "The length of $name ($(length(v))) does not match the number of rows of cr ($(N))"))
end

_checklengthmatch(v::Nothing, name::String, N::Int) = nothing

"""
    post!(gl, gr, gd, ::StataPostHDF, cr::ContrastResult, left::Int=2, right::Int=3; kwargs...)

Export the least-square weights for coefficients indexed by `left` and `right`
from `cr` for Stata module [`posthdf`](https://github.com/junyuan-chen/posthdf).
The contribution of each cell to the difference between two coefficients
are computed and also exported.
The weights and contributions are stored as coefficient estimates
in three groups `gl`, `gr` and `gd` respectively.
The groups can be `HDF5.Group`s or objects that can be indexed by strings.

# Keywords
- `lefttag::String=string(left)`: name to be used as `depvar` in Stata after being prefixed by `"l_"` for the coefficient indexed by `left`.
- `righttag::String=string(right)`: name to be used as `depvar` in Stata after being prefixed by `"r_"` for the coefficient indexed by `right`.
- `model::String="InteractionWeightedDIDs.ContrastResult"`: name of the model.
- `eqnames::Union{AbstractVector, Nothing}=nothing`: equation names prefixed to coefficient names in Stata.
- `colnames::Union{AbstractVector, Nothing}=nothing`: column names used as coefficient names in Stata.
- `at::Union{AbstractVector{<:Real}, Nothing}=nothing`: the `at` vector in Stata.
"""
function post!(gl, gr, gd, ::StataPostHDF, cr::ContrastResult, left::Int=2, right::Int=3;
        lefttag::String=string(left), righttag::String=string(right),
        model::String="InteractionWeightedDIDs.ContrastResult",
        eqnames::Union{AbstractVector, Nothing}=nothing,
        colnames::Union{AbstractVector, Nothing}=nothing,
        at::Union{AbstractVector{<:Real}, Nothing}=nothing)
    N = size(cr.m, 1)
    _checklengthmatch(eqnames, "eqnames", N)
    _checklengthmatch(colnames, "colnames", N)
    _checklengthmatch(at, "at", N)
    gl["depvar"] = string("l_", lefttag)
    wtl = view(cr.m, :, left)[:]
    gl["b"] = wtl
    gr["depvar"] = string("r_", righttag)
    wtr = view(cr.m, :, right)[:]
    gr["b"] = wtr
    gd["depvar"] = string("d_", lefttag, "_", righttag)
    diff = (wtl.-wtr).*view(cr.m,:,1)[:]
    gd["b"] = diff
    colnames === nothing && (colnames = 1:N)
    cnames = eqnames === nothing ? string.(colnames) : string.(eqnames, ":", colnames)
    for g in (gl, gr, gd)
        g["model"] = model
        g["coefnames"] = cnames
        at === nothing || (g["at"] = at)
    end
end


text:
#ENV["CPLEX_STUDIO_BINARIES"]="/home/software/cplex/128/cplex/bin/x86-64_linux"
#ENV["MOSEKBINDIR"]="/home/software/mosek/9/tools/platform/linux64x86/bin"
#ENV["MOSEKLM_LICENSE_FILE"]="/home/software/mosek/mosek.lic"
include("src/core.jl")
function getProblemDataFG_minreturn(pathName, pathName2, pathName3, probSize, type, instance, theK)
    # Size: 200, 300, 400
    # Type: Pard
    # Instance: 1-10, converted to equivalent letter.
    alphabet=["a","b","c","d","e","f","g","h","i","j"]

    # Adjust file path appropriately
    if probSize==200
        if type==:Pard
            pathName=pathName*string(probSize)*"/pard"*string(probSize)*"_"*alphabet[instance]
        end
    else
        if type==:Pard
            pathName=pathName*string(probSize)*"/pard"*string(probSize)*"_"*alphabet[instance]
        end
    end

    # Read in .txt file to find n, μ
    returnData=readdlm(pathName*".txt")
    n=returnData[1,1]
    μ=returnData[2:n+1, 1]
    # Read in .rho file to find r_min
    minReturn=readdlm(pathName*".rho")
    r_bar=minReturn[1,1]
    # Read in .bds file to find min_investment, u
    bounds=readdlm(pathName*".bds")
    min_investment=bounds[1:n, 1]
    ub=bounds[1:n, 2]
    # Read in .mat to find sigma
    cov=readdlm(pathName*".mat")
    Σ=cov[2:n+1, 1:n]
    Σ=convert(Array{Float64,2}, Σ)

    # Read in diagonal vector precomputed by Zheng et al
    if probSize==200
        if type==:Pard
            pathName2=pathName2*"/pard"*string(probSize)*"_"*alphabet[instance]
        end
    else
        if type==:Pard
            pathName2=pathName2*"/pard"*string(probSize)*"_"*alphabet[instance]
        end
    end

    if probSize==200
        if type==:Pard
            pathName3=pathName3*"/pard"*string(probSize)*"_"*alphabet[instance]
        end
    else
        if type==:Pard
            pathName3=pathName3*"/pard"*string(probSize)*"_"*alphabet[instance]
        end
    end
    γ=(1000/n)*ones(n)

    # Purpose of these lines was to attempt to take convex combination of diagonal matrices, as proposed in Zheng et. al. (2014). Unfortunately, just using the Frangioni/Gentile is faster
    dominance=readdlm(pathName2*".diag")
    theν=dominance[2:n+1,1]
    dominance=readdlm(pathName3*".diag")
    theν=0.7*theν+0.3*dominance[2:n+1,1]

    for i=1:n
        γ[i]=((1.0/γ[i])+theν[i])^(-1) #adjust γ to have more diagonal dominance, Σ to be smaller
    end
    Σ=Σ-Diagonal(theν)
    X=sqrt(Σ)
    k=theK

    A=vcat(μ', Matrix(1.0I, n,n)) # Can also explicitly create an upper bound on x, if creating a diagonal matrix is too slow/cumbersome.
    l=vcat([r_bar], zeros(n))
    u=vcat([1e12], ub)
    f=size(X, 2)
    m=size(A, 1)
    μ=zeros(n) #Since constraining μ rather than putting in objective
    Y=zeros(f)
    d=zeros(n)

    return SparsePortfolioData(μ, Y, d, X, A, l, u, k, n, m, f, min_investment, γ)
end

for portfolioNumber in [8]

    for theK in [10]
     thePath=pwd()
     pathName = thePath*"/data/MV/size"
     pathName2 = thePath*"/data/diagonals/s/"
     pathName3 = thePath*"/data/diagonals/s/"
     thePortfolioData=getProblemDataFG_minreturn(pathName,pathName2, pathName3, 200, :Pard, portfolioNumber, theK)

     sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=true, minReturnConstraint=true, useHeuristic=false, useWarmStart=false, useCopyOfVariables=true);
     @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
     println("Variance is:", sparse_Portfolio.portfolioVariance)
     println("Return is:", sparse_Portfolio.expectedReturn)
     println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
     open("results/Table6ResultsOutput_noddterm.csv","a") do fp
     println(fp, "Instance, cardinality, regularizer, method, solvetime, nodecount, iters")
     println(fp, "pard400_"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                       "1000/(n)", ",", "OA", ",",  sparse_Portfolio.solveTime, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
     end

     sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=true, minReturnConstraint=true, useHeuristic=false, useWarmStart=false, useCopyOfVariables=true, usingKelleyPrimal=true, maxCutCallbacks=0);
     @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
     println("Variance is:", sparse_Portfolio.portfolioVariance)
     println("Return is:", sparse_Portfolio.expectedReturn)
     println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
     open("results/Table6ResultsOutput_noddterm.csv","a") do fp
     println(fp, "Instance, cardinality, regularizer, method, solvetime, nodecount, iters")
     println(fp, "pard400_"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                       "1000/(n)", ",", "OA Kelley", ",",  sparse_Portfolio.solveTime, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
     end

     sparse_Portfolio=cutting_planes_portfolios(thePortfolioData; useSOCPLB=true, minReturnConstraint=true, useHeuristic=false, useWarmStart=false, useCopyOfVariables=true, usingKelleyPrimal=true, maxCutCallbacks=200);
     @show @test(abs.(sum(sparse_Portfolio.w)-1.0)<1e-4)
     println("Variance is:", sparse_Portfolio.portfolioVariance)
     println("Return is:", sparse_Portfolio.expectedReturn)
     println("Reg term is:", 0.5*sum(sparse_Portfolio.w[i]^2/thePortfolioData.γ[sparse_Portfolio.indices[i]] for i=1:size(sparse_Portfolio.indices,1)))
     open("results/Table6ResultsOutput_noddterm.csv","a") do fp
     println(fp, "Instance, cardinality, regularizer, method, solvetime, nodecount, iters")
     println(fp, "pard400_"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                       "1000/(n)", ",", "OA Kelley 50", ",",  sparse_Portfolio.solveTime, ",", sparse_Portfolio.nodeCount, ",", sparse_Portfolio.numIters)
     end

    sparse_Portfolio_Cplex=cplex_raw_MISOCP(thePortfolioData, 600.0, 1e-8, portfolioNumber)
    @show @test(abs.(sum(sparse_Portfolio_Cplex.w)-1.0)<1e-4)
    println("Variance is:", sparse_Portfolio_Cplex.portfolioVariance)
    println("Return is:", sparse_Portfolio_Cplex.expectedReturn)
    println("Reg term is:", 0.5*sum(sparse_Portfolio_Cplex.w[i]^2/thePortfolioData.γ[sparse_Portfolio_Cplex.indices[i]] for i=1:size(sparse_Portfolio_Cplex.indices,1)))
    open("results/Table6ResultsOutput_noddterm.csv","a") do fp
    println(fp, "Instance, cardinality, regularizer, method, solvetime")
    println(fp, "pard400_"*string(portfolioNumber), ",", thePortfolioData.k, ",",
                        "1000/sqrt(n)", ",", "CPLEX MISOCP", ",",  sparse_Portfolio_Cplex.solveTime, ",", sparse_Portfolio_Cplex.nodeCount)
        end
    end
end

ext: jl


text:
function post(alg::TMJets20{N}, ivp::IVP{<:AbstractContinuousSystem}, timespan;
              Δt0::TimeInterval=zeroI,
              external::Bool=false,    # if `true`, use the external solver defined in TaylorModels.jl
              kwargs...) where {N}

    @unpack orderQ, orderT, abstol, maxsteps, adaptive, minabstol, disjointness = alg

    # initial time and final time
    t0 = tstart(timespan)
    T = tend(timespan)

    # vector field
    if islinear(ivp) || isaffine(ivp)
        f! = inplace_field!(ivp)
    else
        f! = VectorField(ivp)
    end
    parse_eqs = get(kwargs, :parse_eqs, false)

    n = statedim(ivp)
    ivp_norm = _normalize(ivp)
    X = stateset(ivp_norm)

    # fix the working variables and maximum order in the global
    # parameters struct (_params_TaylorN_)
    set_variables("x", numvars=n, order=2*orderQ)

    # initial set
    X0 = initial_state(ivp_norm)

    if external
        return _solve_external(f!, X0, t0, T, orderQ, orderT, abstol, maxsteps, Δt0; parse_eqs=parse_eqs, kwargs...)
    end

    X0tm = _initialize(X0, orderQ, orderT)

    # optionally absorb initial remainder
    shrink_wrapping = get(kwargs, :shrink_wrapping, true)
    if shrink_wrapping && isa(X0tm, TaylorModelReachSet)
        if !all(iszero, remainder(X0tm))
            X0tm = _shrink_wrapping(X0tm)
        end
    end

    # preallocate output flowpipe
    F = Vector{TaylorModelReachSet{N}}()
    sizehint!(F, maxsteps)

    F, tv, xv, xTM1v, success, _t0 = validated_integ!(F, f!, X0tm, t0, T, orderQ, orderT,
                                                      abstol, maxsteps, X, disjointness, Δt0, adaptive; parse_eqs=parse_eqs)

    if success || !adaptive
        ext = Dict{Symbol, Any}(:tv => tv, :xv => xv, :xTM1v => xTM1v, :actual_abs_tol=>[abstol])
        return Flowpipe(F, ext)
    end

    # save extra data, one vector per iteration
    tv_vec = Vector{typeof(tv)}()
    xv_vec = Vector{typeof(xv)}()
    xTM1v_vec = Vector{typeof(xTM1v)}()

    push!(tv_vec, tv)
    push!(xv_vec, xv)
    push!(xTM1v_vec, xTM1v)

    abstol_vec = [abstol]

    while !success
            # adapt the absolute tolerance
            if abstol > minabstol
                abstol = abstol / 10
                push!(abstol_vec, abstol)
            else
                @warn("Minimum absolute tolerance, $minabstol reached.")
                ext = Dict{Symbol, Any}(:tv => tv, :xv => xv, :xTM1v => xTM1v)
                return Flowpipe(F, ext)
            end

            # new initial states
            if !isempty(F)
                # here we pass the box overapproximation of the final reach-set
                # X0_end_box = set(overapproximate(F[end], Hyperrectangle, orderQ=orderQ, orderT=orderT))
                # X0tm = _initialize(X0_end_box, orderQ, orderT)
                X0tm = xv_vec[end][end]
            end

            # new flowpipe
            Fk = Vector{TaylorModelReachSet{N}}()
            sizehint!(Fk, maxsteps)

            Fk, tv, xv, xTM1v, success, _t0 = validated_integ!(Fk, f!, X0tm, _t0, T, orderQ, orderT,
                                                               abstol, maxsteps, X, disjointness, Δt0, adaptive; parse_eqs=parse_eqs)

            # append the new flowpipe to the accumulated flowpipe and extra data
            append!(F, Fk)

            push!(tv_vec, tv)
            push!(xv_vec, xv)
            push!(xTM1v_vec, xTM1v)
    end

    ext = Dict{Symbol, Any}(:tv => tv_vec, :xv => xv_vec, :xTM1v => xTM1v_vec, :actual_abs_tol=>abstol_vec)
    return Flowpipe(F, ext)
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_PT1

Simulate PT1 model.
"""
module Simulate_PT1

include(joinpath("models", "PT1.jl"))

import .PT1
import ModiaMath

model  = PT1.Model(T=2.0, x0=1.5)
result = ModiaMath.simulate!(model, stopTime=10.0; log=true)

ModiaMath.plot(result, ["x[1]", "der(x[1])"], heading="PT1 model (T*der(x) + x = 0; T = " * string(model.T) * ")")

end


text:
using Pkg
pkg"activate ."
using Revise
using Layered
using Colors

function petals()

    c, tl = canvas(4, 4)

    l = rectlayer!(tl, c.rect, :w, margin=30)

    degrees = range(0, 360, length=21)[1:end-1]

    petals = path!.(ang -> begin
        endpoint = P(ang)
        Path(true, Arc(O, endpoint, 0.2), Arc(endpoint, O, 0.2))
    end, l, deg.(degrees)) .+ Stroke.(LCHuv.(70, 50, degrees)) .+ Linewidth(10) .+
        Fill.(LCHuvA.(70, 50, degrees, 0.5)) .+ Operator(:mult)

    circ = circle!(l, O, 0.25) + Visible(false)

    l + Clip(circ, c.rect)

    draw_svg(c, "./examples/03_petals/petals.svg")
    # c

end; petals()


text:
export degrees, degree, mindegree_vertex, maxdegree_vertex, minmaxdegree_vertex

"""Get degrees of all vertices in a graph."""
degrees(eg::EliminateGraph) = degree.(Ref(eg), vertices(eg))

"""
    neighbors(eg::EliminateGraph, i::Int)

Get degree of vertex `i`.
"""
@inline function degree(eg::EliminateGraph, vi::Int)
    sum(vj->isconnected(eg,vi,vj), vertices(eg))
end

"""find `(vertex, degree)` with minimum degree of a graph."""
function mindegree_vertex(eg::EliminateGraph)
    N = nv(eg)
    dmin = 999999
    vmin = 0
    for vj in vertices(eg)
        dj = 0
        for vi in vertices(eg)
            dj += isconnected(eg,vi,vj)
        end
        if dj < dmin
            dmin = dj
            vmin = vj
        end
    end
    return vmin, dmin
end

"""find `(min_degree_vertex, max_degree_vertex, min_degree, max_degree)` of a graph."""
function minmaxdegree_vertex(eg::EliminateGraph)
    N = nv(eg)
    dmin = 999999
    dmax = 0
    vmin = 0
    vmax = 0
    for vj in vertices(eg)
        dj = 0
        for vi in vertices(eg)
            dj += isconnected(eg,vi,vj)
        end
        if dj < dmin
            dmin = dj
            vmin = vj
        end
        if dj > dmax
            dmax = dj
            vmax = vj
        end
    end
    return vmin, vmax, dmin, dmax
end

"""find `(vertex, degree)` with maximum degree of a graph."""
function maxdegree_vertex(eg::EliminateGraph)
    N = nv(eg)
    dmax = 0
    vmax = 0
    for vj in vertices(eg)
        dj = 0
        for vi in vertices(eg)
            dj += isconnected(eg,vi,vj)
        end
        if dj > dmax
            dmax = dj
            vmax = vj
        end
    end
    return vmax, dmax
end


text:
abstract type Tile end

struct PolyTile <: Tile
	points
	style
end

polytile(points) = PolyTile(points, x->strokepath())
polytile(points, style) = PolyTile(points, style)

struct NgonTile <: Tile
	n
	radius
	style
end

ngontile(n,r) = NgonTile(n,r, x->strokepath())
ngontile(n,r,style) = NgonTile(n,r, style)

function draw(tile::PolyTile)
	@layer begin
		poly(tile.points, :path, close=true)
		tile.style(0)
	end
end

function fill(tile::PolyTile)
	@layer begin
		poly(tile.points, :fill, close=true)
	end
end

function draw(tile::NgonTile)
	@layer begin
		ngon(0,0,tile.radius,tile.n,0, :path)
		tile.style(0)
	end
end

function hexagonaltile(tile::Tile, radius, range)
	for j in -range:range
		@layer begin
			translate(0, j*sqrt(3)*radius)
			for i in -range:range
				@layer begin
					translate(i*1.5*radius,i*sqrt(3)*0.5*radius)
					if (abs(i+j) <= range) draw(tile) end
					# draw(tile)
				end
			end
		end
	end
end


text:
let

using Atom.FuzzyCompletions

comps(line; mod = "Main", context = "", row = 1, column = 1, force = false) =
    Atom.fuzzycompletionadapter(line, mod, context, row, column, force)
# emurate `getSuggestionDetailsOnSelect` API
add_details!(c) = Atom.completiondetail!(JSON.parse(json(c)))

@testset "module completion" begin
    ## basic
    let
        m = Core.eval(@__MODULE__, :(module Foo end))
        Core.eval(m, quote
            macro foo end
        end)
        cs = comps("@", mod = "Main.Foo")
        @test any(cs) do c
            occursin(c.text, "@foo") &&
            c.type == "snippet" &&
            c.icon == "icon-mention" &&
            c.rightLabel == string(m)
        end
    end

    ## advanced
    let cs = comps("match(code", mod = "Atom")
        @test filter(cs) do c
            c.text == "codeblock_regex" &&
            c.rightLabel == "Atom"
        end |> !isempty
    end

    ## don't error for `ModuleCompletion` of modules where `@doc` isn't defined
    let line = "Core.Compiler.type"
        @test (comps(line); true)
    end
end

# method completion
@testset "method completion" begin
    line = "push!("

    ## basic
    let cs = comps(line)
        # test no error occurs in completion processing
        @test length(cs) == length(FuzzyCompletions.completions(line, lastindex(line))[1])
        # test detecting all the `push!` methods available
        @test filter(cs) do c
            c.type == "method"
        end |> length == length(methods(push!))
    end

    ## advances
    # HACK: this allows testing this in a same julia runtime
    let mod = Main.eval(:(module $(gensym("tempmod")) end))
        # initial state check
        let cs = comps(line)
            filter!(c -> c.text == "push!(::$(mod).Singleton)", cs)
            @test isempty(cs)
        end

        # advanced features
        # - dynamic method addition
        # - shows module where the method defined in right label
        # - shows the infered return type in left label from method signature
        @eval mod begin
            struct Singleton end
            import Base: push!
            push!(::Singleton) = ""
        end
        let cs = comps(line)
            filter!(c -> c.text == "push!(::$(mod).Singleton)", cs)
            # dynamic method addition
            @test length(cs) === 1
            c = cs[1]
            @test c.type == "method"
            # shows module where the method defined in right label
            @test c.rightLabel == string(mod)

            c = add_details!(c)
            @test c["leftLabel"] == "String" # shows the infered return type in left label from method signature
        end

        # shows the infered return type in left label from input types
        @eval mod f(a) = a
        let cs = comps("f("; mod = string(mod))
            filter!(c -> c.text == "f(a)", cs)
            @test length(cs) === 1
            c = cs[1]
            @test c.text == "f(a)"
            @test c.type == "method"
            @test c.rightLabel == string(mod)

            c = add_details!(c)
            @test c["leftLabel"] == "" # don't show a return type when it can't be inferred
        end
        let cs = comps("f(1, "; mod = string(mod))
            filter!(c -> c.text == "f(a)", cs)
            @test length(cs) === 1
            c = cs[1]
            @test c.text == "f(a)"
            @test c.type == "method"
            @test c.rightLabel == string(mod)

            # infer return type with input types
            c = add_details!(c)
            @static if VERSION ≥ v"1.1"
                @test c["leftLabel"] == "$(Int)"
            else
                @test_broken c["leftLabel"] == "$(Int)"
            end
        end

        # don't segfault on generate function
        @eval mod begin
            @generated g(x) = :(return $x)
        end
        @test_nowarn let cs = comps("g("; mod = string(mod))
            @test !isempty(cs)
            @test cs[1].type == "method"
            @test cs[1].text == "g(x)"
        end
    end
end

# HACK: this allows testing these while preventing Main from getting polluted
let m = Main.eval(:(module $(gensym("tempmod")) end))
    ms = string(m)
    @eval m dict = Dict(:a => 1, :b => 2)

    @testset "property completion" begin
        props = string.(propertynames(m.dict))
        foreach(comps("dict.", mod = ms)) do c
            @test c.text in props
            @test c.type == "property"
            @test c.rightLabel == string(typeof(m.dict))
        end
    end

    @testset "field completion" begin
        line = "split(\"\", \"\")[1]."
        fields = string.(fieldnames(SubString{String}))
        foreach(comps(line, mod = ms)) do c
            @test c.text in fields
            @test c.type == "property"
            @test c.rightLabel == "SubString{String}" # the type of object who has those fields
            @test c.leftLabel == string(fieldtype(SubString{String}, Symbol(c.text)))
        end
    end

    @testset "dictionary completion" begin
        ks = repr.(keys(m.dict))
        foreach(comps("dict[", mod = ms)) do c
            @test c.text in ks
            @test c.type == "property"
            @test c.icon == "icon-key"
            @test c.rightLabel == ms
            @test c.leftLabel == string(Int)
        end
    end
end

@testset "keyword completion" begin
    for line in ["begin", "beg", "bgin", "bg"]
        cs = comps(line)
        i = findfirst(c -> c.text == "begin", cs)
        @test i !== nothing
        c = cs[i]
        @test c.type == "keyword"
        @test isempty(c.rightLabel)
        c = add_details!(c)
        @test !isempty(c["description"])
    end
end

@testset "path completion" begin
    line = "\""
    cs = comps(line)
    # XXX: currently `withpath` isn't functional
    # once our path completions are able to return paths relative to each editor's one,
    # we can test against `readdir(@__DIR__)`
    @test length(cs) ==
          length(REPLCompletions.completions(line, lastindex(line))[1]) ==
          length(readdir(pwd()))
    foreach(cs) do c
        @test c.type == "path"
        @test c.icon == "icon-file"
        @test c.rightLabel |> isempty
    end
end

# completion suppressing
troublemakers = [" ", "(", "[", "\$", "=", "&"]
@testset "suppressing: $t" for t in troublemakers
    @test length(comps(t, force = false)) === 0 # surpressed
    @test length(comps(t, force = true)) === Atom.MAX_COMPLETIONS # invoked
end

# don't error on fallback case
@test (Atom.fuzzycompletionadapter(""); true)

@testset "local completions" begin
    type = "attribute"
    icon = "icon-chevron-right"

    let context = """
        function foo(x)
            aaa = 3
            a # C3R6

            foo = x # C5R13
            a # C6R6
            abc = 3
            a # C8R6
        end
        """
        comps′(prefix, row, column, force = false) =
            comps(prefix; context = context, row = row, column = column, force = force)

        # local completions ordered by proximity
        prefix = "a"
        let c = comps′(prefix, 3, 6)[1]
            @test c.text == "aaa"
            @test c.rightLabel == "foo"
            @test c.type == type
            @test c.icon == icon
            @test c.description == "aaa = 3"
        end
        let c = comps′(prefix, 6, 6)[1]
            @test c.text == "aaa"
            @test c.rightLabel == "foo"
            @test c.type == type
            @test c.icon == icon
            @test c.description == "aaa = 3"
        end
        let c = comps′(prefix, 8, 6)[1]
            @test c.text == "abc"
            @test c.rightLabel == "foo"
            @test c.type == type
            @test c.icon == icon
            @test c.description == "abc = 3"
        end

        # show a bound line as is if binding verbatim is not so useful
        let c = comps′("x", 5, 13)[1]
            @test c.text == "x"
            @test c.rightLabel == "foo"
            @test c.type == type
            @test c.icon == icon
            @test c.description == "function foo(x)"
        end

        # show all the local bindings in order when forcibly invoked
        let cs = comps′("", 6, 6, true)
            inds = findall(c -> c.type == type && c.icon == icon, cs)
            @test inds == [1, 2, 3, 4]
            @test map(c -> c.text, cs[inds]) == ["foo", "aaa", "x", "abc"]
        end

        # suppress local completions when in "troublemaker" cases
        for t in troublemakers
            cs = comps′(t, 1, 1)
            @test isempty(cs)
        end
    end

    # don't error on fallback case
    @test (Atom.localcompletions("", 1, 1, ""); true)
end

end


text:
""" Re-zeros the array `a`. """
function rezero!(a)
    @inbounds for i in eachindex(a)
        a[i] = 0.0
    end
end

""" Fill the elements of vector `from` into vector `into`. """
function fill_into!(into, from)
    @inbounds for i in eachindex(into)
        into[i] = from[i]
    end
end

""" Fill vector `into` at indices `inds` with the elements at `from[inds]`. """
function fill_at_inds!(into, from, inds)
    @inbounds for i in 1:length(inds)
        into[inds[i]] = from[i]
    end
end

""" Check whether a point is contained within the simplex currently being checked """
function contained!(signs, s_arr, sx, point, dim)
    # Redefine the temporary simplex. This is in-place, so we don't allocate
    # memory. We could also have re-initialised `signs`, but since we're never
    # comparing more than two consecutive signs, this is not necessary.
    rezero!(s_arr)
    rezero!(signs)
    fill_into!(s_arr, sx)

    #Replace first vertex with the point
    fill_at_inds!(s_arr, point, 1:dim)

    # Signed volume
    signs[1] = sign(det(reshape(s_arr, dim + 1, dim + 1)))

    rezero!(s_arr)
    fill_into!(s_arr, sx) #reset

    for κ = 2:dim # Check remaining signs and stop if sign changes
        # Replace the ith vertex with the point we're cheking (leaving the
        # 1 appended to Vi intact.)
        idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1)
        fill_at_inds!(s_arr, point, idxs) # ith change

        signs[κ] = sign(det(reshape(s_arr, dim + 1, dim + 1)))

        if !(signs[κ-1] == signs[κ])
            return false
        end

        rezero!(s_arr)
        fill_into!(s_arr, sx)
    end

    # Last the last vertex with the point in question
    idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1)
    fill_at_inds!(s_arr, point, idxs)

    signs[end] = sign(det(reshape(s_arr, dim + 1, dim + 1)))

    if !(signs[end-1] == signs[end])
       return false
    else
        return true
    end
end

""" The inner loop in the computation of the transfer matrix. """
function innerloop!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int)
    for j in 1:length(inds)
        if contained!(signs, s_arr, Sj[j], pt, dim)
            M[inds[j], i] += 1.0
        end
    end
end

""" Get the simplices at the given indices. """
function get_simplices_at_inds!(simps, inds::Vector{Int}, simplices)
    for i in 1:length(inds)
        simps[i] = simplices[:, inds[i]]
    end
end


text:
export Curve

"""
`Curve` Represents a domain defined by the image of a Fun
"""

immutable IntervalCurve{S<:Space,T} <: IntervalDomain{T}
    curve::Fun{S,T}
end

immutable PeriodicCurve{S<:Space,T} <: PeriodicDomain{T}
    curve::Fun{S,T}
end

typealias Curve{S,T} Union{IntervalCurve{S,T},PeriodicCurve{S,T}}

==(a::Curve,b::Curve)=a.curve==b.curve

for TYP in (:IntervalCurve,:PeriodicCurve)
    @eval points(c::$TYP,n::Integer)=c.curve(points(domain(c.curve),n))
end

checkpoints(d::Curve) = fromcanonical(d,checkpoints(domain(d.curve)))

for op in (:(Base.first),:(Base.last),:(Base.rand))
    @eval $op(c::Curve)=c.curve($op(domain(c.curve)))
end

canonicaldomain(c::Curve)=domain(c.curve)

fromcanonical{S<:Space,T<:Number}(d::Curve{S,T},v::AbstractArray)=eltype(d)[fromcanonical(d,vk) for vk in v]
fromcanonical{S<:Space,T<:Number}(c::Curve{S,T},x) = c.curve(x)
function tocanonical(c::Curve,x)
    rts=roots(c.curve-x)
    @assert length(rts)==1
    first(rts)
end

fromcanonicalD(c::Curve,x)=differentiate(c.curve)(x)

Base.in(x,d::Curve)=in(tocanonical(d,x),canonicaldomain(d))

Base.reverse(d::Curve)=Curve(reverseorientation(d.curve))

isambiguous(d::Curve)=ncoefficients(d.curve)==0 && isambiguous(domain(d.curve))
Base.convert{S,T}(::Type{IntervalCurve{S,T}},::AnyDomain)=Fun([NaN],S(AnyDomain()))
Base.convert{S,T}(::Type{PeriodicCurve{S,T}},::AnyDomain)=Fun([NaN],S(AnyDomain()))


text:
# # Simple diffusion example
#
# This is Oceananigans.jl's simplest example:
# the diffusion of a one-dimensional Gaussian. This example demonstrates
#
#   * How to load `Oceananigans.jl`.
#   * How to instantiate an `Oceananigans.jl` model.
#   * How to create simple `Oceananigans.jl` output.
#   * How to set an initial condition with a function.
#   * How to time-step a model forward.
#   * How to look at results.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## Using `Oceananigans.jl`
#
# Write

using Oceananigans

# to load Oceananigans functions and objects into our script.
#
# ## Instantiating and configuring a model
#
# A core Oceananigans type is `NonhydrostaticModel`. We build an `NonhydrostaticModel`
# by passing it a `grid`, plus information about the equations we would like to solve.
#
# Below, we build a regular rectilinear grid with 128 grid points in the `z`-direction,
# where `z` spans from `z = -0.5` to `z = 0.5`,

grid = RegularRectilinearGrid(size=128, z=(-0.5, 0.5), topology=(Flat, Flat, Bounded))

# The default topology is `(Periodic, Periodic, Bounded)` but since we only want to solve
# a one-dimensional problem, we assign the `x` and `y` dimensions to `Flat`.
# We could specify each of them to be either `Periodic` or `Bounded` but that will define
# a halo in each of those directions, and that is numerically more costly.
# Note that we only specify the extent and size for the `Bounded` dimension.
#
# We next specify a model with an `IsotropicDiffusivity`, which models either
# molecular or turbulent diffusion,

closure = IsotropicDiffusivity(κ=1.0)

# We finally pass these two ingredients to `NonhydrostaticModel`,

model = NonhydrostaticModel(grid=grid, closure=closure)

# Our simple `grid` and `model` use a number of defaults:
#
#   * The default `grid` topology is periodic in `x, y` and bounded in `z`.
#   * The default `Model` has no-flux (insulating and stress-free) boundary conditions on
#     non-periodic boundaries for velocities `u, v, w` and tracers.
#   * The default `Model` has two tracers: temperature `T`, and salinity `S`.
#   * The default `Model` uses a `SeawaterBuoyancy` model with a `LinearEquationOfState`.
#     However, buoyancy is not active in the simulation we run below.
#
# Next, we `set!` an initial condition on the temperature field,
# `model.tracers.T`. Our objective is to observe the diffusion of a Gaussian.

width = 0.1

initial_temperature(x, y, z) = exp(-z^2 / (2width^2))

set!(model, T=initial_temperature)

# ## Visualizing model data
#
# Calling `set!` above changes the data contained in `model.tracers.T`,
# which was initialized as `0`'s when the model was created.
# To see the new data in `model.tracers.T`, we plot it:

using Plots

z = znodes(model.tracers.T)

T_plot = plot(interior(model.tracers.T)[1, 1, :], z,
              linewidth = 2,
              label = "t = 0",
              xlabel = "Temperature (ᵒC)",
              ylabel = "z")

# The function `interior` above extracts a `view` of the physical interior points
# of `model.tracers.T`. This is useful because `model.tracers.T` also contains "halo" points
# that lie outside the physical domain (halo points are used to set boundary conditions
# during time-stepping).

# ## Running a `Simulation`
#
# Next we set-up a `Simulation` that time-steps the model forward and manages output.

## Time-scale for diffusion across a grid cell
diffusion_time_scale = model.grid.Δz^2 / model.closure.κ.T

simulation = Simulation(model, Δt = 0.1 * diffusion_time_scale, stop_iteration = 1000)

# `simulation` will run for 1000 iterations with a time-step that resolves the time-scale
# at which our temperature field diffuses. All that's left is to

run!(simulation)

# ## Visualizing the results
#
# Let's look at how `model.tracers.T` changed during the simulation.

using Printf

plot!(T_plot, interior(model.tracers.T)[1, 1, :], z, linewidth=2,
      label=@sprintf("t = %.3f", model.clock.time))

# Very interesting! Next, we run the simulation a bit longer and make an animation.
# For this, we use the `JLD2OutputWriter` to write data to disk as the simulation progresses.

using Oceananigans.OutputWriters: JLD2OutputWriter, IterationInterval

simulation.output_writers[:temperature] =
    JLD2OutputWriter(model, model.tracers, prefix = "one_dimensional_diffusion",
                     schedule=IterationInterval(100), force = true)

# We run the simulation for 10,000 more iterations,

simulation.stop_iteration += 10000

run!(simulation)

# Finally, we animate the results by opening the JLD2 file, extract the
# iterations we ended up saving at, and plot the evolution of the
# temperature profile in a loop over the iterations.

using JLD2

file = jldopen(simulation.output_writers[:temperature].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

anim = @animate for (i, iter) in enumerate(iterations)

    T = file["timeseries/T/$iter"][1, 1, :]
    t = file["timeseries/t/$iter"]

    plot(T, z, linewidth=2, title=@sprintf("t = %.3f", t),
         label="", xlabel="Temperature", ylabel="z", xlims=(0, 1))
end

mp4(anim, "one_dimensional_diffusion.mp4", fps = 15) # hide


text:
module SegyIO

    myRoot=dirname(dirname(pathof(SegyIO)))
    CHUNKSIZE=2048
    MB2B=1024^2

    # what's being used
    using Distributed
    using Printf

    #Types
    include("types/IBMFloat32.jl")
    include("types/BinaryFileHeader.jl")
    include("types/BinaryTraceHeader.jl")
    include("types/FileHeader.jl")
    include("types/SeisBlock.jl")
    include("types/BlockScan.jl")
    include("types/SeisCon.jl")

    #Reader
    include("read/read_fileheader.jl")
    include("read/read_traceheader.jl")
    include("read/read_trace.jl")
    include("read/read_file.jl")
    include("read/segy_read.jl")
    include("read/read_block.jl")
    include("read/read_block_headers.jl")
    include("read/read_con.jl")
    include("read/read_con_headers.jl")
    include("read/extract_con_headers.jl")

    # Writer
    include("write/write_fileheader.jl")
    include("write/segy_write.jl")
    include("write/segy_write_append.jl")
    include("write/check_fileheader.jl")
    include("write/write_trace.jl")

    # Scanner
    include("scan/scan_file.jl")
    include("scan/segy_scan.jl")
    include("scan/scan_chunk.jl")
    include("scan/scan_block.jl")
    include("scan/scan_shots.jl")
    include("scan/delim_vector.jl")
    include("scan/find_next_delim.jl")

    # Workspace
    th_b2s = SegyIO.th_byte2sample()
    blank_th = BinaryTraceHeader()

    # Methods
    include("methods/ordered_pmap.jl")
    include("methods/merge.jl")
    include("methods/split.jl")
    include("methods/set_header.jl")
    include("methods/get_header.jl")
    include("methods/get_sources.jl")

end # module


text:
#  test MultipleFeatures2D
using RoME, Distributions
using IncrementalInference   # for evalFactor2
using Base.Test

using Optim

# using Gadfly

x1 = [0.0;0;0]
x2 = [1.0;-1.0;pi/2]

l1 = [0.5;0.5]
l2 = [1.0;0.0]
l3 = [1.5;-0.5]
l3b = [1.5;0.0]
L = [l1';l2';l3']

# plot(x=L[:,1],y=L[:,2],Geom.point)

the11 = atan2(l1[2],l1[1])
the12 = atan2(l2[2],l2[1])
the13 = atan2(l3[2],l3[1])

Xi = SE2(x1)
Xj = SE2(x2)
Δxij = se2vee(Xi\Xj)

tl1 = se2vee(Xj \ SE2([l1;0.0]))[1:2]
the21 = atan2(tl1[2],tl1[1])

tl2 = se2vee(Xj \ SE2([l2;0.0]))[1:2]
the22 = atan2(tl2[2],tl2[1])

tl3 = se2vee(Xj \ SE2([l3;0.0]))[1:2]
the23 = atan2(tl3[2],tl3[1])

tl3b = se2vee(Xj \ SE2([l3b;0.0]))[1:2]
the23b = atan2(tl3b[2],tl3b[1])

xir1 = Normal(the11,1e-4)
xir2 = Normal(the12,1e-4)
xir3 = Normal(the13,1e-4)
xjr1 = Normal(the21,1e-4)
xjr2 = Normal(the22,1e-4)
xjr3 = Normal(the23,1e-4)
xjr3b = Normal(the23b,1e-4)

# Not using second feature
mm2 = MultipleFeatures2D(
  xir1,
  xir2,
  xir3,
  xjr1,
  xjr2,
  xjr3,
  xjr3b,
  Categorical([1.0;0.0]),
  SE2(zeros(3))  )

res = randn(3)
idx = 1
meas = getSample(mm2)

wAbi = reshape(x1,3,1)
wAbj = reshape(x2,3,1)
wAo1 = reshape(l1,2,1)
wAo2 = reshape(l2,2,1)
wAo3 = reshape(l3,2,1)

println("test MultipleFeatures unit vector functions")
res, rhat, resid = zeros(3),zeros(2), zeros(2)

for i in 1:3
  r1, rhathat1, α1 = getUvecScaleFeature2D(Xi, SE2(zeros(3)), SE2([L[i,:];0.0]))
  rhat[1] = cos(meas[1][i,idx,1])
  rhat[2] = sin(meas[1][i,idx,1])
  # @show rhat, rhathat1
  resid[1:2] = rhat-rhathat1
  res[1] += dot(resid, resid)
end

@test 0.0 <= res[1] < 1e-5

println("test residual function")

res = zeros(3)

wAbi[:] = 0.0
# @show wAbi

mm2(res, idx, meas,
  wAbi,
  wAbj,
  wAo1,
  wAo2,
  wAo3  )

# @show res

function gg2(x::Float64,y::Float64)
  wAbi[2] = x
  wAbi[3] = y

  res = zeros(3)

  mm2(res, idx, meas,
    wAbi,
    wAbj,
    wAo1,
    wAo2,
    wAo3  )
  return sum(res)[1]
end
XX = linspace(-0.5,0.5,200)
YY = linspace(-pi/4,pi/4,100)
gg2(0.0,0.0)
# Gadfly.plot(z=gg2,x=XX,y=YY, Geom.contour)

# ar99
println("test as minimization problem to pose")
function minmickey(x::Array{Float64,1})
  wAbi[1:3,1] = x[1:3]
  res = zeros(3)
  mm2(res, idx, meas,
    wAbi,
    wAbj,
    wAo1,
    wAo2,
    wAo3  )

  return sum(res[1:2])
end

gg = (x,y) -> minmickey([x;y;0.0])
# gg = (x,y) -> x*y
# plot(z=gg, x=linspace(-2,2,100),y=linspace(-2,2,100), Geom.contour)

gg(0.0,0.1)

@time r = optimize(minmickey, randn(3))
@show r.minimizer

@test norm(res) < 1e-3

println("test as minimization problem to landmark")
wAbi[:] = 0.0
function minmickey(x::Array{Float64,1})
  wAo1[1:2,1] = x[1:2]
  res = zeros(3)
  mm2(res, idx, meas,
    wAbi,
    wAbj,
    wAo1,
    wAo2,
    wAo3  )

  return sum(res[1:2])
end

gg = (x,y) -> minmickey([x;y])
#plot(z=gg, x=linspace(1,2,100),y=linspace(-1,1,100), Geom.contour(levels=50))

# now build in factor graph form for further testing

N=50
fg = initfg()

initCov = 0.05*eye(3)
initCov[3,3] = 0.001
odoCov = deepcopy(initCov)

ipp = PriorPose2(zeros(3,1), initCov^2,[1.0])

v1 = addNode!(fg, :x1, 0.01*randn(3,N), diagm([1.0;1.0;0.1]), N=N)

f1  = addFactor!(fg,[:x1], ipp)

pts2 = [0.014*randn(1,N)+1;  0.014*randn(1,N)-1; 0.002*randn(1,N)+pi/2]

v2 = addNode!(fg, :x2, pts2, diagm([1.0;1.0;0.1]), N=N)

px2 = zeros(3,1)
px2[1:3,1] = [1.1;-1.0;pi/2]
ipp2 = PriorPose2(px2, initCov^2,[1.0])
f1  = addFactor!(fg,[:x2], ipp2)

vl1 = addNode!(fg, :l1, rand(MvNormal(l1,0.001*eye(2)),N), diagm([1.0;1.0]), N=N)
vl2 = addNode!(fg, :l2, rand(MvNormal(l2,0.001*eye(2)),N), diagm([1.0;1.0]), N=N)
vl3 = addNode!(fg, :l3, rand(MvNormal(l3,0.001*eye(2)),N), diagm([1.0;1.0]), N=N)

f2 = addFactor!(fg, [v1;v2;vl1;vl2;vl3], mm2)

ef2pts = evalFactor2(fg, f2, fg.IDs[:l2])

tree = wipeBuildNewTree!(fg)
# spyCliqMat(tree.cliques[1])

[inferOverTreeR!(fg,tree, N=N, dbg=true) for i in 1:1]

# lets see what is happening during MCMC runs
# plotMCMC(tree, :l2, show=true, levels=3) # make true to see pictures, false for testing

# spyCliqMat(tree, :x2)

#
# pts, = getSample(ipp2, 100);
#
# plotKDE( kde!(pts[1:2,:]) )
#
# @show ipp2.Cov
# pts = rand(MvNormal(ipp2.Zi[:,1],ipp2.Cov),N);

# using KernelDensityEstimate
# plotKDE(kde!(ef2pts[1,:]))

# pl = plotKDE([kde!(randn(2,100));kde!(2+randn(2,100))],levels=1,c=["red";"green"],legend=["1";"2"])

# arrp = [getVertKDE(fg, :l1);getVertKDE(fg, :l2);getVertKDE(fg, :l3)]
# arrp1 = deepcopy(arrp)
# plotKDE([arrp;arrp1],levels=1,c=["red";"red";"red";"blue";"blue";"blue"])

# arrp=[marginal(getVertKDE(fg, :x1),[1;2]);marginal(getVertKDE(fg, :x2),[1;2]) ]
# plotKDE(arrp,levels=1)

# plotKDE(marginal(getVertKDE(fg,:x1),[3]))

# pp,arr,parts = localProduct(fg,:x2)
# cc = plotKDE([pp;arr],c=["red";"black";"blue"],levels=1);

# cc = drawLbl(fg,:l1);
# Gadfly.draw(PDF("/home/dehann/Desktop/test.pdf",20cm,15cm),cc)

println("test ambiguous bi-modal multifeature constraint operation")

mm2 = MultipleFeatures2D(
  xir1,
  xir2,
  xir3,
  xjr1,
  xjr2,
  xjr3,
  xjr3b,
  Categorical([0.5;0.5]),
  SE2(zeros(3))  )

#


text:
function clusterfit(
    X::AbstractMatrix{<:Real},
    centers::AbstractMatrix{<:AbstractFloat}
)
    d, n = size(X)

    assignments = Vector{Int}(undef, n)
    dist = pairwise(Euclidean(), centers, X, dims = 2)
    D = eltype(dist)

    costs = Vector{D}(undef, n)
    counts = Vector{Int}(undef, size(centers, 2))

    update_assignments!(dist, assignments, costs, counts)

    return assignments
end


text:
using Behavior
using Behavior.Gherkin

rootpath = joinpath(@__DIR__, "test/behaviour")
featurepath = joinpath(@__DIR__, "test/behaviour/features/connection")
stepspath = joinpath(@__DIR__,"test/behaviour/connection")
configpath = joinpath(@__DIR__,"test/behaviour/config/ConfigEnvironment.jl")

p = ParseOptions(allow_any_step_order = true)

runspec(rootpath; featurepath = featurepath, stepspath = stepspath,  parseoptions=p, execenvpath = configpath, tags="not @ignore-typedb-core")

runspec(rootpath; featurepath = featurepath, stepspath = stepspath,  parseoptions=p, execenvpath = configpath, tags="@actual")


text:
struct DistanceSphereToLine{M,R1,S,m} <: TaskMap{M,R1,S}
    radius::S
    point1::SVector{m,S}
    point2::SVector{m,S}
end
DistanceSphereToLine{M,R1}(r::S, p1, p2) where {M,R1,S} =
    DistanceSphereToLine{M,R1,S,dim(M)}(r, p1, p2)

struct DistanceSphereToLineT{M,R1,S,m} <: TaskMapT{M,R1,S}
    base_map::DistanceSphereToLine{M,R1,S,m}
end

DistanceSphereToLineT{M,R1}(r::S, p1, p2) where {M,R1,S} =
    DistanceSphereToLineT{M,R1,S,m}(DistanceSphereToLine{M,R1}(r, p1, p2))

function task_map_emb(::EmbRep, ::EmbRep, xme, task_map::DistanceSphereToLine{ℝ{m},R1,S}) where {m,S}
    point1, point2 = task_map.point1, task_map.point2
    p1p2 = point2 - point1
    p1x = xme - point1
    p2x = xme - point2
    if dot(p1x, p1p2) > 0.
        if dot(p2x, p1p2) < 0.
            normal_vec = p1x - p1p2*dot(p1x, p1p2)/norm(p1p2)^2
            return SA[norm(normal_vec) - task_map.radius]
        else
            return SA[norm(p2x) - task_map.radius]
        end
    else
        return SA[norm(p1x) - task_map.radius]
    end
end


text:
## Based on the work of Brunton, et al
## B. W. Brunton, S. L. Brunton, J. L. Proctor, J. N. Kutz. Sparse sensor placement optimization for classification.  *SIAM Journal on Applied Mathematics*, 76(5):2099-2122, 2016.
##

module SSPOC

const DIR = @__DIR__
using Pkg
Pkg.activate(DIR*"/../")

using DelimitedFiles, LinearAlgebra, Statistics, Convex, COSMO, StatsBase, Flux
using Flux.Data:DataLoader
using Flux.Losses: logitbinarycrossentropy
using Flux: onecold, onehot
using Flux: onehotbatch
using StatsBase, PyCall

plt = pyimport("matplotlib.pyplot")

cats = Float32.(readdlm(DIR*"/../data/cats.csv", ',')) #|> x -> x ./ 255
dogs = Float32.(readdlm(DIR*"/../data/dogs.csv", ',')) #|> x -> x ./ 255

function run_demo(cats, dogs, n, r)
    out=zeros(n)
    for i = 1:n
        @label start
        cats = deepcopy(cats[:, sample(1:size(cats, 2), size(cats, 2), replace=false)]);
        dogs = deepcopy(dogs[:, sample(1:size(dogs, 2), size(dogs, 2), replace=false)]);

        # create training and testing sets (.74, .25)
        train1 = [cats[:,1:60] dogs[:,1:60]];
        test1 = [cats[:,61:end] dogs[:,61:end]];

        # compute PCA and project data to get eigenrepresentation X̃
        n, m = size(train1)
        ## centering matrix
        C = I(m) - (1/m) * ones(m,m)
        train_centered = train1 * C
        U, S, V = svd(train_centered)
        # plt.plot(s)
        # plt.show()
        Ψ = U[:,1:r];
        X̃ = Ψ' * train1;

        # compute maximally discriminating component using LDA
        ## between-classes scatter Sb
        μ1 = mean(X̃[:, 1:60], dims=2);
        μ2 = mean(X̃[:,61:end], dims=2);
        Sb = (μ2 - μ1) * (μ2 - μ1)';

        function scatter_matrix(data, mu)
            (data - mu .* ones(size(data))) * (data - mu .* ones(size(data)))'
        end
        ## within-classes scatter Sw
        Σ1 = scatter_matrix(X̃[:,1:60], μ1);
        Σ2 = scatter_matrix(X̃[:,61:end], μ2);
        Sw = Σ1 + Σ2;

        # generalized eigenvalue decomposition to get projection matrix W
        Λ, W = eigen(Sb, Sw);
        ## choose eigenvector corresponding to largest eigenvalue w (this represents the maximally discriminating projection vector in the new subspace)
        Λ = reverse(real.(Λ));
        W = reverse(W, dims=2);
        w = real(W[:,1])

        # solve for S (this is the maximally sparse vector in original space that when projected to subspace approximates w)
        n = size(train1, 1)
        s = Variable(n)
        objective = norm(s,1)
        constraint = Ψ' * s == w;
        solver = COSMO.Optimizer;
        problem = minimize(objective, constraint);
        solve!(problem, solver, silent_solver=true)
        ## check linear separability imposed by sparse vector S
        if isnothing(s.value)
            @warn "Encountered a numerical error during optimization!  Repeating sample..."
            num_err += 1
            @goto start
        end
        S = reshape(s.value, n)
        sep = (Ψ' * S)' * (Ψ' * train1)
        # plt.bar(1:length(sep), sep)
        # plt.show()

        # identify indices above threshold
        ## see how close to sparse the sparse array is
        #StatsBase.describe(S)
        ## choose threshold
        # norm(S, 2) / (2*2*60)
        # C = (abs.(S) .> 0.35)
        # n = norm(S,2) / (2 * 2 * 60)
        # C = (abs.(S) .> n)
        thresh = abs(S[reverse(sortperm(abs.(S)))[21]])
        C = abs.(S) .> thresh;
        #C = zeros(length(S))
        #sum(C)

        # create measurement matrix
        Φ = Diagonal(C);

        # scale data (z score transform)
        dt = StatsBase.fit(ZScoreTransform, train1, dims=2);
        StatsBase.transform!(dt, train1);

        # sparsify original data with measurement matrix
        #train2 = train1 + .1 * randn(size(train1))
        y = Φ * train1;
        #y = Φ * hcat(train1, train2)

        ## see what pixels will be used
        # function pl(i)
        # cat_sparse_pixels = deepcopy(cats[:, i])
        # dog_sparse_pixels = deepcopy(dogs[:, i])
        # eig_sparse_pixels = deepcopy(Ψ[:,1])
        # cat_sparse_pixels[C] .= NaN
        # eig_sparse_pixels[C] .= NaN
        # dog_sparse_pixels[C] .= NaN
        # fig, (ax1, ax2, ax3) = plt.subplots(1,3)
        # cmap = plt.cm.gray
        # cmap.set_bad((1, 0, 0, 1))
        # ax1.imshow(reshape(cat_sparse_pixels, (64,64)), cmap=cmap)
        # ax3.imshow(reshape(eig_sparse_pixels, (64,64)), cmap=cmap)
        # ax2.imshow(reshape(dog_sparse_pixels, (64,64)), cmap=cmap)
        # plt.show()
        # end

        # pl(36)

        # create labels
        #labels = Int.(vcat(ones(60), zeros(60)));
        labels = Int.(vcat(ones(60), zeros(60)))
        l = onehotbatch(labels, [0, 1]);

        # prepare data for NN
        train = DataLoader((y, l), batchsize=32);

        # build logistic regression model with L1 regularization

        function loss(x, y)
            logitbinarycrossentropy(m(x), y) + .1 * norm(params(m), 1)
            #logitbinarycrossentropy(m(x), y)
        end

        # optimizer
        #opt = ADAM(1e-3)
        # opt = Descent(1e-2)
        opt = ADADelta(.5)
        m = Dense(4096, 2)
        # train for 100 epochs
        for epoch in 1:100
            Flux.train!(loss, params(m), train, opt)
        end

        # prepare test data
        test1 = [cats[:,61:end] dogs[:,61:end]];
        StatsBase.transform!(dt, test1);
        y_test = Φ * test1;
        labels_test = Int.(vcat(ones(20), zeros(20)));

        out[i] = mean(onecold(sigmoid.(m(y_test)), [0,1]) .== labels_test)
        @info "Fold "*string(i)*" accuracy: "*string(out[i])
    end
mean(out)
end #run_demo

end # module


text:
"""
Test intances collected mainly from various papers about 3D Dubins paths.
"""

rhomin = 40.
pitchlims = pi * [-15., 20.] / 180.
deg2rad(x) = pi * x / 180.

LONG = [
  [
    [200, 500, 200, deg2rad(180), deg2rad(-5)],
    [500, 350, 100, deg2rad(0), deg2rad(-5)],
    [467.70, 449.56],
    "Long 1"
  ],[
    [100, -400, 100, deg2rad(30), deg2rad(0)],
    [500, -700, 0, deg2rad(150), deg2rad(0)],
    [649.52, 636.12],
    "Long 2"
  ],[
    [-200, 200, 250, deg2rad(240), deg2rad(15)],
    [ 500, 800,   0, deg2rad(45), deg2rad(15)],
    [1088.10, 1063.41],
    "Long 3"
  ],[
    [-300, 1200, 350, deg2rad(160), deg2rad(0)],
    [1000,  200,   0, deg2rad(30), deg2rad(0)],
    [1802.60, 1789.21],
    "Long 4"
  ],[
    [-500, -300, 600, deg2rad(150), deg2rad(10)],
    [1200,  900, 100, deg2rad(300), deg2rad(10)],
    [2245.14, 2216.40],
    "Long 5"
  ]
]

SHORT = [
  [
    [120, -30, 250, deg2rad(100), deg2rad(-10)],
    [220, 150, 100, deg2rad(300), deg2rad(-10)],
    [588.60, 583.47],
    "Short 1"
  ],[
    [380, 230, 200, deg2rad(30), deg2rad(0)],
    [280, 150,  30, deg2rad(200), deg2rad(0)],
    [667.71, 658.53],
    "Short 2"
  ],[
    [-80, 10, 250, deg2rad(20), deg2rad(0)],
    [ 50, 70,   0, deg2rad(240), deg2rad(0)],
    [979.34, 968.25],
    "Short 3"
  ],[
    [400, -250, 600, deg2rad(350), deg2rad(0)],
    [600, -150, 300, deg2rad(150), deg2rad(0)],
    [1169.73, 1161.55],
    "Short 4"
  ],[
    [-200, -200, 450, deg2rad(340), deg2rad(0)],
    [-300,  -80, 100, deg2rad(100), deg2rad(0)],
    [1367.56, 1354.12],
    "Short 5"
  ]
]

BEARD = [
  [
    [   0,   0, 350, deg2rad(180)+pi/2,  0],
    [-100, 100, 100, deg2rad(180)+pi/2,  0],
    [NaN],
    "Beard 1"
  ], [
    [  0,   0,  100, deg2rad(70)+pi/2,  0],
    [100, 100,  125, deg2rad(70)+pi/2,  0],
    [NaN],
    "Beard 2"
  ], [
    [  0, 0, 100, deg2rad(0)+pi/2,   0,  0],
    [200, 0, 200, deg2rad(270)-pi/2, 0,  0],
    [NaN],
    "Beard 3"
  ]
]

OTHERS = [
  [
    [500, 100, 300, deg2rad(240), deg2rad(15)],
    [-100, 400, 0, deg2rad(45), deg2rad(15)],
    [NaN],
    "Others 1"
  ], [
    [0., 0, 0, 0, 0],
    [0., 0, 5, 0, 0],
    [NaN],
    "Others 2"
  ], [
    [0., 0, 0, 0, 0],
    [0., 0, 10, 0, 0],
    [NaN],
    "Others 3"
  ], [
    [0., 0, 0, 0, 0],
    [0., 0, 20, 0, 0],
    [NaN],
    "Others 4"
  ], [
    [0., 0, 0, 0, 0],
    [0., 0, 30, 0, 0],
    [NaN],
    "Others 5"
  ], [
    [0., 0, 0, 0, 0],
    [0., 0, 50, 0, 0],
    [NaN],
    "Others 6"
  ], [
    [0., 0, 0, 0, deg2rad(10)],
    [292.4, 0, 100, 0, deg2rad(10)],
    [NaN],
    "Others 7"
  ], [
    [0., 0, 200, 0, 0],
    [0., 0, 20, 0, 0],
    [NaN],
    "Others 8"
  ]
]

HOTA = [
  [
    [100, 10, 200, deg2rad(180), deg2rad(20)],
    [500, 40, 100, deg2rad(0), deg2rad(20)],
    [NaN],
    "Hota 1"
  ], [
    [100, 10, 200, deg2rad(180), deg2rad(20)],
    [120, 20, 210, deg2rad(90), deg2rad(20)],
    [NaN],
    "Hota 2"
  ]
]


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

using ArgParse

s = ArgParseSettings()
@add_arg_table! s begin
    "--conf"
        help = "configuration setup"
        arg_type = String
        default = "2a"
end
parsed_args = parse_args(ARGS, s)
setup = parsed_args["conf"]

using HMMGradients, Flux, Zygote
using Random, Statistics, LinearAlgebra
using FiniteStateTransducers
using DSP, MFCC
using BSON, JLD2, LibSndFile, FileIO, UUIDs
include("WFSTs.jl")
include("Models.jl")
include("Utils.jl")

include("conf/$(setup)/feat_conf.jl")

# get transition matrix
lexicon, ilexicon = get_lexicon()
H, L = get_HL(lexicon)
a, A = get_aA(H)
ippsym = get_iisym(H)

BSON.@load "models/$setup/best_modely_final.bson" best_modely
Flux.testmode!(best_modely)

feat_dir = get_feat_dir(setup)
data = load(joinpath(feat_dir,"train.jld2"))
uttID2feats, uttID2text, uttID2phones = data["uttID2feats"], data["uttID2text"], data["uttID2phones"]

uttID2text_dec   = Dict()
uttID2phones_dec = Dict()
min_dist=2

for uttID in keys(uttID2feats)
  x = uttID2feats[uttID]
  x = feats_post(x)
  y = best_modely(Flux.unsqueeze(x,3))
  gamma = logposterior(size(y,1),a,A,y[:,:])
  ps  = posterior2phones(ippsym,gamma)
  ws  = phones2words_greedy(ilexicon,ps; min_dist=min_dist)
  uttID2phones_dec[uttID] = ps
  uttID2text_dec[uttID] = strip(prod(ws.*" "))
end

uttID2err_textdec = Dict{String,Tuple{String,String}}()
for uttID in keys(uttID2text)
  text, dec = uttID2text[uttID], uttID2text_dec[uttID]
  if text != dec
    uttID2err_textdec[uttID] = (text,dec)
  end
end

accuracy = 1-length(uttID2err_textdec) / length(uttID2text)
wer = get_error_rate(uttID2text, uttID2text_dec)
per = get_error_rate(uttID2phones, uttID2phones_dec; is_split=true)
println("# Setup $setup")
println("* Phone Error Rate (PER): $(round(per * 100, digits=3)) %")
println("* Word Error Rate (WER) : $(round(wer * 100, digits=3)) %")
println("* Accuracy: $(round(accuracy,digits=3))")


text:
# Digital signal plots
function PlotDigital(x::AbstractArray{T}, y::AbstractArray{T}; label_id::String="",
                     count::Integer=min(length(x), length(y)), offset::Integer=0,
                     stride::Integer=1) where {T<:ImPlotData}
    return PlotDigital(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotDigital(x::AbstractArray{T1}, y::AbstractArray{T2};
                     kwargs...) where {T1<:Real,T2<:Real}
    return PlotDigital(promote(x, y)...; kwargs...)
end

# xfield, yfield should be propertynames of eltype(structvec)
function PlotDigital(structvec::Vector{T}, xfield::Symbol, yfield::Symbol;
                     count::Integer=length(structvec), offset::Integer=0, stride::Integer=1,
                     label_id::String="") where {T}
    Tx = fieldtype(T, xfield)
    Ty = fieldtype(T, yfield)
    x_offset = fieldoffset(T, Base.fieldindex(T, xfield))
    y_offset = fieldoffset(T, Base.fieldindex(T, yfield))
    x_ptr = Ptr{Tx}((pointer(structvec, 1) + x_offset))
    y_ptr = Ptr{Ty}((pointer(structvec, 1) + y_offset))

    if !T.mutable
        # this is somewhat illegal and is used only to pass a pointer through AbstractArray argument into ccall
        x = unsafe_wrap(Vector{Tx}, x_ptr, size(structvec); own=false)
        y = unsafe_wrap(Vector{Ty}, y_ptr, size(structvec); own=false)
        stride = stride * sizeof(T)
    else # two new vectors every 1/60 second...
        x = Vector{Tx}(undef, length(structvec))
        y = Vector{Ty}(undef, length(structvec))
        for (i, val) in enumerate(structvec)
            x[i] = getproperty(val, xfield)
            y[i] = getproperty(val, yfield)
        end
        if Tx !== Ty
            x, y = promote(x, y)
        end
        stride = stride * sizeof(eltype(x))
    end

    return PlotDigital(label_id, x, y, count, offset, stride)
end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

using StaticArrays

const _x1 = Grids._x1
const _x2 = Grids._x2
const _x3 = Grids._x3

"""
    knl_min_neighbor_distance!(::Val{N}, ::Val{dim}, direction,
                             min_neighbor_distance, vgeo, topology.realelems)

Computational kernel: Computes the minimum physical distance between node
neighbors within an element.

The `direction` in the reference element controls which nodes are considered
neighbors.
"""
function knl_min_neighbor_distance!(::Val{N}, ::Val{dim}, direction,
                                    min_neighbor_distance, vgeo, elems
                                   ) where {N, dim}

  FT = eltype(min_neighbor_distance)
  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  if direction isa EveryDirection
    mininξ = (true, true, true)
  elseif direction isa HorizontalDirection
    mininξ = (true, dim == 2 ? false : true, false)
  elseif direction isa VerticalDirection
    mininξ = (false, dim == 2 ? true : false, dim == 2 ? false : true)
  end

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          md = typemax(FT)

          ijk = i + Nq * (j - 1) + Nq * Nq * (k - 1)
          x = SVector(vgeo[ijk, _x1, e], vgeo[ijk, _x2, e], vgeo[ijk, _x3, e])

          if mininξ[1]
            @unroll for î = (i-1, i+1)
              if 1 ≤ î ≤ Nq
                îjk = î + Nq * (j-1) + Nq * Nq * (k-1)
                x̂ = SVector(vgeo[îjk, _x1, e],
                            vgeo[îjk, _x2, e],
                            vgeo[îjk, _x3, e])
                md = min(md, norm(x - x̂))
              end
            end
          end

          if mininξ[2]
            @unroll for ĵ = (j-1, j+1)
              if 1 ≤ ĵ ≤ Nq
                iĵk = i + Nq * (ĵ-1) + Nq * Nq * (k-1)
                x̂ = SVector(vgeo[iĵk, _x1, e],
                            vgeo[iĵk, _x2, e],
                            vgeo[iĵk, _x3, e])
                md = min(md, norm(x - x̂))
              end
            end
          end

          if mininξ[3]
            @unroll for k̂ = (k-1, k+1)
              if 1 ≤ k̂ ≤ Nqk
                ijk̂ = i + Nq * (j-1) + Nq * Nq * (k̂-1)
                x̂ = SVector(vgeo[ijk̂, _x1, e],
                            vgeo[ijk̂, _x2, e],
                            vgeo[ijk̂, _x3, e])
                md = min(md, norm(x - x̂))
              end
            end
          end

          min_neighbor_distance[ijk, e] = md
        end
      end
    end
  end
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
function meantoecc(M::FloatingPoint, ecc::FloatingPoint)
    kepler(E) = E - ecc*sin(E) - M
    kepler_der(E) = 1 - ecc*cos(E)
    return newton(M, kepler, kepler_der)
end

function ecctomean(E::FloatingPoint, ecc::FloatingPoint)
    return E - ecc*sin(E)
end

function ecctotrue(E::FloatingPoint, ecc::FloatingPoint)
    return 2*atan2(sqrt(1 + ecc)*sin(E/2), sqrt(1 - ecc)*cos(E/2))
end

function truetoecc(T::FloatingPoint, ecc::FloatingPoint)
    return 2*atan2(sqrt(1 - ecc)*sin(T/2), sqrt(1 + ecc)*cos(T/2))
end

function period(a::FloatingPoint, mu::FloatingPoint)
    return sqrt(4*a^3*pi^2/mu)
end

function period(s::State)
    mu = planets[s.body]["mu"]
    ele = elements(s)
    return period(ele[1], mu)
end

function kepler(ele::Vector, dt::FloatingPoint, mu::FloatingPoint)
    E0 = truetoecc(ele[6], ele[2])
    M0 = ecctomean(E0, ele[2])
    n = 2*pi/period(ele[1], mu)
    M = M0 + n*dt
    E = meantoecc(M, ele[2])
    T = ecctotrue(E, ele[2])
    return [ele[1:5], T]
end

function kepler(ele::Vector, dt::Vector, mu::FloatingPoint)
    out = Array(FloatingPoint, length(dt), 6)
    for i = 1:length(dt)
        out[i,:] = kepler(ele, dt[i], mu)
    end
    return out
end


text:
using ASE
using JuLIP
using BenchmarkTools
using Main.Threads
using DataDeps
using Profile
using ProfileView

register(DataDep(
    "JuLIP_testpots",
    "A few EAM potentials for testing",
    "https://www.dropbox.com/s/leub1c9ft1mm9fg/JuLIP_data.zip?dl=1",
    post_fetch_method = file -> run(`unzip $file`)
    ))
test_pots = joinpath(datadep"JuLIP_testpots", "JuLIP_data") * "/"
Ni = test_pots * "Ni.eam.fs"

eam0 = EAM(Ni)
eam1 = EAM(Ni; s = 1e-3)
display(eam0.ϕ[1].f.spl)
display(eam1.ϕ[1].f.spl)

eam = eam1

@show Main.Threads.nthreads()
@show JuLIP.nthreads()

function evaln(eam, at, N)
    E = zeros(Main.Threads.nthreads())
    F = [ forces(eam, at) * 0  for _ = 1:Main.Threads.nthreads() ]
    @threads for n = 1:N
        tid = threadid()
        at1 = deepcopy(at)
        rattle!(at, 0.1)
        E[tid] += energy(eam, at1) / N
        F[tid] += forces(eam, at1) / N
    end
    return sum(E), sum(F)
end

at = bulk(:Ni, cubic=true) * 7

@info("Energy EAM0")
@btime energy($eam0, $at)

@info("Energy EAM1 (low acc)")
@btime energy($eam1, $at)

@info("Forces")
@btime forces($eam0, $at)
@btime forces($eam1, $at)

##

Profile.clear()
@profile energy(eam, at)

##

tmp = [JuLIP.alloc_temp(eam, at) for i in 1:JuLIP.nthreads()]
@time JuLIP.energy!(tmp, eam, at)

@code_warntype JuLIP.energy!(tmp, eam, at)

##

nlist = neighbourlist(at, cutoff(eam))
Js, Rs = JuLIP.neigs(nlist, 1)
Zs = at.Z[Js]
z0 = at.Z[1]

JuLIP.evaluate!(tmp[1], eam, Rs, Zs, z0)
@btime JuLIP.evaluate!($(tmp[1]), $eam, $Rs, $Zs, $z0)
@code_warntype JuLIP.evaluate!(tmp[1], eam, Rs, Zs, z0)

##

function runn(N, f!, args...)
    s = 0.0
    for n = 1:N
        s += 0.0001 * f!(args...)
    end
    return s
end

##

Profile.clear()
@profile runn(1000, JuLIP.evaluate!, tmp[1], eam, Rs, Zs, z0)

##

spl = eam.ϕ[1].f.spl
r = rnn(:Ni)
@btime $spl($r)

Profile.clear()

@profile runn(100_000, spl, r);

##

ProfileView.view()


text:
@testset "217.contains-duplicate.jl" begin
    @test contains_duplicate([1, 2, 3, 1]) == true
    @test contains_duplicate([1, 2, 3, 4]) == false
    @test contains_duplicate([1, 1, 1, 3, 3, 4, 3, 2, 4, 2]) == true
end


text:
list_test_str = "target area: x=20..30, y=-10..-5"

list_str = "target area: x=156..202, y=-110..-69"

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

function part1(list)
    c = parse(Int, match(r"y=(-?\d*)", list).captures[1])
    return div(c * (c + 1), 2)
end

testit("Test Part 1", part1(list_test_str), ==(45))
testit("Challenge Part 1", part1(list_str), ==(5995))

# Part 2

function part2(list)
    a, b, c, d = parse.(Int, match(r"(-?\d*)\.\.(-?\d*), y=(-?\d*)\.\.(-?\d*)", list).captures)
    good = Matrix{Int}(undef, 0, 3)
    for vx in 0:b
        for vy in c:-c
            xmax = div(vx * (vx + 1), 2)
            y = 0
            n = 0
            while y ≥ c
                n += 1
                # iterating vx, vy and then x and y is probably faster, though
                x = n < vx ? n * vx - div(n * (n - 1), 2) : xmax
                y = n * vy - div(n * (n - 1), 2)
                if x ∈ a:b && y ∈ c:d
                    good = vcat(good, [vx vy n])
                end
            end
        end
    end
    goodmaxima = [maximum(n * g[2] - div(n * (n - 1), 2) for n in 1:g[3]) for g in eachrow(good)]
    return good, size(unique(good[:, 1:2], dims=1), 1), maximum(goodmaxima)
end

testit("Test Part 2", part2(list_test_str)[2], ==(112))
testit("Challenge Part 2", part2(list_str)[2], ==(3202))


text:
### Fecth the require packages ###
using SafeTestsets

### Run the tests ###
@time begin

# Tests all features realted to constructing a model
@time @safetestset "Model Construction" begin include("make_model.jl") end
@time @safetestset "Custom Functions" begin include("custom_functions.jl") end
@time @safetestset "Model Modification" begin include("model_modification.jl") end

# Test api
@time @safetestset "API" begin include("api.jl") end

# Tests various core properties of the package.
@time @safetestset "Higher Order" begin include("higher_order_reactions.jl") end
@time @safetestset "U0 and Parameters Input Types" begin include("u0_n_parameter_inputs.jl") end

# Tests related to solving Ordinary Differential Equations.
@time @safetestset "ODE System Solving" begin include("solve_ODEs.jl") end
@time @safetestset "Make Jacobian" begin include("make_jacobian.jl") end
@time @safetestset "DiffEq Steady State Solving" begin include("steady_state_problems.jl") end

# Tests related to solving Stochastic Differential Equations.
@time @safetestset "SDE System Solving" begin include("solve_SDEs.jl") end

# Tests related to solvingJump Systems.
@time @safetestset "Jump System Solving" begin include("solve_jumps.jl") end

# Miscellaneous tests
#@time @safetestset "Basic Plotting" begin include("plotting.jl") end
@time @safetestset "Latexify" begin include("latexify.jl") end

end # @time


text:
# set to true to support intel fortran compiler
useIntelFortran = false

# construct absolute path to the dependency directory
depsdir  = splitdir(Base.source_path())[1]
builddir = joinpath(depsdir,"builds")
srcdir   = joinpath(depsdir,"src")

# print the directory
println("=== Building spmatvec ===")
println("depsdir  = $depsdir")
println("builddir = $builddir")
println("srcdir   = $srcdir")
println("useIntel = $useIntelFortran")

# create directory to save the binary
if !isdir(builddir)
	 println("creating build directory")
	 mkdir(builddir)
	 if !isdir(builddir)
		  error("Could not create build directory")
	 end
end

# fortran source code
path_src = joinpath(srcdir,"spmatvec.f90")

# the compiled object
path_bin = joinpath(builddir,"spmatvec.so")

# compile source code for linux machine
if Sys.islinux()
	 if useIntelFortran
			run(`ifort -O3 -xHost -fPIC -fpp -openmp -integer-size 64 -diag-disable=7841 -shared  $path_src -o $path_bin`)
	 else
			println("fortran version")
			run(`gfortran --version`)
			run(`gfortran -O3 -fPIC -cpp -fopenmp -fdefault-integer-8 -shared  $path_src -o $path_bin`)
	 end
end

# compile source code for Mac
if Sys.isapple()
	 if useIntelFortran
			run(`ifort -O3 -xHost -fPIC -fpp -openmp -integer-size 64 -diag-disable=7841 -shared  $path_src -o $path_bin`)
	 else
			println("fortran version")
			run(`gfortran --version`)
			run(`gfortran -O3 -fPIC -cpp -fopenmp -fdefault-integer-8 -shared  $path_src -o $path_bin`)
	 end
end

println("The fortran code is compiled successfully")


text:
###################################
### orthogonal matching pursuit ###
###################################

"A container object for all temporary arrays used in `omp`."
type OMPVariables{T <: Float}
    r    :: DenseVector{T}   # vector of residuals
    idxs :: DenseVector{Int} # vector of active indices
    a    :: DenseVector{T}   # temp vector x*b
    b    :: DenseVector{T}   # (sparse) model b
    dots :: DenseVector{T}   # dot products x'*r
end

"A funtion to create an `OMPVariables` object from a matrix `x`, a vector `y`, and a desired sparsity level `k`."
function OMPVariables{T <: Float}(
    x :: DenseMatrix{T},
    y :: DenseVector{T},
    k :: Int
)
    n,p  = size(x)
    r    = copy(y)
    idxs = ones(Int, k)
    a    = zeros(T, n)
    b    = zeros(T, p)
    dots = zeros(T, p)
    OMPVariables{T}(r, idxs, a, b, dots)
end

"""
    omp!(w::OMPVariables, x, y, k [, quiet::Bool = true])

This function computes **one step** of OMP; that is, it greedily computes the vector `b` of sparsity level `k` that minimizes `norm(y - x*b)`.
As a result, `omp!` **requires** computation of previous 1, 2, ..., k-1 indices in order to work correctly!
Thus, `omp!` is best used within a loop or a function such as `omp`.
"""
function omp!{T <: Float}(
    w :: OMPVariables{T},
    x :: DenseMatrix{T},
    y :: DenseVector{T},
    k :: Int;
)
    # dot_p = x' * r
    # need this for basis pursuit, since OMP selects dot product of largest magnitude
    At_mul_B!(w.dots, x, w.r)

    # compute index of next index to add to support
    λ = indmax_abs(w.dots)

    # expand active set indices Λ = [Λ i]
    w.idxs[k] = λ
    idx       = w.idxs[1:k]

    # get subset of x corresponding to active set
    temp = x[:, idx]
    #temp = view(x, :, idx) ### this variant burns CPU to avoid allocating memory; only good for biiiiig matrices

    # z = argmin_{b} ( norm(y - temp*b) )
    z, = lsqr(temp, y)

    # r = y - temp*z
    #BLAS.gemv!('N', -one(T), temp, z, zero(T), w.r)
    #BLAS.axpy!(one(T), y, w.r)
    A_mul_B!(w.r, temp, z)
    BLAS.axpy!(-one(T), y, w.r)

    # save current model to b
    # no need to erase w.b beforehand since OMP builds models stepwise,
    # and previous model (sparsity level k-1) is subset of current model of size k
    w.b[idx] = z

    return nothing
end

"""
    omp(x,y,k) -> SparseMatrixCSC

Perform *o*rthogonal *m*atching *p*ursuit using a matrix `x`, a response vector `y`, and a sparsity level `k`.
This function will compute in greedy fashion all sparsity levels `1, 2, ..., k`.
The result is a matrix of type `SparseMatrixCSC` containing the sparse vector that minimizes `sumsq(y - x*b)` for each model size.

Arguments:
- `x` is the design matrix (dictionary)
- `y` is the response vector (signal)
- `k` is the desired sparsity level

Output:
- `B` is the sparse matrix whose columns provide estimates of `y` with sparisty level `1, 2, ..., k`.
"""
function omp{T <: Float}(
    x :: DenseMatrix{T},
    y :: DenseVector{T},
    k :: Int;
    quiet::Bool = true
)
    # size of problem?
    n,p = size(x)

    # initialize all temporary arrays
    w = OMPVariables(x, y, k)

    # initialize sparse matrix of models
    # will fill this iteratively as we grow our beta
    # result is a "path" of betas
    B = spzeros(T, p, k)

    # compute models 1, 2, ..., k
    for i = 1:k

        # this computes model size i
        omp!(w, x, y, i)

        # output progress if desired
        quiet || @printf("Sparsity level: %d, sum(residuals): %f\n", i, norm(w.r));

        # save model for sparsity level i
        B[:,i] = sparsevec(w.b)
    end

    # return path
    return B
end


text:
# Copyright (c) 2014: SCS.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

struct IndirectSolver <: LinearSolver end

scsint_t(::Type{IndirectSolver}) = Clonglong

function scs_set_default_settings(
    ::Type{IndirectSolver},
    stgs::ScsSettings{I},
) where {I<:Clonglong}
    return @ccall(
        indirect.scs_set_default_settings(stgs::Ref{ScsSettings{I}})::Cvoid,
    )
end

function scs_init(
    ::Type{IndirectSolver},
    data::ScsData{I},
    cone::ScsCone{I},
    stgs::ScsSettings{I},
) where {I<:Clonglong}
    return @ccall indirect.scs_init(
        data::Ref{ScsData{I}},
        cone::Ref{ScsCone{I}},
        stgs::Ref{ScsSettings{I}},
    )::Ptr{Cvoid}
end

function scs_update(
    ::Type{IndirectSolver},
    work::Ptr{Cvoid},
    b::Vector{Float64},
    c::Vector{Float64},
)
    return @ccall direct.scs_update(
        work::Ptr{Cvoid},
        b::Ref{Float64},
        c::Ref{Float64},
    )::Clonglong
end

function scs_solve(
    ::Type{IndirectSolver},
    work::Ptr{Cvoid},
    solution::ScsSolution,
    info::ScsInfo{I},
    warm_start::Integer,
) where {I<:Clonglong}
    return @ccall indirect.scs_solve(
        work::Ptr{Cvoid},
        solution::Ref{ScsSolution},
        info::Ref{ScsInfo{I}},
        warm_start::Clonglong,
    )::Clonglong
end

function scs_finish(::Type{IndirectSolver}, work::Ptr{Cvoid})
    return @ccall indirect.scs_finish(work::Ptr{Cvoid})::Cvoid
end


text:
function inputlines(txtfile, iochannel)
    for line in readlines(txtfile)
        Base.put!(iochannel, line)
    end
    Base.put!(iochannel, nothing)
    println("The other task printed $(take!(iochannel)) lines.")
end

function outputlines(iochannel)
    totallines = 0
    while (line = Base.take!(iochannel)) != nothing
        totallines += 1
        println(line)
    end
    Base.put!(iochannel, totallines)
end

c = Channel(0)
@async inputlines("filename.txt", c)
outputlines(c)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test

o4 = UncertainValue(Uniform, rand(100))
o5 = UncertainValue(Beta, 4, 5)
o6 = UncertainValue(Gamma, 4, 5)
o7 = UncertainValue(rand(1000))
o8 = UncertainValue(Normal, 0, 2)

# 2-tuples
t1 = (o4, o5)

# n-tuples (a 5-tuple, including scalars)
t2 = (o4, o5, o6, o7, o8, 6);

@test resample(t1) isa NTuple{2, T} where T <: Number
@test resample(t2) isa NTuple{6, T} where T <: Number
@test resample(t1, 5) isa Vector{NTuple{2, T}} where {N, T <: Number}
@test resample(t2, 5) isa Vector{NTuple{6, T}} where {N, T <: Number}


text:
# ---
# title: 1048. Longest String Chain
# id: problem1048
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Dynamic Programming
# link: <https://leetcode.com/problems/longest-string-chain/description/>
# hidden: true
# ---
#
# Given a list of words, each word consists of English lowercase letters.
#
# Let's say `word1` is a predecessor of `word2` if and only if we can add
# exactly one letter anywhere in `word1` to make it equal to `word2`.  For
# example, `"abc"` is a predecessor of `"abac"`.
#
# A _word chain  _is a sequence of words `[word_1, word_2, ..., word_k]` with `k
# >= 1`, where `word_1` is a predecessor of `word_2`, `word_2` is a predecessor
# of `word_3`, and so on.
#
# Return the longest possible length of a word chain with words chosen from the
# given list of `words`.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["a","b","ba","bca","bda","bdca"]
#     Output: 4
#     **Explanation** : One of the longest word chain is "a","ba","bda","bdca".
#
#
# **Example 2:**
#
#
#
#     Input: words = ["xbc","pcxbcf","xb","cxbc","pcxbc"]
#     Output: 5
#
#
#
#
# **Constraints:**
#
#   * `1 <= words.length <= 1000`
#   * `1 <= words[i].length <= 16`
#   * `words[i]` only consists of English lowercase letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end