text:
# Simple graph types

module Graphs

using .Iterators: repeated
using ..Dorothy.Utils

export
		AbstractGraph, Graph, GraphView, neighbors!, pairs, connected,
		connected!, follow, follow!, isisolated, pair!, unpair!, isolate!

abstract type AbstractGraph <: AbstractVector{Vector{Int}} end

struct Graph <: AbstractGraph
	neighbors::Vector{Vector{Int}}

	function Graph(n::Integer)
		n >= 0 || error("expected positive number of indices")
		new(fill(emptyindices, n))
	end
end

struct GraphView{T<:AbstractVector{<:Integer}} <: AbstractGraph
	G::Graph
	I::T

	@inline function GraphView{T}(G::Graph, I::T) where
			{T<:AbstractVector{<:Integer}}
		@boundscheck checkindexseries(G, I)
		new(G, I)
	end
end

function Graph(G::AbstractGraph)
	G2 = Graph(length(G))
	for (i, j) in pairs(G)
		pair!(G2, (i, j))
	end
	G2
end

@inline GraphView(G::Graph, I::T) where {T<:AbstractVector{<:Integer}} =
		GraphView{T}(G, I)

@inline function GraphView(G::GraphView, I::AbstractVector{<:Integer})
	@boundscheck checkindexseries(G, I)
	@inbounds GraphView(G.G, G.I[I])
end

Base.size(G::Graph) = (length(G.neighbors),)

Base.size(G::GraphView) = (length(G.I),)

Base.getindex(G::AbstractGraph, i::Int) = neighbors!(Int[], G, i)

Base.getindex(G::Graph, i::Int) = copy(G.neighbors[i])

Base.IndexStyle(::Type{<:AbstractGraph}) = IndexLinear()

Base.parent(G::Graph) = G

Base.parent(G::GraphView) = G.G

Base.parentindices(G::Graph) = 1:length(G)

Base.parentindices(G::GraphView) = G.I

Base.view(G::AbstractGraph, I::AbstractVector{<:Integer}) = GraphView(G, I)

Base.view(G::AbstractGraph, i::Integer) = GraphView(G, [i])

Base.view(G::Graph, I::AbstractVector{<:Integer}) = GraphView(G, copy(I))

@inline function Base.in((i, j)::Tuple{Integer,Integer}, G::Graph)
	@boundscheck begin
		checkbounds(G, i)
		checkbounds(G, j)
	end
	j in G.neighbors[i]
end

@inline function Base.in((i, j)::Tuple{Integer,Integer}, G::GraphView)
	@boundscheck begin
		checkbounds(G, i)
		checkbounds(G, j)
	end
	(G.I[i], G.I[j]) in G.G
end

@inline function neighbors!(dest::AbstractArray{<:Integer}, G::Graph,
		i::Integer)
	@boundscheck checkbounds(G, i)
	empty!(dest)
	append!(dest, G.neighbors[i])
end

@inline function neighbors!(dest::AbstractArray{<:Integer}, G::GraphView,
		i::Integer)
	@boundscheck checkbounds(G, i)
	neighbors!(dest, G.G, G.I[i])
	for j in reverse(eachindex(dest))
		r = searchsorted(G.I, dest[j])
		if isempty(r)
			deleteat!(dest, j)
		else
			dest[j] = r[1]
		end
	end
	dest
end

function Base.pairs(G::Graph)
	pairs = sizehint!(Tuple{Int,Int}[], length(G))
	for i in eachindex(G.neighbors)
		for j in G.neighbors[i]
			if j > i
				push!(pairs, (i, j))
			end
		end
	end
	pairs
end

function Base.pairs(G::GraphView)
	pairs = sizehint!(Tuple{Int,Int}[], length(G))
	tmp = Int[]
	for i in 1:length(G)
		for j in @inbounds neighbors!(tmp, G, i)
			if j > i
				push!(pairs, (i, j))
			end
		end
	end
	pairs
end

connected(G::AbstractGraph, i::Integer) = connected!(Int[], G, i)

function connected!(dest::AbstractVector{<:Integer}, G::AbstractGraph,
		i::Integer, reachable::AbstractVector{Bool} = trues(length(G)))
	@boundscheck checkbounds(G, i)
	empty!(dest)
	tmp = Int[]
	stack = [i]
	while ! isempty(stack)
		j = pop!(stack)
		if reachable[j]
			push!(dest, j)
			append!(stack, neighbors!(tmp, G, j))
			reachable[j] = false
		end
	end
	dest
end

follow(G::AbstractGraph, (i, j)::Tuple{Integer,Integer}) = follow!([], (i, j))

function follow!(dest::AbstractVector{<:Integer}, G::AbstractGraph,
		(i, j)::Tuple{Integer,Integer})
	reachable = trues(length(G))
	reachable[i] = false
	connected!(dest, G, j, reachable)
end

function isisolated(G::AbstractGraph, I::AbstractArray{<:Integer})
	@boundscheck checkbounds(G, I)
	J = Int[]
	for i in I
		for j in neighbors!(J, G, i)
			if ! (j in I)
				return false
			end
		end
	end
	true
end

@inline function isisolated(G::Graph, i::Integer)
	@boundscheck checkbounds(G, i)
	isempty(G.neighbors[i])
end

@inline function isisolated(G::GraphView, i::Integer)
	@boundscheck checkbounds(G, i)
	@inbounds isisolated(G.G, G.I[i])
end

function pair!(G::Graph, pairs::Tuple{Integer,Integer}...)
	@boundscheck begin
		for (i, j) in pairs
			checkbounds(G, i)
			checkbounds(G, j)
			i != j || error("expected non-identical indices")
		end
	end
	for (i, j) in pairs
		J = G.neighbors[i]
		r = searchsorted(J, j)
		if isempty(r)
			if isempty(J)
				G.neighbors[i] = [j]
			else
				insert!(J, first(r), j)
			end
			I = G.neighbors[j]
			if isempty(I)
				G.neighbors[j] = [i]
			else
				insert!(I, searchsortedfirst(I, i), i)
			end
		end
	end
	G
end

function pair!(G::GraphView, pairs::Tuple{Integer,Integer}...)
	@boundscheck begin
		for (i, j) in pairs
			checkbounds(G, i)
			checkbounds(G, j)
			i != j || error("expected non-identical indices")
		end
	end
	@inbounds pair!(G.G, [(G.I[i], G.I[j]) for (i, j) in pairs]...)
end

function unpair!(G::Graph, pairs::Tuple{Integer,Integer}...)
	@boundscheck begin
		for (i, j) in pairs
			checkbounds(G, i)
			checkbounds(G, j)
			i != j || error("expected non-identical indices")
		end
	end
	for (i, j) in pairs
		J = G.neighbors[i]
		r = searchsorted(J, j)
		if ! isempty(r)
			if length(J) == 1
				G.neighbors[i] = emptyindices
			else
				deleteat!(J, first(r))
			end
			I = G.neighbors[j]
			if length(I) == 1
				G.neighbors[j] = emptyindices
			else
				deleteat!(I, searchsortedfirst(I, i))
			end
		end
	end
	G
end

function unpair!(G::GraphView, pairs::Tuple{Integer,Integer}...)
	@boundscheck begin
		for (i, j) in pairs
			checkbounds(G, i)
			checkbounds(G, j)
			i != j || error("expected non-identical indices")
		end
	end
	@inbounds unpair!(G.G, [(G.I[i], G.I[j]) for (i, j) in pairs]...)
end

function isolate!(G::Graph, i::Integer)
	@boundscheck checkbounds(G, i)
	for j in G.neighbors[i]
		I = G.neighbors[j]
		if length(I) == 1
			G.neighbors[j] = emptyindices
		else
			deleteat!(I, searchsortedfirst(I, i))
		end
	end
	G.neighbors[i] = emptyindices
	G
end

@inline function isolate!(G::GraphView, i::Integer)
	@boundscheck checkbounds(G, i)
	for j in @inbounds G[i]
		unpair!(G, (i, j))
	end
	G
end

function isolate!(G::AbstractGraph, I::AbstractVector{<:Integer})
	@boundscheck checkbounds(G, I)
	J = Int[]
	for i in I
		for j in neighbors!(J, G, i)
			if ! (j in I)
				unpair!(G, (i, j))
			end
		end
	end
	G
end

function Base.:(==)(G1::AbstractGraph, G2::AbstractGraph)
	if length(G1) == length(G2)
		J1 = Int[]
		J2 = Int[]
		for i in eachindex(G1)
			J1 = neighbors!(G1, i)
			J2 = neighbors!(G2, i)
			if J1 != J2
				return false
			end
		end
		true
	else
		false
	end
end

function Base.merge!(G::AbstractGraph, src::AbstractGraph)
	@boundscheck length(src) <= length(G) ||
			error("cannot merge source graph larger than destination")
	pair!(G, pairs(src)...)
	G
end

function Base.resize!(G::Graph, n::Integer)
	@boundscheck n >= 0 || error("expected positive length")
	nprev = length(G)
	dn = n - nprev
	if dn > 0
		resize!(G.neighbors, n)
		for i = nprev+1:n
			G.neighbors[i] = emptyindices
		end
	elseif dn < 0
		for i = n+1:nprev
			for j in G.neighbors[i]
				I = G.neighbors[j]
				if length(I) == 1
					G.neighbors[j] = emptyindices
				else
					deleteat!(I, searchsortedfirst(I, i))
				end
			end
			G.neighbors[i] = emptyindices
		end
		resize!(G.neighbors, n)
	end
	G
end

function Base.empty!(G::Graph)
	empty!(G.neighbors)
	G
end

function remap!(G::Graph, map::AbstractVector{<:Integer})
	newneighbors = Int[]
	for (i0, i1) in enumerate(map)
		if i1 != 0
			empty!(newneighbors)
			for j0 in G.neighbors[i0]
				j1 = map[j0]
				if j1 != 0
					push!(newneighbors, j1)
				end
			end
			if isempty(newneighbors)
				G.neighbors[i0] = emptyindices
			else
				G.neighbors[i0] = copy(newneighbors)
			end
		end
	end
end

Base.deleteat!(G::Graph, i::Integer) = deleteat!(G, [i])

@inline function Base.deleteat!(G::Graph, I::AbstractArray{<:Integer})
	@boundscheck checkindexseries(G, I)
	remap!(G, deleteatmap(length(G), I))
	deleteat!(G.neighbors, I)
	G
end

Base.splice!(G::Graph, i::Integer) = deleteat!(G, i)

Base.splice!(G::Graph, i::Integer, replacement::AbstractGraph) =
		splice!(G, i:i, replacement)

Base.splice!(G::Graph, range::AbstractUnitRange{<:Integer}) =
		deleteat!(G, range)

function Base.splice!(G::Graph, range::AbstractUnitRange{<:Integer},
		replacement::AbstractGraph)
	@boundscheck begin
		checkbounds(G, range)
		parent(replacement) !== G ||
				error("replacement must belong to a different graph")
	end
	spliced = Graph(view(G, range))
	nreplacement = length(replacement)
	if nreplacement == 0
		@inbounds deleteat!(G, range)
	else
		remap!(G, splicemap(length(G), range, nreplacement))
		splice!(G.neighbors, range, repeated(emptyindices, nreplacement))
		offset = first(range) - 1
		for (i, j) in pairs(replacement)
			pair!(G, (i+offset, j+offset))
		end
	end
	spliced
end

function Base.append!(G::Graph, src::AbstractGraph)
	n = length(G)
	splice!(G, n+1:n, src)
	G
end

function Base.prepend!(G::Graph, src::AbstractGraph)
	splice!(G, 1:0, src)
	G
end

function Base.insert!(G::Graph, i::Integer, src::AbstractGraph)
	splice!(G, i:i-1, src)
	G
end

end # module


text:
################################################################################
#
#  Submodules
#
################################################################################

include("../EllCrv/FieldsRings.jl")
include("../EllCrv/EllCrv.jl")
include("../EllCrv/EllCrvDiv.jl")
include("../EllCrv/EllCrvMap.jl")
include("../EllCrv/EllCrvModel.jl")

################################################################################
#
#  Imports
#
################################################################################

import AbstractAlgebra

# An aside from the proof of Proposition 4.4, verifying
# the non-vanishing of a constructed determinant given
# a divisor D1 on E

function verify(prime::UInt128)

    if prime == 2

        # Fix a field of prime order.
        Z = AbstractAlgebra.GF(prime)

        # Fix homogeneous coordinates [x:y:z] for an elliptic curve,
        # where x, y, z are in Z.
        coords = ("x1", "y1", "z1")

        # Fix coefficients A, B in Z for the Weierstrass form of the
        # elliptic curve, i.e. in homogenous coordinates [x:y:z],
        # EllCrv := x^3 + Ax^2 + B - x*y.

        coeff = (A, B)

        # Instantiate Polynomial Rings.
        R0, (A, B, xQ, yQ, coords) = AbstractAlgebra.PolynomialRing(Z,
            ["A",
            "B",
            "xQ",
            "yQ",
            "x1",
            "y1",
            "z1"])

            # R1, aq1, aq, a1, a0, r = AbstractAlgebra.PolynomialRing(R0,
        #    "aq1",
        #    "aq",
        #    "a1",
        #    "a0",
        #    "r")

        # Instantiate Short-form Weierstrass Elliptic Curve.
        check = true
        EllCrvWeier = EllipticCurve(coeff, check)

    end

