text:
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v3+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v3+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v4+json")}) = true
Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v5+json")}) = true

function _maybe_fallback_to(mime, d, x, payload)
    if sizeof(payload) > d.config.max_json_bytes
        @warn string(
            "The size of JSON representation (", sizeof(payload), ")",
            " exceeds `max_json_bytes = ", d.config.max_json_bytes, "`.",
            " Falling back to $mime."
        )
        return display(d, mime, x)
    end
    return nothing
end

function _display_vegalite(d, major_version_vegalite, major_version_vega, x)
    payload = stringmime(MIME("application/vnd.vegalite.v$major_version_vegalite+json"), x)
    ans = _maybe_fallback_to(MIME"image/png"(), d, x, payload)
    ans === nothing || return ans

    html_page = """

    """

    displayhtml(d, html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

function _display_vega(d, major_version, x)
    payload = stringmime(MIME("application/vnd.vega.v$major_version+json"), x)
    ans = _maybe_fallback_to(MIME"image/png"(), d, x, payload)
    ans === nothing || return ans

    html_page = """

    """

    displayhtml(d, html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    return _display_vegalite(d, "2", "3", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v3+json")}, x)
    return _display_vegalite(d, "3", "5", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v3+json")}, x)
    return _display_vega(d, "3", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v4+json")}, x)
    return _display_vega(d, "4", x)
end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v5+json")}, x)
    return _display_vega(d, "5", x)
end


text:
function create_MPC(N)
    nx = 6
    nu = 3
    # N= 20
    xg = [-1;-2;-3;0;0;0]
    nz=(N*nx) + (N-1)*nu
    Ac = Array([zeros(3,3) I(3);zeros(3,6)])
    Bc = Array([zeros(3,3);I(3)])

    dt = 0.2
    H = exp(dt*[Ac Bc; zeros(3,9)])
    Ad = sparse(H[1:nx,1:nx])
    Bd = sparse(H[1:nx,nx+1:end])

    xi = [(nu+nx)*(i-1) .+ (1:nx)  for i = 1:N]

    ui = [(nu+nx)*(i-1) .+ (nx .+ (1:nu))  for i = 1:N-1]
    ci = [(nx)*(i-1) .+ (1:nx)  for i = 1:N]

    # cost function
    Qc = 10*I(nx)
    Rc = I(nu)
    Q = spzeros(nz,nz)
    q = zeros(nz)
    for i = 1:N
        if i<N
            Q[xi[i],xi[i]] = Qc
            q[xi[i]] = -Qc*xg
            Q[ui[i],ui[i]] = Rc
        else
            Q[xi[i],xi[i]] = Qc
            q[xi[i]] = -Qc*xg
        end
    end

    A = spzeros(N*nx,nz)
    for i = 1:N-1
        # xkp1 - axk - buk
        A[ci[i],xi[i+1]] = I(nx)
        A[ci[i],xi[i]] = -Ad
        A[ci[i],ui[i]] = -Bd
    end
    A[ci[end],xi[1]] = I(nx)
    b = zeros((N*nx))
    b[ci[end]] = [1;4;-8;.1;.2;-.3]

    G_up = spzeros(nz,nz)
    gi = [(nu)*(i-1) .+ (1:nu)  for i = 1:N-1]
    for i = 1:N-1
        G_up[gi[i],ui[i]] = I(nu)
    end
    G = [G_up;-G_up]
    h = 2*ones(2*nz)

    return Q,q,A,b,G,h
end

# here is a feasible QP
for N = 5:5:150
    local z, iters = CPEG.quadprog(create_MPC(N)...)
    @test typeof(z) == Vector{Float64}
    @test iters >= 1
    @test iters <= ev.qp_solver_opts.max_iters
end

# for N = 10:10:100
#     Q,q,A,b,G,h = create_MPC(N)
#     nx = 6
#     nu = 3
#     nz=(N*nx) + (N-1)*nu
#     h = [2*ones(nz);-3*ones(nz)]  # set this to zero and make it impossible to solve
#     @test_throws ErrorException CPEG.quadprog(Q,q,A,b,G,h)
# end

ext: jl


text:
# Check Arguments Macro ====================================================================

macro check_args(K, param, cond, desc=string(cond))
    quote
        if !($(esc(cond)))
            throw(ArgumentError(string(
                $(string(K)), ": ", $(string(param)), " = ", $(esc(param)), " does not ",
                "satisfy the constraint ", $(string(desc)), ".")))
        end
    end
end

#macro default_kernel_convert(K)
#    θ = fieldnames(eval(K))
#    K_name = esc(K)
#    conversion_call = Expr(:call, :($K_name{T}), [:(κ.$θₖ) for θₖ in θ]...)
#    return quote
#        function test(::Type{K}, κ::$K_name) where {K>:$K_name{T}} where T
#            return $conversion_call
#        end
#    end
#end

# Type Rules ===============================================================================

function promote_float(Tₖ::DataType...)
    if length(Tₖ) == 0
        return Float64
    end
    T = promote_type(Tₖ...)
    return T <: AbstractFloat ? T : Float64
end

# Common Functions =========================================================================

for orientation in (:row, :col)

    row_oriented = orientation == :row
    dim_obs, dim_param = row_oriented ? (1, 2) : (2, 1)
    NT, TN = row_oriented ? ('N', 'T') : ('T', 'N')

    @eval begin

        function dotvectors!(
                ::Val{$(Meta.quot(orientation))},
                xᵀx::Vector{T},
                X::Matrix{T}
            ) where {T<:AbstractFloat}
            if !(size(X, $dim_obs) == length(xᵀx))
                errorstring = string("Dimension mismatch on dimension ", $dim_obs)
                throw(DimensionMismatch(errorstring))
            end
            fill!(xᵀx, zero(T))
            for I in Base.Cartesian.CartesianIndices(size(X))
                xᵀx[I.I[$dim_obs]] += X[I]^2
            end
            xᵀx
        end

        @inline function dotvectors(σ::Val{$(Meta.quot(orientation))}, X::Matrix{T}) where {T<:AbstractFloat}
            dotvectors!(σ, Array{T}(undef, size(X,$dim_obs)), X)
        end

        function gramian!(
                ::Val{$(Meta.quot(orientation))},
                G::Matrix{T},
                X::Matrix{T},
                symmetrize::Bool
            ) where {T<:LinearAlgebra.BLAS.BlasReal}
            LinearAlgebra.BLAS.syrk!('U', $NT, one(T), X, zero(T), G)
            symmetrize ? LinearAlgebra.copytri!(G, 'U') : G
        end

        @inline function gramian!(
                ::Val{$(Meta.quot(orientation))},
                G::Matrix{T},
                X::Matrix{T},
                Y::Matrix{T}
            ) where {T<:LinearAlgebra.BLAS.BlasReal}
            LinearAlgebra.BLAS.gemm!($NT, $TN, one(T), X, Y, zero(T), G)
        end
    end
end


text:
#Logistic map
T(x) = mod(3.92*x*(1-x), 1)

#Arnold's maps
#epsilon = 0.7
#tau = 1.4
#T(x) = mod(x-(epsilon/(2*pi))*sin(2*pi*x)+tau, 1)

#Tsuda's map (from Chaotic Itinerancy in Random Dynamical System Related to Associative Memory Models)
#C = 0.1
#A = 0.07
#T(x) = mod(x+C+A*sin(4*pi*x), 1)

P = NoiseFourier.Fourier1D.assemble_matrix(T, 128; x_0 = 0, x_1 = 1)
D = NoiseFourier.Fourier1D.noise_matrix(0.005, 128)
M = D*P
v2, vinf = NoiseFourier.Fourier1D.rigorous_norm(M)

print(v2)
print(vinf)

function convert_to_float(M, n)
    Q = zeros(Complex{Float64},n,n)
    ze = 0
    zer = []
    max_nnz_rows = 0
    a = 0
    b = 0
    for i in 1:n
        for j in 1:n
            a = Float64(real((M[i,j])))
            b = Float64(imag((M[i,j])))
            Q[i,j] = a + im*b
            if Q[i,j] == 0
                ze = ze +1
            end
        end
        push!(zer, ze)
        ze = 0
    end
    max_nnz_rows = maximum(zer)
    return Q#, max_nnz_rows
end

function idx_to_freq(j,k)
    # k is the size of the matrix
    if j>=0 && j <=k ÷ 2+1
        return j-1
    end
    if j >k ÷ 2+1
        return j-k-1
    end
end

function measure(x,v)
    t = real(v[1])
    for i in 2:128
        a = 2*real(v[i])
        b = 2*imag(v[i])
        k = idx_to_freq(i,128)
        t += a*cos(2*pi*x*k) - b*sin(2*pi*x*k)
    end
    return t
end

Z = convert_to_float(M, 128)
v = eigvecs(Z)[:,128]
p = 50000
x = range(0,1, length = p)
y = zeros(p)
for i in 1:p
    y[i] = measure(x[i],v)
end
plot(x,y,ylims=(0,5),fmt = :png)#,xticks=25:5:75)


text:
# Update position with Mulitnomial Distribution
# Each step is taken with a hopping rate r
# Probability of hopping
# Mulitnomial steps
# Hopping rate
hopping_rate(r::Float64) = r

# Probability
prob(hopping_rate) = [hopping_rate,hopping_rate, 1-4*hopping_rate, hopping_rate,hopping_rate]
prob_no_flux(hopping_rate) = [hopping_rate,hopping_rate, 1-3*hopping_rate, hopping_rate]

# Multinomial
multidist(rate) = Multinomial(1,prob(hopping_rate(rate)))
multidist_no_flux(rate) = Multinomial(1,prob_no_flux(hopping_rate(rate)))

# Stepping functions
function stepping(rate::Float64)
    steps = [(-1, 0),(1, 0),(0, 0),(0, -1),(0, 1)]
    s = steps[findall(!iszero,rand(multidist(rate),1))]
    return Step2D(s[1][1],s[1][2])
end

function stepping_no_north(rate)
    steps = [(-1, 0),(1, 0),(0, 0),(0, -1)]
    s = steps[findall(!iszero,rand(multidist_no_flux(rate),1))]
    return Step2D(s[1][1],s[1][2])
end

function stepping_no_south(rate)
    steps = [(-1, 0),(1, 0),(0, 0),(0, 1)]
    s = steps[findall(!iszero,rand(multidist_no_flux(rate),1))]
    return Step2D(s[1][1],s[1][2])
end

function stepper(x::Int64,rate::Float64)
    if x == 0
        return stepping(rate)
    elseif x == 1
        return stepping(0.0)
    else
        error("input meant to be only 0 or 1 not $(x)")
    end
end

function stepper(x::Int64,s,rate::Float64)
    if x == 0
        return s(rate)
    elseif x == 1
        return s(0.0)
    else
        error("input meant to be only 0 or 1 not $(x)")
    end
end

Base.:+(w::Walker2D, s::Step2D) = Walker2D(w.x+s.x,w.y+s.y)
Base.:+(w1::Walker2D, w2::Walker2D) = Walker2D(w1.x+w2.x,w1.y+w2.y)
Base.:+(s₁::Step2D,s₂::Step2D) = Step2D(s₁.x+s₂.x,s₁.y+s₂.y)

function step(r::Float64;Δx = 1::Real)
    u = rand()
    if u < r
        s = 1
    elseif u < 2*r
        s = -1
    else
        s = 0
    end
    return Step1D(s*Δx)
end

Base.:+(w::Walker1D, s::Step1D) = Walker1D(w.x+s.x)
Base.:+(w1::Walker1D, w2::Walker1D) = Walker1D(w1.x+w2.x)
Base.:+(s₁::Step1D,s₂::Step1D) = Step1D(s₁.x+s₂.x)


text:
facts("PARAFAC2") do

srand(1)
r = 2
A = randn(r, 10)
T = Array{Float64, r}[(randn(4 * (i + 1), r) * A) for i in 1:5]
@time factors = parafac2(T, 2, compute_error=true)
@fact map(size, factors.factors) --> ([(size(t, 1), r) for t in T]...)
@fact length(factors.D) --> length(T)
@fact map(size, factors.D) --> fill((1, r), length(T))
@fact size(factors.A) --> (size(T[1], 2), r)
@fact rel_residue(factors) --> less_than(0.05)

end


text:
using Pkg.TOML, StaticArrays

@testset "ExcitationFields.jl" begin
    params = TOML.parsefile("ExcitationFields.toml")

    Hs = MPISimulations.initialize(MPISimulations.StaticGradientField,"StaticGradientField",params,Float64)
    H1 = MPISimulations.initialize(MPISimulations.SinusoidalField,"SinusoidalField1",params,Float64)
    H2 = MPISimulations.initialize(MPISimulations.SinusoidalField,"SinusoidalField2",params,Float64)
    H3 = MPISimulations.initialize(MPISimulations.SinusoidalField,"SinusoidalField3",params,Float64)
    r = randn(SVector{3,Float64})
    t = randn(Float64)

    excitation = MPISimulations.ExcitationField(Hs,H1)
    @test MPISimulations.magneticFieldStrength(excitation,r,t) ≈ MPISimulations.magneticFieldStrength(Hs,r,t) + MPISimulations.magneticFieldStrength(H1,r,t)

    excitation = MPISimulations.ExcitationField(Hs,H1,H2)
    @test MPISimulations.magneticFieldStrength(excitation,r,t) ≈ MPISimulations.magneticFieldStrength(Hs,r,t) + MPISimulations.magneticFieldStrength(H1,r,t) + MPISimulations.magneticFieldStrength(H2,r,t)

    excitation = MPISimulations.ExcitationField(Hs,H1,H2,H3)
    @test MPISimulations.magneticFieldStrength(excitation,r,t) ≈ MPISimulations.magneticFieldStrength(Hs,r,t) + MPISimulations.magneticFieldStrength(H1,r,t) + MPISimulations.magneticFieldStrength(H2,r,t) + MPISimulations.magneticFieldStrength(H3,r,t)
end


text:
"""
    PySet{T=PyObject}([o])

Wrap the Python set `o` (or anything satisfying the set interface) as a Julia set with elements of type `T`.

If `o` is not given, an empty set is created.
"""
struct PySet{T} <: AbstractSet{T}
    ref::PyRef
    PySet{T}(o) where {T} = new{T}(ispyref(o) ? PyRef(o) : pyset(PyRef, o))
    PySet{T}() where {T} = new{T}(PyRef())
end
PySet(o) = PySet{PyObject}(o)
PySet() = PySet{PyObject}()
export PySet

ispyreftype(::Type{<:PySet}) = true
pyptr(x::PySet) = begin
    ptr = x.ref.ptr
    if isnull(ptr)
        ptr = x.ref.ptr = C.PySet_New(C_NULL)
    end
    ptr
end
Base.unsafe_convert(::Type{CPyPtr}, x::PySet) = checknull(pyptr(x))
C.PyObject_TryConvert__initial(o, ::Type{T}) where {T<:PySet} =
    C.putresult(T(pyborrowedref(o)))

Base.iterate(x::PySet{T}, it::PyRef) where {T} = begin
    ptr = C.PyIter_Next(it)
    if ptr != C_NULL
        r = C.PyObject_Convert(ptr, T)
        C.Py_DecRef(ptr)
        checkm1(r)
        (C.takeresult(T), it)
    elseif C.PyErr_IsSet()
        pythrow()
    else
        nothing
    end
end
Base.iterate(x::PySet) = begin
    it = C.PyObject_GetIter(x)
    isnull(it) && pythrow()
    iterate(x, pynewref(it))
end

# Base.length(x::PySet) = @pyv `len($x)`::Int
Base.length(x::PySet) = Int(checkm1(C.PyObject_Length(x)))

Base.in(_v, x::PySet{T}) where {T} = begin
    v = tryconvertref(T, _v)
    v === PYERR() && pythrow()
    v === NOTIMPLEMENTED() && return false
    pycontains(x, v)
end

Base.push!(x::PySet{T}, v) where {T} = (@py `$x.add($(convertref(T, v)))`; x)

Base.delete!(x::PySet{T}, _v) where {T} = begin
    v = tryconvertref(T, _v)
    v === PYERR() && pythrow()
    v === NOTIMPLEMENTED() && return x
    @py `$x.discard($v)`
    x
end

Base.pop!(x::PySet{T}) where {T} = @pyv `$x.pop()`::T

Base.pop!(x::PySet{T}, _v) where {T} = begin
    v = tryconvert(T, _v)
    v === PYERR() && pythrow()
    (v !== NOTIMPLEMENTED() && v in x) ? (delete!(x, v); v) : error("not an element")
end

Base.pop!(x::PySet{T}, _v, d) where {T} = begin
    v = tryconvert(T, _v)
    v === PYERR() && pythrow()
    (v !== NOTIMPLEMENTED() && v in x) ? (delete!(x, v); v) : d
end

Base.empty!(x::PySet) = (@py `$x.clear()`; x)

Base.copy(x::PySet) = @pyv `$x.copy()`::typeof(x)


text:
module MAXENT
using Random
using Statistics
using LinearAlgebra
using Optim

export maxent

function entropy_like_term( entropy_intermediate_term::Array{Float64,1},
                            dsf::Array{Float64,1},
                            dsf_default::Array{Float64,1},
                            dfrequency1::Array{Float64,1},
                            dfrequency2::Array{Float64,1}
                            )
    #calculate Kullback-Leibler divergence using trapezoidal rule
    #see https://arxiv.org/pdf/1507.01012.pdf
    #set the intermediate result
    broadcast!((x,y) -> (abs(x)*log(abs(x)/y)),entropy_intermediate_term,dsf,dsf_default)
    #return Kullback-Leibler divergence
    -(dot(entropy_intermediate_term,dfrequency1) + dot(entropy_intermediate_term,dfrequency2))
end

function quality_of_fit( isf_m::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1},
                         dsf::Array{Float64,1},dsf_default::Array{Float64,1},
                         entropy_intermediate_term::Array{Float64,1},
                         dfrequency1::Array{Float64,1},
                         dfrequency2::Array{Float64,1},
                         regularization_constant::Float64)
    chi_squared_term = mean(broadcast((x,y,z) -> ((x - y)/z)^2,isf,isf_m,isf_error));
    entropy_like_term_fit = regularization_constant*entropy_like_term( entropy_intermediate_term,
                                                                       dsf, dsf_default, dfrequency1,
                                                                       dfrequency2);
    chi_squared_term_fit/2 - entropy_like_term_fit
end

function quality_of_fit!( isf_m::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1},
                          dsf::Array{Float64,1},dsf_default::Array{Float64,1},
                          entropy_intermediate_term::Array{Float64,1},
                          dfrequency1::Array{Float64,1},
                          dfrequency2::Array{Float64,1},
                          regularization_constant::Float64)
    chi_squared_term_fit = mean(broadcast!((x,y,z) -> ((x - y)/z)^2,isf_m,isf,isf_m,isf_error));
    entropy_like_term_fit = regularization_constant*entropy_like_term( entropy_intermediate_term,
                                                                       dsf, dsf_default, dfrequency1,
                                                                       dfrequency2);
    chi_squared_term_fit/2 - entropy_like_term_fit
end

function set_isf_term(imaginary_time::Array{Float64,1},frequency_bins::Array{Float64,1},beta::Float64)
    b = beta;
    isf_term = Array{Float64,2}(undef,size(imaginary_time,1),size(frequency_bins,1));
    for i in 1:size(imaginary_time,1)
        t = imaginary_time[i];
        for j in 1:size(frequency_bins,1)
            f = frequency_bins[j];
            isf_term[i,j] = (exp(-t*f) + exp(-(b - t)*f));
        end
    end
    isf_term
end

function set_isf_trapezoidal!( isf_m::Array{Float64,1}, isf_m2::Array{Float64,1},
                               isf_term::Array{Float64,2}, isf_term2::Array{Float64,2},
                               dsf::Array{Float64,1}
                               )
    mul!(isf_m,isf_term,abs.(dsf));
    mul!(isf_m2,isf_term2,abs.(dsf));
    isf_m .+= isf_m2;
    isf_m
end

function maxent(dsf::Array{Float64,1},dsf_default::Array{Float64,1},isf::Array{Float64,1},isf_error::Array{Float64,1},
                frequency::Array{Float64,1},imaginary_time::Array{Float64,1};
                temperature::Float64 = 1.2,
                regularization_constant::Float64 = 0.0,
                number_of_iterations::Int64 = 1000,
                allow_f_increases::Bool = false,
                number_of_blas_threads::Int64 = 0)
    #_bts = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
    #println("Number of BLAS threads: $(_bts)");
    if number_of_blas_threads > 0
        LinearAlgebra.BLAS.set_num_threads(number_of_blas_threads)
    end

    moment0 = isf[1];
    beta = 1/temperature;

    isf_term = set_isf_term(imaginary_time,frequency,beta);
    isf_term2 = copy(isf_term);

    df = frequency[2:end] .- frequency[1:size(frequency,1) - 1];
    dfrequency1 = zeros(size(frequency,1));
    dfrequency2 = zeros(size(frequency,1));
    for i in 1:(size(frequency,1) - 1)
        dfrequency1[i] = df[i]/2
        dfrequency2[i+1] = df[i]/2
    end
    isf_term .*= dfrequency1';
    isf_term2 .*= dfrequency2';

    isf_m = Array{Float64,1}(undef,size(isf,1));
    isf_m2 = Array{Float64,1}(undef,size(isf,1));

    entropy_intermediate_term = Array{Float64,1}(undef,size(dsf,1));

    result = Optim.optimize(
                 x -> quality_of_fit!(
                     set_isf_trapezoidal!(
                         isf_m, isf_m2, isf_term, isf_term2, x ),
                     isf, isf_error, x, dsf_default, entropy_intermediate_term,
                     dfrequency1, dfrequency2, regularization_constant ),
                 dsf, BFGS(), Optim.Options(
                     iterations=number_of_iterations,
                     allow_f_increases=allow_f_increases));

    set_isf_trapezoidal!(isf_m, isf_m2, isf_term, isf_term2, result.minimizer );
    chi_squared_term = mean(broadcast((x,y,z) -> ((x - y)/z)^2,isf,isf_m,isf_error));
    result, chi_squared_term
end
end


text:
# ---
# title: 1382. Balance a Binary Search Tree
# id: problem1382
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search Tree
# link: <https://leetcode.com/problems/balance-a-binary-search-tree/description/>
# hidden: true
# ---
#
# Given a binary search tree, return a **balanced** binary search tree with the
# same node values.
#
# A binary search tree is _balanced_ if and only if the depth of the two
# subtrees of every node never differ by more than 1.
#
# If there is more than one answer, return any of them.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/08/22/1515_ex1.png)![](https://assets.leetcode.com/uploads/2019/08/22/1515_ex1_out.png)**
#
#
#
#     Input: root = [1,null,2,null,3,null,4,null,null]
#     Output: [2,1,3,null,null,null,4]
#     Explanation: This is not the only correct answer, [3,1,4,null,2,null,null] is also correct.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is between `1` and `10^4`.
#   * The tree nodes will have distinct values between `1` and `10^5`.
#
#
## @lc code=start
using LeetCode

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


text:
"""

Посчитать число всех горизонтальных прямолинейных перегородок (вертикальных - нет)

"""
#= пока нет границы сверху
    пока нет границы на Ost
        пока сверху перегородка
            идти на Ost
        перегородка += 1
        пока сверху нет перегородки
        идти на Ost

    идти на Nord
    повторить для West
=#

using HorizonSideRobots
originRobot = Robot(animate = true, "20example.sit")
include("../AbstractType.jl")
include("../functions.jl")

function mod_moveToBorder!(robot::Robot, moveSide::HorizonSide, borderSide::HorizonSide)::Integer
    bordersInLine = 0
    state = 0
    while (!isborder(robot, moveSide))
        move!(robot, moveSide)
        if !isborder(robot, borderSide)
            if state == 1
                state = 0
                bordersInLine += 1
            end
        else
            state = 1
        end
    end
    return bordersInLine
end

function mod_moveLikeSnake!(robot::Robot, moveSides::NTuple{2,HorizonSide}, borderSide::HorizonSide)::Integer
    bordersInLine = 0
    for side in moveSides
        bordersInLine += mod_moveToBorder!(robot, side, borderSide)
        println(bordersInLine)
        if !(isborder(robot, borderSide))
            move!(robot, borderSide)
        else
            break
        end
    end
    return bordersInLine
end

function master20!(robot::Robot)
    numberOfBorders = 0
    backPath = moving_in_angle!(robot, West, Sud)
    while !isborder(robot, Nord)
        numberOfBorders += mod_moveLikeSnake!(robot, (Ost, West), Nord)
    end

    moving_in_angle!(robot, West, Sud)
    moving_back_to_start!(robot, backPath)

    return numberOfBorders
end

master20!(originRobot)

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
#=
@author : Spencer Lyon
          Victoria Gregory

@date: 07/09/2014
=#

using Plots
pyplot()

const r = 0.05
const beta = 1.0 / (1.0 + r)
const T = 20  # Time horizon
const S = 5   # Impulse date
const sigma1 = 0.15
const sigma2 = 0.15

function time_path(permanent=false)
    w1 = zeros(T+1)
    w2 = zeros(T+1)
    b = zeros(T+1)
    c = zeros(T+1)

    if permanent === false
        w2[S+2] = 1.0
    else
        w1[S+2] = 1.0
    end

    for t=2:T
        b[t+1] = b[t] - sigma2 * w2[t]
        c[t+1] = c[t] + sigma1 * w1[t+1] + (1 - beta) * sigma2 * w2[t+1]
    end

    return b, c
end

function main()
    L = 0.175

    b1, c1 = time_path(false)
    b2, c2 = time_path(true)
    p = plot(0:T, [c1 c2 b1 b2], layout=(2, 1),
             color=[:green :green :blue :blue],
             label=["consumption" "consumption" "debt" "debt"])
    t = ["impulse-response, transitory income shock"
         "impulse-response, permanent income shock"]
    plot!(title=t', xlabel="Time", ylims=(-L, L), legend=[:topright :bottomright])
    vline!([S S], color=:black, layout=(2, 1), label="")

    return p
end

max_stars_repo_name: mwhchen/quantecon

max_issues_repo_name: mwhchen/quantecon

max_forks_repo_name: mwhchen/quantecon


text:
using MultiModalMuSig
using LinearAlgebra
using SpecialFunctions
using Test

K = [2, 3]
α = [0.1, 0.1]
features = Matrix{Int}[
    # modality 1
    [
        1 1;
        1 2;
        2 1;
        2 2
    ],
    # modality 2
    [
        1 1;
        1 2;
        2 1;
        2 2
    ]
]
X = Vector{Matrix{Int}}[
    # document 1
    [
        # modality 1
        [
            1 5;
            2 8
        ],
        # modality 2
        [
            1 2;
            2 5
        ]
    ],
    # document 2
    [
        # modality 1
        [
            3 4;
            4 9
        ],
        # modality 2
        [
            3 4;
            4 6
        ]
    ],
]

@testset "constructor" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    @test model.D == 2
    @test model.N == [[13, 7], [13, 10]]
    @test model.M == 2
    @test model.I == [2, 2]
    @test model.J == [[2, 2], [2, 2]]
    @test model.V == [4, 4]

    @test length(model.μ) == 5
    @test size(model.Σ) == (5, 5)
    @test size(model.invΣ) == (5, 5)

    @test length(model.ζ) == 2
    @test length(model.ζ[1]) == 2
    @test sum(model.θ[1][1], dims=1) ≈ ones(2)'
    @test length(model.λ[1]) == 5
    @test model.ν[1] == ones(5)

    @test length(model.γ) == 2
    @test length(model.γ[1]) == 2
    @test length(model.γ[1][1]) == 2
    @test length(model.γ[1][1][1]) == 2
    @test all(model.γ[1][1][1] .> 0) == true
end

@testset "calculate_Ndivζ" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.ζ = [[2, 3], [4,5]]

    Ndivζ = [
        sum(X[1][1][:, 2]) / model.ζ[1][1],
        sum(X[1][1][:, 2]) / model.ζ[1][1],
        sum(X[1][2][:, 2]) / model.ζ[1][2],
        sum(X[1][2][:, 2]) / model.ζ[1][2],
        sum(X[1][2][:, 2]) / model.ζ[1][2]
    ]

    @test MultiModalMuSig.calculate_Ndivζ(model, 1) ≈ Ndivζ
end

@testset "calculate_sumθ" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    θ = Matrix{Float64}[]
    push!(θ, [0.4 0.1; 0.6 0.9])
    push!(θ, [0.3 0.4; 0.3 0.5; 0.4 0.1])
    model.θ[1] = θ

    sumθ = [
        X[1][1][1, 2] * θ[1][1, 1] + X[1][1][2, 2] * θ[1][1, 2],
        X[1][1][1, 2] * θ[1][2, 1] + X[1][1][2, 2] * θ[1][2, 2],
        X[1][2][1, 2] * θ[2][1, 1] + X[1][2][2, 2] * θ[2][1, 2],
        X[1][2][1, 2] * θ[2][2, 1] + X[1][2][2, 2] * θ[2][2, 2],
        X[1][2][1, 2] * θ[2][3, 1] + X[1][2][2, 2] * θ[2][3, 2],
    ]

    @test MultiModalMuSig.calculate_sumθ(model, 1) ≈ sumθ
end

@testset "update_λ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)

    λ = Float64[1, 2, 3, 4, 1]
    model.λ[1] = copy(λ)

    MultiModalMuSig.update_λ!(model, 1)
    @test model.λ[1] ≉ λ
end

@testset "ν_objective" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)

    μ = Float64[1, 1, 2, 2, 1]
    invΣ = Matrix{Float64}(I, sum(K), sum(K))
    λ = Float64[1, 2, 3, 4, 1]
    ν = Float64[1, 1, 1, 2, 1]
    ζ = Float64[2, 1]

    L = (
        -0.5 * tr(diagm(0 => ν) * invΣ) -
        sum(X[1][1][:, 2]) / ζ[1] *
            (exp(λ[1] + 0.5ν[1]) + exp(λ[2] + 0.5ν[2])) -
        sum(X[1][2][:, 2]) / ζ[2] *
            (exp(λ[3] + 0.5ν[3]) + exp(λ[4] + 0.5ν[4]) + exp(λ[5] + 0.5ν[5])) +
        0.5 * (log(ν[1]) + log(ν[2]) + log(ν[3]) + log(ν[4]) + log(ν[5]))
    )
    ∇ν = Array{Float64}(undef, sum(K))
    Ndivζ = vcat([fill(model.N[1][m] / ζ[m], K[m]) for m in 1:model.M]...)
    objective = MultiModalMuSig.ν_objective(ν, ∇ν, λ, Ndivζ, μ, invΣ)
    @test objective ≈ L

    grad = -0.5 * diag(invΣ)
    grad[1] += (
        -sum(X[1][1][:, 2]) / (2 * ζ[1]) * exp(λ[1] + 0.5ν[1]) + 1 / (2ν[1])
    )
    grad[2] += (
        -sum(X[1][1][:, 2]) / (2 * ζ[1]) * exp(λ[2] + 0.5ν[2]) + 1 / (2ν[2])
    )
    grad[3] += (
        -sum(X[1][2][:, 2]) / (2 * ζ[2]) * exp(λ[3] + 0.5ν[3]) + 1 / (2ν[3])
    )
    grad[4] += (
        -sum(X[1][2][:, 2]) / (2 * ζ[2]) * exp(λ[4] + 0.5ν[4]) + 1 / (2ν[4])
    )
    grad[5] += (
        -sum(X[1][2][:, 2]) / (2 * ζ[2]) * exp(λ[5] + 0.5ν[5]) + 1 / (2ν[5])
    )
    @test ∇ν ≈ grad

    model.μ = μ
    model.λ[1] = λ
    model.ν[1] = ν
    model.ζ[1] = ζ
    MultiModalMuSig.update_ν!(model, 1)
    @test (model.ν[1] .> 0.0) == fill(true, sum(K))
    @test (model.λ[1] .< 100.0) == fill(true, sum(K))
end

@testset "update_ζ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]
    model.ν = [[1,1,1,2,1], [1,3,1,2,1]]
    MultiModalMuSig.update_ζ!(model, 1)

    ζ = Float64[exp(1.5) + exp(2.5), exp(3.5) + exp(5) + exp(1.5)]
    @test model.ζ[1] ≈ ζ
end

@testset "update_θ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]
    model.γ = [
        [
            [[0.1, 0.2],[0.1, 1.0]],
            [[0.1, 0.1],[1.0, 1.0]]
        ],
        [
            [[0.5, 0.5],[1.0, 1.5]],
            [[1.0, 2.0],[2.0, 3.0]],
            [[1.0, 5.0],[5.0, 2.0]]
        ]
    ]
    MultiModalMuSig.update_Elnϕ!(model)
    MultiModalMuSig.update_θ!(model, 1)

    θ = Array{Float64}(undef, 2, 2)
    θ[1, 1] = exp(1 + digamma(0.1) - digamma(0.3) + digamma(0.1) - digamma(1.1))
    θ[2, 1] = exp(2 + digamma(0.1) - digamma(0.2) + digamma(1.0) - digamma(2.0))
    θ[1, 2] = exp(1 + digamma(0.1) - digamma(0.3) + digamma(1.0) - digamma(1.1))
    θ[2, 2] = exp(2 + digamma(0.1) - digamma(0.2) + digamma(1.0) - digamma(2.0))
    θ[:, 1] ./= sum(θ[:, 1])
    θ[:, 2] ./= sum(θ[:, 2])

    @test sum(model.θ[1][1], dims=1) ≈ ones(size(X[1][1])[1])'
    @test model.θ[1][1] ≈ θ
    @test any(model.θ[1][1] .< 0.0) == false

    MultiModalMuSig.update_θ!(model, 2)
    θ = Array{Float64}(undef, 3, 2)
    θ[1, 1] = exp(1 + digamma(0.5) - digamma(1.0) + digamma(1.0) - digamma(2.5))
    θ[2, 1] = exp(4 + digamma(2.0) - digamma(3.0) + digamma(2.0) - digamma(5.0))
    θ[3, 1] = exp(2 + digamma(5.0) - digamma(6.0) + digamma(5.0) - digamma(7.0))
    θ[1, 2] = exp(1 + digamma(0.5) - digamma(1.0) + digamma(1.5) - digamma(2.5))
    θ[2, 2] = exp(4 + digamma(2.0) - digamma(3.0) + digamma(3.0) - digamma(5.0))
    θ[3, 2] = exp(2 + digamma(5.0) - digamma(6.0) + digamma(2.0) - digamma(7.0))
    θ[:, 1] ./= sum(θ[:, 1])
    θ[:, 2] ./= sum(θ[:, 2])

    @test model.θ[2][2] ≈ θ
end

@testset "update_μ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]

    MultiModalMuSig.update_μ!(model)

    @test model.μ ≈ [1.5, 2.5, 2.0, 4.0, 1.5]
end

@testset "update_Σ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.λ = [[1,2,3,4,1], [2,3,1,4,2]]
    model.ν = [[1,1,1,2,1], [1,3,1,2,1]]
    model.μ = [1, 1, 2, 2, 1]

    MultiModalMuSig.update_Σ!(model)

    diff1 = model.λ[1] .- model.μ
    diff2 = model.λ[2] .- model.μ
    Σ = 0.5 * (
        diagm(0 => model.ν[1]) + diagm(0 => model.ν[2]) +
        (diff1 * diff1') + (diff2 * diff2')
    )
    @test model.Σ ≈ Σ
    @test model.invΣ ≈ inv(Σ)
end

@testset "update_γ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.θ[1][1] = [0.4 0.1; 0.6 0.9]
    model.θ[2][1] = [0.3 0.5; 0.7 0.5]
    MultiModalMuSig.update_γ!(model)

    γ1 = [0.1 + 5 * 0.4 + 8 * 0.1, 0.1 + 4 * 0.3 + 9 * 0.5]
    γ2 = [0.1 + 5 * 0.4 + 4 * 0.3, 0.1 + 8 * 0.1 + 9 * 0.5]
    @test model.γ[1][1][1] ≈ γ1
    @test model.γ[1][1][2] ≈ γ2
end

@testset "update_Elnϕ!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    model.γ[1][1][1] .= [1, 2]

    MultiModalMuSig.update_Elnϕ!(model)

    @test model.Elnϕ[1][1][1][1] ≈ digamma(1) - digamma(3)
end

# TODO
@testset "update_α!" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)

    sum_Elnϕ = sum(model.Elnϕ[1][1][1]) + sum(model.Elnϕ[1][2][1])
    L = K[1] * (lgamma(2α[1]) - 2lgamma(α[1])) + α[1] * sum_Elnϕ
    grad = 2K[1] * (digamma(2α[1]) - digamma(α[1])) + sum_Elnϕ

    ∇α = Array{Float64}(undef, 1)
    @test MultiModalMuSig.α_objective(model.α[1][1:1], ∇α, sum_Elnϕ, K[1], 2) ≈ L
    @test ∇α[1] ≈ grad

    sum_Elnϕ = sum(model.Elnϕ[2][1][2]) + sum(model.Elnϕ[2][2][2])
    L_before = K[2] * (lgamma(2α[2]) - 2lgamma(α[2])) + α[2] * sum_Elnϕ
    MultiModalMuSig.update_α!(model)
    L_after = (
        K[2] * (lgamma(2model.α[2][2]) - 2lgamma(model.α[2][2])) +
        model.α[2][2] * sum_Elnϕ
    )
    @test model.α[1] ≉ fill(α[1], 2)
    @test model.α[2] ≉ fill(α[2], 2)
    @test L_after > L_before
end

@testset "calculate_ElnPϕ" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnPϕ = 0.0
    #@test MultiModalMuSig.calculate_ElnPϕ(model) ≈ ElnPϕ
end

@testset "calculate_ElnPη" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnPη = 0.0
    #@test MultiModalMuSig.calculate_ElnPη(model) ≈ ElnPη
end

@testset "calculate_ElnPZ" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnPZ = 0.0
    #@test MultiModalMuSig.calculate_ElnPZ(model) ≈ ElnPZ
end

@testset "calculate_ElnPX" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnPX = 0.0
    #@test MultiModalMuSig.calculate_ElnPX(model) ≈ ElnPX
end

@testset "calculate_ElnQϕ" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnQϕ = 0.0
    #@test MultiModalMuSig.calculate_ElnQϕ(model) ≈ ElnQϕ
end

@testset "calculate_ElnQη" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnQη = 0.0
    #@test MultiModalMuSig.calculate_ElnQη(model) ≈ ElnQη
end

@testset "calculate_ElnQZ" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ElnQZ = 0.0
    #@test MultiModalMuSig.calculate_ElnQZ(model) ≈ ElnQZ
end

@testset "calculate_elbo" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    @test MultiModalMuSig.calculate_elbo(model) < 0.0
end

@testset "fit" begin
    model = MultiModalMuSig.IMMCTM(K, α, features, X)
    ll = MultiModalMuSig.fit!(model, maxiter=1, verbose=false)
    @test length(ll) == 1
    @test length(ll[1]) == 2
end

@testset "loglikelihoods"  begin
    η = [[1.0, 2.0], [2.0, 3.0]]
    θ = [exp.(η[d]) ./ sum(exp.(η[d])) for d in 1:2]

    γ = [
        [[0.1, 0.2],[0.1, 1.0]],
        [[0.1, 0.1],[1.0, 1.0]]
    ]
    ϕ = [[γ[k][i] ./ sum(γ[k][i]) for i in 1:2] for k in 1:2]

    Xm1 = [X[d][1] for d in 1:2]

    sum_ll = (
        Xm1[1][1, 2] * log(
            θ[1][1] * ϕ[1][1][1] * ϕ[1][2][1] +
            θ[1][2] * ϕ[2][1][1] * ϕ[2][2][1]
        ) +
        Xm1[1][2, 2] * log(
            θ[1][1] * ϕ[1][1][1] * ϕ[1][2][2] +
            θ[1][2] * ϕ[2][1][1] * ϕ[2][2][2]
        ) +
        Xm1[2][1, 2] * log(
            θ[2][1] * ϕ[1][1][2] * ϕ[1][2][1] +
            θ[2][2] * ϕ[2][1][2] * ϕ[2][2][1]
        ) +
        Xm1[2][2, 2] * log(
            θ[2][1] * ϕ[1][1][2] * ϕ[1][2][2] +
            θ[2][2] * ϕ[2][1][2] * ϕ[2][2][2]
        )
    )
    ll = sum_ll / sum(sum(X[d][1][:, 2]) for d in 1:2)
    res = MultiModalMuSig.calculate_modality_loglikelihood(
        Xm1, η, ϕ, features[1]
    )

    @test res ≈ ll
end


text:
module Agent
using Distributions
export AbstractAgent, QLearningAgent, spawn, update, policy

abstract type AbstractAgent end

mutable struct QLearningAgent <: AbstractAgent
    q::Array{Real, 1}
    α::Real
    β::Real
    function QLearningAgent(α::Real, β::Real, k::Int64)
        q = Array{Real, 1}(zeros(k))
        new(q, α, β)
    end
end

function update(agent::QLearningAgent, action::Real, reward::Real)
    agent.q[action] += agent.α * (reward - agent.q[action])
end

function policy(agent::QLearningAgent)
    qmax = maximum(agent.q)
    qexp = exp.(agent.β .* (agent.q .- qmax))
    return qexp ./ sum(qexp)
end

mutable struct ParticleFilterAgent <: AbstractAgent
    npartcles::Real
    particles::Array{Real, 2}
    q::Array{Real, 1}
    σ2::Real
    β::Real
    function ParticleFilterAgent(npartcles::Real, σ2::Real, β::Real, k::Real)
        particles = rand(Normal(), (k, npartcles))
        probs = sigmoid.(particles)
        q = sum(probs, dims=2)[:, 1] / size(probs, 2)
        new(npartcles, particles, q, σ2, β)
    end
end

function sigmoid(x::Real)
    1 / (1 + exp(-x))
end

function choice(x::Array{Real, 1}, p::Array{Real, 1}, n::Real)
    ret = Array{Real, 1}(undef, n)
    p = cumsum(p)
    for i in 1:n
        th = rand()
        for (x_, p_) in zip(x, p)
           if th <= p_
               ret[i] = x_
               break
           end
       end
    end
    return ret
end

function update(agent::ParticleFilterAgent, action::Real, reward::Real)
    noiz = Array{Real, 1}(rand(Normal(0, agent.σ2), agent.npartcles))
    particles = Array{Real, 1}(agent.particles[action, :] + noiz)
    probs = sigmoid.(particles)
    loglikelihood = pdf.(Bernoulli.(probs), reward)
    lnormalize = Array{Real, 1}(loglikelihood ./ sum(loglikelihood))
    agent.q[action] = sum(lnormalize .* probs)
    agent.particles[action, :] = choice(particles, lnormalize, agent.npartcles)
end

function policy(agent::ParticleFilterAgent)
    qmax = maximum(agent.q)
    qexp = exp.(agent.β .* (agent.q .- qmax))
    p = qexp ./ sum(qexp)
    p ./ sum(p)
end

end#modulie


text:
module ChemostatPlots

import Chemostat
const Ch = Chemostat
const ChU = Chemostat.Utils
const ChLP = Chemostat.LP
import LinearAlgebra: normalize
using Plots
using Distributions

include("plot_marginal.jl")
include("plot_projection2D.jl")

end


text:
# Autogenerated wrapper script for Shaderc_jll for x86_64-w64-mingw32-cxx11
export glslangValidator, glslc, libSPIRV_Tools_shared, libshaderc_shared, spirv_as, spirv_cfg, spirv_dis, spirv_link, spirv_lint, spirv_opt, spirv_reduce, spirv_remap, spirv_val

JLLWrappers.@generate_wrapper_header("Shaderc")
JLLWrappers.@declare_library_product(libSPIRV_Tools_shared, "libSPIRV-Tools-shared.dll")
JLLWrappers.@declare_library_product(libshaderc_shared, "libshaderc_shared.dll")
JLLWrappers.@declare_executable_product(glslangValidator)
JLLWrappers.@declare_executable_product(glslc)
JLLWrappers.@declare_executable_product(spirv_as)
JLLWrappers.@declare_executable_product(spirv_cfg)
JLLWrappers.@declare_executable_product(spirv_dis)
JLLWrappers.@declare_executable_product(spirv_link)
JLLWrappers.@declare_executable_product(spirv_lint)
JLLWrappers.@declare_executable_product(spirv_opt)
JLLWrappers.@declare_executable_product(spirv_reduce)
JLLWrappers.@declare_executable_product(spirv_remap)
JLLWrappers.@declare_executable_product(spirv_val)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libSPIRV_Tools_shared,
        "bin\\libSPIRV-Tools-shared.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libshaderc_shared,
        "bin\\libshaderc_shared.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        glslangValidator,
        "bin\\glslangValidator.exe",
    )

    JLLWrappers.@init_executable_product(
        glslc,
        "bin\\glslc.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_as,
        "bin\\spirv-as.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_cfg,
        "bin\\spirv-cfg.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_dis,
        "bin\\spirv-dis.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_link,
        "bin\\spirv-link.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_lint,
        "bin\\spirv-lint.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_opt,
        "bin\\spirv-opt.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_reduce,
        "bin\\spirv-reduce.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_remap,
        "bin\\spirv-remap.exe",
    )

    JLLWrappers.@init_executable_product(
        spirv_val,
        "bin\\spirv-val.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
immutable Laplace <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function Laplace(l::Real, s::Real)
        s > zero(s) || error("scale must be positive")
        new(float64(l), float64(s))
    end
end

Laplace(location::Real) = Laplace(location, 1.0)
Laplace() = Laplace(0.0, 1.0)

typealias Biexponential Laplace

## Support
@distr_support Laplace -Inf Inf

## Properties
mean(d::Laplace) = d.location
median(d::Laplace) = d.location
mode(d::Laplace) = d.location

std(d::Laplace) = sqrt2 * d.scale
var(d::Laplace) = 2.0 * d.scale^2
skewness(d::Laplace) = 0.0
kurtosis(d::Laplace) = 3.0

entropy(d::Laplace) = log(2.0 * d.scale) + 1.0

## Functions
pdf(d::Laplace, x::Float64) = 0.5exp(-abs(x - d.location)/d.scale) / d.scale
logpdf(d::Laplace, x::Float64) = -log(2.0 * d.scale) - abs(x - d.location) / d.scale

cdf(d::Laplace, x::Float64) = x < d.location ? 0.5*exp((x-d.location)/d.scale) : 0.5*(2.0-exp((d.location-x)/d.scale))
ccdf(d::Laplace, x::Float64) = x < d.location ? 0.5-0.5*expm1((x-d.location)/d.scale) : 0.5*exp((d.location-x)/d.scale)
logcdf(d::Laplace, x::Float64) = x < d.location ? loghalf + ((x-d.location)/d.scale) : loghalf + log2mexp((d.location-x)/d.scale)
logccdf(d::Laplace, x::Float64) = x < d.location ? loghalf + log2mexp((x-d.location)/d.scale) : loghalf + ((d.location-x)/d.scale)

quantile(d::Laplace, p::Float64) = p < 0.5 ? d.location + d.scale*log(2.0*p) : d.location - d.scale*log(2.0*(1.0-p))
cquantile(d::Laplace, p::Float64) = p >= 0.5 ? d.location + d.scale*log(2.0*(1.0-p)) : d.location - d.scale*log(2.0*p)
invlogcdf(d::Laplace, lp::Float64) = lp < loghalf ? d.location + d.scale*(logtwo + lp) : d.location - d.scale*(logtwo + log1mexp(lp))
invlogccdf(d::Laplace, lp::Float64) = lp >= loghalf ? d.location + d.scale*(logtwo + log1mexp(lp)) : d.location - d.scale*(logtwo + lp)

function gradlogpdf(d::Laplace, x::Float64)
    d.location != x || error("Gradient is undefined at the location point")
    x > d.location ? - 1.0 / d.scale : 1.0 / d.scale
end

function mgf(d::Laplace, t::Real)
    st = d.scale*t
    exp(t * d.location) / ((1.0-st)*(1.0+st))
end
function cf(d::Laplace, t::Real)
    st = d.scale*t
    exp(im * t * d.location) / ((1.0-st)*(1.0+st))
end

## Sampling
function rand(d::Laplace)
    er = Base.Random.randmtzig_exprnd()
    randbool() ? d.location - d.scale * er : d.location + d.scale * er
end

## Fitting
function fit_mle(::Type{Laplace}, x::Array)
    a = median(x)
    Laplace(a, mad(x, a))
end


text:
module ModuleTransform
using MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleDataNormalization
using MultipleViewGeometry.ModuleTypes
using StaticArrays, LinearAlgebra
export transform
include("transform.jl")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewaySku <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    tier::Any # spec type: Union{ Nothing, String } # spec name: tier
    capacity::Any # spec type: Union{ Nothing, Int32 } # spec name: capacity

    function ApplicationGatewaySku(;name=nothing, tier=nothing, capacity=nothing)
        o = new()
        validate_property(ApplicationGatewaySku, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ApplicationGatewaySku, Symbol("tier"), tier)
        setfield!(o, Symbol("tier"), tier)
        validate_property(ApplicationGatewaySku, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        o
    end
end # type ApplicationGatewaySku

const _property_map_ApplicationGatewaySku = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("tier")=>Symbol("tier"), Symbol("capacity")=>Symbol("capacity"))
const _property_types_ApplicationGatewaySku = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("tier")=>"String", Symbol("capacity")=>"Int32")
Base.propertynames(::Type{ ApplicationGatewaySku }) = collect(keys(_property_map_ApplicationGatewaySku))
Swagger.property_type(::Type{ ApplicationGatewaySku }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewaySku[name]))}
Swagger.field_name(::Type{ ApplicationGatewaySku }, property_name::Symbol) =  _property_map_ApplicationGatewaySku[property_name]

