text:
using Random

mutable struct GridWorld{Int,Int} <: AbstractEnvironment
    state::Int
    reward::Float64
    done::Bool
    size::Int
    pos::Tuple{Int, Int}#Array{Int, 1}
    time_step::Int
	stochastic::Bool
	randchance::Float64
    function GridWorld(size::Int, stochastic::Bool=false, chance::Float64=0.2)
        new{Int,Int}(1, 0., false, size, (1,1), 0, stochastic, chance)
    end
end

function get_gamma(env::GridWorld)
	return 1.0
end

function randomize_env(::Type{T}, rng::AbstractRNG, size::Int, stochastic::Bool) where {T <: GridWorld}
	ϵ = rand(rng, Uniform(0., 0.25))
	return GridWorld(size, stochastic, ϵ)
end

function get_state_desc(env::GridWorld)::Int
    return env.size^2
end

function get_action_desc(env::GridWorld)::Int
    return 4
end

function step!(env::GridWorld, action::Int, rng::AbstractRNG)::Float64
    env.time_step += 1
    if action <= 0 || action > 4
        error("Action needs to be an integer in [1, 4]")
    end
    x,y = env.pos

    if env.stochastic
		temp = rand(rng)
		noeffect = 1. - env.randchance  # e.g., 0.8
		stay = 1. - (env.randchance / 2.) # e.g., 0.9
		side = 1. - (env.randchance / 4.) # e.g., 0.95
		if (temp < noeffect)  # Take proposed action
			# do nothing
		elseif (temp < stay)
			action = 0    # Stay put. Encode as action -1
		elseif (temp < side)
			action += 1
			if (action == 5)
				action = 1
			end
		else
			action -= 1
			if (action == 0)
				action = 4
			end
		end
	end

    if action == 1
        y -= 1
    elseif action == 2
        y += 1
    elseif action == 3
        x -= 1
    elseif action == 4
        x += 1
    end
    x = clamp(x, 1, env.size)
    y = clamp(y, 1, env.size)
    env.pos = (x, y)
    env.reward = -1.0
    env.state = (x-1)*env.size + y
	max_steps = 20 * env.size^2
    env.done = (env.pos == (env.size, env.size)) | (env.time_step ≥ max_steps)

    return env.reward
end

function get_state(env::GridWorld)::Int
    return env.state
end

function get_reward(env::GridWorld)::Float64
    return env.reward
end

function is_terminal(env::GridWorld)::Bool
    return env.done
end

function env_reset!(env::GridWorld, rng)
    env.pos = (1,1)
    env.time_step = 0
    env.state = (env.pos[1]-1)*env.size + (env.pos[2])
    env.done = false
end

function clone(env::GridWorld)::GridWorld
    env2 = GridWorld(env.size, env.stochastic, env.randchance)
    env2.pos = deepcopy(env.pos)
    env2.time_step = deepcopy(env.time_step)
    env2.state = deepcopy(env.state)
    env2.done = deepcopy(env.done)
    env2.reward = deepcopy(env.reward)
    return env2
end

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
using Base: splat

test1 = """
    class: 1-3 or 5-7
    row: 6-11 or 33-44
    seat: 13-40 or 45-50

    your ticket:
    7,1,14

    nearby tickets:
    7,3,47
    40,4,50
    55,2,20
    38,6,12
    """

test2 = """
    class: 0-1 or 4-19
    row: 0-5 or 8-19
    seat: 0-13 or 16-19

    your ticket:
    11,12,13

    nearby tickets:
    3,9,18
    15,1,5
    5,14,9
    """

function parse_input(input::AbstractString)
    (specs, yours, nearby) = strip.(split(input, r"(your|nearby) tickets?:\n"))
    specs = [let (field, ranges...) = split(spec, r": | or "),
                 ranges = [let (a, b) = parse.(Int, split(range, '-'))
                               a:b
                           end for range in ranges]
                 field => union(ranges...)
             end for spec in split(specs, '\n')]
    yours = parse.(Int, split(yours, ','))
    nearby = [parse.(Int, split(nb, ',')) for nb in split(nearby, '\n')]
    return (specs, yours, nearby)
end

function findinvalid(ticket::AbstractVector, specs::AbstractVector)::Vector{Int}
    valid = union(last.(specs)...)
    return findall(∉(valid), ticket)
end

q1(specs, _, nearby) = sum(nb -> nb[findinvalid(nb, specs)] |> sum, nearby)

isvalid(specs) = ticket -> isempty(findinvalid(ticket, specs))

function q2(specs, yours, nearby)::Int
    fields = find_fields(specs, nearby)
    deps = findall(startswith("departure"), first.(specs))
    @assert length(deps) == 6
    return prod(yours[indexin(deps, fields)])
end

function find_fields(specs, nearby)::Vector{Int}
    nearby = filter(isvalid(specs), nearby) |> splat(hcat)
    ranges = last.(specs)
    fields = map(eachrow(nearby)) do vs
        return findall(r -> all(vs .∈ Ref(r)), ranges)
    end
    while any(>(1) ∘ length, fields)
        unamb_idx = findall(==(1) ∘ length, fields)
        unamb_vs = union(fields[unamb_idx]...)
        setdiff!.(fields[setdiff(1:end, unamb_idx)], Ref(unamb_vs))
    end
    return only.(fields)
end

let input1 = parse_input(test1),
    (specs2, _, nearby2) = parse_input(test2)

    @assert q1(input1...) == 71
    @assert find_fields(specs2, nearby2) == [2, 1, 3]
end

input = parse_input(read("day16.in", String))
println("Q1: ", q1(input...))
println("Q2: ", q2(input...))

ext: jl

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
function padding(p)

    θlim, slim = size(p);
    θ_left = Int(θlim/2+1);
    θ_right = Int(θlim/2+θlim);
    s_left = Int(slim/2+1);
    s_right = Int(slim/2+slim);

    p_Padding = zeros(θlim*2, slim*2);
    p_Padding[θ_left:θ_right, s_left:s_right] = p;

    return p_Padding;
end

function conical_radon(p, ϕ, μ)

    ratio = 2;
    p = padding(p);
    θlim, slim = size(p);
    θ = LinRange(-2, 2, θlim);
    U = LinRange(-2, 2, θlim);
    V = LinRange(2, -2, slim);
    C = zeros(length(V), length(U));
    interval = 2*ratio / (θlim - 1);

    for Cx = 1:length(U), Cy=1:length(V)

        U_left = 1:Cx;
        y_left = tan(π / 2 + ϕ) * (θ[U_left] .- U[Cx]) .+ V[Cy];
        bit_left = abs.(y_left) .<= 1;

        if sum(bit_left)>0
            y_left = y_left[bit_left];
            height_left = ratio .- y_left;
            nearest_row_left = round.(Int64, ((height_left) - mod.(height_left, interval)) / interval) .+ 1;
            U_left = U_left[bit_left];
            for i=1:length(nearest_row_left)
                value_left = p[nearest_row_left[i], U_left[i]];
                r_left = sqrt.((value_left-U[Cx])^2 + (value_left-V[Cy])^2);
               C[Cy, Cx]+= value_left*exp(-μ*r_left);
            end
        end

        U_right = Cx:length(U);
        y_right = tan(π / 2 - ϕ) * (θ[U_right] .- U[Cx]) .+ V[Cy];
        bit_right = abs.(y_right) .<= 1;

        if sum(bit_right)>0
            y_right = y_right[bit_right];
            height_right = ratio .- y_right;
            nearest_row_right = round.(Int64, ((height_right) - mod.(height_right, interval)) / interval) .+ 1;
            U_right = U_right[bit_right];
            for i=1:length(nearest_row_right)
                value_right = p[nearest_row_right[i], U_right[i]];
                r_right = sqrt.((value_right-U[Cx])^2 + (value_right-V[Cy])^2);
               C[Cy, Cx]+= value_right*exp(-μ*r_right);
            end
        end
    end

    return C
end


text:
# This file was generated, do not modify it. # hide
#hideall
using CSV, DataFrames, Latexify, Formatting
D = CSV.read("emoji15perf.tsv", DataFrame)
println(latexify(D, latex=false, env=:mdtable, fmt = p -> format(p, commas=true )), " ")


text:
data = split(
"""
.#..##.###...#######
##.############..##.
.#.######.########.#
.###.#######.####.#.
#####.##.#.##.###.##
..#####..#.#########
####################
#.####....###.#.#.##
##.#################
#####.##.###..####..
..######..##.#######
####.##.####...##..#
.#####..#.######.###
##...#.##########...
#.##########.#######
.####.#.###.###.#.##
....##.##.###..#####
.#.#.###########.###
#.#.#.#####.####.###
###.##.####.##.#..##
""", "\n")

data = split(
"""
......#.#.
#..#.#....
..#######.
.#.#.###..
.#..#.....
..#....#.#
#..#....#.
.##.#..###
##...#..#.
.#....####
""", "\n")

function maketuples(data)
    a = Vector{Tuple{Int, Int}}()
    for (i, line) in enumerate(data)
        for j in findall(==("#"), split(line, ""))
            push!(a, (i,j))
        end
    end
    a
end

function countseen(a::Vector, i)
    seen = Set{Tuple{Float64, Float64}}()
    for asts in a
        asts == i && continue
        dir = asts .- i
        dirp = dir[1] == 0 ? dir ./ abs(dir[2]) : dir ./ abs(dir[1])
        push!(seen, dirp)
    end
    length(seen)
end

data = readlines("data/Dec10_data.txt")
a = maketuples(data)
t = countseen(a, a[1])
findmax(map(i->countseen(a, i), a))

using DataStructures
function seenDict(a::Vector, i)
    seen = MultiDict{Tuple{Float64, Float64}, Tuple{Int, Int}}()
    for asts in a
        asts == i && continue
        dir = asts .- i
        dirp = dir[1] == 0 ? dir ./ abs(dir[2]) : dir ./ abs(dir[1])
        insert!(seen, dirp, dir)
    end
    seen
end

function shootstars(myast, a)
    ret = Tuple{Int, Int}[]
    sd = seenDict(a, myast)
    for (k,v) in sd
        sort!(v, by = x->sum((x .- myast).^2), rev = true)
    end
    sortedkeys = sort(collect(keys(sd)), by = x->atan(x[2], x[1]), rev = true)
    counter = 0
    for key in sortedkeys
        if !isempty(sd[key])
            p = pop!(sd[key])
            push!(ret, p .+ myast)
            (counter += 1) == 200 && return p .+ myast .- (1,1)
        end
    end
end
res = shootstars(a[353], a)
res[2]*100 + res[1]


text:
module TestSuite
    using Test, AssociatedLegendrePolynomials
    const LMAX = 5
    const NumTypes = (Float32, Float64, BigFloat)

    function runtests(norm)
        @testset "Normalization $(typeof(norm))" begin
            interface(norm)
        end
    end

    function interface(norm)
        # Test all the aspects of the interface that each normalization type must
        # implement
        @testset "Legendre Interface" begin
            @test norm isa AbstractLegendreNorm
            # Initial condition
            @test isfinite(AssociatedLegendrePolynomials.initcond(norm, Float64))
            # 1-term recurrence coefficients
            @test isfinite(AssociatedLegendrePolynomials.coeff_μ(norm, Float64, 1))
            @test isfinite(AssociatedLegendrePolynomials.coeff_ν(norm, Float64, 1))
            # 2-term recurrence coefficients
            @test isfinite(AssociatedLegendrePolynomials.coeff_α(norm, Float64, 1, 0))
            @test isfinite(AssociatedLegendrePolynomials.coeff_β(norm, Float64, 1, 0))

            # Check that the boundscheck_hook() returns `nothing`
            @test @inferred(AssociatedLegendrePolynomials.boundscheck_hook(norm, 0, 0)) === nothing

            # Make sure object calls have not been shadowed
            Λ₁ = zeros(LMAX+1, LMAX+1)
            Λ₂ = zeros(LMAX+1, LMAX+1)
            @test norm(    LMAX, LMAX, 0.5) == legendre( norm,     LMAX, LMAX, 0.5)
            @test norm(Λ₁, LMAX, LMAX, 0.5) == legendre!(norm, Λ₁, LMAX, LMAX, 0.5)
        end
    end
end


text:
using .Bodies: G_year_AU, Body

struct LeapFrog end
struct PEFRL end
struct FR end

