text:
# A `Table` presents itself as an `AbstractArray` of `NamedTuples`s

"""
    FlexTable(name1 = array1, ...)

Create a column-storage-based `FlexTable` with column names `name1`, etc, from arrays `array1`,
etc. The input arrays `array1`, etc, must share the same dimensionality and indices.

`FlexTable` itself is an `AbstractArray` whose elements are `NamedTuple`s of the form
`(name1 = first(array1), ...)`, etc. Rows of the table are obtained via standard array
indexing `table[i]`, and columns via `table.name`.

`FlexTable` differs from `Table` in that the columns are mutable - you may add, remove,
rename and replace entire columns of a `FlexTable`, but not a `Table`. However, `Table` can
access and iterate rows in local scope with fast, fully type-inferred code while `FlexTable`
will be more efficient with a higher-order interface.
"""
mutable struct FlexTable{N} <: AbstractArray{NamedTuple, N}
    data::NamedTuple{<:Any, <:Tuple{Vararg{AbstractArray{<:Any, N}}}}

    # Inner constructor, to compare axes?
    function FlexTable{N}(d::NamedTuple{<:Any, <:Tuple{Vararg{AbstractArray{<:Any, N}}}}) where N
        new(d)
    end
end

FlexTable(ts...; kwargs...) = _flextable(removenothings(merge(_columns(ts...), kwargs.data)))
FlexTable{N}(ts...; kwargs...) where {N} = _flextable(removenothings(merge(_columns(ts...), kwargs.data)))::FlexTable{N}

_flextable(nt::NamedTuple) = FlexTable{_ndims(nt)}(nt)

FlexTable(t::Table) = FlexTable(columns(t))
Table(t::FlexTable) = Table(columns(t))

_columns(t::FlexTable) = columns(t)

Tables.istable(::Type{<:FlexTable}) = true
Tables.rowaccess(::Type{<:FlexTable}) = true
Tables.columnaccess(::Type{<:FlexTable}) = true
Tables.schema(t::FlexTable) = Tables.Schema(_eltypes(columns(t)))
Tables.materializer(::FlexTable) = FlexTable

"""
    columns(dataframe::FlexTable)

Convert a `FlexTable` into a `NamedTuple` of its columns.
"""
@inline Tables.columns(t::FlexTable) = getfield(t, :data)

@inline Tables.rows(t::FlexTable) = Table(columns(t))

# Simple column access via `table.columnname`
@inline Base.getproperty(t::FlexTable, name::Symbol) = getproperty(columns(t), name)

function Base.setproperty!(t::FlexTable, name::Symbol, a)
    setfield!(t, :data, merge(columns(t), NamedTuple{(name,)}((convert(AbstractArray, a),))))
    return t
end

function Base.setproperty!(t::FlexTable, name::Symbol, ::Nothing)
    setfield!(t, :data, removenothings(merge(columns(t), NamedTuple{(name,)}((nothing,)))))
    return t
end

propertytype(::FlexTable{N}) where {N} = FlexTable{N}

"""
    columnnames(table)

Return a tuple of the column names of a `Table`.
"""
columnnames(t::FlexTable) = keys(columns(t))

# show
Base.show(io::IO, ::MIME"text/plain", t::FlexTable) = showtable(io, t)
Base.show(io::IO, t::FlexTable) = showtable(io, t)

# Basic AbstractArray interface

@inline Base.size(t::FlexTable{N}) where {N} = size(first(columns(t)))::NTuple{N, Integer}
@inline Base.axes(t::FlexTable{N}) where {N} = axes(first(columns(t)))::NTuple{N, Any}
@inline Base.IndexStyle(t::FlexTable) = IndexStyle(first(columns(t)))

function Base.checkbounds(::Type{Bool}, t::FlexTable, i...)
    # Make sure we are in bounds for *every* column. Only safe to do
    # here because each column might mutate size independently of the others!
    all(col -> checkbounds(Bool, col, i...), columns(t))
end

@inline function Base.getindex(t::FlexTable, i::Int)
    @boundscheck checkbounds(t, i)
    map(col -> @inbounds(getindex(col, i)), columns(t))::NamedTuple
end

@inline function Base.getindex(t::FlexTable, i::Int...)
    @boundscheck checkbounds(t, i...)
    map(col -> @inbounds(getindex(col, i...)), columns(t))::NamedTuple
end

@inline function Base.setindex!(t::FlexTable, v::NamedTuple, i::Int)
    @boundscheck begin
        checkbounds(t, i)
        @assert keys(v) === keys(columns(t))
    end
    map((val, col) -> @inbounds(setindex!(col, val, i)), v, columns(t))
    return t
end

@inline function Base.setindex!(t::FlexTable, v::NamedTuple, i::Int...)
    @boundscheck begin
        checkbounds(t, i)
        @assert keys(v) === keys(columns(t))
    end
    map((val, col) -> @inbounds(setindex!(col, val, i...)), v, columns(t))
    return t
end

# similar
@inline Base.similar(t::FlexTable{N}) where {N} = FlexTable{N}(similar(Table(t)))::FlexTable{N}
@inline Base.similar(t::FlexTable{N}, ::Type{NamedTuple}) where {N} = FlexTable{N}(similar(Table(t)))
@inline Base.similar(t::FlexTable{N}, ::Type{NamedTuple{names,T}}) where {N, names, T} = FlexTable{N}(similar(Table(t, NamedTuple{names, T})))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims) = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names,T}}, dims) where {names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))

# Ambiguities...
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims::Union{Integer, AbstractUnitRange}...) = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims::Tuple{Vararg{Int, N}}) where {N} = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple}, dims::Tuple{Union{Integer, OneTo},Vararg{Union{Integer, OneTo}}}) = FlexTable{__ndims(dims)}(similar(Table(t), _eltypes(columns(t)), dims))

@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names, T}}, dims::Union{Integer, AbstractUnitRange}...) where {names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names, T}}, dims::Tuple{Vararg{Int, N}}) where {N, names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))
@inline Base.similar(t::FlexTable, ::Type{NamedTuple{names, T}}, dims::Tuple{Union{Integer, OneTo},Vararg{Union{Integer, OneTo}}}) where {names, T} = FlexTable{__ndims(dims)}(similar(Table(t), NamedTuple{names, T}, dims))

@inline __ndims(::Integer) = 1
@inline __ndims(::AbstractUnitRange) = 1
@inline __ndims(::NTuple{N, Any}) where {N} = N

# empty
Base.empty(t::FlexTable) = empty(t, _eltypes(columns(t)))
function Base.empty(t::FlexTable, ::Type{NamedTuple{names, T}}) where {names, T <: Tuple}
    FlexTable(empty(Table(t)))
end

# Support Vector / deque interface (mutable-length vectors)

function Base.empty!(t::FlexTable)
    map(empty!, columns(t))
    return t
end

function Base.pop!(t::FlexTable)
    return map(pop!, columns(t))::NamedTuple
end

function Base.push!(t::FlexTable, v::NamedTuple)
    map(push!, columns(t), v)
    return t
end

function Base.append!(t::Union{FlexTable, Table}, t2::Union{FlexTable, Table})
    map(append!, columns(t), columns(t2))
    return t
end

function Base.popfirst!(t::FlexTable)
    return map(popfirst!, columns(t))::NamedTuple
end

function Base.pushfirst!(t::FlexTable, v::NamedTuple)
    map(pushfirst!, columns(t), v)
    return t
end

function Base.prepend!(t::Union{FlexTable, Table}, t2::Union{FlexTable, Table})
    map(prepend!, columns(t), columns(t2))
    return t
end

function Base.deleteat!(t::FlexTable, i)
    map(col -> deleteat!(col, i), columns(t))
    return t
end

function Base.insert!(t::FlexTable, i::Integer, v::NamedTuple)
    map((col, val) -> insert!(col, i, val), columns(t), v)
    return t
end

function Base.splice!(t::FlexTable, inds::Integer)
    return map(col -> splice!(col, inds), columns(t))::NamedTuple
end

function Base.splice!(t::FlexTable, inds::AbstractVector)
    return FlexTable{1}(map(col -> splice!(col, inds), columns(t)))
end

function Base.splice!(t::FlexTable, inds::Integer, ins::NamedTuple)
    return map((col, vals) -> splice!(col, inds, vals), columns(t), ins)::NamedTuple
end

function Base.splice!(t::FlexTable, inds::AbstractVector, ins::NamedTuple)
    return FlexTable{1}(map((col, vals) -> splice!(col, inds, vals), columns(t), ins))
end

function Base.splice!(t::FlexTable, inds::Integer, ins::AbstractVector{<:NamedTuple})
    return map((col, vals) -> splice!(col, inds, vals), columns(t), columns(ins))::NamedTuple
end

function Base.splice!(t::FlexTable, inds::AbstractVector, ins::AbstractVector{<:NamedTuple})
    return FlexTable{1}(map((col, vals) -> splice!(col, inds, vals), columns(t), columns(ins)))
end

# TODO splicing in an `AbstractArray{<:NamedTuple}` should be possible...

function Base.resize!(t::FlexTable, i::Integer)
    foreach(col -> resize!(col, i), columns(t))
    return t
end

# Speedups for column-based storage

function Base.getindex(t::FlexTable, inds::Union{AbstractArray, Colon}...)
    return FlexTable{_getindex_dims(inds)}(map(col -> getindex(col, inds...), columns(t)))
end

function Base.view(t::FlexTable, inds::Union{AbstractArray, Colon}...)
    return FlexTable{_getindex_dims(inds)}(map(col -> view(col, inds...), columns(t)))
end

@inline _getindex_dims(inds) = __getindex_dims(0, inds...)
@inline __getindex_dims(n::Int) = n
@inline __getindex_dims(n::Int, ::Int, inds...) = __getindex_dims(n, inds...)
@inline __getindex_dims(n::Int, ::AbstractArray{<:Any, m}, inds...) where {m} = __getindex_dims(n + m, inds...)
@inline __getindex_dims(n::Int, ::Colon, inds...) = __getindex_dims(n + 1, inds...)

# Deprecated for .= syntax (via Base.Broadcast.materialize!)
# It seems `Ref` might be the new cool here. Could also consider `AbstractArray{<:NamedTuple, 0}`?
#function Base.setindex!(t::Table{<:NamedTuple{names}}, v::NamedTuple{names}, inds::Union{AbstractArray, Colon}...) where {names}
#    map((col, val) -> setindex!(col, val, inds...), columns(t), v)
#    return t
#end

function Base.setindex!(t::FlexTable, t2::Union{FlexTable, Table}, inds::Union{AbstractArray, Colon}...)
    map((col, col2) -> setindex!(col, col2, inds...), columns(t), columns(t2))
    return t
end

# Private fields are never exposed since they can conflict with column names
Base.propertynames(t::FlexTable, private::Bool=false) = columnnames(t)

@inline function Base.vcat(ts::Union{FlexTable, Table}...)
    return FlexTable{_vcat_ndims(map(ndims, ts)...)}(map(vcat, map(columns, ts)...))
end

@inline function Base.hcat(ts::Union{FlexTable, Table}...)
    return FlexTable{_hcat_ndims(map(ndims, ts)...)}(map(hcat, map(columns, ts)...))
end

function Base.hvcat(rows::Tuple{Vararg{Int}}, ts::Union{FlexTable, Table}...)
    return FlexTable(map((cols...,) -> hvcat(rows, cols...), map(columns, ts)...))
end

@pure function _vcat_ndims(i::Int...)
    max(1, i...)
end

@pure function _hcat_ndims(i::Int...)
    max(2, i...)
end

function Base.vec(t::FlexTable)
    return FlexTable{1}(map(vec, columns(t)))
end

# "Bulk" operations on FlexTables should generally first unrwap to Tables
_flex(t::Table{<:Any, N}) where {N} = FlexTable(columns(t))
_flex(t) = t

Broadcast.broadcastable(t::FlexTable) = Table(t)

Base.map(f, t::FlexTable{N}) where {N} = _flex(map(f, rows(t)))::AbstractArray{<:Any, N}
Base.map(f, t::FlexTable{N}, t2) where {N} = _flex(map(f, rows(t), t2))::AbstractArray{<:Any, N}
Base.map(f, t, t2::FlexTable{N}) where {N} = _flex(map(f, t, rows(t2)))::AbstractArray{<:Any, N}
Base.map(f, t::FlexTable{N}, t2::FlexTable{N}) where {N} = _flex(map(f, rows(t), rows(t2)))::AbstractArray{<:Any, N}

Base.mapreduce(f, op, t::FlexTable; kwargs...) = mapreduce(f, op, rows(t); kwargs...)

Base.filter(f, t::FlexTable{N}) where {N} = FlexTable(filter(f, rows(t)))::FlexTable{N}

SplitApplyCombine.mapview(f::Union{Function, Type}, t::FlexTable{N}) where {N} = _flex(mapview(f, rows(t)))::AbstractArray{<:Any, N}
SplitApplyCombine.mapview(f::Union{Function, Type}, t::FlexTable{N}, t2) where {N} = _flex(mapview(f, rows(t), t2))::AbstractArray{<:Any, N}
SplitApplyCombine.mapview(f::Union{Function, Type}, t, t2::FlexTable{N}) where {N} = _flex(mapview(f, t, rows(t2)))::AbstractArray{<:Any, N}
SplitApplyCombine.mapview(f::Union{Function, Type}, t::FlexTable{N}, t2::FlexTable{N}) where {N} = _flex(mapview(f, rows(t), rows(t2)))::AbstractArray{<:Any, N}

