text:
using Test
using ModifiedHankelFunctionsOfOrderOneThird
using SpecialFunctions

# SCL1945:
# "Tables of the modified Hankel function of order one-third and of their derivatives"
# Harvard University Press, 1945.

@testset "Power series: zeros of h₂" begin
    h1, h2, h1p, h2p = modifiedhankel(complex(-1.16905371, 2.02486041))
    @test h1 ≈ complex(-0.34342750, 0.59483387)
    @test h1p ≈ complex(-0.06957489, -1.06099210)
    @test h2p ≈ complex(-1.83769223, 1.06099210)

    h1, h2, h1p, h2p = modifiedhankel(complex(-2.04397472, 3.54026807))
    @test h1 ≈ complex(0.29985265, -0.51936003)
    @test h1p ≈ complex(0.03621447, 1.21517637)
    @test h2p ≈ complex(2.10474721, -1.21517637)

    h1, h2, h1p, h2p = modifiedhankel(complex(-2.76027991, 4.78094505))
    @test h1 ≈ complex(-0.27833328, 0.48208739) atol=1e-7
    @test h1p ≈ complex(-0.02507370, -1.30912788)
    @test h2p ≈ complex(-2.26747600, 1.30912788)
end

@testset "Power series: zeros of h₂′" begin
    h1, h2, h1p, h2p = modifiedhankel(complex(-0.50939649, 0.88230059))
    @test h1 ≈ complex(-0.63166599, -0.52691131)
    @test h2 ≈ complex(0, 1.62098891)
    @test h1p ≈ complex(0.89913968, 0)

    h1, h2, h1p, h2p = modifiedhankel(complex(-1.62409879, 2.81302162))
    @test h1 ≈ complex(0.53589957, 0.33980741)
    @test h2 ≈ complex(0, -1.26801270)
    @test h1p ≈ complex(-1.14943284, 0)

    h1, h2, h1p, h2p = modifiedhankel(complex(-2.41004961, 4.17432837))
    @test h1 ≈ complex(-0.49173751, -0.29946086) atol=1e-7
    @test h2 ≈ complex(0, 1.15117521)
    @test h1p ≈ complex(1.26609349, 0)
end

@testset "Power series: spot checks" begin
    h1, h2, h1p, h2p = modifiedhankel(complex(2., 0.))
    @test h1 ≈ complex(0.60991262, 0.36822576)
    @test h1p ≈ complex(-0.59922801, 0.83306447)
    @test h2 ≈ complex(0.60991262, -0.36822576)
    @test h2p ≈ complex(-0.59922801, -0.83306447)

    h1, h2, h1p, h2p = modifiedhankel(complex(4.5, 2.2))
    @test h1 ≈ complex(-0.00170177, -0.00480825) atol=1e-7
    @test h1p ≈ complex(0.01152459, -0.00111977) atol=1e-7
    @test h2 ≈ complex(-7.26047008, 63.40667850)
    @test h2p ≈ complex(133.64890371, 44.66625637)

    h1, h2, h1p, h2p = modifiedhankel(complex(1.1, 4.1))
    @test h1 ≈ complex(-0.00202396, 0.00177285) atol=1e-7
    @test h1p ≈ complex(-0.00042308, -0.00566719) atol=1e-7
    @test h2 ≈ complex(-131.37854842, -8.34958396)
    @test h2p ≈ complex(-176.03846956, 196.31112124)
end

@testset "Comparison with Airy and Bessel functions" begin
    z = complex(1.345,-0.245)
    T = 2/3*z^(3/2)
    k = (3/2)^(2/3)*(1 - im*sqrt(3)/3)
    β = 2^(4/3)*3^(1/6)

    h1, h2, h1p, h2p = modifiedhankel(z)
    h1m, h2m, h1mp, h2mp = modifiedhankel(-z)
    h1c, h2c, h1cp, h2cp = modifiedhankel(conj(z))

    ai, bi = airyai(z), airybi(z)
    aim, bim = airyai(-z), airybi(-z)
    aip, bip = airyaiprime(z), airybiprime(z)

    # See SCL1945 (22), (25)
    @test h2 ≈ conj(h1c)
    @test h2p ≈ conj(h1cp)

    # SCL1945 (18), (19)
    @test h1 ≈ k*(aim - im*bim)
    @test h2 ≈ conj(k)*(aim + im*bim)

    # SCL1945 (20), (21)
    @test ai ≈ 1/2k*h1m + 1/(2*conj(k))*h2m
    @test bi ≈ im/2k*h1m - im/(2*conj(k))*h2m

    # SCL1945 (23), (24)
    @test aip ≈ -1/(2k)*h1mp - 1/(2*conj(k))*h2mp
    @test bip ≈ -im/(2k)*h1mp + im/(2*conj(k))*h2mp

    # SCL1945 (26)
    h1t = modifiedhankel(cis(-2π/3)*z)[1]
    h2t = modifiedhankel(cis(2π/3)*z)[2]
    @test -im*β*aim ≈ h1t
    @test -im*β*aim ≈ -h2t
    @test -im*β*aim ≈ -cis(π/3)*(cis(π/3)*h1 + h2)

    # SCL1945 (27)
    h1t = modifiedhankel(conj(z)*exp(2π*im/3))[1]
    @test h1t ≈ -h2c
    @test h1t ≈ -conj(h1)

    # SCL1945 (28)
    h1tp = modifiedhankel(z*exp(2π*im/3))[3]
    @test h1tp ≈ cis(π/3)*h2p
    @test h1tp ≈ cis(π/3)*conj(h1cp)

    # SCL1945
    @test besselj(1/3,T) ≈ 1/2*T^(-1/3)*(h1 + h2)  # (36)
    @test bessely(1/3,T) ≈ 1/2im*T^(-1/3)*(h1 - h2)  # (37)
    @test hankelh1(1/3,T) ≈ T^(-1/3)*h1  # (38)
    @test hankelh2(1/3,T) ≈ T^(-1/3)*h2  # (39)
    @test besselj(-1/3,T) ≈ 1/2*T^(-1/3)*(cis(π/3)*h1 + cis(-π/3)*h2)  # (40)
    @test bessely(-1/3,T) ≈ 1/2im*T^(-1/3)*(cis(π/3)*h1 - cis(-π/3)*h2)  # (41)
    @test hankelh1(-1/3,T) ≈ cis(π/3)*T^(-1/3)*h1  # (42)
    @test hankelh2(-1/3,T) ≈ cis(-π/3)*T^(-1/3)*h2  # (43)
    @test besselj(2/3,T) ≈ 1/2*(2/3)^(1/3)*T^(-2/3)*(cis(-2π/3)*h1p + cis(2π/3)*h2p)  # (44)
    @test bessely(2/3,T) ≈ 1/2im*(2/3)^(1/3)*T^(-2/3)*(cis(-2π/3)*h1p - cis(2π/3)*h2p)  # (45)
    @test hankelh1(2/3,T) ≈ (2/3)^(1/3)*cis(-2π/3)*T^(-2/3)*h1p  # (46)
    @test hankelh2(2/3,T) ≈ (2/3)^(1/3)*cis(2π/3)*T^(-2/3)*h2p  # (47)
    @test besselj(-2/3,T) ≈ 1/2*(2/3)^(1/3)*T^(-2/3)*(h1p + h2p)  # (48)
    @test bessely(-2/3,T) ≈ 1/2im*(2/3)^(1/3)*T^(-2/3)*(h1p - h2p)  # (49)
    @test hankelh1(-2/3,T) ≈ (2/3)^(1/3)*T^(-2/3)*h1p  # (50)
    @test hankelh2(-2/3,T) ≈ (2/3)^(1/3)*T^(-2/3)*h2p  # (51)
end

@testset "Asymptotic expansion" begin
    # From Table 2, SCL1945

    @test ModifiedHankelFunctionsOfOrderOneThird.C(1) == 0.10416666666666667
    @test ModifiedHankelFunctionsOfOrderOneThird.C(2) == 0.08355034722222222
    @test ModifiedHankelFunctionsOfOrderOneThird.C(3) == 0.12822657455632716
    @test ModifiedHankelFunctionsOfOrderOneThird.C(4) == 0.29184902646414046
    @test ModifiedHankelFunctionsOfOrderOneThird.C(5) == 0.88162726744375765

    @test ModifiedHankelFunctionsOfOrderOneThird.C(6) ≈ 3.321408281862768 atol=1e-14
    @test ModifiedHankelFunctionsOfOrderOneThird.C(7) ≈ 14.9957629868626 atol=1e-12
    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(8)) ≈ 78.923013011587 atol=1e-11
    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(9)) ≈ 474.451538868 atol=1e-8
    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(10)) ≈ 3207.490091 atol=1e-5

    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(11)) ≈ 24086.5496 atol=1e-3
    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(12)) ≈ 198923.12 atol=1e-1
    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(13)) ≈ 1791902.0 atol=1
    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(14)) ≈ 17484377. atol=1

    @test ModifiedHankelFunctionsOfOrderOneThird.C(BigInt(25)) ≈ 7.3900049415704853993e+19
end


text:
## Types

## A container for our counters
mutable struct AMVW_Counter
    zero_index::Int
    start_index::Int
    stop_index::Int
    it_count::Int
    tr::Int
end

## Rotators Our rotators have field names c, s where c nad s are
## either T or Complex{T} It proved to be faster to have immutable
## rotators, rather than mutable ones so there are no "setter" methods
## anymore, a new instance (via `Rotator(c,s,i)`) should be used

abstract type CoreTransform{T} end
abstract type AbstractRotator{T} <: CoreTransform{T} end

is_diagonal(r::AbstractRotator{T}) where {T} = norm(r.s) <= eps(T)

Base.copy(a::AbstractRotator) = AbstractRotator(a.c, a.s, a.i)
## get values
vals(r::AbstractRotator{T}) where {T} = (r.c, r.s)
idx(r::AbstractRotator) = r.i

#the index is superflous for now, and a bit of a hassle to keep immutable
#but might be of help later if twisting is approached. Shouldn't effect speed, but does mean 3N storage (Q, Ct, B)
#so may be
#
struct RealRotator{T} <: AbstractRotator{T}
c::T
s::T
i::Int
RealRotator(c::T, s::T, i::Int) where {T} = new{T}(c,s,i)
RealRotator{T}() where {T} = new()
end

function adjoint(r::RealRotator)
    RealRotator(r.c, -r.s, r.i)
end

Base.one(::Type{RealRotator{T}}) where {T} = RealRotator(one(T), zero(T), 0)# end

##################################################
### Okay, now try with complex C, real s

struct ComplexRealRotator{T} <: AbstractRotator{T}
c::Complex{T}
s::T
i::Int
end

function adjoint(r::ComplexRealRotator)
    ComplexRealRotator(conj(r.c), -r.s, r.i)
end

Base.one(::Type{ComplexRealRotator{T}}) where {T} = ComplexRealRotator(complex(one(T), zero(T)), zero(T), 0)

Base.copy(a::ComplexRealRotator) = ComplexRealRotator(a.c, a.s, a.i)

##
Rotator(c::Complex{T}, s::Complex{T}, i::Int) where {T <: Real} =
    Rotator(c, real(s), i)
Rotator(c::Complex{T}, s::T, i::Int) where {T <: Real} = ComplexRealRotator(c,s,i)
Rotator(c::T, s::T, i::Int) where {T <: Real} = RealRotator(c,s,i)

# ##################################################
# ## We use two complex, rather than 3 reals here.
# ## Will be basically the ame storage, as we don't need to include a D, but not quite (12N, not 11N)

struct ComplexComplexRotator{T} <: AbstractRotator{T}
c::Complex{T}
s::Complex{T}
i::Int
end

function adjoint(r::ComplexComplexRotator)
    ComplexComplexRotator(conj(r.c), -r.s, r.i)
end

Base.one(::Type{ComplexComplexRotator{T}}) where {T} = ComplexComplexRotator(complex(one(T), zero(T)), complex(zero(T), zero(T)), 0)

##################################################
## Factorization Types ##
##################################################

# N -- degree of poly
# POLY -- [p0, p1, p2, ..., pn] (not reversed, same as p.a for Poly type
# reverse in the init_state

abstract type FactorizationType{T, ShiftType, Pencil, Twisted} end

## SingleShift(:SingleShift) -- Q, D, U, Ut
## DoubleShift(:DoubleShift) -- Q, U, V, Ut, Vt, W

## :NoPencil -- Ct, B
## :HasPencil -- [D], Ct, B, Ct1, B1, [D1]

## :NotTwisted --
## :IsTwisted -- sigma

## All have: N, POLY, REIGS, IEIGS,
## A, Bk, R, e1, e2, ctr

############### No Pencil, Not Twisted ###################################

## RDS, no pencil, not twisted
#struct Real_DoubleShift_NoPencil_NotTwisted{T} <: FactorizationType{T, Val{:DoubleShift}, Val{:NoPencil}, Val{:NotTwisted}
mutable struct Real_DoubleShift_NoPencil_NotTwisted{T} <: FactorizationType{T, Val{:DoubleShift}, Val{:NoPencil}, Val{:NotTwisted}}

N::Int
POLY::Vector{T}
##
Q::Vector{RealRotator{T}}
Ct::Vector{RealRotator{T}}  # We use C', not C here
B::Vector{RealRotator{T}}
##
REIGS::Vector{T}
IEIGS::Vector{T}
## reusable storage
U::RealRotator{T}
Ut::RealRotator{T}
V::RealRotator{T}
Vt::RealRotator{T}
W::RealRotator{T}
A::Matrix{T}    # for parts of A = QR
R::Matrix{T}    # temp storage, sometimes R part of QR
e1::Vector{T}   # eigen values e1, e2
e2::Vector{T}
ctrs::AMVW_Counter
end

