text:
using Distributions, DiffResults, ForwardDiff, DiffResults, AdvancedHMC

# Evaluates target distribution π(q)
function π_q(q)
    q₁ = q[1]
    q₂ = q[2]
    in_bound = (-π/2 < q₁ < π/2) & (-π/2 < q₂ < π/2)
    return in_bound*(2/π * exp(-2*(q₁^2 + q₂^2)) + (sin(q₁*q₂)*sin(q₁)*cos(q₂))^2)
end

# Initial sample
q₀ = [0.0, 0.0]

# State space size
Ω_size = length(q₀)

# Number of samples desired
num_samples = 25000

# Number of ϵ updates
num_adapts = 2000

# Evaluates log of target distribution
function ℓπ(q)
    return log(π_q(q))
end

# Evaluates gradient of log of target distribution
# This implementation of NUTS requires the gradient function to be defined this
# way
function ∂ℓπ_∂q(q)
    result = DiffResults.GradientResult(q)
    ForwardDiff.gradient!(result, ℓπ, q)
    return (DiffResults.value(result), DiffResults.gradient(result))
end

# Use Euclidean matrix. 2 here is dimension of sample space
metric = DiagEuclideanMetric(Ω_size)

# Define Hamiltonian function
hamiltonian = Hamiltonian(metric, ℓπ, ∂ℓπ_∂q)

# Initial step size
ϵ₀ = find_good_eps(hamiltonian, q₀)

# Define leapfrog symplectic integrator
lf_integrator = Leapfrog(ϵ₀)

# Define NUTS sampler
sampler = NUTS{MultinomialTS,GeneralisedNoUTurn}(lf_integrator)

# Define ϵ adapter. Used to find apatively determine optimal step size ϵ
adaptor = StanHMCAdaptor(num_adapts, Preconditioner(metric), NesterovDualAveraging(0.8, lf_integrator))

# Get samples. stats stores information for each sample (e.g., acceptance rate)
samples, stats = sample(hamiltonian, sampler, q₀, num_samples, adaptor, num_adapts; progress=true)

samples_mat = zeros(num_samples, 2)
for i=1:num_samples
    samples_mat[i,1] = samples[i][1]
    samples_mat[i,2] = samples[i][2]
end

using PyCall, Plots, Pandas
@pyimport seaborn as sns
using PyPlot
pygui(true)
data = Pandas.DataFrame(Dict(:q₁=>samples_mat[:,1], :q₂=>samples_mat[:,2]))
rc("axes",facecolor = "#f7ece2")
rc("figure",facecolor = "#f7ece2")
rc("savefig",facecolor = "#f7ece2")
sns.jointplot(x = "q₁", y = "q₂", data = data, s = 1, alpha = 0.18)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
using IRTools
using IRTools: IR, CFG, Variable, isexpr, stmt, argument!, return!, xcall, block!,
  branch!, blocks, insertafter!, arguments, argtypes, isreturn, stackify, isconditional

struct WTuple
  parts::Vector{WType}
end

WTuple(Ts::WType...) = WTuple([Ts...])

function Base.show(io::IO, t::WTuple)
  print(io, "(")
  join(io, t.parts, ", ")
  print(io, ")")
end

# TODO: simpler to treat all variables as potential tuples
# Or better, strip out tuples before register allocation
function locals!(ir::IR)
  locals = WType[]
  ret = WType[]
  env = Dict()
  tuples = Dict()
  rename(x::Variable) = env[x]
  rename(x::Real) = Const(x)
  rename(x::Union{Const,Local}) = x
  ltype(x) = rename(x) isa Const ? rename(x).typ : locals[rename(x).id+1]
  local!(T) = (push!(locals, T); Local(length(locals)-1))
  local!(v, T) = get!(() -> local!(T), env, v)
  local!(v, T::WTuple) = tuples[v] = local!.(T.parts)
  for (arg, T) in zip(arguments(ir), argtypes(ir))
    local!(arg, T)
  end
  for b in blocks(ir)
    for (v, st) in b
      ex = st.expr
      if ex isa Variable
        delete!(ir, v)
        _env = haskey(tuples, ex) ? tuples : env
        _env[v] = _env[ex]
      elseif ex == unreachable
        # leave it alone
      elseif !isexpr(ex)
        delete!(ir, v)
        env[v] = Const(ex)
      elseif isexpr(ex, :call)
        for arg in ex.args[2:end]
          haskey(tuples, arg) ?
            [insert!(ir, v, tuples[arg][i]) for i in 1:length(tuples[arg])] :
            insert!(ir, v, rename(arg))
        end
        ir[v] = ex.args[1]::Instruction
        if st.type isa WTuple
          tuples[v] = [(l = local!(T); insertafter!(ir, v, SetLocal(false, l.id)); l)
                       for T in st.type.parts]
        else
          insertafter!(ir, v, SetLocal(false, local!(v, st.type).id))
        end
      elseif isexpr(ex, :tuple)
        tuples[v] = rename.(ex.args)
        delete!(ir, v)
      elseif isexpr(ex, :ref)
        env[v] = tuples[ex.args[1]][ex.args[2]]
        delete!(ir, v)
      else
        error("Unrecognised wasm expression $ex")
      end
    end
    for br in IRTools.branches(b)
      if isreturn(br)
        args = haskey(tuples, arguments(br)[1]) ? tuples[arguments(br)[1]] : arguments(br)
        ret = ltype.(args)
        for arg in args
          push!(b, rename(arg))
        end
        push!(b, Return())
      else
        for (x, y) in zip(arguments(br), arguments(IRTools.block(ir, br.block)))
          push!(b, rename(x))
          push!(b, SetLocal(false, local!(y, ltype(x)).id))
        end
        isconditional(br) && push!(b, rename(br.condition))
        push!(b, Branch(isconditional(br), br.block))
      end
    end
    empty!(IRTools.branches(b))
  end
  return ir, locals, ret
end

isbackedge((from, to)) = to <= from

function reloop(ir, cfg)
  scopes = []
  targets = []
  forw, back = stackify(cfg)
  push!(scopes, Block([]))
  block!() = (bl = Block([]); push!(scopes[1].body, bl); pushfirst!(scopes, bl))
  loop!() = (bl = Loop([]); push!(scopes[1].body, bl); pushfirst!(scopes, bl))
  for b in blocks(ir)
    any(br -> br[2] == b.id, forw) && (popfirst!(scopes); popfirst!(targets))
    brs = vcat(filter(br -> br[1] == b.id, forw), filter(br -> br[2] == b.id, back))
    sort!(brs, by = br -> isbackedge(br) ? (br[1], 1) : (br[2], -1), rev = true)
    for br in brs
      isbackedge(br) ? loop!() : block!()
      pushfirst!(targets, br[2])
    end
    for (v, st) in b
      if st.expr isa Branch
        st.expr.cond && push!(scopes[1].body, i32.eqz)
        target = findfirst(b -> b == st.expr.level, targets)-1
        push!(scopes[1].body, Branch(st.expr.cond, target))
      else
        push!(scopes[1].body, st.expr)
      end
    end
    any(br -> br[1] == b.id, back) && (popfirst!(scopes); popfirst!(targets))
  end
  return scopes[end]
end

flattentype(Ts) = vcat(flattentype.(Ts)...)
flattentype(T::WType) = [T]
flattentype(T::WTuple) = T.parts

function irfunc(name, ir)
  # @show name
  cfg = CFG(ir)
  ir, locals, ret = locals!(ir)
  params = flattentype(argtypes(ir))
  locals = locals[length(params)+1:end]
  Func(name, params, ret, locals, reloop(ir, cfg))
end

ext: jl


text:
exprs = replace.(readlines("day18/input.txt"), '*' => '⨣')

⨣(a, b) = a * b

println(reduce((v, l) -> v + eval(Meta.parse(l)), exprs; init=0))


text:
module ModuleDataNormalization
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleOperators, MultipleViewGeometry.ModuleMathAliases
using StaticArrays, LinearAlgebra
export hartley_transformation, hartley_normalization!, hartley_normalization
include("hartley_transformation.jl")
end


text:
function numlines(buf)
    s = String(take!(buf))
    sum(!isempty(strip(c)) for c in split(s, '\n'))
end

@testset "printtree length" for t in TEST_TREES
    buf = IOBuffer()
    printtree(buf, t)
    @test numlines(buf) == nnodes(t)
end

@testset "printtree horizontal truncation for infinite trees" begin
    t1 = VectorVertex(1, AbstractVertex[])
    push!(t1.chs, t1)
    t2 = Dict()
    t2["ch"] = t2
    buf1 = IOBuffer()
    printtree(buf1, t1; htrunc=10)
    buf2 = IOBuffer()
    printtree(buf2, t2; htrunc=20)
    @test numlines(buf1) == 10 + 1 # one line for ellipsis ⋮
    @test numlines(buf2) == 20 + 1 # one line for ellipsis ⋮
end

@testset "printtree vertical truncation for many children" for l in 0:20, vtrunc in 0:20
    v = 1:l |> collect
    buf = IOBuffer()
    printtree(buf, v; vtrunc=vtrunc)
    @test numlines(buf) == min(l, vtrunc + 1) + 1 # one line for ellipsis ⋮ and one for the header
end

@testset "printtree labelled children sorted" begin
    t1 = Dict("a" => Leaf(1), "b" => Leaf(2), "c" => Leaf(3))
    t2 = Dict("c" => Leaf(3), "a" => Leaf(1), "b" => Leaf(2))
    t3 = Dict("b" => Leaf(2), "c" => Leaf(3), "a" => Leaf(1))
    @test all([t1, t2, t3]) do t
        buff = IOBuffer()
        printtree(buff, t)
        String(take!(buff)) ==
        """
        Dict of
          ├── a: Leaf (1) comm
          ├── b: Leaf (2) comm
          └── c: Leaf (3) comm
        """
    end
end

@testset "printtree traversal encoding" begin
    t = COMPLETE_BINARY_TREE_1
    buf = IOBuffer()
    printtree(buf, t)
    @test String(take!(buf)) ==
        """
        BinaryVertex (1) comm
          ├── BinaryVertex (2) comm
          │     ├── Leaf (4) comm
          │     └── Leaf (5) comm
          └── BinaryVertex (3) comm
                ├── Leaf (6) comm
                └── Leaf (7) comm
        """
    buf = IOBuffer()
    printtree(buf, t; trav=true)
    @test String(take!(buf)) ==
        """
        BinaryVertex (1) [""] comm
          ├── BinaryVertex (2) ["E"] comm
          │     ├── Leaf (4) ["I"] comm
          │     └── Leaf (5) ["M"] comm
          └── BinaryVertex (3) ["U"] comm
                ├── Leaf (6) ["Y"] comm
                └── Leaf (7) ["c"] comm
        """
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Complex{T<:Real} <: Number

Complex number type with real and imaginary part of type `T`.

`Complex32`, `Complex64` and `Complex128` are aliases for
`Complex{Float16}`, `Complex{Float32}` and `Complex{Float64}` respectively.
"""
struct Complex{T<:Real} <: Number
    re::T
    im::T
end
Complex(x::Real, y::Real) = Complex(promote(x,y)...)
Complex(x::Real) = Complex(x, zero(x))

"""
    im

The imaginary unit.

# Examples
```jldoctest
julia> im * im
-1 + 0im
```
"""
const im = Complex(false, true)

const Complex128 = Complex{Float64}
const Complex64  = Complex{Float32}
const Complex32  = Complex{Float16}

convert(::Type{Complex{T}}, x::Real) where {T<:Real} = Complex{T}(x,0)
convert(::Type{Complex{T}}, z::Complex) where {T<:Real} = Complex{T}(real(z),imag(z))
convert(::Type{T}, z::Complex) where {T<:Real} =
    isreal(z) ? convert(T,real(z)) : throw(InexactError())

convert(::Type{Complex}, z::Complex) = z
convert(::Type{Complex}, x::Real) = Complex(x)

promote_rule(::Type{Complex{T}}, ::Type{S}) where {T<:Real,S<:Real} =
    Complex{promote_type(T,S)}
promote_rule(::Type{Complex{T}}, ::Type{Complex{S}}) where {T<:Real,S<:Real} =
    Complex{promote_type(T,S)}

widen(::Type{Complex{T}}) where {T} = Complex{widen(T)}

float(::Type{Complex{T}}) where {T<:AbstractFloat} = Complex{T}
float(::Type{Complex{T}}) where {T} = Complex{float(T)}

"""
    real(z)

Return the real part of the complex number `z`.

# Examples
```jldoctest
julia> real(1 + 3im)
1
```
"""
real(z::Complex) = z.re

"""
    imag(z)

Return the imaginary part of the complex number `z`.

# Examples
```jldoctest
julia> imag(1 + 3im)
3
```
"""
imag(z::Complex) = z.im
real(x::Real) = x
imag(x::Real) = zero(x)

"""
    reim(z)

Return both the real and imaginary parts of the complex number `z`.

# Examples
```jldoctest
julia> reim(1 + 3im)
(1, 3)
```
"""
reim(z) = (real(z), imag(z))

"""
    real(T::Type)

Returns the type that represents the real part of a value of type `T`.
e.g: for `T == Complex{R}`, returns `R`.
Equivalent to `typeof(real(zero(T)))`.

```jldoctest
julia> real(Complex{Int})
Int64

julia> real(Float64)
Float64
```
"""
real(T::Type) = typeof(real(zero(T)))
real(::Type{T}) where {T<:Real} = T
real(::Type{Complex{T}}) where {T<:Real} = T

"""
    isreal(x) -> Bool

Test whether `x` or all its elements are numerically equal to some real number.