SplitApplyCombine.group(by, f, t::FlexTable) = group(by, f, rows(t))
SplitApplyCombine.groupview(by, f, t::FlexTable) = groupview(by, f, rows(t))
SplitApplyCombine.groupinds(by, t::FlexTable) = groupinds(by, rows(t))
SplitApplyCombine.groupreduce(by, f, op, t::FlexTable; kwargs...) = groupreduce(by, f, op, rows(t); kwargs...)

SplitApplyCombine.innerjoin(lkey::Union{Function, Type}, rkey::Union{Function, Type}, f::Union{Function, Type}, cmp::Union{Function, Type}, t1::FlexTable, t2) = _flex(innerjoin(lkey, rkey, f, cmp, rows(t1), t2))
SplitApplyCombine.innerjoin(lkey::Union{Function, Type}, rkey::Union{Function, Type}, f::Union{Function, Type}, cmp::Union{Function, Type}, t1, t2::FlexTable) = _flex(innerjoin(lkey, rkey, f, cmp, t1, rows(t2)))
SplitApplyCombine.innerjoin(lkey::Union{Function, Type}, rkey::Union{Function, Type}, f::Union{Function, Type}, cmp::Union{Function, Type}, t1::FlexTable, t2::FlexTable) = _flex(innerjoin(lkey, rkey, f, cmp, rows(t1), rows(t2)))

Base.:(==)(t1::FlexTable{N}, t2::AbstractArray{<:Any,N}) where {N} = (rows(t1) == t2)
Base.:(==)(t1::AbstractArray{<:Any,N}, t2::FlexTable{N}) where {N} = (t1 == rows(t2))
Base.:(==)(t1::FlexTable{N}, t2::FlexTable{N}) where {N} = (rows(t1) == rows(t2))

Base.isequal(t1::FlexTable{N}, t2::AbstractArray{<:Any,N}) where {N} = isequal(rows(t1), t2)
Base.isequal(t1::AbstractArray{<:Any,N}, t2::FlexTable{N}) where {N} = isequal(t1, rows(t2))
Base.isequal(t1::FlexTable{N}, t2::FlexTable{N}) where {N} = isequal(rows(t1), rows(t2))

Base.isless(t1::FlexTable{1}, t2::AbstractVector) = isless(rows(t1), t2)
Base.isless(t1::AbstractVector, t2::FlexTable{1}) = isless(t1, rows(t2))
Base.isless(t1::FlexTable{1}, t2::FlexTable{1}) = isless(rows(t1), rows(t2))

Base.hash(t::FlexTable, h::UInt) = hash(rows(t), h)

Adapt.adapt_structure(to, t::FlexTable{N}) where N = FlexTable{N}(; Adapt.adapt(to, getfield(t, :data))...)


text:
module GridPortionsTests

using Test
using Gridap

oldgrid = CartesianGrid(partition=(3,3))

n = ncells(oldgrid)
nini = ceil(n/3)
newcell_to_oldcell = collect(Int,nini:n)
grid = GridPortion(oldgrid,newcell_to_oldcell)

test_grid(grid,16,7)

#writevtk(grid,"grid")

model = CartesianDiscreteModel(partition=(3,4,2))

labels = labels_on_dim(model,2)
oldgrid = Grid(model,2)

mask = [ (label == 23 || label == 24 || label == 25) for label in labels ]

newcell_to_oldcell = findall(mask)
grid = GridPortion(oldgrid,newcell_to_oldcell)

test_grid(grid,60,20)

#writevtk(grid,"grid",celldata=["oldcells"=>newcell_to_oldcell])

end # module

ext: jl


text:
immutable SkierObsDist{G<:Function}
    p::PowseekerPOMDP{G}
    gps::Bool
    s::SkierState
    sp::SkierState
end

function rand(rng::AbstractRNG, d::SkierObsDist)
    if d.gps
        return SkierObs(d.sp.time, d.sp.pos + d.p.gps_std*randn(rng, 2), d.sp.psi + d.p.compass_std*randn(rng))
    else
        true_dist = norm(d.sp.pos-d.s.pos)
        dist = true_dist*(1.0 + d.p.dist_std_frac*randn(rng))
        grad_local = true_grad_local(d)
        grad_obs = grad_local + d.p.grad_std*randn(rng, 2)
        return SkierObs(d.sp.time, dist, grad_obs)
    end
end

function true_grad_local(d::SkierObsDist)
    grad_global = mdp(d.p).gradient(d.sp.pos)
    psi = d.sp.psi
    C = @SMatrix [cos(psi) sin(psi); -sin(psi) cos(psi)]
    grad_local = C*grad_global
end

function pdf(d::SkierObsDist, o::SkierObs)
    if d.gps
        meas = get(o.pos) # should error if something is inconsistent
        meas_pdf = exp(-sum((meas - d.sp.pos).^2)/(2*d.p.gps_std^2))
        dir = get(o.dir)
        dir_diff = abs(dir - d.sp.psi)
        while dir_diff > pi
            dir_diff -= 2*pi
        end
        dir_diff = abs(dir_diff)
        dir_pdf = exp(-dir_diff^2/(2*d.p.compass_std^2))
        return meas_pdf*dir_pdf
    else
        true_dist = norm(d.sp.pos-d.s.pos)
        dist_pdf = exp(-(o.dist-true_dist)^2/(2*(true_dist*d.p.dist_std_frac)^2))
        true_grad = true_grad_local(d)
        grad_pdf = exp(-sum((o.grad-true_grad).^2)/(2*d.p.grad_std^2))
        return dist_pdf*grad_pdf
    end
end


text:
function mdl_freq_estimator_ode!(
    device_states,
    output_ode,
    f0,
    ω_sys,
    dynamic_device::PSY.DynamicInverter{C, O, IC, DC, PSY.KauraPLL, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    F <: PSY.Filter,
}

    #Obtain external states inputs for component
    external_ix = get_input_port_ix(dynamic_device, PSY.KauraPLL)
    Vr_filter = device_states[external_ix[1]]
    Vi_filter = device_states[external_ix[2]]

    #V_tR = get_inner_vars(dynamic_device)[VR_inv_var]
    #V_tI = get_inner_vars(dynamic_device)[VI_inv_var]

    #Get parameters
    pll_control = PSY.get_freq_estimator(dynamic_device)
    ω_lp = PSY.get_ω_lp(pll_control)
    kp_pll = PSY.get_kp_pll(pll_control)
    ki_pll = PSY.get_ki_pll(pll_control)
    ωb = 2.0 * pi * f0

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

    #Define internal states for frequency estimator
    internal_states = @view device_states[local_ix]
    vpll_d = internal_states[1]
    vpll_q = internal_states[2]
    ϵ_pll = internal_states[3]
    θ_pll = internal_states[4]

    #Transform to internal dq-PLL reference frame
    V_dq_pll = ri_dq(θ_pll + pi / 2) * [Vr_filter; Vi_filter]

    #Inputs (control signals)

    #Compute 6 states ODEs (D'Arco EPSR122 Model)
    #Output Voltage LPF (internal state)
    #𝜕vpll_d/𝜕t, D'Arco ESPR122 eqn. 12
    output_ode[local_ix[1]] = ω_lp * (V_dq_pll[d] - vpll_d)
    #𝜕vpll_q/𝜕t, D'Arco ESPR122 eqn. 12
    output_ode[local_ix[2]] = ω_lp * (V_dq_pll[q] - vpll_q)
    #PI Integrator (internal state)
    #𝜕dϵ_pll/𝜕t, D'Arco ESPR122 eqn. 13
    output_ode[local_ix[3]] = atan(vpll_q / vpll_d)
    #PLL Frequency Deviation (internal state)
    #𝜕θ_pll/𝜕t, D'Arco ESPR122 eqn. 15
    output_ode[local_ix[4]] = (ωb * kp_pll * atan(vpll_q / vpll_d) + ωb * ki_pll * ϵ_pll)

    #Update inner_vars
    #PLL frequency, D'Arco EPSR122 eqn. 16
    get_inner_vars(dynamic_device)[ω_freq_estimator_var] =
        (kp_pll * atan(vpll_q / vpll_d) + ki_pll * ϵ_pll + ω_sys)
end

function mdl_freq_estimator_ode!(
    device_states,
    output_ode,
    f0,
    ω_sys,
    dynamic_device::PSY.DynamicInverter{C, O, IC, DC, PSY.FixedFrequency, F},
) where {
    C <: PSY.Converter,
    O <: PSY.OuterControl,
    IC <: PSY.InnerControl,
    DC <: PSY.DCSource,
    F <: PSY.Filter,
}

    #Get parameters
    pll_control = PSY.get_freq_estimator(dynamic_device)
    frequency = PSY.get_frequency(pll_control)

    #Update inner_vars
    #PLL frequency
    get_inner_vars(dynamic_device)[ω_freq_estimator_var] = frequency
end


text:
include("../Utils/hist2D.jl")
function OT_sampler(plot_flag,savedata_flag)
    #### Test Examples
    ## Example 1: Beta distribution
    alp = 2
    bet = 5
    B = gamma(alp)*gamma(bet)/gamma(alp+bet)
    f1(x) = x^(alp-1)*(1-x)^(bet-1)/B
    range1 = [0 1]
    M1 = 100
    ## Example 2: Normal Distribution
    mu = 0
    sig2 = 1
    f2(x) = (1/sqrt(2*pi*sig2))*exp((-1/2)*((x-mu)^2)/sig2)
    range2 = [-6 6]
    M2 = 200
    ## Example 3: Cauchy Distribution
    gam = 0.5
    x_0 = 0
    f3(x) = (1/(pi*gam))*(gam^2/((x-x_0)^2+gam^2))
    range3 = [-4 4]
    M3 = 100
    ## Example 4: Bimodal Gaussian distribution
    x_01 = -10
    sigma21 = 2
    x_02 = 10
    sigma22 = 2
    p1 = 0.1
    f4(x) = p1*(1/sqrt(2*pi*sigma21))*exp((-1/2)*((x-x_01)^2)/sigma21) + (1-p1)*(1/sqrt(2*pi*sigma22))*exp((-1/2)*((x-x_02)^2)/sigma22)
    range4 = [-20 20]
    M4 = 200
    ## Example 5: Von Mises Distribution
    k = 8
    I0_k = besseli(0,k)
    f5(x) = exp(k*cos(x))/(2*pi*I0_k)
    range5 = [-1*pi pi]
    M5 = 100
    ##-------------------------------------------------------------------------
    f = Function[f1,f2,f3,f4,f5]
    Mlist = [M1,M2,M3,M4,M5]
    range = [range1,range2,range3,range4,range5]
    for i = 1:length(range)
        # First generate uniform samples from [a,b]
        rangei = range[i]
        a = rangei[1]
        b = rangei[2]
        M = Mlist[i]
        uni_samples = [a+(b-a)*i0/(M+1) for i0 in 1:M] # inline for function
        fi = f[i]
        weight_gen = map(fi,uni_samples)
        w = weight_gen/sum(weight_gen) # weights normalized
        D = abs.(repeat(uni_samples,1,M)-repeat(uni_samples',M,1))
        D_vec = reshape(D,M*M,1) # read along the column
        A1 = zeros(M,M*M)
        A2 = zeros(M,M*M)
        for i1 = 1:M
            A1[i1,M*(i1-1)+1:M*i1] = ones(1,M)
            for j1 = 1:M
                A2[1+j1-1,M*(i1-1)+1+j1-1] = 1
            end
        end
        b2 = ones(M)
        myModel = Model(with_optimizer(GLPK.Optimizer)) #MODEL CONSTRUCTION
        @variable(myModel, p[i2=1:(M*M)] >= 0) # Models x >=0
        @objective(myModel, Min, sum(D_vec[j2]*p[j2] for j2=1:M*M)) # Sets the objective to be minimized. For maximization use Max
        for i3=1:M # for all rows do the following
            @constraint(myModel, sum(A1[i3,j3]*p[j3] for j3=1:M*M) == w[i3]*M)
            @constraint(myModel, sum(A2[i3,j3]*p[j3] for j3=1:M*M) == b2[i3])
        end
        status = @time optimize!(myModel) # solves the model
        p_mat = reshape(value.(p),M,M)
        obj_val = getobjectivevalue(myModel)
        ##
        x_a = zeros(M)
        for i4= 1:M
            x_a[i4] = sum(p_mat[i4,j4]*uni_samples[j4] for j4=1:M)
        end
        if plot_flag == 1
            plotly()
            Plots.histogram(x_a, normed=true,
                legend=:none,
                titlefont=font(9),
                color=:red,
                bins=60)
            Plots.plot!(fi,a,b,color = "blue",show = true)#,reuse = false,show = true) # PGF Plotting
        end
        if savedata_flag == 1
            # Draw a 2D histogram
            (bins , xpoints) = hist2D(x_a,80)
            xmidpts = (xpoints[1:(end-1),1]+xpoints[2:end,1])./2
            data2d  = [xmidpts'; bins']'
            binlength = (xpoints[2]-xpoints[1])
            binnorm = binlength*sum(data2d[:,2])
            data2d[:,2] = data2d[:,2]./binnorm
            str1 = @sprintf("../data/out%i.dat",Int(i))
            writedlm(str1, data2d, " ")
        end
    end
    return 1
end

max_stars_repo_name: niladridas/UKF

max_issues_repo_name: niladridas/UKF

max_forks_repo_name: niladridas/UKF


text:
#=
- An example of using Luxor.jl
- https://cormullion.github.io/pub/2019-06-23-slackmojif.html
=#

using Luxor

function frame(scene, framenumber)
    background("white")
    juliacircles(framenumber)
    sethue("black")
    fontsize(50)
    text(string(framenumber), halign = :center)
end

simpleanimation = Movie(400, 400, "animation")

destination  = "simple.gif"

animate(
    simpleanimation,
    [Scene(simpleanimation, frame)],
    creategif = true,
    pathname = destination,
)

include("viewer.jl")

Viewer(destination)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

using BoundingSphere
const MB = BoundingSphere
using StaticArrays

if VERSION < v"0.7-"
    using Base.Test
    const stdout = STDOUT
    seed!(x) = srand(x)
else
    using Random
    using Random: seed!
    using LinearAlgebra
    using Test
    using InteractiveUtils: subtypes
end

include("helpers.jl")
include("test_broken.jl")
include("test_util.jl")
include("test_welzl.jl")
include("test_degenerate_examples.jl")
include("perf.jl")


text:
using ExactWrightFisher

@testset "test adaptive precomputation approx filtering functions" begin

    function simulate_WF3_data()
        K = 4
        α = ones(K)
        sα = sum(α)
        Pop_size_WF3 = 15
        Ntimes_WF3 = 3
        time_step_WF3 = 0.1
        time_grid_WF3 = [k*time_step_WF3 for k in 0:(Ntimes_WF3-1)]
        Random.seed!(4)
        wfchain_WF3 = Wright_Fisher_K_dim_exact_trajectory(rand(Dirichlet(K,0.3)), time_grid_WF3[1:(end-1)], α)
        wfobs_WF3 = [rand(Multinomial(Pop_size_WF3, wfchain_WF3[:,k])) for k in 1:size(wfchain_WF3,2)] |> l -> hcat(l...)
        data_WF3 = Dict(zip(time_grid_WF3 , [wfobs_WF3[:,t] for t in 1:size(wfobs_WF3,2)]))
        return data_WF3, wfchain_WF3, α
    end
    data, wfchain_WF3, α = simulate_WF3_data()

    @test_nowarn DualOptimalFiltering.filter_WF_adaptive_precomputation_keep_fixed_number(α, data |> DualOptimalFiltering.prepare_WF_dat_1D_2D |> last, 5; silence = false)

    @test_nowarn DualOptimalFiltering.filter_WF_adaptive_precomputation_keep_above_threshold(α, data |> DualOptimalFiltering.prepare_WF_dat_1D_2D |> last, 0.01; silence = false)

    @test_nowarn DualOptimalFiltering.filter_WF_adaptive_precomputation_keep_fixed_fraction(α, data |> DualOptimalFiltering.prepare_WF_dat_1D_2D |> last, 0.95; silence = false)
end;

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
import FileIO
import Observables
import AbstractPlotting

## getsurface and buildsurface

scene = AbstractPlotting.Scene()
transparency = false
segments = 36
radius = 1 + rand()
value1 = constructsphere(q, radius, segments = segments)
color = fill(AbstractPlotting.RGBAf0(rand(4)...), segments, segments)
observable = buildsurface(scene, value1, color, transparency = transparency)
value2 = getsurface(observable, segments)

@test typeof(observable) == Tuple{Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}}}
@test typeof(value2) == Array{ℝ³,2}
@test isapprox(value1, value2, atol = TOLERANCE)

## getsurface with different segments in the coordinate chart axes

number = rand(5:10)
points = [Geographic(rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
top = U1(rand() * 2pi - pi)
bottom = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
value1 = constructwhirl(points,
                        top = top,
                        bottom = bottom,
                        s3rotation = s3rotation,
                        config = config,
                        segments = segments)

color = fill(AbstractPlotting.RGBAf0(rand(4)...), segments, number)
transparency = rand(1:2) == 1 ? true : false
observable = buildsurface(scene, value1, color, transparency = transparency)
value2 = getsurface(observable, segments, number)

@test typeof(observable) == Tuple{Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}}}
@test typeof(value2) == Array{ℝ³,2}
@test isapprox(value1, value2, atol = TOLERANCE)

## builsurface with "observable" color array argument

observablecolor = Observables.Observable(fill(AbstractPlotting.RGBAf0(rand(4)...),
                                              segments,
                                              segments))
segments = rand(5:10)
radius = 1 + rand()
value1 = constructsphere(q, radius, segments = segments)
observable = buildsurface(scene, value1, observablecolor, transparency = transparency)
value2 = getsurface(observable, segments)

@test typeof(observable) == Tuple{Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}}}
@test typeof(value2) == Array{ℝ³,2}
@test isapprox(value1, value2, atol = TOLERANCE)