const _allowed_ApplicationGatewaySku_name = ["Standard_Small", "Standard_Medium", "Standard_Large", "WAF_Medium", "WAF_Large"]

const _allowed_ApplicationGatewaySku_tier = ["Standard", "WAF"]

function check_required(o::ApplicationGatewaySku)
    true
end

function validate_property(::Type{ ApplicationGatewaySku }, name::Symbol, val)
    if name === Symbol("name")
        Swagger.validate_param(name, "ApplicationGatewaySku", :enum, val, _allowed_ApplicationGatewaySku_name)
    end
    if name === Symbol("tier")
        Swagger.validate_param(name, "ApplicationGatewaySku", :enum, val, _allowed_ApplicationGatewaySku_tier)
    end
end


text:
# test for naiveunify.jl
using Test

#include("naiveunify.jl")

@testset "disagreement" begin
 @test disagreement( :(x), :(x)) == ()
 @test disagreement( :(x), :(y)) == (:x,:y)
 @test disagreement( :(x), :(y)) == (:x,:y)
 @test disagreement( :(x), :(y)) == (:x,:y)
 @test disagreement( :(f(x)), :(f(x))) == ()
 @test disagreement( :(f(x)), :(g(x))) == (:(f(x)), :(g(x)))
 @test disagreement( :(f(x)), :(f(g(x)))) == (:x, :(g(x)))
 @test disagreement( :(f(x,y,z)), :(f(x,y,z))) == ()
 @test disagreement( :(f(x,y,x)), :(f(x,y,z))) == (:x,:z)
 @test disagreement( :(f(x,z,x)), :(f(x,y,z))) == (:z,:y)
 @test disagreement( :(f(x,y,y)), :(f(z,y,z))) == (:x,:z)
 @test disagreement( :(f(g(x),y,x)), :(f(x,y,z))) == (:(g(x)),:x)
 @test disagreement( :(f(g(x),y,a)), :(f(g(x),d,z))) == (:y,:d)
 @test disagreement( :(f(g(x),y,a)), :(f(g(h(a)),z,z))) == (:x,:(h(a)))
 @test disagreement( :(f(g(x),h(y,u),x)), :(f(g(x),h(y,hen(i)),z))) == (:u,:(hen(i)))

end

@testset "naiveunify" begin
 @test_throws ICMP unify([],:(P()),:(Q()))

 @test unify([:x],:(),:())== [:x]
 @test unify([:x],:(P()),:(P())) == [:x]
 @test_throws ICMP unify([:x],:(P()),:(Q()))

 @test unify([:x],:(P(x)), :(P(x))) == [:x]
 @test unify([:y],:(P(x)), :(P(x))) == [:y]

 @test unify([:x],:(P(x)), :(P(y))) == [:y]
 @test unify([:y],:(P(x)), :(P(y))) == [:x]

 @testset "with number" begin
  @test unify([], :(f(2)), :(f(2))) == []
  @test unify([:x], :(f(x)), :(f(2))) == [2]
  @test unify([:x], :(f(3)), :(f(x))) == [3]
  @test unify([:x,:y], :(f(3,y)), :(f(x,23))) == [3,23]
 end

# in this case, x and y are constants
 @test_throws ICMP unify([:z],:(P(x)), :(P(y)))

 @test unify([],:(P(a)), :(P(a))) == []
 @test_throws ICMP unify([],:(P(a)), :(P(b)))

 (σ=unify([:x,:y],:(P(x)),:(P(y))); @test σ==[:x,:x] || σ==[:y,:y])

 @test unify([:x,:y],:(P(x)),:(P(a))) == [:a,:y]
 @test unify([:x,:y],:(P(f(x))),:(P(f(a)))) == [:a,:y]

 @test unify([:x,:y],:(P(x)),:(P(f(a)))) == [:(f(a)),:y]

 (σ = unify([:x,:y],:(P(x,f(a,x))),:(P(y,f(a,y)))); @test σ == [:y,:y] || σ == [:y,:y])

 @test unify([:x,:y],:(P(x,f(x))),:(P(a,f(y)))) == [:a,:a]
 @test unify([:x,:y],:(P(x,y)),:(P(y,f(a)))) == [:(f(a)),:(f(a))]

 @test unify([:x,:y], :(f(g(b,y),h(a))), :(f(x,y))) == [:(g(b,h(a))),:(h(a))]

 @testset "unify intervention" begin
  @test unify([:x,:y],:(P(f(x),x)),:(P(y,a))) == [:a,:(f(a))]

  @test unify([:x,:y,:w],:(P(f(w,x),w,x)),:(P(y,h(b),g(a)))) == [:(g(a)),:(f(h(b),g(a))),:(h(b))]

  @test unify([:x,:y,:w,:u,:z,:n],:(P(f(w,z),w,g(n),z)),:(P(y,h(b),u,h(u)))) == [:x,:(f(h(b),h(g(n)))),:(h(b)),:(g(n)),:(h(g(n))),:n]

  @test unify([:x,:y,:z,:w,:u,:n,:v],:((P(f(x),y,h(y,z),z,k(n)))),:((P(w,g(w),u,q(v),v)))) == [:x,:(g(f(x))),:(q(k(n))),:(f(x)),:(h(g(f(x)),q(k(n)))),:n,:(k(n))]

  @test unify([:x,:y,:z,:w,:u,:n,:v],:((P(f(y,z),y,h(z),z,k(n)))),:((P(w,g(u,v),u,q(v),v)))) == [:x,:(g(h(q(k(n))),k(n))),:(q(k(n))),:(f(g(h(q(k(n))),k(n)),q(k(n)))),:(h(q(k(n)))),:n,:(k(n))]

 end

 @test unify([:x,:y], :(f(a=x,b=x)), :(f(a=p,b=y))) == [:p, :p]

 @testset "unify loop test fail" begin
  @test_throws Loop unify([:x,:y],:(P(x)),:(P(f(x))))
  @test_throws Loop unify([:x,:y],:(P(x,y)),:(P(y,f(y))))
 end

end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
#ДАНО: Робот находится в произвольной клетке ограниченного прямоугольного поля без внутренних перегородок и маркеров.
#РЕЗУЛЬТАТ: Робот — в исходном положении в центре прямого креста из маркеров, расставленных вплоть до внешней рамки

function glavnaia(r::Robot) # - главная функция
    for side in (Nord,West,Sud,Ost) # - перебор всех возможных направлений
        putmarkers!(r,side)
        move_markers(r,inverse(side))
    end
    putmarker!(r)
end

#Cтавит маркеры до перегородки, но в исходной клетке маркера нет
putmarkers!(r::Robot,side::HorizonSide) =
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end

# Перемещает в заданном направлении, пока он в клетке с маркером (в итоге робот окажется в клетке без маркера)
move_markers(r::Robot,side::HorizonSide) =
    while ismarker(r)==true
        move!(r,side)
    end

# Возвращает направление, противоположное заданному
inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2, 4))

max_stars_repo_name: Aminulik/mirea-progs

max_issues_repo_name: Aminulik/mirea-progs

max_forks_repo_name: Aminulik/mirea-progs


text:
using IterTools

export get_subseqs, get_subseqs!, get_subseq_counts!, get_subseq_counts, get_subseq_proportions
export get_subpaths!, get_subpaths, get_subpath_counts!, get_subpath_counts!, get_subpath_proportions
export subseq_isin, subpath_isin
export dict_rank

function get_subseqs!(
    out::Set{Vector{T}},
    S::InteractionSequence{T},
    r::Int
    ) where {T<:Union{Int,String}}
    for path in S
        n = length(path)
        for ind in subsets(1:n, r)
            push!(out, path[ind])
        end
    end
end

function get_subseqs(S::InteractionSequence{T}, r::Int) where {T<:Union{Int,String}}
    out = Set{Vector{T}}()
    get_subseqs!(out, S, r)
    return out
end

function get_subseq_counts!(
    out::Dict{Vector{T},S},
    data::InteractionSequenceSample{T},
    r::Int
    ) where {T<:Union{Int,String},S<:Real}

    uniq_paths = Set{Vector{T}}()
    for x in data
        empty!(uniq_paths)
        get_subseqs!(uniq_paths, x, r)
        for path in uniq_paths
            out[path] = get(out, path, 0) + 1
        end
    end
end

function get_subseq_counts(data::InteractionSequenceSample{T}, r::Int) where {T<:Union{Int,String}}
    out = Dict{Vector{T},Int}()
    get_subseq_counts!(out, data, r)
    return out
end

function get_subseq_proportions(data::InteractionSequenceSample{T}, r::Int) where {T<:Union{Int,String}}
    out = Dict{Vector{T},Float64}()
    n = length(data)
    get_subseq_counts!(out, data, r)
    for key in keys(out)
        out[key] /= n
    end
    return out
end

function subseq_isin(S::InteractionSequence{T}, x::Vector{T}) where {T<:Union{Int,String}}
    r = length(x)
    for path in S
        n = length(path)
        for ind in subsets(1:n, r)
            if x == path[ind]
                return true
            end
        end
    end
    return false
end

function subseq_isin(data::InteractionSequenceSample{T}, x::Vector{T}) where {T<:Union{Int,String}}
    out = Vector{Bool}()
    for S in data
        push!(out, subseq_isin(S, x))
    end
    return out
end

# Subpaths
# --------

function get_subpaths!(
    out::Set{Vector{T}},
    S::InteractionSequence{T},
    r::Int
    ) where {T<:Union{Int,String}}
    for path in S
        ind = 1:r
        n = length(path)
        for i in 1:(n-r+1)
            push!(out, path[ind])
            ind = ind .+ 1
        end
    end
end

function get_subpaths(
    S::InteractionSequence{T},
    r::Int
    ) where {T<:Union{Int,String}}
    out = Set{Vector{T}}()
    get_subpaths!(out, S, r)
    return out
end

function get_subpath_counts!(
    out::Dict{Vector{T},S},
    data::InteractionSequenceSample{T},
    r::Int
    ) where {T<:Union{Int,String},S<:Real}

    uniq_paths = Set{Vector{T}}()
    for x in data
        empty!(uniq_paths)
        get_subpaths!(uniq_paths, x, r)
        for path in uniq_paths
            out[path] = get(out, path, 0) + 1
        end
    end
end

function get_subpath_counts(data::InteractionSequenceSample{T}, r::Int) where {T<:Union{Int,String}}
    out = Dict{Vector{T},Int}()
    get_subpath_counts!(out, data, r)
    return out
end

function get_subpath_proportions(data::InteractionSequenceSample{T}, r::Int) where {T<:Union{Int,String}}
    out = Dict{Vector{T},Float64}()
    n = length(data)
    get_subpath_counts!(out, data, r)
    for key in keys(out)
        out[key] /= n
    end
    return out
end

function subpath_isin(S::InteractionSequence{T}, x::Vector{T}) where {T<:Union{Int,String}}
    r = length(x)
    for path in S
        n = length(path)
        ind = 1:r
        for i in 1:(n-r+1)
            if path[ind] == x
                return true
            end
            ind = ind .+ 1
        end
    end
    return false
end

function subpath_isin(data::InteractionSequenceSample{T}, x::Vector{T}) where {T<:Union{Int,String}}
    out = Vector{Bool}()
    for S in data
        push!(out, subpath_isin(S, x))
    end
    return out
end

function dict_rank(d::Dict{T, S}) where {T<:Any, S<:Real}
    pairs = collect(d) # Makes vector of pairs
    vals = collect(values(d))
    ind = sortperm(vals, rev=true)
    return pairs[ind]
end


text:
# Example Configuration for nearly neutral infinite sites simuation
# Command line to run:  julia run.jl examples/is_example2
# Creates output file examples/in_example2.csv
function dfe_disadv( x::Int64 )
  dfe_deleterious( x::Int64; alpha=dfe_disadv_alpha, beta=dfe_disadv_beta )
end

const type_str = "D4"
const mu_list_flag = false
const nn_simtype = 2    # nearly neutral infinite sites
const N_list     = [100,200]
const N_mu_list= [1.0,2.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci # not used in infinite alleles
const ngens = 1000      # Generations after burn-in
const burn_in= 2.0    # generations of burn_in as a multiple of N
dfe = dfe_disadv
dfe_str = "dfe_disadv alpha_disadv=0.2 beta_disadv=1.0"
dfe_disadv_alpha=0.2
dfe_disadv_beta=1.0

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
function getTableauSPARKGLRK(s)
    g = getCoefficientsGLRK(s)

    ω = zeros(s, s+1)
    δ = zeros(0, s)

    for i in 1:s-1
        for j in 1:s
            ω[i,j] = g.b[j] * g.c[j]^(i-1)
        end
    end

    ω[s,s+1] = 1

    return TableauSPARK(Symbol("sparkglrk", s), g.o,
                        g.a, g.a, g.a, g.a,
                        g.a, g.a, g.a, g.a,
                        g.b, g.b, g.b, g.b,
                        g.c, g.c, g.c, g.c,
                        ω, δ)
end


text:
#    Copyright 2015 Raytheon BBN Technologies
#
#     Licensed under the Apache License, Version 2.0 (the "License");
#     you may not use this file except in compliance with the License.
#     You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#     Unless required by applicable law or agreed to in writing, software
#     distributed under the License is distributed on an "AS IS" BASIS,
#     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#     See the License for the specific language governing permissions and
#     limitations under the License.

__precompile__(false)
module QuantumTomography

import Distributions.fit

import Convex, Distributions, SCS, QuantumInfo
import SCS
import SparseArrays
import LinearAlgebra

using LinearAlgebra: Diagonal, norm, ishermitian, isdiag, tr
using Convex: Variable, minimize, maximize, square, norm, isposdef, solve!, tr

include("state-tomo.jl")

include("process-tomo.jl")

end # module


text:
function update_parameters!(net::Net{GPUBackend}, method::SGD, learning_rate, momentum, param_blob, hist_blob, gradient, data_type)
  # hist_blob = momentum * hist_blob
  CuBLAS.scal(net.backend.cublas_ctx, length(hist_blob), convert(data_type, momentum),
      hist_blob.ptr, 1)
  # hist_blob = -learning_rate * gradient + hist_blob
  CuBLAS.axpy(net.backend.cublas_ctx, length(hist_blob), convert(data_type, -learning_rate),
      gradient.ptr, 1, hist_blob.ptr, 1)

  # update parameter
  # param_blob = param_blob + hist_blob
  CuBLAS.axpy(net.backend.cublas_ctx, length(hist_blob), convert(data_type, 1),
      hist_blob.ptr, 1, param_blob.ptr, 1)
end


text:
using Test, LinearAlgebra, MPStates

@testset "Initialization of the Mps class" begin

"""Test whether a tensor product is left canonical."""
function is_left_canonical(A::Vector{Array{T, 3}}) where T
    output = true
    for i=1:length(A)-1
        Ai = A[i]
        rA = reshape(Ai, (size(Ai, 1)*size(Ai, 2), size(Ai, 3)))
        M = rA'*rA
        Id = Matrix{T}(I, size(Ai, 3), size(Ai, 3))
        output = isapprox(M, Id)
    end
    return output
end

"""Test whether a tensor product is right canonical."""
function is_right_canonical(A::Vector{Array{T, 3}}) where T
    output = true
    for i=length(A):-1:2
        Ai = A[i]
        rA = reshape(Ai, (size(Ai, 1), size(Ai, 2)*size(Ai, 3)))
        M = rA*rA'
        Id = Matrix{T}(I, size(Ai, 1), size(Ai, 1))
        output = isapprox(M, Id)
    end
    return output
end

@testset "make tensor state left & right canonical" begin
    # Test if state is initialized in left canonical form if it is normalized.
    rtest1 = MPStates.testMps("rtest1")
    @test is_left_canonical(rtest1.M)
    @test MPStates.norm(rtest1) ≈ 1.
    MPStates.make_right_canonical!(rtest1)
    @test is_right_canonical(rtest1.M)
    @test MPStates.norm(rtest1) ≈ 1.
end
end # @testset "Initialization of the Mps class"


text:
using AWSS3
using Documenter

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

makedocs(;
    modules=[AWSS3],
    sitename="AWSS3.jl",
    format=Documenter.HTML(;
        canonical="https://juliacloud.github.io/AWSS3.jl/stable/", edit_link="main"
    ),
    pages=["Home" => "index.md", "API" => "api.md"],
)

deploydocs(; repo="github.com/JuliaCloud/AWSS3.jl.git", push_preview=true)


text:
@testset "Constructors and basic properties" begin
    let F = FullBinner()
        @test typeof(F) <: AbstractVector{Float64}
        @test eltype(F) == Float64
        @test ndims(F) == 1
        @test length(F) == 0
        @test size(F) == (0,)
        @test lastindex(F) == 0
        @test axes(F) == (Base.OneTo(0),)
        @test isempty(F)
    end

    let F = FullBinner(ComplexF64)
        @test typeof(F) <: AbstractVector{ComplexF64}
        @test eltype(F) == ComplexF64
    end

    let x = [0.84, 0.381169, 0.34174, 0.888868, 0.0973183, 0.722725, 0.0957878, 0.432443, 0.755033, 0.864252]
        F = FullBinner(x)
        @test length(F) == 10
        @test lastindex(F) == 10
        @test axes(F) == (Base.OneTo(10),)
        @test !isempty(F)
    end

    let x = [rand(2,3) for _ in 1:100]
        F = FullBinner(x)
        @test length(F) == 100
        @test eltype(F) == Array{Float64,2}
    end
end

@testset "Scalars statistics" begin
    # Real
    let F = FullBinner(1:10_000)
        @test isapprox(std_error(F), 361.4079699881821)
        bs, stds, cum_stds = all_binning_errors(F)
        @test bs == 1:312
        @test isapprox(sum(stds), 106377.96306621947) # take sum as approx. hash
        @test isapprox(sum(cum_stds), 75541.44622415205)
        @test isapprox(tau(F), 77.86159630295694)

        # beta: convergence
        @test !BinningAnalysis.isconverged(F)
    end

    # Test 0/0 bug in R_value
    @test std_error(FullBinner(fill(1.0, 100))) == 0.0

    # Complex
    let F = FullBinner((1:10_000) .+ ((10_000:-1:1) .* im))
        @test isapprox(std_error(F), 511.10805270701564)
        bs, stds, cum_stds = all_binning_errors(F)
        @test bs == 1:312
        @test isapprox(sum(stds), 150441.1581058718) # take sum as approx. hash
        @test isapprox(sum(cum_stds), 106831.73777147368)
        @test isapprox(tau(F), 77.86159630295694)

        # beta: convergence
        @test !BinningAnalysis.isconverged(F)
    end

    # R -> tau conversion
    @test BinningAnalysis._tau(2.4) == 0.7
end

@testset "Arrays statistics" begin
    # Real
    Random.seed!(123)
    let F = FullBinner([rand(2,3) for _ in 1:100])
        @test length(F) == 100
        @test isapprox(std_error(F), [0.029184472105069394 0.029581605926346424 0.027793717502753976; 0.029105387394205307 0.02741415651581391 0.029933054433434834])
        @test isapprox(tau(F), [0.030183772076860294 0.027610002544459222 -0.047279410457739646; -0.003567404598109447 0.03282685243862249 0.06928989962602228])
    end

    # Complex
    Random.seed!(123)
    let F = FullBinner([rand(ComplexF64, 2,3) for _ in 1:100])
        @test length(F) == 100
        @test isapprox(std_error(F), [0.04445492633322362 0.04004496543964919 0.039737207226072296; 0.04099334255252945 0.039004215520294906 0.0409503504806149])
        @test isapprox(tau(F), [0.062160208644724047 -0.005039096051545122 0.037487751473977315; 0.03850424600452085 -0.032173894672710424 0.0320190059417359])
    end
end

@testset "Cosmetics (show, print, etc.)" begin
    F = FullBinner();
    # empty binner
    oldstdout = stdout
    (read_pipe, write_pipe) = redirect_stdout()
    println(F) # compact
    show(write_pipe, MIME"text/plain"(), F) # full
    redirect_stdout(oldstdout);
    close(write_pipe);

    # compact
    @test readline(read_pipe) == "FullBinner{Float64,Array{Float64,1}}()"
    # full
    @test readline(read_pipe) == "FullBinner{Float64,Array{Float64,1}}"
    @test readline(read_pipe) == "| Count: 0"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);

    # filled binner
    Random.seed!(1234)
    append!(F, rand(1000))
    (read_pipe, write_pipe) = redirect_stdout()
    show(write_pipe, MIME"text/plain"(), F)
    redirect_stdout(oldstdout);
    close(write_pipe);
    @test readline(read_pipe) == "FullBinner{Float64,Array{Float64,1}}"
    @test readline(read_pipe) == "| Count: 1000"
    @test readline(read_pipe) == "| Mean: 0.49685"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);
end


text:
##############################################################################
#
# REFERENCES: Forbes et al. "Statistical Distributions"
#
##############################################################################

immutable BetaPrime <: ContinuousUnivariateDistribution
    alpha::Float64
    beta::Float64
    function BetaPrime(a::Float64, b::Float64)
        (a > zero(a) && b > zero(b)) || error("Alpha and beta must be positive")
        new(float64(a), float64(b))
    end
end

BetaPrime() = BetaPrime(1.0, 1.0)

@distr_support BetaPrime 0.0 Inf

function mean(d::BetaPrime)
    d.beta > 1.0 ? d.alpha / (d.beta - 1.0) : NaN
end

mode(d::BetaPrime) = d.alpha > 1.0 ? (d.alpha - 1.0) / (d.beta + 1.0) : 0.0

function pdf(d::BetaPrime, x::Float64)
    α, β = d.alpha, d.beta
    (x^(α - 1.0) * (1.0 + x)^(-(α + β))) / beta(α, β)
end

cdf(d::BetaPrime, x::Float64) = cdf(Beta(d.alpha, d.beta), x / (one(x) + x))
ccdf(d::BetaPrime, x::Float64) = ccdf(Beta(d.alpha, d.beta), x / (one(x) + x))
logcdf(d::BetaPrime, x::Float64) = logcdf(Beta(d.alpha, d.beta), x / (one(x) + x))
logccdf(d::BetaPrime, x::Float64) = logccdf(Beta(d.alpha, d.beta), x / (one(x) + x))

quantile(d::BetaPrime, p::Float64) = (x = quantile(Beta(d.alpha,d.beta),p); x / (1.0-x))
cquantile(d::BetaPrime, p::Float64) = (x = cquantile(Beta(d.alpha,d.beta),p); x / (1.0-x))
invlogcdf(d::BetaPrime, p::Float64) = (x = invlogcdf(Beta(d.alpha,d.beta),p); x / (1.0-x))
invlogccdf(d::BetaPrime, p::Float64) = (x = invlogccdf(Beta(d.alpha,d.beta),p); x / (1.0-x))

function rand(d::BetaPrime)
    x = rand(Gamma(d.alpha, 1))
    y = rand(Gamma(d.beta, 1))
    x/y
end

function skewness(d::BetaPrime)
    α, β = d.alpha, d.beta
    β > 3.0 ? (2.0 * (2.0 * α + β - 1))/(β - 3.0) * sqrt((β - 2.0)/(α * (α + β - 1.0))) : NaN
end

function var(d::BetaPrime)
    α, β = d.alpha, d.beta
    β > 2.0 ? (α * (α + β - 1.0)) / ((β - 2.0) * (β - 1.0)^2) : NaN
end


text:
using NgSpice
# using Plots
n = NgSpice

netpath = joinpath(@__DIR__, "..", "inputs", "mosfet.cir") |> normpath

n.init()
n.source(netpath)
n.run()
n.display()        # Returns a dataframe of current vectors
n.listallvecs()    # Lists vectors of both active and inactive plots
n.curplot()        # Current active plot
n.listcurvecs()    # Vectors in current active plot

# For plotting the signals, pass `n.graph` an object of `NgSpiceGraphs` type.
# And then, specify the get-vector method.
# Pass "time" or "frequency"
# List of signals to plot
# Finally, the title of the plot
# plot(n.graph, n.getimagvec, "time", ["emit", "vcc"], "Simple MOSFET")
# plot(n.graph, n.getrealvec, "time", ["emit", "vcc"], "Simple MOSFET")

max_stars_repo_name: JuliaComputing/Ngspice

max_issues_repo_name: JuliaComputing/Ngspice

max_forks_repo_name: JuliaComputing/Ngspice


text:
#=
    Contents: Example: Display last visit with LoggedDict as database.
=#
using HttpServer
using AccessControl

# LoggedDict as data store for sessions
using LoggedDicts
sessions = LoggedDict("sessions", "sessions.log", true)    # Logging turned off
AccessControl.update_config!(session = Dict(:datastore => sessions))

# Handler
function home!(req, res)
    session_id = read_sessionid(req)
    if session_id == ""                             # "id" cookie does not exist...session hasn't started...start a new session.
        session_id = session_create!(res, "")
        res.data   = "This is your first visit."
        session_set!(session_id, "lastvisit", string(now()))
    else
        last_visit = session_get(session_id, "lastvisit")
        res.data   = "Welcome back. Your last visit was at $last_visit."
        session_set!(session_id, "lastvisit", string(now()))
    end
end

# App
function app(req::Request)
    res = Response()
    if req.resource == "/home"
        home!(req, res)
    else
        notfound!(res)
    end
    res
end

### Run the app under HTTPS
# Generate keys/server.key and keys/server.crt if they don't already exist
rel(filename::AbstractString, p::AbstractString) = joinpath(dirname(filename), p)
if !isfile("keys/server.crt")
    @unix_only begin
	run(`mkdir -p $(rel(@__FILE__, "keys"))`)
	run(`openssl req -x509 -nodes -days 365 -newkey rsa:2048 -keyout $(rel(@__FILE__, "keys/server.key")) -out $(rel(@__FILE__, "keys/server.crt"))`)
    end
end

# Define and run the server
server = Server((req, res) -> app(req))
cert   = MbedTLS.crt_parse_file(rel(@__FILE__, "keys/server.crt"))
key    = MbedTLS.parse_keyfile(rel(@__FILE__, "keys/server.key"))
run(server, port = 8000, ssl = (cert, key))


text:
"""
    inpolygon(points, poly)

Check if a point or points are inside a polygon.

This algorithm is very efficient for many points, less so a single point.

# Arguments

- `points`: an `AbstractVector` or a `Tuple` or `Real`, an `AbstractVector` of these,
    or a `GeoInterface.AbstractGeometry`.
- `poly`: an `AbstractVector` or nested `AbstractVector` with an inner
    `AbstractVector` or `Tuple` of `Real`, or a `GeoInterface.AbstractPolygon`.

Returns a `Bool` or `AbstractVector{Bool}`.
"""
function inpolygon(point::Union{Tuple,AbstractVector{<:Real}}, poly; kw...)
    inpolygon([point], poly; kw...)[1]
end
function inpolygon(points::GI.AbstractGeometry, poly; kw...)
    inpolygon(collect(_flat_nodes(GI.coordinates(points))), poly, kw...)
end
function inpolygon(points, poly; kw...)
    vecs = _node_vecs(poly)
    edges, nodes = to_edges_and_nodes(vecs)
    inpoly_matrix = inpoly2(points, nodes, edges; kw...)
    return view(inpoly_matrix, :, 1)
end


text:
#------------------------------------------------------------------------------------------------------------
#   Define functions needed to solve the covariance estimation problem:
#           min_{X \in \R^{d x r}}   F(X)   :=    E_{a,b} | <XX^T, aa^T> - b |
#   and it's elastic net variant with
#           F(X) :=  E_{a,b} | <XX^T, aa^T> - b | + (<XX^T, aa^T> - b )^2
#
#   Solve via two methods:
#           1) stochastic subgradient (Φ = 1/2 ||⋅||^2)
#           2) stochastic mirror descent (Φ = a polynomial)
#
#   Note that stochastic subgradients are
#           ∂ (  | <XX^T, a a^T> - b |  )   =    sign(res) * 2 aa^T X
#           ∂ (  | <XX^T, a a^T> - b |  )   =    (2 * res + sign(res) ) * 2 aa^T X
#   where (a, b) ∈ Rᵐ x R are chosen randomly, and
#           res = <XX^T, a a^T> - b = || aTX ||₂² - b  (residual)
#------------------------------------------------------------------------------------------------------------

using LinearAlgebra
using Polynomials
using Random

#------------------------------------------------------------------------------------------
# Draw stochastic a's (gaussian) and corresponding b's
# Returns: maxIter of them, in format:
#                       A = matrix with a's in columns
#                       B = vector with b's in entries
#------------------------------------------------------------------------------------------
function get_ab(XTtrue, stochErr::Float64, maxIter::Int64; noiseType::String="gaussian")
    (r,d) = size(XTtrue);
    A = randn(d, maxIter);  # gaussian entries in a's
    bError = stochErr * randn(1, maxIter);  # normal(0,stochErr²) errors
    if noiseType=="sparse"
        p = 0.25; # probability of corrupting a measurement
        noiseless_idxs = randsubseq(collect(1:maxIter), 1-p);
        bError[noiseless_idxs] = zeros(length(noiseless_idxs));
    end
    B = sum(abs2, XTtrue*A, dims=1) + bError;
    return(A,B)
end

#---------------------------------------------------------------------------------------------------
#   Compute empirical function values, based on (A,B) as above
#   Inputs:     X = point to evaluate function (d x r matrix)
#                   A = matrix with columns as the a's
#                   B = vector with b's as entries
#---------------------------------------------------------------------------------------------------
# vanilla covariance estimation function value
function compute_empirical_function(X, A, B)
    residuals = sum(abs2, X'*A, dims=1) - B;
    return sum(abs, residuals)
end

# elastic net empirical function value
function compute_enet_empirical_function(X, A, B)
    residuals = sum(abs2, X'*A, dims=1) - B;
    return (sum(abs, residuals) + sum(abs2, residuals)) / size(A,2)
end

#----------------------------------------------------------------------------
# Update the subgradient G, in place
#----------------------------------------------------------------------------
function subgrad!(G, XTa, a, b)
    res = sum(abs2, XTa) - b;
    lmul!(0.0,G)  # reset G to zero
    BLAS.ger!(2*sign(res), a, XTa, G);  # rank one update
end

function enet_subgrad!(G, XTa, a, b)
    res = sum(abs2, XTa) - b;
    lmul!(0.0,G)  # reset G to zero
    BLAS.ger!(2*sign.(res)+2*res, a, XTa, G);  # rank one update
end

#------------------------------------------------------------------------------------------------------------
# Find the smallest real, positive root of a polynomial
# input : coefficients, as a vector by increasing degree (e.g. [a0, a1, a2, a3])
#------------------------------------------------------------------------------------------------------------
function get_root(coeffs)
    all_roots = roots(Poly(coeffs));

    # reduce to positive real roots
    pos_real_idxs = map(x->(imag(x)==0.0 && real(x)>0), all_roots);
    pos_real_roots = real(all_roots[pos_real_idxs]);

    # return smallest positive real root
    return(minimum(pos_real_roots))
end

max_stars_repo_name: kelliemac/stochastic-model-based

max_issues_repo_name: kelliemac/stochastic-model-based

max_forks_repo_name: kelliemac/stochastic-model-based


text:
module Sobol
using Random
export SobolSeq, ScaledSobolSeq, next!

include("soboldata.jl") #loads `sobol_a` and `sobol_minit`

abstract type AbstractSobolSeq{N} end

# N iis the dimension of sequence being generated
mutable struct SobolSeq{N} <: AbstractSobolSeq{N}
    m::Array{UInt32,2} #array of size (sdim, 32)
    x::Array{UInt32,1} #previous x = x_n, array of length sdim
    b::Array{UInt32,1} #position of fixed point in x[i] is after bit b[i]
    n::UInt32 #number of x's generated so far
end

ndims(s::AbstractSobolSeq{N}) where {N} = N::Int

function SobolSeq(N::Int)
    (N < 0 || N > (length(sobol_a) + 1)) && error("invalid Sobol dimension")

    m = ones(UInt32, (N, 32))

    #special cases
    N == 0 && return(SobolSeq{0})(m,UInt32[],UInt32[],zero(UInt32))
    #special cases 1
    N == 1 && return(SobolSeq{N}(m,UInt32[0],UInt32[0],zero(UInt32)))

    for i = 2:N
        a = sobol_a[i-1]
        d = floor(Int, log2(a)) #degree of poly

        #set initial values of m from table
        m[i, 1:d] = sobol_minit[1:d, i - 1]
        #fill in remaining values using recurrence
        for j = (d+1):32
            ac = a
            m[i,j] = m[i,j-d]
            for k = 0:d-1
                @inbounds m[i,j] = m[i,j] ⊻ (((ac & one(UInt32)) * m[i, j-d+k]) << (d-k))
                ac >>= 1
            end
        end
    end
    SobolSeq{N}(m,zeros(UInt32,N),zeros(UInt32,N),zero(UInt32))
end
SobolSeq(N::Integer) = SobolSeq(Int(N))

function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat})
    length(x) != ndims(s) && throw(BoundsError())

    if s.n == typemax(s.n)
        return rand!(x)
    end

    s.n += one(s.n)
    c = UInt32(trailing_zeros(s.n))
    sb = s.b
    sx = s.x
    sm = s.m
    for i=1:ndims(s)
        @inbounds b = sb[i]
        # note: ldexp on Float64(sx[i]) is exact, independent of precision of x[i]
        if b >= c
            @inbounds sx[i] = sx[i] ⊻ (sm[i,c+1] << (b-c))
            @inbounds x[i] = ldexp(Float64(sx[i]), ((~b) % Int32))
        else
            @inbounds sx[i] = (sx[i] << (c-b)) ⊻ sm[i,c+1]
            @inbounds sb[i] = c
            @inbounds x[i] = ldexp(Float64(sx[i]), ((~c) % Int32))
        end
    end
    return x
end
next!(s::SobolSeq) = next!(s, Array{Float64,1}(undef, ndims(s)))

# if we know in advance how many points (n) we want to compute, then
# adopt the suggestion of the Joe and Kuo paper, which in turn
# is taken from Acworth et al (1998), of skipping a number of
# points equal to the largest power of 2 smaller than n
# if exactly n points are to be skipped, use the keyword exact=true.
#
# skip!(s, n) skips m such that 2^m < n < 2^(m+1)
# skip!(s, n, exact=true) skips m = n

function skip!(s::SobolSeq, n::Integer, x; exact=false)
    nskip = exact ? n : 1 << floor(Int,log2(n))
    for unused=1:nskip; next!(s,x); end
    return nothing
end
Base.skip(s::SobolSeq, n::Integer; exact=false) = skip!(s, n, Array{Float64,1}(undef, ndims(s)); exact=exact)

function Base.show(io::IO, s::SobolSeq)
    print(io, "$(ndims(s))-dimensional Sobol sequence on [0,1]^$(ndims(s))")
end
function Base.show(io::IO, ::MIME"text/html", s::SobolSeq)
    print(io, "$(ndims(s))-dimensional Sobol sequence on [0,1]$(ndims(s))")
end

# Make an iterator so that we can do "for x in SobolSeq(...)".
# Technically, the Sobol sequence ends after 2^32-1 points, but it
# falls back on pseudorandom numbers after this.  In practice, one is
# unlikely to reach that point.
Base.iterate(s::AbstractSobolSeq, state=nothing) = (next!(s), state)
Base.eltype(::Type{<:AbstractSobolSeq}) = Vector{Float64}
Base.IteratorSize(::Type{<:AbstractSobolSeq}) = Base.IsInfinite()
Base.IteratorEltype(::Type{<:AbstractSobolSeq}) = Base.HasEltype()

# Convenience wrapper for scaled Sobol sequences

struct ScaledSobolSeq{N,T} <: AbstractSobolSeq{N}
    s::SobolSeq{N}
    lb::Vector{T}
    ub::Vector{T}
    function ScaledSobolSeq{N,T}(lb::Vector{T}, ub::Vector{T}) where {N,T}
        length(lb)==length(ub)==N || throw(DimensionMismatch("lb and ub do not have length $N"))
        new(SobolSeq(N), lb, ub)
    end
end
function SobolSeq(N::Integer, lb, ub)
    T = typeof(sum(ub) - sum(lb))
    ScaledSobolSeq{Int(N),T}(copyto!(Vector{T}(undef,N), lb), copyto!(Vector{T}(undef,N), ub))
end
SobolSeq(lb, ub) = SobolSeq(length(lb), lb, ub)

function next!(s::SobolSeq, x::AbstractVector{<:AbstractFloat},
               lb::AbstractVector, ub::AbstractVector)
    length(x) < ndims(s) && throw(BoundsError())
    next!(s,x)
    for i=1:ndims(s)
        x[i] = lb[i] + (ub[i]-lb[i]) * x[i]
    end
    return x
end
function next!(s::SobolSeq{N}, lb::AbstractVector, ub::AbstractVector) where {N}
    T = typeof(float((zero(eltype(ub)) - zero(eltype(lb)))))
    next!(s, Vector{T}(undef, N), lb, ub)
end

next!(s::ScaledSobolSeq, x::AbstractVector{<:AbstractFloat}) = next!(s.s, x, s.lb, s.ub)
next!(s::ScaledSobolSeq{N,T}) where {N,T} = next!(s.s, Vector{float(T)}(undef, N), s.lb, s.ub)
Base.eltype(::Type{ScaledSobolSeq{N,T}}) where {N,T} = Vector{float(T)}

Base.skip(s::ScaledSobolSeq, n) = skip(s.s, n)

function Base.show(io::IO, s::ScaledSobolSeq{N,T}) where {N,T}
    lb = s.lb; ub = s.ub
    print(io, "$N-dimensional scaled $(float(T)) Sobol sequence on [$(lb[1]),$(ub[1])]")
    cnt = 1
    for i = 2:N
        if lb[i] == lb[i-1] && ub[i] == ub[i-1]
            cnt += 1
        else
            if cnt > 1
                print(io, "^", cnt)
            end
            print(io, " x [", lb[i], ",", ub[i], "]")
            cnt = 1
        end
    end
    if cnt > 1
        print(io, "^", cnt)
    end
end

function Base.show(io::IO, ::MIME"text/html", s::ScaledSobolSeq{N,T}) where {N,T}
    lb = s.lb; ub = s.ub
    print(io, "$N-dimensional scaled $(float(T)) Sobol sequence on [$(lb[1]),$(ub[1])]")
    cnt = 1
    for i = 2:N
        if lb[i] == lb[i-1] && ub[i] == ub[i-1]
            cnt += 1
        else
            if cnt > 1
                print(io, "", cnt, "")
            end
            print(io, " × [", lb[i], ",", ub[i], "]")
            cnt = 1
        end
    end
    if cnt > 1
        print(io, "", cnt, "")
    end
end

end # module


text:
using SgdlDev
using Documenter

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

makedocs(;
    modules=[SgdlDev],
    authors="Laurent Daniel <laurent.daniel@sgdl.com> and contributors",
    repo="https://github.com/laurent-daniel/SgdlDev.jl/blob/{commit}{path}#{line}",
    sitename="SgdlDev.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://laurent-daniel.github.io/SgdlDev.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/laurent-daniel/SgdlDev.jl",
    devbranch="master",
)


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/11/2014

#Piggyback and add coordination to SimpleTCAS

module CoordSimpleTCASImpl

export
    addObserver,

    initialize,
    update,

    testThreat,
    selectRA,

    CoordSimpleTCAS,
    CoordSimpleTCASInput,
    CoordSimpleTCASResolutionAdvisory

using AbstractCollisionAvoidanceSystemImpl
using AbstractCollisionAvoidanceSystemInterfaces
using CommonInterfaces

using Util
using RLESUtils, Observers

using AbstractCASCoordImpl

using CASCoordination
using SimpleTCASImpl

import CommonInterfaces.addObserver
import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractCollisionAvoidanceSystemInterfaces.testThreat
import AbstractCollisionAvoidanceSystemInterfaces.selectRA

type CoordSimpleTCAS <: AbstractCollisionAvoidanceSystem
  aircraft_number::Int64
  simpleTCAS::SimpleTCAS
  coordination::AbstractCASCoord

  CoordSimpleTCAS(aircraft_number::Int,coordination::AbstractCASCoord) = new(aircraft_number,SimpleTCAS(),coordination)
end

addObserver(cas::CoordSimpleTCAS, f::Function) = add_observer(cas.observer, f)
addObserver(cas::CoordSimpleTCAS, tag::AbstractString, f::Function) = add_observer(cas.observer, tag, f)

testThreat(cas::CoordSimpleTCAS, input::SimpleTCASInput) = testThreat(cas.simpleTCAS, input)

function selectRA(cas::CoordSimpleTCAS, input::SimpleTCASInput)

  #Make a two aircraft assumption for now
  if cas.aircraft_number == 1
    intruder_number = 2
  else
    intruder_number = 1
  end

  intruderRA = getRecord(cas.coordination,intruder_number) #The coordination object is ownship's RA

  if intruderRA != nothing
    h_d = intruderRA.h_d < 0 ? 1500 / 60 : -1500 / 60  #Do opposite sense, default strength
  else
    h_d = selectRA(cas.simpleTCAS, input)
  end

  return h_d
end

update(cas::CoordSimpleTCAS, input) = update(cas,convert(SimpleTCASInput, input))

function update(cas::CoordSimpleTCAS, input::SimpleTCASInput)

  if cas.simpleTCAS.b_TCAS_activated == false
    cas.simpleTCAS.b_TCAS_activated = testThreat(cas, input)

    if cas.simpleTCAS.b_TCAS_activated
      cas.simpleTCAS.RA = SimpleTCASResolutionAdvisory(selectRA(cas, input))
      setRecord(cas.coordination,cas.aircraft_number,cas.simpleTCAS.RA)
    end
  end

  return cas.simpleTCAS.RA
end

function initialize(cas::CoordSimpleTCAS)
    initialize(cas.simpleTCAS)
end

end


text:
struct PipelineDynamicStateCreateInfo
    mHandleRef::Ref{vk.VkPipelineDynamicStateCreateInfo}
    mPreserve::Vector{Any}

    function PipelineDynamicStateCreateInfo(;
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineDynamicStateCreateFlags = vk.VkFlags(0),
        states::Vector{vk.VkDynamicState} = Vector{vk.VkDynamicState}()
    )

        this = new(Ref(vk.VkPipelineDynamicStateCreateInfo(
            vk.VK_STRUCTURE_TYPE_PIPELINE_DYNAMIC_STATE_CREATE_INFO, #sType::VkStructureType
            pNext, #::Ptr{Cvoid}
            flags, #::VkPipelineDynamicStateCreateFlags
            length(states), #dynamicStateCount::UInt32
            pointer(states) #pDynamicStates::Ptr{VkDynamicState}
        )), [states])
    end
end

function handleRef(this::PipelineDynamicStateCreateInfo)::Ref{vk.VkPipelineDynamicStateCreateInfo}
    this.mHandleRef
end


text:
## Adapted from sundialsTB/kinsol/examples_ser/mkinTest_nds.m

## %mkinTest_dns - KINSOL example problem (serial, dense)
## %   Simple test problem for the Dense linear solver in KINSOL
## %   This example solves the system
## %       y(1)^2 + y(2)^2 = 1
## %       y(2) = y(1)^2
## %

## % Radu Serban <radu@llnl.gov>

using Sundials, Base.Test

## function to be optimized
function sysfn(y_nv, fy_nv, a_in)
    y = convert(Vector, y_nv)
    fy = convert(Vector, fy_nv)
    fy[1] = y[1]^2 + y[2]^2 - 1.0
    fy[2] = y[2] - y[1]^2
    return Sundials.KIN_SUCCESS
end

## Initialize problem
neq = 2
kmem = Sundials.KINCreate()
Sundials.@checkflag Sundials.KINSetFuncNormTol(kmem, 1.0e-5)
Sundials.@checkflag Sundials.KINSetScaledStepTol(kmem, 1.0e-4)
Sundials.@checkflag Sundials.KINSetMaxSetupCalls(kmem, 1)
y = ones(neq)
Sundials.@checkflag Sundials.KINInit(kmem, sysfn, y)
A = Sundials.SUNDenseMatrix(length(y),length(y))
LS = Sundials.SUNDenseLinearSolver(y,A)
Sundials.@checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, A)
## Solve problem
scale = ones(neq)
Sundials.@checkflag Sundials.KINSol(kmem, y, Sundials.KIN_LINESEARCH,
                                    scale, scale)

println("Solution: ", y)
residual = ones(2)
sysfn(y, residual, [1, 2])
println("Residual: ", residual)

@test abs(minimum(residual)) < 1e-5


text:
parse_kw(ps::ParseState, ::Type{Val{Tokens.IMPORT}}) = parse_imports(ps)
parse_kw(ps::ParseState, ::Type{Val{Tokens.IMPORTALL}}) = parse_imports(ps)
parse_kw(ps::ParseState, ::Type{Val{Tokens.USING}}) = parse_imports(ps)

parse_kw(ps::ParseState, ::Type{Val{Tokens.MODULE}}) = parse_module(ps)
parse_kw(ps::ParseState, ::Type{Val{Tokens.BAREMODULE}}) = parse_module(ps)

function parse_module(ps::ParseState)
    startbyte = ps.t.startbyte

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    if ps.nt.kind == Tokens.IDENTIFIER
        next(ps)
        arg = INSTANCE(ps)
    else
        @catcherror ps startbyte arg = @precedence ps 15 @closer ps block @closer ps ws parse_expression(ps)
    end

    block = EXPR{Block}(EXPR[], -ps.nt.startbyte, Variable[], "")
    @scope ps Scope{Tokens.MODULE} @default ps while ps.nt.kind !== Tokens.END
        @catcherror ps startbyte a = @closer ps block parse_doc(ps)
        push!(block.args, a)
    end

    # Construction
    block.span += ps.nt.startbyte
    next(ps)
    ret = EXPR{(kw isa EXPR{KEYWORD{Tokens.MODULE}} ? ModuleH : BareModule)}(EXPR[kw, arg, block, INSTANCE(ps)], ps.nt.startbyte - startbyte, Variable[], "")
    ret.defs = [Variable(Expr(arg), :module, ret)]
    return ret
end

function parse_dot_mod(ps::ParseState, colon = false)
    startbyte = ps.nt.startbyte
    args = EXPR[]

    while ps.nt.kind == Tokens.DOT || ps.nt.kind == Tokens.DDOT || ps.nt.kind == Tokens.DDDOT
        next(ps)
        d = INSTANCE(ps)
        if d isa EXPR{OPERATOR{DotOp,Tokens.DOT,false}}
            push!(args, EXPR{OPERATOR{DotOp,Tokens.DOT,false}}(EXPR[], 1, Variable[], ""))
        elseif d isa EXPR{OPERATOR{ColonOp,Tokens.DDOT,false}}
            push!(args, EXPR{OPERATOR{DotOp,Tokens.DOT,false}}(EXPR[], 1, Variable[], ""))
            push!(args, EXPR{OPERATOR{DotOp,Tokens.DOT,false}}(EXPR[], 1, Variable[], ""))
        elseif d isa EXPR{OPERATOR{DddotOp,Tokens.DDDOT,false}}
            push!(args, EXPR{OPERATOR{DotOp,Tokens.DOT,false}}(EXPR[], 1, Variable[], ""))
            push!(args, EXPR{OPERATOR{DotOp,Tokens.DOT,false}}(EXPR[], 1, Variable[], ""))
            push!(args, EXPR{OPERATOR{DotOp,Tokens.DOT,false}}(EXPR[], 1, Variable[], ""))
        end
    end

    # import/export ..
    if ps.nt.kind == Tokens.COMMA || ps.ws.kind == NewLineWS || ps.nt.kind == Tokens.ENDMARKER
        if length(args) == 2
            return EXPR[INSTANCE(ps)]
        end
    end

    while true
        if ps.nt.kind == Tokens.AT_SIGN
            next(ps)
            next(ps)
            a = INSTANCE(ps)
            a = EXPR{IDENTIFIER}(EXPR[], a.span + 1, Variable[], string("@", a.val))
            push!(args, a)
        elseif ps.nt.kind == Tokens.LPAREN
            next(ps)
            a = EXPR{InvisBrackets}(EXPR[INSTANCE(ps)], -ps.t.startbyte, Variable[], "")
            @catcherror ps startbyte push!(a.args, @default ps @closer ps paren parse_expression(ps))
            next(ps)
            push!(a.args, INSTANCE(ps))
            a.span += ps.nt.startbyte
            push!(args, a)
        elseif ps.nt.kind == Tokens.EX_OR
            @catcherror ps startbyte a = @closer ps comma parse_expression(ps)
            push!(args, a)
        elseif !colon && isoperator(ps.nt) && ps.ndot
            next(ps)
            push!(args, EXPR{OPERATOR{precedence(ps.t),ps.t.kind,false}}(EXPR[], ps.nt.startbyte - ps.t.startbyte - 1, Variable[], ""))
        else
            next(ps)
            push!(args, INSTANCE(ps))
        end

        if ps.nt.kind == Tokens.DOT
            next(ps)
            push!(args, INSTANCE(ps))
        elseif isoperator(ps.nt) && ps.ndot
            push!(args, EXPR{PUNCTUATION{Tokens.DOT}}(EXPR[], 1, Variable[], ""))
        else
            break
        end
        # if ps.nt.kind != Tokens.DOT
        #     break
        # else
        #     next(ps)
        #     push!(puncs, INSTANCE(ps))
        # end
    end
    args
end

function parse_imports(ps::ParseState)
    startbyte = ps.t.startbyte
    kw = INSTANCE(ps)
    kwt = kw isa EXPR{KEYWORD{Tokens.IMPORT}} ? Import :
          kw isa EXPR{KEYWORD{Tokens.IMPORTALL}} ? ImportAll :
          Using
    format_kw(ps)
    tk = ps.t.kind

    arg = parse_dot_mod(ps)

    if ps.nt.kind != Tokens.COMMA && ps.nt.kind != Tokens.COLON
        ret = EXPR{kwt}(EXPR[kw; arg], ps.nt.startbyte - startbyte, Variable[], "")
        ret.defs = [Variable(Expr(ret), :IMPORTS, ret)]
    elseif ps.nt.kind == Tokens.COLON

        ret = EXPR{kwt}(EXPR[kw;arg], 0, Variable[], "")
        t = 0

        next(ps)
        push!(ret.args, INSTANCE(ps))

        @catcherror ps startbyte arg = parse_dot_mod(ps, true)
        append!(ret.args, arg)
        while ps.nt.kind == Tokens.COMMA
            next(ps)
            push!(ret.args, INSTANCE(ps))
            @catcherror ps startbyte arg = parse_dot_mod(ps, true)
            append!(ret.args, arg)
        end
        if Expr(ret).head == :toplevel
            ret.defs = [Variable(d, :IMPORTS, ret) for d in Expr(ret).args]
        else
            ret.defs = [Variable(Expr(ret), :IMPORTS, ret)]
        end
    else
        ret = EXPR{kwt}(EXPR[kw;arg], 0, Variable[], "")
        while ps.nt.kind == Tokens.COMMA
            next(ps)
            push!(ret.args, INSTANCE(ps))
            @catcherror ps startbyte arg = parse_dot_mod(ps)
            append!(ret.args, arg)
        end
        ret.defs = [Variable(d, :IMPORTS, ret) for d in Expr(ret).args]
    end

    ret.span = ps.nt.startbyte - startbyte

    return ret
end

function parse_kw(ps::ParseState, ::Type{Val{Tokens.EXPORT}})
    startbyte = ps.t.startbyte

    # Parsing
    kw = INSTANCE(ps)
    format_kw(ps)
    ret = EXPR{Export}(EXPR[kw; parse_dot_mod(ps)], 0, Variable[], "")

    while ps.nt.kind == Tokens.COMMA
        next(ps)
        push!(ret.args, INSTANCE(ps))
        @catcherror ps startbyte arg = parse_dot_mod(ps)[1]
        push!(ret.args, arg)
    end
    ret.span = ps.nt.startbyte - startbyte

    return ret
end


text:
using Base.Broadcast: broadcastable
using BenchmarkTools
using ChainRulesCore
using LinearAlgebra: Diagonal
using Test

@testset "ChainRulesCore" begin
    @testset "differentials" begin
        include("differentials/abstract_zero.jl")
        include("differentials/one.jl")
        include("differentials/thunks.jl")
        include("differentials/composite.jl")
    end

    include("rules.jl")
end


text:
function extract_info end

# Extract info from FastqInfo
function extract_info(readsinfo::Array{Union{FastqInfo,Missing},1})
	quals = map(e -> round(e.quality, digits=1), readsinfo)
	lengths = map(e -> e.length, readsinfo)
	gcs = map(e -> e.gc, readsinfo)
	DataFrame(quality = quals, length = lengths, gc_content = gcs)
end

# Extract info from BAMInfo
function extract_info(readsinfo::Array{BAMInfo,1})
	quals = map(e -> round(e.quality, digits=1), readsinfo)
	lengths = map(e -> e.length, readsinfo)
	identities = map(e -> round(e.identity, digits=1), readsinfo)
	gcs = map(e -> e.gc, readsinfo)
	DataFrame(quality = quals, length = lengths, identity = identities, gc_content = gcs)
end


text:
#=

# TODO: refactor
function getbatch(dataloader::DataLoader, batchsize = dataloader.batchsize; offset = 0)
    idxs = 1 + offset:offset + batchsize
    samples = [DataLoaders.getsample(dataloader, idx) for idx in idxs]
    batch = dataloader.collatefn(samples)

    return batch
end
getbatch(databunch::DataBunch, train = true; kwargs...) = getbatch(train ? databunch.traindl : databunch.valdl; kwargs...)
getbatch(learner::Learner, train = true; kwargs...) = getbatch(learner.databunch, train; kwargs...)

function getoutputs(model, batch)
    x, y = gpu(batch)
    y_pred = gpu(model)(x)

    return y_pred, y
end

getoutputs(learner::Learner, batch = getbatch(learner)) = getoutputs(learner.model, batch)

function getloss(learner, batch = getbatch(learner))
    y_pred, y = getoutputs(learner, batch)
    learner.lossfn(y_pred, y)
end

=#


text:
# Stricker Klaus 2020 - License is MIT: http://julialang.org/license
# ReadWriteDlm2 - rwd2_dateregex - https://github.com/strickek/ReadWriteDlm2.jl

"""

    dfregex(df::AbstractString, locale::AbstractString=\"english\")

Create a regex string `r\"^...\$\"` for the given `Date` or `DateTime`
`format`string `df`.

The regex groups are named according to the `format`string codes. `locale` is
used to calculate min and max length of month and day names (for codes: UuEe).
"""
function dfregex(df::AbstractString, locale::AbstractString="english")
    # calculate min and max string lengths of months and day_of_weeks names
    Ule = try extrema([length(Dates.monthname(i;locale=locale)) for i in 1:12])catch; (3, 9) end
    ule = try extrema([length(Dates.monthabbr(i;locale=locale)) for i in 1:12])catch; (3, 3) end
    Ele = try extrema([length(Dates.dayname(i;locale=locale)) for i in 1:7])catch; (6, 9) end
    ele = try extrema([length(Dates.dayabbr(i;locale=locale)) for i in 1:7])catch; (3, 3) end

    codechars = 'y', 'Y', 'm', 'u', 'e', 'U', 'E', 'd', 'H', 'M', 'S', 's', 'Z', 'z', '\\'
    r = "^ *"; repeat_count = 1; ldf = length(df); dotsec = false
    for i = 1:ldf
        repeat_next = ((i < ldf) && (df[(i + 1)] == df[i])) ? true : false
        ((df[i] == '.') && (i < ldf) && (df[(i + 1)] == 's')) && (dotsec = true)
        repeat_count = (((i > 2) && (df[(i - 2)] != '\\') && (df[(i - 1)] == df[i])) ||
                        ((i == 2) && (df[1] == df[2]))) ? (repeat_count + 1) : 1
        r = r * (
        ((i > 1) && (df[(i - 1)] == '\\')) ? string(df[i]) :
        ((df[i] == 'y') && (repeat_count < 5) && !repeat_next) ? "(?<y>\\d{1,4})" :
        ((df[i] == 'y') && (repeat_count > 4) && !repeat_next) ? "(?<y>\\d{1,$repeat_count})" :
        ((df[i] == 'Y') && (repeat_count < 5) && !repeat_next) ? "(?<y>\\d{1,4})" :
        ((df[i] == 'Y') && (repeat_count > 4) && !repeat_next) ? "(?<y>\\d{1,$repeat_count})" :
        ((df[i] == 'm') && (repeat_count == 1) && !repeat_next) ? "(?<m>0?[1-9]|1[012])" :
        ((df[i] == 'm') && (repeat_count == 2) && !repeat_next) ? "(?<m>0[1-9]|1[012])" :
        ((df[i] == 'm') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<m>0[1-9]|1[012])" :
        ((df[i] == 'u') && (repeat_count == 1)) ? "(?<u>[A-Za-z\u00C0-\u017F]{$(ule[1]),$(ule[2])})" :
        ((df[i] == 'U') && (repeat_count == 1)) ? "(?<U>[A-Za-z\u00C0-\u017F]{$(Ule[1]),$(Ule[2])})" :
        ((df[i] == 'e') && (repeat_count == 1)) ? "(?<e>[A-Za-z\u00C0-\u017F]{$(ele[1]),$(ele[2])})" :
        ((df[i] == 'E') && (repeat_count == 1)) ? "(?<E>[A-Za-z\u00C0-\u017F]{$(Ele[1]),$(Ele[2])})" :
        ((df[i] == 'd') && (repeat_count == 1) && !repeat_next) ? "(?<d>0?[1-9]|[12]\\d|3[01])" :
        ((df[i] == 'd') && (repeat_count == 2) && !repeat_next) ? "(?<d>0[1-9]|[12]\\d|3[01])" :
        ((df[i] == 'd') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<d>0[1-9]|[12]\\d|3[01])" :
        ((df[i] == 'H') && (repeat_count == 1) && !repeat_next) ? "(?<H>0?\\d|1\\d|2[0-3])" :
        ((df[i] == 'H') && (repeat_count == 2) && !repeat_next) ? "(?<H>0\\d|1\\d|2[0-3])" :
        ((df[i] == 'H') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<H>0\\d|1\\d|2[0-3])" :
        ((df[i] == 'M') && (repeat_count == 1) && !repeat_next) ? "(?<M>\\d|[0-5]\\d)" :
        ((df[i] == 'M') && (repeat_count == 2) && !repeat_next) ? "(?<M>[0-5]\\d)" :
        ((df[i] == 'M') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<M>[0-5]\\d)" :
        ((df[i] == 'S') && (repeat_count == 1) && !repeat_next) ? "(?<S>\\d|[0-5]\\d)" :
        ((df[i] == 'S') && (repeat_count == 2) && !repeat_next) ? "(?<S>[0-5]\\d)" :
        ((df[i] == 'S') && (repeat_count > 2) && !repeat_next) ? "0{$(repeat_count-2)}(?<S>[0-5]\\d)" :
        ((df[i] == '.') && dotsec) ? "" :
        ((df[i] == '.')) ? "\\." :
        ((df[i] == 's') && (dotsec == true) && (repeat_count < 4) && !repeat_next) ? "(\\.(?<s>\\d{0,3}0{0,6}))?" :
        ((df[i] == 's') && (dotsec == true) && (repeat_count > 3) && !repeat_next) ? "(\\.(?<s>\\d{$(repeat_count)}))?" :
        ((df[i] == 's') && (dotsec == false) && (repeat_count < 4) && !repeat_next) ? "(?<s>\\d{3})?" :
        ((df[i] == 's') && (dotsec == false) && (repeat_count > 3) && !repeat_next) ? "(?<s>\\d{$(repeat_count)})?" :
        ((df[i] == 'z') && !repeat_next) ? "(?<z>[\\+|\\-]?(0\\d|1\\d|2[0-3]):?[0-5]\\d)" :
        ((df[i] == 'Z') && !repeat_next) ? "(?<Z>[A-Z]{3,14})" :
        in(df[i], codechars) ? "" : string(df[i])
        )
    end
    return Regex(r * " *" * string('$'))
end


text:
"Structure for single observation"
struct Observation{S<:Real,T<:AbstractVector{S}} <: AbstractVector{S}
    X::T
    dt::S
    npoints
end

Observation(X,dt) = Observation(X,dt,length(X))

size(a::Observation) = (a.npoints,)

function getindex(x::Observation, i::Integer)
    checkbounds(x, i)
    return Base.getindex(x.X, i)
end
getindex(x::Observation, I::AbstractVector{Bool}) = x.X[findall(I)]
getindex(x::Observation, I::AbstractArray{Bool}) = x.X[LinearIndices(I)[findall(I)]]
getindex(x::Observation, ::Colon) = copy(x.X)

"Structure for ensemble observation"
struct EnsembleObservation{S<:Real,T<:AbstractMatrix{S}} <: AbstractMatrix{S}
    X::T
    dt::S
    npoints
    nsample
end

EnsembleObservation(X,dt) = EnsembleObservation(X,dt,size(X,1),size(X,2))

size(a::EnsembleObservation) = (a.npoints,a.nsample)

function getindex(x::EnsembleObservation, i::Integer)
    checkbounds(x, i)
    return Base.getindex(x.X, i)
end

function getindex(x::EnsembleObservation, i::Integer, j::Integer)
    checkbounds(x, i, j)
    return Base.getindex(x.X, i, j)
end
getindex(x::EnsembleObservation, ::Colon) = copy(x.X)


text:
module Donquicote

using PrettyTables, Statistics, StatsBase, ForwardDiff

D(f) = x -> ForwardDiff.derivative(f,float(x))

function newton(f, ∂f, x0; xatol = 1e-7)

    diff = 1
    while diff >= xatol
        x1 = x0 - f(x0) / ∂f(x0);
        diff = abs(x1 - x0);
        x0 = x1;
    end
    return x0
end

function haley(f, ∂f, ∂2f, x0; xatol = 1e-7)
    diff = 1
    while diff >= xatol
        x1 = x0 - (2*f(x0)*∂f(x0)) / (2*∂f(x0)^2 - f(x0)*∂2f(x0));
        diff = abs(x1 - x0);
        x0 = x1;
    end
    return x0

end

function schroeder(f, ∂f, ∂2f, x0; xatol = 1e-7)
    diff = 1
    while diff >= xatol
        x1 = x0 - f(x0) / ∂f(x0)  - (∂2f(x0)*f(x0)^2)/(2*∂f(x0)^3);
        diff = abs(x1 - x0);
        x0 = x1;
    end
    return x0

end

function weightedmean(x, w)
    wxsum = wsum = 0
    for (x,w) in zip(x,w)
        wx = w*x
        if !ismissing(wx) | isfinite(wx)
            wxsum += wx
            wsum += w
        end
    end
    return wxsum / wsum
end

function weightedstd(x, w)
    μ = weightedmean(x, w)
    wxsum = wsum = 0
    for (x,w) in zip(x,w)
        wx = w*(x - μ)^2
        if !ismissing(wx) | isfinite(wx)
            wxsum += wx
            wsum += w
        end
    end
    σ = sqrt((1/wsum)*wxsum)
    return σ
end

function gini(x::Vector{Float64})
    n = length(x)
    sarray = sort(x)
    return 2*(sum(collect(1:n).*sarray))/(n*sum(sarray))-1
end

logm(x) = x > 0.0 ? log(x) : missing

lvar(x) = var(logm.(x))

nanmean(x) = mean(filter(!isnan,x))

format(x) = round(x, digits = 3)

function savehtml(filename, data)
    open("$filename.html", "w") do f
        pretty_table(f, data, backend = :html)
    end
end

export newton, haley, schroeder
export gini, nanmean, weightedmean, weightedstd, format, logm, lvar
export ones, zeros
export savehtml

end