text:
# ------------------------------------------------------------------------------------------
# # 数据结构
#
# 一旦我们开始一次性处理大量数据，采用数组（array）或者字典（dictionaries）这样的数据结构会比单纯依赖普通变量要方便得多。
#
# 要讲到的数据结构：
# 1. 元组 Tuples
# 2. 字典 Dictionaries
# 3. 数组 Arrays
#
# 总的来说，元组和数组都是元素的有序序列（所以可以直接索引元素）。字典和数组都是可变类型。
# 下面进行详细一点的介绍！
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 元组 Tuples
#
# 通过用小括号 `( )` 将一个有序集合括起来可以创建一个元组。
#
# 语法：
# ```julia
# (item1, item2, ...)
# ```
# ------------------------------------------------------------------------------------------

myfavoriteanimals = ("penguins", "cats", "sugargliders")

# ------------------------------------------------------------------------------------------
# 我们可以在这个元组里进行索引，
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# 但因为元组是不可变对象，我们不可以改变它
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1] = "otters"

# ------------------------------------------------------------------------------------------
# ## 命名元组（Named Tuples）
#
# 顾名思义，命名元组 `NamedTuple` 应该和元组 `Tuple` 差不多，
# 只不过每个元素都有名字！命名元组的创建语法是在元组里使用 `=`：
#
# ```julia
# (name1 = item1, name2 = item2, ...)
# ```
# ------------------------------------------------------------------------------------------

myfavoriteanimals = (bird = "penguins", mammal = "cats", marsupial = "sugargliders")

# ------------------------------------------------------------------------------------------
# 就像元组 `Tuples`，命名元组 `NamedTuples` 是有序的，所以也可以通过索引获得元素：
# ------------------------------------------------------------------------------------------

myfavoriteanimals[1]

# ------------------------------------------------------------------------------------------
# 命名元组还有特殊的通过元素名来访问元素的方法：
# ------------------------------------------------------------------------------------------

myfavoriteanimals.bird

# ------------------------------------------------------------------------------------------
# ## 字典 Dictionaries
#
# 如果你的数据集有映射关系，可以考虑使用字典来存储数据。通过 `Dict()` 函数初始化字典，可以为空或者有键值对。
#
# 语法：
# ```julia
# Dict(key1 => value1, key2 => value2, ...)
# ```
#
# 典型的例子如通讯列表，它就是相关联的姓名和电话号码。
# ------------------------------------------------------------------------------------------

myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")

# ------------------------------------------------------------------------------------------
# 在这个例子中，每对姓名和号码就是一个键值对。我们可以通过对应的键来获取Jenny的电话号码（值）
# ------------------------------------------------------------------------------------------

myphonebook["Jenny"]

# ------------------------------------------------------------------------------------------
# 我们可以在这字典中继续增加条目
# ------------------------------------------------------------------------------------------

myphonebook["Kramer"] = "555-FILK"

# ------------------------------------------------------------------------------------------
# 来瞧瞧现在这个电话本长什么样：
# ------------------------------------------------------------------------------------------

myphonebook

# ------------------------------------------------------------------------------------------
# 使用 `pop!` 函数可以从通讯列表中删除 Kramer，并同时获得他的电话号码
# ------------------------------------------------------------------------------------------

pop!(myphonebook, "Kramer")

myphonebook

# ------------------------------------------------------------------------------------------
# 与元组和数组不同，字典是无序的。所以不可通过索引获得元素。
# ------------------------------------------------------------------------------------------

myphonebook[1]

# ------------------------------------------------------------------------------------------
# 上面这个例子中，Julia 会认为你想获得键为 `1` 的值。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 数组 Arrays
#
# 不同于元组，数组是可变的。不同于字典，数组是有序的。
#
# 使用方括号 `[ ]` 括起元素集合来创建数组。
#
# 语法：
# ```julia
# [item1, item2, ...]
# ```
#
# 举个例子，我们可以创建一个数组来记录我的好友
# ------------------------------------------------------------------------------------------

myfriends = ["Ted", "Robyn", "Barney", "Lily", "Marshall"]

# ------------------------------------------------------------------------------------------
# `Array{String,1}` 中的 `1` 表示这是一个一维的向量。
# 一个 `Array{String,2}` 是个二维矩阵，等等。
#
# `Array{String,1}` 中的 `String` 是每个元素的类型。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 要存储一系列数字的话
# ------------------------------------------------------------------------------------------

fibonacci = [1, 1, 2, 3, 5, 8, 13]

mixture = [1, 1, 2, 3, "Ted", "Robyn"]

# ------------------------------------------------------------------------------------------
# 获得数组后，我们可以通过索引来获得数组中的元素。如，想获得 `myfriends` 中的第三个成员，就可以写成
# ------------------------------------------------------------------------------------------

myfriends[3]

# ------------------------------------------------------------------------------------------
# 我们可以通过索引来修个数组中的元素
# ------------------------------------------------------------------------------------------

myfriends[3] = "Baby Bop"

# ------------------------------------------------------------------------------------------
#  是的，Julia 的索引是从 1 开始的，而不像 Python 是从 0 开始的。
#  这种小问题已经有很多争论了。
#
#  我有个拥有 Solomon 智慧的朋友建议使用从 ½ 开始计数来一劳永逸地解决这个争端。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 函数 `push!` 和 `pop!` 也可以改变数组。
# `push!` 在数组最后添加元素，而 `pop!` 删除数组的最后一个元素。
#
# 给 `fibonnaci` 添加一个数
# ------------------------------------------------------------------------------------------

push!(fibonacci, 21)

# ------------------------------------------------------------------------------------------
# 再删掉它
# ------------------------------------------------------------------------------------------

pop!(fibonacci)

fibonacci

# ------------------------------------------------------------------------------------------
# 到目前为止，我给的例子都是一维的标量，但其实数组可以有任意个维度、可以包含其他数组。
#

# 举个例子，下面这个数组就包含了数组：
# ------------------------------------------------------------------------------------------

favorites = [["koobideh", "chocolate", "eggs"],["penguins", "cats", "sugargliders"]]

numbers = [[1, 2, 3], [4, 5], [6, 7, 8, 9]]

# ------------------------------------------------------------------------------------------
# 下面的例子是包含随机数的二维和三维数组。
# ------------------------------------------------------------------------------------------

rand(4, 3)

rand(4, 3, 2)

# ------------------------------------------------------------------------------------------
# 当你复制数组的时候千万要小心！
# ------------------------------------------------------------------------------------------

fibonacci

somenumbers = fibonacci

somenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# 修改 `somenumbers` 同时也修改了 `fibonacci`！
#
# 在上面那个例子中，实际上我们并没有获得 `fibonacci` 的副本。
# 我们只是创建了一个新的获取 `fibonacci` 所绑定数组的元素的途径。
#
# 如果想获得 `fibonacci` 所绑定数组的副本，可以使用 `copy` 函数。
# ------------------------------------------------------------------------------------------

# 先复原 fibonacci
fibonacci[1] = 1
fibonacci

somemorenumbers = copy(fibonacci)

somemorenumbers[1] = 404

fibonacci

# ------------------------------------------------------------------------------------------
# 上面这个例子中，fibonacci 就没有受影响。
# 因此可以判断 `somemorenumbers` 和 `fibonacci` 所绑定的数组是不同的。
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 练习
#
# ### 3.1
# 用下面的代码创建一个数组 `a_ray`：
#
# ```julia
# a_ray = [1, 2, 3]
# ```
#
# 在这个数组的末尾添加数字 `4` 再删掉它。
# ------------------------------------------------------------------------------------------

# 在此作答

@assert a_ray == [1, 2, 3]

# ------------------------------------------------------------------------------------------
# ### 3.2
# 试着通过下面的代码给 `myphonebook` 添加一个键为 "Emergency" 值为 `string(911)` 条目
# ```julia
# myphonebook["Emergency"] = 911
# ```
#
# 为啥不好使呢？
# ------------------------------------------------------------------------------------------

# 在此作答。更正代码以消除错误
myphonebook["Emergency"] = 911

ans = myphonebook["Emergency"]
@assert  ans == "$ans"

# ------------------------------------------------------------------------------------------
# ### 3.3
# 创建一个名为 `flexible_phonebook` 的字典，其中 Jenny 的号码为整型，而 Ghostbusters 的号码为字符串
#
# ```julia
# (Jenny = 8675309, Ghostbusters = 555-2368)
# ```
# ------------------------------------------------------------------------------------------

flexible_phonebook = Dict()

@assert flexible_phonebook == Dict("Jenny" => 8675309, "Ghostbusters" => "555-2368")

# ------------------------------------------------------------------------------------------
# ### 3.4
# 给`flexible_phonebook`添加一个键为"Emergency"值为整数`911`的条目。
# ------------------------------------------------------------------------------------------

# 在此作答

@assert haskey(flexible_phonebook, "Emergency")

@assert flexible_phonebook["Emergency"] == 911

# ------------------------------------------------------------------------------------------
# ### 3.5
# 思考题：为什么`flexible_phonebook`可以添加值为整型的条目而`myphonebook`不行？
# 该如何初始化`myphonebook`才能添加值为整型的条目呢？
# ------------------------------------------------------------------------------------------

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

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
# plot first 6,000 time steps (before stabilization)
include("../Model_scripts/Functions.jl");
# plot the number of infected & uninfected hosts
plot_population_numbers(
    hcat(labels...),
    "c=4.0\n \nNumber of infected and uninfected hosts",
    5,
    6000,
    "$figure_directory/supplemental/supp_3b.png"
)
# plot the average weighted R0
plot_reproductive_rate(5, 6000, "$figure_directory/supplemental/supp_3d.png")
# plot the average V0
plot_vertical_reproductive_ratio(1, 6000, "$figure_directory/supplemental/supp_3f.png")
# plot the average horizontal transmission
plot_horizontal_and_virulence(1, 6000, "$figure_directory/supplemental/supp_3h.png")
# plot the evenness through time
plot_evenness(1, 6100, "$figure_directory/supplemental/supp_3j.png")

ext: jl


text:
include("../src/Transit.jl")

function run_test(path)
    println("Running $path")
    include(path)
end

run_test("tset.jl")
run_test("decoder.jl")
run_test("encoder.jl")
run_test("roundtrip.jl")
run_test("exemplar.jl")


text:
using Test
using HybridSystems

function test_state_prop(automaton)
    prop = HybridSystems.state_property(automaton, Float64)
    @test typeof(prop) == HybridSystems.state_property_type(typeof(automaton), Float64)
    prop[1] = 0.5
    @test prop[1] == 0.5
    iprop = HybridSystems.typed_map(Int, x -> convert(Int, ceil(x)), prop)
    @test typeof(iprop) == HybridSystems.state_property_type(typeof(automaton), Int)
    @test iprop[1] == 1
end

function test_trans_prop(automaton, t1, t2, t3, t4, t5)
    prop = HybridSystems.transition_property(automaton, Int)
    @test typeof(prop) == HybridSystems.transition_property_type(typeof(automaton), Int)
    prop[t2] = 2
    prop[t1] = 3
    prop[t4] = 4
    prop[t3] = 5
    prop[t5] = 1
    @test prop[t5] == 1
    @test prop[t1] == 3
    @test prop[t3] == 5
    @test prop[t2] == 2
    @test prop[t4] == 4
    fprop = HybridSystems.typed_map(Float64, x -> x / 2, prop)
    @test typeof(fprop) == HybridSystems.transition_property_type(typeof(automaton), Float64)
    @test fprop[t4] == 2.0
    @test fprop[t5] == 1/2
    @test fprop[t2] == 1.0
    @test fprop[t3] == 5/2
    @test fprop[t1] == 3/2
end