function leapfrog(::Type{LF}, planets::AbstractVector{Body{T}}; G=G_year_AU, n, dt, keep_history=false) where {T, LF}
    nplanets = length(planets)
    m = getfield.(planets, :m)
    @inbounds if keep_history
        v = zeros(V3{T}, (nplanets,n+1))
        r = zeros(V3{T}, (nplanets,n+1))
        v[:,1] .= getfield.(planets, :v)
        r[:,1] .= getfield.(planets, :r)
    else
        v = zeros(V3{T}, nplanets)
        r = zeros(V3{T}, nplanets)
        v .= getfield.(planets, :v)
        r .= getfield.(planets, :r)
    end
    leapfrog_looper!(LF, r, v, m; G=G, dt=dt, n=n)
end

function leapfrog_looper!(::Type{LF}, r::AbstractVector{V3{T}}, v::AbstractVector{V3{T}},
         m::AbstractVector; G=G_year_AU, n, dt) where {T, LF}
    for i=1:n
        steper(LF)(r, v, m; G=G, dt=dt)
    end
    return r, v
end

function leapfrog_looper!(::Type{LF}, r::AbstractMatrix{V3{T}}, v::AbstractMatrix{V3{T}},
         m::AbstractVector; G=G_year_AU, n, dt) where {T,LF}
    for i=1:n
        v[:,i+1] .= view(v,:,i)
        r[:,i+1] .= view(r,:,i)
        steper(LF)(view(r,:,i+1), view(v,:,i+1), m; G=G, dt=dt)
    end
    return r, v
end

@inline function acceleration(ra, rb, mb, G)
    d = distance(ra, rb)
    (G*mb/d^3) * (rb - ra)
end

function update_acceleration!(a::AbstractVector{V3{T}}, r, m, G) where T
    @assert length(a) == length(m)
    @assert length(a) == length(r)
    @inbounds for j=1:length(m)
        a[j] = zero(V3{T})
        for k=1:length(m)
            if j != k
                a[j] += acceleration(r[j], r[k], m[k], G)
            end
        end
    end
    return a
end


text:
# This file is a part of the ExamplePackage.

using ExamplePackage
using Base.Test

@testset "ExamplePackage.jl" begin
    include("test_sort_numbers.jl")
end


text:
@symbol_func function kappa_x(cur_reactor::AbstractReactor)
  cur_kappa = cur_reactor.kappa_95

  cur_kappa *= 1.12

  cur_kappa
end


text:
import jInv.ForwardShare.getSensTMatVec

function getSensTMatVec(v::Vector,m::Vector,pFor::BasicFWIparam)

    # extract pointers
    Mesh  = pFor.Mesh
    omega = pFor.omega
    gamma = pFor.gamma
    Q     = pFor.Sources
    P     = pFor.Receivers
    LU    = pFor.Ainv
    U     = pFor.Fields

    nsrc = size(Q,2); nfreq = length(omega); nrec = size(P,2)

    # reshape v by receivers x sources x frequencies
    v    = reshape(v,nrec,nsrc,nfreq)

    # allocate space for result
    JTv  = zeros(length(m))
    # derivative of mass matrix
    An2cc = getNodalAverageMatrix(Mesh)
    dM    = repeat(1 .- 1im*vec(gamma),1,nsrc)
    for i=1:nfreq
        tf1=time_ns();
        Lam  = LU[i]\(P*v[:,:,i])
        # println(size(An2cc))
        # println(size(Lam))
        # println(size(U[:,:,i]))
        # println(size(dM))
        # println(size(omega[i]))
        JTvi =  omega[i]^2*dM.*(An2cc*(U[:,:,i].*Lam))
        JTv +=  sum(real(JTvi),dims=2)
        tf2=time_ns();
        # println("Runtime of sensT freq:");
        # println((tf1 - tf2)/1.0e9);
    end
    return vec(JTv)
end

function getSensTMatVecZs(v::Vector,m::Vector,pFor::BasicFWIparam)

        # extract pointers
        Mesh  = pFor.Mesh
        omega = pFor.omega
        gamma = pFor.gamma
        Q     = pFor.Sources
        P     = pFor.Receivers
        LU    = pFor.Ainv
        U     = pFor.Fields

        nsrc = size(Q,2); nfreq = length(omega); nrec = size(P,2)

        # reshape v by receivers x sources x frequencies
        v    = reshape(v,nrec,nsrc,nfreq)

        # allocate space for result
        JTv  = zeros(length(m))
        # derivative of mass matrix
        An2cc = getNodalAverageMatrix(Mesh)
        dM    = repeat(1 .- 1im*vec(gamma),1,1)
        for i=1:nfreq
            tf1 = time_ns();
            for s=1:nsrc
                t1 = time_ns();
                Lam  = LU[i]\(P*v[:,s,i])
                JTv += omega[i]^2*dM.*(An2cc*((LU[i]\Vector(Q[:,s])).*Lam))[:,1];
                t2 = time_ns();
                # println("Runtime of sensT iteration:");
                # println((t1 - t2)/1.0e9);
            end
            tf2=time_ns();
            # println("Runtime of sensT freq:");
            # println((tf1 - tf2)/1.0e9);

        end
        return vec(JTv)
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "integration" begin
    include("test_forwarddiff.jl")
    include("test_zygote.jl")
end


text:
# # Quadratic Program
#
# We want to solve the following quadratic program with decision variable `x`:
# $$
# \begin{array}{ll} \text{minimize} &  1/2 x^\top P x + q^\top x \\
# \text{subject to} &  l \leq A x \leq u
# \end{array}
# $$
# The problem can be solved with `COSMO` in the following way. Start by defining the problem data

using COSMO, SparseArrays, LinearAlgebra, Test

q = [1; 1.];
P = sparse([4. 1; 1 2]);
A = [1. 1; 1 0; 0 1];
l = [1.; 0; 0];
u = [1; 0.7; 0.7];

# First, we decide to solve the problem with two one-sided constraints using `COSMO.Nonnegatives` as the convex set:
Aa = [-A; A]
ba = [u; -l]
constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);

# Next, we define the settings object, the model and then assemble everything:
settings = COSMO.Settings(verbose=true);
model = COSMO.Model();
assemble!(model, P, q, constraint1, settings = settings);
res = COSMO.optimize!(model);

# Alternatively, we can also use two-sided constraints with `COSMO.Box`:
constraint1 = COSMO.Constraint(A, zeros(3), COSMO.Box(l, u));

model = COSMO.Model();
assemble!(model, P, q, constraint1, settings = settings);
res_box = COSMO.optimize!(model);

# Let's check that the solution is correct:
@testset "QP Problem" begin
  @test norm(res.x[1:2] - [0.3; 0.7], Inf) < 1e-3
  @test norm(res_box.x[1:2] - [0.3; 0.7], Inf) < 1e-3
  @test abs(res.obj_val - 1.88) < 1e-3
  @test abs(res_box.obj_val - 1.88) < 1e-3
end
nothing


text:
@testset "402.remove-k-digits.jl" begin
    @test remove_kdigits("1432219", 3) == "1219"
    @test remove_kdigits("10200", 1) == "200"
    @test remove_kdigits("10", 1) == "0"
    @test remove_kdigits("112", 1) == "11"
end


text:
using CovarianceMatrices
using Random, DataFrames, GLM
using StanSample
using StatisticalRethinking

Random.seed!(1)
n = 500
x = randn(n, 5)
u = zeros(2*n)
u[1] = rand()
for j in 2:2*n
    u[j] = 0.78*u[j-1] + randn()
end
u = u[n+1:2*n]
y = 0.1 .+ x*[0.2, 0.3, 0.0, 0.0, 0.5] + u

df = DataFrame()
for i in 1:5
    df[!, Symbol("x$i")] = x[:, i]
end
df[!,:y] = y

#Using the data in df, the coefficient of the regression can be estimated using GLM

lm1 = glm(@formula(y~x1+x2+x3+x4+x5), df, Normal(), IdentityLink())
lm1 |> display
println()

vcov(QuadraticSpectralKernel{Andrews}(), lm1, prewhite = false) |> display
println()

# For the previous example:

stderror(QuadraticSpectralKernel{Andrews}(), lm1, prewhite = false) |> display
println()

stderror(QuadraticSpectralKernel{NeweyWest}(), lm1, prewhite = false) |> display
println()

# Sometime is useful to access the bandwidth selected by the automatic procedures.
# This can be done using the optimalbandwidth method

optimalbandwidth(QuadraticSpectralKernel{NeweyWest}(), lm1; prewhite = false) |> display
println()

optimalbandwidth(QuadraticSpectralKernel{Andrews}(), lm1; prewhite = false) |> display
println()


text:
# Op_fixed initializes all the c^dagger c operators
# in the fixed particle number subspace
struct Op_fixed
    dim::Int
    nume::Int
    cdctot::SparseMatrixCSC{Float64,Int64}
    le::Int
    basis::SparseMatrixCSC{Float64,Int64}
    Op_fixed(dim,nume) = new(dim, nume, operators_fixed(dim,nume)...)
end

dim(o::Op_fixed) = o.dim
nume(o::Op_fixed) = o.nume
cdctot(o::Op_fixed) = o.cdctot
le(o::Op_fixed) = o.le
basis(o::Op_fixed) = o.basis

ada(o::Op_fixed,i,j) = cdctot(o)[(1+(i-1)*le(o)):i*le(o),(1+(j-1)*le(o)):j*le(o)]
aad(o::Op_fixed,i,j) = ada(o,j,i)

# Here we create the full matrix with all
# fixed particle operators definitions
function operators_fixed(n::Int,m::Int)
    l1 = binomial(n,m)
    l2 = n*l1
    b,ind = basis_m(n,m)
    z = spzeros(l2,l2)
    for i in 1:n
        for j in i:n
            if i==j
                z[(1+(i-1)*l1):i*l1,(1+(j-1)*l1):j*l1] = 1/2*ada(b,ind,i,j)
            else
                z[(1+(i-1)*l1):i*l1,(1+(j-1)*l1):j*l1] = ada(b,ind,i,j)
            end
        end
    end
    op_general = z+z'
    return op_general, l1, b
end

#operator c^dagger c for a system with fixed particle number.
# Much faster than using fiexed over the whole operators

function ada(n::Int64, m::Int64, i::Int64, j::Int64)
    base, ind = basis_m(n,m)
    indice = ind .-1
    l = binomial(n,m)
    op = spzeros(l,l)
    if i==j
        for k in 1:l
            if base[k,:][j] == 1
                op[k,k] = 1
            end
        end
    else
        for k in 1:l
            if base[k,:][j] == 1 && base[k,:][i] == 0
                ill = Int(indice[k]-2^(n-j)+2^(n-i))
                ill2 = myfind(indice, ill)[1]
                sign = (-1)^(sum(base[k,1:(j-1)])+sum(base[ill2,1:(i-1)]))
                op[ill2,k] = sign
            end
        end
    end
    return op
end

# When using many operators, it is easier to predefine the basis outside
# the function
function ada(base::SparseArrays.SparseMatrixCSC{Float64,Int64}, ind::SparseArrays.SparseVector{Float64,Int64}, i::Int64, j::Int64)
    l = size(base)[1]
    n = size(base)[2]
    op = spzeros(l,l)
    indice = ind .-1
    if i==j
        for k in 1:l
            if base[k,:][j] == 1
                op[k,k] = 1
            end
        end
    else
        for k in 1:l
            if base[k,:][j] == 1 && base[k,:][i] == 0
                ill = Int(indice[k]-2^(n-j)+2^(n-i))
                ill2 = myfind(indice, ill)[1]
                sign = (-1)^(sum(base[k,1:(j-1)])+sum(base[ill2,1:(i-1)]))
                op[ill2,k] = sign
            end
        end
    end
    return op
end

# the c c^dagger operator, is the
function aad(n::Int64, m::Int64, i::Int64, j::Int64)
    return -ada(n,m,j,i)
end

function aad(base::SparseArrays.SparseMatrixCSC{Float64,Int64}, indice::SparseArrays.SparseVector{Float64,Int64}, i::Int64, j::Int64)
    return -ada(base,indice,j,i)
end

