text:
"""
    ro = IdOffsetRange(r::AbstractUnitRange, offset=0)

Construct an "identity offset range". Numerically, `collect(ro) == collect(r) .+ offset`,
with the additional property that `axes(ro, 1) = axes(r, 1) .+ offset`.
When `r` starts at 1, then `ro[i] == i` and even `ro[ro] == ro`,
i.e., it's the "identity," which is the origin of the "Id" in `IdOffsetRange`.

# Examples

The most common case is shifting a range that starts at 1 (either `1:n` or `Base.OneTo(n)`):
```jldoctest; setup=:(import OffsetArrays)
julia> ro = OffsetArrays.IdOffsetRange(1:3, -2)
OffsetArrays.IdOffsetRange(-1:1)

julia> axes(ro, 1)
OffsetArrays.IdOffsetRange(-1:1)

julia> ro[-1]
-1

julia> ro[3]
ERROR: BoundsError: attempt to access 3-element UnitRange{$Int} at index [5]
```

If the range doesn't start at 1, the values may be different from the indices:
```jldoctest; setup=:(import OffsetArrays)
julia> ro = OffsetArrays.IdOffsetRange(11:13, -2)
OffsetArrays.IdOffsetRange(9:11)

julia> axes(ro, 1)     # 11:13 is indexed by 1:3, and the offset is also applied to the axes
OffsetArrays.IdOffsetRange(-1:1)

julia> ro[-1]
9

julia> ro[3]
ERROR: BoundsError: attempt to access 3-element UnitRange{$Int} at index [5]
```

# Extended help

Construction/coercion preserves the (shifted) values of the input range, but may modify
the indexes if required by the specified types. For example,

    r = OffsetArrays.IdOffsetRange{Int,UnitRange{Int}}(3:4)

has `r[1] == 3` and `r[2] == 4`, whereas

    r = OffsetArrays.IdOffsetRange{Int,Base.OneTo{Int}}(3:4)

has `r[3] == 3` and `r[4] == 4`, and `r[1]` would throw a `BoundsError`.
In this latter case, a shift in the axes was needed because `Base.OneTo` ranges
must start with value 1.

!!! warning

    In the future, *conversion* will preserve both the values and
    the indices, throwing an error when this is not achievable. For instance,

        r = convert(OffsetArrays.IdOffsetRange{Int,UnitRange{Int}}, 3:4)

    has `r[1] == 3` and `r[2] == 4` and would satisfy `r == 3:4`, whereas

    ```julia
    julia> convert(OffsetArrays.IdOffsetRange{Int,Base.OneTo{Int}}, 3:4)    # future behavior, not present behavior
    ERROR: ArgumentError: first element must be 1, got 3
    ```

    where the error will arise because the result could not have the same axes as the input.

    An important corollary is that `typeof(r1)(r2)` and `oftype(r1, r2)` will behave differently:
    the first coerces `r2` to be of the type of `r1`, whereas the second converts.
    Developers are urged to future-proof their code by choosing the behavior appropriate for each usage.
"""
struct IdOffsetRange{T<:Integer,I<:AbstractUnitRange{T}} <: AbstractUnitRange{T}
    parent::I
    offset::T

    IdOffsetRange{T,I}(r::I, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}} = new{T,I}(r, offset)

    #= This method is necessary to avoid a StackOverflowError in IdOffsetRange{T,I}(r::IdOffsetRange, offset::Integer).
    The type signature in that method is more specific than IdOffsetRange{T,I}(r::I, offset::T),
    so it ends up calling itself if I <: IdOffsetRange.
    =#
    function IdOffsetRange{T,IdOffsetRange{T,I}}(r::IdOffsetRange{T,I}, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}}
        new{T,IdOffsetRange{T,I}}(r, offset)
    end
end

# Construction/coercion from arbitrary AbstractUnitRanges
function IdOffsetRange{T,I}(r::AbstractUnitRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, o = offset_coerce(I, r)
    return IdOffsetRange{T,I}(rc, convert(T, o+offset)::T)
end
function IdOffsetRange{T}(r::AbstractUnitRange, offset::Integer = 0) where T<:Integer
    rc = convert(AbstractUnitRange{T}, r)::AbstractUnitRange{T}
    return IdOffsetRange{T,typeof(rc)}(rc, convert(T, offset)::T)
end
IdOffsetRange(r::AbstractUnitRange{T}, offset::Integer = 0) where T<:Integer =
    IdOffsetRange{T,typeof(r)}(r, convert(T, offset)::T)

# Coercion from other IdOffsetRanges
IdOffsetRange{T,I}(r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
function IdOffsetRange{T,I}(r::IdOffsetRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, offset_rc = offset_coerce(I, r.parent)
    return IdOffsetRange{T,I}(rc, convert(T, r.offset + offset + offset_rc)::T)
end
function IdOffsetRange{T}(r::IdOffsetRange, offset::Integer = 0) where T<:Integer
    return IdOffsetRange{T}(r.parent, r.offset + offset)
end
IdOffsetRange(r::IdOffsetRange) = r

# TODO: uncomment these when Julia is ready
# # Conversion preserves both the values and the indexes, throwing an InexactError if this
# # is not possible.
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r.parent), r.offset)
# Base.convert(::Type{IdOffsetRange{T,I}}, r::AbstractUnitRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r), 0)

offset_coerce(::Type{Base.OneTo{T}}, r::Base.OneTo) where T<:Integer = convert(Base.OneTo{T}, r), 0
function offset_coerce(::Type{Base.OneTo{T}}, r::AbstractUnitRange) where T<:Integer
    o = first(r) - 1
    return Base.OneTo{T}(last(r) - o), o
end
# function offset_coerce(::Type{Base.OneTo{T}}, r::IdOffsetRange) where T<:Integer
#     rc, o = offset_coerce(Base.OneTo{T}, r.parent)

# Fallback, specialze this method if `convert(I, r)` doesn't do what you need
offset_coerce(::Type{I}, r::AbstractUnitRange) where I<:AbstractUnitRange =
    convert(I, r)::I, 0

@inline Base.parent(r::IdOffsetRange) = r.parent
@inline Base.axes(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.axes1(r::IdOffsetRange) = IdOffsetRange(Base.axes1(r.parent), r.offset)
@inline Base.unsafe_indices(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.length(r::IdOffsetRange) = length(r.parent)
Base.reduced_index(i::IdOffsetRange) = typeof(i)(first(i):first(i))
# Workaround for #92 on Julia < 1.4
Base.reduced_index(i::IdentityUnitRange{<:IdOffsetRange}) = typeof(i)(first(i):first(i))
for f in [:firstindex, :lastindex, :first, :last]
    @eval @inline Base.$f(r::IdOffsetRange) = $f(r.parent) + r.offset
end

@inline function Base.iterate(r::IdOffsetRange)
    ret = iterate(r.parent)
    ret === nothing && return nothing
    return (ret[1] + r.offset, ret[2])
end
@inline function Base.iterate(r::IdOffsetRange, i)
    ret = iterate(r.parent, i)
    ret === nothing && return nothing
    return (ret[1] + r.offset, ret[2])
end

@propagate_inbounds Base.getindex(r::IdOffsetRange, i::Integer) = r.parent[i - r.offset] + r.offset
@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::AbstractUnitRange{<:Integer})
    return r.parent[s .- r.offset] .+ r.offset
end
@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::IdentityUnitRange)
    return IdOffsetRange(r.parent[s .- r.offset], r.offset)
end
@propagate_inbounds function Base.getindex(r::IdOffsetRange, s::IdOffsetRange)
    return IdOffsetRange(r.parent[s.parent .+ (s.offset - r.offset)] .+ (r.offset - s.offset), s.offset)
end

# offset-preserve broadcasting
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(-), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .- x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .+ x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), x::Integer, r::IdOffsetRange{T}) where T =
    IdOffsetRange{T}(x .+ r.parent, r.offset)

Base.show(io::IO, r::IdOffsetRange) = print(io, "OffsetArrays.IdOffsetRange(",first(r), ':', last(r),")")

# Optimizations
@inline Base.checkindex(::Type{Bool}, inds::IdOffsetRange, i::Real) = Base.checkindex(Bool, inds.parent, i - inds.offset)

if VERSION < v"1.5.2"
    # issue 100, 133: IdOffsetRange as another index-preserving case shouldn't comtribute offsets
    # fixed by https://github.com/JuliaLang/julia/pull/37204
    @inline Base.compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{IdOffsetRange}, I::Tuple) =
        Base.compute_linindex(parent, I) - stride1*first(Base.axes1(inds[1]))
end

# This was deemed "too private" to extend: see issue #184
# # Fixes an inference failure in Base.mapfirst!
# # Test: A = OffsetArray(rand(4,4), (-3,5)); R = similar(A, (1:1, 6:9)); maximum!(R, A)
# if isdefined(Base, :_firstslice)
#     Base._firstslice(i::IdOffsetRange) = IdOffsetRange(Base._firstslice(i.parent), i.offset)
# end


text:
abstract type PersistentSet{T} <: AbstractSet{T} end

struct PersistentHashSet{T}
    dict::PersistentHashDict{T, Nothing}
    # TODO: this constructor is inconsistent with everything else
    # and with Set in base; probably good to deprecate.
    #PersistentSet{T}(d::PersistentHashDict{T, Nothing}) where {T} = new{T}(d)
    PersistentSet{T}() where {T} = new{T}(PersistentHashDict{T, Nothing}())
    PersistentSet{T}(els...) where {T} = new{T}(pmap{T,Nothing}([el=>nothing for el in els]...)
end
PersistentSet{T}(itr) where {T} = _union(PersistentSet{T}(), itr)
PersistentSet() = PersistentSet{Any}()
PersistentSet(itr) = PersistentSet{eltype(itr)}(itr)

Base.@deprecate(PersistentSet(x1::T, x2::T, xs::T...) where {T},
                PersistentSet(T[x1, x2, xs...]))

Base.hash(s::PersistentSet,h::UInt) =
    hash(s.dict, h+(UInt(0xf7dca1a5fd7090be)))

Base.conj(s::PersistentSet{T}, val) where {T} =
    PersistentSet{T}(assoc(s.dict, val, nothing))

push(s::PersistentSet, val) = conj(s, val)

disj(s::PersistentSet{T}, val) where {T} =
    PersistentSet{T}(dissoc(s.dict, val))

Base.length(s::PersistentSet) = length(s.dict)
Base.eltype(s::PersistentSet{T}) where {T} = T

Base.isequal(s1::PersistentSet, s2::PersistentSet) = isequal(s1.dict, s2.dict)
==(s1::PersistentSet, s2::PersistentSet) = s1.dict == s2.dict

Base.in(val, s::PersistentSet) = haskey(s.dict, val)

Base.iterate(s::PersistentSet) = iterate(s, iterate(s.dict))
Base.iterate(s::PersistentSet, ::Nothing) = nothing
Base.iterate(s::PersistentSet, (kv, dict_state)) = (kv[1], iterate(s.dict, dict_state))

function Base.filter(f::Function, s::PersistentSet{T}) where T
    filtered = filter((kv) -> f(kv[1]), s.dict)
    PersistentSet{T}(keys(filtered))
end

function Base.setdiff(l::PersistentSet, r::( Union{PersistentSet, Set}))
    notinr(el) = !(el in r)
    filter(notinr, l)
end

import Base.-
-(l::PersistentSet, r::( Union{PersistentSet, Set})) = setdiff(l, r)

Base.isempty(s::PersistentSet) = length(s.dict) == 0

function _union(s::PersistentSet, xs)
    for x in xs
        s = conj(s, x)
    end
    s
end

join_eltype() =  Union{}
join_eltype(v1, vs...) = typejoin(eltype(v1), join_eltype(vs...))

Base.union(s::PersistentSet) = s
Base.union(x::PersistentSet, xs::PersistentSet...) =
    reduce(_union, [x, xs...], init=PersistentSet{join_eltype(x, xs...)}())

function Base.isless(s1::PersistentSet, s2::PersistentSet)
    length(s1) < length(s2) &&
    all(el -> el in s2, s1)
end

import Base.<=
<=(s1::PersistentSet, s2::PersistentSet) = all(el -> el in s2, s1)

function Base.show(io::IO, s::PersistentSet{T}) where T
    print(io, "PersistentSet{$T}(")
    print(io, join([k for (k, v) in s.dict], ", "))
    print(io, ")")
end

ext: jl

max_issues_repo_name: noahbenson/FunctionalUtilities

max_forks_repo_name: noahbenson/FunctionalUtilities


text:
p1 = PetriCospanOb(1)
p2 = PetriCospanOb(2)
pc1 = id(p1)
pc2 = id(p2)
I = munit(PetriCospanOb)

@test dom(pc1 ⊗ pc2) == dom(pc1) ⊗ dom(pc2)
@test codom(pc1 ⊗ pc2) == codom(pc1) ⊗ codom(pc2)
@test dom(σ(p1,p2)) == p1 ⊗ p2
@test codom(σ(p1,p2)) == p2 ⊗ p1
compare_petricospan(σ(p1, p2), braid(p1,p2))

@test dom(Δ(p1)) == p1
@test codom(Δ(p1)) == p1 ⊗ p1
@test dom(◊(p1)) == p1
@test codom(◊(p1)) == I

compare_petricospan(pair(pc1,pc1), Δ(p1) ⋅ (pc1 ⊗ pc1))
compare_petricospan(proj1(p1,p2), pc1 ⊗ ◊(p2))
compare_petricospan(proj2(p1,p2), ◊(p1) ⊗ pc2)


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f5f69190-5f08-11ec-1a42-0952046242a1
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="HypertextLiteral", version="0.6"),
        Pkg.PackageSpec(name="JSON", version="0.21"),
    ])
    using PlutoUI, HypertextLiteral, JSON
end

# ╔═╡ b2bbc3d2-0f4c-4bee-a82e-495c678d2d37
ClickCounter(text="Click") = @htl("""

<button id="add_button">$(text)</button>
<button id="reset_button">reset</button>

""")

# ╔═╡ 2165c042-d7e9-4282-a536-168736e42b5a
@bind num_clicks ClickCounter()

# ╔═╡ 3c6d479f-59f0-4933-a07d-0a866c214fc3
num_clicks

# ╔═╡ 9adf09be-ba1c-4688-9c2f-66f917ecd9b5
x = ["hello pluton!", 232000,2,2,12 ,12,2,21,1,2, 120000]

# ╔═╡ e4379e02-9738-4852-86e3-2467fbe6dce4
state = Dict(
	:x => x
	)

# ╔═╡ 2aca634c-298b-416d-9ab9-835f1a5e022a
@htl("""
<link href="https://cdn.bootcdn.net/ajax/libs/tailwindcss/2.2.7/utilities.min.css" rel="stylesheet">

""")

# ╔═╡ Cell order:
# ╠═f5f69190-5f08-11ec-1a42-0952046242a1
# ╠═b2bbc3d2-0f4c-4bee-a82e-495c678d2d37
# ╠═2165c042-d7e9-4282-a536-168736e42b5a
# ╠═3c6d479f-59f0-4933-a07d-0a866c214fc3
# ╠═2aca634c-298b-416d-9ab9-835f1a5e022a
# ╠═9adf09be-ba1c-4688-9c2f-66f917ecd9b5
# ╠═e4379e02-9738-4852-86e3-2467fbe6dce4

max_stars_repo_name: BioMedAI/Pluto-Note

max_issues_repo_name: BioMedAI/Pluto-Note

max_forks_repo_name: BioMedAI/Pluto-Note


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SeqTests

using Test, OEISUtils, SeqUtils
export SeqTest, is_oeis_installed

const ShowTest = false

# Directory of oeis local data.
srcdir = realpath(joinpath(dirname(@__FILE__)))
ROOTDIR = dirname(srcdir)
datadir = joinpath(ROOTDIR, "data")

"""
Returns the path where the oeis data is expected.
"""
oeis_path() = joinpath(datadir, "stripped")

"""
Indicates if the local copy of the OEIS data (the so-called 'stripped' file) is installed (in ../data).
"""
is_oeis_installed() = isfile(oeis_path())

"""
Indicates if the local copy of the OEIS data (the so-called 'stripped' file) is not installed and warns.
"""
function oeis_notinstalled()
    if !is_oeis_installed()
        @warn("OEIS data not installed! Download stripped.gz from oeis.org,")
        @warn("expand it and put it in the directory ../data.")
        return true
    end
    return false
end

function SeqTest(seqarray, kind, offset=0)
    if kind == 'V' return SeqVTest(seqarray, offset) end
    if kind == 'B' return SeqBTest(seqarray) end
    if kind == 'L' return SeqLTest(seqarray) end
    if kind == 'T' return SeqTTest(seqarray) end
    if kind == 'P' return SeqPTest(seqarray) end

    error("Test function not found!")
end

function SeqVTest(seq, offset=0)
    name = SeqName(seq)
    O = oeis_local(name, 10)
    S = ZArray(10, seq, offset)
    if ShowTest
        println("V --> ", name)
        println(O); println(S)
    end
    @test all(S[1:10] .== O[1:10])
end

function SeqBTest(seqarray)
    for seq in seqarray
        name = SeqName(seq)
        # the parameter is not 'length' but 'search bound'.
        O = oeis_local(name, 12)
        S = seq(300)
        if ShowTest
            println("B --> ", name);
            println(O); println(S)
        end
        @test all(S[0:11] .== O[0:11])
    end
end

function SeqLTest(seqarray)
    for seq in seqarray
        name = SeqName(seq)
        O = oeis_local(name, 12)
        S = seq(12)
        if ShowTest
            println("L --> ", name);
            println(O); println(S)
        end
        @test all(S .== O)
    end
end

function SeqTTest(seqarray)
    for seq in seqarray
        name = SeqName(seq)
        O = oeis_local(name, 21)
        S = seq(6)
        if ShowTest
            println("T --> ", name);
            println(O); println(S)
            # ShowAsΔ(O); ShowAsΔ(S)
        end
        @test all(S .== O)
    end
end

function SeqPTest(seqarray)
    for seq in seqarray
        name = SeqName(seq)
        O = oeis_local(name, 28)
        S = seq(7)
        if ShowTest
            println("P --> ", name);
            ShowAsΔ(O); ShowAsΔ(S)
        end
    end
end

end # module


text:
function Base.parse(::Type{Rational{T}}, x::AbstractString) where {T <: Integer}
    ms, ns = split(x, '/', keepempty = false)
    m = parse(T, ms)
    n = parse(T, ns)
    return m // n
end

Base.parse(::Type{Rational}, x::AbstractString) = parse(Rational{Int32}, x)

function get_ratio(x::AbstractString)
    r = tryparse(Float64, x)
    return r === nothing ? convert(Float64, parse(Rational, x)) : r
end

get_ratio(x::Float64) = x


text:
# Tests for Stationary Arm Model
using InteractiveUtils

@testset "Stationary Arm Models" begin

# Get list of stationary arm models
list_of_arms    = subtypes( Arms.StationaryArmBase )

# For each of them, check whether the basic interface is implemented
# In case of stationary arms, only pull has to be implemented
basefn__pull!   = which( Arms.pull!, Tuple{Arms.StationaryArmBase} )
basefn__tick!   = which( Arms.tick!, Tuple{Arms.StationaryArmBase} )
basefn_reset!   = which( Arms.reset!, Tuple{Arms.StationaryArmBase} )

@testset "Sanity Check" begin
    for arm_model ∈ list_of_arms
        @test basefn__pull! != which( Arms.pull!, Tuple{arm_model} )
        @test basefn__tick! == which( Arms.tick!, Tuple{arm_model} )
        @test basefn_reset! == which( Arms.reset!, Tuple{arm_model} )
    end
end

@testset "Bernoulli Arm" begin
    arm = Arms.Bernoulli( 0.3 )
    @test Arms.mean( arm )      == 0.3
    @test Arms.pull!( arm )     ∈ [0,1]
    @test Arms.tick!( arm )     == nothing
    @test Arms.reset!( arm )    == nothing
    @test Arms.mean( arm )      == 0.3           # Test after reseting
end

@testset "Beta Arm" begin
    arm = Arms.Beta( 2.0, 3.0 )
    @test Arms.mean( arm )      == 0.40
    @test 0.0 <= Arms.pull!( arm ) <= 1.0
    @test Arms.tick!( arm )     == nothing
    @test Arms.reset!( arm )    == nothing
    @test Arms.mean( arm )      == 0.40
end

@testset "Expoenential Arm" begin
    arm = Arms.Exponential( 2.0 )
    @test Arms.mean( arm )      == 2.0
    @test Arms.pull!( arm )     >= 0.0
    @test Arms.tick!( arm )     == nothing
    @test Arms.reset!( arm)     == nothing
    @test Arms.mean( arm )      == 2.0
end

@testset "Normal Arm" begin
    arm = Arms.Normal( 2.0, 1.0 )
    @test Arms.mean( arm )      == 2.0
    @test typeof( Arms.pull!(arm) ) <: Real
    @test Arms.tick!( arm )     == nothing
    @test Arms.reset!( arm )    == nothing
    @test Arms.mean( arm )      == 2.0
end

end


text:
using Primes
import Primes.factor
import Random: rand, SamplerType, AbstractRNG
export factor

"""
    num_irreducibles(::Type{<:UnivariatePolynomial{F}}, r)

Number of irreducible polynomials over `F` of degree `r`.
"""
function num_irreducibles(::Type{<:UnivariatePolynomial{G}}, r::Integer) where G
    k = order(G)
    iszero(k) && throw(ArgumentError("order of base type is zero"))
    T = mintype_for(k, r, false)
    necklace(T(k), r)
end
function num_irreducibles(::Type{G}, r) where G<:Ring
    num_irreducibles(G[:x], r)
end

"""
    isirreducible(p::F[X])

Returns iff `p` is an irreducible (prime) polynomial over field `F`. See also `factor`.
"""
function isirreducible(p::UnivariatePolynomial{<:QuotientRing})
    deg(p) <= 1 && return true
    iszero(p.coeff[1]) && return false
    pp = gcd(p, derive(p))
    deg(pp) > 0 && return false
    isddf(p)
end
"""
    isreducible(p::F[X])

Returns iff `p` is a reducible (prime) polynomial over field `F`. See also `factor`.
"""
isreducible(p::UnivariatePolynomial{<:QuotientRing}) = !isirreducible(p)

import Base.Iterators: Filter, take, drop
"""::UnivariatePolynomial{<:QuotientRing}
    irreducible(P, n, nr=0)

Returns an irreducible polynomial with in `P` with degree `n`. Skip first `nr` hits.
"""
irreducible(::Type{P}, n) where P<:UnivariatePolynomial = first(irreducibles(P, n))
function irreducible(::Type{P}, n, nr::Integer) where P<:UnivariatePolynomial
    first(drop(irreducibles(P, n), nr))
end
"""
    reducible(P, n, nr=0)

Returns a reducible polynomial with in `P` with degree `n`. Skip first `nr` hits.
"""
reducible(::Type{P}, n) where P<:UnivariatePolynomial = first(reducibles(P, n))
function reducible(::Type{P}, n, nr::Integer) where P<:UnivariatePolynomial
    first(drop(reducibles(P, n), nr))
end
"""
    irreducibles(P, n)

Returns iterator of all irreducible monic polynomials in `P` with degree `n`.
"""
function irreducibles(::Type{P}, n) where P<:UnivariatePolynomial{<:QuotientRing}
    Base.Iterators.Filter(isirreducible, Monic(P, n))
end
function irreducibles(::Type{G}, n) where G<:Ring
    irreducibles(G[:x], n)
end
"""
reducibles(P, n)

Returns iterator of all reducible monic polynomials in `P` with degree `n`.
"""
function reducibles(::Type{P}, n) where P<:UnivariatePolynomial{<:QuotientRing}
    Base.Iterators.Filter(!isirreducible, Monic(P, n))
end
function reducibles(::Type{G}, n) where G<:Ring
    reducibles(G[:x], n)
end

"""
    factor(p::F[:x])

Factorize polynomial in `F[X]` where `F` is a field
(`ZZ/p` or `GF(p,m)` with `p` prime number).
"""
function factor(p::P) where P<:UnivariatePolynomial{<:QuotientRing}
    res = Pair{P,Int}[]
    u = lcunit(p)
    if !isone(u)
        p /= u
        push!(res, P(u) => 1)
    end
    if deg(p) <= 1
        if !isone(p) || isempty(res)
            push!(res, p => 1)
        end
        return res
    end
    pp = sff(p)
    for (q, k) in pp
        qq = ddf(q)
        for (r, l) in qq
            rr = edf(r, l)
            for s in rr
                push!(res, s => k)
            end
        end
    end
    sort!(res)
end

"""
    sff(p)

`Square-free factorization`.
Algorithm to split polynomial `p` into a product of powers of squarefree factors.
Return an array of pairs of squarefree factors and corresponding powers.
"""
function sff(f::P) where {Z<:QuotientRing,P<:UnivariatePolynomial{Z}}
    q = order(Z)
    p = characteristic(Z)
    i = 1
    R = Pair{P,Int}[]
    fs = derive(f)
    c = gcd(f, fs) # c contains all multiple factors of f
    w = f / c # w is square-free
    while !isunit(w)
        y = gcd(w, c)
        z = w / y
        if deg(z) > 0
            push!(R, Pair(z, i))
        end
        c /= y
        w = y
        i += 1
    end

    if deg(c) > 0
        c = shrink(c, p)
        for (g, i) in sff(c)
            push!(R, Pair(g, i * p))
        end
    end
    R
end

# assuming p(x) = q(x^p), return q. Formally q(x) = p(x^(1/p)).
function shrink(a::P, p::Integer) where P<:UnivariatePolynomial
    n = deg(a)
    c = similar(a.coeff, (n÷p)+1)
    for k = 0:n÷p
        c[k+1] = a.coeff[p*k+1]
    end
    P(c)
end

"""
    ddf(p)

`Distinct-degree factorization`.
Input is a squarefree polynomial.
Returns a list of pairs `g_i => d_i` of polynomials g_i, each of which is a product of
all irreducible monic polynomials of equal degree `d_i`. The product of all `g_i == p`.
"""
function ddf(f::P) where {Z<:QuotientRing, P<:UnivariatePolynomial{Z}}
    q = order(Z)
    S = Pair{P, Int}[]
    x = monom(typeof(f), 1)
    i = 1
    fs = f
    xqi = x
    while deg(fs) >= 2i
        xqi = powermod(xqi, q, fs)
        g = gcd(fs, xqi - x)
        if deg(g) > 0
            push!(S, Pair(g, i))
            fs /= g
        end
        i += 1
    end
    if deg(fs) > 0
        push!(S, Pair(fs, deg(fs)))
    end
    if isempty(S)
        push!(S, Pair(f, 1))
    end
    S
end

function isddf(f::P) where {Z<:QuotientRing, P<:UnivariatePolynomial{Z}}
    q = order(Z)
    x = monom(typeof(f), 1)
    i = 1
    fs = f
    xqi = x
    while deg(fs) >= 2i
        xqi = powermod(xqi, q, fs)
        g = gcd(fs, xqi - x)
        deg(g) > 0 && return false
        i += 1
    end
    return true
end

"""
    edf(p::Polynomial, d::Integer)

`Equal-degree factorization`.
Algorithm of Cantor-Zassenhaus to find the factors of `p`, a product of monomials of
degree `d`. (Such polynomials are in the output of `ddf`).
The base type for `p` must be a finite field. Odd charcteristic is a covered special case.
"""
function edf(f::P, d::Integer) where {Z<:QuotientRing,P<:UnivariatePolynomial{Z}}
    q = order(Z)
    n = deg(f)
    S = [f]
    n == d && return S
    rem(n, d) == 0 || throw(DomainError((n, d), "degree of f must be multiple of d = $d"))
    ex = big(q)^d ÷ 2 # isodd(q) ? (q^d - 1) ÷ 2 : (q^d ÷ 2)
    r = div(n, d)
    power = isodd(q) ? powermod : powersum
    while length(S) < r
        h = P(rand(Z, n))
        s = length(S)
        g = power(h, ex, f) - 1
        for k = 1:s
            u = S[k]
            gu = gcd(g, u)
            if 0 < deg(gu) < deg(u)
                S[k] = gu
                push!(S, u / gu)
            end
        end
    end
    S
end

"""
    powersum(h, ex, f)

Calculate the sum `h + h^2 + h^4 + h^8 + ... + h^ex mod f`
"""
function powersum(h, ex, f)
    s = h
    n = 1
    while n < ex
        h = rem(h * h, f)
        n *= 2
        s += h
    end
    s
end

# random samplers for finite rings
function rand(r::AbstractRNG, ::SamplerType{Z}) where {Z<:ZZmod}
    m = modulus(Z)
    Z(rand(r, 0:m-1))
end
# Random field element of `Q = P / (polynomial)`, whith `basetype(P) <: ZZmod`.
function rand(r::AbstractRNG, ::SamplerType{Q}) where {Z,P<:UnivariatePolynomial{Z},Q<:Quotient{P}}

    m = deg(modulus(Q))
    r = Q(P(rand(r, Z, m)))
end

function Base.isless(p::T, q::T) where T<:Pair{<:Ring,<:Integer}
    p.first < q.first ||
    p.first == q.first && p.second == q.second
end

import Base: prod
function Base.prod(ff::Vector{<:Pair{T,<:Integer}}) where T<:Ring
    res = one(T)
    for p in ff
        res *= p.first^p.second
    end
    res
end


text:
push!(LOAD_PATH, "..")

using Documenter
using Bibliography
using Literate
using Plots  # to avoid capturing precompilation output by Literate
using Oceananigans
using Oceananigans.Operators
using Oceananigans.Grids
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.TurbulenceClosures
using Oceananigans.TimeSteppers
using Oceananigans.AbstractOperations

#=
bib_filepath = joinpath(dirname(@__FILE__), "oceananigans.bib")
const BIBLIOGRAPHY = import_bibtex(bib_filepath)
@info "Bibliography: found $(length(BIBLIOGRAPHY)) entries."

include("bibliography.jl")
include("citations.jl")
=#

#####
##### Generate examples
#####

# Gotta set this environment variable when using the GR run-time on Travis CI.
# This happens as examples will use Plots.jl to make plots and movies.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = "100"

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")

examples = [
           # "internal_wave.jl",
            "eady_turbulence.jl"
           ]

for example in examples
    example_filepath = joinpath(EXAMPLES_DIR, example)
    Literate.markdown(example_filepath, OUTPUT_DIR, documenter=true)
end

#####
##### Organize page hierarchies
#####

example_pages = [
                 #"Internal wave"                    => "generated/internal_wave.md",
                 "Eady turbulence"                  => "generated/eady_turbulence.md"
                ]

pages = [
         "Home" => "index.md",
         "Examples" => example_pages
        ]

#####
##### Build and deploy docs
#####

format = Documenter.HTML(collapselevel = 1,
                            prettyurls = false
                        )

makedocs(sitename = "Oceananigans.jl",
          authors = "Ali Ramadhan, Gregory Wagner, John Marshall, Jean-Michel Campin, Chris Hill",
           format = format,
            pages = pages,
          modules = [Oceananigans],
          doctest = false,
           strict = false,
            clean = true,
        checkdocs = :none  # Should fix our docstring so we can use checkdocs=:exports with strict=true.
        )


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Qt6Base"
version = v"6.3.0"

# Set this to true first when updating the version. It will build only for the host (linux musl).
# After that JLL is in the registyry, set this to false to build for the other platforms, using
# this same package as host build dependency.
const host_build = false

# Collection of sources required to build qt6
sources = [
    ArchiveSource("https://download.qt.io/official_releases/qt/$(version.major).$(version.minor)/$version/submodules/qtbase-everywhere-src-$version.tar.xz",
                  "b865aae43357f792b3b0a162899d9bf6a1393a55c4e5e4ede5316b157b1a0f99"),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/11.0-11.1/MacOSX11.1.sdk.tar.xz",
                  "9b86eab03176c56bb526de30daa50fa819937c54b280364784ce431885341bf6"),
    ArchiveSource("https://sourceforge.net/projects/mingw-w64/files/mingw-w64/mingw-w64-release/mingw-w64-v10.0.0.tar.bz2",
                  "ba6b430aed72c63a3768531f6a3ffc2b0fde2c57a3b251450dcf489a894f0894")
]

script = raw"""
cd $WORKSPACE/srcdir

BIN_DIR="/opt/bin/${bb_full_target}"

mkdir build
cd build/

qtsrcdir=`ls -d ../qtbase-everywhere-src-*`

commonoptions=" \
-opensource -confirm-license \
-openssl-linked  -nomake examples -release \
"

commoncmakeoptions="-DCMAKE_PREFIX_PATH=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DQT_HOST_PATH=$host_prefix -DQT_FEATURE_openssl_linked=ON"

export LD_LIBRARY_PATH=$host_libdir:$LD_LIBRARY_PATH
export OPENSSL_LIBS="-L${libdir} -lssl -lcrypto"

sed -i 's/"-march=haswell"/"-mavx2" "-mf16c"/' $qtsrcdir/cmake/QtCompilerOptimization.cmake

case "$target" in

    x86_64-linux-musl*)
        export LD_LIBRARY_PATH=$WORKSPACE/srcdir/build/lib:$host_libdir:$LD_LIBRARY_PATH
        ../qtbase-everywhere-src-*/configure -prefix $prefix $commonoptions -fontconfig -- -DCMAKE_PREFIX_PATH=${prefix} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_HOST_TOOLCHAIN}
    ;;

    *mingw*)
        cd $WORKSPACE/srcdir/mingw*/mingw-w64-headers
        ./configure --prefix=/opt/$target/$target/sys-root --enable-sdk=all --host=$target
        make install

        cd ../mingw-w64-crt/
        if [ ${target} == "i686-w64-mingw32" ]; then
            _crt_configure_args="--disable-lib64 --enable-lib32"
        elif [ ${target} == "x86_64-w64-mingw32" ]; then
            _crt_configure_args="--disable-lib32 --enable-lib64"
        fi
        ./configure --prefix=/opt/$target/$target/sys-root --enable-sdk=all --host=$target --enable-wildcard ${_crt_configure_args}
        make -j${nproc}
        make install

        cd ../mingw-w64-libraries/winpthreads
        ./configure --prefix=/opt/$target/$target/sys-root --host=$target --enable-static --enable-shared
        make -j${nproc}
        make install

        cd $WORKSPACE/srcdir/build
        ../qtbase-everywhere-src-*/configure -prefix $prefix $commonoptions -opengl dynamic -- $commoncmakeoptions
    ;;

    *apple-darwin*)
        apple_sdk_root=$WORKSPACE/srcdir/MacOSX11.1.sdk
        sed -i "s!/opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root!$apple_sdk_root!" $CMAKE_TARGET_TOOLCHAIN
        ../qtbase-everywhere-src-*/configure -prefix $prefix $commonoptions -- $commoncmakeoptions -DCMAKE_SYSROOT=$apple_sdk_root -DCMAKE_FRAMEWORK_PATH=$apple_sdk_root/System/Library/Frameworks -DCMAKE_OSX_DEPLOYMENT_TARGET=10.15 -DCUPS_INCLUDE_DIR=$apple_sdk_root/usr/include -DCUPS_LIBRARIES=$apple_sdk_root/usr/lib/libcups.tbd
    ;;

    *freebsd*)
        sed -i 's/-Wl,--no-undefined//' $qtsrcdir/mkspecs/freebsd-clang/qmake.conf
        sed -i 's/-Wl,--no-undefined//' $qtsrcdir/cmake/QtFlagHandlingHelpers.cmake
        ../qtbase-everywhere-src-*/configure -prefix $prefix $commonoptions -fontconfig -- $commoncmakeoptions -DQT_PLATFORM_DEFINITION_DIR=$host_prefix/mkspecs/freebsd-clang
    ;;

    *)
        if [[ "$target" != *"aarch64"* ]]; then
            sed -i 's/ELFOSABI_GNU/ELFOSABI_LINUX/' $qtsrcdir/src/corelib/plugin/qelfparser_p.cpp
            sed -i 's/.*EM_AARCH64.*//' $qtsrcdir/src/corelib/plugin/qelfparser_p.cpp
        fi
        sed -i 's/.*EM_BLACKFIN.*//' $qtsrcdir/src/corelib/plugin/qelfparser_p.cpp
        ../qtbase-everywhere-src-*/configure -prefix $prefix $commonoptions -fontconfig -- $commoncmakeoptions
    ;;
esac

cmake --build . --parallel ${nproc}
cmake --install .
install_license $WORKSPACE/srcdir/qtbase-everywhere-src-*/LICENSE.LGPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
if host_build
    platforms = [Platform("x86_64", "linux",cxxstring_abi=:cxx11,libc="musl")]
else
    platforms = expand_cxxstring_abis(filter(!Sys.isapple, supported_platforms()))
    filter!(p -> arch(p) != "armv6l", platforms) # No OpenGL on armv6
    platforms_macos = [ Platform("x86_64", "macos"), Platform("aarch64", "macos") ]
end

# The products that we will ensure are always built
products = [
    LibraryProduct(["Qt6Concurrent", "libQt6Concurrent", "QtConcurrent"], :libqt6concurrent),
    LibraryProduct(["Qt6Core", "libQt6Core", "QtCore"], :libqt6core),
    LibraryProduct(["Qt6DBus", "libQt6DBus", "QtDBus"], :libqt6dbus),
    LibraryProduct(["Qt6Gui", "libQt6Gui", "QtGui"], :libqt6gui),
    LibraryProduct(["Qt6Network", "libQt6Network", "QtNetwork"], :libqt6network),
    LibraryProduct(["Qt6OpenGL", "libQt6OpenGL", "QtOpenGL"], :libqt6opengl),
    LibraryProduct(["Qt6PrintSupport", "libQt6PrintSupport", "QtPrintSupport"], :libqt6printsupport),
    LibraryProduct(["Qt6Sql", "libQt6Sql", "QtSql"], :libqt6sql),
    LibraryProduct(["Qt6Test", "libQt6Test", "QtTest"], :libqt6test),
    LibraryProduct(["Qt6Widgets", "libQt6Widgets", "QtWidgets"], :libqt6widgets),
    LibraryProduct(["Qt6Xml", "libQt6Xml", "QtXml"], :libqt6xml),
]

products_macos = [
    FrameworkProduct("QtConcurrent", :libqt6concurrent),
    FrameworkProduct("QtCore", :libqt6core),
    FrameworkProduct("QtDBus", :libqt6dbus),
    FrameworkProduct("QtGui", :libqt6gui),
    FrameworkProduct("QtNetwork", :libqt6network),
    FrameworkProduct("QtOpenGL", :libqt6opengl),
    FrameworkProduct("QtPrintSupport", :libqt6printsupport),
    FrameworkProduct("QtSql", :libqt6sql),
    FrameworkProduct("QtTest", :libqt6test),
    FrameworkProduct("QtWidgets", :libqt6widgets),
    FrameworkProduct("QtXml", :libqt6xml),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Xorg_libX11_jll"),
    Dependency("Xorg_libXext_jll"),
    BuildDependency("Xorg_xproto_jll"),
    BuildDependency("Xorg_glproto_jll"),
    Dependency("Xorg_libxcb_jll"),
    Dependency("Xorg_xcb_util_wm_jll"),
    Dependency("Xorg_xcb_util_image_jll"),
    Dependency("Xorg_xcb_util_keysyms_jll"),
    Dependency("Xorg_xcb_util_renderutil_jll"),
    BuildDependency("Xorg_kbproto_jll"),
    BuildDependency("Xorg_renderproto_jll"),
    Dependency("Xorg_libXrender_jll"),
    Dependency("xkbcommon_jll"),
    Dependency("Libglvnd_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("Glib_jll", v"2.59.0"; compat="2.59.0"),
    Dependency("Zlib_jll"),
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenSSL_jll"),
    BuildDependency(PackageSpec(name="LLVM_full_jll", version=v"11.0.1")),
]

if !host_build
    push!(dependencies, HostBuildDependency("Qt6Base_jll"))
end

include("../../fancy_toys.jl")

@static if !host_build
    if any(should_build_platform.(triplet.(platforms_macos)))
        build_tarballs(ARGS, name, version, sources, script, platforms_macos, products_macos, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
    end
end
if any(should_build_platform.(triplet.(platforms)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"9", julia_compat="1.6")
end


text:
## Initialization Filters
using LinearAlgebra, FFTW

"""
dims is the size of the filters
e.g., dims=(3,2) will give you 2 dimensional filters, each of size 3x2

These routines always return prod(dims) number of filters

:DCT filters are ordered from low to high frequency
"""
generatefilters(type,dims;form=:list) = _generatefilters(Val(type),dims,Val(form))

# DCT
_dctmatrix(n) = dct(Matrix(I,(n,n)),1)'
function _generatefilters(::Val{:DCT},dims,::Val{:list})
    H = _generatefilters(Val(:DCT),dims,Val(:matrix))
    return Array.(_filterlist(H,dims))
end
function _generatefilters(::Val{:DCT},dims,::Val{:matrix})
    @assert length(dims) >= 1  "Must provide at least one dimention for filter initialization."
    @assert all(dims .>= 1) "All filter dimensions must have positive size"
    H = _dctmatrix(dims[1])
    for i in 2:length(dims)
        H = kron(_dctmatrix(dims[i]), H)
    end
    @assert H'H ≈ I "DCT filter generation led to non-orthogonal filters"
    h = _reordersnake(Array.(_filterlist(H,dims)), dims)
    return  _filtermatrix(h)[1] ./ sqrt(prod(dims))
end

# Helper functions
function _reordersnake(h, dims)
    freqmeasure = zeros(dims)
    for i in 1:length(dims)
        orderdims = [i; 1:(i-1); (i+1):length(dims)]
        reorderdims = [2:i; 1; (i+1):length(dims)]
        repdims = [1; collect(dims[orderdims][2:end])]
        temp = repeat(1:dims[i], outer=repdims)
        freqmeasure += permutedims(temp, reorderdims)
    end
    perminds = sortperm(vec(freqmeasure))
    return h[perminds]
end


text:
using QuantumESPRESSOBase
using Documenter

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

makedocs(;
    modules=[QuantumESPRESSOBase],
    authors="Qi Zhang <singularitti@outlook.com>",
    repo="https://github.com/MineralsCloud/QuantumESPRESSOBase.jl/blob/{commit}{path}#{line}",
    sitename="QuantumESPRESSOBase.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://MineralsCloud.github.io/QuantumESPRESSOBase.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => ["Installation" => "install.md", "Development" => "develop.md"],
        "API by module" => [
            "`QuantumESPRESSOBase` module" => "api/api.md",
            "`Inputs` module" => "api/Inputs/Inputs.md",
            "`Inputs.PWscf` module" => "api/Inputs/PWscf.md",
            "`Inputs.PHonon` module" => "api/Inputs/PHonon.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/MineralsCloud/QuantumESPRESSOBase.jl",
)


text:
module Mustache

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

using Tables

include("utils.jl")
include("tokens.jl")
include("context.jl")
include("writer.jl")
include("parse.jl")
include("render.jl")

export @mt_str, @jmt_str, render, render_from_file

end


text:
"""
Calculate global coordinates of solution points

    global_sp(xi, r)

Line elements
"""
function global_sp(xi::AbstractVector{T}, r::AbstractVector{T}) where {T<:Real} # line elements
    xsp = similar(xi, first(eachindex(xi)):last(eachindex(xi))-1, length(r))
    for i in axes(xsp, 1), j in axes(r, 1)
        xsp[i, j] = r_x(r[j], xi[i], xi[i+1])
    end

    return xsp
end

"""
    global_sp(xi, yi, r)

Quadrilateral elements
"""
function global_sp(
    xi::AbstractArray{<:Real,2},
    yi::AbstractArray{<:Real,2},
    r::AbstractArray{<:Real,1},
)
    xp = similar(
        xi,
        first(axes(xi, 1)):last(axes(xi, 1))-1,
        first(axes(xi, 2)):last(axes(xi, 2)),
        length(r),
        length(r),
    )
    yp = similar(xp)
    for i in axes(xp, 1), j in axes(xp, 2), k in axes(r, 1), l in axes(r, 1)
        xp[i, j, k, l] = r_x(r[k], xi[i, j], xi[i+1, j])
        yp[i, j, k, l] = r_x(r[l], yi[i, j], yi[i, j+1])
    end

    return xp, yp
end

"""
    global_sp(points, cellid, N)

Triangle elements

- @arg points: vetices of elements
- @arg cellid: point ids of elements
- @arg N: polynomial degree

The right triangle is used in the reference space
- 1, 2: bottom points
- 3: top point
"""
function global_sp(
    points::AbstractMatrix{T1},
    cellid::AbstractMatrix{T2},
    N::Integer,
) where {T1<:Real,T2<:Integer}
    pl, wl = tri_quadrature(N)
    Np = size(wl, 1)

    spg = zeros(eltype(points), size(cellid, 1), Np, 2)
    for i in axes(spg, 1), j in axes(spg, 2)
        id1, id2, id3 = cellid[i, :]
        spg[i, j, :] .=
            rs_xy(pl[j, :], points[id1, 1:2], points[id2, 1:2], points[id3, 1:2])
    end

    return spg
end

#--- deprecated codes for equilateral triangle element ---#
#=function global_sp(
    points::AbstractMatrix{T1},
    cellid::AbstractMatrix{T2},
    rs::AbstractMatrix{T3},
) where {T1<:Real,T2<:Integer,T3<:Real}
    r = rs[:, 1]
    s = rs[:, 2]
    xsp = similar(points, size(cellid, 1), size(rs, 1), size(rs, 2))

    for i in axes(xsp, 1), j in axes(xsp, 2)
        xsp[i, j, :] =
            (-3.0 * r[j] + 2.0 - √3 * s[j]) / 6.0 .* points[cellid[i, 1], 1:2] +
            (3.0 * r[j] + 2.0 - √3 * s[j]) / 6.0 .* points[cellid[i, 2], 1:2] +
            (2.0 + 2.0 * √3 * s[j]) / 6.0 .* points[cellid[i, 3], 1:2]
    end

    return xsp
end=#

"""
Calculate global coordinates of flux points

    global_fp(points, cellid, N)

Triangle elements
"""
function global_fp(points, cellid, N)
    pf, wf = triface_quadrature(N)

    fpg = zeros(size(cellid, 1), 3, N + 1, 2)
    for i in axes(fpg, 1)
        id1, id2, id3 = cellid[i, :]
        for j in axes(fpg, 2), k in axes(fpg, 3)
            fpg[i, j, k, :] .=
                rs_xy(pf[j, k, :], points[id1, 1:2], points[id2, 1:2], points[id3, 1:2])
        end
    end

    return fpg
end

#--- deprecated face-indexed codes ---#
#=function global_fp(points, cellid, faceCells, facePoints, N)
    pf, wf = triface_quadrature(N)

    fpg = zeros(size(faceCells, 1), N+1, 2)
    for i in axes(fpg, 1), j in axes(fpg, 2)
        idc = ifelse(faceCells[i, 1] != -1, faceCells[i, 1], faceCells[i, 2])
        id1, id2, id3 = cellid[idc, :]

        if !(id3 in facePoints[i, :])
            idf = 1
        elseif !(id1 in facePoints[i, :])
            idf = 2
        elseif !(id2 in facePoints[i, :])
            idf = 3
        end

        fpg[i, j, :] .= rs_xy(pf[idf, j, :], points[id1, 1:2], points[id2, 1:2], points[id3, 1:2])
    end

    return fpg
end=#


text:
using Statistics
using PyCall
pygui(:tk)
import PyPlot; const plt = PyPlot
import InteratomicPotentials as Potentials

n=50
########################################################################################
########################################################################################
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 2, 1)
rdf_true = zeros(100, n)
rdf_snap = zeros(100, n)
rdf_lj = zeros(100, n)
for j = 1:n
    # println(j)
    # println("True")
    bins, rdf = Potentials.get_rdf("examples/Argon13/DATA_training/true_ensemble/seed_$j/tmp.rdf")
    rdf_true[:, j] = rdf
    if j == 1
        ax.plot(bins, rdf, color="red", linewidth = 1.0, linestyle="--", label= "True LJ")
    else
        ax.plot(bins, rdf, color="red", linewidth = 0.1, linestyle="--")
    end

    # println("LJ")
    bins, rdf = Potentials.get_rdf("examples/Argon13/fit_lj/samples/sample_1/seed_$j/tmp.rdf")
    rdf_lj[:, j] = rdf
    if j == 1
        ax.plot(bins, rdf, color="blue", linewidth = 1.0, linestyle="--", label= "Fitted LJ")
    else
        ax.plot(bins, rdf, color="blue", linewidth = 0.1, linestyle="--")
    end

    # println("SNAP")
    bins, rdf = Potentials.get_rdf("examples/Argon13/fit_snap/samples/sample_1/seed_$j/tmp.rdf")
    if j == 1
        ax.plot(bins, rdf, color="black", linewidth = 1.0, linestyle="--", label= "Fitted SNAP")
    else
        ax.plot(bins, rdf, color="black", linewidth = 0.1, linestyle="--")
    end
    rdf_snap[:, j] = rdf

end
max_var_true = round(maximum(std(rdf_true, dims = 2)), digits = 2)
max_var_lj = round(maximum(std(rdf_lj, dims = 2)), digits = 2)
max_var_snap = round(maximum(std(rdf_snap, dims = 2)), digits = 2)
ax.text(1.5, 20.0, "Maximum STD: \n True LJ:         $max_var_true \n Fitted LJ:       $max_var_lj \n Fitted SNAP: $max_var_snap")
ax.set_title("Variability of RDF over random seeds")
ax.set_ylabel("RDF")
ax.set_xlabel("r (LJ Units)")
ax.legend()
ax.set_xlim([0.75, 2.5])
plt.grid("on")

# ax = fig.add_subplot(1, 3, 2)
# for j = 1:n
#     bins, rdf = Potentials.get_rdf("examples/Argon13/DATA_training/true_ensemble/seed_$j/tmp.rdf")
#     ax.plot(bins, rdf, color="red", linewidth = 0.1, linestyle="--")
# end

# rdf_snap = zeros(100, 48)
# rdf_lj = zeros(100, 48)
# count = 1
# for j = 1:4
#     for i = 1:12
#     bins, rdf = Potentials.get_rdf("examples/Argon13/fit_snap/samples/sample_$i/seed_$j/tmp.rdf")
#     rdf_snap[:, count] = rdf
#     ax.plot(bins, rdf, color = "black", linewidth = 0.1, linestyle = "--")

#     bins, rdf = Potentials.get_rdf("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.rdf")
#     rdf_lj[:, count] = rdf
#     ax.plot(bins, rdf, color="blue", linewidth = 0.1, linestyle="--")
#     plt.title("Variability of RDF over random samples")
#     count += 1
#     end
# end
# max_var_lj = round(maximum(std(rdf_lj, dims = 2)), digits = 2)
# max_var_snap = round(maximum(std(rdf_snap, dims = 2)), digits = 2)
# ax.text(1.5, 20.0, "Maximum STD: \n \n Fitted LJ:       $max_var_lj \n Fitted SNAP: $max_var_snap")
# ax.legend()
# ax.set_xlabel("r (LJ Units)")
# ax.set_xlim([0.75, 2.5])
# plt.grid("on")

ax = fig.add_subplot(1, 2,2)
for j = 1:n
    bins, rdf = Potentials.get_rdf("examples/Argon13/DATA_training/true_ensemble/seed_$j/tmp.rdf")
    ax.plot(bins, rdf, color="red", linewidth = 0.1, linestyle="--")
end

rdf_snap = zeros(100, 40*12)
rdf_lj = zeros(100, 40*12)
count = 1
for j = 1:50
    # println(j)
    for i = 1:12
        # println(i)
        try Potentials.get_rdf("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.rdf")
            # println("SNAP")
            bins, rdf = Potentials.get_rdf("examples/Argon13/fit_snap/samples/sample_$i/seed_$j/tmp.rdf")
            rdf_snap[:, count] = rdf
            ax.plot(bins, rdf, color = "black", linewidth = 0.1, linestyle = "--")

            # println("LJ")
            bins, rdf = Potentials.get_rdf("examples/Argon13/fit_lj/samples/sample_$i/seed_$j/tmp.rdf")
            rdf_lj[:, count] = rdf
            ax.plot(bins, rdf, color="blue", linewidth = 0.1, linestyle="--")
            plt.title("Variability of RDF over seeds and samples")
            count += 1
        catch
            continue
        end

        if count > 40*12
            break
        end
    end

end
max_var_lj = round(maximum(std(rdf_lj, dims = 2)), digits = 2)
max_var_snap = round(maximum(std(rdf_snap, dims = 2)), digits = 2)
ax.text(1.5, 20.0, "Maximum STD: \n \n Fitted LJ:       $max_var_lj \n Fitted SNAP: $max_var_snap")
ax.legend()
ax.set_xlabel("r (LJ Units)")
ax.set_xlim([0.75, 2.5])
plt.grid("on")

plt.gcf()
plt.savefig("variation_of_rdf_two.png", dpi=600, bbox_inches = "tight")


text:
export Disk,𝔻

# T is Real or Complex
# TT is (Real,Real) or Complex
"""
    Disk(c,r)

represents the disk centred at `c` with radius `r`.
"""
struct Disk{T} <: Domain{T}
    center::T
    radius::Float64
end

Disk(c,r) = Disk{typeof(c)}(c,r)
Disk(c::Real,r) = Disk(complex(c),r)
Disk(c::Integer,r) = Disk(float(c),r)
Disk(c::Complex{<:Integer},r) = Disk(float(c),r)

Disk() = Disk(Vec(0.,0.),1.)
Disk(::AnyDomain) = Disk(NaN,(NaN,NaN))

const 𝔻 = Disk()

isambiguous(d::Disk) = isnan(d.radius) && all(isnan,d.center)

polar(r,θ) = Vec(r*cos(θ),r*sin(θ))
ipolar(x,y) = Vec(sqrt(abs2(x)+abs2(y)), atan(y,x))
polar(rθ::Vec) = polar(rθ...)
ipolar(xy::Vec) = ipolar(xy...)

fromcanonical(D::Disk{T},xy::Vec) where {T<:Vec} = Vec((D.radius.*xy .+ D.center)...)
tocanonical(D::Disk{T},x,y) where {T<:Vec} = Vec((xy .- D.center)./D.radius)

checkpoints(d::Disk) = [fromcanonical(d,Vec(.1,.2243)),fromcanonical(d,Vec(-.212423,-.3))]

# function points(d::Disk,n,m,k)
#     ptsx=0.5*(1-gaussjacobi(n,1.,0.)[1])
#     ptst=points(PeriodicSegment(),m)
#
#     Float64[fromcanonical(d,x,t)[k] for x in ptsx, t in ptst]
# end

boundary(d::Disk) = Circle(d.center,d.radius)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.DLR_Visualization (Modia3D/renderer/DLR_Visualization/_module.jl)
#

function Composition.initializeVisualization(renderer::DummyRenderer, velements::Vector{Composition.Object3D})::Nothing
end

function Composition.visualize!(renderer::DummyRenderer, time::Float64)
end

function Composition.closeVisualization(renderer::DummyRenderer)
end

Composition.isVisible(feature::Modia3D.AbstractVisualElement, renderer::DummyRenderer) = false


text:
const KFOLDCROSSVALIDATION_DEFAULT = false
const NUMFOLDS_DEFAULT = 5
const TESTSETSHARE_DEFAULT = 0.15
const VALID_PRELIMINARYSELECTION = [:lasso]
const EXPORTCSV_DEFAULT = nothing
const EXPORTSUMMARY_DEFAULT = nothing
const EXPORTLATEX_DEFAULT = nothing


text:
"""
    RotationUtilModule

Module to handle two-dimensional and three-dimensional rotations.
"""
module RotationUtilModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import LinearAlgebra: norm

_I3 = [i==j ? one(FFlt) : zero(FFlt) for i=1:3, j=1:3]

"""
    rotmat3!(Rmout::FFltMat, a::T) where {T}

Compute a 3D rotation matrix in-place.

`a` = array, vector, or tuple with three floating-point numbers
"""
function rotmat3!(Rmout::FFltMat, a::T) where {T}
    m,n = size(Rmout);
    @assert (m == n) && (m == 3) && (length(a) == m)
    na = norm(a);
    if (na == 0.)
        for i = 1:size(Rmout, 1), j = 1:size(Rmout, 2)
            Rmout[i, j] = (i == j ? 1.0 : 0.0)
        end
    else
        # thetatilde = zeros(3,3);
        # skewmat!(thetatilde, a);
        # a = a/na;
        # ca = cos(na);
        # sa = sin(na);
        # aa = a*a';
        # copyto!(Rmout, ca * (_I3-aa) + sa/na*thetatilde + aa);
        as = (a[1]/na, a[2]/na, a[3]/na)
        ca = cos(na);
        sa = sin(na);
        oca = (1.0 - ca)
        for j in 1:3, i in 1:3
            Rmout[i, j] = oca * as[i] * as[j]
        end
        Rmout[1, 2] += sa * -as[3]
        Rmout[1, 3] += sa * as[2];
        Rmout[2, 1] += sa * as[3]
        Rmout[2, 3] += sa * -as[1];
        Rmout[3, 1] += sa * -as[2];
        Rmout[3, 2] += sa * as[1];
        Rmout[1, 1] += ca
        Rmout[2, 2] += ca
        Rmout[3, 3] += ca
    end
    return Rmout
end

"""
    rotmat3(a::T) where {T}

Prepare a rotation matrix from a rotation vector
"""
function rotmat3(a::T) where {T}
    Rmout = fill(0.0, 3, 3)
    return rotmat3!(Rmout, a)
end

"""
    skewmat!(S, theta)

Compute skew-symmetric matrix.
"""
function skewmat!(S, theta)
    @assert  length(theta)== 3 "Input must be a 3-vector"
    # S[:,:].=[0.0     -theta[3]  theta[2];
    #          theta[3]   0.0    -theta[1];
    #         -theta[2] theta[1]  0.0];
    S[1, 1] = S[2, 2] = S[3, 3] = 0.0
    S[1, 2] = -theta[3]
    S[1, 3] = theta[2];
    S[2, 1] = theta[3]
    S[2, 3] = -theta[1];
    S[3, 1] = -theta[2]
    S[3, 2] = theta[1]
    return S;
end

# function cross3(theta::AbstractVector{T}, v::AbstractVector{T}) where T
#     @assert (length(theta)== 3) && (length(v)== 3) "Inputs must be 3-vectors"
#     return [ 0.0-theta[3]*v[2]+theta[2]*v[3];
#             theta[3]*v[1]+0.0-theta[1]*v[3];
#             -theta[2]*v[1]+theta[1]*v[2]+0.0];
# end

function cross3!(result::AbstractVector{T}, theta::AbstractVector{T}, v::AbstractVector{T}) where T
    @assert (length(theta)== 3) && (length(v)== 3) "Inputs must be 3-vectors"
    result[1] = -theta[3]*v[2]+theta[2]*v[3];
    result[2] = theta[3]*v[1]-theta[1]*v[3];
    result[3] = -theta[2]*v[1]+theta[1]*v[2];
    return result
end

function cross2(theta::AbstractVector{T}, v::AbstractVector{T}) where T
    @assert (length(theta)== 2) && (length(v)== 2) "Inputs must be 2-vectors"
    return -theta[2]*v[1] +theta[1]*v[2];
end

# """
#     rotmat(theta::FFltMat)
#
# Compute rotation matrix from a rotation vector or from the associated skew
# matrix.
#
# Example:
# ```
# theta=[0.8147;  0.9058;       0.1270]
# julia>        RotationUtilModule.rotmat(theta')
# 3x3 Array{Float64,2}:
#   0.631452  0.22755    0.741276
#   0.422633  0.700499  -0.57505
#  -0.650116  0.676404   0.346161
# ```
# """
# function rotmat(theta::FFltMat)
#     m,n=size(theta);
#     if (m == 3) && (n == 3)         # in this case we are assuming  the input is a skew symmetric matrix
#         # if (norm(theta' + theta) > 1e-9)
#         #     error("theta must be a skew-symmetric matrix");
#         # end
#         a=[ -theta[2,3],  theta[1,3],  -theta[1,2]];
#         # R = expm(thetatilde); # but we are not going to do that, we will use the formula below
#     else                            # we are going to assume that the input is the rotation vector
#         a= theta[:];
#     end
#     na = norm(a);
#     if (na == 0)
#         R = eye(3,3);
#     else
#         thetatilde=zeros(3,3); skewmat!(thetatilde, a);
#         a=a/na;
#         ca = cos(na);
#         sa = sin(na);
#         aa = a*a';
#         R = ca * (eye(3,3)-aa) + sa/na*thetatilde + aa;
#     end
#     return R;
# end

end


text:
using PyPlot
using QuantEcon: mc_compute_stationary

P =[0.971 0.029 0.000
    0.145 0.778 0.077
    0.000 0.508 0.492]

psi = [0.0 0.2 0.8]

fig = figure()
ax = fig[:gca](projection="3d")
ax[:set_xlim](0, 1)
ax[:set_ylim](0, 1)
ax[:set_zlim](0, 1)
ax[:set_xticks]((0.25, 0.5, 0.75))
ax[:set_yticks]((0.25, 0.5, 0.75))
ax[:set_zticks]((0.25, 0.5, 0.75))

t = 20
x_vals = Array(Float64, t)
y_vals = Array(Float64, t)
z_vals = Array(Float64, t)

for i=1:t
    x_vals[i] = psi[1]
    y_vals[i] = psi[2]
    z_vals[i] = psi[3]
    psi = psi*P
end

ax[:scatter](x_vals, y_vals, z_vals, c="r", s = 60)

P = MarkovChain(P)
psi_star = mc_compute_stationary(P)
ax[:scatter](psi_star[1], psi_star[2], psi_star[3], c="k", s = 60)


text:
module Transformation

using Images
using Interpolations

using DeformableRegistration.ImageProcessing

export getCellCenteredGrid, getNodalGrid, stg2cen, cen2stg, checkForOddNumberOfGridPoints, getCellCenteredGridRanges
export transformGridAffine, transformWorldCoordinate, scaledArray, + , - , *

struct scaledArray
           data::Array{Float64, 1}
           dimensions::Tuple{Int64, Int64}
           voxelsize::Array{Float64, 1}
           shift::Array{Float64, 1}
end

import Base.+
import Base.-
import Base.*

function +(a::scaledArray, b::scaledArray)
    @assert a.voxelsize == b.voxelsize "Array's world matrices must match"
    @assert a.shift == b.shift "Array's world matrices must match"
    return(scaledArray(a.data+b.data, a.dimensions, a.voxelsize, a.shift))
end

function -(a::scaledArray, b::scaledArray)
    return(a + (-1 * b))
end

function *(s::Number, a::scaledArray)
    b = deepcopy(a)
    b.data[:] = s * b.data[:]
    return b
end

function getCellCenteredGrid(I::regImage)
  return getCellCenteredGrid(I.voxelsize,I.shift,size(I.data))
end

function getCellCenteredGrid(I::scaledArray)
  return getCellCenteredGrid(I.voxelsize,I.shift,I.dimensions)
end

function getCellCenteredGrid(voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64})
    x = Array(shift[1]+voxelsize[1]/2:voxelsize[1]:shift[1]+voxelsize[1]*(gridSize[1]))
    y = Array(shift[2]+voxelsize[2]/2:voxelsize[2]:shift[2]+voxelsize[2]*(gridSize[2]))
    xgrid = repmat(x,1,gridSize[2])[:]
    ygrid = repmat(y',gridSize[1],1)[:]
    return scaledArray(vcat(xgrid,ygrid), gridSize, voxelsize, shift)
end

function getNodalGrid(image::regImage)
  return getNodalGrid(I.voxelsize,I.shift, size(I.data))
end

function getNodalGrid(voxelsize::Array{Float64,1},shift::Array{Float64,1}, gridSize::Tuple{Int64,Int64})
    x = Array(shift[1]:voxelsize[1]:shift[1]+voxelsize[1]*(gridSize[1]))
    y = Array(shift[2]:voxelsize[2]:shift[2]+voxelsize[2]*(gridSize[2]))
    xgrid = repmat(x,1,length(y))[:]
    ygrid = repmat(y',length(x),1)[:]
    return scaledArray(vcat(xgrid,ygrid), gridSize, voxelsize, shift)
end

# # staggered grid to Cellcentered grid
# function stg2cen(staggeredGrid::Array{Float64,1},gridSize::Tuple{Int,Int})
#     N = prod(gridSize)
#     M = gridSize[1]*(gridSize[2]+1)
#     centeredGrid = zeros(2*N)
#     for j=1:gridSize[2]
#         for i=1:gridSize[1]
#             centeredGrid[  (j-1)*gridSize[1]+i] = (staggeredGrid[(j-1)* (gridSize[1]+1)+i]      + staggeredGrid[j*gridSize[1]+i]            ) / 2
#             centeredGrid[N+(j-1)*gridSize[1]+i] = (staggeredGrid[M+(j-1)*(gridSize[1]+1)+i] + staggeredGrid[M+(j-1)*(gridSize[1]+1)+i+1]) / 2
#         end
#     end
#     return centeredGrid
# end
#
# # Cellcentered grid to staggeredGrid
# function cen2stg(centeredGrid::Array{Float64,1},gridSize::Array{Int64,1})
#     NX = (gridSize[2]+1)*gridSize[1]
#     NY = (gridSize[1]+1)*gridSize[2]
#     M = prod(gridSize)
#     staggeredGrid = zeros(NX+NY)
#
#     # staggered grid in horizontal direction
#     # inner part
#     for i=gridSize[1]+1:M
#         staggeredGrid[i] = ( centeredGrid[i-gridSize[1]] + centeredGrid[i] ) / 2
#     end
#     # outer part (edge area)
#     for i=1:gridSize[1]
#         staggeredGrid[i] = centeredGrid[i] / 2
#     end
#     for i=M+1:NX
#         staggeredGrid[i] = centeredGrid[i-gridSize[1]] / 2
#     end
#
#     # staggered grid in vertical direction
#     # inner part
#     for j=1:gridSize[2]
#         staggeredGrid[NX+(j-1)*(gridSize[1]+1)+1] = centeredGrid[M+1+(j-1)*gridSize[1]] / 2
#         for i=2:gridSize[1]
#             staggeredGrid[NX+(j-1)*(gridSize[1]+1)+i] = (centeredGrid[M+i+(j-1)*gridSize[1]-1] + centeredGrid[M+i+(j-1)*gridSize[1]]) / 2
#         end
#         staggeredGrid[NX+j*(gridSize[1]+1)] = centeredGrid[M+j*gridSize[1]] / 2
#     end
#
#     return staggeredGrid
# end
#
# function cen2stg(distanceOutput::Tuple{Number,Array{Float64,1},Function,Tuple{Array{Float64,1},Array{Float64,1}}}, refImg::ImageMeta)
#   D = distanceOutput[1]
#   dD = cen2stg(distanceOutput[2],getSize(refImg))
#   d2D(grid) = cen2stg(distanceOutput[3](stg2cen(grid,getSize(refImg))),getSize(refImg))
#   dTransformedImage = distanceOutput[4]
#   return(D, dD, d2D, dTransformedImage)
# end
#
# function cen2stg(distanceOutput::Tuple{Number,Array{Float64,1},Function}, refImg::ImageMeta)
#   D = distanceOutput[1]
#   dD = cen2stg(distanceOutput[2],getSize(refImg))
#   d2D(grid) = cen2stg(distanceOutput[3](stg2cen(grid,getSize(refImg))),getSize(refImg))
# 	return(D, dD, d2D)
# end

# function cen2stg(distanceOutput::Tuple{Number,Array{Float64,1},SparseMatrixCSC{Float64,Int64},Tuple{Array{Float64,1},Array{Float64,1}}}, refImg::ImageMeta)
#   D = distanceOutput[1]
#   dD = cen2stg(distanceOutput[2],getSize(refImg))
#   d2D(grid) = cen2stg(distanceOutput[3]*stg2cen(grid,getSize(refImg)),getSize(refImg))
#   dTransformedImage = distanceOutput[4]
#   return(D, dD, d2D, dTransformedImage)
# end
#
# function cen2stg(distanceOutput::Tuple{Number,Array{Float64,1},SparseMatrixCSC{Float64,Int64}}, refImg::ImageMeta)
#   D = distanceOutput[1]
#   dD = cen2stg(distanceOutput[2],getSize(refImg))
#   d2D(grid) = cen2stg(distanceOutput[3]*stg2cen(grid,getSize(refImg)),getSize(refImg))
#   return(D, dD, d2D)
# end
#
# function cen2stg(distanceOutput::Tuple{Number,Any,Any,Tuple{Array{Float64,1},Array{Float64,1}}}, refImg::ImageMeta)
#   D = distanceOutput[1]
#   dD = distanceOutput[2]
#   d2D = distanceOutput[3]
#   dTransformedImage = distanceOutput[4]
#   return(D, dD, d2D, dTransformedImage)
# end
#
# function cen2stg(distanceOutput::Tuple{Number,Number,Number}, refImg::ImageMeta)
#   D = distanceOutput[1]
#   dD = distanceOutput[2]
#   d2D = distanceOutput[3]
#   return(D, dD, d2D)
# end
#
function checkForOddNumberOfGridPoints(I::ImageMeta,grid::Array{Float64,1})
  if( 2*prod(size(I)) == size(grid,1) )
      return false
  else
      return true
  end
end

function getCellCenteredGridRanges(I::regImage, gridSize)
    x = Array(I.shift[1]+I.voxelsize[1]/2:I.voxelsize[1]:I.shift[1]+I.voxelsize[1]*(gridSize[1]))
    y = Array(I.shift[2]+I.voxelsize[2]/2:I.voxelsize[2]:I.shift[2]+I.voxelsize[2]*(gridSize[1]))
    return x,y
end

function transformWorldCoordinate(I::regImage,y::Float64,x::Float64)
    xNew = (x - I.shift[1]) / I.voxelsize[1] + .5
    yNew = (y - I.shift[2]) / I.voxelsize[2] + .5
    return yNew,xNew
end

function transformWorldCoordinate(I::regImage,grid::Array{Float64,1})
    N = size(grid,1)
    newGrid = zeros(N)
    for i=1:N
      newGrid[i] = (grid[i] - I.shift[1]) / I.voxelsize[1] + .5
      newGrid[N+i] = (grid[N+i] - I.shift[2]) / I.voxelsize[2] + .5
    end
    return newGrid
end

function transformGridAffine(grid::scaledArray,affineParameters)
    transformationMatrix = [affineParameters[1] affineParameters[2]; affineParameters[4] affineParameters[5]]
    translationParameters = [affineParameters[3]; affineParameters[6]]
    return transformGridAffine(transformationMatrix, grid, translationParameters)
end

function transformGridAffine(transformationMatrix::Matrix,gridOriginal::scaledArray,translationParameters::Array{Float64,1})
    N = prod(gridOriginal.dimensions)
    grid = zeros(2*N)
    @simd for i=1:N
        @inbounds grid[i]   = transformationMatrix[1,1]*gridOriginal.data[i]+transformationMatrix[1,2]*gridOriginal.data[N+i]+translationParameters[1]
        @inbounds grid[N+i] = transformationMatrix[2,1]*gridOriginal.data[i]+transformationMatrix[2,2]*gridOriginal.data[N+i]+translationParameters[2]
    end
    return scaledArray(grid, gridOriginal.dimensions, gridOriginal.voxelsize, gridOriginal.shift)
end

end


text:
# Models.jl
"""
Type of IRT models. User-defined model can be used.
A struct which hyper type is `IRTmodel` must contain these containers below:

- Model parameter(s)
- group
- fixed

Then, some methods should be implemented:

- `_take` Extract model item parameter to be estimated in EM cycle.
- `_distribute` and `_distribute!` Attach new value to contaiers of parameters.
- `_irf` Item response function. 2 type of return value, one is Float64,
anothoer is eltype(par), which is used in ForwardDiff, should be implemented.
- `_init` and `_init!` Initialize item parameters, which are used for the starting value in the estimation.
- `_rescale!` Rescaling item parameters at each EM cycles.
- `_checkitem` Check whether estimated parameters are taking appropriate values.

See "IRTmodels/graded.jl" for an example.

"""
abstract type IRTmodel end

"""
Take model paraemters.
"""
function take(par::IRTmodel)
    _take(par)
end

"""
Take model parameters based on `q` matrix(vector).
"""
function take(par::IRTmodel, q)
    _take(par, q)
end

"""
Distribute `new` value to `old` as model parameters.
"""
function distribute!(new, old::IRTmodel)
    _distribute!(new, old)
end

function distribute!(new, old::IRTmodel, q)
    _distribute!(new, old, q)
end

function distribute(new, old::IRTmodel)::IRTmodel
    _distribute(new, old)
end

function distribute(new, old::IRTmodel, q)::IRTmodel
    _distribute(new, old, q)
end

function copy(x::IRTmodel)
    _copy(x)
end

function rescale!(pars::IRTmodel, μ, Σ, q)
    _rescale!(pars, μ, Σ, q)
end

"""
    irf(p::graded, θ, u, q)
Item respose function

# Arguments

`p` Probability model contains its parameters.
`θ` A latent trait vector
`q` Q matrix(vector).
`u` An item response.

# Example
```julia
# Graded response model
pars = graded(a = 1.0, d = [2, 0, -1.5])
[irf.(Ref(pars), -4:0.1:4, k, 1) for k in 0:3] |> plot
```
"""
function irf(p::IRTmodel, θ, u, q)
    _irf(p, θ, u, q)
end

# For the internal use in ForwardDiff
function irf(p, m::IRTmodel, θ, u, q)
    _irf(p, m, θ, u, q)
end


text:
immutable Multinomial <: DiscreteMultivariateDistribution
    n::Int
    prob::Vector{Float64}
    drawtable::DiscreteDistributionTable
    function Multinomial{T <: Real}(n::Integer, p::Vector{T})
        if n <= 0
            error("Multinomial: n must be positive")
        end
        sump = 0.0
        for i in 1:length(p)
            if p[i] < 0.0
                error("Multinomial: probabilities must be non-negative")
            end
            sump += p[i]
        end
        for i in 1:length(p)
            p[i] /= sump
        end
        new(int(n), p, DiscreteDistributionTable(p))
    end
end

function Multinomial(n::Integer, d::Integer)
    if d <= 1
        error("d must be greater than 1")
    end
    prob = Array(Float64, d)
    fill!(prob, 1.0 / d)
    Multinomial(n, prob)
end

Multinomial(d::Integer) = Multinomial(1, d)

entropy(d::Multinomial) = pventropy(d.prob)

function insupport{T <: Real}(d::Multinomial, x::Vector{T})
    n = length(x)
    if length(d.prob) != n
        return false
    end
    s = 0.0
    for i in 1:n
        if x[i] < 0.0 || !isinteger(x[i])
            return false
        end
        s += x[i]
    end
    if abs(s - d.n) > 1e-8
        return false
    end
    return true
end

mean(d::Multinomial) = d.n .* d.prob

function mgf(d::Multinomial, t::AbstractVector)
    p, n = d.prob, d.n
    k = length(p)
    s = 0.0
    for i in 1:k
        s += p[i] * exp(t[i])
    end
    return s^n
end

function cf(d::Multinomial, t::AbstractVector)
    p, n = d.prob, d.n
    k = length(p)
    s = 0.0 + 0.0im
    for i in 1:k
        s += p[i] * exp(im * t[i])
    end
    return s^n
end

pdf{T <: Real}(d::Multinomial, x::Vector{T}) = exp(logpdf(d, x))

function logpdf{T <: Real}(d::Multinomial, x::Vector{T})
    if !insupport(d, x)
        return -Inf
    else
        s = lgamma(d.n + 1.0)
        for i in 1:length(x)
            s -= lgamma(x[i] + 1.0)
            s += x[i] * log(d.prob[i])
        end
        return s
    end
end

function rand!(d::Multinomial, x::Vector)
    for itr in 1:d.n
        i = draw(d.drawtable)
        x[i] += 1
    end
    return x
end

function rand(d::Multinomial)
    x = zeros(Int, length(d.prob))
    return rand!(d, x)
end

function var(d::Multinomial)
    n = length(d.prob)
    S = Array(Float64, n, n)
    for j in 1:n
        for i in 1:n
            if i == j
                S[i, j] = d.n * d.prob[i] * (1.0 - d.prob[i])
            else
                S[i, j] = -d.n * d.prob[i] * d.prob[j]
            end
        end
    end
    return S
end


text:
literal(i::Int) = @bf [right!, repeated(:inc!, i), right!, inc!]

@bf swap = [step!(-2), dec!, step!(-1), move!(1),
            step!(2), move!(-2),
            step!(-1), move!(1), inc!, step!(2)]

@bf dup = [dec!, step!(-1), move!(1),
           step!(1), move!(-1, 1), inc!,
           step!(2), inc!]

@bf empty! = [[drop], while!]

@bf halt = [[rstack!, Flip(:empty!)], call]

@bf ! = [left!, [right!, dec!, left!, reset!], while!, # Set flag, reset
         right!, [left!, inc!, right!, dec!], while!, inc!]

@bf + = [dec!, left!, move!(-2), left!]
@bf - = [dec!, left!, move!(-2, mode = :dec!), left!]

@bf * = [dec!, step!(-2), dec!, left!, move!(1), right!,
         [dec!, right!, move!(1, -2), right!, move!(-1), step!(-2)], while!,
         right!, reset!, left!, inc!]

@bf < = [left!, move!(-1), left!, dec!, left!, move!(2), right!,
         [right!, right!,
          if!(:pass, @bf [step!(-3), inc!, right!, reset!, step!(2)]),
          left!, dec!, left!, dec!], while!,
         inc!, step!(2), drop]

# End of unsafe operations

@bf != = [-]
@bf == = [-, !]
@bf or = [+]
@bf and = [*]
@bf <= = [[==], [<], bi2, or]
@bf > = [swap, <]
@bf >= = [swap, <=]

@bf sq = [dup, *]

# x f -- x
@bf dip = [swap, rpush, [rpop], rpush, call]
#  x y f -- x y
@bf dip2 = [swap, [dip], dip]

# x y -- y
@bf nip = [swap, drop]
# x y -- x x y
@bf dupd = [[dup], dip]
# x y -- x y x
@bf over = [dupd, swap]
# x y -- x y x y
@bf dup2 = [over, over]
# x y -- y x y
@bf tuck = [swap, over]
# x y z -- x y z x
@bf pick = [[over], dip, swap]
# x y z -- y x z
@bf swapd = [[swap], dip]
# x y z -- y z x
@bf rotl = [swapd, swap]
# x y z -- z x y
@bf rotr = [swap, swapd]

# x f -- x
@bf keep = [dupd, dip]
# x y f -- x y
@bf keep2 = [[dup2], dip, dip2]
# x f g -- fx gx
@bf bi = [[keep], dip, call]
# x y f g -- fxy gxy
@bf bi2 = [[keep2], dip, call]
# x y f g -- fx gy
@bf bi_ = [[dip], dip, call] # bi*
# x y f -- fx fy
@bf bia = [dup, bi_] # bi@

@bf iff = [rotl, [], [swap], if!, drop, call]
@bf when = [[], iff]
@bf unless = [[], swap, iff]

@bf loop = [over, dip2, rotl,
            [dup, dip2, loop], [drop, drop], iff]

@bf doloop = [dup, dip2, loop]

max_stars_repo_name: JuliaTagBot/BrainForth

max_issues_repo_name: JuliaTagBot/BrainForth

max_forks_repo_name: JuliaTagBot/BrainForth


text:
mutable struct whole_body_data_t <: LCMType
    timestamp::Int64

    num_positions::Int32
    q_des::Vector{Float32}

    spline::piecewise_polynomial_t

    num_constrained_dofs::Int32
    constrained_dofs::Vector{Int32}
end

@lcmtypesetup(whole_body_data_t,
    q_des => (num_positions,),
    constrained_dofs => (num_constrained_dofs,),
)


text:
@testset "SubSVDD" begin
    dummy_data, labels = generate_mvn_with_outliers(4, 100, 42, true, true)
    pools = fill(:U, size(dummy_data, 2))
    subspaces = [[1, 2], [3, 4]]

    model = SubSVDD(dummy_data, subspaces, pools)

    @testset "create" begin
        @test model.state == SVDD.model_created
        @test_throws SVDD.ModelStateException SVDD.predict(model, 1)
    end
    C = 1.0
    gamma = 0.5
    gamma_strategy = FixedGammaStrategy(MLKernels.GaussianKernel(gamma))
    C_strategy = FixedCStrategy(C)
    init_strategy = SVDD.SimpleSubspaceStrategy(gamma_strategy, C_strategy, gamma_scope=Val(:Global))

    @testset "initialize" begin
        SVDD.initialize!(model, init_strategy)
        @test length(model.kernel_fct) == 2
        @test !(all(model.K[1] .≈ model.K[2]))
        @test model.state == SVDD.model_initialized
        @test model.data == dummy_data
        @test model.C ≈ C
        @test length(model.K) == length(subspaces)
        @test size(model.K[1]) == (size(dummy_data, 2), size(dummy_data, 2))
        @test all(map(x -> all(diag(x) .≈ 1.0), model.K))
    end

    @testset "fit" begin
        SVDD.fit!(model, TEST_SOLVER)
        @test model.state == SVDD.model_fitted
        @test length(model.alpha_values) == length(subspaces)
        @test all(length.(model.alpha_values) .== size(model.data, 2))
        @test all(model.R .> 0.0)
        @test all(model.const_term .> 0.0)
    end

    @testset "calculate_c_delta" begin
        @test SVDD.calculate_upper_limit(model.alpha_values, model.C, model.v, 1) ≈ model.C .- model.alpha_values[2]
        @test SVDD.calculate_upper_limit(model.alpha_values, model.C, model.v, 2) ≈ model.C .- model.alpha_values[1]
        model.v .= 0.5
        @test SVDD.calculate_upper_limit(model.alpha_values, model.C, model.v, 2) ≈ 0.5*model.C .- model.alpha_values[1]
    end

    @testset "find_support_vectors" begin
        dummy_model = SubSVDD(dummy_data, [[1,2]], pools)
        dummy_model.alpha_values = [[0.1, 0.1-1e-10, 0.05, 0.0, 1e-10]]
        dummy_model.v = ones(Float64, 5)
        dummy_model.C = 0.1
        @test SVDD.find_support_vectors(dummy_model, 1) == [3]

        dummy_model.C = 0.2
        dummy_model.subspaces = [[1,2], [3,4]]
        push!(dummy_model.alpha_values, [0.09, 0.0, 0.0, 0.0, 0.0])
        @test SVDD.find_support_vectors(dummy_model, 1) == [1,2,3]
    end

    @testset "@each_subspace" begin
        @testset "generic function" begin
            g(m, subspace_idx) = subspace_idx
            expected = [1, 2]
            actual = @eachsubspace g(model)
            @test expected == actual
        end
        @testset "predict" begin
            expected = @eachsubspace SVDD.predict(model)
            actual = [SVDD.predict(model, k) for k in eachindex(model.subspaces)]
            @test all(expected .≈ actual)
        end
    end

    @testset "predict" begin
        predictions = @eachsubspace SVDD.predict(model)
        @test length(predictions) == 2
        @test all(length.(predictions) .== size(model.data, 2))
        @test all(map(x -> all(x .< 1e-5), predictions))
    end

    @testset "update_with_feedback" begin
        model.v .= 1.0
        expected_pools = labelmap2vec(model.pools)
        dummy_remaining_indices = collect(eachindex(expected_pools))

        @assert model.weight_update_strategy === nothing
        @test_throws ErrorException update_with_feedback!(model, model.data, expected_pools, [1], dummy_remaining_indices, dummy_remaining_indices)

        update_strategy = SVDD.FixedWeightStrategy(1.1, 0.9)
        set_param!(model, Dict(:weight_update_strategy => update_strategy))
        @test model.weight_update_strategy == update_strategy

        update_with_feedback!(model, model.data, expected_pools, Int[], dummy_remaining_indices, dummy_remaining_indices)
        @test labelmap2vec(model.pools) == expected_pools

        query_ids = [1,2]
        query_pool_labels = [:Lin, :Lout]
        expected_pools[query_ids] .= query_pool_labels
        update_with_feedback!(model, model.data, expected_pools, query_ids, dummy_remaining_indices, dummy_remaining_indices)
        @test model.v[1] ≈ 1.1
        @test model.v[2] ≈ 0.9
        @test all(model.v[3:end] .≈ 1.0)
        @test labelmap2vec(model.pools) == expected_pools
    end
end


text:
using Plots
using LinearAlgebra
using ForwardDiff
using ProgressMeter
using NearestNeighbors

using ContactAlgorithms

# This is a one-file project. Simplicity is key.

# Setup

# walls
width = 10.
height = 10.

w1 = (pos=[0., 0.], dir=[0.,1.])
w2 = (pos=[0., height], dir=[0.,-1.])
w3 = (pos=[0., 0.], dir=[1.,0.])
w4 = (pos=[width, 0.], dir=[-1.,0.])
walls = (w1,w2,w3,w4)

# initial state
N = 30
sqN = Int(round(sqrt(N)))
N = sqN * sqN
d = 2
x0 = hcat( [[x,y]
            for x in range(0.25*width, stop=0.75*width, length=sqN)
            for y in range(0.15*height, stop=0.75*height, length=sqN)]...  )
x0 += 0.1 * randn(d,N)
v0 = zeros(d,N)
R = 0.3

M = ones(N)' # mass
M_inv = M.^(-1) # inverse mass
μ = 0.1

# small steps stuff
d0 = zeros(1,N)
d0wall = zeros(1,N)
v_max = 5
# collision constraints
function c_coll(x, y, R)
    return sqrt(sum((x-y).^2)) - 2*R
end

function c_coll_vmax(x, y, R, i)
    return sqrt(sum((x-y).^2)) - 2*R+max(d0[i]-v_max*dtₛ,0)
end
# wall constraints
function c_wall_vmax(x, wall,i)
    return dot( x - wall.pos, wall.dir )+max(d0wall[i]-v_max*dtₛ,0)
end

# forces
F(x) = [0., -5.]
#F(x) = [5., 5.] .- x

scatter(x0[1,:], x0[2,:], markersize=R* (11*3.3),
        xlims=(-0.5,10.5), ylims=(-0.5,10.5), aspect_ratio=:equal)
savefig("initial_plot.png")

anim = Animation()
show_frames = true

dt = 0.5
#n_stab = 2
n_steps = 10 # internal steps
dtₛ = dt / n_steps
t_end = 10.
α = 0.0
γ = 0.3
R_collision_detection = 4. * R

# simulate
let x0=x0, v0=v0
    global t, X, V, d0, d0wall
    t = 0
    X = [x0]
    V = [v0]

    p = Progress(Int(t_end/dt))
    while t < t_end
        t += dt

        x = X[end]
        v = V[end]
        # detect collisions
        balltree = KDTree(x, reorder=false)

        idxs = inrange(balltree, x, R_collision_detection, false)
        for i = 1:N
            for j = idxs[i]
                if j < i
                    c = c_coll(x[:,i], x[:,j], R)
                    if c < 0
                        d0[i] = -c
                    end
                end
            end
        end
        for wall in (w1,w2,w3,w4)
            for i = 1:N
                c = c_wall(x[:,i], wall)
                if c < 0
                    d0wall[i] = -c
                end
            end
        end

        # solve collisions
        for n = 1:n_steps
            # predict next step
            x_old, x = x, x + dtₛ * v + dtₛ^2 * M_inv .* ( F(x) .- γ * v)
            #x_old, x = x, x + dtₛ * M_inv .* F(x)

            # wall constraints
            for wall in (w1,w2,w3,w4)
                for i = 1:N
                    c = c_wall_vmax(x[:,i], wall,i)
                    if c < 0
                        Dc = ForwardDiff.gradient(z -> c_wall_vmax(z,wall,i), x[:,i])
                        Δλ = -c / ( Dc' * M_inv[i] * Dc + α )
                        Δx = M_inv[i] * Dc * Δλ
                        x[:,i] += Δx
                    end
                end
            end

            # particles constraints
            for i = 1:N
                for j = idxs[i]
                    if j < i
                        c = c_coll_vmax(x[:,i], x[:,j], R,i)
                        if c < 0
                            Dcᵢ = ForwardDiff.gradient(z -> c_coll_vmax(z, x[:,j], R, i), x[:,i])
                            Dcⱼ = ForwardDiff.gradient(z -> c_coll_vmax(x[:,i], z, R, i), x[:,j])
                            denom = (Dcᵢ)' * M_inv[i] * Dcᵢ + (Dcⱼ)' * M_inv[j] * Dcⱼ + α
                            Δλᵢ = -c / denom
                            Δλⱼ = -c / denom
                            Δxᵢ = M_inv[i] * Dcᵢ * Δλᵢ
                            Δxⱼ = M_inv[j] * Dcⱼ * Δλⱼ
                            x[:,i] += Δxᵢ
                            x[:,j] += Δxⱼ
                        end
                    end
                end
            end

            # update velocity
            v = (x - x_old) / dtₛ
            #v += dtₛ * M_inv .* F(x);
        end

        push!(X, x)
        push!(V, v)

        if show_frames
            scatter(X[end][1,:], X[end][2,:], markersize=R*(11*3.3),
                    xlims=(-0.5,10.5), ylims=(-0.5,10.5), aspect_ratio=:equal)
            frame(anim)
            display(current())
        end

        next!(p)
    end
end

if !show_frames
    anim = @animate for i = 1:length(X)
        scatter(X[i][1,:], X[i][2,:], markersize=R* (11*3.3),
                xlims=(-0.5,10.5), ylims=(-0.5,10.5), aspect_ratio=:equal)
    end
end

gif(anim, "output.gif", fps=10)

max_stars_repo_name: SteffenPL/ContactAlgorithms

max_issues_repo_name: SteffenPL/ContactAlgorithms


text:
function bitcat(::Type{T}, arr)::T where T<:Number
    out = zero(T)
    for x in arr
        out <<= sizeof(x)*8
        out |= convert(T, x)
    end
    out
end

function geteventdata(TYPE::Symbol, data::Vector{UInt8}, pairs::Pair{Symbol,DataType}...)
    e = Event{TYPE}()
    i = 1
    for pair in pairs
        s = sizeof(pair.second)
        setindex!(getfield(e, :fields), bitcat(pair.second, data[i+s-1:-1:i]), pair.first)
        i += s
    end
    return e
end

eventtype(data::Vector{UInt8}) = bitcat(UInt32, data[4:-1:1])

parseevent(window::Window, data::Vector{UInt8}) = parseevent(window, Val(eventtype(data)), data)
parseevent(window::Window, ::Val, data::Vector{UInt8}) = Event{:notsupported}()

function parseevent(window::Window, ::Val{SDL.QUIT}, data::Vector{UInt8})
    return geteventdata(:quit, data, :type => UInt32, :timestamp => UInt32)
end


text:
export TikhonovReg, getTikhonovReg, regularizer

"""
Tikhnov Regularizer R(x) = 0.5*alpha*(H*(x-xref), x-xref), with H spd
"""
mutable struct TikhonovReg{T}
    H
    alpha::T
    xref::Array{T}
end

getTikhonovReg(TYPE;H=I,alpha=1.0,xref=zeros(0))=TikhonovReg{TYPE}(H,convert(TYPE,alpha),convert.(TYPE,vec(xref)))

function regularizer(this::TikhonovReg,x)
    u = vec(copy(x))
    if !isempty(this.xref)
        u.-=this.xref
    end
    d2S = this.alpha*this.H
    dS = d2S*u
    S  = 0.5*dot(dS,u)
    return S,dS,d2S
end


text:
# addprocs(3)
using Test
# using Compat
# using IncrementalInference

include("testApproxConv.jl")

include("testHasPriors913.jl")

include("testInitVariableOrder.jl")

include("testExpXstroke.jl")

include("testBasicRecycling.jl")

include("TestModuleFunctions.jl")

include("testStateMachine.jl")

include("testCompareVariablesFactors.jl")

include("typeReturnMemRef.jl")

include("basicGraphsOperations.jl")

include("testMixturePrior.jl")

include("testPartialFactors.jl")

include("testBayesTreeiSAM2Example.jl")

include("testSpecialSampler.jl")

include("testSaveLoadDFG.jl")

include("testJunctionTreeConstruction.jl")

#FIXME fails on MetaBayesTree
include("testTreeSaveLoad.jl")

include("saveconvertertypes.jl")

include("testgraphpackingconverters.jl")

include("testNLsolve.jl")

# Randomized roots no longer supported, see PR #1075
# include("testNumericRootGenericRandomized.jl")

include("testCommonConvWrapper.jl")

include("testBasicForwardConvolve.jl")

include("testFactorMetadata.jl")

include("testBasicCSM.jl")

include("testCliqueFactors.jl")

include("testCcolamdOrdering.jl")

include("testBasicGraphs.jl")

include("testDefaultDeconv.jl")

include("testJointEnforcement.jl")

include("testlocalconstraintexamples.jl")

include("testBasicTreeInit.jl")

include("testSolveOrphanedFG.jl")

include("testSolveSetPPE.jl")

include("testEuclidDistance.jl")

include("priorusetest.jl")

include("testVariousNSolveSize.jl")

include("testExplicitMultihypo.jl")

include("TestCSMMultihypo.jl")

include("testMultihypoFMD.jl")

include("testMultiHypo2Door.jl")

include("testMultimodal1D.jl")

include("testMultihypoAndChain.jl")

include("testMultithreaded.jl")

include("testpartialconstraint.jl")

include("testnullhypothesis.jl")

include("testmultihypothesisapi.jl")

include("fourdoortest.jl")

include("testSphere1.jl")

include("testMixtureLinearConditional.jl")

include("testFluxModelsDistribution.jl")

include("testDERelative.jl")

include("testAnalysisTools.jl")

include("testBasicParametric.jl")

# dont run test on ARM, as per issue #527
if Base.Sys.ARCH in [:x86_64;]
  include("testTexTreeIllustration.jl")
end

include("testMultiprocess.jl")

include("testDeadReckoningTether.jl")

include("testCSMMonitor.jl")

include("testSkipUpDown.jl")

include("testTreeFunctions.jl")

#