## builsurface with image as color

color = FileIO.load("data/basemap90grid.png")
segments = rand(5:10)
radius = 1 + rand()
value1 = constructsphere(q, radius, segments = segments)
observable = buildsurface(scene, value1, color, transparency = transparency)
value2 = getsurface(observable, segments)

@test typeof(observable) == Tuple{Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}},
                                  Observables.Observable{Array{Float64,2}}}
@test typeof(value2) == Array{ℝ³,2}
@test isapprox(value1, value2, atol = TOLERANCE)

## updatesurfce

value2 = map(x -> 2x, value1)
updatesurface(value2, observable)
value3 = getsurface(observable, segments)

@test isapprox(value1, value2, atol = TOLERANCE) == false
@test isapprox(value2, value3, atol = TOLERANCE)


text:
using Test
using Statistics
using Random
using TrackerFlux
using Tracker
using Flux
using Pkg
using UUIDs

Random.seed!(1234)

@testset "TrackerFlux.jl" begin

x = randn(Float32, 10, 1, 100)
y = mean(x, dims = 1)

θ0, re = Flux.destructure(Dense(2, 1))
@test Tracker.gradient(θ0) do θ
    sum(re(θ)([0.5, 2.0]))
end[1].data == [0.5, 2.0, 1]

model = Chain(LSTM(10, 100), LSTM(100, 1)) |> TrackerFlux.track
θ0, re = Flux.destructure(model)
@test Flux.destructure(re(θ0))[1] == θ0

model(x[:, :, 1])
@test Tracker.istracked(model[1].state[1])
Flux.truncate!(model)
ver = Pkg.dependencies()[UUID("587475ba-b771-5e3f-ad9e-33799f191a9c")].version
ver < v"0.12" && @test !Tracker.istracked(model[1].state[1])
@test !Tracker.istracked(TrackerFlux.untrack(model)[1].cell.Wh)

function loss(x, y)
    xs = Flux.unstack(x, 3)
    ys = Flux.unstack(y, 3)
    ŷs = model.(xs)
    l = 0f0
    for t in 1:length(ŷs)
        l += Flux.mse(ys[t], ŷs[t])
    end
    return l / length(ŷs)
end
ps = Flux.params(model)
data = repeat([(x, y)], 100)
opt = ADAMW(1e-3, (0.9, 0.999), 1e-4)
cb = () -> Flux.reset!(model)
TrackerFlux.overload_gradient()
Flux.@epochs 10 Flux.train!(loss, ps, data, opt, cb = cb)
@test loss(x, y) < 0.02

end


text:
using RomanNumerals
using Compat, Compat.Test
using Primes

# Constructor tests
@test RomanNumeral(46) == RomanNumeral("XLVI")
@test RomanNumeral(1999) == RomanNumeral("MCMXCIX")
@test rn"XXXX" == rn"XL"
@test RomanNumeral(1) == rn"I"
@test rn"xx" == rn"XX"
@test parse(RomanNumeral, "xx") == rn"XX"

# Invalid numeral tests
@test_throws Meta.ParseError RomanNumeral("nope")
@test_throws Meta.ParseError RomanNumeral("XLX")

# arithmetic tests
@test rn"I" + rn"IX" == rn"X"
@test rn"L" - RomanNumeral(1) == rn"XLIX"
@test rn"X" > rn"IX"
@test isqrt(rn"MXXIV") == rn"XXXII"
@test rn"XX" ^ rn"II" == rn"CD"

# primus numeri et factorii (quis scit?)
@test factor(rn"XX")[rn"V"]    == rn"I"
@test length(primes(rn"XXIX")) == rn"X"


text:
__precompile__(true)
module JuliaOpModule
    using ADCME
    using PyCall
    global do_it_op
    function __init__()
        global do_it_op = load_op("$(@__DIR__)/DoItOp/build/libDoItOp", "do_it_op")
    end

    export DoIt!, do_it
    function DoIt!(x::Array{Float64}, y::Array{Float64})
        x[:] = y
    end
    function do_it(o::PyObject)
        do_it_op(o)
    end
end


text:
module Ranges

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template, register_normal_element

using Stipple

export range, RangeData, slider

register_normal_element("q__range", context = @__MODULE__)
register_normal_element("q__slider", context = @__MODULE__)

Base.@kwdef mutable struct RangeData{T}
  range::UnitRange{T}
end