@testset "Automaton" begin
    @testset "OneStateAutomaton" begin
        automaton = OneStateAutomaton(5)
        t = map(σ -> begin
            t = HybridSystems.OneStateTransition(σ)
            @test source(automaton, t) == 1
            @test target(automaton, t) == 1
            @test event(automaton, t) == σ
            return t
        end, 1:5)
        test_state_prop(automaton)
        test_trans_prop(automaton, t...)
        @test collect(states(automaton)) == [1]
        @test nstates(automaton) == 1
        @test collect(transitions(automaton)) == t
        @test collect(in_transitions(automaton, 1)) == t
        @test collect(out_transitions(automaton, 1)) == t
        @test has_transition(automaton, 1, 1) == true
        @test length(transitions(automaton)) == 5
        @test ntransitions(automaton) == 5
        @test !has_transition(automaton, HybridSystems.OneStateTransition(0))
        @test  has_transition(automaton, HybridSystems.OneStateTransition(1))
        @test  has_transition(automaton, HybridSystems.OneStateTransition(2))
        @test  has_transition(automaton, HybridSystems.OneStateTransition(3))
        @test  has_transition(automaton, HybridSystems.OneStateTransition(4))
        @test  has_transition(automaton, HybridSystems.OneStateTransition(5))
        @test !has_transition(automaton, HybridSystems.OneStateTransition(6))
        @test !has_transition(automaton, HybridSystems.OneStateTransition(7))
    end

    @testset "LightAutomaton" begin
        automaton = LightAutomaton(2)
        t1 = add_transition!(automaton, 1, 2, 1)
        @test has_transition(automaton, 1, 2)
        @test !has_transition(automaton, 1, 1)
        t2 = add_transition!(automaton, 1, 1, 1)
        t3 = add_transition!(automaton, 2, 1, 2)
        @test !has_transition(automaton, 2, 2)
        t4 = add_transition!(automaton, 2, 2, 1)
        t5 = add_transition!(automaton, 1, 2, 2)
        @test has_transition(automaton, t1)
        @test has_transition(automaton, t5)
        @test has_transition(automaton, t3)
        @test has_transition(automaton, t2)
        @test has_transition(automaton, t4)
        @test has_transition(automaton, 1, 2)
        @test has_transition(automaton, 1, 1)
        @test has_transition(automaton, 2, 1)
        @test has_transition(automaton, 2, 2)
        @test length(out_transitions(automaton, 1)) == 3
        @test length(out_transitions(automaton, 2)) == 2
        @test length(in_transitions(automaton, 1)) == 2
        @test length(in_transitions(automaton, 2)) == 3
        @test eltype(transitions(automaton)) == typeof(first(transitions(automaton)))
        test_state_prop(automaton)
        test_trans_prop(automaton, t1, t2, t3, t4, t5)
        @test ntransitions(automaton) == 5
        @test length(transitions(automaton)) == 5
        rem_transition!(automaton, t5)
        @test has_transition(automaton, t1)
        @test !has_transition(automaton, t5)
        @test has_transition(automaton, t2)
        @test ntransitions(automaton) == 4
        @test length(transitions(automaton)) == 4
        rem_transition!(automaton, t1)
        @test !has_transition(automaton, t1)
        @test !has_transition(automaton, t5)
        @test has_transition(automaton, t2)
        @test ntransitions(automaton) == 3
        @test !has_transition(automaton, 1, 2)
        rem_transition!(automaton, t2)
        @test !has_transition(automaton, t1)
        @test !has_transition(automaton, t5)
        @test !has_transition(automaton, t2)
        @test ntransitions(automaton) == 2
        t1 = add_transition!(automaton, 1, 2, 1)
        @test has_transition(automaton, 1, 2)
        @test ntransitions(automaton) == 3
        t2 = add_transition!(automaton, 1, 1, 1)
        @test ntransitions(automaton) == 4
        t5 = add_transition!(automaton, 1, 2, 2)
        @test ntransitions(automaton) == 5
        t5 = rem_state!(automaton, 2)
        @test ntransitions(automaton) == 1
    end
end


text:
"""
    reduce(step, xf, reducible; [init, simd, basesize]) :: T

Thread-based parallelization of [`foldl`](@ref).  The "bottom"
reduction function `step(::T, ::T) :: T` must be associative and
`init` must be its identity element.

Transducers composing `xf` must be stateless (e.g., [`Map`](@ref),
[`Filter`](@ref), [`Cat`](@ref), etc.) except for [`ScanEmit`](@ref).
Note that [`Scan`](@ref) is not supported (although possible in
theory).  Early termination requires Julia ≥ 1.3.

See [`foldl`](@ref), [`dreduce`](@ref).

# Keyword Arguments
- `basesize::Integer = length(reducible) ÷ nthreads()`: A size of
  chunk in `reducible` that is processed by each worker.  A smaller
  size may be required when:
    * computation time for processing each item fluctuates a lot
    * computation can be terminated by [`reduced`](@ref) or
      transducers using it, such as [`ReduceIf`](@ref)
- For other keyword arguments, see [`foldl`](@ref).
"""
Base.reduce

"""
    mapreduce(xf, step, reducible; init, simd)

!!! warning

    `mapreduce` exists primary for backward compatibility.  It is
    recommended to use `reduce`.

Like [`reduce`](@ref) but `step` is _not_ automatically wrapped by
[`Completing`](@ref).
"""
Base.mapreduce

struct SizedReducible{T} <: Reducible
    reducible::T
    basesize::Int
end

foldable(reducible::SizedReducible) = reducible.reducible

issmall(reducible::SizedReducible) =
    length(reducible.reducible) <= max(reducible.basesize, 1)

function halve(reducible::SizedReducible)
    left, right = halve(reducible.reducible)
    return (
        SizedReducible(left, reducible.basesize),
        SizedReducible(right, reducible.basesize),
    )
end

function halve(arr::AbstractArray)
    # TODO: support "slow" arrays
    mid = length(arr) ÷ 2
    left = @view arr[firstindex(arr):firstindex(arr) - 1 + mid]
    right = @view arr[firstindex(arr) + mid:end]
    return (left, right)
end

struct TaskContext
    listening::Vector{Threads.Atomic{Bool}}
    cancellables::Vector{Threads.Atomic{Bool}}
end

TaskContext() = TaskContext([], [])

function splitcontext(ctx::TaskContext)
    c = Threads.Atomic{Bool}(false)
    return (
        fg = TaskContext(ctx.listening, vcat(ctx.cancellables, c)),
        bg = TaskContext(vcat(ctx.listening, c), ctx.cancellables),
    )
end

function should_abort(ctx::TaskContext)
    for c in ctx.listening
        c[] && return true
    end
    return false
end

function cancel!(ctx::TaskContext)
    for c in ctx.cancellables
        c[] = true
    end
end

function transduce_assoc(
    xform::Transducer, step, init, coll;
    simd::SIMDFlag = Val(false),
    basesize::Integer = length(coll) ÷ Threads.nthreads(),
)
    reducible = SizedReducible(coll, basesize)
    rf = maybe_usesimd(Reduction(xform, step), simd)
    @static if VERSION >= v"1.3-alpha"
        acc = @return_if_reduced _reduce(TaskContext(), rf, init, reducible)
    else
        acc = @return_if_reduced _reduce_threads_for(rf, init, reducible)
    end
    return complete(rf, acc)
end

function _reduce(ctx, rf, init, reducible::Reducible)
    should_abort(ctx) && return init
    if issmall(reducible)
        acc = foldl_nocomplete(rf, _start_init(rf, init), foldable(reducible))
        if acc isa Reduced
            cancel!(ctx)
        end
        return acc
    else
        left, right = halve(reducible)
        fg, bg = splitcontext(ctx)
        task = @spawn _reduce(bg, rf, init, right)
        a0 = _reduce(fg, rf, init, left)
        b0 = fetch(task)
        a = @return_if_reduced a0
        b = @return_if_reduced b0
        should_abort(ctx) && return a  # slight optimization
        return combine(rf, a, b)
    end
end

function _reduce_threads_for(rf, init, reducible::SizedReducible{<:AbstractArray})
    arr = reducible.reducible
    basesize = reducible.basesize
    nthreads = max(
        1,
        basesize <= 1 ? length(arr) : length(arr) ÷ basesize
    )
    if nthreads == 1
        return foldl_nocomplete(rf, _start_init(rf, init), arr)
    else
        w = length(arr) ÷ nthreads
        results = Vector{Any}(undef, nthreads)
        Threads.@threads for i in 1:nthreads
            if i == nthreads
                chunk = @view arr[(i - 1) * w + 1:end]
            else
                chunk = @view arr[(i - 1) * w + 1:i * w]
            end
            results[i] = foldl_nocomplete(rf, _start_init(rf, init), chunk)
        end
        i = findfirst(isreduced, results)
        i === nothing || return results[i]
        # It can be done in `log2(n)` for loops but it's not clear if
        # `combine` is compute-intensive enough so that launching
        # threads is worth enough.  Let's merge the `results`
        # sequentially for now.
        c = foldl(results) do a, b
            combine(rf, a, b)
        end
        return c
    end
end

# AbstractArray for disambiguation
Base.mapreduce(xform::Transducer, step, itr::AbstractArray;
               init = MissingInit(), kwargs...) =
    unreduced(transduce_assoc(xform, step, init, itr; kwargs...))

Base.mapreduce(xform::Transducer, step, itr;
               init = MissingInit(), kwargs...) =
    unreduced(transduce_assoc(xform, step, init, itr; kwargs...))

Base.reduce(step, xform::Transducer, itr; kwargs...) =
    mapreduce(xform, Completing(step), itr; kwargs...)


text:
###
### Janzen-Connell Model
### (c) Daniel Vedder, MIT license
###

struct Species
    id::UInt8 # < 256
    max_age::Int16 # < 32768
    max_size::Int8 # < 128m
    growth_rate::Int8
    seed_production::Int16
    dispersal_distance::Int16
    pathogen_resistance::Float16
end

#The default species
Species(id) = Species(id, 150, 25, 2, 10, 200, 0)

mutable struct Pathogen
    host::UInt8
    infectious::Bool
    infection_rate::Float16
    infection_radius::Int16
    lethality::Float16
end

#The default pathogen
Pathogen(host) = Pathogen(host, false, 0.8, settings["transmission"], 0.05)

mutable struct Tree
    uid::Int64 #XXX should be large enough - hopefully...
    species::Species
    age::Int16
    size::Int8
    mature::Bool
    infection::Union{Pathogen,Nothing}
    position::NamedTuple{(:x, :y), Tuple{Int16,Int16}}
end

# Store and administrate a running counter for tree UID values
let idcounter::Int64 = 0
    global function nextid()
        idcounter += 1
        return idcounter
    end

    global currentcounter() = idcounter
    global resetcounter() = idcounter = 0
end

#Create a seed
Tree(sp, xpos, ypos) = Tree(nextid(), sp, 0, 1, false, nothing, (x=xpos, y=ypos))

"""
Vary a number i randomly  by up to +/- p% (utility function)
"""
function vary(i::Number; p::Int=50)
    i == 0 && return 0
    v = (p/100) * i
    s = i/100
    if isinteger(i)
        s = 1
        v = round(typeof(i), v)
    end
    n = i + rand(-v:s:v)
    return n
end

let species::Vector{Species} = Vector{Species}(undef,settings["species"])

    """
    Initialise all species. If `default` is true, use the standard trait values,
    otherwise create variable species.
    """
    global function createspecies()
        for n in 1:settings["species"]
            s = Species(n)
            if settings["neutral"]
                species[n] = s
            else
                max_age = vary(s.max_age)
                max_size = vary(s.max_size)
                growth_rate = vary(s.growth_rate)
                seed_production = vary(s.seed_production)
                dispersal_distance = vary(s.dispersal_distance)
                pathogen_resistance = vary(s.pathogen_resistance)
                species[n] = Species(n, max_age, max_size, growth_rate,
                                     seed_production, dispersal_distance, pathogen_resistance)
            end
        end
    end

    """
    Return the species object with this ID.
    """
    global function getspecies(id)
        return species[id]
    end

    #XXX Do I need a reset() function?
end

max_stars_repo_name: veddox/jcm

max_issues_repo_name: veddox/jcm

max_forks_repo_name: veddox/jcm


text:
using MemPool
import MemPool: mmwrite, mmread, MMSer

#### DataValueArray

function mmwrite(io::AbstractSerializer, xs::DataValueArray)
    Base.serialize_type(io, MMSer{DataValueArray})

    mmwrite(io, BitArray(xs.isnull))
    mmwrite(io, xs.values)
end

function mmread(::Type{DataValueArray}, io, mmap)
    isnull = deserialize(io)
    vals = deserialize(io)
    DataValueArray(vals, isnull)
end

using PooledArrays

function mmwrite(io::AbstractSerializer, xs::PooledArray)
    Base.serialize_type(io, MMSer{PooledArray})

    mmwrite(io, xs.pool)
    mmwrite(io, xs.refs)
end

function mmread(::Type{PooledArray}, io, mmap)
    pool = deserialize(io)
    refs = deserialize(io)
    PooledArray(PooledArrays.RefArray(refs), pool)
end

# Columns, NDSparse

function mmwrite(io::AbstractSerializer, xs::Columns)
    Base.serialize_type(io, MMSer{Columns})

    if eltype(xs) <: NamedTuple
        fnames = fieldnames(eltype(xs))
    else
        fnames = length(columns(xs))
    end

    serialize(io, fnames)
    for c in columns(xs)
        mmwrite(io, c)
    end
end

function mmread(::Type{Columns}, io, mmap)
    fnames = deserialize(io)
    if isa(fnames, Int)
        cols = [deserialize(io) for i=1:fnames]
        Columns(cols...)
    else
        cols = [deserialize(io) for i=1:length(fnames)]
        Columns(cols...; names=fnames)
    end
end

function mmwrite(io::AbstractSerializer, xs::NDSparse)
    Base.serialize_type(io, MMSer{NDSparse})

    flush!(xs)
    mmwrite(io, xs.index)
    mmwrite(io, xs.data)
end

function mmread(::Type{NDSparse}, io, mmap)
    idx = deserialize(io)
    data = deserialize(io)
    NDSparse(idx, data, presorted=true, copy=false)
end

function mmwrite(io::AbstractSerializer, xs::NextTable)
    Base.serialize_type(io, MMSer{NextTable})

    #flush!(xs)
    mmwrite(io, rows(xs))
    mmwrite(io, xs.pkey)
    mmwrite(io, xs.perms)
end

function mmread(::Type{NextTable}, io, mmap)
    data = deserialize(io)
    pkey = deserialize(io)
    perms = deserialize(io)
    table(data, pkey=pkey, perms=perms, presorted=true, copy=false)
end


text:
# SAM Flags
# =========
#
# Bitwise flags (or FLAG).

for (name, bits, doc) in [
        (:PAIRED,        UInt16(0x001), "the read is paired in sequencing, no matter whether it is mapped in a pair"),
        (:PROPER_PAIR,   UInt16(0x002), "the read is mapped in a proper pair"                                       ),
        (:UNMAP,         UInt16(0x004), "the read itself is unmapped; conflictive with SAM.FLAG_PROPER_PAIR"        ),
        (:MUNMAP,        UInt16(0x008), "the mate is unmapped"                                                      ),
        (:REVERSE,       UInt16(0x010), "the read is mapped to the reverse strand"                                  ),
        (:MREVERSE,      UInt16(0x020), "the mate is mapped to the reverse strand"                                  ),
        (:READ1,         UInt16(0x040), "this is read1"                                                             ),
        (:READ2,         UInt16(0x080), "this is read2"                                                             ),
        (:SECONDARY,     UInt16(0x100), "not primary alignment"                                                     ),
        (:QCFAIL,        UInt16(0x200), "QC failure"                                                                ),
        (:DUP,           UInt16(0x400), "optical or PCR duplicate"                                                  ),
        (:SUPPLEMENTARY, UInt16(0x800), "supplementary alignment"                                                   ),]
    @assert bits isa UInt16
    sym = Symbol("FLAG_", name)
    doc = string(@sprintf("0x%04x: ", bits), doc)
    @eval begin
        @doc $(doc) const $(sym) = $(bits)
    end
end


text:
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/N2Oforcing.jl")

    add_comp!(m, n2oforcing)

    set_param!(m, :n2oforcing, :c_N2Oconcentration, readpagedata(m, "test/validationdata/$valdir/c_n2oconcentration.csv"))
    set_param!(m, :n2oforcing, :c_CH4concentration, readpagedata(m, "test/validationdata/$valdir/c_ch4concentration.csv"))

    ##running Model
    run(m)

    forcing = m[:n2oforcing,:f_N2Oforcing]
    forcing_compare = readpagedata(m, "test/validationdata/$valdir/f_n2oforcing.csv")

    @test forcing ≈ forcing_compare rtol = 1e-3
end


text:
using ProgressMeter
using MATLAB
using Printf
include("PSS_function.jl")

H=range(.11,stop=0.145,length=1001)
@showprogress 1 "Computing..." for Energy in H

mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;pbaspect([1.3 1 1]);"
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/PSS_color_movie/"
h=replace(@sprintf("%.15f",Energy),"."=>"_")
file_name=location*h*".fig"
println(file_name)
mat"axis([ -.025,.025,-.3,.3 ])"
location="/mnt/bdd38f66-9ece-451a-b915-952523c139d2/"

n_iter_Q=51#50
Q_end=.3
n_iter_P=n_iter_Q;#51
P_end=.025
t_end=1e3
global max_hit=1000
#create an empty list to store dH
H_differences=zeros(0)
#julia's version of linspace
ArrP=range(0,stop=P_end,length=n_iter_P)
ArrQ=range(0,stop=Q_end,length=n_iter_Q)

#defining colors for PSS
COLOR=["#393b79" ,"#5254a3","#6b6ecf","#9c9ede" ,"#637939","#8ca252" ,"#b5cf6b" ,"#cedb9c" ,"#8c6d31","#bd9e39" ,"#e7ba52","#e7cb94","#843c39","#ad494a" ,"#d6616b","#e7969c" ,"#7b4173" ,"#a55194","#ce6dbd" ,"#de9ed6"]
#
for j=1:n_iter_P
    Q1,P1,dH=PSS_function(ArrQ[j],0, Energy, t_end, max_hit)
    Q2,P2,dH=PSS_function(-ArrQ[j],0, Energy, t_end, max_hit)
    Q3,P3,dH=PSS_function(0,ArrP[j], Energy, t_end, max_hit)
    Q4,P4,dH=PSS_function(0,-ArrP[j], Energy, t_end, max_hit)
    Q=vcat(Q1,Q2)
    P=vcat(P1,P2)
    current_color=COLOR[mod(2*j,length(COLOR))+1]
    Q_2=vcat(Q3,Q4)
    P_2=vcat(P3,P4)
    current_color_2=COLOR[mod(j+5,length(COLOR))+1]
    mat"plot($P_2,$Q_2,'.','MarkerSize',.1,'color',$current_color_2); hold on;"
    mat"plot($P,$Q,'.','MarkerSize',.1,'color',$current_color); hold on;"
end

n_iter_Q=13;#50
Q_start=-.3
Q_end=.3
n_iter_P=n_iter_Q;#51
P_start=-.025
P_end=.025
#create an empty list to store dH
H_differences=zeros(0)
#julia's version of linspace
ArrP=range(0,stop=P_end,length=n_iter_P)
ArrQ=range(0,stop=Q_end,length=n_iter_Q)

for j=1:n_iter_P
    for k=1:n_iter_Q
        Q1,P1,dH=PSS_function(ArrQ[k], ArrP[j], Energy, t_end, max_hit)
        Q2,P2,dH=PSS_function(-ArrQ[k], ArrP[j], Energy, t_end, max_hit)
        Q3,P3,dH=PSS_function(ArrQ[k], -ArrP[j], Energy, t_end, max_hit)
        Q4,P4,dH=PSS_function(-ArrQ[k], -ArrP[j], Energy, t_end, max_hit)
        Q=vcat(Q1,Q2)
        P=vcat(P1,P2)
        Q_2=vcat(Q3,Q4)
        P_2=vcat(P3,P4)
        current_color=COLOR[mod(j+k,length(COLOR))+1]
        current_color_2=COLOR[mod(j+5,length(COLOR))+1]
        #
        mat"plot($P,$Q,'.','MarkerSize',.1,'color',$current_color); hold on;"
        mat"plot($P_2,$Q_2,'.','MarkerSize',.1,'color',$current_color_2); hold on;"

        # mat"plot($P,$Q,'k.','MarkerSize',.1)"

    end
end
#
mat"savefig($file_name)"
mat"close"
end


text:
using Documenter, DiffEqOperators

include("pages.jl")

makedocs(
    sitename="DiffEqOperators.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[DiffEqOperators],

    format=Documenter.HTML(analytics = "UA-90474609-3",
                           assets=["assets/favicon.ico"],
                           canonical="https://diffeqoperators.sciml.ai/stable/"),

    pages=pages
)

deploydocs(
    repo="github.com/SciML/DiffEqOperators.jl";
    push_preview=true
)


text:
function create_system_data(; with_forecasts = false, time_series_in_memory = false)
    data = IS.SystemData(; time_series_in_memory = time_series_in_memory)

    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(data, component)

    if with_forecasts
        file = joinpath(FORECASTS_DIR, "ComponentsAsColumnsNoTime.json")
        IS.add_forecasts!(IS.InfrastructureSystemsType, data, file)

        forecasts = get_all_forecasts(data)
        @assert length(forecasts) > 0
    end

    return data
end

function create_system_data_shared_forecasts(; time_series_in_memory = false)
    data = IS.SystemData(; time_series_in_memory = time_series_in_memory)

    name1 = "Component1"
    name2 = "Component2"
    component1 = IS.TestComponent(name1, 5)
    component2 = IS.TestComponent(name2, 6)
    IS.add_component!(data, component1)
    IS.add_component!(data, component2)

    ts_data = create_time_series_data()
    forecast = IS.DeterministicInternal("get_val", ts_data)
    IS.add_forecast!(data, component1, forecast, ts_data)
    IS.add_forecast!(data, component2, forecast, ts_data)

    return data
end

function get_all_forecasts(data)
    return collect(IS.iterate_forecasts(data))
end

function create_time_series_data()
    dates = collect(
        Dates.DateTime("1/1/2020 00:00:00", "d/m/y H:M:S"):Dates.Hour(1):Dates.DateTime(
            "1/1/2020 23:00:00",
            "d/m/y H:M:S",
        ),
    )
    data = collect(1:24)
    component_name = "gen"
    ta = TimeSeries.TimeArray(dates, data, [component_name])
    return IS.TimeSeriesData(ta)
end


text:
#using NeuralQuantum
using Test

@testset "NeuralQuantum" begin
    println("Testing machines...")
    @testset "Machines" begin
        include("Machines/test_grad.jl")
        include("Machines/test_ndmcomplex.jl")
        include("Machines/test_lut.jl")
    end

    println("Testing algorithm caches...")
    @testset "Algorithm" begin
        include("Algorithms_caches/Caches_test.jl")
    end

    println("Testing Problems and algorithms...")
    @testset "Problem & algs" begin
        include("Problems/problem_algs.jl")
        include("Problems/problem_operators.jl")
        include("Problems/observables.jl")
        include("Problems/hamiltonian.jl")
    end

    println("Testing States...")
    @testset "States" begin
        include("States/states.jl")
        include("States/modified_states.jl")
    end

    println("Testing Operators...")
    @testset "Operators" begin
        include("operators.jl")
    end

end


text:
#=  Here I replicate a tiny example given in Azevedo and Leshno (2016). Representing each
    student profile as a unit square and highlighting the portions admitted to each school
    gives us insight into the one-to-one correspondence between market-clearing cutoffs
    and stable matches.                                                         =#

using Plots
using Plots.PlotMeasures
using DeferredAcceptance

function argsort(vec::AbstractArray{<:Real, 1})::AbstractArray{<:Real, 1}
    return invperm(sortperm(vec))
end

students = [1 2; 2 1]
students_dist = [1., 1.]

α = 0.75

capacities = rand(2)
capacities *= 2 * α / sum(capacities)

assn, rdist, cutoffs = nonatomicdeferredacceptance(students, students_dist, nothing, capacities,
                   verbose=true, return_cutoffs=true)

p = plot(xlims=(0, 1),
         ylims=(0, 1),
         xlabel="Score at school 1",
         ylabel="Score at school 2",
         layout=2,
         size=(800, 470),
         bottom_margin=70px,
         legend=nothing)

rect(x, y) = Shape([x, 1, 1, x], [y, y, 1, 1])

plot!([Shape([cutoffs[1], 1, 1, cutoffs[1]], [0, 0, 1, 1]),
       Shape([0, cutoffs[1], cutoffs[1], 0], [cutoffs[2], cutoffs[2], 1, 1])],
      color=[:dodgerblue, :olivedrab],
      title=students[:, 1],
      subplot=1)

plot!([Shape([cutoffs[1], 1, 1, cutoffs[1]], [0, 0, cutoffs[2], cutoffs[2]]),
       Shape([0, 1, 1, 0], [cutoffs[2], cutoffs[2], 1, 1])],
      color=[:dodgerblue, :olivedrab],
      stroke=0,
      title=students[:, 2],
      subplot=2)

annotate!(1.08, -.25, text("Correspondence between school cutoffs and\n"*
                           "stable assignments, after Azevedo and Leshno (2016)"), subplot=1)

# Equivalent discrete problem
n = 100     # Number of students in each profile
students = hcat(repeat([1, 2], 1, n), repeat([2, 1], 1, n))
scores = rand(n * 2, 2)
schools = mapslices(argsort, -scores, dims=1)
capacities = round.(Int, capacities .*= n)
@assert sum(capacities) == 3 * n / 2 "oops"

assn, dist = deferredacceptance(students, schools, capacities)

colors = [:dodgerblue, :olivedrab, :crimson]
markers = [:utriangle, :hexagon, :+]

q = plot(xlims=(0, 1),
         ylims=(0, 1),
         xlabel="Score at school 1",
         ylabel="Score at school 2",
         layout=2,
         size=(800, 470),
         bottom_margin=70px,
         legend=nothing)

scatter!(scores[1:n, 1],
         scores[1:n, 2],
         title=students[:, 1],
         color=colors[assn[1:n]],
         markershape=markers[assn[1:n]],
         msw=0,
         marker=markers[assn[1:n]],
         subplot=1)

scatter!(scores[n + 1:end, 1],
         scores[n + 1:end, 2],
         title=students[:, n + 1],
         color=colors[assn[n + 1:end]],
         markershape=markers[assn[n + 1:end]],
         msw=0,
         subplot=2)

annotate!(1.08, -.25, text("Correspondence between school cutoffs\n"*
                           "and stable assignments, discrete form"), subplot=1)

# savefig(p, string("examples/nonatomic/plot-cutoffs.pdf"))
# savefig(p, string("examples/nonatomic/plot-cutoffs.png"))
#
# savefig(q, string("examples/nonatomic/plot-discrete.pdf"))
# savefig(q, string("examples/nonatomic/plot-discrete.png"))

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
using Distributed
using OMEinsumContractionOrders, OMEinsum, CUDA

println("find $(length(devices())) GPU devices")
const procs = addprocs(length(devices())-nprocs()+1)
const gpus = collect(devices())
const process_device_map = Dict(zip(procs, gpus))
@show process_device_map

@everywhere begin  # these packages/functions should be accessible on all processes
    using OMEinsumContractionOrders, OMEinsum, CUDA
    CUDA.allowscalar(false)

    function do_work(f, jobs, results) # define work function everywhere
        while true
            job = take!(jobs)
            @info "running $job on device $(Distributed.myid())"
            res = f(job)
            put!(results, res)
        end
    end
end

"""
    multiprocess_run(func, inputs::AbstractVector)

Run `func` in parallel for a vector f `inputs`.
Returns a vector of results.
"""
function multiprocess_run(func, inputs::AbstractVector{T}) where T
    n = length(inputs)
    jobs = RemoteChannel(()->Channel{T}(n));
    results = RemoteChannel(()->Channel{Any}(n));
    for i in 1:n
        put!(jobs, inputs[i])
    end
    for p in workers() # start tasks on the workers to process requests in parallel
        remote_do(do_work, p, func, jobs, results)
    end
    return Any[take!(results) for i=1:n]
end