function Base.convert(::Type{FactorizationType{T, Val{:DoubleShift}, Val{:NoPencil}, Val{:NotTwisted}}}, ps::Vector{T}) where {T <: AbstractFloat}

    N = length(ps) - 1

    Real_DoubleShift_NoPencil_NotTwisted(N, ps,
                                         Vector{AMVW.RealRotator{T}}(undef, N), #Q
                                         Vector{AMVW.RealRotator{T}}(undef, N), #Ct '
                                         Vector{AMVW.RealRotator{T}}(undef, N), # B
                                         zeros(T, N),  zeros(T, N), #EIGS
                                         one(RealRotator{T}), one(RealRotator{T}), #U, Ut
                                         one(RealRotator{T}), one(RealRotator{T}), #V, Vt
    one(RealRotator{T}), # W
    zeros(T, 2, 2),zeros(T, 3, 2), # A, R
    zeros(T,2), zeros(T,2),
    AMVW_Counter(0,1,N-1, 0, N-2)
    )
end

# ComplexReal Double Shift, no pencil, not twisted
mutable struct ComplexReal_SingleShift_NoPencil_NotTwisted{T} <: FactorizationType{T, Val{:SingleShift}, Val{:NoPencil}, Val{:NotTwisted}}

N::Int
POLY::Vector{Complex{T}}
Q::Vector{ComplexRealRotator{T}}
D::Vector{Complex{T}}
Ct::Vector{ComplexRealRotator{T}}  # We use C', not C here
B::Vector{ComplexRealRotator{T}}
#
REIGS::Vector{T}
IEIGS::Vector{T}
# reusable storage
U::ComplexRealRotator{T}
Ut::ComplexRealRotator{T}
A::Matrix{Complex{T}}    # for parts of A = QR
R::Matrix{Complex{T}}    # temp storage, sometimes R part of QR
e1::Vector{T}   # eigen values e1, e2, store as (re,imag)
e2::Vector{T}
ray::Bool
ctrs::AMVW_Counter
end

function Base.convert(::Type{FactorizationType{T, Val{:SingleShift}, Val{:NoPencil}, Val{:NotTwisted}}}, ps::Vector{Complex{T}}; ray=true) where {T}

    N = length(ps) - 1

    ComplexReal_SingleShift_NoPencil_NotTwisted(N, ps,
                                                Vector{ComplexRealRotator{T}}(undef,N), #Q
                                                ones(Complex{T}, N+1), # D
                                                Vector{ComplexRealRotator{T}}(undef,N), #Ct
                                                Vector{ComplexRealRotator{T}}(undef,N), #B
                                                zeros(T, N),  zeros(T, N), #EIGS
    one(ComplexRealRotator{T}), one(ComplexRealRotator{T}), #U, Ut
    zeros(Complex{T}, 2, 2),zeros(Complex{T}, 3, 2), # A R
    zeros(T,2), zeros(T,2),
    #    true,  # true for Wilkinson, false for Rayleigh.Make adjustable!
    ray,
    AMVW_Counter(0,1,N-1, 0, N-2)
    )
end

############## Has Pencil, Not twisted ####################################

## RDS, no pencil, not twisted
mutable struct Real_DoubleShift_HasPencil_NotTwisted{T} <: FactorizationType{T, Val{:DoubleShift}, Val{:HasPencil}, Val{:NotTwisted}}

N::Int
POLY::Vector{T}
##
Q::Vector{RealRotator{T}}
Ct::Vector{RealRotator{T}}  # We use C', not C here
B::Vector{RealRotator{T}}
#
Ct1::Vector{RealRotator{T}}  # W = Ct * B
B1::Vector{RealRotator{T}}   #
##
REIGS::Vector{T}
IEIGS::Vector{T}
## reusable storage
U::RealRotator{T}
Ut::RealRotator{T}
V::RealRotator{T}
Vt::RealRotator{T}
W::RealRotator{T}
#
A::Matrix{T}    # for parts of A = QR
R::Matrix{T}    # temp storage, sometimes R part of QR
e1::Vector{T}   # eigen values e1, e2
e2::Vector{T}
ctrs::AMVW_Counter
end

function Base.convert(::Type{FactorizationType{T, Val{:DoubleShift}, Val{:HasPencil}, Val{:NotTwisted}}}, ps::Vector{T}) where {T}

    N = length(ps) - 1

    Real_DoubleShift_HasPencil_NotTwisted(N, ps,
                                          Vector{RealRotator{T}}(undef, N), #Q
                                          Vector{RealRotator{T}}(undef, N), #Ct
                                          Vector{RealRotator{T}}(undef, N), #B
                                          Vector{RealRotator{T}}(undef, N), #Ct1
                                          Vector{RealRotator{T}}(undef, N), #B1
    zeros(T, N),  zeros(T, N), #EIGS
    one(RealRotator{T}), one(RealRotator{T}),
    one(RealRotator{T}), one(RealRotator{T}),
    one(RealRotator{T}), #U,U',V,V',W
    zeros(T, 2, 2), zeros(T, 3, 2), # A R
    zeros(T,2), zeros(T,2),
    AMVW_Counter(0,1,N-1, 0, N-2)
    )
end

# ComplexReal Double Shift, has pencil, not twisted

mutable struct ComplexReal_SingleShift_HasPencil_NotTwisted{T} <: FactorizationType{T, Val{:SingleShift}, Val{:HasPencil}, Val{:NotTwisted}}

N::Int
POLY::Vector{Complex{T}}
Q::Vector{ComplexRealRotator{T}}
D::Vector{Complex{T}}
Ct::Vector{ComplexRealRotator{T}}  # We use C', not C here
B::Vector{ComplexRealRotator{T}}
# pencil part
D1::Vector{Complex{T}}
Ct1::Vector{ComplexRealRotator{T}}  # not inverses
B1::Vector{ComplexRealRotator{T}}
#
REIGS::Vector{T}
IEIGS::Vector{T}
## reusable storage
U::ComplexRealRotator{T}
Ut::ComplexRealRotator{T}
#
A::Matrix{Complex{T}}    # for parts of A = QR
R::Matrix{Complex{T}}    # temp storage, sometimes R part of QR
e1::Vector{T}   # eigen values e1, e2, store as (re,imag)
e2::Vector{T}
#
ray::Bool
ctrs::AMVW_Counter
end

function Base.convert(::Type{FactorizationType{T, Val{:SingleShift}, Val{:HasPencil}, Val{:NotTwisted}}}, ps::Vector{Complex{T}}; ray=true) where {T}

    N = length(ps) - 1

    ComplexReal_SingleShift_HasPencil_NotTwisted(N, ps,
                                                 Vector{ComplexRealRotator{T}}(undef, N), #Q
                                                 ones(Complex{T}, N+1), # D
                                                 Vector{ComplexRealRotator{T}}(undef, N), #Ct
    Vector{ComplexRealRotator{T}}(undef, N), #B
    ones(Complex{T}, N+1), # D1
    Vector{ComplexRealRotator{T}}(undef, N), #Ct1
    Vector{ComplexRealRotator{T}}(undef, N), #B1
    zeros(T, N),  zeros(T, N), #EIGS
    one(ComplexRealRotator{T}), one(ComplexRealRotator{T}), #U, Ut
    zeros(Complex{T}, 2, 2),zeros(Complex{T}, 3, 2), # A , R
    zeros(T,2), zeros(T,2),
    #    true,  # true for Wilkinson, false for Rayleigh.Make adjustable!
    ray,
    AMVW_Counter(0,1,N-1, 0, N-2)
    )
end

ext: jl


text:
module Validation

    include("betacv.jl")

    export betacv,
    intra_cluster_weights,
    inter_cluster_weights,
    weights

    export betacv_fused

    export betacv_pairwise,
    intra_cluster_weights_pairwise,
    inter_cluster_weights_pairwise,
    weights_pairwise

end # module

max_stars_repo_name: swiesend/FluxMO

max_issues_repo_name: swiesend/FluxMO

max_forks_repo_name: swiesend/FluxMO


text:
# Contains helper functions
# Functions for inserting text are based on
# https://discourse.julialang.org/t/write-to-a-particular-line-in-a-file/50179

############# Project Structure Helper Function ###############
# Collection of functions creating the latex files and
# directories. Inside the "target directory" which is equivalent
# the `path` varible, there are some files:
# `main.tex` - Main file to run build the pdf;
# `julia_font.tex` - Font information;
# `julia_listings.tex` - Listings configurations;
# `julia_listings_unicode.tex` - Listings unicodes;
# `preface.tex` - Preface when using a book latex template.
# There are four folders:
# `figures/` - Location where plots are saved;
# `fonts/` - Location of the `.ttf` files for the Julia-Mono font;
# `frontmatter/`- Files used in the book latex template;
# `notebooks/`  - Contains the latex files for each converted notebook.

"""
    createfolders(targetdir="./")
Creates the directories for the latex project.
"""
function createfolders(targetdir="./")
    if !isdir(targetdir)
        mkpath(targetdir * "/notebooks")
        mkpath(targetdir * "/figures")
        mkpath(targetdir * "/frontmatter")
    else
        if !isdir(targetdir * "/notebooks")
            mkpath(targetdir * "/notebooks")
        end
        if !isdir(targetdir * "/figures")
            mkpath(targetdir * "/figures")
        end
        if !isdir(targetdir * "/frontmatter")
            mkpath(targetdir * "/frontmatter")
        end
    end
end

"""
    downloadfonts(targetdir="./"; fontpath=nothing)
Downloads the Julia-Mono fonts necessary. If the user already
has the fonts installed, he can provide `fontpath` to
be used, instead of downloading `.ttf` files.
Note that in such case, the `fonts/` dictory will be empty.
If you use a font manager, the Julia-Mono font might
be found in a path such as the one below:
```
fontpath = "/home/username/.local/share/fonts/Unknown Vendor/TrueType/JuliaMono/"
```
"""
function downloadfonts(targetdir="./"; fontpath=nothing)
    if fontpath === nothing
        if !isdir(targetdir * "/fonts")
            mkpath(targetdir * "/fonts")
            fontsourcepath = joinpath(@__DIR__, "../templates/fonts/")
            for font in readdir(fontsourcepath)
                cp(joinpath(fontsourcepath, font),
                   joinpath(targetdir * "/fonts/", font))
            end
        end
        # the code below sets font path in the julia_font.tex file
        julia_font_tex = targetdir * "/julia_font.tex"
        writetext(julia_font_tex, "./fonts/,", 6)
    else
        julia_font_tex = targetdir * "/julia_font.tex"
        writetext(julia_font_tex, fontpath * "/,", 6)
    end
end

"""
    createproject(targetdir="./", template=:book, fontpath=nothing)
This function just call all the other auxiliary functions
in order to create the proper files and directories.
Note that first the folders must be created, then
the main latex template files, followed by the auxiliary
latex listings, and finally, the fonts.
"""
function createproject(targetdir="./", template=:book, fontpath=nothing)
    createfolders(targetdir)
    createtemplate(targetdir, template)
    createauxiliarytex(targetdir)
    downloadfonts(targetdir, fontpath=fontpath)
end

############# Text Helper Function ###############
# Collection of functions for inserting text in
# specific locations of a file

"""
    skiplines(io::IO, n)
Helper function for skipping lines when reading a file.
"""
function skiplines(io::IO, n)
    i = 1
    while i <= n
        eof(io) && error("File contains less than $n lines")
        i += read(io, Char) === '\n'
    end
end

"""
    writetext(file::String, text::String, linenumber::Integer, endline=true)
Writes a `text` to a `file` in an specific `linenumber`. By default,
the text is appended at the end of the line (`endline = true`). If
`endline` is set to false, then the text is actually written in
the beggining of the line. The text is appended, and does not overwrite
what is already written in the file.
"""
function writetext(file::String, text::String, linenumber::Integer, endline=true)
    f = open(file, "r+");
    if endline
        skiplines(f, linenumber);
        skip(f, -1)
    else
        skiplines(f, linenumber - 1);
    end
    mark(f)
    buf = IOBuffer()
    write(buf, f)
    seekstart(buf)
    reset(f)
    print(f, text);
    write(f, buf)
    close(f)
end

"""
    insertlineabove(file::String, text::String, linenumber::Integer)
Inserts a line of `text` in a `file` above the `linenumber`.
"""
function insertlineabove(file::String, text::String, linenumber::Integer)
    if linenumber == 1
        writetext(file, text * "\n", linenumber, false)
    else
        writetext(file, "\n" * text, linenumber - 1)
    end
end

"""
    insertlinebelow(file::String, text::String, linenumber::Integer)
Inserts a line of `text` in a `file` below the `linenumber`.
"""
function insertlinebelow(file::String, text::String, linenumber::Integer)
    writetext(file, "\n" * text, linenumber)
end

################ Parsing Nested Dictionary #####################
"""
    nestedget(dict::Dict, keys::AbstractVector, default = nothing)
Checks whether `dict["key1"]["key2"]...["keyn"]` exists,
returing it's value or the `default`.
"""
function nestedget(dict::Dict, keys::AbstractVector, default=nothing)
    if isempty(keys)
        return dict
    end
    key = popfirst!(keys)
    if haskey(dict, key)
        return nestedget(dict[key], keys)
    end
    return default
end

"""
    nestedget(non_dict, keys::AbstractVector, default = nothing)
This function is an auxiliary used by `nestedget` above
in order to extrac the value from the last key.
"""
function nestedget(non_dict, keys::AbstractVector, default = nothing)
    if isempty(keys)
        return non_dict
    end
    return default
end


text:
testEventTimes = [1.0,2.0,3.0]