# This is the creation operator adagger
# it serves as ancilla creation or in the same
# sp space. m is the number of particles in The
# starting state
function ad(n::Int64, m::Int64, i::Int64)
    base1, ind1 = basis_m(n,m)
    l2 = binomial(n,m)
    # if I am in the original sp space
    if i<=n
        base2, ind2 = basis_m(n,m+1)
        l1 = binomial(n,m+1)
        op = spzeros(l1,l2)
        for k in 1:l2
            if base1[k,:][i] == 0
                j = myfind(ind2,(ind1[k]+2^(n-i)))[1]
                sign = (-1)^(sum(base1[k,1:(i-1)]))
                op[j,k] = sign
            end
        end
    # else, create an ancilla
    else
        base2, ind2 = basis_m(i,m+1)
        l1 = binomial(i,m+1)
        op = spzeros(l1,l2)
        for k in 1:l2
            x = join(Int.(append!(Array(base1[k,:]),1)))
            x = parse(Int, x; base=2) + 1
            j = myfind(ind2,x)[1]
            sign = (-1)^m
            op[j,k] = sign
        end
    end
    return op
end

# This is the destruction operator.
# m is the number of particles in The
# starting state. It can also be built from ad
# computing the transpose
function a(n::Int64, m::Int64, i::Int64)
    if i <= n
        base1, ind1 = basis_m(n,m)
        base2, ind2 = basis_m(n,m-1)
        l1 = binomial(n,m-1)
        l2 = binomial(n,m)
        op = spzeros(l1,l2)
        for k in 1:l2
            if base1[k,:][i] != 0
                j = myfind(ind2,(ind1[k]-2^(n-i)))[1]
                sign = (-1)^(sum(base1[k,1:(i-1)]))
                op[j,k] = sign
            end
        end
        return op
    else
        return spzeros(binomial(n,m))
    end
end

# This serves for applying several creations op
# at once. Modes will be a vector with the chosen
# modes
function ad(n::Int64, m::Int64, modes::Array{Int64,1})
    l = length(modes)
    op  = ad(n,m,modes[1])
    for k in 2:l
        op = ad(n,m+k-1,modes[k])*op
    end
    return op
end

# This serves for applying several destruction op
# at once. Modes will be a vector with the chosen
# modes
function a(n::Int64, m::Int64, modes::Array{Int64,1})
    l = length(modes)
    op  = a(n,m,modes[1])
    for k in 2:l
        op = a(n,m-k+1,modes[k])*op
    end
    return op
end

# The following function initializes the creation
# operators for each mode given a fixed dimension d
# and the starting number of particles m
function adf(d::Int, m::Int)
    adtot = [ad(d,m,i) for i in 1:d]
    return adtot
end

# The following function initializes the destruction
# operators for each mode given a fixed dimension d
# and the starting number of particles m
function af(d::Int, m::Int)
    atot = [a(d,m,i) for i in 1:d]
    return atot
end

#This functions trimms the operators matrix
#for fixed particle number subspace
function fixed(o, nume::Int64)
    n = Int(log(2,size(o)[1]))
    posn = ffilter(n, nume)
    reduced = o[setdiff(1:end,posn), setdiff(1:end,posn)]
    return reduced
end

function ffilter(n::Int64, m::Int64)
    _,_,baso,_ = integer_digits(n)
    posn = zeros(Int64, Int(2^n-binomial(n,m)))
    countern = 1
    for i in 1:2^n
        if sum(baso[i,:]) != m
            posn[countern] = Int(i)
            countern = countern + 1
        end
    end
    return posn
end

function myfind(c,j)
    a = similar(c, Int)
    count = 1
    @inbounds for i in eachindex(c)
        a[count] = i
        count += (c[i] == j)
    end
    return a[1:(count-1)]
end

## What follows is the op_semifixed struct
# that I previously used for initalizing all
# creation and destruction op in fixed
# particle number subspaces, switching
# from one to the other.
# I prefer now initalizing only adf and af
# as it is clearer and spanns most cases

# You should add LazyStack to the dependancies:
#LazyStack = "1fad7336-0346-5a1a-a56f-a06ba010965b"
# in Fermionic.jl:
# using LazyStack

# # The following struct initialized the creation
# # and destruction operators, given a diemension
# struct Op_semifixed
#     d::Int
#     adtot::Array{Array{Float64,3},1}
#     Op_semifixed(d) = new(d, opnfixed(d))
# end
#
# ad(o::Op_semifixed, m::Int, i::Int) = sparse(o.adtot[m][:,:,i])
# # a(o::Op_semifixed, m::Int, i::Int) = sparse(o.adtot[m-1][:,:,i]')
#
# function a(o::Op_semifixed, m::Int, i::Int)
#     if i <= o.d
#         return sparse(o.adtot[m-1][:,:,i]')
#     else
#         return spzeros(binomial(o.d,m))
#     end
# end

# # This function is the core of Op_semifixed
# # initialization. It uses rstack from the package
# # Lazystack. It looses the sparse structure in doing so
# # A pending task is to improve this mechanism
# # We only initialize ad, and compute a as the transpose
# # Changing this function would allow deleting
# # a package from the dependancies of the package
# function opnfixed(d::Int)
#     r = [rstack([ad(d,j,i) for i in 1:d];fill=0.) for i in 1:d for j in 1:(d-1)]
#     return r
# end

# # This serves for applying several creations op
# # at once. Modes will be a vector with the chosen
# # modes
# function ad(o::Op_semifixed, m::Int64, modes::Array{Int64,1})
#     l = length(modes)
#     op  = ad(o,m,modes[1])
#     for k in 2:l
#         op = ad(o,m+k-1,modes[k])*op
#     end
#     return op
# end

# # This serves for applying several destruction op
# # at once. Modes will be a vector with the chosen
# # modes
# function a(o::Op_semifixed, m::Int64, modes::Array{Int64,1})
#     l = length(modes)
#     op  = a(o,m,modes[1])
#     for k in 2:l
#         op = a(o,m-k+1,modes[k])*op
#     end
#     return op
# end

# The documentation:
# - **Op_semifixed(d::Int)**: initializes the fermionic operators $a_i$ and $a_i^\dagger$ corresponding to dimension **d** in the fixed particle subspace. Once initizialized, the following functions are made avaiable.
#     - **a(o::Op_semifixed, n::Int, i::Int)**: returns the fermionic destruction operator in mode i, i.e. $a_i$, corresponding to the previously defined **o** in the fixed particle subspace with $n$ starting particles.
#     - **a(o::Op_semifixed,  n::Int, mod::Array{Int,1})**: returns the consecutive fermionic destruction operator of modes in vecor **mod**, i.e. $a_{mod_m}...a_{mod_1}$, corresponding to the previously defined **o** in the fixed particle subspace with $n$ starting particles.
#     - **ad(o::Op_semifixed,  n::Int, i::Int)**: returns the fermionic creation operator in mode **i**, i.e. $a^\dagger_i$, corresponding to the previously defined **o** in the fixed particle subspace with $n$ starting particles.
#      - **ad(o::Op_semifixed,  n::Int, mod::Array{Int,1})**: returns the consecutive fermionic creation operator of modes in vecor **mod**, i.e. $a^\dagger_{mod_m}...a^\dagger_{mod_1}$, corresponding to the previously defined **o** in the fixed particle subspace with $n$ starting particles.


text:
using REPL.TerminalMenus
export interactive_guess

function interactive_guess(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: DefaultRecipe
    cnames = keys(model)
    menu = RadioMenu(string.(cnames), pagesize=20)
    choice = request("Choose a component:", menu)
    (choice == -1)  &&  (return nothing)
    interactive_guess(source, pspec, model, cnames[choice])
end

function interactive_guess(source::QSO{T}, pspec::PreparedSpectrum, model::Model, cname::Symbol) where T <: DefaultRecipe
    pnames = collect(keys(GFit.getparams(model[cname])))
    choices = [string(p.name) * (p.index > 0 ? string(p.index) : "") for p in pnames]
    menu = RadioMenu(choices, pagesize=20)
    choice = request("Choose a parameter of component $cname:", menu)
    (choice == -1)  &&  (return nothing)
    interactive_guess(source, pspec, model, cname, pnames[choice])
end

interactive_guess(source::QSO{T}, pspec::PreparedSpectrum, model::Model, cname::Symbol, pname::Symbol) where T <: DefaultRecipe =
    interactive_guess(source, pspec, model, cname, GFit.ParamID(pname, 0))

function interactive_guess(source::QSO{T}, pspec::PreparedSpectrum, model::Model, cname::Symbol, pid::GFit.ParamID) where T <: DefaultRecipe
    if pid.index == 0
        par = getproperty(model[cname], pid.name)
    else
        par = getproperty(model[cname], pid.name)[pid.index]
    end
    orig_val = par.val
    println("Component: $cname, Parameter: $(pid.name)", (pid.index > 0  ?  "[" * string(pid.index) * "]"  :  ""))
    println("Initial value: ", par.val)
    println("('q' to accept current value, 'c' to restore initial value)")
    while true
        evaluate!(model)
        @gp    :interactive_guess title=(source.name * ", " * pspec.orig.label) "set grid" :-
        @gp :- :interactive_guess reverse([Gnuplot.recipe(model)..., Gnuplot.recipe(pspec.data)])
        print("Insert new value: ")
        input = readline()
        if input == "q"
            println("Consider adding the following code to your recipe:")
            println()
            println("function guess_emission_lines!(source::QSO{T}, pspec::PreparedSpectrum, model::Model) where T <: ", string(typeof(source))[5:end-1])
            println("    guess_emission_lines!(parent_recipe(source), pspec, model)")
            println("    if (source.name == \"", source.name, "\")  &&  (pspec.orig.label == \"", pspec.orig.label, "\")")
            println("        model[:$(cname)].$(pid.name)", (pid.index > 0  ?  "[" * string(pid.index) * "]"  :  ""), ".val = ", par.val)
            println("    end")
            println("end")
            println()
            println("(press ENTER to continue)"); readline()
            break
        elseif input == "c"
            par.val = orig_val
            break
        else
            par.val = Meta.parse(input)
        end
    end
    println("Using value $(par.val)")
    evaluate!(model)
    Gnuplot.quit(:interactive_guess)
end


text:
# Methods for the Nonlinear interface

abstract AbstractNonlinearModel <: AbstractMathProgModel
export AbstractNonlinearModel

abstract AbstractNLPEvaluator
export AbstractNLPEvaluator

# methods for AbstractNLPEvaluator
@define_interface begin
    NonlinearModel
    initialize
    features_available
    eval_f
    eval_g
    eval_grad_f
    jac_structure
    hesslag_structure
    eval_jac_g
    eval_jac_prod
    eval_jac_prod_t
    eval_hesslag_prod
    eval_hesslag
    isobjlinear
    isobjquadratic
    isconstrlinear
    obj_expr
    constr_expr
end

# fallback options
isobjlinear(::AbstractNLPEvaluator) = false
isobjquadratic(::AbstractNLPEvaluator) = false
isconstrlinear(::AbstractNLPEvaluator, i::Integer) = false

# getreducedcosts and getconstrduals already have stubs


text:
__precompile__()

module RobustStats

using Compat
using DataFrames
using Printf
using Rmath
using Statistics
using StatsBase
import Base.show

export
    tmean,
    winval,
    winmean,
    winvar,
    wincov,
    trimse,
    trimci,
    idealf,
    pbvar,
    bivar,
    tauloc,
    tauvar,
    outbox,
    msmedse,
    binomci,
    acbinomci,
    sint,
    hpsi,
    onestep,
    bootstrapci,
    bootstrapse,
    mom,
    momci,
    contam_randn,
    trimpb,
    pcorb,
    yuend,
    bisquareWM,
    huberWM,
    trimean,
    iqrn,
    shorthrange,
    shorthrange!,
    scaleS,
    scaleS!,
    scaleQ,
    scaleQ!

include("types.jl")
include("utils.jl")
include("show.jl")
include("functions.jl")
include("location_estimators.jl")
include("dispersion_estimators.jl")
end


text:
"""
    make_adjacency_matrix(i2f)

Make full adjacency matrix based only on the form of n-grams regardless of whether
they are seen in the training data. This usually takes hours for large datasets,
as all possible combinations are considered.

# Obligatory Arguments
- `i2f::Dict`: the dictionary returning features given indices

# Optional Arguments
- `tokenized::Bool=false`:if true, the dataset target is assumed to be tokenized
- `sep_token::Union{Nothing, String, Char}=nothing`: separator token
- `verbose::Bool=false`: if true, more information will be printed

# Examples
```julia
# without tokenization
i2f = Dict([(1, "#ab"), (2, "abc"), (3, "bc#"), (4, "#bc"), (5, "ab#")])
JudiLing.make_adjacency_matrix(i2f)

# with tokenization
i2f = Dict([(1, "#-a-b"), (2, "a-b-c"), (3, "b-c-#"), (4, "#-b-c"), (5, "a-b-#")])
JudiLing.make_adjacency_matrix(
    i2f,
    tokenized=true,
    sep_token="-")
```
"""
function make_full_adjacency_matrix(
    i2f;
    tokenized = false,
    sep_token = nothing,
    verbose = false,
)

    ngrams = [i2f[i] for i = 1:length(i2f)]

    if tokenized && !isnothing(sep_token)
        words = split.(ngrams, sep_token)
    else
        words = split.(ngrams, "")
    end

    n_ngrams = length(ngrams)

    I = Int64[]
    J = Int64[]
    V = Int64[]

    iter = 1:n_ngrams
    if verbose
        pb = Progress(n_ngrams)
    end

    for i in iter
        for j = 1:n_ngrams
            if isattachable(words[i], words[j])
                push!(I, i)
                push!(J, j)
                push!(V, 1)
            end
        end
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    sparse(I, J, V, n_ngrams, n_ngrams, *)
end

"""
    make_combined_adjacency_matrix(data_train, data_val)

Make combined adjacency matrix.

# Obligatory Arguments
- `data_train::DataFrame`: training dataset
- `data_val::DataFrame`: validation dataset

# Optional Arguments
- `grams=3`: the number of grams for cues
- `target_col=:Words`: the column name for target strings
- `tokenized=false`:if true, the dataset target is assumed to be tokenized
- `sep_token=nothing`: separator
- `keep_sep=false`: if true, keep separators in cues
- `start_end_token="#"`: start and end token in boundary cues
- `verbose=false`: if true, more information is printed

# Examples
```julia
JudiLing.make_combined_adjacency_matrix(
    latin_train,
    latin_val,
    grams=3,
    target_col=:Word,
    tokenized=false,
    keep_sep=false
    )
```
"""
function make_combined_adjacency_matrix(
    data_train,
    data_val;
    grams = 3,
    target_col = :Words,
    tokenized = false,
    sep_token = nothing,
    keep_sep = false,
    start_end_token = "#",
    verbose = false,
)

    t, v = make_combined_cue_matrix(
        data_train,
        data_val;
        grams = grams,
        target_col = target_col,
        tokenized = tokenized,
        sep_token = sep_token,
        keep_sep = keep_sep,
        start_end_token = start_end_token,
        verbose = verbose,
    )

    t.A
end


text:
function generate_random_unit_quaternion()

	q0::Float64 = randn()
	q1::Float64 = randn()
	q2::Float64 = randn()
	q3::Float64 = randn()

	q_norm::Float64 = sqrt( q0^2 + q1^2 + q2^2 + q3^2 )

	q0 /= q_norm
	q1 /= q_norm
	q2 /= q_norm
	q3 /= q_norm

	return (q0, q1, q2, q3)
end


text:
"""
Nodal admittance matrix (Ybus) is an N x N matrix describing a power system with N buses. It represents the nodal admittance of the buses in a power system.

The Ybus Struct is indexed using the Bus Numbers, no need for them to be sequential
"""
struct Ybus{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{ComplexF64}
    data::SparseArrays.SparseMatrixCSC{ComplexF64, Int}
    axes::Ax
    lookup::L
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    b::ACBranch,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]
    Y_l = (1 / (get_r(b) + get_x(b) * 1im))
    Y11 = Y_l + (1im * get_b(b).from)
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = -Y_l
    ybus[bus_from_no, bus_to_no] += Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += Y12
    Y22 = Y_l + (1im * get_b(b).to)
    ybus[bus_to_no, bus_to_no] += Y22
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    b::Transformer2W,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]
    Y_t = 1 / (get_r(b) + get_x(b) * 1im)
    Y11 = Y_t
    b = get_primary_shunt(b)
    ybus[bus_from_no, bus_from_no] += Y11
    ybus[bus_from_no, bus_to_no] += -Y_t
    ybus[bus_to_no, bus_from_no] += -Y_t
    ybus[bus_to_no, bus_to_no] += Y_t + (1im * b)
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    b::TapTransformer,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]

    Y_t = 1 / (get_r(b) + get_x(b) * 1im)
    c = 1 / get_tap(b)
    b = get_primary_shunt(b)

    Y11 = (Y_t * c^2)
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = (-Y_t * c)
    ybus[bus_from_no, bus_to_no] += Y12
    #Y21 = Y12
    ybus[bus_to_no, bus_from_no] += Y12
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += Y22 + (1im * b)
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    b::PhaseShiftingTransformer,
    num_bus::Dict{Int, Int},
)
    arc = get_arc(b)
    bus_from_no = num_bus[arc.from.number]
    bus_to_no = num_bus[arc.to.number]
    Y_t = 1 / (get_r(b) + get_x(b) * 1im)
    tap = (get_tap(b) * exp(get_α(b) * 1im))
    c_tap = (get_tap(b) * exp(-1 * get_α(b) * 1im))
    b = get_primary_shunt(b)
    Y11 = (Y_t / abs(tap)^2)
    ybus[bus_from_no, bus_from_no] += Y11
    Y12 = (-Y_t / c_tap)
    ybus[bus_from_no, bus_to_no] += Y12
    Y21 = (-Y_t / tap)
    ybus[bus_to_no, bus_from_no] += Y21
    Y22 = Y_t
    ybus[bus_to_no, bus_to_no] += Y22 + (1im * b)
    return