"""
    multigpu_einsum(code::SlicedEinsum, xs::AbstractArray...; size_info = nothing, process_device_map::Dict)

Multi-GPU contraction of a sliced einsum specified by `code`.
Each time, the program take the slice and upload them to a specific GPU device and do the contraction.
Other arguments are

* `xs` are input tensors allocated in **main memory**,
* `size_info` specifies extra size information,
* `process_device_map` is a map between processes and GPU devices.
"""
function multigpu_einsum(se::SlicedEinsum{LT,ET}, @nospecialize(xs::AbstractArray...); size_info = nothing, process_device_map::Dict) where {LT, ET}
    length(se.slicing) == 0 && return se.eins(xs...; size_info=size_info)
    size_dict = size_info===nothing ? Dict{OMEinsum.labeltype(se),Int}() : copy(size_info)
    OMEinsum.get_size_dict!(se, xs, size_dict)

    it = OMEinsumContractionOrders.SliceIterator(se, size_dict)
    res = OMEinsum.get_output_array(xs, getindex.(Ref(size_dict), it.iyv))
    eins_sliced = OMEinsumContractionOrders.drop_slicedim(se.eins, se.slicing)
    inputs = collect(enumerate([copy(x) for x in it]))
    @info "start multiple process contraction!"
    results = multiprocess_run(inputs) do (k, slicemap)
        @info "computing slice $k/$(length(it))"
        device!(process_device_map[Distributed.myid()])
        xsi = ntuple(i->CuArray(OMEinsumContractionOrders.take_slice(xs[i], it.ixsv[i], slicemap)), length(xs))
        Array(einsum(eins_sliced, xsi, it.size_dict_sliced))
    end
    # accumulate results to `res`
    for (resi, (k, slicemap)) in zip(results, inputs)
        OMEinsumContractionOrders.fill_slice!(res, it.iyv, resi, slicemap)
    end
    return res
end

# A using case
# ---------------------------------------
using Yao, YaoToEinsum, Yao.EasyBuild

# I. create a quantum circuit
nbit = 20
c = Yao.dispatch!(variational_circuit(nbit, 10), :random)

# II. convert a tensor network
# 1. specify input and output states as product states,
prod_state = Dict(zip(1:nbit, zeros(Int, nbit)))
# 2. convert the circuit to einsum code,
code, xs = YaoToEinsum.yao2einsum(c; initial_state=prod_state, final_state=prod_state)
# 3. optimize the contraction order
size_dict = OMEinsum.get_size_dict(getixsv(code), xs)
slicedcode = optimize_code(code, size_dict, TreeSA(; nslices=5), MergeGreedy())

# III. do the contraction on multiple GPUs in parallel
@info "time/space complexity is $(timespace_complexity(slicedcode, size_dict))"
multigpu_einsum(slicedcode, xs...; process_device_map=process_device_map)


text:
include("utils.jl")

@testset "Graph Manifold" begin
    M = Euclidean(2)
    x = [[1.0, 4.0], [2.0, 5.0], [3.0, 6.0]]
    y = [[4.0, 5.0], [6.0, 7.0], [8.0, 9.0]]
    z = [[6.0, 5.0], [4.0, 3.0], [2.0, 8.0]]
    @testset "Plain Graph" begin
        G = SimpleGraph(3)
        add_edge!(G, 1, 2)
        add_edge!(G, 2, 3)
        N = GraphManifold(G, M, VertexManifold())
        @test manifold_dimension(N) == manifold_dimension(M) * nv(G)
        @test manifold_dimension(GraphManifold(G, M, EdgeManifold())) ==
              manifold_dimension(M) * ne(G)
        @test is_manifold_point(N, x)
        @test !is_manifold_point(N, [x..., [0.0, 0.0]]) # an entry too much
        @test_throws DomainError is_manifold_point(N, [x..., [0.0, 0.0]], true)
        @test is_tangent_vector(N, x, log(N, x, y))
        @test !is_tangent_vector(N, x[1:2], log(N, x, y))
        @test_throws DomainError is_tangent_vector(N, x[1:2], log(N, x, y), true)
        @test !is_tangent_vector(N, x[1:2], log(N, x, y)[1:2])
        @test_throws DomainError is_tangent_vector(N, x, log(N, x, y)[1:2], true)
        @test incident_log(N, x) == [x[2] - x[1], x[1] - x[2] + x[3] - x[2], x[2] - x[3]]

        pts = [x, y, z]
        test_manifold(
            N,
            pts;
            test_representation_size=false,
            test_reverse_diff=false, #VERSION > v"1.2",
        )
        @test sprint(show, "text/plain", N) == """
        GraphManifold
        Graph:
         {3, 2} undirected simple Int64 graph
        Manifold on vertices:
         Euclidean(2; field = ℝ)"""

        NE = GraphManifold(G, M, EdgeManifold())
        @test is_manifold_point(NE, x[1:2])
        @test !is_manifold_point(NE, x) # an entry too much
        @test_throws DomainError is_manifold_point(NE, x, true)
        @test is_tangent_vector(NE, x[1:2], log(N, x, y)[1:2])
        @test !is_tangent_vector(NE, x, log(N, x, y))
        @test_throws DomainError is_tangent_vector(NE, x, log(N, x, y), true)
        @test !is_tangent_vector(N, x[1:2], log(N, x, y))
        @test_throws DomainError is_tangent_vector(NE, x[1:2], log(N, x, y), true)

        test_manifold(
            NE,
            [x[1:2], y[1:2], z[1:2]];
            test_representation_size=false,
            test_reverse_diff=false, #VERSION > v"1.2",
        )
        @test sprint(show, "text/plain", NE) == """
        GraphManifold
        Graph:
         {3, 2} undirected simple Int64 graph
        Manifold on edges:
         Euclidean(2; field = ℝ)"""

        G2 = SimpleDiGraph(3)
        add_edge!(G2, 1, 2)
        add_edge!(G2, 2, 3)
        N2 = GraphManifold(G2, M, VertexManifold())
        @test incident_log(N2, x) == [x[2] - x[1], x[3] - x[2], [0.0, 0.0]]
    end
    @testset "Weighted Graph" begin
        G3 = SimpleWeightedGraph(3)
        add_edge!(G3, 1, 2, 1.5)
        add_edge!(G3, 2, 3, 0.5)
        N3 = GraphManifold(G3, M, VertexManifold())
        @test incident_log(N3, x) == [
            1.5 * (x[2] - x[1]),
            1.5 * (x[1] - x[2]) + 0.5 * (x[3] - x[2]),
            0.5 * (x[2] - x[3]),
        ]
    end
end


text:
# Autogenerated wrapper script for Calcium_jll for i686-linux-musl
export libcalcium

