text:
"""
$(SIGNATURES)

Move node at `index` to `new_location`.
"""
@propagate_inbounds function assign(sa_struct::SAStruct, index, new_location)
    node = sa_struct.nodes[index]
    # Assign this location to the new node.
    assign(node, new_location)
    # Update the grid to point to the node.
    return sa_struct.grid[location(node)] = index
end

@propagate_inbounds function unsafe_assign(sa_struct::SAStruct, node::SANode, new_location)
    return assign(node, new_location)
end

"""
$(SIGNATURES)

Move node at `index` from its current location to `new_location`.
"""
@propagate_inbounds function move(sa_struct::SAStruct, index, new_location)
    node = sa_struct.nodes[index]
    # Zero out this node's original location in the grid, then assign the new
    # location.
    sa_struct.grid[location(node)] = 0
    sa_struct.grid[new_location] = index
    return unsafe_assign(sa_struct, node, new_location)
end

"""
$(SIGNATURES)

Swap the locations of two nodes with indices `node1_idx` and `node2_idx`.
"""
@propagate_inbounds function swap(sa_struct::SAStruct, node1_idx, node2_idx)
    # Get references to these objects to make life easier.
    n1 = sa_struct.nodes[node1_idx]
    n2 = sa_struct.nodes[node2_idx]
    # Swap address/component assignments
    s = location(n1)
    t = location(n2)

    unsafe_assign(sa_struct, n1, t)
    unsafe_assign(sa_struct, n2, s)

    # Swap grid.
    sa_struct.grid[t] = node1_idx
    sa_struct.grid[s] = node2_idx

    return nothing
end

################################################################################
# DEFAULT METRIC FUNCTIONS
################################################################################
#=
Right now, this first function is basically a dispatch method to handle the
case when we have multiple channel types.

Ideas include:

1. Allow type dispatch to do everything.
2. Make this a "generated" function based on all the types of the channels and
    build a "case" statement out of it.
=#
@propagate_inbounds function channel_cost(sa_struct::SAStruct, idx::Int)
    return channel_cost(sa_struct, sa_struct.channels[idx])
end

"""
    channel_cost(sa_struct{T}, channel :: SAChannel) where {T <: RuleSet}

Return the cost of `channel`. Default implementation accumulates the distances
between the source and sink addresses of the channel using `sa_struct.distance`.

Method List
-----------
$(METHODLIST)
"""
@propagate_inbounds function channel_cost(sa_struct::SAStruct, channel::TwoChannel)
    a = sa_struct.nodes[channel.source]
    b = sa_struct.nodes[channel.sink]
    return Float64(getdistance(sa_struct.distance, a, b))
end

@propagate_inbounds function channel_cost(sa_struct::SAStruct, channel::MultiChannel)
    cost = 0.0
    for src in channel.sources, snk in channel.sinks
        a = sa_struct.nodes[src]
        b = sa_struct.nodes[snk]
        cost += getdistance(sa_struct.distance, a, b)
    end
    return cost
end

"""
    address_cost(sa_struct{T}, node :: SANode, address_data::AddressData) where {T <: RuleSet}

Return the address cost for `node` for RuleSet `T`. Default return value
is `zero(Float64)`.

Called by default during [`node_cost`](@ref) and [`node_pair_cost`](@ref)

Method List
-----------
$(METHODLIST)
"""
address_cost(sa_struct::SAStruct, node::SANode, address_data::AddressData) = zero(Float64)
function address_cost(sa_struct::SAStruct, node::SANode, address_data::DefaultAddressData)
    return address_data[getlocation(node)]
end

# Convenience wrapper.
function address_cost(sa_struct::SAStruct, node::SANode)
    return address_cost(sa_struct, node, sa_struct.address_data)
end

"""
    aux_cost(sa_struct{T}) where {T <: RuleSet}

Return an auxiliary cost associated with the entire mapping of the `sa_struct`.
May use any field of `sa_struct` but may only mutate `sa_struct.aux`.

Default: `zero(Float64)`

Method List
-----------
$(METHODLIST)
"""
aux_cost(sa_struct::SAStruct) = zero(Float64)

function map_cost(sa_struct::SAStruct)
    cost = aux_cost(sa_struct)
    for i in eachindex(sa_struct.channels)
        cost += channel_cost(sa_struct, i)
    end
    for n in sa_struct.nodes
        cost += address_cost(sa_struct, n)
    end
    return cost
end

"""
    node_cost(sa_struct{T}, idx) where {T <: RuleSet}

Return the cost of the node with index `idx` in ruleset `T`.

Default implementation sums the node's:
* incoming channels
* outgoing channels
* address cost
* auxiliary cost.

Method List
-----------
$(METHODLIST)
"""
@propagate_inbounds function node_cost(sa_struct::SAStruct, idx::Integer)
    # Unpack node data type.
    n = sa_struct.nodes[idx]
    cost = aux_cost(sa_struct)
    for channel in n.outchannels
        cost += channel_cost(sa_struct, channel)
    end
    for channel in n.inchannels
        cost += channel_cost(sa_struct, channel)
    end
    cost += address_cost(sa_struct, n)

    return cost
end

#=
Node pair cost is slightly more subtle than just each independent node cost if

    1. The communication resources in the array are asymmetric. Then if two
        nodes are swapped and are connected by a channel, the double cost of
        that channel is not cancelled out correctly.

    2. For multi-pin nets, if a source and sink are swapped, then the the
        objective function is calculated incorrectly due to counting the channel
        twice.
=#
"""
    node_pair_cost(sa_struct{T}, idx1, idx2) where {T <: RuleSet}

Compute the cost of the pair of nodes with indices `idx1` and `idx2`. Call this
function when computing the cost of two nodes because in general, the total
cost of two nodes is not the sum of the individual nodes' costs.

Method List
-----------
$(METHODLIST)
"""
@propagate_inbounds function node_pair_cost(sa_struct::SAStruct, idx1, idx2)
    cost = node_cost(sa_struct, idx1)
    # Get the two node types for calculating the cost of the second node.
    a = sa_struct.nodes[idx1]
    b = sa_struct.nodes[idx2]

    for channel in b.outchannels
        if !in(channel, a.inchannels)
            cost += channel_cost(sa_struct, channel)
        end
    end
    for channel in b.inchannels
        if !in(channel, a.outchannels)
            cost += channel_cost(sa_struct, channel)
        end
    end
    cost += address_cost(sa_struct, b)
    return cost
end


text:
# we assume that kspan is created by a command like
#    [(i,j) for i=range(kxmin,kxmax,lengtH=Nkx), j=range(kymin,kymax,length=Nky)]
# i.e. that kspan is supplied as an Nkx×Nky array of 2-element tuples
function sdos(eigvals::AbstractVector{<:Real},
              eigvecs::AbstractMatrix{<:Number},
              tb::TBModel,
              kspan::AbstractArray = (ks=range(-3π,3π,length=50); [(kx,ky) for kx in ks, ky in ks]),
              enespan::AbstractVector{<:Real} = range(extrema(eigvals)...,length=100),
              eta::Real = 0.1)

    eigvals_c = eigvals .- eta*im # add finite loss

    eigvec_k_abs2 = zeros(Float64,    size(eigvecs, 1))
    expfac        = zeros(ComplexF64, size(tb.pos, 1))
    term          = zeros(Float64,    size(eigvals))
    sdosvals      = zeros(Float64,    (size(kspan)..., size(enespan)...))

    @showprogress 2 "Calculating sdos ... " for kidx in CartesianIndices(kspan)
        expfac .= cis.( kspan[kidx][1].*tb.pos[:,1] .+ kspan[kidx][2].*tb.pos[:,2] )

        eigvec_k_abs2 .= abs2.(dropdims(sum(repeat(expfac, tb.basisdim) .* eigvecs, dims=1),
                                        dims=1)) # TODO: ... probably allocates...

        for eneidx = eachindex(enespan)
            term .= imag.( eigvec_k_abs2 ./ (enespan[eneidx] .- eigvals_c) )
            sdosvals[kidx,eneidx] = -sum(term)
        end
    end

    return sdosvals
end

ext: jl

max_issues_repo_name: thchr/super-tb

max_forks_repo_name: thchr/super-tb


text:
module Math

export sin, cos, tan, sinh, cosh, tanh, asin, acos, atan,
       asinh, acosh, atanh, sec, csc, cot, asec, acsc, acot,
       sech, csch, coth, asech, acsch, acoth,
       sinpi, cospi, sinc, cosc,
       cosd, cotd, cscd, secd, sind, tand,
       acosd, acotd, acscd, asecd, asind, atand, atan2,
       rad2deg, deg2rad,
       log, log2, log10, log1p, exponent, exp, exp2, exp10, expm1,
       cbrt, sqrt, erf, erfc, erfcx, erfi, dawson,
       ceil, floor, trunc, round, significand,
       lgamma, hypot, gamma, lfact, max, min, minmax, ldexp, frexp,
       clamp, modf, ^, mod2pi,
       airy, airyai, airyprime, airyaiprime, airybi, airybiprime, airyx,
       besselj0, besselj1, besselj, besseljx,
       bessely0, bessely1, bessely, besselyx,
       hankelh1, hankelh2, hankelh1x, hankelh2x,
       besseli, besselix, besselk, besselkx, besselh,
       beta, lbeta, eta, zeta, polygamma, invdigamma, digamma, trigamma,
       erfinv, erfcinv, @evalpoly

import Base: log, exp, sin, cos, tan, sinh, cosh, tanh, asin,
             acos, atan, asinh, acosh, atanh, sqrt, log2, log10,
             max, min, minmax, ceil, floor, trunc, round, ^, exp2,
             exp10, expm1, log1p

import Core.Intrinsics: nan_dom_err, sqrt_llvm, box, unbox, powi_llvm

# non-type specific math functions

clamp{X,L,H}(x::X, lo::L, hi::H) =
    ifelse(x > hi, convert(promote_type(X,L,H), hi),
           ifelse(x < lo,
                  convert(promote_type(X,L,H), lo),
                  convert(promote_type(X,L,H), x)))

clamp{T}(x::AbstractArray{T,1}, lo, hi) = [clamp(xx, lo, hi) for xx in x]
clamp{T}(x::AbstractArray{T,2}, lo, hi) =
    [clamp(x[i,j], lo, hi) for i in 1:size(x,1), j in 1:size(x,2)]
clamp{T}(x::AbstractArray{T}, lo, hi) =
    reshape([clamp(xx, lo, hi) for xx in x], size(x))

# evaluate p[1] + x * (p[2] + x * (....)), i.e. a polynomial via Horner's rule
macro horner(x, p...)
    ex = esc(p[end])
    for i = length(p)-1:-1:1
        ex = :($(esc(p[i])) + t * $ex)
    end
    Expr(:block, :(t = $(esc(x))), ex)
end

# Evaluate p[1] + z*p[2] + z^2*p[3] + ... + z^(n-1)*p[n].  This uses
# Horner's method if z is real, but for complex z it uses a more
# efficient algorithm described in Knuth, TAOCP vol. 2, section 4.6.4,
# equation (3).
macro evalpoly(z, p...)
    a = :($(esc(p[end])))
    b = :($(esc(p[end-1])))
    as = {}
    for i = length(p)-2:-1:1
        ai = symbol(string("a", i))
        push!(as, :($ai = $a))
        a = :($b + r*$ai)
        b = :($(esc(p[i])) - s * $ai)
    end
    ai = :a0
    push!(as, :($ai = $a))
    C = Expr(:block,
             :(t = $(esc(z))),
             :(x = real(t)),
             :(y = imag(t)),
             :(r = x + x),
             :(s = x*x + y*y),
             as...,
             :($ai * t + $b))
    R = Expr(:macrocall, symbol("@horner"), esc(z), p...)
    :(isa($(esc(z)), Complex) ? $C : $R)
end

rad2deg(z::Real) = oftype(z, 57.29577951308232*z)
deg2rad(z::Real) = oftype(z, 0.017453292519943295*z)
rad2deg(z::Integer) = rad2deg(float(z))
deg2rad(z::Integer) = deg2rad(float(z))
@vectorize_1arg Real rad2deg
@vectorize_1arg Real deg2rad

log(b,x) = log(x)./log(b)

# type specific math functions

const libm = Base.libm_name
const openspecfun = "libopenspecfun"

# functions with no domain error
for f in (:cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp, :erf, :erfc, :exp2, :expm1)
    @eval begin
        ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x)
        ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x)
        ($f)(x::Real) = ($f)(float(x))
        @vectorize_1arg Number $f
    end
end

# fallback definitions to prevent infinite loop from $f(x::Real) def above
cbrt(x::FloatingPoint) = x^(1//3)
exp2(x::FloatingPoint) = 2^x
for f in (:sinh, :cosh, :tanh, :atan, :asinh, :exp, :erf, :erfc, :expm1)
    @eval ($f)(x::FloatingPoint) = error("not implemented for ", typeof(x))
end

# TODO: GNU libc has exp10 as an extension; should openlibm?
exp10(x::Float64) = 10.0^x
exp10(x::Float32) = 10.0f0^x
exp10(x::Integer) = exp10(float(x))
@vectorize_1arg Number exp10

# functions that return NaN on non-NaN argument for domain error
for f in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
          :lgamma, :log1p)
    @eval begin
        ($f)(x::Float64) = nan_dom_err(ccall(($(string(f)),libm), Float64, (Float64,), x), x)
        ($f)(x::Float32) = nan_dom_err(ccall(($(string(f,"f")),libm), Float32, (Float32,), x), x)
        ($f)(x::Real) = ($f)(float(x))
        @vectorize_1arg Number $f
    end
end

sqrt(x::Float64) = box(Float64,sqrt_llvm(unbox(Float64,x)))
sqrt(x::Float32) = box(Float32,sqrt_llvm(unbox(Float32,x)))
sqrt(x::Real) = sqrt(float(x))
@vectorize_1arg Number sqrt

for f in (:ceil, :trunc, :significand) # :rint, :nearbyint
    @eval begin
        ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x)
        ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x)
        @vectorize_1arg Real $f
    end
end

round(x::Float32) = ccall((:roundf, libm), Float32, (Float32,), x)
@vectorize_1arg Real round

floor(x::Float32) = ccall((:floorf, libm), Float32, (Float32,), x)
@vectorize_1arg Real floor

hypot(x::Real, y::Real) = hypot(promote(float(x), float(y))...)
function hypot{T<:FloatingPoint}(x::T, y::T)
    x = abs(x)
    y = abs(y)
    if x < y
        x, y = y, x
    end
    if x == 0
        r = y/one(x)
    else
        r = y/x
        if isnan(r)
            isinf(x) && return x
            isinf(y) && return y
            return r
        end
    end
    x * sqrt(one(r)+r*r)
end

atan2(x::Real, y::Real) = atan2(promote(float(x),float(y))...)
atan2{T<:FloatingPoint}(x::T, y::T) = Base.no_op_err("atan2", T)

for f in (:atan2, :hypot)
    @eval begin
        ($f)(x::Float64, y::Float64) = ccall(($(string(f)),libm), Float64, (Float64, Float64,), x, y)
        ($f)(x::Float32, y::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32, Float32), x, y)
        @vectorize_2arg Number $f
    end
end

max{T<:FloatingPoint}(x::T, y::T) = ifelse((y > x) | (x != x), y, x)
@vectorize_2arg Real max

min{T<:FloatingPoint}(x::T, y::T) = ifelse((y < x) | (x != x), y, x)
@vectorize_2arg Real min

minmax{T<:FloatingPoint}(x::T, y::T) =  x <= y ? (x, y) :
                                        x >  y ? (y, x) :
                                        x == x ? (x, x) : (y, y)

function exponent(x::Float64)
    if x==0 || !isfinite(x)
        throw(DomainError())
    end
    int(ccall((:ilogb,libm), Int32, (Float64,), x))
end
function exponent(x::Float32)
    if x==0 || !isfinite(x)
        throw(DomainError())
    end
    int(ccall((:ilogbf,libm), Int32, (Float32,), x))
end
@vectorize_1arg Real exponent

ldexp(x::Float64,e::Int) = ccall((:scalbn,libm),  Float64, (Float64,Int32), x, int32(e))
ldexp(x::Float32,e::Int) = ccall((:scalbnf,libm), Float32, (Float32,Int32), x, int32(e))
# TODO: vectorize ldexp

begin
    local exp::Array{Int32,1} = zeros(Int32,1)
    global frexp
    function frexp(x::Float64)
        s = ccall((:frexp,libm), Float64, (Float64, Ptr{Int32}), x, exp)
        (s, int(exp[1]))
    end
    function frexp(x::Float32)
        s = ccall((:frexpf,libm), Float32, (Float32, Ptr{Int32}), x, exp)
        (s, int(exp[1]))
    end
    function frexp(A::Array{Float64})
        f = similar(A)
        e = Array(Int, size(A))
        for i = 1:length(A)
            f[i] = ccall((:frexp,libm), Float64, (Float64, Ptr{Int32}), A[i], exp)
            e[i] = exp[1]
        end
        return (f, e)
    end
    function frexp(A::Array{Float32})
        f = similar(A)
        e = Array(Int, size(A))
        for i = 1:length(A)
            f[i] = ccall((:frexpf,libm), Float32, (Float32, Ptr{Int32}), A[i], exp)
            e[i] = exp[1]
        end
        return (f, e)
    end
end

modf(x) = rem(x,one(x)), trunc(x)

const _modff_temp = Float32[0]
function modf(x::Float32)
    f = ccall((:modff,libm), Float32, (Float32,Ptr{Float32}), x, _modff_temp)
    f, _modff_temp[1]
end

const _modf_temp = Float64[0]
function modf(x::Float64)
    f = ccall((:modf,libm), Float64, (Float64,Ptr{Float64}), x, _modf_temp)
    f, _modf_temp[1]
end

^(x::Float64, y::Float64) = nan_dom_err(ccall((:pow,libm),  Float64, (Float64,Float64), x, y), x+y)
^(x::Float32, y::Float32) = nan_dom_err(ccall((:powf,libm), Float32, (Float32,Float32), x, y), x+y)

^(x::Float64, y::Integer) =
    box(Float64, powi_llvm(unbox(Float64,x), unbox(Int32,int32(y))))
^(x::Float32, y::Integer) =
    box(Float32, powi_llvm(unbox(Float32,x), unbox(Int32,int32(y))))

function angle_restrict_symm(theta)
    const P1 = 4 * 7.8539812564849853515625e-01
    const P2 = 4 * 3.7748947079307981766760e-08
    const P3 = 4 * 2.6951514290790594840552e-15

    y = 2*floor(theta/(2*pi))
    r = ((theta - y*P1) - y*P2) - y*P3
    if (r > pi)
        r -= (2*pi)
    end
    return r
end

## mod2pi-related calculations ##

function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64)
    # as above, but only compute and return high double
    r = xh+yh
    s = (abs(xh) > abs(yh)) ? (xh-r+yh+yl+xl) : (yh-r+xh+xl+yl)
    zh = r+s
    return zh
end

function ieee754_rem_pio2(x::Float64)
    # rem_pio2 essentially computes x mod pi/2 (ie within a quarter circle)
    # and returns the result as
    # y between + and - pi/4 (for maximal accuracy (as the sign bit is exploited)), and
    # n, where n specifies the integer part of the division, or, at any rate,
    # in which quadrant we are.
    # The invariant fulfilled by the returned values seems to be
    #  x = y + n*pi/2 (where y = y1+y2 is a double-double and y2 is the "tail" of y).
    # Note: for very large x (thus n), the invariant might hold only modulo 2pi
    # (in other words, n might be off by a multiple of 4, or a multiple of 100)

    # this is just wrapping up
    # https://github.com/JuliaLang/openlibm/blob/master/src/e_rem_pio2.c

    y = [0.0,0.0]
    n = ccall((:__ieee754_rem_pio2, libm), Cint, (Float64,Ptr{Float64}), x, y)
    return (n,y)
end

# multiples of pi/2, as double-double (ie with "tail")
const pi1o2_h  = 1.5707963267948966     # convert(Float64, pi * BigFloat(1/2))
const pi1o2_l  = 6.123233995736766e-17  # convert(Float64, pi * BigFloat(1/2) - pi1o2_h)

const pi2o2_h  = 3.141592653589793      # convert(Float64, pi * BigFloat(1))
const pi2o2_l  = 1.2246467991473532e-16 # convert(Float64, pi * BigFloat(1) - pi2o2_h)

const pi3o2_h  = 4.71238898038469       # convert(Float64, pi * BigFloat(3/2))
const pi3o2_l  = 1.8369701987210297e-16 # convert(Float64, pi * BigFloat(3/2) - pi3o2_h)

const pi4o2_h  = 6.283185307179586      # convert(Float64, pi * BigFloat(2))
const pi4o2_l  = 2.4492935982947064e-16 # convert(Float64, pi * BigFloat(2) - pi4o2_h)

function mod2pi(x::Float64) # or modtau(x)
# with r = mod2pi(x)
# a) 0 <= r < 2π  (note: boundary open or closed - a bit fuzzy, due to rem_pio2 implementation)
# b) r-x = k*2π with k integer

# note: mod(n,4) is 0,1,2,3; while mod(n-1,4)+1 is 1,2,3,4.
# We use the latter to push negative y in quadrant 0 into the positive (one revolution, + 4*pi/2)

    if x < pi4o2_h
        if 0.0 <= x return x end
        if x > -pi4o2_h
            return add22condh(x,0.0,pi4o2_h,pi4o2_l)
        end
    end

    (n,y) = ieee754_rem_pio2(x)

    if iseven(n)
        if n & 2 == 2 # add pi
            return add22condh(y[1],y[2],pi2o2_h,pi2o2_l)
        else # add 0 or 2pi
            if y[1] > 0.0
                return y[1]
            else # else add 2pi
                return add22condh(y[1],y[2],pi4o2_h,pi4o2_l)
            end
        end
    else # add pi/2 or 3pi/2
        if n & 2 == 2 # add 3pi/2
            return add22condh(y[1],y[2],pi3o2_h,pi3o2_l)
        else # add pi/2
            return add22condh(y[1],y[2],pi1o2_h,pi1o2_l)
        end
    end
end

mod2pi(x::Float32) = float32(mod2pi(float64(x)))
mod2pi(x::Int32) = mod2pi(float64(x))
function mod2pi(x::Int64)
  fx = float64(x)
  fx == x || error("Integer argument to mod2pi is too large: $x")
  mod2pi(fx)
end

# More special functions

include("special/trig.jl")
include("special/bessel.jl")
include("special/erf.jl")
include("special/gamma.jl")

end # module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Plots, StatsBase, Statistics, LaTeXStrings, JLD

function InitialNetwork(dim, P)
    Network = sample([-1,0], Weights([1-P, P]),(dim,dim))
    return Network
end

function FindNeighbors(i, j, Network, dim)
    Numbereds = 0
    colorList = []
    for neighbor in ([0,-1], [-1,0])
        if j + neighbor[2] >= 1 && i + neighbor[1] >= 1 && Network[([i, j] + neighbor)...] > 0
            push!(colorList,Network[([i, j] + neighbor)...])
            Numbereds += 1
        end
    end
    return colorList, Numbereds #colorList contain the node and its neighbors
end

function RecursionMapping(x, L)
    while L[x] != L[L[x]]
        x = L[x]
    end
    return L[x]
end

@inline function CheckPercolation(dim, L, Network)
    fside = []
    lside = []
    for i in 1:dim
        if Network[i,dim] != -1
            push!(lside, RecursionMapping(Network[i,dim],L))
        end
        if Network[i,1] != -1
            push!(fside, RecursionMapping(Network[i,1],L))
        end
    end
    if length(intersect(fside,lside)) > 0
        return 1
    else
        return 0
    end
end

@inline function HKNetworkDynamic(dim, P)
    Network = InitialNetwork(dim, P)
    S = []
    L = []
    k = 1

    for j in 1:dim
        for i in 1:dim
            if Network[i,j] != 0
                continue
            end
            NColor, NNum = FindNeighbors(i,j,Network,dim)
            if NNum == 0
                Network[i,j] = k
                push!(L, k)
                push!(S, 1)
                k += 1
            elseif NNum == 1
                Network[i,j] = RecursionMapping(NColor[1],L)
                S[RecursionMapping(NColor[1],L)] += 1
            elseif NNum == 2 && RecursionMapping(NColor[1],L) == RecursionMapping(NColor[2],L)
                Network[i,j] = RecursionMapping(NColor[1],L)
                S[RecursionMapping(NColor[1],L)] += 1
            else
                S[RecursionMapping(NColor[1],L)] += 1 + S[RecursionMapping(NColor[2],L)]
                S[RecursionMapping(NColor[2],L)] = 0
                Network[i,j] = RecursionMapping(NColor[1],L)
                L[NColor[2]] = RecursionMapping(NColor[1],L)
            end
        end
    end
    # return Network
    return CheckPercolation(dim, L, Network)
end

@time HKNetworkDynamic(100, 0.6)

Plist = hcat(0.55:0.002:0.65)
Avglist = []
STDlist = []
TotData = []
# PercolationCheck(40,0.55,1)
for P in Plist
    print("\r$P")
    AllRuns = []
    for i in 1:1000
        push!(AllRuns, HKNetworkDynamic(200, P))
    end
    push!(TotData, AllRuns)
    push!(STDlist, std(AllRuns))
    push!(Avglist, mean(AllRuns))
end

save("200x200-HK.jld", "data", TotData)

scatter(Plist,
    Avglist,
    yerr = STDlist,
    legend = nothing,
    xlabel = L"P",
    ylabel = L"Avg\ Q",
    title = L"Percolation \ for \ 200\times200 \ Network\ (1000 \ runs)")
savefig("../../Figs/Q5/200x200.pdf")


text:
###################################################################
#                TEST Metric estimations
###################################################################
using Test
#include("../src/OptimizedNeuralNetwork.jl")
include("../src/metric_estimation.jl")

@testset "Regression metrics" begin
    f(x)=x.*x
    X = [1,-1,-2,2]
    Y = map(f, X)
    # Comprobamos que devuelve que para este caso en concreto son correctas:
    # la media de error, mediana de error, la desviación media y el coeficiente de correlación
    @test regression(X,Y,f) == (0,0,0,1)
end


text:
export AbstractPolicy,
       update!

"""
    AbstractPolicy

Supertype for policies. A policy takes in a state and returns an action.
"""
abstract type AbstractPolicy end

function update! end


text:
using DashBootstrapComponents

text_input = dbc_formgroup([
    dbc_label("Text"),
    dbc_input(placeholder="Input goes here...", type="text"),
    dbc_formtext("Type something in the box above"),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
module TestUncertainties

using Modia
@usingModiaPlot
using Modia.Measurements

TestModel = Model(
    T = 0.2 ± 0.02, # Cannot convert a particle distribution to a float if not all particles are the same.
    u0 = 10 ± 1, # Cannot convert a particle distribution to a float if not all particles are the same.
    k = 1 ± 0.1,
    x = Var(init = 5 ± 0.5),
    v = Var(init = 1 ± 0.1),

    m = 10 ± 1,
    F = 20 ± 2,
    equations = :[
    T*der(x) = u0
    y = k*x

    F = m * a
    der(v) = a # Cannot convert a particle distribution to a float if not all particles are the same.
    ]
)

model = @instantiateModel(TestModel, FloatType = Measurement{Float64})

simulate!(model, stopTime = 1.0, log=false,
          requiredFinalStates = Measurements.Measurement{Float64}[54.99999999999998 ± 7.08872343937891, 2.999999999999999 ± 0.3])

plot(model, ["T", "x", "der(x)", "y", "a", "der(v)"], figure=1)

include("../models/Electric.jl")

FilterCircuit = Model(
    R = Resistor        | Map(R = (100±10)u"Ω"),
    C = Capacitor       | Map(C = (0.01±0.001)u"F", v=Var(init=(0.0±1.0)u"V")),
    V = ConstantVoltage | Map(V = (10.0±1.0)u"V"),
    ground = Ground,
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

filterCircuit = @instantiateModel(FilterCircuit, FloatType=Measurements.Measurement{Float64})
simulate!(filterCircuit, Tsit5(), stopTime = 3.0, logParameters = true, logStates = true)
plot(filterCircuit, ("V.v", "C.v"), figure=2)

Pendulum = Model(
   L = (0.8±0.1)u"m",
   m = (1.0±0.1)u"kg",
   d = (0.5±0.05)u"N*m*s/rad",
   g = 9.81u"m/s^2",
   phi = Var(init = (pi/2±0.1)*u"rad"),
   w   = Var(init = 0u"rad/s"),
   equations = :[
          w = der(phi)
        0.0 = m*L^2*der(w) + d*w + m*g*L*sin(phi)
   ]
)

pendulum = @instantiateModel(Pendulum, FloatType=Measurements.Measurement{Float64})
simulate!(pendulum, Tsit5(), stopTime = 10.0)
plot(pendulum, "phi", figure=4)

end


text:
module PackageCompiler

using Libdl, SnoopCompile
Sys.iswindows() && using WinRPM

include("compiler_flags.jl")
include("static_julia.jl")
include("api.jl")
include("snooping.jl")
include("system_image.jl")

const sysimage_binaries = ("sys.$(Libdl.dlext)",)

function copy_system_image(src, dest, ignore_missing = false)
    for file in sysimage_binaries
        # backup
        srcfile = joinpath(src, file)
        destfile = joinpath(dest, file)
        if !isfile(srcfile)
            ignore_missing && continue
            error("No file: $srcfile")
        end
        if isfile(destfile)
            if isfile(destfile * ".backup")
                rm(destfile * ".backup", force = true)
            end
            mv(destfile, destfile * ".backup", force = true)
        end
        @info "Copying system image: $srcfile to $destfile"
        cp(srcfile, destfile, force = true)
    end
end

julia_cpu_target(x) = error("CPU target needs to be a string or `nothing`")
julia_cpu_target(x::String) = x # TODO: match against available targets
function julia_cpu_target(::Nothing)
    replace(Base.julia_cmd().exec[2], "-C" => "")
end

"""
Reverts a forced compilation of the system image.
This will restore any previously backed up system image files, or
build a new, clean system image.
"""
function revert(debug = false)
    syspath = default_sysimg_path(debug)
    sysimg_backup = dirname(get_backup!(debug))
    copy_system_image(sysimg_backup, syspath)
end

function get_root_dir(path)
    path, name = splitdir(path)
    if isempty(name)
        return splitdir(path)[2]
    else
        name
    end
end

function sysimg_folder(files...)
    base_path = normpath(abspath(joinpath(@__DIR__, "..", "sysimg")))
    isdir(base_path) || mkpath(base_path)
    normpath(abspath(joinpath(base_path, files...)))
end

function sysimgbackup_folder(files...)
    backup = sysimg_folder("backup")
    isdir(backup) || mkpath(backup)
    sysimg_folder("backup", files...)
end

function package_folder(package...)
    packages = normpath(abspath(joinpath(@__DIR__, "..", "packages")))
    isdir(packages) || mkpath(packages)
    normpath(abspath(joinpath(packages, package...)))
end

"""
    compile_package(packages...; kw_args...)

with packages being either a string naming a package, or a tuple `(package_name, precompile_file)`.
If no precompile file is given, it will use the packages `runtests.jl`, which is a good canditate
for figuring out what functions to compile!
"""
function compile_package(packages...; kw_args...)
    args = map(packages) do package
        # If no explicit path to a seperate precompile file, use runtests
        isa(package, String) && return (package, "test/runtests.jl")
        isa(package, Tuple{String, String}) && return package
        error("Unrecognized package. Use `packagename::String`, or `(packagename::String, rel_path_to_testfile::String)`. Found: `$package`")
    end
    compile_package(args...; kw_args...)
end

"""
    compile_package(packages::Tuple{String, String}...; force = false, reuse = false, debug = false, cpu_target = nothing)

Compile a list of packages. Each package comes as a tuple of `(package_name, precompile_file)`
where the precompile file should contain all function calls, that should get compiled into the system image.
Usually the `runtests.jl` file is a good candidate, since it should run all important functions of a package.
"""
function compile_package(packages::Tuple{String, String}...; force = false, reuse = false, debug = false, cpu_target = nothing)
    userimg = sysimg_folder("precompile.jl")
    if !reuse
        snoop_userimg(userimg, packages...)
    end
    !isfile(userimg) && reuse && error("Nothing to reuse. Please run `compile_package(reuse = true)`")
    image_path = sysimg_folder()
    build_sysimg(image_path, userimg, cpu_target=cpu_target)
    imgfile = joinpath(image_path, "sys.$(Libdl.dlext)")
    syspath = joinpath(default_sysimg_path(debug), "sys.$(Libdl.dlext)")
    if force
        try
            backup = syspath * ".packagecompiler_backup"
            isfile(backup) || mv(syspath, backup)
            cp(imgfile, syspath)
            @info """Replaced system image successfully. Next start of julia will load the newly compiled system image.
                If you encounter any errors with the new julia image, try `PackageCompiler.revert([debug = false])`."""
        catch e
            @warn "An error occured while replacing sysimg files:" error=e
            @info "Recovering old system image from backup"
            # if any file is missing in default system image, revert!
            if !isfile(syspath)
                @info "$syspath missing. Reverting!"
                revert(debug)
            end
        end
    else
        @info """Not replacing system image.
            You can start julia with $(`julia -J $imgfile`) at a posix shell to load the compiled files."""
    end
    imgfile
end

function __init__()
    if Base.julia_cmd().exec[2] != "-Cnative"
        @warn "Your Julia system image is not compiled natively for this CPU architecture.\n" *
            "Please run `PackageCompiler.force_native_image!()` for optimal Julia performance."
    end
end

export compile_package, revert, force_native_image!, executable_ext, build_executable, build_shared_lib, static_julia

end # module


text:
# Real Space

export RealSpace

const RealSpaceType{T} = Union{T, Array{T, N} where N}

"""
    RealSpace{T, DType, S} <: AbstractSpace{T, S}

A real space contains element of type `DType` (e.g Array or a Number), the
numeric content of each element has type `T` (`eltype(DType) = T`)

A naive example is the 1-D real space

```julia-repl
julia> RealSpace(min=-1, max=1)
QMTK.RealSpace{Float64,Float64,QMTK.Randomized}(-1.0, 1.0, -0.04293376352914002)
```
"""
mutable struct RealSpace{T<:Real, DType<:RealSpaceType{T}, S} <: AbstractSpace{T, S}
    min::T
    max::T
    data::DType
end

# Utils
_uniform(min, max) = _uniform(Float64, min, max)
_uniform(shape::Tuple, min, max) = _uniform(Float64, shape, min, max)
_uniform(::Type{T}, min, max) where T = T(max - min) * rand(T) + T(min)
_uniform(::Type{T}, shape::Tuple, min, max) where T = T(max-min) * rand(T, shape) + T(min)

# Constructors
RealSpace(;min=0, max=1) = RealSpace(Float64, min=min, max=max)
RealSpace(shape::Integer...; min=0, max=1) = RealSpace(Float64, shape, min=min, max=max)
RealSpace(::Type{T}, shape::Integer...; min=0, max=1) where T = RealSpace(T, shape, min=min, max=max)

# Basic Constructors
RealSpace(::Type{T}; min=0, max=1) where T = RealSpace{T, T, Randomized}(min, max, _uniform(T, min, max))
RealSpace(::Type{T}, shape::Tuple; min=0, max=1) where T =
    RealSpace{T, Array{T, length(shape)}, Randomized}(min, max, _uniform(T, shape, min, max))

RealSpace(::Type{S}, space::RealSpace{T, D}) where {T, D, S} = RealSpace{T, D, S}(space.min, space.max, space.data)

function reset!(space::RealSpace{T, D, Randomized}) where {T, D <: AbstractArray}
    fill!(space.data, space.min)
    return RealSpace(UnRandomized, space)
end

function reset!(space::RealSpace{T, D, Randomized}) where {T, D <: Real}
    space.data = space.min
    return RealSpace(UnRandomized, space)
end

copy(space::RealSpace{T, D, S}) where {T, D <: Real, S} =
    RealSpace{T, D, S}(space.min, space.max, space.data)
copy(space::RealSpace{T, D, S}) where {T, D <: AbstractArray, S} =
    RealSpace{T, D, S}(space.min, space.max, copy(space.data))

acquire(space::RealSpace{T, D, S}) where {T, D <: Real, S} = space.data
acquire(space::RealSpace{T, D, S}) where {T, D <: AbstractArray, S} = copy(space.data)

traverse(space::RealSpace; step::Real=1e-2) = RealTraverser(space, step)

import Base: start, next, done, length, eltype

struct RealTraverser{T, D}
    data::RealSpace{T, D}
    step::T
end

# TODO: traverse array field
length(itr::RealTraverser) = floor(Int, (itr.data.max - itr.data.min) / itr.step) + 1
eltype(itr::RealTraverser{T, D}) where {T, D <: Real} = D

start(itr::RealTraverser{T, D}) where {T, D <: Real} = 0
next(itr::RealTraverser{T, D}, state) where {T, D <: Real} = (itr.data.min + state * itr.step, state + 1)
done(itr::RealTraverser{T, D}, state) where {T, D <: Real} = itr.data.min + state * itr.step > itr.data.max

# Properties
data(space::RealSpace) = space.data

# Randomize methods
# Note: shake! will be dispatched, we use a temperory function

@inline function _shake!(space::RealSpace{T, DType}) where {T, DType <: Real}
    space.data = (space.max - space.min) * rand() + space.min
end

@inline function _shake!(space::RealSpace{T, DType}) where {T, DType <: AbstractArray}
    rand!(space.data)
    space.data.*= (space.max - space.min)
    space.data.+= space.min
    return space.data
end

shake!(space::RealSpace{T, D, Randomized}) where {T, D} = _shake!(space)

function randomize!(space::RealSpace{T, D, UnRandomized}) where {T, D}
    _shake!(space)
    return RealSpace(Randomized, space)
end


text:
using SimBioReactor
using Test
using Gtk

@testset "mainGUI" begin
    w =SimBioReactorGUI()
    destroy(w)
end


text:
###################################################################################################
# CONSTANTS
###################################################################################################
const ETA_MAX_ABS=5.0
const CI_THRESH=1.0
###################################################################################################
# FUNCTIONS
###################################################################################################
"""
    `create_Ux([N1,...,NK],[α1,...,αK],β)`

Function that creates a function to compute potential energy for a region with [N1,...,NK],
parameters [α1,...,αK], and correlation β.

# Examples
```julia-repl
julia> Ux_fun = CpelAsm.create_Ux([2,2],[1.0,1.0],1.0);
julia> Ux_fun([1,1,1,1])
-7.0
```
"""
function create_Ux(n::Vector{Int64},a::Vector{Float64},b::Float64)

    # Define function given n, α, and β
    function Ux(x::Vector{Int64})::Float64
        u = a[1] * sum(x[1:n[1]]) + b * sum(x[1:(end-1)] .* x[2:end])
        @inbounds for i in 2:length(n)
            u += a[i] * sum(x[(sum(n[1:(i-1)])+1):(sum(n[1:i]))])
        end
        return -u
    end

    # Return function
    return Ux

end # end create_Ux
"""
    `get_W(N,α,β)`

Function that returns W^{N-1} computed efficiently via sum of rank-1 matrices
assuming α and β.

# Examples
```julia-repl
julia> CpelAsm.get_W(4,0.0,0.0)
2×2 Array{Float64,2}:
 4.0  4.0
 4.0  4.0
```
"""
function get_W(n::Int64, a::Float64, b::Float64)::Array{Float64,2}

    # if n=1 then return identity
    n==1 && return [1.0 0.0;0.0 1.0]

    # Compute relevant quantities
    exp_a = exp(a)
    exp_b = exp(b)
    cosh_a = 0.5*(exp_a+1.0/exp_a)
    sinh_a = exp_a-cosh_a

    # Eigenvalues
    aux1 = exp_b * cosh_a
    aux2 = sqrt(1.0 + exp_b^4*sinh_a^2)
    lambda1N = (aux1-1.0/exp_b*aux2)^(n-1)
    lambda2N = (aux1+1.0/exp_b*aux2)^(n-1)

    # Eigenvectors
    aux1 = -exp_b^2 * sinh_a
    e1 = [aux1-aux2; 1.0]
    e1 /= sqrt(e1'*e1)
    e2 = [aux1+aux2; 1.0]
    e2 /= sqrt(e2'*e2)

    # Return W^{N-1}
    return e1*e1'*lambda1N + e2*e2'*lambda2N

end
"""
    `get_V([α1,α2],β)`

Function that returns V matrix assuming [α1,α2] and β.

# Examples
```julia-repl
julia> CpelAsm.get_V([1.0,1.0],1.0)
2×2 Array{Float64,2}:
 1.0       0.367879
 0.367879  7.38906
```
"""
function get_V(a::Vector{Float64},b::Float64)::Array{Float64,2}

    # Compute auxiliary vars
    exp_b = exp(b)
    exp_a_p = exp(0.5*sum(a))
    exp_a_m = exp(0.5*(a[2]-a[1]))

    # Return V
    return [exp_b/exp_a_p exp_a_m/exp_b;1.0/(exp_b*exp_a_m) exp_b*exp_a_p]

end
"""
    `get_u(α)`

Function that boundary vector for Z computation.

# Examples
```julia-repl
julia> CpelAsm.get_u(0.0)
2-element Vector{Float64}:
 1.0
 1.0
```
"""
function get_u(a::Float64)::Vector{Float64}

    # Compute auxiliary
    exp_aux = exp(a/2.0)

    # Return u
    return [1.0/exp_aux; exp_aux]

end
"""
    `comp_Z([N1,...,NK],[α1,...,αK],β)`

Compute partition function of a model with [N1,...,NK] CpG cites and with
parameters [α1,...,αK] and β.

# Examples
```julia-repl
julia> CpelAsm.comp_Z([1,1,1],[1.0,1.0,1.0],1.0)
155.37102759254836
```
"""
function comp_Z(n::Vector{Int64},a::Vector{Float64},b::Float64)::Float64

    # Boundary conditions.
    y = get_u(a[1])'*get_W(n[1],a[1],b)
    if length(n)>1
        y *= prod([get_V(a[(i-1):i],b)*get_W(n[i],a[i],b) for i in 2:length(n)])
    end
    y *= get_u(a[end])

    # Return Z
    return max(1e-100,y)

end
"""
    `get_grad_logZ([N1,...,NK],θhat)`

Numerically computes the gradient of the log partition function of a model with [N1,...,NK]
CpG cites and estimated parameter vector θhat. This is equivalent to computing the expected
value of sufficient statistis (SS).

# Examples
```julia-repl
julia> CpelAsm.get_grad_logZ([1,1,1],[1.0,-1.0,1.0,0.0])
4-element Array{Float64,1}:
  0.7615941559702503
 -0.7615941559702503
  0.7615941559335818
 -1.1600513167692226
```
"""
function get_grad_logZ(n::Vector{Int64},θhat::Vector{Float64})::Vector{Float64}

    # Define function
    function f(θ::Vector{Float64})
        log(comp_Z(n,θ[1:(end-1)],θ[end]))
    end

    # Return ∇A(θ)
    return Calculus.gradient(f,θhat)

end
"""
    `check_boundary(θhat)`

Function that returns a bool indicating whether model with parameter estimate vector θhat is on the
boundary of the parameter space in any of its components.

# Examples
```julia-repl
julia> CpelAsm.check_boundary([1.0,1.0,1.0,1.0])
false
julia> CpelAsm.check_boundary([1.0,5.0,1.0,1.0])
true
```
"""
function check_boundary(θhat::Vector{Float64})::Bool

    # Return true θhat on boundary.
    return any(isapprox.(abs.(θhat),ETA_MAX_ABS;atol=5e-2)) || any(abs.(θhat).>ETA_MAX_ABS)

end
"""
    `comp_g([R1,...,RK],[α1,...,αK],β,αp1,αp2)`

Compute scaling factor in a model with [R1,...,RK] unobserved CpG cites from each block, with
parameters [α1,...,αK] and β. The transfer-matrix method is used to obtain a computationally
efficient expression without the need of recursive expressions. αp1 and αp2 are determined by the
respective kind of boundary.

# Examples
```julia-repl
julia> CpelAsm.comp_g([1],[1.0],1.0,3.0,3.0)
20.135323991555527
```
"""
function comp_g(r::Vector{Int64},a::Vector{Float64},b::Float64,ap1::Float64,ap2::Float64)::Float64

    # Return one if r=0
    r==0 && return 1.0

    # Boundary conditions
    y = get_u(ap1)'*get_W(r[1],a[1],b)
    if length(r)>1
        y *= prod([get_V(a[(i-1):i],b)*get_W(r[i],a[i],b) for i in 2:length(r)])
    end
    y *= get_u(ap2)

    # Return scaling factor
    return max(1e-100,y)

end
"""
    `comp_lkhd(X,[N1,...,NK],[α1,...,αK],β)`

Compute likelihood of a partial/full observation X that can be missing values anywhere (given by
0's).

# Examples
```julia-repl
julia> CpelAsm.comp_lkhd([1,1,0,1,1],[5],[1.0],1.0)
0.953646032691218
```
"""
function comp_lkhd(x::Vector{Int64},n::Vector{Int64},a::Vector{Float64},b::Float64)::Float64

    # Avoid the rest if N=1
    length(x)==1 && return 0.5 * exp(x[1]*a[1]) / cosh(a[1])

    # Get partition function and energy function
    Z = comp_Z(n,a,b)
    Ux = create_Ux(n,a,b)

    # Find changes to/from 0 in x vector
    zerost = findall(isequal(-1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1
    zeroend = findall(isequal(1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1

    # Determine whether it starts/finishes with 0
    x[1]==0 && pushfirst!(zerost,1)
    x[end]==0 && push!(zeroend,length(x)+1)

    # Find subregion label for each CpG site
    subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)

    # Get overall scaling factor as product of individual factors
    sf = 1.0
    @inbounds for i in 1:length(zerost)
        # Find boundary conditions (other X's or boundary of X)
        ap1 = zerost[i]==1 ? a[1] : 2.0*x[zerost[i]-1]*b+a[subid[zerost[i]]]
        ap2 = zeroend[i]==sum(n)+1 ? a[end] : 2.0*x[zeroend[i]]*b+a[subid[zeroend[i]-1]]

        # Figure out b (block IDs of missing) and r (α indices)
        b_id = subid[zerost[i]:(zeroend[i]-1)]
        n_miss = [count(x->x==n,b_id) for n in unique(b_id)]

        # Call scaling factor function
        sf *= comp_g(n_miss,a[unique(b_id)],b,ap1,ap2)

    end

    # Return energy function properly scaled.
    return exp(-Ux(x)) * sf / Z

end
"""
    `create_Llkhd([N1,...,NK],XOBS)`

Create function to compute the minus log-likelihood function for a region with N CpG sites given
the M partial observations XOBS.

# Examples
```julia-repl
julia> n=[2,2];
julia> xobs=CpelAsm.gen_ising_full_data(20,n);
julia> LogLike=CpelAsm.create_Llkhd(n,xobs);
julia> LogLike([0.0,0.0,0.0])
55.45177444479562
```
"""
function create_Llkhd(n::Vector{Int64},xobs::Array{Vector{Int64},1})

    # Define minus log-likelihood function
    function Llkhd_fun(theta::Vector{Float64})::Float64
        # Get parameters
        aux = 0.0
        a = theta[1:(end-1)]
        b = theta[end]

        # Get energy function and partition function
        Ux = create_Ux(n,a,b)
        logZ = log(comp_Z(n,a,b))

        # Initialize variables used in for loop
        ap1 = ap2 = 0.0                 # Boundary values left and right
        b_id = n_miss = []              # Block IDs of missing CpG; number missing per ID
        zerost = zeroend = [] # Indices of zeros
        subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)    # ID's of CpGs

        # Contribution of each observation
        for x in xobs

            # Find changes to/from 0 in x vector
            zerost = findall(isequal(-1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1
            zeroend = findall(isequal(1),abs.(x[2:end]) - abs.(x[1:(end-1)])) .+ 1

            # Determine whether it starts/finishes with 0
            x[1]==0 && pushfirst!(zerost,1)
            x[end]==0 && push!(zeroend,length(x)+1)

            # Get scaling factors due to all missing data
            @inbounds for i in 1:length(zerost)
                # Find boundary conditions (other X's or boundary of X)
                ap1 = zerost[i]==1 ? a[1] : 2.0*x[zerost[i]-1]*b+a[subid[zerost[i]]]
                ap2 = zeroend[i]==sum(n)+1 ? a[end] : 2.0*x[zeroend[i]]*b+a[subid[zeroend[i]-1]]

                # Figure out b (block IDs of missing) and r (α indices)
                b_id = subid[zerost[i]:(zeroend[i]-1)]
                n_miss = [count(x->x==n,b_id) for n in unique(b_id)]

                # Call scaling factor function
                aux += log(comp_g(n_miss,a[unique(b_id)],b,ap1,ap2))
            end

            # Add log of it to Llkhd
            aux += -Ux(x)
        end

        # Return MINUS log-likelihood. Add as many logZ as samples we have
        -aux + length(xobs) * logZ
    end

    # Return function
    return Llkhd_fun

end # end create_Llkhd
"""
    `est_alpha(XOBS)`

Estimate parameter α in N=1 case (first entry in returned vector).

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> n=[1];
julia> xobs=CpelAsm.gen_ising_full_data(100,n);
julia> CpelAsm.est_alpha(xobs)
2-element Array{Float64,1}:
 -0.020002667306849575
  0.0
```
"""
function est_alpha(xobs::Array{Vector{Int64},1})::Vector{Float64}
    # Derivation of estimate:
        # log(p)+log(2) = α-log(cosh(α)); log(1-p)+log(2) = -α-log(cosh(α))
        # α = (log(p)-log(1-p))/2

    # Proportion of X=1
    phat = length(findall(x->x==[1],xobs)) / length(xobs)
    a = 0.5*(log(phat)-log(1.0-phat))

    # Return estimate
    return [min(max(-ETA_MAX_ABS,a),ETA_MAX_ABS),0.0]

end # end est_alpha
"""
    `est_theta_sa([N1,...,NK],XOBS)`

Estimate parameter vector [α1,...,αK,β] using simulated annealing.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> n=[4]
julia> xobs=CpelAsm.gen_ising_full_data(100,n);
julia> CpelAsm.est_theta_sa(n,xobs)
2-element Vector{Float64}:
 -0.051768715945913854
 -0.022598858872278995
```
"""
function est_theta_sa(n::Vector{Int64},xobs::Array{Vector{Int64},1})::Vector{Float64}

    # If N=1, then estimate α
    sum(n)==1 && return est_alpha(xobs)

    # Define boundaries and initialization
    L = create_Llkhd(n,xobs)
    init = zeros(Float64,length(n)+1)
    lower = -ETA_MAX_ABS * ones(Float64,length(n)+1)
    upper = ETA_MAX_ABS * ones(Float64,length(n)+1)
    opts = Optim.Options(iterations=10^6,show_trace=false,store_trace=false)

    # Boxed Simulated Annealing (SAMI)
    optim = Optim.optimize(L,lower,upper,init,SAMIN(rt=1e-4;f_tol=1e-3,verbosity=0),opts)

    # Return estimate
    return optim.minimizer

end # end est_theta_sa
###################################################################################################
# COMPETING MODELS
###################################################################################################
"""
    `mle_mult(XOBS)`

Estimate parameter vector p=[p_1,...,p_{2^N}] based on a multinomial model and full observations.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=CpelAsm.gen_mult_full_data(100);
julia> CpelAsm.mle_mult(xobs)
16-element Vector{Float64}:
 0.07
 0.08
 0.06
 0.06
 0.08
 0.08
 0.04
 0.05
 0.04
 0.07
 0.07
 0.02
 0.04
 0.12
 0.07
 0.05
```
"""
function mle_mult(xobs::Array{Vector{Int64},1})::Vector{Float64}

    # Get xcal
    N = size(xobs[1])[1]
    xcal = generate_xcal(N)

    # Get empirical estimates
    phat = zeros(Float64,2^N)
    for x in xobs
        phat[findfirst(y-> x==y, xcal)] += 1
    end
    phat = phat ./ sum(phat)

    # Return phat
    return phat

end # end mle_mult
"""
    `mle_bin_semi(XOBS)`

Estimate parameter vector `p=[p_1,...,p_{N}]`, based on a binomial model that assumes a different
probability of methylation at each CpG site, from potentially partial observations. This results
in a semi-parametric model that grows in size proportionally to the number of CpG sites.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=CpelAsm.gen_mult_full_data(100);
julia> CpelAsm.mle_bin_semi(xobs)
4-element Array{Float64,1}:
 0.49
 0.51
 0.46
 0.44
```
"""
function mle_bin_semi(xobs::Array{Vector{Int64},1})::Vector{Float64}

    # Get xcal
    N = size(xobs[1])[1]

    # Get empirical estimate
    θ = [length(findall(x->x[i]==1,xobs)) / (length(findall(x->x[i]==1,xobs))+
         length(findall(x->x[i]==-1,xobs))) for i in 1:N]

    # Return phat
    return θ

end # end mle_bin_semi
"""
    `mle_bin_param(XOBS)`

Estimate parameter vector `p=[p_1,...,p_{2^N}]` based on a binomial model
that assumes that each CpG site has the same probability of being methylated,
as well as fully observed vectors. This results in a parametric model that does
not grow with the number of CpG sites considered.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> xobs=CpelAsm.gen_mult_full_data(100);
julia> CpelAsm.mle_bin_param(xobs)
16-element Array{Float64,1}:
 0.07596914062500001
 0.068733984375
 0.068733984375
 0.062187890625
 0.068733984375
 0.062187890625
 0.062187890625
 0.056265234374999994
 0.068733984375
 0.062187890625
 0.062187890625
 0.056265234374999994
 0.062187890625
 0.056265234374999994
 0.056265234374999994
 0.05090664062499999
```
"""
function mle_bin_param(xobs::Array{Vector{Int64},1})::Vector{Float64}

    # Get xcal
    N = size(xobs[1])[1]
    M = length(xobs)
    xcal = generate_xcal(N)

    # Get empirical estimates θ=[θ_1,...,θ_N]
    θ = sum([length(findall(x->x==1,x)) / N for x in xobs]) / M

    # Compute probabilities of each pattern
    i = 1
    phat = zeros(Float64,2^N)
    for x in xcal
        k = length(findall(y->y==1,x))
        phat[i] = θ^(k) * (1-θ)^(N-k)
        i+=1
    end

    # Return phat
    return phat

end # end mle_bin_param
###################################################################################################
# NOT CURRENTLY USED (EM Algorithm)
###################################################################################################
"""
    `keep_region(M,[N1,...,NK],θhat)`

Function that returns a bool indicating whether model with [N1,...,NK] CpG cites and with parameter
estimate vector θhat, from M observations, should be kept for downstream analysis. Not currently
used.

# Examples
```julia-repl
julia> CpelAsm.keep_region(10,[1,1,1],[1.0,1.0,1.0,1.0])
false
```
"""
function keep_region(m::Int64,n::Vector{Int64},θhat::Vector{Float64})::Bool

    # Compute FIM considering special case N=1
    fim = sum(n)==1 ? 1-(tanh(θhat[1]))^2 : get_fim(n,θhat)
    sum(n)==1 && return 4.0/sqrt(fim*m)<=CI_THRESH

    # Return true if all CIs are below threshold width
    return det(fim)!=0 ? all(4*sqrt.(max.(0.0,diag(inv(fim)))/m).<=CI_THRESH) : false

end
"""
    `get_fim([N1,...,NK],θhat)`

Estimate Fisher information matrix of a model with [N1,...,NK] CpG cites and evaluated at estimated
parameter vector θhat.

# Examples
```julia-repl
julia> CpelAsm.get_fim([1,1,1],[1.0,-1.0,1.0,0.0])
4×4 Array{Float64,2}:
  0.419975    -3.02772e-6  -3.02772e-6  -0.319849
 -3.02772e-6   0.419978     0.0          0.639701
 -3.02772e-6   0.0          0.419975    -0.319852
 -0.319849     0.639701    -0.319852     1.8143
```
"""
function get_fim(n::Vector{Int64},θhat::Vector{Float64})::Array{Float64,2}

    # Define function
    function f(θ::Vector{Float64})
        log(comp_Z(n,θ[1:(end-1)],θ[end]))
    end

    # Return I(θ)
    return Calculus.hessian(f,θhat)

end
"""
    `comp_ES(𝑋,[N1,...,N_K],θ)`

Function that computes expected value of the complete vector of sufficient statistics 𝐄[S(X)|𝑋;θ],
where X is the full methylation state, given the observed methylation vector 𝑋 and parameter
vector θ.

# Examples
```julia-repl
julia> CpelAsm.comp_ES([1,-1,0,-1],[2,2],[0.0,0.0,0.0])
3-element Array{Float64,1}:
  0.0
 -1.0
 -1.0
```
"""
function comp_ES(x::Vector{Int64},n::Vector{Int64},θ::Vector{Float64})::Vector{Float64}

    # Compute marginal of X=x
    p_x = comp_lkhd(x,n,θ[1:(end-1)],θ[end])

    # Compute E[W_n|X;θ]
    w = zeros(Int64,length(x))
    E_W_X = zeros(Float64,length(x))
    @inbounds for i=1:length(x)
        w[i] = 1
        E_W_X[i] = x[i]==0 ? (2.0*comp_lkhd(x+w,n,θ[1:(end-1)],θ[end])/p_x)-1.0 : x[i]
        w[i] = 0
    end

    # Compute E[W_{n}W_{n+1}|X;θ]
    w1 = zeros(Int64,length(x))
    w2 = zeros(Int64,length(x))
    E_WW_X = zeros(Float64,length(x)-1)
    @inbounds for i=1:(length(x)-1)
        if (x[i]!=0)&(x[i+1]!=0)
            E_WW_X[i] += x[i]*x[i+1]
            continue
        end
        if (x[i]!=0)&(x[i+1]==0)
            E_WW_X[i] += x[i]*E_W_X[i+1]
            continue
        end
        if (x[i]==0)&(x[i+1]!=0)
            E_WW_X[i] += E_W_X[i]*x[i+1]
            continue
        end
        if (x[i]==0)&(x[i+1]==0)
            w1[i] = -1
            w2[i+1] = -1
            E_WW_X[i] += comp_lkhd(x+w1+w2,n,θ[1:(end-1)],θ[end])
            w1[i] = 1
            w2[i+1] = 1
            E_WW_X[i] += comp_lkhd(x+w1+w2,n,θ[1:(end-1)],θ[end])
            w1[i] = -1
            E_WW_X[i] += -comp_lkhd(x+w1+w2,n,θ[1:(end-1)],θ[end])
            w1[i] = 0
            w2[i+1] = 0
            E_WW_X[i] *= 2.0/p_x
            E_WW_X[i] += -1.0
        end
    end

    # Get subregion ID for each CpG site
    subid = vcat([i*ones(Int64,n[i]) for i in 1:length(n)]...)

    # Initialize with observed part for k=1,...,K
    Sout = [sum(E_W_X[subid.==id]) for id in unique(subid)]

    # Append K+1 component
    push!(Sout,sum(E_WW_X))

    # Return vector E[S(W)|X;θ]
    return Sout

end # end comp_ES
"""
    `comp_ET(XOBS,[N1,...,N_K],θ)`

Function that computes expectation

    𝐄[T(X_1,...,X_M)|𝑋_1,...,𝑋_M;θ],

where X_i is the complete methylation state that corresponds to the i-th observation, 𝑋_i is the
corresponding observed vector, and θ is the parameter vector the expectation is computed wrt.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> n=[2,2];
julia> xobs=CpelAsm.gen_ising_full_data(100,n);
julia> CpelAsm.comp_ET(xobs,n,[0.0,0.0,0.0])
3-element Array{Float64,1}:
   0.0
 -20.0
  -6.0
```
"""
function comp_ET(xobs::Array{Vector{Int64},1},n::Vector{Int64},θ::Vector{Float64})::Vector{Float64}

    # Return vector E[T({W1,...,Wm})|{X1,...,Xm};θ]
    return sum(map(x -> comp_ES(x,n,θ),xobs))

end # end comp_ET
"""
    `em_alg([N1,...,N_K],XOBS)`

Function that runs a single instance of the EM algorithm given observations in XOBS. Returns a pair
`θhat` and `convergence` flag.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> n=[2,2];
julia> xobs=CpelAsm.gen_ising_full_data(100,n);
julia> CpelAsm.em_alg(n,xobs)
([-0.00114078, -0.10281, -0.0235107], true)
```
"""
function em_alg(n::Vector{Int64},xobs::Array{Vector{Int64},1})::Tuple{Vector{Float64},Bool}

    # Check if we have complete data
    complete = isnothing(findfirst(x->any(x.==0),xobs))

    # Initialize θhat
    θhat = rand(Uniform(-1.5,1.5),length(n)+1)

    # Iterate until ||θ-θhat||<ϵ
    ϵ = 5e-3
    em_convergence = false
    @inbounds for i=1:10
        # (Re)-compute E[T|X1,...,XM;θ^{(i-1)}]
        ET = comp_ET(xobs,n,θhat)
        # Set up system
        function f!(U,θ)
            ∇logZ = get_grad_logZ(n,θ)
            @inbounds for j=1:(length(n)+1)
                U[j] =  ∇logZ[j] - ET[j] / length(xobs)
            end
        end
        #  Try to solve NL system
	    sol =
        try
	        # Call solver
	        nlsolve(f!,θhat;iterations=20,ftol=1e-3)
        catch x
	        # Report x error if found
	        print_log("Issue with NL solver.")
	        break
        end
        # Leave if no convergence
        converged(sol) || break
        # Check Euclidean distance
        em_convergence = complete || (sqrt(sum((sol.zero-θhat).^2)/sum(θhat.^2))<=ϵ)
        # Update θhat
        θhat = sol.zero
        # Break if converged
        em_convergence && break
    end

    # Return θhat
    return (θhat,em_convergence)

end # end em_alg
"""
    `est_theta_em([N1,...,N_K],XOBS)`

Function that runs multiple instances of the EM algorithm given observations in XOBS. Returns a pair
(θhat,convergence flag). This is done so as to try different initializations for the EM algorithm.
The θhat that maximizes the log-likelihood is returned as long as at least one instance of the EM
algorithm has converged.

# Examples
```julia-repl
julia> Random.seed!(1234);
julia> n=[2,2];
julia> xobs=CpelAsm.gen_ising_full_data(100,n);
julia> CpelAsm.est_theta_em(n,xobs)
([-0.00114078, -0.10281, -0.0235107], true)
```
"""
function est_theta_em(n::Vector{Int64},xobs::Array{Vector{Int64},1})::Tuple{Vector{Float64},Bool}

    # If N=1, then estimate α
    sum(n)==1 && return (est_alpha(xobs),true)

    # Check if we have complete data
    complete = isnothing(findfirst(x->any(x.==0),xobs))

    # Initialize
    min_LogLike = Inf
    convergence = false
    LogLike = create_Llkhd(n,xobs)
    θhat = zeros(Float64,length(n)+1)

    # Run EM with a number of different initializations
    @inbounds for i=1:15
        # EM Algorithm: obtain candidate and convergence flag
        θcand,converged = em_alg(n,xobs)
        # If did not converge move onto new initialization
        converged || continue
        # If it converged and minimizes minus log-likelihood, then keep estimate
        if LogLike(θcand)<min_LogLike
            θhat = θcand
            min_LogLike = min_LogLike
        end
        # Set convergence true since at least one init worked
        convergence = true
        # If full data then no need to re-run EM
        complete && break
    end

    # Return θhat
    return (θhat,convergence)

end # end est_theta_em


text:
"""
Test code for aerodynamic coefficients of a straight wing
"""

using ForwardDiff

@test begin
	NW = 65
	NM = 10
	NF = 21

	missile_xs = collect(LinRange(-0.5, 0.0, NM))
	fus_xs = collect(LinRange(- 3.5, 3.5, NF))

	fus_Rs = sqrt.(1.0 .- (fus_xs ./ 3.5) .^ 2) .* 0.2
	fus_As = fus_Rs .^ 2 .* π

	fus_z = 0.1

	pts = hcat(
		[
			[0.0, y, 0.0] for y in LinRange(- 3.5, 3.5, NW)
		]...
	)

	pts = hcat(
		pts,
		[
			[missile_xs[id], - 3.5, 0.0] for id = 1:(NM - 1)
		]...,
		[
			[missile_xs[id], 3.5, 0.0] for id = 1:(NM - 1)
		]...,
		[
			[fus_xs[id], 0.0, fus_z] for id = 1:NF
		]...
	)

	acft = Aircraft(pts)

	for i = 1:(NW - 1)
		add_wing_strip!(
			acft,
			i,
			i + 1,
			1.0;
			e = 0.25,
			m = 0.2
		)
	end

	AM = 0.01 ^ 2 * π

	for i = 1:(NM - 1)
		add_fuselage_cut!(
			acft,
			i + NW,
			(i == (NM - 1) ? 1 : i + NW + 1),
			(i == 1 ? 0.0 : AM),
			(i == NM - 1 ? 0.0 : AM)
		)

		add_fuselage_cut!(
			acft,
			i + NW + (NM - 1),
			(i == (NM - 1) ? NW : i + NW + NM),
			(i == 1 ? 0.0 : AM),
			(i == NM - 1 ? 0.0 : AM)
		)
	end

	for i = 1:(NF - 1)
		add_fuselage_cut!(
			acft,
			NW + 2 * (NM - 1) + i,
			NW + 2 * (NM - 1) + i + 1,
			fus_As[i],
			fus_As[i + 1]
		)
	end

	add_link!(
		acft,
		NW ÷ 2 + 1,
		NW + 2 * (NM - 1) + NF ÷ 2 + 1
	)

	add_point_drag!(
		acft,
		NW ÷ 2 + 1,
		4 * π * 0.025 ^ 3 / 3
	)

	add_engine!(
		acft,
		NW + 1,
		1.0
	)
	add_engine!(
		acft,
		NW + NM,
		1.0
	)

	fcon = get_state(
		acft,
		100.0;
		α = 5.0
	)

	FM = [
		100.0 100.0
		(- 10.0) 10.0
	]

	@btime fcon_derivative, coefs = state_space($acft, $fcon, 100.0; α = 5.0, Sref = 7.0, bref = 7.0, engine_loads = $FM)
	fcon_derivative, coefs = state_space(acft, fcon, 100.0; α = 5.0, Sref = 7.0, bref = 7.0, engine_loads = FM)

	# plot_aircraft(acft)

	true
end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
module Tmp

using Test
using Mimi
import Mimi:
    compdefs, compdef, external_param_conns

@defcomp X begin
    x = Parameter(index = [time])
    y = Variable(index = [time])
    function run_timestep(p, v, d, t)
        v.y[t] = 1
    end
end

@defcomp X_repl begin
    x = Parameter(index = [time])
    y = Variable(index = [time])
    function run_timestep(p, v, d, t)
        v.y[t] = 2
    end
end

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X, exports=[:x => :z])          # Original component X
add_comp!(m, X_repl)
set_param!(m, :X, :x, zeros(6))

if false
    run(m)
    @test m[:X, :y] == ones(6)

    replace_comp!(m, X_repl, :X)
    run(m)

    @test length(components(m)) == 1        # Only one component exists in the model
    @test m[:X, :y] == 2 * ones(6)          # Successfully ran the run_timestep function from X_repl
end

end # module

using Mimi
m = Tmp.m


text:
using Documenter, MendelIHT

makedocs(
    doctest = false, #will set to true later
    format = Documenter.HTML(),
    sitename = "MendelIHT",
    authors = "Benjamin Chu, Kevin Keys",
    clean = true,
    pages = [
        "Home"            => "index.md",
        "Getting Started" => "man/getting_started.md",
        "Examples"        => "man/examples.md",
        "Mathematical Details" => "man/math.md",
        "Contributing"    => "man/contributing.md",
        "FAQ"             => "man/FAQ.md",
        "API"             => "man/api.md",
    ]
)

deploydocs(
    repo   = "github.com/OpenMendel/MendelIHT.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to print the tables.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export pretty_table

################################################################################
#                               Public Functions
################################################################################

"""
    pretty_table([io::IO | String | HTML,] table;  kwargs...)

Print to `io` the `table`.

If `io` is omitted, then it defaults to `stdout`. If `String` is passed in the
place of `io`, then a `String` with the printed table will be returned by the
function. If `HTML` is passed in the place of `io`, then an `HTML` object is
returned with the printed table.

When printing, it will be verified if `table` complies with **Tables.jl** API.
If it is compliant, then this interface will be used to print the table. If it
is not compliant, then only the following types are supported:

1. `AbstractVector`: any vector can be printed.
2. `AbstractMatrix`: any matrix can be printed.
3. `Dict`: any `Dict` can be printed. In this case, the special keyword
    `sortkeys` can be used to select whether or not the user wants to print the
    dictionary with the keys sorted. If it is `false`, then the elements will be
    printed on the same order returned by the functions `keys` and `values`.
    Notice that this assumes that the keys are sortable, if they are not, then
    an error will be thrown.

# Keywords

- `alignment::Union{Symbol, Vector{Symbol}}`: Select the alignment of the
    columns (see the section `Alignment`).
- `backend::Union{Symbol, T_BACKENDS}`: Select which back-end will be used to
    print the table (see the section `Backend`). Notice that the
    additional configuration in `kwargs...` depends on the selected backend.
- `cell_alignment::Union{Nothing, Dict{Tuple{Int, Int}, Symbol}, Function, Tuple}`:
    A tuple of functions with the signature `f(data, i, j)` that overrides the
    alignment of the cell `(i, j)` to the value returned by `f`. It can also be a
    single function, when it is assumed that only one alignment function is
    required, or `nothing`, when no cell alignment modification will be
    performed. If the function `f` does not return a valid alignment symbol as
    shown in section `Alignment`, then it will be discarded. For
    convenience, it can also be a dictionary of type `(i, j) => a` that
    overrides the alignment of the cell `(i, j)` to `a`. `a` must be a symbol
    like specified in the section `Alignment`. (**Default** = `nothing`)

!!! note
    If more than one alignment function is passed to `cell_alignment`, then
    the functions will be evaluated in the same order of the tuple. The
    first one that returns a valid alignment symbol for each cell is applied,
    and the rest is discarded.

- `cell_first_line_only::Bool`: If `true`, then only the first line of each cell
    will be printed. (**Default** = `false`)
- `compact_printing::Bool`: Select if the option `:compact` will be used when
    printing the data. (**Default** = `true`)
- `filters_row::Union{Nothing, Tuple}`: Filters for the rows (see the section
    `Filters`).
- `filters_col::Union{Nothing, Tuple}`: Filters for the columns (see the section
    `Filters`).
- `formatters::Union{Nothing, Function, Tuple}`: See the section
    `Formatters`.
- `header::Union{Symbol, Vector{Symbol}}`: The header must be a tuple of
    vectors. Each one must have the number of elements equal to the number of
    columns in the table. The first vector is considered the header and the
    others are the subheaders. If it is `nothing`, then a default value based on
    the type will be used. If a single vector is passed, then it will be
    considered the header. (**Default** = `nothing`)
- `header_alignment::Union{Symbol, Vector{Symbol}}`: Select the alignment of the
    header columns (see the section `Alignment`). If the symbol that
    specifies the alignment is `:s` for a specific column, then the same
    alignment in the keyword `alignment` for that column will be used.
    (**Default** = `:s`)
- `header_cell_alignment::Union{Nothing, Dict{Tuple{Int, Int}, Symbol}, Function, Tuple}`:
    This keyword has the same structure of `cell_alignment` but in this case it
    operates in the header. Thus, `(i, j)` will be a cell in the header matrix
    that contains the header and sub-headers. This means that the `data` field
    in the functions will be the same value passed in the keyword `header`.
    (**Default** = `nothing`)

!!! note
      If more than one alignment function is passed to `header_cell_alignment`,
      then the functions will be evaluated in the same order of the tuple. The
      first one that returns a valid alignment symbol for each cell is applied,
      and the rest is discarded.

- `limit_printing::Bool`: If `true`, then the cells will be converted using the
    property `:limit => true` of `IOContext`. (**Default** = `true`)
- `renderer::Symbol`: A symbol that indicates which function should be used to
    convert an object to a string. It can be `:print` to use the function
    `print` or `:show` to use the function `show`. Notice that this selection is
    applicable only to the table data. Headers, sub-headers, and row name column
    are always rendered with print. (**Default** = `:print`)
- `row_names::Union{Nothing, AbstractVector}`: A vector containing the row names
    that will be appended to the left of the table. If it is `nothing`, then the
    column with the row names will not be shown. Notice that the size of this
    vector must match the number of rows in the table. (**Default** = `nothing`)
- `row_name_alignment::Symbol`: Alignment of the column with the rows name (see
    the section `Alignment`).
- `row_name_column_title::AbstractString`: Title of the column with the row
    names. (**Default** = "")
- `row_number_column_title::AbstractString`: Title of the column with the row
    numbers. (**Default** = "Row")
- `show_row_number::Bool`: If `true`, then a new column will be printed showing
    the row number. (**Default** = `false`)
- `title::AbstractString`: The title of the table. If it is empty, then no title
    will be printed. (**Default** = "")
- `title_alignment::Symbol`: Alignment of the title, which must be a symbol as
    explained in the section `Alignment`. This argument is ignored in the
    LaTeX backend. (**Default** = :l)

!!! note
    Notice that all back-ends have the keyword `tf` to specify the table
    printing format. Thus, if the keyword `backend` is not present or if it is
    `nothing`, then the back-end will be automatically inferred from the type of
    the keyword `tf`. In this case, if `tf` is also not present, then it just
    fall-back to the text back-end unless `HTML` is passed as the first
    argument. In this case, the default back-end is set to HTML.

# Alignment

The keyword `alignment` can be a `Symbol` or a vector of `Symbol`.

If it is a symbol, we have the following behavior:

- `:l` or `:L`: the text of all columns will be left-aligned;
- `:c` or `:C`: the text of all columns will be center-aligned;
- `:r` or `:R`: the text of all columns will be right-aligned;
- Otherwise it defaults to `:r`.

If it is a vector, then it must have the same number of symbols as the number of
columns in `data`. The *i*-th symbol in the vector specify the alignment of the
-i*-th column using the same symbols as described previously.

!!! note
    In HTML backend, the user can select `:n` ou `:N` to print the cell without
    any alignment annotation.

# Filters

It is possible to specify filters to filter the data that will be printed. There
are two types of filters: the row filters, which are specified by the keyword
`filters_row`, and the column filters, which are specified by the keyword
`filters_col`.

The filters are a tuple of functions that must have the following signature:

```julia
f(data,i)::Bool
```

in which `data` is a pointer to the matrix that is being printed and `i` is the
i-th row in the case of the row filters or the i-th column in the case of column
filters. If this function returns `true` for `i`, then the i-th row (in case of
`filters_row`) or the i-th column (in case of `filters_col`) will be printed.
Otherwise, it will be omitted.

A set of filters can be passed inside of a tuple. Notice that, in this case,
-*all filters** for a specific row or column must be return `true` so that it
can be printed, *i.e* the set of filters has an `AND` logic.

If the keyword is set to `nothing`, which is the default, then no filtering will
be applied to the data.

!!! note
    The filters do not change the row and column numbering for the others
    modifiers such as column width specification, formatters, and highlighters.
    Thus, for example, if only the 4-th row is printed, then it will also be
    referenced inside the formatters and highlighters as 4 instead of 1.

---

# Pretty table text back-end

This back-end produces text tables. This back-end can be used by selecting
`back-end = :text`.

# Keywords

- `border_crayon::Crayon`: Crayon to print the border.
- `header_crayon::Union{Crayon, Vector{Crayon}}`: Crayon to print the header.
- `subheader_crayon::Union{Crayon, Vector{Crayon}}`: Crayon to print
    sub-headers.
- `rownum_header_crayon::Crayon`: Crayon for the header of the column with the
    row numbers.
- `text_crayon::Crayon`: Crayon to print default text.
- `omitted_cell_summary_crayon::Crayon`: Crayon used to print the omitted cell
    summary.
- `alignment_anchor_fallback::Symbol`: This keyword controls the line alignment
    when using the regex alignment anchors if a match is not found. If it is
    `:l`, then the left of the line will be aligned with the anchor. If it is
    `:c`, then the line center will be aligned with the anchor. Otherwise, the
    end of the line will be aligned with the anchor. (**Default** = `:l`)
- `alignment_anchor_fallback_override::Dict{Int, Symbol}`: A `Dict{Int, Symbol}`
    to override the behavior of `fallback_alignment_anchor` for a specific
    column. Example: `Dict(3 => :c)` changes the fallback alignment anchor
    behavior for `:c` only for the column 3.
- `alignment_anchor_regex::Dict{Int, AbstractVector{Regex}}`: A dictionary
    `Dict{Int, AbstractVector{Regex}}` with a set of regexes that is used to
    align the values in the columns (keys). The characters at the first regex
    match (or anchor) of each line in every cell of the column will be aligned.
    The regex match is searched in the same order as the regexes appear on the
    vector. The regex matching is applied after the cell conversion to string,
    which includes the formatters. If no match is found for a specific line,
    then the alignment of this line depends on the options
    `alignment_anchor_fallback` and `alignment_anchor_fallback_override`. If the
    key `0` is present, then the related regexes will be used to align all the
    columns. In this case, all the other keys will be neglected. Example:
    `Dict(2 => [r"\\."])` aligns the decimal point of the cells in the second
    column. (**Default** = `Dict{Int, Vector{Regex}}()`)
- `autowrap::Bool`: If `true`, then the text will be wrapped on spaces to fit
    the column. Notice that this function requires `linebreaks = true` and the
    column must have a fixed size (see `columns_width`).
- `body_hlines::Vector{Int}`: A vector of `Int` indicating row numbers in which
    an additional horizontal line should be drawn after the row. Notice that
    numbers lower than 1 and equal or higher than the number of printed rows
    will be neglected. This vector will be appended to the one in `hlines`, but
    the indices here are related to the printed rows of the body. Thus, if `1`
    is added to `body_hlines`, then a horizontal line will be drawn after the
    first data row. (**Default** = `Int[]`)
- `body_hlines_format::Union{Nothing, NTuple{4, Char}}`: A tuple of 4 characters
    specifying the format of the horizontal lines that will be drawn by
    `body_hlines`. The characters must be the left intersection, the middle
    intersection, the right intersection, and the row. If it is `nothing`, then
    it will use the same format specified in `tf`. (**Default** = `nothing`)
- `columns_width::Union{Int, AbstractVector{Int}}`: A set of integers specifying
    the width of each column. If the width is equal or lower than 0, then it
    will be automatically computed to fit the large cell in the column. If it is
    a single integer, then this number will be used as the size of all columns.
    (**Default** = 0)
- `crop::Symbol`: Select the printing behavior when the data is bigger than the
    available display size (see `display_size`). It can be `:both` to crop on
    vertical and horizontal direction, `:horizontal` to crop only on horizontal
    direction, `:vertical` to crop only on vertical direction, or `:none` to do
    not crop the data at all. If the `io` has `:limit => true`, then `crop` is
    set to `:both` by default. Otherwise, it is set to `:none` by default.
- `crop_num_lines_at_beginning::Int`: Number of lines to be left at the
    beginning of the printing when vertically cropping the output. Notice that
    the lines required to show the title are automatically computed.
    (**Default** = 0)
- `crop_subheader::Bool`: If `true`, then the sub-header size will not be taken
    into account when computing the column size. Hence, the print algorithm can
    crop it to save space. This has no effect if the user selects a fixed column
    width. (**Default** = `false`)
- `continuation_row_alignment::Symbol`: A symbol that defines the alignment of
    the cells in the continuation row. This row is printed if the table is
    vertically cropped. (**Default** = `:c`)
- `display_size::Tuple{Int, Int}`: A tuple of two integers that defines the
    display size (num. of rows, num. of columns) that is available to print the
    table. It is used to crop the data depending on the value of the keyword
    `crop`. Notice that if a dimension is not positive, then it will be treated
    as unlimited. (**Default** = `displaysize(io)`)
- `ellipsis_line_skip::Integer`: An integer defining how many lines will be
    skipped from showing the ellipsis that indicates the text was cropped.
    (**Default** = 0)
- `equal_columns_width::Bool`: If `true`, then all the columns will have the
    same width. (**Default** = `false`)
- `highlighters::Union{Highlighter, Tuple}`: An instance of `Highlighter` or a
    tuple with a list of text highlighters (see the section `Text
    highlighters`).
- `hlines::Union{Nothing, Symbol, AbstractVector}`: This variable controls where
    the horizontal lines will be drawn. It can be `nothing`, `:all`, `:none` or
    a vector of integers. (**Default** = `nothing`)
    - If it is `nothing`, which is the default, then the configuration will be
        obtained from the table format in the variable `tf` (see
        [`TextFormat`](@ref)).
    - If it is `:all`, then all horizontal lines will be drawn.
    - If it is `:none`, then no horizontal line will be drawn.
    - If it is a vector of integers, then the horizontal lines will be drawn
        only after the rows in the vector. Notice that the top line will be
        drawn if `0` is in `hlines`, and the header and subheaders are
        considered as only 1 row. Furthermore, it is important to mention that
        the row number in this variable is related to the **printed rows**.
        Thus, it is affected by filters, and by the option to suppress the
        header `noheader`. Finally, for convenience, the top and bottom lines
        can be drawn by adding the symbols `:begin` and `:end` to this vector,
        respectively, and the line after the header can be drawn by adding the
        symbol `:header`.

!!! info
    The values of `body_hlines` will be appended to this vector. Thus,
    horizontal lines can be drawn even if `hlines` is `:none`.

- `linebreaks::Bool`: If `true`, then `\\n` will break the line inside the
    cells. (**Default** = `false`)
- `maximum_columns_width::Union{Int, AbstractVector{Int}}`: A set of integers
    specifying the maximum width of each column. If the width is equal or lower
    than 0, then it will be ignored. If it is a single integer, then this number
    will be used as the maximum width of all columns. Notice that the parameter
    `columns_width` has precedence over this one. (**Default** = 0)
- `minimum_columns_width::Union{Int, AbstractVector{Int}}`: A set of integers
    specifying the minimum width of each column. If the width is equal or lower
    than 0, then it will be ignored. If it is a single integer, then this number
    will be used as the minimum width of all columns. Notice that the parameter
    `columns_width` has precedence over this one. (**Default** = 0)
- `newline_at_end::Bool`: If `false`, then the table will not end with a newline
    character. (**Default** = `true`)
- `noheader::Bool`: If `true`, then the header will not be printed. Notice that
    all keywords and parameters related to the header and sub-headers will be
    ignored. (**Default** = `false`)
- `nosubheader::Bool`: If `true`, then the sub-header will not be printed,
    *i.e.* the header will contain only one line. Notice that this option has no
    effect if `noheader = true`. (**Default** = `false`)
- `overwrite::Bool`: If `true`, then the same number of lines in the printed
    table will be deleted from the output `io`. This can be used to update the
    table in the display continuously. (**Default** = `false`)
- `row_number_alignment::Symbol`: Select the alignment of the row number column
    (see the section `Alignment`). (**Default** = `:r`)
- `show_omitted_cell_summary::Bool`: If `true`, then a summary will be printed
    after the table with the number of columns and rows that were omitted.
    (**Default** = `true`)
- `tf::TextFormat`: Table format used to print the table (see
    [`TextFormat`](@ref)). (**Default** = `tf_unicode`)
- `title_autowrap::Bool`: If `true`, then the title text will be wrapped
    considering the title size. Otherwise, lines larger than the title size will
    be cropped. (**Default** = `false`)
- `title_crayon::Crayon`: Crayon to print the title.
- `title_same_width_as_table::Bool`: If `true`, then the title width will match
    that of the table. Otherwise, the title size will be equal to the display
    width. (**Default** = `false`)
- `vcrop_mode::Symbol`: This variable defines the vertical crop behavior. If it
    is `:bottom`, then the data, if required, will be cropped in the bottom. On
    the other hand, if it is `:middle`, then the data will be cropped in the
    middle if necessary. (**Default** = `:bottom`)
- `vlines::Union{Nothing, Symbol, AbstractVector}`: This variable controls where
    the vertical lines will be drawn. It can be `nothing`, `:all`, `:none` or a
    vector of integers. (**Default** = `nothing`)
    - If it is `nothing`, which is the default, then the configuration will be
        obtained from the table format in the variable `tf` (see
        [`TextFormat`](@ref)).
    - If it is `:all`, then all vertical lines will be drawn.
    - If it is `:none`, then no vertical line will be drawn.
    - If it is a vector of integers, then the vertical lines will be drawn only
        after the columns in the vector. Notice that the top line will be drawn
        if `0` is in `vlines`. Furthermore, it is important to mention that the
        column number in this variable is related to the **printed column**.
        Thus, it is affected by filters, and by the options `row_names` and
        `show_row_number`. Finally, for convenience, the left and right vertical
        lines can be drawn by adding the symbols `:begin` and `:end` to this
        vector, respectively, and the line after the header can be drawn by
        adding the symbol `:header`.

The keywords `header_crayon` and `subheader_crayon` can be a `Crayon` or a
`Vector{Crayon}`. In the first case, the `Crayon` will be applied to all the
elements. In the second, each element can have its own crayon, but the length of
the vector must be equal to the number of columns in the data.

## Crayons

A `Crayon` is an object that handles a style for text printed on terminals. It
is defined in the package
[Crayons.jl](https://github.com/KristofferC/Crayons.jl). There are many options
available to customize the style, such as foreground color, background color,
bold text, etc.

A `Crayon` can be created in two different ways:

```julia-repl
julia> Crayon(foreground = :blue, background = :black, bold = :true)

julia> crayon"blue bg:black bold"
```

For more information, see the package documentation.

## Text highlighters

A set of highlighters can be passed as a `Tuple` to the `highlighters` keyword.
Each highlighter is an instance of the structure `Highlighter` that contains
three fields:

- `f::Function`: Function with the signature `f(data, i, j)` in which should
    return `true` if the element `(i, j)` in `data` must be highlighter, or
    `false` otherwise.
- `fd::Function`: Function with the signature `f(h,data,i,j)` in which `h` is
    the highlighter. This function must return the `Crayon` to be applied to the
    cell that must be highlighted.
- `crayon::Crayon`: The `Crayon` to be applied to the highlighted cell if the
    default `fd` is used.

The function `f` has the following signature:

    f(data, i, j)

in which `data` is a reference to the data that is being printed, and `i` and
`j` are the element coordinates that are being tested. If this function returns
`true`, then the cell `(i, j)` will be highlighted.

If the function `f` returns true, then the function `fd(h, data, i, j)` will be
called and must return a `Crayon` that will be applied to the cell.

A highlighter can be constructed using three helpers:

    Highlighter(f::Function; kwargs...)

where it will construct a `Crayon` using the keywords in `kwargs` and apply it
to the highlighted cell,

    Highlighter(f::Function, crayon::Crayon)

where it will apply the `crayon` to the highlighted cell, and

    Highlighter(f::Function, fd::Function)

where it will apply the `Crayon` returned by the function `fd` to the
highlighted cell.

!!! info
    If only a single highlighter is wanted, then it can be passed directly to
    the keyword `highlighter` without being inside a `Tuple`.

!!! note
    If multiple highlighters are valid for the element `(i, j)`, then the
    applied style will be equal to the first match considering the order in the
    tuple `highlighters`.

!!! note
    If the highlighters are used together with [Formatters](@ref), then the
    change in the format **will not** affect the parameter `data` passed to the
    highlighter function `f`. It will always receive the original, unformatted
    value.

---

# Pretty table HTML backend

This backend produces HTML tables. This backend can be used by selecting
`backend = Val(:html)`.

# Keywords

- `allow_html_in_cells::Bool`: By default, special characters like `<`, `>`,
    `"`, etc. are replaced in HTML backend to generate valid code. However, this
    algorithm blocks the usage of HTML code inside of the cells. If this keyword
    is `true`, then the escape algorithm **will not** be applied, allowing HTML
    code inside all the cells. In this case, the user must ensure that the
    output code is valid. (**Default** = `false`)
- `highlighters::Union{HTMLHighlighter, Tuple}`: An instance of
    [`HTMLHighlighter`](@ref) or a tuple with a list of HTML highlighters (see
    the section `HTML highlighters`).
- `linebreaks::Bool`: If `true`, then `\\n` will be replaced by `<br>`.
    (**Default** = `false`)
- `minify::Bool`: If `true`, then output will be displayed minified, *i.e.*
    without unnecessary indentation or newlines. (**Default** = `false`)
- `noheader::Bool`: If `true`, then the header will not be printed. Notice that
    all keywords and parameters related to the header and sub-headers will be
    ignored. (**Default** = `false`)
- `nosubheader::Bool`: If `true`, then the sub-header will not be printed,
    *i.e.* the header will contain only one line. Notice that this option has no
    effect if `noheader = true`. (**Default** = `false`)
- `standalone::Bool`: If `true`, then a complete HTML page will be generated.
    Otherwise, only the content between the tags `<table>` and `</table>` will
    be printed (with the tags included). (**Default** = `true`)
- `tf::HTMLTableFormat`: An instance of the structure [`HTMLTableFormat`](@ref)
    that defines the general format of the HTML table.

## HTML highlighters

A set of highlighters can be passed as a `Tuple` to the `highlighters` keyword.
Each highlighter is an instance of the structure [`HTMLHighlighter`](@ref). It
contains the following two public fields:

- `f::Function`: Function with the signature `f(data, i, j)` in which should
    return `true` if the element `(i,j)` in `data` must be highlighted, or
    `false` otherwise.
- `fd::Function`: Function with the signature `f(h, data, i, j)` in which `h` is
    the highlighter. This function must return the `HTMLDecoration` to be
    applied to the cell that must be highlighted.

The function `f` has the following signature:

    f(data, i, j)

in which `data` is a reference to the data that is being printed, and `i` and
`j` are the element coordinates that are being tested. If this function returns
`true`, then the highlight style will be applied to the `(i, j)` element.
Otherwise, the default style will be used.

If the function `f` returns true, then the function `fd(h, data, i, j)` will be
called and must return an element of type [`HTMLDecoration`](@ref) that contains
the decoration to be applied to the cell.

A HTML highlighter can be constructed using two helpers:

    HTMLHighlighter(f::Function, decoration::HTMLDecoration)

    HTMLHighlighter(f::Function, fd::Function)

The first will apply a fixed decoration to the highlighted cell specified in
`decoration` whereas the second let the user select the desired decoration by
specifying the function `fd`.

!!! info
    If only a single highlighter is wanted, then it can be passed directly to
    the keyword `highlighter` without being inside a `Tuple`.

!!! note
    If multiple highlighters are valid for the element `(i, j)`, then the
    applied style will be equal to the first match considering the order in the
    tuple `highlighters`.

!!! note
    If the highlighters are used together with [Formatters](@ref), then the
    change in the format **will not** affect the parameter `data` passed to the
    highlighter function `f`. It will always receive the original, unformatted
    value.

---

# Pretty table LaTeX backend

This backend produces LaTeX tables. This backend can be used by selecting
`backend = Val(:latex)`.

# Keywords

- `body_hlines::Vector{Int}`: A vector of `Int` indicating row numbers in which
    an additional horizontal line should be drawn after the row. Notice that
    numbers lower than 1 and equal or higher than the number of printed rows
    will be neglected. This vector will be appended to the one in `hlines`, but
    the indices here are related to the printed rows of the body. Thus, if `1`
    is added to `body_hlines`, then a horizontal line will be drawn after the
    first data row. (**Default** = `Int[]`)
- `highlighters::Union{LatexHighlighter, Tuple}`: An instance of
    `LatexHighlighter` or a tuple with a list of LaTeX highlighters (see the
    section `LaTeX highlighters`).
- `hlines::Union{Nothing, Symbol, AbstractVector}`: This variable controls where
    the horizontal lines will be drawn. It can be `nothing`, `:all`, `:none` or
    a vector of integers. (**Default** = `nothing`)
    - If it is `nothing`, which is the default, then the configuration will be
        obtained from the table format in the variable `tf` (see
        [`LatexTableFormat`](@ref)).
    - If it is `:all`, then all horizontal lines will be drawn.
    - If it is `:none`, then no horizontal line will be drawn.
    - If it is a vector of integers, then the horizontal lines will be drawn
        only after the rows in the vector. Notice that the top line will be
        drawn if `0` is in `hlines`, and the header and subheaders are
        considered as only 1 row. Furthermore, it is important to mention that
        the row number in this variable is related to the **printed rows**.
        Thus, it is affected by filters, and by the option to suppress the
        header `noheader`. Finally, for convenience, the top and bottom lines
        can be drawn by adding the symbols `:begin` and `:end` to this vector,
        respectively, and the line after the header can be drawn by adding the
        symbol `:header`.

!!! info
    The values of `body_hlines` will be appended to this vector. Thus,
    horizontal lines can be drawn even if `hlines` is `:none`.

- `label::AbstractString`: The label of the table. If empty, then no label will
    be added. (**Default** = "")
- `longtable_footer::Union{Nothing, AbstractString}`: The string that will be
    drawn in the footer of the tables before a page break. This only works if
    `table_type` is `:longtable`. If it is `nothing`, then no footer will be
    used. (**Default** = `nothing`)
- `noheader::Bool`: If `true`, then the header will not be printed. Notice that
    all keywords and parameters related to the header and sub-headers will be
    ignored. (**Default** = `false`)
- `nosubheader::Bool`: If `true`, then the sub-header will not be printed,
    *i.e.* the header will contain only one line. Notice that this option has no
    effect if `noheader = true`. (**Default** = `false`)
- `row_number_alignment::Symbol`: Select the alignment of the row number column
    (see the section `Alignment`). (**Default** = `:r`)
- `table_type::Union{Nothing, Symbol}`: Select which LaTeX environment will be
    used to print the table. Currently supported options are `:tabular` for
    `tabular` or `:longtable` for `longtable`. If it is `nothing` then the
    default option of the table format will be used. (**Default** = `nothing`)
- `tf::LatexTableFormat`: An instance of the structure
    [`LatexTableFormat`](@ref) that defines the general format of the LaTeX table.
- `vlines::Union{Nothing, Symbol, AbstractVector}`: This variable controls where
    the vertical lines will be drawn. It can be `:all`, `:none` or a vector of
    integers. In the first case (the default behavior), all vertical lines will
    be drawn. In the second case, no vertical line will be drawn. In the third
    case, the vertical lines will be drawn only after the columns in the vector.
    Notice that the left border will be drawn if `0` is in `vlines`.
    Furthermore, it is important to mention that the column number in this
    variable is related to the **printed columns**. Thus, it is affected by
    filters, and by the columns added using the variable `show_row_number`.
    Finally, for convenience, the left and right border can be drawn by adding
    the symbols `:begin` and `:end` to this vector, respectively.
    (**Default** = `:none`)
- `wrap_table::Union{Nothing, String}`: This variable controls whether to wrap
    the table in a environment defined by the variable `wrap_table_environment`.
    Defaults to `true`. When `false`, the printed table begins with
    `\\begin{tabular}`. This option does not work with `:longtable`. If it is
    `nothing` then the default option of the table format will be used.
    (**Default** = `nothing`)
- `wrap_table_environment::Union{Nothing, String}`: Environment that will be
    used to wrap the table if the option `wrap_table` is `true`. If it is
    `nothing` then the default option of the table format will be used.
    (**Default** = `nothing`)

## LaTeX highlighters

A set of highlighters can be passed as a `Tuple` to the `highlighters` keyword.
Each highlighter is an instance of the structure [`LatexHighlighter`](@ref). It
contains the following two fields:

- `f::Function`: Function with the signature `f(data, i, j)` in which should
    return `true` if the element `(i, j)` in `data` must be highlighted, or
    `false` otherwise.
- `fd::Functions`: A function with the signature `f(data, i, j, str)::String` in
    which `data` is the matrix, `(i, j)` is the element position in the table,
    and `str` is the data converted to string. This function must return a
    string that will be placed in the cell.

The function `f` has the following signature:

    f(data, i, j)

in which `data` is a reference to the data that is being printed, `i` and `j`
are the element coordinates that are being tested. If this function returns
`true`, then the highlight style will be applied to the `(i, j)` element.
Otherwise, the default style will be used.

If the function `f` returns true, then the function `fd(data, i, j, str)` will
be called and must return the LaTeX string that will be placed in the cell.

There are two helpers that can be used to create LaTeX highlighters:

    LatexHighlighter(f::Function, envs::Union{String,Vector{String}})

    LatexHighlighter(f::Function, fd::Function)

The first will apply recursively all the LaTeX environments in `envs` to the
highlighted text whereas the second let the user select the desired decoration
by specifying the function `fd`.

Thus, for example:

    LatexHighlighter((data, i, j) -> true, ["textbf", "small"])

will wrap all the cells in the table in the following environment:

    \\textbf{\\small{<Cell text>}}

!!! info
    If only a single highlighter is wanted, then it can be passed directly to
    the keyword `highlighter` without being inside a `Tuple`.

!!! note
    If multiple highlighters are valid for the element `(i, j)`, then the
    applied style will be equal to the first match considering the order in the
    tuple `highlighters`.

!!! note
    If the highlighters are used together with [Formatters](@ref), then the
    change in the format **will not** affect the parameter `data` passed to the
    highlighter function `f`. It will always receive the original, unformatted
    value.

---

# Formatters

The keyword `formatters` can be used to pass functions to format the values in
the columns. It must be a tuple of functions in which each function has the
following signature:

    f(v, i, j)

where `v` is the value in the cell, `i` is the row number, and `j` is the column
number. Thus, it must return the formatted value of the cell `(i, j)` that has
the value `v`. Notice that the returned value will be converted to string after
using the function `sprint`.

This keyword can also be a single function, meaning that only one formatter is
available, or `nothing`, meaning that no formatter will be used.

For example, if we want to multiply all values in odd rows of the column 2 by π,
then the formatter should look like:

    formatters = (v, i, j) -> (j == 2 && isodd(i)) ? v * π : v

If multiple formatters are available, then they will be applied in the same
order as they are located in the tuple. Thus, for the following `formatters`:

    formatters = (f1, f2, f3)

each element `v` in the table (i-th row and j-th column) will be formatted by:

    v = f1(v,i,j)
    v = f2(v,i,j)
    v = f3(v,i,j)

Thus, the user must be ensure that the type of `v` between the calls are
compatible.

"""
@inline function pretty_table(data; kwargs...)
    io = stdout isa Base.TTY ? IOContext(stdout, :limit => true) : stdout
    _pretty_table(io, data; kwargs...)
end

pretty_table(io::IO, data; kwargs...) = _pretty_table(io, data; kwargs...)

function pretty_table(::Type{String}, data; kwargs...)
    io = IOBuffer()
    _pretty_table(io, data; kwargs...)
    return String(take!(io))
end

function pretty_table(::Type{HTML}, data; kwargs...)
    # If the keywords does not set the backend or the table format, use the HTML
    # backend by default.
    if !haskey(kwargs, :backend) && !haskey(kwargs, :tf)
        str = pretty_table(String, data; backend = Val(:html), kwargs...)
    else
        str = pretty_table(String, data; kwargs...)
    end

    return HTML(str)
end


text:
using Pkg

Pkg.activate(@__DIR__)

using Artifacts
using FastAI
using Publish

cd(@__DIR__) do
    @eval Publish.Themes.default() = artifact"flux-theme"
end

serve(FastAI, port = 8001)


text:
# Depending on how the internals of `@patch` and `apply` are setup there could be issues
# with reusing a patch for a variety of reasons. These were discovered with experimenting
# with using Cassette with Mocking.

@testset "patch reuse" begin
    @testset "identical signature" begin
        # Declare the original function
        f() = 0

        # Create a patch `p1` and use it
        p1 = @patch f() = 1
        apply(p1) do
            @test (@mock f()) == 1
        end

        # Create another patch `p2` with the same function signature as `p1`
        p2 = @patch f() = 2
        apply(p2) do
            @test (@mock f()) == 2
        end

        # Verify we can re-use the `p1` patch
        apply(p1) do
            @test (@mock f()) == 1
        end
    end

    @testset "optional arguments" begin
        f(args::T...) where T = zero(T)

        # Create a patch function that contains two methods: `f(x)` and `f(x, y)`
        p1 = @patch f(x, y=0) = x + y
        apply(p1) do
            @test (@mock f(5)) == 5
        end

        # Create a second patch which could possibly interfere with `f(x)` from `p1` later
        p2 = @patch f(x) = -x
        apply(p2) do
            @test (@mock f(5)) == -5
        end

        # Validate that executing `f(x)` executes `f(x, 0)`
        apply(p1) do
            @test (@mock f(5)) == 5
        end
    end

    @testset "specificity" begin
        f(::Real) = Real

        # Create our initial patch
        p1 = @patch f(::Number) = Number
        apply(p1) do
            @test (@mock f(1)) == Number
        end

        # Create a second patch which is more specific than the first
        p2 = @patch f(::Integer) = Integer
        apply(p2) do
            @test (@mock f(1)) == Integer
        end

        # Validate that first patch is called and not the more specific second patch
        apply(p1) do
            @test (@mock f(1)) == Number
        end
    end
end