# Examples
```jldoctest
julia> isreal(5.)
true

julia> isreal([4.; complex(0,1)])
false
```
"""
isreal(x::Real) = true
isreal(z::Complex) = iszero(imag(z))
isinteger(z::Complex) = isreal(z) & isinteger(real(z))
isfinite(z::Complex) = isfinite(real(z)) & isfinite(imag(z))
isnan(z::Complex) = isnan(real(z)) | isnan(imag(z))
isinf(z::Complex) = isinf(real(z)) | isinf(imag(z))
iszero(z::Complex) = iszero(real(z)) & iszero(imag(z))

"""
    complex(r, [i])

Convert real numbers or arrays to complex. `i` defaults to zero.
"""
complex(z::Complex) = z
complex(x::Real) = Complex(x)
complex(x::Real, y::Real) = Complex(x, y)

"""
    complex(T::Type)

Returns an appropriate type which can represent a value of type `T` as a complex number.
Equivalent to `typeof(complex(zero(T)))`.

# Examples
```jldoctest
julia> complex(Complex{Int})
Complex{Int64}

julia> complex(Int)
Complex{Int64}
```
"""
complex(::Type{T}) where {T<:Real} = Complex{T}
complex(::Type{Complex{T}}) where {T<:Real} = Complex{T}

flipsign(x::Complex, y::Real) = ifelse(signbit(y), -x, x)

function show(io::IO, z::Complex)
    r, i = reim(z)
    compact = get(io, :compact, false)
    show(io, r)
    if signbit(i) && !isnan(i)
        i = -i
        print(io, compact ? "-" : " - ")
    else
        print(io, compact ? "+" : " + ")
    end
    show(io, i)
    if !(isa(i,Integer) && !isa(i,Bool) || isa(i,AbstractFloat) && isfinite(i))
        print(io, "*")
    end
    print(io, "im")
end
show(io::IO, z::Complex{Bool}) =
    print(io, z == im ? "im" : "Complex($(z.re),$(z.im))")

function read(s::IO, ::Type{Complex{T}}) where T<:Real
    r = read(s,T)
    i = read(s,T)
    Complex{T}(r,i)
end
function write(s::IO, z::Complex)
    write(s,real(z),imag(z))
end

## byte order swaps: real and imaginary part are swapped individually
bswap(z::Complex) = Complex(bswap(real(z)), bswap(imag(z)))

## equality and hashing of complex numbers ##

==(z::Complex, w::Complex) = (real(z) == real(w)) & (imag(z) == imag(w))
==(z::Complex, x::Real) = isreal(z) && real(z) == x
==(x::Real, z::Complex) = isreal(z) && real(z) == x

isequal(z::Complex, w::Complex) = isequal(real(z),real(w)) & isequal(imag(z),imag(w))

if UInt === UInt64
    const h_imag = 0x32a7a07f3e7cd1f9
else
    const h_imag = 0x3e7cd1f9
end
const hash_0_imag = hash(0, h_imag)

function hash(z::Complex, h::UInt)
    # TODO: with default argument specialization, this would be better:
    # hash(real(z), h ⊻ hash(imag(z), h ⊻ h_imag) ⊻ hash(0, h ⊻ h_imag))
    hash(real(z), h ⊻ hash(imag(z), h_imag) ⊻ hash_0_imag)
end

## generic functions of complex numbers ##

"""
    conj(z)

Compute the complex conjugate of a complex number `z`.

# Examples
```jldoctest
julia> conj(1 + 3im)
1 - 3im
```
"""
conj(z::Complex) = Complex(real(z),-imag(z))
abs(z::Complex)  = hypot(real(z), imag(z))
abs2(z::Complex) = real(z)*real(z) + imag(z)*imag(z)
inv(z::Complex)  = conj(z)/abs2(z)
inv(z::Complex{<:Integer}) = inv(float(z))

-(z::Complex) = Complex(-real(z), -imag(z))
+(z::Complex, w::Complex) = Complex(real(z) + real(w), imag(z) + imag(w))
-(z::Complex, w::Complex) = Complex(real(z) - real(w), imag(z) - imag(w))
*(z::Complex, w::Complex) = Complex(real(z) * real(w) - imag(z) * imag(w),
                                    real(z) * imag(w) + imag(z) * real(w))

muladd(z::Complex, w::Complex, x::Complex) =
    Complex(muladd(real(z), real(w), real(x)) - imag(z)*imag(w), # TODO: use mulsub given #15985
            muladd(real(z), imag(w), muladd(imag(z), real(w), imag(x))))

# handle Bool and Complex{Bool}
# avoid type signature ambiguity warnings
+(x::Bool, z::Complex{Bool}) = Complex(x + real(z), imag(z))
+(z::Complex{Bool}, x::Bool) = Complex(real(z) + x, imag(z))
-(x::Bool, z::Complex{Bool}) = Complex(x - real(z), - imag(z))
-(z::Complex{Bool}, x::Bool) = Complex(real(z) - x, imag(z))
*(x::Bool, z::Complex{Bool}) = Complex(x * real(z), x * imag(z))
*(z::Complex{Bool}, x::Bool) = Complex(real(z) * x, imag(z) * x)

+(x::Bool, z::Complex) = Complex(x + real(z), imag(z))
+(z::Complex, x::Bool) = Complex(real(z) + x, imag(z))
-(x::Bool, z::Complex) = Complex(x - real(z), - imag(z))
-(z::Complex, x::Bool) = Complex(real(z) - x, imag(z))
*(x::Bool, z::Complex) = Complex(x * real(z), x * imag(z))
*(z::Complex, x::Bool) = Complex(real(z) * x, imag(z) * x)

+(x::Real, z::Complex{Bool}) = Complex(x + real(z), imag(z))
+(z::Complex{Bool}, x::Real) = Complex(real(z) + x, imag(z))
function -(x::Real, z::Complex{Bool})
    # we don't want the default type for -(Bool)
    re = x-real(z)
    Complex(re, - oftype(re, imag(z)))
end
-(z::Complex{Bool}, x::Real) = Complex(real(z) - x, imag(z))
*(x::Real, z::Complex{Bool}) = Complex(x * real(z), x * imag(z))
*(z::Complex{Bool}, x::Real) = Complex(real(z) * x, imag(z) * x)

# adding or multiplying real & complex is common
+(x::Real, z::Complex) = Complex(x + real(z), imag(z))
+(z::Complex, x::Real) = Complex(x + real(z), imag(z))
function -(x::Real, z::Complex)
    # we don't want the default type for -(Bool)
    re = x - real(z)
    Complex(re, - oftype(re, imag(z)))
end
-(z::Complex, x::Real) = Complex(real(z) - x, imag(z))
*(x::Real, z::Complex) = Complex(x * real(z), x * imag(z))
*(z::Complex, x::Real) = Complex(x * real(z), x * imag(z))

muladd(x::Real, z::Complex, y::Number) = muladd(z, x, y)
muladd(z::Complex, x::Real, y::Real) = Complex(muladd(real(z),x,y), imag(z)*x)
muladd(z::Complex, x::Real, w::Complex) =
    Complex(muladd(real(z),x,real(w)), muladd(imag(z),x,imag(w)))
muladd(x::Real, y::Real, z::Complex) = Complex(muladd(x,y,real(z)), imag(z))
muladd(z::Complex, w::Complex, x::Real) =
    Complex(muladd(real(z), real(w), x) - imag(z)*imag(w), # TODO: use mulsub given #15985
            muladd(real(z), imag(w), imag(z) * real(w)))

/(a::R, z::S) where {R<:Real,S<:Complex} = (T = promote_type(R,S); a*inv(T(z)))
/(z::Complex, x::Real) = Complex(real(z)/x, imag(z)/x)

function /(a::Complex{T}, b::Complex{T}) where T<:Real
    are = real(a); aim = imag(a); bre = real(b); bim = imag(b)
    if abs(bre) <= abs(bim)
        if isinf(bre) && isinf(bim)
            r = sign(bre)/sign(bim)
        else
            r = bre / bim
        end
        den = bim + r*bre
        Complex((are*r + aim)/den, (aim*r - are)/den)
    else
        if isinf(bre) && isinf(bim)
            r = sign(bim)/sign(bre)
        else
            r = bim / bre
        end
        den = bre + r*bim
        Complex((are + aim*r)/den, (aim - are*r)/den)
    end
end

inv(z::Complex{<:Union{Float16,Float32}}) =
    oftype(z, conj(widen(z))/abs2(widen(z)))

/(z::Complex{T}, w::Complex{T}) where {T<:Union{Float16,Float32}} =
    oftype(z, widen(z)*inv(widen(w)))

# robust complex division for double precision
# the first step is to scale variables if appropriate ,then do calculations
# in a way that avoids over/underflow (subfuncs 1 and 2), then undo the scaling.
# scaling variable s and other techniques
# based on arxiv.1210.4539
#             a + i*b
#  p + i*q = ---------
#             c + i*d
function /(z::Complex128, w::Complex128)
    a, b = reim(z); c, d = reim(w)
    half = 0.5
    two = 2.0
    ab = max(abs(a), abs(b))
    cd = max(abs(c), abs(d))
    ov = realmax(a)
    un = realmin(a)
    ϵ = eps(Float64)
    bs = two/(ϵ*ϵ)
    s = 1.0
    ab >= half*ov  && (a=half*a; b=half*b; s=two*s ) # scale down a,b
    cd >= half*ov  && (c=half*c; d=half*d; s=s*half) # scale down c,d
    ab <= un*two/ϵ && (a=a*bs; b=b*bs; s=s/bs      ) # scale up a,b
    cd <= un*two/ϵ && (c=c*bs; d=d*bs; s=s*bs      ) # scale up c,d
    abs(d)<=abs(c) ? ((p,q)=robust_cdiv1(a,b,c,d)  ) : ((p,q)=robust_cdiv1(b,a,d,c); q=-q)
    return Complex128(p*s,q*s) # undo scaling
end
function robust_cdiv1(a::Float64, b::Float64, c::Float64, d::Float64)
    r = d/c
    t = 1.0/(c+d*r)
    p = robust_cdiv2(a,b,c,d,r,t)
    q = robust_cdiv2(b,-a,c,d,r,t)
    return p,q
end
function robust_cdiv2(a::Float64, b::Float64, c::Float64, d::Float64, r::Float64, t::Float64)
    if r != 0
        br = b*r
        return (br != 0 ? (a+br)*t : a*t + (b*t)*r)
    else
        return (a + d*(b/c)) * t
    end
end

function inv(w::Complex128)
    c, d = reim(w)
    half = 0.5
    two = 2.0
    cd = max(abs(c), abs(d))
    ov = realmax(c)
    un = realmin(c)
    ϵ = eps(Float64)
    bs = two/(ϵ*ϵ)
    s = 1.0
    cd >= half*ov  && (c=half*c; d=half*d; s=s*half) # scale down c,d
    cd <= un*two/ϵ && (c=c*bs; d=d*bs; s=s*bs      ) # scale up c,d
    if abs(d)<=abs(c)
        r = d/c
        t = 1.0/(c+d*r)
        p = t
        q = -r * t
    else
        c, d = d, c
        r = d/c
        t = 1.0/(c+d*r)
        p = r * t
        q = -t
    end
    return Complex128(p*s,q*s) # undo scaling
end

function ssqs(x::T, y::T) where T<:AbstractFloat
    k::Int = 0
    ρ = x*x + y*y
    if !isfinite(ρ) && (isinf(x) || isinf(y))
        ρ = convert(T, Inf)
    elseif isinf(ρ) || (ρ==0 && (x!=0 || y!=0)) || ρ<nextfloat(zero(T))/(2*eps(T)^2)
        m::T = max(abs(x), abs(y))
        k = m==0 ? m : exponent(m)
        xk, yk = ldexp(x,-k), ldexp(y,-k)
        ρ = xk*xk + yk*yk
    end
    ρ, k
end

function sqrt(z::Complex{<:AbstractFloat})
    x, y = reim(z)
    if x==y==0
        return Complex(zero(x),y)
    end
    ρ, k::Int = ssqs(x, y)
    if isfinite(x) ρ=ldexp(abs(x),-k)+sqrt(ρ) end
    if isodd(k)
        k = div(k-1,2)
    else
        k = div(k,2)-1
        ρ += ρ
    end
    ρ = ldexp(sqrt(ρ),k) #sqrt((abs(z)+abs(x))/2) without over/underflow
    ξ = ρ
    η = y
    if ρ != 0
        if isfinite(η) η=(η/ρ)/2 end
        if x<0
            ξ = abs(η)
            η = copysign(ρ,y)
        end
    end
    Complex(ξ,η)
end
sqrt(z::Complex) = sqrt(float(z))

# function sqrt(z::Complex)
#     rz = float(real(z))
#     iz = float(imag(z))
#     r = sqrt((hypot(rz,iz)+abs(rz))/2)
#     if r == 0
#         return Complex(zero(iz), iz)
#     end
#     if rz >= 0
#         return Complex(r, iz/r/2)
#     end
#     return Complex(abs(iz)/r/2, copysign(r,iz))
# end

# compute exp(im*theta)
function cis(theta::Real)
    s, c = sincos(theta)
    Complex(c, s)
end

"""
    cis(z)

Return ``\\exp(iz)``.

# Examples
```jldoctest
julia> cis(π) ≈ -1
true
```
"""
function cis(z::Complex)
    v = exp(-imag(z))
    s, c = sincos(real(z))
    Complex(v * c, v * s)
end

"""
    angle(z)

Compute the phase angle in radians of a complex number `z`.