function Base.range(
                range::AbstractRange{T} where T <: Real,
                fieldname::Union{Symbol,Nothing} = nothing,
                args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                lazy = false,
                kwargs...)

  wrap() do
    q__range( "", lazy ? @on(:change," val => { $fieldname } ") : "" , args...;
              attributes(
                        [ Symbol(":min") => first(range),
                          Symbol(":max") => last(range),
                          Symbol(":step") => step(range),
                          ( lazy ? () : (:fieldname => fieldname,) )...,
                          kwargs...
                        ], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

function slider(range::AbstractRange{T} where T <: Real,
                fieldname::Union{Symbol, Nothing} = nothing,
                args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                lazy = false,
                kwargs...)

  wrap() do
    q__slider("", lazy ? @on(:change," val => { $fieldname } ") : "", args...;
    attributes(
      [Symbol(":min") => first(range), Symbol(":max") => last(range), Symbol(":step") => step(range),
      ( lazy ? () : (:fieldname => fieldname,) )..., kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

function Stipple.render(rd::RangeData{T}, fieldname::Union{Symbol,Nothing} = nothing) where {T,R}
  Dict(:min => rd.range.start, :max => rd.range.stop)
end

function Base.parse(::Type{RangeData{T}}, d::Dict{X,Y}) where {T,X,Y}
  RangeData(d["min"]:d["max"])
end

end


text:
using KeywordCalls
using Documenter

DocMeta.setdocmeta!(KeywordCalls, :DocTestSetup, :(using KeywordCalls); recursive=true)

makedocs(;
    modules=[KeywordCalls],
    authors="Chad Scherrer <chad.scherrer@gmail.com> and contributors",
    repo="https://github.com/cscherrer/KeywordCalls.jl/blob/{commit}{path}#{line}",
    sitename="KeywordCalls.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://cscherrer.github.io/KeywordCalls.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/cscherrer/KeywordCalls.jl",
)


text:
using COSMO, Test, LinearAlgebra, SparseArrays, Random, QDLDL, Pkg
include("COSMOTestUtils.jl")

# check if optional dependencies are available
test_iterative_solvers = pkg_installed("IterativeSolvers", "42fd0dbc-a981-5370-80f2-aaf504508153") && pkg_installed("LinearMaps", "7a12625a-238d-50fd-b39a-03d52299707e")
test_pardiso = pkg_installed("Pardiso", "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2")
test_iterative_solvers = false
test_pardiso = false

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

function make_test_kkt(P::AbstractMatrix{T}, A::AbstractMatrix{T}, sigma, rho) where {T <: AbstractFloat}

    m, n = size(A)
    R = length(rho)   == 1 ? diagm(0 => ones(T, m) ./ rho[1]) : Diagonal( one(T) ./ rho)
    S = length(sigma) == 1 ? (sigma[1]) * I : Diagonal(sigma)

    #compute the full KKT matrix
    K    = [P+S A'; A -R]
    return K
end

function add_kwargs(array; kwargs...)
    push!(array, kwargs)
end

@testset "KKT solver" begin

    @testset "KKT solver (T = $(T))" for T in UnitTestFloats
        rng = Random.MersenneTwister(1)

        solver_types = []
        solver_tols = []
        params = []
        push!(solver_types, COSMO.QdldlKKTSolver)
        T <: Float32 ? push!(solver_tols, T(1e-5)) : push!(solver_tols, 1e-10)
        add_kwargs(params)
        T == Float64 && push!(solver_types, COSMO.CholmodKKTSolver) #SuiteSparse Cholmod only supports doubles
        T <: Float32 ? push!(solver_tols, T(1e-5)) : push!(solver_tols, 1e-10)
        add_kwargs(params)

        # optional dependencies
        if test_pardiso && T <: Float64
            using Pardiso
            if Pardiso.LOCAL_MKL_FOUND
                push!(solver_types, COSMO.MKLPardisoKKTSolver)
                push!(solver_tols, 1e-10)
                    add_kwargs(params)
                end
            if Pardiso.PARDISO_LIB_FOUND
                push!(solver_types, COSMO.PardisoDirectKKTSolver)
                push!(solver_tols, 1e-10)
                add_kwargs(params)
            end
            if Pardiso.PARDISO_LIB_FOUND
                push!(solver_types, COSMO.PardisoIndirectKKTSolver)
                push!(solver_tols, 5e-5)
                add_kwargs(params)
            end
        end

        if test_iterative_solvers && T <: Float64
            using IterativeSolvers, LinearMaps
            push!(solver_types, COSMO.IndirectReducedKKTSolver)
            push!(solver_tols, 1e-3)
            add_kwargs(params, solver_type=:CG)

            push!(solver_types, COSMO.IndirectReducedKKTSolver)
            push!(solver_tols, 1e-3)
            add_kwargs(params, solver_type=:MINRES)

            push!(solver_types, COSMO.IndirectKKTSolver)
            push!(solver_tols, 1e-3)
            add_kwargs(params, solver_type=:MINRES)
        end

            @testset "$(solver_types[i]): KKT solver tests"  for i = 1:length(solver_types)
                m = 10
                n = 20

                for rho1 in [rand(T), rand(T, m)],
                    rho2 in [rand(T), rand(T, m)],
                    sigma in [rand(T)]

                    Pf = generate_pos_def_matrix(rng, n, MT = T)
                    P  = sparse(Pf)
                    A  = sprand(T, m, n, 0.2)
                    b = rand(T, m + n)

                    F = solver_types[i](P, A, sigma, rho1; params[i]...)
                    if test_iterative_solvers
                        if isa(F, COSMO.IndirectReducedKKTSolver) || isa(F, COSMO.IndirectKKTSolver)
                            # Technically, we should have been able to set the tolerance even lower
                            # But, currently, issues in IterativeSolvers.jl do not allow this
                            # See: https://github.com/JuliaMath/IterativeSolvers.jl/pull/244
                            F.iteration_counter = 10^4 # This forces CG's/MINRES tolerance to be 1/10^4
                        end
                    end
                    J = make_test_kkt(P, A, sigma, rho1)
                    x = copy(b)

                    #test a single solve
                    COSMO.solve!(F, x, b)
                    @test norm(x - Matrix(J) \ b) <= solver_tols[i]

                    # Check that warm starting works
                    # Invoking again an indirect solver should result in the solution with only
                    # one matrix vector multiplication
                    if test_iterative_solvers
                        if isa(F, COSMO.IndirectReducedKKTSolver) || isa(F, COSMO.IndirectKKTSolver)
                            # The calculation of the residual, and thus the termination criterion, of
                            # MINRES is approximate. Thus warm started solutions won't necessarily finish in one step
                            # For this reason we don't check warm starting with MINRES for now :(
                            if F.solver_type != :MINRES
                                COSMO.solve!(F, x, b)
                                # @test F.multiplications[end] <= 1 uncommented as test causes non-deterministic behaviour
                                @test norm(x - J \ b) <= solver_tols[i]
                            end
                        end
                    end

                    #test a rho update and solve
                    J = make_test_kkt(P, A, sigma, rho2)
                    x = copy(b)
                    COSMO.update_rho!(F, rho2)
                    COSMO.solve!(F, x, b)
                    @test norm(x - Matrix(J) \ b) <= solver_tols[i]

                 end

            end

        if T == Float64
            # Try every solver on the same problem and compare results
            solver_types = Array{Union{Type{<: COSMO.AbstractKKTSolver}, COSMO.OptionsFactory{<: COSMO.AbstractKKTSolver}}}(undef, 0)
            solver_names = Array{String}(undef, 0)
            push!(solver_types, COSMO.QdldlKKTSolver, COSMO.CholmodKKTSolver)
            push!(solver_names, "QDLDL", "CHOLMOD")
            if test_pardiso
                Pardiso.LOCAL_MKL_FOUND && push!(solver_types, COSMO.MKLPardisoKKTSolver)
                if Pardiso.PARDISO_LIB_FOUND
                    push!(solver_types, COSMO.PardisoDirectKKTSolver)
                    push!(solver_types, COSMO.PardisoIndirectKKTSolver)
                    push!(solver_names, "PardisoDirectKKTSolver", "PardisoIndirectKKTSolver")

                end
            end
            if test_iterative_solvers
                push!(solver_types, COSMO.CGIndirectKKTSolver)
                push!(solver_types, COSMO.MINRESIndirectKKTSolver)
                push!(solver_names, "CGIndirectKKTSolver", "MINRESIndirectKKTSolver")
                push!(solver_types, with_options(COSMO.CGIndirectKKTSolver))
                push!(solver_types, with_options(COSMO.MINRESIndirectKKTSolver))
                push!(solver_names, "CGIndirectKKTSolver (options)", "MINRESIndirectKKTSolver (options)")

            end

             P = T[4. 1;1 2]
             q = T[1; 1.]
             A = T[1. 1;1 0; 0 1]
             l = T[1.;0;0]
             u = T[1.;0.7;0.7]

             constraint1 = COSMO.Constraint(-A, u, COSMO.Nonnegatives)
             constraint2 = COSMO.Constraint(A, -l, COSMO.Nonnegatives)
             constraints = [constraint1; constraint2]
             @testset "$(solver_names[i]): simple problem" for i = 1:length(solver_types)
                 model = COSMO.Model{T}()
                 settings = COSMO.Settings{T}(kkt_solver = solver_types[i])
                 assemble!(model, P, q, constraints, settings = settings)
                 res = COSMO.optimize!(model);
                 @test isapprox(norm(res.x - [0.3; 0.7]), 0., atol=1e-3)
                 @test isapprox(res.obj_val, 1.88, atol=1e-3)
             end
         end

         # Check the KKT triangle assembly utility functions
         @testset "KKT matrix assembly" begin
             n = 3
             m = 2
             P = sparse([1. 0 1; 0 0 1; 1 1 1])
             A = sparse([2. 2 2; 0 3 0])
             # known number of nonzero entries in upper and lower triangle of KKT matrix K
             Kcolnz_upper_sol = [1; 1; 3; 4; 2]
             Kcolnz_lower_sol = [3; 4; 2; 1; 1]
             Kcolnz_upper = zeros(Int64, m + n)
             COSMO._count_upper_triangle!(Kcolnz_upper, P, A, n)
             @test Kcolnz_upper == Kcolnz_upper_sol
             Kcolnz_lower = zeros(Int64, m + n)
             COSMO._count_lower_triangle!(Kcolnz_lower, P, A, n)
             @test Kcolnz_lower == Kcolnz_lower_sol

             # make sure all the assembly function return the same matrix K
             m = 11
             n = 15
             P  = sparse(generate_pos_def_matrix(rng, n, MT = T))
             A  = sprand(rng, T, m, n, 0.2)
             sigma = T(1e-6)
             rho = T(0.1) * ones(T, m)
             K_full = COSMO._assemble_kkt_full(P, A, sigma, rho)
             K_tril = COSMO.assemble_kkt_triangle(P, A, sigma, rho, :L)
             K_triu = COSMO.assemble_kkt_triangle(P, A, sigma, rho, :U)
             @test triu(K_full) == K_triu
             @test tril(K_full) == K_tril
         end
    end
end
nothing


text:
if VERSION >= v"0.7"
  using DelimitedFiles
end

include("../../../src/transit_poly_struct.jl")

function profile_transit_poly!(trans,flux,b)
    t1 = time_ns()
    for k=1:10
        for i=1:length(b)
          trans.b = b[i]
          flux[i] = transit_poly_g(trans)
        end
    end
    elapsed = 1e-10 * (time_ns() - t1)
    return elapsed
end

b = readdlm("b.txt")
u = readdlm("u.txt")
u_n = zeros(length(u))
for i = 1:length(u)
    u_n[i] = u[i]
end
trans = transit_init(0.1, 0.0, u_n, false)
flux = zeros(length(b))
println(profile_transit_poly!(trans, flux, b))
writedlm("flux.txt", flux)

# Multiprecision
trans_big = transit_init(big(0.1), big(0.0), big.(u_n), false)
flux_big = zeros(BigFloat, length(b))
for i=1:length(b)
  trans_big.b = big(b[i])
  flux_big[i] = transit_poly_g(trans_big)
end
writedlm("flux_multi.txt", flux_big)


text:
module LowFrequencyEconometrics

#Import your packages
using KernelDensity,
		Distributions,
		Dierckx,
		StatsBase,
		DataFrames,
		PyPlot
#import Base: func1 #Any function you add dispatches to need to be imported directly

#abstract AbType #Define abstract types before the types they abstract!
include("util\\psi_compute.jl")
include("util\\simul_.jl")
include("util\\den_invariant.jl")
include("util\\frac_cov.jl")
include("util\\acv2cov.jl")
include("util\\ar1_cov.jl")
include("bcdModel\\bcdModel.jl") #Include all the functionality
include("bcdModel\\bcd_spectrum.jl") #Include all the functionality
include("bcdModel\\Sigma_Compute.jl") #Include all the functionality
include("tests\\lfH_test.jl")
include("tests\\lfS_test.jl")
include("tests\\lfst_test.jl")
include("tests\\lfur_test.jl")

export
#from "util\\simul_.jl"
simul_arma,
simul_llm,
#from "util\\frac_cov.jl"
frac_cov,
#from "util\\psi_compute.jl"
psi_compute,
#from "util\\den_invariant.jl"
den_invariant,
#from "util\\acv2cov.jl"
acv2cov,
#from "util\\ar1_cov.jl"
ar1_cov,
#from "bcdModel\\bcd_spectrum.jl"
plot_locS, #Export the functions you want users to use
#from "bcdModel\\bcdModel.jl"
bcdModel,
#from "bcdModel\\Sigma_Compute.jl"
q_Sigma_cd_grid,
interp_q_Sigma_cd_grid,
estimate_d,
Sigma_Compute2,
#from "tests\\lfH_test.jl"
#from "tests\\lfS_test.jl"
#from "tests\\lfst_test.jl"
lfst_test,
lfst_pv,
mw_lfst_pv
#from "tests\\lfur_test.jl"

end


text:
using DataStructures: OrderedDict
using Flux, CuArrays
export Game, showAll, validMoves, randomPlace, place, getGameState, int, PLAYERS

## Game setting
# board size
BOARD_SIZE = Int8(9)
ACTION_SIZE = Int8(9)
HISTORY = Int8(4)

# trainsformation of matrix mapping
ALL_MOVES = (
    original_dont_move  = Int8[1, 2, 3, 4, 5, 6, 7, 8, 9],
    original_dim1_flip  = Int8[7, 8, 9, 4, 5, 6, 1, 2, 3],
    right_rotate_1x90°  = Int8[7, 4, 1, 8, 5, 2, 9, 6, 3],
    right_rot1x90_flip  = Int8[9, 6, 3, 8, 5, 2, 7, 4, 1],
    right_rotate_2x90°  = Int8[9, 8, 7, 6, 5, 4, 3, 2, 1],
    right_rot2x90_flip  = Int8[3, 2, 1, 6, 5, 4, 9, 8, 7],
    right_rotate_3x90°  = Int8[3, 6, 9, 2, 5, 8, 1, 4, 7],
    right_rot3x90_flip  = Int8[1, 4, 7, 2, 5, 8, 3, 6, 9]
)

# transition from players
PLAYERS = Dict{Union{String, Int8}, Union{String, Int8}}("O"=>Int8(1), Int8(1)=>"O", "X"=>Int8(-1), Int8(-1)=>"X", zero(Int8)=>".")
# NEXT_TURN = Dict{Union{String, Int8}, Union{String, Int8}}("O"=>"X", Int8(1)=>Int8(-1), "X"=>"O", Int8(-1)=>Int8(1))

# The state we need to evaluate winning situtaions
STATE = Array{Array{Int8,1}}([
        [1, 2, 3], [4, 5, 6], [7, 8, 9],
        [1, 4, 7], [2, 5, 8], [3, 6, 9],
        [1, 5, 9], [3, 5, 7]
        ])

# give an action, get all possible winning situations
SLOT = Array{Array{Array{Int8,1},1},1}([
        [[1, 2, 3], [1, 4, 7], [1, 5, 9]],
        [[1, 2, 3], [2, 5, 8]],
        [[1, 2, 3], [3, 5, 7], [3, 6, 9]],
        [[1, 4, 7], [4, 5, 6]],
        [[4, 5, 6], [2, 5, 8], [1, 5, 9], [3, 5, 7]],
        [[4, 5, 6], [3, 6, 9]],
        [[1, 4, 7], [7, 8, 9], [3, 5, 7]],
        [[7, 8, 9], [2, 5, 8]],
        [[7, 8, 9], [3, 6, 9], [1, 5, 9]]
        ])

# offset move -> actual board index
OFFSET = Dict{Tuple{Int8,Int8},Tuple{Int8,Int8}}(
  (3, 5) => (2, 8), (3, 6) => (2, 9), (6, 9) => (6, 9), (4, 4) => (5, 1), (9, 4) => (8, 7), (9, 1) => (7, 7), (8, 9) => (9, 6), (3, 1) => (1, 7), (4, 5) => (5, 2),
  (9, 9) => (9, 9), (2, 4) => (2, 4), (6, 5) => (5, 8), (4, 9) => (6, 3), (8, 4) => (8, 4), (1, 2) => (1, 2), (9, 2) => (7, 8), (5, 1) => (4, 4), (9, 8) => (9, 8),
  (3, 4) => (2, 7), (9, 6) => (8, 9), (8, 2) => (7, 5), (2, 3) => (1, 6), (2, 6) => (2, 6), (5, 4) => (5, 4), (2, 1) => (1, 4), (2, 5) => (2, 5), (8, 5) => (8, 5),
  (5, 9) => (6, 6), (1, 4) => (2, 1), (3, 9) => (3, 9), (2, 8) => (3, 5), (5, 8) => (6, 5), (4, 8) => (6, 2), (7, 8) => (9, 2), (1, 1) => (1, 1), (7, 7) => (9, 1),
  (7, 9) => (9, 3), (5, 7) => (6, 4), (7, 2) => (7, 2), (2, 2) => (1, 5), (6, 3) => (4, 9), (1, 7) => (3, 1), (4, 7) => (6, 1), (6, 1) => (4, 7), (6, 8) => (6, 8),
  (8, 8) => (9, 5), (7, 1) => (7, 1), (1, 3) => (1, 3), (7, 5) => (8, 2), (2, 9) => (3, 6), (3, 2) => (1, 8), (7, 6) => (8, 3), (9, 3) => (7, 9), (3, 7) => (3, 7),
  (6, 6) => (5, 9), (5, 5) => (5, 5), (2, 7) => (3, 4), (4, 1) => (4, 1), (8, 3) => (7, 6), (4, 2) => (4, 2), (1, 6) => (2, 3), (4, 3) => (4, 3), (9, 7) => (9, 7),
  (5, 6) => (5, 6), (8, 6) => (8, 6), (4, 6) => (5, 3), (8, 1) => (7, 4), (3, 3) => (1, 9), (6, 2) => (4, 8), (7, 4) => (8, 1), (1, 9) => (3, 3), (8, 7) => (9, 4),
  (6, 4) => (5, 7), (5, 2) => (4, 5), (1, 8) => (3, 2), (6, 7) => (6, 7), (7, 3) => (7, 3), (5, 3) => (4, 6), (3, 8) => (3, 8), (1, 5) => (2, 2), (9, 5) => (8, 8)
)

int(x::Union{Char, String}, T=Int) = parse(T, x)
convertState(state) = [int(i, Int8) for i in state]

mutable struct Game
    boards::AbstractArray{Float32,3}
    vmoves::Dict{Int8, AbstractArray}
    nextPlayer::Int8
    steps::Int8
    winner::Int8
    state::String
    function Game()
        boards = zeros(Float32, (BOARD_SIZE, BOARD_SIZE, 2HISTORY+1))
        vmoves = Dict(Int8(1)=>trues(9), Int8(-1)=>trues(9), Int8(0)=>zeros(Int8, 9,9))
        new(boards, vmoves, Int8(-1),  Int8(1), Int8(0), "")
    end
end

function nozero(x::AbstractArray)
    for i in x
        i==0 && return false
    end
    return true
end

@inline function rmInvalid!(self::Game, x::Integer, y::Integer)
    if @inbounds self.vmoves[0][x,y] != 0
        self.winner = -self.nextPlayer
        return
    end
    @inbounds self.vmoves[0][x,y] = self.nextPlayer

    @inbounds for triple in SLOT[y]
        @inbounds code = sum(self.vmoves[0][x, triple])
        if  code == 2self.nextPlayer
            @inbounds self.vmoves[-self.nextPlayer][x] = false
        end
        if @inbounds code ==-self.nextPlayer && nozero(self.vmoves[0][x, triple])
            count = 0
            for t in STATE
                @inbounds count+= (sum(self.vmoves[0][x, t]) == -2self.nextPlayer)
                count>0 && break
            end
            count >0 && continue
            @inbounds self.vmoves[self.nextPlayer][x] = true
        end
        if @inbounds code==3self.nextPlayer
            self.winner = self.nextPlayer
            return
        end
    end
    return
end

function Game(state::String)
    self = Game()
    isempty(state) && return self
    moves = convertState(state)
    num_moves =  length(moves)
    subboard = moves[1]
    update_bound = max(2, num_moves - HISTORY+1)
    upp = max(1, HISTORY+2-num_moves)
    @sync for i=2:num_moves
        @inbounds rmInvalid!(self, subboard, moves[i])
        @inbounds x,y = OFFSET[(subboard, moves[i])]
        @inbounds self.boards[x, y, upp:HISTORY] .= self.nextPlayer
        @inbounds self.boards[x, y, HISTORY+1:end-(upp-1)] .= -self.nextPlayer
        upp += (i ≥ update_bound);
        self.nextPlayer *= -1;
        @inbounds subboard = moves[i]
    end
    self.steps = num_moves
    @inbounds self.boards[:,:,HISTORY+1] .= self.nextPlayer
    self.state = state
    return self
end

function Base.show(io::IO, self::Game)
    "print out the board"
    compact = get(io, :compact, false)
    if !compact
        board = self.boards[:,:,HISTORY]
        for i=1:9
            print(io, PLAYERS[board[i,1]]," ", PLAYERS[board[i,2]]," ",PLAYERS[board[i,3]]," | ")
            print(io, PLAYERS[board[i,4]]," ", PLAYERS[board[i,5]]," ",PLAYERS[board[i,6]]," | ")
            print(io, PLAYERS[board[i,7]]," ", PLAYERS[board[i,8]]," ",PLAYERS[board[i,9]],"\n")
            if i in (3, 6)
                print(io, "------+-------+------\n")
            end
        end
        showState = isempty(self.state) ? "0" : self.state
        print(io, "\nSubboard: ",showState[end],"; next Player: ", PLAYERS[self.nextPlayer], ";\n")
        print(io, "\nGame States: ",showState,"; Steps: ",self.steps-1,"; ")
        iszero(self.winner) ? print(io, "Game continue...\n") : print(io, PLAYERS[self.winner], " wins\n")
    end
end

function checkWin(self::Game, move::Integer)
    @inbounds x = int(self.state[end])
    @inbounds for triple in SLOT[move]
        @inbounds if sum(self.vmoves[0][x, triple]) == 2self.nextPlayer
            return true
        end
    end
    return false
end

function validMoves(self::Game)
    "Returns all the legal moves for current subboard"
    isempty(self.state) && return keys(OFFSET)
    @inbounds x = int(self.state[end])
    @inbounds available_moves = Array{Int8}(findall(iszero, self.vmoves[0][x,:]))
    if @inbounds !self.vmoves[-self.nextPlayer][x]
        return filter(available_moves) do move
            checkWin(self, move)
        end
    end
    validmoves = filter(available_moves) do move
        @inbounds self.vmoves[self.nextPlayer][move]
    end
    @inbounds if isempty(validmoves)
        return isempty(available_moves) ? [Int8(1)] : available_moves
    end
    return validmoves
end

function place(self::Game, move::Union{Tuple, Integer})
    !iszero(self.winner) && return false
    if isempty(self.state)
        @inbounds rmInvalid!(self, move[1], move[2])
        @inbounds self.state = "$(move[1])$(move[2])"
    else
        @inbounds rmInvalid!(self, int(self.state[end]), move)
        x, y = OFFSET[(int(self.state[end]), move)]
        if self.boards[x, y, HISTORY] != 0
            self.winner = -self.nextPlayer
            return false
        end
        self.state = "$(self.state)$move"
    end
    num_moves = length(self.state)
    iters = max(1, HISTORY-num_moves+2)
    player = self.nextPlayer
    for i=HISTORY:-1:iters
        @inbounds subboard = int(self.state[end+i-HISTORY-1])
        @inbounds x, y = OFFSET[(subboard, int(self.state[end+i-HISTORY]))]
        @inbounds self.boards[x, y, i] = player
        @inbounds self.boards[x, y, end-i+1] = -player
        player *= -1
    end
    self.steps += 1
    self.nextPlayer *= -1
    @inbounds self.boards[:,:,HISTORY+1] .= self.nextPlayer
    return true
end

function randomPlace(self::Game, times::Integer=1)
    for i = 1:times
        !iszero(self.winner) && break
        moves = validMoves(self)
        place(self, rand(moves))
    end
end

function conv(state::String, typ::Int)
    result = ""
    for (i,e) in enumerate(state)
        result *= string(ALL_MOVES[typ][int(e)])
    end
    return result
end

function getGameState(self::Game)
    return reshape(self.boards, (BOARD_SIZE, BOARD_SIZE, 2HISTORY+1,1))
end

function getAugmentation(self::Game, p::Array{Float32, 1})
    """ The order is
        {game.state => M1, flip1 => M2,
         rotr1x90   => M3, flip2 => M4
         rotr2x90   => M5, flip3 => M6
         rotr3x90   => M7, flip4 => M8}
    """
    # strings = String[];
    # push!(strings, conv(self.state, 2))
    data = Array{Array{Float32, 3}, Array{Float32, 1}, Integer}[]
    gameState1 = self.boards
    push!(data, [gameState1, p, game.nextPlayer])
    # flip along with dim1
    # states = cat(dims=4, reverse(gameState1, dims=1))
    @inbounds push!(data, [reverse(gameState1, dims=1), p[ALL_MOVES[2]], game.nextPlayer])
    for i=3:2:8
        # right rotate 90°
        # push!(strings, conv(self.state, i))
        gameState1 = mapslices(rotr90, gameState1, dims=[1,2])
        @inbounds push!(data, [gameState1, p[ALL_MOVES[i]], game.nextPlayer])
        # states = cat(dims=4, states, gameState1)
        # flip along with dim1
        # push!(strings, conv(self.state, i+1))
        @inbounds push!(data, [reverse(gameState1, dims=1), p[ALL_MOVES[i]], game.nextPlayer])
        # states = cat(dims=4, states, reverse(gameState1, dims=1))
    end
    return data
end


text:
module InterfaceTests

using Test
using InplaceArrays.Arrays

a = rand(20,12)

test_array(a,a)
test_array(a,a,≈)

end # module


text:
# vortex particle definition

# import modules
import LinearAlgebra
la = LinearAlgebra

"""
`Particle(args)`

Description here.
"""
mutable struct Particle

    # properties here
    gamma::Array{Number, 1}       # vortex strength vector
    x::Array{Number, 1}           # particle position
    v::Array{Number, 1}           # particle velocity

    Particle(gamma,x,v) = new(gamma,x,v)

end

"""
`vinduced(particle::Particle,X)`

Calculates the induced velocity produced by the vortex particle argument at the specified location.
"""
function vinduced(particle,x)
    r = x .- particle.x
    vinduced = -la.cross(r/4/pi/la.norm(r)^3, particle.gamma)

    return vinduced
end

"""
`place(coordinates::Array{Array{Number,1},1},gammas::Array{Array{Number,1},1},vs::Array{Array{Number,1},1})`

Places vortex particles at the specified locations.
"""
function place(coordinates, gammas, vs)
    particles = Particle[]
    for i = range(1,length=length(coordinates))
        push!(particles,Particle(gammas[i],coordinates[i],vs[i]))
    end

    return particles
end

"""
`netv(particles::Array{Particle,1},x::Array{Particle,1})`

Returns the net velocity induced by the vector of Particles at location x.
"""
function netv(particles, x; err = 0.001)
    # find velocity induced at a point due to all vortex particles
    velocity = [0.0, 0.0, 0.0]
    for particle in particles
        if la.norm(particle.x .- x) < err
            continue
        else
            velocity .+= vinduced(particle, x)
        end
    end

    return velocity
end

"""
`step(particles::Array{Particle,1},timestep::Number,Uinf)`

Calculates induced velocities and advances forward the specified timestep.
"""
function advance(particles,timestep,Uinf)
    # update velocity of each particle
    for particle in particles
        # yogi = netv(particles, particle.x) .+ Uinf(particle.x)
        particle.v = netv(particles, particle.x) .+ Uinf(particle.x)#yogi # WHY DO I NEED YOGI????
    end

    # move each particle
    for particle in particles
        particle.x .= particle.x + particle.v * timestep
    end
end

"""
`vfield(particles::Array{Particle,1},xs::Array{Number,1},ys::Array{Number,1},zs::Array{Number,1})`

Evaluates the velocity field at discrete points on a cartesian grid defined by the arguments:

    * `xs` defines the X values of the cartesian grid
    * `ys` defines the Y values of the cartesian grid
    * `zs` defines the Z values of the cartesian grid

Returns six vectors:

    * `gridx` defines the X coordinate of each grid point
    * `gridy` defines the Y coordinate of each grid point
    * `gridz` defines the Z coordinate of each grid point
    * `vx` defines the X component of the vorticity-induced velocity at each grid point
    * `vy` defines the Y component of the vorticity-induced velocity at each grid point
    * `vz` defines the Z component of the vorticity-induced velocity at each grid point
"""
function vfield(particles,xs, ys, zs)
    # build grid
    gridx = Number[]
    gridy = Number[]
    gridz = Number[]
    vx = Number[]
    vy = Number[]
    vz = Number[]

    for x in xs
        for y in ys
            for z in zs
                # add coordinate
                push!(gridx,x)
                push!(gridy,y)
                push!(gridz,z)
                v = netv(particles,[x,y,z]) # get induced velocity
                if la.norm(v) > 2.0
                    v = v./la.norm(v)*2.0
                end
                # unpack velocity components
                push!(vx,v[1])
                push!(vy,v[2])
                push!(vz,v[3])
            end
        end
    end

    return gridx, gridy, gridz, vx, vy, vz
end

"""
`gfield(particles::Array{Particle,1})`

Useful for plotting the circulation vector of each particle

Returns six vectors:

    * `gridx` defines the X coordinate of each particle
    * `gridy` defines the Y coordinate of each particle
    * `gridz` defines the Z coordinate of each particle
    * `gx` defines the X component of the circulation of each particle
    * `gy` defines the Y component of the circulation of each particle
    * `gz` defines the Z component of the circulation of each particle
"""
function gfield(particles)
    # set up vectors
    gridx = Number[]
    gridy = Number[]
    gridz = Number[]
    gx = Number[]
    gy = Number[]
    gz = Number[]
    # get coordinates and circulations
    for particle in particles
        push!(gridx,particle.x[1])
        push!(gridy,particle.x[2])
        push!(gridz,particle.x[3])
        push!(gx,particle.gamma[1])
        push!(gy,particle.gamma[2])
        push!(gz,particle.gamma[3])
    end

    return gridx, gridy, gridz, gx, gy, gz
end

"""
`vorticity(particles::Array{Particle,1}, x::Vector{Number})`

Returns the particle-induced vorticity at coordinates x.
"""
function vorticity(particles, x)
    # update velocity of each particle
    for particle in particles

    end

end


text:
import BSON
import JSON
import Random

Random.seed!(0)

struct Bob
  x::String
  y::Vector{Float64}
end

bob = Bob("Arthur", randn(3))

save_path = "results/julia_bob.bson"
BSON.bson(save_path, bob=bob)

bob_in = BSON.load(save_path)

# TODO: Try JSON?
# jbob = JSON.json(bob)
# jbob


text:
export Epinions

"""
This is a who-trust-whom online social network of a a general consumer review site Epinions.com.
Members of the site can decide whether to ''trust'' each other. All the trust relationships interact
and form the Web of Trust which is then combined with review ratings to determine which reviews
are shown to the user.

This network has:
Vertices      : 75890
Edges         : 508837
Vertex Props  : -
Edge Props    : -
"""
abstract Epinions <: Dataset

function download(::Type{Epinions})
   dpath = filepath("Epinions.txt")
   tpath = "/tmp/soc-Epinions1.txt"

   if isfile(dpath)
      return dpath
   end

   if !isfile(tpath)
      download("http://snap.stanford.edu/data/soc-Epinions1.txt.gz", "/tmp/soc-Epinions1.txt.gz")
      run(`gzip -d /tmp/soc-Epinions1.txt.gz`)
   end

   Nv = 75890
   Ne = 508837

   fin = open(tpath)
   vs  = readdlm(fin, Int)
   close(fin)

   fout = open(dpath, "w")
   write(fout, "$Nv,$Ne\n")
   write(fout, "\n")
   write(fout, "\n")
   write(fout, "\n")
   write(fout, "\n")
   writedlm(fout, collect(0:(Nv-1)))
   writecsv(fout, vs)
   close(fout)

   dpath
end


text:
# The following is a particular combination of generic parameters such that
# the goes close to a singularity. An increase in the precision of the residual
# is only necessary for 1 or 2 steps
@testset "Steiner High Precision" begin
    p = Complex{Float64}[
        -0.8760865081446109-0.8245161978967308im,
        -0.371989298448952+1.175807247192591im,
        -0.5889429565516365-0.6815249624614181im,
        0.05893625668644283-0.15242494424272904im,
        -0.04017905091825077-0.32922740646951476im,
        -1.1620912695313244-0.880462208384591im,
        -0.6128442732846304+0.2299959308120725im,
        0.04451305089395737-0.5613666038122553im,
        0.4275591895528981-0.4560517627610251im,
        -1.44959201510603+0.9148102094587954im,
        0.9642383181691287-1.028506672030121im,
        -0.04531547209868653+1.1658486088983617im,
        -0.7706473388214675+1.6213385645828058im,
        0.3905889030007126-0.837367437811382im,
        0.521976803491032+0.8715564019511688im,
        0.9913423747115939-0.9251461816473662im,
        -1.5976395911608239-0.6869570695244973im,
        0.3812500672931469-0.2640531551530478im,
        0.6050593029483629+1.0735839629840869im,
        0.24229108273519206+0.31651798241470175im,
        0.2343873046310603+0.9209493301878944im,
        -0.04243469809131875+0.2930962153305209im,
        -0.05095304883782437-0.07163754176755009im,
        0.055623959461651266+0.11458865587499502im,
        0.4749235468690426+1.3139464606178932im,
        -0.5340685054046019-0.05295858535911429im,
        0.6838868206371786-0.5368790063173585im,
        -0.6012500786841708-0.41545534317815064im,
        0.14315666179481693+0.7966975827842401im,
        0.06217289908817691-1.268704027512965im,
    ]

    q = [
        0.18802370367009363,
        -1.9260285280414342,
        0.6835970218915922,
        -0.5817725614568748,
        -0.5033365002810842,
        0.6894548688621689,
        0.03748633251693155,
        -0.12168863868778622,
        -1.6197914296308757,
        -2.2192790882018585,
        1.4133360740976828,
        1.6335108803709772,
        0.13316023698203722,
        0.6415334044239753,
        0.36899573610352765,
        -0.6026231774364099,
        1.086473021454091,
        -0.7385757989526022,
        1.4247229002869566,
        -0.2974094335738924,
        -0.2651970300878602,
        -1.0758341565830896,
        -0.2676109895616729,
        -1.0203047093858557,
        -0.2897673462238505,
        1.950855352655958,
        -1.3290748919968962,
        -0.6377519368264665,
        0.850257024568716,
        0.2844129100094346,
    ]

    s_p = Complex{Float64}[
        5.7580258277373275-4.533830475896814im,
        -4.471642985901334+11.790406705747543im,
        -3.1685796570857536-4.127973850160156im,
        -3.5677280267248253-1.4277434015840338im,
        1.0995526935530473-3.7157339191392107im,
        0.05304893282894325-8.310443430899932im,
        14.954165680284154-4.813664506429587im,
        -1.0940288937042102+1.525825504520944im,
        -1.008339359846099+0.3508156442116616im,
        0.07785841864101696-0.09487678203891686im,
        0.04803901724013503-0.35755534225201563im,
        -0.16148337320880576+0.09850690772683601im,
        -0.02859806047287574-0.3486320723561723im,
        -0.5931267530551303-0.49056820640225385im,
        -0.1700776249816018-1.1761079973154294im,
    ]

    s_q = Complex{Float64}[
        0.7971076145172316+0.5372402576809401im,
        0.2071466227304207+1.259009278293226im,
        -0.3398858196752808+0.3751960862136744im,
        -1.3675207743723174-0.09914031977945449im,
        -0.9246609053707422-1.25038329814984im,
        0.391715711840247+0.2299495066551252im,
        0.29477984269056495-0.26961698034518383im,
        1.0814832524360534-0.0564333571343061im,
        0.46580859060989904+0.5449535876793804im,
        -2.0668707771828707+4.099238639488696im,
        -1.2045014977210438-6.995994517223501im,
        0.6499685548019725+0.22568940175977828im,
        -1.094860871004778-2.07591025239193im,
        0.14092679017444423-0.6261846990179475im,
        0.8690167435198312-0.5157825273396847im,
    ]
    F = let
        @var x[1:2] a[1:5] c[1:6] y[1:2, 1:5]

        #tangential conics
        f =
            a[1] * x[1]^2 +
            a[2] * x[1] * x[2] +
            a[3] * x[2]^2 +
            a[4] * x[1] +
            a[5] * x[2] +
            1
        ∇ = differentiate(f, x)
        #5 conics
        g =
            c[1] * x[1]^2 +
            c[2] * x[1] * x[2] +
            c[3] * x[2]^2 +
            c[4] * x[1] +
            c[5] * x[2] +
            c[6]
        ∇_2 = differentiate(g, x)
        #the general system
        #f_a_0 is tangent to g_b₀ at x₀
        function Incidence(f, a₀, g, b₀, x₀)
            fᵢ = f(x => x₀, a => a₀)
            ∇ᵢ = ∇(x => x₀, a => a₀)
            Cᵢ = g(x => x₀, c => b₀)
            ∇_Cᵢ = ∇_2(x => x₀, c => b₀)
            [fᵢ; Cᵢ; det([∇ᵢ ∇_Cᵢ])]
        end
        @var v[1:6, 1:5]
        F = vcat(map(i -> Incidence(f, a, g, v[:, i], y[:, i]), 1:5)...)
        System(F, [a; vec(y)], vec(v))
    end

    @testset "Steiner - AD: $AD" for AD = 1:3
        tracker = Tracker(
            ParameterHomotopy(F, p, q),
            options = TrackerOptions(automatic_differentiation = AD),
        )
        r_q = track(tracker, s_p, 1, 0)
        @test is_success(r_q)
        @test solution(r_q) ≈ s_q
        r_p = track(tracker, r_q, 0, 1)
        @test is_success(r_p)
        @test solution(r_p) ≈ s_p
    end
end


text:
using IonBase
using IonBase.Templates
using Configurations

include(joinpath(pkgdir(IonBase), "test", "utils.jl"))
with_test_ion() do
    IonBase.copy_templates()
end

with_test_ion() do
    IonBase.create("Foo"; user="me", template="package", authors="A,B")
end

t = Template("basic", [
    CompatHelper(),
    Git(),
    License(),
    ProjectFile(),
    Readme(),
    SrcDir(),
    TagBot(),
    Tests(),
])

toml("templates/basic.toml", t)

t = Template("academic", [
    CompatHelper(),
    Git(),
    License(),
    ProjectFile(),
    Readme(),
    SrcDir(),
    TagBot(),
    Tests(),
    Citation(),
])

toml("templates/academic.toml", t)

t = Template("package", [
    CompatHelper(),
    Git(),
    License(),
    ProjectFile(),
    Readme(),
    SrcDir(),
    TagBot(),
    Tests(),
    GitHubActions(),
    Codecov(),
    Documenter(;depoly="GitHubActions"),
])

toml("templates/package.toml", t)

t = Template("comonicon", [
    CompatHelper(),
    Git(),
    License(),
    ProjectFile(),
    Readme(),
    SrcDir(),
    TagBot(),
    Tests(),
    GitHubActions(),
    Codecov(),
    Documenter(;depoly="GitHubActions"),
    Comonicon(),
])

toml("templates/comonicon.toml", t)

t = Template("comonicon-sysimg", [
    CompatHelper(),
    Git(),
    License(),
    ProjectFile(),
    Readme(),
    SrcDir(),
    TagBot(),
    Tests(),
    GitHubActions(),
    Codecov(),
    Documenter(;depoly="GitHubActions"),
    Comonicon(),
    SystemImage(),
])

toml("templates/comonicon-sysimg.toml", t)


text:
using DataStructures
using JSON

#TODO: Apagar as linhas para remover após consertar o matching
@enum(Tokens,
      CT_VALUE=1, # ??
      CT_INT,       # 123123
      CT_CHAR,      # 'c'
      CT_FLOAT,     # 123123.123123f
      #CT_VEC,       # REMOVER
      CT_STRING,    # "foo bar"

      #SMCL,         # ??
      EPSILON,          # epsilonILON
      ID,           # 'id'
      CONST,        #  'const'
      #EXPR,         # REMOVER

      OPR_PM,       # '+ -'
      OPR_DM,       # ' * /'

      OPRLN,        # ??
      #OPRLR_EQ,     # REMOVER
      #OPRLR_LG,     # REMOVER
      OPRLR_LGEQ,   # '== >= <= !='
      OPR_ATR,      # '='
      #OPRLR_UN_NEG, # REMOVER
      OPRLR_AND,    # 'and'
      OPRLR_OR,     # 'or'
      OPRL_NOT,     # 'not'
      FN_PRINT,     # "print"
      #FN_READ,      # REMOVER
      DRED,         # COMENTARIO ??

      IF,       # 'if'
      ELSE,      # 'else'
      FOR,      # 'for'
      WHILE,    # 'while'
      CONTINUE,     # 'continue'
      BREAK,        # 'break'
      RETURN,       # 'return'

      COMMA,        # ','

      O_BRCKT,      # '('
      C_BRCKT,      # ')'
      O_SQRBRCKT,   # '['
      C_SQRBRCKT,   # ']'
      O_C_BRCKT,    # '{'
      C_C_BRCKT,    # '}'
      #O_PRTSIS,    # REMOVER
      #C_PRTSIS,    # REMOVER

      IDT_INT,      # 'int'
      IDT_FLOAT,    # 'float'
      IDT_CHAR,     # 'char'
      IDT_STRING,   # 'string'
      IDT_BOOL,     # 'bool'
      IDT_VOID,     # 'void'
      CT_FALSE,     # 'false'
      CT_TRUE,      # 'true'

      FN_MAIN,      # "main"
      VEC_IN,       # '::'
      LEX_ERR,      # erro lexico
      #=
      OPR_SUM,      # REMOVER
      OPR_SUB,      # REMOVER
      CTN,          # REMOVER
      PARAMS_Tk,    # REMOVER
      EXP_BOOL,     # REMOVER

      FN_DCLR,      # REMOVER
      FN_CALL_TK,   # REMOVER
      OPRLR_GT,     # REMOVER
      OPRLR_GEQ,    # REMOVER
      OPR_CONCAT,   # REMOVER
      =#
      OPRLR_REL,     #Operador lógico relacional
      CT_B,         # true ou false
      CMNT_LN,      # COMENTARIO DE LINHA
      #OPRLR_DIF,    # REMOVER
      EOF)

tks_names = String[]

for tk in instances(Tokens)
	push!(tks_names,string(tk))
end
for ti =1:length(tks_names)
      tks_names[ti] = lowercase(tks_names[ti])
end
#map!(lowercase,tks_names);

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
@enum mjtWarning begin             # warning types
   WARN_INERTIA      = 0         # (near) singular inertia matrix
   WARN_CONTACTFULL              # too many contacts in contact list
   WARN_CNSTRFULL                # too many constraints
   WARN_VGEOMFULL                # too many visual geoms
   WARN_BADQPOS                  # bad number in qpos
   WARN_BADQVEL                  # bad number in qvel
   WARN_BADQACC                  # bad number in qacc
   WARN_BADCTRL                  # bad number in ctrl

   NWARNING                      # number of warnings
end

@enum mjtTimer begin
   # main api
   TIMER_STEP        = 0         # step
   TIMER_FORWARD                 # forward
   TIMER_INVERSE                 # inverse

   # breakdown of step/forward
   TIMER_POSITION                # fwdPosition
   TIMER_VELOCITY                # fwdVelocity
   TIMER_ACTUATION               # fwdActuation
   TIMER_ACCELERATION            # fwdAcceleration
   TIMER_CONSTRAINT              # fwdConstraint

   # breakdown of fwdPosition
   TIMER_POS_KINEMATICS          # kinematics, com, tendon, transmission
   TIMER_POS_INERTIA             # inertia computations
   TIMER_POS_COLLISION           # collision detection
   TIMER_POS_MAKE                # make constraints
   TIMER_POS_PROJECT             # project constraints

   NTIMER                        # number of timers
end

#------------------------------ mjContact ----------------------------------------------
struct mjContact
   dist::mjtNum
   pos::SVector{3, mjtNum}
   frame::SVector{9, mjtNum}
   includemargin::mjtNum
   friction::SVector{5, mjtNum}
   solref::SVector{NREF, mjtNum}
   solimp::SVector{NIMP, mjtNum}
   mu::mjtNum
   H::SVector{36, mjtNum}
   dim::Cint
   geom1::Cint
   geom2::Cint
   exclude::Cint
   efc_address::Cint
end

#------------------------------ diagnostics --------------------------------------------

struct mjWarningStat
   lastinfo::Cint
   number::Cint
end

struct mjTimerStat
   duration::mjtNum
   number::Cint
end

struct mjSolverStat
   improvement::mjtNum
   gradient::mjtNum
   lineslope::mjtNum
   nactive::Cint
   nchange::Cint
   neval::Cint
   nupdate::Cint
end

#---------------------------------- mjData ---------------------------------------------

mutable struct mjData
   nstack::Cint
   nbuffer::Cint

   pstack::Cint

   maxuse_stack::Cint
   maxuse_con::Cint
   maxuse_efc::Cint

   warning::SVector{Int(NWARNING), mjWarningStat}
   timer::SVector{Int(NTIMER), mjTimerStat}
   solver::SVector{NSOLVER, mjSolverStat}
   solver_iter::Cint
   solver_nnz::Cint
   solver_fwdinv::SVector{2, mjtNum}

   ne::Cint
   nf::Cint
   nefc::Cint
   ncon::Cint
   time::mjtNum
   energy::SVector{2, mjtNum}

   buffer::Ptr{Cvoid}
   stack::Ptr{mjtNum}

   qpos::Ptr{mjtNum}
   qvel::Ptr{mjtNum}
   act::Ptr{mjtNum}
   qacc_warmstart::Ptr{mjtNum}

   ctrl::Ptr{mjtNum}
   qfrc_applied::Ptr{mjtNum}
   xfrc_applied::Ptr{mjtNum}

   qacc::Ptr{mjtNum}
   act_dot::Ptr{mjtNum}

   mocap_pos::Ptr{mjtNum}
   mocap_quat::Ptr{mjtNum}

   userdata::Ptr{mjtNum}

   sensordata::Ptr{mjtNum}

   xpos::Ptr{mjtNum}
   xquat::Ptr{mjtNum}
   xmat::Ptr{mjtNum}
   xipos::Ptr{mjtNum}
   ximat::Ptr{mjtNum}
   xanchor::Ptr{mjtNum}
   xaxis::Ptr{mjtNum}
   geom_xpos::Ptr{mjtNum}
   geom_xmat::Ptr{mjtNum}
   site_xpos::Ptr{mjtNum}
   site_xmat::Ptr{mjtNum}
   cam_xpos::Ptr{mjtNum}
   cam_xmat::Ptr{mjtNum}
   light_xpos::Ptr{mjtNum}
   light_xdir::Ptr{mjtNum}

   subtree_com::Ptr{mjtNum}
   cdof::Ptr{mjtNum}
   cinert::Ptr{mjtNum}

   ten_wrapadr::Ptr{Cint}
   ten_wrapnum::Ptr{Cint}
   ten_J_rownnz::Ptr{Cint}
   ten_J_rowadr::Ptr{Cint}
   ten_J_colind::Ptr{Cint}
   ten_length::Ptr{mjtNum}
   ten_J::Ptr{mjtNum}
   wrap_obj::Ptr{Cint}
   wrap_xpos::Ptr{mjtNum}

   actuator_length::Ptr{mjtNum}
   actuator_moment::Ptr{mjtNum}

   crb::Ptr{mjtNum}
   qM::Ptr{mjtNum}
   qLD::Ptr{mjtNum}
   qLDiagInv::Ptr{mjtNum}
   qLDiagSqrtInv::Ptr{mjtNum}

   contact::Ptr{mjContact}

   efc_type::Ptr{Cint}
   efc_id::Ptr{Cint}
   efc_J_rownnz::Ptr{Cint}
   efc_J_rowadr::Ptr{Cint}
   efc_J_rowsuper::Ptr{Cint}
   efc_J_colind::Ptr{Cint}
   efc_JT_rownnz::Ptr{Cint}
   efc_JT_rowadr::Ptr{Cint}
   efc_JT_rowsuper::Ptr{Cint}
   efc_JT_colind::Ptr{Cint}
   efc_J::Ptr{mjtNum}
   efc_JT::Ptr{mjtNum}
   efc_pos::Ptr{mjtNum}
   efc_margin::Ptr{mjtNum}
   efc_frictionloss::Ptr{mjtNum}
   efc_diagApprox::Ptr{mjtNum}
   efc_KBIP::Ptr{mjtNum}
   efc_D::Ptr{mjtNum}
   efc_R::Ptr{mjtNum}

   efc_AR_rownnz::Ptr{Cint}
   efc_AR_rowadr::Ptr{Cint}
   efc_AR_colind::Ptr{Cint}
   efc_AR::Ptr{mjtNum}

   ten_velocity::Ptr{mjtNum}
   actuator_velocity::Ptr{mjtNum}
   cvel::Ptr{mjtNum}
   cdof_dot::Ptr{mjtNum}
   qfrc_bias::Ptr{mjtNum}
   qfrc_passive::Ptr{mjtNum}
   efc_vel::Ptr{mjtNum}
   efc_aref::Ptr{mjtNum}
   subtree_linvel::Ptr{mjtNum}
   subtree_angmom::Ptr{mjtNum}
   actuator_force::Ptr{mjtNum}
   qfrc_actuator::Ptr{mjtNum}
   qfrc_unc::Ptr{mjtNum}
   qacc_unc::Ptr{mjtNum}

   efc_b::Ptr{mjtNum}
   efc_force::Ptr{mjtNum}
   efc_state::Ptr{Cint}
   qfrc_constraint::Ptr{mjtNum}

   qfrc_inverse::Ptr{mjtNum}

   cacc::Ptr{mjtNum}
   cfrc_int::Ptr{mjtNum}
   cfrc_ext::Ptr{mjtNum}
end

# Callback function types TODO

const mjfGeneric = Ptr{Cvoid}
const mjfConFilt = Ptr{Cvoid}
const mjfSensor = Ptr{Cvoid}
const mjfTime = Ptr{Cvoid}
const mjfAct = Ptr{Cvoid}
const mjfCollision = Ptr{Cvoid}


text:
using Printf, ProgressMeter
@everywhere using MDToolbox

using BSON: @save, @load
@load "msm_hmm_afm.bson" afms models state qs afmize_config probe_radius sigma_noise

configs = [AfmizeConfig(10.0 * (pi / 180),
    r,
    MDToolbox.Point2D(-250, -200),
    MDToolbox.Point2D(250, 200),
    MDToolbox.Point2D(6.25, 6.25),
    MDToolbox.defaultParameters())
    for r in [10.0, 15.0, 20.0, 25.0, 30.0]]

# all frames with parallel map
@time res = @showprogress pmap(x -> getafmposterior(x, models, qs, configs), afms)

@save "msm_hmm_afm.posterior.bson" res configs

max_stars_repo_name: matsunagalab/notebooks

max_issues_repo_name: matsunagalab/notebooks

max_forks_repo_name: matsunagalab/notebooks


text:
X, _ = make_blobs(500, 3, rng=555)

@testset "AffinityPropagation" begin
    m = AffinityPropagation()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    if fp.labels !== nothing
        p  = predict(mach, X)
        @test p isa CategoricalArray
    end
    @test keys(fp) == (:cluster_centers_indices, :cluster_centers, :labels, :affinity_matrix)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
    # first fix #163 on MLJBase
    @test_broken infos[:target_scitype] == AbstractVector{Multiclass}
end

@testset "AgglomerativeClustering" begin
    m = AgglomerativeClustering()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    @test keys(fp) == (:n_clusters, :labels, :n_leaves, :n_connected_components, :children)
    @test fp.labels isa CategoricalArray
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
end

@testset "Birch" begin
    m = Birch()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    p = predict(mach, X)
    Xt = transform(mach, X)
    @test p isa CategoricalArray
    @test Xt isa Tables.MatrixTable
    @test keys(fp) == (:root, :dummy_leaf, :subcluster_centers, :subcluster_labels, :labels)
    @test fp.labels isa CategoricalArray
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
    @test infos[:output_scitype] == Table(Continuous)
    # first fix #163 on MLJBase
    @test_broken infos[:target_scitype] == AbstractVector{Multiclass}
end

@testset "DBSCAN" begin
    m = DBSCAN()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    @test keys(fp) == (:core_sample_indices,  :components, :labels)
    @test fp.labels isa CategoricalArray
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
end

@testset "FeatureAgglomeration" begin
    m = FeatureAgglomeration()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    @test keys(fp) == (:n_clusters, :labels, :n_leaves, :n_connected_components,  :children, :distances)
    @test fp.distances === nothing
    @test fp.labels isa CategoricalArray
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:output_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])

    # NOTE: they're not equal (not sure why)
    Xt  = transform(mach, X)
    Xit = inverse_transform(mach, Xt)
    @test Xit isa Tables.MatrixTable
end

@testset "KMeans" begin
    m = KMeans(n_clusters=4)
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    Xt = transform(mach, X)
    p  = predict(mach, X)
    @test size(fp.cluster_centers) == (4, 3)
    @test keys(fp) == (:cluster_centers, :labels, :inertia)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:output_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
    # first fix #163 on MLJBase
    @test_broken infos[:target_scitype] == AbstractVector{Multiclass}
end

@testset "MBKMeans" begin
    m = MiniBatchKMeans(n_clusters=4)
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    Xt = transform(mach, X)
    p  = predict(mach, X)
    @test size(fp.cluster_centers) == (4, 3)
    @test keys(fp) == (:cluster_centers, :labels, :inertia)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test infos[:output_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
    # first fix #163 on MLJBase
    @test_broken infos[:target_scitype] == AbstractVector{Multiclass}
end

@testset "MeanShift" begin
    m = MeanShift()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    p  = predict(mach, X)
    @test p isa CategoricalArray
    @test keys(fp) == (:cluster_centers, :labels)
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
    # first fix #163 on MLJBase
    @test_broken infos[:target_scitype] == AbstractVector{Multiclass}
end

@testset "OPTICS" begin
    m = OPTICS()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    @test keys(fp) == (:labels, :reachability, :ordering, :core_distances, :predecessor, :cluster_hierarchy)
    @test fp.labels isa CategoricalArray
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
end

@testset "SpectralClustering" begin
    m = SpectralClustering()
    mach = machine(m, X)
    fit!(mach)
    fp = fitted_params(mach)
    @test keys(fp) == (:labels, :affinity_matrix)
    @test fp.labels isa CategoricalArray
    infos = info_dict(m)
    @test infos[:input_scitype] == Table(Continuous)
    @test !isempty(infos[:docstring])
end


text:
#!/usr/bin/env julia

"""
# Problem 1: Multiples of 3 and 5

If we list all the natural numbers below 10 that are multiples of 3 or 5, we
get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of `m` or `n` below `limit`.
"""
module Problem001

include("Sequences.jl")

"""
    solve(; m::Integer = 3, n::Integer = 5, limit::Integer = 1000)

Returns the solution for [`Problem001`](@ref) with the given parameters.

## Preconditions
- `m > 0`
- `n > 0`
- `limit > 0`
"""
function solve(; m::Integer = 3, n::Integer = 5, limit::Integer = 1000)
    m_sum = sum_divisible_by(m, limit)
    n_sum = sum_divisible_by(n, limit)
    lcm_sum = sum_divisible_by(lcm(m, n), limit)
    return m_sum + n_sum - lcm_sum
end

"""
    sum_divisible_by(n::Integer, limit::Integer)

Returns the sum of natural numbers below `limit` that are divisible by `n`.

## Preconditions
- `n > 0`
- `limit > 0`

## Examples
```jldoctest
julia> sum_divisible_by(3, 10)
18

julia> sum_divisible_by(5, 10)
5

julia> sum_divisible_by(5, 5)
0
```
"""
function sum_divisible_by(n::Integer, limit::Integer)
    return Sequences.arithmetic_series(n, (limit - 1) ÷ n, n)
end

end

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "gmsh"
version = v"4.4.1"

# Collection of sources required to build Gmsh
sources = [
    "http://gmsh.info/src/gmsh-$version-source.tgz" =>
    "853c6438fc4e4b765206e66a514b09182c56377bb4b73f1d0d26eda7eb8af0dc",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd gmsh-*
mkdir build
cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DDEFAULT=0 -DENABLE_BUILD_SHARED=1 -DENABLE_MESH=1 -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    MacOS(:x86_64),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libgmsh", :libgmsh),
#   ExecutableProduct(prefix, "gmsh", :gmsh)
]

# Dependencies that must be installed before this package can be built
dependencies = [
   "https://github.com/JuliaLinearAlgebra/OpenBLASBuilder/releases/download/v0.3.0-3/build_OpenBLAS.v0.3.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
function compute_wavelet_filter(Type, par)
    """
        compute_wavelet_filter - Generate Orthonormal QMF Filter for Wavelet Transform

           [h,g] = compute_wavelet_filter(Type,Par)

         Inputs
           Type   string, 'Haar', 'Beylkin', 'Coiflet', 'Daubechies',
                  'Symmlet', 'Vaidyanathan','Battle'
           Par    integer, it is a parameter related to the support and vanishing
                  moments of the wavelets, explained below for each wavelet.

        Outputs
          h   low pass quadrature mirror filter
          g   high pass

         Description
           The Haar filter (which could be considered a Daubechies-2) was the
           first wavelet, though not called as such, and is discontinuous.

           The Beylkin filter places roots for the frequency response function
           close to the Nyquist frequency on the real axis.

           The Coiflet filters are designed to give both the mother and father
           wavelets 2*Par vanishing moments; here Par may be one of 1,2,3,4 or 5.

           The Daubechies filters are minimal phase filters that generate wavelets
           which have a minimal support for a given number of vanishing moments.
           They are indexed by their length, Par, which may be one of
           2,4,6,8,10,12,14,16,18 or 20. The number of vanishing moments is par/2.

           Symmlets are also wavelets within a minimum size support for a given
           number of vanishing moments, but they are as symmetrical as possible,
           as opposed to the Daubechies filters which are highly asymmetrical.
           They are indexed by Par, which specifies the number of vanishing
           moments and is equal to half the size of the support. It ranges
           from 4 to 10.

           The Vaidyanathan filter gives an exact reconstruction, but does not
           satisfy any moment condition.  The filter has been optimized for
           speech coding.

           The Battle-Lemarie filter generate spline orthogonal wavelet basis.
           The parameter Par gives the degree of the spline. The number of
           vanishing moments is Par+1.

        See Also
           FWT_PO, IWT_PO, FWT2_PO, IWT2_PO, WPAnalysis

        References
            The books by Daubechies and Wickerhauser.

        Warning : only Daubechies implemented for the moment !
    """

    if Type == "Daubechies"

        if par == 1
            f = [1,1]/sqrt(2)
        end

        if par == 4
            f = [.482962913145,.836516303738,
                .224143868042,-.129409522551]
        end

        if par == 6
            f = [.332670552950,.806891509311,
            .459877502118,-.135011020010,
                -.085441273882,.035226291882]
        end

        if par == 8
            f = [ .230377813309,.714846570553,
                .630880767930,-.027983769417,
                -.187034811719,.030841381836,
                .032883011667,-.010597401785]
        end

        if par == 10
            f = [.160102397974,.603829269797,.724308528438,
                .138428145901,-.242294887066,-.032244869585,
                .077571493840,-.006241490213,-.012580751999,
                .003335725285]
        end

        if par == 12
            f = [.111540743350,.494623890398,.751133908021,
                .315250351709,-.226264693965,-.129766867567,
                .097501605587,.027522865530,-.031582039317,
                .000553842201,.004777257511,-.001077301085]
        end

        if par == 14
            f = [.077852054085,.396539319482,.729132090846,
                .469782287405,-.143906003929,-.224036184994,
                .071309219267,.080612609151,-.038029936935,
                -.016574541631,.012550998556,.000429577973,
                -.001801640704,.000353713800]
        end

        if par == 16
            f = [.054415842243,.312871590914,.675630736297,
                .585354683654,-.015829105256,-.284015542962,
                .000472484574,.128747426620,-.017369301002,
                -.044088253931,.013981027917,.008746094047,
                -.004870352993,-.000391740373,.000675449406,
                -.000117476784]
        end

        if par == 18
            f = [.038077947364,.243834674613,.604823123690,
                .657288078051,.133197385825,-.293273783279,
                -.096840783223,.148540749338,.030725681479,
                -.067632829061,.000250947115,.022361662124,
                -.004723204758,-.004281503682,.001847646883,
                .000230385764,-.000251963189,.000039347320]
        end

        if par == 20
            f = [.026670057901,.188176800078,.527201188932,
                .688459039454,.281172343661,-.249846424327,
                -.195946274377,.127369340336,.093057364604,
                -.071394147166,-.029457536822,.033212674059,
                .003606553567,-.010733175483,.001395351747,
                .001992405295,-.000685856695,-.000116466855,
                .000093588670,-.000013264203]
        end

    else
        error("Wrong arguments, see comments for acceptable values")
    end

    f = f/norm(f) #list?

    if length(f) % 2 == 0
        f = [[0]; f]
    end

    return f
end


text:
module Learning

"""
    Optimizer
WIP - Currently LearningUpdate and Optimizer are haphazardly similar....
"""
abstract type Optimizer end
"""
    LearningUpdate
WIP - Currently LearningUpdate and Optimizer are haphazardly similar....
"""
abstract type LearningUpdate end

"""
    AbstractValueFunctionn
Abstract type definition for a value function
"""
abstract type AbstractValueFunction end

weights(q::AbstractValueFunction) = nothing
update!(value::AbstractValueFunction, Δθ) = throw("Implement weight update for $(typeof(value))")
update!(value::AbstractValueFunction, s, δ) = throw("Implement weight update for $(typeof(value))")

feature_type(value::AbstractValueFunction) = Float64

"""
    AbstractVFunction
Abstract type definition for a state value function
"""
abstract type AbstractVFunction <: AbstractValueFunction end

"""
    update!(value::ValueFunction, opt::Optimizer, ρ, s_t, s_tp1, reward, γ, terminal)

# Arguments:
`value::ValueFunction`:
`opt::Optimizer`:
`ρ`: Importance sampling ratios (Array of Floats)
`s_t`: States at time t
`s_tp1`: States at time t + 1
`reward`: cumulant or reward for value function
`γ`: discount factor
`terminal`: Determining termination of the episode (if applicable).
"""
function update!(value::AbstractVFunction, lu::LearningUpdate, ϕ_t, ϕ_tp1, reward, γ, ρ, terminal)
    throw(ErrorException("Implement update for $(typeof(opt)) and $(typeof(value))"))
end

"""
    AbstractQFunction
Abstract type definition for a state-action value function
"""
abstract type AbstractQFunction <: AbstractValueFunction end

update!(value::AbstractQFunction, s, a, δ) = throw("Implement weight update for $(typeof(value))")

"""
    update!(value::ValueFunction, opt::Optimizer, ρ, s_t, s_tp1, reward, γ, terminal)

# Arguments:
`value::ValueFunction`:
`opt::Optimizer`:
`ρ`: Importance sampling ratios (Array of Floats)
`s_t`: States at time t
`s_tp1`: States at time t + 1
`reward`: cumulant or reward for value function
`γ`: discount factor
`terminal`: Determining termination of the episode (if applicable).
`a_t`: Action at time t
`a_tp1`: Action at time t + 1
`target_policy`: Action at time t
"""
function update!(value::AbstractQFunction, lu::LearningUpdate, ϕ_t, ϕ_tp1, reward, γ, ρ, terminal, a_t, a_tp1, target_policy)
    throw(ErrorException("Implement update for $(typeof(opt)) and $(typeof(value))"))
end

# include("learning/Updates.jl")

export LinearRL
include("learning/LinearRL.jl")

export TabularRL
include("learning/TabularRL.jl")

# include("learning/FluxRL.jl")

end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
"""
Hard sphere fluid in N dimensions
"""
mutable struct HardSphereFluid{N,T}
	box::RectangularBox{N,T}
	balls::Vector{MovableBall{N,T}}
end

#
# function HardSphereFluid(box::RectangularBox{N,T}, balls::Vector{MovableBall{N,T}}) where {N,T}
# 	partner1, partner2, collision_type, min_collision_time = find_collision(balls, box)
# 	current_time = zero(T)
#
# 	return HardSphereFluid{N,T}(box, balls, current_time, min_collision_time, partner1, partner2, collision_type)
#
# end

function HardSphereFluid{N,T}(box::RectangularBox{N,T}, num_balls, r) where {N,T}

	balls = [MovableBall(zero(SVector{N,T}), zero(SVector{N,T}), r) for i in 1:num_balls]

	return HardSphereFluid{N,T}(box, balls)
end

HardSphereFluid(N, num_balls, r) = HardSphereFluid{N,Float64}(unit_hypercube(N, Float64), num_balls, r)

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
function update_p!(state::State, data::Data, prior::Prior)
  a, b = params(prior.p)
  state.p = rand(Beta(a + state.beta, b + 1 - state.beta))
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
#=
MIT License

Optimal Testing and Containment Strategies for Universities in Mexico amid COVID-19

Copyright © 2021 Test and Contain. Luis Benavides-Vázquez, Héctor Alonso Guzmán-Gutiérrez, Jakob Jonnerby, Philip Lazos, Edwin Lock, Francisco J. Marmolejo-Cossío, Ninad Rajgopal,and José Roberto Tello-Ayala. https://www.testandcontain.com/

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#
struct GroupIterator2
    options::Array{Int,1}  # list of possible group sizes
    k::Int  # number of categories
end
function Base.iterate(iter::GroupIterator2, state::Int=0)
    state >= length(iter.options)^iter.k && return nothing
    k = iter.k
    o = length(iter.options)
    alloc = zeros(Int, k)
    temp = 1
    for i in k:-1:1
        @inbounds alloc[i] = iter.options[state ÷ temp % o + 1]
        temp *= o
    end
    return alloc, state+1
end
Base.length(iter::GroupIterator2) = length(iter.options)^iter.k
Base.eltype(::Type{GroupIterator2}) = Array{Int,1}

struct GroupIterator
    options::Array{Int,1}
    k::Int  # number of categories
end
function Base.iterate(iter::GroupIterator)
    state = ones(Int, iter.k)
    item = map(i->iter.options[i], state)
    return item, (item, state)
end
function Base.iterate(iter::GroupIterator, state)
    state === nothing && return nothing
    _giterate(iter, state) === nothing && return nothing
    return copy(state[1]), state
end
function _giterate(iter::GroupIterator, next, i::Int=1)
    i > iter.k && return nothing
    item, state = next
    @inbounds state[i] += 1
    if state[i] > length(iter.options)  # carry over
        out = _giterate(iter, next, i+1)
        out === nothing && return nothing
        @inbounds state[i] = 1  # reset
    end
    @inbounds item[i] = iter.options[state[i]]
    return true
end
Base.length(iter::GroupIterator) = length(iter.options)^iter.k
Base.eltype(::Type{GroupIterator}) = Array{Int,1}
Base.firstindex(::GroupIterator) = 1

struct TestIterator
    budget::Int
    step::Int
    lower::Vector{Int}
    upper::Vector{Int}
    k::Int  # number of categories
end
function Base.iterate(iter::TestIterator)
    k = iter.k
    t = copy(iter.lower)
    t[k] = iter.budget - sum(t[1:end-1])
    while !(iter.lower[k] <= t[k] <= iter.upper[k])
        !_titerate(iter, t) && return nothing
    end
    return t, t
end
function Base.iterate(iter::TestIterator, t)
    t === nothing && return nothing
    k = iter.k
    !_titerate(iter, t) && return nothing
    while !(iter.lower[k] <= t[k] <= iter.upper[k])
        !_titerate(iter, t) && return nothing
    end
    return t, t
end
"""Iterates to next element. Returns false if spent, true otherwise."""
function _titerate(iter::TestIterator, t)
    # Invariant: t[k] always contains remaining budget
    k = iter.k
    i = k-1
    while i >= 1 && (t[i] >= iter.upper[i] || t[k] < 0)
        # Reset allocation to category i
        t[k] += t[i] - iter.lower[i]  # reclaim budget
        t[i] = iter.lower[i]
        i -= 1
    end
    i === 0 && return false
    t[i] += iter.step
    t[k] -= iter.step
    return true
end
Base.eltype(::Type{TestIterator}) = Vector{Int}
Base.firstindex(::TestIterator) = 1

max_stars_repo_name: hguzmang/eaamo

max_issues_repo_name: hguzmang/eaamo

max_forks_repo_name: hguzmang/eaamo