@testset "Count Events" begin

    @testset "Basic" begin

        testParentVector = [0,1,0]
        testCountEvents = HawkesProcesses.countEvent(testParentVector, testEventTimes)
        inferTests && @inferred HawkesProcesses.countEvent(testParentVector, testEventTimes)
        @test isa(testCountEvents, Tuple)
        @test length(testCountEvents) == 3
        @test testCountEvents[1] == 2
        @test testCountEvents[2] == [1,0,0]
        @test testCountEvents[3] == [1]
    end

    @testset "No Children" begin

        testCountEvents_nochildren = HawkesProcesses.countEvent([0,0,0], testEventTimes)
        inferTests && @inferred HawkesProcesses.countEvent([0,0,0], testEventTimes)
        @test isa(testCountEvents_nochildren, Tuple)
        @test length(testCountEvents_nochildren) == 3
        @test testCountEvents_nochildren[1] == 3
        @test testCountEvents_nochildren[2] == [0,0,0]
        @test isempty(testCountEvents_nochildren[3])
    end

end

# testEventTimes2 = [1,2,3] + .5
# testParentVector2 = [0,0,1]
#
# testCountEvents_hierarchical = map(
#     (x,y) -> HawkesProcesseses.countEvent(x,y),
#     [testParentVector, testParentVector2],
#     [testEventTimes, testEventTimes2],)