using FLINT_jll
using Arb_jll
using Antic_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Calcium")
JLLWrappers.@declare_library_product(libcalcium, "libcalcium.so.0")
function __init__()
    JLLWrappers.@generate_init_header(FLINT_jll, Arb_jll, Antic_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libcalcium,
        "lib/libcalcium.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module VertexSafeGraphs

import LightGraphs
const LG = LightGraphs

export VSafeGraph

"""
    VSafeGraph{T, G<:LG.AbstractGraph{T}, V<:AbstractVector{Int}}

A `LightGraphs.AbstractGraph` type maintaining vertex numbering, even after vertex removal.
"""
struct VSafeGraph{T, G<:LG.AbstractGraph{T}, V<:AbstractVector{Int}} <: LG.AbstractGraph{T}
    g::G
    deleted_vertices::V
    VSafeGraph(g::G, v::V) where {T, G<:LG.AbstractGraph{T}, V<:AbstractVector{Int}} = new{T, G, V}(g, v)
end

VSafeGraph(g::G) where {G<:LG.AbstractGraph} = VSafeGraph(g, Vector{Int}())
VSafeGraph(nv::Integer) = VSafeGraph(LG.SimpleGraph(nv))

LG.edges(g::VSafeGraph) = LG.edges(g.g)
LG.edgetype(g::VSafeGraph) = LG.edgetype(g.g)

LG.is_directed(g::VSafeGraph) = LG.is_directed(g.g)
LG.is_directed(::Type{VSafeGraph}) = false
LG.is_directed(::Type{<:VSafeGraph{T,G}}) where {T, G} = LG.is_directed(G)

LG.ne(g::VSafeGraph) = LG.ne(g.g)
LG.nv(g::VSafeGraph) = LG.nv(g.g) - length(g.deleted_vertices)
LG.vertices(g::VSafeGraph) = (v for v in LG.vertices(g.g) if !(v in g.deleted_vertices))

LG.has_vertex(g::VSafeGraph, v) = LG.has_vertex(g.g, v) && !(v in g.deleted_vertices)

LG.has_edge(g::VSafeGraph, src, dst) = LG.has_edge(g.g, src, dst)
LG.has_edge(g::VSafeGraph, edge::LG.AbstractEdge) = LG.has_edge(g.g, LG.src(edge), LG.dst(edge))

LG.add_vertex!(g::VSafeGraph) = LG.add_vertex!(g.g)

LG.rem_edge!(g::VSafeGraph, v1, v2) = LG.rem_edge!(g.g, v1, v2)

Base.copy(g::VSafeGraph) = VSafeGraph(copy(g.g), copy(g.deleted_vertices))

function LG.outneighbors(g::VSafeGraph, v)
    if LG.has_vertex(g, v)
        LG.outneighbors(g.g, v)
    else
        throw(ArgumentError("$v is not a valid vertex in graph."))
    end
end

function LG.inneighbors(g::VSafeGraph, v)
    if LG.has_vertex(g, v)
        LG.inneighbors(g.g, v)
    else
        throw(ArgumentError("$v is not a valid vertex in graph."))
    end
end

function LG.add_edge!(g::VSafeGraph, v1, v2)
    if !LG.has_vertex(g, v1) || !LG.has_vertex(g, v2)
        return false
    end
    return LG.add_edge!(g.g, v1, v2)
end

LG.add_edge!(g::VSafeGraph, edge::LG.AbstractEdge) = LG.add_edge!(g, LG.src(edge), LG.dst(edge))

function LG.rem_vertex!(g::VSafeGraph, v1)
    if !LG.has_vertex(g, v1) || v1 in g.deleted_vertices
        return false
    end
    for v2 in LG.outneighbors(g, v1)
        LG.rem_edge!(g, v1, v2)
    end
    for v2 in LG.inneighbors(g, v1)
        LG.rem_edge!(g, v2, v1)
    end
    push!(g.deleted_vertices, v1)
    return true
end

end # module


text:
module Tests

include("utilities.jl")

const linear_tests = Dict{String, Function}()
const soc_tests = Dict{String, Function}()
const sd_tests = Dict{String, Function}()

include("term.jl")
include("term_fixed.jl")
include("quartic_constant.jl")
include("quadratic.jl")
include("quartic_ideal.jl")
include("univariate_sum.jl")
include("rearrangement.jl")
include("choi.jl")
include("simple_matrix.jl")
include("horn.jl")
include("concave_then_convex_cubic.jl")
include("lyapunov_switched_system.jl")
include("motzkin.jl")
include("BPT12e399.jl")
include("maxcut.jl")
include("chebyshev.jl")
include("quartic_comparison.jl")
include("options_pricing.jl")
include("sosdemo5.jl")
include("sosdemo9.jl")
include("sosdemo10.jl")

@test_suite linear
@test_suite soc
@test_suite sd

end


text:
const NUCLEOTIDE_BITS = Dict('A' => unsigned(0),
                       'C' => unsigned(1),
                       'G' => unsigned(2),
                       'T' => unsigned(3))

const KmerType = Array{UInt32, 1}

"""
    kmer_count(str::String, k::Int)

Count kmers of size `k` in string, return array with bins from "A...A" to "T...T"
and value in each bin corresponding to number of occurences of that kmer within `str`.
"""
function kmer_count(str::String, k::Int)
    # TODO: could directly encode `str` as 2-bit BioSequence
    bins = zeros(eltype(KmerType), 4^k)

    mask = unsigned(4^k - 1)  # all ones
    kmer = unsigned(0)
    for c in str[1:k-1]
        kmer = (kmer << 2) + NUCLEOTIDE_BITS[c]
    end
    for c in str[k:end]
        kmer = ((kmer << 2) & mask) + NUCLEOTIDE_BITS[c]
        bins[kmer + 1] += 1
    end
    return bins
end

const AA_DICT = Dict('M'=>11,'I'=>8,'Y'=>20,'L'=>10,'*'=>0,'F'=>5,'Q'=>14,
    'D'=>3,'V'=>18,'E'=>4,'T'=>17,'H'=>7,'P'=>13,'G'=>6,'N'=>12,'K'=>9,'C'=>2,
    'R'=>15,'W'=>19,'A'=>1,'S'=>16,'X'=>22)

"""
    sparse_aa_kmer_count(str::String, k::Int)

Counts amino acid kmers in string (in all reading frames).
k = length of kmer in amino acids.
This is sparse and kinda slow.
"""
function sparse_aa_kmer_count(str::String, k::Int)
    bins = spzeros(UInt64, 21^k)
    kmer = unsigned(0)
    mask = unsigned(21) ^ k
    straas = generate_aa_seqs(str)
    for straa in straas
        for aa in straa[1:k-1]
            kmer = kmer * 21 + AA_DICT[aa]
        end
        for aa in straa[k:end]
            kmer = ((kmer * 21) % mask) + AA_DICT[aa]
            bins[kmer+1] += 1
        end
    end
    return bins
end

"""
    corrected_kmer_dist(kmers1::Array, kmers2::Array; k = nothing)

Compute distance between kmer vectors that is corrected towards edit distance for small differences.
The default for `k` is `k = Int(log(4, length(kmers1)))`
"""
function corrected_kmer_dist(kmers1::Array, kmers2::Array; k = nothing)
    if k == nothing
        k = Int(log(4, length(kmers1)))
    end
    return sqeuclidean(kmers1, kmers2)/ (k*(sum(kmers1) + sum(kmers2)))
end

"""
    corrected_kmer_dist(k::Int)

Returns function that computes distance between `k`-mer vectors
that is corrected towards edit distance for small differences.
"""
function corrected_kmer_dist(k::Int)
    function kdist(kmers1, kmers2)
        return corrected_kmer_dist(kmers1, kmers2, k=k)
    end
    return kdist
end
"""
    corrected_kmer_dist_full(kmers1::Array, kmers2::Array; k = 6)
"""
function corrected_kmer_dist_full(kmers1::Array, kmers2::Array; k = 6)
    return sqeuclidean(kmers1, kmers2)/(2*k)
end


text:
#
# Advent of Code
# Day 5: How About a Nice Game of Chess?
#

# Part 1 completes in: 12.37s
# Part 2 completes in: 31.63s
#
# In Python 2.7.10 (32-bit)
# Part 1 takes ~23.54s
# Part 2 takes ~63.59s
#
# In Python 2.7.12 (64-bit)
# Part 1 takes ~20.71s
# Part 2 takes ~57.06s
#
# In Python 3.5.2 (64-bit)
# Part 1 takes ~28.07s
# Part 2 takes ~75.11s
#
# Using DV's C++ implementation,
# Part 1 takes 66.097s / 38.344s with /Ox
# Part 2 takes 187.31s / 106.322s with /Ox
#
# Using MZ's Python implementation (2.7.10, 32-bit)
# Part 1 takes ~27.52s
# Part 2 takes ~76.02s
#
#
using Nettle

println("Advent of Code")
println("Day 5: How About a Nice Game of Chess?")
println("")
println("Part 1")

#DOOR_ID = "abc"
DOOR_ID = "ojvtpuvg"

notFound = true
counter = 0
foundCounter = 0
password = ""

tic()

while foundCounter < 8
  testHash = hexdigest("md5",DOOR_ID * string(counter))
  if testHash[1:5] == "00000"
    foundCounter += 1
    password = password * string(testHash[6])
    println("$password")
  end

  counter += 1
end

toc()

println("The part 1 password for $DOOR_ID is $password ($counter iterations).")

# part 2
println("\nPart2")

notFound = true
counter = 0
foundCounter = 0
password = "--------"

tic()

while foundCounter < 8
  testHash = hexdigest("md5",DOOR_ID * string(counter))
  if testHash[1:5] == "00000"
    if isnumber(testHash[6])
      position = parse(string(testHash[6])) + 1
      if position <= 8
        if password[position] == '-'
          password = password[1:position-1] * string(testHash[7]) * password[position+1:end]
          foundCounter += 1
          println("$password")
        end
      end
    end
  end
  counter += 1
end
toc()

println("The part 2 password for $DOOR_ID is $password ($counter iterations).")


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03td(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    NMAX = 100
    LDA  = NMAX
    LDG  = NMAX
    LDRES  = NMAX
    LDU1 = NMAX
    LDU2 = NMAX
    LDWORK = 8*NMAX
    SELECT = Array{Bool,1}(undef, NMAX)
    LOWER = Array{Bool,1}(undef, NMAX)
    A = Array{Float64,2}(undef, LDA,NMAX)
    G = Array{Float64,2}(undef, LDG,NMAX)
    U1 = Array{Float64,2}(undef, LDU1,NMAX)
    U2 = Array{Float64,2}(undef, LDU2,NMAX)
    WR = Array{Float64,1}(undef, NMAX)
    WI = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    RES = Array{Float64,2}(undef, LDRES,NMAX)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    TYP =  vs[2][1]
    COMPU =  vs[3][1]
if ( N<=0 || N>NMAX )
else

    vs = String[]
    _isz = N
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    SELECT[1:_isz] .= parsex.(Bool, vs)

    vs = String[]
    _isz = N
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    LOWER[1:_isz] .= parsex.(Bool, vs)

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       G[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
if ( LSAME( COMPU, 'U' ) )

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U1[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U2[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
end # if
# interp call 1
    select = convert(Vector{BlasBool},SELECT)
    lower = convert(Vector{BlasBool},LOWER)
    M, INFO = SLICOT.mb03td!(TYP, COMPU, select, lower, N, A, G, U1, U2, WR, WI, LDWORK)
    println(io, "M = $M")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
if ( LSAME( COMPU, 'U' ) )
# interp output 1
    println(io, "U1:")
    _nc = N
    _nr = N
    show(io, "text/plain", U1[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "U2:")
    _nc = N
    _nr = N
    show(io, "text/plain", U2[1:_nr,1:_nc])
    println(io)

# unable to translate write statement:
#  write MA02JD( .FALSE., .FALSE., N, U1, LDU1, U2, LDU2, RES, LDRES )
end # if
# interp output 3
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

    Gfull = zeros(N,N)
if ( LSAME( TYP, 'S' ) )
    # interp output 4
    for i in 1:N
        for j in 1:i-1
            Gfull[i,j] = -G[j,i]
        end
        for j in i+1:N
            Gfull[i,j] = G[i,j]
        end
    end
else
    for i in 1:N
        for j in 1:i-1
            Gfull[i,j] = G[j,i]
        end
        for j in i:N
            Gfull[i,j] = G[i,j]
        end
    end
end # if
    println(io, "G:")
    _nc = N
    _nr = N
    show(io, "text/plain", Gfull[1:_nr,1:_nc])
    println(io)

end # if
end # if
    close(f)
end # run_mb03td()


text:
module Invariants

using InlineTest
using TextWrap: wrap
import AbstractTrees
import InvariantsCore
import InvariantsCore: AbstractInvariant, InvariantList, AnyInvariant, AllInvariant,
                       AsMarkdown, errormessage, satisfies, invariant, check, check_throw,
                       title, description, showdescription, showtitle

include("wrapio.jl")
include("tree.jl")
include("show.jl")

include("invariants/hasmethod.jl")

function exampleinvariant(symbol = :n)
    return Invariant("`$symbol` is positive",
        description = "The number `$symbol` should be larger than `0`.") do x
        if !(x isa Number)
            return "`$symbol` has type $(typeof(x)), but it should be a `Number` type."
        else
            x > 0 && return nothing
            return "`$symbol` is not a positive number, got value `$x`. Please pass a number larger than 0."
        end
    end
end

export invariant, check, check_throw

end


text:
# benchmark of deque

using DataStructures

# push_back

function batch_pushback!(v::Container, n::Int, e::T) where {Container,T}
    for i = 1 : n
        push!(v, e)
    end
end

v = Int[]
q = deque(Int)

batch_pushback!(v, 10, 0)
t1 = @elapsed batch_pushback!(v, 10^7, 0)

batch_pushback!(q, 10, 0)
t2 = @elapsed batch_pushback!(q, 10^7, 0)

println("push back 10^7 integers:")
@printf("    Vector:   elapsed = %8.4fs\n", t1)
@printf("    Deque:    elapsed = %8.4fs\n", t2)

# push_front

function batch_pushfront!(v::Container, n::Int, e::T) where {Container,T}
    for i = 1 : n
        unshift!(v, e)
    end
end

v = Int[]
q = deque(Int)

batch_pushfront!(v, 10, 0)
t1 = @elapsed batch_pushfront!(v, 10^7, 0)

batch_pushfront!(q, 10, 0)
t2 = @elapsed batch_pushfront!(q, 10^7, 0)

println("push front 10^7 integers:")
@printf("    Vector:   elapsed = %8.4fs\n", t1)
@printf("    Deque:    elapsed = %8.4fs\n", t2)

# traverse

function traverse(container)
    for e in container
    end
end

traverse(v)
t1 = @elapsed traverse(v)

traverse(q)
t2 = @elapsed traverse(q)

println("traverse 10^7 integers:")
@printf("    Vector:   elapsed = %8.4fs\n", t1)
@printf("    Deque:    elapsed = %8.4fs\n", t2)


text:
using SolveDSGE

filename = "model1a.txt"
path = joinpath(@__DIR__,filename)

process_model(path)

processed_filename = "model1a_processed.txt"
processed_path =  joinpath(@__DIR__,processed_filename)

dsge = retrieve_processed_model(processed_path)

x = [0.05, 21.4, 1.9, 2.5]

tol = 1e-8
maxiters = 1000
ss_obj = compute_steady_state(dsge,x,tol,maxiters)
ss = ss_obj.zero

N   = PerturbationScheme(ss,1.0,"first")
NN  = PerturbationScheme(ss,1.0,"second")
NNN = PerturbationScheme(ss,1.0,"third")

soln_fo = solve_model(dsge,N)
soln_so = solve_model(dsge,NN)
soln_to = solve_model(dsge,NNN)

P = ChebyshevSchemeStoch(ss,chebyshev_nodes,[21,21],9,3,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
PP = ChebyshevSchemeStoch(ss,chebyshev_nodes,[21,21],9,4,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
PPP = ChebyshevSchemeStoch(ss,chebyshev_nodes,[71,71],9,6,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)

soln_nla = solve_model(dsge,P)
soln_nlb = solve_model(dsge,soln_fo,PP)
soln_nlc = solve_model(dsge,soln_so,PPP)
soln_nld = solve_model(dsge,soln_to,P)
soln_nle = solve_model(dsge,soln_nld,PP)
soln_nlf = solve_model(dsge,soln_nle,PPP)

L = SmolyakSchemeStoch(ss,chebyshev_gauss_lobatto,9,4,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
LL = SmolyakSchemeStoch(ss,clenshaw_curtis_equidistant,9,4,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
LLL = SmolyakSchemeStoch(ss,chebyshev_gauss_lobatto,9,5,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)

soln_nlg = solve_model(dsge,L)
soln_nlh = solve_model(dsge,soln_to,LL)
soln_nli = solve_model(dsge,soln_nlh,LLL)
soln_nlj = solve_model(dsge,soln_nlf,LL)

M = PiecewiseLinearSchemeStoch(ss,[21,21],9,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)
MM = PiecewiseLinearSchemeStoch(ss,[31,31],9,[0.0960769 26.0; -0.0960769 8.0],tol,1e-6,maxiters)

soln_nlk = solve_model(dsge,M)
soln_nll = solve_model(dsge,soln_to,MM)
soln_nlm = solve_model(dsge,soln_nlf,MM)
soln_nln = solve_model(dsge,soln_nlj,MM)
soln_nlo = solve_model(dsge,soln_nln,MM)

simulated_data1 = simulate(soln_nlf,ss[1:2],100000)
pos_imps, neg_imps = impulses(soln_nlf,50,[1],10000)

nodesf, f = approximate_density(simulated_data1[3,:],10,1.95,1.99)
nodesF, F = approximate_distribution(simulated_data1[3,:],10,1.95,1.99)


text:
using AugmentedGPLikelihoods
using AugmentedGPLikelihoods.SpecialDistributions
using AugmentedGPLikelihoods.TestUtils: test_auglik
const AGPL = AugmentedGPLikelihoods
using Distributions
using GPLikelihoods
using GPLikelihoods.TestInterface: test_interface
using MeasureBase
using LogExpFunctions
using Random
using SpecialFunctions
using Test
using TupleVectors

@testset "AugmentedGPLikelihoods.jl" begin
    @info "Testing likelihoods"
    @testset "Likelihoods" begin
        include("likelihoods/bernoulli.jl")
        include("likelihoods/laplace.jl")
        include("likelihoods/negativebinomial.jl")
        include("likelihoods/poisson.jl")
        include("likelihoods/studentt.jl")
    end

    @info "Testing SpecialDistributions"
    @testset "SpecialDistributions" begin
        include("SpecialDistributions/ntdist.jl")
        include("SpecialDistributions/polyagamma.jl")
        include("SpecialDistributions/polyagammapoisson.jl")
    end

    include("utils.jl")
end


text:
using StatsBase
using Statistics
using Distributions
using Random

# N.B. See Mimi/WIP/load_empirical_dist.jl for helper functions.

struct EmpiricalDistribution{T}
    values::Vector{T}
    weights::ProbabilityWeights
    dist::Distribution

    # Create an empirical distribution from a vector of values and an optional
    # vector of probabilities for each value. If not provided, the values are
    # assumed to be equally likely.
    # N.B. This doesn't copy the values vector, so caller must, if required
    function EmpiricalDistribution(values::Vector{T}, probs::Union{Nothing, Vector{Float64}}=nothing) where T
        n = length(values)
        if probs === nothing
            probs = Vector{Float64}(n)
            probs[:] = 1/n
        elseif length(probs) != n
            error("Vectors of values and probabilities must be equal lengths")
        end

        # If the probs don't exactly sum to 1, but are close, tweak
        # the final value so the sum is 1.
        total = sum(probs)
        if 0 < (1 - total) < 1E-5
            probs[end] += 1 - total
        end

        weights = ProbabilityWeights(probs)

        return new{T}(values, weights, Categorical(probs))
    end
end

#
# Delegate a few functions that we require in our application.
# No need to be exhaustive here.
#
function Statistics.mean(d::EmpiricalDistribution)
    return mean(d.values, d.weights)
end

function Statistics.std(d::EmpiricalDistribution)
    return std(d.values, d.weights, corrected=true)
end

function Statistics.var(d::EmpiricalDistribution)
    return var(d.values, d.weights, corrected=true)
end

function Statistics.quantile(d::EmpiricalDistribution, args...)
    indices = quantile.(d.dist, args...)
    return d.values[indices]
end

function Statistics.rand(d::EmpiricalDistribution, args::Vararg{Integer,N}) where {N}
    indices = rand(d.dist, args...)
    return d.values[indices]
end

function Random.rand!(d::EmpiricalDistribution, args::Vararg{Integer,N}) where {N}
    indices = rand!(d.dist, args...)
    return d.values[indices]
end


text:
##file to include in ComplexModeler.jl which contains all files associated to Matpower

include("RTEMatpowerGenerator.jl")
include("MatpowerLoad.jl")
include("MatpowerROPFShunt.jl")
include("MatpowerLine_π.jl")
include("MatpowerVolt.jl")
#include("read_matpower_simple_RTE.jl")
include("read_matpower_simple_MAT.jl")


text:
ofeltype(x, y) = convert(float(eltype(x)), y)

"""
    reduce_nodes(aggr, g, x)

For a batched graph `g`, return the graph-wise aggregation of the node
features `x`. The aggregation operator `aggr` can be `+`, `mean`, `max`, or `min`.
The returned array will have last dimension `g.num_graphs`.
"""
function reduce_nodes(aggr, g::GNNGraph, x)
    @assert size(x)[end] == g.num_nodes
    indexes = graph_indicator(g)
    return NNlib.scatter(aggr, x, indexes)
end

"""
    reduce_edges(aggr, g, e)

For a batched graph `g`, return the graph-wise aggregation of the edge
features `e`. The aggregation operator `aggr` can be `+`, `mean`, `max`, or `min`.
The returned array will have last dimension `g.num_graphs`.
"""
function reduce_edges(aggr, g::GNNGraph, e)
    @assert size(e)[end] == g.num_edges
    s, t = edge_index(g)
    indexes = graph_indicator(g)[s]
    return NNlib.scatter(aggr, e, indexes)
end

"""
    softmax_nodes(g, x)

Graph-wise softmax of the node features `x`.
"""
function softmax_nodes(g::GNNGraph, x)
    @assert size(x)[end] == g.num_nodes
    gi = graph_indicator(g)
    max_ = gather(scatter(max, x, gi), gi)
    num = exp.(x .- max_)
    den = reduce_nodes(+, g, num)
    den = gather(den, gi)
    return num ./ den
end

"""
    softmax_edges(g, e)

Graph-wise softmax of the edge features `e`.
"""
function softmax_edges(g::GNNGraph, e)
    @assert size(e)[end] == g.num_edges
    gi = graph_indicator(g, edges=true)
    max_ = gather(scatter(max, e, gi), gi)
    num = exp.(e .- max_)
    den = reduce_edges(+, g, num)
    den = gather(den, gi)
    return num ./ den
end

@doc raw"""
    softmax_edge_neighbors(g, e)

Softmax over each node's neighborhood of the edge features `e`.

```math
\mathbf{e}'_{j\to i} = \frac{e^{\mathbf{e}_{j\to i}}}
                    {\sum_{j'\in N(i)} e^{\mathbf{e}_{j'\to i}}}.
```
"""
function softmax_edge_neighbors(g::GNNGraph, e)
    @assert size(e)[end] == g.num_edges
    s, t = edge_index(g)
    max_ = gather(scatter(max, e, t), t)
    num = exp.(e .- max_)
    den = gather(scatter(+, num, t), t)
    return num ./ den
end

"""
    broadcast_nodes(g, x)

Graph-wise broadcast array `x` of size `(*, g.num_graphs)`
to size `(*, g.num_nodes)`.
"""
function broadcast_nodes(g::GNNGraph, x)
    @assert size(x)[end] == g.num_graphs
    gi = graph_indicator(g)
    return gather(x, gi)
end

"""
    broadcast_edges(g, x)

Graph-wise broadcast array `x` of size `(*, g.num_graphs)`
to size `(*, g.num_edges)`.
"""
function broadcast_edges(g::GNNGraph, x)
    @assert size(x)[end] == g.num_graphs
    gi = graph_indicator(g, edges=true)
    return gather(x, gi)
end


text:
using ProximalSEM
using Documenter

DocMeta.setdocmeta!(ProximalSEM, :DocTestSetup, :(using ProximalSEM); recursive=true)

makedocs(;
    modules=[ProximalSEM],
    authors="Maximilian S. Ernst, Aaron Peikert",
    repo="https://github.com/StructuralEquationModels/ProximalSEM.jl/blob/{commit}{path}#{line}",
    sitename="ProximalSEM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://StructuralEquationModels.github.io/ProximalSEM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/StructuralEquationModels/ProximalSEM.jl",
    devbranch="main",
)


text:
export slice

function slice end

function (sim::Simulation{Dim})(::typeof(slice), dir::Direction{D}, p::Real) where {D,Dim}

dirs = directions(sim.grid)
maskdir =  dirs .!= fill(dir,Dim)

ext =  extent(sim.grid)[maskdir]
spac=  spacing(sim.grid)[maskdir]

sim_wg = Simulation(λ₀ = sim.λ₀ ;  grid = Grid(extent = ext , spacing =  spac ))

ind = round(Int,p / spacing(sim.grid,dir))

sim_wg.J_x = selectdim(sim.J_x, D , ind)
sim_wg.J_y = selectdim(sim.J_y, D , ind)
sim_wg.J_z = selectdim(sim.J_z, D , ind)

sim_wg.Q = selectdim(sim.Q, D , ind)
sim_wg.S_x = selectdim(sim.S_x, D , ind)
sim_wg.S_y = selectdim(sim.S_y, D , ind)
sim_wg.S_z = selectdim(sim.S_z, D , ind)

sim_wg.ϵᵣ_xx = selectdim(sim.ϵᵣ_xx, D , ind)
sim_wg.ϵᵣ_yy = selectdim(sim.ϵᵣ_yy, D , ind)
sim_wg.ϵᵣ_zz = selectdim(sim.ϵᵣ_zz, D , ind)

sim_wg.μᵣ_xx = selectdim(sim.μᵣ_xx, D , ind)
sim_wg.μᵣ_yy = selectdim(sim.μᵣ_yy, D , ind)
sim_wg.μᵣ_zz = selectdim(sim.μᵣ_zz, D , ind)

sim_wg.e⁻ⁱᵏᴸ = sim.e⁻ⁱᵏᴸ[maskdir]

return sim_wg

end


text:
using FileIO
import ParallelTrainingNN
using Plots
using StatPlots
import CSV
using DataFrames
using LaTeXStrings
pgfplots()
nn = ParallelTrainingNN

# Import data
df = CSV.read("processed_data/training_set.csv");
ti, ui, yi = Vector(df[:time]), Vector(df[:input]), Vector(df[:output]);
Ts = (ti[2] - ti[1])/3600;
identification_data = nn.IdData(yi, ui, Ts);
df = CSV.read("processed_data/test_set.csv");
tv, uv, yv = Vector(df[:time]), Vector(df[:input]), Vector(df[:output]);
Ts = (tv[2] - tv[1])/3600;
validation_data = nn.IdData(yv, uv, Ts);

# Define neural network model
mdl = nn.FeedforwardNetwork(2, 1, [10])
yterms = [[1]];
uterms = [[1]];
mdl = nn.learn_normalization(mdl, yterms, uterms, identification_data);
freerun = nn.free_run_simulation(mdl, yterms, uterms, validation_data);
# Get data
ys = nn.get_slice(freerun, validation_data.y);

# Compute MSE
seeds = 1:100
shoot_len_list = [3]
θ_est = Matrix{Vector{Float64}}(length(seeds),  length(shoot_len_list))
MSE = Matrix{Float64}(length(seeds), length(shoot_len_list))
execution_time = Matrix{Float64}(length(seeds), length(shoot_len_list))
nfev = Matrix{Float64}(length(seeds), length(shoot_len_list))
res = Matrix(length(seeds), length(shoot_len_list))
for (i, s) in enumerate(seeds)
    for (j, shoot_len) in enumerate(shoot_len_list)
        file_name = "solutions/sol"*string(s)*"_"*string(shoot_len)*".jld2"
        res[i, j] = try
            load(file_name, "res")
        catch
            false
        end
        if isa(res[i, j], Dict)
            θ_est[i, j] = res[i, j]["x"][1:41]
            execution_time[i, j] = res[i, j]["execution_time"]
            nfev[i, j] = res[i, j]["nfev"]
            MSE[i, j] = mean((ys - nn.predict(freerun, θ_est[i, j])).^2 )
        end
    end
end

# Read MSE file
ptch = CSV.read("MSE_parallel_considered_harmfull.csv")

# Empirical cumulative distribution function
prob = 0.01:0.01:1
mse_min = 3e2
mse_max = 2e5
a = 0.05
n = length(seeds)
c = sqrt(1/(2n)*log(2/(1-a)))
ARXmse = 1489.5
# add some points at the extremity to make it easier to vizualize
#prob = vcat(0, prob, 1)
oness = ones(length(prob))
p2 = plot(xlim=[mse_min, mse_max], ylims = [0, 1], yticks=0:0.1:1,
          xlabel="MSE", legend=:bottomright)
cdf = sort(ptch[:series_parallel])
plot!(p2, vcat(mse_min, cdf, mse_max),
      vcat(0, prob, 1),
      label= "NN ARX",
      xscale=:log10,
      color=:green)
plot!(p2, cdf, [prob, prob], fillrange=[max.(prob-c, 0), min.(prob+c, 1)],
     fillalpha=0.1, label= "", xscale=:log10, color=:green)
cdf = sort(ptch[:parallel])
plot!(p2,vcat(mse_min, cdf, mse_max),
      vcat(0, prob, 1),
      label= "NN OE - SS",
     color=:blue)
plot!(p2, cdf, [prob, prob], fillrange=[max.(prob-c, 0), min.(prob+c, 1)],
      fillalpha=0.1, label= "", xscale=:log10, color=:blue)
plot!(p2, cdf[cdf .> ARXmse],
      prob[cdf .> ARXmse],
      label= "",
      xscale=:log10,
      marker=:circle,
      color=:blue)
cdf = sort(MSE[:, 1])
plot!(p2, vcat(mse_min, cdf, mse_max),
      vcat(0, prob, 1),
      label= "NN OE - MS", xscale=:log10,
      color=:red)
plot!(p2, cdf, [prob, prob], fillrange=[max.(prob-c, 0), min.(prob+c, 1)],
      fillalpha=0.1, label= "", xscale=:log10, color=:red)
plot!(p2, cdf[cdf .> ARXmse],
      prob[cdf .> ARXmse],
      label= "",
      xscale=:log10,
      marker=:circle,
      color=:red)
plot!(p2, 1489.5*oness, prob, color=:black, lw=1, s=:dash,
      label= "Linear ARX")

savefig("neuralnet_mse.tex")

println(mean(execution_time[:, 1]));
println(mean(nfev[:, 1]));


text:
# Create STO Problem
function stoproblem(
  x0::Array{Float64,1},                 # Initial State
  A::Array{Float64,3};                  # Linear Dynamics
  ngrid::Int64=2,                       # Number of Linearization points in the fixed grid (2 for linear case. We do not need them by default)
  t0::Float64=0.0,                      # Initial Time
  tf::Float64=1.0,                      # Final Time
  Q::Array{Float64, 2}=emptyfmat,       # Cost Matrix
  E::Array{Float64, 2}=emptyfmat,      # Final Cost Matrix
  lb::Array{Float64, 1}=emptyfvec,      # Lower Bound on Intervals
  ub::Array{Float64, 1}=emptyfvec,      # Upper Bound on Intervals
  tau0ws::Array{Float64,1}=emptyfvec,   # Warm Start tau vector
  solver::MathProgBase.AbstractMathProgSolver=Ipopt.IpoptSolver(maxiter = 30))

  # Get Dimensions
  nx = size(A, 1)      # State Dimension
  N = size(A, 3) - 1   # Get number of switches

  # Adjust variables which have not been initalized
  if isempty(Q)
    Q = eye(nx)
  end

  if isempty(E)
    E = zeros(nx, nx)
  end

  if isempty(lb)
    lb = zeros(N+1)
  end

  if isempty(ub)
    ub = Inf*ones(N+1)
  end

  if isempty(tau0ws)
    tau0ws = collect(linspace(t0, tf, N+2))
    tau0ws = tau0ws[2:end-1]
  end

  # Define warm starting delta0
  delta0ws = tau2delta(tau0ws, t0, tf)

  # Create Discretization Grid
  tgrid = collect(linspace(t0, tf, ngrid))

  # Initialize time vectors
  tvec = Array(Float64, N + ngrid)    # Complete grid
  tauIdx = Array(Int, N + 2)      # Indeces of switching times in the complete grid
  tgridIdx = Array(Int, ngrid)      # Indeces of switching times in the complete grid
  deltacomplete = Array(Float64, N + ngrid - 1)   # Intervals over the whole grid

  ### Initialize NLP Evaluator
  # Preallocate arrays
  deltafun_prev = Array(Float64, N+1)
  deltagrad_prev = Array(Float64, N+1)
  deltahess_prev = Array(Float64, N+1)
  xpts = Array(Float64, nx, N+ngrid); xpts[:, 1] = x0   # Set Initial State
  expMat = Array(Float64, nx, nx, N+ngrid-1)
  Phi = Array(Float64, nx, nx, N+2, N+2)
  M = Array(Float64, nx, nx, N+ngrid-1)
  S = Array(Float64, nx, nx, N+ngrid)
  C = Array(Float64, nx, nx, N+1)

  # Decompose Dynamics Matrices
  V = Array(Complex{Float64}, 2*nx, 2*nx, N+1)
  invV = Array(Complex{Float64}, 2*nx, 2*nx, N+1)
  D =  Array(Complex{Float64}, 2*nx, N+1)
  isDiag = Array(Bool, N+1)

  for i = 1:N+1
    D[:, i], V[:, :, i] = eig([-A[:, :, i]'  Q;
                             zeros(nx, nx) A[:, :, i]])
    if cond(V[:, :, i]) == Inf  # Non diagonalizable matrix
      isDiag[i] = false
    else
      invV[:, :, i] = inv(V[:, :, i])
      isDiag[i] = true
    end
  end

  # Construct Matrix of Indeces for lower triangular Matrix (Hessian)
  IndTril = find(tril(ones(N+1, N+1)))
  Itril, Jtril, _ = findnz(tril(ones(N+1, N+1)))

  # Construct Constraints Matrix (Vector)
  Ag = ones(1, N+1)
  Ig, Jg, Vg = findnz(Ag)
  bg = [tf]   # Only one constraints for the sum of the switching intervals

  # Initialize objective evaluator
  obj = Array(Float64, 0)
  deltaval = Array(Float64, N+1, 0)
  nobjeval = 0                           # Number of objective function evaluations
  ngradeval = 0                           # Number of gradient evaluations
  nhesseval = 0                           # Number of hessian evaluations

  # Construct NLP evaluator
  STOev = linSTOev(x0, nx, A, N, t0, tf, tf, Q, E, ngrid, tgrid, tvec, tauIdx, tgridIdx, deltacomplete, V, invV, D, isDiag, IndTril, Itril, Jtril, Ag, Ig, Jg, Vg, bg, deltafun_prev, deltagrad_prev, deltahess_prev, xpts, expMat, Phi, M, S, C,
  obj, deltaval, nobjeval, ngradeval, nhesseval)

  # Generate Model
  m = MathProgBase.NonlinearModel(solver)

  ### Load NLP Program into the model
  MathProgBase.loadproblem!(m, N+1, length(bg), lb, ub, bg, bg, :Min, STOev)

  # Propagate dynamic for new switching times
  propagateDynamics!(STOev, delta0ws)

  ### Add Warm Starting Point
  MathProgBase.setwarmstart!(m, delta0ws)

  # Create STO
  STOproblem = linSTO(m, STOev, delta0ws)

  return STOproblem  # Return STO

end

# Create STO Problem for Nonlinear STO
function stoproblem(
  x0::Array{Float64,1},             # Initial State
  nonlin_dyn::Function,             # Nonlinear Dynamics
  nonlin_dyn_deriv::Function,       # Nonlinear Dynamics Derivative
  uvec::Array{Float64, 2};          # Vector of integer Inputs per switching combination
  ngrid::Int64=10,                  # Number of Linearization points in the fixed grid
  t0::Float64=0.0,                  # Initial Time
  tf::Float64=1.0,                  # Final Time
  Q::Array{Float64, 2}=emptyfmat,   # Cost Matrix
  E::Array{Float64, 2}=emptyfmat,  # Final Cost Matrix
  lb::Array{Float64, 1}=emptyfvec,  # Lower Bound on Intervals
  ub::Array{Float64, 1}=emptyfvec,  # Upper Bound on Intervals
  tau0ws::Array{Float64,1}=emptyfvec, # Warm Start tau vector
  solver::MathProgBase.AbstractMathProgSolver=Ipopt.IpoptSolver(maxiter = 30))

  # Get Dimensions
  nx = length(x0)         # State Dimension

  # Generate artificial switches repeating the input sequence
  N = size(uvec, 2) - 1   # Get total number of switches

  # Adjust variables which have not been initalized
  if isempty(Q)
    Q = eye(nx)
  end

  if isempty(E)
    E = zeros(nx, nx)
  end

  if isempty(lb)
    lb = zeros(N+1)
  end

  if isempty(ub)
    ub = Inf*ones(N+1)
  end

  if isempty(tau0ws)
    tau0ws = collect(linspace(t0, tf, N+2))  # Currently counting tau_0 and tau_{N+1}. They are removed below after xpts initialization.
    tau0ws = tau0ws[2:end-1]  # Include only switching instants

  end

  # Create Discretization grid
  tgrid = collect(linspace(t0, tf, ngrid))

  # Initialize time vectors
  tvec = Array(Float64, N + ngrid)    # Complete grid
  tauIdx = Array(Int, N + 2)      # Indeces of switching times in the complete grid
  tgridIdx = Array(Int, ngrid)      # Indeces of switching times in the complete grid
  deltacomplete = Array(Float64, N + ngrid - 1)   # Intervals over the whole grid

  ### Initialize NLP Evaluator
  # Extend System State To get Affine to Linear Dynamics
  nx = nx + 1

  # Extend Initial State and Cost Matrix
  x0 = [x0; 1]
  spz = Array(Float64, 1,1); spz[1,1] = 0.0; spz = sparse(spz)  # Sparse scalar
  Q = full(blkdiag(sparse(Q), spz))
  E = full(blkdiag(sparse(E), spz))

  # Define Required Matrices and Switching Instants
  A = Array(Float64, nx, nx, N+ngrid-1)

  # Initialize Switching Instants
  xpts = Array(Float64, nx, N+ngrid)
  xpts[:, 1] = x0

  # Define warm starting delta0ws from tau0ws
  delta0ws = tau2delta(tau0ws, t0, tf)

  deltafun_prev = Array(Float64, N+1)
  deltagrad_prev = Array(Float64, N+1)
  deltahess_prev = Array(Float64, N+1)
  expMat = Array(Float64, nx, nx, N+ngrid-1)
  Phi = Array(Float64, nx, nx, N+2, N+2)
  M = Array(Float64, nx, nx, N+ngrid-1)
  S = Array(Float64, nx, nx, N+ngrid)
  C = Array(Float64, nx, nx, N+1)  # Only at sw times (including initial time)

  # Construct Matrix of Indeces for lower triangular Matrix (Hessian)
  IndTril = find(tril(ones(N+1, N+1)))
  Itril, Jtril, _ = findnz(tril(ones(N+1, N+1)))

  # Construct Constraints Matrix (Vector)
  Ag = ones(1, N+1)
  Ig, Jg, Vg = findnz(Ag)
  bg = [tf]   # Only one constraints for the sum of the switching intervals

  # Initialize objective evaluator
  obj = Array(Float64, 0)
  deltaval = Array(Float64, N+1, 0)
  nobjeval = 0                           # Number of objective function evaluations
  ngradeval = 0                           # Number of gradient evaluations
  nhesseval = 0                           # Number of hessian evaluations

  # Construct NLPEvaluator
  STOev = nlinSTOev(x0, nx, A, N, t0, tf, tf, Q, E, uvec, ngrid, tgrid, tvec, tauIdx, tgridIdx, deltacomplete, nonlin_dyn, nonlin_dyn_deriv, IndTril, Itril, Jtril, Ag, Ig, Jg, Vg, bg, deltafun_prev, deltagrad_prev, deltahess_prev, xpts, expMat, Phi, M, S, C, obj, deltaval, nobjeval, ngradeval, nhesseval)

  # Propagate Dynamics to compute matrix exponentials and states at the switching times
  propagateDynamics!(STOev, delta0ws)

  # Generate Model
  m = MathProgBase.NonlinearModel(solver)

  ### Load NLP Program into the model
  MathProgBase.loadproblem!(m, N+1, length(bg), lb, ub, bg, bg, :Min, STOev)

  ### Add Warm Starting Point for the solver
  MathProgBase.setwarmstart!(m, delta0ws)

  # Create STO
  STOproblem = nlinSTO(m, STOev, delta0ws)

  return STOproblem  # Return STO

end

# Solve Optimization for Linear System
function solve!(m::STO)

  # Perform STO
  m.soltime = @elapsed MathProgBase.optimize!(m.model)
  m.stat = MathProgBase.status(m.model)
  m.delta = MathProgBase.getsolution(m.model)
  m.tau, _ = delta2tau(m.delta, m.STOev.t0)
  m.objval = MathProgBase.getobjval(m.model)

end

"Set warm starting point"
function setwarmstart!(m::STO, tau0ws::Array{Float64,1})

  # Define warm starting delta0
  delta0ws = tau2delta(tau0ws, m.STOev.t0, m.STOev.tf)

  # Set Warm Starting Point for Nonlinear Solver
  MathProgBase.setwarmstart!(m.model, delta0ws)

  # Reset Cost Function Iterates
  m.STOev.obj = Array(Float64,0)
  m.STOev.deltaval = Array(Float64, m.STOev.N+1, 0)

end

"Set initial state x0"
function setx0!(m::STO, x0::Array{Float64,1})

  # Define initial point
  m.STOev.x0 = x0

end

# Return Variables from STO
gettau(m::STO) = m.tau
getdelta(m::STO) = m.delta
getobjval(m::STO) = m.objval
getstat(m::STO) = m.stat
getsoltime(m::STO) = m.soltime
getnobjeval(m::STO) = m.STOev.nobjeval
getngradeval(m::STO) = m.STOev.ngradeval
getnhesseval(m::STO) = m.STOev.nhesseval

# Convert from Switching Times to Intervals
function tau2delta(tau::Array{Float64, 1}, t0::Float64, tf::Float64)

  # Extend vector of switching times with initial and final time
  tauexp = [t0; tau; tf]  # Extend tau vector to simplify numbering

  # Return Vector of differences
  return  diff(tauexp)

end

# Convert from Intervals to Switching Times
function delta2tau(delta::Array{Float64, 1}, t0::Float64)

  # Define tau vector
  tau = Array(Float64, length(delta)-1)

  # Initialize first tau element
  tau[1] = t0 + delta[1]

  # Iteratively add all deltas to create switching times vector
  for i = 2:length(tau)
    tau[i] = tau[i-1] + delta[i]
  end

  tfdelta = tau[end] + delta[end]

  return tau, tfdelta

end


text:
#X(u)_function.jl

"""
Implementation of FHIP's susceptibility formula χ(Ω) for the Optical Absorption of Polarons and uses parts of Devresse's, Sitter's & Goovaerts paper that expands upon calculations in the FHIP paper (particularly for ℜχ(Ω)).

See FHIP 1962, equations (35a), (44), (47a, b, c):
https://link.aps.org/doi/10.1103/PhysRev.127.1004

and Devreese's et al. 1971 paper, equations (15), (16), and Appendix A:
https://link.aps.org/doi/10.1103/PhysRevB.5.2367.
"""

using QuadGK # Have to include this as this script somehow doesn't pick it up in the main PolaronMakie.jl module.

"""
----------------------------------------------------------------------
Simple implementation for finite temperatures. Useful for most values of arguments.
----------------------------------------------------------------------
"""

"""
ℜχ(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64)

    Calculate the real part of χ(Ω) (equation (35a) in FHIP) for a given frequency Ω. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. β is the reduced thermodynamical beta = ħω/kT with k being Boltzman's constant. This attempts to just solve the integral directly using adaptive gaussian quadrature methods.
"""
function ℜχ(Ω, β, α, v, w)
    integrand(x) = (1 - cos(Ω * x)) * ℑS(x, 0.0, v, w, β, α)
    return QuadGK.quadgk(x -> integrand(x), 0.0, Inf; order = 14)[1]
end

"""
ℑχ(Ω::Float64, β::Float64, α::Float64, v::Float64, w::Float64)

    Calculate the imaginary part of χ(Ω) (equations (35a) and (44) in FHIP) for a given frequency Ω. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. β is the reduced thermodynamical beta = ħω/kT with k being Boltzman's constant. This attempts to just solve the integral directly using adaptive gaussian quadrature methods.
"""
function ℑχ(Ω, β, α, v, w)
    integrand(x) = -sin(Ω * x) * ℑS(x, 0.0, v, w, β, α)
    return QuadGK.quadgk(x -> integrand(x), 0.0, Inf; order = 14)[1]
end

# Create a complex number from the real and imaginary parts of χ(Ω).

χ(Ω, β, α, v, w) = ℜχ(Ω, β, α, v, w) + 1im * ℑχ(Ω, β, α, v, w)

"""
----------------------------------------------------------------------
Low temperature implementation as β -> ∞.
----------------------------------------------------------------------
"""

"""
ℜχ_0(Ω::Float64, α::Float64, v::Float64, w::Float64, N::Int)

    Calculate the real part of χ(Ω) in a zero temperature approximation (equation (16) in Devreese's et al.) for a given frequency Ω. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. N is the upper limit of a sum that is analytically infinite, however most of the information is encapuslated by N < 50.
"""
function ℜχ_0(Ω, α, v, w, N = 10) # Large N for greater accuracy.
    R = (v^2 - w^2) / (w^2 * v)
    integrand(x, n) =
        ((n + 1 / 2) * x^(n - 1 / 2) * exp(-R * x) - R * x^(n + 1 / 2) * exp(-R * x)) *
        log(abs((1 + n * v + x)^2 / (Ω^2 - (1 + n * v + x)^2))^(1 / 2))
    total_sum = 0.0
    for n = 0:Int(N)
        total_sum += -(1 / gamma(n + 3 / 2)) * quadgk(x -> integrand(x, n), 0, Inf)[1]
    end
    return total_sum
end

"""
ℑχ_0(Ω::Float64, α::Float64, v::Float64, w::Float64, N::Int)

    Calculate the imaginary part of χ(Ω) in a zero temperature approximation (equation (16) in Devreese's et al.) for a given frequency Ω. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. N is the upper limit of a sum that is analytically infinite, however most of the information is encapuslated by N < 50.
"""
function ℑχ_0(Ω, α, v, w, N = 10)

    R = (v^2 - w^2) / (w^2 * v)
    coefficient = 2 / 3 * α * (v / w)^3

    total_sum = 0.0
    for n = 0:Int(N)
        total_sum +=
            -coefficient * sqrt(π) / (gamma(-n - 1 / 2) * gamma(n + 1)) * (-2 * R)^n /
            BigCombinatorics.doublefactorial(2 * n + 1) *
            abs(Ω - 1 - n * v)^(n + 1 / 2) *
            exp(-R * abs(Ω - 1 - n * v)) *
            (1 + sign(Ω - 1 - n * v))
    end
    return total_sum
end

"""
----------------------------------------------------------------------
Finite temperature implementation for ℑχ using BesselK functions.
----------------------------------------------------------------------

Details come from Appendix A of Devreese's et al. paper, although they take the limit of β -> ∞, whereas here we do not. Likewise, we provide proper treatment of any expansions without immediate approximatons (i.e. in Appendix A they ignore the even part of the cos^n(vx) expansion without any apparent justification).
"""

"""
ℑχ_expansion(Ω::Float64, α::Float64, v::Float64, w::Float64, N::Int)

    Calculate the imaginary part of χ(Ω) at finite temperature using an infinite expansion of BesselK functions (see Appendix A in Devreese's et al.), for a given frequency Ω. v and w are the variational Polaron parameters that minimise the free energy, for the supplied α Frohlich coupling. N is the upper limit of a sum that is analytically infinite, however most of the information is encapuslated by N < 50.

    Function can either implement the BesselK functions to FLoat64 precision using the SpecialFunctions.jl package (which eventually Overflows), or to arbitrary precision with the ArbNumerics.jl package. Change comments where appropriate. Set the precision of the ArbReal type with 'setworkingprecision(ArbReal, bits=128)'. This needs to be the same precision as BigFloat, set with 'setprecision(BigFloat, 128)'.
"""
function ℑχ_expansion(Ω, β, α, v, w, N = 10) # larger N increases accuracy.

    # Set arguments to BigFloat precision. Without this the calculations break down due to large values of hyperbolic functions.

    Ω = BigFloat(Ω)
    β = BigFloat(β)
    α = BigFloat(α)
    v = BigFloat(v)
    w = BigFloat(w)

    # Initialise constants.

    R = (v^2 - w^2) / (w^2 * v)
    a = sqrt(β^2 / 4 + R * β * coth(β * v / 2))
    b = R * β / sinh(β * v / 2)

    # Coefficient independent of summation variables n & k.

    coefficient =
        2 * α / (3 * sqrt(π)) * (v / w)^3 * β^(3 / 2) * sinh(β * Ω / 2) / sinh(β / 2)

    # Initialise total value of double summation as a BigFloat.

    total_sum = BigFloat(0.0)

    for n = 0:Int(N)   # Infinite summation from the Binomial expansion of (x^2 + a^2 - b * cos(v * x))^(-3/2). |b * cos(v * x) / (x^2 + a^2)| < 1 for v > 0 and β > 0.

        # Coefficient that depends on n summation variable.

        n_coefficient =
            -2 * beta(-1 / 2 - n, n + 1)^(-1) * (-1)^n * b^n * (4 * a)^(-n - 1) * sqrt(π) /
            gamma(n + 3 / 2)

        # Summation over expansion of (cos(v * x))^n splits into sums over odd n and even n.

        if isodd(n) # For odd n.

            # Finite sum over k for odd n from (cos(v * x))^n expansion.

            for k = -1:Int((n - 1) / 2)

                # Coefficient that depends on k summation variable.

                k_coefficient_inverse = ((n + 1) * beta(n - k + 1, k + 1))

                # Arguments of resultant cosines obtained from cos^n(vx) * cos(x) * cos(Ωx).

                y = [
                    (Ω + 1 + v * (n - 2 * k)),
                    (Ω + 1 - v * (n - 2 * k)),
                    (Ω - 1 + v * (n - 2 * k)),
                    (Ω - 1 - v * (n - 2 * k)),
                ]

                for x in y # Iterate over cosine arguments for brevity.

                    """
                    For low-ish values of β (~ < 40) and Ω (~ < 50). Integral over cos(yx) / (x^2 + a^2)^(n + 3/2) can be identified with BesselK functions.
                    """

                    total_sum +=
                        coefficient *
                        n_coefficient *
                        besselk(Int(n + 1), abs(Float64(x)) * Float64(a)) *
                        abs(x)^(n + 1) / k_coefficient_inverse

                    """
                    For larger values of β (~ > 40) and Ω (~ > 50). This is an arbitrary precision implementation of BesselK functions by ArbNumerics.jl that is not available in the SpecialFunctions.jl package. Comment above and uncomment below for higher precision calculations.
                    """

                    # Note: May be type issues somewhere resulting in Any types and slowing down code.

                    # total_sum += coefficient * n_coefficient * k_coefficient * besselk(ArbReal(n + 1), ArbReal(abs(x) * a)) * abs(x)^(n + 1)
                end
            end
        end

        if iseven(n) # For even n.

            # Finite sum over k for even n from (cos(v * x))^n expansion.

            for k = -1:Int(n / 2 - 1)

                # Coefficient that depends on k summation variable.

                k_coefficient_inverse = ((n + 1) * beta(n - k + 1, k + 1))

                # Arguments of resultant cosines obtained from cos^n(vx) * cos(x) * cos(Ωx).

                y = [
                    (Ω + 1 + v * (n - 2 * k)),
                    (Ω + 1 - v * (n - 2 * k)),
                    (Ω - 1 + v * (n - 2 * k)),
                    (Ω - 1 - v * (n - 2 * k)),
                ]

                for x in y # Iterate over cosine arguments for brevity.

                    """
                    For low-ish values of β (~ < 40) and Ω (~ < 50). Integral over cos(yx) / (x^2 + a^2)^(n + 3/2) can be identified with BesselK functions.
                    """

                    total_sum +=
                        coefficient *
                        n_coefficient *
                        besselk(Int(n + 1), abs(Float64(x)) * Float64(a)) *
                        abs(x)^(n + 1) / k_coefficient_inverse

                    """
                    For larger values of β (~ > 40) and Ω (~ > 50). This is an arbitrary precision implementation of BesselK functions by ArbNumerics.jl that is not available in the SpecialFunctions.jl package. Comment above and uncomment below for higher precision calculations.
                    """

                    # Note: May be type issues somewhere resulting in Any types and slowing down code.

                    # total_sum += coefficient * n_coefficient * k_coefficient * besselk(ArbReal(n + 1), ArbReal(abs(x) * a)) * abs(x)^(n + 1)
                end

                # Coefficient that depends on k summation variable. This term is an extra term that arises for the even part of the cos^n(vx) expansion.

                k_even_coefficient_inverse = ((n + 1) * beta(n / 2 + 1, n / 2 + 1))

                # Arguments of resultant cosines obtained from just cos(x) * cos(Ωx).

                y_even = [(Ω + 1), (Ω - 1)]

                for x in y_even # Iterate over cosine arguments for brevity.

                    """
                    For low-ish values of β (~ < 40) and Ω (~ < 50). Integral over cos(y_even x) / (x^2 + a^2)^(n + 3/2) can be identified with BesselK functions.
                    """

                    total_sum +=
                        coefficient *
                        n_coefficient *
                        besselk(Int(n + 1), abs(Float64(x)) * Float64(a)) *
                        abs(x)^(n + 1) / k_even_coefficient_inverse

                    """
                    For larger values of β (~ > 40) and Ω (~ > 50). This is an arbitrary precision implementation of BesselK functions by ArbNumerics.jl that is not available in the SpecialFunctions.jl package. Comment above and uncomment below for higher precision calculations.
                    """

                    # Note: May be type issues somewhere resulting in Any types and slowing down code.

                    # total_sum += coefficient * n_coefficient * k_even_coefficient * besselk(ArbReal(n + 1), ArbReal(abs(x) * a)) * abs(x)^(n + 1)
                end
            end
        end
    end

    # Return final value obtained from double summation.
    return total_sum
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

download("https://www.ngdc.noaa.gov/geomag/data/poles/NP.xy", "NP.xy")


text:
module PointsetRegistration

export
    # accelerated_routines.jl
    neg_XtRY, wYPXt, cost, partial_cost,

    # clustering.jl
    kernel_kmeans, max_k_cut, factor_k_cut,

    # grid.jl
    grid, grid_coords,

    # GPM.jl
    GPM_orthogonal, GPM_permutation,

    # density_grid.jl
    parallel_density_grid, density_grid, density_grid_stats,

    # register.jl
    register, exact_register, gaussian_register, register_to_reference,

    # synchronization.jl
    parallel_pairwise_register, stiefel_to_permutations,

    # pointsets.jl
    trim, random_perturb!, random_permute!, random_rotate!, random_rotate_reflect!, pointsets_transform, pointsets_rmsd,

    # group_representations
    irrep_O3!, irrep_O3

include("accelerated_routines.jl")
include("clustering.jl")
include("grid.jl")
include("GPM.jl")
include("density_grid.jl")
include("register.jl")
include("synchronization.jl")
include("pointsets.jl")
include("irreps_O3.jl")

end


text:
using Luxor
using OffsetArrays

function make_drawing(width, height, img_path, bkg_color, origin_p)
    d = Drawing(width, height, img_path)
    background(bkg_color)
    origin(origin_p)
    return d
end

width = 500
height = 500
path = "voronoi.png"
color = "black"

my_draw = make_drawing(width, height, path, color, O)

slope(x₁, y₁, x₂, y₂) = (y₂ - y₁) / (x₂ - x₁)
slope(P1, P2) = (P2.y - P1.y) / (P2.x - P1.x)
midpoint(x₁, y₁, x₂, y₂) = ((x₁ + x₂) / 2, (y₁ + y₂) / 2)
midpoint(P1, P2) = ((P1.x + P2.x) / 2, (P2.y + P2.y) / 2)
point_slope(x, x₁, y₁, m) = (x, m * (x - x₁) + y₁)

mat = zeros(width, height)

p1 = (x = 125, y = 125)
p2 = (x = 375, y = 375)

sp = slope(p1, p2)
mp = midpoint(p1, p2)

d(x1, y1, x2, y2, x, y) = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1)

h1_plane = []
h2_plane = []
h3_plane = []

for x in 1:500, y in 1:500
	if d(500, 0, mp[1], mp[2], x, y) < 0
		push!(h1_plane, Point(x, y))
	elseif d(500, 0, mp[1], mp[2], x, y) > 0
		push!(h2_plane, Point(x, y))
	else
		push!(h3_plane, Point(x, y))
	end
end

# for p in vals
	# mat[Int(p.x), Int(p.y)] = 1
# end

# function half_plane(p1, p2, width, height)

    # mp = midpoint(p1.x, p1.y, p2.x, p2.y) |> Point
    # sp = slope(p1.x, p1.y, p2.x, p2.y)
    # first_point = point_slope(-250, mp.x, mp.y, sp == 0 ? 0 : -1 / sp) |> Point
    # second_point = point_slope(250, mp.x, mp.y, sp == 0 ? 0 : -1 / sp) |> Point

    # half_1 = [first_point, second_point, Point(width / 2, height / 2)]
    # half_2 = [first_point, second_point, Point(-width / 2, height / 2)]

    # return half_1, half_2
# end

setcolor(1, 0, 0, 1)
for c in h1_plane
	circle(c, 1; action = :stroke)
end

setcolor(0, 0, 1, 1)
for c in h2_plane
	circle(c, 1; action = :stroke)
end

setcolor(1, 1, 0, 1)
for c in h3_plane
	circle(c, 1; action = :stroke)
end

setcolor(1, 0, 0, 1)
circle(p1.x, p1.y, 15, action = :fill)

setcolor(0, 0, 1, 1)
circle(p2.x, p2.y, 15, action = :fill)

finish()

max_stars_repo_name: TheCedarPrince/DataExplorers

max_issues_repo_name: TheCedarPrince/DataExplorers

max_forks_repo_name: TheCedarPrince/DataExplorers


text:
@testset "7.reverse-integer.jl" begin
    @test reverse_integer(Int32(123)) == 321
    @test reverse_integer(Int32(-123)) == -321
end


text:
module JungleHelperTheoStatueGate

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/TheoStatueGate" TheoStatueGate(x::Integer, y::Integer, height::Integer=48, sprite::String="default")

const placements = Ahorn.PlacementDict(
    "Green Crystal Gate (Jungle Helper)" => Ahorn.EntityPlacement(
        TheoStatueGate,
        "point",
        Dict{String, Any}(
            "sprite" => "theo"
        )
    )
)

const textures = String["default", "mirror", "theo"]

Ahorn.editingOptions(entity::TheoStatueGate) = Dict{String, Any}(
    "sprite" => textures
)

function Ahorn.selection(entity::TheoStatueGate)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x - 4, y, 15, height)
end

const sprites = Dict{String, String}(
    "default" => "objects/door/TempleDoor00",
    "mirror" => "objects/door/TempleDoorB00",
    "theo" => "objects/door/TempleDoorC00"
)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::TheoStatueGate, room::Maple.Room)
    sprite = get(entity.data, "sprite", "default")

    if haskey(sprites, sprite)
        Ahorn.drawImage(ctx, sprites[sprite], -4, 0)
    end
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper