text:
import ComplexOptInterface
const COI = ComplexOptInterface

function matrix_cone(::Type{COI.HermitianPositiveSemidefiniteConeTriangle}, d)
    return COI.HermitianPositiveSemidefiniteConeTriangle(d)
end

vectorized_matrix(Q, basis, ::Type, ::Type) = MultivariateMoments.SymMatrix(Q, basis)
function vectorized_matrix(Q, basis, ::Type{COI.HermitianPositiveSemidefiniteConeTriangle}, T::Type)
    return MultivariateMoments.VectorizedHermitianMatrix{eltype(Q), T}(Q, basis)
end
function vectorized_matrix(Q, basis, ::Type{COI.HermitianPositiveSemidefiniteConeTriangle}, ::Type{Complex{T}}) where T
    return vectorized_matrix(Q, basis, COI.HermitianPositiveSemidefiniteConeTriangle, T)
end

# Need these two methods to avoid ambiguity
function build_gram_matrix(q::Vector{MOI.VariableIndex},
                           basis::AbstractPolynomialBasis, matrix_cone_type, T::Type)
    return build_gram_matrix([MOI.SingleVariable(vi) for vi in q], basis, matrix_cone_type, T)
end
#function build_gram_matrix(q::Vector{MOI.VariableIndex},
#                           basis::AbstractPolynomialBasis, T::Type, vectorization::HermitianVectorized)
#    return build_gram_matrix([MOI.SingleVariable(vi) for vi in q], basis, T, vectorization)
#end
function build_gram_matrix(q::Vector,
                           basis::AbstractPolynomialBasis, matrix_cone_type, T::Type)
    n = length(basis)
    Q = vectorized_matrix(q, length(basis), matrix_cone_type, T)
    U = _promote_sum(eltype(Q), T)
#    N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
#    if length(q) != N
#        throw(DimensionMismatch("Length of vectorization $(length(q)) does not match the side dimension $(n) of the symmetric matrix, it should be $N."))
#    end
    return GramMatrix{eltype(Q), typeof(basis), U}(Q, basis)
end
#function build_gram_matrix(q::Vector,
#                           basis::AbstractPolynomialBasis, T::Type, ::HermitianVectorized)
#    n = length(basis)
#    N = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
#    M = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n - 1))
#    U = _promote_sum(eltype(q), T)
#    if length(q) != N + M
#        throw(DimensionMismatch("Length of real vectorization $(length(q)) does not match the side dimension $(n) of the hermitian matrix, it should be $(N + M) for complex hermitian matrices."))
#    end
#    # Complex case, the first part of `q` is the real part,
#    # the second part is the complex part of the upper off-diagonal entries.
#    C = [zero(U) for i in 1:N]
#    k_real = 0
#    k_imag = 0
#    for j in 1:n
#        for i in 1:j
#            k_real += 1
#            C[k_real] = MA.operate!(MA.add_mul, C[k_real], one(T), q[k_real])
#            if i != j
#                k_imag += 1
#                C[k_real] = MA.operate!(MA.add_mul, C[k_real], one(T) * im, q[N + k_imag])
#            end
#        end
#    end
#    return build_gram_matrix(C, basis, T, SymmetricVectorized())
#end

function build_matrix(Q::Function, bases::Vector{<:AbstractPolynomialBasis}, f::Function)
    return map(eachindex(bases)) do i
        return f(Q(i), bases[i])
    end
end
function build_matrix(Q::Function, bases::Vector{<:Vector{<:AbstractPolynomialBasis}}, f::Function)
    return [
        f(Q(i), bases[i][j])
        for i in eachindex(bases) for j in eachindex(bases[i])
    ]
end

function build_gram_matrix(Q::Function, bases::Vector, matrix_cone_type, T)
    return SparseGramMatrix(build_matrix(Q, bases, (Q, b) -> build_gram_matrix(Q, b, matrix_cone_type, T)))
end

function build_moment_matrix(Q::Function, bases::Vector)
    return MultivariateMoments.SparseMomentMatrix(build_matrix(Q, bases, build_moment_matrix))
end

function union_constraint_indices_types(MCT)
    return Union{MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(matrix_cone(MCT, 0))},
                 MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(matrix_cone(MCT, 1))},
                 MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(matrix_cone(MCT, 2))},
                 MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(matrix_cone(MCT, 3))}}
end

function add_gram_matrix(model::MOI.ModelLike, matrix_cone_type::Type,
                         basis::AbstractPolynomialBasis, T::Type)
    Q, cQ = MOI.add_constrained_variables(model, matrix_cone(matrix_cone_type, length(basis)))
    q = build_gram_matrix(Q, basis, matrix_cone_type, T)
    return q, Q, cQ
end
_first(b::AbstractPolynomialBasis) = b
_first(b::Vector) = first(b)
function add_gram_matrix(model::MOI.ModelLike, matrix_cone_type::Type,
                         bases::Vector, T::Type)
    Qs = Vector{Vector{MOI.VariableIndex}}(undef, length(bases))
    cQs = Vector{union_constraint_indices_types(matrix_cone_type)}(undef, length(bases))
    cur = 0
    g = build_gram_matrix(bases, matrix_cone_type, T) do i
        # If `bases` is a vector of vector, the function
        # is called several times for the same `i` but we only need to
        # add it once.
        if i != cur
            @assert i == cur + 1
            cur = i
            Qs[i], cQs[i] = MOI.add_constrained_variables(model, matrix_cone(matrix_cone_type, length(_first(bases[i]))))
        end
        return Qs[i]
    end
    return g, Qs, cQs
end

function build_moment_matrix(q::Vector,
                             basis::AbstractPolynomialBasis)
    return MomentMatrix(MultivariateMoments.SymMatrix(q, length(basis)),
                        basis)
end

struct SOSPolynomialSet{DT <: AbstractSemialgebraicSet,
                        MT <: MP.AbstractMonomial,
                        MVT <: AbstractVector{MT},
                        CT <: Certificate.AbstractCertificate} <: MOI.AbstractVectorSet
    domain::DT
    monomials::MVT
    certificate::CT
end
MOI.dimension(set::SOSPolynomialSet) = length(set.monomials)
Base.copy(set::SOSPolynomialSet) = set


text:
struct SolutionCache{M <: MOI.ModelLike} <: MOI.ModelLike
    modattr::Dict{MOI.AbstractModelAttribute, Any}
    varattr::Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}
    conattr::Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}

    function SolutionCache(model::MOI.ModelLike)
        M = typeof(model)
        cache = new{M}(Dict{MOI.AbstractModelAttribute, Any}(),
                       Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}(),
                       Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}())
        return cache
    end

    function SolutionCache(model::MOI.ModelLike, variables::Vector{<:VI}, constraints::Vector{<:CI})
        M = typeof(model)
        cache = new{M}(Dict{MOI.AbstractModelAttribute, Any}(),
                       Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}(),
                       Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}())
        cache_model_attributes!(cache, model)
        cache_variable_attributes!(identity, cache, model, variables)
        cache_constraint_attributes!(identity, cache, model, constraints)
        return cache
    end
end

function cache_model_attributes!(cache::SolutionCache, src::MOI.ModelLike)
    attributes = [MOI.ObjectiveValue(),
                  MOI.DualObjectiveValue(),
                  MOI.ObjectiveBound(),
                  MOI.RelativeGap(),
                  MOI.SolveTime(),
                  MOI.SimplexIterations(),
                  MOI.BarrierIterations(),
                  MOI.NodeCount(),
                  MOI.TerminationStatus(),
                  MOI.RawStatusString(),
                  MOI.PrimalStatus(),
                  MOI.DualStatus()]
    for attr in attributes
        try
            value = MOI.get(src, attr)
            MOI.set(cache, attr, value)
        catch
        end
    end
    return nothing
end

function cache_model_attributes!(cache::SolutionCache, src::MOI.ModelLike, stage::Integer, scenario_index::Integer)
    attributes = [MOI.ObjectiveValue(),
                  MOI.DualObjectiveValue(),
                  MOI.ObjectiveBound(),
                  MOI.RelativeGap(),
                  MOI.SolveTime(),
                  MOI.SimplexIterations(),
                  MOI.BarrierIterations(),
                  MOI.NodeCount(),
                  MOI.TerminationStatus(),
                  MOI.RawStatusString(),
                  MOI.PrimalStatus(),
                  MOI.DualStatus()]
    for attr in attributes
        try
            value = MOI.get(src, ScenarioDependentModelAttribute(stage, scenario_index, attr))
            MOI.set(cache, attr, value)
        catch
        end
    end
    return nothing
end

function cache_variable_attributes!(index_map::Function, cache::SolutionCache, src::MOI.ModelLike, variables::Vector{<:VI})
    attr = MOI.VariablePrimal()
    for vi in variables
        try
            value = MOI.get(src, attr, index_map(vi))
            MOI.set(cache, attr, vi, value)
        catch
        end
    end
    return nothing
end
function cache_variable_attributes!(cache::SolutionCache, src::MOI.ModelLike)
    variables = MOI.get(src, MOI.ListOfVariableIndices())
    cache_variable_attributes!(identity, cache, src, variables)
    return nothing
end
function cache_variable_attributes!(cache::SolutionCache, src::MOI.ModelLike, variables::Vector{<:VI}, stage::Integer, scenario_index::Integer)
    attr = MOI.VariablePrimal()
    for vi in variables
        try
            value = MOI.get(src, ScenarioDependentVariableAttribute(stage, scenario_index, attr), vi)
            MOI.set(cache, attr, vi, value)
        catch
        end
    end
    return nothing
end
function cache_variable_attributes!(cache::SolutionCache, src::MOI.ModelLike, stage::Integer, scenario_index::Integer)
    variables = MOI.get(src, ScenarioDependentModelAttribute(stage, scenario_index, MOI.ListOfVariableIndices()))
    cache_variable_attributes!(cache, src, variables, stage, scenario_index)
    return nothing
end

function cache_constraint_attributes!(index_map::Function, cache::SolutionCache, src::MOI.ModelLike, constraints::Vector{<:CI})
    conattrs = [MOI.ConstraintPrimal(), MOI.ConstraintDual(), MOI.ConstraintBasisStatus()]
    for ci in constraints
        for attr in conattrs
            try
                value = MOI.get(src, attr, index_map(ci))
                MOI.set(cache, attr, ci, value)
            catch
            end
        end
    end
    return nothing
end
function cache_constraint_attributes!(cache::SolutionCache, src::MOI.ModelLike)
    ctypes = filter(t -> is_decision_type(t[1]), MOI.get(src, MOI.ListOfConstraints()))
    constraints = mapreduce(vcat, ctypes) do (F, S)
        return MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
    end
    cache_constraint_attributes!(identity, cache, src, constraints)
    return nothing
end
function cache_constraint_attributes!(cache::SolutionCache, src::MOI.ModelLike, constraints::Vector{<:CI}, stage::Integer, scenario_index::Integer)
    conattrs = [MOI.ConstraintPrimal(), MOI.ConstraintDual(), MOI.ConstraintBasisStatus()]
    for ci in constraints
        for attr in conattrs
            try
                value = MOI.get(src, ScenarioDependentConstraintAttribute(stage, scenario_index, attr), ci)
                MOI.set(cache, attr, ci, value)
            catch
            end
        end
    end
    return nothing
end
function cache_constraint_attributes!(cache::SolutionCache, src::MOI.ModelLike, stage::Integer, scenario_index::Integer)
    ctypes = MOI.get(src, ScenarioDependentModelAttribute(stage, scenario_index, MOI.ListOfConstraints()))
    constraints = mapreduce(vcat, filter(t -> is_decision_type(t[1]), ctypes)) do (F, S)
        return MOI.get(src, ScenarioDependentModelAttribute(stage, scenario_index, MOI.ListOfConstraintIndices{F,S}()))
    end
    cache_constraint_attributes!(cache, src, constraints, stage, scenario_index)
    return nothing
end

function MOI.get(cache::SolutionCache{M},
                 attr::MOI.AbstractModelAttribute) where M
    if !haskey(cache.modattr, attr)
        throw(ArgumentError("Origin model type $M did not support accessing the attribute $attr"))
    end
    return cache.modattr[attr]
end

function MOI.get(cache::SolutionCache{M},
                 attr::MOI.AbstractVariableAttribute,
                 vi::VI) where M
    attribute_dict = get(cache.varattr, attr, nothing)
    if attribute_dict === nothing
        throw(ArgumentError("Origin model type $M did not support accessing the attribute $attr"))
    end
    if !haskey(attribute_dict, vi)
        throw(ArgumentError("Origin model type $M did not support accessing the attribute $attr"))
    end
    return get(attribute_dict, vi, nothing)
end

function MOI.get(cache::SolutionCache{M},
                 attr::MOI.AbstractConstraintAttribute,
                 ci::CI) where M
    attribute_dict = get(cache.conattr, attr, nothing)
    if attribute_dict === nothing
        throw(ArgumentError("Origin model type $M did not support accessing the attribute $attr"))
    end
    if !haskey(attribute_dict, ci)
        throw(ArgumentError("Origin model type $M did not support accessing the attribute $attr"))
    end
    return get(attribute_dict, ci, nothing)
end

function MOI.set(cache::SolutionCache,
                 attr::MOI.AbstractModelAttribute,
                 value)
    cache.modattr[attr] = value
    return nothing
end

function MOI.set(cache::SolutionCache,
                 attr::MOI.AbstractVariableAttribute,
                 vi::VI,
                 value)
    if !haskey(cache.varattr, attr)
        cache.varattr[attr] = Dict{VI, Any}()
    end
    cache.varattr[attr][vi] = value
    return nothing
end

function MOI.set(cache::SolutionCache,
                 attr::MOI.AbstractConstraintAttribute,
                 ci::CI,
                 value)
    if !haskey(cache.conattr, attr)
        cache.conattr[attr] = Dict{VI, Any}()
    end
    cache.conattr[attr][ci] = value
    return nothing
end

MOI.supports(::SolutionCache, ::MOI.AbstractModelAttribute) = true
MOI.supports(::SolutionCache, ::Union{MOI.AbstractVariableAttribute, MOI.AbstractConstraintAttribute}, ::Type{<:MOI.Index}) = true

ext: jl


text:
using LinearAlgebra
using Distributions

"""
    multinomial_hyper(α::AbstractArray{Float32,1})
[Dirichlet Distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution)
"""
struct multinomial_dist <: distibution_sample
    α::AbstractArray{Float32,1}
end

function log_likelihood!(r::AbstractArray,x::AbstractArray, distibution_sample::multinomial_dist , group::Int64 = -1)
    r .= (distibution_sample.α' * x)[1,:]
end


text:
# star.jl : read star data

const STAR_P = Matrix{Float64}([
-0.90097 -0.43388
-0.90097 0.43388
-0.86451 -0.39443
-0.86051 0.3901
-0.84371 -0.42971
-0.84321 0.42967
-0.83242 -0.3597
-0.82 0.34626
-0.80909 0.38907
-0.79777 -0.32222
-0.79647 -0.42626
-0.78564 0.30909
-0.78504 0.42543
-0.78119 -0.37694
-0.76462 0.36736
-0.76098 -0.2824
-0.74803 -0.33296
-0.74789 0.26824
-0.74516 -0.42252
-0.73603 0.3214
-0.73588 0.42185
-0.7251 -0.24357
-0.72388 -0.3783
-0.71338 0.37276
-0.70751 -0.29252
-0.70711 0.22412
-0.7019 0.27417
-0.69493 -0.33601
-0.69103 -0.20671
-0.68995 -0.4185
-0.68588 0.32592
-0.68306 0.41799
-0.67067 -0.26099
-0.67052 0.18452
-0.66532 -0.37173
-0.65929 0.37056
-0.65747 0.27921
-0.65548 -0.16825
-0.64934 0.23178
-0.64892 -0.31472
-0.63935 -0.22039
-0.63575 -0.41454
-0.6342 0.14522
-0.6327 0.32338
-0.62951 0.41409
-0.61904 -0.12881
-0.61827 -0.27111
-0.61587 0.19149
-0.61573 -0.36309
-0.60899 0.27269
-0.60537 0.36672
-0.60508 -0.17994
-0.59709 0.10506
-0.5928 -0.31758
-0.58773 0.22856
-0.58637 -0.22848
-0.58409 -0.41078
-0.58363 0.1492
-0.58166 -0.088365
-0.58099 0.31782
-0.57615 0.4102
-0.57004 -0.13886
-0.56538 0.070752
-0.56303 -0.27407
-0.56217 -0.36309
-0.55741 0.18381
-0.55577 0.27079
-0.55277 -0.18618
-0.55266 0.36222
-0.54288 -0.0464
-0.54058 0.1249
-0.53676 -0.31795
-0.53485 -0.096479
-0.53262 0.035302
-0.53164 -0.23095
-0.5314 0.22464
-0.53051 -0.40687
-0.52792 0.31485
-0.52333 0.40634
-0.51889 -0.14291
-0.51102 0.08396
-0.50924 0.1724
-0.50734 -0.27411
-0.5072 -0.36071
-0.50253 0.26766
-0.50173 -0.049168
-0.5 6.1232e-17
-0.49973 0.35886
-0.49927 -0.18727
-0.48658 -0.096198
-0.48413 0.12923
-0.4807 -0.3161
-0.47988 0.21857
-0.47758 0.045037
-0.47671 -0.23021
-0.47618 -0.40291
-0.47496 0.3116
-0.47046 0.40249
-0.46772 -0.14124
-0.45927 -0.044641
-0.45451 0.174
-0.45375 0.088745
-0.45196 -0.35764
-0.45172 -0.27206
-0.45062 0.2637
-0.44684 0.35529
-0.44669 0.0016115
-0.44615 -0.18499
-0.4386 -0.09242
-0.42706 0.13181
-0.42561 0.21846
-0.42503 -0.31329
-0.42304 0.045802
-0.42241 0.30807
-0.42221 -0.2276
-0.42156 -0.39892
-0.41765 0.39864
-0.41682 -0.13817
-0.41227 -0.04467
-0.39878 0.17514
-0.39793 0.089063
-0.39694 0.26251
-0.3969 -0.35431
-0.39651 -0.26935
-0.39412 0.35173
-0.39366 0.001052
-0.39317 -0.18244
-0.3886 -0.091603
-0.3712 0.13202
-0.37057 0.21877
-0.36967 -0.31059
-0.36898 0.0449
-0.36884 0.30604
-0.3678 -0.22567
-0.36684 -0.39493
-0.36489 0.39479
-0.36487 -0.13701
-0.36285 -0.045128
-0.3436 0.1753
-0.34315 0.088222
-0.34228 0.26249
-0.34186 -0.3514
-0.34117 -0.26838
-0.34081 0.34897
-0.34068 -0.00013538
-0.33987 -0.1813
-0.33757 -0.091348
-0.31642 0.13155
-0.31585 0.21896
-0.3153 0.043673
-0.31408 0.30628
-0.31368 -0.22498
-0.31365 -0.31104
-0.3125 -0.13653
-0.31188 -0.045899
-0.31174 -0.39092
-0.31174 0.39092
-0.30442 0.43887
-0.30393 -0.44209
-0.29743 0.48462
-0.29655 -0.49034
-0.28935 0.53747
-0.28921 0.17521
-0.28916 0.087303
-0.28812 -0.54553
-0.28792 0.26298
-0.28778 -0.0013406
-0.28669 -0.26862
-0.28663 -0.18102
-0.28584 -0.091567
-0.28572 0.35044
-0.285 -0.35404
-0.28158 0.58835
-0.28038 -0.59624
-0.27409 0.63736
-0.27238 -0.64859
-0.26531 0.69482
-0.26446 -0.40743
-0.26368 -0.70555
-0.26338 0.40202
-0.2625 0.13109
-0.26202 0.042633
-0.2618 0.21931
-0.26011 0.30774
-0.2601 -0.046717
-0.26008 -0.22536
-0.25998 -0.13664
-0.25928 -0.31295
-0.25528 -0.76052
-0.25446 -0.51103
-0.25408 0.5011
-0.25388 0.7697
-0.24941 -0.45752
-0.2481 0.45088
-0.24713 -0.81383
-0.24581 0.82249
-0.24032 0.55182
-0.23919 -0.86583
-0.23815 0.87264
-0.23697 -0.56078
-0.23609 -0.36312
-0.23579 0.086624
-0.23557 0.17524
-0.23553 0.35662
-0.23504 -0.002297
-0.23441 0.26415
-0.23376 -0.092042
-0.23374 -0.18147
-0.23333 -0.27031
-0.23325 0.73325
-0.23252 0.65288
-0.23187 -0.91373
-0.23135 -0.61057
-0.23101 0.91936
-0.22887 -0.67028
-0.22531 0.60088
-0.22252 -0.97493
-0.22252 0.97493
-0.21581 -0.72508
-0.21535 -0.41209
-0.21367 0.40522
-0.21344 -0.50208
-0.21123 0.49821
-0.20927 0.13084
-0.20913 0.041902
-0.20884 -0.31844
-0.20874 0.31136
-0.20861 0.22001
-0.20801 -0.047429
-0.20761 -0.13716
-0.20736 -0.22668
-0.20668 0.78864
-0.20535 0.69446
-0.20304 -0.77766
-0.19828 -0.85007
-0.19566 -0.93569
-0.19378 -0.45624
-0.19335 0.54695
-0.19329 -0.54705
-0.19299 -0.63772
-0.19277 0.45063
-0.19265 0.84414
-0.19134 0.92939
-0.18583 -0.36612
-0.18479 0.64517
-0.1847 0.35875
-0.18293 0.086203
-0.18267 0.1755
-0.18251 -0.0029727
-0.1825 0.26618
-0.18219 -0.27359
-0.18166 -0.092601
-0.18144 -0.18245
-0.18083 0.74443
-0.17865 -0.68578
-0.17712 -0.59117
-0.17285 0.59105
-0.17014 -0.81263
-0.16952 -0.89752
-0.16476 -0.73334
-0.16447 0.89015
-0.16393 -0.50114
-0.1631 0.4968
-0.16304 -0.41166
-0.16172 0.40501
-0.1584 0.7994
-0.15796 -0.32039
-0.15734 0.31272
-0.15662 0.13077
-0.15656 0.041421
-0.15648 0.22102
-0.15592 -0.22862
-0.15591 -0.047964
-0.15557 -0.13786
-0.15483 0.69682
-0.14568 -0.63798
-0.14194 0.54211
-0.14192 -0.54781
-0.14166 -0.85683
-0.14128 -0.77297
-0.14108 0.63205
-0.13926 -0.45594
-0.13846 0.4501
-0.13543 0.84774
-0.13418 -0.3662
-0.13304 0.3589
-0.13183 0.75321
-0.13099 -0.2748
-0.13081 0.26691
-0.13049 0.17585
-0.13047 0.085931
-0.13019 -0.0033918
-0.12987 -0.18356
-0.12973 -0.093085
-0.12882 -0.68588
-0.12138 -0.59269
-0.11781 0.58427
-0.11549 0.6697
-0.11364 -0.5016
-0.1132 -0.81527
-0.11156 0.49489
-0.1103 -0.41122
-0.10925 0.4046
-0.10889 -0.72971
-0.10647 -0.32053
-0.10571 0.31275
-0.10534 0.8038
-0.10466 0.22121
-0.10455 -0.22931
-0.10444 0.13069
-0.10422 0.041171
-0.10392 -0.13842
-0.1039 -0.048317
-0.10016 0.71561
-0.095607 -0.63891
-0.089738 -0.54773
-0.087782 0.53858
-0.087699 0.62453
-0.085715 -0.45578
-0.084288 0.44933
-0.0841 -0.77277
-0.082161 -0.3658
-0.081172 0.35844
-0.079537 -0.27486
-0.079153 0.2666
-0.078577 0.17557
-0.078407 -0.18374
-0.078215 0.08586
-0.078032 -0.0036225
-0.077956 -0.093385
-0.076428 0.76157
-0.075891 -0.68438
-0.067551 -0.59284
-0.065297 0.67144
-0.063239 -0.50133
-0.061583 0.58009
-0.057965 0.49334
-0.0576 -0.41064
-0.056719 0.40375
-0.054748 -0.32021
-0.054084 0.31192
-0.053766 -0.72847
-0.053162 -0.22889
-0.052935 0.22033
-0.052269 -0.13843
-0.052248 0.13064
-0.052056 0.041098
-0.051973 -0.04852
-0.046193 0.71741
-0.043757 -0.63832
-0.037033 -0.54741
-0.035993 0.62533
-0.032291 -0.45497
-0.032165 0.53626
-0.030027 0.44817
-0.02997 -0.36538
-0.029393 0.35704
-0.028162 -0.2738
-0.027632 0.26471
-0.026778 -0.18334
-0.026327 0.17539
-0.026123 -0.093441
-0.026048 0.085884
-0.026009 -0.0037078
-0.022437 -0.68272
-0.014591 0.67126
-0.014242 -0.59313
-0.014037 -0.50051
-0.0072757 0.58106
-0.0046829 -0.4105
-0.0046299 0.40196
-0.0034727 0.49154
-0.0033654 -0.31831
-0.0029089 0.30833
-0.0014897 -0.22801
-0.00040952 -0.1383
-0.0004079 0.21995
-6.0753e-05 -0.048552
-7.7757e-06 0.1307
1.6295e-06 0.041118
0.009434 -0.63617
0.01588 -0.55032
0.016596 0.62571
0.020993 0.35079
0.021249 -0.36227
0.021936 -0.45567
0.02211 0.53798
0.023678 -0.27229
0.024397 0.44559
0.025221 -0.18305
0.025465 0.26411
0.025783 -0.093389
0.025956 -0.003685
0.026059 0.085988
0.026133 0.17549
0.030418 -0.50773
0.041342 -0.58957
0.043864 0.39161
0.046026 -0.4054
0.047708 0.58028
0.048706 0.49581
0.04874 -0.3159
0.050804 -0.22763
0.051056 0.30769
0.051577 -0.13819
0.051863 -0.048508
0.051998 0.041203
0.052218 0.1309
0.052529 0.22012
0.068437 -0.4902
0.073579 -0.35847
0.073981 -0.54191
0.074378 0.42014
0.074718 -0.44356
0.076152 0.35076
0.076588 -0.27192
0.076704 0.46184
0.077391 -0.18292
0.077741 -0.093325
0.077899 -0.0036247
0.078102 0.086145
0.078611 0.17582
0.079154 0.26442
0.079351 0.53406
0.10022 -0.39761
0.1027 -0.31568
0.10326 -0.22756
0.10363 -0.13812
0.10381 -0.048464
0.10395 0.041289
0.10435 0.13116
0.10552 0.30836
0.10555 0.22067
0.10619 0.38702
0.11126 -0.48746
0.11126 0.48746
0.12631 -0.43772
0.12668 0.43406
0.12961 -0.18289
0.1297 -0.27201
0.12974 -0.093294
0.12984 -0.0036
0.12996 -0.35761
0.13007 0.086268
0.13094 0.17628
0.13299 0.26521
0.13454 0.34881
0.15567 -0.22765
0.15572 -0.13811
0.1558 -0.048482
0.15588 0.041283
0.15633 0.13137
0.15721 -0.40006
0.15737 -0.31546
0.15778 0.46918
0.15844 0.22141
0.15897 -0.51488
0.15919 -0.46995
0.15936 0.51511
0.15971 0.39543
0.16154 0.30783
0.18181 -0.093335
0.18183 -0.003708
0.18185 -0.18291
0.18195 0.086236
0.18294 0.17668
0.18297 -0.2721
0.1851 -0.359
0.18764 0.26544
0.18778 0.43678
0.18819 -0.43823
0.18824 0.35365
0.20344 -0.54044
0.20398 0.54075
0.20659 0.48958
0.20696 -0.49083
0.20783 0.041052
0.20788 -0.13814
0.20788 -0.048658
0.20796 -0.22771
0.20809 0.13136
0.21133 0.2218
0.21149 -0.31652
0.21467 -0.40095
0.21528 0.3107
0.21563 0.39831
0.23377 0.085875
0.23391 -0.0040934
0.23403 -0.093519
0.23416 -0.18283
0.23456 0.17672
0.23564 -0.27294
0.2377 -0.45236
0.23805 0.45033
0.24026 -0.36015
0.24133 0.26726
0.24227 0.35594
0.244 -0.56375
0.24473 0.56417
0.25337 -0.50665
0.25381 0.50523
0.25959 -0.22797
0.25966 0.13082
0.25983 0.040411
0.26015 -0.049131
0.26024 -0.13815
0.2631 0.22287
0.26517 -0.40919
0.26559 -0.31784
0.26578 0.40561
0.26912 0.31201
0.2851 0.17665
0.28559 0.084787
0.28571 -0.4598
0.28616 -0.0049461
0.28653 -0.18277
0.2866 -0.093714
0.28664 0.45681
0.28748 -0.27389
0.2884 -0.54183
0.28846 0.54119
0.28964 -0.58997
0.29094 0.59072
0.29204 -0.36508
0.29362 0.35887
0.29549 0.26811
0.31033 -0.2283
0.31076 0.12958
0.311 -0.50228
0.31201 0.038627
0.31214 0.50018
0.31282 -0.13809
0.31288 -0.049423
0.31379 0.2232
0.31507 0.40688
0.31545 -0.41294
0.31814 -0.32086
0.32239 0.31264
0.33051 -0.61346
0.33207 0.61436
0.3344 0.17551
0.33637 -0.55602
0.33716 0.082063
0.33719 0.55324
0.33754 -0.27605
0.33889 -0.45788
0.33906 -0.18266
0.3393 -0.0062857
0.33943 0.45187
0.33951 -0.093042
0.34411 -0.36656
0.34484 0.35812
0.3505 0.26774
0.35927 -0.2292
0.35996 -0.50815
0.361 0.12586
0.36211 0.4981
0.36249 0.2223
0.36497 0.036381
0.36517 -0.13624
0.36621 0.40462
0.36685 -0.049946
0.36781 -0.41073
0.37191 -0.32199
0.37424 -0.59139
0.37451 0.58857
0.37624 0.31017
0.37639 -0.63983
0.37834 0.64095
0.38138 0.17035
0.38633 -0.27703
0.38812 -0.45858
0.38908 0.53488
0.38947 0.079671
0.39066 -0.17912
0.39272 0.45156
0.3929 -0.55149
0.39354 -0.0078401
0.39389 -0.092765
0.3957 0.3566
0.39825 -0.36427
0.40414 -0.22276
0.40663 0.21532
0.40755 0.26274
0.40999 -0.50491
0.41241 0.12291
0.41629 0.40545
0.41791 -0.66369
0.41852 0.035961
0.41856 -0.41178
0.41917 -0.13542
0.42029 0.49168
0.42049 0.56891
0.42211 -0.051552
0.42233 -0.60669
0.42416 0.61745
0.42503 0.66778
0.42837 0.30823
0.42843 -0.31834
0.42903 -0.26547
0.42933 0.1674
0.43909 -0.46004
0.44105 -0.55522
0.44325 0.07943
0.44474 -0.17693
0.44476 0.35705
0.44743 -0.0063528
0.4494 0.52837
0.44989 -0.096828
0.45017 -0.36762
0.45048 -0.21694
0.45048 0.21694
0.45098 0.45008
0.45746 -0.51029
0.45938 -0.68752
0.4648 0.2637
0.46501 0.12265
0.46502 0.40598
0.46725 0.58084
0.46739 0.69212
0.4682 -0.63389
0.46872 -0.2765
0.46943 -0.41992
0.4731 0.039431
0.47328 0.63303
0.47457 -0.58663
0.47492 -0.050802
0.47492 0.16872
0.47591 -0.14395
0.47953 0.31179
0.48182 0.48582
0.48424 -0.32717
0.48919 -0.47473
0.4941 0.35935
0.49848 -0.53425
0.49908 0.085807
0.49991 0.53645
0.50022 -0.0047965
0.50058 -0.38051
0.50197 -0.095032
0.5022 -0.19653
0.50484 -0.71365
0.50517 -0.6667
0.50642 0.66989
0.50658 0.1948
0.50897 0.4079
0.51369 0.58833
0.51646 0.72032
0.51735 0.13658
0.51741 -0.43547
0.51859 -0.58406
0.52209 0.45074
0.52229 -0.14157
0.52562 0.040312
0.52804 -0.047559
0.52986 -0.63206
0.53186 0.63371
0.53478 -0.49217
0.53653 0.49789
0.54707 -0.73792
0.54997 0.087741
0.55077 -0.54439
0.55151 0.54681
0.55293 0.68943
0.55333 -0.68708
0.55383 -0.090598
0.5545 -0.0031196
0.55463 0.17582
0.55477 -0.17577
0.56336 0.74728
0.5662 0.59477
0.56741 -0.59871
0.568 0.1289
0.57323 -0.1287
0.57751 0.047381
0.5815 0.64473
0.58307 -0.050436
0.58365 -0.65176
0.5863 -0.76046
0.59702 0.6954
0.60147 -0.70994
0.60429 -0.00060389
0.60766 0.15489
0.60832 0.098195
0.60994 -0.15399
0.61055 0.73957
0.61463 -0.096944
0.62349 -0.78183
0.62349 0.78183
0.62485 0.049113
0.63357 -0.04583
0.65437 0.13645
0.65517 0.0032417
0.65801 -0.13501
0.66114 0.090053
0.66771 -0.082749
0.67516 0.045181
0.68649 -0.030527
0.70693 0.1157
0.70879 -0.11497
0.71456 0.018926
0.71778 0.070637
0.72067 -0.065609
0.73873 -0.016611
0.76121 -0.09427
0.76204 0.093945
0.76516 0.042919
0.77593 -0.047034
0.78493 0.0029256
0.81349 0.073631
0.81492 -0.073069
0.82671 0.028963
0.8287 -0.022385
0.86619 -0.052826
0.87162 0.050683
0.87833 0.001808
0.91453 -0.033742
0.92781 0.028498
0.95816 -0.016518
1 -2.4493e-16
])

const STAR_T = Matrix{Int64}([
44 36 31
2 4 6
32 36 45
657 676 664
664 644 657
689 670 686
1 5 3
9 4 8
13 6 9
9 6 4
18 20 12
24 20 31
31 36 24
24 36 32
32 21 24
680 664 676
670 664 680
680 686 670
720 718 719
421 407 430
451 465 444
444 430 451
444 452 431
465 452 444
431 421 444
421 430 444
165 190 161
200 190 165
311 270 292
292 329 311
294 330 313
329 292 313
348 330 362
329 313 348
348 313 330
688 681 679
615 642 643
665 681 678
665 643 642
665 678 656
656 643 665
5 11 7
7 3 5
15 12 20
20 24 15
15 9 8
8 12 15
15 24 21
15 21 13
13 9 15
19 30 23
241 265 283
283 265 303
173 197 216
257 216 238
238 216 197
192 177 210
628 645 612
649 645 628
713 711 714
671 676 657
657 647 671
717 719 718
599 622 600
329 348 364
512 503 528
466 452 481
481 452 465
408 421 431
431 422 408
291 270 311
174 200 165
217 236 212
212 236 259
37 44 31
203 181 224
224 181 202
221 180 204
221 265 241
227 204 184
660 665 642
660 679 681
681 665 660
615 643 621
621 643 656
656 651 621
606 618 587
587 618 613
487 481 465
503 481 487
556 574 586
586 570 556
530 504 515
14 17 10
10 7 14
14 7 11
14 23 17
14 11 19
19 23 14
17 23 28
173 216 175
162 160 191
162 197 173
191 197 162
191 160 194
194 221 241
180 221 194
223 197 191
223 238 197
223 194 241
191 194 223
277 297 257
257 238 277
648 664 670
644 664 648
648 620 644
214 243 218
214 261 243
214 199 261
284 266 307
287 266 254
331 307 287
287 307 266
261 199 242
284 261 242
242 266 284
199 196 242
245 216 257
606 598 630
630 598 616
705 711 708
708 713 709
711 713 708
646 652 613
613 618 646
712 707 709
712 710 707
715 710 712
709 713 712
613 652 631
663 671 647
647 637 663
715 712 716
716 712 713
718 715 716
716 717 718
716 713 714
714 717 716
615 621 588
588 621 595
614 603 635
635 651 662
622 599 632
534 553 526
558 553 534
503 512 491
491 481 503
466 481 491
491 482 466
394 407 421
421 408 394
240 255 276
57 77 65
27 37 31
31 20 27
27 20 18
96 116 103
84 77 96
96 103 84
84 103 92
84 65 77
181 203 163
224 202 247
247 269 224
270 291 247
247 291 269
393 364 378
378 364 348
378 348 362
362 392 378
393 378 406
406 378 392
420 430 407
420 407 393
393 406 420
249 292 270
167 205 182
362 330 345
141 151 133
206 227 184
204 227 246
246 221 204
265 221 246
626 621 651
635 603 626
626 651 635
626 603 595
595 621 626
435 437 414
558 570 583
526 553 538
538 553 567
511 496 526
485 496 511
526 538 511
511 538 525
552 538 567
525 538 552
513 487 504
504 530 513
528 503 513
503 487 513
574 556 543
543 513 530
543 556 528
528 513 543
649 628 617
617 636 649
35 23 30
35 28 23
35 30 42
98 79 88
88 106 98
241 283 263
263 223 241
238 223 263
263 277 238
331 287 314
349 331 314
266 242 232
232 242 196
232 196 192
254 266 232
192 210 232
232 210 254
254 210 233
233 210 177
211 175 216
216 245 211
245 233 211
177 175 211
211 233 177
257 297 281
281 245 257
684 690 695
675 690 684
682 690 675
687 694 696
629 605 641
668 645 649
649 673 668
644 620 625
657 644 625
625 647 657
593 596 566
597 625 620
596 625 597
566 596 569
569 540 566
596 597 569
634 618 606
634 646 618
606 630 634
659 663 637
659 631 652
637 631 659
595 603 577
577 565 595
555 542 577
577 542 565
382 351 367
397 382 367
412 382 397
595 565 568
568 588 595
696 694 700
672 662 651
656 678 672
672 651 656
685 687 674
694 687 685
610 632 599
526 496 507
507 534 526
512 534 507
507 496 482
482 491 507
507 491 512
541 534 512
541 512 528
541 570 558
558 534 541
528 556 541
541 556 570
380 364 393
393 407 380
407 394 380
379 394 408
346 379 361
359 325 344
256 240 276
239 200 256
219 179 189
176 179 215
255 240 215
179 219 215
215 219 255
47 41 33
10 17 16
29 33 41
142 135 156
44 60 51
51 61 45
36 44 51
51 45 36
44 37 50
50 60 44
114 106 97
97 88 78
106 88 97
101 111 93
116 135 123
123 103 116
135 142 123
120 111 101
248 203 224
224 269 248
225 247 202
202 182 225
225 182 205
205 249 225
270 247 225
225 249 270
590 610 581
581 610 599
571 599 600
492 480 467
467 480 448
502 527 516
502 480 492
492 519 502
502 519 527
428 441 419
419 405 428
419 406 392
392 405 419
430 420 442
451 430 442
442 463 451
504 487 477
451 463 477
477 465 451
477 487 465
515 504 488
463 479 488
504 477 488
488 477 463
292 249 273
294 313 273
273 313 292
164 202 181
164 182 202
362 345 376
376 392 362
376 405 392
390 405 376
312 330 294
312 345 330
439 427 448
426 443 454
358 375 343
390 375 403
345 312 327
360 375 390
390 376 360
360 376 345
345 327 360
343 375 360
360 327 343
141 133 122
122 133 114
180 194 158
158 194 160
98 106 117
204 180 171
184 204 171
171 151 184
114 133 125
136 117 125
125 106 114
125 117 106
183 163 203
227 206 250
250 271 289
268 246 227
268 250 289
227 250 268
195 198 235
259 279 235
235 212 259
198 212 235
410 412 435
435 414 410
425 414 437
601 583 570
601 570 586
601 628 612
612 583 601
586 617 601
601 617 628
567 553 584
584 598 567
584 553 558
558 583 584
616 598 584
612 616 584
584 583 612
341 325 359
606 587 580
580 552 567
580 598 606
567 598 580
576 587 613
525 552 535
54 64 47
69 88 79
69 51 60
69 60 78
78 88 69
69 79 61
61 51 69
254 233 275
298 314 275
275 233 245
275 287 254
275 314 287
275 281 298
245 281 275
684 673 661
661 675 684
649 636 661
661 673 649
701 707 702
696 698 691
691 687 696
705 708 703
703 700 705
703 698 696
696 700 703
629 641 653
611 641 605
600 622 611
654 661 636
675 661 654
636 624 654
654 624 638
677 691 682
687 691 677
699 695 702
619 625 596
596 593 619
647 625 619
619 637 647
613 631 592
592 576 613
557 576 592
683 673 684
683 668 673
603 614 585
585 577 603
555 577 585
614 602 585
381 397 367
381 350 365
367 350 381
332 342 365
365 350 332
332 350 315
539 568 565
539 542 522
565 542 539
590 602 623
623 610 590
633 602 614
633 614 635
650 623 633
633 623 602
572 602 590
572 585 602
572 546 555
555 585 572
650 633 658
635 662 658
658 633 635
669 653 641
669 685 674
674 653 669
363 380 394
394 379 363
379 346 363
308 344 325
308 325 289
289 271 308
315 350 333
367 351 333
333 350 367
213 256 200
240 256 213
200 174 213
213 174 176
176 215 213
213 215 240
189 195 234
234 219 189
195 235 234
65 84 72
72 54 65
64 54 72
72 84 92
56 41 47
47 64 56
64 75 56
127 108 115
37 27 39
39 50 37
55 50 39
17 28 25
25 16 17
41 56 52
46 52 62
101 93 82
105 93 111
114 97 105
105 122 114
111 122 105
76 55 66
66 82 76
76 82 93
131 123 142
92 103 112
103 123 112
123 131 112
119 138 126
145 155 167
145 138 155
132 126 145
145 126 138
228 248 271
271 250 228
228 250 206
203 248 228
228 183 203
206 183 228
508 514 536
536 514 546
524 551 537
524 495 508
508 536 524
524 536 551
590 581 563
581 551 563
563 572 590
551 536 563
546 572 563
563 536 546
564 581 599
599 571 564
564 551 581
537 551 564
454 443 468
537 519 509
509 524 537
495 524 509
509 519 492
545 519 537
537 564 545
545 564 571
571 554 545
527 519 545
545 554 527
448 480 464
464 439 448
480 502 490
490 464 480
478 464 490
490 502 516
429 419 441
429 442 420
429 420 406
406 419 429
463 442 450
441 462 450
450 429 441
442 429 450
450 479 463
450 462 479
274 312 294
172 142 156
231 272 253
448 427 440
440 426 454
440 467 448
454 467 440
388 403 375
375 358 388
402 426 416
416 388 402
403 388 416
427 403 416
416 440 427
426 440 416
418 403 427
428 405 418
418 405 390
390 403 418
418 439 428
427 439 418
343 327 309
272 288 309
253 272 293
293 327 312
272 309 293
293 309 327
293 274 253
312 274 293
151 171 144
144 133 151
144 125 133
136 125 144
166 151 141
184 151 166
166 206 184
166 183 206
246 268 286
303 265 286
265 246 286
396 382 412
412 410 396
395 408 422
395 379 408
361 379 395
423 391 409
422 433 409
409 433 423
409 395 422
445 452 466
431 452 445
445 422 431
445 433 422
466 482 456
456 445 466
433 445 456
549 576 557
65 54 49
49 35 42
42 57 49
57 65 49
698 701 693
693 690 682
682 691 693
693 691 698
709 707 704
707 701 704
704 708 709
704 703 708
704 701 698
698 703 704
629 653 640
636 617 604
604 624 636
604 589 624
574 589 604
586 574 604
604 617 586
667 677 682
667 654 638
682 675 667
675 654 667
706 699 702
702 707 706
707 710 706
557 592 573
573 593 566
637 619 608
608 619 593
608 631 637
608 592 631
593 573 608
608 573 592
695 699 692
684 695 692
692 683 684
276 255 295
295 315 276
295 332 315
632 610 639
610 623 639
639 623 650
328 363 346
328 291 311
290 308 271
290 248 269
271 248 290
276 315 296
315 333 296
296 256 276
321 342 304
342 332 304
332 295 304
258 235 279
258 234 235
258 279 300
115 108 95
75 95 89
108 99 89
89 95 108
146 127 134
143 152 134
134 152 146
134 127 115
118 99 108
108 127 118
186 231 208
208 231 253
185 205 167
167 155 185
71 82 66
16 25 22
33 29 22
22 25 33
52 46 38
38 29 41
41 52 38
67 50 55
55 76 67
78 60 67
60 50 67
93 105 85
85 76 93
85 105 97
85 67 76
85 97 78
78 67 85
124 112 131
124 134 115
124 131 143
143 134 124
94 81 74
128 138 119
146 152 169
169 154 146
169 152 186
187 154 169
186 208 169
169 208 187
132 145 150
150 145 167
167 182 150
182 164 150
149 139 163
149 166 141
163 183 149
183 166 149
120 139 130
111 120 130
130 122 111
141 122 130
130 149 141
139 149 130
163 139 148
181 163 148
148 164 181
110 120 101
129 139 120
120 110 129
129 110 121
129 148 139
555 546 529
546 514 529
529 514 500
529 542 555
457 476 458
458 437 435
435 457 458
500 514 493
493 476 500
493 514 508
473 498 500
500 476 473
473 476 457
453 468 443
453 425 437
443 425 453
483 467 454
454 468 483
492 467 483
483 468 495
483 509 492
495 509 483
547 554 575
516 527 547
527 554 547
575 554 582
605 575 582
582 611 605
600 611 582
582 571 600
582 554 571
439 464 449
449 462 441
449 464 478
478 462 449
449 441 428
428 439 449
479 462 489
489 462 478
142 172 153
143 131 153
153 131 142
222 200 239
190 200 222
302 282 319
251 288 272
272 231 251
251 226 267
267 288 251
402 385 411
411 426 402
443 426 411
411 425 443
157 144 171
180 158 157
157 171 180
136 144 157
306 286 268
306 268 289
289 325 306
325 341 306
351 382 368
382 396 368
415 384 404
404 391 423
404 432 415
423 432 404
371 339 357
357 384 371
359 391 374
374 341 359
391 404 374
374 404 384
357 341 374
374 384 357
361 395 377
377 344 361
377 409 391
395 409 377
359 344 377
377 391 359
469 456 482
485 461 469
469 496 485
482 496 469
461 432 446
446 469 461
456 469 446
433 456 446
423 433 446
446 432 423
349 314 334
334 314 298
297 277 317
561 549 535
561 535 552
587 576 561
576 549 561
561 580 587
552 580 561
383 369 400
352 369 383
531 549 557
40 49 54
40 54 47
28 35 40
35 49 40
47 33 40
40 25 28
33 25 40
697 701 702
697 693 701
702 695 697
697 695 690
690 693 697
629 640 609
591 575 605
591 609 579
591 605 629
629 609 591
677 667 655
655 667 638
559 589 574
574 543 559
559 543 530
363 328 347
329 364 347
311 329 347
347 328 311
364 380 347
380 363 347
291 328 310
310 328 346
269 291 310
310 290 269
344 308 326
308 290 326
290 310 326
326 310 346
326 346 361
361 344 326
319 282 299
219 234 260
255 219 260
260 295 255
260 304 295
83 95 75
83 72 92
64 72 83
83 75 64
92 112 104
104 83 92
95 83 104
115 95 104
104 124 115
112 124 104
75 89 68
62 52 68
68 56 75
68 52 56
170 154 187
187 207 170
170 207 185
170 185 155
229 249 205
205 185 229
185 207 229
229 273 249
86 73 90
90 100 86
80 89 99
99 90 80
80 90 73
80 73 62
62 68 80
80 68 89
59 46 62
62 73 59
58 71 66
119 126 107
107 100 119
155 138 147
138 128 147
147 170 155
154 170 147
137 128 118
137 127 146
137 118 127
146 154 137
154 147 137
137 147 128
91 110 101
101 82 91
91 71 81
82 71 91
148 129 140
132 150 140
140 121 132
140 129 121
140 150 164
164 148 140
500 498 520
520 529 500
520 498 522
522 542 520
542 529 520
471 458 476
476 493 471
437 458 471
471 453 437
516 547 532
517 489 506
506 490 516
478 490 506
506 489 478
516 532 506
506 532 517
505 488 479
479 489 505
515 488 505
505 489 517
505 533 515
517 533 505
178 172 156
178 201 172
178 156 159
267 226 244
226 201 244
168 152 143
143 153 168
186 152 168
286 306 323
323 339 303
303 286 323
357 339 323
323 341 357
323 306 341
339 371 355
398 384 415
398 371 384
415 432 447
447 432 461
472 461 485
472 447 461
460 447 472
352 334 318
318 334 298
318 281 297
298 281 318
277 263 301
301 317 277
301 263 283
349 334 366
366 334 352
352 383 366
549 531 518
535 549 518
501 499 474
494 518 501
501 518 531
557 573 544
544 531 557
544 573 566
566 540 544
540 523 544
627 609 640
627 655 638
640 655 627
677 655 666
674 687 666
687 677 666
666 655 640
666 653 674
666 640 653
548 559 530
548 530 515
515 533 548
262 222 239
262 299 282
351 368 335
319 299 335
316 296 333
316 335 299
316 333 351
351 335 316
234 258 280
280 260 234
280 258 300
304 260 280
280 300 321
321 304 280
187 208 230
230 207 187
253 274 230
230 208 253
99 118 109
109 90 99
118 128 109
100 90 109
119 100 109
109 128 119
74 81 63
81 71 63
53 58 43
71 58 53
53 63 71
70 73 86
70 59 73
43 58 48
55 39 48
66 55 48
48 58 66
39 34 48
48 34 43
26 27 18
26 39 27
26 34 39
94 107 113
113 107 126
113 126 132
132 121 113
87 107 94
87 94 74
86 100 87
100 107 87
87 70 86
110 91 102
94 113 102
81 94 102
102 91 81
121 110 102
102 113 121
484 493 508
484 471 493
508 495 484
495 468 484
468 453 484
453 471 484
562 548 533
579 562 550
550 562 533
550 533 517
517 532 550
560 547 575
560 532 547
560 550 532
579 550 560
575 591 560
560 591 579
190 222 193
193 178 159
159 161 193
193 161 190
267 244 285
282 302 264
302 285 264
264 285 244
226 251 209
209 251 231
209 231 186
186 168 209
188 153 172
188 168 153
188 209 168
226 209 188
188 201 226
172 201 188
358 340 373
373 385 402
402 388 373
373 388 358
324 340 358
324 309 288
324 358 343
343 309 324
319 335 353
353 335 368
385 370 399
399 411 385
414 425 399
425 411 399
322 285 302
497 472 485
485 511 497
497 511 525
474 459 475
475 501 474
494 501 475
470 494 475
372 401 387
387 354 372
369 354 387
400 369 387
387 424 400
387 401 424
389 401 372
372 355 389
389 355 371
371 398 389
434 398 415
415 447 434
434 447 460
424 401 436
352 318 336
317 354 336
297 317 336
336 318 297
336 369 352
336 354 369
337 355 372
372 354 337
317 301 337
337 354 317
521 544 523
531 544 521
523 499 521
499 501 521
521 501 531
607 627 638
638 624 607
624 589 607
579 609 594
609 627 594
627 607 594
594 562 579
296 316 278
278 262 239
299 262 278
278 316 299
239 256 278
256 296 278
207 230 252
273 229 252
252 229 207
294 273 252
252 274 294
252 230 274
178 193 220
201 178 220
220 244 201
220 264 244
370 353 386
386 368 396
386 353 368
386 410 414
386 396 410
414 399 386
386 399 370
340 322 356
356 370 385
385 373 356
356 373 340
338 353 370
370 356 338
338 356 322
338 322 302
338 302 319
319 353 338
267 285 305
285 322 305
305 322 340
305 288 267
305 324 288
340 324 305
535 518 510
510 518 494
525 535 510
510 497 525
472 497 486
470 460 486
460 472 486
497 510 486
486 494 470
486 510 494
459 436 455
470 475 455
455 475 459
401 389 417
417 436 401
339 355 320
355 337 320
303 339 320
320 337 301
320 283 303
320 301 283
578 607 589
578 594 607
562 594 578
589 559 578
559 548 578
548 562 578
237 193 222
237 220 193
237 262 282
222 262 237
282 264 237
264 220 237
413 389 398
413 417 389
398 434 413
438 455 436
436 417 438
417 413 438
438 413 434
438 434 460
438 460 470
470 455 438
])

star() = (STAR_P,STAR_T)


text:
using AssetRegistry
using Base64: stringmime

const blinksetup = joinpath(
    dirname(@__FILE__), "..", "..",
    "packages", "blink-provider", "dist", "blink.js"
)

using Sockets

struct BlinkConnection <: WebIO.AbstractConnection
    page::Blink.Page
end

Blink.body!(p::Blink.Page, x::AbstractWidget) = Blink.body!(p, Widgets.render(x))
function Blink.body!(p::Blink.Page, x::Union{Node, Scope})
    wait(p)

    bs = AssetRegistry.register(blinksetup)
    Blink.loadjs!(p, bs)

    conn = BlinkConnection(p)
    Blink.handle(p, "webio") do msg
        WebIO.dispatch(conn, msg)
    end

    Blink.body!(p, stringmime(MIME"text/html"(), x))
end

Blink.body!(p::Blink.Window, x::AbstractWidget) = Blink.body!(p, Widgets.render(x))
function Blink.body!(p::Blink.Window, x::Union{Node, Scope})
    Blink.body!(p.content, x)
end

function Sockets.send(b::BlinkConnection, data)
    Blink.msg(b.page, Dict(:type=>"webio", :data=>data))
end

Base.isopen(b::BlinkConnection) = Blink.active(b.page)

function WebIO.register_renderable(T::Type, ::Val{:blink})
    eval(:(Blink.body!(p::Union{Blink.Window, Blink.Page}, x::$T) =
           Blink.body!(p, WebIO.render(x))))
end

WebIO.setup_provider(::Val{:blink}) = nothing  # blink setup has no side-effects
WebIO.setup(:blink)


text:
module samplers

using Distributions, LinearAlgebra
using StatsBase
using Printf
using CSV
using DataFrames

include("outFiles.jl")
include("misc.jl")

export runSampler

#main sampler
function runSampler(rowID,Y,X,Z,chainLength,burnIn,outputFreq,priorVCV,varM_prior,M,paths2maps,rS) ##varE will be fixed for now

	println("varM_prior $(varM_prior)")
	#output settings
	these2Keep  = collect((burnIn+outputFreq):outputFreq:chainLength) #print these iterations

        #some info
	##This is not really nFix, but the "blocks" of fixed effects
        nFix  = length(X)
	nRand = length(Z)
	nData = length(Y)
	nMarkerSets = length(M)
	nMarkers    = [size(M[m],2) for m in keys(M)]
	println("nMarkers: $nMarkers")

        #initial computations and settings
	ycorr = deepcopy(Y)

	##make iXpX, Z', zpz (for uncor)
        iXpX = deepcopy(X)
        for x in keys(X)
                iXpX[x] = inv(X[x]'X[x])
        end

	Zp  = deepcopy(Z) #similar(Z')
#	zpz = Array{Array{Float64, 1},1}(undef,0)
	zpz = Dict{Any,Any}()
	for z in keys(Z)
                zpz[z] = diag(Z[z]'Z[z])
		Zp[z]  = Z[z]'
        end

        #key positions for speed
        XKeyPos = Dict{Any,Any}()
        for xSet in keys(X)
              pos = findall(xSet.==collect(keys(X)))[]
              XKeyPos[xSet] = pos
        end
        println("XKeyPos: $XKeyPos")

        ##make b and u arrays
        b = Array{Array{Float64, 1},1}(undef,0)
        ##counts columns per effect
        nColEachX = []
        for x in keys(X)
                println(x)
                nCol = size(X[x],2)
                push!(b,fill(0.0,nCol))
                nColEachX = push!(nColEachX,nCol)
        end

        u = Array{Array{Float64, 1},1}(undef,0)
        ##counts columns per effect
        nColEachZ = []
	##get priors per effect
	iVarStr = [] #inverses will be computed
	varU_prior = []
        for z in keys(Z)
                println(z)
                nCol = size(Z[z],2)
                push!(u,fill(0.0,nCol))
                nColEachZ = push!(nColEachZ,nCol)
		#var structures and priors
		if isempty(priorVCV[z][1])
			println("priorVCV $z is empty, an identity matrix will be used")
			push!(iVarStr,Matrix(1.0I,nCol,nCol))
		else 	push!(iVarStr,inv(priorVCV[z][1]))
		end
		push!(varU_prior,priorVCV[z][2])
        end
	println("prior variances $(varU_prior)")

	#set up for E
	isempty(last(priorVCV,1)[1][1]) ? strE = Matrix(1.0I,nData,nData) : strE = last(priorVCV,1)[1][1]
	varE_prior = last(priorVCV,1)[1][2] #since last returns a tupple

	#parameters for priors
        dfE = 4.0
	dfDefault = 4.0

	dfM = Array{Float64}(undef,0)
        for m in 1:nMarkerSets
		size(varM_prior[m],1)>1 ? println("multivariate prior for marker set $m df=$(3+size(varM_prior[m],1))") : println("univariate prior for marker set $m df=$(3+size(varM_prior[m],1))")
                push!(dfM,3.0+size(varM_prior[m],1))
        end

	println("dfM $dfM")

	if varE_prior==0.0
		varE_prior  = 0.0005
       		scaleE     = 0.0005
        else
       		scaleE    = varE_prior*(dfE-2.0)/dfE
   	end

	##no 0.0005 prior adapted here yet, also no correlated random effects
	scaleU = zeros(nRand)
	for z in 1:nRand
		scaleU[z] = varU_prior[z]*(dfDefault-2.0)/dfDefault
	end

	scaleM = Array{Any,1}(undef,nMarkerSets)
	for m in 1:nMarkerSets
		nMComp = size(varM_prior[m],1)
                nMComp > 1 ? scaleM[m] = varM_prior[m].*(dfM[m]-nMComp-1.0)  : scaleM[m] = varM_prior[m]*(dfM[m]-2.0)/dfM[m] #I make float and array of float
        end

	#pre-computations using priors, not relevant for correlated random effects
   	νS_E = scaleE*dfE
	νS_U = zeros(nRand)
	for z in 1:nRand
                νS_U[z] = scaleU[z]*dfDefault
        end

	#ADD MARKERS
		# read map file and make regions
	regionArray =  Array{Array{UnitRange{Int64},1},1}(undef,0)
	for m in 1:nMarkerSets
		theseRegions = prep2RegionData(paths2maps[m],rS[m]) ###first data
		push!(regionArray,theseRegions)
	end
	println("size regionArray: $(length(regionArray))")
	println("size regionArray: $(length.(regionArray))")

	nRegions  = length.(regionArray) #per component

		#make mpm
#		Mp  = deepcopy(M) #not needed coz I use BLAS.dot
       		mpm = Dict{Any,Any}()
       		for m in keys(M)
               		mpm[m] = diag(M[m]'M[m]) #will not work for large matrices!!!!
#                	Mp[m]  = M[m]'
        	end
		#key positions for speed
		MKeyPos = Dict{Any,Any}()
		for mSet in keys(M)
			pos = findall(mSet.==collect(keys(M)))[]
			MKeyPos[mSet] = pos
		end
		println("MKeyPos: $MKeyPos")
		#storage

	varU = varU_prior #for storage

	###FIXED FOR NOW
	mpmMat = 0

	beta = zeros(Float64,nMarkerSets,maximum(nMarkers)) #can allow unequal length! Remove tail zeros for printing....
#	vcovBeta = fill(Matrix(Diagonal(varM)),maximum(nRegions)) #can allow unequal length! Remove tail zeros for printing....

	varBeta = Array{Any,1}(undef,0)
        for m in 1:nMarkerSets
                push!(varBeta,fill(varM_prior[m],nRegions[m]))
        end

	varBeta2 = Dict{Any,Any}()
	for mSet in keys(M)
		varBeta2[mSet] = fill(varM_prior[MKeyPos[mSet]],nRegions[MKeyPos[mSet]]) #later, direct reference to key when varM_prior is a dictionary
	end
	println("keys of varBeta2: $(keys(varBeta2))")

	#Start McMC
        for iter in 1:chainLength
		#sample residual variance
	       	varE = sampleVarE(νS_E,ycorr,dfE,nData)

		#sample fixed effects
        	#always returns corrected Y and new b
	        sampleX!(X,b,iXpX,nFix,nColEachX,XKeyPos,ycorr,varE)

		#sample random effects
		# always returns corrected Y and new u
		sampleZ!(iVarStr,Z,Zp,zpz,nRand,varE,varU,u,ycorr)

		#sample variances
		sampleRanVar!(varU,nRand,νS_U,u,dfDefault,iVarStr)

		#sample marker effects
@time	        sampleM!(M,beta,mpm,nMarkerSets,MKeyPos,regionArray,nRegions,ycorr,varE,varBeta)

		#sample marker variances
		sampleMarkerVar!(beta,varBeta,nMarkerSets,regionArray,scaleM,dfM)
           sampleMarkerVar2!(beta,varBeta2,nMarkerSets,MKeyPos,nRegions,regionArray,scaleM,dfM)

        	#print
		if iter in these2Keep
			IO.outMCMC(pwd(),"b",vcat(b...)') ### currently no path is provided!!!!
			IO.outMCMC(pwd(),"u",vcat(u...)')
			IO.outMCMC(pwd(),"varE",varE)
			IO.outMCMC(pwd(),"varU",varU')
			for markers in 1:nMarkerSets
				IO.outMCMC(pwd(),"varM$(markers)",varBeta[markers]')
			end
	#		if onScreen==true
            			println("b, $(vcat(b...))") #i always vectorize b. maybe better to make it vector initially
				println("vG, $(var(M["M1"]*beta[1,:]))")
        #		end
		end
	end
end

#Sampling fixed effects

function sampleX!(X,b,iXpX,nFix,nColEachX,keyX,ycorr,varE)
        #block for each effect
        for x in keys(X)
                ycorr    .+= X[x]*b[keyX[x]]
                rhs      = X[x]'*ycorr
                meanMu   = iXpX[x]*rhs
                if nColEachX[keyX[x]] == 1
                        b[keyX[x]] .= rand(Normal(meanMu[],sqrt((iXpX[x]*varE))[]))
                else b[keyX[x]] .= rand(MvNormal(vec(meanMu),convert(Array,Symmetric(iXpX[x]*varE))))
                end
                ycorr    .-= X[x]*b[keyX[x]]
        end
end

#Sampling random effects
function sampleZ!(iStrMat,Zmat,ZpMat,zpzMat,nRand,varE,varU,u,ycorr)
	#block for each effect
	for z in keys(Zmat)
		pos = findall(z.==collect(keys(Zmat)))[]
		uVec = deepcopy(u[pos])
		iMat = iStrMat[z]
		tempzpz = zpzMat[z] ###added
		λz = varE/(varU[pos])
	        ycorr .+= Zmat[z]*uVec
	        Yi = ZpMat[z]*ycorr #computation of Z'ycorr for ALL  rhsU
		nCol = length(zpzMat[z])
	        for i in 1:nCol
        	        uVec[i] = 0.0 #also excludes individual from iMat! Nice trick.
			rhsU = Yi[i] - λz*dot(iMat[:,i],uVec)
                	lhsU = tempzpz[i] + (iMat[i,i]*λz)[1]
			invLhsU = 1.0/lhsU
                	meanU = invLhsU*rhsU
                	uVec[i] = rand(Normal(meanU,sqrt(invLhsU*varE)))
        	end
		u[pos] = uVec
        	ycorr .-= Zmat[z]*uVec
	end
end

#Sampling marker effects
function sampleM!(MMat,beta,mpmMat,nMSet,keyM,regionsMat,regions,ycorr,varE,varBeta)
        #for each marker set
        for mSet in keys(MMat)
                for r in 1:regions[keyM[mSet]]
                        theseLoci = regionsMat[keyM[mSet]][r]
                        regionSize = length(theseLoci)
                        lambda = varE/(varBeta[keyM[mSet]][r])
                        for locus in theseLoci
                                BLAS.axpy!(beta[keyM[mSet],locus],MMat[mSet][:,locus],ycorr)
                                rhs = BLAS.dot(MMat[mSet][:,locus],ycorr)
                                lhs = mpmMat[mSet][locus] + lambda
                                meanBeta = lhs\rhs
                                beta[keyM[mSet],locus] = sampleBeta(meanBeta, lhs, varE)
                                BLAS.axpy!(-1.0*beta[keyM[mSet],locus],MMat[mSet][:,locus],ycorr)
                        end
                end
        end
end

#sample marker effects
function sampleBeta(meanBeta, lhs, varE)
    return rand(Normal(meanBeta,sqrt(lhs\varE)))
end

#sample random effects' variances
function sampleRanVar!(varU,nRand,νS_ranVar,effVec,df_ranVar,iStrMat)
	for z in 1:nRand
		n = size(iStrMat[z],2)
		varU[z] = (νS_ranVar[z] + effVec[z]'*iStrMat[z]*effVec[z])/rand(Chisq(df_ranVar + n))
	end
end

#Sampling marker effects' variances
function sampleMarkerVar!(beta,varBeta,nMSet,regionsMat,scaleM,dfM)
        #for each marker set
        for mSet in 1:nMSet
                for r in 1:length(regionsMat[mSet]) #dont have to compute 1000000 times, take it out
                        theseLoci = regionsMat[mSet][r]
                        regionSize = length(theseLoci)
@time			varBeta[mSet][r] = sampleVarBeta(scaleM[mSet],dfM[mSet],beta[mSet,theseLoci],regionSize)
                end
        end
end

#function sampleMarkerVar2!(beta,varBeta2,nMSet,keyM,regions,regionsMat,scaleM,dfM)
#        #for each marker set
#        for mSet in keys(varBeta2)
#                for r in 1:regions[keyM[mSet]] #dont have to compute 1000000 times, take it out
#                        theseLoci = regionsMat[keyM[mSet]][r]
#                        regionSize = length(theseLoci)
#@time                        varBeta2[mSet][r] = sampleVarBeta(scaleM[keyM[mSet]],dfM[keyM[mSet]],beta[keyM[mSet],theseLoci],regionSize)
#                end
#        end
#end

function sampleMarkerVar2!(beta,varBeta2,nMSet,keyM,regions,regionsMat,scaleM,dfM)
        #for each marker set
        for mSet in keys(varBeta2)
                for r in regionsMat[keyM[mSet]] #dont have to compute 1000000 times, take it out
                       	println("r : $r")
			print(r)
			theseLoci = collect(r)
                        regionSize = length(theseLoci)
@time                        varBeta2[mSet][r] = sampleVarBeta(scaleM[keyM[mSet]],dfM[keyM[mSet]],beta[keyM[mSet],theseLoci],regionSize)
                end
        end
end

#sample marker variances
function sampleVarBeta(scalem,dfm,whichLoci,regionSize)
	return (scalem*dfm + dot(whichLoci,whichLoci)) / rand(Chisq(dfm + regionSize))
end

#Sample residual variance
function sampleVarE(νS_e,yCorVec,df_e,nRecords)
    return((νS_e + dot(yCorVec,yCorVec))/rand(Chisq(df_e + nRecords)))
end

end

max_stars_repo_name: animBreedTools/NextGP

max_issues_repo_name: animBreedTools/NextGP

max_forks_repo_name: animBreedTools/NextGP


text:
"""
    EconometricModel(f::FormulaTerm, data::AbstractDataFrame;
                     contrasts::Dict{Symbol} = Dict{Symbol,Union{<:AbstractContrasts,<:AbstractTerm}})

    Use fit(EconometricModel, f, data, contrasts = contrasts)
    Formula has syntax: @formula(response ~ exogenous + (endogenous ~ instruments) +
                                 weights(wts))
    For absorbing categorical features use the term `absorb(features)`
    For the between estimator use the term `between(features)`
    For the one-way random effects model use the terms `PID(pid) + TID(tid)`
"""
mutable struct EconometricModel{E<:ModelEstimator,
                                F<:FormulaTerm,
                                Y<:AbstractVecOrMat{<:Number},
                                W<:FrequencyWeights,
                                Ŷ<:AbstractVecOrMat{<:Float64},
                                N<:Tuple{<:Union{<:AbstractVector{<:AbstractString},
                                                 <:AbstractString},
                                         <:AbstractVector{<:AbstractString}}} <: EconometricsModel
    estimator::E
    f::F
    data::DataFrame
    X::Matrix{Float64}
    y::Y
    w::W
    β::Vector{Float64}
    Ψ::Hermitian{Float64,Matrix{Float64}}
    ŷ::Ŷ
    vars::N
    iv::Int
end
function show(io::IO, obj::EconometricModel{<:LinearModelEstimators})
    show(io, obj.estimator)
    ℓℓ = loglikelihood(obj)
    println(io, @sprintf("Number of observations: %i", nobs(obj)))
    println(io, @sprintf("Loglikelihood: %.2f", ℓℓ))
    println(io, @sprintf("R-squared: %.4f", r2(obj)))
    W, F, p = wald(obj)
    if !isnan(p)
        println(io, @sprintf("Wald: %.2f ∼ F(%i, %i) ⟹  Pr > F = %.4f", W, params(F)..., p))
    end
    println(io, string("Formula: ", obj.f))
    show(io, coeftable(obj))
end
function show(io::IO, obj::EconometricModel{<:ContinuousResponse})
    show(io, obj.estimator)
    ℓℓ₀ = nullloglikelihood(obj)
    ℓℓ = loglikelihood(obj)
    absorbed = !isempty(obj.estimator.groups)
    println(io, @sprintf("Number of observations: %i", nobs(obj)))
    absorbed || println(io, @sprintf("Null Loglikelihood: %.2f", ℓℓ₀))
    println(io, @sprintf("Loglikelihood: %.2f", ℓℓ))
    println(io, @sprintf("R-squared: %.4f", r2(obj)))
    if !absorbed
        lr = 2(ℓℓ - ℓℓ₀)
        vars = coefnames(obj)
        k = count(x -> !occursin("(Intercept)", x), vars[2])
        if k > zero(k)
            χ² = Chisq(k)
            p = ccdf(χ², lr)
            println(io, @sprintf("LR Test: %.2f ∼ χ²(%i) ⟹  Pr > χ² = %.4f", lr, k, p))
        end
    else
        W, F, p = wald(obj)
        if !isnan(p)
            println(io, @sprintf("Wald: %.2f ∼ F(%i, %i) ⟹  Pr > F = %.4f", W, params(F)..., p))
        end
    end
    println(io, string("Formula: ", obj.f))
    show(io, coeftable(obj))
end
function show(io::IO, obj::EconometricModel{<:NominalResponse})
    show(io, obj.estimator)
    ℓℓ₀ = nullloglikelihood(obj)
    ℓℓ = loglikelihood(obj)
    println(io, @sprintf("Number of observations: %i", nobs(obj)))
    println(io, @sprintf("Null Loglikelihood: %.2f", ℓℓ₀))
    println(io, @sprintf("Loglikelihood: %.2f", ℓℓ))
    println(io, @sprintf("R-squared: %.4f", r2(obj)))
    lr = 2(ℓℓ - ℓℓ₀)
    vars = coefnames(obj)
    k = length(vars[1]) * count(x -> !occursin("(Intercept)", x), vars[2])
    if k > zero(k)
        χ² = Chisq(k)
        p = ccdf(χ², lr)
        println(io, @sprintf("LR Test: %.2f ∼ χ²(%i) ⟹  Pr > χ² = %.4f", lr, k, p))
    end
    println(io, string("Formula: ", obj.f))
    show(io, coeftable(obj))
end
function show(io::IO, obj::EconometricModel{<:OrdinalResponse})
    show(io, obj.estimator)
    ℓℓ₀ = nullloglikelihood(obj)
    ℓℓ = loglikelihood(obj)
    println(io, @sprintf("Number of observations: %i", nobs(obj)))
    println(io, @sprintf("Null Loglikelihood: %.2f", ℓℓ₀))
    println(io, @sprintf("Loglikelihood: %.2f", ℓℓ))
    println(io, @sprintf("R-squared: %.4f", r2(obj)))
    lr = 2(ℓℓ - ℓℓ₀)
    vars = coefnames(obj)
    k = length(vars[2])
    if k > zero(k)
        χ² = Chisq(k)
        p = ccdf(χ², lr)
        println(io, @sprintf("LR Test: %.2f ∼ χ²(%i) ⟹  Pr > χ² = %.4f", lr, k, p))
    end
    println(io, string("Formula: ", obj.f))
    show(io, coeftable(obj))
end
function fit(::Type{<:EconometricModel},
             f::FormulaTerm,
             data::AbstractDataFrame;
             contrasts::Dict{Symbol} = Dict{Symbol,Union{<:AbstractContrasts,<:AbstractTerm}}())
    data, f, exogenous, iv, absorbed, pid, tid, wts, effect, y, X, z, Z =
        decompose(deepcopy(f), data, contrasts)
    ispanel = !isempty(pid[1])
    istime = !isempty(tid[1])
    hdf = !isempty(absorbed[1])
    isbetween = !isempty(effect[1])
    instrumental = !isa(iv.lhs, InterceptTerm)
    if isa(y, AbstractCategoricalVector)
        @assert !isbetween "Between Estimator only defined for continous response"
        @assert !hdf "Absorbing covariates only is only defined for continous response"
        @assert !ispanel "Panel is reserved for the random effects estimator"
        @assert !instrumental "Only exogenous variables are supported for categorical responses"
        estimator = isordered(y) ?
            OrdinalResponse(exogenous.lhs.contrasts) :
            NominalResponse(exogenous.lhs.contrasts)
    elseif isbetween
        @assert !ispanel "Panel ID not required for the between estimator"
        @assert !hdf "Absorbing covariates not implemented with the between estimator"
        estimator = BetweenEstimator(effect[1][1], effect[2])
    elseif ispanel && istime
        @assert !hdf "Absorbing covariates not implemented with the between estimator"
        @assert hasintercept(exogenous) "Random Effects Models require an InterceptTerm"
        estimator = RandomEffectEstimator(pid, tid, X, y, z, Z, wts)
    else
        estimator = ContinuousResponse(absorbed[2])
    end
    X, y, β, Ψ, ŷ, wts, piv = solve(estimator, X, y, z, Z, wts)
    vars = (coefnames(exogenous.lhs),
            convert(Vector{String}, vcat(coefnames(exogenous.rhs), coefnames(iv.lhs))[piv]))
    EconometricModel(estimator, f, data, X, y, wts, β, Ψ, ŷ, vars, size(Z, 2))
end


text:
using KitBase, FluxRC, OrdinaryDiffEq, Plots
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

begin
    x0 = -1
    x1 = 1
    nx = 100
    nface = nx + 1
    dx = (x1 - x0) / nx
    deg = 2 # polynomial degree
    nsp = deg + 1
    u0 = -5
    u1 = 5
    nu = 28
    cfl = 0.1
    dt = cfl * dx
    t = 0.0
    a = 1.0
end

pspace = FluxRC.FRPSpace1D(x0, x1, nx, deg)
vspace = VSpace1D(u0, u1, nu)
δ = heaviside.(vspace.u)

begin
    xFace = collect(x0:dx:x1)
    xGauss = FluxRC.legendre_point(deg)
    xsp = FluxRC.global_sp(xFace, xGauss)
    ll = FluxRC.lagrange_point(xGauss, -1.0)
    lr = FluxRC.lagrange_point(xGauss, 1.0)
    lpdm = FluxRC.∂lagrange(xGauss)
    dgl, dgr = FluxRC.∂radau(deg, xGauss)
end

u = zeros(nx, nsp)
f = zeros(nx, nu, nsp)
for i = 1:nx, ppp1 = 1:nsp
    #u[i, ppp1] = exp(-20.0 * pspace.xp[i, ppp1]^2)
    u[i, ppp1] = 1.0 - sin(π * pspace.xp[i, ppp1])
    prim = conserve_prim(u[i, ppp1], a)
    f[i, :, ppp1] .= maxwellian(vspace.u, prim)
end

e2f = zeros(Int, nx, 2)
for i = 1:nx
    if i == 1
        e2f[i, 2] = nface
        e2f[i, 1] = i + 1
    elseif i == nx
        e2f[i, 2] = i
        e2f[i, 1] = 1
    else
        e2f[i, 2] = i
        e2f[i, 1] = i + 1
    end
end
f2e = zeros(Int, nface, 2)
for i = 1:nface
    if i == 1
        f2e[i, 1] = i
        f2e[i, 2] = nx
    elseif i == nface
        f2e[i, 1] = 1
        f2e[i, 2] = i - 1
    else
        f2e[i, 1] = i
        f2e[i, 2] = i - 1
    end
end

function mol!(du, u, p, t) # method of lines
    dx, e2f, f2e, a, velo, weights, δ, deg, ll, lr, lpdm, dgl, dgr = p

    ncell = size(u, 1)
    nu = size(u, 2)
    nsp = size(u, 3)

    ρ = similar(u, ncell, nsp)
    M = similar(u, ncell, nu, nsp)
    for i = 1:ncell, k = 1:nsp
        #ρ[i, k] = discrete_moments(u[i, :, k], weights)
        ρ[i, k] = sum(u[i, :, k] .* weights)
        #prim = conserve_prim(ρ[i, k], a)
        prim = [ρ[i, k], a, 1.0]
        #M[i, :, k] .= maxwellian(velo, prim)
        @. M[i, :, k] = prim[1] * sqrt(prim[3] / π) * exp(-prim[3] * (velo - prim[2])^2)
    end
    τ = 2.0 * 0.001

    f = similar(u)
    for i = 1:ncell, j = 1:nu, k = 1:nsp
        J = 0.5 * dx[i]
        f[i, j, k] = velo[j] * u[i, j, k] / J
    end

    u_face = zeros(eltype(u), ncell, nu, 2)
    f_face = zeros(eltype(u), ncell, nu, 2)
    for i = 1:ncell, j = 1:nu, k = 1:nsp
        # right face of element i
        u_face[i, j, 1] += u[i, j, k] * lr[k]
        f_face[i, j, 1] += f[i, j, k] * lr[k]

        # left face of element i
        u_face[i, j, 2] += u[i, j, k] * ll[k]
        f_face[i, j, 2] += f[i, j, k] * ll[k]
    end

    f_interaction = similar(u, nface, nu)
    for i = 1:nface
        @. f_interaction[i, :] =
            f_face[f2e[i, 1], :, 2] * (1.0 - δ) + f_face[f2e[i, 2], :, 1] * (δ)
    end

    rhs1 = zeros(eltype(u), ncell, nu, nsp)
    for i = 1:ncell, j = 1:nu, ppp1 = 1:nsp, k = 1:nsp
        rhs1[i, j, ppp1] += f[i, j, k] * lpdm[ppp1, k]
    end

    for i = 1:ncell, j = 1:nu, ppp1 = 1:nsp
        du[i, j, ppp1] =
            -(
                rhs1[i, j, ppp1] +
                (f_interaction[e2f[i, 2], j] - f_face[i, j, 2]) * dgl[ppp1] +
                (f_interaction[e2f[i, 1], j] - f_face[i, j, 1]) * dgr[ppp1]
            ) + (M[i, j, ppp1] - u[i, j, ppp1]) / τ
    end
end

tspan = (0.0, 0.5)
p = (pspace.dx, e2f, f2e, a, vspace.u, vspace.weights, δ, deg, ll, lr, lpdm, dgl, dgr)
prob = ODEProblem(mol!, f, tspan, p)
sol = solve(
    prob,
    ROCK4(),
    saveat = tspan[2],
    #reltol = 1e-8,
    #abstol = 1e-8,
    adaptive = false,
    dt = 0.0005,
    progress = true,
    progress_steps = 10,
    progress_name = "frode",
    #autodiff = false,
)
prob = remake(prob, u0 = sol.u[end], p = p, t = tspan)

#--- post process ---#
plot(xsp[:, 2], u[:, 2])
begin
    ρ = zeros(nx, nsp)
    for i = 1:nx, j = 1:nsp
        ρ[i, j] = moments_conserve(sol.u[end][i, :, j], vspace.u, vspace.weights)[1]
    end
    plot!(xsp[:, 2], ρ[:, 2])
end


text:
@testset "UTILS" begin
    l_a = String[]
    l_t = Tuple[]
    l_v = Vector{String}()
    l_d = Dict()
    @test RediSearch.isa_list(l_a) == true
    @test RediSearch.isa_list(l_t) == true
    @test RediSearch.isa_list(l_v) == true
    @test RediSearch.isa_list(l_d) == false
end


text:
#
# Use this script to plot the evolution of the boundary layer depth
#

## Plotting mixed layer depth
using Plots, PyPlot

using OceanTurbulenceParameterEstimation
using Oceananigans.Fields: interior

# les analysis
# derivative function
function δ(z, Φ)
    m = length(Φ)-1
    Φz = ones(m)
    for i in 1:m
        Φz[i] = (Φ[i+1]-Φ[i])/(z[i+1]-z[i])
    end
    return Φz
end

"""
Approximates the mixed layer depth for pure convection.
"""
function approximate_mixed_layer_depth(model_time_series, data::OneDimensionalTimeSeries, targets)

        data.constants[:αg]
        Qᵇ = data.boundary_conditions.Qᵇ
        N² = data.boundary_conditions.dbdz_bottom

        Nz = data.grid.Nz
        z = data.grid.zC
        Lz = data.grid.Lz
        Nt = length(data.t)

        # For the LES solution
        h2_les = randn(Nt)
        for i in targets

            b = interior(data.b[i]) # remove halos

            Bz = δ([z...], [b...])
            mBz = maximum(Bz)
            tt = (2*N² + mBz)/3
            bools = Bz .> tt
            zA = (z[1:(Nz-1)] .+ z[2:Nz] ) ./ 2

            h2_les[i] = any(bools) ? -minimum(zA[bools]) : Lz
        end

        # For the model solution
        h2_model = randn(Nt)
        for i in 1:Nt-1
            B = interior(model_time_series.b[i])
            Bz = δ([z...], [B...])
            mBz = maximum(Bz)
            tt = (2*N² + mBz)/3
            bools = Bz .> tt
            zA = (z[1:(Nz-1)] + z[2:Nz] )./2

            h2_model[i] = any(bools) ? -minimum(zA[bools]) : Lz # mixed layer reached the bottom
        end

    return [h2_les, h2_model]
end

LESdata = GeneralStrat

ParametersToOptimize = CATKEParametersRiDependent
RelevantParameters = ParametersToOptimize

params = Parameters(RelevantParameters = RelevantParameters,
               ParametersToOptimize = ParametersToOptimize)

function make_all_the_plots(params)

    directory = pwd() * "/Results/plotGeneralStrat_default_parameters/$(RelevantParameters)/"
    isdir(directory) || mkpath(directory)

    for (i, LEScase) in enumerate(values(LESdata))

        td = OneDimensionalTimeSeries(LEScase.filename; grid_type=OneDimensionalEnsembleGrid, Nz=32)

        relative_weights = Dict(:b => 1.0, :u => 1.0, :v => 1.0, :e => 1.0)

        # Single simulation
        case_loss, default_parameters = get_loss(LEScase, td, params, relative_weights; Δt=10.0, fields=relevant_fields(LEScase),
                                                parameter_specific_kwargs[params.RelevantParameters]...)

        targets = case_loss.loss.targets

        ℱ = model_time_series(case_loss.simulation, case_loss.data, default_parameters)
        h2_les, h2_model = approximate_mixed_layer_depth(ℱ, td, targets)

        days = td.t[targets] ./ 86400
        first = LEScase.first
        last = isnothing(LEScase.last) ? length(td.t) : LEScase.last

        Plots.plot(days, h2_les, label = "LES mixed layer depth", linewidth = 3 , color = :red, grid = true, gridstyle = :dash, gridalpha = 0.25, framestyle = :box , title=td.name, ylabel = "Depth [meters]", xlabel = "days", legend = :topleft)
        p = plot!(days, h2_model, color = :purple, label = "TKE mixed layer depth", linewidth = 3 )
        Plots.savefig(p, directory*td.name*"_mixed_layer_depth.pdf")

        Plots.plot(days, h2_les, label = "LES mixed layer depth", linewidth = 3 , color = :red, grid = true, gridstyle = :dash, gridalpha = 0.25, framestyle = :box , title=td.name, ylabel = "Depth [meters]", xlabel = "days", yscale = :log10, xscale = :log10, legend = :topleft)
        p = plot!(days, h2_model, color = :purple, label = "TKE mixed layer depth", linewidth = 3 )
        Plots.savefig(p, directory*td.name*"_mixed_layer_depth_log10.pdf")

        # p = visualize!(md, td, 1:180:length(td), best_parameters)
        # PyPlot.savefig(directory*td.name*".png")
    end
end

make_all_the_plots(params)


text:
#=
Derived from https://github.com/aparamon/HDF5Plugin-Zstandard, zstd_h5plugin.c
Licensed under Apache License Version 2.0, see licenses/H5Zzstd_LICENSE.txt

The following license applies to the Julia port.
Copyright (c) 2021 Mark Kittisopikul and Howard Hughes Medical Institute. License MIT, see LICENSE.txt
=#
module H5Zzstd

using CodecZstd
import CodecZstd.LibZstd
using HDF5.API
import HDF5.Filters: Filter, filterid, register_filter, filterid, filtername, filter_func, filter_cfunc

const H5Z_FILTER_ZSTD = API.H5Z_filter_t(32015)
const zstd_name = "Zstandard compression: http://www.zstd.net"

export H5Z_filter_zstd, H5Z_FILTER_ZSTD, ZstdFilter

# cd_values First optional value is the compressor aggression
#           Default is CodecZstd.LibZstd.ZSTD_CLEVEL_DEFAULT
function H5Z_filter_zstd(flags::Cuint, cd_nelmts::Csize_t,
                        cd_values::Ptr{Cuint}, nbytes::Csize_t,
                        buf_size::Ptr{Csize_t}, buf::Ptr{Ptr{Cvoid}})::Csize_t
    inbuf = unsafe_load(buf)
    outbuf = C_NULL
    origSize = nbytes
    ret_value = Csize_t(0)

    try

    if flags & API.H5Z_FLAG_REVERSE != 0
        #decompresssion

        decompSize = LibZstd.ZSTD_getDecompressedSize(inbuf, origSize)
        outbuf = Libc.malloc(decompSize)
        if outbuf == C_NULL
            error("zstd_h5plugin: Cannot allocate memory for outbuf during decompression.")
        end
        decompSize = LibZstd.ZSTD_decompress(outbuf, decompSize, inbuf, origSize)
        Libc.free(inbuf)
        unsafe_store!(buf, outbuf)
        outbuf = C_NULL
        ret_value = Csize_t(decompSize)
    else
        # compression

        if cd_nelmts > 0
            aggression = Cint(unsafe_load(cd_values))
        else
            aggression = CodecZstd.LibZstd.ZSTD_CLEVEL_DEFAULT
        end

        if aggression < 1
            aggression = 1 # ZSTD_minCLevel()
        elseif aggression > LibZstd.ZSTD_maxCLevel()
            aggression = LibZstd.ZSTD_maxCLevel()
        end

        compSize = LibZstd.ZSTD_compressBound(origSize)
        outbuf = Libc.malloc(compSize)
        if outbuf == C_NULL
            error("zstd_h5plugin: Cannot allocate memory for outbuf during compression.")
        end

        compSize = LibZstd.ZSTD_compress(outbuf, compSize, inbuf, origSize, aggression)

        Libc.free(unsafe_load(buf))
        unsafe_store!(buf, outbuf)
        unsafe_store!(buf_size, compSize)
        outbuf = C_NULL
        ret_value = compSize
    end
    catch e
        #  "In the case of failure, the return value is 0 (zero) and all pointer arguments are left unchanged."
        ret_value = Csize_t(0)
        @error "H5Zzstd Non-Fatal ERROR: " err
        display(stacktrace(catch_backtrace()))
    finally

        if outbuf != C_NULL
            free(outbuf)
        end

    end # try catch finally
    return Csize_t(ret_value)
end

# Filters Module

"""
    ZstdFilter(clevel)

Zstandard compression filter. `clevel` determines the compression level.

# External Links
* [Zstandard HDF5 Filter ID 32015](https://portal.hdfgroup.org/display/support/Filters#Filters-32015)
* [Zstandard HDF5 Plugin Repository (C code)](https://github.com/aparamon/HDF5Plugin-Zstandard)
"""
struct ZstdFilter <: Filter
    clevel::Cuint
end
ZstdFilter() = ZstdFilter(CodecZstd.LibZstd.ZSTD_CLEVEL_DEFAULT)

filterid(::Type{ZstdFilter}) = H5Z_FILTER_ZSTD
filtername(::Type{ZstdFilter}) = zstd_name
filter_func(::Type{ZstdFilter}) = H5Z_filter_zstd
filter_cfunc(::Type{ZstdFilter}) = @cfunction(H5Z_filter_zstd, Csize_t,
                                             (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                                             Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))

function __init__()
    register_filter(ZstdFilter)
end

end # module H5Zzstd


text:
abstract type SpecialArray{T,n} <: AbstractArray{T,n} end

Base.size(a::SpecialArray) = size(a.values)
Base.IndexStyle(::Type{<:SpecialArray}) = IndexLinear()
#Base.linearindexing{T<:SpecialArray}(::Type{T}) = Base.LinearFast()
Base.getindex(a::SpecialArray, i::Int) = a.values[i]
Base.convert(::Type{Array{T,n}}, x::SpecialArray) where {T,n} = x.values
Base.promote_rule(::Type{AbstractArray{T,n}},::SpecialArray{T,n}) where {T,n} =
  Type{Array{T,n}}

"""
    ParameterValuesArray{T,n}

Same as Array, but tracks if user has modified values.
"""
mutable struct ParameterValuesArray{T,n} <: SpecialArray{T,n}
  values :: Array{T,n}
  ismodified :: Bool
end
function Base.setindex!(a::ParameterValuesArray, v, i::Int)
  if ~isfinite(v)
    throw(DomainError("$v is not valid LTspice parameter value"))
  end
  a.ismodified = true
  a.values[i] = v
end
Base.convert(::Type{ParameterValuesArray{T,n}}, x::AbstractArray{S,n}) where {T,n,S} =
  ParameterValuesArray{T,n}(Array{T,n}(x),true)
Base.convert(::Type{LTspice.ParameterValuesArray{T,n}}, x::LTspice.ParameterValuesArray{T,n}) where {T,n} =
  ParameterValuesArray{T,n}(Array{T,n}(x),true)
"""
    MeasurementValuesArray{T,n}

Same as Array, but setindex! returns error.
"""
mutable struct MeasurementValuesArray{T,n} <: SpecialArray{T,n}
  values :: Array{T,n}
end
function Base.setindex!(a::MeasurementValuesArray, ::Any, ::Any)
  throw(ErrorException("Measurments may not be set.  Run simulation to update."))
end


text:
# * Diagonal matrices
DiagonalBlockDiagonal(A::AbstractMatrix, rows, cols) =
    BandedBlockBandedMatrix(A, rows,cols, (0,0), (0,0))

DiagonalBlockDiagonal(A::AbstractMatrix, rows) =
    DiagonalBlockDiagonal(A, rows, rows)

function (B::FEDVROrRestricted)(D::Diagonal)
    n = size(B,2)
    @assert size(D) == (n,n)
    all(order(B) .== 2) ? D : DiagonalBlockDiagonal(D, block_structure(B))
end

# * Dense operators
# ** Matrix construction

function Matrix(::UndefInitializer, B::FEDVROrRestricted{T}, ::Type{U}=T) where {T,U}
    if all(order(B) .== 2)
        n = size(B,2)
        dl = Vector{U}(undef, n-1)
        d = Vector{U}(undef, n)
        du = Vector{U}(undef, n-1)
        Tridiagonal(dl, d, du)
    else
        rows = block_structure(B)
        l,u = block_bandwidths(B,rows)

        BlockSkylineMatrix{U}(undef, rows, rows, (l,u))
    end
end

function Base.zeros(B::FEDVROrRestricted{T}, ::Type{U}=T) where {T,U}
    M = Matrix(undef, B, U)
    M .= zero(U)
    M
end

# ** Individual elements access

function set_element_common!(fun::Function, A, B::FEDVROrRestricted, i)
        els = elements(B)
    i ∉ els && throw(ArgumentError("Invalid element $(i) ∉ $(els)"))
    el1,eln = first(els),last(els)

    is = indices(B,2)

    elb = element_boundaries(parent(B))
    # If the first element only contributes through the bridge
    # function, also consider the second element as the "first"
    # element, since it will be exactly in the north-west corner of
    # the matrix.
    el1′ = el1 + (first(is) == elb[el1+1] ? 1 : 0)
    isfirst = i == el1 || i == el1′
    # If the last element only contributes through the bridge
    # function, also consider the next-to-last element as the "last"
    # element, since it will be exactly in the south-east corner of
    # the matrix.
    eln′ = eln - (last(is) == elb[eln] ? 1 : 0)
    islast = i == eln || i == eln′

    b = fun(i)
    o = order(B,i)
    @assert size(b,1) == size(b,2) == o

    # Find start and end of the interior block; the first and last
    # elements (within the restriction) have their interior blocks as
    # their outer blocks, which are cropped as the restriction
    # prescribes.
    s = 1 + (i == el1 ? first(is)-elb[el1] : (i == el1′ ? first(is)-elb[el1′] : 1))
    e = o - (i == eln ? elb[eln+1]-last(is) : (i == eln′ ? elb[eln]-last(is) : 1))

    # First element only contributes through the bridge function
    if isfirst && s == o
        A[1,1] += b[end,end]
        return
    end

    # Last element only contributes through the bridge function
    if islast && e == 1
        A[end,end] += b[1,1]
        return
    end

    b,o,s:e,isfirst,islast,el1′
end

function set_element!(fun::Function, A::BlockSkylineMatrix, B::FEDVROrRestricted, i)
    res = set_element_common!(fun, A, B, i)
    isnothing(res) && return
    b,o,sel,isfirst,islast,el1′ = res

    # Find block index of the north-west block of element i; this is a
    # linear search which is not optimal for large number of elements.
    j = isfirst ? 1 : 0
    if !isfirst
        for el in el1′:i-1
            j += 1 + (order(B,el) > 2 || el == el1′ ? 1 : 0)
        end
    end

    if o > 2
        ji = j+!isfirst
        # Interior block
        A[Block(ji,ji)] += b[sel,sel]

        if !isfirst
            # North-west corner
            A[Block(j,j)] .+= b[1,1]
            # West edge
            A[Block(j+1,j)] = b[sel,1:1]
            # North edge
            A[Block(j,ji)] = b[1:1,sel]
            if !islast
                # Nort-east corner
                A[Block(j,ji+1)] = b[1:1,end:end]
                # South-west corner
                A[Block(j+2,j)] = b[end:end,1:1]
            end
        end
        if !islast
            # South-east corner
            A[Block(ji+1,ji+1)] .+= b[end,end]
            # East edge
            A[Block(ji,ji+1)] = b[sel,end:end]
            # South edge
            A[Block(ji+1,ji)] = b[end:end, sel]
        end
    else
        A[Block(j,j)] .+= b[1,1]
        A[Block(j+1,j+1)] .+= b[2,2]
        A[Block(j,j+1)] .= b[1,2]
        A[Block(j+1,j)] .= b[1,2]
    end
end

function set_element!(fun::Function, A::Tridiagonal, B::FEDVROrRestricted, i)
    res = set_element_common!(fun, A, B, i)
    isnothing(res) && return
    b,o,sel,isfirst,islast,el1′ = res

    ii = i - el1′ + 1

    A.du[ii] = b[1,2]
    A.d[ii:ii+1] .+= diag(b)
    A.dl[ii] = b[2,1]
end

function set_elements!(fun::Function, A, B::FEDVROrRestricted)
    A .= false
    for i in elements(B)
        set_element!(fun, A, B, i)
    end
    A
end

# * Scalar operators

@simplify function *(Ac::QuasiAdjoint{<:Any,<:FEDVR},
                     D::QuasiDiagonal,
                     B::FEDVR)
    A = parent(Ac)
    A == B || throw(ArgumentError("Cannot multiply functions on different grids"))

    Diagonal(getindex.(Ref(D.diag), B.x))
end

# A & B restricted
@simplify function *(Ac::AdjointRestrictedFEDVR,
                     D::QuasiDiagonal,
                     B::RestrictedFEDVR)
    # This is mainly for type-stability; it would be trivial to
    # generate the proper banded matrix with one off-diagonal, from
    # the combination of two differently restricted bases, but we
    # would like to have Diagonal as the result if they are equal, and
    # this has higher priority. On the other hand, you typically only
    # compute scalar operators in the beginning of the calculation,
    # and thus type-instability is not a big problem, so this
    # behaviour may change in the future.
    reverse(axes(Ac)) == axes(B) || throw(DimensionMismatch("axes must be same"))
    A = parent(Ac)
    parent(A) == parent(B) || throw(ArgumentError("Cannot multiply incompatible FEDVR expansions"))

    a,b = restriction_extents(B)
    Diagonal(getindex.(Ref(D.diag), parent(B).x[1+a:end-b]))
end


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Unk-token pre-processing.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-06-07

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

using ArgParse

using CoNLLX

argparse = ArgParseSettings()
@add_arg_table argparse begin
    "--input", "-i"
        help = "input source (default: STDIN)"
        default = "/dev/stdin"
    "--output", "-o"
        help = "output sink (default: STDOUT)"
        default = "/dev/stdout"
    "--cut-off", "-c"
        help = "frequency cut-off (default: 1)"
        arg_type = Uint
        range_tester = (x) -> x > 0
        default = uint(1)
    "--unk-tok", "-u"
        help = "token representation for UNKs (default: <UNK>)"
        arg_type = String
        default = "<UNK>"
    "vocab_tsv"
        help = "token/occurences pair TSV"
        required = true
end

pargs = parse_args(argparse)

vocab = Set{String}()
open(pargs["vocab_tsv"], "r") do vocab_f
    for line in eachline(vocab_f)
        tok, occs_s = split(line, '\t', 2)
        occs = uint(occs_s)
        if occs > pargs["cut-off"]
            push!(vocab, tok)
        end
    end
end

open(pargs["input"], "r") do input_f
    open(pargs["output"], "w") do output_f
        for sent in conllxparse(input_f)
            for tok in sent
                if !(tok.form in vocab)
                    tok.form = pargs["unk-tok"]
                end
            end
            print(output_f, sent)
            print(output_f, "\n\n")
        end
    end
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
module FindSteadyStates

import DifferentialEquations
import ModelingToolkit

using DifferentialEquations: SteadyStateProblem, solve, Tsit5,AutoTsit5, Rosenbrock23, EnsembleProblem, EnsembleThreads, remake, ODEProblem, solve, SSRootfind, EnsembleSolution

using Random: shuffle!, rand
using LinearAlgebra: eigvals

using Parameters: @with_kw, @unpack

include("types.jl")
include("utils.jl")

# Paramger Generator
include("param_gen/types.jl")
include("param_gen/utils.jl")
include("param_gen/sampling.jl")
include("param_gen/grid_search.jl")

# Solvers
include("solve.jl")

# Stability
include("jacobian.jl")

end


text:
using Luna

a = 9e-6
gas = :Ar
pres = 5
flength = 25e-2

τfwhm = 30e-15
λ0 = 1500e-9
energy = 1.7e-6

grid = Grid.RealGrid(flength, λ0, (200e-9, 3000e-9), 2e-12)

m = Capillary.MarcatiliMode(a, gas, pres, loss=false)
aeff = let m=m
    z -> Modes.Aeff(m, z=z)
end

energyfun, energyfunω = Fields.energyfuncs(grid)

dens0 = PhysData.density(gas, pres)
densityfun(z) = dens0

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

plasma = Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot)
responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),
             plasma)

linop, βfun!, β1, αfun = LinearOps.make_const_linop(grid, m, λ0)

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, βfun!, aeff)

statsfun = Stats.default(grid, Eω, m, linop, transform; gas=gas, windows=((150e-9, 300e-9),))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)

Luna.run(Eω, grid, linop, transform, FT, output)

Plotting.pygui(true)
Plotting.stats(output)
Plotting.prop_2D(output, :λ, λrange=(500e-9,1800e-9), trange=(-500e-15,100e-15), dBmin=-30.0)
Plotting.time_1D(output, [5e-2, 15e-2, 25e-2])
Plotting.spec_1D(output, [5e-2, 15e-2, 25e-2])


text:
include("header.jl")

@testset "cnn" begin

    arrtype = KnetLayers.arrtype

    @testset "conv" begin
        m = Conv(height=3, width=3, inout=3=>5, stride=1, padding=1, mode=1)
        x = arrtype(zeros(10,10,3,2))
        y = m(x)
        @test size(y) == (10,10,5,2)

        m = Conv(height=3, width=3, inout=3=>5, stride=1, padding=1, mode=1, binit=nothing, activation=nothing)
        x = arrtype(zeros(10,10,3,2))
        y = m(x)
        @test size(y) == (10,10,5,2)

        m = DeConv(height=3, width=3, inout=5=>3, stride=2, padding=1, mode=1)
        x = arrtype(zeros(10,10,5,2))
        y = m(x)
        @test size(y) == (19,19,3,2)

    end
end


text:
using MolecularGraph.MolecularGraphModel: lgnodematcher, lgedgematcher
using MolecularGraph:
    fastidentityfilter, fastsubstrfilter,
    atommatch, bondmatch

@testset "substructure" begin

@testset "filter" begin
    BuOH = smilestomol("CCCO")
    butane = smilestomol("CCCC")
    propane = smilestomol("CCC")
    @test fastidentityfilter(butane, BuOH)
    @test !fastidentityfilter(propane, butane)
    @test fastsubstrfilter(butane, propane)
    @test !fastsubstrfilter(propane, butane)
end

@testset "substruct" begin
    null = smilestomol("")
    @test !isstructmatch(null, null, prefilter=false)
    @test !issubstructmatch(null, null, prefilter=false)

    hexane = smilestomol("CCCCCC")
    iso = smilestomol("CC(C)CCC")
    iso2 = smilestomol("CCCC(C)C")
    cyclohexane = smilestomol("C1CCCCC1")
    @test !isstructmatch(hexane, iso, prefilter=false)
    @test isstructmatch(iso, iso2, prefilter=false)
    @test !isstructmatch(cyclohexane, hexane, prefilter=false)
    @test issubstructmatch(cyclohexane, hexane, prefilter=false)

    tms = smilestomol("C[Si](C)(C)C")
    tsm = smilestomol("[Si]C([Si])([Si])[Si]")
    @test !isstructmatch(tms, tsm, prefilter=false)

    sulfide = smilestomol("CSSC")
    disconn = smilestomol("CS.SC")
    @test !isstructmatch(sulfide, disconn, prefilter=false)
    @test issubstructmatch(sulfide, disconn, prefilter=false)

    # Invalid (no edges)
    nacl = smilestomol("[Na+].[Cl-]")
    na = smilestomol("[Na]")
    @test !issubstructmatch(nacl, na, prefilter=false)

    tetrahedrane = smilestomol("C12C3C1C23")
    fused = smilestomol("C1C2C1C2")
    spiro = smilestomol("C1CC12CC2")
    @test issubstructmatch(tetrahedrane, fused, prefilter=false)
    @test !issubstructmatch(tetrahedrane, spiro, prefilter=false)
end

@testset "connectedquery" begin
    anyatom = parse(SMARTS, "*")
    hexane = smilestomol("CCCCCC")
    @test isquerymatch(hexane, anyatom)

    priamine = parse(SMARTS, "[NX3;H2]")
    aniline = smilestomol("C1=CC=CC=C1N")
    dieamine = smilestomol("CCNCC")
    @test isquerymatch(aniline, priamine)
    @test !isquerymatch(dieamine, priamine)

    alcohol = parse(SMARTS, "[#6][OD]")
    diether = smilestomol("CCOCC")
    phenol = smilestomol("c1ccccc1O")
    glycerol = smilestomol("OCC(O)CO")
    @test !isquerymatch(diether, alcohol)
    @test isquerymatch(phenol, alcohol)
    @test isquerymatch(glycerol, alcohol)

    aliphring = parse(SMARTS, "*@;!:*")
    cyclopentane = smilestomol("C1CCCC1")
    pyrrole = smilestomol("N1C=CC=C1")
    @test isquerymatch(cyclopentane, aliphring)
    @test !isquerymatch(pyrrole, aliphring)

    notamide = parse(SMARTS, raw"[NX3;!$(NC=O)]")
    triamine = smilestomol("CCN(CC)CC")
    acetamide = smilestomol("CC(=O)N")
    @test isquerymatch(triamine, notamide)
    @test !isquerymatch(acetamide, notamide)

    peroxide = parse(SMARTS, "[OX2][OX2]")
    po1 = smilestomol("COOC")
    npo1 = smilestomol("COCOCOC")
    @test isquerymatch(po1, peroxide)
    @test !isquerymatch(npo1, peroxide)

    sixmem = parse(SMARTS, "[*r6]1[*r6][*r6][*r6][*r6][*r6]1")
    pyridine = smilestomol("n1ccccc1")
    pyrrole = smilestomol("n1cccc1")
    @test isquerymatch(pyridine, sixmem)
    @test !isquerymatch(pyrrole, sixmem)
end

@testset "disconnectedquery" begin
    disconn = parse(SMARTS, "[#7,#8].[!#6].N")
    hetero1 = smilestomol("CCN(O)[O-]")
    hetero2 = smilestomol("CCN=N")
    hetero3 = smilestomol("BrCCOC#N")
    hetero4 = smilestomol("CCS(=O)(=O)O")
    @test isquerymatch(hetero1, disconn)
    @test !isquerymatch(hetero2, disconn)
    @test isquerymatch(hetero3, disconn)
    @test !isquerymatch(hetero4, disconn)
end

end # substructure


text:
@reexport module ITensorNetworks

using ITensors

using ITensors: data, contract

include("ITensors.jl")
include("networks/lattices.jl")
include("networks/inds_network.jl")
include("networks/itensor_network.jl")
include("MPSTensor/MPSTensor.jl")
include("TreeTensor/TreeTensor.jl")
include("models/models.jl")
include("approximations/approximations.jl")
include("peps/peps.jl")
include("interfaces/sweep_contractor.jl")

end


text:
# FIXME this thing eats up macro calls!
"""
Remove LineNumberNode from quoted blocks of code
"""
rmlines(e::Expr) = Expr(e.head, map(rmlines, filter(x -> !(x isa LineNumberNode), e.args))...)
rmlines(a) = a

# TODO binarize block?
"""
HARD FIX of n-arity of operators in `Expr` trees.
"""
function binarize!(e, ops::Vector{Symbol})
    if !(e isa Expr)
        return e
    end

    start = isexpr(e, :call) ? 2 : 1
    n = length(e.args)

    for i ∈ start:n
        e.args[i] = binarize!(e.args[i], ops)
    end

    if isexpr(e, :call)
        op = e.args[1]
        if op ∈ ops && length(e.args) > 3
            return foldl((x,y) -> Expr(:call, op, x, y), @view e.args[2:end])
        end
    end
    return e
end

# function binarize!(e, ops::Vector{Symbol})
#     df_walk!(binarize_step, e, ops)
# end

function clean_block_step(e)
    if isexpr(e, :block)
        if length(e.args) == 1
            return e.args[1]
        elseif length(e.args) > 3
            return foldl((x,y) -> Expr(:block, x, y), @view e.args[2:end])
        end
    end
    return e
end

function cleanblock(e)
    df_walk!(clean_block_step, e)
end

const binarize_ops = [:(+), :(*)]

"""
Binarize n-ary operators (`+` and `*`) and call [`rmlines`](@ref)
"""
cleanast_rec(ex) = rmlines(ex)  |>
    x -> binarize!(x, binarize_ops)

function cleanast(e::Expr)
    # TODO better line removal
    if isexpr(e, :block)
        return Expr(e.head, filter(x -> !(x isa LineNumberNode), e.args)...)
    end

    # Binarize
    if isexpr(e, :call)
        op = e.args[1]
        if op ∈ binarize_ops && length(e.args) > 3
            return foldl((x,y) -> Expr(:call, op, x, y), @view e.args[2:end])
        end
    end
    return e
end

function cleanast!(e::Expr)
    # TODO better line removal
    if isexpr(e, :block)
        e.args = filter(x -> !(x isa LineNumberNode), e.args)
    end

    # Binarize
    if isexpr(e, :call)
        op = e.args[1]
        if op ∈ binarize_ops && length(e.args) > 3
            ne = foldl((x,y) -> Expr(:call, op, x, y), @view e.args[2:end])
            e.args = ne.args
        end
    end
    return e
end

remove_assertions(e) = df_walk( x -> (isexpr(x, :(::)) ? x.args[1] : x),
    e; skip_call=true )

unquote_sym(e) = df_walk( x -> (x isa QuoteNode && x.value isa Symbol ? x.value : x),
    e; skip_call=true )

function collect_symbols(ex)
    syms = Set{Symbol}()
    df_walk( x -> (if x isa Symbol; push!(syms, x); end; x), ex; skip_call=true )
    syms
end


text:
"""
将三角切分的结果再进行切分
"""
module Refine

using ..Basics

export split_triangle, refine_triangles, find_adjs_by_adjoint

"""
将三角形重新切分成四个，可以用在直角三角形上
依旧保证出来的四个是直角三角形，且第一个角是直角
"""
@generated function split_triangle(rtri::Basics.AbstractTriangle{T}) where T
    type2func = Dict(
        :COM => Triangle,
        :RT => RtTriangle,
        :EQ => EqTriangle
    )
    func = type2func[T]
    if T == :RT
        return quote
            #新的四个直角三角形的四个直角
            rtvexs = Vector{Point2D}(undef, 4)
            rtvexs[1] = rtri.vertex[1]
            rtvexs[2] = middle_point(rtri.vertex[1], rtri.vertex[2])
            rtvexs[3] = middle_point(rtri.vertex[2], rtri.vertex[3])
            rtvexs[4] = middle_point(rtri.vertex[3], rtri.vertex[1])
            #
            newtris = Vector{Basics.AbstractTriangle{T}}(undef, 4)
            newtris[1] = RtTriangle(rtvexs[2], rtvexs[3], rtvexs[1])
            newtris[2] = RtTriangle(rtvexs[2], rtri.vertex[2], rtvexs[3])
            newtris[3] = RtTriangle(rtvexs[4], rtvexs[3], rtri.vertex[3])
            newtris[4] = RtTriangle(rtvexs[4], rtvexs[3], rtvexs[1])
            return newtris
        end
    end
    return quote
        #新的四个直角三角形的四个直角
        rtvexs = Vector{Point2D}(undef, 4)
        rtvexs[1] = rtri.vertex[1]
        rtvexs[2] = middle_point(rtri.vertex[1], rtri.vertex[2])
        rtvexs[3] = middle_point(rtri.vertex[2], rtri.vertex[3])
        rtvexs[4] = middle_point(rtri.vertex[3], rtri.vertex[1])
        #
        newtris = Vector{Basics.AbstractTriangle{T}}(undef, 4)
        newtris[1] = ($func)(rtvexs[1], rtvexs[2], rtvexs[4])
        newtris[2] = ($func)(rtvexs[2], rtri.vertex[2], rtvexs[3])
        newtris[3] = ($func)(rtvexs[3], rtvexs[4], rtvexs[2])
        newtris[4] = ($func)(rtvexs[4], rtvexs[3], rtri.vertex[3])
        return newtris
    end # end quote
end

"""
将一组的三角形细化
"""
function refine_triangles(ltris::Vector{T}) where T <: Basics.AbstractTriangle
    newtris = Matrix{T}(undef, length(ltris), 4)
    for (idx, tri) in enumerate(ltris)
        newtris[idx, :] = split_triangle(tri)
    end
    newtris = reshape(newtris, length(ltris)*4)
    return newtris
end

"""
判断一条边是不是某个三角形的
"""
macro isedge(edge, tri)
    return esc(quote
        topt1 = [vex - ($edge).pt1 for vex in ($tri).vertex]
        topt1 = [sqrt(pt.x^2 + pt.y^2) for pt in topt1]
        haspt1 = isapprox(minimum(topt1), 0., atol=1e-6)
        topt2 = [vex - ($edge).pt2 for vex in ($tri).vertex]
        topt2 = [sqrt(pt.x^2 + pt.y^2) for pt in topt2]
        haspt2 = isapprox(minimum(topt2), 0., atol=1e-6)
        haspt1 && haspt2
    end)
end

"""
通过判断边是否重合来判断两个三角形是否挨在一起
"""
function find_adjs_by_adjoint(ltris::Vector{T}) where T <: Basics.AbstractTriangle
    ladjs = Vector{Tuple{
        Union{Missing, T}, Union{Missing, T}, Union{Missing, T}
    }}(undef, length(ltris))
    for idx in 1:1:length(ltris)
        ladjs[idx] = (missing, missing, missing)
    end
    @Threads.threads for idx in 1:1:length(ltris)
        tri = ltris[idx]
        fadjs = Vector{Union{Missing, T}}(undef, 3)
        fadjs[:] = [missing, missing, missing]
        for (adjidx, adjtri) in enumerate(ltris)
            if idx == adjidx
                continue
            end
            if @isedge tri.edges[1] adjtri
                fadjs[1] = adjtri
                continue
            end
            if @isedge tri.edges[2] adjtri
                fadjs[2] = adjtri
                continue
            end
            if @isedge tri.edges[3] adjtri
                fadjs[3] = adjtri
                continue
            end
        end
        ladjs[idx] = Tuple(fadjs)
    end
    return ladjs
end

end # end module


text:
# Test weighted sum

using NumericExtensions
using Base.Test

xi = rand(1:100, 6)
x = rand(6)
y = rand(6)
z = rand(6)
w = rand(6)

# full reduction

@test_approx_eq wsum(w, xi) sum(w .* xi)
@test_approx_eq wsum(w, x) sum(w .* x)
@test_approx_eq wsum(w, Abs2Fun(), x) sum(w .* abs2(x))
@test_approx_eq wsum(w, Multiply(), x, y) sum(w .* (x .* y))
@test_approx_eq wsum(w, FMA(), x, y, z) sum(w .* (x + y .* z))
@test_approx_eq wsumfdiff(w, Abs2Fun(), x, y) sum(w .* abs2(x - y))

# partial reduction on matrices

xi = rand(1:100, 5, 6)
x = rand(5, 6)
y = rand(5, 6)
w1 = rand(5)
w2 = rand(1, 6)

@test_approx_eq wsum(w1, xi, 1) sum(w1 .* xi, 1)
@test_approx_eq wsum(w2, xi, 2) sum(w2 .* xi, 2)

@test_approx_eq wsum(w1, x, 1) sum(w1 .* x, 1)
@test_approx_eq wsum(w2, x, 2) sum(w2 .* x, 2)

@test_approx_eq wsum(w1, Abs2Fun(), x, 1) sum(w1 .* abs2(x), 1)
@test_approx_eq wsum(w2, Abs2Fun(), x, 2) sum(w2 .* abs2(x), 2)

@test_approx_eq wsum(w1, Multiply(), x, y, 1) sum(w1 .* (x .* y), 1)
@test_approx_eq wsum(w2, Multiply(), x, y, 2) sum(w2 .* (x .* y), 2)

@test_approx_eq wsumfdiff(w1, Abs2Fun(), x, y, 1) sum(w1 .* abs2(x - y), 1)
@test_approx_eq wsumfdiff(w2, Abs2Fun(), x, y, 2) sum(w2 .* abs2(x - y), 2)

# partial reduction on cubes

xi = rand(1:100, 5, 6, 3)
x = rand(5, 6, 3)
y = rand(5, 6, 3)
w1 = rand(5)
w2 = rand(1, 6)
w3 = rand(1, 1, 3)

@test_approx_eq wsum(w1, xi, 1) sum(w1 .* xi, 1)
@test_approx_eq wsum(w2, xi, 2) sum(w2 .* xi, 2)
@test_approx_eq wsum(w3, xi, 3) sum(w3 .* xi, 3)

@test_approx_eq wsum(w1, x, 1) sum(w1 .* x, 1)
@test_approx_eq wsum(w2, x, 2) sum(w2 .* x, 2)
@test_approx_eq wsum(w3, x, 3) sum(w3 .* x, 3)

@test_approx_eq wsum(w1, Abs2Fun(), x, 1) sum(w1 .* abs2(x), 1)
@test_approx_eq wsum(w2, Abs2Fun(), x, 2) sum(w2 .* abs2(x), 2)
@test_approx_eq wsum(w3, Abs2Fun(), x, 3) sum(w3 .* abs2(x), 3)

@test_approx_eq wsum(w1, Multiply(), x, y, 1) sum(w1 .* (x .* y), 1)
@test_approx_eq wsum(w2, Multiply(), x, y, 2) sum(w2 .* (x .* y), 2)
@test_approx_eq wsum(w3, Multiply(), x, y, 3) sum(w3 .* (x .* y), 3)

@test_approx_eq wsumfdiff(w1, Abs2Fun(), x, y, 1) sum(w1 .* abs2(x - y), 1)
@test_approx_eq wsumfdiff(w2, Abs2Fun(), x, y, 2) sum(w2 .* abs2(x - y), 2)
@test_approx_eq wsumfdiff(w3, Abs2Fun(), x, y, 3) sum(w3 .* abs2(x - y), 3)

# convenience functions

x = randn(6)
y = rand(6)
w = rand(6)

@test_approx_eq wsumabs(w, x) sum(w .* abs(x))
@test_approx_eq wsumabsdiff(w, x, y) sum(w .* abs(x - y))
@test_approx_eq wsumabsdiff(w, x, 1) sum(w .* abs(x - 1))

@test_approx_eq wsumsq(w, x) sum(w .* abs2(x))
@test_approx_eq wsumsqdiff(w, x, y) sum(w .* abs2(x - y))
@test_approx_eq wsumsqdiff(w, x, 1) sum(w .* abs2(x - 1))

x = randn(5, 6)
y = rand(5, 6)
w1 = rand(5)
w2 = rand(1, 6)

@test_approx_eq wsumabs(w1, x, 1) sum(w1 .* abs(x), 1)
@test_approx_eq wsumabs(w2, x, 2) sum(w2 .* abs(x), 2)

r = zeros(6); wsumabs!(r, w1, x, 1)
@test_approx_eq r vec(wsumabs(w1, x, 1))
r = zeros(5); wsumabs!(r, w2, x, 2)
@test_approx_eq r vec(wsumabs(w2, x, 2))

@test_approx_eq wsumabsdiff(w1, x, y, 1) sum(w1 .* abs(x - y), 1)
@test_approx_eq wsumabsdiff(w2, x, y, 2) sum(w2 .* abs(x - y), 2)

r = zeros(6); wsumabsdiff!(r, w1, x, y, 1)
@test_approx_eq r vec(wsumabsdiff(w1, x, y, 1))
r = zeros(5); wsumabsdiff!(r, w2, x, y, 2)
@test_approx_eq r vec(wsumabsdiff(w2, x, y, 2))

@test_approx_eq wsumsq(w1, x, 1) sum(w1 .* abs2(x), 1)
@test_approx_eq wsumsq(w2, x, 2) sum(w2 .* abs2(x), 2)

r = zeros(6); wsumsq!(r, w1, x, 1)
@test_approx_eq r vec(wsumsq(w1, x, 1))
r = zeros(5); wsumsq!(r, w2, x, 2)
@test_approx_eq r vec(wsumsq(w2, x, 2))

@test_approx_eq wsumsqdiff(w1, x, y, 1) sum(w1 .* abs2(x - y), 1)
@test_approx_eq wsumsqdiff(w2, x, y, 2) sum(w2 .* abs2(x - y), 2)

r = zeros(6); wsumsqdiff!(r, w1, x, y, 1)
@test_approx_eq r vec(wsumsqdiff(w1, x, y, 1))
r = zeros(5); wsumsqdiff!(r, w2, x, y, 2)
@test_approx_eq r vec(wsumsqdiff(w2, x, y, 2))