text:
const ħoec = ustrip(u"T*nm^2",ħ/e)
const μB = ustrip(u"T^-1*meV",μ_B)
const σ0τz = @SMatrix[1 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1]
const σ0τ0 = @SMatrix[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
const σzτ0 = @SMatrix[1 0 0 0; 0 -1 0 0; 0 0 1 0; 0 0 0 -1]
const σzτz = @SMatrix[1 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1]
const σyτy = @SMatrix[0 0 0 -1; 0 0 1 0; 0 1 0 0; -1 0 0 0]
const σyτz = @SMatrix[0 -im 0 0; im 0 0 0; 0 0 0 im; 0 0 -im 0]
const σyτ0 = @SMatrix[0 -im 0 0; im 0 0 0; 0 0 0 -im; 0 0 im 0]
const σxτz = @SMatrix[0 1 0 0; 1 0 0 0; 0 0 0 -1; 0 0 -1 0]
const σxτ0 = @SMatrix[0 1 0 0; 1 0 0 0; 0 0 0 1; 0 0 1 0]
const σ0 = SA[1 0; 0 1]
const σx = SA[0 1; 1 0]
const σy = SA[0 -im; im 0]
const σz = SA[1 0; 0 -1]

# Default presets
@with_kw struct Params @deftype Float64
    # Units: nm, meV
    # Params AB bernal bilayer obtained from DFT calculations:
    # ref -> https://arxiv.org/pdf/1511.06706.pdf
    scale = 1                  #scaling parameter
    a0 = 0.246 * scale
    dinter = 1.36 * a0/scale  #Note that we leave dinter unaltered by scale
    dintra = a0/sqrt(3)
    dinter_warp = sqrt(dintra^2 + dinter^2)
    Ls = 0                      #units of a0
    Ws = 0
    Ln = 30                     #units of a0
    W = 30                      #units of a0
    d = 0                       # transition length of NS contact

    t0 =  1e3 * 2.6 / scale     #in meV
    t1 = 1e3 * .36              #in meV
    t3 = 1e3 * .28 / scale      #in meV
    μN = 0 * t0
    μS = t0 / 3.25              # t0/3 > μS > 2 * t1
    τ = 1                       # NS transparency
    barrier = 0                 # NS barrier

    U_dimer = 0 * 0.015         #Energy difference between dimer and non dimer sites:
                                # https://arxiv.org/pdf/1901.01332.pdf
    U = 0                       #Interlayer bias
    Δ = 1.3                     #That of SC MoRe: ref ->  https://arxiv.org/abs/2006.05522
    EZ::SVector{3,Float64} = [0,20,0]
    E::SVector{3,Float64} = [0,0,0]
    λ = 30                      #Between [-2.5,2.5] meV:
                                #ref -> https://arxiv.org/abs/1905.08688
    α = 20                      #Rashba coef
    g = 2
end

"""
    latBLG_unbounded(p = Params()
creates an unbounded lattice for the BLG in Bernal stacking
see: Params()
"""
function latBLG_unbounded(p = Params())
    @unpack a0, dinter = p
    sAbot = sublat((0.0,-1.0a0/sqrt(3.0), - dinter/2); name = :Ab)
    sBbot = sublat((0.0, 0.0a0/sqrt(3.0), - dinter/2); name = :Bb)
    sAtop = sublat((0.0, 0.0a0/sqrt(3.0), + dinter/2); name = :At)
    sBtop = sublat((0.0, 1.0a0/sqrt(3.0), + dinter/2); name = :Bt)
    br = a0 * SA[cos(pi/3) sin(pi/3) 0; -cos(pi/3) sin(pi/3) 0]'
    lat = lattice(sAtop, sBtop, sAbot, sBbot; bravais = br)
    return lat
end

"""
    latBLG(p, θ)
creates an unbounded lattice for the BLG in Bernal stacking with bravais vectors rotated by
θ degrees wrt latBLG_unbounded(::Params()).
    see: latBLG_unbounded(::Params()), Params()
"""
function latBLG(p, θ)
    (; a0, dinter, Ln, W) = p
    lat0_slg = LP.honeycomb(; a0, dim = 3)
    function rotated_region(r)
        r´ = SA[r[1], floor(abs(r[2])/(0.5*√3 * a0)) * sign(r[2])*(0.5*√3 * a0), r[3]]
        r´ = SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * r´
        return abs(r´[1]) <= Ln/2 && abs(r´[2]) <= W/2
    end
    lat_slg = unitcell(lat0_slg, region = rotated_region)
    Quantica.transform!(r -> SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * r, lat_slg)
    lat_bot = lattice(Quantica.transform!(r -> r +
        SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * SA[0, -0.5*a0/sqrt(3.0), -dinter/2],
        copy(lat_slg)); names = (:Ab, :Bb), type = Float64)
    lat_top = lattice(Quantica.transform!(r -> r +
        SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * SA[0, 0.5*a0/sqrt(3.0), dinter/2],
        copy(lat_slg)); names = (:At, :Bt), type = Float64)
    return lat_top, lat_bot
end

latSLG(p = Params()) = LP.honeycomb(; a0 = p.a0, dim = Val(3))

"""
    `modelN(p = Params())`
Hamiltonian terms of the normal system (without induced superconductivity):
    see: `Params()`
"""
function modelN(p = Params())
    @unpack a0, dinter, dintra, dinter_warp, Ls, Ln, W, μN, μS, t0, t1,
            t3, Δ, λ, g, U, U_dimer, α, scale, EZ, E = p

    modelintra = hopping(-t0 * σ0, range = dintra,
            sublats = (:At=>:Bt,:Ab=>:Bb,:Bb=>:Ab,:Bt=>:At, :A=>:B, :B=>:A))

    modelinter = hopping(t1 * σ0, range = dinter, sublats = (:At=>:Bb,:Bb=>:At)) +
                 hopping(t3 * σ0, range = dinter_warp, sublats = (:Bt=>:Ab, :Ab=>:Bt))

    hopI(λ, dr) = λ/(2*3*sqrt(3)) * im * cos(3*atan(dr[2],dr[1])) * σz

    modelIsing = hopping((r, dr) -> hopI(λ,dr) * sign(r[3]),
        sublats = (:At, :Bt ,:Ab, :Bb) .=> (:At, :Bt ,:Ab, :Bb), range = a0)

    modelKaneMele =
        hopping((r, dr) -> hopI(λ, dr), sublats = :A=>:A, range = a0) -
        hopping((r, dr) -> hopI(λ, dr), sublats = :B=>:B, range = a0)

    hopR(α, dr) = 2α/3 * im * (dr[2]/dintra * σx - dr[1]/dintra * σy)

    modelRashba =
        hopping((r,dr) -> hopR(α, dr), range = a0/√3,
            sublats = (:At=>:Bt, :Bt=>:At, :A=>:B, :B=>:A)) -
        hopping((r,dr) -> hopR(α, dr), range = a0/√3, sublats = (:Ab=>:Bb, :Bb=>:Ab))

    modelonsiteN = onsite(-μN * σ0 + EZ'*SA[σx, σy, σz])

    model0 = modelonsiteN + modelintra + modelinter + modelIsing +
        modelKaneMele + modelRashba

    field! = @onsite!((o, r) -> o + E'r * σ0)

    return (; model0, field!)
end

"""
    `modelS(p = Params())`
Hamiltonian terms of the proximitized system:
    see: `Params()`
"""
function modelS(p = Params())
    @unpack a0, dinter, dintra, dinter_warp, Ls, Ln, W, μN,μS , t0, t1,
            t3, Δ, λ, g, U, U_dimer, α, scale, EZ, E = p

    modelintra = hopping(-t0 * σ0τz, range = dintra,
            sublats = (:At=>:Bt,:Ab=>:Bb,:Bb=>:Ab,:Bt=>:At, :A=>:B, :B=>:A))

    modelinter = hopping(t1 * σ0τz, range = dinter, sublats = (:At=>:Bb, :Bb=>:At)) +
                 hopping(t3 * σ0τz, range = dinter_warp, sublats = (:Bt=>:Ab, :Ab=>:Bt))

    hopI(λ, dr) = λ/(2*3*sqrt(3)) * im * cos(3*atan(dr[2],dr[1])) * σzτ0

    modelIsing = hopping((r, dr) -> hopI(λ,dr) * sign(r[3]),
        sublats = (:At, :Bt ,:Ab, :Bb) .=> (:At, :Bt ,:Ab, :Bb), range = a0)

    modelKaneMele =
        hopping((r, dr) -> 0*hopI(λ, dr), sublats = :A=> :A, range = a0) -
        hopping((r, dr) -> 0*hopI(λ, dr), sublats = :B=> :B, range = a0)

    hopR(α, dr) = 2α/3 * im * (dr[2]/dintra * σxτ0 - dr[1]/dintra * σyτz)

    modelRashba =
        hopping((r,dr) -> hopR(α, dr), range = a0/√3,
            sublats = (:At=>:Bt, :Bt=>:At, :A=>:B, :B=>:A)) -
        hopping((r,dr) -> hopR(α, dr), range = a0/√3, sublats = (:Ab=>:Bb, :Bb=>:Ab))

    modelonsiteN = onsite(-μN * σ0τz + EZ' * SA[σxτz, σyτ0, σzτz])

    model0 = modelonsiteN + modelintra + modelinter + modelIsing +
        modelKaneMele + modelRashba

    field! = @onsite!((o, r) -> o + E'r * σ0τz)

    return (; model0, field!, modelinter)
end

############
## Deprecated
############

# function latBLG_unbounded(p, θ)
#     (; a0, dinter, Ln, W) = p
#     lat0_slg = LP.honeycomb(; a0, dim = 3)
#     function rotated_region(r)
#         r´ = SA[r[1], floor(abs(r[2])/(0.5*√3 * a0)) * sign(r[2])*(0.5*√3 * a0), r[3]]
#         r´ = SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * r´
#         return abs(r´[1]) <= Ln/2 && abs(r´[2]) <= W/2
#     end
#     lat_slg = unitcell(lat0_slg, region = rotated_region)
#     Quantica.transform!(r -> SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * r, lat_slg)
#     lat_bot = lattice(Quantica.transform!(r ->
#    r + SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] *
#          SA[0, -0.5*a0/sqrt(3.0), -dinter/2], copy(lat_slg));
#           names = (:Ab, :Bb), bravais = br)
#     lat_top = lattice(Quantica.transform!(r -> r +
#         SA[cos(θ) -sin(θ) 0; sin(θ) cos(θ) 0; 0 0 1] * SA[0, 0.5*a0/sqrt(3.0), dinter/2],
#         copy(lat_slg)); names = (:At, :Bt), bravais = br)
#     return lat_top, lat_bot
# end

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
module MedicalImagingUtils

using Images

include("./calciumscoring.jl")

export
    # Export calciumscoring.jl functions
    compute_calcium_scores

end


text:
@testset "TimeSeries" begin
    @testset "constructor" begin
        timestamps = collect(DateTime(2000):Hour(1):DateTime(2000, 2, 1))
        vals       = rand(length(timestamps))
        ts         = TimeSeries("teste", timestamps, vals)
        @test isa(ts, TimeSeries)

        timestamps = collect(DateTime(2000):Day(1):DateTime(2000, 2, 1))
        vals       = rand(length(timestamps))
        ts         = TimeSeries("teste", timestamps, vals)
        @test isa(ts, TimeSeries)

        timestamps = collect(DateTime(2000):Second(1):DateTime(2000, 1, 1, 1, 1, 1))
        vals       = rand(length(timestamps))
        ts         = TimeSeries("teste", timestamps, vals)
        @test isa(ts, TimeSeries)

        timestamps = collect(DateTime(2000):Day(1):DateTime(2000, 2, 1))
        vals       = rand(length(timestamps))[1:end - 1]
        @test_throws DimensionMismatch TimeSeries("teste", timestamps, vals)

        timestamps = [DateTime(2001); DateTime(2000, 2, 1)]
        vals       = [1, 2]
        @test_throws ErrorException TimeSeries("teste", timestamps, vals)

        timestamps = [DateTime(2001); DateTime(2001, 1, 1)]
        vals       = [1, 2]
        @test_throws ErrorException TimeSeries("teste", timestamps, vals)
    end

    @testset "arithmetic methods (+ - *)" begin
        stamps1  = [DateTime(2000, 1, i) for i = 1:2]
        vals1    = collect(1:length(stamps1))
        ts1      = TimeSeries("teste 1", stamps1, vals1)

        stamps2  = [DateTime(2000, 1, i) for i = 2:4]
        vals2    = collect(1:length(stamps2))
        ts2      = TimeSeries("teste 2", stamps2, vals2)

        stamps3  = [DateTime(2000, 1, i) for i = 1:4]
        vals3    = collect(1:length(stamps3))
        ts3      = TimeSeries("teste 3", stamps3, vals3)

        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 + ts2, TimeSeries("teste 1 + teste 2", [DateTime(2000, 1, i) for i = 1:4], [1, 3, 2, 3]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 + ts2 + ts3, TimeSeries("teste 1 + teste 2 + teste 3", [DateTime(2000, 1, i) for i = 1:4], [2, 5, 5, 7]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(- ts1, TimeSeries(string("- ", ts1.name), ts1.timestamps, -ts1.vals))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 - ts2, TimeSeries("teste 1 - teste 2", [DateTime(2000, 1, i) for i = 1:4], [1, 1, -2, -3]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(+ ts1 - ts2, TimeSeries("teste 1 - teste 2", [DateTime(2000, 1, i) for i = 1:4], [1, 1, -2, -3]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 - ts2, TimeSeries("teste 1 - teste 2", [DateTime(2000, 1, i) for i = 1:4], [1, 1, -2, -3]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 - ts2 - ts3, TimeSeries("teste 1 - teste 2 - teste 3", [DateTime(2000, 1, i) for i = 1:4], [0, -1, -5, -7]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 - ts1 + ts1, TimeSeries("teste 1 - teste 1 + teste 1", [DateTime(2000, 1, i) for i = 1:2], [1, 2]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 * ts2, TimeSeries("teste 1 * teste 2", [DateTime(2000, 1, i) for i = 1:4], [1, 2, 2, 3]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 * ts1 * ts1, TimeSeries("teste 1 * teste 1 * teste 1", [DateTime(2000, 1, i) for i = 1:2], [1, 8]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 * ts2 * ts3, TimeSeries("teste 1 * teste 2 * teste 3", [DateTime(2000, 1, i) for i = 1:4], [1, 4, 6, 12]))
        @test TimeSeriesInterface.verify_if_ts_are_equal(ts1 + ts2 * ts3, TimeSeries("teste 1 + teste 2 * teste 3", [DateTime(2000, 1, i) for i = 1:4], [2, 4, 6, 12]))
        @test TimeSeriesInterface.verify_if_ts_are_equal((ts1 + ts2) * ts3, TimeSeries("teste 1 + teste 2 * teste 3", [DateTime(2000, 1, i) for i = 1:4], [1, 6, 6, 12]))

    end

    @testset "aggregation" begin
        timestamps = DateTime.(["2020-02-23T18:00", "2020-02-23T19:00", "2020-02-24T17:00", "2020-02-24T18:00",
                                "2020-04-23T18:00", "2020-05-23T18:00", "2020-07-13T13:00"])
        vals = Vector(1:length(timestamps))
        ts = TimeSeries("teste", timestamps, vals)
        ts = hourly_to_monthly(ts)
        @test ts.name == "teste"
        @test all(ts.timestamps .== DateTime.(["2020-02-01T00:00", "2020-04-01T00:00", "2020-05-01T00:00", "2020-07-01T00:00",]))
        @test all(ts.vals .≈  [2.5, 5, 6, 7])
    end

    @testset "disaggregation" begin
        timestamps = DateTime.(["2020-02-01T00:00", "2020-04-01T00:00", "2020-05-01T00:00", "2020-07-01T00:00",])
        vals = Vector(1:length(timestamps))
        ts = TimeSeries("teste", timestamps, vals)
        ts = monthly_to_hourly(ts)
        @test ts.name == "teste"
        expected_time_stamps = vcat(collect(DateTime(2020, 2, 1):Hour(1):DateTime(2020, 2, 29, 23)),
                    collect(DateTime(2020, 4, 1):Hour(1):DateTime(2020, 4, 30, 23)),
                    collect(DateTime(2020, 5, 1):Hour(1):DateTime(2020, 5, 31, 23)),
                    collect(DateTime(2020, 7, 1):Hour(1):DateTime(2020, 7, 31, 23)))
        @test all(ts.timestamps .== expected_time_stamps)
        expected_vals = ones(length(expected_time_stamps))
        expected_vals[1:696] .= 1
        expected_vals[697:1416] .= 2
        expected_vals[1417:2160] .= 3
        expected_vals[2161:end] .= 4
        @test all(ts.vals .≈ expected_vals)
    end
end


text:
include("../ode_test_functions/ode_test_functions.jl")
include("../src/ode_module.jl")

module otest
using ..odesample
using ..ode
using PyPlot
PyPlot.pygui(true)

time_interval   = ode.TimeInterval(5.0)
initial_values  = [0.0, 1.0]
options         = ode.OdeOptions(step_size = 1e-3)

params = odesample.HarmonicOscillatorOptions(
    m = 1.0,
    d = 0.5,
    c = 10.0)
fcn(t, q) = odesample.harmonic_oscillator(t, q, params)

ode_res = ode.ode_solver_implicit_euler(fcn, time_interval, initial_values, options)

PyPlot.figure()
PyPlot.subplot(2,1,1)
PyPlot.grid()
PyPlot.plot( ode_res.t, ode_res.q[:, 1] )

PyPlot.subplot(2,1,2)
PyPlot.grid()
PyPlot.plot( ode_res.t, ode_res.q[:, 2] )

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Aliases and short functions go here
# ===========================================================================

# Length of data vector in each channel in S
Lx(S::GphysData) = [length(S.x[i]) for i in 1:S.n]

# Change a string separator to a path separator
change_sep(S::Array{String,1}) = [replace(i, "/" => Base.Filesystem.pathsep()) for i in S]

# Check that fields are preserved from S1 to S2
test_fields_preserved(S1::GphysData, S2::GphysData, x::Int, y::Int) =
  @test(minimum([getfield(S1,f)[x]==getfield(S2,f)[y] for f in datafields]))
test_fields_preserved(S1::SeisChannel, S2::GphysData, y::Int) =
  @test(minimum([getfield(S1,f)==getfield(S2,f)[y] for f in datafields]))

# Randomly colored printing
printcol(r::Float64) = (r ≥ 1.00 ? 1 : r ≥ 0.75 ? 202 : r ≥ 0.50 ? 190 : r ≥ 0.25 ? 148 : 10)

# This function exists because Appveyor can't decide on its own write permissions
function safe_rm(file::String)
  try
    rm(file)
  catch err
    @warn(string("Can't remove ", file, ": throws error ", err))
  end
  return nothing
end

# test that each field has the right number of entries
function sizetest(S::GphysData, nt::Integer)
  @test ≈(S.n, nt)
  @test ≈(maximum([length(getfield(S,i)) for i in datafields]), nt)
  @test ≈(minimum([length(getfield(S,i)) for i in datafields]), nt)
  return nothing
end

# Basic sanity checks
function basic_checks(S::GphysData; allow_empty::Bool=false)
  # are we in scoping hell again?
  @test isempty(S) == false

  # does each field have the right number of entries?
  sizetest(S, S.n)

  for i = 1:S.n
    if allow_empty == false
      # were data read in?
      @test isempty(S.x[i]) == false

      # does :t indexing match :x?
      @test length(S.x[i]) == (S.fs[i] == 0.0 ? size(S.t[i],1) : length(S.x[i]))

    elseif isempty(S.x[i]) == false
      # does :t indexing match :x?
      @test length(S.x[i]) == (S.fs[i] == 0.0 ? size(S.t[i],1) : length(S.x[i]))
    else
      # an empty :x should never occur without an empty :t
      @test isempty(S.t[i])
    end
  end
  return nothing
end

# Get the start and end times of each channel in S
function get_edge_times(S::GphysData)
  ts = [S.t[i][1,2] for i=1:S.n]
  te = copy(ts)
  for i=1:S.n
    if S.fs[i] == 0.0
      te[i] = S.t[i][end,2]
    else
      te[i] = endtime(S.t[i], S.fs[i])
    end
  end
  return ts, te
end

# Convert lat, lon to x, y
function latlon2xy(xlat::Float64, xlon::Float64)
  s = sign(xlon)
  c = 111194.6976
  y = c*xlat
  d = acosd(cosd(xlon*s)*cosd(xlat))
  x = sqrt(c^2*d^2-y^2)
  return [round(Int32, s*x), round(Int32, y)]
end

# A simple time counting loop
function loop_time(ts::Int64, te::Int64; ti::Int64=86400000000)
  t1 = deepcopy(ts)
  j = 0
  while t1 < te
    j += 1
    t1 = min(ts + ti, te)
    s_str = int2tstr(ts + 1)
    t_str = int2tstr(t1)
    ts += ti
  end
  return j
end

# Remove low-gain seismic data channels
function remove_low_gain!(S::GphysData)
  i_low = findall([occursin(r".EL?", S.id[i]) for i=1:S.n])
  if !isempty(i_low)
    for k = length(i_low):-1:1
      @warn(join(["Low-gain, low-fs channel removed: ", S.id[i_low[k]]]))
      S -= S.id[i_low[k]]
    end
  end
  return nothing
end

# Test that data are time synched correctly within a SeisData structure
function sync_test!(S::GphysData)
  local L = [length(S.x[i])/S.fs[i] for i = 1:S.n]
  local t = [S.t[i][1,2] for i = 1:S.n]
  @test maximum(L) - minimum(L) ≤ maximum(2.0./S.fs)
  @test maximum(t) - minimum(t) ≤ maximum(2.0./S.fs)
  return nothing
end

function test_del_ranges(x_del::Array{UnitRange, 1}, ii::Array{Int64, 1})
  rr = [collect(i) for i in x_del]
  jj = vcat(rr...)
  @test jj == ii
end


text:
##############################
# Mapping functions on Dicts #
##############################
function mapvals!(f, d)
    for (k,v) in d
        d[k] = f(v)
    end
    return d
end
mapvals(f, d) = Dict(zip(collect(keys(d)), map(f, collect(values(d)))))

mapkeys(f, d) = Dict(zip(map(f, collect(keys(d))), collect(values(d))))
mapkv(f, d) = Dict(map(f, collect(d)))

#############
# Filtering #
#############
Base.filter!(f::Function, kt::Ket) = (filter!(f, dict(kt)); return kt)
Base.filter!(f::Function, br::Bra) = (filter!((k,v)->f(k,v'), dict(br)); return br)
Base.filter!(f::Function, op::OpSum) = (filter!(f, dict(op)); return op)
Base.filter!(f::Function, opc::DualOpSum) = (filter!((k,v)->f(k',v'), dict(opc)); return opc)

Base.filter(f::Function, kt::Ket) = similar(kt, filter(f, dict(kt)))
Base.filter(f::Function, br::Bra) = similar(br, filter((k,v)->f(k,v'), dict(br)))
Base.filter(f::Function, op::OpSum) = similar(op, filter(f, dict(op)))
Base.filter(f::Function, opc::DualOpSum) = similar(opc, filter((k,v)->f(k',v'), dict(opc)))
Base.filter(f::Function, op::OuterProduct) = filter(f, convert(OpSum, op))

########################
# mapcoeffs!/mapcoeffs #
########################
mapcoeffs!(f::Union(Function,DataType), k::Ket) = (mapvals!(f, dict(k)); return k)
mapcoeffs!(f::Union(Function,DataType), b::Bra) = (mapvals!(v->f(v')', dict(b)); return b)
mapcoeffs!(f::Union(Function,DataType), op::OpSum) = (mapvals!(f, dict(op)); return op)
mapcoeffs!(f::Union(Function,DataType), opc::DualOpSum) = (mapvals!(v->f(v')', dict(opc)); return opc)

mapcoeffs(f::Union(Function,DataType), kt::Ket) = similar(kt, mapvals(f, dict(kt)))
mapcoeffs(f::Union(Function,DataType), br::Bra) = similar(br, mapvals(v->f(v')', dict(br)))
mapcoeffs(f::Union(Function,DataType), op::OpSum) = similar(op, mapvals(f, dict(op)))
mapcoeffs(f::Union(Function,DataType), opc::DualOpSum) = similar(opc, mapvals(v->f(v')', dict(opc)))
mapcoeffs(f::Union(Function,DataType), op::OuterProduct) = mapcoeffs(f, convert(OpSum, op))

########################
# maplabels!/maplabels #
########################
maplabels(f::Union(Function,DataType), s::DiracState) = similar(s, mapkeys(f, dict(s)))
maplabels(f::Union(Function,DataType), op::OpSum) = similar(op, mapkeys(f, dict(op)))
maplabels(f::Union(Function,DataType), opc::DualOpSum) = OpSum(ptype(opc), mapkeys(label->f(label'), dict(op)))
maplabels(f::Union(Function,DataType), op::OuterProduct) = maplabels(f, convert(OpSum, op))

#######
# map #
#######
Base.map(f::Union(Function,DataType), obj::AbstractDirac) = similar(obj, mapkv(kv->f(kv[1], kv[2]), dict(obj)))
Base.map(f::Union(Function,DataType), br::Bra) = similar(br, mapkv(kv->br_tup(f(kv[1], kv[2]')), dict(d)))
Base.map(f::Union(Function,DataType), opc::DualOpSum) = OpSum(ptype(opc), mapkv(kv->f(kv[1]', kv[2]'), dict(d)))
Base.map(f::Union(Function,DataType), op::OuterProduct) = map(f, convert(OpSum, op))

br_tup(tup) = (tup[1], tup[2]')

export maplabels!,
    mapcoeffs!,
    maplabels,
    mapcoeffs


text:
struct BoundaryPaddedMesh
    mesh::Any
    gridsize::Any
    nfmside::Any
    numghostlayers::Any
    bottomghostcoords::Any
    rightghostcoords::Any
    topghostcoords::Any
    leftghostcoords::Any
end

function Base.show(io::IO, mesh::BoundaryPaddedMesh)
    dim = dimension(mesh)
    numghostlayers = number_of_ghost_layers(mesh)
    gridsize = grid_size(mesh)
    nfmside = nodes_per_mesh_side(mesh)

    str =
        "BoundaryPaddedMesh\n\tDimension : $dim\n\tGrid Size : $gridsize" *
        "\n\tNodes per side : $nfmside" *
        "\n\tNum. Ghost Layers : $numghostlayers"
    print(io, str)
end

function BoundaryPaddedMesh(mesh, numghostlayers)
    x0 = reference_corner(mesh)
    nfmside = nodes_per_mesh_side(mesh)
    meshwidths = mesh_widths(mesh)

    gridsize = meshwidths ./ (nfmside .- 1)

    bottomghostcoords = bottom_ghost_coordinates(
        x0,
        meshwidths[1],
        nfmside[1],
        gridsize[2],
        numghostlayers,
    )
    rightghostcoords = right_ghost_coordinates(
        x0,
        meshwidths,
        nfmside[2],
        gridsize[1],
        numghostlayers,
    )
    topghostcoords = top_ghost_coordinates(
        x0,
        meshwidths,
        nfmside[1],
        gridsize[2],
        numghostlayers,
    )
    leftghostcoords = left_ghost_coordinates(
        x0,
        meshwidths[2],
        nfmside[2],
        gridsize[1],
        numghostlayers,
    )

    return BoundaryPaddedMesh(
        mesh,
        gridsize,
        nfmside,
        numghostlayers,
        bottomghostcoords,
        rightghostcoords,
        topghostcoords,
        leftghostcoords,
    )
end

function background_mesh(mesh::BoundaryPaddedMesh)
    return mesh.mesh
end

function dimension(mesh::BoundaryPaddedMesh)
    return dimension(background_mesh(mesh))
end

function nodes_per_mesh_side(paddedmesh::BoundaryPaddedMesh)
    return paddedmesh.nfmside
end

function grid_size(mesh::BoundaryPaddedMesh)
    return mesh.gridsize
end

function number_of_ghost_layers(mesh::BoundaryPaddedMesh)
    return mesh.numghostlayers
end

function bottom_ghost_coordinates(mesh::BoundaryPaddedMesh)
    return mesh.bottomghostcoords
end

function right_ghost_coordinates(mesh::BoundaryPaddedMesh)
    return mesh.rightghostcoords
end

function top_ghost_coordinates(mesh::BoundaryPaddedMesh)
    return mesh.topghostcoords
end

function left_ghost_coordinates(mesh::BoundaryPaddedMesh)
    return mesh.leftghostcoords
end

function uniform_grid_coordinates(x0, widths, npts)
    xrange = range(x0[1], stop = x0[1] + widths[1], length = npts[1])
    yrange = range(x0[2], stop = x0[2] + widths[2], length = npts[2])

    xcoords = repeat(xrange, inner = npts[2])
    ycoords = repeat(yrange, outer = npts[1])

    return vcat(xcoords', ycoords')
end

function bottom_ghost_coordinates(x0, xwidth, npts, dy, numlayers)

    return uniform_grid_coordinates(
        x0 - [0.0, numlayers * dy],
        [xwidth, (numlayers - 1) * dy],
        [npts, numlayers],
    )
end

function right_ghost_coordinates(x0, widths, npts, dx, numlayers)

    return uniform_grid_coordinates(
        x0 + [widths[1] + dx, 0.0],
        [(numlayers - 1) * dx, widths[2]],
        [numlayers, npts],
    )
end

function top_ghost_coordinates(x0, widths, npts, dy, numlayers)

    return uniform_grid_coordinates(
        x0 + [0.0, widths[2] + dy],
        [widths[1], (numlayers - 1) * dy],
        [npts, numlayers],
    )
end

function left_ghost_coordinates(x0, ywidth, npts, dx, numlayers)

    return uniform_grid_coordinates(
        x0 - [numlayers * dx, 0.0],
        [(numlayers - 1) * dx, ywidth],
        [numlayers, npts],
    )
end

struct BoundaryPaddedLevelSet
    interiordist::Any
    bottomghostdist::Any
    rightghostdist::Any
    topghostdist::Any
    leftghostdist::Any
    paddedmesh::Any
end

function BoundaryPaddedLevelSet(
    paddedmesh::BoundaryPaddedMesh,
    spatialseedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius,
)

    bottomghostdist = distance_to_zero_levelset(
        bottom_ghost_coordinates(paddedmesh),
        spatialseedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )
    rightghostdist = distance_to_zero_levelset(
        right_ghost_coordinates(paddedmesh),
        spatialseedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )
    topghostdist = distance_to_zero_levelset(
        top_ghost_coordinates(paddedmesh),
        spatialseedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )
    leftghostdist = distance_to_zero_levelset(
        left_ghost_coordinates(paddedmesh),
        spatialseedpoints,
        seedcellids,
        levelset,
        tol,
        boundingradius,
    )

    return BoundaryPaddedLevelSet(
        coefficients(levelset),
        bottomghostdist,
        rightghostdist,
        topghostdist,
        leftghostdist,
        paddedmesh,
    )
end

function background_mesh(p::BoundaryPaddedLevelSet)
    return p.paddedmesh
end

function interior_values(paddedlevelset::BoundaryPaddedLevelSet)
    return paddedlevelset.interiordist
end

function left_ghost_values(paddedlevelset::BoundaryPaddedLevelSet)
    return paddedlevelset.leftghostdist
end

function right_ghost_values(paddedlevelset::BoundaryPaddedLevelSet)
    return paddedlevelset.rightghostdist
end

function bottom_ghost_values(paddedlevelset::BoundaryPaddedLevelSet)
    return paddedlevelset.bottomghostdist
end

function top_ghost_values(paddedlevelset::BoundaryPaddedLevelSet)
    return paddedlevelset.topghostdist
end

function first_order_horizontal_backward_difference(paddedlevelset)
    interior = interior_values(paddedlevelset)
    left = left_ghost_values(paddedlevelset)

    cols, rows = nodes_per_mesh_side(background_mesh(paddedlevelset))
    dx = grid_size(background_mesh(paddedlevelset))[1]

    npts = length(interior)
    derivative = zeros(npts)

    derivative[1:rows] = (interior[1:rows] - left) / dx

    for col = 2:cols
        prevstart = (col - 2) * rows + 1
        prevstop = prevstart + rows - 1

        start = (col - 1) * rows + 1
        stop = start + rows - 1

        derivative[start:stop] =
            (interior[start:stop] - interior[prevstart:prevstop]) / dx
    end

    return derivative
end

function first_order_horizontal_forward_difference(paddedlevelset)

    interior = interior_values(paddedlevelset)
    right = right_ghost_values(paddedlevelset)
    cols, rows = nodes_per_mesh_side(background_mesh(paddedlevelset))
    dx = grid_size(background_mesh(paddedlevelset))[1]

    npts = length(interior)
    derivative = zeros(npts)

    for col = 1:cols-1
        start = (col - 1) * rows + 1
        stop = start + rows - 1
        nextstart = stop + 1
        nextstop = nextstart + rows - 1

        derivative[start:stop] =
            (interior[nextstart:nextstop] - interior[start:stop]) / dx
    end

    start = (cols - 1) * rows + 1
    stop = start + rows - 1
    derivative[start:stop] = (right - interior[start:stop]) / dx

    return derivative
end

function first_order_vertical_backward_difference(paddedlevelset)

    interior = interior_values(paddedlevelset)
    bottom = bottom_ghost_values(paddedlevelset)
    cols, rows = nodes_per_mesh_side(background_mesh(paddedlevelset))
    dy = grid_size(background_mesh(paddedlevelset))[2]
    numghostlayers = number_of_ghost_layers(background_mesh(paddedlevelset))

    npts = length(interior)
    derivative = zeros(npts)

    interioridx = range(1, step = rows, length = cols)
    ghostidx = range(numghostlayers, step = numghostlayers, length = cols)

    derivative[interioridx] = (interior[interioridx] - bottom[ghostidx]) / dy

    for row = 2:rows
        interioridx = range(row, step = rows, length = cols)
        previnterioridx = range(row - 1, step = rows, length = cols)

        derivative[interioridx] =
            (interior[interioridx] - interior[previnterioridx]) / dy
    end

    return derivative
end

function first_order_vertical_forward_difference(paddedlevelset)
    interior = interior_values(paddedlevelset)
    top = top_ghost_values(paddedlevelset)
    cols, rows = nodes_per_mesh_side(background_mesh(paddedlevelset))
    dy = grid_size(background_mesh(paddedlevelset))[2]
    numghostlayers = number_of_ghost_layers(background_mesh(paddedlevelset))

    npts = length(interior)
    derivative = zeros(npts)

    for row = 1:rows-1
        interioridx = range(row, step = rows, length = cols)
        nextinterioridx = range(row + 1, step = rows, length = cols)

        derivative[interioridx] =
            (interior[nextinterioridx] - interior[interioridx]) / dy
    end

    interioridx = range(rows, step = rows, length = cols)
    ghostidx = range(1, step = numghostlayers, length = cols)

    derivative[interioridx] = (top[ghostidx] - interior[interioridx]) / dy

    return derivative
end

function first_order_nabla(u1, u2, v1, v2)
    return sqrt.(
        (max.(u1, 0.0)) .^ 2 +
        (min.(u2, 0.0)) .^ 2 +
        (max.(v1, 0.0)) .^ 2 +
        (min.(v2, 0.0)) .^ 2,
    )
end

function step_first_order_levelset(paddedlevelset, speed, dt)
    Dmx = first_order_horizontal_backward_difference(paddedlevelset)
    Dpx = first_order_horizontal_forward_difference(paddedlevelset)
    Dmy = first_order_vertical_backward_difference(paddedlevelset)
    Dpy = first_order_vertical_forward_difference(paddedlevelset)

    delplus = first_order_nabla(Dmx, Dpx, Dmy, Dpy)
    delminus = first_order_nabla(Dpx, Dmx, Dpy, Dmy)

    return paddedlevelset.interiordist -
           dt * (max.(speed, 0.0) .* delplus + min.(speed, 0.0) .* delminus)
end


text:
function compare_to(obj::Enum, arg0::Enum)
    return jcall(obj, "compareTo", jint, (Enum,), arg0)
end

function describe_constable(obj::Enum)
    return jcall(obj, "describeConstable", Optional, ())
end

function equals(obj::Enum, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_capitalized_abbreviation(obj::Month)
    return jcall(obj, "getCapitalizedAbbreviation", JString, ())
end

function get_capitalized_name(obj::Month)
    return jcall(obj, "getCapitalizedName", JString, ())
end

function get_declaring_class(obj::Enum)
    return jcall(obj, "getDeclaringClass", Class, ())
end

function get_lower_case_abbreviation(obj::Month)
    return jcall(obj, "getLowerCaseAbbreviation", JString, ())
end

function get_lower_case_name(obj::Month)
    return jcall(obj, "getLowerCaseName", JString, ())
end

function get_month(::Type{Month}, arg0::jint)
    return jcall(Month, "getMonth", Month, (jint,), arg0)
end

function get_number(obj::Month)
    return jcall(obj, "getNumber", jint, ())
end

function get_upper_case_abbreviation(obj::Month)
    return jcall(obj, "getUpperCaseAbbreviation", JString, ())
end

function get_upper_case_name(obj::Month)
    return jcall(obj, "getUpperCaseName", JString, ())
end

function hash_code(obj::Enum)
    return jcall(obj, "hashCode", jint, ())
end

function name(obj::Enum)
    return jcall(obj, "name", JString, ())
end

function ordinal(obj::Enum)
    return jcall(obj, "ordinal", jint, ())
end

function parse_month(::Type{Month}, arg0::JString)
    return jcall(Month, "parseMonth", Month, (JString,), arg0)
end

function to_string(obj::Enum)
    return jcall(obj, "toString", JString, ())
end

function value_of(::Type{Enum}, arg0::Class, arg1::JString)
    return jcall(Enum, "valueOf", Enum, (Class, JString), arg0, arg1)
end

function value_of(::Type{Month}, arg0::JString)
    return jcall(Month, "valueOf", Month, (JString,), arg0)
end

function values(::Type{Month})
    return jcall(Month, "values", Vector{Month}, ())
end


text:
""" Julia program to find the Length of Longest Increasing Subsequence

In this problem, given an array we have to find the length of the longest increasing subsequence that array can make.
The problem can be solved using Dynamic Programming
"""

function length_longest_increasing_subsequence(arr, n)
    max_len = 0
    # Initialize the dp array with the 1 as value, as the maximum length
    # at each point is atleast 1, by including that value in the sequence
    dp = ones(Int, n)

    """ Now Lets Fill the dp array in Bottom-Up manner
    Compare Each i'th element to its previous elements from 0 to i-1,
    If arr[i] > arr[j](where j = 0 to i-1), then it qualifies for increasing subsequence and
    If dp[i] < dp[j] + 1, then that subsequence  qualifies for being the longest one"""
    for i in 1:n
        for j in 1:(i-1)
            if(arr[i] > arr[j] && dp[i] < dp[j] + 1)
                dp[i] = dp[j] + 1
            end
        end
    end

    # Now Find the largest element in the dp array
    max_len = findmax(dp)[1]
    return max_len
end

print("What is the length of the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("No numbers present in the array!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
res = length_longest_increasing_subsequence(arr, n)
print("The length of the longest increasing subsequence of the given array is $res")

"""
Time Complexity - O(n^2), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
What is the length of the array? 5
Enter the numbers: 5 4 3 2 1
The length of the longest increasing subsequence of the given array is 1

"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
const DEFAULT_PRECISION = Ref{Clong}(256)

"""
    precision(<:Union{Arf,Arb,Acb})
Get the default precision (in bits) currently used for `Arblib` arithmetic.
"""
Base.precision(::Type{<:Union{Arf,Arb,Acb}}) = Int(DEFAULT_PRECISION[])

"""
    setprecision(::Type{<:Union{Arf, Arb, Acb}}, precision::Int)
Set the precision (in bits) to be used for `Arblib` arithmetic.
!!! warning
    This function is not thread-safe. It will affect code running on all threads, but
    its behavior is undefined if called concurrently with computations that use the
    setting.
"""
function Base.setprecision(::Type{<:Union{Arf,Arb,Acb}}, precision::Integer)
    if precision < 2
        throw(DomainError(precision, "`precision` cannot be less than 2."))
    end
    DEFAULT_PRECISION[] = precision
    return precision
end


text:
module NeuroEntities

export Entity

@enum Entity::UInt8 begin
    acquisition
    label
    space
    flip
    reconstruction
    subject
    sample
    processing
    part
    split
    description
    tracer
    recording
    resolution
    direction
    session
    task
    mtransfer
    run
    modality
    density
    inversion
    ceagent
    echo
end

function Base.String(x::Entity)
    if x === acquisition
        return "acq"
    elseif x === label
        return "label"
    elseif x === space
        return "space"
    elseif x === flip
        return "flip"
    elseif x === reconstruction
        return "rec"
    elseif x === subject
        return "sub"
    elseif x === sample
        return "sample"
    elseif x === processing
        return "proc"
    elseif x === part
        return "part"
    elseif x === split
        return "split"
    elseif x === description
        return "desc"
    elseif x === tracer
        return "trc"
    elseif x === recording
        return "recording"
    elseif x === resolution
        return "res"
    elseif x === direction
        return "dir"
    elseif x === session
        return "ses"
    elseif x === task
        return "task"
    elseif x === mtransfer
        return "mt"
    elseif x === run
        return "run"
    elseif x === modality
        return "mod"
    elseif x === density
        return "den"
    elseif x === inversion
        return "inv"
    elseif x === ceagent
        return "ce"
    elseif x === echo
        return "echo"
    end
end

function Entity(x::AbstractString)
    if x == "acq"
        return acquisition
    elseif x == "label"
        return label
    elseif x == "space"
        return space
    elseif x == "flip"
        return flip
    elseif x == "rec"
        return reconstruction
    elseif x == "sub"
        return subject
    elseif x == "sample"
        return sample
    elseif x == "proc"
        return processing
    elseif x == "part"
        return part
    elseif x == "split"
        return split
    elseif x == "desc"
        return description
    elseif x == "trc"
        return tracer
    elseif x == "recording"
        return recording
    elseif x == "res"
        return resolution
    elseif x == "dir"
        return direction
    elseif x == "ses"
        return session
    elseif x == "task"
        return task
    elseif x == "mt"
        return mtransfer
    elseif x == "run"
        return run
    elseif x == "mod"
        return modality
    elseif x == "den"
        return density
    elseif x == "inv"
        return inversion
    elseif x == "ce"
        return ceagent
    elseif x == "echo"
        return echo
    else
        error(x * " cannot be converted to an Entity type.")    end
end

function Entity(x::Symbol)
    if x == :acquisition
        return acquisition
    elseif x == :label
        return label
    elseif x == :space
        return space
    elseif x == :flip
        return flip
    elseif x == :reconstruction
        return reconstruction
    elseif x == :subject
        return subject
    elseif x == :sample
        return sample
    elseif x == :processing
        return processing
    elseif x == :part
        return part
    elseif x == :split
        return split
    elseif x == :description
        return description
    elseif x == :tracer
        return tracer
    elseif x == :recording
        return recording
    elseif x == :resolution
        return resolution
    elseif x == :direction
        return direction
    elseif x == :session
        return session
    elseif x == :task
        return task
    elseif x == :mtransfer
        return mtransfer
    elseif x == :run
        return run
    elseif x == :modality
        return modality
    elseif x == :density
        return density
    elseif x == :inversion
        return inversion
    elseif x == :ceagent
        return ceagent
    elseif x == :echo
        return echo
    else
        error("$(x) cannot be converted to an Entity type.")
    end
end

end


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

## dummy stub for https://github.com/JuliaBinaryWrappers/LibUnwind_jll.jl

module LibUnwind_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libunwind

# These get calculated in __init__()
PATH = Ref("")
LIBPATH = Ref("")
artifact_dir = ""
libunwind_handle = C_NULL
libunwind_path = ""

const libunwind = "libunwind.so.8"

function __init__()
    global artifact_dir = dirname(Sys.BINDIR)
    global LIBPATH[] = joinpath(Sys.BINDIR, Base.LIBDIR, "julia")
    # We only do something on Linux/FreeBSD
    @static if Sys.islinux() || Sys.isfreebsd()
        global libunwind_handle = dlopen(libunwind)
        global libunwind_path = dlpath(libunwind_handle)
    end
end

is_available() = @static (Sys.islinux() || Sys.isfreebsd()) ? true : false

end  # module LibUnwind_jll

max_stars_repo_name: saberdy/julia

max_issues_repo_name: saberdy/julia

max_forks_repo_name: saberdy/julia


text:
module OpenTelemetry

using Reexport

@reexport using OpenTelemetryAPI
@reexport using OpenTelemetrySDK
@reexport using OpenTelemetryExporterOtlpProtoGrpc
@reexport using OpenTelemetryExporterPrometheus
@reexport using OpenTelemetryProto

function __init__()
    if global_tracer_provider() isa OpenTelemetryAPI.DummyTracerProvider
        global_tracer_provider!(
            TracerProvider(
                span_processor = CompositSpanProcessor(
                    SimpleSpanProcessor(ConsoleExporter()),
                ),
            ),
        )
    end

    if global_meter_provider() isa OpenTelemetryAPI.DummyMeterProvider
        global_meter_provider!(MeterProvider())
    end
end

end


text:
using Turing
using Turing.RandomMeasures
using Distributions
using Random
using Libtask
using AutoGibbs

const data = [-1.48, -1.40, -1.16, -1.08, -1.02, 0.14, 0.51, 0.53, 0.78]

function update_histogram!(histogram, bin)
    if bin > length(histogram)
        push!(histogram, 1)
    else
        histogram[bin] += 1
    end

    return histogram
end

function histogram(d, z)
    nk = Int[]
    K = 0

    for n in 1:d.N
        nk = update_histogram!(nk, z[n])
        K = max(K, z[n])
    end

    return nk, K
end

mutable struct CRP <: Distributions.DiscreteMultivariateDistribution
    N::Int
    α::Float64
end

CRP(N, α, G₀) = CRP(N, α)

function Distributions.rand(rng::AbstractRNG, d::CRP)
    nk = Int[]
    D = ChineseRestaurantProcess(DirichletProcess(d.α), nk)
    z = rand(rng, D, d.N)
end

function Distributions.logpdf(d::CRP, z::AbstractVector{<:Int})
    nk, K = histogram(d, z)
    D = ChineseRestaurantProcess(DirichletProcess(d.α), nk)
    return mapreduce(zi -> logpdf(D, zi), +, z; init=0.0)
end

function (c::GibbsConditional{V, L})(θ) where {
    V<:VarName, L<:LogLikelihood{<:CRP}}

    Ω = c.base.dist.N

    conditionals = similar(Ωs, DiscreteNonParametric)

    for index in eachindex(Ωs, independent_distributions, conditionals)
        sub_vn = DynamicPPL.VarName(c.vn, (DynamicPPL.getindexing(c.vn)..., (index,)))
        θs_on_support = fixvalues(θ, sub_vn => Ωs[index])
        logtable = map(θs_on_support) do θ′
            c.base(θ′) + reduce(+, (β(θ′) for (vn, β) in c.blanket if vn == sub_vn), init=0.0)
        end
        conditionals[index] = DiscreteNonParametric(Ωs[index], softmax!(vec(logtable)))
    end

    return Product(conditionals)
end

@model function imm_oneshot(y, α, ::Type{T}=Array{Float64, 1}) where {T}
    N = length(y)

    G₀ = Normal()
    z ~ CRP(N, α, G₀)
    K = maximum(z)

    μ = T(undef, K)
    for k = 1:K
        μ[k] ~ G₀
    end

    for n = 1:N
        y[n] ~ Normal(μ[z[n]], 1.0)
    end
end

model_imm_oneshot = imm_oneshot(data[5:7], 10.0)

sample(model_imm_oneshot, Gibbs(AutoConditional(:z), HMC(0.01, 10, :μ)), 2)


text:
"""
Script for simulating a delta robot when accounting for the Base and platform radii
    - Lagrangian Formulation
        - position and velocity initial conditions
        - robot parameters (lengths and masses)
        - time span to simulate over
    - Hamiltonian Formulation
        - position and velocity (velocities are internally converted to momentas)
        - robot parameters (lengths and masses)
        - time span to simulate over
"""
using Revise
using Plots
using DifferentialEquations
using LinearAlgebra
using Sundials
using DASSL
using DeltaRobot

# Set initial conditions
u0, du0, diff_vars = delta_robot_ics(
    [pi/4; 3*pi/4; 0.0; pi/4; 3*pi/4; 0.0; pi/4; 3*pi/4; 0.0],
    [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0],
    [0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
)
@assert length(u0) == length(du0) # make sure position and velocity initial condition vectors have same dimensionality
@assert length(u0) == length(diff_vars) # make sure position and vector denoting which equations are differential equations have the same dimensionality
@assert u0[10:18] == du0[1:9] # make sure initial conditions are consistent
delta_robot_params = DeltaRobotParams(2.5, 3.0, 4.0, 0.312, 0.410, 9.8, 0.114, 0.051, 0.1, 0.1, 0.1, 100.0, 5.0); # Actual robot parameters with random d-params
tspan = (0.0, 20.0) # time span to simulate the system over
dt = 0.05 # time step

# Setup & solve the DAE under Lagrangian formalism
L_prob = DAEProblem(DeltaRobot.L_delta_robot_base!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars) # statement of the DAE problem
L_sol = solve(L_prob, IDA(), saveat=dt) # try :BCG

# Animate the solution to the differential equation
anim_L = DeltaRobot.animate_robot(L_sol, delta_robot_params)
gif(anim_L, "robot_animation.gif", fps= 1 / dt)

# Second set of initial conditions
phi0 = deg2rad.([20., 20., 40.0])
x0 = DeltaRobot.for_kin(phi0, [3*pi/4, 0.0, 3*pi/4, 0.0, 3*pi/4, 0.0], delta_robot_params)
u0, du0, diff_vars = delta_robot_ics(x0[1], zeros(9), zeros(9))
@assert length(u0) == length(du0) # make sure position and velocity initial condition vectors have same dimensionality
@assert length(u0) == length(diff_vars) # make sure position and vector denoting which equations are differential equations have the same dimensionality
@assert u0[10:18] == du0[1:9] # make sure initial conditions are consistent
DeltaRobot.holonomic_constraint(x0[1], delta_robot_params)

# Setup & solve the DAE under Lagrangian formalism
L_prob = DAEProblem(DeltaRobot.L_delta_robot_base!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars) # statement of the DAE problem
L_sol = solve(L_prob, IDA(), saveat=dt)

# Animate the solution to the differential equation
anim_L = DeltaRobot.animate_robot(L_sol, delta_robot_params)
gif(anim_L, "demo.gif", fps=1/dt)

# Gravity compensation
phi0 = deg2rad.([20., 20., 20.0])
x0 = DeltaRobot.for_kin(phi0, [3*pi/4, 0.0, 3*pi/4, 0.0, 3*pi/4, 0.0], delta_robot_params)
u0, du0, diff_vars = delta_robot_ics(x0[1], zeros(9), zeros(9))
@assert length(u0) == length(du0) # make sure position and velocity initial condition vectors have same dimensionality
@assert length(u0) == length(diff_vars) # make sure position and vector denoting which equations are differential equations have the same dimensionality
@assert u0[10:18] == du0[1:9] # make sure initial conditions are consistent
DeltaRobot.holonomic_constraint(x0[1], delta_robot_params)
grav_comp_prob = DAEProblem(DeltaRobot.gravity_compensated_model!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars) # statement of the DAE problem
grav_comp_sol = solve(grav_comp_prob,  IDA(linear_solver=:BCG), saveat=dt, reltol=1e-6,abstol=1e-6)

# Animate the solution to the differential equation
grav_comp_anim = DeltaRobot.animate_robot(grav_comp_sol, delta_robot_params)
gif(grav_comp_anim, fps=1/dt)

phi0 = deg2rad.([80., 80.0, 40.0])
x0 = DeltaRobot.for_kin(phi0, [3*pi/4, 0.0, 3*pi/4, 0.0, 3*pi/4, 0.0], delta_robot_params)
u0, du0, diff_vars = delta_robot_ics(x0[1], zeros(9), zeros(9))
@assert length(u0) == length(du0) # make sure position and velocity initial condition vectors have same dimensionality
@assert length(u0) == length(diff_vars) # make sure position and vector denoting which equations are differential equations have the same dimensionality
@assert u0[10:18] == du0[1:9] # make sure initial conditions are consistent
DeltaRobot.holonomic_constraint(x0[1], delta_robot_params)
grav_comp_prob = DAEProblem(DeltaRobot.gravity_compensated_model!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars) # statement of the DAE problem
grav_comp_sol = solve(grav_comp_prob,  IDA(linear_solver=:BCG), saveat=dt, reltol=1e-6,abstol=1e-6)

# Animate the solution to the differential equation
grav_comp_anim = DeltaRobot.animate_robot(grav_comp_sol, delta_robot_params)
gif(grav_comp_anim, fps=1/dt)

prob = DAEProblem(DeltaRobot.L_delta_robot_base_pdff!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars)
sol = solve(prob, IDA(linear_solver=:BCG), saveat=dt, reltol=1e-6, abstol=1e-6)
anim = DeltaRobot.animate_robot(sol, delta_robot_params)
gif(anim, fps=1/dt)

# # PD 1
# PD_prob = DAEProblem(DeltaRobot.gravity_compensated_model_PD!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars)
# # IDA(), IDA(linear_solver=:BCG), DABDF2(), DImplicitEuler(), dassl()
# PD_sol = solve(PD_prob, IDA(linear_solver=:GMRES), saveat=dt, reltol=1e-2*2, abstol=1e-2*2)
# PD_anim = DeltaRobot.animate_robot(PD_sol, delta_robot_params)
# gif(PD_anim, fps=1/dt)
#
# # PD
# phi0 = deg2rad.([0., 0., 0.])
# x0 = DeltaRobot.for_kin(phi0, [3*pi/4, 0.0, 3*pi/4, 0.0, 3*pi/4, 0.0], delta_robot_params)
# u0, du0, diff_vars = delta_robot_ics(x0[1], zeros(9), zeros(9))
# @assert length(u0) == length(du0) # make sure position and velocity initial condition vectors have same dimensionality
# @assert length(u0) == length(diff_vars) # make sure position and vector denoting which equations are differential equations have the same dimensionality
# @assert u0[10:18] == du0[1:9] # make sure initial conditions are consistent
# DeltaRobot.holonomic_constraint(x0[1], delta_robot_params)
# PD_prob = DAEProblem(DeltaRobot.gravity_compensated_model_PD!, du0, u0, tspan, delta_robot_params, differential_vars=diff_vars)
# # IDA(), IDA(linear_solver=:BCG), DABDF2(), DImplicitEuler(), dassl()
# PD_sol = solve(PD_prob,  IDA(linear_solver=:BCG), saveat=dt, reltol=1e-4, abstol=1e-4)
# PD_anim = DeltaRobot.animate_robot(PD_sol, delta_robot_params)
# gif(PD_anim, fps=1/dt)

max_stars_repo_name: zborffs/Delta

max_issues_repo_name: zborffs/Delta

max_forks_repo_name: zborffs/Delta


text:
using PSCOPF

using Test
using Dates

@testset "test_check_branch" begin

    @testset "definition_example" begin
        branch = PSCOPF.Networks.Branch("branch_1_2", 1000.)
        @test PSCOPF.check(branch)
    end

    @testset "limit>0" begin
        branch = PSCOPF.Networks.Branch("branch_1_2", -1000.)
        @test !PSCOPF.check(branch)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
p = Path()
@test length(p) == 0
@test start(p) == 0.0
@test finish(p) == 0.0
@test duration(p) == 0.0

p = Path([Jump(0.5,"a","c"),Jump(1.5,"c","b"),Jump(2.2,"b","d")])
@test length(p)==3
@test start(p) == 0.5
@test finish(p) == 2.2
@test duration(p)==2.2-0.5
@test is_valid(p)

p2 = Path([Jump(0.5,"a","c"),Jump(1.5,"d","b"),Jump(2.2,"b","d")])
@test !is_valid(p2)

dp = DurationPath([DurationJump(0.5,"a","c",0.6),DurationJump(1.5,"c","b",0.5),DurationJump(2.2,"b","d",0.3)])
@test length(dp) == 3
@test start(dp) == 0.5
@test finish(dp) == 2.5
@test duration(dp) == 2.0
@test is_valid(dp)

p2 = Path([Jump(2.5,"d","c"),Jump(3.1,"c","e")])
@test is_valid(p+p2)
@test length(p+p2) == 5
@test start(p+p2) == 0.5
@test finish(p+p2) == 3.1
@test duration(p+p2) == 3.1-0.5
#@test Path([Jump(2.2,"b","d"),Jump(2.5,"d","c")]) ⊆ (p+p2)

dp2 = DurationPath([DurationJump(2.5,"d","e",1.6),DurationJump(4.5,"e","b",0.01)])
@test is_valid(dp+dp2)
@test length(dp+dp2) == 5
@test start(dp+dp2) == 0.5
@test finish(dp+dp2) == 4.51
@test duration(dp+dp2) == 4.51-0.5
#@test DurationPath([DurationJump(1.5,"c","b",0.5),DurationJump(2.2,"b","d",0.3),DurationJump(2.5,"d","e",1.6)]) ⊆ (dp+dp2)


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Script for the baseline model
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

using DataFrames, RData, LinearAlgebra, Statistics, XLSX

include("transform_WIOD_2016_3.jl") # Script with functions to import and transform raw data
include("price_hat.jl") # Script with function to obtain the price index
include("wage_hat.jl") # Script with function to obtain the wages and gross output

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

years = 2000:2014 # vector with years covered by WIOD Rev. 2016
N = 41 # number of countries
S = 35 # number of industries
dir = "C:/Users/u0148308/Desktop/raw/" # location of raw data

# -------------------------------------------------------------------------------------------------------------------------------------------------------------

function baseline(dir::String, year::Integer, N::Integer, S::Integer)

    # iteration parameters
    vfactor = 0.4
    tolerance = 1e-6
    max_iteration = 100
    iteration = 0
    max_error = 1e7

    # sectoral trade elasticity
    θ = 5 # assumption
    θ = fill(θ, N*S) # NS×1, work with country-industry elasticities

    # trade costs
    τ_hat_Z = ones(N*S, N*S) # NS×NS
    τ_hat_F = ones(N*S, N) # NS×N

    # initialize wages and price indices
    w_hat = ones(N) # N×1

    # ------------

    Z, F, Y, F_ctry, TB_ctry, VA_ctry, VA_coeff, γ, α, π_Z, π_F = transform_WIOD_2016(dir, "excel", year, N, S)

    #TB_ctry .= 0.0 # (if active => adjustments such that there is no trade deficit, i.e. counterfactual in itself)

    # ------------

    while max_error > tolerance && iteration <= max_iteration

        P_hat_Z, P_hat_F, π_hat_Z, π_hat_F, cost_hat = create_price_index_hat(w_hat, τ_hat_Z, τ_hat_F, VA_coeff, γ, π_Z, π_F, θ)

        # compute counterfactual trade shares
        π_hat_Z = ifelse.(isinf.(π_hat_Z), 0.0, π_hat_Z) # remove Inf
        π_hat_F = ifelse.(isinf.(π_hat_F), 0.0, π_hat_F)

        global π_prime_Z = π_Z .* π_hat_Z # NS×NS
        global π_prime_F = π_F .* π_hat_F # NS×N

        # store last wage in case new optimization obtains negative wages
        w_hat_prev = copy(w_hat)

        w_hat, Y_prime, ETB_ctry = create_wages_hat(w_hat, vfactor, π_prime_Z, π_prime_F, VA_ctry, TB_ctry, γ, α)

        # update iteration parameters
        error = abs.(w_hat .- w_hat_prev)
        max_error = maximum(error) # update error
        iteration += 1 # update iteration count

        if minimum(w_hat) < 0.0
            w_hat[:] = copy(w_hat_prev)
            println("Iteration $iteration completed with error $max_error (wage negative, rerun with previous estimate)")
        else
            println("Iteration $iteration completed with error $max_error")
        end

    end

    # ------------

    return w_hat, Y_prime, ETB_ctry, π_prime_Z, π_prime_F, α, γ

end

w_hat, Y_prime, ETB_ctry, π_prime_Z, π_prime_F, α, γ = baseline(dir, 2011, N, S)


text:
module BaseModel
export AbstractModel

abstract type AbstractModel end

end


text:
### MCRange

abstract type MCRange{T<:Integer} <: Range{T} end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
This is the reference implmentation of the BLAKE3 algorithm done in Julia.  No attempt has
been made to improve the performance it's primarly purpose is to be easy to read and used
to debug/verify faster implmentations.
"""
module RefImpl

using Base.Iterators

const OUT_LEN = UInt(32)
const KEY_LEN = UInt(32)
const BLOCK_LEN = UInt(64)
const CHUNK_LEN = UInt(1024)

const CHUNK_START = UInt32(1) << 0
const CHUNK_END   = UInt32(1) << 1
const PARENT      = UInt32(1) << 2
const ROOT        = UInt32(1) << 3
const KEYED_HASH  = UInt32(1) << 4
const DERIVE_KEY_CONTEXT  = UInt32(1) << 5
const DERIVE_KEY_MATERIAL = UInt32(1) << 6

const IV = [
    UInt32(0x6A09E667), UInt32(0xBB67AE85), UInt32(0x3C6EF372), UInt32(0xA54FF53A),
    UInt32(0x510E527F), UInt32(0x9B05688C), UInt32(0x1F83D9AB), UInt32(0x5BE0CD19)
]

const MSG_PERMUTATION = [ 2, 6, 3, 10, 7, 0, 4, 13, 1, 11, 12, 5, 9, 14, 15, 8 ]

function g(state::AbstractVector{UInt32}, a::Int, b::Int, c::Int, d::Int, mx::UInt32, my::UInt32)::Nothing
    state[a] = state[a] + state[b] + mx;
    state[d] = bitrotate(state[d] ⊻ state[a], -16)
    state[c] = state[c] + state[d]
    state[b] = bitrotate(state[b] ⊻ state[c], -12)
    state[a] = state[a] + state[b] + my
    state[d] = bitrotate(state[d] ⊻ state[a], -8)
    state[c] = state[c] + state[d]
    state[b] = bitrotate(state[b] ⊻ state[c], -7)
    nothing
end

function round(state::AbstractVector{UInt32}, m::AbstractVector{UInt32})::Nothing
    # Mix the columns
    g(state, 1, 5,  9, 13, m[1], m[2])
    g(state, 2, 6, 10, 14, m[3], m[4])
    g(state, 3, 7, 11, 15, m[5], m[6])
    g(state, 4, 8, 12, 16, m[7], m[8])
    # Mix the diagonals
    g(state, 1, 6, 11, 16, m[9],  m[10])
    g(state, 2, 7, 12, 13, m[11], m[12])
    g(state, 3, 8,  9, 14, m[13], m[14])
    g(state, 4, 5, 10, 15, m[15], m[16])
    nothing
end

function permute(m::AbstractVector{UInt32})::Nothing
    local permuted = Vector{UInt32}(undef, 16)

    for i in 1:16
        permuted[i] = m[MSG_PERMUTATION[i]+1]
    end
    copyto!(m, permuted)
    nothing
end

function compress(
            chaining_value::AbstractVector{UInt32},
            block_words::AbstractVector{UInt32},
            counter::UInt64,
            block_len::UInt32,
            flags::UInt32
        )::AbstractVector{UInt32}
    local state = [
        chaining_value[1], chaining_value[2], chaining_value[3], chaining_value[4],
        chaining_value[5], chaining_value[6], chaining_value[7], chaining_value[8],
        IV[1], IV[2], IV[3], IV[4],
        UInt32(counter & 0xffffffff), UInt32(counter >> 32),
        block_len, flags
    ]
    local block = copy(block_words)

    round(state, block) # round 1
    permute(block)
    round(state, block) # round 2
    permute(block)
    round(state, block) # round 3
    permute(block)
    round(state, block) # round 4
    permute(block)
    round(state, block) # round 5
    permute(block)
    round(state, block) # round 6
    permute(block)
    round(state, block) # round 7

    for i in 1:8
        state[i] ⊻= state[i + 8]
        state[i + 8] ⊻= chaining_value[i]
    end

    return state
end

function first_8_words(compression_output::AbstractVector{UInt32})::Vector{UInt32}
    return [
        compression_output[1], compression_output[2],
        compression_output[3], compression_output[4],
        compression_output[5], compression_output[6],
        compression_output[7], compression_output[8]
    ]
end

function words_from_little_endian_bytes(bytes::AbstractVector{UInt8}, words::Vector{UInt32})
    for (byte_block, i) in zip(partition(bytes, 4), eachindex(words))
        words[i] = ltoh(reinterpret(UInt32, byte_block)[1])
    end
end

struct Output
    input_chaining_value::Vector{UInt32}
    block_words::Vector{UInt32}
    counter::UInt64
    block_len::UInt32
    flags::UInt32
end

function chaining_value(self::Output)::Vector{UInt32}
    return first_8_words(compress(
        self.input_chaining_value,
        self.block_words,
        self.counter,
        self.block_len,
        self.flags
    ))
end

function root_output_bytes(self::Output, out_slice::AbstractVector{UInt8})::Nothing
    local output_block_counter = UInt64(0)

    for out_block in partition(out_slice, 2 * OUT_LEN)
        local words = compress(
            self.input_chaining_value,
            self.block_words,
            output_block_counter,
            self.block_len,
            self.flags | ROOT
        )
        for (word, out_word) in zip(words, partition(out_block, 4))
            out_word[1] = UInt8(word & 0xff)
            if length(out_word) >= 2
                out_word[2] = UInt8((word >> 8) & 0xff)
            end
            if length(out_word) >= 3
                 out_word[3] = UInt8((word >> 16) & 0xff)
            end
            if length(out_word) >= 4
                out_word[4] = UInt8((word >> 24) & 0xff)
            end
        end
        output_block_counter += 1
    end
    nothing
end

mutable struct ChunkState
    chaining_value::Vector{UInt32}
    chunk_counter::UInt64
    block::Vector{UInt8}
    block_len::UInt8
    blocks_compressed::UInt8
    flags::UInt32
end

function ChunkState(key::AbstractVector{UInt32}, chunk_counter::UInt64, flags::UInt32)::ChunkState
    return ChunkState(key, chunk_counter, zeros(UInt8, BLOCK_LEN), 0, 0, flags)
end

function len(self::ChunkState)::UInt
    return BLOCK_LEN * UInt(self.blocks_compressed) + UInt(self.block_len)
end

function start_flag(self::ChunkState)::UInt32
    if self.blocks_compressed == 0
        return CHUNK_START
    end

    return 0
end

function update(self::ChunkState, input::AbstractVector{UInt8})::Nothing
    local data = view(input, 1:length(input))

    while isempty(data) == false
        # If the block buffer is full, compress it and clear it. More
        # input is coming, so this compression is not CHUNK_END.

        if UInt(self.block_len) == BLOCK_LEN
            local block_words = zeros(UInt32, 16)
            words_from_little_endian_bytes(self.block, block_words)
            self.chaining_value = first_8_words(compress(
                self.chaining_value,
                block_words,
                self.chunk_counter,
                UInt32(BLOCK_LEN),
                self.flags | start_flag(self)
            ))
            self.blocks_compressed += UInt8(1)
            self.block = zeros(UInt8, BLOCK_LEN)
            self.block_len = 0
        end

        local want = BLOCK_LEN - UInt(self.block_len)
        local take = min(want, length(data))
        copyto!(self.block, self.block_len+1, data, 1, take)
        self.block_len += UInt8(take)
        data = view(data, (take+1):length(data))
    end
end

function output(self::ChunkState)::Output
    local block_words = zeros(UInt32, 16)

    words_from_little_endian_bytes(self.block, block_words)

    return Output(
        self.chaining_value,
        block_words,
        self.chunk_counter,
        self.block_len,
        self.flags | start_flag(self) | CHUNK_END
    )
end

function parent_output(
            left_child_cv::AbstractVector{UInt32},
            right_child_cv::AbstractVector{UInt32},
            key::AbstractVector{UInt32},
            flags::UInt32
        )::Output
    local block_words = vcat(left_child_cv, right_child_cv)

    return Output(key, block_words, 0, UInt32(BLOCK_LEN), PARENT | flags)
end

function parent_cv(
            left_child_cv::AbstractVector{UInt32},
            right_child_cv::AbstractVector{UInt32},
            key::AbstractVector{UInt32},
            flags::UInt32
        )::Vector{UInt32}
    return chaining_value(parent_output(left_child_cv, right_child_cv, key, flags))
end

mutable struct Blake3Ctx
    chunk_state::ChunkState
    key::Vector{UInt32}
    cv_stack::Array{UInt32, 2}
    cv_stack_len::UInt8
    flags::UInt32
    Blake3Ctx(key::Vector{UInt32}=IV, flags::UInt32=UInt32(0)) = new(
        ChunkState(key, UInt64(0), flags), key, Array{UInt32}(undef, 8, 54), 0, flags
    )
    Blake3Ctx(key::Vector{UInt8}) = begin
        local key_words = zeros(UInt32, 8)
        words_from_little_endian_bytes(key, key_words)
        return Blake3Ctx(key_words, KEYED_HASH)
    end
    Blake3Ctx(content::AbstractString) = begin
        local context_hasher = Blake3Ctx(IV, DERIVE_KEY_CONTEXT)
        update!(context_hasher, Vector{UInt8}(content))
        local context_key = zeros(UInt8, KEY_LEN)
        digest(context_hasher, context_key)
        local context_key_words = zeros(UInt32, 8)
        words_from_little_endian_bytes(context_key, context_key_words)
        return Blake3Ctx(context_key_words, DERIVE_KEY_MATERIAL)
    end
end

function push_stack(self::Blake3Ctx, cv::AbstractVector{UInt32})::Nothing
    self.cv_stack_len += 1
    self.cv_stack[:, self.cv_stack_len] .= cv
    nothing
end

function pop_stack(self::Blake3Ctx)::Vector{UInt32}
    self.cv_stack_len -= 1
    return self.cv_stack[:, self.cv_stack_len+1]
end

# Section 5.1.2 of the BLAKE3 spec explains this algorithm in more detail.
function add_chunk_chaining_value(
            self::Blake3Ctx,
            new_cv::AbstractVector{UInt32},
            total_chunks::UInt64
        )::Nothing
    # This chunk might complete some subtrees. For each completed subtree,
    # its left child will be the current top entry in the CV stack, and
    # its right child will be the current value of `new_cv`. Pop each left
    # child off the stack, merge it with `new_cv`, and overwrite `new_cv`
    # with the result. After all these merges, push the final value of
    # `new_cv` onto the stack. The number of completed subtrees is given
    # by the number of trailing 0-bits in the new total number of chunks.
    while total_chunks & 1 == 0
        new_cv = parent_cv(pop_stack(self), new_cv, self.key, self.flags)
        total_chunks = total_chunks >> 1;
    end
    push_stack(self, new_cv)
    nothing
end

"""
    update!(self::Blake3Ctx, input::AbstractVector{UInt8}, datalen=length(input))::Nothing

Hash more data.  This method can be called multiple times until all the data has been
hashed.
"""
function update!(self::Blake3Ctx, input::AbstractVector{UInt8}, datalen=length(input))::Nothing
    local data = view(input, 1:datalen)

    while isempty(data) == false
        # If the current chunk is complete, finalize it and reset the
        # chunk state. More input is coming, so this chunk is not ROOT.
        if len(self.chunk_state) == CHUNK_LEN
            local chunk_cv     = chaining_value(output(self.chunk_state))
            local total_chunks = self.chunk_state.chunk_counter + 1
            add_chunk_chaining_value(self, chunk_cv, total_chunks)
            self.chunk_state = ChunkState(self.key, total_chunks, self.flags)
        end

        # Compress input bytes into the current chunk state.
        local want = CHUNK_LEN - len(self.chunk_state)
        local take = min(want, length(data))
        update(self.chunk_state, view(data, 1:take))
        data = view(data, (take+1):length(data))
    end
end

"""
    digest(self::Blake3Ctx, out_slice::AbstractArray{UInt8})::Nothing

Generate the hash and write it to the provided array.  The contents of out_slice will be
overwritten.

__Details:__

Normally you would pass in an array of 32 bytes which the method will fill with the hash.
However if you want more (or less) bytes you can pass in a an array of a different size.
"""
function digest(self::Blake3Ctx, out_slice::AbstractArray{UInt8})::Nothing
    # Starting with the Output from the current chunk, compute all the
    # parent chaining values along the right edge of the tree, until we
    # have the root Output.
    local out                    = output(self.chunk_state)
    local parent_nodes_remaining = UInt(self.cv_stack_len)

    while parent_nodes_remaining > 0
        parent_nodes_remaining -= 1
        out = parent_output(
            self.cv_stack[:, parent_nodes_remaining+1],
            chaining_value(out),
            self.key,
            self.flags
        )
    end
    root_output_bytes(out, out_slice)
end

"""
    digest(self::Blake3Ctx, bytes=32)::Vector{UInt8}

Generate the hash and return it.

__Details:__

Normally you would want 32 bytes back (256bits) but you can have it return more or less
bytes based on your need.
"""
function digest(self::Blake3Ctx, bytes=32)::Vector{UInt8}
    local output = Vector{UInt8}(undef, bytes)
    digest(self, output)
    return output
end

end


text:
using OmegaCore
using Test

function ptest1()
  x = 1 ~ Normal(0, 1)
  y = 2 ~ Normal(0, 1)
  z(ω) = x(ω) + y(ω)
  xy(ω) = (x(ω), y(ω))
  xy |ᶜ z ==ₚ z_

  # Defines pararametric representatio nover x and y
  proposal(::typeof(z), z_, θ) =
    (x => z_ - θ, y => θ)

  # defines distribution over x and y
  function proposal(::typeof(z), z_, ω)
    θ = 3 ~ Normal(0, 1)(ω)
    (x => z_ - θ, y => θ)
  end

end

function ptest2()
  β = 1 ~ Beta(0.5)
  x = 2 ~ Bernoulli(β)
  function proposal(::typeof(x), x_, ω)
    if x_
      r = uniform(0, p)
  end

end

def bernoulli(p):
  u = sample(rand(), "u")
  return u < p
with proposal(b):
  if b:
    r = sample(uniform(0, p), "r")
  else:
    r = sample(uniform(p, 1), "r")
  return {"u": r}
with proposal(t):
  return {"r": t["u"]}

  # Issues:

# How will this be implemented?
Conceptually what we want to do is something parametric relational programming, where we take the values we know, e.g. the conditioned outputs,
and use those to construct proposals on the rest.
The difficulties are how to fit this into a forward simulation framework, which Julia is

In posthook of funciton, we evaluate the proposal
-- e.g. after computing z(w), we do the proposal to update values of x and y (WHERE?)
This seems a bit late, to compute z we have to have already computed x and y.  We could overwrite it of course.

How about in prehook, z depends on x and y so in some sense contains x and y, so __before__ we evaluate z(w) we do the proposal for x and y

# How we will determine "who" chooses value for "who"

# What if no pinv given

# is typeof(z) specific enough, i.e. different randvars can have same type?

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# right hand side function in our system
function f!(du, u, p, t)
    du[1] = u[1] * (1.0-u[2]);
    du[2] = u[2] * (u[1]-1.0);
    du
end

# set up ODE
function solve_ode(u0, Tmax)
    # define the ODE problem
    prob = ODEProblem(f!, u0, (0., Tmax))
    # solve the problem
    sol = solve(prob)

    return sol
end


text:
struct TypeListItem{TThis, TNext} end
struct EmptyTypeList end
const TypeList = Union{EmptyTypeList, TypeListItem}

encode() = EmptyTypeList
encode(t) = TypeListItem{t, EmptyTypeList}
encode(t, ts...) = TypeListItem{t, encode(ts...)}

decode(t::Tuple) = t
decode(::Type{EmptyTypeList}) = ()
decode(::Type{TypeListItem{TThis, EmptyTypeList}}) where TThis = (TThis,)
decode(::Type{TypeListItem{TThis, TNext}}) where {TThis, TNext} = (TThis, decode(TNext)...)

Base.length(::Type{EmptyTypeList}) = 0
function Base.length(::Type{TypeListItem{TThis, TNext}}) where {TThis, TNext}
    return 1 + length(TNext)
end

Base.findfirst(::Type, ::Type{EmptyTypeList}) = nothing
function Base.findfirst(t::Type, ::Type{TypeListItem{TThis, TNext}}) where {TThis, TNext}
    t == TThis && return 1
    tailidx = findfirst(t, TNext)
    return isnothing(tailidx) ? nothing : tailidx + 1
end

pretty(::Type{TypeListItem{TThis, TNext}}) where {TThis, TNext} = begin
    return "$(TThis), " * pretty(TNext)
end

pretty(::Type{EmptyTypeList}) = ""


text:
module ComradeGalactic

using Comrade
using Reexport

@reexport using GalacticOptim

function GalacticOptim.OptimizationFunction(post::Comrade.Posterior, args...; kwargs...)
    throw("Transform the posterior first using `asflat` or `ascube`")
end

function GalacticOptim.OptimizationFunction(post::Comrade.TransformedPosterior, args...; kwargs...)
    ℓ(x,p) = -logdensityof(post, x)
    return OptimizationFunction(ℓ, args...; kwargs...)
end

# """
#     solve(prob, opt, args...; transform=nothing, kwargs...)
# Find the optimum of the problem.

# Optional keyword argument `transform``
# will transform the solution to parameter space, and the return will be
# a tuple with the galactic optim solution in the first element and the optimum
# location in model space in the second argument
# """
# function GalacticOptim.solve(prob::OptimizationProblem, opt, transform::Union{Nothing, HypercubeTransform.AbstractHypercubeTransform, HypercubeTransform.TransformVariables.AbstractTransform, HypercubeTransform.NamedFlatTransform}, args...; kwargs...)
#     sol = solve(prob, opt, args...; kwargs...)
#     if isnothing(transform)
#         return sol
#     else
#         p = HypercubeTransform.transform(transform, sol.u)
#         return sol, p
#     end
# end

end


text:
@noinline function simulate_unit!(simprim::SimulationPrimitives, i, do_r=true)

    sim = SimulationDrawsVector(simprim, i)
    simtmp = SimulationTmp(simprim)
    theta_d, theta_r, theta_q = split_thetas(simprim)
    unit, obs_r, wells = getindex.(data(simprim), i)

    qm = _qm(sim)
    fill!(qm,0)

    M = _num_sim(sim)

    # royalty posterior
    if do_r
        simloglik_royalty!(first(obs_r), theta_r, sim, false)
        logsumexp!(qm)
    else
        invM = 1/M
        qm .= invM
    end

    model = _model(DataDrill(unit))
    solve_vf_and_update_itp!(model, theta_d, ichars(unit), false)
    @assert any( data(EVobj(value_function(model))) .!= 0 )

    # do update
    for m = OneTo(M)
        simm = getindex(sim, m)

        if num_initial_leases(unit) > 0
            for lease in InitialDrilling(unit)
                weight = qm[m] * j1chars(lease)
                simulate_lease!(simprim, lease, simm, weight)
            end
        else
            DD = DevelopmentDrilling(unit)
            for lease in DD
                weight = qm[m]
                simulate_lease!(simprim, lease, simm, weight)
            end
        end

    end
end

function simulate_unit!(i, do_r=true)
     simprim = get_g_SimulationPrimitives()
     simulate_unit!(simprim, i, do_r)
end


text:
srcillum_helper(field_ginsu::Array{T}, nthreads) where {T<:Float32} =
    ccall((:srcillum_helper_float,libillumination), Cvoid,
        (Ptr{Float32},Csize_t,Cint), field_ginsu, length(field_ginsu), nthreads)

srcillum_helper(field_ginsu::Array{T}, nthreads) where {T<:Float64} =
    ccall((:srcillum_helper_double,libillumination), Cvoid,
        (Ptr{Float64},Csize_t,Cint), field_ginsu, length(field_ginsu), nthreads)

illum_accumulate(field_ginsu_accum::Array{T}, field_ginsu::Array{T}, nthreads) where {T<:Float32} =
    ccall((:illum_accumulate_float,libillumination), Cvoid,
        (Ptr{Float32},Ptr{Float32},Csize_t,Cint), field_ginsu_accum, field_ginsu, length(field_ginsu), nthreads)

illum_accumulate(field_ginsu_accum::Array{T}, field_ginsu::Array{T}, nthreads) where {T<:Float64} =
    ccall((:illum_accumulate_double,libillumination), Cvoid,
        (Ptr{Float64},Ptr{Float64},Csize_t,Cint), field_ginsu_accum, field_ginsu, length(field_ginsu), nthreads)

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD