text:
import Base.cross, Base.dot, Base.-, Base.+, Base.*, Base./

type Vector
	x::Real
	y::Real
	z::Real
end

# Outer Constructors

Vector() = Vector(0, 0, 0)

# Methods on Vector

function toList(a::Vector)
	return [a.x, a.y, a.z]
end

function toVec(a::Array)
	return Vector(a[1], a[2], a[3])
end

function +(a::Vector, b::Vector)
	return toVec(toList(a) .+ toList(b))
end

function -(a::Vector, b::Vector)
	return toVec(toList(a) .- toList(b))
end

function *(a::Vector, b::Real)
	return toVec(toList(a) .* b)
end

function magnitude(a::Vector)
	return sqrt(a.x ^ 2 + a.y ^ 2 + a.z ^ 2)
end

function normal(a::Vector)
	mag = magnitude(a)
	return toVec(toList(a) ./ mag)
end

function dot(a::Vector, b::Vector)
	return dot(toList(a), toList(b))
end

function cross(a::Vector, b::Vector)
	return toVec(cross(toList(a), toList(b)))
end


text:
# ------------------------------------------------------------------------------------------
# # Julia 中的线性代数
# > Based on work by Andreas Noack Jensen (MIT) (http://www.econ.ku.dk/phdstudent/noack/)
# > with edits from Jane Herriman
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 基础线性代数操作
# 定义一个随机矩阵
# ------------------------------------------------------------------------------------------

A = rand(1:4,3,3)

# ------------------------------------------------------------------------------------------
# 定义一个元素全为 1 的向量
# ------------------------------------------------------------------------------------------

x = fill(1.0, (3,)) # = fill(1.0, 3)

# ------------------------------------------------------------------------------------------
# 注意 $A$ 的类型为 `Array{Int64,2}`，而 $x$ 的类型为 `Array{Float64,1}`。
#
# Julia 定义向量 `Vector{Type}` 的别名为 `Array{Type,1}`，矩阵`Matrix{Type}` 的别名为 `Array{Type,2}`。
#
# 许多线性代数的基础操作和其他语言一样
#
# ### 乘法
# ------------------------------------------------------------------------------------------

b = A*x

# ------------------------------------------------------------------------------------------
# ### 转置
# 就像在其他语言中 `A'` 表示对 `A` 进行共轭转置
# ------------------------------------------------------------------------------------------

A'

# ------------------------------------------------------------------------------------------
# 我们还可以通过 `transpose` 函数获得转置矩阵
# ------------------------------------------------------------------------------------------

transpose(A)

# ------------------------------------------------------------------------------------------
# ### 转置的乘法
# Julia 中某些情况下可以省略 `*` 号
# ------------------------------------------------------------------------------------------

A'A

# ------------------------------------------------------------------------------------------
# ### 解线性方程组
# 用方阵 $A$ 表示的线性方程组 $Ax=b$ 可以用左除运算符（函数）`\` 求解
# ------------------------------------------------------------------------------------------

A\b

# ------------------------------------------------------------------------------------------
# 超定线性方程组（overdetermined linear system）的情况下 `A\b` 给出*最小二乘解（least squares solution）*
# ------------------------------------------------------------------------------------------

Atall = rand(3, 2)

Atall\b

# ------------------------------------------------------------------------------------------
# 缺秩最小二乘问题（a rank-deficient least squares problem）情况下，会给出*极小范数最小二乘解（minimum norm least
# squares solution）*
# ------------------------------------------------------------------------------------------

v = rand(3)
rankdef = hcat(v, v)

rankdef\b

# ------------------------------------------------------------------------------------------
# 欠定问题（an underdetermined solution）下，会给出 最小范数解（minimum norm solution）
# ------------------------------------------------------------------------------------------

bshort = rand(2)
Ashort = rand(2, 3)

Ashort\bshort

# ------------------------------------------------------------------------------------------
# ## `LinearAlgebra` 包
#
# 如上所示，Julia默认自带了很多线性代数操作，而标准库 `LinearAlgebra` 提供了更多的相关名词和函数。
# 特别是它提供了因子分解和一些结构化矩阵类型。
# 和其他包一样，可以通过 `using LinearAlgebra` 引入这些额外的特征。
# ------------------------------------------------------------------------------------------

using LinearAlgebra

# ------------------------------------------------------------------------------------------
# ## 练习
#
# ### 10.1
# 计算向量 `v` 和它自己的内积，并赋值给 `dot_v`。
# ------------------------------------------------------------------------------------------

v = [1,2,3]

# 在此作答
dot_v = v

@assert dot_v == 14

# ------------------------------------------------------------------------------------------
# ### 10.2
# 计算向量 v 和它自己的外积(Outer product)，并赋值给 `outer_v`。
#
# 提示：外积是一种克罗内克积(Kronecker product)
#
# 计算叉积(cross product)，赋值给 `cross_v`
# ------------------------------------------------------------------------------------------

# 在此作答
outer_v = v
cross_v = v

@assert outer_v == [1 2 3
                    2 4 6
                    3 6 9]

@assert cross_v == [0, 0, 0]

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

ext: jl

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
"""
"""
type CellLocation{T<:Real}
  x::Array{T,1}
  y::Array{T,1}
  z::Array{T,1}
end

"""
"""
type CellSize{T<:Real}
  x::Array{T,1}
  y::Array{T,1}
  z::Array{T,1}
end

"""
"""
type FaceLocation{T<:Real}
  x::Array{T,1}
  y::Array{T,1}
  z::Array{T,1}
end

"""
"""
type MeshStructure
  dimension::Real
  dims::Array{Int,1}
  cellsize::CellSize
  cellcenters::CellLocation
  facecenters::FaceLocation
  corner::Array{Int,1}
  edge::Array{Int,1}
end

"""
"""
type CellValue{T<:Real}
  domain::MeshStructure
  value::Array{T}
end

"""
"""
type CellVector{T<:Real}
  domain::MeshStructure
  xvalue::Array{T}
  yvalue::Array{T}
  zvalue::Array{T}
end

"""
"""
type FaceValue{T<:Real}
  domain::MeshStructure
  xvalue::Array{T}
  yvalue::Array{T}
  zvalue::Array{T}
end

"""
"""
type BorderValue{T<:Real}
  a::Array{T}
  b::Array{T}
  c::Array{T}
  periodic::Bool
end

"""
"""
type BoundaryCondition
  domain::MeshStructure
  left::BorderValue
  right::BorderValue
  bottom::BorderValue
  top::BorderValue
  back::BorderValue
  front::BorderValue
end


text:
euclidex = quote
    while a != b
        if a > b
            a = a - b
            p = p - q
            r = r - s
        else
            q = q - p
            b = b - a
            s = s - r
        end
    end
end

fermat = quote
    while r != 0
        if r > 0
            r = r - v
            v = v + 2
        else
            r = r + u
            u = u + 2
        end
    end
end

wensley = quote
    while d>= E
        if P < a+b
            b = b/2
            d = d/2
        else
            a = a+b
            y = y+d/2
            b = b/2
            d = d/2
        end
    end
end

lcm = quote
    while x != y
        if x > y
            x = x - y
            v = v + u
        else
            y = y - x
            u = u + v
        end
    end
end

knuth = quote
    while (s >= d) && (r != 0)
        if 2*r-rp+q < 0
        t  = r
        r  = 2*r-rp+q+d+2
        rp = t
        q  = q+4
        d  = d+2
        elseif (2*r-rp+q >= 0) && (2*r-rp+q < d+2)
        t  = r
        r  = 2*r-rp+q
        rp = t
        d  = d+2
        elseif (2*r-rp+q >= 0) && (2*r-rp+q >= d+2) && (2*r-rp+q < 2*d+4)
        t  = r
        r  = 2*r-rp+q-d-2
        rp = t
        q  = q-4
        d  = d+2
        else # ((2*r-rp+q >= 0) && (2*r-rp+q >= 2*d+4))
        t  = r
        r  = 2*r-rp+q-2*d-4
        rp = t
        q  = q-8
        d  = d+2
        end
    end
end

divbin = quote
    while b != B
        x = 2*x
        b = b/2
        if r >= b
            x = x+1
            r = r-b
        end
    end
end

psolv2m = quote
    n = 1
    while true
        if true
            a = 2*(n+1)*(n+3/2)*a
            b = 4*(n+1)*b
            c = 1/2*(n+3/2)*c
            n = n+1
        else
            a = 2*a
            b = 4*b
            c = 1/2*c
        end
    end
end

psolv3m = quote
    n1, n2 = 1, 1
    while true
        if true
            a = 2*(n1+1)*(n1+3/2)*a
            b = 4*(n1+1)*b
            c = 1/2*(n1+3/2)*c
            n1 = n1+1
        elseif c>0
            a = 2*a
            b = 4*b
            c = 1/2*c
        else
            a = 2*(n2+1)^3*(n2+3/2)^3*a
            b = 4*(n2+1)^3*b
            c = 1/2*(n2+3/2)^3*c
            n2 = n2+1
        end
    end
end

psolv4m = quote
    n1, n2 = 1, 1
    while true
        if true
            a = 2*(n1+1)*(n1+3/2)*a
            b = 4*(n1+1)*b
            c = 1/2*(n1+3/2)*c
            n1 = n1+1
        elseif c>0
            a = 2*a
            b = 4*b
            c = 1/2*c
        elseif c<0
            a = 2*(n2+1)^3*(n2+3/2)^3*a
            b = 4*(n2+1)^3*b
            c = 1/2*(n2+3/2)^3*c
            n2 = n2+2
        else
            a = 6*a
            b = 36*b
            c = 1/6*c
        end
    end
end

psolv10m = quote
    n1, n2, n3, n4, n5 = 1, 1, 1, 1, 1
    while true
        if true
            a = 2*(n1+1)*(n1+3/2)*a
            b = 4*(n1+1)*b
            c = 1/2*(n1+3/2)*c
            n1 = n1+1
        elseif c>0
            a = 2*a
            b = 4*b
            c = 1/2*c
        elseif c<0
            a = 2*(n2+1)^3*(n2+3/2)^3*a
            b = 4*(n2+1)^3*b
            c = 1/2*(n2+3/2)^3*c
            n2 = n2+1
        elseif c<0
            a = 6*a
            b = 36*b
            c = 1/6*c
        elseif c<0
            a = 2*(n3+1)^5*(n3+3/2)^5*a
            b = 4*(n3+1)^5*b
            c = 1/2*(n3+3/2)^5*c
            n3 = n3+1
        elseif c<0
            a = 3*a
            b = 9*b
            c = 1/3*c
        elseif c<0
            a = 2*(n4+1)^7*(n4+3/2)^7*a
            b = 4*(n4+1)^7*b
            c = 1/2*(n4+3/2)^7*c
            n4 = n4+1
        elseif c<0
            a = 5*a
            b = 25*b
            c = 1/5*c
        elseif c<0
            a = 2*(n5+1)^9*(n5+3/2)^9*a
            b = 4*(n5+1)^9*b
            c = 1/2*(n5+3/2)^9*c
            n5 = n5+1
        else
            a = 4*a
            b = 16*b
            c = 1/4*c
        end
    end
end

prodbin = quote
    while y!=0
        if y % 2 ==1
            z = z+x;
            y = y-1;
        end
        x = 2*x;
        y = y/2;
    end
end

dijkstra = quote
    while q!=1
        q=q/4;
        h=p+q;
        p=p/2;
        if r>=h
            p=p+q;
            r=r-h;
        end
    end
end

z3sqrt = quote
    while 2*p*r >= e
        if 2*r-2*q-p >= 0
            r = 2*r-2*q-p;
            q = q+p;
            p = p/2;
        else
            r = 2*r;
            p = p/2;
        end
    end
end

writers = quote
    r = 0;
    w = 0;
    k = 14;
    # c1 = 3;
    # c2 = 2;
    while true
        global r, w, k, c1, c2
        if w == 0
            r = r+1;
            k = k-c1;
        elseif r == 0
            w = w+1;
            k = k-c2;
        elseif w==0
            r = r-1;
            k = k+c1;
        elseif r==0
            w = w-1;
            k = k+c2;
        end
        @info  2*w+k+3*r-14==0
    end
end


text:
using LinearAlgebra
@testset "BBCS Chebyshev exponential" begin

    # Basic implementation
    for (k, θ) in enumerate([1.38e-5, 2.92e-3, 0.1295, 0.636, 2.212])
        Aθ = θ * I * 1im + [0 1e-8; -1e-8 0]
        (graph, cref) = graph_bbcs_cheb_exp(k, T = Complex{BigFloat})
        tv = range(-θ, θ, length = 41)
        @test all(
            abs.(eval_graph(graph, big.(tv)) - exp.(-1im * big.(tv))) .<
            eps() / 2,
        )
        (graph, cref) = graph_bbcs_cheb_exp(k, T = ComplexF64)
        @test eval_graph(graph, Aθ) ≈ exp(-1im * Aθ)
        @test sum(values(graph.operations) .== :mult) == k
    end

    # Scalign-and-squaring phase
    for (k, θ) in
        enumerate([2 * 2.212, 4 * 2.212, 8 * 2.212, 16 * 2.212, 32 * 2.212])
        kk = k + 5
        Aθ = θ * I * 1im + [0 1e-8; -1e-8 0]
        (graph, cref) = graph_bbcs_cheb_exp(kk, T = Complex{BigFloat})
        tv = range(-θ, θ, length = 41)
        @test all(
            abs.(eval_graph(graph, big.(tv)) - exp.(-1im * big.(tv))) .<
            eps() * 2^(k - 1),
        )
        (graph, cref) = graph_bbcs_cheb_exp(kk, T = ComplexF64)
        @test eval_graph(graph, Aθ) ≈ exp(-1im * Aθ)
        @test sum(values(graph.operations) .== :mult) == kk
    end
end


text:
using Test

@testset "TESTSETS: Test read of TCX string" begin
    @testset "CASE: Test string is not XML" begin
        err, _ = TCX.parse_tcx_str("This is not a TCX string\nNot even a little")
        @test err == TCX.CLIENT_ERROR
    end

    @testset "CASE: Test a TCX records run." begin
        path_to_file = joinpath(@__DIR__, tcx_sample_file)
        err, data = TCX.parse_tcx_str(read(path_to_file, String))
        @test (err == TCX.OK) & (getActivityType(data) == "Running")
    end

    @testset "CASE: Test string is not TCX" begin
        err, _ = @test_logs(
            (:warn, "Invalid TCX string: <NotTCX>Oh no! Blorp!</NotTCX>"),
            TCX.parse_tcx_str("<NotTCX>Oh no! Blorp!</NotTCX>")
        )
        @test err == TCX.CLIENT_TCX_ERROR
    end
end


text:
# Testing:
#
# - consistency of first-passage time densities

using DiffModels
using Test

function testpdfequal(d::AbstractDrift, b::AbstractBounds, tmax,
                      g1ref::Vector{Float64}, g2ref::Vector{Float64}, tol)
    g1, g2 = pdf(d, b, tmax)
    @test length(g1) == length(g2)
    @test length(g1) == length(g1ref)
    @test length(g2) == length(g2ref)
    for n = 1:length(g1)
        @test g1[n] ≈ g1ref[n] atol=tol
        @test g2[n] ≈ g2ref[n] atol=tol
    end
end

# symmetric bounds

mu = 1.2
bound = 1.1
t = 0.5
dt = 0.005
tmax = 3.0
maxn = length(0:dt:tmax)

@test (pdful(ConstDrift(mu, dt), ConstSymBounds(bound, dt), t)[1] ≈
       pdfu(ConstDrift(mu, dt), ConstSymBounds(bound, dt), t))
@test (pdful(ConstDrift(mu, dt), ConstSymBounds(bound, dt), t)[2] ≈
       pdfl(ConstDrift(mu, dt), ConstSymBounds(bound, dt), t))
@test (pdful(ConstDrift(mu, dt), ConstAsymBounds(bound, -bound, dt), t)[1] ≈
       pdfu(ConstDrift(mu, dt), ConstAsymBounds(bound, -bound, dt), t))
@test (pdful(ConstDrift(mu, dt), ConstAsymBounds(bound, -bound, dt), t)[2] ≈
       pdfl(ConstDrift(mu, dt), ConstAsymBounds(bound, -bound, dt), t))

g1, g2 = pdf(VarDrift(fill(mu, maxn), dt),
             VarAsymBounds(VarBound(fill(bound, maxn), dt),
                           VarBound(fill(bound, maxn), dt)), tmax)
testpdfequal(ConstDrift(mu, dt), ConstSymBounds(bound, dt), tmax, g1, g2, 1e-4)
testpdfequal(ConstDrift(mu, dt), ConstAsymBounds(bound, -bound, dt), tmax, g1, g2, 1e-4)
testpdfequal(ConstDrift(mu, dt), LinearSymBounds(bound, 0.0, dt), tmax, g1, g2, 1e-4)
testpdfequal(ConstDrift(mu, dt),
             AsymBounds{LinearBound, LinearBound}(LinearBound(bound, 0.0, dt),
                                                  LinearBound(bound, 0.0, dt)),
             tmax, g1, g2, 1e-4)

# asymmetric bounds

mu = 1.2
bounds = [1.1, 0.5]
bound_slopes = [-0.02, -0.01]
t = 0.5
dt = 0.005
tmax = 3.0
maxn = length(0:dt:tmax)

@test (pdful(ConstDrift(mu, dt), ConstAsymBounds(bounds[1], -bounds[2], dt), t)[1] ≈
       pdfu(ConstDrift(mu, dt), ConstAsymBounds(bounds[1], -bounds[2], dt), t))
@test (pdful(ConstDrift(mu, dt), ConstAsymBounds(bounds[1], -bounds[2], dt), t)[2] ≈
       pdfl(ConstDrift(mu, dt), ConstAsymBounds(bounds[1], -bounds[2], dt), t))

g1, g2 = pdf(VarDrift(fill(mu, maxn), dt),
             VarAsymBounds(VarBound(fill(bounds[1], maxn), dt),
                           VarBound(fill(bounds[2], maxn), dt)), tmax)
testpdfequal(ConstDrift(mu, dt), ConstAsymBounds(bounds[1], -bounds[2], dt), tmax, g1, g2, 1e-4)
testpdfequal(ConstDrift(mu, dt),
             AsymBounds{LinearBound, LinearBound}(LinearBound(bounds[1], 0.0, dt),
                                                  LinearBound(bounds[2], 0.0, dt)),
             tmax, g1, g2, 1e-4)

g1, g2 = pdf(VarDrift(fill(mu, maxn), dt),
             VarAsymBounds(VarBound(bounds[1] .+ bound_slopes[1] * dt * (0:maxn-1), dt),
                           VarBound(bounds[2] .+ bound_slopes[2] * dt * (0:maxn-1), dt)),
             tmax)
testpdfequal(ConstDrift(mu, dt),
             AsymBounds{LinearBound, LinearBound}(LinearBound(bounds[1], bound_slopes[1], dt),
                                                  LinearBound(bounds[2], bound_slopes[2], dt)),
             tmax, g1, g2, 1e-4)


text:
const RegSolveAlg = Cuint

const REG_SOLVE_FGMRES = RegSolveAlg(0)
const REG_SOLVE_LGMRES = RegSolveAlg(1)

struct RegSolveCtrl{T<:ElFloatType}
    alg::RegSolveAlg
    relTol::T
    relTolRefine::T
    maxIts::ElInt
    maxRefineIts::ElInt
    restart::ElInt
    progress::ElBool
    time::ElBool
end

function RegSolveCtrl(::Type{T};
                    alg = REG_SOLVE_FGMRES,
                    relTol = eps(T)^0.5,
                    relTolRefine = eps(T)^0.8,
                    maxIts = 4,
                    maxRefineIts = 2,
                    restart = 4,
                    progress = false,
                    time = false) where {T<:ElFloatType}

    RegSolveCtrl{T}(RegSolveAlg(alg),
        T(relTol),
        T(relTolRefine),
        ElInt(maxIts),
        ElInt(maxRefineIts),
        ElInt(restart),
        ElBool(progress),
        ElBool(time))
end

mutable struct Permutation
    obj::Ptr{Cvoid}
end

function destroy(P::Permutation)
    ElError(ccall(("ElPermutationDestroy", libEl), Cuint,
        (Ptr{Cvoid},),
        P.obj))
    return nothing
end

function Permutation()
    obj = Ref{Ptr{Cvoid}}(0)
    ElError(ccall(("ElPermutationCreate", libEl), Cuint,
        (Ref{Ptr{Cvoid}},),
        obj))
    P = Permutation(obj[])
    finalizer(destroy, P)
    return P
end

mutable struct DistPermutation
    obj::Ptr{Cvoid}
end

function destroy(P::DistPermutation)
    ElError(ccall(("ElDistPermutationDestroy", libEl), Cuint,
        (Ptr{Cvoid},),
        P.obj))
    return nothing
end

function DistPermutation(grid::Grid = DefaultGrid[])
    obj = Ref{Ptr{Cvoid}}(0)
    ElError(ccall(("ElDistPermutationCreate", libEl), Cuint,
        (Ref{Ptr{Cvoid}}, Ptr{Cvoid}),
        obj, grid.obj))
    P = DistPermutation(obj[])
    finalizer(destroy, P)
    return P
end

for mattype in ("", "Dist")
    mat = Symbol(mattype, "Matrix")
    _p  = Symbol(mattype, "Permutation")

    for (elty, ext) in ((:Float32, :s),
                        (:Float64, :d),
                        (:ComplexF32, :c),
                        (:ComplexF64, :z))
        @eval begin

            function _cholesky!(uplo::UpperOrLower, A::$mat{$elty})
                ElError(ccall(($(string("ElCholesky", mattype, "_", ext)), libEl), Cuint,
                    (UpperOrLower, Ptr{Cvoid}),
                    uplo, A.obj))
                return A
            end

            function _lu!(A::$mat{$elty})
                p = $_p()
                ElError(ccall(($(string("ElLUPartialPiv", mattype, "_", ext)), libEl), Cuint,
                    (Ptr{Cvoid}, Ptr{Cvoid}),
                    A.obj, p.obj))
                return A, p
            end
        end
    end
end


text:
"""
The string part of JSCode.
"""
struct JSString
    source::String
end

"""
Javascript code that supports interpolation of Julia Objects.
Construction of JSCode via string macro:
```julia
jsc = js"console.log(\$(some_julia_variable))"
```
This will decompose into:
```julia
jsc.source == [JSString("console.log("), some_julia_variable, JSString("\"")]
```
"""
struct JSCode
    source::Vector{Union{JSString, Any}}
end

"""
Represent an asset stored at an URL.
We try to always have online & local files for assets
If one gives an online resource, it will be downloaded, to host it locally.
"""
struct Asset
    media_type::Symbol
    # We try to always have online & local files for assets
    # If you only give an online resource, we will download it
    # to also be able to host it locally
    online_path::String
    local_path::String
    onload::Union{Nothing, JSCode}
end

"""
Encapsulates frontend dependencies. Can be used in the following way:

```Julia
const noUiSlider = Dependency(
    :noUiSlider,
    # js & css dependencies are supported
    [
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.js",
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.css"
    ]
)
# use the dependency on the frontend:
evaljs(session, js"\$(noUiSlider).some_function(...)")
```
jsrender will make sure that all dependencies get loaded.
"""
struct Dependency
    name::Symbol # The JS Module name that will get loaded
    assets::Vector{Asset}
end

"""
A web session with a user
"""
struct Session
    # indicates, whether session is in fuse mode, so it won't
    # send any messages, until fuse ends and send them all at once
    fusing::Base.RefValue{Bool}
    connections::Vector{WebSocket}
    observables::Dict{String, Tuple{Bool, Observable}} # Bool -> if already registered with Frontend
    message_queue::Vector{Dict{Symbol, Any}}
    dependencies::Set{Asset}
    on_document_load::Vector{JSCode}
end

struct Routes
    table::Vector{Pair{Any, Any}}
end

function Routes(pairs::Pair...)
    return Routes([pairs...])
end

pattern_priority(x::Pair) = pattern_priority(x[1])
pattern_priority(x::String) = 1
pattern_priority(x::Tuple) = 2
pattern_priority(x::Regex) = 3

function Base.setindex!(routes::Routes, f, pattern)
    idx = findfirst(pair-> pair[1] == pattern, routes.table)
    if idx !== nothing
        routes.table[idx] = pattern => f
    else
        push!(routes.table, pattern => f)
    end
    # Sort for priority so that exact string matches come first
    sort!(routes.table, by = pattern_priority)
    return
end

function apply_handler(f, args...)
    f(args...)
end

function apply_handler(chain::Tuple, context, args...)
    f = first(chain)
    result = f(args...)
    apply_handler(Base.tail(chain), context, result...)
end

function delegate(routes::Routes, application, request::Request, args...)
    for (pattern, f) in routes.table
        match = match_request(pattern, request)
        if match !== nothing
            context = (
                routes = routes,
                application = application,
                request = request,
                match = match
            )
            return apply_handler(f, context, args...)
        end
    end
    # If no route is found we have a classic case of 404!
    # What a classic this response!
    return response_404("Didn't find route for $(request.target)")
end

function match_request(pattern::String, request)
    return request.target == pattern ? pattern : nothing
end

function match_request(pattern::Regex, request)
    return match(pattern, request.target)
end

function match_request(pattern::Tuple, request)
    return Matcha.matchat(pattern, request.target)
end

"""
The application one serves
"""
struct Application
    url::String
    port::Int
    sessions::Dict{String, Dict{String, Session}}
    server_task::Ref{Task}
    server_connection::Ref{TCPServer}
    routes::Routes
    websocket_routes::Routes
end

WebSockets.getrawstream(io::IO) = io

function websocket_request()
    headers = [
        "Host" => "127.0.0.1",
        "User-Agent" => "Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:68.0) Gecko/20100101 Firefox/68.0",
        "Accept" => "*/*",
        "Accept-Encoding" => "gzip, deflate, br",
        "Accept-Language" => "de,en-US;q=0.7,en;q=0.3",
        "Cache-Control" => "no-cache",
        "Connection" => "keep-alive, Upgrade",
        "Dnt" => "1",
        "Origin" => "https://localhost",
        "Pragma" => "no-cache",
        "Sec-Websocket-Extensions" => "permessage-deflate",
        "Sec-Websocket-Key" => "BL3d8I8KC5faPjubRM0riA==",
        "Sec-Websocket-Version" => "13",
        "Upgrade" => "websocket",
    ]
    msg = HTTP.Request(
        "GET",
        "/",
        headers,
        UInt8[],
        parent = nothing,
        version = v"1.1.0"
    )
    return Stream(msg, IOBuffer())
end

"""
warmup(application::Application)

Warms up the application, by sending a couple of request.
"""
function warmup(application::Application)
    yield() # yield to server task to give it a chance to get started
    task = application.server_task[]
    if Base.istaskdone(task)
        error("Webserver doesn't serve! Error: $(fetch(task))")
    end
    # Make a websocket request
    stream = websocket_request()
    try
        stream_handler(application, stream)
    catch e
        # TODO make it not error so we can test this properly
        # This will error, since its not a propper websocket request
        @debug "Error in stream_handler" exception=e
    end

    bundle_url = url(JSCallLibLocal) # online URL
    target = AssetRegistry.register(JSCallLibLocal.local_path) # http target part
    request = Request("GET", AssetRegistry.register(JSCallLibLocal.local_path))

    delegate(application.routes, application, request)

    wait_time = 10.0; start = time() # wait for max 10 s
    success = false
    while time() - start < wait_time
        # Block as long as our server doesn't actually serve the bundle
        if Base.istaskdone(task)
            error("Webserver doesn't serve! Error: $(fetch(task))")
            break
        end
        # async + fetch allows the server to do work while the request gets sent!
        resp = fetch(@async WebSockets.HTTP.get(bundle_url))
        if resp.status == 200
            success = true
            break
        end
        sleep(0.1)
    end
    if !success
        error("Waited $(wait_time)s for webserver to start up, and didn't receive any response")
    end
end

function stream_handler(application::Application, stream::Stream)
    if WebSockets.is_upgrade(stream)
        WebSockets.upgrade(stream) do request, websocket
            delegate(
                application.websocket_routes, application, request, websocket
            )
        end
        return
    end
    f = HTTP.RequestHandlerFunction() do request
        delegate(
            application.routes, application, request,
        )
    end
    HTTP.handle(f, stream)
end

const MATCH_HEX = r"[\da-f]"
const MATCH_UUID4 = MATCH_HEX^8 * r"-" * (MATCH_HEX^4 * r"-")^3 * MATCH_HEX^12
const MATCH_SESSION_ID = MATCH_UUID4 * r"/" * MATCH_HEX^4

function serve_dom(context, dom)
    application = context.application
    session_id = string(uuid4())
    session = Session()
    application.sessions[session_id] = Dict("base" => session)
    return html(dom2html(session, session_id, dom(session, context.request)))
end

"""
Application(
        dom, url::String, port::Int;
        verbose = false
    )

Creates an application that manages the global server state!
"""
function Application(
        dom, url::String, port::Int;
        verbose = false,
        routes = Routes(
            "/" => ctx-> serve_dom(ctx, dom),
            r"/assetserver/" * MATCH_HEX^40 * r"-.*" => file_server,
            r".*" => (context)-> response_404()
        ),
        websocket_routes = Routes(
            r"/" * MATCH_SESSION_ID => websocket_handler
        )
    )

    application = Application(
        url, port, Dict{String, Dict{String, Session}}(),
        Ref{Task}(), Ref{TCPServer}(),
        routes,
        websocket_routes
    )
    start(application)
    # warmup server!
    warmup(application)

    return application
end

function isrunning(application::Application)
    return (isassigned(application.server_task) &&
        isassigned(application.server_connection) &&
        !istaskdone(application.server_task[]) &&
        isopen(application.server_connection[]))
end

function Base.close(application::Application)
    # Closing the io connection should shut down the HTTP listen loop
    for (id, clients) in application.sessions
        for (id, session) in clients
            close(session)
        end
    end

    close(application.server_connection[])
    # For good measures, wait until the task finishes!
    try
        wait(application.server_task[])
    catch e
        # the wait will throw with the below exception
        if !(e isa TaskFailedException && e.task.exception.code == -4079)
            rethrow(e)
        end
    end
end

function start(application::Application)
    isrunning(application) && return
    address = Sockets.InetAddr(
        parse(Sockets.IPAddr, application.url), application.port
    )
    ioserver = Sockets.listen(address)
    application.server_connection[] = ioserver
    # pass tcp connection to listen, so that we can close the server
    application.server_task[] = @async HTTP.listen(
            application.url, application.port, server = ioserver
        ) do stream::Stream
        Base.invokelatest(stream_handler, application, stream)
    end
    return
end

function route!(application::Application, pattern_f::Pair)
    application.routes[pattern_f[1]] = pattern_f[2]
end

function route!(f, application::Application, pattern)
    route!(application, pattern => f)
end

function websocket_route!(application::Application, pattern_f::Pair)
    application.webscoket_routes[pattern_f[1]] = pattern_f[2]
end


text:
function setindex_multithreads_worker( channel::Channel{Tuple}, buf::Array{T,N}, ba::BigArray{D,T} ) where {D,T,N}
    C = get_encoding(ba)
    mipLevelName = get_mip_level_name(ba)
    for (blockID, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in channel
        # println("global range of chunk: $(cartesian_range2string(chunkGlobalRange))")
        chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer = adjust_volume_boundary(ba, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer)
        delay = 0.05
        for t in 1:4
            try
                chk = buf[rangeInBuffer]
                key = joinpath(mipLevelName, cartesian_range2string( chunkGlobalRange ) )
                ba.kvStore[ key ] = encode( chk, C)
                # for the GSDicts backend, the haskey function will really download the object!
                #@assert haskey(ba.kvStore, key)
                break
            catch err
                println("catch an error while saving in BigArray: $err")
                @show typeof(err)
                @show stacktrace()
                if t==4
                    println("rethrow the error: $err")
                    rethrow()
                else
                    @warn("retry for the $(t)'s time.")
                end
                sleep(delay*(0.8+(0.4*rand())))
                delay *= 10
                println("retry for the $(t)'s time: $(string(chunkGlobalRange))")
            end
        end
    end
end

"""
    put array in RAM to a BigArray backend
this version uses channel to control the number of asynchronized request
"""
function setindex_multithreads!( ba::BigArray{D,T}, buf::Array{T,N},
                       idxes::Union{UnitRange, Int, Colon} ... ) where {D,T,N}
    idxes = colon2unit_range(buf, idxes)
    offset = get_offset(ba)
    chunkSize = BigArrays.get_chunk_size(ba)
    # check alignment
    @assert all(map((x,y,z)->mod(first(x) - 1 - y, z),
                    idxes, offset.I, chunkSize).==0)
                    "the start of index should align with BigArray chunk size"
    t1 = time()
    baIter = ChunkIterator(idxes, chunkSize; offset=offset)
    @sync begin
        channel = Channel{Tuple}( CHUNK_CHANNEL_SIZE )
        @async begin
            for iter in baIter
                put!(channel, iter)
            end
            close(channel)
        end
        for i in 1:TASK_NUM
            @async setindex_multithreads_worker(channel, buf, ba)
        end
    end
    elapsed = time() - t1 # sec
    println("saving speed: $(sizeof(buf)/1024/1024/elapsed) MB/s")
end

function getindex_multithreads_worker!(chan::Channel{Tuple}, buf::Array{T,N}, ba::BigArray{D,T}) where {D,T,N}
    baRange = CartesianIndices(ba)
    C = get_encoding(ba)
    mipLevelName = get_mip_level_name(ba)
    for (blockId, chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer) in chan
        if any(map((x,y)->x>y, first(globalRange).I, last(baRange).I)) ||
            any(map((x,y)->x<y, last(globalRange).I, first(baRange).I))
            @warn("out of volume range, keep it as zeros")
            continue
        end
        chunkGlobalRange, globalRange, rangeInChunk, rangeInBuffer =
            adjust_volume_boundary(ba, chunkGlobalRange, globalRange,
                                   rangeInChunk, rangeInBuffer)
        chunkSize = (last(chunkGlobalRange) - first(chunkGlobalRange) +
                     one(CartesianIndex{N})).I
        try
            #println("global range of chunk: $(cartesian_range2string(chunkGlobalRange))")
            key = cartesian_range2string(chunkGlobalRange)
            v = ba.kvStore[joinpath(mipLevelName, cartesian_range2string(chunkGlobalRange))]
            v = Codings.decode(v, C)
            chk = reinterpret(T, v)
            chk = reshape(chk, chunkSize)
            @inbounds buf[rangeInBuffer] = chk[rangeInChunk]
        catch err
            if isa(err, KeyError) && ba.fillMissing
                println("no suck key in kvstore: $(err), will fill this block as zeros")
                break
            else
                println("catch an error while getindex in BigArray: $err with type of $(typeof(err))")
                rethrow()
            end
        end
    end
end

function getindex_multithreads( ba::BigArray, idxes::Union{UnitRange, Int}...)
    t1 = time()
    sz = map(length, idxes)
    buf = zeros(eltype(ba), sz)
    baIter = ChunkIterator(idxes, get_chunk_size(ba); offset=get_offset(ba))

    @sync begin
        channel = Channel{Tuple}( CHUNK_CHANNEL_SIZE )
        @async begin
            for iter in baIter
                put!(channel, iter)
            end
            close(channel)
        end
        # control the number of concurrent requests here
        for i in 1:TASK_NUM
            @async getindex_multithreads_worker!(channel, buf, ba)
        end
    end
    elapsed = time() - t1 # seconds
    println("cutout speed: $(sizeof(buf)/1024/1024/elapsed) MB/s")
    OffsetArray(buf, idxes...)
end


text:
module AbstractPlotting

using Observables, GeometryTypes, StaticArrays, ColorTypes, Colors, IntervalSets, PlotUtils
using ColorBrewer, FixedPointNumbers, Packing, SignedDistanceFields
using Markdown # documentation
using Serialization # serialize events
using StructArrays
# Text related packages
using FreeType, FreeTypeAbstraction, UnicodeFun
using LinearAlgebra, Statistics
import ImageMagick, FileIO
import FileIO: save
using Printf: @sprintf

using Base: RefValue
using Base.Iterators: repeated, drop
import Base: getindex, setindex!, push!, append!, parent, get, get!, delete!, haskey

module ContoursHygiene
    import Contour
end
using .ContoursHygiene
const Contours = ContoursHygiene.Contour

include("utilities/quaternions.jl")
include("types.jl")
include("utilities/utilities.jl")
include("utilities/logging.jl")
include("utilities/texture_atlas.jl")
include("interaction/nodes.jl")

# Basic scene/plot/recipe interfaces + types
include("scenes.jl")
include("theming.jl")
include("recipes.jl")
include("interfaces.jl")
include("conversions.jl")

# camera types + functions
include("camera/projection_math.jl")
include("camera/camera.jl")
include("camera/camera2d.jl")
include("camera/camera3d.jl")

# some default recipes
include("basic_recipes/basic_recipes.jl")
include("basic_recipes/multiple.jl")
# layouting of plots
include("layouting/transformation.jl")
include("layouting/data_limits.jl")
include("layouting/layouting.jl")
include("layouting/boundingbox.jl")
# more default recipes
include("basic_recipes/buffers.jl")
include("basic_recipes/axis.jl")
include("basic_recipes/legend.jl")
include("basic_recipes/title.jl")

# Interactiveness
include("interaction/events.jl")
include("interaction/gui.jl")
include("interaction/interactive_api.jl")

# documentation and help functions
include("documentation/documentation.jl")
include("display.jl")

# help functions and supporting functions
export help, help_attributes, help_arguments

# Abstract/Concrete scene + plot types
export AbstractScene, SceneLike, Scene, AbstractScreen
export AbstractPlot, Combined, Atomic, Axis

# Theming, working with Plots
export Attributes, Theme, attributes, default_theme, theme, set_theme!
export title

# Node/Signal related
export Node, node, lift, map_once, to_value, on

# utilities and macros
export @recipe, @extract, @extractvalue, @key_str, @get_attribute
export broadcast_foreach, to_vector, replace_automatic!

# conversion infrastructure
export @key_str, convert_attribute, convert_arguments
export to_color, to_colormap, to_rotation, to_font, to_align, to_textsize
export to_ndim, Reverse

# Transformations
export translated, translate!, transform!, scale!, rotate!, grid, Accum, Absolute
export boundingbox, insertplots!, center!, translation
export hbox, vbox

# camera related
export AbstractCamera, EmptyCamera, Camera, Camera2D, Camera3D, cam2d!, cam2d
export campixel!, campixel, cam3d!, update_cam!, rotate_cam!, translate_cam!
export pixelarea, plots, cameracontrols, cameracontrols!, camera, events
export to_world

# picking + interactive use cases + events
export mouseover, ispressed, onpick, pick, Events, Keyboard, Mouse, mouse_selection
export register_callbacks
export window_area
export window_open
export mouse_buttons
export mouse_position
export mousedrag
export scroll
export keyboard_buttons
export unicode_input
export dropped_files
export hasfocus
export entered_window
export disconnect!, must_update, force_update!, update!, update_limits!

# currently special-cased functions (`textslider`) for example
export textslider

# gui
export slider, button, playbutton

# Raymarching algorithms
export RaymarchAlgorithm, IsoValue, Absorption, MaximumIntensityProjection, AbsorptionRGBA, IndexedAbsorptionRGBA
export Billboard

# Reexports of
# Color/Vector types convenient for 3d/2d graphics
export Quaternion, Quaternionf0, qrotation
export RGBAf0, RGBf0, VecTypes, RealVector, FRect, FRect2D, IRect2D
export FRect3D, IRect3D, Rect3D, Transformation
export IRect, FRect, Rect, Sphere, Circle
export Vec4f0, Vec3f0, Vec2f0, Point4f0, Point3f0, Point2f0
export Vec, Vec2, Vec3, Vec4, Point, Point2, Point3, Point4
export (..), GLNormalUVMesh

# conflicting identifiers
using GeometryTypes: widths
export widths, decompose

# building blocks for series recipes
export PlotList, PlotSpec

export plot!, plot

export Stepper, step!, replay_events, record_events, RecordEvents, record, VideoStream
export VideoStream, recordframe!, record
export save

# colormap stuff from PlotUtils, and showlibrary, showgradients
export clibraries, cgradients, clibrary, showlibrary, showgradients

# default icon for Makie
function icon()
    path = joinpath(dirname(pathof(AbstractPlotting)), "..", "assets", "icons")
    icons = FileIO.load.(joinpath.(path, readdir(path)))
    icons = reinterpret.(NTuple{4,UInt8}, icons)
end

const cached_logo = Ref{String}()

function logo()
    if !isassigned(cached_logo)
        cached_logo[] = download(
            "https://raw.githubusercontent.com/JuliaPlots/Makie.jl/sd/abstract/docs/src/assets/logo.png",
            joinpath(@__DIR__, "logo.png")
        )
    end
    FileIO.load(cached_logo[])
end

const has_ffmpeg = Ref(false)

const config_file = "theme.jl"
const config_path = joinpath(homedir(), ".config", "makie", config_file)

function __init__()
    pushdisplay(PlotDisplay())
    has_ffmpeg[] = try
        success(`ffmpeg -h`)
    catch
        false
    end
    cfg_path = config_path
    if isfile(cfg_path)
        theme = include(cfg_path)
        if theme isa Attributes
            set_theme!(theme)
        else
            @warn("Found config file in $(cfg_path), which doesn't return an instance of Attributes. Ignoring faulty file!")
        end
    end
end

end # module


text:
# no longer used for besseli but could be used in future for Jn, Yn
#=
@inline function down_recurrence(x, in, inp1, nu, branch)
    # this prevents us from looping through large values of nu when the loop will always return zero
    (iszero(in) || iszero(inp1)) && return zero(x)
    (isinf(inp1) || isinf(inp1)) && return in

    inm1 = in
    x2 = 2 / x
    for n in branch:-1:nu+1
        a = x2 * n
        inm1 = muladd(a, in, inp1)
        inp1 = in
        in = inm1
    end
    return inm1
end
=#
@inline function up_recurrence(x, k0, k1, nu)
    nu == 0 && return k0
    nu == 1 && return k1

    # this prevents us from looping through large values of nu when the loop will always return zero
    (iszero(k0) || iszero(k1)) && return zero(x)

    k2 = k0
    x2 = 2 / x
    for n in 1:nu-1
        a = x2 * n
        k2 = muladd(a, k1, k0)
        k0 = k1
        k1 = k2
    end
    return k2, k0
end


text:
module HerbLanguage

using Base: Symbol, vect, tail, Callable, uv_status_string, replaceuserpath, RangeStepRegular
include("structs.jl")
include("context.jl")
include("parse.jl")
include("operators.jl")
include("utils.jl")
include("structs-prob.jl")

export Type, thing, Constant, Variable, Functor, Structure, List, LPair, Predicate, Literal, Proposition, Negation, Conj, Clause, Disjunction, Recursion, Program
export Context, c_const!, c_const, c_var!, c_var, c_functor!, c_functor, c_pred!, c_pred, c_type!, c_type!, c_prop!, c_prop
export from_string, @prolog, @problog, expr_to_prolog, expr_to_problog
export is_ground, get_vars, get_arity, get_functor_name, subs
export ProbabilisticLiteral, AnnotatedDisjunction

end # module


text:
using Augmentor
using ImageCore
using Random
using ImageCore: GenericImage

"""
    make_gif(img, pl, num_sample; random_seed=1337, kwargs...)

Augment `img` with pipeline `pl` for num_sample times, and concatenate them into a 3 dimensional
image.

# Examples

The main purpose of this function is to generate a 3-dimensional image so that we could save a gif
cover using `ImageMagick.save`.

```julia
using Augmentor, ImageMagick
cover = make_gif(testpattern(RGB), FlipX(), 2)
ImageMagick.save("flipx.gif", cover; fps=1)
```

`img` can be a list of images, too. In this case, additional `kwargs` are passed to `mosaicview` so
that you could control how images are ordered.

```julia
pl = ElasticDistortion(6, scale=0.3, border=true) |>
     Rotate([10, -5, -3, 0, 3, 5, 10]) |>
     ShearX(-10:10) * ShearY(-10:10) |>
     CropSize(28, 28) |>
     Zoom(0.9:0.1:1.2)

n_samples, n_frames = 24, 10
imgs = [MNIST.convert2image(MNIST.traintensor(i)) for i in 1:n_samples]
preview = make_gif(imgs, pl, n_frames; nrow=1)
```
"""
function make_gif(img::GenericImage, pl, num_sample; post_op=center_pad ∘ drawborder, random_seed=1337)
    Random.seed!(random_seed)

    fillvalue = oneunit(eltype(img))
    frames = sym_paddedviews(
        fillvalue,
        post_op(img),
        [post_op(augment(img, pl)) for _ in 1:num_sample-1]...
    )
    cat(frames..., dims=3)
end

function make_gif(img, pl, num_sample; post_op=drawborder, random_seed=1337, kwargs...)
    fillvalue = oneunit(eltype(img[1]))

    init_frame = mosaicview(post_op.(img); kwargs...)
    frames = map(1:num_sample-1) do _
        mosaicview(map(x->post_op(augment(x, pl)), img)...; kwargs...)
    end

    frames = sym_paddedviews(fillvalue, init_frame, frames...)
    cat(frames..., dims=3)
end

"""
    center_pad(img, sz=(240, 200))

Pad img with white pixels to height:width ratio `sz[1]:sz[2]`.

Note that `sz` here is not the output size.
"""
function center_pad(img::AbstractMatrix, sz=(240, 200))
    # the default size (240, 200) is used in DemoCards
    fillvalue = oneunit(eltype(img))

    # make sure we don't shrink the image
    h, w = size(img)
    ratio = sz[1]/sz[2]
    pad_sz = h/w > ratio ? (h, round(Int, w / ratio)) : (round(Int, h * ratio), w)
    pad_sz = max.(size(img), pad_sz)

    offset = (pad_sz .- size(img)) .÷ 2
    PaddedView(fillvalue, img, ntuple(i -> -offset[i]:pad_sz[i]-offset[i]+1, ndims(img)))
end

function drawborder(img, fillvalue=colorant"pink")
    img = copy(img)
    img[1, 1:end] .= fillvalue
    img[1:end, 1] .= fillvalue
    img[end, 1:end] .= fillvalue
    img[1:end, end] .= fillvalue
    img
end


text:
module TestTreeSet

using Phylo
using DataFrames
using JuliaDB

using Test

species = ["Dog", "Cat", "Human"]
ntips = 10
df = DataFrame(species = species, count=[10, 20, 3])
observations = ["Dog", "Cat", "Dog", "Dog"]
jdb = table((species = observations, count = [1, 2, 3, 4]))

@testset "TreeSet" begin
    @test length(rand(Ultrametric(ntips), 10)) ==
        length(rand(Nonultrametric(ntips), 10))
    ts = rand(Ultrametric(ntips), 1:10)
    @test length(gettrees(ts)) == length(gettreenames(ts)) == ntrees(ts) == 10
    @test length(getleafnames(ts)) == ntips
    for name in gettreenames(ts)
        @test name ∈ 1:10
        @test length(getbranchnames(ts[name])) == ntips * 2 - 2
        @test length(getbranchnames(ts)[name]) == ntips * 2 - 2
        @test length(getleafnames(ts[name])) == ntips
    end
    for tree in gettrees(ts)
        @test length(getnodenames(tree)) == ntips * 2 - 1
        @test length(getbranchnames(tree)) == ntips * 2 - 2
        @test nleaves(tree) == ntips
    end
    @test nleaves(ts) == ntips
    @test_throws ArgumentError getleafnames(TreeSet(Dict{String, NamedTree}(),
                                                    Dict{String, Dict{String,
                                                                      Any}}()))
    @test nleaves(rand(Ultrametric(df), 20)) == length(species)
    @test nleaves(rand(Ultrametric(jdb), 10)) ==
        length(unique(observations))
    @test getleafinfo(rand(Ultrametric(df))) == df
    @test getleafinfo(rand(Ultrametric(jdb), 2)) == jdb
    @test nleaves(rand(Nonultrametric(df), 20)) == length(species)
    @test nleaves(rand(Nonultrametric(jdb), 10)) ==
        length(unique(observations))
    @test getleafinfo(rand(Nonultrametric(df))) == df
    @test getleafinfo(rand(Nonultrametric(jdb), 2)) == jdb

    @test nodetype(typeof(ts)) == nodetype(typeof(first(gettrees(ts))))
    @test branchtype(typeof(ts)) == branchtype(typeof(first(gettrees(ts))))
end
end


text:
using Revise
using Plots
using Basins
using DynamicalSystems

# C. Grebogi, S. W. McDonald, E. Ott, J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983
function grebogi_map(dz,z, p, n)
    θ = z[1]; x = z[2]
    J₀=0.3; a=1.32; b=0.9;
    dz[1]= θ + a*sin(2*θ) - b*sin(4*θ) -x*sin(θ)
    dz[1] = mod(dz[1],2π) # to avoid problems with attracto at θ=π
    dz[2]=-J₀*cos(θ)
    return
end

# dummy function to keep the initializator happy
function grebogi_map_J(J,z0, p, n)
   return
end
ds = DiscreteDynamicalSystem(grebogi_map,[1., -1.], [] , grebogi_map_J)
integ  = integrator(ds)

θg=range(0,2π,length=600)
xg=range(-0.5,0.5,length=600)

@time bsn=Basins.basins_map2D(θg, xg, integ)
#@time bns2=ChaosTools.basin_map(θg, xg, integ)

#plot(θg,xg,bsn.basin',seriestype=:heatmap)

# Estimation using box counting
@show bd = box_counting_dim(θg,xg, bsn.basin)
α = 2 - bd

# Estimation using original method
D,e,f = uncertainty_exponent(bsn, integ)
@show 2-D

D2,e,f = Basins.static_estimate(θg,xg,bsn)
@show 2-D2


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    fit(estimator, data)

Fit estimator to geospatial `data` and return a fitted estimator.
"""
function fit end

"""
    predict(estimator, var, uₒ)

Predict variable `var` with the `estimator` at point or geometry `uₒ`.
"""
function predict end

"""
    status(fitted)

Return the status of the `fitted` estimator. (e.g. the
factorization of the linear system was successful)
"""
function status end


text:
"""
    txarray,xarray,tqarray,qarray = qss1(f,x₀,tspan,Δq,tuarray,uarray)

Solve the initial value problem (IVP) for a given first-order ODE with inputs using the method of hysteretically quantized state system (QSS) or order 1.

For a single first-order explicit ODEs (aka state equation) `ẋ = f(x,u)`, with the initial value `x₀` specified at an initial time `t₀`, and possibly a piecewise (control) input `u` given by a pair of vectors `tuarray` and `uarray`of times and values, find the solution on the time span `tspan=(t₀,t₁)`. The only parameter for the method is the quantum `Δq`, while the hysteresis band has the same width as the quantum.

# Arguments
- `f`: function defining the right hand side of the differential equation.
- `x₀`: initial condition.
- `tspan`: time span.
- `Δq`: quantum.
- `tuarray`: array of times at which the the input changes to a new constant value.
- `uarray`: array of values of the piecewise constant input.

# Outputs
- `txarray`: a vector of times when the value of `x` is computed, between these values it evolves linearly.
- `xarray`: a vector of values of `x` corresponding to the times.
- `tqarray`: a vector of times when the value of piecewise constant `q` changes.
- `qarray`: a vector of values of the piecewise constant quantized `q`.

Note that for a first-order system with no external events the `x` and `q` arrays are identical up to the initial values (`x₀` before vs. after quantization).
"""
function qss1(f,x₀,tspan,Δq,tuarray,uarray)
    (t₀,t₁) = tspan             # Initial and final times for the simulation.
    if ~isempty(tuarray) && tuarray[1]==t₀  # Setting the initial value of the input:
        u = uarray[1]           # - either to the first value of the (piecewise) input,
        popfirst!(uarray)
        popfirst!(tuarray)
    else
        u = 0.0                 # - or to zero.
    end
    txarray = [t₀,]             # Initializing an array of times. The length not known apriori.
    xarray = [x₀,]              # Initializing an array of state variables.
    t = t₀                      # Initializing the current time.
    x = x₀                      # Initializing the current state.
    q = Δq*floor(x₀/Δq)         # Setting the quantized initial state.
    tqarray = [t,]              # Initializing an array of times for the quantized state variable.
    qarray = [q,]               # Initializing an array of (samples of) the quant. state variable.
    while t < t₁
        k = f(q,u,t)            # Derivative corresponding to the quantized state.
        Δt = abs(Δq/k)          # Time increment based only on f and not u.
        if isempty(tuarray) && t+Δt > t₁ # Predicted crossing of the quant. level beyond time span.
            Δt = t₁-t           # The time from now till the end.
            t = t₁              # The last time x is evaluated is at the very end of the time span.
            x += k*Δt           # The state x updated at the very end of the timespan.
            push!(txarray,t)    # Store the time when x experiences a change in derivative and
            push!(xarray,x)     # the actual value of x when that happens. No change to q and its t.
        elseif (~isempty(tuarray) && t+Δt < tuarray[1]) || (isempty(tuarray) && t+Δt <= t₁)   # If the next crossing of quant. level is before the change of u or before the end.
            t += Δt             # Set the next time.
            q += sign(k)*Δq     # Set the next quantized state.
            x = q               # State variable x is following q at the moment q changes.
            push!(txarray,t)    # Append the time to the array of times.
            push!(xarray,q)     # Append the quantized state to the array of (samples of) states.
            push!(tqarray,t)
            push!(qarray,q)     # The same for the quantized states (up to the first value).
        elseif ~isempty(tuarray) && t+Δt >= tuarray[1] # If the change in u arrives before x crosses the quant. level.
            Δt = tuarray[1]-t   # The time from now till the moment u changes.
            t = tuarray[1]      # The moment when u changes.
            u = uarray[1]       # The new value of u.
            popfirst!(tuarray)  # Discard the time so that it is not used in the next iteration.
            popfirst!(uarray)   # Discard also the value of control for the same reason.
            x += k*Δt           # The state x gets updated but q does not.
            push!(txarray,t)    # Store the time when x experiences a change in derivative and
            push!(xarray,x)     # the actual value of x when that happens. No change to q and its t.
        end
    end
    return (txarray,xarray,tqarray,qarray)
end

max_stars_repo_name: hurak/QuantizedStateSystems

max_issues_repo_name: hurak/QuantizedStateSystems

max_forks_repo_name: hurak/QuantizedStateSystems


text:
export iauPn00a
"""
Precession-nutation, IAU 2000A model:  a multi-purpose function,
supporting classical (equinox-based) use directly and CIO-based
use indirectly.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
    date1,date2  double          TT as a 2-part Julian Date (Note 1)

Returned:
    dpsi,deps    double          nutation (Note 2)
    epsa         double          mean obliquity (Note 3)
    rb           double[3][3]    frame bias matrix (Note 4)
    rp           double[3][3]    precession matrix (Note 5)
    rbp          double[3][3]    bias-precession matrix (Note 6)
    rn           double[3][3]    nutation matrix (Note 7)
    rbpn         double[3][3]    GCRS-to-true matrix (Notes 8,9)

Notes:

    1.  The TT date date1+date2 is a Julian Date, apportioned in any
    convenient way between the two arguments.  For example,
    JD(TT)=2450123.7 could be expressed in any of these ways,
    among others:

            date1          date2

        2450123.7           0.0       (JD method)
        2451545.0       -1421.3       (J2000 method)
        2400000.5       50123.2       (MJD method)
        2450123.5           0.2       (date & time method)

    The JD method is the most natural and convenient to use in
    cases where the loss of several decimal digits of resolution
    is acceptable.  The J2000 method is best matched to the way
    the argument is handled internally and will deliver the
    optimum resolution.  The MJD method and the date & time methods
    are both good compromises between resolution and convenience.

    2.  The nutation components (luni-solar + planetary, IAU 2000A) in
    longitude and obliquity are in radians and with respect to the
    equinox and ecliptic of date.  Free core nutation is omitted;
    for the utmost accuracy, use the iauPn00  function, where the
    nutation components are caller-specified.  For faster but
    slightly less accurate results, use the iauPn00b function.

    3.  The mean obliquity is consistent with the IAU 2000 precession.

    4.  The matrix rb transforms vectors from GCRS to J2000.0 mean
    equator and equinox by applying frame bias.

    5.  The matrix rp transforms vectors from J2000.0 mean equator and
    equinox to mean equator and equinox of date by applying
    precession.

    6.  The matrix rbp transforms vectors from GCRS to mean equator and
    equinox of date by applying frame bias then precession.  It is
    the product rp x rb.

    7.  The matrix rn transforms vectors from mean equator and equinox
    of date to true equator and equinox of date by applying the
    nutation (luni-solar + planetary).

    8.  The matrix rbpn transforms vectors from GCRS to true equator and
    equinox of date.  It is the product rn x rbp, applying frame
    bias, precession and nutation in that order.

    9.  The X,Y,Z coordinates of the IAU 2000A Celestial Intermediate
    Pole are elements (3,1-3) of the GCRS-to-true matrix,
    i.e. rbpn[2][0-2].

    10. It is permissible to re-use the same array in the returned
    arguments.  The arrays are filled in the order given.

Called:
    iauNut00a    nutation, IAU 2000A
    iauPn00      bias/precession/nutation results, IAU 2000

Reference:

    Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

    n.b. The celestial ephemeris origin (CEO) was renamed "celestial
        intermediate origin" (CIO) by IAU 2006 Resolution 2.

This revision:  2013 November 14

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPn00a(double date1, double date2,
#               double *dpsi, double *deps, double *epsa,
#               double rb[3][3], double rp[3][3], double rbp[3][3],
#               double rn[3][3], double rbpn[3][3])

function iauPn00a(date1::Real, date2::Real)
    # Allocate return values
    ref_dpsi = Ref{Float64}(0.0)
    ref_deps = Ref{Float64}(0.0)
    ref_epsa = Ref{Float64}(0.0)
    rb       = zeros(Float64, 3, 3)
    rp       = zeros(Float64, 3, 3)
    rbp      = zeros(Float64, 3, 3)
    rn       = zeros(Float64, 3, 3)
    rbpn     = zeros(Float64, 3, 3)

    ccall((:iauPn00a, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble},
          Ref{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble},
          Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
          convert(Float64, date1),
          convert(Float64, date2),
          ref_dpsi, ref_deps, ref_epsa, rb, rp, rbp, rn, rbpn)

    return ref_dpsi[], ref_deps[], ref_epsa[], SMatrix{3,3}(rb'), SMatrix{3,3}(rp'),
       SMatrix{3,3}(rbp'), SMatrix{3,3}(rn'), SMatrix{3,3}(rbpn')
 end


text:
module JpegTurbo

include("libjpeg.jl")

end


text:
println("Loading packages...")
using Cytof5, Random
using JLD2, FileIO

dirContains(d::String, file::String) = file in split(read(`ls $d`, String))

RESULTS_PATH = "results/flowSearch/"
DATA_PATHS = map(x -> "$RESULTS_PATH/$x", split(read(`ls $RESULTS_PATH`, String)))
DATA_PATHS = filter(d -> dirContains(d, "dat.jld2"), DATA_PATHS)

MCMC_ITER = 1000
BURN = 10000
K_MCMC = 10
L_MCMC = 5
b0PriorSd = 0.1
b1PriorScale = 0.1

function sim(path)
  @load "$(path)/dat.jld2" dat
  d = Cytof5.Model.Data(dat[:y])

  println("Generating Priors ...")
  @time c = Cytof5.Model.defaultConstants(d, K_MCMC, L_MCMC,
                                          b0PriorSd=b0PriorSd, b1PriorScale=b1PriorScale)

  println("Generating initial state ...")
  @time init = Cytof5.Model.genInitialState(c, d)

  println("Fitting Model ...");
  @time out, lastState, ll, metrics =
    Cytof5.Model.cytof5_fit(init, c, d,
                            monitors=[[:Z, :lam, :W,
                                       :b0, :b1, :v,
                                       :sig2, :mus,
                                       :alpha, :v,
                                       :eta],
                                      [:y_imputed]],
                            thins=[1, 100],
                            nmcmc=MCMC_ITER, nburn=BURN,
                            computeDIC=true, computeLPML=true,
                            printFreq=10,
                            flushOutput=true)

  println("Saving Data ...")
  @save "$(path)/mcmcout.jld2" out ll lastState c dat metrics;

  println("MCMC Completed.")
end

for path in DATA_PATHS
  println(path)
  sim(path)
end


text:
function make_gui_widget(cur_binding::Expr)
  if cur_binding.head == :escape && length(cur_binding.args) == 1
    cur_binding = cur_binding.args[1]
  end

  strip_escape!(cur_binding)

  @assert cur_binding.head == :(=)
  cur_symbol, cur_range = cur_binding.args

  Expr(
    :(=), esc(cur_symbol),
    Expr(
      :call, widget, esc(cur_range), Expr(
        :kw, :label, string(cur_symbol)
      )
    )
  )
end


text:
######### Jags program example  ###########

using StatsPlots, Jags, Statistics

ProjDir = joinpath(dirname(@__FILE__))
cd(ProjDir) do

  bones = "
  model {
     for (i in 1:nChild) {
        theta[i] ~ dnorm(0.0, 0.001);

        for (j in 1:nInd) {
           # Cumulative probability of > grade k given theta
           for (k in 1:(ncat[j]-1)) {
              logit(Q[i,j,k]) <- delta[j]*(theta[i] - gamma[j,k]);
           }
           Q[i,j,ncat[j]] <- 0;
        }

        for (j in 1:nInd) {
           # Probability of observing grade k given theta
           p[i,j,1] <- 1 - Q[i,j,1];
           for (k in 2:ncat[j]) {
              p[i,j,k] <- Q[i,j,(k-1)] - Q[i,j,k];
           }
           grade[i,j] ~ dcat(p[i,j,1:ncat[j]]);
        }
     }
  }
  "

  data = Dict{String, Any}()
  data["nChild"] = 13
  data["nInd"] = 34
  data["gamma"] = reshape([
    0.7425, 10.267, 10.5215, 9.3877, 0.2593, -0.5998,
    10.5891, 6.6701, 8.8921, 12.4275, 12.4788, 13.7778, 5.8374, 6.9485,
    13.7184, 14.3476, 4.8066, 9.1037, 10.7483, 0.3887, 3.2573, 11.6273,
    15.8842, 14.8926, 15.5487, 15.4091, 3.9216, 15.475, 0.4927, 1.3059,
    1.5012, 0.8021, 5.0022, 4.0168, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1.0153, 7.0421, 14.4242,
    17.4685, 16.7409, 16.872, 17.0061, 5.2099, 16.9406, 1.3556, 1.8793,
    1.8902, 2.3873, 6.3704, 5.1537, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, 17.4944, 2.3016, 2.497, 2.3689, 3.9525, 8.2832, 7.1053,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 3.2535, 3.2306,
    2.9495, 5.3198, 10.4988, 10.3038], 34, 4)
  data["delta"] = [
    2.9541, 0.6603, 0.7965, 1.0495, 5.7874, 3.8376, 0.6324, 0.8272,
    0.6968, 0.8747, 0.8136, 0.8246, 0.6711, 0.978, 1.1528, 1.6923,
    1.0331, 0.5381, 1.0688, 8.1123, 0.9974, 1.2656, 1.1802, 1.368,
    1.5435, 1.5006, 1.6766, 1.4297, 3.385, 3.3085, 3.4007, 2.0906,
    1.0954, 1.5329]
  data["ncat"] = [
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3,
    3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 5, 5]
  data["grade"] = reshape([
    1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 1, 2, NaN,
    2, 2, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2,
    2, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1,
    1, 1, 1, 1, 1, 1, 2, 2, 2, NaN, 2, 2, 1, 1, 1, 1, 1, 2, 1, 2,
    2, 2, 2, 2, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2, 2, 2, 2, 2, 1, 1, 1,
    1, 1, 1, NaN, NaN, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 1, NaN, NaN, 1,
    1, NaN, 2, NaN, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, NaN, 2, 2, 1, 1, 1,
    NaN, 1, 1, 1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2,
    2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, NaN, 2, 2, 1, 1, 1, 1,
    1, 1, 1, 1, 1, 1, NaN, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2,
    2, 2, 1, 1, 1, 1, 1, NaN, NaN, 2, 1, NaN, 1, 2, 2, 1, 1, 1, 1, 1,
    1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
    1, 1, 1, 1, 2, 2, 3, 2, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, NaN, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, NaN, 1, 1,
    1, 1, 1, 1, 1, 1, 1, 1, NaN, NaN, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    1, 1, 2, NaN, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1,
    1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
    2, 4, 2, 3, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 3, 5, 5, 5,
    5, 5, 5, 5, 5, 5, 5, 1, 1, 3, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 2,
    2, 3, 3, 4, 5, 5, 5, 5, 5, 5, 5, 5, 1, 1, 1, 1, 2, 3, 3, 3, 4,
    5, 5, 5, 5, 1, 1, 1, 1, 3, 3, 3, 4, 4, 5, 5, 5, 5], 13, 34)

  inits = Dict{String, Any}()
  inits["theta"] = [0.5, 1, 2, 3, 5, 6, 7, 8, 9, 12, 13, 16, 18]
  inits["grade"] = reshape([
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, 1,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, 1, NaN, NaN, 1,
    NaN, 1, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, NaN, NaN, NaN,
    NaN, NaN, 1, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, 1, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    1, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, 1, 1, NaN, NaN, 1, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, 1, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN,
    NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN], 13, 34)

  monitors = Dict{String, Bool}("theta" => true)

  jagsmodel = Jagsmodel(name="bones1",
  model=bones,
    monitor=monitors,
    #ncommands=4, nchains=1,
    #adapt=1000, nsamples=10000, thin=1,
    #deviance=true, dic=true, popt=true,
    #updatedatafile=true, updateinitfiles=true,
    #pdir=ProjDir
    );

  println("\nJagsmodel that will be used:")
  jagsmodel |> display

  @time sim = jags(jagsmodel, data, [inits], ProjDir)
  println()

  ## Plotting
  p = plot(sim)

end #cd


text:
# Copyright 2017-2018 Gandalf Software, Inc. (Scott Paul Jones)
# Licensed under MIT License, see LICENSE.md

ValidatedStyle(::Type{<:AbstractString}) = UnknownValidity()
ValidatedStyle(::Type{<:Str}) = AlwaysValid()

ValidatedStyle(A::T) where {T<:Union{AbsChar,AbstractString}} = ValidatedStyle(T)

CharSetStyle(::Type{<:Str{C}}) where {C<:CSE} = CharSetStyle(C)

# Different character sets
is_valid(::Type{Str{<:CSE{CS1}}}, str::T) where {T<:Str{<:CSE{CS2}},CS1} where {CS2} =
    _isvalid(ValidatedStyle(T), CS1, CS2, str)

# Same character set
is_valid(::Type{Str{<:CSE{CS}}}, str::T) where {T<:Str{<:CSE{CS}}} where {CS} =
    _isvalid(ValidatedStyle(T), str)

_isvalid(::UnknownValidity, str::T) where {T<:Str} = _isvalid(T, pointer(str), ncodeunits(str))

is_valid(::Type{T}, str::T) where {T<:Str}       = _isvalid(ValidatedStyle(T), str)
is_valid(str::T) where {T<:Str}       = _isvalid(ValidatedStyle(T), str)

# For now, there is only support for immutable `Str`s, when mutable `Str`s are added.

"""
    MutableStyle(A)
    MutableStyle(typeof(A))

`MutableStyle` specifies the whether a string type is mutable or not

    MutableStyle(::Type{<:MyString}) = MutableStr()

The default is `ImmutableStr`
"""
abstract type MutableStyle end
struct ImmutableStr <: MutableStyle end
struct MutableStr   <: MutableStyle end

MutableStyle(A::AbstractString) = MutableStyle(typeof(A))

MutableStyle(A::Str) = ImmutableStr()

_ismutable(::ImmutableStr, str::Type{T}) where {T<:Str} = false
_ismutable(::MutableStr, str::Type{T}) where {T<:Str} = true
is_mutable(::Type{T}) where {T<:Str} = _ismutable(MutableStyle(T))

isimmutable(str::T) where {T<:Str} = !is_mutable(T)

"""
    CanContain(Union{A, typeof(A)}, Union{B, typeof(B)})

`CanContainStyle` specifies whether the first string can contain a substring of the second type,
    and if so, what is the most efficient method of comparison

    CanContain(::Type{<:MyString}, ::Type{String}) = ByteCompare()

Returns an instance of type `CompareStyle`, default `CodePointCompare`
"""
CanContain(::Type{<:CSE}, ::Type{<:CSE}) = CodePointCompare()

CanContain(::Type{C}, ::Type{C}) where {C<:CSE} = ByteCompare()

CanContain(::Type{<:Binary_CSEs}, ::Type{<:Union{_UCS2CSE,_UTF32CSE}}) =
    NoCompare()
CanContain(::Type{<:Binary_CSEs}, ::Type{<:Byte_CSEs}) =
    ByteCompare()
CanContain(::Type{<:Binary_CSEs}, ::Type{<:Union{Word_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{ASCIICSE}, ::Type{<:SubSet_CSEs}) =
    NoCompare()
CanContain(::Type{ASCIICSE}, ::Type{<:Union{Binary_CSEs, LatinCSE, UTF8_CSEs}}) =
    ByteCompare()
CanContain(::Type{ASCIICSE}, ::Type{<:Union{Word_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{<:Latin_CSEs}, ::Type{<:Union{_UCS2CSE,_UTF32CSE}}) =
    NoCompare()
CanContain(::Type{<:Latin_CSEs}, ::Type{<:Union{Binary_CSEs,ASCIICSE,Latin_CSEs}}) =
    ByteCompare()
CanContain(::Type{<:Latin_CSEs}, ::Type{<:UTF8_CSEs}) =
    ASCIICompare()
CanContain(::Type{<:Latin_CSEs}, ::Type{<:Union{Word_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{<:UTF8_CSEs}, ::Type{<:Union{ASCIICSE,Binary_CSEs}}) =
    ByteCompare()
CanContain(::Type{<:UTF8_CSEs}, ::Type{<:Latin_CSEs}) =
    ASCIICompare()

CanContain(::Type{<:Union{Text2CSE,UCS2CSE}}, ::Type{_UTF32CSE}) =
    NoCompare()
CanContain(::Type{<:Union{Text2CSE,UCS2_CSEs}}, ::Type{<:Word_CSEs}) =
    ByteCompare()
CanContain(::Type{<:Union{Text2CSE,UCS2_CSEs}},
           ::Type{<:Union{ASCIICSE, Binary_CSEs, Latin_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{UTF16CSE}, ::Type{<:Union{ASCIICSE,Binary_CSEs,Latin_CSEs}}) =
    WidenCompare()
CanContain(::Type{UTF16CSE}, ::Type{<:Union{Text2CSE, UCS2_CSEs}}) =
    ByteCompare()

CanContain(::Type{<:Quad_CSEs}, ::Type{<:Quad_CSEs}) =
    ByteCompare()
CanContain(::Type{<:Quad_CSEs},
           ::Type{<:Union{Binary_CSEs,ASCIICSE,Latin_CSEs,Text2CSE,UCS2_CSEs}}) =
    WidenCompare()

CanContain(::Type{S}, ::Type{T}) where {S<:AbstractString, T<:AbstractString} =
    CanContain(cse(S), cse(T))
CanContain(::Type{T}, ::Type{T}) where {T<:AbstractString} = ByteCompare()

CanContain(A::AbstractString, B::AbstractString) = CanContain(typeof(A), typeof(B))

"""
    EqualsStyle(Union{A, typeof(A)}, Union{B, typeof(B)})

`EqualsStyle` specifies whether the first string can equal a substring of the second type,
    and if so, what is the most efficient method of comparison

This is determined by the CanContain trait
Returns an instance of type `CompareStyle`, default `CodePointCompare`
"""
EqualsStyle(::Type{S}, ::Type{T}) where {S<:CSE,T<:CSE} =
    Base.promote_typeof(CanContain(S, T), CanContain(T, S))()
EqualsStyle(::Type{T}, ::Type{T}) where {T<:CSE} = ByteCompare()

EqualsStyle(::Type{S}, ::Type{T}) where {S<:AbstractString, T<:AbstractString} =
    EqualsStyle(cse(S), cse(T))
EqualsStyle(::Type{T}, ::Type{T}) where {T<:AbstractString} = ByteCompare()

EqualsStyle(A::AbstractString, B::AbstractString) =  EqualsStyle(typeof(A), typeof(B))


text:
module Kucoin

using HTTP: Response, URI, request, iserror
using HTTP.WebSockets: WebSocket, WebSockets  # open
using SHA: hmac_sha256
using JSON3: JSON3  # Array, Object, read, write
using Base64: base64encode
using UUIDs: uuid1
import Base: print, take!, fetch, isready

# REST API
export ApiData
export get_symbols_list, get_all_tickers
export list_accounts, get_trade_fees
export place_limit_order, place_market_order, cancel_all_orders
export get_ws_token

# WebSocket API
export WebSocketClient
export subscribe, unsubscribe, take!, fetch, isready

"""
    Kucoin.ApiData(key, secret, passphrase)

Holds API data susch as API key, secret, passphrase, etc.

# Fields
- `key::String`: kucoin api key
- `secret::String`: kucoin api secret
- `passphrase::String`: kucoin api passphrase
- `encrypted_passphrase::String`: encrypted passpharase
"""
Base.@kwdef struct ApiData
    key::String
    secret::String
    passphrase::String
    encrypted_passphrase::String = _sign(secret, passphrase)
end

include("utils.jl")
include("market.jl")
include("user.jl")
include("trade.jl")
include("websockets.jl")

end


text:
"""
Common interface for all modes and key lengths
"""
function encrypt end

encrypt(plain, cipher::Cipher{CBC,T,U}; iv=rand(UInt8, 16)) where {T,U} = encrypt(transcode(UInt8, plain), cipher, iv)
encrypt(plain, cipher::Cipher{CTR,T,U}; iv=rand(UInt8, 16)) where {T,U} = encrypt(transcode(UInt8, plain), cipher, iv)
encrypt(plain, cipher::Cipher) = encrypt(transcode(UInt8, plain), cipher)

encrypt(plaintext::AbstractArray{UInt8}, cipher::Cipher{CBC,T,U}, iv) where {T,U} = CipherText(AESCBC(plaintext, iv, cipher.key, cipher.cache), iv, get_key_length(cipher), CBC)
encrypt(plaintext::AbstractArray{UInt8}, cipher::Cipher{CTR,T,U}, iv) where {T,U} = CipherText(AESCTR(plaintext, iv, cipher.key, cipher.cache), iv, get_key_length(cipher), CTR)
encrypt(plaintext::AbstractArray{UInt8}, cipher::Cipher{ECB,T,U}) where {T,U} = CipherText(AESECB(plaintext, cipher.key, cipher.cache), nothing, get_key_length(cipher), ECB)

encrypt!(ciphertext, plaintext, cipher::Cipher{CBC,T,U}; iv=rand(UInt8, 16)) where {T,U} = encrypt!(ciphertext, transcode(UInt8, plaintext), cipher, iv)
encrypt!(ciphertext, plaintext, cipher::Cipher{CTR,T,U}; iv=rand(UInt8, 16)) where {T,U} = encrypt!(ciphertext, transcode(UInt8, plaintext), cipher, iv)
encrypt!(ciphertext, plaintext, cipher::Cipher) = encrypt!(ciphertext, transcode(UInt8, plaintext), cipher)

encrypt!(ciphertext, plaintext::AbstractArray{UInt8}, cipher::Cipher{CBC,T,U}, iv) where {T,U} = CipherText(AESCBC!(ciphertext, plaintext, iv, cipher.key, cipher.cache), iv, get_key_length(cipher), CBC)
encrypt!(ciphertext, plaintext::AbstractArray{UInt8}, cipher::Cipher{CTR,T,U}, iv) where {T,U} = CipherText(AESCTR!(ciphertext, plaintext, iv, cipher.key, cipher.cache), iv, get_key_length(cipher), CTR)
encrypt!(ciphertext, plaintext::AbstractArray{UInt8}, cipher::Cipher{ECB,T,U}) where {T,U} = CipherText(AESECB!(ciphertext, plaintext, cipher.key, cipher.cache), nothing, get_key_length(cipher), ECB)


text:
import Statistics
import StatsBase

"""
    r2_score(ytrue, ypred)

Computes coefficient of determination. Higher values are better. Best value
is 1.
"""
function r2_score(
        ytrue::AbstractVector{<:Real},
        ypred::AbstractVector{<:Real},
        )
    if length(ytrue) != length(ypred)
        error("length(ytrue) != length(ypred)")
    end
    if length(ytrue) == 0
        error("length(ytrue) == 0")
    end
    # ybar = mean of the true y values
    ybar = Statistics.mean(ytrue)
    # SStot = total sum of squares
    SStot = sum( (ytrue .- ybar).^2 )
    # SSres = sum of squares of residuals
    residuals = ytrue .- ypred
    SSres = sum( residuals.^2 )
    R2 = 1 - SSres/SStot
    return R2
end


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

# TODO: remove the type stubs below between 0.7 and 1.0

"""
    RowVector(vector)

A lazy-view wrapper of an [`AbstractVector`](@ref), which turns a length-`n` vector into a `1×n`
shaped row vector and represents the transpose of a vector (the elements are also transposed
recursively).

By convention, a vector can be multiplied by a matrix on its left (`A * v`) whereas a row
vector can be multiplied by a matrix on its right (such that `RowVector(v) * A = RowVector(Transpose(A) * v)`). It
differs from a `1×n`-sized matrix by the facts that its transpose returns a vector and the
inner product `RowVector(v1) * v2` returns a scalar, but will otherwise behave similarly.

# Examples
```jldoctest
julia> a = [1; 2; 3; 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> RowVector(a)
1×4 RowVector{Int64,Array{Int64,1}}:
 1  2  3  4

julia> RowVector(a)[3]
3

julia> RowVector(a)[1,3]
3

julia> RowVector(a)[3,1]
ERROR: BoundsError: attempt to access 1×4 RowVector{Int64,Array{Int64,1}} at index [3, 1]
[...]

julia> RowVector(a)*a
30

julia> B = [1 2; 3 4; 5 6; 7 8]
4×2 Array{Int64,2}:
 1  2
 3  4
 5  6
 7  8

julia> RowVector(a)*B
1×2 RowVector{Int64,Array{Int64,1}}:
 50  60
```
"""
struct RowVector{T,V<:AbstractVector} <: AbstractMatrix{T}
    vec::V
    function RowVector{T,V}(v::V) where V<:AbstractVector where T
        check_types(T,v)
        new(v)
    end
end
const ConjRowVector{T,CV<:ConjVector} = RowVector{T,CV}

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
"""
    tofs(s, T1, T2)

Time change mapping `t` in `[T1, T2]` (``X``-time) to `s` in `[T1, T2]` (`U`-time).
"""
tofs(s, T1, T2) = T1 .+ (s .- T1).*(2.0 .- (s .- T1)/(T2-T1))

"""
    soft(t, T1, T2)

Time change mapping `s` in `[T1, T2]`` (``U``-time) to `t`` in `[T1, T2]` (`X`-time).
"""
soft(t, T1, T2, T=T2-T1) = T2 .- sqrt.(T*(T2 .- t))

xofu(s, u, T1, T2, v, P) = Vs(s, T1, T2, v, P) .- (T2-s)*u
uofx(s, x, T1, T2, v, P) = (Vs(s, T1, T2, v, P) .- x)/(T2-s)
txofsu(s, u, T1, T2, v, P) = (tofs(s, T1, T2), xofu(s, u, T1, T2, v, P))

"""
    Vs(s, T1, T2, v, P)

Time changed `V` for generation of `U`.
"""
function Vs(s, T1, T2, v, P::LinPro, phim = exp(-P.B*(T2-T1)*(1.0 - (s-T1)/(T2-T1))^2))
    phim*( v .- P.μ) .-  P.μ
end
Vs(s, T1, T2, v, P::Ptilde) = V(tofs(s, T1, T2), T2, v, P)

"""
    dotVs (s, T1, T2, v, P)

Time changed time derivative of `V` for generation of `U`.
"""
function dotVs(s, T1, T2, v, P::LinPro, phim = exp(-P.B*(T2-T1)*(1.0 - (s-T1)/(T2-T1))^2))
    phim*( P.B*v .+ P.beta)
end
dotVs(s, T1, T2, v, P::Ptilde) = dotV(tofs(s, T1, T2), T2, v, P)

function Ju(s, T1, T2, P::LinPro, x, phim = exp(-P.B*(T2-T1)*(1.0 - (s-T1)/(T2-T1))^2))
     sl = P.lambda*(T2-T1)/(T2-s)^2
    ( phim*sl*phim'-sl)\x
end

function Ju(s, T1, T2, P::Ptilde, x)
    H(tofs(s, T1, T2), T2, P, x)*(T2-s)^2/(T2-T1)
end
function J(s, T1, T2, P::Ptilde)
    H(tofs(s, T1, T2), T2, P)*(T2-s)^2/(T2-T1)
end

ubridge(W, Po) = ubridge!(copy(W), W, Po)
function ubridge!(X, W::SamplePath{T}, Po) where T
    T1 = Po.tt[1]
    T2 = Po.tt[end]
    v = Po.v[2]
    Pt = ptilde(Po)
    N = length(W)
    N != length(X) && error("X and W differ in length.")

    ss = W.tt
    ww = W.yy
    tt = X.tt
    xx = X.yy

    u = uofx(T1, Po.v[1], T1, T2, v, Pt)

    for i in 1:N-1
        s = ss[i]
        t, x = txofsu(s, u, T1, T2, v, Pt)
        tt[i], xx[.., i] = t, x
        bU = 2/(T2-T1)*dotVs(s, T1, T2, v, Pt) - 2/(T2-T1)*b(t, x, target(Po)) +   1/(T2-s)*(u - 2.0*a(t, x, Po)*Ju(s, T1, T2, Pt, u) )
        σUdW = (-sqrt(2.0/((T2-T1)*(T2-s))))*(σ(t, x, Po)*(ww[.., i+1]-ww[..,i]))
        u += bU*(ss[i+1]-s) + σUdW
    end
    xx[.., N] = v
    X
end

uthetamethod(W, P, theta=0.5) = uthetamethod!(copy(W), W, Po, theta)
function uthetamethod!(Y, u, W::SamplePath, Po, theta=0.5)

    T1 = Po.t0
    T2 = Po.t1
    v = Po.v1
    Pt = ptilde(Po)
    N = length(W)
    N != length(X) && error("X and W differ in length.")

    @assert(constdiff(P))

    ss = W.tt
    ww = W.yy
    tt = X.tt
    xx = X.yy

    u = uofx(T1, Po.v0, T1, T2, v, Pt)

    for i in 1:N-2 # fix me
        s = ss[i]
        t, x = txofsu(s, u, T1, T2, v, Pt)
        tt[i], xx[.., i] = t, x
        x2 = x
        dw = (ww[.., i+1]-ww[..,i])
        delta1 = 2/(T2-T1)*dotVs(s, T1, T2, v, Pt) - 2/(T2-T1)*b(t, x, target(Po)) +   1/(T2-s)*(u - 2.0*a(t, x, Po)*Ju(s, T1, T2, Pt, u) )*(ss[i+1]-s)
        local delta2

        eps2 = 5e-6
        MM = 8
        for mm in 1:MM

            delta2 = 2/(T2-T1)*dotVs(s, T1, T2, v, Pt) - 2/(T2-T1)*b(t, x2, target(Po)) +   1/(T2-s)*(u2 - 2.0*a(t, x2, Po)*Ju(s, T1, T2, Pt, u2) )*(ss[i+1]-s)
            buderiv = -2/(T2-T1)*bderiv(t, x2, target(Po)) + 1/(T2-s)*(I - 2.0*a(t, x, Po)*J(s, T1, T2, Pt))*(ss[i+1]-s)
            bderiv(tt[i+1], y2, P)
            dy2 = -inv(I - theta*(buderiv*(s[i+1]-s)))*(u2 - u - (1-theta)*delta1 - theta*delta2 - σ(tt[i], y, P)*dw)

            x2 += dx2
            x2 = txofsu(s, u2, T1, T2, v, Pt)
            if  maximum(abs(dx2)) < eps2
                break;
            end

            if mm == MM
                warn("thetamethod: no convergence $i $y $y2  $dy2")
            end
        end
        u = u + (1-theta)*delta1 + theta*delta2 + σ(tt[i], x, P)*dw
    end
    yy[.., N-1] = y
    Y
end

# using left approximation
function ullikelihood(Y::SamplePath{T}, Po) where T
    yy = Y.yy
    tt = Y.tt
    T1 = Po.tt[1]
    T2 = Po.tt[end]
    v = Po.v[2]
    P = target(Po)
    Pt = ptilde(Po)
    s2 = soft(tt[1], T1, T2)
    som::Float64 = 0.
    for i in 1:length(tt)-1 #skip last value, summing over n-1 elements
        t = tt[i]
        x = yy[.., i]
        s = s2
        s2 = soft(tt[i+1], T1, T2)

        j = J(s, T1, T2, Pt)
        ju = j*uofx(s, yy[.., i], T1, T2, v, Pt)
        som += 2.0*dot(b(t, x, P)  - b(t,x, Pt),ju)*(s2-s)

        if !constdiff(Po)
            ad = a(t,x, P) - a(t,x, Pt)
            som += -1.0/(T2-s)*(tr(j*ad) - T*dot(ju,ad*ju))*(s2-s)
        end
    end
    som
end
function ullikelihoodtrapez(Y::SamplePath{T}, Po) where T
    yy = Y.yy
    tt = Y.tt
    T1 = Po.tt[1]
    T2 = Po.tt[end]
    v = Po.v[2]
    P = target(Po)
    Pt = ptilde(Po)
    ss = soft(tt, T1, T2)
    som::Float64 = 0.
    for i in [1]
        t = tt[i]
        x = yy[.., i]
        j = J(ss[i], T1, T2, Pt)
        ju = j*uofx(ss[i], yy[.., i], T1, T2, v, Pt)
        som += dot(b(t, x, P)  - b(t,x, Pt),ju)*(ss[2]-ss[1])
    end

    for i in 2:length(tt)-1 #skip last value, summing over n-1 elements
        t = tt[i]
        x = yy[.., i]
        j = J(ss[i], T1, T2, Pt)
        ju = j*uofx(ss[i], yy[.., i], T1, T2, v, Pt)
        som += dot(b(t, x, P)  - b(t,x, Pt),ju)*(ss[i+1]-ss[i-1])

        if !constdiff(Po)
            error("not implemented")
        end
    end
    som
end

uinnovations(Y, Po) = uinnovations!(copy(Y), Y, Po)
function uinnovations!(W, Y::SamplePath{T}, Po) where T

    N = length(W)
    N != length(Y) && error("Y and W differ in length.")

    yy = Y.yy
    tt = Y.tt
    ww = W.yy
    ss = W.tt
    Pt = ptilde(Po)

    T1 = Po.tt[1]
    T2 = Po.tt[end]
    v = Po.v[2]

    w = zero(ww[.., 1])
    s = s2 = soft(tt[1], T1, T2)
    u2 = uofx(s2, yy[.., 1], T1, T2, v, Pt)

    for i in 1:N-1
        t = tt[i]
        s, u = s2, u2
        ww[.., i] = w
        ss[i] = s

        s2 = soft(tt[i+1], T1, T2)
        u2 = uofx(s2, yy[.., i+1], T1, T2, v, Pt)

        bU = 2/(T2-T1)*dotVs(s, T1, T2, v, Pt) - 2/(T2-T1)*b(t,  yy[.., i], target(Po)) +   1/(T2-s)*(u - 2.0*a(t,  yy[.., i], Po)*Ju(s, T1, T2, Pt, u) )
        σU = -sqrt(2.0/((T2-T1)*(T2-s)))*σ(t, yy[.., i], Po)

        w = w + inv(σU)*(u2 - u - bU*(s2 - s))
    end
    ww[.., N] = ww[.., N-1] + randn(typeof(ww[.., N]))*sqrt(s2 - s)
    SamplePath{T}(ss, ww)
end


text:
function uview_col(x::SparseMatrixCSC, j::Integer)
    vals = nonzeros(x)
    rng = nzrange(x, j)
    return view(vals, rng)
end

"""
for each period `t`, we need to compute
    1. State transitions: Pr(sₜ₊₁ | sₜ)
    2. Expected value of Type I extreme value shocks
"""
# function update_sparse_state_transition!(simtmp::SimulationTmp,  uv::NTuple{2,<:Real}, z::NTuple{NZ,Real}, itypidx::NTuple{NI,Real}, state_if_no_drilling::Integer) where {NZ,NI}

function update_sparse_state_transition!(simtmp::SimulationTmp,
    obs::ObservationDrill, sim::SimulationDraw, theta)

    model = _model(obs)
    wp = statespace(model)
    P = Pprime(simtmp)

    vals = nonzeros(P)
    fill!(vals, 0)

    states_after_drilling = end_lrn(wp)+1 : end_inf(wp)
    states_can_start_from = flatten( (_x(obs), states_after_drilling, )  )

    for si in states_can_start_from

        obs_cf = ObservationDrill(model, ichars(obs), zchars(obs), 0, si)
        actions = actionspace(obs_cf)
        dp1sp = actions .+ 1

        local ubV = view(vdfull(simtmp), dp1sp)
        local actionprobs = uview_col(P, si)

        # if we can't drill
        if end_ex0(wp) < si <= end_lrn(wp) || si == end_inf(wp)
            @assert length(actions) == 1
            fill!(ubV, 0)
            fill!(actionprobs, 1)
        else
            for d in actions
                ubV[d+1] = full_payoff!(theta, d, obs_cf, theta, sim, false)
            end
        end

        # this is ok to do b/c logsumexp(x::Float64) = x
        lse = logsumexp!(actionprobs, ubV)
        isfinite(lse) || throw(error("lse = $lse with actionprobs = $actionprobs and ubV = $ubV"))
        eeps = lse - dot(actionprobs, ubV)
        isfinite(eeps) || throw(error("Eeps = $eeps with actionprobs = $actionprobs and ubV = $ubV and lse=$lse"))
        setindex!(Eeps(simtmp), eeps, si)
      end
end

struct LeaseCounterfactual{L<:DrillLease} <: AbstractObservationGroup
    data::L
end

zchars( x::LeaseCounterfactual) = zchars(DataDrill(x))

@delegate LeaseCounterfactual.data [jtstart, tstart, uniti, _regime, ichars]

firstindex(x::LeaseCounterfactual) = tstart(_data(x))   # tptr[j]
length(    x::LeaseCounterfactual) = length(zchars(x)) - jtstart(_data(x)) + 1
lastindex( x::LeaseCounterfactual) =  tstart(_data(x)) + length(x)-1
eachindex( x::LeaseCounterfactual) = firstindex(x) : lastindex(x)

first_state(x::LeaseCounterfactual) = _x(DataDrill(x), tstart(x))

function state_if_never_drilled(x::LeaseCounterfactual, t)
    model = _model(DataDrill(x))
    wp = statespace(model)
    return state_if_never_drilled(wp, first_state(x), t)
end

function getindex(x::LeaseCounterfactual, t)
    time_since_start = t - tstart(x)
    zt = jtstart(x) + time_since_start
    d = DataDrill(x)
    x0 = state_if_never_drilled(x,time_since_start)
    return ObservationDrill(_model(d), ichars(x), zchars(d,zt), 0, x0), zt
end


text:
#
# Random initial set of parameters
#
function initP!(p0, options :: Options, lower, upper)
  for i in 1:length(p0)
    pmin = max(options.p0_range[1],lower[i])
    pmax = min(options.p0_range[2],upper[i])
    p0[i] = pmin + rand()*(pmax-pmin)
  end
  return p0
end


text:
# This is a simple script which implements the "graphics hello world" described in chapter 2

# constants
rows = 200;
cols = 100;
max_color = 255; # colors go from 0-255
fio = open("outputs\\helloWorld.ppm", "w");
# format meta information, P3 format name, second line tells rows by column
println(fio, "P3\n$rows $cols\n$max_color");
# actual pixels
rgb = [
    [i / rows, j / cols, 0.2]
    for
    i in range(0, rows - 1, step = 1), j in range(cols - 1, 0, step = -1)
]
[
    println(fio, "$(c[1]) $(c[2]) $(c[3])")
    for c in map(x -> round.(Int, x * max_color), rgb)
]
close(fio);


text:
#
# git clone https://github.com/m3g/M3GTools
#

# Load M3GTools
using M3GTools

# Load simulation data
mysim = M3GTools.init(psf="../simulation_files/structure.psf",
                  pdb="../simulation_files/structure.pdb",
                  parfiles = [ "/home/leandro/programs/toppar/charmm/par_all36_prot.prm",
                               "/home/leandro/programs/toppar/charmm/toppar_water_ions.str",
                               "/home/leandro/programs/toppar/charmm/par_all36_lipid.prm" ],
                  vmd="vmd")

# Select atoms of residue number 33 (be careful with the fact that VMD starts counting from zero):
sel1 = M3GTools.select(mysim,"protein and residue 33")

# Select atoms of residue number 29
sel2 = M3GTools.select(mysim,"protein and residue 29")

# Compute electrostatic interaction between these two selections
elecenergy = M3GTools.coulomb(mysim.atoms,sel1,sel2);

# Compute vdw interaction between these two selections
vdwenergy = M3GTools.vdw(mysim.atoms,sel1,sel2);

println("Electrostatic interaction between residues 30 and 34 = $elecenergy")
println("vdW interaction between residues 30 and 34 = $vdwenergy")

# Compute the number of hydrogen bonds between two selections
nhb = M3GTools.nhbonds(mysim.atoms,mysim.bonds,sel1,sel2)
println("Number of H-bonds between residues 30 and 34 = $nhb")


text:
# Copyright (c) 2021-2022, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon the Sudoku module being defined first.
using Test

function valid_subarray()
    @test Sudoku.valid_subarray(3,[1,2,3,4,5,6,7,8,9])
    @test Sudoku.valid_subarray(3,[1,2,3,7,8,9,4,5,6])
    @test !Sudoku.valid_subarray(3,[2,2,3,4,5,6,7,8,9])
end

function valid_puzzle_2()
    p = Sudoku.puzzle_n(2)
    @test Sudoku.valid_puzzle(p)
    p[1,1] = p[1,2] # invalidate puzzle
    @test !Sudoku.valid_puzzle(p)
end

function satisfies_puzzle_2()
    p = Sudoku.puzzle_n(2)
    s = similar(p)
    s .= p
    p[1,1] = 0
    @test Sudoku.satisfies(p,s)
    Sudoku.row_permute!(p,2,1,2)
    @test !Sudoku.satisfies(p,s)
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
function add_tracing_callbacks(b::Gtk.GtkBuilder,handles::Tracker_Handles)

    signal_connect(combine_cb,b["tracing_combine_button"],"clicked",Void,(),false,(handles,))

    nothing
end

function combine_cb(w::Ptr,user_data::Tuple{Tracker_Handles})

    han, = user_data

    han.combine_mode = getproperty(han.b["tracing_combine_button"],:active,Bool)

    nothing
end

function combine_start(han,x,y)

    println("start")
    han.partial = deepcopy(han.wt.whiskers[han.woi_id])
    han.combine_mode = 2

    nothing
end

ccw(x1,x2,x3,y1,y2,y3)=(y3-y1) * (x2-x1) > (y2-y1) * (x3-x1)

function intersect(x1,x2,x3,x4,y1,y2,y3,y4)
    (ccw(x1,x3,x4,y1,y3,y4) != ccw(x2,x3,x4,y2,y3,y4))&&(ccw(x1,x2,x3,y1,y2,y3) != ccw(x1,x2,x4,y1,y2,y4))
end

function combine_end(han,x,y)

    out1=1
    out2=1

    for i=2:han.partial.len
        for j=2:han.wt.whiskers[han.woi_id].len
            if intersect(han.partial.x[i-1],han.partial.x[i],han.wt.whiskers[han.woi_id].x[j-1],
                han.wt.whiskers[han.woi_id].x[j],han.partial.y[i-1],han.partial.y[i],
                han.wt.whiskers[han.woi_id].y[j-1],han.wt.whiskers[han.woi_id].y[j])
                out1=i
                out2=j
                break
            end
        end
    end

    if out1==1
        println("No intersection found, looking for closest match")
        d=1000.0
        for i=2:han.partial.len
            for j=2:han.wt.whiskers[han.woi_id].len
                this_d = sqrt((han.partial.x[i]-han.wt.whiskers[han.woi_id].x[j]).^2+(han.partial.y[i]-han.wt.whiskers[han.woi_id].y[j]).^2)
                if this_d<d
                    d = this_d
                    out1=i
                    out2=j
                end
            end
        end
    end

    println("Segments combined")
    new_x = [han.wt.whiskers[han.woi_id].x[1:out2]; han.partial.x[out1:end]]
    new_y = [han.wt.whiskers[han.woi_id].y[1:out2]; han.partial.y[out1:end]]
    new_scores = [han.wt.whiskers[han.woi_id].scores[1:out2]; han.partial.scores[out1:end]]
    new_thick = [han.wt.whiskers[han.woi_id].thick[1:out2]; han.partial.thick[out1:end]]
    han.woi[han.displayed_frame].x=new_x
    han.woi[han.displayed_frame].y=new_y
    han.woi[han.displayed_frame].thick=new_thick
    han.woi[han.displayed_frame].scores=new_scores
    han.woi[han.displayed_frame].len = length(new_thick)

    han.combine_mode = 1
    redraw_all(han)

    nothing
end


text:
# DiffEqSolution to Array
function diffeqsol_to_array(x)
    xarr = cpu(x)
    return reshape(xarr, size(xarr)[1:2])
end

function diffeqsol_to_cuarray(x)
    xarr = gpu(x)
    return reshape(xarr, size(xarr)[1:2])
end

function diffeqsol_to_trackedarray(x)
    xarr = x.u[end]
    return reshape(xarr, size(xarr)[1:2])
end

function diffeqsol_to_3dtrackedarray(x)
    return hcat(map(_x -> reshape(_x, size(_x, 1), 1, size(_x, 2)), x.u)...)
end

_convert_tspan(tspan, p) = eltype(p).(tspan)

_convert_tspan(tspan, p::TrackedArray) = Tracker.collect(eltype(p).(tspan))

norm_batched(x::AbstractArray) = sqrt.(sum(x .^ 2, dims = 1))

# Running Average Meter
mutable struct AverageMeter{T}
    last_value::T
    sum::T
    count::Int

    AverageMeter(T = Float32) = new{T}(T(0), T(0), 0)
end

function reset!(am::AverageMeter{T}) where {T}
    am.last_value = T(0)
    am.sum = T(0)
    am.count = 0
    return am
end

function update!(am::AverageMeter{T}, val::T) where {T}
    am.last_value = val
    am.sum += val
    am.count += 1
    return am.sum / am.count
end

(am::AverageMeter)() = am.sum / am.count

# Pretty Logger
function table_logger(header::Vector{String}, record::Vector{String} = [])
    n = length(header) + length(record)
    ind_lens = vcat(length.(header), length.(record))
    span = sum(ind_lens .+ 3) + 1
    println("="^span)
    for h in vcat(header, record)
        print("| $h ")
    end
    println("|")
    println("="^span)

    avg_meters = Dict{String,AverageMeter}(rec => AverageMeter() for rec in record)

    patterns = ["%$l.4f" for l in ind_lens]
    fmtrfuncs = generate_formatter.(patterns)
    function internal_logger(last::Bool, records::Dict, args::Vararg)
        if length(records) > 0
            for (rec, val) in records
                update!(avg_meters[rec], val)
            end
            return
        end
        if last
            println("="^span)
            return
        end
        for h in [
            fmtrfunc(arg) for (fmtrfunc, arg) in
            zip(fmtrfuncs, vcat([args...], [avg_meters[rec]() for rec in record]))
        ]
            print("| $h ")
        end
        println("|")
    end
end

# MvNormal for FFJORD
struct BatchedMultiVariateNormal{M,N,O,P,Q}
    μ::M
    cov::N
    inv_cov::N
    det_cov::O
    lt_decom::P
    k::Q
end

Flux.@functor BatchedMultiVariateNormal

function BatchedMultiVariateNormal(μ::AbstractVector, cov::AbstractMatrix)
    inv_cov = inv(cov)
    det_cov = det(cov)
    lt_decom = Array(cholesky(cov))
    return BatchedMultiVariateNormal(
        reshape(μ, :, 1),
        cov,
        inv_cov,
        det_cov,
        lt_decom,
        length(μ),
    )
end

function (mvnorm::BatchedMultiVariateNormal)(x::AbstractMatrix{T}) where {T}
    denom = sqrt((T(2π)^mvnorm.k) * T(mvnorm.det_cov))
    diff = x .- mvnorm.μ
    return exp.(diag(-permutedims(diff, (2, 1)) * (mvnorm.inv_cov * diff) / 2)) / denom
end

function sample(mvnorm::BatchedMultiVariateNormal{CuArray{T,2}}, nsamples::Int) where {T}
    samples = CUDA.randn(T, mvnorm.k, nsamples)::CuArray{T,2}
    return mvnorm.μ .+ mvnorm.lt_decom * samples
end

function sample(
    mvnorm::BatchedMultiVariateNormal{TrackedArray{T,2,CuArray{T,2}}},
    nsamples::Int,
) where {T}
    samples = CUDA.randn(T, mvnorm.k, nsamples)::CuArray{T,2}
    return mvnorm.μ .+ mvnorm.lt_decom * samples
end

function sample(mvnorm::BatchedMultiVariateNormal{Array{T,2}}, nsamples::Int) where {T}
    samples = randn(T, mvnorm.k, nsamples)
    return mvnorm.μ .+ mvnorm.lt_decom * samples
end

function sample(
    mvnorm::BatchedMultiVariateNormal{TrackedArray{T,2,Array{T,2}}},
    nsamples::Int,
) where {T}
    samples = randn(T, mvnorm.k, nsamples)
    return mvnorm.μ .+ mvnorm.lt_decom * samples
end

# Gradient Updates
function update_parameters!(ps, gs, opt)
    for (p, g) in zip(ps, gs)
        length(p) == 0 && continue
        # any(isnan.(g)) && error("NaN Gradient Detected. Terminating...")
        Flux.Optimise.update!(opt, data(p), data(g))
        p.tracker.f = Tracker.Call()
    end
end


text:
module SummaryWrapper

using JavaCall

const AbstractStorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.AbstractStorelessUnivariateStatistic
const AggregatableStatistic = @jimport org.hipparchus.stat.descriptive.AggregatableStatistic
const Class = @jimport java.lang.Class
const DoubleConsumer = @jimport java.util.function.DoubleConsumer
const Iterable = @jimport java.lang.Iterable
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const Product = @jimport org.hipparchus.stat.descriptive.summary.Product
const StorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.StorelessUnivariateStatistic
const Sum = @jimport org.hipparchus.stat.descriptive.summary.Sum
const SumOfLogs = @jimport org.hipparchus.stat.descriptive.summary.SumOfLogs
const SumOfSquares = @jimport org.hipparchus.stat.descriptive.summary.SumOfSquares
const UnivariateStatistic = @jimport org.hipparchus.stat.descriptive.UnivariateStatistic
const WeightedEvaluation = @jimport org.hipparchus.stat.descriptive.WeightedEvaluation

export Product
export Sum
export SumOfLogs
export SumOfSquares
export accept
export aggregate
export and_then
export clear
export copy
export equals
export evaluate
export get_class
export get_n
export get_result
export hash_code
export increment
export increment_all
export notify
export notify_all
export to_string
export wait

include("product.jl")
include("sum.jl")
include("sum_of_logs.jl")
include("sum_of_squares.jl")

end


text:
using ACE
using ACE: alloc_temp, alloc_B, AtomState

species = [:Ti, :Mo]
Bμ0 = ACE.Species1PBasisCtr(species)
Bμ1 = ACE.Species1PBasisNeig(species)
Rn = ACE.Rn1pBasis(ACE.Utils.radial_basis())
Ylm = ACE.Ylm1pBasis(10)
B1p = ACE.Product1pBasis( (Bμ0, Bμ1, Rn, Ylm) )

ACE.symbols(B1p)
ACE.indexrange(B1p)

# TODO : automate this...
ACE.init1pspec!(B1p, maxdeg=5)
length(B1p)

φ = ACE.Invariant()
ace = PIBasis(B1p, ACE.One1pBasis(), 3, 5) |> length #  length ~ 62,000
ace = PIBasis(B1p, ACE.One1pBasis(), 3, 5; property = φ)   # filtered: length ~ 7,000  !!!

# full specification of the 𝑨 Basis
ACE.get_spec(ace)

symace = ACE.SymmetricBasis(ace, φ)

function rand_state()
   𝒓 = ACE.rand_radial(Rn.R) * ACE.Random.rand_sphere()
   μ = rand(species)
   return AtomState(μ, 𝒓)
end

A = alloc_B(B1p)
tmp = alloc_temp(B1p)

Xs = [ rand_state() for _ = 1:30 ]
X0 = rand_state()

ACE.evaluate!(A, tmp, B1p, Xs, X0)
using BenchmarkTools
@btime ACE.evaluate!($A, $tmp, $B1p, $Xs, $X0)

AA = alloc_B(ace)
tmpace = alloc_temp(ace)
ACE.evaluate!(AA, tmpace, ace, Xs, X0)

B = alloc_B(symace)
tmp = alloc_temp(symace)
ACE.evaluate!(B, tmp, symace, Xs, X0)

#---

using StaticArrays, Random

function rand_rot(Xs::AbstractVector{<: AtomState})
   K = (@SMatrix rand(3,3)) .- 0.5
   Q = exp(K - K')
   return [ AtomState(X.mu, Q * X.rr) for X in Xs ]
end

rand_refl(Xs::AbstractVector{<: AtomState}) = (σ = rand([-1, 1]);
      [ AtomState(X.mu, σ * X.rr)  for X in Xs ] )

rand_sym(Xs) = shuffle(rand_refl(rand_rot(Xs)))

for ntest = 1:10
   Xs1 = rand_sym(Xs)
   B1 = alloc_B(symace)
   B1 = ACE.evaluate!(B1, tmp, symace, Xs1, X0)
   @show norm(B1 - B)
   # isapprox(B, B1; rtol = 1e-10)  # use this intead of ≈
end


text:
@testset "Trapping" begin
    # TODO: Remove dummy argument
    function trapkern(x)
        AMDGPU.trap()
        nothing
    end
    function debugtrapkern(x)
        AMDGPU.debugtrap()
        nothing
    end

    iob = IOBuffer()
    AMDGPU.code_gcn(iob, trapkern, Tuple{Int}; kernel=true)
    @test occursin("s_trap 2", String(take!(iob)))
    iob = IOBuffer()
    AMDGPU.code_gcn(iob, debugtrapkern, Tuple{Int}; kernel=true)
    @test occursin("s_trap 3", String(take!(iob)))
end