end

ext: jl


text:
# using GLMakie
# using CairoMakie
using Plots

function field_plot(x; title = "", kwargs...)
    l = size(x,4)
    plot(
        [
            plot(
        x,
                title = "$title\ncomponent $i",
                c = cgrad([:blue, :red], [0.0, 0.5, 1.0]),
                kwargs...,
            ) for x in eachslice(x,dims=4)
        ]...,
        layout = l,
    )
    # fig = Figure()
    # fig = Figure(
    #     backgroundcolor = RGBf0(0.98, 0.98, 0.98),
    #     resolution = (1000, 700),
    # )
    # for i = 1:length(X)
    #     volume(fig[1, i], X[i])
    # end
    # plots = [volume(X[i]) for i = 1:length(X)]
    # plot(plots..., layout = length(X))
end

function vis(x::EquivConv)
    r=LinRange(0.,x.rmax,64)
    plots = [Plots.plot(r,p.op.radfunc.(r)) for p in x.paths]
    Plots.plot(plots..., layout = length(plots))
end


text:
using Pages
using HTTP
using JSON
using Dates

Endpoint("/hello") do request::HTTP.Request
"Hello world"
end

Endpoint("/echo") do request::HTTP.Request
        params = HTTP.queryparams(HTTP.URI(request.target).query)
        println("Body: $(params)")
        open("logfile.txt", "a") do io
            println(io, params["id"])
        end
        response = JSON.json(params)
end

Endpoint("/echo",POST) do request::HTTP.Request
        data = String(request.body)
        open("logfile.txt","a") do io
		println("Parameters: $(data)")
		print(io, data)
        println(io, Dates.format(now(), " dd/mm/yyyy HH:MM:ss"))
        end
        response = JSON.json(Dict(:data => data, "time" => Dates.format(now(), " dd/mm/yyyy HH:MM:ss")))
end

Endpoint("/") do request::HTTP.Request
    read(joinpath(@__DIR__,"index.html"),String)
end

open("logfile.txt", "r") do io
    a = readlines(io)
    # print([parse(Int,x[1]) for x in a])
    print(a)
end

#@async Pages.start(8000);
Pages.start(8000);


text:
using Test, LinearMaps, LinearAlgebra, BlockArrays

@testset "Non-traditional axes" begin

    A = rand(ComplexF64,2,4)
    B = PseudoBlockMatrix{ComplexF64}(undef, [2,3], [3,4])

    ax1 = axes(B)[1]
    ax2 = axes(B)[2]
    fill!(B,0)
    B[Block(1),Block(2)] .= A

    N = @inferred LinearMap(B)
    @test axes(N) == (ax1,ax2)
    @test axes(N, 1) == ax1
    @test axes(N, 2) == ax2
    @test_throws ErrorException axes(N, 3)

    @test eltype(N) == eltype(B)

    u = similar(Array{ComplexF64}, ax2)
    v = PseudoBlockVector{ComplexF64}(undef, [3,5])
    w = similar(Array{ComplexF64}, ax1)

    for i in eachindex(u) u[i] = rand(ComplexF64) end
    for i in eachindex(v) v[i] = rand(ComplexF64) end
    for i in eachindex(w) w[i] = rand(ComplexF64) end

    @test B*u == N*u
    @test_throws DimensionMismatch N*v

    # Lu = L*u
    Nu = N*u

    @test axes(Nu)[1] == axes(N)[1] == axes(B)[1]
    @test blocklengths(axes(Nu)[1]) == blocklengths(axes(N)[1]) == blocklengths(axes(B)[1]) == [2,3]

    for trans in (adjoint, transpose)
        Nt = trans(LinearMap(N))
        Bt = trans(B)
        Ntw = Nt*w
        Btw = Bt*w
        @test Ntw ≈ Btw
        @test axes(Ntw)[1] == axes(Nt)[1] == axes(Bt)[1]
        @test blocklengths(axes(Ntw)[1]) == blocklengths(axes(Nt)[1]) == blocklengths(axes(Bt)[1]) == [3,4]
    end

    C = B + 2N
    @test axes(C) === axes(B) === axes(N)
    @test C*u ≈ 3*Nu

    Cu = C*u
    @test axes(C)[1] == ax1
    @test blocklengths(axes(C)[1]) == blocklengths(ax1)

    A = rand(ComplexF64,2,2)
    B = PseudoBlockMatrix{ComplexF64}(undef, [2,2], [2,2])
    ax1 = axes(B)[1]
    ax2 = axes(B)[2]
    fill!(B,0)
    B[Block(1),Block(2)] .= A
    L = LinearMap(B)
    L2 = L*L
    B2 = B*B
    @test axes(L2) == axes(B2)
    B2 = B*Matrix(B)
    L2 = L*LinearMap(Matrix(B))
    @test axes(L2) == axes(B2)
    u = similar(Array{ComplexF64}, ax2)
    B2u = B2*u; L2u = L2*u
    @test axes(B2u)[1] == axes(L2u)[1] == axes(B2)[1] == axes(L2)[1]
    @test blocklengths(axes(B2u)[1]) == blocklengths(axes(L2u)[1]) == [2,2]
end


text:
function EuropeanCallMC(K, T, S, sig, r, div, N, M)

## setting parameters

dt = T/N
nudt = (r - div - 0.5 * sig^2) * dt
sigsdt = sig * sqrt(dt)
lnS = log(S)

## variables to store the sums of the option values

sum_CT = 0
sum_CT2 = 0

## outer loop is the number of simulations bigger is better

for j = 1:M

  lnSt = lnS

    ## inner loop is simulating each price path of length N
  for i = 1:N

    lnSt = lnSt + nudt + sigsdt * randn() ## note, it only add the changes, does not store the path

  end

    ## wraps up and stores all the values of each path
  ST = exp(lnSt)
  CT = max(0, ST - K)
  sum_CT = sum_CT + CT
  sum_CT2 = sum_CT2 + CT^2

end

  ## discounting the average call value
call_value = sum_CT / M * exp(-r*T)

SD = (sqrt(sum_CT2 - 1/M * (sum_CT^2)) * exp(-2*r*T)) / (M-1)
SE = SD / sqrt(M)

println("The Call Value is $call_value")
println("The SD is $SD")
println("The SE is $SE")

MC = [call_value, SD, SE]

end


text:
@doc raw"""
    UnweightedParticleEnsemble{T}

An ensemble of `N` particles, each of dimension `n`.
"""
struct UnweightedParticleEnsemble{T} <: UnweightedParticleRepresentation{Vector{T}}
    positions::Matrix{T}
end

UnweightedParticleEnsemble(vec::Vector{T}, N::Int) where T<:Number = UnweightedParticleEnsemble(repeat(vec, 1, N))

particle_dim(ens::UnweightedParticleEnsemble)    = size(ens.positions, 1)

no_of_particles(ens::UnweightedParticleEnsemble) = size(ens.positions, 2)

get_pos(ens::UnweightedParticleEnsemble, i) = view(ens.positions, :, i)

Base.show(io::IO, ens::UnweightedParticleEnsemble) = print(io, "Unweighted particle ensemble
    # of particles: ", no_of_particles(ens),"
    particle type:  ", particle_dim(ens),"-dimensional ", eltype(ens))

function UnweightedParticleEnsemble(model::HiddenStateModel, N::Int)
    return UnweightedParticleEnsemble(hcat([initialize(model) for i in 1:N]...))
end

function propagate!(ens::UnweightedParticleEnsemble{S}, model::HiddenStateModel{Vector{S}, ContinuousTime}, dt) where S
    propagate!(ens.positions, model, dt)
end

mean(ens::UnweightedParticleEnsemble) = Statistics.mean(ens.positions, dims=2)
cov(ens::UnweightedParticleEnsemble)  = Statistics.cov(ens.positions, dims=2, corrected=false)
var(ens::UnweightedParticleEnsemble)  = Statistics.var(ens.positions, dims=2, corrected=false)

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to create a printing recipe for the lines.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Convert `Decotration` to string.
function convert(::Type{String}, d::Decoration)

    # Check if we have a reset.
    _is_default(d) && return ""
    d.reset && return "$(CSI)0m"

    str = ""
    str *= isempty(d.foreground) ? "" : "$(CSI)" * d.foreground * "m"
    str *= isempty(d.background) ? "" : "$(CSI)" * d.background * "m"
    str *= "$(CSI)" * (d.bold ? "1" : "22") * "m"
    str *= "$(CSI)" * (d.underline ? "4" : "24") * "m"
    str *= "$(CSI)" * (d.reversed ? "7" : "27") * "m"

    return str
end

string(d::Decoration) = convert(String, d)

"""
    _printing_recipe(str::AbstractString, start_char::Int, max_chars::Int, highlight_matches::SearchMatches, active_match::SearchMatch)

Create a printing recipe of the line `str`. The recipe is composed of:

* A vector with the strings;
* A vector with objects of type `Decoration` describing the decoration of each
  string.

This function also return the number of cropped characters in this line.

"""
function _printing_recipe(
    str::String,
    start_char::Int,
    max_chars::Int,
    highlight_matches::SearchMatches,
    active_match::SearchMatch,
    freeze_columns::Int
)
    # Current state.
    decoration::Decoration = Decoration()

    # Vector containing the decorations.
    d = Vector{Decoration}(undef, 0)

    # Vetor containing the strings.
    s = Vector{String}(undef, 0)

    # Initial state for the state machine.
    state = :string

    # Indicate that we reached end of line.
    eol = false

    # Variables to handle highlighting.
    num_highlights = length(highlight_matches)
    hl_state = :normal
    old_decoration::Decoration = Decoration()

    # Vectors with the beginning and end of the highlights.
    hl_beg = map(x->x[2], highlight_matches)
    hl_end = hl_beg .+ map(x->x[3] - 1, highlight_matches)
    hl_active_id = 0

    for i = 1:num_highlights
        highlight_matches[i] === active_match && (hl_active_id = i)
    end

    # Current string.
    str_i = ""

    # Inform that we have a new decoration.
    new_decoration = false

    # Number of processed characters in the string.
    num_processed_chars = 0

    # Number of printed characters in the display.
    num_printed_chars = 0

    # State of string processing.
    string_state = :initial_cropping

    # Current escape sequence code.
    code = ""

    # Make sure the parameters have correct values.
    start_char < 1 && (start_char = 1)

    # Variable to compute the number of cropped chars.
    cropped_chars = 0

    # Process the freeze columns.
    if freeze_columns > 0
        freeze_columns > max_chars && (freeze_columns = max_chars)
        start_char < freeze_columns && (start_char = freeze_columns)
        string_state = :freeze
    end

    for c in str
        # Check if we have a escape charater. In this case, change the state to
        # `:escape_seq`. We also need to add the string assembled so far to the
        # output vector.
        if c == '\e'
            state = :escape_seq_beg

            # If we are not highlighting somehting, we have at least one space
            # in the screen, and the current string is not empty, then we need
            # to flush the string and decoration to start a new segment.
            if (hl_state == :normal) &&
               ((string_state == :view) || (string_state == :freeze))&&
               !isempty(str_i)
                push!(s, str_i)
                push!(d, decoration)
                str_i = ""
                new_decoration = false
            end

            continue
        end

        # If we are not in the state `:escape_seq`, then just add the character
        # to the string.
        if state == :string
            cw = textwidth(c)

            # If this character matches the beginning of highlight, we need to
            # flush the string until now and start a new segment.
            if (hl_state == :normal)

                hl_id = 0

                for i = 1:num_highlights
                    hl_beg[i] == (num_processed_chars + cw) && (hl_id = i)
                end

                if hl_id > 0
                    hl_state = :highlight

                    push!(s, str_i)
                    push!(d, decoration)
                    str_i = ""

                    old_decoration = decoration
                    decoration = _search_highlighting[hl_id == hl_active_id]
                end
            end

            num_processed_chars += cw

            # Check if we are in the freeze columns.
            if string_state == :freeze
                if num_printed_chars + cw > freeze_columns
                    string_state = :initial_cropping

                    push!(s, str_i)
                    push!(d, decoration)
                    str_i = ""
                    new_decoration = false
                end
            end

            # Check if we finished the initial cropping with this character.
            if string_state == :initial_cropping
                num_processed_chars ≥ start_char && (string_state = :view)
            end

            # Check if we reached the final of the viewable area.
            if string_state == :view
                num_printed_chars + cw > max_chars && (string_state = :final_cropping)
            end

            if (string_state == :view) || (string_state == :freeze)
                str_i *= string(c)
                num_printed_chars += cw
            end

            if hl_state == :highlight
                hl_id = 0

                for i = 1:num_highlights
                    hl_end[i] == num_processed_chars && (hl_id = i)
                end

                # TODO: If there is a highlight that starts inside the old one,
                # we could save a push here.
                if hl_id > 0
                    push!(s, str_i)
                    push!(d, decoration)
                    str_i = ""
                    decoration = old_decoration

                    # Add an empty segment to reset all formatting created by
                    # the highlight.
                    push!(s, "")
                    push!(d, Decoration(reset = true))

                    hl_state = :normal
                end
            end

            # Check if we are in the final character to be printed.
            if string_state == :final_cropping
                if !eol
                    push!(s, str_i)
                    push!(d, decoration)
                    str_i = ""
                    eol = true
                end

                cropped_chars += cw
            end

        elseif state == :escape_seq_beg
            state = c == '[' ? :escape_seq : :string

        elseif state == :escape_seq
            if isdigit(c) || c == ';'
                code *= string(c)

            elseif c == 'm'
                state = :string

                if hl_state == :highlight
                    old_decoration = _parse_ansi_code(old_decoration, code)
                else
                    decoration = _parse_ansi_code(decoration, code)
                end

                code = ""
                new_decoration = true

            else
                state = :string
                code = ""
            end
        end
    end

    if !isempty(str_i)
        push!(s, str_i)
        push!(d, decoration)
    elseif new_decoration
        push!(d, decoration)
    end

    return s, d, cropped_chars
end

################################################################################
#                              Private functions
################################################################################

function _parse_ansi_code(decoration::Decoration, code::String)
    tokens = split(code, ';')

    # Unpack fields.
    foreground = decoration.foreground
    background = decoration.background
    bold       = decoration.bold
    underline  = decoration.underline
    reset      = decoration.reset
    reversed   = decoration.reversed
    force      = decoration.force

    # `reset` must not be copied to other decorations. Hence, we need to reset
    # it here.
    reset = false
    force = false

    i = 1
    while i ≤ length(tokens)
        code_i = tryparse(Int, tokens[i], base = 10)

        if code_i == 0
            # If we have a reset, neglect all the other configurations.
            return Decoration(reset = true)
        elseif code_i == 1
            bold = true
        elseif code_i == 4
            underline  = true
        elseif code_i == 7
            reversed = true
        elseif code_i == 22
            bold = false
            force = true
        elseif code_i == 24
            underline = false
            force = true
        elseif code_i == 27
            reversed = false
            force = true
        elseif 30 <= code_i <= 37
            foreground = "$code_i"
        # 256-color support for foreground.
        elseif code_i == 38
            # In this case, we can have an extended color code. To check this,
            # we must have at least two more codes.
            if i+2 ≤ length(tokens)
                code_i_1 = tryparse(Int, tokens[i+1], base = 10)
                code_i_2 = tryparse(Int, tokens[i+2], base = 10)

                if code_i_1 == 5
                    foreground = "38;5;$code_i_2"
                end

                i += 2
            end
        elseif code_i == 39
            foreground = "39"
        elseif 40 <= code_i <= 47
            background = "$code_i"
        # 256-color support for background.
        elseif code_i == 48
            # In this case, we can have an extended color code. To check this,
            # we must have at least two more codes.
            if i+2 ≤ length(tokens)
                code_i_1 = tryparse(Int, tokens[i+1], base = 10)
                code_i_2 = tryparse(Int, tokens[i+2], base = 10)

                if code_i_1 == 5
                    background = "48;5;$code_i_2"
                end

                i += 2
            end
        elseif code_i == 49
            background = "49"
        # Bright foreground colors defined by Aixterm.
        elseif 90 <= code_i <= 97
            foreground = "$code_i"
        # Bright background colors defined by Aixterm.
        elseif 100 <= code_i <= 107
            background = "$code_i"
        end

        i += 1
    end

    return Decoration(
        foreground,
        background,
        bold,
        underline,
        reset,
        reversed,
        force
    )
end

_is_default(d::Decoration) = d === Decoration()


text:
using Test, Random, AdvancedHMC, ForwardDiff, AbstractMCMC
using Statistics: mean
include("common.jl")

@testset "`gdemo`" begin
    rng = MersenneTwister(0)

    n_samples = 5_000
    n_adapts = 5_000

    θ_init = randn(rng, 2)

    model = AdvancedHMC.DifferentiableDensityModel(ℓπ_gdemo, ForwardDiff)
    init_eps = Leapfrog(1e-3)
    κ = NUTS(init_eps)
    metric = DiagEuclideanMetric(2)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, κ.τ.integrator))

    samples = AbstractMCMC.sample(
        rng, model, κ, metric, adaptor, n_adapts + n_samples;
        nadapts = n_adapts,
        init_params = θ_init,
        progress=false,
        verbose=false
    );

    # Transform back to original space.
    # NOTE: We're not correcting for the `logabsdetjac` here since, but
    # we're only interested in the mean it doesn't matter.
    for t in samples
        t.z.θ .= invlink_gdemo(t.z.θ)
    end
    m_est = mean(samples[n_adapts + 1:end]) do t
        t.z.θ
    end

    @test m_est ≈ [49 / 24, 7 / 6] atol=RNDATOL
end


text:
using MakieGallery, AbstractPlotting, GLMakie, Test
using MakieGallery: @block, @cell
empty!(MakieGallery.plotting_backends)
push!(MakieGallery.plotting_backends, "Makie")
database = MakieGallery.load_database()

exclude = (
    "Cobweb plot", # has some weird scaling issue on CI
    "Colormap collection", # has one size different...
    # doesn't match 0.035520551315007046 <= 0.032. Looked at the artifacts and it looks fairly similar
    # so blaming video compression
    "Interaction with Mouse"
)
# Download is broken on CI
filter!(entry-> !("download" in entry.tags) && !(entry.title in exclude), database)

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end
recordings = MakieGallery.record_examples(test_record_path)
@test length(recordings) == length(database)
MakieGallery.run_comparison(test_record_path, tested_diff_path)


text:
# read the output power
export psp_readpower

const READ_POWER = [0x57, 0x0d]

"""
    psp_readpower(io)

Read the output power.
"""
function psp_readpower(io_psp::IO)
  write(io_psp, READ_POWER)
  powerstring = bytestring(read(io_psp,UInt8, 9))
  if  powerstring[1] != 'W' ||
      powerstring[5] != '.' ||
      powerstring[7:9] != "\r\r\n"
        error("PSP Invalid Response: ",powerstring)
  end
  parse(Float64,powerstring[2:6])
end


text:
# resample(uvals::Vector{AbstractUncertainValue})
@test resample(uvals) isa Vector
@test length(resample(uvals)) == length(uvals)

# resample(uvals::Vector{AbstractUncertainValue}, n::Int)
n = 2
@test resample(uvals, n) isa Vector{Vector{T}} where T<:Real
# TODO: make sure all values are promoted to floats, otherwise performance loss is
# significant
#@test resample(uvals, n) isa Vector{Vector{Float}}

@test length(resample(uvals, n)) == n

# resample_elwise(uvals::Vector{AbstractUncertainValue}, n::Int)
n = 2
@test resample_elwise(uvals, n) isa Vector{Vector{T}} where T<:Real
@test length(resample_elwise(uvals, n)) == length(uvals)
@test length(resample_elwise(uvals, n)[1]) == n


text:
using SpectralGM
using Test

@testset "SpectralGM.jl" begin
    # Write your tests here.
end


text:
#region Uncertainties filtered parameters
# using FiniteDiff
# Base.eps(Real) = eps(Float64)

function A0_B0_est_for_white_cov_mat_obj_SD_filter_time_seq(model::SdErgm, N, obsT, vEstSdResParAll, ftot_0)

    T = length(obsT)
    nPar = length(vEstSdResParAll)
    gradT = zeros(nPar, T)
    hessT = zeros(nPar, nPar, T)

    vecUnParAll = unrestrict_all_par(model, vEstSdResParAll)

    function obj_fun_t(xUn, t)
        logLikeVecT = seq_loglike_sd_filter(model, N, obsT[1:t], xUn, ftot_0)
        return - logLikeVecT[end]
    end

    for t = 1:T

        obj_fun(vecUnParAll) = obj_fun_t(vecUnParAll, t)

        gradT[:,t] = deepcopy(ForwardDiff.gradient(obj_fun, vecUnParAll))
        hessT[:,:,t] =  deepcopy(ForwardDiff.hessian(obj_fun, vecUnParAll))
        # gradT[:,t] = deepcopy(FiniteDiff.finite_difference_gradient(obj_fun, vecUnParAll))
        # hessT[:,:,t] =  deepcopy(FiniteDiff.finite_difference_hessian(obj_fun, vecUnParAll))
    end

    OPGradSum = sum([gt * gt' for gt in eachcol(gradT[:,2:end])] )
    HessSum = dropdims(sum(hessT[:,:,2:end], dims=3 ), dims=3)

    return OPGradSum, HessSum
end

function white_estimate_cov_mat_static_sd_par(model::SdErgm,  N,obsT, indTvPar, ftot_0, vEstSdResParAll; returnAllMats=false, enforcePosDef = true)

    T = length(obsT)
    nErgmPar = number_ergm_par(model)
    errorFlag = false

    OPGradSum, HessSum = A0_B0_est_for_white_cov_mat_obj_SD_filter_time_seq(model, N, obsT, vEstSdResParAll, ftot_0)

    parCovHat = pinv(HessSum) * OPGradSum * pinv(HessSum)

    if enforcePosDef
        parCovHatPosDef, minEigenVal = make_pos_def(parCovHat)
    else
        parCovHatPosDef = parCovHat
        minEigenVal = 0
    end

    if minEigenVal < 0
        minEiegOPGrad = minimum(eigen(OPGradSum).values)
        minEiegHess = minimum(eigen(HessSum).values)
        Logging.@info("Outer Prod minimum eigenvalue $(minEiegOPGrad) , hessian minimum eigenvalue $(minEiegHess)")

        # if the negative eigenvalue of the cov mat is due to a negative eigenvalue of the hessian, do not use that estimate
        if minEiegHess < 0
            errorFlag = true
        end
    end

    mvSDUnParEstCov = Symmetric(parCovHatPosDef)
    if returnAllMats
        return mvSDUnParEstCov, errorFlag, OPGradSum, HessSum
    else
        return mvSDUnParEstCov, errorFlag
    end
end

function get_B_A_mats_for_TV(model::SdErgm, indTvPar, vEstSdResPar)
    # @debug "[get_B_A_mats_for_TV][begin]"
    nTvPar = sum(indTvPar)
    nErgmPar = number_ergm_par(model)
    B = zeros(nTvPar)
    A = zeros(nTvPar)

    lastIndRead = 0
    lastIndWrite = 0
    for e in 1:nErgmPar
        if indTvPar[e]
            bInd = lastIndRead+2
            aInd = lastIndRead+3
            B[lastIndWrite+1] = vEstSdResPar[bInd]
            A[lastIndWrite+1] = vEstSdResPar[aInd]
            lastIndRead += 3
            lastIndWrite += 1
        else
            lastIndRead += 1
        end

    end

    # @debug "[get_B_A_mats_for_TV][end]"
    return B, A
end

function distrib_filtered_par_from_sample_static_par(model::SdErgm, N, obsT, indTvPar, ftot_0, sampleUnParAll)

    T = length(obsT)
    nErgmPar = number_ergm_par(model)
    nTvPar = sum(indTvPar)
    nSample = size(sampleUnParAll)[2]

    # any(.!indTvPar) ? error("assuming all parameters TV. otherwise need to set filter unc to zero for static ones") : ()

    distribFilteredSD = zeros(nSample, nErgmPar,T)
    filtCovHatSample = zeros(nTvPar, T, nSample)
    flagIntegrated = get_option(model,  "integrated")
    for n=1:nSample
        vResPar = restrict_all_par(model, sampleUnParAll[:,n])

        vecSDParRe, vConstPar = divide_SD_par_from_const(model, vResPar)

        distribFilteredSD[n, :, :] , ~, ~, invScalMatT = DynNets.score_driven_filter( model, N, obsT, vecSDParRe, indTvPar;ftot_0=ftot_0, vConstPar=vConstPar)

        BMatSD, AMatSD = get_B_A_mats_for_TV(model, indTvPar, vResPar)

        constFiltUncCoeff = (BMatSD.^(-1)).*AMatSD

        if model.scoreScalingType[end] == 'D'
            for t in 1:T
                filtCovHatSample[:,t,n] = constFiltUncCoeff./ invScalMatT[:,:,t][I(nErgmPar)][indTvPar]
            end
        else
            error()
        end

    end

    # remove samples that resulted in inconsistent P_t
    indBadSamples = [any(filtCovHatSample[:,:,n].<0) | any(.!isfinite.(filtCovHatSample[:,:,n])) for n in 1:nSample]

    if mean(indBadSamples) > 0.5
        Logging.@error("Sampling from mvNormal resulted in $(mean(indBadSamples)*100) % of bad samples. Rejected")
        errFlag = true
    else
        errFlag = false
    end

    filtCovHatSample = filtCovHatSample[:,:,.!indBadSamples]

    return distribFilteredSD, filtCovHatSample, errFlag
end

function distrib_static_par_from_mv_normal(model::SdErgm, N, obsT, indTvPar, ftot_0, vEstSdUnPar, mvSDUnParEstCov; nSample = 500)

    T = length(obsT)
    nErgmPar = number_ergm_par(model)
    nTvPar = sum(indTvPar)
    # any(.!indTvPar) ? error("assuming all parameters TV. otherwise need to set filtering unc to zero for static ones") : ()

    zeroMeanSample =  rand(MvNormal(mvSDUnParEstCov), nSample )
    sampleUnParAll = zeroMeanSample .+ vEstSdUnPar

    return sampleUnParAll
end

"""
Use parametric bootstrap to sample from the distribution of static parameters of the SD filter
"""
function par_bootstrap_distrib_filtered_par(model::SdErgm, N, obsT, indTvPar, ftot_0, vEstSdResPar; nSample = 250, logFlag = false)

    T = length(obsT)
    nErgmPar = number_ergm_par(model)

    nStaticSdPar = length(vEstSdResPar)
    vUnParEstDistrib = SharedArray(zeros(length(vEstSdResPar), nSample))
    distribFilteredSD = SharedArray(zeros(nSample, nErgmPar,T))
    filtCovHatSample = SharedArray(zeros(nErgmPar, nSample))
    errFlagVec = SharedArray{Bool}((nSample))
    flagIntegrated = get_option(model,  "integrated")

    Threads.@threads for n=1:nSample

        # try
            # sample SD dgp
            fVecTDgp, A_T_dgp, ~ = DynNets.score_driven_dgp( model, N, vEstSdResPar, indTvPar, (N,T))
            obsTNew = [stats_from_mat(model, A_T_dgp[:,:,t]) for t in 1:T ]

            # Estimate SD static params from SD DGP
            arrayAllParHat, conv_flag,UM, ftot_0 = estimate(model, N, obsTNew; indTvPar=indTvPar, indTargPar=falses(2), ftot_0=ftot_0)

            vResPar = DynNets.array_2_vec_all_par(model, arrayAllParHat, indTvPar)

            vUnPar = unrestrict_all_par(model, deepcopy(vResPar))

            vUnParEstDistrib[:,n] =  vUnPar

            vecSDParRe, vConstPar = divide_SD_par_from_const(model, vResPar)

            distribFilteredSD[n, :, :] , ~, ~ = DynNets.score_driven_filter( model, N, obsT, vecSDParRe, indTvPar; ftot_0=ftot_0, vConstPar=vConstPar)

            BMatSD, AMatSD = get_B_A_mats_for_TV(model, indTvPar, vResPar)

            filtCovHat = (BMatSD.^(-1)).*AMatSD
            filtCovHat[.!indTvPar] .= 0
            filtCovHatSample[:,n] = filtCovHat

        # catch
        #     errFlagVec[n] = true
        # end

    end

   return distribFilteredSD, filtCovHatSample, errFlagVec, vUnParEstDistrib
end

function non_par_bootstrap_distrib_filtered_par(model::SdErgm, N, obsT, indTvPar, ftot_0; nBootStrap = 100)

    vEstSdUnParBootDist = SharedArray(zeros(3*sum(model.indTvPar), nBootStrap))
    T = length(obsT)
    @time Threads.@threads for k=1:nBootStrap
        vEstSdUnParBootDist[:, k] = rand(1:T, T) |> (inds->(inds |> x-> DynNets.estimate(model, N, obsT; indTvPar=indTvPar, ftot_0 = ftot_0, shuffleObsInds = x) |> x-> getindex(x, 1) |> x -> DynNets.array_2_vec_all_par(model, x, model.indTvPar))) |> x -> DynNets.unrestrict_all_par(model, x)
    end

    return vEstSdUnParBootDist
end

function conf_bands_buccheri(model::SdErgm, N, obsT, indTvPar, fVecT_filt, distribFilteredSD, filtCovHatSample, quantilesVals::Vector{Vector{Float64}}, dropOutliers::Bool; nSample = 250, winsorProp=0.005)
    @debug "[conf_bands_buccheri][begin]"

    T = length(obsT)
    nErgmPar = number_ergm_par(model)

    parUncVarianceT = zeros(nErgmPar,T)
    filtUncVarianceT = zeros(nErgmPar,T)

    distribFilteredSD[isnan.(distribFilteredSD)] .= 0
    fVecT_filt[isnan.(fVecT_filt)] .= 0

    any(isnan.(filtCovHatSample)) ? error() : ()

    filtCovDiagHatMean = dropdims(mean(filtCovHatSample, dims=3), dims=3)
    @debug "[conf_bands_buccheri][size(filtCovDiagHatMean)=$(size(filtCovDiagHatMean))]"

    #for each time compute the variance of the filtered par under the normal distrib of static par
    indAmongTV = 0
    for k=1:nErgmPar
        indTvPar[k] ? indAmongTV += 1 : ()

        for t=1:T
            a_t_vec = distribFilteredSD[:,k,t]
            aHat_t = fVecT_filt[k,t]

            # add filtering and parameter unc
            diff = a_t_vec[isfinite.(a_t_vec)] .- aHat_t
            extr = quantile(diff, [0.01, 0.99])
            if dropOutliers & (extr[1] != extr[2])
                diffNoOutl = filter(x->extr[1] < x < extr[2], diff)
            else
                diffNoOutl = diff
            end

            # parUncVarianceT[k, t] = MinCovDet().fit(diffNoOutl).covariance_

            parUncVarianceT[k, t] = var(winsor(diffNoOutl, prop=winsorProp))

            isnan(parUncVarianceT[k, t]) ? (@show a_t_vec; @show aHat_t; @show diff; @show extr; @show t; @show k; error()) : ()

            if indTvPar[k]
                filtUncVarianceT[k, t] = filtCovDiagHatMean[indAmongTV,t]
            end
        end
    end

    nBands = length(quantilesVals)
    nQuant = 2*nBands

    confBandsPar = repeat(fVecT_filt, outer=(1,1,nBands, 2))
    confBandsParFilt = repeat(fVecT_filt, outer=(1,1,nBands, 2))

    for p = 1:number_ergm_par(model)
        for t=1:T
            sigmaPar = sqrt(parUncVarianceT[p,t])
            sigmaFiltPar = sqrt(parUncVarianceT[p,t] + filtUncVarianceT[p,t])
            for b=1:nBands
                if sigmaPar !== 0
                    confBandsPar[p, t, b,  :] = quantile.(Normal(fVecT_filt[p, t], sigmaPar), quantilesVals[b])
                end
                # confBandsPar[p, t, b,  :] = fVecT_filt[p, t] .+  [1, -1] .* (1.96 * sigma)
                confBandsParFilt[p, t, b,  :] = quantile.(Normal(fVecT_filt[p, t], sigmaFiltPar), quantilesVals[b])
                # confBandsParFilt[p, t, b,  :] = fVecT_filt[p, t] .+  [1, -1] .* (1.96 * sigma)

                # catch e
                #     Logging.@error( (parUncVarianceT[p,t] , filtUncVarianceT[p,t], filtCovDiagHatMean))

                #     error(e)
                # end
            end
        end
    end
    return confBandsParFilt, confBandsPar
end

function conf_bands_par_uncertainty_blasques(model::SdErgm, obsT, fVecT_filt, distribFilteredSD, quantilesVals::Vector{Vector{Float64}}; nSample = 500)

    T = length(obsT)
    nErgmPar = number_ergm_par(model)

    parUncVarianceT = zeros(nErgmPar,T)

    nBands = length(quantilesVals)
    nQuant = 2*nBands
    confBandsPar = repeat(fVecT_filt, outer=(1,1,nBands, 2))

    distribFilteredSD[isnan.(distribFilteredSD)] .= 0
    fVecT_filt[isnan.(fVecT_filt)] .= 0

    #for each time compute the variance of the filtered par under the normal distrib of static par

    for k=1:nErgmPar
        for t=1:T
            for b=1:nBands
                length(quantilesVals[b]) ==2 ? () : error()

                a_t_vec = distribFilteredSD[:,k,t]

                # add filtering and parameter unc
                confBandsPar[k, t, b, :] = Statistics.quantile(a_t_vec[.! isnan.(a_t_vec)], quantilesVals[b])
            end
        end
    end

    return confBandsPar
end

function conf_bands_coverage(parDgpTIn, confBandsIn; offset=0 )

    nErgmPar, T, nBands, nQuant  = size(confBandsIn)

    T = T-offset
    parDgpT = parDgpTIn[:, 1:end-offset]
    confBands = confBandsIn[:, 1+offset:end, :, :]

    nQuant ==2 ? () : error()
    nParDgp, TDgp = size(parDgpT)
    nParDgp == nErgmPar ? () : error()
    T == TDgp ? () : error()

    isCovered = falses(nErgmPar, T, nBands)
    for b in 1:nBands
        for p in 1:nErgmPar
            for t in 1:T
                ub = confBands[p, t, b, 1]
                lb = confBands[p, t, b, 2]
                ub<lb ? error("wrong bands ordering") : ()
                isCov = lb <= parDgpT[p, t] <= ub
                isCovered[p, t, b] = isCov
            end
        end
    end
    return isCovered
end

function conf_bands_given_SD_estimates(model::SdErgm, N, obsT, vEstSdUnParAll, ftot_0, quantilesVals::Vector{Vector{Float64}}; indTvPar = model.indTvPar, parDgpT=zeros(2,2), plotFlag=false, parUncMethod = "WHITE-MLE", dropOutliers = false, offset = 1,  nSample = 100 , mvSDUnParEstCov = Symmetric(zeros(3,3)), sampleStaticUnPar = zeros(3,3), winsorProp=0, xval=nothing)

    T = length(obsT)
    nStaticPar = length(indTvPar) + 2*sum(indTvPar)
    nErgmPar = number_ergm_par(model)

    vEstSdResParAll = restrict_all_par(model, vEstSdUnParAll)
    vEstSdResPar, vConstPar = divide_SD_par_from_const(model, vEstSdResParAll)

    fVecT_filt , target_fun_val_T, sVecT_filt = score_driven_filter(model, N, obsT,  vEstSdResPar, indTvPar; ftot_0 = ftot_0, vConstPar=vConstPar)

    isDistribFilteredSDAvail = false
    # are we providing a distribution of static parameters ?
    if sampleStaticUnPar != zeros(3,3)

        errFlagEstCov = false
        errFlagSample = false

    else

        # if an external covariance matrix is given, use it to get a  static mv normal sample of stati parameters
        if mvSDUnParEstCov != Symmetric(zeros(3,3))
            errFlagEstCov = false
            Logging.@info("[conf_bands_given_SD_estimates][using external gaussian covariance matrix]")

            sampleMvNormal = true

        elseif parUncMethod[1:2] == "PB"
            @debug "[conf_bands_given_SD_estimates][PB][vEstSdResParAll=$vEstSdResParAll]"
            # Parametric Bootstrap
            distribFilteredSD, filtCovHatSample, _,  sampleStaticUnPar = par_bootstrap_distrib_filtered_par(model, N, obsT, indTvPar, ftot_0, vEstSdResParAll)
            mean(errFlagVec) > 0.2 ? errFlagEstCov=true : errFlagEstCov = false

            if contains(parUncMethod, "SAMPLE")
                @debug "[conf_bands_given_SD_estimates][PB][SAMPLE]"

                    error("need to clean the sample of parametric bootstrapped paths for outliers(e.g. static estimates)")
                isDistribFilteredSDAvail = true
                errFlagSample = false
                sampleMvNormal = false
            elseif contains(parUncMethod, "MVN")
                @debug "[conf_bands_given_SD_estimates][PB][MVN]"
                sampleMvNormal = true
                mvSDUnParEstCov = MinCovDet().fit(sampleStaticUnPar').covariance_
            end

        elseif contains(parUncMethod, "WHITE-MLE")
            # Huber- White robust estimator
            try
                @debug "[conf_bands_given_SD_estimates][WHITE-MLE][vEstSdResParAll=$vEstSdResParAll]"
                mvSDUnParEstCov, errFlagEstCov = white_estimate_cov_mat_static_sd_par(model, N, obsT, indTvPar, ftot_0, vEstSdResParAll)
                    sampleMvNormal = true
            catch
                Logging.@error("Error in computing white estimator")
                errFlagEstCov = true
            end

        elseif parUncMethod == "HESS"
                @debug "[conf_bands_given_SD_estimates][HESS][vEstSdResParAll=$vEstSdResParAll]"
            # Hessian of the objective function
            OPGradSum, HessSum = A0_B0_est_for_white_cov_mat_obj_SD_filter_time_seq(model, N, obsT, vEstSdResParAll, ftot_0)

            errFlagEstCov = false
            sampleMvNormal = true

            hess = HessSum
            mvSDUnParEstCov, minEigenVal = make_pos_def(Symmetric(pinv(hess)))

        elseif contains(parUncMethod, "NPB")
            # non parametric bootstrap
            @debug "[conf_bands_given_SD_estimates][NPB][vEstSdResParAll=$vEstSdResParAll]"

            vEstSdUnParBootDist = non_par_bootstrap_distrib_filtered_par(model, N, obsT, indTvPar, ftot_0; nBootStrap = 100)

            if contains(parUncMethod, "SAMPLE")
                @debug "[conf_bands_given_SD_estimates][NPB][SAMPLE]"
                error("need to clean the sample of non parametric bootstrapped paths for outliers(e.g. static estimates)")
                sampleMvNormal = false

            elseif contains(parUncMethod, "MVN")
                @debug "[conf_bands_given_SD_estimates][NPB][MVN]"
                sampleMvNormal = true
                mvSDUnParEstCov = MinCovDet().fit(sampleStaticUnPar').covariance_
            end

        end

        if sampleMvNormal
            @debug "[conf_bands_given_SD_estimates][samplingMvNormal]"

            if errFlagEstCov
                @debug "[conf_bands_given_SD_estimates][error in estimating the covariance matrix]"
                mvSDUnParEstCov = Symmetric(zeros(nStaticPar, nStaticPar))
                distribFilteredSD = zeros(nSample, nErgmPar,T)
                filtCovHatSample = zeros(nErgmPar, T, nSample)
                errFlagSample = true
            else
                vEstSdUnPar = unrestrict_all_par(model, vEstSdResParAll)
                mvSDUnParEstCov = Symmetric(mvSDUnParEstCov)
                @debug "covMAt =$mvSDUnParEstCov"
                sampleStaticUnPar = distrib_static_par_from_mv_normal(model, N, obsT, indTvPar, ftot_0, vEstSdUnPar, mvSDUnParEstCov; nSample = nSample)

                errFlagSample = false
            end

        end
    end

    if !isDistribFilteredSDAvail
        @debug "[conf_bands_given_SD_estimates][obtaining distribution of filtered latent]"
        distribFilteredSD, filtCovHatSample, errFlagSample = distrib_filtered_par_from_sample_static_par(model, N, obsT, indTvPar, ftot_0, sampleStaticUnPar)
    end

    # Compute confidence bands given distribution of filtered paths
    errFlagBands=false
    nBands = length(quantilesVals)
    confBandsPar = repeat(fVecT_filt, outer=(1,1,nBands, 2))
    confBandsFiltPar = repeat(fVecT_filt, outer=(1,1,nBands, 2))
    # try
        confBandsFiltPar,  confBandsPar = conf_bands_buccheri(model, N, obsT, indTvPar, fVecT_filt, distribFilteredSD, filtCovHatSample, quantilesVals, dropOutliers, winsorProp=winsorProp)
    # catch
        # errFlagBands = true
        # Logging.@error("Error in producing confidence bands")
    # end

    if plotFlag
        @debug "[conf_bands_given_SD_estimates][plotting]"
        fig, ax = plot_filtered_and_conf_bands(model, N, fVecT_filt, confBandsFiltPar ; parDgpTIn=parDgpT, nameConfBand1= "$parUncMethod - Filter+Par", nameConfBand2= "Par", confBands2In=confBandsPar, offset=offset, xval=xval)

    end

    errFlag = errFlagEstCov | errFlagSample | errFlagBands

    @debug "[conf_bands_given_SD_estimates][end]"
    fVecT_filt::Array{<:Real,2}
    confBandsFiltPar::Array{<:Real,4}
    confBandsPar::Array{<:Real,4}
    errFlag::Bool
    mvSDUnParEstCov::LinearAlgebra.Symmetric{Float64,Array{Float64,2}}
    distribFilteredSD::Array{Float64,3}

    if plotFlag
        return (; fVecT_filt, confBandsFiltPar, confBandsPar, errFlag, mvSDUnParEstCov, distribFilteredSD, fig, ax )
    else
        return (; fVecT_filt, confBandsFiltPar, confBandsPar, errFlag, mvSDUnParEstCov, distribFilteredSD )
    end
end

function estimate_filter_and_conf_bands(model::SdErgm, A_T, quantilesVals::Vector{Vector{Float64}}; indTvPar = model.indTvPar, parDgpT=zeros(2,2), plotFlag=false, parUncMethod = "WHITE-MLE",show_trace = false)

    N = size(A_T)[1]

    obsT = seq_of_obs_from_seq_of_mats(model, A_T)

    obsT, vEstSdResParAll, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0 = estimate_and_filter(model, N, obsT; indTvPar = indTvPar, show_trace = show_trace )
    Logging.@debug("[estimate_filter_and_conf_bands][vEstSdResParAll = $vEstSdResParAll] ")

    vEstSdUnPar = unrestrict_all_par(model, vEstSdResParAll)

    fVecT_filt, confBandsFiltPar, confBandsPar, errFlag, mvSDUnParEstCov, distribFilteredSD = conf_bands_given_SD_estimates(model, N, obsT, vEstSdUnPar, ftot_0, quantilesVals; indTvPar = indTvPar, parDgpT=parDgpT, plotFlag=plotFlag, parUncMethod = parUncMethod)

    return obsT, vEstSdResParAll, fVecT_filt, target_fun_val_T, sVecT_filt, conv_flag, ftot_0, confBandsFiltPar, confBandsPar, errFlag

end


text:
# ========== ACWT Utilities ==========
"""
    autocorr(f::OrthoFilter)

Generates the autocorrelation filter for a given wavelet filter.
"""
function autocorr(f::OrthoFilter)
    H = WT.qmf(f)
    l = length(H)
    result = zeros(l - 1)
    for k in 1:(l - 1)
        for i in 1:(l - k)
            @inbounds result[k] += H[i] * H[i + k]
        end
        result[k] *= 2
    end
    return result
end

"""
    pfilter(f::OrthoFilter)

Generates the high-pass autocorrelation filter

**See also:** [`qfilter`](@ref), [`autocorr`](@ref)
"""
function pfilter(f::OrthoFilter)
    a = autocorr(f)
    c1 = 1 / sqrt(2)
    c2 = c1 / 2
    b = c2 * a
    return vcat(reverse(b), c1, b)
end

"""
    qfilter(f::OrthoFilter)

Generates the low-pass autocorrelation filter.

**See also:** [`pfilter`](@ref), [`autocorr`](@ref)
"""
function qfilter(f::OrthoFilter)
    a = autocorr(f)
    c1 = 1 / sqrt(2)
    c2 = c1 / 2
    b = -c2 * a
    return vcat(reverse(b), c1, b)
end

"""
    make_acqmfpair(f::OrthoFilter)

Generates the autocorrelation quadratic mirror filters.

**See also:** [`make_acreverseqmfpair`](@ref), [`pfilter`](@ref), [`qfilter`](@ref)
"""
function make_acqmfpair(f::OrthoFilter)
    pmfilter, qmfilter = pfilter(f), qfilter(f)
    return pmfilter, qmfilter
end

"""
    make_acreverseqmfpair(f::OrthoFilter)

Generates the reverse autocorrelation quadratic mirror filters.

**See also:** [`make_acqmfpair`](@ref), [`pfilter`](@ref), [`qfilter`](@ref)
"""
function make_acreverseqmfpair(f::OrthoFilter)
    pmf, qmf = make_acqmfpair(f)
    return reverse(pmf), reverse(qmf)
end


text:
using Match

required = ["byr:" "iyr:" "eyr:" "hgt:" "hcl:" "ecl:" "pid:"]

function solve(input)
	ret = 0
	lines = collect(Iterators.takewhile(line -> line != "", input))
	while length(lines) > 0
		passport = join(lines, ' ')
		if all(contains.(passport, required))
			ret += 1
		end
		lines = collect(Iterators.takewhile(line -> line != "", input))
	end
	ret
end

function solve_with_validations(input)
	ret = 0
	lines = collect(Iterators.takewhile(line -> line != "", input))
	while length(lines) > 0
		passport = join(lines, ' ')
		fields = split(passport)
		if all(validate.(passport, required))
			ret += 1
		end
		lines = collect(Iterators.takewhile(line -> line != "", input))
	end
	ret
end

function validate(passport, req)
	@match req begin
		"byr:" => begin
			m = match(r"byr:\d{4}\b", passport)
			if m === nothing
				false
			else
				n = parse(Int16, m.match[5:end])
				1920 <= n && n <= 2002
			end
		end
		"iyr:" => begin
			m = match(r"iyr:\d{4}\b", passport)
			if m === nothing
				false
			else
				n = parse(Int16, m.match[5:end])
				2010 <= n && n <= 2020
			end
		end
		"eyr:" => begin
			m = match(r"eyr:\d{4}\b", passport)
			if m === nothing
				false
			else
				n = parse(Int16, m.match[5:end])
				2020 <= n && n <= 2030
			end
		end
		"hgt:" => begin
			m = match(r"hgt:(\d{3}cm|\d{2}in)\b", passport)
			if m === nothing
				false
			else
				if contains(m.match, "cm")
					n = parse(Int16, m.match[5:7])
					150 <= n && n <= 193
				else
					n = parse(Int16, m.match[5:6])
					59 <= n && n <= 76
				end
			end
		end
		"hcl:" => occursin(r"hcl:#[\dabcdef]{6}\b", passport)
		"ecl:" => occursin(r"ecl:(amb|blu|brn|gry|grn|hzl|oth)\b", passport)
		"pid:" => occursin(r"pid:\d{9}\b", passport)
		_ => true
	end
end

@assert solve(eachline("input/test/test_day4.txt")) == 2 ["test case 1 failed"]

println("First star: " * string(solve(eachline("input/day4.txt"))))
println("Second star: " * string(solve_with_validations(eachline("input/day4.txt"))))


text:
include("macros.jl")

using StaticArrays

# old style call
function eval_UC_spline(a, b, orders, C::Array{Float64}, S::Matrix{Float64})
    d,N = size(S)

    n_x = size(C,1)
    dims = tuple(size(C)[2:end]...)
    out = zeros(n_x,N)
    CC = reshape(reinterpret(SVector{n_x,Float64},vec(C)),dims)
    SS = reshape(reinterpret(SVector{d,Float64},vec(S)),(N,))
    V  = reshape(reinterpret(SVector{n_x,Float64},vec(out)),(N,))
    eval_UC_spline!(a, b, orders, CC, SS, V)
    return out
end

function eval_UC_spline(a, b, orders, C::Array{T, d}, S::Vector{SVector{d,Float64}}) where T where d
    N = length(S)
    vals = zeros(T,N)
    eval_UC_spline!(a, b, orders, C, S, vals)
    return vals
end

@generated function eval_UC_spline!(a, b, orders, C::Array{T, d}, S::Vector{SVector{d, Float64}}, V::Vector{T}) where d where T
    # d = C.parameters[2]-1 # first dimension of C indexes the splines
    # the speed penalty of extrapolating points when iterating over point
    # seems very small so this is the default
    fun = (create_function(d,"natural"))
    return fun.args[2].args[2]
end


text:
module SampledDomains

export CartesianDomain2D, make_centered_domain2D

#TODO rewrite as parametric type
struct CartesianDomain2D
    xrange::AbstractRange
    yrange::AbstractRange
end

function make_centered_domain2D(xlength, ylength, pixelsize)
    xrange = ((1:xlength) .- (1+xlength)/2) .* pixelsize
    yrange = ((1:ylength) .- (1+ylength)/2) .* pixelsize
    CartesianDomain2D(xrange, yrange)
end

end


text:
# BAI
# ===
#
# Index for BAM files.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# An index type for the BAM file format.
struct BAI
    # BGZF file index
    index::GenomicFeatures.Indexes.BGZFIndex

    # number of unmapped reads
    n_no_coors::Union{Nothing, Int}
end

"""
    BAI(filename::AbstractString)

Load a BAI index from `filename`.
"""
function BAI(filename::AbstractString)
    return open(read_bai, filename)
end

"""
    BAI(input::IO)

Load a BAI index from `input`.
"""
function BAI(input::IO)
    return read_bai(input)
end

# Read a BAI object from `input`.
function read_bai(input::IO)
    # check magic bytes
    B = read(input, UInt8)
    A = read(input, UInt8)
    I = read(input, UInt8)
    x = read(input, UInt8)
    if B != UInt8('B') || A != UInt8('A') || I != UInt8('I') || x != 0x01
        error("input is not a valid BAI file")
    end

    # read contents
    n_refs = read(input, Int32)
    index = GenomicFeatures.Indexes.read_bgzfindex(input, n_refs)
    if !eof(input)
        n_no_coors = read(input, UInt64)
    else
        n_no_coors = nothing
    end

    return BAI(index, n_no_coors)
end


text:
# state_updates.jl
# David Merrell
#
# A collection of functions for updating variables in a trace,
# or storing updates from MCMC.

import Base: lock, ReentrantLock
import Base.Threads: @threads

function old_update_loop(tr, proposal_param_vec::Vector{Float64})

    lambda_prop_std = proposal_param_vec[end]

    V = length(proposal_param_vec) - 1
    lambda_acc = false
    tr, lambda_acc = Gen.mh(tr, lambda_proposal, (lambda_prop_std,))

    acc_vec = zeros(V)
    for i=1:V
        tr, acc_vec[i] = Gen.mh(tr, smart_proposal,
                                (i, proposal_param_vec[i], V),
                                smart_involution)
    end

    return tr, (lambda_acc, acc_vec)
end

"""
Loop through the model variables and perform
Metropolis-Hastings updates on them.
It uses the parent set proposal distribution.
"""
function vertex_lambda_update_loop(tr, proposal_param_vec::Vector{Float64})

    lambda_prop_std = proposal_param_vec[end]

    V = length(proposal_param_vec) - 1

    acc_vec = zeros(V)
    lambda_acc = false
    for i=1:V
        tr, lambda_acc = Gen.mh(tr, lambda_vec_proposal, (i, lambda_prop_std))
        tr, acc_vec[i] = Gen.mh(tr, smart_proposal,
				(i, proposal_param_vec[i], V),
				smart_involution)
    end

    return tr, (lambda_acc, acc_vec)
end

"""
Equivalent to `vertex_lambda_update_loop`, except that it uses
multithread parallelism to update the parent sets.
"""
function multithread_update_loop(tr, proposal_param_vec::Vector{Float64})

    lambda_prop_std = proposal_param_vec[end]
    cmap = Gen.choicemap()
    r = ReentrantLock()

    V = length(proposal_param_vec) - 1

    acc_vec = zeros(Bool, V)
    lambda_acc = zeros(Bool, V)
    # Parallelize the loop over vertices
    @threads for i=1:V

        # Update the inverse temperature for this vertex
        new_tr, lambda_acc[i] = Gen.mh(tr, lambda_vec_proposal, (i, lambda_prop_std))
        if lambda_acc[i]
            lock(r) do
                cmap[:lambda_vec => i => :lambda] = new_tr[:lambda_vec => i => :lambda]
            end
        end

        # Update the parent set for this vertex
        new_tr, acc_vec[i] = Gen.mh(new_tr, smart_proposal,
				(i, proposal_param_vec[i], V),
				smart_involution)
        if acc_vec[i]
            # Update the choice map in a thread-safe fashion
            lock(r) do
                cmap[:parent_sets => i => :parents] = new_tr[:parent_sets => i => :parents]
            end
        end
    end

    tr, _, _, _ = Gen.update(tr, Gen.get_args(tr), (), cmap)
    return tr, (true, acc_vec)
end

function uniform_update_loop(tr, proposal_param_vec::Vector{Float64})

    lambda_prop_std = proposal_param_vec[end]

    V = length(proposal_param_vec) - 1

    acc_vec = zeros(V)
    lambda_acc = false
    for i=1:V
        tr, lambda_acc = Gen.mh(tr, lambda_vec_proposal, (i, lambda_prop_std))
        tr, acc_vec[i] = Gen.mh(tr, uniform_proposal, (V,), uniform_involution)
    end

    return tr, (lambda_acc, acc_vec)
end

import DataStructures: DefaultDict

"""
Aggregate summary statistics for the quantities of interest
while the markov chain runs, in an online fashion.
"""
function update_results_summary(results, tr, args)

    V = args[1]
    if results == nothing
        results = Dict("parent_sets" => [DefaultDict{Int64,Int64}(0) for i=1:V],
                       "lambda_mean" => 0.0,
                       "lambda_var" => 0.0,
                       "n" => 0
                       )
    end

    increment_counts!(results["parent_sets"], tr)

    results["lambda_mean"],
    results["lambda_var"],
    results["n"] = update_moments(results["lambda_mean"],
                                  results["lambda_var"],
                                  tr[:lambda],results["n"])

    return results
end

function update_moments(cur_mean, cur_var, new_x, cur_N)
    new_N = cur_N + 1
    diff = new_x - cur_mean
    new_mean = cur_mean + diff/new_N
    diff2 = new_x - new_mean
    if cur_N == 0
        new_var = 0.0
    else
        new_var = ((cur_N-1)*cur_var + diff*diff2)/cur_N
    end
    return new_mean, new_var, new_N
end

"""
Store a compact representation of all samples by
tracking *changes* in the model variables.
"""
function update_results_storediff(results, tr, args)

    V = args[1]
    vector_lambda = Bool(args[2])
    # Initialize the state
    if results == nothing
        results = Dict("parent_sets" => [DefaultDict{Int64,Vector{Tuple{Int64,Int64}}}([]) for i=1:V],
                       "n" => 0,
                       "prev_parents" => [Vector{Int64}() for i=1:V]
                       )
        if vector_lambda
            results["lambda"] = [Vector{Tuple{Int64,Float64}}() for i=1:V]
            for i=1:V
                push!(results["lambda"][i], (0, tr[:lambda_vec => i => :lambda]))
            end
        else
            results["lambda"] = Vector{Tuple{Int64,Float64}}()
            push!(results["lambda"], (0, tr[:lambda]))
        end

        for i=1:V
            for parent in tr[:parent_sets => i => :parents]
                push!(results["parent_sets"][i][parent], (0,1))
            end
            results["prev_parents"][i] = tr[:parent_sets => i => :parents]
        end

    else # Update the state

        n = results["n"] + 1
        # Lambda diffs
        if vector_lambda
            for i=1:V
                if results["lambda"][i][end][2] != tr[:lambda_vec => i => :lambda]
                    push!(results["lambda"][i], (n, tr[:lambda_vec => i => :lambda]))
                end
            end
        else
            if results["lambda"][end][2] != tr[:lambda]
                push!(results["lambda"], (n, tr[:lambda]))
            end
        end
        # Parent set diffs
        for i=1:V
            prev_ps = results["prev_parents"][i]
            cur_ps = tr[:parent_sets => i => :parents]
            # Has a new parent been added?
            new_parents = setdiff(cur_ps, prev_ps)
            for new_parent in new_parents
                push!(results["parent_sets"][i][new_parent], (n, 1))
            end

            # Has a parent been removed?
            removed_parents = setdiff(prev_ps, cur_ps)
            for rem_parent in removed_parents
                push!(results["parent_sets"][i][rem_parent], (n, 0))
            end

            results["prev_parents"][i] = copy(cur_ps)
        end
        # Finished. update n.
        results["n"] = n

    end

    return results
end

"""
Helper function for updating edge counts during inference.

The entry
tr[:parent_sets => i => :parents]
indicates existence of edge j --> i
(kind of backward from what you might expect.)

So this matrix should be read "row = child; column = parent"
"""
function increment_counts!(edge_counts, tr)

    for i=1:length(edge_counts) # for each vertex
        # increment the parents
        ps = tr[:parent_sets => i => :parents]
        for p in ps
            edge_counts[i][p] += 1
	end
    end
end

max_stars_repo_name: gitter-lab/ssps

max_issues_repo_name: gitter-lab/ssps

max_forks_repo_name: gitter-lab/ssps


text:
"""
Roulettewheel `IndividualsSelector`, with chance relative to fitness.

"""
struct RouletteWheelSelector <: IndividualsSelector end

#function select(::RouletteWheelSelector, population::FitPopulation, n::Integer)
#    #TODO implement
#end


text:
# Advent of Code, Day 20
#
# Transform an image by successively replacing pixels with values
# looked up from a translation table, the index being the value
# of the 9 cells surrounding each pixel, converted from binary to
# decimal. Much complicated by the fact that the input data has a
# 1 in the first position of the translation table, meaning that
# empty areas are filled with 1s, which muck up the pixel count.
#
# AK, 20/12/2021

# Read data into list of rows
f = "sample.txt"
f = "input.txt"
d = split(strip(read(f, String)), "\n")

# First row is "algorithm", convert # to 1 and . to 0
code = map(x -> x == '#' ? 1 : 0, [c for c in d[1]])
#println(code)

# Remaining rows are the image, convert to matrix
rows = map(r -> [c == '#' ? 1 : 0 for c in r], d[3:end])
image = Matrix(transpose(reshape(vcat(rows...), (:, length(rows)))))

# Transform an image by iteratively applying translation function
function transform(img, iterations)

    for iter in 1:iterations

        # If the code starts with a 1, it causes the whole image to
        # flip on each second iteration. To undo this, we neeed to
        # fill the border with 0 on the first iteration, 1 on the
        # second iteration, and so on.
        if code[1] == 1 && iseven(iter)
            bfill = 1
        else
            bfill = 0
        end

        # Enhance the image and display number of on bits
        img = enhance(img, bfill)
        #display(img)
        println("Iteration ", iter, ": bits set = ", sum(img))

    end
end

# Enhance the image, by converting 3x3 areas around a pixel
# to 9-digit binary number, converting this to decimal, and
# then looking up replacement bit in the "algorithm" vector
function enhance(img, padWith)

    # Make a copy of the image with a 2-thick border
    # around it, with the value of the padding alternating
    # to get around the flashing image problem when there
    # is a 1 at the beginning of the translation code.
    padded = addPadding(img, 2, padWith)

    # Make an empty array for the output, one bigger
    # on each side then the original image
    nrows, ncols = size(img)
    out = zeros(Int, nrows+2, ncols+2)

    # Go though the original image, and replace pixels
    # based on algorithm lookup
    for r in 1:nrows+2
        for c in 1:ncols+2
            n = numAt(padded, r+1, c+1)
            out[r,c] = code[Int(n+1)]
        end
    end

    #Return the enhanced  image
    return out

end

# Return the "number" centered at the given coordinates, i.e.,
# take the 3x3 cells around that number, string them along
# row-wise, and convert that binary representation to a decimal
# number
function numAt(img, r, c)

    # Get the 3x3 bits around the given coordinates
    # and turn them into a vector
    bits = reshape(transpose(img[r-1:r+1,c-1:c+1]), 1, :)

    # Convert this binary sequence to a decimal number
    #return bin2dec(bits)
    n = 0
    for b in bits
        n = n * 2 + b
    end
    return n

end

# Add a given number of layer of padding around a matrix,
# using the provide fill value (1/0)
function addPadding(img, n, padWith)

    # Add padding on top and bottom
    nr, nc = size(img)
    v = padWith == 1 ? ones(Int, n, nc) : zeros(Int, n, nc)
    img = vcat(v, img, v)

    # Add padding on sides
    nr, nc = size(img)
    v = padWith == 1 ? ones(Int, nr, n) : zeros(Int, nr, n)
    img = hcat(v, img, v)

    return img

end

# Perform the two transformations (Part 1 has 2 iterations, Part 2 has 50).
# For sample data, should be 35 and 3351, respectively
println("PART 1:")
transform(image, 2)

println("\nPART 2:")
transform(image, 50)


text:
include("spline_kernel.jl")

# Removing t = 0, such that Σ is invertible
t = Vector(0.1:0.1:1)
n = length(t);

# Creating a test matrix Σ = tril(UV') + triu(VU',1) that is PSD
p = 2;
U, V = spline_kernel(t, p);
Σ    = spline_kernel_matrix(U, V) + I;
chol = cholesky(Σ)

# Creating a symmetric exended generator representable semiseperable matrix
K  = DiaSymSemiseparable(U,V,ones(n))
# Calculating its Cholesky factorization
L = DiaSymSemiseparableChol(K)
# Creating a test vector
x = randn(n);

# Testing multiplication
@test isapprox(L*x, chol.L*x, atol=1e-6)
@test isapprox(L'*x, chol.U*x, atol=1e-6)

# Testing inverses (Using Cholesky factorizations)
@test isapprox(L'\(L\x), chol.U\(chol.L\x), atol=1e-6)
@test isapprox(inv(L), chol.U\(chol.L\Diagonal(ones(K.n))), atol=1e-6)
@test isapprox(inv(L,x), chol.U\(chol.L\x), atol=1e-6)
@test isapprox(K*(K\x), x)
@test isapprox(K'*(K'\x), x)

# Testing logdet
@test isapprox(logdet(L), logdet(chol.L), atol=1e-10)
@test isapprox(logdet(L'), logdet(chol.U), atol=1e-10)

# Testing traces
@test isapprox(tr(L), tr(chol.L)) # Trace of L
@test isapprox(trinv(L), tr(chol\Diagonal(ones(K.n)))) # Trace of (K+D)^(-1)
D = SymSemiseparable(U,V);
@test isapprox(tr(L,D), tr(chol\spline_kernel_matrix(U, V))) # Trace of (K+D)^(-1)K

# Testing using the SymSemiseparableChol struct as a lower triangular matrix
C = DiaSymSemiseparableChol(U,V,ones(K.n));
@test isapprox(C*x, (tril(U*V',-1) + I)*x)


text:
function reinterpret(a::Array{V}) where V<:MultiValue{S,T,N,L} where {S,T,N,L}
  b = reinterpret(T,a)
  sa = size(a)
  s = (L,sa...)
  reshape(b,s)
end


text:
export PerformanceEvent, Performance
export getnotesequence, setlength!

using ..NoteSequences: SeqNote, DEFAULT_QPM, DEFAULT_TPQ
using ..NoteSequences: MIN_MIDI_VELOCITY, MAX_MIDI_VELOCITY, MIN_MIDI_PITCH, MAX_MIDI_PITCH
using MacroTools: @forward

"""
    PerformanceEvent <: Any

Event-based performance representation of music data from Oore et al.
In this representation, midi data is encoded as a 388-length one hot vector.
It corresponds to NOTE_ON, NOTE_OFF events for each of the 128 midi pitches,
32 bins for the 128 midi velocities and 100 TIME_SHIFT events
(it can represent a time shift from 10 ms up to 1 second).
`PerformanceEvent` is the base representation.

## Fields
* `event_type::Int`  :  Type of the event. One of {NOTE_ON, NOTE_OFF, TIME_SHIFT, VELOCITY}.
* `event_value::Int` :  Value of the event corresponding to its type.
"""
struct PerformanceEvent
    event_type::Int
    event_value::Int

    function PerformanceEvent(type::Int, value::Int)
        if (type == NOTE_ON || type == NOTE_OFF)
            if !(MIN_MIDI_PITCH <= value <= MAX_MIDI_PITCH)
                error("The event has an invalid pitch value")
            end
        elseif type == TIME_SHIFT
            if !(value >= 0)
                error("The event has an invalid time shift value")
            end
        elseif type == VELOCITY
            if !(MIN_MIDI_VELOCITY <= value <= MAX_MIDI_VELOCITY)
                error("The event has an invalid velocity value")
            end
        else
            error("PerformanceEvent has invalid type")
        end

        new(type, value)
    end
end

function Base.show(io::IO, a::PerformanceEvent)
    if a.event_type == NOTE_ON
        s = "NOTE-ON"
    elseif a.event_type == NOTE_OFF
        s = "NOTE-OFF"
    elseif a.event_type == TIME_SHIFT
        s = "TIME-SHIFT"
    elseif a.event_type == VELOCITY
        s = "VELOCITY"
    else
        error("PerformanceEvent has invalid type")
    end

    s *= " $(a.event_value)"
    print(io, s)
end

"""
    Performance <: Any

`Performance` represents a polyphonic sequence
as a stream of `PerformanceEvent`s.

See [`PerformanceEvent`](@ref)

## Fields
* `events::Vector{PerformanceEvent}`: The performance in a vector of `PerformanceEvent`s.
* `program::Int`: MIDI program number for this `Performance`
* `startstep::Int`: The start step of the performance relative to the source sequence.
* `velocity_bins::Int`: Number of bins for velocity values.
* `steps_per_second::Int`: Number of steps per second for quantization.
* `max_shift_steps::Int`: Maximum number of steps in a `TIME_SHIFT` event.
"""
mutable struct Performance
    events::Vector{PerformanceEvent}
    program::Int
    startstep::Int
    velocity_bins::Int
    steps_per_second::Int
    max_shift_steps::Int

    function Performance(startstep::Int,
                         velocity_bins::Int,
                         steps_per_second::Int,
                         max_shift_steps::Int;
                         program::Int=-1,
                         events::Vector{PerformanceEvent}=Vector{PerformanceEvent}())

        new(events,
            program,
            startstep,
            velocity_bins,
            steps_per_second,
            max_shift_steps)
    end
end

function Performance(quantizedns::NoteSequence;
                     startstep::Int = 0,
                     velocity_bins::Int = 0,
                     max_shift_steps::Int = DEFAULT_MAX_SHIFT_STEPS,
                     instrument::Int = -1)

    if !quantizedns.isquantized
        throw(ArgumentError("The `NoteSequence` must be quantized."))
    end

    steps_per_second = quantizedns.sps
    events = getperfevents(quantizedns, startstep, velocity_bins, max_shift_steps, instrument)

    Performance(startstep, velocity_bins, steps_per_second, max_shift_steps, events=events)
end

function Performance(steps_per_second::Int;
                     startstep::Int = 0,
                     velocity_bins::Int = 0,
                     max_shift_steps::Int = DEFAULT_MAX_SHIFT_STEPS,
                     program::Int = -1)

    Performance(startstep, velocity_bins, steps_per_second, max_shift_steps, program=program)
end

function Base.getproperty(performance::Performance, sym::Symbol)
    if sym === :numsteps
        steps = 0
        for event in performance
            if event.event_type == TIME_SHIFT
                steps += event.event_value
            end
        end
        return steps
    else
        getfield(performance, sym)
    end
end

@forward Performance.events Base.length, Base.getindex, Base.lastindex, Base.pop!

Base.setindex!(p::Performance, event::PerformanceEvent, idx::Int) = setindex!(p.events, event, idx)
Base.iterate(p::Performance, state=1) = iterate(p.events, state)
Base.push!(p::Performance, event::PerformanceEvent) = push!(p.events, event)
Base.append!(p::Performance, events::Vector{PerformanceEvent}) = append!(p.events, events)

function Base.append!(p1::Performance, p2::Performance)
    append!(p1, p2.events)
end

function Base.copy(p::Performance)
    Performance(
        copy(p.events),
        p.program,
        p.startstep,
        p.velocity_bins,
        p.steps_per_second,
        p.max_shift_steps)
end

"""
    getperfevents(quantizedns::NoteSequence, startstep::Int, velocity_bins::Int,
                  max_shift_steps::Int, instrument::Int=-1)

Extract performance events from a quantized [`NoteSequence`](@ref).

The extraction starts from `startstep`. The number of velocity bins to use
is given by `velocity_bins`. If it's 0, velocity events will not be included.
`max_shift_steps` is the maximum number of steps for a single `TIME_SHIFT` event.

If `instrument` is -1, performance events will be extracted for all instruments in the `NoteSequence`
Otherwise, performance events will be extracted only for the given instrument.
"""
function getperfevents(quantizedns::NoteSequence,
                       startstep::Int,
                       velocity_bins::Int,
                       max_shift_steps::Int,
                       instrument::Int = -1)

    if !quantizedns.isquantized
        throw(ArgumentError("The `NoteSequence` is not quantized."))
    end

    notes = [note for note in quantizedns.notes
             if note.start_time >= startstep &&
             (instrument == -1 || note.instrument == instrument)]

    sort!(notes, by=note->(note.start_time, note.pitch))

    onsets = [(note.start_time, idx, false) for (idx, note) in enumerate(notes)]
    offsets = [(note.end_time, idx, true) for (idx, note) in enumerate(notes)]
    noteevents = sort(vcat(onsets, offsets))

    currentstep = startstep
    currentvelocitybin = 0
    perfevents = Vector{PerformanceEvent}()

    for (step, idx, isoffset) in noteevents
        if step > currentstep
            while step > (currentstep + max_shift_steps)
                push!(perfevents, PerformanceEvent(TIME_SHIFT, max_shift_steps))
                currentstep += max_shift_steps
            end
            push!(perfevents, PerformanceEvent(TIME_SHIFT, step - currentstep))
            currentstep = step
        end

        if velocity_bins > 0
            velocity = velocity2bin(notes[idx].velocity, velocity_bins)
            if !isoffset && velocity != currentvelocitybin
                currentvelocitybin = velocity
                push!(perfevents, PerformanceEvent(VELOCITY, currentvelocitybin))
            end
        end

        push!(perfevents, PerformanceEvent(ifelse(isoffset, NOTE_OFF, NOTE_ON), notes[idx].pitch))
    end

    perfevents
end

"""
    getnotesequence(performance::Performance, velocity::Int=100, instrument::Int=1, program::Int=-1)

Return a `NoteSequence` from the `Performance`.

All of the notes in the `NoteSequence` will have the same `velocity`, `program` number
and `instrument` number.
"""
function getnotesequence(performance::Performance,
                         velocity::Int = 100,
                         instrument::Int = 1,
                         program::Int = -1)

    ticks_per_step = second2tick(1) / performance.steps_per_second
    ns = tosequence(performance, ticks_per_step, velocity, instrument, program)
    push!(ns.tempos, NoteSequences.Tempo(0, 120.0))
    push!(ns.timesignatures, MIDI.TimeSignatureEvent(0, 4, 4, 24, 8))

    ns
end

function tosequence(performance::Performance,
                    ticks_per_step::Float64,
                    velocity::Int,
                    instrument::Int,
                    program::Int)

    sequence = NoteSequence(DEFAULT_TPQ, isquantized=false)
    seqstart_time = performance.startstep * ticks_per_step

    DEFAULT_PROGRAM = 0
    if program == -1
        program = ifelse(performance.program != -1, performance.program, DEFAULT_PROGRAM)
    end

    # Map the pitch of a note to a list of start steps and velocities
    # (since may be active multiple times)
    pitchmap = DefaultDict{Int, Vector{Tuple{Int, Int}}}(Vector{Tuple{Int, Int}})
    step = 0

    for event in performance
        if event.event_type == NOTE_ON
            push!(pitchmap[event.event_value], (step, velocity))
        elseif event.event_type == NOTE_OFF
            if event.event_value in keys(pitchmap)
                pitchstartstep, pitchvelocity = popfirst!(pitchmap[event.event_value])

                if isempty(pitchmap[event.event_value])
                    delete!(pitchmap, event.event_value)
                end

                # If start step and end step are the same, ignore
                if step == pitchstartstep
                    continue
                end

                start_time = round(ticks_per_step * pitchstartstep) + seqstart_time
                end_time = round(ticks_per_step * step) + seqstart_time
                note = SeqNote(event.event_value, pitchvelocity, start_time, end_time, instrument, program)
                push!(sequence.notes, note)
                if note.end_time > sequence.total_time
                    sequence.total_time = note.end_time
                end
            end
        elseif event.event_type == TIME_SHIFT
            step += event.event_value
        elseif event.event_type == VELOCITY
            if event.event_value != velocity
                velocity = bin2velocity(event.event_value, performance.velocity_bins)
            end
        else
            throw(error("Unkown event type $(event.event_type)"))
        end
    end

    # End all the notes that weren't ended
    for pitch in keys(pitchmap)
        for (pitchstartstep, pitchvelocity) in pitchmap[pitch]
            if step == pitchstartstep
                continue
            end

            start_time = round(ticks_per_step * pitchstartstep) + seqstart_time
            end_time = round(ticks_per_step * step) + seqstart_time
            # Maybe end after 5 seconds?
            # end_time = start_time + round(ticks_per_step * 5 * performance.steps_per_second)
            note = SeqNote(pitch, pitchvelocity, start_time, end_time, instrument, program)
            push!(sequence.notes, note)
            if note.end_time > sequence.total_time
                sequence.total_time = note.end_time
            end
        end
    end

    sequence
end

"""
    truncate(performance::Performance, numevents::Int)

Truncate the performance to exactly `numevents` events.
"""
function Base.truncate(performance::Performance, numevents::Int)
    performance.events = performance.events[1:numevents]
end

function append_steps(performance::Performance, numsteps::Int)
    max_shift_steps = performance.max_shift_steps
    if (!isempty(performance) &&
        performance[end].event_type == TIME_SHIFT &&
        performance[end].event_value < max_shift_steps)

        steps = min(numsteps, max_shift_steps - performance[end].event_value)
        performance[end] = PerformanceEvent(TIME_SHIFT, performance[end].event_value + steps)
        numsteps -= steps
    end

    while numsteps >= max_shift_steps
        push!(performance, PerformanceEvent(TIME_SHIFT, max_shift_steps))
        numsteps -= max_shift_steps
    end

    if numsteps > 0
        push!(performance, PerformanceEvent(TIME_SHIFT, numsteps))
    end
end

function trim_steps(performance::Performance, numsteps::Int)
    trimmed = 0
    while !isempty(performance) && trimmed < numsteps
        if performance[end].event_type == TIME_SHIFT
            if trimmed + performance[end].event_value > numsteps
                performance[end] = PerformanceEvent(TIME_SHIFT, performance[end].event_value - numsteps + trimmed)
                trimmed = numsteps
            else
                trimmed += performance[end].event_value
                pop!(performance)
            end
        else
            pop!(performance)
        end
    end
end

"""
    setlength!(performance::Performance, steps::Int)

Sets the total length of the `Performance` based on the length given by `steps`.

If the length of the performance is greater than the given steps, it is trimmed.
Otherwise, the performance is padded with `TIME_SHIFT`s.
"""
function setlength!(performance::Performance, steps::Int)
    if performance.numsteps < steps
        append_steps(performance, steps - performance.numsteps)
    elseif performance.numsteps > steps
        trim_steps(performance, performance.numsteps - steps)
    end

    @assert performance.numsteps == steps
end

"""     binsize(velocity_bins)
Returns the size of a bin given the total number of bins.
"""
binsize(velocity_bins::Int) = Int(ceil(
        (MAX_MIDI_VELOCITY - MIN_MIDI_VELOCITY + 1) / velocity_bins))

velocity2bin(velocity::Int, velocity_bins::Int) = ((velocity - MIN_MIDI_VELOCITY) ÷ binsize(velocity_bins)) + 1

"""     bin2velocity(bin, velocity_bins)
Returns a velocity value given a bin and the total number of velocity bins.
"""
bin2velocity(bin::Int, velocity_bins::Int) = MIN_MIDI_VELOCITY + (bin - 1) * binsize(velocity_bins)


text:
export ARCSpectral_abs

function ARCSpectral_abs(nlp 		:: AbstractNLPModel,
               		     nlpstop 	:: NLPStopping;
					     kwargs...
               		     )

	return TRARC2(nlp,
				  nlpstop;
				  TR = TrustRegion(10.0),
				  c = Combi(ARCTR.hessian_dense, ARCTR.PDataSpectral{eltype(nlp.meta.x0)}, ARCTR.solve_modelARCDiagAbs, ARCTR.preprocessSpectral, ARCTR.decreaseFact, ARCTR.Tparam{eltype(nlp.meta.x0)}()),
				  kwargs...
				  )
end


text:
####################################################################
# Using Meek's Rules to Update PDAG
####################################################################

#Revert a graph to undirected edges and unshielded colliders (i.e. parents not adjacent)
function graphVStructure!(g::PDAG)

    #loop through all vertices
    for x in vertices(g)
        #get the parents of current vertex
        parentsX = parents(g,x)
        #if all the parents are adjacent, undirect the edges
        #(if there is only 1 parent, then it is still a clique)
        if isclique(g, parentsX)
            for p in parentsX
                addedge!(g, x, p, directed=false)
            end
        end
    end
end

#Thoughts on improvement
    #Would it be cleaner to find all unique triples where one edge is undirected?
    #Can this update be more local? Most edges will remain unchanged.

#This is set up so we only loop through all the vertices and edges once
function meekRules!(g::PDAG)

    #Loop through all the edges in the graph (x-y)
    for edge in edges(g)
        #We only need to update undirected edges
        if !edge.directed

            #For clarity extract the edge vertices
            x, y = edge.parent, edge.child

            #Label the neighbors of the nodes that comprise the current edge
            #Neighbors can be labelled "parent", "child", or "undirected"
            #e.g. v₁→x-v₂, v₁ is a parent and v₂ is undirected
            xNeighbors, yNeighbors =  categorizeNeighbors(g, x, y)

            #Rule 1: If x or y have a unique parent, direct edge away from it
            #Check if x has a unique parent
            if R1(xNeighbors, yNeighbors)
                #Add x→y
                addedge!(g, x, y)
            #Check if y has a unique parent
            elseif R1(yNeighbors, xNeighbors)
                #Add y→x
                addedge!(g, y, x)
            end

            #Rule 2: Direct the edge away from a potential cycle
            #Check if x has a child that is a parent of y
            if R2(xNeighbors, yNeighbors)
                #Add x→y
                addedge!(g, x, y)
            #Check if y has a child that is a parent of x
            elseif R2(yNeighbors, xNeighbors)
                #Add y→x
                addedge!(g, y, x)
            end

            #Rule 3: Double Triangle, diagonal
            if R3(xNeighbors, yNeighbors)
                #Add x→y
                addedge!(g, x, y)
            elseif R3(yNeighbors, xNeighbors)
                #Add y→x
                addedge!(g, y, x)
            end

            #Rule 4: Double Triangle, side
            #Requires an additional test that can't be check with the neighbor dictionaries
            #So, we need to pass the graph through as well
            if R4(xNeighbors, yNeighbors, g)
                #Add x→y
                addedge!(g, x, y)
            elseif R4(yNeighbors, xNeighbors, g)
                #Add y→x
                addedge!(g, y, x)
            end

        end
    end

    return nothing
end

function R1(neighborSet1, neighborSet2)
    #given x-y, look for patterns that match v₁→x and not(v₁→y)
    for (v₁, category) in neighborSet1
        if category == :parent && !haskey(neighborSet2,v₁)
            return true
        end
    end
    return false
end

function R2(neighborSet1, neighborSet2)
    #given x-y, look for patterns that match x→v₁→y
    for (v₁, category) in neighborSet1
        if category == :child && haskey(neighborSet2,v₁) && neighborSet2[v₁]==:parent
            return true
        end
    end
    return false
end

function R3(neighborSet1, neighborSet2)

    #given x-y, count the number of patterns that match x-v₁→y
    counter = 0
    for (v₁, category) in neighborSet1
        if category == :undirected && haskey(neighborSet2,v₁) && neighborSet2[v₁]==:parent
            counter += 1
            #If we find two such paths, the Rule 3 pattern is matched
            if counter == 2
                return true
            end
        end
    end

    return false
end

function R4(neighborSet1, neighborSet2, g)

    #given x-y, look for patterns that match x-v₁→v₂→y and x-v₂
    for (v₂, category) in neighborSet1
        #first look for x-v₂→y
        if category == :undirected && haskey(neighborSet2,v₂) && neighborSet2[v₂]==:parent
            #check for x-v₁ and v₁→v₂
            for (v₁, category) in neighborSet1
                if category == :undirected && isparent(g, v₁, v₂)
                    return true
                end
            end
        end
    end

    return false
end

#Categorize neighboring edges as "parent", "child", or "undirected"
function categorizeNeighbors(g::PDAG, x, y)

    #Create a dictionary of neighbors (e.g. vertex 5 => :parent)
    xNeighbors = Dict{Int, Symbol}()
    yNeighbors = Dict{Int, Symbol}()

    #Loop through all the vertices in the graph
    for vᵢ in vertices(g)

        #If vᵢ is adjacent, give it a category
        if isadjacent(g, x, vᵢ)
            xNeighbors[vᵢ] = setCategory(g, vᵢ, x)
        end

        #Same procedure to the other vertex in the current edge
        if isadjacent(g, y, vᵢ)
            yNeighbors[vᵢ] = setCategory(g, vᵢ, y)
        end
    end

    return (xNeighbors, yNeighbors)
end

function setCategory(g::PDAG, v₁, v₂)

    #Test if v₁ → v₂
    if isparent(g,v₁,v₂)
        return :parent
    #Test if v₁ ← v₂
    elseif ischild(g,v₁,v₂)
        return :child
    #If not then we must have an undirected edge
    else
        return :undirected
    end
end


text:
########################################################################################################
#HERE STARTS THE FUNCTIONS FOR THE HUGGETT MODEL
#######################################################################################################

#find q using bisection (thats why I need update params here - its lbub vector (upper and lower bound))
function UpdateGEq_Hugget_bisection(Params,SolParam_old, ass_sum,lbub, iter)

    lb = lbub[1]
    ub = lbub[2]

    if(iter ==1 )
        ub = Params[7]
        return(ub), [lb, ub]
    elseif iter==2
        ub = Params[7]
        lb = Params[6]
        return  (ub+lb)/2,[lb, ub]
    else
        if(ass_sum<=0.0)
            ub = SolParam_old[iter]
            return (ub+lb)/2,[lb, ub]
        else
            lb = SolParam_old[iter]
            return (ub+lb)/2,[lb, ub]
        end
    end

end

function SaveHAEq_Hugget(Policy,Distr, q)
    #ok this might be better, for now just plot policy functions and print q (but q is ok)
    policy_c = Policy[1]
    plot(Policy[3], policy_c[1,:], label = "income =1.0")
    plot!(Policy[3], policy_c[2,:], label = "income =0.1")
    png("plot_policies")
    println(" equilibrium q is ", q)
end
#check convergence
function ConvSolParam_Hugget_bisection(SolParam_old, SolParam_new, ass_sum, maxiter, iter)
    return abs(ass_sum)>1e-6 && iter<=maxiter
end
#

function SolveAgP_Huggett_EGM(Params, NumParams, q; maxiter=2000)
    """
    Find consumption and asset policies using Endogenous grid method
    """
    #read Params
    β = Params[1]
    σ = Params[2]
    ϵ_vals = Params[3]
    P = Params[4]
    ϵ_n = size(P,1)

    a_n = NumParams[1]
    a_grid  = NumParams[2]
    a_min = NumParams[3]

    #initialize policies
    policy_c_old = zeros(ϵ_n,a_n)
    policy_c_new = zeros(ϵ_n,a_n)
    policy_a = zeros(ϵ_n,a_n)
    EGrid = zeros(a_n)
    c_prev_vals = zeros(a_n)

    #initialize consumption function
    for j in 1:ϵ_n
        for i in 1:a_n
            policy_c_old[j,i] = max(a_grid[i]+ϵ_vals[j],1e-7)
        end
    end

    iter = 0

    #convergence check
    con_measure = ones(ϵ_n,a_n)

    #ITERATIONS STARTS HERE
    while( maximum(abs.(con_measure))>=1e-8 && iter<=maxiter )
        if(iter!=0)
            policy_c_old = copy(policy_c_new)
        end
        for j in 1:ϵ_n
            a_prev_old = a_min+0.001
            for i in 1:a_n
                Futur_cons = 0.0
                for jj in 1:ϵ_n
                    Futur_cons = Futur_cons+P[j,jj]*max(policy_c_old[jj,i], 1e-8)^(-σ)
                end
                #calculate todays consumption if future assets are given by gird
                c_prev = max((β/q* Futur_cons)^(-1/σ), 1e-7)
                a_prev = maximum([c_prev+q*a_grid[i]-ϵ_vals[j], a_min])

                #update endogenous grid
                EGrid[i] = a_prev

            end

            #linearly interpolate, to get policy function of exogenous grid
            Policy_EGM = LinearInterpolation(EGrid, a_grid, extrapolation_bc = Flat())
            #find next period assets and consumption, updating policy functions
            policy_a[j,:] = Policy_EGM.(a_grid)
            policy_c_new[j,:] = a_grid.+ϵ_vals[j] - q.*policy_a[j,:]

        end
        #check if consumption function converged
        con_measure = policy_c_old-policy_c_new
        iter =iter+1

    end

    return (policy_c_new, policy_a, a_grid)
end

#Additional function for monte calrlo simulation, might be helpful , for now not used
function MarkovChain_sim(start_val, M, P)
    ϵ_n = size(P,1)
    rand_choices = rand(M)
    P_increasing = zeros(ϵ_n,ϵ_n)
    P_increasing[:,1] = P[:,1]
    shock_history = zeros(M)
    for j = 1: ϵ_n
        for jj = 2:ϵ_n
            P_increasing[j,jj]=  P_increasing[j,jj-1]+P[j,jj]
        end
    end

    for i  in 1:M
        sim_val = 0
        for jj in 1:ϵ_n
            if P_increasing[start_val,jj] >= rand_choices[i]
                shock_history[i] = jj
                break
            end

        end
        if sim_val == 0
            shock_history[i] = ϵ_n
        end
        start_val = shock_history[i]
    end
    return shock_history
end

function SolveDistr_Hugget_Iter(Params, NumParams,q, Policy; maxiter = 1000 )
    """
    This function finds the distribution using the updaing over the grid (instead of Monte Carlo). I just start with some initial probability measure, update it for all possible events (rounding the future assets to grid points)
    is preatty fast, but might not work that well with huge state space

    """
    #read params
    ϵ_vals = Params[3]
    P = Params[4]
    ϵ_n = size(P,1)

    a_n = NumParams[1]
    a_grid  = NumParams[2]
    a_min = NumParams[3]

    #read policies
    policy_c = Policy[1]
    #here I build array of interpolations to not construct it each time
    Policy_func_c =  Array{Interpolations.Extrapolation}(undef,ϵ_n)
    for j in 1:ϵ_n
        Policy_func_c[j] = LinearInterpolation(a_grid, policy_c[j,:] , extrapolation_bc = Flat())
    end

    #define starting value
    start_ϵ = 1
    start_a = searchsortedlast(a_grid, 0.0)

    #initialize measures
    measure_old = zeros(ϵ_n,a_n)
    measure_new = zeros(ϵ_n,a_n)
    measure_new[start_ϵ, start_a] = 1.0
    #convergence params
    con_measure =ones(ϵ_n, a_n)

    iter =0

    #ITERATIONS STARTS HERE
    while(maximum(abs.(con_measure))>=1e-16 && iter<=maxiter)
        measure_old = copy(measure_new)
        measure_new = zeros(ϵ_n,a_n)
        for j in 1:ϵ_n
            for i in 1:a_n
                if measure_old[j,i]>0.0
                    c = Policy_func_c[j](a_grid[i])
                    #find new asset
                    a_new = 1/q*(ϵ_vals[j] + a_grid[i]-c)
                    #Now define the nearest neighbords on the gird
                    a_new_ongrid_1 = searchsortedlast(a_grid, a_new)
                    a_new_ongrid_2 = min(a_new_ongrid_1+1, a_n)
                    if a_new_ongrid_1 ==0
                        a_new_ongrid_1 = 1
                        weight = 1.0
                    elseif a_new_ongrid_1==a_new_ongrid_2
                        weight = 1.0
                    else
                        weight = (a_new-a_grid[a_new_ongrid_1])/(a_grid[a_new_ongrid_2]-a_grid[a_new_ongrid_1])
                    end
                    #Update measure
                    for jj in 1:ϵ_n

                        measure_new[jj,a_new_ongrid_1] =  P[j,jj]*(1.0-weight)*measure_old[j,i] + measure_new[jj,a_new_ongrid_1]

                        measure_new[jj,a_new_ongrid_2] = P[j,jj]*weight* measure_old[j,i]+measure_new[jj,a_new_ongrid_2]
                    end

                end
            end

        end
        #check convergence
        con_measure = measure_new-measure_old
        iter =iter+1
    end
    Asset_vals = zeros(ϵ_n, a_n)
    #find sum of assets in the economy (in equilibrium 0), should be
    for j = 1:ϵ_n
        Asset_vals[j,:] = a_grid #1/q*(a_grid.-Policy_func_c[j].(a_grid).+ϵ_vals[j])
    end

    return measure_new, sum(Asset_vals.*measure_new)
end

function Find_eq_Hugget(Params, NumParams)
    """
    These functions solve hugget model for the parameters tuple, using the genral function SolveHAEq and defined above specific functions for Huggett model

    """
    start_q = Params[6] #it should be lower bound
    flag_init = 10.0 #this for sure is not 0, so loop will not stop at first iteration
    #params for bisection
    lbub = [Params[6], Params[7]]
    #Now add the special functions for Hugget model to the general function
    SolveHAEq(Params, NumParams, start_q, SolveAgP_Huggett_EGM, SolveDistr_Hugget_Iter, UpdateGEq_Hugget_bisection, ConvSolParam_Hugget_bisection, SaveHAEq_Hugget, flag_init,lbub)
    #Done if we done
    println("DONE")

end


text:
using Documenter, POMCP

makedocs(modules=[POMCP])

deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/JuliaPOMDP/POMCP.jl.git",
    julia = "0.5",
    osname = "linux"
)


text:
println("F1) Tests of the semi-empirical computation of ADK rates.")

@warn("\n\n !!! This example does not work properly at present !!! \n\n")
#
wa = 1