text:
include("cpa.jl")
include("dpa.jl")
include("ta.jl")
using PyPlot
plt = PyPlot

function display_cpa(results)
  #to get the keymax, possible key strats from 0
  keymax = length(results.ρ[1,:,1]) - 1
  #to display the graphs for each key unit guessed
  for i in 1:length(results.keyGraphs)
    fig = plt.figure()
    plt.plot(results.keyGraphs[i], color="maroon")
    plt.xlim(0,keymax)
    xlabel("key byte values")
    ylabel("maximum correlation values")
    title("max correlation vs key values for unit$(i)")
    plt.savefig(join([results.res_dir, "/max_cc_keyvals_for_byte$(i)"]))
    plt.close(fig)
  end
  f = open(join([results.res_dir, "/log.txt"]),"w")
  write(f,join(["The recovered key is: ", results.key]))
  close(f)
  println(join(["The resulting cpa plots are saved at ", results.res_dir]))
  #can add more functions to manipulate ρ and show details
end

function display_dpa(results)
  num_unit = length(results.diff_of_means)
  for i in 1:num_unit
    dir = join([results.res_dir, "/byte$(i)"])
    mkdir(dir)
    diff_vecs = results.diff_of_means[i]
    for j in 1:length(diff_vecs)
      fig = plt.figure()
      plt.plot(diff_vecs[i])
      plt.xlim(0,length(diff_vecs[i]))
      xlabel("time")
      ylabel("difference mean power")
      t = "DPA subkey=$(j-1) for byte$(i)"
      title(t)
      plt.savefig(join([dir,"/subkey$(j-1).png"]))
      plt.close(fig)
    end
  end
  println(join(["The resulting dpa plots are saved at ", results.res_dir]))
end

function display_ta(results)
  num_unit = length(results.templates)
  templates = results.templates
  f = open(join([results.res_dir, "/log.txt"]), "w")
  write(f, join(["The recoverd key: ", results.recovered_key, "\n\n"]))

  for i in 1:num_unit
    template = templates[i]
    write(f,"_____________________\n")
    write(f,"byte$(template.this_unit)\n")
    write(f,"POI: ")
    for point in map(string,template.POI)
      write(f,join([point, " "]))
    end
    write(f,"\n\n")

    fig = plt.figure()
    plt.plot(template.tempSumDiff)
    xlabel("time")
    ylabel("sum of diff of means trace")
    title("sum of diff of means for byte=$(i)")
    plt.savefig(join([results.res_dir, "/sumdiffmean_$(i)"]))
    plt.close(fig)
  end
  close(f)
  println(join(["The results of TA are saved at ", results.res_dir]))
end

max_stars_repo_name: udion/zypher

max_issues_repo_name: udion/zypher

max_forks_repo_name: udion/zypher


text:
"""
	This is a script alternative to the julia_get_proteins_under_200_aa.ipynb
	notebook for those who don't have IJulia or would like to run it as a script

	Notebook to preprocess the raw data file and
	handle it properly.
	Will prune the unnecessary data for now.
	Reducing data file from 600mb to 170mb.

	Select only proteins under L aminoacids (AAs).
"""

L = 200										# Set maximum AA length
N = 995										# Set maximum number of proteins
RAW_DATA_PATH = "../data/training_30.txt"	# Path to raw data file
DESTIN_PATH = "../data/full_under_200.txt"	# Path to destin file

# Open the file and read content

f = try open(RAW_DATA_PATH) catch
	println("File not found. Check it's there. Instructions in the readme.")
	exit(0)
	end
lines = readlines(f)

function coords_split(lister, splice)
    # Split all passed sequences by "splice" and return an array of them
    # Convert string fragments to float
    coords = []
    for c in lister
        push!(coords, [parse(Float64, a) for a in split(c, splice)])
    end
    return coords
end

function norm(vector)
	# Could use "Using LinearAlgebra + built-in norm()" but gotta learn Julia
    return sqrt(sum([v*v for v in vector]))
end

# Scan first n proteins
names = []
seqs = []
coords = []
pssms = []

try
    # Record names, seqs and coords for each protein btwn 1-n
    for i in 1:length(lines)
        if length(coords) == N
            break
        end

        # Start recording
        if lines[i] == "[ID]"
            push!(names, lines[i+1])
        elseif lines[i] == "[PRIMARY]"
            push!(seqs, lines[i+1])
        elseif lines[i] == "[TERTIARY]"
            push!(coords, coords_split(lines[i+1:i+3], "\t"))
        elseif lines[i] == "[EVOLUTIONARY]"
            push!(pssms, coords_split(lines[i+1:i+21], "\t"))
            # Progress control
            if length(names)%50 == 0
                println("Currently @ ", length(names), " out of n: ", N)
            end
        end
    end
catch
    println("Error while reading file. Check it's complete or download again.")
    exit(0)
end

# Check proteins w/ length under L
println("\n\nTotal number of proteins: ", length(seqs))
under = []
for i in 1:length(seqs)
    if length(seqs[i])<L
        push!(under, i)
        # Uncomment for debugging purposes
        # println("Seelected with: ", length(seqs[i]), " number: ", i)
    end
end
println("Number of proteins under ", L, " : ", length(under), "\n\n")

# Get distances btwn pairs of AAs - only for prots under 200
dists = []
try
    for k in under
        # Get distances from coordinates
        dist = []
        for i in 1:length(coords[k][1])
            # Only pick coords for C-alpha carbons! - position (1/3 of total data)
            # i%3 == 2 Because juia arrays start at 1 - Python: i%3 == 1
            if i%3 == 2
                aad = [] # Distance to every AA from a given AA
                for j in 1:length(coords[k][1])
                    if j%3 == 2
                        push!(aad, norm([coords[k][1][i],coords[k][2][i],coords[k][3][i]]
                    	    			-
                    		    		[coords[k][1][j],coords[k][2][j],coords[k][3][j]]))
                    end
                end
                push!(dist, aad)
            end
        end
        push!(dists, dist)

        # Progress control
        if length(dists)%50 == 0
            println("Dists Currently @ ", length(dists), " out of n: ", N)
        end
    end
catch
	println("Error while calculating distances. Set N to smaller value.")
end

# Check everything's alright
n = 2
println("\n\nSample protein data (example)")
println("id: ", names[n])
println("seq: ", seqs[n])
println("sample coord: ", coords[n][1][1])
println("sample dist: ", dists[n][1][5])

# Data is OK. Save it to a file.
using DelimitedFiles
open(DESTIN_PATH, "w") do f
    aux = [0]
    for k in under
        push!(aux, aux[length(aux)]+1)
        # ID
        write(f, "\n[ID]\n")
        write(f, names[k])
        # Seq
        write(f, "\n[PRIMARY]\n")
        write(f, seqs[k])
        # PSSMS
        write(f, "\n[EVOLUTIONARY]\n")
        writedlm(f, pssms[k])
        # Coords
        write(f, "\n[TERTIARY]\n")
        writedlm(f, coords[k])
        # Dists
        write(f, "\n[DIST]\n")
        # Check that saved proteins are less than 200 AAs
        if length(dists[aux[length(aux)]][1])>L
            println("error when checking protein in dists n: ",
            		 aux[length(aux)], " length: ", length(dists[aux[length(aux)]][1]))
            break
        else
            writedlm(f, dists[aux[length(aux)]])
        end
    end
end

println("\n\nScript execution went fine. Data is ready at: ", DESTIN_PATH)
exit(0)

ext: jl

max_stars_repo_name: entruv/MiniFold

max_issues_repo_name: entruv/MiniFold

max_forks_repo_name: entruv/MiniFold


text:
module ThreadLocalCountersTests

include("test_counters.jl")

end  # module ThreadLocalCountersTests


text:
using QContinuousMeasurement

using Test
using Random

@testset "Dicke basis" begin
    @testset "ω = $ω" for ω in range(0,stop=1.,length=3)
        @testset "Nj = $Nj" for Nj in [1, 4, 5]
            Ntraj = 1
            Tfinal = 1.0
            dt = 0.001
            κ_ind = 1.
            κ_coll = 1.
            θ = 0. # Angle of the independent noise (θ = 0 : parallel)
            η = 1.0

            seed = 20

            Random.seed!(seed)
            @time res_dicke = Eff_QFI_HD_Dicke(Nj, Ntraj,# Number of trajectories
                Tfinal,                              # Final time
                dt;                                  # Time step
                κ = κ_ind,
                κcoll = κ_coll,
                ω = ω,
                η = η)

            Random.seed!(seed)
            @time res_sup = Eff_QFI_HD(Nj, Ntraj,# Number of trajectories
                            Tfinal,                              # Final time
                            dt;                                  # Time step
                            κ = κ_ind,
                            κcoll = κ_coll,
                            ω = ω,
                            η = η)      # Initial state

            er_FI = maximum(abs.(res_dicke.FI - res_sup.FI) ./ res_sup.FI)

            #println(er_FI)
            if er_FI > 1e-8
                @warn "Relative error FI high" er_FI
            end

            er_QFI = maximum(abs.(res_dicke.QFI - res_sup.QFI) ./ res_sup.QFI)
            if er_QFI > 1e-8
                @warn "Relative error QFI high" er_QFI
            end

            #println(er_QFI)

            rtol = 1e-6
            @test res_dicke.FI ≈ res_sup.FI rtol=rtol atol=dt^2
            @test res_dicke.QFI ≈ res_sup.QFI rtol=rtol atol=dt^2
            @test res_dicke.jx ≈ res_sup.jx rtol=rtol atol=dt^2
            @test res_dicke.jy ≈ res_sup.jy rtol=rtol atol=dt^2
            @test res_dicke.jz ≈ res_sup.jz rtol=rtol atol=dt^2
        end
    end
end


text:
"""
    densify(S::SparseMatrix, threshold=0.3)

Convert sparse `S` to `Diagonal` if `S` is diagonal or to `full(S)` if
the proportion of nonzeros exceeds `threshold`.
"""
function densify(A::SparseMatrixCSC, threshold::Real = 0.3)
    S = sparse(A)
    m, n = size(S)
    if m == n && isdiag(S)  # convert diagonal sparse to Diagonal
        Diagonal(diag(S))
    elseif nnz(S)/(S.m * S.n) ≤ threshold
        A
    else
        Array(S)
    end
end
densify(A::AbstractMatrix, threshold::Real = 0.3) = A
function densify(A::BlockedSparse, threshold::Real=0.3)
    Asp = A.cscmat
    Ad = densify(Asp)
    Ad === Asp ? A : Ad
end

function LinearMixedModel(f, trms, wts)
    n = size(trms[1], 1)
    T = eltype(trms[end])
    sqrtwts = sqrt.(wts)
    sz = size.(trms, 2)
    if !isempty(wts)
        reweight!.(trms, Vector[sqrtwts])
    end
    nt = length(trms)
    A = BlockArrays._BlockArray(AbstractMatrix{T}, sz, sz)
    L = BlockArrays._BlockArray(AbstractMatrix{T}, sz, sz)
    for j in 1:nt, i in j:nt
        Lij = L[Block(i,j)] = densify(trms[i]'trms[j])
        A[Block(i,j)] = deepcopy(isa(Lij, BlockedSparse) ? Lij.cscmat : Lij)
    end
                  # check for fill-in due to non-nested grouping factors
    for i in 2:nt
        ti = trms[i]
        if isa(ti, AbstractFactorReTerm)
            for j in 1:(i - 1)
                tj = trms[j]
                if isa(tj, AbstractFactorReTerm) && !isnested(tj, ti)
                    for k in i:nt
                        L[Block(k, i)] = Matrix(L[Block(k, i)])
                    end
                    break
                end
            end
        end
    end
    optsum = OptSummary(getθ(trms), lowerbd(trms), :LN_BOBYQA;
        ftol_rel = T(1.0e-12), ftol_abs = T(1.0e-8))
    fill!(optsum.xtol_abs, 1.0e-10)
    LinearMixedModel(f, trms, sqrtwts, A, LowerTriangular(L), optsum)
end

function LinearMixedModel(f::Formula, fr::AbstractDataFrame;
    weights::Vector = [], contrasts=Dict(), rdist::Distribution=Normal())
    mf = ModelFrame(f, fr, contrasts=contrasts)
    X = ModelMatrix(mf).m
    n = size(X, 1)
    T = eltype(X)
    y = model_response(mf, rdist)
    tdict = Dict{Symbol, Vector{Any}}()
    for t in filter(x -> Meta.isexpr(x, :call) && x.args[1] == :|, mf.terms.terms)
        fnm = t.args[3]
        isa(fnm, Symbol) || throw(ArgumentError("rhs of $t must be a symbol"))
        tdict[fnm] = haskey(tdict, fnm) ? push!(tdict[fnm], t.args[2]) : [t.args[2]]
    end
    isempty(tdict) && throw(ArgumentError("No random-effects terms found in $f"))
    trms = AbstractTerm{T}[]
    for (grp, lhs) in tdict
        gr = compress(categorical(getindex(mf.df, grp)))
        if (length(lhs) == 1 && lhs[1] == 1)
            push!(trms, ScalarFactorReTerm(gr, grp))
        else
            m = T[]
            coefnms = String[]
            trsize = Int[]
            for l in lhs
                if l == 1
                    append!(m, ones(T, n))
                    push!(coefnms, "(Intercept)")
                    push!(trsize, 1)
                else
                    fr = ModelFrame(@eval(@formula($nothing ~ $l)), mf.df)
                    append!(m, ModelMatrix(fr).m)
                    cnms = coefnames(fr)
                    append!(coefnms, cnms)
                    push!(trsize, length(cnms))
                end
            end
            push!(trms,
                  length(coefnms) == 1 ? ScalarFactorReTerm(gr, m, grp, coefnms) :
                  VectorFactorReTerm(gr, collect(adjoint(reshape(m, (n, sum(trsize))))), grp,
                      coefnms,  trsize))
        end
    end
    sort!(trms, by = nrandomeff, rev = true)
    push!(trms, MatrixTerm(X, coefnames(mf)))
    push!(trms, MatrixTerm(y))
    LinearMixedModel(f, trms, oftype(y, weights))
end

StatsBase.model_response(mf::ModelFrame, d::Distribution) =
    model_response(mf.df[mf.terms.eterms[1]], d)

StatsBase.model_response(v::AbstractVector, d::Distribution) = Vector{partype(d)}(v)

function StatsBase.model_response(v::CategoricalVector, d::Bernoulli)
    levs = levels(v)
    nlevs = length(levs)
    @argcheck(nlevs ≤ 2)
    nlevs < 2 ? zeros(v, partype(d)) : partype(d)[cv == levs[2] for cv in v]
end

StatsBase.fit(::Type{LinearMixedModel},
              f::Formula,
              fr::AbstractDataFrame;
              weights=[],
              contrasts=Dict(),
              verbose=false,
              REML=false) =
    fit!(LinearMixedModel(f, fr, weights=weights, contrasts=contrasts),
         verbose=verbose,
         REML=REML)

"""
    updateL!(m::LinearMixedModel)

Update the blocked lower Cholesky factor, `m.L`, from `m.A` and `m.trms` (used for Λ only)

This is the crucial step in evaluating the objective, given a new parameter value.
"""
function updateL!(m::LinearMixedModel{T}) where {T}
    A = m.A
    Ldat = m.L.data
    trms = m.trms
    nblk = length(trms)
    for (j, trm) in enumerate(trms)
        Ljj = scaleInflate!(Ldat[Block(j, j)], A[Block(j, j)], trm)
        LjjH = isa(Ljj, Diagonal) ? Ljj : Hermitian(Ljj, :L)
        LjjT = isa(Ljj, Diagonal) ? Ljj : LowerTriangular(Ljj)
        for jj in 1:(j - 1)
            rankUpdate!(LjjH, Ldat[Block(j, jj)], -one(T))
        end
        cholUnblocked!(Ljj, Val{:L})
        for i in (j + 1):nblk
            Lij = lmul!(Λ(trms[i])', rmul!(copyto!(Ldat[Block(i, j)], A[Block(i, j)]), Λ(trm)))
            for jj in 1:(j - 1)
                mulαβ!(Lij, Ldat[Block(i, jj)], Ldat[Block(j, jj)]', -one(T), one(T))
            end
            rdiv!(Lij, LjjT')
        end
    end
    m
end

StatsBase.coef(m::MixedModel) = fixef(m, false)

"""
    fit!(m::LinearMixedModel; verbose=false, REML=false)

Optimize the objective of a `LinearMixedModel`.  When `verbose` is `true` the values of the
objective and the parameters are printed on stdout at each function evaluation.  The
objective is negative twice the log-likelihood when `REML` is `false` (the default) or
the REML criterion.
"""
function StatsBase.fit!(m::LinearMixedModel{T}; verbose=false, REML=nothing) where {T}
    optsum = m.optsum
    opt = Opt(optsum)
    feval = 0
    if isa(REML, Bool)
        optsum.REML = REML
    end
    function obj(x, g)
        isempty(g) || error("gradient not defined")
        feval += 1
        val = objective(updateL!(setθ!(m, x)))
        feval == 1 && (optsum.finitial = val)
        verbose && println("f_", feval, ": ", round(val, digits=5), " ", x)
        val
    end
    NLopt.min_objective!(opt, obj)
    fmin, xmin, ret = NLopt.optimize!(opt, copyto!(optsum.final, optsum.initial))
    ## check if small non-negative parameter values can be set to zero
    xmin_ = copy(xmin)
    lb = optsum.lowerbd
    for i in eachindex(xmin_)
        if iszero(lb[i]) && zero(T) < xmin_[i] < T(0.001)
            xmin_[i] = zero(T)
        end
    end
    if xmin_ ≠ xmin
        if (zeroobj = obj(xmin_, T[])) ≤ (fmin + 1.e-5)
            fmin = zeroobj
            copyto!(xmin, xmin_)
        end
    end
    ## ensure that the parameter values saved in m are xmin
    updateL!(setθ!(m, xmin))

    optsum.feval = feval
    optsum.final = xmin
    optsum.fmin = fmin
    optsum.returnvalue = ret
    ret == :ROUNDOFF_LIMITED && @warn("NLopt was roundoff limited")
    if ret ∈ [:FAILURE, :INVALID_ARGS, :OUT_OF_MEMORY, :FORCED_STOP, :MAXFEVAL_REACHED]
        @warn("NLopt optimization failure: $ret")
    end
    m
end

function fitted!(v::AbstractArray{T}, m::LinearMixedModel{T}) where {T}
    ## FIXME: Create and use `effects(m) -> β, b` w/o calculating β twice
    trms = m.trms
    mul!(vec(v), trms[end - 1], fixef(m))
    b = ranef(m)
    for j in eachindex(b)
        unscaledre!(vec(v), trms[j], b[j])
    end
    v
end

StatsBase.fitted(m::LinearMixedModel{T}) where {T} = fitted!(Vector{T}(undef, nobs(m)), m)

StatsBase.predict(m::LinearMixedModel) = fitted(m)

StatsBase.residuals(m::LinearMixedModel) = model_response(m) .- fitted(m)

"""
    lowerbd(m::LinearMixedModel)

Return the vector of lower bounds on the covariance parameter vector `θ`
"""
lowerbd(m::LinearMixedModel) = lowerbd(m.trms)

StatsBase.model_response(m::LinearMixedModel) = vec(m.trms[end].x)

"""
    objective(m::LinearMixedModel)

Return negative twice the log-likelihood of model `m` or the REML criterion,
according to the value of `m.optsum.REML`
"""
function objective(m::LinearMixedModel)
    wts = m.sqrtwts
    logdet(m) + varest_denom(m)*(1 + log2π + log(varest(m))) -
        (isempty(wts) ? 0 : 2sum(log, wts))
end

"""
    fixef!(v::Vector{T}, m::LinearMixedModel{T})

Overwrite `v` with the pivoted and, possibly, truncated fixed-effects coefficients of model `m`
"""
function fixef!(v::AbstractVector{T}, m::LinearMixedModel{T}) where {T}
    L = feL(m)
    @argcheck(length(v) == size(L, 1), DimensionMismatch)
    ldiv!(adjoint(L), copyto!(v, m.L.data.blocks[end, end - 1]))
end

"""
    fixef(m::MixedModel, permuted=true)

Return the fixed-effects parameter vector estimate of `m`.

If `permuted` is `true` the vector elements are permuted according to
`m.trms[end - 1].piv` and truncated to the rank of that term.
"""
function fixef(m::LinearMixedModel{T}, permuted=true) where {T}
    permuted && return fixef!(Vector{T}(undef, size(m)[2]), m)
    Xtrm = m.trms[end - 1]
    piv = Xtrm.piv
    v = fill(-zero(T), size(piv))
    fixef!(view(v, 1:Xtrm.rank), m)
    invpermute!(v, piv)
end

StatsBase.dof(m::LinearMixedModel) = size(m)[2] + sum(nθ, m.trms) + 1

function StatsBase.dof_residual(m::LinearMixedModel)
    (n, p, q, k) = size(m)
    n - m.optsum.REML * p
end

StatsBase.loglikelihood(m::LinearMixedModel) = -objective(m)/2

StatsBase.nobs(m::LinearMixedModel) = length(m.trms[end].wtx)

function Base.size(m::LinearMixedModel)
    trms = m.trms
    k = length(trms) - 2
    q = sum(size(trms[j], 2) for j in 1:k)
    length(trms[end]), trms[k+1].rank, q, k
end

"""
    sdest(m::LinearMixedModel)

Return the estimate of σ, the standard deviation of the per-observation noise.
"""
sdest(m::LinearMixedModel) = √varest(m)

"""
    setθ!{T}(m::LinearMixedModel{T}, v::Vector{T})

Install `v` as the θ parameters in `m`.
"""
function setθ!(m::LinearMixedModel, v)
    offset = 0
    for trm in m.trms
        if isa(trm, ScalarFactorReTerm)
            offset += 1
            setθ!(trm, v[offset])
        else
            k = nθ(trm)
            iszero(k) || setθ!(trm, view(v, (1:k) .+ offset))
            offset += k
        end
    end
    m
end

"""
    sqrtpwrss(m::LinearMixedModel)

Return the square root of the penalized, weighted residual sum-of-squares (pwrss).

This value is the contents of the `1 × 1` bottom right block of `m.L`
"""
sqrtpwrss(m::LinearMixedModel) = @views m.L.data.blocks[end, end][1]

function varest_denom(m::LinearMixedModel)
    (n, p, q, k) = size(m)
    n - m.optsum.REML * p
end

"""
    varest(m::LinearMixedModel)

Returns the estimate of σ², the variance of the conditional distribution of Y given B.
"""
varest(m::LinearMixedModel) = pwrss(m) / varest_denom(m)

"""
    pwrss(m::LinearMixedModel)

The penalized residual sum-of-squares.
"""
pwrss(m::LinearMixedModel) = abs2(sqrtpwrss(m))

function StatsBase.deviance(m::LinearMixedModel)
    isfit(m) || error("Model has not been fit")
    objective(m)
end

"""
    isfit(m::LinearMixedModel)

Return a `Bool` indicating if the model has been fit.
"""
isfit(m::LinearMixedModel) = m.optsum.fmin < Inf

"""
    lrt(mods::LinearMixedModel...)

Perform sequential likelihood ratio tests on a sequence of models.

The returned value is a `DataFrame` containing information on the likelihood ratio tests.
"""
function lrt(mods::LinearMixedModel...) # not tested
    (nm = length(mods)) > 1 || throw(ArgumentError("at least two models are required for a likelihood ratio test"))
    m1 = mods[1]
    n = nobs(m1)
    for i in 2:nm
        nobs(mods[i]) == n || throw(ArgumentError("number of observations must be constant across models"))
    end
    mods = mods[sortperm([dof(m)::Int for m in mods])]
    degf = Int[dof(m) for m in mods]
    dev = [deviance(m)::Float64 for m in mods]
    csqr = pushfirst!([(dev[i-1]-dev[i])::Float64 for i in 2:nm],NaN)
    pval = pushfirst!([ccdf(Chisq(degf[i]-degf[i-1]),csqr[i])::Float64 for i in 2:nm],NaN)
    DataFrame(Df = degf, Deviance = dev, Chisq=csqr,pval=pval)
end

"""
    reweight!{T}(m::LinearMixedModel{T}, wts::Vector{T})

Update `m.sqrtwts` from `wts` and reweight each term.  Recompute `m.A` and `m.L`.
"""
function reweight!(m::LinearMixedModel, weights)
    trms = m.trms
    m.sqrtwts .= sqrt.(weights)
    reweight!.(trms, Vector[m.sqrtwts])
    ntrm = length(trms)
    A = m.A
    for j in 1:ntrm, i in j:ntrm
        mul!(A[Block(i, j)], trms[i]', trms[j])
    end
    updateL!(m)
end

function Base.show(io::IO, m::LinearMixedModel)
    if !isfit(m)
        @warn("Model has not been fit")
        return nothing
    end
    n, p, q, k = size(m)
    REML = m.optsum.REML
    println(io, "Linear mixed model fit by ", REML ? "REML" : "maximum likelihood")
    println(io, " ", m.formula)
    oo = objective(m)
    if REML
        println(io, " REML criterion at convergence: ", oo)
    else
        nums = showoff([-oo/ 2, oo, aic(m), bic(m)])
        fieldwd = max(maximum(textwidth.(nums)) + 1, 11)
        for label in [" logLik", "-2 logLik", "AIC", "BIC"]
            print(io, rpad(lpad(label, (fieldwd + textwidth(label)) >> 1), fieldwd))
        end
        println(io)
        for num in nums
            print(io, lpad(num, fieldwd))
        end
        println(io)
    end
    println(io)

    show(io,VarCorr(m))

    gl = grplevels(m)
    @printf(io," Number of obs: %d; levels of grouping factors: %d", n, gl[1])
    for l in gl[2:end] @printf(io, ", %d", l) end
    println(io)
    @printf(io,"\n  Fixed-effects parameters:\n")
    show(io,coeftable(m))
end


text:
@testset "NumFieldOrd" begin
  x = PolynomialRing(FlintQQ, "x", cached = false)[2]
  K, a = number_field(x^2+1, check = false, cached = false)
  Kns, gKns = number_field([x^2+1], check = false, cached = false)
  Kt, t = PolynomialRing(K, "t", cached = false)
  L, b = number_field(t^2+3, "b", cached = false, check = false)
  Lns, gLns = number_field([t^2+3], cached = false, check = false)

  OK = maximal_order(K)
  @inferred nf(OK)
  @test is_simple(OK)
  @test is_commutative(OK)
  @test @inferred is_equation_order(OK)
  @test @inferred is_maximal(OK)
  @test @inferred degree(OK) == 2
  @test @inferred absolute_degree(OK) == 2
  BOK = absolute_basis(OK)
  @test @inferred discriminant(BOK) == @inferred discriminant(OK)
  @test @inferred signature(OK) == (0, 1)
  @test @inferred discriminant(OK, FlintQQ) == discriminant(OK)

  OKns = maximal_order(Kns)
  @inferred nf(OKns)
  @test !is_simple(OKns)
  @test is_commutative(OKns)
  @test @inferred !is_equation_order(OKns)
  @test @inferred is_maximal(OKns)
  @test @inferred degree(OKns) == 2
  @test @inferred absolute_degree(OKns) == 2
  BOKns = absolute_basis(OKns)
  @test @inferred discriminant(BOKns) == @inferred discriminant(OKns)
  @test @inferred signature(OKns) == (0, 1)
  @test @inferred discriminant(OKns, FlintQQ) == discriminant(OKns)

  OL = maximal_order(L)
  @inferred nf(OL)
  @test is_simple(OL)
  @test is_commutative(OL)
  @test @inferred !is_equation_order(OL)
  @test @inferred is_maximal(OL)
  @test @inferred degree(OL) == 2
  @test @inferred absolute_degree(OL) == 4
  BOL = absolute_basis(OL)
  @test numerator(det(trace_matrix(map(elem_in_nf, BOL)))) == @inferred absolute_discriminant(OL)
  @test @inferred signature(OL) == (0, 2)
  @test @inferred discriminant(OL, FlintQQ) == absolute_discriminant(OL)

  OLns = maximal_order(Lns)
  @inferred nf(OLns)
  @test !is_simple(OLns)
  @test is_commutative(OLns)
  @test @inferred !is_equation_order(OLns)
  @test @inferred is_maximal(OLns)
  @test @inferred degree(OLns) == 2
  @test @inferred absolute_degree(OLns) == 4
  BOLns = absolute_basis(OLns)
  @test numerator(det(trace_matrix(map(elem_in_nf, BOLns)))) == @inferred absolute_discriminant(OLns)
  @test @inferred signature(OLns) == (0, 2)
  @test @inferred discriminant(OLns, FlintQQ) == absolute_discriminant(OLns)

end

@testset "RelResidueField" begin
  Qs, s = QQ["s"]
  K,a = number_field(s^4+s+1, "a")
  Kt, t = K["t"]
  E, b = number_field(t^2-a*t+1,"b")
  Eu, u = E["u"]
  L, c = number_field(u^3-7, "c")
  OK = maximal_order(K)
  OE = maximal_order(E)
  OL = maximal_order(L)
  p = prime_decomposition(OK, 11)[1][1]
  P = prime_decomposition(OE,p)[1][1]
  PP = prime_decomposition(OL, P)[1][1]
  FL, projL = relative_residue_field(OL, PP)
  mL = extend(projL, L)

  e,f = PP.splitting_type
  @test degree(defining_polynomial(FL)) == f
  e,f = P.splitting_type
  @test degree(defining_polynomial(base_field(FL))) == f
  e, f = p.splitting_type
  @test degree(defining_polynomial(base_field(base_field(FL)))) == f

  @test domain(projL.map_subfield) === OE
  @test mL(gen(L)//2) == gen(FL)//2
  @test degree(Hecke.absolute_field(FL)[1]) == 6
  @test characteristic(Hecke.prime_field(FL)) == 11

  K,a = rationals_as_number_field()
  Kt, t = K["t"]
  E,b = number_field(t^4-5*t^3-6*t^2+5*t+1,"b")
  OK = maximal_order(K)
  OE = maximal_order(E)
  p = prime_decomposition(OK, 2)[1][1]
  P = prime_decomposition(OE, p)[1][1]
  @test is_index_divisor(OE, p)
  FE, projE = relative_residue_field(OE, P)
  _, f = P.splitting_type
  @test degree(defining_polynomial(FE)) == f
  mE = extend(projE, E)
  @test !iszero(mE(gen(E)))
  @test get_attribute(P, :rel_residue_field_map) !== nothing

end


text:
using HypergeometricFunctions
using SpecialFunctions
using StaticArrays # To have access to static arrays
using Interpolations # To have access to interpolation functions

##################################################
# Distribution function in (E,L)
##################################################

"""
    _H(x,[q=qCalc])

Computes the `H(x,q)` function used in the distribution function.

# Remarks
- Uses pre-computed interpolated functions when using `₂F₁`.

# Arguments
- `x::Float64`: `x` parameter.
- `q ::Float64`: Anisotropy parameter.
"""
function _H(x::Float64, q::Float64=qCalc)
    if (x <= 1)
        pref = 1/(GAMMA_ca)
        HG   = H_1(x)
        return pref*HG
    else
        pref = 1/(GAMMA_db*GAMMA_bc)
        HG   = H_3(1/x)
        return pref*x^(-q/2)*HG
    end
end

"""
    _tF(tE,tL)

Non-dimensional distribution function of a Plummer sphere.

# Arguments
- `tE::Float64`: Reduced energy `tE = E/E0`.
- `tL::Float64`: Reduced energy `tL = L/L0`.
"""
function _tF(tE::Float64, tL::Float64)

    if (tE < 0.0 || tL < 0.0) # If E or L are negative, the DF vanishes
        return 0.0
    else # If E and L are positive

        if (qCalc == 0.0) # Isotropic case
             return 3.0/(7.0*PI^3) * (2.0*tE)^(7/2)
        elseif (qCalc == 2.0) # Anisotropic case
            x = tL^2/(2*tE)
            if (x <= 1)
                return 6.0/(2.0*PI)^3 * (2.0*tE - tL^2)^(3/2)
            else
                return 0.0
            end
        else # Case q<2
            x = tL^2/(2*tE)
            return (3.0*GAMMA_6q/(2.0*(2.0*PI)^(5/2)) *
                       tE^(7/2-qCalc) * _H(x,qCalc))
        end

    end
end

"""
    _F(E,L)

Distribution function of a Plummer sphere.
Normalized to the total mass of the cluster `_M`.

# Arguments
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
"""
function _F(E::Float64, L::Float64)
    tE = _tE(E)
    tL = _tL(L)
    DF  = _tF(tE,tL)
    return _M*_F0*DF
end

##################################################
# Distribution function in (E,L)
# For flux computation
##################################################

"""
    _HdH(x,[q=qCalc])

Computes the `H(x,q)` function used in the distribution function and its derivative.
Returns a tuple `(H, dH/dx)`.

# Remarks
- Uses pre-computed interpolated functions when using `₂F₁`.

# Arguments
- `x::Float64`: `x` parameter.
- `q ::Float64`: Anisotropy parameter.
"""
function _HdH(x::Float64, q::Float64=qCalc)
    if (x <= 1)
        pref = 1/(GAMMA_ca)
        HG   = H_1(x)
        HGp  = H_2(x)
        H = pref*HG
        dH = pref*x^(-1)*((q/2)*(q-3.5)*x*HGp)
        return H, dH
    else
        pref = 1/(GAMMA_db*GAMMA_bc)
        HG   = H_3(1/x)
        HGp  = H_4(1/x)
        H = pref*x^(-q/2)*HG
        dH = pref*x^(-q/2-1)*((-q/2)*HG - (1/x)*HGp*(q/2)*(q/2)/(4.5-q/2))
        return H, dH
    end
end

"""
    _tFdF(tE,tL)

Non-dimensional distribution function of a Plummer sphere and its gradient.
Returns a tuple `(F, dtF/dtE, dtF/dtL)`.

# Arguments
- `tE::Float64`: Reduced energy `tE = E/E0`.
- `tL::Float64`: Reduced energy `tL = L/L0`.
"""
function _tFdF(tE::Float64, tL::Float64, q::Float64=qCalc)
    if (tE <= 0.0 || tL <= 0.0) # If E or L are negative, the DF vanishes
        DF = 0.0
        dE = 0.0
        dL = 0.0
        return DF, dE, dL
    else
        x = tL^2/(2*tE)
        if (q == 0.0) # Isotropic case
            pref = 3/(PI^3) * (2*tE)^(5/2)
            DF = pref/7 * (2*tE)
            dE = pref
            dL = 0.0
            return DF, dE, dL
        else
            if (q == 2.0)
                if (x <= 1.0)
                    pref = 18/(2*PI)^3*(2*tE-tL^2)^(1/2)
                    DF = (1/3)*pref* (2*tE - tL^2)
                    dE = pref
                    dL = -pref*tL
                    return DF, dE, dL
                else
                    DF = 0.0
                    dE = 0.0
                    dL = 0.0
                    return DF, dE, dL
                end
            else
                pref = 3*GAMMA_6q/(2*(2*PI)^(5/2))*tE^(5/2-q)
                pref0 = pref*tE
                pref1 = pref
                H, dH = _HdH(x,q)

                DF = pref0 * H
                dE = pref1*((7/2-q)*H - x*dH)
                dL = pref1*(tL*dH)
                return DF, dE, dL
            end
        end
    end
end

"""
    _FdF(E,L)

Distribution function of a Plummer sphere and its gradient.
Normalized to the total mass of the cluster `_M`.
Returns a tuple `(F, dF/dE, dF/dL)`.

# Arguments
- `E::Float64`: Energy.
- `L::Float64`: Angular momentum.
"""
function _FdF(E::Float64, L::Float64, q::Float64=qCalc)
    tE = _tE(E)
    tL = _tL(L)
    tDF, tdE, tdL = _tFdF(tE,tL,q)
    DF = _M*_F0*tDF
    dE = (_M*_F0/_E0)*tdE
    dL = (_M*_F0/_L0)*tdL
    return DF, dE, dL
end

##################################################
# Interpolate the hypergeometric function
# HG1   = _₂F₁(qCalc/2,qCalc-3.5,1.0,x)
# HGp1  = _₂F₁(qCalc/2+1,qCalc-2.5,2.0,x)
# HG2   = _₂F₁(qCalc/2,qCalc/2,4.5-qCalc/2,1/x)
# HGp2  = _₂F₁(qCalc/2+1,qCalc/2+1,5.5-qCalc/2,1/x)

# valeur en x=1: https://homepage.tudelft.nl/11r49/documents/wi4006/hyper.pdf equation 4
# 2F1(a,b,c,1) = gamma(c)gamma(c-a-b)/(gamma(c-a)gamma(c-b))

# Valeur en x=1
# HG1   = _₂F₁(qCalc/2,qCalc-3.5,1.0,1) = gamma(4.5-1.5*qCalc)/(gamma(1-qCalc/2)gamma(4.5-qCalc))
# HGp1  = _₂F₁(qCalc/2+1,qCalc-2.5,2.0,1) = gamma(3.5-1.5*qCalc)/(gamma(1-qCalc/2)gamma(4.5-qCalc))
# HG2   = _₂F₁(qCalc/2,qCalc/2,4.5-qCalc/2,1) = gamma(4.5-qCalc/2)gamma(4.5-1.5*qCalc)/(gamma(4.5-qCalc)gamma(4.5-qCalc))
# HGp2  = _₂F₁(qCalc/2+1,qCalc/2+1,5.5-qCalc/2,1) = gamma(5.5-qCalc/2)gamma(3.5-1.5*qCalc)/(gamma(4.5-qCalc)gamma(4.5-qCalc))
##################################################

"""
    getHGInt(nbxInt=1000)

Interpolate the following hypergeometric function :
- `HG1   = _₂F₁(qCalc/2,qCalc-3.5,1.0,x)`
- `HGp1  = _₂F₁(qCalc/2+1,qCalc-2.5,2.0,x)`
- `HG2   = _₂F₁(qCalc/2,qCalc/2,4.5-qCalc/2,1/x)`
- `HGp2  = _₂F₁(qCalc/2+1,qCalc/2+1,5.5-qCalc/2,1/x)`

# Remarks
- Manually sets the values at x=0 and x=1.
- `₂F₁(a,b,c,0) = 1.0` .
- `₂F₁(a,b,c,1) = gamma(c)gamma(c-a-b)/(gamma(c-a)gamma(c-b))`.

# Arguments
- `nbxInt::Int64`: Number of interpolation sampling points.
"""
function getHGInt(nbxInt::Int64=1000)
    xminInt = 0.0
    xmaxInt = 1.0
    #####
    rangexInt = range(xminInt,length=nbxInt,xmaxInt)
    tabxInt = collect(rangexInt)
    tabHG1Int = zeros(Float64,nbxInt)
    tabHGp1Int = zeros(Float64,nbxInt)
    tabHG2Int = zeros(Float64,nbxInt)
    tabHGp2Int = zeros(Float64,nbxInt)

    #####
    for indx=2:nbxInt-1
        xloc = tabxInt[indx]
        hg1loc = _₂F₁(qCalc/2,qCalc-3.5,1.0,xloc)
        hgp1loc = _₂F₁(qCalc/2+1,qCalc-2.5,2.0,xloc)
        hg2loc = _₂F₁(qCalc/2,qCalc/2,4.5-qCalc/2,xloc)
        hgp2loc = _₂F₁(qCalc/2+1,qCalc/2+1,5.5-qCalc/2,xloc)
        tabHG1Int[indx] = hg1loc
        tabHGp1Int[indx] = hgp1loc
        tabHG2Int[indx] = hg2loc
        tabHGp2Int[indx] = hgp2loc
    end
    # x=0
    tabHG1Int[1] = 1.0
    tabHGp1Int[1] = 1.0
    tabHG2Int[1] = 1.0
    tabHGp2Int[1] = 1.0

    # x=1
    tabHG1Int[nbxInt] = gamma(4.5-1.5*qCalc)/(gamma(1-qCalc/2)*gamma(4.5-qCalc))
    tabHGp1Int[nbxInt] = gamma(3.5-1.5*qCalc)/(gamma(1-qCalc/2)*gamma(4.5-qCalc))
    tabHG2Int[nbxInt] = gamma(4.5-qCalc/2)*gamma(4.5-1.5*qCalc)/(gamma(4.5-qCalc)*gamma(4.5-qCalc))
    tabHGp2Int[nbxInt] = gamma(5.5-qCalc/2)*gamma(3.5-1.5*qCalc)/(gamma(4.5-qCalc)*gamma(4.5-qCalc))
    #####
    intHG1 = Interpolations.scale(interpolate(tabHG1Int, BSpline(Cubic(Line(OnGrid())))),rangexInt)
    intHGp1 = Interpolations.scale(interpolate(tabHGp1Int, BSpline(Cubic(Line(OnGrid())))),rangexInt)
    intHG2 = Interpolations.scale(interpolate(tabHG2Int, BSpline(Cubic(Line(OnGrid())))),rangexInt)
    intHGp2 = Interpolations.scale(interpolate(tabHGp2Int, BSpline(Cubic(Line(OnGrid())))),rangexInt)
    #####
    return [intHG1, intHGp1, intHG2, intHGp2]
end

"""
    tabHyperGeoInt

Interpolate the following hypergeometric function :
- `HG1   = _₂F₁(qCalc/2,qCalc-3.5,1.0,x)`
- `HGp1  = _₂F₁(qCalc/2+1,qCalc-2.5,2.0,x)`
- `HG2   = _₂F₁(qCalc/2,qCalc/2,4.5-qCalc/2,1/x)`
- `HGp2  = _₂F₁(qCalc/2+1,qCalc/2+1,5.5-qCalc/2,1/x)`

Interpolations are put inside a StaticArrays of size 4.
"""
const tabHyperGeoInt = SVector{4}(getHGInt())

"""
    H_1

Interpolate the following hypergeometric function :
- `HG1   = _₂F₁(qCalc/2,qCalc-3.5,1.0,x)`
"""
const H_1 = tabHyperGeoInt[1]

"""
    H_2

Interpolate the following hypergeometric function :
- `HGp1  = _₂F₁(qCalc/2+1,qCalc-2.5,2.0,x)`
"""
const H_2 = tabHyperGeoInt[2]

"""
    H_3

Interpolate the following hypergeometric function :
- `HG2   = _₂F₁(qCalc/2,qCalc/2,4.5-qCalc/2,1/x)`
"""
const H_3 = tabHyperGeoInt[3]

"""
    H_4

Interpolate the following hypergeometric function :
- `HGp2  = _₂F₁(qCalc/2+1,qCalc/2+1,5.5-qCalc/2,1/x)`
"""
const H_4 = tabHyperGeoInt[4]

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
@testset "1297.maximum-number-of-occurrences-of-a-substring.jl" begin
    @test furthest_building([4, 2, 7, 6, 9, 14, 12], 5, 1) == 4
    @test furthest_building([4, 12, 2, 7, 3, 18, 20, 3, 19], 10, 2) == 7
    @test furthest_building([14, 3, 19, 3], 17, 0) == 3
end


text:
# import Pkg; Pkg.activate(@__DIR__);
# Pkg.instantiate()
using RigidBodyDynamics
using StaticArrays
using RobotDynamics
using BenchmarkTools

# Defining Arthur model using RigidBodyDynamics
struct Arthur{C} <: AbstractModel
    mechanism::Mechanism{Float64}
    statecache::C
    dyncache::DynamicsResultCache{Float64}
    function Arthur(mechanism::Mechanism)
        statecache = StateCache(mechanism)
        rescache = DynamicsResultCache(mechanism)
        new{typeof(statecache)}(mechanism, statecache, rescache)
    end
end

#TODO: Change Path
Arthur(; mechanism=RigidBodyDynamics.URDF.parse_urdf(joinpath(@__DIR__,"../../../kortex_description/arms/gen3/7dof/urdf/Gen3_URDF_v13.urdf"), remove_fixed_tree_joints = false)) = Arthur(mechanism)

# State, s, is [q q̇ x ẋ F]
# x will be input from the camera
# q, q̇, ẋ will be taken or derived from the arm
# F will be input from the F/T Sensor
# Input, u, is Torque (τ)
function RobotDynamics.dynamics(model::Arthur, x::AbstractVector{T1}, u::AbstractVector{T2}) where {T1,T2}
    # Create a state of the mechanism model and a result struct for the dynamics
    T = promote_type(T1,T2)
    state = model.statecache[T]
    res = model.dyncache[T]

    # Get states and constants of system not dependent on model state
    # num_q = RigidBodyDynamics.num_positions(model.mechanism)
    # q = x[SA[1, 2, 3, 4, 5, 6, 7]]
    # q̇ = x[SA[8, 9, 10, 11, 12, 13, 14]]
    # p = x[2*num_q + 1:2*num_q + 6]
    # ṗ = x[2*num_q + 7:2*num_q + 12]
    # F = x[2*num_q + 13:2*num_q + 18]
    # Be = zeros(T, 6, 6)

    # if (norm(ṗ) > 1e-5)
    #     for k = 1:6
    #         Be[k,k] = norm(F) / norm(ṗ)
    #     end
    # end

    # Set mechanism state to current state
    # copyto!(parent(state.q), 1, x, 1, num_positions(state))
    # copyto!(parent(state.v), 1, x, num_positions(state) + 1, num_velocities(state))
    # setdirty!(state)
    copyto!(state, x[SA[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]])

    # w = Wrench{T}(default_frame(bodies(model.mechanism)[end-1]), F[4:6], F[1:3])
    # wrenches = BodyDict{Wrench{T}}(b => zero(Wrench{T}, root_frame(model.mechanism)) for b in bodies(model.mechanism))
    # wrenches[bodies(model.mechanism)[end-1].id] = transform(w, transform_to_root(mechanismState, bodies(model.mechanism)[end-1]))    dynamics!(dynamicsResult, mechanismState, u, wrenches)
    # dynamics!(res, state, u, wrenches)
    dynamics!(res, state, u)

    # p̈ = [res.accelerations[bodies(model.mechanism)[end].id].linear; res.accelerations[bodies(model.mechanism)[end].id].angular]
    # Ḟ = Be*p̈
    # return SVector{32}([q̇; q̈; ṗ; p̈; Ḟ])
    # return SVector{14}([x[SA[8, 9, 10, 11, 12, 13, 14]]; res.v̇])
    return [x[SA[8, 9, 10, 11, 12, 13, 14]]; res.v̇]
end

RobotDynamics.state_dim(::Arthur) = 14
RobotDynamics.control_dim(::Arthur) = 7

# cd(joinpath(@__DIR__,"../../../kortex_description/arms/gen3/7dof/urdf/GEN3_URDF_V12 copy.urdf"))
# readdir(joinpath(@__DIR__,"../../../kortex_description/arms/gen3/7dof/urdf/"))
# Arthur()

"""
Integrates the dynamics ODE 1 dt forward, x_{k+1} = rk4(x_k,u_k,dt).

returns x_{k+1}
"""
function rk4(model::Arthur,x::AbstractVector,u::AbstractVector,dt::Float64)
    # rk4 for integration
    k1 = dt*dynamics(model,x,u)
    k2 = dt*dynamics(model,x + k1/2,u)
    k3 = dt*dynamics(model,x + k2/2,u)
    k4 = dt*dynamics(model,x + k3,u)
    return x + (1/6)*(k1 + 2*k2 + 2*k3 + k4)
end

# model = Arthur()
# n,m = size(model)
# x = rand(n)
# u = rand(m)
# T1 = Float64
# T2 = Float64
# T = promote_type(T1,T2)
# state = model.statecache[T]
# res = model.dyncache[T]
# # copyto!(state, x)
# copyto!(parent(state.q), 1, x, 1, num_positions(state))
# copyto!(parent(state.v), 1, x, num_positions(state) + 1, num_velocities(state))
# @btime begin
#     setdirty!(state)
#     dynamics_bias!($res, $state)
#     mass_matrix!($res, $state)
#     # [x[SA[8, 9, 10, 11, 12, 13, 14]]; res.v̇]
# end
# @btime dynamics($model, $x, $u);
# suite = BenchmarkGroup()
# suite["test"] = @benchmarkable(begin

# end,
# # setup = begin

# # end,
# evals=1)
# overhead = BenchmarkTools.estimate_overhead()
# results = run(suite, verbose=true, overhead=overhead, gctrial=false)
# for result in results
#     println("$(first(result)):")
#     display(last(result))
#     println()
# end
# """
#     simulate(model, x0, ctrl; [kwargs...])

# Simulate `model` starting from `x0` using the `get_control(ctrl, x, t)` method to get the
# closed-loop feedback command.

# # Keyword Arguments
# * `tf`: final time
# * `dt`: simulation time step
# """
# function simulate(model::Arthur, x0, ctrl; tf=ctrl.times[end], dt=1e-2)
#     n,m = size(model)
#     times = range(0, tf, step=dt)
#     N = length(times)
#     X = [@SVector zeros(n) for k = 1:N]
#     U = [@SVector zeros(m) for k = 1:N-1]
#     X[1] = x0

#     tstart = time_ns()
#     for k = 1:N-1
#         U[k] = get_control(ctrl, X[k], times[k])
# #         u = clamp(U[k], umin, umax)
#         X[k+1] = discrete_dynamics(RK4, model, X[k], U[k], times[k], dt)
#     end
#     tend = time_ns()
#     rate = N / (tend - tstart) * 1e9
#     println("Controller ran at $rate Hz")
#     return X,U,times
# end


text:
# Use PCA to split points. SVD implementation.

import Random
import PyPlot
import Statistics

using LinearAlgebra
#import Interpolations
#using Revise

include("../src/RKHSRegularization.jl")
import .RKHSRegularization # https://github.com/RoyCCWang/RKHSRegularization

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])
import Random
import PyPlot
import Statistics

using LinearAlgebra
#import Interpolations
#using Revise

include("../src/RKHSRegularization.jl")
import .RKHSRegularization # https://github.com/RoyCCWang/RKHSRegularization

PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

PyPlot.close("all")

Random.seed!(25)

fig_num = 1

D = 2
N = 200
#X = randn(D,N)

import Distributions

μ = randn(D)
Σ = randn(D,D)
Σ = Σ'*Σ
dist = Distributions.MvNormal(μ, Σ)
X = rand(dist, N)

# ceneterd data matrix: zero mean and transposed.
Z = collect( X[:,n] - μ for n = 1:size(X,2) )
Z_mat = RKHSRegularization.array2matrix(Z)
Z_mat = Z_mat'

U, s, V = svd(Z_mat)

u = V[:,1]

X_set = collect( X[:,n] for n = 1:size(X,2) )
flags, functional_evals, c = RKHSRegularization.splitpoints(u, X_set)
X1 = X[:, flags]
X2 = X[:, .!flags]

# dot(u,x)+c = 0 to y = m*x + b.
function get2Dline(u::Vector{T}, c::T) where T

    m = -u[1]/u[2]
    b = c/u[2]

    return m, b
end

m, b = get2Dline(u, c)
t = LinRange(-5, 5, 300)
y = m .* t .+ b

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.scatter(X1[1,:], X1[2,:], label = "X1")
PyPlot.scatter(X2[1,:], X2[2,:], label = "X2")
#PyPlot.plot(y, t)
PyPlot.plot(t,y)

title_string = "split points"
PyPlot.title(title_string)
PyPlot.legend()
PyPlot.axis("equal")

max_stars_repo_name: RoyCCWang/RKHSRegularization

max_issues_repo_name: RoyCCWang/RKHSRegularization

max_forks_repo_name: RoyCCWang/RKHSRegularization


text:
using AstroLib: kepler_solver

struct KeplerianOrbit <: AbstractOrbit
    a
    ecc
    period
    ρ_star
    r_star
    m_star
    m_planet
    m_total
    n
    a_star
    a_planet
    K0
    M0
    t_periastron
    t0
    tref
    incl
    sin_incl
    cos_incl
end

const G_grav = 2942.2062175044193

function KeplerianOrbit(; period, m_planet=0, r_star=1, m_star=1, ρ_star = 3 * m_star / (4 * π * r_star^3), ecc=nothing, t_periastron=0, t0 = 0 + t_periastron)
    a = cbrt(G_grav * (m_star + m_planet) * period^2 / (4 * π^2))
    m_total = m_star + m_planet
    n = 2 * π / period
    a_star = a * m_planet / m_total
    a_planet = -a * m_star / m_total
    K0 = n * a / m_total
    if ecc === nothing
        M0 = 0.5 * π
        incl_factor = 1
    else
        error("not implemented")
    end
    incl = 0.5 * π
    sin_incl, cos_incl = sincos(incl)

    tref = t_periastron - t0
    return KeplerianOrbit(a,
                          ecc,
                          period,
                          ρ_star,
                          r_star,
                          m_star,
                          m_planet,
                          m_total,
                          n,
                          a_star,
                          a_planet,
                          K0,
                          M0,
                          t_periastron,
                          t0,
                          tref,
                          incl,
                          sin_incl,
                          cos_incl)
end

function relative_position(orbit::KeplerianOrbit, t)
    sinf, cosf = get_true_anomaly(orbit, t)
    if orbit.ecc === nothing
        r = -orbit.a
    else
        r = -orbit.a * (1 - orbit.ecc^2) / (1 + orbit.ecc * cosf)
    end
    return rotate_vector(orbit, r * cosf, r * sinf)
end

function get_true_anomaly(orbit::KeplerianOrbit, t)
    M = orbit.n * ((t - orbit.t0) - orbit.tref)
    return kepler(M, orbit.ecc)
end

function kepler(M, ecc)
    E = kepler_solver(M, ecc)
    return sincos(E)
end
kepler(M, ::Nothing) = sincos(M)

function rotate_vector(orbit::KeplerianOrbit, x, y)
    # rotate about z0 axis by ω
    if orbit.ecc === nothing
        x1, y1 = x, y
    else
        # x1 = orbit.cosω * x - orbit.sinω * y
        # y1 = orbit.sinω * x + orbit.cosω * y
    end

    # rotate about x1 axis by -incl
    x2 = x1
    y2 = orbit.cos_incl * y1
    Z = -orbit.sin_incl * y1

    return SA[x2, y2, Z]

    # rotate about z2 axis by Ω
    # if orbit.Ω === nothing
    #     return SA[x2, y2, Z]
    # end
    # X = orbit.cosΩ * x2 - orbit.sinΩ * y2
    # Y = orbit.sinΩ * x2 + orbit.cosΩ * y2
    # return SA[X, Y, Z]
end


text:
# -- Dimensions ---------------------------------------------------------------

type PartialDimension
    name::Symbol
end

typealias Dimension Union(PartialDimension, Int)

#show(d::PartialDimension) = print("dimension(:$(d.name))")
show(d::PartialDimension) = print(d.name)

function broadcast_dim(n, m)
    if (n == m) || (n == 1)
        return m
    elseif m == 1
        return n
    end
    error("Cannot broadcast together dimensions: $n != $m")
end

function broadcast_shape(s::(Dimension...), t::(Dimension...))
    (ns, nt) = (length(s), length(t))
    if ns < nt
        return broadcast_shape(t, s)
    elseif ns > nt
        return tuple(broadcast_shape(s[1:nt], t)..., s[(nt+1):ns]...)
    end
    return map(broadcast_dim, s, t)
end

## Factory ##
dimension(x::Int) = x
dimension(name::Symbol) = PartialDimension(name)

# == PartialArray =============================================================

# Array depending on yet unknown values
abstract PartialArray{T,N} <: AbstractArray{T,N}

# required fields: shape

size(A::PartialArray) = A.shape
size{T,N}(A::PartialArray{T,N}, d) = (d <= N ? A.shape[d] : 1)

show{T,N}(A::PartialArray{T,N}) = print(A.name)

peel(args::Tuple) = map(peel, args)
peel(args...) = peel(args)

#uinds = gensym(32) # just some number of dims that should be enough
uinds = (:_i, :_j, :_k, :_l) # todo: use gensyms instead

# -- InputArray ---------------------------------------------------------------

type InputArray{T,N} <: PartialArray{T,N}
    name::Symbol
    shape::(Dimension...)

    #todo: check length(shape) == N!
end

## Factory ##
input_array(name, T, shape) = InputArray{T,length(shape)}(name, shape)
input_array(name, T, shape...) = input_array(name, T, shape)

peel(A::InputArray) = expr(:ref, A.name, uinds[1:ndims(A)]...)

# -- OuputArray ---------------------------------------------------------------

#type OutputArray{T,N} <: PartialArray{T,N}
type OutputArray{T,N}  <: AbstractArray{T,N} # don't want outputs as inputs
    name::Symbol
    shape::(Dimension...)

    writes::Vector

    OutputArray(name, shape) = new(name, shape, {})
end

## Factory ##
output_array(name, T, shape) = OutputArray{T,length(shape)}(name, shape)
output_array(name, T, shape...) = output_array(name, T, shape)

function assign{T,N}(dest::OutputArray{T,N}, A::PartialArray{T,N})
    ref = expr(:ref, dest.name, uinds[1:ndims(dest)]...)
    a = peel(A)
    write_op = :( ($ref) = $a )
    append!(dest.writes, {write_op})
end

# -- KernelArray --------------------------------------------------------------

# partial array composed of partial elements
type KernelArray{T,N} <: PartialArray{T,N}
    # inds?
    ex::Expr

    shape::(Dimension...)
end

KernelArray{T}(ex, ::Type{T}, shape) = KernelArray{T,length(shape)}(ex, shape)

peel(A::KernelArray) = A.ex

# -- PartialArray math --------------------------------------------------------

function .*{TA,TB}(A::PartialArray{TA}, B::PartialArray{TB})
    shape = broadcast_shape(size(A), size(B))
    T = promote_type(TA, TB)

    (a, b) = peel(A, B)
    ex = :($a*$b)  # ($a).*($b) instead?
    KernelArray(ex, T, shape)
end

# -- Test code ----------------------------------------------------------------

if false
    VType = Float32

    n = dimension(:n)
    m = dimension(:m)

    A = input_array(:A, VType, (n,m))
    B = input_array(:B, VType, (n,m))
    dest = output_array(:dest, VType, (n,m))

    C = A.*B
    dest[] = C
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
include("roadways.jl")
include("vehicles.jl")
include("overlays.jl")


text:
## This file contains direct Julia implementations of 2d and 3d convolutions
using Base.Threads

# Helper functions for restricting x/w overreach
function clamp_lo(x, w)
    idx = 1
    while idx <= length(x) && x[idx] <= 0
        idx += 1
    end
    return (x[idx:end], w[idx:end])
end
function clamp_hi(x, w, L)
    idx = length(x)
    while idx >= 1 && x[idx] > L
        idx -= 1
    end
    return (x[1:idx], w[1:idx])
end

"""
    conv_direct!(y, x, w, cdims; alpha=1, beta=0)

Direct convolution implementation; used for debugging, tests, and mixing/matching of
strange datatypes within a single convolution.  Uses naive nested for loop implementation
and does not attempt to optimize performance.  Rather, this implementation is intended to
be maximally understandable and debuggable, to aid in testing other, more performant
implementations.  We also explicitly support mixing and matching of strange datatypes,
so that if the user really wants to convolve an image of `UInt8`'s with a `Float16`
kernel, storing the result in a `Float32` output, there is at least a function call
for that madness.

The keyword arguments `alpha` and `beta` control accumulation behavior; this function
calculates `y = alpha * x * w + beta * y`, therefore by setting `beta` to a nonzero
value, the user is able to accumulate values into a preallocated `y` buffer, or by
setting `alpha` to a nonunitary value, an arbitrary gain factor can be applied.

By defaulting `beta` to `false`, we make use of the Bradbury promotion trick to override
`NaN`'s that may pre-exist within our output buffer, as `false*NaN == 0.0`, whereas
`0.0*NaN == NaN`.  Only set `beta` if you are certain that none of the elements within
`y` are `NaN`.

The basic implementation performs 3-dimensional convolution; 1-dimensional and 2-
dimensional casesa are supported by simply reshaping `y`, `x` and `w`, for which
wrapper methods are available.
"""
conv_direct!

function conv_direct!(y::AbstractArray{yT,5}, x::AbstractArray{xT,5},
                      w::AbstractArray{wT,5}, cdims::DenseConvDims;
                      alpha::yT = yT(1), beta = false) where {yT, xT, wT}
    check_dims(size(x), size(w), size(y), cdims)

    width, height, depth = input_size(cdims)
    kernel_w, kernel_h, kernel_d = kernel_size(cdims)
    out_c = channels_out(cdims)
    pad_w_lo, pad_w_hi, pad_h_lo, pad_h_hi, pad_d_lo, pad_d_hi = padding(cdims)
    dil_w, dil_h, dil_d = dilation(cdims)
    stride_w, stride_h, stride_d = stride(cdims)
    out_width, out_height, out_depth = output_size(cdims)

    # Create a method that, at compile-time, determines how we're going to index into `w`
    kproj(k, M, cdims::ConvDims{N,S,P,D,true}) where {N, S, P, D} = k
    kproj(k, M, cdims::ConvDims{N,S,P,D,false}) where {N, S, P, D} = M - k + 1

    # A helper function to project from output (w, h) to input (input_w, input_h)
    project(idx, stride, pad) = (idx - 1)*stride - pad + 1

    # Use `calc_padding_regions` to determine where we do or don't need to worry about padding
    padded_regions, central_region = calc_padding_regions(cdims)

    # Start with the central region
    w_region, h_region, d_region = central_region
    @inbounds for batch in 1:size(x, 5),
        c_out in 1:out_c,
        d_idx in d_region,
        h_idx in h_region,
        w_idx in w_region

        # Since we're in the central region, we don't need to worry about clamping
        dotprod = yT(0)
        for c_in in 1:channels_in(cdims),
            kd in 1:kernel_d,
            kh in 1:kernel_h,
            kw in 1:kernel_w

            # Hoist me, you coward.
            x_d = project(d_idx, stride_d, pad_d_lo) + (kd - 1)*dil_d
            x_h = project(h_idx, stride_h, pad_h_lo) + (kh - 1)*dil_h
            x_w = project(w_idx, stride_w, pad_w_lo) + (kw - 1)*dil_w

            x_val = x[x_w, x_h, x_d, c_in, batch]
            w_val = w[kproj(kw, kernel_w, cdims),
                    kproj(kh, kernel_h, cdims),
                    kproj(kd, kernel_d, cdims),
                    c_in, c_out]
            dotprod = muladd(x_val, w_val, dotprod)
        end
        y[w_idx, h_idx, d_idx, c_out, batch] = alpha*dotprod + beta*y[w_idx, h_idx, d_idx, c_out, batch]
    end

    # Next, do potentially-padded regions:
    @inbounds for (w_region, h_region, d_region) in padded_regions,
        batch in 1:size(x, 5),
        c_out in 1:out_c,
        d_idx in d_region,
        h_idx in h_region,
        w_idx in w_region

        # Probe for out-of-bounds accesses on `x` and `continue` if we hit one
        dotprod = yT(0)
        for c_in in 1:channels_in(cdims),
            kd in 1:kernel_d

            x_d = project(d_idx, stride_d, pad_d_lo) + (kd - 1)*dil_d
            if x_d <= 0 || x_d > depth
                continue
            end

            for kh in 1:kernel_h
                x_h = project(h_idx, stride_h, pad_h_lo) + (kh - 1)*dil_h
                if x_h <= 0 || x_h > height
                    continue
                end

                for kw in 1:kernel_w
                    x_w = project(w_idx, stride_w, pad_w_lo) + (kw - 1)*dil_w
                    if x_w <= 0 || x_w > width
                        continue
                    end

                    x_val = x[x_w, x_h, x_d, c_in, batch]
                    w_val = w[kproj(kw, kernel_w, cdims),
                            kproj(kh, kernel_h, cdims),
                            kproj(kd, kernel_d, cdims),
                            c_in, c_out]
                    dotprod = muladd(x_val, w_val, dotprod)
                end
            end
        end

        y[w_idx, h_idx, d_idx, c_out, batch] = alpha*dotprod + beta*y[w_idx, h_idx, d_idx, c_out, batch]
    end

    return y
end

## Gradient definitions
"""
    ∇conv_data_direct!(dx, dy, w, cdims; alpha=1, beta=0)

Calculate the gradient imposed upon `x` in the convolution `y = x * w`.
"""
∇conv_data_direct!

function ∇conv_data_direct!(dx::AbstractArray{xT,5}, dy::AbstractArray{yT,5},
                            w::AbstractArray{wT,5}, cdims::DenseConvDims;
                            alpha::xT=xT(1), beta=false) where {xT, yT, wT}
    w = transpose_swapbatch(w[end:-1:1, end:-1:1, end:-1:1, :, :])
    dy = predilate(dy, stride(cdims))
    ctdims = DenseConvDims(dy, w; padding=transpose_pad(cdims),
                                  dilation=dilation(cdims),
                                  flipkernel=flipkernel(cdims))
    dx = conv_direct!(dx, dy, w, ctdims; alpha=alpha, beta=beta)
    return transpose_swapbatch(dx)
end

"""
    ∇conv_filter_direct!(dw, x, dy, cdims; alpha=1, beta=0)

Calculate the gradient imposed upon `w` in the convolution `y = x * w`.
"""
∇conv_filter_direct!

function ∇conv_filter_direct!(dw::AbstractArray{wT,5}, x::AbstractArray{xT,5},
                              dy::AbstractArray{yT,5}, cdims::DenseConvDims;
                              alpha::wT=wT(1), beta=false) where {xT, yT, wT}
    x = transpose_swapbatch(x[end:-1:1, end:-1:1, end:-1:1, :, :])
    dy = transpose_swapbatch(predilate(dy, stride(cdims)))
    ctdims = DenseConvDims(dy, x; padding=transpose_pad(cdims),
                                    stride=dilation(cdims))
    conv_direct!(dw, dy, x, ctdims; alpha=alpha, beta=beta)
    if flipkernel(cdims)
        dw .= dw[end:-1:1, end:-1:1, end:-1:1, :, :]
    end
    return dw
end


text:
using Base.Meta

"""
    AbstractVertex{P}

A supertype for all vertex types with a payload of type `P`.
"""
abstract type AbstractVertex{P} end

"""
    payload(v)

Get a vertex's payload
"""
payload(v::AbstractVertex) = v.payload

zero!(v::AbstractVertex{P}) where P = v.payload = zero(P)

"""
    id(v)

Get the id of a vertex.
"""
id(v::AbstractVertex) = v.id

"""
    above(v)

Get an array of all vertices above `v`.
"""
above(v::AbstractVertex) = v.above

"""
    below(v)

Get an array of all vertices below `v`.
"""
below(v::AbstractVertex) = v.below

Base.length(v::AbstractVertex) = length(id(v))

Base.eachindex(v::AbstractVertex) = eachindex(id(v))

Base.getindex(v::AbstractVertex, idx...) = getindex(id(v), idx...)

bookend(x::AbstractString, left::AbstractString="[", right::AbstractString="]") = left * x * right

prettyname(v::AbstractVertex) = bookend(join(repr.(id(v)), ", "))

function Base.show(io::IO, v::AbstractVertex)
    ex = Meta.parse(string(typeof(v)))
    t = if ex isa Symbol
        string(ex)
    else
        string(ex.args[1])
    end
    print(io, t, "(", prettyname(v), ", ", payload(v), ")")
end

function Base.show(io::IO, ::MIME"text/dot", v::AbstractVertex)
    print(io, "[label=\"", prettyname(v), "\\n")
    try
        show(io, MIME("text/dot"), payload(v))
    catch
        print(io, payload(v))
    end
    print(io, "\"]")
end

abstract type AbstractUnnamedVertex{P} <: AbstractVertex{P} end

abstract type AbstractNamedVertex{N,P} <: AbstractVertex{P} end

name(v::AbstractNamedVertex) = v.name

function prettyname(v::AbstractNamedVertex)
    ns = String[]
    for n in name(v)
        push!(ns, if length(n) == 1
            join(repr.(n), ", ")
        else
            bookend(join(repr.(n), ", "), "{", "}")
        end)
    end
    bookend(join(ns, ", "))
end

"""
    VertexID = Vector{Int}

A type alias for integer arrays used as ids of vertices.
"""
const VertexID = Vector{Int}

"""
    UnnamedVertex{P} <: AbstractVertex{P}

A standard vertex containing a `id`, `payload`, and arrays of vertices
`above` and `below` it.
"""
mutable struct UnnamedVertex{P, V <: AbstractUnnamedVertex{P}} <: AbstractUnnamedVertex{P}
    id::VertexID
    payload::P
    above::Vector{V}
    below::Vector{V}
end

"""
    UnnamedVertex(id, p)

Construct a `Vertex` with a given `id` and payload (`p`), with no vertices
above or below.
"""
function UnnamedVertex(id::VertexID, p::P) where P
    UnnamedVertex(id, p, UnnamedVertex{P}[], UnnamedVertex{P}[])
end

"""
    UnnamedVertex{P}(id)

Construct a `Vertex` with a given `id`, a "zeroed" payload of type `P` and no
vertices above or below.
"""
UnnamedVertex{P}(id::VertexID) where P = UnnamedVertex(id, zero(P))

clone(v::UnnamedVertex) = UnnamedVertex(id(v), payload(v))

mutable struct Vertex{N, P, V <: AbstractNamedVertex{N,P}} <: AbstractNamedVertex{N,P}
    id::VertexID
    name::Vector{Vector{N}}
    payload::P
    above::Vector{V}
    below::Vector{V}
end

function Vertex(id::VertexID, name::AbstractVector{Vector{N}}, p::P) where {N,P}
    Vertex(id, name, p, Vertex{N,P}[], Vertex{N,P}[])
end

function Vertex{N,P}(id::VertexID, name::AbstractVector{Vector{N}}) where {N,P}
    Vertex(id, name, zero(P))
end

clone(v::Vertex) = Vertex(id(v), name(v), payload(v))

function Base.convert(::Type{UnnamedVertex{T}}, v::Vertex{N,T}) where {N,T}
    a = convert(Vector{UnnamedVertex{T}}, above(v))
    b = convert(Vector{UnnamedVertex{T}}, below(v))
    UnnamedVertex(id(v), payload(v), a, b)
end

"""
    isbelow(a, b)

Determine whether or not `a` is below `b`.
"""
@inline function isbelow(xs::VertexID, ys::VertexID)
    for i in eachindex(ys)
        is_valid = false
        for j in eachindex(xs)
            if xs[j] & ys[i] == xs[j]
                is_valid = true
                break
            end
        end
        if !is_valid
            return false
        end
    end
    true
end

isbelow(a::AbstractVertex, b::AbstractVertex) = isbelow(id(a), id(b))

"""
    genvertices(::Type{V}, n) where {V <: AbstractUnnamedVertex}
"""
function genvertices(::Type{V}, n::Int) where {V <: AbstractUnnamedVertex}
    if n < 1
        throw(DomainError(n, "at least one node is required"))
    end
    vs = V[]
    m = (1 << n) - 1;
    for i in 1:m
        push!(vs, V([i]))
        genvertices!(vs, i + 1, m, [i])
    end
    vs
end

function genvertices!(vs::AbstractVector{V}, i, m, c) where {V <: AbstractUnnamedVertex}
    if i < m
        genvertices!(vs, i+1, m, c[:])
    end

    if i <= m
        z = 0
        for j in 1:length(c)
            z = i & c[j]
            if (z == i || z == c[j])
                return vs
            end
        end
        push!(c, i)
        push!(vs, V(c))
        genvertices!(vs, i+1, m, c[:])
    end

    vs
end

function extractname(nodenames::AbstractVector{N}, id::Int) where N
    name = N[]
    for i in eachindex(nodenames)
        if id & (1 << (i-1)) != 0
            push!(name, nodenames[i])
        end
    end
    name
end

extractname(nodenames::AbstractVector, id::VertexID) = map(i -> extractname(nodenames, i), id)

function genvertices(::Type{V}, names::AbstractVector{N}) where {N, V <: AbstractNamedVertex{N}}
    if isempty(names)
        throw(DomainError(n, "at least one node is required"))
    end
    n = length(names)
    vs = V[]
    m = (1 << n) - 1;
    for i in 1:m
        id = [i]
        push!(vs, V(id, extractname(names, id)))
        genvertices!(vs, names, i + 1, m, [i])
    end
    vs
end

function genvertices!(vs::AbstractVector{V}, names, i, m, c) where {V <: AbstractNamedVertex}
    if i < m
        genvertices!(vs, names, i+1, m, c[:])
    end

    if i <= m
        z = 0
        for j in 1:length(c)
            z = i & c[j]
            if (z == i || z == c[j])
                return vs
            end
        end
        push!(c, i)
        push!(vs, V(c, extractname(names, c)))
        genvertices!(vs, names, i+1, m, c[:])
    end

    vs
end

"""
    toposort!(vs)

Topologically sort vertices according to the `isbelow` relation.
"""
function toposort!(vs::AbstractVector{V}) where {V <: AbstractVertex}
    n, u, v = length(vs), 1, 1
    while v < n - 1
        u = v
        for i in u:n
            is_bottom = true
            for j in u:n
                if i != j && isbelow(vs[j], vs[i])
                    is_bottom = false
                    break
                end
            end
            if is_bottom
                vs[v], vs[i] = vs[i], vs[v]
                v += 1
            end
        end
    end
    vs
end

"""
    Hasse{V}

A Hasse diagram (lattice) with vertices of type `V`
"""
mutable struct Hasse{V <: AbstractVertex}
    top::V
    bottom::V
    vertices::Vector{V}
end

function Hasse(::Type{P}, n::Int) where P
    if n > 5
        @warn """Hasse diagrams are computationally expensive to build for n > 5.
        This is gonna take a while and lot of memory""" n
    end
    Hasse(genvertices(UnnamedVertex{P}, n))
end

function Hasse(::Type{P}, names::AbstractVector{N}) where {N,P}
    if length(names) > 5
        @warn """Hasse diagrams are computationally expensive to build for n > 5.
        This is gonna take a while and lot of memory""" n=length(names)
    end
    Hasse(genvertices(Vertex{N,P}, names))
end

function Hasse(vs::AbstractVector{V}; sort::Bool=true) where {V <: AbstractVertex}
    if isempty(vs)
        throw(ArgumentError("vertex list is empty"))
    end

    if sort
        toposort!(vs)
    end

    for i in 1:length(vs)
        for j in i+1:length(vs)
            if isbelow(vs[i], vs[j])
                stop = false
                for v in above(vs[i])
                    if isbelow(v, vs[j])
                        stop = true
                        break
                    end
                end
                if stop
                    break
                else
                    push!(above(vs[i]), vs[j])
                    push!(below(vs[j]), vs[i])
                end
            end
        end
    end

    Hasse(vs[end], vs[1], vs)
end

"""
    top(h)

Get the top vertex of the Hasse diagram `h`.
"""
top(h::Hasse) = h.top

"""
    bottom(h)

Get the bottom vertex of the Hasse diagram `h`.
"""
bottom(h::Hasse) = h.bottom

"""
    vertices(h)

Get the array of vertices of the Hasse diagram `h`.
"""
vertices(h::Hasse) = h.vertices

Base.length(h::Hasse) = length(vertices(h))

Base.eachindex(h::Hasse) = eachindex(vertices(h))

Base.getindex(h::Hasse, idx...) = getindex(vertices(h), idx...)

zero!(h::Hasse) = zero!.(vertices(h))

"""
    edgelist(h)

Get the array of edges — pairs of vertex ids — of the Hasse diagram `h`.
"""
function edgelist(h::Hasse)
    edges = NTuple{2, VertexID}[]
    vs = vertices(h)
    for v in vs, w in v.above
        push!(edges, (id(v), id(w)))
    end
    edges
end

function prune(f::Function, h::Hasse)
    vs = clone.(filter(f, vertices(h)))
    if isempty(vs)
        throw(ErrorException("pruned diagram has no vertices"))
    end
    Hasse(vs)
end
prune(h::Hasse) = prune(!iszero ∘ payload, h)

Base.istextmime(::MIME"text/dot") = true

function Base.show(io::IO, ::MIME"text/dot", h::Hasse; names=nothing)
    println(io, "digraph Hasse {")
    println(io, "  rankdir=\"BT\";")
    vs = vertices(h)
    for i in eachindex(vs)
        print(io, "  ", i, " ")
        show(io, MIME("text/dot"), vs[i])
        println(io, ";")
    end
    println(io)
    for i in eachindex(vs)
        for v in above(vs[i])
            j = findfirst(w -> v == w, vs)
            println(io, "  ", i, " -> ", j, ";")
        end
    end
    println(io, "}")
end

function graphviz(filename::AbstractString, h::Hasse;
                  format=lowercase(last(splitext(filename)))[2:end])
    open(`dot -T$format -o$filename`, "w", stdout) do io
        show(io, MIME("text/dot"), h)
    end
end


text:
# edges = [[2000, [35, 40], [355,500,600], 1000],
#         [2001, [35, 50], [100,500,600], 2000],
#         [2002, [37, 40], [100,400,600], 3000]];
# edges = hcat(edges...)';

# create_graph_dicts(edges)

# # pre_to_post, post_to_pre = create_post_pre_dicts(segs)
# @test length(pre_to_post[35]) == 2
# @test pre_to_post[35][1] == 40
# @test pre_to_post[35][2] == 50
# @test pre_to_post[37][1] == 40
# @test post_to_pre[40][1] == 35
# @test post_to_pre[40][2] == 37
# @test post_to_pre[50][1] == 35


text:
using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")
using JLD2
using NNFEM

tid = 1
@load "../Data/domain$tid.jld2" domain

nodes = domain.nodes

plot_line = (i, j)->begin
    x = [nodes[i,1]+u[i];nodes[j,1]+u[j]]
    y = [nodes[i,2]+v[i];nodes[j,2]+v[j]]
    plot(x, y, color="grey",alpha=0.5)
end

close("all")
scatter(nodes[:, 1], nodes[:,2], color="red")
u,v = domain.state[1:domain.nnodes], domain.state[domain.nnodes+1:end]
scatter(nodes[:, 1] + u, nodes[:,2] + v, color="blue")
for i = 1:3
    plot_line(2i, 2i+2)
    plot_line(2i, 2i-1)
    plot_line(2i, 2i+1)
    plot_line(2i-1, 2i+1)
    plot_line(2i-1, 2i+2)
    plot_line(2i+1, 2i+2)
end
ylim(-0.2,1.5)
xlabel("x")
ylabel("y")
# mpl.save("truss2d_loc.tex")
savefig("truss2d_loc.pdf")

close("all")
X, Y = prepare_strain_stress_data1D(domain)
scatter(X[:,1], Y, marker=".",s=5)
xlabel("Strain")
ylabel("Stress")
grid("on")
# mpl.save("truss2d_stress.tex")
savefig("truss2d_stress.pdf")


text:
const TLS_BUSY_READ_SECS = 1
const TLS_BUSY_READ_YIELD_SECS = 0.001
const TLS_READBUFF_SIZE = MbedTLS.MBEDTLS_SSL_MAX_CONTENT_LEN * 5
const TLS_MIN_WRITEBUFF_SIZE = MbedTLS.MBEDTLS_SSL_MAX_CONTENT_LEN
const TCP_MAX_WRITEBUFF_SIZE = 1024*512
const TCP_MIN_WRITEBUFF_SIZE = 1024*64

struct BufferedTLSSocket <: IO
    in::IOBuffer                # no read lock, single task reads socket and distributes messages to channels
    out::IOBuffer
    sock::MbedTLS.SSLContext
    readbuff::Vector{UInt8}
    out_lck::ReentrantLock      # protect out::IOBuffer when there are multiple channels on the connection

    function BufferedTLSSocket(sock::MbedTLS.SSLContext; readbuff_size::Int=TLS_READBUFF_SIZE)
        new(PipeBuffer(), PipeBuffer(), sock, Vector{UInt8}(undef, readbuff_size), ReentrantLock())
    end
end

isopen(bio::BufferedTLSSocket) = isopen(bio.sock)
close(bio::BufferedTLSSocket) = close(bio.sock)

function read(bio::BufferedTLSSocket, ::Type{UInt8})
    fill_in(bio, 1)
    read(bio.in, UInt8)
end

function read(bio::BufferedTLSSocket, T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}})
    fill_in(bio, sizeof(T))
    read(bio.in, T)
end

function read!(bio::BufferedTLSSocket, buff::Vector{UInt8})
    fill_in(bio, length(buff))
    read!(bio.in, buff)
end

function peek(bio::BufferedTLSSocket, T::Union{Type{Int16},Type{UInt16},Type{Int32},Type{UInt32},Type{Int64},Type{UInt64},Type{Int128},Type{UInt128},Type{Float16},Type{Float32},Type{Float64}})
    fill_in(bio, sizeof(T))
    peek(bio.in, T)
end

function fill_in(bio::BufferedTLSSocket, atleast::Int)
    avail = bytesavailable(bio.in)
    if atleast > avail
        while (atleast > avail) && isopen(bio.sock)
            bytes_read = isreadable(bio.sock) ? readbytes!(bio.sock, bio.readbuff; all=false) : 0
            if bytes_read > 0
                avail += Base.write(bio.in, first(bio.readbuff, bytes_read))
            else
                MbedTLS.wait_for_decrypted_data(bio.sock)
            end
        end
    end
end

function write(bio::BufferedTLSSocket, data::UInt8)
    lock(bio.out_lck) do
        write(bio.out, data)
    end
end
function write(bio::BufferedTLSSocket, data::Union{Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128,Float16,Float32,Float64})
    lock(bio.out_lck) do
        write(bio.out, data)
    end
end
function write(bio::BufferedTLSSocket, data::Array)
    lock(bio.out_lck) do
        write(bio.out, data)
    end
end
function flush(bio::BufferedTLSSocket)
    lock(bio.out_lck) do
        write(bio.sock, take!(bio.out))
    end
    nothing
end


text:
module testPQPolygonSDK

# include -
include("Include.jl")

# methods and types that we export -
export model
export url
export api

# types -
export AbstractPolygonEndpointModel
export PQPolygonSDKUserModel

# market endpoints -
export PolygonPreviousCloseEndpointModel
export PolygonAggregatesEndpointModel
export PolygonOptionsContractReferenceEndpoint
export PolygonGroupedDailyEndpointModel
export PolygonDailyOpenCloseEndpointModel

# reference endpoints -
export PolygonTickerNewsEndpointModel
export PolygonTickerDetailsEndpointModel
export PolygonMarketHolidaysEndpointModel #ycpan
export PolygonExchangesEndpointModel #ycpan
export PolygonStockSplitsEndpointModel  #ycpan
export PolygonMarketStatusEndpointModel #ycpan
export PolygonDividendsEndpointModel #ycpan
export PolygonTickersEndpointModel #ycpan
export PolygonConditionsEndpointModel #ycpan
export PolygonStockFinancialsEndpointModel #ycpan
export PolygonTickerTypesEndpointModel #ycpan

end # module


text:
"""
    struct Affine
        w
        b
        σ
    end

Affine layer of neural network

"""
struct Affine{T1<:AbstractArray,T2}
    w::T1
    b::T2
    σ::Function
end

function Affine(
    in::Integer,
    out::Integer,
    σ = identity::Function;
    fw = randn::Function,
    fb = zeros::Function,
    isBias = true::Bool,
    precision = Float32,
)
    if isBias
        return Affine(fw(precision, out, in), fb(precision, out), σ)
    else
        return Affine(fw(precision, out, in), Flux.Zeros(precision, out), σ)
    end
end

Flux.@functor Affine # works with Flux.params

(L::Affine)(x::AbstractVector) = L.w * x .+ L.b

"""
    struct FastAffine{I,F,F2} <: DiffEqFlux.FastLayer
        out::I
        in::I
        σ::F
        initial_params::F2
    end
Equivalent FastDense layer with controllable type
"""
struct FastAffine{I,F,F2} <: DiffEqFlux.FastLayer
    out::I
    in::I
    σ::F
    initial_params::F2
end

function FastAffine(
    in::Integer,
    out::Integer,
    σ = identity;
    fw = randn,
    fb = Flux.zeros,
    precision = Float32,
)
    initial_params() = vcat(vec(fw(precision, out, in)), fb(precision, out))
    return FastAffine(out, in, σ, initial_params)
end

(f::FastAffine)(x, p) =
    f.σ.(reshape(p[1:(f.out*f.in)], f.out, f.in) * x .+ p[(f.out*f.in+1):end])

DiffEqFlux.paramlength(f::FastAffine) = f.out * (f.in + 1)
DiffEqFlux.initial_params(f::FastAffine) = f.initial_params()


text:
struct RPMS{AT} <: IRMeasurement
    sig::AT
    prepad::Int
    gain::Float64
end

"""
    rpms(samples; options...)

Generate periodic noise using the random phase multisine medhod [1], by starting
with a flat spectrum and randomizing the phases. Note that this is,
nondeterministic, so make sure you use the same measurement instance to generate
and analyze your stimulus.

## Options
$optiondoc_gain
$optiondoc_prepad

[1]: Mateljan, I., Signal selection for the room acoustics measurement, 1999,
IEEE Workshop on Applications of Signal Processing to Audio and Acoustics,
pp. 199–202.
"""
function rpms(L; prepad=L, gain=rpms_gain)
    N = L÷2+1
    X = ones(ComplexF64, N)
    # set phases to be uniform random in [0,2π]
    @. X *= exp(im*2π*rand())
    # DC component should be zero
    X[1] = 0.0
    if iseven(L)
        # Nyquist component should be real (it needs to equal its complex
        # conjugate)
        X[N] = 1.0
    end
    x = irfft(X, L) * sqrt(L)
    RPMS(x, prepad, gain)
end

# TODO: add trials support
stimulus(s::RPMS) = [zeros(s.prepad); s.sig * s.gain; zeros(length(s.sig))]
prepadding(s::RPMS) = s.prepad

function _analyze(s::RPMS, response::AbstractArray)
    mapslices(response, dims=1) do v
        xcorr(v[s.prepad+1:end], s.sig, padmode=:longest)[end÷2+1:end] / length(s.sig) / s.gain
    end
end