# Examples
```jldoctest
julia> rad2deg(angle(1 + im))
45.0

julia> rad2deg(angle(1 - im))
-45.0

julia> rad2deg(angle(-1 - im))
-135.0
```
"""
angle(z::Complex) = atan2(imag(z), real(z))

function log(z::Complex{T}) where T<:AbstractFloat
    const T1::T  = 1.25
    const T2::T  = 3
    const ln2::T = log(convert(T,2))  #0.6931471805599453
    x, y = reim(z)
    ρ, k = ssqs(x,y)
    ax = abs(x)
    ay = abs(y)
    if ax < ay
        θ, β = ax, ay
    else
        θ, β = ay, ax
    end
    if k==0 && (0.5 < β*β) && (β <= T1 || ρ < T2)
        ρρ = log1p((β-1)*(β+1)+θ*θ)/2
    else
        ρρ = log(ρ)/2 + k*ln2
    end
    Complex(ρρ, angle(z))
end
log(z::Complex) = log(float(z))

# function log(z::Complex)
#     ar = abs(real(z))
#     ai = abs(imag(z))
#     if ar < ai
#         r = ar/ai
#         re = log(ai) + log1p(r*r)/2
#     else
#         if ar == 0
#             re = isnan(ai) ? ai : -inv(ar)
#         elseif isinf(ai)
#             re = oftype(ar,Inf)
#         else
#             r = ai/ar
#             re = log(ar) + log1p(r*r)/2
#         end
#     end
#     Complex(re, angle(z))
# end

function log10(z::Complex)
    a = log(z)
    a/log(oftype(real(a),10))
end
function log2(z::Complex)
    a = log(z)
    a/log(oftype(real(a),2))
end

function exp(z::Complex)
    zr, zi = reim(z)
    if isnan(zr)
        Complex(zr, zi==0 ? zi : zr)
    elseif !isfinite(zi)
        if zr == Inf
            Complex(-zr, oftype(zr,NaN))
        elseif zr == -Inf
            Complex(-zero(zr), copysign(zero(zi), zi))
        else
            Complex(oftype(zr,NaN), oftype(zi,NaN))
        end
    else
        er = exp(zr)
        if iszero(zi)
            Complex(er, zi)
        else
            s, c = sincos(zi)
            Complex(er * c, er * s)
        end
    end
end

function expm1(z::Complex{T}) where T<:Real
    Tf = float(T)
    zr,zi = reim(z)
    if isnan(zr)
        Complex(zr, zi==0 ? zi : zr)
    elseif !isfinite(zi)
        if zr == Inf
            Complex(-zr, oftype(zr,NaN))
        elseif zr == -Inf
            Complex(-one(zr), copysign(zero(zi), zi))
        else
            Complex(oftype(zr,NaN), oftype(zi,NaN))
        end
    else
        erm1 = expm1(zr)
        if zi == 0
            Complex(erm1, zi)
        else
            er = erm1+one(erm1)
            if isfinite(er)
                wr = erm1 - 2 * er * (sin(convert(Tf, 0.5) * zi))^2
                return Complex(wr, er * sin(zi))
            else
                s, c = sincos(zi)
                return Complex(er * c, er * s)
            end
        end
    end
end

function log1p(z::Complex{T}) where T
    zr,zi = reim(z)
    if isfinite(zr)
        isinf(zi) && return log(z)
        # This is based on a well-known trick for log1p of real z,
        # allegedly due to Kahan, only modified to handle real(u) <= 0
        # differently to avoid inaccuracy near z==-2 and for correct branch cut
        u = one(float(T)) + z
        u == 1 ? convert(typeof(u), z) : real(u) <= 0 ? log(u) : log(u)*z/(u-1)
    elseif isnan(zr)
        Complex(zr, zr)
    elseif isfinite(zi)
        Complex(T(Inf), copysign(zr > 0 ? zero(T) : convert(T, pi), zi))
    else
        Complex(T(Inf), T(NaN))
    end
end

function ^(z::Complex{T}, p::Complex{T})::Complex{T} where T<:AbstractFloat
    if p == 2 #square
        zr, zi = reim(z)
        x = (zr-zi)*(zr+zi)
        y = 2*zr*zi
        if isnan(x)
            if isinf(y)
                x = copysign(zero(T),zr)
            elseif isinf(zi)
                x = convert(T,-Inf)
            elseif isinf(zr)
                x = convert(T,Inf)
            end
        elseif isnan(y) && isinf(x)
            y = copysign(zero(T), y)
        end
        Complex(x,y)
    elseif z!=0
        if p!=0 && isinteger(p)
            rp = real(p)
            if rp < 0
                return power_by_squaring(inv(z), convert(Integer, -rp))
            else
                return power_by_squaring(z, convert(Integer, rp))
            end
        end
        exp(p*log(z))
    elseif p!=0 #0^p
        zero(z) #CHECK SIGNS
    else #0^0
        zer = copysign(zero(T),real(p))*copysign(zero(T),imag(z))
        Complex(one(T), zer)
    end
end

function exp2(z::Complex{T}) where T<:AbstractFloat
    er = exp2(real(z))
    theta = imag(z) * log(convert(T, 2))
    s, c = sincos(theta)
    Complex(er * c, er * s)
end
exp2(z::Complex) = exp2(float(z))

function exp10(z::Complex{T}) where T<:AbstractFloat
    er = exp10(real(z))
    theta = imag(z) * log(convert(T, 10))
    s, c = sincos(theta)
    Complex(er * c, er * s)
end
exp10(z::Complex) = exp10(float(z))

function ^(z::T, p::T) where T<:Complex
    if isinteger(p)
        rp = real(p)
        if rp < 0
            return power_by_squaring(inv(float(z)), convert(Integer, -rp))
        else
            return power_by_squaring(float(z), convert(Integer, rp))
        end
    end
    pr, pim = reim(p)
    zr, zi = reim(z)
    r = abs(z)
    rp = r^pr
    theta = atan2(zi, zr)
    ntheta = pr*theta
    if pim != 0 && r != 0
        rp = rp*exp(-pim*theta)
        ntheta = ntheta + pim*log(r)
    end
    sinntheta, cosntheta = sincos(ntheta)
    re, im = rp*cosntheta, rp*sinntheta
    if isinf(rp)
        if isnan(re)
            re = copysign(zero(re), cosntheta)
        end
        if isnan(im)
            im = copysign(zero(im), sinntheta)
        end
    end

    # apply some corrections to force known zeros
    if pim == 0
        if isinteger(pr)
            if zi == 0
                im = copysign(zero(im), im)
            elseif zr == 0
                if isinteger(0.5*pr) # pr is even
                    im = copysign(zero(im), im)
                else
                    re = copysign(zero(re), re)
                end
            end
        else
            dr = pr*2
            if isinteger(dr) && zi == 0
                if zr < 0
                    re = copysign(zero(re), re)
                else
                    im = copysign(zero(im), im)
                end
            end
        end
    end

    Complex(re, im)
end

^(z::Complex, n::Bool) = n ? z : one(z)
^(z::Complex, n::Integer) = z^Complex(n)

^(z::Complex{<:AbstractFloat}, n::Bool) = n ? z : one(z)  # to resolve ambiguity
^(z::Complex{<:Integer}, n::Bool) = n ? z : one(z)        # to resolve ambiguity

^(z::Complex{<:AbstractFloat}, n::Integer) =
    n>=0 ? power_by_squaring(z,n) : power_by_squaring(inv(z),-n)
^(z::Complex{<:Integer}, n::Integer) = power_by_squaring(z,n) # DomainError for n<0

function sin(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if zr == 0
        Complex(F(zr), sinh(zi))
    elseif !isfinite(zr)
        if zi == 0 || isinf(zi)
            Complex(F(NaN), F(zi))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        s, c = sincos(zr)
        Complex(s * cosh(zi), c * sinh(zi))
    end
end

function cos(z::Complex{T}) where T
    F = float(T)
    zr, zi = reim(z)
    if zr == 0
        Complex(cosh(zi), isnan(zi) ? F(zr) : -flipsign(F(zr),zi))
    elseif !isfinite(zr)
        if zi == 0
            Complex(F(NaN), isnan(zr) ? zero(F) : -flipsign(F(zi),zr))
        elseif isinf(zi)
            Complex(F(Inf), F(NaN))
        else
            Complex(F(NaN), F(NaN))
        end
    else
        s, c = sincos(zr)
        Complex(c * cosh(zi), -s * sinh(zi))
    end
end

function tan(z::Complex)
    zr, zi = reim(z)
    w = tanh(Complex(-zi, zr))
    Complex(imag(w), -real(w))
end

function asin(z::Complex)
    zr, zi = reim(z)
    if isinf(zr) && isinf(zi)
        return Complex(copysign(oftype(zr,pi)/4, zr),zi)
    elseif isnan(zi) && isinf(zr)
        return Complex(zi, oftype(zr, Inf))
    end
    ξ = zr == 0       ? zr :
        !isfinite(zr) ? oftype(zr,pi)/2 * sign(zr) :
        atan2(zr, real(sqrt(1-z)*sqrt(1+z)))
    η = asinh(copysign(imag(sqrt(conj(1-z))*sqrt(1+z)), imag(z)))
    Complex(ξ,η)
end

function acos(z::Complex{<:AbstractFloat})
    zr, zi = reim(z)
    if isnan(zr)
        if isinf(zi) return Complex(zr, -zi)
        else         return Complex(zr, zr) end
    elseif isnan(zi)
        if isinf(zr) return Complex(zi, abs(zr))
        elseif zr==0 return Complex(oftype(zr,pi)/2, zi)
        else         return Complex(zi, zi) end
    elseif zr==zi==0
        return Complex(oftype(zr,pi)/2, -zi)
    elseif zr==Inf && zi===0.0
        return Complex(zi, -zr)
    elseif zr==-Inf && zi===-0.0
        return Complex(oftype(zi,pi), -zr)
    end
    ξ = 2*atan2(real(sqrt(1-z)), real(sqrt(1+z)))
    η = asinh(imag(sqrt(conj(1+z))*sqrt(1-z)))
    if isinf(zr) && isinf(zi) ξ -= oftype(η,pi)/4 * sign(zr) end
    Complex(ξ,η)
end
acos(z::Complex) = acos(float(z))

function atan(z::Complex)
    w = atanh(Complex(-imag(z),real(z)))
    Complex(imag(w),-real(w))
end

function sinh(z::Complex)
    zr, zi = reim(z)
    w = sin(Complex(zi, zr))
    Complex(imag(w),real(w))
end

function cosh(z::Complex)
    zr, zi = reim(z)
    cos(Complex(zi,-zr))
end

function tanh(z::Complex{T}) where T<:AbstractFloat
    const Ω = prevfloat(typemax(T))
    ξ, η = reim(z)
    if isnan(ξ) && η==0 return Complex(ξ, η) end
    if 4*abs(ξ) > asinh(Ω) #Overflow?
        Complex(copysign(one(T),ξ),
                copysign(zero(T),η*(isfinite(η) ? sin(2*abs(η)) : one(η))))
    else
        t = tan(η)
        β = 1+t*t #sec(η)^2
        s = sinh(ξ)
        ρ = sqrt(1 + s*s) #cosh(ξ)
        if isinf(t)
            Complex(ρ/s,1/t)
        else
            Complex(β*ρ*s,t)/(1+β*s*s)
        end
    end
end
tanh(z::Complex) = tanh(float(z))

function asinh(z::Complex)
    w = asin(Complex(-imag(z),real(z)))
    Complex(imag(w),-real(w))
end

function acosh(z::Complex)
    zr, zi = reim(z)
    if isnan(zr) || isnan(zi)
        if isinf(zr) || isinf(zi)
            return Complex(oftype(zr, Inf), oftype(zi, NaN))
        else
            return Complex(oftype(zr, NaN), oftype(zi, NaN))
        end
    elseif zr==-Inf && zi===-0.0 #Edge case is wrong - WHY?
        return Complex(oftype(zr,Inf), oftype(zi, -pi))
    end
    ξ = asinh(real(sqrt(conj(z-1))*sqrt(z+1)))
    η = 2atan2(imag(sqrt(z-1)),real(sqrt(z+1)))
    if isinf(zr) && isinf(zi)
        η -= oftype(η,pi)/4 * sign(zi) * sign(zr)
    end
    Complex(ξ, η)
end

function atanh(z::Complex{T}) where T<:AbstractFloat
    const Ω = prevfloat(typemax(T))
    const θ = sqrt(Ω)/4
    const ρ = 1/θ
    x, y = reim(z)
    ax = abs(x)
    ay = abs(y)
    if ax > θ || ay > θ #Prevent overflow
        if isnan(y)
            if isinf(x)
                return Complex(copysign(zero(x),x), y)
            else
                return Complex(real(1/z), y)
            end
        end
        if isinf(y)
            return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
        end
        return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
    elseif ax==1
        if y == 0
            ξ = copysign(oftype(x,Inf),x)
            η = zero(y)
        else
            ym = ay+ρ
            ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym))
            η = copysign(oftype(y,pi)/2 + atan(ym/2), y)/2
        end
    else #Normal case
        ysq = (ay+ρ)^2
        if x == 0
            ξ = x
        else
            ξ = log1p(4x/((1-x)^2 + ysq))/4
        end
        η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2
    end
    Complex(ξ, η)
end
atanh(z::Complex) = atanh(float(z))

function lexcmp(a::Complex, b::Complex)
    c = cmp(real(a), real(b))
    c == 0 || return c
    cmp(imag(a), imag(b))
end

#Rounding complex numbers
#Requires two different RoundingModes for the real and imaginary components
"""
    round(z, RoundingModeReal, RoundingModeImaginary)