end

function _ybus!(
    ybus::SparseArrays.SparseMatrixCSC{ComplexF64, Int},
    fa::FixedAdmittance,
    num_bus::Dict{Int, Int},
)
    bus = get_bus(fa)
    bus_no = num_bus[get_number(bus)]
    ybus[bus_no, bus_no] += fa.Y
    return
end

function _buildybus(branches, nodes, fixed_admittances)
    buscount = length(nodes)
    num_bus = Dict{Int, Int}()

    for (ix, b) in enumerate(nodes)
        num_bus[get_number(b)] = ix
    end
    ybus = SparseArrays.spzeros(ComplexF64, buscount, buscount)
    for b in branches
        if get_name(b) == "init"
            throw(DataFormatError("The data in Branch is invalid"))
        end
        get_available(b) && _ybus!(ybus, b, num_bus)
    end
    for fa in fixed_admittances
        get_available(fa) && _ybus!(ybus, fa, num_bus)
    end
    return SparseArrays.dropzeros!(ybus)
end

"""
Builds a Ybus from a collection of buses and branches. The return is a Ybus Array indexed with the bus numbers and the branch names.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
"""
function Ybus(
    branches,
    nodes,
    fixed_admittances = Vector{FixedAdmittance}();
    check_connectivity::Bool = true,
)
    nodes = sort!(collect(nodes), by = x -> get_number(x))
    bus_ax = get_number.(nodes)
    axes = (bus_ax, bus_ax)
    bus_lookup = _make_ax_ref(bus_ax)
    look_up = (bus_lookup, bus_lookup)
    ybus = _buildybus(branches, nodes, fixed_admittances)
    if check_connectivity
        connected = validate_connectivity(ybus, nodes, bus_lookup)
        !connected && throw(DataFormatError("Network not connected"))
    end
    return Ybus(ybus, axes, look_up)
end

"""
Builds a Ybus from the system. The return is a Ybus Array indexed with the bus numbers and the branch names.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
"""
function Ybus(sys::System; check_connectivity::Bool = true)
    branches = get_components(ACBranch, sys)
    nodes = get_components(Bus, sys)
    fixed_admittances = get_components(FixedAdmittance, sys)
    return Ybus(branches, nodes, fixed_admittances; check_connectivity = check_connectivity)
end

"""
Nodal incidence matrix (Adjacency) is an N x N matrix describing a power system with N buses. It represents the directed connectivity of the buses in a power system.

The Adjacency Struct is indexed using the Bus Numbers, no need for them to be sequential
"""
struct Adjacency{Ax, L <: NTuple{2, Dict}} <: PowerNetworkMatrix{Int}
    data::SparseArrays.SparseMatrixCSC{Int, Int}
    axes::Ax
    lookup::L
end

"""
Builds a Adjacency from a collection of buses and branches. The return is an N x N Adjacency Array indexed with the bus numbers.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
"""
function Adjacency(branches, nodes; check_connectivity::Bool = true)
    buscount = length(nodes)
    bus_ax = get_number.(nodes)
    axes = (bus_ax, bus_ax)
    bus_lookup = _make_ax_ref(bus_ax)
    look_up = (bus_lookup, bus_lookup)

    a = SparseArrays.spzeros(Int, buscount, buscount)

    for b in branches
        (fr_b, to_b) = get_bus_indices(b, bus_lookup)
        a[fr_b, to_b] = 1
        a[to_b, fr_b] = -1
        a[fr_b, fr_b] = 1
        a[to_b, to_b] = 1
    end

    if check_connectivity
        connected = validate_connectivity(a, nodes, bus_lookup)
        !connected && throw(DataFormatError("Network not connected"))
    end

    return Adjacency(a, axes, look_up)
end

"""
Builds a Adjacency from the system. The return is an N x N Adjacency Array indexed with the bus numbers.

# Keyword arguments
- `check_connectivity::Bool`: Checks connectivity of the network using Goderya's algorithm
"""
function Adjacency(sys::System; check_connectivity::Bool = true)
    nodes = sort!(collect(get_components(Bus, sys)), by = x -> get_number(x))
    branches = get_components(Branch, sys, get_available)
    return Adjacency(branches, nodes; check_connectivity = check_connectivity)
end

function _goderya(ybus::SparseArrays.SparseMatrixCSC)
    node_count = size(ybus)[1]
    max_I = node_count^2
    I, J, val = SparseArrays.findnz(ybus)
    T = SparseArrays.sparse(I, J, ones(Int, length(val)))
    T_ = T * T
    for n in 1:(node_count - 1)
        I, _, _ = SparseArrays.findnz(T_)
        if length(I) == max_I
            @info "The System has no islands"
            break
        elseif length(I) < max_I
            temp = T_ * T
            I_temp, _, _ = SparseArrays.findnz(temp)
            if all(I_temp == I)
                @warn "The system contains islands" maxlog = 1
            end
            T_ = temp
        else
            @assert false
        end
        #@assert n < node_count - 1
    end
    return I
end

function validate_connectivity(sys::System)
    a = Adjacency(sys; check_connectivity = false)
    return validate_connectivity(a.data, collect(get_components(Bus, sys)), a.lookup[1])
end

function validate_connectivity(M, nodes, bus_lookup)
    I = _goderya(M)

    node_count = length(nodes)
    connections = Dict([i => count(x -> x == i, I) for i in Set(I)])

    if length(Set(I)) == node_count
        connected = true
        if any(values(connections) .!= node_count)
            cc = Set(values(connections))
            @warn "Network has at least $(length(cc)) connected components with $cc nodes"
            connected = false
        end
    else
        disconnected_nodes = get_name.(nodes[setdiff(values(bus_lookup), I)])
        @warn "Principal connected component does not contain:" disconnected_nodes
        connected = false
    end
    return connected
end


text:
"""
```
shock_decompositions(m, system, histshocks)

shock_decompositions(system, forecast_horizons, histshocks, start_index,
    end_index)
```

### Inputs

- `system::System{S}`: state-space system matrices
- `forecast_horizons::Int`: number of periods ahead to forecast
- `histshocks::Matrix{S}`: matrix of size `nshocks` x `hist_periods` of
  historical smoothed shocks
- `start_index::Int`: first index from which to return computed shock
  decompositions
- `end_index::Int`: last index for which to return computed shock decompositions

where `S<:AbstractFloat`.

### Outputs

- `states::Array{S, 3}`: matrix of size `nstates` x `nperiods` x `nshocks` of
  state shock decompositions
- `obs::Array{S, 3}`: matrix of size `nobs` x `nperiods` x `nshocks` of
  observable shock decompositions
- `pseudo::Array{S, 3}`: matrix of size `npseudo` x `nperiods` x `nshocks` of
  pseudo-observable shock decompositions

where `nperiods = `end_index - start_index + 1`.
"""
function shock_decompositions{S<:AbstractFloat}(m::AbstractModel,
    system::System{S}, histshocks::Matrix{S})

    horizon   = forecast_horizons(m)
    start_ind = index_shockdec_start(m)
    end_ind   = index_shockdec_end(m)

    shock_decompositions(system, horizon, histshocks, start_ind, end_ind)
end

function shock_decompositions{S<:AbstractFloat}(system::System{S},
    forecast_horizons::Int, histshocks::Matrix{S},
    start_index::Int, end_index::Int)

    # Setup
    nshocks      = size(system[:RRR], 2)
    nstates      = size(system[:TTT], 2)
    nobs         = size(system[:ZZ], 1)
    npseudo      = size(system[:ZZ_pseudo], 1)
    histperiods  = size(histshocks, 2)
    allperiods   = histperiods + forecast_horizons

    states = zeros(S, nstates, allperiods, nshocks)
    obs    = zeros(S, nobs, allperiods, nshocks)
    pseudo = zeros(S, npseudo, allperiods, nshocks)

    # Check dates
    if forecast_horizons <= 0 || start_index < 1 || end_index > allperiods
        throw(DomainError())
    end

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    z0 = zeros(S, nstates)

    for i = 1:nshocks
        # Isolate single shock
        shocks = zeros(S, nshocks, allperiods)
        shocks[i, 1:histperiods] = histshocks[i, :]

        # Iterate state space forward
        states[:, :, i], obs[:, :, i], pseudo[:, :, i], _ = forecast(system, z0, shocks)
    end

    # Return shock decompositions in appropriate range
    if start_index == 1 && end_index == allperiods
        return states, obs, pseudo
    else
        range = start_index:end_index
        return states[:, range, :], obs[:, range, :], pseudo[:, range, :]
    end
end

"""
```
deterministic_trends(m, system, z0)

deterministic_trends(system, z0, nperiods, start_index, end_index)
```

Compute deterministic trend values of states, observables, and
pseudo-observables, given a model object and system matrices. The deterministic
trend for a single draw is simply the series that would be obtained by iterating
the state-space system forward, beginning from a state vector `z0` in the last
presample period.

### Inputs

- `m::AbstractModel`: model object
- `system::System{S}`: state-space system matrices
- `z0`::Vector{S}: initial state vector

where `S<:AbstractFloat`.

### Outputs

- `states::Matrix{S}`: matrix of size `nstates` x `nperiods` of state
  steady-state values
- `obs::Matrix{S}`: matrix of size `nobs` x `nperiods` of observable
  steady-state values
- `pseudo::Matrix{S}`: matrix of size `npseudo` x `nperiods` of
  pseudo-observable steady-state values

where `nperiods` is the number of quarters between
`date_shockdec_start(m)` and `date_shockdec_end(m)`, inclusive.
"""
function deterministic_trends{S<:AbstractFloat}(m::AbstractModel{S},
    system::System{S}, z0::Vector{S})

    # Dates: We compute the deterministic trend starting from the
    # first historical period.  However, since it is only used to
    # compute shock decompositions, we truncate and only store
    # results for periods corresponding to the shockdec period.
    nperiods    = subtract_quarters(date_forecast_end(m), date_mainsample_start(m)) + 1
    start_index = index_shockdec_start(m)
    end_index   = index_shockdec_end(m)

    deterministic_trends(system, z0, nperiods, start_index, end_index)
end

function deterministic_trends{S<:AbstractFloat}(system::System{S}, z0::Vector{S}, nperiods::Int,
    start_index::Int, end_index::Int)

    # Set constant system matrices to 0
    system = zero_system_constants(system)

    # Construct matrix of 0 shocks for entire history and forecast horizon
    nshocks  = size(system[:RRR], 2)
    shocks   = zeros(S, nshocks, nperiods)

    # Use forecast to iterate state-space system forward without shocks or the ZLB procedure
    states, obs, pseudo, _ = forecast(system, z0, shocks)

    # Return deterministic trends in appropriate range
    if start_index == 1 && end_index == nperiods
        return states, obs, pseudo
    else
        range = start_index:end_index
        return states[:, range], obs[:, range], pseudo[:, range]
    end
end

"""
```
trends{S<:AbstractFloat}(system::System{S})
```

Compute trend (steady-state) states, observables, and pseudo-observables. The
trend is used for plotting shock decompositions.
"""
function trends{S<:AbstractFloat}(system::System{S})

    state_trend  = system[:CCC]
    obs_trend    = system[:ZZ]*system[:CCC] + system[:DD]
    pseudo_trend = system[:ZZ_pseudo]*system[:CCC] + system[:DD_pseudo]

    return state_trend, obs_trend, pseudo_trend
end


text:
module SantasLittleHelpers

using Statistics, Dates, LinearAlgebra
using Missings
using StaticKernels
using ShiftedArrays: lag
using DataFrames

export skew, kurt, autocorrelate, varianceratio, vr2, vr4, vr6, vr8, vr12, loggrowth
export isnanmissing, nantomissing!, showallrows
export applyrolling, makekernel, logdiff
export yrqtrfun, simiterator, perioditerator, monthtoquarter
export conditionalcor, downsidecor, upsidecor

include("rollingkernels.jl")
include("statistics.jl")
include("labelling.jl")
include("datautils.jl")

end # module


text:
estimates = collect(KalmanFilter(lssm, lobs), measures)
mu = [e[1] for e in estimates]
emp_var = empirical_var(mu-states,n,k)
@test emp_var < n*2.


text:
include("ops.jl")
using Test

mpi_init()
A = Float64[1 0 2 0 3
     0 0 1 2 3
     0 3 4 1 0
     0 4 2 1 0
     1 1 2 0 1]
A = A * A'
B = SparseTensor(A)
rhs = rand(5)
ilower = 0
iupper = 4
solver = "GMRES"
printlevel = 2
rows, ncols, cols, out =  mpi_create_matrix(B.o.indices,B.o.values,ilower,iupper)
# TODO: specify your input parameters

# error()
u = mpi_tensor_solve(rows,ncols,cols,out,rhs,ilower,iupper,solver,printlevel)
sess = Session(); init(sess)
u_out = run(sess, u)
u_ref = A\rhs

@show u_out - u_ref

# uncomment it for testing gradients
error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    return sum(mpi_tensor_solve(rows,ncols,cols,values,rhs,ilower,iupper,solver,printlevel)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(10,20))
v_ = rand(10,20)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
using Distributed

@everywhere module CoordTrans

    using NCDatasets
    using SharedArrays
    using Formatting
    using Distributed

    mutable struct WeightInfo

        NNN_max    :: Int64

        s_N        :: Int64
        d_N        :: Int64

        NN_idx     :: AbstractArray{Int64, 2}
        s_gc_lat   :: AbstractArray{Float64, 1}
        s_gc_lon   :: AbstractArray{Float64, 1}
        d_gc_lat   :: AbstractArray{Float64, 1}
        d_gc_lon   :: AbstractArray{Float64, 1}

        s_dims     :: AbstractArray{Int64, 1}
        d_dims     :: AbstractArray{Int64, 1}

        s_wgt      :: AbstractArray{Float64, 1}

    end

    mutable struct GridInfo

        N      :: Int64

        gc_lon :: AbstractArray{Float64}
        gc_lat :: AbstractArray{Float64}

        area   :: AbstractArray{Float64}
        mask   :: AbstractArray{Float64}

        unit_of_angle :: Symbol

        dims   :: AbstractArray{Int64}

        function GridInfo(;
            gc_lon :: AbstractArray{Float64,1},
            gc_lat :: AbstractArray{Float64,1},
            area   :: AbstractArray{Float64,1},
            mask   :: AbstractArray{Float64,1},
            unit_of_angle :: Symbol,
            dims   = nothing,
        )

            N = length(gc_lon)

            for var in [gc_lat, area, mask]
                if length(var) != N
                    throw(ErrorException("Not all input has the same length."))
                end
            end

            if unit_of_angle == :deg

                gc_lon .*= π / 180.0
                gc_lat .*= π / 180.0

            elseif unit_of_angle == :rad
                # do nothing

            else
                throw(ErrorException("`unit_of_angle` must be `:deg` or `:rad`."))
            end

            if dims == nothing
                dims = (N,)
            end

            dims = convert(Array{Int64}, dims)

            if reduce(*, dims) != N
                throw(ErrorException("Dims does not match the number of elements."))
            end

            return new(N, gc_lon, gc_lat, area, mask, unit_of_angle, dims)
        end

    end

    function genWeight_NearestNeighbors(
        filename :: AbstractString,
        gi_s     :: GridInfo,
        gi_d     :: GridInfo,
        NNN_max  :: Integer;
    )

        if length(gi_s.dims) != length(gi_d.dims)
            throw(ErrorException("Source and destination grid does not have same number dimensions"))
        end

        trans = SharedArray{Float64}((NNN_max, gi_d.N))

        # s_coord and d_coord are the coordinates of grid points
        # in 3-dimensional cartesian coordinate

        s_coord = SharedArray{Float64}(3, gi_s.N)
        d_coord = SharedArray{Float64}(3, gi_d.N)

        s_NaN_idx = (gi_s.mask .== 0)

        @sync @distributed for i = 1:gi_s.N

            s_coord[1, i] = cos(gi_s.gc_lat[i]) * cos(gi_s.gc_lon[i])
            s_coord[2, i] = cos(gi_s.gc_lat[i]) * sin(gi_s.gc_lon[i])
            s_coord[3, i] = sin(gi_s.gc_lat[i])

        end

        @sync @distributed for i = 1:gi_d.N

            d_coord[1, i] = cos(gi_d.gc_lat[i]) * cos(gi_d.gc_lon[i])
            d_coord[2, i] = cos(gi_d.gc_lat[i]) * sin(gi_d.gc_lon[i])
            d_coord[3, i] = sin(gi_d.gc_lat[i])

        end

        #s_NaN_idx = (s_mask .== 0)

        println("Start making transform matrix... ")

        @time @sync @distributed for i = 1:gi_d.N

            # For every point find its nearest-neighbors

            #print("\r", i, "/", d_N)

            if gi_d.mask[i] == 0
                trans[:, i] .= 0
                continue
            end

            dist2 = (  (s_coord[1, :] .- d_coord[1, i]).^2
                     + (s_coord[2, :] .- d_coord[2, i]).^2
                     + (s_coord[3, :] .- d_coord[3, i]).^2 )

            # Decided not to apply this condition because in
            # extreme cases there might be a small area of water
            # that is surrounded by lands.

            dist2[s_NaN_idx] .= NaN

            idx_arr = collect(1:gi_s.N)
            sort!(idx_arr; by=(k)->dist2[k])
            trans[:, i] = idx_arr[1:NNN_max]

        end
        if any(isnan.(trans))
            throw(ErrorException("Weird!"))
        end
        trans = convert(Array{Int64}, trans)
        println(typeof(gi_s.dims))

        wi = WeightInfo(
            NNN_max,
            gi_s.N,
            gi_d.N,
            trans,
            gi_s.gc_lat,
            gi_s.gc_lon,
            gi_d.gc_lat,
            gi_d.gc_lon,
            gi_s.dims,
            gi_d.dims,
            gi_s.area,
        )

        writeWeightInfo(wi, filename)
    end

    function writeWeightInfo(
        wi::WeightInfo,
        filename :: AbstractString;
        missing_value = 1e20,
    )

        Dataset(filename, "c") do ds

            defDim(ds, "s_N", wi.s_N)
            defDim(ds, "d_N", wi.d_N)
            defDim(ds, "NNN_max", size(wi.NN_idx)[1])
            defDim(ds, "dims", length(wi.s_dims))

            for (varname, vardata, vardims) in (
                ("NN_idx",    wi.NN_idx, ("NNN_max", "d_N")),
                ("s_gc_lat",  wi.s_gc_lat, ("s_N",)),
                ("s_gc_lon",  wi.s_gc_lon, ("s_N",)),
                ("d_gc_lat",  wi.d_gc_lat, ("d_N",)),
                ("d_gc_lon",  wi.d_gc_lon, ("d_N",)),
                ("s_dims",    wi.s_dims, ("dims",)),
                ("d_dims",    wi.d_dims, ("dims",)),
                ("s_wgt",     wi.s_wgt, ("s_N",)),
            )

                print(format("Output data: {} ...", varname))

                dtype = eltype(vardata)

                v = defVar(ds, varname, eltype(vardata), vardims)

                if dtype <: AbstractFloat
                    v.attrib["_FillValue"] = missing_value
                end

                v[:] = vardata
                println("done.")
            end

        end

    end

    function readWeightInfo(
        filename :: AbstractString
    )
        local wi

        Dataset(filename, "r") do ds

            wi = WeightInfo(
                ds.dim["NNN_max"],
                ds.dim["s_N"],
                ds.dim["d_N"],
                replace(ds["NN_idx"], missing=>0),
                replace(ds["s_gc_lat"][:], missing=>NaN),
                replace(ds["d_gc_lon"][:], missing=>NaN),
                replace(ds["s_gc_lat"][:], missing=>NaN),
                replace(ds["d_gc_lon"][:], missing=>NaN),
                replace(ds["s_dims"][:], missing=>0),
                replace(ds["d_dims"][:], missing=>0),
                replace(ds["s_wgt"][:], missing=>NaN),
            )

        end

        return wi

    end

    function convertData!(
        wi      :: WeightInfo,
        s_data  :: AbstractArray{Float64, 1},
        d_data  :: AbstractArray{Float64, 1},
    )

        NN_idx = wi.NN_idx
        s_wgt  = wi.s_wgt

        NNN = size(NN_idx)[1]
        for i = 1 : length(d_data)

            if NN_idx[1, i] == 0
                d_data[i] = NaN
                continue
            end

            wgt_sum = 0.0
            d_data[i] = 0.0

            for j = 1:NNN

                idx = NN_idx[j, i]
                data = s_data[idx]

                if isfinite(data)
                    wgt_sum += s_wgt[idx]
                    d_data[i] += data * s_wgt[idx]
                else
                    break
                end
            end

            d_data[i] = (wgt_sum == 0) ? NaN : d_data[i] / wgt_sum

        end

    end

    function convertFile(
        in_filename   :: AbstractString,
        out_filename  :: AbstractString,
        wgt_filename  :: Union{AbstractString, WeightInfo};
        varnames      :: Union{Nothing, Array} = nothing,
        xydim         :: AbstractString = "grid",
        xdim          :: AbstractString = "lon",
        ydim          :: AbstractString = "lat",
        zdim          :: Union{AbstractString, Nothing} = "lev",
        tdim          :: AbstractString = "time",
        xydim_val     :: Union{AbstractArray{Float64, 1}, Nothing} = nothing,
        xdim_val      :: Union{AbstractArray{Float64, 1}, Nothing} = nothing,
        ydim_val      :: Union{AbstractArray{Float64, 1}, Nothing} = nothing,
        zdim_val      :: Union{AbstractArray{Float64, 1}, Nothing} = nothing,
    )

        if typeof(wgt_filename) <: AbstractString
            wi = readWeightInfo(wgt_filename)
        else
            wi = wgt_filename
        end

        dim_len = length(wi.d_dims)
        d_data_tmp = zeros(Float64, wi.d_N)

        ds_in  = Dataset(in_filename, "r")
        ds_out = Dataset(out_filename, "c")

        if ! (dim_len in (1, 2))
            throw(ErrorException("Weird"))
        end

        # Defining dimension
        if dim_len == 1
            defDim(ds_out, xydim, wi.d_dims[1])
        elseif dim_len == 2
            defDim(ds_out, xdim, wi.d_dims[1])
            defDim(ds_out, ydim, wi.d_dims[2])
        end

        defDim(ds_out, tdim, Inf)

        if zdim != nothing && zdim in ds_in.dim
            defDim(ds_out, zdim, ds_in.dim[zdim])
        end

        # Fill in dimension values

        if dim_len == 1

            if xydim_val != nothing
                defVar(ds_out, xydim, xydim_val, (xydim,))
            end

        elseif dim_len == 2

            if xdim_val != nothing
                defVar(ds_out, xdim, xdim_val, (xdim,))
            end

            if ydim_val != nothing
                defVar(ds_out, ydim, ydim_val, (ydim,))
            end
        end

        if zdim_val != nothing
            defVar(ds_out, zdim, zdim_val, (zdim,))
        elseif zdim in ds_in
            defVar(ds_out, zdim, replace(ds_in[zdim][:], missing=>0), (zdim,))
        end

        if varnames == nothing

            varnames = Array{Any}(undef, 0)

            for varname in keys(ds_in)
                println("varname:", varname, "; Dimnames: ", dimnames(ds_in[varname]))
                if dim_len == 1
                    if dimnames(ds_in[varname])[1] == xydim
                        push!(varnames, varname)
                    end
                elseif dim_len == 2
                    if dimnames(ds_in[varname])[1:2] == (xdim, ydim,)
                        push!(varnames, varname)
                    end
                end
            end
        end

        println("Defined dimensions: ", keys(ds_out.dim))

        # Converting variables
        for varname in varnames

            println("Dealing with varname: ", varname)

            if varname in ds_out
                println(format("Varname: {} already exists. Skip.", varname))
                continue
            end

            if ! (varname in ds_in)
                println(format("Cannot find var: {}. Skip.", varname))
                continue
            end

            cf_var = ds_in[varname]
            cf_var_dimnames = dimnames(cf_var)

            # Cannot read data yet. Some data are really large
            #s_data = replace(cf_var[:], missing=>NaN)

            s_data_dims = size(cf_var)
            s_data_dims_len = length(s_data_dims)

            d1 = reduce(*, s_data_dims[1:dim_len])
            d2 = (s_data_dims_len > dim_len) ? reduce(*, s_data_dims[dim_len+1:end]) : 1

            #s_data = reshape(s_data, d1, d2)

            attrib = Dict()
            for (k,v) in cf_var.attrib
                if typeof(v) <: AbstractFloat
                    attrib[k] = convert(Float64, v)
                else
                    attrib[k] = v
                end
            end

            println(attrib)
            println(cf_var_dimnames)
            v = defVar(ds_out, varname, Float64, cf_var_dimnames, attrib=attrib)

            for k = 1:d2

                # Construct reading / writing shape
                idx = Array{Any}(undef, 0)
                for i=1:dim_len
                    push!(idx, :)
                end

                if s_data_dims_len == dim_len + 1
                    push!(idx, k)
                elseif s_data_dims_len == dim_len + 2
                    zidx = mod(k-1, ds_in.dim[zdim]) + 1
                    tidx = floor(Integer, (k - 1) / ds_in.dim[zdim]) + 1
                    push!(idx, zidx)
                    push!(idx, tidx)
                end
                #println(idx)

                convertData!(wi, reshape( replace( cf_var[idx...], missing=>NaN ), d1 ), d_data_tmp)
                v[idx...] = reshape(d_data_tmp, wi.d_dims[1], wi.d_dims[2])
            end
        end

        close(ds_in)
        close(ds_out)
    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# continuously drive in a circle to demonstrate fixed lag

using RoME

fg = initfg()
getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = true
getSolverParams(fg).dbg = true

defaultFixedLagOnTree!(fg, 5, limitfixeddown=true)

tree = BayesTree()

generateCanonicalFG_Circle(6, fg=fg, offsetPoses=0, stopEarly=3, cyclePoses=6)
tree = solveTree!(fg, tree, recordcliqs=ls(fg));
drawTree(tree,filepath=joinLogPath(fg,"tree3.pdf"))

generateCanonicalFG_Circle(6, fg=fg, offsetPoses=3, stopEarly=6, cyclePoses=6)
tree = solveTree!(fg, tree);
drawTree(tree,filepath=joinLogPath(fg,"tree6.pdf"))

generateCanonicalFG_Circle(12, fg=fg, offsetPoses=6, stopEarly=9, cyclePoses=6)
tree = solveTree!(fg, tree);
drawTree(tree,filepath=joinLogPath(fg,"tree9.pdf"))

generateCanonicalFG_Circle(12, fg=fg, offsetPoses=9, stopEarly=12, cyclePoses=6)
tree = solveTree!(fg, tree);
drawTree(tree,filepath=joinLogPath(fg,"tree12.pdf"))

generateCanonicalFG_Circle(18, fg=fg, offsetPoses=12, stopEarly=15, cyclePoses=6)
tree = solveTree!(fg, tree);
drawTree(tree,filepath=joinLogPath(fg,"tree15.pdf"))

generateCanonicalFG_Circle(18, fg=fg, offsetPoses=15, stopEarly=18, cyclePoses=6)
tree = solveTree!(fg, tree);
drawTree(tree,filepath=joinLogPath(fg,"tree18.pdf"))

drawGraph(fg, show=true)

## visualization of numerical results

using RoMEPlotting

Gadfly.set_default_plot_size(35cm,20cm)
plotSLAM2D(fg)


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: https://julialang.org/license

prefix = joinpath(homedir(), ".julia", "registries", "Uncurated")

write_toml(prefix, "Registry") do io
    repo = "https://github.com/JuliaRegistries/Uncurated.git"
    uuid = string(uuid5(uuid_registry, repo))
    println(io, "name = ", repr("Uncurated"))
    println(io, "uuid = ", repr(uuid))
    println(io, "repo = ", repr(repo))
    println(io, "\ndescription = \"\"\"")
    print(io, """
        Official uncurated Julia package registry where people can
        register any package they want without too much debate about
        naming and without enforced standards on documentation or
        testing. We nevertheless encourage documentation, testing and
        some amount of consideration when choosing package names.
        """)
    println(io, "\"\"\"")
    println(io, "\n[packages]")
    for (pkg, p) in sort!(collect(pkgs), by=(p->p.uuid.value)∘last)
        bucket = string(uppercase(first(pkg)))
        path = joinpath(bucket, pkg)
        println(io, p.uuid, " = { name = ", repr(pkg), ", path = ", repr(path), " }")
    end
end

buckets = Dict()
for (pkg, p) in pkgs
    bucket = string(uppercase(first(pkg)))
    push!(get!(buckets, bucket, []), (pkg, p))
end

const trees, stdlibs = gitmeta(pkgs)

for pkg in STDLIBS
    tree = stdlib_trees[pkg]
    deps = Dict(dep => Require(VersionInterval()) for dep in stdlib_deps[pkg])
    pkgs[pkg] = Package(
        UUID(stdlib_uuids[pkg]),
        "https://github.com/JuliaLang/julia.git",
        Dict(VersionNumber(0,7,0,("DEV",),("r"*tree[1:8],)) => Version(tree, deps)),
    )
end

for (pkg, p) in pkgs
    uuid = string(p.uuid)
    haskey(stdlibs, uuid) || continue
    for (ver, v) in p.versions
        n = get(stdlibs[uuid], v.sha1, 0)
        n == 0 && continue
        for lib in STDLIBS
            if n & 1 != 0
                v.requires[lib] = Require(VersionInterval())
            end
            n >>>= 1
        end
    end
end

for (bucket, b_pkgs) in buckets, (pkg, p) in b_pkgs
    haskey(stdlibs, pkg) && continue
    url = p.url
    uuid = string(p.uuid)
    startswith(url, "git://github.com") && (url = "https"*url[4:end])

    # Package.toml
    write_toml(prefix, bucket, pkg, "Package") do io
        println(io, "name = ", repr(pkg))
        println(io, "uuid = ", repr(uuid))
        println(io, "repo = ", repr(url))
    end

    # Versions.toml
    write_toml(prefix, bucket, pkg, "Versions") do io
        for (i, (ver, v)) in enumerate(sort!(collect(p.versions), by=first))
            i > 1 && println(io)
            println(io, "[", toml_key(string(ver)), "]")
            println(io, "git-tree-sha1 = ", repr(trees[uuid][v.sha1]))
        end
    end
    versions = sort!(collect(keys(p.versions)))

    function write_versions_data(f::Function, name::String; lt::Function=isless)
        data = Dict{VersionNumber,Dict{String,String}}()
        for (ver, v) in p.versions, (dep, d) in v.requires
            val = f(dep, d)
            val == nothing && continue
            haskey(data, ver) || (data[ver] = Dict{String,String}())
            # BinDeps injects a dependency on Libdl
            if name == "Deps" && dep == "BinDeps"
                data[ver]["Libdl"] = "\"8f399da3-3557-5675-b5ff-fb832c97cbdb\""
            end
            data[ver][dep] = val
        end
        compressed = compress_versions_data(data, versions)
        !isempty(compressed) && write_toml(prefix, bucket, pkg, name) do io
            vers = unique(getindex.(compressed, 1))
            keys = sort!(unique(getindex.(compressed, 2)), lt=lt)
            what = (vers, keys)
            ord = (1, 2)
            for (i, x) in enumerate(what[ord[1]])
                i > 1 && println(io)
                println(io, "[", toml_key(x), "]")
                for y in what[ord[2]]
                    for t in compressed
                        t[ord[1]] == x && t[ord[2]] == y || continue
                        println(io, toml_key(y), " = ", t[3])
                    end
                end
            end
        end
    end

    # Deps.toml
    write_versions_data("Deps") do dep, d
        dep == "julia" ? nothing : repr(string(pkgs[dep].uuid))
    end

    # Compat.toml
    write_versions_data("Compat", lt=packagelt) do dep, d
        dep in STDLIBS ? nothing : versions_repr(compress_versions(
            d.versions, collect(keys(pkgs[dep].versions))
        ))
    end
end

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
module Jevo

include("fitness.jl")
include("populations.jl")
include("initiate.jl")
include("wright_fisher.jl")
include("mutations.jl")
include("utils.jl")
end # module


text:
"""
Transient vibration analysis of free-floating thin plate

A homogeneous free-floating(unsupported) square plate from the NAFEMS
Benchmark, Test No. FV12.

The plate is discretized with shell elements. Because no displacements are
prevented, the structure has six rigid body modes (six zero vibration
frequencies).

"""

using LinearAlgebra
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using Gnuplot

function solve()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.05*phun("m");
    L = 10.0*phun("m");

    n = 4*16
    @info "Mesh: $n elements per side"

    tolerance = L/n/1000
    fens, fes = T3block(L,L,n,n);
    fens.xyz[:, 1] .-= L/2
    fens.xyz[:, 2] .-= L/2
    fens.xyz = xyz3(fens)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), mater)

# Construct the requisite fields, geometry and displacement
# Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

# Apply EBC's
    l1 = collect(1:count(fens))

    applyebc!(dchi)
    numberdofs!(dchi);

# Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    K = FEMMShellT3FFModule.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, geom0, dchi);
# @show sum(sum(M, dims = 1))/3,.

# Solve
    evals, evecs, nconv = eigs(K, M; nev=1, which=:LM, explicittransform=:none)
    omega_max = sqrt(evals[1])

    dt = 0.99* 2/omega_max

    U0 = gathersysvec(dchi)
    U1 = deepcopy(U0)
    V0 = deepcopy(U0)
    V1 = deepcopy(U0)
    A0 = deepcopy(U0)
    A1 = deepcopy(U0)
    F0 = deepcopy(U0)

    qpoint = selectnode(fens; nearestto=[L/4 L/4 0])[1]
    cpoint = selectnode(fens; nearestto=[0 0 0])[1]
    applyebc!(dchi)
    numberdofs!(dchi);
    qpointdof = dchi.dofnums[qpoint, 3]
    cpointdof = dchi.dofnums[cpoint, 3]

    V0[qpointdof] = 100.0

    A0 .= M\F0;

    nsteps = 100000
    cdeflections = fill(0.0, nsteps+1)
    displacements = []
    nbtw = 1000
    t = 0.0
    for step in 1:nsteps
        cdeflections[step] = U0[cpointdof]
        # displacement update
        @. U1 = U0 + dt*V0 + (dt^2)/2*A0;
# Compute the new acceleration.
        A1 .= M\(F0-K*U1);
# Update the velocity
        @. V1 = V0 + (dt/2)* (A0+A1);
        t = t + dt
        @. U0 = U1
        @. V0 = V1
        @. A0 = A1
        if rem(step, nbtw) == 0
            push!(displacements, deepcopy(U1))
        end
    end
    cdeflections[end] = U0[cpointdof]

    @gp  "set terminal windows 0 "  :-

    @gp  :- collect(0.0:dt:(nsteps*dt)) cdeflections " lw 2 lc rgb 'red' with p title 'Deflection at the center' "  :-

    @gp  :- "set xlabel 'Time'" :-
    @gp  :- "set ylabel 'Deflection'" :-
    @gp  :- "set title 'Free-floating plate'"

# Visualization
    @show length(displacements)
    tbox = plot_space_box(reshape(inflatebox!(boundingbox(fens.xyz), 0.1 * L), 2, 3))
    tenv0 = plot_midsurface(fens, fes; x=geom0.values, u=0.0 .* dchi.values[:, 1:3], R=Rfield0.values, facecolor="rgb(125, 155, 125)", opacity=0.3);
    plots = cat(tbox, tenv0; dims=1)
    layout = default_layout_3d(;width=600, height=600)
    layout[:scene][:aspectmode] = "data"
    pl = render(plots; layout=layout, title = "Step ")
    sleep(0.5)

    scale = 500;
    for i in 1:length(displacements)
        scattersysvec!(dchi, scale .* displacements[i])
        update_rotation_field!(Rfield0, dchi)
        tenv1 = plot_midsurface(fens, fes; x=geom0.values, u=dchi.values[:, 1:3], R=Rfield1.values, facecolor="rgb(50, 155, 225)");
        plots = cat(tbox, tenv0, tenv1; dims=1)
        pl.plot.layout[:title] = "Step $(i)"
        react!(pl, plots, pl.plot.layout)
        sleep(0.115)
    end

end

solve()

nothing


text:
using LinearAlgebra
using PyPlot

function gauss_quad(p)
    n = ceil((p+1)/2)
    b = 1:n-1
    b = @. b / sqrt(4*b^2 - 1)
    eval, evec = eigen(diagm(1 => b, -1 => b))
    return eval, 2*evec[1,:].^2
end

function legendre_poly(x, p)
    z = zeros(size(x))
    o = ones(size(x))
    y = hcat(o, x, repeat(z, 1, p-1))
    dy = hcat(z, o, repeat(z, 1, p-1))
    for i = 1:p-1
        @. y[:,i+2] = ((2i+1)*x*y[:,i+1] - i*y[:,i]) / (i+1)
        @. dy[:,i+2] = ((2i+1)*(x*dy[:,i+1] + y[:,i+1]) - i*dy[:,i]) / (i+1)
    end
    y, dy
end

function dgconvect(; n=10, p=1, T=1.0, dt=1e-3)
    # Discretization
    h = 1 / n
    s = @.cos((pi*collect(0:p))/p)
    s1 = 0.5*(-s.+1)
    x = @. h*s1 + (0:h:1-h)'

    # Gaussian initial condition (and exact solution if shifted)
    uinit(x) = exp(-(x - 0.5)^2 / 0.1^2)

    # Coefficient Matrix
    lhsc = legendre_poly(-s,p)[1]
    cmx = lhsc \ Matrix{Float64}(I,p+1,p+1)

    # Mass and Stiffness Matrices

    Mel = zeros(p+1,p+1)
    Kel = zeros(p+1,p+1)
    gx, gw = gauss_quad(2p)
    #gx = (h/2).*(gx.+1)
    for i = 0:p, j = 0:p
	    mdvec = zeros(p+1)
	    kdvec = zeros(p+1)
	    for k = 0:p
	         lpval, dlpval = legendre_poly([gx[k+1]], p)
		 mdvec[k+1] = dot(lpval,cmx[:,i+1]) * dot(lpval,cmx[:,j+1])
		 kdvec[k+1] = dot(dlpval,cmx[:,i+1]) * dot(lpval,cmx[:,j+1])
	    end
	    Mel[i+1,j+1] = dot(mdvec, gw)
	    Kel[i+1,j+1] = dot(kdvec, gw)
    end

    Mel = Mel.*(h/2)
    Kel = Kel

    # RHS function in semi-discrete system
    function rhs(u)
        r = Kel * u
        r[end,:] = r[end,:] - u[end,:]
        r[1,:] = r[1,:] + u[end, [end; 1:end-1]]
        r = Mel \ r
    end

    # Setup
    u = uinit.(x)
    nsteps = round(Int, T/dt)

    # Setup plotting
    xx = collect(0:0.01:1) # For exact solution
    clf(); axis([0, 1, -0.1, 1.1]); grid(true);
    ph1 = plot(x, u, "b")   # Many line object (since discontinuous)
    ph2, = plot(xx, uinit.(xx), "k") # Single line object

    # Main loop
    for it = 1:nsteps
        # Runge-Kutta 4
        k1 = dt * rhs(u)
        k2 = dt * rhs(u + k1/2)
        k3 = dt * rhs(u + k2/2)
        k4 = dt * rhs(u + k3)
        u += (k1 + 2*k2 + 2*k3 + k4) / 6

        # Plotting
        if mod(it, round(nsteps/100)) == 0
            uexact = @. uinit(mod(xx - dt*it, 1.0))
            for line in ph1 line.remove(); end
            ph1 = plot(x, u, "b")
            ph2.set_data(xx, uexact)
            pause(1e-3)
        end
    end
    ua = uinit.(x)
    xw = zeros(size(x)[1],size(x)[2])
    for j = 1:p
	xw[p-j+2,:] = x[p-j+2,:] - x[p-j+1,:]
    end
    eu = (sum((u.^2).*xw))^0.5
    eua	= (sum((ua.^2).*xw))^0.5
    return u, abs(eua-eu)
end

function dgconvect_convergence()
    errors = zeros(5,5)
    for i = 0:4, j = 4:8
        p1 = ceil(Int,2^i)
        n1 = ceil(Int,(2^j)/p1)
        println(p1, " ", n1)
	errors[j-3,i+1] = dgconvect( n=n1, p=p1, dt=2e-4)[2]
    end
    clf()
    loglog([16 32 64 128 256]', errors)
    rates = @. log2(errors[end-1,:]) - log2(errors[end,:])
    return errors, rates
end

function dgconvdiff(; n=10, p=1, T=1.0, dt=1e-3, k=1e-3)
    # Discretization
    h = 1 / n
    s = @.cos((pi*collect(0:p))/p)
    s1 = 0.5*(-s.+1)
    x = @. h*s1 + (0:h:1-h)'

    # Gaussian initial condition (and exact solution if shifted)
    uinit(x) = exp(-(x - 0.5)^2 / 0.1^2)

    # Coefficient Matrix
    lhsc = legendre_poly(-s,p)[1]
    cmx = lhsc \ Matrix{Float64}(I,p+1,p+1)

    # Mass and Stiffness Matrices

    Mel = zeros(p+1,p+1)
    Kel = zeros(p+1,p+1)
    gx, gw = gauss_quad(2p)
    #gx = (h/2).*(gx.+1)
    for i = 0:p, j = 0:p
	    mdvec = zeros(p+1)
	    kdvec = zeros(p+1)
	    for k = 0:p
	         lpval, dlpval = legendre_poly([gx[k+1]], p)
		 mdvec[k+1] = dot(lpval,cmx[:,i+1]) * dot(lpval,cmx[:,j+1])
		 kdvec[k+1] = dot(dlpval,cmx[:,i+1]) * dot(lpval,cmx[:,j+1])
	    end
	    Mel[i+1,j+1] = dot(mdvec, gw)
	    Kel[i+1,j+1] = dot(kdvec, gw)
    end

    Mel = Mel.*(h/2)
    Kel = Kel

    u = uinit.(x)
    us = uinit.(x)

    # RHS function in semi-discrete system
    function rhs(u)
        up = vcat(u[end,:]' , u[1:end-1,:])
        un = vcat(u[2:end,:] , u[1,:]')
        xrev = abs.(x.-vcat(x[2,:]' , x[1:end-1,:]))
        u2 = (2u - un - up)./(2xrev)

        r = -Kel * (u)
        r[end,:] = r[end,:] + u[end,:]
        r[1,:] = r[1,:] - u[end, [end; 1:end-1]]
        r = Mel \ r
        sigma = r
        r = Kel * (u - k*(r))
        r[end,:] = r[end,:] - u[end,:] + k*sigma[end,:]
        r[1,:] = r[1,:] + u[end, [end; 1:end-1]] - k*sigma[end, [end; 1:end-1]]
        r = Mel \ r
    end

    # Setup

    nsteps = round(Int, T/dt)

    # Setup plotting
    xx = collect(0:0.01:1) # For exact solution
    uex(x,t,i,k) = exp(-100((x-t-0.5+i)^2/(1+400*k*t)))/sqrt(1+400*k*t)
    clf(); axis([0, 1, -0.1, 1.1]); grid(true);
    ph1 = plot(x, u, "b")   # Many line object (since discontinuous)
    ph2, = plot(xx, uinit.(xx), "k") # Single line object

    # Main loop
    for it = 1:nsteps
        # Runge-Kutta 4
        k1 = dt * rhs(u)
        k2 = dt * rhs(u + k1/2)
        k3 = dt * rhs(u + k2/2)
        k4 = dt * rhs(u + k3)
        u += (k1 + 2*k2 + 2*k3 + k4) / 6

        # Plotting
        if mod(it, round(nsteps/100)) == 0
            uexact = @. uex(xx,dt*it,0,k)+uex(xx,dt*it,1,k)+uex(xx,dt*it,-1,k)
            for line in ph1 line.remove(); end
            ph1 = plot(x, u, "b")
            ph2.set_data(xx, uexact)
            pause(1e-3)
        end
    end
    ua = @. uex(x,1,0,k)+uex(x,1,1,k) + uex(x,1,-1,k)
    xw = zeros(size(x)[1],size(x)[2])
    for j = 1:p
	xw[p-j+2,:] = x[p-j+2,:] - x[p-j+1,:]
    end
    eu = (sum((u.^2).*xw))^0.5
    eua	= (sum((ua.^2).*xw))^0.5
    return u, abs(eua-eu)
end

function dgconvdiff_convergence()
    errors = zeros(4,5)
    for i = 0:4, j = 4:7
        p1 = ceil(Int,2^i)
        n1 = ceil(Int,(2^j)/p1)
        println(p1, " ", n1)
	errors[j-3,i+1] = dgconvdiff( n=n1, p=p1, dt=1e-3, k = 1e-3)[2]
    end
    clf()
    loglog([16 32 64 128]', errors)
    rates = @. log2(errors[end-1,:]) - log2(errors[end,:])
    return errors, rates
end


text:
function particalsort!(A, b)
    l = 1
    r = length(A)
    while l < r
        while A[r] > b
            r -= 1
        end

        while A[l] <= b
            l += 1
        end
        if r > l
            A[l], A[r] = A[r], A[l]
        end
    end
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 7bba9ce0-1efe-11ec-2557-b9c38d2fbc38
md"""
# Chapter 17
#### Nicetak
 $(import Dates; Dates.format(Dates.today(), Dates.DateFormat("U d, Y")))
"""

# ╔═╡ d7594aeb-a318-4cf1-9ee0-e22512aa4681
md"""
## 17.1
ハミルトンサイクル問題がNPに属することを示してください.
"""

# ╔═╡ 645f0723-2329-4f5b-ae46-f1709b86ac38
md"""
**解答**

ハミルトンサイクルの解の候補 $p$ が存在するとして, これがハミルトンサイクルであることを確かめるには

- すべての頂点を通ること
- 合法的な経路であること

を確かめる必要がある. 解候補が$|V|$より長いときは, 自明にハミルトンサイクルでないことが確かめられるので,　たかだか $|V|$ の長さの解候補を見れば良い. このとき, 解候補の順番にたどっていけば, すべての頂点を通ることは自明に確かめられる. 合法的な経路であることは, 解候補の示すある頂点に対して, 次の頂点への辺があることを確かめればよいので, 一つの頂点に対して最大で $|V|$ の長さを見ればよいので, $O(|V|^2)$.

"""

# ╔═╡ 7063faec-9747-443a-8339-dd792109ccb5
md"""
## 17.2
グラフ $G = (V, E)$ と正の整数 $k$ があたえられたときに, $G$ がサイズ $k$ 以上の完全グラフを部分グラフとして含むかどうかを判定する問題を**クリーク問題**とよびます. クリーク問題がNPに属することを示してください.
"""

# ╔═╡ 6bae3ae3-6488-43ae-9548-d9e6a74647ff
md"""
**解答**

解候補は部分グラフである必要があるため, 入力は多項式長である必要があり, これは多項式時間で判定できる.
今, 完全グラフの解候補 $G' = (V', E')$ が与えられているとする. これが部分グラフであることを確かめるためには, 各 $e \in E'$ が $E$ に含まれていることを確かめればよいので, $O(|E|^2)$.
また, これが完全グラフであることを確かめるためには, すべての頂点がすべての他頂点に対して辺を持つことを確かめればよい. 一つの頂点に対して最大で $|V| - 1$ の辺が $E'$ に含まれているか確かめればよいので, $O(|V|^2|E|)$.
"""

# ╔═╡ d7cf1c75-c8de-4f91-9eb4-2c43c0a6df4b
md"""
## 17.3
安定集合問題がNP完全であることを利用して, クリーク問題もNP完全であることを示してください.
"""

# ╔═╡ 67708d53-4285-4025-9ac3-aea1902cafe6
md"""
**解答**

グラフ $G$ の頂点の部分集合 $S$ が安定集合であることは, 補グラフ $\bar{G}$ において, $S$ が完全グラフであることと同値である. したがって, サイズ $k$ 以上の安定集合が存在することと, 補グラフ $\bar{G}$ においてサイズ $k$ 以上の完全グラフが存在することは同値である. 補グラフの構成は自明に多項式時間 $O(|V|^2)$ で可能なので, クリーク問題は安定集合問題に多項式還元される.
"""

# ╔═╡ Cell order:
# ╟─7bba9ce0-1efe-11ec-2557-b9c38d2fbc38
# ╟─d7594aeb-a318-4cf1-9ee0-e22512aa4681
# ╟─645f0723-2329-4f5b-ae46-f1709b86ac38
# ╟─7063faec-9747-443a-8339-dd792109ccb5
# ╟─6bae3ae3-6488-43ae-9548-d9e6a74647ff
# ╟─d7cf1c75-c8de-4f91-9eb4-2c43c0a6df4b
# ╟─67708d53-4285-4025-9ac3-aea1902cafe6


text:
export size, ztr, copy, statespace_recursion

"""
    size(model::StateSpaceModel)

Return the dimensions `n`, `p`, `m` and `r` of the `StateSpaceModel`
"""
function size(model::StateSpaceModel)
    return model.dim.n, model.dim.p, model.dim.m, model.dim.r
end

"""
    ztr(model::StateSpaceModel)

Return the state space model arrays `Z`, `T` and `R` of the `StateSpaceModel`
"""
function ztr(model::StateSpaceModel)
    return model.Z, model.T, model.R
end

"""
    statespace_recursion(model::StateSpaceModel{Typ}, N::Int, initial_a::Matrix{Typ}) where Typ

Runs the recursion of a state space model ignoring the y.
"""
function statespace_recursion(model::StateSpaceModel{Typ}, initial_a::Matrix{Typ}) where Typ

    # If model has any unknows is not possible to perform the recursion
    unknowns = Unknowns(model)
    if unknowns.n_unknowns != 0
        error("StateSpaceModel has unknown parameters.")
    end

    # Save the distributions H and Q
    dist_H = MvNormal(model.H)
    dist_Q = MvNormal(model.Q)

    p, m, n = size(model.Z)

    if size(initial_a, 2) != m
        error("intial_a must be a 1 by $m matrix.")
    end

    y = Matrix{Typ}(undef, n, p)
    α = Matrix{Typ}(undef, n, m)
    α[1, :] = initial_a

    for t = 1:n-1
        y[t, :]   = model.Z[:, :, t]*α[t, :] + model.d[t, :] + rand(dist_H)
        α[t+1, :] = model.T*α[t, :] + model.c[t, :] + model.R*rand(dist_Q)
    end
    y[n, :] = model.Z[:, :, n]*α[n, :] +  model.d[n, :] + rand(dist_H)

    return y, α
end

# Linar Algebra wrappers
function gram_in_time(mat::Array{T, 3}) where T
    gram_in_time = similar(mat)
    @inbounds @views for t = 1:size(gram_in_time, 3)
        gram_in_time[:, :, t] = gram(mat[:, :, t])
    end
    return gram_in_time
end

function gram(mat::AbstractArray{T}) where T
    if size(mat, 1) == 1
        gram_mat = Matrix{T}(undef, 1, 1)
        gram_mat[1, 1] = mat[1, 1]^2
        return gram_mat
    else
        return LinearAlgebra.BLAS.gemm('N', 'T', mat, mat) # mat*mat'
    end
end

"""
    check_steady_state(P_t1::Matrix{T}, P_t::Matrix{T}, tol::T) where T

Return `true` if steady state was attained with respect to tolerance `tol`, `false` otherwise.
The steady state is checked by the following equation `maximum(abs.((P_t1 - P_t)./P_t1)) < tol`.
"""
function check_steady_state(P_t1::Matrix{T}, P_t::Matrix{T}, tol::T) where T
    return maximum(abs.((P_t1 - P_t)./P_t1)) < tol ? true : false
end

function check_steady_state(P::AbstractArray{T}, t::Int, tol::T) where T
    @inbounds for j in axes(P, 2), i in axes(P, 1)
        if abs((P[i, j, t+1] - P[i, j, t])/P[i, j, t+1]) > tol
            return false
        end
    end
    return true
end

"""
    ensure_pos_sym!(M::Matrix{T}; ϵ::T = 1e-8) where T

Ensure that matrix `M` is positive and symmetric to avoid numerical errors when numbers are small by doing `(M + M')/2 + ϵ*I`
"""
function ensure_pos_sym!(M::AbstractArray{T}, t::Int; ϵ::T = T(1e-8)) where T
    @inbounds for j in axes(M, 2), i in 1:j
        if i == j
            M[i, i, t] = (M[i, i, t] + M[i, i, t])/2 + ϵ
        else
            M[i, j, t] = (M[i, j, t] + M[j, i, t])/2
            M[j, i, t] = M[i, j, t]
        end
    end
    return
end

function ensure_pos_sym!(M::AbstractArray{T}; ϵ::T = T(1e-8)) where T
    @inbounds for j in axes(M, 2), i in 1:j
        if i == j
            M[i, i] = (M[i, i] + M[i, i])/2 + ϵ
        else
            M[i, j] = (M[i, j] + M[j, i])/2
            M[j, i] = M[i, j]
        end
    end
    return
end

function sum_matrix!(mat_prin::AbstractArray{T}, mat_sum::AbstractMatrix{T}, t::Int, offset::Int) where T
    @inbounds for j in axes(mat_prin, 2), i in axes(mat_prin, 1)
        mat_prin[i, j, t + offset] = mat_prin[i, j, t + offset] + mat_sum[i, j]
    end
    return
end

function sum_matrix!(mat_prin::AbstractArray{T}, mat_sum::AbstractArray{T}, t::Int, offset::Int) where T
    @inbounds for j in axes(mat_prin, 2), i in axes(mat_prin, 1)
        mat_prin[i, j, t + offset] = mat_prin[i, j, t + offset] + mat_sum[i, j, t]
    end
    return
end

function invertF(F::Array{T, 3}, t::Int) where T
    return @inbounds @views size(F, 1) == 1 ? 1/(F[1, 1, t]) : inv(F[:, :, t])
end

function invertF(F::Vector{T}, t::Int) where T
    return 1/F[t]
end

function logdetF(F::Array{T, 3}, t::Int) where T
    return @inbounds @views logdet(F[:, :, t])
end

function find_missing_observations(y::Matrix{T}) where T
    missing_obs = Vector{Int}(undef, 0)
    for i in axes(y, 1), j in axes(y, 2)
        if isnan(y[i, j])
            # If is empty or if the index i is not already pushed
            (isempty(missing_obs) || (missing_obs[end] != i)) && push!(missing_obs, i)
        end
    end
    return missing_obs
end

function Base.show(io::IO, ss::StateSpace)
    println("An estimated state-space model with")
    println("filter_type = $(ss.filter_type),")
    println("optimization_method = $(typeof(ss.opt_method)).")
    return nothing
end

function Base.show(io::IO, model::StateSpaceModel)
    println("A $(model.mode) state-space model with")
    println("n = $(model.dim.n),")
    println("p = $(model.dim.p),")
    println("m = $(model.dim.m),")
    println("r = $(model.dim.r).")
    return nothing
end

function build_H(p::Int, T)
    H = fill(NaN, p, p)
    return T.(H)
end

function build_Q(r::Int, p::Int, T)
    Q = kron(ones(p, p), Matrix{T}(I, Int(r/p), Int(r/p)))
    Q[findall(isequal(1), Q)] .= NaN
    return Q
end

function ensure_is_matrix(y::Vector{T}) where T
    return y[:, :]
end
function ensure_is_matrix(y::Matrix{T}) where T
    return y
end

function build_ss_dim(y::Matrix{Typ}, Z::Array{Typ, 3}, T::Matrix{Typ}, R::Matrix{Typ}) where Typ
    ny, py = size(y)
    pz, mz, nz = size(Z)
    mt1, mt2 = size(T)
    mr, rr = size(R)
    if !((mz == mt1 == mt2 == mr) && (pz == py))
        error("StateSpaceModel dimension mismatch")
    end
    return StateSpaceDimensions(ny, py, mr, rr)
end

function build_ss_dim(y::Matrix{Typ}, Z::Matrix{Typ}, T::Matrix{Typ}, R::Matrix{Typ}) where Typ
    ny, py = size(y)
    pz, mz = size(Z)
    mt, mt = size(T)
    mr, rr = size(R)
    if !((mz == mt == mr) && (pz == py))
        error("StateSpaceModel dimension mismatch")
    end
    return StateSpaceDimensions(ny, py, mr, rr)
end

function has_unknowns(model::StateSpaceModel{T}) where T
    unknowns = Unknowns(model)
    return unknowns.n_unknowns == 0 ? false : true
end