Returns the nearest integral value of the same type as the complex-valued `z` to `z`,
breaking ties using the specified [`RoundingMode`](@ref)s. The first
[`RoundingMode`](@ref) is used for rounding the real components while the
second is used for rounding the imaginary components.
"""
function round(z::Complex{<:AbstractFloat}, ::RoundingMode{MR}, ::RoundingMode{MI}) where {MR,MI}
    Complex(round(real(z), RoundingMode{MR}()),
            round(imag(z), RoundingMode{MI}()))
end
round(z::Complex) = Complex(round(real(z)), round(imag(z)))

function round(z::Complex, digits::Integer, base::Integer=10)
    Complex(round(real(z), digits, base),
            round(imag(z), digits, base))
end

float(z::Complex{<:AbstractFloat}) = z
float(z::Complex) = Complex(float(real(z)), float(imag(z)))

big(::Type{Complex{T}}) where {T<:Real} = Complex{big(T)}
big(z::Complex{T}) where {T<:Real} = Complex{big(T)}(z)

## Array operations on complex numbers ##

complex(A::AbstractArray{<:Complex}) = A

function complex(A::AbstractArray{T}) where T
    if !isleaftype(T)
        error("`complex` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(complex(zero(T)))}, A)
end

## promotion to complex ##

_default_type(T::Type{Complex}) = Complex{Int}

max_stars_repo_name: omus/julia

max_issues_repo_name: omus/julia

max_forks_repo_name: omus/julia


text:
using DiffEqBase: DEProblem, remake
using Distributions: Normal, cquantile

dimension(prob::DEProblem) = length(prob.u0)

PhaseRelaxer(prob::LEProblem,
             ::LEProblem,
             ::LESolution,
             ::Terminator,
             integrator_options) =
    PhaseRelaxer(get_integrator(prob.phase_prob; integrator_options...),
                 prob.t_tran)

function step!(relaxer::PhaseRelaxer)
    step!(relaxer.integrator, relaxer.t_tran)
    relaxer.i += 1
end
Base.length(relaxer::PhaseRelaxer) = 1  # FIXME

function get_tangent_prob(prob;
                          x0 = prob.phase_prob.u0,
                          kwargs...)
    u0 = copy(prob.tangent_prob.u0)
    u0[:, 1] = x0
    return remake(prob.tangent_prob; u0=u0, kwargs...)
end

current_state(relaxer::Union{PhaseRelaxer, AbstractRenormalizer}) =
    current_state(relaxer.integrator)

function get_tangent_integrator(prob::LEProblem, relaxer;
                                integrator_options...)
    # Carry on simulated time to tangent integrator [*]:
    tspan = if relaxer.integrator isa ODEIntegrator
        (relaxer.integrator.t, Inf)
    else
        # TODO: Support this in discrete problem
        prob.phase_prob.tspan
    end
    x0 = current_state(relaxer)
    tangent_prob = get_tangent_prob(prob; x0=x0, tspan=tspan)
    return get_integrator(tangent_prob; integrator_options...)
end
# [*] Another way to propagate simulated time to tangent integrator is
# to use set_t!() after init().  However, init() calls phase_prob.f
# with phase_prob.tspan[1] which is less than relaxer.integrator.t.
# It causes problems in some rare cases with stateful problem
# definitions.
# See: [[./clv/clv_problem.jl::get_integrator]]

TangentRenormalizer(relaxer::PhaseRelaxer, prob::LEProblem, sol::LESolution,
                    tmnr::Terminator, integrator_options) =
    TangentRenormalizer(get_tangent_integrator(prob, relaxer;
                                               integrator_options...),
                        prob.t_attr, prob.t_renorm, sol, tmnr)

MLERenormalizer(relaxer::PhaseRelaxer, prob::LEProblem, sol::LESolution,
                tmnr::Terminator, integrator_options) =
    MLERenormalizer(get_tangent_integrator(prob, relaxer;
                                           integrator_options...),
                    prob.t_attr, prob.t_renorm, sol, tmnr)

Base.length(stage::AbstractRenormalizer) =
    ceil(Int, stage.t_attr / stage.t_renorm)

@inline function keepgoing!(stage::AbstractRenormalizer)
    u0 = current_state(stage)
    u0[:, 1] = stage.phase_state
    u0[:, 2:end] = stage.tangent_state
    set_u!(stage.integrator, u0)
    step!(stage.integrator, stage.t_renorm, true)
    assert_success(stage.integrator)
    # TODO: maybe use step!(solver.integrator, dt) aka "inexact step"
    u0 = current_state(stage)
    stage.phase_state[:] = @view u0[:, 1]
    stage.tangent_state[:] = @view u0[:, 2:end]
end

"""
S_n = ((n-1)/n) S_{n-1} + r_n / n
"""
@inline function lyap_add_R!(n, lyap, R)
    for i = 1:length(lyap)
        lyap[i] = lyap_add_r(n, lyap[i], R[i, i])
    end
end

@inline function lyap_add_r(n, lyap, r)
    return (n - 1) / n * lyap + log(r) / n
end

@inline function eye!(A)
    A .= 0
    for i in 1:min(size(A)...)
        @inbounds A[i, i] = 1
    end
    A
end

""" A = A * diag(sgn) """
@inline function A_mul_sign!(A, sgn)
    for i = 1:size(A)[2]
        if !sgn[i]
            A[:, i] *= -1
        end
    end
    A
end

""" A = diag(sgn) * A """
@inline function sign_mul_A!(sgn, A)
    for i = 1:size(A)[1]
        if !sgn[i]
            A[i, :] *= -1
        end
    end
    A
end

""" sgn = sign(diag(A))  """
@inline function sign_diag!(sgn, A)
    for i = 1:size(A)[1]
        sgn[i] = A[i, i] >= 0
    end
    sgn
end

"""
    step!(stage::AbstractRenormalizer)

Evolve the dynamics and then do an orthonormalization.
"""
function step!(stage::AbstractRenormalizer)
    keepgoing!(stage)
    post_evolve!(stage)
    record!(stage, Val{:OS})
    record!(stage, Val{:FTLE})
    stage.is_finished = should_terminate!(stage.tmnr, stage.sol)
end

is_finished(stage::AbstractRenormalizer) =
    stage.is_finished || (stage.i >= length(stage))

function post_evolve!(stage::TangentRenormalizer)
    dim_lyap = length(stage.inst_exponents)
    P = stage.tangent_state
    F = qrfact!(P)
    Q = stage.Q
    A_mul_B!(F[:Q], eye!(Q))  # Q = Matrix(F[:Q])[...]; but lesser allocation
    R = stage.R = F[:R]

    sign_R = stage.sign_R
    sign_diag!(sign_R, R)        # sign_R = diagm(sign(diag(R)))
    A_mul_sign!(Q, sign_R)       # Q = Q * sign_R
    sign_mul_A!(sign_R, R)       # R = signR * R

    stage.i += 1
    stage.tangent_state, stage.Q = Q, stage.tangent_state
    # At this point:
    # stage.tangent_state === Q == 𝑮ₙ₊ₖ
    # stage.Q is a mess, as qrfact! use it as a "buffer"
end

function post_evolve!(stage::MLERenormalizer)
    v = stage.tangent_state[:, 1]
    r = norm(v)
    n = (stage.i += 1)
    stage.exponent = lyap_add_r(n, stage.exponent, r)
    stage.inst_exponent = log(r) / t_chunk(stage)
    # TODO: don't re-calculate log(r)
    stage.tangent_state[:, 1] .= v ./ r
end

t_chunk(stage) = stage.t_renorm

function record!(stage::TangentRenormalizer{<: LESolRecOS}, ::Type{Val{:OS}})
    dim_lyap = length(stage.inst_exponents)
    R = stage.R
    dt = t_chunk(stage)
    @assert size(R) == (dim_lyap, dim_lyap)
    for i in 1:dim_lyap
        @inbounds stage.inst_exponents[i] = log(R[i, i]) / dt
    end
    OnlineStats.fit!(stage.sol.series, stage.inst_exponents)
end

function record!(stage::MLERenormalizer{<: LESolRecOS}, ::Type{Val{:OS}})
    # FIXME: Assuming stage.sol.main is VecMean and alike
    OnlineStats.fit!(stage.sol.series, [stage.inst_exponent])
end

function record!(stage::AbstractRenormalizer{<: LESolRecFTLE},
                 ::Type{Val{:FTLE}})
    stage.sol.ftle_history[stage.i] .= ftle(stage)
    stage.sol.num_orth = stage.i
end

current_result(stage::TangentRenormalizer) = stage.inst_exponents
current_result(stage::MLERenormalizer) = stage.inst_exponent

"""
    lyapunov_exponents(solver)

Get the result of Lyapunov exponents calculation stored in `solver`.
"""
lyapunov_exponents(solver::Union{LESolverRecOS,
                                 AbstractRenormalizer{<: LESolRecOS}}) =
    lyapunov_exponents(solver.sol)

lyapunov_exponents(sol::LESolRecOS) = mean(sol.main_stat)

@inline function lyapunov_exponents(stage::MLERenormalizer)
    [stage.exponent / t_chunk(stage)]
end

"""Get finite-time Lyapunov exponents (FTLE)"""
ftle(stage::TangentRenormalizer) = stage.inst_exponents
ftle(stage::MLERenormalizer) = [stage.inst_exponent]

ftle_history(solver::Union{LESolverRecFTLE,
                           AbstractRenormalizer{<: LESolRecFTLE}},
             args...) =
    ftle_history(solver.sol, args...)

ftle_history(sol::LESolRecFTLE) = @view sol.ftle_history[1:sol.num_orth]
ftle_history(sol::LESolRecFTLE, i) =
    [x[i] for x in sol.ftle_history[1:sol.num_orth]]

function ftle_history(::Type{A}, sol::LESolRecFTLE) where {A <: AbstractMatrix}
    dim_lyap = length(sol.ftle_history[1])
    num_orth = sol.num_orth
    arr = A(dim_lyap, num_orth)
    for (i, ftle) in enumerate(ftle_history(sol))
        arr[:, i] .= ftle
    end
    return arr
end

const Vecs = AbstractVector{<: AbstractVector}

exponents_history(sol) = exponents_history(ftle_history(sol))

function exponents_history(ftle_history::Vecs)
    t_attr = length(ftle_history)
    dim_lyap = length(ftle_history[1])

    m = VecMean(dim_lyap)
    s = OnlineStats.Series(m)
    le_hist = similar(ftle_history[1], (dim_lyap, t_attr))
    for i in 1:t_attr
        OnlineStats.fit!(s, ftle_history[i])
        le_hist[:, i] .= mean(m)
    end
    return le_hist
end

exponents_stat_history(sol::Union{LESolRecFTLE,
                                  AbstractRenormalizer{<: LESolRecFTLE},
                                  LESolverRecFTLE}) =
    exponents_stat_history(ftle_history(sol))

function exponents_stat_history(ftle_history::Vecs, coverageprob = 0.95)
    t_attr = length(ftle_history)
    dim_lyap = length(ftle_history[1])
    c = cquantile(Normal(), (1 - coverageprob) / 2)
    # Usually t_attr is big so using TDist does not change much here.

    v = VecVariance(dim_lyap)
    s = OnlineStats.Series(v)
    le_hist = similar(ftle_history[1], (dim_lyap, t_attr))
    ci_hist = similar(le_hist)

    OnlineStats.fit!(s, ftle_history[1])
    le_hist[:, 1] .= mean(v)
    ci_hist[:, 1] .= NaN
    for i in 2:t_attr
        OnlineStats.fit!(s, ftle_history[i])
        # c = cquantile(TDist(i - 1), (1 - coverageprob) / 2)
        le_hist[:, i] .= mean(v)
        ci_hist[:, i] .= std(v) / sqrt(i) * c
    end
    return le_hist, ci_hist
end
# MAYBE: use OnlineStats.Bootstrap
# http://joshday.github.io/OnlineStats.jl/latest/api.html#OnlineStats.Bootstrap


text:
function priors_twogroup_newvariant_cities(θ)
	@unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    LP = 0.
    LP += logpdf(Beta(8,2),ct_min1)
    LP += logpdf(Gamma(3,2.5/3),R₀)
    LP += logpdf(Beta(45,5),ϵ)
    LP += logpdf(Gamma(10,4.5/10),χ₁)
    LP += logpdf(Gamma(10,1.5/10),χ₂)
    LP += logpdf(Gamma(3,1e-4/3),p_test₁)#In the region of 1% detection rate
    LP += logpdf(Gamma(3,5e-4/3),p_test₂)
    LP += logpdf(Beta(40,60),p_choose1)
	LP += logpdf(Beta(35,15),P_eff)
    LP += logpdf(Gamma(100,1.5/100),extra_transmissibility)
    LP += logpdf(Gamma(5,100/5),influx_exposed_new_variant)
	LP += logpdf(Normal(0,0.1),log(χ_boost))
	LP += logpdf(Normal(0,0.1),log(p_test_boost))
	LP += logpdf(Gamma(3,100/3),E₀)

    return LP
end

function priors_twogroup_newvariant_semi_urbanrural(θ)
	@unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    LP = 0.
    LP += logpdf(Beta(8,2),ct_min1)
    LP += logpdf(Gamma(3,2/3),R₀)
    LP += logpdf(Beta(45,5),ϵ)
    LP += logpdf(Gamma(5,4.5/5),χ₁)
    LP += logpdf(Gamma(10,1.5/10),χ₂)
    LP += logpdf(Gamma(3,2e-5/3),p_test₁)
    LP += logpdf(Gamma(3,1e-4/3),p_test₂)
    LP += logpdf(Beta(80,20),p_choose1)
	LP += logpdf(Beta(95,5),P_eff)
    LP += logpdf(Gamma(100,1.5/100),extra_transmissibility)
    LP += logpdf(Gamma(5,100/5),influx_exposed_new_variant)
	LP += logpdf(Normal(0,0.1),log(χ_boost))
	LP += logpdf(Normal(0,0.1),log(p_test_boost))
	LP += logpdf(Gamma(3,1/3),E₀)

    return LP
end

function priors_twogroup_newvariant_rural(θ)
	@unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    LP = 0.
    LP += logpdf(Beta(8,2),ct_min1)
    LP += logpdf(Gamma(3,1.5/3),R₀)
    LP += logpdf(Beta(45,5),ϵ)
    LP += logpdf(Gamma(3,4.5/3),χ₁)
    LP += logpdf(Gamma(10,1.5/10),χ₂)
    LP += logpdf(Gamma(3,1e-5/3),p_test₁)
    LP += logpdf(Gamma(3,5e-5/3),p_test₂)
    LP += logpdf(Beta(90,10),p_choose1)
	LP += logpdf(Beta(95,5),P_eff)
    LP += logpdf(Gamma(100,1.5/100),extra_transmissibility)
    LP += logpdf(Gamma(5,100/5),influx_exposed_new_variant)
	LP += logpdf(Normal(0,0.1),log(χ_boost))
	LP += logpdf(Normal(0,0.1),log(p_test_boost))
	LP += logpdf(Gamma(3,0.1/3),E₀)

    return LP
end

function priors_twogroup_newvariant_semi_urbanrural_fitted_detection(θ)
	@unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    LP = 0.
    LP += logpdf(Beta(8,2),ct_min1)
    LP += logpdf(Gamma(3,2/3),R₀)
    LP += logpdf(Beta(45,5),ϵ)
    LP += logpdf(semi_urban_rural_fits.d_χ1,χ₁)
    LP += logpdf(semi_urban_rural_fits.d_χ2,χ₂)
    LP += logpdf(semi_urban_rural_fits.d_ptest1,p_test₁)
    LP += logpdf(semi_urban_rural_fits.d_ptest2,p_test₂)
    LP += logpdf(Beta(80,20),p_choose1)
	LP += logpdf(Beta(95,5),P_eff)
    LP += logpdf(Gamma(100,1.5/100),extra_transmissibility)
    LP += logpdf(Gamma(5,100/5),influx_exposed_new_variant)
	LP += logpdf(Normal(0,0.1),log(χ_boost))
	LP += logpdf(Normal(0,0.1),log(p_test_boost))
	LP += logpdf(Gamma(3,1/3),E₀)

    return LP
end

function priors_twogroup_newvariant_rural_fitted_detection(θ)
	@unpack ct_min1,R₀,ϵ,χ₁,χ₂,p_test₁,p_test₂,p_choose1,P_eff,schooleffect,extra_transmissibility,influx_exposed_new_variant,p_test_boost,χ_boost,E₀ = θ

    LP = 0.
    LP += logpdf(Beta(8,2),ct_min1)
    LP += logpdf(Gamma(3,1.5/3),R₀)
    LP += logpdf(Beta(45,5),ϵ)
    LP += logpdf(rural_fits.d_χ1,χ₁)
    LP += logpdf(rural_fits.d_χ2,χ₂)
    LP += logpdf(rural_fits.d_ptest1,p_test₁)
    LP += logpdf(rural_fits.d_ptest2,p_test₂)
    LP += logpdf(Beta(90,10),p_choose1)
	LP += logpdf(Beta(95,5),P_eff)
    LP += logpdf(Gamma(100,1.5/100),extra_transmissibility)
    LP += logpdf(Gamma(5,100/5),influx_exposed_new_variant)
	LP += logpdf(Normal(0,0.1),log(χ_boost))
	LP += logpdf(Normal(0,0.1),log(p_test_boost))
	LP += logpdf(Gamma(3,0.1/3),E₀)

    return LP
end


text:
using Test
using LightGraphs

include("../../src/Generate/filteredges.jl")
include("../../src/Generate/addedges.jl")

@testset "addedges" begin
    @testset "check_options" begin
        options = Dict(:degree_max => Dict("degree" => 3),
                       :degree_min => Dict("degree" => 4))
        @test_throws ErrorException(":degree_max cannot be strictly inferior to :degree_min") check_options_src(options)
        options = Dict(:size_biggest_clique_max => Dict("size" => 3),
                       :size_biggest_clique_min => Dict("size" => 4))
        @test_throws ErrorException(":size_biggest_clique_max cannot be strictly inferior to :size_biggest_clique_min") check_options_src(options)
        options = Dict(:size_smallest_clique_max => Dict("size" => 3),
                       :size_smallest_clique_min => Dict("size" => 4))
        @test_throws ErrorException(":size_smallest_clique_max cannot be strictly inferior to :size_smallest_clique_min") check_options_src(options)

        options = Dict(:distance_max => Dict("dist" => 3),
                       :distance_min => Dict("dist" => 4))
        @test_throws ErrorException(":distance_max cannot be strictly inferior to :distance_min") check_options_dst(options)
        options = Dict(:same_clique => nothing,
                       :no_same_clique => nothing)
        @test_throws ErrorException(":same_clique and :no_same_clique cannot be set together") check_options_dst(options)
        options = Dict(:has_edge => nothing,
                       :no_edge => nothing)
        @test_throws ErrorException(":has_edge and :no_edge cannot be set together") check_options_dst(options)
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
# using Dates

@testset "Sid.jl" begin
    filename = "data\\test.sid"
    dc = DDR2import.Sid.read(filename)
    @test dc["EEEI"][1] == "BALTI"
    @test dc["EEEI"][end] == "TLL"
    @test dc["BIBL"][1] == "RH"
end


text:
module GrupoBeta
include("Algoritmo_QR.jl")
include("Metodo_potencias.jl")

export autovalores
export autovaloresRef
export autovetores
export autovetoresRef
export iteracoes_QR
export tudo_QR
export tudo_QR_Ref
export autovalor_dominante
export autovetor_dominante
export iteracoes_potencias
export tudo_potencias

# A função deverá ser chamada no terminal com a matriz que se
# deseja os respectivos resultados

# Chamada de exemplo: autovalores([1 1 1; 1 2 1; 1 1 3])

function autovalores(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)
  return M
end

function autovaloresRef(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)
  m, n = size(M)
  for i = 1:m
    for j = 1:n
      if M[i,j] < 10e-10
        M[i,j] = 0.0
      end
    end
    println(M[i,:])
  end
end

function autovetores(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)
  return X
end

function autovetoresRef(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)
  m, n = size(X)
  for i = 1:m
    for j = 1:n
      if X[i,j] < 10e-10
        X[i,j] = 0.0
      end
    end
    println(X[i,:])
  end
end

function iteracoes_QR(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)
  return iter
end

function tudo_QR(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)
  return E, M, X, iter
end

function tudo_QR_Ref(E::Matrix)
  E, M, X, iter = algoritmo_QR(E)

  mE, nE = size(E)
  for i = 1:mE
    println(E[i,:])
  end

  println("")
  mM, nM = size(M)
  for i = 1:mM
    for j = 1:nM
      if M[i,j] < 10e-10
        M[i,j] = 0.0
      end
    end
    println(M[i,:])
  end

  println("")
  mX, nX = size(X)
  for i = 1:mX
    for j = 1:nX
      if X[i,j] < 10e-10
        X[i,j] = 0.0
      end
    end
    println(X[i,:])
  end

  println("")
  println(iter)
end

function autovalor_dominante(E::Matrix)
  λ, y0, iter = metodo_potencias(E)
  return λ
end

function autovetor_dominante(E::Matrix)
  λ, y0, iter = metodo_potencias(E)
  return y0
end

function iteracoes_potencias(E::Matrix)
  λ, y0, iter = metodo_potencias(E)
  return iter
end

function tudo_potencias(E::Matrix)
  return metodo_potencias(E)
end
end

max_stars_repo_name: EleonoraAvello/GrupoBetaPacote

max_issues_repo_name: EleonoraAvello/GrupoBetaPacote

max_forks_repo_name: EleonoraAvello/GrupoBetaPacote


text:
# module MountainCar

using Random

module MountainCarConst
const vel_limit = (-0.07, 0.07)
const pos_limit = (-1.2, 0.5)
const pos_initial_range = (-0.6, 0.4)
end

"""
MountainCar(rng)

# Arguments:
`rng`: The random number generator.
"""
mutable struct MountainCar <: AbstractEnvironment
    pos::Float64
    vel::Float64
    actions::AbstractSet
    normalized::Bool
    function MountainCar(pos, vel, normalized::Bool=false)
        mcc = MountainCarConst
        @boundscheck (pos >= mcc.pos_limit[1] && pos <= mcc.pos_limit[2])
        @boundscheck (vel >= mcc.vel_limit[1] && vel <= mcc.vel_limit[2])
        new(pos, vel, Set(1:3), normalized)
    end
end

MountainCar(normalized::Bool=false) = MountainCar(0.0, 0.0, normalized)
MountainCar(rng::AbstractRNG, normalized::Bool=false) =
    MountainCar((rand(rng)*(MountainCarConst.pos_initial_range[2]
                            - MountainCarConst.pos_initial_range[1])
                 + MountainCarConst.pos_initial_range[1]),
                0.0,
                normalized)

function reset!(env::MountainCar; rng = Random.GLOBAL_RNG, kwargs...)
    # throw("Implement reset! for environment $(typeof(env))")
    env.pos = (rand(rng)*(MountainCarConst.pos_initial_range[2]
                          - MountainCarConst.pos_initial_range[1])
               + MountainCarConst.pos_initial_range[1])
    env.vel = 0
end

get_actions(env::MountainCar) = env.actions
valid_action(env::MountainCar, action) = action in env.actions

function environment_step!(env::MountainCar, action; rng=Random.GLOBAL_RNG, kwargs...)
    @boundscheck valid_action(env, action)
    env.vel = clamp(env.vel + (action - 2)*0.001 - 0.0025*cos(3*env.pos), MountainCarConst.vel_limit...)
    env.pos = clamp(env.pos + env.vel, MountainCarConst.pos_limit...)
end

function get_reward(env::MountainCar) # -> determines if the agent_state is terminal
    if env.pos >= MountainCarConst.pos_limit[2]
        return 0
    end
    return -1
end

function is_terminal(env::MountainCar) # -> determines if the agent_state is terminal
    return env.pos >= MountainCarConst.pos_limit[2]
end

function get_state(env::MountainCar)
    if env.normalized
        return get_normalized_state(env)
    else
        return [env.pos, env.vel]
    end
end

function get_normalized_state(env::MountainCar)
    pos_limit = MountainCarConst.pos_limit
    vel_limit = MountainCarConst.vel_limit
    return [(env.pos - pos_limit[1])/(pos_limit[2] - pos_limit[1]), (env.vel - vel_limit[1])/(vel_limit[2] - vel_limit[1])]
end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
export defaultfont

Txt(pos::Point,
    text::String,
    size::Real,
    halign::Symbol,
    valign::Symbol=:c,
    angle::Angle=deg(0),
    font::String=defaultfont()) = Txt(pos, text, size, halign, valign, angle, font, nothing)

Txt(t::Txt, extent::TextExtent) = Txt(t.pos, t.text, t.size, t.halign, t.valign, t.angle, t.font, extent)

function Txt(pos::Point, text::String, size::Real; ha=:c, va=:c, angle=deg(0), font=defaultfont())
    Txt(pos, text, size, ha, va, angle, font)
end

needed_attributes(::Type{Txt}) = (Visible, Textfill)

function Base.:*(r::Real, te::TextExtent)
    TextExtent(
        te.xbearing * r,
        te.ybearing * r,
        te.width * r,
        te.height * r,
        te.xadvance * r,
        te.yadvance * r
    )
end

_defaultfont = "Helvetica Neue"

defaultfont() = _defaultfont
function defaultfont(font::String)
    global _defaultfont = font
end


text:
const number_of_chebyshevs_to_compile_into_binaries = 20

function next_chebyshev(previous::Sum_Of_Functions, two_previous::Sum_Of_Functions)
    return PE_Function(2.0,0.0,0.0,1) * previous - two_previous
end

first_kind_chebyshevs = Array{Sum_Of_Functions}(undef, number_of_chebyshevs_to_compile_into_binaries)
first_kind_chebyshevs[1] = Sum_Of_Functions([PE_Function(1.0,0.0,0.0,0)])
first_kind_chebyshevs[2] = Sum_Of_Functions([PE_Function(1.0,0.0,0.0,1)])
for i in 3:number_of_chebyshevs_to_compile_into_binaries
    first_kind_chebyshevs[i] = next_chebyshev(first_kind_chebyshevs[i-1], first_kind_chebyshevs[i-2])
end

second_kind_chebyshevs = Array{Sum_Of_Functions}(undef, number_of_chebyshevs_to_compile_into_binaries)
second_kind_chebyshevs[1] = Sum_Of_Functions([PE_Function(1.0,0.0,0.0,0)])
second_kind_chebyshevs[2] = Sum_Of_Functions([PE_Function(2.0,0.0,0.0,1)])
for i in 3:number_of_chebyshevs_to_compile_into_binaries
    second_kind_chebyshevs[i] = next_chebyshev(second_kind_chebyshevs[i-1], second_kind_chebyshevs[i-2])
end

"""
    get_chevyshevs_up_to(N::Integer, first_kind::Bool = true)

Get the first N chebyshev polynomials returned as a vector of `UnivariateFunction`s.
The first 20 polynomials of each are precompiled into the binaries for speed. If
you need more than that they will be calculated at runtime.

These can be from either the first kind or second kind polynomial sequence.
### Inputs
* `N` - How many chebyshev polynomials do you want.
* `first_kind` - A Bool. If true you get first kind polynomials. If false you get second kind.
### Returns
* A `Vector` of `UnivariateFunction`s for each polynomial.
"""
function get_chevyshevs_up_to(N::Integer, first_kind::Bool = true)
    chebyshevs = Array{Sum_Of_Functions}(undef, N)
    if N >= number_of_chebyshevs_to_compile_into_binaries
        if first_kind
            chebyshevs[1:number_of_chebyshevs_to_compile_into_binaries] = first_kind_chebyshevs
        else
            chebyshevs[1:number_of_chebyshevs_to_compile_into_binaries] = second_kind_chebyshevs
        end
        for i in (number_of_chebyshevs_to_compile_into_binaries+1):N
            chebyshevs[i] = next_chebyshev(chebyshevs[i-1], chebyshevs[i-2])
        end
    else
        if first_kind
            chebyshevs[1:N] = first_kind_chebyshevs[1:N]
        else
            chebyshevs[1:N] = second_kind_chebyshevs[1:N]
        end
    end
    return chebyshevs
end


text:
srand(0)
# Column-stochastic x
x_cs = rand(10, 10, 3)
x_cs = x_cs./sum(x_cs, 1)
# Row stochastic x
x_rs = similar(x_cs);
for k = 1:3
  x_rs[:, :, k] = x_cs[:, :, k]';
end

r = rand(10)
q = Array(Float64, 10, 3)

v = rand(10)
my_v = copy(v)

function normal_max(q, v, x, r)
  for k = 1:3
    q[:, k] = r + x[:, :, k] * v
  end
  [maximum(q[r, :]) for r = 1:10]
end

v = normal_max(q, v, x_rs, r)

function running_max(newv, oldv, x, r, q)
  newv[:] = copy(oldv)
  for k = 1:3
    for j = 1:10
      myq = r[j] + dot(x[:, j, k], oldv)
      #println( abs(myq - q[j, k]))
      if newv[j] < myq
        newv[j] = myq
      end
    end
  end
end

new_v = similar(my_v)
running_max(new_v, my_v, x_cs, r, q)

# myv and v should be equal
abs(v .- new_v) .< eps(max(maximum(abs(v)), maximum(abs(myv))))

# Indmax

x = rand(10, 3)

v1 = [maximum(x[r, :]) for r in 1:10]
p1 = [indmax(x[r, :]) for r in 1:10]

p2 = ones(Int, 10)
v2 = zeros(10)
for j = 1:3
  for i = 1:10
    if v2[i] < x[i, j]
      v2[i] = x[i, j]
      p2[i] = j
    end
  end
end


text:
using BivariateCopulas
using StatsBase
using Test

include("archimedean/clayton.jl")


text:
@doc raw"""
    IsoRectangle2

An object `r` of the data type [`IsoRectangle2`](@ref) is a rectangle in the
Euclidean plane ``Ε^2`` with sides parallel to the ``x`` and ``y`` axis
of the coordinate system.

Although they are represented in a canonical form by only two vertices, namely
the lower left and the upper right vertex, we provide functions for "accessing"
the other vertices as well. The vertices are returned in counterclockwise order.

Iso-oriented rectangles and bounding boxes are quite similar. The difference
however is that bounding boxes have always double coordinates, whereas the
coordinate type of an iso-oriented rectangle is chosen by the user.
"""
IsoRectangle2

"""
    IsoRectangle2(p::Point2, q::Point2)

Introduces an iso-oriented rectangle `r` with diagonal opposite vertices `p` and
`q`.

Note that the object is brought in the canonical form.
"""
IsoRectangle2(p::Point2, q::Point2)

"""
    IsoRectangle2(p::Point2, q::Point2, ::Integer)

Introduces an iso-oriented rectangle `r` with diagonal opposite vertices `p` and
`q`.

Note that the object is brought in the canonical form.
"""
IsoRectangle2(p::Point2, q::Point2, ::Integer)

@doc raw"""
    IsoRectangle2(left::Point2, right::Point2, bottom::Point2, top::Point2)

Introduces an iso-oriented rectangle `r` whose minimal ``x`` coordinate is the
one of `left`, the maximal ``x`` coordinate is the one of `right`, the minimal
``y`` coordinate is the one of `bottom`, the maximal ``y`` coordinate is the one
of ``top``.
"""
IsoRectangle2(left::Point2, right::Point2, bottom::Point2, top::Point2)

"""
    IsoRectangle2(min_hx::Real, min_hy::Real,
                  max_hx::Real, max_hy::Real, hw::Real = 1)

Introduces an iso-oriented rectangle `r` with diagonal opposite vertices
`(min_hx/hw, min_hy/hw)` and `(max_hx/hw, max_hy/hw)`.

!!! info "Precondition"

    `hw` ``≠`` `0`.
"""
IsoRectangle2(min_hx::Real, min_hy::Real,
              max_hx::Real, max_hy::Real, hw::Real = 1) =
    IsoRectangle2(convert.(FT, (min_hx, min_hy, max_hx, max_hy, hw))...)

"""
    IsoRectangle2(bbox::Bbox2)

If [`RT`](@ref) is constructible from `Float64`, introduces an iso-oriented
rectangle from `bbox`.
"""
IsoRectangle2(bbox::Bbox2)

"""
    ==(r::IsoRectangle2, r2::IsoRectangle2)

Test for equality: two iso-oriented rectangles are equal, iff their lower left
and upper right vertices are equal.
"""
==(r::IsoRectangle2, r2::IsoRectangle2)

"""
    vertex(r::IsoRectangle2, i::Integer)

Returns the `i`ᵗʰ vertex modulo 4 of `r` in counterclockwise order, starting
with the lower left vertex.
"""
vertex(r::IsoRectangle2, i::Integer)

"""
    min(r::IsoRectangle2)

Returns the lower left vertex of `r` (= `vertex(r, 0)`).
"""
min(r::IsoRectangle2)

"""
    max(r::IsoRectangle2)

Returns the upper right vertex of `r` (= `vertex(r, 2)`).
"""
max(r::IsoRectangle2)

"""
    xmin(r::IsoRectangle2)

Returns the ``x`` coordinate of the lower left vertex of `r`.
"""
xmin(r::IsoRectangle2)

"""
    ymin(r::IsoRectangle2)

Returns the ``y`` coordinate of the lower left vertex of `r`.
"""
ymin(r::IsoRectangle2)

"""
    xmax(r::IsoRectangle2)

Returns the ``x`` coordinate of the upper right vertex of `r`.
"""
xmax(r::IsoRectangle2)

"""
    ymax(r::IsoRectangle2)

Returns the ``y`` coordinate of the upper right vertex of `r`.
"""
ymax(r::IsoRectangle2)

"""
    min_coord(r::IsoRectangle2, i::Integer)

Returns the `i`ᵗʰ Cartesian coordinate of the lower left vertex of `r`.

!!! info "Precondition"

    ``0 ≤ i ≤ 1``
"""
min_coord(r::IsoRectangle2, i::Integer)

"""
    max_coord(r::IsoRectangle2, i::Integer)

Returns the `i`ᵗʰ Cartesian coordinate of the upper right vertex of `r`.

!!! info "Precondition"

    ``0 ≤ i ≤ 1``
"""
max_coord(r::IsoRectangle2, i::Integer)

"""
    is_degenerate(r::IsoRectangle2)

`r` is degenerate, if all vertices are collinear.
"""
is_degenerate(r::IsoRectangle2)

"""
    bounded_side(r::IsoRectangle2, p::Point2)

Returns either [`ON_UNBOUNDED_SIDE`](@ref), [`ON_BOUNDED_SIDE`](@ref), or the
constant [`ON_BOUNDARY`](@ref), depending on where point `p` is.
"""
bounded_side(r::IsoRectangle2, p::Point2)

"""
    has_on_boundary(r::IsoRectangle2, p::Point2)

Returns `true`, iff `p` lies on the boundary of `r`.
"""
has_on_boundary(c::IsoRectangle2, p::Point2)

"""
    has_on_bounded_side(r::IsoRectangle2, p::Point2)

Returns `true`, iff `p` lies properly inside of `r`.
"""
has_on_bounded_side(c::IsoRectangle2, p::Point2)

"""
    has_on_unbounded_side(r::IsoRectangle2, p::Point2)

Returns `true`, iff `p` lies properly outside of `r`.
"""
has_on_unbounded_side(c::IsoRectangle2, p::Point2)

"""
    area(r::IsoRectangle2)

Returns the area of `r`.
"""
area(r::IsoRectangle2)

"""
    bbox(r::IsoRectangle2)

Returns a bounding box containing `r`.
"""
bbox(r::IsoRectangle2)

@doc raw"""
    transform(r::IsoRectangle2, t::AffTransformation2)

Returns the iso-oriented rectangle obtained by applying `t` on the lower left
and upper right corner of `r`.

!!! info "Precondition"

    The angle at a rotation must be a multiple of ``π \over 2``, otherwise the
    resulting rectangle does not have the same side length. Note that rotating
    about an arbitrary angle can even result in a degenerate iso-oriented
    rectangle.
"""
transform(r::IsoRectangle2, t::AffTransformation2)


text:
module MPIExtensions

using MPI
import MPI: Buffer, Comm, MPI_Comm, MPIPtr, MPI_Request, Request, RBuffer, libmpi, Op, MPI_Op, MPI_Datatype, @mpichk,
            free, Allreduce!, IN_PLACE
import Base: wait
import MPI: Wait!, Waitall!
import Flux: cpu, gpu
import CUDA: CuArray

const mpi_is_cuda_aware = Ref(false)

function __init__()
    # If using `mpi_is_cuda_aware` anywhere use the @static macro
    # since we anyways need to recompile the code when MPI
    # implementation changes
    mpi_is_cuda_aware[] = MPI.has_cuda()
    if !mpi_is_cuda_aware[]
        @warn "MPI Implementation is not CUDA Aware"
    end
end

"""
    Allreduce!(v, op, comm)

Perform `MPI.Allreduce!` ensuring that CuArrays are safely transfered to
CPU if CUDA-aware MPI is unavailable.
"""
function Allreduce!(v::CuArray, op::Any, comm::MPI_Comm)
    @static if !mpi_is_cuda_aware[]
        v = cpu(v)
    end
    MPI.Allreduce!(v, op, comm)
    @static if !mpi_is_cuda_aware[]
        v = gpu(v)
    end
    return v
end

"""
    Bcast!(v, op, comm)

Perform `MPI.Bcast!` ensuring that CuArrays are safely transfered to
CPU if CUDA-aware MPI is unavailable.
"""
function Bcast!(v::CuArray, root::Integer, comm::MPI_Comm)
    @static if !mpi_is_cuda_aware[]
        v = cpu(v)
    end
    MPI.Bcast!(v, root, comm)
    @static if !mpi_is_cuda_aware[]
        v = gpu(v)
    end
    return v
end

"""
    Reduce!(v, op, comm)

Perform `MPI.Reduce!` ensuring that CuArrays are safely transfered to
CPU if CUDA-aware MPI is unavailable.
"""
function Reduce!(v::CuArray, op::Any, root::Integer, comm::MPI_Comm)
    @static if !mpi_is_cuda_aware[]
        v = cpu(v)
    end
    MPI.Reduce!(v, op, root, comm)
    @static if !mpi_is_cuda_aware[]
        v = gpu(v)
    end
    return v
end

"""
    JuliaTaskRequest

Mostly a hack to work around asynchronous MPI calls using CuArray.

!!! warn
    This is mostly an internal implementation detail and might be removed
    without any deprecation warning
"""
mutable struct JuliaTaskRequest
    task::Any
end

wait(req::JuliaTaskRequest) = wait(req.task)
Wait!(req::JuliaTaskRequest) = wait(req)
Waitall!(reqs::Vector{JuliaTaskRequest}) = wait.(reqs)

"""
    Iallreduce!(sendbuf, recvbuf, op, comm)
    Iallreduce!(sendrecvbuf, op, comm)

Performs non-blocking elementwise reduction using the operator `op` on the buffer `sendbuf`.
`sendbuf` can also be a scalar, in which case `recvbuf` will be a value of the same type.

`recvbuf` and an MPI_Request object are returned. The value in `recvbuf` is only valid after
the request has been completed. (`MPI.Wait!`)
"""
function Iallreduce!(rbuf::RBuffer, op::Union{Op,MPI_Op}, comm::Comm)
    if rbuf.recvdata isa CuArray
        # Iallreduce! is segfaulting for CuArray. So a hack to circumvent it.
        task = @async Allreduce!(rbuf, op, comm)
        return rbuf.recvdata, JuliaTaskRequest(task)
    else
        req = Request()
        # int MPI_Iallreduce(const void *sendbuf, void *recvbuf, int count,
        #                    MPI_Datatype datatype, MPI_Op op, MPI_Comm comm,
        #                    MPI_Request * request)
        @mpichk ccall((:MPI_Iallreduce, libmpi), Cint,
                      (MPIPtr, MPIPtr, Cint, MPI_Datatype, MPI_Op, MPI_Comm, Ptr{MPI_Request}), rbuf.senddata,
                      rbuf.recvdata, rbuf.count, rbuf.datatype, op, comm, req)
        req.buffer = rbuf
        finalizer(free, req)
        return rbuf.recvdata, req
    end
end

Iallreduce!(rbuf::RBuffer, op, comm::Comm) = Iallreduce!(rbuf, Op(op, eltype(rbuf)), comm)

Iallreduce!(sendbuf, recvbuf, op, comm::Comm) = Iallreduce!(RBuffer(sendbuf, recvbuf), op, comm)

Iallreduce!(buf, op, comm::Comm) = Iallreduce!(IN_PLACE, buf, op, comm)

"""
    Ibcast!(buf, root, comm)

Non-blocking broadcast of the buffer `buf` to all processes in `comm`.

`recvbuf` and an MPI_Request object are returned. The value in `recvbuf` is only valid after
the request has been completed. (`MPI.Wait!`)
"""
function Ibcast!(buf::Buffer, root::Integer, comm::Comm)
    req = Request()
    # int MPI_Ibcast(void *buffer, int count, MPI_Datatype datatype, int root,
    #                MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Ibcast, libmpi), Cint, (MPIPtr, Cint, MPI_Datatype, Cint, MPI_Comm, Ptr{MPI_Request}), buf.data,
                  buf.count, buf.datatype, root, comm, req)
    req.buffer = buf
    finalizer(free, req)
    return buf.data, req
end

Ibcast!(buf, root::Integer, comm::Comm) = Ibcast!(Buffer(buf), root, comm)

end


text:
cuda1 = [
"sqrt",
"exp",
"log",
"sin",
"tanh",
("neg", "-", "-xi"),
("invx", "invx", "1/xi"),
("relu", "relu", "(xi>0?xi:0)"),
("sigm", "sigm", "1/(1+exp(-xi))"),
("abs", "abs", "(xi<0?-xi:xi)"),
("abs2", "abs2", "(xi*xi)"),
("sign", "sign", "(xi>0?1:xi<0?-1:0)"),
]

function cuda1src(f, j=f, ex="$f(xi)"; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
            print(s,
"""
__global__ void _$F(int n, $T *x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    $T xi = x[i];
    y[i] = $ex;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void $F(int n, $T *x, $T *y) {
    _$F<<<$BLK,$THR>>>(n,x,y);
  }
}
""")
        end
    end
end

for a in cuda1
    isa(a,Tuple) || (a=(a,))
    print(cuda1src(a...))
end

function cuda1fill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _fill_$F(int n, $T x, $T *y) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (i < n) {
    y[i] = x;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void fill_$F(int n, $T x, $T *y) {
    _fill_$F<<<$BLK,$THR>>>(n,x,y);
  }
}
""")
        end
    end
end

print(cuda1fill())

function cuda1xfill(; BLK=256, THR=256)
    sprint() do s
        for (T,F) in [("float","32"),("double","64")]
            print(s,
"""
__global__ void _xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
  int row, col, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    yidx = row + col * incy;
    y[yidx] = x;
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void xfill_$F(int nrows, int ncols, $T x, $T *y, int incy) {
    _xfill_$F<<<$BLK,$THR>>>(nrows, ncols, x, y, incy);
  }
}
""")
        end
    end
end

print(cuda1xfill())

# copy a block from x into y
function cuda1xcopy(; BLK=256, THR=256)
"""
__global__ void _xcopy(int nrows, int ncols, const char *x, int incx, char *y, int incy) {
  int row, col, xidx, yidx;
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  while (1) {
    row = i % nrows;
    col = i / nrows;
    if (col >= ncols) break;
    xidx = row + col * incx;
    yidx = row + col * incy;
    y[yidx] = x[xidx];
    i += blockDim.x * gridDim.x;
  }
}
extern "C" {
  void xcopy(int nrows, int ncols, const void *x, int incx, void *y, int incy) {
    _xcopy<<<$BLK,$THR>>>(nrows,ncols,(char*)x,incx,(char*)y,incy);
  }
}
"""
end

print(cuda1xcopy())


text:
module Ranger

# Example:
# Ranger.print_trace(stacktrace(), "run")
function print_trace(trace, stop_at::String)
    println("------- trace ---------")
    for tr in trace
        if string(tr.func) == stop_at
            break
        end
        println(tr)
    end
    println("-----------------------")
end

include("geometry/geometry.jl")
include("math/math.jl")
include("world.jl")
include("nodes/abstracts.jl")
include("events/events.jl")
include("animation/animation.jl")
include("rendering/rendering.jl")
include("nodes/nodes.jl")
include("particles/particles.jl")

include("nodes/custom/custom.jl")

include("engine/engine.jl")

end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# BigWig Overlap
# ==============

struct OverlapIterator
    reader::Reader
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32
end

function Base.eltype(::Type{OverlapIterator})
    return Record
end

function Base.IteratorSize(::Type{OverlapIterator})
    return Base.SizeUnknown()
end

function GenomicFeatures.eachoverlap(reader::Reader, interval::Interval)
    if haskey(reader.chroms, seqname(interval))
        id, _ = reader.chroms[seqname(interval)]
    else
        id = typemax(UInt32)
    end
    return OverlapIterator(reader, id, leftposition(interval) - 1, rightposition(interval))
end

mutable struct OverlapIteratorState
    # inflating data stream
    stream::IOBuffer
    data::Vector{UInt8}
    done::Bool
    header::SectionHeader
    record::Record
    blocks::Vector{BBI.Block}
    current_block::Int
    n_records::UInt16
    current_record::UInt16
end

function Base.iterate(iter::OverlapIterator)
    data = Vector{UInt8}(undef, iter.reader.header.uncompress_buf_size)
    blocks = BBI.find_overlapping_blocks(iter.reader.index, iter.chromid, iter.chromstart, iter.chromend)
    # dummy header
    header = SectionHeader(0, 0, 0, 0, 0, 0, 0, 0)
    state = OverlapIteratorState(IOBuffer(), data, false, header, Record(), blocks, 1, 0, 0)
    return iterate(iter, state)
end

function Base.iterate(iter::OverlapIterator, state::OverlapIteratorState)
    advance!(iter, state)
    if state.done
        return nothing
    end
    return copy(state.record), state
end

function advance!(iter::OverlapIterator, state::OverlapIteratorState)
    while true
        # find a section that has at least one record
        while state.current_record == state.n_records && state.current_block ≤ lastindex(state.blocks)
            block = state.blocks[state.current_block]
            seek(iter.reader.stream, block.offset)
            size = BBI.uncompress!(state.data, read(iter.reader.stream, block.size))
            state.stream = IOBuffer(state.data[1:size])
            state.header = read(state.stream, SectionHeader)
            state.current_block += 1
            state.n_records = state.header.itemcount
            state.current_record = 0
        end
        if state.current_record == state.n_records && state.current_block > lastindex(state.blocks)
            state.done = true
            return state
        end

        # read a new record
        _read!(iter.reader, state, state.record)
        if overlaps(state.record, iter.chromid, iter.chromstart, iter.chromend)
            return state
        end
    end
end

function overlaps(record::Record, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    return record.header.chromid == chromid && !(record.chromend ≤ chromstart || record.chromstart ≥ chromend)
end


text:
using Traitor
using Base.Test

abstract Size

immutable Big <: Size; end
immutable Medium <: Size; end
immutable Small <: Size; end

Size(::Union{Type{Int32},Type{Int64}}) = Small
Size(::Type{BigInt}) = Big
Size(::Type{Int128}) = Medium

@traitor function howbig(x::Any::Big)
    "Huge!"
end
@traitor function howbig(x::::Medium)
    "So-so"
end
@test howbig(Int128(1)) == "So-so"
@test howbig(BigInt(1)) == "Huge!"

@traitor function howbig(x::Integer::Small)
    "Teensy..."
end

@test howbig(1) == "Teensy..."
# These next ones throw...
# @test howbig(Int128(1)) == "So-so"
# @test howbig(BigInt(1)) == "Huge!"
# ...until we define
@traitor howbig(x::Integer::Big) = "Huge!"
@traitor howbig(x::Integer::Medium) = "So-so"

@test howbig(1) == "Teensy..."
@test howbig(Int128(1)) == "So-so"
@test howbig(BigInt(1)) == "Huge!"

@test @inferred(supertrait(Small)) == Size
@test @inferred(supertrait(Union{Small,Medium})) == Size

abstract Fooness

immutable FooA <: Fooness ; end
immutable FooB <: Fooness ; end

Fooness(::Type{Any}) = FooA
Fooness(::Type{Int16}) = FooB

@test_throws ErrorException supertrait(Union{Small,FooB})


text:
module FMRIModels

using CategoricalArrays, DataArrays, NIfTI

import MixedModels: ReMat
export ExpandedReMat

include.(["remat.jl"])

end


text:
using OrdinaryDiffEq, DiffEqBase, DiffEqCallbacks, Test
using Random
Random.seed!(213)
CACHE_TEST_ALGS = [Euler(),Midpoint(),RK4(),SSPRK22(),SSPRK33(),
  CarpenterKennedy2N54(), HSLDDRK64(),
  CFRLDDRK64(), TSLDDRK74(),
  CKLLSRK43_2(),
  ParsaniKetchesonDeconinck3S32(),
  BS3(),BS5(),DP5(),DP8(),Feagin10(),Feagin12(),Feagin14(),TanYam7(),
  Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9(),OwrenZen3(),OwrenZen4(),OwrenZen5(),
  AutoTsit5(Rosenbrock23())]
broken_CACHE_TEST_ALGS = [ORK256(), DGLDDRK73_C(),KenCarp4()]

using InteractiveUtils

NON_IMPLICIT_ALGS = filter((x)->isconcretetype(x) && !OrdinaryDiffEq.isimplicit(x()),union(subtypes(OrdinaryDiffEq.OrdinaryDiffEqAlgorithm),subtypes(OrdinaryDiffEq.OrdinaryDiffEqAdaptiveAlgorithm)))

f = function (du,u,p,t)
  for i in 1:length(u)
    du[i] = (0.3/length(u))*u[i]
  end
end

condition = function (u,t,integrator)
  1-maximum(u)
end

affect! = function (integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()/5 + 0.25
  u[maxidx] = Θ
  u[end] = 1-Θ
  nothing
end

callback = ContinuousCallback(condition,affect!)

u0 = [0.2]
tspan = (0.0,10.0)
prob = ODEProblem(f,u0,tspan)

println("Check for stochastic errors")
for i in 1:10
  @test_nowarn sol = solve(prob,Tsit5(),callback=callback)
end

println("Check some other integrators")
sol = solve(prob,Rosenbrock23(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1
sol = solve(prob,Rosenbrock32(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1
@test_broken sol = solve(prob,KenCarp4(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1
@test_broken sol = solve(prob,TRBDF2(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1

for alg in CACHE_TEST_ALGS
  @show alg
  sol = solve(prob,alg,callback=callback,dt=1/2)
  @test length(sol[end]) > 1
end

for alg in broken_CACHE_TEST_ALGS
  @show alg
  @test_broken length(solve(prob,alg,callback=callback,dt=1/2)[end]) > 1
end

sol = solve(prob,Rodas4(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1
sol = solve(prob,Rodas5(chunk_size=1),callback=callback,dt=1/2)
@test length(sol[end]) > 1

# Force switching

function f2(du,u,p,t)
  @assert length(u) == length(du) "length(u) = $(length(u)), length(du) = $(length(du)) at time $(t)"
  for i in 1:length(u)
    if t > 10
      du[i] = -10000*u[i]
    else
      du[i] = 0.3*u[i]
    end
  end
  return du
end

function condition2(u, t, integrator)
  1-maximum(u)
end

function affect2!(integrator)
  u = integrator.u
  resize!(integrator,length(u)+1)
  maxidx = findmax(u)[2]
  Θ = rand()
  u[maxidx] = Θ
  u[end] = 1-Θ
  nothing
end

callback = ContinuousCallback(condition2,affect2!)
u0 = [0.2]
tspan = (0.0,20.0)
prob = ODEProblem(f2,u0,tspan)
sol = solve(prob, AutoTsit5(Rosenbrock23()), callback=callback)
@test length(sol[end]) > 1


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://discourse.julialang.org/t/handling-of-package-inter-dependency/64568/24

# %%
module SCI
export fun1,fun2
fun1(x) = x
fun2(x) = x^2
end

module COM
using ..SCI
"""function to initialize a container"""
init_container() = []
"""function to read a message from a pipe, assume it reads from the argument for simplicity"""
read_message(x) = x
# function to check if instruction is authorized (omitted for simplicity)

"""function doing something based on the message received"""
function compute!(container, msg)
    if msg == :fun1
        push!(container, fun1(rand()))
    elseif msg ==:fun2
        push!(container, fun2(rand()))
    end
    return container
end
end

# Load COM+SCI utilities
using .COM

@show container = COM.init_container()
@show msg = COM.read_message(:fun1)
@show container = COM.compute!(container, msg)
@show msg = COM.read_message(:fun2)
@show container = COM.compute!(container, msg);

# %%

max_stars_repo_path: 0012/BambOoxX's MWE with minor corrections.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module ChordalExtensionGraph

export LabelledGraph, add_node!, add_edge!, add_clique!, chordal_extension

# With a `Vector{Vector{Int}}` with unsorted neighbors, computing fill-in
# would be inefficient.
# With a `Vector{Vector{Int}}` with sorted neighbors, adding an edge would be
# inefficient.
# With `Vector{Set{Int}}` both adding edges and computing fill-in is efficient.

struct Graph
    neighbors::Vector{Set{Int}}
    disabled::BitSet
end
Graph() = Graph(Set{Int}[], BitSet())
Base.broadcastable(g::Graph) = Ref(g)
Base.copy(G::Graph) = Graph(deepcopy(G.neighbors), copy(G.disabled))
function add_node!(G::Graph)
    push!(G.neighbors, Set{Int}())
    return length(G.neighbors)
end
num_nodes(G::Graph) = length(G.neighbors)
function add_edge!(G::Graph, i::Int, j::Int)
    if !(j in G.neighbors[i])
        push!(G.neighbors[i], j)
        push!(G.neighbors[j], i)
    end
    return (i, j)

end
function add_clique!(G::Graph, nodes::Vector{Int})
    for i in 1:(length(nodes) - 1)
        for j in (i + 1):length(nodes)
            add_edge!(G, nodes[i], nodes[j])
        end
    end
end
disable_node!(G::Graph, node::Int) = push!(G.disabled, node)
is_enabled(G::Graph, node::Int) = !(node in G.disabled)
enable_all_nodes!(G::Graph) = empty!(G.disabled)

"""
    neighbors(G::Graph, node::Int}

Return neighbors of `node` in `G`.
"""
function neighbors(G::Graph, node::Int)
    return G.neighbors[node]
end

function _num_edges_subgraph(G::Graph, nodes::Union{Vector{Int}, Set{Int}}, node::Int)
    neighs = neighbors(G, node)
    return count(nodes) do node
        is_enabled(G, node) && node in neighs
    end
end
function num_edges_subgraph(G::Graph, nodes::Union{Vector{Int}, Set{Int}})
    return mapreduce(+, nodes; init = 0) do node
        is_enabled(G, node) ? _num_edges_subgraph(G, nodes, node) : 0
    end
end

function num_missing_edges_subgraph(G::Graph, nodes::Union{Vector{Int}, Set{Int}})
    n = count(node -> is_enabled(G, node), nodes)
    # A clique is a completely connected graph. As such it has n*(n-1)/2 undirected
    # or equivalently n*(n-1) directed edges.
    return div(n * (n - 1) - num_edges_subgraph(G, nodes), 2)
end

"""
    fill_in(G::Graph{T}, i::T}

Return number of edges that need to be added to make the neighbors of `i` a clique.
"""
function fill_in(G::Graph, node::Int)
    return num_missing_edges_subgraph(G, neighbors(G, node))
end

"""
    is_clique(G::Graph{T}, x::Vector{T})

Return a `Bool` indication whether `x` is a clique in `G`.
"""
function is_clique(G::Graph, nodes::Vector{Int})
    return iszero(num_missing_edges_subgraph(G, nodes))
end

struct FillInCache
    graph::Graph
    fill_in::Vector{Int}
end
FillInCache(graph::Graph) = FillInCache(graph, [fill_in(graph, i) for i in 1:num_nodes(graph)])
Base.copy(G::FillInCache) = FillInCache(copy(G.graph), copy(G.fill_in))

num_nodes(G::FillInCache) = num_nodes(G.graph)
neighbors(G::FillInCache, node::Int) = neighbors(G.graph, node)
function add_edge!(G::FillInCache, i::Int, j::Int)
    ni = neighbors(G, i)
    nj = neighbors(G, j)
    if i in nj
        @assert j in ni
        return
    end
    @assert !(j in ni)
    for node in ni
        @assert node != i
        @assert node != j
        if node in nj
            G.fill_in[node] -= 1
        end
    end
    G.fill_in[i] += count(k -> is_enabled(G.graph, k), ni) - _num_edges_subgraph(G.graph, ni, j)
    G.fill_in[j] += count(k -> is_enabled(G.graph, k), nj) - _num_edges_subgraph(G.graph, nj, i)
    add_edge!(G.graph, i, j)
end
fill_in(G::FillInCache, node::Int) = G.fill_in[node]
function disable_node!(G::FillInCache, node::Int)
    for neighbor in neighbors(G, node)
        nodes = neighbors(G, neighbor)
        G.fill_in[neighbor] -= (length(nodes) - 1) - _num_edges_subgraph(G.graph, nodes, node)
    end
    disable_node!(G.graph, node)
end
is_enabled(G::FillInCache, node::Int) = is_enabled(G.graph, node)

"""
    struct LabelledGraph{T}
        n2int::Dict{T,Int}
        int2n::Vector{T}
        graph::Graph
    end

Type to represend a graph with nodes of type T.
"""
struct LabelledGraph{T}
    n2int::Dict{T, Int}
    int2n::Vector{T}
    graph::Graph
end

Base.broadcastable(g::LabelledGraph) = Ref(g)
Base.copy(G::LabelledGraph{T}) where T = LabelledGraph(copy(G.n2int), copy(G.int2n), copy(G.graph))

function LabelledGraph{T}() where T
    return LabelledGraph(Dict{T, Int}(), T[], Graph())
end

function LabelledGraph()
    error("`LabelledGraph()` constructor is not valid, you need to specify the type of nodes as type parameter, e.g. `LabelledGraph{Int}()`.")
end

function Base.show(io::IO, G::LabelledGraph{T}) where T
    println(io, "$T-LabelledGraph with")
    println(io, "Nodes:")
    for node in G.int2n
        print(io, "$node, ")
    end
    if !isempty(G.int2n)
        println(io, "")
    end
    println(io, "Edges:")
    for (x, i) in G.n2int
        for j in G.graph.neighbors[i]
            println(io, "( $x, $(G.int2n[j]) )")
        end
    end
end

"""
    add_node!(G::LabelledGraph{T}, i::T) where T

Add the node i to graph G.
If i is already a node of G, only return the reference.
"""
function add_node!(G::LabelledGraph{T}, i::T) where T
    if haskey(G.n2int, i)
        idx = G.n2int[i]
    else
        push!(G.int2n, i)
        idx = add_node!(G.graph)
        @assert idx == length(G.int2n)
        G.n2int[i] = idx
    end
    return idx
end

"""
    add_edge!(G::LabelledGraph{T}, i::T, j::T) where T

Add the unweighted edge (i, j) to graph G. Duplicate edges are not taken into account.
"""
function add_edge!(G::LabelledGraph{T}, xi::T, xj::T) where T
    add_edge!(G.graph, add_node!(G, xi), add_node!(G, xj))
    return (xi, xj)
end

"""
    add_edge!(G::LabelledGraph{T}, e::Tuple{T,T}) where T

Add the unweighted edge e to graph G. Duplicate edges are not taken into account.
"""
function add_edge!(G::LabelledGraph{T}, e::Tuple{T,T}) where T
    add_edge!(G, first(e), last(e))
end

"""
    add_clique!(G::LabelledGraph{T}, x::Vector{T}) where T

Add all elements of x as nodes to G and add edges such that x is fully
connected in G.
"""
function add_clique!(G::LabelledGraph{T}, x::Vector{T}) where T
    add_clique!(G.graph, add_node!.(G, x))
end

abstract type AbstractGreedyAlgorithm end

struct GreedyFillIn <: AbstractGreedyAlgorithm end
cache(G::Graph, ::GreedyFillIn) = FillInCache(G)
heuristic_value(G::FillInCache, node::Int, ::GreedyFillIn) = fill_in(G, node)

function _greedy_triangulation!(G, algo::AbstractGreedyAlgorithm)
    candidate_cliques = Vector{Int}[]
    for i in 1:num_nodes(G)
        node = argmin(map(1:num_nodes(G)) do node
            if is_enabled(G, node)
                heuristic_value(G, node, algo)
            else
                typemax(Int)
            end
        end)
#=
        # look at all its neighbors. In G we want the neighbors to form a clique.
        neighbor_nodes = collect(neighbors(G, node))

        # add neighbors and node as a potentially maximal clique
        candidate_clique = [neighbor for neighbor in neighbor_nodes if is_enabled(G, neighbor)]
        push!(candidate_clique, node)
        push!(candidate_cliques, candidate_clique)

        # add edges to G to make the new candidate_clique at least potentially a clique
        for i in eachindex(neighbor_nodes)
            if is_enabled(G, i)
                for j in (i + 1):length(neighbor_nodes)
                    if is_enabled(G, j)
                        add_edge!(G, neighbor_nodes[i], neighbor_nodes[j])
                    end
                end
            end
        end
=#
	neighbor_nodes = [neighbor for neighbor in neighbors(G, node) if is_enabled(G, neighbor)]
	push!(candidate_cliques, [node, neighbor_nodes...])
	for i in eachindex(neighbor_nodes)
		for j in (i+1):length(neighbor_nodes)
			add_edge!(G, neighbor_nodes[i], neighbor_nodes[j])
		end
	end

    disable_node!(G, node)
    end
    return candidate_cliques
end

function chordal_extension(G::Graph, algo::AbstractGreedyAlgorithm)
    H = copy(G)

    candidate_cliques = _greedy_triangulation!(cache(H, algo), algo)
    enable_all_nodes!(H)

    sort!.(candidate_cliques)
    unique!(candidate_cliques)

    # check whether candidate cliques are actually cliques
    candidate_cliques = candidate_cliques[[is_clique(H, clique) for clique in candidate_cliques]]
    sort!(candidate_cliques, by = x -> length(x))
    reverse!(candidate_cliques) # TODO use `rev=true` in `sort!`.

    maximal_cliques = [first(candidate_cliques)]
    for clique in Iterators.drop(candidate_cliques, 1)
        if all(other_clique -> !(clique ⊆ other_clique), maximal_cliques)
            push!(maximal_cliques, clique)
        end
    end

    if length(Set(Iterators.flatten(maximal_cliques))) != num_nodes(H)
        error("Maximal cliques do not cover all nodes.")
    end

    return H, maximal_cliques
end

"""
    chordal_extension(G::LabelledGraph{T})

Return a chordal extension of G and the corresponding maximal cliques.

The algoritm is Algorithm 3 in [BA10] with the GreedyFillIn heuristic of Table I.

[BA10] Bodlaender, Hans L., and Arie MCA Koster.
Treewidth computations I. Upper bounds.
Information and Computation 208, no. 3 (2010): 259-275.
Utrecht University, Utrecht, The Netherlands www.cs.uu.nl
"""
function chordal_extension(G::LabelledGraph{T}, algo::AbstractGreedyAlgorithm) where T
    H, cliques = chordal_extension(G.graph, algo)
    return LabelledGraph{T}(G.n2int, G.int2n, H), [[G.int2n[i] for i in clique] for clique in cliques]
end
end #module


text:
module Models

export IncompressibleModel, NonDimensionalModel, Clock, tick!, state

using Adapt

using Oceananigans.Architectures
using Oceananigans.Fields
using Oceananigans.Coriolis
using Oceananigans.Buoyancy
using Oceananigans.TurbulenceClosures
using Oceananigans.BoundaryConditions
using Oceananigans.Solvers
using Oceananigans.Forcing
using Oceananigans.Utils

"""
    AbstractModel

Abstract supertype for models.
"""
abstract type AbstractModel end

"""
    state(model)

Returns a `NamedTuple` with fields `velocities, tracers, diffusivities, tendencies`
corresponding to `NamedTuple`s of `OffsetArray`s that reference each of the field's data.
"""
@inline state(model) = (   velocities = datatuple(model.velocities),
                              tracers = datatuple(model.tracers),
                        diffusivities = datatuple(model.diffusivities))

include("clock.jl")
include("incompressible_model.jl")
include("non_dimensional_model.jl")
include("show_models.jl")

end


text:
export JuliaContext, create_function, call_function, isboxed, isghosttype

"""
    JuliaContext()

Returns the (session-bound) LLVM context used by the Julia compiler.
"""
function JuliaContext()
    ptr = cglobal(:jl_LLVMContext)
    ref = if VERSION >= v"1.3.0-DEV.550"
        ptr = convert(Ptr{LLVM.API.LLVMContextRef}, ptr)
        unsafe_load(ptr)
    else
        convert(LLVM.API.LLVMContextRef, ptr)
    end
    LLVM.Context(ref)
end

"""
    create_function(rettyp::LLVMType, argtyp::Vector{LLVMType}, [name::String])

Create an LLVM function, given its return type `rettyp` and a vector of argument types
`argtyp`. The function is marked for inlining, to be embedded in the caller's body.
Returns both the newly created function, and its type.
"""
function create_function(rettyp::LLVMType=LLVM.VoidType(JuliaContext()),
                         argtyp::Vector{<:LLVMType}=LLVMType[],
                         name::String="")
    ctx = context(rettyp)
    mod = LLVM.Module("llvmcall", ctx)

    ft = LLVM.FunctionType(rettyp, argtyp)
    f = LLVM.Function(mod, name, ft)
    push!(function_attributes(f), EnumAttribute("alwaysinline", 0, ctx))
    linkage!(f, LLVM.API.LLVMPrivateLinkage)

    return f, ft
end

"""
    call_function(f::LLVM.Function, rettyp::Type, argtyp::Type, args::Expr)

Generate a call to an LLVM function `f`, given its return type `rettyp` and a tuple-type for
the arguments. The arguments should be passed as an expression yielding a tuple of the
argument values (eg. `:((1,2))`), which will be splatted into the call to the function.
"""
function call_function(llvmf::LLVM.Function, rettyp::Type=Nothing, argtyp::Type=Tuple{},
                       args::Expr=:())
    ref = LLVM.ref(llvmf)
    quote
        Base.@_inline_meta
        Base.llvmcall($(convert(Ptr{Cvoid},ref)), $rettyp, $argtyp, $args...)
    end
end

"""
    isboxed(typ::Type)

Return if a type would be boxed when instantiated in the code generator.
"""
function isboxed(typ::Type)
    isboxed_ref = Ref{Bool}()
    if VERSION >= v"1.5.0-DEV.393"
        ccall(:jl_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Ptr{Bool}), typ, isboxed_ref)
    else
        ccall(:julia_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Ptr{Bool}), typ, isboxed_ref)
    end
    return isboxed_ref[]
end

"""
    convert(LLVMType, typ::Type)

Convert a Julia type `typ` to its LLVM representation. Fails if the type would be boxed.
"""
function Base.convert(::Type{LLVMType}, typ::Type, allow_boxed::Bool=false)
    isboxed_ref = Ref{Bool}()
    if VERSION >= v"1.5.0-DEV.393"
        llvmtyp = LLVMType(ccall(:jl_type_to_llvm, LLVM.API.LLVMTypeRef,
                                 (Any, Ptr{Bool}), typ, isboxed_ref))
    else
        llvmtyp = LLVMType(ccall(:julia_type_to_llvm, LLVM.API.LLVMTypeRef,
                                 (Any, Ptr{Bool}), typ, isboxed_ref))
    end
    if !allow_boxed && isboxed_ref[]
        error("Conversion of boxed type $typ is not allowed")
    end
    return llvmtyp
end

"""
    isghosttype(t::Type)
    isghosttype(T::LLVMType)

Check if a type is a ghost type, implying it would not be emitted by the Julia compiler.
This only works for types created by the Julia compiler (living in its LLVM context).
"""
isghosttype(@nospecialize(T::LLVMType)) = T == LLVM.VoidType(JuliaContext()) || isempty(T)
isghosttype(@nospecialize(t::Type)) = isghosttype(convert(LLVMType, t, true))